/[lmdze]/trunk/dyn3d/advyp.f
ViewVC logotype

Diff of /trunk/dyn3d/advyp.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

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

Legend:
Removed from v.76  
changed lines
  Added in v.81

  ViewVC Help
Powered by ViewVC 1.1.21