/[lmdze]/trunk/libf/dyn3d/bilan_dyn.f
ViewVC logotype

Annotation of /trunk/libf/dyn3d/bilan_dyn.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
File size: 16608 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/bilan_dyn.F,v 1.5 2005/03/16 10:12:17 fairhead Exp $
3     !
4     SUBROUTINE bilan_dyn (ntrac,dt_app,dt_cum,
5     s ps,masse,pk,flux_u,flux_v,teta,phi,ucov,vcov,trac)
6    
7     c AFAIRE
8     c Prevoir en champ nq+1 le diagnostique de l'energie
9     c en faisant Qzon=Cv T + L * ...
10     c vQ..A=Cp T + L * ...
11    
12     USE IOIPSL
13    
14     use dimens_m
15     use paramet_m
16     use comconst
17     use comvert
18     use comgeom, only: constang_2d, cu_2d, cv_2d, rlatv
19     use temps
20     use iniprint
21     use inigrads_m, only: inigrads
22    
23     IMPLICIT NONE
24    
25    
26     c====================================================================
27     c
28     c Sous-programme consacre à des diagnostics dynamiques de base
29     c
30     c
31     c De facon generale, les moyennes des scalaires Q sont ponderees par
32     c la masse.
33     c
34     c Les flux de masse sont eux simplement moyennes.
35     c
36     c====================================================================
37    
38     c Arguments :
39     c ===========
40    
41     integer ntrac
42     real dt_app,dt_cum
43     real ps(iip1,jjp1)
44     real masse(iip1,jjp1,llm),pk(iip1,jjp1,llm)
45     real flux_u(iip1,jjp1,llm)
46     real flux_v(iip1,jjm,llm)
47     real teta(iip1,jjp1,llm)
48     real phi(iip1,jjp1,llm)
49     real ucov(iip1,jjp1,llm)
50     real vcov(iip1,jjm,llm)
51     real trac(iip1,jjp1,llm,ntrac)
52    
53     c Local :
54     c =======
55    
56     integer icum,ncum
57     logical first
58     real zz,zqy,zfactv(jjm,llm)
59    
60     integer nQ
61     parameter (nQ=7)
62    
63    
64     cym character*6 nom(nQ)
65     cym character*6 unites(nQ)
66     character*6,save :: nom(nQ)
67     character*6,save :: unites(nQ)
68    
69     character*10 file
70     integer ifile
71     parameter (ifile=4)
72    
73     integer itemp,igeop,iecin,iang,iu,iovap,iun
74     integer i_sortie
75    
76     save first,icum,ncum
77     save itemp,igeop,iecin,iang,iu,iovap,iun
78     save i_sortie
79    
80     real time
81     integer itau
82     save time,itau
83     data time,itau/0.,0/
84    
85     data first/.true./
86     data itemp,igeop,iecin,iang,iu,iovap,iun/1,2,3,4,5,6,7/
87     data i_sortie/1/
88    
89     real ww
90    
91     c variables dynamiques intermédiaires
92     REAL vcont(iip1,jjm,llm),ucont(iip1,jjp1,llm)
93     REAL ang(iip1,jjp1,llm),unat(iip1,jjp1,llm)
94     REAL massebx(iip1,jjp1,llm),masseby(iip1,jjm,llm)
95     REAL vorpot(iip1,jjm,llm)
96     REAL w(iip1,jjp1,llm),ecin(iip1,jjp1,llm),convm(iip1,jjp1,llm)
97     REAL bern(iip1,jjp1,llm)
98    
99     c champ contenant les scalaires advectés.
100     real Q(iip1,jjp1,llm,nQ)
101    
102     c champs cumulés
103     real ps_cum(iip1,jjp1)
104     real masse_cum(iip1,jjp1,llm)
105     real flux_u_cum(iip1,jjp1,llm)
106     real flux_v_cum(iip1,jjm,llm)
107     real Q_cum(iip1,jjp1,llm,nQ)
108     real flux_uQ_cum(iip1,jjp1,llm,nQ)
109     real flux_vQ_cum(iip1,jjm,llm,nQ)
110     real flux_wQ_cum(iip1,jjp1,llm,nQ)
111     real dQ(iip1,jjp1,llm,nQ)
112    
113     save ps_cum,masse_cum,flux_u_cum,flux_v_cum
114     save Q_cum,flux_uQ_cum,flux_vQ_cum
115    
116     c champs de tansport en moyenne zonale
117     integer ntr,itr
118     parameter (ntr=5)
119    
120     cym character*10 znom(ntr,nQ)
121     cym character*20 znoml(ntr,nQ)
122     cym character*10 zunites(ntr,nQ)
123     character*10,save :: znom(ntr,nQ)
124     character*20,save :: znoml(ntr,nQ)
125     character*10,save :: zunites(ntr,nQ)
126    
127     integer iave,itot,immc,itrs,istn
128     data iave,itot,immc,itrs,istn/1,2,3,4,5/
129     character*3 ctrs(ntr)
130     data ctrs/' ','TOT','MMC','TRS','STN'/
131    
132     real zvQ(jjm,llm,ntr,nQ),zvQtmp(jjm,llm)
133     real zavQ(jjm,ntr,nQ),psiQ(jjm,llm+1,nQ)
134     real zmasse(jjm,llm),zamasse(jjm)
135    
136     real zv(jjm,llm),psi(jjm,llm+1)
137    
138     integer i,j,l,iQ
139    
140    
141     c Initialisation du fichier contenant les moyennes zonales.
142     c ---------------------------------------------------------
143    
144     character*10 infile
145    
146     integer fileid
147     integer thoriid, zvertiid
148     save fileid
149    
150     integer ndex3d(jjm*llm)
151    
152     C Variables locales
153     C
154     integer tau0
155     real zjulian
156     character*3 str
157     character*10 ctrac
158     integer ii,jj
159     integer zan, dayref
160     C
161     real rlong(jjm),rlatg(jjm)
162    
163     !!print *, "Call sequence information: bilan_dyn"
164    
165     c=====================================================================
166     c Initialisation
167     c=====================================================================
168    
169     time=time+dt_app
170     itau=itau+1
171    
172     if (first) then
173    
174    
175     icum=0
176     c initialisation des fichiers
177     first=.false.
178     c ncum est la frequence de stokage en pas de temps
179     ncum=dt_cum/dt_app
180     if (abs(ncum*dt_app-dt_cum).gt.1.e-5*dt_app) then
181 guez 12 print *,
182 guez 3 . 'Pb : le pas de cumule doit etre multiple du pas'
183 guez 12 print *,'dt_app=',dt_app
184     print *,'dt_cum=',dt_cum
185 guez 3 stop
186     endif
187    
188     if (i_sortie.eq.1) then
189     file='dynzon'
190     call inigrads(ifile
191     s ,(/0./),180./pi,0.,0.,rlatv,-90.,90.,180./pi
192     s ,presnivs,1.
193     s ,dt_cum,file,'dyn_zon ')
194     endif
195    
196     nom(itemp)='T'
197     nom(igeop)='gz'
198     nom(iecin)='K'
199     nom(iang)='ang'
200     nom(iu)='u'
201     nom(iovap)='ovap'
202     nom(iun)='un'
203    
204     unites(itemp)='K'
205     unites(igeop)='m2/s2'
206     unites(iecin)='m2/s2'
207     unites(iang)='ang'
208     unites(iu)='m/s'
209     unites(iovap)='kg/kg'
210     unites(iun)='un'
211    
212    
213     c Initialisation du fichier contenant les moyennes zonales.
214     c ---------------------------------------------------------
215    
216     infile='dynzon'
217    
218     zan = annee_ref
219     dayref = day_ref
220     CALL ymds2ju(zan, 1, dayref, 0.0, zjulian)
221     tau0 = itau_dyn
222    
223     rlong=0.
224     rlatg=rlatv*180./pi
225    
226     call histbeg_totreg(infile, 1, rlong(:1), jjm, rlatg,
227     . 1, 1, 1, jjm,
228     . tau0, zjulian, dt_cum, thoriid, fileid)
229    
230     C
231     C Appel a histvert pour la grille verticale
232     C
233     call histvert(fileid, 'presnivs', 'Niveaux sigma','mb',
234     . llm, presnivs, zvertiid)
235     C
236     C Appels a histdef pour la definition des variables a sauvegarder
237     do iQ=1,nQ
238     do itr=1,ntr
239     if(itr.eq.1) then
240     znom(itr,iQ)=nom(iQ)
241     znoml(itr,iQ)=nom(iQ)
242     zunites(itr,iQ)=unites(iQ)
243     else
244     znom(itr,iQ)=ctrs(itr)//'v'//nom(iQ)
245     znoml(itr,iQ)='transport : v * '//nom(iQ)//' '//ctrs(itr)
246     zunites(itr,iQ)='m/s * '//unites(iQ)
247     endif
248     enddo
249     enddo
250    
251     c Declarations des champs avec dimension verticale
252     c print*,'1HISTDEF'
253     do iQ=1,nQ
254     do itr=1,ntr
255     IF (prt_level > 5)
256 guez 12 . print *,'var ',itr,iQ
257 guez 3 . ,znom(itr,iQ),znoml(itr,iQ),zunites(itr,iQ)
258     call histdef(fileid,znom(itr,iQ),znoml(itr,iQ),
259     . zunites(itr,iQ),1,jjm,thoriid,llm,1,llm,zvertiid,
260     . 32,'ave(X)',dt_cum,dt_cum)
261     enddo
262     c Declarations pour les fonctions de courant
263     c print*,'2HISTDEF'
264     call histdef(fileid,'psi'//nom(iQ)
265     . ,'stream fn. '//znoml(itot,iQ),
266     . zunites(itot,iQ),1,jjm,thoriid,llm,1,llm,zvertiid,
267     . 32,'ave(X)',dt_cum,dt_cum)
268     enddo
269    
270    
271     c Declarations pour les champs de transport d'air
272     c print*,'3HISTDEF'
273     call histdef(fileid, 'masse', 'masse',
274     . 'kg', 1, jjm, thoriid, llm, 1, llm, zvertiid,
275     . 32, 'ave(X)', dt_cum, dt_cum)
276     call histdef(fileid, 'v', 'v',
277     . 'm/s', 1, jjm, thoriid, llm, 1, llm, zvertiid,
278     . 32, 'ave(X)', dt_cum, dt_cum)
279     c Declarations pour les fonctions de courant
280     c print*,'4HISTDEF'
281     call histdef(fileid,'psi','stream fn. MMC ','mega t/s',
282     . 1,jjm,thoriid,llm,1,llm,zvertiid,
283     . 32,'ave(X)',dt_cum,dt_cum)
284    
285    
286     c Declaration des champs 1D de transport en latitude
287     c print*,'5HISTDEF'
288     do iQ=1,nQ
289     do itr=2,ntr
290     call histdef(fileid,'a'//znom(itr,iQ),znoml(itr,iQ),
291     . zunites(itr,iQ),1,jjm,thoriid,1,1,1,-99,
292     . 32,'ave(X)',dt_cum,dt_cum)
293     enddo
294     enddo
295    
296    
297     c print*,'8HISTDEF'
298     CALL histend(fileid)
299    
300    
301     endif
302    
303    
304     c=====================================================================
305     c Calcul des champs dynamiques
306     c ----------------------------
307    
308     c énergie cinétique
309     ucont(:,:,:)=0
310     CALL covcont(llm,ucov,vcov,ucont,vcont)
311     CALL enercin(vcov,ucov,vcont,ucont,ecin)
312    
313     c moment cinétique
314     do l=1,llm
315     ang(:,:,l)=ucov(:,:,l)+constang_2d(:,:)
316     unat(:,:,l)=ucont(:,:,l)*cu_2d(:,:)
317     enddo
318    
319     Q(:,:,:,itemp)=teta(:,:,:)*pk(:,:,:)/cpp
320     Q(:,:,:,igeop)=phi(:,:,:)
321     Q(:,:,:,iecin)=ecin(:,:,:)
322     Q(:,:,:,iang)=ang(:,:,:)
323     Q(:,:,:,iu)=unat(:,:,:)
324     Q(:,:,:,iovap)=trac(:,:,:,1)
325     Q(:,:,:,iun)=1.
326    
327    
328     c=====================================================================
329     c Cumul
330     c=====================================================================
331     c
332     if(icum.EQ.0) then
333     ps_cum=0.
334     masse_cum=0.
335     flux_u_cum=0.
336     flux_v_cum=0.
337     Q_cum=0.
338     flux_vQ_cum=0.
339     flux_uQ_cum=0.
340     endif
341    
342     IF (prt_level > 5)
343 guez 12 . print *,'dans bilan_dyn ',icum,'->',icum+1
344 guez 3 icum=icum+1
345    
346     c accumulation des flux de masse horizontaux
347     ps_cum=ps_cum+ps
348     masse_cum=masse_cum+masse
349     flux_u_cum=flux_u_cum+flux_u
350     flux_v_cum=flux_v_cum+flux_v
351     do iQ=1,nQ
352     Q_cum(:,:,:,iQ)=Q_cum(:,:,:,iQ)+Q(:,:,:,iQ)*masse(:,:,:)
353     enddo
354    
355     c=====================================================================
356     c FLUX ET TENDANCES
357     c=====================================================================
358    
359     c Flux longitudinal
360     c -----------------
361     do iQ=1,nQ
362     do l=1,llm
363     do j=1,jjp1
364     do i=1,iim
365     flux_uQ_cum(i,j,l,iQ)=flux_uQ_cum(i,j,l,iQ)
366     s +flux_u(i,j,l)*0.5*(Q(i,j,l,iQ)+Q(i+1,j,l,iQ))
367     enddo
368     flux_uQ_cum(iip1,j,l,iQ)=flux_uQ_cum(1,j,l,iQ)
369     enddo
370     enddo
371     enddo
372    
373     c flux méridien
374     c -------------
375     do iQ=1,nQ
376     do l=1,llm
377     do j=1,jjm
378     do i=1,iip1
379     flux_vQ_cum(i,j,l,iQ)=flux_vQ_cum(i,j,l,iQ)
380     s +flux_v(i,j,l)*0.5*(Q(i,j,l,iQ)+Q(i,j+1,l,iQ))
381     enddo
382     enddo
383     enddo
384     enddo
385    
386    
387     c tendances
388     c ---------
389    
390     c convergence horizontale
391     call convflu(flux_uQ_cum,flux_vQ_cum,llm*nQ,dQ)
392    
393     c calcul de la vitesse verticale
394     call convmas(flux_u_cum,flux_v_cum,convm)
395     CALL vitvert(convm,w)
396    
397     do iQ=1,nQ
398     do l=1,llm-1
399     do j=1,jjp1
400     do i=1,iip1
401     ww=-0.5*w(i,j,l+1)*(Q(i,j,l,iQ)+Q(i,j,l+1,iQ))
402     dQ(i,j,l ,iQ)=dQ(i,j,l ,iQ)-ww
403     dQ(i,j,l+1,iQ)=dQ(i,j,l+1,iQ)+ww
404     enddo
405     enddo
406     enddo
407     enddo
408     IF (prt_level > 5)
409 guez 12 . print *,'Apres les calculs fait a chaque pas'
410 guez 3 c=====================================================================
411     c PAS DE TEMPS D'ECRITURE
412     c=====================================================================
413     if (icum.eq.ncum) then
414     c=====================================================================
415    
416     IF (prt_level > 5)
417 guez 12 . print *,'Pas d ecriture'
418 guez 3
419     c Normalisation
420     do iQ=1,nQ
421     Q_cum(:,:,:,iQ)=Q_cum(:,:,:,iQ)/masse_cum(:,:,:)
422     enddo
423     zz=1./float(ncum)
424     ps_cum=ps_cum*zz
425     masse_cum=masse_cum*zz
426     flux_u_cum=flux_u_cum*zz
427     flux_v_cum=flux_v_cum*zz
428     flux_uQ_cum=flux_uQ_cum*zz
429     flux_vQ_cum=flux_vQ_cum*zz
430     dQ=dQ*zz
431    
432    
433     c A retravailler eventuellement
434     c division de dQ par la masse pour revenir aux bonnes grandeurs
435     do iQ=1,nQ
436     dQ(:,:,:,iQ)=dQ(:,:,:,iQ)/masse_cum(:,:,:)
437     enddo
438    
439     c=====================================================================
440     c Transport méridien
441     c=====================================================================
442    
443     c cumul zonal des masses des mailles
444     c ----------------------------------
445     zv=0.
446     zmasse=0.
447     call massbar(masse_cum,massebx,masseby)
448     do l=1,llm
449     do j=1,jjm
450     do i=1,iim
451     zmasse(j,l)=zmasse(j,l)+masseby(i,j,l)
452     zv(j,l)=zv(j,l)+flux_v_cum(i,j,l)
453     enddo
454     zfactv(j,l)=cv_2d(1,j)/zmasse(j,l)
455     enddo
456     enddo
457    
458     c print*,'3OK'
459     c --------------------------------------------------------------
460     c calcul de la moyenne zonale du transport :
461     c ------------------------------------------
462     c
463     c --
464     c TOT : la circulation totale [ vq ]
465     c
466     c - -
467     c MMC : mean meridional circulation [ v ] [ q ]
468     c
469     c ---- -- - -
470     c TRS : transitoires [ v'q'] = [ vq ] - [ v q ]
471     c
472     c - * - * - - - -
473     c STT : stationaires [ v q ] = [ v q ] - [ v ] [ q ]
474     c
475     c - -
476     c on utilise aussi l'intermediaire TMP : [ v q ]
477     c
478     c la variable zfactv transforme un transport meridien cumule
479     c en kg/s * unte-du-champ-transporte en m/s * unite-du-champ-transporte
480     c
481     c --------------------------------------------------------------
482    
483    
484     c ----------------------------------------
485     c Transport dans le plan latitude-altitude
486     c ----------------------------------------
487    
488     zvQ=0.
489     psiQ=0.
490     do iQ=1,nQ
491     zvQtmp=0.
492     do l=1,llm
493     do j=1,jjm
494     c print*,'j,l,iQ=',j,l,iQ
495     c Calcul des moyennes zonales du transort total et de zvQtmp
496     do i=1,iim
497     zvQ(j,l,itot,iQ)=zvQ(j,l,itot,iQ)
498     s +flux_vQ_cum(i,j,l,iQ)
499     zqy= 0.5*(Q_cum(i,j,l,iQ)*masse_cum(i,j,l)+
500     s Q_cum(i,j+1,l,iQ)*masse_cum(i,j+1,l))
501     zvQtmp(j,l)=zvQtmp(j,l)+flux_v_cum(i,j,l)*zqy
502     s /(0.5*(masse_cum(i,j,l)+masse_cum(i,j+1,l)))
503     zvQ(j,l,iave,iQ)=zvQ(j,l,iave,iQ)+zqy
504     enddo
505     c print*,'aOK'
506     c Decomposition
507     zvQ(j,l,iave,iQ)=zvQ(j,l,iave,iQ)/zmasse(j,l)
508     zvQ(j,l,itot,iQ)=zvQ(j,l,itot,iQ)*zfactv(j,l)
509     zvQtmp(j,l)=zvQtmp(j,l)*zfactv(j,l)
510     zvQ(j,l,immc,iQ)=zv(j,l)*zvQ(j,l,iave,iQ)*zfactv(j,l)
511     zvQ(j,l,itrs,iQ)=zvQ(j,l,itot,iQ)-zvQtmp(j,l)
512     zvQ(j,l,istn,iQ)=zvQtmp(j,l)-zvQ(j,l,immc,iQ)
513     enddo
514     enddo
515     c fonction de courant meridienne pour la quantite Q
516     do l=llm,1,-1
517     do j=1,jjm
518     psiQ(j,l,iQ)=psiQ(j,l+1,iQ)+zvQ(j,l,itot,iQ)
519     enddo
520     enddo
521     enddo
522    
523     c fonction de courant pour la circulation meridienne moyenne
524     psi=0.
525     do l=llm,1,-1
526     do j=1,jjm
527     psi(j,l)=psi(j,l+1)+zv(j,l)
528     zv(j,l)=zv(j,l)*zfactv(j,l)
529     enddo
530     enddo
531    
532     c print*,'4OK'
533     c sorties proprement dites
534     if (i_sortie.eq.1) then
535     do iQ=1,nQ
536     do itr=1,ntr
537     call histwrite(fileid,znom(itr,iQ),itau,zvQ(:,:,itr,iQ)
538     s ,jjm*llm,ndex3d)
539     enddo
540     call histwrite(fileid,'psi'//nom(iQ),itau,psiQ(:,1:llm,iQ)
541     s ,jjm*llm,ndex3d)
542     enddo
543    
544     call histwrite(fileid,'masse',itau,zmasse
545     s ,jjm*llm,ndex3d)
546     call histwrite(fileid,'v',itau,zv
547     s ,jjm*llm,ndex3d)
548     psi=psi*1.e-9
549     call histwrite(fileid,'psi',itau,psi(:,1:llm),jjm*llm,ndex3d)
550    
551     endif
552    
553    
554     c -----------------
555     c Moyenne verticale
556     c -----------------
557    
558     zamasse=0.
559     do l=1,llm
560     zamasse(:)=zamasse(:)+zmasse(:,l)
561     enddo
562     zavQ=0.
563     do iQ=1,nQ
564     do itr=2,ntr
565     do l=1,llm
566     zavQ(:,itr,iQ)=zavQ(:,itr,iQ)+zvQ(:,l,itr,iQ)*zmasse(:,l)
567     enddo
568     zavQ(:,itr,iQ)=zavQ(:,itr,iQ)/zamasse(:)
569     call histwrite(fileid,'a'//znom(itr,iQ),itau,zavQ(:,itr,iQ)
570     s ,jjm*llm,ndex3d)
571     enddo
572     enddo
573    
574     c on doit pouvoir tracer systematiquement la fonction de courant.
575    
576     c=====================================================================
577     c/////////////////////////////////////////////////////////////////////
578     icum=0 !///////////////////////////////////////
579     endif ! icum.eq.ncum !///////////////////////////////////////
580     c/////////////////////////////////////////////////////////////////////
581     c=====================================================================
582    
583     return
584     end

  ViewVC Help
Powered by ViewVC 1.1.21