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