/[lmdze]/trunk/dyn3d/advxp.f
ViewVC logotype

Annotation of /trunk/dyn3d/advxp.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 82 - (hide annotations)
Wed Mar 5 14:57:53 2014 UTC (10 years, 3 months ago) by guez
File size: 18284 byte(s)
Changed all ".f90" suffixes to ".f".
1 guez 3
2 guez 81 ! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/advxp.F,v 1.1.1.1 2004/05/19
3     ! 12:53:06 lmdzadmin Exp $
4 guez 3
5 guez 81 SUBROUTINE advxp(limit, dtx, pbaru, sm, s0, ssx, sy, sz, ssxx, ssxy, ssxz, &
6     syy, syz, szz, ntra)
7     USE dimens_m
8     USE paramet_m
9     USE comconst
10     USE disvert_m
11     IMPLICIT NONE
12     ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
13     ! C
14     ! second-order moments (SOM) advection of tracer in X direction C
15     ! C
16     ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
17 guez 3
18 guez 81 ! parametres principaux du modele
19 guez 3
20    
21 guez 81 INTEGER ntra
22     ! PARAMETER (ntra = 1)
23 guez 3
24 guez 81 ! definition de la grille du modele
25 guez 3
26 guez 81 REAL dtx
27     REAL, INTENT (IN) :: pbaru(iip1, jjp1, llm)
28 guez 3
29 guez 81 ! moments: SM total mass in each grid box
30     ! S0 mass of tracer in each grid box
31     ! Si 1rst order moment in i direction
32     ! Sij 2nd order moment in i and j directions
33    
34     REAL sm(iip1, jjp1, llm), s0(iip1, jjp1, llm, ntra)
35     REAL ssx(iip1, jjp1, llm, ntra), sy(iip1, jjp1, llm, ntra), &
36     sz(iip1, jjp1, llm, ntra)
37     REAL ssxx(iip1, jjp1, llm, ntra), ssxy(iip1, jjp1, llm, ntra), &
38     ssxz(iip1, jjp1, llm, ntra), syy(iip1, jjp1, llm, ntra), &
39     syz(iip1, jjp1, llm, ntra), szz(iip1, jjp1, llm, ntra)
40    
41     ! Local :
42     ! -------
43    
44     ! mass fluxes across the boundaries (UGRI,VGRI,WGRI)
45     ! mass fluxes in kg
46     ! declaration :
47    
48     REAL ugri(iip1, jjp1, llm)
49    
50     ! Rem : VGRI et WGRI ne sont pas utilises dans
51     ! cette subroutine ( advection en x uniquement )
52    
53    
54     ! Tij are the moments for the current latitude and level
55    
56     REAL tm(iim)
57     REAL t0(iim, ntra), tx(iim, ntra)
58     REAL ty(iim, ntra), tz(iim, ntra)
59     REAL txx(iim, ntra), txy(iim, ntra)
60     REAL txz(iim, ntra), tyy(iim, ntra)
61     REAL tyz(iim, ntra), tzz(iim, ntra)
62    
63     ! the moments F are similarly defined and used as temporary
64     ! storage for portions of the grid boxes in transit
65    
66     REAL fm(iim)
67     REAL f0(iim, ntra), fx(iim, ntra)
68     REAL fy(iim, ntra), fz(iim, ntra)
69     REAL fxx(iim, ntra), fxy(iim, ntra)
70     REAL fxz(iim, ntra), fyy(iim, ntra)
71     REAL fyz(iim, ntra), fzz(iim, ntra)
72    
73     ! work arrays
74    
75     REAL alf(iim), alf1(iim), alfq(iim), alf1q(iim)
76     REAL alf2(iim), alf3(iim), alf4(iim)
77    
78     REAL smnew(iim), uext(iim)
79     REAL sqi, sqf
80     REAL temptm
81     REAL slpmax
82     REAL s1max, s1new, s2new
83    
84     LOGICAL limit
85     INTEGER num(jjp1), lonk, numk
86     INTEGER lon, lati, latf, niv
87     INTEGER i, i2, i3, j, jv, l, k, iter
88    
89     lon = iim
90     lati = 2
91     latf = jjm
92     niv = llm
93    
94     ! *** Test : diagnostique de la qtite totale de traceur
95     ! dans l'atmosphere avant l'advection
96    
97     sqi = 0.
98     sqf = 0.
99    
100     DO l = 1, llm
101     DO j = 1, jjp1
102 guez 3 DO i = 1, iim
103 guez 81 sqi = sqi + s0(i, j, l, ntra)
104 guez 3 END DO
105 guez 81 END DO
106     END DO
107     PRINT *, '------ DIAG DANS ADVX2 - ENTREE -----'
108     PRINT *, 'sqi=', sqi
109     ! test
110     ! -------------------------------------
111     DO j = 1, jjp1
112     num(j) = 1
113     END DO
114     ! DO l=1,llm
115     ! NUM(2,l)=6
116     ! NUM(3,l)=6
117     ! NUM(jjm-1,l)=6
118     ! NUM(jjm,l)=6
119     ! ENDDO
120     ! DO j=2,6
121     ! NUM(j)=12
122     ! ENDDO
123     ! DO j=jjm-5,jjm-1
124     ! NUM(j)=12
125     ! ENDDO
126    
127     ! Interface : adaptation nouveau modele
128     ! -------------------------------------
129    
130     ! ---------------------------------------------------------
131     ! Conversion des flux de masses en kg/s
132     ! pbaru est en N/s d'ou :
133     ! ugri est en kg/s
134    
135     DO l = 1, llm
136     DO j = 1, jjp1
137     DO i = 1, iip1
138     ugri(i, j, llm+1-l) = pbaru(i, j, l)
139 guez 3 END DO
140 guez 81 END DO
141     END DO
142    
143     ! ---------------------------------------------------------
144     ! start here
145    
146     ! boucle principale sur les niveaux et les latitudes
147    
148     DO l = 1, niv
149     DO k = lati, latf
150    
151    
152     ! initialisation
153    
154     ! program assumes periodic boundaries in X
155    
156     DO i = 2, lon
157     smnew(i) = sm(i, k, l) + (ugri(i-1,k,l)-ugri(i,k,l))*dtx
158 guez 3 END DO
159 guez 81 smnew(1) = sm(1, k, l) + (ugri(lon,k,l)-ugri(1,k,l))*dtx
160 guez 3
161 guez 81 ! modifications for extended polar zones
162 guez 3
163 guez 81 numk = num(k)
164     lonk = lon/numk
165 guez 3
166 guez 81 IF (numk>1) THEN
167 guez 3
168 guez 81 DO i = 1, lon
169     tm(i) = 0.
170     END DO
171     DO jv = 1, ntra
172     DO i = 1, lon
173     t0(i, jv) = 0.
174     tx(i, jv) = 0.
175     ty(i, jv) = 0.
176     tz(i, jv) = 0.
177     txx(i, jv) = 0.
178     txy(i, jv) = 0.
179     txz(i, jv) = 0.
180     tyy(i, jv) = 0.
181     tyz(i, jv) = 0.
182     tzz(i, jv) = 0.
183     END DO
184     END DO
185    
186     DO i2 = 1, numk
187    
188     DO i = 1, lonk
189     i3 = (i-1)*numk + i2
190     tm(i) = tm(i) + sm(i3, k, l)
191     alf(i) = sm(i3, k, l)/tm(i)
192     alf1(i) = 1. - alf(i)
193     alfq(i) = alf(i)*alf(i)
194     alf1q(i) = alf1(i)*alf1(i)
195     alf2(i) = alf1(i) - alf(i)
196     alf3(i) = alf(i)*alf1(i)
197     END DO
198    
199     DO jv = 1, ntra
200     DO i = 1, lonk
201     i3 = (i-1)*numk + i2
202     temptm = -alf(i)*t0(i, jv) + alf1(i)*s0(i3, k, l, jv)
203     t0(i, jv) = t0(i, jv) + s0(i3, k, l, jv)
204     txx(i, jv) = alfq(i)*ssxx(i3, k, l, jv) + alf1q(i)*txx(i, jv) + &
205     5.*(alf3(i)*(ssx(i3,k,l,jv)-tx(i,jv))+alf2(i)*temptm)
206     tx(i, jv) = alf(i)*ssx(i3, k, l, jv) + alf1(i)*tx(i, jv) + &
207     3.*temptm
208     txy(i, jv) = alf(i)*ssxy(i3, k, l, jv) + alf1(i)*txy(i, jv) + &
209     3.*(alf1(i)*sy(i3,k,l,jv)-alf(i)*ty(i,jv))
210     txz(i, jv) = alf(i)*ssxz(i3, k, l, jv) + alf1(i)*txz(i, jv) + &
211     3.*(alf1(i)*sz(i3,k,l,jv)-alf(i)*tz(i,jv))
212     ty(i, jv) = ty(i, jv) + sy(i3, k, l, jv)
213     tz(i, jv) = tz(i, jv) + sz(i3, k, l, jv)
214     tyy(i, jv) = tyy(i, jv) + syy(i3, k, l, jv)
215     tyz(i, jv) = tyz(i, jv) + syz(i3, k, l, jv)
216     tzz(i, jv) = tzz(i, jv) + szz(i3, k, l, jv)
217     END DO
218     END DO
219    
220     END DO
221    
222 guez 3 ELSE
223    
224 guez 81 DO i = 1, lon
225     tm(i) = sm(i, k, l)
226     END DO
227     DO jv = 1, ntra
228     DO i = 1, lon
229     t0(i, jv) = s0(i, k, l, jv)
230     tx(i, jv) = ssx(i, k, l, jv)
231     ty(i, jv) = sy(i, k, l, jv)
232     tz(i, jv) = sz(i, k, l, jv)
233     txx(i, jv) = ssxx(i, k, l, jv)
234     txy(i, jv) = ssxy(i, k, l, jv)
235     txz(i, jv) = ssxz(i, k, l, jv)
236     tyy(i, jv) = syy(i, k, l, jv)
237     tyz(i, jv) = syz(i, k, l, jv)
238     tzz(i, jv) = szz(i, k, l, jv)
239     END DO
240     END DO
241 guez 3
242 guez 81 END IF
243    
244     DO i = 1, lonk
245     uext(i) = ugri(i*numk, k, l)
246 guez 3 END DO
247 guez 81
248     ! place limits on appropriate moments before transport
249     ! (if flux-limiting is to be applied)
250    
251     IF (.NOT. limit) GO TO 13
252    
253     DO jv = 1, ntra
254     DO i = 1, lonk
255     IF (t0(i,jv)>0.) THEN
256     slpmax = t0(i, jv)
257     s1max = 1.5*slpmax
258     s1new = amin1(s1max, amax1(-s1max,tx(i,jv)))
259     s2new = amin1(2.*slpmax-abs(s1new)/3., amax1(abs( &
260     s1new)-slpmax,txx(i,jv)))
261     tx(i, jv) = s1new
262     txx(i, jv) = s2new
263     txy(i, jv) = amin1(slpmax, amax1(-slpmax,txy(i,jv)))
264     txz(i, jv) = amin1(slpmax, amax1(-slpmax,txz(i,jv)))
265     ELSE
266     tx(i, jv) = 0.
267     txx(i, jv) = 0.
268     txy(i, jv) = 0.
269     txz(i, jv) = 0.
270     END IF
271     END DO
272 guez 3 END DO
273    
274 guez 81 13 CONTINUE
275    
276     ! calculate flux and moments between adjacent boxes
277     ! 1- create temporary moments/masses for partial boxes in transit
278     ! 2- reajusts moments remaining in the box
279    
280     ! flux from IP to I if U(I).lt.0
281    
282     DO i = 1, lonk - 1
283     IF (uext(i)<0.) THEN
284     fm(i) = -uext(i)*dtx
285     alf(i) = fm(i)/tm(i+1)
286     tm(i+1) = tm(i+1) - fm(i)
287     END IF
288 guez 3 END DO
289 guez 81
290     i = lonk
291     IF (uext(i)<0.) THEN
292     fm(i) = -uext(i)*dtx
293     alf(i) = fm(i)/tm(1)
294     tm(1) = tm(1) - fm(i)
295     END IF
296    
297     ! flux from I to IP if U(I).gt.0
298    
299     DO i = 1, lonk
300     IF (uext(i)>=0.) THEN
301     fm(i) = uext(i)*dtx
302     alf(i) = fm(i)/tm(i)
303     tm(i) = tm(i) - fm(i)
304     END IF
305 guez 3 END DO
306 guez 81
307     DO i = 1, lonk
308     alfq(i) = alf(i)*alf(i)
309     alf1(i) = 1. - alf(i)
310     alf1q(i) = alf1(i)*alf1(i)
311     alf2(i) = alf1(i) - alf(i)
312     alf3(i) = alf(i)*alfq(i)
313     alf4(i) = alf1(i)*alf1q(i)
314 guez 3 END DO
315    
316 guez 81 DO jv = 1, ntra
317     DO i = 1, lonk - 1
318    
319     IF (uext(i)<0.) THEN
320    
321     f0(i, jv) = alf(i)*(t0(i+1,jv)-alf1(i)*(tx(i+1, &
322     jv)-alf2(i)*txx(i+1,jv)))
323     fx(i, jv) = alfq(i)*(tx(i+1,jv)-3.*alf1(i)*txx(i+1,jv))
324     fxx(i, jv) = alf3(i)*txx(i+1, jv)
325     fy(i, jv) = alf(i)*(ty(i+1,jv)-alf1(i)*txy(i+1,jv))
326     fz(i, jv) = alf(i)*(tz(i+1,jv)-alf1(i)*txz(i+1,jv))
327     fxy(i, jv) = alfq(i)*txy(i+1, jv)
328     fxz(i, jv) = alfq(i)*txz(i+1, jv)
329     fyy(i, jv) = alf(i)*tyy(i+1, jv)
330     fyz(i, jv) = alf(i)*tyz(i+1, jv)
331     fzz(i, jv) = alf(i)*tzz(i+1, jv)
332    
333     t0(i+1, jv) = t0(i+1, jv) - f0(i, jv)
334     tx(i+1, jv) = alf1q(i)*(tx(i+1,jv)+3.*alf(i)*txx(i+1,jv))
335     txx(i+1, jv) = alf4(i)*txx(i+1, jv)
336     ty(i+1, jv) = ty(i+1, jv) - fy(i, jv)
337     tz(i+1, jv) = tz(i+1, jv) - fz(i, jv)
338     tyy(i+1, jv) = tyy(i+1, jv) - fyy(i, jv)
339     tyz(i+1, jv) = tyz(i+1, jv) - fyz(i, jv)
340     tzz(i+1, jv) = tzz(i+1, jv) - fzz(i, jv)
341     txy(i+1, jv) = alf1q(i)*txy(i+1, jv)
342     txz(i+1, jv) = alf1q(i)*txz(i+1, jv)
343    
344     END IF
345    
346     END DO
347     END DO
348    
349     i = lonk
350     IF (uext(i)<0.) THEN
351    
352     DO jv = 1, ntra
353    
354     f0(i, jv) = alf(i)*(t0(1,jv)-alf1(i)*(tx(1,jv)-alf2(i)*txx(1,jv)))
355     fx(i, jv) = alfq(i)*(tx(1,jv)-3.*alf1(i)*txx(1,jv))
356     fxx(i, jv) = alf3(i)*txx(1, jv)
357     fy(i, jv) = alf(i)*(ty(1,jv)-alf1(i)*txy(1,jv))
358     fz(i, jv) = alf(i)*(tz(1,jv)-alf1(i)*txz(1,jv))
359     fxy(i, jv) = alfq(i)*txy(1, jv)
360     fxz(i, jv) = alfq(i)*txz(1, jv)
361     fyy(i, jv) = alf(i)*tyy(1, jv)
362     fyz(i, jv) = alf(i)*tyz(1, jv)
363     fzz(i, jv) = alf(i)*tzz(1, jv)
364    
365     t0(1, jv) = t0(1, jv) - f0(i, jv)
366     tx(1, jv) = alf1q(i)*(tx(1,jv)+3.*alf(i)*txx(1,jv))
367     txx(1, jv) = alf4(i)*txx(1, jv)
368     ty(1, jv) = ty(1, jv) - fy(i, jv)
369     tz(1, jv) = tz(1, jv) - fz(i, jv)
370     tyy(1, jv) = tyy(1, jv) - fyy(i, jv)
371     tyz(1, jv) = tyz(1, jv) - fyz(i, jv)
372     tzz(1, jv) = tzz(1, jv) - fzz(i, jv)
373     txy(1, jv) = alf1q(i)*txy(1, jv)
374     txz(1, jv) = alf1q(i)*txz(1, jv)
375    
376     END DO
377    
378     END IF
379    
380     DO jv = 1, ntra
381     DO i = 1, lonk
382    
383     IF (uext(i)>=0.) THEN
384    
385     f0(i, jv) = alf(i)*(t0(i,jv)+alf1(i)*(tx(i,jv)+alf2(i)*txx(i, &
386     jv)))
387     fx(i, jv) = alfq(i)*(tx(i,jv)+3.*alf1(i)*txx(i,jv))
388     fxx(i, jv) = alf3(i)*txx(i, jv)
389     fy(i, jv) = alf(i)*(ty(i,jv)+alf1(i)*txy(i,jv))
390     fz(i, jv) = alf(i)*(tz(i,jv)+alf1(i)*txz(i,jv))
391     fxy(i, jv) = alfq(i)*txy(i, jv)
392     fxz(i, jv) = alfq(i)*txz(i, jv)
393     fyy(i, jv) = alf(i)*tyy(i, jv)
394     fyz(i, jv) = alf(i)*tyz(i, jv)
395     fzz(i, jv) = alf(i)*tzz(i, jv)
396    
397     t0(i, jv) = t0(i, jv) - f0(i, jv)
398     tx(i, jv) = alf1q(i)*(tx(i,jv)-3.*alf(i)*txx(i,jv))
399     txx(i, jv) = alf4(i)*txx(i, jv)
400     ty(i, jv) = ty(i, jv) - fy(i, jv)
401     tz(i, jv) = tz(i, jv) - fz(i, jv)
402     tyy(i, jv) = tyy(i, jv) - fyy(i, jv)
403     tyz(i, jv) = tyz(i, jv) - fyz(i, jv)
404     tzz(i, jv) = tzz(i, jv) - fzz(i, jv)
405     txy(i, jv) = alf1q(i)*txy(i, jv)
406     txz(i, jv) = alf1q(i)*txz(i, jv)
407    
408     END IF
409    
410     END DO
411     END DO
412    
413     ! puts the temporary moments Fi into appropriate neighboring boxes
414    
415     DO i = 1, lonk
416     IF (uext(i)<0.) THEN
417     tm(i) = tm(i) + fm(i)
418     alf(i) = fm(i)/tm(i)
419     END IF
420     END DO
421    
422     DO i = 1, lonk - 1
423     IF (uext(i)>=0.) THEN
424     tm(i+1) = tm(i+1) + fm(i)
425     alf(i) = fm(i)/tm(i+1)
426     END IF
427     END DO
428    
429     i = lonk
430     IF (uext(i)>=0.) THEN
431     tm(1) = tm(1) + fm(i)
432     alf(i) = fm(i)/tm(1)
433     END IF
434    
435     DO i = 1, lonk
436     alf1(i) = 1. - alf(i)
437     alfq(i) = alf(i)*alf(i)
438     alf1q(i) = alf1(i)*alf1(i)
439     alf2(i) = alf1(i) - alf(i)
440     alf3(i) = alf(i)*alf1(i)
441     END DO
442    
443     DO jv = 1, ntra
444     DO i = 1, lonk
445    
446     IF (uext(i)<0.) THEN
447    
448     temptm = -alf(i)*t0(i, jv) + alf1(i)*f0(i, jv)
449     t0(i, jv) = t0(i, jv) + f0(i, jv)
450     txx(i, jv) = alfq(i)*fxx(i, jv) + alf1q(i)*txx(i, jv) + &
451     5.*(alf3(i)*(fx(i,jv)-tx(i,jv))+alf2(i)*temptm)
452     tx(i, jv) = alf(i)*fx(i, jv) + alf1(i)*tx(i, jv) + 3.*temptm
453     txy(i, jv) = alf(i)*fxy(i, jv) + alf1(i)*txy(i, jv) + &
454     3.*(alf1(i)*fy(i,jv)-alf(i)*ty(i,jv))
455     txz(i, jv) = alf(i)*fxz(i, jv) + alf1(i)*txz(i, jv) + &
456     3.*(alf1(i)*fz(i,jv)-alf(i)*tz(i,jv))
457     ty(i, jv) = ty(i, jv) + fy(i, jv)
458     tz(i, jv) = tz(i, jv) + fz(i, jv)
459     tyy(i, jv) = tyy(i, jv) + fyy(i, jv)
460     tyz(i, jv) = tyz(i, jv) + fyz(i, jv)
461     tzz(i, jv) = tzz(i, jv) + fzz(i, jv)
462    
463     END IF
464    
465     END DO
466     END DO
467    
468     DO jv = 1, ntra
469     DO i = 1, lonk - 1
470    
471     IF (uext(i)>=0.) THEN
472    
473     temptm = alf(i)*t0(i+1, jv) - alf1(i)*f0(i, jv)
474     t0(i+1, jv) = t0(i+1, jv) + f0(i, jv)
475     txx(i+1, jv) = alfq(i)*fxx(i, jv) + alf1q(i)*txx(i+1, jv) + &
476     5.*(alf3(i)*(tx(i+1,jv)-fx(i,jv))-alf2(i)*temptm)
477     tx(i+1, jv) = alf(i)*fx(i, jv) + alf1(i)*tx(i+1, jv) + 3.*temptm
478     txy(i+1, jv) = alf(i)*fxy(i, jv) + alf1(i)*txy(i+1, jv) + &
479     3.*(alf(i)*ty(i+1,jv)-alf1(i)*fy(i,jv))
480     txz(i+1, jv) = alf(i)*fxz(i, jv) + alf1(i)*txz(i+1, jv) + &
481     3.*(alf(i)*tz(i+1,jv)-alf1(i)*fz(i,jv))
482     ty(i+1, jv) = ty(i+1, jv) + fy(i, jv)
483     tz(i+1, jv) = tz(i+1, jv) + fz(i, jv)
484     tyy(i+1, jv) = tyy(i+1, jv) + fyy(i, jv)
485     tyz(i+1, jv) = tyz(i+1, jv) + fyz(i, jv)
486     tzz(i+1, jv) = tzz(i+1, jv) + fzz(i, jv)
487    
488     END IF
489    
490     END DO
491     END DO
492    
493     i = lonk
494     IF (uext(i)>=0.) THEN
495     DO jv = 1, ntra
496     temptm = alf(i)*t0(1, jv) - alf1(i)*f0(i, jv)
497     t0(1, jv) = t0(1, jv) + f0(i, jv)
498     txx(1, jv) = alfq(i)*fxx(i, jv) + alf1q(i)*txx(1, jv) + &
499     5.*(alf3(i)*(tx(1,jv)-fx(i,jv))-alf2(i)*temptm)
500     tx(1, jv) = alf(i)*fx(i, jv) + alf1(i)*tx(1, jv) + 3.*temptm
501     txy(1, jv) = alf(i)*fxy(i, jv) + alf1(i)*txy(1, jv) + &
502     3.*(alf(i)*ty(1,jv)-alf1(i)*fy(i,jv))
503     txz(1, jv) = alf(i)*fxz(i, jv) + alf1(i)*txz(1, jv) + &
504     3.*(alf(i)*tz(1,jv)-alf1(i)*fz(i,jv))
505     ty(1, jv) = ty(1, jv) + fy(i, jv)
506     tz(1, jv) = tz(1, jv) + fz(i, jv)
507     tyy(1, jv) = tyy(1, jv) + fyy(i, jv)
508     tyz(1, jv) = tyz(1, jv) + fyz(i, jv)
509     tzz(1, jv) = tzz(1, jv) + fzz(i, jv)
510     END DO
511     END IF
512    
513     ! retour aux mailles d'origine (passage des Tij aux Sij)
514    
515     IF (numk>1) THEN
516    
517     DO i2 = 1, numk
518    
519     DO i = 1, lonk
520    
521     i3 = i2 + (i-1)*numk
522     sm(i3, k, l) = smnew(i3)
523     alf(i) = smnew(i3)/tm(i)
524     tm(i) = tm(i) - smnew(i3)
525    
526     alfq(i) = alf(i)*alf(i)
527     alf1(i) = 1. - alf(i)
528     alf1q(i) = alf1(i)*alf1(i)
529     alf2(i) = alf1(i) - alf(i)
530     alf3(i) = alf(i)*alfq(i)
531     alf4(i) = alf1(i)*alf1q(i)
532    
533     END DO
534    
535     DO jv = 1, ntra
536     DO i = 1, lonk
537    
538     i3 = i2 + (i-1)*numk
539     s0(i3, k, l, jv) = alf(i)*(t0(i,jv)-alf1(i)*(tx(i, &
540     jv)-alf2(i)*txx(i,jv)))
541     ssx(i3, k, l, jv) = alfq(i)*(tx(i,jv)-3.*alf1(i)*txx(i,jv))
542     ssxx(i3, k, l, jv) = alf3(i)*txx(i, jv)
543     sy(i3, k, l, jv) = alf(i)*(ty(i,jv)-alf1(i)*txy(i,jv))
544     sz(i3, k, l, jv) = alf(i)*(tz(i,jv)-alf1(i)*txz(i,jv))
545     ssxy(i3, k, l, jv) = alfq(i)*txy(i, jv)
546     ssxz(i3, k, l, jv) = alfq(i)*txz(i, jv)
547     syy(i3, k, l, jv) = alf(i)*tyy(i, jv)
548     syz(i3, k, l, jv) = alf(i)*tyz(i, jv)
549     szz(i3, k, l, jv) = alf(i)*tzz(i, jv)
550    
551     ! reajusts moments remaining in the box
552    
553     t0(i, jv) = t0(i, jv) - s0(i3, k, l, jv)
554     tx(i, jv) = alf1q(i)*(tx(i,jv)+3.*alf(i)*txx(i,jv))
555     txx(i, jv) = alf4(i)*txx(i, jv)
556     ty(i, jv) = ty(i, jv) - sy(i3, k, l, jv)
557     tz(i, jv) = tz(i, jv) - sz(i3, k, l, jv)
558     tyy(i, jv) = tyy(i, jv) - syy(i3, k, l, jv)
559     tyz(i, jv) = tyz(i, jv) - syz(i3, k, l, jv)
560     tzz(i, jv) = tzz(i, jv) - szz(i3, k, l, jv)
561     txy(i, jv) = alf1q(i)*txy(i, jv)
562     txz(i, jv) = alf1q(i)*txz(i, jv)
563    
564     END DO
565     END DO
566    
567     END DO
568    
569     ELSE
570    
571     DO i = 1, lon
572     sm(i, k, l) = tm(i)
573     END DO
574     DO jv = 1, ntra
575     DO i = 1, lon
576     s0(i, k, l, jv) = t0(i, jv)
577     ssx(i, k, l, jv) = tx(i, jv)
578     sy(i, k, l, jv) = ty(i, jv)
579     sz(i, k, l, jv) = tz(i, jv)
580     ssxx(i, k, l, jv) = txx(i, jv)
581     ssxy(i, k, l, jv) = txy(i, jv)
582     ssxz(i, k, l, jv) = txz(i, jv)
583     syy(i, k, l, jv) = tyy(i, jv)
584     syz(i, k, l, jv) = tyz(i, jv)
585     szz(i, k, l, jv) = tzz(i, jv)
586     END DO
587     END DO
588    
589     END IF
590    
591     END DO
592     END DO
593    
594     ! ---------- bouclage cyclique
595    
596     DO l = 1, llm
597     DO j = 1, jjp1
598     sm(iip1, j, l) = sm(1, j, l)
599     s0(iip1, j, l, ntra) = s0(1, j, l, ntra)
600     ssx(iip1, j, l, ntra) = ssx(1, j, l, ntra)
601     sy(iip1, j, l, ntra) = sy(1, j, l, ntra)
602     sz(iip1, j, l, ntra) = sz(1, j, l, ntra)
603     END DO
604     END DO
605    
606     ! ----------- qqtite totale de traceur dans tte l'atmosphere
607     DO l = 1, llm
608     DO j = 1, jjp1
609     DO i = 1, iim
610     sqf = sqf + s0(i, j, l, ntra)
611     END DO
612     END DO
613     END DO
614    
615     PRINT *, '------ DIAG DANS ADVX2 - SORTIE -----'
616     PRINT *, 'sqf=', sqf
617     ! -------------------------------------------------------------
618     RETURN
619     END SUBROUTINE advxp

  ViewVC Help
Powered by ViewVC 1.1.21