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