source: LMDZ6/trunk/libf/dyn3d_common/advyp.f90 @ 5248

Last change on this file since 5248 was 5246, checked in by abarral, 20 hours ago

Convert fixed-form to free-form sources .F -> .{f,F}90
(WIP: some .F remain, will be handled in subsequent commits)

  • 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: 17.8 KB
Line 
1!
2! $Header$
3!
4SUBROUTINE ADVYP(LIMIT,DTY,PBARV,SM,S0,SSX,SY,SZ &
5        ,SSXX,SSXY,SSXZ,SYY,SYZ,SZZ,ntra )
6  IMPLICIT NONE
7  !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
8  !                                                                 C
9  !  second-order moments (SOM) advection of tracer in Y direction  C
10  !                                                                 C
11  !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
12                                                             ! C
13  !  Source : Pascal Simon ( Meteo, CNRM )                         C
14  !  Adaptation : A.A. (LGGE)                                      C
15  !  Derniere Modif : 19/10/95 LAST
16                                                             ! C
17  !  sont les arguments d'entree pour le s-pg                      C
18  !                                                                C
19  !  argument de sortie du s-pg                                    C
20  !                                                                C
21  !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
22  !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
23  !
24  !  Rem : Probleme aux poles il faut reecrire ce cas specifique
25  !    Attention au sens de l'indexation
26  !
27  !  parametres principaux du modele
28  !
29  !
30  include "dimensions.h"
31  include "paramet.h"
32  include "comgeom.h"
33
34  !  Arguments :
35  !  ----------
36  !  dty : frequence fictive d'appel du transport
37  !  parbu,pbarv : flux de masse en x et y en Pa.m2.s-1
38
39  INTEGER :: lon,lat,niv
40  INTEGER :: i,j,jv,k,kp,l
41  INTEGER :: ntra
42   ! PARAMETER (ntra = 1)
43
44  REAL :: dty
45  REAL :: pbarv ( iip1,jjm, llm )
46
47  !  moments: SM  total mass in each grid box
48        ! S0  mass of tracer in each grid box
49        ! Si  1rst order moment in i direction
50  !
51  REAL :: SM(iip1,jjp1,llm) &
52        ,S0(iip1,jjp1,llm,ntra)
53  REAL :: SSX(iip1,jjp1,llm,ntra) &
54        ,SY(iip1,jjp1,llm,ntra) &
55        ,SZ(iip1,jjp1,llm,ntra) &
56        ,SSXX(iip1,jjp1,llm,ntra) &
57        ,SSXY(iip1,jjp1,llm,ntra) &
58        ,SSXZ(iip1,jjp1,llm,ntra) &
59        ,SYY(iip1,jjp1,llm,ntra) &
60        ,SYZ(iip1,jjp1,llm,ntra) &
61        ,SZZ(iip1,jjp1,llm,ntra)
62  !
63  !  Local :
64  !  -------
65
66  !  mass fluxes across the boundaries (UGRI,VGRI,WGRI)
67  !  mass fluxes in kg
68  !  declaration :
69
70  REAL :: VGRI(iip1,0:jjp1,llm)
71
72  !  Rem : UGRI et WGRI ne sont pas utilises dans
73  !  cette subroutine ( advection en y uniquement )
74  !  Rem 2 :le dimensionnement de VGRI depend de celui de pbarv
75  !
76  !  the moments F are similarly defined and used as temporary
77  !  storage for portions of the grid boxes in transit
78  !
79  !  the moments Fij are used as temporary storage for
80  !  portions of the grid boxes in transit at the current level
81  !
82  !  work arrays
83  !
84  !
85  REAL :: F0(iim,0:jjp1,ntra),FM(iim,0:jjp1)
86  REAL :: FX(iim,jjm,ntra),FY(iim,jjm,ntra)
87  REAL :: FZ(iim,jjm,ntra)
88  REAL :: FXX(iim,jjm,ntra),FXY(iim,jjm,ntra)
89  REAL :: FXZ(iim,jjm,ntra),FYY(iim,jjm,ntra)
90  REAL :: FYZ(iim,jjm,ntra),FZZ(iim,jjm,ntra)
91  REAL :: S00(ntra)
92  REAL :: SM0             ! Just temporal variable
93  !
94  !  work arrays
95  !
96  REAL :: ALF(iim,0:jjp1),ALF1(iim,0:jjp1)
97  REAL :: ALFQ(iim,0:jjp1),ALF1Q(iim,0:jjp1)
98  REAL :: ALF2(iim,0:jjp1),ALF3(iim,0:jjp1)
99  REAL :: ALF4(iim,0:jjp1)
100  REAL :: TEMPTM          ! Just temporal variable
101  REAL :: SLPMAX,S1MAX,S1NEW,S2NEW
102  !
103  !  Special pour poles
104  !
105  REAL :: sbms,sfms,sfzs,sbmn,sfmn,sfzn
106  REAL :: sns0(ntra),snsz(ntra),snsm
107  REAL :: qy1(iim,llm,ntra),qylat(iim,llm,ntra)
108  REAL :: cx1(llm,ntra), cxLAT(llm,ntra)
109  REAL :: cy1(llm,ntra), cyLAT(llm,ntra)
110  REAL :: z1(iim), zcos(iim), zsin(iim)
111  REAL :: SSUM
112  EXTERNAL SSUM
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).GT.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).LE.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).LE.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).GT.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).GT.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).GT.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).GT.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).LT.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).LT.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).LT.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).LT.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).GE.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).GE.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).LT.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).LT.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).LT.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).LT.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.