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

Annotation of /trunk/dyn3d/advn.f

Parent Directory Parent Directory | Revision Log Revision Log

Revision 12 - (hide annotations)
Mon Jul 21 16:05:07 2008 UTC (15 years, 10 months ago) by guez
Original Path: trunk/libf/dyn3d/advn.f
File size: 26057 byte(s) 
-- Minor modification of input/output:

Created procedure "read_logic". Variables of module "logic" are read
by "read_logic" instead of "conf_gcm". Variable "offline" of module
"conf_gcm" is read from namelist instead of "*.def".

Deleted arguments "dtime", "co2_ppm_etat0", "solaire_etat0",
"tabcntr0" and local variables "radpas", "tab_cntrl" of
"phyetat0". "phyetat0" does not read "controle" in "startphy.nc" any
longer. "phyetat0" now reads global attribute "itau_phy" from
"startphy.nc". "phyredem" does not create variable "controle" in
"startphy.nc" any longer. "phyredem" now writes global attribute
"itau_phy" of "startphy.nc". Deleted argument "tabcntr0" of
"printflag". Removed diagnostic messages written by "printflag" for
comparison of the variable "controle" of "startphy.nc" and the
variables read from "*.def" or namelist input.

-- Removing unwanted functionality:

Removed variable "lunout" from module "iniprint", replaced everywhere
by standard output.

Removed case "ocean == 'couple'" in "clmain", "interfsurf_hq" and
"physiq". Removed procedure "interfoce_cpl".

-- Should not change anything at run time:

Automated creation of graphs in documentation. More documentation on
input files.

Converted Fortran files to free format: "phyredem.f90", "printflag.f90".

Split module "clesphy" into "clesphys" and "clesphys2".

Removed variables "conser", "leapf", "forward", "apphys", "apdiss" and
"statcl" from module "logic". Added arguments "conser" to "advect",
"leapf" to "integrd". Added local variables "forward", "leapf",
"apphys", "conser", "apdiss" in "leapfrog".

Added intent attributes.

Deleted arguments "dtime" of "phyredem", "pdtime" of "flxdtdq", "sh"
of "phytrac", "dt" of "yamada".

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