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

Annotation of /trunk/dyn3d/advn.f

Parent Directory Parent Directory | Revision Log Revision Log

Revision 112 - (hide annotations)
Thu Sep 18 13:36:51 2014 UTC (9 years, 8 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 guez 3
2 guez 81 ! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/advn.F,v 1.1.1.1 2004/05/19
3     ! 12:53:06 lmdzadmin Exp $
4 guez 3
5 guez 81 SUBROUTINE advn(q, masse, w, pbaru, pbarv, pdt, mode)
6 guez 3
7 guez 81 ! Auteur : F. Hourdin
8 guez 3
9 guez 81 ! ********************************************************************
10     ! Shema d'advection " pseudo amont " .
11     ! ********************************************************************
12     ! q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg ....
13 guez 3
14 guez 81 ! pbaru,pbarv,w flux de masse en u ,v ,w
15     ! pdt pas de temps
16 guez 3
17 guez 81 ! --------------------------------------------------------------------
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 guez 3
26    
27    
28 guez 81 ! 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 guez 3
36 guez 81 ! Local
37     ! ---------
38 guez 3
39 guez 112 INTEGER ij, l
40     REAL zm(ip1jmp1, llm)
41 guez 81 REAL mu(ip1jmp1, llm)
42     REAL mv(ip1jm, llm)
43     REAL mw(ip1jmp1, llm+1)
44 guez 112 REAL zq(ip1jmp1, llm), qpn, qps
45 guez 81 REAL zqg(ip1jmp1, llm), zqd(ip1jmp1, llm)
46     REAL zqs(ip1jmp1, llm), zqn(ip1jmp1, llm)
47     REAL zqh(ip1jmp1, llm), zqb(ip1jmp1, llm)
48 guez 112 REAL temps1, temps2, temps3
49     REAL ssum
50 guez 81 LOGICAL testcpu
51     SAVE testcpu
52     SAVE temps1, temps2, temps3
53     REAL zzpbar, zzw
54 guez 3
55 guez 81 REAL qmin, qmax
56     DATA qmin, qmax/0., 1./
57     DATA testcpu/.FALSE./
58     DATA temps1, temps2, temps3/0., 0., 0./
59 guez 3
60 guez 81 zzpbar = 0.5*pdt
61     zzw = pdt
62 guez 3
63 guez 81 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 guez 3
75 guez 81 DO ij = 1, ip1jmp1
76     mw(ij, llm+1) = 0.
77     END DO
78 guez 3
79 guez 81 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 guez 3
94 guez 81 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 guez 3
104 guez 81 ! 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 guez 3
118 guez 81 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 guez 3
127 guez 81 RETURN
128     END SUBROUTINE advn
129 guez 3
130 guez 81 SUBROUTINE advnqx(q, qg, qd)
131 guez 3
132 guez 81 ! Auteurs: Calcul des valeurs de q aux point u.
133 guez 3
134 guez 81 ! --------------------------------------------------------------------
135     USE dimens_m
136     USE paramet_m
137     USE conf_gcm_m
138     IMPLICIT NONE
139 guez 3
140    
141    
142 guez 81 ! Arguments:
143     ! ----------
144     REAL q(ip1jmp1, llm), qg(ip1jmp1, llm), qd(ip1jmp1, llm)
145 guez 3
146 guez 81 ! Local
147     ! ---------
148 guez 3
149 guez 81 INTEGER ij, l
150 guez 3
151 guez 81 REAL dxqu(ip1jmp1), zqu(ip1jmp1)
152     REAL zqmax(ip1jmp1), zqmin(ip1jmp1)
153     LOGICAL extremum(ip1jmp1)
154 guez 3
155 guez 81 INTEGER mode
156     SAVE mode
157     DATA mode/1/
158 guez 3
159 guez 81 ! 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 guez 3
191 guez 81 ! calcul des valeurs max et min acceptees aux interfaces
192 guez 3
193 guez 81 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 guez 3
224 guez 81 GO TO 8888
225 guez 3
226 guez 81 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 guez 3
230 guez 81 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 guez 3
237 guez 81 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 guez 3
247 guez 81 ! Auteurs: Calcul des valeurs de q aux point v.
248 guez 3
249 guez 81 ! --------------------------------------------------------------------
250     USE dimens_m
251     USE paramet_m
252     USE conf_gcm_m
253     IMPLICIT NONE
254 guez 3
255    
256    
257 guez 81 ! Arguments:
258     ! ----------
259     REAL q(ip1jmp1, llm), qs(ip1jmp1, llm), qn(ip1jmp1, llm)
260 guez 3
261 guez 81 ! Local
262     ! ---------
263 guez 3
264 guez 81 INTEGER ij, l
265 guez 3
266 guez 81 REAL dyqv(ip1jm), zqv(ip1jm, llm)
267     REAL zqmax(ip1jm), zqmin(ip1jm)
268     LOGICAL extremum(ip1jmp1)
269 guez 3
270 guez 81 INTEGER mode
271     SAVE mode
272     DATA mode/1/
273 guez 3
274 guez 81 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 guez 3
283 guez 81 ! 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 guez 3
290 guez 81 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 guez 3
295 guez 81 DO ij = iip2, ip1jm
296     extremum(ij) = dyqv(ij)*dyqv(ij-iip1) <= 0.
297     END DO
298 guez 3
299 guez 81 ! 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 guez 3
307 guez 81 ! 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 guez 3
313 guez 81 DO ij = 1, ip1jm
314     zqv(ij, l) = min(max(zqmin(ij),zqv(ij,l)), zqmax(ij))
315     END DO
316 guez 3
317 guez 81 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 guez 3
329 guez 81 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 guez 3
336 guez 81 END DO
337     END IF
338     RETURN
339     END SUBROUTINE advnqy
340 guez 3
341 guez 81 SUBROUTINE advnqz(q, qh, qb)
342 guez 3
343 guez 81 ! Auteurs: Calcul des valeurs de q aux point v.
344 guez 3
345 guez 81 ! --------------------------------------------------------------------
346     USE dimens_m
347     USE paramet_m
348     USE conf_gcm_m
349     IMPLICIT NONE
350 guez 3
351    
352    
353 guez 81 ! Arguments:
354     ! ----------
355     REAL q(ip1jmp1, llm), qh(ip1jmp1, llm), qb(ip1jmp1, llm)
356 guez 3
357 guez 81 ! Local
358     ! ---------
359 guez 3
360 guez 81 INTEGER ij, l
361 guez 3
362 guez 81 REAL dzqw(ip1jmp1, llm+1), zqw(ip1jmp1, llm+1)
363     REAL zqmax(ip1jmp1, llm), zqmin(ip1jmp1, llm)
364     LOGICAL extremum(ip1jmp1, llm)
365 guez 3
366 guez 81 INTEGER mode
367     SAVE mode
368 guez 3
369 guez 81 DATA mode/1/
370 guez 3
371 guez 81 ! calcul des pentes en u:
372     ! -----------------------
373 guez 3
374 guez 81 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 guez 3
403 guez 81 ! 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 guez 32
411 guez 81 ! 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