GCC Code Coverage Report


Directory: ./
File: phys/diagphy.f90
Date: 2022-01-11 19:19:34
Exec Total Coverage
Lines: 0 122 0.0%
Branches: 0 26 0.0%

Line Branch Exec Source
1
2 ! $Header$
3
4 SUBROUTINE diagphy(airephy, tit, iprt, tops, topl, sols, soll, sens, evap, &
5 rain_fall, snow_fall, ts, d_etp_tot, d_qt_tot, d_ec_tot, fs_bound, &
6 fq_bound)
7 ! ======================================================================
8
9 ! Purpose:
10 ! Compute the thermal flux and the watter mass flux at the atmosphere
11 ! boundaries. Print them and also the atmospheric enthalpy change and
12 ! the atmospheric mass change.
13
14 ! Arguments:
15 ! airephy-------input-R- grid area
16 ! tit---------input-A15- Comment to be added in PRINT (CHARACTER*15)
17 ! iprt--------input-I- PRINT level ( <=0 : no PRINT)
18 ! tops(klon)--input-R- SW rad. at TOA (W/m2), positive up.
19 ! topl(klon)--input-R- LW rad. at TOA (W/m2), positive down
20 ! sols(klon)--input-R- Net SW flux above surface (W/m2), positive up
21 ! (i.e. -1 * flux absorbed by the surface)
22 ! soll(klon)--input-R- Net LW flux above surface (W/m2), positive up
23 ! (i.e. flux emited - flux absorbed by the surface)
24 ! sens(klon)--input-R- Sensible Flux at surface (W/m2), positive down
25 ! evap(klon)--input-R- Evaporation + sublimation watter vapour mass flux
26 ! (kg/m2/s), positive up
27 ! rain_fall(klon)
28 ! --input-R- Liquid watter mass flux (kg/m2/s), positive down
29 ! snow_fall(klon)
30 ! --input-R- Solid watter mass flux (kg/m2/s), positive down
31 ! ts(klon)----input-R- Surface temperature (K)
32 ! d_etp_tot---input-R- Heat flux equivalent to atmospheric enthalpy
33 ! change (W/m2)
34 ! d_qt_tot----input-R- Mass flux equivalent to atmospheric watter mass
35 ! change (kg/m2/s)
36 ! d_ec_tot----input-R- Flux equivalent to atmospheric cinetic energy
37 ! change (W/m2)
38
39 ! fs_bound---output-R- Thermal flux at the atmosphere boundaries (W/m2)
40 ! fq_bound---output-R- Watter mass flux at the atmosphere boundaries
41 ! (kg/m2/s)
42
43 ! J.L. Dufresne, July 2002
44 ! Version prise sur
45 ! ~rlmd833/LMDZOR_201102/modipsl/modeles/LMDZ.3.3/libf/phylmd
46 ! le 25 Novembre 2002.
47 ! ======================================================================
48
49 USE dimphy
50 IMPLICIT NONE
51
52 include "YOMCST.h"
53 include "YOETHF.h"
54
55 ! Input variables
56 REAL airephy(klon)
57 CHARACTER *15 tit
58 INTEGER iprt
59 REAL tops(klon), topl(klon), sols(klon), soll(klon)
60 REAL sens(klon), evap(klon), rain_fall(klon), snow_fall(klon)
61 REAL ts(klon)
62 REAL d_etp_tot, d_qt_tot, d_ec_tot
63 ! Output variables
64 REAL fs_bound, fq_bound
65
66 ! Local variables
67 REAL stops, stopl, ssols, ssoll
68 REAL ssens, sfront, slat
69 REAL airetot, zcpvap, zcwat, zcice
70 REAL rain_fall_tot, snow_fall_tot, evap_tot
71
72 INTEGER i
73
74 INTEGER pas
75 SAVE pas
76 DATA pas/0/
77 !$OMP THREADPRIVATE(pas)
78
79 pas = pas + 1
80 stops = 0.
81 stopl = 0.
82 ssols = 0.
83 ssoll = 0.
84 ssens = 0.
85 sfront = 0.
86 evap_tot = 0.
87 rain_fall_tot = 0.
88 snow_fall_tot = 0.
89 airetot = 0.
90
91 ! Pour les chaleur specifiques de la vapeur d'eau, de l'eau et de
92 ! la glace, on travaille par difference a la chaleur specifique de l'
93 ! air sec. En effet, comme on travaille a niveau de pression donne,
94 ! toute variation de la masse d'un constituant est totalement
95 ! compense par une variation de masse d'air.
96
97 zcpvap = rcpv - rcpd
98 zcwat = rcw - rcpd
99 zcice = rcs - rcpd
100
101 DO i = 1, klon
102 stops = stops + tops(i)*airephy(i)
103 stopl = stopl + topl(i)*airephy(i)
104 ssols = ssols + sols(i)*airephy(i)
105 ssoll = ssoll + soll(i)*airephy(i)
106 ssens = ssens + sens(i)*airephy(i)
107 sfront = sfront + (evap(i)*zcpvap-rain_fall(i)*zcwat-snow_fall(i)*zcice)* &
108 ts(i)*airephy(i)
109 evap_tot = evap_tot + evap(i)*airephy(i)
110 rain_fall_tot = rain_fall_tot + rain_fall(i)*airephy(i)
111 snow_fall_tot = snow_fall_tot + snow_fall(i)*airephy(i)
112 airetot = airetot + airephy(i)
113 END DO
114 stops = stops/airetot
115 stopl = stopl/airetot
116 ssols = ssols/airetot
117 ssoll = ssoll/airetot
118 ssens = ssens/airetot
119 sfront = sfront/airetot
120 evap_tot = evap_tot/airetot
121 rain_fall_tot = rain_fall_tot/airetot
122 snow_fall_tot = snow_fall_tot/airetot
123
124 slat = rlvtt*rain_fall_tot + rlstt*snow_fall_tot
125 ! Heat flux at atm. boundaries
126 fs_bound = stops - stopl - (ssols+ssoll) + ssens + sfront + slat
127 ! Watter flux at atm. boundaries
128 fq_bound = evap_tot - rain_fall_tot - snow_fall_tot
129
130 IF (iprt>=1) WRITE (6, 6666) tit, pas, fs_bound, d_etp_tot, fq_bound, &
131 d_qt_tot
132
133 IF (iprt>=1) WRITE (6, 6668) tit, pas, d_etp_tot + d_ec_tot - fs_bound, &
134 d_qt_tot - fq_bound
135
136 IF (iprt>=2) WRITE (6, 6667) tit, pas, stops, stopl, ssols, ssoll, ssens, &
137 slat, evap_tot, rain_fall_tot + snow_fall_tot
138
139 RETURN
140
141 6666 FORMAT ('Phys. Flux Budget ', A15, 1I6, 2F8.2, 2(1PE13.5))
142 6667 FORMAT ('Phys. Boundary Flux ', A15, 1I6, 6F8.2, 2(1PE13.5))
143 6668 FORMAT ('Phys. Total Budget ', A15, 1I6, F8.2, 2(1PE13.5))
144
145 END SUBROUTINE diagphy
146
147 ! ======================================================================
148 SUBROUTINE diagetpq(airephy, tit, iprt, idiag, idiag2, dtime, t, q, ql, qs, &
149 u, v, paprs, pplay, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
150 ! ======================================================================
151
152 ! Purpose:
153 ! Calcul la difference d'enthalpie et de masse d'eau entre 2 appels,
154 ! et calcul le flux de chaleur et le flux d'eau necessaire a ces
155 ! changements. Ces valeurs sont moyennees sur la surface de tout
156 ! le globe et sont exprime en W/2 et kg/s/m2
157 ! Outil pour diagnostiquer la conservation de l'energie
158 ! et de la masse dans la physique. Suppose que les niveau de
159 ! pression entre couche ne varie pas entre 2 appels.
160
161 ! Plusieurs de ces diagnostics peuvent etre fait en parallele: les
162 ! bilans sont sauvegardes dans des tableaux indices. On parlera
163 ! "d'indice de diagnostic"
164
165
166 ! ======================================================================
167 ! Arguments:
168 ! airephy-------input-R- grid area
169 ! tit-----imput-A15- Comment added in PRINT (CHARACTER*15)
170 ! iprt----input-I- PRINT level ( <=1 : no PRINT)
171 ! idiag---input-I- indice dans lequel sera range les nouveaux
172 ! bilans d' entalpie et de masse
173 ! idiag2--input-I-les nouveaux bilans d'entalpie et de masse
174 ! sont compare au bilan de d'enthalpie de masse de
175 ! l'indice numero idiag2
176 ! Cas parriculier : si idiag2=0, pas de comparaison, on
177 ! sort directement les bilans d'enthalpie et de masse
178 ! dtime----input-R- time step (s)
179 ! t--------input-R- temperature (K)
180 ! q--------input-R- vapeur d'eau (kg/kg)
181 ! ql-------input-R- liquid watter (kg/kg)
182 ! qs-------input-R- solid watter (kg/kg)
183 ! u--------input-R- vitesse u
184 ! v--------input-R- vitesse v
185 ! paprs----input-R- pression a intercouche (Pa)
186 ! pplay----input-R- pression au milieu de couche (Pa)
187
188 ! the following total value are computed by UNIT of earth surface
189
190 ! d_h_vcol--output-R- Heat flux (W/m2) define as the Enthalpy
191 ! change (J/m2) during one time step (dtime) for the whole
192 ! atmosphere (air, watter vapour, liquid and solid)
193 ! d_qt------output-R- total water mass flux (kg/m2/s) defined as the
194 ! total watter (kg/m2) change during one time step (dtime),
195 ! d_qw------output-R- same, for the watter vapour only (kg/m2/s)
196 ! d_ql------output-R- same, for the liquid watter only (kg/m2/s)
197 ! d_qs------output-R- same, for the solid watter only (kg/m2/s)
198 ! d_ec------output-R- Cinetic Energy Budget (W/m2) for vertical air column
199
200 ! other (COMMON...)
201 ! RCPD, RCPV, ....
202
203 ! J.L. Dufresne, July 2002
204 ! ======================================================================
205
206 USE dimphy
207 IMPLICIT NONE
208
209 include "YOMCST.h"
210 include "YOETHF.h"
211
212 ! Input variables
213 REAL airephy(klon)
214 CHARACTER *15 tit
215 INTEGER iprt, idiag, idiag2
216 REAL dtime
217 REAL t(klon, klev), q(klon, klev), ql(klon, klev), qs(klon, klev)
218 REAL u(klon, klev), v(klon, klev)
219 REAL paprs(klon, klev+1), pplay(klon, klev)
220 ! Output variables
221 REAL d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec
222
223 ! Local variables
224
225 REAL h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot, h_qs_tot, qw_tot, ql_tot, &
226 qs_tot, ec_tot
227 ! h_vcol_tot-- total enthalpy of vertical air column
228 ! (air with watter vapour, liquid and solid) (J/m2)
229 ! h_dair_tot-- total enthalpy of dry air (J/m2)
230 ! h_qw_tot---- total enthalpy of watter vapour (J/m2)
231 ! h_ql_tot---- total enthalpy of liquid watter (J/m2)
232 ! h_qs_tot---- total enthalpy of solid watter (J/m2)
233 ! qw_tot------ total mass of watter vapour (kg/m2)
234 ! ql_tot------ total mass of liquid watter (kg/m2)
235 ! qs_tot------ total mass of solid watter (kg/m2)
236 ! ec_tot------ total cinetic energy (kg/m2)
237
238 REAL zairm(klon, klev) ! layer air mass (kg/m2)
239 REAL zqw_col(klon)
240 REAL zql_col(klon)
241 REAL zqs_col(klon)
242 REAL zec_col(klon)
243 REAL zh_dair_col(klon)
244 REAL zh_qw_col(klon), zh_ql_col(klon), zh_qs_col(klon)
245
246 REAL d_h_dair, d_h_qw, d_h_ql, d_h_qs
247
248 REAL airetot, zcpvap, zcwat, zcice
249
250 INTEGER i, k
251
252 INTEGER ndiag ! max number of diagnostic in parallel
253 PARAMETER (ndiag=10)
254 INTEGER pas(ndiag)
255 SAVE pas
256 DATA pas/ndiag*0/
257 !$OMP THREADPRIVATE(pas)
258
259 REAL h_vcol_pre(ndiag), h_dair_pre(ndiag), h_qw_pre(ndiag), &
260 h_ql_pre(ndiag), h_qs_pre(ndiag), qw_pre(ndiag), ql_pre(ndiag), &
261 qs_pre(ndiag), ec_pre(ndiag)
262 SAVE h_vcol_pre, h_dair_pre, h_qw_pre, h_ql_pre, h_qs_pre, qw_pre, ql_pre, &
263 qs_pre, ec_pre
264 !$OMP THREADPRIVATE(h_vcol_pre, h_dair_pre, h_qw_pre, h_ql_pre)
265 !$OMP THREADPRIVATE(h_qs_pre, qw_pre, ql_pre, qs_pre , ec_pre)
266 ! ======================================================================
267
268 DO k = 1, klev
269 DO i = 1, klon
270 ! layer air mass
271 zairm(i, k) = (paprs(i,k)-paprs(i,k+1))/rg
272 END DO
273 END DO
274
275 ! Reset variables
276 DO i = 1, klon
277 zqw_col(i) = 0.
278 zql_col(i) = 0.
279 zqs_col(i) = 0.
280 zec_col(i) = 0.
281 zh_dair_col(i) = 0.
282 zh_qw_col(i) = 0.
283 zh_ql_col(i) = 0.
284 zh_qs_col(i) = 0.
285 END DO
286
287 zcpvap = rcpv
288 zcwat = rcw
289 zcice = rcs
290
291 ! Compute vertical sum for each atmospheric column
292 ! ================================================
293 DO k = 1, klev
294 DO i = 1, klon
295 ! Watter mass
296 zqw_col(i) = zqw_col(i) + q(i, k)*zairm(i, k)
297 zql_col(i) = zql_col(i) + ql(i, k)*zairm(i, k)
298 zqs_col(i) = zqs_col(i) + qs(i, k)*zairm(i, k)
299 ! Cinetic Energy
300 zec_col(i) = zec_col(i) + 0.5*(u(i,k)**2+v(i,k)**2)*zairm(i, k)
301 ! Air enthalpy
302 zh_dair_col(i) = zh_dair_col(i) + rcpd*(1.-q(i,k)-ql(i,k)-qs(i,k))* &
303 zairm(i, k)*t(i, k)
304 zh_qw_col(i) = zh_qw_col(i) + zcpvap*q(i, k)*zairm(i, k)*t(i, k)
305 zh_ql_col(i) = zh_ql_col(i) + zcwat*ql(i, k)*zairm(i, k)*t(i, k) - &
306 rlvtt*ql(i, k)*zairm(i, k)
307 zh_qs_col(i) = zh_qs_col(i) + zcice*qs(i, k)*zairm(i, k)*t(i, k) - &
308 rlstt*qs(i, k)*zairm(i, k)
309
310 END DO
311 END DO
312
313 ! Mean over the planete surface
314 ! =============================
315 qw_tot = 0.
316 ql_tot = 0.
317 qs_tot = 0.
318 ec_tot = 0.
319 h_vcol_tot = 0.
320 h_dair_tot = 0.
321 h_qw_tot = 0.
322 h_ql_tot = 0.
323 h_qs_tot = 0.
324 airetot = 0.
325
326 DO i = 1, klon
327 qw_tot = qw_tot + zqw_col(i)*airephy(i)
328 ql_tot = ql_tot + zql_col(i)*airephy(i)
329 qs_tot = qs_tot + zqs_col(i)*airephy(i)
330 ec_tot = ec_tot + zec_col(i)*airephy(i)
331 h_dair_tot = h_dair_tot + zh_dair_col(i)*airephy(i)
332 h_qw_tot = h_qw_tot + zh_qw_col(i)*airephy(i)
333 h_ql_tot = h_ql_tot + zh_ql_col(i)*airephy(i)
334 h_qs_tot = h_qs_tot + zh_qs_col(i)*airephy(i)
335 airetot = airetot + airephy(i)
336 END DO
337
338 qw_tot = qw_tot/airetot
339 ql_tot = ql_tot/airetot
340 qs_tot = qs_tot/airetot
341 ec_tot = ec_tot/airetot
342 h_dair_tot = h_dair_tot/airetot
343 h_qw_tot = h_qw_tot/airetot
344 h_ql_tot = h_ql_tot/airetot
345 h_qs_tot = h_qs_tot/airetot
346
347 h_vcol_tot = h_dair_tot + h_qw_tot + h_ql_tot + h_qs_tot
348
349 ! Compute the change of the atmospheric state compare to the one
350 ! stored in "idiag2", and convert it in flux. THis computation
351 ! is performed IF idiag2 /= 0 and IF it is not the first CALL
352 ! for "idiag"
353 ! ===================================
354
355 IF ((idiag2>0) .AND. (pas(idiag2)/=0)) THEN
356 d_h_vcol = (h_vcol_tot-h_vcol_pre(idiag2))/dtime
357 d_h_dair = (h_dair_tot-h_dair_pre(idiag2))/dtime
358 d_h_qw = (h_qw_tot-h_qw_pre(idiag2))/dtime
359 d_h_ql = (h_ql_tot-h_ql_pre(idiag2))/dtime
360 d_h_qs = (h_qs_tot-h_qs_pre(idiag2))/dtime
361 d_qw = (qw_tot-qw_pre(idiag2))/dtime
362 d_ql = (ql_tot-ql_pre(idiag2))/dtime
363 d_qs = (qs_tot-qs_pre(idiag2))/dtime
364 d_ec = (ec_tot-ec_pre(idiag2))/dtime
365 d_qt = d_qw + d_ql + d_qs
366 ELSE
367 d_h_vcol = 0.
368 d_h_dair = 0.
369 d_h_qw = 0.
370 d_h_ql = 0.
371 d_h_qs = 0.
372 d_qw = 0.
373 d_ql = 0.
374 d_qs = 0.
375 d_ec = 0.
376 d_qt = 0.
377 END IF
378
379 IF (iprt>=2) THEN
380 WRITE (6, 9000) tit, pas(idiag), d_qt, d_qw, d_ql, d_qs
381 9000 FORMAT ('Phys. Watter Mass Budget (kg/m2/s)', A15, 1I6, 10(1PE14.6))
382 WRITE (6, 9001) tit, pas(idiag), d_h_vcol
383 9001 FORMAT ('Phys. Enthalpy Budget (W/m2) ', A15, 1I6, 10(F8.2))
384 WRITE (6, 9002) tit, pas(idiag), d_ec
385 9002 FORMAT ('Phys. Cinetic Energy Budget (W/m2) ', A15, 1I6, 10(F8.2))
386 END IF
387
388 ! Store the new atmospheric state in "idiag"
389
390 pas(idiag) = pas(idiag) + 1
391 h_vcol_pre(idiag) = h_vcol_tot
392 h_dair_pre(idiag) = h_dair_tot
393 h_qw_pre(idiag) = h_qw_tot
394 h_ql_pre(idiag) = h_ql_tot
395 h_qs_pre(idiag) = h_qs_tot
396 qw_pre(idiag) = qw_tot
397 ql_pre(idiag) = ql_tot
398 qs_pre(idiag) = qs_tot
399 ec_pre(idiag) = ec_tot
400
401 RETURN
402 END SUBROUTINE diagetpq
403