Sources/dyn3d/bilan_dyn.f

module bilan_dyn_m

  IMPLICIT NONE

contains

  SUBROUTINE bilan_dyn(ps, masse, pk, flux_u, flux_v, teta, phi, ucov, vcov, &
       trac, dt_app, dt_cum)

    ! From LMDZ4/libf/dyn3d/bilan_dyn.F, version 1.5 2005/03/16
    ! 10:12:17 fairhead

    ! Sous-programme consacré à des diagnostics dynamiques de base.
    ! De façon générale, les moyennes des scalaires Q sont pondérées
    ! par la masse. Les flux de masse sont, eux, simplement moyennés.

    USE histcom, ONLY: histbeg_totreg, histdef, histend, histvert
    USE calendar, ONLY: ymds2ju
    USE histwrite_m, ONLY: histwrite
    USE dimens_m, ONLY: iim, jjm, llm
    USE paramet_m, ONLY: iip1, jjp1
    USE comconst, ONLY: cpp
    USE comvert, ONLY: presnivs
    USE comgeom, ONLY: constang_2d, cu_2d, cv_2d, rlatv
    USE temps, ONLY: annee_ref, day_ref, itau_dyn
    USE inigrads_m, ONLY: inigrads
    USE nr_util, ONLY: pi

    ! Arguments:

    real, intent(in):: dt_app, dt_cum
    real ps(iip1, jjp1)
    real masse(iip1, jjp1, llm), pk(iip1, jjp1, llm)
    real flux_u(iip1, jjp1, llm)
    real flux_v(iip1, jjm, llm)
    real, intent(in):: teta(iip1, jjp1, llm)
    real phi(iip1, jjp1, llm)
    real ucov(iip1, jjp1, llm)
    real vcov(iip1, jjm, llm)
    real, intent(in):: trac(:, :, :) ! (iim + 1, jjm + 1, llm)

    ! Local:

    integer:: icum  = 0
    integer, save:: ncum
    logical:: first = .true.
    real zz, zqy, zfactv(jjm, llm)

    integer, parameter:: nQ = 7
    character(len=4), parameter:: nom(nQ) = (/'T   ', 'gz  ', 'K   ', 'ang ', &
         'u   ', 'ovap', 'un  '/)
    character(len=5), parameter:: unites(nQ) = (/'K    ', 'm2/s2', 'm2/s2', &
         'ang  ', 'm/s  ', 'kg/kg', 'un   '/)

    real:: time = 0.
    integer:: itau = 0
    real ww

    ! Variables dynamiques intermédiaires
    REAL vcont(iip1, jjm, llm), ucont(iip1, jjp1, llm)
    REAL ang(iip1, jjp1, llm), unat(iip1, jjp1, llm)
    REAL massebx(iip1, jjp1, llm), masseby(iip1, jjm, llm)
    REAL w(iip1, jjp1, llm), ecin(iip1, jjp1, llm), convm(iip1, jjp1, llm)

    ! Champ contenant les scalaires advectés
    real Q(iip1, jjp1, llm, nQ)

    ! Champs cumulés
    real, save:: ps_cum(iip1, jjp1)
    real, save:: masse_cum(iip1, jjp1, llm)
    real, save:: flux_u_cum(iip1, jjp1, llm)
    real, save:: flux_v_cum(iip1, jjm, llm)
    real, save:: Q_cum(iip1, jjp1, llm, nQ)
    real, save:: flux_uQ_cum(iip1, jjp1, llm, nQ)
    real, save:: flux_vQ_cum(iip1, jjm, llm, nQ)
    real dQ(iip1, jjp1, llm, nQ)

    ! champs de tansport en moyenne zonale
    integer itr
    integer, parameter:: ntr = 5

    character(len=10), save:: znom(ntr, nQ)
    character(len=26), save:: znoml(ntr, nQ)
    character(len=12), save:: zunites(ntr, nQ)

    integer, parameter:: iave = 1, itot = 2, immc = 3, itrs = 4, istn = 5
    character(len=3), parameter:: ctrs(ntr) = (/'   ', 'TOT', 'MMC', 'TRS', &
         'STN'/)

    real zvQ(jjm, llm, ntr, nQ), zvQtmp(jjm, llm)
    real zavQ(jjm, 2: ntr, nQ), psiQ(jjm, llm + 1, nQ)
    real zmasse(jjm, llm)

    real zv(jjm, llm), psi(jjm, llm + 1)

    integer i, j, l, iQ

    ! Initialisation du fichier contenant les moyennes zonales.

    integer, save:: fileid
    integer thoriid, zvertiid

    real zjulian
    integer zan, dayref

    real rlong(jjm), rlatg(jjm)

    !-----------------------------------------------------------------

    !!print *, "Call sequence information: bilan_dyn"

    ! Initialisation

    time = time + dt_app
    itau = itau + 1

    first_call: if (first) then
       ! initialisation des fichiers
       first = .false.
       ! ncum est la frequence de stokage en pas de temps
       ncum = dt_cum / dt_app
       if (abs(ncum * dt_app - dt_cum) > 1e-5 * dt_app) then
          print *, 'Problème : le pas de cumul doit être multiple du pas'
          print *, 'dt_app = ', dt_app
          print *, 'dt_cum = ', dt_cum
          stop 1
       endif

       call inigrads(i_f=4, x=(/0./), fx=180./pi, xmin=0., xmax=0., y=rlatv, &
            ymin=-90., ymax=90., fy=180./pi, z=presnivs, fz=1., dt=dt_cum, &
            file='dynzon', titlel='dyn_zon ')

       ! Initialisation du fichier contenant les moyennes zonales

       zan = annee_ref
       dayref = day_ref
       CALL ymds2ju(zan, 1, dayref, 0.0, zjulian)

       rlong = 0.
       rlatg = rlatv*180./pi

       call histbeg_totreg('dynzon', rlong(:1), rlatg, 1, 1, 1, jjm, itau_dyn, &
            zjulian, dt_cum, thoriid, fileid)

       ! Appel à histvert pour la grille verticale

       call histvert(fileid, 'presnivs', 'Niveaux sigma', 'mb', llm, presnivs, &
            zvertiid)

       ! Appels à histdef pour la définition des variables à sauvegarder
       do iQ = 1, nQ
          do itr = 1, ntr
             if (itr == 1) then
                znom(itr, iQ) = nom(iQ)
                znoml(itr, iQ) = nom(iQ)
                zunites(itr, iQ) = unites(iQ)
             else
                znom(itr, iQ) = ctrs(itr)//'v'//nom(iQ)
                znoml(itr, iQ) = 'transport : v * '//nom(iQ)//' '//ctrs(itr)
                zunites(itr, iQ) = 'm/s * '//unites(iQ)
             endif
          enddo
       enddo

       ! Déclarations des champs avec dimension verticale
       do iQ = 1, nQ
          do itr = 1, ntr
             call histdef(fileid, znom(itr, iQ), znoml(itr, iQ), &
                  zunites(itr, iQ), 1, jjm, thoriid, llm, 1, llm, zvertiid, &
                  'ave(X)', dt_cum, dt_cum)
          enddo
          ! Declarations pour les fonctions de courant
          call histdef(fileid, 'psi'//nom(iQ), 'stream fn. '//znoml(itot, iQ), &
               zunites(itot, iQ), 1, jjm, thoriid, llm, 1, llm, zvertiid, &
               'ave(X)', dt_cum, dt_cum)
       enddo

       ! Declarations pour les champs de transport d'air
       call histdef(fileid, 'masse', 'masse', &
            'kg', 1, jjm, thoriid, llm, 1, llm, zvertiid, &
            'ave(X)', dt_cum, dt_cum)
       call histdef(fileid, 'v', 'v', &
            'm/s', 1, jjm, thoriid, llm, 1, llm, zvertiid, &
            'ave(X)', dt_cum, dt_cum)
       ! Declarations pour les fonctions de courant
       call histdef(fileid, 'psi', 'stream fn. MMC ', 'mega t/s', &
            1, jjm, thoriid, llm, 1, llm, zvertiid, &
            'ave(X)', dt_cum, dt_cum)

       ! Declaration des champs 1D de transport en latitude
       do iQ = 1, nQ
          do itr = 2, ntr
             call histdef(fileid, 'a'//znom(itr, iQ), znoml(itr, iQ), &
                  zunites(itr, iQ), 1, jjm, thoriid, 1, 1, 1, -99, &
                  'ave(X)', dt_cum, dt_cum)
          enddo
       enddo

       CALL histend(fileid)
    endif first_call

    ! Calcul des champs dynamiques

    ! Énergie cinétique
    ucont = 0
    CALL covcont(llm, ucov, vcov, ucont, vcont)
    CALL enercin(vcov, ucov, vcont, ucont, ecin)

    ! moment cinétique
    do l = 1, llm
       ang(:, :, l) = ucov(:, :, l) + constang_2d
       unat(:, :, l) = ucont(:, :, l)*cu_2d
    enddo

    Q(:, :, :, 1) = teta * pk / cpp
    Q(:, :, :, 2) = phi
    Q(:, :, :, 3) = ecin
    Q(:, :, :, 4) = ang
    Q(:, :, :, 5) = unat
    Q(:, :, :, 6) = trac
    Q(:, :, :, 7) = 1.

    ! Cumul

    if (icum == 0) then
       ps_cum = 0.
       masse_cum = 0.
       flux_u_cum = 0.
       flux_v_cum = 0.
       Q_cum = 0.
       flux_vQ_cum = 0.
       flux_uQ_cum = 0.
    endif

    icum = icum + 1

    ! Accumulation des flux de masse horizontaux
    ps_cum = ps_cum + ps
    masse_cum = masse_cum + masse
    flux_u_cum = flux_u_cum + flux_u
    flux_v_cum = flux_v_cum + flux_v
    do iQ = 1, nQ
       Q_cum(:, :, :, iQ) = Q_cum(:, :, :, iQ) + Q(:, :, :, iQ)*masse
    enddo

    ! FLUX ET TENDANCES

    ! Flux longitudinal
    forall (iQ = 1: nQ, i = 1: iim) flux_uQ_cum(i, :, :, iQ) &
         = flux_uQ_cum(i, :, :, iQ) &
         + flux_u(i, :, :) * 0.5 * (Q(i, :, :, iQ) + Q(i + 1, :, :, iQ))
    flux_uQ_cum(iip1, :, :, :) = flux_uQ_cum(1, :, :, :)

    ! Flux méridien
    forall (iQ = 1: nQ, j = 1: jjm) flux_vQ_cum(:, j, :, iQ) &
         = flux_vQ_cum(:, j, :, iQ) &
         + flux_v(:, j, :) * 0.5 * (Q(:, j, :, iQ) + Q(:, j + 1, :, iQ))

    ! tendances

    ! convergence horizontale
    call convflu(flux_uQ_cum, flux_vQ_cum, llm*nQ, dQ)

    ! calcul de la vitesse verticale
    call convmas(flux_u_cum, flux_v_cum, convm)
    CALL vitvert(convm, w)

    do iQ = 1, nQ
       do l = 1, llm-1
          do j = 1, jjp1
             do i = 1, iip1
                ww = -0.5*w(i, j, l + 1)*(Q(i, j, l, iQ) + Q(i, j, l + 1, iQ))
                dQ(i, j, l, iQ) = dQ(i, j, l, iQ)-ww
                dQ(i, j, l + 1, iQ) = dQ(i, j, l + 1, iQ) + ww
             enddo
          enddo
       enddo
    enddo

    ! PAS DE TEMPS D'ECRITURE

    writing_step: if (icum == ncum) then
       ! Normalisation
       do iQ = 1, nQ
          Q_cum(:, :, :, iQ) = Q_cum(:, :, :, iQ)/masse_cum
       enddo
       zz = 1. / real(ncum)
       ps_cum = ps_cum*zz
       masse_cum = masse_cum*zz
       flux_u_cum = flux_u_cum*zz
       flux_v_cum = flux_v_cum*zz
       flux_uQ_cum = flux_uQ_cum*zz
       flux_vQ_cum = flux_vQ_cum*zz
       dQ = dQ*zz

       ! A retravailler eventuellement
       ! division de dQ par la masse pour revenir aux bonnes grandeurs
       do iQ = 1, nQ
          dQ(:, :, :, iQ) = dQ(:, :, :, iQ)/masse_cum
       enddo

       ! Transport méridien

       ! cumul zonal des masses des mailles

       zv = 0.
       zmasse = 0.
       call massbar(masse_cum, massebx, masseby)
       do l = 1, llm
          do j = 1, jjm
             do i = 1, iim
                zmasse(j, l) = zmasse(j, l) + masseby(i, j, l)
                zv(j, l) = zv(j, l) + flux_v_cum(i, j, l)
             enddo
             zfactv(j, l) = cv_2d(1, j)/zmasse(j, l)
          enddo
       enddo

       ! Transport dans le plan latitude-altitude

       zvQ = 0.
       psiQ = 0.
       do iQ = 1, nQ
          zvQtmp = 0.
          do l = 1, llm
             do j = 1, jjm
                ! Calcul des moyennes zonales du transort total et de zvQtmp
                do i = 1, iim
                   zvQ(j, l, itot, iQ) = zvQ(j, l, itot, iQ) &
                        + flux_vQ_cum(i, j, l, iQ)
                   zqy =  0.5 * (Q_cum(i, j, l, iQ) * masse_cum(i, j, l) &
                        + Q_cum(i, j + 1, l, iQ) * masse_cum(i, j + 1, l))
                   zvQtmp(j, l) = zvQtmp(j, l) + flux_v_cum(i, j, l) * zqy &
                        / (0.5 * (masse_cum(i, j, l) + masse_cum(i, j + 1, l)))
                   zvQ(j, l, iave, iQ) = zvQ(j, l, iave, iQ) + zqy
                enddo
                ! Decomposition
                zvQ(j, l, iave, iQ) = zvQ(j, l, iave, iQ)/zmasse(j, l)
                zvQ(j, l, itot, iQ) = zvQ(j, l, itot, iQ)*zfactv(j, l)
                zvQtmp(j, l) = zvQtmp(j, l)*zfactv(j, l)
                zvQ(j, l, immc, iQ) = zv(j, l)*zvQ(j, l, iave, iQ)*zfactv(j, l)
                zvQ(j, l, itrs, iQ) = zvQ(j, l, itot, iQ)-zvQtmp(j, l)
                zvQ(j, l, istn, iQ) = zvQtmp(j, l)-zvQ(j, l, immc, iQ)
             enddo
          enddo
          ! fonction de courant meridienne pour la quantite Q
          do l = llm, 1, -1
             do j = 1, jjm
                psiQ(j, l, iQ) = psiQ(j, l + 1, iQ) + zvQ(j, l, itot, iQ)
             enddo
          enddo
       enddo

       ! fonction de courant pour la circulation meridienne moyenne
       psi = 0.
       do l = llm, 1, -1
          do j = 1, jjm
             psi(j, l) = psi(j, l + 1) + zv(j, l)
             zv(j, l) = zv(j, l)*zfactv(j, l)
          enddo
       enddo

       ! sorties proprement dites
       do iQ = 1, nQ
          do itr = 1, ntr
             call histwrite(fileid, znom(itr, iQ), itau, zvQ(:, :, itr, iQ))
          enddo
          call histwrite(fileid, 'psi'//nom(iQ), itau, psiQ(:, :llm, iQ))
       enddo

       call histwrite(fileid, 'masse', itau, zmasse)
       call histwrite(fileid, 'v', itau, zv)
       psi = psi*1.e-9
       call histwrite(fileid, 'psi', itau, psi(:, :llm))

       ! Intégrale verticale

       forall (iQ = 1: nQ, itr = 2: ntr) zavQ(:, itr, iQ) &
            = sum(zvQ(:, :, itr, iQ) * zmasse, dim=2) / cv_2d(1, :)

       do iQ = 1, nQ
          do itr = 2, ntr
             call histwrite(fileid, 'a'//znom(itr, iQ), itau, zavQ(:, itr, iQ))
          enddo
       enddo

       ! On doit pouvoir tracer systematiquement la fonction de courant.
       icum = 0
    endif writing_step

  end SUBROUTINE bilan_dyn

end module bilan_dyn_m
1	guez	40	module bilan_dyn_m
2	guez	3
3	guez	40	IMPLICIT NONE
4	guez	3
5	guez	40	contains
6	guez	3
7	guez	40	SUBROUTINE bilan_dyn(ps, masse, pk, flux_u, flux_v, teta, phi, ucov, vcov, &
8			trac, dt_app, dt_cum)
9	guez	3
10	guez	40	! From LMDZ4/libf/dyn3d/bilan_dyn.F, version 1.5 2005/03/16
11			! 10:12:17 fairhead
12	guez	3
13	guez	55	! Sous-programme consacré à des diagnostics dynamiques de base.
14			! De façon générale, les moyennes des scalaires Q sont pondérées
15			! par la masse. Les flux de masse sont, eux, simplement moyennés.
16	guez	3
17	guez	40	USE histcom, ONLY: histbeg_totreg, histdef, histend, histvert
18			USE calendar, ONLY: ymds2ju
19			USE histwrite_m, ONLY: histwrite
20			USE dimens_m, ONLY: iim, jjm, llm
21			USE paramet_m, ONLY: iip1, jjp1
22			USE comconst, ONLY: cpp
23			USE comvert, ONLY: presnivs
24			USE comgeom, ONLY: constang_2d, cu_2d, cv_2d, rlatv
25			USE temps, ONLY: annee_ref, day_ref, itau_dyn
26			USE inigrads_m, ONLY: inigrads
27			USE nr_util, ONLY: pi
28	guez	3
29	guez	40	! Arguments:
30	guez	3
31	guez	40	real, intent(in):: dt_app, dt_cum
32			real ps(iip1, jjp1)
33			real masse(iip1, jjp1, llm), pk(iip1, jjp1, llm)
34			real flux_u(iip1, jjp1, llm)
35			real flux_v(iip1, jjm, llm)
36	guez	44	real, intent(in):: teta(iip1, jjp1, llm)
37	guez	40	real phi(iip1, jjp1, llm)
38			real ucov(iip1, jjp1, llm)
39			real vcov(iip1, jjm, llm)
40			real, intent(in):: trac(:, :, :) ! (iim + 1, jjm + 1, llm)
41	guez	3
42	guez	40	! Local:
43	guez	3
44	guez	54	integer:: icum = 0
45			integer, save:: ncum
46	guez	40	logical:: first = .true.
47			real zz, zqy, zfactv(jjm, llm)
48	guez	3
49	guez	54	integer, parameter:: nQ = 7
50			character(len=4), parameter:: nom(nQ) = (/'T ', 'gz ', 'K ', 'ang ', &
51			'u ', 'ovap', 'un '/)
52			character(len=5), parameter:: unites(nQ) = (/'K ', 'm2/s2', 'm2/s2', &
53			'ang ', 'm/s ', 'kg/kg', 'un '/)
54	guez	3
55	guez	40	real:: time = 0.
56			integer:: itau = 0
57			real ww
58	guez	3
59	guez	40	! Variables dynamiques intermédiaires
60			REAL vcont(iip1, jjm, llm), ucont(iip1, jjp1, llm)
61			REAL ang(iip1, jjp1, llm), unat(iip1, jjp1, llm)
62			REAL massebx(iip1, jjp1, llm), masseby(iip1, jjm, llm)
63			REAL w(iip1, jjp1, llm), ecin(iip1, jjp1, llm), convm(iip1, jjp1, llm)
64	guez	3
65	guez	40	! Champ contenant les scalaires advectés
66			real Q(iip1, jjp1, llm, nQ)
67	guez	3
68	guez	40	! Champs cumulés
69			real, save:: ps_cum(iip1, jjp1)
70			real, save:: masse_cum(iip1, jjp1, llm)
71			real, save:: flux_u_cum(iip1, jjp1, llm)
72			real, save:: flux_v_cum(iip1, jjm, llm)
73			real, save:: Q_cum(iip1, jjp1, llm, nQ)
74			real, save:: flux_uQ_cum(iip1, jjp1, llm, nQ)
75			real, save:: flux_vQ_cum(iip1, jjm, llm, nQ)
76			real dQ(iip1, jjp1, llm, nQ)
77	guez	3
78	guez	40	! champs de tansport en moyenne zonale
79			integer itr
80	guez	54	integer, parameter:: ntr = 5
81	guez	3
82	guez	40	character(len=10), save:: znom(ntr, nQ)
83	guez	54	character(len=26), save:: znoml(ntr, nQ)
84			character(len=12), save:: zunites(ntr, nQ)
85	guez	3
86	guez	54	integer, parameter:: iave = 1, itot = 2, immc = 3, itrs = 4, istn = 5
87			character(len=3), parameter:: ctrs(ntr) = (/' ', 'TOT', 'MMC', 'TRS', &
88			'STN'/)
89	guez	3
90	guez	40	real zvQ(jjm, llm, ntr, nQ), zvQtmp(jjm, llm)
91	guez	54	real zavQ(jjm, 2: ntr, nQ), psiQ(jjm, llm + 1, nQ)
92			real zmasse(jjm, llm)
93	guez	3
94	guez	54	real zv(jjm, llm), psi(jjm, llm + 1)
95	guez	3
96	guez	40	integer i, j, l, iQ
97	guez	3
98	guez	40	! Initialisation du fichier contenant les moyennes zonales.
99	guez	3
100	guez	40	integer, save:: fileid
101			integer thoriid, zvertiid
102	guez	3
103	guez	40	real zjulian
104			integer zan, dayref
105	guez	3
106	guez	40	real rlong(jjm), rlatg(jjm)
107	guez	3
108	guez	40	!-----------------------------------------------------------------
109	guez	3
110	guez	40	!!print *, "Call sequence information: bilan_dyn"
111	guez	3
112	guez	40	! Initialisation
113	guez	3
114	guez	54	time = time + dt_app
115			itau = itau + 1
116	guez	3
117	guez	54	first_call: if (first) then
118	guez	40	! initialisation des fichiers
119	guez	54	first = .false.
120	guez	40	! ncum est la frequence de stokage en pas de temps
121	guez	54	ncum = dt_cum / dt_app
122	guez	40	if (abs(ncum * dt_app - dt_cum) > 1e-5 * dt_app) then
123			print *, 'Problème : le pas de cumul doit être multiple du pas'
124	guez	54	print *, 'dt_app = ', dt_app
125			print *, 'dt_cum = ', dt_cum
126	guez	40	stop 1
127			endif
128	guez	3
129	guez	54	call inigrads(i_f=4, x=(/0./), fx=180./pi, xmin=0., xmax=0., y=rlatv, &
130			ymin=-90., ymax=90., fy=180./pi, z=presnivs, fz=1., dt=dt_cum, &
131			file='dynzon', titlel='dyn_zon ')
132	guez	3
133	guez	40	! Initialisation du fichier contenant les moyennes zonales
134	guez	3
135	guez	40	zan = annee_ref
136			dayref = day_ref
137			CALL ymds2ju(zan, 1, dayref, 0.0, zjulian)
138	guez	3
139	guez	54	rlong = 0.
140			rlatg = rlatv*180./pi
141	guez	3
142	guez	40	call histbeg_totreg('dynzon', rlong(:1), rlatg, 1, 1, 1, jjm, itau_dyn, &
143			zjulian, dt_cum, thoriid, fileid)
144	guez	3
145	guez	40	! Appel à histvert pour la grille verticale
146	guez	3
147	guez	40	call histvert(fileid, 'presnivs', 'Niveaux sigma', 'mb', llm, presnivs, &
148			zvertiid)
149	guez	3
150	guez	40	! Appels à histdef pour la définition des variables à sauvegarder
151	guez	54	do iQ = 1, nQ
152			do itr = 1, ntr
153			if (itr == 1) then
154			znom(itr, iQ) = nom(iQ)
155			znoml(itr, iQ) = nom(iQ)
156			zunites(itr, iQ) = unites(iQ)
157	guez	40	else
158	guez	54	znom(itr, iQ) = ctrs(itr)//'v'//nom(iQ)
159			znoml(itr, iQ) = 'transport : v * '//nom(iQ)//' '//ctrs(itr)
160			zunites(itr, iQ) = 'm/s * '//unites(iQ)
161	guez	40	endif
162			enddo
163			enddo
164	guez	3
165	guez	40	! Déclarations des champs avec dimension verticale
166	guez	54	do iQ = 1, nQ
167			do itr = 1, ntr
168	guez	40	call histdef(fileid, znom(itr, iQ), znoml(itr, iQ), &
169			zunites(itr, iQ), 1, jjm, thoriid, llm, 1, llm, zvertiid, &
170			'ave(X)', dt_cum, dt_cum)
171			enddo
172			! Declarations pour les fonctions de courant
173			call histdef(fileid, 'psi'//nom(iQ), 'stream fn. '//znoml(itot, iQ), &
174			zunites(itot, iQ), 1, jjm, thoriid, llm, 1, llm, zvertiid, &
175			'ave(X)', dt_cum, dt_cum)
176			enddo
177	guez	3
178	guez	40	! Declarations pour les champs de transport d'air
179			call histdef(fileid, 'masse', 'masse', &
180			'kg', 1, jjm, thoriid, llm, 1, llm, zvertiid, &
181			'ave(X)', dt_cum, dt_cum)
182			call histdef(fileid, 'v', 'v', &
183			'm/s', 1, jjm, thoriid, llm, 1, llm, zvertiid, &
184			'ave(X)', dt_cum, dt_cum)
185			! Declarations pour les fonctions de courant
186			call histdef(fileid, 'psi', 'stream fn. MMC ', 'mega t/s', &
187			1, jjm, thoriid, llm, 1, llm, zvertiid, &
188			'ave(X)', dt_cum, dt_cum)
189	guez	3
190	guez	40	! Declaration des champs 1D de transport en latitude
191	guez	54	do iQ = 1, nQ
192			do itr = 2, ntr
193	guez	40	call histdef(fileid, 'a'//znom(itr, iQ), znoml(itr, iQ), &
194			zunites(itr, iQ), 1, jjm, thoriid, 1, 1, 1, -99, &
195			'ave(X)', dt_cum, dt_cum)
196			enddo
197			enddo
198	guez	3
199	guez	40	CALL histend(fileid)
200	guez	54	endif first_call
201	guez	3
202	guez	40	! Calcul des champs dynamiques
203	guez	3
204	guez	40	! Énergie cinétique
205			ucont = 0
206			CALL covcont(llm, ucov, vcov, ucont, vcont)
207			CALL enercin(vcov, ucov, vcont, ucont, ecin)
208	guez	3
209	guez	40	! moment cinétique
210	guez	54	do l = 1, llm
211			ang(:, :, l) = ucov(:, :, l) + constang_2d
212			unat(:, :, l) = ucont(:, :, l)*cu_2d
213	guez	40	enddo
214	guez	3
215	guez	54	Q(:, :, :, 1) = teta * pk / cpp
216			Q(:, :, :, 2) = phi
217			Q(:, :, :, 3) = ecin
218			Q(:, :, :, 4) = ang
219			Q(:, :, :, 5) = unat
220			Q(:, :, :, 6) = trac
221			Q(:, :, :, 7) = 1.
222	guez	3
223	guez	40	! Cumul
224	guez	3
225	guez	54	if (icum == 0) then
226			ps_cum = 0.
227			masse_cum = 0.
228			flux_u_cum = 0.
229			flux_v_cum = 0.
230			Q_cum = 0.
231			flux_vQ_cum = 0.
232			flux_uQ_cum = 0.
233	guez	40	endif
234	guez	3
235	guez	54	icum = icum + 1
236	guez	3
237	guez	40	! Accumulation des flux de masse horizontaux
238	guez	54	ps_cum = ps_cum + ps
239			masse_cum = masse_cum + masse
240			flux_u_cum = flux_u_cum + flux_u
241			flux_v_cum = flux_v_cum + flux_v
242			do iQ = 1, nQ
243			Q_cum(:, :, :, iQ) = Q_cum(:, :, :, iQ) + Q(:, :, :, iQ)*masse
244	guez	40	enddo
245	guez	3
246	guez	40	! FLUX ET TENDANCES
247	guez	3
248	guez	40	! Flux longitudinal
249	guez	54	forall (iQ = 1: nQ, i = 1: iim) flux_uQ_cum(i, :, :, iQ) &
250			= flux_uQ_cum(i, :, :, iQ) &
251			+ flux_u(i, :, :) * 0.5 * (Q(i, :, :, iQ) + Q(i + 1, :, :, iQ))
252			flux_uQ_cum(iip1, :, :, :) = flux_uQ_cum(1, :, :, :)
253	guez	3
254	guez	54	! Flux méridien
255			forall (iQ = 1: nQ, j = 1: jjm) flux_vQ_cum(:, j, :, iQ) &
256			= flux_vQ_cum(:, j, :, iQ) &
257			+ flux_v(:, j, :) * 0.5 * (Q(:, j, :, iQ) + Q(:, j + 1, :, iQ))
258	guez	3
259	guez	40	! tendances
260	guez	3
261	guez	40	! convergence horizontale
262			call convflu(flux_uQ_cum, flux_vQ_cum, llm*nQ, dQ)
263	guez	3
264	guez	40	! calcul de la vitesse verticale
265			call convmas(flux_u_cum, flux_v_cum, convm)
266			CALL vitvert(convm, w)
267	guez	3
268	guez	54	do iQ = 1, nQ
269			do l = 1, llm-1
270			do j = 1, jjp1
271			do i = 1, iip1
272			ww = -0.5w(i, j, l + 1)(Q(i, j, l, iQ) + Q(i, j, l + 1, iQ))
273			dQ(i, j, l, iQ) = dQ(i, j, l, iQ)-ww
274			dQ(i, j, l + 1, iQ) = dQ(i, j, l + 1, iQ) + ww
275	guez	40	enddo
276			enddo
277			enddo
278			enddo
279	guez	3
280	guez	40	! PAS DE TEMPS D'ECRITURE
281	guez	3
282	guez	40	writing_step: if (icum == ncum) then
283			! Normalisation
284	guez	54	do iQ = 1, nQ
285			Q_cum(:, :, :, iQ) = Q_cum(:, :, :, iQ)/masse_cum
286	guez	40	enddo
287	guez	54	zz = 1. / real(ncum)
288			ps_cum = ps_cum*zz
289			masse_cum = masse_cum*zz
290			flux_u_cum = flux_u_cum*zz
291			flux_v_cum = flux_v_cum*zz
292			flux_uQ_cum = flux_uQ_cum*zz
293			flux_vQ_cum = flux_vQ_cum*zz
294			dQ = dQ*zz
295	guez	3
296	guez	40	! A retravailler eventuellement
297			! division de dQ par la masse pour revenir aux bonnes grandeurs
298	guez	54	do iQ = 1, nQ
299			dQ(:, :, :, iQ) = dQ(:, :, :, iQ)/masse_cum
300	guez	40	enddo
301	guez	3
302	guez	40	! Transport méridien
303	guez	3
304	guez	40	! cumul zonal des masses des mailles
305	guez	3
306	guez	54	zv = 0.
307			zmasse = 0.
308	guez	40	call massbar(masse_cum, massebx, masseby)
309	guez	54	do l = 1, llm
310			do j = 1, jjm
311			do i = 1, iim
312			zmasse(j, l) = zmasse(j, l) + masseby(i, j, l)
313			zv(j, l) = zv(j, l) + flux_v_cum(i, j, l)
314	guez	40	enddo
315	guez	54	zfactv(j, l) = cv_2d(1, j)/zmasse(j, l)
316	guez	40	enddo
317			enddo
318	guez	3
319	guez	40	! Transport dans le plan latitude-altitude
320	guez	3
321	guez	54	zvQ = 0.
322			psiQ = 0.
323			do iQ = 1, nQ
324			zvQtmp = 0.
325			do l = 1, llm
326			do j = 1, jjm
327	guez	40	! Calcul des moyennes zonales du transort total et de zvQtmp
328	guez	54	do i = 1, iim
329			zvQ(j, l, itot, iQ) = zvQ(j, l, itot, iQ) &
330			+ flux_vQ_cum(i, j, l, iQ)
331			zqy = 0.5 * (Q_cum(i, j, l, iQ) * masse_cum(i, j, l) &
332			+ Q_cum(i, j + 1, l, iQ) * masse_cum(i, j + 1, l))
333			zvQtmp(j, l) = zvQtmp(j, l) + flux_v_cum(i, j, l) * zqy &
334			/ (0.5 * (masse_cum(i, j, l) + masse_cum(i, j + 1, l)))
335			zvQ(j, l, iave, iQ) = zvQ(j, l, iave, iQ) + zqy
336	guez	40	enddo
337			! Decomposition
338	guez	54	zvQ(j, l, iave, iQ) = zvQ(j, l, iave, iQ)/zmasse(j, l)
339			zvQ(j, l, itot, iQ) = zvQ(j, l, itot, iQ)*zfactv(j, l)
340			zvQtmp(j, l) = zvQtmp(j, l)*zfactv(j, l)
341			zvQ(j, l, immc, iQ) = zv(j, l)zvQ(j, l, iave, iQ)zfactv(j, l)
342			zvQ(j, l, itrs, iQ) = zvQ(j, l, itot, iQ)-zvQtmp(j, l)
343			zvQ(j, l, istn, iQ) = zvQtmp(j, l)-zvQ(j, l, immc, iQ)
344	guez	40	enddo
345			enddo
346			! fonction de courant meridienne pour la quantite Q
347	guez	54	do l = llm, 1, -1
348			do j = 1, jjm
349			psiQ(j, l, iQ) = psiQ(j, l + 1, iQ) + zvQ(j, l, itot, iQ)
350	guez	40	enddo
351			enddo
352			enddo
353	guez	3
354	guez	40	! fonction de courant pour la circulation meridienne moyenne
355	guez	54	psi = 0.
356			do l = llm, 1, -1
357			do j = 1, jjm
358			psi(j, l) = psi(j, l + 1) + zv(j, l)
359			zv(j, l) = zv(j, l)*zfactv(j, l)
360	guez	40	enddo
361			enddo
362	guez	3
363	guez	40	! sorties proprement dites
364	guez	54	do iQ = 1, nQ
365			do itr = 1, ntr
366			call histwrite(fileid, znom(itr, iQ), itau, zvQ(:, :, itr, iQ))
367	guez	40	enddo
368	guez	54	call histwrite(fileid, 'psi'//nom(iQ), itau, psiQ(:, :llm, iQ))
369			enddo
370	guez	3
371	guez	54	call histwrite(fileid, 'masse', itau, zmasse)
372			call histwrite(fileid, 'v', itau, zv)
373			psi = psi*1.e-9
374			call histwrite(fileid, 'psi', itau, psi(:, :llm))
375	guez	3
376	guez	55	! Intégrale verticale
377	guez	3
378	guez	54	forall (iQ = 1: nQ, itr = 2: ntr) zavQ(:, itr, iQ) &
379	guez	55	= sum(zvQ(:, :, itr, iQ) * zmasse, dim=2) / cv_2d(1, :)
380	guez	54
381			do iQ = 1, nQ
382			do itr = 2, ntr
383	guez	40	call histwrite(fileid, 'a'//znom(itr, iQ), itau, zavQ(:, itr, iQ))
384			enddo
385			enddo
386	guez	3
387	guez	54	! On doit pouvoir tracer systematiquement la fonction de courant.
388			icum = 0
389	guez	40	endif writing_step
390	guez	3
391	guez	40	end SUBROUTINE bilan_dyn
392	guez	3
393	guez	40	end module bilan_dyn_m