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