/[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 12 - (show annotations)
Mon Jul 21 16:05:07 2008 UTC (15 years, 9 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 !
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
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 . print *,'var ',itr,iQ
257 . ,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 . print *,'dans bilan_dyn ',icum,'->',icum+1
344 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 . print *,'Apres les calculs fait a chaque pas'
410 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 . print *,'Pas d ecriture'
418
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