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