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, 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, save:: icum, ncum
    logical:: first = .true.
    real zz, zqy, zfactv(jjm, llm)

    integer, parameter:: nQ=7

    character(len=6), save:: nom(nQ)
    character(len=6), save:: unites(nQ)

    character(len=10) file
    integer, parameter:: ifile=4

    integer itemp, igeop, iecin, iang, iu, iovap, iun
    integer:: i_sortie = 1

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

    data itemp, igeop, iecin, iang, iu, iovap, iun/1, 2, 3, 4, 5, 6, 7/

    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 vorpot(iip1, jjm, llm)
    REAL w(iip1, jjp1, llm), ecin(iip1, jjp1, llm), convm(iip1, jjp1, llm)
    REAL bern(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 flux_wQ_cum(iip1, jjp1, 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=20), save:: znoml(ntr, nQ)
    character(len=10), save:: zunites(ntr, nQ)

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

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

    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
    integer ndex3d(jjm*llm)

    ! Variables locales

    real zjulian
    character(len=3) str
    character(len=10) ctrac
    integer ii, jj
    integer zan, dayref

    real rlong(jjm), rlatg(jjm)

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

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

    ! Initialisation

    time=time+dt_app
    itau=itau+1

    if (first) then
       icum=0
       ! 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

       if (i_sortie == 1) then
          file='dynzon'
          call inigrads(ifile , (/0./), 180./pi, 0., 0., rlatv, -90., 90., &
               180./pi , presnivs, 1. , dt_cum, file, 'dyn_zon ')
       endif

       nom(itemp)='T'
       nom(igeop)='gz'
       nom(iecin)='K'
       nom(iang)='ang'
       nom(iu)='u'
       nom(iovap)='ovap'
       nom(iun)='un'

       unites(itemp)='K'
       unites(igeop)='m2/s2'
       unites(iecin)='m2/s2'
       unites(iang)='ang'
       unites(iu)='m/s'
       unites(iovap)='kg/kg'
       unites(iun)='un'

       ! 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

    ! 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(:, :, :, itemp)=teta*pk/cpp
    Q(:, :, :, igeop)=phi
    Q(:, :, :, iecin)=ecin
    Q(:, :, :, iang)=ang
    Q(:, :, :, iu)=unat
    Q(:, :, :, iovap)=trac
    Q(:, :, :, iun)=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
    do iQ=1, nQ
       do l=1, llm
          do j=1, jjp1
             do i=1, iim
                flux_uQ_cum(i, j, l, iQ)=flux_uQ_cum(i, j, l, iQ) &
                     +flux_u(i, j, l)*0.5*(Q(i, j, l, iQ)+Q(i+1, j, l, iQ))
             enddo
             flux_uQ_cum(iip1, j, l, iQ)=flux_uQ_cum(1, j, l, iQ)
          enddo
       enddo
    enddo

    ! flux méridien
    do iQ=1, nQ
       do l=1, llm
          do j=1, jjm
             do i=1, iip1
                flux_vQ_cum(i, j, l, iQ)=flux_vQ_cum(i, j, l, iQ) &
                     +flux_v(i, j, l)*0.5*(Q(i, j, l, iQ)+Q(i, j+1, l, iQ))
             enddo
          enddo
       enddo
    enddo

    ! 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./float(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
       if (i_sortie == 1) then
          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(:, 1: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(:, 1:llm))
       endif

       ! Moyenne verticale

       zamasse=0.
       do l=1, llm
          zamasse(:)=zamasse(:)+zmasse(:, l)
       enddo
       zavQ=0.
       do iQ=1, nQ
          do itr=2, ntr
             do l=1, llm
                zavQ(:, itr, iQ) = zavQ(:, itr, iQ) &
                     + zvQ(:, l, itr, iQ) * zmasse(:, l)
             enddo
             zavQ(:, itr, iQ)=zavQ(:, itr, iQ)/zamasse(:)
             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	40	! Sous-programme consacré à des diagnostics dynamiques de base
14			! De façon générale, les moyennes des scalaires Q sont pondérées par
15			! 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	40	integer, save:: icum, ncum
45			logical:: first = .true.
46			real zz, zqy, zfactv(jjm, llm)
47	guez	3
48	guez	40	integer, parameter:: nQ=7
49	guez	3
50	guez	40	character(len=6), save:: nom(nQ)
51			character(len=6), save:: unites(nQ)
52	guez	3
53	guez	40	character(len=10) file
54			integer, parameter:: ifile=4
55	guez	3
56	guez	40	integer itemp, igeop, iecin, iang, iu, iovap, iun
57			integer:: i_sortie = 1
58	guez	3
59	guez	40	real:: time = 0.
60			integer:: itau = 0
61	guez	3
62	guez	40	data itemp, igeop, iecin, iang, iu, iovap, iun/1, 2, 3, 4, 5, 6, 7/
63	guez	3
64	guez	40	real ww
65	guez	3
66	guez	40	! Variables dynamiques intermédiaires
67			REAL vcont(iip1, jjm, llm), ucont(iip1, jjp1, llm)
68			REAL ang(iip1, jjp1, llm), unat(iip1, jjp1, llm)
69			REAL massebx(iip1, jjp1, llm), masseby(iip1, jjm, llm)
70			REAL vorpot(iip1, jjm, llm)
71			REAL w(iip1, jjp1, llm), ecin(iip1, jjp1, llm), convm(iip1, jjp1, llm)
72			REAL bern(iip1, jjp1, llm)
73	guez	3
74	guez	40	! Champ contenant les scalaires advectés
75			real Q(iip1, jjp1, llm, nQ)
76	guez	3
77	guez	40	! Champs cumulés
78			real, save:: ps_cum(iip1, jjp1)
79			real, save:: masse_cum(iip1, jjp1, llm)
80			real, save:: flux_u_cum(iip1, jjp1, llm)
81			real, save:: flux_v_cum(iip1, jjm, llm)
82			real, save:: Q_cum(iip1, jjp1, llm, nQ)
83			real, save:: flux_uQ_cum(iip1, jjp1, llm, nQ)
84			real, save:: flux_vQ_cum(iip1, jjm, llm, nQ)
85			real flux_wQ_cum(iip1, jjp1, llm, nQ)
86			real dQ(iip1, jjp1, llm, nQ)
87	guez	3
88	guez	40	! champs de tansport en moyenne zonale
89			integer itr
90			integer, parameter:: ntr=5
91	guez	3
92	guez	40	character(len=10), save:: znom(ntr, nQ)
93			character(len=20), save:: znoml(ntr, nQ)
94			character(len=10), save:: zunites(ntr, nQ)
95	guez	3
96	guez	40	integer iave, itot, immc, itrs, istn
97			data iave, itot, immc, itrs, istn/1, 2, 3, 4, 5/
98			character(len=3) ctrs(ntr)
99			data ctrs/' ', 'TOT', 'MMC', 'TRS', 'STN'/
100	guez	3
101	guez	40	real zvQ(jjm, llm, ntr, nQ), zvQtmp(jjm, llm)
102			real zavQ(jjm, ntr, nQ), psiQ(jjm, llm+1, nQ)
103			real zmasse(jjm, llm), zamasse(jjm)
104	guez	3
105	guez	40	real zv(jjm, llm), psi(jjm, llm+1)
106	guez	3
107	guez	40	integer i, j, l, iQ
108	guez	3
109	guez	40	! Initialisation du fichier contenant les moyennes zonales.
110	guez	3
111	guez	40	integer, save:: fileid
112			integer thoriid, zvertiid
113			integer ndex3d(jjm*llm)
114	guez	3
115	guez	40	! Variables locales
116	guez	3
117	guez	40	real zjulian
118			character(len=3) str
119			character(len=10) ctrac
120			integer ii, jj
121			integer zan, dayref
122	guez	3
123	guez	40	real rlong(jjm), rlatg(jjm)
124	guez	3
125	guez	40	!-----------------------------------------------------------------
126	guez	3
127	guez	40	!!print *, "Call sequence information: bilan_dyn"
128	guez	3
129	guez	40	! Initialisation
130	guez	3
131	guez	40	time=time+dt_app
132			itau=itau+1
133	guez	3
134	guez	40	if (first) then
135			icum=0
136			! initialisation des fichiers
137			first=.false.
138			! ncum est la frequence de stokage en pas de temps
139			ncum=dt_cum/dt_app
140			if (abs(ncum * dt_app - dt_cum) > 1e-5 * dt_app) then
141			print *, 'Problème : le pas de cumul doit être multiple du pas'
142			print *, 'dt_app=', dt_app
143			print *, 'dt_cum=', dt_cum
144			stop 1
145			endif
146	guez	3
147	guez	40	if (i_sortie == 1) then
148			file='dynzon'
149			call inigrads(ifile , (/0./), 180./pi, 0., 0., rlatv, -90., 90., &
150			180./pi , presnivs, 1. , dt_cum, file, 'dyn_zon ')
151			endif
152	guez	3
153	guez	40	nom(itemp)='T'
154			nom(igeop)='gz'
155			nom(iecin)='K'
156			nom(iang)='ang'
157			nom(iu)='u'
158			nom(iovap)='ovap'
159			nom(iun)='un'
160	guez	3
161	guez	40	unites(itemp)='K'
162			unites(igeop)='m2/s2'
163			unites(iecin)='m2/s2'
164			unites(iang)='ang'
165			unites(iu)='m/s'
166			unites(iovap)='kg/kg'
167			unites(iun)='un'
168	guez	3
169	guez	40	! Initialisation du fichier contenant les moyennes zonales
170	guez	3
171	guez	40	zan = annee_ref
172			dayref = day_ref
173			CALL ymds2ju(zan, 1, dayref, 0.0, zjulian)
174	guez	3
175	guez	40	rlong=0.
176			rlatg=rlatv*180./pi
177	guez	3
178	guez	40	call histbeg_totreg('dynzon', rlong(:1), rlatg, 1, 1, 1, jjm, itau_dyn, &
179			zjulian, dt_cum, thoriid, fileid)
180	guez	3
181	guez	40	! Appel à histvert pour la grille verticale
182	guez	3
183	guez	40	call histvert(fileid, 'presnivs', 'Niveaux sigma', 'mb', llm, presnivs, &
184			zvertiid)
185	guez	3
186	guez	40	! Appels à histdef pour la définition des variables à sauvegarder
187			do iQ=1, nQ
188			do itr=1, ntr
189			if(itr == 1) then
190			znom(itr, iQ)=nom(iQ)
191			znoml(itr, iQ)=nom(iQ)
192			zunites(itr, iQ)=unites(iQ)
193			else
194			znom(itr, iQ)=ctrs(itr)//'v'//nom(iQ)
195			znoml(itr, iQ)='transport : v * '//nom(iQ)//' '//ctrs(itr)
196			zunites(itr, iQ)='m/s * '//unites(iQ)
197			endif
198			enddo
199			enddo
200	guez	3
201	guez	40	! Déclarations des champs avec dimension verticale
202			do iQ=1, nQ
203			do itr=1, ntr
204			call histdef(fileid, znom(itr, iQ), znoml(itr, iQ), &
205			zunites(itr, iQ), 1, jjm, thoriid, llm, 1, llm, zvertiid, &
206			'ave(X)', dt_cum, dt_cum)
207			enddo
208			! Declarations pour les fonctions de courant
209			call histdef(fileid, 'psi'//nom(iQ), 'stream fn. '//znoml(itot, iQ), &
210			zunites(itot, iQ), 1, jjm, thoriid, llm, 1, llm, zvertiid, &
211			'ave(X)', dt_cum, dt_cum)
212			enddo
213	guez	3
214	guez	40	! Declarations pour les champs de transport d'air
215			call histdef(fileid, 'masse', 'masse', &
216			'kg', 1, jjm, thoriid, llm, 1, llm, zvertiid, &
217			'ave(X)', dt_cum, dt_cum)
218			call histdef(fileid, 'v', 'v', &
219			'm/s', 1, jjm, thoriid, llm, 1, llm, zvertiid, &
220			'ave(X)', dt_cum, dt_cum)
221			! Declarations pour les fonctions de courant
222			call histdef(fileid, 'psi', 'stream fn. MMC ', 'mega t/s', &
223			1, jjm, thoriid, llm, 1, llm, zvertiid, &
224			'ave(X)', dt_cum, dt_cum)
225	guez	3
226	guez	40	! Declaration des champs 1D de transport en latitude
227			do iQ=1, nQ
228			do itr=2, ntr
229			call histdef(fileid, 'a'//znom(itr, iQ), znoml(itr, iQ), &
230			zunites(itr, iQ), 1, jjm, thoriid, 1, 1, 1, -99, &
231			'ave(X)', dt_cum, dt_cum)
232			enddo
233			enddo
234	guez	3
235	guez	40	CALL histend(fileid)
236			endif
237	guez	3
238	guez	40	! Calcul des champs dynamiques
239	guez	3
240	guez	40	! Énergie cinétique
241			ucont = 0
242			CALL covcont(llm, ucov, vcov, ucont, vcont)
243			CALL enercin(vcov, ucov, vcont, ucont, ecin)
244	guez	3
245	guez	40	! moment cinétique
246			do l=1, llm
247			ang(:, :, l)=ucov(:, :, l)+constang_2d
248			unat(:, :, l)=ucont(:, :, l)*cu_2d
249			enddo
250	guez	3
251	guez	40	Q(:, :, :, itemp)=teta*pk/cpp
252			Q(:, :, :, igeop)=phi
253			Q(:, :, :, iecin)=ecin
254			Q(:, :, :, iang)=ang
255			Q(:, :, :, iu)=unat
256			Q(:, :, :, iovap)=trac
257			Q(:, :, :, iun)=1.
258	guez	3
259	guez	40	! Cumul
260	guez	3
261	guez	40	if(icum == 0) then
262			ps_cum=0.
263			masse_cum=0.
264			flux_u_cum=0.
265			flux_v_cum=0.
266			Q_cum=0.
267			flux_vQ_cum=0.
268			flux_uQ_cum=0.
269			endif
270	guez	3
271	guez	40	icum=icum+1
272	guez	3
273	guez	40	! Accumulation des flux de masse horizontaux
274			ps_cum=ps_cum+ps
275			masse_cum=masse_cum+masse
276			flux_u_cum=flux_u_cum+flux_u
277			flux_v_cum=flux_v_cum+flux_v
278			do iQ=1, nQ
279			Q_cum(:, :, :, iQ)=Q_cum(:, :, :, iQ)+Q(:, :, :, iQ)*masse
280			enddo
281	guez	3
282	guez	40	! FLUX ET TENDANCES
283	guez	3
284	guez	40	! Flux longitudinal
285			do iQ=1, nQ
286			do l=1, llm
287			do j=1, jjp1
288			do i=1, iim
289			flux_uQ_cum(i, j, l, iQ)=flux_uQ_cum(i, j, l, iQ) &
290			+flux_u(i, j, l)0.5(Q(i, j, l, iQ)+Q(i+1, j, l, iQ))
291			enddo
292			flux_uQ_cum(iip1, j, l, iQ)=flux_uQ_cum(1, j, l, iQ)
293			enddo
294			enddo
295			enddo
296	guez	3
297	guez	40	! flux méridien
298			do iQ=1, nQ
299			do l=1, llm
300			do j=1, jjm
301			do i=1, iip1
302			flux_vQ_cum(i, j, l, iQ)=flux_vQ_cum(i, j, l, iQ) &
303			+flux_v(i, j, l)0.5(Q(i, j, l, iQ)+Q(i, j+1, l, iQ))
304			enddo
305			enddo
306			enddo
307			enddo
308	guez	3
309	guez	40	! tendances
310	guez	3
311	guez	40	! convergence horizontale
312			call convflu(flux_uQ_cum, flux_vQ_cum, llm*nQ, dQ)
313	guez	3
314	guez	40	! calcul de la vitesse verticale
315			call convmas(flux_u_cum, flux_v_cum, convm)
316			CALL vitvert(convm, w)
317	guez	3
318	guez	40	do iQ=1, nQ
319			do l=1, llm-1
320			do j=1, jjp1
321			do i=1, iip1
322			ww=-0.5w(i, j, l+1)(Q(i, j, l, iQ)+Q(i, j, l+1, iQ))
323			dQ(i, j, l , iQ)=dQ(i, j, l , iQ)-ww
324			dQ(i, j, l+1, iQ)=dQ(i, j, l+1, iQ)+ww
325			enddo
326			enddo
327			enddo
328			enddo
329	guez	3
330	guez	40	! PAS DE TEMPS D'ECRITURE
331	guez	3
332	guez	40	writing_step: if (icum == ncum) then
333			! Normalisation
334			do iQ=1, nQ
335			Q_cum(:, :, :, iQ)=Q_cum(:, :, :, iQ)/masse_cum
336			enddo
337			zz=1./float(ncum)
338			ps_cum=ps_cum*zz
339			masse_cum=masse_cum*zz
340			flux_u_cum=flux_u_cum*zz
341			flux_v_cum=flux_v_cum*zz
342			flux_uQ_cum=flux_uQ_cum*zz
343			flux_vQ_cum=flux_vQ_cum*zz
344			dQ=dQ*zz
345	guez	3
346	guez	40	! A retravailler eventuellement
347			! division de dQ par la masse pour revenir aux bonnes grandeurs
348			do iQ=1, nQ
349			dQ(:, :, :, iQ)=dQ(:, :, :, iQ)/masse_cum
350			enddo
351	guez	3
352	guez	40	! Transport méridien
353	guez	3
354	guez	40	! cumul zonal des masses des mailles
355	guez	3
356	guez	40	zv=0.
357			zmasse=0.
358			call massbar(masse_cum, massebx, masseby)
359			do l=1, llm
360			do j=1, jjm
361			do i=1, iim
362			zmasse(j, l)=zmasse(j, l)+masseby(i, j, l)
363			zv(j, l)=zv(j, l)+flux_v_cum(i, j, l)
364			enddo
365			zfactv(j, l)=cv_2d(1, j)/zmasse(j, l)
366			enddo
367			enddo
368	guez	3
369	guez	40	! Transport dans le plan latitude-altitude
370	guez	3
371	guez	40	zvQ=0.
372			psiQ=0.
373			do iQ=1, nQ
374			zvQtmp=0.
375			do l=1, llm
376			do j=1, jjm
377			! Calcul des moyennes zonales du transort total et de zvQtmp
378			do i=1, iim
379			zvQ(j, l, itot, iQ)=zvQ(j, l, itot, iQ) &
380			+flux_vQ_cum(i, j, l, iQ)
381			zqy= 0.5(Q_cum(i, j, l, iQ)masse_cum(i, j, l)+ &
382			Q_cum(i, j+1, l, iQ)*masse_cum(i, j+1, l))
383			zvQtmp(j, l)=zvQtmp(j, l)+flux_v_cum(i, j, l)*zqy &
384			/(0.5*(masse_cum(i, j, l)+masse_cum(i, j+1, l)))
385			zvQ(j, l, iave, iQ)=zvQ(j, l, iave, iQ)+zqy
386			enddo
387			! Decomposition
388			zvQ(j, l, iave, iQ)=zvQ(j, l, iave, iQ)/zmasse(j, l)
389			zvQ(j, l, itot, iQ)=zvQ(j, l, itot, iQ)*zfactv(j, l)
390			zvQtmp(j, l)=zvQtmp(j, l)*zfactv(j, l)
391			zvQ(j, l, immc, iQ)=zv(j, l)zvQ(j, l, iave, iQ)zfactv(j, l)
392			zvQ(j, l, itrs, iQ)=zvQ(j, l, itot, iQ)-zvQtmp(j, l)
393			zvQ(j, l, istn, iQ)=zvQtmp(j, l)-zvQ(j, l, immc, iQ)
394			enddo
395			enddo
396			! fonction de courant meridienne pour la quantite Q
397			do l=llm, 1, -1
398			do j=1, jjm
399			psiQ(j, l, iQ)=psiQ(j, l+1, iQ)+zvQ(j, l, itot, iQ)
400			enddo
401			enddo
402			enddo
403	guez	3
404	guez	40	! fonction de courant pour la circulation meridienne moyenne
405			psi=0.
406			do l=llm, 1, -1
407			do j=1, jjm
408			psi(j, l)=psi(j, l+1)+zv(j, l)
409			zv(j, l)=zv(j, l)*zfactv(j, l)
410			enddo
411			enddo
412	guez	3
413	guez	40	! sorties proprement dites
414			if (i_sortie == 1) then
415			do iQ=1, nQ
416			do itr=1, ntr
417			call histwrite(fileid, znom(itr, iQ), itau, zvQ(:, :, itr, iQ))
418			enddo
419			call histwrite(fileid, 'psi'//nom(iQ), itau, psiQ(:, 1:llm, iQ))
420			enddo
421	guez	3
422	guez	40	call histwrite(fileid, 'masse', itau, zmasse)
423			call histwrite(fileid, 'v', itau, zv)
424			psi=psi*1.e-9
425			call histwrite(fileid, 'psi', itau, psi(:, 1:llm))
426			endif
427	guez	3
428	guez	40	! Moyenne verticale
429	guez	3
430	guez	40	zamasse=0.
431			do l=1, llm
432			zamasse(:)=zamasse(:)+zmasse(:, l)
433			enddo
434			zavQ=0.
435			do iQ=1, nQ
436			do itr=2, ntr
437			do l=1, llm
438			zavQ(:, itr, iQ) = zavQ(:, itr, iQ) &
439			+ zvQ(:, l, itr, iQ) * zmasse(:, l)
440			enddo
441			zavQ(:, itr, iQ)=zavQ(:, itr, iQ)/zamasse(:)
442			call histwrite(fileid, 'a'//znom(itr, iQ), itau, zavQ(:, itr, iQ))
443			enddo
444			enddo
445	guez	3
446	guez	40	! on doit pouvoir tracer systematiquement la fonction de courant.
447			icum=0
448			endif writing_step
449	guez	3
450	guez	40	end SUBROUTINE bilan_dyn
451	guez	3
452	guez	40	end module bilan_dyn_m