source: LMDZ.3.3/trunk/libf/dyn3d/bilan_dyn.F @ 3730

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

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