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

Last change on this file since 706 was 630, checked in by Laurent Fairhead, 20 years ago

Import d'une version parallele de la dynamique YM
LF

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