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