libf/dyn3d/bilan_dyn.f90

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	module bilan_dyn_m
2
3	IMPLICIT NONE
4
5	contains
6
7	SUBROUTINE bilan_dyn(ps, masse, pk, flux_u, flux_v, teta, phi, ucov, vcov, &
8	trac)
9
10	! From LMDZ4/libf/dyn3d/bilan_dyn.F, version 1.5 2005/03/16 10:12:17
11
12	! 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
16	USE comconst, ONLY: cpp
17	USE comgeom, ONLY: constang_2d, cu_2d, cv_2d
18	USE dimens_m, ONLY: iim, jjm, llm
19	USE histwrite_m, ONLY: histwrite
20	use init_dynzon_m, only: ncum, fileid, znom, ntr, nq, nom
21	USE paramet_m, ONLY: iip1, jjp1
22
23	! Arguments:
24
25	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	real, intent(in):: teta(iip1, jjp1, llm)
30	real, intent(in):: phi(iip1, jjp1, llm)
31	real, intent(in):: ucov(iip1, jjp1, llm)
32	real, intent(in):: vcov(iip1, jjm, llm)
33	real, intent(in):: trac(:, :, :) ! (iim + 1, jjm + 1, llm)
34
35	! Local:
36
37	integer:: icum = 0
38	integer:: itau = 0
39	real zqy, zfactv(jjm, llm)
40
41	real ww
42
43	! 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
49	! Champ contenant les scalaires advectés
50	real Q(iip1, jjp1, llm, nQ)
51
52	! 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
62	! champs de tansport en moyenne zonale
63	integer itr
64	integer, parameter:: iave = 1, itot = 2, immc = 3, itrs = 4, istn = 5
65
66	real zvQ(jjm, llm, ntr, nQ), zvQtmp(jjm, llm)
67	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	integer i, j, l, iQ
71
72	!-----------------------------------------------------------------
73
74	! Calcul des champs dynamiques
75
76	! Énergie cinétique
77	ucont = 0
78	CALL covcont(llm, ucov, vcov, ucont, vcont)
79	CALL enercin(vcov, ucov, vcont, ucont, ecin)
80
81	! moment cinétique
82	do l = 1, llm
83	ang(:, :, l) = ucov(:, :, l) + constang_2d
84	unat(:, :, l) = ucont(:, :, l)*cu_2d
85	enddo
86
87	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
95	! Cumul
96
97	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	endif
106
107	itau = itau + 1
108	icum = icum + 1
109
110	! Accumulation des flux de masse horizontaux
111	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	enddo
118
119	! FLUX ET TENDANCES
120
121	! Flux longitudinal
122	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
127	! 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
132	! tendances
133
134	! convergence horizontale
135	call convflu(flux_uQ_cum, flux_vQ_cum, llm*nQ, dQ)
136
137	! calcul de la vitesse verticale
138	call convmas(flux_u_cum, flux_v_cum, convm)
139	CALL vitvert(convm, w)
140
141	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	enddo
149	enddo
150	enddo
151	enddo
152
153	! PAS DE TEMPS D'ECRITURE
154
155	writing_step: if (icum == ncum) then
156	! Normalisation
157	do iQ = 1, nQ
158	Q_cum(:, :, :, iQ) = Q_cum(:, :, :, iQ)/masse_cum
159	enddo
160	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
168	! A retravailler eventuellement
169	! division de dQ par la masse pour revenir aux bonnes grandeurs
170	do iQ = 1, nQ
171	dQ(:, :, :, iQ) = dQ(:, :, :, iQ)/masse_cum
172	enddo
173
174	! Transport méridien
175
176	! cumul zonal des masses des mailles
177
178	zv = 0.
179	zmasse = 0.
180	call massbar(masse_cum, massebx, masseby)
181	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	enddo
187	zfactv(j, l) = cv_2d(1, j)/zmasse(j, l)
188	enddo
189	enddo
190
191	! Transport dans le plan latitude-altitude
192
193	zvQ = 0.
194	psiQ = 0.
195	do iQ = 1, nQ
196	zvQtmp = 0.
197	do l = 1, llm
198	do j = 1, jjm
199	! Calcul des moyennes zonales du transort total et de zvQtmp
200	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	enddo
209	! Decomposition
210	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	enddo
217	enddo
218	! fonction de courant meridienne pour la quantite Q
219	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	enddo
223	enddo
224	enddo
225
226	! fonction de courant pour la circulation meridienne moyenne
227	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	enddo
233	enddo
234
235	! sorties proprement dites
236	do iQ = 1, nQ
237	do itr = 1, ntr
238	call histwrite(fileid, znom(itr, iQ), itau, zvQ(:, :, itr, iQ))
239	enddo
240	call histwrite(fileid, 'psi'//nom(iQ), itau, psiQ(:, :llm, iQ))
241	enddo
242
243	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
248	! Intégrale verticale
249
250	forall (iQ = 1: nQ, itr = 2: ntr) zavQ(:, itr, iQ) &
251	= sum(zvQ(:, :, itr, iQ) * zmasse, dim=2) / cv_2d(1, :)
252
253	do iQ = 1, nQ
254	do itr = 2, ntr
255	call histwrite(fileid, 'a'//znom(itr, iQ), itau, zavQ(:, itr, iQ))
256	enddo
257	enddo
258
259	! On doit pouvoir tracer systematiquement la fonction de courant.
260	icum = 0
261	endif writing_step
262
263	end SUBROUTINE bilan_dyn
264
265	end module bilan_dyn_m