dyn3d/PPM3d/ppm3d.f

module ppm3d_m

  implicit none

contains

  SUBROUTINE ppm3d(igd, q, ps1, ps2, u, v, w, ndt, iord, jord, kord, nc, imr, &
       jnp, j1, nlay, ap, bp, pt, ae, fill, umax)

    ! 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 number of constituents
    ! IMR: first dimension (E-W); number of Grid intervals in E-W is IMR
    ! JNP: 2nd dimension (N-S); number of Grid intervals in N-S is JNP-1
    ! NLAY: 3rd dimension (number 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 number 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

    integer imr
    INTEGER igd, iord, jord, kord, nc, jnp, j1, nlay, ae
    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 imrd2
    REAL pt
    LOGICAL cross, fill

    ! 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)
    integer ndt0, nstep
    DATA ndt0, nstep/0, 0/
    DATA cross/.TRUE./
    SAVE dtdy, dtdy5, rcap, js0, jn0, iml, dtdx, dtdx5, acosp, cosp, cose, dap, dbk
    real cr1, d5, dl, dp, dt, dtdy, dtdy5, pi, rcap, ru, sum1, sum2, ztc
    integer i, iad, ic, iiu, imh, imjm, iml, iu, j, j2, jad, jmp, jmr, jn, jn0
    integer js, js0, jt, k, krd, l, maxdt

    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
       END IF
       IF (jnp<nlay) THEN
          WRITE (6, *) 'JNP must be >= NLAY'
          STOP
       END IF
       imrd2 = mod(imr, 2)
       IF (j1==2 .AND. imrd2/=0) THEN
          WRITE (6, *) 'if j1=2 IMR must be an even integer'
          STOP
       END IF


       IF (jmax<jnp .OR. kmax<nlay) THEN
          WRITE (6, *) 'Jmax or Kmax is too small'
          STOP
       END IF

       DO k = 1, nlay
          dap(k) = (ap(k+1)-ap(k))*pt
          dbk(k) = bp(k+1) - bp(k)
       END DO

       pi = 4.*atan(1.)
       dl = 2.*pi/float(imr)
       dp = pi/float(jmr)

       IF (igd==0) THEN
          ! Compute analytic cosine at cell edges
          CALL cosa(cosp, cose, jnp, pi, dp)
       ELSE
          ! Define cosine consistent with GEOS-GCM (using dycore2.0 or later)
          CALL cosc(cosp, cose, jnp, pi, dp)
       END IF

       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
    END IF

    IF (ndt0/=ndt) THEN
       dt = ndt
       ndt0 = ndt

       IF (umax<180.) THEN
          WRITE (6, *) 'Umax may be too small!'
       END IF
       cr1 = abs(umax*dt)/(dl*ae)
       maxdt = dp*ae/abs(umax) + 0.5
       WRITE (6, *) 'Largest time step for max(V)=', umax, ' is ', maxdt
       IF (maxdt<abs(ndt)) THEN
          WRITE (6, *) 'Warning!!! NDT maybe too large!'
       END IF

       IF (cr1>=0.95) THEN
          js0 = 0
          jn0 = 0
          iml = imr - 2
          ztc = 0.
       ELSE
          ztc = acos(cr1)*(180./pi)

          js0 = float(jmr)*(90.-ztc)/180. + 2
          js0 = max(js0, j1+1)
          iml = min(6*js0/(j1-1)+2, 4*imr/5)
          jn0 = jnp - js0 + 1
       END IF


       DO j = 2, jmr
          dtdx(j) = dt/(dl*ae*cosp(j))

          dtdx5(j) = 0.5*dtdx(j)
       END DO


       dtdy = dt/(ae*dp)
       dtdy5 = 0.5*dtdy

    END IF

    ! *********** 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)
          END DO
       END DO
    END DO


    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
    END IF

    ! 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 number 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 number
          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
       END IF

       ! 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
             END IF
          END DO
       END DO

2222   CONTINUE
       DO j = jn0, j2 - 1
          DO i = 1, imr
             IF (abs(crx(i,j))>1.) THEN
                jn = j
                GO TO 2233
             END IF
          END DO
       END DO
2233   CONTINUE

       IF (j1/=2) THEN ! Enlarged polar cap.
          DO i = 1, imr
             dpi(i, 2, k) = 0.
             dpi(i, jmr, k) = 0.
          END DO
       END IF

       ! ******* 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))
          END DO
       END DO

       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)
       END DO

       sum1 = -sum1*rcap
       sum2 = sum2*rcap
       DO i = 1, imr
          dpi(i, 1, k) = sum1
          dpi(i, jnp, k) = sum2
       END DO

       ! 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))
          END DO
       END DO

       DO j = j1, j2
          pu(1, j) = 0.5*(delp2(1,j,k)+delp2(imr,j,k))
       END DO

       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))
          END DO
       END DO

       DO j = j1, j2
          ua(imr, j) = 0.5*(crx(imr,j)+crx(1,j))
       END DO
       ! cccccccccccccccccccccccccccccccccccccccccccccccccccccc
       ! Rajouts pour LMDZ.3.3
       ! cccccccccccccccccccccccccccccccccccccccccccccccccccccc
       DO i = 1, imr
          DO j = 1, jnp
             va(i, j) = 0.
          END DO
       END DO

       DO i = 1, imr*(jmr-1)
          va(i, 2) = 0.5*(cry(i,2)+cry(i,3))
       END DO

       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)
          END DO
          va(imr, 1) = va(1, 1)
          va(imr, jnp) = va(1, jnp)
       END IF

       ! ****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.
          END DO

          ! E-W advective cross term
          DO j = j1, j2
             IF (j>js .AND. j<jn) cycle

             DO i = 1, imr
                qtmp(i) = q(i, j, k, ic)
             END DO

             DO i = -iml, 0
                qtmp(i) = q(imr+i, j, k, ic)
                qtmp(imr+1-i) = q(1-i, j, k, ic)
             END DO

             DO i = 1, imr
                iu = ua(i, j)
                ru = ua(i, j) - iu
                iiu = i - iu
                IF (ua(i,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))
                END IF
                wk1(i, j, 1) = wk1(i, j, 1) - qtmp(i)
             END DO
          END DO

          IF (jn/=0) THEN
             DO j = js + 1, jn - 1

                DO i = 1, imr
                   qtmp(i) = q(i, j, k, ic)
                END DO

                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))
                END DO
             END DO
          END IF
          ! ****6***0*********0*********0*********0*********0*********0**********72
          ! Contribution from the N-S advection
          DO i = 1, imr*(j2-j1+1)
             jt = float(j1) - va(i, j1)
             wk1(i, j1, 2) = va(i, j1)*(q(i,jt,k,ic)-q(i,jt+1,k,ic))
          END DO

          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)
          END DO

          IF (cross) THEN
             ! Add cross terms in the vertical direction.
             IF (iord>=2) THEN
                iad = 2
             ELSE
                iad = 1
             END IF

             IF (jord>=2) THEN
                jad = 2
             ELSE
                jad = 1
             END IF
             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)
                END DO
             END DO
          END IF

          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)
             END DO
          END DO
       END DO

       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, jmp, k, ic)
             END DO
          END DO
       END IF
    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
    END IF

  END SUBROUTINE ppm3d

end module ppm3d_m
1	module ppm3d_m
2
3	implicit none
4
5	contains
6
7	SUBROUTINE ppm3d(igd, q, ps1, ps2, u, v, w, ndt, iord, jord, kord, nc, imr, &
8	jnp, j1, nlay, ap, bp, pt, ae, fill, umax)
9
10	! rajout de déclarations
11	! integer Jmax,kmax,ndt0,nstep,k,j,i,ic,l,js,jn,imh,iad,jad,krd
12	! integer iu,iiu,j2,jmr,js0,jt
13	! real dtdy,dtdy5,rcap,iml,jn0,imjm,pi,dl,dp
14	! real dt,cr1,maxdt,ztc,d5,sum1,sum2,ru
15
16	! ********************************************************************
17
18	! =============
19	! INPUT:
20	! =============
21
22	! Q(IMR,JNP,NLAY,NC): mixing ratios at current time (t)
23	! NC: total number of constituents
24	! IMR: first dimension (E-W); number of Grid intervals in E-W is IMR
25	! JNP: 2nd dimension (N-S); number of Grid intervals in N-S is JNP-1
26	! NLAY: 3rd dimension (number of layers); vertical index increases from 1
27	! at
28	! the model top to NLAY near the surface (see fig. below).
29	! It is assumed that 6 <= NLAY <= JNP (for dynamic memory allocation)
30
31	! PS1(IMR,JNP): surface pressure at current time (t)
32	! PS2(IMR,JNP): surface pressure at mid-time-level (t+NDT/2)
33	! PS2 is replaced by the predicted PS (at t+NDT) on output.
34	! Note: surface pressure can have any unit or can be multiplied by any
35	! const.
36
37	! The pressure at layer edges are defined as follows:
38
39	! p(i,j,k) = AP(k)PT + BP(k)PS(i,j) (1)
40
41	! Where PT is a constant having the same unit as PS.
42	! AP and BP are unitless constants given at layer edges
43	! defining the vertical coordinate.
44	! BP(1) = 0., BP(NLAY+1) = 1.
45	! The pressure at the model top is PTOP = AP(1)*PT
46
47	! For pure sigma system set AP(k) = 1 for all k, PT = PTOP,
48	! BP(k) = sige(k) (sigma at edges), PS = Psfc - PTOP.
49
50	! Note: the sigma-P coordinate is a subset of Eq. 1, which in turn
51	! is a subset of the following even more general sigma-P-thelta coord.
52	! currently under development.
53	! p(i,j,k) = (AP(k)PT + BP(k)PS(i,j))/(D(k)-C(k)TE*(-1/kapa))
54
55	! /////////////////////////////////
56	! / \ ------------- PTOP -------------- AP(1), BP(1)
57	! \|
58	! delp(1) \| ........... Q(i,j,1) ............
59	! \|
60	! W(1) \ / --------------------------------- AP(2), BP(2)
61
62
63
64	! W(k-1) / \ --------------------------------- AP(k), BP(k)
65	! \|
66	! delp(K) \| ........... Q(i,j,k) ............
67	! \|
68	! W(k) \ / --------------------------------- AP(k+1), BP(k+1)
69
70
71
72	! / \ --------------------------------- AP(NLAY), BP(NLAY)
73	! \|
74	! delp(NLAY) \| ........... Q(i,j,NLAY) .........
75	! \|
76	! W(NLAY)=0 \ / ------------- surface ----------- AP(NLAY+1), BP(NLAY+1)
77	! //////////////////////////////////
78
79	! U(IMR,JNP,NLAY) & V(IMR,JNP,NLAY):winds (m/s) at mid-time-level (t+NDT/2)
80	! U and V may need to be polar filtered in advance in some cases.
81
82	! IGD: grid type on which winds are defined.
83	! IGD = 0: A-Grid [all variables defined at the same point from south
84	! pole (j=1) to north pole (j=JNP) ]
85
86	! IGD = 1 GEOS-GCM C-Grid
87	! [North]
88
89	! V(i,j)
90	! \|
91	! \|
92	! \|
93	! U(i-1,j)---Q(i,j)---U(i,j) [EAST]
94	! \|
95	! \|
96	! \|
97	! V(i,j-1)
98
99	! U(i, 1) is defined at South Pole.
100	! V(i, 1) is half grid north of the South Pole.
101	! V(i,JMR) is half grid south of the North Pole.
102
103	! V must be defined at j=1 and j=JMR if IGD=1
104	! V at JNP need not be given.
105
106	! NDT: time step in seconds (need not be constant during the course of
107	! the integration). Suggested value: 30 min. for 4x5, 15 min. for 2x2.5
108	! (Lat-Lon) resolution. Smaller values are recommanded if the model
109	! has a well-resolved stratosphere.
110
111	! J1 defines the size of the polar cap:
112	! South polar cap edge is located at -90 + (j1-1.5)*180/(JNP-1) deg.
113	! North polar cap edge is located at 90 - (j1-1.5)*180/(JNP-1) deg.
114	! There are currently only two choices (j1=2 or 3).
115	! IMR must be an even integer if j1 = 2. Recommended value: J1=3.
116
117	! IORD, JORD, and KORD are integers controlling various options in E-W,
118	! N-S,
119	! and vertical transport, respectively. Recommended values for positive
120	! definite scalars: IORD=JORD=3, KORD=5. Use KORD=3 for non-
121	! positive definite scalars or when linear correlation between constituents
122	! is to be maintained.
123
124	! _ORD=
125	! 1: 1st order upstream scheme (too diffusive, not a useful option; it
126	! can be used for debugging purposes; this is THE only known "linear"
127	! monotonic advection scheme.).
128	! 2: 2nd order van Leer (full monotonicity constraint;
129	! see Lin et al 1994, MWR)
130	! 3: monotonic PPM* (slightly improved PPM of Collela & Woodward 1984)
131	! 4: semi-monotonic PPM (same as 3, but overshoots are allowed)
132	! 5: positive-definite PPM (constraint on the subgrid distribution is
133	! only strong enough to prevent generation of negative values;
134	! both overshoots & undershoots are possible).
135	! 6: un-constrained PPM (nearly diffusion free; slightly faster but
136	! positivity not quaranteed. Use this option only when the fields
137	! and winds are very smooth).
138
139	! *PPM: Piece-wise Parabolic Method
140
141	! Note that KORD <=2 options are no longer supported. DO not use option 4
142	! or 5.
143	! for non-positive definite scalars (such as Ertel Potential Vorticity).
144
145	! The implicit numerical diffusion decreases as _ORD increases.
146	! The last two options (ORDER=5, 6) should only be used when there is
147	! significant explicit diffusion (such as a turbulence parameterization).
148	! You
149	! might get dispersive results otherwise.
150	! No filter of any kind is applied to the constituent fields here.
151
152	! AE: Radius of the sphere (meters).
153	! Recommended value for the planet earth: 6.371E6
154
155	! fill(logical): flag to do filling for negatives (see note below).
156
157	! Umax: Estimate (upper limit) of the maximum U-wind speed (m/s).
158	! (220 m/s is a good value for troposphere model; 280 m/s otherwise)
159
160	! =============
161	! Output
162	! =============
163
164	! Q: mixing ratios at future time (t+NDT) (original values are
165	! over-written)
166	! W(NLAY): large-scale vertical mass flux as diagnosed from the hydrostatic
167	! relationship. W will have the same unit as PS1 and PS2 (eg, mb).
168	! W must be divided by NDT to get the correct mass-flux unit.
169	! The vertical Courant number C = W/delp_UPWIND, where delp_UPWIND
170	! is the pressure thickness in the "upwind" direction. For example,
171	! C(k) = W(k)/delp(k) if W(k) > 0;
172	! C(k) = W(k)/delp(k+1) if W(k) < 0.
173	! ( W > 0 is downward, ie, toward surface)
174	! PS2: predicted PS at t+NDT (original values are over-written)
175
176	! ********************************************************************
177	! NOTES:
178	! This forward-in-time upstream-biased transport scheme reduces to
179	! the 2nd order center-in-time center-in-space mass continuity eqn.
180	! if Q = 1 (constant fields will remain constant). This also ensures
181	! that the computed vertical velocity to be identical to GEOS-1 GCM
182	! for on-line transport.
183
184	! A larger polar cap is used if j1=3 (recommended for C-Grid winds or when
185	! winds are noisy near poles).
186
187	! Flux-Form Semi-Lagrangian transport in the East-West direction is used
188	! when and where Courant number is greater than one.
189
190	! The user needs to change the parameter Jmax or Kmax if the resolution
191	! is greater than 0.5 deg in N-S or 150 layers in the vertical direction.
192	! (this TransPort Core is otherwise resolution independent and can be used
193	! as a library routine).
194
195	! PPM is 4th order accurate when grid spacing is uniform (x & y); 3rd
196	! order accurate for non-uniform grid (vertical sigma coord.).
197
198	! Time step is limitted only by transport in the meridional direction.
199	! (the FFSL scheme is not implemented in the meridional direction).
200
201	! Since only 1-D limiters are applied, negative values could
202	! potentially be generated when large time step is used and when the
203	! initial fields contain discontinuities.
204	! This does not necessarily imply the integration is unstable.
205	! These negatives are typically very small. A filling algorithm is
206	! activated if the user set "fill" to be true.
207
208	! The van Leer scheme used here is nearly as accurate as the original PPM
209	! due to the use of a 4th order accurate reference slope. The PPM imple-
210	! mented here is an improvement over the original and is also based on
211	! the 4th order reference slope.
212
213	! **60******0*****0*****0*****0*****0********72
214
215	! User modifiable parameters
216
217	integer, PARAMETER:: jmax=361, kmax=150
218
219	! **60******0*****0*****0*****0*****0********72
220
221	! Input-Output arrays
222
223	integer imr
224	INTEGER igd, iord, jord, kord, nc, jnp, j1, nlay, ae
225	REAL q(imr, jnp, nlay, nc), ps1(imr, jnp), ps2(imr, jnp), &
226	u(imr, jnp, nlay), v(imr, jnp, nlay), ap(nlay+1), bp(nlay+1), &
227	w(imr, jnp, nlay), ndt, val(nlay), umax
228	INTEGER imrd2
229	REAL pt
230	LOGICAL cross, fill
231
232	! Local dynamic arrays
233
234	REAL crx(imr, jnp), cry(imr, jnp), xmass(imr, jnp), ymass(imr, jnp), &
235	fx1(imr+1), dpi(imr, jnp, nlay), delp1(imr, jnp, nlay), &
236	wk1(imr, jnp, nlay), pu(imr, jnp), pv(imr, jnp), dc2(imr, jnp), &
237	delp2(imr, jnp, nlay), dq(imr, jnp, nlay, nc), va(imr, jnp), &
238	ua(imr, jnp), qtmp(-imr:2*imr)
239
240	! Local static arrays
241
242	REAL dtdx(jmax), dtdx5(jmax), acosp(jmax), cosp(jmax), cose(jmax), &
243	dap(kmax), dbk(kmax)
244	integer ndt0, nstep
245	DATA ndt0, nstep/0, 0/
246	DATA cross/.TRUE./
247	SAVE dtdy, dtdy5, rcap, js0, jn0, iml, dtdx, dtdx5, acosp, cosp, cose, dap, dbk
248	real cr1, d5, dl, dp, dt, dtdy, dtdy5, pi, rcap, ru, sum1, sum2, ztc
249	integer i, iad, ic, iiu, imh, imjm, iml, iu, j, j2, jad, jmp, jmr, jn, jn0
250	integer js, js0, jt, k, krd, l, maxdt
251
252	jmr = jnp - 1
253	imjm = imr*jnp
254	j2 = jnp - j1 + 1
255	nstep = nstep + 1
256
257	! ********* Initialization ********************
258	IF (nstep==1) THEN
259
260	WRITE (6, *) '------------------------------------ '
261	WRITE (6, *) 'NASA/GSFC Transport Core Version 4.5'
262	WRITE (6, *) '------------------------------------ '
263
264	WRITE (6, *) 'IMR=', imr, ' JNP=', jnp, ' NLAY=', nlay, ' j1=', j1
265	WRITE (6, *) 'NC=', nc, iord, jord, kord, ndt
266
267	! controles sur les parametres
268	IF (nlay<6) THEN
269	WRITE (6, *) 'NLAY must be >= 6'
270	STOP
271	END IF
272	IF (jnp<nlay) THEN
273	WRITE (6, *) 'JNP must be >= NLAY'
274	STOP
275	END IF
276	imrd2 = mod(imr, 2)
277	IF (j1==2 .AND. imrd2/=0) THEN
278	WRITE (6, *) 'if j1=2 IMR must be an even integer'
279	STOP
280	END IF
281
282
283	IF (jmax<jnp .OR. kmax<nlay) THEN
284	WRITE (6, *) 'Jmax or Kmax is too small'
285	STOP
286	END IF
287
288	DO k = 1, nlay
289	dap(k) = (ap(k+1)-ap(k))*pt
290	dbk(k) = bp(k+1) - bp(k)
291	END DO
292
293	pi = 4.*atan(1.)
294	dl = 2.*pi/float(imr)
295	dp = pi/float(jmr)
296
297	IF (igd==0) THEN
298	! Compute analytic cosine at cell edges
299	CALL cosa(cosp, cose, jnp, pi, dp)
300	ELSE
301	! Define cosine consistent with GEOS-GCM (using dycore2.0 or later)
302	CALL cosc(cosp, cose, jnp, pi, dp)
303	END IF
304
305	DO j = 2, jmr
306	acosp(j) = 1./cosp(j)
307	END DO
308
309	! Inverse of the Scaled polar cap area.
310
311	rcap = dp/(imr(1.-cos((j1-1.5)dp)))
312	acosp(1) = rcap
313	acosp(jnp) = rcap
314	END IF
315
316	IF (ndt0/=ndt) THEN
317	dt = ndt
318	ndt0 = ndt
319
320	IF (umax<180.) THEN
321	WRITE (6, *) 'Umax may be too small!'
322	END IF
323	cr1 = abs(umaxdt)/(dlae)
324	maxdt = dp*ae/abs(umax) + 0.5
325	WRITE (6, *) 'Largest time step for max(V)=', umax, ' is ', maxdt
326	IF (maxdt<abs(ndt)) THEN
327	WRITE (6, *) 'Warning!!! NDT maybe too large!'
328	END IF
329
330	IF (cr1>=0.95) THEN
331	js0 = 0
332	jn0 = 0
333	iml = imr - 2
334	ztc = 0.
335	ELSE
336	ztc = acos(cr1)*(180./pi)
337
338	js0 = float(jmr)*(90.-ztc)/180. + 2
339	js0 = max(js0, j1+1)
340	iml = min(6js0/(j1-1)+2, 4imr/5)
341	jn0 = jnp - js0 + 1
342	END IF
343
344
345	DO j = 2, jmr
346	dtdx(j) = dt/(dlaecosp(j))
347
348	dtdx5(j) = 0.5*dtdx(j)
349	END DO
350
351
352	dtdy = dt/(ae*dp)
353	dtdy5 = 0.5*dtdy
354
355	END IF
356
357	! ********* End Initialization ********************
358
359	! delp = pressure thickness: the psudo-density in a hydrostatic system.
360	DO k = 1, nlay
361	DO j = 1, jnp
362	DO i = 1, imr
363	delp1(i, j, k) = dap(k) + dbk(k)*ps1(i, j)
364	delp2(i, j, k) = dap(k) + dbk(k)*ps2(i, j)
365	END DO
366	END DO
367	END DO
368
369
370	IF (j1/=2) THEN
371	DO ic = 1, nc
372	DO l = 1, nlay
373	DO i = 1, imr
374	q(i, 2, l, ic) = q(i, 1, l, ic)
375	q(i, jmr, l, ic) = q(i, jnp, l, ic)
376	END DO
377	END DO
378	END DO
379	END IF
380
381	! Compute "tracer density"
382	DO ic = 1, nc
383	DO k = 1, nlay
384	DO j = 1, jnp
385	DO i = 1, imr
386	dq(i, j, k, ic) = q(i, j, k, ic)*delp1(i, j, k)
387	END DO
388	END DO
389	END DO
390	END DO
391
392	DO k = 1, nlay
393
394	IF (igd==0) THEN
395	! Convert winds on A-Grid to Courant number on C-Grid.
396	CALL a2c(u(1,1,k), v(1,1,k), imr, jmr, j1, j2, crx, cry, dtdx5, dtdy5)
397	ELSE
398	! Convert winds on C-grid to Courant number
399	DO j = j1, j2
400	DO i = 2, imr
401	crx(i, j) = dtdx(j)*u(i-1, j, k)
402	END DO
403	END DO
404
405
406	DO j = j1, j2
407	crx(1, j) = dtdx(j)*u(imr, j, k)
408	END DO
409
410	DO i = 1, imr*jmr
411	cry(i, 2) = dtdy*v(i, 1, k)
412	END DO
413	END IF
414
415	! Determine JS and JN
416	js = j1
417	jn = j2
418
419	DO j = js0, j1 + 1, -1
420	DO i = 1, imr
421	IF (abs(crx(i,j))>1.) THEN
422	js = j
423	GO TO 2222
424	END IF
425	END DO
426	END DO
427
428	2222 CONTINUE
429	DO j = jn0, j2 - 1
430	DO i = 1, imr
431	IF (abs(crx(i,j))>1.) THEN
432	jn = j
433	GO TO 2233
434	END IF
435	END DO
436	END DO
437	2233 CONTINUE
438
439	IF (j1/=2) THEN ! Enlarged polar cap.
440	DO i = 1, imr
441	dpi(i, 2, k) = 0.
442	dpi(i, jmr, k) = 0.
443	END DO
444	END IF
445
446	! ***** Compute horizontal mass fluxes **********
447
448	! N-S component
449	DO j = j1, j2 + 1
450	d5 = 0.5*cose(j)
451	DO i = 1, imr
452	ymass(i, j) = cry(i, j)d5(delp2(i,j,k)+delp2(i,j-1,k))
453	END DO
454	END DO
455
456	DO j = j1, j2
457	DO i = 1, imr
458	dpi(i, j, k) = (ymass(i,j)-ymass(i,j+1))*acosp(j)
459	END DO
460	END DO
461
462	! Poles
463	sum1 = ymass(imr, j1)
464	sum2 = ymass(imr, j2+1)
465	DO i = 1, imr - 1
466	sum1 = sum1 + ymass(i, j1)
467	sum2 = sum2 + ymass(i, j2+1)
468	END DO
469
470	sum1 = -sum1*rcap
471	sum2 = sum2*rcap
472	DO i = 1, imr
473	dpi(i, 1, k) = sum1
474	dpi(i, jnp, k) = sum2
475	END DO
476
477	! E-W component
478
479	DO j = j1, j2
480	DO i = 2, imr
481	pu(i, j) = 0.5*(delp2(i,j,k)+delp2(i-1,j,k))
482	END DO
483	END DO
484
485	DO j = j1, j2
486	pu(1, j) = 0.5*(delp2(1,j,k)+delp2(imr,j,k))
487	END DO
488
489	DO j = j1, j2
490	DO i = 1, imr
491	xmass(i, j) = pu(i, j)*crx(i, j)
492	END DO
493	END DO
494
495	DO j = j1, j2
496	DO i = 1, imr - 1
497	dpi(i, j, k) = dpi(i, j, k) + xmass(i, j) - xmass(i+1, j)
498	END DO
499	END DO
500
501	DO j = j1, j2
502	dpi(imr, j, k) = dpi(imr, j, k) + xmass(imr, j) - xmass(1, j)
503	END DO
504
505	DO j = j1, j2
506	DO i = 1, imr - 1
507	ua(i, j) = 0.5*(crx(i,j)+crx(i+1,j))
508	END DO
509	END DO
510
511	DO j = j1, j2
512	ua(imr, j) = 0.5*(crx(imr,j)+crx(1,j))
513	END DO
514	! cccccccccccccccccccccccccccccccccccccccccccccccccccccc
515	! Rajouts pour LMDZ.3.3
516	! cccccccccccccccccccccccccccccccccccccccccccccccccccccc
517	DO i = 1, imr
518	DO j = 1, jnp
519	va(i, j) = 0.
520	END DO
521	END DO
522
523	DO i = 1, imr*(jmr-1)
524	va(i, 2) = 0.5*(cry(i,2)+cry(i,3))
525	END DO
526
527	IF (j1==2) THEN
528	imh = imr/2
529	DO i = 1, imh
530	va(i, 1) = 0.5*(cry(i,2)-cry(i+imh,2))
531	va(i+imh, 1) = -va(i, 1)
532	va(i, jnp) = 0.5*(cry(i,jnp)-cry(i+imh,jmr))
533	va(i+imh, jnp) = -va(i, jnp)
534	END DO
535	va(imr, 1) = va(1, 1)
536	va(imr, jnp) = va(1, jnp)
537	END IF
538
539	! **60******0*****0*****0*****0*****0********72
540	DO ic = 1, nc
541
542	DO i = 1, imjm
543	wk1(i, 1, 1) = 0.
544	wk1(i, 1, 2) = 0.
545	END DO
546
547	! E-W advective cross term
548	DO j = j1, j2
549	IF (j>js .AND. j<jn) cycle
550
551	DO i = 1, imr
552	qtmp(i) = q(i, j, k, ic)
553	END DO
554
555	DO i = -iml, 0
556	qtmp(i) = q(imr+i, j, k, ic)
557	qtmp(imr+1-i) = q(1-i, j, k, ic)
558	END DO
559
560	DO i = 1, imr
561	iu = ua(i, j)
562	ru = ua(i, j) - iu
563	iiu = i - iu
564	IF (ua(i,j)>=0.) THEN
565	wk1(i, j, 1) = qtmp(iiu) + ru*(qtmp(iiu-1)-qtmp(iiu))
566	ELSE
567	wk1(i, j, 1) = qtmp(iiu) + ru*(qtmp(iiu)-qtmp(iiu+1))
568	END IF
569	wk1(i, j, 1) = wk1(i, j, 1) - qtmp(i)
570	END DO
571	END DO
572
573	IF (jn/=0) THEN
574	DO j = js + 1, jn - 1
575
576	DO i = 1, imr
577	qtmp(i) = q(i, j, k, ic)
578	END DO
579
580	qtmp(0) = q(imr, j, k, ic)
581	qtmp(imr+1) = q(1, j, k, ic)
582
583	DO i = 1, imr
584	iu = i - ua(i, j)
585	wk1(i, j, 1) = ua(i, j)*(qtmp(iu)-qtmp(iu+1))
586	END DO
587	END DO
588	END IF
589	! **60******0*****0*****0*****0*****0********72
590	! Contribution from the N-S advection
591	DO i = 1, imr*(j2-j1+1)
592	jt = float(j1) - va(i, j1)
593	wk1(i, j1, 2) = va(i, j1)*(q(i,jt,k,ic)-q(i,jt+1,k,ic))
594	END DO
595
596	DO i = 1, imjm
597	wk1(i, 1, 1) = q(i, 1, k, ic) + 0.5*wk1(i, 1, 1)
598	wk1(i, 1, 2) = q(i, 1, k, ic) + 0.5*wk1(i, 1, 2)
599	END DO
600
601	IF (cross) THEN
602	! Add cross terms in the vertical direction.
603	IF (iord>=2) THEN
604	iad = 2
605	ELSE
606	iad = 1
607	END IF
608
609	IF (jord>=2) THEN
610	jad = 2
611	ELSE
612	jad = 1
613	END IF
614	CALL xadv(imr, jnp, j1, j2, wk1(1,1,2), ua, js, jn, iml, dc2, iad)
615	CALL yadv(imr, jnp, j1, j2, wk1(1,1,1), va, pv, w, jad)
616	DO j = 1, jnp
617	DO i = 1, imr
618	q(i, j, k, ic) = q(i, j, k, ic) + dc2(i, j) + pv(i, j)
619	END DO
620	END DO
621	END IF
622
623	CALL xtp(imr, jnp, iml, j1, j2, jn, js, pu, dq(1,1,k,ic), wk1(1,1,2), &
624	crx, fx1, xmass, iord)
625
626	CALL ytp(imr, jnp, j1, j2, acosp, rcap, dq(1,1,k,ic), wk1(1,1,1), cry, &
627	dc2, ymass, wk1(1,1,3), wk1(1,1,4), wk1(1,1,5), wk1(1,1,6), jord)
628
629	END DO
630	END DO
631
632	! ***** Compute vertical mass flux (same unit as PS) *********
633
634	! 1st step: compute total column mass CONVERGENCE.
635
636	DO j = 1, jnp
637	DO i = 1, imr
638	cry(i, j) = dpi(i, j, 1)
639	END DO
640	END DO
641
642	DO k = 2, nlay
643	DO j = 1, jnp
644	DO i = 1, imr
645	cry(i, j) = cry(i, j) + dpi(i, j, k)
646	END DO
647	END DO
648	END DO
649
650	DO j = 1, jnp
651	DO i = 1, imr
652
653	! 2nd step: compute PS2 (PS at n+1) using the hydrostatic assumption.
654	! Changes (increases) to surface pressure = total column mass
655	! convergence
656
657	ps2(i, j) = ps1(i, j) + cry(i, j)
658
659	! 3rd step: compute vertical mass flux from mass conservation
660	! principle.
661
662	w(i, j, 1) = dpi(i, j, 1) - dbk(1)*cry(i, j)
663	w(i, j, nlay) = 0.
664	END DO
665	END DO
666
667	DO k = 2, nlay - 1
668	DO j = 1, jnp
669	DO i = 1, imr
670	w(i, j, k) = w(i, j, k-1) + dpi(i, j, k) - dbk(k)*cry(i, j)
671	END DO
672	END DO
673	END DO
674
675	DO k = 1, nlay
676	DO j = 1, jnp
677	DO i = 1, imr
678	delp2(i, j, k) = dap(k) + dbk(k)*ps2(i, j)
679	END DO
680	END DO
681	END DO
682
683	krd = max(3, kord)
684	DO ic = 1, nc
685
686	! **60******0*****0*****0*****0*****0********72
687
688	CALL fzppm(imr, jnp, nlay, j1, dq(1,1,1,ic), w, q(1,1,1,ic), wk1, dpi, &
689	dc2, crx, cry, pu, pv, xmass, ymass, delp1, krd)
690
691
692	IF (fill) CALL qckxyz(dq(1,1,1,ic), dc2, imr, jnp, nlay, j1, j2, cosp, &
693	acosp, .FALSE., ic, nstep)
694
695	! Recover tracer mixing ratio from "density" using predicted
696	! "air density" (pressure thickness) at time-level n+1
697
698	DO k = 1, nlay
699	DO j = 1, jnp
700	DO i = 1, imr
701	q(i, j, k, ic) = dq(i, j, k, ic)/delp2(i, j, k)
702	END DO
703	END DO
704	END DO
705
706	IF (j1/=2) THEN
707	DO k = 1, nlay
708	DO i = 1, imr
709	! j=1 c'est le pôle Sud, j=JNP c'est le pôle Nord
710	q(i, 2, k, ic) = q(i, 1, k, ic)
711	q(i, jmr, k, ic) = q(i, jmp, k, ic)
712	END DO
713	END DO
714	END IF
715	END DO
716
717	IF (j1/=2) THEN
718	DO k = 1, nlay
719	DO i = 1, imr
720	w(i, 2, k) = w(i, 1, k)
721	w(i, jmr, k) = w(i, jnp, k)
722	END DO
723	END DO
724	END IF
725
726	END SUBROUTINE ppm3d
727
728	end module ppm3d_m