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

Annotation of /trunk/dyn3d/advyp.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 82 - (hide annotations)
Wed Mar 5 14:57:53 2014 UTC (10 years, 2 months ago) by guez
File size: 19953 byte(s)
Changed all ".f90" suffixes to ".f".
1 guez 3
2 guez 81 ! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/advyp.F,v 1.1.1.1 2004/05/19
3     ! 12:53:06 lmdzadmin Exp $
4 guez 3
5 guez 81 SUBROUTINE advyp(limit, dty, pbarv, sm, s0, ssx, sy, sz, ssxx, ssxy, ssxz, &
6     syy, syz, szz, ntra)
7     USE dimens_m
8     USE comconst
9     USE paramet_m
10     USE disvert_m
11     USE comgeom
12     IMPLICIT NONE
13     ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
14     ! C
15     ! second-order moments (SOM) advection of tracer in Y direction C
16     ! C
17     ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
18     ! C
19     ! Source : Pascal Simon ( Meteo, CNRM ) C
20     ! Adaptation : A.A. (LGGE) C
21     ! Derniere Modif : 19/10/95 LAST
22     ! C
23     ! sont les arguments d'entree pour le s-pg C
24     ! C
25     ! argument de sortie du s-pg C
26     ! C
27     ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
28     ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
29 guez 3
30 guez 81 ! Rem : Probleme aux poles il faut reecrire ce cas specifique
31     ! Attention au sens de l'indexation
32 guez 3
33 guez 81 ! parametres principaux du modele
34 guez 3
35    
36    
37 guez 81 ! Arguments :
38     ! ----------
39     ! dty : frequence fictive d'appel du transport
40     ! parbu,pbarv : flux de masse en x et y en Pa.m2.s-1
41 guez 3
42 guez 81 INTEGER lon, lat, niv
43     INTEGER i, j, jv, k, kp, l
44     INTEGER ntra
45     ! PARAMETER (ntra = 1)
46 guez 3
47 guez 81 REAL dty
48     REAL, INTENT (IN) :: pbarv(iip1, jjm, llm)
49 guez 3
50 guez 81 ! moments: SM total mass in each grid box
51     ! S0 mass of tracer in each grid box
52     ! Si 1rst order moment in i direction
53    
54     REAL sm(iip1, jjp1, llm), s0(iip1, jjp1, llm, ntra)
55     REAL ssx(iip1, jjp1, llm, ntra), sy(iip1, jjp1, llm, ntra), &
56     sz(iip1, jjp1, llm, ntra), ssxx(iip1, jjp1, llm, ntra), &
57     ssxy(iip1, jjp1, llm, ntra), ssxz(iip1, jjp1, llm, ntra), &
58     syy(iip1, jjp1, llm, ntra), syz(iip1, jjp1, llm, ntra), &
59     szz(iip1, jjp1, llm, ntra)
60    
61     ! Local :
62     ! -------
63    
64     ! mass fluxes across the boundaries (UGRI,VGRI,WGRI)
65     ! mass fluxes in kg
66     ! declaration :
67    
68     REAL vgri(iip1, 0:jjp1, llm)
69    
70     ! Rem : UGRI et WGRI ne sont pas utilises dans
71     ! cette subroutine ( advection en y uniquement )
72     ! Rem 2 :le dimensionnement de VGRI depend de celui de pbarv
73    
74     ! the moments F are similarly defined and used as temporary
75     ! storage for portions of the grid boxes in transit
76    
77     ! the moments Fij are used as temporary storage for
78     ! portions of the grid boxes in transit at the current level
79    
80     ! work arrays
81    
82    
83     REAL f0(iim, 0:jjp1, ntra), fm(iim, 0:jjp1)
84     REAL fx(iim, jjm, ntra), fy(iim, jjm, ntra)
85     REAL fz(iim, jjm, ntra)
86     REAL fxx(iim, jjm, ntra), fxy(iim, jjm, ntra)
87     REAL fxz(iim, jjm, ntra), fyy(iim, jjm, ntra)
88     REAL fyz(iim, jjm, ntra), fzz(iim, jjm, ntra)
89     REAL s00(ntra)
90     REAL sm0 ! Just temporal variable
91    
92     ! work arrays
93    
94     REAL alf(iim, 0:jjp1), alf1(iim, 0:jjp1)
95     REAL alfq(iim, 0:jjp1), alf1q(iim, 0:jjp1)
96     REAL alf2(iim, 0:jjp1), alf3(iim, 0:jjp1)
97     REAL alf4(iim, 0:jjp1)
98     REAL temptm ! Just temporal variable
99     REAL slpmax, s1max, s1new, s2new
100    
101     ! Special pour poles
102    
103     REAL sbms, sfms, sfzs, sbmn, sfmn, sfzn
104     REAL sns0(ntra), snsz(ntra), snsm
105     REAL qy1(iim, llm, ntra), qylat(iim, llm, ntra)
106     REAL cx1(llm, ntra), cxlat(llm, ntra)
107     REAL cy1(llm, ntra), cylat(llm, ntra)
108     REAL z1(iim), zcos(iim), zsin(iim)
109     REAL ssum
110     EXTERNAL ssum
111    
112     REAL sqi, sqf
113     LOGICAL limit
114    
115     lon = iim ! rem : Il est possible qu'un pbl. arrive ici
116     lat = jjp1 ! a cause des dim. differentes entre les
117     niv = llm ! tab. S et VGRI
118    
119     ! -----------------------------------------------------------------
120     ! initialisations
121    
122     sbms = 0.
123     sfms = 0.
124     sfzs = 0.
125     sbmn = 0.
126     sfmn = 0.
127     sfzn = 0.
128    
129     ! -----------------------------------------------------------------
130     ! *** Test : diag de la qtite totale de traceur dans
131     ! l'atmosphere avant l'advection en Y
132    
133     sqi = 0.
134     sqf = 0.
135    
136     DO l = 1, llm
137     DO j = 1, jjp1
138     DO i = 1, iim
139     sqi = sqi + s0(i, j, l, ntra)
140 guez 3 END DO
141 guez 81 END DO
142     END DO
143     PRINT *, '---------- DIAG DANS ADVY - ENTREE --------'
144     PRINT *, 'sqi=', sqi
145 guez 3
146 guez 81 ! -----------------------------------------------------------------
147     ! Interface : adaptation nouveau modele
148     ! -------------------------------------
149 guez 3
150 guez 81 ! Conversion des flux de masses en kg
151     ! -AA 20/10/94 le signe -1 est necessaire car indexation opposee
152 guez 3
153 guez 81 DO l = 1, llm
154     DO j = 1, jjm
155     DO i = 1, iip1
156     vgri(i, j, llm+1-l) = -1.*pbarv(i, j, l)
157     END DO
158     END DO
159     END DO
160 guez 3
161 guez 81 ! AA Initialisation de flux fictifs aux bords sup. des boites pol.
162 guez 3
163 guez 81 DO l = 1, llm
164     DO i = 1, iip1
165     vgri(i, 0, l) = 0.
166     vgri(i, jjp1, l) = 0.
167     END DO
168     END DO
169    
170     ! ----------------- START HERE -----------------------
171     ! boucle sur les niveaux
172    
173     DO l = 1, niv
174    
175     ! place limits on appropriate moments before transport
176     ! (if flux-limiting is to be applied)
177    
178     IF (.NOT. limit) GO TO 11
179    
180     DO jv = 1, ntra
181     DO k = 1, lat
182     DO i = 1, lon
183     IF (s0(i,k,l,jv)>0.) THEN
184     slpmax = amax1(s0(i,k,l,jv), 0.)
185     s1max = 1.5*slpmax
186     s1new = amin1(s1max, amax1(-s1max,sy(i,k,l,jv)))
187     s2new = amin1(2.*slpmax-abs(s1new)/3., amax1(abs( &
188     s1new)-slpmax,syy(i,k,l,jv)))
189     sy(i, k, l, jv) = s1new
190     syy(i, k, l, jv) = s2new
191     ssxy(i, k, l, jv) = amin1(slpmax, amax1(-slpmax,ssxy(i,k,l,jv)))
192     syz(i, k, l, jv) = amin1(slpmax, amax1(-slpmax,syz(i,k,l,jv)))
193     ELSE
194     sy(i, k, l, jv) = 0.
195     syy(i, k, l, jv) = 0.
196     ssxy(i, k, l, jv) = 0.
197     syz(i, k, l, jv) = 0.
198     END IF
199     END DO
200 guez 3 END DO
201 guez 81 END DO
202 guez 3
203 guez 81 11 CONTINUE
204 guez 3
205 guez 81 ! le flux a travers le pole Nord est traite separement
206 guez 3
207 guez 81 sm0 = 0.
208     DO jv = 1, ntra
209     s00(jv) = 0.
210     END DO
211 guez 3
212 guez 81 DO i = 1, lon
213 guez 3
214 guez 81 IF (vgri(i,0,l)<=0.) THEN
215     fm(i, 0) = -vgri(i, 0, l)*dty
216     alf(i, 0) = fm(i, 0)/sm(i, 1, l)
217     sm(i, 1, l) = sm(i, 1, l) - fm(i, 0)
218     sm0 = sm0 + fm(i, 0)
219     END IF
220 guez 3
221 guez 81 alfq(i, 0) = alf(i, 0)*alf(i, 0)
222     alf1(i, 0) = 1. - alf(i, 0)
223     alf1q(i, 0) = alf1(i, 0)*alf1(i, 0)
224     alf2(i, 0) = alf1(i, 0) - alf(i, 0)
225     alf3(i, 0) = alf(i, 0)*alfq(i, 0)
226     alf4(i, 0) = alf1(i, 0)*alf1q(i, 0)
227 guez 3
228 guez 81 END DO
229     ! print*,'ADVYP 21'
230 guez 3
231 guez 81 DO jv = 1, ntra
232     DO i = 1, lon
233 guez 3
234 guez 81 IF (vgri(i,0,l)<=0.) THEN
235 guez 3
236 guez 81 f0(i, 0, jv) = alf(i, 0)*(s0(i,1,l,jv)-alf1(i,0)*(sy(i,1,l, &
237     jv)-alf2(i,0)*syy(i,1,l,jv)))
238 guez 3
239 guez 81 s00(jv) = s00(jv) + f0(i, 0, jv)
240     s0(i, 1, l, jv) = s0(i, 1, l, jv) - f0(i, 0, jv)
241     sy(i, 1, l, jv) = alf1q(i, 0)*(sy(i,1,l,jv)+3.*alf(i,0)*syy(i,1,l, &
242     jv))
243     syy(i, 1, l, jv) = alf4(i, 0)*syy(i, 1, l, jv)
244     ssx(i, 1, l, jv) = alf1(i, 0)*(ssx(i,1,l,jv)+alf(i,0)*ssxy(i,1,l,jv &
245     ))
246     sz(i, 1, l, jv) = alf1(i, 0)*(sz(i,1,l,jv)+alf(i,0)*ssxz(i,1,l,jv))
247     ssxx(i, 1, l, jv) = alf1(i, 0)*ssxx(i, 1, l, jv)
248     ssxz(i, 1, l, jv) = alf1(i, 0)*ssxz(i, 1, l, jv)
249     szz(i, 1, l, jv) = alf1(i, 0)*szz(i, 1, l, jv)
250     ssxy(i, 1, l, jv) = alf1q(i, 0)*ssxy(i, 1, l, jv)
251     syz(i, 1, l, jv) = alf1q(i, 0)*syz(i, 1, l, jv)
252 guez 3
253 guez 81 END IF
254 guez 3
255 guez 81 END DO
256     END DO
257 guez 3
258 guez 81 DO i = 1, lon
259     IF (vgri(i,0,l)>0.) THEN
260     fm(i, 0) = vgri(i, 0, l)*dty
261     alf(i, 0) = fm(i, 0)/sm0
262     END IF
263     END DO
264 guez 3
265 guez 81 DO jv = 1, ntra
266     DO i = 1, lon
267     IF (vgri(i,0,l)>0.) THEN
268     f0(i, 0, jv) = alf(i, 0)*s00(jv)
269     END IF
270     END DO
271     END DO
272    
273     ! puts the temporary moments Fi into appropriate neighboring boxes
274    
275     ! print*,'av ADVYP 25'
276     DO i = 1, lon
277    
278     IF (vgri(i,0,l)>0.) THEN
279     sm(i, 1, l) = sm(i, 1, l) + fm(i, 0)
280     alf(i, 0) = fm(i, 0)/sm(i, 1, l)
281     END IF
282    
283     alfq(i, 0) = alf(i, 0)*alf(i, 0)
284     alf1(i, 0) = 1. - alf(i, 0)
285     alf1q(i, 0) = alf1(i, 0)*alf1(i, 0)
286     alf2(i, 0) = alf1(i, 0) - alf(i, 0)
287     alf3(i, 0) = alf1(i, 0)*alf(i, 0)
288    
289     END DO
290     ! print*,'av ADVYP 25'
291    
292     DO jv = 1, ntra
293     DO i = 1, lon
294    
295     IF (vgri(i,0,l)>0.) THEN
296    
297     temptm = alf(i, 0)*s0(i, 1, l, jv) - alf1(i, 0)*f0(i, 0, jv)
298     s0(i, 1, l, jv) = s0(i, 1, l, jv) + f0(i, 0, jv)
299     syy(i, 1, l, jv) = alf1q(i, 0)*syy(i, 1, l, jv) + &
300     5.*(alf3(i,0)*sy(i,1,l,jv)-alf2(i,0)*temptm)
301     sy(i, 1, l, jv) = alf1(i, 0)*sy(i, 1, l, jv) + 3.*temptm
302     ssxy(i, 1, l, jv) = alf1(i, 0)*ssxy(i, 1, l, jv) + &
303     3.*alf(i, 0)*ssx(i, 1, l, jv)
304     syz(i, 1, l, jv) = alf1(i, 0)*syz(i, 1, l, jv) + &
305     3.*alf(i, 0)*sz(i, 1, l, jv)
306    
307     END IF
308    
309     END DO
310     END DO
311    
312     ! calculate flux and moments between adjacent boxes
313     ! 1- create temporary moments/masses for partial boxes in transit
314     ! 2- reajusts moments remaining in the box
315    
316     ! flux from KP to K if V(K).lt.0 and from K to KP if V(K).gt.0
317    
318     ! print*,'av ADVYP 30'
319     DO k = 1, lat - 1
320     kp = k + 1
321     DO i = 1, lon
322    
323     IF (vgri(i,k,l)<0.) THEN
324     fm(i, k) = -vgri(i, k, l)*dty
325     alf(i, k) = fm(i, k)/sm(i, kp, l)
326     sm(i, kp, l) = sm(i, kp, l) - fm(i, k)
327     ELSE
328     fm(i, k) = vgri(i, k, l)*dty
329     alf(i, k) = fm(i, k)/sm(i, k, l)
330     sm(i, k, l) = sm(i, k, l) - fm(i, k)
331     END IF
332    
333     alfq(i, k) = alf(i, k)*alf(i, k)
334     alf1(i, k) = 1. - alf(i, k)
335     alf1q(i, k) = alf1(i, k)*alf1(i, k)
336     alf2(i, k) = alf1(i, k) - alf(i, k)
337     alf3(i, k) = alf(i, k)*alfq(i, k)
338     alf4(i, k) = alf1(i, k)*alf1q(i, k)
339    
340     END DO
341     END DO
342     ! print*,'ap ADVYP 30'
343    
344     DO jv = 1, ntra
345     DO k = 1, lat - 1
346     kp = k + 1
347     DO i = 1, lon
348    
349     IF (vgri(i,k,l)<0.) THEN
350    
351     f0(i, k, jv) = alf(i, k)*(s0(i,kp,l,jv)-alf1(i,k)*(sy(i,kp,l, &
352     jv)-alf2(i,k)*syy(i,kp,l,jv)))
353     fy(i, k, jv) = alfq(i, k)*(sy(i,kp,l,jv)-3.*alf1(i,k)*syy(i,kp,l, &
354     jv))
355     fyy(i, k, jv) = alf3(i, k)*syy(i, kp, l, jv)
356     fx(i, k, jv) = alf(i, k)*(ssx(i,kp,l,jv)-alf1(i,k)*ssxy(i,kp,l,jv &
357     ))
358     fz(i, k, jv) = alf(i, k)*(sz(i,kp,l,jv)-alf1(i,k)*syz(i,kp,l,jv))
359     fxy(i, k, jv) = alfq(i, k)*ssxy(i, kp, l, jv)
360     fyz(i, k, jv) = alfq(i, k)*syz(i, kp, l, jv)
361     fxx(i, k, jv) = alf(i, k)*ssxx(i, kp, l, jv)
362     fxz(i, k, jv) = alf(i, k)*ssxz(i, kp, l, jv)
363     fzz(i, k, jv) = alf(i, k)*szz(i, kp, l, jv)
364    
365     s0(i, kp, l, jv) = s0(i, kp, l, jv) - f0(i, k, jv)
366     sy(i, kp, l, jv) = alf1q(i, k)*(sy(i,kp,l,jv)+3.*alf(i,k)*syy(i, &
367     kp,l,jv))
368     syy(i, kp, l, jv) = alf4(i, k)*syy(i, kp, l, jv)
369     ssx(i, kp, l, jv) = ssx(i, kp, l, jv) - fx(i, k, jv)
370     sz(i, kp, l, jv) = sz(i, kp, l, jv) - fz(i, k, jv)
371     ssxx(i, kp, l, jv) = ssxx(i, kp, l, jv) - fxx(i, k, jv)
372     ssxz(i, kp, l, jv) = ssxz(i, kp, l, jv) - fxz(i, k, jv)
373     szz(i, kp, l, jv) = szz(i, kp, l, jv) - fzz(i, k, jv)
374     ssxy(i, kp, l, jv) = alf1q(i, k)*ssxy(i, kp, l, jv)
375     syz(i, kp, l, jv) = alf1q(i, k)*syz(i, kp, l, jv)
376    
377     ELSE
378    
379     f0(i, k, jv) = alf(i, k)*(s0(i,k,l,jv)+alf1(i,k)*(sy(i,k,l, &
380     jv)+alf2(i,k)*syy(i,k,l,jv)))
381     fy(i, k, jv) = alfq(i, k)*(sy(i,k,l,jv)+3.*alf1(i,k)*syy(i,k,l,jv &
382     ))
383     fyy(i, k, jv) = alf3(i, k)*syy(i, k, l, jv)
384     fx(i, k, jv) = alf(i, k)*(ssx(i,k,l,jv)+alf1(i,k)*ssxy(i,k,l,jv))
385     fz(i, k, jv) = alf(i, k)*(sz(i,k,l,jv)+alf1(i,k)*syz(i,k,l,jv))
386     fxy(i, k, jv) = alfq(i, k)*ssxy(i, k, l, jv)
387     fyz(i, k, jv) = alfq(i, k)*syz(i, k, l, jv)
388     fxx(i, k, jv) = alf(i, k)*ssxx(i, k, l, jv)
389     fxz(i, k, jv) = alf(i, k)*ssxz(i, k, l, jv)
390     fzz(i, k, jv) = alf(i, k)*szz(i, k, l, jv)
391    
392     s0(i, k, l, jv) = s0(i, k, l, jv) - f0(i, k, jv)
393     sy(i, k, l, jv) = alf1q(i, k)*(sy(i,k,l,jv)-3.*alf(i,k)*syy(i,k,l &
394     ,jv))
395     syy(i, k, l, jv) = alf4(i, k)*syy(i, k, l, jv)
396     ssx(i, k, l, jv) = ssx(i, k, l, jv) - fx(i, k, jv)
397     sz(i, k, l, jv) = sz(i, k, l, jv) - fz(i, k, jv)
398     ssxx(i, k, l, jv) = ssxx(i, k, l, jv) - fxx(i, k, jv)
399     ssxz(i, k, l, jv) = ssxz(i, k, l, jv) - fxz(i, k, jv)
400     szz(i, k, l, jv) = szz(i, k, l, jv) - fzz(i, k, jv)
401     ssxy(i, k, l, jv) = alf1q(i, k)*ssxy(i, k, l, jv)
402     syz(i, k, l, jv) = alf1q(i, k)*syz(i, k, l, jv)
403    
404     END IF
405    
406     END DO
407     END DO
408     END DO
409     ! print*,'ap ADVYP 31'
410    
411     ! puts the temporary moments Fi into appropriate neighboring boxes
412    
413     DO k = 1, lat - 1
414     kp = k + 1
415     DO i = 1, lon
416    
417     IF (vgri(i,k,l)<0.) THEN
418     sm(i, k, l) = sm(i, k, l) + fm(i, k)
419     alf(i, k) = fm(i, k)/sm(i, k, l)
420     ELSE
421     sm(i, kp, l) = sm(i, kp, l) + fm(i, k)
422     alf(i, k) = fm(i, k)/sm(i, kp, l)
423     END IF
424    
425     alfq(i, k) = alf(i, k)*alf(i, k)
426     alf1(i, k) = 1. - alf(i, k)
427     alf1q(i, k) = alf1(i, k)*alf1(i, k)
428     alf2(i, k) = alf1(i, k) - alf(i, k)
429     alf3(i, k) = alf1(i, k)*alf(i, k)
430    
431     END DO
432     END DO
433     ! print*,'ap ADVYP 32'
434    
435     DO jv = 1, ntra
436     DO k = 1, lat - 1
437     kp = k + 1
438     DO i = 1, lon
439    
440     IF (vgri(i,k,l)<0.) THEN
441    
442     temptm = -alf(i, k)*s0(i, k, l, jv) + alf1(i, k)*f0(i, k, jv)
443     s0(i, k, l, jv) = s0(i, k, l, jv) + f0(i, k, jv)
444     syy(i, k, l, jv) = alfq(i, k)*fyy(i, k, jv) + &
445     alf1q(i, k)*syy(i, k, l, jv) + 5.*(alf3(i,k)*(fy(i,k,jv)-sy(i, &
446     k,l,jv))+alf2(i,k)*temptm)
447     sy(i, k, l, jv) = alf(i, k)*fy(i, k, jv) + &
448     alf1(i, k)*sy(i, k, l, jv) + 3.*temptm
449     ssxy(i, k, l, jv) = alf(i, k)*fxy(i, k, jv) + &
450     alf1(i, k)*ssxy(i, k, l, jv) + 3.*(alf1(i,k)*fx(i,k,jv)-alf(i,k &
451     )*ssx(i,k,l,jv))
452     syz(i, k, l, jv) = alf(i, k)*fyz(i, k, jv) + &
453     alf1(i, k)*syz(i, k, l, jv) + 3.*(alf1(i,k)*fz(i,k,jv)-alf(i,k) &
454     *sz(i,k,l,jv))
455     ssx(i, k, l, jv) = ssx(i, k, l, jv) + fx(i, k, jv)
456     sz(i, k, l, jv) = sz(i, k, l, jv) + fz(i, k, jv)
457     ssxx(i, k, l, jv) = ssxx(i, k, l, jv) + fxx(i, k, jv)
458     ssxz(i, k, l, jv) = ssxz(i, k, l, jv) + fxz(i, k, jv)
459     szz(i, k, l, jv) = szz(i, k, l, jv) + fzz(i, k, jv)
460    
461     ELSE
462    
463     temptm = alf(i, k)*s0(i, kp, l, jv) - alf1(i, k)*f0(i, k, jv)
464     s0(i, kp, l, jv) = s0(i, kp, l, jv) + f0(i, k, jv)
465     syy(i, kp, l, jv) = alfq(i, k)*fyy(i, k, jv) + &
466     alf1q(i, k)*syy(i, kp, l, jv) + 5.*(alf3(i,k)*(sy(i,kp,l, &
467     jv)-fy(i,k,jv))-alf2(i,k)*temptm)
468     sy(i, kp, l, jv) = alf(i, k)*fy(i, k, jv) + &
469     alf1(i, k)*sy(i, kp, l, jv) + 3.*temptm
470     ssxy(i, kp, l, jv) = alf(i, k)*fxy(i, k, jv) + &
471     alf1(i, k)*ssxy(i, kp, l, jv) + 3.*(alf(i,k)*ssx(i,kp,l,jv)- &
472     alf1(i,k)*fx(i,k,jv))
473     syz(i, kp, l, jv) = alf(i, k)*fyz(i, k, jv) + &
474     alf1(i, k)*syz(i, kp, l, jv) + 3.*(alf(i,k)*sz(i,kp,l,jv)-alf1( &
475     i,k)*fz(i,k,jv))
476     ssx(i, kp, l, jv) = ssx(i, kp, l, jv) + fx(i, k, jv)
477     sz(i, kp, l, jv) = sz(i, kp, l, jv) + fz(i, k, jv)
478     ssxx(i, kp, l, jv) = ssxx(i, kp, l, jv) + fxx(i, k, jv)
479     ssxz(i, kp, l, jv) = ssxz(i, kp, l, jv) + fxz(i, k, jv)
480     szz(i, kp, l, jv) = szz(i, kp, l, jv) + fzz(i, k, jv)
481    
482     END IF
483    
484     END DO
485     END DO
486     END DO
487     ! print*,'ap ADVYP 33'
488    
489     ! traitement special pour le pole Sud (idem pole Nord)
490    
491     k = lat
492    
493     sm0 = 0.
494     DO jv = 1, ntra
495     s00(jv) = 0.
496     END DO
497    
498     DO i = 1, lon
499    
500     IF (vgri(i,k,l)>=0.) THEN
501     fm(i, k) = vgri(i, k, l)*dty
502     alf(i, k) = fm(i, k)/sm(i, k, l)
503     sm(i, k, l) = sm(i, k, l) - fm(i, k)
504     sm0 = sm0 + fm(i, k)
505     END IF
506    
507     alfq(i, k) = alf(i, k)*alf(i, k)
508     alf1(i, k) = 1. - alf(i, k)
509     alf1q(i, k) = alf1(i, k)*alf1(i, k)
510     alf2(i, k) = alf1(i, k) - alf(i, k)
511     alf3(i, k) = alf(i, k)*alfq(i, k)
512     alf4(i, k) = alf1(i, k)*alf1q(i, k)
513    
514     END DO
515     ! print*,'ap ADVYP 41'
516    
517     DO jv = 1, ntra
518     DO i = 1, lon
519    
520     IF (vgri(i,k,l)>=0.) THEN
521     f0(i, k, jv) = alf(i, k)*(s0(i,k,l,jv)+alf1(i,k)*(sy(i,k,l, &
522     jv)+alf2(i,k)*syy(i,k,l,jv)))
523     s00(jv) = s00(jv) + f0(i, k, jv)
524    
525     s0(i, k, l, jv) = s0(i, k, l, jv) - f0(i, k, jv)
526     sy(i, k, l, jv) = alf1q(i, k)*(sy(i,k,l,jv)-3.*alf(i,k)*syy(i,k,l, &
527     jv))
528     syy(i, k, l, jv) = alf4(i, k)*syy(i, k, l, jv)
529     ssx(i, k, l, jv) = alf1(i, k)*(ssx(i,k,l,jv)-alf(i,k)*ssxy(i,k,l,jv &
530     ))
531     sz(i, k, l, jv) = alf1(i, k)*(sz(i,k,l,jv)-alf(i,k)*syz(i,k,l,jv))
532     ssxx(i, k, l, jv) = alf1(i, k)*ssxx(i, k, l, jv)
533     ssxz(i, k, l, jv) = alf1(i, k)*ssxz(i, k, l, jv)
534     szz(i, k, l, jv) = alf1(i, k)*szz(i, k, l, jv)
535     ssxy(i, k, l, jv) = alf1q(i, k)*ssxy(i, k, l, jv)
536     syz(i, k, l, jv) = alf1q(i, k)*syz(i, k, l, jv)
537     END IF
538    
539     END DO
540     END DO
541     ! print*,'ap ADVYP 42'
542    
543     DO i = 1, lon
544     IF (vgri(i,k,l)<0.) THEN
545     fm(i, k) = -vgri(i, k, l)*dty
546     alf(i, k) = fm(i, k)/sm0
547     END IF
548     END DO
549     ! print*,'ap ADVYP 43'
550    
551     DO jv = 1, ntra
552     DO i = 1, lon
553     IF (vgri(i,k,l)<0.) THEN
554     f0(i, k, jv) = alf(i, k)*s00(jv)
555     END IF
556     END DO
557     END DO
558    
559     ! puts the temporary moments Fi into appropriate neighboring boxes
560    
561     DO i = 1, lon
562    
563     IF (vgri(i,k,l)<0.) THEN
564     sm(i, k, l) = sm(i, k, l) + fm(i, k)
565     alf(i, k) = fm(i, k)/sm(i, k, l)
566     END IF
567    
568     alfq(i, k) = alf(i, k)*alf(i, k)
569     alf1(i, k) = 1. - alf(i, k)
570     alf1q(i, k) = alf1(i, k)*alf1(i, k)
571     alf2(i, k) = alf1(i, k) - alf(i, k)
572     alf3(i, k) = alf1(i, k)*alf(i, k)
573    
574     END DO
575     ! print*,'ap ADVYP 45'
576    
577     DO jv = 1, ntra
578     DO i = 1, lon
579    
580     IF (vgri(i,k,l)<0.) THEN
581    
582     temptm = -alf(i, k)*s0(i, k, l, jv) + alf1(i, k)*f0(i, k, jv)
583     s0(i, k, l, jv) = s0(i, k, l, jv) + f0(i, k, jv)
584     syy(i, k, l, jv) = alf1q(i, k)*syy(i, k, l, jv) + &
585     5.*(-alf3(i,k)*sy(i,k,l,jv)+alf2(i,k)*temptm)
586     sy(i, k, l, jv) = alf1(i, k)*sy(i, k, l, jv) + 3.*temptm
587     ssxy(i, k, l, jv) = alf1(i, k)*ssxy(i, k, l, jv) - &
588     3.*alf(i, k)*ssx(i, k, l, jv)
589     syz(i, k, l, jv) = alf1(i, k)*syz(i, k, l, jv) - &
590     3.*alf(i, k)*sz(i, k, l, jv)
591    
592     END IF
593    
594     END DO
595     END DO
596     ! print*,'ap ADVYP 46'
597    
598     END DO
599    
600     ! --------------------------------------------------
601     ! bouclage cyclique horizontal .
602    
603     DO l = 1, llm
604     DO jv = 1, ntra
605     DO j = 1, jjp1
606     sm(iip1, j, l) = sm(1, j, l)
607     s0(iip1, j, l, jv) = s0(1, j, l, jv)
608     ssx(iip1, j, l, jv) = ssx(1, j, l, jv)
609     sy(iip1, j, l, jv) = sy(1, j, l, jv)
610     sz(iip1, j, l, jv) = sz(1, j, l, jv)
611     END DO
612     END DO
613     END DO
614    
615     ! -------------------------------------------------------------------
616     ! *** Test negativite:
617    
618     ! DO jv = 1,ntra
619     ! DO l = 1,llm
620     ! DO j = 1,jjp1
621     ! DO i = 1,iip1
622     ! IF (s0( i,j,l,jv ).lt.0.) THEN
623     ! PRINT*, '------ S0 < 0 en FIN ADVYP ---'
624     ! PRINT*, 'S0(',i,j,l,jv,')=', S0(i,j,l,jv)
625     ! c STOP
626     ! ENDIF
627     ! ENDDO
628     ! ENDDO
629     ! ENDDO
630     ! ENDDO
631    
632    
633     ! -------------------------------------------------------------------
634     ! *** Test : diag de la qtite totale de traceur dans
635     ! l'atmosphere avant l'advection en Y
636    
637     DO l = 1, llm
638     DO j = 1, jjp1
639     DO i = 1, iim
640     sqf = sqf + s0(i, j, l, ntra)
641     END DO
642     END DO
643     END DO
644     PRINT *, '---------- DIAG DANS ADVY - SORTIE --------'
645     PRINT *, 'sqf=', sqf
646     ! print*,'ap ADVYP fin'
647    
648     ! -----------------------------------------------------------------
649    
650     RETURN
651     END SUBROUTINE advyp
652    
653    
654    
655    
656    
657    
658    
659    
660    
661    
662    
663    

  ViewVC Help
Powered by ViewVC 1.1.21