trunk/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)

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

    ! 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 calendar, ONLY: ymds2ju
    USE conf_gcm_m, ONLY: day_step, iperiod, periodav
    USE comconst, ONLY: cpp
    USE comvert, ONLY: presnivs
    USE comgeom, ONLY: constang_2d, cu_2d, cv_2d, rlatv
    USE dimens_m, ONLY: iim, jjm, llm
    USE histcom, ONLY: histbeg_totreg, histdef, histend, histvert
    USE histwrite_m, ONLY: histwrite
    USE nr_util, ONLY: pi
    USE paramet_m, ONLY: iip1, jjp1
    USE temps, ONLY: annee_ref, day_ref, itau_dyn

    ! Arguments:

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

    ! Local:

    real dt_cum
    integer:: icum  = 0
    integer, save:: ncum
    logical:: first = .true.
    real 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   '/)

    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"

    first_call: if (first) then
       ! initialisation des fichiers
       first = .false.
       ! ncum est la frequence de stokage en pas de temps
       ncum = day_step / iperiod * periodav
       dt_cum = ncum * dt_app

       ! 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
          ! Déclarations 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

       ! Déclarations 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)
       ! Déclarations 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)

       ! Déclaration 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

    itau = itau + 1

    ! 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
       ps_cum = ps_cum / ncum
       masse_cum = masse_cum / ncum
       flux_u_cum = flux_u_cum / ncum
       flux_v_cum = flux_v_cum / ncum
       flux_uQ_cum = flux_uQ_cum / ncum
       flux_vQ_cum = flux_vQ_cum / ncum
       dQ = dQ / ncum

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