source: LMDZ4/branches/IPSL-CM4_IPCC_branch/libf/dyn3d/bilan_dyn.F @ 3536

Last change on this file since 3536 was 524, checked in by lmdzadmin, 21 years ago

Initial revision

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 15.8 KB
Line 
1!
2! $Header$
3!
4c
5c $Header$
6c
7      SUBROUTINE bilan_dyn (ntrac,dt_app,dt_cum,
8     s  ps,masse,pk,flux_u,flux_v,teta,phi,ucov,vcov,trac)
9
10c   AFAIRE
11c   Prevoir en champ nq+1 le diagnostique de l'energie
12c   en faisant Qzon=Cv T + L * ...
13c             vQ..A=Cp T + L * ...
14
15      USE IOIPSL
16
17      IMPLICIT NONE
18
19#include "dimensions.h"
20#include "paramet.h"
21#include "comconst.h"
22#include "comvert.h"
23#include "comgeom2.h"
24#include "temps.h"
25#include "iniprint.h"
26
27c====================================================================
28c
29c   Sous-programme consacre à des diagnostics dynamiques de base
30c
31c
32c   De facon generale, les moyennes des scalaires Q sont ponderees par
33c   la masse.
34c
35c   Les flux de masse sont eux simplement moyennes.
36c
37c====================================================================
38
39c   Arguments :
40c   ===========
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
54c   Local :
55c   =======
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      character*6 nom(nQ)
66      character*6 unites(nQ)
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
118      character*10 znom(ntr,nQ)
119      character*20 znoml(ntr,nQ)
120      character*10 zunites(ntr,nQ)
121      integer iave,itot,immc,itrs,istn
122      data iave,itot,immc,itrs,istn/1,2,3,4,5/
123      character*3 ctrs(ntr)
124      data ctrs/'  ','TOT','MMC','TRS','STN'/
125
126      real zvQ(jjm,llm,ntr,nQ),zvQtmp(jjm,llm)
127      real zavQ(jjm,ntr,nQ),psiQ(jjm,llm+1,nQ)
128      real zmasse(jjm,llm),zamasse(jjm)
129
130      real zv(jjm,llm),psi(jjm,llm+1)
131
132      integer i,j,l,iQ
133
134
135c   Initialisation du fichier contenant les moyennes zonales.
136c   ---------------------------------------------------------
137
138      character*10 infile
139
140      integer fileid
141      integer thoriid, zvertiid
142      save fileid
143
144      integer ndex3d(jjm*llm)
145
146C   Variables locales
147C
148      integer tau0
149      real zjulian
150      character*3 str
151      character*10 ctrac
152      integer ii,jj
153      integer zan, dayref
154C
155      real rlong(jjm),rlatg(jjm)
156
157
158
159c=====================================================================
160c   Initialisation
161c=====================================================================
162
163      time=time+dt_app
164      itau=itau+1
165
166      if (first) then
167
168
169        icum=0
170c       initialisation des fichiers
171        first=.false.
172c   ncum est la frequence de stokage en pas de temps
173        ncum=dt_cum/dt_app
174        if (abs(ncum*dt_app-dt_cum).gt.1.e-5*dt_app) then
175           WRITE(lunout,*)
176     .            'Pb : le pas de cumule doit etre multiple du pas'
177           WRITE(lunout,*)'dt_app=',dt_app
178           WRITE(lunout,*)'dt_cum=',dt_cum
179           stop
180        endif
181
182        if (i_sortie.eq.1) then
183         file='dynzon'
184         call inigrads(ifile,1
185     s  ,0.,180./pi,0.,0.,jjm,rlatv,-90.,90.,180./pi
186     s  ,llm,presnivs,1.
187     s  ,dt_cum,file,'dyn_zon ')
188        endif
189
190        nom(itemp)='T'
191        nom(igeop)='gz'
192        nom(iecin)='K'
193        nom(iang)='ang'
194        nom(iu)='u'
195        nom(iovap)='ovap'
196        nom(iun)='un'
197
198        unites(itemp)='K'
199        unites(igeop)='m2/s2'
200        unites(iecin)='m2/s2'
201        unites(iang)='ang'
202        unites(iu)='m/s'
203        unites(iovap)='kg/kg'
204        unites(iun)='un'
205
206
207c   Initialisation du fichier contenant les moyennes zonales.
208c   ---------------------------------------------------------
209
210      infile='dynzon'
211
212      zan = annee_ref
213      dayref = day_ref
214      CALL ymds2ju(zan, 1, dayref, 0.0, zjulian)
215      tau0 = itau_dyn
216     
217      rlong=0.
218      rlatg=rlatv*180./pi
219       
220      call histbeg(infile, 1, rlong, jjm, rlatg,
221     .             1, 1, 1, jjm,
222     .             tau0, zjulian, dt_cum, thoriid, fileid)
223
224C
225C  Appel a histvert pour la grille verticale
226C
227      call histvert(fileid, 'presnivs', 'Niveaux sigma','mb',
228     .              llm, presnivs, zvertiid)
229C
230C  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
245c   Declarations des champs avec dimension verticale
246c      print*,'1HISTDEF'
247      do iQ=1,nQ
248         do itr=1,ntr
249      IF (prt_level > 5)
250     . WRITE(lunout,*)'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     .        32,'ave(X)',dt_cum,dt_cum)
255         enddo
256c   Declarations pour les fonctions de courant
257c      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     .      32,'ave(X)',dt_cum,dt_cum)
262      enddo
263
264
265c   Declarations pour les champs de transport d'air
266c      print*,'3HISTDEF'
267      call histdef(fileid, 'masse', 'masse',
268     .             'kg', 1, jjm, thoriid, llm, 1, llm, zvertiid,
269     .             32, 'ave(X)', dt_cum, dt_cum)
270      call histdef(fileid, 'v', 'v',
271     .             'm/s', 1, jjm, thoriid, llm, 1, llm, zvertiid,
272     .             32, 'ave(X)', dt_cum, dt_cum)
273c   Declarations pour les fonctions de courant
274c      print*,'4HISTDEF'
275          call histdef(fileid,'psi','stream fn. MMC ','mega t/s',
276     .      1,jjm,thoriid,llm,1,llm,zvertiid,
277     .      32,'ave(X)',dt_cum,dt_cum)
278
279
280c   Declaration des champs 1D de transport en latitude
281c      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     .        32,'ave(X)',dt_cum,dt_cum)
287         enddo
288      enddo
289
290
291c      print*,'8HISTDEF'
292               CALL histend(fileid)
293
294
295      endif
296
297
298c=====================================================================
299c   Calcul des champs dynamiques
300c   ----------------------------
301
302c   énergie cinétique
303      CALL covcont(llm,ucov,vcov,ucont,vcont)
304      CALL enercin(vcov,ucov,vcont,ucont,ecin)
305
306c   moment cinétique
307      do l=1,llm
308         ang(:,:,l)=ucov(:,:,l)+constang(:,:)
309         unat(:,:,l)=ucont(:,:,l)*cu(:,:)
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)=q(:,:,:,1)
318      Q(:,:,:,iun)=1.
319
320
321c=====================================================================
322c   Cumul
323c=====================================================================
324c
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     . WRITE(lunout,*)'dans bilan_dyn ',icum,'->',icum+1
337      icum=icum+1
338
339c   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
348c=====================================================================
349c  FLUX ET TENDANCES
350c=====================================================================
351
352c   Flux longitudinal
353c   -----------------
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
366c    flux méridien
367c    -------------
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
380c    tendances
381c    ---------
382
383c   convergence horizontale
384      call  convflu(flux_uQ_cum,flux_vQ_cum,llm*nQ,dQ)
385
386c   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     . WRITE(lunout,*)'Apres les calculs fait a chaque pas'
403c=====================================================================
404c   PAS DE TEMPS D'ECRITURE
405c=====================================================================
406      if (icum.eq.ncum) then
407c=====================================================================
408
409      IF (prt_level > 5)
410     . WRITE(lunout,*)'Pas d ecriture'
411
412c   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
426c   A retravailler eventuellement
427c   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 
432c=====================================================================
433c   Transport méridien
434c=====================================================================
435
436c   cumul zonal des masses des mailles
437c   ----------------------------------
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(1,j)/zmasse(j,l)
448         enddo
449      enddo
450
451c     print*,'3OK'
452c   --------------------------------------------------------------
453c   calcul de la moyenne zonale du transport :
454c   ------------------------------------------
455c
456c                                     --
457c TOT : la circulation totale       [ vq ]
458c
459c                                      -     -
460c MMC : mean meridional circulation [ v ] [ q ]
461c
462c                                     ----      --       - -
463c TRS : transitoires                [ v'q'] = [ vq ] - [ v q ]
464c
465c                                     - * - *       - -       -     -
466c STT : stationaires                [ v   q   ] = [ v q ] - [ v ] [ q ]
467c
468c                                              - -
469c    on utilise aussi l'intermediaire TMP :  [ v q ]
470c
471c    la variable zfactv transforme un transport meridien cumule
472c    en kg/s * unte-du-champ-transporte en m/s * unite-du-champ-transporte
473c
474c   --------------------------------------------------------------
475
476
477c   ----------------------------------------
478c   Transport dans le plan latitude-altitude
479c   ----------------------------------------
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
487c              print*,'j,l,iQ=',j,l,iQ
488c   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
498c              print*,'aOK'
499c   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
508c   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
516c   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
525c     print*,'4OK'
526c   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     s      ,jjm*llm,ndex3d)
532         enddo
533         call histwrite(fileid,'psi'//nom(iQ),itau,psiQ(:,1:llm,iQ)
534     s      ,jjm*llm,ndex3d)
535      enddo
536
537      call histwrite(fileid,'masse',itau,zmasse
538     s   ,jjm*llm,ndex3d)
539      call histwrite(fileid,'v',itau,zv
540     s   ,jjm*llm,ndex3d)
541      psi=psi*1.e-9
542      call histwrite(fileid,'psi',itau,psi(:,1:llm),jjm*llm,ndex3d)
543
544      endif
545
546
547c   -----------------
548c   Moyenne verticale
549c   -----------------
550
551      zamasse=0.
552      do l=1,llm
553         zamasse(:)=zamasse(:)+zmasse(:,l)
554      enddo
555      zavQ=0.
556      do iQ=1,nQ
557         do itr=2,ntr
558            do l=1,llm
559               zavQ(:,itr,iQ)=zavQ(:,itr,iQ)+zvQ(:,l,itr,iQ)*zmasse(:,l)
560            enddo
561            zavQ(:,itr,iQ)=zavQ(:,itr,iQ)/zamasse(:)
562            call histwrite(fileid,'a'//znom(itr,iQ),itau,zavQ(:,itr,iQ)
563     s      ,jjm*llm,ndex3d)
564         enddo
565      enddo
566
567c     on doit pouvoir tracer systematiquement la fonction de courant.
568
569c=====================================================================
570c/////////////////////////////////////////////////////////////////////
571      icum=0                  !///////////////////////////////////////
572      endif ! icum.eq.ncum    !///////////////////////////////////////
573c/////////////////////////////////////////////////////////////////////
574c=====================================================================
575
576      return
577      end
Note: See TracBrowser for help on using the repository browser.