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	module ppm3d_m
2
3	implicit none
4
5	contains
6
7	SUBROUTINE ppm3d(igd, q, ps1, ps2, u, v, w, ndt, iord, jord, kord, nc, imr, &
8	jnp, j1, nlay, ap, bp, pt, ae, fill, umax)
9
10	! INPUT:
11	! =============
12
13	! 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
22	! 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
28	! The pressure at layer edges are defined as follows:
29
30	! p(i,j,k) = AP(k)PT + BP(k)PS(i,j) (1)
31
32	! 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
38	! 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
41	! 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
46	! Cf. ppm3d.txt.
47
48	! 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
51	! 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
55	! IGD = 1 GEOS-GCM C-Grid
56	! Cf. ppm3d.txt.
57
58	! 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
62	! V must be defined at j=1 and j=JMR if IGD=1
63	! V at JNP need not be given.
64
65	! 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
70	! 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
76	! 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
83	! _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
98	! *PPM: Piece-wise Parabolic Method
99
100	! 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
104	! 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
111	! AE: Radius of the sphere (meters).
112	! Recommended value for the planet earth: 6.371E6
113
114	! fill(logical): flag to do filling for negatives (see note below).
115
116	! 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
119	! =============
120	! Output
121	! =============
122
123	! 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
135	! ********************************************************************
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
143	! A larger polar cap is used if j1=3 (recommended for C-Grid winds or when
144	! winds are noisy near poles).
145
146	! Flux-Form Semi-Lagrangian transport in the East-West direction is used
147	! when and where Courant number is greater than one.
148
149	! 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
154	! PPM is 4th order accurate when grid spacing is uniform (x & y); 3rd
155	! order accurate for non-uniform grid (vertical sigma coord.).
156
157	! Time step is limitted only by transport in the meridional direction.
158	! (the FFSL scheme is not implemented in the meridional direction).
159
160	! 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
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
172	! **60******0*****0*****0*****0*****0********72
173
174	! User modifiable parameters
175
176	integer, PARAMETER:: jmax=361, kmax=150
177
178	! **60******0*****0*****0*****0*****0********72
179
180	! Input-Output arrays
181
182	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	w(imr, jnp, nlay), ndt, umax
187	INTEGER imrd2
188	REAL pt
189	LOGICAL cross, fill
190
191	! Local dynamic arrays
192
193	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
199	! Local static arrays
200
201	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
211	jmr = jnp - 1
212	imjm = imr*jnp
213	j2 = jnp - j1 + 1
214	nstep = nstep + 1
215
216	! ********* Initialization ********************
217	IF (nstep==1) THEN
218
219	WRITE (6, *) '------------------------------------ '
220	WRITE (6, *) 'NASA/GSFC Transport Core Version 4.5'
221	WRITE (6, *) '------------------------------------ '
222
223	WRITE (6, *) 'IMR=', imr, ' JNP=', jnp, ' NLAY=', nlay, ' j1=', j1
224	WRITE (6, *) 'NC=', nc, iord, jord, kord, ndt
225
226	! 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
241
242	IF (jmax<jnp .OR. kmax<nlay) THEN
243	WRITE (6, *) 'Jmax or Kmax is too small'
244	STOP
245	END IF
246
247	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
252	pi = 4.*atan(1.)
253	dl = 2.*pi/float(imr)
254	dp = pi/float(jmr)
255
256	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
264	DO j = 2, jmr
265	acosp(j) = 1./cosp(j)
266	END DO
267
268	! Inverse of the Scaled polar cap area.
269
270	rcap = dp/(imr(1.-cos((j1-1.5)dp)))
271	acosp(1) = rcap
272	acosp(jnp) = rcap
273	END IF
274
275	IF (ndt0/=ndt) THEN
276	dt = ndt
277	ndt0 = ndt
278
279	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
289	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
297	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
303
304	DO j = 2, jmr
305	dtdx(j) = dt/(dlaecosp(j))
306
307	dtdx5(j) = 0.5*dtdx(j)
308	END DO
309
310
311	dtdy = dt/(ae*dp)
312	dtdy5 = 0.5*dtdy
313
314	END IF
315
316	! ********* End Initialization ********************
317
318	! 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	END DO
327
328
329	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	DO ic = 1, nc
342	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	END DO
350
351	DO k = 1, nlay
352
353	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
364
365	DO j = j1, j2
366	crx(1, j) = dtdx(j)*u(imr, j, k)
367	END DO
368
369	DO i = 1, imr*jmr
370	cry(i, 2) = dtdy*v(i, 1, k)
371	END DO
372	END IF
373
374	! Determine JS and JN
375	js = j1
376	jn = j2
377
378	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
387	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
398	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
405	! ***** Compute horizontal mass fluxes **********
406
407	! 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
415	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
421	! 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
429	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
436	! E-W component
437
438	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
444	DO j = j1, j2
445	pu(1, j) = 0.5*(delp2(1,j,k)+delp2(imr,j,k))
446	END DO
447
448	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
454	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
460	DO j = j1, j2
461	dpi(imr, j, k) = dpi(imr, j, k) + xmass(imr, j) - xmass(1, j)
462	END DO
463
464	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
470	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
482	DO i = 1, imr*(jmr-1)
483	va(i, 2) = 0.5*(cry(i,2)+cry(i,3))
484	END DO
485
486	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
498	! **60******0*****0*****0*****0*****0********72
499	DO ic = 1, nc
500
501	DO i = 1, imjm
502	wk1(i, 1, 1) = 0.
503	wk1(i, 1, 2) = 0.
504	END DO
505
506	! E-W advective cross term
507	DO j = j1, j2
508	IF (j>js .AND. j<jn) cycle
509
510	DO i = 1, imr
511	qtmp(i) = q(i, j, k, ic)
512	END DO
513
514	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
519	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
532	IF (jn/=0) THEN
533	DO j = js + 1, jn - 1
534
535	DO i = 1, imr
536	qtmp(i) = q(i, j, k, ic)
537	END DO
538
539	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	END DO
554
555	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	END DO
559
560	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
568	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
582	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
585	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
588	END DO
589	END DO
590
591	! ***** Compute vertical mass flux (same unit as PS) *********
592
593	! 1st step: compute total column mass CONVERGENCE.
594
595	DO j = 1, jnp
596	DO i = 1, imr
597	cry(i, j) = dpi(i, j, 1)
598	END DO
599	END DO
600
601	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	END DO
608
609	DO j = 1, jnp
610	DO i = 1, imr
611
612	! 2nd step: compute PS2 (PS at n+1) using the hydrostatic assumption.
613	! Changes (increases) to surface pressure = total column mass
614	! convergence
615
616	ps2(i, j) = ps1(i, j) + cry(i, j)
617
618	! 3rd step: compute vertical mass flux from mass conservation
619	! principle.
620
621	w(i, j, 1) = dpi(i, j, 1) - dbk(1)*cry(i, j)
622	w(i, j, nlay) = 0.
623	END DO
624	END DO
625
626	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	END DO
633
634	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	END DO
641
642	krd = max(3, kord)
643	DO ic = 1, nc
644
645	! **60******0*****0*****0*****0*****0********72
646
647	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
650
651	IF (fill) CALL qckxyz(dq(1,1,1,ic), dc2, imr, jnp, nlay, j1, j2, cosp, &
652	acosp, .FALSE., ic, nstep)
653
654	! Recover tracer mixing ratio from "density" using predicted
655	! "air density" (pressure thickness) at time-level n+1
656
657	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	! j=1 c'est le p\^ole Sud, j=JNP c'est le p\^ole Nord
669	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	END DO
675
676	IF (j1/=2) THEN
677	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	END IF
684
685	END SUBROUTINE ppm3d
686
687	end module ppm3d_m