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

Annotation of /trunk/dyn3d/bilan_dyn.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 20 - (hide annotations)
Wed Oct 15 16:19:57 2008 UTC (15 years, 7 months ago) by guez
Original Path: trunk/libf/dyn3d/bilan_dyn.f
File size: 16422 byte(s)
Deleted argument "presnivs" of "physiq", "ini_histhf", "ini_histhf3d",
"ini_histday", "ini_histins", "ini_histrac", "phytrac". Access it from
"comvert" instead.

Replaced calls to NetCDF Fortran 77 interface by calls to Fortran 90
interface or to NetCDF95.

Procedure "gr_phy_write_3d" now works with a variable of arbitrary
size in the second dimension.

Annotated use statements with "only" clause.

Replaced calls to NetCDF interface version 2 by calls to Fortran 90
interface in "guide.f90" and "read_reanalyse.f".

In "write_histrac", replaced calls to "gr_fi_ecrit" by calls to
"gr_phy_write_2d" and "gr_phy_write_3d".

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

  ViewVC Help
Powered by ViewVC 1.1.21