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

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