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

Contents of /trunk/dyn3d/advyp.f

Parent Directory Parent Directory | Revision Log Revision Log

Revision 113 - (show annotations)
Thu Sep 18 19:56:46 2014 UTC (9 years, 7 months ago) by guez
File size: 19750 byte(s) 
Moved the call to read_serre out of conf_gcm so that it can be called
only in the program ce0l, not in gcm. In gcm, variables of module
serre are read from start file. Added reading of dzoomx, dzoomy, taux,
tauy from start file, in dynetat0. Those variables were written by
dynredem0 but not read.

Removed possibility fxyhypb = false, because the geometric part of the
program is such a mess. Could then remove variables transx, transy,
alphax, alphay, pxo, pyo of module serre.

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