! +-======-+ ! Copyright (c) 2003-2018 United States Government as represented by ! the Admistrator of the National Aeronautics and Space Administration. ! All Rights Reserved. ! ! THIS OPEN SOURCE AGREEMENT ("AGREEMENT") DEFINES THE RIGHTS OF USE, ! REPRODUCTION, DISTRIBUTION, MODIFICATION AND REDISTRIBUTION OF CERTAIN ! COMPUTER SOFTWARE ORIGINALLY RELEASED BY THE UNITED STATES GOVERNMENT AS ! REPRESENTED BY THE GOVERNMENT AGENCY LISTED BELOW ("GOVERNMENT AGENCY"). ! THE UNITED STATES GOVERNMENT, AS REPRESENTED BY GOVERNMENT AGENCY, IS AN ! INTENDED THIRD-PARTY BENEFICIARY OF ALL SUBSEQUENT DISTRIBUTIONS OR ! REDISTRIBUTIONS OF THE SUBJECT SOFTWARE. ANYONE WHO USES, REPRODUCES, ! DISTRIBUTES, MODIFIES OR REDISTRIBUTES THE SUBJECT SOFTWARE, AS DEFINED ! HEREIN, OR ANY PART THEREOF, IS, BY THAT ACTION, ACCEPTING IN FULL THE ! RESPONSIBILITIES AND OBLIGATIONS CONTAINED IN THIS AGREEMENT. ! ! Government Agency: National Aeronautics and Space Administration ! Government Agency Original Software Designation: GSC-15354-1 ! Government Agency Original Software Title: GEOS-5 GCM Modeling Software ! User Registration Requested. Please Visit http://opensource.gsfc.nasa.gov ! Government Agency Point of Contact for Original Software: ! Dale Hithon, SRA Assistant, (301) 286-2691 ! ! +-======-+ #include "MAPL_Generic.h" subroutine windfix ( ua,va,plea, & ub,vb,pleb,im,jm,lm,VM,GRIDana,method, & vintdiva,vintdivb,vintdivc ) use ESMF use MAPL_Mod use MAPL_CommsMod implicit none type (ESMF_VM) :: VM type (ESMF_Grid) :: GRIDana integer :: DIMS(ESMF_MAXGRIDDIM) integer :: RC,STATUS character(len=ESMF_MAXSTR) :: IAm 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) real vintdiva(im,jm) real vintdivb(im,jm) real vintdivc(im,jm) real, pointer :: uglo(:,:,:) => null() real, pointer :: vglo(:,:,:) => null() real, pointer :: dglo(:,:,:) => null() real, pointer :: dpglo(:,:,:) => null() real, allocatable :: sum1(:,:) real, allocatable :: sum2(:,:) real, allocatable :: sum3(:,:) real, allocatable :: lambda(:,:) real, allocatable :: chi(:,:) real, allocatable :: uchi(:,:) real, allocatable :: vchi(:,:) real, allocatable :: dumdiv(:,:) real, allocatable :: dumchi(:,:) integer L, comm, myid, npes integer img, jmg, imjmg integer msgn integer k, kmax real undef Iam = "WINDFIX" msgn = 0 ! Scalar Flag undef = 1e15 call ESMF_VMGet (VM, localpet=myid, petcount=NPES, RC=STATUS) VERIFY_(STATUS) call ESMF_VmGet (VM, mpicommunicator=comm, rc=status) VERIFY_(STATUS) call MAPL_GridGet(GRIDana, globalCellCountPerDim=DIMS, RC=STATUS) img = DIMS(1) ! global grid dim jmg = DIMS(2) ! global grid dim imjmg = img*jmg 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 ! --------------------------- call MAPL_CollectiveGather3D(gridAna,ub,uglo,rc=status) VERIFY_(status) call MAPL_CollectiveGather3D(gridAna,vb,vglo,rc=status) VERIFY_(status) call MAPL_CollectiveGather3D(gridAna,dpb,dpglo,rc=status) VERIFY_(status) if (size(dpglo)>1) then allocate(dglo(img,jmg,size(dpglo,3)),stat=status) VERIFY_(status) end if ! Compute Vorticity and Divergence ! -------------------------------- if (size(uglo)>1) then do l=1,size(uglo,3) call getdiv (uglo(1,1,l),vglo(1,1,l),dpglo(1,1,l),dglo(1,1,l),img,jmg ) enddo end if call MAPL_CollectiveScatter3d(gridAna,dglo,divb,rc=status) VERIFY_(status) deallocate(uglo,vglo,dpglo) nullify(uglo,vglo,dpglo) if (associated(dglo)) then deallocate(dglo) nullify(dglo) end if ! Gather Winds for Analysis ! ------------------------- call MAPL_CollectiveGather3D(gridAna,ua,uglo,rc=status) VERIFY_(status) call MAPL_CollectiveGather3D(gridAna,va,vglo,rc=status) VERIFY_(status) call MAPL_CollectiveGather3D(gridAna,dpa,dpglo,rc=status) VERIFY_(status) if (size(dpglo)>1) then allocate(dglo(img,jmg,size(dpglo,3)),stat=status) VERIFY_(status) end if ! Compute Vorticity and Divergence ! -------------------------------- if (size(dpglo)>1) then do l=1,size(uglo,3) call getdiv (uglo(1,1,l),vglo(1,1,l),dpglo(1,1,l),dglo(1,1,l),img,jmg ) enddo end if call MAPL_CollectiveScatter3d(gridAna,dglo,diva,rc=status) VERIFY_(status) if (associated(dglo)) then deallocate(dglo) nullify(dglo) end if ! Compute Divergence Increment (to force vanishing vertical integral) ! ------------------------------------------------------------------- ! Minimize Relative Change ! ------------------------ 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 ! Minimize Absolute Change ! ------------------------ 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 vintdiva = sum1 vintdivb = sum2 vintdivc = sum3 ! call writit ( sum1,im,jm,1,55,lattice ) ! call writit ( sum2,im,jm,1,55,lattice ) ! call writit ( sum3,im,jm,1,55,lattice ) ! Gather and Broadcast Divergence Increment ! ----------------------------------------- call MAPL_CollectiveGather3D(gridAna,diva,dglo,rc=status) VERIFY_(status) ! Compute Wind Increments Associated with Divergence Increment ! ------------------------------------------------------------ if (size(dpglo)>1) then do l=1,size(dglo,3) call VELPOT_SP (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) end do end if ! Scatter Winds ! ------------- call MAPL_CollectiveScatter3d(gridAna,uglo,ua,rc=status) VERIFY_(status) call MAPL_CollectiveScatter3d(gridAna,vglo,va,rc=status) VERIFY_(status) if (associated(dglo)) then deallocate(dglo,uglo,vglo,dpglo) nullify(dglo,uglo,vglo,dpglo) end if deallocate ( sum1,sum2,sum3,lambda ) deallocate ( chi,uchi,vchi ) return end SUBROUTINE GETDIV ( U,V,DP,DIV,IM,JM ) ! ******************************************************************** ! **** **** ! **** THIS PROGRAM CALCULATES DIVERGENCE **** ! **** AT EACH LEVEL FOR A NON-STAGGERED A-GRID **** ! **** **** ! **** INPUT: **** ! **** U ....... ZONAL WIND **** ! **** V ....... MERIDIONAL WIND **** ! **** IM ...... NUMBER OF LONGITUDE POINTS **** ! **** JM ...... NUMBER OF LATITUDE POINTS **** ! **** **** ! **** OUTPUT: **** ! **** DIV (IM,JM) .... DIVERGENCE **** ! **** **** ! ******************************************************************** 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/ ! ********************************************************* ! **** INITIALIZATION FOR DIVERGENCE **** ! ********************************************************* 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 ! ******************************************************** ! **** CALCULATE AVERAGE QUANTITIES **** ! ********************************************************* 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 ! ********************************************************* ! **** CALCULATE HORIZONTAL DIVERGENCE **** ! ********************************************************* 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 ! ********************************************************* ! **** CALCULATE HORIZONTAL DIVERGENCE AT POLES **** ! ********************************************************* 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) ! ********************************************************* ! **** **** ! **** THIS PROGRAM CALCULATES THE HORIZONTAL **** ! **** GRADIENT OF THE INPUT FIELD Q **** ! **** **** ! **** ARGUMENTS: **** ! **** Q ....... FIELD TO BE DIFFERENTIATED **** ! **** DQDX .... LONGITUDINAL Q-DERIVATIVE **** ! **** DQDY .... MERIDIONAL Q-DERIVATIVE **** ! **** IM ...... NUMBER OF LONGITUDINAL POINTS **** ! **** JM ...... NUMBER OF LATITUDINAL POINTS **** ! **** **** ! ********************************************************* 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 ! ********************************************************* ! **** INITIALIZATION **** ! ********************************************************* 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 ! ********************************************************* ! **** CALCULATE AVERAGE QUANTITIES **** ! ********************************************************* 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 ! ********************************************************* ! **** CALCULATE Q-GRADIENTS **** ! ********************************************************* 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 ! ********************************************************* ! **** CALCULATE Q-GRADIENTS (POLES) **** ! ********************************************************* 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 ! APPLY BOUNDARY CONDITIONS AT THE POLES ! ========================================== 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_SP (DIV,VELP,im,jnp) implicit none integer IM,JNP real DIV(IM,JNP) real VELP(IM,JNP) real*8, allocatable :: VP(:,:) real*8, allocatable :: even(:,:) real*8, allocatable :: odd(:,:) real*8 :: xlmbda(1) real*8 :: pertrb(1) real*8, allocatable :: wsha(:) real*8, allocatable :: wshs(:) real*8, allocatable :: work(:) integer :: i,j,imp integer :: mtrunc,ntrunc integer :: lwsha, lwshs, lwork integer :: l1, l2, maxnt integer :: isym, nt, ierror real*8 PI real*8 RAD imp = im+1 mtrunc = (im+2)/2 ntrunc = jnp l1 = jnp l2 = (jnp+1)/2 maxnt = 1 lwsha = 2*jnp*l2+3*((l1-2)*(jnp+jnp-l1-1))/2+imp+15 lwshs = 2*jnp*l2+3*((l1-2)*(jnp+jnp-l1-1))/2+imp+15 lwork = jnp*(2*maxnt*imp+max0(6*l2,imp)+2*maxnt*l1+1) allocate ( vp(jnp,im) ) allocate ( even(mtrunc,ntrunc) ) allocate ( odd(mtrunc,ntrunc) ) allocate ( wsha(lwsha) ) allocate ( wshs(lwshs) ) allocate ( work(lwork) ) ! Transpose the input array ! ------------------------- do j=1,jnp do i=1,im vp(j,i) = div(i,j) enddo enddo ! do j=1,jnp ! do i=1,im/2 ! vp(j,i) = div(i+im/2,jnp-j+1) ! enddo ! do i=im/2+1,im ! vp(j,i) = div(i-im/2,jnp-j+1) ! enddo ! enddo ! === SET THE INPUT VARIABLES RAD = 6371000.0 PI = 3.14159265358979D0 xlmbda(1) = 0.0 isym = 0 nt = 1 call SHAECI(jnp,im,wsha,lwsha,work,lwork,ierror) call SHSECI(jnp,im,wshs,lwshs,work,lwork,ierror) call SHAEC(jnp,im,isym,nt,vp,jnp,imp,even,odd,mtrunc,ntrunc,wsha,lwsha,work,lwork,ierror) call ISLAPEC(jnp,im,isym,nt,xlmbda,vp,jnp,imp,even,odd,mtrunc,ntrunc,wshs,lwshs,work,lwork,pertrb,ierror) ! Scale by earth radius ! --------------------- do j=1,jnp do i=1,im VELP(I,J) = VP(J,I) * RAD * RAD enddo enddo ! do j=1,jnp ! do i=1,im/2 ! VELP(I+im/2,jnp-J+1) = VP(J,I) * RAD * RAD ! enddo ! do i=im/2+1,im ! VELP(I-im/2,jnp-J+1) = VP(J,I) * RAD * RAD ! enddo ! enddo ! Remove global mean ! ------------------ CALL ZEROG (VELP,IM,JNP) deallocate ( vp ) deallocate ( even ) deallocate ( odd ) deallocate ( wshs ) deallocate ( wsha ) deallocate ( work ) 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 ! Transpose the input array ! ------------------------- do j=1,jnp do i=1,im vp(j,i) = div(i,j) enddo vp(j,imp) = vp(j,1) enddo ! === 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 ! Scale by earth radius ! --------------------- do j=1,jnp do i=1,im VELP(I,J) = VP(J,I) * RAD * RAD enddo enddo ! Remove global mean ! ------------------ 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) ! Ensure unique pole values ! ------------------------- 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 ! Compute global average ! ---------------------- 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 ! print *, 'Remove Global Average: ', qave RETURN END subroutine writit ( q,im,jm,lm,ku,lattice ) use G3_MPI_Util_Mod implicit none type ( dynamics_lattice_type ) lattice integer im,jm,lm,L,ku,img,jmg real q(im,jm,lm) real, allocatable :: glo(:,:) real*4, allocatable :: a(:,:) img = lattice%imglobal jmg = lattice%jmglobal allocate ( glo(img,jmg) ) allocate ( a(img,jmg) ) do L=1,lm call G3_GATHER ( glo,q(:,:,L),lattice ) if( lattice%myid.eq.0 ) then a = glo write(ku) a endif enddo deallocate ( glo ) deallocate ( a ) return end subroutine regrid ( qin,im_in,jm_in,qout,im_out,jm_out,undef,msgn ) implicit none integer im_in ,jm_in ,msgn integer im_out,jm_out real undef real qin(im_in ,jm_in ) real qout(im_out,jm_out) real q10x10(360*6,180*6) ! Parse Arbitray Field (im,jm) to 10'x10' Variable ! ------------------------------------------------ call bin_10x10 ( qin,im_in,jm_in,q10x10 ) ! Bin 10'x10' Variable to Output Field (im_out,jm_out) ! ---------------------------------------------------- call averaged_10x10 ( q10x10,qout,im_out,jm_out,undef,msgn ) return end subroutine averaged_10x10 ( z10x10,z,im,jm,undef,msgn ) !*********************************************************************** ! ! PURPOSE: ! ======== ! Average a (10m X 10m) input array to an output array (im,jm) ! ! INPUT: ! ====== ! z10x10 ..... Input array(360*6,180*6) ! msgn ....... Integer Flag for scalar (0) or vector (1) ! ! OUTPUT: ! ======= ! z .......... Output array(im,jm) ! im ......... Longitudinal dimension of z ! jm ......... Latitudinal dimension of z ! ! NOTES: ! ====== ! Input array z10x10 represents values within a 10min X 10min grid-box. ! Each box is referenced by the latitude and longitude of ! its southwest corner, not its center point. Thus, ! the height associated with a coordinate actually ! represents the heights centered to the northeast of that point. ! ! Output array z(im,jm) is assumed to be on an A-grid. ! z(i,j) represents the value at the center of the grid-box. ! z(1,j) is located at lon=-180. ! z(i,1) is located at lat=-90. ! z(i,jm) is located at lat=+90. ! !*********************************************************************** !* GODDARD LABORATORY FOR ATMOSPHERES * !*********************************************************************** implicit none integer im,jm,msgn real z(im,jm) real dlam(im), dphi(jm) real z10x10(360*6,180*6) integer i,j,ibeg,iend,jbeg,jend integer ii,jj,itmp real sum1,sum2 real zlat,zlon real lon1,lon2,wx real lat1,lat2,wy real lonbeg,lonend,lat,coslat real latbeg,latend real undef real pi,dz real lon_cmp(im) real lat_cmp(jm) logical defined pi = 4.*atan(1.) dz = pi/(6.*180) dlam = 2*pi/ im dphi = pi/(jm-1) ! Compute Computational Lambda's and Phi's ! ---------------------------------------- lon_cmp(1) = -pi do i=2,im lon_cmp(i) = lon_cmp(i-1) + dlam(i-1) enddo lat_cmp(1) = -pi*0.5 do j=2,jm-1 lat_cmp(j) = lat_cmp(j-1) + dphi(j-1) enddo lat_cmp(jm) = pi*0.5 ! Compute average away from poles ! ------------------------------- do j=2,jm-1 do i=1,im zlat = lat_cmp(j) zlon = lon_cmp(i) latbeg = zlat-dphi(j-1)/2 latend = zlat+dphi(j) /2 if( i.eq.1 ) then lonbeg = zlon-dlam(im) /2 else lonbeg = zlon-dlam(i-1)/2 endif lonend = zlon+dlam(i) /2 ibeg = 1.+(lonbeg+pi) /dz iend = 1.+(lonend+pi) /dz jbeg = 1.+(latbeg+pi/2)/dz jend = 1.+(latend+pi/2)/dz sum1 = 0 sum2 = 0 do jj=jbeg,jend lat = -pi/2+(jj-0.5)*dz coslat = cos(lat) lat1 = -pi/2 + (jj-1)*dz lat2 = -pi/2 + jj *dz wy = 1.0 if( lat1.lt.latbeg ) wy = (lat2-latbeg)/dz if( lat2.gt.latend ) wy = (latend-lat1)/dz if(ibeg.ge.1) then do ii=ibeg,iend if( defined(z10x10(ii,jj),undef) ) then lon1 = -pi + (ii-1)*dz lon2 = -pi + ii *dz wx = 1.0 if( lon1.lt.lonbeg ) wx = (lon2-lonbeg)/dz if( lon2.gt.lonend ) wx = (lonend-lon1)/dz sum1 = sum1 + z10x10(ii,jj)*coslat*wx*wy sum2 = sum2 + coslat*wx*wy endif enddo else itmp = 1.+(lonbeg+0.1*dz+3*pi)/dz do ii=itmp,360*6 if( defined(z10x10(ii,jj),undef) ) then lon1 = -pi + (ii-1)*dz lon2 = -pi + ii *dz wx = 1.0 if( lon1.lt.lonbeg+2*pi ) wx = (lon2-lonbeg-2*pi)/dz if( lon2.gt.lonend+2*pi ) wx = (2*pi+lonend-lon1)/dz sum1 = sum1 + z10x10(ii,jj)*coslat*wx*wy sum2 = sum2 + coslat*wx*wy endif enddo do ii=1,iend if( defined(z10x10(ii,jj),undef) ) then lon1 = -pi + (ii-1)*dz lon2 = -pi + ii *dz wx = 1.0 if( lon1.lt.lonbeg ) wx = (lon2-lonbeg)/dz if( lon2.gt.lonend ) wx = (lonend-lon1)/dz sum1 = sum1 + z10x10(ii,jj)*coslat*wx*wy sum2 = sum2 + coslat*wx*wy endif enddo endif enddo if( sum2.ne.0.0 ) then z(i,j) = sum1/sum2 else z(i,j) = undef endif enddo enddo ! Compute average at South Pole ! ----------------------------- j=1 do i=1,im zlat = lat_cmp(j) zlon = lon_cmp(i) latbeg = zlat latend = zlat+dphi(j) /2 if( i.eq.1 ) then lonbeg = zlon-dlam(im) /2 else lonbeg = zlon-dlam(i-1)/2 endif lonend = zlon+dlam(i) /2 ibeg = 1.+(lonbeg+pi) /dz iend = 1.+(lonend+pi) /dz jbeg = 1 jend = 1.+(latend+pi/2)/dz sum1 = 0 sum2 = 0 do jj=jbeg,jend lat = -pi/2+(jj-0.5)*dz coslat = cos(lat) lat1 = -pi/2 + (jj-1)*dz lat2 = -pi/2 + jj *dz wy = 1.0 if( lat1.lt.latbeg ) wy = (lat2-latbeg)/dz if( lat2.gt.latend ) wy = (latend-lat1)/dz if(ibeg.ge.1) then do ii=ibeg,iend if( defined(z10x10(ii,jj),undef) ) then lon1 = -pi + (ii-1)*dz lon2 = -pi + ii *dz wx = 1.0 if( lon1.lt.lonbeg ) wx = (lon2-lonbeg)/dz if( lon2.gt.lonend ) wx = (lonend-lon1)/dz sum1 = sum1 + z10x10(ii,jj)*coslat*wx*wy sum2 = sum2 + coslat*wx*wy endif enddo else itmp = 1.+(lonbeg+0.1*dz+3*pi)/dz do ii=itmp,360*6 if( defined(z10x10(ii,jj),undef) ) then lon1 = -pi + (ii-1)*dz lon2 = -pi + ii *dz wx = 1.0 if( lon1.lt.lonbeg+2*pi ) wx = (lon2-lonbeg-2*pi)/dz if( lon2.gt.lonend+2*pi ) wx = (2*pi+lonend-lon1)/dz sum1 = sum1 + z10x10(ii,jj)*coslat*wx*wy sum2 = sum2 + coslat*wx*wy endif enddo do ii=1,iend if( defined(z10x10(ii,jj),undef) ) then lon1 = -pi + (ii-1)*dz lon2 = -pi + ii *dz wx = 1.0 if( lon1.lt.lonbeg ) wx = (lon2-lonbeg)/dz if( lon2.gt.lonend ) wx = (lonend-lon1)/dz sum1 = sum1 + z10x10(ii,jj)*coslat*wx*wy sum2 = sum2 + coslat*wx*wy endif enddo endif enddo if( sum2.ne.0.0 ) then z(i,j) = sum1/sum2 else z(i,j) = undef endif enddo ! Compute average at North Pole ! ----------------------------- j=jm do i=1,im zlat = lat_cmp(j) zlon = lon_cmp(i) latbeg = zlat-dphi(j-1)/2 latend = zlat if( i.eq.1 ) then lonbeg = zlon-dlam(im) /2 else lonbeg = zlon-dlam(i-1)/2 endif lonend = zlon+dlam(i) /2 ibeg = 1.+(lonbeg+pi) /dz iend = 1.+(lonend+pi) /dz jbeg = 1.+(latbeg+pi/2)/dz jend = 1080 sum1 = 0 sum2 = 0 do jj=jbeg,jend lat = -pi/2+(jj-0.5)*dz coslat = cos(lat) lat1 = -pi/2 + (jj-1)*dz lat2 = -pi/2 + jj *dz wy = 1.0 if( lat1.lt.latbeg ) wy = (lat2-latbeg)/dz if( lat2.gt.latend ) wy = (latend-lat1)/dz if(ibeg.ge.1) then do ii=ibeg,iend if( defined(z10x10(ii,jj),undef) ) then lon1 = -pi + (ii-1)*dz lon2 = -pi + ii *dz wx = 1.0 if( lon1.lt.lonbeg ) wx = (lon2-lonbeg)/dz if( lon2.gt.lonend ) wx = (lonend-lon1)/dz sum1 = sum1 + z10x10(ii,jj)*coslat*wx*wy sum2 = sum2 + coslat*wx*wy endif enddo else itmp = 1.+(lonbeg+0.1*dz+3*pi)/dz do ii=itmp,360*6 if( defined(z10x10(ii,jj),undef) ) then lon1 = -pi + (ii-1)*dz lon2 = -pi + ii *dz wx = 1.0 if( lon1.lt.lonbeg+2*pi ) wx = (lon2-lonbeg-2*pi)/dz if( lon2.gt.lonend+2*pi ) wx = (2*pi+lonend-lon1)/dz sum1 = sum1 + z10x10(ii,jj)*coslat*wx*wy sum2 = sum2 + coslat*wx*wy endif enddo do ii=1,iend if( defined(z10x10(ii,jj),undef) ) then lon1 = -pi + (ii-1)*dz lon2 = -pi + ii *dz wx = 1.0 if( lon1.lt.lonbeg ) wx = (lon2-lonbeg)/dz if( lon2.gt.lonend ) wx = (lonend-lon1)/dz sum1 = sum1 + z10x10(ii,jj)*coslat*wx*wy sum2 = sum2 + coslat*wx*wy endif enddo endif enddo if( sum2.ne.0.0 ) then z(i,j) = sum1/sum2 else z(i,j) = undef endif enddo ! Average Pole Values ! ------------------- if( msgn.eq.0 ) then sum1 = 0 j = 0 do i=1,im if( defined(z(i,1),undef) ) then sum1 = sum1 + z(i,1) j = j + 1 endif enddo if( j.ne.0 ) then z(:,1) = sum1/j else z(:,1) = undef endif sum2 = 0 j = 0 do i=1,im if( defined(z(i,jm),undef) ) then sum2 = sum2 + z(i,jm) j = j + 1 endif enddo if( j.ne.0 ) then z(:,jm) = sum2/j else z(:,jm) = undef endif endif return end subroutine bin_10x10 ( z,im,jm,z10x10 ) !*********************************************************************** ! ! PURPOSE: ! ======== ! Compute a (10m X 10m) array binned from an input array (im,jm) ! ! INPUT: ! ====== ! z .......... Input array(im,jm) ! im ......... Longitudinal dimension of z ! jm ......... Latitudinal dimension of z ! ! OUTPUT: ! ======= ! z10x10 ..... Output array(360*6,180*6) ! ! NOTES: ! ====== ! Input array z(im,jm) is assumed to be on an A-grid. ! z(i,j) represents the value at the center of the grid-box. ! z(1,j) is located at lon=-180. ! z(i,1) is located at lat=-90. ! z(i,jm) is located at lat=+90. ! ! Output array z10x10 represents values within a 10min X 10min grid-box. ! Each box is referenced by the latitude and longitude of ! its southwest corner, not its center point. Thus, ! the height associated with a coordinate actually ! represents the heights centered to the northeast of that point. ! !*********************************************************************** !* GODDARD LABORATORY FOR ATMOSPHERES * !*********************************************************************** implicit none integer im,jm real z(im,jm) real z10x10(360*6,180*6) integer i,j,ii,jj,ibeg,iend,jbeg,jend real zlatc,zlonc real lonbeg,lonend,lat real latbeg,latend real pi,dl,dp,dz pi = 4.*atan(1.) dl = 2*pi/im dp = pi/(jm-1) dz = pi/(6.*180) do j=1,180*6 do i=1,360*6 zlatc = -pi/2+(j-0.5)*dz ! Latitude at center of 10x10 box zlonc = -pi +(i-0.5)*dz ! Longitude at center of 10x10 box ! Find bounding lat and lon on IMxJM grid ! --------------------------------------- iend = nint( 1.+(zlonc+pi)/dl ) lonend = -pi + (iend-1)*dl if( lonend.ge.zlonc ) then lonbeg = -pi + (iend-2)*dl else iend = iend+1 lonbeg = lonend lonend = -pi + (iend-1)*dl endif ibeg = iend-1 jend = nint( 1.+(zlatc+pi/2)/dp ) latend = -pi/2 + (jend-1)*dp if( latend.ge.zlatc ) then latbeg = -pi/2 + (jend-2)*dp else jend = jend+1 latbeg = latend latend = -pi/2 + (jend-1)*dp endif jbeg = jend-1 if(iend.gt.im) iend=iend-im if( zlonc.le.lonbeg+0.5*dl ) then ii = ibeg else ii = iend endif if( zlatc.le.latbeg+0.5*dp ) then jj = jbeg else jj = jend endif z10x10(i,j) = z(ii,jj) enddo enddo return end function defined ( q,undef ) implicit none logical defined real q,undef defined = abs(q-undef).gt.0.1*abs(undef) return end