source: LMDZ4/branches/LMDZ4-dev-20091210/libf/dyn3d/bilan_dyn.F @ 2063

Last change on this file since 2063 was 1186, checked in by Ehouarn Millour, 15 years ago

Cleanup around IOIPSL, so that LMDZ dynamics may be used without IOIPSL.

  • moved ersatz IOIPSL routines (ioipsl_* , taken from IOIPSLv2_1_8, so that 'getin' function may be used even if not using the IOIPSL library) from dyn3d/dyn3dpar to bibio.
  • enclosed 'use ioipsl' instruction with #ifdef CPP_IOIPSL cpp keys.

EM

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 16.1 KB
Line 
1!
2! $Id: bilan_dyn.F 1186 2009-06-18 09:20:44Z 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
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#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
26c====================================================================
27c
28c   Sous-programme consacre à des diagnostics dynamiques de base
29c
30c
31c   De facon generale, les moyennes des scalaires Q sont ponderees par
32c   la masse.
33c
34c   Les flux de masse sont eux simplement moyennes.
35c
36c====================================================================
37
38c   Arguments :
39c   ===========
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
53c   Local :
54c   =======
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
64cym      character*6 nom(nQ)
65cym      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
91c   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
99c   champ contenant les scalaires advectés.
100      real Q(iip1,jjp1,llm,nQ)
101   
102c   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
116c   champs de tansport en moyenne zonale
117      integer ntr,itr
118      parameter (ntr=5)
119
120cym      character*10 znom(ntr,nQ)
121cym      character*20 znoml(ntr,nQ)
122cym      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
141c   Initialisation du fichier contenant les moyennes zonales.
142c   ---------------------------------------------------------
143
144      character*10 infile
145
146      integer fileid
147      integer thoriid, zvertiid
148      save fileid
149
150      integer ndex3d(jjm*llm)
151
152C   Variables locales
153C
154      integer tau0
155      real zjulian
156      character*3 str
157      character*10 ctrac
158      integer ii,jj
159      integer zan, dayref
160C
161      real rlong(jjm),rlatg(jjm)
162
163
164
165c=====================================================================
166c   Initialisation
167c=====================================================================
168
169      time=time+dt_app
170      itau=itau+1
171cIM
172      ndex3d=0
173
174      if (first) then
175
176
177        icum=0
178c       initialisation des fichiers
179        first=.false.
180c   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
215c   Initialisation du fichier contenant les moyennes zonales.
216c   ---------------------------------------------------------
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
232C
233C  Appel a histvert pour la grille verticale
234C
235      call histvert(fileid, 'presnivs', 'Niveaux sigma','mb',
236     .              llm, presnivs, zvertiid)
237C
238C  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
253c   Declarations des champs avec dimension verticale
254c      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
264c   Declarations pour les fonctions de courant
265c      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
273c   Declarations pour les champs de transport d'air
274c      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)
281c   Declarations pour les fonctions de courant
282c      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
288c   Declaration des champs 1D de transport en latitude
289c      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
299c      print*,'8HISTDEF'
300               CALL histend(fileid)
301
302
303      endif
304
305
306c=====================================================================
307c   Calcul des champs dynamiques
308c   ----------------------------
309
310c   énergie cinétique
311      ucont(:,:,:)=0
312      CALL covcont(llm,ucov,vcov,ucont,vcont)
313      CALL enercin(vcov,ucov,vcont,ucont,ecin)
314
315c   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
330c=====================================================================
331c   Cumul
332c=====================================================================
333c
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
348c   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
357c=====================================================================
358c  FLUX ET TENDANCES
359c=====================================================================
360
361c   Flux longitudinal
362c   -----------------
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
375c    flux méridien
376c    -------------
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
389c    tendances
390c    ---------
391
392c   convergence horizontale
393      call  convflu(flux_uQ_cum,flux_vQ_cum,llm*nQ,dQ)
394
395c   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'
412c=====================================================================
413c   PAS DE TEMPS D'ECRITURE
414c=====================================================================
415      if (icum.eq.ncum) then
416c=====================================================================
417
418      IF (prt_level > 5)
419     . WRITE(lunout,*)'Pas d ecriture'
420
421c   Normalisation
422      do iQ=1,nQ
423         Q_cum(:,:,:,iQ)=Q_cum(:,:,:,iQ)/masse_cum(:,:,:)
424      enddo
425      zz=1./float(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
435c   A retravailler eventuellement
436c   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 
441c=====================================================================
442c   Transport méridien
443c=====================================================================
444
445c   cumul zonal des masses des mailles
446c   ----------------------------------
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
460c     print*,'3OK'
461c   --------------------------------------------------------------
462c   calcul de la moyenne zonale du transport :
463c   ------------------------------------------
464c
465c                                     --
466c TOT : la circulation totale       [ vq ]
467c
468c                                      -     -
469c MMC : mean meridional circulation [ v ] [ q ]
470c
471c                                     ----      --       - -
472c TRS : transitoires                [ v'q'] = [ vq ] - [ v q ]
473c
474c                                     - * - *       - -       -     -
475c STT : stationaires                [ v   q   ] = [ v q ] - [ v ] [ q ]
476c
477c                                              - -
478c    on utilise aussi l'intermediaire TMP :  [ v q ]
479c
480c    la variable zfactv transforme un transport meridien cumule
481c    en kg/s * unte-du-champ-transporte en m/s * unite-du-champ-transporte
482c
483c   --------------------------------------------------------------
484
485
486c   ----------------------------------------
487c   Transport dans le plan latitude-altitude
488c   ----------------------------------------
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
496c              print*,'j,l,iQ=',j,l,iQ
497c   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
507c              print*,'aOK'
508c   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
517c   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
525c   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
534c     print*,'4OK'
535c   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
556c   -----------------
557c   Moyenne verticale
558c   -----------------
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
576c     on doit pouvoir tracer systematiquement la fonction de courant.
577
578c=====================================================================
579c/////////////////////////////////////////////////////////////////////
580      icum=0                  !///////////////////////////////////////
581      endif ! icum.eq.ncum    !///////////////////////////////////////
582c/////////////////////////////////////////////////////////////////////
583c=====================================================================
584
585      return
586      end
Note: See TracBrowser for help on using the repository browser.