trunk/dyn3d/advz.f


! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/advz.F,v 1.2 2005/05/25 13:10:09
! fairhead Exp $

SUBROUTINE advz(limit, dtz, w, sm, s0, sx, sy, sz)
  USE dimens_m
  USE paramet_m
  USE comconst
  USE disvert_m
  IMPLICIT NONE

  ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  ! C
  ! first-order moments (FOM) advection of tracer in Z direction  C
  ! C
  ! Source : Pascal Simon (Meteo,CNRM)                            C
  ! Adaptation : A.Armengaud (LGGE) juin 94                       C
  ! C
  ! C
  ! sont des arguments d'entree pour le s-pg...                   C
  ! C
  ! dq est l'argument de sortie pour le s-pg                      C
  ! C
  ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

  ! parametres principaux du modele


  ! Arguments :
  ! -----------
  ! dtz : frequence fictive d'appel du transport
  ! w : flux de masse en z en Pa.m2.s-1

  INTEGER ntra
  PARAMETER (ntra=1)

  REAL, INTENT (IN) :: dtz
  REAL w(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), &
    sz(iip1, jjp1, llm, ntra)


  ! Local :
  ! -------

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

  REAL wgri(iip1, jjp1, 0:llm)


  ! the moments F are used as temporary  storage for
  ! portions of grid boxes in transit at the current latitude

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

  ! work arrays

  REAL alf(iim), alf1(iim), alfq(iim), alf1q(iim)
  REAL temptm ! Just temporal variable
  REAL sqi, sqf

  LOGICAL limit
  INTEGER lon, lat, niv
  INTEGER i, j, jv, k, l, lp

  lon = iim
  lat = jjp1
  niv = llm

  ! *** Test : diag de la qqtite totale de traceur
  ! dans l'atmosphere avant l'advection en z
  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 ADVZ - ENTREE ---------'
  PRINT *, 'sqi=', sqi

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

  ! Conversion du flux de masse en kg.s-1

  DO l = 1, llm
    DO j = 1, jjp1
      DO i = 1, iip1
        ! wgri (i,j,llm+1-l) =  w (i,j,l) / g
        wgri(i, j, llm+1-l) = w(i, j, l)
        ! wgri (i,j,0) = 0.                ! a detruire ult.
        ! wgri (i,j,l) = 0.1               !    w (i,j,l)
        ! wgri (i,j,llm) = 0.              ! a detruire ult.
      END DO
    END DO
  END DO
  DO j = 1, jjp1
    DO i = 1, iip1
      wgri(i, j, 0) = 0.
    END DO
  END DO

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

  ! start here
  ! boucle sur les latitudes

  DO k = 1, lat

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

    IF (.NOT. limit) GO TO 101

    DO jv = 1, ntra
      DO l = 1, niv
        DO i = 1, lon
          sz(i, k, l, jv) = sign(amin1(amax1(s0(i,k,l,jv), &
            0.),abs(sz(i,k,l,jv))), sz(i,k,l,jv))
        END DO
      END DO
    END DO

101 CONTINUE

    ! boucle sur les niveaux intercouches de 1 a NIV-1
    ! (flux nul au sommet L=0 et a la base L=NIV)

    ! calculate flux and moments between adjacent boxes
    ! (flux from LP to L if WGRI(L).lt.0, from L to LP if WGRI(L).gt.0)
    ! 1- create temporary moments/masses for partial boxes in transit
    ! 2- reajusts moments remaining in the box

    DO l = 1, niv - 1
      lp = l + 1

      DO i = 1, lon

        IF (wgri(i,k,l)<0.) THEN
          fm(i, l) = -wgri(i, k, l)*dtz
          alf(i) = fm(i, l)/sm(i, k, lp)
          sm(i, k, lp) = sm(i, k, lp) - fm(i, l)
        ELSE
          fm(i, l) = wgri(i, k, l)*dtz
          alf(i) = fm(i, l)/sm(i, k, l)
          sm(i, k, l) = sm(i, k, l) - fm(i, l)
        END IF

        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, lon

          IF (wgri(i,k,l)<0.) THEN

            f0(i, l, jv) = alf(i)*(s0(i,k,lp,jv)-alf1(i)*sz(i,k,lp,jv))
            fz(i, l, jv) = alfq(i)*sz(i, k, lp, jv)
            fx(i, l, jv) = alf(i)*sx(i, k, lp, jv)
            fy(i, l, jv) = alf(i)*sy(i, k, lp, jv)

            s0(i, k, lp, jv) = s0(i, k, lp, jv) - f0(i, l, jv)
            sz(i, k, lp, jv) = alf1q(i)*sz(i, k, lp, jv)
            sx(i, k, lp, jv) = sx(i, k, lp, jv) - fx(i, l, jv)
            sy(i, k, lp, jv) = sy(i, k, lp, jv) - fy(i, l, jv)

          ELSE

            f0(i, l, jv) = alf(i)*(s0(i,k,l,jv)+alf1(i)*sz(i,k,l,jv))
            fz(i, l, jv) = alfq(i)*sz(i, k, l, jv)
            fx(i, l, jv) = alf(i)*sx(i, k, l, jv)
            fy(i, l, jv) = alf(i)*sy(i, k, l, jv)

            s0(i, k, l, jv) = s0(i, k, l, jv) - f0(i, l, jv)
            sz(i, k, l, jv) = alf1q(i)*sz(i, k, l, jv)
            sx(i, k, l, jv) = sx(i, k, l, jv) - fx(i, l, jv)
            sy(i, k, l, jv) = sy(i, k, l, jv) - fy(i, l, jv)

          END IF

        END DO
      END DO

    END DO

    ! puts the temporary moments Fi into appropriate neighboring boxes

    DO l = 1, niv - 1
      lp = l + 1

      DO i = 1, lon

        IF (wgri(i,k,l)<0.) THEN
          sm(i, k, l) = sm(i, k, l) + fm(i, l)
          alf(i) = fm(i, l)/sm(i, k, l)
        ELSE
          sm(i, k, lp) = sm(i, k, lp) + fm(i, l)
          alf(i) = fm(i, l)/sm(i, k, lp)
        END IF

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

      END DO

      DO jv = 1, ntra
        DO i = 1, lon

          IF (wgri(i,k,l)<0.) THEN

            temptm = -alf(i)*s0(i, k, l, jv) + alf1(i)*f0(i, l, jv)
            s0(i, k, l, jv) = s0(i, k, l, jv) + f0(i, l, jv)
            sz(i, k, l, jv) = alf(i)*fz(i, l, jv) + alf1(i)*sz(i, k, l, jv) + &
              3.*temptm
            sx(i, k, l, jv) = sx(i, k, l, jv) + fx(i, l, jv)
            sy(i, k, l, jv) = sy(i, k, l, jv) + fy(i, l, jv)

          ELSE

            temptm = alf(i)*s0(i, k, lp, jv) - alf1(i)*f0(i, l, jv)
            s0(i, k, lp, jv) = s0(i, k, lp, jv) + f0(i, l, jv)
            sz(i, k, lp, jv) = alf(i)*fz(i, l, jv) + &
              alf1(i)*sz(i, k, lp, jv) + 3.*temptm
            sx(i, k, lp, jv) = sx(i, k, lp, jv) + fx(i, l, jv)
            sy(i, k, lp, jv) = sy(i, k, lp, jv) + fy(i, l, jv)

          END IF

        END DO
      END DO

    END DO

    ! fin de la boucle principale sur les latitudes

  END DO

  ! *** ------------------- bouclage cyclique  en X ------------

  ! 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)
  ! sx(iip1,j,l,ntra) = sx(1,j,l,ntra)
  ! sy(iip1,j,l,ntra) = sy(1,j,l,ntra)
  ! sz(iip1,j,l,ntra) = sz(1,j,l,ntra)
  ! ENDDO
  ! ENDDO

  ! -------------------------------------------------------------
  ! *** Test : diag de la qqtite totale de traceur
  ! dans l'atmosphere avant l'advection en z
  DO l = 1, llm
    DO j = 1, jjp1
      DO i = 1, iim
        ! IM 240305            sqf = sqf + S0(i,j,l,9)
        sqf = sqf + s0(i, j, l, ntra)
      END DO
    END DO
  END DO
  PRINT *, '-------- DIAG DANS ADVZ - SORTIE ---------'
  PRINT *, 'sqf=', sqf

  ! -------------------------------------------------------------
  RETURN
END SUBROUTINE advz
! _______________________________________________________________
! _______________________________________________________________
1	guez	3
2	guez	81	! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/advz.F,v 1.2 2005/05/25 13:10:09
3			! fairhead Exp $
4	guez	3
5	guez	81	SUBROUTINE advz(limit, dtz, w, sm, s0, sx, sy, sz)
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 Z direction C
15			! C
16			! Source : Pascal Simon (Meteo,CNRM) C
17			! Adaptation : A.Armengaud (LGGE) juin 94 C
18			! C
19			! C
20			! sont des arguments d'entree pour le s-pg... C
21			! C
22			! dq est l'argument de sortie pour le s-pg C
23			! C
24			! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
25	guez	3
26	guez	81	! parametres principaux du modele
27	guez	3
28
29	guez	81	! Arguments :
30			! -----------
31			! dtz : frequence fictive d'appel du transport
32			! w : flux de masse en z en Pa.m2.s-1
33	guez	3
34	guez	81	INTEGER ntra
35			PARAMETER (ntra=1)
36	guez	3
37	guez	81	REAL, INTENT (IN) :: dtz
38			REAL w(iip1, jjp1, llm)
39	guez	3
40	guez	81	! moments: SM total mass in each grid box
41			! S0 mass of tracer in each grid box
42			! Si 1rst order moment in i direction
43	guez	3
44	guez	81	REAL sm(iip1, jjp1, llm), s0(iip1, jjp1, llm, ntra)
45			REAL sx(iip1, jjp1, llm, ntra), sy(iip1, jjp1, llm, ntra), &
46			sz(iip1, jjp1, llm, ntra)
47	guez	3
48
49	guez	81	! Local :
50			! -------
51	guez	3
52	guez	81	! mass fluxes across the boundaries (UGRI,VGRI,WGRI)
53			! mass fluxes in kg
54			! declaration :
55	guez	3
56	guez	81	REAL wgri(iip1, jjp1, 0:llm)
57	guez	3
58
59	guez	81	! the moments F are used as temporary storage for
60			! portions of grid boxes in transit at the current latitude
61	guez	3
62	guez	81	REAL fm(iim, llm)
63			REAL f0(iim, llm, ntra), fx(iim, llm, ntra)
64			REAL fy(iim, llm, ntra), fz(iim, llm, ntra)
65
66			! work arrays
67
68			REAL alf(iim), alf1(iim), alfq(iim), alf1q(iim)
69			REAL temptm ! Just temporal variable
70			REAL sqi, sqf
71
72			LOGICAL limit
73			INTEGER lon, lat, niv
74			INTEGER i, j, jv, k, l, lp
75
76			lon = iim
77			lat = jjp1
78			niv = llm
79
80			! *** Test : diag de la qqtite totale de traceur
81			! dans l'atmosphere avant l'advection en z
82			sqi = 0.
83			sqf = 0.
84
85			DO l = 1, llm
86			DO j = 1, jjp1
87			DO i = 1, iim
88			! IM 240305 sqi = sqi + S0(i,j,l,9)
89			sqi = sqi + s0(i, j, l, ntra)
90			END DO
91			END DO
92			END DO
93			PRINT *, '-------- DIAG DANS ADVZ - ENTREE ---------'
94			PRINT *, 'sqi=', sqi
95
96			! -----------------------------------------------------------------
97			! Interface : adaptation nouveau modele
98			! -------------------------------------
99
100			! Conversion du flux de masse en kg.s-1
101
102			DO l = 1, llm
103			DO j = 1, jjp1
104			DO i = 1, iip1
105			! wgri (i,j,llm+1-l) = w (i,j,l) / g
106			wgri(i, j, llm+1-l) = w(i, j, l)
107			! wgri (i,j,0) = 0. ! a detruire ult.
108			! wgri (i,j,l) = 0.1 ! w (i,j,l)
109			! wgri (i,j,llm) = 0. ! a detruire ult.
110			END DO
111			END DO
112			END DO
113			DO j = 1, jjp1
114			DO i = 1, iip1
115			wgri(i, j, 0) = 0.
116			END DO
117			END DO
118
119			! -----------------------------------------------------------------
120
121			! start here
122			! boucle sur les latitudes
123
124			DO k = 1, lat
125
126			! place limits on appropriate moments before transport
127			! (if flux-limiting is to be applied)
128
129			IF (.NOT. limit) GO TO 101
130
131			DO jv = 1, ntra
132			DO l = 1, niv
133			DO i = 1, lon
134			sz(i, k, l, jv) = sign(amin1(amax1(s0(i,k,l,jv), &
135			0.),abs(sz(i,k,l,jv))), sz(i,k,l,jv))
136			END DO
137			END DO
138			END DO
139
140			101 CONTINUE
141
142			! boucle sur les niveaux intercouches de 1 a NIV-1
143			! (flux nul au sommet L=0 et a la base L=NIV)
144
145			! calculate flux and moments between adjacent boxes
146			! (flux from LP to L if WGRI(L).lt.0, from L to LP if WGRI(L).gt.0)
147			! 1- create temporary moments/masses for partial boxes in transit
148			! 2- reajusts moments remaining in the box
149
150			DO l = 1, niv - 1
151			lp = l + 1
152
153			DO i = 1, lon
154
155			IF (wgri(i,k,l)<0.) THEN
156			fm(i, l) = -wgri(i, k, l)*dtz
157			alf(i) = fm(i, l)/sm(i, k, lp)
158			sm(i, k, lp) = sm(i, k, lp) - fm(i, l)
159			ELSE
160			fm(i, l) = wgri(i, k, l)*dtz
161			alf(i) = fm(i, l)/sm(i, k, l)
162			sm(i, k, l) = sm(i, k, l) - fm(i, l)
163			END IF
164
165			alfq(i) = alf(i)*alf(i)
166			alf1(i) = 1. - alf(i)
167			alf1q(i) = alf1(i)*alf1(i)
168
169			END DO
170
171			DO jv = 1, ntra
172			DO i = 1, lon
173
174			IF (wgri(i,k,l)<0.) THEN
175
176			f0(i, l, jv) = alf(i)(s0(i,k,lp,jv)-alf1(i)sz(i,k,lp,jv))
177			fz(i, l, jv) = alfq(i)*sz(i, k, lp, jv)
178			fx(i, l, jv) = alf(i)*sx(i, k, lp, jv)
179			fy(i, l, jv) = alf(i)*sy(i, k, lp, jv)
180
181			s0(i, k, lp, jv) = s0(i, k, lp, jv) - f0(i, l, jv)
182			sz(i, k, lp, jv) = alf1q(i)*sz(i, k, lp, jv)
183			sx(i, k, lp, jv) = sx(i, k, lp, jv) - fx(i, l, jv)
184			sy(i, k, lp, jv) = sy(i, k, lp, jv) - fy(i, l, jv)
185
186			ELSE
187
188			f0(i, l, jv) = alf(i)(s0(i,k,l,jv)+alf1(i)sz(i,k,l,jv))
189			fz(i, l, jv) = alfq(i)*sz(i, k, l, jv)
190			fx(i, l, jv) = alf(i)*sx(i, k, l, jv)
191			fy(i, l, jv) = alf(i)*sy(i, k, l, jv)
192
193			s0(i, k, l, jv) = s0(i, k, l, jv) - f0(i, l, jv)
194			sz(i, k, l, jv) = alf1q(i)*sz(i, k, l, jv)
195			sx(i, k, l, jv) = sx(i, k, l, jv) - fx(i, l, jv)
196			sy(i, k, l, jv) = sy(i, k, l, jv) - fy(i, l, jv)
197
198			END IF
199
200			END DO
201			END DO
202
203			END DO
204
205			! puts the temporary moments Fi into appropriate neighboring boxes
206
207			DO l = 1, niv - 1
208			lp = l + 1
209
210			DO i = 1, lon
211
212			IF (wgri(i,k,l)<0.) THEN
213			sm(i, k, l) = sm(i, k, l) + fm(i, l)
214			alf(i) = fm(i, l)/sm(i, k, l)
215			ELSE
216			sm(i, k, lp) = sm(i, k, lp) + fm(i, l)
217			alf(i) = fm(i, l)/sm(i, k, lp)
218			END IF
219
220			alf1(i) = 1. - alf(i)
221			alfq(i) = alf(i)*alf(i)
222			alf1q(i) = alf1(i)*alf1(i)
223
224			END DO
225
226			DO jv = 1, ntra
227			DO i = 1, lon
228
229			IF (wgri(i,k,l)<0.) THEN
230
231			temptm = -alf(i)s0(i, k, l, jv) + alf1(i)f0(i, l, jv)
232			s0(i, k, l, jv) = s0(i, k, l, jv) + f0(i, l, jv)
233			sz(i, k, l, jv) = alf(i)fz(i, l, jv) + alf1(i)sz(i, k, l, jv) + &
234			3.*temptm
235			sx(i, k, l, jv) = sx(i, k, l, jv) + fx(i, l, jv)
236			sy(i, k, l, jv) = sy(i, k, l, jv) + fy(i, l, jv)
237
238			ELSE
239
240			temptm = alf(i)s0(i, k, lp, jv) - alf1(i)f0(i, l, jv)
241			s0(i, k, lp, jv) = s0(i, k, lp, jv) + f0(i, l, jv)
242			sz(i, k, lp, jv) = alf(i)*fz(i, l, jv) + &
243			alf1(i)sz(i, k, lp, jv) + 3.temptm
244			sx(i, k, lp, jv) = sx(i, k, lp, jv) + fx(i, l, jv)
245			sy(i, k, lp, jv) = sy(i, k, lp, jv) + fy(i, l, jv)
246
247			END IF
248
249			END DO
250			END DO
251
252			END DO
253
254			! fin de la boucle principale sur les latitudes
255
256			END DO
257
258			! *** ------------------- bouclage cyclique en X ------------
259
260			! DO l = 1,llm
261			! DO j = 1,jjp1
262			! SM(iip1,j,l) = SM(1,j,l)
263			! S0(iip1,j,l,ntra) = S0(1,j,l,ntra)
264			! sx(iip1,j,l,ntra) = sx(1,j,l,ntra)
265			! sy(iip1,j,l,ntra) = sy(1,j,l,ntra)
266			! sz(iip1,j,l,ntra) = sz(1,j,l,ntra)
267			! ENDDO
268			! ENDDO
269
270			! -------------------------------------------------------------
271			! *** Test : diag de la qqtite totale de traceur
272			! dans l'atmosphere avant l'advection en z
273			DO l = 1, llm
274			DO j = 1, jjp1
275			DO i = 1, iim
276			! IM 240305 sqf = sqf + S0(i,j,l,9)
277			sqf = sqf + s0(i, j, l, ntra)
278			END DO
279			END DO
280			END DO
281			PRINT *, '-------- DIAG DANS ADVZ - SORTIE ---------'
282			PRINT *, 'sqf=', sqf
283
284			! -------------------------------------------------------------
285			RETURN
286			END SUBROUTINE advz
287			! _______________________________________________________________
288			! _______________________________________________________________