source: LMDZ6/trunk/libf/dyn3d_common/ppm3d.f90 @ 5248

!
! $Id: ppm3d.f90 5246 2024-10-21 12:58:45Z abarral $
!

!From lin@explorer.gsfc.nasa.gov Wed Apr 15 17:44:44 1998
!Date: Wed, 15 Apr 1998 11:37:03 -0400
!From: lin@explorer.gsfc.nasa.gov
!To: Frederic.Hourdin@lmd.jussieu.fr
!Subject: 3D transport module of the GSFC CTM and GEOS GCM


!This code is sent to you by S-J Lin, DAO, NASA-GSFC

!Note: this version is intended for machines like CRAY
!-90. No multitasking directives implemented.


! ********************************************************************
!
! TransPort Core for Goddard Chemistry Transport Model (G-CTM), Goddard
! Earth Observing System General Circulation Model (GEOS-GCM), and Data
! Assimilation System (GEOS-DAS).
!
! ********************************************************************
!
! Purpose: given horizontal winds on  a hybrid sigma-p surfaces,
!      one call to tpcore updates the 3-D mixing ratio
!      fields one time step (NDT). [vertical mass flux is computed
!      internally consistent with the discretized hydrostatic mass
!      continuity equation of the C-Grid GEOS-GCM (for IGD=1)].
!
! Schemes: Multi-dimensional Flux Form Semi-Lagrangian (FFSL) scheme based
!      on the van Leer or PPM.
!      (see Lin and Rood 1996).
! Version 4.5
! Last modified: Dec. 5, 1996
! Major changes from version 4.0: a more general vertical hybrid sigma-
! pressure coordinate.
! Subroutines modified: xtp, ytp, fzppm, qckxyz
! Subroutines deleted: vanz
!
! Author: Shian-Jiann Lin
! mail address:
!             Shian-Jiann Lin*
!             Code 910.3, NASA/GSFC, Greenbelt, MD 20771
!             Phone: 301-286-9540
!             E-mail: lin@dao.gsfc.nasa.gov
!
! *affiliation:
!             Joint Center for Earth Systems Technology
!             The University of Maryland Baltimore County
!             NASA - Goddard Space Flight Center
! References:
!
! 1. Lin, S.-J., and R. B. Rood, 1996: Multidimensional flux form semi-
!    Lagrangian transport schemes. Mon. Wea. Rev., 124, 2046-2070.
!
! 2. Lin, S.-J., W. C. Chao, Y. C. Sud, and G. K. Walker, 1994: A class of
!    the van Leer-type transport schemes and its applications to the moist-
!    ure transport in a General Circulation Model. Mon. Wea. Rev., 122,
!    1575-1593.
!
! ****6***0*********0*********0*********0*********0*********0**********72
!
subroutine ppm3d(IGD,Q,PS1,PS2,U,V,W,NDT,IORD,JORD,KORD,NC,IMR, &
        JNP,j1,NLAY,AP,BP,PT,AE,fill,dum,Umax)

  implicit none

  ! rajout de d�clarations
  !  integer Jmax,kmax,ndt0,nstep,k,j,i,ic,l,js,jn,imh,iad,jad,krd
  !  integer iu,iiu,j2,jmr,js0,jt
  !  real dtdy,dtdy5,rcap,iml,jn0,imjm,pi,dl,dp
  !  real dt,cr1,maxdt,ztc,d5,sum1,sum2,ru
  !
  ! ********************************************************************
  !
  ! =============
  ! INPUT:
  ! =============
  !
  ! Q(IMR,JNP,NLAY,NC): mixing ratios at current time (t)
  ! NC: total # of constituents
  ! IMR: first dimension (E-W); # of Grid intervals in E-W is IMR
  ! JNP: 2nd dimension (N-S); # of Grid intervals in N-S is JNP-1
  ! NLAY: 3rd dimension (# of layers); vertical index increases from 1 at
  !   the model top to NLAY near the surface (see fig. below).
  !   It is assumed that 6 <= NLAY <= JNP (for dynamic memory allocation)
  !
  ! PS1(IMR,JNP): surface pressure at current time (t)
  ! PS2(IMR,JNP): surface pressure at mid-time-level (t+NDT/2)
  ! PS2 is replaced by the predicted PS (at t+NDT) on output.
  ! Note: surface pressure can have any unit or can be multiplied by any
  !   const.
  !
  ! The pressure at layer edges are defined as follows:
  !
  !    p(i,j,k) = AP(k)*PT  +  BP(k)*PS(i,j)          (1)
  !
  ! Where PT is a constant having the same unit as PS.
  ! AP and BP are unitless constants given at layer edges
  ! defining the vertical coordinate.
  ! BP(1) = 0., BP(NLAY+1) = 1.
  ! The pressure at the model top is PTOP = AP(1)*PT
  !
  ! For pure sigma system set AP(k) = 1 for all k, PT = PTOP,
  ! BP(k) = sige(k) (sigma at edges), PS = Psfc - PTOP.
  !
  ! Note: the sigma-P coordinate is a subset of Eq. 1, which in turn
  ! is a subset of the following even more general sigma-P-thelta coord.
  ! currently under development.
  !  p(i,j,k) = (AP(k)*PT + BP(k)*PS(i,j))/(D(k)-C(k)*TE**(-1/kapa))
  !
  !              /////////////////////////////////
  !          / \ ------------- PTOP --------------  AP(1), BP(1)
  !           |
  !    delp(1)    |  ........... Q(i,j,1) ............
  !           |
  !  W(1)    \ / ---------------------------------  AP(2), BP(2)
  !
  !
  !
  ! W(k-1)   / \ ---------------------------------  AP(k), BP(k)
  !           |
  !    delp(K)    |  ........... Q(i,j,k) ............
  !           |
  !  W(k)    \ / ---------------------------------  AP(k+1), BP(k+1)
  !
  !
  !
  !          / \ ---------------------------------  AP(NLAY), BP(NLAY)
  !           |
  !  delp(NLAY)   |  ........... Q(i,j,NLAY) .........
  !           |
  !   W(NLAY)=0  \ / ------------- surface ----------- AP(NLAY+1), BP(NLAY+1)
  !             //////////////////////////////////
  !
  ! U(IMR,JNP,NLAY) & V(IMR,JNP,NLAY):winds (m/s) at mid-time-level (t+NDT/2)
  ! U and V may need to be polar filtered in advance in some cases.
  !
  ! IGD:      grid type on which winds are defined.
  ! IGD = 0:  A-Grid  [all variables defined at the same point from south
  !               pole (j=1) to north pole (j=JNP) ]
  !
  ! IGD = 1  GEOS-GCM C-Grid
  !                                  [North]
  !
  !                                   V(i,j)
  !                                      |
  !                                      |
  !                                      |
  !                         U(i-1,j)---Q(i,j)---U(i,j) [EAST]
  !                                      |
  !                                      |
  !                                      |
  !                                   V(i,j-1)
  !
  !     U(i,  1) is defined at South Pole.
  !     V(i,  1) is half grid north of the South Pole.
  !     V(i,JMR) is half grid south of the North Pole.
  !
  !     V must be defined at j=1 and j=JMR if IGD=1
  !     V at JNP need not be given.
  !
  ! NDT: time step in seconds (need not be constant during the course of
  !  the integration). Suggested value: 30 min. for 4x5, 15 min. for 2x2.5
  !  (Lat-Lon) resolution. Smaller values are recommanded if the model
  !  has a well-resolved stratosphere.
  !
  ! J1 defines the size of the polar cap:
  ! South polar cap edge is located at -90 + (j1-1.5)*180/(JNP-1) deg.
  ! North polar cap edge is located at  90 - (j1-1.5)*180/(JNP-1) deg.
  ! There are currently only two choices (j1=2 or 3).
  ! IMR must be an even integer if j1 = 2. Recommended value: J1=3.
  !
  ! IORD, JORD, and KORD are integers controlling various options in E-W, N-S,
  ! and vertical transport, respectively. Recommended values for positive
  ! definite scalars: IORD=JORD=3, KORD=5. Use KORD=3 for non-
  ! positive definite scalars or when linear correlation between constituents
  ! is to be maintained.
  !
  !  _ORD=
  !    1: 1st order upstream scheme (too diffusive, not a useful option; it
  !       can be used for debugging purposes; this is THE only known "linear"
  !       monotonic advection scheme.).
  !    2: 2nd order van Leer (full monotonicity constraint;
  !       see Lin et al 1994, MWR)
  !    3: monotonic PPM* (slightly improved PPM of Collela & Woodward 1984)
  !    4: semi-monotonic PPM (same as 3, but overshoots are allowed)
  !    5: positive-definite PPM (constraint on the subgrid distribution is
  !       only strong enough to prevent generation of negative values;
  !       both overshoots & undershoots are possible).
  !    6: un-constrained PPM (nearly diffusion free; slightly faster but
  !       positivity not quaranteed. Use this option only when the fields
  !       and winds are very smooth).
  !
  ! *PPM: Piece-wise Parabolic Method
  !
  ! Note that KORD <=2 options are no longer supported. DO not use option 4 or 5.
  ! for non-positive definite scalars (such as Ertel Potential Vorticity).
  !
  ! The implicit numerical diffusion decreases as _ORD increases.
  ! The last two options (ORDER=5, 6) should only be used when there is
  ! significant explicit diffusion (such as a turbulence parameterization). You
  ! might get dispersive results otherwise.
  ! No filter of any kind is applied to the constituent fields here.
  !
  ! AE: Radius of the sphere (meters).
  ! Recommended value for the planet earth: 6.371E6
  !
  ! fill(logical):   flag to do filling for negatives (see note below).
  !
  ! Umax: Estimate (upper limit) of the maximum U-wind speed (m/s).
  ! (220 m/s is a good value for troposphere model; 280 m/s otherwise)
  !
  ! =============
  ! Output
  ! =============
  !
  ! Q: mixing ratios at future time (t+NDT) (original values are over-written)
  ! W(NLAY): large-scale vertical mass flux as diagnosed from the hydrostatic
  !      relationship. W will have the same unit as PS1 and PS2 (eg, mb).
  !      W must be divided by NDT to get the correct mass-flux unit.
  !      The vertical Courant number C = W/delp_UPWIND, where delp_UPWIND
  !      is the pressure thickness in the "upwind" direction. For example,
  !      C(k) = W(k)/delp(k)   if W(k) > 0;
  !      C(k) = W(k)/delp(k+1) if W(k) < 0.
  !          ( W > 0 is downward, ie, toward surface)
  ! PS2: predicted PS at t+NDT (original values are over-written)
  !
  ! ********************************************************************
  ! NOTES:
  ! This forward-in-time upstream-biased transport scheme reduces to
  ! the 2nd order center-in-time center-in-space mass continuity eqn.
  ! if Q = 1 (constant fields will remain constant). This also ensures
  ! that the computed vertical velocity to be identical to GEOS-1 GCM
  ! for on-line transport.
  !
  ! A larger polar cap is used if j1=3 (recommended for C-Grid winds or when
  ! winds are noisy near poles).
  !
  ! Flux-Form Semi-Lagrangian transport in the East-West direction is used
  ! when and where Courant # 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
  !

  real :: Q(IMR,JNP,NLAY,NC),PS1(IMR,JNP),PS2(IMR,JNP), &
        U(IMR,JNP,NLAY),V(IMR,JNP,NLAY),AP(NLAY+1), &
        BP(NLAY+1),W(IMR,JNP,NLAY),NDT,val(NLAY),Umax
  integer :: IGD,IORD,JORD,KORD,NC,IMR,JNP,j1,NLAY,AE
  integer :: IMRD2
  real :: PT
  logical :: cross, fill, dum
  !
  ! Local dynamic arrays
  !
  real :: CRX(IMR,JNP),CRY(IMR,JNP),xmass(IMR,JNP),ymass(IMR,JNP), &
        fx1(IMR+1),DPI(IMR,JNP,NLAY),delp1(IMR,JNP,NLAY), &
        WK1(IMR,JNP,NLAY),PU(IMR,JNP),PV(IMR,JNP),DC2(IMR,JNP), &
        delp2(IMR,JNP,NLAY),DQ(IMR,JNP,NLAY,NC),VA(IMR,JNP), &
        UA(IMR,JNP),qtmp(-IMR:2*IMR)
  !
  ! Local static  arrays
  !
  real :: DTDX(Jmax), DTDX5(Jmax), acosp(Jmax), &
        cosp(Jmax), cose(Jmax), DAP(kmax),DBK(Kmax)
  data NDT0, NSTEP /0, 0/
  data cross /.true./
  REAL :: DTDY, DTDY5, RCAP
  INTEGER :: JS0, JN0, IML, JMR, IMJM
  SAVE DTDY, DTDY5, RCAP, JS0, JN0, IML, &
        DTDX, DTDX5, ACOSP, COSP, COSE, DAP,DBK
  !
  INTEGER :: NDT0, NSTEP, j2, k,j,i,ic,l,JS,JN,IMH
  INTEGER :: IU,IIU,JT,iad,jad,krd
  REAL :: r23,r3,PI,DL,DP,DT,CR1,MAXDT,ZTC,D5
  REAL :: sum1,sum2,ru

  JMR = JNP -1
  IMJM  = IMR*JNP
  j2 = JNP - j1 + 1
  NSTEP = NSTEP + 1
  !
  ! *********** Initialization **********************
  if(NSTEP.eq.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.LT.6) then
    write(6,*) 'NLAY must be >= 6'
    stop
  endif
  if (JNP.LT.NLAY) then
     write(6,*) 'JNP must be >= NLAY'
    stop
  endif
  IMRD2=mod(IMR,2)
  if (j1.eq.2.and.IMRD2.NE.0) then
     write(6,*) 'if j1=2 IMR must be an even integer'
    stop
  endif

  !
  if(Jmax.lt.JNP .or. Kmax.lt.NLAY) then
    write(6,*) 'Jmax or Kmax is too small'
    stop
  endif
  !
  DO k=1,NLAY
  DAP(k) = (AP(k+1) - AP(k))*PT
  DBK(k) =  BP(k+1) - BP(k)
  ENDDO
  !
  PI = 4. * ATAN(1.)
  DL = 2.*PI / REAL(IMR)
  DP =    PI / REAL(JMR)
  !
  if(IGD.eq.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)
  endif
  !
  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
  endif
  !
  if(NDT0 .ne. NDT) then
  DT   = NDT
  NDT0 = NDT

    if(Umax .lt. 180.) then
     write(6,*) 'Umax may be too small!'
    endif
  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 .lt. abs(NDT)) then
        write(6,*) 'Warning!!! NDT maybe too large!'
  endif
  !
  if(CR1.ge.0.95) then
  JS0 = 0
  JN0 = 0
  IML = IMR-2
  ZTC = 0.
  else
  ZTC  = acos(CR1) * (180./PI)
  !
  JS0 = REAL(JMR)*(90.-ZTC)/180. + 2
  JS0 = max(JS0, J1+1)
  IML = min(6*JS0/(J1-1)+2, 4*IMR/5)
  JN0 = JNP-JS0+1
  endif
  !
  !
  do J=2,JMR
  DTDX(j)  = DT / ( DL*AE*COSP(J) )

  ! print*,'dtdx=',dtdx(j)
  DTDX5(j) = 0.5*DTDX(j)
  enddo
  !

  DTDY  = DT /(AE*DP)
   ! print*,'dtdy=',dtdy
  DTDY5 = 0.5*DTDY
  !
  !  write(6,*) 'J1=',J1,' J2=', J2
  endif
  !
  ! *********** 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)
        enddo
     enddo
  enddo

  !
  if(j1.ne.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
  endif
  !
  ! 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.eq.0) then
  ! Convert winds on A-Grid to Courant # 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 #
  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)
50  continue
  end do
  !
  do i=1,IMR*JMR
    CRY(i,2) = DTDY*V(i,1,k)
  end do
  endif
  !
  ! Determine JS and JN
  JS = j1
  JN = j2
  !
  do j=JS0,j1+1,-1
  do i=1,IMR
  if(abs(CRX(i,j)).GT.1.) then
        JS = j
        go to 2222
  endif
  enddo
  enddo
  !
2222   continue
  do j=JN0,j2-1
  do i=1,IMR
  if(abs(CRX(i,j)).GT.1.) then
        JN = j
        go to 2233
  endif
  enddo
  enddo
2233   continue
  !
  if(j1.ne.2) then           ! Enlarged polar cap.
  do i=1,IMR
  DPI(i,  2,k) = 0.
  DPI(i,JMR,k) = 0.
  enddo
  endif
  !
  ! ******* 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))
  enddo
  enddo
  !
  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)
  enddo
  !
  sum1 = - sum1 * RCAP
  sum2 =   sum2 * RCAP
  do i=1,IMR
  DPI(i,  1,k) = sum1
  DPI(i,JNP,k) = sum2
  enddo
  !
  ! 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))
  enddo
  enddo
  !
  do j=j1,j2
  PU(1,j) = 0.5 * (delp2(1,j,k) + delp2(IMR,j,k))
  enddo
  !
  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))
  enddo
  enddo
  !
  DO j=j1,j2
  UA(imr,j) = 0.5 * (CRX(imr,j)+CRX(1,j))
  enddo
  !cccccccccccccccccccccccccccccccccccccccccccccccccccccc
  ! Rajouts pour LMDZ.3.3
  !cccccccccccccccccccccccccccccccccccccccccccccccccccccc
  do i=1,IMR
     do j=1,JNP
         VA(i,j)=0.
     enddo
  enddo

  do i=1,imr*(JMR-1)
  VA(i,2) = 0.5*(CRY(i,2)+CRY(i,3))
  enddo
  !
  if(j1.eq.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)
  enddo
  VA(IMR,1)=VA(1,1)
  VA(IMR,JNP)=VA(1,JNP)
  endif
  !
  ! ****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.
  enddo
  !
  ! E-W advective cross term
  do j=J1,J2
  if(J.GT.JS  .and. J.LT.JN) GO TO 250
  !
  do i=1,IMR
  qtmp(i) = q(i,j,k,IC)
  enddo
  !
  do i=-IML,0
  qtmp(i)       = q(IMR+i,j,k,IC)
  qtmp(IMR+1-i) = q(1-i,j,k,IC)
  enddo
  !
  DO i=1,IMR
  iu = UA(i,j)
  ru = UA(i,j) - iu
  iiu = i-iu
  if(UA(i,j).GE.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))
  endif
  wk1(i,j,1) = wk1(i,j,1) - qtmp(i)
  END DO
250 continue
  end do
  !
  if(JN.ne.0) then
  do j=JS+1,JN-1
  !
  do i=1,IMR
  qtmp(i) = q(i,j,k,IC)
  enddo
  !
  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))
  enddo
  enddo
  endif
  ! ****6***0*********0*********0*********0*********0*********0**********72
  ! Contribution from the N-S advection
  do i=1,imr*(j2-j1+1)
  JT = REAL(J1) - VA(i,j1)
  wk1(i,j1,2) = VA(i,j1) * (q(i,jt,k,IC) - q(i,jt+1,k,IC))
  enddo
  !
  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)
  enddo
  !
    if(cross) then
  ! Add cross terms in the vertical direction.
    if(IORD .GE. 2) then
            iad = 2
    else
            iad = 1
    endif
  !
    if(JORD .GE. 2) then
            jad = 2
    else
            jad = 1
    endif
  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)
  enddo
  enddo
  endif
  !
  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)
         ! print*,'i=',i,'j=',j,'k=',k,'Q(i,j,k,IC)=',Q(i,j,k,IC)
  enddo
  enddo
  enddo
  !
  if(j1.ne.2) then
  DO k=1,NLAY
  DO I=1,IMR
  ! j=1 c'est le p�le Sud, j=JNP c'est le p�le Nord
  Q(I,  2,k,IC) = Q(I,  1,k,IC)
  Q(I,JMR,k,IC) = Q(I,JNP,k,IC)
  END DO
  END DO
  endif
  end do
  !
  if(j1.ne.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
  endif
  !
  RETURN
END SUBROUTINE ppm3d
!
!****6***0*********0*********0*********0*********0*********0**********72
subroutine FZPPM(IMR,JNP,NLAY,j1,DQ,WZ,P,DC,DQDT,AR,AL,A6, &
        flux,wk1,wk2,wz2,delp,KORD)
  implicit none
  integer,parameter :: kmax = 150
  real,parameter :: R23 = 2./3., R3 = 1./3.
  integer :: IMR,JNP,NLAY,J1,KORD
  real :: WZ(IMR,JNP,NLAY),P(IMR,JNP,NLAY),DC(IMR,JNP,NLAY), &
        wk1(IMR,*),delp(IMR,JNP,NLAY),DQ(IMR,JNP,NLAY), &
        DQDT(IMR,JNP,NLAY)
  ! Assuming JNP >= NLAY
  real :: AR(IMR,*),AL(IMR,*),A6(IMR,*),flux(IMR,*),wk2(IMR,*), &
        wz2(IMR,*)
  integer :: JMR,IMJM,NLAYM1,LMT,K,I,J
  real :: c0,c1,c2,tmp,qmax,qmin,a,b,fct,a1,a2,cm,cp
  !
  JMR = JNP - 1
  IMJM = IMR*JNP
  NLAYM1 = NLAY - 1
  !
  LMT = KORD - 3
  !
  ! ****6***0*********0*********0*********0*********0*********0**********72
  ! Compute DC for PPM
  ! ****6***0*********0*********0*********0*********0*********0**********72
  !
  do k=1,NLAYM1
  do i=1,IMJM
  DQDT(i,1,k) = P(i,1,k+1) - P(i,1,k)
  end do
  end do
  !
  DO k=2,NLAYM1
  DO I=1,IMJM
   c0 =  delp(i,1,k) / (delp(i,1,k-1)+delp(i,1,k)+delp(i,1,k+1))
   c1 = (delp(i,1,k-1)+0.5*delp(i,1,k))/(delp(i,1,k+1)+delp(i,1,k))
   c2 = (delp(i,1,k+1)+0.5*delp(i,1,k))/(delp(i,1,k-1)+delp(i,1,k))
  tmp = c0*(c1*DQDT(i,1,k) + c2*DQDT(i,1,k-1))
  Qmax = max(P(i,1,k-1),P(i,1,k),P(i,1,k+1)) - P(i,1,k)
  Qmin = P(i,1,k) - min(P(i,1,k-1),P(i,1,k),P(i,1,k+1))
  DC(i,1,k) = sign(min(abs(tmp),Qmax,Qmin), tmp)
  END DO
  END DO

  !
  ! ****6***0*********0*********0*********0*********0*********0**********72
  ! Loop over latitudes  (to save memory)
  ! ****6***0*********0*********0*********0*********0*********0**********72
  !
  DO j=1,JNP
  if((j.eq.2 .or. j.eq.JMR) .and. j1.ne.2) goto 2000
  !
  DO k=1,NLAY
  DO i=1,IMR
  wz2(i,k) =   WZ(i,j,k)
  wk1(i,k) =    P(i,j,k)
  wk2(i,k) = delp(i,j,k)
  flux(i,k) = DC(i,j,k)  !this flux is actually the monotone slope
  enddo
  enddo
  !
  !****6***0*********0*********0*********0*********0*********0**********72
  ! Compute first guesses at cell interfaces
  ! First guesses are required to be continuous.
  ! ****6***0*********0*********0*********0*********0*********0**********72
  !
  ! three-cell parabolic subgrid distribution at model top
  ! two-cell parabolic with zero gradient subgrid distribution
  ! at the surface.
  !
  ! First guess top edge value
  DO i=1,IMR
  ! three-cell PPM
  ! Compute a,b, and c of q = aP**2 + bP + c using cell averages and delp
  a = 3.*( DQDT(i,j,2) - DQDT(i,j,1)*(wk2(i,2)+wk2(i,3))/ &
        (wk2(i,1)+wk2(i,2)) ) / &
        ( (wk2(i,2)+wk2(i,3))*(wk2(i,1)+wk2(i,2)+wk2(i,3)) )
  b = 2.*DQDT(i,j,1)/(wk2(i,1)+wk2(i,2)) - &
        R23*a*(2.*wk2(i,1)+wk2(i,2))
  AL(i,1) =  wk1(i,1) - wk2(i,1)*(R3*a*wk2(i,1) + 0.5*b)
  AL(i,2) =  wk2(i,1)*(a*wk2(i,1) + b) + AL(i,1)
  !
  ! Check if change sign
  if(wk1(i,1)*AL(i,1).le.0.) then
             AL(i,1) = 0.
         flux(i,1) = 0.
    else
         flux(i,1) =  wk1(i,1) - AL(i,1)
    endif
  END DO
  !
  ! Bottom
  DO i=1,IMR
  ! 2-cell PPM with zero gradient right at the surface
  !
  fct = DQDT(i,j,NLAYM1)*wk2(i,NLAY)**2 / &
        ( (wk2(i,NLAY)+wk2(i,NLAYM1))*(2.*wk2(i,NLAY)+wk2(i,NLAYM1)))
  AR(i,NLAY) = wk1(i,NLAY) + fct
  AL(i,NLAY) = wk1(i,NLAY) - (fct+fct)
  if(wk1(i,NLAY)*AR(i,NLAY).le.0.) AR(i,NLAY) = 0.
  flux(i,NLAY) = AR(i,NLAY) -  wk1(i,NLAY)
  END DO

  !
  !****6***0*********0*********0*********0*********0*********0**********72
  ! 4th order interpolation in the interior.
  !****6***0*********0*********0*********0*********0*********0**********72
  !
  DO k=3,NLAYM1
  DO i=1,IMR
  c1 =  DQDT(i,j,k-1)*wk2(i,k-1) / (wk2(i,k-1)+wk2(i,k))
  c2 =  2. / (wk2(i,k-2)+wk2(i,k-1)+wk2(i,k)+wk2(i,k+1))
  A1   =  (wk2(i,k-2)+wk2(i,k-1)) / (2.*wk2(i,k-1)+wk2(i,k))
  A2   =  (wk2(i,k  )+wk2(i,k+1)) / (2.*wk2(i,k)+wk2(i,k-1))
  AL(i,k) = wk1(i,k-1) + c1 + c2 * &
        ( wk2(i,k  )*(c1*(A1 - A2)+A2*flux(i,k-1)) - &
        wk2(i,k-1)*A1*flux(i,k)  )
   ! print *,'AL1',i,k, AL(i,k)
  END DO
  END DO
  !
  do i=1,IMR*NLAYM1
  AR(i,1) = AL(i,2)
   ! print *,'AR1',i,AR(i,1)
  end do
  !
  do i=1,IMR*NLAY
  A6(i,1) = 3.*(wk1(i,1)+wk1(i,1) - (AL(i,1)+AR(i,1)))
   ! print *,'A61',i,A6(i,1)
  end do
  !
  !****6***0*********0*********0*********0*********0*********0**********72
  ! Top & Bot always monotonic
  call lmtppm(flux(1,1),A6(1,1),AR(1,1),AL(1,1),wk1(1,1),IMR,0)
  call lmtppm(flux(1,NLAY),A6(1,NLAY),AR(1,NLAY),AL(1,NLAY), &
        wk1(1,NLAY),IMR,0)
  !
  ! Interior depending on KORD
  if(LMT.LE.2) &
        call lmtppm(flux(1,2),A6(1,2),AR(1,2),AL(1,2),wk1(1,2), &
        IMR*(NLAY-2),LMT)
  !
  !****6***0*********0*********0*********0*********0*********0**********72
  !
  DO i=1,IMR*NLAYM1
  IF(wz2(i,1).GT.0.) then
    CM = wz2(i,1) / wk2(i,1)
    flux(i,2) = AR(i,1)+0.5*CM*(AL(i,1)-AR(i,1)+A6(i,1)*(1.-R23*CM))
  else
     ! print *,'test2-0',i,j,wz2(i,1),wk2(i,2)
    CP= wz2(i,1) / wk2(i,2)
     ! print *,'testCP',CP
    flux(i,2) = AL(i,2)+0.5*CP*(AL(i,2)-AR(i,2)-A6(i,2)*(1.+R23*CP))
     ! print *,'test2',i, AL(i,2),AR(i,2),A6(i,2),R23
  endif
  END DO
  !
  DO i=1,IMR*NLAYM1
  flux(i,2) = wz2(i,1) * flux(i,2)
250 CONTINUE
  END DO
  !
  do i=1,IMR
  DQ(i,j,   1) = DQ(i,j,   1) - flux(i,   2)
  DQ(i,j,NLAY) = DQ(i,j,NLAY) + flux(i,NLAY)
  end do
  !
  do k=2,NLAYM1
  do i=1,IMR
    DQ(i,j,k) = DQ(i,j,k) + flux(i,k) - flux(i,k+1)
  end do
  end do
2000 CONTINUE
  END DO
  return
end subroutine fzppm
!
subroutine xtp(IMR,JNP,IML,j1,j2,JN,JS,PU,DQ,Q,UC, &
        fx1,xmass,IORD)
  implicit none
  integer :: IMR,JNP,IML,j1,j2,JN,JS,IORD
  real :: PU,DQ,Q,UC,fx1,xmass
  real :: dc,qtmp
  integer :: ISAVE(IMR)
  dimension UC(IMR,*),DC(-IML:IMR+IML+1),xmass(IMR,JNP) &
        ,fx1(IMR+1),DQ(IMR,JNP),qtmp(-IML:IMR+1+IML)
  dimension PU(IMR,JNP),Q(IMR,JNP)
  integer :: jvan,j1vl,j2vl,j,i,iu,itmp,ist,imp
  real :: rut
  !
  IMP = IMR + 1
  !
  ! van Leer at high latitudes
  jvan = max(1,JNP/18)
  j1vl = j1+jvan
  j2vl = j2-jvan
  !
  do j=j1,j2
  !
  do i=1,IMR
  qtmp(i) = q(i,j)
  enddo
  !
  if(j.ge.JN .or. j.le.JS) goto 2222
  ! ************* Eulerian **********
  !
  qtmp(0)     = q(IMR,J)
  qtmp(-1)    = q(IMR-1,J)
  qtmp(IMP)   = q(1,J)
  qtmp(IMP+1) = q(2,J)
  !
  IF(IORD.eq.1 .or. j.eq.j1.OR. j.eq.j2) THEN
  DO i=1,IMR
  iu = REAL(i) - uc(i,j)
    fx1(i) = qtmp(iu)
  END DO
  ELSE
  call xmist(IMR,IML,Qtmp,DC)
  DC(0) = DC(IMR)
  !
  if(IORD.eq.2 .or. j.le.j1vl .or. j.ge.j2vl) then
  DO i=1,IMR
  iu = REAL(i) - uc(i,j)
    fx1(i) = qtmp(iu) + DC(iu)*(sign(1.,uc(i,j))-uc(i,j))
  END DO
  else
  call fxppm(IMR,IML,UC(1,j),Qtmp,DC,fx1,IORD)
  endif
  !
  ENDIF
  !
  DO i=1,IMR
    fx1(i) = fx1(i)*xmass(i,j)
  END DO
  !
  goto 1309
  !
  ! ***** Conservative (flux-form) Semi-Lagrangian transport *****
  !
2222   continue
  !
  do i=-IML,0
  qtmp(i)     = q(IMR+i,j)
  qtmp(IMP-i) = q(1-i,j)
  enddo
  !
  IF(IORD.eq.1 .or. j.eq.j1.OR. j.eq.j2) THEN
  DO i=1,IMR
  itmp = INT(uc(i,j))
  ISAVE(i) = i - itmp
  iu = i - uc(i,j)
    fx1(i) = (uc(i,j) - itmp)*qtmp(iu)
  END DO
  ELSE
  call xmist(IMR,IML,Qtmp,DC)
  !
  do i=-IML,0
  DC(i)     = DC(IMR+i)
  DC(IMP-i) = DC(1-i)
  enddo
  !
  DO i=1,IMR
  itmp = INT(uc(i,j))
  rut  = uc(i,j) - itmp
  ISAVE(i) = i - itmp
  iu = i - uc(i,j)
    fx1(i) = rut*(qtmp(iu) + DC(iu)*(sign(1.,rut) - rut))
  END DO
  ENDIF
  !
  do i=1,IMR
  IF(uc(i,j).GT.1.) then
  !DIR$ NOVECTOR
    do ist = ISAVE(i),i-1
    fx1(i) = fx1(i) + qtmp(ist)
    enddo
  elseIF(uc(i,j).LT.-1.) then
    do ist = i,ISAVE(i)-1
    fx1(i) = fx1(i) - qtmp(ist)
    enddo
  !DIR$ VECTOR
  endif
  end do
  do i=1,IMR
  fx1(i) = PU(i,j)*fx1(i)
  enddo
  !
  ! ***************************************
  !
1309   fx1(IMP) = fx1(1)
  DO i=1,IMR
    DQ(i,j) =  DQ(i,j) + fx1(i)-fx1(i+1)
  END DO
  !
  ! ***************************************
  !
  end do
  return
end subroutine xtp
!
subroutine fxppm(IMR,IML,UT,P,DC,flux,IORD)
  implicit none
  integer :: IMR,IML,IORD
  real :: UT,P,DC,flux
  real,parameter ::  R3 = 1./3., R23 = 2./3.
  DIMENSION UT(*),flux(*),P(-IML:IMR+IML+1),DC(-IML:IMR+IML+1)
  REAL :: AR(0:IMR),AL(0:IMR),A6(0:IMR)
  integer :: LMT,IMP,JLVL,i
   ! logical first
   ! data first /.true./
   ! SAVE LMT
   ! if(first) then
  !
  ! correction calcul de LMT a chaque passage pour pouvoir choisir
  ! plusieurs schemas PPM pour differents traceurs
  !  IF (IORD.LE.0) then
  !        if(IMR.GE.144) then
  !              LMT = 0
  !        elseif(IMR.GE.72) then
  !              LMT = 1
  !        else
  !              LMT = 2
  !        endif
  !  else
  !        LMT = IORD - 3
  !  endif
  !
  LMT = IORD - 3
   ! write(6,*) 'PPM option in E-W direction = ', LMT
   ! first = .false.
   ! endif
  !
  DO i=1,IMR
    AL(i) = 0.5*(p(i-1)+p(i)) + (DC(i-1) - DC(i))*R3
  END DO
  !
  do i=1,IMR-1
    AR(i) = AL(i+1)
  end do
  AR(IMR) = AL(1)
  !
  do i=1,IMR
    A6(i) = 3.*(p(i)+p(i)  - (AL(i)+AR(i)))
  end do
  !
  if(LMT.LE.2) call lmtppm(DC(1),A6(1),AR(1),AL(1),P(1),IMR,LMT)
  !
  AL(0) = AL(IMR)
  AR(0) = AR(IMR)
  A6(0) = A6(IMR)
  !
  DO i=1,IMR
  IF(UT(i).GT.0.) then
  flux(i) = AR(i-1) + 0.5*UT(i)*(AL(i-1) - AR(i-1) + &
        A6(i-1)*(1.-R23*UT(i)) )
  else
  flux(i) = AL(i) - 0.5*UT(i)*(AR(i) - AL(i) + &
        A6(i)*(1.+R23*UT(i)))
  endif
  enddo
  return
end subroutine fxppm
!
subroutine xmist(IMR,IML,P,DC)
  implicit none
  integer :: IMR,IML
  real,parameter :: R24 = 1./24.
  real :: P(-IML:IMR+1+IML),DC(-IML:IMR+1+IML)
  integer :: i
  real :: tmp,pmax,pmin
  !
  do  i=1,IMR
  tmp = R24*(8.*(p(i+1) - p(i-1)) + p(i-2) - p(i+2))
  Pmax = max(P(i-1), p(i), p(i+1)) - p(i)
  Pmin = p(i) - min(P(i-1), p(i), p(i+1))
    DC(i) = sign(min(abs(tmp),Pmax,Pmin), tmp)
  end do
  return
end subroutine xmist
!
subroutine ytp(IMR,JNP,j1,j2,acosp,RCAP,DQ,P,VC,DC2 &
        ,ymass,fx,A6,AR,AL,JORD)
  implicit none
  integer :: IMR,JNP,j1,j2,JORD
  real :: acosp,RCAP,DQ,P,VC,DC2,ymass,fx,A6,AR,AL
  dimension P(IMR,JNP),VC(IMR,JNP),ymass(IMR,JNP) &
        ,DC2(IMR,JNP),DQ(IMR,JNP),acosp(JNP)
  ! Work array
  DIMENSION fx(IMR,JNP),AR(IMR,JNP),AL(IMR,JNP),A6(IMR,JNP)
  integer :: JMR,len,i,jt,j
  real :: sum1,sum2
  !
  JMR = JNP - 1
  len = IMR*(J2-J1+2)
  !
  if(JORD.eq.1) then
  DO i=1,len
  JT = REAL(J1) - VC(i,J1)
    fx(i,j1) = p(i,JT)
  END DO
  else

  call ymist(IMR,JNP,j1,P,DC2,4)
  !
  if(JORD.LE.0 .or. JORD.GE.3) then

  call fyppm(VC,P,DC2,fx,IMR,JNP,j1,j2,A6,AR,AL,JORD)

  else
  DO i=1,len
  JT = REAL(J1) - VC(i,J1)
    fx(i,j1) = p(i,JT) + (sign(1.,VC(i,j1))-VC(i,j1))*DC2(i,JT)
  END DO
  endif
  endif
  !
  DO i=1,len
    fx(i,j1) = fx(i,j1)*ymass(i,j1)
  END DO
  !
  DO j=j1,j2
  DO i=1,IMR
    DQ(i,j) = DQ(i,j) + (fx(i,j) - fx(i,j+1)) * acosp(j)
  END DO
  END DO
  !
  ! Poles
  sum1 = fx(IMR,j1  )
  sum2 = fx(IMR,J2+1)
  do i=1,IMR-1
  sum1 = sum1 + fx(i,j1  )
  sum2 = sum2 + fx(i,J2+1)
  enddo
  !
  sum1 = DQ(1,  1) - sum1 * RCAP
  sum2 = DQ(1,JNP) + sum2 * RCAP
  do i=1,IMR
  DQ(i,  1) = sum1
  DQ(i,JNP) = sum2
  enddo
  !
  if(j1.ne.2) then
  do i=1,IMR
  DQ(i,  2) = sum1
  DQ(i,JMR) = sum2
  enddo
  endif
  !
  return
end subroutine ytp
!
subroutine  ymist(IMR,JNP,j1,P,DC,ID)
  implicit none
  integer :: IMR,JNP,j1,ID
  real,parameter :: R24 = 1./24.
  real :: P(IMR,JNP),DC(IMR,JNP)
  integer :: iimh,jmr,ijm3,imh,i
  real :: pmax,pmin,tmp
  !
  IMH = IMR / 2
  JMR = JNP - 1
  IJM3 = IMR*(JMR-3)
  !
  IF(ID.EQ.2) THEN
  do i=1,IMR*(JMR-1)
  tmp = 0.25*(p(i,3) - p(i,1))
  Pmax = max(p(i,1),p(i,2),p(i,3)) - p(i,2)
  Pmin = p(i,2) - min(p(i,1),p(i,2),p(i,3))
  DC(i,2) = sign(min(abs(tmp),Pmin,Pmax),tmp)
  end do
  ELSE
  do i=1,IMH
  ! J=2
  tmp = (8.*(p(i,3) - p(i,1)) + p(i+IMH,2) - p(i,4))*R24
  Pmax = max(p(i,1),p(i,2),p(i,3)) - p(i,2)
  Pmin = p(i,2) - min(p(i,1),p(i,2),p(i,3))
  DC(i,2) = sign(min(abs(tmp),Pmin,Pmax),tmp)
  ! J=JMR
  tmp=(8.*(p(i,JNP)-p(i,JMR-1))+p(i,JMR-2)-p(i+IMH,JMR))*R24
  Pmax = max(p(i,JMR-1),p(i,JMR),p(i,JNP)) - p(i,JMR)
  Pmin = p(i,JMR) - min(p(i,JMR-1),p(i,JMR),p(i,JNP))
  DC(i,JMR) = sign(min(abs(tmp),Pmin,Pmax),tmp)
  end do
  do i=IMH+1,IMR
  ! J=2
  tmp = (8.*(p(i,3) - p(i,1)) + p(i-IMH,2) - p(i,4))*R24
  Pmax = max(p(i,1),p(i,2),p(i,3)) - p(i,2)
  Pmin = p(i,2) - min(p(i,1),p(i,2),p(i,3))
  DC(i,2) = sign(min(abs(tmp),Pmin,Pmax),tmp)
  ! J=JMR
  tmp=(8.*(p(i,JNP)-p(i,JMR-1))+p(i,JMR-2)-p(i-IMH,JMR))*R24
  Pmax = max(p(i,JMR-1),p(i,JMR),p(i,JNP)) - p(i,JMR)
  Pmin = p(i,JMR) - min(p(i,JMR-1),p(i,JMR),p(i,JNP))
  DC(i,JMR) = sign(min(abs(tmp),Pmin,Pmax),tmp)
  end do
  !
  do i=1,IJM3
  tmp = (8.*(p(i,4) - p(i,2)) + p(i,1) - p(i,5))*R24
  Pmax = max(p(i,2),p(i,3),p(i,4)) - p(i,3)
  Pmin = p(i,3) - min(p(i,2),p(i,3),p(i,4))
  DC(i,3) = sign(min(abs(tmp),Pmin,Pmax),tmp)
  end do
  ENDIF
  !
  if(j1.ne.2) then
  do i=1,IMR
  DC(i,1) = 0.
  DC(i,JNP) = 0.
  enddo
  else
  ! Determine slopes in polar caps for scalars!
  !
  do i=1,IMH
  ! South
  tmp = 0.25*(p(i,2) - p(i+imh,2))
  Pmax = max(p(i,2),p(i,1), p(i+imh,2)) - p(i,1)
  Pmin = p(i,1) - min(p(i,2),p(i,1), p(i+imh,2))
  DC(i,1)=sign(min(abs(tmp),Pmax,Pmin),tmp)
  ! North.
  tmp = 0.25*(p(i+imh,JMR) - p(i,JMR))
  Pmax = max(p(i+imh,JMR),p(i,jnp), p(i,JMR)) - p(i,JNP)
  Pmin = p(i,JNP) - min(p(i+imh,JMR),p(i,jnp), p(i,JMR))
  DC(i,JNP) = sign(min(abs(tmp),Pmax,pmin),tmp)
  end do
  !
  do i=imh+1,IMR
  DC(i,  1) =  - DC(i-imh,  1)
  DC(i,JNP) =  - DC(i-imh,JNP)
  end do
  endif
  return
end subroutine ymist
!
subroutine fyppm(VC,P,DC,flux,IMR,JNP,j1,j2,A6,AR,AL,JORD)
  implicit none
  integer :: IMR,JNP,j1,j2,JORD
  real,parameter :: R3 = 1./3., R23 = 2./3.
  real :: VC(IMR,*),flux(IMR,*),P(IMR,*),DC(IMR,*)
  ! Local work arrays.
  real :: AR(IMR,JNP),AL(IMR,JNP),A6(IMR,JNP)
  integer :: LMT,i
  integer :: IMH,JMR,j11,IMJM1,len
   ! logical first
   ! data first /.true./
   ! SAVE LMT
  !
  IMH = IMR / 2
  JMR = JNP - 1
  j11 = j1-1
  IMJM1 = IMR*(J2-J1+2)
  len   = IMR*(J2-J1+3)
   ! if(first) then
   ! IF(JORD.LE.0) then
   !       if(JMR.GE.90) then
   !             LMT = 0
   !       elseif(JMR.GE.45) then
   !             LMT = 1
   !       else
   !             LMT = 2
   !       endif
   ! else
   !       LMT = JORD - 3
   ! endif
  !
  !  first = .false.
  !  endif
  !
  ! modifs pour pouvoir choisir plusieurs schemas PPM
  LMT = JORD - 3
  !
  DO i=1,IMR*JMR
  AL(i,2) = 0.5*(p(i,1)+p(i,2)) + (DC(i,1) - DC(i,2))*R3
  AR(i,1) = AL(i,2)
  END DO
  !
  !Poles:
  !
  DO i=1,IMH
  AL(i,1) = AL(i+IMH,2)
  AL(i+IMH,1) = AL(i,2)
  !
  AR(i,JNP) = AR(i+IMH,JMR)
  AR(i+IMH,JNP) = AR(i,JMR)
  ENDDO

  !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
  !   Rajout pour LMDZ.3.3
  !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
  AR(IMR,1)=AL(1,1)
  AR(IMR,JNP)=AL(1,JNP)
  !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc


  do i=1,len
    A6(i,j11) = 3.*(p(i,j11)+p(i,j11)  - (AL(i,j11)+AR(i,j11)))
  end do
  !
  if(LMT.le.2) call lmtppm(DC(1,j11),A6(1,j11),AR(1,j11) &
        ,AL(1,j11),P(1,j11),len,LMT)
  !

  DO i=1,IMJM1
  IF(VC(i,j1).GT.0.) then
  flux(i,j1) = AR(i,j11) + 0.5*VC(i,j1)*(AL(i,j11) - AR(i,j11) + &
        A6(i,j11)*(1.-R23*VC(i,j1)) )
  else
  flux(i,j1) = AL(i,j1) - 0.5*VC(i,j1)*(AR(i,j1) - AL(i,j1) + &
        A6(i,j1)*(1.+R23*VC(i,j1)))
  endif
  END DO
  return
end subroutine fyppm
!
  subroutine yadv(IMR,JNP,j1,j2,p,VA,ady,wk,IAD)
    implicit none
    integer :: IMR,JNP,j1,j2,IAD
    REAL :: p(IMR,JNP),ady(IMR,JNP),VA(IMR,JNP)
    REAL :: WK(IMR,-1:JNP+2)
    INTEGER :: JMR,IMH,i,j,jp
    REAL :: rv,a1,b1,sum1,sum2
  !
    JMR = JNP-1
    IMH = IMR/2
    do j=1,JNP
    do i=1,IMR
    wk(i,j) = p(i,j)
    enddo
    enddo
  ! Poles:
    do i=1,IMH
    wk(i,   -1) = p(i+IMH,3)
    wk(i+IMH,-1) = p(i,3)
    wk(i,    0) = p(i+IMH,2)
    wk(i+IMH,0) = p(i,2)
    wk(i,JNP+1) = p(i+IMH,JMR)
    wk(i+IMH,JNP+1) = p(i,JMR)
    wk(i,JNP+2) = p(i+IMH,JNP-2)
    wk(i+IMH,JNP+2) = p(i,JNP-2)
    enddo
     ! write(*,*) 'toto 1'
  ! --------------------------------
  IF(IAD.eq.2) then
  do j=j1-1,j2+1
  do i=1,IMR
   ! write(*,*) 'avt NINT','i=',i,'j=',j
  JP = NINT(VA(i,j))
  rv = JP - VA(i,j)
   ! write(*,*) 'VA=',VA(i,j), 'JP1=',JP,'rv=',rv
  JP = j - JP
   ! write(*,*) 'JP2=',JP
  a1 = 0.5*(wk(i,jp+1)+wk(i,jp-1)) - wk(i,jp)
  b1 = 0.5*(wk(i,jp+1)-wk(i,jp-1))
   ! write(*,*) 'a1=',a1,'b1=',b1
  ady(i,j) = wk(i,jp) + rv*(a1*rv + b1) - wk(i,j)
  enddo
  enddo
   ! write(*,*) 'toto 2'
  !
  ELSEIF(IAD.eq.1) then
    do j=j1-1,j2+1
  do i=1,imr
  JP = REAL(j)-VA(i,j)
  ady(i,j) = VA(i,j)*(wk(i,jp)-wk(i,jp+1))
  enddo
  enddo
  ENDIF
  !
    if(j1.ne.2) then
    sum1 = 0.
    sum2 = 0.
  do i=1,imr
  sum1 = sum1 + ady(i,2)
  sum2 = sum2 + ady(i,JMR)
  enddo
    sum1 = sum1 / IMR
    sum2 = sum2 / IMR
  !
  do i=1,imr
  ady(i,  2) =  sum1
  ady(i,JMR) =  sum2
  ady(i,  1) =  sum1
  ady(i,JNP) =  sum2
  enddo
    else
  ! Poles:
    sum1 = 0.
    sum2 = 0.
  do i=1,imr
  sum1 = sum1 + ady(i,1)
  sum2 = sum2 + ady(i,JNP)
  enddo
    sum1 = sum1 / IMR
    sum2 = sum2 / IMR
  !
  do i=1,imr
  ady(i,  1) =  sum1
  ady(i,JNP) =  sum2
  enddo
    endif
  !
    return
end subroutine yadv
!
  subroutine xadv(IMR,JNP,j1,j2,p,UA,JS,JN,IML,adx,IAD)
    implicit none
    INTEGER :: IMR,JNP,j1,j2,JS,JN,IML,IAD
    REAL :: p(IMR,JNP),adx(IMR,JNP),qtmp(-IMR:IMR+IMR),UA(IMR,JNP)
    INTEGER :: JMR,j,i,ip,iu,iiu
    REAL :: ru,a1,b1
  !
    JMR = JNP-1
  do j=j1,j2
  if(J.GT.JS  .and. J.LT.JN) GO TO 1309
  !
  do i=1,IMR
  qtmp(i) = p(i,j)
  enddo
  !
  do i=-IML,0
  qtmp(i)       = p(IMR+i,j)
  qtmp(IMR+1-i) = p(1-i,j)
  enddo
  !
  IF(IAD.eq.2) THEN
  DO i=1,IMR
  IP = NINT(UA(i,j))
  ru = IP - UA(i,j)
  IP = i - IP
  a1 = 0.5*(qtmp(ip+1)+qtmp(ip-1)) - qtmp(ip)
  b1 = 0.5*(qtmp(ip+1)-qtmp(ip-1))
  adx(i,j) = qtmp(ip) + ru*(a1*ru + b1)
  enddo
  ELSEIF(IAD.eq.1) then
  DO i=1,IMR
  iu = UA(i,j)
  ru = UA(i,j) - iu
  iiu = i-iu
  if(UA(i,j).GE.0.) then
  adx(i,j) = qtmp(iiu)+ru*(qtmp(iiu-1)-qtmp(iiu))
  else
  adx(i,j) = qtmp(iiu)+ru*(qtmp(iiu)-qtmp(iiu+1))
  endif
  enddo
  ENDIF
  !
  do i=1,IMR
  adx(i,j) = adx(i,j) - p(i,j)
  enddo
1309 continue
  end do
  !
  ! Eulerian upwind
  !
  do j=JS+1,JN-1
  !
  do i=1,IMR
  qtmp(i) = p(i,j)
  enddo
  !
  qtmp(0)     = p(IMR,J)
  qtmp(IMR+1) = p(1,J)
  !
  IF(IAD.eq.2) THEN
  qtmp(-1)     = p(IMR-1,J)
  qtmp(IMR+2) = p(2,J)
  do i=1,imr
  IP = NINT(UA(i,j))
  ru = IP - UA(i,j)
  IP = i - IP
  a1 = 0.5*(qtmp(ip+1)+qtmp(ip-1)) - qtmp(ip)
  b1 = 0.5*(qtmp(ip+1)-qtmp(ip-1))
  adx(i,j) = qtmp(ip)- p(i,j) + ru*(a1*ru + b1)
  enddo
  ELSEIF(IAD.eq.1) then
  ! 1st order
  DO i=1,IMR
  IP = i - UA(i,j)
  adx(i,j) = UA(i,j)*(qtmp(ip)-qtmp(ip+1))
  enddo
  ENDIF
  enddo
  !
    if(j1.ne.2) then
  do i=1,IMR
  adx(i,  2) = 0.
  adx(i,JMR) = 0.
  enddo
    endif
  ! set cross term due to x-adv at the poles to zero.
  do i=1,IMR
  adx(i,  1) = 0.
  adx(i,JNP) = 0.
  enddo
    return
end subroutine xadv
!
subroutine lmtppm(DC,A6,AR,AL,P,IM,LMT)
  implicit none
  !
  ! A6 =  CURVATURE OF THE TEST PARABOLA
  ! AR =  RIGHT EDGE VALUE OF THE TEST PARABOLA
  ! AL =  LEFT  EDGE VALUE OF THE TEST PARABOLA
  ! DC =  0.5 * MISMATCH
  ! P  =  CELL-AVERAGED VALUE
  ! IM =  VECTOR LENGTH
  !
  ! OPTIONS:
  !
  ! LMT = 0: FULL MONOTONICITY
  ! LMT = 1: SEMI-MONOTONIC CONSTRAINT (NO UNDERSHOOTS)
  ! LMT = 2: POSITIVE-DEFINITE CONSTRAINT
  !
  real,parameter :: R12 = 1./12.
  real :: A6(IM),AR(IM),AL(IM),P(IM),DC(IM)
  integer :: IM,LMT
  INTEGER :: i
  REAL :: da1,da2,a6da,fmin
  !
  if(LMT.eq.0) then
  ! Full constraint
  do i=1,IM
  if(DC(i).eq.0.) then
        AR(i) = p(i)
        AL(i) = p(i)
        A6(i) = 0.
  else
  da1  = AR(i) - AL(i)
  da2  = da1**2
  A6DA = A6(i)*da1
  if(A6DA .lt. -da2) then
        A6(i) = 3.*(AL(i)-p(i))
        AR(i) = AL(i) - A6(i)
  elseif(A6DA .gt. da2) then
        A6(i) = 3.*(AR(i)-p(i))
        AL(i) = AR(i) - A6(i)
  endif
  endif
  end do
  elseif(LMT.eq.1) then
  ! Semi-monotonic constraint
  do i=1,IM
  if(abs(AR(i)-AL(i)) .GE. -A6(i)) go to 150
  if(p(i).lt.AR(i) .and. p(i).lt.AL(i)) then
        AR(i) = p(i)
        AL(i) = p(i)
        A6(i) = 0.
  elseif(AR(i) .gt. AL(i)) then
        A6(i) = 3.*(AL(i)-p(i))
        AR(i) = AL(i) - A6(i)
  else
        A6(i) = 3.*(AR(i)-p(i))
        AL(i) = AR(i) - A6(i)
  endif
150 continue
  end do
  elseif(LMT.eq.2) then
  do i=1,IM
  if(abs(AR(i)-AL(i)) .GE. -A6(i)) go to 250
  fmin = p(i) + 0.25*(AR(i)-AL(i))**2/A6(i) + A6(i)*R12
  if(fmin.ge.0.) go to 250
  if(p(i).lt.AR(i) .and. p(i).lt.AL(i)) then
        AR(i) = p(i)
        AL(i) = p(i)
        A6(i) = 0.
  elseif(AR(i) .gt. AL(i)) then
        A6(i) = 3.*(AL(i)-p(i))
        AR(i) = AL(i) - A6(i)
  else
        A6(i) = 3.*(AR(i)-p(i))
        AL(i) = AR(i) - A6(i)
  endif
250 continue
  end do
  endif
  return
end subroutine lmtppm
!
subroutine A2C(U,V,IMR,JMR,j1,j2,CRX,CRY,dtdx5,DTDY5)
  implicit none
  integer :: IMR,JMR,j1,j2
  real :: U(IMR,*),V(IMR,*),CRX(IMR,*),CRY(IMR,*),DTDX5(*),DTDY5
  integer :: i,j
  !
  do j=j1,j2
  do i=2,IMR
    CRX(i,J) = dtdx5(j)*(U(i,j)+U(i-1,j))
  end do
  end do
  !
  do j=j1,j2
    CRX(1,J) = dtdx5(j)*(U(1,j)+U(IMR,j))
  end do
  !
  do i=1,IMR*JMR
    CRY(i,2) = DTDY5*(V(i,2)+V(i,1))
  end do
  return
end subroutine a2c
!
subroutine cosa(cosp,cose,JNP,PI,DP)
  implicit none
  integer :: JNP
  real :: cosp(*),cose(*),PI,DP
  integer :: JMR,j,jeq
  real :: ph5
  JMR = JNP-1
  do j=2,JNP
    ph5  =  -0.5*PI + (REAL(J-1)-0.5)*DP
    cose(j) = cos(ph5)
  end do
  !
  JEQ = (JNP+1) / 2
  if(JMR .eq. 2*(JMR/2) ) then
  do j=JNP, JEQ+1, -1
   cose(j) =  cose(JNP+2-j)
  enddo
  else
  ! cell edge at equator.
   cose(JEQ+1) =  1.
  do j=JNP, JEQ+2, -1
   cose(j) =  cose(JNP+2-j)
   enddo
  endif
  !
  do j=2,JMR
    cosp(j) = 0.5*(cose(j)+cose(j+1))
  end do
  cosp(1) = 0.
  cosp(JNP) = 0.
  return
end subroutine cosa
!
subroutine cosc(cosp,cose,JNP,PI,DP)
  implicit none
  integer :: JNP
  real :: cosp(*),cose(*),PI,DP
  real :: phi
  integer :: j
  !
  phi = -0.5*PI
  do j=2,JNP-1
  phi  =  phi + DP
    cosp(j) = cos(phi)
  end do
    cosp(  1) = 0.
    cosp(JNP) = 0.
  !
  do j=2,JNP
    cose(j) = 0.5*(cosp(j)+cosp(j-1))
  end do
  !
  do j=2,JNP-1
   cosp(j) = 0.5*(cose(j)+cose(j+1))
  end do
  return
end subroutine cosc
!
SUBROUTINE qckxyz (Q,qtmp,IMR,JNP,NLAY,j1,j2,cosp,acosp, &
        cross,IC,NSTEP)
  !
  real,parameter :: tiny = 1.E-60
  INTEGER :: IMR,JNP,NLAY,j1,j2,IC,NSTEP
  REAL :: Q(IMR,JNP,NLAY),qtmp(IMR,JNP),cosp(*),acosp(*)
  logical :: cross
  INTEGER :: NLAYM1,len,ip,L,icr,ipy,ipx,i
  real :: qup,qly,dup,sum
  !
  NLAYM1 = NLAY-1
  len = IMR*(j2-j1+1)
  ip = 0
  !
  ! Top layer
  L = 1
    icr = 1
  call filns(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,ipy,tiny)
  if(ipy.eq.0) goto 50
  call filew(q(1,1,L),qtmp,IMR,JNP,j1,j2,ipx,tiny)
  if(ipx.eq.0) goto 50
  !
  if(cross) then
  call filcr(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,icr,tiny)
  endif
  if(icr.eq.0) goto 50
  !
  ! Vertical filling...
  do i=1,len
  IF( Q(i,j1,1).LT.0.) THEN
  ip = ip + 1
      Q(i,j1,2) = Q(i,j1,2) + Q(i,j1,1)
      Q(i,j1,1) = 0.
  endif
  enddo
  !
50   continue
  DO L = 2,NLAYM1
  icr = 1
  !
  call filns(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,ipy,tiny)
  if(ipy.eq.0) goto 225
  call filew(q(1,1,L),qtmp,IMR,JNP,j1,j2,ipx,tiny)
  if(ipx.eq.0) go to 225
  if(cross) then
  call filcr(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,icr,tiny)
  endif
  if(icr.eq.0) goto 225
  !
  do i=1,len
  IF( Q(I,j1,L).LT.0.) THEN
  !
  ip = ip + 1
  ! From above
      qup =  Q(I,j1,L-1)
      qly = -Q(I,j1,L)
      dup  = min(qly,qup)
      Q(I,j1,L-1) = qup - dup
      Q(I,j1,L  ) = dup-qly
  ! Below
      Q(I,j1,L+1) = Q(I,j1,L+1) + Q(I,j1,L)
      Q(I,j1,L)   = 0.
  ENDIF
  ENDDO
225 CONTINUE
  END DO
  !
  ! BOTTOM LAYER
  sum = 0.
  L = NLAY
  !
  call filns(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,ipy,tiny)
  if(ipy.eq.0) goto 911
  call filew(q(1,1,L),qtmp,IMR,JNP,j1,j2,ipx,tiny)
  if(ipx.eq.0) goto 911
  !
  call filcr(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,icr,tiny)
  if(icr.eq.0) goto 911
  !
  DO  I=1,len
  IF( Q(I,j1,L).LT.0.) THEN
  ip = ip + 1
  !
  ! From above
  !
      qup = Q(I,j1,NLAYM1)
      qly = -Q(I,j1,L)
      dup = min(qly,qup)
      Q(I,j1,NLAYM1) = qup - dup
  ! From "below" the surface.
      sum = sum + qly-dup
      Q(I,j1,L) = 0.
   ENDIF
  ENDDO
  !
911   continue
  !
  if(ip.gt.IMR) then
  write(6,*) 'IC=',IC,' STEP=',NSTEP, &
        ' Vertical filling pts=',ip
  endif
  !
  if(sum.gt.1.e-25) then
  write(6,*) IC,NSTEP,' Mass source from the ground=',sum
  endif
  RETURN
END SUBROUTINE qckxyz
!
subroutine filcr(q,IMR,JNP,j1,j2,cosp,acosp,icr,tiny)
  implicit none
  integer :: IMR,JNP,j1,j2,icr
  real :: q(IMR,*),cosp(*),acosp(*),tiny
  integer :: i,j
  real :: dq,dn,d0,d1,ds,d2
  icr = 0
  do j=j1+1,j2-1
  DO i=1,IMR-1
  IF(q(i,j).LT.0.) THEN
  icr =  1
  dq  = - q(i,j)*cosp(j)
  ! N-E
  dn = q(i+1,j+1)*cosp(j+1)
  d0 = max(0.,dn)
  d1 = min(dq,d0)
  q(i+1,j+1) = (dn - d1)*acosp(j+1)
  dq = dq - d1
  ! S-E
  ds = q(i+1,j-1)*cosp(j-1)
  d0 = max(0.,ds)
  d2 = min(dq,d0)
  q(i+1,j-1) = (ds - d2)*acosp(j-1)
  q(i,j) = (d2 - dq)*acosp(j) + tiny
  endif
50  CONTINUE
  END DO
  if(icr.eq.0 .and. q(IMR,j).ge.0.) goto 65
  DO i=2,IMR
  IF(q(i,j).LT.0.) THEN
  icr =  1
  dq  = - q(i,j)*cosp(j)
  ! N-W
  dn = q(i-1,j+1)*cosp(j+1)
  d0 = max(0.,dn)
  d1 = min(dq,d0)
  q(i-1,j+1) = (dn - d1)*acosp(j+1)
  dq = dq - d1
  ! S-W
  ds = q(i-1,j-1)*cosp(j-1)
  d0 = max(0.,ds)
  d2 = min(dq,d0)
  q(i-1,j-1) = (ds - d2)*acosp(j-1)
  q(i,j) = (d2 - dq)*acosp(j) + tiny
  endif
  END DO
  ! *****************************************
  ! i=1
  i=1
  IF(q(i,j).LT.0.) THEN
  icr =  1
  dq  = - q(i,j)*cosp(j)
  ! N-W
  dn = q(IMR,j+1)*cosp(j+1)
  d0 = max(0.,dn)
  d1 = min(dq,d0)
  q(IMR,j+1) = (dn - d1)*acosp(j+1)
  dq = dq - d1
  ! S-W
  ds = q(IMR,j-1)*cosp(j-1)
  d0 = max(0.,ds)
  d2 = min(dq,d0)
  q(IMR,j-1) = (ds - d2)*acosp(j-1)
  q(i,j) = (d2 - dq)*acosp(j) + tiny
  endif
  ! *****************************************
  ! i=IMR
  i=IMR
  IF(q(i,j).LT.0.) THEN
  icr =  1
  dq  = - q(i,j)*cosp(j)
  ! N-E
  dn = q(1,j+1)*cosp(j+1)
  d0 = max(0.,dn)
  d1 = min(dq,d0)
  q(1,j+1) = (dn - d1)*acosp(j+1)
  dq = dq - d1
  ! S-E
  ds = q(1,j-1)*cosp(j-1)
  d0 = max(0.,ds)
  d2 = min(dq,d0)
  q(1,j-1) = (ds - d2)*acosp(j-1)
  q(i,j) = (d2 - dq)*acosp(j) + tiny
  endif
  ! *****************************************
65  continue
  end do
  !
  do i=1,IMR
  if(q(i,j1).lt.0. .or. q(i,j2).lt.0.) then
  icr = 1
  goto 80
  endif
  enddo
  !
80   continue
  !
  if(q(1,1).lt.0. .or. q(1,jnp).lt.0.) then
  icr = 1
  endif
  !
  return
end subroutine filcr
!
subroutine filns(q,IMR,JNP,j1,j2,cosp,acosp,ipy,tiny)
  implicit none
  integer :: IMR,JNP,j1,j2,ipy
  real :: q(IMR,*),cosp(*),acosp(*),tiny
  real :: DP,CAP1,dq,dn,d0,d1,ds,d2
  INTEGER :: i,j
   ! logical first
   ! data first /.true./
   ! save cap1
  !
  !  if(first) then
  DP = 4.*ATAN(1.)/REAL(JNP-1)
  CAP1 = IMR*(1.-COS((j1-1.5)*DP))/DP
   ! first = .false.
   ! endif
  !
  ipy = 0
  do j=j1+1,j2-1
  DO i=1,IMR
  IF(q(i,j).LT.0.) THEN
  ipy =  1
  dq  = - q(i,j)*cosp(j)
  ! North
  dn = q(i,j+1)*cosp(j+1)
  d0 = max(0.,dn)
  d1 = min(dq,d0)
  q(i,j+1) = (dn - d1)*acosp(j+1)
  dq = dq - d1
  ! South
  ds = q(i,j-1)*cosp(j-1)
  d0 = max(0.,ds)
  d2 = min(dq,d0)
  q(i,j-1) = (ds - d2)*acosp(j-1)
  q(i,j) = (d2 - dq)*acosp(j) + tiny
  endif
  END DO
  end do
  !
  do i=1,imr
  IF(q(i,j1).LT.0.) THEN
  ipy =  1
  dq  = - q(i,j1)*cosp(j1)
  ! North
  dn = q(i,j1+1)*cosp(j1+1)
  d0 = max(0.,dn)
  d1 = min(dq,d0)
  q(i,j1+1) = (dn - d1)*acosp(j1+1)
  q(i,j1) = (d1 - dq)*acosp(j1) + tiny
  endif
  enddo
  !
  j = j2
  do i=1,imr
  IF(q(i,j).LT.0.) THEN
  ipy =  1
  dq  = - q(i,j)*cosp(j)
  ! South
  ds = q(i,j-1)*cosp(j-1)
  d0 = max(0.,ds)
  d2 = min(dq,d0)
  q(i,j-1) = (ds - d2)*acosp(j-1)
  q(i,j) = (d2 - dq)*acosp(j) + tiny
  endif
  enddo
  !
  ! Check Poles.
  if(q(1,1).lt.0.) then
  dq = q(1,1)*cap1/REAL(IMR)*acosp(j1)
  do i=1,imr
  q(i,1) = 0.
  q(i,j1) = q(i,j1) + dq
  if(q(i,j1).lt.0.) ipy = 1
  enddo
  endif
  !
  if(q(1,JNP).lt.0.) then
  dq = q(1,JNP)*cap1/REAL(IMR)*acosp(j2)
  do i=1,imr
  q(i,JNP) = 0.
  q(i,j2) = q(i,j2) + dq
  if(q(i,j2).lt.0.) ipy = 1
  enddo
  endif
  !
  return
end subroutine filns
!
subroutine filew(q,qtmp,IMR,JNP,j1,j2,ipx,tiny)
  implicit none
  integer :: IMR,JNP,j1,j2,ipx
  real :: q(IMR,*),qtmp(JNP,IMR),tiny
  integer :: i,j
  real :: d0,d1,d2
  !
  ipx = 0
  ! Copy & swap direction for vectorization.
  do i=1,imr
  do j=j1,j2
    qtmp(j,i) = q(i,j)
  end do
  end do
  !
  do i=2,imr-1
  do j=j1,j2
  if(qtmp(j,i).lt.0.) then
  ipx =  1
  ! west
  d0 = max(0.,qtmp(j,i-1))
  d1 = min(-qtmp(j,i),d0)
  qtmp(j,i-1) = qtmp(j,i-1) - d1
  qtmp(j,i) = qtmp(j,i) + d1
  ! east
  d0 = max(0.,qtmp(j,i+1))
  d2 = min(-qtmp(j,i),d0)
  qtmp(j,i+1) = qtmp(j,i+1) - d2
  qtmp(j,i) = qtmp(j,i) + d2 + tiny
  endif
  end do
  end do
  !
  i=1
  do j=j1,j2
  if(qtmp(j,i).lt.0.) then
  ipx =  1
  ! west
  d0 = max(0.,qtmp(j,imr))
  d1 = min(-qtmp(j,i),d0)
  qtmp(j,imr) = qtmp(j,imr) - d1
  qtmp(j,i) = qtmp(j,i) + d1
  ! east
  d0 = max(0.,qtmp(j,i+1))
  d2 = min(-qtmp(j,i),d0)
  qtmp(j,i+1) = qtmp(j,i+1) - d2
  !
  qtmp(j,i) = qtmp(j,i) + d2 + tiny
  endif
65  continue
  end do
  i=IMR
  do j=j1,j2
  if(qtmp(j,i).lt.0.) then
  ipx =  1
  ! west
  d0 = max(0.,qtmp(j,i-1))
  d1 = min(-qtmp(j,i),d0)
  qtmp(j,i-1) = qtmp(j,i-1) - d1
  qtmp(j,i) = qtmp(j,i) + d1
  ! east
  d0 = max(0.,qtmp(j,1))
  d2 = min(-qtmp(j,i),d0)
  qtmp(j,1) = qtmp(j,1) - d2
  !
  qtmp(j,i) = qtmp(j,i) + d2 + tiny
  endif
  end do
  !
  if(ipx.ne.0) then
  do j=j1,j2
  do i=1,imr
    q(i,j) = qtmp(j,i)
  end do
  end do
  else
  !
  ! Poles.
  if(q(1,1).lt.0.OR. q(1,JNP).lt.0.) ipx = 1
  endif
  return
end subroutine filew
!
subroutine zflip(q,im,km,nc)
  implicit none
  ! This routine flip the array q (in the vertical).
  integer :: im,km,nc
  real :: q(im,km,nc)
  ! local dynamic array
  real :: qtmp(im,km)
  integer :: IC,k,i
  !
  do IC = 1, nc
  !
  do k=1,km
  do i=1,im
  qtmp(i,k) = q(i,km+1-k,IC)
  end do
  end do
  !
  do i=1,im*km
    q(i,1,IC) = qtmp(i,1)
2000 continue
  end do
  end do
  return
end subroutine zflip