/[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.f revision 76 by guez, Fri Nov 15 18:45:49 2013 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  ! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/ppm3d.F,v 1.1.1.1 2004/05/19 12:53:07 lmdzadmin Exp $  
3  !    implicit none
4    
5  cFrom lin@explorer.gsfc.nasa.gov Wed Apr 15 17:44:44 1998  contains
6  cDate: Wed, 15 Apr 1998 11:37:03 -0400  
7  cFrom: lin@explorer.gsfc.nasa.gov    SUBROUTINE ppm3d(igd, q, ps1, ps2, u, v, w, ndt, iord, jord, kord, nc, imr, &
8  cTo: Frederic.Hourdin@lmd.jussieu.fr         jnp, j1, nlay, ap, bp, pt, ae, fill, umax)
9  cSubject: 3D transport module of the GSFC CTM and GEOS GCM  
10        ! INPUT:
11        ! =============
12  cThis code is sent to you by S-J Lin, DAO, NASA-GSFC  
13        ! Q(IMR,JNP,NLAY,NC): mixing ratios at current time (t)
14  cNote: this version is intended for machines like CRAY      ! NC: total number of constituents
15  C-90. No multitasking directives implemented.      ! IMR: first dimension (E-W); number of Grid intervals in E-W is IMR
16        ! JNP: 2nd dimension (N-S); number of Grid intervals in N-S is JNP-1
17              ! NLAY: 3rd dimension (number of layers); vertical index increases from 1
18  C ********************************************************************      ! at
19  C      ! the model top to NLAY near the surface (see fig. below).
20  C TransPort Core for Goddard Chemistry Transport Model (G-CTM), Goddard      ! It is assumed that 6 <= NLAY <= JNP (for dynamic memory allocation)
21  C Earth Observing System General Circulation Model (GEOS-GCM), and Data  
22  C Assimilation System (GEOS-DAS).      ! PS1(IMR,JNP): surface pressure at current time (t)
23  C      ! PS2(IMR,JNP): surface pressure at mid-time-level (t+NDT/2)
24  C ********************************************************************      ! PS2 is replaced by the predicted PS (at t+NDT) on output.
25  C      ! Note: surface pressure can have any unit or can be multiplied by any
26  C Purpose: given horizontal winds on  a hybrid sigma-p surfaces,      ! const.
27  C          one call to tpcore updates the 3-D mixing ratio  
28  C          fields one time step (NDT). [vertical mass flux is computed      ! The pressure at layer edges are defined as follows:
29  C          internally consistent with the discretized hydrostatic mass  
30  C          continuity equation of the C-Grid GEOS-GCM (for IGD=1)].      ! p(i,j,k) = AP(k)*PT  +  BP(k)*PS(i,j)          (1)
31  C  
32  C Schemes: Multi-dimensional Flux Form Semi-Lagrangian (FFSL) scheme based      ! Where PT is a constant having the same unit as PS.
33  C          on the van Leer or PPM.      ! AP and BP are unitless constants given at layer edges
34  C          (see Lin and Rood 1996).      ! defining the vertical coordinate.
35  C Version 4.5      ! BP(1) = 0., BP(NLAY+1) = 1.
36  C Last modified: Dec. 5, 1996      ! The pressure at the model top is PTOP = AP(1)*PT
37  C Major changes from version 4.0: a more general vertical hybrid sigma-  
38  C pressure coordinate.      ! For pure sigma system set AP(k) = 1 for all k, PT = PTOP,
39  C Subroutines modified: xtp, ytp, fzppm, qckxyz      ! BP(k) = sige(k) (sigma at edges), PS = Psfc - PTOP.
40  C Subroutines deleted: vanz  
41  C      ! Note: the sigma-P coordinate is a subset of Eq. 1, which in turn
42  C Author: Shian-Jiann Lin      ! is a subset of the following even more general sigma-P-thelta coord.
43  C mail address:      ! currently under development.
44  C                 Shian-Jiann Lin*      ! p(i,j,k) = (AP(k)*PT + BP(k)*PS(i,j))/(D(k)-C(k)*TE**(-1/kapa))
45  C                 Code 910.3, NASA/GSFC, Greenbelt, MD 20771  
46  C                 Phone: 301-286-9540      ! Cf. ppm3d.txt.
47  C                 E-mail: lin@dao.gsfc.nasa.gov  
48  C      ! U(IMR,JNP,NLAY) & V(IMR,JNP,NLAY):winds (m/s) at mid-time-level (t+NDT/2)
49  C *affiliation:      ! U and V may need to be polar filtered in advance in some cases.
50  C                 Joint Center for Earth Systems Technology  
51  C                 The University of Maryland Baltimore County      ! IGD:      grid type on which winds are defined.
52  C                 NASA - Goddard Space Flight Center      ! IGD = 0:  A-Grid  [all variables defined at the same point from south
53  C References:      ! pole (j=1) to north pole (j=JNP) ]
54  C  
55  C 1. Lin, S.-J., and R. B. Rood, 1996: Multidimensional flux form semi-      ! IGD = 1  GEOS-GCM C-Grid
56  C    Lagrangian transport schemes. Mon. Wea. Rev., 124, 2046-2070.      ! Cf. ppm3d.txt.
57  C  
58  C 2. Lin, S.-J., W. C. Chao, Y. C. Sud, and G. K. Walker, 1994: A class of      ! U(i,  1) is defined at South Pole.
59  C    the van Leer-type transport schemes and its applications to the moist-      ! V(i,  1) is half grid north of the South Pole.
60  C    ure transport in a General Circulation Model. Mon. Wea. Rev., 122,      ! V(i,JMR) is half grid south of the North Pole.
61  C    1575-1593.  
62  C      ! V must be defined at j=1 and j=JMR if IGD=1
63  C ****6***0*********0*********0*********0*********0*********0**********72      ! V at JNP need not be given.
64  C  
65        subroutine ppm3d(IGD,Q,PS1,PS2,U,V,W,NDT,IORD,JORD,KORD,NC,IMR,      ! NDT: time step in seconds (need not be constant during the course of
66       &                  JNP,j1,NLAY,AP,BP,PT,AE,fill,dum,Umax)      ! 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  c      implicit none      ! has a well-resolved stratosphere.
69    
70  c     rajout de déclarations      ! J1 defines the size of the polar cap:
71  c      integer Jmax,kmax,ndt0,nstep,k,j,i,ic,l,js,jn,imh,iad,jad,krd      ! South polar cap edge is located at -90 + (j1-1.5)*180/(JNP-1) deg.
72  c      integer iu,iiu,j2,jmr,js0,jt      ! North polar cap edge is located at  90 - (j1-1.5)*180/(JNP-1) deg.
73  c      real dtdy,dtdy5,rcap,iml,jn0,imjm,pi,dl,dp      ! There are currently only two choices (j1=2 or 3).
74  c      real dt,cr1,maxdt,ztc,d5,sum1,sum2,ru      ! IMR must be an even integer if j1 = 2. Recommended value: J1=3.
75  C  
76  C ********************************************************************      ! IORD, JORD, and KORD are integers controlling various options in E-W,
77  C      ! N-S,
78  C =============      ! and vertical transport, respectively. Recommended values for positive
79  C INPUT:      ! definite scalars: IORD=JORD=3, KORD=5. Use KORD=3 for non-
80  C =============      ! positive definite scalars or when linear correlation between constituents
81  C      ! is to be maintained.
82  C Q(IMR,JNP,NLAY,NC): mixing ratios at current time (t)  
83  C NC: total number of constituents      ! _ORD=
84  C IMR: first dimension (E-W); number of Grid intervals in E-W is IMR      ! 1: 1st order upstream scheme (too diffusive, not a useful option; it
85  C JNP: 2nd dimension (N-S); number of Grid intervals in N-S is JNP-1      ! can be used for debugging purposes; this is THE only known "linear"
86  C NLAY: 3rd dimension (number of layers); vertical index increases from 1 at      ! monotonic advection scheme.).
87  C       the model top to NLAY near the surface (see fig. below).      ! 2: 2nd order van Leer (full monotonicity constraint;
88  C       It is assumed that 6 <= NLAY <= JNP (for dynamic memory allocation)      ! see Lin et al 1994, MWR)
89  C      ! 3: monotonic PPM* (slightly improved PPM of Collela & Woodward 1984)
90  C PS1(IMR,JNP): surface pressure at current time (t)      ! 4: semi-monotonic PPM (same as 3, but overshoots are allowed)
91  C PS2(IMR,JNP): surface pressure at mid-time-level (t+NDT/2)      ! 5: positive-definite PPM (constraint on the subgrid distribution is
92  C PS2 is replaced by the predicted PS (at t+NDT) on output.      ! only strong enough to prevent generation of negative values;
93  C Note: surface pressure can have any unit or can be multiplied by any      ! both overshoots & undershoots are possible).
94  C       const.      ! 6: un-constrained PPM (nearly diffusion free; slightly faster but
95  C      ! positivity not quaranteed. Use this option only when the fields
96  C The pressure at layer edges are defined as follows:      ! and winds are very smooth).
97  C  
98  C        p(i,j,k) = AP(k)*PT  +  BP(k)*PS(i,j)          (1)      ! *PPM: Piece-wise Parabolic Method
99  C  
100  C Where PT is a constant having the same unit as PS.      ! Note that KORD <=2 options are no longer supported. DO not use option 4
101  C AP and BP are unitless constants given at layer edges      ! or 5.
102  C defining the vertical coordinate.      ! for non-positive definite scalars (such as Ertel Potential Vorticity).
103  C BP(1) = 0., BP(NLAY+1) = 1.  
104  C The pressure at the model top is PTOP = AP(1)*PT      ! The implicit numerical diffusion decreases as _ORD increases.
105  C      ! The last two options (ORDER=5, 6) should only be used when there is
106  C For pure sigma system set AP(k) = 1 for all k, PT = PTOP,      ! significant explicit diffusion (such as a turbulence parameterization).
107  C BP(k) = sige(k) (sigma at edges), PS = Psfc - PTOP.      ! You
108  C      ! might get dispersive results otherwise.
109  C Note: the sigma-P coordinate is a subset of Eq. 1, which in turn      ! No filter of any kind is applied to the constituent fields here.
110  C is a subset of the following even more general sigma-P-thelta coord.  
111  C currently under development.      ! AE: Radius of the sphere (meters).
112  C  p(i,j,k) = (AP(k)*PT + BP(k)*PS(i,j))/(D(k)-C(k)*TE**(-1/kapa))      ! Recommended value for the planet earth: 6.371E6
113  C  
114  C                  /////////////////////////////////      ! fill(logical):   flag to do filling for negatives (see note below).
115  C              / \ ------------- PTOP --------------  AP(1), BP(1)  
116  C               |      ! Umax: Estimate (upper limit) of the maximum U-wind speed (m/s).
117  C    delp(1)    |  ........... Q(i,j,1) ............        ! (220 m/s is a good value for troposphere model; 280 m/s otherwise)
118  C               |  
119  C      W(1)    \ / ---------------------------------  AP(2), BP(2)      ! =============
120  C      ! Output
121  C      ! =============
122  C  
123  C     W(k-1)   / \ ---------------------------------  AP(k), BP(k)      ! Q: mixing ratios at future time (t+NDT) (original values are
124  C               |      ! over-written)
125  C    delp(K)    |  ........... Q(i,j,k) ............      ! W(NLAY): large-scale vertical mass flux as diagnosed from the hydrostatic
126  C               |      ! relationship. W will have the same unit as PS1 and PS2 (eg, mb).
127  C      W(k)    \ / ---------------------------------  AP(k+1), BP(k+1)      ! W must be divided by NDT to get the correct mass-flux unit.
128  C      ! The vertical Courant number C = W/delp_UPWIND, where delp_UPWIND
129  C      ! is the pressure thickness in the "upwind" direction. For example,
130  C      ! C(k) = W(k)/delp(k)   if W(k) > 0;
131  C              / \ ---------------------------------  AP(NLAY), BP(NLAY)      ! C(k) = W(k)/delp(k+1) if W(k) < 0.
132  C               |      ! ( W > 0 is downward, ie, toward surface)
133  C  delp(NLAY)   |  ........... Q(i,j,NLAY) .........        ! PS2: predicted PS at t+NDT (original values are over-written)
134  C               |  
135  C   W(NLAY)=0  \ / ------------- surface ----------- AP(NLAY+1), BP(NLAY+1)      ! ********************************************************************
136  C                 //////////////////////////////////      ! NOTES:
137  C      ! This forward-in-time upstream-biased transport scheme reduces to
138  C U(IMR,JNP,NLAY) & V(IMR,JNP,NLAY):winds (m/s) at mid-time-level (t+NDT/2)      ! the 2nd order center-in-time center-in-space mass continuity eqn.
139  C U and V may need to be polar filtered in advance in some cases.      ! if Q = 1 (constant fields will remain constant). This also ensures
140  C      ! that the computed vertical velocity to be identical to GEOS-1 GCM
141  C IGD:      grid type on which winds are defined.      ! for on-line transport.
142  C IGD = 0:  A-Grid  [all variables defined at the same point from south  
143  C                   pole (j=1) to north pole (j=JNP) ]      ! A larger polar cap is used if j1=3 (recommended for C-Grid winds or when
144  C      ! winds are noisy near poles).
145  C IGD = 1  GEOS-GCM C-Grid  
146  C                                      [North]      ! Flux-Form Semi-Lagrangian transport in the East-West direction is used
147  C      ! when and where Courant number is greater than one.
148  C                                       V(i,j)  
149  C                                          |      ! The user needs to change the parameter Jmax or Kmax if the resolution
150  C                                          |      ! is greater than 0.5 deg in N-S or 150 layers in the vertical direction.
151  C                                          |      ! (this TransPort Core is otherwise resolution independent and can be used
152  C                             U(i-1,j)---Q(i,j)---U(i,j) [EAST]      ! as a library routine).
153  C                                          |  
154  C                                          |      ! PPM is 4th order accurate when grid spacing is uniform (x & y); 3rd
155  C                                          |      ! order accurate for non-uniform grid (vertical sigma coord.).
156  C                                       V(i,j-1)  
157  C      ! Time step is limitted only by transport in the meridional direction.
158  C         U(i,  1) is defined at South Pole.      ! (the FFSL scheme is not implemented in the meridional direction).
159  C         V(i,  1) is half grid north of the South Pole.  
160  C         V(i,JMR) is half grid south of the North Pole.      ! Since only 1-D limiters are applied, negative values could
161  C      ! potentially be generated when large time step is used and when the
162  C         V must be defined at j=1 and j=JMR if IGD=1      ! initial fields contain discontinuities.
163  C         V at JNP need not be given.      ! This does not necessarily imply the integration is unstable.
164  C      ! These negatives are typically very small. A filling algorithm is
165  C NDT: time step in seconds (need not be constant during the course of      ! activated if the user set "fill" to be true.
166  C      the integration). Suggested value: 30 min. for 4x5, 15 min. for 2x2.5  
167  C      (Lat-Lon) resolution. Smaller values are recommanded if the model      ! The van Leer scheme used here is nearly as accurate as the original PPM
168  C      has a well-resolved stratosphere.      ! due to the use of a 4th order accurate reference slope. The PPM imple-
169  C      ! mented here is an improvement over the original and is also based on
170  C J1 defines the size of the polar cap:      ! the 4th order reference slope.
171  C South polar cap edge is located at -90 + (j1-1.5)*180/(JNP-1) deg.  
172  C North polar cap edge is located at  90 - (j1-1.5)*180/(JNP-1) deg.      ! ****6***0*********0*********0*********0*********0*********0**********72
173  C There are currently only two choices (j1=2 or 3).  
174  C IMR must be an even integer if j1 = 2. Recommended value: J1=3.      ! User modifiable parameters
175  C  
176  C IORD, JORD, and KORD are integers controlling various options in E-W, N-S,      integer, PARAMETER:: jmax=361, kmax=150
177  C and vertical transport, respectively. Recommended values for positive  
178  C definite scalars: IORD=JORD=3, KORD=5. Use KORD=3 for non-      ! ****6***0*********0*********0*********0*********0*********0**********72
179  C positive definite scalars or when linear correlation between constituents  
180  C is to be maintained.      ! Input-Output arrays
181  C  
182  C  _ORD=      integer imr
183  C        1: 1st order upstream scheme (too diffusive, not a useful option; it      INTEGER igd, iord, jord, kord, nc, jnp, j1, nlay, ae
184  C           can be used for debugging purposes; this is THE only known "linear"      REAL q(imr, jnp, nlay, nc), ps1(imr, jnp), ps2(imr, jnp), &
185  C           monotonic advection scheme.).           u(imr, jnp, nlay), v(imr, jnp, nlay), ap(nlay+1), bp(nlay+1), &
186  C        2: 2nd order van Leer (full monotonicity constraint;           w(imr, jnp, nlay), ndt, umax
187  C           see Lin et al 1994, MWR)      INTEGER imrd2
188  C        3: monotonic PPM* (slightly improved PPM of Collela & Woodward 1984)      REAL pt
189  C        4: semi-monotonic PPM (same as 3, but overshoots are allowed)      LOGICAL cross, fill
190  C        5: positive-definite PPM (constraint on the subgrid distribution is  
191  C           only strong enough to prevent generation of negative values;      ! Local dynamic arrays
192  C           both overshoots & undershoots are possible).  
193  C        6: un-constrained PPM (nearly diffusion free; slightly faster but      REAL crx(imr, jnp), cry(imr, jnp), xmass(imr, jnp), ymass(imr, jnp), &
194  C           positivity not quaranteed. Use this option only when the fields           fx1(imr+1), dpi(imr, jnp, nlay), delp1(imr, jnp, nlay), &
195  C           and winds are very smooth).           wk1(imr, jnp, nlay), pu(imr, jnp), pv(imr, jnp), dc2(imr, jnp), &
196  C           delp2(imr, jnp, nlay), dq(imr, jnp, nlay, nc), va(imr, jnp), &
197  C *PPM: Piece-wise Parabolic Method           ua(imr, jnp), qtmp(-imr:2*imr)
198  C  
199  C Note that KORD <=2 options are no longer supported. DO not use option 4 or 5.      ! Local static  arrays
200  C for non-positive definite scalars (such as Ertel Potential Vorticity).  
201  C      REAL dtdx(jmax), dtdx5(jmax), acosp(jmax), cosp(jmax), cose(jmax), &
202  C The implicit numerical diffusion decreases as _ORD increases.           dap(kmax), dbk(kmax)
203  C The last two options (ORDER=5, 6) should only be used when there is      integer ndt0, nstep
204  C significant explicit diffusion (such as a turbulence parameterization). You      DATA ndt0, nstep/0, 0/
205  C might get dispersive results otherwise.      DATA cross/.TRUE./
206  C No filter of any kind is applied to the constituent fields here.      SAVE dtdy, dtdy5, rcap, js0, jn0, iml, dtdx, dtdx5, acosp, cosp, cose, dap, dbk
207  C      real cr1, d5, dl, dp, dt, dtdy, dtdy5, pi, rcap, ru, sum1, sum2, ztc
208  C AE: Radius of the sphere (meters).      integer i, iad, ic, iiu, imh, imjm, iml, iu, j, j2, jad, jmp, jmr, jn, jn0
209  C     Recommended value for the planet earth: 6.371E6      integer js, js0, jt, k, krd, l, maxdt
210  C  
211  C fill(logical):   flag to do filling for negatives (see note below).      jmr = jnp - 1
212  C      imjm = imr*jnp
213  C Umax: Estimate (upper limit) of the maximum U-wind speed (m/s).      j2 = jnp - j1 + 1
214  C (220 m/s is a good value for troposphere model; 280 m/s otherwise)      nstep = nstep + 1
215  C  
216  C =============      ! *********** Initialization **********************
217  C Output      IF (nstep==1) THEN
218  C =============  
219  C         WRITE (6, *) '------------------------------------ '
220  C Q: mixing ratios at future time (t+NDT) (original values are over-written)         WRITE (6, *) 'NASA/GSFC Transport Core Version 4.5'
221  C W(NLAY): large-scale vertical mass flux as diagnosed from the hydrostatic         WRITE (6, *) '------------------------------------ '
222  C          relationship. W will have the same unit as PS1 and PS2 (eg, mb).  
223  C          W must be divided by NDT to get the correct mass-flux unit.         WRITE (6, *) 'IMR=', imr, ' JNP=', jnp, ' NLAY=', nlay, ' j1=', j1
224  C          The vertical Courant number C = W/delp_UPWIND, where delp_UPWIND         WRITE (6, *) 'NC=', nc, iord, jord, kord, ndt
225  C          is the pressure thickness in the "upwind" direction. For example,  
226  C          C(k) = W(k)/delp(k)   if W(k) > 0;         ! controles sur les parametres
227  C          C(k) = W(k)/delp(k+1) if W(k) < 0.         IF (nlay<6) THEN
228  C              ( W > 0 is downward, ie, toward surface)            WRITE (6, *) 'NLAY must be >= 6'
229  C PS2: predicted PS at t+NDT (original values are over-written)            STOP
230  C         END IF
231  C ********************************************************************         IF (jnp<nlay) THEN
232  C NOTES:            WRITE (6, *) 'JNP must be >= NLAY'
233  C This forward-in-time upstream-biased transport scheme reduces to            STOP
234  C the 2nd order center-in-time center-in-space mass continuity eqn.         END IF
235  C if Q = 1 (constant fields will remain constant). This also ensures         imrd2 = mod(imr, 2)
236  C that the computed vertical velocity to be identical to GEOS-1 GCM         IF (j1==2 .AND. imrd2/=0) THEN
237  C for on-line transport.            WRITE (6, *) 'if j1=2 IMR must be an even integer'
238  C            STOP
239  C A larger polar cap is used if j1=3 (recommended for C-Grid winds or when         END IF
240  C winds are noisy near poles).  
241  C  
242  C Flux-Form Semi-Lagrangian transport in the East-West direction is used         IF (jmax<jnp .OR. kmax<nlay) THEN
243  C when and where Courant number is greater than one.            WRITE (6, *) 'Jmax or Kmax is too small'
244  C            STOP
245  C The user needs to change the parameter Jmax or Kmax if the resolution         END IF
246  C is greater than 0.5 deg in N-S or 150 layers in the vertical direction.  
247  C (this TransPort Core is otherwise resolution independent and can be used         DO k = 1, nlay
248  C as a library routine).            dap(k) = (ap(k+1)-ap(k))*pt
249  C            dbk(k) = bp(k+1) - bp(k)
250  C PPM is 4th order accurate when grid spacing is uniform (x & y); 3rd         END DO
251  C order accurate for non-uniform grid (vertical sigma coord.).  
252  C         pi = 4.*atan(1.)
253  C Time step is limitted only by transport in the meridional direction.         dl = 2.*pi/float(imr)
254  C (the FFSL scheme is not implemented in the meridional direction).         dp = pi/float(jmr)
255  C  
256  C Since only 1-D limiters are applied, negative values could         IF (igd==0) THEN
257  C potentially be generated when large time step is used and when the            ! Compute analytic cosine at cell edges
258  C initial fields contain discontinuities.            CALL cosa(cosp, cose, jnp, pi, dp)
259  C This does not necessarily imply the integration is unstable.         ELSE
260  C These negatives are typically very small. A filling algorithm is            ! Define cosine consistent with GEOS-GCM (using dycore2.0 or later)
261  C activated if the user set "fill" to be true.            CALL cosc(cosp, cose, jnp, pi, dp)
262  C         END IF
263  C The van Leer scheme used here is nearly as accurate as the original PPM  
264  C due to the use of a 4th order accurate reference slope. The PPM imple-         DO j = 2, jmr
265  C mented here is an improvement over the original and is also based on            acosp(j) = 1./cosp(j)
266  C the 4th order reference slope.         END DO
267  C  
268  C ****6***0*********0*********0*********0*********0*********0**********72         ! Inverse of the Scaled polar cap area.
269  C  
270  C     User modifiable parameters         rcap = dp/(imr*(1.-cos((j1-1.5)*dp)))
271  C         acosp(1) = rcap
272        parameter (Jmax = 361, kmax = 150)         acosp(jnp) = rcap
273  C      END IF
274  C ****6***0*********0*********0*********0*********0*********0**********72  
275  C      IF (ndt0/=ndt) THEN
276  C Input-Output arrays         dt = ndt
277  C         ndt0 = ndt
278          
279        real Q(IMR,JNP,NLAY,NC),PS1(IMR,JNP),PS2(IMR,JNP),         IF (umax<180.) THEN
280       &     U(IMR,JNP,NLAY),V(IMR,JNP,NLAY),AP(NLAY+1),            WRITE (6, *) 'Umax may be too small!'
281       &     BP(NLAY+1),W(IMR,JNP,NLAY),NDT,val(NLAY),Umax         END IF
282        integer IGD,IORD,JORD,KORD,NC,IMR,JNP,j1,NLAY,AE         cr1 = abs(umax*dt)/(dl*ae)
283        integer IMRD2         maxdt = dp*ae/abs(umax) + 0.5
284        real    PT               WRITE (6, *) 'Largest time step for max(V)=', umax, ' is ', maxdt
285        logical  cross, fill, dum         IF (maxdt<abs(ndt)) THEN
286  C            WRITE (6, *) 'Warning!!! NDT maybe too large!'
287  C Local dynamic arrays         END IF
288  C  
289        real CRX(IMR,JNP),CRY(IMR,JNP),xmass(IMR,JNP),ymass(IMR,JNP),         IF (cr1>=0.95) THEN
290       &     fx1(IMR+1),DPI(IMR,JNP,NLAY),delp1(IMR,JNP,NLAY),            js0 = 0
291       &     WK1(IMR,JNP,NLAY),PU(IMR,JNP),PV(IMR,JNP),DC2(IMR,JNP),            jn0 = 0
292       &     delp2(IMR,JNP,NLAY),DQ(IMR,JNP,NLAY,NC),VA(IMR,JNP),            iml = imr - 2
293       &     UA(IMR,JNP),qtmp(-IMR:2*IMR)            ztc = 0.
294  C         ELSE
295  C Local static  arrays            ztc = acos(cr1)*(180./pi)
296  C  
297        real DTDX(Jmax), DTDX5(Jmax), acosp(Jmax),            js0 = float(jmr)*(90.-ztc)/180. + 2
298       &     cosp(Jmax), cose(Jmax), DAP(kmax),DBK(Kmax)            js0 = max(js0, j1+1)
299        data NDT0, NSTEP /0, 0/            iml = min(6*js0/(j1-1)+2, 4*imr/5)
300        data cross /.true./            jn0 = jnp - js0 + 1
301        SAVE DTDY, DTDY5, RCAP, JS0, JN0, IML,         END IF
302       &     DTDX, DTDX5, ACOSP, COSP, COSE, DAP,DBK  
303  C  
304                       DO j = 2, jmr
305        JMR = JNP -1            dtdx(j) = dt/(dl*ae*cosp(j))
306        IMJM  = IMR*JNP  
307        j2 = JNP - j1 + 1            dtdx5(j) = 0.5*dtdx(j)
308        NSTEP = NSTEP + 1         END DO
309  C  
310  C *********** Initialization **********************  
311        if(NSTEP.eq.1) then         dtdy = dt/(ae*dp)
312  c         dtdy5 = 0.5*dtdy
313        write(6,*) '------------------------------------ '  
314        write(6,*) 'NASA/GSFC Transport Core Version 4.5'      END IF
315        write(6,*) '------------------------------------ '  
316  c      ! *********** End Initialization **********************
317        WRITE(6,*) 'IMR=',IMR,' JNP=',JNP,' NLAY=',NLAY,' j1=',j1  
318        WRITE(6,*) 'NC=',NC,IORD,JORD,KORD,NDT      ! delp = pressure thickness: the psudo-density in a hydrostatic system.
319  C      DO k = 1, nlay
320  C controles sur les parametres         DO j = 1, jnp
321        if(NLAY.LT.6) then            DO i = 1, imr
322          write(6,*) 'NLAY must be >= 6'               delp1(i, j, k) = dap(k) + dbk(k)*ps1(i, j)
323          stop               delp2(i, j, k) = dap(k) + dbk(k)*ps2(i, j)
324        endif            END DO
325        if (JNP.LT.NLAY) then         END DO
326           write(6,*) 'JNP must be >= NLAY'      END DO
327          stop  
328        endif  
329        IMRD2=mod(IMR,2)      IF (j1/=2) THEN
330        if (j1.eq.2.and.IMRD2.NE.0) then         DO ic = 1, nc
331           write(6,*) 'if j1=2 IMR must be an even integer'            DO l = 1, nlay
332          stop               DO i = 1, imr
333        endif                  q(i, 2, l, ic) = q(i, 1, l, ic)
334                    q(i, jmr, l, ic) = q(i, jnp, l, ic)
335  C               END DO
336        if(Jmax.lt.JNP .or. Kmax.lt.NLAY) then            END DO
337          write(6,*) 'Jmax or Kmax is too small'         END DO
338          stop      END IF
339        endif  
340  C      ! Compute "tracer density"
341        DO k=1,NLAY      DO ic = 1, nc
342        DAP(k) = (AP(k+1) - AP(k))*PT         DO k = 1, nlay
343        DBK(k) =  BP(k+1) - BP(k)            DO j = 1, jnp
344        ENDDO                   DO i = 1, imr
345  C                  dq(i, j, k, ic) = q(i, j, k, ic)*delp1(i, j, k)
346        PI = 4. * ATAN(1.)               END DO
347        DL = 2.*PI / float(IMR)            END DO
348        DP =    PI / float(JMR)         END DO
349  C      END DO
350        if(IGD.eq.0) then  
351  C Compute analytic cosine at cell edges      DO k = 1, nlay
352              call cosa(cosp,cose,JNP,PI,DP)  
353        else         IF (igd==0) THEN
354  C Define cosine consistent with GEOS-GCM (using dycore2.0 or later)            ! Convert winds on A-Grid to Courant number on C-Grid.
355              call cosc(cosp,cose,JNP,PI,DP)            CALL a2c(u(1,1,k), v(1,1,k), imr, jmr, j1, j2, crx, cry, dtdx5, dtdy5)
356        endif         ELSE
357  C            ! Convert winds on C-grid to Courant number
358        do 15 J=2,JMR            DO j = j1, j2
359  15    acosp(j) = 1. / cosp(j)               DO i = 2, imr
360  C                  crx(i, j) = dtdx(j)*u(i-1, j, k)
361  C Inverse of the Scaled polar cap area.               END DO
362  C            END DO
363        RCAP  = DP / (IMR*(1.-COS((j1-1.5)*DP)))  
364        acosp(1)   = RCAP  
365        acosp(JNP) = RCAP            DO j = j1, j2
366        endif               crx(1, j) = dtdx(j)*u(imr, j, k)
367  C            END DO
368        if(NDT0 .ne. NDT) then  
369        DT   = NDT            DO i = 1, imr*jmr
370        NDT0 = NDT               cry(i, 2) = dtdy*v(i, 1, k)
371              END DO
372          if(Umax .lt. 180.) then         END IF
373           write(6,*) 'Umax may be too small!'  
374          endif         ! Determine JS and JN
375        CR1  = abs(Umax*DT)/(DL*AE)         js = j1
376        MaxDT = DP*AE / abs(Umax) + 0.5         jn = j2
377        write(6,*)'Largest time step for max(V)=',Umax,' is ',MaxDT  
378        if(MaxDT .lt. abs(NDT)) then         DO j = js0, j1 + 1, -1
379              write(6,*) 'Warning!!! NDT maybe too large!'            DO i = 1, imr
380        endif               IF (abs(crx(i,j))>1.) THEN
381  C                  js = j
382        if(CR1.ge.0.95) then                  GO TO 2222
383        JS0 = 0               END IF
384        JN0 = 0            END DO
385        IML = IMR-2         END DO
386        ZTC = 0.  
387        else  2222   CONTINUE
388        ZTC  = acos(CR1) * (180./PI)         DO j = jn0, j2 - 1
389  C            DO i = 1, imr
390        JS0 = float(JMR)*(90.-ZTC)/180. + 2               IF (abs(crx(i,j))>1.) THEN
391        JS0 = max(JS0, J1+1)                  jn = j
392        IML = min(6*JS0/(J1-1)+2, 4*IMR/5)                  GO TO 2233
393        JN0 = JNP-JS0+1               END IF
394        endif            END DO
395  C             END DO
396  C  2233   CONTINUE
397        do J=2,JMR  
398        DTDX(j)  = DT / ( DL*AE*COSP(J) )         IF (j1/=2) THEN ! Enlarged polar cap.
399              DO i = 1, imr
400        DTDX5(j) = 0.5*DTDX(j)               dpi(i, 2, k) = 0.
401        enddo               dpi(i, jmr, k) = 0.
402  C            END DO
403                 END IF
404        DTDY  = DT /(AE*DP)  
405        DTDY5 = 0.5*DTDY         ! ******* Compute horizontal mass fluxes ************
406  C  
407        endif         ! N-S component
408  C         DO j = j1, j2 + 1
409  C *********** End Initialization **********************            d5 = 0.5*cose(j)
410  C            DO i = 1, imr
411  C delp = pressure thickness: the psudo-density in a hydrostatic system.               ymass(i, j) = cry(i, j)*d5*(delp2(i,j,k)+delp2(i,j-1,k))
412        do  k=1,NLAY            END DO
413           do  j=1,JNP         END DO
414              do  i=1,IMR  
415                 delp1(i,j,k)=DAP(k)+DBK(k)*PS1(i,j)         DO j = j1, j2
416                 delp2(i,j,k)=DAP(k)+DBK(k)*PS2(i,j)                  DO i = 1, imr
417              enddo               dpi(i, j, k) = (ymass(i,j)-ymass(i,j+1))*acosp(j)
418           enddo            END DO
419        enddo         END DO
420              
421  C         ! Poles
422        if(j1.ne.2) then         sum1 = ymass(imr, j1)
423        DO 40 IC=1,NC         sum2 = ymass(imr, j2+1)
424        DO 40 L=1,NLAY         DO i = 1, imr - 1
425        DO 40 I=1,IMR            sum1 = sum1 + ymass(i, j1)
426        Q(I,  2,L,IC) = Q(I,  1,L,IC)            sum2 = sum2 + ymass(i, j2+1)
427  40    Q(I,JMR,L,IC) = Q(I,JNP,L,IC)         END DO
428        endif  
429  C         sum1 = -sum1*rcap
430  C Compute "tracer density"         sum2 = sum2*rcap
431        DO 550 IC=1,NC         DO i = 1, imr
432        DO 44 k=1,NLAY            dpi(i, 1, k) = sum1
433        DO 44 j=1,JNP            dpi(i, jnp, k) = sum2
434        DO 44 i=1,IMR         END DO
435  44    DQ(i,j,k,IC) = Q(i,j,k,IC)*delp1(i,j,k)  
436  550     continue         ! E-W component
437  C  
438        do 1500 k=1,NLAY         DO j = j1, j2
439  C            DO i = 2, imr
440        if(IGD.eq.0) then               pu(i, j) = 0.5*(delp2(i,j,k)+delp2(i-1,j,k))
441  C Convert winds on A-Grid to Courant number on C-Grid.            END DO
442        call A2C(U(1,1,k),V(1,1,k),IMR,JMR,j1,j2,CRX,CRY,dtdx5,DTDY5)         END DO
443        else  
444  C Convert winds on C-grid to Courant number         DO j = j1, j2
445        do 45 j=j1,j2            pu(1, j) = 0.5*(delp2(1,j,k)+delp2(imr,j,k))
446        do 45 i=2,IMR         END DO
447  45    CRX(i,J) = dtdx(j)*U(i-1,j,k)  
448             DO j = j1, j2
449  C            DO i = 1, imr
450        do 50 j=j1,j2               xmass(i, j) = pu(i, j)*crx(i, j)
451  50    CRX(1,J) = dtdx(j)*U(IMR,j,k)            END DO
452  C         END DO
453        do 55 i=1,IMR*JMR  
454  55    CRY(i,2) = DTDY*V(i,1,k)         DO j = j1, j2
455        endif            DO i = 1, imr - 1
456  C                   dpi(i, j, k) = dpi(i, j, k) + xmass(i, j) - xmass(i+1, j)
457  C Determine JS and JN            END DO
458        JS = j1         END DO
459        JN = j2  
460  C         DO j = j1, j2
461        do j=JS0,j1+1,-1            dpi(imr, j, k) = dpi(imr, j, k) + xmass(imr, j) - xmass(1, j)
462        do i=1,IMR         END DO
463        if(abs(CRX(i,j)).GT.1.) then  
464              JS = j         DO j = j1, j2
465              go to 2222            DO i = 1, imr - 1
466        endif               ua(i, j) = 0.5*(crx(i,j)+crx(i+1,j))
467        enddo            END DO
468        enddo         END DO
469  C  
470  2222  continue         DO j = j1, j2
471        do j=JN0,j2-1            ua(imr, j) = 0.5*(crx(imr,j)+crx(1,j))
472        do i=1,IMR         END DO
473        if(abs(CRX(i,j)).GT.1.) then         ! cccccccccccccccccccccccccccccccccccccccccccccccccccccc
474              JN = j         ! Rajouts pour LMDZ.3.3
475              go to 2233         ! cccccccccccccccccccccccccccccccccccccccccccccccccccccc
476        endif         DO i = 1, imr
477        enddo            DO j = 1, jnp
478        enddo               va(i, j) = 0.
479  2233  continue            END DO
480  C         END DO
481        if(j1.ne.2) then           ! Enlarged polar cap.  
482        do i=1,IMR         DO i = 1, imr*(jmr-1)
483        DPI(i,  2,k) = 0.            va(i, 2) = 0.5*(cry(i,2)+cry(i,3))
484        DPI(i,JMR,k) = 0.         END DO
485        enddo  
486        endif         IF (j1==2) THEN
487  C            imh = imr/2
488  C ******* Compute horizontal mass fluxes ************            DO i = 1, imh
489  C               va(i, 1) = 0.5*(cry(i,2)-cry(i+imh,2))
490  C N-S component               va(i+imh, 1) = -va(i, 1)
491        do j=j1,j2+1               va(i, jnp) = 0.5*(cry(i,jnp)-cry(i+imh,jmr))
492        D5 = 0.5 * COSE(j)               va(i+imh, jnp) = -va(i, jnp)
493        do i=1,IMR            END DO
494        ymass(i,j) = CRY(i,j)*D5*(delp2(i,j,k) + delp2(i,j-1,k))            va(imr, 1) = va(1, 1)
495        enddo            va(imr, jnp) = va(1, jnp)
496        enddo         END IF
497  C  
498        do 95 j=j1,j2         ! ****6***0*********0*********0*********0*********0*********0**********72
499        DO 95 i=1,IMR         DO ic = 1, nc
500  95    DPI(i,j,k) = (ymass(i,j) - ymass(i,j+1)) * acosp(j)  
501  C            DO i = 1, imjm
502  C Poles               wk1(i, 1, 1) = 0.
503        sum1 = ymass(IMR,j1  )               wk1(i, 1, 2) = 0.
504        sum2 = ymass(IMR,J2+1)            END DO
505        do i=1,IMR-1  
506        sum1 = sum1 + ymass(i,j1  )            ! E-W advective cross term
507        sum2 = sum2 + ymass(i,J2+1)            DO j = j1, j2
508        enddo               IF (j>js .AND. j<jn) cycle
509  C  
510        sum1 = - sum1 * RCAP               DO i = 1, imr
511        sum2 =   sum2 * RCAP                  qtmp(i) = q(i, j, k, ic)
512        do i=1,IMR               END DO
513        DPI(i,  1,k) = sum1  
514        DPI(i,JNP,k) = sum2               DO i = -iml, 0
515        enddo                  qtmp(i) = q(imr+i, j, k, ic)
516  C                  qtmp(imr+1-i) = q(1-i, j, k, ic)
517  C E-W component               END DO
518  C  
519        do j=j1,j2               DO i = 1, imr
520        do i=2,IMR                  iu = ua(i, j)
521        PU(i,j) = 0.5 * (delp2(i,j,k) + delp2(i-1,j,k))                  ru = ua(i, j) - iu
522        enddo                  iiu = i - iu
523        enddo                  IF (ua(i,j)>=0.) THEN
524  C                     wk1(i, j, 1) = qtmp(iiu) + ru*(qtmp(iiu-1)-qtmp(iiu))
525        do j=j1,j2                  ELSE
526        PU(1,j) = 0.5 * (delp2(1,j,k) + delp2(IMR,j,k))                     wk1(i, j, 1) = qtmp(iiu) + ru*(qtmp(iiu)-qtmp(iiu+1))
527        enddo                  END IF
528  C                  wk1(i, j, 1) = wk1(i, j, 1) - qtmp(i)
529        do 110 j=j1,j2               END DO
530        DO 110 i=1,IMR            END DO
531  110   xmass(i,j) = PU(i,j)*CRX(i,j)  
532  C            IF (jn/=0) THEN
533        DO 120 j=j1,j2               DO j = js + 1, jn - 1
534        DO 120 i=1,IMR-1  
535  120   DPI(i,j,k) = DPI(i,j,k) + xmass(i,j) - xmass(i+1,j)                  DO i = 1, imr
536  C                     qtmp(i) = q(i, j, k, ic)
537        DO 130 j=j1,j2                  END DO
538  130   DPI(IMR,j,k) = DPI(IMR,j,k) + xmass(IMR,j) - xmass(1,j)  
539  C                  qtmp(0) = q(imr, j, k, ic)
540        DO j=j1,j2                  qtmp(imr+1) = q(1, j, k, ic)
541        do i=1,IMR-1  
542        UA(i,j) = 0.5 * (CRX(i,j)+CRX(i+1,j))                  DO i = 1, imr
543        enddo                     iu = i - ua(i, j)
544        enddo                     wk1(i, j, 1) = ua(i, j)*(qtmp(iu)-qtmp(iu+1))
545  C                  END DO
546        DO j=j1,j2               END DO
547        UA(imr,j) = 0.5 * (CRX(imr,j)+CRX(1,j))            END IF
548        enddo            ! ****6***0*********0*********0*********0*********0*********0**********72
549  ccccccccccccccccccccccccccccccccccccccccccccccccccccccc            ! Contribution from the N-S advection
550  c Rajouts pour LMDZ.3.3            DO i = 1, imr*(j2-j1+1)
551  ccccccccccccccccccccccccccccccccccccccccccccccccccccccc               jt = float(j1) - va(i, j1)
552        do i=1,IMR               wk1(i, j1, 2) = va(i, j1)*(q(i,jt,k,ic)-q(i,jt+1,k,ic))
553           do j=1,JNP            END DO
554               VA(i,j)=0.  
555           enddo            DO i = 1, imjm
556        enddo               wk1(i, 1, 1) = q(i, 1, k, ic) + 0.5*wk1(i, 1, 1)
557                 wk1(i, 1, 2) = q(i, 1, k, ic) + 0.5*wk1(i, 1, 2)
558        do i=1,imr*(JMR-1)            END DO
559        VA(i,2) = 0.5*(CRY(i,2)+CRY(i,3))  
560        enddo            IF (cross) THEN
561  C               ! Add cross terms in the vertical direction.
562        if(j1.eq.2) then               IF (iord>=2) THEN
563          IMH = IMR/2                  iad = 2
564        do i=1,IMH               ELSE
565        VA(i,      1) = 0.5*(CRY(i,2)-CRY(i+IMH,2))                  iad = 1
566        VA(i+IMH,  1) = -VA(i,1)               END IF
567        VA(i,    JNP) = 0.5*(CRY(i,JNP)-CRY(i+IMH,JMR))  
568        VA(i+IMH,JNP) = -VA(i,JNP)               IF (jord>=2) THEN
569        enddo                  jad = 2
570        VA(IMR,1)=VA(1,1)               ELSE
571        VA(IMR,JNP)=VA(1,JNP)                  jad = 1
572        endif               END IF
573  C               CALL xadv(imr, jnp, j1, j2, wk1(1,1,2), ua, js, jn, iml, dc2, iad)
574  C ****6***0*********0*********0*********0*********0*********0**********72               CALL yadv(imr, jnp, j1, j2, wk1(1,1,1), va, pv, w, jad)
575        do 1000 IC=1,NC               DO j = 1, jnp
576  C                  DO i = 1, imr
577        do i=1,IMJM                     q(i, j, k, ic) = q(i, j, k, ic) + dc2(i, j) + pv(i, j)
578        wk1(i,1,1) = 0.                  END DO
579        wk1(i,1,2) = 0.               END DO
580        enddo            END IF
581  C  
582  C E-W advective cross term            CALL xtp(imr, jnp, iml, j1, j2, jn, js, pu, dq(1,1,k,ic), wk1(1,1,2), &
583        do 250 j=J1,J2                 crx, fx1, xmass, iord)
584        if(J.GT.JS  .and. J.LT.JN) GO TO 250  
585  C            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)
587        qtmp(i) = q(i,j,k,IC)  
588        enddo         END DO
589  C      END DO
590        do i=-IML,0  
591        qtmp(i)       = q(IMR+i,j,k,IC)      ! ******* Compute vertical mass flux (same unit as PS) ***********
592        qtmp(IMR+1-i) = q(1-i,j,k,IC)  
593        enddo      ! 1st step: compute total column mass CONVERGENCE.
594  C  
595        DO 230 i=1,IMR      DO j = 1, jnp
596        iu = UA(i,j)         DO i = 1, imr
597        ru = UA(i,j) - iu            cry(i, j) = dpi(i, j, 1)
598        iiu = i-iu         END DO
599        if(UA(i,j).GE.0.) then      END DO
600        wk1(i,j,1) = qtmp(iiu)+ru*(qtmp(iiu-1)-qtmp(iiu))  
601        else      DO k = 2, nlay
602        wk1(i,j,1) = qtmp(iiu)+ru*(qtmp(iiu)-qtmp(iiu+1))         DO j = 1, jnp
603        endif            DO i = 1, imr
604        wk1(i,j,1) = wk1(i,j,1) - qtmp(i)               cry(i, j) = cry(i, j) + dpi(i, j, k)
605  230   continue            END DO
606  250   continue         END DO
607  C      END DO
608        if(JN.ne.0) then  
609        do j=JS+1,JN-1      DO j = 1, jnp
610  C         DO i = 1, imr
611        do i=1,IMR  
612        qtmp(i) = q(i,j,k,IC)            ! 2nd step: compute PS2 (PS at n+1) using the hydrostatic assumption.
613        enddo            ! Changes (increases) to surface pressure = total column mass
614  C            ! convergence
615        qtmp(0)     = q(IMR,J,k,IC)  
616        qtmp(IMR+1) = q(  1,J,k,IC)            ps2(i, j) = ps1(i, j) + cry(i, j)
617  C  
618        do i=1,imr            ! 3rd step: compute vertical mass flux from mass conservation
619        iu = i - UA(i,j)            ! principle.
620        wk1(i,j,1) = UA(i,j)*(qtmp(iu) - qtmp(iu+1))  
621        enddo            w(i, j, 1) = dpi(i, j, 1) - dbk(1)*cry(i, j)
622        enddo            w(i, j, nlay) = 0.
623        endif         END DO
624  C ****6***0*********0*********0*********0*********0*********0**********72      END DO
625  C Contribution from the N-S advection  
626        do i=1,imr*(j2-j1+1)      DO k = 2, nlay - 1
627        JT = float(J1) - VA(i,j1)         DO j = 1, jnp
628        wk1(i,j1,2) = VA(i,j1) * (q(i,jt,k,IC) - q(i,jt+1,k,IC))            DO i = 1, imr
629        enddo               w(i, j, k) = w(i, j, k-1) + dpi(i, j, k) - dbk(k)*cry(i, j)
630  C            END DO
631        do i=1,IMJM         END DO
632        wk1(i,1,1) = q(i,1,k,IC) + 0.5*wk1(i,1,1)      END DO
633        wk1(i,1,2) = q(i,1,k,IC) + 0.5*wk1(i,1,2)  
634        enddo      DO k = 1, nlay
635  C         DO j = 1, jnp
636          if(cross) then            DO i = 1, imr
637  C Add cross terms in the vertical direction.               delp2(i, j, k) = dap(k) + dbk(k)*ps2(i, j)
638          if(IORD .GE. 2) then            END DO
639                  iad = 2         END DO
640          else      END DO
641                  iad = 1  
642          endif      krd = max(3, kord)
643  C      DO ic = 1, nc
644          if(JORD .GE. 2) then  
645                  jad = 2         ! ****6***0*********0*********0*********0*********0*********0**********72
646          else  
647                  jad = 1         CALL fzppm(imr, jnp, nlay, j1, dq(1,1,1,ic), w, q(1,1,1,ic), wk1, dpi, &
648          endif              dc2, crx, cry, pu, pv, xmass, ymass, delp1, krd)
649        call xadv(IMR,JNP,j1,j2,wk1(1,1,2),UA,JS,JN,IML,DC2,iad)  
650        call yadv(IMR,JNP,j1,j2,wk1(1,1,1),VA,PV,W,jad)  
651        do j=1,JNP         IF (fill) CALL qckxyz(dq(1,1,1,ic), dc2, imr, jnp, nlay, j1, j2, cosp, &
652        do i=1,IMR              acosp, .FALSE., ic, nstep)
653        q(i,j,k,IC) = q(i,j,k,IC) + DC2(i,j) + PV(i,j)  
654        enddo         ! Recover tracer mixing ratio from "density" using predicted
655        enddo         ! "air density" (pressure thickness) at time-level n+1
656        endif  
657  C         DO k = 1, nlay
658        call xtp(IMR,JNP,IML,j1,j2,JN,JS,PU,DQ(1,1,k,IC),wk1(1,1,2)            DO j = 1, jnp
659       &        ,CRX,fx1,xmass,IORD)               DO i = 1, imr
660                    q(i, j, k, ic) = dq(i, j, k, ic)/delp2(i, j, k)
661        call ytp(IMR,JNP,j1,j2,acosp,RCAP,DQ(1,1,k,IC),wk1(1,1,1),CRY,               END DO
662       &  DC2,ymass,WK1(1,1,3),wk1(1,1,4),WK1(1,1,5),WK1(1,1,6),JORD)            END DO
663  C         END DO
664  1000  continue  
665  1500  continue         IF (j1/=2) THEN
666  C            DO k = 1, nlay
667  C ******* Compute vertical mass flux (same unit as PS) ***********               DO i = 1, imr
668  C                  ! j=1 c'est le p\^ole Sud, j=JNP c'est le p\^ole Nord
669  C 1st step: compute total column mass CONVERGENCE.                  q(i, 2, k, ic) = q(i, 1, k, ic)
670  C                  q(i, jmr, k, ic) = q(i, jmp, k, ic)
671        do 320 j=1,JNP               END DO
672        do 320 i=1,IMR            END DO
673  320   CRY(i,j) = DPI(i,j,1)         END IF
674  C      END DO
675        do 330 k=2,NLAY  
676        do 330 j=1,JNP      IF (j1/=2) THEN
677        do 330 i=1,IMR         DO k = 1, nlay
678        CRY(i,j)  = CRY(i,j) + DPI(i,j,k)            DO i = 1, imr
679  330   continue               w(i, 2, k) = w(i, 1, k)
680  C               w(i, jmr, k) = w(i, jnp, k)
681        do 360 j=1,JNP            END DO
682        do 360 i=1,IMR         END DO
683  C      END IF
684  C 2nd step: compute PS2 (PS at n+1) using the hydrostatic assumption.  
685  C Changes (increases) to surface pressure = total column mass convergence    END SUBROUTINE ppm3d
686  C  
687        PS2(i,j)  = PS1(i,j) + CRY(i,j)  end module ppm3d_m
 C  
 C 3rd step: compute vertical mass flux from mass conservation principle.  
 C  
       W(i,j,1) = DPI(i,j,1) - DBK(1)*CRY(i,j)  
       W(i,j,NLAY) = 0.  
 360   continue  
 C  
       do 370 k=2,NLAY-1  
       do 370 j=1,JNP  
       do 370 i=1,IMR  
       W(i,j,k) = W(i,j,k-1) + DPI(i,j,k) - DBK(k)*CRY(i,j)  
 370   continue  
 C  
       DO 380 k=1,NLAY  
       DO 380 j=1,JNP  
       DO 380 i=1,IMR  
       delp2(i,j,k) = DAP(k) + DBK(k)*PS2(i,j)  
 380   continue  
 C  
         KRD = max(3, KORD)  
       do 4000 IC=1,NC  
 C  
 C****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)  
 C  
       
       if(fill) call qckxyz(DQ(1,1,1,IC),DC2,IMR,JNP,NLAY,j1,j2,  
      &                     cosp,acosp,.false.,IC,NSTEP)  
 C  
 C Recover tracer mixing ratio from "density" using predicted  
 C "air density" (pressure thickness) at time-level n+1  
 C  
       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)  
       enddo  
       enddo  
       enddo  
 C      
       if(j1.ne.2) then  
       DO 400 k=1,NLAY  
       DO 400 I=1,IMR  
 c     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)  
 400   CONTINUE  
       endif  
 4000  continue  
 C  
       if(j1.ne.2) then  
       DO 5000 k=1,NLAY  
       DO 5000 i=1,IMR  
       W(i,  2,k) = W(i,  1,k)  
       W(i,JMR,k) = W(i,JNP,k)  
 5000  continue  
       endif  
 C  
       RETURN  
       END  
 C  
 C****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)  
 C Assuming JNP >= NLAY  
       real AR(IMR,*),AL(IMR,*),A6(IMR,*),flux(IMR,*),wk2(IMR,*),  
      &     wz2(IMR,*)  
 C  
       JMR = JNP - 1  
       IMJM = IMR*JNP  
       NLAYM1 = NLAY - 1  
 C  
       LMT = KORD - 3  
 C  
 C ****6***0*********0*********0*********0*********0*********0**********72  
 C Compute DC for PPM  
 C ****6***0*********0*********0*********0*********0*********0**********72  
 C  
       do 1000 k=1,NLAYM1  
       do 1000 i=1,IMJM  
       DQDT(i,1,k) = P(i,1,k+1) - P(i,1,k)  
 1000  continue  
 C  
       DO 1220 k=2,NLAYM1  
       DO 1220 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)    
 1220  CONTINUE  
       
 C      
 C ****6***0*********0*********0*********0*********0*********0**********72  
 C Loop over latitudes  (to save memory)  
 C ****6***0*********0*********0*********0*********0*********0**********72  
 C  
       DO 2000 j=1,JNP  
       if((j.eq.2 .or. j.eq.JMR) .and. j1.ne.2) goto 2000  
 C  
       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  
 C  
 C****6***0*********0*********0*********0*********0*********0**********72  
 C Compute first guesses at cell interfaces  
 C First guesses are required to be continuous.  
 C ****6***0*********0*********0*********0*********0*********0**********72  
 C  
 C three-cell parabolic subgrid distribution at model top  
 C two-cell parabolic with zero gradient subgrid distribution  
 C at the surface.  
 C  
 C First guess top edge value  
       DO 10 i=1,IMR  
 C three-cell PPM  
 C 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)  
 C  
 C 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  
 10    continue  
 C  
 C Bottom  
       DO 15 i=1,IMR  
 C 2-cell PPM with zero gradient right at the surface  
 C  
       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)  
 15    continue  
       
 C  
 C****6***0*********0*********0*********0*********0*********0**********72  
 C 4th order interpolation in the interior.  
 C****6***0*********0*********0*********0*********0*********0**********72  
 C  
       DO 14 k=3,NLAYM1  
       DO 12 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)  )  
 12    CONTINUE  
 14    continue  
 C  
       do 20 i=1,IMR*NLAYM1  
       AR(i,1) = AL(i,2)  
 20    continue  
 C  
       do 30 i=1,IMR*NLAY  
       A6(i,1) = 3.*(wk1(i,1)+wk1(i,1) - (AL(i,1)+AR(i,1)))  
 30    continue  
 C  
 C****6***0*********0*********0*********0*********0*********0**********72  
 C 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)  
 C  
 C 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)  
 C  
 C****6***0*********0*********0*********0*********0*********0**********72  
 C  
       DO 140 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  
         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))  
       endif  
 140   continue  
 C  
       DO 250 i=1,IMR*NLAYM1  
       flux(i,2) = wz2(i,1) * flux(i,2)  
 250   continue  
 C  
       do 350 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)  
 350   continue  
 C  
       do 360 k=2,NLAYM1  
       do 360 i=1,IMR  
 360   DQ(i,j,k) = DQ(i,j,k) + flux(i,k) - flux(i,k+1)  
 2000  continue  
       return  
       end  
 C  
       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)  
 C  
       IMP = IMR + 1  
 C  
 C van Leer at high latitudes  
       jvan = max(1,JNP/18)  
       j1vl = j1+jvan  
       j2vl = j2-jvan  
 C  
       do 1310 j=j1,j2  
 C  
       do i=1,IMR  
       qtmp(i) = q(i,j)  
       enddo  
 C  
       if(j.ge.JN .or. j.le.JS) goto 2222  
 C ************* Eulerian **********  
 C  
       qtmp(0)     = q(IMR,J)  
       qtmp(-1)    = q(IMR-1,J)  
       qtmp(IMP)   = q(1,J)  
       qtmp(IMP+1) = q(2,J)  
 C  
       IF(IORD.eq.1 .or. j.eq.j1. or. j.eq.j2) THEN  
       DO 1406 i=1,IMR  
       iu = float(i) - uc(i,j)  
 1406  fx1(i) = qtmp(iu)  
       ELSE  
       call xmist(IMR,IML,Qtmp,DC)  
       DC(0) = DC(IMR)  
 C  
       if(IORD.eq.2 .or. j.le.j1vl .or. j.ge.j2vl) then  
       DO 1408 i=1,IMR  
       iu = float(i) - uc(i,j)  
 1408  fx1(i) = qtmp(iu) + DC(iu)*(sign(1.,uc(i,j))-uc(i,j))  
       else  
       call fxppm(IMR,IML,UC(1,j),Qtmp,DC,fx1,IORD)  
       endif  
 C  
       ENDIF  
 C  
       DO 1506 i=1,IMR  
 1506  fx1(i) = fx1(i)*xmass(i,j)  
 C  
       goto 1309  
 C  
 C ***** Conservative (flux-form) Semi-Lagrangian transport *****  
 C  
 2222  continue  
 C  
       do i=-IML,0  
       qtmp(i)     = q(IMR+i,j)  
       qtmp(IMP-i) = q(1-i,j)  
       enddo  
 C  
       IF(IORD.eq.1 .or. j.eq.j1. or. j.eq.j2) THEN  
       DO 1306 i=1,IMR  
       itmp = INT(uc(i,j))  
       ISAVE(i) = i - itmp  
       iu = i - uc(i,j)  
 1306  fx1(i) = (uc(i,j) - itmp)*qtmp(iu)  
       ELSE  
       call xmist(IMR,IML,Qtmp,DC)  
 C  
       do i=-IML,0  
       DC(i)     = DC(IMR+i)  
       DC(IMP-i) = DC(1-i)  
       enddo  
 C  
       DO 1307 i=1,IMR  
       itmp = INT(uc(i,j))  
       rut  = uc(i,j) - itmp  
       ISAVE(i) = i - itmp  
       iu = i - uc(i,j)  
 1307  fx1(i) = rut*(qtmp(iu) + DC(iu)*(sign(1.,rut) - rut))  
       ENDIF  
 C  
       do 1308 i=1,IMR  
       IF(uc(i,j).GT.1.) then  
 CDIR$ 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  
 CDIR$ VECTOR  
       endif  
 1308  continue  
       do i=1,IMR  
       fx1(i) = PU(i,j)*fx1(i)  
       enddo  
 C  
 C ***************************************  
 C  
 1309  fx1(IMP) = fx1(1)  
       DO 1215 i=1,IMR  
 1215  DQ(i,j) =  DQ(i,j) + fx1(i)-fx1(i+1)  
 C  
 C ***************************************  
 C  
 1310  continue  
       return  
       end  
 C  
       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  
 c      logical first  
 c      data first /.true./  
 c      SAVE LMT  
 c      if(first) then  
 C  
 C correction calcul de LMT a chaque passage pour pouvoir choisir  
 c plusieurs schemas PPM pour differents traceurs  
 c      IF (IORD.LE.0) then  
 c            if(IMR.GE.144) then  
 c                  LMT = 0  
 c            elseif(IMR.GE.72) then  
 c                  LMT = 1  
 c            else  
 c                  LMT = 2  
 c            endif  
 c      else  
 c            LMT = IORD - 3  
 c      endif  
 C  
       LMT = IORD - 3  
   
       DO 10 i=1,IMR  
 10    AL(i) = 0.5*(p(i-1)+p(i)) + (DC(i-1) - DC(i))*R3  
 C  
       do 20 i=1,IMR-1  
 20    AR(i) = AL(i+1)  
       AR(IMR) = AL(1)  
 C  
       do 30 i=1,IMR  
 30    A6(i) = 3.*(p(i)+p(i)  - (AL(i)+AR(i)))  
 C  
       if(LMT.LE.2) call lmtppm(DC(1),A6(1),AR(1),AL(1),P(1),IMR,LMT)  
 C  
       AL(0) = AL(IMR)  
       AR(0) = AR(IMR)  
       A6(0) = A6(IMR)  
 C  
       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  
 C  
       subroutine xmist(IMR,IML,P,DC)  
       parameter( R24 = 1./24.)  
       dimension P(-IML:IMR+1+IML),DC(-IML:IMR+1+IML)  
 C  
       do 10  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))  
 10    DC(i) = sign(min(abs(tmp),Pmax,Pmin), tmp)  
       return  
       end  
 C  
       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)  
 C Work array  
       DIMENSION fx(IMR,JNP),AR(IMR,JNP),AL(IMR,JNP),A6(IMR,JNP)  
 C  
       JMR = JNP - 1  
       len = IMR*(J2-J1+2)  
 C  
       if(JORD.eq.1) then  
       DO 1000 i=1,len  
       JT = float(J1) - VC(i,J1)  
 1000  fx(i,j1) = p(i,JT)  
       else  
     
       call ymist(IMR,JNP,j1,P,DC2,4)  
 C  
       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 1200 i=1,len  
       JT = float(J1) - VC(i,J1)  
 1200  fx(i,j1) = p(i,JT) + (sign(1.,VC(i,j1))-VC(i,j1))*DC2(i,JT)  
       endif  
       endif  
 C  
       DO 1300 i=1,len  
 1300  fx(i,j1) = fx(i,j1)*ymass(i,j1)  
 C  
       DO 1400 j=j1,j2  
       DO 1400 i=1,IMR  
 1400  DQ(i,j) = DQ(i,j) + (fx(i,j) - fx(i,j+1)) * acosp(j)  
 C  
 C 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  
 C  
       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  
 C  
       if(j1.ne.2) then  
       do i=1,IMR  
       DQ(i,  2) = sum1  
       DQ(i,JMR) = sum2  
       enddo  
       endif  
 C  
       return  
       end  
 C  
       subroutine  ymist(IMR,JNP,j1,P,DC,ID)  
       parameter ( R24 = 1./24. )  
       dimension P(IMR,JNP),DC(IMR,JNP)  
 C  
       IMH = IMR / 2  
       JMR = JNP - 1  
       IJM3 = IMR*(JMR-3)  
 C  
       IF(ID.EQ.2) THEN  
       do 10 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)  
 10    CONTINUE  
       ELSE  
       do 12 i=1,IMH  
 C 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)  
 C 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)  
 12    CONTINUE  
       do 14 i=IMH+1,IMR  
 C 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)  
 C 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)  
 14    CONTINUE  
 C  
       do 15 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)  
 15    CONTINUE  
       ENDIF  
 C  
       if(j1.ne.2) then  
       do i=1,IMR  
       DC(i,1) = 0.  
       DC(i,JNP) = 0.  
       enddo  
       else  
 C Determine slopes in polar caps for scalars!  
 C  
       do 13 i=1,IMH  
 C 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)  
 C 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)  
 13    continue  
 C  
       do 25 i=imh+1,IMR  
       DC(i,  1) =  - DC(i-imh,  1)  
       DC(i,JNP) =  - DC(i-imh,JNP)  
 25    continue  
       endif  
       return  
       end  
 C  
       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,*)  
 C Local work arrays.  
       real AR(IMR,JNP),AL(IMR,JNP),A6(IMR,JNP)  
       integer LMT  
 c      logical first  
 C      data first /.true./  
 C      SAVE LMT  
 C  
       IMH = IMR / 2  
       JMR = JNP - 1  
       j11 = j1-1  
       IMJM1 = IMR*(J2-J1+2)  
       len   = IMR*(J2-J1+3)  
 C      if(first) then  
 C      IF(JORD.LE.0) then  
 C            if(JMR.GE.90) then  
 C                  LMT = 0  
 C            elseif(JMR.GE.45) then  
 C                  LMT = 1  
 C            else  
 C                  LMT = 2  
 C            endif  
 C      else  
 C            LMT = JORD - 3  
 C      endif  
 C  
 C      first = .false.  
 C      endif  
 C      
 c modifs pour pouvoir choisir plusieurs schemas PPM  
       LMT = JORD - 3        
 C  
       DO 10 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)  
 10    CONTINUE  
 C  
 CPoles:  
 C  
       DO i=1,IMH  
       AL(i,1) = AL(i+IMH,2)  
       AL(i+IMH,1) = AL(i,2)  
 C  
       AR(i,JNP) = AR(i+IMH,JMR)  
       AR(i+IMH,JNP) = AR(i,JMR)  
       ENDDO  
   
 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc  
 c   Rajout pour LMDZ.3.3  
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc  
       AR(IMR,1)=AL(1,1)  
       AR(IMR,JNP)=AL(1,JNP)  
 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc  
         
             
       do 30 i=1,len  
 30    A6(i,j11) = 3.*(p(i,j11)+p(i,j11)  - (AL(i,j11)+AR(i,j11)))  
 C  
       if(LMT.le.2) call lmtppm(DC(1,j11),A6(1,j11),AR(1,j11)  
      &                       ,AL(1,j11),P(1,j11),len,LMT)  
 C  
       
       DO 140 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  
 140   continue  
       return  
       end  
 C  
         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)  
 C  
         JMR = JNP-1  
         IMH = IMR/2  
         do j=1,JNP  
         do i=1,IMR  
         wk(i,j) = p(i,j)  
         enddo  
         enddo  
 C 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  
   
       IF(IAD.eq.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)  
       enddo  
       enddo  
   
       ELSEIF(IAD.eq.1) then  
         do j=j1-1,j2+1  
       do i=1,imr  
       JP = float(j)-VA(i,j)  
       ady(i,j) = VA(i,j)*(wk(i,jp)-wk(i,jp+1))  
       enddo  
       enddo  
       ENDIF  
 C  
         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  
 C  
       do i=1,imr  
       ady(i,  2) =  sum1  
       ady(i,JMR) =  sum2  
       ady(i,  1) =  sum1  
       ady(i,JNP) =  sum2  
       enddo  
         else  
 C 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  
 C  
       do i=1,imr  
       ady(i,  1) =  sum1  
       ady(i,JNP) =  sum2  
       enddo  
         endif  
 C  
         return  
         end  
 C  
         subroutine xadv(IMR,JNP,j1,j2,p,UA,JS,JN,IML,adx,IAD)  
         REAL p(IMR,JNP),adx(IMR,JNP),qtmp(-IMR:IMR+IMR),UA(IMR,JNP)  
 C  
         JMR = JNP-1  
       do 1309 j=j1,j2  
       if(J.GT.JS  .and. J.LT.JN) GO TO 1309  
 C  
       do i=1,IMR  
       qtmp(i) = p(i,j)  
       enddo  
 C  
       do i=-IML,0  
       qtmp(i)       = p(IMR+i,j)  
       qtmp(IMR+1-i) = p(1-i,j)  
       enddo  
 C  
       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  
 C  
       do i=1,IMR  
       adx(i,j) = adx(i,j) - p(i,j)  
       enddo  
 1309  continue  
 C  
 C Eulerian upwind  
 C  
       do j=JS+1,JN-1  
 C  
       do i=1,IMR  
       qtmp(i) = p(i,j)  
       enddo  
 C  
       qtmp(0)     = p(IMR,J)  
       qtmp(IMR+1) = p(1,J)  
 C  
       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  
 C 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  
 C  
         if(j1.ne.2) then  
       do i=1,IMR  
       adx(i,  2) = 0.  
       adx(i,JMR) = 0.  
       enddo  
         endif  
 C 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  
 C  
       subroutine lmtppm(DC,A6,AR,AL,P,IM,LMT)  
 C  
 C A6 =  CURVATURE OF THE TEST PARABOLA  
 C AR =  RIGHT EDGE VALUE OF THE TEST PARABOLA  
 C AL =  LEFT  EDGE VALUE OF THE TEST PARABOLA  
 C DC =  0.5 * MISMATCH  
 C P  =  CELL-AVERAGED VALUE  
 C IM =  VECTOR LENGTH  
 C  
 C OPTIONS:  
 C  
 C LMT = 0: FULL MONOTONICITY  
 C LMT = 1: SEMI-MONOTONIC CONSTRAINT (NO UNDERSHOOTS)  
 C LMT = 2: POSITIVE-DEFINITE CONSTRAINT  
 C  
       parameter ( R12 = 1./12. )  
       dimension A6(IM),AR(IM),AL(IM),P(IM),DC(IM)  
 C  
       if(LMT.eq.0) then  
 C Full constraint  
       do 100 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  
 100   continue  
       elseif(LMT.eq.1) then  
 C Semi-monotonic constraint  
       do 150 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  
       elseif(LMT.eq.2) then  
       do 250 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  
       endif  
       return  
       end  
 C  
       subroutine A2C(U,V,IMR,JMR,j1,j2,CRX,CRY,dtdx5,DTDY5)  
       dimension U(IMR,*),V(IMR,*),CRX(IMR,*),CRY(IMR,*),DTDX5(*)  
 C  
       do 35 j=j1,j2  
       do 35 i=2,IMR  
 35    CRX(i,J) = dtdx5(j)*(U(i,j)+U(i-1,j))  
 C  
       do 45 j=j1,j2  
 45    CRX(1,J) = dtdx5(j)*(U(1,j)+U(IMR,j))  
 C  
       do 55 i=1,IMR*JMR  
 55    CRY(i,2) = DTDY5*(V(i,2)+V(i,1))  
       return  
       end  
 C  
       subroutine cosa(cosp,cose,JNP,PI,DP)  
       dimension cosp(*),cose(*)  
       JMR = JNP-1  
       do 55 j=2,JNP  
         ph5  =  -0.5*PI + (FLOAT(J-1)-0.5)*DP  
 55      cose(j) = cos(ph5)  
 C  
       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  
 C cell edge at equator.  
        cose(JEQ+1) =  1.  
       do j=JNP, JEQ+2, -1  
        cose(j) =  cose(JNP+2-j)  
        enddo  
       endif  
 C  
       do 66 j=2,JMR  
 66    cosp(j) = 0.5*(cose(j)+cose(j+1))  
       cosp(1) = 0.  
       cosp(JNP) = 0.  
       return  
       end  
 C  
       subroutine cosc(cosp,cose,JNP,PI,DP)  
       dimension cosp(*),cose(*)  
 C  
       phi = -0.5*PI  
       do 55 j=2,JNP-1  
       phi  =  phi + DP  
 55    cosp(j) = cos(phi)  
         cosp(  1) = 0.  
         cosp(JNP) = 0.  
 C  
       do 66 j=2,JNP  
         cose(j) = 0.5*(cosp(j)+cosp(j-1))  
 66    CONTINUE  
 C  
       do 77 j=2,JNP-1  
        cosp(j) = 0.5*(cose(j)+cose(j+1))  
 77    CONTINUE  
       return  
       end  
 C  
       SUBROUTINE qckxyz (Q,qtmp,IMR,JNP,NLAY,j1,j2,cosp,acosp,  
      &                   cross,IC,NSTEP)  
 C  
       parameter( tiny = 1.E-60 )  
       DIMENSION Q(IMR,JNP,NLAY),qtmp(IMR,JNP),cosp(*),acosp(*)  
       logical cross  
 C  
       NLAYM1 = NLAY-1  
       len = IMR*(j2-j1+1)  
       ip = 0  
 C  
 C 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  
 C  
       if(cross) then  
       call filcr(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,icr,tiny)  
       endif  
       if(icr.eq.0) goto 50  
 C  
 C 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  
 C  
 50    continue  
       DO 225 L = 2,NLAYM1  
       icr = 1  
 C  
       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  
 C  
       do i=1,len  
       IF( Q(I,j1,L).LT.0.) THEN  
 C  
       ip = ip + 1  
 C 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  
 C Below  
           Q(I,j1,L+1) = Q(I,j1,L+1) + Q(I,j1,L)  
           Q(I,j1,L)   = 0.  
       ENDIF  
       ENDDO  
 225   CONTINUE  
 C  
 C BOTTOM LAYER  
       sum = 0.  
       L = NLAY  
 C  
       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  
 C  
       call filcr(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,icr,tiny)  
       if(icr.eq.0) goto 911  
 C  
       DO  I=1,len  
       IF( Q(I,j1,L).LT.0.) THEN  
       ip = ip + 1  
 c  
 C From above  
 C  
           qup = Q(I,j1,NLAYM1)  
           qly = -Q(I,j1,L)  
           dup = min(qly,qup)  
           Q(I,j1,NLAYM1) = qup - dup  
 C From "below" the surface.  
           sum = sum + qly-dup  
           Q(I,j1,L) = 0.  
        ENDIF  
       ENDDO  
 C  
 911   continue  
 C  
       if(ip.gt.IMR) then  
       write(6,*) 'IC=',IC,' STEP=',NSTEP,  
      &           ' Vertical filling pts=',ip  
       endif  
 C  
       if(sum.gt.1.e-25) then  
       write(6,*) IC,NSTEP,' Mass source from the ground=',sum  
       endif  
       RETURN  
       END  
 C  
       subroutine filcr(q,IMR,JNP,j1,j2,cosp,acosp,icr,tiny)  
       dimension q(IMR,*),cosp(*),acosp(*)  
       icr = 0  
       do 65 j=j1+1,j2-1  
       DO 50 i=1,IMR-1  
       IF(q(i,j).LT.0.) THEN  
       icr =  1  
       dq  = - q(i,j)*cosp(j)  
 C 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  
 C 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  
       if(icr.eq.0 .and. q(IMR,j).ge.0.) goto 65  
       DO 55 i=2,IMR  
       IF(q(i,j).LT.0.) THEN  
       icr =  1  
       dq  = - q(i,j)*cosp(j)  
 C 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  
 C 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  
 55    continue  
 C *****************************************  
 C i=1  
       i=1  
       IF(q(i,j).LT.0.) THEN  
       icr =  1  
       dq  = - q(i,j)*cosp(j)  
 C 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  
 C 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  
 C *****************************************  
 C i=IMR  
       i=IMR  
       IF(q(i,j).LT.0.) THEN  
       icr =  1  
       dq  = - q(i,j)*cosp(j)  
 C 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  
 C 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  
 C *****************************************  
 65    continue  
 C  
       do i=1,IMR  
       if(q(i,j1).lt.0. .or. q(i,j2).lt.0.) then  
       icr = 1  
       goto 80  
       endif  
       enddo  
 C  
 80    continue  
 C  
       if(q(1,1).lt.0. .or. q(1,jnp).lt.0.) then  
       icr = 1  
       endif  
 C  
       return  
       end  
 C  
       subroutine filns(q,IMR,JNP,j1,j2,cosp,acosp,ipy,tiny)  
       dimension q(IMR,*),cosp(*),acosp(*)  
 c      logical first  
 c      data first /.true./  
 c      save cap1  
 C  
 c      if(first) then  
       DP = 4.*ATAN(1.)/float(JNP-1)  
       CAP1 = IMR*(1.-COS((j1-1.5)*DP))/DP  
 c      first = .false.  
 c      endif  
 C  
       ipy = 0  
       do 55 j=j1+1,j2-1  
       DO 55 i=1,IMR  
       IF(q(i,j).LT.0.) THEN  
       ipy =  1  
       dq  = - q(i,j)*cosp(j)  
 C 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  
 C 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  
 55    continue  
 C  
       do i=1,imr  
       IF(q(i,j1).LT.0.) THEN  
       ipy =  1  
       dq  = - q(i,j1)*cosp(j1)  
 C 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  
 C  
       j = j2  
       do i=1,imr  
       IF(q(i,j).LT.0.) THEN  
       ipy =  1  
       dq  = - q(i,j)*cosp(j)  
 C 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  
 C  
 C Check Poles.  
       if(q(1,1).lt.0.) then  
       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).lt.0.) ipy = 1  
       enddo  
       endif  
 C  
       if(q(1,JNP).lt.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).lt.0.) ipy = 1  
       enddo  
       endif  
 C  
       return  
       end  
 C  
       subroutine filew(q,qtmp,IMR,JNP,j1,j2,ipx,tiny)  
       dimension q(IMR,*),qtmp(JNP,IMR)  
 C  
       ipx = 0  
 C Copy & swap direction for vectorization.  
       do 25 i=1,imr  
       do 25 j=j1,j2  
 25    qtmp(j,i) = q(i,j)  
 C  
       do 55 i=2,imr-1  
       do 55 j=j1,j2  
       if(qtmp(j,i).lt.0.) then  
       ipx =  1  
 c 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  
 c 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  
 55    continue  
 c  
       i=1  
       do 65 j=j1,j2  
       if(qtmp(j,i).lt.0.) then  
       ipx =  1  
 c 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  
 c east  
       d0 = max(0.,qtmp(j,i+1))  
       d2 = min(-qtmp(j,i),d0)  
       qtmp(j,i+1) = qtmp(j,i+1) - d2  
 c  
       qtmp(j,i) = qtmp(j,i) + d2 + tiny  
       endif  
 65    continue  
       i=IMR  
       do 75 j=j1,j2  
       if(qtmp(j,i).lt.0.) then  
       ipx =  1  
 c 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  
 c east  
       d0 = max(0.,qtmp(j,1))  
       d2 = min(-qtmp(j,i),d0)  
       qtmp(j,1) = qtmp(j,1) - d2  
 c  
       qtmp(j,i) = qtmp(j,i) + d2 + tiny  
       endif  
 75    continue  
 C  
       if(ipx.ne.0) then  
       do 85 j=j1,j2  
       do 85 i=1,imr  
 85    q(i,j) = qtmp(j,i)  
       else  
 C  
 C Poles.  
       if(q(1,1).lt.0. or. q(1,JNP).lt.0.) ipx = 1  
       endif  
       return  
       end  
Legend:
Removed from v.76  
changed lines
  Added in v.169
  ViewVC Help
Powered by ViewVC 1.1.21