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

Last change on this file since 5258 was 5246, checked in by abarral, 6 weeks ago

Convert fixed-form to free-form sources .F -> .{f,F}90
(WIP: some .F remain, will be handled in subsequent commits)

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