source: trunk/LMDZ.COMMON/libf/dyn3d_common/advyp.F @ 3493

Last change on this file since 3493 was 1422, checked in by milmd, 10 years ago

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