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

Last change on this file since 5159 was 5159, checked in by abarral, 7 weeks ago

Put dimensions.h and paramet.h into modules

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