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 massbar_m, only: massbar
    USE paramet_m, ONLY: iip1, jjp1

    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 qy, factv(jjm, llm)

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

    ! champs de tansport en moyenne zonale
    integer itr
    integer, parameter:: iave = 1, itot = 2, immc = 3, itrs = 4, istn = 5

    real vq(jjm, llm, ntr, nQ), vqtmp(jjm, llm)
    real avq(jjm, 2: ntr, nQ), psiQ(jjm, llm + 1, nQ)
    real zmasse(jjm, llm)
    real v(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
    forall (l = 1: llm)
       ang(:, :, l) = ucov(:, :, l) + constang_2d
       unat(:, :, l) = ucont(:, :, l) * cu_2d
    end forall

    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
    forall (iQ = 1: nQ) Q_cum(:, :, :, iQ) = Q_cum(:, :, :, iQ) &
         + Q(:, :, :, iQ) * masse

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

    writing_step: if (icum == ncum) then
       ! Normalisation
       forall (iQ = 1: nQ) Q_cum(:, :, :, iQ) = Q_cum(:, :, :, iQ) / masse_cum
       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

       ! Transport méridien

       ! Cumul zonal des masses des mailles

       v = 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)
                v(j, l) = v(j, l) + flux_v_cum(i, j, l)
             enddo
             factv(j, l) = cv_2d(1, j) / zmasse(j, l)
          enddo
       enddo

       ! Transport dans le plan latitude-altitude

       vq = 0.
       psiQ = 0.
       do iQ = 1, nQ
          vqtmp = 0.
          do l = 1, llm
             do j = 1, jjm
                ! Calcul des moyennes zonales du transport total et de vqtmp
                do i = 1, iim
                   vq(j, l, itot, iQ) = vq(j, l, itot, iQ) &
                        + flux_vQ_cum(i, j, l, iQ)
                   qy =  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))
                   vqtmp(j, l) = vqtmp(j, l) + flux_v_cum(i, j, l) * qy &
                        / (0.5 * (masse_cum(i, j, l) + masse_cum(i, j + 1, l)))
                   vq(j, l, iave, iQ) = vq(j, l, iave, iQ) + qy
                enddo
                ! Decomposition
                vq(j, l, iave, iQ) = vq(j, l, iave, iQ) / zmasse(j, l)
                vq(j, l, itot, iQ) = vq(j, l, itot, iQ) * factv(j, l)
                vqtmp(j, l) = vqtmp(j, l) * factv(j, l)
                vq(j, l, immc, iQ) = v(j, l) * vq(j, l, iave, iQ) * factv(j, l)
                vq(j, l, itrs, iQ) = vq(j, l, itot, iQ) - vqtmp(j, l)
                vq(j, l, istn, iQ) = vqtmp(j, l) - vq(j, l, immc, iQ)
             enddo
          enddo
          ! Fonction de courant méridienne pour la quantité Q
          do l = llm, 1, -1
             do j = 1, jjm
                psiQ(j, l, iQ) = psiQ(j, l + 1, iQ) + vq(j, l, itot, iQ)
             enddo
          enddo
       enddo

       ! Fonction de courant pour la circulation méridienne moyenne
       psi = 0.
       do l = llm, 1, -1
          do j = 1, jjm
             psi(j, l) = psi(j, l + 1) + v(j, l)
             v(j, l) = v(j, l) * factv(j, l)
          enddo
       enddo

       ! Sorties proprement dites
       do iQ = 1, nQ
          do itr = 1, ntr
             call histwrite(fileid, znom(itr, iQ), itau, vq(:, :, 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, v)
       psi = psi * 1e-9
       call histwrite(fileid, 'psi', itau, psi(:, :llm))

       ! Intégrale verticale

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

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

       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	91	use massbar_m, only: massbar
22	guez	56	USE paramet_m, ONLY: iip1, jjp1
23	guez	3
24	guez	56	real, intent(in):: ps(iip1, jjp1)
25			real, intent(in):: masse(iip1, jjp1, llm), pk(iip1, jjp1, llm)
26			real, intent(in):: flux_u(iip1, jjp1, llm)
27			real, intent(in):: flux_v(iip1, jjm, llm)
28	guez	44	real, intent(in):: teta(iip1, jjp1, llm)
29	guez	56	real, intent(in):: phi(iip1, jjp1, llm)
30	guez	62	real, intent(in):: ucov(:, :, :) ! (iip1, jjp1, llm)
31	guez	56	real, intent(in):: vcov(iip1, jjm, llm)
32	guez	40	real, intent(in):: trac(:, :, :) ! (iim + 1, jjm + 1, llm)
33	guez	3
34	guez	40	! Local:
35	guez	3
36	guez	54	integer:: icum = 0
37	guez	57	integer:: itau = 0
38	guez	62	real qy, factv(jjm, llm)
39	guez	3
40	guez	40	! Variables dynamiques intermédiaires
41			REAL vcont(iip1, jjm, llm), ucont(iip1, jjp1, llm)
42			REAL ang(iip1, jjp1, llm), unat(iip1, jjp1, llm)
43			REAL massebx(iip1, jjp1, llm), masseby(iip1, jjm, llm)
44	guez	62	REAL ecin(iip1, jjp1, llm)
45	guez	3
46	guez	40	! Champ contenant les scalaires advectés
47			real Q(iip1, jjp1, llm, nQ)
48	guez	3
49	guez	40	! Champs cumulés
50			real, save:: ps_cum(iip1, jjp1)
51			real, save:: masse_cum(iip1, jjp1, llm)
52			real, save:: flux_u_cum(iip1, jjp1, llm)
53			real, save:: flux_v_cum(iip1, jjm, llm)
54			real, save:: Q_cum(iip1, jjp1, llm, nQ)
55			real, save:: flux_uQ_cum(iip1, jjp1, llm, nQ)
56			real, save:: flux_vQ_cum(iip1, jjm, llm, nQ)
57	guez	3
58	guez	40	! champs de tansport en moyenne zonale
59			integer itr
60	guez	54	integer, parameter:: iave = 1, itot = 2, immc = 3, itrs = 4, istn = 5
61	guez	3
62	guez	62	real vq(jjm, llm, ntr, nQ), vqtmp(jjm, llm)
63			real avq(jjm, 2: ntr, nQ), psiQ(jjm, llm + 1, nQ)
64	guez	54	real zmasse(jjm, llm)
65	guez	62	real v(jjm, llm), psi(jjm, llm + 1)
66	guez	40	integer i, j, l, iQ
67	guez	3
68	guez	40	!-----------------------------------------------------------------
69	guez	3
70	guez	40	! Calcul des champs dynamiques
71	guez	3
72	guez	40	! Énergie cinétique
73			ucont = 0
74			CALL covcont(llm, ucov, vcov, ucont, vcont)
75			CALL enercin(vcov, ucov, vcont, ucont, ecin)
76	guez	3
77	guez	40	! moment cinétique
78	guez	62	forall (l = 1: llm)
79	guez	54	ang(:, :, l) = ucov(:, :, l) + constang_2d
80	guez	62	unat(:, :, l) = ucont(:, :, l) * cu_2d
81			end forall
82	guez	3
83	guez	54	Q(:, :, :, 1) = teta * pk / cpp
84			Q(:, :, :, 2) = phi
85			Q(:, :, :, 3) = ecin
86			Q(:, :, :, 4) = ang
87			Q(:, :, :, 5) = unat
88			Q(:, :, :, 6) = trac
89			Q(:, :, :, 7) = 1.
90	guez	3
91	guez	40	! Cumul
92	guez	3
93	guez	54	if (icum == 0) then
94			ps_cum = 0.
95			masse_cum = 0.
96			flux_u_cum = 0.
97			flux_v_cum = 0.
98			Q_cum = 0.
99			flux_vQ_cum = 0.
100			flux_uQ_cum = 0.
101	guez	40	endif
102	guez	3
103	guez	57	itau = itau + 1
104	guez	54	icum = icum + 1
105	guez	3
106	guez	40	! Accumulation des flux de masse horizontaux
107	guez	54	ps_cum = ps_cum + ps
108			masse_cum = masse_cum + masse
109			flux_u_cum = flux_u_cum + flux_u
110			flux_v_cum = flux_v_cum + flux_v
111	guez	62	forall (iQ = 1: nQ) Q_cum(:, :, :, iQ) = Q_cum(:, :, :, iQ) &
112			+ Q(:, :, :, iQ) * masse
113	guez	3
114	guez	40	! Flux longitudinal
115	guez	54	forall (iQ = 1: nQ, i = 1: iim) flux_uQ_cum(i, :, :, iQ) &
116			= flux_uQ_cum(i, :, :, iQ) &
117			+ flux_u(i, :, :) * 0.5 * (Q(i, :, :, iQ) + Q(i + 1, :, :, iQ))
118			flux_uQ_cum(iip1, :, :, :) = flux_uQ_cum(1, :, :, :)
119	guez	3
120	guez	54	! Flux méridien
121			forall (iQ = 1: nQ, j = 1: jjm) flux_vQ_cum(:, j, :, iQ) &
122			= flux_vQ_cum(:, j, :, iQ) &
123			+ flux_v(:, j, :) * 0.5 * (Q(:, j, :, iQ) + Q(:, j + 1, :, iQ))
124	guez	3
125	guez	40	writing_step: if (icum == ncum) then
126			! Normalisation
127	guez	62	forall (iQ = 1: nQ) Q_cum(:, :, :, iQ) = Q_cum(:, :, :, iQ) / masse_cum
128	guez	56	ps_cum = ps_cum / ncum
129			masse_cum = masse_cum / ncum
130			flux_u_cum = flux_u_cum / ncum
131			flux_v_cum = flux_v_cum / ncum
132			flux_uQ_cum = flux_uQ_cum / ncum
133			flux_vQ_cum = flux_vQ_cum / ncum
134	guez	3
135	guez	40	! Transport méridien
136	guez	3
137	guez	62	! Cumul zonal des masses des mailles
138	guez	3
139	guez	62	v = 0.
140	guez	54	zmasse = 0.
141	guez	40	call massbar(masse_cum, massebx, masseby)
142	guez	54	do l = 1, llm
143			do j = 1, jjm
144			do i = 1, iim
145			zmasse(j, l) = zmasse(j, l) + masseby(i, j, l)
146	guez	62	v(j, l) = v(j, l) + flux_v_cum(i, j, l)
147	guez	40	enddo
148	guez	62	factv(j, l) = cv_2d(1, j) / zmasse(j, l)
149	guez	40	enddo
150			enddo
151	guez	3
152	guez	40	! Transport dans le plan latitude-altitude
153	guez	3
154	guez	62	vq = 0.
155	guez	54	psiQ = 0.
156			do iQ = 1, nQ
157	guez	62	vqtmp = 0.
158	guez	54	do l = 1, llm
159			do j = 1, jjm
160	guez	62	! Calcul des moyennes zonales du transport total et de vqtmp
161	guez	54	do i = 1, iim
162	guez	62	vq(j, l, itot, iQ) = vq(j, l, itot, iQ) &
163	guez	54	+ flux_vQ_cum(i, j, l, iQ)
164	guez	62	qy = 0.5 * (Q_cum(i, j, l, iQ) * masse_cum(i, j, l) &
165	guez	54	+ Q_cum(i, j + 1, l, iQ) * masse_cum(i, j + 1, l))
166	guez	62	vqtmp(j, l) = vqtmp(j, l) + flux_v_cum(i, j, l) * qy &
167	guez	54	/ (0.5 * (masse_cum(i, j, l) + masse_cum(i, j + 1, l)))
168	guez	62	vq(j, l, iave, iQ) = vq(j, l, iave, iQ) + qy
169	guez	40	enddo
170			! Decomposition
171	guez	62	vq(j, l, iave, iQ) = vq(j, l, iave, iQ) / zmasse(j, l)
172			vq(j, l, itot, iQ) = vq(j, l, itot, iQ) * factv(j, l)
173			vqtmp(j, l) = vqtmp(j, l) * factv(j, l)
174			vq(j, l, immc, iQ) = v(j, l) * vq(j, l, iave, iQ) * factv(j, l)
175			vq(j, l, itrs, iQ) = vq(j, l, itot, iQ) - vqtmp(j, l)
176			vq(j, l, istn, iQ) = vqtmp(j, l) - vq(j, l, immc, iQ)
177	guez	40	enddo
178			enddo
179	guez	62	! Fonction de courant méridienne pour la quantité Q
180	guez	54	do l = llm, 1, -1
181			do j = 1, jjm
182	guez	62	psiQ(j, l, iQ) = psiQ(j, l + 1, iQ) + vq(j, l, itot, iQ)
183	guez	40	enddo
184			enddo
185			enddo
186	guez	3
187	guez	62	! Fonction de courant pour la circulation méridienne moyenne
188	guez	54	psi = 0.
189			do l = llm, 1, -1
190			do j = 1, jjm
191	guez	62	psi(j, l) = psi(j, l + 1) + v(j, l)
192			v(j, l) = v(j, l) * factv(j, l)
193	guez	40	enddo
194			enddo
195	guez	3
196	guez	62	! Sorties proprement dites
197	guez	54	do iQ = 1, nQ
198			do itr = 1, ntr
199	guez	62	call histwrite(fileid, znom(itr, iQ), itau, vq(:, :, itr, iQ))
200	guez	40	enddo
201	guez	62	call histwrite(fileid, 'psi' // nom(iQ), itau, psiQ(:, :llm, iQ))
202	guez	54	enddo
203	guez	3
204	guez	54	call histwrite(fileid, 'masse', itau, zmasse)
205	guez	62	call histwrite(fileid, 'v', itau, v)
206			psi = psi * 1e-9
207	guez	54	call histwrite(fileid, 'psi', itau, psi(:, :llm))
208	guez	3
209	guez	55	! Intégrale verticale
210	guez	3
211	guez	62	forall (iQ = 1: nQ, itr = 2: ntr) avq(:, itr, iQ) &
212			= sum(vq(:, :, itr, iQ) * zmasse, dim=2) / cv_2d(1, :)
213	guez	54
214			do iQ = 1, nQ
215			do itr = 2, ntr
216	guez	62	call histwrite(fileid, 'a' // znom(itr, iQ), itau, avq(:, itr, iQ))
217	guez	40	enddo
218			enddo
219	guez	3
220	guez	54	icum = 0
221	guez	40	endif writing_step
222	guez	3
223	guez	40	end SUBROUTINE bilan_dyn
224	guez	3
225	guez	40	end module bilan_dyn_m