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

Diff of /trunk/dyn3d/advxp.f

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

trunk/dyn3d/advxp.f revision 76 by guez, Fri Nov 15 18:45:49 2013 UTC trunk/dyn3d/advxp.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/advxp.F,v 1.1.1.1 2004/05/19 12:53:06 lmdzadmin Exp $  ! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/advxp.F,v 1.1.1.1 2004/05/19
3  !  ! 12:53:06 lmdzadmin Exp $
4         SUBROUTINE ADVXP(LIMIT,DTX,PBARU,SM,S0,SSX,SY,SZ  
5       .                ,SSXX,SSXY,SSXZ,SYY,SYZ,SZZ,ntra)  SUBROUTINE advxp(limit, dtx, pbaru, sm, s0, ssx, sy, sz, ssxx, ssxy, ssxz, &
6         use dimens_m      syy, syz, szz, ntra)
7        use paramet_m    USE dimens_m
8        use comconst    USE paramet_m
9        use disvert_m    USE comconst
10         IMPLICIT NONE    USE disvert_m
11  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC    IMPLICIT NONE
12  C                                                                 C    ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
13  C  second-order moments (SOM) advection of tracer in X direction  C    ! C
14  C                                                                 C    ! second-order moments (SOM) advection of tracer in X direction  C
15  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC    ! C
16  C    ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
17  C  parametres principaux du modele  
18  C    ! parametres principaux du modele
19    
20         INTEGER ntra  
21  c      PARAMETER (ntra = 1)    INTEGER ntra
22  C    ! PARAMETER (ntra = 1)
23  C  definition de la grille du modele  
24  C    ! definition de la grille du modele
25        REAL dtx  
26        REAL, intent(in):: pbaru ( iip1,jjp1,llm )    REAL dtx
27  C    REAL, INTENT (IN) :: pbaru(iip1, jjp1, llm)
28  C  moments: SM  total mass in each grid box  
29  C           S0  mass of tracer in each grid box    ! moments: SM  total mass in each grid box
30  C           Si  1rst order moment in i direction    ! S0  mass of tracer in each grid box
31  C           Sij 2nd  order moment in i and j directions    ! Si  1rst order moment in i direction
32  C    ! Sij 2nd  order moment in i and j directions
33        REAL SM(iip1,jjp1,llm)  
34       +    ,S0(iip1,jjp1,llm,ntra)    REAL sm(iip1, jjp1, llm), s0(iip1, jjp1, llm, ntra)
35        REAL SSX(iip1,jjp1,llm,ntra)    REAL ssx(iip1, jjp1, llm, ntra), sy(iip1, jjp1, llm, ntra), &
36       +    ,SY(iip1,jjp1,llm,ntra)      sz(iip1, jjp1, llm, ntra)
37       +    ,SZ(iip1,jjp1,llm,ntra)    REAL ssxx(iip1, jjp1, llm, ntra), ssxy(iip1, jjp1, llm, ntra), &
38        REAL SSXX(iip1,jjp1,llm,ntra)      ssxz(iip1, jjp1, llm, ntra), syy(iip1, jjp1, llm, ntra), &
39       +    ,SSXY(iip1,jjp1,llm,ntra)      syz(iip1, jjp1, llm, ntra), szz(iip1, jjp1, llm, ntra)
40       +    ,SSXZ(iip1,jjp1,llm,ntra)  
41       +    ,SYY(iip1,jjp1,llm,ntra)    ! Local :
42       +    ,SYZ(iip1,jjp1,llm,ntra)    ! -------
43       +    ,SZZ(iip1,jjp1,llm,ntra)  
44      ! mass fluxes across the boundaries (UGRI,VGRI,WGRI)
45  C  Local :    ! mass fluxes in kg
46  C  -------    ! declaration :
47    
48  C  mass fluxes across the boundaries (UGRI,VGRI,WGRI)    REAL ugri(iip1, jjp1, llm)
49  C  mass fluxes in kg  
50  C  declaration :    ! Rem : VGRI et WGRI ne sont pas utilises dans
51      ! cette subroutine ( advection en x uniquement )
52         REAL UGRI(iip1,jjp1,llm)  
53    
54  C  Rem : VGRI et WGRI ne sont pas utilises dans    ! Tij are the moments for the current latitude and level
55  C  cette subroutine ( advection en x uniquement )  
56  C    REAL tm(iim)
57  C    REAL t0(iim, ntra), tx(iim, ntra)
58  C  Tij are the moments for the current latitude and level    REAL ty(iim, ntra), tz(iim, ntra)
59  C    REAL txx(iim, ntra), txy(iim, ntra)
60        REAL TM (iim)    REAL txz(iim, ntra), tyy(iim, ntra)
61        REAL T0 (iim,NTRA),TX (iim,NTRA)    REAL tyz(iim, ntra), tzz(iim, ntra)
62        REAL TY (iim,NTRA),TZ (iim,NTRA)  
63        REAL TXX(iim,NTRA),TXY(iim,NTRA)    ! the moments F are similarly defined and used as temporary
64        REAL TXZ(iim,NTRA),TYY(iim,NTRA)    ! storage for portions of the grid boxes in transit
65        REAL TYZ(iim,NTRA),TZZ(iim,NTRA)  
66  C    REAL fm(iim)
67  C  the moments F are similarly defined and used as temporary    REAL f0(iim, ntra), fx(iim, ntra)
68  C  storage for portions of the grid boxes in transit    REAL fy(iim, ntra), fz(iim, ntra)
69  C    REAL fxx(iim, ntra), fxy(iim, ntra)
70        REAL FM (iim)    REAL fxz(iim, ntra), fyy(iim, ntra)
71        REAL F0 (iim,NTRA),FX (iim,NTRA)    REAL fyz(iim, ntra), fzz(iim, ntra)
72        REAL FY (iim,NTRA),FZ (iim,NTRA)  
73        REAL FXX(iim,NTRA),FXY(iim,NTRA)    ! work arrays
74        REAL FXZ(iim,NTRA),FYY(iim,NTRA)  
75        REAL FYZ(iim,NTRA),FZZ(iim,NTRA)    REAL alf(iim), alf1(iim), alfq(iim), alf1q(iim)
76  C    REAL alf2(iim), alf3(iim), alf4(iim)
77  C  work arrays  
78  C    REAL smnew(iim), uext(iim)
79        REAL ALF (iim),ALF1(iim),ALFQ(iim),ALF1Q(iim)    REAL sqi, sqf
80        REAL ALF2(iim),ALF3(iim),ALF4(iim)    REAL temptm
81  C    REAL slpmax
82        REAL SMNEW(iim),UEXT(iim)    REAL s1max, s1new, s2new
83        REAL sqi,sqf  
84        REAL TEMPTM    LOGICAL limit
85        REAL SLPMAX    INTEGER num(jjp1), lonk, numk
86        REAL S1MAX,S1NEW,S2NEW    INTEGER lon, lati, latf, niv
87      INTEGER i, i2, i3, j, jv, l, k, iter
88        LOGICAL LIMIT  
89        INTEGER NUM(jjp1),LONK,NUMK    lon = iim
90        INTEGER lon,lati,latf,niv    lati = 2
91        INTEGER i,i2,i3,j,jv,l,k,iter    latf = jjm
92      niv = llm
93        lon = iim  
94        lati=2    ! *** Test : diagnostique de la qtite totale de traceur
95        latf = jjm    ! dans l'atmosphere avant l'advection
96        niv = llm  
97      sqi = 0.
98  C *** Test : diagnostique de la qtite totale de traceur    sqf = 0.
99  C            dans l'atmosphere avant l'advection  
100  c    DO l = 1, llm
101        sqi =0.      DO j = 1, jjp1
       sqf =0.  
 c  
       DO l = 1, llm  
       DO j = 1, jjp1  
102        DO i = 1, iim        DO i = 1, iim
103           sqi = sqi + S0(i,j,l,ntra)          sqi = sqi + s0(i, j, l, ntra)
104        END DO        END DO
105        END DO
106      END DO
107      PRINT *, '------ DIAG DANS ADVX2 - ENTREE -----'
108      PRINT *, 'sqi=', sqi
109      ! test
110      ! -------------------------------------
111      DO j = 1, jjp1
112        num(j) = 1
113      END DO
114      ! DO l=1,llm
115      ! NUM(2,l)=6
116      ! NUM(3,l)=6
117      ! NUM(jjm-1,l)=6
118      ! NUM(jjm,l)=6
119      ! ENDDO
120      ! DO j=2,6
121      ! NUM(j)=12
122      ! ENDDO
123      ! DO j=jjm-5,jjm-1
124      ! NUM(j)=12
125      ! ENDDO
126    
127      ! Interface : adaptation nouveau modele
128      ! -------------------------------------
129    
130      ! ---------------------------------------------------------
131      ! Conversion des flux de masses en kg/s
132      ! pbaru est en N/s d'ou :
133      ! ugri est en kg/s
134    
135      DO l = 1, llm
136        DO j = 1, jjp1
137          DO i = 1, iip1
138            ugri(i, j, llm+1-l) = pbaru(i, j, l)
139        END DO        END DO
140        END DO
141      END DO
142    
143      ! ---------------------------------------------------------
144      ! start here
145    
146      ! boucle principale sur les niveaux et les latitudes
147    
148      DO l = 1, niv
149        DO k = lati, latf
150    
151    
152          ! initialisation
153    
154          ! program assumes periodic boundaries in X
155    
156          DO i = 2, lon
157            smnew(i) = sm(i, k, l) + (ugri(i-1,k,l)-ugri(i,k,l))*dtx
158        END DO        END DO
159        PRINT*,'------ DIAG DANS ADVX2 - ENTREE -----'        smnew(1) = sm(1, k, l) + (ugri(lon,k,l)-ugri(1,k,l))*dtx
160        PRINT*,'sqi=',sqi  
161  c test        ! modifications for extended polar zones
162  c  -------------------------------------  
163          DO 300 j =1,jjp1        numk = num(k)
164           NUM(j) =1        lonk = lon/numk
165   300  CONTINUE  
166  c       DO l=1,llm        IF (numk>1) THEN
167  c      NUM(2,l)=6  
168  c      NUM(3,l)=6          DO i = 1, lon
169  c      NUM(jjm-1,l)=6              tm(i) = 0.
170  c      NUM(jjm,l)=6          END DO
171  c      ENDDO          DO jv = 1, ntra
172  c        DO j=2,6            DO i = 1, lon
173  c       NUM(j)=12              t0(i, jv) = 0.
174  c       ENDDO              tx(i, jv) = 0.
175  c       DO j=jjm-5,jjm-1              ty(i, jv) = 0.
176  c       NUM(j)=12              tz(i, jv) = 0.
177  c       ENDDO              txx(i, jv) = 0.
178                txy(i, jv) = 0.
179  C  Interface : adaptation nouveau modele              txz(i, jv) = 0.
180  C  -------------------------------------              tyy(i, jv) = 0.
181  C              tyz(i, jv) = 0.
182  C  ---------------------------------------------------------              tzz(i, jv) = 0.
183  C  Conversion des flux de masses en kg/s            END DO
184  C  pbaru est en N/s d'ou :          END DO
185  C  ugri est en kg/s  
186            DO i2 = 1, numk
187         DO 500 l = 1,llm  
188         DO 500 j = 1,jjp1            DO i = 1, lonk
189         DO 500 i = 1,iip1              i3 = (i-1)*numk + i2
190         ugri (i,j,llm+1-l) =pbaru (i,j,l)              tm(i) = tm(i) + sm(i3, k, l)
191   500   CONTINUE              alf(i) = sm(i3, k, l)/tm(i)
192                alf1(i) = 1. - alf(i)
193  C  ---------------------------------------------------------              alfq(i) = alf(i)*alf(i)
194  C  start here              alf1q(i) = alf1(i)*alf1(i)
195  C              alf2(i) = alf1(i) - alf(i)
196  C  boucle principale sur les niveaux et les latitudes              alf3(i) = alf(i)*alf1(i)
197  C                END DO
198        DO 1 L=1,NIV  
199        DO 1 K=lati,latf            DO jv = 1, ntra
200                DO i = 1, lonk
201  C                i3 = (i-1)*numk + i2
202  C  initialisation                temptm = -alf(i)*t0(i, jv) + alf1(i)*s0(i3, k, l, jv)
203  C                t0(i, jv) = t0(i, jv) + s0(i3, k, l, jv)
204  C  program assumes periodic boundaries in X                txx(i, jv) = alfq(i)*ssxx(i3, k, l, jv) + alf1q(i)*txx(i, jv) + &
205  C                  5.*(alf3(i)*(ssx(i3,k,l,jv)-tx(i,jv))+alf2(i)*temptm)
206        DO 10 I=2,LON                tx(i, jv) = alf(i)*ssx(i3, k, l, jv) + alf1(i)*tx(i, jv) + &
207           SMNEW(I)=SM(I,K,L)+(UGRI(I-1,K,L)-UGRI(I,K,L))*DTX                  3.*temptm
208   10   CONTINUE                txy(i, jv) = alf(i)*ssxy(i3, k, l, jv) + alf1(i)*txy(i, jv) + &
209        SMNEW(1)=SM(1,K,L)+(UGRI(LON,K,L)-UGRI(1,K,L))*DTX                  3.*(alf1(i)*sy(i3,k,l,jv)-alf(i)*ty(i,jv))
210  C                txz(i, jv) = alf(i)*ssxz(i3, k, l, jv) + alf1(i)*txz(i, jv) + &
211  C  modifications for extended polar zones                  3.*(alf1(i)*sz(i3,k,l,jv)-alf(i)*tz(i,jv))
212  C                ty(i, jv) = ty(i, jv) + sy(i3, k, l, jv)
213        NUMK=NUM(K)                tz(i, jv) = tz(i, jv) + sz(i3, k, l, jv)
214        LONK=LON/NUMK                tyy(i, jv) = tyy(i, jv) + syy(i3, k, l, jv)
215  C                tyz(i, jv) = tyz(i, jv) + syz(i3, k, l, jv)
216        IF(NUMK.GT.1) THEN                tzz(i, jv) = tzz(i, jv) + szz(i3, k, l, jv)
217  C              END DO
218        DO 111 I=1,LON            END DO
219           TM(I)=0.  
220   111  CONTINUE          END DO
221        DO 112 JV=1,NTRA  
       DO 1120 I=1,LON  
          T0 (I,JV)=0.  
          TX (I,JV)=0.  
          TY (I,JV)=0.  
          TZ (I,JV)=0.  
          TXX(I,JV)=0.  
          TXY(I,JV)=0.  
          TXZ(I,JV)=0.  
          TYY(I,JV)=0.  
          TYZ(I,JV)=0.  
          TZZ(I,JV)=0.  
  1120 CONTINUE  
  112  CONTINUE  
 C  
       DO 11 I2=1,NUMK  
 C  
          DO 113 I=1,LONK  
             I3=(I-1)*NUMK+I2  
             TM(I)=TM(I)+SM(I3,K,L)  
             ALF(I)=SM(I3,K,L)/TM(I)  
             ALF1(I)=1.-ALF(I)  
             ALFQ(I)=ALF(I)*ALF(I)  
             ALF1Q(I)=ALF1(I)*ALF1(I)  
             ALF2(I)=ALF1(I)-ALF(I)  
             ALF3(I)=ALF(I)*ALF1(I)  
  113     CONTINUE  
 C  
          DO 114 JV=1,NTRA  
          DO 1140 I=1,LONK  
             I3=(I-1)*NUMK+I2  
             TEMPTM=-ALF(I)*T0(I,JV)+ALF1(I)*S0(I3,K,L,JV)  
             T0 (I,JV)=T0(I,JV)+S0(I3,K,L,JV)  
             TXX(I,JV)=ALFQ(I)*SSXX(I3,K,L,JV)+ALF1Q(I)*TXX(I,JV)  
      +        +5.*( ALF3(I)*(SSX(I3,K,L,JV)-TX(I,JV))+ALF2(I)*TEMPTM )  
             TX (I,JV)=ALF(I)*SSX(I3,K,L,JV)+ALF1(I)*TX(I,JV)+3.*TEMPTM  
             TXY(I,JV)=ALF (I)*SSXY(I3,K,L,JV)+ALF1(I)*TXY(I,JV)  
      +           +3.*(ALF1(I)*SY (I3,K,L,JV)-ALF (I)*TY (I,JV))  
             TXZ(I,JV)=ALF (I)*SSXZ(I3,K,L,JV)+ALF1(I)*TXZ(I,JV)  
      +           +3.*(ALF1(I)*SZ (I3,K,L,JV)-ALF (I)*TZ (I,JV))  
             TY (I,JV)=TY (I,JV)+SY (I3,K,L,JV)  
             TZ (I,JV)=TZ (I,JV)+SZ (I3,K,L,JV)  
             TYY(I,JV)=TYY(I,JV)+SYY(I3,K,L,JV)  
             TYZ(I,JV)=TYZ(I,JV)+SYZ(I3,K,L,JV)  
             TZZ(I,JV)=TZZ(I,JV)+SZZ(I3,K,L,JV)  
  1140    CONTINUE  
  114     CONTINUE  
 C  
  11   CONTINUE  
 C  
       ELSE  
 C  
       DO 115 I=1,LON  
          TM(I)=SM(I,K,L)  
  115  CONTINUE  
       DO 116 JV=1,NTRA  
       DO 1160 I=1,LON  
          T0 (I,JV)=S0 (I,K,L,JV)  
          TX (I,JV)=SSX (I,K,L,JV)  
          TY (I,JV)=SY (I,K,L,JV)  
          TZ (I,JV)=SZ (I,K,L,JV)  
          TXX(I,JV)=SSXX(I,K,L,JV)  
          TXY(I,JV)=SSXY(I,K,L,JV)  
          TXZ(I,JV)=SSXZ(I,K,L,JV)  
          TYY(I,JV)=SYY(I,K,L,JV)  
          TYZ(I,JV)=SYZ(I,K,L,JV)  
          TZZ(I,JV)=SZZ(I,K,L,JV)  
  1160 CONTINUE  
  116  CONTINUE  
 C  
       ENDIF  
 C  
       DO 117 I=1,LONK  
          UEXT(I)=UGRI(I*NUMK,K,L)  
  117  CONTINUE  
 C  
 C  place limits on appropriate moments before transport  
 C      (if flux-limiting is to be applied)  
 C  
       IF(.NOT.LIMIT) GO TO 13  
 C  
       DO 12 JV=1,NTRA  
       DO 120 I=1,LONK  
         IF(T0(I,JV).GT.0.) THEN  
           SLPMAX=T0(I,JV)  
           S1MAX=1.5*SLPMAX  
           S1NEW=AMIN1(S1MAX,AMAX1(-S1MAX,TX(I,JV)))  
           S2NEW=AMIN1( 2.*SLPMAX-ABS(S1NEW)/3. ,  
      +                 AMAX1(ABS(S1NEW)-SLPMAX,TXX(I,JV)) )  
           TX (I,JV)=S1NEW  
           TXX(I,JV)=S2NEW  
           TXY(I,JV)=AMIN1(SLPMAX,AMAX1(-SLPMAX,TXY(I,JV)))  
           TXZ(I,JV)=AMIN1(SLPMAX,AMAX1(-SLPMAX,TXZ(I,JV)))  
         ELSE  
           TX (I,JV)=0.  
           TXX(I,JV)=0.  
           TXY(I,JV)=0.  
           TXZ(I,JV)=0.  
         ENDIF  
  120  CONTINUE  
  12   CONTINUE  
 C  
  13   CONTINUE  
 C  
 C  calculate flux and moments between adjacent boxes  
 C  1- create temporary moments/masses for partial boxes in transit  
 C  2- reajusts moments remaining in the box  
 C  
 C  flux from IP to I if U(I).lt.0  
 C  
       DO 140 I=1,LONK-1  
          IF(UEXT(I).LT.0.) THEN  
            FM(I)=-UEXT(I)*DTX  
            ALF(I)=FM(I)/TM(I+1)  
            TM(I+1)=TM(I+1)-FM(I)  
          ENDIF  
  140  CONTINUE  
 C  
       I=LONK  
       IF(UEXT(I).LT.0.) THEN  
         FM(I)=-UEXT(I)*DTX  
         ALF(I)=FM(I)/TM(1)  
         TM(1)=TM(1)-FM(I)  
       ENDIF  
 C  
 C  flux from I to IP if U(I).gt.0  
 C  
       DO 141 I=1,LONK  
          IF(UEXT(I).GE.0.) THEN  
            FM(I)=UEXT(I)*DTX  
            ALF(I)=FM(I)/TM(I)  
            TM(I)=TM(I)-FM(I)  
          ENDIF  
  141  CONTINUE  
 C  
       DO 142 I=1,LONK  
          ALFQ(I)=ALF(I)*ALF(I)  
          ALF1(I)=1.-ALF(I)  
          ALF1Q(I)=ALF1(I)*ALF1(I)  
          ALF2(I)=ALF1(I)-ALF(I)  
          ALF3(I)=ALF(I)*ALFQ(I)  
          ALF4(I)=ALF1(I)*ALF1Q(I)  
  142  CONTINUE  
 C  
       DO 150 JV=1,NTRA  
       DO 1500 I=1,LONK-1  
 C  
          IF(UEXT(I).LT.0.) THEN  
 C  
            F0 (I,JV)=ALF (I)* ( T0(I+1,JV)-ALF1(I)*  
      +             ( TX(I+1,JV)-ALF2(I)*TXX(I+1,JV) ) )  
            FX (I,JV)=ALFQ(I)*(TX(I+1,JV)-3.*ALF1(I)*TXX(I+1,JV))  
            FXX(I,JV)=ALF3(I)*TXX(I+1,JV)  
            FY (I,JV)=ALF (I)*(TY(I+1,JV)-ALF1(I)*TXY(I+1,JV))  
            FZ (I,JV)=ALF (I)*(TZ(I+1,JV)-ALF1(I)*TXZ(I+1,JV))  
            FXY(I,JV)=ALFQ(I)*TXY(I+1,JV)  
            FXZ(I,JV)=ALFQ(I)*TXZ(I+1,JV)  
            FYY(I,JV)=ALF (I)*TYY(I+1,JV)  
            FYZ(I,JV)=ALF (I)*TYZ(I+1,JV)  
            FZZ(I,JV)=ALF (I)*TZZ(I+1,JV)  
 C  
            T0 (I+1,JV)=T0(I+1,JV)-F0(I,JV)  
            TX (I+1,JV)=ALF1Q(I)*(TX(I+1,JV)+3.*ALF(I)*TXX(I+1,JV))  
            TXX(I+1,JV)=ALF4(I)*TXX(I+1,JV)  
            TY (I+1,JV)=TY (I+1,JV)-FY (I,JV)  
            TZ (I+1,JV)=TZ (I+1,JV)-FZ (I,JV)  
            TYY(I+1,JV)=TYY(I+1,JV)-FYY(I,JV)  
            TYZ(I+1,JV)=TYZ(I+1,JV)-FYZ(I,JV)  
            TZZ(I+1,JV)=TZZ(I+1,JV)-FZZ(I,JV)  
            TXY(I+1,JV)=ALF1Q(I)*TXY(I+1,JV)  
            TXZ(I+1,JV)=ALF1Q(I)*TXZ(I+1,JV)  
 C  
          ENDIF  
 C  
  1500 CONTINUE  
  150  CONTINUE  
 C  
       I=LONK  
       IF(UEXT(I).LT.0.) THEN  
 C  
         DO 151 JV=1,NTRA  
 C  
            F0 (I,JV)=ALF (I)* ( T0(1,JV)-ALF1(I)*  
      +             ( TX(1,JV)-ALF2(I)*TXX(1,JV) ) )  
            FX (I,JV)=ALFQ(I)*(TX(1,JV)-3.*ALF1(I)*TXX(1,JV))  
            FXX(I,JV)=ALF3(I)*TXX(1,JV)  
            FY (I,JV)=ALF (I)*(TY(1,JV)-ALF1(I)*TXY(1,JV))  
            FZ (I,JV)=ALF (I)*(TZ(1,JV)-ALF1(I)*TXZ(1,JV))  
            FXY(I,JV)=ALFQ(I)*TXY(1,JV)  
            FXZ(I,JV)=ALFQ(I)*TXZ(1,JV)  
            FYY(I,JV)=ALF (I)*TYY(1,JV)  
            FYZ(I,JV)=ALF (I)*TYZ(1,JV)  
            FZZ(I,JV)=ALF (I)*TZZ(1,JV)  
 C  
            T0 (1,JV)=T0(1,JV)-F0(I,JV)  
            TX (1,JV)=ALF1Q(I)*(TX(1,JV)+3.*ALF(I)*TXX(1,JV))  
            TXX(1,JV)=ALF4(I)*TXX(1,JV)  
            TY (1,JV)=TY (1,JV)-FY (I,JV)  
            TZ (1,JV)=TZ (1,JV)-FZ (I,JV)  
            TYY(1,JV)=TYY(1,JV)-FYY(I,JV)  
            TYZ(1,JV)=TYZ(1,JV)-FYZ(I,JV)  
            TZZ(1,JV)=TZZ(1,JV)-FZZ(I,JV)  
            TXY(1,JV)=ALF1Q(I)*TXY(1,JV)  
            TXZ(1,JV)=ALF1Q(I)*TXZ(1,JV)  
 C  
  151    CONTINUE  
 C  
       ENDIF  
 C  
       DO 152 JV=1,NTRA  
       DO 1520 I=1,LONK  
 C  
          IF(UEXT(I).GE.0.) THEN  
 C  
            F0 (I,JV)=ALF (I)* ( T0(I,JV)+ALF1(I)*  
      +             ( TX(I,JV)+ALF2(I)*TXX(I,JV) ) )  
            FX (I,JV)=ALFQ(I)*(TX(I,JV)+3.*ALF1(I)*TXX(I,JV))  
            FXX(I,JV)=ALF3(I)*TXX(I,JV)  
            FY (I,JV)=ALF (I)*(TY(I,JV)+ALF1(I)*TXY(I,JV))  
            FZ (I,JV)=ALF (I)*(TZ(I,JV)+ALF1(I)*TXZ(I,JV))  
            FXY(I,JV)=ALFQ(I)*TXY(I,JV)  
            FXZ(I,JV)=ALFQ(I)*TXZ(I,JV)  
            FYY(I,JV)=ALF (I)*TYY(I,JV)  
            FYZ(I,JV)=ALF (I)*TYZ(I,JV)  
            FZZ(I,JV)=ALF (I)*TZZ(I,JV)  
 C  
            T0 (I,JV)=T0(I,JV)-F0(I,JV)  
            TX (I,JV)=ALF1Q(I)*(TX(I,JV)-3.*ALF(I)*TXX(I,JV))  
            TXX(I,JV)=ALF4(I)*TXX(I,JV)  
            TY (I,JV)=TY (I,JV)-FY (I,JV)  
            TZ (I,JV)=TZ (I,JV)-FZ (I,JV)  
            TYY(I,JV)=TYY(I,JV)-FYY(I,JV)  
            TYZ(I,JV)=TYZ(I,JV)-FYZ(I,JV)  
            TZZ(I,JV)=TZZ(I,JV)-FZZ(I,JV)  
            TXY(I,JV)=ALF1Q(I)*TXY(I,JV)  
            TXZ(I,JV)=ALF1Q(I)*TXZ(I,JV)  
 C  
          ENDIF  
 C  
  1520 CONTINUE  
  152  CONTINUE  
 C  
 C  puts the temporary moments Fi into appropriate neighboring boxes  
 C  
       DO 160 I=1,LONK  
          IF(UEXT(I).LT.0.) THEN  
            TM(I)=TM(I)+FM(I)  
            ALF(I)=FM(I)/TM(I)  
          ENDIF  
  160  CONTINUE  
 C  
       DO 161 I=1,LONK-1  
          IF(UEXT(I).GE.0.) THEN  
            TM(I+1)=TM(I+1)+FM(I)  
            ALF(I)=FM(I)/TM(I+1)  
          ENDIF  
  161  CONTINUE  
 C  
       I=LONK  
       IF(UEXT(I).GE.0.) THEN  
         TM(1)=TM(1)+FM(I)  
         ALF(I)=FM(I)/TM(1)  
       ENDIF  
 C  
       DO 162 I=1,LONK  
          ALF1(I)=1.-ALF(I)  
          ALFQ(I)=ALF(I)*ALF(I)  
          ALF1Q(I)=ALF1(I)*ALF1(I)  
          ALF2(I)=ALF1(I)-ALF(I)  
          ALF3(I)=ALF(I)*ALF1(I)  
  162  CONTINUE  
 C  
       DO 170 JV=1,NTRA  
       DO 1700 I=1,LONK  
 C  
          IF(UEXT(I).LT.0.) THEN  
 C  
            TEMPTM=-ALF(I)*T0(I,JV)+ALF1(I)*F0(I,JV)  
            T0 (I,JV)=T0(I,JV)+F0(I,JV)  
            TXX(I,JV)=ALFQ(I)*FXX(I,JV)+ALF1Q(I)*TXX(I,JV)  
      +          +5.*( ALF3(I)*(FX(I,JV)-TX(I,JV))+ALF2(I)*TEMPTM )  
            TX (I,JV)=ALF (I)*FX (I,JV)+ALF1(I)*TX (I,JV)+3.*TEMPTM  
            TXY(I,JV)=ALF (I)*FXY(I,JV)+ALF1(I)*TXY(I,JV)  
      +          +3.*(ALF1(I)*FY (I,JV)-ALF (I)*TY (I,JV))  
            TXZ(I,JV)=ALF (I)*FXZ(I,JV)+ALF1(I)*TXZ(I,JV)  
      +          +3.*(ALF1(I)*FZ (I,JV)-ALF (I)*TZ (I,JV))  
            TY (I,JV)=TY (I,JV)+FY (I,JV)  
            TZ (I,JV)=TZ (I,JV)+FZ (I,JV)  
            TYY(I,JV)=TYY(I,JV)+FYY(I,JV)  
            TYZ(I,JV)=TYZ(I,JV)+FYZ(I,JV)  
            TZZ(I,JV)=TZZ(I,JV)+FZZ(I,JV)  
 C  
          ENDIF  
 C  
  1700 CONTINUE  
  170  CONTINUE  
 C  
       DO 171 JV=1,NTRA  
       DO 1710 I=1,LONK-1  
 C  
          IF(UEXT(I).GE.0.) THEN  
 C  
            TEMPTM=ALF(I)*T0(I+1,JV)-ALF1(I)*F0(I,JV)  
            T0 (I+1,JV)=T0(I+1,JV)+F0(I,JV)  
            TXX(I+1,JV)=ALFQ(I)*FXX(I,JV)+ALF1Q(I)*TXX(I+1,JV)  
      +           +5.*( ALF3(I)*(TX(I+1,JV)-FX(I,JV))-ALF2(I)*TEMPTM )  
            TX (I+1,JV)=ALF(I)*FX (I  ,JV)+ALF1(I)*TX (I+1,JV)+3.*TEMPTM  
            TXY(I+1,JV)=ALF(I)*FXY(I  ,JV)+ALF1(I)*TXY(I+1,JV)  
      +            +3.*(ALF(I)*TY (I+1,JV)-ALF1(I)*FY (I  ,JV))  
            TXZ(I+1,JV)=ALF(I)*FXZ(I  ,JV)+ALF1(I)*TXZ(I+1,JV)  
      +            +3.*(ALF(I)*TZ (I+1,JV)-ALF1(I)*FZ (I  ,JV))  
            TY (I+1,JV)=TY (I+1,JV)+FY (I,JV)  
            TZ (I+1,JV)=TZ (I+1,JV)+FZ (I,JV)  
            TYY(I+1,JV)=TYY(I+1,JV)+FYY(I,JV)  
            TYZ(I+1,JV)=TYZ(I+1,JV)+FYZ(I,JV)  
            TZZ(I+1,JV)=TZZ(I+1,JV)+FZZ(I,JV)  
 C  
          ENDIF  
 C  
  1710 CONTINUE  
  171  CONTINUE  
 C  
       I=LONK  
       IF(UEXT(I).GE.0.) THEN  
         DO 172 JV=1,NTRA  
            TEMPTM=ALF(I)*T0(1,JV)-ALF1(I)*F0(I,JV)  
            T0 (1,JV)=T0(1,JV)+F0(I,JV)  
            TXX(1,JV)=ALFQ(I)*FXX(I,JV)+ALF1Q(I)*TXX(1,JV)  
      +         +5.*( ALF3(I)*(TX(1,JV)-FX(I,JV))-ALF2(I)*TEMPTM )  
            TX (1,JV)=ALF(I)*FX(I,JV)+ALF1(I)*TX(1,JV)+3.*TEMPTM  
            TXY(1,JV)=ALF(I)*FXY(I,JV)+ALF1(I)*TXY(1,JV)  
      +          +3.*(ALF(I)*TY (1,JV)-ALF1(I)*FY (I,JV))  
            TXZ(1,JV)=ALF(I)*FXZ(I,JV)+ALF1(I)*TXZ(1,JV)  
      +          +3.*(ALF(I)*TZ (1,JV)-ALF1(I)*FZ (I,JV))  
            TY (1,JV)=TY (1,JV)+FY (I,JV)  
            TZ (1,JV)=TZ (1,JV)+FZ (I,JV)  
            TYY(1,JV)=TYY(1,JV)+FYY(I,JV)  
            TYZ(1,JV)=TYZ(1,JV)+FYZ(I,JV)  
            TZZ(1,JV)=TZZ(1,JV)+FZZ(I,JV)  
  172    CONTINUE  
       ENDIF  
 C  
 C  retour aux mailles d'origine (passage des Tij aux Sij)  
 C  
       IF(NUMK.GT.1) THEN  
 C  
       DO 18 I2=1,NUMK  
 C  
          DO 180 I=1,LONK  
 C  
             I3=I2+(I-1)*NUMK  
             SM(I3,K,L)=SMNEW(I3)  
             ALF(I)=SMNEW(I3)/TM(I)  
             TM(I)=TM(I)-SMNEW(I3)  
 C  
             ALFQ(I)=ALF(I)*ALF(I)  
             ALF1(I)=1.-ALF(I)  
             ALF1Q(I)=ALF1(I)*ALF1(I)  
             ALF2(I)=ALF1(I)-ALF(I)  
             ALF3(I)=ALF(I)*ALFQ(I)  
             ALF4(I)=ALF1(I)*ALF1Q(I)  
 C  
  180     CONTINUE  
 C  
          DO 181 JV=1,NTRA  
          DO 181 I=1,LONK  
 C  
             I3=I2+(I-1)*NUMK  
             S0 (I3,K,L,JV)=ALF (I)* ( T0(I,JV)-ALF1(I)*  
      +              ( TX(I,JV)-ALF2(I)*TXX(I,JV) ) )  
             SSX (I3,K,L,JV)=ALFQ(I)*(TX(I,JV)-3.*ALF1(I)*TXX(I,JV))  
             SSXX(I3,K,L,JV)=ALF3(I)*TXX(I,JV)  
             SY (I3,K,L,JV)=ALF (I)*(TY(I,JV)-ALF1(I)*TXY(I,JV))  
             SZ (I3,K,L,JV)=ALF (I)*(TZ(I,JV)-ALF1(I)*TXZ(I,JV))  
             SSXY(I3,K,L,JV)=ALFQ(I)*TXY(I,JV)  
             SSXZ(I3,K,L,JV)=ALFQ(I)*TXZ(I,JV)  
             SYY(I3,K,L,JV)=ALF (I)*TYY(I,JV)  
             SYZ(I3,K,L,JV)=ALF (I)*TYZ(I,JV)  
             SZZ(I3,K,L,JV)=ALF (I)*TZZ(I,JV)  
 C  
 C   reajusts moments remaining in the box  
 C  
             T0 (I,JV)=T0(I,JV)-S0(I3,K,L,JV)  
             TX (I,JV)=ALF1Q(I)*(TX(I,JV)+3.*ALF(I)*TXX(I,JV))  
             TXX(I,JV)=ALF4 (I)*TXX(I,JV)  
             TY (I,JV)=TY (I,JV)-SY (I3,K,L,JV)  
             TZ (I,JV)=TZ (I,JV)-SZ (I3,K,L,JV)  
             TYY(I,JV)=TYY(I,JV)-SYY(I3,K,L,JV)  
             TYZ(I,JV)=TYZ(I,JV)-SYZ(I3,K,L,JV)  
             TZZ(I,JV)=TZZ(I,JV)-SZZ(I3,K,L,JV)  
             TXY(I,JV)=ALF1Q(I)*TXY(I,JV)  
             TXZ(I,JV)=ALF1Q(I)*TXZ(I,JV)  
 C  
  181     CONTINUE  
 C  
  18   CONTINUE  
 C  
222        ELSE        ELSE
223  C  
224        DO 190 I=1,LON          DO i = 1, lon
225           SM(I,K,L)=TM(I)            tm(i) = sm(i, k, l)
226   190  CONTINUE          END DO
227        DO 191 JV=1,NTRA          DO jv = 1, ntra
228        DO 1910 I=1,LON            DO i = 1, lon
229           S0 (I,K,L,JV)=T0 (I,JV)              t0(i, jv) = s0(i, k, l, jv)
230           SSX (I,K,L,JV)=TX (I,JV)              tx(i, jv) = ssx(i, k, l, jv)
231           SY (I,K,L,JV)=TY (I,JV)              ty(i, jv) = sy(i, k, l, jv)
232           SZ (I,K,L,JV)=TZ (I,JV)              tz(i, jv) = sz(i, k, l, jv)
233           SSXX(I,K,L,JV)=TXX(I,JV)              txx(i, jv) = ssxx(i, k, l, jv)
234           SSXY(I,K,L,JV)=TXY(I,JV)              txy(i, jv) = ssxy(i, k, l, jv)
235           SSXZ(I,K,L,JV)=TXZ(I,JV)              txz(i, jv) = ssxz(i, k, l, jv)
236           SYY(I,K,L,JV)=TYY(I,JV)              tyy(i, jv) = syy(i, k, l, jv)
237           SYZ(I,K,L,JV)=TYZ(I,JV)              tyz(i, jv) = syz(i, k, l, jv)
238           SZZ(I,K,L,JV)=TZZ(I,JV)              tzz(i, jv) = szz(i, k, l, jv)
239   1910 CONTINUE            END DO
240   191  CONTINUE          END DO
241  C  
242        ENDIF        END IF
243  C  
244   1    CONTINUE        DO i = 1, lonk
245            uext(i) = ugri(i*numk, k, l)
 c ---------- bouclage cyclique  
   
       DO l = 1,llm  
       DO j = 1,jjp1  
          SM(iip1,j,l) = SM(1,j,l)  
          S0(iip1,j,l,ntra) = S0(1,j,l,ntra)  
          SSX(iip1,j,l,ntra) = SSX(1,j,l,ntra)  
          SY(iip1,j,l,ntra) = SY(1,j,l,ntra)  
          SZ(iip1,j,l,ntra) = SZ(1,j,l,ntra)  
246        END DO        END DO
247    
248          ! place limits on appropriate moments before transport
249          ! (if flux-limiting is to be applied)
250    
251          IF (.NOT. limit) GO TO 13
252    
253          DO jv = 1, ntra
254            DO i = 1, lonk
255              IF (t0(i,jv)>0.) THEN
256                slpmax = t0(i, jv)
257                s1max = 1.5*slpmax
258                s1new = amin1(s1max, amax1(-s1max,tx(i,jv)))
259                s2new = amin1(2.*slpmax-abs(s1new)/3., amax1(abs( &
260                  s1new)-slpmax,txx(i,jv)))
261                tx(i, jv) = s1new
262                txx(i, jv) = s2new
263                txy(i, jv) = amin1(slpmax, amax1(-slpmax,txy(i,jv)))
264                txz(i, jv) = amin1(slpmax, amax1(-slpmax,txz(i,jv)))
265              ELSE
266                tx(i, jv) = 0.
267                txx(i, jv) = 0.
268                txy(i, jv) = 0.
269                txz(i, jv) = 0.
270              END IF
271            END DO
272        END DO        END DO
273    
274  C ----------- qqtite totale de traceur dans tte l'atmosphere  13    CONTINUE
275        DO l = 1, llm  
276        DO j = 1, jjp1        ! calculate flux and moments between adjacent boxes
277        DO i = 1, iim        ! 1- create temporary moments/masses for partial boxes in transit
278          sqf = sqf + S0(i,j,l,ntra)        ! 2- reajusts moments remaining in the box
279    
280          ! flux from IP to I if U(I).lt.0
281    
282          DO i = 1, lonk - 1
283            IF (uext(i)<0.) THEN
284              fm(i) = -uext(i)*dtx
285              alf(i) = fm(i)/tm(i+1)
286              tm(i+1) = tm(i+1) - fm(i)
287            END IF
288          END DO
289    
290          i = lonk
291          IF (uext(i)<0.) THEN
292            fm(i) = -uext(i)*dtx
293            alf(i) = fm(i)/tm(1)
294            tm(1) = tm(1) - fm(i)
295          END IF
296    
297          ! flux from I to IP if U(I).gt.0
298    
299          DO i = 1, lonk
300            IF (uext(i)>=0.) THEN
301              fm(i) = uext(i)*dtx
302              alf(i) = fm(i)/tm(i)
303              tm(i) = tm(i) - fm(i)
304            END IF
305          END DO
306    
307          DO i = 1, lonk
308            alfq(i) = alf(i)*alf(i)
309            alf1(i) = 1. - alf(i)
310            alf1q(i) = alf1(i)*alf1(i)
311            alf2(i) = alf1(i) - alf(i)
312            alf3(i) = alf(i)*alfq(i)
313            alf4(i) = alf1(i)*alf1q(i)
314          END DO
315    
316          DO jv = 1, ntra
317            DO i = 1, lonk - 1
318    
319              IF (uext(i)<0.) THEN
320    
321                f0(i, jv) = alf(i)*(t0(i+1,jv)-alf1(i)*(tx(i+1, &
322                  jv)-alf2(i)*txx(i+1,jv)))
323                fx(i, jv) = alfq(i)*(tx(i+1,jv)-3.*alf1(i)*txx(i+1,jv))
324                fxx(i, jv) = alf3(i)*txx(i+1, jv)
325                fy(i, jv) = alf(i)*(ty(i+1,jv)-alf1(i)*txy(i+1,jv))
326                fz(i, jv) = alf(i)*(tz(i+1,jv)-alf1(i)*txz(i+1,jv))
327                fxy(i, jv) = alfq(i)*txy(i+1, jv)
328                fxz(i, jv) = alfq(i)*txz(i+1, jv)
329                fyy(i, jv) = alf(i)*tyy(i+1, jv)
330                fyz(i, jv) = alf(i)*tyz(i+1, jv)
331                fzz(i, jv) = alf(i)*tzz(i+1, jv)
332    
333                t0(i+1, jv) = t0(i+1, jv) - f0(i, jv)
334                tx(i+1, jv) = alf1q(i)*(tx(i+1,jv)+3.*alf(i)*txx(i+1,jv))
335                txx(i+1, jv) = alf4(i)*txx(i+1, jv)
336                ty(i+1, jv) = ty(i+1, jv) - fy(i, jv)
337                tz(i+1, jv) = tz(i+1, jv) - fz(i, jv)
338                tyy(i+1, jv) = tyy(i+1, jv) - fyy(i, jv)
339                tyz(i+1, jv) = tyz(i+1, jv) - fyz(i, jv)
340                tzz(i+1, jv) = tzz(i+1, jv) - fzz(i, jv)
341                txy(i+1, jv) = alf1q(i)*txy(i+1, jv)
342                txz(i+1, jv) = alf1q(i)*txz(i+1, jv)
343    
344              END IF
345    
346            END DO
347          END DO
348    
349          i = lonk
350          IF (uext(i)<0.) THEN
351    
352            DO jv = 1, ntra
353    
354              f0(i, jv) = alf(i)*(t0(1,jv)-alf1(i)*(tx(1,jv)-alf2(i)*txx(1,jv)))
355              fx(i, jv) = alfq(i)*(tx(1,jv)-3.*alf1(i)*txx(1,jv))
356              fxx(i, jv) = alf3(i)*txx(1, jv)
357              fy(i, jv) = alf(i)*(ty(1,jv)-alf1(i)*txy(1,jv))
358              fz(i, jv) = alf(i)*(tz(1,jv)-alf1(i)*txz(1,jv))
359              fxy(i, jv) = alfq(i)*txy(1, jv)
360              fxz(i, jv) = alfq(i)*txz(1, jv)
361              fyy(i, jv) = alf(i)*tyy(1, jv)
362              fyz(i, jv) = alf(i)*tyz(1, jv)
363              fzz(i, jv) = alf(i)*tzz(1, jv)
364    
365              t0(1, jv) = t0(1, jv) - f0(i, jv)
366              tx(1, jv) = alf1q(i)*(tx(1,jv)+3.*alf(i)*txx(1,jv))
367              txx(1, jv) = alf4(i)*txx(1, jv)
368              ty(1, jv) = ty(1, jv) - fy(i, jv)
369              tz(1, jv) = tz(1, jv) - fz(i, jv)
370              tyy(1, jv) = tyy(1, jv) - fyy(i, jv)
371              tyz(1, jv) = tyz(1, jv) - fyz(i, jv)
372              tzz(1, jv) = tzz(1, jv) - fzz(i, jv)
373              txy(1, jv) = alf1q(i)*txy(1, jv)
374              txz(1, jv) = alf1q(i)*txz(1, jv)
375    
376            END DO
377    
378          END IF
379    
380          DO jv = 1, ntra
381            DO i = 1, lonk
382    
383              IF (uext(i)>=0.) THEN
384    
385                f0(i, jv) = alf(i)*(t0(i,jv)+alf1(i)*(tx(i,jv)+alf2(i)*txx(i, &
386                  jv)))
387                fx(i, jv) = alfq(i)*(tx(i,jv)+3.*alf1(i)*txx(i,jv))
388                fxx(i, jv) = alf3(i)*txx(i, jv)
389                fy(i, jv) = alf(i)*(ty(i,jv)+alf1(i)*txy(i,jv))
390                fz(i, jv) = alf(i)*(tz(i,jv)+alf1(i)*txz(i,jv))
391                fxy(i, jv) = alfq(i)*txy(i, jv)
392                fxz(i, jv) = alfq(i)*txz(i, jv)
393                fyy(i, jv) = alf(i)*tyy(i, jv)
394                fyz(i, jv) = alf(i)*tyz(i, jv)
395                fzz(i, jv) = alf(i)*tzz(i, jv)
396    
397                t0(i, jv) = t0(i, jv) - f0(i, jv)
398                tx(i, jv) = alf1q(i)*(tx(i,jv)-3.*alf(i)*txx(i,jv))
399                txx(i, jv) = alf4(i)*txx(i, jv)
400                ty(i, jv) = ty(i, jv) - fy(i, jv)
401                tz(i, jv) = tz(i, jv) - fz(i, jv)
402                tyy(i, jv) = tyy(i, jv) - fyy(i, jv)
403                tyz(i, jv) = tyz(i, jv) - fyz(i, jv)
404                tzz(i, jv) = tzz(i, jv) - fzz(i, jv)
405                txy(i, jv) = alf1q(i)*txy(i, jv)
406                txz(i, jv) = alf1q(i)*txz(i, jv)
407    
408              END IF
409    
410            END DO
411          END DO
412    
413          ! puts the temporary moments Fi into appropriate neighboring boxes
414    
415          DO i = 1, lonk
416            IF (uext(i)<0.) THEN
417              tm(i) = tm(i) + fm(i)
418              alf(i) = fm(i)/tm(i)
419            END IF
420          END DO
421    
422          DO i = 1, lonk - 1
423            IF (uext(i)>=0.) THEN
424              tm(i+1) = tm(i+1) + fm(i)
425              alf(i) = fm(i)/tm(i+1)
426            END IF
427          END DO
428    
429          i = lonk
430          IF (uext(i)>=0.) THEN
431            tm(1) = tm(1) + fm(i)
432            alf(i) = fm(i)/tm(1)
433          END IF
434    
435          DO i = 1, lonk
436            alf1(i) = 1. - alf(i)
437            alfq(i) = alf(i)*alf(i)
438            alf1q(i) = alf1(i)*alf1(i)
439            alf2(i) = alf1(i) - alf(i)
440            alf3(i) = alf(i)*alf1(i)
441          END DO
442    
443          DO jv = 1, ntra
444            DO i = 1, lonk
445    
446              IF (uext(i)<0.) THEN
447    
448                temptm = -alf(i)*t0(i, jv) + alf1(i)*f0(i, jv)
449                t0(i, jv) = t0(i, jv) + f0(i, jv)
450                txx(i, jv) = alfq(i)*fxx(i, jv) + alf1q(i)*txx(i, jv) + &
451                  5.*(alf3(i)*(fx(i,jv)-tx(i,jv))+alf2(i)*temptm)
452                tx(i, jv) = alf(i)*fx(i, jv) + alf1(i)*tx(i, jv) + 3.*temptm
453                txy(i, jv) = alf(i)*fxy(i, jv) + alf1(i)*txy(i, jv) + &
454                  3.*(alf1(i)*fy(i,jv)-alf(i)*ty(i,jv))
455                txz(i, jv) = alf(i)*fxz(i, jv) + alf1(i)*txz(i, jv) + &
456                  3.*(alf1(i)*fz(i,jv)-alf(i)*tz(i,jv))
457                ty(i, jv) = ty(i, jv) + fy(i, jv)
458                tz(i, jv) = tz(i, jv) + fz(i, jv)
459                tyy(i, jv) = tyy(i, jv) + fyy(i, jv)
460                tyz(i, jv) = tyz(i, jv) + fyz(i, jv)
461                tzz(i, jv) = tzz(i, jv) + fzz(i, jv)
462    
463              END IF
464    
465            END DO
466        END DO        END DO
467    
468          DO jv = 1, ntra
469            DO i = 1, lonk - 1
470    
471              IF (uext(i)>=0.) THEN
472    
473                temptm = alf(i)*t0(i+1, jv) - alf1(i)*f0(i, jv)
474                t0(i+1, jv) = t0(i+1, jv) + f0(i, jv)
475                txx(i+1, jv) = alfq(i)*fxx(i, jv) + alf1q(i)*txx(i+1, jv) + &
476                  5.*(alf3(i)*(tx(i+1,jv)-fx(i,jv))-alf2(i)*temptm)
477                tx(i+1, jv) = alf(i)*fx(i, jv) + alf1(i)*tx(i+1, jv) + 3.*temptm
478                txy(i+1, jv) = alf(i)*fxy(i, jv) + alf1(i)*txy(i+1, jv) + &
479                  3.*(alf(i)*ty(i+1,jv)-alf1(i)*fy(i,jv))
480                txz(i+1, jv) = alf(i)*fxz(i, jv) + alf1(i)*txz(i+1, jv) + &
481                  3.*(alf(i)*tz(i+1,jv)-alf1(i)*fz(i,jv))
482                ty(i+1, jv) = ty(i+1, jv) + fy(i, jv)
483                tz(i+1, jv) = tz(i+1, jv) + fz(i, jv)
484                tyy(i+1, jv) = tyy(i+1, jv) + fyy(i, jv)
485                tyz(i+1, jv) = tyz(i+1, jv) + fyz(i, jv)
486                tzz(i+1, jv) = tzz(i+1, jv) + fzz(i, jv)
487    
488              END IF
489    
490            END DO
491        END DO        END DO
492    
493          i = lonk
494          IF (uext(i)>=0.) THEN
495            DO jv = 1, ntra
496              temptm = alf(i)*t0(1, jv) - alf1(i)*f0(i, jv)
497              t0(1, jv) = t0(1, jv) + f0(i, jv)
498              txx(1, jv) = alfq(i)*fxx(i, jv) + alf1q(i)*txx(1, jv) + &
499                5.*(alf3(i)*(tx(1,jv)-fx(i,jv))-alf2(i)*temptm)
500              tx(1, jv) = alf(i)*fx(i, jv) + alf1(i)*tx(1, jv) + 3.*temptm
501              txy(1, jv) = alf(i)*fxy(i, jv) + alf1(i)*txy(1, jv) + &
502                3.*(alf(i)*ty(1,jv)-alf1(i)*fy(i,jv))
503              txz(1, jv) = alf(i)*fxz(i, jv) + alf1(i)*txz(1, jv) + &
504                3.*(alf(i)*tz(1,jv)-alf1(i)*fz(i,jv))
505              ty(1, jv) = ty(1, jv) + fy(i, jv)
506              tz(1, jv) = tz(1, jv) + fz(i, jv)
507              tyy(1, jv) = tyy(1, jv) + fyy(i, jv)
508              tyz(1, jv) = tyz(1, jv) + fyz(i, jv)
509              tzz(1, jv) = tzz(1, jv) + fzz(i, jv)
510            END DO
511          END IF
512    
513          ! retour aux mailles d'origine (passage des Tij aux Sij)
514    
515          IF (numk>1) THEN
516    
517            DO i2 = 1, numk
518    
519              DO i = 1, lonk
520    
521                i3 = i2 + (i-1)*numk
522                sm(i3, k, l) = smnew(i3)
523                alf(i) = smnew(i3)/tm(i)
524                tm(i) = tm(i) - smnew(i3)
525    
526                alfq(i) = alf(i)*alf(i)
527                alf1(i) = 1. - alf(i)
528                alf1q(i) = alf1(i)*alf1(i)
529                alf2(i) = alf1(i) - alf(i)
530                alf3(i) = alf(i)*alfq(i)
531                alf4(i) = alf1(i)*alf1q(i)
532    
533              END DO
534    
535              DO jv = 1, ntra
536                DO i = 1, lonk
537    
538                  i3 = i2 + (i-1)*numk
539                  s0(i3, k, l, jv) = alf(i)*(t0(i,jv)-alf1(i)*(tx(i, &
540                    jv)-alf2(i)*txx(i,jv)))
541                  ssx(i3, k, l, jv) = alfq(i)*(tx(i,jv)-3.*alf1(i)*txx(i,jv))
542                  ssxx(i3, k, l, jv) = alf3(i)*txx(i, jv)
543                  sy(i3, k, l, jv) = alf(i)*(ty(i,jv)-alf1(i)*txy(i,jv))
544                  sz(i3, k, l, jv) = alf(i)*(tz(i,jv)-alf1(i)*txz(i,jv))
545                  ssxy(i3, k, l, jv) = alfq(i)*txy(i, jv)
546                  ssxz(i3, k, l, jv) = alfq(i)*txz(i, jv)
547                  syy(i3, k, l, jv) = alf(i)*tyy(i, jv)
548                  syz(i3, k, l, jv) = alf(i)*tyz(i, jv)
549                  szz(i3, k, l, jv) = alf(i)*tzz(i, jv)
550    
551                  ! reajusts moments remaining in the box
552    
553                  t0(i, jv) = t0(i, jv) - s0(i3, k, l, jv)
554                  tx(i, jv) = alf1q(i)*(tx(i,jv)+3.*alf(i)*txx(i,jv))
555                  txx(i, jv) = alf4(i)*txx(i, jv)
556                  ty(i, jv) = ty(i, jv) - sy(i3, k, l, jv)
557                  tz(i, jv) = tz(i, jv) - sz(i3, k, l, jv)
558                  tyy(i, jv) = tyy(i, jv) - syy(i3, k, l, jv)
559                  tyz(i, jv) = tyz(i, jv) - syz(i3, k, l, jv)
560                  tzz(i, jv) = tzz(i, jv) - szz(i3, k, l, jv)
561                  txy(i, jv) = alf1q(i)*txy(i, jv)
562                  txz(i, jv) = alf1q(i)*txz(i, jv)
563    
564                END DO
565              END DO
566    
567            END DO
568    
569          ELSE
570    
571            DO i = 1, lon
572              sm(i, k, l) = tm(i)
573            END DO
574            DO jv = 1, ntra
575              DO i = 1, lon
576                s0(i, k, l, jv) = t0(i, jv)
577                ssx(i, k, l, jv) = tx(i, jv)
578                sy(i, k, l, jv) = ty(i, jv)
579                sz(i, k, l, jv) = tz(i, jv)
580                ssxx(i, k, l, jv) = txx(i, jv)
581                ssxy(i, k, l, jv) = txy(i, jv)
582                ssxz(i, k, l, jv) = txz(i, jv)
583                syy(i, k, l, jv) = tyy(i, jv)
584                syz(i, k, l, jv) = tyz(i, jv)
585                szz(i, k, l, jv) = tzz(i, jv)
586              END DO
587            END DO
588    
589          END IF
590    
591        END DO
592      END DO
593    
594      ! ---------- bouclage cyclique
595    
596      DO l = 1, llm
597        DO j = 1, jjp1
598          sm(iip1, j, l) = sm(1, j, l)
599          s0(iip1, j, l, ntra) = s0(1, j, l, ntra)
600          ssx(iip1, j, l, ntra) = ssx(1, j, l, ntra)
601          sy(iip1, j, l, ntra) = sy(1, j, l, ntra)
602          sz(iip1, j, l, ntra) = sz(1, j, l, ntra)
603        END DO
604      END DO
605    
606      ! ----------- qqtite totale de traceur dans tte l'atmosphere
607      DO l = 1, llm
608        DO j = 1, jjp1
609          DO i = 1, iim
610            sqf = sqf + s0(i, j, l, ntra)
611        END DO        END DO
612        END DO
613      END DO
614    
615        PRINT*,'------ DIAG DANS ADVX2 - SORTIE -----'    PRINT *, '------ DIAG DANS ADVX2 - SORTIE -----'
616        PRINT*,'sqf=',sqf    PRINT *, 'sqf=', sqf
617  c-------------------------------------------------------------    ! -------------------------------------------------------------
618        RETURN    RETURN
619        END  END SUBROUTINE advxp
Legend:
Removed from v.76  
changed lines
  Added in v.81
  ViewVC Help
Powered by ViewVC 1.1.21