/[lmdze]/trunk/dyn3d/ppm3d.f90
ViewVC logotype

Diff of /trunk/dyn3d/ppm3d.f90

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

trunk/dyn3d/ppm3d.f revision 80 by guez, Fri Nov 15 18:45:49 2013 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 number of constituents    ! Q(IMR,JNP,NLAY,NC): mixing ratios at current time (t)
84  C IMR: first dimension (E-W); number of Grid intervals in E-W is IMR    ! NC: total number of constituents
85  C JNP: 2nd dimension (N-S); number 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 (number 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 number 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        DTDX5(j) = 0.5*DTDX(j)        jn0 = jnp - js0 + 1
401        enddo      END IF
402  C  
403          
404        DTDY  = DT /(AE*DP)      DO j = 2, jmr
405        DTDY5 = 0.5*DTDY        dtdx(j) = dt/(dl*ae*cosp(j))
406  C  
407        endif        dtdx5(j) = 0.5*dtdx(j)
408  C      END DO
409  C *********** End Initialization **********************  
410  C  
411  C delp = pressure thickness: the psudo-density in a hydrostatic system.      dtdy = dt/(ae*dp)
412        do  k=1,NLAY      dtdy5 = 0.5*dtdy
413           do  j=1,JNP  
414              do  i=1,IMR    END IF
415                 delp1(i,j,k)=DAP(k)+DBK(k)*PS1(i,j)  
416                 delp2(i,j,k)=DAP(k)+DBK(k)*PS2(i,j)          ! *********** End Initialization **********************
417              enddo  
418           enddo    ! delp = pressure thickness: the psudo-density in a hydrostatic system.
419        enddo    DO k = 1, nlay
420                  DO j = 1, jnp
421  C        DO i = 1, imr
422        if(j1.ne.2) then          delp1(i, j, k) = dap(k) + dbk(k)*ps1(i, j)
423        DO 40 IC=1,NC          delp2(i, j, k) = dap(k) + dbk(k)*ps2(i, j)
424        DO 40 L=1,NLAY        END DO
425        DO 40 I=1,IMR      END DO
426        Q(I,  2,L,IC) = Q(I,  1,L,IC)    END DO
427  40    Q(I,JMR,L,IC) = Q(I,JNP,L,IC)  
428        endif  
429  C    IF (j1/=2) THEN
430  C Compute "tracer density"      DO ic = 1, nc
431        DO 550 IC=1,NC        DO l = 1, nlay
432        DO 44 k=1,NLAY          DO i = 1, imr
433        DO 44 j=1,JNP            q(i, 2, l, ic) = q(i, 1, l, ic)
434        DO 44 i=1,IMR            q(i, jmr, l, ic) = q(i, jnp, l, ic)
435  44    DQ(i,j,k,IC) = Q(i,j,k,IC)*delp1(i,j,k)          END DO
436  550     continue        END DO
437  C      END DO
438        do 1500 k=1,NLAY    END IF
439  C  
440        if(IGD.eq.0) then    ! Compute "tracer density"
441  C Convert winds on A-Grid to Courant number on C-Grid.    DO ic = 1, nc
442        call A2C(U(1,1,k),V(1,1,k),IMR,JMR,j1,j2,CRX,CRY,dtdx5,DTDY5)      DO k = 1, nlay
443        else        DO j = 1, jnp
444  C Convert winds on C-grid to Courant number          DO i = 1, imr
445        do 45 j=j1,j2            dq(i, j, k, ic) = q(i, j, k, ic)*delp1(i, j, k)
446        do 45 i=2,IMR          END DO
447  45    CRX(i,J) = dtdx(j)*U(i-1,j,k)        END DO
448          END DO
449  C    END DO
450        do 50 j=j1,j2  
451  50    CRX(1,J) = dtdx(j)*U(IMR,j,k)    DO k = 1, nlay
452  C  
453        do 55 i=1,IMR*JMR      IF (igd==0) THEN
454  55    CRY(i,2) = DTDY*V(i,1,k)        ! Convert winds on A-Grid to Courant number on C-Grid.
455        endif        CALL a2c(u(1,1,k), v(1,1,k), imr, jmr, j1, j2, crx, cry, dtdx5, dtdy5)
456  C          ELSE
457  C Determine JS and JN        ! Convert winds on C-grid to Courant number
458        JS = j1        DO j = j1, j2
459        JN = j2          DO i = 2, imr
460  C            crx(i, j) = dtdx(j)*u(i-1, j, k)
461        do j=JS0,j1+1,-1          END DO
462        do i=1,IMR        END DO
463        if(abs(CRX(i,j)).GT.1.) then  
464              JS = j  
465              go to 2222        DO j = j1, j2
466        endif          crx(1, j) = dtdx(j)*u(imr, j, k)
467        enddo        END DO
468        enddo  
469  C        DO i = 1, imr*jmr
470  2222  continue          cry(i, 2) = dtdy*v(i, 1, k)
471        do j=JN0,j2-1        END DO
472        do i=1,IMR      END IF
473        if(abs(CRX(i,j)).GT.1.) then  
474              JN = j      ! Determine JS and JN
475              go to 2233      js = j1
476        endif      jn = j2
477        enddo  
478        enddo      DO j = js0, j1 + 1, -1
479  2233  continue        DO i = 1, imr
480  C          IF (abs(crx(i,j))>1.) THEN
481        if(j1.ne.2) then           ! Enlarged polar cap.            js = j
482        do i=1,IMR            GO TO 2222
483        DPI(i,  2,k) = 0.          END IF
484        DPI(i,JMR,k) = 0.        END DO
485        enddo      END DO
486        endif  
487  C  2222 CONTINUE
488  C ******* Compute horizontal mass fluxes ************      DO j = jn0, j2 - 1
489  C        DO i = 1, imr
490  C N-S component          IF (abs(crx(i,j))>1.) THEN
491        do j=j1,j2+1            jn = j
492        D5 = 0.5 * COSE(j)            GO TO 2233
493        do i=1,IMR          END IF
494        ymass(i,j) = CRY(i,j)*D5*(delp2(i,j,k) + delp2(i,j-1,k))        END DO
495        enddo      END DO
496        enddo  2233 CONTINUE
497  C  
498        do 95 j=j1,j2      IF (j1/=2) THEN ! Enlarged polar cap.
499        DO 95 i=1,IMR        DO i = 1, imr
500  95    DPI(i,j,k) = (ymass(i,j) - ymass(i,j+1)) * acosp(j)          dpi(i, 2, k) = 0.
501  C          dpi(i, jmr, k) = 0.
502  C Poles        END DO
503        sum1 = ymass(IMR,j1  )      END IF
504        sum2 = ymass(IMR,J2+1)  
505        do i=1,IMR-1      ! ******* Compute horizontal mass fluxes ************
506        sum1 = sum1 + ymass(i,j1  )  
507        sum2 = sum2 + ymass(i,J2+1)      ! N-S component
508        enddo      DO j = j1, j2 + 1
509  C        d5 = 0.5*cose(j)
510        sum1 = - sum1 * RCAP        DO i = 1, imr
511        sum2 =   sum2 * RCAP          ymass(i, j) = cry(i, j)*d5*(delp2(i,j,k)+delp2(i,j-1,k))
512        do i=1,IMR        END DO
513        DPI(i,  1,k) = sum1      END DO
514        DPI(i,JNP,k) = sum2  
515        enddo      DO j = j1, j2
516  C        DO i = 1, imr
517  C E-W component          dpi(i, j, k) = (ymass(i,j)-ymass(i,j+1))*acosp(j)
518  C        END DO
519        do j=j1,j2      END DO
520        do i=2,IMR  
521        PU(i,j) = 0.5 * (delp2(i,j,k) + delp2(i-1,j,k))      ! Poles
522        enddo      sum1 = ymass(imr, j1)
523        enddo      sum2 = ymass(imr, j2+1)
524  C      DO i = 1, imr - 1
525        do j=j1,j2        sum1 = sum1 + ymass(i, j1)
526        PU(1,j) = 0.5 * (delp2(1,j,k) + delp2(IMR,j,k))        sum2 = sum2 + ymass(i, j2+1)
527        enddo      END DO
528  C  
529        do 110 j=j1,j2      sum1 = -sum1*rcap
530        DO 110 i=1,IMR      sum2 = sum2*rcap
531  110   xmass(i,j) = PU(i,j)*CRX(i,j)      DO i = 1, imr
532  C        dpi(i, 1, k) = sum1
533        DO 120 j=j1,j2        dpi(i, jnp, k) = sum2
534        DO 120 i=1,IMR-1      END DO
535  120   DPI(i,j,k) = DPI(i,j,k) + xmass(i,j) - xmass(i+1,j)  
536  C      ! E-W component
537        DO 130 j=j1,j2  
538  130   DPI(IMR,j,k) = DPI(IMR,j,k) + xmass(IMR,j) - xmass(1,j)      DO j = j1, j2
539  C        DO i = 2, imr
540        DO j=j1,j2          pu(i, j) = 0.5*(delp2(i,j,k)+delp2(i-1,j,k))
541        do i=1,IMR-1        END DO
542        UA(i,j) = 0.5 * (CRX(i,j)+CRX(i+1,j))      END DO
543        enddo  
544        enddo      DO j = j1, j2
545  C        pu(1, j) = 0.5*(delp2(1,j,k)+delp2(imr,j,k))
546        DO j=j1,j2      END DO
547        UA(imr,j) = 0.5 * (CRX(imr,j)+CRX(1,j))  
548        enddo      DO j = j1, j2
549  ccccccccccccccccccccccccccccccccccccccccccccccccccccccc        DO i = 1, imr
550  c Rajouts pour LMDZ.3.3          xmass(i, j) = pu(i, j)*crx(i, j)
551  ccccccccccccccccccccccccccccccccccccccccccccccccccccccc        END DO
552        do i=1,IMR      END DO
553           do j=1,JNP  
554               VA(i,j)=0.      DO j = j1, j2
555           enddo        DO i = 1, imr - 1
556        enddo          dpi(i, j, k) = dpi(i, j, k) + xmass(i, j) - xmass(i+1, j)
557          END DO
558        do i=1,imr*(JMR-1)      END DO
559        VA(i,2) = 0.5*(CRY(i,2)+CRY(i,3))  
560        enddo      DO j = j1, j2
561  C        dpi(imr, j, k) = dpi(imr, j, k) + xmass(imr, j) - xmass(1, j)
562        if(j1.eq.2) then      END DO
563          IMH = IMR/2  
564        do i=1,IMH      DO j = j1, j2
565        VA(i,      1) = 0.5*(CRY(i,2)-CRY(i+IMH,2))        DO i = 1, imr - 1
566        VA(i+IMH,  1) = -VA(i,1)          ua(i, j) = 0.5*(crx(i,j)+crx(i+1,j))
567        VA(i,    JNP) = 0.5*(CRY(i,JNP)-CRY(i+IMH,JMR))        END DO
568        VA(i+IMH,JNP) = -VA(i,JNP)      END DO
569        enddo  
570        VA(IMR,1)=VA(1,1)      DO j = j1, j2
571        VA(IMR,JNP)=VA(1,JNP)        ua(imr, j) = 0.5*(crx(imr,j)+crx(1,j))
572        endif      END DO
573  C      ! cccccccccccccccccccccccccccccccccccccccccccccccccccccc
574  C ****6***0*********0*********0*********0*********0*********0**********72      ! Rajouts pour LMDZ.3.3
575        do 1000 IC=1,NC      ! cccccccccccccccccccccccccccccccccccccccccccccccccccccc
576  C      DO i = 1, imr
577        do i=1,IMJM        DO j = 1, jnp
578        wk1(i,1,1) = 0.          va(i, j) = 0.
579        wk1(i,1,2) = 0.        END DO
580        enddo      END DO
581  C  
582  C E-W advective cross term      DO i = 1, imr*(jmr-1)
583        do 250 j=J1,J2        va(i, 2) = 0.5*(cry(i,2)+cry(i,3))
584        if(J.GT.JS  .and. J.LT.JN) GO TO 250      END DO
585  C  
586        do i=1,IMR      IF (j1==2) THEN
587        qtmp(i) = q(i,j,k,IC)        imh = imr/2
588        enddo        DO i = 1, imh
589  C          va(i, 1) = 0.5*(cry(i,2)-cry(i+imh,2))
590        do i=-IML,0          va(i+imh, 1) = -va(i, 1)
591        qtmp(i)       = q(IMR+i,j,k,IC)          va(i, jnp) = 0.5*(cry(i,jnp)-cry(i+imh,jmr))
592        qtmp(IMR+1-i) = q(1-i,j,k,IC)          va(i+imh, jnp) = -va(i, jnp)
593        enddo        END DO
594  C        va(imr, 1) = va(1, 1)
595        DO 230 i=1,IMR        va(imr, jnp) = va(1, jnp)
596        iu = UA(i,j)      END IF
597        ru = UA(i,j) - iu  
598        iiu = i-iu      ! ****6***0*********0*********0*********0*********0*********0**********72
599        if(UA(i,j).GE.0.) then      DO ic = 1, nc
600        wk1(i,j,1) = qtmp(iiu)+ru*(qtmp(iiu-1)-qtmp(iiu))  
601        else        DO i = 1, imjm
602        wk1(i,j,1) = qtmp(iiu)+ru*(qtmp(iiu)-qtmp(iiu+1))          wk1(i, 1, 1) = 0.
603        endif          wk1(i, 1, 2) = 0.
604        wk1(i,j,1) = wk1(i,j,1) - qtmp(i)        END DO
605  230   continue  
606  250   continue        ! E-W advective cross term
607  C        DO j = j1, j2
608        if(JN.ne.0) then          IF (j>js .AND. j<jn) GO TO 250
609        do j=JS+1,JN-1  
610  C          DO i = 1, imr
611        do i=1,IMR            qtmp(i) = q(i, j, k, ic)
612        qtmp(i) = q(i,j,k,IC)          END DO
613        enddo  
614  C          DO i = -iml, 0
615        qtmp(0)     = q(IMR,J,k,IC)            qtmp(i) = q(imr+i, j, k, ic)
616        qtmp(IMR+1) = q(  1,J,k,IC)            qtmp(imr+1-i) = q(1-i, j, k, ic)
617  C          END DO
618        do i=1,imr  
619        iu = i - UA(i,j)          DO i = 1, imr
620        wk1(i,j,1) = UA(i,j)*(qtmp(iu) - qtmp(iu+1))            iu = ua(i, j)
621        enddo            ru = ua(i, j) - iu
622        enddo            iiu = i - iu
623        endif            IF (ua(i,j)>=0.) THEN
624  C ****6***0*********0*********0*********0*********0*********0**********72              wk1(i, j, 1) = qtmp(iiu) + ru*(qtmp(iiu-1)-qtmp(iiu))
625  C Contribution from the N-S advection            ELSE
626        do i=1,imr*(j2-j1+1)              wk1(i, j, 1) = qtmp(iiu) + ru*(qtmp(iiu)-qtmp(iiu+1))
627        JT = float(J1) - VA(i,j1)            END IF
628        wk1(i,j1,2) = VA(i,j1) * (q(i,jt,k,IC) - q(i,jt+1,k,IC))            wk1(i, j, 1) = wk1(i, j, 1) - qtmp(i)
629        enddo          END DO
630  C  250   END DO
631        do i=1,IMJM  
632        wk1(i,1,1) = q(i,1,k,IC) + 0.5*wk1(i,1,1)        IF (jn/=0) THEN
633        wk1(i,1,2) = q(i,1,k,IC) + 0.5*wk1(i,1,2)          DO j = js + 1, jn - 1
634        enddo  
635  C            DO i = 1, imr
636          if(cross) then              qtmp(i) = q(i, j, k, ic)
637  C Add cross terms in the vertical direction.            END DO
638          if(IORD .GE. 2) then  
639                  iad = 2            qtmp(0) = q(imr, j, k, ic)
640          else            qtmp(imr+1) = q(1, j, k, ic)
641                  iad = 1  
642          endif            DO i = 1, imr
643  C              iu = i - ua(i, j)
644          if(JORD .GE. 2) then              wk1(i, j, 1) = ua(i, j)*(qtmp(iu)-qtmp(iu+1))
645                  jad = 2            END DO
646          else          END DO
647                  jad = 1        END IF
648          endif        ! ****6***0*********0*********0*********0*********0*********0**********72
649        call xadv(IMR,JNP,j1,j2,wk1(1,1,2),UA,JS,JN,IML,DC2,iad)        ! Contribution from the N-S advection
650        call yadv(IMR,JNP,j1,j2,wk1(1,1,1),VA,PV,W,jad)        DO i = 1, imr*(j2-j1+1)
651        do j=1,JNP          jt = float(j1) - va(i, j1)
652        do i=1,IMR          wk1(i, j1, 2) = va(i, j1)*(q(i,jt,k,ic)-q(i,jt+1,k,ic))
653        q(i,j,k,IC) = q(i,j,k,IC) + DC2(i,j) + PV(i,j)        END DO
654        enddo  
655        enddo        DO i = 1, imjm
656        endif          wk1(i, 1, 1) = q(i, 1, k, ic) + 0.5*wk1(i, 1, 1)
657  C          wk1(i, 1, 2) = q(i, 1, k, ic) + 0.5*wk1(i, 1, 2)
658        call xtp(IMR,JNP,IML,j1,j2,JN,JS,PU,DQ(1,1,k,IC),wk1(1,1,2)        END DO
659       &        ,CRX,fx1,xmass,IORD)  
660          IF (cross) THEN
661        call ytp(IMR,JNP,j1,j2,acosp,RCAP,DQ(1,1,k,IC),wk1(1,1,1),CRY,          ! Add cross terms in the vertical direction.
662       &  DC2,ymass,WK1(1,1,3),wk1(1,1,4),WK1(1,1,5),WK1(1,1,6),JORD)          IF (iord>=2) THEN
663  C            iad = 2
664  1000  continue          ELSE
665  1500  continue            iad = 1
666  C          END IF
667  C ******* Compute vertical mass flux (same unit as PS) ***********  
668  C          IF (jord>=2) THEN
669  C 1st step: compute total column mass CONVERGENCE.            jad = 2
670  C          ELSE
671        do 320 j=1,JNP            jad = 1
672        do 320 i=1,IMR          END IF
673  320   CRY(i,j) = DPI(i,j,1)          CALL xadv(imr, jnp, j1, j2, wk1(1,1,2), ua, js, jn, iml, dc2, iad)
674  C          CALL yadv(imr, jnp, j1, j2, wk1(1,1,1), va, pv, w, jad)
675        do 330 k=2,NLAY          DO j = 1, jnp
676        do 330 j=1,JNP            DO i = 1, imr
677        do 330 i=1,IMR              q(i, j, k, ic) = q(i, j, k, ic) + dc2(i, j) + pv(i, j)
678        CRY(i,j)  = CRY(i,j) + DPI(i,j,k)            END DO
679  330   continue          END DO
680  C        END IF
681        do 360 j=1,JNP  
682        do 360 i=1,IMR        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  C 2nd step: compute PS2 (PS at n+1) using the hydrostatic assumption.  
685  C Changes (increases) to surface pressure = total column mass convergence        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        PS2(i,j)  = PS1(i,j) + CRY(i,j)  
688  C      END DO
689  C 3rd step: compute vertical mass flux from mass conservation principle.    END DO
690  C  
691        W(i,j,1) = DPI(i,j,1) - DBK(1)*CRY(i,j)    ! ******* Compute vertical mass flux (same unit as PS) ***********
692        W(i,j,NLAY) = 0.  
693  360   continue    ! 1st step: compute total column mass CONVERGENCE.
694  C  
695        do 370 k=2,NLAY-1    DO j = 1, jnp
696        do 370 j=1,JNP      DO i = 1, imr
697        do 370 i=1,IMR        cry(i, j) = dpi(i, j, 1)
698        W(i,j,k) = W(i,j,k-1) + DPI(i,j,k) - DBK(k)*CRY(i,j)      END DO
699  370   continue    END DO
700  C  
701        DO 380 k=1,NLAY    DO k = 2, nlay
702        DO 380 j=1,JNP      DO j = 1, jnp
703        DO 380 i=1,IMR        DO i = 1, imr
704        delp2(i,j,k) = DAP(k) + DBK(k)*PS2(i,j)          cry(i, j) = cry(i, j) + dpi(i, j, k)
705  380   continue        END DO
706  C      END DO
707          KRD = max(3, KORD)    END DO
708        do 4000 IC=1,NC  
709  C    DO j = 1, jnp
710  C****6***0*********0*********0*********0*********0*********0**********72      DO i = 1, imr
711      
712        call FZPPM(IMR,JNP,NLAY,j1,DQ(1,1,1,IC),W,Q(1,1,1,IC),WK1,DPI,        ! 2nd step: compute PS2 (PS at n+1) using the hydrostatic assumption.
713       &           DC2,CRX,CRY,PU,PV,xmass,ymass,delp1,KRD)        ! Changes (increases) to surface pressure = total column mass
714  C        ! convergence
715        
716        if(fill) call qckxyz(DQ(1,1,1,IC),DC2,IMR,JNP,NLAY,j1,j2,        ps2(i, j) = ps1(i, j) + cry(i, j)
717       &                     cosp,acosp,.false.,IC,NSTEP)  
718  C        ! 3rd step: compute vertical mass flux from mass conservation
719  C Recover tracer mixing ratio from "density" using predicted        ! principle.
720  C "air density" (pressure thickness) at time-level n+1  
721  C        w(i, j, 1) = dpi(i, j, 1) - dbk(1)*cry(i, j)
722        DO k=1,NLAY        w(i, j, nlay) = 0.
723        DO j=1,JNP      END DO
724        DO i=1,IMR    END DO
725              Q(i,j,k,IC) = DQ(i,j,k,IC) / delp2(i,j,k)  
726        enddo    DO k = 2, nlay - 1
727        enddo      DO j = 1, jnp
728        enddo        DO i = 1, imr
729  C              w(i, j, k) = w(i, j, k-1) + dpi(i, j, k) - dbk(k)*cry(i, j)
730        if(j1.ne.2) then        END DO
731        DO 400 k=1,NLAY      END DO
732        DO 400 I=1,IMR    END DO
733  c     j=1 c'est le pôle Sud, j=JNP c'est le pôle Nord  
734        Q(I,  2,k,IC) = Q(I,  1,k,IC)    DO k = 1, nlay
735        Q(I,JMR,k,IC) = Q(I,JMP,k,IC)      DO j = 1, jnp
736  400   CONTINUE        DO i = 1, imr
737        endif          delp2(i, j, k) = dap(k) + dbk(k)*ps2(i, j)
738  4000  continue        END DO
739  C      END DO
740        if(j1.ne.2) then    END DO
741        DO 5000 k=1,NLAY  
742        DO 5000 i=1,IMR    krd = max(3, kord)
743        W(i,  2,k) = W(i,  1,k)    DO ic = 1, nc
744        W(i,JMR,k) = W(i,JNP,k)  
745  5000  continue      ! ****6***0*********0*********0*********0*********0*********0**********72
746        endif  
747  C      CALL fzppm(imr, jnp, nlay, j1, dq(1,1,1,ic), w, q(1,1,1,ic), wk1, dpi, &
748        RETURN        dc2, crx, cry, pu, pv, xmass, ymass, delp1, krd)
749        END  
750  C  
751  C****6***0*********0*********0*********0*********0*********0**********72      IF (fill) CALL qckxyz(dq(1,1,1,ic), dc2, imr, jnp, nlay, j1, j2, cosp, &
752        subroutine FZPPM(IMR,JNP,NLAY,j1,DQ,WZ,P,DC,DQDT,AR,AL,A6,        acosp, .FALSE., ic, nstep)
753       &                 flux,wk1,wk2,wz2,delp,KORD)  
754        parameter ( kmax = 150 )      ! Recover tracer mixing ratio from "density" using predicted
755        parameter ( R23 = 2./3., R3 = 1./3.)      ! "air density" (pressure thickness) at time-level n+1
756        real WZ(IMR,JNP,NLAY),P(IMR,JNP,NLAY),DC(IMR,JNP,NLAY),  
757       &     wk1(IMR,*),delp(IMR,JNP,NLAY),DQ(IMR,JNP,NLAY),      DO k = 1, nlay
758       &     DQDT(IMR,JNP,NLAY)        DO j = 1, jnp
759  C Assuming JNP >= NLAY          DO i = 1, imr
760        real AR(IMR,*),AL(IMR,*),A6(IMR,*),flux(IMR,*),wk2(IMR,*),            q(i, j, k, ic) = dq(i, j, k, ic)/delp2(i, j, k)
761       &     wz2(IMR,*)          END DO
762  C        END DO
763        JMR = JNP - 1      END DO
764        IMJM = IMR*JNP  
765        NLAYM1 = NLAY - 1      IF (j1/=2) THEN
766  C        DO k = 1, nlay
767        LMT = KORD - 3          DO i = 1, imr
768  C            ! j=1 c'est le pôle Sud, j=JNP c'est le pôle Nord
769  C ****6***0*********0*********0*********0*********0*********0**********72            q(i, 2, k, ic) = q(i, 1, k, ic)
770  C Compute DC for PPM            q(i, jmr, k, ic) = q(i, jmp, k, ic)
771  C ****6***0*********0*********0*********0*********0*********0**********72          END DO
772  C        END DO
773        do 1000 k=1,NLAYM1      END IF
774        do 1000 i=1,IMJM    END DO
775        DQDT(i,1,k) = P(i,1,k+1) - P(i,1,k)  
776  1000  continue    IF (j1/=2) THEN
777  C      DO k = 1, nlay
778        DO 1220 k=2,NLAYM1        DO i = 1, imr
779        DO 1220 I=1,IMJM              w(i, 2, k) = w(i, 1, k)
780         c0 =  delp(i,1,k) / (delp(i,1,k-1)+delp(i,1,k)+delp(i,1,k+1))          w(i, jmr, k) = w(i, jnp, k)
781         c1 = (delp(i,1,k-1)+0.5*delp(i,1,k))/(delp(i,1,k+1)+delp(i,1,k))            END DO
782         c2 = (delp(i,1,k+1)+0.5*delp(i,1,k))/(delp(i,1,k-1)+delp(i,1,k))      END DO
783        tmp = c0*(c1*DQDT(i,1,k) + c2*DQDT(i,1,k-1))    END IF
784        Qmax = max(P(i,1,k-1),P(i,1,k),P(i,1,k+1)) - P(i,1,k)  
785        Qmin = P(i,1,k) - min(P(i,1,k-1),P(i,1,k),P(i,1,k+1))    RETURN
786        DC(i,1,k) = sign(min(abs(tmp),Qmax,Qmin), tmp)    END SUBROUTINE ppm3d
787  1220  CONTINUE  
788        ! ****6***0*********0*********0*********0*********0*********0**********72
789  C      SUBROUTINE fzppm(imr, jnp, nlay, j1, dq, wz, p, dc, dqdt, ar, al, a6, flux, &
790  C ****6***0*********0*********0*********0*********0*********0**********72      wk1, wk2, wz2, delp, kord)
791  C Loop over latitudes  (to save memory)    PARAMETER (kmax=150)
792  C ****6***0*********0*********0*********0*********0*********0**********72    PARAMETER (r23=2./3., r3=1./3.)
793  C    REAL wz(imr, jnp, nlay), p(imr, jnp, nlay), dc(imr, jnp, nlay), &
794        DO 2000 j=1,JNP      wk1(imr, *), delp(imr, jnp, nlay), dq(imr, jnp, nlay), &
795        if((j.eq.2 .or. j.eq.JMR) .and. j1.ne.2) goto 2000      dqdt(imr, jnp, nlay)
796  C    ! Assuming JNP >= NLAY
797        DO k=1,NLAY    REAL ar(imr, *), al(imr, *), a6(imr, *), flux(imr, *), wk2(imr, *), &
798        DO i=1,IMR      wz2(imr, *)
799        wz2(i,k) =   WZ(i,j,k)  
800        wk1(i,k) =    P(i,j,k)    jmr = jnp - 1
801        wk2(i,k) = delp(i,j,k)    imjm = imr*jnp
802        flux(i,k) = DC(i,j,k)  !this flux is actually the monotone slope    nlaym1 = nlay - 1
803        enddo  
804        enddo    lmt = kord - 3
805  C  
806  C****6***0*********0*********0*********0*********0*********0**********72    ! ****6***0*********0*********0*********0*********0*********0**********72
807  C Compute first guesses at cell interfaces    ! Compute DC for PPM
808  C First guesses are required to be continuous.    ! ****6***0*********0*********0*********0*********0*********0**********72
809  C ****6***0*********0*********0*********0*********0*********0**********72  
810  C    DO k = 1, nlaym1
811  C three-cell parabolic subgrid distribution at model top      DO i = 1, imjm
812  C two-cell parabolic with zero gradient subgrid distribution        dqdt(i, 1, k) = p(i, 1, k+1) - p(i, 1, k)
813  C at the surface.      END DO
814  C    END DO
815  C First guess top edge value  
816        DO 10 i=1,IMR    DO k = 2, nlaym1
817  C three-cell PPM      DO i = 1, imjm
818  C Compute a,b, and c of q = aP**2 + bP + c using cell averages and delp        c0 = delp(i, 1, k)/(delp(i,1,k-1)+delp(i,1,k)+delp(i,1,k+1))
819        a = 3.*( DQDT(i,j,2) - DQDT(i,j,1)*(wk2(i,2)+wk2(i,3))/        c1 = (delp(i,1,k-1)+0.5*delp(i,1,k))/(delp(i,1,k+1)+delp(i,1,k))
820       &         (wk2(i,1)+wk2(i,2)) ) /        c2 = (delp(i,1,k+1)+0.5*delp(i,1,k))/(delp(i,1,k-1)+delp(i,1,k))
821       &       ( (wk2(i,2)+wk2(i,3))*(wk2(i,1)+wk2(i,2)+wk2(i,3)) )        tmp = c0*(c1*dqdt(i,1,k)+c2*dqdt(i,1,k-1))
822        b = 2.*DQDT(i,j,1)/(wk2(i,1)+wk2(i,2)) -        qmax = max(p(i,1,k-1), p(i,1,k), p(i,1,k+1)) - p(i, 1, k)
823       &    R23*a*(2.*wk2(i,1)+wk2(i,2))        qmin = p(i, 1, k) - min(p(i,1,k-1), p(i,1,k), p(i,1,k+1))
824        AL(i,1) =  wk1(i,1) - wk2(i,1)*(R3*a*wk2(i,1) + 0.5*b)        dc(i, 1, k) = sign(min(abs(tmp),qmax,qmin), tmp)
825        AL(i,2) =  wk2(i,1)*(a*wk2(i,1) + b) + AL(i,1)      END DO
826  C    END DO
827  C Check if change sign  
828        if(wk1(i,1)*AL(i,1).le.0.) then  
829                   AL(i,1) = 0.    ! ****6***0*********0*********0*********0*********0*********0**********72
830               flux(i,1) = 0.    ! Loop over latitudes  (to save memory)
831          else    ! ****6***0*********0*********0*********0*********0*********0**********72
832               flux(i,1) =  wk1(i,1) - AL(i,1)  
833          endif    DO j = 1, jnp
834  10    continue      IF ((j==2 .OR. j==jmr) .AND. j1/=2) GO TO 2000
835  C  
836  C Bottom      DO k = 1, nlay
837        DO 15 i=1,IMR        DO i = 1, imr
838  C 2-cell PPM with zero gradient right at the surface          wz2(i, k) = wz(i, j, k)
839  C          wk1(i, k) = p(i, j, k)
840        fct = DQDT(i,j,NLAYM1)*wk2(i,NLAY)**2 /          wk2(i, k) = delp(i, j, k)
841       & ( (wk2(i,NLAY)+wk2(i,NLAYM1))*(2.*wk2(i,NLAY)+wk2(i,NLAYM1)))          flux(i, k) = dc(i, j, k) !this flux is actually the monotone slope
842        AR(i,NLAY) = wk1(i,NLAY) + fct        END DO
843        AL(i,NLAY) = wk1(i,NLAY) - (fct+fct)      END DO
844        if(wk1(i,NLAY)*AR(i,NLAY).le.0.) AR(i,NLAY) = 0.  
845        flux(i,NLAY) = AR(i,NLAY) -  wk1(i,NLAY)      ! ****6***0*********0*********0*********0*********0*********0**********72
846  15    continue      ! Compute first guesses at cell interfaces
847            ! First guesses are required to be continuous.
848  C      ! ****6***0*********0*********0*********0*********0*********0**********72
849  C****6***0*********0*********0*********0*********0*********0**********72  
850  C 4th order interpolation in the interior.      ! three-cell parabolic subgrid distribution at model top
851  C****6***0*********0*********0*********0*********0*********0**********72      ! two-cell parabolic with zero gradient subgrid distribution
852  C      ! at the surface.
853        DO 14 k=3,NLAYM1  
854        DO 12 i=1,IMR      ! First guess top edge value
855        c1 =  DQDT(i,j,k-1)*wk2(i,k-1) / (wk2(i,k-1)+wk2(i,k))      DO i = 1, imr
856        c2 =  2. / (wk2(i,k-2)+wk2(i,k-1)+wk2(i,k)+wk2(i,k+1))        ! three-cell PPM
857        A1   =  (wk2(i,k-2)+wk2(i,k-1)) / (2.*wk2(i,k-1)+wk2(i,k))        ! Compute a,b, and c of q = aP**2 + bP + c using cell averages and delp
858        A2   =  (wk2(i,k  )+wk2(i,k+1)) / (2.*wk2(i,k)+wk2(i,k-1))        a = 3.*(dqdt(i,j,2)-dqdt(i,j,1)*(wk2(i,2)+wk2(i,3))/(wk2(i,1)+wk2(i, &
859        AL(i,k) = wk1(i,k-1) + c1 + c2 *          2)))/((wk2(i,2)+wk2(i,3))*(wk2(i,1)+wk2(i,2)+wk2(i,3)))
860       &        ( wk2(i,k  )*(c1*(A1 - A2)+A2*flux(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       &          wk2(i,k-1)*A1*flux(i,k)  )        al(i, 1) = wk1(i, 1) - wk2(i, 1)*(r3*a*wk2(i,1)+0.5*b)
862  12    CONTINUE        al(i, 2) = wk2(i, 1)*(a*wk2(i,1)+b) + al(i, 1)
863  14    continue  
864  C        ! Check if change sign
865        do 20 i=1,IMR*NLAYM1        IF (wk1(i,1)*al(i,1)<=0.) THEN
866        AR(i,1) = AL(i,2)          al(i, 1) = 0.
867  20    continue          flux(i, 1) = 0.
 C  
       do 30 i=1,IMR*NLAY  
       A6(i,1) = 3.*(wk1(i,1)+wk1(i,1) - (AL(i,1)+AR(i,1)))  
 30    continue  
 C  
 C****6***0*********0*********0*********0*********0*********0**********72  
 C Top & Bot always monotonic  
       call lmtppm(flux(1,1),A6(1,1),AR(1,1),AL(1,1),wk1(1,1),IMR,0)  
       call lmtppm(flux(1,NLAY),A6(1,NLAY),AR(1,NLAY),AL(1,NLAY),  
      &            wk1(1,NLAY),IMR,0)  
 C  
 C Interior depending on KORD  
       if(LMT.LE.2)  
      &  call lmtppm(flux(1,2),A6(1,2),AR(1,2),AL(1,2),wk1(1,2),  
      &              IMR*(NLAY-2),LMT)  
 C  
 C****6***0*********0*********0*********0*********0*********0**********72  
 C  
       DO 140 i=1,IMR*NLAYM1  
       IF(wz2(i,1).GT.0.) then  
         CM = wz2(i,1) / wk2(i,1)  
         flux(i,2) = AR(i,1)+0.5*CM*(AL(i,1)-AR(i,1)+A6(i,1)*(1.-R23*CM))  
       else  
         CP= wz2(i,1) / wk2(i,2)          
         flux(i,2) = AL(i,2)+0.5*CP*(AL(i,2)-AR(i,2)-A6(i,2)*(1.+R23*CP))  
       endif  
 140   continue  
 C  
       DO 250 i=1,IMR*NLAYM1  
       flux(i,2) = wz2(i,1) * flux(i,2)  
 250   continue  
 C  
       do 350 i=1,IMR  
       DQ(i,j,   1) = DQ(i,j,   1) - flux(i,   2)  
       DQ(i,j,NLAY) = DQ(i,j,NLAY) + flux(i,NLAY)  
 350   continue  
 C  
       do 360 k=2,NLAYM1  
       do 360 i=1,IMR  
 360   DQ(i,j,k) = DQ(i,j,k) + flux(i,k) - flux(i,k+1)  
 2000  continue  
       return  
       end  
 C  
       subroutine xtp(IMR,JNP,IML,j1,j2,JN,JS,PU,DQ,Q,UC,  
      &               fx1,xmass,IORD)  
       dimension UC(IMR,*),DC(-IML:IMR+IML+1),xmass(IMR,JNP)  
      &    ,fx1(IMR+1),DQ(IMR,JNP),qtmp(-IML:IMR+1+IML)  
       dimension PU(IMR,JNP),Q(IMR,JNP),ISAVE(IMR)  
 C  
       IMP = IMR + 1  
 C  
 C van Leer at high latitudes  
       jvan = max(1,JNP/18)  
       j1vl = j1+jvan  
       j2vl = j2-jvan  
 C  
       do 1310 j=j1,j2  
 C  
       do i=1,IMR  
       qtmp(i) = q(i,j)  
       enddo  
 C  
       if(j.ge.JN .or. j.le.JS) goto 2222  
 C ************* Eulerian **********  
 C  
       qtmp(0)     = q(IMR,J)  
       qtmp(-1)    = q(IMR-1,J)  
       qtmp(IMP)   = q(1,J)  
       qtmp(IMP+1) = q(2,J)  
 C  
       IF(IORD.eq.1 .or. j.eq.j1. or. j.eq.j2) THEN  
       DO 1406 i=1,IMR  
       iu = float(i) - uc(i,j)  
 1406  fx1(i) = qtmp(iu)  
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  
   
       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    
1152        IF(IAD.eq.2) then      IF (jord<=0 .OR. jord>=3) THEN
1153        do j=j1-1,j2+1  
1154        do i=1,IMR        CALL fyppm(vc, p, dc2, fx, imr, jnp, j1, j2, a6, ar, al, jord)
1155        JP = NINT(VA(i,j))        
1156        rv = JP - VA(i,j)      ELSE
1157        JP = j - JP        DO i = 1, len
1158        a1 = 0.5*(wk(i,jp+1)+wk(i,jp-1)) - wk(i,jp)          jt = float(j1) - vc(i, j1)
1159        b1 = 0.5*(wk(i,jp+1)-wk(i,jp-1))          fx(i, j1) = p(i, jt) + (sign(1.,vc(i,j1))-vc(i,j1))*dc2(i, jt)
1160        ady(i,j) = wk(i,jp) + rv*(a1*rv + b1) - wk(i,j)        END DO
1161        enddo      END IF
1162        enddo    END IF
1163    
1164        ELSEIF(IAD.eq.1) then    DO i = 1, len
1165          do j=j1-1,j2+1      fx(i, j1) = fx(i, j1)*ymass(i, j1)
1166        do i=1,imr    END DO
1167        JP = float(j)-VA(i,j)  
1168        ady(i,j) = VA(i,j)*(wk(i,jp)-wk(i,jp+1))    DO j = j1, j2
1169        enddo      DO i = 1, imr
1170        enddo        dq(i, j) = dq(i, j) + (fx(i,j)-fx(i,j+1))*acosp(j)
1171        ENDIF      END DO
1172  C    END DO
1173          if(j1.ne.2) then  
1174          sum1 = 0.    ! Poles
1175          sum2 = 0.    sum1 = fx(imr, j1)
1176        do i=1,imr    sum2 = fx(imr, j2+1)
1177        sum1 = sum1 + ady(i,2)    DO i = 1, imr - 1
1178        sum2 = sum2 + ady(i,JMR)      sum1 = sum1 + fx(i, j1)
1179        enddo      sum2 = sum2 + fx(i, j2+1)
1180          sum1 = sum1 / IMR    END DO
1181          sum2 = sum2 / IMR  
1182  C    sum1 = dq(1, 1) - sum1*rcap
1183        do i=1,imr    sum2 = dq(1, jnp) + sum2*rcap
1184        ady(i,  2) =  sum1    DO i = 1, imr
1185        ady(i,JMR) =  sum2      dq(i, 1) = sum1
1186        ady(i,  1) =  sum1      dq(i, jnp) = sum2
1187        ady(i,JNP) =  sum2    END DO
1188        enddo  
1189          else    IF (j1/=2) THEN
1190  C Poles:      DO i = 1, imr
1191          sum1 = 0.        dq(i, 2) = sum1
1192          sum2 = 0.        dq(i, jmr) = sum2
1193        do i=1,imr      END DO
1194        sum1 = sum1 + ady(i,1)    END IF
1195        sum2 = sum2 + ady(i,JNP)  
1196        enddo    RETURN
1197          sum1 = sum1 / IMR  END SUBROUTINE ytp
1198          sum2 = sum2 / IMR  
1199  C  SUBROUTINE ymist(imr, jnp, j1, p, dc, id)
1200        do i=1,imr    PARAMETER (r24=1./24.)
1201        ady(i,  1) =  sum1    DIMENSION p(imr, jnp), dc(imr, jnp)
1202        ady(i,JNP) =  sum2  
1203        enddo    imh = imr/2
1204          endif    jmr = jnp - 1
1205  C    ijm3 = imr*(jmr-3)
1206          return  
1207          end    IF (id==2) THEN
1208  C      DO i = 1, imr*(jmr-1)
1209          subroutine xadv(IMR,JNP,j1,j2,p,UA,JS,JN,IML,adx,IAD)        tmp = 0.25*(p(i,3)-p(i,1))
1210          REAL p(IMR,JNP),adx(IMR,JNP),qtmp(-IMR:IMR+IMR),UA(IMR,JNP)        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          JMR = JNP-1        dc(i, 2) = sign(min(abs(tmp),pmin,pmax), tmp)
1213        do 1309 j=j1,j2      END DO
1214        if(J.GT.JS  .and. J.LT.JN) GO TO 1309    ELSE
1215  C      DO i = 1, imh
1216        do i=1,IMR        ! J=2
1217        qtmp(i) = p(i,j)        tmp = (8.*(p(i,3)-p(i,1))+p(i+imh,2)-p(i,4))*r24
1218        enddo        pmax = max(p(i,1), p(i,2), p(i,3)) - p(i, 2)
1219  C        pmin = p(i, 2) - min(p(i,1), p(i,2), p(i,3))
1220        do i=-IML,0        dc(i, 2) = sign(min(abs(tmp),pmin,pmax), tmp)
1221        qtmp(i)       = p(IMR+i,j)        ! J=JMR
1222        qtmp(IMR+1-i) = p(1-i,j)        tmp = (8.*(p(i,jnp)-p(i,jmr-1))+p(i,jmr-2)-p(i+imh,jmr))*r24
1223        enddo        pmax = max(p(i,jmr-1), p(i,jmr), p(i,jnp)) - p(i, jmr)
1224  C        pmin = p(i, jmr) - min(p(i,jmr-1), p(i,jmr), p(i,jnp))
1225        IF(IAD.eq.2) THEN        dc(i, jmr) = sign(min(abs(tmp),pmin,pmax), tmp)
1226        DO i=1,IMR      END DO
1227        IP = NINT(UA(i,j))      DO i = imh + 1, imr
1228        ru = IP - UA(i,j)        ! J=2
1229        IP = i - IP        tmp = (8.*(p(i,3)-p(i,1))+p(i-imh,2)-p(i,4))*r24
1230        a1 = 0.5*(qtmp(ip+1)+qtmp(ip-1)) - qtmp(ip)        pmax = max(p(i,1), p(i,2), p(i,3)) - p(i, 2)
1231        b1 = 0.5*(qtmp(ip+1)-qtmp(ip-1))        pmin = p(i, 2) - min(p(i,1), p(i,2), p(i,3))
1232        adx(i,j) = qtmp(ip) + ru*(a1*ru + b1)        dc(i, 2) = sign(min(abs(tmp),pmin,pmax), tmp)
1233        enddo        ! J=JMR
1234        ELSEIF(IAD.eq.1) then        tmp = (8.*(p(i,jnp)-p(i,jmr-1))+p(i,jmr-2)-p(i-imh,jmr))*r24
1235        DO i=1,IMR        pmax = max(p(i,jmr-1), p(i,jmr), p(i,jnp)) - p(i, jmr)
1236        iu = UA(i,j)        pmin = p(i, jmr) - min(p(i,jmr-1), p(i,jmr), p(i,jnp))
1237        ru = UA(i,j) - iu        dc(i, jmr) = sign(min(abs(tmp),pmin,pmax), tmp)
1238        iiu = i-iu      END DO
1239        if(UA(i,j).GE.0.) then  
1240        adx(i,j) = qtmp(iiu)+ru*(qtmp(iiu-1)-qtmp(iiu))      DO i = 1, ijm3
1241        else        tmp = (8.*(p(i,4)-p(i,2))+p(i,1)-p(i,5))*r24
1242        adx(i,j) = qtmp(iiu)+ru*(qtmp(iiu)-qtmp(iiu+1))        pmax = max(p(i,2), p(i,3), p(i,4)) - p(i, 3)
1243        endif        pmin = p(i, 3) - min(p(i,2), p(i,3), p(i,4))
1244        enddo        dc(i, 3) = sign(min(abs(tmp),pmin,pmax), tmp)
1245        ENDIF      END DO
1246  C    END IF
1247        do i=1,IMR  
1248        adx(i,j) = adx(i,j) - p(i,j)    IF (j1/=2) THEN
1249        enddo      DO i = 1, imr
1250  1309  continue        dc(i, 1) = 0.
1251  C        dc(i, jnp) = 0.
1252  C Eulerian upwind      END DO
1253  C    ELSE
1254        do j=JS+1,JN-1      ! Determine slopes in polar caps for scalars!
1255  C  
1256        do i=1,IMR      DO i = 1, imh
1257        qtmp(i) = p(i,j)        ! South
1258        enddo        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        qtmp(0)     = p(IMR,J)        pmin = p(i, 1) - min(p(i,2), p(i,1), p(i+imh,2))
1261        qtmp(IMR+1) = p(1,J)        dc(i, 1) = sign(min(abs(tmp),pmax,pmin), tmp)
1262  C        ! North.
1263        IF(IAD.eq.2) THEN        tmp = 0.25*(p(i+imh,jmr)-p(i,jmr))
1264        qtmp(-1)     = p(IMR-1,J)        pmax = max(p(i+imh,jmr), p(i,jnp), p(i,jmr)) - p(i, jnp)
1265        qtmp(IMR+2) = p(2,J)        pmin = p(i, jnp) - min(p(i+imh,jmr), p(i,jnp), p(i,jmr))
1266        do i=1,imr        dc(i, jnp) = sign(min(abs(tmp),pmax,pmin), tmp)
1267        IP = NINT(UA(i,j))      END DO
1268        ru = IP - UA(i,j)  
1269        IP = i - IP      DO i = imh + 1, imr
1270        a1 = 0.5*(qtmp(ip+1)+qtmp(ip-1)) - qtmp(ip)        dc(i, 1) = -dc(i-imh, 1)
1271        b1 = 0.5*(qtmp(ip+1)-qtmp(ip-1))        dc(i, jnp) = -dc(i-imh, jnp)
1272        adx(i,j) = qtmp(ip)- p(i,j) + ru*(a1*ru + b1)      END DO
1273        enddo    END IF
1274        ELSEIF(IAD.eq.1) then    RETURN
1275  C 1st order  END SUBROUTINE ymist
1276        DO i=1,IMR  
1277        IP = i - UA(i,j)  SUBROUTINE fyppm(vc, p, dc, flux, imr, jnp, j1, j2, a6, ar, al, jord)
1278        adx(i,j) = UA(i,j)*(qtmp(ip)-qtmp(ip+1))    PARAMETER (r3=1./3., r23=2./3.)
1279        enddo    REAL vc(imr, *), flux(imr, *), p(imr, *), dc(imr, *)
1280        ENDIF    ! Local work arrays.
1281        enddo    REAL ar(imr, jnp), al(imr, jnp), a6(imr, jnp)
1282  C    INTEGER lmt
1283          if(j1.ne.2) then    ! logical first
1284        do i=1,IMR    ! data first /.true./
1285        adx(i,  2) = 0.    ! SAVE LMT
1286        adx(i,JMR) = 0.  
1287        enddo    imh = imr/2
1288          endif    jmr = jnp - 1
1289  C set cross term due to x-adv at the poles to zero.    j11 = j1 - 1
1290        do i=1,IMR    imjm1 = imr*(j2-j1+2)
1291        adx(i,  1) = 0.    len = imr*(j2-j1+3)
1292        adx(i,JNP) = 0.    ! if(first) then
1293        enddo    ! IF(JORD.LE.0) then
1294          return    ! if(JMR.GE.90) then
1295          end    ! LMT = 0
1296  C    ! elseif(JMR.GE.45) then
1297        subroutine lmtppm(DC,A6,AR,AL,P,IM,LMT)    ! LMT = 1
1298  C    ! else
1299  C A6 =  CURVATURE OF THE TEST PARABOLA    ! LMT = 2
1300  C AR =  RIGHT EDGE VALUE OF THE TEST PARABOLA    ! endif
1301  C AL =  LEFT  EDGE VALUE OF THE TEST PARABOLA    ! else
1302  C DC =  0.5 * MISMATCH    ! LMT = JORD - 3
1303  C P  =  CELL-AVERAGED VALUE    ! endif
1304  C IM =  VECTOR LENGTH  
1305  C    ! first = .false.
1306  C OPTIONS:    ! endif
1307  C  
1308  C LMT = 0: FULL MONOTONICITY    ! modifs pour pouvoir choisir plusieurs schemas PPM
1309  C LMT = 1: SEMI-MONOTONIC CONSTRAINT (NO UNDERSHOOTS)    lmt = jord - 3
1310  C LMT = 2: POSITIVE-DEFINITE CONSTRAINT  
1311  C    DO i = 1, imr*jmr
1312        parameter ( R12 = 1./12. )      al(i, 2) = 0.5*(p(i,1)+p(i,2)) + (dc(i,1)-dc(i,2))*r3
1313        dimension A6(IM),AR(IM),AL(IM),P(IM),DC(IM)      ar(i, 1) = al(i, 2)
1314  C    END DO
1315        if(LMT.eq.0) then  
1316  C Full constraint    ! Poles:
1317        do 100 i=1,IM  
1318        if(DC(i).eq.0.) then    DO i = 1, imh
1319              AR(i) = p(i)      al(i, 1) = al(i+imh, 2)
1320              AL(i) = p(i)      al(i+imh, 1) = al(i, 2)
1321              A6(i) = 0.  
1322        else      ar(i, jnp) = ar(i+imh, jmr)
1323        da1  = AR(i) - AL(i)      ar(i+imh, jnp) = ar(i, jmr)
1324        da2  = da1**2    END DO
1325        A6DA = A6(i)*da1  
1326        if(A6DA .lt. -da2) then    ! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1327              A6(i) = 3.*(AL(i)-p(i))    ! Rajout pour LMDZ.3.3
1328              AR(i) = AL(i) - A6(i)    ! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1329        elseif(A6DA .gt. da2) then    ar(imr, 1) = al(1, 1)
1330              A6(i) = 3.*(AR(i)-p(i))    ar(imr, jnp) = al(1, jnp)
1331              AL(i) = AR(i) - A6(i)    ! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1332        endif  
1333        endif  
1334  100   continue    DO i = 1, len
1335        elseif(LMT.eq.1) then      a6(i, j11) = 3.*(p(i,j11)+p(i,j11)-(al(i,j11)+ar(i,j11)))
1336  C Semi-monotonic constraint    END DO
1337        do 150 i=1,IM  
1338        if(abs(AR(i)-AL(i)) .GE. -A6(i)) go to 150    IF (lmt<=2) CALL lmtppm(dc(1,j11), a6(1,j11), ar(1,j11), al(1,j11), &
1339        if(p(i).lt.AR(i) .and. p(i).lt.AL(i)) then      p(1,j11), len, lmt)
1340              AR(i) = p(i)  
1341              AL(i) = p(i)  
1342              A6(i) = 0.    DO i = 1, imjm1
1343        elseif(AR(i) .gt. AL(i)) then      IF (vc(i,j1)>0.) THEN
1344              A6(i) = 3.*(AL(i)-p(i))        flux(i, j1) = ar(i, j11) + 0.5*vc(i, j1)*(al(i,j11)-ar(i,j11)+a6(i,j11) &
1345              AR(i) = AL(i) - A6(i)          *(1.-r23*vc(i,j1)))
1346        else      ELSE
1347              A6(i) = 3.*(AR(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              AL(i) = AR(i) - A6(i)          +r23*vc(i,j1)))
1349        endif      END IF
1350  150   continue    END DO
1351        elseif(LMT.eq.2) then    RETURN
1352        do 250 i=1,IM  END SUBROUTINE fyppm
1353        if(abs(AR(i)-AL(i)) .GE. -A6(i)) go to 250  
1354        fmin = p(i) + 0.25*(AR(i)-AL(i))**2/A6(i) + A6(i)*R12  SUBROUTINE yadv(imr, jnp, j1, j2, p, va, ady, wk, iad)
1355        if(fmin.ge.0.) go to 250    REAL p(imr, jnp), ady(imr, jnp), va(imr, jnp)
1356        if(p(i).lt.AR(i) .and. p(i).lt.AL(i)) then    REAL wk(imr, -1:jnp+2)
1357              AR(i) = p(i)  
1358              AL(i) = p(i)    jmr = jnp - 1
1359              A6(i) = 0.    imh = imr/2
1360        elseif(AR(i) .gt. AL(i)) then    DO j = 1, jnp
1361              A6(i) = 3.*(AL(i)-p(i))      DO i = 1, imr
1362              AR(i) = AL(i) - A6(i)        wk(i, j) = p(i, j)
1363        else      END DO
1364              A6(i) = 3.*(AR(i)-p(i))    END DO
1365              AL(i) = AR(i) - A6(i)    ! Poles:
1366        endif    DO i = 1, imh
1367  250   continue      wk(i, -1) = p(i+imh, 3)
1368        endif      wk(i+imh, -1) = p(i, 3)
1369        return      wk(i, 0) = p(i+imh, 2)
1370        end      wk(i+imh, 0) = p(i, 2)
1371  C      wk(i, jnp+1) = p(i+imh, jmr)
1372        subroutine A2C(U,V,IMR,JMR,j1,j2,CRX,CRY,dtdx5,DTDY5)      wk(i+imh, jnp+1) = p(i, jmr)
1373        dimension U(IMR,*),V(IMR,*),CRX(IMR,*),CRY(IMR,*),DTDX5(*)      wk(i, jnp+2) = p(i+imh, jnp-2)
1374  C      wk(i+imh, jnp+2) = p(i, jnp-2)
1375        do 35 j=j1,j2    END DO
1376        do 35 i=2,IMR  
1377  35    CRX(i,J) = dtdx5(j)*(U(i,j)+U(i-1,j))    IF (iad==2) THEN
1378  C      DO j = j1 - 1, j2 + 1
1379        do 45 j=j1,j2        DO i = 1, imr
1380  45    CRX(1,J) = dtdx5(j)*(U(1,j)+U(IMR,j))          jp = nint(va(i,j))
1381  C          rv = jp - va(i, j)
1382        do 55 i=1,IMR*JMR          jp = j - jp
1383  55    CRY(i,2) = DTDY5*(V(i,2)+V(i,1))          a1 = 0.5*(wk(i,jp+1)+wk(i,jp-1)) - wk(i, jp)
1384        return          b1 = 0.5*(wk(i,jp+1)-wk(i,jp-1))
1385        end          ady(i, j) = wk(i, jp) + rv*(a1*rv+b1) - wk(i, j)
1386  C        END DO
1387        subroutine cosa(cosp,cose,JNP,PI,DP)      END DO
1388        dimension cosp(*),cose(*)  
1389        JMR = JNP-1    ELSE IF (iad==1) THEN
1390        do 55 j=2,JNP      DO j = j1 - 1, j2 + 1
1391          ph5  =  -0.5*PI + (FLOAT(J-1)-0.5)*DP        DO i = 1, imr
1392  55      cose(j) = cos(ph5)          jp = float(j) - va(i, j)
1393  C          ady(i, j) = va(i, j)*(wk(i,jp)-wk(i,jp+1))
1394        JEQ = (JNP+1) / 2        END DO
1395        if(JMR .eq. 2*(JMR/2) ) then      END DO
1396        do j=JNP, JEQ+1, -1    END IF
1397         cose(j) =  cose(JNP+2-j)  
1398        enddo    IF (j1/=2) THEN
1399        else      sum1 = 0.
1400  C cell edge at equator.      sum2 = 0.
1401         cose(JEQ+1) =  1.      DO i = 1, imr
1402        do j=JNP, JEQ+2, -1        sum1 = sum1 + ady(i, 2)
1403         cose(j) =  cose(JNP+2-j)        sum2 = sum2 + ady(i, jmr)
1404         enddo      END DO
1405        endif      sum1 = sum1/imr
1406  C      sum2 = sum2/imr
1407        do 66 j=2,JMR  
1408  66    cosp(j) = 0.5*(cose(j)+cose(j+1))      DO i = 1, imr
1409        cosp(1) = 0.        ady(i, 2) = sum1
1410        cosp(JNP) = 0.        ady(i, jmr) = sum2
1411        return        ady(i, 1) = sum1
1412        end        ady(i, jnp) = sum2
1413  C      END DO
1414        subroutine cosc(cosp,cose,JNP,PI,DP)    ELSE
1415        dimension cosp(*),cose(*)      ! Poles:
1416  C      sum1 = 0.
1417        phi = -0.5*PI      sum2 = 0.
1418        do 55 j=2,JNP-1      DO i = 1, imr
1419        phi  =  phi + DP        sum1 = sum1 + ady(i, 1)
1420  55    cosp(j) = cos(phi)        sum2 = sum2 + ady(i, jnp)
1421          cosp(  1) = 0.      END DO
1422          cosp(JNP) = 0.      sum1 = sum1/imr
1423  C      sum2 = sum2/imr
1424        do 66 j=2,JNP  
1425          cose(j) = 0.5*(cosp(j)+cosp(j-1))      DO i = 1, imr
1426  66    CONTINUE        ady(i, 1) = sum1
1427  C        ady(i, jnp) = sum2
1428        do 77 j=2,JNP-1      END DO
1429         cosp(j) = 0.5*(cose(j)+cose(j+1))    END IF
1430  77    CONTINUE  
1431        return    RETURN
1432        end  END SUBROUTINE yadv
1433  C  
1434        SUBROUTINE qckxyz (Q,qtmp,IMR,JNP,NLAY,j1,j2,cosp,acosp,  SUBROUTINE xadv(imr, jnp, j1, j2, p, ua, js, jn, iml, adx, iad)
1435       &                   cross,IC,NSTEP)    REAL p(imr, jnp), adx(imr, jnp), qtmp(-imr:imr+imr), ua(imr, jnp)
1436  C  
1437        parameter( tiny = 1.E-60 )    jmr = jnp - 1
1438        DIMENSION Q(IMR,JNP,NLAY),qtmp(IMR,JNP),cosp(*),acosp(*)    DO j = j1, j2
1439        logical cross      IF (j>js .AND. j<jn) GO TO 1309
1440  C  
1441        NLAYM1 = NLAY-1      DO i = 1, imr
1442        len = IMR*(j2-j1+1)        qtmp(i) = p(i, j)
1443        ip = 0      END DO
1444  C  
1445  C Top layer      DO i = -iml, 0
1446        L = 1        qtmp(i) = p(imr+i, j)
1447          icr = 1        qtmp(imr+1-i) = p(1-i, j)
1448        call filns(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,ipy,tiny)      END DO
1449        if(ipy.eq.0) goto 50  
1450        call filew(q(1,1,L),qtmp,IMR,JNP,j1,j2,ipx,tiny)      IF (iad==2) THEN
1451        if(ipx.eq.0) goto 50        DO i = 1, imr
1452  C          ip = nint(ua(i,j))
1453        if(cross) then          ru = ip - ua(i, j)
1454        call filcr(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,icr,tiny)          ip = i - ip
1455        endif          a1 = 0.5*(qtmp(ip+1)+qtmp(ip-1)) - qtmp(ip)
1456        if(icr.eq.0) goto 50          b1 = 0.5*(qtmp(ip+1)-qtmp(ip-1))
1457  C          adx(i, j) = qtmp(ip) + ru*(a1*ru+b1)
1458  C Vertical filling...        END DO
1459        do i=1,len      ELSE IF (iad==1) THEN
1460        IF( Q(i,j1,1).LT.0.) THEN        DO i = 1, imr
1461        ip = ip + 1          iu = ua(i, j)
1462            Q(i,j1,2) = Q(i,j1,2) + Q(i,j1,1)          ru = ua(i, j) - iu
1463            Q(i,j1,1) = 0.          iiu = i - iu
1464        endif          IF (ua(i,j)>=0.) THEN
1465        enddo            adx(i, j) = qtmp(iiu) + ru*(qtmp(iiu-1)-qtmp(iiu))
1466  C          ELSE
1467  50    continue            adx(i, j) = qtmp(iiu) + ru*(qtmp(iiu)-qtmp(iiu+1))
1468        DO 225 L = 2,NLAYM1          END IF
1469        icr = 1        END DO
1470  C      END IF
1471        call filns(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,ipy,tiny)  
1472        if(ipy.eq.0) goto 225      DO i = 1, imr
1473        call filew(q(1,1,L),qtmp,IMR,JNP,j1,j2,ipx,tiny)        adx(i, j) = adx(i, j) - p(i, j)
1474        if(ipx.eq.0) go to 225      END DO
1475        if(cross) then  1309 END DO
1476        call filcr(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,icr,tiny)  
1477        endif    ! Eulerian upwind
1478        if(icr.eq.0) goto 225  
1479  C    DO j = js + 1, jn - 1
1480        do i=1,len  
1481        IF( Q(I,j1,L).LT.0.) THEN      DO i = 1, imr
1482  C        qtmp(i) = p(i, j)
1483        END DO
1484    
1485        qtmp(0) = p(imr, j)
1486        qtmp(imr+1) = p(1, j)
1487    
1488        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            d0 = max(0., qtmp(j,i+1))
1974            d2 = min(-qtmp(j,i), d0)
1975            qtmp(j, i+1) = qtmp(j, i+1) - d2
1976            qtmp(j, i) = qtmp(j, i) + d2 + tiny
1977          END IF
1978        END DO
1979      END DO
1980    
1981      i = 1
1982      DO j = j1, j2
1983        IF (qtmp(j,i)<0.) THEN
1984          ipx = 1
1985          ! west
1986          d0 = max(0., qtmp(j,imr))
1987          d1 = min(-qtmp(j,i), d0)
1988          qtmp(j, imr) = qtmp(j, imr) - d1
1989          qtmp(j, i) = qtmp(j, i) + d1
1990          ! east
1991          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.80  
changed lines
  Added in v.81
  ViewVC Help
Powered by ViewVC 1.1.21