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

Contents of /trunk/dyn3d/advn.f

Parent Directory Parent Directory | Revision Log Revision Log

Revision 112 - (show annotations)
Thu Sep 18 13:36:51 2014 UTC (9 years, 7 months ago) by guez
File size: 22469 byte(s) 
Removed 8 first arguments of fxyhyper, use variables of module serre
instead.

Moved reading of variables of module serre from procedure conf_gcm to
new procedure read_serre.

In guide, added conditions to avoid useless calls to tau2alpha and
writefield. Bugfix: offline corresponds to alpha = 1. Open only one
NetCDF file to read number of vertical levels.

In tau2alpha, added conditions to avoid useless computations of dxdyu
and dxdyv. gamma is not needed for a regular grid.
1
2 ! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/advn.F,v 1.1.1.1 2004/05/19
3 ! 12:53:06 lmdzadmin Exp $
4
5 SUBROUTINE advn(q, masse, w, pbaru, pbarv, pdt, mode)
6
7 ! Auteur : F. Hourdin
8
9 ! ********************************************************************
10 ! Shema d'advection " pseudo amont " .
11 ! ********************************************************************
12 ! q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg ....
13
14 ! pbaru,pbarv,w flux de masse en u ,v ,w
15 ! pdt pas de temps
16
17 ! --------------------------------------------------------------------
18 USE dimens_m
19 USE paramet_m
20 USE comconst
21 USE disvert_m
22 USE conf_gcm_m
23 USE comgeom
24 IMPLICIT NONE
25
26
27
28 ! Arguments:
29 ! ----------
30 INTEGER mode
31 REAL masse(ip1jmp1, llm)
32 REAL, INTENT (IN) :: pbaru(ip1jmp1, llm), pbarv(ip1jm, llm)
33 REAL q(ip1jmp1, llm)
34 REAL w(ip1jmp1, llm), pdt
35
36 ! Local
37 ! ---------
38
39 INTEGER ij, l
40 REAL zm(ip1jmp1, llm)
41 REAL mu(ip1jmp1, llm)
42 REAL mv(ip1jm, llm)
43 REAL mw(ip1jmp1, llm+1)
44 REAL zq(ip1jmp1, llm), qpn, qps
45 REAL zqg(ip1jmp1, llm), zqd(ip1jmp1, llm)
46 REAL zqs(ip1jmp1, llm), zqn(ip1jmp1, llm)
47 REAL zqh(ip1jmp1, llm), zqb(ip1jmp1, llm)
48 REAL temps1, temps2, temps3
49 REAL ssum
50 LOGICAL testcpu
51 SAVE testcpu
52 SAVE temps1, temps2, temps3
53 REAL zzpbar, zzw
54
55 REAL qmin, qmax
56 DATA qmin, qmax/0., 1./
57 DATA testcpu/.FALSE./
58 DATA temps1, temps2, temps3/0., 0., 0./
59
60 zzpbar = 0.5*pdt
61 zzw = pdt
62
63 DO l = 1, llm
64 DO ij = iip2, ip1jm
65 mu(ij, l) = pbaru(ij, l)*zzpbar
66 END DO
67 DO ij = 1, ip1jm
68 mv(ij, l) = pbarv(ij, l)*zzpbar
69 END DO
70 DO ij = 1, ip1jmp1
71 mw(ij, l) = w(ij, l)*zzw
72 END DO
73 END DO
74
75 DO ij = 1, ip1jmp1
76 mw(ij, llm+1) = 0.
77 END DO
78
79 DO l = 1, llm
80 qpn = 0.
81 qps = 0.
82 DO ij = 1, iim
83 qpn = qpn + q(ij, l)*masse(ij, l)
84 qps = qps + q(ip1jm+ij, l)*masse(ip1jm+ij, l)
85 END DO
86 qpn = qpn/ssum(iim, masse(1,l), 1)
87 qps = qps/ssum(iim, masse(ip1jm+1,l), 1)
88 DO ij = 1, iip1
89 q(ij, l) = qpn
90 q(ip1jm+ij, l) = qps
91 END DO
92 END DO
93
94 DO ij = 1, ip1jmp1
95 mw(ij, llm+1) = 0.
96 END DO
97 DO l = 1, llm
98 DO ij = 1, ip1jmp1
99 zq(ij, l) = q(ij, l)
100 zm(ij, l) = masse(ij, l)
101 END DO
102 END DO
103
104 ! call minmaxq(zq,qmin,qmax,'avant vlx ')
105 CALL advnqx(zq, zqg, zqd)
106 CALL advnx(zq, zqg, zqd, zm, mu, mode)
107 CALL advnqy(zq, zqs, zqn)
108 CALL advny(zq, zqs, zqn, zm, mv)
109 CALL advnqz(zq, zqh, zqb)
110 CALL advnz(zq, zqh, zqb, zm, mw)
111 ! call vlz(zq,0.,zm,mw)
112 CALL advnqy(zq, zqs, zqn)
113 CALL advny(zq, zqs, zqn, zm, mv)
114 CALL advnqx(zq, zqg, zqd)
115 CALL advnx(zq, zqg, zqd, zm, mu, mode)
116 ! call minmaxq(zq,qmin,qmax,'apres vlx ')
117
118 DO l = 1, llm
119 DO ij = 1, ip1jmp1
120 q(ij, l) = zq(ij, l)
121 END DO
122 DO ij = 1, ip1jm + 1, iip1
123 q(ij+iim, l) = q(ij, l)
124 END DO
125 END DO
126
127 RETURN
128 END SUBROUTINE advn
129
130 SUBROUTINE advnqx(q, qg, qd)
131
132 ! Auteurs: Calcul des valeurs de q aux point u.
133
134 ! --------------------------------------------------------------------
135 USE dimens_m
136 USE paramet_m
137 USE conf_gcm_m
138 IMPLICIT NONE
139
140
141
142 ! Arguments:
143 ! ----------
144 REAL q(ip1jmp1, llm), qg(ip1jmp1, llm), qd(ip1jmp1, llm)
145
146 ! Local
147 ! ---------
148
149 INTEGER ij, l
150
151 REAL dxqu(ip1jmp1), zqu(ip1jmp1)
152 REAL zqmax(ip1jmp1), zqmin(ip1jmp1)
153 LOGICAL extremum(ip1jmp1)
154
155 INTEGER mode
156 SAVE mode
157 DATA mode/1/
158
159 ! calcul des pentes en u:
160 ! -----------------------
161 IF (mode==0) THEN
162 DO l = 1, llm
163 DO ij = 1, ip1jm
164 qd(ij, l) = q(ij, l)
165 qg(ij, l) = q(ij, l)
166 END DO
167 END DO
168 ELSE
169 DO l = 1, llm
170 DO ij = iip2, ip1jm - 1
171 dxqu(ij) = q(ij+1, l) - q(ij, l)
172 zqu(ij) = 0.5*(q(ij+1,l)+q(ij,l))
173 END DO
174 DO ij = iip1 + iip1, ip1jm, iip1
175 dxqu(ij) = dxqu(ij-iim)
176 zqu(ij) = zqu(ij-iim)
177 END DO
178 DO ij = iip2, ip1jm - 1
179 zqu(ij) = zqu(ij) - dxqu(ij+1)/12.
180 END DO
181 DO ij = iip1 + iip1, ip1jm, iip1
182 zqu(ij) = zqu(ij-iim)
183 END DO
184 DO ij = iip2 + 1, ip1jm
185 zqu(ij) = zqu(ij) + dxqu(ij-1)/12.
186 END DO
187 DO ij = iip1 + iip1, ip1jm, iip1
188 zqu(ij-iim) = zqu(ij)
189 END DO
190
191 ! calcul des valeurs max et min acceptees aux interfaces
192
193 DO ij = iip2, ip1jm - 1
194 zqmax(ij) = max(q(ij+1,l), q(ij,l))
195 zqmin(ij) = min(q(ij+1,l), q(ij,l))
196 END DO
197 DO ij = iip1 + iip1, ip1jm, iip1
198 zqmax(ij) = zqmax(ij-iim)
199 zqmin(ij) = zqmin(ij-iim)
200 END DO
201 DO ij = iip2 + 1, ip1jm
202 extremum(ij) = dxqu(ij)*dxqu(ij-1) <= 0.
203 END DO
204 DO ij = iip1 + iip1, ip1jm, iip1
205 extremum(ij-iim) = extremum(ij)
206 END DO
207 DO ij = iip2, ip1jm
208 zqu(ij) = min(max(zqmin(ij),zqu(ij)), zqmax(ij))
209 END DO
210 DO ij = iip2 + 1, ip1jm
211 IF (extremum(ij)) THEN
212 qg(ij, l) = q(ij, l)
213 qd(ij, l) = q(ij, l)
214 ELSE
215 qd(ij, l) = zqu(ij)
216 qg(ij, l) = zqu(ij-1)
217 END IF
218 END DO
219 DO ij = iip1 + iip1, ip1jm, iip1
220 qd(ij-iim, l) = qd(ij, l)
221 qg(ij-iim, l) = qg(ij, l)
222 END DO
223
224 GO TO 8888
225
226 DO ij = iip2 + 1, ip1jm
227 IF (extremum(ij) .AND. .NOT. extremum(ij-1)) qd(ij-1, l) = q(ij, l)
228 END DO
229
230 DO ij = iip1 + iip1, ip1jm, iip1
231 qd(ij-iim, l) = qd(ij, l)
232 END DO
233 DO ij = iip2, ip1jm - 1
234 IF (extremum(ij) .AND. .NOT. extremum(ij+1)) qg(ij+1, l) = q(ij, l)
235 END DO
236
237 DO ij = iip1 + iip1, ip1jm, iip1
238 qg(ij, l) = qg(ij-iim, l)
239 END DO
240 8888 CONTINUE
241 END DO
242 END IF
243 RETURN
244 END SUBROUTINE advnqx
245 SUBROUTINE advnqy(q, qs, qn)
246
247 ! Auteurs: Calcul des valeurs de q aux point v.
248
249 ! --------------------------------------------------------------------
250 USE dimens_m
251 USE paramet_m
252 USE conf_gcm_m
253 IMPLICIT NONE
254
255
256
257 ! Arguments:
258 ! ----------
259 REAL q(ip1jmp1, llm), qs(ip1jmp1, llm), qn(ip1jmp1, llm)
260
261 ! Local
262 ! ---------
263
264 INTEGER ij, l
265
266 REAL dyqv(ip1jm), zqv(ip1jm, llm)
267 REAL zqmax(ip1jm), zqmin(ip1jm)
268 LOGICAL extremum(ip1jmp1)
269
270 INTEGER mode
271 SAVE mode
272 DATA mode/1/
273
274 IF (mode==0) THEN
275 DO l = 1, llm
276 DO ij = 1, ip1jmp1
277 qn(ij, l) = q(ij, l)
278 qs(ij, l) = q(ij, l)
279 END DO
280 END DO
281 ELSE
282
283 ! calcul des pentes en u:
284 ! -----------------------
285 DO l = 1, llm
286 DO ij = 1, ip1jm
287 dyqv(ij) = q(ij, l) - q(ij+iip1, l)
288 END DO
289
290 DO ij = iip2, ip1jm - iip1
291 zqv(ij, l) = 0.5*(q(ij+iip1,l)+q(ij,l))
292 zqv(ij, l) = zqv(ij, l) + (dyqv(ij+iip1)-dyqv(ij-iip1))/12.
293 END DO
294
295 DO ij = iip2, ip1jm
296 extremum(ij) = dyqv(ij)*dyqv(ij-iip1) <= 0.
297 END DO
298
299 ! Pas de pentes aux poles
300 DO ij = 1, iip1
301 zqv(ij, l) = q(ij, l)
302 zqv(ip1jm-iip1+ij, l) = q(ip1jm+ij, l)
303 extremum(ij) = .TRUE.
304 extremum(ip1jmp1-iip1+ij) = .TRUE.
305 END DO
306
307 ! calcul des valeurs max et min acceptees aux interfaces
308 DO ij = 1, ip1jm
309 zqmax(ij) = max(q(ij+iip1,l), q(ij,l))
310 zqmin(ij) = min(q(ij+iip1,l), q(ij,l))
311 END DO
312
313 DO ij = 1, ip1jm
314 zqv(ij, l) = min(max(zqmin(ij),zqv(ij,l)), zqmax(ij))
315 END DO
316
317 DO ij = iip2, ip1jm
318 IF (extremum(ij)) THEN
319 qs(ij, l) = q(ij, l)
320 qn(ij, l) = q(ij, l)
321 ! if (.not.extremum(ij-iip1)) qs(ij-iip1,l)=q(ij,l)
322 ! if (.not.extremum(ij+iip1)) qn(ij+iip1,l)=q(ij,l)
323 ELSE
324 qs(ij, l) = zqv(ij, l)
325 qn(ij, l) = zqv(ij-iip1, l)
326 END IF
327 END DO
328
329 DO ij = 1, iip1
330 qs(ij, l) = q(ij, l)
331 qn(ij, l) = q(ij, l)
332 qs(ip1jm+ij, l) = q(ip1jm+ij, l)
333 qn(ip1jm+ij, l) = q(ip1jm+ij, l)
334 END DO
335
336 END DO
337 END IF
338 RETURN
339 END SUBROUTINE advnqy
340
341 SUBROUTINE advnqz(q, qh, qb)
342
343 ! Auteurs: Calcul des valeurs de q aux point v.
344
345 ! --------------------------------------------------------------------
346 USE dimens_m
347 USE paramet_m
348 USE conf_gcm_m
349 IMPLICIT NONE
350
351
352
353 ! Arguments:
354 ! ----------
355 REAL q(ip1jmp1, llm), qh(ip1jmp1, llm), qb(ip1jmp1, llm)
356
357 ! Local
358 ! ---------
359
360 INTEGER ij, l
361
362 REAL dzqw(ip1jmp1, llm+1), zqw(ip1jmp1, llm+1)
363 REAL zqmax(ip1jmp1, llm), zqmin(ip1jmp1, llm)
364 LOGICAL extremum(ip1jmp1, llm)
365
366 INTEGER mode
367 SAVE mode
368
369 DATA mode/1/
370
371 ! calcul des pentes en u:
372 ! -----------------------
373
374 IF (mode==0) THEN
375 DO l = 1, llm
376 DO ij = 1, ip1jmp1
377 qb(ij, l) = q(ij, l)
378 qh(ij, l) = q(ij, l)
379 END DO
380 END DO
381 ELSE
382 DO l = 2, llm
383 DO ij = 1, ip1jmp1
384 dzqw(ij, l) = q(ij, l-1) - q(ij, l)
385 zqw(ij, l) = 0.5*(q(ij,l-1)+q(ij,l))
386 END DO
387 END DO
388 DO ij = 1, ip1jmp1
389 dzqw(ij, 1) = 0.
390 dzqw(ij, llm+1) = 0.
391 END DO
392 DO l = 2, llm
393 DO ij = 1, ip1jmp1
394 zqw(ij, l) = zqw(ij, l) + (dzqw(ij,l+1)-dzqw(ij,l-1))/12.
395 END DO
396 END DO
397 DO l = 2, llm - 1
398 DO ij = 1, ip1jmp1
399 extremum(ij, l) = dzqw(ij, l)*dzqw(ij, l+1) <= 0.
400 END DO
401 END DO
402
403 ! Pas de pentes en bas et en haut
404 DO ij = 1, ip1jmp1
405 zqw(ij, 2) = q(ij, 1)
406 zqw(ij, llm) = q(ij, llm)
407 extremum(ij, 1) = .TRUE.
408 extremum(ij, llm) = .TRUE.
409 END DO
410
411 ! calcul des valeurs max et min acceptees aux interfaces
412 DO l = 2, llm
413 DO ij = 1, ip1jmp1
414 zqmax(ij, l) = max(q(ij,l-1), q(ij,l))
415 zqmin(ij, l) = min(q(ij,l-1), q(ij,l))
416 END DO
417 END DO
418
419 DO l = 2, llm
420 DO ij = 1, ip1jmp1
421 zqw(ij, l) = min(max(zqmin(ij,l),zqw(ij,l)), zqmax(ij,l))
422 END DO
423 END DO
424
425 DO l = 2, llm - 1
426 DO ij = 1, ip1jmp1
427 IF (extremum(ij,l)) THEN
428 qh(ij, l) = q(ij, l)
429 qb(ij, l) = q(ij, l)
430 ELSE
431 qh(ij, l) = zqw(ij, l+1)
432 qb(ij, l) = zqw(ij, l)
433 END IF
434 END DO
435 END DO
436 ! do l=2,llm-1
437 ! do ij=1,ip1jmp1
438 ! if(extremum(ij,l)) then
439 ! if (.not.extremum(ij,l-1)) qh(ij,l-1)=q(ij,l)
440 ! if (.not.extremum(ij,l+1)) qb(ij,l+1)=q(ij,l)
441 ! endif
442 ! enddo
443 ! enddo
444
445 DO ij = 1, ip1jmp1
446 qb(ij, 1) = q(ij, 1)
447 qh(ij, 1) = q(ij, 1)
448 qb(ij, llm) = q(ij, llm)
449 qh(ij, llm) = q(ij, llm)
450 END DO
451
452 END IF
453
454 RETURN
455 END SUBROUTINE advnqz
456
457 SUBROUTINE advnx(q, qg, qd, masse, u_m, mode)
458
459 ! Auteur : F. Hourdin
460
461 ! ********************************************************************
462 ! Shema d'advection " pseudo amont " .
463 ! ********************************************************************
464 ! nq,iq,q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg ....
465
466
467 ! --------------------------------------------------------------------
468 USE dimens_m
469 USE paramet_m
470 USE comconst
471 USE disvert_m
472 USE conf_gcm_m
473 IMPLICIT NONE
474
475
476
477 ! Arguments:
478 ! ----------
479 INTEGER mode
480 REAL masse(ip1jmp1, llm)
481 REAL u_m(ip1jmp1, llm)
482 REAL q(ip1jmp1, llm), qd(ip1jmp1, llm), qg(ip1jmp1, llm)
483
484 ! Local
485 ! ---------
486
487 INTEGER i, j, ij, l, indu(ip1jmp1), niju, iju, ijq
488 INTEGER n0, nl(llm)
489
490 REAL new_m, zu_m, zdq, zz
491 REAL zsigg(ip1jmp1, llm), zsigd(ip1jmp1, llm), zsig
492 REAL u_mq(ip1jmp1, llm)
493
494 REAL zm, zq, zsigm, zsigp, zqm, zqp, zu
495
496 LOGICAL ladvplus(ip1jmp1, llm)
497
498 REAL prec
499 SAVE prec
500
501 DATA prec/1.E-15/
502
503 DO l = 1, llm
504 DO ij = iip2, ip1jm
505 zdq = qd(ij, l) - qg(ij, l)
506 IF (abs(zdq)>prec) THEN
507 zsigd(ij, l) = (q(ij,l)-qg(ij,l))/zdq
508 zsigg(ij, l) = 1. - zsigd(ij, l)
509 ELSE
510 zsigd(ij, l) = 0.5
511 zsigg(ij, l) = 0.5
512 qd(ij, l) = q(ij, l)
513 qg(ij, l) = q(ij, l)
514 END IF
515 END DO
516 END DO
517
518 ! calcul de la pente maximum dans la maille en valeur absolue
519
520 DO l = 1, llm
521 DO ij = iip2, ip1jm - 1
522 IF (u_m(ij,l)>=0.) THEN
523 zsigp = zsigd(ij, l)
524 zsigm = zsigg(ij, l)
525 zqp = qd(ij, l)
526 zqm = qg(ij, l)
527 zm = masse(ij, l)
528 zq = q(ij, l)
529 ELSE
530 zsigm = zsigd(ij+1, l)
531 zsigp = zsigg(ij+1, l)
532 zqm = qd(ij+1, l)
533 zqp = qg(ij+1, l)
534 zm = masse(ij+1, l)
535 zq = q(ij+1, l)
536 END IF
537 zu = abs(u_m(ij,l))
538 ladvplus(ij, l) = zu > zm
539 zsig = zu/zm
540 IF (zsig==0.) zsigp = 0.1
541 IF (mode==1) THEN
542 IF (zsig<=zsigp) THEN
543 u_mq(ij, l) = u_m(ij, l)*zqp
544 ELSE IF (mode==1) THEN
545 u_mq(ij, l) = sign(zm, u_m(ij,l))*(zsigp*zqp+(zsig-zsigp)*zqm)
546 END IF
547 ELSE
548 IF (zsig<=zsigp) THEN
549 u_mq(ij, l) = u_m(ij, l)*(zqp-0.5*zsig/zsigp*(zqp-zq))
550 ELSE
551 zz = 0.5*(zsig-zsigp)/zsigm
552 u_mq(ij, l) = sign(zm, u_m(ij,l))*(0.5*(zq+zqp)*zsigp+(zsig-zsigp)* &
553 (zq+zz*(zqm-zq)))
554 END IF
555 END IF
556 END DO
557 END DO
558
559 DO l = 1, llm
560 DO ij = iip1 + iip1, ip1jm, iip1
561 u_mq(ij, l) = u_mq(ij-iim, l)
562 ladvplus(ij, l) = ladvplus(ij-iim, l)
563 END DO
564 END DO
565
566 ! =================================================================
567 ! SCHEMA SEMI-LAGRAGIEN EN X DANS LES REGIONS POLAIRES
568 ! =================================================================
569 ! tris des regions a traiter
570 n0 = 0
571 DO l = 1, llm
572 nl(l) = 0
573 DO ij = iip2, ip1jm
574 IF (ladvplus(ij,l)) THEN
575 nl(l) = nl(l) + 1
576 u_mq(ij, l) = 0.
577 END IF
578 END DO
579 n0 = n0 + nl(l)
580 END DO
581
582 IF (n0>1) THEN
583 IF (prt_level>9) PRINT *, &
584 'Nombre de points pour lesquels on advect plus que le', &
585 'contenu de la maille : ', n0
586
587 DO l = 1, llm
588 IF (nl(l)>0) THEN
589 iju = 0
590 ! indicage des mailles concernees par le traitement special
591 DO ij = iip2, ip1jm
592 IF (ladvplus(ij,l) .AND. mod(ij,iip1)/=0) THEN
593 iju = iju + 1
594 indu(iju) = ij
595 END IF
596 END DO
597 niju = iju
598
599 ! traitement des mailles
600 DO iju = 1, niju
601 ij = indu(iju)
602 j = (ij-1)/iip1 + 1
603 zu_m = u_m(ij, l)
604 u_mq(ij, l) = 0.
605 IF (zu_m>0.) THEN
606 ijq = ij
607 i = ijq - (j-1)*iip1
608 ! accumulation pour les mailles completements advectees
609 DO WHILE (zu_m>masse(ijq,l))
610 u_mq(ij, l) = u_mq(ij, l) + q(ijq, l)*masse(ijq, l)
611 zu_m = zu_m - masse(ijq, l)
612 i = mod(i-2+iim, iim) + 1
613 ijq = (j-1)*iip1 + i
614 END DO
615 ! MODIFS SPECIFIQUES DU SCHEMA
616 ! ajout de la maille non completement advectee
617 zsig = zu_m/masse(ijq, l)
618 IF (zsig<=zsigd(ijq,l)) THEN
619 u_mq(ij, l) = u_mq(ij, l) + zu_m*(qd(ijq,l)-0.5*zsig/zsigd(ijq, &
620 l)*(qd(ijq,l)-q(ijq,l)))
621 ELSE
622 ! u_mq(ij,l)=u_mq(ij,l)+zu_m*q(ijq,l)
623 ! goto 8888
624 zz = 0.5*(zsig-zsigd(ijq,l))/zsigg(ijq, l)
625 IF (.NOT. (zz>0. .AND. zz<=0.5)) THEN
626 PRINT *, 'probleme2 au point ij=', ij, ' l=', l
627 PRINT *, 'zz=', zz
628 STOP
629 END IF
630 u_mq(ij, l) = u_mq(ij, l) + masse(ijq, l)*(0.5*(q(ijq, &
631 l)+qd(ijq,l))*zsigd(ijq,l)+(zsig-zsigd(ijq,l))*(q(ijq, &
632 l)+zz*(qg(ijq,l)-q(ijq,l))))
633 END IF
634 ELSE
635 ijq = ij + 1
636 i = ijq - (j-1)*iip1
637 ! accumulation pour les mailles completements advectees
638 DO WHILE (-zu_m>masse(ijq,l))
639 u_mq(ij, l) = u_mq(ij, l) - q(ijq, l)*masse(ijq, l)
640 zu_m = zu_m + masse(ijq, l)
641 i = mod(i, iim) + 1
642 ijq = (j-1)*iip1 + i
643 END DO
644 ! ajout de la maille non completement advectee
645 ! 2eme MODIF SPECIFIQUE
646 zsig = -zu_m/masse(ij+1, l)
647 IF (zsig<=zsigg(ijq,l)) THEN
648 u_mq(ij, l) = u_mq(ij, l) + zu_m*(qg(ijq,l)-0.5*zsig/zsigg(ijq, &
649 l)*(qg(ijq,l)-q(ijq,l)))
650 ELSE
651 ! u_mq(ij,l)=u_mq(ij,l)+zu_m*q(ijq,l)
652 ! goto 9999
653 zz = 0.5*(zsig-zsigg(ijq,l))/zsigd(ijq, l)
654 IF (.NOT. (zz>0. .AND. zz<=0.5)) THEN
655 PRINT *, 'probleme22 au point ij=', ij, ' l=', l
656 PRINT *, 'zz=', zz
657 STOP
658 END IF
659 u_mq(ij, l) = u_mq(ij, l) - masse(ijq, l)*(0.5*(q(ijq, &
660 l)+qg(ijq,l))*zsigg(ijq,l)+(zsig-zsigg(ijq,l))*(q(ijq, &
661 l)+zz*(qd(ijq,l)-q(ijq,l))))
662 END IF
663 ! fin de la modif
664 END IF
665 END DO
666 END IF
667 END DO
668 END IF ! n0.gt.0
669
670 ! bouclage en latitude
671 DO l = 1, llm
672 DO ij = iip1 + iip1, ip1jm, iip1
673 u_mq(ij, l) = u_mq(ij-iim, l)
674 END DO
675 END DO
676
677 ! =================================================================
678 ! CALCUL DE LA CONVERGENCE DES FLUX
679 ! =================================================================
680
681 DO l = 1, llm
682 DO ij = iip2 + 1, ip1jm
683 new_m = masse(ij, l) + u_m(ij-1, l) - u_m(ij, l)
684 q(ij, l) = (q(ij,l)*masse(ij,l)+u_mq(ij-1,l)-u_mq(ij,l))/new_m
685 masse(ij, l) = new_m
686 END DO
687 ! Modif Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
688 DO ij = iip1 + iip1, ip1jm, iip1
689 q(ij-iim, l) = q(ij, l)
690 masse(ij-iim, l) = masse(ij, l)
691 END DO
692 END DO
693
694 RETURN
695 END SUBROUTINE advnx
696 SUBROUTINE advny(q, qs, qn, masse, v_m)
697
698 ! Auteur : F. Hourdin
699
700 ! ********************************************************************
701 ! Shema d'advection " pseudo amont " .
702 ! ********************************************************************
703 ! nq,iq,q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg ....
704
705
706 ! --------------------------------------------------------------------
707 USE dimens_m
708 USE paramet_m
709 USE comgeom
710 USE conf_gcm_m
711 IMPLICIT NONE
712
713
714
715 ! Arguments:
716 ! ----------
717 REAL masse(ip1jmp1, llm)
718 REAL v_m(ip1jm, llm)
719 REAL q(ip1jmp1, llm), qn(ip1jmp1, llm), qs(ip1jmp1, llm)
720
721 ! Local
722 ! ---------
723
724 INTEGER ij, l
725
726 REAL new_m, zdq, zz
727 REAL zsigs(ip1jmp1), zsign(ip1jmp1), zsig
728 REAL v_mq(ip1jm, llm)
729 REAL convpn, convps, convmpn, convmps, massen, masses
730 REAL zm, zq, zsigm, zsigp, zqm, zqp
731 REAL ssum
732 REAL prec
733 SAVE prec
734
735 DATA prec/1.E-15/
736
737 DO l = 1, llm
738 DO ij = 1, ip1jmp1
739 zdq = qn(ij, l) - qs(ij, l)
740 IF (abs(zdq)>prec) THEN
741 zsign(ij) = (q(ij,l)-qs(ij,l))/zdq
742 zsigs(ij) = 1. - zsign(ij)
743 ELSE
744 zsign(ij) = 0.5
745 zsigs(ij) = 0.5
746 END IF
747 END DO
748
749 ! calcul de la pente maximum dans la maille en valeur absolue
750
751 DO ij = 1, ip1jm
752 IF (v_m(ij,l)>=0.) THEN
753 zsigp = zsign(ij+iip1)
754 zsigm = zsigs(ij+iip1)
755 zqp = qn(ij+iip1, l)
756 zqm = qs(ij+iip1, l)
757 zm = masse(ij+iip1, l)
758 zq = q(ij+iip1, l)
759 ELSE
760 zsigm = zsign(ij)
761 zsigp = zsigs(ij)
762 zqm = qn(ij, l)
763 zqp = qs(ij, l)
764 zm = masse(ij, l)
765 zq = q(ij, l)
766 END IF
767 zsig = abs(v_m(ij,l))/zm
768 IF (zsig==0.) zsigp = 0.1
769 IF (zsig<=zsigp) THEN
770 v_mq(ij, l) = v_m(ij, l)*(zqp-0.5*zsig/zsigp*(zqp-zq))
771 ELSE
772 zz = 0.5*(zsig-zsigp)/zsigm
773 v_mq(ij, l) = sign(zm, v_m(ij,l))*(0.5*(zq+zqp)*zsigp+(zsig-zsigp)*( &
774 zq+zz*(zqm-zq)))
775 END IF
776 END DO
777 END DO
778
779 DO l = 1, llm
780 DO ij = iip2, ip1jm
781 new_m = masse(ij, l) + v_m(ij, l) - v_m(ij-iip1, l)
782 q(ij, l) = (q(ij,l)*masse(ij,l)+v_mq(ij,l)-v_mq(ij-iip1,l))/new_m
783 masse(ij, l) = new_m
784 END DO
785 ! .-. ancienne version
786 convpn = ssum(iim, v_mq(1,l), 1)
787 convmpn = ssum(iim, v_m(1,l), 1)
788 massen = ssum(iim, masse(1,l), 1)
789 new_m = massen + convmpn
790 q(1, l) = (q(1,l)*massen+convpn)/new_m
791 DO ij = 1, iip1
792 q(ij, l) = q(1, l)
793 masse(ij, l) = new_m*aire(ij)/apoln
794 END DO
795
796 convps = -ssum(iim, v_mq(ip1jm-iim,l), 1)
797 convmps = -ssum(iim, v_m(ip1jm-iim,l), 1)
798 masses = ssum(iim, masse(ip1jm+1,l), 1)
799 new_m = masses + convmps
800 q(ip1jm+1, l) = (q(ip1jm+1,l)*masses+convps)/new_m
801 DO ij = ip1jm + 1, ip1jmp1
802 q(ij, l) = q(ip1jm+1, l)
803 masse(ij, l) = new_m*aire(ij)/apols
804 END DO
805 END DO
806
807 RETURN
808 END SUBROUTINE advny
809 SUBROUTINE advnz(q, qh, qb, masse, w_m)
810
811 ! Auteurs: F.Hourdin
812
813 ! ********************************************************************
814 ! Shema d'advection " pseudo amont " .
815 ! b designe le bas et h le haut
816 ! il y a une correspondance entre le b en z et le d en x
817 ! ********************************************************************
818
819
820 ! --------------------------------------------------------------------
821 USE dimens_m
822 USE paramet_m
823 USE comgeom
824 USE conf_gcm_m
825 IMPLICIT NONE
826
827
828
829 ! Arguments:
830 ! ----------
831 REAL masse(ip1jmp1, llm)
832 REAL w_m(ip1jmp1, llm+1)
833 REAL q(ip1jmp1, llm), qb(ip1jmp1, llm), qh(ip1jmp1, llm)
834
835
836 ! Local
837 ! ---------
838
839 INTEGER ij, l
840
841 REAL new_m, zdq, zz
842 REAL zsigh(ip1jmp1, llm), zsigb(ip1jmp1, llm), zsig
843 REAL w_mq(ip1jmp1, llm+1)
844 REAL zm, zq, zsigm, zsigp, zqm, zqp
845 REAL prec
846 SAVE prec
847
848 DATA prec/1.E-13/
849
850 DO l = 1, llm
851 DO ij = 1, ip1jmp1
852 zdq = qb(ij, l) - qh(ij, l)
853 IF (abs(zdq)>prec) THEN
854 zsigb(ij, l) = (q(ij,l)-qh(ij,l))/zdq
855 zsigh(ij, l) = 1. - zsigb(ij, l)
856 zsigb(ij, l) = min(max(zsigb(ij,l),0.), 1.)
857 ELSE
858 zsigb(ij, l) = 0.5
859 zsigh(ij, l) = 0.5
860 END IF
861 END DO
862 END DO
863
864 ! calcul de la pente maximum dans la maille en valeur absolue
865 DO l = 2, llm
866 DO ij = 1, ip1jmp1
867 IF (w_m(ij,l)>=0.) THEN
868 zsigp = zsigb(ij, l)
869 zsigm = zsigh(ij, l)
870 zqp = qb(ij, l)
871 zqm = qh(ij, l)
872 zm = masse(ij, l)
873 zq = q(ij, l)
874 ELSE
875 zsigm = zsigb(ij, l-1)
876 zsigp = zsigh(ij, l-1)
877 zqm = qb(ij, l-1)
878 zqp = qh(ij, l-1)
879 zm = masse(ij, l-1)
880 zq = q(ij, l-1)
881 END IF
882 zsig = abs(w_m(ij,l))/zm
883 IF (zsig==0.) zsigp = 0.1
884 IF (zsig<=zsigp) THEN
885 w_mq(ij, l) = w_m(ij, l)*(zqp-0.5*zsig/zsigp*(zqp-zq))
886 ELSE
887 zz = 0.5*(zsig-zsigp)/zsigm
888 w_mq(ij, l) = sign(zm, w_m(ij,l))*(0.5*(zq+zqp)*zsigp+(zsig-zsigp)*( &
889 zq+zz*(zqm-zq)))
890 END IF
891 END DO
892 END DO
893
894 DO ij = 1, ip1jmp1
895 w_mq(ij, llm+1) = 0.
896 w_mq(ij, 1) = 0.
897 END DO
898
899 DO l = 1, llm
900 DO ij = 1, ip1jmp1
901 new_m = masse(ij, l) + w_m(ij, l+1) - w_m(ij, l)
902 q(ij, l) = (q(ij,l)*masse(ij,l)+w_mq(ij,l+1)-w_mq(ij,l))/new_m
903 masse(ij, l) = new_m
904 END DO
905 END DO
906
907 END SUBROUTINE advnz
  ViewVC Help
Powered by ViewVC 1.1.21