dyn3d/PPM3d/ppm3d.f

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

  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./
  SAVE dtdy, dtdy5, rcap, js0, jn0, iml, dtdx, dtdx5, acosp, cosp, cose, dap, &
    dbk


  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) GO TO 250

        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
250   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

  RETURN
END SUBROUTINE ppm3d

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