Line data Source code
1 : module qneg_module
2 :
3 : use shr_kind_mod, only: r8 => shr_kind_r8, CS => SHR_KIND_CS
4 : use perf_mod, only: t_startf, t_stopf
5 : use cam_logfile, only: iulog
6 : use cam_abortutils, only: endrun
7 : use shr_sys_mod, only: shr_sys_flush
8 : use cam_history_support, only: max_fieldname_len
9 : use ppgrid, only: pcols
10 : use constituents, only: pcnst, cnst_name
11 :
12 : implicit none
13 : private
14 : save
15 :
16 : ! Public interface.
17 :
18 : public :: qneg_readnl
19 : public :: qneg_init
20 : public :: qneg3
21 : public :: qneg4
22 : public :: qneg_print_summary
23 :
24 : ! Private module variables
25 : character(len=8) :: print_qneg_warn
26 : logical :: log_warnings = .false.
27 : logical :: collect_stats = .false.
28 : logical :: timestep_reset = .false.
29 :
30 : real(r8), parameter :: tol = 1.e-12_r8
31 : real(r8), parameter :: worst_reset = 1.e35_r8
32 :
33 : ! Diagnostic field names
34 : integer, parameter :: num_diag_fields = (2 * pcnst) + 1
35 : character(len=max_fieldname_len) :: diag_names(num_diag_fields)
36 : logical :: cnst_out_calc = .false.
37 : logical :: cnst_outfld(num_diag_fields) = .false.
38 :
39 : ! Summary buffers
40 : integer, parameter :: num3_bins = 24
41 : integer, parameter :: num4_bins = 4
42 : character(len=CS) :: qneg3_warn_labels(num3_bins) = ''
43 : character(len=CS) :: qneg4_warn_labels(num4_bins) = ''
44 : integer :: qneg3_warn_num(pcnst, num3_bins) = 0
45 : integer :: qneg4_warn_num(num4_bins) = 0
46 : real(r8) :: qneg3_warn_worst(pcnst, num3_bins) = worst_reset
47 : real(r8) :: qneg4_warn_worst(num4_bins) = worst_reset
48 :
49 : private :: calc_cnst_out
50 : private :: find_index3
51 : private :: find_index4
52 : interface reset_stats
53 : module procedure reset_stats_scalar
54 : module procedure reset_stats_array
55 : end interface reset_stats
56 :
57 : contains
58 :
59 1536 : subroutine qneg_readnl(nlfile)
60 : use namelist_utils, only: find_group_name
61 : use units, only: getunit, freeunit
62 : use spmd_utils, only: mpicom, mstrid=>masterprocid, mpi_character, masterproc
63 : ! File containing namelist input.
64 : character(len=*), intent(in) :: nlfile
65 :
66 : ! Local variables
67 : integer :: unitn, ierr
68 : character(len=*), parameter :: sub = 'qneg_readnl'
69 :
70 : namelist /qneg_nl/ print_qneg_warn
71 :
72 1536 : print_qneg_warn = ''
73 :
74 1536 : if (masterproc) then
75 2 : unitn = getunit()
76 2 : open( unitn, file=trim(nlfile), status='old' )
77 2 : call find_group_name(unitn, 'qneg_nl', status=ierr)
78 2 : if (ierr == 0) then
79 2 : read(unitn, qneg_nl, iostat=ierr)
80 2 : if (ierr /= 0) then
81 0 : call endrun(sub // ':: ERROR reading namelist qneg_nl')
82 : end if
83 : end if
84 2 : close(unitn)
85 2 : call freeunit(unitn)
86 : end if
87 :
88 1536 : call mpi_bcast(print_qneg_warn, len(print_qneg_warn), mpi_character, mstrid, mpicom, ierr)
89 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: print_qneg_warn")
90 :
91 3072 : select case(trim(print_qneg_warn))
92 : case('summary')
93 1536 : collect_stats = .true.
94 1536 : timestep_reset = .false.
95 : case('timestep')
96 0 : collect_stats = .true.
97 0 : timestep_reset = .true.
98 : case('off')
99 0 : collect_stats = .false.
100 0 : timestep_reset = .false.
101 : case default
102 3072 : call endrun(sub//" FATAL: '"//trim(print_qneg_warn)//"' is not a valid value for print_qneg_warn")
103 : end select
104 :
105 1536 : if (masterproc) then
106 2 : if (collect_stats) then
107 2 : if (timestep_reset) then
108 0 : write(iulog, *) sub, ": QNEG statistics will be collected and printed for each timestep"
109 : else
110 2 : write(iulog, *) sub, ": QNEG statistics will be collected and printed at the end of the run"
111 : end if
112 : else
113 0 : write(iulog, *) sub, ": QNEG statistics will not be collected"
114 : end if
115 : end if
116 :
117 1536 : end subroutine qneg_readnl
118 :
119 1536 : subroutine qneg_init()
120 : use cam_history, only: addfld, horiz_only
121 : use constituents, only: cnst_longname
122 :
123 : integer :: index
124 :
125 64512 : do index = 1, pcnst
126 62976 : diag_names(index) = trim(cnst_name(index))//'_qneg3'
127 : call addfld(diag_names(index), (/ 'lev' /), 'I', 'kg/kg', &
128 125952 : trim(cnst_longname(index))//' QNEG3 error (cell)')
129 62976 : diag_names(pcnst+index) = trim(cnst_name(index))//'_qneg3_col'
130 : call addfld(diag_names(pcnst+index), horiz_only, 'I', 'kg/kg', &
131 64512 : trim(cnst_longname(index))//' QNEG3 error (column)')
132 : end do
133 1536 : diag_names((2*pcnst) + 1) = 'qflux_exceeded'
134 : call addfld(diag_names((2*pcnst) + 1), horiz_only, 'I', 'kg/m^2/s', &
135 1536 : 'qflux excess (QNEG4)')
136 :
137 1536 : end subroutine qneg_init
138 :
139 1536 : subroutine calc_cnst_out()
140 1536 : use cam_history, only: hist_fld_active, history_initialized
141 : integer :: index
142 :
143 1536 : if (history_initialized()) then
144 : ! to protect against routines that call qneg3 too early
145 129024 : do index = 1, num_diag_fields
146 129024 : cnst_outfld(index) = hist_fld_active(trim(diag_names(index)))
147 : end do
148 1536 : cnst_out_calc = .true.
149 : end if
150 :
151 1536 : end subroutine calc_cnst_out
152 :
153 743173128 : integer function find_index3(nam) result(index)
154 : ! Find a valid or new index for 'nam' entries
155 : character(len=*), intent(in) :: nam
156 :
157 : integer :: i
158 :
159 743173128 : index = -1
160 9419391144 : do i = 1, num3_bins
161 9419391144 : if (trim(nam) == trim(qneg3_warn_labels(i))) then
162 : ! We found this entry, return its index
163 : index = i
164 743173128 : exit
165 8676248736 : else if (len_trim(qneg3_warn_labels(i)) == 0) then
166 : ! We have run out of known entries, use a new one and reset its stats
167 30720 : qneg3_warn_labels(i) = nam
168 30720 : index = i
169 30720 : call reset_stats(qneg3_warn_num(:, index), qneg3_warn_worst(:,index))
170 30720 : exit
171 : end if
172 : end do
173 743174664 : end function find_index3
174 :
175 1489176 : integer function find_index4(nam) result(index)
176 : ! Find a valid or new index for 'nam' entries
177 : character(len=*), intent(in) :: nam
178 :
179 : integer :: i
180 :
181 1489176 : index = -1
182 1489176 : do i = 1, num4_bins
183 1489176 : if (trim(nam) == trim(qneg4_warn_labels(i))) then
184 : ! We found this entry, return its index
185 : index = i
186 1489176 : exit
187 1536 : else if (len_trim(qneg4_warn_labels(i)) == 0) then
188 : ! We have run out of known entries, use a new one and reset its stats
189 1536 : qneg4_warn_labels(i) = nam
190 1536 : index = i
191 : call reset_stats(qneg4_warn_num(index), qneg4_warn_worst(index))
192 1536 : exit
193 : end if
194 : end do
195 1489176 : end function find_index4
196 :
197 743173128 : subroutine qneg3 (subnam, idx, ncol, ncold, lver, lconst_beg, &
198 743173128 : lconst_end, qmin, q)
199 : !-----------------------------------------------------------------------
200 : !
201 : ! Purpose:
202 : ! Check moisture and tracers for minimum value, reset any below
203 : ! minimum value to minimum value and return information to allow
204 : ! warning message to be printed. The global average is NOT preserved.
205 : !
206 : ! Method:
207 : ! <Describe the algorithm(s) used in the routine.>
208 : ! <Also include any applicable external references.>
209 : !
210 : ! Author: J. Rosinski
211 : !
212 : !-----------------------------------------------------------------------
213 : use cam_history, only: outfld
214 :
215 : !------------------------------Arguments--------------------------------
216 : !
217 : ! Input arguments
218 : !
219 : character(len=*), intent(in) :: subnam ! name of calling routine
220 :
221 : integer, intent(in) :: idx ! chunk/latitude index
222 : integer, intent(in) :: ncol ! number of atmospheric columns
223 : integer, intent(in) :: ncold ! declared number of atmospheric columns
224 : integer, intent(in) :: lver ! number of vertical levels in column
225 : integer, intent(in) :: lconst_beg ! beginning constituent
226 : integer, intent(in) :: lconst_end ! ending constituent
227 :
228 : real(r8), intent(in) :: qmin(lconst_beg:lconst_end) ! Global minimum constituent concentration
229 :
230 : !
231 : ! Input/Output arguments
232 : !
233 : real(r8), intent(inout) :: q(ncold,lver,lconst_beg:lconst_end) ! moisture/tracer field
234 : !
235 : !---------------------------Local workspace-----------------------------
236 : !
237 : integer :: nvals ! number of values found < qmin
238 : integer :: i, k ! longitude, level indices
239 : integer :: index ! For storing stats
240 : integer :: m ! constituent index
241 : integer :: iw,kw ! i,k indices of worst violator
242 :
243 : logical :: found ! true => at least 1 minimum violator found
244 :
245 1486346256 : real(r8) :: badvals(ncold, lver) ! Collector for outfld calls
246 1486346256 : real(r8) :: badcols(ncold) ! Column sum for outfld
247 : real(r8) :: worst ! biggest violator
248 : !
249 : !-----------------------------------------------------------------------
250 : !
251 :
252 743173128 : call t_startf ('qneg3')
253 : ! The first time we call this, we need to determine whether to call outfld
254 743173128 : if (.not. cnst_out_calc) then
255 1536 : call calc_cnst_out()
256 : end if
257 :
258 743173128 : if (collect_stats) then
259 743173128 : index = find_index3(trim(subnam))
260 : else
261 : index = -1
262 : end if
263 :
264 1605975696 : do m = lconst_beg, lconst_end
265 : nvals = 0
266 81103441392 : found = .false.
267 81103441392 : worst = worst_reset
268 >13649*10^8 : badvals(:,:) = 0.0_r8
269 81103441392 : iw = -1
270 81103441392 : kw = -1
271 : !
272 : ! Test all field values for being less than minimum value. Set q = qmin
273 : ! for all such points. Trace offenders and identify worst one.
274 : !
275 81103441392 : do k = 1, lver
276 >13406*10^8 : do i = 1, ncol
277 >13398*10^8 : if (q(i,k,m) < qmin(m)) then
278 7550970901 : found = .true.
279 : nvals = nvals + 1
280 7550970901 : badvals(i, k) = q(i, k, m)
281 7550970901 : if (index > 0) then
282 7550970901 : qneg3_warn_num(m, index) = qneg3_warn_num(m, index) + 1
283 : end if
284 7550970901 : if (q(i,k,m) < worst) then
285 486955362 : worst = q(i,k,m)
286 486955362 : iw = i
287 486955362 : kw = k
288 486955362 : if (index > 0) then
289 486955362 : qneg3_warn_worst(m, index) = worst
290 : end if
291 : end if
292 7550970901 : q(i,k,m) = qmin(m)
293 : end if
294 : end do
295 : end do
296 : ! Maybe output bad values
297 862802568 : if ((cnst_outfld(m)) .and. (worst < worst_reset)) then
298 0 : call outfld(trim(diag_names(m)), badvals, pcols, idx)
299 : end if
300 1605975696 : if ((cnst_outfld(pcnst+m)) .and. (worst < worst_reset)) then
301 0 : do i = 1, pcols
302 0 : badcols(i) = SUM(badvals(i,:))
303 : end do
304 0 : call outfld(trim(diag_names(pcnst+m)), badcols, pcols, idx)
305 : end if
306 : end do
307 743173128 : call t_stopf ('qneg3')
308 :
309 743173128 : end subroutine qneg3
310 :
311 1489176 : subroutine qneg4 (subnam, lchnk, ncol, ztodt, &
312 1489176 : qbot, srfrpdel, shflx, lhflx, qflx)
313 : !-----------------------------------------------------------------------
314 : !
315 : ! Purpose:
316 : ! Check if moisture flux into the ground is exceeding the total
317 : ! moisture content of the lowest model layer (creating negative moisture
318 : ! values). If so, then subtract the excess from the moisture and
319 : ! latent heat fluxes and add it to the sensible heat flux.
320 : !
321 : ! Method:
322 : ! <Describe the algorithm(s) used in the routine.>
323 : ! <Also include any applicable external references.>
324 : !
325 : ! Author: J. Olson
326 : !
327 : !-----------------------------------------------------------------------
328 743173128 : use physconst, only: gravit, latvap
329 : use constituents, only: qmin
330 : use cam_history, only: outfld
331 :
332 : !
333 : ! Input arguments
334 : !
335 : character(len=*), intent(in) :: subnam ! name of calling routine
336 : !
337 : integer, intent(in) :: lchnk ! chunk index
338 : integer, intent(in) :: ncol ! number of atmospheric columns
339 : !
340 : real(r8), intent(in) :: ztodt ! two times model timestep (2 delta-t)
341 : real(r8), intent(in) :: qbot(ncol,pcnst) ! moisture at lowest model level
342 : real(r8), intent(in) :: srfrpdel(ncol) ! 1./(pint(K+1)-pint(K))
343 : !
344 : ! Input/Output arguments
345 : !
346 : real(r8), intent(inout) :: shflx(ncol) ! Surface sensible heat flux (J/m2/s)
347 : real(r8), intent(inout) :: lhflx(ncol) ! Surface latent heat flux (J/m2/s)
348 : real(r8), intent(inout) :: qflx (ncol,pcnst) ! surface water flux (kg/m^2/s)
349 : !
350 : !---------------------------Local workspace-----------------------------
351 : !
352 : integer :: i ! column index
353 : integer :: iw ! i index of worst violator
354 : integer :: index ! caller bin index
355 : !
356 : real(r8):: worst ! biggest violator
357 2978352 : real(r8):: excess(ncol) ! Excess downward sfc latent heat flux
358 : !
359 : !-----------------------------------------------------------------------
360 :
361 1489176 : call t_startf ('qneg4')
362 : ! The first time we call this, we need to determine whether to call outfld
363 1489176 : if (.not. cnst_out_calc) then
364 0 : call calc_cnst_out()
365 : end if
366 :
367 1489176 : if (collect_stats) then
368 1489176 : index = find_index4(trim(subnam))
369 : else
370 : index = -1
371 : end if
372 :
373 : !
374 : ! Compute excess downward (negative) q flux compared to a theoretical
375 : ! maximum downward q flux. The theoretical max is based upon the
376 : ! given moisture content of lowest level of the model atmosphere.
377 : !
378 1489176 : worst = worst_reset
379 24865776 : do i = 1, ncol
380 23376600 : excess(i) = qflx(i,1) - (qmin(1) - qbot(i,1))/(ztodt*gravit*srfrpdel(i))
381 : !
382 : ! If there is an excess downward (negative) q flux, then subtract
383 : ! excess from "qflx" and "lhflx" and add to "shflx".
384 : !
385 24865776 : if (excess(i) < 0._r8) then
386 2 : if (excess(i) < worst) then
387 2 : iw = i
388 2 : worst = excess(i)
389 : end if
390 2 : qflx (i,1) = qflx (i,1) - excess(i)
391 2 : lhflx(i) = lhflx(i) - excess(i)*latvap
392 2 : shflx(i) = shflx(i) + excess(i)*latvap
393 2 : if (index > 0) then
394 2 : qneg4_warn_num(index) = qneg4_warn_num(index) + 1
395 : end if
396 : end if
397 : end do
398 : ! Maybe output bad values
399 1489176 : if ((cnst_outfld((2*pcnst)+1)) .and. (worst < worst_reset)) then
400 0 : do i = 1, ncol
401 0 : if (excess(i) > 0.0_r8) then
402 0 : excess(i) = 0.0_r8
403 : end if
404 : end do
405 0 : call outfld(trim(diag_names((2*pcnst)+1)), excess(1:ncol), ncol, lchnk)
406 : end if
407 1489176 : call t_stopf ('qneg4')
408 :
409 1489176 : end subroutine qneg4
410 :
411 369408 : subroutine qneg_print_summary(end_of_run)
412 1489176 : use spmd_utils, only: mpicom, masterprocid, masterproc
413 : use spmd_utils, only: MPI_MIN, MPI_SUM, MPI_INTEGER, MPI_REAL8
414 :
415 : logical, intent(in) :: end_of_run
416 :
417 : integer :: global_warn_num(pcnst)
418 : real(r8) :: global_warn_worst(pcnst)
419 : integer :: index, m
420 : integer :: ierr
421 :
422 369408 : if (collect_stats) then
423 369408 : if (timestep_reset .or. end_of_run) then
424 38400 : do index = 1, num3_bins
425 : ! QNEG3
426 36864 : call reset_stats(global_warn_num(:), global_warn_worst(:))
427 : call MPI_REDUCE(qneg3_warn_num(:, index), global_warn_num(:), &
428 36864 : pcnst, MPI_INTEGER, MPI_SUM, masterprocid, mpicom, ierr)
429 : call MPI_REDUCE(qneg3_warn_worst(:, index), global_warn_worst(:),&
430 36864 : pcnst, MPI_REAL8, MPI_MIN, masterprocid, mpicom, ierr)
431 36864 : if (masterproc) then
432 2016 : do m = 1, pcnst
433 1968 : if ( (global_warn_num(m) > 0) .and. &
434 : (abs(global_warn_worst(m)) > tol)) then
435 67 : write(iulog, 9100) trim(qneg3_warn_labels(index)), &
436 67 : trim(cnst_name(m)), global_warn_num(m), &
437 134 : global_warn_worst(m)
438 : end if
439 2016 : call shr_sys_flush(iulog)
440 : end do
441 : end if
442 38400 : call reset_stats(qneg3_warn_num(:,index), qneg3_warn_worst(:,index))
443 : end do
444 7680 : do index = 1, num4_bins
445 : ! QNEG4
446 6144 : call reset_stats(qneg4_warn_num(index), qneg4_warn_worst(index))
447 6144 : call reset_stats(global_warn_num(1), global_warn_worst(1))
448 : call MPI_REDUCE(qneg4_warn_num(index), global_warn_num(1), &
449 6144 : 1, MPI_INTEGER, MPI_SUM, masterprocid, mpicom, ierr)
450 : call MPI_REDUCE(qneg4_warn_worst(index), global_warn_worst(1), &
451 6144 : 1, MPI_REAL8, MPI_MIN, masterprocid, mpicom, ierr)
452 6144 : if (masterproc) then
453 8 : if ( (global_warn_num(1) > 0) .and. &
454 : (abs(global_warn_worst(1)) > tol)) then
455 0 : write(iulog, 9101) trim(qneg4_warn_labels(index)), &
456 0 : global_warn_num(1), global_warn_worst(1)
457 : end if
458 8 : call shr_sys_flush(iulog)
459 : end if
460 13824 : call reset_stats(qneg4_warn_num(index), qneg4_warn_worst(index))
461 : end do
462 : end if
463 : end if
464 :
465 369408 : return
466 : 9100 format(' QNEG3 from ', a, ':', a, &
467 : ' Min. mixing ratio violated at ', i9, ' points. Worst = ', e10.1)
468 : 9101 format(' QNEG4 from ',a,': moisture flux exceeded at', &
469 : i9, ' points. Worst = ', e10.1)
470 : end subroutine qneg_print_summary
471 :
472 104448 : subroutine reset_stats_array(num_array, worst_array)
473 : ! Private routine to reset statistics
474 : integer, intent(inout) :: num_array(:)
475 : real(r8), intent(inout) :: worst_array(:)
476 :
477 4386816 : num_array(:) = 0
478 4386816 : worst_array(:) = worst_reset
479 104448 : end subroutine reset_stats_array
480 :
481 19968 : subroutine reset_stats_scalar(num, worst)
482 : ! Private routine to reset statistics
483 : integer, intent(inout) :: num
484 : real(r8), intent(inout) :: worst
485 :
486 19968 : num = 0
487 19968 : worst = worst_reset
488 19968 : end subroutine reset_stats_scalar
489 :
490 : end module qneg_module
|
