! $Id: ppm3d.f90 5159 2024-08-02 19:58:25Z abarral $ !From lin@explorer.gsfc.nasa.gov Wed Apr 15 17:44:44 1998 !Date: Wed, 15 Apr 1998 11:37:03 -0400 !From: lin@explorer.gsfc.nasa.gov !To: Frederic.Hourdin@lmd.jussieu.fr !Subject: 3D transport module of the GSFC CTM and GEOS GCM !This code is sent to you by S-J Lin, DAO, NASA-GSFC !Note: this version is intended for machines like CRAY !-90. No multitasking directives implemented. ! ******************************************************************** ! TransPort Core for Goddard Chemistry Transport Model (G-CTM), Goddard ! Earth Observing System General Circulation Model (GEOS-GCM), and Data ! Assimilation System (GEOS-DAS). ! ******************************************************************** ! Purpose: given horizontal winds on a hybrid sigma-p surfaces, ! one CALL to tpcore updates the 3-D mixing ratio ! fields one time step (NDT). [vertical mass flux is computed ! internally consistent with the discretized hydrostatic mass ! continuity equation of the C-Grid GEOS-GCM (for IGD=1)]. ! Schemes: Multi-dimensional Flux Form Semi-Lagrangian (FFSL) scheme based ! on the van Leer or PPM. ! (see Lin and Rood 1996). ! Version 4.5 ! Last modified: Dec. 5, 1996 ! Major changes from version 4.0: a more general vertical hybrid sigma- ! pressure coordinate. ! Subroutines modified: xtp, ytp, fzppm, qckxyz ! Subroutines deleted: vanz ! Author: Shian-Jiann Lin ! mail address: ! Shian-Jiann Lin* ! Code 910.3, NASA/GSFC, Greenbelt, MD 20771 ! Phone: 301-286-9540 ! E-mail: lin@dao.gsfc.nasa.gov ! *affiliation: ! Joint Center for Earth Systems Technology ! The University of Maryland Baltimore County ! NASA - Goddard Space Flight Center ! References: ! 1. Lin, S.-J., and R. B. Rood, 1996: Multidimensional flux form semi- ! Lagrangian transport schemes. Mon. Wea. Rev., 124, 2046-2070. ! 2. Lin, S.-J., W. C. Chao, Y. C. Sud, and G. K. Walker, 1994: A class of ! the van Leer-type transport schemes and its applications to the moist- ! ure transport in a General Circulation Model. Mon. Wea. Rev., 122, ! 1575-1593. ! ****6***0*********0*********0*********0*********0*********0**********72 SUBROUTINE ppm3d(IGD,Q,PS1,PS2,U,V,W,NDT,IORD,JORD,KORD,NC,IMR, & JNP,j1,NLAY,AP,BP,PT,AE,fill,dum,Umax) IMPLICIT NONE ! rajout de déclarations ! integer Jmax,kmax,ndt0,nstep,k,j,i,ic,l,js,jn,imh,iad,jad,krd ! integer iu,iiu,j2,jmr,js0,jt ! real dtdy,dtdy5,rcap,iml,jn0,imjm,pi,dl,dp ! real dt,cr1,maxdt,ztc,d5,sum1,sum2,ru ! ******************************************************************** ! ============= ! INPUT: ! ============= ! Q(IMR,JNP,NLAY,NC): mixing ratios at current time (t) ! NC: total # of constituents ! IMR: first dimension (E-W); # of Grid intervals in E-W is IMR ! JNP: 2nd dimension (N-S); # of Grid intervals in N-S is JNP-1 ! NLAY: 3rd dimension (# of layers); vertical index increases from 1 at ! the model top to NLAY near the surface (see fig. below). ! It is assumed that 6 <= NLAY <= JNP (for dynamic memory allocation) ! PS1(IMR,JNP): surface pressure at current time (t) ! PS2(IMR,JNP): surface pressure at mid-time-level (t+NDT/2) ! PS2 is replaced by the predicted PS (at t+NDT) on output. ! Note: surface pressure can have any unit or can be multiplied by any ! const. ! The pressure at layer edges are defined as follows: ! p(i,j,k) = AP(k)*PT + BP(k)*PS(i,j) (1) ! Where PT is a constant having the same unit as PS. ! AP and BP are unitless constants given at layer edges ! defining the vertical coordinate. ! BP(1) = 0., BP(NLAY+1) = 1. ! The pressure at the model top is PTOP = AP(1)*PT ! For pure sigma system set AP(k) = 1 for all k, PT = PTOP, ! BP(k) = sige(k) (sigma at edges), PS = Psfc - PTOP. ! Note: the sigma-P coordinate is a subset of Eq. 1, which in turn ! is a subset of the following even more general sigma-P-thelta coord. ! currently under development. ! p(i,j,k) = (AP(k)*PT + BP(k)*PS(i,j))/(D(k)-C(k)*TE**(-1/kapa)) ! ///////////////////////////////// ! / \ ------------- PTOP -------------- AP(1), BP(1) ! | ! delp(1) | ........... Q(i,j,1) ............ ! | ! W(1) \ / --------------------------------- AP(2), BP(2) ! W(k-1) / \ --------------------------------- AP(k), BP(k) ! | ! delp(K) | ........... Q(i,j,k) ............ ! | ! W(k) \ / --------------------------------- AP(k+1), BP(k+1) ! / \ --------------------------------- AP(NLAY), BP(NLAY) ! | ! delp(NLAY) | ........... Q(i,j,NLAY) ......... ! | ! W(NLAY)=0 \ / ------------- surface ----------- AP(NLAY+1), BP(NLAY+1) ! ////////////////////////////////// ! U(IMR,JNP,NLAY) & V(IMR,JNP,NLAY):winds (m/s) at mid-time-level (t+NDT/2) ! U and V may need to be polar filtered in advance in some cases. ! IGD: grid type on which winds are defined. ! IGD = 0: A-Grid [all variables defined at the same point from south ! pole (j=1) to north pole (j=JNP) ] ! IGD = 1 GEOS-GCM C-Grid ! [North] ! V(i,j) ! | ! | ! | ! U(i-1,j)---Q(i,j)---U(i,j) [EAST] ! | ! | ! | ! V(i,j-1) ! U(i, 1) is defined at South Pole. ! V(i, 1) is half grid north of the South Pole. ! V(i,JMR) is half grid south of the North Pole. ! V must be defined at j=1 and j=JMR if IGD=1 ! V at JNP need not be given. ! NDT: time step in seconds (need not be constant during the course of ! the integration). Suggested value: 30 min. for 4x5, 15 min. for 2x2.5 ! (Lat-Lon) resolution. Smaller values are recommanded if the model ! has a well-resolved stratosphere. ! J1 defines the size of the polar cap: ! South polar cap edge is located at -90 + (j1-1.5)*180/(JNP-1) deg. ! North polar cap edge is located at 90 - (j1-1.5)*180/(JNP-1) deg. ! There are currently only two choices (j1=2 or 3). ! IMR must be an even integer if j1 = 2. Recommended value: J1=3. ! IORD, JORD, and KORD are integers controlling various options in E-W, N-S, ! and vertical transport, respectively. Recommended values for positive ! definite scalars: IORD=JORD=3, KORD=5. Use KORD=3 for non- ! positive definite scalars or when linear correlation between constituents ! is to be maintained. ! _ORD= ! 1: 1st order upstream scheme (too diffusive, not a useful option; it ! can be used for debugging purposes; this is THE only known "linear" ! monotonic advection scheme.). ! 2: 2nd order van Leer (full monotonicity constraint; ! see Lin et al 1994, MWR) ! 3: monotonic PPM* (slightly improved PPM of Collela & Woodward 1984) ! 4: semi-monotonic PPM (same as 3, but overshoots are allowed) ! 5: positive-definite PPM (constraint on the subgrid distribution is ! only strong enough to prevent generation of negative values; ! both overshoots & undershoots are possible). ! 6: un-constrained PPM (nearly diffusion free; slightly faster but ! positivity not quaranteed. Use this option only when the fields ! and winds are very smooth). ! *PPM: Piece-wise Parabolic Method ! Note that KORD <=2 options are no longer supported. DO not use option 4 or 5. ! for non-positive definite scalars (such as Ertel Potential Vorticity). ! The implicit numerical diffusion decreases as _ORD increases. ! The last two options (ORDER=5, 6) should only be used when there is ! significant explicit diffusion (such as a turbulence parameterization). You ! might get dispersive results otherwise. ! No filter of any kind is applied to the constituent fields here. ! AE: Radius of the sphere (meters). ! Recommended value for the planet earth: 6.371E6 ! fill(logical): flag to do filling for negatives (see note below). ! Umax: Estimate (upper limit) of the maximum U-wind speed (m/s). ! (220 m/s is a good value for troposphere model; 280 m/s otherwise) ! ============= ! Output ! ============= ! Q: mixing ratios at future time (t+NDT) (original values are over-written) ! W(NLAY): large-scale vertical mass flux as diagnosed from the hydrostatic ! relationship. W will have the same unit as PS1 and PS2 (eg, mb). ! W must be divided by NDT to get the correct mass-flux unit. ! The vertical Courant number C = W/delp_UPWIND, where delp_UPWIND ! is the pressure thickness in the "upwind" direction. For example, ! C(k) = W(k)/delp(k) if W(k) > 0; ! C(k) = W(k)/delp(k+1) if W(k) < 0. ! ( W > 0 is downward, ie, toward surface) ! PS2: predicted PS at t+NDT (original values are over-written) ! ******************************************************************** ! NOTES: ! This forward-in-time upstream-biased transport scheme reduces to ! the 2nd order center-in-time center-in-space mass continuity eqn. ! if Q = 1 (constant fields will remain constant). This also ensures ! that the computed vertical velocity to be identical to GEOS-1 GCM ! for on-line transport. ! A larger polar cap is used if j1=3 (recommended for C-Grid winds or when ! winds are noisy near poles). ! Flux-Form Semi-Lagrangian transport in the East-West direction is used ! when and where Courant # is greater than one. ! The user needs to change the parameter Jmax or Kmax if the resolution ! is greater than 0.5 deg in N-S or 150 layers in the vertical direction. ! (this TransPort Core is otherwise resolution independent and can be used ! as a library routine). ! PPM is 4th order accurate when grid spacing is uniform (x & y); 3rd ! order accurate for non-uniform grid (vertical sigma coord.). ! Time step is limitted only by transport in the meridional direction. ! (the FFSL scheme is not implemented in the meridional direction). ! Since only 1-D limiters are applied, negative values could ! potentially be generated when large time step is used and when the ! initial fields contain discontinuities. ! This does not necessarily imply the integration is unstable. ! These negatives are typically very small. A filling algorithm is ! activated if the user set "fill" to be true. ! The van Leer scheme used here is nearly as accurate as the original PPM ! due to the use of a 4th order accurate reference slope. The PPM imple- ! mented here is an improvement over the original and is also based on ! the 4th order reference slope. ! ****6***0*********0*********0*********0*********0*********0**********72 ! User modifiable parameters INTEGER,parameter :: Jmax = 361, kmax = 150 ! ****6***0*********0*********0*********0*********0*********0**********72 ! Input-Output arrays ! REAL :: Q(IMR,JNP,NLAY,NC),PS1(IMR,JNP),PS2(IMR,JNP), & U(IMR,JNP,NLAY),V(IMR,JNP,NLAY),AP(NLAY+1), & BP(NLAY+1),W(IMR,JNP,NLAY),NDT,val(NLAY),Umax INTEGER :: IGD,IORD,JORD,KORD,NC,IMR,JNP,j1,NLAY,AE INTEGER :: IMRD2 REAL :: PT LOGICAL :: cross, fill, dum ! Local dynamic arrays REAL :: CRX(IMR,JNP),CRY(IMR,JNP),xmass(IMR,JNP),ymass(IMR,JNP), & fx1(IMR+1),DPI(IMR,JNP,NLAY),delp1(IMR,JNP,NLAY), & WK1(IMR,JNP,NLAY),PU(IMR,JNP),PV(IMR,JNP),DC2(IMR,JNP), & delp2(IMR,JNP,NLAY),DQ(IMR,JNP,NLAY,NC),VA(IMR,JNP), & UA(IMR,JNP),qtmp(-IMR:2*IMR) ! Local static arrays REAL :: DTDX(Jmax), DTDX5(Jmax), acosp(Jmax), & cosp(Jmax), cose(Jmax), DAP(kmax),DBK(Kmax) data NDT0, NSTEP /0, 0/ data cross /.TRUE./ REAL :: DTDY, DTDY5, RCAP INTEGER :: JS0, JN0, IML, JMR, IMJM SAVE DTDY, DTDY5, RCAP, JS0, JN0, IML, & DTDX, DTDX5, ACOSP, COSP, COSE, DAP,DBK INTEGER :: NDT0, NSTEP, j2, k,j,i,ic,l,JS,JN,IMH INTEGER :: IU,IIU,JT,iad,jad,krd REAL :: r23,r3,PI,DL,DP,DT,CR1,MAXDT,ZTC,D5 REAL :: sum1,sum2,ru JMR = JNP -1 IMJM = IMR*JNP j2 = JNP - j1 + 1 NSTEP = NSTEP + 1 ! *********** Initialization ********************** IF(NSTEP==1) THEN WRITE(6,*) '------------------------------------ ' WRITE(6,*) 'NASA/GSFC Transport Core Version 4.5' WRITE(6,*) '------------------------------------ ' WRITE(6,*) 'IMR=',IMR,' JNP=',JNP,' NLAY=',NLAY,' j1=',j1 WRITE(6,*) 'NC=',NC,IORD,JORD,KORD,NDT ! controles sur les parametres IF(NLAY<6) THEN WRITE(6,*) 'NLAY must be >= 6' stop ENDIF IF (JNP= NLAY' stop ENDIF IMRD2=mod(IMR,2) IF (j1==2.AND.IMRD2/=0) THEN WRITE(6,*) 'if j1=2 IMR must be an even integer' stop ENDIF IF(Jmax=0.95) THEN JS0 = 0 JN0 = 0 IML = IMR-2 ZTC = 0. else ZTC = acos(CR1) * (180./PI) JS0 = REAL(JMR)*(90.-ZTC)/180. + 2 JS0 = max(JS0, J1+1) IML = min(6*JS0/(J1-1)+2, 4*IMR/5) JN0 = JNP-JS0+1 ENDIF DO J=2,JMR DTDX(j) = DT / ( DL*AE*COSP(J) ) ! PRINT*,'dtdx=',dtdx(j) DTDX5(j) = 0.5*DTDX(j) enddo ! DTDY = DT /(AE*DP) ! PRINT*,'dtdy=',dtdy DTDY5 = 0.5*DTDY ! WRITE(6,*) 'J1=',J1,' J2=', J2 ENDIF ! *********** End Initialization ********************** ! delp = pressure thickness: the psudo-density in a hydrostatic system. DO k=1,NLAY DO j=1,JNP DO i=1,IMR delp1(i,j,k)=DAP(k)+DBK(k)*PS1(i,j) delp2(i,j,k)=DAP(k)+DBK(k)*PS2(i,j) enddo enddo enddo IF(j1/=2) THEN DO IC=1,NC DO L=1,NLAY DO I=1,IMR Q(I, 2,L,IC) = Q(I, 1,L,IC) Q(I,JMR,L,IC) = Q(I,JNP,L,IC) END DO END DO END DO ENDIF ! Compute "tracer density" DO IC=1,NC DO k=1,NLAY DO j=1,JNP DO i=1,IMR DQ(i,j,k,IC) = Q(i,j,k,IC)*delp1(i,j,k) END DO END DO END DO END DO DO k=1,NLAY IF(IGD==0) THEN ! Convert winds on A-Grid to Courant # on C-Grid. CALL A2C(U(1,1,k),V(1,1,k),IMR,JMR,j1,j2,CRX,CRY,dtdx5,DTDY5) else ! Convert winds on C-grid to Courant # DO j=j1,j2 DO i=2,IMR CRX(i,J) = dtdx(j)*U(i-1,j,k) END DO END DO DO j=j1,j2 CRX(1,J) = dtdx(j)*U(IMR,j,k) END DO DO i=1,IMR*JMR CRY(i,2) = DTDY*V(i,1,k) END DO ENDIF ! Determine JS and JN JS = j1 JN = j2 DO j=JS0,j1+1,-1 DO i=1,IMR IF(abs(CRX(i,j))>1.) THEN JS = j go to 2222 ENDIF enddo enddo 2222 continue DO j=JN0,j2-1 DO i=1,IMR IF(abs(CRX(i,j))>1.) THEN JN = j go to 2233 ENDIF enddo enddo 2233 continue IF(j1/=2) then ! Enlarged polar cap. DO i=1,IMR DPI(i, 2,k) = 0. DPI(i,JMR,k) = 0. enddo ENDIF ! ******* Compute horizontal mass fluxes ************ ! N-S component DO j=j1,j2+1 D5 = 0.5 * COSE(j) DO i=1,IMR ymass(i,j) = CRY(i,j)*D5*(delp2(i,j,k) + delp2(i,j-1,k)) enddo enddo DO j=j1,j2 DO i=1,IMR DPI(i,j,k) = (ymass(i,j) - ymass(i,j+1)) * acosp(j) END DO END DO ! Poles sum1 = ymass(IMR,j1 ) sum2 = ymass(IMR,J2+1) DO i=1,IMR-1 sum1 = sum1 + ymass(i,j1 ) sum2 = sum2 + ymass(i,J2+1) enddo sum1 = - sum1 * RCAP sum2 = sum2 * RCAP DO i=1,IMR DPI(i, 1,k) = sum1 DPI(i,JNP,k) = sum2 enddo ! E-W component DO j=j1,j2 DO i=2,IMR PU(i,j) = 0.5 * (delp2(i,j,k) + delp2(i-1,j,k)) enddo enddo DO j=j1,j2 PU(1,j) = 0.5 * (delp2(1,j,k) + delp2(IMR,j,k)) enddo DO j=j1,j2 DO i=1,IMR xmass(i,j) = PU(i,j)*CRX(i,j) END DO END DO DO j=j1,j2 DO i=1,IMR-1 DPI(i,j,k) = DPI(i,j,k) + xmass(i,j) - xmass(i+1,j) END DO END DO DO j=j1,j2 DPI(IMR,j,k) = DPI(IMR,j,k) + xmass(IMR,j) - xmass(1,j) END DO DO j=j1,j2 DO i=1,IMR-1 UA(i,j) = 0.5 * (CRX(i,j)+CRX(i+1,j)) enddo enddo DO j=j1,j2 UA(imr,j) = 0.5 * (CRX(imr,j)+CRX(1,j)) enddo !cccccccccccccccccccccccccccccccccccccccccccccccccccccc ! Rajouts pour LMDZ.3.3 !cccccccccccccccccccccccccccccccccccccccccccccccccccccc DO i=1,IMR DO j=1,JNP VA(i,j)=0. enddo enddo DO i=1,imr*(JMR-1) VA(i,2) = 0.5*(CRY(i,2)+CRY(i,3)) enddo IF(j1==2) THEN IMH = IMR/2 DO i=1,IMH VA(i, 1) = 0.5*(CRY(i,2)-CRY(i+IMH,2)) VA(i+IMH, 1) = -VA(i,1) VA(i, JNP) = 0.5*(CRY(i,JNP)-CRY(i+IMH,JMR)) VA(i+IMH,JNP) = -VA(i,JNP) enddo VA(IMR,1)=VA(1,1) VA(IMR,JNP)=VA(1,JNP) ENDIF ! ****6***0*********0*********0*********0*********0*********0**********72 DO IC=1,NC DO i=1,IMJM wk1(i,1,1) = 0. wk1(i,1,2) = 0. enddo ! E-W advective cross term DO j=J1,J2 IF(J>JS .AND. J=0.) THEN wk1(i,j,1) = qtmp(iiu)+ru*(qtmp(iiu-1)-qtmp(iiu)) else wk1(i,j,1) = qtmp(iiu)+ru*(qtmp(iiu)-qtmp(iiu+1)) ENDIF wk1(i,j,1) = wk1(i,j,1) - qtmp(i) END DO 250 continue END DO IF(JN/=0) THEN DO j=JS+1,JN-1 DO i=1,IMR qtmp(i) = q(i,j,k,IC) enddo qtmp(0) = q(IMR,J,k,IC) qtmp(IMR+1) = q( 1,J,k,IC) DO i=1,imr iu = i - UA(i,j) wk1(i,j,1) = UA(i,j)*(qtmp(iu) - qtmp(iu+1)) enddo enddo ENDIF ! ****6***0*********0*********0*********0*********0*********0**********72 ! Contribution from the N-S advection DO i=1,imr*(j2-j1+1) JT = REAL(J1) - VA(i,j1) wk1(i,j1,2) = VA(i,j1) * (q(i,jt,k,IC) - q(i,jt+1,k,IC)) enddo DO i=1,IMJM wk1(i,1,1) = q(i,1,k,IC) + 0.5*wk1(i,1,1) wk1(i,1,2) = q(i,1,k,IC) + 0.5*wk1(i,1,2) enddo IF(cross) THEN ! Add cross terms in the vertical direction. IF(IORD >= 2) THEN iad = 2 else iad = 1 endif IF(JORD >= 2) THEN jad = 2 else jad = 1 endif CALL xadv(IMR,JNP,j1,j2,wk1(1,1,2),UA,JS,JN,IML,DC2,iad) CALL yadv(IMR,JNP,j1,j2,wk1(1,1,1),VA,PV,W,jad) DO j=1,JNP DO i=1,IMR q(i,j,k,IC) = q(i,j,k,IC) + DC2(i,j) + PV(i,j) enddo enddo ENDIF CALL xtp(IMR,JNP,IML,j1,j2,JN,JS,PU,DQ(1,1,k,IC),wk1(1,1,2) & ,CRX,fx1,xmass,IORD) CALL ytp(IMR,JNP,j1,j2,acosp,RCAP,DQ(1,1,k,IC),wk1(1,1,1),CRY, & DC2,ymass,WK1(1,1,3),wk1(1,1,4),WK1(1,1,5),WK1(1,1,6),JORD) END DO END DO ! ******* Compute vertical mass flux (same unit as PS) *********** ! 1st step: compute total column mass CONVERGENCE. DO j=1,JNP DO i=1,IMR CRY(i,j) = DPI(i,j,1) END DO END DO DO k=2,NLAY DO j=1,JNP DO i=1,IMR CRY(i,j) = CRY(i,j) + DPI(i,j,k) END DO END DO END DO DO j=1,JNP DO i=1,IMR ! 2nd step: compute PS2 (PS at n+1) using the hydrostatic assumption. ! Changes (increases) to surface pressure = total column mass convergence PS2(i,j) = PS1(i,j) + CRY(i,j) ! 3rd step: compute vertical mass flux from mass conservation principle. W(i,j,1) = DPI(i,j,1) - DBK(1)*CRY(i,j) W(i,j,NLAY) = 0. END DO END DO DO k=2,NLAY-1 DO j=1,JNP DO i=1,IMR W(i,j,k) = W(i,j,k-1) + DPI(i,j,k) - DBK(k)*CRY(i,j) END DO END DO END DO DO k=1,NLAY DO j=1,JNP DO i=1,IMR delp2(i,j,k) = DAP(k) + DBK(k)*PS2(i,j) END DO END DO END DO KRD = max(3, KORD) DO IC=1,NC !****6***0*********0*********0*********0*********0*********0**********72 CALL FZPPM(IMR,JNP,NLAY,j1,DQ(1,1,1,IC),W,Q(1,1,1,IC),WK1,DPI, & DC2,CRX,CRY,PU,PV,xmass,ymass,delp1,KRD) ! IF(fill) CALL qckxyz(DQ(1,1,1,IC),DC2,IMR,JNP,NLAY,j1,j2, & cosp,acosp,.FALSE.,IC,NSTEP) ! Recover tracer mixing ratio from "density" using predicted ! "air density" (pressure thickness) at time-level n+1 DO k=1,NLAY DO j=1,JNP DO i=1,IMR Q(i,j,k,IC) = DQ(i,j,k,IC) / delp2(i,j,k) ! PRINT*,'i=',i,'j=',j,'k=',k,'Q(i,j,k,IC)=',Q(i,j,k,IC) enddo enddo enddo IF(j1/=2) THEN DO k=1,NLAY DO I=1,IMR ! j=1 c'est le pôle Sud, j=JNP c'est le pôle Nord Q(I, 2,k,IC) = Q(I, 1,k,IC) Q(I,JMR,k,IC) = Q(I,JNP,k,IC) END DO END DO ENDIF END DO IF(j1/=2) THEN DO k=1,NLAY DO i=1,IMR W(i, 2,k) = W(i, 1,k) W(i,JMR,k) = W(i,JNP,k) END DO END DO ENDIF RETURN END SUBROUTINE ppm3d !****6***0*********0*********0*********0*********0*********0**********72 SUBROUTINE FZPPM(IMR,JNP,NLAY,j1,DQ,WZ,P,DC,DQDT,AR,AL,A6, & flux,wk1,wk2,wz2,delp,KORD) IMPLICIT NONE INTEGER,parameter :: kmax = 150 REAL,parameter :: R23 = 2./3., R3 = 1./3. INTEGER :: IMR,JNP,NLAY,J1,KORD REAL :: WZ(IMR,JNP,NLAY),P(IMR,JNP,NLAY),DC(IMR,JNP,NLAY), & wk1(IMR,*),delp(IMR,JNP,NLAY),DQ(IMR,JNP,NLAY), & DQDT(IMR,JNP,NLAY) ! Assuming JNP >= NLAY REAL :: AR(IMR,*),AL(IMR,*),A6(IMR,*),flux(IMR,*),wk2(IMR,*), & wz2(IMR,*) INTEGER :: JMR,IMJM,NLAYM1,LMT,K,I,J REAL :: c0,c1,c2,tmp,qmax,qmin,a,b,fct,a1,a2,cm,cp JMR = JNP - 1 IMJM = IMR*JNP NLAYM1 = NLAY - 1 LMT = KORD - 3 ! ****6***0*********0*********0*********0*********0*********0**********72 ! Compute DC for PPM ! ****6***0*********0*********0*********0*********0*********0**********72 DO k=1,NLAYM1 DO i=1,IMJM DQDT(i,1,k) = P(i,1,k+1) - P(i,1,k) END DO END DO DO k=2,NLAYM1 DO I=1,IMJM c0 = delp(i,1,k) / (delp(i,1,k-1)+delp(i,1,k)+delp(i,1,k+1)) c1 = (delp(i,1,k-1)+0.5*delp(i,1,k))/(delp(i,1,k+1)+delp(i,1,k)) c2 = (delp(i,1,k+1)+0.5*delp(i,1,k))/(delp(i,1,k-1)+delp(i,1,k)) tmp = c0*(c1*DQDT(i,1,k) + c2*DQDT(i,1,k-1)) Qmax = max(P(i,1,k-1),P(i,1,k),P(i,1,k+1)) - P(i,1,k) Qmin = P(i,1,k) - min(P(i,1,k-1),P(i,1,k),P(i,1,k+1)) DC(i,1,k) = sign(min(abs(tmp),Qmax,Qmin), tmp) END DO END DO ! ****6***0*********0*********0*********0*********0*********0**********72 ! Loop over latitudes (to save memory) ! ****6***0*********0*********0*********0*********0*********0**********72 DO j=1,JNP IF((j==2 .OR. j==JMR) .AND. j1/=2) goto 2000 DO k=1,NLAY DO i=1,IMR wz2(i,k) = WZ(i,j,k) wk1(i,k) = P(i,j,k) wk2(i,k) = delp(i,j,k) flux(i,k) = DC(i,j,k) !this flux is actually the monotone slope enddo enddo !****6***0*********0*********0*********0*********0*********0**********72 ! Compute first guesses at cell interfaces ! First guesses are required to be continuous. ! ****6***0*********0*********0*********0*********0*********0**********72 ! three-cell parabolic subgrid distribution at model top ! two-cell parabolic with zero gradient subgrid distribution ! at the surface. ! First guess top edge value DO i=1,IMR ! three-cell PPM ! Compute a,b, and c of q = aP**2 + bP + c using cell averages and delp a = 3.*( DQDT(i,j,2) - DQDT(i,j,1)*(wk2(i,2)+wk2(i,3))/ & (wk2(i,1)+wk2(i,2)) ) / & ( (wk2(i,2)+wk2(i,3))*(wk2(i,1)+wk2(i,2)+wk2(i,3)) ) b = 2.*DQDT(i,j,1)/(wk2(i,1)+wk2(i,2)) - & R23*a*(2.*wk2(i,1)+wk2(i,2)) AL(i,1) = wk1(i,1) - wk2(i,1)*(R3*a*wk2(i,1) + 0.5*b) AL(i,2) = wk2(i,1)*(a*wk2(i,1) + b) + AL(i,1) ! Check if change sign IF(wk1(i,1)*AL(i,1)<=0.) THEN AL(i,1) = 0. flux(i,1) = 0. else flux(i,1) = wk1(i,1) - AL(i,1) endif END DO ! Bottom DO i=1,IMR ! 2-cell PPM with zero gradient right at the surface fct = DQDT(i,j,NLAYM1)*wk2(i,NLAY)**2 / & ( (wk2(i,NLAY)+wk2(i,NLAYM1))*(2.*wk2(i,NLAY)+wk2(i,NLAYM1))) AR(i,NLAY) = wk1(i,NLAY) + fct AL(i,NLAY) = wk1(i,NLAY) - (fct+fct) IF(wk1(i,NLAY)*AR(i,NLAY)<=0.) AR(i,NLAY) = 0. flux(i,NLAY) = AR(i,NLAY) - wk1(i,NLAY) END DO !****6***0*********0*********0*********0*********0*********0**********72 ! 4th order interpolation in the interior. !****6***0*********0*********0*********0*********0*********0**********72 DO k=3,NLAYM1 DO i=1,IMR c1 = DQDT(i,j,k-1)*wk2(i,k-1) / (wk2(i,k-1)+wk2(i,k)) c2 = 2. / (wk2(i,k-2)+wk2(i,k-1)+wk2(i,k)+wk2(i,k+1)) A1 = (wk2(i,k-2)+wk2(i,k-1)) / (2.*wk2(i,k-1)+wk2(i,k)) A2 = (wk2(i,k )+wk2(i,k+1)) / (2.*wk2(i,k)+wk2(i,k-1)) AL(i,k) = wk1(i,k-1) + c1 + c2 * & ( wk2(i,k )*(c1*(A1 - A2)+A2*flux(i,k-1)) - & wk2(i,k-1)*A1*flux(i,k) ) ! print *,'AL1',i,k, AL(i,k) END DO END DO DO i=1,IMR*NLAYM1 AR(i,1) = AL(i,2) ! print *,'AR1',i,AR(i,1) END DO DO i=1,IMR*NLAY A6(i,1) = 3.*(wk1(i,1)+wk1(i,1) - (AL(i,1)+AR(i,1))) ! print *,'A61',i,A6(i,1) END DO !****6***0*********0*********0*********0*********0*********0**********72 ! Top & Bot always monotonic CALL lmtppm(flux(1,1),A6(1,1),AR(1,1),AL(1,1),wk1(1,1),IMR,0) CALL lmtppm(flux(1,NLAY),A6(1,NLAY),AR(1,NLAY),AL(1,NLAY), & wk1(1,NLAY),IMR,0) ! Interior depending on KORD IF(LMT<=2) & CALL lmtppm(flux(1,2),A6(1,2),AR(1,2),AL(1,2),wk1(1,2), & IMR*(NLAY-2),LMT) !****6***0*********0*********0*********0*********0*********0**********72 DO i=1,IMR*NLAYM1 IF(wz2(i,1)>0.) THEN CM = wz2(i,1) / wk2(i,1) flux(i,2) = AR(i,1)+0.5*CM*(AL(i,1)-AR(i,1)+A6(i,1)*(1.-R23*CM)) else ! print *,'test2-0',i,j,wz2(i,1),wk2(i,2) CP= wz2(i,1) / wk2(i,2) ! print *,'testCP',CP flux(i,2) = AL(i,2)+0.5*CP*(AL(i,2)-AR(i,2)-A6(i,2)*(1.+R23*CP)) ! print *,'test2',i, AL(i,2),AR(i,2),A6(i,2),R23 ENDIF END DO DO i=1,IMR*NLAYM1 flux(i,2) = wz2(i,1) * flux(i,2) END DO DO i=1,IMR DQ(i,j, 1) = DQ(i,j, 1) - flux(i, 2) DQ(i,j,NLAY) = DQ(i,j,NLAY) + flux(i,NLAY) END DO DO k=2,NLAYM1 DO i=1,IMR DQ(i,j,k) = DQ(i,j,k) + flux(i,k) - flux(i,k+1) END DO END DO 2000 continue END DO RETURN END SUBROUTINE fzppm SUBROUTINE xtp(IMR,JNP,IML,j1,j2,JN,JS,PU,DQ,Q,UC, & fx1,xmass,IORD) IMPLICIT NONE INTEGER :: IMR,JNP,IML,j1,j2,JN,JS,IORD REAL :: PU,DQ,Q,UC,fx1,xmass REAL :: dc,qtmp INTEGER :: ISAVE(IMR) dimension UC(IMR,*),DC(-IML:IMR+IML+1),xmass(IMR,JNP) & ,fx1(IMR+1),DQ(IMR,JNP),qtmp(-IML:IMR+1+IML) dimension PU(IMR,JNP),Q(IMR,JNP) INTEGER :: jvan,j1vl,j2vl,j,i,iu,itmp,ist,imp REAL :: rut IMP = IMR + 1 ! van Leer at high latitudes jvan = max(1,JNP/18) j1vl = j1+jvan j2vl = j2-jvan DO j=j1,j2 DO i=1,IMR qtmp(i) = q(i,j) enddo IF(j>=JN .OR. j<=JS) goto 2222 ! ************* Eulerian ********** qtmp(0) = q(IMR,J) qtmp(-1) = q(IMR-1,J) qtmp(IMP) = q(1,J) qtmp(IMP+1) = q(2,J) IF(IORD==1 .OR. j==j1 .OR. j==j2) THEN DO i=1,IMR iu = REAL(i) - uc(i,j) fx1(i) = qtmp(iu) END DO ELSE CALL xmist(IMR,IML,Qtmp,DC) DC(0) = DC(IMR) IF(IORD==2 .OR. j<=j1vl .OR. j>=j2vl) THEN DO i=1,IMR iu = REAL(i) - uc(i,j) fx1(i) = qtmp(iu) + DC(iu)*(sign(1.,uc(i,j))-uc(i,j)) END DO else CALL fxppm(IMR,IML,UC(1,j),Qtmp,DC,fx1,IORD) ENDIF ENDIF DO i=1,IMR fx1(i) = fx1(i)*xmass(i,j) END DO goto 1309 ! ***** Conservative (flux-form) Semi-Lagrangian transport ***** 2222 continue DO i=-IML,0 qtmp(i) = q(IMR+i,j) qtmp(IMP-i) = q(1-i,j) enddo IF(IORD==1 .OR. j==j1 .OR. j==j2) THEN DO i=1,IMR itmp = INT(uc(i,j)) ISAVE(i) = i - itmp iu = i - uc(i,j) fx1(i) = (uc(i,j) - itmp)*qtmp(iu) END DO ELSE CALL xmist(IMR,IML,Qtmp,DC) DO i=-IML,0 DC(i) = DC(IMR+i) DC(IMP-i) = DC(1-i) enddo DO i=1,IMR itmp = INT(uc(i,j)) rut = uc(i,j) - itmp ISAVE(i) = i - itmp iu = i - uc(i,j) fx1(i) = rut*(qtmp(iu) + DC(iu)*(sign(1.,rut) - rut)) END DO ENDIF DO i=1,IMR IF(uc(i,j)>1.) THEN !DIR$ NOVECTOR DO ist = ISAVE(i),i-1 fx1(i) = fx1(i) + qtmp(ist) enddo elseIF(uc(i,j)<-1.) THEN DO ist = i,ISAVE(i)-1 fx1(i) = fx1(i) - qtmp(ist) enddo !DIR$ VECTOR ENDIF END DO DO i=1,IMR fx1(i) = PU(i,j)*fx1(i) enddo ! *************************************** 1309 fx1(IMP) = fx1(1) DO i=1,IMR DQ(i,j) = DQ(i,j) + fx1(i)-fx1(i+1) END DO ! *************************************** END DO RETURN END SUBROUTINE xtp SUBROUTINE fxppm(IMR,IML,UT,P,DC,flux,IORD) IMPLICIT NONE INTEGER :: IMR,IML,IORD REAL :: UT,P,DC,flux REAL,parameter :: R3 = 1./3., R23 = 2./3. DIMENSION UT(*),flux(*),P(-IML:IMR+IML+1),DC(-IML:IMR+IML+1) REAL :: AR(0:IMR),AL(0:IMR),A6(0:IMR) INTEGER :: LMT,IMP,JLVL,i ! LOGICAL first ! data first /.TRUE./ ! SAVE LMT ! IF(first) THEN ! correction calcul de LMT a chaque passage pour pouvoir choisir ! plusieurs schemas PPM pour differents traceurs ! IF (IORD.LE.0) THEN ! IF(IMR.GE.144) THEN ! LMT = 0 ! elseif(IMR.GE.72) THEN ! LMT = 1 ! else ! LMT = 2 ! endif ! else ! LMT = IORD - 3 ! ENDIF LMT = IORD - 3 ! WRITE(6,*) 'PPM option in E-W direction = ', LMT ! first = .FALSE. ! END IF DO i=1,IMR AL(i) = 0.5*(p(i-1)+p(i)) + (DC(i-1) - DC(i))*R3 END DO DO i=1,IMR-1 AR(i) = AL(i+1) END DO AR(IMR) = AL(1) DO i=1,IMR A6(i) = 3.*(p(i)+p(i) - (AL(i)+AR(i))) END DO IF(LMT<=2) CALL lmtppm(DC(1),A6(1),AR(1),AL(1),P(1),IMR,LMT) AL(0) = AL(IMR) AR(0) = AR(IMR) A6(0) = A6(IMR) DO i=1,IMR IF(UT(i)>0.) THEN flux(i) = AR(i-1) + 0.5*UT(i)*(AL(i-1) - AR(i-1) + & A6(i-1)*(1.-R23*UT(i)) ) else flux(i) = AL(i) - 0.5*UT(i)*(AR(i) - AL(i) + & A6(i)*(1.+R23*UT(i))) ENDIF enddo RETURN END SUBROUTINE fxppm SUBROUTINE xmist(IMR,IML,P,DC) IMPLICIT NONE INTEGER :: IMR,IML REAL,parameter :: R24 = 1./24. REAL :: P(-IML:IMR+1+IML),DC(-IML:IMR+1+IML) INTEGER :: i REAL :: tmp,pmax,pmin DO i=1,IMR tmp = R24*(8.*(p(i+1) - p(i-1)) + p(i-2) - p(i+2)) Pmax = max(P(i-1), p(i), p(i+1)) - p(i) Pmin = p(i) - min(P(i-1), p(i), p(i+1)) DC(i) = sign(min(abs(tmp),Pmax,Pmin), tmp) END DO RETURN END SUBROUTINE xmist SUBROUTINE ytp(IMR,JNP,j1,j2,acosp,RCAP,DQ,P,VC,DC2 & ,ymass,fx,A6,AR,AL,JORD) IMPLICIT NONE INTEGER :: IMR,JNP,j1,j2,JORD REAL :: acosp,RCAP,DQ,P,VC,DC2,ymass,fx,A6,AR,AL dimension P(IMR,JNP),VC(IMR,JNP),ymass(IMR,JNP) & ,DC2(IMR,JNP),DQ(IMR,JNP),acosp(JNP) ! Work array DIMENSION fx(IMR,JNP),AR(IMR,JNP),AL(IMR,JNP),A6(IMR,JNP) INTEGER :: JMR,len,i,jt,j REAL :: sum1,sum2 JMR = JNP - 1 len = IMR*(J2-J1+2) IF(JORD==1) THEN DO i=1,len JT = REAL(J1) - VC(i,J1) fx(i,j1) = p(i,JT) END DO else CALL ymist(IMR,JNP,j1,P,DC2,4) IF(JORD<=0 .OR. JORD>=3) THEN CALL fyppm(VC,P,DC2,fx,IMR,JNP,j1,j2,A6,AR,AL,JORD) else DO i=1,len JT = REAL(J1) - VC(i,J1) fx(i,j1) = p(i,JT) + (sign(1.,VC(i,j1))-VC(i,j1))*DC2(i,JT) END DO ENDIF ENDIF DO i=1,len fx(i,j1) = fx(i,j1)*ymass(i,j1) END DO DO j=j1,j2 DO i=1,IMR DQ(i,j) = DQ(i,j) + (fx(i,j) - fx(i,j+1)) * acosp(j) END DO END DO ! Poles sum1 = fx(IMR,j1 ) sum2 = fx(IMR,J2+1) DO i=1,IMR-1 sum1 = sum1 + fx(i,j1 ) sum2 = sum2 + fx(i,J2+1) enddo sum1 = DQ(1, 1) - sum1 * RCAP sum2 = DQ(1,JNP) + sum2 * RCAP DO i=1,IMR DQ(i, 1) = sum1 DQ(i,JNP) = sum2 enddo IF(j1/=2) THEN DO i=1,IMR DQ(i, 2) = sum1 DQ(i,JMR) = sum2 enddo ENDIF RETURN END SUBROUTINE ytp subroutine ymist(IMR,JNP,j1,P,DC,ID) IMPLICIT NONE INTEGER :: IMR,JNP,j1,ID REAL,parameter :: R24 = 1./24. REAL :: P(IMR,JNP),DC(IMR,JNP) INTEGER :: iimh,jmr,ijm3,imh,i REAL :: pmax,pmin,tmp IMH = IMR / 2 JMR = JNP - 1 IJM3 = IMR*(JMR-3) IF(ID==2) THEN DO i=1,IMR*(JMR-1) tmp = 0.25*(p(i,3) - p(i,1)) Pmax = max(p(i,1),p(i,2),p(i,3)) - p(i,2) Pmin = p(i,2) - min(p(i,1),p(i,2),p(i,3)) DC(i,2) = sign(min(abs(tmp),Pmin,Pmax),tmp) END DO ELSE DO i=1,IMH ! J=2 tmp = (8.*(p(i,3) - p(i,1)) + p(i+IMH,2) - p(i,4))*R24 Pmax = max(p(i,1),p(i,2),p(i,3)) - p(i,2) Pmin = p(i,2) - min(p(i,1),p(i,2),p(i,3)) DC(i,2) = sign(min(abs(tmp),Pmin,Pmax),tmp) ! J=JMR tmp=(8.*(p(i,JNP)-p(i,JMR-1))+p(i,JMR-2)-p(i+IMH,JMR))*R24 Pmax = max(p(i,JMR-1),p(i,JMR),p(i,JNP)) - p(i,JMR) Pmin = p(i,JMR) - min(p(i,JMR-1),p(i,JMR),p(i,JNP)) DC(i,JMR) = sign(min(abs(tmp),Pmin,Pmax),tmp) END DO DO i=IMH+1,IMR ! J=2 tmp = (8.*(p(i,3) - p(i,1)) + p(i-IMH,2) - p(i,4))*R24 Pmax = max(p(i,1),p(i,2),p(i,3)) - p(i,2) Pmin = p(i,2) - min(p(i,1),p(i,2),p(i,3)) DC(i,2) = sign(min(abs(tmp),Pmin,Pmax),tmp) ! J=JMR tmp=(8.*(p(i,JNP)-p(i,JMR-1))+p(i,JMR-2)-p(i-IMH,JMR))*R24 Pmax = max(p(i,JMR-1),p(i,JMR),p(i,JNP)) - p(i,JMR) Pmin = p(i,JMR) - min(p(i,JMR-1),p(i,JMR),p(i,JNP)) DC(i,JMR) = sign(min(abs(tmp),Pmin,Pmax),tmp) END DO DO i=1,IJM3 tmp = (8.*(p(i,4) - p(i,2)) + p(i,1) - p(i,5))*R24 Pmax = max(p(i,2),p(i,3),p(i,4)) - p(i,3) Pmin = p(i,3) - min(p(i,2),p(i,3),p(i,4)) DC(i,3) = sign(min(abs(tmp),Pmin,Pmax),tmp) END DO ENDIF IF(j1/=2) THEN DO i=1,IMR DC(i,1) = 0. DC(i,JNP) = 0. enddo else ! Determine slopes in polar caps for scalars! DO i=1,IMH ! South tmp = 0.25*(p(i,2) - p(i+imh,2)) Pmax = max(p(i,2),p(i,1), p(i+imh,2)) - p(i,1) Pmin = p(i,1) - min(p(i,2),p(i,1), p(i+imh,2)) DC(i,1)=sign(min(abs(tmp),Pmax,Pmin),tmp) ! North. tmp = 0.25*(p(i+imh,JMR) - p(i,JMR)) Pmax = max(p(i+imh,JMR),p(i,jnp), p(i,JMR)) - p(i,JNP) Pmin = p(i,JNP) - min(p(i+imh,JMR),p(i,jnp), p(i,JMR)) DC(i,JNP) = sign(min(abs(tmp),Pmax,pmin),tmp) END DO DO i=imh+1,IMR DC(i, 1) = - DC(i-imh, 1) DC(i,JNP) = - DC(i-imh,JNP) END DO ENDIF RETURN END SUBROUTINE ymist SUBROUTINE fyppm(VC,P,DC,flux,IMR,JNP,j1,j2,A6,AR,AL,JORD) IMPLICIT NONE INTEGER :: IMR,JNP,j1,j2,JORD REAL,parameter :: R3 = 1./3., R23 = 2./3. REAL :: VC(IMR,*),flux(IMR,*),P(IMR,*),DC(IMR,*) ! Local work arrays. REAL :: AR(IMR,JNP),AL(IMR,JNP),A6(IMR,JNP) INTEGER :: LMT,i INTEGER :: IMH,JMR,j11,IMJM1,len ! LOGICAL first ! data first /.TRUE./ ! SAVE LMT IMH = IMR / 2 JMR = JNP - 1 j11 = j1-1 IMJM1 = IMR*(J2-J1+2) len = IMR*(J2-J1+3) ! IF(first) THEN ! IF(JORD.LE.0) THEN ! IF(JMR.GE.90) THEN ! LMT = 0 ! elseif(JMR.GE.45) THEN ! LMT = 1 ! else ! LMT = 2 ! endif ! else ! LMT = JORD - 3 ! END IF ! first = .FALSE. ! ENDIF ! modifs pour pouvoir choisir plusieurs schemas PPM LMT = JORD - 3 DO i=1,IMR*JMR AL(i,2) = 0.5*(p(i,1)+p(i,2)) + (DC(i,1) - DC(i,2))*R3 AR(i,1) = AL(i,2) END DO !Poles: DO i=1,IMH AL(i,1) = AL(i+IMH,2) AL(i+IMH,1) = AL(i,2) AR(i,JNP) = AR(i+IMH,JMR) AR(i+IMH,JNP) = AR(i,JMR) ENDDO !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! Rajout pour LMDZ.3.3 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc AR(IMR,1)=AL(1,1) AR(IMR,JNP)=AL(1,JNP) !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc DO i=1,len A6(i,j11) = 3.*(p(i,j11)+p(i,j11) - (AL(i,j11)+AR(i,j11))) END DO IF(LMT<=2) CALL lmtppm(DC(1,j11),A6(1,j11),AR(1,j11) & ,AL(1,j11),P(1,j11),len,LMT) ! DO i=1,IMJM1 IF(VC(i,j1)>0.) THEN flux(i,j1) = AR(i,j11) + 0.5*VC(i,j1)*(AL(i,j11) - AR(i,j11) + & A6(i,j11)*(1.-R23*VC(i,j1)) ) else flux(i,j1) = AL(i,j1) - 0.5*VC(i,j1)*(AR(i,j1) - AL(i,j1) + & A6(i,j1)*(1.+R23*VC(i,j1))) ENDIF END DO RETURN END SUBROUTINE fyppm SUBROUTINE yadv(IMR,JNP,j1,j2,p,VA,ady,wk,IAD) IMPLICIT NONE INTEGER :: IMR,JNP,j1,j2,IAD REAL :: p(IMR,JNP),ady(IMR,JNP),VA(IMR,JNP) REAL :: WK(IMR,-1:JNP+2) INTEGER :: JMR,IMH,i,j,jp REAL :: rv,a1,b1,sum1,sum2 JMR = JNP-1 IMH = IMR/2 DO j=1,JNP DO i=1,IMR wk(i,j) = p(i,j) enddo enddo ! Poles: DO i=1,IMH wk(i, -1) = p(i+IMH,3) wk(i+IMH,-1) = p(i,3) wk(i, 0) = p(i+IMH,2) wk(i+IMH,0) = p(i,2) wk(i,JNP+1) = p(i+IMH,JMR) wk(i+IMH,JNP+1) = p(i,JMR) wk(i,JNP+2) = p(i+IMH,JNP-2) wk(i+IMH,JNP+2) = p(i,JNP-2) enddo ! WRITE(*,*) 'toto 1' ! -------------------------------- IF(IAD==2) THEN DO j=j1-1,j2+1 DO i=1,IMR ! WRITE(*,*) 'avt NINT','i=',i,'j=',j JP = NINT(VA(i,j)) rv = JP - VA(i,j) ! WRITE(*,*) 'VA=',VA(i,j), 'JP1=',JP,'rv=',rv JP = j - JP ! WRITE(*,*) 'JP2=',JP a1 = 0.5*(wk(i,jp+1)+wk(i,jp-1)) - wk(i,jp) b1 = 0.5*(wk(i,jp+1)-wk(i,jp-1)) ! WRITE(*,*) 'a1=',a1,'b1=',b1 ady(i,j) = wk(i,jp) + rv*(a1*rv + b1) - wk(i,j) enddo enddo ! WRITE(*,*) 'toto 2' ELSEIF(IAD==1) THEN DO j=j1-1,j2+1 DO i=1,imr JP = REAL(j)-VA(i,j) ady(i,j) = VA(i,j)*(wk(i,jp)-wk(i,jp+1)) enddo enddo ENDIF IF(j1/=2) THEN sum1 = 0. sum2 = 0. DO i=1,imr sum1 = sum1 + ady(i,2) sum2 = sum2 + ady(i,JMR) enddo sum1 = sum1 / IMR sum2 = sum2 / IMR DO i=1,imr ady(i, 2) = sum1 ady(i,JMR) = sum2 ady(i, 1) = sum1 ady(i,JNP) = sum2 enddo else ! Poles: sum1 = 0. sum2 = 0. DO i=1,imr sum1 = sum1 + ady(i,1) sum2 = sum2 + ady(i,JNP) enddo sum1 = sum1 / IMR sum2 = sum2 / IMR DO i=1,imr ady(i, 1) = sum1 ady(i,JNP) = sum2 enddo endif RETURN END SUBROUTINE yadv SUBROUTINE xadv(IMR,JNP,j1,j2,p,UA,JS,JN,IML,adx,IAD) IMPLICIT NONE INTEGER :: IMR,JNP,j1,j2,JS,JN,IML,IAD REAL :: p(IMR,JNP),adx(IMR,JNP),qtmp(-IMR:IMR+IMR),UA(IMR,JNP) INTEGER :: JMR,j,i,ip,iu,iiu REAL :: ru,a1,b1 JMR = JNP-1 DO j=j1,j2 IF(J>JS .AND. J=0.) THEN adx(i,j) = qtmp(iiu)+ru*(qtmp(iiu-1)-qtmp(iiu)) else adx(i,j) = qtmp(iiu)+ru*(qtmp(iiu)-qtmp(iiu+1)) ENDIF enddo ENDIF DO i=1,IMR adx(i,j) = adx(i,j) - p(i,j) enddo 1309 continue END DO ! Eulerian upwind DO j=JS+1,JN-1 DO i=1,IMR qtmp(i) = p(i,j) enddo qtmp(0) = p(IMR,J) qtmp(IMR+1) = p(1,J) IF(IAD==2) THEN qtmp(-1) = p(IMR-1,J) qtmp(IMR+2) = p(2,J) DO i=1,imr IP = NINT(UA(i,j)) ru = IP - UA(i,j) IP = i - IP a1 = 0.5*(qtmp(ip+1)+qtmp(ip-1)) - qtmp(ip) b1 = 0.5*(qtmp(ip+1)-qtmp(ip-1)) adx(i,j) = qtmp(ip)- p(i,j) + ru*(a1*ru + b1) enddo ELSEIF(IAD==1) THEN ! 1st order DO i=1,IMR IP = i - UA(i,j) adx(i,j) = UA(i,j)*(qtmp(ip)-qtmp(ip+1)) enddo ENDIF enddo IF(j1/=2) THEN DO i=1,IMR adx(i, 2) = 0. adx(i,JMR) = 0. enddo endif ! set cross term due to x-adv at the poles to zero. DO i=1,IMR adx(i, 1) = 0. adx(i,JNP) = 0. enddo RETURN END SUBROUTINE xadv SUBROUTINE lmtppm(DC,A6,AR,AL,P,IM,LMT) IMPLICIT NONE ! A6 = CURVATURE OF THE TEST PARABOLA ! AR = RIGHT EDGE VALUE OF THE TEST PARABOLA ! AL = LEFT EDGE VALUE OF THE TEST PARABOLA ! DC = 0.5 * MISMATCH ! P = CELL-AVERAGED VALUE ! IM = VECTOR LENGTH ! OPTIONS: ! LMT = 0: FULL MONOTONICITY ! LMT = 1: SEMI-MONOTONIC CONSTRAINT (NO UNDERSHOOTS) ! LMT = 2: POSITIVE-DEFINITE CONSTRAINT REAL,parameter :: R12 = 1./12. REAL :: A6(IM),AR(IM),AL(IM),P(IM),DC(IM) INTEGER :: IM,LMT INTEGER :: i REAL :: da1,da2,a6da,fmin IF(LMT==0) THEN ! Full constraint DO i=1,IM IF(DC(i)==0.) THEN AR(i) = p(i) AL(i) = p(i) A6(i) = 0. else da1 = AR(i) - AL(i) da2 = da1**2 A6DA = A6(i)*da1 IF(A6DA < -da2) THEN A6(i) = 3.*(AL(i)-p(i)) AR(i) = AL(i) - A6(i) elseif(A6DA > da2) THEN A6(i) = 3.*(AR(i)-p(i)) AL(i) = AR(i) - A6(i) ENDIF ENDIF END DO elseif(LMT==1) THEN ! Semi-monotonic constraint DO i=1,IM IF(abs(AR(i)-AL(i)) >= -A6(i)) go to 150 IF(p(i) AL(i)) THEN A6(i) = 3.*(AL(i)-p(i)) AR(i) = AL(i) - A6(i) else A6(i) = 3.*(AR(i)-p(i)) AL(i) = AR(i) - A6(i) ENDIF 150 continue END DO elseif(LMT==2) THEN DO i=1,IM IF(abs(AR(i)-AL(i)) >= -A6(i)) go to 250 fmin = p(i) + 0.25*(AR(i)-AL(i))**2/A6(i) + A6(i)*R12 IF(fmin>=0.) go to 250 IF(p(i) AL(i)) THEN A6(i) = 3.*(AL(i)-p(i)) AR(i) = AL(i) - A6(i) else A6(i) = 3.*(AR(i)-p(i)) AL(i) = AR(i) - A6(i) ENDIF 250 continue END DO ENDIF RETURN END SUBROUTINE lmtppm SUBROUTINE A2C(U,V,IMR,JMR,j1,j2,CRX,CRY,dtdx5,DTDY5) IMPLICIT NONE INTEGER :: IMR,JMR,j1,j2 REAL :: U(IMR,*),V(IMR,*),CRX(IMR,*),CRY(IMR,*),DTDX5(*),DTDY5 INTEGER :: i,j DO j=j1,j2 DO i=2,IMR CRX(i,J) = dtdx5(j)*(U(i,j)+U(i-1,j)) END DO END DO DO j=j1,j2 CRX(1,J) = dtdx5(j)*(U(1,j)+U(IMR,j)) END DO DO i=1,IMR*JMR CRY(i,2) = DTDY5*(V(i,2)+V(i,1)) END DO RETURN END SUBROUTINE a2c SUBROUTINE cosa(cosp,cose,JNP,PI,DP) IMPLICIT NONE INTEGER :: JNP REAL :: cosp(*),cose(*),PI,DP INTEGER :: JMR,j,jeq REAL :: ph5 JMR = JNP-1 DO j=2,JNP ph5 = -0.5*PI + (REAL(J-1)-0.5)*DP cose(j) = cos(ph5) END DO JEQ = (JNP+1) / 2 IF(JMR == 2*(JMR/2) ) THEN DO j=JNP, JEQ+1, -1 cose(j) = cose(JNP+2-j) enddo else ! cell edge at equator. cose(JEQ+1) = 1. DO j=JNP, JEQ+2, -1 cose(j) = cose(JNP+2-j) enddo ENDIF DO j=2,JMR cosp(j) = 0.5*(cose(j)+cose(j+1)) END DO cosp(1) = 0. cosp(JNP) = 0. RETURN END SUBROUTINE cosa SUBROUTINE cosc(cosp,cose,JNP,PI,DP) IMPLICIT NONE INTEGER :: JNP REAL :: cosp(*),cose(*),PI,DP REAL :: phi INTEGER :: j phi = -0.5*PI DO j=2,JNP-1 phi = phi + DP cosp(j) = cos(phi) END DO cosp( 1) = 0. cosp(JNP) = 0. DO j=2,JNP cose(j) = 0.5*(cosp(j)+cosp(j-1)) END DO DO j=2,JNP-1 cosp(j) = 0.5*(cose(j)+cose(j+1)) END DO RETURN END SUBROUTINE cosc SUBROUTINE qckxyz(Q,qtmp,IMR,JNP,NLAY,j1,j2,cosp,acosp, & cross,IC,NSTEP) REAL,parameter :: tiny = 1.E-60 INTEGER :: IMR,JNP,NLAY,j1,j2,IC,NSTEP REAL :: Q(IMR,JNP,NLAY),qtmp(IMR,JNP),cosp(*),acosp(*) LOGICAL :: cross INTEGER :: NLAYM1,len,ip,L,icr,ipy,ipx,i REAL :: qup,qly,dup,sum NLAYM1 = NLAY-1 len = IMR*(j2-j1+1) ip = 0 ! Top layer L = 1 icr = 1 CALL filns(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,ipy,tiny) IF(ipy==0) goto 50 CALL filew(q(1,1,L),qtmp,IMR,JNP,j1,j2,ipx,tiny) IF(ipx==0) goto 50 IF(cross) THEN CALL filcr(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,icr,tiny) ENDIF IF(icr==0) goto 50 ! Vertical filling... DO i=1,len IF( Q(i,j1,1)<0.) THEN ip = ip + 1 Q(i,j1,2) = Q(i,j1,2) + Q(i,j1,1) Q(i,j1,1) = 0. ENDIF enddo 50 continue DO L = 2,NLAYM1 icr = 1 CALL filns(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,ipy,tiny) IF(ipy==0) goto 225 CALL filew(q(1,1,L),qtmp,IMR,JNP,j1,j2,ipx,tiny) IF(ipx==0) go to 225 IF(cross) THEN CALL filcr(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,icr,tiny) ENDIF IF(icr==0) goto 225 DO i=1,len IF( Q(I,j1,L)<0.) THEN ip = ip + 1 ! From above qup = Q(I,j1,L-1) qly = -Q(I,j1,L) dup = min(qly,qup) Q(I,j1,L-1) = qup - dup Q(I,j1,L ) = dup-qly ! Below Q(I,j1,L+1) = Q(I,j1,L+1) + Q(I,j1,L) Q(I,j1,L) = 0. ENDIF ENDDO 225 CONTINUE END DO ! BOTTOM LAYER sum = 0. L = NLAY CALL filns(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,ipy,tiny) IF(ipy==0) goto 911 CALL filew(q(1,1,L),qtmp,IMR,JNP,j1,j2,ipx,tiny) IF(ipx==0) goto 911 CALL filcr(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,icr,tiny) IF(icr==0) goto 911 DO I=1,len IF( Q(I,j1,L)<0.) THEN ip = ip + 1 ! From above qup = Q(I,j1,NLAYM1) qly = -Q(I,j1,L) dup = min(qly,qup) Q(I,j1,NLAYM1) = qup - dup ! From "below" the surface. sum = sum + qly-dup Q(I,j1,L) = 0. ENDIF ENDDO 911 continue IF(ip>IMR) THEN WRITE(6,*) 'IC=',IC,' STEP=',NSTEP, & ' Vertical filling pts=',ip ENDIF IF(sum>1.e-25) THEN WRITE(6,*) IC,NSTEP,' Mass source from the ground=',sum ENDIF RETURN END SUBROUTINE qckxyz SUBROUTINE filcr(q,IMR,JNP,j1,j2,cosp,acosp,icr,tiny) IMPLICIT NONE INTEGER :: IMR,JNP,j1,j2,icr REAL :: q(IMR,*),cosp(*),acosp(*),tiny INTEGER :: i,j REAL :: dq,dn,d0,d1,ds,d2 icr = 0 DO j=j1+1,j2-1 DO i=1,IMR-1 IF(q(i,j)<0.) THEN icr = 1 dq = - q(i,j)*cosp(j) ! N-E dn = q(i+1,j+1)*cosp(j+1) d0 = max(0.,dn) d1 = min(dq,d0) q(i+1,j+1) = (dn - d1)*acosp(j+1) dq = dq - d1 ! S-E ds = q(i+1,j-1)*cosp(j-1) d0 = max(0.,ds) d2 = min(dq,d0) q(i+1,j-1) = (ds - d2)*acosp(j-1) q(i,j) = (d2 - dq)*acosp(j) + tiny ENDIF END DO IF(icr==0 .AND. q(IMR,j)>=0.) goto 65 DO i=2,IMR IF(q(i,j)<0.) THEN icr = 1 dq = - q(i,j)*cosp(j) ! N-W dn = q(i-1,j+1)*cosp(j+1) d0 = max(0.,dn) d1 = min(dq,d0) q(i-1,j+1) = (dn - d1)*acosp(j+1) dq = dq - d1 ! S-W ds = q(i-1,j-1)*cosp(j-1) d0 = max(0.,ds) d2 = min(dq,d0) q(i-1,j-1) = (ds - d2)*acosp(j-1) q(i,j) = (d2 - dq)*acosp(j) + tiny ENDIF END DO ! ***************************************** ! i=1 i=1 IF(q(i,j)<0.) THEN icr = 1 dq = - q(i,j)*cosp(j) ! N-W dn = q(IMR,j+1)*cosp(j+1) d0 = max(0.,dn) d1 = min(dq,d0) q(IMR,j+1) = (dn - d1)*acosp(j+1) dq = dq - d1 ! S-W ds = q(IMR,j-1)*cosp(j-1) d0 = max(0.,ds) d2 = min(dq,d0) q(IMR,j-1) = (ds - d2)*acosp(j-1) q(i,j) = (d2 - dq)*acosp(j) + tiny ENDIF ! ***************************************** ! i=IMR i=IMR IF(q(i,j)<0.) THEN icr = 1 dq = - q(i,j)*cosp(j) ! N-E dn = q(1,j+1)*cosp(j+1) d0 = max(0.,dn) d1 = min(dq,d0) q(1,j+1) = (dn - d1)*acosp(j+1) dq = dq - d1 ! S-E ds = q(1,j-1)*cosp(j-1) d0 = max(0.,ds) d2 = min(dq,d0) q(1,j-1) = (ds - d2)*acosp(j-1) q(i,j) = (d2 - dq)*acosp(j) + tiny ENDIF ! ***************************************** 65 continue END DO DO i=1,IMR IF(q(i,j1)<0. .OR. q(i,j2)<0.) THEN icr = 1 goto 80 ENDIF enddo 80 continue IF(q(1,1)<0. .OR. q(1,jnp)<0.) THEN icr = 1 ENDIF RETURN END SUBROUTINE filcr SUBROUTINE filns(q,IMR,JNP,j1,j2,cosp,acosp,ipy,tiny) IMPLICIT NONE INTEGER :: IMR,JNP,j1,j2,ipy REAL :: q(IMR,*),cosp(*),acosp(*),tiny REAL :: DP,CAP1,dq,dn,d0,d1,ds,d2 INTEGER :: i,j ! LOGICAL first ! data first /.TRUE./ ! save cap1 ! IF(first) THEN DP = 4.*ATAN(1.)/REAL(JNP-1) CAP1 = IMR*(1.-COS((j1-1.5)*DP))/DP ! first = .FALSE. ! END IF ipy = 0 DO j=j1+1,j2-1 DO i=1,IMR IF(q(i,j)<0.) THEN ipy = 1 dq = - q(i,j)*cosp(j) ! North dn = q(i,j+1)*cosp(j+1) d0 = max(0.,dn) d1 = min(dq,d0) q(i,j+1) = (dn - d1)*acosp(j+1) dq = dq - d1 ! South ds = q(i,j-1)*cosp(j-1) d0 = max(0.,ds) d2 = min(dq,d0) q(i,j-1) = (ds - d2)*acosp(j-1) q(i,j) = (d2 - dq)*acosp(j) + tiny ENDIF END DO END DO DO i=1,imr IF(q(i,j1)<0.) THEN ipy = 1 dq = - q(i,j1)*cosp(j1) ! North dn = q(i,j1+1)*cosp(j1+1) d0 = max(0.,dn) d1 = min(dq,d0) q(i,j1+1) = (dn - d1)*acosp(j1+1) q(i,j1) = (d1 - dq)*acosp(j1) + tiny ENDIF enddo j = j2 DO i=1,imr IF(q(i,j)<0.) THEN ipy = 1 dq = - q(i,j)*cosp(j) ! South ds = q(i,j-1)*cosp(j-1) d0 = max(0.,ds) d2 = min(dq,d0) q(i,j-1) = (ds - d2)*acosp(j-1) q(i,j) = (d2 - dq)*acosp(j) + tiny ENDIF enddo ! Check Poles. IF(q(1,1)<0.) THEN dq = q(1,1)*cap1/REAL(IMR)*acosp(j1) DO i=1,imr q(i,1) = 0. q(i,j1) = q(i,j1) + dq IF(q(i,j1)<0.) ipy = 1 enddo ENDIF IF(q(1,JNP)<0.) THEN dq = q(1,JNP)*cap1/REAL(IMR)*acosp(j2) DO i=1,imr q(i,JNP) = 0. q(i,j2) = q(i,j2) + dq IF(q(i,j2)<0.) ipy = 1 enddo ENDIF RETURN END SUBROUTINE filns SUBROUTINE filew(q,qtmp,IMR,JNP,j1,j2,ipx,tiny) IMPLICIT NONE INTEGER :: IMR,JNP,j1,j2,ipx REAL :: q(IMR,*),qtmp(JNP,IMR),tiny INTEGER :: i,j REAL :: d0,d1,d2 ipx = 0 ! Copy & swap direction for vectorization. DO i=1,imr DO j=j1,j2 qtmp(j,i) = q(i,j) END DO END DO DO i=2,imr-1 DO j=j1,j2 IF(qtmp(j,i)<0.) THEN ipx = 1 ! west d0 = max(0.,qtmp(j,i-1)) d1 = min(-qtmp(j,i),d0) qtmp(j,i-1) = qtmp(j,i-1) - d1 qtmp(j,i) = qtmp(j,i) + d1 ! east d0 = max(0.,qtmp(j,i+1)) d2 = min(-qtmp(j,i),d0) qtmp(j,i+1) = qtmp(j,i+1) - d2 qtmp(j,i) = qtmp(j,i) + d2 + tiny ENDIF END DO END DO i=1 DO j=j1,j2 IF(qtmp(j,i)<0.) THEN ipx = 1 ! west d0 = max(0.,qtmp(j,imr)) d1 = min(-qtmp(j,i),d0) qtmp(j,imr) = qtmp(j,imr) - d1 qtmp(j,i) = qtmp(j,i) + d1 ! east d0 = max(0.,qtmp(j,i+1)) d2 = min(-qtmp(j,i),d0) qtmp(j,i+1) = qtmp(j,i+1) - d2 qtmp(j,i) = qtmp(j,i) + d2 + tiny ENDIF END DO i=IMR DO j=j1,j2 IF(qtmp(j,i)<0.) THEN ipx = 1 ! west d0 = max(0.,qtmp(j,i-1)) d1 = min(-qtmp(j,i),d0) qtmp(j,i-1) = qtmp(j,i-1) - d1 qtmp(j,i) = qtmp(j,i) + d1 ! east d0 = max(0.,qtmp(j,1)) d2 = min(-qtmp(j,i),d0) qtmp(j,1) = qtmp(j,1) - d2 qtmp(j,i) = qtmp(j,i) + d2 + tiny ENDIF END DO IF(ipx/=0) THEN DO j=j1,j2 DO i=1,imr q(i,j) = qtmp(j,i) END DO END DO else ! Poles. IF(q(1,1)<0 .OR. q(1,JNP)<0.) ipx = 1 ENDIF RETURN END SUBROUTINE filew SUBROUTINE zflip(q,im,km,nc) IMPLICIT NONE ! This routine flip the array q (in the vertical). INTEGER :: im,km,nc REAL :: q(im,km,nc) ! local dynamic array REAL :: qtmp(im,km) INTEGER :: IC,k,i DO IC = 1, nc DO k=1,km DO i=1,im qtmp(i,k) = q(i,km+1-k,IC) END DO END DO DO i=1,im*km q(i,1,IC) = qtmp(i,1) END DO END DO RETURN END SUBROUTINE zflip