source: LMDZ4/branches/LMDZ4V5.0-LF/libf/dyn3dpar/ppm3d.F @ 4778

Last change on this file since 4778 was 1299, checked in by Laurent Fairhead, 15 years ago

Nettoyage general pour se rapprocher des normes et éviter des erreurs a la
compilation:

  • tous les FLOAT() sont remplacés par des REAL()
  • tous les STOP dans phylmd sont remplacés par des appels à abort_gcm
  • le common control défini dans le fichier control.h est remplacé par le module control_mod pour éviter des messages sur l'alignement des variables dans les déclarations
  • des $Header$ remplacés par des $Id$ pour svn

Quelques remplacements à faire ont pu m'échapper

General cleanup of the code to try and adhere to norms and to prevent some
compilation errors:

  • all FLOAT() instructions have been replaced by REAL() instructions
  • all STOP instructions in phylmd have been replaced by calls to abort_gcm
  • the common block control defined in the control.h file has been replaced by the control_mod to prevent compilation warnings on the alignement of declared variables
  • $Header$ replaced by $Id$ for svn

Some changes which should have been made might have escaped me
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 52.3 KB
Line 
1!
2! $Id: ppm3d.F 1299 2010-01-20 14:27:21Z fhourdin $
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
68c      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      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      SAVE DTDY, DTDY5, RCAP, JS0, JN0, IML,
302     &     DTDX, DTDX5, ACOSP, COSP, COSE, DAP,DBK
303C
304           
305      JMR = JNP -1
306      IMJM  = IMR*JNP
307      j2 = JNP - j1 + 1
308      NSTEP = NSTEP + 1
309C
310C *********** Initialization **********************
311      if(NSTEP.eq.1) then
312c
313      write(6,*) '------------------------------------ '
314      write(6,*) 'NASA/GSFC Transport Core Version 4.5'
315      write(6,*) '------------------------------------ '
316c
317      WRITE(6,*) 'IMR=',IMR,' JNP=',JNP,' NLAY=',NLAY,' j1=',j1
318      WRITE(6,*) 'NC=',NC,IORD,JORD,KORD,NDT
319C
320C 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
335C
336      if(Jmax.lt.JNP .or. Kmax.lt.NLAY) then
337        write(6,*) 'Jmax or Kmax is too small'
338        stop
339      endif
340C
341      DO k=1,NLAY
342      DAP(k) = (AP(k+1) - AP(k))*PT
343      DBK(k) =  BP(k+1) - BP(k)
344      ENDDO     
345C
346      PI = 4. * ATAN(1.)
347      DL = 2.*PI / REAL(IMR)
348      DP =    PI / REAL(JMR)
349C
350      if(IGD.eq.0) then
351C Compute analytic cosine at cell edges
352            call cosa(cosp,cose,JNP,PI,DP)
353      else
354C Define cosine consistent with GEOS-GCM (using dycore2.0 or later)
355            call cosc(cosp,cose,JNP,PI,DP)
356      endif
357C
358      do 15 J=2,JMR
35915    acosp(j) = 1. / cosp(j)
360C
361C Inverse of the Scaled polar cap area.
362C
363      RCAP  = DP / (IMR*(1.-COS((j1-1.5)*DP)))
364      acosp(1)   = RCAP
365      acosp(JNP) = RCAP
366      endif
367C
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
381C
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)
389C
390      JS0 = REAL(JMR)*(90.-ZTC)/180. + 2
391      JS0 = max(JS0, J1+1)
392      IML = min(6*JS0/(J1-1)+2, 4*IMR/5)
393      JN0 = JNP-JS0+1
394      endif
395C     
396C
397      do J=2,JMR
398      DTDX(j)  = DT / ( DL*AE*COSP(J) )
399
400c     print*,'dtdx=',dtdx(j)
401      DTDX5(j) = 0.5*DTDX(j)
402      enddo
403C
404     
405      DTDY  = DT /(AE*DP)
406c      print*,'dtdy=',dtdy
407      DTDY5 = 0.5*DTDY
408C
409c      write(6,*) 'J1=',J1,' J2=', J2
410      endif
411C
412C *********** End Initialization **********************
413C
414C 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         
424C
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)
43040    Q(I,JMR,L,IC) = Q(I,JNP,L,IC)
431      endif
432C
433C 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
43844    DQ(i,j,k,IC) = Q(i,j,k,IC)*delp1(i,j,k)
439550     continue
440C
441      do 1500 k=1,NLAY
442C
443      if(IGD.eq.0) then
444C Convert winds on A-Grid to Courant # on C-Grid.
445      call A2C(U(1,1,k),V(1,1,k),IMR,JMR,j1,j2,CRX,CRY,dtdx5,DTDY5)
446      else
447C Convert winds on C-grid to Courant #
448      do 45 j=j1,j2
449      do 45 i=2,IMR
45045    CRX(i,J) = dtdx(j)*U(i-1,j,k)
451   
452C
453      do 50 j=j1,j2
45450    CRX(1,J) = dtdx(j)*U(IMR,j,k)
455C
456      do 55 i=1,IMR*JMR
45755    CRY(i,2) = DTDY*V(i,1,k)
458      endif
459C     
460C Determine JS and JN
461      JS = j1
462      JN = j2
463C
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
472C
4732222  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
4822233  continue
483C
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
490C
491C ******* Compute horizontal mass fluxes ************
492C
493C 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
500C
501      do 95 j=j1,j2
502      DO 95 i=1,IMR
50395    DPI(i,j,k) = (ymass(i,j) - ymass(i,j+1)) * acosp(j)
504C
505C 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
512C
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
519C
520C E-W component
521C
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
527C
528      do j=j1,j2
529      PU(1,j) = 0.5 * (delp2(1,j,k) + delp2(IMR,j,k))
530      enddo
531C
532      do 110 j=j1,j2
533      DO 110 i=1,IMR
534110   xmass(i,j) = PU(i,j)*CRX(i,j)
535C
536      DO 120 j=j1,j2
537      DO 120 i=1,IMR-1
538120   DPI(i,j,k) = DPI(i,j,k) + xmass(i,j) - xmass(i+1,j)
539C
540      DO 130 j=j1,j2
541130   DPI(IMR,j,k) = DPI(IMR,j,k) + xmass(IMR,j) - xmass(1,j)
542C
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
548C
549      DO j=j1,j2
550      UA(imr,j) = 0.5 * (CRX(imr,j)+CRX(1,j))
551      enddo
552ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
553c Rajouts pour LMDZ.3.3
554ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
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
564C
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
576C
577C ****6***0*********0*********0*********0*********0*********0**********72
578      do 1000 IC=1,NC
579C
580      do i=1,IMJM
581      wk1(i,1,1) = 0.
582      wk1(i,1,2) = 0.
583      enddo
584C
585C E-W advective cross term
586      do 250 j=J1,J2
587      if(J.GT.JS  .and. J.LT.JN) GO TO 250
588C
589      do i=1,IMR
590      qtmp(i) = q(i,j,k,IC)
591      enddo
592C
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
597C
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)
608230   continue
609250   continue
610C
611      if(JN.ne.0) then
612      do j=JS+1,JN-1
613C
614      do i=1,IMR
615      qtmp(i) = q(i,j,k,IC)
616      enddo
617C
618      qtmp(0)     = q(IMR,J,k,IC)
619      qtmp(IMR+1) = q(  1,J,k,IC)
620C
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
627C ****6***0*********0*********0*********0*********0*********0**********72
628C Contribution from the N-S advection
629      do i=1,imr*(j2-j1+1)
630      JT = REAL(J1) - VA(i,j1)
631      wk1(i,j1,2) = VA(i,j1) * (q(i,jt,k,IC) - q(i,jt+1,k,IC))
632      enddo
633C
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
638C
639        if(cross) then
640C Add cross terms in the vertical direction.
641        if(IORD .GE. 2) then
642                iad = 2
643        else
644                iad = 1
645        endif
646C
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
660C
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)
666C
6671000  continue
6681500  continue
669C
670C ******* Compute vertical mass flux (same unit as PS) ***********
671C
672C 1st step: compute total column mass CONVERGENCE.
673C
674      do 320 j=1,JNP
675      do 320 i=1,IMR
676320   CRY(i,j) = DPI(i,j,1)
677C
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)
682330   continue
683C
684      do 360 j=1,JNP
685      do 360 i=1,IMR
686C
687C 2nd step: compute PS2 (PS at n+1) using the hydrostatic assumption.
688C Changes (increases) to surface pressure = total column mass convergence
689C
690      PS2(i,j)  = PS1(i,j) + CRY(i,j)
691C
692C 3rd step: compute vertical mass flux from mass conservation principle.
693C
694      W(i,j,1) = DPI(i,j,1) - DBK(1)*CRY(i,j)
695      W(i,j,NLAY) = 0.
696360   continue
697C
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)
702370   continue
703C
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)
708380   continue
709C
710        KRD = max(3, KORD)
711      do 4000 IC=1,NC
712C
713C****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)
717C
718   
719      if(fill) call qckxyz(DQ(1,1,1,IC),DC2,IMR,JNP,NLAY,j1,j2,
720     &                     cosp,acosp,.false.,IC,NSTEP)
721C
722C Recover tracer mixing ratio from "density" using predicted
723C "air density" (pressure thickness) at time-level n+1
724C
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)
729c            print*,'i=',i,'j=',j,'k=',k,'Q(i,j,k,IC)=',Q(i,j,k,IC)
730      enddo
731      enddo
732      enddo
733C     
734      if(j1.ne.2) then
735      DO 400 k=1,NLAY
736      DO 400 I=1,IMR
737c     j=1 c'est le pôle Sud, j=JNP c'est le pôle Nord
738      Q(I,  2,k,IC) = Q(I,  1,k,IC)
739      Q(I,JMR,k,IC) = Q(I,JNP,k,IC)
740400   CONTINUE
741      endif
7424000  continue
743C
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)
7495000  continue
750      endif
751C
752      RETURN
753      END
754C
755C****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)
763C Assuming JNP >= NLAY
764      real AR(IMR,*),AL(IMR,*),A6(IMR,*),flux(IMR,*),wk2(IMR,*),
765     &     wz2(IMR,*)
766C
767      JMR = JNP - 1
768      IMJM = IMR*JNP
769      NLAYM1 = NLAY - 1
770C
771      LMT = KORD - 3
772C
773C ****6***0*********0*********0*********0*********0*********0**********72
774C Compute DC for PPM
775C ****6***0*********0*********0*********0*********0*********0**********72
776C
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)
7801000  continue
781C
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)   
7911220  CONTINUE
792     
793C     
794C ****6***0*********0*********0*********0*********0*********0**********72
795C Loop over latitudes  (to save memory)
796C ****6***0*********0*********0*********0*********0*********0**********72
797C
798      DO 2000 j=1,JNP
799      if((j.eq.2 .or. j.eq.JMR) .and. j1.ne.2) goto 2000
800C
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
809C
810C****6***0*********0*********0*********0*********0*********0**********72
811C Compute first guesses at cell interfaces
812C First guesses are required to be continuous.
813C ****6***0*********0*********0*********0*********0*********0**********72
814C
815C three-cell parabolic subgrid distribution at model top
816C two-cell parabolic with zero gradient subgrid distribution
817C at the surface.
818C
819C First guess top edge value
820      DO 10 i=1,IMR
821C three-cell PPM
822C 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)
830C
831C 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
83810    continue
839C
840C Bottom
841      DO 15 i=1,IMR
842C 2-cell PPM with zero gradient right at the surface
843C
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)
85015    continue
851     
852C
853C****6***0*********0*********0*********0*********0*********0**********72
854C 4th order interpolation in the interior.
855C****6***0*********0*********0*********0*********0*********0**********72
856C
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)  )
866C      print *,'AL1',i,k, AL(i,k)
86712    CONTINUE
86814    continue
869C
870      do 20 i=1,IMR*NLAYM1
871      AR(i,1) = AL(i,2)
872C      print *,'AR1',i,AR(i,1)
87320    continue
874C
875      do 30 i=1,IMR*NLAY
876      A6(i,1) = 3.*(wk1(i,1)+wk1(i,1) - (AL(i,1)+AR(i,1)))
877C      print *,'A61',i,A6(i,1)
87830    continue
879C
880C****6***0*********0*********0*********0*********0*********0**********72
881C 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)
885C
886C 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)
890C
891C****6***0*********0*********0*********0*********0*********0**********72
892C
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
898C        print *,'test2-0',i,j,wz2(i,1),wk2(i,2)
899        CP= wz2(i,1) / wk2(i,2)       
900C        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))
902C        print *,'test2',i, AL(i,2),AR(i,2),A6(i,2),R23
903      endif
904140   continue
905C
906      DO 250 i=1,IMR*NLAYM1
907      flux(i,2) = wz2(i,1) * flux(i,2)
908250   continue
909C
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)
913350   continue
914C
915      do 360 k=2,NLAYM1
916      do 360 i=1,IMR
917360   DQ(i,j,k) = DQ(i,j,k) + flux(i,k) - flux(i,k+1)
9182000  continue
919      return
920      end
921C
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)
927C
928      IMP = IMR + 1
929C
930C van Leer at high latitudes
931      jvan = max(1,JNP/18)
932      j1vl = j1+jvan
933      j2vl = j2-jvan
934C
935      do 1310 j=j1,j2
936C
937      do i=1,IMR
938      qtmp(i) = q(i,j)
939      enddo
940C
941      if(j.ge.JN .or. j.le.JS) goto 2222
942C ************* Eulerian **********
943C
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)
948C
949      IF(IORD.eq.1 .or. j.eq.j1. or. j.eq.j2) THEN
950      DO 1406 i=1,IMR
951      iu = REAL(i) - uc(i,j)
9521406  fx1(i) = qtmp(iu)
953      ELSE
954      call xmist(IMR,IML,Qtmp,DC)
955      DC(0) = DC(IMR)
956C
957      if(IORD.eq.2 .or. j.le.j1vl .or. j.ge.j2vl) then
958      DO 1408 i=1,IMR
959      iu = REAL(i) - uc(i,j)
9601408  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
964C
965      ENDIF
966C
967      DO 1506 i=1,IMR
9681506  fx1(i) = fx1(i)*xmass(i,j)
969C
970      goto 1309
971C
972C ***** Conservative (flux-form) Semi-Lagrangian transport *****
973C
9742222  continue
975C
976      do i=-IML,0
977      qtmp(i)     = q(IMR+i,j)
978      qtmp(IMP-i) = q(1-i,j)
979      enddo
980C
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)
9861306  fx1(i) = (uc(i,j) - itmp)*qtmp(iu)
987      ELSE
988      call xmist(IMR,IML,Qtmp,DC)
989C
990      do i=-IML,0
991      DC(i)     = DC(IMR+i)
992      DC(IMP-i) = DC(1-i)
993      enddo
994C
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)
10001307  fx1(i) = rut*(qtmp(iu) + DC(iu)*(sign(1.,rut) - rut))
1001      ENDIF
1002C
1003      do 1308 i=1,IMR
1004      IF(uc(i,j).GT.1.) then
1005CDIR$ 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
1013CDIR$ VECTOR
1014      endif
10151308  continue
1016      do i=1,IMR
1017      fx1(i) = PU(i,j)*fx1(i)
1018      enddo
1019C
1020C ***************************************
1021C
10221309  fx1(IMP) = fx1(1)
1023      DO 1215 i=1,IMR
10241215  DQ(i,j) =  DQ(i,j) + fx1(i)-fx1(i+1)
1025C
1026C ***************************************
1027C
10281310  continue
1029      return
1030      end
1031C
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
1037c      logical first
1038c      data first /.true./
1039c      SAVE LMT
1040c      if(first) then
1041C
1042C correction calcul de LMT a chaque passage pour pouvoir choisir
1043c plusieurs schemas PPM pour differents traceurs
1044c      IF (IORD.LE.0) then
1045c            if(IMR.GE.144) then
1046c                  LMT = 0
1047c            elseif(IMR.GE.72) then
1048c                  LMT = 1
1049c            else
1050c                  LMT = 2
1051c            endif
1052c      else
1053c            LMT = IORD - 3
1054c      endif
1055C
1056      LMT = IORD - 3
1057c      write(6,*) 'PPM option in E-W direction = ', LMT
1058c      first = .false.
1059C      endif
1060C
1061      DO 10 i=1,IMR
106210    AL(i) = 0.5*(p(i-1)+p(i)) + (DC(i-1) - DC(i))*R3
1063C
1064      do 20 i=1,IMR-1
106520    AR(i) = AL(i+1)
1066      AR(IMR) = AL(1)
1067C
1068      do 30 i=1,IMR
106930    A6(i) = 3.*(p(i)+p(i)  - (AL(i)+AR(i)))
1070C
1071      if(LMT.LE.2) call lmtppm(DC(1),A6(1),AR(1),AL(1),P(1),IMR,LMT)
1072C
1073      AL(0) = AL(IMR)
1074      AR(0) = AR(IMR)
1075      A6(0) = A6(IMR)
1076C
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
1088C
1089      subroutine xmist(IMR,IML,P,DC)
1090      parameter( R24 = 1./24.)
1091      dimension P(-IML:IMR+1+IML),DC(-IML:IMR+1+IML)
1092C
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))
109710    DC(i) = sign(min(abs(tmp),Pmax,Pmin), tmp)
1098      return
1099      end
1100C
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)
1105C Work array
1106      DIMENSION fx(IMR,JNP),AR(IMR,JNP),AL(IMR,JNP),A6(IMR,JNP)
1107C
1108      JMR = JNP - 1
1109      len = IMR*(J2-J1+2)
1110C
1111      if(JORD.eq.1) then
1112      DO 1000 i=1,len
1113      JT = REAL(J1) - VC(i,J1)
11141000  fx(i,j1) = p(i,JT)
1115      else
1116   
1117      call ymist(IMR,JNP,j1,P,DC2,4)
1118C
1119      if(JORD.LE.0 .or. JORD.GE.3) then
1120   
1121      call fyppm(VC,P,DC2,fx,IMR,JNP,j1,j2,A6,AR,AL,JORD)
1122   
1123      else
1124      DO 1200 i=1,len
1125      JT = REAL(J1) - VC(i,J1)
11261200  fx(i,j1) = p(i,JT) + (sign(1.,VC(i,j1))-VC(i,j1))*DC2(i,JT)
1127      endif
1128      endif
1129C
1130      DO 1300 i=1,len
11311300  fx(i,j1) = fx(i,j1)*ymass(i,j1)
1132C
1133      DO 1400 j=j1,j2
1134      DO 1400 i=1,IMR
11351400  DQ(i,j) = DQ(i,j) + (fx(i,j) - fx(i,j+1)) * acosp(j)
1136C
1137C 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
1144C
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
1151C
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
1158C
1159      return
1160      end
1161C
1162      subroutine  ymist(IMR,JNP,j1,P,DC,ID)
1163      parameter ( R24 = 1./24. )
1164      dimension P(IMR,JNP),DC(IMR,JNP)
1165C
1166      IMH = IMR / 2
1167      JMR = JNP - 1
1168      IJM3 = IMR*(JMR-3)
1169C
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)
117610    CONTINUE
1177      ELSE
1178      do 12 i=1,IMH
1179C 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)
1184C 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)
118912    CONTINUE
1190      do 14 i=IMH+1,IMR
1191C 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)
1196C 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)
120114    CONTINUE
1202C
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)
120815    CONTINUE
1209      ENDIF
1210C
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
1217C Determine slopes in polar caps for scalars!
1218C
1219      do 13 i=1,IMH
1220C 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)
1225C 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)
123013    continue
1231C
1232      do 25 i=imh+1,IMR
1233      DC(i,  1) =  - DC(i-imh,  1)
1234      DC(i,JNP) =  - DC(i-imh,JNP)
123525    continue
1236      endif
1237      return
1238      end
1239C
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,*)
1243C Local work arrays.
1244      real AR(IMR,JNP),AL(IMR,JNP),A6(IMR,JNP)
1245      integer LMT
1246c      logical first
1247C      data first /.true./
1248C      SAVE LMT
1249C
1250      IMH = IMR / 2
1251      JMR = JNP - 1
1252      j11 = j1-1
1253      IMJM1 = IMR*(J2-J1+2)
1254      len   = IMR*(J2-J1+3)
1255C      if(first) then
1256C      IF(JORD.LE.0) then
1257C            if(JMR.GE.90) then
1258C                  LMT = 0
1259C            elseif(JMR.GE.45) then
1260C                  LMT = 1
1261C            else
1262C                  LMT = 2
1263C            endif
1264C      else
1265C            LMT = JORD - 3
1266C      endif
1267C
1268C      first = .false.
1269C      endif
1270C     
1271c modifs pour pouvoir choisir plusieurs schemas PPM
1272      LMT = JORD - 3     
1273C
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)
127710    CONTINUE
1278C
1279CPoles:
1280C
1281      DO i=1,IMH
1282      AL(i,1) = AL(i+IMH,2)
1283      AL(i+IMH,1) = AL(i,2)
1284C
1285      AR(i,JNP) = AR(i+IMH,JMR)
1286      AR(i+IMH,JNP) = AR(i,JMR)
1287      ENDDO
1288
1289ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1290c   Rajout pour LMDZ.3.3
1291cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1292      AR(IMR,1)=AL(1,1)
1293      AR(IMR,JNP)=AL(1,JNP)
1294ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1295     
1296           
1297      do 30 i=1,len
129830    A6(i,j11) = 3.*(p(i,j11)+p(i,j11)  - (AL(i,j11)+AR(i,j11)))
1299C
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)
1302C
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
1312140   continue
1313      return
1314      end
1315C
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)
1319C
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
1327C 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
1338c        write(*,*) 'toto 1'
1339C --------------------------------
1340      IF(IAD.eq.2) then
1341      do j=j1-1,j2+1
1342      do i=1,IMR
1343c      write(*,*) 'avt NINT','i=',i,'j=',j
1344      JP = NINT(VA(i,j))     
1345      rv = JP - VA(i,j)
1346c      write(*,*) 'VA=',VA(i,j), 'JP1=',JP,'rv=',rv
1347      JP = j - JP
1348c      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))
1351c      write(*,*) 'a1=',a1,'b1=',b1
1352      ady(i,j) = wk(i,jp) + rv*(a1*rv + b1) - wk(i,j)
1353      enddo
1354      enddo
1355c      write(*,*) 'toto 2'
1356C
1357      ELSEIF(IAD.eq.1) then
1358        do j=j1-1,j2+1
1359      do i=1,imr
1360      JP = REAL(j)-VA(i,j)
1361      ady(i,j) = VA(i,j)*(wk(i,jp)-wk(i,jp+1))
1362      enddo
1363      enddo
1364      ENDIF
1365C
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
1375C
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
1383C 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
1392C
1393      do i=1,imr
1394      ady(i,  1) =  sum1
1395      ady(i,JNP) =  sum2
1396      enddo
1397        endif
1398C
1399        return
1400        end
1401C
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)
1404C
1405        JMR = JNP-1
1406      do 1309 j=j1,j2
1407      if(J.GT.JS  .and. J.LT.JN) GO TO 1309
1408C
1409      do i=1,IMR
1410      qtmp(i) = p(i,j)
1411      enddo
1412C
1413      do i=-IML,0
1414      qtmp(i)       = p(IMR+i,j)
1415      qtmp(IMR+1-i) = p(1-i,j)
1416      enddo
1417C
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
1439C
1440      do i=1,IMR
1441      adx(i,j) = adx(i,j) - p(i,j)
1442      enddo
14431309  continue
1444C
1445C Eulerian upwind
1446C
1447      do j=JS+1,JN-1
1448C
1449      do i=1,IMR
1450      qtmp(i) = p(i,j)
1451      enddo
1452C
1453      qtmp(0)     = p(IMR,J)
1454      qtmp(IMR+1) = p(1,J)
1455C
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
1468C 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
1475C
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
1482C 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
1489C
1490      subroutine lmtppm(DC,A6,AR,AL,P,IM,LMT)
1491C
1492C A6 =  CURVATURE OF THE TEST PARABOLA
1493C AR =  RIGHT EDGE VALUE OF THE TEST PARABOLA
1494C AL =  LEFT  EDGE VALUE OF THE TEST PARABOLA
1495C DC =  0.5 * MISMATCH
1496C P  =  CELL-AVERAGED VALUE
1497C IM =  VECTOR LENGTH
1498C
1499C OPTIONS:
1500C
1501C LMT = 0: FULL MONOTONICITY
1502C LMT = 1: SEMI-MONOTONIC CONSTRAINT (NO UNDERSHOOTS)
1503C LMT = 2: POSITIVE-DEFINITE CONSTRAINT
1504C
1505      parameter ( R12 = 1./12. )
1506      dimension A6(IM),AR(IM),AL(IM),P(IM),DC(IM)
1507C
1508      if(LMT.eq.0) then
1509C 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
1527100   continue
1528      elseif(LMT.eq.1) then
1529C 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
1543150   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
1560250   continue
1561      endif
1562      return
1563      end
1564C
1565      subroutine A2C(U,V,IMR,JMR,j1,j2,CRX,CRY,dtdx5,DTDY5)
1566      dimension U(IMR,*),V(IMR,*),CRX(IMR,*),CRY(IMR,*),DTDX5(*)
1567C
1568      do 35 j=j1,j2
1569      do 35 i=2,IMR
157035    CRX(i,J) = dtdx5(j)*(U(i,j)+U(i-1,j))
1571C
1572      do 45 j=j1,j2
157345    CRX(1,J) = dtdx5(j)*(U(1,j)+U(IMR,j))
1574C
1575      do 55 i=1,IMR*JMR
157655    CRY(i,2) = DTDY5*(V(i,2)+V(i,1))
1577      return
1578      end
1579C
1580      subroutine cosa(cosp,cose,JNP,PI,DP)
1581      dimension cosp(*),cose(*)
1582      JMR = JNP-1
1583      do 55 j=2,JNP
1584        ph5  =  -0.5*PI + (REAL(J-1)-0.5)*DP
158555      cose(j) = cos(ph5)
1586C
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
1593C 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
1599C
1600      do 66 j=2,JMR
160166    cosp(j) = 0.5*(cose(j)+cose(j+1))
1602      cosp(1) = 0.
1603      cosp(JNP) = 0.
1604      return
1605      end
1606C
1607      subroutine cosc(cosp,cose,JNP,PI,DP)
1608      dimension cosp(*),cose(*)
1609C
1610      phi = -0.5*PI
1611      do 55 j=2,JNP-1
1612      phi  =  phi + DP
161355    cosp(j) = cos(phi)
1614        cosp(  1) = 0.
1615        cosp(JNP) = 0.
1616C
1617      do 66 j=2,JNP
1618        cose(j) = 0.5*(cosp(j)+cosp(j-1))
161966    CONTINUE
1620C
1621      do 77 j=2,JNP-1
1622       cosp(j) = 0.5*(cose(j)+cose(j+1))
162377    CONTINUE
1624      return
1625      end
1626C
1627      SUBROUTINE qckxyz (Q,qtmp,IMR,JNP,NLAY,j1,j2,cosp,acosp,
1628     &                   cross,IC,NSTEP)
1629C
1630      parameter( tiny = 1.E-60 )
1631      DIMENSION Q(IMR,JNP,NLAY),qtmp(IMR,JNP),cosp(*),acosp(*)
1632      logical cross
1633C
1634      NLAYM1 = NLAY-1
1635      len = IMR*(j2-j1+1)
1636      ip = 0
1637C
1638C 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
1645C
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
1650C
1651C 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
1659C
166050    continue
1661      DO 225 L = 2,NLAYM1
1662      icr = 1
1663C
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
1672C
1673      do i=1,len
1674      IF( Q(I,j1,L).LT.0.) THEN
1675C
1676      ip = ip + 1
1677C 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
1683C 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
1688225   CONTINUE
1689C
1690C BOTTOM LAYER
1691      sum = 0.
1692      L = NLAY
1693C
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
1698C
1699      call filcr(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,icr,tiny)
1700      if(icr.eq.0) goto 911
1701C
1702      DO  I=1,len
1703      IF( Q(I,j1,L).LT.0.) THEN
1704      ip = ip + 1
1705c
1706C From above
1707C
1708          qup = Q(I,j1,NLAYM1)
1709          qly = -Q(I,j1,L)
1710          dup = min(qly,qup)
1711          Q(I,j1,NLAYM1) = qup - dup
1712C From "below" the surface.
1713          sum = sum + qly-dup
1714          Q(I,j1,L) = 0.
1715       ENDIF
1716      ENDDO
1717C
1718911   continue
1719C
1720      if(ip.gt.IMR) then
1721      write(6,*) 'IC=',IC,' STEP=',NSTEP,
1722     &           ' Vertical filling pts=',ip
1723      endif
1724C
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
1730C
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)
1739C 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
1745C 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
175250    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)
1758C 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
1764C 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
177155    continue
1772C *****************************************
1773C i=1
1774      i=1
1775      IF(q(i,j).LT.0.) THEN
1776      icr =  1
1777      dq  = - q(i,j)*cosp(j)
1778C 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
1784C 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
1791C *****************************************
1792C i=IMR
1793      i=IMR
1794      IF(q(i,j).LT.0.) THEN
1795      icr =  1
1796      dq  = - q(i,j)*cosp(j)
1797C 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
1803C 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
1810C *****************************************
181165    continue
1812C
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
1819C
182080    continue
1821C
1822      if(q(1,1).lt.0. .or. q(1,jnp).lt.0.) then
1823      icr = 1
1824      endif
1825C
1826      return
1827      end
1828C
1829      subroutine filns(q,IMR,JNP,j1,j2,cosp,acosp,ipy,tiny)
1830      dimension q(IMR,*),cosp(*),acosp(*)
1831c      logical first
1832c      data first /.true./
1833c      save cap1
1834C
1835c      if(first) then
1836      DP = 4.*ATAN(1.)/REAL(JNP-1)
1837      CAP1 = IMR*(1.-COS((j1-1.5)*DP))/DP
1838c      first = .false.
1839c      endif
1840C
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)
1847C 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
1853C 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
186055    continue
1861C
1862      do i=1,imr
1863      IF(q(i,j1).LT.0.) THEN
1864      ipy =  1
1865      dq  = - q(i,j1)*cosp(j1)
1866C 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
1874C
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)
1880C 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
1888C
1889C Check Poles.
1890      if(q(1,1).lt.0.) then
1891      dq = q(1,1)*cap1/REAL(IMR)*acosp(j1)
1892      do i=1,imr
1893      q(i,1) = 0.
1894      q(i,j1) = q(i,j1) + dq
1895      if(q(i,j1).lt.0.) ipy = 1
1896      enddo
1897      endif
1898C
1899      if(q(1,JNP).lt.0.) then
1900      dq = q(1,JNP)*cap1/REAL(IMR)*acosp(j2)
1901      do i=1,imr
1902      q(i,JNP) = 0.
1903      q(i,j2) = q(i,j2) + dq
1904      if(q(i,j2).lt.0.) ipy = 1
1905      enddo
1906      endif
1907C
1908      return
1909      end
1910C
1911      subroutine filew(q,qtmp,IMR,JNP,j1,j2,ipx,tiny)
1912      dimension q(IMR,*),qtmp(JNP,IMR)
1913C
1914      ipx = 0
1915C Copy & swap direction for vectorization.
1916      do 25 i=1,imr
1917      do 25 j=j1,j2
191825    qtmp(j,i) = q(i,j)
1919C
1920      do 55 i=2,imr-1
1921      do 55 j=j1,j2
1922      if(qtmp(j,i).lt.0.) then
1923      ipx =  1
1924c 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
1929c 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
193555    continue
1936c
1937      i=1
1938      do 65 j=j1,j2
1939      if(qtmp(j,i).lt.0.) then
1940      ipx =  1
1941c 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
1946c 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
1950c
1951      qtmp(j,i) = qtmp(j,i) + d2 + tiny
1952      endif
195365    continue
1954      i=IMR
1955      do 75 j=j1,j2
1956      if(qtmp(j,i).lt.0.) then
1957      ipx =  1
1958c 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
1963c east
1964      d0 = max(0.,qtmp(j,1))
1965      d2 = min(-qtmp(j,i),d0)
1966      qtmp(j,1) = qtmp(j,1) - d2
1967c
1968      qtmp(j,i) = qtmp(j,i) + d2 + tiny
1969      endif
197075    continue
1971C
1972      if(ipx.ne.0) then
1973      do 85 j=j1,j2
1974      do 85 i=1,imr
197585    q(i,j) = qtmp(j,i)
1976      else
1977C
1978C Poles.
1979      if(q(1,1).lt.0. or. q(1,JNP).lt.0.) ipx = 1
1980      endif
1981      return
1982      end
1983C
1984      subroutine zflip(q,im,km,nc)
1985C This routine flip the array q (in the vertical).
1986      real q(im,km,nc)
1987C local dynamic array
1988      real qtmp(im,km)
1989C
1990      do 4000 IC = 1, nc
1991C
1992      do 1000 k=1,km
1993      do 1000 i=1,im
1994      qtmp(i,k) = q(i,km+1-k,IC)
19951000  continue
1996C
1997      do 2000 i=1,im*km
19982000  q(i,1,IC) = qtmp(i,1)
19994000  continue
2000      return
2001      end
Note: See TracBrowser for help on using the repository browser.