trunk/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 disvert_m
  IMPLICIT NONE

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

  ! parametres principaux du modele


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

  INTEGER ntra
  PARAMETER (ntra=1)

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

  REAL dtx
  REAL, INTENT (IN) :: pbaru(iip1, jjp1, llm)

  ! moments: SM  total mass in each grid box
  ! S0  mass of tracer in each grid box
  ! Si  1rst order moment in i direction

  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)

  ! Local :
  ! -------

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

  REAL ugri(iip1, jjp1, llm)

  ! Rem : VGRI et WGRI ne sont pas utilises dans
  ! cette subroutine ( advection en x uniquement )

  ! Ti are the moments for the current latitude and level

  REAL tm(iim)
  REAL t0(iim, ntra), tx(iim, ntra)
  REAL ty(iim, ntra), tz(iim, ntra)
  REAL temptm ! just a temporary variable

  ! the moments F are similarly defined and used as temporary
  ! storage for portions of the grid boxes in transit

  REAL fm(iim)
  REAL f0(iim, ntra), fx(iim, ntra)
  REAL fy(iim, ntra), fz(iim, ntra)

  ! work arrays

  REAL alf(iim), alf1(iim), alfq(iim), alf1q(iim)

  REAL smnew(iim), uext(iim)

  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

  ! *** Test de passage d'arguments ******


  ! -------------------------------------
  DO j = 1, jjp1
    num(j) = 1
  END DO
  sqi = 0.
  sqf = 0.

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


  ! Interface : adaptation nouveau modele
  ! -------------------------------------

  ! ---------------------------------------------------------
  ! Conversion des flux de masses en kg/s
  ! pbaru est en N/s d'ou :
  ! ugri est en kg/s

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


  ! ---------------------------------------------------------
  ! ---------------------------------------------------------
  ! ---------------------------------------------------------

  ! start here

  ! boucle principale sur les niveaux et les latitudes

  DO l = 1, niv
    DO k = lati, latf

      ! initialisation

      ! program assumes periodic boundaries in X

      DO i = 2, lon
        smnew(i) = sm(i, k, l) + (ugri(i-1,k,l)-ugri(i,k,l))*dtx
      END DO
      smnew(1) = sm(1, k, l) + (ugri(lon,k,l)-ugri(1,k,l))*dtx

      ! modifications for extended polar zones

      numk = num(k)
      lonk = lon/numk

      IF (numk>1) THEN

        DO i = 1, lon
          tm(i) = 0.
        END DO
        DO jv = 1, ntra
          DO i = 1, lon
            t0(i, jv) = 0.
            tx(i, jv) = 0.
            ty(i, jv) = 0.
            tz(i, jv) = 0.
          END DO
        END DO

        DO i2 = 1, numk

          DO 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)
          END DO

          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)
            END DO
          END DO

        END DO

      ELSE

        DO i = 1, lon
          tm(i) = sm(i, k, l)
        END DO
        DO jv = 1, ntra
          DO 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)
          END DO
        END DO

      END IF

      DO i = 1, lonk
        uext(i) = ugri(i*numk, k, l)
      END DO

      ! place limits on appropriate moments before transport
      ! (if flux-limiting is to be applied)

      IF (.NOT. limit) GO TO 13

      DO jv = 1, ntra
        DO i = 1, lonk
          tx(i, jv) = sign(amin1(amax1(t0(i,jv),0.),abs(tx(i,jv))), tx(i,jv))
        END DO
      END DO

13    CONTINUE

      ! calculate flux and moments between adjacent boxes
      ! 1- create temporary moments/masses for partial boxes in transit
      ! 2- reajusts moments remaining in the box

      ! flux from IP to I if U(I).lt.0

      DO i = 1, lonk - 1
        IF (uext(i)<0.) THEN
          fm(i) = -uext(i)*dtx
          alf(i) = fm(i)/tm(i+1)
          tm(i+1) = tm(i+1) - fm(i)
        END IF
      END DO

      i = lonk
      IF (uext(i)<0.) THEN
        fm(i) = -uext(i)*dtx
        alf(i) = fm(i)/tm(1)
        tm(1) = tm(1) - fm(i)
      END IF

      ! flux from I to IP if U(I).gt.0

      DO i = 1, lonk
        IF (uext(i)>=0.) THEN
          fm(i) = uext(i)*dtx
          alf(i) = fm(i)/tm(i)
          tm(i) = tm(i) - fm(i)
        END IF
      END DO

      DO i = 1, lonk
        alfq(i) = alf(i)*alf(i)
        alf1(i) = 1. - alf(i)
        alf1q(i) = alf1(i)*alf1(i)
      END DO

      DO jv = 1, ntra
        DO i = 1, lonk - 1

          IF (uext(i)<0.) THEN

            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)

            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)

          END IF

        END DO
      END DO

      i = lonk
      IF (uext(i)<0.) THEN

        DO jv = 1, ntra

          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)

          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)

        END DO

      END IF

      DO jv = 1, ntra
        DO i = 1, lonk

          IF (uext(i)>=0.) THEN

            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)

            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)

          END IF

        END DO
      END DO

      ! puts the temporary moments Fi into appropriate neighboring boxes

      DO i = 1, lonk
        IF (uext(i)<0.) THEN
          tm(i) = tm(i) + fm(i)
          alf(i) = fm(i)/tm(i)
        END IF
      END DO

      DO i = 1, lonk - 1
        IF (uext(i)>=0.) THEN
          tm(i+1) = tm(i+1) + fm(i)
          alf(i) = fm(i)/tm(i+1)
        END IF
      END DO

      i = lonk
      IF (uext(i)>=0.) THEN
        tm(1) = tm(1) + fm(i)
        alf(i) = fm(i)/tm(1)
      END IF

      DO i = 1, lonk
        alf1(i) = 1. - alf(i)
      END DO

      DO jv = 1, ntra
        DO i = 1, lonk

          IF (uext(i)<0.) THEN

            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)

          END IF

        END DO
      END DO

      DO jv = 1, ntra
        DO i = 1, lonk - 1

          IF (uext(i)>=0.) THEN

            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)

          END IF

        END DO
      END DO

      i = lonk
      IF (uext(i)>=0.) THEN
        DO 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)
        END DO
      END IF

      ! retour aux mailles d'origine (passage des Tij aux Sij)

      IF (numk>1) THEN

        DO i2 = 1, numk

          DO i = 1, lonk

            i3 = i2 + (i-1)*numk
            sm(i3, k, l) = smnew(i3)
            alf(i) = smnew(i3)/tm(i)
            tm(i) = tm(i) - smnew(i3)

            alfq(i) = alf(i)*alf(i)
            alf1(i) = 1. - alf(i)
            alf1q(i) = alf1(i)*alf1(i)

          END DO
        END DO

        DO jv = 1, ntra
          DO i = 1, lonk

            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)

            ! reajusts moments remaining in the box

            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)
          END DO
        END DO


      ELSE

        DO i = 1, lon
          sm(i, k, l) = tm(i)
        END DO
        DO jv = 1, ntra
          DO 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)
          END DO
        END DO

      END IF

    END DO
  END DO

  ! ---------- 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
  END DO

  ! ----------- qqtite totale de traceur dans tte l'atmosphere
  DO l = 1, llm
    DO j = 1, jjp1
      DO i = 1, iim
        ! IM 240405          sqf = sqf + S0(i,j,l,9)
        sqf = sqf + s0(i, j, l, ntra)
      END DO
    END DO
  END DO

  PRINT *, '------ DIAG DANS ADVX - SORTIE -----'
  PRINT *, 'sqf=', sqf
  ! -------------

  RETURN
END SUBROUTINE advx
! _________________________________________________________________
! _________________________________________________________________
1	guez	3
2	guez	81	! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/advx.F,v 1.2 2005/05/25 13:10:09
3			! fairhead Exp $
4	guez	3
5	guez	81	SUBROUTINE advx(limit, dtx, pbaru, sm, s0, sx, sy, sz, lati, latf)
6			USE dimens_m
7			USE paramet_m
8			USE comconst
9			USE disvert_m
10			IMPLICIT NONE
11	guez	3
12	guez	81	! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
13			! C
14			! first-order moments (FOM) advection of tracer in X direction C
15			! C
16			! Source : Pascal Simon (Meteo,CNRM) C
17			! Adaptation : A.Armengaud (LGGE) juin 94 C
18			! C
19			! limit,dtx,pbaru,pbarv,sm,s0,sx,sy,sz C
20			! sont des arguments d'entree pour le s-pg... C
21			! C
22			! sm,s0,sx,sy,sz C
23			! sont les arguments de sortie pour le s-pg C
24			! C
25			! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
26	guez	3
27	guez	81	! parametres principaux du modele
28	guez	3
29
30	guez	81	! Arguments :
31			! -----------
32			! dtx : frequence fictive d'appel du transport
33			! pbaru, pbarv : flux de masse en x et y en Pa.m2.s-1
34	guez	3
35	guez	81	INTEGER ntra
36			PARAMETER (ntra=1)
37	guez	3
38	guez	81	! ATTENTION partout ou on trouve ntra, insertion de boucle
39			! possible dans l'avenir.
40	guez	3
41	guez	81	REAL dtx
42			REAL, INTENT (IN) :: pbaru(iip1, jjp1, llm)
43	guez	3
44	guez	81	! moments: SM total mass in each grid box
45			! S0 mass of tracer in each grid box
46			! Si 1rst order moment in i direction
47	guez	3
48	guez	81	REAL sm(iip1, jjp1, llm), s0(iip1, jjp1, llm, ntra)
49			REAL sx(iip1, jjp1, llm, ntra), sy(iip1, jjp1, llm, ntra)
50			REAL sz(iip1, jjp1, llm, ntra)
51	guez	3
52	guez	81	! Local :
53			! -------
54	guez	3
55	guez	81	! mass fluxes across the boundaries (UGRI,VGRI,WGRI)
56			! mass fluxes in kg
57			! declaration :
58	guez	3
59	guez	81	REAL ugri(iip1, jjp1, llm)
60	guez	3
61	guez	81	! Rem : VGRI et WGRI ne sont pas utilises dans
62			! cette subroutine ( advection en x uniquement )
63	guez	3
64	guez	81	! Ti are the moments for the current latitude and level
65	guez	3
66	guez	81	REAL tm(iim)
67			REAL t0(iim, ntra), tx(iim, ntra)
68			REAL ty(iim, ntra), tz(iim, ntra)
69			REAL temptm ! just a temporary variable
70	guez	3
71	guez	81	! the moments F are similarly defined and used as temporary
72			! storage for portions of the grid boxes in transit
73	guez	3
74	guez	81	REAL fm(iim)
75			REAL f0(iim, ntra), fx(iim, ntra)
76			REAL fy(iim, ntra), fz(iim, ntra)
77	guez	3
78	guez	81	! work arrays
79	guez	3
80	guez	81	REAL alf(iim), alf1(iim), alfq(iim), alf1q(iim)
81
82			REAL smnew(iim), uext(iim)
83
84			REAL sqi, sqf
85
86			LOGICAL limit
87			INTEGER num(jjp1), lonk, numk
88			INTEGER lon, lati, latf, niv
89			INTEGER i, i2, i3, j, jv, l, k, itrac
90
91			lon = iim
92			niv = llm
93
94			! * Test de passage d'arguments ****
95
96
97			! -------------------------------------
98			DO j = 1, jjp1
99			num(j) = 1
100			END DO
101			sqi = 0.
102			sqf = 0.
103
104			DO l = 1, llm
105			DO j = 1, jjp1
106			DO i = 1, iim
107			! IM 240305 sqi = sqi + S0(i,j,l,9)
108			sqi = sqi + s0(i, j, l, ntra)
109			END DO
110			END DO
111			END DO
112			PRINT *, '-------- DIAG DANS ADVX - ENTREE ---------'
113			PRINT *, 'sqi=', sqi
114
115
116			! Interface : adaptation nouveau modele
117			! -------------------------------------
118
119			! ---------------------------------------------------------
120			! Conversion des flux de masses en kg/s
121			! pbaru est en N/s d'ou :
122			! ugri est en kg/s
123
124			DO l = 1, llm
125			DO j = 1, jjm + 1
126			DO i = 1, iip1
127			! ugri (i,j,llm+1-l) = pbaru (i,j,l) * ( dsig(l) / g )
128			ugri(i, j, llm+1-l) = pbaru(i, j, l)
129			END DO
130			END DO
131			END DO
132
133
134			! ---------------------------------------------------------
135			! ---------------------------------------------------------
136			! ---------------------------------------------------------
137
138			! start here
139
140			! boucle principale sur les niveaux et les latitudes
141
142			DO l = 1, niv
143			DO k = lati, latf
144
145			! initialisation
146
147			! program assumes periodic boundaries in X
148
149			DO i = 2, lon
150			smnew(i) = sm(i, k, l) + (ugri(i-1,k,l)-ugri(i,k,l))*dtx
151			END DO
152			smnew(1) = sm(1, k, l) + (ugri(lon,k,l)-ugri(1,k,l))*dtx
153
154			! modifications for extended polar zones
155
156			numk = num(k)
157			lonk = lon/numk
158
159			IF (numk>1) THEN
160
161			DO i = 1, lon
162			tm(i) = 0.
163			END DO
164			DO jv = 1, ntra
165			DO i = 1, lon
166			t0(i, jv) = 0.
167			tx(i, jv) = 0.
168			ty(i, jv) = 0.
169			tz(i, jv) = 0.
170			END DO
171			END DO
172
173			DO i2 = 1, numk
174
175			DO i = 1, lonk
176			i3 = (i-1)*numk + i2
177			tm(i) = tm(i) + sm(i3, k, l)
178			alf(i) = sm(i3, k, l)/tm(i)
179			alf1(i) = 1. - alf(i)
180			END DO
181
182			DO jv = 1, ntra
183			DO i = 1, lonk
184			i3 = (i-1)*numk + i2
185			temptm = -alf(i)t0(i, jv) + alf1(i)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) + alf1(i)tx(i, jv) + &
188			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			END DO
192			END DO
193
194			END DO
195
196	guez	3	ELSE
197	guez	81
198			DO i = 1, lon
199			tm(i) = sm(i, k, l)
200	guez	3	END DO
201	guez	81	DO jv = 1, ntra
202			DO 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			END DO
208			END DO
209
210			END IF
211
212			DO i = 1, lonk
213			uext(i) = ugri(i*numk, k, l)
214	guez	3	END DO
215
216	guez	81	! place limits on appropriate moments before transport
217			! (if flux-limiting is to be applied)
218
219			IF (.NOT. limit) GO TO 13
220
221			DO jv = 1, ntra
222			DO i = 1, lonk
223			tx(i, jv) = sign(amin1(amax1(t0(i,jv),0.),abs(tx(i,jv))), tx(i,jv))
224	guez	3	END DO
225			END DO
226
227	guez	81	13 CONTINUE
228
229			! calculate flux and moments between adjacent boxes
230			! 1- create temporary moments/masses for partial boxes in transit
231			! 2- reajusts moments remaining in the box
232
233			! flux from IP to I if U(I).lt.0
234
235			DO i = 1, lonk - 1
236			IF (uext(i)<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			END IF
241			END DO
242
243			i = lonk
244			IF (uext(i)<0.) THEN
245			fm(i) = -uext(i)*dtx
246			alf(i) = fm(i)/tm(1)
247			tm(1) = tm(1) - fm(i)
248			END IF
249
250			! flux from I to IP if U(I).gt.0
251
252			DO i = 1, lonk
253			IF (uext(i)>=0.) THEN
254			fm(i) = uext(i)*dtx
255			alf(i) = fm(i)/tm(i)
256			tm(i) = tm(i) - fm(i)
257			END IF
258			END DO
259
260			DO i = 1, lonk
261			alfq(i) = alf(i)*alf(i)
262			alf1(i) = 1. - alf(i)
263			alf1q(i) = alf1(i)*alf1(i)
264			END DO
265
266			DO jv = 1, ntra
267			DO i = 1, lonk - 1
268
269			IF (uext(i)<0.) THEN
270
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
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
281			END IF
282
283			END DO
284			END DO
285
286			i = lonk
287			IF (uext(i)<0.) THEN
288
289			DO jv = 1, ntra
290
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
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
301			END DO
302
303			END IF
304
305			DO jv = 1, ntra
306			DO i = 1, lonk
307
308			IF (uext(i)>=0.) THEN
309
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
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
320			END IF
321
322			END DO
323			END DO
324
325			! puts the temporary moments Fi into appropriate neighboring boxes
326
327			DO i = 1, lonk
328			IF (uext(i)<0.) THEN
329			tm(i) = tm(i) + fm(i)
330			alf(i) = fm(i)/tm(i)
331			END IF
332			END DO
333
334			DO i = 1, lonk - 1
335			IF (uext(i)>=0.) THEN
336			tm(i+1) = tm(i+1) + fm(i)
337			alf(i) = fm(i)/tm(i+1)
338			END IF
339			END DO
340
341			i = lonk
342			IF (uext(i)>=0.) THEN
343			tm(1) = tm(1) + fm(i)
344			alf(i) = fm(i)/tm(1)
345			END IF
346
347			DO i = 1, lonk
348			alf1(i) = 1. - alf(i)
349			END DO
350
351			DO jv = 1, ntra
352			DO i = 1, lonk
353
354			IF (uext(i)<0.) THEN
355
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
362			END IF
363
364			END DO
365			END DO
366
367			DO jv = 1, ntra
368			DO i = 1, lonk - 1
369
370			IF (uext(i)>=0.) THEN
371
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
378			END IF
379
380			END DO
381			END DO
382
383			i = lonk
384			IF (uext(i)>=0.) THEN
385			DO 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			END DO
392			END IF
393
394			! retour aux mailles d'origine (passage des Tij aux Sij)
395
396			IF (numk>1) THEN
397
398			DO i2 = 1, numk
399
400			DO i = 1, lonk
401
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
407			alfq(i) = alf(i)*alf(i)
408			alf1(i) = 1. - alf(i)
409			alf1q(i) = alf1(i)*alf1(i)
410
411			END DO
412			END DO
413
414			DO jv = 1, ntra
415			DO i = 1, lonk
416
417			i3 = i2 + (i-1)*numk
418			s0(i3, k, l, jv) = alf(i)(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
423			! reajusts moments remaining in the box
424
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			END DO
430			END DO
431
432
433			ELSE
434
435			DO i = 1, lon
436			sm(i, k, l) = tm(i)
437			END DO
438			DO jv = 1, ntra
439			DO 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			END DO
445			END DO
446
447			END IF
448
449			END DO
450			END DO
451
452			! ---------- bouclage cyclique
453			DO itrac = 1, ntra
454			DO l = 1, llm
455			DO j = lati, latf
456			sm(iip1, j, l) = sm(1, j, l)
457			s0(iip1, j, l, itrac) = s0(1, j, l, itrac)
458			sx(iip1, j, l, itrac) = sx(1, j, l, itrac)
459			sy(iip1, j, l, itrac) = sy(1, j, l, itrac)
460			sz(iip1, j, l, itrac) = sz(1, j, l, itrac)
461			END DO
462			END DO
463			END DO
464
465			! ----------- qqtite totale de traceur dans tte l'atmosphere
466			DO l = 1, llm
467			DO j = 1, jjp1
468			DO i = 1, iim
469			! IM 240405 sqf = sqf + S0(i,j,l,9)
470			sqf = sqf + s0(i, j, l, ntra)
471			END DO
472			END DO
473			END DO
474
475			PRINT *, '------ DIAG DANS ADVX - SORTIE -----'
476			PRINT *, 'sqf=', sqf
477			! -------------
478
479			RETURN
480			END SUBROUTINE advx
481			! _________________________________________________________________
482			! _________________________________________________________________