source: LMDZ4/branches/V3_test/libf/dyn3d/advyp.F @ 1270

Last change on this file since 1270 was 524, checked in by lmdzadmin, 21 years ago

Initial revision

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