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
353	guez	167	IF (igd==0) THEN
354			! Convert winds on A-Grid to Courant number on C-Grid.
355			CALL a2c(u(1,1,k), v(1,1,k), imr, jmr, j1, j2, crx, cry, dtdx5, dtdy5)
356			ELSE
357			! Convert winds on C-grid to Courant number
358			DO j = j1, j2
359			DO i = 2, imr
360			crx(i, j) = dtdx(j)*u(i-1, j, k)
361			END DO
362			END DO
363	guez	81
364
365	guez	167	DO j = j1, j2
366			crx(1, j) = dtdx(j)*u(imr, j, k)
367			END DO
368	guez	81
369	guez	167	DO i = 1, imr*jmr
370			cry(i, 2) = dtdy*v(i, 1, k)
371			END DO
372			END IF
373	guez	81
374	guez	167	! Determine JS and JN
375			js = j1
376			jn = j2
377	guez	81
378	guez	167	DO j = js0, j1 + 1, -1
379			DO i = 1, imr
380			IF (abs(crx(i,j))>1.) THEN
381			js = j
382			GO TO 2222
383			END IF
384			END DO
385			END DO
386	guez	81
387	guez	167	2222 CONTINUE
388			DO j = jn0, j2 - 1
389			DO i = 1, imr
390			IF (abs(crx(i,j))>1.) THEN
391			jn = j
392			GO TO 2233
393			END IF
394			END DO
395			END DO
396			2233 CONTINUE
397	guez	81
398	guez	167	IF (j1/=2) THEN ! Enlarged polar cap.
399			DO i = 1, imr
400			dpi(i, 2, k) = 0.
401			dpi(i, jmr, k) = 0.
402			END DO
403			END IF
404	guez	81
405	guez	167	! ***** Compute horizontal mass fluxes **********
406	guez	81
407	guez	167	! N-S component
408			DO j = j1, j2 + 1
409			d5 = 0.5*cose(j)
410			DO i = 1, imr
411			ymass(i, j) = cry(i, j)d5(delp2(i,j,k)+delp2(i,j-1,k))
412			END DO
413			END DO
414	guez	81
415	guez	167	DO j = j1, j2
416			DO i = 1, imr
417			dpi(i, j, k) = (ymass(i,j)-ymass(i,j+1))*acosp(j)
418			END DO
419			END DO
420	guez	81
421	guez	167	! Poles
422			sum1 = ymass(imr, j1)
423			sum2 = ymass(imr, j2+1)
424			DO i = 1, imr - 1
425			sum1 = sum1 + ymass(i, j1)
426			sum2 = sum2 + ymass(i, j2+1)
427			END DO
428	guez	81
429	guez	167	sum1 = -sum1*rcap
430			sum2 = sum2*rcap
431			DO i = 1, imr
432			dpi(i, 1, k) = sum1
433			dpi(i, jnp, k) = sum2
434			END DO
435	guez	81
436	guez	167	! E-W component
437	guez	81
438	guez	167	DO j = j1, j2
439			DO i = 2, imr
440			pu(i, j) = 0.5*(delp2(i,j,k)+delp2(i-1,j,k))
441			END DO
442			END DO
443	guez	81
444	guez	167	DO j = j1, j2
445			pu(1, j) = 0.5*(delp2(1,j,k)+delp2(imr,j,k))
446			END DO
447	guez	81
448	guez	167	DO j = j1, j2
449			DO i = 1, imr
450			xmass(i, j) = pu(i, j)*crx(i, j)
451			END DO
452			END DO
453	guez	81
454	guez	167	DO j = j1, j2
455			DO i = 1, imr - 1
456			dpi(i, j, k) = dpi(i, j, k) + xmass(i, j) - xmass(i+1, j)
457			END DO
458			END DO
459	guez	81
460	guez	167	DO j = j1, j2
461			dpi(imr, j, k) = dpi(imr, j, k) + xmass(imr, j) - xmass(1, j)
462			END DO
463	guez	81
464	guez	167	DO j = j1, j2
465			DO i = 1, imr - 1
466			ua(i, j) = 0.5*(crx(i,j)+crx(i+1,j))
467			END DO
468			END DO
469	guez	81
470	guez	167	DO j = j1, j2
471			ua(imr, j) = 0.5*(crx(imr,j)+crx(1,j))
472			END DO
473			! cccccccccccccccccccccccccccccccccccccccccccccccccccccc
474			! Rajouts pour LMDZ.3.3
475			! cccccccccccccccccccccccccccccccccccccccccccccccccccccc
476			DO i = 1, imr
477			DO j = 1, jnp
478			va(i, j) = 0.
479			END DO
480			END DO
481	guez	81
482	guez	167	DO i = 1, imr*(jmr-1)
483			va(i, 2) = 0.5*(cry(i,2)+cry(i,3))
484			END DO
485	guez	81
486	guez	167	IF (j1==2) THEN
487			imh = imr/2
488			DO i = 1, imh
489			va(i, 1) = 0.5*(cry(i,2)-cry(i+imh,2))
490			va(i+imh, 1) = -va(i, 1)
491			va(i, jnp) = 0.5*(cry(i,jnp)-cry(i+imh,jmr))
492			va(i+imh, jnp) = -va(i, jnp)
493			END DO
494			va(imr, 1) = va(1, 1)
495			va(imr, jnp) = va(1, jnp)
496			END IF
497	guez	81
498	guez	167	! **60******0*****0*****0*****0*****0********72
499			DO ic = 1, nc
500	guez	81
501	guez	167	DO i = 1, imjm
502			wk1(i, 1, 1) = 0.
503			wk1(i, 1, 2) = 0.
504			END DO
505	guez	81
506	guez	167	! E-W advective cross term
507			DO j = j1, j2
508			IF (j>js .AND. j<jn) cycle
509	guez	81
510	guez	167	DO i = 1, imr
511			qtmp(i) = q(i, j, k, ic)
512			END DO
513	guez	81
514	guez	167	DO i = -iml, 0
515			qtmp(i) = q(imr+i, j, k, ic)
516			qtmp(imr+1-i) = q(1-i, j, k, ic)
517			END DO
518	guez	81
519	guez	167	DO i = 1, imr
520			iu = ua(i, j)
521			ru = ua(i, j) - iu
522			iiu = i - iu
523			IF (ua(i,j)>=0.) THEN
524			wk1(i, j, 1) = qtmp(iiu) + ru*(qtmp(iiu-1)-qtmp(iiu))
525			ELSE
526			wk1(i, j, 1) = qtmp(iiu) + ru*(qtmp(iiu)-qtmp(iiu+1))
527			END IF
528			wk1(i, j, 1) = wk1(i, j, 1) - qtmp(i)
529			END DO
530			END DO
531	guez	81
532	guez	167	IF (jn/=0) THEN
533			DO j = js + 1, jn - 1
534	guez	81
535	guez	167	DO i = 1, imr
536			qtmp(i) = q(i, j, k, ic)
537			END DO
538	guez	81
539	guez	167	qtmp(0) = q(imr, j, k, ic)
540			qtmp(imr+1) = q(1, j, k, ic)
541
542			DO i = 1, imr
543			iu = i - ua(i, j)
544			wk1(i, j, 1) = ua(i, j)*(qtmp(iu)-qtmp(iu+1))
545			END DO
546			END DO
547			END IF
548			! **60******0*****0*****0*****0*****0********72
549			! Contribution from the N-S advection
550			DO i = 1, imr*(j2-j1+1)
551			jt = float(j1) - va(i, j1)
552			wk1(i, j1, 2) = va(i, j1)*(q(i,jt,k,ic)-q(i,jt+1,k,ic))
553	guez	81	END DO
554
555	guez	167	DO i = 1, imjm
556			wk1(i, 1, 1) = q(i, 1, k, ic) + 0.5*wk1(i, 1, 1)
557			wk1(i, 1, 2) = q(i, 1, k, ic) + 0.5*wk1(i, 1, 2)
558	guez	81	END DO
559
560	guez	167	IF (cross) THEN
561			! Add cross terms in the vertical direction.
562			IF (iord>=2) THEN
563			iad = 2
564			ELSE
565			iad = 1
566			END IF
567	guez	81
568	guez	167	IF (jord>=2) THEN
569			jad = 2
570			ELSE
571			jad = 1
572			END IF
573			CALL xadv(imr, jnp, j1, j2, wk1(1,1,2), ua, js, jn, iml, dc2, iad)
574			CALL yadv(imr, jnp, j1, j2, wk1(1,1,1), va, pv, w, jad)
575			DO j = 1, jnp
576			DO i = 1, imr
577			q(i, j, k, ic) = q(i, j, k, ic) + dc2(i, j) + pv(i, j)
578			END DO
579			END DO
580			END IF
581	guez	81
582	guez	167	CALL xtp(imr, jnp, iml, j1, j2, jn, js, pu, dq(1,1,k,ic), wk1(1,1,2), &
583			crx, fx1, xmass, iord)
584	guez	81
585	guez	167	CALL ytp(imr, jnp, j1, j2, acosp, rcap, dq(1,1,k,ic), wk1(1,1,1), cry, &
586			dc2, ymass, wk1(1,1,3), wk1(1,1,4), wk1(1,1,5), wk1(1,1,6), jord)
587	guez	81
588	guez	167	END DO
589	guez	81	END DO
590
591	guez	167	! ***** Compute vertical mass flux (same unit as PS) *********
592	guez	81
593	guez	167	! 1st step: compute total column mass CONVERGENCE.
594	guez	81
595	guez	167	DO j = 1, jnp
596			DO i = 1, imr
597			cry(i, j) = dpi(i, j, 1)
598			END DO
599	guez	81	END DO
600
601	guez	167	DO k = 2, nlay
602			DO j = 1, jnp
603			DO i = 1, imr
604			cry(i, j) = cry(i, j) + dpi(i, j, k)
605			END DO
606			END DO
607	guez	81	END DO
608
609	guez	167	DO j = 1, jnp
610			DO i = 1, imr
611	guez	81
612	guez	167	! 2nd step: compute PS2 (PS at n+1) using the hydrostatic assumption.
613			! Changes (increases) to surface pressure = total column mass
614			! convergence
615	guez	81
616	guez	167	ps2(i, j) = ps1(i, j) + cry(i, j)
617	guez	81
618	guez	167	! 3rd step: compute vertical mass flux from mass conservation
619			! principle.
620	guez	81
621	guez	167	w(i, j, 1) = dpi(i, j, 1) - dbk(1)*cry(i, j)
622			w(i, j, nlay) = 0.
623			END DO
624	guez	81	END DO
625
626	guez	167	DO k = 2, nlay - 1
627			DO j = 1, jnp
628			DO i = 1, imr
629			w(i, j, k) = w(i, j, k-1) + dpi(i, j, k) - dbk(k)*cry(i, j)
630			END DO
631			END DO
632	guez	81	END DO
633
634	guez	167	DO k = 1, nlay
635			DO j = 1, jnp
636			DO i = 1, imr
637			delp2(i, j, k) = dap(k) + dbk(k)*ps2(i, j)
638			END DO
639			END DO
640	guez	81	END DO
641
642	guez	167	krd = max(3, kord)
643			DO ic = 1, nc
644	guez	81
645	guez	167	! **60******0*****0*****0*****0*****0********72
646	guez	81
647	guez	167	CALL fzppm(imr, jnp, nlay, j1, dq(1,1,1,ic), w, q(1,1,1,ic), wk1, dpi, &
648			dc2, crx, cry, pu, pv, xmass, ymass, delp1, krd)
649	guez	81
650
651	guez	167	IF (fill) CALL qckxyz(dq(1,1,1,ic), dc2, imr, jnp, nlay, j1, j2, cosp, &
652			acosp, .FALSE., ic, nstep)
653	guez	81
654	guez	167	! Recover tracer mixing ratio from "density" using predicted
655			! "air density" (pressure thickness) at time-level n+1
656	guez	81
657	guez	167	DO k = 1, nlay
658			DO j = 1, jnp
659			DO i = 1, imr
660			q(i, j, k, ic) = dq(i, j, k, ic)/delp2(i, j, k)
661			END DO
662			END DO
663			END DO
664
665			IF (j1/=2) THEN
666			DO k = 1, nlay
667			DO i = 1, imr
668	guez	168	! j=1 c'est le p\^ole Sud, j=JNP c'est le p\^ole Nord
669	guez	167	q(i, 2, k, ic) = q(i, 1, k, ic)
670			q(i, jmr, k, ic) = q(i, jmp, k, ic)
671			END DO
672			END DO
673			END IF
674	guez	81	END DO
675
676			IF (j1/=2) THEN
677	guez	167	DO k = 1, nlay
678			DO i = 1, imr
679			w(i, 2, k) = w(i, 1, k)
680			w(i, jmr, k) = w(i, jnp, k)
681			END DO
682			END DO
683	guez	81	END IF
684
685	guez	167	END SUBROUTINE ppm3d
686	guez	81
687	guez	167	end module ppm3d_m