source: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/ppm3d.f90 @ 5221

Last change on this file since 5221 was 5160, checked in by abarral, 5 months ago

Put .h into modules

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