/[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 27 - (show annotations)
Thu Mar 25 14:29:07 2010 UTC (14 years, 2 months ago) by guez
File size: 16337 byte(s)
"dyn3d" and "filtrez" do not contain any included file so make rules
have been updated.

"comdissip.f90" was useless, removed it.

"dynredem0" wrote undefined value in "controle(31)", that was
overwritten by "dynredem1". Now "dynredem0" just writes 0 to
"controle(31)".

Removed arguments of "inidissip". "inidissip" now accesses the
variables by use association.

In program "etat0_lim", "itaufin" is not defined so "dynredem1" wrote
undefined value to "controle(31)". Added argument "itau" of
"dynredem1" to correct that.

"itaufin" does not need to be a module variable (of "temps"), made it
a local variable of "leapfrog".

Removed calls to "diagedyn" from "leapfrog".

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

  ViewVC Help
Powered by ViewVC 1.1.21