source: LMDZ6/branches/Amaury_dev/libf/dyn3d/bilan_dyn.F90 @ 5105

Last change on this file since 5105 was 5105, checked in by abarral, 8 weeks ago

Replace 1DUTILS.h by module lmdz_1dutils.f90
Replace 1DConv.h by module lmdz_old_1dconv.f90 (it's only used by old_* files)
Convert *.F to *.f90
Fix gradsdef.h formatting
Remove unnecessary "RETURN" at the end of functions/subroutines

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