source: LMDZ6/branches/Amaury_dev/libf/dyn3d/lmdz_bilan_dyn.f90 @ 5186

Last change on this file since 5186 was 5186, checked in by abarral, 2 months ago

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