source: LMDZ4/trunk/libf/dyn3dpar/bilan_dyn_p.F @ 885

Last change on this file since 885 was 764, checked in by Laurent Fairhead, 17 years ago

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