source: LMDZ6/trunk/libf/dyn3d/bilan_dyn.F90 @ 5267

Last change on this file since 5267 was 5267, checked in by abarral, 2 days ago

Remove CPP_IOIPSL cpp keys uses

  • 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.3 KB
Line 
1!
2! $Id: bilan_dyn.F90 5267 2024-10-23 15:34:46Z 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  IMPLICIT NONE
18
19  include "dimensions.h"
20  include "paramet.h"
21  include "comgeom2.h"
22  include "iniprint.h"
23
24  !====================================================================
25  !
26  !   Sous-programme consacre � des diagnostics dynamiques de base
27  !
28  !
29  !   De facon generale, les moyennes des scalaires Q sont ponderees par
30  !   la masse.
31  !
32  !   Les flux de masse sont eux simplement moyennes.
33  !
34  !====================================================================
35
36  !   Arguments :
37  !   ===========
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
51  !   Local :
52  !   =======
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
62  !ym      character*6 nom(nQ)
63  !ym      character*6 unites(nQ)
64  character*6,save :: nom(nQ)
65  character*6,save :: unites(nQ)
66
67  character(len=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
89  !   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
97  !   champ contenant les scalaires advect�s.
98  real :: Q(iip1,jjp1,llm,nQ)
99
100  !   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
114  !   champs de tansport en moyenne zonale
115  integer :: ntr,itr
116  parameter (ntr=5)
117
118  !ym      character*10 znom(ntr,nQ)
119  !ym      character*20 znoml(ntr,nQ)
120  !ym      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(len=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
139  !   Initialisation du fichier contenant les moyennes zonales.
140  !   ---------------------------------------------------------
141
142  character(len=10) :: infile
143
144  integer :: fileid
145  integer :: thoriid, zvertiid
146  save fileid
147
148  integer :: ndex3d(jjm*llm)
149
150  !   Variables locales
151  !
152  integer :: tau0
153  real :: zjulian
154  character(len=3) :: str
155  character(len=10) :: ctrac
156  integer :: ii,jj
157  integer :: zan, dayref
158  !
159  real :: rlong(jjm),rlatg(jjm)
160
161
162
163  !=====================================================================
164  !   Initialisation
165  !=====================================================================
166
167  time=time+dt_app
168  itau=itau+1
169  !IM
170  ndex3d=0
171
172  if (first) then
173
174
175    icum=0
176    ! initialisation des fichiers
177    first=.false.
178  !   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       call abort_gcm('bilan_dyn','stopped',1)
186    endif
187
188    if (i_sortie.eq.1) then
189     file='dynzon'
190     call inigrads(ifile,1 &
191           ,0.,180./pi,0.,0.,jjm,rlatv,-90.,90.,180./pi &
192           ,llm,presnivs,1. &
193           ,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
213  !   Initialisation du fichier contenant les moyennes zonales.
214  !   ---------------------------------------------------------
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
230  !
231  !  Appel a histvert pour la grille verticale
232  !
233  call histvert(fileid, 'presnivs', 'Niveaux sigma','mb', &
234        llm, presnivs, zvertiid)
235  !
236  !  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
251  !   Declarations des champs avec dimension verticale
252   ! 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
262  !   Declarations pour les fonctions de courant
263   ! 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
271  !   Declarations pour les champs de transport d'air
272   ! 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)
279  !   Declarations pour les fonctions de courant
280   ! 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
286  !   Declaration des champs 1D de transport en latitude
287   ! 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
297   ! print*,'8HISTDEF'
298           CALL histend(fileid)
299
300
301  endif
302
303
304  !=====================================================================
305  !   Calcul des champs dynamiques
306  !   ----------------------------
307
308  !   �nergie cin�tique
309  ucont(:,:,:)=0
310  CALL covcont(llm,ucov,vcov,ucont,vcont)
311  CALL enercin(vcov,ucov,vcont,ucont,ecin)
312
313  !   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
328  !=====================================================================
329  !   Cumul
330  !=====================================================================
331  !
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
346  !   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
355  !=====================================================================
356  !  FLUX ET TENDANCES
357  !=====================================================================
358
359  !   Flux longitudinal
360  !   -----------------
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                    +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
373  !    flux m�ridien
374  !    -------------
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                    +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
387  !    tendances
388  !    ---------
389
390  !   convergence horizontale
391  call  convflu(flux_uQ_cum,flux_vQ_cum,llm*nQ,dQ)
392
393  !   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'
410  !=====================================================================
411  !   PAS DE TEMPS D'ECRITURE
412  !=====================================================================
413  if (icum.eq.ncum) then
414  !=====================================================================
415
416  IF (prt_level > 5) &
417        WRITE(lunout,*)'Pas d ecriture'
418
419  !   Normalisation
420  do iQ=1,nQ
421     Q_cum(:,:,:,iQ)=Q_cum(:,:,:,iQ)/masse_cum(:,:,:)
422  enddo
423  zz=1./REAL(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
433  !   A retravailler eventuellement
434  !   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
439  !=====================================================================
440  !   Transport m�ridien
441  !=====================================================================
442
443  !   cumul zonal des masses des mailles
444  !   ----------------------------------
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
458  ! print*,'3OK'
459  !   --------------------------------------------------------------
460  !   calcul de la moyenne zonale du transport :
461  !   ------------------------------------------
462  !
463  !                                 --
464  ! TOT : la circulation totale       [ vq ]
465  !
466  !                                  -     -
467  ! MMC : mean meridional circulation [ v ] [ q ]
468  !
469  !                                 ----      --       - -
470  ! TRS : transitoires                [ v'q'] = [ vq ] - [ v q ]
471  !
472  !                                 - * - *       - -       -     -
473  ! STT : stationaires                [ v   q   ] = [ v q ] - [ v ] [ q ]
474  !
475  !                                          - -
476  !    on utilise aussi l'intermediaire TMP :  [ v q ]
477  !
478  !    la variable zfactv transforme un transport meridien cumule
479  !    en kg/s * unte-du-champ-transporte en m/s * unite-du-champ-transporte
480  !
481  !   --------------------------------------------------------------
482
483
484  !   ----------------------------------------
485  !   Transport dans le plan latitude-altitude
486  !   ----------------------------------------
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
494           ! print*,'j,l,iQ=',j,l,iQ
495  !   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                    +flux_vQ_cum(i,j,l,iQ)
499              zqy=      0.5*(Q_cum(i,j,l,iQ)*masse_cum(i,j,l)+ &
500                    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                    /(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
505           ! print*,'aOK'
506  !   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
515  !   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
523  !   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
532  ! print*,'4OK'
533  !   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              ,jjm*llm,ndex3d)
539     enddo
540     call histwrite(fileid,'psi'//nom(iQ),itau,psiQ(:,1:llm,iQ) &
541           ,jjm*llm,ndex3d)
542  enddo
543
544  call histwrite(fileid,'masse',itau,zmasse &
545        ,jjm*llm,ndex3d)
546  call histwrite(fileid,'v',itau,zv &
547        ,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
554  !   -----------------
555  !   Moyenne verticale
556  !   -----------------
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              ,jjm*llm,ndex3d)
571     enddo
572  enddo
573
574  ! on doit pouvoir tracer systematiquement la fonction de courant.
575
576  !=====================================================================
577  !/////////////////////////////////////////////////////////////////////
578  icum=0                  !///////////////////////////////////////
579  endif ! icum.eq.ncum    !///////////////////////////////////////
580  !/////////////////////////////////////////////////////////////////////
581  !=====================================================================
582
583  return
584end subroutine bilan_dyn
Note: See TracBrowser for help on using the repository browser.