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

Contents of /trunk/dyn3d/bilan_dyn.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 39 - (show 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 !
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 histcom
13 use calendar
14 use histwrite_m
15 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 use nr_util, only: pi
24
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 print *,
181 . 'Pb : le pas de cumule doit etre multiple du pas'
182 print *,'dt_app=',dt_app
183 print *,'dt_cum=',dt_cum
184 stop
185 endif
186
187 if (i_sortie.eq.1) then
188 file='dynzon'
189 call inigrads(ifile ,(/0./),180./pi,0.,0.,rlatv,-90.,90.,
190 $ 180./pi ,presnivs,1. ,dt_cum,file,'dyn_zon ')
191 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 call histbeg_totreg('dynzon', rlong(:1), rlatg,
221 . 1, 1, 1, jjm,
222 . itau_dyn, zjulian, dt_cum, thoriid, fileid)
223
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 . print *,'var ',itr,iQ
251 . ,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 . 'ave(X)',dt_cum,dt_cum)
255 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 . 'ave(X)',dt_cum,dt_cum)
262 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 . 'ave(X)', dt_cum, dt_cum)
270 call histdef(fileid, 'v', 'v',
271 . 'm/s', 1, jjm, thoriid, llm, 1, llm, zvertiid,
272 . 'ave(X)', dt_cum, dt_cum)
273 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 . 'ave(X)',dt_cum,dt_cum)
278
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 . 'ave(X)',dt_cum,dt_cum)
287 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 . print *,'dans bilan_dyn ',icum,'->',icum+1
338 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 . print *,'Apres les calculs fait a chaque pas'
404 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 . print *,'Pas d ecriture'
412
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 call histwrite(fileid,znom(itr,iQ),itau,zvQ(:,:,itr,iQ))
532 enddo
533 call histwrite(fileid,'psi'//nom(iQ),itau,psiQ(:,1:llm,iQ))
534 enddo
535
536 call histwrite(fileid,'masse',itau,zmasse)
537 call histwrite(fileid,'v',itau,zv)
538 psi=psi*1.e-9
539 call histwrite(fileid,'psi',itau,psi(:,1:llm))
540
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 call histwrite(fileid,'a'//znom(itr,iQ),itau,zavQ(:,itr,iQ))
560 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