GCC Code Coverage Report

Directory: ./
File: phys/add_phys_tend_mod.f90
Date: 2022-01-11 19:19:34
Exec Total Coverage
Lines: 60 291 20.6%
Branches: 54 348 15.5%
Line Branch Exec Source
1 !
2 ! $Id$
3 !
4 !
5 MODULE add_phys_tend_mod
6
7 IMPLICIT NONE
8 ! flag to compute diagnostics to check energy conservation. If fl_ebil==0, no check
9 INTEGER, SAVE :: fl_ebil
10 !$OMP THREADPRIVATE(fl_ebil)
11 ! flag to include modifcations to ensure energy conservation. If fl_cor_ebil==0, no corrections
12 ! Note that with time, all these modifications should be included by default
13 INTEGER, SAVE :: fl_cor_ebil
14 !$OMP THREADPRIVATE(fl_cor_ebil)
15
16 CONTAINS
17
18 ✗ SUBROUTINE add_pbl_tend(zdu, zdv, zdt, zdq, zdql, zdqi, paprs, text,abortphy,flag_inhib_tend, itap)
19 ! ======================================================================
20 ! Ajoute les tendances de couche limite, soit determinees par la
21 ! parametrisation
22 ! physique, soit forcees, aux variables d etat de la dynamique t_seri,
23 ! q_seri ...
24 ! ======================================================================
25
26
27 ! ======================================================================
28 ! Declarations
29 ! ======================================================================
30
31 USE dimphy, ONLY: klon, klev
32 ! USE dimphy
33 USE phys_local_var_mod
34 USE phys_state_var_mod
35 USE mod_grid_phy_lmdz, ONLY: nbp_lev
36 IMPLICIT NONE
37 REAL,SAVE,ALLOCATABLE :: hthturb_gcssold(:)
38 REAL,SAVE,ALLOCATABLE :: hqturb_gcssold(:)
39 !$OMP THREADPRIVATE(hthturb_gcssold,hqturb_gcssold)
40 REAL,SAVE :: dtime_frcg
41 LOGICAL,SAVE :: turb_fcg_gcssold
42 LOGICAL,SAVE :: firstcall=.true.
43 !$OMP THREADPRIVATE(firstcall,turb_fcg_gcssold,dtime_frcg)
44 INTEGER abortphy
45 ! COMMON /turb_forcing/dtime_frcg, hthturb_gcssold, hqturb_gcssold, &
46 ! turb_fcg_gcssold
47
48 ! Arguments :
49 ! ------------
50 REAL zdu(klon, klev), zdv(klon, klev)
51 REAL zdt(klon, klev), zdq(klon, klev), zdql(klon, klev), zdqi(klon, klev)
52 CHARACTER *(*) text
53 REAL paprs(klon,klev+1)
54 INTEGER flag_inhib_tend ! if flag_inhib_tend != 0, tendencies are not added
55 INTEGER itap ! time step number
56
57 ! Local :
58 ! --------
59 ✗ REAL zzdt(klon, klev), zzdq(klon, klev)
60 INTEGER i, k
61
62 ✗ IF (firstcall) THEN
63 ✗ ALLOCATE(hthturb_gcssold(nbp_lev))
64 ✗ ALLOCATE(hqturb_gcssold(nbp_lev))
65 ✗ firstcall=.false.
66 ENDIF
67
68 ✗ IF (turb_fcg_gcssold) THEN
69 ✗ DO k = 1, klev
70 ✗ DO i = 1, klon
71 ✗ zzdt(i, k) = hthturb_gcssold(k)*dtime_frcg
72 ✗ zzdq(i, k) = hqturb_gcssold(k)*dtime_frcg
73 END DO
74 END DO
75 ✗ PRINT *, ' add_pbl_tend, dtime_frcg ', dtime_frcg
76 ✗ PRINT *, ' add_pbl_tend, zzdt ', zzdt
77 ✗ PRINT *, ' add_pbl_tend, zzdq ', zzdq
78 ✗ CALL add_phys_tend(zdu, zdv, zzdt, zzdq, zdql, zdqi, paprs, text,abortphy,flag_inhib_tend, itap, 0)
79 ELSE
80 ✗ CALL add_phys_tend(zdu, zdv, zdt, zdq, zdql, zdqi, paprs, text,abortphy,flag_inhib_tend, itap, 0)
81 END IF
82
83
84 ✗ RETURN
85 END SUBROUTINE add_pbl_tend
86 !
87 ! $Id: add_phys_tend.F90 2611 2016-08-03 15:41:26Z jyg $
88 !
89 6240 SUBROUTINE add_phys_tend (zdu,zdv,zdt,zdq,zdql,zdqi,paprs,text, &
90 abortphy,flag_inhib_tend, itap, diag_mode)
91 !======================================================================
92 ! Ajoute les tendances des variables physiques aux variables
93 ! d'etat de la dynamique t_seri, q_seri ...
94 ! On en profite pour faire des tests sur les tendances en question.
95 !======================================================================
96
97
98 !======================================================================
99 ! Declarations
100 !======================================================================
101
102 USE dimphy, ONLY: klon, klev
103 USE phys_state_var_mod, ONLY : phys_tstep
104 USE phys_local_var_mod, ONLY: u_seri, v_seri, ql_seri, qs_seri, q_seri, t_seri
105 USE phys_state_var_mod, ONLY: ftsol
106 USE geometry_mod, ONLY: longitude_deg, latitude_deg
107 USE print_control_mod, ONLY: prt_level
108 USE cmp_seri_mod
109 USE phys_output_var_mod, ONLY : d_qw_col, d_ql_col, d_qs_col, d_qt_col, d_ek_col, d_h_dair_col &
110 & , d_h_qw_col, d_h_ql_col, d_h_qs_col, d_h_col
111 IMPLICIT none
112 include "YOMCST.h"
113 include "clesphys.h"
114
115 ! Arguments :
116 !------------
117 REAL, DIMENSION(klon,klev), INTENT(IN) :: zdu, zdv
118 REAL, DIMENSION(klon,klev), INTENT(IN) :: zdt, zdql, zdqi
119 REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs
120 CHARACTER*(*), INTENT(IN) :: text
121 INTEGER, INTENT(IN) :: abortphy
122 INTEGER, INTENT(IN) :: flag_inhib_tend ! if not 0, tendencies are not added
123 INTEGER, INTENT(IN) :: itap ! time step number
124 INTEGER, INTENT(IN) :: diag_mode ! 0 -> normal effective mode
125 ! 1 -> only conservation stats are computed
126 !
127 REAL, DIMENSION(klon,klev), INTENT(INOUT) :: zdq
128
129 ! Local :
130 !--------
131 REAL zt,zq
132 REAL zq_int, zqp_int, zq_new
133
134 12480 REAL zqp(klev)
135
136 ! Save variables, used in diagnostic mode (diag_mode=1).
137 12480 REAL, DIMENSION(klon,klev) :: sav_u_seri, sav_v_seri
138 12480 REAL, DIMENSION(klon,klev) :: sav_ql_seri, sav_qs_seri, sav_q_seri
139 12480 REAL, DIMENSION(klon,klev) :: sav_t_seri
140 12480 REAL, DIMENSION(klon,klev) :: sav_zdq
141 !
142 INTEGER i, k,j, n
143 12480 INTEGER jadrs(klon*klev), jbad
144 12480 INTEGER jqadrs(klon*klev), jqbad
145 12480 INTEGER kadrs(klon*klev)
146 12480 INTEGER kqadrs(klon*klev)
147
148 12480 LOGICAL done(klon)
149
150 integer debug_level
151 logical, save :: first=.true.
152 !$OMP THREADPRIVATE(first)
153 !
154 !======================================================================
155 ! Variables for energy conservation tests
156 !======================================================================
157 !
158
159 ! zh_col------- total enthalpy of vertical air column
160 ! (air with watter vapour, liquid and solid) (J/m2)
161 ! zh_dair_col--- total enthalpy of dry air (J/m2)
162 ! zh_qw_col---- total enthalpy of watter vapour (J/m2)
163 ! zh_ql_col---- total enthalpy of liquid watter (J/m2)
164 ! zh_qs_col---- total enthalpy of solid watter (J/m2)
165 ! zqw_col------ total mass of watter vapour (kg/m2)
166 ! zql_col------ total mass of liquid watter (kg/m2)
167 ! zqs_col------ total mass of solid watter (kg/m2)
168 ! zek_col------ total kinetic energy (kg/m2)
169 !
170 12480 REAL zairm(klon, klev) ! layer air mass (kg/m2)
171 12480 REAL zqw_col(klon,2)
172 12480 REAL zql_col(klon,2)
173 12480 REAL zqs_col(klon,2)
174 12480 REAL zek_col(klon,2)
175 12480 REAL zh_dair_col(klon,2)
176 12480 REAL zh_qw_col(klon,2), zh_ql_col(klon,2), zh_qs_col(klon,2)
177 12480 REAL zh_col(klon,2)
178
179 REAL zcpvap, zcwat, zcice
180
181 !======================================================================
182 ! Initialisations
183
184
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 6240 times.
 6240 IF (prt_level >= 5) then
185 ✗ write (*,*) "In add_phys_tend, after ",text
186 ✗ call flush
187 end if
188
189 ! if flag_inhib_tend != 0, tendencies are not added
190
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 6240 times.
 6240 IF (flag_inhib_tend /= 0) then
191 ! If requiered, diagnostics are shown
192 ✗ IF (flag_inhib_tend > 0) then
193 ! print some diagnostics if xxx_seri have changed
194 ✗ call cmp_seri(flag_inhib_tend,text)
195 END IF
196 ✗ RETURN ! on n ajoute pas les tendance
197 END IF
198
199
1/2
✓ Branch 0 taken 6240 times.
✗ Branch 1 not taken.
 6240 IF (abortphy==1) RETURN ! on n ajoute pas les tendance si le modele
200 ! a deja plante.
201
202 6240 debug_level=10
203
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 6239 times.
 6240 if (first) then
204 1 print *,"TestJLD rcpv, rcw, rcs",rcpv, rcw, rcs
205 1 first=.false.
206 endif
207 !
208 ! print *,'add_phys_tend: paprs ',paprs
209 ! When in diagnostic mode, save initial values of out variables
210
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 6240 times.
 6240 IF (diag_mode == 1) THEN
211 ✗ sav_u_seri(:,:) = u_seri(:,:)
212 ✗ sav_v_seri(:,:) = v_seri(:,:)
213 ✗ sav_ql_seri(:,:) = ql_seri(:,:)
214 ✗ sav_qs_seri(:,:) = qs_seri(:,:)
215 ✗ sav_q_seri(:,:) = q_seri(:,:)
216 ✗ sav_t_seri(:,:) = t_seri(:,:)
217 ✗ sav_zdq(:,:) = zdq(:,:)
218 ENDIF ! (diag_mode == 1)
219 !======================================================================
220 ! Diagnostics for energy conservation tests
221 !======================================================================
222
2/2
✓ Branch 0 taken 243360 times.
✓ Branch 1 taken 6240 times.
 249600 DO k = 1, klev
223 ! layer air mass
224
2/2
✓ Branch 0 taken 241899840 times.
✓ Branch 1 taken 243360 times.
 242149440 zairm(:, k) = (paprs(:,k)-paprs(:,k+1))/rg
225 END DO
226
227
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 6240 times.
 6240 if (fl_ebil .GT. 0) then
228 ! ------------------------------------------------
229 ! Compute vertical sum for each atmospheric column
230 ! ------------------------------------------------
231 n=1 ! begining of time step
232
233 ✗ zcpvap = rcpv
234 zcwat = rcw
235 zcice = rcs
236
237 CALL integr_v(klon, klev, zcpvap, &
238 t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, zairm, &
239 zqw_col(:,n), zql_col(:,n), zqs_col(:,n), zek_col(:,n), zh_dair_col(:,n), &
240 ✗ zh_qw_col(:,n), zh_ql_col(:,n), zh_qs_col(:,n), zh_col(:,n))
241
242 end if ! end if (fl_ebil .GT. 0)
243
244 !======================================================================
245 ! Ajout des tendances sur le vent et l'eau liquide
246 !======================================================================
247
248
4/4
✓ Branch 0 taken 243360 times.
✓ Branch 1 taken 6240 times.
✓ Branch 2 taken 241899840 times.
✓ Branch 3 taken 243360 times.
 242149440 u_seri(:,:)=u_seri(:,:)+zdu(:,:)
249
4/4
✓ Branch 0 taken 243360 times.
✓ Branch 1 taken 6240 times.
✓ Branch 2 taken 241899840 times.
✓ Branch 3 taken 243360 times.
 242149440 v_seri(:,:)=v_seri(:,:)+zdv(:,:)
250
4/4
✓ Branch 0 taken 243360 times.
✓ Branch 1 taken 6240 times.
✓ Branch 2 taken 241899840 times.
✓ Branch 3 taken 243360 times.
 242149440 ql_seri(:,:)=ql_seri(:,:)+zdql(:,:)
251
4/4
✓ Branch 0 taken 243360 times.
✓ Branch 1 taken 6240 times.
✓ Branch 2 taken 241899840 times.
✓ Branch 3 taken 243360 times.
 242149440 qs_seri(:,:)=qs_seri(:,:)+zdqi(:,:)
252
253 !======================================================================
254 ! On ajoute les tendances de la temperature et de la vapeur d'eau
255 ! en verifiant que ca ne part pas dans les choux
256 !======================================================================
257
258 6240 jbad=0
259 6240 jqbad=0
260
2/2
✓ Branch 0 taken 243360 times.
✓ Branch 1 taken 6240 times.
 249600 DO k = 1, klev
261
2/2
✓ Branch 0 taken 241899840 times.
✓ Branch 1 taken 243360 times.
 242149440 DO i = 1, klon
262 241899840 zt=t_seri(i,k)+zdt(i,k)
263 241899840 zq=q_seri(i,k)+zdq(i,k)
264
3/6
✓ Branch 0 taken 241899840 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 241899840 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 241899840 times.
 241899840 IF ( zt>370. .or. zt<130. .or. abs(zdt(i,k))>50. ) then
265 ✗ jbad = jbad + 1
266 ✗ jadrs(jbad) = i
267 ✗ kadrs(jbad) = k
268 ENDIF
269
3/6
✓ Branch 0 taken 241899840 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 241899840 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 241899840 times.
 241899840 IF ( zq<0. .or. zq>0.1 .or. abs(zdq(i,k))>1.e-2 ) then
270 ✗ jqbad = jqbad + 1
271 ✗ jqadrs(jqbad) = i
272 ✗ kqadrs(jqbad) = k
273 ENDIF
274 241899840 t_seri(i,k)=zt
275 242143200 q_seri(i,k)=zq
276 ENDDO
277 ENDDO
278
279 !=====================================================================================
280 ! Impression et stop en cas de probleme important
281 !=====================================================================================
282
283
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 6240 times.
 6240 IF (jbad .GT. 0) THEN
284 ✗ DO j = 1, jbad
285 ✗ i=jadrs(j)
286 ✗ if(prt_level.ge.debug_level) THEN
287 ✗ print*,'PLANTAGE POUR LE POINT i lon lat =',&
288 ✗ i,longitude_deg(i),latitude_deg(i),text
289 ✗ print*,'l T dT Q dQ '
290 ✗ DO k = 1, klev
291 ✗ write(*,'(i3,2f14.4,2e14.2)') k,t_seri(i,k),zdt(i,k),q_seri(i,k),zdq(i,k)
292 ENDDO
293 ✗ call print_debug_phys(i,debug_level,text)
294 endif
295 ENDDO
296 ENDIF
297 !
298 !=====================================================================================
299 ! Impression, warning et correction en cas de probleme moins important
300 !=====================================================================================
301
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 6240 times.
 6240 IF (jqbad .GT. 0) THEN
302 ✗ done(:) = .false. !jyg
303 ✗ DO j = 1, jqbad
304 ✗ i=jqadrs(j)
305 ✗ if(prt_level.ge.debug_level) THEN
306 ✗ print*,'WARNING : EAU POUR LE POINT i lon lat =',&
307 ✗ i,longitude_deg(i),latitude_deg(i),text
308 ✗ print*,'l T dT Q dQ '
309 ✗ DO k = 1, klev
310 ✗ write(*,'(i3,2f14.4,2e14.2)') k,t_seri(i,k),zdt(i,k),q_seri(i,k),zdq(i,k)
311 ENDDO
312 endif
313 ✗ IF (ok_conserv_q) THEN
314 !jyg<20140228 Corrections pour conservation de l'eau
315 ✗ IF (.NOT.done(i)) THEN !jyg
316 ✗ DO k = 1, klev
317 ✗ zqp(k) = max(q_seri(i,k),1.e-15)
318 ENDDO
319 ✗ zq_int = 0.
320 ✗ zqp_int = 0.
321 ✗ DO k = 1, klev
322 ✗ zq_int = zq_int + q_seri(i,k)*(paprs(i,k)-paprs(i,k+1))/Rg
323 ✗ zqp_int = zqp_int + zqp(k) *(paprs(i,k)-paprs(i,k+1))/Rg
324 ENDDO
325 ✗ if(prt_level.ge.debug_level) THEN
326 ✗ print*,' cas q_seri<1.e-15 i k zq_int zqp_int zq_int/zqp_int :', &
327 ✗ i, kqadrs(j), zq_int, zqp_int, zq_int/zqp_int
328 endif
329 ✗ DO k = 1, klev
330 ✗ zq_new = zqp(k)*zq_int/zqp_int
331 ✗ zdq(i,k) = zdq(i,k) + zq_new - q_seri(i,k)
332 ✗ q_seri(i,k) = zq_new
333 ENDDO
334 ✗ done(i) = .true.
335 ENDIF !(.NOT.done(i))
336 ELSE
337 !jyg>
338 ✗ DO k = 1, klev
339 ✗ zq=q_seri(i,k)+zdq(i,k)
340 ✗ if (zq.lt.1.e-15) then
341 ✗ if (q_seri(i,k).lt.1.e-15) then
342 ✗ if(prt_level.ge.debug_level) THEN
343 ✗ print*,' cas q_seri<1.e-15 i k q_seri zq zdq :',i,k,q_seri(i,k),zq,zdq(i,k)
344 endif
345 ✗ q_seri(i,k)=1.e-15
346 ✗ zdq(i,k)=(1.e-15-q_seri(i,k))
347 endif
348 endif
349 ! zq=q_seri(i,k)+zdq(i,k)
350 ! if (zq.lt.1.e-15) then
351 ! zdq(i,k)=(1.e-15-q_seri(i,k))
352 ! endif
353 ENDDO
354 !jyg<20140228
355 ENDIF ! (ok_conserv_q)
356 !jyg>
357 ENDDO ! j = 1, jqbad
358 ENDIF
359 !
360
361 !IM ajout memes tests pour reverifier les jbad, jqbad beg
362 6240 jbad=0
363 6240 jqbad=0
364
2/2
✓ Branch 0 taken 243360 times.
✓ Branch 1 taken 6240 times.
 249600 DO k = 1, klev
365
2/2
✓ Branch 0 taken 241899840 times.
✓ Branch 1 taken 243360 times.
 242149440 DO i = 1, klon
366
3/6
✓ Branch 0 taken 241899840 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 241899840 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 241899840 times.
 241899840 IF ( t_seri(i,k)>370. .or. t_seri(i,k)<130. .or. abs(zdt(i,k))>50. ) then
367 ✗ jbad = jbad + 1
368 ✗ jadrs(jbad) = i
369 ! if(prt_level.ge.debug_level) THEN
370 ! print*,'cas2 i k t_seri zdt',i,k,t_seri(i,k),zdt(i,k)
371 ! endif
372 ENDIF
373
3/6
✓ Branch 0 taken 241899840 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 241899840 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 241899840 times.
 242143200 IF ( q_seri(i,k)<0. .or. q_seri(i,k)>0.1 .or. abs(zdq(i,k))>1.e-2 ) then
374 ✗ jqbad = jqbad + 1
375 ✗ jqadrs(jqbad) = i
376 ✗ kqadrs(jqbad) = k
377 ! if(prt_level.ge.debug_level) THEN
378 ! print*,'cas2 i k q_seri zdq',i,k,q_seri(i,k),zdq(i,k)
379 ! endif
380 ENDIF
381 ENDDO
382 ENDDO
383
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 6240 times.
 6240 IF (jbad .GT. 0) THEN
384 ✗ DO j = 1, jbad
385 ✗ i=jadrs(j)
386 ✗ k=kadrs(j)
387 ✗ if(prt_level.ge.debug_level) THEN
388 ✗ print*,'PLANTAGE2 POUR LE POINT i itap lon lat txt jbad zdt t',&
389 ✗ i,itap,longitude_deg(i),latitude_deg(i),text,jbad, &
390 ✗ & zdt(i,k),t_seri(i,k)-zdt(i,k)
391 !!! if(prt_level.ge.10.and.itap.GE.229.and.i.EQ.3027) THEN
392 ✗ print*,'l T dT Q dQ '
393 ✗ DO k = 1, klev
394 ✗ write(*,'(i3,2f14.4,2e14.2)') k,t_seri(i,k),zdt(i,k),q_seri(i,k),zdq(i,k)
395 ENDDO
396 ✗ call print_debug_phys(i,debug_level,text)
397 endif
398 ENDDO
399 ENDIF
400 !
401
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 6240 times.
 6240 IF (jqbad .GT. 0) THEN
402 ✗ DO j = 1, jqbad
403 ✗ i=jqadrs(j)
404 ✗ k=kqadrs(j)
405 ✗ if(prt_level.ge.debug_level) THEN
406 ✗ print*,'WARNING : EAU2 POUR LE POINT i itap lon lat txt jqbad zdq q zdql ql',&
407 ✗ i,itap,longitude_deg(i),latitude_deg(i),text,jqbad,&
408 ✗ & zdq(i,k), q_seri(i,k)-zdq(i,k), zdql(i,k), ql_seri(i,k)-zdql(i,k)
409 !!! if(prt_level.ge.10.and.itap.GE.229.and.i.EQ.3027) THEN
410 ✗ print*,'l T dT Q dQ '
411 ✗ DO k = 1, klev
412 ✗ write(*,'(i3,2f14.4,2e14.2)') k,t_seri(i,k),zdt(i,k),q_seri(i,k),zdq(i,k)
413 ENDDO
414 ✗ call print_debug_phys(i,debug_level,text)
415 endif
416 ENDDO
417 ENDIF
418
419 !======================================================================
420 ! Contrôle des min/max pour arrêt du modèle
421 ! Si le modele est en mode abortphy, on retire les tendances qu'on
422 ! vient d'ajouter. Pas exactement parce qu'on ne tient pas compte des
423 ! seuils.
424 !======================================================================
425
426 6240 CALL hgardfou(t_seri,ftsol,text,abortphy)
427 IF (abortphy==1) THEN
428 Print*,'ERROR ABORT hgardfou dans ',text
429 ! JLD pourquoi on ne modifie pas de meme t_seri et q_seri ?
430 u_seri(:,:)=u_seri(:,:)-zdu(:,:)
431 v_seri(:,:)=v_seri(:,:)-zdv(:,:)
432 ql_seri(:,:)=ql_seri(:,:)-zdql(:,:)
433 qs_seri(:,:)=qs_seri(:,:)-zdqi(:,:)
434 ENDIF
435
436 !======================================================================
437 ! Diagnostics for energy conservation tests
438 !======================================================================
439
440
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 6240 times.
 6240 if (fl_ebil .GT. 0) then
441
442 ! ------------------------------------------------
443 ! Compute vertical sum for each atmospheric column
444 ! ------------------------------------------------
445 n=2 ! end of time step
446
447 CALL integr_v(klon, klev, zcpvap, &
448 t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, zairm, &
449 zqw_col(:,n), zql_col(:,n), zqs_col(:,n), zek_col(:,n), zh_dair_col(:,n), &
450 ✗ zh_qw_col(:,n), zh_ql_col(:,n), zh_qs_col(:,n), zh_col(:,n))
451
452 ! ------------------------------------------------
453 ! Compute the changes by unit of time
454 ! ------------------------------------------------
455
456 ✗ d_qw_col(:) = (zqw_col(:,2)-zqw_col(:,1))/phys_tstep
457 ✗ d_ql_col(:) = (zql_col(:,2)-zql_col(:,1))/phys_tstep
458 ✗ d_qs_col(:) = (zqs_col(:,2)-zqs_col(:,1))/phys_tstep
459 ✗ d_qt_col(:) = d_qw_col(:) + d_ql_col(:) + d_qs_col(:)
460
461 ✗ d_ek_col(:) = (zek_col(:,2)-zek_col(:,1))/phys_tstep
462
463 ✗ d_h_dair_col(:) = (zh_dair_col(:,2)-zh_dair_col(:,1))/phys_tstep
464 ✗ d_h_qw_col(:) = (zh_qw_col(:,2)-zh_qw_col(:,1))/phys_tstep
465 ✗ d_h_ql_col(:) = (zh_ql_col(:,2)-zh_ql_col(:,1))/phys_tstep
466 ✗ d_h_qs_col(:) = (zh_qs_col(:,2)-zh_qs_col(:,1))/phys_tstep
467
468 ✗ d_h_col = (zh_col(:,2)-zh_col(:,1))/phys_tstep
469
470 end if ! end if (fl_ebil .GT. 0)
471 !
472 ! When in diagnostic mode, restore "out" variables to initial values.
473
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 6240 times.
 6240 IF (diag_mode == 1) THEN
474 ✗ u_seri(:,:) = sav_u_seri(:,:)
475 ✗ v_seri(:,:) = sav_v_seri(:,:)
476 ✗ ql_seri(:,:) = sav_ql_seri(:,:)
477 ✗ qs_seri(:,:) = sav_qs_seri(:,:)
478 ✗ q_seri(:,:) = sav_q_seri(:,:)
479 ✗ t_seri(:,:) = sav_t_seri(:,:)
480 ✗ zdq(:,:) = sav_zdq(:,:)
481 ENDIF ! (mode == 1)
482
483 RETURN
484 END SUBROUTINE add_phys_tend
485
486 ✗ SUBROUTINE diag_phys_tend (nlon, nlev, uu, vv, temp, qv, ql, qs, &
487 ✗ zdu,zdv,zdt,zdq,zdql,zdqs,paprs,text)
488 !======================================================================
489 ! Ajoute les tendances des variables physiques aux variables
490 ! d'etat de la dynamique t_seri, q_seri ...
491 ! On en profite pour faire des tests sur les tendances en question.
492 !======================================================================
493
494
495 !======================================================================
496 ! Declarations
497 !======================================================================
498
499 ✗ USE phys_state_var_mod, ONLY : phys_tstep, ftsol
500 USE geometry_mod, ONLY: longitude_deg, latitude_deg
501 USE print_control_mod, ONLY: prt_level
502 USE cmp_seri_mod
503 USE phys_output_var_mod, ONLY : d_qw_col, d_ql_col, d_qs_col, d_qt_col, d_ek_col, d_h_dair_col &
504 & , d_h_qw_col, d_h_ql_col, d_h_qs_col, d_h_col
505 IMPLICIT none
506 include "YOMCST.h"
507 include "clesphys.h"
508
509 ! Arguments :
510 !------------
511 INTEGER, INTENT(IN) :: nlon, nlev
512 REAL, DIMENSION(nlon,nlev), INTENT(IN) :: uu, vv
513 REAL, DIMENSION(nlon,nlev), INTENT(IN) :: temp, qv, ql, qs
514 REAL, DIMENSION(nlon,nlev), INTENT(IN) :: zdu, zdv
515 REAL, DIMENSION(nlon,nlev), INTENT(IN) :: zdt, zdq, zdql, zdqs
516 REAL, DIMENSION(nlon,nlev+1), INTENT(IN) :: paprs
517 CHARACTER*(*), INTENT(IN) :: text
518
519 ! Local :
520 !--------
521 ✗ REAL, DIMENSION(nlon,nlev) :: uu_n, vv_n
522 ✗ REAL, DIMENSION(nlon,nlev) :: temp_n, qv_n, ql_n, qs_n
523
524
525 !
526 INTEGER k, n
527
528 integer debug_level
529 logical, save :: first=.true.
530 !$OMP THREADPRIVATE(first)
531 !
532 !======================================================================
533 ! Variables for energy conservation tests
534 !======================================================================
535 !
536
537 ! zh_col------- total enthalpy of vertical air column
538 ! (air with watter vapour, liquid and solid) (J/m2)
539 ! zh_dair_col--- total enthalpy of dry air (J/m2)
540 ! zh_qw_col---- total enthalpy of watter vapour (J/m2)
541 ! zh_ql_col---- total enthalpy of liquid watter (J/m2)
542 ! zh_qs_col---- total enthalpy of solid watter (J/m2)
543 ! zqw_col------ total mass of watter vapour (kg/m2)
544 ! zql_col------ total mass of liquid watter (kg/m2)
545 ! zqs_col------ total mass of solid watter (kg/m2)
546 ! zek_col------ total kinetic energy (kg/m2)
547 !
548 ✗ REAL zairm(nlon, nlev) ! layer air mass (kg/m2)
549 ✗ REAL zqw_col(nlon,2)
550 ✗ REAL zql_col(nlon,2)
551 ✗ REAL zqs_col(nlon,2)
552 ✗ REAL zek_col(nlon,2)
553 ✗ REAL zh_dair_col(nlon,2)
554 ✗ REAL zh_qw_col(nlon,2), zh_ql_col(nlon,2), zh_qs_col(nlon,2)
555 ✗ REAL zh_col(nlon,2)
556
557 !======================================================================
558 ! Initialisations
559
560 ✗ IF (prt_level >= 5) then
561 ✗ write (*,*) "In diag_phys_tend, after ",text
562 ✗ call flush
563 end if
564
565 debug_level=10
566 ✗ if (first) then
567 ✗ print *,"TestJLD rcpv, rcw, rcs",rcpv, rcw, rcs
568 ✗ first=.false.
569 endif
570 !
571 ! print *,'add_phys_tend: paprs ',paprs
572 !======================================================================
573 ! Diagnostics for energy conservation tests
574 !======================================================================
575 ✗ DO k = 1, nlev
576 ! layer air mass
577 ✗ zairm(:, k) = (paprs(:,k)-paprs(:,k+1))/rg
578 END DO
579
580 ✗ if (fl_ebil .GT. 0) then
581 ! ------------------------------------------------
582 ! Compute vertical sum for each atmospheric column
583 ! ------------------------------------------------
584 n=1 ! begining of time step
585
586 CALL integr_v(nlon, nlev, rcpv, &
587 temp, qv, ql, qs, uu, vv, zairm, &
588 zqw_col(:,n), zql_col(:,n), zqs_col(:,n), zek_col(:,n), zh_dair_col(:,n), &
589 ✗ zh_qw_col(:,n), zh_ql_col(:,n), zh_qs_col(:,n), zh_col(:,n))
590
591 end if ! end if (fl_ebil .GT. 0)
592
593 !======================================================================
594 ! Ajout des tendances sur le vent, la temperature et les diverses phases de l'eau
595 !======================================================================
596
597 ✗ uu_n(:,:)=uu(:,:)+zdu(:,:)
598 ✗ vv_n(:,:)=vv(:,:)+zdv(:,:)
599 ✗ qv_n(:,:)=qv(:,:)+zdq(:,:)
600 ✗ ql_n(:,:)=ql(:,:)+zdql(:,:)
601 ✗ qs_n(:,:)=qs(:,:)+zdqs(:,:)
602 ✗ temp_n(:,:)=temp(:,:)+zdt(:,:)
603
604
605
606 !======================================================================
607 ! Diagnostics for energy conservation tests
608 !======================================================================
609
610 ✗ if (fl_ebil .GT. 0) then
611
612 ! ------------------------------------------------
613 ! Compute vertical sum for each atmospheric column
614 ! ------------------------------------------------
615 n=2 ! end of time step
616
617 CALL integr_v(nlon, nlev, rcpv, &
618 temp_n, qv_n, ql_n, qs_n, uu_n, vv_n, zairm, &
619 zqw_col(:,n), zql_col(:,n), zqs_col(:,n), zek_col(:,n), zh_dair_col(:,n), &
620 ✗ zh_qw_col(:,n), zh_ql_col(:,n), zh_qs_col(:,n), zh_col(:,n))
621
622 ! ------------------------------------------------
623 ! Compute the changes by unit of time
624 ! ------------------------------------------------
625
626 ✗ d_qw_col(:) = (zqw_col(:,2)-zqw_col(:,1))/phys_tstep
627 ✗ d_ql_col(:) = (zql_col(:,2)-zql_col(:,1))/phys_tstep
628 ✗ d_qs_col(:) = (zqs_col(:,2)-zqs_col(:,1))/phys_tstep
629 ✗ d_qt_col(:) = d_qw_col(:) + d_ql_col(:) + d_qs_col(:)
630
631 ✗ d_ek_col(:) = (zek_col(:,2)-zek_col(:,1))/phys_tstep
632
633 ✗ print *,'zdu ', zdu
634 ✗ print *,'zdv ', zdv
635 ✗ print *,'d_ek_col, zek_col(2), zek_col(1) ',d_ek_col(1), zek_col(1,2), zek_col(1,1)
636
637 ✗ d_h_dair_col(:) = (zh_dair_col(:,2)-zh_dair_col(:,1))/phys_tstep
638 ✗ d_h_qw_col(:) = (zh_qw_col(:,2)-zh_qw_col(:,1))/phys_tstep
639 ✗ d_h_ql_col(:) = (zh_ql_col(:,2)-zh_ql_col(:,1))/phys_tstep
640 ✗ d_h_qs_col(:) = (zh_qs_col(:,2)-zh_qs_col(:,1))/phys_tstep
641
642 ✗ d_h_col = (zh_col(:,2)-zh_col(:,1))/phys_tstep
643
644 end if ! end if (fl_ebil .GT. 0)
645 !
646
647 ✗ RETURN
648 END SUBROUTINE diag_phys_tend
649
650 ✗ SUBROUTINE integr_v(nlon, nlev, zcpvap, &
651 ✗ temp, qv, ql, qs, uu, vv, zairm, &
652 ✗ zqw_col, zql_col, zqs_col, zek_col, zh_dair_col, &
653 zh_qw_col, zh_ql_col, zh_qs_col, zh_col)
654
655 IMPLICIT none
656 include "YOMCST.h"
657
658 INTEGER, INTENT(IN) :: nlon,nlev
659 REAL, INTENT(IN) :: zcpvap
660 REAL, DIMENSION(nlon,nlev), INTENT(IN) :: temp, qv, ql, qs, uu, vv
661 REAL, DIMENSION(nlon,nlev), INTENT(IN) :: zairm
662 REAL, DIMENSION(nlon), INTENT(OUT) :: zqw_col
663 REAL, DIMENSION(nlon), INTENT(OUT) :: zql_col
664 REAL, DIMENSION(nlon), INTENT(OUT) :: zqs_col
665 REAL, DIMENSION(nlon), INTENT(OUT) :: zek_col
666 REAL, DIMENSION(nlon), INTENT(OUT) :: zh_dair_col
667 REAL, DIMENSION(nlon), INTENT(OUT) :: zh_qw_col
668 REAL, DIMENSION(nlon), INTENT(OUT) :: zh_ql_col
669 REAL, DIMENSION(nlon), INTENT(OUT) :: zh_qs_col
670 REAL, DIMENSION(nlon), INTENT(OUT) :: zh_col
671
672 INTEGER :: i, k
673
674
675 ! Reset variables
676 ✗ zqw_col(:) = 0.
677 ✗ zql_col(:) = 0.
678 ✗ zqs_col(:) = 0.
679 ✗ zek_col(:) = 0.
680 ✗ zh_dair_col(:) = 0.
681 ✗ zh_qw_col(:) = 0.
682 ✗ zh_ql_col(:) = 0.
683 ✗ zh_qs_col(:) = 0.
684
685 !JLD write (*,*) "rcpd, zcpvap, zcwat, zcice ",rcpd, zcpvap, zcwat, zcice
686
687 ! ------------------------------------------------
688 ! Compute vertical sum for each atmospheric column
689 ! ------------------------------------------------
690 ✗ DO k = 1, nlev
691 ✗ DO i = 1, nlon
692 ! Watter mass
693 ✗ zqw_col(i) = zqw_col(i) + qv(i, k)*zairm(i, k)
694 ✗ zql_col(i) = zql_col(i) + ql(i, k)*zairm(i, k)
695 ✗ zqs_col(i) = zqs_col(i) + qs(i, k)*zairm(i, k)
696 ! Kinetic Energy
697 ✗ zek_col(i) = zek_col(i) + 0.5*(uu(i,k)**2+vv(i,k)**2)*zairm(i, k)
698 ! Air enthalpy : dry air, water vapour, liquid, solid
699 zh_dair_col(i) = zh_dair_col(i) + rcpd*(1.-qv(i,k)-ql(i,k)-qs(i,k))* &
700 ✗ zairm(i, k)*temp(i, k)
701 ✗ zh_qw_col(i) = zh_qw_col(i) + zcpvap*temp(i, k) *qv(i, k)*zairm(i, k) !jyg
702 ✗ zh_ql_col(i) = zh_ql_col(i) + (zcpvap*temp(i, k) - rlvtt)*ql(i, k)*zairm(i, k) !jyg
703 ✗ zh_qs_col(i) = zh_qs_col(i) + (zcpvap*temp(i, k) - rlstt)*qs(i, k)*zairm(i, k) !jyg
704 END DO
705 END DO
706 ! compute total air enthalpy
707 ✗ zh_col(:) = zh_dair_col(:) + zh_qw_col(:) + zh_ql_col(:) + zh_qs_col(:)
708
709 ✗ END SUBROUTINE integr_v
710
711 6720 SUBROUTINE prt_enerbil (text, itap)
712 !======================================================================
713 ! Print enenrgy budget diagnotics for the 1D case
714 !======================================================================
715
716 !======================================================================
717 ! Declarations
718 !======================================================================
719
720 ✗ USE dimphy, ONLY: klon, klev
721 USE phys_state_var_mod, ONLY : phys_tstep
722 USE phys_state_var_mod, ONLY : topsw, toplw, solsw, sollw, rain_con, snow_con
723 USE geometry_mod, ONLY: longitude_deg, latitude_deg
724 USE print_control_mod, ONLY: prt_level
725 USE cmp_seri_mod
726 USE phys_output_var_mod, ONLY : d_qw_col, d_ql_col, d_qs_col, d_qt_col, d_ek_col, d_h_dair_col &
727 & , d_h_qw_col, d_h_ql_col, d_h_qs_col, d_h_col
728 USE phys_local_var_mod, ONLY: evap, sens
729 USE phys_local_var_mod, ONLY: u_seri, v_seri, ql_seri, qs_seri, q_seri, t_seri &
730 & , rain_lsc, snow_lsc
731 USE climb_hq_mod, ONLY : d_h_col_vdf, f_h_bnd
732 IMPLICIT none
733 include "YOMCST.h"
734
735 ! Arguments :
736 !------------
737 CHARACTER*(*) text ! text specifing the involved parametrization
738 integer itap ! time step number
739 ! local variables
740 ! ---------------
741 real bilq_seuil, bilh_seuil ! thresold on error in Q and H budget
742 real bilq_error, bilh_error ! erros in Q and H budget
743 real bilq_bnd, bilh_bnd ! Q and H budget due to exchange with boundaries
744 integer bilq_ok, bilh_ok
745 CHARACTER*(12) status
746
747 bilq_seuil = 1.E-10
748 bilh_seuil = 1.E-1
749 bilq_ok=0
750 bilh_ok=0
751
752 !!print *,'prt_level:',prt_level,' fl_ebil:',fl_ebil,' fl_cor_ebil:',fl_cor_ebil
753
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 6720 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 6720 if ( (fl_ebil .GT. 0) .and. (klon .EQ. 1)) then
754
755 ✗ bilq_bnd = 0.
756 ✗ bilh_bnd = 0.
757
758 ✗ param: SELECT CASE (text)
759 CASE("vdf") param
760 ✗ bilq_bnd = evap(1)
761 ✗ bilh_bnd = sens(1)+(rcpv-rcpd)*evap(1)*t_seri(1,1)
762 CASE("lsc") param
763 ✗ bilq_bnd = - rain_lsc(1) - snow_lsc(1)
764 bilh_bnd = (-(rcw-rcpd)*t_seri(1,1) + rlvtt) * rain_lsc(1) &
765 ✗ & + (-(rcs-rcpd)*t_seri(1,1) + rlstt) * snow_lsc(1)
766 CASE("convection") param
767 ✗ bilq_bnd = - rain_con(1) - snow_con(1)
768 bilh_bnd = (-(rcw-rcpd)*t_seri(1,1) + rlvtt) * rain_con(1) &
769 ✗ & + (-(rcs-rcpd)*t_seri(1,1) + rlstt) * snow_con(1)
770 CASE("SW") param
771 ✗ bilh_bnd = topsw(1) - solsw(1)
772 CASE("LW") param
773 ✗ bilh_bnd = -(toplw(1) + sollw(1))
774 CASE DEFAULT param
775 bilq_bnd = 0.
776 ✗ bilh_bnd = 0.
777 END SELECT param
778
779 ✗ bilq_error = d_qt_col(1) - bilq_bnd
780 ✗ bilh_error = d_h_col(1) - bilh_bnd
781 ! are the errors too large?
782 ✗ if ( abs(bilq_error) .gt. bilq_seuil) bilq_ok=1
783 ✗ if ( abs(bilh_error) .gt. bilh_seuil) bilh_ok=1
784 !
785 ! Print diagnostics
786 ! =================
787 ✗ if ( (bilq_ok .eq. 0).and.(bilh_ok .eq. 0) ) then
788 ✗ status="enerbil-OK"
789 else
790 ✗ status="enerbil-PB"
791 end if
792
793 ✗ if ( prt_level .GE. 3) then
794 ✗ write(*,9010) text,status," itap:",itap,"enerbilERROR: Q", bilq_error," H", bilh_error
795 9010 format (1x,A8,2x,A12,A6,I4,A18,E15.6,A5,E15.6)
796 end if
797 ✗ if ( prt_level .GE. 3) then
798 ✗ write(*,9000) text,"enerbil: Q,H,KE budget", d_qt_col(1), d_h_col(1),d_ek_col(1)
799 end if
800 ✗ if ( prt_level .GE. 5) then
801 ✗ write(*,9000) text,"enerbil at boundaries: Q, H",bilq_bnd, bilh_bnd
802 ✗ write(*,9000) text,"enerbil: water budget",d_qt_col(1),d_qw_col(1),d_ql_col(1),d_qs_col(1)
803 ✗ write(*,9000) text,"enerbil: enthalpy budget",d_h_col(1),d_h_dair_col(1),d_h_qw_col(1),d_h_ql_col(1),d_h_qs_col(1)
804 end if
805
806 ✗ specific_diag: SELECT CASE (text)
807 CASE("vdf") specific_diag
808 ✗ if ( prt_level .GE. 5) then
809 ✗ write(*,9000) text,"enerbil: d_h, bilh, sens,t_seri", d_h_col(1), bilh_bnd, sens(1), t_seri(1,1)
810 ✗ write(*,9000) text,"enerbil: d_h_col_vdf, f_h, diff",d_h_col_vdf, f_h_bnd, bilh_bnd-sens(1)
811 end if
812 CASE("lsc") specific_diag
813 ✗ if ( prt_level .GE. 5) then
814 ✗ write(*,9000) text,"enerbil: rain, bil_lat, bil_sens", rain_lsc(1), rlvtt * rain_lsc(1), -(rcw-rcpd)*t_seri(1,1) * rain_lsc(1)
815 ✗ write(*,9000) text,"enerbil: snow, bil_lat, bil_sens", snow_lsc(1), rlstt * snow_lsc(1), -(rcs-rcpd)*t_seri(1,1) * snow_lsc(1)
816 end if
817 CASE("convection") specific_diag
818 ✗ if ( prt_level .GE. 5) then
819 ✗ write(*,9000) text,"enerbil: rain, bil_lat, bil_sens", rain_con(1), rlvtt * rain_con(1), -(rcw-rcpd)*t_seri(1,1) * rain_con(1)
820 ✗ write(*,9000) text,"enerbil: snow, bil_lat, bil_sens", snow_con(1), rlstt * snow_con(1), -(rcs-rcpd)*t_seri(1,1) * snow_con(1)
821 end if
822 END SELECT specific_diag
823
824 9000 format (1x,A8,2x,A35,10E15.6)
825
826 end if ! end if (fl_ebil .GT. 0)
827
828 6720 END SUBROUTINE prt_enerbil
829
830 END MODULE add_phys_tend_mod
831