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

Annotation of /trunk/dyn3d/bilan_dyn.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 39 - (hide annotations)
Tue Jan 25 15:11:05 2011 UTC (13 years, 4 months ago) by guez
Original Path: trunk/libf/dyn3d/bilan_dyn.f
File size: 16406 byte(s)
"pi" comes from "nr_util". Removed subroutine "initialize" in module
"comconst".

Copied the content of "fxy_sin.h" into "fxysinus", instead of getting
it from an "include" line. Removed file "fxy_sin.h".

"ps" has rank 2 in "gcm" and "dynetat0".

Assumed-shape for argument "q" of "integrd".

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

  ViewVC Help
Powered by ViewVC 1.1.21