1 |
|
2 |
! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/advx.F,v 1.2 2005/05/25 13:10:09 |
3 |
! fairhead Exp $ |
4 |
|
5 |
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 |
|
12 |
! 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 |
|
27 |
! parametres principaux du modele |
28 |
|
29 |
|
30 |
! 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 |
|
35 |
INTEGER ntra |
36 |
PARAMETER (ntra=1) |
37 |
|
38 |
! ATTENTION partout ou on trouve ntra, insertion de boucle |
39 |
! possible dans l'avenir. |
40 |
|
41 |
REAL dtx |
42 |
REAL, INTENT (IN) :: pbaru(iip1, jjp1, llm) |
43 |
|
44 |
! 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 |
|
48 |
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 |
|
52 |
! Local : |
53 |
! ------- |
54 |
|
55 |
! mass fluxes across the boundaries (UGRI,VGRI,WGRI) |
56 |
! mass fluxes in kg |
57 |
! declaration : |
58 |
|
59 |
REAL ugri(iip1, jjp1, llm) |
60 |
|
61 |
! Rem : VGRI et WGRI ne sont pas utilises dans |
62 |
! cette subroutine ( advection en x uniquement ) |
63 |
|
64 |
! Ti are the moments for the current latitude and level |
65 |
|
66 |
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 |
|
71 |
! the moments F are similarly defined and used as temporary |
72 |
! storage for portions of the grid boxes in transit |
73 |
|
74 |
REAL fm(iim) |
75 |
REAL f0(iim, ntra), fx(iim, ntra) |
76 |
REAL fy(iim, ntra), fz(iim, ntra) |
77 |
|
78 |
! work arrays |
79 |
|
80 |
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 |
ELSE |
197 |
|
198 |
DO i = 1, lon |
199 |
tm(i) = sm(i, k, l) |
200 |
END DO |
201 |
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 |
END DO |
215 |
|
216 |
! 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 |
END DO |
225 |
END DO |
226 |
|
227 |
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 |
! _________________________________________________________________ |