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

Contents of /trunk/libf/dyn3d/bilan_dyn.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 20 - (show annotations)
Wed Oct 15 16:19:57 2008 UTC (15 years, 7 months ago) by guez
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 !
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 print *,
182 . 'Pb : le pas de cumule doit etre multiple du pas'
183 print *,'dt_app=',dt_app
184 print *,'dt_cum=',dt_cum
185 stop
186 endif
187
188 if (i_sortie.eq.1) then
189 file='dynzon'
190 call inigrads(ifile ,(/0./),180./pi,0.,0.,rlatv,-90.,90.,
191 $ 180./pi ,presnivs,1. ,dt_cum,file,'dyn_zon ')
192 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 call histbeg_totreg(infile, rlong(:1), rlatg,
225 . 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 . print *,'var ',itr,iQ
255 . ,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 . 'ave(X)',dt_cum,dt_cum)
259 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 . 'ave(X)',dt_cum,dt_cum)
266 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 . 'ave(X)', dt_cum, dt_cum)
274 call histdef(fileid, 'v', 'v',
275 . 'm/s', 1, jjm, thoriid, llm, 1, llm, zvertiid,
276 . 'ave(X)', dt_cum, dt_cum)
277 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 . 'ave(X)',dt_cum,dt_cum)
282
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 . 'ave(X)',dt_cum,dt_cum)
291 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 . print *,'dans bilan_dyn ',icum,'->',icum+1
342 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 . print *,'Apres les calculs fait a chaque pas'
408 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 . print *,'Pas d ecriture'
416
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 call histwrite(fileid,znom(itr,iQ),itau,zvQ(:,:,itr,iQ))
536 enddo
537 call histwrite(fileid,'psi'//nom(iQ),itau,psiQ(:,1:llm,iQ))
538 enddo
539
540 call histwrite(fileid,'masse',itau,zmasse)
541 call histwrite(fileid,'v',itau,zv)
542 psi=psi*1.e-9
543 call histwrite(fileid,'psi',itau,psi(:,1:llm))
544
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 call histwrite(fileid,'a'//znom(itr,iQ),itau,zavQ(:,itr,iQ))
564 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