/[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 30 - (show annotations)
Thu Apr 1 09:07:28 2010 UTC (14 years, 1 month ago) by guez
File size: 16378 byte(s)
Imported Source files of the external library "IOIPSL_Lionel" into
"libf/IOIPSL".

Split "cray.f90" into "scopy.f90" and "ssum.f90".

Rewrote "leapfrog" in order to have a clearer algorithmic structure.

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

  ViewVC Help
Powered by ViewVC 1.1.21