My Project
 All Classes Files Functions Variables Macros
ppm3d.F
Go to the documentation of this file.
1 !
2 ! $Id: ppm3d.F 1403 2010-07-01 09:02:53Z fairhead $
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 c 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  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  SAVE dtdy, dtdy5, rcap, js0, jn0, iml,
302  & dtdx, dtdx5, acosp, cosp, cose, dap,dbk
303 C
304 
305  jmr = jnp -1
306  imjm = imr*jnp
307  j2 = jnp - j1 + 1
308  nstep = nstep + 1
309 C
310 C *********** Initialization **********************
311  if(nstep.eq.1) then
312 c
313  write(6,*) '------------------------------------ '
314  write(6,*) 'NASA/GSFC Transport Core Version 4.5'
315  write(6,*) '------------------------------------ '
316 c
317  WRITE(6,*) 'IMR=',imr,' JNP=',jnp,' NLAY=',nlay,' j1=',j1
318  WRITE(6,*) 'NC=',nc,iord,jord,kord,ndt
319 C
320 C controles sur les parametres
321  if(nlay.LT.6) then
322  write(6,*) 'NLAY must be >= 6'
323  stop
324  endif
325  if (jnp.LT.nlay) then
326  write(6,*) 'JNP must be >= NLAY'
327  stop
328  endif
329  imrd2=mod(imr,2)
330  if (j1.eq.2.and.imrd2.NE.0) then
331  write(6,*) 'if j1=2 IMR must be an even integer'
332  stop
333  endif
334 
335 C
336  if(jmax.lt.jnp .or. kmax.lt.nlay) then
337  write(6,*) 'Jmax or Kmax is too small'
338  stop
339  endif
340 C
341  DO k=1,nlay
342  dap(k) = (ap(k+1) - ap(k))*pt
343  dbk(k) = bp(k+1) - bp(k)
344  ENDDO
345 C
346  pi = 4. * atan(1.)
347  dl = 2.*pi / REAL(imr)
348  dp = pi / REAL(jmr)
349 C
350  if(igd.eq.0) then
351 C Compute analytic cosine at cell edges
352  call cosa(cosp,cose,jnp,pi,dp)
353  else
354 C Define cosine consistent with GEOS-GCM (using dycore2.0 or later)
355  call cosc(cosp,cose,jnp,pi,dp)
356  endif
357 C
358  do 15 j=2,jmr
359 15 acosp(j) = 1. / cosp(j)
360 C
361 C Inverse of the Scaled polar cap area.
362 C
363  rcap = dp / (imr*(1.-cos((j1-1.5)*dp)))
364  acosp(1) = rcap
365  acosp(jnp) = rcap
366  endif
367 C
368  if(ndt0 .ne. ndt) then
369  dt = ndt
370  ndt0 = ndt
371 
372  if(umax .lt. 180.) then
373  write(6,*) 'Umax may be too small!'
374  endif
375  cr1 = abs(umax*dt)/(dl*ae)
376  maxdt = dp*ae / abs(umax) + 0.5
377  write(6,*)'Largest time step for max(V)=',umax,' is ',maxdt
378  if(maxdt .lt. abs(ndt)) then
379  write(6,*) 'Warning!!! NDT maybe too large!'
380  endif
381 C
382  if(cr1.ge.0.95) then
383  js0 = 0
384  jn0 = 0
385  iml = imr-2
386  ztc = 0.
387  else
388  ztc = acos(cr1) * (180./pi)
389 C
390  js0 = REAL(jmr)*(90.-ztc)/180. + 2
391  js0 = max(js0, j1+1)
392  iml = min(6*js0/(j1-1)+2, 4*imr/5)
393  jn0 = jnp-js0+1
394  endif
395 C
396 C
397  do j=2,jmr
398  dtdx(j) = dt / ( dl*ae*cosp(j) )
399 
400 c print*,'dtdx=',dtdx(j)
401  dtdx5(j) = 0.5*dtdx(j)
402  enddo
403 C
404 
405  dtdy = dt /(ae*dp)
406 c print*,'dtdy=',dtdy
407  dtdy5 = 0.5*dtdy
408 C
409 c write(6,*) 'J1=',J1,' J2=', J2
410  endif
411 C
412 C *********** End Initialization **********************
413 C
414 C delp = pressure thickness: the psudo-density in a hydrostatic system.
415  do k=1,nlay
416  do j=1,jnp
417  do i=1,imr
418  delp1(i,j,k)=dap(k)+dbk(k)*ps1(i,j)
419  delp2(i,j,k)=dap(k)+dbk(k)*ps2(i,j)
420  enddo
421  enddo
422  enddo
423 
424 C
425  if(j1.ne.2) then
426  DO 40 ic=1,nc
427  DO 40 l=1,nlay
428  DO 40 i=1,imr
429  q(i, 2,l,ic) = q(i, 1,l,ic)
430 40 q(i,jmr,l,ic) = q(i,jnp,l,ic)
431  endif
432 C
433 C Compute "tracer density"
434  DO 550 ic=1,nc
435  DO 44 k=1,nlay
436  DO 44 j=1,jnp
437  DO 44 i=1,imr
438 44 dq(i,j,k,ic) = q(i,j,k,ic)*delp1(i,j,k)
439 550 continue
440 C
441  do 1500 k=1,nlay
442 C
443  if(igd.eq.0) then
444 C Convert winds on A-Grid to Courant # on C-Grid.
445  call a2c(u(1,1,k),v(1,1,k),imr,jmr,j1,j2,crx,cry,dtdx5,dtdy5)
446  else
447 C Convert winds on C-grid to Courant #
448  do 45 j=j1,j2
449  do 45 i=2,imr
450 45 crx(i,j) = dtdx(j)*u(i-1,j,k)
451 
452 C
453  do 50 j=j1,j2
454 50 crx(1,j) = dtdx(j)*u(imr,j,k)
455 C
456  do 55 i=1,imr*jmr
457 55 cry(i,2) = dtdy*v(i,1,k)
458  endif
459 C
460 C Determine JS and JN
461  js = j1
462  jn = j2
463 C
464  do j=js0,j1+1,-1
465  do i=1,imr
466  if(abs(crx(i,j)).GT.1.) then
467  js = j
468  go to 2222
469  endif
470  enddo
471  enddo
472 C
473 2222 continue
474  do j=jn0,j2-1
475  do i=1,imr
476  if(abs(crx(i,j)).GT.1.) then
477  jn = j
478  go to 2233
479  endif
480  enddo
481  enddo
482 2233 continue
483 C
484  if(j1.ne.2) then ! Enlarged polar cap.
485  do i=1,imr
486  dpi(i, 2,k) = 0.
487  dpi(i,jmr,k) = 0.
488  enddo
489  endif
490 C
491 C ******* Compute horizontal mass fluxes ************
492 C
493 C N-S component
494  do j=j1,j2+1
495  d5 = 0.5 * cose(j)
496  do i=1,imr
497  ymass(i,j) = cry(i,j)*d5*(delp2(i,j,k) + delp2(i,j-1,k))
498  enddo
499  enddo
500 C
501  do 95 j=j1,j2
502  DO 95 i=1,imr
503 95 dpi(i,j,k) = (ymass(i,j) - ymass(i,j+1)) * acosp(j)
504 C
505 C Poles
506  sum1 = ymass(imr,j1 )
507  sum2 = ymass(imr,j2+1)
508  do i=1,imr-1
509  sum1 = sum1 + ymass(i,j1 )
510  sum2 = sum2 + ymass(i,j2+1)
511  enddo
512 C
513  sum1 = - sum1 * rcap
514  sum2 = sum2 * rcap
515  do i=1,imr
516  dpi(i, 1,k) = sum1
517  dpi(i,jnp,k) = sum2
518  enddo
519 C
520 C E-W component
521 C
522  do j=j1,j2
523  do i=2,imr
524  pu(i,j) = 0.5 * (delp2(i,j,k) + delp2(i-1,j,k))
525  enddo
526  enddo
527 C
528  do j=j1,j2
529  pu(1,j) = 0.5 * (delp2(1,j,k) + delp2(imr,j,k))
530  enddo
531 C
532  do 110 j=j1,j2
533  DO 110 i=1,imr
534 110 xmass(i,j) = pu(i,j)*crx(i,j)
535 C
536  DO 120 j=j1,j2
537  DO 120 i=1,imr-1
538 120 dpi(i,j,k) = dpi(i,j,k) + xmass(i,j) - xmass(i+1,j)
539 C
540  DO 130 j=j1,j2
541 130 dpi(imr,j,k) = dpi(imr,j,k) + xmass(imr,j) - xmass(1,j)
542 C
543  DO j=j1,j2
544  do i=1,imr-1
545  ua(i,j) = 0.5 * (crx(i,j)+crx(i+1,j))
546  enddo
547  enddo
548 C
549  DO j=j1,j2
550  ua(imr,j) = 0.5 * (crx(imr,j)+crx(1,j))
551  enddo
552 ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
553 c Rajouts pour LMDZ.3.3
554 ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
555  do i=1,imr
556  do j=1,jnp
557  va(i,j)=0.
558  enddo
559  enddo
560 
561  do i=1,imr*(jmr-1)
562  va(i,2) = 0.5*(cry(i,2)+cry(i,3))
563  enddo
564 C
565  if(j1.eq.2) then
566  imh = imr/2
567  do i=1,imh
568  va(i, 1) = 0.5*(cry(i,2)-cry(i+imh,2))
569  va(i+imh, 1) = -va(i,1)
570  va(i, jnp) = 0.5*(cry(i,jnp)-cry(i+imh,jmr))
571  va(i+imh,jnp) = -va(i,jnp)
572  enddo
573  va(imr,1)=va(1,1)
574  va(imr,jnp)=va(1,jnp)
575  endif
576 C
577 C ****6***0*********0*********0*********0*********0*********0**********72
578  do 1000 ic=1,nc
579 C
580  do i=1,imjm
581  wk1(i,1,1) = 0.
582  wk1(i,1,2) = 0.
583  enddo
584 C
585 C E-W advective cross term
586  do 250 j=j1,j2
587  if(j.GT.js .and. j.LT.jn) go to 250
588 C
589  do i=1,imr
590  qtmp(i) = q(i,j,k,ic)
591  enddo
592 C
593  do i=-iml,0
594  qtmp(i) = q(imr+i,j,k,ic)
595  qtmp(imr+1-i) = q(1-i,j,k,ic)
596  enddo
597 C
598  DO 230 i=1,imr
599  iu = ua(i,j)
600  ru = ua(i,j) - iu
601  iiu = i-iu
602  if(ua(i,j).GE.0.) then
603  wk1(i,j,1) = qtmp(iiu)+ru*(qtmp(iiu-1)-qtmp(iiu))
604  else
605  wk1(i,j,1) = qtmp(iiu)+ru*(qtmp(iiu)-qtmp(iiu+1))
606  endif
607  wk1(i,j,1) = wk1(i,j,1) - qtmp(i)
608 230 continue
609 250 continue
610 C
611  if(jn.ne.0) then
612  do j=js+1,jn-1
613 C
614  do i=1,imr
615  qtmp(i) = q(i,j,k,ic)
616  enddo
617 C
618  qtmp(0) = q(imr,j,k,ic)
619  qtmp(imr+1) = q( 1,j,k,ic)
620 C
621  do i=1,imr
622  iu = i - ua(i,j)
623  wk1(i,j,1) = ua(i,j)*(qtmp(iu) - qtmp(iu+1))
624  enddo
625  enddo
626  endif
627 C ****6***0*********0*********0*********0*********0*********0**********72
628 C Contribution from the N-S advection
629  do i=1,imr*(j2-j1+1)
630  jt = REAL(J1) - va(i,j1)
631  wk1(i,j1,2) = va(i,j1) * (q(i,jt,k,ic) - q(i,jt+1,k,ic))
632  enddo
633 C
634  do i=1,imjm
635  wk1(i,1,1) = q(i,1,k,ic) + 0.5*wk1(i,1,1)
636  wk1(i,1,2) = q(i,1,k,ic) + 0.5*wk1(i,1,2)
637  enddo
638 C
639  if(cross) then
640 C Add cross terms in the vertical direction.
641  if(iord .GE. 2) then
642  iad = 2
643  else
644  iad = 1
645  endif
646 C
647  if(jord .GE. 2) then
648  jad = 2
649  else
650  jad = 1
651  endif
652  call xadv(imr,jnp,j1,j2,wk1(1,1,2),ua,js,jn,iml,dc2,iad)
653  call yadv(imr,jnp,j1,j2,wk1(1,1,1),va,pv,w,jad)
654  do j=1,jnp
655  do i=1,imr
656  q(i,j,k,ic) = q(i,j,k,ic) + dc2(i,j) + pv(i,j)
657  enddo
658  enddo
659  endif
660 C
661  call xtp(imr,jnp,iml,j1,j2,jn,js,pu,dq(1,1,k,ic),wk1(1,1,2)
662  & ,crx,fx1,xmass,iord)
663 
664  call ytp(imr,jnp,j1,j2,acosp,rcap,dq(1,1,k,ic),wk1(1,1,1),cry,
665  & dc2,ymass,wk1(1,1,3),wk1(1,1,4),wk1(1,1,5),wk1(1,1,6),jord)
666 C
667 1000 continue
668 1500 continue
669 C
670 C ******* Compute vertical mass flux (same unit as PS) ***********
671 C
672 C 1st step: compute total column mass CONVERGENCE.
673 C
674  do 320 j=1,jnp
675  do 320 i=1,imr
676 320 cry(i,j) = dpi(i,j,1)
677 C
678  do 330 k=2,nlay
679  do 330 j=1,jnp
680  do 330 i=1,imr
681  cry(i,j) = cry(i,j) + dpi(i,j,k)
682 330 continue
683 C
684  do 360 j=1,jnp
685  do 360 i=1,imr
686 C
687 C 2nd step: compute PS2 (PS at n+1) using the hydrostatic assumption.
688 C Changes (increases) to surface pressure = total column mass convergence
689 C
690  ps2(i,j) = ps1(i,j) + cry(i,j)
691 C
692 C 3rd step: compute vertical mass flux from mass conservation principle.
693 C
694  w(i,j,1) = dpi(i,j,1) - dbk(1)*cry(i,j)
695  w(i,j,nlay) = 0.
696 360 continue
697 C
698  do 370 k=2,nlay-1
699  do 370 j=1,jnp
700  do 370 i=1,imr
701  w(i,j,k) = w(i,j,k-1) + dpi(i,j,k) - dbk(k)*cry(i,j)
702 370 continue
703 C
704  DO 380 k=1,nlay
705  DO 380 j=1,jnp
706  DO 380 i=1,imr
707  delp2(i,j,k) = dap(k) + dbk(k)*ps2(i,j)
708 380 continue
709 C
710  krd = max(3, kord)
711  do 4000 ic=1,nc
712 C
713 C****6***0*********0*********0*********0*********0*********0**********72
714 
715  call fzppm(imr,jnp,nlay,j1,dq(1,1,1,ic),w,q(1,1,1,ic),wk1,dpi,
716  & dc2,crx,cry,pu,pv,xmass,ymass,delp1,krd)
717 C
718 
719  if(fill) call qckxyz(dq(1,1,1,ic),dc2,imr,jnp,nlay,j1,j2,
720  & cosp,acosp,.false.,ic,nstep)
721 C
722 C Recover tracer mixing ratio from "density" using predicted
723 C "air density" (pressure thickness) at time-level n+1
724 C
725  DO k=1,nlay
726  DO j=1,jnp
727  DO i=1,imr
728  q(i,j,k,ic) = dq(i,j,k,ic) / delp2(i,j,k)
729 c print*,'i=',i,'j=',j,'k=',k,'Q(i,j,k,IC)=',Q(i,j,k,IC)
730  enddo
731  enddo
732  enddo
733 C
734  if(j1.ne.2) then
735  DO 400 k=1,nlay
736  DO 400 i=1,imr
737 c j=1 c'est le pôle Sud, j=JNP c'est le pôle Nord
738  q(i, 2,k,ic) = q(i, 1,k,ic)
739  q(i,jmr,k,ic) = q(i,jnp,k,ic)
740 400 CONTINUE
741  endif
742 4000 continue
743 C
744  if(j1.ne.2) then
745  DO 5000 k=1,nlay
746  DO 5000 i=1,imr
747  w(i, 2,k) = w(i, 1,k)
748  w(i,jmr,k) = w(i,jnp,k)
749 5000 continue
750  endif
751 C
752  RETURN
753  END
754 C
755 C****6***0*********0*********0*********0*********0*********0**********72
756  subroutine fzppm(IMR,JNP,NLAY,j1,DQ,WZ,P,DC,DQDT,AR,AL,A6,
757  & flux,wk1,wk2,wz2,delp,kord)
758  parameter( kmax = 150 )
759  parameter( r23 = 2./3., r3 = 1./3.)
760  real wz(imr,jnp,nlay),p(imr,jnp,nlay),dc(imr,jnp,nlay),
761  & wk1(imr,*),delp(imr,jnp,nlay),dq(imr,jnp,nlay),
762  & dqdt(imr,jnp,nlay)
763 C Assuming JNP >= NLAY
764  real ar(imr,*),al(imr,*),a6(imr,*),flux(imr,*),wk2(imr,*),
765  & wz2(imr,*)
766 C
767  jmr = jnp - 1
768  imjm = imr*jnp
769  nlaym1 = nlay - 1
770 C
771  lmt = kord - 3
772 C
773 C ****6***0*********0*********0*********0*********0*********0**********72
774 C Compute DC for PPM
775 C ****6***0*********0*********0*********0*********0*********0**********72
776 C
777  do 1000 k=1,nlaym1
778  do 1000 i=1,imjm
779  dqdt(i,1,k) = p(i,1,k+1) - p(i,1,k)
780 1000 continue
781 C
782  DO 1220 k=2,nlaym1
783  DO 1220 i=1,imjm
784  c0 = delp(i,1,k) / (delp(i,1,k-1)+delp(i,1,k)+delp(i,1,k+1))
785  c1 = (delp(i,1,k-1)+0.5*delp(i,1,k))/(delp(i,1,k+1)+delp(i,1,k))
786  c2 = (delp(i,1,k+1)+0.5*delp(i,1,k))/(delp(i,1,k-1)+delp(i,1,k))
787  tmp = c0*(c1*dqdt(i,1,k) + c2*dqdt(i,1,k-1))
788  qmax = max(p(i,1,k-1),p(i,1,k),p(i,1,k+1)) - p(i,1,k)
789  qmin = p(i,1,k) - min(p(i,1,k-1),p(i,1,k),p(i,1,k+1))
790  dc(i,1,k) = sign(min(abs(tmp),qmax,qmin), tmp)
791 1220 CONTINUE
792 
793 C
794 C ****6***0*********0*********0*********0*********0*********0**********72
795 C Loop over latitudes (to save memory)
796 C ****6***0*********0*********0*********0*********0*********0**********72
797 C
798  DO 2000 j=1,jnp
799  if((j.eq.2 .or. j.eq.jmr) .and. j1.ne.2) goto 2000
800 C
801  DO k=1,nlay
802  DO i=1,imr
803  wz2(i,k) = wz(i,j,k)
804  wk1(i,k) = p(i,j,k)
805  wk2(i,k) = delp(i,j,k)
806  flux(i,k) = dc(i,j,k) !this flux is actually the monotone slope
807  enddo
808  enddo
809 C
810 C****6***0*********0*********0*********0*********0*********0**********72
811 C Compute first guesses at cell interfaces
812 C First guesses are required to be continuous.
813 C ****6***0*********0*********0*********0*********0*********0**********72
814 C
815 C three-cell parabolic subgrid distribution at model top
816 C two-cell parabolic with zero gradient subgrid distribution
817 C at the surface.
818 C
819 C First guess top edge value
820  DO 10 i=1,imr
821 C three-cell PPM
822 C Compute a,b, and c of q = aP**2 + bP + c using cell averages and delp
823  a = 3.*( dqdt(i,j,2) - dqdt(i,j,1)*(wk2(i,2)+wk2(i,3))/
824  & (wk2(i,1)+wk2(i,2)) ) /
825  & ( (wk2(i,2)+wk2(i,3))*(wk2(i,1)+wk2(i,2)+wk2(i,3)) )
826  b = 2.*dqdt(i,j,1)/(wk2(i,1)+wk2(i,2)) -
827  & r23*a*(2.*wk2(i,1)+wk2(i,2))
828  al(i,1) = wk1(i,1) - wk2(i,1)*(r3*a*wk2(i,1) + 0.5*b)
829  al(i,2) = wk2(i,1)*(a*wk2(i,1) + b) + al(i,1)
830 C
831 C Check if change sign
832  if(wk1(i,1)*al(i,1).le.0.) then
833  al(i,1) = 0.
834  flux(i,1) = 0.
835  else
836  flux(i,1) = wk1(i,1) - al(i,1)
837  endif
838 10 continue
839 C
840 C Bottom
841  DO 15 i=1,imr
842 C 2-cell PPM with zero gradient right at the surface
843 C
844  fct = dqdt(i,j,nlaym1)*wk2(i,nlay)**2 /
845  & ( (wk2(i,nlay)+wk2(i,nlaym1))*(2.*wk2(i,nlay)+wk2(i,nlaym1)))
846  ar(i,nlay) = wk1(i,nlay) + fct
847  al(i,nlay) = wk1(i,nlay) - (fct+fct)
848  if(wk1(i,nlay)*ar(i,nlay).le.0.) ar(i,nlay) = 0.
849  flux(i,nlay) = ar(i,nlay) - wk1(i,nlay)
850 15 continue
851 
852 C
853 C****6***0*********0*********0*********0*********0*********0**********72
854 C 4th order interpolation in the interior.
855 C****6***0*********0*********0*********0*********0*********0**********72
856 C
857  DO 14 k=3,nlaym1
858  DO 12 i=1,imr
859  c1 = dqdt(i,j,k-1)*wk2(i,k-1) / (wk2(i,k-1)+wk2(i,k))
860  c2 = 2. / (wk2(i,k-2)+wk2(i,k-1)+wk2(i,k)+wk2(i,k+1))
861  a1 = (wk2(i,k-2)+wk2(i,k-1)) / (2.*wk2(i,k-1)+wk2(i,k))
862  a2 = (wk2(i,k )+wk2(i,k+1)) / (2.*wk2(i,k)+wk2(i,k-1))
863  al(i,k) = wk1(i,k-1) + c1 + c2 *
864  & ( wk2(i,k )*(c1*(a1 - a2)+a2*flux(i,k-1)) -
865  & wk2(i,k-1)*a1*flux(i,k) )
866 C print *,'AL1',i,k, AL(i,k)
867 12 CONTINUE
868 14 continue
869 C
870  do 20 i=1,imr*nlaym1
871  ar(i,1) = al(i,2)
872 C print *,'AR1',i,AR(i,1)
873 20 continue
874 C
875  do 30 i=1,imr*nlay
876  a6(i,1) = 3.*(wk1(i,1)+wk1(i,1) - (al(i,1)+ar(i,1)))
877 C print *,'A61',i,A6(i,1)
878 30 continue
879 C
880 C****6***0*********0*********0*********0*********0*********0**********72
881 C Top & Bot always monotonic
882  call lmtppm(flux(1,1),a6(1,1),ar(1,1),al(1,1),wk1(1,1),imr,0)
883  call lmtppm(flux(1,nlay),a6(1,nlay),ar(1,nlay),al(1,nlay),
884  & wk1(1,nlay),imr,0)
885 C
886 C Interior depending on KORD
887  if(lmt.LE.2)
888  & call lmtppm(flux(1,2),a6(1,2),ar(1,2),al(1,2),wk1(1,2),
889  & imr*(nlay-2),lmt)
890 C
891 C****6***0*********0*********0*********0*********0*********0**********72
892 C
893  DO 140 i=1,imr*nlaym1
894  IF(wz2(i,1).GT.0.) then
895  cm = wz2(i,1) / wk2(i,1)
896  flux(i,2) = ar(i,1)+0.5*cm*(al(i,1)-ar(i,1)+a6(i,1)*(1.-r23*cm))
897  else
898 C print *,'test2-0',i,j,wz2(i,1),wk2(i,2)
899  cp= wz2(i,1) / wk2(i,2)
900 C print *,'testCP',CP
901  flux(i,2) = al(i,2)+0.5*cp*(al(i,2)-ar(i,2)-a6(i,2)*(1.+r23*cp))
902 C print *,'test2',i, AL(i,2),AR(i,2),A6(i,2),R23
903  endif
904 140 continue
905 C
906  DO 250 i=1,imr*nlaym1
907  flux(i,2) = wz2(i,1) * flux(i,2)
908 250 continue
909 C
910  do 350 i=1,imr
911  dq(i,j, 1) = dq(i,j, 1) - flux(i, 2)
912  dq(i,j,nlay) = dq(i,j,nlay) + flux(i,nlay)
913 350 continue
914 C
915  do 360 k=2,nlaym1
916  do 360 i=1,imr
917 360 dq(i,j,k) = dq(i,j,k) + flux(i,k) - flux(i,k+1)
918 2000 continue
919  return
920  end
921 C
922  subroutine xtp(IMR,JNP,IML,j1,j2,JN,JS,PU,DQ,Q,UC,
923  & fx1,xmass,iord)
924  dimension uc(imr,*),dc(-iml:imr+iml+1),xmass(imr,jnp)
925  & ,fx1(imr+1),dq(imr,jnp),qtmp(-iml:imr+1+iml)
926  dimension pu(imr,jnp),q(imr,jnp),isave(imr)
927 C
928  imp = imr + 1
929 C
930 C van Leer at high latitudes
931  jvan = max(1,jnp/18)
932  j1vl = j1+jvan
933  j2vl = j2-jvan
934 C
935  do 1310 j=j1,j2
936 C
937  do i=1,imr
938  qtmp(i) = q(i,j)
939  enddo
940 C
941  if(j.ge.jn .or. j.le.js) goto 2222
942 C ************* Eulerian **********
943 C
944  qtmp(0) = q(imr,j)
945  qtmp(-1) = q(imr-1,j)
946  qtmp(imp) = q(1,j)
947  qtmp(imp+1) = q(2,j)
948 C
949  IF(iord.eq.1 .or. j.eq.j1. or. j.eq.j2) THEN
950  DO 1406 i=1,imr
951  iu = REAL(i) - uc(i,j)
952 1406 fx1(i) = qtmp(iu)
953  ELSE
954  call xmist(imr,iml,qtmp,dc)
955  dc(0) = dc(imr)
956 C
957  if(iord.eq.2 .or. j.le.j1vl .or. j.ge.j2vl) then
958  DO 1408 i=1,imr
959  iu = REAL(i) - uc(i,j)
960 1408 fx1(i) = qtmp(iu) + dc(iu)*(sign(1.,uc(i,j))-uc(i,j))
961  else
962  call fxppm(imr,iml,uc(1,j),qtmp,dc,fx1,iord)
963  endif
964 C
965  ENDIF
966 C
967  DO 1506 i=1,imr
968 1506 fx1(i) = fx1(i)*xmass(i,j)
969 C
970  goto 1309
971 C
972 C ***** Conservative (flux-form) Semi-Lagrangian transport *****
973 C
974 2222 continue
975 C
976  do i=-iml,0
977  qtmp(i) = q(imr+i,j)
978  qtmp(imp-i) = q(1-i,j)
979  enddo
980 C
981  IF(iord.eq.1 .or. j.eq.j1. or. j.eq.j2) THEN
982  DO 1306 i=1,imr
983  itmp = int(uc(i,j))
984  isave(i) = i - itmp
985  iu = i - uc(i,j)
986 1306 fx1(i) = (uc(i,j) - itmp)*qtmp(iu)
987  ELSE
988  call xmist(imr,iml,qtmp,dc)
989 C
990  do i=-iml,0
991  dc(i) = dc(imr+i)
992  dc(imp-i) = dc(1-i)
993  enddo
994 C
995  DO 1307 i=1,imr
996  itmp = int(uc(i,j))
997  rut = uc(i,j) - itmp
998  isave(i) = i - itmp
999  iu = i - uc(i,j)
1000 1307 fx1(i) = rut*(qtmp(iu) + dc(iu)*(sign(1.,rut) - rut))
1001  ENDIF
1002 C
1003  do 1308 i=1,imr
1004  IF(uc(i,j).GT.1.) then
1005 CDIR$ NOVECTOR
1006  do ist = isave(i),i-1
1007  fx1(i) = fx1(i) + qtmp(ist)
1008  enddo
1009  elseIF(uc(i,j).LT.-1.) then
1010  do ist = i,isave(i)-1
1011  fx1(i) = fx1(i) - qtmp(ist)
1012  enddo
1013 CDIR$ VECTOR
1014  endif
1015 1308 continue
1016  do i=1,imr
1017  fx1(i) = pu(i,j)*fx1(i)
1018  enddo
1019 C
1020 C ***************************************
1021 C
1022 1309 fx1(imp) = fx1(1)
1023  DO 1215 i=1,imr
1024 1215 dq(i,j) = dq(i,j) + fx1(i)-fx1(i+1)
1025 C
1026 C ***************************************
1027 C
1028 1310 continue
1029  return
1030  end
1031 C
1032  subroutine fxppm(IMR,IML,UT,P,DC,flux,IORD)
1033  parameter( r3 = 1./3., r23 = 2./3. )
1034  dimension ut(*),flux(*),p(-iml:imr+iml+1),dc(-iml:imr+iml+1)
1035  dimension ar(0:imr),al(0:imr),a6(0:imr)
1036  integer lmt
1037 c logical first
1038 c data first /.true./
1039 c SAVE LMT
1040 c if(first) then
1041 C
1042 C correction calcul de LMT a chaque passage pour pouvoir choisir
1043 c plusieurs schemas PPM pour differents traceurs
1044 c IF (IORD.LE.0) then
1045 c if(IMR.GE.144) then
1046 c LMT = 0
1047 c elseif(IMR.GE.72) then
1048 c LMT = 1
1049 c else
1050 c LMT = 2
1051 c endif
1052 c else
1053 c LMT = IORD - 3
1054 c endif
1055 C
1056  lmt = iord - 3
1057 c write(6,*) 'PPM option in E-W direction = ', LMT
1058 c first = .false.
1059 C endif
1060 C
1061  DO 10 i=1,imr
1062 10 al(i) = 0.5*(p(i-1)+p(i)) + (dc(i-1) - dc(i))*r3
1063 C
1064  do 20 i=1,imr-1
1065 20 ar(i) = al(i+1)
1066  ar(imr) = al(1)
1067 C
1068  do 30 i=1,imr
1069 30 a6(i) = 3.*(p(i)+p(i) - (al(i)+ar(i)))
1070 C
1071  if(lmt.LE.2) call lmtppm(dc(1),a6(1),ar(1),al(1),p(1),imr,lmt)
1072 C
1073  al(0) = al(imr)
1074  ar(0) = ar(imr)
1075  a6(0) = a6(imr)
1076 C
1077  DO i=1,imr
1078  IF(ut(i).GT.0.) then
1079  flux(i) = ar(i-1) + 0.5*ut(i)*(al(i-1) - ar(i-1) +
1080  & a6(i-1)*(1.-r23*ut(i)) )
1081  else
1082  flux(i) = al(i) - 0.5*ut(i)*(ar(i) - al(i) +
1083  & a6(i)*(1.+r23*ut(i)))
1084  endif
1085  enddo
1086  return
1087  end
1088 C
1089  subroutine xmist(IMR,IML,P,DC)
1090  parameter( r24 = 1./24.)
1091  dimension p(-iml:imr+1+iml),dc(-iml:imr+1+iml)
1092 C
1093  do 10 i=1,imr
1094  tmp = r24*(8.*(p(i+1) - p(i-1)) + p(i-2) - p(i+2))
1095  pmax = max(p(i-1), p(i), p(i+1)) - p(i)
1096  pmin = p(i) - min(p(i-1), p(i), p(i+1))
1097 10 dc(i) = sign(min(abs(tmp),pmax,pmin), tmp)
1098  return
1099  end
1100 C
1101  subroutine ytp(IMR,JNP,j1,j2,acosp,RCAP,DQ,P,VC,DC2
1102  & ,ymass,fx,a6,ar,al,jord)
1103  dimension p(imr,jnp),vc(imr,jnp),ymass(imr,jnp)
1104  & ,dc2(imr,jnp),dq(imr,jnp),acosp(jnp)
1105 C Work array
1106  dimension fx(imr,jnp),ar(imr,jnp),al(imr,jnp),a6(imr,jnp)
1107 C
1108  jmr = jnp - 1
1109  len = imr*(j2-j1+2)
1110 C
1111  if(jord.eq.1) then
1112  DO 1000 i=1,len
1113  jt = REAL(J1) - vc(i,j1)
1114 1000 fx(i,j1) = p(i,jt)
1115  else
1116 
1117  call ymist(imr,jnp,j1,p,dc2,4)
1118 C
1119  if(jord.LE.0 .or. jord.GE.3) then
1120 
1121  call fyppm(vc,p,dc2,fx,imr,jnp,j1,j2,a6,ar,al,jord)
1122 
1123  else
1124  DO 1200 i=1,len
1125  jt = REAL(J1) - vc(i,j1)
1126 1200 fx(i,j1) = p(i,jt) + (sign(1.,vc(i,j1))-vc(i,j1))*dc2(i,jt)
1127  endif
1128  endif
1129 C
1130  DO 1300 i=1,len
1131 1300 fx(i,j1) = fx(i,j1)*ymass(i,j1)
1132 C
1133  DO 1400 j=j1,j2
1134  DO 1400 i=1,imr
1135 1400 dq(i,j) = dq(i,j) + (fx(i,j) - fx(i,j+1)) * acosp(j)
1136 C
1137 C Poles
1138  sum1 = fx(imr,j1 )
1139  sum2 = fx(imr,j2+1)
1140  do i=1,imr-1
1141  sum1 = sum1 + fx(i,j1 )
1142  sum2 = sum2 + fx(i,j2+1)
1143  enddo
1144 C
1145  sum1 = dq(1, 1) - sum1 * rcap
1146  sum2 = dq(1,jnp) + sum2 * rcap
1147  do i=1,imr
1148  dq(i, 1) = sum1
1149  dq(i,jnp) = sum2
1150  enddo
1151 C
1152  if(j1.ne.2) then
1153  do i=1,imr
1154  dq(i, 2) = sum1
1155  dq(i,jmr) = sum2
1156  enddo
1157  endif
1158 C
1159  return
1160  end
1161 C
1162  subroutine ymist(IMR,JNP,j1,P,DC,ID)
1163  parameter( r24 = 1./24. )
1164  dimension p(imr,jnp),dc(imr,jnp)
1165 C
1166  imh = imr / 2
1167  jmr = jnp - 1
1168  ijm3 = imr*(jmr-3)
1169 C
1170  IF(id.EQ.2) THEN
1171  do 10 i=1,imr*(jmr-1)
1172  tmp = 0.25*(p(i,3) - p(i,1))
1173  pmax = max(p(i,1),p(i,2),p(i,3)) - p(i,2)
1174  pmin = p(i,2) - min(p(i,1),p(i,2),p(i,3))
1175  dc(i,2) = sign(min(abs(tmp),pmin,pmax),tmp)
1176 10 CONTINUE
1177  ELSE
1178  do 12 i=1,imh
1179 C J=2
1180  tmp = (8.*(p(i,3) - p(i,1)) + p(i+imh,2) - p(i,4))*r24
1181  pmax = max(p(i,1),p(i,2),p(i,3)) - p(i,2)
1182  pmin = p(i,2) - min(p(i,1),p(i,2),p(i,3))
1183  dc(i,2) = sign(min(abs(tmp),pmin,pmax),tmp)
1184 C J=JMR
1185  tmp=(8.*(p(i,jnp)-p(i,jmr-1))+p(i,jmr-2)-p(i+imh,jmr))*r24
1186  pmax = max(p(i,jmr-1),p(i,jmr),p(i,jnp)) - p(i,jmr)
1187  pmin = p(i,jmr) - min(p(i,jmr-1),p(i,jmr),p(i,jnp))
1188  dc(i,jmr) = sign(min(abs(tmp),pmin,pmax),tmp)
1189 12 CONTINUE
1190  do 14 i=imh+1,imr
1191 C J=2
1192  tmp = (8.*(p(i,3) - p(i,1)) + p(i-imh,2) - p(i,4))*r24
1193  pmax = max(p(i,1),p(i,2),p(i,3)) - p(i,2)
1194  pmin = p(i,2) - min(p(i,1),p(i,2),p(i,3))
1195  dc(i,2) = sign(min(abs(tmp),pmin,pmax),tmp)
1196 C J=JMR
1197  tmp=(8.*(p(i,jnp)-p(i,jmr-1))+p(i,jmr-2)-p(i-imh,jmr))*r24
1198  pmax = max(p(i,jmr-1),p(i,jmr),p(i,jnp)) - p(i,jmr)
1199  pmin = p(i,jmr) - min(p(i,jmr-1),p(i,jmr),p(i,jnp))
1200  dc(i,jmr) = sign(min(abs(tmp),pmin,pmax),tmp)
1201 14 CONTINUE
1202 C
1203  do 15 i=1,ijm3
1204  tmp = (8.*(p(i,4) - p(i,2)) + p(i,1) - p(i,5))*r24
1205  pmax = max(p(i,2),p(i,3),p(i,4)) - p(i,3)
1206  pmin = p(i,3) - min(p(i,2),p(i,3),p(i,4))
1207  dc(i,3) = sign(min(abs(tmp),pmin,pmax),tmp)
1208 15 CONTINUE
1209  ENDIF
1210 C
1211  if(j1.ne.2) then
1212  do i=1,imr
1213  dc(i,1) = 0.
1214  dc(i,jnp) = 0.
1215  enddo
1216  else
1217 C Determine slopes in polar caps for scalars!
1218 C
1219  do 13 i=1,imh
1220 C South
1221  tmp = 0.25*(p(i,2) - p(i+imh,2))
1222  pmax = max(p(i,2),p(i,1), p(i+imh,2)) - p(i,1)
1223  pmin = p(i,1) - min(p(i,2),p(i,1), p(i+imh,2))
1224  dc(i,1)=sign(min(abs(tmp),pmax,pmin),tmp)
1225 C North.
1226  tmp = 0.25*(p(i+imh,jmr) - p(i,jmr))
1227  pmax = max(p(i+imh,jmr),p(i,jnp), p(i,jmr)) - p(i,jnp)
1228  pmin = p(i,jnp) - min(p(i+imh,jmr),p(i,jnp), p(i,jmr))
1229  dc(i,jnp) = sign(min(abs(tmp),pmax,pmin),tmp)
1230 13 continue
1231 C
1232  do 25 i=imh+1,imr
1233  dc(i, 1) = - dc(i-imh, 1)
1234  dc(i,jnp) = - dc(i-imh,jnp)
1235 25 continue
1236  endif
1237  return
1238  end
1239 C
1240  subroutine fyppm(VC,P,DC,flux,IMR,JNP,j1,j2,A6,AR,AL,JORD)
1241  parameter( r3 = 1./3., r23 = 2./3. )
1242  real vc(imr,*),flux(imr,*),p(imr,*),dc(imr,*)
1243 C Local work arrays.
1244  real ar(imr,jnp),al(imr,jnp),a6(imr,jnp)
1245  integer lmt
1246 c logical first
1247 C data first /.true./
1248 C SAVE LMT
1249 C
1250  imh = imr / 2
1251  jmr = jnp - 1
1252  j11 = j1-1
1253  imjm1 = imr*(j2-j1+2)
1254  len = imr*(j2-j1+3)
1255 C if(first) then
1256 C IF(JORD.LE.0) then
1257 C if(JMR.GE.90) then
1258 C LMT = 0
1259 C elseif(JMR.GE.45) then
1260 C LMT = 1
1261 C else
1262 C LMT = 2
1263 C endif
1264 C else
1265 C LMT = JORD - 3
1266 C endif
1267 C
1268 C first = .false.
1269 C endif
1270 C
1271 c modifs pour pouvoir choisir plusieurs schemas PPM
1272  lmt = jord - 3
1273 C
1274  DO 10 i=1,imr*jmr
1275  al(i,2) = 0.5*(p(i,1)+p(i,2)) + (dc(i,1) - dc(i,2))*r3
1276  ar(i,1) = al(i,2)
1277 10 CONTINUE
1278 C
1279 CPoles:
1280 C
1281  DO i=1,imh
1282  al(i,1) = al(i+imh,2)
1283  al(i+imh,1) = al(i,2)
1284 C
1285  ar(i,jnp) = ar(i+imh,jmr)
1286  ar(i+imh,jnp) = ar(i,jmr)
1287  ENDDO
1288 
1289 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1290 c Rajout pour LMDZ.3.3
1291 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1292  ar(imr,1)=al(1,1)
1293  ar(imr,jnp)=al(1,jnp)
1294 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1295 
1296 
1297  do 30 i=1,len
1298 30 a6(i,j11) = 3.*(p(i,j11)+p(i,j11) - (al(i,j11)+ar(i,j11)))
1299 C
1300  if(lmt.le.2) call lmtppm(dc(1,j11),a6(1,j11),ar(1,j11)
1301  & ,al(1,j11),p(1,j11),len,lmt)
1302 C
1303 
1304  DO 140 i=1,imjm1
1305  IF(vc(i,j1).GT.0.) then
1306  flux(i,j1) = ar(i,j11) + 0.5*vc(i,j1)*(al(i,j11) - ar(i,j11) +
1307  & a6(i,j11)*(1.-r23*vc(i,j1)) )
1308  else
1309  flux(i,j1) = al(i,j1) - 0.5*vc(i,j1)*(ar(i,j1) - al(i,j1) +
1310  & a6(i,j1)*(1.+r23*vc(i,j1)))
1311  endif
1312 140 continue
1313  return
1314  end
1315 C
1316  subroutine yadv(IMR,JNP,j1,j2,p,VA,ady,wk,IAD)
1317  REAL p(imr,jnp),ady(imr,jnp),va(imr,jnp)
1318  REAL wk(imr,-1:jnp+2)
1319 C
1320  jmr = jnp-1
1321  imh = imr/2
1322  do j=1,jnp
1323  do i=1,imr
1324  wk(i,j) = p(i,j)
1325  enddo
1326  enddo
1327 C Poles:
1328  do i=1,imh
1329  wk(i, -1) = p(i+imh,3)
1330  wk(i+imh,-1) = p(i,3)
1331  wk(i, 0) = p(i+imh,2)
1332  wk(i+imh,0) = p(i,2)
1333  wk(i,jnp+1) = p(i+imh,jmr)
1334  wk(i+imh,jnp+1) = p(i,jmr)
1335  wk(i,jnp+2) = p(i+imh,jnp-2)
1336  wk(i+imh,jnp+2) = p(i,jnp-2)
1337  enddo
1338 c write(*,*) 'toto 1'
1339 C --------------------------------
1340  IF(iad.eq.2) then
1341  do j=j1-1,j2+1
1342  do i=1,imr
1343 c write(*,*) 'avt NINT','i=',i,'j=',j
1344  jp = nint(va(i,j))
1345  rv = jp - va(i,j)
1346 c write(*,*) 'VA=',VA(i,j), 'JP1=',JP,'rv=',rv
1347  jp = j - jp
1348 c write(*,*) 'JP2=',JP
1349  a1 = 0.5*(wk(i,jp+1)+wk(i,jp-1)) - wk(i,jp)
1350  b1 = 0.5*(wk(i,jp+1)-wk(i,jp-1))
1351 c write(*,*) 'a1=',a1,'b1=',b1
1352  ady(i,j) = wk(i,jp) + rv*(a1*rv + b1) - wk(i,j)
1353  enddo
1354  enddo
1355 c write(*,*) 'toto 2'
1356 C
1357  ELSEIF(iad.eq.1) then
1358  do j=j1-1,j2+1
1359  do i=1,imr
1360  jp = REAL(j)-va(i,j)
1361  ady(i,j) = va(i,j)*(wk(i,jp)-wk(i,jp+1))
1362  enddo
1363  enddo
1364  ENDIF
1365 C
1366  if(j1.ne.2) then
1367  sum1 = 0.
1368  sum2 = 0.
1369  do i=1,imr
1370  sum1 = sum1 + ady(i,2)
1371  sum2 = sum2 + ady(i,jmr)
1372  enddo
1373  sum1 = sum1 / imr
1374  sum2 = sum2 / imr
1375 C
1376  do i=1,imr
1377  ady(i, 2) = sum1
1378  ady(i,jmr) = sum2
1379  ady(i, 1) = sum1
1380  ady(i,jnp) = sum2
1381  enddo
1382  else
1383 C Poles:
1384  sum1 = 0.
1385  sum2 = 0.
1386  do i=1,imr
1387  sum1 = sum1 + ady(i,1)
1388  sum2 = sum2 + ady(i,jnp)
1389  enddo
1390  sum1 = sum1 / imr
1391  sum2 = sum2 / imr
1392 C
1393  do i=1,imr
1394  ady(i, 1) = sum1
1395  ady(i,jnp) = sum2
1396  enddo
1397  endif
1398 C
1399  return
1400  end
1401 C
1402  subroutine xadv(IMR,JNP,j1,j2,p,UA,JS,JN,IML,adx,IAD)
1403  REAL p(imr,jnp),adx(imr,jnp),qtmp(-imr:imr+imr),ua(imr,jnp)
1404 C
1405  jmr = jnp-1
1406  do 1309 j=j1,j2
1407  if(j.GT.js .and. j.LT.jn) go to 1309
1408 C
1409  do i=1,imr
1410  qtmp(i) = p(i,j)
1411  enddo
1412 C
1413  do i=-iml,0
1414  qtmp(i) = p(imr+i,j)
1415  qtmp(imr+1-i) = p(1-i,j)
1416  enddo
1417 C
1418  IF(iad.eq.2) THEN
1419  DO i=1,imr
1420  ip = nint(ua(i,j))
1421  ru = ip - ua(i,j)
1422  ip = i - ip
1423  a1 = 0.5*(qtmp(ip+1)+qtmp(ip-1)) - qtmp(ip)
1424  b1 = 0.5*(qtmp(ip+1)-qtmp(ip-1))
1425  adx(i,j) = qtmp(ip) + ru*(a1*ru + b1)
1426  enddo
1427  ELSEIF(iad.eq.1) then
1428  DO i=1,imr
1429  iu = ua(i,j)
1430  ru = ua(i,j) - iu
1431  iiu = i-iu
1432  if(ua(i,j).GE.0.) then
1433  adx(i,j) = qtmp(iiu)+ru*(qtmp(iiu-1)-qtmp(iiu))
1434  else
1435  adx(i,j) = qtmp(iiu)+ru*(qtmp(iiu)-qtmp(iiu+1))
1436  endif
1437  enddo
1438  ENDIF
1439 C
1440  do i=1,imr
1441  adx(i,j) = adx(i,j) - p(i,j)
1442  enddo
1443 1309 continue
1444 C
1445 C Eulerian upwind
1446 C
1447  do j=js+1,jn-1
1448 C
1449  do i=1,imr
1450  qtmp(i) = p(i,j)
1451  enddo
1452 C
1453  qtmp(0) = p(imr,j)
1454  qtmp(imr+1) = p(1,j)
1455 C
1456  IF(iad.eq.2) THEN
1457  qtmp(-1) = p(imr-1,j)
1458  qtmp(imr+2) = p(2,j)
1459  do i=1,imr
1460  ip = nint(ua(i,j))
1461  ru = ip - ua(i,j)
1462  ip = i - ip
1463  a1 = 0.5*(qtmp(ip+1)+qtmp(ip-1)) - qtmp(ip)
1464  b1 = 0.5*(qtmp(ip+1)-qtmp(ip-1))
1465  adx(i,j) = qtmp(ip)- p(i,j) + ru*(a1*ru + b1)
1466  enddo
1467  ELSEIF(iad.eq.1) then
1468 C 1st order
1469  DO i=1,imr
1470  ip = i - ua(i,j)
1471  adx(i,j) = ua(i,j)*(qtmp(ip)-qtmp(ip+1))
1472  enddo
1473  ENDIF
1474  enddo
1475 C
1476  if(j1.ne.2) then
1477  do i=1,imr
1478  adx(i, 2) = 0.
1479  adx(i,jmr) = 0.
1480  enddo
1481  endif
1482 C set cross term due to x-adv at the poles to zero.
1483  do i=1,imr
1484  adx(i, 1) = 0.
1485  adx(i,jnp) = 0.
1486  enddo
1487  return
1488  end
1489 C
1490  subroutine lmtppm(DC,A6,AR,AL,P,IM,LMT)
1491 C
1492 C A6 = CURVATURE OF THE TEST PARABOLA
1493 C AR = RIGHT EDGE VALUE OF THE TEST PARABOLA
1494 C AL = LEFT EDGE VALUE OF THE TEST PARABOLA
1495 C DC = 0.5 * MISMATCH
1496 C P = CELL-AVERAGED VALUE
1497 C IM = VECTOR LENGTH
1498 C
1499 C OPTIONS:
1500 C
1501 C LMT = 0: FULL MONOTONICITY
1502 C LMT = 1: SEMI-MONOTONIC CONSTRAINT (NO UNDERSHOOTS)
1503 C LMT = 2: POSITIVE-DEFINITE CONSTRAINT
1504 C
1505  parameter( r12 = 1./12. )
1506  dimension a6(im),ar(im),al(im),p(im),dc(im)
1507 C
1508  if(lmt.eq.0) then
1509 C Full constraint
1510  do 100 i=1,im
1511  if(dc(i).eq.0.) then
1512  ar(i) = p(i)
1513  al(i) = p(i)
1514  a6(i) = 0.
1515  else
1516  da1 = ar(i) - al(i)
1517  da2 = da1**2
1518  a6da = a6(i)*da1
1519  if(a6da .lt. -da2) then
1520  a6(i) = 3.*(al(i)-p(i))
1521  ar(i) = al(i) - a6(i)
1522  elseif(a6da .gt. da2) then
1523  a6(i) = 3.*(ar(i)-p(i))
1524  al(i) = ar(i) - a6(i)
1525  endif
1526  endif
1527 100 continue
1528  elseif(lmt.eq.1) then
1529 C Semi-monotonic constraint
1530  do 150 i=1,im
1531  if(abs(ar(i)-al(i)) .GE. -a6(i)) go to 150
1532  if(p(i).lt.ar(i) .and. p(i).lt.al(i)) then
1533  ar(i) = p(i)
1534  al(i) = p(i)
1535  a6(i) = 0.
1536  elseif(ar(i) .gt. al(i)) then
1537  a6(i) = 3.*(al(i)-p(i))
1538  ar(i) = al(i) - a6(i)
1539  else
1540  a6(i) = 3.*(ar(i)-p(i))
1541  al(i) = ar(i) - a6(i)
1542  endif
1543 150 continue
1544  elseif(lmt.eq.2) then
1545  do 250 i=1,im
1546  if(abs(ar(i)-al(i)) .GE. -a6(i)) go to 250
1547  fmin = p(i) + 0.25*(ar(i)-al(i))**2/a6(i) + a6(i)*r12
1548  if(fmin.ge.0.) go to 250
1549  if(p(i).lt.ar(i) .and. p(i).lt.al(i)) then
1550  ar(i) = p(i)
1551  al(i) = p(i)
1552  a6(i) = 0.
1553  elseif(ar(i) .gt. al(i)) then
1554  a6(i) = 3.*(al(i)-p(i))
1555  ar(i) = al(i) - a6(i)
1556  else
1557  a6(i) = 3.*(ar(i)-p(i))
1558  al(i) = ar(i) - a6(i)
1559  endif
1560 250 continue
1561  endif
1562  return
1563  end
1564 C
1565  subroutine a2c(U,V,IMR,JMR,j1,j2,CRX,CRY,dtdx5,DTDY5)
1566  dimension u(imr,*),v(imr,*),crx(imr,*),cry(imr,*),dtdx5(*)
1567 C
1568  do 35 j=j1,j2
1569  do 35 i=2,imr
1570 35 crx(i,j) = dtdx5(j)*(u(i,j)+u(i-1,j))
1571 C
1572  do 45 j=j1,j2
1573 45 crx(1,j) = dtdx5(j)*(u(1,j)+u(imr,j))
1574 C
1575  do 55 i=1,imr*jmr
1576 55 cry(i,2) = dtdy5*(v(i,2)+v(i,1))
1577  return
1578  end
1579 C
1580  subroutine cosa(cosp,cose,JNP,PI,DP)
1581  dimension cosp(*),cose(*)
1582  jmr = jnp-1
1583  do 55 j=2,jnp
1584  ph5 = -0.5*pi + (REAL(j-1)-0.5)*dp
1585 55 cose(j) = cos(ph5)
1586 C
1587  jeq = (jnp+1) / 2
1588  if(jmr .eq. 2*(jmr/2) ) then
1589  do j=jnp, jeq+1, -1
1590  cose(j) = cose(jnp+2-j)
1591  enddo
1592  else
1593 C cell edge at equator.
1594  cose(jeq+1) = 1.
1595  do j=jnp, jeq+2, -1
1596  cose(j) = cose(jnp+2-j)
1597  enddo
1598  endif
1599 C
1600  do 66 j=2,jmr
1601 66 cosp(j) = 0.5*(cose(j)+cose(j+1))
1602  cosp(1) = 0.
1603  cosp(jnp) = 0.
1604  return
1605  end
1606 C
1607  subroutine cosc(cosp,cose,JNP,PI,DP)
1608  dimension cosp(*),cose(*)
1609 C
1610  phi = -0.5*pi
1611  do 55 j=2,jnp-1
1612  phi = phi + dp
1613 55 cosp(j) = cos(phi)
1614  cosp( 1) = 0.
1615  cosp(jnp) = 0.
1616 C
1617  do 66 j=2,jnp
1618  cose(j) = 0.5*(cosp(j)+cosp(j-1))
1619 66 CONTINUE
1620 C
1621  do 77 j=2,jnp-1
1622  cosp(j) = 0.5*(cose(j)+cose(j+1))
1623 77 CONTINUE
1624  return
1625  end
1626 C
1627  SUBROUTINE qckxyz (Q,qtmp,IMR,JNP,NLAY,j1,j2,cosp,acosp,
1628  & cross,ic,nstep)
1629 C
1630  parameter( tiny = 1.e-60 )
1631  dimension q(imr,jnp,nlay),qtmp(imr,jnp),cosp(*),acosp(*)
1632  logical cross
1633 C
1634  nlaym1 = nlay-1
1635  len = imr*(j2-j1+1)
1636  ip = 0
1637 C
1638 C Top layer
1639  l = 1
1640  icr = 1
1641  call filns(q(1,1,l),imr,jnp,j1,j2,cosp,acosp,ipy,tiny)
1642  if(ipy.eq.0) goto 50
1643  call filew(q(1,1,l),qtmp,imr,jnp,j1,j2,ipx,tiny)
1644  if(ipx.eq.0) goto 50
1645 C
1646  if(cross) then
1647  call filcr(q(1,1,l),imr,jnp,j1,j2,cosp,acosp,icr,tiny)
1648  endif
1649  if(icr.eq.0) goto 50
1650 C
1651 C Vertical filling...
1652  do i=1,len
1653  IF( q(i,j1,1).LT.0.) THEN
1654  ip = ip + 1
1655  q(i,j1,2) = q(i,j1,2) + q(i,j1,1)
1656  q(i,j1,1) = 0.
1657  endif
1658  enddo
1659 C
1660 50 continue
1661  DO 225 l = 2,nlaym1
1662  icr = 1
1663 C
1664  call filns(q(1,1,l),imr,jnp,j1,j2,cosp,acosp,ipy,tiny)
1665  if(ipy.eq.0) goto 225
1666  call filew(q(1,1,l),qtmp,imr,jnp,j1,j2,ipx,tiny)
1667  if(ipx.eq.0) go to 225
1668  if(cross) then
1669  call filcr(q(1,1,l),imr,jnp,j1,j2,cosp,acosp,icr,tiny)
1670  endif
1671  if(icr.eq.0) goto 225
1672 C
1673  do i=1,len
1674  IF( q(i,j1,l).LT.0.) THEN
1675 C
1676  ip = ip + 1
1677 C From above
1678  qup = q(i,j1,l-1)
1679  qly = -q(i,j1,l)
1680  dup = min(qly,qup)
1681  q(i,j1,l-1) = qup - dup
1682  q(i,j1,l ) = dup-qly
1683 C Below
1684  q(i,j1,l+1) = q(i,j1,l+1) + q(i,j1,l)
1685  q(i,j1,l) = 0.
1686  ENDIF
1687  ENDDO
1688 225 CONTINUE
1689 C
1690 C BOTTOM LAYER
1691  sum = 0.
1692  l = nlay
1693 C
1694  call filns(q(1,1,l),imr,jnp,j1,j2,cosp,acosp,ipy,tiny)
1695  if(ipy.eq.0) goto 911
1696  call filew(q(1,1,l),qtmp,imr,jnp,j1,j2,ipx,tiny)
1697  if(ipx.eq.0) goto 911
1698 C
1699  call filcr(q(1,1,l),imr,jnp,j1,j2,cosp,acosp,icr,tiny)
1700  if(icr.eq.0) goto 911
1701 C
1702  DO i=1,len
1703  IF( q(i,j1,l).LT.0.) THEN
1704  ip = ip + 1
1705 c
1706 C From above
1707 C
1708  qup = q(i,j1,nlaym1)
1709  qly = -q(i,j1,l)
1710  dup = min(qly,qup)
1711  q(i,j1,nlaym1) = qup - dup
1712 C From "below" the surface.
1713  sum = sum + qly-dup
1714  q(i,j1,l) = 0.
1715  ENDIF
1716  ENDDO
1717 C
1718 911 continue
1719 C
1720  if(ip.gt.imr) then
1721  write(6,*) 'IC=',ic,' STEP=',nstep,
1722  & ' Vertical filling pts=',ip
1723  endif
1724 C
1725  if(sum.gt.1.e-25) then
1726  write(6,*) ic,nstep,' Mass source from the ground=',sum
1727  endif
1728  RETURN
1729  END
1730 C
1731  subroutine filcr(q,IMR,JNP,j1,j2,cosp,acosp,icr,tiny)
1732  dimension q(imr,*),cosp(*),acosp(*)
1733  icr = 0
1734  do 65 j=j1+1,j2-1
1735  DO 50 i=1,imr-1
1736  IF(q(i,j).LT.0.) THEN
1737  icr = 1
1738  dq = - q(i,j)*cosp(j)
1739 C N-E
1740  dn = q(i+1,j+1)*cosp(j+1)
1741  d0 = max(0.,dn)
1742  d1 = min(dq,d0)
1743  q(i+1,j+1) = (dn - d1)*acosp(j+1)
1744  dq = dq - d1
1745 C S-E
1746  ds = q(i+1,j-1)*cosp(j-1)
1747  d0 = max(0.,ds)
1748  d2 = min(dq,d0)
1749  q(i+1,j-1) = (ds - d2)*acosp(j-1)
1750  q(i,j) = (d2 - dq)*acosp(j) + tiny
1751  endif
1752 50 continue
1753  if(icr.eq.0 .and. q(imr,j).ge.0.) goto 65
1754  DO 55 i=2,imr
1755  IF(q(i,j).LT.0.) THEN
1756  icr = 1
1757  dq = - q(i,j)*cosp(j)
1758 C N-W
1759  dn = q(i-1,j+1)*cosp(j+1)
1760  d0 = max(0.,dn)
1761  d1 = min(dq,d0)
1762  q(i-1,j+1) = (dn - d1)*acosp(j+1)
1763  dq = dq - d1
1764 C S-W
1765  ds = q(i-1,j-1)*cosp(j-1)
1766  d0 = max(0.,ds)
1767  d2 = min(dq,d0)
1768  q(i-1,j-1) = (ds - d2)*acosp(j-1)
1769  q(i,j) = (d2 - dq)*acosp(j) + tiny
1770  endif
1771 55 continue
1772 C *****************************************
1773 C i=1
1774  i=1
1775  IF(q(i,j).LT.0.) THEN
1776  icr = 1
1777  dq = - q(i,j)*cosp(j)
1778 C N-W
1779  dn = q(imr,j+1)*cosp(j+1)
1780  d0 = max(0.,dn)
1781  d1 = min(dq,d0)
1782  q(imr,j+1) = (dn - d1)*acosp(j+1)
1783  dq = dq - d1
1784 C S-W
1785  ds = q(imr,j-1)*cosp(j-1)
1786  d0 = max(0.,ds)
1787  d2 = min(dq,d0)
1788  q(imr,j-1) = (ds - d2)*acosp(j-1)
1789  q(i,j) = (d2 - dq)*acosp(j) + tiny
1790  endif
1791 C *****************************************
1792 C i=IMR
1793  i=imr
1794  IF(q(i,j).LT.0.) THEN
1795  icr = 1
1796  dq = - q(i,j)*cosp(j)
1797 C N-E
1798  dn = q(1,j+1)*cosp(j+1)
1799  d0 = max(0.,dn)
1800  d1 = min(dq,d0)
1801  q(1,j+1) = (dn - d1)*acosp(j+1)
1802  dq = dq - d1
1803 C S-E
1804  ds = q(1,j-1)*cosp(j-1)
1805  d0 = max(0.,ds)
1806  d2 = min(dq,d0)
1807  q(1,j-1) = (ds - d2)*acosp(j-1)
1808  q(i,j) = (d2 - dq)*acosp(j) + tiny
1809  endif
1810 C *****************************************
1811 65 continue
1812 C
1813  do i=1,imr
1814  if(q(i,j1).lt.0. .or. q(i,j2).lt.0.) then
1815  icr = 1
1816  goto 80
1817  endif
1818  enddo
1819 C
1820 80 continue
1821 C
1822  if(q(1,1).lt.0. .or. q(1,jnp).lt.0.) then
1823  icr = 1
1824  endif
1825 C
1826  return
1827  end
1828 C
1829  subroutine filns(q,IMR,JNP,j1,j2,cosp,acosp,ipy,tiny)
1830  dimension q(imr,*),cosp(*),acosp(*)
1831 c logical first
1832 c data first /.true./
1833 c save cap1
1834 C
1835 c if(first) then
1836  dp = 4.*atan(1.)/REAL(jnp-1)
1837  cap1 = imr*(1.-cos((j1-1.5)*dp))/dp
1838 c first = .false.
1839 c endif
1840 C
1841  ipy = 0
1842  do 55 j=j1+1,j2-1
1843  DO 55 i=1,imr
1844  IF(q(i,j).LT.0.) THEN
1845  ipy = 1
1846  dq = - q(i,j)*cosp(j)
1847 C North
1848  dn = q(i,j+1)*cosp(j+1)
1849  d0 = max(0.,dn)
1850  d1 = min(dq,d0)
1851  q(i,j+1) = (dn - d1)*acosp(j+1)
1852  dq = dq - d1
1853 C South
1854  ds = q(i,j-1)*cosp(j-1)
1855  d0 = max(0.,ds)
1856  d2 = min(dq,d0)
1857  q(i,j-1) = (ds - d2)*acosp(j-1)
1858  q(i,j) = (d2 - dq)*acosp(j) + tiny
1859  endif
1860 55 continue
1861 C
1862  do i=1,imr
1863  IF(q(i,j1).LT.0.) THEN
1864  ipy = 1
1865  dq = - q(i,j1)*cosp(j1)
1866 C North
1867  dn = q(i,j1+1)*cosp(j1+1)
1868  d0 = max(0.,dn)
1869  d1 = min(dq,d0)
1870  q(i,j1+1) = (dn - d1)*acosp(j1+1)
1871  q(i,j1) = (d1 - dq)*acosp(j1) + tiny
1872  endif
1873  enddo
1874 C
1875  j = j2
1876  do i=1,imr
1877  IF(q(i,j).LT.0.) THEN
1878  ipy = 1
1879  dq = - q(i,j)*cosp(j)
1880 C South
1881  ds = q(i,j-1)*cosp(j-1)
1882  d0 = max(0.,ds)
1883  d2 = min(dq,d0)
1884  q(i,j-1) = (ds - d2)*acosp(j-1)
1885  q(i,j) = (d2 - dq)*acosp(j) + tiny
1886  endif
1887  enddo
1888 C
1889 C Check Poles.
1890  if(q(1,1).lt.0.) then
1891  dq = q(1,1)*cap1/REAL(imr)*acosp(j1)
1892  do i=1,imr
1893  q(i,1) = 0.
1894  q(i,j1) = q(i,j1) + dq
1895  if(q(i,j1).lt.0.) ipy = 1
1896  enddo
1897  endif
1898 C
1899  if(q(1,jnp).lt.0.) then
1900  dq = q(1,jnp)*cap1/REAL(imr)*acosp(j2)
1901  do i=1,imr
1902  q(i,jnp) = 0.
1903  q(i,j2) = q(i,j2) + dq
1904  if(q(i,j2).lt.0.) ipy = 1
1905  enddo
1906  endif
1907 C
1908  return
1909  end
1910 C
1911  subroutine filew(q,qtmp,IMR,JNP,j1,j2,ipx,tiny)
1912  dimension q(imr,*),qtmp(jnp,imr)
1913 C
1914  ipx = 0
1915 C Copy & swap direction for vectorization.
1916  do 25 i=1,imr
1917  do 25 j=j1,j2
1918 25 qtmp(j,i) = q(i,j)
1919 C
1920  do 55 i=2,imr-1
1921  do 55 j=j1,j2
1922  if(qtmp(j,i).lt.0.) then
1923  ipx = 1
1924 c west
1925  d0 = max(0.,qtmp(j,i-1))
1926  d1 = min(-qtmp(j,i),d0)
1927  qtmp(j,i-1) = qtmp(j,i-1) - d1
1928  qtmp(j,i) = qtmp(j,i) + d1
1929 c east
1930  d0 = max(0.,qtmp(j,i+1))
1931  d2 = min(-qtmp(j,i),d0)
1932  qtmp(j,i+1) = qtmp(j,i+1) - d2
1933  qtmp(j,i) = qtmp(j,i) + d2 + tiny
1934  endif
1935 55 continue
1936 c
1937  i=1
1938  do 65 j=j1,j2
1939  if(qtmp(j,i).lt.0.) then
1940  ipx = 1
1941 c west
1942  d0 = max(0.,qtmp(j,imr))
1943  d1 = min(-qtmp(j,i),d0)
1944  qtmp(j,imr) = qtmp(j,imr) - d1
1945  qtmp(j,i) = qtmp(j,i) + d1
1946 c east
1947  d0 = max(0.,qtmp(j,i+1))
1948  d2 = min(-qtmp(j,i),d0)
1949  qtmp(j,i+1) = qtmp(j,i+1) - d2
1950 c
1951  qtmp(j,i) = qtmp(j,i) + d2 + tiny
1952  endif
1953 65 continue
1954  i=imr
1955  do 75 j=j1,j2
1956  if(qtmp(j,i).lt.0.) then
1957  ipx = 1
1958 c west
1959  d0 = max(0.,qtmp(j,i-1))
1960  d1 = min(-qtmp(j,i),d0)
1961  qtmp(j,i-1) = qtmp(j,i-1) - d1
1962  qtmp(j,i) = qtmp(j,i) + d1
1963 c east
1964  d0 = max(0.,qtmp(j,1))
1965  d2 = min(-qtmp(j,i),d0)
1966  qtmp(j,1) = qtmp(j,1) - d2
1967 c
1968  qtmp(j,i) = qtmp(j,i) + d2 + tiny
1969  endif
1970 75 continue
1971 C
1972  if(ipx.ne.0) then
1973  do 85 j=j1,j2
1974  do 85 i=1,imr
1975 85 q(i,j) = qtmp(j,i)
1976  else
1977 C
1978 C Poles.
1979  if(q(1,1).lt.0. or. q(1,jnp).lt.0.) ipx = 1
1980  endif
1981  return
1982  end
1983 C
1984  subroutine zflip(q,im,km,nc)
1985 C This routine flip the array q (in the vertical).
1986  real q(im,km,nc)
1987 C local dynamic array
1988  real qtmp(im,km)
1989 C
1990  do 4000 ic = 1, nc
1991 C
1992  do 1000 k=1,km
1993  do 1000 i=1,im
1994  qtmp(i,k) = q(i,km+1-k,ic)
1995 1000 continue
1996 C
1997  do 2000 i=1,im*km
1998 2000 q(i,1,ic) = qtmp(i,1)
1999 4000 continue
2000  return
2001  end