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

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

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

trunk/libf/dyn3d/ppm3d.f revision 3 by guez, Wed Feb 27 13:16:39 2008 UTC trunk/dyn3d/ppm3d.f90 revision 81 by guez, Wed Mar 5 14:38:41 2014 UTC
# Line 1  Line 1 
1  !  
2  ! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/ppm3d.F,v 1.1.1.1 2004/05/19 12:53:07 lmdzadmin Exp $  ! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/ppm3d.F,v 1.1.1.1 2004/05/19
3  !  ! 12:53:07 lmdzadmin Exp $
4    
5  cFrom lin@explorer.gsfc.nasa.gov Wed Apr 15 17:44:44 1998  
6  cDate: Wed, 15 Apr 1998 11:37:03 -0400  ! From lin@explorer.gsfc.nasa.gov Wed Apr 15 17:44:44 1998
7  cFrom: lin@explorer.gsfc.nasa.gov  ! Date: Wed, 15 Apr 1998 11:37:03 -0400
8  cTo: Frederic.Hourdin@lmd.jussieu.fr  ! From: lin@explorer.gsfc.nasa.gov
9  cSubject: 3D transport module of the GSFC CTM and GEOS GCM  ! To: Frederic.Hourdin@lmd.jussieu.fr
10    ! Subject: 3D transport module of the GSFC CTM and GEOS GCM
11    
12  cThis code is sent to you by S-J Lin, DAO, NASA-GSFC  
13    ! This code is sent to you by S-J Lin, DAO, NASA-GSFC
14  cNote: this version is intended for machines like CRAY  
15  C-90. No multitasking directives implemented.  ! Note: this version is intended for machines like CRAY
16    ! -90. No multitasking directives implemented.
17          
18  C ********************************************************************  
19  C  ! ********************************************************************
20  C TransPort Core for Goddard Chemistry Transport Model (G-CTM), Goddard  
21  C Earth Observing System General Circulation Model (GEOS-GCM), and Data  ! TransPort Core for Goddard Chemistry Transport Model (G-CTM), Goddard
22  C Assimilation System (GEOS-DAS).  ! Earth Observing System General Circulation Model (GEOS-GCM), and Data
23  C  ! Assimilation System (GEOS-DAS).
24  C ********************************************************************  
25  C  ! ********************************************************************
26  C Purpose: given horizontal winds on  a hybrid sigma-p surfaces,  
27  C          one call to tpcore updates the 3-D mixing ratio  ! Purpose: given horizontal winds on  a hybrid sigma-p surfaces,
28  C          fields one time step (NDT). [vertical mass flux is computed  ! one call to tpcore updates the 3-D mixing ratio
29  C          internally consistent with the discretized hydrostatic mass  ! fields one time step (NDT). [vertical mass flux is computed
30  C          continuity equation of the C-Grid GEOS-GCM (for IGD=1)].  ! internally consistent with the discretized hydrostatic mass
31  C  ! continuity equation of the C-Grid GEOS-GCM (for IGD=1)].
32  C Schemes: Multi-dimensional Flux Form Semi-Lagrangian (FFSL) scheme based  
33  C          on the van Leer or PPM.  ! Schemes: Multi-dimensional Flux Form Semi-Lagrangian (FFSL) scheme based
34  C          (see Lin and Rood 1996).  ! on the van Leer or PPM.
35  C Version 4.5  ! (see Lin and Rood 1996).
36  C Last modified: Dec. 5, 1996  ! Version 4.5
37  C Major changes from version 4.0: a more general vertical hybrid sigma-  ! Last modified: Dec. 5, 1996
38  C pressure coordinate.  ! Major changes from version 4.0: a more general vertical hybrid sigma-
39  C Subroutines modified: xtp, ytp, fzppm, qckxyz  ! pressure coordinate.
40  C Subroutines deleted: vanz  ! Subroutines modified: xtp, ytp, fzppm, qckxyz
41  C  ! Subroutines deleted: vanz
42  C Author: Shian-Jiann Lin  
43  C mail address:  ! Author: Shian-Jiann Lin
44  C                 Shian-Jiann Lin*  ! mail address:
45  C                 Code 910.3, NASA/GSFC, Greenbelt, MD 20771  ! Shian-Jiann Lin*
46  C                 Phone: 301-286-9540  ! Code 910.3, NASA/GSFC, Greenbelt, MD 20771
47  C                 E-mail: lin@dao.gsfc.nasa.gov  ! Phone: 301-286-9540
48  C  ! E-mail: lin@dao.gsfc.nasa.gov
49  C *affiliation:  
50  C                 Joint Center for Earth Systems Technology  ! *affiliation:
51  C                 The University of Maryland Baltimore County  ! Joint Center for Earth Systems Technology
52  C                 NASA - Goddard Space Flight Center  ! The University of Maryland Baltimore County
53  C References:  ! NASA - Goddard Space Flight Center
54  C  ! References:
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.  ! 1. Lin, S.-J., and R. B. Rood, 1996: Multidimensional flux form semi-
57  C  ! Lagrangian transport schemes. Mon. Wea. Rev., 124, 2046-2070.
58  C 2. Lin, S.-J., W. C. Chao, Y. C. Sud, and G. K. Walker, 1994: A class of  
59  C    the van Leer-type transport schemes and its applications to the moist-  ! 2. Lin, S.-J., W. C. Chao, Y. C. Sud, and G. K. Walker, 1994: A class of
60  C    ure transport in a General Circulation Model. Mon. Wea. Rev., 122,  ! the van Leer-type transport schemes and its applications to the moist-
61  C    1575-1593.  ! ure transport in a General Circulation Model. Mon. Wea. Rev., 122,
62  C  ! 1575-1593.
63  C ****6***0*********0*********0*********0*********0*********0**********72  
64  C  ! ****6***0*********0*********0*********0*********0*********0**********72
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)  SUBROUTINE ppm3d(igd, q, ps1, ps2, u, v, w, ndt, iord, jord, kord, nc, imr, &
67        jnp, j1, nlay, ap, bp, pt, ae, fill, dum, umax)
68  c      implicit none  
69      ! implicit none
70  c     rajout de déclarations  
71  c      integer Jmax,kmax,ndt0,nstep,k,j,i,ic,l,js,jn,imh,iad,jad,krd    ! rajout de déclarations
72  c      integer iu,iiu,j2,jmr,js0,jt    ! integer Jmax,kmax,ndt0,nstep,k,j,i,ic,l,js,jn,imh,iad,jad,krd
73  c      real dtdy,dtdy5,rcap,iml,jn0,imjm,pi,dl,dp    ! integer iu,iiu,j2,jmr,js0,jt
74  c      real dt,cr1,maxdt,ztc,d5,sum1,sum2,ru    ! real dtdy,dtdy5,rcap,iml,jn0,imjm,pi,dl,dp
75  C    ! real dt,cr1,maxdt,ztc,d5,sum1,sum2,ru
76  C ********************************************************************  
77  C    ! ********************************************************************
78  C =============  
79  C INPUT:    ! =============
80  C =============    ! INPUT:
81  C    ! =============
82  C Q(IMR,JNP,NLAY,NC): mixing ratios at current time (t)  
83  C NC: total # of constituents    ! Q(IMR,JNP,NLAY,NC): mixing ratios at current time (t)
84  C IMR: first dimension (E-W); # of Grid intervals in E-W is IMR    ! NC: total number of constituents
85  C JNP: 2nd dimension (N-S); # of Grid intervals in N-S is JNP-1    ! IMR: first dimension (E-W); number of Grid intervals in E-W is IMR
86  C NLAY: 3rd dimension (# of layers); vertical index increases from 1 at    ! JNP: 2nd dimension (N-S); number of Grid intervals in N-S is JNP-1
87  C       the model top to NLAY near the surface (see fig. below).    ! NLAY: 3rd dimension (number of layers); vertical index increases from 1
88  C       It is assumed that 6 <= NLAY <= JNP (for dynamic memory allocation)    ! at
89  C    ! the model top to NLAY near the surface (see fig. below).
90  C PS1(IMR,JNP): surface pressure at current time (t)    ! It is assumed that 6 <= NLAY <= JNP (for dynamic memory allocation)
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.    ! PS1(IMR,JNP): surface pressure at current time (t)
93  C Note: surface pressure can have any unit or can be multiplied by any    ! PS2(IMR,JNP): surface pressure at mid-time-level (t+NDT/2)
94  C       const.    ! PS2 is replaced by the predicted PS (at t+NDT) on output.
95  C    ! Note: surface pressure can have any unit or can be multiplied by any
96  C The pressure at layer edges are defined as follows:    ! const.
97  C  
98  C        p(i,j,k) = AP(k)*PT  +  BP(k)*PS(i,j)          (1)    ! The pressure at layer edges are defined as follows:
99  C  
100  C Where PT is a constant having the same unit as PS.    ! p(i,j,k) = AP(k)*PT  +  BP(k)*PS(i,j)          (1)
101  C AP and BP are unitless constants given at layer edges  
102  C defining the vertical coordinate.    ! Where PT is a constant having the same unit as PS.
103  C BP(1) = 0., BP(NLAY+1) = 1.    ! AP and BP are unitless constants given at layer edges
104  C The pressure at the model top is PTOP = AP(1)*PT    ! defining the vertical coordinate.
105  C    ! BP(1) = 0., BP(NLAY+1) = 1.
106  C For pure sigma system set AP(k) = 1 for all k, PT = PTOP,    ! The pressure at the model top is PTOP = AP(1)*PT
107  C BP(k) = sige(k) (sigma at edges), PS = Psfc - PTOP.  
108  C    ! For pure sigma system set AP(k) = 1 for all k, PT = PTOP,
109  C Note: the sigma-P coordinate is a subset of Eq. 1, which in turn    ! BP(k) = sige(k) (sigma at edges), PS = Psfc - PTOP.
110  C is a subset of the following even more general sigma-P-thelta coord.  
111  C currently under development.    ! Note: the sigma-P coordinate is a subset of Eq. 1, which in turn
112  C  p(i,j,k) = (AP(k)*PT + BP(k)*PS(i,j))/(D(k)-C(k)*TE**(-1/kapa))    ! is a subset of the following even more general sigma-P-thelta coord.
113  C    ! currently under development.
114  C                  /////////////////////////////////    ! p(i,j,k) = (AP(k)*PT + BP(k)*PS(i,j))/(D(k)-C(k)*TE**(-1/kapa))
115  C              / \ ------------- PTOP --------------  AP(1), BP(1)  
116  C               |    ! /////////////////////////////////
117  C    delp(1)    |  ........... Q(i,j,1) ............      ! / \ ------------- PTOP --------------  AP(1), BP(1)
118  C               |    ! |
119  C      W(1)    \ / ---------------------------------  AP(2), BP(2)    ! delp(1)    |  ........... Q(i,j,1) ............
120  C    ! |
121  C    ! W(1)    \ / ---------------------------------  AP(2), BP(2)
122  C  
123  C     W(k-1)   / \ ---------------------------------  AP(k), BP(k)  
124  C               |  
125  C    delp(K)    |  ........... Q(i,j,k) ............    ! W(k-1)   / \ ---------------------------------  AP(k), BP(k)
126  C               |    ! |
127  C      W(k)    \ / ---------------------------------  AP(k+1), BP(k+1)    ! delp(K)    |  ........... Q(i,j,k) ............
128  C    ! |
129  C    ! W(k)    \ / ---------------------------------  AP(k+1), BP(k+1)
130  C  
131  C              / \ ---------------------------------  AP(NLAY), BP(NLAY)  
132  C               |  
133  C  delp(NLAY)   |  ........... Q(i,j,NLAY) .........      ! / \ ---------------------------------  AP(NLAY), BP(NLAY)
134  C               |    ! |
135  C   W(NLAY)=0  \ / ------------- surface ----------- AP(NLAY+1), BP(NLAY+1)    ! delp(NLAY)   |  ........... Q(i,j,NLAY) .........
136  C                 //////////////////////////////////    ! |
137  C    ! W(NLAY)=0  \ / ------------- surface ----------- AP(NLAY+1), BP(NLAY+1)
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.  
140  C    ! U(IMR,JNP,NLAY) & V(IMR,JNP,NLAY):winds (m/s) at mid-time-level (t+NDT/2)
141  C IGD:      grid type on which winds are defined.    ! U and V may need to be polar filtered in advance in some cases.
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) ]    ! IGD:      grid type on which winds are defined.
144  C    ! IGD = 0:  A-Grid  [all variables defined at the same point from south
145  C IGD = 1  GEOS-GCM C-Grid    ! pole (j=1) to north pole (j=JNP) ]
146  C                                      [North]  
147  C    ! IGD = 1  GEOS-GCM C-Grid
148  C                                       V(i,j)    ! [North]
149  C                                          |  
150  C                                          |    ! V(i,j)
151  C                                          |    ! |
152  C                             U(i-1,j)---Q(i,j)---U(i,j) [EAST]    ! |
153  C                                          |    ! |
154  C                                          |    ! U(i-1,j)---Q(i,j)---U(i,j) [EAST]
155  C                                          |    ! |
156  C                                       V(i,j-1)    ! |
157  C    ! |
158  C         U(i,  1) is defined at South Pole.    ! V(i,j-1)
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.    ! U(i,  1) is defined at South Pole.
161  C    ! V(i,  1) is half grid north of the South Pole.
162  C         V must be defined at j=1 and j=JMR if IGD=1    ! V(i,JMR) is half grid south of the North Pole.
163  C         V at JNP need not be given.  
164  C    ! V must be defined at j=1 and j=JMR if IGD=1
165  C NDT: time step in seconds (need not be constant during the course of    ! V at JNP need not be given.
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    ! NDT: time step in seconds (need not be constant during the course of
168  C      has a well-resolved stratosphere.    ! the integration). Suggested value: 30 min. for 4x5, 15 min. for 2x2.5
169  C    ! (Lat-Lon) resolution. Smaller values are recommanded if the model
170  C J1 defines the size of the polar cap:    ! has a well-resolved stratosphere.
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.    ! J1 defines the size of the polar cap:
173  C There are currently only two choices (j1=2 or 3).    ! South polar cap edge is located at -90 + (j1-1.5)*180/(JNP-1) deg.
174  C IMR must be an even integer if j1 = 2. Recommended value: J1=3.    ! North polar cap edge is located at  90 - (j1-1.5)*180/(JNP-1) deg.
175  C    ! There are currently only two choices (j1=2 or 3).
176  C IORD, JORD, and KORD are integers controlling various options in E-W, N-S,    ! IMR must be an even integer if j1 = 2. Recommended value: J1=3.
177  C and vertical transport, respectively. Recommended values for positive  
178  C definite scalars: IORD=JORD=3, KORD=5. Use KORD=3 for non-    ! IORD, JORD, and KORD are integers controlling various options in E-W,
179  C positive definite scalars or when linear correlation between constituents    ! N-S,
180  C is to be maintained.    ! and vertical transport, respectively. Recommended values for positive
181  C    ! definite scalars: IORD=JORD=3, KORD=5. Use KORD=3 for non-
182  C  _ORD=    ! positive definite scalars or when linear correlation between constituents
183  C        1: 1st order upstream scheme (too diffusive, not a useful option; it    ! is to be maintained.
184  C           can be used for debugging purposes; this is THE only known "linear"  
185  C           monotonic advection scheme.).    ! _ORD=
186  C        2: 2nd order van Leer (full monotonicity constraint;    ! 1: 1st order upstream scheme (too diffusive, not a useful option; it
187  C           see Lin et al 1994, MWR)    ! can be used for debugging purposes; this is THE only known "linear"
188  C        3: monotonic PPM* (slightly improved PPM of Collela & Woodward 1984)    ! monotonic advection scheme.).
189  C        4: semi-monotonic PPM (same as 3, but overshoots are allowed)    ! 2: 2nd order van Leer (full monotonicity constraint;
190  C        5: positive-definite PPM (constraint on the subgrid distribution is    ! see Lin et al 1994, MWR)
191  C           only strong enough to prevent generation of negative values;    ! 3: monotonic PPM* (slightly improved PPM of Collela & Woodward 1984)
192  C           both overshoots & undershoots are possible).    ! 4: semi-monotonic PPM (same as 3, but overshoots are allowed)
193  C        6: un-constrained PPM (nearly diffusion free; slightly faster but    ! 5: positive-definite PPM (constraint on the subgrid distribution is
194  C           positivity not quaranteed. Use this option only when the fields    ! only strong enough to prevent generation of negative values;
195  C           and winds are very smooth).    ! both overshoots & undershoots are possible).
196  C    ! 6: un-constrained PPM (nearly diffusion free; slightly faster but
197  C *PPM: Piece-wise Parabolic Method    ! positivity not quaranteed. Use this option only when the fields
198  C    ! and winds are very smooth).
199  C Note that KORD <=2 options are no longer supported. DO not use option 4 or 5.  
200  C for non-positive definite scalars (such as Ertel Potential Vorticity).    ! *PPM: Piece-wise Parabolic Method
201  C  
202  C The implicit numerical diffusion decreases as _ORD increases.    ! Note that KORD <=2 options are no longer supported. DO not use option 4
203  C The last two options (ORDER=5, 6) should only be used when there is    ! or 5.
204  C significant explicit diffusion (such as a turbulence parameterization). You    ! for non-positive definite scalars (such as Ertel Potential Vorticity).
205  C might get dispersive results otherwise.  
206  C No filter of any kind is applied to the constituent fields here.    ! The implicit numerical diffusion decreases as _ORD increases.
207  C    ! The last two options (ORDER=5, 6) should only be used when there is
208  C AE: Radius of the sphere (meters).    ! significant explicit diffusion (such as a turbulence parameterization).
209  C     Recommended value for the planet earth: 6.371E6    ! You
210  C    ! might get dispersive results otherwise.
211  C fill(logical):   flag to do filling for negatives (see note below).    ! No filter of any kind is applied to the constituent fields here.
212  C  
213  C Umax: Estimate (upper limit) of the maximum U-wind speed (m/s).    ! AE: Radius of the sphere (meters).
214  C (220 m/s is a good value for troposphere model; 280 m/s otherwise)    ! Recommended value for the planet earth: 6.371E6
215  C  
216  C =============    ! fill(logical):   flag to do filling for negatives (see note below).
217  C Output  
218  C =============    ! Umax: Estimate (upper limit) of the maximum U-wind speed (m/s).
219  C    ! (220 m/s is a good value for troposphere model; 280 m/s otherwise)
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    ! =============
222  C          relationship. W will have the same unit as PS1 and PS2 (eg, mb).    ! Output
223  C          W must be divided by NDT to get the correct mass-flux unit.    ! =============
224  C          The vertical Courant number C = W/delp_UPWIND, where delp_UPWIND  
225  C          is the pressure thickness in the "upwind" direction. For example,    ! Q: mixing ratios at future time (t+NDT) (original values are
226  C          C(k) = W(k)/delp(k)   if W(k) > 0;    ! over-written)
227  C          C(k) = W(k)/delp(k+1) if W(k) < 0.    ! W(NLAY): large-scale vertical mass flux as diagnosed from the hydrostatic
228  C              ( W > 0 is downward, ie, toward surface)    ! relationship. W will have the same unit as PS1 and PS2 (eg, mb).
229  C PS2: predicted PS at t+NDT (original values are over-written)    ! W must be divided by NDT to get the correct mass-flux unit.
230  C    ! The vertical Courant number C = W/delp_UPWIND, where delp_UPWIND
231  C ********************************************************************    ! is the pressure thickness in the "upwind" direction. For example,
232  C NOTES:    ! C(k) = W(k)/delp(k)   if W(k) > 0;
233  C This forward-in-time upstream-biased transport scheme reduces to    ! C(k) = W(k)/delp(k+1) if W(k) < 0.
234  C the 2nd order center-in-time center-in-space mass continuity eqn.    ! ( W > 0 is downward, ie, toward surface)
235  C if Q = 1 (constant fields will remain constant). This also ensures    ! PS2: predicted PS at t+NDT (original values are over-written)
236  C that the computed vertical velocity to be identical to GEOS-1 GCM  
237  C for on-line transport.    ! ********************************************************************
238  C    ! NOTES:
239  C A larger polar cap is used if j1=3 (recommended for C-Grid winds or when    ! This forward-in-time upstream-biased transport scheme reduces to
240  C winds are noisy near poles).    ! the 2nd order center-in-time center-in-space mass continuity eqn.
241  C    ! if Q = 1 (constant fields will remain constant). This also ensures
242  C Flux-Form Semi-Lagrangian transport in the East-West direction is used    ! that the computed vertical velocity to be identical to GEOS-1 GCM
243  C when and where Courant # is greater than one.    ! for on-line transport.
244  C  
245  C The user needs to change the parameter Jmax or Kmax if the resolution    ! A larger polar cap is used if j1=3 (recommended for C-Grid winds or when
246  C is greater than 0.5 deg in N-S or 150 layers in the vertical direction.    ! winds are noisy near poles).
247  C (this TransPort Core is otherwise resolution independent and can be used  
248  C as a library routine).    ! Flux-Form Semi-Lagrangian transport in the East-West direction is used
249  C    ! when and where Courant number is greater than one.
250  C PPM is 4th order accurate when grid spacing is uniform (x & y); 3rd  
251  C order accurate for non-uniform grid (vertical sigma coord.).    ! The user needs to change the parameter Jmax or Kmax if the resolution
252  C    ! is greater than 0.5 deg in N-S or 150 layers in the vertical direction.
253  C Time step is limitted only by transport in the meridional direction.    ! (this TransPort Core is otherwise resolution independent and can be used
254  C (the FFSL scheme is not implemented in the meridional direction).    ! as a library routine).
255  C  
256  C Since only 1-D limiters are applied, negative values could    ! PPM is 4th order accurate when grid spacing is uniform (x & y); 3rd
257  C potentially be generated when large time step is used and when the    ! order accurate for non-uniform grid (vertical sigma coord.).
258  C initial fields contain discontinuities.  
259  C This does not necessarily imply the integration is unstable.    ! Time step is limitted only by transport in the meridional direction.
260  C These negatives are typically very small. A filling algorithm is    ! (the FFSL scheme is not implemented in the meridional direction).
261  C activated if the user set "fill" to be true.  
262  C    ! Since only 1-D limiters are applied, negative values could
263  C The van Leer scheme used here is nearly as accurate as the original PPM    ! potentially be generated when large time step is used and when the
264  C due to the use of a 4th order accurate reference slope. The PPM imple-    ! initial fields contain discontinuities.
265  C mented here is an improvement over the original and is also based on    ! This does not necessarily imply the integration is unstable.
266  C the 4th order reference slope.    ! These negatives are typically very small. A filling algorithm is
267  C    ! activated if the user set "fill" to be true.
268  C ****6***0*********0*********0*********0*********0*********0**********72  
269  C    ! The van Leer scheme used here is nearly as accurate as the original PPM
270  C     User modifiable parameters    ! due to the use of a 4th order accurate reference slope. The PPM imple-
271  C    ! mented here is an improvement over the original and is also based on
272        parameter (Jmax = 361, kmax = 150)    ! the 4th order reference slope.
273  C  
274  C ****6***0*********0*********0*********0*********0*********0**********72    ! ****6***0*********0*********0*********0*********0*********0**********72
275  C  
276  C Input-Output arrays    ! User modifiable parameters
277  C  
278            PARAMETER (jmax=361, kmax=150)
279        real Q(IMR,JNP,NLAY,NC),PS1(IMR,JNP),PS2(IMR,JNP),  
280       &     U(IMR,JNP,NLAY),V(IMR,JNP,NLAY),AP(NLAY+1),    ! ****6***0*********0*********0*********0*********0*********0**********72
281       &     BP(NLAY+1),W(IMR,JNP,NLAY),NDT,val(NLAY),Umax  
282        integer IGD,IORD,JORD,KORD,NC,IMR,JNP,j1,NLAY,AE    ! Input-Output arrays
283        integer IMRD2  
284        real    PT        
285        logical  cross, fill, dum    REAL q(imr, jnp, nlay, nc), ps1(imr, jnp), ps2(imr, jnp), &
286  C      u(imr, jnp, nlay), v(imr, jnp, nlay), ap(nlay+1), bp(nlay+1), &
287  C Local dynamic arrays      w(imr, jnp, nlay), ndt, val(nlay), umax
288  C    INTEGER igd, iord, jord, kord, nc, imr, jnp, j1, nlay, ae
289        real CRX(IMR,JNP),CRY(IMR,JNP),xmass(IMR,JNP),ymass(IMR,JNP),    INTEGER imrd2
290       &     fx1(IMR+1),DPI(IMR,JNP,NLAY),delp1(IMR,JNP,NLAY),    REAL pt
291       &     WK1(IMR,JNP,NLAY),PU(IMR,JNP),PV(IMR,JNP),DC2(IMR,JNP),    LOGICAL cross, fill, dum
292       &     delp2(IMR,JNP,NLAY),DQ(IMR,JNP,NLAY,NC),VA(IMR,JNP),  
293       &     UA(IMR,JNP),qtmp(-IMR:2*IMR)    ! Local dynamic arrays
294  C  
295  C Local static  arrays    REAL crx(imr, jnp), cry(imr, jnp), xmass(imr, jnp), ymass(imr, jnp), &
296  C      fx1(imr+1), dpi(imr, jnp, nlay), delp1(imr, jnp, nlay), &
297        real DTDX(Jmax), DTDX5(Jmax), acosp(Jmax),      wk1(imr, jnp, nlay), pu(imr, jnp), pv(imr, jnp), dc2(imr, jnp), &
298       &     cosp(Jmax), cose(Jmax), DAP(kmax),DBK(Kmax)      delp2(imr, jnp, nlay), dq(imr, jnp, nlay, nc), va(imr, jnp), &
299        data NDT0, NSTEP /0, 0/      ua(imr, jnp), qtmp(-imr:2*imr)
300        data cross /.true./  
301        SAVE DTDY, DTDY5, RCAP, JS0, JN0, IML,    ! Local static  arrays
302       &     DTDX, DTDX5, ACOSP, COSP, COSE, DAP,DBK  
303  C    REAL dtdx(jmax), dtdx5(jmax), acosp(jmax), cosp(jmax), cose(jmax), &
304                    dap(kmax), dbk(kmax)
305        JMR = JNP -1    DATA ndt0, nstep/0, 0/
306        IMJM  = IMR*JNP    DATA cross/.TRUE./
307        j2 = JNP - j1 + 1    SAVE dtdy, dtdy5, rcap, js0, jn0, iml, dtdx, dtdx5, acosp, cosp, cose, dap, &
308        NSTEP = NSTEP + 1      dbk
309  C  
310  C *********** Initialization **********************  
311        if(NSTEP.eq.1) then    jmr = jnp - 1
312  c    imjm = imr*jnp
313        write(6,*) '------------------------------------ '    j2 = jnp - j1 + 1
314        write(6,*) 'NASA/GSFC Transport Core Version 4.5'    nstep = nstep + 1
315        write(6,*) '------------------------------------ '  
316  c    ! *********** Initialization **********************
317        WRITE(6,*) 'IMR=',IMR,' JNP=',JNP,' NLAY=',NLAY,' j1=',j1    IF (nstep==1) THEN
318        WRITE(6,*) 'NC=',NC,IORD,JORD,KORD,NDT  
319  C      WRITE (6, *) '------------------------------------ '
320  C controles sur les parametres      WRITE (6, *) 'NASA/GSFC Transport Core Version 4.5'
321        if(NLAY.LT.6) then      WRITE (6, *) '------------------------------------ '
322          write(6,*) 'NLAY must be >= 6'  
323          stop      WRITE (6, *) 'IMR=', imr, ' JNP=', jnp, ' NLAY=', nlay, ' j1=', j1
324        endif      WRITE (6, *) 'NC=', nc, iord, jord, kord, ndt
325        if (JNP.LT.NLAY) then  
326           write(6,*) 'JNP must be >= NLAY'      ! controles sur les parametres
327          stop      IF (nlay<6) THEN
328        endif        WRITE (6, *) 'NLAY must be >= 6'
329        IMRD2=mod(IMR,2)        STOP
330        if (j1.eq.2.and.IMRD2.NE.0) then      END IF
331           write(6,*) 'if j1=2 IMR must be an even integer'      IF (jnp<nlay) THEN
332          stop        WRITE (6, *) 'JNP must be >= NLAY'
333        endif        STOP
334        END IF
335  C      imrd2 = mod(imr, 2)
336        if(Jmax.lt.JNP .or. Kmax.lt.NLAY) then      IF (j1==2 .AND. imrd2/=0) THEN
337          write(6,*) 'Jmax or Kmax is too small'        WRITE (6, *) 'if j1=2 IMR must be an even integer'
338          stop        STOP
339        endif      END IF
340  C  
341        DO k=1,NLAY  
342        DAP(k) = (AP(k+1) - AP(k))*PT      IF (jmax<jnp .OR. kmax<nlay) THEN
343        DBK(k) =  BP(k+1) - BP(k)        WRITE (6, *) 'Jmax or Kmax is too small'
344        ENDDO            STOP
345  C      END IF
346        PI = 4. * ATAN(1.)  
347        DL = 2.*PI / float(IMR)      DO k = 1, nlay
348        DP =    PI / float(JMR)        dap(k) = (ap(k+1)-ap(k))*pt
349  C        dbk(k) = bp(k+1) - bp(k)
350        if(IGD.eq.0) then      END DO
351  C Compute analytic cosine at cell edges  
352              call cosa(cosp,cose,JNP,PI,DP)      pi = 4.*atan(1.)
353        else      dl = 2.*pi/float(imr)
354  C Define cosine consistent with GEOS-GCM (using dycore2.0 or later)      dp = pi/float(jmr)
355              call cosc(cosp,cose,JNP,PI,DP)  
356        endif      IF (igd==0) THEN
357  C        ! Compute analytic cosine at cell edges
358        do 15 J=2,JMR        CALL cosa(cosp, cose, jnp, pi, dp)
359  15    acosp(j) = 1. / cosp(j)      ELSE
360  C        ! Define cosine consistent with GEOS-GCM (using dycore2.0 or later)
361  C Inverse of the Scaled polar cap area.        CALL cosc(cosp, cose, jnp, pi, dp)
362  C      END IF
363        RCAP  = DP / (IMR*(1.-COS((j1-1.5)*DP)))  
364        acosp(1)   = RCAP      DO j = 2, jmr
365        acosp(JNP) = RCAP        acosp(j) = 1./cosp(j)
366        endif      END DO
367  C  
368        if(NDT0 .ne. NDT) then      ! Inverse of the Scaled polar cap area.
369        DT   = NDT  
370        NDT0 = NDT      rcap = dp/(imr*(1.-cos((j1-1.5)*dp)))
371        acosp(1) = rcap
372          if(Umax .lt. 180.) then      acosp(jnp) = rcap
373           write(6,*) 'Umax may be too small!'    END IF
374          endif  
375        CR1  = abs(Umax*DT)/(DL*AE)    IF (ndt0/=ndt) THEN
376        MaxDT = DP*AE / abs(Umax) + 0.5      dt = ndt
377        write(6,*)'Largest time step for max(V)=',Umax,' is ',MaxDT      ndt0 = ndt
378        if(MaxDT .lt. abs(NDT)) then  
379              write(6,*) 'Warning!!! NDT maybe too large!'      IF (umax<180.) THEN
380        endif        WRITE (6, *) 'Umax may be too small!'
381  C      END IF
382        if(CR1.ge.0.95) then      cr1 = abs(umax*dt)/(dl*ae)
383        JS0 = 0      maxdt = dp*ae/abs(umax) + 0.5
384        JN0 = 0      WRITE (6, *) 'Largest time step for max(V)=', umax, ' is ', maxdt
385        IML = IMR-2      IF (maxdt<abs(ndt)) THEN
386        ZTC = 0.        WRITE (6, *) 'Warning!!! NDT maybe too large!'
387        else      END IF
388        ZTC  = acos(CR1) * (180./PI)  
389  C      IF (cr1>=0.95) THEN
390        JS0 = float(JMR)*(90.-ZTC)/180. + 2        js0 = 0
391        JS0 = max(JS0, J1+1)        jn0 = 0
392        IML = min(6*JS0/(J1-1)+2, 4*IMR/5)        iml = imr - 2
393        JN0 = JNP-JS0+1        ztc = 0.
394        endif      ELSE
395  C            ztc = acos(cr1)*(180./pi)
396  C  
397        do J=2,JMR        js0 = float(jmr)*(90.-ztc)/180. + 2
398        DTDX(j)  = DT / ( DL*AE*COSP(J) )        js0 = max(js0, j1+1)
399          iml = min(6*js0/(j1-1)+2, 4*imr/5)
400  c     print*,'dtdx=',dtdx(j)        jn0 = jnp - js0 + 1
401        DTDX5(j) = 0.5*DTDX(j)      END IF
402        enddo  
403  C  
404              DO j = 2, jmr
405        DTDY  = DT /(AE*DP)        dtdx(j) = dt/(dl*ae*cosp(j))
406  c      print*,'dtdy=',dtdy  
407        DTDY5 = 0.5*DTDY        dtdx5(j) = 0.5*dtdx(j)
408  C      END DO
409  c      write(6,*) 'J1=',J1,' J2=', J2  
410        endif  
411  C      dtdy = dt/(ae*dp)
412  C *********** End Initialization **********************      dtdy5 = 0.5*dtdy
413  C  
414  C delp = pressure thickness: the psudo-density in a hydrostatic system.    END IF
415        do  k=1,NLAY  
416           do  j=1,JNP    ! *********** End Initialization **********************
417              do  i=1,IMR  
418                 delp1(i,j,k)=DAP(k)+DBK(k)*PS1(i,j)    ! delp = pressure thickness: the psudo-density in a hydrostatic system.
419                 delp2(i,j,k)=DAP(k)+DBK(k)*PS2(i,j)          DO k = 1, nlay
420              enddo      DO j = 1, jnp
421           enddo        DO i = 1, imr
422        enddo          delp1(i, j, k) = dap(k) + dbk(k)*ps1(i, j)
423                      delp2(i, j, k) = dap(k) + dbk(k)*ps2(i, j)
424  C        END DO
425        if(j1.ne.2) then      END DO
426        DO 40 IC=1,NC    END DO
427        DO 40 L=1,NLAY  
428        DO 40 I=1,IMR  
429        Q(I,  2,L,IC) = Q(I,  1,L,IC)    IF (j1/=2) THEN
430  40    Q(I,JMR,L,IC) = Q(I,JNP,L,IC)      DO ic = 1, nc
431        endif        DO l = 1, nlay
432  C          DO i = 1, imr
433  C Compute "tracer density"            q(i, 2, l, ic) = q(i, 1, l, ic)
434        DO 550 IC=1,NC            q(i, jmr, l, ic) = q(i, jnp, l, ic)
435        DO 44 k=1,NLAY          END DO
436        DO 44 j=1,JNP        END DO
437        DO 44 i=1,IMR      END DO
438  44    DQ(i,j,k,IC) = Q(i,j,k,IC)*delp1(i,j,k)    END IF
439  550     continue  
440  C    ! Compute "tracer density"
441        do 1500 k=1,NLAY    DO ic = 1, nc
442  C      DO k = 1, nlay
443        if(IGD.eq.0) then        DO j = 1, jnp
444  C Convert winds on A-Grid to Courant # on C-Grid.          DO i = 1, imr
445        call A2C(U(1,1,k),V(1,1,k),IMR,JMR,j1,j2,CRX,CRY,dtdx5,DTDY5)            dq(i, j, k, ic) = q(i, j, k, ic)*delp1(i, j, k)
446        else          END DO
447  C Convert winds on C-grid to Courant #        END DO
448        do 45 j=j1,j2      END DO
449        do 45 i=2,IMR    END DO
450  45    CRX(i,J) = dtdx(j)*U(i-1,j,k)  
451        DO k = 1, nlay
452  C  
453        do 50 j=j1,j2      IF (igd==0) THEN
454  50    CRX(1,J) = dtdx(j)*U(IMR,j,k)        ! Convert winds on A-Grid to Courant number on C-Grid.
455  C        CALL a2c(u(1,1,k), v(1,1,k), imr, jmr, j1, j2, crx, cry, dtdx5, dtdy5)
456        do 55 i=1,IMR*JMR      ELSE
457  55    CRY(i,2) = DTDY*V(i,1,k)        ! Convert winds on C-grid to Courant number
458        endif        DO j = j1, j2
459  C              DO i = 2, imr
460  C Determine JS and JN            crx(i, j) = dtdx(j)*u(i-1, j, k)
461        JS = j1          END DO
462        JN = j2        END DO
463  C  
464        do j=JS0,j1+1,-1  
465        do i=1,IMR        DO j = j1, j2
466        if(abs(CRX(i,j)).GT.1.) then          crx(1, j) = dtdx(j)*u(imr, j, k)
467              JS = j        END DO
468              go to 2222  
469        endif        DO i = 1, imr*jmr
470        enddo          cry(i, 2) = dtdy*v(i, 1, k)
471        enddo        END DO
472  C      END IF
473  2222  continue  
474        do j=JN0,j2-1      ! Determine JS and JN
475        do i=1,IMR      js = j1
476        if(abs(CRX(i,j)).GT.1.) then      jn = j2
477              JN = j  
478              go to 2233      DO j = js0, j1 + 1, -1
479        endif        DO i = 1, imr
480        enddo          IF (abs(crx(i,j))>1.) THEN
481        enddo            js = j
482  2233  continue            GO TO 2222
483  C          END IF
484        if(j1.ne.2) then           ! Enlarged polar cap.        END DO
485        do i=1,IMR      END DO
486        DPI(i,  2,k) = 0.  
487        DPI(i,JMR,k) = 0.  2222 CONTINUE
488        enddo      DO j = jn0, j2 - 1
489        endif        DO i = 1, imr
490  C          IF (abs(crx(i,j))>1.) THEN
491  C ******* Compute horizontal mass fluxes ************            jn = j
492  C            GO TO 2233
493  C N-S component          END IF
494        do j=j1,j2+1        END DO
495        D5 = 0.5 * COSE(j)      END DO
496        do i=1,IMR  2233 CONTINUE
497        ymass(i,j) = CRY(i,j)*D5*(delp2(i,j,k) + delp2(i,j-1,k))  
498        enddo      IF (j1/=2) THEN ! Enlarged polar cap.
499        enddo        DO i = 1, imr
500  C          dpi(i, 2, k) = 0.
501        do 95 j=j1,j2          dpi(i, jmr, k) = 0.
502        DO 95 i=1,IMR        END DO
503  95    DPI(i,j,k) = (ymass(i,j) - ymass(i,j+1)) * acosp(j)      END IF
504  C  
505  C Poles      ! ******* Compute horizontal mass fluxes ************
506        sum1 = ymass(IMR,j1  )  
507        sum2 = ymass(IMR,J2+1)      ! N-S component
508        do i=1,IMR-1      DO j = j1, j2 + 1
509        sum1 = sum1 + ymass(i,j1  )        d5 = 0.5*cose(j)
510        sum2 = sum2 + ymass(i,J2+1)        DO i = 1, imr
511        enddo          ymass(i, j) = cry(i, j)*d5*(delp2(i,j,k)+delp2(i,j-1,k))
512  C        END DO
513        sum1 = - sum1 * RCAP      END DO
514        sum2 =   sum2 * RCAP  
515        do i=1,IMR      DO j = j1, j2
516        DPI(i,  1,k) = sum1        DO i = 1, imr
517        DPI(i,JNP,k) = sum2          dpi(i, j, k) = (ymass(i,j)-ymass(i,j+1))*acosp(j)
518        enddo        END DO
519  C      END DO
520  C E-W component  
521  C      ! Poles
522        do j=j1,j2      sum1 = ymass(imr, j1)
523        do i=2,IMR      sum2 = ymass(imr, j2+1)
524        PU(i,j) = 0.5 * (delp2(i,j,k) + delp2(i-1,j,k))      DO i = 1, imr - 1
525        enddo        sum1 = sum1 + ymass(i, j1)
526        enddo        sum2 = sum2 + ymass(i, j2+1)
527  C      END DO
528        do j=j1,j2  
529        PU(1,j) = 0.5 * (delp2(1,j,k) + delp2(IMR,j,k))      sum1 = -sum1*rcap
530        enddo      sum2 = sum2*rcap
531  C      DO i = 1, imr
532        do 110 j=j1,j2        dpi(i, 1, k) = sum1
533        DO 110 i=1,IMR        dpi(i, jnp, k) = sum2
534  110   xmass(i,j) = PU(i,j)*CRX(i,j)      END DO
535  C  
536        DO 120 j=j1,j2      ! E-W component
537        DO 120 i=1,IMR-1  
538  120   DPI(i,j,k) = DPI(i,j,k) + xmass(i,j) - xmass(i+1,j)      DO j = j1, j2
539  C        DO i = 2, imr
540        DO 130 j=j1,j2          pu(i, j) = 0.5*(delp2(i,j,k)+delp2(i-1,j,k))
541  130   DPI(IMR,j,k) = DPI(IMR,j,k) + xmass(IMR,j) - xmass(1,j)        END DO
542  C      END DO
543        DO j=j1,j2  
544        do i=1,IMR-1      DO j = j1, j2
545        UA(i,j) = 0.5 * (CRX(i,j)+CRX(i+1,j))        pu(1, j) = 0.5*(delp2(1,j,k)+delp2(imr,j,k))
546        enddo      END DO
547        enddo  
548  C      DO j = j1, j2
549        DO j=j1,j2        DO i = 1, imr
550        UA(imr,j) = 0.5 * (CRX(imr,j)+CRX(1,j))          xmass(i, j) = pu(i, j)*crx(i, j)
551        enddo        END DO
552  ccccccccccccccccccccccccccccccccccccccccccccccccccccccc      END DO
553  c Rajouts pour LMDZ.3.3  
554  ccccccccccccccccccccccccccccccccccccccccccccccccccccccc      DO j = j1, j2
555        do i=1,IMR        DO i = 1, imr - 1
556           do j=1,JNP          dpi(i, j, k) = dpi(i, j, k) + xmass(i, j) - xmass(i+1, j)
557               VA(i,j)=0.        END DO
558           enddo      END DO
559        enddo  
560        DO j = j1, j2
561        do i=1,imr*(JMR-1)        dpi(imr, j, k) = dpi(imr, j, k) + xmass(imr, j) - xmass(1, j)
562        VA(i,2) = 0.5*(CRY(i,2)+CRY(i,3))      END DO
563        enddo  
564  C      DO j = j1, j2
565        if(j1.eq.2) then        DO i = 1, imr - 1
566          IMH = IMR/2          ua(i, j) = 0.5*(crx(i,j)+crx(i+1,j))
567        do i=1,IMH        END DO
568        VA(i,      1) = 0.5*(CRY(i,2)-CRY(i+IMH,2))      END DO
569        VA(i+IMH,  1) = -VA(i,1)  
570        VA(i,    JNP) = 0.5*(CRY(i,JNP)-CRY(i+IMH,JMR))      DO j = j1, j2
571        VA(i+IMH,JNP) = -VA(i,JNP)        ua(imr, j) = 0.5*(crx(imr,j)+crx(1,j))
572        enddo      END DO
573        VA(IMR,1)=VA(1,1)      ! cccccccccccccccccccccccccccccccccccccccccccccccccccccc
574        VA(IMR,JNP)=VA(1,JNP)      ! Rajouts pour LMDZ.3.3
575        endif      ! cccccccccccccccccccccccccccccccccccccccccccccccccccccc
576  C      DO i = 1, imr
577  C ****6***0*********0*********0*********0*********0*********0**********72        DO j = 1, jnp
578        do 1000 IC=1,NC          va(i, j) = 0.
579  C        END DO
580        do i=1,IMJM      END DO
581        wk1(i,1,1) = 0.  
582        wk1(i,1,2) = 0.      DO i = 1, imr*(jmr-1)
583        enddo        va(i, 2) = 0.5*(cry(i,2)+cry(i,3))
584  C      END DO
585  C E-W advective cross term  
586        do 250 j=J1,J2      IF (j1==2) THEN
587        if(J.GT.JS  .and. J.LT.JN) GO TO 250        imh = imr/2
588  C        DO i = 1, imh
589        do i=1,IMR          va(i, 1) = 0.5*(cry(i,2)-cry(i+imh,2))
590        qtmp(i) = q(i,j,k,IC)          va(i+imh, 1) = -va(i, 1)
591        enddo          va(i, jnp) = 0.5*(cry(i,jnp)-cry(i+imh,jmr))
592  C          va(i+imh, jnp) = -va(i, jnp)
593        do i=-IML,0        END DO
594        qtmp(i)       = q(IMR+i,j,k,IC)        va(imr, 1) = va(1, 1)
595        qtmp(IMR+1-i) = q(1-i,j,k,IC)        va(imr, jnp) = va(1, jnp)
596        enddo      END IF
597  C  
598        DO 230 i=1,IMR      ! ****6***0*********0*********0*********0*********0*********0**********72
599        iu = UA(i,j)      DO ic = 1, nc
600        ru = UA(i,j) - iu  
601        iiu = i-iu        DO i = 1, imjm
602        if(UA(i,j).GE.0.) then          wk1(i, 1, 1) = 0.
603        wk1(i,j,1) = qtmp(iiu)+ru*(qtmp(iiu-1)-qtmp(iiu))          wk1(i, 1, 2) = 0.
604        else        END DO
605        wk1(i,j,1) = qtmp(iiu)+ru*(qtmp(iiu)-qtmp(iiu+1))  
606        endif        ! E-W advective cross term
607        wk1(i,j,1) = wk1(i,j,1) - qtmp(i)        DO j = j1, j2
608  230   continue          IF (j>js .AND. j<jn) GO TO 250
609  250   continue  
610  C          DO i = 1, imr
611        if(JN.ne.0) then            qtmp(i) = q(i, j, k, ic)
612        do j=JS+1,JN-1          END DO
613  C  
614        do i=1,IMR          DO i = -iml, 0
615        qtmp(i) = q(i,j,k,IC)            qtmp(i) = q(imr+i, j, k, ic)
616        enddo            qtmp(imr+1-i) = q(1-i, j, k, ic)
617  C          END DO
618        qtmp(0)     = q(IMR,J,k,IC)  
619        qtmp(IMR+1) = q(  1,J,k,IC)          DO i = 1, imr
620  C            iu = ua(i, j)
621        do i=1,imr            ru = ua(i, j) - iu
622        iu = i - UA(i,j)            iiu = i - iu
623        wk1(i,j,1) = UA(i,j)*(qtmp(iu) - qtmp(iu+1))            IF (ua(i,j)>=0.) THEN
624        enddo              wk1(i, j, 1) = qtmp(iiu) + ru*(qtmp(iiu-1)-qtmp(iiu))
625        enddo            ELSE
626        endif              wk1(i, j, 1) = qtmp(iiu) + ru*(qtmp(iiu)-qtmp(iiu+1))
627  C ****6***0*********0*********0*********0*********0*********0**********72            END IF
628  C Contribution from the N-S advection            wk1(i, j, 1) = wk1(i, j, 1) - qtmp(i)
629        do i=1,imr*(j2-j1+1)          END DO
630        JT = float(J1) - VA(i,j1)  250   END DO
631        wk1(i,j1,2) = VA(i,j1) * (q(i,jt,k,IC) - q(i,jt+1,k,IC))  
632        enddo        IF (jn/=0) THEN
633  C          DO j = js + 1, jn - 1
634        do i=1,IMJM  
635        wk1(i,1,1) = q(i,1,k,IC) + 0.5*wk1(i,1,1)            DO i = 1, imr
636        wk1(i,1,2) = q(i,1,k,IC) + 0.5*wk1(i,1,2)              qtmp(i) = q(i, j, k, ic)
637        enddo            END DO
638  C  
639          if(cross) then            qtmp(0) = q(imr, j, k, ic)
640  C Add cross terms in the vertical direction.            qtmp(imr+1) = q(1, j, k, ic)
641          if(IORD .GE. 2) then  
642                  iad = 2            DO i = 1, imr
643          else              iu = i - ua(i, j)
644                  iad = 1              wk1(i, j, 1) = ua(i, j)*(qtmp(iu)-qtmp(iu+1))
645          endif            END DO
646  C          END DO
647          if(JORD .GE. 2) then        END IF
648                  jad = 2        ! ****6***0*********0*********0*********0*********0*********0**********72
649          else        ! Contribution from the N-S advection
650                  jad = 1        DO i = 1, imr*(j2-j1+1)
651          endif          jt = float(j1) - va(i, j1)
652        call xadv(IMR,JNP,j1,j2,wk1(1,1,2),UA,JS,JN,IML,DC2,iad)          wk1(i, j1, 2) = va(i, j1)*(q(i,jt,k,ic)-q(i,jt+1,k,ic))
653        call yadv(IMR,JNP,j1,j2,wk1(1,1,1),VA,PV,W,jad)        END DO
654        do j=1,JNP  
655        do i=1,IMR        DO i = 1, imjm
656        q(i,j,k,IC) = q(i,j,k,IC) + DC2(i,j) + PV(i,j)          wk1(i, 1, 1) = q(i, 1, k, ic) + 0.5*wk1(i, 1, 1)
657        enddo          wk1(i, 1, 2) = q(i, 1, k, ic) + 0.5*wk1(i, 1, 2)
658        enddo        END DO
659        endif  
660  C        IF (cross) THEN
661        call xtp(IMR,JNP,IML,j1,j2,JN,JS,PU,DQ(1,1,k,IC),wk1(1,1,2)          ! Add cross terms in the vertical direction.
662       &        ,CRX,fx1,xmass,IORD)          IF (iord>=2) THEN
663              iad = 2
664        call ytp(IMR,JNP,j1,j2,acosp,RCAP,DQ(1,1,k,IC),wk1(1,1,1),CRY,          ELSE
665       &  DC2,ymass,WK1(1,1,3),wk1(1,1,4),WK1(1,1,5),WK1(1,1,6),JORD)            iad = 1
666  C          END IF
667  1000  continue  
668  1500  continue          IF (jord>=2) THEN
669  C            jad = 2
670  C ******* Compute vertical mass flux (same unit as PS) ***********          ELSE
671  C            jad = 1
672  C 1st step: compute total column mass CONVERGENCE.          END IF
673  C          CALL xadv(imr, jnp, j1, j2, wk1(1,1,2), ua, js, jn, iml, dc2, iad)
674        do 320 j=1,JNP          CALL yadv(imr, jnp, j1, j2, wk1(1,1,1), va, pv, w, jad)
675        do 320 i=1,IMR          DO j = 1, jnp
676  320   CRY(i,j) = DPI(i,j,1)            DO i = 1, imr
677  C              q(i, j, k, ic) = q(i, j, k, ic) + dc2(i, j) + pv(i, j)
678        do 330 k=2,NLAY            END DO
679        do 330 j=1,JNP          END DO
680        do 330 i=1,IMR        END IF
681        CRY(i,j)  = CRY(i,j) + DPI(i,j,k)  
682  330   continue        CALL xtp(imr, jnp, iml, j1, j2, jn, js, pu, dq(1,1,k,ic), wk1(1,1,2), &
683  C          crx, fx1, xmass, iord)
684        do 360 j=1,JNP  
685        do 360 i=1,IMR        CALL ytp(imr, jnp, j1, j2, acosp, rcap, dq(1,1,k,ic), wk1(1,1,1), cry, &
686  C          dc2, ymass, wk1(1,1,3), wk1(1,1,4), wk1(1,1,5), wk1(1,1,6), jord)
687  C 2nd step: compute PS2 (PS at n+1) using the hydrostatic assumption.  
688  C Changes (increases) to surface pressure = total column mass convergence      END DO
689  C    END DO
690        PS2(i,j)  = PS1(i,j) + CRY(i,j)  
691  C    ! ******* Compute vertical mass flux (same unit as PS) ***********
692  C 3rd step: compute vertical mass flux from mass conservation principle.  
693  C    ! 1st step: compute total column mass CONVERGENCE.
694        W(i,j,1) = DPI(i,j,1) - DBK(1)*CRY(i,j)  
695        W(i,j,NLAY) = 0.    DO j = 1, jnp
696  360   continue      DO i = 1, imr
697  C        cry(i, j) = dpi(i, j, 1)
698        do 370 k=2,NLAY-1      END DO
699        do 370 j=1,JNP    END DO
700        do 370 i=1,IMR  
701        W(i,j,k) = W(i,j,k-1) + DPI(i,j,k) - DBK(k)*CRY(i,j)    DO k = 2, nlay
702  370   continue      DO j = 1, jnp
703  C        DO i = 1, imr
704        DO 380 k=1,NLAY          cry(i, j) = cry(i, j) + dpi(i, j, k)
705        DO 380 j=1,JNP        END DO
706        DO 380 i=1,IMR      END DO
707        delp2(i,j,k) = DAP(k) + DBK(k)*PS2(i,j)    END DO
708  380   continue  
709  C    DO j = 1, jnp
710          KRD = max(3, KORD)      DO i = 1, imr
711        do 4000 IC=1,NC  
712  C        ! 2nd step: compute PS2 (PS at n+1) using the hydrostatic assumption.
713  C****6***0*********0*********0*********0*********0*********0**********72        ! Changes (increases) to surface pressure = total column mass
714            ! convergence
715        call FZPPM(IMR,JNP,NLAY,j1,DQ(1,1,1,IC),W,Q(1,1,1,IC),WK1,DPI,  
716       &           DC2,CRX,CRY,PU,PV,xmass,ymass,delp1,KRD)        ps2(i, j) = ps1(i, j) + cry(i, j)
717  C  
718              ! 3rd step: compute vertical mass flux from mass conservation
719        if(fill) call qckxyz(DQ(1,1,1,IC),DC2,IMR,JNP,NLAY,j1,j2,        ! principle.
720       &                     cosp,acosp,.false.,IC,NSTEP)  
721  C        w(i, j, 1) = dpi(i, j, 1) - dbk(1)*cry(i, j)
722  C Recover tracer mixing ratio from "density" using predicted        w(i, j, nlay) = 0.
723  C "air density" (pressure thickness) at time-level n+1      END DO
724  C    END DO
725        DO k=1,NLAY  
726        DO j=1,JNP    DO k = 2, nlay - 1
727        DO i=1,IMR      DO j = 1, jnp
728              Q(i,j,k,IC) = DQ(i,j,k,IC) / delp2(i,j,k)        DO i = 1, imr
729  c            print*,'i=',i,'j=',j,'k=',k,'Q(i,j,k,IC)=',Q(i,j,k,IC)          w(i, j, k) = w(i, j, k-1) + dpi(i, j, k) - dbk(k)*cry(i, j)
730        enddo        END DO
731        enddo      END DO
732        enddo    END DO
733  C      
734        if(j1.ne.2) then    DO k = 1, nlay
735        DO 400 k=1,NLAY      DO j = 1, jnp
736        DO 400 I=1,IMR        DO i = 1, imr
737  c     j=1 c'est le pôle Sud, j=JNP c'est le pôle Nord          delp2(i, j, k) = dap(k) + dbk(k)*ps2(i, j)
738        Q(I,  2,k,IC) = Q(I,  1,k,IC)        END DO
739        Q(I,JMR,k,IC) = Q(I,JMP,k,IC)      END DO
740  400   CONTINUE    END DO
741        endif  
742  4000  continue    krd = max(3, kord)
743  C    DO ic = 1, nc
744        if(j1.ne.2) then  
745        DO 5000 k=1,NLAY      ! ****6***0*********0*********0*********0*********0*********0**********72
746        DO 5000 i=1,IMR  
747        W(i,  2,k) = W(i,  1,k)      CALL fzppm(imr, jnp, nlay, j1, dq(1,1,1,ic), w, q(1,1,1,ic), wk1, dpi, &
748        W(i,JMR,k) = W(i,JNP,k)        dc2, crx, cry, pu, pv, xmass, ymass, delp1, krd)
749  5000  continue  
750        endif  
751  C      IF (fill) CALL qckxyz(dq(1,1,1,ic), dc2, imr, jnp, nlay, j1, j2, cosp, &
752        RETURN        acosp, .FALSE., ic, nstep)
753        END  
754  C      ! Recover tracer mixing ratio from "density" using predicted
755  C****6***0*********0*********0*********0*********0*********0**********72      ! "air density" (pressure thickness) at time-level n+1
756        subroutine FZPPM(IMR,JNP,NLAY,j1,DQ,WZ,P,DC,DQDT,AR,AL,A6,  
757       &                 flux,wk1,wk2,wz2,delp,KORD)      DO k = 1, nlay
758        parameter ( kmax = 150 )        DO j = 1, jnp
759        parameter ( R23 = 2./3., R3 = 1./3.)          DO i = 1, imr
760        real WZ(IMR,JNP,NLAY),P(IMR,JNP,NLAY),DC(IMR,JNP,NLAY),            q(i, j, k, ic) = dq(i, j, k, ic)/delp2(i, j, k)
761       &     wk1(IMR,*),delp(IMR,JNP,NLAY),DQ(IMR,JNP,NLAY),          END DO
762       &     DQDT(IMR,JNP,NLAY)        END DO
763  C Assuming JNP >= NLAY      END DO
764        real AR(IMR,*),AL(IMR,*),A6(IMR,*),flux(IMR,*),wk2(IMR,*),  
765       &     wz2(IMR,*)      IF (j1/=2) THEN
766  C        DO k = 1, nlay
767        JMR = JNP - 1          DO i = 1, imr
768        IMJM = IMR*JNP            ! j=1 c'est le pôle Sud, j=JNP c'est le pôle Nord
769        NLAYM1 = NLAY - 1            q(i, 2, k, ic) = q(i, 1, k, ic)
770  C            q(i, jmr, k, ic) = q(i, jmp, k, ic)
771        LMT = KORD - 3          END DO
772  C        END DO
773  C ****6***0*********0*********0*********0*********0*********0**********72      END IF
774  C Compute DC for PPM    END DO
775  C ****6***0*********0*********0*********0*********0*********0**********72  
776  C    IF (j1/=2) THEN
777        do 1000 k=1,NLAYM1      DO k = 1, nlay
778        do 1000 i=1,IMJM        DO i = 1, imr
779        DQDT(i,1,k) = P(i,1,k+1) - P(i,1,k)          w(i, 2, k) = w(i, 1, k)
780  1000  continue          w(i, jmr, k) = w(i, jnp, k)
781  C        END DO
782        DO 1220 k=2,NLAYM1      END DO
783        DO 1220 I=1,IMJM        END IF
784         c0 =  delp(i,1,k) / (delp(i,1,k-1)+delp(i,1,k)+delp(i,1,k+1))  
785         c1 = (delp(i,1,k-1)+0.5*delp(i,1,k))/(delp(i,1,k+1)+delp(i,1,k))        RETURN
786         c2 = (delp(i,1,k+1)+0.5*delp(i,1,k))/(delp(i,1,k-1)+delp(i,1,k))  END SUBROUTINE ppm3d
787        tmp = c0*(c1*DQDT(i,1,k) + c2*DQDT(i,1,k-1))  
788        Qmax = max(P(i,1,k-1),P(i,1,k),P(i,1,k+1)) - P(i,1,k)  ! ****6***0*********0*********0*********0*********0*********0**********72
789        Qmin = P(i,1,k) - min(P(i,1,k-1),P(i,1,k),P(i,1,k+1))  SUBROUTINE fzppm(imr, jnp, nlay, j1, dq, wz, p, dc, dqdt, ar, al, a6, flux, &
790        DC(i,1,k) = sign(min(abs(tmp),Qmax,Qmin), tmp)        wk1, wk2, wz2, delp, kord)
791  1220  CONTINUE    PARAMETER (kmax=150)
792          PARAMETER (r23=2./3., r3=1./3.)
793  C        REAL wz(imr, jnp, nlay), p(imr, jnp, nlay), dc(imr, jnp, nlay), &
794  C ****6***0*********0*********0*********0*********0*********0**********72      wk1(imr, *), delp(imr, jnp, nlay), dq(imr, jnp, nlay), &
795  C Loop over latitudes  (to save memory)      dqdt(imr, jnp, nlay)
796  C ****6***0*********0*********0*********0*********0*********0**********72    ! Assuming JNP >= NLAY
797  C    REAL ar(imr, *), al(imr, *), a6(imr, *), flux(imr, *), wk2(imr, *), &
798        DO 2000 j=1,JNP      wz2(imr, *)
799        if((j.eq.2 .or. j.eq.JMR) .and. j1.ne.2) goto 2000  
800  C    jmr = jnp - 1
801        DO k=1,NLAY    imjm = imr*jnp
802        DO i=1,IMR    nlaym1 = nlay - 1
803        wz2(i,k) =   WZ(i,j,k)  
804        wk1(i,k) =    P(i,j,k)    lmt = kord - 3
805        wk2(i,k) = delp(i,j,k)  
806        flux(i,k) = DC(i,j,k)  !this flux is actually the monotone slope    ! ****6***0*********0*********0*********0*********0*********0**********72
807        enddo    ! Compute DC for PPM
808        enddo    ! ****6***0*********0*********0*********0*********0*********0**********72
809  C  
810  C****6***0*********0*********0*********0*********0*********0**********72    DO k = 1, nlaym1
811  C Compute first guesses at cell interfaces      DO i = 1, imjm
812  C First guesses are required to be continuous.        dqdt(i, 1, k) = p(i, 1, k+1) - p(i, 1, k)
813  C ****6***0*********0*********0*********0*********0*********0**********72      END DO
814  C    END DO
815  C three-cell parabolic subgrid distribution at model top  
816  C two-cell parabolic with zero gradient subgrid distribution    DO k = 2, nlaym1
817  C at the surface.      DO i = 1, imjm
818  C        c0 = delp(i, 1, k)/(delp(i,1,k-1)+delp(i,1,k)+delp(i,1,k+1))
819  C First guess top edge value        c1 = (delp(i,1,k-1)+0.5*delp(i,1,k))/(delp(i,1,k+1)+delp(i,1,k))
820        DO 10 i=1,IMR        c2 = (delp(i,1,k+1)+0.5*delp(i,1,k))/(delp(i,1,k-1)+delp(i,1,k))
821  C three-cell PPM        tmp = c0*(c1*dqdt(i,1,k)+c2*dqdt(i,1,k-1))
822  C Compute a,b, and c of q = aP**2 + bP + c using cell averages and delp        qmax = max(p(i,1,k-1), p(i,1,k), p(i,1,k+1)) - p(i, 1, k)
823        a = 3.*( DQDT(i,j,2) - DQDT(i,j,1)*(wk2(i,2)+wk2(i,3))/        qmin = p(i, 1, k) - min(p(i,1,k-1), p(i,1,k), p(i,1,k+1))
824       &         (wk2(i,1)+wk2(i,2)) ) /        dc(i, 1, k) = sign(min(abs(tmp),qmax,qmin), tmp)
825       &       ( (wk2(i,2)+wk2(i,3))*(wk2(i,1)+wk2(i,2)+wk2(i,3)) )      END DO
826        b = 2.*DQDT(i,j,1)/(wk2(i,1)+wk2(i,2)) -    END DO
827       &    R23*a*(2.*wk2(i,1)+wk2(i,2))  
828        AL(i,1) =  wk1(i,1) - wk2(i,1)*(R3*a*wk2(i,1) + 0.5*b)  
829        AL(i,2) =  wk2(i,1)*(a*wk2(i,1) + b) + AL(i,1)    ! ****6***0*********0*********0*********0*********0*********0**********72
830  C    ! Loop over latitudes  (to save memory)
831  C Check if change sign    ! ****6***0*********0*********0*********0*********0*********0**********72
832        if(wk1(i,1)*AL(i,1).le.0.) then  
833                   AL(i,1) = 0.    DO j = 1, jnp
834               flux(i,1) = 0.      IF ((j==2 .OR. j==jmr) .AND. j1/=2) GO TO 2000
835          else  
836               flux(i,1) =  wk1(i,1) - AL(i,1)      DO k = 1, nlay
837          endif        DO i = 1, imr
838  10    continue          wz2(i, k) = wz(i, j, k)
839  C          wk1(i, k) = p(i, j, k)
840  C Bottom          wk2(i, k) = delp(i, j, k)
841        DO 15 i=1,IMR          flux(i, k) = dc(i, j, k) !this flux is actually the monotone slope
842  C 2-cell PPM with zero gradient right at the surface        END DO
843  C      END DO
844        fct = DQDT(i,j,NLAYM1)*wk2(i,NLAY)**2 /  
845       & ( (wk2(i,NLAY)+wk2(i,NLAYM1))*(2.*wk2(i,NLAY)+wk2(i,NLAYM1)))      ! ****6***0*********0*********0*********0*********0*********0**********72
846        AR(i,NLAY) = wk1(i,NLAY) + fct      ! Compute first guesses at cell interfaces
847        AL(i,NLAY) = wk1(i,NLAY) - (fct+fct)      ! First guesses are required to be continuous.
848        if(wk1(i,NLAY)*AR(i,NLAY).le.0.) AR(i,NLAY) = 0.      ! ****6***0*********0*********0*********0*********0*********0**********72
849        flux(i,NLAY) = AR(i,NLAY) -  wk1(i,NLAY)  
850  15    continue      ! three-cell parabolic subgrid distribution at model top
851            ! two-cell parabolic with zero gradient subgrid distribution
852  C      ! at the surface.
853  C****6***0*********0*********0*********0*********0*********0**********72  
854  C 4th order interpolation in the interior.      ! First guess top edge value
855  C****6***0*********0*********0*********0*********0*********0**********72      DO i = 1, imr
856  C        ! three-cell PPM
857        DO 14 k=3,NLAYM1        ! Compute a,b, and c of q = aP**2 + bP + c using cell averages and delp
858        DO 12 i=1,IMR        a = 3.*(dqdt(i,j,2)-dqdt(i,j,1)*(wk2(i,2)+wk2(i,3))/(wk2(i,1)+wk2(i, &
859        c1 =  DQDT(i,j,k-1)*wk2(i,k-1) / (wk2(i,k-1)+wk2(i,k))          2)))/((wk2(i,2)+wk2(i,3))*(wk2(i,1)+wk2(i,2)+wk2(i,3)))
860        c2 =  2. / (wk2(i,k-2)+wk2(i,k-1)+wk2(i,k)+wk2(i,k+1))        b = 2.*dqdt(i, j, 1)/(wk2(i,1)+wk2(i,2)) - r23*a*(2.*wk2(i,1)+wk2(i,2))
861        A1   =  (wk2(i,k-2)+wk2(i,k-1)) / (2.*wk2(i,k-1)+wk2(i,k))        al(i, 1) = wk1(i, 1) - wk2(i, 1)*(r3*a*wk2(i,1)+0.5*b)
862        A2   =  (wk2(i,k  )+wk2(i,k+1)) / (2.*wk2(i,k)+wk2(i,k-1))        al(i, 2) = wk2(i, 1)*(a*wk2(i,1)+b) + al(i, 1)
863        AL(i,k) = wk1(i,k-1) + c1 + c2 *  
864       &        ( wk2(i,k  )*(c1*(A1 - A2)+A2*flux(i,k-1)) -        ! Check if change sign
865       &          wk2(i,k-1)*A1*flux(i,k)  )        IF (wk1(i,1)*al(i,1)<=0.) THEN
866  C      print *,'AL1',i,k, AL(i,k)          al(i, 1) = 0.
867  12    CONTINUE          flux(i, 1) = 0.
 14    continue  
 C  
       do 20 i=1,IMR*NLAYM1  
       AR(i,1) = AL(i,2)  
 C      print *,'AR1',i,AR(i,1)  
 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)))  
 C      print *,'A61',i,A6(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  
 C        print *,'test2-0',i,j,wz2(i,1),wk2(i,2)  
         CP= wz2(i,1) / wk2(i,2)          
 C        print *,'testCP',CP  
         flux(i,2) = AL(i,2)+0.5*CP*(AL(i,2)-AR(i,2)-A6(i,2)*(1.+R23*CP))  
 C        print *,'test2',i, AL(i,2),AR(i,2),A6(i,2),R23  
       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)  
868        ELSE        ELSE
869        call xmist(IMR,IML,Qtmp,DC)          flux(i, 1) = wk1(i, 1) - al(i, 1)
870        DC(0) = DC(IMR)        END IF
871  C      END DO
872        if(IORD.eq.2 .or. j.le.j1vl .or. j.ge.j2vl) then  
873        DO 1408 i=1,IMR      ! Bottom
874        iu = float(i) - uc(i,j)      DO i = 1, imr
875  1408  fx1(i) = qtmp(iu) + DC(iu)*(sign(1.,uc(i,j))-uc(i,j))        ! 2-cell PPM with zero gradient right at the surface
876        else  
877        call fxppm(IMR,IML,UC(1,j),Qtmp,DC,fx1,IORD)        fct = dqdt(i, j, nlaym1)*wk2(i, nlay)**2/((wk2(i,nlay)+wk2(i, &
878        endif          nlaym1))*(2.*wk2(i,nlay)+wk2(i,nlaym1)))
879  C        ar(i, nlay) = wk1(i, nlay) + fct
880        ENDIF        al(i, nlay) = wk1(i, nlay) - (fct+fct)
881  C        IF (wk1(i,nlay)*ar(i,nlay)<=0.) ar(i, nlay) = 0.
882        DO 1506 i=1,IMR        flux(i, nlay) = ar(i, nlay) - wk1(i, nlay)
883  1506  fx1(i) = fx1(i)*xmass(i,j)      END DO
884  C  
885        goto 1309  
886  C      ! ****6***0*********0*********0*********0*********0*********0**********72
887  C ***** Conservative (flux-form) Semi-Lagrangian transport *****      ! 4th order interpolation in the interior.
888  C      ! ****6***0*********0*********0*********0*********0*********0**********72
889  2222  continue  
890  C      DO k = 3, nlaym1
891        do i=-IML,0        DO i = 1, imr
892        qtmp(i)     = q(IMR+i,j)          c1 = dqdt(i, j, k-1)*wk2(i, k-1)/(wk2(i,k-1)+wk2(i,k))
893        qtmp(IMP-i) = q(1-i,j)          c2 = 2./(wk2(i,k-2)+wk2(i,k-1)+wk2(i,k)+wk2(i,k+1))
894        enddo          a1 = (wk2(i,k-2)+wk2(i,k-1))/(2.*wk2(i,k-1)+wk2(i,k))
895  C          a2 = (wk2(i,k)+wk2(i,k+1))/(2.*wk2(i,k)+wk2(i,k-1))
896        IF(IORD.eq.1 .or. j.eq.j1. or. j.eq.j2) THEN          al(i, k) = wk1(i, k-1) + c1 + c2*(wk2(i,k)*(c1*(a1-a2)+a2*flux(i, &
897        DO 1306 i=1,IMR            k-1))-wk2(i,k-1)*a1*flux(i,k))
898        itmp = INT(uc(i,j))        END DO
899        ISAVE(i) = i - itmp      END DO
900        iu = i - uc(i,j)  
901  1306  fx1(i) = (uc(i,j) - itmp)*qtmp(iu)      DO i = 1, imr*nlaym1
902          ar(i, 1) = al(i, 2)
903        END DO
904    
905        DO i = 1, imr*nlay
906          a6(i, 1) = 3.*(wk1(i,1)+wk1(i,1)-(al(i,1)+ar(i,1)))
907        END DO
908    
909        ! ****6***0*********0*********0*********0*********0*********0**********72
910        ! Top & Bot always monotonic
911        CALL lmtppm(flux(1,1), a6(1,1), ar(1,1), al(1,1), wk1(1,1), imr, 0)
912        CALL lmtppm(flux(1,nlay), a6(1,nlay), ar(1,nlay), al(1,nlay), &
913          wk1(1,nlay), imr, 0)
914    
915        ! Interior depending on KORD
916        IF (lmt<=2) CALL lmtppm(flux(1,2), a6(1,2), ar(1,2), al(1,2), wk1(1,2), &
917          imr*(nlay-2), lmt)
918    
919        ! ****6***0*********0*********0*********0*********0*********0**********72
920    
921        DO i = 1, imr*nlaym1
922          IF (wz2(i,1)>0.) THEN
923            cm = wz2(i, 1)/wk2(i, 1)
924            flux(i, 2) = ar(i, 1) + 0.5*cm*(al(i,1)-ar(i,1)+a6(i,1)*(1.-r23*cm))
925        ELSE        ELSE
926        call xmist(IMR,IML,Qtmp,DC)          cp = wz2(i, 1)/wk2(i, 2)
927  C          flux(i, 2) = al(i, 2) + 0.5*cp*(al(i,2)-ar(i,2)-a6(i,2)*(1.+r23*cp))
928        do i=-IML,0        END IF
929        DC(i)     = DC(IMR+i)      END DO
930        DC(IMP-i) = DC(1-i)  
931        enddo      DO i = 1, imr*nlaym1
932  C        flux(i, 2) = wz2(i, 1)*flux(i, 2)
933        DO 1307 i=1,IMR      END DO
934        itmp = INT(uc(i,j))  
935        rut  = uc(i,j) - itmp      DO i = 1, imr
936        ISAVE(i) = i - itmp        dq(i, j, 1) = dq(i, j, 1) - flux(i, 2)
937        iu = i - uc(i,j)        dq(i, j, nlay) = dq(i, j, nlay) + flux(i, nlay)
938  1307  fx1(i) = rut*(qtmp(iu) + DC(iu)*(sign(1.,rut) - rut))      END DO
939        ENDIF  
940  C      DO k = 2, nlaym1
941        do 1308 i=1,IMR        DO i = 1, imr
942        IF(uc(i,j).GT.1.) then          dq(i, j, k) = dq(i, j, k) + flux(i, k) - flux(i, k+1)
943  CDIR$ NOVECTOR        END DO
944          do ist = ISAVE(i),i-1      END DO
945          fx1(i) = fx1(i) + qtmp(ist)  2000 END DO
946          enddo    RETURN
947        elseIF(uc(i,j).LT.-1.) then  END SUBROUTINE fzppm
948          do ist = i,ISAVE(i)-1  
949          fx1(i) = fx1(i) - qtmp(ist)  SUBROUTINE xtp(imr, jnp, iml, j1, j2, jn, js, pu, dq, q, uc, fx1, xmass, &
950          enddo      iord)
951  CDIR$ VECTOR    DIMENSION uc(imr, *), dc(-iml:imr+iml+1), xmass(imr, jnp), fx1(imr+1), &
952        endif      dq(imr, jnp), qtmp(-iml:imr+1+iml)
953  1308  continue    DIMENSION pu(imr, jnp), q(imr, jnp), isave(imr)
954        do i=1,IMR  
955        fx1(i) = PU(i,j)*fx1(i)    imp = imr + 1
956        enddo  
957  C    ! van Leer at high latitudes
958  C ***************************************    jvan = max(1, jnp/18)
959  C    j1vl = j1 + jvan
960  1309  fx1(IMP) = fx1(1)    j2vl = j2 - jvan
961        DO 1215 i=1,IMR  
962  1215  DQ(i,j) =  DQ(i,j) + fx1(i)-fx1(i+1)    DO j = j1, j2
963  C  
964  C ***************************************      DO i = 1, imr
965  C        qtmp(i) = q(i, j)
966  1310  continue      END DO
967        return  
968        end      IF (j>=jn .OR. j<=js) GO TO 2222
969  C      ! ************* Eulerian **********
970        subroutine fxppm(IMR,IML,UT,P,DC,flux,IORD)  
971        parameter ( R3 = 1./3., R23 = 2./3. )      qtmp(0) = q(imr, j)
972        DIMENSION UT(*),flux(*),P(-IML:IMR+IML+1),DC(-IML:IMR+IML+1)      qtmp(-1) = q(imr-1, j)
973        DIMENSION AR(0:IMR),AL(0:IMR),A6(0:IMR)      qtmp(imp) = q(1, j)
974        integer LMT      qtmp(imp+1) = q(2, j)
975  c      logical first  
976  c      data first /.true./      IF (iord==1 .OR. j==j1 .OR. j==j2) THEN
977  c      SAVE LMT        DO i = 1, imr
978  c      if(first) then          iu = float(i) - uc(i, j)
979  C          fx1(i) = qtmp(iu)
980  C correction calcul de LMT a chaque passage pour pouvoir choisir        END DO
981  c plusieurs schemas PPM pour differents traceurs      ELSE
982  c      IF (IORD.LE.0) then        CALL xmist(imr, iml, qtmp, dc)
983  c            if(IMR.GE.144) then        dc(0) = dc(imr)
984  c                  LMT = 0  
985  c            elseif(IMR.GE.72) then        IF (iord==2 .OR. j<=j1vl .OR. j>=j2vl) THEN
986  c                  LMT = 1          DO i = 1, imr
987  c            else            iu = float(i) - uc(i, j)
988  c                  LMT = 2            fx1(i) = qtmp(iu) + dc(iu)*(sign(1.,uc(i,j))-uc(i,j))
989  c            endif          END DO
 c      else  
 c            LMT = IORD - 3  
 c      endif  
 C  
       LMT = IORD - 3  
 c      write(6,*) 'PPM option in E-W direction = ', LMT  
 c      first = .false.  
 C      endif  
 C  
       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  
990        ELSE        ELSE
991        do 12 i=1,IMH          CALL fxppm(imr, iml, uc(1,j), qtmp, dc, fx1, iord)
992  C J=2        END IF
993        tmp = (8.*(p(i,3) - p(i,1)) + p(i+IMH,2) - p(i,4))*R24  
994        Pmax = max(p(i,1),p(i,2),p(i,3)) - p(i,2)      END IF
995        Pmin = p(i,2) - min(p(i,1),p(i,2),p(i,3))  
996        DC(i,2) = sign(min(abs(tmp),Pmin,Pmax),tmp)      DO i = 1, imr
997  C J=JMR        fx1(i) = fx1(i)*xmass(i, j)
998        tmp=(8.*(p(i,JNP)-p(i,JMR-1))+p(i,JMR-2)-p(i+IMH,JMR))*R24      END DO
999        Pmax = max(p(i,JMR-1),p(i,JMR),p(i,JNP)) - p(i,JMR)  
1000        Pmin = p(i,JMR) - min(p(i,JMR-1),p(i,JMR),p(i,JNP))      GO TO 1309
1001        DC(i,JMR) = sign(min(abs(tmp),Pmin,Pmax),tmp)  
1002  12    CONTINUE      ! ***** Conservative (flux-form) Semi-Lagrangian transport *****
1003        do 14 i=IMH+1,IMR  
1004  C J=2  2222 CONTINUE
1005        tmp = (8.*(p(i,3) - p(i,1)) + p(i-IMH,2) - p(i,4))*R24  
1006        Pmax = max(p(i,1),p(i,2),p(i,3)) - p(i,2)      DO i = -iml, 0
1007        Pmin = p(i,2) - min(p(i,1),p(i,2),p(i,3))        qtmp(i) = q(imr+i, j)
1008        DC(i,2) = sign(min(abs(tmp),Pmin,Pmax),tmp)        qtmp(imp-i) = q(1-i, j)
1009  C J=JMR      END DO
1010        tmp=(8.*(p(i,JNP)-p(i,JMR-1))+p(i,JMR-2)-p(i-IMH,JMR))*R24  
1011        Pmax = max(p(i,JMR-1),p(i,JMR),p(i,JNP)) - p(i,JMR)      IF (iord==1 .OR. j==j1 .OR. j==j2) THEN
1012        Pmin = p(i,JMR) - min(p(i,JMR-1),p(i,JMR),p(i,JNP))        DO i = 1, imr
1013        DC(i,JMR) = sign(min(abs(tmp),Pmin,Pmax),tmp)          itmp = int(uc(i,j))
1014  14    CONTINUE          isave(i) = i - itmp
1015  C          iu = i - uc(i, j)
1016        do 15 i=1,IJM3          fx1(i) = (uc(i,j)-itmp)*qtmp(iu)
1017        tmp = (8.*(p(i,4) - p(i,2)) + p(i,1) - p(i,5))*R24        END DO
1018        Pmax = max(p(i,2),p(i,3),p(i,4)) - p(i,3)      ELSE
1019        Pmin = p(i,3) - min(p(i,2),p(i,3),p(i,4))        CALL xmist(imr, iml, qtmp, dc)
1020        DC(i,3) = sign(min(abs(tmp),Pmin,Pmax),tmp)  
1021  15    CONTINUE        DO i = -iml, 0
1022        ENDIF          dc(i) = dc(imr+i)
1023  C          dc(imp-i) = dc(1-i)
1024        if(j1.ne.2) then        END DO
1025        do i=1,IMR  
1026        DC(i,1) = 0.        DO i = 1, imr
1027        DC(i,JNP) = 0.          itmp = int(uc(i,j))
1028        enddo          rut = uc(i, j) - itmp
1029        else          isave(i) = i - itmp
1030  C Determine slopes in polar caps for scalars!          iu = i - uc(i, j)
1031  C          fx1(i) = rut*(qtmp(iu)+dc(iu)*(sign(1.,rut)-rut))
1032        do 13 i=1,IMH        END DO
1033  C South      END IF
1034        tmp = 0.25*(p(i,2) - p(i+imh,2))  
1035        Pmax = max(p(i,2),p(i,1), p(i+imh,2)) - p(i,1)      DO i = 1, imr
1036        Pmin = p(i,1) - min(p(i,2),p(i,1), p(i+imh,2))        IF (uc(i,j)>1.) THEN
1037        DC(i,1)=sign(min(abs(tmp),Pmax,Pmin),tmp)          ! DIR$ NOVECTOR
1038  C North.          DO ist = isave(i), i - 1
1039        tmp = 0.25*(p(i+imh,JMR) - p(i,JMR))            fx1(i) = fx1(i) + qtmp(ist)
1040        Pmax = max(p(i+imh,JMR),p(i,jnp), p(i,JMR)) - p(i,JNP)          END DO
1041        Pmin = p(i,JNP) - min(p(i+imh,JMR),p(i,jnp), p(i,JMR))        ELSE IF (uc(i,j)<-1.) THEN
1042        DC(i,JNP) = sign(min(abs(tmp),Pmax,pmin),tmp)          DO ist = i, isave(i) - 1
1043  13    continue            fx1(i) = fx1(i) - qtmp(ist)
1044  C          END DO
1045        do 25 i=imh+1,IMR          ! DIR$ VECTOR
1046        DC(i,  1) =  - DC(i-imh,  1)        END IF
1047        DC(i,JNP) =  - DC(i-imh,JNP)      END DO
1048  25    continue      DO i = 1, imr
1049        endif        fx1(i) = pu(i, j)*fx1(i)
1050        return      END DO
1051        end  
1052  C      ! ***************************************
1053        subroutine fyppm(VC,P,DC,flux,IMR,JNP,j1,j2,A6,AR,AL,JORD)  
1054        parameter ( R3 = 1./3., R23 = 2./3. )  1309 fx1(imp) = fx1(1)
1055        real VC(IMR,*),flux(IMR,*),P(IMR,*),DC(IMR,*)      DO i = 1, imr
1056  C Local work arrays.        dq(i, j) = dq(i, j) + fx1(i) - fx1(i+1)
1057        real AR(IMR,JNP),AL(IMR,JNP),A6(IMR,JNP)      END DO
1058        integer LMT  
1059  c      logical first      ! ***************************************
1060  C      data first /.true./  
1061  C      SAVE LMT    END DO
1062  C    RETURN
1063        IMH = IMR / 2  END SUBROUTINE xtp
1064        JMR = JNP - 1  
1065        j11 = j1-1  SUBROUTINE fxppm(imr, iml, ut, p, dc, flux, iord)
1066        IMJM1 = IMR*(J2-J1+2)    PARAMETER (r3=1./3., r23=2./3.)
1067        len   = IMR*(J2-J1+3)    DIMENSION ut(*), flux(*), p(-iml:imr+iml+1), dc(-iml:imr+iml+1)
1068  C      if(first) then    DIMENSION ar(0:imr), al(0:imr), a6(0:imr)
1069  C      IF(JORD.LE.0) then    INTEGER lmt
1070  C            if(JMR.GE.90) then    ! logical first
1071  C                  LMT = 0    ! data first /.true./
1072  C            elseif(JMR.GE.45) then    ! SAVE LMT
1073  C                  LMT = 1    ! if(first) then
1074  C            else  
1075  C                  LMT = 2    ! correction calcul de LMT a chaque passage pour pouvoir choisir
1076  C            endif    ! plusieurs schemas PPM pour differents traceurs
1077  C      else    ! IF (IORD.LE.0) then
1078  C            LMT = JORD - 3    ! if(IMR.GE.144) then
1079  C      endif    ! LMT = 0
1080  C    ! elseif(IMR.GE.72) then
1081  C      first = .false.    ! LMT = 1
1082  C      endif    ! else
1083  C        ! LMT = 2
1084  c modifs pour pouvoir choisir plusieurs schemas PPM    ! endif
1085        LMT = JORD - 3          ! else
1086  C    ! LMT = IORD - 3
1087        DO 10 i=1,IMR*JMR            ! endif
1088        AL(i,2) = 0.5*(p(i,1)+p(i,2)) + (DC(i,1) - DC(i,2))*R3  
1089        AR(i,1) = AL(i,2)    lmt = iord - 3
1090  10    CONTINUE  
1091  C    DO i = 1, imr
1092  CPoles:      al(i) = 0.5*(p(i-1)+p(i)) + (dc(i-1)-dc(i))*r3
1093  C    END DO
1094        DO i=1,IMH  
1095        AL(i,1) = AL(i+IMH,2)    DO i = 1, imr - 1
1096        AL(i+IMH,1) = AL(i,2)      ar(i) = al(i+1)
1097  C    END DO
1098        AR(i,JNP) = AR(i+IMH,JMR)    ar(imr) = al(1)
1099        AR(i+IMH,JNP) = AR(i,JMR)  
1100        ENDDO    DO i = 1, imr
1101        a6(i) = 3.*(p(i)+p(i)-(al(i)+ar(i)))
1102  ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc    END DO
1103  c   Rajout pour LMDZ.3.3  
1104  cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc    IF (lmt<=2) CALL lmtppm(dc(1), a6(1), ar(1), al(1), p(1), imr, lmt)
1105        AR(IMR,1)=AL(1,1)  
1106        AR(IMR,JNP)=AL(1,JNP)    al(0) = al(imr)
1107  ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc    ar(0) = ar(imr)
1108            a6(0) = a6(imr)
1109              
1110        do 30 i=1,len    DO i = 1, imr
1111  30    A6(i,j11) = 3.*(p(i,j11)+p(i,j11)  - (AL(i,j11)+AR(i,j11)))      IF (ut(i)>0.) THEN
1112  C        flux(i) = ar(i-1) + 0.5*ut(i)*(al(i-1)-ar(i-1)+a6(i-1)*(1.-r23*ut(i)))
1113        if(LMT.le.2) call lmtppm(DC(1,j11),A6(1,j11),AR(1,j11)      ELSE
1114       &                       ,AL(1,j11),P(1,j11),len,LMT)        flux(i) = al(i) - 0.5*ut(i)*(ar(i)-al(i)+a6(i)*(1.+r23*ut(i)))
1115  C      END IF
1116          END DO
1117        DO 140 i=1,IMJM1    RETURN
1118        IF(VC(i,j1).GT.0.) then  END SUBROUTINE fxppm
1119        flux(i,j1) = AR(i,j11) + 0.5*VC(i,j1)*(AL(i,j11) - AR(i,j11) +  
1120       &                         A6(i,j11)*(1.-R23*VC(i,j1)) )  SUBROUTINE xmist(imr, iml, p, dc)
1121        else    PARAMETER (r24=1./24.)
1122        flux(i,j1) = AL(i,j1) - 0.5*VC(i,j1)*(AR(i,j1) - AL(i,j1) +    DIMENSION p(-iml:imr+1+iml), dc(-iml:imr+1+iml)
1123       &                        A6(i,j1)*(1.+R23*VC(i,j1)))  
1124        endif    DO i = 1, imr
1125  140   continue      tmp = r24*(8.*(p(i+1)-p(i-1))+p(i-2)-p(i+2))
1126        return      pmax = max(p(i-1), p(i), p(i+1)) - p(i)
1127        end      pmin = p(i) - min(p(i-1), p(i), p(i+1))
1128  C      dc(i) = sign(min(abs(tmp),pmax,pmin), tmp)
1129          subroutine yadv(IMR,JNP,j1,j2,p,VA,ady,wk,IAD)    END DO
1130          REAL p(IMR,JNP),ady(IMR,JNP),VA(IMR,JNP)    RETURN
1131          REAL WK(IMR,-1:JNP+2)  END SUBROUTINE xmist
1132  C  
1133          JMR = JNP-1  SUBROUTINE ytp(imr, jnp, j1, j2, acosp, rcap, dq, p, vc, dc2, ymass, fx, a6, &
1134          IMH = IMR/2      ar, al, jord)
1135          do j=1,JNP    DIMENSION p(imr, jnp), vc(imr, jnp), ymass(imr, jnp), dc2(imr, jnp), &
1136          do i=1,IMR      dq(imr, jnp), acosp(jnp)
1137          wk(i,j) = p(i,j)    ! Work array
1138          enddo    DIMENSION fx(imr, jnp), ar(imr, jnp), al(imr, jnp), a6(imr, jnp)
1139          enddo  
1140  C Poles:    jmr = jnp - 1
1141          do i=1,IMH    len = imr*(j2-j1+2)
1142          wk(i,   -1) = p(i+IMH,3)  
1143          wk(i+IMH,-1) = p(i,3)    IF (jord==1) THEN
1144          wk(i,    0) = p(i+IMH,2)      DO i = 1, len
1145          wk(i+IMH,0) = p(i,2)        jt = float(j1) - vc(i, j1)
1146          wk(i,JNP+1) = p(i+IMH,JMR)        fx(i, j1) = p(i, jt)
1147          wk(i+IMH,JNP+1) = p(i,JMR)      END DO
1148          wk(i,JNP+2) = p(i+IMH,JNP-2)    ELSE
1149          wk(i+IMH,JNP+2) = p(i,JNP-2)  
1150          enddo      CALL ymist(imr, jnp, j1, p, dc2, 4)
1151  c        write(*,*) 'toto 1'  
1152  C --------------------------------      IF (jord<=0 .OR. jord>=3) THEN
1153        IF(IAD.eq.2) then  
1154        do j=j1-1,j2+1        CALL fyppm(vc, p, dc2, fx, imr, jnp, j1, j2, a6, ar, al, jord)
1155        do i=1,IMR  
1156  c      write(*,*) 'avt NINT','i=',i,'j=',j      ELSE
1157        JP = NINT(VA(i,j))              DO i = 1, len
1158        rv = JP - VA(i,j)          jt = float(j1) - vc(i, j1)
1159  c      write(*,*) 'VA=',VA(i,j), 'JP1=',JP,'rv=',rv          fx(i, j1) = p(i, jt) + (sign(1.,vc(i,j1))-vc(i,j1))*dc2(i, jt)
1160        JP = j - JP        END DO
1161  c      write(*,*) 'JP2=',JP      END IF
1162        a1 = 0.5*(wk(i,jp+1)+wk(i,jp-1)) - wk(i,jp)    END IF
1163        b1 = 0.5*(wk(i,jp+1)-wk(i,jp-1))  
1164  c      write(*,*) 'a1=',a1,'b1=',b1    DO i = 1, len
1165        ady(i,j) = wk(i,jp) + rv*(a1*rv + b1) - wk(i,j)      fx(i, j1) = fx(i, j1)*ymass(i, j1)
1166        enddo    END DO
1167        enddo  
1168  c      write(*,*) 'toto 2'    DO j = j1, j2
1169  C      DO i = 1, imr
1170        ELSEIF(IAD.eq.1) then        dq(i, j) = dq(i, j) + (fx(i,j)-fx(i,j+1))*acosp(j)
1171          do j=j1-1,j2+1      END DO
1172        do i=1,imr    END DO
1173        JP = float(j)-VA(i,j)  
1174        ady(i,j) = VA(i,j)*(wk(i,jp)-wk(i,jp+1))    ! Poles
1175        enddo    sum1 = fx(imr, j1)
1176        enddo    sum2 = fx(imr, j2+1)
1177        ENDIF    DO i = 1, imr - 1
1178  C      sum1 = sum1 + fx(i, j1)
1179          if(j1.ne.2) then      sum2 = sum2 + fx(i, j2+1)
1180          sum1 = 0.    END DO
1181          sum2 = 0.  
1182        do i=1,imr    sum1 = dq(1, 1) - sum1*rcap
1183        sum1 = sum1 + ady(i,2)    sum2 = dq(1, jnp) + sum2*rcap
1184        sum2 = sum2 + ady(i,JMR)    DO i = 1, imr
1185        enddo      dq(i, 1) = sum1
1186          sum1 = sum1 / IMR      dq(i, jnp) = sum2
1187          sum2 = sum2 / IMR    END DO
1188  C  
1189        do i=1,imr    IF (j1/=2) THEN
1190        ady(i,  2) =  sum1      DO i = 1, imr
1191        ady(i,JMR) =  sum2        dq(i, 2) = sum1
1192        ady(i,  1) =  sum1        dq(i, jmr) = sum2
1193        ady(i,JNP) =  sum2      END DO
1194        enddo    END IF
1195          else  
1196  C Poles:    RETURN
1197          sum1 = 0.  END SUBROUTINE ytp
1198          sum2 = 0.  
1199        do i=1,imr  SUBROUTINE ymist(imr, jnp, j1, p, dc, id)
1200        sum1 = sum1 + ady(i,1)    PARAMETER (r24=1./24.)
1201        sum2 = sum2 + ady(i,JNP)    DIMENSION p(imr, jnp), dc(imr, jnp)
1202        enddo  
1203          sum1 = sum1 / IMR    imh = imr/2
1204          sum2 = sum2 / IMR    jmr = jnp - 1
1205  C    ijm3 = imr*(jmr-3)
1206        do i=1,imr  
1207        ady(i,  1) =  sum1    IF (id==2) THEN
1208        ady(i,JNP) =  sum2      DO i = 1, imr*(jmr-1)
1209        enddo        tmp = 0.25*(p(i,3)-p(i,1))
1210          endif        pmax = max(p(i,1), p(i,2), p(i,3)) - p(i, 2)
1211  C        pmin = p(i, 2) - min(p(i,1), p(i,2), p(i,3))
1212          return        dc(i, 2) = sign(min(abs(tmp),pmin,pmax), tmp)
1213          end      END DO
1214  C    ELSE
1215          subroutine xadv(IMR,JNP,j1,j2,p,UA,JS,JN,IML,adx,IAD)      DO i = 1, imh
1216          REAL p(IMR,JNP),adx(IMR,JNP),qtmp(-IMR:IMR+IMR),UA(IMR,JNP)        ! J=2
1217  C        tmp = (8.*(p(i,3)-p(i,1))+p(i+imh,2)-p(i,4))*r24
1218          JMR = JNP-1        pmax = max(p(i,1), p(i,2), p(i,3)) - p(i, 2)
1219        do 1309 j=j1,j2        pmin = p(i, 2) - min(p(i,1), p(i,2), p(i,3))
1220        if(J.GT.JS  .and. J.LT.JN) GO TO 1309        dc(i, 2) = sign(min(abs(tmp),pmin,pmax), tmp)
1221  C        ! J=JMR
1222        do i=1,IMR        tmp = (8.*(p(i,jnp)-p(i,jmr-1))+p(i,jmr-2)-p(i+imh,jmr))*r24
1223        qtmp(i) = p(i,j)        pmax = max(p(i,jmr-1), p(i,jmr), p(i,jnp)) - p(i, jmr)
1224        enddo        pmin = p(i, jmr) - min(p(i,jmr-1), p(i,jmr), p(i,jnp))
1225  C        dc(i, jmr) = sign(min(abs(tmp),pmin,pmax), tmp)
1226        do i=-IML,0      END DO
1227        qtmp(i)       = p(IMR+i,j)      DO i = imh + 1, imr
1228        qtmp(IMR+1-i) = p(1-i,j)        ! J=2
1229        enddo        tmp = (8.*(p(i,3)-p(i,1))+p(i-imh,2)-p(i,4))*r24
1230  C        pmax = max(p(i,1), p(i,2), p(i,3)) - p(i, 2)
1231        IF(IAD.eq.2) THEN        pmin = p(i, 2) - min(p(i,1), p(i,2), p(i,3))
1232        DO i=1,IMR        dc(i, 2) = sign(min(abs(tmp),pmin,pmax), tmp)
1233        IP = NINT(UA(i,j))        ! J=JMR
1234        ru = IP - UA(i,j)        tmp = (8.*(p(i,jnp)-p(i,jmr-1))+p(i,jmr-2)-p(i-imh,jmr))*r24
1235        IP = i - IP        pmax = max(p(i,jmr-1), p(i,jmr), p(i,jnp)) - p(i, jmr)
1236        a1 = 0.5*(qtmp(ip+1)+qtmp(ip-1)) - qtmp(ip)        pmin = p(i, jmr) - min(p(i,jmr-1), p(i,jmr), p(i,jnp))
1237        b1 = 0.5*(qtmp(ip+1)-qtmp(ip-1))        dc(i, jmr) = sign(min(abs(tmp),pmin,pmax), tmp)
1238        adx(i,j) = qtmp(ip) + ru*(a1*ru + b1)      END DO
1239        enddo  
1240        ELSEIF(IAD.eq.1) then      DO i = 1, ijm3
1241        DO i=1,IMR        tmp = (8.*(p(i,4)-p(i,2))+p(i,1)-p(i,5))*r24
1242        iu = UA(i,j)        pmax = max(p(i,2), p(i,3), p(i,4)) - p(i, 3)
1243        ru = UA(i,j) - iu        pmin = p(i, 3) - min(p(i,2), p(i,3), p(i,4))
1244        iiu = i-iu        dc(i, 3) = sign(min(abs(tmp),pmin,pmax), tmp)
1245        if(UA(i,j).GE.0.) then      END DO
1246        adx(i,j) = qtmp(iiu)+ru*(qtmp(iiu-1)-qtmp(iiu))    END IF
1247        else  
1248        adx(i,j) = qtmp(iiu)+ru*(qtmp(iiu)-qtmp(iiu+1))    IF (j1/=2) THEN
1249        endif      DO i = 1, imr
1250        enddo        dc(i, 1) = 0.
1251        ENDIF        dc(i, jnp) = 0.
1252  C      END DO
1253        do i=1,IMR    ELSE
1254        adx(i,j) = adx(i,j) - p(i,j)      ! Determine slopes in polar caps for scalars!
1255        enddo  
1256  1309  continue      DO i = 1, imh
1257  C        ! South
1258  C Eulerian upwind        tmp = 0.25*(p(i,2)-p(i+imh,2))
1259  C        pmax = max(p(i,2), p(i,1), p(i+imh,2)) - p(i, 1)
1260        do j=JS+1,JN-1        pmin = p(i, 1) - min(p(i,2), p(i,1), p(i+imh,2))
1261  C        dc(i, 1) = sign(min(abs(tmp),pmax,pmin), tmp)
1262        do i=1,IMR        ! North.
1263        qtmp(i) = p(i,j)        tmp = 0.25*(p(i+imh,jmr)-p(i,jmr))
1264        enddo        pmax = max(p(i+imh,jmr), p(i,jnp), p(i,jmr)) - p(i, jnp)
1265  C        pmin = p(i, jnp) - min(p(i+imh,jmr), p(i,jnp), p(i,jmr))
1266        qtmp(0)     = p(IMR,J)        dc(i, jnp) = sign(min(abs(tmp),pmax,pmin), tmp)
1267        qtmp(IMR+1) = p(1,J)      END DO
1268  C  
1269        IF(IAD.eq.2) THEN      DO i = imh + 1, imr
1270        qtmp(-1)     = p(IMR-1,J)        dc(i, 1) = -dc(i-imh, 1)
1271        qtmp(IMR+2) = p(2,J)        dc(i, jnp) = -dc(i-imh, jnp)
1272        do i=1,imr      END DO
1273        IP = NINT(UA(i,j))    END IF
1274        ru = IP - UA(i,j)    RETURN
1275        IP = i - IP  END SUBROUTINE ymist
1276        a1 = 0.5*(qtmp(ip+1)+qtmp(ip-1)) - qtmp(ip)  
1277        b1 = 0.5*(qtmp(ip+1)-qtmp(ip-1))  SUBROUTINE fyppm(vc, p, dc, flux, imr, jnp, j1, j2, a6, ar, al, jord)
1278        adx(i,j) = qtmp(ip)- p(i,j) + ru*(a1*ru + b1)    PARAMETER (r3=1./3., r23=2./3.)
1279        enddo    REAL vc(imr, *), flux(imr, *), p(imr, *), dc(imr, *)
1280        ELSEIF(IAD.eq.1) then    ! Local work arrays.
1281  C 1st order    REAL ar(imr, jnp), al(imr, jnp), a6(imr, jnp)
1282        DO i=1,IMR    INTEGER lmt
1283        IP = i - UA(i,j)    ! logical first
1284        adx(i,j) = UA(i,j)*(qtmp(ip)-qtmp(ip+1))    ! data first /.true./
1285        enddo    ! SAVE LMT
1286        ENDIF  
1287        enddo    imh = imr/2
1288  C    jmr = jnp - 1
1289          if(j1.ne.2) then    j11 = j1 - 1
1290        do i=1,IMR    imjm1 = imr*(j2-j1+2)
1291        adx(i,  2) = 0.    len = imr*(j2-j1+3)
1292        adx(i,JMR) = 0.    ! if(first) then
1293        enddo    ! IF(JORD.LE.0) then
1294          endif    ! if(JMR.GE.90) then
1295  C set cross term due to x-adv at the poles to zero.    ! LMT = 0
1296        do i=1,IMR    ! elseif(JMR.GE.45) then
1297        adx(i,  1) = 0.    ! LMT = 1
1298        adx(i,JNP) = 0.    ! else
1299        enddo    ! LMT = 2
1300          return    ! endif
1301          end    ! else
1302  C    ! LMT = JORD - 3
1303        subroutine lmtppm(DC,A6,AR,AL,P,IM,LMT)    ! endif
1304  C  
1305  C A6 =  CURVATURE OF THE TEST PARABOLA    ! first = .false.
1306  C AR =  RIGHT EDGE VALUE OF THE TEST PARABOLA    ! endif
1307  C AL =  LEFT  EDGE VALUE OF THE TEST PARABOLA  
1308  C DC =  0.5 * MISMATCH    ! modifs pour pouvoir choisir plusieurs schemas PPM
1309  C P  =  CELL-AVERAGED VALUE    lmt = jord - 3
1310  C IM =  VECTOR LENGTH  
1311  C    DO i = 1, imr*jmr
1312  C OPTIONS:      al(i, 2) = 0.5*(p(i,1)+p(i,2)) + (dc(i,1)-dc(i,2))*r3
1313  C      ar(i, 1) = al(i, 2)
1314  C LMT = 0: FULL MONOTONICITY    END DO
1315  C LMT = 1: SEMI-MONOTONIC CONSTRAINT (NO UNDERSHOOTS)  
1316  C LMT = 2: POSITIVE-DEFINITE CONSTRAINT    ! Poles:
1317  C  
1318        parameter ( R12 = 1./12. )    DO i = 1, imh
1319        dimension A6(IM),AR(IM),AL(IM),P(IM),DC(IM)      al(i, 1) = al(i+imh, 2)
1320  C      al(i+imh, 1) = al(i, 2)
1321        if(LMT.eq.0) then  
1322  C Full constraint      ar(i, jnp) = ar(i+imh, jmr)
1323        do 100 i=1,IM      ar(i+imh, jnp) = ar(i, jmr)
1324        if(DC(i).eq.0.) then    END DO
1325              AR(i) = p(i)  
1326              AL(i) = p(i)    ! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1327              A6(i) = 0.    ! Rajout pour LMDZ.3.3
1328        else    ! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1329        da1  = AR(i) - AL(i)    ar(imr, 1) = al(1, 1)
1330        da2  = da1**2    ar(imr, jnp) = al(1, jnp)
1331        A6DA = A6(i)*da1    ! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1332        if(A6DA .lt. -da2) then  
1333              A6(i) = 3.*(AL(i)-p(i))  
1334              AR(i) = AL(i) - A6(i)    DO i = 1, len
1335        elseif(A6DA .gt. da2) then      a6(i, j11) = 3.*(p(i,j11)+p(i,j11)-(al(i,j11)+ar(i,j11)))
1336              A6(i) = 3.*(AR(i)-p(i))    END DO
1337              AL(i) = AR(i) - A6(i)  
1338        endif    IF (lmt<=2) CALL lmtppm(dc(1,j11), a6(1,j11), ar(1,j11), al(1,j11), &
1339        endif      p(1,j11), len, lmt)
1340  100   continue  
1341        elseif(LMT.eq.1) then  
1342  C Semi-monotonic constraint    DO i = 1, imjm1
1343        do 150 i=1,IM      IF (vc(i,j1)>0.) THEN
1344        if(abs(AR(i)-AL(i)) .GE. -A6(i)) go to 150        flux(i, j1) = ar(i, j11) + 0.5*vc(i, j1)*(al(i,j11)-ar(i,j11)+a6(i,j11) &
1345        if(p(i).lt.AR(i) .and. p(i).lt.AL(i)) then          *(1.-r23*vc(i,j1)))
1346              AR(i) = p(i)      ELSE
1347              AL(i) = p(i)        flux(i, j1) = al(i, j1) - 0.5*vc(i, j1)*(ar(i,j1)-al(i,j1)+a6(i,j1)*(1. &
1348              A6(i) = 0.          +r23*vc(i,j1)))
1349        elseif(AR(i) .gt. AL(i)) then      END IF
1350              A6(i) = 3.*(AL(i)-p(i))    END DO
1351              AR(i) = AL(i) - A6(i)    RETURN
1352        else  END SUBROUTINE fyppm
1353              A6(i) = 3.*(AR(i)-p(i))  
1354              AL(i) = AR(i) - A6(i)  SUBROUTINE yadv(imr, jnp, j1, j2, p, va, ady, wk, iad)
1355        endif    REAL p(imr, jnp), ady(imr, jnp), va(imr, jnp)
1356  150   continue    REAL wk(imr, -1:jnp+2)
1357        elseif(LMT.eq.2) then  
1358        do 250 i=1,IM    jmr = jnp - 1
1359        if(abs(AR(i)-AL(i)) .GE. -A6(i)) go to 250    imh = imr/2
1360        fmin = p(i) + 0.25*(AR(i)-AL(i))**2/A6(i) + A6(i)*R12    DO j = 1, jnp
1361        if(fmin.ge.0.) go to 250      DO i = 1, imr
1362        if(p(i).lt.AR(i) .and. p(i).lt.AL(i)) then        wk(i, j) = p(i, j)
1363              AR(i) = p(i)      END DO
1364              AL(i) = p(i)    END DO
1365              A6(i) = 0.    ! Poles:
1366        elseif(AR(i) .gt. AL(i)) then    DO i = 1, imh
1367              A6(i) = 3.*(AL(i)-p(i))      wk(i, -1) = p(i+imh, 3)
1368              AR(i) = AL(i) - A6(i)      wk(i+imh, -1) = p(i, 3)
1369        else      wk(i, 0) = p(i+imh, 2)
1370              A6(i) = 3.*(AR(i)-p(i))      wk(i+imh, 0) = p(i, 2)
1371              AL(i) = AR(i) - A6(i)      wk(i, jnp+1) = p(i+imh, jmr)
1372        endif      wk(i+imh, jnp+1) = p(i, jmr)
1373  250   continue      wk(i, jnp+2) = p(i+imh, jnp-2)
1374        endif      wk(i+imh, jnp+2) = p(i, jnp-2)
1375        return    END DO
1376        end  
1377  C    IF (iad==2) THEN
1378        subroutine A2C(U,V,IMR,JMR,j1,j2,CRX,CRY,dtdx5,DTDY5)      DO j = j1 - 1, j2 + 1
1379        dimension U(IMR,*),V(IMR,*),CRX(IMR,*),CRY(IMR,*),DTDX5(*)        DO i = 1, imr
1380  C          jp = nint(va(i,j))
1381        do 35 j=j1,j2          rv = jp - va(i, j)
1382        do 35 i=2,IMR          jp = j - jp
1383  35    CRX(i,J) = dtdx5(j)*(U(i,j)+U(i-1,j))          a1 = 0.5*(wk(i,jp+1)+wk(i,jp-1)) - wk(i, jp)
1384  C          b1 = 0.5*(wk(i,jp+1)-wk(i,jp-1))
1385        do 45 j=j1,j2          ady(i, j) = wk(i, jp) + rv*(a1*rv+b1) - wk(i, j)
1386  45    CRX(1,J) = dtdx5(j)*(U(1,j)+U(IMR,j))        END DO
1387  C      END DO
1388        do 55 i=1,IMR*JMR  
1389  55    CRY(i,2) = DTDY5*(V(i,2)+V(i,1))    ELSE IF (iad==1) THEN
1390        return      DO j = j1 - 1, j2 + 1
1391        end        DO i = 1, imr
1392  C          jp = float(j) - va(i, j)
1393        subroutine cosa(cosp,cose,JNP,PI,DP)          ady(i, j) = va(i, j)*(wk(i,jp)-wk(i,jp+1))
1394        dimension cosp(*),cose(*)        END DO
1395        JMR = JNP-1      END DO
1396        do 55 j=2,JNP    END IF
1397          ph5  =  -0.5*PI + (FLOAT(J-1)-0.5)*DP  
1398  55      cose(j) = cos(ph5)    IF (j1/=2) THEN
1399  C      sum1 = 0.
1400        JEQ = (JNP+1) / 2      sum2 = 0.
1401        if(JMR .eq. 2*(JMR/2) ) then      DO i = 1, imr
1402        do j=JNP, JEQ+1, -1        sum1 = sum1 + ady(i, 2)
1403         cose(j) =  cose(JNP+2-j)        sum2 = sum2 + ady(i, jmr)
1404        enddo      END DO
1405        else      sum1 = sum1/imr
1406  C cell edge at equator.      sum2 = sum2/imr
1407         cose(JEQ+1) =  1.  
1408        do j=JNP, JEQ+2, -1      DO i = 1, imr
1409         cose(j) =  cose(JNP+2-j)        ady(i, 2) = sum1
1410         enddo        ady(i, jmr) = sum2
1411        endif        ady(i, 1) = sum1
1412  C        ady(i, jnp) = sum2
1413        do 66 j=2,JMR      END DO
1414  66    cosp(j) = 0.5*(cose(j)+cose(j+1))    ELSE
1415        cosp(1) = 0.      ! Poles:
1416        cosp(JNP) = 0.      sum1 = 0.
1417        return      sum2 = 0.
1418        end      DO i = 1, imr
1419  C        sum1 = sum1 + ady(i, 1)
1420        subroutine cosc(cosp,cose,JNP,PI,DP)        sum2 = sum2 + ady(i, jnp)
1421        dimension cosp(*),cose(*)      END DO
1422  C      sum1 = sum1/imr
1423        phi = -0.5*PI      sum2 = sum2/imr
1424        do 55 j=2,JNP-1  
1425        phi  =  phi + DP      DO i = 1, imr
1426  55    cosp(j) = cos(phi)        ady(i, 1) = sum1
1427          cosp(  1) = 0.        ady(i, jnp) = sum2
1428          cosp(JNP) = 0.      END DO
1429  C    END IF
1430        do 66 j=2,JNP  
1431          cose(j) = 0.5*(cosp(j)+cosp(j-1))    RETURN
1432  66    CONTINUE  END SUBROUTINE yadv
1433  C  
1434        do 77 j=2,JNP-1  SUBROUTINE xadv(imr, jnp, j1, j2, p, ua, js, jn, iml, adx, iad)
1435         cosp(j) = 0.5*(cose(j)+cose(j+1))    REAL p(imr, jnp), adx(imr, jnp), qtmp(-imr:imr+imr), ua(imr, jnp)
1436  77    CONTINUE  
1437        return    jmr = jnp - 1
1438        end    DO j = j1, j2
1439  C      IF (j>js .AND. j<jn) GO TO 1309
1440        SUBROUTINE qckxyz (Q,qtmp,IMR,JNP,NLAY,j1,j2,cosp,acosp,  
1441       &                   cross,IC,NSTEP)      DO i = 1, imr
1442  C        qtmp(i) = p(i, j)
1443        parameter( tiny = 1.E-60 )      END DO
1444        DIMENSION Q(IMR,JNP,NLAY),qtmp(IMR,JNP),cosp(*),acosp(*)  
1445        logical cross      DO i = -iml, 0
1446  C        qtmp(i) = p(imr+i, j)
1447        NLAYM1 = NLAY-1        qtmp(imr+1-i) = p(1-i, j)
1448        len = IMR*(j2-j1+1)      END DO
1449        ip = 0  
1450  C      IF (iad==2) THEN
1451  C Top layer        DO i = 1, imr
1452        L = 1          ip = nint(ua(i,j))
1453          icr = 1          ru = ip - ua(i, j)
1454        call filns(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,ipy,tiny)          ip = i - ip
1455        if(ipy.eq.0) goto 50          a1 = 0.5*(qtmp(ip+1)+qtmp(ip-1)) - qtmp(ip)
1456        call filew(q(1,1,L),qtmp,IMR,JNP,j1,j2,ipx,tiny)          b1 = 0.5*(qtmp(ip+1)-qtmp(ip-1))
1457        if(ipx.eq.0) goto 50          adx(i, j) = qtmp(ip) + ru*(a1*ru+b1)
1458  C        END DO
1459        if(cross) then      ELSE IF (iad==1) THEN
1460        call filcr(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,icr,tiny)        DO i = 1, imr
1461        endif          iu = ua(i, j)
1462        if(icr.eq.0) goto 50          ru = ua(i, j) - iu
1463  C          iiu = i - iu
1464  C Vertical filling...          IF (ua(i,j)>=0.) THEN
1465        do i=1,len            adx(i, j) = qtmp(iiu) + ru*(qtmp(iiu-1)-qtmp(iiu))
1466        IF( Q(i,j1,1).LT.0.) THEN          ELSE
1467        ip = ip + 1            adx(i, j) = qtmp(iiu) + ru*(qtmp(iiu)-qtmp(iiu+1))
1468            Q(i,j1,2) = Q(i,j1,2) + Q(i,j1,1)          END IF
1469            Q(i,j1,1) = 0.        END DO
1470        endif      END IF
1471        enddo  
1472  C      DO i = 1, imr
1473  50    continue        adx(i, j) = adx(i, j) - p(i, j)
1474        DO 225 L = 2,NLAYM1      END DO
1475        icr = 1  1309 END DO
1476  C  
1477        call filns(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,ipy,tiny)    ! Eulerian upwind
1478        if(ipy.eq.0) goto 225  
1479        call filew(q(1,1,L),qtmp,IMR,JNP,j1,j2,ipx,tiny)    DO j = js + 1, jn - 1
1480        if(ipx.eq.0) go to 225  
1481        if(cross) then      DO i = 1, imr
1482        call filcr(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,icr,tiny)        qtmp(i) = p(i, j)
1483        endif      END DO
1484        if(icr.eq.0) goto 225  
1485  C      qtmp(0) = p(imr, j)
1486        do i=1,len      qtmp(imr+1) = p(1, j)
1487        IF( Q(I,j1,L).LT.0.) THEN  
1488  C      IF (iad==2) THEN
1489          qtmp(-1) = p(imr-1, j)
1490          qtmp(imr+2) = p(2, j)
1491          DO i = 1, imr
1492            ip = nint(ua(i,j))
1493            ru = ip - ua(i, j)
1494            ip = i - ip
1495            a1 = 0.5*(qtmp(ip+1)+qtmp(ip-1)) - qtmp(ip)
1496            b1 = 0.5*(qtmp(ip+1)-qtmp(ip-1))
1497            adx(i, j) = qtmp(ip) - p(i, j) + ru*(a1*ru+b1)
1498          END DO
1499        ELSE IF (iad==1) THEN
1500          ! 1st order
1501          DO i = 1, imr
1502            ip = i - ua(i, j)
1503            adx(i, j) = ua(i, j)*(qtmp(ip)-qtmp(ip+1))
1504          END DO
1505        END IF
1506      END DO
1507    
1508      IF (j1/=2) THEN
1509        DO i = 1, imr
1510          adx(i, 2) = 0.
1511          adx(i, jmr) = 0.
1512        END DO
1513      END IF
1514      ! set cross term due to x-adv at the poles to zero.
1515      DO i = 1, imr
1516        adx(i, 1) = 0.
1517        adx(i, jnp) = 0.
1518      END DO
1519      RETURN
1520    END SUBROUTINE xadv
1521    
1522    SUBROUTINE lmtppm(dc, a6, ar, al, p, im, lmt)
1523    
1524      ! A6 =  CURVATURE OF THE TEST PARABOLA
1525      ! AR =  RIGHT EDGE VALUE OF THE TEST PARABOLA
1526      ! AL =  LEFT  EDGE VALUE OF THE TEST PARABOLA
1527      ! DC =  0.5 * MISMATCH
1528      ! P  =  CELL-AVERAGED VALUE
1529      ! IM =  VECTOR LENGTH
1530    
1531      ! OPTIONS:
1532    
1533      ! LMT = 0: FULL MONOTONICITY
1534      ! LMT = 1: SEMI-MONOTONIC CONSTRAINT (NO UNDERSHOOTS)
1535      ! LMT = 2: POSITIVE-DEFINITE CONSTRAINT
1536    
1537      PARAMETER (r12=1./12.)
1538      DIMENSION a6(im), ar(im), al(im), p(im), dc(im)
1539    
1540      IF (lmt==0) THEN
1541        ! Full constraint
1542        DO i = 1, im
1543          IF (dc(i)==0.) THEN
1544            ar(i) = p(i)
1545            al(i) = p(i)
1546            a6(i) = 0.
1547          ELSE
1548            da1 = ar(i) - al(i)
1549            da2 = da1**2
1550            a6da = a6(i)*da1
1551            IF (a6da<-da2) THEN
1552              a6(i) = 3.*(al(i)-p(i))
1553              ar(i) = al(i) - a6(i)
1554            ELSE IF (a6da>da2) THEN
1555              a6(i) = 3.*(ar(i)-p(i))
1556              al(i) = ar(i) - a6(i)
1557            END IF
1558          END IF
1559        END DO
1560      ELSE IF (lmt==1) THEN
1561        ! Semi-monotonic constraint
1562        DO i = 1, im
1563          IF (abs(ar(i)-al(i))>=-a6(i)) GO TO 150
1564          IF (p(i)<ar(i) .AND. p(i)<al(i)) THEN
1565            ar(i) = p(i)
1566            al(i) = p(i)
1567            a6(i) = 0.
1568          ELSE IF (ar(i)>al(i)) THEN
1569            a6(i) = 3.*(al(i)-p(i))
1570            ar(i) = al(i) - a6(i)
1571          ELSE
1572            a6(i) = 3.*(ar(i)-p(i))
1573            al(i) = ar(i) - a6(i)
1574          END IF
1575    150 END DO
1576      ELSE IF (lmt==2) THEN
1577        DO i = 1, im
1578          IF (abs(ar(i)-al(i))>=-a6(i)) GO TO 250
1579          fmin = p(i) + 0.25*(ar(i)-al(i))**2/a6(i) + a6(i)*r12
1580          IF (fmin>=0.) GO TO 250
1581          IF (p(i)<ar(i) .AND. p(i)<al(i)) THEN
1582            ar(i) = p(i)
1583            al(i) = p(i)
1584            a6(i) = 0.
1585          ELSE IF (ar(i)>al(i)) THEN
1586            a6(i) = 3.*(al(i)-p(i))
1587            ar(i) = al(i) - a6(i)
1588          ELSE
1589            a6(i) = 3.*(ar(i)-p(i))
1590            al(i) = ar(i) - a6(i)
1591          END IF
1592    250 END DO
1593      END IF
1594      RETURN
1595    END SUBROUTINE lmtppm
1596    
1597    SUBROUTINE a2c(u, v, imr, jmr, j1, j2, crx, cry, dtdx5, dtdy5)
1598      DIMENSION u(imr, *), v(imr, *), crx(imr, *), cry(imr, *), dtdx5(*)
1599    
1600      DO j = j1, j2
1601        DO i = 2, imr
1602          crx(i, j) = dtdx5(j)*(u(i,j)+u(i-1,j))
1603        END DO
1604      END DO
1605    
1606      DO j = j1, j2
1607        crx(1, j) = dtdx5(j)*(u(1,j)+u(imr,j))
1608      END DO
1609    
1610      DO i = 1, imr*jmr
1611        cry(i, 2) = dtdy5*(v(i,2)+v(i,1))
1612      END DO
1613      RETURN
1614    END SUBROUTINE a2c
1615    
1616    SUBROUTINE cosa(cosp, cose, jnp, pi, dp)
1617      DIMENSION cosp(*), cose(*)
1618    
1619      jmr = jnp - 1
1620      DO j = 2, jnp
1621        ph5 = -0.5*pi + (float(j-1)-0.5)*dp
1622        cose(j) = cos(ph5)
1623      END DO
1624    
1625      jeq = (jnp+1)/2
1626      IF (jmr==2*(jmr/2)) THEN
1627        DO j = jnp, jeq + 1, -1
1628          cose(j) = cose(jnp+2-j)
1629        END DO
1630      ELSE
1631        ! cell edge at equator.
1632        cose(jeq+1) = 1.
1633        DO j = jnp, jeq + 2, -1
1634          cose(j) = cose(jnp+2-j)
1635        END DO
1636      END IF
1637    
1638      DO j = 2, jmr
1639        cosp(j) = 0.5*(cose(j)+cose(j+1))
1640      END DO
1641      cosp(1) = 0.
1642      cosp(jnp) = 0.
1643      RETURN
1644    END SUBROUTINE cosa
1645    
1646    SUBROUTINE cosc(cosp, cose, jnp, pi, dp)
1647      DIMENSION cosp(*), cose(*)
1648    
1649      phi = -0.5*pi
1650      DO j = 2, jnp - 1
1651        phi = phi + dp
1652        cosp(j) = cos(phi)
1653      END DO
1654      cosp(1) = 0.
1655      cosp(jnp) = 0.
1656    
1657      DO j = 2, jnp
1658        cose(j) = 0.5*(cosp(j)+cosp(j-1))
1659      END DO
1660    
1661      DO j = 2, jnp - 1
1662        cosp(j) = 0.5*(cose(j)+cose(j+1))
1663      END DO
1664      RETURN
1665    END SUBROUTINE cosc
1666    
1667    SUBROUTINE qckxyz(q, qtmp, imr, jnp, nlay, j1, j2, cosp, acosp, cross, ic, &
1668        nstep)
1669    
1670      PARAMETER (tiny=1.E-60)
1671      DIMENSION q(imr, jnp, nlay), qtmp(imr, jnp), cosp(*), acosp(*)
1672      LOGICAL cross
1673    
1674      nlaym1 = nlay - 1
1675      len = imr*(j2-j1+1)
1676      ip = 0
1677    
1678      ! Top layer
1679      l = 1
1680      icr = 1
1681      CALL filns(q(1,1,l), imr, jnp, j1, j2, cosp, acosp, ipy, tiny)
1682      IF (ipy==0) GO TO 50
1683      CALL filew(q(1,1,l), qtmp, imr, jnp, j1, j2, ipx, tiny)
1684      IF (ipx==0) GO TO 50
1685    
1686      IF (cross) THEN
1687        CALL filcr(q(1,1,l), imr, jnp, j1, j2, cosp, acosp, icr, tiny)
1688      END IF
1689      IF (icr==0) GO TO 50
1690    
1691      ! Vertical filling...
1692      DO i = 1, len
1693        IF (q(i,j1,1)<0.) THEN
1694        ip = ip + 1        ip = ip + 1
1695  C From above        q(i, j1, 2) = q(i, j1, 2) + q(i, j1, 1)
1696            qup =  Q(I,j1,L-1)        q(i, j1, 1) = 0.
1697            qly = -Q(I,j1,L)      END IF
1698            dup  = min(qly,qup)    END DO
1699            Q(I,j1,L-1) = qup - dup  
1700            Q(I,j1,L  ) = dup-qly  50 CONTINUE
1701  C Below    DO l = 2, nlaym1
1702            Q(I,j1,L+1) = Q(I,j1,L+1) + Q(I,j1,L)      icr = 1
1703            Q(I,j1,L)   = 0.  
1704        ENDIF      CALL filns(q(1,1,l), imr, jnp, j1, j2, cosp, acosp, ipy, tiny)
1705        ENDDO      IF (ipy==0) GO TO 225
1706  225   CONTINUE      CALL filew(q(1,1,l), qtmp, imr, jnp, j1, j2, ipx, tiny)
1707  C      IF (ipx==0) GO TO 225
1708  C BOTTOM LAYER      IF (cross) THEN
1709        sum = 0.        CALL filcr(q(1,1,l), imr, jnp, j1, j2, cosp, acosp, icr, tiny)
1710        L = NLAY      END IF
1711  C      IF (icr==0) GO TO 225
1712        call filns(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,ipy,tiny)  
1713        if(ipy.eq.0) goto 911      DO i = 1, len
1714        call filew(q(1,1,L),qtmp,IMR,JNP,j1,j2,ipx,tiny)        IF (q(i,j1,l)<0.) THEN
1715        if(ipx.eq.0) goto 911  
1716  C          ip = ip + 1
1717        call filcr(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,icr,tiny)          ! From above
1718        if(icr.eq.0) goto 911          qup = q(i, j1, l-1)
1719  C          qly = -q(i, j1, l)
1720        DO  I=1,len          dup = min(qly, qup)
1721        IF( Q(I,j1,L).LT.0.) THEN          q(i, j1, l-1) = qup - dup
1722            q(i, j1, l) = dup - qly
1723            ! Below
1724            q(i, j1, l+1) = q(i, j1, l+1) + q(i, j1, l)
1725            q(i, j1, l) = 0.
1726          END IF
1727        END DO
1728    225 END DO
1729    
1730      ! BOTTOM LAYER
1731      sum = 0.
1732      l = nlay
1733    
1734      CALL filns(q(1,1,l), imr, jnp, j1, j2, cosp, acosp, ipy, tiny)
1735      IF (ipy==0) GO TO 911
1736      CALL filew(q(1,1,l), qtmp, imr, jnp, j1, j2, ipx, tiny)
1737      IF (ipx==0) GO TO 911
1738    
1739      CALL filcr(q(1,1,l), imr, jnp, j1, j2, cosp, acosp, icr, tiny)
1740      IF (icr==0) GO TO 911
1741    
1742      DO i = 1, len
1743        IF (q(i,j1,l)<0.) THEN
1744        ip = ip + 1        ip = ip + 1
1745  c  
1746  C From above        ! From above
1747  C  
1748            qup = Q(I,j1,NLAYM1)        qup = q(i, j1, nlaym1)
1749            qly = -Q(I,j1,L)        qly = -q(i, j1, l)
1750            dup = min(qly,qup)        dup = min(qly, qup)
1751            Q(I,j1,NLAYM1) = qup - dup        q(i, j1, nlaym1) = qup - dup
1752  C From "below" the surface.        ! From "below" the surface.
1753            sum = sum + qly-dup        sum = sum + qly - dup
1754            Q(I,j1,L) = 0.        q(i, j1, l) = 0.
1755         ENDIF      END IF
1756        ENDDO    END DO
1757  C  
1758  911   continue  911 CONTINUE
1759  C  
1760        if(ip.gt.IMR) then    IF (ip>imr) THEN
1761        write(6,*) 'IC=',IC,' STEP=',NSTEP,      WRITE (6, *) 'IC=', ic, ' STEP=', nstep, ' Vertical filling pts=', ip
1762       &           ' Vertical filling pts=',ip    END IF
1763        endif  
1764  C    IF (sum>1.E-25) THEN
1765        if(sum.gt.1.e-25) then      WRITE (6, *) ic, nstep, ' Mass source from the ground=', sum
1766        write(6,*) IC,NSTEP,' Mass source from the ground=',sum    END IF
1767        endif    RETURN
1768        RETURN  END SUBROUTINE qckxyz
1769        END  
1770  C  SUBROUTINE filcr(q, imr, jnp, j1, j2, cosp, acosp, icr, tiny)
1771        subroutine filcr(q,IMR,JNP,j1,j2,cosp,acosp,icr,tiny)    DIMENSION q(imr, *), cosp(*), acosp(*)
1772        dimension q(IMR,*),cosp(*),acosp(*)  
1773        icr = 0    icr = 0
1774        do 65 j=j1+1,j2-1    DO j = j1 + 1, j2 - 1
1775        DO 50 i=1,IMR-1      DO i = 1, imr - 1
1776        IF(q(i,j).LT.0.) THEN        IF (q(i,j)<0.) THEN
1777        icr =  1          icr = 1
1778        dq  = - q(i,j)*cosp(j)          dq = -q(i, j)*cosp(j)
1779  C N-E          ! N-E
1780        dn = q(i+1,j+1)*cosp(j+1)          dn = q(i+1, j+1)*cosp(j+1)
1781        d0 = max(0.,dn)          d0 = max(0., dn)
1782        d1 = min(dq,d0)          d1 = min(dq, d0)
1783        q(i+1,j+1) = (dn - d1)*acosp(j+1)          q(i+1, j+1) = (dn-d1)*acosp(j+1)
1784        dq = dq - d1          dq = dq - d1
1785  C S-E          ! S-E
1786        ds = q(i+1,j-1)*cosp(j-1)          ds = q(i+1, j-1)*cosp(j-1)
1787        d0 = max(0.,ds)          d0 = max(0., ds)
1788        d2 = min(dq,d0)          d2 = min(dq, d0)
1789        q(i+1,j-1) = (ds - d2)*acosp(j-1)          q(i+1, j-1) = (ds-d2)*acosp(j-1)
1790        q(i,j) = (d2 - dq)*acosp(j) + tiny          q(i, j) = (d2-dq)*acosp(j) + tiny
1791        endif        END IF
1792  50    continue      END DO
1793        if(icr.eq.0 .and. q(IMR,j).ge.0.) goto 65      IF (icr==0 .AND. q(imr,j)>=0.) GO TO 65
1794        DO 55 i=2,IMR      DO i = 2, imr
1795        IF(q(i,j).LT.0.) THEN        IF (q(i,j)<0.) THEN
1796        icr =  1          icr = 1
1797        dq  = - q(i,j)*cosp(j)          dq = -q(i, j)*cosp(j)
1798  C N-W          ! N-W
1799        dn = q(i-1,j+1)*cosp(j+1)          dn = q(i-1, j+1)*cosp(j+1)
1800        d0 = max(0.,dn)          d0 = max(0., dn)
1801        d1 = min(dq,d0)          d1 = min(dq, d0)
1802        q(i-1,j+1) = (dn - d1)*acosp(j+1)          q(i-1, j+1) = (dn-d1)*acosp(j+1)
1803        dq = dq - d1          dq = dq - d1
1804  C S-W          ! S-W
1805        ds = q(i-1,j-1)*cosp(j-1)          ds = q(i-1, j-1)*cosp(j-1)
1806        d0 = max(0.,ds)          d0 = max(0., ds)
1807        d2 = min(dq,d0)          d2 = min(dq, d0)
1808        q(i-1,j-1) = (ds - d2)*acosp(j-1)          q(i-1, j-1) = (ds-d2)*acosp(j-1)
1809        q(i,j) = (d2 - dq)*acosp(j) + tiny          q(i, j) = (d2-dq)*acosp(j) + tiny
1810        endif        END IF
1811  55    continue      END DO
1812  C *****************************************      ! *****************************************
1813  C i=1      ! i=1
1814        i=1      i = 1
1815        IF(q(i,j).LT.0.) THEN      IF (q(i,j)<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  
1816        icr = 1        icr = 1
1817        goto 80        dq = -q(i, j)*cosp(j)
1818        endif        ! N-W
1819        enddo        dn = q(imr, j+1)*cosp(j+1)
1820  C        d0 = max(0., dn)
1821  80    continue        d1 = min(dq, d0)
1822  C        q(imr, j+1) = (dn-d1)*acosp(j+1)
1823        if(q(1,1).lt.0. .or. q(1,jnp).lt.0.) then        dq = dq - d1
1824          ! S-W
1825          ds = q(imr, j-1)*cosp(j-1)
1826          d0 = max(0., ds)
1827          d2 = min(dq, d0)
1828          q(imr, j-1) = (ds-d2)*acosp(j-1)
1829          q(i, j) = (d2-dq)*acosp(j) + tiny
1830        END IF
1831        ! *****************************************
1832        ! i=IMR
1833        i = imr
1834        IF (q(i,j)<0.) THEN
1835        icr = 1        icr = 1
1836        endif        dq = -q(i, j)*cosp(j)
1837  C        ! N-E
1838        return        dn = q(1, j+1)*cosp(j+1)
1839        end        d0 = max(0., dn)
1840  C        d1 = min(dq, d0)
1841        subroutine filns(q,IMR,JNP,j1,j2,cosp,acosp,ipy,tiny)        q(1, j+1) = (dn-d1)*acosp(j+1)
       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)  
1842        dq = dq - d1        dq = dq - d1
1843  C South        ! S-E
1844        ds = q(i,j-1)*cosp(j-1)        ds = q(1, j-1)*cosp(j-1)
1845        d0 = max(0.,ds)        d0 = max(0., ds)
1846        d2 = min(dq,d0)        d2 = min(dq, d0)
1847        q(i,j-1) = (ds - d2)*acosp(j-1)        q(1, j-1) = (ds-d2)*acosp(j-1)
1848        q(i,j) = (d2 - dq)*acosp(j) + tiny        q(i, j) = (d2-dq)*acosp(j) + tiny
1849        endif      END IF
1850  55    continue      ! *****************************************
1851  C  65 END DO
1852        do i=1,imr  
1853        IF(q(i,j1).LT.0.) THEN    DO i = 1, imr
1854        ipy =  1      IF (q(i,j1)<0. .OR. q(i,j2)<0.) THEN
1855        dq  = - q(i,j1)*cosp(j1)        icr = 1
1856  C North        GO TO 80
1857        dn = q(i,j1+1)*cosp(j1+1)      END IF
1858        d0 = max(0.,dn)    END DO
1859        d1 = min(dq,d0)  
1860        q(i,j1+1) = (dn - d1)*acosp(j1+1)  80 CONTINUE
1861        q(i,j1) = (d1 - dq)*acosp(j1) + tiny  
1862        endif    IF (q(1,1)<0. .OR. q(1,jnp)<0.) THEN
1863        enddo      icr = 1
1864  C    END IF
1865        j = j2  
1866        do i=1,imr    RETURN
1867        IF(q(i,j).LT.0.) THEN  END SUBROUTINE filcr
1868        ipy =  1  
1869        dq  = - q(i,j)*cosp(j)  SUBROUTINE filns(q, imr, jnp, j1, j2, cosp, acosp, ipy, tiny)
1870  C South    DIMENSION q(imr, *), cosp(*), acosp(*)
1871        ds = q(i,j-1)*cosp(j-1)    ! logical first
1872        d0 = max(0.,ds)    ! data first /.true./
1873        d2 = min(dq,d0)    ! save cap1
1874        q(i,j-1) = (ds - d2)*acosp(j-1)  
1875        q(i,j) = (d2 - dq)*acosp(j) + tiny    ! if(first) then
1876        endif    dp = 4.*atan(1.)/float(jnp-1)
1877        enddo    cap1 = imr*(1.-cos((j1-1.5)*dp))/dp
1878  C    ! first = .false.
1879  C Check Poles.    ! endif
1880        if(q(1,1).lt.0.) then  
1881        dq = q(1,1)*cap1/float(IMR)*acosp(j1)    ipy = 0
1882        do i=1,imr    DO j = j1 + 1, j2 - 1
1883        q(i,1) = 0.      DO i = 1, imr
1884        q(i,j1) = q(i,j1) + dq        IF (q(i,j)<0.) THEN
1885        if(q(i,j1).lt.0.) ipy = 1          ipy = 1
1886        enddo          dq = -q(i, j)*cosp(j)
1887        endif          ! North
1888  C          dn = q(i, j+1)*cosp(j+1)
1889        if(q(1,JNP).lt.0.) then          d0 = max(0., dn)
1890        dq = q(1,JNP)*cap1/float(IMR)*acosp(j2)          d1 = min(dq, d0)
1891        do i=1,imr          q(i, j+1) = (dn-d1)*acosp(j+1)
1892        q(i,JNP) = 0.          dq = dq - d1
1893        q(i,j2) = q(i,j2) + dq          ! South
1894        if(q(i,j2).lt.0.) ipy = 1          ds = q(i, j-1)*cosp(j-1)
1895        enddo          d0 = max(0., ds)
1896        endif          d2 = min(dq, d0)
1897  C          q(i, j-1) = (ds-d2)*acosp(j-1)
1898        return          q(i, j) = (d2-dq)*acosp(j) + tiny
1899        end        END IF
1900  C      END DO
1901        subroutine filew(q,qtmp,IMR,JNP,j1,j2,ipx,tiny)    END DO
1902        dimension q(IMR,*),qtmp(JNP,IMR)  
1903  C    DO i = 1, imr
1904        ipx = 0      IF (q(i,j1)<0.) THEN
1905  C Copy & swap direction for vectorization.        ipy = 1
1906        do 25 i=1,imr        dq = -q(i, j1)*cosp(j1)
1907        do 25 j=j1,j2        ! North
1908  25    qtmp(j,i) = q(i,j)        dn = q(i, j1+1)*cosp(j1+1)
1909  C        d0 = max(0., dn)
1910        do 55 i=2,imr-1        d1 = min(dq, d0)
1911        do 55 j=j1,j2        q(i, j1+1) = (dn-d1)*acosp(j1+1)
1912        if(qtmp(j,i).lt.0.) then        q(i, j1) = (d1-dq)*acosp(j1) + tiny
1913        ipx =  1      END IF
1914  c west    END DO
1915        d0 = max(0.,qtmp(j,i-1))  
1916        d1 = min(-qtmp(j,i),d0)    j = j2
1917        qtmp(j,i-1) = qtmp(j,i-1) - d1    DO i = 1, imr
1918        qtmp(j,i) = qtmp(j,i) + d1      IF (q(i,j)<0.) THEN
1919  c east        ipy = 1
1920        d0 = max(0.,qtmp(j,i+1))        dq = -q(i, j)*cosp(j)
1921        d2 = min(-qtmp(j,i),d0)        ! South
1922        qtmp(j,i+1) = qtmp(j,i+1) - d2        ds = q(i, j-1)*cosp(j-1)
1923        qtmp(j,i) = qtmp(j,i) + d2 + tiny        d0 = max(0., ds)
1924        endif        d2 = min(dq, d0)
1925  55    continue        q(i, j-1) = (ds-d2)*acosp(j-1)
1926  c        q(i, j) = (d2-dq)*acosp(j) + tiny
1927        i=1      END IF
1928        do 65 j=j1,j2    END DO
1929        if(qtmp(j,i).lt.0.) then  
1930        ipx =  1    ! Check Poles.
1931  c west    IF (q(1,1)<0.) THEN
1932        d0 = max(0.,qtmp(j,imr))      dq = q(1, 1)*cap1/float(imr)*acosp(j1)
1933        d1 = min(-qtmp(j,i),d0)      DO i = 1, imr
1934        qtmp(j,imr) = qtmp(j,imr) - d1        q(i, 1) = 0.
1935        qtmp(j,i) = qtmp(j,i) + d1        q(i, j1) = q(i, j1) + dq
1936  c east        IF (q(i,j1)<0.) ipy = 1
1937        d0 = max(0.,qtmp(j,i+1))      END DO
1938        d2 = min(-qtmp(j,i),d0)    END IF
1939        qtmp(j,i+1) = qtmp(j,i+1) - d2  
1940  c    IF (q(1,jnp)<0.) THEN
1941        qtmp(j,i) = qtmp(j,i) + d2 + tiny      dq = q(1, jnp)*cap1/float(imr)*acosp(j2)
1942        endif      DO i = 1, imr
1943  65    continue        q(i, jnp) = 0.
1944        i=IMR        q(i, j2) = q(i, j2) + dq
1945        do 75 j=j1,j2        IF (q(i,j2)<0.) ipy = 1
1946        if(qtmp(j,i).lt.0.) then      END DO
1947        ipx =  1    END IF
1948  c west  
1949        d0 = max(0.,qtmp(j,i-1))    RETURN
1950        d1 = min(-qtmp(j,i),d0)  END SUBROUTINE filns
1951        qtmp(j,i-1) = qtmp(j,i-1) - d1  
1952        qtmp(j,i) = qtmp(j,i) + d1  SUBROUTINE filew(q, qtmp, imr, jnp, j1, j2, ipx, tiny)
1953  c east    DIMENSION q(imr, *), qtmp(jnp, imr)
1954        d0 = max(0.,qtmp(j,1))  
1955        d2 = min(-qtmp(j,i),d0)    ipx = 0
1956        qtmp(j,1) = qtmp(j,1) - d2    ! Copy & swap direction for vectorization.
1957  c    DO i = 1, imr
1958        qtmp(j,i) = qtmp(j,i) + d2 + tiny      DO j = j1, j2
1959        endif        qtmp(j, i) = q(i, j)
1960  75    continue      END DO
1961  C    END DO
1962        if(ipx.ne.0) then  
1963        do 85 j=j1,j2    DO i = 2, imr - 1
1964        do 85 i=1,imr      DO j = j1, j2
1965  85    q(i,j) = qtmp(j,i)        IF (qtmp(j,i)<0.) THEN
1966        else          ipx = 1
1967  C          ! west
1968  C Poles.          d0 = max(0., qtmp(j,i-1))
1969        if(q(1,1).lt.0. or. q(1,JNP).lt.0.) ipx = 1          d1 = min(-qtmp(j,i), d0)
1970        endif          qtmp(j, i-1) = qtmp(j, i-1) - d1
1971        return          qtmp(j, i) = qtmp(j, i) + d1
1972        end          ! east
1973  C          d0 = max(0., qtmp(j,i+1))
1974        subroutine zflip(q,im,km,nc)          d2 = min(-qtmp(j,i), d0)
1975  C This routine flip the array q (in the vertical).          qtmp(j, i+1) = qtmp(j, i+1) - d2
1976        real q(im,km,nc)          qtmp(j, i) = qtmp(j, i) + d2 + tiny
1977  C local dynamic array        END IF
1978        real qtmp(im,km)      END DO
1979  C    END DO
1980        do 4000 IC = 1, nc  
1981  C    i = 1
1982        do 1000 k=1,km    DO j = j1, j2
1983        do 1000 i=1,im      IF (qtmp(j,i)<0.) THEN
1984        qtmp(i,k) = q(i,km+1-k,IC)        ipx = 1
1985  1000  continue        ! west
1986  C        d0 = max(0., qtmp(j,imr))
1987        do 2000 i=1,im*km        d1 = min(-qtmp(j,i), d0)
1988  2000  q(i,1,IC) = qtmp(i,1)        qtmp(j, imr) = qtmp(j, imr) - d1
1989  4000  continue        qtmp(j, i) = qtmp(j, i) + d1
1990        return        ! east
1991        end        d0 = max(0., qtmp(j,i+1))
1992          d2 = min(-qtmp(j,i), d0)
1993          qtmp(j, i+1) = qtmp(j, i+1) - d2
1994    
1995          qtmp(j, i) = qtmp(j, i) + d2 + tiny
1996        END IF
1997      END DO
1998      i = imr
1999      DO j = j1, j2
2000        IF (qtmp(j,i)<0.) THEN
2001          ipx = 1
2002          ! west
2003          d0 = max(0., qtmp(j,i-1))
2004          d1 = min(-qtmp(j,i), d0)
2005          qtmp(j, i-1) = qtmp(j, i-1) - d1
2006          qtmp(j, i) = qtmp(j, i) + d1
2007          ! east
2008          d0 = max(0., qtmp(j,1))
2009          d2 = min(-qtmp(j,i), d0)
2010          qtmp(j, 1) = qtmp(j, 1) - d2
2011    
2012          qtmp(j, i) = qtmp(j, i) + d2 + tiny
2013        END IF
2014      END DO
2015    
2016      IF (ipx/=0) THEN
2017        DO j = j1, j2
2018          DO i = 1, imr
2019            q(i, j) = qtmp(j, i)
2020          END DO
2021        END DO
2022      ELSE
2023    
2024        ! Poles.
2025        IF (q(1,1)<0 .OR. q(1,jnp)<0.) ipx = 1
2026      END IF
2027      RETURN
2028    END SUBROUTINE filew
Legend:
Removed from v.3  
changed lines
  Added in v.81
  ViewVC Help
Powered by ViewVC 1.1.21