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

Last change on this file since 559 was 559, checked in by lmdzadmin, 20 years ago

Initialisations diverses YM
LF

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