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)

    ! 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 comconst, ONLY: cpp
    USE comgeom, ONLY: constang_2d, cu_2d, cv_2d
    USE dimens_m, ONLY: iim, jjm, llm
    USE histwrite_m, ONLY: histwrite
    use init_dynzon_m, only: ncum, fileid, znom, ntr, nq, nom
    USE paramet_m, ONLY: iip1, jjp1

    ! 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)

    ! Local:

    integer:: icum  = 0
    integer:: itau = 0
    real zqy, zfactv(jjm, llm)

    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:: iave = 1, itot = 2, immc = 3, itrs = 4, istn = 5

    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

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

    ! 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

    itau = itau + 1
    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	57	trac)
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 comconst, ONLY: cpp
17	guez	57	USE comgeom, ONLY: constang_2d, cu_2d, cv_2d
18	guez	56	USE dimens_m, ONLY: iim, jjm, llm
19			USE histwrite_m, ONLY: histwrite
20	guez	57	use init_dynzon_m, only: ncum, fileid, znom, ntr, nq, nom
21	guez	56	USE paramet_m, ONLY: iip1, jjp1
22	guez	3
23	guez	40	! Arguments:
24	guez	3
25	guez	56	real, intent(in):: ps(iip1, jjp1)
26			real, intent(in):: masse(iip1, jjp1, llm), pk(iip1, jjp1, llm)
27			real, intent(in):: flux_u(iip1, jjp1, llm)
28			real, intent(in):: flux_v(iip1, jjm, llm)
29	guez	44	real, intent(in):: teta(iip1, jjp1, llm)
30	guez	56	real, intent(in):: phi(iip1, jjp1, llm)
31			real, intent(in):: ucov(iip1, jjp1, llm)
32			real, intent(in):: vcov(iip1, jjm, llm)
33	guez	40	real, intent(in):: trac(:, :, :) ! (iim + 1, jjm + 1, llm)
34	guez	3
35	guez	40	! Local:
36	guez	3
37	guez	54	integer:: icum = 0
38	guez	57	integer:: itau = 0
39	guez	56	real zqy, zfactv(jjm, llm)
40	guez	3
41	guez	40	real ww
42	guez	3
43	guez	40	! Variables dynamiques intermédiaires
44			REAL vcont(iip1, jjm, llm), ucont(iip1, jjp1, llm)
45			REAL ang(iip1, jjp1, llm), unat(iip1, jjp1, llm)
46			REAL massebx(iip1, jjp1, llm), masseby(iip1, jjm, llm)
47			REAL w(iip1, jjp1, llm), ecin(iip1, jjp1, llm), convm(iip1, jjp1, llm)
48	guez	3
49	guez	40	! Champ contenant les scalaires advectés
50			real Q(iip1, jjp1, llm, nQ)
51	guez	3
52	guez	40	! Champs cumulés
53			real, save:: ps_cum(iip1, jjp1)
54			real, save:: masse_cum(iip1, jjp1, llm)
55			real, save:: flux_u_cum(iip1, jjp1, llm)
56			real, save:: flux_v_cum(iip1, jjm, llm)
57			real, save:: Q_cum(iip1, jjp1, llm, nQ)
58			real, save:: flux_uQ_cum(iip1, jjp1, llm, nQ)
59			real, save:: flux_vQ_cum(iip1, jjm, llm, nQ)
60			real dQ(iip1, jjp1, llm, nQ)
61	guez	3
62	guez	40	! champs de tansport en moyenne zonale
63			integer itr
64	guez	54	integer, parameter:: iave = 1, itot = 2, immc = 3, itrs = 4, istn = 5
65	guez	3
66	guez	40	real zvQ(jjm, llm, ntr, nQ), zvQtmp(jjm, llm)
67	guez	54	real zavQ(jjm, 2: ntr, nQ), psiQ(jjm, llm + 1, nQ)
68			real zmasse(jjm, llm)
69			real zv(jjm, llm), psi(jjm, llm + 1)
70	guez	40	integer i, j, l, iQ
71	guez	3
72	guez	40	!-----------------------------------------------------------------
73	guez	3
74	guez	40	! Calcul des champs dynamiques
75	guez	3
76	guez	40	! Énergie cinétique
77			ucont = 0
78			CALL covcont(llm, ucov, vcov, ucont, vcont)
79			CALL enercin(vcov, ucov, vcont, ucont, ecin)
80	guez	3
81	guez	40	! moment cinétique
82	guez	54	do l = 1, llm
83			ang(:, :, l) = ucov(:, :, l) + constang_2d
84			unat(:, :, l) = ucont(:, :, l)*cu_2d
85	guez	40	enddo
86	guez	3
87	guez	54	Q(:, :, :, 1) = teta * pk / cpp
88			Q(:, :, :, 2) = phi
89			Q(:, :, :, 3) = ecin
90			Q(:, :, :, 4) = ang
91			Q(:, :, :, 5) = unat
92			Q(:, :, :, 6) = trac
93			Q(:, :, :, 7) = 1.
94	guez	3
95	guez	40	! Cumul
96	guez	3
97	guez	54	if (icum == 0) then
98			ps_cum = 0.
99			masse_cum = 0.
100			flux_u_cum = 0.
101			flux_v_cum = 0.
102			Q_cum = 0.
103			flux_vQ_cum = 0.
104			flux_uQ_cum = 0.
105	guez	40	endif
106	guez	3
107	guez	57	itau = itau + 1
108	guez	54	icum = icum + 1
109	guez	3
110	guez	40	! Accumulation des flux de masse horizontaux
111	guez	54	ps_cum = ps_cum + ps
112			masse_cum = masse_cum + masse
113			flux_u_cum = flux_u_cum + flux_u
114			flux_v_cum = flux_v_cum + flux_v
115			do iQ = 1, nQ
116			Q_cum(:, :, :, iQ) = Q_cum(:, :, :, iQ) + Q(:, :, :, iQ)*masse
117	guez	40	enddo
118	guez	3
119	guez	40	! FLUX ET TENDANCES
120	guez	3
121	guez	40	! Flux longitudinal
122	guez	54	forall (iQ = 1: nQ, i = 1: iim) flux_uQ_cum(i, :, :, iQ) &
123			= flux_uQ_cum(i, :, :, iQ) &
124			+ flux_u(i, :, :) * 0.5 * (Q(i, :, :, iQ) + Q(i + 1, :, :, iQ))
125			flux_uQ_cum(iip1, :, :, :) = flux_uQ_cum(1, :, :, :)
126	guez	3
127	guez	54	! Flux méridien
128			forall (iQ = 1: nQ, j = 1: jjm) flux_vQ_cum(:, j, :, iQ) &
129			= flux_vQ_cum(:, j, :, iQ) &
130			+ flux_v(:, j, :) * 0.5 * (Q(:, j, :, iQ) + Q(:, j + 1, :, iQ))
131	guez	3
132	guez	40	! tendances
133	guez	3
134	guez	40	! convergence horizontale
135			call convflu(flux_uQ_cum, flux_vQ_cum, llm*nQ, dQ)
136	guez	3
137	guez	40	! calcul de la vitesse verticale
138			call convmas(flux_u_cum, flux_v_cum, convm)
139			CALL vitvert(convm, w)
140	guez	3
141	guez	54	do iQ = 1, nQ
142			do l = 1, llm-1
143			do j = 1, jjp1
144			do i = 1, iip1
145			ww = -0.5w(i, j, l + 1)(Q(i, j, l, iQ) + Q(i, j, l + 1, iQ))
146			dQ(i, j, l, iQ) = dQ(i, j, l, iQ)-ww
147			dQ(i, j, l + 1, iQ) = dQ(i, j, l + 1, iQ) + ww
148	guez	40	enddo
149			enddo
150			enddo
151			enddo
152	guez	3
153	guez	40	! PAS DE TEMPS D'ECRITURE
154	guez	3
155	guez	40	writing_step: if (icum == ncum) then
156			! Normalisation
157	guez	54	do iQ = 1, nQ
158			Q_cum(:, :, :, iQ) = Q_cum(:, :, :, iQ)/masse_cum
159	guez	40	enddo
160	guez	56	ps_cum = ps_cum / ncum
161			masse_cum = masse_cum / ncum
162			flux_u_cum = flux_u_cum / ncum
163			flux_v_cum = flux_v_cum / ncum
164			flux_uQ_cum = flux_uQ_cum / ncum
165			flux_vQ_cum = flux_vQ_cum / ncum
166			dQ = dQ / ncum
167	guez	3
168	guez	40	! A retravailler eventuellement
169			! division de dQ par la masse pour revenir aux bonnes grandeurs
170	guez	54	do iQ = 1, nQ
171			dQ(:, :, :, iQ) = dQ(:, :, :, iQ)/masse_cum
172	guez	40	enddo
173	guez	3
174	guez	40	! Transport méridien
175	guez	3
176	guez	40	! cumul zonal des masses des mailles
177	guez	3
178	guez	54	zv = 0.
179			zmasse = 0.
180	guez	40	call massbar(masse_cum, massebx, masseby)
181	guez	54	do l = 1, llm
182			do j = 1, jjm
183			do i = 1, iim
184			zmasse(j, l) = zmasse(j, l) + masseby(i, j, l)
185			zv(j, l) = zv(j, l) + flux_v_cum(i, j, l)
186	guez	40	enddo
187	guez	54	zfactv(j, l) = cv_2d(1, j)/zmasse(j, l)
188	guez	40	enddo
189			enddo
190	guez	3
191	guez	40	! Transport dans le plan latitude-altitude
192	guez	3
193	guez	54	zvQ = 0.
194			psiQ = 0.
195			do iQ = 1, nQ
196			zvQtmp = 0.
197			do l = 1, llm
198			do j = 1, jjm
199	guez	40	! Calcul des moyennes zonales du transort total et de zvQtmp
200	guez	54	do i = 1, iim
201			zvQ(j, l, itot, iQ) = zvQ(j, l, itot, iQ) &
202			+ flux_vQ_cum(i, j, l, iQ)
203			zqy = 0.5 * (Q_cum(i, j, l, iQ) * masse_cum(i, j, l) &
204			+ Q_cum(i, j + 1, l, iQ) * masse_cum(i, j + 1, l))
205			zvQtmp(j, l) = zvQtmp(j, l) + flux_v_cum(i, j, l) * zqy &
206			/ (0.5 * (masse_cum(i, j, l) + masse_cum(i, j + 1, l)))
207			zvQ(j, l, iave, iQ) = zvQ(j, l, iave, iQ) + zqy
208	guez	40	enddo
209			! Decomposition
210	guez	54	zvQ(j, l, iave, iQ) = zvQ(j, l, iave, iQ)/zmasse(j, l)
211			zvQ(j, l, itot, iQ) = zvQ(j, l, itot, iQ)*zfactv(j, l)
212			zvQtmp(j, l) = zvQtmp(j, l)*zfactv(j, l)
213			zvQ(j, l, immc, iQ) = zv(j, l)zvQ(j, l, iave, iQ)zfactv(j, l)
214			zvQ(j, l, itrs, iQ) = zvQ(j, l, itot, iQ)-zvQtmp(j, l)
215			zvQ(j, l, istn, iQ) = zvQtmp(j, l)-zvQ(j, l, immc, iQ)
216	guez	40	enddo
217			enddo
218			! fonction de courant meridienne pour la quantite Q
219	guez	54	do l = llm, 1, -1
220			do j = 1, jjm
221			psiQ(j, l, iQ) = psiQ(j, l + 1, iQ) + zvQ(j, l, itot, iQ)
222	guez	40	enddo
223			enddo
224			enddo
225	guez	3
226	guez	40	! fonction de courant pour la circulation meridienne moyenne
227	guez	54	psi = 0.
228			do l = llm, 1, -1
229			do j = 1, jjm
230			psi(j, l) = psi(j, l + 1) + zv(j, l)
231			zv(j, l) = zv(j, l)*zfactv(j, l)
232	guez	40	enddo
233			enddo
234	guez	3
235	guez	40	! sorties proprement dites
236	guez	54	do iQ = 1, nQ
237			do itr = 1, ntr
238			call histwrite(fileid, znom(itr, iQ), itau, zvQ(:, :, itr, iQ))
239	guez	40	enddo
240	guez	54	call histwrite(fileid, 'psi'//nom(iQ), itau, psiQ(:, :llm, iQ))
241			enddo
242	guez	3
243	guez	54	call histwrite(fileid, 'masse', itau, zmasse)
244			call histwrite(fileid, 'v', itau, zv)
245			psi = psi*1.e-9
246			call histwrite(fileid, 'psi', itau, psi(:, :llm))
247	guez	3
248	guez	55	! Intégrale verticale
249	guez	3
250	guez	54	forall (iQ = 1: nQ, itr = 2: ntr) zavQ(:, itr, iQ) &
251	guez	55	= sum(zvQ(:, :, itr, iQ) * zmasse, dim=2) / cv_2d(1, :)
252	guez	54
253			do iQ = 1, nQ
254			do itr = 2, ntr
255	guez	40	call histwrite(fileid, 'a'//znom(itr, iQ), itau, zavQ(:, itr, iQ))
256			enddo
257			enddo
258	guez	3
259	guez	54	! On doit pouvoir tracer systematiquement la fonction de courant.
260			icum = 0
261	guez	40	endif writing_step
262	guez	3
263	guez	40	end SUBROUTINE bilan_dyn
264	guez	3
265	guez	40	end module bilan_dyn_m