Sources/dyn3d/advx.f

!
! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/advx.F,v 1.2 2005/05/25 13:10:09 fairhead Exp $
!
      SUBROUTINE  advx(limit,dtx,pbaru,sm,s0,
     $     sx,sy,sz,lati,latf)
      use dimens_m
      use paramet_m
      use comconst
      use comvert
      IMPLICIT NONE

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                C
C  first-order moments (FOM) advection of tracer in X direction  C
C                                                                C
C  Source : Pascal Simon (Meteo,CNRM)                            C
C  Adaptation : A.Armengaud (LGGE) juin 94                       C
C                                                                C
C  limit,dtx,pbaru,pbarv,sm,s0,sx,sy,sz                       C
C  sont des arguments d'entree pour le s-pg...                   C
C                                                                C
C  sm,s0,sx,sy,sz                                                C
C  sont les arguments de sortie pour le s-pg                     C
C                                                                C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C  parametres principaux du modele
C

C  Arguments :
C  -----------
C  dtx : frequence fictive d'appel du transport 
C  pbaru, pbarv : flux de masse en x et y en Pa.m2.s-1

       INTEGER ntra
       PARAMETER (ntra = 1)

C ATTENTION partout ou on trouve ntra, insertion de boucle
C           possible dans l'avenir.

      REAL dtx
      REAL, intent(in):: pbaru ( iip1,jjp1,llm )

C  moments: SM  total mass in each grid box
C           S0  mass of tracer in each grid box
C           Si  1rst order moment in i direction
C
      REAL SM(iip1,jjp1,llm),S0(iip1,jjp1,llm,ntra)
      REAL sx(iip1,jjp1,llm,ntra)
     $    ,sy(iip1,jjp1,llm,ntra)
      REAL sz(iip1,jjp1,llm,ntra)

C  Local :
C  ------- 

C  mass fluxes across the boundaries (UGRI,VGRI,WGRI)
C  mass fluxes in kg
C  declaration :

      REAL UGRI(iip1,jjp1,llm)

C  Rem : VGRI et WGRI ne sont pas utilises dans 
C  cette subroutine ( advection en x uniquement )
C
C  Ti are the moments for the current latitude and level
C
      REAL TM(iim)
      REAL T0(iim,ntra),TX(iim,ntra)
      REAL TY(iim,ntra),TZ(iim,ntra)
      REAL TEMPTM                ! just a temporary variable
C
C  the moments F are similarly defined and used as temporary
C  storage for portions of the grid boxes in transit
C
      REAL FM(iim)
      REAL F0(iim,ntra),FX(iim,ntra)
      REAL FY(iim,ntra),FZ(iim,ntra)
C
C  work arrays
C
      REAL ALF(iim),ALF1(iim),ALFQ(iim),ALF1Q(iim)
C
      REAL SMNEW(iim),UEXT(iim)
C
      REAL sqi,sqf

      LOGICAL LIMIT
      INTEGER NUM(jjp1),LONK,NUMK
      INTEGER lon,lati,latf,niv
      INTEGER i,i2,i3,j,jv,l,k,itrac 

      lon = iim 
      niv = llm 

C *** Test de passage d'arguments ******


C  -------------------------------------
      DO 300 j = 1,jjp1 
         NUM(j) = 1
  300 CONTINUE
      sqi = 0.
      sqf = 0.

      DO l = 1,llm
         DO j = 1,jjp1
            DO i = 1,iim
cIM 240305            sqi = sqi + S0(i,j,l,9)
               sqi = sqi + S0(i,j,l,ntra)
            ENDDO
         ENDDO
      ENDDO
      PRINT*,'-------- DIAG DANS ADVX - ENTREE ---------'
      PRINT*,'sqi=',sqi


C  Interface : adaptation nouveau modele
C  -------------------------------------
C
C  ---------------------------------------------------------
C  Conversion des flux de masses en kg/s
C  pbaru est en N/s d'ou :
C  ugri est en kg/s

      DO 500 l = 1,llm
         DO 500 j = 1,jjm+1
            DO 500 i = 1,iip1  
C            ugri (i,j,llm+1-l) = pbaru (i,j,l) * ( dsig(l) / g )
             ugri (i,j,llm+1-l) = pbaru (i,j,l)
  500 CONTINUE


C  ---------------------------------------------------------
C  ---------------------------------------------------------
C  ---------------------------------------------------------
  
C  start here          
C
C  boucle principale sur les niveaux et les latitudes
C
      DO 1 L=1,NIV
      DO 1 K=lati,latf
C
C  initialisation
C
C  program assumes periodic boundaries in X
C
      DO 10 I=2,LON
         SMNEW(I)=SM(I,K,L)+(UGRI(I-1,K,L)-UGRI(I,K,L))*DTX
 10   CONTINUE
      SMNEW(1)=SM(1,K,L)+(UGRI(LON,K,L)-UGRI(1,K,L))*DTX
C
C  modifications for extended polar zones
C
      NUMK=NUM(K)
      LONK=LON/NUMK
C
      IF(NUMK.GT.1) THEN
C
      DO 111 I=1,LON
         TM(I)=0.
 111  CONTINUE
      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.
 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)
 113     CONTINUE
C
         DO  JV=1,NTRA
         DO  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)
            TX(I,JV)=ALF(I)  *sx(I3,K,L,JV)+
     $       ALF1(I)*TX(I,JV) +3.*TEMPTM
            TY(I,JV)=TY(I,JV)+sy(I3,K,L,JV)
            TZ(I,JV)=TZ(I,JV)+sz(I3,K,L,JV)
         ENDDO 
         ENDDO
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)=sx(I,K,L,JV)
         TY(I,JV)=sy(I,K,L,JV)
         TZ(I,JV)=sz(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
        TX(I,JV)=SIGN(AMIN1(AMAX1(T0(I,JV),0.),ABS(TX(I,JV))),TX(I,JV))
 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)
 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) )
           FX(I,JV)=ALFQ(I)*TX(I+1,JV)
           FY(I,JV)=ALF (I)*TY(I+1,JV)
           FZ(I,JV)=ALF (I)*TZ(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)
           TY(I+1,JV)=TY(I+1,JV)-FY(I,JV)
           TZ(I+1,JV)=TZ(I+1,JV)-FZ(I,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) )
           FX (I,JV)=ALFQ(I)*TX(1,JV)
           FY (I,JV)=ALF (I)*TY(1,JV)
           FZ (I,JV)=ALF (I)*TZ(1,JV)
C
           T0(1,JV)=T0(1,JV)-F0(I,JV)
           TX(1,JV)=ALF1Q(I)*TX(1,JV)
           TY(1,JV)=TY(1,JV)-FY(I,JV)
           TZ(1,JV)=TZ(1,JV)-FZ(I,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) )
           FX(I,JV)=ALFQ(I)*TX(I,JV)
           FY(I,JV)=ALF (I)*TY(I,JV)
           FZ(I,JV)=ALF (I)*TZ(I,JV)
C
           T0(I,JV)=T0(I,JV)-F0(I,JV)
           TX(I,JV)=ALF1Q(I)*TX(I,JV)
           TY(I,JV)=TY(I,JV)-FY(I,JV)
           TZ(I,JV)=TZ(I,JV)-FZ(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)
 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)
           TX(I,JV)=ALF(I)*FX(I,JV)+ALF1(I)*TX(I,JV)+3.*TEMPTM
           TY(I,JV)=TY(I,JV)+FY(I,JV)
           TZ(I,JV)=TZ(I,JV)+FZ(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)
           TX(I+1,JV)=ALF(I)*FX(I,JV)+ALF1(I)*TX(I+1,JV)+3.*TEMPTM
           TY(I+1,JV)=TY(I+1,JV)+FY(I,JV)
           TZ(I+1,JV)=TZ(I+1,JV)+FZ(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)
           TX(1,JV)=ALF(I)*FX(I,JV)+ALF1(I)*TX(1,JV)+3.*TEMPTM
           TY(1,JV)=TY(1,JV)+FY(I,JV)
           TZ(1,JV)=TZ(1,JV)+FZ(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 180 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)
C
 180     CONTINUE
C
         DO  JV=1,NTRA
         DO  I=1,LONK
C
            I3=I2+(I-1)*NUMK
            S0(I3,K,L,JV)=ALF (I)
     $       * (T0(I,JV)-ALF1(I)*TX(I,JV))
            sx(I3,K,L,JV)=ALFQ(I)*TX(I,JV)
            sy(I3,K,L,JV)=ALF (I)*TY(I,JV)
            sz(I3,K,L,JV)=ALF (I)*TZ(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)
            TY(I,JV)=TY(I,JV)-sy(I3,K,L,JV)
            TZ(I,JV)=TZ(I,JV)-sz(I3,K,L,JV)
          ENDDO
          ENDDO
C
C
      ELSE
C
      DO 190 I=1,LON
         SM(I,K,L)=TM(I)
 190  CONTINUE
      DO 191 JV=1,NTRA
      DO 1910 I=1,LON
         S0(I,K,L,JV)=T0(I,JV)
         sx(I,K,L,JV)=TX(I,JV)
         sy(I,K,L,JV)=TY(I,JV)
         sz(I,K,L,JV)=TZ(I,JV)
 1910 CONTINUE
 191  CONTINUE
C
      ENDIF
C
 1    CONTINUE
C
C ---------- bouclage cyclique 
      DO itrac=1,ntra
      DO l = 1,llm
        DO j = lati,latf
           SM(iip1,j,l) = SM(1,j,l)
           S0(iip1,j,l,itrac) = S0(1,j,l,itrac)
           sx(iip1,j,l,itrac) = sx(1,j,l,itrac)
           sy(iip1,j,l,itrac) = sy(1,j,l,itrac)
           sz(iip1,j,l,itrac) = sz(1,j,l,itrac)
        END DO
      END DO
      ENDDO 

c ----------- qqtite totale de traceur dans tte l'atmosphere
      DO l = 1, llm
        DO j = 1, jjp1
          DO i = 1, iim
cIM 240405          sqf = sqf + S0(i,j,l,9)
             sqf = sqf + S0(i,j,l,ntra)
          END DO  
        END DO
      END DO
c
      PRINT*,'------ DIAG DANS ADVX - SORTIE -----'
      PRINT*,'sqf=',sqf
c-------------

      RETURN
      END
C_________________________________________________________________
C_________________________________________________________________
1	guez	3	!
2			! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/advx.F,v 1.2 2005/05/25 13:10:09 fairhead Exp $
3			!
4			SUBROUTINE advx(limit,dtx,pbaru,sm,s0,
5			$ sx,sy,sz,lati,latf)
6			use dimens_m
7			use paramet_m
8			use comconst
9			use comvert
10			IMPLICIT NONE
11
12			CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
13			C C
14			C first-order moments (FOM) advection of tracer in X direction C
15			C C
16			C Source : Pascal Simon (Meteo,CNRM) C
17			C Adaptation : A.Armengaud (LGGE) juin 94 C
18			C C
19			C limit,dtx,pbaru,pbarv,sm,s0,sx,sy,sz C
20			C sont des arguments d'entree pour le s-pg... C
21			C C
22			C sm,s0,sx,sy,sz C
23			C sont les arguments de sortie pour le s-pg C
24			C C
25			CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
26			C
27			C parametres principaux du modele
28			C
29
30			C Arguments :
31			C -----------
32			C dtx : frequence fictive d'appel du transport
33			C pbaru, pbarv : flux de masse en x et y en Pa.m2.s-1
34
35			INTEGER ntra
36			PARAMETER (ntra = 1)
37
38			C ATTENTION partout ou on trouve ntra, insertion de boucle
39			C possible dans l'avenir.
40
41			REAL dtx
42	guez	31	REAL, intent(in):: pbaru ( iip1,jjp1,llm )
43	guez	3
44			C moments: SM total mass in each grid box
45			C S0 mass of tracer in each grid box
46			C Si 1rst order moment in i direction
47			C
48			REAL SM(iip1,jjp1,llm),S0(iip1,jjp1,llm,ntra)
49			REAL sx(iip1,jjp1,llm,ntra)
50			$ ,sy(iip1,jjp1,llm,ntra)
51			REAL sz(iip1,jjp1,llm,ntra)
52
53			C Local :
54			C -------
55
56			C mass fluxes across the boundaries (UGRI,VGRI,WGRI)
57			C mass fluxes in kg
58			C declaration :
59
60			REAL UGRI(iip1,jjp1,llm)
61
62			C Rem : VGRI et WGRI ne sont pas utilises dans
63			C cette subroutine ( advection en x uniquement )
64			C
65			C Ti are the moments for the current latitude and level
66			C
67			REAL TM(iim)
68			REAL T0(iim,ntra),TX(iim,ntra)
69			REAL TY(iim,ntra),TZ(iim,ntra)
70			REAL TEMPTM ! just a temporary variable
71			C
72			C the moments F are similarly defined and used as temporary
73			C storage for portions of the grid boxes in transit
74			C
75			REAL FM(iim)
76			REAL F0(iim,ntra),FX(iim,ntra)
77			REAL FY(iim,ntra),FZ(iim,ntra)
78			C
79			C work arrays
80			C
81			REAL ALF(iim),ALF1(iim),ALFQ(iim),ALF1Q(iim)
82			C
83			REAL SMNEW(iim),UEXT(iim)
84			C
85			REAL sqi,sqf
86
87			LOGICAL LIMIT
88			INTEGER NUM(jjp1),LONK,NUMK
89			INTEGER lon,lati,latf,niv
90			INTEGER i,i2,i3,j,jv,l,k,itrac
91
92			lon = iim
93			niv = llm
94
95			C * Test de passage d'arguments ****
96
97
98			C -------------------------------------
99			DO 300 j = 1,jjp1
100			NUM(j) = 1
101			300 CONTINUE
102			sqi = 0.
103			sqf = 0.
104
105			DO l = 1,llm
106			DO j = 1,jjp1
107			DO i = 1,iim
108			cIM 240305 sqi = sqi + S0(i,j,l,9)
109			sqi = sqi + S0(i,j,l,ntra)
110			ENDDO
111			ENDDO
112			ENDDO
113			PRINT*,'-------- DIAG DANS ADVX - ENTREE ---------'
114			PRINT*,'sqi=',sqi
115
116
117			C Interface : adaptation nouveau modele
118			C -------------------------------------
119			C
120			C ---------------------------------------------------------
121			C Conversion des flux de masses en kg/s
122			C pbaru est en N/s d'ou :
123			C ugri est en kg/s
124
125			DO 500 l = 1,llm
126			DO 500 j = 1,jjm+1
127			DO 500 i = 1,iip1
128			C ugri (i,j,llm+1-l) = pbaru (i,j,l) * ( dsig(l) / g )
129			ugri (i,j,llm+1-l) = pbaru (i,j,l)
130			500 CONTINUE
131
132
133			C ---------------------------------------------------------
134			C ---------------------------------------------------------
135			C ---------------------------------------------------------
136
137			C start here
138			C
139			C boucle principale sur les niveaux et les latitudes
140			C
141			DO 1 L=1,NIV
142			DO 1 K=lati,latf
143			C
144			C initialisation
145			C
146			C program assumes periodic boundaries in X
147			C
148			DO 10 I=2,LON
149			SMNEW(I)=SM(I,K,L)+(UGRI(I-1,K,L)-UGRI(I,K,L))*DTX
150			10 CONTINUE
151			SMNEW(1)=SM(1,K,L)+(UGRI(LON,K,L)-UGRI(1,K,L))*DTX
152			C
153			C modifications for extended polar zones
154			C
155			NUMK=NUM(K)
156			LONK=LON/NUMK
157			C
158			IF(NUMK.GT.1) THEN
159			C
160			DO 111 I=1,LON
161			TM(I)=0.
162			111 CONTINUE
163			DO 112 JV=1,NTRA
164			DO 1120 I=1,LON
165			T0(I,JV)=0.
166			TX(I,JV)=0.
167			TY(I,JV)=0.
168			TZ(I,JV)=0.
169			1120 CONTINUE
170			112 CONTINUE
171			C
172			DO 11 I2=1,NUMK
173			C
174			DO 113 I=1,LONK
175			I3=(I-1)*NUMK+I2
176			TM(I)=TM(I)+SM(I3,K,L)
177			ALF(I)=SM(I3,K,L)/TM(I)
178			ALF1(I)=1.-ALF(I)
179			113 CONTINUE
180			C
181			DO JV=1,NTRA
182			DO I=1,LONK
183			I3=(I-1)*NUMK+I2
184			TEMPTM=-ALF(I)*T0(I,JV)+ALF1(I)
185			$ *S0(I3,K,L,JV)
186			T0(I,JV)=T0(I,JV)+S0(I3,K,L,JV)
187			TX(I,JV)=ALF(I) *sx(I3,K,L,JV)+
188			$ ALF1(I)TX(I,JV) +3.TEMPTM
189			TY(I,JV)=TY(I,JV)+sy(I3,K,L,JV)
190			TZ(I,JV)=TZ(I,JV)+sz(I3,K,L,JV)
191			ENDDO
192			ENDDO
193			C
194			11 CONTINUE
195			C
196			ELSE
197			C
198			DO 115 I=1,LON
199			TM(I)=SM(I,K,L)
200			115 CONTINUE
201			DO 116 JV=1,NTRA
202			DO 1160 I=1,LON
203			T0(I,JV)=S0(I,K,L,JV)
204			TX(I,JV)=sx(I,K,L,JV)
205			TY(I,JV)=sy(I,K,L,JV)
206			TZ(I,JV)=sz(I,K,L,JV)
207			1160 CONTINUE
208			116 CONTINUE
209			C
210			ENDIF
211			C
212			DO 117 I=1,LONK
213			UEXT(I)=UGRI(I*NUMK,K,L)
214			117 CONTINUE
215			C
216			C place limits on appropriate moments before transport
217			C (if flux-limiting is to be applied)
218			C
219			IF(.NOT.LIMIT) GO TO 13
220			C
221			DO 12 JV=1,NTRA
222			DO 120 I=1,LONK
223			TX(I,JV)=SIGN(AMIN1(AMAX1(T0(I,JV),0.),ABS(TX(I,JV))),TX(I,JV))
224			120 CONTINUE
225			12 CONTINUE
226			C
227			13 CONTINUE
228			C
229			C calculate flux and moments between adjacent boxes
230			C 1- create temporary moments/masses for partial boxes in transit
231			C 2- reajusts moments remaining in the box
232			C
233			C flux from IP to I if U(I).lt.0
234			C
235			DO 140 I=1,LONK-1
236			IF(UEXT(I).LT.0.) THEN
237			FM(I)=-UEXT(I)*DTX
238			ALF(I)=FM(I)/TM(I+1)
239			TM(I+1)=TM(I+1)-FM(I)
240			ENDIF
241			140 CONTINUE
242			C
243			I=LONK
244			IF(UEXT(I).LT.0.) THEN
245			FM(I)=-UEXT(I)*DTX
246			ALF(I)=FM(I)/TM(1)
247			TM(1)=TM(1)-FM(I)
248			ENDIF
249			C
250			C flux from I to IP if U(I).gt.0
251			C
252			DO 141 I=1,LONK
253			IF(UEXT(I).GE.0.) THEN
254			FM(I)=UEXT(I)*DTX
255			ALF(I)=FM(I)/TM(I)
256			TM(I)=TM(I)-FM(I)
257			ENDIF
258			141 CONTINUE
259			C
260			DO 142 I=1,LONK
261			ALFQ(I)=ALF(I)*ALF(I)
262			ALF1(I)=1.-ALF(I)
263			ALF1Q(I)=ALF1(I)*ALF1(I)
264			142 CONTINUE
265			C
266			DO 150 JV=1,NTRA
267			DO 1500 I=1,LONK-1
268			C
269			IF(UEXT(I).LT.0.) THEN
270			C
271			F0(I,JV)=ALF (I)* ( T0(I+1,JV)-ALF1(I)*TX(I+1,JV) )
272			FX(I,JV)=ALFQ(I)*TX(I+1,JV)
273			FY(I,JV)=ALF (I)*TY(I+1,JV)
274			FZ(I,JV)=ALF (I)*TZ(I+1,JV)
275			C
276			T0(I+1,JV)=T0(I+1,JV)-F0(I,JV)
277			TX(I+1,JV)=ALF1Q(I)*TX(I+1,JV)
278			TY(I+1,JV)=TY(I+1,JV)-FY(I,JV)
279			TZ(I+1,JV)=TZ(I+1,JV)-FZ(I,JV)
280			C
281			ENDIF
282			C
283			1500 CONTINUE
284			150 CONTINUE
285			C
286			I=LONK
287			IF(UEXT(I).LT.0.) THEN
288			C
289			DO 151 JV=1,NTRA
290			C
291			F0 (I,JV)=ALF (I)* ( T0(1,JV)-ALF1(I)*TX(1,JV) )
292			FX (I,JV)=ALFQ(I)*TX(1,JV)
293			FY (I,JV)=ALF (I)*TY(1,JV)
294			FZ (I,JV)=ALF (I)*TZ(1,JV)
295			C
296			T0(1,JV)=T0(1,JV)-F0(I,JV)
297			TX(1,JV)=ALF1Q(I)*TX(1,JV)
298			TY(1,JV)=TY(1,JV)-FY(I,JV)
299			TZ(1,JV)=TZ(1,JV)-FZ(I,JV)
300			C
301			151 CONTINUE
302			C
303			ENDIF
304			C
305			DO 152 JV=1,NTRA
306			DO 1520 I=1,LONK
307			C
308			IF(UEXT(I).GE.0.) THEN
309			C
310			F0(I,JV)=ALF (I)* ( T0(I,JV)+ALF1(I)*TX(I,JV) )
311			FX(I,JV)=ALFQ(I)*TX(I,JV)
312			FY(I,JV)=ALF (I)*TY(I,JV)
313			FZ(I,JV)=ALF (I)*TZ(I,JV)
314			C
315			T0(I,JV)=T0(I,JV)-F0(I,JV)
316			TX(I,JV)=ALF1Q(I)*TX(I,JV)
317			TY(I,JV)=TY(I,JV)-FY(I,JV)
318			TZ(I,JV)=TZ(I,JV)-FZ(I,JV)
319			C
320			ENDIF
321			C
322			1520 CONTINUE
323			152 CONTINUE
324			C
325			C puts the temporary moments Fi into appropriate neighboring boxes
326			C
327			DO 160 I=1,LONK
328			IF(UEXT(I).LT.0.) THEN
329			TM(I)=TM(I)+FM(I)
330			ALF(I)=FM(I)/TM(I)
331			ENDIF
332			160 CONTINUE
333			C
334			DO 161 I=1,LONK-1
335			IF(UEXT(I).GE.0.) THEN
336			TM(I+1)=TM(I+1)+FM(I)
337			ALF(I)=FM(I)/TM(I+1)
338			ENDIF
339			161 CONTINUE
340			C
341			I=LONK
342			IF(UEXT(I).GE.0.) THEN
343			TM(1)=TM(1)+FM(I)
344			ALF(I)=FM(I)/TM(1)
345			ENDIF
346			C
347			DO 162 I=1,LONK
348			ALF1(I)=1.-ALF(I)
349			162 CONTINUE
350			C
351			DO 170 JV=1,NTRA
352			DO 1700 I=1,LONK
353			C
354			IF(UEXT(I).LT.0.) THEN
355			C
356			TEMPTM=-ALF(I)T0(I,JV)+ALF1(I)F0(I,JV)
357			T0(I,JV)=T0(I,JV)+F0(I,JV)
358			TX(I,JV)=ALF(I)FX(I,JV)+ALF1(I)TX(I,JV)+3.*TEMPTM
359			TY(I,JV)=TY(I,JV)+FY(I,JV)
360			TZ(I,JV)=TZ(I,JV)+FZ(I,JV)
361			C
362			ENDIF
363			C
364			1700 CONTINUE
365			170 CONTINUE
366			C
367			DO 171 JV=1,NTRA
368			DO 1710 I=1,LONK-1
369			C
370			IF(UEXT(I).GE.0.) THEN
371			C
372			TEMPTM=ALF(I)T0(I+1,JV)-ALF1(I)F0(I,JV)
373			T0(I+1,JV)=T0(I+1,JV)+F0(I,JV)
374			TX(I+1,JV)=ALF(I)FX(I,JV)+ALF1(I)TX(I+1,JV)+3.*TEMPTM
375			TY(I+1,JV)=TY(I+1,JV)+FY(I,JV)
376			TZ(I+1,JV)=TZ(I+1,JV)+FZ(I,JV)
377			C
378			ENDIF
379			C
380			1710 CONTINUE
381			171 CONTINUE
382			C
383			I=LONK
384			IF(UEXT(I).GE.0.) THEN
385			DO 172 JV=1,NTRA
386			TEMPTM=ALF(I)T0(1,JV)-ALF1(I)F0(I,JV)
387			T0(1,JV)=T0(1,JV)+F0(I,JV)
388			TX(1,JV)=ALF(I)FX(I,JV)+ALF1(I)TX(1,JV)+3.*TEMPTM
389			TY(1,JV)=TY(1,JV)+FY(I,JV)
390			TZ(1,JV)=TZ(1,JV)+FZ(I,JV)
391			172 CONTINUE
392			ENDIF
393			C
394			C retour aux mailles d'origine (passage des Tij aux Sij)
395			C
396			IF(NUMK.GT.1) THEN
397			C
398			DO 180 I2=1,NUMK
399			C
400			DO 180 I=1,LONK
401			C
402			I3=I2+(I-1)*NUMK
403			SM(I3,K,L)=SMNEW(I3)
404			ALF(I)=SMNEW(I3)/TM(I)
405			TM(I)=TM(I)-SMNEW(I3)
406			C
407			ALFQ(I)=ALF(I)*ALF(I)
408			ALF1(I)=1.-ALF(I)
409			ALF1Q(I)=ALF1(I)*ALF1(I)
410			C
411			180 CONTINUE
412			C
413			DO JV=1,NTRA
414			DO I=1,LONK
415			C
416			I3=I2+(I-1)*NUMK
417			S0(I3,K,L,JV)=ALF (I)
418			$ * (T0(I,JV)-ALF1(I)*TX(I,JV))
419			sx(I3,K,L,JV)=ALFQ(I)*TX(I,JV)
420			sy(I3,K,L,JV)=ALF (I)*TY(I,JV)
421			sz(I3,K,L,JV)=ALF (I)*TZ(I,JV)
422			C
423			C reajusts moments remaining in the box
424			C
425			T0(I,JV)=T0(I,JV)-S0(I3,K,L,JV)
426			TX(I,JV)=ALF1Q(I)*TX(I,JV)
427			TY(I,JV)=TY(I,JV)-sy(I3,K,L,JV)
428			TZ(I,JV)=TZ(I,JV)-sz(I3,K,L,JV)
429			ENDDO
430			ENDDO
431			C
432			C
433			ELSE
434			C
435			DO 190 I=1,LON
436			SM(I,K,L)=TM(I)
437			190 CONTINUE
438			DO 191 JV=1,NTRA
439			DO 1910 I=1,LON
440			S0(I,K,L,JV)=T0(I,JV)
441			sx(I,K,L,JV)=TX(I,JV)
442			sy(I,K,L,JV)=TY(I,JV)
443			sz(I,K,L,JV)=TZ(I,JV)
444			1910 CONTINUE
445			191 CONTINUE
446			C
447			ENDIF
448			C
449			1 CONTINUE
450			C
451			C ---------- bouclage cyclique
452			DO itrac=1,ntra
453			DO l = 1,llm
454			DO j = lati,latf
455			SM(iip1,j,l) = SM(1,j,l)
456			S0(iip1,j,l,itrac) = S0(1,j,l,itrac)
457			sx(iip1,j,l,itrac) = sx(1,j,l,itrac)
458			sy(iip1,j,l,itrac) = sy(1,j,l,itrac)
459			sz(iip1,j,l,itrac) = sz(1,j,l,itrac)
460			END DO
461			END DO
462			ENDDO
463
464			c ----------- qqtite totale de traceur dans tte l'atmosphere
465			DO l = 1, llm
466			DO j = 1, jjp1
467			DO i = 1, iim
468			cIM 240405 sqf = sqf + S0(i,j,l,9)
469			sqf = sqf + S0(i,j,l,ntra)
470			END DO
471			END DO
472			END DO
473			c
474			PRINT*,'------ DIAG DANS ADVX - SORTIE -----'
475			PRINT*,'sqf=',sqf
476			c-------------
477
478			RETURN
479			END
480			C_________________________________________________________________
481			C_________________________________________________________________