source: LMDZ6/trunk/libf/dyn3d/bilan_dyn.f90

Last change on this file was 5272, checked in by abarral, 31 hours ago

Turn paramet.h into a module

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 15.5 KB
Line 
1!
2! $Id: bilan_dyn.f90 5272 2024-10-24 15:53:15Z abarral $
3!
4SUBROUTINE bilan_dyn (ntrac,dt_app,dt_cum, &
5        ps,masse,pk,flux_u,flux_v,teta,phi,ucov,vcov,trac)
6
7  !   AFAIRE
8  !   Prevoir en champ nq+1 le diagnostique de l'energie
9  !   en faisant Qzon=Cv T + L * ...
10  !             vQ..A=Cp T + L * ...
11
12  USE IOIPSL
13  USE comconst_mod, ONLY: pi, cpp
14  USE comvert_mod, ONLY: presnivs
15  USE temps_mod, ONLY: annee_ref, day_ref, itau_dyn
16
17  USE dimensions_mod, ONLY: iim, jjm, llm, ndm
18USE paramet_mod_h, ONLY: iip1, iip2, iip3, jjp1, llmp1, llmp2, llmm1, kftd, ip1jm, ip1jmp1, &
19          ip1jmi1, ijp1llm, ijmllm, mvar, jcfil, jcfllm
20IMPLICIT NONE
21
22
23
24  include "comgeom2.h"
25  include "iniprint.h"
26
27  !====================================================================
28  !
29  !   Sous-programme consacre � des diagnostics dynamiques de base
30  !
31  !
32  !   De facon generale, les moyennes des scalaires Q sont ponderees par
33  !   la masse.
34  !
35  !   Les flux de masse sont eux simplement moyennes.
36  !
37  !====================================================================
38
39  !   Arguments :
40  !   ===========
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
54  !   Local :
55  !   =======
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
65  !ym      character*6 nom(nQ)
66  !ym      character*6 unites(nQ)
67  character*6,save :: nom(nQ)
68  character*6,save :: unites(nQ)
69
70  character(len=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
92  !   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
100  !   champ contenant les scalaires advect�s.
101  real :: Q(iip1,jjp1,llm,nQ)
102
103  !   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
117  !   champs de tansport en moyenne zonale
118  integer :: ntr,itr
119  parameter (ntr=5)
120
121  !ym      character*10 znom(ntr,nQ)
122  !ym      character*20 znoml(ntr,nQ)
123  !ym      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(len=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
142  !   Initialisation du fichier contenant les moyennes zonales.
143  !   ---------------------------------------------------------
144
145  character(len=10) :: infile
146
147  integer :: fileid
148  integer :: thoriid, zvertiid
149  save fileid
150
151  integer :: ndex3d(jjm*llm)
152
153  !   Variables locales
154  !
155  integer :: tau0
156  real :: zjulian
157  character(len=3) :: str
158  character(len=10) :: ctrac
159  integer :: ii,jj
160  integer :: zan, dayref
161  !
162  real :: rlong(jjm),rlatg(jjm)
163
164
165
166  !=====================================================================
167  !   Initialisation
168  !=====================================================================
169
170  time=time+dt_app
171  itau=itau+1
172  !IM
173  ndex3d=0
174
175  if (first) then
176
177
178    icum=0
179    ! initialisation des fichiers
180    first=.false.
181  !   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       call abort_gcm('bilan_dyn','stopped',1)
189    endif
190
191    if (i_sortie.eq.1) then
192     file='dynzon'
193     call inigrads(ifile,1 &
194           ,0.,180./pi,0.,0.,jjm,rlatv,-90.,90.,180./pi &
195           ,llm,presnivs,1. &
196           ,dt_cum,file,'dyn_zon ')
197    endif
198
199    nom(itemp)='T'
200    nom(igeop)='gz'
201    nom(iecin)='K'
202    nom(iang)='ang'
203    nom(iu)='u'
204    nom(iovap)='ovap'
205    nom(iun)='un'
206
207    unites(itemp)='K'
208    unites(igeop)='m2/s2'
209    unites(iecin)='m2/s2'
210    unites(iang)='ang'
211    unites(iu)='m/s'
212    unites(iovap)='kg/kg'
213    unites(iun)='un'
214
215
216  !   Initialisation du fichier contenant les moyennes zonales.
217  !   ---------------------------------------------------------
218
219  infile='dynzon'
220
221  zan = annee_ref
222  dayref = day_ref
223  CALL ymds2ju(zan, 1, dayref, 0.0, zjulian)
224  tau0 = itau_dyn
225
226  rlong=0.
227  rlatg=rlatv*180./pi
228
229  call histbeg(infile, 1, rlong, jjm, rlatg, &
230        1, 1, 1, jjm, &
231        tau0, zjulian, dt_cum, thoriid, fileid)
232
233  !
234  !  Appel a histvert pour la grille verticale
235  !
236  call histvert(fileid, 'presnivs', 'Niveaux sigma','mb', &
237        llm, presnivs, zvertiid)
238  !
239  !  Appels a histdef pour la definition des variables a sauvegarder
240  do iQ=1,nQ
241     do itr=1,ntr
242        if(itr.eq.1) then
243           znom(itr,iQ)=nom(iQ)
244           znoml(itr,iQ)=nom(iQ)
245           zunites(itr,iQ)=unites(iQ)
246        else
247           znom(itr,iQ)=ctrs(itr)//'v'//nom(iQ)
248           znoml(itr,iQ)='transport : v * '//nom(iQ)//' '//ctrs(itr)
249           zunites(itr,iQ)='m/s * '//unites(iQ)
250        endif
251     enddo
252  enddo
253
254  !   Declarations des champs avec dimension verticale
255   ! print*,'1HISTDEF'
256  do iQ=1,nQ
257     do itr=1,ntr
258  IF (prt_level > 5) &
259        WRITE(lunout,*)'var ',itr,iQ &
260        ,znom(itr,iQ),znoml(itr,iQ),zunites(itr,iQ)
261        call histdef(fileid,znom(itr,iQ),znoml(itr,iQ), &
262              zunites(itr,iQ),1,jjm,thoriid,llm,1,llm,zvertiid, &
263              32,'ave(X)',dt_cum,dt_cum)
264     enddo
265  !   Declarations pour les fonctions de courant
266   ! print*,'2HISTDEF'
267      call histdef(fileid,'psi'//nom(iQ) &
268            ,'stream fn. '//znoml(itot,iQ), &
269            zunites(itot,iQ),1,jjm,thoriid,llm,1,llm,zvertiid, &
270            32,'ave(X)',dt_cum,dt_cum)
271  enddo
272
273
274  !   Declarations pour les champs de transport d'air
275   ! print*,'3HISTDEF'
276  call histdef(fileid, 'masse', 'masse', &
277        'kg', 1, jjm, thoriid, llm, 1, llm, zvertiid, &
278        32, 'ave(X)', dt_cum, dt_cum)
279  call histdef(fileid, 'v', 'v', &
280        'm/s', 1, jjm, thoriid, llm, 1, llm, zvertiid, &
281        32, 'ave(X)', dt_cum, dt_cum)
282  !   Declarations pour les fonctions de courant
283   ! print*,'4HISTDEF'
284      call histdef(fileid,'psi','stream fn. MMC ','mega t/s', &
285            1,jjm,thoriid,llm,1,llm,zvertiid, &
286            32,'ave(X)',dt_cum,dt_cum)
287
288
289  !   Declaration des champs 1D de transport en latitude
290   ! print*,'5HISTDEF'
291  do iQ=1,nQ
292     do itr=2,ntr
293        call histdef(fileid,'a'//znom(itr,iQ),znoml(itr,iQ), &
294              zunites(itr,iQ),1,jjm,thoriid,1,1,1,-99, &
295              32,'ave(X)',dt_cum,dt_cum)
296     enddo
297  enddo
298
299
300   ! print*,'8HISTDEF'
301           CALL histend(fileid)
302
303
304  endif
305
306
307  !=====================================================================
308  !   Calcul des champs dynamiques
309  !   ----------------------------
310
311  !   �nergie cin�tique
312  ucont(:,:,:)=0
313  CALL covcont(llm,ucov,vcov,ucont,vcont)
314  CALL enercin(vcov,ucov,vcont,ucont,ecin)
315
316  !   moment cin�tique
317  do l=1,llm
318     ang(:,:,l)=ucov(:,:,l)+constang(:,:)
319     unat(:,:,l)=ucont(:,:,l)*cu(:,:)
320  enddo
321
322  Q(:,:,:,itemp)=teta(:,:,:)*pk(:,:,:)/cpp
323  Q(:,:,:,igeop)=phi(:,:,:)
324  Q(:,:,:,iecin)=ecin(:,:,:)
325  Q(:,:,:,iang)=ang(:,:,:)
326  Q(:,:,:,iu)=unat(:,:,:)
327  Q(:,:,:,iovap)=trac(:,:,:,1)
328  Q(:,:,:,iun)=1.
329
330
331  !=====================================================================
332  !   Cumul
333  !=====================================================================
334  !
335  if(icum.EQ.0) then
336     ps_cum=0.
337     masse_cum=0.
338     flux_u_cum=0.
339     flux_v_cum=0.
340     Q_cum=0.
341     flux_vQ_cum=0.
342     flux_uQ_cum=0.
343  endif
344
345  IF (prt_level > 5) &
346        WRITE(lunout,*)'dans bilan_dyn ',icum,'->',icum+1
347  icum=icum+1
348
349  !   accumulation des flux de masse horizontaux
350  ps_cum=ps_cum+ps
351  masse_cum=masse_cum+masse
352  flux_u_cum=flux_u_cum+flux_u
353  flux_v_cum=flux_v_cum+flux_v
354  do iQ=1,nQ
355  Q_cum(:,:,:,iQ)=Q_cum(:,:,:,iQ)+Q(:,:,:,iQ)*masse(:,:,:)
356  enddo
357
358  !=====================================================================
359  !  FLUX ET TENDANCES
360  !=====================================================================
361
362  !   Flux longitudinal
363  !   -----------------
364  do iQ=1,nQ
365     do l=1,llm
366        do j=1,jjp1
367           do i=1,iim
368              flux_uQ_cum(i,j,l,iQ)=flux_uQ_cum(i,j,l,iQ) &
369                    +flux_u(i,j,l)*0.5*(Q(i,j,l,iQ)+Q(i+1,j,l,iQ))
370           enddo
371           flux_uQ_cum(iip1,j,l,iQ)=flux_uQ_cum(1,j,l,iQ)
372        enddo
373     enddo
374  enddo
375
376  !    flux m�ridien
377  !    -------------
378  do iQ=1,nQ
379     do l=1,llm
380        do j=1,jjm
381           do i=1,iip1
382              flux_vQ_cum(i,j,l,iQ)=flux_vQ_cum(i,j,l,iQ) &
383                    +flux_v(i,j,l)*0.5*(Q(i,j,l,iQ)+Q(i,j+1,l,iQ))
384           enddo
385        enddo
386     enddo
387  enddo
388
389
390  !    tendances
391  !    ---------
392
393  !   convergence horizontale
394  call  convflu(flux_uQ_cum,flux_vQ_cum,llm*nQ,dQ)
395
396  !   calcul de la vitesse verticale
397  call convmas(flux_u_cum,flux_v_cum,convm)
398  CALL vitvert(convm,w)
399
400  do iQ=1,nQ
401     do l=1,llm-1
402        do j=1,jjp1
403           do i=1,iip1
404              ww=-0.5*w(i,j,l+1)*(Q(i,j,l,iQ)+Q(i,j,l+1,iQ))
405              dQ(i,j,l  ,iQ)=dQ(i,j,l  ,iQ)-ww
406              dQ(i,j,l+1,iQ)=dQ(i,j,l+1,iQ)+ww
407           enddo
408        enddo
409     enddo
410  enddo
411  IF (prt_level > 5) &
412        WRITE(lunout,*)'Apres les calculs fait a chaque pas'
413  !=====================================================================
414  !   PAS DE TEMPS D'ECRITURE
415  !=====================================================================
416  if (icum.eq.ncum) then
417  !=====================================================================
418
419  IF (prt_level > 5) &
420        WRITE(lunout,*)'Pas d ecriture'
421
422  !   Normalisation
423  do iQ=1,nQ
424     Q_cum(:,:,:,iQ)=Q_cum(:,:,:,iQ)/masse_cum(:,:,:)
425  enddo
426  zz=1./REAL(ncum)
427  ps_cum=ps_cum*zz
428  masse_cum=masse_cum*zz
429  flux_u_cum=flux_u_cum*zz
430  flux_v_cum=flux_v_cum*zz
431  flux_uQ_cum=flux_uQ_cum*zz
432  flux_vQ_cum=flux_vQ_cum*zz
433  dQ=dQ*zz
434
435
436  !   A retravailler eventuellement
437  !   division de dQ par la masse pour revenir aux bonnes grandeurs
438  do iQ=1,nQ
439     dQ(:,:,:,iQ)=dQ(:,:,:,iQ)/masse_cum(:,:,:)
440  enddo
441
442  !=====================================================================
443  !   Transport m�ridien
444  !=====================================================================
445
446  !   cumul zonal des masses des mailles
447  !   ----------------------------------
448  zv=0.
449  zmasse=0.
450  call massbar(masse_cum,massebx,masseby)
451  do l=1,llm
452     do j=1,jjm
453        do i=1,iim
454           zmasse(j,l)=zmasse(j,l)+masseby(i,j,l)
455           zv(j,l)=zv(j,l)+flux_v_cum(i,j,l)
456        enddo
457        zfactv(j,l)=cv(1,j)/zmasse(j,l)
458     enddo
459  enddo
460
461  ! print*,'3OK'
462  !   --------------------------------------------------------------
463  !   calcul de la moyenne zonale du transport :
464  !   ------------------------------------------
465  !
466  !                                 --
467  ! TOT : la circulation totale       [ vq ]
468  !
469  !                                  -     -
470  ! MMC : mean meridional circulation [ v ] [ q ]
471  !
472  !                                 ----      --       - -
473  ! TRS : transitoires                [ v'q'] = [ vq ] - [ v q ]
474  !
475  !                                 - * - *       - -       -     -
476  ! STT : stationaires                [ v   q   ] = [ v q ] - [ v ] [ q ]
477  !
478  !                                          - -
479  !    on utilise aussi l'intermediaire TMP :  [ v q ]
480  !
481  !    la variable zfactv transforme un transport meridien cumule
482  !    en kg/s * unte-du-champ-transporte en m/s * unite-du-champ-transporte
483  !
484  !   --------------------------------------------------------------
485
486
487  !   ----------------------------------------
488  !   Transport dans le plan latitude-altitude
489  !   ----------------------------------------
490
491  zvQ=0.
492  psiQ=0.
493  do iQ=1,nQ
494     zvQtmp=0.
495     do l=1,llm
496        do j=1,jjm
497           ! print*,'j,l,iQ=',j,l,iQ
498  !   Calcul des moyennes zonales du transort total et de zvQtmp
499           do i=1,iim
500              zvQ(j,l,itot,iQ)=zvQ(j,l,itot,iQ) &
501                    +flux_vQ_cum(i,j,l,iQ)
502              zqy=      0.5*(Q_cum(i,j,l,iQ)*masse_cum(i,j,l)+ &
503                    Q_cum(i,j+1,l,iQ)*masse_cum(i,j+1,l))
504              zvQtmp(j,l)=zvQtmp(j,l)+flux_v_cum(i,j,l)*zqy &
505                    /(0.5*(masse_cum(i,j,l)+masse_cum(i,j+1,l)))
506              zvQ(j,l,iave,iQ)=zvQ(j,l,iave,iQ)+zqy
507           enddo
508           ! print*,'aOK'
509  !   Decomposition
510           zvQ(j,l,iave,iQ)=zvQ(j,l,iave,iQ)/zmasse(j,l)
511           zvQ(j,l,itot,iQ)=zvQ(j,l,itot,iQ)*zfactv(j,l)
512           zvQtmp(j,l)=zvQtmp(j,l)*zfactv(j,l)
513           zvQ(j,l,immc,iQ)=zv(j,l)*zvQ(j,l,iave,iQ)*zfactv(j,l)
514           zvQ(j,l,itrs,iQ)=zvQ(j,l,itot,iQ)-zvQtmp(j,l)
515           zvQ(j,l,istn,iQ)=zvQtmp(j,l)-zvQ(j,l,immc,iQ)
516        enddo
517     enddo
518  !   fonction de courant meridienne pour la quantite Q
519     do l=llm,1,-1
520        do j=1,jjm
521           psiQ(j,l,iQ)=psiQ(j,l+1,iQ)+zvQ(j,l,itot,iQ)
522        enddo
523     enddo
524  enddo
525
526  !   fonction de courant pour la circulation meridienne moyenne
527  psi=0.
528  do l=llm,1,-1
529     do j=1,jjm
530        psi(j,l)=psi(j,l+1)+zv(j,l)
531        zv(j,l)=zv(j,l)*zfactv(j,l)
532     enddo
533  enddo
534
535  ! print*,'4OK'
536  !   sorties proprement dites
537  if (i_sortie.eq.1) then
538  do iQ=1,nQ
539     do itr=1,ntr
540        call histwrite(fileid,znom(itr,iQ),itau,zvQ(:,:,itr,iQ) &
541              ,jjm*llm,ndex3d)
542     enddo
543     call histwrite(fileid,'psi'//nom(iQ),itau,psiQ(:,1:llm,iQ) &
544           ,jjm*llm,ndex3d)
545  enddo
546
547  call histwrite(fileid,'masse',itau,zmasse &
548        ,jjm*llm,ndex3d)
549  call histwrite(fileid,'v',itau,zv &
550        ,jjm*llm,ndex3d)
551  psi=psi*1.e-9
552  call histwrite(fileid,'psi',itau,psi(:,1:llm),jjm*llm,ndex3d)
553
554  endif
555
556
557  !   -----------------
558  !   Moyenne verticale
559  !   -----------------
560
561  zamasse=0.
562  do l=1,llm
563     zamasse(:)=zamasse(:)+zmasse(:,l)
564  enddo
565  zavQ=0.
566  do iQ=1,nQ
567     do itr=2,ntr
568        do l=1,llm
569           zavQ(:,itr,iQ)=zavQ(:,itr,iQ)+zvQ(:,l,itr,iQ)*zmasse(:,l)
570        enddo
571        zavQ(:,itr,iQ)=zavQ(:,itr,iQ)/zamasse(:)
572        call histwrite(fileid,'a'//znom(itr,iQ),itau,zavQ(:,itr,iQ) &
573              ,jjm*llm,ndex3d)
574     enddo
575  enddo
576
577  ! on doit pouvoir tracer systematiquement la fonction de courant.
578
579  !=====================================================================
580  !/////////////////////////////////////////////////////////////////////
581  icum=0                  !///////////////////////////////////////
582  endif ! icum.eq.ncum    !///////////////////////////////////////
583  !/////////////////////////////////////////////////////////////////////
584  !=====================================================================
585
586  return
587end subroutine bilan_dyn
Note: See TracBrowser for help on using the repository browser.