source: LMDZ4/trunk/libf/dyn3d/bilan_dyn.F @ 604

Last change on this file since 604 was 604, checked in by Laurent Fairhead, 19 years ago

Repere par GG
LF

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