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)

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

    ! Cf. ppm3d.txt.

    ! 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
    ! Cf. ppm3d.txt.

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