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

Annotation of /trunk/Sources/dyn3d/PPM3d/ppm3d.f

Parent Directory Parent Directory | Revision Log Revision Log

Revision 10 - (hide annotations)
Fri Apr 18 14:45:53 2008 UTC (16 years, 1 month ago) by guez
Original Path: trunk/libf/dyn3d/ppm3d.f
File size: 53742 byte(s) 
Added NetCDF directory "/home/guez/include" in "g95.mk" and
"nag_tools.mk".

Added some "intent" attributes in "PVtheta", "advtrac", "caladvtrac",
"calfis", "diagedyn", "dissip", "vlspltqs", "aeropt", "ajsec",
"calltherm", "clmain", "cltrac", "cltracrn", "concvl", "conema3",
"conflx", "fisrtilp", "newmicro", "nuage", "diagcld1", "diagcld2",
"drag_noro", "lift_noro", "SUGWD", "physiq", "phytrac", "radlwsw", "thermcell".

Removed the case "ierr == 0" in "abort_gcm"; moved call to "histclo"
and messages for end of run from "abort_gcm" to "gcm"; replaced call
to "abort_gcm" in "leapfrog" by exit from outer loop.

In "calfis": removed argument "pp" and variable "unskap"; changed
"pksurcp" from scalar to rank 2; use "pressure_var"; rewrote
computation of "zplev", "zplay", "ztfi", "pcvgt" using "dyn_phy";
added computation of "pls".

Removed unused variable in "dynredem0".

In "exner_hyb": changed "dellta" from scalar to rank 1; replaced call
to "ssum" by call to "sum"; removed variables "xpn" and "xps";
replaced some loops by array expressions.

In "leapfrog": use "pressure_var"; deleted variables "p", "longcles".

Converted common blocks "YOECUMF", "YOEGWD" to modules.

Removed argument "pplay" in "cvltr", "diagetpq", "nflxtr".

Created module "raddimlw" from include file "raddimlw.h".

Corrected call to "new_unit" in "test_disvert".
1 guez 3 !
2     ! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/ppm3d.F,v 1.1.1.1 2004/05/19 12:53:07 lmdzadmin Exp $
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 guez 10 C NC: total number of constituents
84     C IMR: first dimension (E-W); number of Grid intervals in E-W is IMR
85     C JNP: 2nd dimension (N-S); number of Grid intervals in N-S is JNP-1
86     C NLAY: 3rd dimension (number of layers); vertical index increases from 1 at
87 guez 3 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 guez 10 C when and where Courant number is greater than one.
244 guez 3 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 / float(IMR)
348     DP = PI / float(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 = float(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 guez 10 C Convert winds on A-Grid to Courant number on C-Grid.
445 guez 3 call A2C(U(1,1,k),V(1,1,k),IMR,JMR,j1,j2,CRX,CRY,dtdx5,DTDY5)
446     else
447 guez 10 C Convert winds on C-grid to Courant number
448 guez 3 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 = float(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,JMP,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 = float(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 = float(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 = float(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 = float(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 = float(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 + (FLOAT(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.)/float(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/float(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/float(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
  ViewVC Help
Powered by ViewVC 1.1.21