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

Contents of /trunk/Sources/dyn3d/advyp.f

Parent Directory Parent Directory | Revision Log Revision Log

Revision 178 - (show annotations)
Fri Mar 11 18:47:26 2016 UTC (8 years, 1 month ago) by guez
File size: 19445 byte(s) 
Moved variables date0, deltat, datasz_max, ncvar_ids, point, buff_pos,
buffer, regular from module histcom_var to modules where they are
defined.

Removed procedure ioipslmpp, useless for a sequential program.

Added argument datasz_max to histwrite_real (to avoid circular
dependency with histwrite).

Removed useless variables and computations everywhere.

Changed real litteral constants from default kind to double precision
in lwb, lwu, lwvn, sw1s, swtt, swtt1, swu.

Removed unused arguments: paer of sw, sw1s, sw2s, swclr; pcldsw of
sw1s, sw2s; pdsig, prayl of swr; co2_ppm of clmain, clqh; tsol of
transp_lay; nsrf of screenp; kcrit and kknu of gwstress; pstd of
orosetup.

Added output of relative humidity.
1
2 ! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/advyp.F,v 1.1.1.1 2004/05/19
3 ! 12:53:06 lmdzadmin Exp $
4
5 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
30 ! Rem : Probleme aux poles il faut reecrire ce cas specifique
31 ! Attention au sens de l'indexation
32
33 ! parametres principaux du modele
34
35
36
37 ! 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
42 INTEGER lon, lat, niv
43 INTEGER i, j, jv, k, kp, l
44 INTEGER ntra
45 ! PARAMETER (ntra = 1)
46
47 REAL dty
48 REAL, INTENT (IN) :: pbarv(iip1, jjm, llm)
49
50 ! 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 sqi, sqf
104 LOGICAL limit
105
106 lon = iim ! rem : Il est possible qu'un pbl. arrive ici
107 lat = jjp1 ! a cause des dim. differentes entre les
108 niv = llm ! tab. S et VGRI
109
110 ! *** Test : diag de la qtite totale de traceur dans
111 ! l'atmosphere avant l'advection en Y
112
113 sqi = 0.
114 sqf = 0.
115
116 DO l = 1, llm
117 DO j = 1, jjp1
118 DO i = 1, iim
119 sqi = sqi + s0(i, j, l, ntra)
120 END DO
121 END DO
122 END DO
123 PRINT *, '---------- DIAG DANS ADVY - ENTREE --------'
124 PRINT *, 'sqi=', sqi
125
126 ! -----------------------------------------------------------------
127 ! Interface : adaptation nouveau modele
128 ! -------------------------------------
129
130 ! Conversion des flux de masses en kg
131 ! -AA 20/10/94 le signe -1 est necessaire car indexation opposee
132
133 DO l = 1, llm
134 DO j = 1, jjm
135 DO i = 1, iip1
136 vgri(i, j, llm+1-l) = -1.*pbarv(i, j, l)
137 END DO
138 END DO
139 END DO
140
141 ! AA Initialisation de flux fictifs aux bords sup. des boites pol.
142
143 DO l = 1, llm
144 DO i = 1, iip1
145 vgri(i, 0, l) = 0.
146 vgri(i, jjp1, l) = 0.
147 END DO
148 END DO
149
150 ! ----------------- START HERE -----------------------
151 ! boucle sur les niveaux
152
153 DO l = 1, niv
154
155 ! place limits on appropriate moments before transport
156 ! (if flux-limiting is to be applied)
157
158 IF (.NOT. limit) GO TO 11
159
160 DO jv = 1, ntra
161 DO k = 1, lat
162 DO i = 1, lon
163 IF (s0(i,k,l,jv)>0.) THEN
164 slpmax = amax1(s0(i,k,l,jv), 0.)
165 s1max = 1.5*slpmax
166 s1new = amin1(s1max, amax1(-s1max,sy(i,k,l,jv)))
167 s2new = amin1(2.*slpmax-abs(s1new)/3., amax1(abs( &
168 s1new)-slpmax,syy(i,k,l,jv)))
169 sy(i, k, l, jv) = s1new
170 syy(i, k, l, jv) = s2new
171 ssxy(i, k, l, jv) = amin1(slpmax, amax1(-slpmax,ssxy(i,k,l,jv)))
172 syz(i, k, l, jv) = amin1(slpmax, amax1(-slpmax,syz(i,k,l,jv)))
173 ELSE
174 sy(i, k, l, jv) = 0.
175 syy(i, k, l, jv) = 0.
176 ssxy(i, k, l, jv) = 0.
177 syz(i, k, l, jv) = 0.
178 END IF
179 END DO
180 END DO
181 END DO
182
183 11 CONTINUE
184
185 ! le flux a travers le pole Nord est traite separement
186
187 sm0 = 0.
188 DO jv = 1, ntra
189 s00(jv) = 0.
190 END DO
191
192 DO i = 1, lon
193
194 IF (vgri(i,0,l)<=0.) THEN
195 fm(i, 0) = -vgri(i, 0, l)*dty
196 alf(i, 0) = fm(i, 0)/sm(i, 1, l)
197 sm(i, 1, l) = sm(i, 1, l) - fm(i, 0)
198 sm0 = sm0 + fm(i, 0)
199 END IF
200
201 alfq(i, 0) = alf(i, 0)*alf(i, 0)
202 alf1(i, 0) = 1. - alf(i, 0)
203 alf1q(i, 0) = alf1(i, 0)*alf1(i, 0)
204 alf2(i, 0) = alf1(i, 0) - alf(i, 0)
205 alf3(i, 0) = alf(i, 0)*alfq(i, 0)
206 alf4(i, 0) = alf1(i, 0)*alf1q(i, 0)
207
208 END DO
209 ! print*,'ADVYP 21'
210
211 DO jv = 1, ntra
212 DO i = 1, lon
213
214 IF (vgri(i,0,l)<=0.) THEN
215
216 f0(i, 0, jv) = alf(i, 0)*(s0(i,1,l,jv)-alf1(i,0)*(sy(i,1,l, &
217 jv)-alf2(i,0)*syy(i,1,l,jv)))
218
219 s00(jv) = s00(jv) + f0(i, 0, jv)
220 s0(i, 1, l, jv) = s0(i, 1, l, jv) - f0(i, 0, jv)
221 sy(i, 1, l, jv) = alf1q(i, 0)*(sy(i,1,l,jv)+3.*alf(i,0)*syy(i,1,l, &
222 jv))
223 syy(i, 1, l, jv) = alf4(i, 0)*syy(i, 1, l, jv)
224 ssx(i, 1, l, jv) = alf1(i, 0)*(ssx(i,1,l,jv)+alf(i,0)*ssxy(i,1,l,jv &
225 ))
226 sz(i, 1, l, jv) = alf1(i, 0)*(sz(i,1,l,jv)+alf(i,0)*ssxz(i,1,l,jv))
227 ssxx(i, 1, l, jv) = alf1(i, 0)*ssxx(i, 1, l, jv)
228 ssxz(i, 1, l, jv) = alf1(i, 0)*ssxz(i, 1, l, jv)
229 szz(i, 1, l, jv) = alf1(i, 0)*szz(i, 1, l, jv)
230 ssxy(i, 1, l, jv) = alf1q(i, 0)*ssxy(i, 1, l, jv)
231 syz(i, 1, l, jv) = alf1q(i, 0)*syz(i, 1, l, jv)
232
233 END IF
234
235 END DO
236 END DO
237
238 DO i = 1, lon
239 IF (vgri(i,0,l)>0.) THEN
240 fm(i, 0) = vgri(i, 0, l)*dty
241 alf(i, 0) = fm(i, 0)/sm0
242 END IF
243 END DO
244
245 DO jv = 1, ntra
246 DO i = 1, lon
247 IF (vgri(i,0,l)>0.) THEN
248 f0(i, 0, jv) = alf(i, 0)*s00(jv)
249 END IF
250 END DO
251 END DO
252
253 ! puts the temporary moments Fi into appropriate neighboring boxes
254
255 ! print*,'av ADVYP 25'
256 DO i = 1, lon
257
258 IF (vgri(i,0,l)>0.) THEN
259 sm(i, 1, l) = sm(i, 1, l) + fm(i, 0)
260 alf(i, 0) = fm(i, 0)/sm(i, 1, l)
261 END IF
262
263 alfq(i, 0) = alf(i, 0)*alf(i, 0)
264 alf1(i, 0) = 1. - alf(i, 0)
265 alf1q(i, 0) = alf1(i, 0)*alf1(i, 0)
266 alf2(i, 0) = alf1(i, 0) - alf(i, 0)
267 alf3(i, 0) = alf1(i, 0)*alf(i, 0)
268
269 END DO
270 ! print*,'av ADVYP 25'
271
272 DO jv = 1, ntra
273 DO i = 1, lon
274
275 IF (vgri(i,0,l)>0.) THEN
276
277 temptm = alf(i, 0)*s0(i, 1, l, jv) - alf1(i, 0)*f0(i, 0, jv)
278 s0(i, 1, l, jv) = s0(i, 1, l, jv) + f0(i, 0, jv)
279 syy(i, 1, l, jv) = alf1q(i, 0)*syy(i, 1, l, jv) + &
280 5.*(alf3(i,0)*sy(i,1,l,jv)-alf2(i,0)*temptm)
281 sy(i, 1, l, jv) = alf1(i, 0)*sy(i, 1, l, jv) + 3.*temptm
282 ssxy(i, 1, l, jv) = alf1(i, 0)*ssxy(i, 1, l, jv) + &
283 3.*alf(i, 0)*ssx(i, 1, l, jv)
284 syz(i, 1, l, jv) = alf1(i, 0)*syz(i, 1, l, jv) + &
285 3.*alf(i, 0)*sz(i, 1, l, jv)
286
287 END IF
288
289 END DO
290 END DO
291
292 ! calculate flux and moments between adjacent boxes
293 ! 1- create temporary moments/masses for partial boxes in transit
294 ! 2- reajusts moments remaining in the box
295
296 ! flux from KP to K if V(K).lt.0 and from K to KP if V(K).gt.0
297
298 ! print*,'av ADVYP 30'
299 DO k = 1, lat - 1
300 kp = k + 1
301 DO i = 1, lon
302
303 IF (vgri(i,k,l)<0.) THEN
304 fm(i, k) = -vgri(i, k, l)*dty
305 alf(i, k) = fm(i, k)/sm(i, kp, l)
306 sm(i, kp, l) = sm(i, kp, l) - fm(i, k)
307 ELSE
308 fm(i, k) = vgri(i, k, l)*dty
309 alf(i, k) = fm(i, k)/sm(i, k, l)
310 sm(i, k, l) = sm(i, k, l) - fm(i, k)
311 END IF
312
313 alfq(i, k) = alf(i, k)*alf(i, k)
314 alf1(i, k) = 1. - alf(i, k)
315 alf1q(i, k) = alf1(i, k)*alf1(i, k)
316 alf2(i, k) = alf1(i, k) - alf(i, k)
317 alf3(i, k) = alf(i, k)*alfq(i, k)
318 alf4(i, k) = alf1(i, k)*alf1q(i, k)
319
320 END DO
321 END DO
322 ! print*,'ap ADVYP 30'
323
324 DO jv = 1, ntra
325 DO k = 1, lat - 1
326 kp = k + 1
327 DO i = 1, lon
328
329 IF (vgri(i,k,l)<0.) THEN
330
331 f0(i, k, jv) = alf(i, k)*(s0(i,kp,l,jv)-alf1(i,k)*(sy(i,kp,l, &
332 jv)-alf2(i,k)*syy(i,kp,l,jv)))
333 fy(i, k, jv) = alfq(i, k)*(sy(i,kp,l,jv)-3.*alf1(i,k)*syy(i,kp,l, &
334 jv))
335 fyy(i, k, jv) = alf3(i, k)*syy(i, kp, l, jv)
336 fx(i, k, jv) = alf(i, k)*(ssx(i,kp,l,jv)-alf1(i,k)*ssxy(i,kp,l,jv &
337 ))
338 fz(i, k, jv) = alf(i, k)*(sz(i,kp,l,jv)-alf1(i,k)*syz(i,kp,l,jv))
339 fxy(i, k, jv) = alfq(i, k)*ssxy(i, kp, l, jv)
340 fyz(i, k, jv) = alfq(i, k)*syz(i, kp, l, jv)
341 fxx(i, k, jv) = alf(i, k)*ssxx(i, kp, l, jv)
342 fxz(i, k, jv) = alf(i, k)*ssxz(i, kp, l, jv)
343 fzz(i, k, jv) = alf(i, k)*szz(i, kp, l, jv)
344
345 s0(i, kp, l, jv) = s0(i, kp, l, jv) - f0(i, k, jv)
346 sy(i, kp, l, jv) = alf1q(i, k)*(sy(i,kp,l,jv)+3.*alf(i,k)*syy(i, &
347 kp,l,jv))
348 syy(i, kp, l, jv) = alf4(i, k)*syy(i, kp, l, jv)
349 ssx(i, kp, l, jv) = ssx(i, kp, l, jv) - fx(i, k, jv)
350 sz(i, kp, l, jv) = sz(i, kp, l, jv) - fz(i, k, jv)
351 ssxx(i, kp, l, jv) = ssxx(i, kp, l, jv) - fxx(i, k, jv)
352 ssxz(i, kp, l, jv) = ssxz(i, kp, l, jv) - fxz(i, k, jv)
353 szz(i, kp, l, jv) = szz(i, kp, l, jv) - fzz(i, k, jv)
354 ssxy(i, kp, l, jv) = alf1q(i, k)*ssxy(i, kp, l, jv)
355 syz(i, kp, l, jv) = alf1q(i, k)*syz(i, kp, l, jv)
356
357 ELSE
358
359 f0(i, k, jv) = alf(i, k)*(s0(i,k,l,jv)+alf1(i,k)*(sy(i,k,l, &
360 jv)+alf2(i,k)*syy(i,k,l,jv)))
361 fy(i, k, jv) = alfq(i, k)*(sy(i,k,l,jv)+3.*alf1(i,k)*syy(i,k,l,jv &
362 ))
363 fyy(i, k, jv) = alf3(i, k)*syy(i, k, l, jv)
364 fx(i, k, jv) = alf(i, k)*(ssx(i,k,l,jv)+alf1(i,k)*ssxy(i,k,l,jv))
365 fz(i, k, jv) = alf(i, k)*(sz(i,k,l,jv)+alf1(i,k)*syz(i,k,l,jv))
366 fxy(i, k, jv) = alfq(i, k)*ssxy(i, k, l, jv)
367 fyz(i, k, jv) = alfq(i, k)*syz(i, k, l, jv)
368 fxx(i, k, jv) = alf(i, k)*ssxx(i, k, l, jv)
369 fxz(i, k, jv) = alf(i, k)*ssxz(i, k, l, jv)
370 fzz(i, k, jv) = alf(i, k)*szz(i, k, l, jv)
371
372 s0(i, k, l, jv) = s0(i, k, l, jv) - f0(i, k, jv)
373 sy(i, k, l, jv) = alf1q(i, k)*(sy(i,k,l,jv)-3.*alf(i,k)*syy(i,k,l &
374 ,jv))
375 syy(i, k, l, jv) = alf4(i, k)*syy(i, k, l, jv)
376 ssx(i, k, l, jv) = ssx(i, k, l, jv) - fx(i, k, jv)
377 sz(i, k, l, jv) = sz(i, k, l, jv) - fz(i, k, jv)
378 ssxx(i, k, l, jv) = ssxx(i, k, l, jv) - fxx(i, k, jv)
379 ssxz(i, k, l, jv) = ssxz(i, k, l, jv) - fxz(i, k, jv)
380 szz(i, k, l, jv) = szz(i, k, l, jv) - fzz(i, k, jv)
381 ssxy(i, k, l, jv) = alf1q(i, k)*ssxy(i, k, l, jv)
382 syz(i, k, l, jv) = alf1q(i, k)*syz(i, k, l, jv)
383
384 END IF
385
386 END DO
387 END DO
388 END DO
389 ! print*,'ap ADVYP 31'
390
391 ! puts the temporary moments Fi into appropriate neighboring boxes
392
393 DO k = 1, lat - 1
394 kp = k + 1
395 DO i = 1, lon
396
397 IF (vgri(i,k,l)<0.) THEN
398 sm(i, k, l) = sm(i, k, l) + fm(i, k)
399 alf(i, k) = fm(i, k)/sm(i, k, l)
400 ELSE
401 sm(i, kp, l) = sm(i, kp, l) + fm(i, k)
402 alf(i, k) = fm(i, k)/sm(i, kp, l)
403 END IF
404
405 alfq(i, k) = alf(i, k)*alf(i, k)
406 alf1(i, k) = 1. - alf(i, k)
407 alf1q(i, k) = alf1(i, k)*alf1(i, k)
408 alf2(i, k) = alf1(i, k) - alf(i, k)
409 alf3(i, k) = alf1(i, k)*alf(i, k)
410
411 END DO
412 END DO
413 ! print*,'ap ADVYP 32'
414
415 DO jv = 1, ntra
416 DO k = 1, lat - 1
417 kp = k + 1
418 DO i = 1, lon
419
420 IF (vgri(i,k,l)<0.) THEN
421
422 temptm = -alf(i, k)*s0(i, k, l, jv) + alf1(i, k)*f0(i, k, jv)
423 s0(i, k, l, jv) = s0(i, k, l, jv) + f0(i, k, jv)
424 syy(i, k, l, jv) = alfq(i, k)*fyy(i, k, jv) + &
425 alf1q(i, k)*syy(i, k, l, jv) + 5.*(alf3(i,k)*(fy(i,k,jv)-sy(i, &
426 k,l,jv))+alf2(i,k)*temptm)
427 sy(i, k, l, jv) = alf(i, k)*fy(i, k, jv) + &
428 alf1(i, k)*sy(i, k, l, jv) + 3.*temptm
429 ssxy(i, k, l, jv) = alf(i, k)*fxy(i, k, jv) + &
430 alf1(i, k)*ssxy(i, k, l, jv) + 3.*(alf1(i,k)*fx(i,k,jv)-alf(i,k &
431 )*ssx(i,k,l,jv))
432 syz(i, k, l, jv) = alf(i, k)*fyz(i, k, jv) + &
433 alf1(i, k)*syz(i, k, l, jv) + 3.*(alf1(i,k)*fz(i,k,jv)-alf(i,k) &
434 *sz(i,k,l,jv))
435 ssx(i, k, l, jv) = ssx(i, k, l, jv) + fx(i, k, jv)
436 sz(i, k, l, jv) = sz(i, k, l, jv) + fz(i, k, jv)
437 ssxx(i, k, l, jv) = ssxx(i, k, l, jv) + fxx(i, k, jv)
438 ssxz(i, k, l, jv) = ssxz(i, k, l, jv) + fxz(i, k, jv)
439 szz(i, k, l, jv) = szz(i, k, l, jv) + fzz(i, k, jv)
440
441 ELSE
442
443 temptm = alf(i, k)*s0(i, kp, l, jv) - alf1(i, k)*f0(i, k, jv)
444 s0(i, kp, l, jv) = s0(i, kp, l, jv) + f0(i, k, jv)
445 syy(i, kp, l, jv) = alfq(i, k)*fyy(i, k, jv) + &
446 alf1q(i, k)*syy(i, kp, l, jv) + 5.*(alf3(i,k)*(sy(i,kp,l, &
447 jv)-fy(i,k,jv))-alf2(i,k)*temptm)
448 sy(i, kp, l, jv) = alf(i, k)*fy(i, k, jv) + &
449 alf1(i, k)*sy(i, kp, l, jv) + 3.*temptm
450 ssxy(i, kp, l, jv) = alf(i, k)*fxy(i, k, jv) + &
451 alf1(i, k)*ssxy(i, kp, l, jv) + 3.*(alf(i,k)*ssx(i,kp,l,jv)- &
452 alf1(i,k)*fx(i,k,jv))
453 syz(i, kp, l, jv) = alf(i, k)*fyz(i, k, jv) + &
454 alf1(i, k)*syz(i, kp, l, jv) + 3.*(alf(i,k)*sz(i,kp,l,jv)-alf1( &
455 i,k)*fz(i,k,jv))
456 ssx(i, kp, l, jv) = ssx(i, kp, l, jv) + fx(i, k, jv)
457 sz(i, kp, l, jv) = sz(i, kp, l, jv) + fz(i, k, jv)
458 ssxx(i, kp, l, jv) = ssxx(i, kp, l, jv) + fxx(i, k, jv)
459 ssxz(i, kp, l, jv) = ssxz(i, kp, l, jv) + fxz(i, k, jv)
460 szz(i, kp, l, jv) = szz(i, kp, l, jv) + fzz(i, k, jv)
461
462 END IF
463
464 END DO
465 END DO
466 END DO
467 ! print*,'ap ADVYP 33'
468
469 ! traitement special pour le pole Sud (idem pole Nord)
470
471 k = lat
472
473 sm0 = 0.
474 DO jv = 1, ntra
475 s00(jv) = 0.
476 END DO
477
478 DO i = 1, lon
479
480 IF (vgri(i,k,l)>=0.) THEN
481 fm(i, k) = vgri(i, k, l)*dty
482 alf(i, k) = fm(i, k)/sm(i, k, l)
483 sm(i, k, l) = sm(i, k, l) - fm(i, k)
484 sm0 = sm0 + fm(i, k)
485 END IF
486
487 alfq(i, k) = alf(i, k)*alf(i, k)
488 alf1(i, k) = 1. - alf(i, k)
489 alf1q(i, k) = alf1(i, k)*alf1(i, k)
490 alf2(i, k) = alf1(i, k) - alf(i, k)
491 alf3(i, k) = alf(i, k)*alfq(i, k)
492 alf4(i, k) = alf1(i, k)*alf1q(i, k)
493
494 END DO
495 ! print*,'ap ADVYP 41'
496
497 DO jv = 1, ntra
498 DO i = 1, lon
499
500 IF (vgri(i,k,l)>=0.) THEN
501 f0(i, k, jv) = alf(i, k)*(s0(i,k,l,jv)+alf1(i,k)*(sy(i,k,l, &
502 jv)+alf2(i,k)*syy(i,k,l,jv)))
503 s00(jv) = s00(jv) + f0(i, k, jv)
504
505 s0(i, k, l, jv) = s0(i, k, l, jv) - f0(i, k, jv)
506 sy(i, k, l, jv) = alf1q(i, k)*(sy(i,k,l,jv)-3.*alf(i,k)*syy(i,k,l, &
507 jv))
508 syy(i, k, l, jv) = alf4(i, k)*syy(i, k, l, jv)
509 ssx(i, k, l, jv) = alf1(i, k)*(ssx(i,k,l,jv)-alf(i,k)*ssxy(i,k,l,jv &
510 ))
511 sz(i, k, l, jv) = alf1(i, k)*(sz(i,k,l,jv)-alf(i,k)*syz(i,k,l,jv))
512 ssxx(i, k, l, jv) = alf1(i, k)*ssxx(i, k, l, jv)
513 ssxz(i, k, l, jv) = alf1(i, k)*ssxz(i, k, l, jv)
514 szz(i, k, l, jv) = alf1(i, k)*szz(i, k, l, jv)
515 ssxy(i, k, l, jv) = alf1q(i, k)*ssxy(i, k, l, jv)
516 syz(i, k, l, jv) = alf1q(i, k)*syz(i, k, l, jv)
517 END IF
518
519 END DO
520 END DO
521 ! print*,'ap ADVYP 42'
522
523 DO i = 1, lon
524 IF (vgri(i,k,l)<0.) THEN
525 fm(i, k) = -vgri(i, k, l)*dty
526 alf(i, k) = fm(i, k)/sm0
527 END IF
528 END DO
529 ! print*,'ap ADVYP 43'
530
531 DO jv = 1, ntra
532 DO i = 1, lon
533 IF (vgri(i,k,l)<0.) THEN
534 f0(i, k, jv) = alf(i, k)*s00(jv)
535 END IF
536 END DO
537 END DO
538
539 ! puts the temporary moments Fi into appropriate neighboring boxes
540
541 DO i = 1, lon
542
543 IF (vgri(i,k,l)<0.) THEN
544 sm(i, k, l) = sm(i, k, l) + fm(i, k)
545 alf(i, k) = fm(i, k)/sm(i, k, l)
546 END IF
547
548 alfq(i, k) = alf(i, k)*alf(i, k)
549 alf1(i, k) = 1. - alf(i, k)
550 alf1q(i, k) = alf1(i, k)*alf1(i, k)
551 alf2(i, k) = alf1(i, k) - alf(i, k)
552 alf3(i, k) = alf1(i, k)*alf(i, k)
553
554 END DO
555 ! print*,'ap ADVYP 45'
556
557 DO jv = 1, ntra
558 DO i = 1, lon
559
560 IF (vgri(i,k,l)<0.) THEN
561
562 temptm = -alf(i, k)*s0(i, k, l, jv) + alf1(i, k)*f0(i, k, jv)
563 s0(i, k, l, jv) = s0(i, k, l, jv) + f0(i, k, jv)
564 syy(i, k, l, jv) = alf1q(i, k)*syy(i, k, l, jv) + &
565 5.*(-alf3(i,k)*sy(i,k,l,jv)+alf2(i,k)*temptm)
566 sy(i, k, l, jv) = alf1(i, k)*sy(i, k, l, jv) + 3.*temptm
567 ssxy(i, k, l, jv) = alf1(i, k)*ssxy(i, k, l, jv) - &
568 3.*alf(i, k)*ssx(i, k, l, jv)
569 syz(i, k, l, jv) = alf1(i, k)*syz(i, k, l, jv) - &
570 3.*alf(i, k)*sz(i, k, l, jv)
571
572 END IF
573
574 END DO
575 END DO
576 ! print*,'ap ADVYP 46'
577
578 END DO
579
580 ! --------------------------------------------------
581 ! bouclage cyclique horizontal .
582
583 DO l = 1, llm
584 DO jv = 1, ntra
585 DO j = 1, jjp1
586 sm(iip1, j, l) = sm(1, j, l)
587 s0(iip1, j, l, jv) = s0(1, j, l, jv)
588 ssx(iip1, j, l, jv) = ssx(1, j, l, jv)
589 sy(iip1, j, l, jv) = sy(1, j, l, jv)
590 sz(iip1, j, l, jv) = sz(1, j, l, jv)
591 END DO
592 END DO
593 END DO
594
595 ! -------------------------------------------------------------------
596 ! *** Test negativite:
597
598 ! DO jv = 1,ntra
599 ! DO l = 1,llm
600 ! DO j = 1,jjp1
601 ! DO i = 1,iip1
602 ! IF (s0( i,j,l,jv ).lt.0.) THEN
603 ! PRINT*, '------ S0 < 0 en FIN ADVYP ---'
604 ! PRINT*, 'S0(',i,j,l,jv,')=', S0(i,j,l,jv)
605 ! c STOP
606 ! ENDIF
607 ! ENDDO
608 ! ENDDO
609 ! ENDDO
610 ! ENDDO
611
612
613 ! -------------------------------------------------------------------
614 ! *** Test : diag de la qtite totale de traceur dans
615 ! l'atmosphere avant l'advection en Y
616
617 DO l = 1, llm
618 DO j = 1, jjp1
619 DO i = 1, iim
620 sqf = sqf + s0(i, j, l, ntra)
621 END DO
622 END DO
623 END DO
624 PRINT *, '---------- DIAG DANS ADVY - SORTIE --------'
625 PRINT *, 'sqf=', sqf
626 ! print*,'ap ADVYP fin'
627
628 ! -----------------------------------------------------------------
629
630 RETURN
631 END SUBROUTINE advyp
632
633
634
635
636
637
638
639
640
641
642
643
  ViewVC Help
Powered by ViewVC 1.1.21