trunk/dyn3d/advz.f90


! $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
2	! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/advz.F,v 1.2 2005/05/25 13:10:09
3	! fairhead Exp $
4
5	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
12	! 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
26	! parametres principaux du modele
27
28
29	! Arguments :
30	! -----------
31	! dtz : frequence fictive d'appel du transport
32	! w : flux de masse en z en Pa.m2.s-1
33
34	INTEGER ntra
35	PARAMETER (ntra=1)
36
37	REAL, INTENT (IN) :: dtz
38	REAL w(iip1, jjp1, llm)
39
40	! 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
44	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
48
49	! Local :
50	! -------
51
52	! mass fluxes across the boundaries (UGRI,VGRI,WGRI)
53	! mass fluxes in kg
54	! declaration :
55
56	REAL wgri(iip1, jjp1, 0:llm)
57
58
59	! the moments F are used as temporary storage for
60	! portions of grid boxes in transit at the current latitude
61
62	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	! _______________________________________________________________