source: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/advyp.f90 @ 5209

Last change on this file since 5209 was 5159, checked in by abarral, 7 weeks ago

Put dimensions.h and paramet.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: 20.8 KB
Line 
1! $Header$
2
3SUBROUTINE ADVYP(LIMIT, DTY, PBARV, SM, S0, SSX, SY, SZ &
4        , SSXX, SSXY, SSXZ, SYY, SYZ, SZZ, ntra)
5  USE lmdz_comgeom
6
7  USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm
8  USE lmdz_paramet
9  IMPLICIT NONE
10  !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
11  !                                                                 C
12  !  second-order moments (SOM) advection of tracer in Y direction  C
13  !                                                                 C
14  !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
15  ! C
16  !  Source : Pascal Simon ( Meteo, CNRM )                         C
17  !  Adaptation : A.A. (LGGE)                                      C
18  !  Derniere Modif : 19/10/95 LAST
19  ! C
20  !  sont les arguments d'entree pour le s-pg                      C
21  !                                                                C
22  !  argument de sortie du s-pg                                    C
23  !                                                                C
24  !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
25  !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
26
27  !  Rem : Probleme aux poles il faut reecrire ce cas specifique
28  !    Attention au sens de l'indexation
29
30  !  parametres principaux du modele
31
32  !
33
34
35
36  !  Arguments :
37  !  ----------
38  !  dty : frequence fictive d'appel du transport
39  !  parbu,pbarv : flux de masse en x et y en Pa.m2.s-1
40
41  INTEGER :: lon, lat, niv
42  INTEGER :: i, j, jv, k, kp, l
43  INTEGER :: ntra
44  ! PARAMETER (ntra = 1)
45
46  REAL :: dty
47  REAL :: pbarv (iip1, jjm, llm)
48
49  !  moments: SM  total mass in each grid box
50  ! S0  mass of tracer in each grid box
51  ! Si  1rst order moment in i direction
52
53  REAL :: SM(iip1, jjp1, llm) &
54          , S0(iip1, jjp1, llm, ntra)
55  REAL :: SSX(iip1, jjp1, llm, ntra) &
56          , SY(iip1, jjp1, llm, ntra) &
57          , SZ(iip1, jjp1, llm, ntra) &
58          , SSXX(iip1, jjp1, llm, ntra) &
59          , SSXY(iip1, jjp1, llm, ntra) &
60          , SSXZ(iip1, jjp1, llm, ntra) &
61          , SYY(iip1, jjp1, llm, ntra) &
62          , SYZ(iip1, jjp1, llm, ntra) &
63          , SZZ(iip1, jjp1, llm, ntra)
64
65  !  Local :
66  !  -------
67
68  !  mass fluxes across the boundaries (UGRI,VGRI,WGRI)
69  !  mass fluxes in kg
70  !  declaration :
71
72  REAL :: VGRI(iip1, 0:jjp1, llm)
73
74  !  Rem : UGRI et WGRI ne sont pas utilises dans
75  !  cette SUBROUTINE ( advection en y uniquement )
76  !  Rem 2 :le dimensionnement de VGRI depend de celui de pbarv
77
78  !  the moments F are similarly defined and used as temporary
79  !  storage for portions of the grid boxes in transit
80
81  !  the moments Fij are used as temporary storage for
82  !  portions of the grid boxes in transit at the current level
83
84  !  work arrays
85
86
87  REAL :: F0(iim, 0:jjp1, ntra), FM(iim, 0:jjp1)
88  REAL :: FX(iim, jjm, ntra), FY(iim, jjm, ntra)
89  REAL :: FZ(iim, jjm, ntra)
90  REAL :: FXX(iim, jjm, ntra), FXY(iim, jjm, ntra)
91  REAL :: FXZ(iim, jjm, ntra), FYY(iim, jjm, ntra)
92  REAL :: FYZ(iim, jjm, ntra), FZZ(iim, jjm, ntra)
93  REAL :: S00(ntra)
94  REAL :: SM0             ! Just temporal variable
95
96  !  work arrays
97
98  REAL :: ALF(iim, 0:jjp1), ALF1(iim, 0:jjp1)
99  REAL :: ALFQ(iim, 0:jjp1), ALF1Q(iim, 0:jjp1)
100  REAL :: ALF2(iim, 0:jjp1), ALF3(iim, 0:jjp1)
101  REAL :: ALF4(iim, 0:jjp1)
102  REAL :: TEMPTM          ! Just temporal variable
103  REAL :: SLPMAX, S1MAX, S1NEW, S2NEW
104
105  !  Special pour poles
106
107  REAL :: sbms, sfms, sfzs, sbmn, sfmn, sfzn
108  REAL :: sns0(ntra), snsz(ntra), snsm
109  REAL :: qy1(iim, llm, ntra), qylat(iim, llm, ntra)
110  REAL :: cx1(llm, ntra), cxLAT(llm, ntra)
111  REAL :: cy1(llm, ntra), cyLAT(llm, ntra)
112  REAL :: z1(iim), zcos(iim), zsin(iim)
113
114  REAL :: sqi, sqf
115  LOGICAL :: LIMIT
116
117  lon = iim         ! rem : Il est possible qu'un pbl. arrive ici
118  lat = jjp1        ! a cause des dim. differentes entre les
119  niv = llm         !       tab. S et VGRI
120
121  !-----------------------------------------------------------------
122  ! initialisations
123
124  sbms = 0.
125  sfms = 0.
126  sfzs = 0.
127  sbmn = 0.
128  sfmn = 0.
129  sfzn = 0.
130
131  !-----------------------------------------------------------------
132  ! *** Test : diag de la qtite totale de traceur dans
133  ! l'atmosphere avant l'advection en Y
134
135  sqi = 0.
136  sqf = 0.
137
138  DO l = 1, llm
139    DO j = 1, jjp1
140      DO i = 1, iim
141        sqi = sqi + S0(i, j, l, ntra)
142      END DO
143    END DO
144  END DO
145  PRINT*, '---------- DIAG DANS ADVY - ENTREE --------'
146  PRINT*, 'sqi=', sqi
147
148  !-----------------------------------------------------------------
149  !  Interface : adaptation nouveau modele
150  !  -------------------------------------
151
152  !  Conversion des flux de masses en kg
153  !-AA 20/10/94  le signe -1 est necessaire car indexation opposee
154
155  DO l = 1, llm
156    DO j = 1, jjm
157      DO i = 1, iip1
158        vgri (i, j, llm + 1 - l) = -1. * pbarv (i, j, l)
159      END DO
160    END DO
161  END DO
162
163  !AA Initialisation de flux fictifs aux bords sup. des boites pol.
164
165  DO l = 1, llm
166    DO i = 1, iip1
167      vgri(i, 0, l) = 0.
168      vgri(i, jjp1, l) = 0.
169    ENDDO
170  ENDDO
171
172  !----------------- START HERE -----------------------
173  !  boucle sur les niveaux
174
175  DO L = 1, NIV
176
177    !  place limits on appropriate moments before transport
178    !  (if flux-limiting is to be applied)
179
180    IF(.NOT.LIMIT) GO TO 11
181
182    DO JV = 1, NTRA
183      DO K = 1, LAT
184        DO I = 1, LON
185          IF(S0(I, K, L, JV)>0.) THEN
186            SLPMAX = AMAX1(S0(I, K, L, JV), 0.)
187            S1MAX = 1.5 * SLPMAX
188            S1NEW = AMIN1(S1MAX, AMAX1(-S1MAX, SY(I, K, L, JV)))
189            S2NEW = AMIN1(2. * SLPMAX - ABS(S1NEW) / 3., &
190                    AMAX1(ABS(S1NEW) - SLPMAX, SYY(I, K, L, JV)))
191            SY (I, K, L, JV) = S1NEW
192            SYY(I, K, L, JV) = S2NEW
193            SSXY(I, K, L, JV) = AMIN1(SLPMAX, AMAX1(-SLPMAX, SSXY(I, K, L, JV)))
194            SYZ(I, K, L, JV) = AMIN1(SLPMAX, AMAX1(-SLPMAX, SYZ(I, K, L, JV)))
195          ELSE
196            SY (I, K, L, JV) = 0.
197            SYY(I, K, L, JV) = 0.
198            SSXY(I, K, L, JV) = 0.
199            SYZ(I, K, L, JV) = 0.
200          ENDIF
201        END DO
202      END DO
203    END DO
204
205    11   CONTINUE
206
207    !  le flux a travers le pole Nord est traite separement
208
209    SM0 = 0.
210    DO JV = 1, NTRA
211      S00(JV) = 0.
212    END DO
213
214    DO I = 1, LON
215
216      IF(VGRI(I, 0, L)<=0.) THEN
217        FM(I, 0) = -VGRI(I, 0, L) * DTY
218        ALF(I, 0) = FM(I, 0) / SM(I, 1, L)
219        SM(I, 1, L) = SM(I, 1, L) - FM(I, 0)
220        SM0 = SM0 + FM(I, 0)
221      ENDIF
222
223      ALFQ(I, 0) = ALF(I, 0) * ALF(I, 0)
224      ALF1(I, 0) = 1. - ALF(I, 0)
225      ALF1Q(I, 0) = ALF1(I, 0) * ALF1(I, 0)
226      ALF2(I, 0) = ALF1(I, 0) - ALF(I, 0)
227      ALF3(I, 0) = ALF(I, 0) * ALFQ(I, 0)
228      ALF4(I, 0) = ALF1(I, 0) * ALF1Q(I, 0)
229
230    END DO
231    ! PRINT*,'ADVYP 21'
232
233    DO JV = 1, NTRA
234      DO I = 1, LON
235
236        IF(VGRI(I, 0, L)<=0.) THEN
237
238          F0(I, 0, JV) = ALF(I, 0) * (S0(I, 1, L, JV) - ALF1(I, 0) * &
239                  (SY(I, 1, L, JV) - ALF2(I, 0) * SYY(I, 1, L, JV)))
240
241          S00(JV) = S00(JV) + F0(I, 0, JV)
242          S0 (I, 1, L, JV) = S0(I, 1, L, JV) - F0(I, 0, JV)
243          SY (I, 1, L, JV) = ALF1Q(I, 0) * &
244                  (SY(I, 1, L, JV) + 3. * ALF(I, 0) * SYY(I, 1, L, JV))
245          SYY(I, 1, L, JV) = ALF4 (I, 0) * SYY(I, 1, L, JV)
246          SSX (I, 1, L, JV) = ALF1 (I, 0) * &
247                  (SSX(I, 1, L, JV) + ALF(I, 0) * SSXY(I, 1, L, JV))
248          SZ (I, 1, L, JV) = ALF1 (I, 0) * &
249                  (SZ(I, 1, L, JV) + ALF(I, 0) * SSXZ(I, 1, L, JV))
250          SSXX(I, 1, L, JV) = ALF1 (I, 0) * SSXX(I, 1, L, JV)
251          SSXZ(I, 1, L, JV) = ALF1 (I, 0) * SSXZ(I, 1, L, JV)
252          SZZ(I, 1, L, JV) = ALF1 (I, 0) * SZZ(I, 1, L, JV)
253          SSXY(I, 1, L, JV) = ALF1Q(I, 0) * SSXY(I, 1, L, JV)
254          SYZ(I, 1, L, JV) = ALF1Q(I, 0) * SYZ(I, 1, L, JV)
255
256        ENDIF
257
258      END DO
259    END DO
260
261    DO I = 1, LON
262      IF(VGRI(I, 0, L)>0.) THEN
263        FM(I, 0) = VGRI(I, 0, L) * DTY
264        ALF(I, 0) = FM(I, 0) / SM0
265      ENDIF
266    END DO
267
268    DO JV = 1, NTRA
269      DO I = 1, LON
270        IF(VGRI(I, 0, L)>0.) THEN
271          F0(I, 0, JV) = ALF(I, 0) * S00(JV)
272        ENDIF
273      END DO
274    END DO
275
276    !  puts the temporary moments Fi into appropriate neighboring boxes
277
278    ! PRINT*,'av ADVYP 25'
279    DO I = 1, LON
280
281      IF(VGRI(I, 0, L)>0.) THEN
282        SM(I, 1, L) = SM(I, 1, L) + FM(I, 0)
283        ALF(I, 0) = FM(I, 0) / SM(I, 1, L)
284      ENDIF
285
286      ALFQ(I, 0) = ALF(I, 0) * ALF(I, 0)
287      ALF1(I, 0) = 1. - ALF(I, 0)
288      ALF1Q(I, 0) = ALF1(I, 0) * ALF1(I, 0)
289      ALF2(I, 0) = ALF1(I, 0) - ALF(I, 0)
290      ALF3(I, 0) = ALF1(I, 0) * ALF(I, 0)
291
292    END DO
293    ! PRINT*,'av ADVYP 25'
294
295    DO JV = 1, NTRA
296      DO I = 1, LON
297
298        IF(VGRI(I, 0, L)>0.) THEN
299
300          TEMPTM = ALF(I, 0) * S0(I, 1, L, JV) - ALF1(I, 0) * F0(I, 0, JV)
301          S0 (I, 1, L, JV) = S0(I, 1, L, JV) + F0(I, 0, JV)
302          SYY(I, 1, L, JV) = ALF1Q(I, 0) * SYY(I, 1, L, JV) &
303                  + 5. * (ALF3 (I, 0) * SY (I, 1, L, JV) - ALF2(I, 0) * TEMPTM)
304          SY (I, 1, L, JV) = ALF1 (I, 0) * SY (I, 1, L, JV) + 3. * TEMPTM
305          SSXY(I, 1, L, JV) = ALF1 (I, 0) * SSXY(I, 1, L, JV) + 3. * ALF(I, 0) * SSX(I, 1, L, JV)
306          SYZ(I, 1, L, JV) = ALF1 (I, 0) * SYZ(I, 1, L, JV) + 3. * ALF(I, 0) * SZ(I, 1, L, JV)
307
308        ENDIF
309
310      END DO
311    END DO
312
313    !  calculate flux and moments between adjacent boxes
314    !  1- create temporary moments/masses for partial boxes in transit
315    !  2- reajusts moments remaining in the box
316
317    !  flux from KP to K if V(K).lt.0 and from K to KP if V(K).gt.0
318
319    ! PRINT*,'av ADVYP 30'
320    DO K = 1, LAT - 1
321      KP = K + 1
322      DO I = 1, LON
323
324        IF(VGRI(I, K, L)<0.) THEN
325          FM(I, K) = -VGRI(I, K, L) * DTY
326          ALF(I, K) = FM(I, K) / SM(I, KP, L)
327          SM(I, KP, L) = SM(I, KP, L) - FM(I, K)
328        ELSE
329          FM(I, K) = VGRI(I, K, L) * DTY
330          ALF(I, K) = FM(I, K) / SM(I, K, L)
331          SM(I, K, L) = SM(I, K, L) - FM(I, K)
332        ENDIF
333
334        ALFQ(I, K) = ALF(I, K) * ALF(I, K)
335        ALF1(I, K) = 1. - ALF(I, K)
336        ALF1Q(I, K) = ALF1(I, K) * ALF1(I, K)
337        ALF2(I, K) = ALF1(I, K) - ALF(I, K)
338        ALF3(I, K) = ALF(I, K) * ALFQ(I, K)
339        ALF4(I, K) = ALF1(I, K) * ALF1Q(I, K)
340
341      END DO
342    END DO
343    ! PRINT*,'ap ADVYP 30'
344
345    DO JV = 1, NTRA
346      DO K = 1, LAT - 1
347        KP = K + 1
348        DO I = 1, LON
349
350          IF(VGRI(I, K, L)<0.) THEN
351
352            F0 (I, K, JV) = ALF (I, K) * (S0(I, KP, L, JV) - ALF1(I, K) * &
353                    (SY(I, KP, L, JV) - ALF2(I, K) * SYY(I, KP, L, JV)))
354            FY (I, K, JV) = ALFQ(I, K) * &
355                    (SY(I, KP, L, JV) - 3. * ALF1(I, K) * SYY(I, KP, L, JV))
356            FYY(I, K, JV) = ALF3(I, K) * SYY(I, KP, L, JV)
357            FX (I, K, JV) = ALF (I, K) * &
358                    (SSX(I, KP, L, JV) - ALF1(I, K) * SSXY(I, KP, L, JV))
359            FZ (I, K, JV) = ALF (I, K) * &
360                    (SZ(I, KP, L, JV) - ALF1(I, K) * SYZ(I, KP, L, JV))
361            FXY(I, K, JV) = ALFQ(I, K) * SSXY(I, KP, L, JV)
362            FYZ(I, K, JV) = ALFQ(I, K) * SYZ(I, KP, L, JV)
363            FXX(I, K, JV) = ALF (I, K) * SSXX(I, KP, L, JV)
364            FXZ(I, K, JV) = ALF (I, K) * SSXZ(I, KP, L, JV)
365            FZZ(I, K, JV) = ALF (I, K) * SZZ(I, KP, L, JV)
366
367            S0 (I, KP, L, JV) = S0(I, KP, L, JV) - F0(I, K, JV)
368            SY (I, KP, L, JV) = ALF1Q(I, K) * &
369                    (SY(I, KP, L, JV) + 3. * ALF(I, K) * SYY(I, KP, L, JV))
370            SYY(I, KP, L, JV) = ALF4(I, K) * SYY(I, KP, L, JV)
371            SSX (I, KP, L, JV) = SSX (I, KP, L, JV) - FX (I, K, JV)
372            SZ (I, KP, L, JV) = SZ (I, KP, L, JV) - FZ (I, K, JV)
373            SSXX(I, KP, L, JV) = SSXX(I, KP, L, JV) - FXX(I, K, JV)
374            SSXZ(I, KP, L, JV) = SSXZ(I, KP, L, JV) - FXZ(I, K, JV)
375            SZZ(I, KP, L, JV) = SZZ(I, KP, L, JV) - FZZ(I, K, JV)
376            SSXY(I, KP, L, JV) = ALF1Q(I, K) * SSXY(I, KP, L, JV)
377            SYZ(I, KP, L, JV) = ALF1Q(I, K) * SYZ(I, KP, L, JV)
378
379          ELSE
380
381            F0 (I, K, JV) = ALF (I, K) * (S0(I, K, L, JV) + ALF1(I, K) * &
382                    (SY(I, K, L, JV) + ALF2(I, K) * SYY(I, K, L, JV)))
383            FY (I, K, JV) = ALFQ(I, K) * &
384                    (SY(I, K, L, JV) + 3. * ALF1(I, K) * SYY(I, K, L, JV))
385            FYY(I, K, JV) = ALF3(I, K) * SYY(I, K, L, JV)
386            FX (I, K, JV) = ALF (I, K) * (SSX(I, K, L, JV) + ALF1(I, K) * SSXY(I, K, L, JV))
387            FZ (I, K, JV) = ALF (I, K) * (SZ(I, K, L, JV) + ALF1(I, K) * SYZ(I, K, L, JV))
388            FXY(I, K, JV) = ALFQ(I, K) * SSXY(I, K, L, JV)
389            FYZ(I, K, JV) = ALFQ(I, K) * SYZ(I, K, L, JV)
390            FXX(I, K, JV) = ALF (I, K) * SSXX(I, K, L, JV)
391            FXZ(I, K, JV) = ALF (I, K) * SSXZ(I, K, L, JV)
392            FZZ(I, K, JV) = ALF (I, K) * SZZ(I, K, L, JV)
393
394            S0 (I, K, L, JV) = S0 (I, K, L, JV) - F0 (I, K, JV)
395            SY (I, K, L, JV) = ALF1Q(I, K) * &
396                    (SY(I, K, L, JV) - 3. * ALF(I, K) * SYY(I, K, L, JV))
397            SYY(I, K, L, JV) = ALF4(I, K) * SYY(I, K, L, JV)
398            SSX (I, K, L, JV) = SSX (I, K, L, JV) - FX (I, K, JV)
399            SZ (I, K, L, JV) = SZ (I, K, L, JV) - FZ (I, K, JV)
400            SSXX(I, K, L, JV) = SSXX(I, K, L, JV) - FXX(I, K, JV)
401            SSXZ(I, K, L, JV) = SSXZ(I, K, L, JV) - FXZ(I, K, JV)
402            SZZ(I, K, L, JV) = SZZ(I, K, L, JV) - FZZ(I, K, JV)
403            SSXY(I, K, L, JV) = ALF1Q(I, K) * SSXY(I, K, L, JV)
404            SYZ(I, K, L, JV) = ALF1Q(I, K) * SYZ(I, K, L, JV)
405
406          ENDIF
407
408        END DO
409      END DO
410    END DO
411    ! PRINT*,'ap ADVYP 31'
412
413    !  puts the temporary moments Fi into appropriate neighboring boxes
414
415    DO K = 1, LAT - 1
416      KP = K + 1
417      DO I = 1, LON
418
419        IF(VGRI(I, K, L)<0.) THEN
420          SM(I, K, L) = SM(I, K, L) + FM(I, K)
421          ALF(I, K) = FM(I, K) / SM(I, K, L)
422        ELSE
423          SM(I, KP, L) = SM(I, KP, L) + FM(I, K)
424          ALF(I, K) = FM(I, K) / SM(I, KP, L)
425        ENDIF
426
427        ALFQ(I, K) = ALF(I, K) * ALF(I, K)
428        ALF1(I, K) = 1. - ALF(I, K)
429        ALF1Q(I, K) = ALF1(I, K) * ALF1(I, K)
430        ALF2(I, K) = ALF1(I, K) - ALF(I, K)
431        ALF3(I, K) = ALF1(I, K) * ALF(I, K)
432
433      END DO
434    END DO
435    ! PRINT*,'ap ADVYP 32'
436
437    DO JV = 1, NTRA
438      DO K = 1, LAT - 1
439        KP = K + 1
440        DO I = 1, LON
441
442          IF(VGRI(I, K, L)<0.) THEN
443
444            TEMPTM = -ALF(I, K) * S0(I, K, L, JV) + ALF1(I, K) * F0(I, K, JV)
445            S0 (I, K, L, JV) = S0(I, K, L, JV) + F0(I, K, JV)
446            SYY(I, K, L, JV) = ALFQ(I, K) * FYY(I, K, JV) + ALF1Q(I, K) * SYY(I, K, L, JV) &
447                    + 5. * (ALF3(I, K) * (FY(I, K, JV) - SY(I, K, L, JV)) + ALF2(I, K) * TEMPTM)
448            SY (I, K, L, JV) = ALF(I, K) * FY(I, K, JV) + ALF1(I, K) * SY(I, K, L, JV) &
449                    + 3. * TEMPTM
450            SSXY(I, K, L, JV) = ALF (I, K) * FXY(I, K, JV) + ALF1(I, K) * SSXY(I, K, L, JV) &
451                    + 3. * (ALF1(I, K) * FX (I, K, JV) - ALF (I, K) * SSX (I, K, L, JV))
452            SYZ(I, K, L, JV) = ALF (I, K) * FYZ(I, K, JV) + ALF1(I, K) * SYZ(I, K, L, JV) &
453                    + 3. * (ALF1(I, K) * FZ (I, K, JV) - ALF (I, K) * SZ (I, K, L, JV))
454            SSX (I, K, L, JV) = SSX (I, K, L, JV) + FX (I, K, JV)
455            SZ (I, K, L, JV) = SZ (I, K, L, JV) + FZ (I, K, JV)
456            SSXX(I, K, L, JV) = SSXX(I, K, L, JV) + FXX(I, K, JV)
457            SSXZ(I, K, L, JV) = SSXZ(I, K, L, JV) + FXZ(I, K, JV)
458            SZZ(I, K, L, JV) = SZZ(I, K, L, JV) + FZZ(I, K, JV)
459
460          ELSE
461
462            TEMPTM = ALF(I, K) * S0(I, KP, L, JV) - ALF1(I, K) * F0(I, K, JV)
463            S0 (I, KP, L, JV) = S0(I, KP, L, JV) + F0(I, K, JV)
464            SYY(I, KP, L, JV) = ALFQ(I, K) * FYY(I, K, JV) + ALF1Q(I, K) * SYY(I, KP, L, JV) &
465                    + 5. * (ALF3(I, K) * (SY(I, KP, L, JV) - FY(I, K, JV)) - ALF2(I, K) * TEMPTM)
466            SY (I, KP, L, JV) = ALF(I, K) * FY(I, K, JV) + ALF1(I, K) * SY(I, KP, L, JV) &
467                    + 3. * TEMPTM
468            SSXY(I, KP, L, JV) = ALF(I, K) * FXY(I, K, JV) + ALF1(I, K) * SSXY(I, KP, L, JV) &
469                    + 3. * (ALF(I, K) * SSX(I, KP, L, JV) - ALF1(I, K) * FX(I, K, JV))
470            SYZ(I, KP, L, JV) = ALF(I, K) * FYZ(I, K, JV) + ALF1(I, K) * SYZ(I, KP, L, JV) &
471                    + 3. * (ALF(I, K) * SZ(I, KP, L, JV) - ALF1(I, K) * FZ(I, K, JV))
472            SSX (I, KP, L, JV) = SSX (I, KP, L, JV) + FX (I, K, JV)
473            SZ (I, KP, L, JV) = SZ (I, KP, L, JV) + FZ (I, K, JV)
474            SSXX(I, KP, L, JV) = SSXX(I, KP, L, JV) + FXX(I, K, JV)
475            SSXZ(I, KP, L, JV) = SSXZ(I, KP, L, JV) + FXZ(I, K, JV)
476            SZZ(I, KP, L, JV) = SZZ(I, KP, L, JV) + FZZ(I, K, JV)
477
478          ENDIF
479
480        END DO
481      END DO
482    END DO
483    ! PRINT*,'ap ADVYP 33'
484
485    !  traitement special pour le pole Sud (idem pole Nord)
486
487    K = LAT
488
489    SM0 = 0.
490    DO JV = 1, NTRA
491      S00(JV) = 0.
492    END DO
493
494    DO I = 1, LON
495
496      IF(VGRI(I, K, L)>=0.) THEN
497        FM(I, K) = VGRI(I, K, L) * DTY
498        ALF(I, K) = FM(I, K) / SM(I, K, L)
499        SM(I, K, L) = SM(I, K, L) - FM(I, K)
500        SM0 = SM0 + FM(I, K)
501      ENDIF
502
503      ALFQ(I, K) = ALF(I, K) * ALF(I, K)
504      ALF1(I, K) = 1. - ALF(I, K)
505      ALF1Q(I, K) = ALF1(I, K) * ALF1(I, K)
506      ALF2(I, K) = ALF1(I, K) - ALF(I, K)
507      ALF3(I, K) = ALF(I, K) * ALFQ(I, K)
508      ALF4(I, K) = ALF1(I, K) * ALF1Q(I, K)
509
510    END DO
511    ! PRINT*,'ap ADVYP 41'
512
513    DO JV = 1, NTRA
514      DO I = 1, LON
515
516        IF(VGRI(I, K, L)>=0.) THEN
517          F0 (I, K, JV) = ALF(I, K) * (S0(I, K, L, JV) + ALF1(I, K) * &
518                  (SY(I, K, L, JV) + ALF2(I, K) * SYY(I, K, L, JV)))
519          S00(JV) = S00(JV) + F0(I, K, JV)
520
521          S0 (I, K, L, JV) = S0 (I, K, L, JV) - F0 (I, K, JV)
522          SY (I, K, L, JV) = ALF1Q(I, K) * &
523                  (SY(I, K, L, JV) - 3. * ALF(I, K) * SYY(I, K, L, JV))
524          SYY(I, K, L, JV) = ALF4 (I, K) * SYY(I, K, L, JV)
525          SSX (I, K, L, JV) = ALF1(I, K) * (SSX(I, K, L, JV) - ALF(I, K) * SSXY(I, K, L, JV))
526          SZ (I, K, L, JV) = ALF1(I, K) * (SZ(I, K, L, JV) - ALF(I, K) * SYZ(I, K, L, JV))
527          SSXX(I, K, L, JV) = ALF1 (I, K) * SSXX(I, K, L, JV)
528          SSXZ(I, K, L, JV) = ALF1 (I, K) * SSXZ(I, K, L, JV)
529          SZZ(I, K, L, JV) = ALF1 (I, K) * SZZ(I, K, L, JV)
530          SSXY(I, K, L, JV) = ALF1Q(I, K) * SSXY(I, K, L, JV)
531          SYZ(I, K, L, JV) = ALF1Q(I, K) * SYZ(I, K, L, JV)
532        ENDIF
533
534      END DO
535    END DO
536    ! PRINT*,'ap ADVYP 42'
537
538    DO I = 1, LON
539      IF(VGRI(I, K, L)<0.) THEN
540        FM(I, K) = -VGRI(I, K, L) * DTY
541        ALF(I, K) = FM(I, K) / SM0
542      ENDIF
543    END DO
544    ! PRINT*,'ap ADVYP 43'
545
546    DO JV = 1, NTRA
547      DO I = 1, LON
548        IF(VGRI(I, K, L)<0.) THEN
549          F0(I, K, JV) = ALF(I, K) * S00(JV)
550        ENDIF
551      END DO
552    END DO
553
554    !  puts the temporary moments Fi into appropriate neighboring boxes
555
556    DO I = 1, LON
557
558      IF(VGRI(I, K, L)<0.) THEN
559        SM(I, K, L) = SM(I, K, L) + FM(I, K)
560        ALF(I, K) = FM(I, K) / SM(I, K, L)
561      ENDIF
562
563      ALFQ(I, K) = ALF(I, K) * ALF(I, K)
564      ALF1(I, K) = 1. - ALF(I, K)
565      ALF1Q(I, K) = ALF1(I, K) * ALF1(I, K)
566      ALF2(I, K) = ALF1(I, K) - ALF(I, K)
567      ALF3(I, K) = ALF1(I, K) * ALF(I, K)
568
569    END DO
570    ! PRINT*,'ap ADVYP 45'
571
572    DO JV = 1, NTRA
573      DO I = 1, LON
574
575        IF(VGRI(I, K, L)<0.) THEN
576
577          TEMPTM = -ALF(I, K) * S0(I, K, L, JV) + ALF1(I, K) * F0(I, K, JV)
578          S0 (I, K, L, JV) = S0(I, K, L, JV) + F0(I, K, JV)
579          SYY(I, K, L, JV) = ALF1Q(I, K) * SYY(I, K, L, JV) &
580                  + 5. * (-ALF3 (I, K) * SY (I, K, L, JV) + ALF2(I, K) * TEMPTM)
581          SY (I, K, L, JV) = ALF1(I, K) * SY (I, K, L, JV) + 3. * TEMPTM
582          SSXY(I, K, L, JV) = ALF1(I, K) * SSXY(I, K, L, JV) - 3. * ALF(I, K) * SSX(I, K, L, JV)
583          SYZ(I, K, L, JV) = ALF1(I, K) * SYZ(I, K, L, JV) - 3. * ALF(I, K) * SZ(I, K, L, JV)
584
585        ENDIF
586
587      END DO
588    END DO
589    ! PRINT*,'ap ADVYP 46'
590
591  END DO
592
593  !--------------------------------------------------
594  ! bouclage cyclique horizontal .
595
596  DO l = 1, llm
597    DO jv = 1, ntra
598      DO j = 1, jjp1
599        SM(iip1, j, l) = SM(1, j, l)
600        S0(iip1, j, l, jv) = S0(1, j, l, jv)
601        SSX(iip1, j, l, jv) = SSX(1, j, l, jv)
602        SY(iip1, j, l, jv) = SY(1, j, l, jv)
603        SZ(iip1, j, l, jv) = SZ(1, j, l, jv)
604      END DO
605    END DO
606  END DO
607
608  ! -------------------------------------------------------------------
609  ! *** Test  negativite:
610
611  ! DO jv = 1,ntra
612  !  DO l = 1,llm
613  !    DO j = 1,jjp1
614  !      DO i = 1,iip1
615  !         IF (s0( i,j,l,jv ).lt.0.) THEN
616  !            PRINT*, '------ S0 < 0 en FIN ADVYP ---'
617  !            PRINT*, 'S0(',i,j,l,jv,')=', S0(i,j,l,jv)
618  !c                 STOP
619  !         ENDIF
620  !      ENDDO
621  !    ENDDO
622  !  ENDDO
623  ! ENDDO
624
625
626  ! -------------------------------------------------------------------
627  ! *** Test : diag de la qtite totale de traceur dans
628  !       l'atmosphere avant l'advection en Y
629
630  DO l = 1, llm
631    DO j = 1, jjp1
632      DO i = 1, iim
633        sqf = sqf + S0(i, j, l, ntra)
634      END DO
635    END DO
636  END DO
637  PRINT*, '---------- DIAG DANS ADVY - SORTIE --------'
638  PRINT*, 'sqf=', sqf
639  ! PRINT*,'ap ADVYP fin'
640
641  !-----------------------------------------------------------------
642
643  RETURN
644END SUBROUTINE ADVYP
645
646
647
648
649
650
651
652
653
654
655
656
Note: See TracBrowser for help on using the repository browser.