source: LMDZ6/trunk/libf/dyn3d_common/ppm3d.F @ 5069

Last change on this file since 5069 was 2197, checked in by Ehouarn Millour, 10 years ago

Added 'implicit none' statements and proper variable definitions where they were missing.
EM

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