source: LMDZ.3.3/branches/rel-LF/libf/dyn3d/bilan_dyn.F @ 1125

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

Repere par GG
LF

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