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

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

Put comgeom.h, comgeom2.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.7 KB
Line 
1! $Id: bilan_dyn.F90 5136 2024-07-28 14:17:54Z 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
18  IMPLICIT NONE
19
20  INCLUDE "dimensions.h"
21  INCLUDE "paramet.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    icum = 0
173    ! initialisation des fichiers
174    first = .FALSE.
175    !   ncum est la frequence de stokage en pas de temps
176    ncum = dt_cum / dt_app
177    IF (abs(ncum * dt_app - dt_cum)>1.e-5 * dt_app) THEN
178      WRITE(lunout, *) &
179              'Pb : le pas de cumule doit etre multiple du pas'
180      WRITE(lunout, *)'dt_app=', dt_app
181      WRITE(lunout, *)'dt_cum=', dt_cum
182      CALL abort_gcm('bilan_dyn', 'stopped', 1)
183    endif
184
185    IF (i_sortie==1) THEN
186      file = 'dynzon'
187      CALL inigrads(ifile, 1 &
188              , 0., 180. / pi, 0., 0., jjm, rlatv, -90., 90., 180. / pi &
189              , llm, presnivs, 1. &
190              , dt_cum, file, 'dyn_zon ')
191    endif
192
193    nom(itemp) = 'T'
194    nom(igeop) = 'gz'
195    nom(iecin) = 'K'
196    nom(iang) = 'ang'
197    nom(iu) = 'u'
198    nom(iovap) = 'ovap'
199    nom(iun) = 'un'
200
201    unites(itemp) = 'K'
202    unites(igeop) = 'm2/s2'
203    unites(iecin) = 'm2/s2'
204    unites(iang) = 'ang'
205    unites(iu) = 'm/s'
206    unites(iovap) = 'kg/kg'
207    unites(iun) = 'un'
208
209
210    !   Initialisation du fichier contenant les moyennes zonales.
211    !   ---------------------------------------------------------
212
213    infile = 'dynzon'
214
215    zan = annee_ref
216    dayref = day_ref
217    CALL ymds2ju(zan, 1, dayref, 0.0, zjulian)
218    tau0 = itau_dyn
219
220    rlong = 0.
221    rlatg = rlatv * 180. / pi
222
223    CALL histbeg(infile, 1, rlong, jjm, rlatg, &
224            1, 1, 1, jjm, &
225            tau0, zjulian, dt_cum, thoriid, fileid)
226
227    !
228    !  Appel a histvert pour la grille verticale
229    !
230    CALL histvert(fileid, 'presnivs', 'Niveaux sigma', 'mb', &
231            llm, presnivs, zvertiid)
232    !
233    !  Appels a histdef pour la definition des variables a sauvegarder
234    do iQ = 1, nQ
235      do itr = 1, ntr
236        IF(itr==1) THEN
237          znom(itr, iQ) = nom(iQ)
238          znoml(itr, iQ) = nom(iQ)
239          zunites(itr, iQ) = unites(iQ)
240        else
241          znom(itr, iQ) = ctrs(itr) // 'v' // nom(iQ)
242          znoml(itr, iQ) = 'transport : v * ' // nom(iQ) // ' ' // ctrs(itr)
243          zunites(itr, iQ) = 'm/s * ' // unites(iQ)
244        endif
245      enddo
246    enddo
247
248    !   Declarations des champs avec dimension verticale
249    ! PRINT*,'1HISTDEF'
250    do iQ = 1, nQ
251      do itr = 1, ntr
252        IF (prt_level > 5) &
253                WRITE(lunout, *)'var ', itr, iQ &
254                        , znom(itr, iQ), znoml(itr, iQ), zunites(itr, iQ)
255        CALL histdef(fileid, znom(itr, iQ), znoml(itr, iQ), &
256                zunites(itr, iQ), 1, jjm, thoriid, llm, 1, llm, zvertiid, &
257                32, 'ave(X)', dt_cum, dt_cum)
258      enddo
259      !   Declarations pour les fonctions de courant
260      ! PRINT*,'2HISTDEF'
261      CALL histdef(fileid, 'psi' // nom(iQ) &
262              , 'stream fn. ' // znoml(itot, iQ), &
263              zunites(itot, iQ), 1, jjm, thoriid, llm, 1, llm, zvertiid, &
264              32, 'ave(X)', dt_cum, dt_cum)
265    enddo
266
267
268    !   Declarations pour les champs de transport d'air
269    ! PRINT*,'3HISTDEF'
270    CALL histdef(fileid, 'masse', 'masse', &
271            'kg', 1, jjm, thoriid, llm, 1, llm, zvertiid, &
272            32, 'ave(X)', dt_cum, dt_cum)
273    CALL histdef(fileid, 'v', 'v', &
274            'm/s', 1, jjm, thoriid, llm, 1, llm, zvertiid, &
275            32, 'ave(X)', dt_cum, dt_cum)
276    !   Declarations pour les fonctions de courant
277    ! PRINT*,'4HISTDEF'
278    CALL histdef(fileid, 'psi', 'stream fn. MMC ', 'mega t/s', &
279            1, jjm, thoriid, llm, 1, llm, zvertiid, &
280            32, 'ave(X)', dt_cum, dt_cum)
281
282
283    !   Declaration des champs 1D de transport en latitude
284    ! PRINT*,'5HISTDEF'
285    do iQ = 1, nQ
286      do itr = 2, ntr
287        CALL histdef(fileid, 'a' // znom(itr, iQ), znoml(itr, iQ), &
288                zunites(itr, iQ), 1, jjm, thoriid, 1, 1, 1, -99, &
289                32, 'ave(X)', dt_cum, dt_cum)
290      enddo
291    enddo
292
293
294    ! PRINT*,'8HISTDEF'
295    CALL histend(fileid)
296
297  ENDIF
298
299
300  !=====================================================================
301  !   Calcul des champs dynamiques
302  !   ----------------------------
303
304  !   énergie cinétique
305  ucont(:, :, :) = 0
306  CALL covcont(llm, ucov, vcov, ucont, vcont)
307  CALL enercin(vcov, ucov, vcont, ucont, ecin)
308
309  !   moment cinétique
310  do l = 1, llm
311    ang(:, :, l) = ucov(:, :, l) + constang(:, :)
312    unat(:, :, l) = ucont(:, :, l) * cu(:, :)
313  enddo
314
315  Q(:, :, :, itemp) = teta(:, :, :) * pk(:, :, :) / cpp
316  Q(:, :, :, igeop) = phi(:, :, :)
317  Q(:, :, :, iecin) = ecin(:, :, :)
318  Q(:, :, :, iang) = ang(:, :, :)
319  Q(:, :, :, iu) = unat(:, :, :)
320  Q(:, :, :, iovap) = trac(:, :, :, 1)
321  Q(:, :, :, iun) = 1.
322
323
324  !=====================================================================
325  !   Cumul
326  !=====================================================================
327  !
328  IF(icum==0) THEN
329    ps_cum = 0.
330    masse_cum = 0.
331    flux_u_cum = 0.
332    flux_v_cum = 0.
333    Q_cum = 0.
334    flux_vQ_cum = 0.
335    flux_uQ_cum = 0.
336  ENDIF
337
338  IF (prt_level > 5) &
339          WRITE(lunout, *)'dans bilan_dyn ', icum, '->', icum + 1
340  icum = icum + 1
341
342  !   accumulation des flux de masse horizontaux
343  ps_cum = ps_cum + ps
344  masse_cum = masse_cum + masse
345  flux_u_cum = flux_u_cum + flux_u
346  flux_v_cum = flux_v_cum + flux_v
347  do iQ = 1, nQ
348    Q_cum(:, :, :, iQ) = Q_cum(:, :, :, iQ) + Q(:, :, :, iQ) * masse(:, :, :)
349  enddo
350
351  !=====================================================================
352  !  FLUX ET TENDANCES
353  !=====================================================================
354
355  !   Flux longitudinal
356  !   -----------------
357  do iQ = 1, nQ
358    do l = 1, llm
359      do j = 1, jjp1
360        do i = 1, iim
361          flux_uQ_cum(i, j, l, iQ) = flux_uQ_cum(i, j, l, iQ) &
362                  + flux_u(i, j, l) * 0.5 * (Q(i, j, l, iQ) + Q(i + 1, j, l, iQ))
363        enddo
364        flux_uQ_cum(iip1, j, l, iQ) = flux_uQ_cum(1, j, l, iQ)
365      enddo
366    enddo
367  enddo
368
369  !    flux méridien
370  !    -------------
371  do iQ = 1, nQ
372    do l = 1, llm
373      do j = 1, jjm
374        do i = 1, iip1
375          flux_vQ_cum(i, j, l, iQ) = flux_vQ_cum(i, j, l, iQ) &
376                  + flux_v(i, j, l) * 0.5 * (Q(i, j, l, iQ) + Q(i, j + 1, l, iQ))
377        enddo
378      enddo
379    enddo
380  enddo
381
382
383  !    tendances
384  !    ---------
385
386  !   convergence horizontale
387  CALL  convflu(flux_uQ_cum, flux_vQ_cum, llm * nQ, dQ)
388
389  !   calcul de la vitesse verticale
390  CALL convmas(flux_u_cum, flux_v_cum, convm)
391  CALL vitvert(convm, w)
392
393  do iQ = 1, nQ
394    do l = 1, llm - 1
395      do j = 1, jjp1
396        do i = 1, iip1
397          ww = -0.5 * w(i, j, l + 1) * (Q(i, j, l, iQ) + Q(i, j, l + 1, iQ))
398          dQ(i, j, l, iQ) = dQ(i, j, l, iQ) - ww
399          dQ(i, j, l + 1, iQ) = dQ(i, j, l + 1, iQ) + ww
400        enddo
401      enddo
402    enddo
403  enddo
404  IF (prt_level > 5) &
405          WRITE(lunout, *)'Apres les calculs fait a chaque pas'
406  !=====================================================================
407  !   PAS DE TEMPS D'ECRITURE
408  !=====================================================================
409  IF (icum==ncum) THEN
410    !=====================================================================
411
412    IF (prt_level > 5) &
413            WRITE(lunout, *)'Pas d ecriture'
414
415    !   Normalisation
416    do iQ = 1, nQ
417      Q_cum(:, :, :, iQ) = Q_cum(:, :, :, iQ) / masse_cum(:, :, :)
418    enddo
419    zz = 1. / REAL(ncum)
420    ps_cum = ps_cum * zz
421    masse_cum = masse_cum * zz
422    flux_u_cum = flux_u_cum * zz
423    flux_v_cum = flux_v_cum * zz
424    flux_uQ_cum = flux_uQ_cum * zz
425    flux_vQ_cum = flux_vQ_cum * zz
426    dQ = dQ * zz
427
428
429    !   A retravailler eventuellement
430    !   division de dQ par la masse pour revenir aux bonnes grandeurs
431    do iQ = 1, nQ
432      dQ(:, :, :, iQ) = dQ(:, :, :, iQ) / masse_cum(:, :, :)
433    enddo
434
435    !=====================================================================
436    !   Transport méridien
437    !=====================================================================
438
439    !   cumul zonal des masses des mailles
440    !   ----------------------------------
441    zv = 0.
442    zmasse = 0.
443    CALL massbar(masse_cum, massebx, masseby)
444    do l = 1, llm
445      do j = 1, jjm
446        do i = 1, iim
447          zmasse(j, l) = zmasse(j, l) + masseby(i, j, l)
448          zv(j, l) = zv(j, l) + flux_v_cum(i, j, l)
449        enddo
450        zfactv(j, l) = cv(1, j) / zmasse(j, l)
451      enddo
452    enddo
453
454    ! PRINT*,'3OK'
455    !   --------------------------------------------------------------
456    !   calcul de la moyenne zonale du transport :
457    !   ------------------------------------------
458    !
459    !                                 --
460    ! TOT : la circulation totale       [ vq ]
461    !
462    !                                  -     -
463    ! MMC : mean meridional circulation [ v ] [ q ]
464    !
465    !                                 ----      --       - -
466    ! TRS : transitoires                [ v'q'] = [ vq ] - [ v q ]
467    !
468    !                                 - * - *       - -       -     -
469    ! STT : stationaires                [ v   q   ] = [ v q ] - [ v ] [ q ]
470    !
471    !                                          - -
472    !    on utilise aussi l'intermediaire TMP :  [ v q ]
473    !
474    !    la variable zfactv transforme un transport meridien cumule
475    !    en kg/s * unte-du-champ-transporte en m/s * unite-du-champ-transporte
476    !
477    !   --------------------------------------------------------------
478
479
480    !   ----------------------------------------
481    !   Transport dans le plan latitude-altitude
482    !   ----------------------------------------
483
484    zvQ = 0.
485    psiQ = 0.
486    do iQ = 1, nQ
487      zvQtmp = 0.
488      do l = 1, llm
489        do j = 1, jjm
490          ! PRINT*,'j,l,iQ=',j,l,iQ
491          !   Calcul des moyennes zonales du transort total et de zvQtmp
492          do i = 1, iim
493            zvQ(j, l, itot, iQ) = zvQ(j, l, itot, iQ) &
494                    + flux_vQ_cum(i, j, l, iQ)
495            zqy = 0.5 * (Q_cum(i, j, l, iQ) * masse_cum(i, j, l) + &
496                    Q_cum(i, j + 1, l, iQ) * masse_cum(i, j + 1, l))
497            zvQtmp(j, l) = zvQtmp(j, l) + flux_v_cum(i, j, l) * zqy &
498                    / (0.5 * (masse_cum(i, j, l) + masse_cum(i, j + 1, l)))
499            zvQ(j, l, iave, iQ) = zvQ(j, l, iave, iQ) + zqy
500          enddo
501          ! PRINT*,'aOK'
502          !   Decomposition
503          zvQ(j, l, iave, iQ) = zvQ(j, l, iave, iQ) / zmasse(j, l)
504          zvQ(j, l, itot, iQ) = zvQ(j, l, itot, iQ) * zfactv(j, l)
505          zvQtmp(j, l) = zvQtmp(j, l) * zfactv(j, l)
506          zvQ(j, l, immc, iQ) = zv(j, l) * zvQ(j, l, iave, iQ) * zfactv(j, l)
507          zvQ(j, l, itrs, iQ) = zvQ(j, l, itot, iQ) - zvQtmp(j, l)
508          zvQ(j, l, istn, iQ) = zvQtmp(j, l) - zvQ(j, l, immc, iQ)
509        enddo
510      enddo
511      !   fonction de courant meridienne pour la quantite Q
512      do l = llm, 1, -1
513        do j = 1, jjm
514          psiQ(j, l, iQ) = psiQ(j, l + 1, iQ) + zvQ(j, l, itot, iQ)
515        enddo
516      enddo
517    enddo
518
519    !   fonction de courant pour la circulation meridienne moyenne
520    psi = 0.
521    do l = llm, 1, -1
522      do j = 1, jjm
523        psi(j, l) = psi(j, l + 1) + zv(j, l)
524        zv(j, l) = zv(j, l) * zfactv(j, l)
525      enddo
526    enddo
527
528    ! PRINT*,'4OK'
529    !   sorties proprement dites
530    IF (i_sortie==1) THEN
531      do iQ = 1, nQ
532        do itr = 1, ntr
533          CALL histwrite(fileid, znom(itr, iQ), itau, zvQ(:, :, itr, iQ) &
534                  , jjm * llm, ndex3d)
535        enddo
536        CALL histwrite(fileid, 'psi' // nom(iQ), itau, psiQ(:, 1:llm, iQ) &
537                , jjm * llm, ndex3d)
538      enddo
539
540      CALL histwrite(fileid, 'masse', itau, zmasse &
541              , jjm * llm, ndex3d)
542      CALL histwrite(fileid, 'v', itau, zv &
543              , jjm * llm, ndex3d)
544      psi = psi * 1.e-9
545      CALL histwrite(fileid, 'psi', itau, psi(:, 1:llm), jjm * llm, ndex3d)
546
547    endif
548
549
550    !   -----------------
551    !   Moyenne verticale
552    !   -----------------
553
554    zamasse = 0.
555    do l = 1, llm
556      zamasse(:) = zamasse(:) + zmasse(:, l)
557    enddo
558    zavQ = 0.
559    do iQ = 1, nQ
560      do itr = 2, ntr
561        do l = 1, llm
562          zavQ(:, itr, iQ) = zavQ(:, itr, iQ) + zvQ(:, l, itr, iQ) * zmasse(:, l)
563        enddo
564        zavQ(:, itr, iQ) = zavQ(:, itr, iQ) / zamasse(:)
565        CALL histwrite(fileid, 'a' // znom(itr, iQ), itau, zavQ(:, itr, iQ) &
566                , jjm * llm, ndex3d)
567      enddo
568    enddo
569
570    ! on doit pouvoir tracer systematiquement la fonction de courant.
571
572    !=====================================================================
573    !/////////////////////////////////////////////////////////////////////
574    icum = 0                  !///////////////////////////////////////
575  ENDIF ! icum.EQ.ncum    !///////////////////////////////////////
576  !/////////////////////////////////////////////////////////////////////
577  !=====================================================================
578
579
580END SUBROUTINE  bilan_dyn
Note: See TracBrowser for help on using the repository browser.