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

Last change on this file since 1032 was 693, checked in by lmdzadmin, 19 years ago

Initialisation ndex3d
IM

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