Changeset 5159 for LMDZ6/branches/Amaury_dev/libf/dyn3d_common/ppm3d.f90

Timestamp:

Aug 2, 2024, 9:58:25 PM (7 weeks ago)

Author:

abarral

Message:

Put dimensions.h and paramet.h into modules

File:

• LMDZ6/branches/Amaury_dev/libf/dyn3d_common/ppm3d.f90 (modified) (147 diffs)

LMDZ6/branches/Amaury_dev/libf/dyn3d_common/ppm3d.f90

@@ -16 +16 @@
 
 ! ********************************************************************
-!
+
 ! 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
@@ -28 +28 @@
 !      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.
@@ -38 +38 @@
 ! Subroutines modified: xtp, ytp, fzppm, qckxyz
 ! Subroutines deleted: vanz
-!
+
 ! Author: Shian-Jiann Lin
 ! mail address:
@@ -45 +45 @@
 !             Phone: 301-286-9540
 !             E-mail: lin@dao.gsfc.nasa.gov
-!
+
 ! *affiliation:
 !             Joint Center for Earth Systems Technology
@@ -51 +51 @@
 !             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)
@@ -72 +72 @@
   !  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
@@ -86 +86 @@
   !   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)
@@ -92 +92 @@
   ! 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
@@ -102 +102 @@
   ! 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)
@@ -117 +117 @@
   !           |
   !  W(1)    \ / ---------------------------------  AP(2), BP(2)
-  !
-  !
-  !
+
+
+
   ! W(k-1)   / \ ---------------------------------  AP(k), BP(k)
   !           |
@@ -125 +125 @@
   !           |
   !  W(k)    \ / ---------------------------------  AP(k+1), BP(k+1)
-  !
-  !
-  !
+
+
+
   !          / \ ---------------------------------  AP(NLAY), BP(NLAY)
   !           |
@@ -134 +134 @@
   !   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)
   !                                      |
@@ -154 +154 @@
   !                                      |
   !                                   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.
@@ -172 +172 @@
   ! 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
@@ -178 +178 @@
   ! 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
@@ -193 +193 @@
   !       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
@@ -204 +204 @@
   ! 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
@@ -227 +227 @@
   !          ( W > 0 is downward, ie, toward surface)
   ! PS2: predicted PS at t+NDT (original values are over-written)
-  !
+
   ! ********************************************************************
   ! NOTES:
@@ -235 +235 @@
   ! 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
@@ -259 +259 @@
   ! 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
   !
@@ -283 +283 @@
   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), &
@@ -291 +291 @@
         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)
@@ -302 +302 @@
   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
@@ -312 +312 @@
   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
@@ -338 +338 @@
   ENDIF
 
-  !
+
   IF(Jmax<JNP .OR. Kmax<NLAY) THEN
     WRITE(6,*) 'Jmax or Kmax is too small'
     stop
   ENDIF
-  !
+
   DO k=1,NLAY
   DAP(k) = (AP(k+1) - AP(k))*PT
   DBK(k) =  BP(k+1) - BP(k)
   ENDDO
-  !
+
   PI = 4. * ATAN(1.)
   DL = 2.*PI / REAL(IMR)
   DP =    PI / REAL(JMR)
-  !
+
   IF(IGD==0) THEN
   ! Compute analytic cosine at cell edges
@@ -360 +360 @@
         CALL cosc(cosp,cose,JNP,PI,DP)
   ENDIF
-  !
+
   DO J=2,JMR
   acosp(j) = 1. / cosp(j)
   END DO
-  !
+
   ! Inverse of the Scaled polar cap area.
-  !
+
   RCAP  = DP / (IMR*(1.-COS((j1-1.5)*DP)))
   acosp(1)   = RCAP
   acosp(JNP) = RCAP
   ENDIF
-  !
+
   IF(NDT0 /= NDT) THEN
   DT   = NDT
@@ -385 +385 @@
         WRITE(6,*) 'Warning!!! NDT maybe too large!'
   ENDIF
-  !
+
   IF(CR1>=0.95) THEN
   JS0 = 0
@@ -393 +393 @@
   else
   ZTC  = acos(CR1) * (180./PI)
-  !
+
   JS0 = REAL(JMR)*(90.-ZTC)/180. + 2
   JS0 = max(JS0, J1+1)
@@ -399 +399 @@
   JN0 = JNP-JS0+1
   ENDIF
-  !
-  !
+
+
   DO J=2,JMR
   DTDX(j)  = DT / ( DL*AE*COSP(J) )
@@ -412 +412 @@
    ! 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
@@ -428 +428 @@
   enddo
 
-  !
+
   IF(j1/=2) THEN
   DO IC=1,NC
@@ -439 +439 @@
   END DO
   ENDIF
-  !
+
   ! Compute "tracer density"
   DO IC=1,NC
@@ -450 +450 @@
   END DO
   END DO
-  !
+
   DO k=1,NLAY
-  !
+
   IF(IGD==0) THEN
   ! Convert winds on A-Grid to Courant # on C-Grid.
@@ -464 +464 @@
   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
@@ -486 +486 @@
   enddo
   enddo
-  !
+
 2222   continue
   DO j=JN0,j2-1
@@ -497 +497 @@
   enddo
 2233   continue
-  !
+
   IF(j1/=2) then           ! Enlarged polar cap.
   DO i=1,IMR
@@ -504 +504 @@
   enddo
   ENDIF
-  !
+
   ! ******* Compute horizontal mass fluxes ************
-  !
+
   ! N-S component
   DO j=j1,j2+1
@@ -514 +514 @@
   enddo
   enddo
-  !
+
   DO j=j1,j2
   DO i=1,IMR
@@ -520 +520 @@
   END DO
   END DO
-  !
+
   ! Poles
   sum1 = ymass(IMR,j1  )
@@ -528 +528 @@
   sum2 = sum2 + ymass(i,J2+1)
   enddo
-  !
+
   sum1 = - sum1 * RCAP
   sum2 =   sum2 * RCAP
@@ -535 +535 @@
   DPI(i,JNP,k) = sum2
   enddo
-  !
+
   ! E-W component
-  !
+
   DO j=j1,j2
   DO i=2,IMR
@@ -543 +543 @@
   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
@@ -553 +553 @@
   END DO
   END DO
-  !
+
   DO j=j1,j2
   DO i=1,IMR-1
@@ -559 +559 @@
   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
@@ -569 +569 @@
   enddo
   enddo
-  !
+
   DO j=j1,j2
   UA(imr,j) = 0.5 * (CRX(imr,j)+CRX(1,j))
@@ -585 +585 @@
   VA(i,2) = 0.5*(CRY(i,2)+CRY(i,3))
   enddo
-  !
+
   IF(j1==2) THEN
     IMH = IMR/2
@@ -597 +597 @@
   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<JN) GO TO 250
-  !
+
   DO i=1,IMR
   qtmp(i) = q(i,j,k,IC)
   enddo
-  !
+
   DO i=-IML,0
   qtmp(i)       = q(IMR+i,j,k,IC)
   qtmp(IMR+1-i) = q(1-i,j,k,IC)
   enddo
-  !
+
   DO i=1,IMR
   iu = UA(i,j)
@@ -632 +632 @@
 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)
@@ -655 +655 @@
   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.
@@ -668 +668 @@
             iad = 1
     endif
-  !
+
     IF(JORD >= 2) THEN
             jad = 2
@@ -682 +682 @@
   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)
@@ -688 +688 @@
   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
-  !
+
+  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
@@ -701 +701 @@
   END DO
   END DO
-  !
+
   DO k=2,NLAY
   DO j=1,JNP
@@ -709 +709 @@
   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
@@ -732 +732 @@
   END DO
   END DO
-  !
+
   DO k=1,NLAY
   DO j=1,JNP
@@ -740 +740 @@
   END DO
   END DO
-  !
+
     KRD = max(3, KORD)
   DO IC=1,NC
-  !
+
   !****6***0*********0*********0*********0*********0*********0**********72
 
@@ -752 +752 @@
   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
@@ -764 +764 @@
   enddo
   enddo
-  !
+
   IF(j1/=2) THEN
   DO k=1,NLAY
@@ -775 +775 @@
   ENDIF
   END DO
-  !
+
   IF(j1/=2) THEN
   DO k=1,NLAY
@@ -784 +784 @@
   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, &
@@ -803 +803 @@
   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
@@ -819 +819 @@
   END DO
   END DO
-  !
+
   DO k=2,NLAYM1
   DO I=1,IMJM
@@ -832 +832 @@
   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
@@ -848 +848 @@
   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
@@ -869 +869 @@
   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
@@ -878 +878 @@
     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)))
@@ -891 +891 @@
   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
@@ -908 +908 @@
   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
@@ -924 +924 @@
   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
@@ -944 +944 @@
   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
@@ -963 +963 @@
   RETURN
 END SUBROUTINE fzppm
-!
+
 SUBROUTINE xtp(IMR,JNP,IML,j1,j2,JN,JS,PU,DQ,Q,UC, &
         fx1,xmass,IORD)
@@ -976 +976 @@
   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
@@ -1006 +1006 @@
   CALL xmist(IMR,IML,Qtmp,DC)
   DC(0) = DC(IMR)
-  !
+
   IF(IORD==2 .OR. j<=j1vl .OR. j>=j2vl) THEN
   DO i=1,IMR
@@ -1015 +1015 @@
   CALL fxppm(IMR,IML,UC(1,j),Qtmp,DC,fx1,IORD)
   ENDIF
-  !
-  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
@@ -1042 +1042 @@
   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))
@@ -1056 +1056 @@
   END DO
   ENDIF
-  !
+
   DO i=1,IMR
   IF(uc(i,j)>1.) THEN
@@ -1073 +1073 @@
   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
@@ -1099 +1099 @@
    ! SAVE LMT
    ! IF(first) THEN
-  !
+
   ! correction calcul de LMT a chaque passage pour pouvoir choisir
   ! plusieurs schemas PPM pour differents traceurs
@@ -1113 +1113 @@
   !        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
@@ -1149 +1149 @@
   RETURN
 END SUBROUTINE fxppm
-!
+
 SUBROUTINE xmist(IMR,IML,P,DC)
   IMPLICIT NONE
@@ -1157 +1157 @@
   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))
@@ -1166 +1166 @@
   RETURN
 END SUBROUTINE xmist
-!
+
 SUBROUTINE ytp(IMR,JNP,j1,j2,acosp,RCAP,DQ,P,VC,DC2 &
         ,ymass,fx,A6,AR,AL,JORD)
@@ -1178 +1178 @@
   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
@@ -1190 +1190 @@
 
   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)
@@ -1201 +1201 @@
   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
@@ -1211 +1211 @@
   END DO
   END DO
-  !
+
   ! Poles
   sum1 = fx(IMR,j1  )
@@ -1219 +1219 @@
   sum2 = sum2 + fx(i,J2+1)
   enddo
-  !
+
   sum1 = DQ(1,  1) - sum1 * RCAP
   sum2 = DQ(1,JNP) + sum2 * RCAP
@@ -1226 +1226 @@
   DQ(i,JNP) = sum2
   enddo
-  !
+
   IF(j1/=2) THEN
   DO i=1,IMR
@@ -1233 +1233 @@
   enddo
   ENDIF
-  !
+
   RETURN
 END SUBROUTINE ytp
-!
+
 subroutine  ymist(IMR,JNP,j1,P,DC,ID)
   IMPLICIT NONE
@@ -1244 +1244 @@
   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)
@@ -1281 +1281 @@
   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
@@ -1289 +1289 @@
   END DO
   ENDIF
-  !
+
   IF(j1/=2) THEN
   DO i=1,IMR
@@ -1297 +1297 @@
   else
   ! Determine slopes in polar caps for scalars!
-  !
+
   DO i=1,IMH
   ! South
@@ -1310 +1310 @@
   DC(i,JNP) = sign(min(abs(tmp),Pmax,pmin),tmp)
   END DO
-  !
+
   DO i=imh+1,IMR
   DC(i,  1) =  - DC(i-imh,  1)
@@ -1318 +1318 @@
   RETURN
 END SUBROUTINE ymist
-!
+
 SUBROUTINE fyppm(VC,P,DC,flux,IMR,JNP,j1,j2,A6,AR,AL,JORD)
   IMPLICIT NONE
@@ -1331 +1331 @@
    ! data first /.TRUE./
    ! SAVE LMT
-  !
+
   IMH = IMR / 2
   JMR = JNP - 1
@@ -1349 +1349 @@
    !       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)
@@ -1382 +1382 @@
   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)
@@ -1398 +1398 @@
   RETURN
 END SUBROUTINE fyppm
-!
+
   SUBROUTINE yadv(IMR,JNP,j1,j2,p,VA,ady,wk,IAD)
     IMPLICIT NONE
@@ -1406 +1406 @@
     INTEGER :: JMR,IMH,i,j,jp
     REAL :: rv,a1,b1,sum1,sum2
-  !
+
     JMR = JNP-1
     IMH = IMR/2
@@ -1443 +1443 @@
   enddo
    ! WRITE(*,*) 'toto 2'
-  !
+
   ELSEIF(IAD==1) THEN
     DO j=j1-1,j2+1
@@ -1452 +1452 @@
   enddo
   ENDIF
-  !
+
     IF(j1/=2) THEN
     sum1 = 0.
@@ -1462 +1462 @@
     sum1 = sum1 / IMR
     sum2 = sum2 / IMR
-  !
+
   DO i=1,imr
   ady(i,  2) =  sum1
@@ -1479 +1479 @@
     sum1 = sum1 / IMR
     sum2 = sum2 / IMR
-  !
+
   DO i=1,imr
   ady(i,  1) =  sum1
@@ -1485 +1485 @@
   enddo
     endif
-  !
+
     RETURN
 END SUBROUTINE yadv
-!
+
   SUBROUTINE xadv(IMR,JNP,j1,j2,p,UA,JS,JN,IML,adx,IAD)
     IMPLICIT NONE
@@ -1495 +1495 @@
     INTEGER :: JMR,j,i,ip,iu,iiu
     REAL :: ru,a1,b1
-  !
+
     JMR = JNP-1
   DO j=j1,j2
   IF(J>JS  .AND. J<JN) GO TO 1309
-  !
+
   DO i=1,IMR
   qtmp(i) = p(i,j)
   enddo
-  !
+
   DO i=-IML,0
   qtmp(i)       = p(IMR+i,j)
   qtmp(IMR+1-i) = p(1-i,j)
   enddo
-  !
+
   IF(IAD==2) THEN
   DO i=1,IMR
@@ -1530 +1530 @@
   enddo
   ENDIF
-  !
+
   DO i=1,IMR
   adx(i,j) = adx(i,j) - p(i,j)
@@ -1536 +1536 @@
 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)
@@ -1567 +1567 @@
   ENDIF
   enddo
-  !
+
     IF(j1/=2) THEN
   DO i=1,IMR
@@ -1581 +1581 @@
     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
@@ -1591 +1591 @@
   ! 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)
@@ -1603 +1603 @@
   INTEGER :: i
   REAL :: da1,da2,a6da,fmin
-  !
+
   IF(LMT==0) THEN
   ! Full constraint
@@ -1662 +1662 @@
   RETURN
 END SUBROUTINE lmtppm
-!
+
 SUBROUTINE A2C(U,V,IMR,JMR,j1,j2,CRX,CRY,dtdx5,DTDY5)
   IMPLICIT NONE
@@ -1668 +1668 @@
   REAL :: U(IMR,*),V(IMR,*),CRX(IMR,*),CRY(IMR,*),DTDX5(*),DTDY5
   INTEGER :: i,j
-  !
+
   DO j=j1,j2
   DO i=2,IMR
@@ -1674 +1674 @@
   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))
@@ -1684 +1684 @@
   RETURN
 END SUBROUTINE a2c
-!
+
 SUBROUTINE cosa(cosp,cose,JNP,PI,DP)
   IMPLICIT NONE
@@ -1696 +1696 @@
     cose(j) = cos(ph5)
   END DO
-  !
+
   JEQ = (JNP+1) / 2
   IF(JMR == 2*(JMR/2) ) THEN
@@ -1709 +1709 @@
    enddo
   ENDIF
-  !
+
   DO j=2,JMR
   cosp(j) = 0.5*(cose(j)+cose(j+1))
@@ -1717 +1717 @@
   RETURN
 END SUBROUTINE cosa
-!
+
 SUBROUTINE cosc(cosp,cose,JNP,PI,DP)
   IMPLICIT NONE
@@ -1724 +1724 @@
   REAL :: phi
   INTEGER :: j
-  !
+
   phi = -0.5*PI
   DO j=2,JNP-1
@@ -1732 +1732 @@
     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))
@@ -1742 +1742 @@
   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
@@ -1752 +1752 @@
   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
@@ -1764 +1764 @@
   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
@@ -1778 +1778 @@
   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
@@ -1791 +1791 @@
   ENDIF
   IF(icr==0) goto 225
-  !
+
   DO i=1,len
   IF( Q(I,j1,L)<0.) THEN
-  !
+
   ip = ip + 1
   ! From above
@@ -1809 +1809 @@
 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)
@@ -1837 +1837 @@
    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
@@ -1850 +1850 @@
   RETURN
 END SUBROUTINE qckxyz
-!
+
 SUBROUTINE filcr(q,IMR,JNP,j1,j2,cosp,acosp,icr,tiny)
   IMPLICIT NONE
@@ -1937 +1937 @@
 65   continue
   END DO
-  !
+
   DO i=1,IMR
   IF(q(i,j1)<0. .OR. q(i,j2)<0.) THEN
@@ -1944 +1944 @@
   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
@@ -1963 +1963 @@
    ! data first /.TRUE./
    ! save cap1
-  !
+
   !  IF(first) THEN
   DP = 4.*ATAN(1.)/REAL(JNP-1)
@@ -1969 +1969 @@
    ! first = .FALSE.
    ! END IF
-  !
+
   ipy = 0
   DO j=j1+1,j2-1
@@ -1991 +1991 @@
   END DO
   END DO
-  !
+
   DO i=1,imr
   IF(q(i,j1)<0.) THEN
@@ -2004 +2004 @@
   ENDIF
   enddo
-  !
+
   j = j2
   DO i=1,imr
@@ -2018 +2018 @@
   ENDIF
   enddo
-  !
+
   ! Check Poles.
   IF(q(1,1)<0.) THEN
@@ -2028 +2028 @@
   enddo
   ENDIF
-  !
+
   IF(q(1,JNP)<0.) THEN
   dq = q(1,JNP)*cap1/REAL(IMR)*acosp(j2)
@@ -2037 +2037 @@
   enddo
   ENDIF
-  !
+
   RETURN
 END SUBROUTINE filns
-!
+
 SUBROUTINE filew(q,qtmp,IMR,JNP,j1,j2,ipx,tiny)
   IMPLICIT NONE
@@ -2047 +2047 @@
   INTEGER :: i,j
   REAL :: d0,d1,d2
-  !
+
   ipx = 0
   ! Copy & swap direction for vectorization.
@@ -2055 +2055 @@
   END DO
   END DO
-  !
+
   DO i=2,imr-1
   DO j=j1,j2
@@ -2073 +2073 @@
   END DO
   END DO
-  !
+
   i=1
   DO j=j1,j2
@@ -2087 +2087 @@
   d2 = min(-qtmp(j,i),d0)
   qtmp(j,i+1) = qtmp(j,i+1) - d2
-  !
+
   qtmp(j,i) = qtmp(j,i) + d2 + tiny
   ENDIF
@@ -2104 +2104 @@
   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
@@ -2116 +2116 @@
   END DO
   else
-  !
+
   ! Poles.
   IF(q(1,1)<0 .OR. q(1,JNP)<0.) ipx = 1
@@ -2122 +2122 @@
   RETURN
 END SUBROUTINE filew
-!
+
 SUBROUTINE zflip(q,im,km,nc)
   IMPLICIT NONE
@@ -2131 +2131 @@
   REAL :: qtmp(im,km)
   INTEGER :: IC,k,i
-  !
+
   DO IC = 1, nc
-  !
+
   DO k=1,km
   DO i=1,im
@@ -2139 +2139 @@
   END DO
   END DO
-  !
+
   DO i=1,im*km
   q(i,1,IC) = qtmp(i,1)