C +-======-+ C Copyright (c) 2003-2018 United States Government as represented by C the Admistrator of the National Aeronautics and Space Administration. C All Rights Reserved. C C THIS OPEN SOURCE AGREEMENT ("AGREEMENT") DEFINES THE RIGHTS OF USE, C REPRODUCTION, DISTRIBUTION, MODIFICATION AND REDISTRIBUTION OF CERTAIN C COMPUTER SOFTWARE ORIGINALLY RELEASED BY THE UNITED STATES GOVERNMENT AS C REPRESENTED BY THE GOVERNMENT AGENCY LISTED BELOW ("GOVERNMENT AGENCY"). C THE UNITED STATES GOVERNMENT, AS REPRESENTED BY GOVERNMENT AGENCY, IS AN C INTENDED THIRD-PARTY BENEFICIARY OF ALL SUBSEQUENT DISTRIBUTIONS OR C REDISTRIBUTIONS OF THE SUBJECT SOFTWARE. ANYONE WHO USES, REPRODUCES, C DISTRIBUTES, MODIFIES OR REDISTRIBUTES THE SUBJECT SOFTWARE, AS DEFINED C HEREIN, OR ANY PART THEREOF, IS, BY THAT ACTION, ACCEPTING IN FULL THE C RESPONSIBILITIES AND OBLIGATIONS CONTAINED IN THIS AGREEMENT. C C Government Agency: National Aeronautics and Space Administration C Government Agency Original Software Designation: GSC-15354-1 C Government Agency Original Software Title: GEOS-5 GCM Modeling Software C User Registration Requested. Please Visit http://opensource.gsfc.nasa.gov C Government Agency Point of Contact for Original Software: C Dale Hithon, SRA Assistant, (301) 286-2691 C C +-======-+ subroutine windfix ( ua,va,plea, . ub,vb,pleb,im,jm,lm,lattice,grid,method ) use dynamics_lattice_module implicit none type ( dynamics_lattice_type ) lattice integer im,jm,lm,method real ua(im,jm,lm) real va(im,jm,lm) real plea(im,jm,lm+1) real ub(im,jm,lm) real vb(im,jm,lm) real pleb(im,jm,lm+1) real u(im,jm,lm) real v(im,jm,lm) real dpa(im,jm,lm) real dpb(im,jm,lm) real diva(im,jm,lm) real divb(im,jm,lm) character(*) grid integer index(lm),ierror real, allocatable :: uglo(:,:,:) real, allocatable :: vglo(:,:,:) real, allocatable :: dglo(:,:,:) real, allocatable :: dpglo(:,:,:) real, allocatable :: sum1(:,:) real, allocatable :: sum2(:,:) real, allocatable :: sum3(:,:) real, allocatable :: lambda(:,:) real, allocatable :: chi(:,:) real, allocatable :: uchi(:,:) real, allocatable :: vchi(:,:) integer L, comm, myid, npes integer img, jmg img = lattice%imglobal jmg = lattice%jmglobal comm = lattice%comm myid = lattice%myid npes = lattice%nx * lattice%ny do L=1,lm index(L) = mod(L-1,npes) enddo allocate ( uglo(img,jmg,lm) ) allocate ( vglo(img,jmg,lm) ) allocate ( dglo(img,jmg,lm) ) allocate ( dpglo(img,jmg,lm) ) allocate ( chi(img,jmg) ) allocate ( uchi(img,jmg) ) allocate ( vchi(img,jmg) ) allocate ( sum1(im,jm) ) allocate ( sum2(im,jm) ) allocate ( sum3(im,jm) ) allocate ( lambda(im,jm) ) ! Compute Pressure Thickness ! -------------------------- do L=1,lm dpa(:,:,L) = ( plea(:,:,L+1)-plea(:,:,L) ) dpb(:,:,L) = ( pleb(:,:,L+1)-pleb(:,:,L) ) enddo ! Gather Winds for Background ! --------------------------- do L=1,lm call gather_2d ( uglo(1,1,L), ub(1,1,L),lattice ) call gather_2d ( vglo(1,1,L), vb(1,1,L),lattice ) call gather_2d ( dpglo(1,1,L),dpb(1,1,L),lattice ) if( trim(grid).eq.'d' .or. . trim(grid).eq.'D' ) then if( lattice%myid.eq.0 ) then call dtoa ( uglo(1,1,L),uglo(1,1,L),img,jmg,1,2 ) call dtoa ( vglo(1,1,L),vglo(1,1,L),img,jmg,1,1 ) endif endif enddo #ifdef mpi call mpi_bcast ( uglo,img*jmg*lm,lattice%mpi_rkind,0,comm,ierror ) call mpi_bcast ( vglo,img*jmg*lm,lattice%mpi_rkind,0,comm,ierror ) call mpi_bcast ( dpglo,img*jmg*lm,lattice%mpi_rkind,0,comm,ierror ) #endif c Compute Vorticity and Divergence c -------------------------------- do L=1,lm if( index(L).eq.myid ) call getdiv (uglo(1,1,L),vglo(1,1,L),dpglo(1,1,L),dglo(1,1,L),img,jmg ) enddo #ifdef mpi call mpi_barrier (comm,ierror) do L=1,lm call mpi_bcast ( dglo(1,1,L),img*jmg,lattice%mpi_rkind,index(L),comm,ierror ) enddo call mpi_barrier (comm,ierror) #endif do L=1,lm call scatter_2d ( dglo(1,1,L),divb(1,1,L),lattice ) enddo ! Gather Winds for Analysis ! ------------------------- do L=1,lm call gather_2d ( uglo(1,1,L), ua(1,1,L),lattice ) call gather_2d ( vglo(1,1,L), va(1,1,L),lattice ) call gather_2d (dpglo(1,1,L),dpa(1,1,L),lattice ) if( trim(grid).eq.'d' .or. . trim(grid).eq.'D' ) then if( lattice%myid.eq.0 ) then call dtoa ( uglo(1,1,L),uglo(1,1,L),img,jmg,1,2 ) call dtoa ( vglo(1,1,L),vglo(1,1,L),img,jmg,1,1 ) endif endif enddo #ifdef mpi call mpi_bcast ( uglo,img*jmg*lm,lattice%mpi_rkind,0,comm,ierror ) call mpi_bcast ( vglo,img*jmg*lm,lattice%mpi_rkind,0,comm,ierror ) call mpi_bcast ( dpglo,img*jmg*lm,lattice%mpi_rkind,0,comm,ierror ) #endif c Compute Vorticity and Divergence c -------------------------------- do L=1,lm if( index(L).eq.myid ) call getdiv (uglo(1,1,L),vglo(1,1,L),dpglo(1,1,L),dglo(1,1,L),img,jmg ) enddo #ifdef mpi call mpi_barrier (comm,ierror) do L=1,lm call mpi_bcast ( dglo(1,1,L),img*jmg,lattice%mpi_rkind,index(L),comm,ierror ) enddo call mpi_barrier (comm,ierror) #endif do L=1,lm call scatter_2d ( dglo(1,1,L),diva(1,1,L),lattice ) enddo c Compute Divergence Increment (to force vanishing vertical integral) c ------------------------------------------------------------------- c Minimize Relative Change c ------------------------ if( method.eq.1 ) then sum1(:,:) = 0.0 sum2(:,:) = 0.0 sum3(:,:) = 0.0 do L=1,lm sum1(:,:) = sum1(:,:) + diva(:,:,L) sum2(:,:) = sum2(:,:) + divb(:,:,L) sum3(:,:) = sum3(:,:) + (diva(:,:,L)-divb(:,:,L))**2 enddo where( sum3 .ne. 0.0 ) lambda = (sum1-sum2) / sum3 elsewhere lambda = 0.0 endwhere do L=1,lm diva(:,:,L) = -lambda(:,:)*( diva(:,:,L)-divb(:,:,L) )**2 enddo sum3(:,:) = sum1(:,:) do L=1,lm sum3(:,:) = sum3(:,:) + diva(:,:,L) enddo endif c Minimize Absolute Change c ------------------------ if( method.eq.2 ) then sum1(:,:) = 0.0 sum2(:,:) = 0.0 sum3(:,:) = 0.0 do L=1,lm sum1(:,:) = sum1(:,:) + diva(:,:,L) sum2(:,:) = sum2(:,:) + divb(:,:,L) sum3(:,:) = sum3(:,:) + dpa(:,:,L)**2 enddo lambda = (sum1-sum2) / sum3 do L=1,lm diva(:,:,L) = -lambda(:,:)*dpa(:,:,L)**2 enddo sum3(:,:) = sum1(:,:) do L=1,lm sum3(:,:) = sum3(:,:) + diva(:,:,L) enddo endif call writit ( sum1,im,jm,1,55,lattice ) call writit ( sum2,im,jm,1,55,lattice ) call writit ( sum3,im,jm,1,55,lattice ) c Gather and Broadcast Divergence Increment c ----------------------------------------- do L=1,lm call gather_2d ( dglo(1,1,L),diva(1,1,L),lattice ) enddo #ifdef mpi call mpi_bcast ( dglo,img*jmg*lm,lattice%mpi_rkind,0,comm,ierror ) #endif c Compute Wind Increments Associated with Divergence Increment c ------------------------------------------------------------ do L=1,lm if( index(L).eq.myid ) then call VELPOT (dglo(1,1,L),chi,img,jmg) call gradq (chi, uchi,vchi,img,jmg) uglo(:,:,L) = uglo(:,:,L) + uchi(:,:)/dpglo(:,:,L) vglo(:,:,L) = vglo(:,:,L) + vchi(:,:)/dpglo(:,:,L) endif enddo #ifdef mpi call mpi_barrier (comm,ierror) do L=1,lm call mpi_bcast ( uglo(1,1,L),img*jmg,lattice%mpi_rkind,index(L),comm,ierror ) call mpi_bcast ( vglo(1,1,L),img*jmg,lattice%mpi_rkind,index(L),comm,ierror ) enddo call mpi_barrier (comm,ierror) #endif ! Scatter Winds ! ------------- do L=1,lm if( trim(grid).eq.'d' .or. . trim(grid).eq.'D' ) then if( lattice%myid.eq.0 ) then call atod ( uglo(1,1,L),uglo(1,1,L),img,jmg,1,2 ) call atod ( vglo(1,1,L),vglo(1,1,L),img,jmg,1,1 ) endif endif call scatter_2d ( uglo(1,1,L),ua(1,1,L),lattice ) call scatter_2d ( vglo(1,1,L),va(1,1,L),lattice ) enddo deallocate ( sum1,sum2,sum3,lambda ) deallocate ( chi,uchi,vchi ) deallocate ( uglo,vglo,dglo,dpglo ) return end SUBROUTINE GETDIV ( U,V,DP,DIV,IM,JM ) C ******************************************************************** C **** **** C **** THIS PROGRAM CALCULATES DIVERGENCE **** C **** AT EACH LEVEL FOR A NON-STAGGERED A-GRID **** C **** **** C **** INPUT: **** C **** U ....... ZONAL WIND **** C **** V ....... MERIDIONAL WIND **** C **** IM ...... NUMBER OF LONGITUDE POINTS **** C **** JM ...... NUMBER OF LATITUDE POINTS **** C **** **** C **** OUTPUT: **** C **** DIV (IM,JM) .... DIVERGENCE **** C **** **** C ******************************************************************** real U(IM,JM) real V(IM,JM) real DP(IM,JM) real DIV(IM,JM) real P1X (IM,JM) real P1Y (IM,JM) real TMP1(IM,JM) real TMP2(IM,JM) real cosij(IM,JM) DIMENSION MSGN(2) DATA MSGN /-1,1/ C ********************************************************* C **** INITIALIZATION FOR DIVERGENCE **** C ********************************************************* A = 6.372e6 pi = 4.*atan(1.) dlon = 2*pi/ im dlat = pi/(jm-1) C11 = 1.0 / (4.0*A*IM*(1.0-COS(0.5*dlat))) CX1 = 1.0 / (4.0*A*dlon) CY1 = 1.0 / (4.0*A*dlat) do j=2,jm-1 phi = -pi/2.+(j-1)*dlat cosphi = cos(phi) do i=1,im cosij(i,j) = cosphi enddo enddo cosij(:,1) = 0.0 cosij(:,jm) = 0.0 C ******************************************************** C **** CALCULATE AVERAGE QUANTITIES **** C ********************************************************* DO j=2,jm-1 i =im DO ip1=1,im P1X(i,j) = ( U(ip1,j)+U(i,j) )*( DP(ip1,j)+DP(i,j) ) i =ip1 ENDDO ENDDO DO j=1,jm-1 DO I=1,im P1Y(I,j) = ( V(I,J+1)*COSIJ(I,J+1)+V(I,j)*COSIJ(I,j) )*( DP(I,J+1)+DP(I,J) ) ENDDO ENDDO C ********************************************************* C **** CALCULATE HORIZONTAL DIVERGENCE **** C ********************************************************* DO j=2,jm-1 im1=im DO i=1,im TMP1(i,j) = ( P1X(i,j)-P1X(im1,j) )*CX1 im1=i ENDDO DO I=1,im TMP2(I,j) = ( P1Y(I,j) -P1Y(I,j-1) )*CY1 DIV (I,j) = ( TMP1(I,j)+TMP2(I,j) )/(cosij(i,j)) ENDDO ENDDO C ********************************************************* C **** CALCULATE HORIZONTAL DIVERGENCE AT POLES **** C ********************************************************* DO 100 M=1,2 JPOLE = 1 + (M-1)*(jm-1) JPH = 1 + (M-1)*(jm-2) SUM11 = 0.0 DO I=1,im SUM11 = SUM11 + P1Y(I,JPH) ENDDO DO I=1,im DIV(I,JPOLE) = - MSGN(M) * C11*SUM11 ENDDO 100 CONTINUE RETURN END SUBROUTINE GRADQ (Q,DQDX,DQDY,IM,JM) C ********************************************************* C **** **** C **** THIS PROGRAM CALCULATES THE HORIZONTAL **** C **** GRADIENT OF THE INPUT FIELD Q **** C **** **** C **** ARGUMENTS: **** C **** Q ....... FIELD TO BE DIFFERENTIATED **** C **** DQDX .... LONGITUDINAL Q-DERIVATIVE **** C **** DQDY .... MERIDIONAL Q-DERIVATIVE **** C **** IM ...... NUMBER OF LONGITUDINAL POINTS **** C **** JM ...... NUMBER OF LATITUDINAL POINTS **** C **** **** C ********************************************************* use MAPL_ConstantsMod implicit none integer im,jm real Q(IM,JM) real DQDX(IM,JM) real DQDY(IM,JM) real Q1X(IM,JM) real Q2X(IM,JM) real Q1Y(IM,JM) real Q2Y(IM,JM) real acos(JM) real sinl(IM) real cosl(IM) real cx1,cx2,cy1,cy2,uc,vc,us,vs real dl,dp,a,pi,fjeq,phi integer i,j,m,ip1,ip2,jpole,msgn C ********************************************************* C **** INITIALIZATION **** C ********************************************************* a = MAPL_RADIUS pi = 4.0*atan(1.0) dl = 2.0*pi/im dp = pi/(jm-1) CX1 = 2.0 / ( 3.0*A*DL) CX2 = 1.0 / (12.0*A*DL) CY1 = 2.0 / ( 3.0*A*DP) CY2 = 1.0 / (12.0*A*DP) Q1X(:,:) = 0.0 Q2X(:,:) = 0.0 Q1Y(:,:) = 0.0 Q2Y(:,:) = 0.0 fjeq = ( jm+1 )*0.5 do j=2,jm-1 phi = dp * (j-fjeq) acos(j) = 1.0/( cos(phi) ) enddo do i=1,im/2 cosl(i) = -cos((i-1)*dl) cosl(i+im/2) = -cosl(i) sinl(i) = -sin((i-1)*dl) sinl(i+im/2) = -sinl(i) enddo C ********************************************************* C **** CALCULATE AVERAGE QUANTITIES **** C ********************************************************* do j = 2,jm-1 i = im-1 ip1 = im do ip2 = 1,im Q1X(i ,j) = Q(ip1,j) + Q(i,j) Q2X(ip1,j) = Q(ip2,j) + Q(i,j) i = ip1 ip1 = ip2 enddo enddo do j=1,jm-1 do i=1,im Q1Y(i,j) = Q(i,j+1) + Q(i,j) enddo enddo do j=2,jm-1 do i=1,im Q2Y(i,j) = Q(i,j+1) + Q(i,j-1) enddo enddo do i=1,im/2 Q2Y(i, 1) = Q(i,2) Q2Y(i,jm) = Q(i,jm-1) enddo do i=1,im/2 Q2Y(i , 1) = Q(i+im/2,2) + Q2Y(i,1) Q2Y(i+im/2, 1) = Q2Y(i,1) Q2Y(i ,jm) = Q(i+im/2,jm-1) + Q2Y(i,jm) Q2Y(i+im/2,jm) = Q2Y(i,jm) enddo C ********************************************************* C **** CALCULATE Q-GRADIENTS **** C ********************************************************* do j = 2,jm-1 i = im-1 ip1 = im do ip2 = 1,im DQDX(ip1,j) = ACOS(j) * ( ( Q1X(ip1,j)-Q1X(i,j) )*CX1 . - ( Q2X(ip2,j)-Q2X(i,j) )*CX2 ) i = ip1 ip1 = ip2 enddo enddo do j=2,jm-1 do i=1,im DQDY(i,j) = ( Q1Y(i,j) -Q1Y(i,j-1) )*CY1 . - ( Q2Y(i,j+1)-Q2Y(i,j-1) )*CY2 enddo enddo C ********************************************************* C **** CALCULATE Q-GRADIENTS (POLES) **** C ********************************************************* do i=1,im/2 Q1Y(i, 2) = Q(i, 1) + Q(i+im/2,2) Q1Y(i+im/2, 2) = Q(i+im/2, 1) + Q(i, 2) Q2Y(i, 1) = Q(i, 1) + Q(i+im/2,3) Q2Y(i+im/2, 1) = Q(i+im/2, 1) + Q(i, 3) Q1Y(i, jm) = Q(i, jm) + Q(i+im/2,jm-1) Q1Y(i+im/2,jm) = Q(i+im/2,jm) + Q(i, jm-1) Q2Y(i, jm) = Q(i, jm) + Q(i+im/2,jm-2) Q2Y(i+im/2,jm) = Q(i+im/2,jm) + Q(i, jm-2) enddo do i=1,im DQDY(i,jm) = ( Q1Y(i,jm)-Q1Y(i,jm-1) )*CY1 . - ( Q2Y(i,jm)-Q2Y(i,jm-1) )*CY2 DQDY(i, 1) = ( Q1Y(i,1)-Q1Y(i,2) )*CY1 . - ( Q2Y(i,2)-Q2Y(i,1) )*CY2 enddo C APPLY BOUNDARY CONDITIONS AT THE POLES C ========================================== DO 170 M=1,2 MSGN = (-1)**M JPOLE = 1 + (M-1)*(jm - 1) VC = 0.0 VS = 0.0 DO 180 I=1,IM VC = VC + DQDY(I,JPOLE)*COSL(I) VS = VS + DQDY(I,JPOLE)*SINL(I) 180 CONTINUE VC = 2.0 * VC / IM VS = 2.0 * VS / IM UC = - MSGN*VS US = MSGN*VC DO 190 I=1,IM DQDX(I,JPOLE) = US*SINL(I) + UC*COSL(I) DQDY(I,JPOLE) = VS*SINL(I) + VC*COSL(I) 190 CONTINUE 170 CONTINUE RETURN END SUBROUTINE VELPOT (DIV,VELP,im,jnp) integer IM,JNP real DIV(IM,JNP) real VELP(IM,JNP) real*8, allocatable :: VP(:,:) real*8, allocatable :: w(:) real*8, allocatable :: bdtf(:) real*8, allocatable :: bdts(:) real*8, allocatable :: bdps(:) real*8, allocatable :: bdpf(:) real*8 ts,tf,ps,pf,elmbda,pertrb,pi imp = im+1 iwk = 11*jnp+6*imp allocate ( vp(jnp,imp) ) allocate ( w(iwk) ) allocate ( bdtf(imp) ) allocate ( bdts(imp) ) allocate ( bdps(jnp) ) allocate ( bdpf(jnp) ) vp(:,:)=0.0 w(:)=0.0 bdtf(:)=0.0 bdts(:)=0.0 bdps(:)=0.0 bdpf(:)=0.0 c Transpose the input array c ------------------------- do j=1,jnp do i=1,im vp(j,i) = div(i,j) enddo vp(j,imp) = vp(j,1) enddo C === SET THE INPUT VARIABLES RAD = 6371000.0 PI = 3.14159265358979D0 INTL=0 TS=0.0 TF=PI M=JNP-1 MBDCND=9 PS=0.0 PF=2*PI N=IM NBDCND=0 ELMBDA=0 PERTRB=0 IDIMF=M+1 CALL PWSSSP (INTL,TS,TF,M,MBDCND,BDTS,BDTF,PS,PF,N,NBDCND,BDPS, * BDPF,ELMBDA,VP,IDIMF,PERTRB,IERROR,W) if( ierror.ne.0 ) then print *, 'PWSSSP IERROR = ',ierror stop endif c Scale by earth radius c --------------------- do j=1,jnp do i=1,im VELP(I,J) = VP(J,I) * RAD * RAD enddo enddo c Remove global mean c ------------------ CALL ZEROG (VELP,IM,JNP) deallocate ( vp ) deallocate ( w ) deallocate ( bdtf ) deallocate ( bdts ) deallocate ( bdps ) deallocate ( bdpf ) RETURN END SUBROUTINE ZEROG (VEL,IM,JNP) integer IM,JNP real VEL(IM,JNP) pi = 4.0*atan(1.0) dl = 2*pi/im dp = pi/(jnp-1) cap = 1-cos(0.5*dp) c Ensure unique pole values c ------------------------- sum1 = 0.0 sum2 = 0.0 do i=1,im sum1 = sum1 + vel(i,1) sum2 = sum2 + vel(i,jnp) enddo sum1 = sum1/im sum2 = sum2/im do i=1,im vel(i,1) = sum1 vel(i,jnp) = sum2 enddo c Compute global average c ---------------------- sum1 = 0.0 sum2 = 0.0 do i=1,im sum1 = sum1 + cap*vel(i,1) sum2 = sum2 + cap enddo do j=2,jnp-1 cosj = cos( -pi/2 + (j-1)*dp ) do i=1,im sum1 = sum1 + cosj*dp*vel(i,j) sum2 = sum2 + cosj*dp enddo enddo do i=1,im sum1 = sum1 + cap*vel(i,jnp) sum2 = sum2 + cap enddo qave = sum1/sum2 do j=1,jnp do i=1,im vel(i,j) = vel(i,j)-qave enddo enddo c print *, 'Remove Global Average: ', qave RETURN END