/[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 15 - (show annotations)
Fri Aug 1 15:24:12 2008 UTC (15 years, 9 months ago) by guez
File size: 16431 byte(s)
-- Minor modification of input/output:

Added variable "Sigma_O3_Royer" to "histday.nc". "ecrit_day" is not
modified in "physiq". Removed variables "pyu1", "pyv1", "ftsol1",
"ftsol2", "ftsol3", "ftsol4", "psrf1", "psrf2", "psrf3", "psrf4"
"mfu", "mfd", "en_u", "en_d", "de_d", "de_u", "coefh" from
"histrac.nc".

Variable "raz_date" of module "conf_gcm_m" has logical type instead of
integer type.

-- Should not change any result at run time:

Modified calls to "IOIPSL_Lionel" procedures because the interfaces of
these procedures have been simplified.

Changed name of variable in module "start_init_orog_m": "masque" to
"mask".

Created a module containing procedure "phyredem".

Removed arguments "punjours", "pdayref" and "ptimestep" of procedure
"iniphysiq".

Renamed procedure "gr_phy_write" to "gr_phy_write_2d". Created
procedure "gr_phy_write_3d".

Removed procedures "ini_undefstd", "moy_undefSTD", "calcul_STDlev",
"calcul_divers".

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, rlong(:1), 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 . '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 . '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 . 'ave(X)', dt_cum, dt_cum)
276 call histdef(fileid, 'v', 'v',
277 . 'm/s', 1, jjm, thoriid, llm, 1, llm, zvertiid,
278 . '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 . '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 . '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 enddo
539 call histwrite(fileid,'psi'//nom(iQ),itau,psiQ(:,1:llm,iQ))
540 enddo
541
542 call histwrite(fileid,'masse',itau,zmasse)
543 call histwrite(fileid,'v',itau,zv)
544 psi=psi*1.e-9
545 call histwrite(fileid,'psi',itau,psi(:,1:llm))
546
547 endif
548
549
550 c -----------------
551 c Moyenne verticale
552 c -----------------
553
554 zamasse=0.
555 do l=1,llm
556 zamasse(:)=zamasse(:)+zmasse(:,l)
557 enddo
558 zavQ=0.
559 do iQ=1,nQ
560 do itr=2,ntr
561 do l=1,llm
562 zavQ(:,itr,iQ)=zavQ(:,itr,iQ)+zvQ(:,l,itr,iQ)*zmasse(:,l)
563 enddo
564 zavQ(:,itr,iQ)=zavQ(:,itr,iQ)/zamasse(:)
565 call histwrite(fileid,'a'//znom(itr,iQ),itau,zavQ(:,itr,iQ))
566 enddo
567 enddo
568
569 c on doit pouvoir tracer systematiquement la fonction de courant.
570
571 c=====================================================================
572 c/////////////////////////////////////////////////////////////////////
573 icum=0 !///////////////////////////////////////
574 endif ! icum.eq.ncum !///////////////////////////////////////
575 c/////////////////////////////////////////////////////////////////////
576 c=====================================================================
577
578 return
579 end

  ViewVC Help
Powered by ViewVC 1.1.21