Sources/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\'e \`a des diagnostics dynamiques de
    ! base.  De fa\c{}con g\'en\'erale, les moyennes des scalaires Q
    ! sont pond\'er\'ees par la masse. Les flux de masse sont, eux,
    ! simplement moyenn\'es.

    USE comconst, ONLY: cpp
    USE comgeom, ONLY: constang_2d, cu_2d, cv_2d
    use covcont_m, only: covcont
    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\'ediaires
    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\'es
    real Q(iip1, jjp1, llm, nQ)

    ! Champs cumul\'es
    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

    ! \'Energie cin\'etique
    ucont = 0
    CALL covcont(llm, ucov, vcov, ucont, vcont)
    CALL enercin(vcov, ucov, vcont, ucont, ecin)

    ! moment cin\'etique
    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\'eridien
    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\'eridien

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