/[lmdze]/trunk/Sources/dyn3d/PPM3d/ppm3d.f
ViewVC logotype

Diff of /trunk/Sources/dyn3d/PPM3d/ppm3d.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/dyn3d/ppm3d.f90 revision 81 by guez, Wed Mar 5 14:38:41 2014 UTC trunk/Sources/dyn3d/PPM3d/ppm3d.f revision 169 by guez, Mon Sep 14 17:13:16 2015 UTC
# Line 1  Line 1 
1    module ppm3d_m
2    
3  ! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/ppm3d.F,v 1.1.1.1 2004/05/19    implicit none
 ! 12:53:07 lmdzadmin Exp $  
4    
5    contains
6    
7  ! From lin@explorer.gsfc.nasa.gov Wed Apr 15 17:44:44 1998    SUBROUTINE ppm3d(igd, q, ps1, ps2, u, v, w, ndt, iord, jord, kord, nc, imr, &
8  ! Date: Wed, 15 Apr 1998 11:37:03 -0400         jnp, j1, nlay, ap, bp, pt, ae, fill, umax)
 ! 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 number of constituents  
   ! IMR: first dimension (E-W); number of Grid intervals in E-W is IMR  
   ! JNP: 2nd dimension (N-S); number of Grid intervals in N-S is JNP-1  
   ! NLAY: 3rd dimension (number of layers); vertical index increases from 1  
   ! at  
   ! the model top to NLAY near the surface (see fig. below).  
   ! It is assumed that 6 <= NLAY <= JNP (for dynamic memory allocation)  
   
   ! PS1(IMR,JNP): surface pressure at current time (t)  
   ! PS2(IMR,JNP): surface pressure at mid-time-level (t+NDT/2)  
   ! PS2 is replaced by the predicted PS (at t+NDT) on output.  
   ! Note: surface pressure can have any unit or can be multiplied by any  
   ! const.  
   
   ! The pressure at layer edges are defined as follows:  
   
   ! p(i,j,k) = AP(k)*PT  +  BP(k)*PS(i,j)          (1)  
   
   ! Where PT is a constant having the same unit as PS.  
   ! AP and BP are unitless constants given at layer edges  
   ! defining the vertical coordinate.  
   ! BP(1) = 0., BP(NLAY+1) = 1.  
   ! The pressure at the model top is PTOP = AP(1)*PT  
   
   ! For pure sigma system set AP(k) = 1 for all k, PT = PTOP,  
   ! BP(k) = sige(k) (sigma at edges), PS = Psfc - PTOP.  
   
   ! Note: the sigma-P coordinate is a subset of Eq. 1, which in turn  
   ! is a subset of the following even more general sigma-P-thelta coord.  
   ! currently under development.  
   ! p(i,j,k) = (AP(k)*PT + BP(k)*PS(i,j))/(D(k)-C(k)*TE**(-1/kapa))  
   
   ! /////////////////////////////////  
   ! / \ ------------- PTOP --------------  AP(1), BP(1)  
   ! |  
   ! delp(1)    |  ........... Q(i,j,1) ............  
   ! |  
   ! W(1)    \ / ---------------------------------  AP(2), BP(2)  
   
   
   
   ! W(k-1)   / \ ---------------------------------  AP(k), BP(k)  
   ! |  
   ! delp(K)    |  ........... Q(i,j,k) ............  
   ! |  
   ! W(k)    \ / ---------------------------------  AP(k+1), BP(k+1)  
   
   
   
   ! / \ ---------------------------------  AP(NLAY), BP(NLAY)  
   ! |  
   ! delp(NLAY)   |  ........... Q(i,j,NLAY) .........  
   ! |  
   ! W(NLAY)=0  \ / ------------- surface ----------- AP(NLAY+1), BP(NLAY+1)  
   ! //////////////////////////////////  
   
   ! U(IMR,JNP,NLAY) & V(IMR,JNP,NLAY):winds (m/s) at mid-time-level (t+NDT/2)  
   ! U and V may need to be polar filtered in advance in some cases.  
   
   ! IGD:      grid type on which winds are defined.  
   ! IGD = 0:  A-Grid  [all variables defined at the same point from south  
   ! pole (j=1) to north pole (j=JNP) ]  
   
   ! IGD = 1  GEOS-GCM C-Grid  
   ! [North]  
   
   ! V(i,j)  
   ! |  
   ! |  
   ! |  
   ! U(i-1,j)---Q(i,j)---U(i,j) [EAST]  
   ! |  
   ! |  
   ! |  
   ! V(i,j-1)  
   
   ! U(i,  1) is defined at South Pole.  
   ! V(i,  1) is half grid north of the South Pole.  
   ! V(i,JMR) is half grid south of the North Pole.  
   
   ! V must be defined at j=1 and j=JMR if IGD=1  
   ! V at JNP need not be given.  
   
   ! NDT: time step in seconds (need not be constant during the course of  
   ! the integration). Suggested value: 30 min. for 4x5, 15 min. for 2x2.5  
   ! (Lat-Lon) resolution. Smaller values are recommanded if the model  
   ! has a well-resolved stratosphere.  
   
   ! J1 defines the size of the polar cap:  
   ! South polar cap edge is located at -90 + (j1-1.5)*180/(JNP-1) deg.  
   ! North polar cap edge is located at  90 - (j1-1.5)*180/(JNP-1) deg.  
   ! There are currently only two choices (j1=2 or 3).  
   ! IMR must be an even integer if j1 = 2. Recommended value: J1=3.  
   
   ! IORD, JORD, and KORD are integers controlling various options in E-W,  
   ! N-S,  
   ! and vertical transport, respectively. Recommended values for positive  
   ! definite scalars: IORD=JORD=3, KORD=5. Use KORD=3 for non-  
   ! positive definite scalars or when linear correlation between constituents  
   ! is to be maintained.  
   
   ! _ORD=  
   ! 1: 1st order upstream scheme (too diffusive, not a useful option; it  
   ! can be used for debugging purposes; this is THE only known "linear"  
   ! monotonic advection scheme.).  
   ! 2: 2nd order van Leer (full monotonicity constraint;  
   ! see Lin et al 1994, MWR)  
   ! 3: monotonic PPM* (slightly improved PPM of Collela & Woodward 1984)  
   ! 4: semi-monotonic PPM (same as 3, but overshoots are allowed)  
   ! 5: positive-definite PPM (constraint on the subgrid distribution is  
   ! only strong enough to prevent generation of negative values;  
   ! both overshoots & undershoots are possible).  
   ! 6: un-constrained PPM (nearly diffusion free; slightly faster but  
   ! positivity not quaranteed. Use this option only when the fields  
   ! and winds are very smooth).  
   
   ! *PPM: Piece-wise Parabolic Method  
   
   ! Note that KORD <=2 options are no longer supported. DO not use option 4  
   ! or 5.  
   ! for non-positive definite scalars (such as Ertel Potential Vorticity).  
   
   ! The implicit numerical diffusion decreases as _ORD increases.  
   ! The last two options (ORDER=5, 6) should only be used when there is  
   ! significant explicit diffusion (such as a turbulence parameterization).  
   ! You  
   ! might get dispersive results otherwise.  
   ! No filter of any kind is applied to the constituent fields here.  
   
   ! AE: Radius of the sphere (meters).  
   ! Recommended value for the planet earth: 6.371E6  
   
   ! fill(logical):   flag to do filling for negatives (see note below).  
   
   ! Umax: Estimate (upper limit) of the maximum U-wind speed (m/s).  
   ! (220 m/s is a good value for troposphere model; 280 m/s otherwise)  
   
   ! =============  
   ! Output  
   ! =============  
   
   ! Q: mixing ratios at future time (t+NDT) (original values are  
   ! over-written)  
   ! W(NLAY): large-scale vertical mass flux as diagnosed from the hydrostatic  
   ! relationship. W will have the same unit as PS1 and PS2 (eg, mb).  
   ! W must be divided by NDT to get the correct mass-flux unit.  
   ! The vertical Courant number C = W/delp_UPWIND, where delp_UPWIND  
   ! is the pressure thickness in the "upwind" direction. For example,  
   ! C(k) = W(k)/delp(k)   if W(k) > 0;  
   ! C(k) = W(k)/delp(k+1) if W(k) < 0.  
   ! ( W > 0 is downward, ie, toward surface)  
   ! PS2: predicted PS at t+NDT (original values are over-written)  
   
   ! ********************************************************************  
   ! NOTES:  
   ! This forward-in-time upstream-biased transport scheme reduces to  
   ! the 2nd order center-in-time center-in-space mass continuity eqn.  
   ! if Q = 1 (constant fields will remain constant). This also ensures  
   ! that the computed vertical velocity to be identical to GEOS-1 GCM  
   ! for on-line transport.  
   
   ! A larger polar cap is used if j1=3 (recommended for C-Grid winds or when  
   ! winds are noisy near poles).  
   
   ! Flux-Form Semi-Lagrangian transport in the East-West direction is used  
   ! when and where Courant number is greater than one.  
   
   ! The user needs to change the parameter Jmax or Kmax if the resolution  
   ! is greater than 0.5 deg in N-S or 150 layers in the vertical direction.  
   ! (this TransPort Core is otherwise resolution independent and can be used  
   ! as a library routine).  
   
   ! PPM is 4th order accurate when grid spacing is uniform (x & y); 3rd  
   ! order accurate for non-uniform grid (vertical sigma coord.).  
   
   ! Time step is limitted only by transport in the meridional direction.  
   ! (the FFSL scheme is not implemented in the meridional direction).  
   
   ! Since only 1-D limiters are applied, negative values could  
   ! potentially be generated when large time step is used and when the  
   ! initial fields contain discontinuities.  
   ! This does not necessarily imply the integration is unstable.  
   ! These negatives are typically very small. A filling algorithm is  
   ! activated if the user set "fill" to be true.  
   
   ! The van Leer scheme used here is nearly as accurate as the original PPM  
   ! due to the use of a 4th order accurate reference slope. The PPM imple-  
   ! mented here is an improvement over the original and is also based on  
   ! the 4th order reference slope.  
   
   ! ****6***0*********0*********0*********0*********0*********0**********72  
   
   ! User modifiable parameters  
   
   PARAMETER (jmax=361, kmax=150)  
   
   ! ****6***0*********0*********0*********0*********0*********0**********72  
   
   ! Input-Output arrays  
   
   
   REAL q(imr, jnp, nlay, nc), ps1(imr, jnp), ps2(imr, jnp), &  
     u(imr, jnp, nlay), v(imr, jnp, nlay), ap(nlay+1), bp(nlay+1), &  
     w(imr, jnp, nlay), ndt, val(nlay), umax  
   INTEGER igd, iord, jord, kord, nc, imr, jnp, j1, nlay, ae  
   INTEGER imrd2  
   REAL pt  
   LOGICAL cross, fill, dum  
   
   ! Local dynamic arrays  
   
   REAL crx(imr, jnp), cry(imr, jnp), xmass(imr, jnp), ymass(imr, jnp), &  
     fx1(imr+1), dpi(imr, jnp, nlay), delp1(imr, jnp, nlay), &  
     wk1(imr, jnp, nlay), pu(imr, jnp), pv(imr, jnp), dc2(imr, jnp), &  
     delp2(imr, jnp, nlay), dq(imr, jnp, nlay, nc), va(imr, jnp), &  
     ua(imr, jnp), qtmp(-imr:2*imr)  
   
   ! Local static  arrays  
   
   REAL dtdx(jmax), dtdx5(jmax), acosp(jmax), cosp(jmax), cose(jmax), &  
     dap(kmax), dbk(kmax)  
   DATA ndt0, nstep/0, 0/  
   DATA cross/.TRUE./  
   SAVE dtdy, dtdy5, rcap, js0, jn0, iml, dtdx, dtdx5, acosp, cosp, cose, dap, &  
     dbk  
   
   
   jmr = jnp - 1  
   imjm = imr*jnp  
   j2 = jnp - j1 + 1  
   nstep = nstep + 1  
   
   ! *********** Initialization **********************  
   IF (nstep==1) THEN  
   
     WRITE (6, *) '------------------------------------ '  
     WRITE (6, *) 'NASA/GSFC Transport Core Version 4.5'  
     WRITE (6, *) '------------------------------------ '  
   
     WRITE (6, *) 'IMR=', imr, ' JNP=', jnp, ' NLAY=', nlay, ' j1=', j1  
     WRITE (6, *) 'NC=', nc, iord, jord, kord, ndt  
   
     ! controles sur les parametres  
     IF (nlay<6) THEN  
       WRITE (6, *) 'NLAY must be >= 6'  
       STOP  
     END IF  
     IF (jnp<nlay) THEN  
       WRITE (6, *) 'JNP must be >= NLAY'  
       STOP  
     END IF  
     imrd2 = mod(imr, 2)  
     IF (j1==2 .AND. imrd2/=0) THEN  
       WRITE (6, *) 'if j1=2 IMR must be an even integer'  
       STOP  
     END IF  
   
   
     IF (jmax<jnp .OR. kmax<nlay) THEN  
       WRITE (6, *) 'Jmax or Kmax is too small'  
       STOP  
     END IF  
   
     DO k = 1, nlay  
       dap(k) = (ap(k+1)-ap(k))*pt  
       dbk(k) = bp(k+1) - bp(k)  
     END DO  
9    
10      pi = 4.*atan(1.)      ! INPUT:
11      dl = 2.*pi/float(imr)      ! =============
12      dp = pi/float(jmr)  
13        ! Q(IMR,JNP,NLAY,NC): mixing ratios at current time (t)
14      IF (igd==0) THEN      ! NC: total number of constituents
15        ! Compute analytic cosine at cell edges      ! IMR: first dimension (E-W); number of Grid intervals in E-W is IMR
16        CALL cosa(cosp, cose, jnp, pi, dp)      ! JNP: 2nd dimension (N-S); number of Grid intervals in N-S is JNP-1
17      ELSE      ! NLAY: 3rd dimension (number of layers); vertical index increases from 1
18        ! Define cosine consistent with GEOS-GCM (using dycore2.0 or later)      ! at
19        CALL cosc(cosp, cose, jnp, pi, dp)      ! the model top to NLAY near the surface (see fig. below).
20      END IF      ! It is assumed that 6 <= NLAY <= JNP (for dynamic memory allocation)
21    
22      DO j = 2, jmr      ! PS1(IMR,JNP): surface pressure at current time (t)
23        acosp(j) = 1./cosp(j)      ! PS2(IMR,JNP): surface pressure at mid-time-level (t+NDT/2)
24      END DO      ! 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      ! Inverse of the Scaled polar cap area.      ! const.
27    
28      rcap = dp/(imr*(1.-cos((j1-1.5)*dp)))      ! The pressure at layer edges are defined as follows:
29      acosp(1) = rcap  
30      acosp(jnp) = rcap      ! p(i,j,k) = AP(k)*PT  +  BP(k)*PS(i,j)          (1)
31    END IF  
32        ! Where PT is a constant having the same unit as PS.
33    IF (ndt0/=ndt) THEN      ! AP and BP are unitless constants given at layer edges
34      dt = ndt      ! defining the vertical coordinate.
35      ndt0 = ndt      ! BP(1) = 0., BP(NLAY+1) = 1.
36        ! The pressure at the model top is PTOP = AP(1)*PT
37      IF (umax<180.) THEN  
38        WRITE (6, *) 'Umax may be too small!'      ! For pure sigma system set AP(k) = 1 for all k, PT = PTOP,
39      END IF      ! BP(k) = sige(k) (sigma at edges), PS = Psfc - PTOP.
40      cr1 = abs(umax*dt)/(dl*ae)  
41      maxdt = dp*ae/abs(umax) + 0.5      ! Note: the sigma-P coordinate is a subset of Eq. 1, which in turn
42      WRITE (6, *) 'Largest time step for max(V)=', umax, ' is ', maxdt      ! is a subset of the following even more general sigma-P-thelta coord.
43      IF (maxdt<abs(ndt)) THEN      ! currently under development.
44        WRITE (6, *) 'Warning!!! NDT maybe too large!'      ! p(i,j,k) = (AP(k)*PT + BP(k)*PS(i,j))/(D(k)-C(k)*TE**(-1/kapa))
45      END IF  
46        ! Cf. ppm3d.txt.
47      IF (cr1>=0.95) THEN  
48        js0 = 0      ! U(IMR,JNP,NLAY) & V(IMR,JNP,NLAY):winds (m/s) at mid-time-level (t+NDT/2)
49        jn0 = 0      ! U and V may need to be polar filtered in advance in some cases.
50        iml = imr - 2  
51        ztc = 0.      ! IGD:      grid type on which winds are defined.
52      ELSE      ! IGD = 0:  A-Grid  [all variables defined at the same point from south
53        ztc = acos(cr1)*(180./pi)      ! pole (j=1) to north pole (j=JNP) ]
54    
55        js0 = float(jmr)*(90.-ztc)/180. + 2      ! IGD = 1  GEOS-GCM C-Grid
56        js0 = max(js0, j1+1)      ! Cf. ppm3d.txt.
57        iml = min(6*js0/(j1-1)+2, 4*imr/5)  
58        jn0 = jnp - js0 + 1      ! U(i,  1) is defined at South Pole.
59      END IF      ! 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        ! ****6***0*********0*********0*********0*********0*********0**********72
173    
174      DO j = 2, jmr      ! User modifiable parameters
       dtdx(j) = dt/(dl*ae*cosp(j))  
175    
176        dtdx5(j) = 0.5*dtdx(j)      integer, PARAMETER:: jmax=361, kmax=150
     END DO  
177    
178        ! ****6***0*********0*********0*********0*********0*********0**********72
179    
180      dtdy = dt/(ae*dp)      ! Input-Output arrays
     dtdy5 = 0.5*dtdy  
181    
182    END IF      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(umax*dt)/(dl*ae)
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(6*js0/(j1-1)+2, 4*imr/5)
300              jn0 = jnp - js0 + 1
301           END IF
302    
303    
304           DO j = 2, jmr
305              dtdx(j) = dt/(dl*ae*cosp(j))
306    
307              dtdx5(j) = 0.5*dtdx(j)
308           END DO
309    
   ! *********** End Initialization **********************  
310    
311    ! delp = pressure thickness: the psudo-density in a hydrostatic system.         dtdy = dt/(ae*dp)
312    DO k = 1, nlay         dtdy5 = 0.5*dtdy
     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  
313    
314        END IF
315    
316    IF (j1/=2) THEN      ! *********** End Initialization **********************
     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  
317    
318    ! Compute "tracer density"      ! delp = pressure thickness: the psudo-density in a hydrostatic system.
   DO ic = 1, nc  
319      DO k = 1, nlay      DO k = 1, nlay
320        DO j = 1, jnp         DO j = 1, jnp
321          DO i = 1, imr            DO i = 1, imr
322            dq(i, j, k, ic) = q(i, j, k, ic)*delp1(i, j, k)               delp1(i, j, k) = dap(k) + dbk(k)*ps1(i, j)
323          END DO               delp2(i, j, k) = dap(k) + dbk(k)*ps2(i, j)
324        END DO            END DO
325      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  
326      END DO      END DO
327    
 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  
328    
329      IF (j1/=2) THEN ! Enlarged polar cap.      IF (j1/=2) THEN
330        DO i = 1, imr         DO ic = 1, nc
331          dpi(i, 2, k) = 0.            DO l = 1, nlay
332          dpi(i, jmr, k) = 0.               DO i = 1, imr
333        END DO                  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      END IF
339    
340      ! ******* Compute horizontal mass fluxes ************      ! Compute "tracer density"
341        DO ic = 1, nc
342      ! N-S component         DO k = 1, nlay
343      DO j = j1, j2 + 1            DO j = 1, jnp
344        d5 = 0.5*cose(j)               DO i = 1, imr
345        DO i = 1, imr                  dq(i, j, k, ic) = q(i, j, k, ic)*delp1(i, j, k)
346          ymass(i, j) = cry(i, j)*d5*(delp2(i,j,k)+delp2(i,j-1,k))               END DO
347        END DO            END DO
348      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  
349      END DO      END DO
350    
351      DO j = j1, j2      DO k = 1, nlay
       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  
352    
353      DO i = 1, imr*(jmr-1)         IF (igd==0) THEN
354        va(i, 2) = 0.5*(cry(i,2)+cry(i,3))            ! Convert winds on A-Grid to Courant number on C-Grid.
355      END DO            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    
     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  
364    
365      ! ****6***0*********0*********0*********0*********0*********0**********72            DO j = j1, j2
366      DO ic = 1, nc               crx(1, j) = dtdx(j)*u(imr, j, k)
367              END DO
368    
369        DO i = 1, imjm            DO i = 1, imr*jmr
370          wk1(i, 1, 1) = 0.               cry(i, 2) = dtdy*v(i, 1, k)
371          wk1(i, 1, 2) = 0.            END DO
372        END DO         END IF
   
       ! E-W advective cross term  
       DO j = j1, j2  
         IF (j>js .AND. j<jn) GO TO 250  
   
         DO i = 1, imr  
           qtmp(i) = q(i, j, k, ic)  
         END DO  
   
         DO i = -iml, 0  
           qtmp(i) = q(imr+i, j, k, ic)  
           qtmp(imr+1-i) = q(1-i, j, k, ic)  
         END DO  
   
         DO i = 1, imr  
           iu = ua(i, j)  
           ru = ua(i, j) - iu  
           iiu = i - iu  
           IF (ua(i,j)>=0.) THEN  
             wk1(i, j, 1) = qtmp(iiu) + ru*(qtmp(iiu-1)-qtmp(iiu))  
           ELSE  
             wk1(i, j, 1) = qtmp(iiu) + ru*(qtmp(iiu)-qtmp(iiu+1))  
           END IF  
           wk1(i, j, 1) = wk1(i, j, 1) - qtmp(i)  
         END DO  
 250   END DO  
373    
374        IF (jn/=0) THEN         ! Determine JS and JN
375          DO j = js + 1, jn - 1         js = j1
376           jn = j2
377    
378           DO j = js0, j1 + 1, -1
379            DO i = 1, imr            DO i = 1, imr
380              qtmp(i) = q(i, j, k, ic)               IF (abs(crx(i,j))>1.) THEN
381                    js = j
382                    GO TO 2222
383                 END IF
384            END DO            END DO
385           END DO
386    
387            qtmp(0) = q(imr, j, k, ic)  2222   CONTINUE
388            qtmp(imr+1) = q(1, j, k, ic)         DO j = jn0, j2 - 1
   
389            DO i = 1, imr            DO i = 1, imr
390              iu = i - ua(i, j)               IF (abs(crx(i,j))>1.) THEN
391              wk1(i, j, 1) = ua(i, j)*(qtmp(iu)-qtmp(iu+1))                  jn = j
392                    GO TO 2233
393                 END IF
394            END DO            END DO
395          END DO         END DO
396        END IF  2233   CONTINUE
397        ! ****6***0*********0*********0*********0*********0*********0**********72  
398        ! Contribution from the N-S advection         IF (j1/=2) THEN ! Enlarged polar cap.
       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  
399            DO i = 1, imr            DO i = 1, imr
400              q(i, j, k, ic) = q(i, j, k, ic) + dc2(i, j) + pv(i, j)               dpi(i, 2, k) = 0.
401                 dpi(i, jmr, k) = 0.
402            END DO            END DO
403          END DO         END IF
       END IF  
   
       CALL xtp(imr, jnp, iml, j1, j2, jn, js, pu, dq(1,1,k,ic), wk1(1,1,2), &  
         crx, fx1, xmass, iord)  
   
       CALL ytp(imr, jnp, j1, j2, acosp, rcap, dq(1,1,k,ic), wk1(1,1,1), cry, &  
         dc2, ymass, wk1(1,1,3), wk1(1,1,4), wk1(1,1,5), wk1(1,1,6), jord)  
   
     END DO  
   END DO  
   
   ! ******* Compute vertical mass flux (same unit as PS) ***********  
   
   ! 1st step: compute total column mass CONVERGENCE.  
   
   DO j = 1, jnp  
     DO i = 1, imr  
       cry(i, j) = dpi(i, j, 1)  
     END DO  
   END DO  
   
   DO k = 2, nlay  
     DO j = 1, jnp  
       DO i = 1, imr  
         cry(i, j) = cry(i, j) + dpi(i, j, k)  
       END DO  
     END DO  
   END DO  
   
   DO j = 1, jnp  
     DO i = 1, imr  
   
       ! 2nd step: compute PS2 (PS at n+1) using the hydrostatic assumption.  
       ! Changes (increases) to surface pressure = total column mass  
       ! convergence  
   
       ps2(i, j) = ps1(i, j) + cry(i, j)  
   
       ! 3rd step: compute vertical mass flux from mass conservation  
       ! principle.  
   
       w(i, j, 1) = dpi(i, j, 1) - dbk(1)*cry(i, j)  
       w(i, j, nlay) = 0.  
     END DO  
   END DO  
   
   DO k = 2, nlay - 1  
     DO j = 1, jnp  
       DO i = 1, imr  
         w(i, j, k) = w(i, j, k-1) + dpi(i, j, k) - dbk(k)*cry(i, j)  
       END DO  
     END DO  
   END DO  
   
   DO k = 1, nlay  
     DO j = 1, jnp  
       DO i = 1, imr  
         delp2(i, j, k) = dap(k) + dbk(k)*ps2(i, j)  
       END DO  
     END DO  
   END DO  
   
   krd = max(3, kord)  
   DO ic = 1, nc  
   
     ! ****6***0*********0*********0*********0*********0*********0**********72  
   
     CALL fzppm(imr, jnp, nlay, j1, dq(1,1,1,ic), w, q(1,1,1,ic), wk1, dpi, &  
       dc2, crx, cry, pu, pv, xmass, ymass, delp1, krd)  
   
   
     IF (fill) CALL qckxyz(dq(1,1,1,ic), dc2, imr, jnp, nlay, j1, j2, cosp, &  
       acosp, .FALSE., ic, nstep)  
   
     ! Recover tracer mixing ratio from "density" using predicted  
     ! "air density" (pressure thickness) at time-level n+1  
   
     DO k = 1, nlay  
       DO j = 1, jnp  
         DO i = 1, imr  
           q(i, j, k, ic) = dq(i, j, k, ic)/delp2(i, j, k)  
         END DO  
       END DO  
     END DO  
   
     IF (j1/=2) THEN  
       DO k = 1, nlay  
         DO i = 1, imr  
           ! j=1 c'est le pôle Sud, j=JNP c'est le pôle Nord  
           q(i, 2, k, ic) = q(i, 1, k, ic)  
           q(i, jmr, k, ic) = q(i, jmp, k, ic)  
         END DO  
       END DO  
     END IF  
   END DO  
   
   IF (j1/=2) THEN  
     DO k = 1, nlay  
       DO i = 1, imr  
         w(i, 2, k) = w(i, 1, k)  
         w(i, jmr, k) = w(i, jnp, k)  
       END DO  
     END DO  
   END IF  
   
   RETURN  
 END SUBROUTINE ppm3d  
   
 ! ****6***0*********0*********0*********0*********0*********0**********72  
 SUBROUTINE fzppm(imr, jnp, nlay, j1, dq, wz, p, dc, dqdt, ar, al, a6, flux, &  
     wk1, wk2, wz2, delp, kord)  
   PARAMETER (kmax=150)  
   PARAMETER (r23=2./3., r3=1./3.)  
   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, *)  
   
   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==2 .OR. j==jmr) .AND. j1/=2) GO TO 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  
       END DO  
     END DO  
   
     ! ****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)<=0.) THEN  
         al(i, 1) = 0.  
         flux(i, 1) = 0.  
       ELSE  
         flux(i, 1) = wk1(i, 1) - al(i, 1)  
       END IF  
     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)<=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))  
       END DO  
     END DO  
   
     DO i = 1, imr*nlaym1  
       ar(i, 1) = al(i, 2)  
     END DO  
   
     DO i = 1, imr*nlay  
       a6(i, 1) = 3.*(wk1(i,1)+wk1(i,1)-(al(i,1)+ar(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<=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)>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  
         cp = wz2(i, 1)/wk2(i, 2)  
         flux(i, 2) = al(i, 2) + 0.5*cp*(al(i,2)-ar(i,2)-a6(i,2)*(1.+r23*cp))  
       END IF  
     END DO  
   
     DO i = 1, imr*nlaym1  
       flux(i, 2) = wz2(i, 1)*flux(i, 2)  
     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 END DO  
   RETURN  
 END SUBROUTINE fzppm  
   
 SUBROUTINE xtp(imr, jnp, iml, j1, j2, jn, js, pu, dq, q, uc, fx1, xmass, &  
     iord)  
   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), isave(imr)  
   
   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)  
     END DO  
   
     IF (j>=jn .OR. j<=js) GO TO 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==1 .OR. j==j1 .OR. j==j2) THEN  
       DO i = 1, imr  
         iu = float(i) - uc(i, j)  
         fx1(i) = qtmp(iu)  
       END DO  
     ELSE  
       CALL xmist(imr, iml, qtmp, dc)  
       dc(0) = dc(imr)  
   
       IF (iord==2 .OR. j<=j1vl .OR. j>=j2vl) THEN  
         DO i = 1, imr  
           iu = float(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)  
       END IF  
   
     END IF  
   
     DO i = 1, imr  
       fx1(i) = fx1(i)*xmass(i, j)  
     END DO  
   
     GO TO 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)  
     END DO  
   
     IF (iord==1 .OR. j==j1 .OR. j==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)  
       END DO  
   
       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  
     END IF  
   
     DO i = 1, imr  
       IF (uc(i,j)>1.) THEN  
         ! DIR$ NOVECTOR  
         DO ist = isave(i), i - 1  
           fx1(i) = fx1(i) + qtmp(ist)  
         END DO  
       ELSE IF (uc(i,j)<-1.) THEN  
         DO ist = i, isave(i) - 1  
           fx1(i) = fx1(i) - qtmp(ist)  
         END DO  
         ! DIR$ VECTOR  
       END IF  
     END DO  
     DO i = 1, imr  
       fx1(i) = pu(i, j)*fx1(i)  
     END DO  
   
     ! ***************************************  
   
 1309 fx1(imp) = fx1(1)  
     DO i = 1, imr  
       dq(i, j) = dq(i, j) + fx1(i) - fx1(i+1)  
     END DO  
404    
405      ! ***************************************         ! ******* Compute horizontal mass fluxes ************
   
   END DO  
   RETURN  
 END SUBROUTINE xtp  
   
 SUBROUTINE fxppm(imr, iml, ut, p, dc, flux, iord)  
   PARAMETER (r3=1./3., r23=2./3.)  
   DIMENSION ut(*), flux(*), p(-iml:imr+iml+1), dc(-iml:imr+iml+1)  
   DIMENSION ar(0:imr), al(0:imr), a6(0:imr)  
   INTEGER lmt  
   ! 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  
   
   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<=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)>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)))  
     END IF  
   END DO  
   RETURN  
 END SUBROUTINE fxppm  
   
 SUBROUTINE xmist(imr, iml, p, dc)  
   PARAMETER (r24=1./24.)  
   DIMENSION p(-iml:imr+1+iml), dc(-iml:imr+1+iml)  
   
   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)  
   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)  
   
   jmr = jnp - 1  
   len = imr*(j2-j1+2)  
   
   IF (jord==1) THEN  
     DO i = 1, len  
       jt = float(j1) - vc(i, j1)  
       fx(i, j1) = p(i, jt)  
     END DO  
   ELSE  
   
     CALL ymist(imr, jnp, j1, p, dc2, 4)  
   
     IF (jord<=0 .OR. jord>=3) THEN  
   
       CALL fyppm(vc, p, dc2, fx, imr, jnp, j1, j2, a6, ar, al, jord)  
   
     ELSE  
       DO i = 1, len  
         jt = float(j1) - vc(i, j1)  
         fx(i, j1) = p(i, jt) + (sign(1.,vc(i,j1))-vc(i,j1))*dc2(i, jt)  
       END DO  
     END IF  
   END IF  
406    
407    DO i = 1, len         ! N-S component
408      fx(i, j1) = fx(i, j1)*ymass(i, j1)         DO j = j1, j2 + 1
409    END DO            d5 = 0.5*cose(j)
410              DO i = 1, imr
411    DO j = j1, j2               ymass(i, j) = cry(i, j)*d5*(delp2(i,j,k)+delp2(i,j-1,k))
412      DO i = 1, imr            END DO
413        dq(i, j) = dq(i, j) + (fx(i,j)-fx(i,j+1))*acosp(j)         END DO
     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)  
   END DO  
   
   sum1 = dq(1, 1) - sum1*rcap  
   sum2 = dq(1, jnp) + sum2*rcap  
   DO i = 1, imr  
     dq(i, 1) = sum1  
     dq(i, jnp) = sum2  
   END DO  
   
   IF (j1/=2) THEN  
     DO i = 1, imr  
       dq(i, 2) = sum1  
       dq(i, jmr) = sum2  
     END DO  
   END IF  
   
   RETURN  
 END SUBROUTINE ytp  
414    
415  SUBROUTINE ymist(imr, jnp, j1, p, dc, id)         DO j = j1, j2
416    PARAMETER (r24=1./24.)            DO i = 1, imr
417    DIMENSION p(imr, jnp), dc(imr, jnp)               dpi(i, j, k) = (ymass(i,j)-ymass(i,j+1))*acosp(j)
418              END DO
419    imh = imr/2         END DO
   jmr = jnp - 1  
   ijm3 = imr*(jmr-3)  
   
   IF (id==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  
420    
421      DO i = 1, ijm3         ! Poles
422        tmp = (8.*(p(i,4)-p(i,2))+p(i,1)-p(i,5))*r24         sum1 = ymass(imr, j1)
423        pmax = max(p(i,2), p(i,3), p(i,4)) - p(i, 3)         sum2 = ymass(imr, j2+1)
424        pmin = p(i, 3) - min(p(i,2), p(i,3), p(i,4))         DO i = 1, imr - 1
425        dc(i, 3) = sign(min(abs(tmp),pmin,pmax), tmp)            sum1 = sum1 + ymass(i, j1)
426      END DO            sum2 = sum2 + ymass(i, j2+1)
427    END IF         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    IF (j1/=2) THEN         DO j = j1, j2
445      DO i = 1, imr            pu(1, j) = 0.5*(delp2(1,j,k)+delp2(imr,j,k))
446        dc(i, 1) = 0.         END DO
       dc(i, jnp) = 0.  
     END DO  
   ELSE  
     ! Determine slopes in polar caps for scalars!  
447    
448      DO i = 1, imh         DO j = j1, j2
449        ! South            DO i = 1, imr
450        tmp = 0.25*(p(i,2)-p(i+imh,2))               xmass(i, j) = pu(i, j)*crx(i, j)
451        pmax = max(p(i,2), p(i,1), p(i+imh,2)) - p(i, 1)            END DO
452        pmin = p(i, 1) - min(p(i,2), p(i,1), p(i+imh,2))         END DO
       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  
453    
454      DO i = imh + 1, imr         DO j = j1, j2
455        dc(i, 1) = -dc(i-imh, 1)            DO i = 1, imr - 1
456        dc(i, jnp) = -dc(i-imh, jnp)               dpi(i, j, k) = dpi(i, j, k) + xmass(i, j) - xmass(i+1, j)
457      END DO            END DO
458    END IF         END DO
   RETURN  
 END SUBROUTINE ymist  
   
 SUBROUTINE fyppm(vc, p, dc, flux, imr, jnp, j1, j2, a6, ar, al, jord)  
   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  
   ! 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)  
   END DO  
   
   ! 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<=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)>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)))  
     END IF  
   END DO  
   RETURN  
 END SUBROUTINE fyppm  
   
 SUBROUTINE yadv(imr, jnp, j1, j2, p, va, ady, wk, iad)  
   REAL p(imr, jnp), ady(imr, jnp), va(imr, jnp)  
   REAL wk(imr, -1:jnp+2)  
   
   jmr = jnp - 1  
   imh = imr/2  
   DO j = 1, jnp  
     DO i = 1, imr  
       wk(i, j) = p(i, j)  
     END DO  
   END DO  
   ! 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)  
   END DO  
   
   IF (iad==2) THEN  
     DO j = j1 - 1, j2 + 1  
       DO i = 1, imr  
         jp = nint(va(i,j))  
         rv = jp - va(i, j)  
         jp = j - 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))  
         ady(i, j) = wk(i, jp) + rv*(a1*rv+b1) - wk(i, j)  
       END DO  
     END DO  
459    
460    ELSE IF (iad==1) THEN         DO j = j1, j2
461      DO j = j1 - 1, j2 + 1            dpi(imr, j, k) = dpi(imr, j, k) + xmass(imr, j) - xmass(1, j)
462        DO i = 1, imr         END DO
463          jp = float(j) - va(i, j)  
464          ady(i, j) = va(i, j)*(wk(i,jp)-wk(i,jp+1))         DO j = j1, j2
465        END DO            DO i = 1, imr - 1
466      END DO               ua(i, j) = 0.5*(crx(i,j)+crx(i+1,j))
467    END IF            END DO
468           END DO
469    
470    IF (j1/=2) THEN         DO j = j1, j2
471      sum1 = 0.            ua(imr, j) = 0.5*(crx(imr,j)+crx(1,j))
472      sum2 = 0.         END DO
473      DO i = 1, imr         ! cccccccccccccccccccccccccccccccccccccccccccccccccccccc
474        sum1 = sum1 + ady(i, 2)         ! Rajouts pour LMDZ.3.3
475        sum2 = sum2 + ady(i, jmr)         ! cccccccccccccccccccccccccccccccccccccccccccccccccccccc
476      END DO         DO i = 1, imr
477      sum1 = sum1/imr            DO j = 1, jnp
478      sum2 = sum2/imr               va(i, j) = 0.
479              END DO
480           END DO
481    
482      DO i = 1, imr         DO i = 1, imr*(jmr-1)
483        ady(i, 2) = sum1            va(i, 2) = 0.5*(cry(i,2)+cry(i,3))
484        ady(i, jmr) = sum2         END DO
485        ady(i, 1) = sum1  
486        ady(i, jnp) = sum2         IF (j1==2) THEN
487      END DO            imh = imr/2
488    ELSE            DO i = 1, imh
489      ! Poles:               va(i, 1) = 0.5*(cry(i,2)-cry(i+imh,2))
490      sum1 = 0.               va(i+imh, 1) = -va(i, 1)
491      sum2 = 0.               va(i, jnp) = 0.5*(cry(i,jnp)-cry(i+imh,jmr))
492      DO i = 1, imr               va(i+imh, jnp) = -va(i, jnp)
493        sum1 = sum1 + ady(i, 1)            END DO
494        sum2 = sum2 + ady(i, jnp)            va(imr, 1) = va(1, 1)
495      END DO            va(imr, jnp) = va(1, jnp)
496      sum1 = sum1/imr         END IF
497      sum2 = sum2/imr  
498           ! ****6***0*********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      DO i = 1, imr            ! E-W advective cross term
507        ady(i, 1) = sum1            DO j = j1, j2
508        ady(i, jnp) = sum2               IF (j>js .AND. j<jn) cycle
509      END DO  
510    END IF               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    RETURN            IF (jn/=0) THEN
533  END SUBROUTINE yadv               DO j = js + 1, jn - 1
534    
535  SUBROUTINE xadv(imr, jnp, j1, j2, p, ua, js, jn, iml, adx, iad)                  DO i = 1, imr
536    REAL p(imr, jnp), adx(imr, jnp), qtmp(-imr:imr+imr), ua(imr, jnp)                     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              ! ****6***0*********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    jmr = jnp - 1            DO i = 1, imjm
556    DO j = j1, j2               wk1(i, 1, 1) = q(i, 1, k, ic) + 0.5*wk1(i, 1, 1)
557      IF (j>js .AND. j<jn) GO TO 1309               wk1(i, 1, 2) = q(i, 1, k, ic) + 0.5*wk1(i, 1, 2)
558              END DO
559    
560      DO i = 1, imr            IF (cross) THEN
561        qtmp(i) = p(i, j)               ! Add cross terms in the vertical direction.
562      END DO               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      DO i = -iml, 0            CALL xtp(imr, jnp, iml, j1, j2, jn, js, pu, dq(1,1,k,ic), wk1(1,1,2), &
583        qtmp(i) = p(imr+i, j)                 crx, fx1, xmass, iord)
       qtmp(imr+1-i) = p(1-i, j)  
     END DO  
584    
585      IF (iad==2) THEN            CALL ytp(imr, jnp, j1, j2, acosp, rcap, dq(1,1,k,ic), wk1(1,1,1), cry, &
586        DO i = 1, imr                 dc2, ymass, wk1(1,1,3), wk1(1,1,4), wk1(1,1,5), wk1(1,1,6), jord)
         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)  
       END DO  
     ELSE IF (iad==1) THEN  
       DO i = 1, imr  
         iu = ua(i, j)  
         ru = ua(i, j) - iu  
         iiu = i - iu  
         IF (ua(i,j)>=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))  
         END IF  
       END DO  
     END IF  
587    
588      DO i = 1, imr         END DO
       adx(i, j) = adx(i, j) - p(i, j)  
589      END DO      END DO
 1309 END DO  
590    
591    ! Eulerian upwind      ! ******* Compute vertical mass flux (same unit as PS) ***********
592    
593    DO j = js + 1, jn - 1      ! 1st step: compute total column mass CONVERGENCE.
594    
595      DO i = 1, imr      DO j = 1, jnp
596        qtmp(i) = p(i, j)         DO i = 1, imr
597              cry(i, j) = dpi(i, j, 1)
598           END DO
599      END DO      END DO
600    
601      qtmp(0) = p(imr, j)      DO k = 2, nlay
602      qtmp(imr+1) = p(1, j)         DO j = 1, jnp
603              DO i = 1, imr
604      IF (iad==2) THEN               cry(i, j) = cry(i, j) + dpi(i, j, k)
605        qtmp(-1) = p(imr-1, j)            END DO
606        qtmp(imr+2) = p(2, j)         END DO
       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)  
       END DO  
     ELSE IF (iad==1) THEN  
       ! 1st order  
       DO i = 1, imr  
         ip = i - ua(i, j)  
         adx(i, j) = ua(i, j)*(qtmp(ip)-qtmp(ip+1))  
       END DO  
     END IF  
   END DO  
   
   IF (j1/=2) THEN  
     DO i = 1, imr  
       adx(i, 2) = 0.  
       adx(i, jmr) = 0.  
     END DO  
   END IF  
   ! set cross term due to x-adv at the poles to zero.  
   DO i = 1, imr  
     adx(i, 1) = 0.  
     adx(i, jnp) = 0.  
   END DO  
   RETURN  
 END SUBROUTINE xadv  
   
 SUBROUTINE lmtppm(dc, a6, ar, al, p, im, lmt)  
   
   ! 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  
   
   PARAMETER (r12=1./12.)  
   DIMENSION a6(im), ar(im), al(im), p(im), dc(im)  
   
   IF (lmt==0) THEN  
     ! Full constraint  
     DO i = 1, im  
       IF (dc(i)==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<-da2) THEN  
           a6(i) = 3.*(al(i)-p(i))  
           ar(i) = al(i) - a6(i)  
         ELSE IF (a6da>da2) THEN  
           a6(i) = 3.*(ar(i)-p(i))  
           al(i) = ar(i) - a6(i)  
         END IF  
       END IF  
     END DO  
   ELSE IF (lmt==1) THEN  
     ! Semi-monotonic constraint  
     DO i = 1, im  
       IF (abs(ar(i)-al(i))>=-a6(i)) GO TO 150  
       IF (p(i)<ar(i) .AND. p(i)<al(i)) THEN  
         ar(i) = p(i)  
         al(i) = p(i)  
         a6(i) = 0.  
       ELSE IF (ar(i)>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)  
       END IF  
 150 END DO  
   ELSE IF (lmt==2) THEN  
     DO i = 1, im  
       IF (abs(ar(i)-al(i))>=-a6(i)) GO TO 250  
       fmin = p(i) + 0.25*(ar(i)-al(i))**2/a6(i) + a6(i)*r12  
       IF (fmin>=0.) GO TO 250  
       IF (p(i)<ar(i) .AND. p(i)<al(i)) THEN  
         ar(i) = p(i)  
         al(i) = p(i)  
         a6(i) = 0.  
       ELSE IF (ar(i)>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)  
       END IF  
 250 END DO  
   END IF  
   RETURN  
 END SUBROUTINE lmtppm  
   
 SUBROUTINE a2c(u, v, imr, jmr, j1, j2, crx, cry, dtdx5, dtdy5)  
   DIMENSION u(imr, *), v(imr, *), crx(imr, *), cry(imr, *), dtdx5(*)  
   
   DO j = j1, j2  
     DO i = 2, imr  
       crx(i, j) = dtdx5(j)*(u(i,j)+u(i-1,j))  
607      END DO      END DO
   END DO  
608    
609    DO j = j1, j2      DO j = 1, jnp
610      crx(1, j) = dtdx5(j)*(u(1,j)+u(imr,j))         DO i = 1, imr
   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)  
   DIMENSION cosp(*), cose(*)  
   
   jmr = jnp - 1  
   DO j = 2, jnp  
     ph5 = -0.5*pi + (float(j-1)-0.5)*dp  
     cose(j) = cos(ph5)  
   END DO  
   
   jeq = (jnp+1)/2  
   IF (jmr==2*(jmr/2)) THEN  
     DO j = jnp, jeq + 1, -1  
       cose(j) = cose(jnp+2-j)  
     END DO  
   ELSE  
     ! cell edge at equator.  
     cose(jeq+1) = 1.  
     DO j = jnp, jeq + 2, -1  
       cose(j) = cose(jnp+2-j)  
     END DO  
   END IF  
611    
612    DO j = 2, jmr            ! 2nd step: compute PS2 (PS at n+1) using the hydrostatic assumption.
613      cosp(j) = 0.5*(cose(j)+cose(j+1))            ! Changes (increases) to surface pressure = total column mass
614    END DO            ! convergence
   cosp(1) = 0.  
   cosp(jnp) = 0.  
   RETURN  
 END SUBROUTINE cosa  
   
 SUBROUTINE cosc(cosp, cose, jnp, pi, dp)  
   DIMENSION cosp(*), cose(*)  
   
   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)  
   
   PARAMETER (tiny=1.E-60)  
   DIMENSION q(imr, jnp, nlay), qtmp(imr, jnp), cosp(*), acosp(*)  
   LOGICAL cross  
   
   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==0) GO TO 50  
   CALL filew(q(1,1,l), qtmp, imr, jnp, j1, j2, ipx, tiny)  
   IF (ipx==0) GO TO 50  
   
   IF (cross) THEN  
     CALL filcr(q(1,1,l), imr, jnp, j1, j2, cosp, acosp, icr, tiny)  
   END IF  
   IF (icr==0) GO TO 50  
   
   ! Vertical filling...  
   DO i = 1, len  
     IF (q(i,j1,1)<0.) THEN  
       ip = ip + 1  
       q(i, j1, 2) = q(i, j1, 2) + q(i, j1, 1)  
       q(i, j1, 1) = 0.  
     END IF  
   END DO  
615    
616  50 CONTINUE            ps2(i, j) = ps1(i, j) + cry(i, j)
   DO l = 2, nlaym1  
     icr = 1  
   
     CALL filns(q(1,1,l), imr, jnp, j1, j2, cosp, acosp, ipy, tiny)  
     IF (ipy==0) GO TO 225  
     CALL filew(q(1,1,l), qtmp, imr, jnp, j1, j2, ipx, tiny)  
     IF (ipx==0) GO TO 225  
     IF (cross) THEN  
       CALL filcr(q(1,1,l), imr, jnp, j1, j2, cosp, acosp, icr, tiny)  
     END IF  
     IF (icr==0) GO TO 225  
617    
618      DO i = 1, len            ! 3rd step: compute vertical mass flux from mass conservation
619        IF (q(i,j1,l)<0.) THEN            ! principle.
620    
621          ip = ip + 1            w(i, j, 1) = dpi(i, j, 1) - dbk(1)*cry(i, j)
622          ! From above            w(i, j, nlay) = 0.
623          qup = q(i, j1, l-1)         END DO
         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.  
       END IF  
624      END DO      END DO
 225 END DO  
   
   ! BOTTOM LAYER  
   sum = 0.  
   l = nlay  
   
   CALL filns(q(1,1,l), imr, jnp, j1, j2, cosp, acosp, ipy, tiny)  
   IF (ipy==0) GO TO 911  
   CALL filew(q(1,1,l), qtmp, imr, jnp, j1, j2, ipx, tiny)  
   IF (ipx==0) GO TO 911  
   
   CALL filcr(q(1,1,l), imr, jnp, j1, j2, cosp, acosp, icr, tiny)  
   IF (icr==0) GO TO 911  
   
   DO i = 1, len  
     IF (q(i,j1,l)<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.  
     END IF  
   END DO  
   
 911 CONTINUE  
625    
626    IF (ip>imr) THEN      DO k = 2, nlay - 1
627      WRITE (6, *) 'IC=', ic, ' STEP=', nstep, ' Vertical filling pts=', ip         DO j = 1, jnp
628    END IF            DO i = 1, imr
629                 w(i, j, k) = w(i, j, k-1) + dpi(i, j, k) - dbk(k)*cry(i, j)
630    IF (sum>1.E-25) THEN            END DO
631      WRITE (6, *) ic, nstep, ' Mass source from the ground=', sum         END DO
   END IF  
   RETURN  
 END SUBROUTINE qckxyz  
   
 SUBROUTINE filcr(q, imr, jnp, j1, j2, cosp, acosp, icr, tiny)  
   DIMENSION q(imr, *), cosp(*), acosp(*)  
   
   icr = 0  
   DO j = j1 + 1, j2 - 1  
     DO i = 1, imr - 1  
       IF (q(i,j)<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  
       END IF  
     END DO  
     IF (icr==0 .AND. q(imr,j)>=0.) GO TO 65  
     DO i = 2, imr  
       IF (q(i,j)<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  
       END IF  
632      END DO      END DO
     ! *****************************************  
     ! i=1  
     i = 1  
     IF (q(i,j)<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  
     END IF  
     ! *****************************************  
     ! i=IMR  
     i = imr  
     IF (q(i,j)<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  
     END IF  
     ! *****************************************  
 65 END DO  
   
   DO i = 1, imr  
     IF (q(i,j1)<0. .OR. q(i,j2)<0.) THEN  
       icr = 1  
       GO TO 80  
     END IF  
   END DO  
   
 80 CONTINUE  
633    
634    IF (q(1,1)<0. .OR. q(1,jnp)<0.) THEN      DO k = 1, nlay
635      icr = 1         DO j = 1, jnp
636    END IF            DO i = 1, imr
637                 delp2(i, j, k) = dap(k) + dbk(k)*ps2(i, j)
638    RETURN            END DO
639  END SUBROUTINE filcr         END DO
   
 SUBROUTINE filns(q, imr, jnp, j1, j2, cosp, acosp, ipy, tiny)  
   DIMENSION q(imr, *), cosp(*), acosp(*)  
   ! logical first  
   ! data first /.true./  
   ! save cap1  
   
   ! if(first) then  
   dp = 4.*atan(1.)/float(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)<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  
       END IF  
640      END DO      END DO
   END DO  
641    
642    DO i = 1, imr      krd = max(3, kord)
643      IF (q(i,j1)<0.) THEN      DO ic = 1, nc
       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  
     END IF  
   END DO  
644    
645    j = j2         ! ****6***0*********0*********0*********0*********0*********0**********72
   DO i = 1, imr  
     IF (q(i,j)<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  
     END IF  
   END DO  
646    
647    ! Check Poles.         CALL fzppm(imr, jnp, nlay, j1, dq(1,1,1,ic), w, q(1,1,1,ic), wk1, dpi, &
648    IF (q(1,1)<0.) THEN              dc2, crx, cry, pu, pv, xmass, ymass, delp1, krd)
     dq = q(1, 1)*cap1/float(imr)*acosp(j1)  
     DO i = 1, imr  
       q(i, 1) = 0.  
       q(i, j1) = q(i, j1) + dq  
       IF (q(i,j1)<0.) ipy = 1  
     END DO  
   END IF  
649    
   IF (q(1,jnp)<0.) THEN  
     dq = q(1, jnp)*cap1/float(imr)*acosp(j2)  
     DO i = 1, imr  
       q(i, jnp) = 0.  
       q(i, j2) = q(i, j2) + dq  
       IF (q(i,j2)<0.) ipy = 1  
     END DO  
   END IF  
650    
651    RETURN         IF (fill) CALL qckxyz(dq(1,1,1,ic), dc2, imr, jnp, nlay, j1, j2, cosp, &
652  END SUBROUTINE filns              acosp, .FALSE., ic, nstep)
653    
654  SUBROUTINE filew(q, qtmp, imr, jnp, j1, j2, ipx, tiny)         ! Recover tracer mixing ratio from "density" using predicted
655    DIMENSION q(imr, *), qtmp(jnp, imr)         ! "air density" (pressure thickness) at time-level n+1
656    
657    ipx = 0         DO k = 1, nlay
658    ! Copy & swap direction for vectorization.            DO j = 1, jnp
659    DO i = 1, imr               DO i = 1, imr
660      DO j = j1, j2                  q(i, j, k, ic) = dq(i, j, k, ic)/delp2(i, j, k)
661        qtmp(j, i) = q(i, j)               END DO
662      END DO            END DO
663    END DO         END DO
664    
665    DO i = 2, imr - 1         IF (j1/=2) THEN
666      DO j = j1, j2            DO k = 1, nlay
667        IF (qtmp(j,i)<0.) THEN               DO i = 1, imr
668          ipx = 1                  ! j=1 c'est le p\^ole Sud, j=JNP c'est le p\^ole Nord
669          ! west                  q(i, 2, k, ic) = q(i, 1, k, ic)
670          d0 = max(0., qtmp(j,i-1))                  q(i, jmr, k, ic) = q(i, jmp, k, ic)
671          d1 = min(-qtmp(j,i), d0)               END DO
672          qtmp(j, i-1) = qtmp(j, i-1) - d1            END DO
673          qtmp(j, i) = qtmp(j, i) + d1         END IF
         ! 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  
       END IF  
674      END DO      END DO
   END DO  
   
   i = 1  
   DO j = j1, j2  
     IF (qtmp(j,i)<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  
     END IF  
   END DO  
   i = imr  
   DO j = j1, j2  
     IF (qtmp(j,i)<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  
675    
676        qtmp(j, i) = qtmp(j, i) + d2 + tiny      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      END IF
   END DO  
684    
685    IF (ipx/=0) THEN    END SUBROUTINE ppm3d
     DO j = j1, j2  
       DO i = 1, imr  
         q(i, j) = qtmp(j, i)  
       END DO  
     END DO  
   ELSE  
686    
687      ! Poles.  end module ppm3d_m
     IF (q(1,1)<0 .OR. q(1,jnp)<0.) ipx = 1  
   END IF  
   RETURN  
 END SUBROUTINE filew  
Legend:
Removed from v.81  
changed lines
  Added in v.169
  ViewVC Help
Powered by ViewVC 1.1.21