/[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 27 - (hide annotations)
Thu Mar 25 14:29:07 2010 UTC (14 years, 2 months ago) by guez
File size: 16337 byte(s)
"dyn3d" and "filtrez" do not contain any included file so make rules
have been updated.

"comdissip.f90" was useless, removed it.

"dynredem0" wrote undefined value in "controle(31)", that was
overwritten by "dynredem1". Now "dynredem0" just writes 0 to
"controle(31)".

Removed arguments of "inidissip". "inidissip" now accesses the
variables by use association.

In program "etat0_lim", "itaufin" is not defined so "dynredem1" wrote
undefined value to "controle(31)". Added argument "itau" of
"dynredem1" to correct that.

"itaufin" does not need to be a module variable (of "temps"), made it
a local variable of "leapfrog".

Removed calls to "diagedyn" from "leapfrog".

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     integer fileid
145     integer thoriid, zvertiid
146     save fileid
147    
148     integer ndex3d(jjm*llm)
149    
150     C Variables locales
151     C
152     real zjulian
153     character*3 str
154     character*10 ctrac
155     integer ii,jj
156     integer zan, dayref
157     C
158     real rlong(jjm),rlatg(jjm)
159    
160     !!print *, "Call sequence information: bilan_dyn"
161    
162     c=====================================================================
163     c Initialisation
164     c=====================================================================
165    
166     time=time+dt_app
167     itau=itau+1
168    
169     if (first) then
170    
171    
172     icum=0
173     c initialisation des fichiers
174     first=.false.
175     c ncum est la frequence de stokage en pas de temps
176     ncum=dt_cum/dt_app
177     if (abs(ncum*dt_app-dt_cum).gt.1.e-5*dt_app) then
178 guez 12 print *,
179 guez 3 . 'Pb : le pas de cumule doit etre multiple du pas'
180 guez 12 print *,'dt_app=',dt_app
181     print *,'dt_cum=',dt_cum
182 guez 3 stop
183     endif
184    
185     if (i_sortie.eq.1) then
186     file='dynzon'
187 guez 20 call inigrads(ifile ,(/0./),180./pi,0.,0.,rlatv,-90.,90.,
188     $ 180./pi ,presnivs,1. ,dt_cum,file,'dyn_zon ')
189 guez 3 endif
190    
191     nom(itemp)='T'
192     nom(igeop)='gz'
193     nom(iecin)='K'
194     nom(iang)='ang'
195     nom(iu)='u'
196     nom(iovap)='ovap'
197     nom(iun)='un'
198    
199     unites(itemp)='K'
200     unites(igeop)='m2/s2'
201     unites(iecin)='m2/s2'
202     unites(iang)='ang'
203     unites(iu)='m/s'
204     unites(iovap)='kg/kg'
205     unites(iun)='un'
206    
207    
208     c Initialisation du fichier contenant les moyennes zonales.
209     c ---------------------------------------------------------
210    
211     zan = annee_ref
212     dayref = day_ref
213     CALL ymds2ju(zan, 1, dayref, 0.0, zjulian)
214    
215     rlong=0.
216     rlatg=rlatv*180./pi
217    
218 guez 27 call histbeg_totreg('dynzon', rlong(:1), rlatg,
219 guez 3 . 1, 1, 1, jjm,
220 guez 27 . itau_dyn, zjulian, dt_cum, thoriid, fileid)
221 guez 3
222     C
223     C Appel a histvert pour la grille verticale
224     C
225     call histvert(fileid, 'presnivs', 'Niveaux sigma','mb',
226     . llm, presnivs, zvertiid)
227     C
228     C Appels a histdef pour la definition des variables a sauvegarder
229     do iQ=1,nQ
230     do itr=1,ntr
231     if(itr.eq.1) then
232     znom(itr,iQ)=nom(iQ)
233     znoml(itr,iQ)=nom(iQ)
234     zunites(itr,iQ)=unites(iQ)
235     else
236     znom(itr,iQ)=ctrs(itr)//'v'//nom(iQ)
237     znoml(itr,iQ)='transport : v * '//nom(iQ)//' '//ctrs(itr)
238     zunites(itr,iQ)='m/s * '//unites(iQ)
239     endif
240     enddo
241     enddo
242    
243     c Declarations des champs avec dimension verticale
244     c print*,'1HISTDEF'
245     do iQ=1,nQ
246     do itr=1,ntr
247     IF (prt_level > 5)
248 guez 12 . print *,'var ',itr,iQ
249 guez 3 . ,znom(itr,iQ),znoml(itr,iQ),zunites(itr,iQ)
250     call histdef(fileid,znom(itr,iQ),znoml(itr,iQ),
251     . zunites(itr,iQ),1,jjm,thoriid,llm,1,llm,zvertiid,
252 guez 15 . 'ave(X)',dt_cum,dt_cum)
253 guez 3 enddo
254     c Declarations pour les fonctions de courant
255     c print*,'2HISTDEF'
256     call histdef(fileid,'psi'//nom(iQ)
257     . ,'stream fn. '//znoml(itot,iQ),
258     . zunites(itot,iQ),1,jjm,thoriid,llm,1,llm,zvertiid,
259 guez 15 . 'ave(X)',dt_cum,dt_cum)
260 guez 3 enddo
261    
262    
263     c Declarations pour les champs de transport d'air
264     c print*,'3HISTDEF'
265     call histdef(fileid, 'masse', 'masse',
266     . 'kg', 1, jjm, thoriid, llm, 1, llm, zvertiid,
267 guez 15 . 'ave(X)', dt_cum, dt_cum)
268 guez 3 call histdef(fileid, 'v', 'v',
269     . 'm/s', 1, jjm, thoriid, llm, 1, llm, zvertiid,
270 guez 15 . 'ave(X)', dt_cum, dt_cum)
271 guez 3 c Declarations pour les fonctions de courant
272     c print*,'4HISTDEF'
273     call histdef(fileid,'psi','stream fn. MMC ','mega t/s',
274     . 1,jjm,thoriid,llm,1,llm,zvertiid,
275 guez 15 . 'ave(X)',dt_cum,dt_cum)
276 guez 3
277    
278     c Declaration des champs 1D de transport en latitude
279     c print*,'5HISTDEF'
280     do iQ=1,nQ
281     do itr=2,ntr
282     call histdef(fileid,'a'//znom(itr,iQ),znoml(itr,iQ),
283     . zunites(itr,iQ),1,jjm,thoriid,1,1,1,-99,
284 guez 15 . 'ave(X)',dt_cum,dt_cum)
285 guez 3 enddo
286     enddo
287    
288    
289     c print*,'8HISTDEF'
290     CALL histend(fileid)
291    
292    
293     endif
294    
295    
296     c=====================================================================
297     c Calcul des champs dynamiques
298     c ----------------------------
299    
300     c énergie cinétique
301     ucont(:,:,:)=0
302     CALL covcont(llm,ucov,vcov,ucont,vcont)
303     CALL enercin(vcov,ucov,vcont,ucont,ecin)
304    
305     c moment cinétique
306     do l=1,llm
307     ang(:,:,l)=ucov(:,:,l)+constang_2d(:,:)
308     unat(:,:,l)=ucont(:,:,l)*cu_2d(:,:)
309     enddo
310    
311     Q(:,:,:,itemp)=teta(:,:,:)*pk(:,:,:)/cpp
312     Q(:,:,:,igeop)=phi(:,:,:)
313     Q(:,:,:,iecin)=ecin(:,:,:)
314     Q(:,:,:,iang)=ang(:,:,:)
315     Q(:,:,:,iu)=unat(:,:,:)
316     Q(:,:,:,iovap)=trac(:,:,:,1)
317     Q(:,:,:,iun)=1.
318    
319    
320     c=====================================================================
321     c Cumul
322     c=====================================================================
323     c
324     if(icum.EQ.0) then
325     ps_cum=0.
326     masse_cum=0.
327     flux_u_cum=0.
328     flux_v_cum=0.
329     Q_cum=0.
330     flux_vQ_cum=0.
331     flux_uQ_cum=0.
332     endif
333    
334     IF (prt_level > 5)
335 guez 12 . print *,'dans bilan_dyn ',icum,'->',icum+1
336 guez 3 icum=icum+1
337    
338     c accumulation des flux de masse horizontaux
339     ps_cum=ps_cum+ps
340     masse_cum=masse_cum+masse
341     flux_u_cum=flux_u_cum+flux_u
342     flux_v_cum=flux_v_cum+flux_v
343     do iQ=1,nQ
344     Q_cum(:,:,:,iQ)=Q_cum(:,:,:,iQ)+Q(:,:,:,iQ)*masse(:,:,:)
345     enddo
346    
347     c=====================================================================
348     c FLUX ET TENDANCES
349     c=====================================================================
350    
351     c Flux longitudinal
352     c -----------------
353     do iQ=1,nQ
354     do l=1,llm
355     do j=1,jjp1
356     do i=1,iim
357     flux_uQ_cum(i,j,l,iQ)=flux_uQ_cum(i,j,l,iQ)
358     s +flux_u(i,j,l)*0.5*(Q(i,j,l,iQ)+Q(i+1,j,l,iQ))
359     enddo
360     flux_uQ_cum(iip1,j,l,iQ)=flux_uQ_cum(1,j,l,iQ)
361     enddo
362     enddo
363     enddo
364    
365     c flux méridien
366     c -------------
367     do iQ=1,nQ
368     do l=1,llm
369     do j=1,jjm
370     do i=1,iip1
371     flux_vQ_cum(i,j,l,iQ)=flux_vQ_cum(i,j,l,iQ)
372     s +flux_v(i,j,l)*0.5*(Q(i,j,l,iQ)+Q(i,j+1,l,iQ))
373     enddo
374     enddo
375     enddo
376     enddo
377    
378    
379     c tendances
380     c ---------
381    
382     c convergence horizontale
383     call convflu(flux_uQ_cum,flux_vQ_cum,llm*nQ,dQ)
384    
385     c calcul de la vitesse verticale
386     call convmas(flux_u_cum,flux_v_cum,convm)
387     CALL vitvert(convm,w)
388    
389     do iQ=1,nQ
390     do l=1,llm-1
391     do j=1,jjp1
392     do i=1,iip1
393     ww=-0.5*w(i,j,l+1)*(Q(i,j,l,iQ)+Q(i,j,l+1,iQ))
394     dQ(i,j,l ,iQ)=dQ(i,j,l ,iQ)-ww
395     dQ(i,j,l+1,iQ)=dQ(i,j,l+1,iQ)+ww
396     enddo
397     enddo
398     enddo
399     enddo
400     IF (prt_level > 5)
401 guez 12 . print *,'Apres les calculs fait a chaque pas'
402 guez 3 c=====================================================================
403     c PAS DE TEMPS D'ECRITURE
404     c=====================================================================
405     if (icum.eq.ncum) then
406     c=====================================================================
407    
408     IF (prt_level > 5)
409 guez 12 . print *,'Pas d ecriture'
410 guez 3
411     c Normalisation
412     do iQ=1,nQ
413     Q_cum(:,:,:,iQ)=Q_cum(:,:,:,iQ)/masse_cum(:,:,:)
414     enddo
415     zz=1./float(ncum)
416     ps_cum=ps_cum*zz
417     masse_cum=masse_cum*zz
418     flux_u_cum=flux_u_cum*zz
419     flux_v_cum=flux_v_cum*zz
420     flux_uQ_cum=flux_uQ_cum*zz
421     flux_vQ_cum=flux_vQ_cum*zz
422     dQ=dQ*zz
423    
424    
425     c A retravailler eventuellement
426     c division de dQ par la masse pour revenir aux bonnes grandeurs
427     do iQ=1,nQ
428     dQ(:,:,:,iQ)=dQ(:,:,:,iQ)/masse_cum(:,:,:)
429     enddo
430    
431     c=====================================================================
432     c Transport méridien
433     c=====================================================================
434    
435     c cumul zonal des masses des mailles
436     c ----------------------------------
437     zv=0.
438     zmasse=0.
439     call massbar(masse_cum,massebx,masseby)
440     do l=1,llm
441     do j=1,jjm
442     do i=1,iim
443     zmasse(j,l)=zmasse(j,l)+masseby(i,j,l)
444     zv(j,l)=zv(j,l)+flux_v_cum(i,j,l)
445     enddo
446     zfactv(j,l)=cv_2d(1,j)/zmasse(j,l)
447     enddo
448     enddo
449    
450     c print*,'3OK'
451     c --------------------------------------------------------------
452     c calcul de la moyenne zonale du transport :
453     c ------------------------------------------
454     c
455     c --
456     c TOT : la circulation totale [ vq ]
457     c
458     c - -
459     c MMC : mean meridional circulation [ v ] [ q ]
460     c
461     c ---- -- - -
462     c TRS : transitoires [ v'q'] = [ vq ] - [ v q ]
463     c
464     c - * - * - - - -
465     c STT : stationaires [ v q ] = [ v q ] - [ v ] [ q ]
466     c
467     c - -
468     c on utilise aussi l'intermediaire TMP : [ v q ]
469     c
470     c la variable zfactv transforme un transport meridien cumule
471     c en kg/s * unte-du-champ-transporte en m/s * unite-du-champ-transporte
472     c
473     c --------------------------------------------------------------
474    
475    
476     c ----------------------------------------
477     c Transport dans le plan latitude-altitude
478     c ----------------------------------------
479    
480     zvQ=0.
481     psiQ=0.
482     do iQ=1,nQ
483     zvQtmp=0.
484     do l=1,llm
485     do j=1,jjm
486     c print*,'j,l,iQ=',j,l,iQ
487     c Calcul des moyennes zonales du transort total et de zvQtmp
488     do i=1,iim
489     zvQ(j,l,itot,iQ)=zvQ(j,l,itot,iQ)
490     s +flux_vQ_cum(i,j,l,iQ)
491     zqy= 0.5*(Q_cum(i,j,l,iQ)*masse_cum(i,j,l)+
492     s Q_cum(i,j+1,l,iQ)*masse_cum(i,j+1,l))
493     zvQtmp(j,l)=zvQtmp(j,l)+flux_v_cum(i,j,l)*zqy
494     s /(0.5*(masse_cum(i,j,l)+masse_cum(i,j+1,l)))
495     zvQ(j,l,iave,iQ)=zvQ(j,l,iave,iQ)+zqy
496     enddo
497     c print*,'aOK'
498     c Decomposition
499     zvQ(j,l,iave,iQ)=zvQ(j,l,iave,iQ)/zmasse(j,l)
500     zvQ(j,l,itot,iQ)=zvQ(j,l,itot,iQ)*zfactv(j,l)
501     zvQtmp(j,l)=zvQtmp(j,l)*zfactv(j,l)
502     zvQ(j,l,immc,iQ)=zv(j,l)*zvQ(j,l,iave,iQ)*zfactv(j,l)
503     zvQ(j,l,itrs,iQ)=zvQ(j,l,itot,iQ)-zvQtmp(j,l)
504     zvQ(j,l,istn,iQ)=zvQtmp(j,l)-zvQ(j,l,immc,iQ)
505     enddo
506     enddo
507     c fonction de courant meridienne pour la quantite Q
508     do l=llm,1,-1
509     do j=1,jjm
510     psiQ(j,l,iQ)=psiQ(j,l+1,iQ)+zvQ(j,l,itot,iQ)
511     enddo
512     enddo
513     enddo
514    
515     c fonction de courant pour la circulation meridienne moyenne
516     psi=0.
517     do l=llm,1,-1
518     do j=1,jjm
519     psi(j,l)=psi(j,l+1)+zv(j,l)
520     zv(j,l)=zv(j,l)*zfactv(j,l)
521     enddo
522     enddo
523    
524     c print*,'4OK'
525     c sorties proprement dites
526     if (i_sortie.eq.1) then
527     do iQ=1,nQ
528     do itr=1,ntr
529 guez 15 call histwrite(fileid,znom(itr,iQ),itau,zvQ(:,:,itr,iQ))
530 guez 3 enddo
531 guez 15 call histwrite(fileid,'psi'//nom(iQ),itau,psiQ(:,1:llm,iQ))
532 guez 3 enddo
533    
534 guez 15 call histwrite(fileid,'masse',itau,zmasse)
535     call histwrite(fileid,'v',itau,zv)
536 guez 3 psi=psi*1.e-9
537 guez 15 call histwrite(fileid,'psi',itau,psi(:,1:llm))
538 guez 3
539     endif
540    
541    
542     c -----------------
543     c Moyenne verticale
544     c -----------------
545    
546     zamasse=0.
547     do l=1,llm
548     zamasse(:)=zamasse(:)+zmasse(:,l)
549     enddo
550     zavQ=0.
551     do iQ=1,nQ
552     do itr=2,ntr
553     do l=1,llm
554     zavQ(:,itr,iQ)=zavQ(:,itr,iQ)+zvQ(:,l,itr,iQ)*zmasse(:,l)
555     enddo
556     zavQ(:,itr,iQ)=zavQ(:,itr,iQ)/zamasse(:)
557 guez 15 call histwrite(fileid,'a'//znom(itr,iQ),itau,zavQ(:,itr,iQ))
558 guez 3 enddo
559     enddo
560    
561     c on doit pouvoir tracer systematiquement la fonction de courant.
562    
563     c=====================================================================
564     c/////////////////////////////////////////////////////////////////////
565     icum=0 !///////////////////////////////////////
566     endif ! icum.eq.ncum !///////////////////////////////////////
567     c/////////////////////////////////////////////////////////////////////
568     c=====================================================================
569    
570     return
571     end

  ViewVC Help
Powered by ViewVC 1.1.21