Line data Source code
1 : !-----------------------------------------------------------------------
2 : ! $Id: fill_holes.F90 8738 2018-07-19 19:58:53Z bmg2@uwm.edu $
3 : !===============================================================================
4 : module fill_holes
5 :
6 : implicit none
7 :
8 : public :: fill_holes_driver, &
9 : fill_holes_vertical, &
10 : hole_filling_hm_one_lev, &
11 : fill_holes_hydromet, &
12 : fill_holes_wv, &
13 : clip_hydromet_conc_mvr, &
14 : setup_stats_indices
15 :
16 : private ! Set Default Scope
17 :
18 : contains
19 :
20 : !=============================================================================
21 178701120 : subroutine fill_holes_vertical( nz, ngrdcol, threshold, upper_hf_level, &
22 178701120 : dz, rho_ds, &
23 178701120 : field )
24 :
25 : ! Description:
26 : ! This subroutine clips values of 'field' that are below 'threshold' as much
27 : ! as possible (i.e. "fills holes"), but conserves the total integrated mass
28 : ! of 'field'. This prevents clipping from acting as a spurious source.
29 : !
30 : ! Mass is conserved by reducing the clipped field everywhere by a constant
31 : ! multiplicative coefficient.
32 : !
33 : ! This subroutine does not guarantee that the clipped field will exceed
34 : ! threshold everywhere; blunt clipping is needed for that.
35 : !
36 : ! The lowest level (k=1) should not be included, as the hole-filling scheme
37 : ! should not alter the set value of 'field' at the surface (for momentum
38 : ! level variables), or consider the value of 'field' at a level below the
39 : ! surface (for thermodynamic level variables).
40 : !
41 : ! For momentum level variables only, the hole-filling scheme should not
42 : ! alter the set value of 'field' at the upper boundary level (k=nz).
43 : ! So for momemtum level variables, call with upper_hf_level=nz-1, and
44 : ! for thermodynamic level variables, call with upper_hf_level=nz.
45 : !
46 : ! References:
47 : ! ``Numerical Methods for Wave Equations in Geophysical Fluid
48 : ! Dynamics'', Durran (1999), p. 292.
49 : !-----------------------------------------------------------------------
50 :
51 : use grid_class, only: &
52 : grid ! Type
53 :
54 : use clubb_precision, only: &
55 : core_rknd ! Variable(s)
56 :
57 : use constants_clubb, only: &
58 : eps, &
59 : one, &
60 : num_hf_draw_points ! The number of points on either side of the hole;
61 : ! Mass is drawn from these points to fill the hole
62 :
63 : implicit none
64 :
65 : ! --------------------- Input variables ---------------------
66 : integer, intent(in) :: &
67 : nz, &
68 : ngrdcol
69 :
70 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
71 : dz ! Spacing between thermodynamic grid levels; centered over
72 : ! momentum grid levels
73 : ! OR
74 : ! Spcaing between momentum grid levels; centered over
75 : ! thermodynamic grid levels
76 :
77 : integer, intent(in) :: &
78 : upper_hf_level ! Upper grid level of global hole-filling range []
79 :
80 : real( kind = core_rknd ), intent(in) :: &
81 : threshold ! A threshold (e.g. w_tol*w_tol) below which field must not
82 : ! fall [Units vary; same as field]
83 :
84 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
85 : rho_ds ! Dry, static density on thermodynamic or momentum levels [kg/m^3]
86 :
87 : ! --------------------- Input/Output variable ---------------------
88 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(inout) :: &
89 : field ! The field (e.g. wp2) that contains holes [Units same as threshold]
90 :
91 : ! --------------------- Local Variables ---------------------
92 : integer :: &
93 : i, & ! Loop index for column []
94 : k, & ! Loop index for absolute grid level []
95 : j, &
96 : k_start, & ! Lower grid level of local hole-filling range []
97 : k_end ! Upper grid level of local hole-filling range []
98 :
99 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
100 357402240 : rho_ds_dz, & ! rho_ds * dz
101 357402240 : invrs_denom_integral, & ! Inverse of the integral in the denominator (see description)
102 357402240 : field_clipped ! The raw field (e.g. wp2) that contains no holes
103 : ! [Units same as threshold]
104 :
105 : real( kind = core_rknd ) :: &
106 : field_avg, & ! Vertical average of field [Units of field]
107 : field_clipped_avg, & ! Vertical average of clipped field [Units of field]
108 : mass_fraction ! Coefficient that multiplies clipped field
109 : ! in order to conserve mass. []
110 :
111 : real( kind = core_rknd ), dimension(ngrdcol) :: &
112 357402240 : denom_integral_global, & ! Integral in the denominator for global filling
113 357402240 : numer_integral_global, & ! Integral in the numerator for global filling
114 357402240 : field_avg_global, & ! Vertical average of field [Units of field]
115 357402240 : mass_fraction_global ! Coefficient that multiplies clipped field
116 : ! in order to conserve mass. []
117 :
118 : logical :: &
119 : l_field_below_threshold
120 :
121 : real( kind = core_rknd ) :: &
122 : rho_k_sum
123 :
124 : ! --------------------- Begin Code ---------------------
125 :
126 178701120 : l_field_below_threshold = .false.
127 :
128 : !$acc parallel loop gang vector collapse(2) default(present) &
129 : !$acc reduction(.or.:l_field_below_threshold)
130 15368296320 : do k = 1, nz
131 >25380*10^7 : do i = 1, ngrdcol
132 >25363*10^7 : if ( field(i,k) < threshold ) then
133 8966831621 : l_field_below_threshold = .true.
134 : end if
135 : end do
136 : end do
137 : !$acc end parallel loop
138 :
139 : ! If all field values are above the specified threshold, no hole filling is required
140 178701120 : if ( .not. l_field_below_threshold ) then
141 : return
142 : end if
143 :
144 : !$acc enter data create( invrs_denom_integral, field_clipped, denom_integral_global, rho_ds_dz, &
145 : !$acc numer_integral_global, field_avg_global, mass_fraction_global )
146 :
147 : !$acc parallel loop gang vector collapse(2) default(present)
148 2380439580 : do k = 1, nz
149 39324789090 : do i = 1, ngrdcol
150 39297109560 : rho_ds_dz(i,k) = rho_ds(i,k) * dz(i,k)
151 : end do
152 : end do
153 : !$acc end parallel loop
154 :
155 :
156 : ! denom_integral does not change throughout the hole filling algorithm
157 : ! so we can calculate it before hand. This results in unneccesary computations,
158 : ! but is parallelizable and reduces the cost of the serial k loop
159 : !$acc parallel loop gang vector collapse(2) default(present)
160 462318936 : do i = 1, ngrdcol
161 34805633527 : do k = 2+num_hf_draw_points, upper_hf_level-num_hf_draw_points
162 :
163 : ! This loop and division could be written more compactly as
164 : ! invrs_denom_integral(i,k) = one / sum(rho_ds_dz(i,k-num_hf_draw_points:k+num_hf_draw_points))
165 : ! but has been manually written in loop form to improve performance
166 : ! when using OpenMP target offloading
167 : ! See: https://github.com/larson-group/clubb/issues/1138#issuecomment-1974918151
168 34343314591 : rho_k_sum = 0.0_core_rknd
169 :
170 >20605*10^7 : do j = k - num_hf_draw_points, k + num_hf_draw_points
171 >20605*10^7 : rho_k_sum = rho_k_sum + rho_ds_dz(i,j)
172 : end do
173 :
174 34777953997 : invrs_denom_integral(i,k) = one / rho_k_sum
175 : end do
176 : end do
177 : !$acc end parallel loop
178 :
179 : !$acc parallel loop gang vector default(present)
180 462318936 : do i = 1, ngrdcol
181 :
182 : ! Make one pass up the profile, filling holes as much as we can using
183 : ! nearby mass, ignoring the first level.
184 : !
185 : ! This loop can only be done in serial due to the field values required for the next
186 : ! iteration potentially changing. We could in theory expose more parallelism in cases
187 : ! where there are large enough gaps between vertical levels which need hole-filling,
188 : ! but levels which require hole-filling are often close or consecutive.
189 34805633527 : do k = 2+num_hf_draw_points, upper_hf_level-num_hf_draw_points
190 :
191 34343314591 : k_start = k - num_hf_draw_points
192 34343314591 : k_end = k + num_hf_draw_points
193 :
194 >16718*10^7 : if ( any( field(i,k_start:k_end) < threshold ) ) then
195 :
196 : ! Compute the field's vertical average cenetered at k, which we must conserve,
197 : ! see description of the vertical_avg function in advance_helper_module
198 : field_avg = sum( rho_ds_dz(i,k_start:k_end) * field(i,k_start:k_end) ) &
199 55769172702 : * invrs_denom_integral(i,k)
200 :
201 : ! Clip small or negative values from field.
202 9294862117 : if ( field_avg >= threshold ) then
203 : ! We know we can fill in holes completely
204 7520717778 : field_clipped(i,k_start:k_end) = max( threshold, field(i,k_start:k_end) )
205 : else
206 : ! We can only fill in holes partly;
207 : ! to do so, we remove all mass above threshold.
208 48248454924 : field_clipped(i,k_start:k_end) = min( threshold, field(i,k_start:k_end) )
209 : endif
210 :
211 : ! Compute the clipped field's vertical integral.
212 : ! clipped_total_mass >= original_total_mass,
213 : ! see description of the vertical_avg function in advance_helper_module
214 : field_clipped_avg = sum( rho_ds_dz(i,k_start:k_end) * field_clipped(i,k_start:k_end) ) &
215 55769172702 : * invrs_denom_integral(i,k)
216 :
217 : ! Avoid divide by zero issues by doing nothing if field_clipped_avg ~= threshold
218 9294862117 : if ( abs(field_clipped_avg-threshold) > abs(field_clipped_avg+threshold)*eps/2) then
219 : ! Compute coefficient that makes the clipped field have the same mass as the
220 : ! original field. We should always have mass_fraction > 0.
221 : mass_fraction = ( field_avg - threshold ) &
222 9294862046 : / ( field_clipped_avg - threshold )
223 :
224 : ! Calculate normalized, filled field
225 : field(i,k_start:k_end) = threshold &
226 55769172276 : + mass_fraction * ( field_clipped(i,k_start:k_end) - threshold )
227 : endif
228 :
229 : endif
230 :
231 : end do
232 :
233 : end do
234 : !$acc end parallel loop
235 :
236 : l_field_below_threshold = .false.
237 :
238 : !$acc parallel loop gang vector collapse(2) default(present) reduction(.or.:l_field_below_threshold)
239 2380439580 : do k = 1, nz
240 39324789090 : do i = 1, ngrdcol
241 39297109560 : if ( field(i,k) < threshold ) then
242 7580956695 : l_field_below_threshold = .true.
243 : end if
244 : end do
245 : end do
246 : !$acc end parallel loop
247 :
248 : ! If all field values are above the threshold, no further hole filling is required
249 27679530 : if ( .not. l_field_below_threshold ) then
250 : !$acc exit data delete( invrs_denom_integral, field_clipped, denom_integral_global, rho_ds_dz, &
251 : !$acc numer_integral_global, field_avg_global, mass_fraction_global )
252 : return
253 : end if
254 :
255 :
256 : ! Now we fill holes globally to maximize the chance that all holes are filled.
257 : ! To improve parallelism we assume that global hole filling needs to be done
258 : ! for each grid column, perform all calculations required, and only check
259 : ! if any holes need filling before the final step of updating the field.
260 :
261 : ! Compute the numerator and denominator integrals
262 : !$acc parallel loop gang vector default(present)
263 187199796 : do i = 1, ngrdcol
264 175993909 : numer_integral_global(i) = 0.0_core_rknd
265 187199796 : denom_integral_global(i) = 0.0_core_rknd
266 : end do
267 : !$acc end parallel loop
268 :
269 : !$acc parallel loop gang vector default(present)
270 187199796 : do i = 1, ngrdcol
271 14794806837 : do k = 2, upper_hf_level
272 14607607041 : numer_integral_global(i) = numer_integral_global(i) + rho_ds_dz(i,k) * field(i,k)
273 :
274 14783600950 : denom_integral_global(i) = denom_integral_global(i) + rho_ds_dz(i,k)
275 : end do
276 : end do
277 : !$acc end parallel loop
278 :
279 :
280 : !$acc parallel loop gang vector default(present)
281 187199796 : do i = 1, ngrdcol
282 :
283 : ! Find the vertical average of field, using the precomputed numerator and denominator,
284 : ! see description of the vertical_avg function in advance_helper_module
285 175993909 : field_avg_global(i) = numer_integral_global(i) / denom_integral_global(i)
286 :
287 : ! Clip small or negative values from field.
288 187199796 : if ( field_avg_global(i) >= threshold ) then
289 : ! We know we can fill in holes completely
290 14741288806 : field_clipped(i,2:upper_hf_level) = max( threshold, field(i,2:upper_hf_level) )
291 : else
292 : ! We can only fill in holes partly;
293 : ! to do so, we remove all mass above threshold.
294 42312144 : field_clipped(i,2:upper_hf_level) = min( threshold, field(i,2:upper_hf_level) )
295 : endif
296 :
297 : end do
298 : !$acc end parallel loop
299 :
300 : ! To compute the clipped field's vertical integral we only need to recompute the numerator
301 : !$acc parallel loop gang vector default(present)
302 187199796 : do i = 1, ngrdcol
303 175993909 : numer_integral_global(i) = 0.0_core_rknd
304 14794806837 : do k = 2, upper_hf_level
305 14783600950 : numer_integral_global(i) = numer_integral_global(i) + rho_ds_dz(i,k) * field_clipped(i,k)
306 : end do
307 : end do
308 : !$acc end parallel loop
309 :
310 : !$acc parallel loop gang vector default(present)
311 187199796 : do i = 1, ngrdcol
312 :
313 : ! Do not complete calculations or update field values for this
314 : ! column if there are no holes that need filling
315 5908083612 : if ( any( field(i,2:upper_hf_level) < threshold ) ) then
316 :
317 : ! Compute the clipped field's vertical integral,
318 : ! see description of the vertical_avg function in advance_helper_module
319 145463245 : field_clipped_avg = numer_integral_global(i) / denom_integral_global(i)
320 :
321 : ! Avoid divide by zero issues by doing nothing if field_clipped_avg ~= threshold
322 145463245 : if ( abs(field_clipped_avg-threshold) > abs(field_clipped_avg+threshold)*eps/2) then
323 :
324 : ! Compute coefficient that makes the clipped field have the same mass as the
325 : ! original field. We should always have mass_fraction > 0.
326 : mass_fraction_global(i) = ( field_avg_global(i) - threshold ) &
327 145463245 : / ( field_clipped_avg - threshold )
328 :
329 : ! Calculate normalized, filled field
330 : field(i,2:upper_hf_level) = threshold + mass_fraction_global(i) &
331 : * ( field_clipped(i,2:upper_hf_level) &
332 12218912580 : - threshold )
333 : end if
334 : end if
335 :
336 : end do
337 : !$acc end parallel loop
338 :
339 : !$acc exit data delete( invrs_denom_integral, field_clipped, denom_integral_global, rho_ds_dz, &
340 : !$acc numer_integral_global, field_avg_global, mass_fraction_global )
341 :
342 : return
343 :
344 : end subroutine fill_holes_vertical
345 :
346 : !===============================================================================
347 0 : subroutine hole_filling_hm_one_lev( num_hm_fill, hm_one_lev, & ! Intent(in)
348 0 : hm_one_lev_filled ) ! Intent(out)
349 :
350 : ! Description:
351 : ! Fills holes between same-phase (i.e. either liquid or frozen) hydrometeors for
352 : ! one height level.
353 : !
354 : ! Warning: Do not input hydrometeors of different phases, e.g. liquid and frozen.
355 : ! Otherwise heat will not be conserved.
356 : !
357 : ! References:
358 : !
359 : ! None
360 : !-----------------------------------------------------------------------
361 :
362 : use constants_clubb, only: &
363 : one, & ! Variable(s)
364 : zero, &
365 : eps
366 :
367 : use clubb_precision, only: &
368 : core_rknd ! Variable(s)
369 :
370 : use error_code, only: &
371 : clubb_at_least_debug_level ! Procedure
372 :
373 : implicit none
374 :
375 : ! Input Variables
376 : integer, intent(in) :: num_hm_fill ! number of hydrometeors involved
377 :
378 : real(kind = core_rknd), dimension(num_hm_fill), intent(in) :: hm_one_lev
379 :
380 : ! Output Variables
381 : real(kind = core_rknd), dimension(num_hm_fill), intent(out) :: hm_one_lev_filled
382 :
383 : ! Local Variables
384 : integer :: num_neg_hm ! number of holes
385 :
386 : real(kind = core_rknd) :: &
387 : total_hole, & ! Size of the hole ( missing mass, less than 0 )
388 : total_mass ! Total mass to fill the hole
389 : ! total mass of water substance = total_mass + total_hole
390 :
391 : integer :: i ! loop iterator
392 :
393 : !-----------------------------------------------------------------------
394 :
395 : !----- Begin Code -----
396 :
397 : ! Initialization
398 0 : hm_one_lev_filled = 0._core_rknd
399 : total_hole = 0._core_rknd
400 : total_mass = 0._core_rknd
401 : num_neg_hm = 0
402 :
403 : ! Determine the total size of the hole and the number of neg. hydrometeors
404 : ! and the total mass of hole filling material
405 0 : do i=1, num_hm_fill
406 : ! print *, "hm_one_lev(",i,") = ", hm_one_lev(i)
407 0 : if ( hm_one_lev(i) < zero ) then
408 0 : total_hole = total_hole + hm_one_lev(i) ! less than zero
409 : num_neg_hm = num_neg_hm + 1
410 : else
411 0 : total_mass = total_mass + hm_one_lev(i)
412 : endif
413 :
414 : enddo
415 :
416 : ! print *, "total_hole = ", total_hole
417 : ! print *, "total_mass = ", total_mass
418 : ! print *, "num_neg_hm = ", num_neg_hm
419 :
420 : ! There is no water substance at all to fill the hole
421 0 : if ( abs(total_mass) < eps ) then
422 :
423 0 : if ( clubb_at_least_debug_level( 2 ) ) then
424 0 : print *, "Warning: One-level hole filling was not successful! total_mass ~= 0"
425 : endif
426 :
427 0 : hm_one_lev_filled = hm_one_lev
428 :
429 : return
430 : endif
431 :
432 : ! Fill the holes and adjust the remaining quantities:
433 : ! hm_filled(i) = 0, if hm(i) < 0
434 : ! or
435 : ! hm_filled(i) = (1 + total_hole/total_mass)*hm(i), if hm(i) > 0
436 0 : do i=1, num_hm_fill
437 :
438 : ! if there is not enough material, fill the holes partially with all the material available
439 0 : if ( abs(total_hole) > total_mass ) then
440 :
441 0 : if ( clubb_at_least_debug_level( 2 ) ) then
442 : print *, "Warning: One-level hole filling was not able to fill holes completely!" // &
443 0 : " The holes were filled partially. |total_hole| > total_mass"
444 : endif
445 :
446 0 : hm_one_lev_filled(i) = min(hm_one_lev(i), zero) * ( one + total_mass / total_hole )
447 :
448 : else ! fill holes completely
449 0 : hm_one_lev_filled(i) = max(hm_one_lev(i), zero) * ( one + total_hole / total_mass )
450 :
451 : endif
452 :
453 : enddo
454 :
455 : ! Assertion checks (water substance conservation, non-negativity)
456 0 : if ( clubb_at_least_debug_level( 2 ) ) then
457 :
458 0 : if ( abs(sum( hm_one_lev ) - sum(hm_one_lev_filled)) > &
459 : abs(sum( hm_one_lev ) + sum(hm_one_lev_filled)) * eps/2 ) then
460 0 : print *, "Warning: Hole filling was not conservative!"
461 : endif
462 :
463 0 : if ( any( hm_one_lev_filled < zero ) ) then
464 0 : print *, "Warning: Hole filling failed! A hole could not be filled."
465 : endif
466 :
467 : endif
468 :
469 : return
470 :
471 : end subroutine hole_filling_hm_one_lev
472 : !-----------------------------------------------------------------------
473 :
474 : !-----------------------------------------------------------------------
475 0 : subroutine fill_holes_hydromet( nz, hydromet_dim, hydromet, & ! Intent(in)
476 0 : l_frozen_hm, l_mix_rat_hm, & ! Intent(in)
477 0 : hydromet_filled ) ! Intent(out)
478 :
479 : ! Description:
480 : ! Fills holes between same-phase hydrometeors(i.e. for frozen hydrometeors).
481 : ! The hole filling conserves water substance between all same-phase (frozen or liquid)
482 : ! hydrometeors at each height level.
483 : !
484 : ! Attention: The hole filling for the liquid phase hydrometeors is not yet implemented
485 : !
486 : ! Attention: l_frozen_hm and l_mix_rat_hm need to be set up before this subroutine is called!
487 : !
488 : ! References:
489 : !
490 : ! None
491 : !-----------------------------------------------------------------------
492 :
493 : use clubb_precision, only: &
494 : core_rknd
495 :
496 : use constants_clubb, only: &
497 : zero
498 :
499 : implicit none
500 :
501 : ! Input Variables
502 : integer, intent(in) :: hydromet_dim, nz
503 :
504 : real( kind = core_rknd ), dimension(nz,hydromet_dim), intent(in) :: &
505 : hydromet
506 :
507 : logical, dimension(hydromet_dim), intent(in) :: &
508 : l_frozen_hm, & ! if true, then the hydrometeor is frozen; otherwise liquid
509 : l_mix_rat_hm ! if true, then the quantity is a hydrometeor mixing ratio
510 :
511 : ! Output Variables
512 : real( kind = core_rknd ), dimension(nz,hydromet_dim), intent(out) :: &
513 : hydromet_filled
514 :
515 : ! Local Variables
516 : integer :: i,j ! Loop iterators
517 :
518 : integer :: num_frozen_hm ! Number of frozen hydrometeor mixing ratios
519 :
520 : real( kind = core_rknd ), dimension(:,:), allocatable :: &
521 0 : hydromet_frozen, & ! Frozen hydrometeor mixing ratios
522 0 : hydromet_frozen_filled ! Frozen hydrometeor mixing ratios after hole filling
523 :
524 : !-----------------------------------------------------------------------
525 :
526 : !----- Begin Code -----
527 :
528 : ! Determine the number of frozen hydrometeor mixing ratios
529 0 : num_frozen_hm = 0
530 0 : do i=1,hydromet_dim
531 0 : if ( l_frozen_hm(i) .and. l_mix_rat_hm(i) ) then
532 0 : num_frozen_hm = num_frozen_hm + 1
533 : endif
534 : enddo
535 :
536 : ! Allocation
537 0 : allocate( hydromet_frozen(nz,num_frozen_hm) )
538 0 : allocate( hydromet_frozen_filled(nz,num_frozen_hm) )
539 :
540 : ! Determine frozen hydrometeor mixing ratios
541 0 : j = 1
542 0 : do i = 1,hydromet_dim
543 0 : if ( l_frozen_hm(i) .and. l_mix_rat_hm(i) ) then
544 0 : hydromet_frozen(:,j) = hydromet(:,i)
545 0 : j = j+1
546 : endif
547 : enddo
548 :
549 : ! Fill holes for the frozen hydrometeors
550 0 : do i=1,nz
551 0 : if ( any( hydromet_frozen(i,:) < zero ) ) then
552 : call hole_filling_hm_one_lev( num_frozen_hm, hydromet_frozen(i,:), & ! Intent(in)
553 0 : hydromet_frozen_filled(i,:) ) ! Intent(out)
554 : else
555 0 : hydromet_frozen_filled(i,:) = hydromet_frozen(i,:)
556 : endif
557 : enddo
558 :
559 : ! Setup the filled hydromet array
560 : j = 1
561 0 : do i=1, hydromet_dim
562 0 : if ( l_frozen_hm(i) .and. l_mix_rat_hm(i) ) then
563 0 : hydromet_filled(:,i) = hydromet_frozen_filled(:,j)
564 0 : j = j+1
565 : else
566 0 : hydromet_filled(:,i) = hydromet(:,i)
567 : endif
568 : enddo
569 :
570 : !!! Here we could do the same hole filling for all the liquid phase hydrometeors
571 :
572 0 : return
573 0 : end subroutine fill_holes_hydromet
574 : !-----------------------------------------------------------------------
575 :
576 : !-----------------------------------------------------------------------
577 0 : subroutine fill_holes_wv( nz, dt, exner, hydromet_name, & ! Intent(in)
578 0 : rvm_mc, thlm_mc, hydromet )! Intent(inout)
579 :
580 : ! Description:
581 : ! Fills holes using the cloud water mixing ratio from the current height level.
582 : !
583 : ! References:
584 : !
585 : ! None
586 : !-----------------------------------------------------------------------
587 :
588 : use clubb_precision, only: &
589 : core_rknd
590 :
591 : use constants_clubb, only: &
592 : zero_threshold, &
593 : Lv, &
594 : Ls, &
595 : Cp
596 :
597 : implicit none
598 :
599 : ! Input Variables
600 : integer, intent(in) :: nz
601 :
602 : real( kind = core_rknd ), intent(in) :: &
603 : dt ! Timestep [s]
604 :
605 : character(len=10), intent(in) :: hydromet_name
606 :
607 : real( kind = core_rknd ), dimension(nz), intent(in) :: &
608 : exner ! Exner function [-]
609 :
610 : ! Input/Output Variables
611 : real( kind = core_rknd ), dimension(nz), intent(inout) :: &
612 : hydromet, & ! Hydrometeor array [units vary]
613 : rvm_mc, &
614 : thlm_mc
615 :
616 : ! Local Variables
617 : integer :: k ! Loop iterator
618 :
619 : real( kind = core_rknd ) :: rvm_clip_tndcy
620 : !-----------------------------------------------------------------------
621 :
622 : !----- Begin Code -----
623 :
624 0 : do k = 2, nz, 1
625 :
626 0 : if ( hydromet(k) < zero_threshold ) then
627 :
628 : ! Set rvm_clip_tndcy to the time tendency applied to vapor and removed
629 : ! from the hydrometeor.
630 0 : rvm_clip_tndcy = hydromet(k) / dt
631 :
632 : ! Adjust the tendency rvm_mc accordingly
633 0 : rvm_mc(k) = rvm_mc(k) + rvm_clip_tndcy
634 :
635 : ! Adjust the tendency of thlm_mc according to whether the
636 : ! effect is an evaporation or sublimation tendency.
637 0 : select case ( trim( hydromet_name ) )
638 : case( "rrm" )
639 0 : thlm_mc(k) = thlm_mc(k) - rvm_clip_tndcy * ( Lv / ( Cp*exner(k) ) )
640 : case( "rim", "rsm", "rgm" )
641 0 : thlm_mc(k) = thlm_mc(k) - rvm_clip_tndcy * ( Ls / ( Cp*exner(k) ) )
642 : case default
643 0 : error stop "Fatal error in microphys_driver"
644 : end select
645 :
646 : ! Set the mixing ratio to 0
647 0 : hydromet(k) = zero_threshold
648 :
649 : endif ! hydromet(k,i) < 0
650 :
651 : enddo ! k = 2..gr%nz
652 :
653 0 : return
654 : end subroutine fill_holes_wv
655 : !-----------------------------------------------------------------------
656 :
657 : !-----------------------------------------------------------------------
658 0 : subroutine fill_holes_driver( gr, nz, dt, hydromet_dim, hm_metadata, & ! Intent(in)
659 : l_fill_holes_hm, & ! Intent(in)
660 0 : rho_ds_zm, rho_ds_zt, exner, & ! Intent(in)
661 : stats_metadata, & ! Intent(in)
662 : stats_zt, & ! intent(inout)
663 0 : thlm_mc, rvm_mc, hydromet ) ! Intent(inout)
664 :
665 : ! Description:
666 : ! Fills holes between same-phase hydrometeors(i.e. for frozen hydrometeors).
667 : ! The hole filling conserves water substance between all same-phase (frozen or liquid)
668 : ! hydrometeors at each height level.
669 : !
670 : ! Attention: The hole filling for the liquid phase hydrometeors is not yet implemented
671 : !
672 : ! Attention: l_frozen_hm and l_mix_rat_hm need to be set up before this subroutine is called!
673 : !
674 : ! References:
675 : !
676 : ! None
677 : !-----------------------------------------------------------------------
678 :
679 : use grid_class, only: &
680 : grid ! Type
681 :
682 : use clubb_precision, only: &
683 : core_rknd ! Variable(s)
684 :
685 : use constants_clubb, only: &
686 : zero, &
687 : zero_threshold, &
688 : Lv, &
689 : Ls, &
690 : Cp, &
691 : fstderr, &
692 : num_hf_draw_points
693 :
694 : use corr_varnce_module, only: &
695 : hm_metadata_type
696 :
697 : use stats_type_utilities, only: &
698 : stat_begin_update, & ! Subroutines
699 : stat_end_update
700 :
701 : use stats_variables, only: &
702 : stats_metadata_type
703 :
704 : use error_code, only: &
705 : clubb_at_least_debug_level ! Procedure
706 :
707 : use stats_type, only: &
708 : stats ! Type
709 :
710 : implicit none
711 :
712 : !----------------------- Input Variables -----------------------
713 : type (grid), target, intent(in) :: &
714 : gr
715 :
716 : integer, intent(in) :: &
717 : hydromet_dim, &
718 : nz
719 :
720 : type (hm_metadata_type), intent(in) :: &
721 : hm_metadata
722 :
723 : logical, intent(in) :: &
724 : l_fill_holes_hm
725 :
726 : real( kind = core_rknd ), intent(in) :: &
727 : dt ! Timestep [s]
728 :
729 : real( kind = core_rknd ), dimension(nz), intent(in) :: &
730 : rho_ds_zm, & ! Dry, static density on momentum levels [kg/m^3]
731 : rho_ds_zt ! Dry, static density on thermo. levels [kg/m^3]
732 :
733 : real( kind = core_rknd ), dimension(nz), intent(in) :: &
734 : exner ! Exner function [-]
735 :
736 : type (stats_metadata_type), intent(in) :: &
737 : stats_metadata
738 :
739 : !----------------------- Input/Output Variables -----------------------
740 : type (stats), target, intent(inout) :: &
741 : stats_zt
742 :
743 : real( kind = core_rknd ), dimension(nz, hydromet_dim), intent(inout) :: &
744 : hydromet ! Mean of hydrometeor fields [units vary]
745 :
746 : real( kind = core_rknd ), dimension(nz), intent(inout) :: &
747 : rvm_mc, & ! Microphysics contributions to vapor water [kg/kg/s]
748 : thlm_mc ! Microphysics contributions to liquid potential temp [K/s]
749 :
750 : !----------------------- Local Variables -----------------------
751 : integer :: i, k ! Loop iterators
752 :
753 : real( kind = core_rknd ), dimension(nz, hydromet_dim) :: &
754 0 : hydromet_filled, & ! Frozen hydrometeor mixing ratios after hole filling
755 0 : hydromet_clipped ! Clipped mean of hydrometeor fields [units vary]
756 :
757 : character( len = 10 ) :: hydromet_name
758 :
759 : real( kind = core_rknd ) :: &
760 : max_velocity ! Maximum sedimentation velocity [m/s]
761 :
762 : integer :: ixrm_hf, ixrm_wvhf, ixrm_cl, &
763 : ixrm_bt, ixrm_mc
764 :
765 : logical :: l_hole_fill = .true.
766 :
767 : !----------------------- Begin Code -----------------------
768 :
769 : ! Start stats output for the _hf variables (changes in the hydromet array
770 : ! due to fill_holes_hydromet and fill_holes_vertical)
771 0 : if ( stats_metadata%l_stats_samp ) then
772 :
773 0 : do i = 1, hydromet_dim
774 :
775 : ! Set up the stats indices for hydrometeor at index i
776 : call setup_stats_indices( i, stats_metadata, hydromet_dim, & ! Intent(in)
777 : hm_metadata%hydromet_list, & ! Intent(in)
778 : ixrm_bt, ixrm_hf, ixrm_wvhf, & ! Intent(inout)
779 : ixrm_cl, ixrm_mc, & ! Intent(inout)
780 0 : max_velocity ) ! Intent(inout)
781 :
782 0 : call stat_begin_update( gr%nz, ixrm_hf, hydromet(:,i) / dt, & ! intent(in)
783 0 : stats_zt ) ! intent(inout)
784 :
785 : enddo ! i = 1, hydromet_dim
786 :
787 : endif ! stats_metadata%l_stats_samp
788 :
789 : ! If we're dealing with negative hydrometeors, we first try to fill the
790 : ! holes proportionally from other same-phase hydrometeors at each height
791 : ! level.
792 0 : if ( any( hydromet < zero_threshold ) .and. l_fill_holes_hm ) then
793 :
794 : call fill_holes_hydromet( nz, hydromet_dim, hydromet, & ! Intent(in)
795 : hm_metadata%l_frozen_hm, hm_metadata%l_mix_rat_hm, & ! Intent(in)
796 0 : hydromet_filled ) ! Intent(out)
797 :
798 0 : hydromet = hydromet_filled
799 :
800 : endif ! any( hydromet < zero ) .and. l_fill_holes_hm
801 :
802 0 : hydromet_filled = zero
803 :
804 0 : do i = 1, hydromet_dim
805 :
806 : ! Set up the stats indices for hydrometeor at index i
807 : call setup_stats_indices( i, stats_metadata, hydromet_dim, & ! Intent(in)
808 : hm_metadata%hydromet_list, & ! Intent(in)
809 : ixrm_bt, ixrm_hf, ixrm_wvhf, & ! Intent(inout)
810 : ixrm_cl, ixrm_mc, & ! Intent(inout)
811 0 : max_velocity ) ! Intent(inout)
812 :
813 0 : hydromet_name = hm_metadata%hydromet_list(i)
814 :
815 : ! Print warning message if any hydrometeor species has a value < 0.
816 0 : if ( clubb_at_least_debug_level( 1 ) ) then
817 0 : if ( any( hydromet(:,i) < zero_threshold ) ) then
818 :
819 0 : do k = 1, nz
820 0 : if ( hydromet(k,i) < zero_threshold ) then
821 0 : write(fstderr,*) trim( hydromet_name ) //" < ", &
822 0 : zero_threshold, &
823 0 : " in fill_holes_driver at k= ", k
824 : endif ! hydromet(k,i) < 0
825 : enddo ! k = 1, nz
826 :
827 : endif ! hydromet(:,i) < 0
828 : endif ! clubb_at_least_debug_level( 1 )
829 :
830 :
831 : ! Store the previous value of the hydrometeor for the effect of the
832 : ! hole-filling scheme.
833 : ! if ( stats_metadata%l_stats_samp ) then
834 : ! call stat_begin_update( ixrm_hf, hydromet(:,i) &
835 : ! / dt, stats_zt )
836 : ! endif
837 :
838 : ! If we're dealing with a mixing ratio and hole filling is enabled,
839 : ! then we apply the hole filling algorithm
840 0 : if ( any( hydromet(:,i) < zero_threshold ) ) then
841 :
842 0 : if ( hydromet_name(1:1) == "r" .and. l_hole_fill ) then
843 :
844 : !$acc data copyin( gr, gr%dzt, rho_ds_zt ) &
845 : !$acc copy( hydromet(:,i) )
846 :
847 : ! Apply the hole filling algorithm
848 : ! upper_hf_level = nz since we are filling the zt levels
849 : call fill_holes_vertical( gr%nz, 1, zero_threshold, gr%nz, & ! In
850 : gr%dzt, rho_ds_zt, & ! In
851 0 : hydromet(:,i) ) ! InOut
852 :
853 : !$acc end data
854 :
855 : endif ! Variable is a mixing ratio and l_hole_fill is true
856 :
857 : endif ! hydromet(:,i) < 0
858 :
859 : ! Enter the new value of the hydrometeor for the effect of the
860 : ! hole-filling scheme.
861 0 : if ( stats_metadata%l_stats_samp ) then
862 0 : call stat_end_update( gr%nz, ixrm_hf, hydromet(:,i) / dt, & ! intent(in)
863 0 : stats_zt ) ! intent(inout)
864 : endif
865 :
866 : ! Store the previous value of the hydrometeor for the effect of the water
867 : ! vapor hole-filling scheme.
868 0 : if ( stats_metadata%l_stats_samp ) then
869 0 : call stat_begin_update( gr%nz, ixrm_wvhf, hydromet(:,i) / dt, & ! intent(in)
870 0 : stats_zt ) ! intent(inout)
871 : endif
872 :
873 0 : if ( any( hydromet(:,i) < zero_threshold ) ) then
874 :
875 0 : if ( hydromet_name(1:1) == "r" .and. l_hole_fill ) then
876 :
877 : ! If the hole filling algorithm failed, then we attempt to fill
878 : ! the missing mass with water vapor mixing ratio.
879 : ! We noticed this is needed for ASEX A209, particularly if Latin
880 : ! hypercube sampling is enabled. -dschanen 11 Nov 2010
881 : call fill_holes_wv( nz, dt, exner, hydromet_name, & ! Intent(in)
882 0 : rvm_mc, thlm_mc, hydromet(:,i) ) ! Intent(out)
883 :
884 : endif ! Variable is a mixing ratio and l_hole_fill is true
885 :
886 : endif ! hydromet(:,i) < 0
887 :
888 : ! Enter the new value of the hydrometeor for the effect of the water vapor
889 : ! hole-filling scheme.
890 0 : if ( stats_metadata%l_stats_samp ) then
891 0 : call stat_end_update( gr%nz, ixrm_wvhf, hydromet(:,i) / dt, & ! intent(in)
892 0 : stats_zt ) ! intent(inout)
893 : endif
894 :
895 : ! Clipping for hydrometeor mixing ratios.
896 0 : if ( hm_metadata%l_mix_rat_hm(i) ) then
897 :
898 : ! Store the previous value of the hydrometeor for the effect of
899 : ! clipping.
900 0 : if ( stats_metadata%l_stats_samp ) then
901 0 : call stat_begin_update( gr%nz, ixrm_cl, hydromet(:,i) / dt, & ! intent(in)
902 0 : stats_zt ) ! intent(inout)
903 : endif
904 :
905 0 : if ( any( hydromet(:,i) < zero_threshold ) ) then
906 :
907 : ! Clip any remaining negative values of precipitating hydrometeor
908 : ! mixing ratios to 0.
909 0 : where ( hydromet(:,i) < zero_threshold )
910 : hydromet(:,i) = zero_threshold
911 : end where
912 :
913 : endif ! hydromet(:,i) < 0
914 :
915 : ! Eliminate very small values of mean precipitating hydrometeor mixing
916 : ! ratios by setting them to 0.
917 0 : do k = 2, gr%nz, 1
918 :
919 0 : if ( hydromet(k,i) <= hm_metadata%hydromet_tol(i) ) then
920 :
921 : rvm_mc(k) &
922 : = rvm_mc(k) &
923 0 : + ( hydromet(k,i) / dt )
924 :
925 0 : if ( .not. hm_metadata%l_frozen_hm(i) ) then
926 :
927 : ! Rain water mixing ratio
928 :
929 : thlm_mc(k) &
930 : = thlm_mc(k) &
931 : - ( Lv / ( Cp * exner(k) ) ) &
932 0 : * ( hydromet(k,i) / dt )
933 :
934 : else ! Frozen hydrometeor mixing ratio
935 :
936 : thlm_mc(k) &
937 : = thlm_mc(k) &
938 : - ( Ls / ( Cp * exner(k) ) ) &
939 0 : * ( hydromet(k,i) / dt )
940 :
941 : endif ! l_frozen_hm(i)
942 :
943 0 : hydromet(k,i) = zero
944 :
945 : endif ! hydromet(k,i) <= hydromet_tol(i)
946 :
947 : enddo ! k = 2, gr%nz, 1
948 :
949 :
950 : ! Enter the new value of the hydrometeor for the effect of clipping.
951 0 : if ( stats_metadata%l_stats_samp ) then
952 : call stat_end_update( gr%nz, ixrm_cl, hydromet(:,i) / dt, & ! intent(in)
953 0 : stats_zt ) ! intent(inout)
954 : endif
955 :
956 : endif ! l_mix_rat_hm(i)
957 :
958 : enddo ! i = 1, hydromet_dim, 1
959 :
960 : ! Calculate clipping for hydrometeor concentrations.
961 : call clip_hydromet_conc_mvr( gr%nz, hydromet_dim, hm_metadata, & ! Intent(in)
962 : hydromet, & ! Intent(in)
963 0 : hydromet_clipped ) ! Intent(out)
964 :
965 : ! Clip hydrometeor concentrations and output stats.
966 0 : do i = 1, hydromet_dim
967 :
968 0 : if ( .not. hm_metadata%l_mix_rat_hm(i) ) then
969 :
970 0 : if ( stats_metadata%l_stats_samp ) then
971 :
972 : ! Set up the stats indices for hydrometeor at index i
973 : call setup_stats_indices( i, stats_metadata, hydromet_dim, & ! Intent(in)
974 : hm_metadata%hydromet_list, & ! Intent(in)
975 : ixrm_bt, ixrm_hf, ixrm_wvhf, & ! In/Out
976 : ixrm_cl, ixrm_mc, & ! In/Out
977 0 : max_velocity ) ! In/Out
978 :
979 : ! Store the previous value of the hydrometeor for the effect of
980 : ! clipping.
981 0 : call stat_begin_update( gr%nz, ixrm_cl, hydromet(:,i) / dt, & ! intent(in)
982 0 : stats_zt ) ! intent(inout)
983 :
984 : endif ! stats_metadata%l_stats_samp
985 :
986 : ! Apply clipping of hydrometeor concentrations.
987 0 : hydromet(:,i) = hydromet_clipped(:,i)
988 :
989 : ! Enter the new value of the hydrometeor for the effect of clipping.
990 0 : if ( stats_metadata%l_stats_samp ) then
991 : call stat_end_update( gr%nz, ixrm_cl, hydromet(:,i) / dt, & ! intent(in)
992 0 : stats_zt ) ! intent(inout)
993 : endif
994 :
995 : endif ! .not. l_mix_rat_hm(i)
996 :
997 : enddo ! i = 1, hydromet_dim, 1
998 :
999 :
1000 0 : return
1001 :
1002 : end subroutine fill_holes_driver
1003 :
1004 : !=============================================================================
1005 0 : subroutine clip_hydromet_conc_mvr( nz, hydromet_dim, hm_metadata, & ! Intent(in)
1006 0 : hydromet, & ! Intent(in)
1007 0 : hydromet_clipped ) ! Intent(out)
1008 :
1009 : ! Description:
1010 : ! Increases the value of a hydrometeor concentration when it is too small,
1011 : ! according to the hydrometeor mixing ratio (which remains unchanged) and
1012 : ! the maximum acceptable drop or particle mean volume radius for that
1013 : ! hydrometeor species.
1014 :
1015 : ! References:
1016 : !-----------------------------------------------------------------------
1017 :
1018 : use grid_class, only: &
1019 : grid ! Type
1020 :
1021 : use corr_varnce_module, only: &
1022 : hm_metadata_type
1023 :
1024 : use constants_clubb, only: &
1025 : pi, & ! Variable(s)
1026 : four_thirds, &
1027 : one, &
1028 : zero, &
1029 : rho_lw, &
1030 : rho_ice
1031 :
1032 : use index_mapping, only: &
1033 : Nx2rx_hm_idx, & ! Procedure(s)
1034 : mvr_hm_max
1035 :
1036 : use clubb_precision, only: &
1037 : core_rknd ! Variable(s)
1038 :
1039 : implicit none
1040 :
1041 : ! Input Variables
1042 : integer, intent(in) :: &
1043 : nz
1044 :
1045 : integer, intent(in) :: &
1046 : hydromet_dim ! Number of hydrometeor fields
1047 :
1048 : type (hm_metadata_type), intent(in) :: &
1049 : hm_metadata
1050 :
1051 : real( kind = core_rknd ), dimension(nz,hydromet_dim), intent(in) :: &
1052 : hydromet ! Mean of hydrometeor fields [units vary]
1053 :
1054 : ! Output Variable
1055 : real( kind = core_rknd ), dimension(nz,hydromet_dim), intent(out) :: &
1056 : hydromet_clipped ! Clipped mean of hydrometeor fields [units vary]
1057 :
1058 : ! Local Variables
1059 : real( kind = core_rknd ) :: &
1060 : Nxm_min_coef ! Coefficient for min mean value of a concentration [1/kg]
1061 :
1062 : integer :: &
1063 : k, & ! Vertical grid index
1064 : idx ! Hydrometeor species index
1065 :
1066 :
1067 : ! Clipping for hydrometeor concentrations.
1068 0 : do idx = 1, hydromet_dim
1069 :
1070 0 : if ( .not. hm_metadata%l_mix_rat_hm(idx) ) then
1071 :
1072 : ! The field is a hydrometeor concentration.
1073 :
1074 0 : if ( .not. hm_metadata%l_frozen_hm(idx) ) then
1075 :
1076 : ! Clipping for mean rain drop concentration, <Nr>.
1077 : !
1078 : ! When mean rain water mixing ratio, <rr>, is found at a grid
1079 : ! level, mean rain drop concentration must be at least a minimum
1080 : ! value so that average rain drop mean volume radius stays within
1081 : ! an upper bound. Otherwise, mean rain drop concentration is 0.
1082 : ! This can also be applied to values at a sample point, rather than
1083 : ! a grid-box mean.
1084 :
1085 : ! The minimum mean rain drop concentration is given by:
1086 : !
1087 : ! <Nr> = <rr> / ( (4/3) * pi * rho_lw * mvr_rain_max^3 );
1088 : !
1089 : ! where mvr_rain_max is the maximum acceptable average rain drop
1090 : ! mean volume radius.
1091 :
1092 : Nxm_min_coef &
1093 0 : = one / ( four_thirds * pi * rho_lw * mvr_hm_max(idx,hm_metadata)**3 )
1094 :
1095 : else ! l_frozen_hm(idx)
1096 :
1097 : ! Clipping for mean frozen hydrometeor concentration, <Nx>, where x
1098 : ! stands for any frozen hydrometeor species.
1099 : !
1100 : ! When mean frozen hydrometeor mixing ratio, <rx>, is found at a
1101 : ! grid level, mean frozen hydrometeor concentration must be at
1102 : ! least a minimum value so that average frozen hydrometeor mean
1103 : ! volume radius stays within an upper bound. Otherwise, mean
1104 : ! frozen hydrometeor concentration is 0. This can also be applied
1105 : ! to values at a sample point, rather than a grid-box mean.
1106 :
1107 : ! The minimum mean frozen hydrometeor concentration is given by:
1108 : !
1109 : ! <Nx> = <rx> / ( (4/3) * pi * rho_ice * mvr_x_max^3 );
1110 : !
1111 : ! where mvr_x_max is the maximum acceptable average frozen
1112 : ! hydrometeor mean volume radius for frozen hydrometeor species, x.
1113 :
1114 : Nxm_min_coef &
1115 0 : = one / ( four_thirds * pi * rho_ice * mvr_hm_max(idx,hm_metadata)**3 )
1116 :
1117 : endif ! .not. l_frozen_hm(idx)
1118 :
1119 : ! Loop over vertical levels and increase hydrometeor concentrations
1120 : ! when necessary.
1121 0 : do k = 2, nz, 1
1122 :
1123 0 : if ( hydromet(k,Nx2rx_hm_idx(idx,hm_metadata)) > zero ) then
1124 :
1125 : ! Hydrometeor mixing ratio, <rx>, is found at the grid level.
1126 0 : hydromet_clipped(k,idx) &
1127 : = max( hydromet(k,idx), &
1128 0 : Nxm_min_coef * hydromet(k,Nx2rx_hm_idx(idx,hm_metadata)) )
1129 :
1130 : else ! <rx> = 0
1131 :
1132 0 : hydromet_clipped(k,idx) = zero
1133 :
1134 : endif ! hydromet(k,Nx2rx_hm_idx(idx)) > 0
1135 :
1136 : enddo ! k = 2, nz, 1
1137 :
1138 : ! The lowest thermodynamic level is below the model's lower boundary.
1139 0 : hydromet_clipped(1,idx) = hydromet(1,idx)
1140 :
1141 : else ! l_mix_rat_hm(idx)
1142 :
1143 : ! The field is a hydrometeor mixing ratio.
1144 0 : hydromet_clipped(:,idx) = hydromet(:,idx)
1145 :
1146 : endif ! .not. l_mix_rat_hm(idx)
1147 :
1148 : enddo ! idx = 1, hydromet_dim, 1
1149 :
1150 :
1151 0 : return
1152 :
1153 : end subroutine clip_hydromet_conc_mvr
1154 :
1155 : !-----------------------------------------------------------------------
1156 0 : subroutine setup_stats_indices( ihm, stats_metadata, hydromet_dim, & ! Intent(in)
1157 0 : hydromet_list, & ! Intent(in)
1158 : ixrm_bt, ixrm_hf, ixrm_wvhf, & ! Intent(inout)
1159 : ixrm_cl, ixrm_mc, & ! Intent(inout)
1160 : max_velocity ) ! Intent(inout)
1161 :
1162 : ! Description:
1163 : !
1164 : ! Determines the stats output indices depending on the hydrometeor.
1165 :
1166 : ! Attention: hydromet_list needs to be set up before this routine is called.
1167 : !
1168 : ! Bogus example
1169 : ! References:
1170 : !
1171 : ! None
1172 : !-----------------------------------------------------------------------
1173 :
1174 : use clubb_precision, only: &
1175 : core_rknd
1176 :
1177 : use constants_clubb, only: &
1178 : zero
1179 :
1180 : use stats_variables, only: &
1181 : stats_metadata_type
1182 :
1183 : implicit none
1184 :
1185 : ! Input Variables
1186 : integer, intent(in) :: ihm, hydromet_dim
1187 :
1188 : type (stats_metadata_type), intent(in) :: &
1189 : stats_metadata
1190 :
1191 : character(len=10), dimension(hydromet_dim), intent(in) :: &
1192 : hydromet_list
1193 :
1194 : ! Input/Output Variables
1195 : real( kind = core_rknd ), intent(inout) :: &
1196 : max_velocity ! Maximum sedimentation velocity [m/s]
1197 :
1198 : integer, intent(inout) :: ixrm_hf, ixrm_wvhf, ixrm_cl, &
1199 : ixrm_bt, ixrm_mc
1200 :
1201 : !-----------------------------------------------------------------------
1202 :
1203 : !----- Begin Code -----
1204 :
1205 : ! Initializing max_velocity in order to avoid a compiler warning.
1206 : ! Regardless of the case, it will be reset in the 'select case'
1207 : ! statement immediately below.
1208 0 : max_velocity = zero
1209 :
1210 0 : select case ( trim( hydromet_list(ihm) ) )
1211 : case ( "rrm" )
1212 0 : ixrm_bt = stats_metadata%irrm_bt
1213 0 : ixrm_hf = stats_metadata%irrm_hf
1214 0 : ixrm_wvhf = stats_metadata%irrm_wvhf
1215 0 : ixrm_cl = stats_metadata%irrm_cl
1216 0 : ixrm_mc = stats_metadata%irrm_mc
1217 :
1218 0 : max_velocity = -9.1_core_rknd ! m/s
1219 :
1220 : case ( "rim" )
1221 0 : ixrm_bt = stats_metadata%irim_bt
1222 0 : ixrm_hf = stats_metadata%irim_hf
1223 0 : ixrm_wvhf = stats_metadata%irim_wvhf
1224 0 : ixrm_cl = stats_metadata%irim_cl
1225 0 : ixrm_mc = stats_metadata%irim_mc
1226 :
1227 0 : max_velocity = -1.2_core_rknd ! m/s
1228 :
1229 : case ( "rsm" )
1230 0 : ixrm_bt = stats_metadata%irsm_bt
1231 0 : ixrm_hf = stats_metadata%irsm_hf
1232 0 : ixrm_wvhf = stats_metadata%irsm_wvhf
1233 0 : ixrm_cl = stats_metadata%irsm_cl
1234 0 : ixrm_mc = stats_metadata%irsm_mc
1235 :
1236 : ! Morrison limit
1237 : ! max_velocity = -1.2_core_rknd ! m/s
1238 : ! Made up limit. The literature suggests that it is quite possible
1239 : ! that snow flake might achieve a terminal velocity of 2 m/s, and this
1240 : ! happens in the COAMPS microphysics -dschanen 29 Sept 2009
1241 0 : max_velocity = -2.0_core_rknd ! m/s
1242 :
1243 : case ( "rgm" )
1244 0 : ixrm_bt = stats_metadata%irgm_bt
1245 0 : ixrm_hf = stats_metadata%irgm_hf
1246 0 : ixrm_wvhf = stats_metadata%irgm_wvhf
1247 0 : ixrm_cl = stats_metadata%irgm_cl
1248 0 : ixrm_mc = stats_metadata%irgm_mc
1249 :
1250 0 : max_velocity = -20._core_rknd ! m/s
1251 :
1252 : case ( "Nrm" )
1253 0 : ixrm_bt = stats_metadata%iNrm_bt
1254 0 : ixrm_hf = 0
1255 0 : ixrm_wvhf = 0
1256 0 : ixrm_cl = stats_metadata%iNrm_cl
1257 0 : ixrm_mc = stats_metadata%iNrm_mc
1258 :
1259 0 : max_velocity = -9.1_core_rknd ! m/s
1260 :
1261 : case ( "Nim" )
1262 0 : ixrm_bt = stats_metadata%iNim_bt
1263 0 : ixrm_hf = 0
1264 0 : ixrm_wvhf = 0
1265 0 : ixrm_cl = stats_metadata%iNim_cl
1266 0 : ixrm_mc = stats_metadata%iNim_mc
1267 :
1268 0 : max_velocity = -1.2_core_rknd ! m/s
1269 :
1270 : case ( "Nsm" )
1271 0 : ixrm_bt = stats_metadata%iNsm_bt
1272 0 : ixrm_hf = 0
1273 0 : ixrm_wvhf = 0
1274 0 : ixrm_cl = stats_metadata%iNsm_cl
1275 0 : ixrm_mc = stats_metadata%iNsm_mc
1276 :
1277 : ! Morrison limit
1278 : ! max_velocity = -1.2_core_rknd ! m/s
1279 : ! Made up limit. The literature suggests that it is quite possible
1280 : ! that snow flake might achieve a terminal velocity of 2 m/s, and this
1281 : ! happens in the COAMPS microphysics -dschanen 29 Sept 2009
1282 0 : max_velocity = -2.0_core_rknd ! m/s
1283 :
1284 : case ( "Ngm" )
1285 0 : ixrm_bt = stats_metadata%iNgm_bt
1286 0 : ixrm_hf = 0
1287 0 : ixrm_wvhf = 0
1288 0 : ixrm_cl = stats_metadata%iNgm_cl
1289 0 : ixrm_mc = stats_metadata%iNgm_mc
1290 :
1291 0 : max_velocity = -20._core_rknd ! m/s
1292 :
1293 : case ( "Ncm" )
1294 0 : ixrm_bt = stats_metadata%iNcm_bt
1295 0 : ixrm_hf = 0
1296 0 : ixrm_wvhf = 0
1297 0 : ixrm_cl = stats_metadata%iNcm_cl
1298 0 : ixrm_mc = stats_metadata%iNcm_mc
1299 :
1300 : ! Use the rain water limit, since Morrison has no explicit limit on
1301 : ! cloud water. Presumably these numbers are never large.
1302 : ! -dschanen 28 Sept 2009
1303 0 : max_velocity = -9.1_core_rknd ! m/s
1304 :
1305 : case default
1306 0 : ixrm_bt = 0
1307 0 : ixrm_hf = 0
1308 0 : ixrm_wvhf = 0
1309 0 : ixrm_cl = 0
1310 0 : ixrm_mc = 0
1311 :
1312 0 : max_velocity = -9.1_core_rknd ! m/s
1313 :
1314 : end select
1315 :
1316 :
1317 0 : return
1318 :
1319 : end subroutine setup_stats_indices
1320 : !-----------------------------------------------------------------------
1321 :
1322 : end module fill_holes
|