Line data Source code
1 : !-------------------------------------------------------------------------
2 : ! $Id$
3 : !===============================================================================
4 : module precipitation_fraction
5 :
6 : ! Description:
7 : ! Sets overall precipitation fraction as well as the precipitation fraction
8 : ! in each PDF component.
9 :
10 : implicit none
11 :
12 : private
13 :
14 : public :: precip_fraction
15 :
16 : private :: component_precip_frac_weighted, &
17 : component_precip_frac_specify, &
18 : precip_frac_assert_check
19 :
20 : integer, parameter, public :: &
21 : precip_frac_calc_type = 2 ! Option used to calculate component precip_frac
22 :
23 : contains
24 :
25 : !=============================================================================
26 0 : subroutine precip_fraction( nz, ngrdcol, hydromet_dim, &
27 0 : hydromet, cloud_frac, cloud_frac_1, &
28 0 : l_mix_rat_hm, l_frozen_hm, hydromet_tol, &
29 0 : cloud_frac_2, ice_supersat_frac, &
30 0 : ice_supersat_frac_1, ice_supersat_frac_2, &
31 0 : mixt_frac, clubb_params, &
32 : stats_metadata, &
33 0 : stats_sfc, &
34 0 : precip_frac, &
35 0 : precip_frac_1, &
36 0 : precip_frac_2, &
37 0 : precip_frac_tol )
38 :
39 : ! Description:
40 : ! Determines (overall) precipitation fraction over the horizontal domain, as
41 : ! well as the precipitation fraction within each PDF component, at every
42 : ! vertical grid level.
43 :
44 : ! References:
45 : !-----------------------------------------------------------------------
46 :
47 : use constants_clubb, only: &
48 : one, & ! Constant(s)
49 : zero, &
50 : cloud_frac_min, &
51 : fstderr
52 :
53 : use parameter_indices, only: &
54 : nparams, & ! Variable(s)
55 : iupsilon_precip_frac_rat
56 :
57 : use stats_type_utilities, only: &
58 : stat_update_var_pt ! Procedure(s)
59 :
60 : use clubb_precision, only: &
61 : core_rknd ! Variable(s)
62 :
63 : use error_code, only: &
64 : clubb_at_least_debug_level, & ! Procedure
65 : err_code, & ! Error Indicator
66 : clubb_fatal_error ! Constant
67 :
68 : use stats_type, only: &
69 : stats ! Type
70 :
71 : use stats_variables, only: &
72 : stats_metadata_type
73 :
74 : implicit none
75 :
76 : !------------------------- Input Variables -------------------------
77 : integer, intent(in) :: &
78 : nz, & ! Number of model vertical grid levels
79 : ngrdcol, & ! Number of grid columns
80 : hydromet_dim ! Number of hydrometeor species
81 :
82 : real( kind = core_rknd ), dimension(ngrdcol,nz,hydromet_dim), intent(in) :: &
83 : hydromet ! Mean of hydrometeor, hm (overall) [units vary]
84 :
85 : logical, dimension(hydromet_dim), intent(in) :: &
86 : l_frozen_hm, & ! if true, then the hydrometeor is frozen; otherwise liquid
87 : l_mix_rat_hm ! if true, then the quantity is a hydrometeor mixing ratio
88 :
89 : real( kind = core_rknd ), dimension(hydromet_dim), intent(in) :: &
90 : hydromet_tol ! Tolerance values for all hydrometeors [units vary]
91 :
92 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
93 : cloud_frac, & ! Cloud fraction (overall) [-]
94 : cloud_frac_1, & ! Cloud fraction (1st PDF component) [-]
95 : cloud_frac_2, & ! Cloud fraction (2nd PDF component) [-]
96 : ice_supersat_frac, & ! Ice supersaturation fraction (overall) [-]
97 : ice_supersat_frac_1, & ! Ice supersaturation fraction (1st PDF comp.) [-]
98 : ice_supersat_frac_2, & ! Ice supersaturation fraction (2nd PDF comp.) [-]
99 : mixt_frac ! Mixture fraction [-]
100 :
101 : real( kind = core_rknd ), dimension(nparams), intent(in) :: &
102 : clubb_params ! Array of CLUBB's tunable parameters [units vary]
103 :
104 : type (stats_metadata_type), intent(in) :: &
105 : stats_metadata
106 :
107 : !------------------------- Inout Variables -------------------------
108 : type (stats), target, dimension(ngrdcol), intent(inout) :: &
109 : stats_sfc
110 :
111 : !------------------------- Output Variables -------------------------
112 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
113 : precip_frac, & ! Precipitation fraction (overall) [-]
114 : precip_frac_1, & ! Precipitation fraction (1st PDF component) [-]
115 : precip_frac_2 ! Precipitation fraction (2nd PDF component) [-]
116 :
117 : real( kind = core_rknd ), dimension(ngrdcol), intent(out) :: &
118 : precip_frac_tol ! Minimum precip. frac. when hydromet. are present [-]
119 :
120 : !------------------------- Local Variables -------------------------
121 :
122 : ! "Maximum allowable" hydrometeor mixing ratio in-precip component mean.
123 : real( kind = core_rknd ), parameter :: &
124 : max_hm_ip_comp_mean = 0.0025_core_rknd ! [kg/kg]
125 :
126 : real( kind = core_rknd ), parameter :: &
127 : precip_frac_tol_coef = 0.1_core_rknd ! Coefficient for precip_frac_tol
128 :
129 : integer :: &
130 : j, k, ivar ! Loop indices
131 :
132 :
133 : ! Initialize the precipitation fraction variables (precip_frac,
134 : ! precip_frac_1, and precip_frac_2) to 0.
135 0 : precip_frac(:,:) = zero
136 0 : precip_frac_1(:,:) = zero
137 0 : precip_frac_2(:,:) = zero
138 :
139 : ! Set the minimum allowable precipitation fraction when hydrometeors are
140 : ! found at a grid level.
141 0 : if ( any( l_frozen_hm ) ) then
142 : ! Ice microphysics included.
143 0 : do j = 1, ngrdcol
144 0 : precip_frac_tol(j) &
145 : = max( precip_frac_tol_coef &
146 : * max( maxval( cloud_frac(j,:) ), maxval( ice_supersat_frac(j,:) ) ), &
147 0 : cloud_frac_min )
148 : end do
149 : else
150 : ! Warm microphysics.
151 0 : do j = 1, ngrdcol
152 0 : precip_frac_tol(j) = max( precip_frac_tol_coef * maxval( cloud_frac(j,:) ), &
153 0 : cloud_frac_min )
154 : end do
155 : endif
156 :
157 : !!! Find overall precipitation fraction.
158 : ! The precipitation fraction is the greatest cloud fraction at or above a
159 : ! vertical level.
160 0 : if ( any( l_frozen_hm ) ) then
161 :
162 : ! Ice microphysics included.
163 0 : precip_frac(:,nz) = max( cloud_frac(:,nz), ice_supersat_frac(:,nz) )
164 :
165 0 : do k = nz-1, 1, -1
166 0 : precip_frac(:,k) = max( precip_frac(:,k+1), cloud_frac(:,k), &
167 0 : ice_supersat_frac(:,k) )
168 : end do
169 :
170 : else
171 :
172 : ! Warm microphysics.
173 0 : precip_frac(:,nz) = cloud_frac(:,nz)
174 :
175 0 : do k = nz-1, 1, -1
176 0 : precip_frac(:,k) = max( precip_frac(:,k+1), cloud_frac(:,k) )
177 : end do
178 :
179 : end if
180 :
181 : !!! Special checks for overall precipitation fraction
182 0 : do k = 1, nz, 1
183 0 : do j = 1, ngrdcol
184 :
185 0 : if ( any( hydromet(j,k,:) >= hydromet_tol(:) ) &
186 0 : .and. precip_frac(j,k) < precip_frac_tol(j) ) then
187 :
188 : ! In a scenario where we find any hydrometeor at this grid level, but
189 : ! no cloud at or above this grid level, set precipitation fraction to
190 : ! a minimum threshold value.
191 0 : precip_frac(j,k) = precip_frac_tol(j)
192 :
193 0 : elseif ( all( hydromet(j,k,:) < hydromet_tol(:) ) ) then
194 :
195 : ! The means (overall) of every precipitating hydrometeor are all less
196 : ! than their respective tolerance amounts. They are all considered to
197 : ! have values of 0. There are not any hydrometeor species found at
198 : ! this grid level. There is also no cloud at or above this grid
199 : ! level, so set precipitation fraction to 0.
200 0 : precip_frac(j,k) = zero
201 :
202 : endif
203 :
204 : end do
205 : enddo ! Special checks for overall precipitation fraction loop: k = 1, nz, 1
206 :
207 :
208 : !!! Find precipitation fraction within each PDF component.
209 : !
210 : ! The overall precipitation fraction, f_p, is given by the equation:
211 : !
212 : ! f_p = a * f_p(1) + ( 1 - a ) * f_p(2);
213 : !
214 : ! where "a" is the mixture fraction (weight of PDF component 1), f_p(1) is
215 : ! the precipitation fraction within PDF component 1, and f_p(2) is the
216 : ! precipitation fraction within PDF component 2. Overall precipitation
217 : ! fraction is found according the method above, and mixture fraction is
218 : ! already determined, leaving f_p(1) and f_p(2) to be solved for. The
219 : ! values for f_p(1) and f_p(2) must satisfy the above equation.
220 0 : if ( precip_frac_calc_type == 1 ) then
221 :
222 : ! Calculatate precip_frac_1 and precip_frac_2 based on the greatest
223 : ! weighted cloud_frac_1 at or above a grid level.
224 : call component_precip_frac_weighted( nz, ngrdcol, hydromet_dim, & ! intent(in)
225 : l_frozen_hm, hydromet_tol, & ! intent(in)
226 : hydromet(:,:,:), precip_frac(:,:), & ! intent(in)
227 : cloud_frac_1(:,:), cloud_frac_2(:,:), & ! intent(in)
228 : ice_supersat_frac_1(:,:), & ! intent(in)
229 : ice_supersat_frac_2(:,:), mixt_frac(:,:), & !intent(in)
230 : precip_frac_tol(:), & ! intent(in)
231 : precip_frac_1(:,:), precip_frac_2(:,:) ) ! intent(out)
232 :
233 : elseif ( precip_frac_calc_type == 2 ) then
234 :
235 : ! Specified method.
236 : call component_precip_frac_specify( nz, ngrdcol, hydromet_dim, hydromet_tol, & ! intent(in)
237 : clubb_params(iupsilon_precip_frac_rat), & ! intent(in)
238 : hydromet(:,:,:), precip_frac(:,:), & ! intent(in)
239 : mixt_frac(:,:), precip_frac_tol(:), & ! intent(in)
240 : precip_frac_1(:,:), precip_frac_2(:,:) ) ! intent(out)
241 :
242 : else ! Invalid option selected.
243 :
244 : write(fstderr,*) "Invalid option to calculate precip_frac_1 " &
245 : // "and precip_frac_2."
246 : err_code = clubb_fatal_error
247 : return
248 :
249 : endif ! precip_frac_calc_type
250 :
251 :
252 : ! Increase Precipiation Fraction under special conditions.
253 : !
254 : ! There are scenarios that sometimes occur that require precipitation
255 : ! fraction to be boosted. Precipitation fraction is calculated from cloud
256 : ! fraction and ice supersaturation fraction. For numerical reasons, CLUBB's
257 : ! PDF may become entirely subsaturated with respect to liquid and ice,
258 : ! resulting in both a cloud fraction of 0 and an ice supersaturation
259 : ! fraction of 0. When this happens, precipitation fraction drops to 0 when
260 : ! there aren't any hydrometeors present at that grid level, or to
261 : ! precip_frac_tol when there is at least one hydrometeor present at that
262 : ! grid level. However, sometimes there are large values of hydrometeors
263 : ! found at that grid level. When this occurs, the PDF component in-precip
264 : ! mean of a hydrometeor can become ridiculously large. This is because the
265 : ! ith PDF component in-precip mean of a hydrometeor, mu_hm_i, is given by
266 : ! the equation:
267 : !
268 : ! mu_hm_i = hm_i / precip_frac_i;
269 : !
270 : ! where hm_i is the overall ith PDF component mean of the hydrometeor, and
271 : ! precip_frac_i is the ith PDF component precipitation fraction. When
272 : ! precip_frac_i has a value of precip_frac_tol and hm_i is large, mu_hm_i
273 : ! can be huge. This can cause enormous microphysical process rates and
274 : ! result in numerical instability. It is also very inaccurate.
275 : !
276 : ! In order to limit this problem, the ith PDF component precipitation
277 : ! fraction is increased in order to decrease mu_hm_i. First, an "upper
278 : ! limit" is set for mu_hm_i when the hydrometeor is a mixing ratio. This is
279 : ! called max_hm_ip_comp_mean. At every vertical level and for every
280 : ! hydrometeor mixing ratio, a check is made to try to prevent mu_hm_i from
281 : ! exceeding the "upper limit". The check is:
282 : !
283 : ! hm_i / precip_frac_i ( which = mu_hm_i ) > max_hm_ip_comp_mean,
284 : !
285 : ! which can be rewritten:
286 : !
287 : ! hm_i > precip_frac_i * max_hm_ip_comp_mean.
288 : !
289 : ! Since hm_i has not been calculated yet, the assumption for this check is
290 : ! that all of the hydrometeor is found in one PDF component, which is the
291 : ! worst-case scenario in violating this limit. The check becomes:
292 : !
293 : ! <hm> / ( mixt_frac * precip_frac_1 ) > max_hm_ip_comp_mean;
294 : ! in PDF comp. 1; and
295 : ! <hm> / ( ( 1 - mixt_frac ) * precip_frac_2 ) > max_hm_ip_comp_mean;
296 : ! in PDF comp. 2.
297 : !
298 : ! These limits can be rewritten as:
299 : !
300 : ! <hm> > mixt_frac * precip_frac_1 * max_hm_ip_comp_mean;
301 : ! in PDF comp. 1; and
302 : ! <hm> > ( 1 - mixt_frac ) * precip_frac_2 * max_hm_ip_comp_mean;
303 : ! in PDF comp. 2.
304 : !
305 : ! When component precipitation fraction is found to be in excess of the
306 : ! limit, precip_frac_i is increased to:
307 : !
308 : ! <hm> / ( mixt_frac * max_hm_ip_comp_mean );
309 : ! when the limit is exceeded in PDF comp. 1; and
310 : ! <hm> / ( ( 1 - mixt_frac ) * max_hm_ip_comp_mean );
311 : ! when the limit is exceeded in PDF comp. 2.
312 : !
313 : ! Of course, precip_frac_i is not allowed to exceed 1, so when
314 : ! <hm> / mixt_frac (or <hm> / ( 1 - mixt_frac )) is already greater than
315 : ! max_hm_ip_comp_mean, mu_hm_i will also have to be greater than
316 : ! max_hm_ip_comp_mean. However, the value of mu_hm_i is still reduced when
317 : ! compared to what it would have been using precip_frac_tol. In the event
318 : ! that multiple hydrometeor mixing ratios violate the check, the code is set
319 : ! up so that precip_frac_i is increased based on the highest hm_i.
320 :
321 :
322 0 : do ivar = 1, hydromet_dim
323 0 : if ( l_mix_rat_hm(ivar) ) then
324 0 : do k = 1, nz
325 0 : do j = 1, ngrdcol
326 :
327 : ! The hydrometeor is a mixing ratio.
328 :
329 0 : if ( hydromet(j,k,ivar) >= hydromet_tol(ivar) .and. &
330 : hydromet(j,k,ivar) > mixt_frac(j,k) * precip_frac_1(j,k) &
331 : * max_hm_ip_comp_mean ) then
332 :
333 : ! Increase precipitation fraction in the 1st PDF component.
334 : precip_frac_1(j,k) &
335 : = min( hydromet(j,k,ivar) &
336 0 : / ( mixt_frac(j,k) * max_hm_ip_comp_mean ), one )
337 :
338 : ! The value of precip_frac_1 must be at least precip_frac_tol
339 : ! when precipitation is found in the 1st PDF component.
340 0 : precip_frac_1(j,k) = max( precip_frac_1(j,k), precip_frac_tol(j) )
341 :
342 : endif ! <hm>/(mixt_frac*precip_frac_1) > max_hm_ip_comp_mean
343 :
344 0 : if ( hydromet(j,k,ivar) >= hydromet_tol(ivar) .and. &
345 : hydromet(j,k,ivar) > ( one - mixt_frac(j,k) ) * precip_frac_2(j,k) &
346 0 : * max_hm_ip_comp_mean ) then
347 :
348 : ! Increase precipitation fraction in the 2nd PDF component.
349 : precip_frac_2(j,k) &
350 : = min( hydromet(j,k,ivar) &
351 0 : / ( ( one - mixt_frac(j,k) ) * max_hm_ip_comp_mean ), one )
352 :
353 : ! The value of precip_frac_2 must be at least precip_frac_tol
354 : ! when precipitation is found in the 2nd PDF component.
355 0 : precip_frac_2(j,k) = max( precip_frac_2(j,k), precip_frac_tol(j) )
356 :
357 : endif ! <hm>/((1-mixt_frac)*precip_frac_2) > max_hm_ip_comp_mean
358 :
359 : end do
360 : end do ! k = 1, nz
361 : end if ! l_mix_rat_hm(ivar)
362 : end do ! ivar = 1, hydromet_dim
363 :
364 : ! Recalculate overall precipitation fraction for consistency.
365 : precip_frac(:,:) = mixt_frac(:,:) * precip_frac_1(:,:) &
366 0 : + ( one - mixt_frac(:,:) ) * precip_frac_2(:,:)
367 :
368 : ! Double check that precip_frac_tol <= precip_frac <= 1 when hydrometeors
369 : ! are found at a grid level.
370 : ! PLEASE DO NOT ALTER precip_frac, precip_frac_1, or precip_frac_2 anymore
371 : ! after this point in the code.
372 0 : do k = 1, nz, 1
373 0 : do j = 1, ngrdcol
374 0 : if ( any( hydromet(j,k,:) >= hydromet_tol(:) ) ) then
375 0 : precip_frac(j,k) = min( max( precip_frac(j,k), precip_frac_tol(j) ), one )
376 : end if ! any( hydromet(k,:) >= hydromet_tol(:) )
377 : end do
378 : end do ! k = 1, nz, 1
379 :
380 :
381 : ! Statistics
382 0 : if ( stats_metadata%l_stats_samp ) then
383 0 : if ( stats_metadata%iprecip_frac_tol > 0 ) then
384 0 : do j = 1, ngrdcol
385 0 : call stat_update_var_pt( stats_metadata%iprecip_frac_tol, 1, precip_frac_tol(j), & ! intent(in)
386 0 : stats_sfc(j) ) ! intent(inout)
387 : end do
388 : end if ! stats_metadata%iprecip_frac_tol
389 : end if ! stats_metadata%l_stats_samp
390 :
391 :
392 : ! Assertion check for precip_frac, precip_frac_1, and precip_frac_2.
393 0 : if ( clubb_at_least_debug_level( 2 ) ) then
394 0 : do j = 1, ngrdcol
395 : call precip_frac_assert_check( nz, hydromet_dim, hydromet_tol, & ! intent(in)
396 0 : hydromet(j,:,:), mixt_frac(j,:), precip_frac(j,:), & ! in
397 0 : precip_frac_1(j,:), precip_frac_2(j,:), & ! intent(in)
398 0 : precip_frac_tol(j) ) ! intent(in)
399 : end do
400 : endif
401 :
402 0 : return
403 :
404 : end subroutine precip_fraction
405 :
406 : !=============================================================================
407 : subroutine component_precip_frac_weighted( nz, ngrdcol, hydromet_dim, &
408 : l_frozen_hm, hydromet_tol, &
409 : hydromet, precip_frac, &
410 : cloud_frac_1, cloud_frac_2, &
411 : ice_supersat_frac_1, &
412 : ice_supersat_frac_2, mixt_frac, &
413 : precip_frac_tol, &
414 : precip_frac_1, precip_frac_2 )
415 :
416 : ! Description:
417 : ! Set precipitation fraction in each component of the PDF. The weighted 1st
418 : ! PDF component precipitation fraction (weighted_pfrac_1) at a grid level is
419 : ! calculated by the greatest value of mixt_frac * cloud_frac_1 at or above
420 : ! the relevant grid level. Likewise, the weighted 2nd PDF component
421 : ! precipitation fraction (weighted_pfrac_2) at a grid level is calculated by
422 : ! the greatest value of ( 1 - mixt_frac ) * cloud_frac_2 at or above the
423 : ! relevant grid level.
424 : !
425 : ! The fraction weighted_pfrac_1 / ( weighted_pfrac_1 + weighted_pfrac_2 ) is
426 : ! the weighted_pfrac_1 fraction. Multiplying this fraction by overall
427 : ! precipitation fraction and then dividing by mixt_frac produces the 1st PDF
428 : ! component precipitation fraction (precip_frac_1). Then, calculate the 2nd
429 : ! PDF component precipitation fraction (precip_frac_2) accordingly.
430 :
431 : ! References:
432 : !-----------------------------------------------------------------------
433 :
434 : use constants_clubb, only: &
435 : one, & ! Constant(s)
436 : zero
437 :
438 : use clubb_precision, only: &
439 : core_rknd ! Variable(s)
440 :
441 : implicit none
442 :
443 : ! Input Variables
444 : integer, intent(in) :: &
445 : nz, & ! Number of model vertical grid levels
446 : ngrdcol, & ! Number of grid columns
447 : hydromet_dim ! Number of hydrometeor species
448 :
449 : logical, dimension(hydromet_dim), intent(in) :: &
450 : l_frozen_hm ! if true, then the hydrometeor is frozen; otherwise liquid
451 :
452 : real( kind = core_rknd ), dimension(hydromet_dim), intent(in) :: &
453 : hydromet_tol ! Tolerance values for all hydrometeors [units vary]
454 :
455 : real( kind = core_rknd ), dimension(ngrdcol,nz,hydromet_dim), intent(in) :: &
456 : hydromet ! Mean of hydrometeor, hm (overall) [units vary]
457 :
458 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
459 : precip_frac, & ! Precipitation fraction (overall) [-]
460 : cloud_frac_1, & ! Cloud fraction (1st PDF component) [-]
461 : cloud_frac_2, & ! Cloud fraction (2nd PDF component) [-]
462 : ice_supersat_frac_1, & ! Ice supersaturation fraction (1st PDF comp.) [-]
463 : ice_supersat_frac_2, & ! Ice supersaturation fraction (2nd PDF comp.) [-]
464 : mixt_frac ! Mixture fraction [-]
465 :
466 : real( kind = core_rknd ), dimension(ngrdcol), intent(in) :: &
467 : precip_frac_tol ! Minimum precip. frac. when hydromet. are present [-]
468 :
469 : ! Output Variables
470 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
471 : precip_frac_1, & ! Precipitation fraction (1st PDF component) [-]
472 : precip_frac_2 ! Precipitation fraction (2nd PDF component) [-]
473 :
474 : ! Local Variables
475 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
476 : weighted_pfrac_1, & ! Product of mixt_frac and cloud_frac_1 [-]
477 : weighted_pfrac_2 ! Product of ( 1 - mixt_frac ) and cloud_frac_2 [-]
478 :
479 : integer :: k, j ! Loop index
480 :
481 :
482 : !!! Find precipitation fraction within PDF component 1.
483 : ! The method used to find overall precipitation fraction will also be to
484 : ! find precipitation fraction within PDF component 1 and PDF component 2.
485 : ! In order to do so, it is assumed (poorly) that PDF component 1 overlaps
486 : ! PDF component 1 and that PDF component 2 overlaps PDF component 2 at every
487 : ! vertical level in the vertical profile.
488 :
489 : ! The weighted precipitation fraction in PDF component 1 is the greatest
490 : ! value of the product of mixture fraction (mixt_frac) and 1st PDF
491 : ! component cloud fraction at or above a vertical level. Likewise, the
492 : ! weighted precipitation fraction in PDF component 2 is the greatest
493 : ! value of the product of ( 1 - mixt_frac ) and 2nd PDF component cloud
494 : ! fraction at or above a vertical level.
495 : if ( any( l_frozen_hm ) ) then
496 :
497 : ! Ice microphysics included.
498 :
499 : ! Weighted precipitation fraction in PDF component 1.
500 : weighted_pfrac_1(:,nz) &
501 : = max( mixt_frac(:,nz) * cloud_frac_1(:,nz), &
502 : mixt_frac(:,nz) * ice_supersat_frac_1(:,nz) )
503 :
504 : ! Weighted precipitation fraction in PDF component 2.
505 : weighted_pfrac_2(:,nz) &
506 : = max( ( one - mixt_frac(:,nz) ) * cloud_frac_2(:,nz), &
507 : ( one - mixt_frac(:,nz) ) * ice_supersat_frac_2(:,nz) )
508 :
509 : do k = nz, 1, -1
510 : do j = 1, ngrdcol
511 : ! Weighted precipitation fraction in PDF component 1.
512 : weighted_pfrac_1(j,k) &
513 : = max( weighted_pfrac_1(j,k+1), &
514 : mixt_frac(j,k) * cloud_frac_1(j,k), &
515 : mixt_frac(j,k) * ice_supersat_frac_1(j,k) )
516 :
517 : ! Weighted precipitation fraction in PDF component 2.
518 : weighted_pfrac_2(j,k) &
519 : = max( weighted_pfrac_2(j,k+1), &
520 : ( one - mixt_frac(j,k) ) * cloud_frac_2(j,k), &
521 : ( one - mixt_frac(j,k) ) * ice_supersat_frac_2(j,k) )
522 : end do
523 : end do
524 :
525 :
526 : else
527 :
528 : ! Warm microphysics.
529 :
530 : ! Weighted precipitation fraction in PDF component 1.
531 : weighted_pfrac_1(:,nz) = mixt_frac(:,nz) * cloud_frac_1(:,nz)
532 :
533 : ! Weighted precipitation fraction in PDF component 2.
534 : weighted_pfrac_2(:,nz) = ( one - mixt_frac(:,nz) ) * cloud_frac_2(:,nz)
535 :
536 : do k = nz, 1, -1
537 : do j = 1, ngrdcol
538 : ! Weighted precipitation fraction in PDF component 1.
539 : weighted_pfrac_1(j,k) &
540 : = max( weighted_pfrac_1(j,k+1), &
541 : mixt_frac(j,k) * cloud_frac_1(j,k) )
542 :
543 : ! Weighted precipitation fraction in PDF component 2.
544 : weighted_pfrac_2(j,k) &
545 : = max( weighted_pfrac_2(j,k+1), &
546 : ( one - mixt_frac(j,k) ) * cloud_frac_2(j,k) )
547 : end do
548 : end do
549 :
550 : end if
551 :
552 :
553 : ! Calculate precip_frac_1 and special cases for precip_frac_1.
554 : do k = 1, nz, 1
555 : do j = 1, ngrdcol
556 : ! Calculate precipitation fraction in the 1st PDF component.
557 : if ( weighted_pfrac_1(j,k) + weighted_pfrac_2(j,k) > zero ) then
558 :
559 : ! Adjust weighted 1st PDF component precipitation fraction by
560 : ! multiplying it by a factor. That factor is overall precipitation
561 : ! fraction divided by the sum of the weighted 1st PDF component
562 : ! precipitation fraction and the weighted 2nd PDF component
563 : ! precipitation fraction. The 1st PDF component precipitation
564 : ! fraction is then found by dividing the adjusted weighted 1st PDF
565 : ! component precipitation fraction by mixture fraction.
566 : precip_frac_1(j,k) &
567 : = weighted_pfrac_1(j,k) &
568 : * ( precip_frac(j,k) &
569 : / ( weighted_pfrac_1(j,k) + weighted_pfrac_2(j,k) ) ) &
570 : / mixt_frac(j,k)
571 : else
572 :
573 : ! Usually, the sum of the weighted 1st PDF component precipitation
574 : ! fraction and the 2nd PDF component precipitation fraction go to 0
575 : ! when overall precipitation fraction goes to 0. Since 1st PDF
576 : ! component weighted precipitation fraction is 0, 1st PDF component
577 : ! precipitation fraction also 0.
578 : precip_frac_1(j,k) = zero
579 :
580 : end if
581 : end do
582 : end do
583 :
584 : do k = 1, nz, 1
585 : do j = 1, ngrdcol
586 : ! Special cases for precip_frac_1.
587 : if ( any( hydromet(j,k,:) >= hydromet_tol(:) ) &
588 : .and. precip_frac_1(j,k) &
589 : > min( one, precip_frac(j,k) / mixt_frac(j,k) ) ) then
590 :
591 : ! Using the above method, it is possible for precip_frac_1 to be
592 : ! greater than 1. For example, the mixture fraction at level k+1 is
593 : ! 0.10 and the cloud_frac_1 at level k+1 is 1, resulting in a
594 : ! weighted_pfrac_1 of 0.10. This product is greater than the product
595 : ! of mixt_frac and cloud_frac_1 at level k. The mixture fraction at
596 : ! level k is 0.05, resulting in a precip_frac_1 of 2. The value of
597 : ! precip_frac_1 is limited at 1. The leftover precipitation fraction
598 : ! (a result of the decreasing weight of PDF component 1 between the
599 : ! levels) is applied to PDF component 2.
600 : ! Additionally, when weighted_pfrac_1 at level k is greater than
601 : ! overall precipitation fraction at level k, the resulting calculation
602 : ! of precip_frac_2 at level k will be negative.
603 : precip_frac_1(j,k) = min( one, precip_frac(j,k) / mixt_frac(j,k) )
604 :
605 : elseif ( any( hydromet(j,k,:) >= hydromet_tol(:) ) &
606 : .and. precip_frac_1(j,k) > zero &
607 : .and. precip_frac_1(j,k) < precip_frac_tol(j) ) then
608 :
609 : ! In a scenario where we find precipitation in the 1st PDF component
610 : ! (it is allowed to have a value of 0 when all precipitation is found
611 : ! in the 2nd PDF component) but it is tiny (less than tolerance
612 : ! level), boost 1st PDF component precipitation fraction to tolerance
613 : ! level.
614 : precip_frac_1(j,k) = precip_frac_tol(j)
615 :
616 : elseif ( all( hydromet(j,k,:) < hydromet_tol(:) ) ) then
617 :
618 : ! The means (overall) of every precipitating hydrometeor are all less
619 : ! than their respective tolerance amounts. They are all considered to
620 : ! have values of 0. There are not any hydrometeor species found at
621 : ! this grid level. There is also no cloud at or above this grid
622 : ! level, so set 1st component precipitation fraction to 0.
623 : precip_frac_1(j,k) = zero
624 :
625 : end if
626 : end do
627 : end do
628 :
629 :
630 : !!! Find precipitation fraction within PDF component 2.
631 : ! The equation for precipitation fraction within PDF component 2 is:
632 : !
633 : ! f_p(2) = ( f_p - a * f_p(1) ) / ( 1 - a );
634 : !
635 : ! given the overall precipitation fraction, f_p (calculated above), the
636 : ! precipitation fraction within PDF component 1, f_p(1) (calculated above),
637 : ! and mixture fraction, a. Any leftover precipitation fraction from
638 : ! precip_frac_1 will be included in this calculation of precip_frac_2.
639 : do k = 1, nz, 1
640 : do j = 1, ngrdcol
641 :
642 : if ( any( hydromet(j,k,:) >= hydromet_tol(:) ) ) then
643 :
644 : ! Calculate precipitation fraction in the 2nd PDF component.
645 : precip_frac_2(j,k) &
646 : = max( ( precip_frac(j,k) - mixt_frac(j,k) * precip_frac_1(j,k) ) &
647 : / ( one - mixt_frac(j,k) ), &
648 : zero )
649 :
650 : ! Special cases for precip_frac_2.
651 : if ( precip_frac_2(j,k) > one ) then
652 :
653 : ! Again, it is possible for precip_frac_2 to be greater than 1.
654 : ! For example, the mixture fraction at level k+1 is 0.10 and the
655 : ! cloud_frac_1 at level k+1 is 1, resulting in a weighted_pfrac_1
656 : ! of 0.10. This product is greater than the product of mixt_frac
657 : ! and cloud_frac_1 at level k. Additionally, precip_frac (overall)
658 : ! is 1 for level k. The mixture fraction at level k is 0.5,
659 : ! resulting in a precip_frac_1 of 0.2. Using the above equation,
660 : ! precip_frac_2 is calculated to be 1.8. The value of
661 : ! precip_frac_2 is limited at 1. The leftover precipitation
662 : ! fraction (as a result of the increasing weight of component 1
663 : ! between the levels) is applied to PDF component 1.
664 : precip_frac_2(j,k) = one
665 :
666 : ! Recalculate precipitation fraction in the 1st PDF component.
667 : precip_frac_1(j,k) &
668 : = ( precip_frac(j,k) - ( one - mixt_frac(j,k) ) ) / mixt_frac(j,k)
669 :
670 : ! Double check precip_frac_1
671 : if ( precip_frac_1(j,k) > one ) then
672 :
673 : precip_frac_1(j,k) = one
674 :
675 : precip_frac_2(j,k) = ( precip_frac(j,k) - mixt_frac(j,k) ) &
676 : / ( one - mixt_frac(j,k) )
677 :
678 : elseif ( precip_frac_1(j,k) > zero .and. precip_frac_1(j,k) < precip_frac_tol(j) ) then
679 :
680 : precip_frac_1(j,k) = precip_frac_tol(j)
681 :
682 : ! fp = a*fp1+(1-a)*fp2 solving for fp2
683 : precip_frac_2(j,k) = precip_frac_1(j,k) &
684 : * ( ( ( precip_frac(j,k) / precip_frac_1(j,k)) &
685 : - mixt_frac(j,k) ) / ( one - mixt_frac(j,k) ) )
686 :
687 : end if
688 :
689 : elseif ( precip_frac_2(j,k) > zero &
690 : .and. precip_frac_2(j,k) < precip_frac_tol(j) ) then
691 :
692 : ! In a scenario where we find precipitation in the 2nd PDF
693 : ! component (it is allowed to have a value of 0 when all
694 : ! precipitation is found in the 1st PDF component) but it is tiny
695 : ! (less than tolerance level), boost 2nd PDF component
696 : ! precipitation fraction to tolerance level.
697 : precip_frac_2(j,k) = precip_frac_tol(j)
698 :
699 : ! Recalculate precipitation fraction in the 1st PDF component.
700 : precip_frac_1(j,k) &
701 : = ( precip_frac(j,k) - ( one - mixt_frac(j,k) ) * precip_frac_2(j,k) ) &
702 : / mixt_frac(j,k)
703 :
704 : ! Double check precip_frac_1
705 : if ( precip_frac_1(j,k) > one ) then
706 :
707 : precip_frac_1(j,k) = one
708 :
709 : precip_frac_2(j,k) = ( precip_frac(j,k) - mixt_frac(j,k) ) &
710 : / ( one - mixt_frac(j,k) )
711 :
712 : elseif ( precip_frac_1(j,k) > zero .and. precip_frac_1(j,k) < precip_frac_tol(j) ) then
713 :
714 : precip_frac_1(j,k) = precip_frac_tol(j)
715 :
716 : ! fp = a*fp1+(1-a)*fp2 solving for fp2
717 : precip_frac_2(j,k) = precip_frac_1(j,k) &
718 : * ( ( ( precip_frac(j,k) / precip_frac_1(j,k)) &
719 : - mixt_frac(j,k) ) / ( one - mixt_frac(j,k) ) )
720 :
721 : end if
722 :
723 : end if ! Special cases for precip_frac_2
724 :
725 : else ! all( hydromet(k,:) < hydromet_tol(:) )
726 :
727 : ! The means (overall) of every precipitating hydrometeor are all less
728 : ! than their respective tolerance amounts. They are all considered to
729 : ! have values of 0. There are not any hydrometeor species found at
730 : ! this grid level. There is also no cloud at or above this grid
731 : ! level, so set 2nd component precipitation fraction to 0.
732 : precip_frac_2(j,k) = zero
733 :
734 : end if ! any( hydromet(k,:) > hydromet_tol(:) )
735 : end do
736 : end do
737 :
738 : return
739 :
740 : end subroutine component_precip_frac_weighted
741 :
742 : !=============================================================================
743 0 : subroutine component_precip_frac_specify( nz, ngrdcol, hydromet_dim, hydromet_tol, &
744 : upsilon_precip_frac_rat, &
745 0 : hydromet, precip_frac, &
746 0 : mixt_frac, precip_frac_tol, &
747 0 : precip_frac_1, precip_frac_2 )
748 :
749 : ! Description:
750 : ! Calculates the precipitation fraction in each PDF component.
751 : !
752 : ! The equation for precipitation fraction is:
753 : !
754 : ! f_p = mixt_frac * f_p(1) + ( 1 - mixt_frac ) * f_p(2);
755 : !
756 : ! where f_p is overall precipitation fraction, f_p(1) is precipitation
757 : ! fraction in the 1st PDF component, f_p(2) is precipitation fraction in the
758 : ! 2nd PDF component, and mixt_frac is the mixture fraction. Using this
759 : ! method, a new specified parameter is introduced, upsilon, where:
760 : !
761 : ! upsilon = mixt_frac * f_p(1) / f_p; and where 0 <= upsilon <= 1.
762 : !
763 : ! In other words, upsilon is the ratio of mixt_frac * f_p(1) to f_p. Since
764 : ! f_p and mixt_frac are calculated previously, and upsilon is specified,
765 : ! f_p(1) can be calculated by:
766 : !
767 : ! f_p(1) = upsilon * f_p / mixt_frac;
768 : !
769 : ! and has an upper limit of 1. The value of f_p(2) can then be calculated
770 : ! by:
771 : !
772 : ! f_p(2) = ( f_p - mixt_frac * f_p(1) ) / ( 1 - mixt_frac );
773 : !
774 : ! and also has an upper limit of 1. When upsilon = 1, all of the
775 : ! precipitation is found in the 1st PDF component (as long as
776 : ! f_p <= mixt_frac, otherwise it would cause f_p(1) to be greater than 1).
777 : ! When upsilon = 0, all of the precipitation is found in the 2nd PDF
778 : ! component (as long as f_p <= 1 - mixt_frac, otherwise it would cause
779 : ! f_p(2) to be greater than 1). When upsilon is between 0 and 1,
780 : ! precipitation is split between the two PDF components accordingly.
781 :
782 : ! References:
783 : !-----------------------------------------------------------------------
784 :
785 : use constants_clubb, only: &
786 : one, & ! Constant(s)
787 : zero, &
788 : eps
789 :
790 : use clubb_precision, only: &
791 : core_rknd ! Variable(s)
792 :
793 : implicit none
794 :
795 : ! Input Variables
796 : integer, intent(in) :: &
797 : nz, & ! Number of model vertical grid levels
798 : ngrdcol, & ! Number of grid columns
799 : hydromet_dim ! Number of hydrometeor species
800 :
801 : real( kind = core_rknd ), dimension(hydromet_dim), intent(in) :: &
802 : hydromet_tol ! Tolerance values for all hydrometeors [units vary]
803 :
804 : real( kind = core_rknd ), intent(in) :: &
805 : upsilon_precip_frac_rat ! ratio mixt_frac*precip_frac_1/precip_frac [-]
806 :
807 : real( kind = core_rknd ), dimension(ngrdcol,nz,hydromet_dim), intent(in) :: &
808 : hydromet ! Mean of hydrometeor, hm (overall) [units vary]
809 :
810 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
811 : precip_frac, & ! Precipitation fraction (overall) [-]
812 : mixt_frac ! Mixture fraction [-]
813 :
814 : real( kind = core_rknd ), dimension(ngrdcol), intent(in) :: &
815 : precip_frac_tol ! Minimum precip. frac. when hydromet. are present [-]
816 :
817 : ! Output Variables
818 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
819 : precip_frac_1, & ! Precipitation fraction (1st PDF component) [-]
820 : precip_frac_2 ! Precipitation fraction (2nd PDF component) [-]
821 :
822 : integer :: k, j ! Loop index.
823 :
824 :
825 : ! There are hydrometeors found at this grid level.
826 0 : if ( abs(upsilon_precip_frac_rat - one) < &
827 : abs(upsilon_precip_frac_rat + one) / 2 * eps ) then
828 :
829 :
830 : ! Loop over all vertical levels.
831 0 : do k = 1, nz, 1
832 0 : do j = 1, ngrdcol
833 :
834 0 : if ( any( hydromet(j,k,:) >= hydromet_tol(:) ) ) then
835 :
836 0 : if ( precip_frac(j,k) <= mixt_frac(j,k) ) then
837 :
838 : ! All the precipitation is found in the 1st PDF component.
839 0 : precip_frac_1(j,k) = precip_frac(j,k) / mixt_frac(j,k)
840 0 : precip_frac_2(j,k) = zero
841 :
842 : else ! precip_frac(k) > mixt_frac(k)
843 :
844 : ! Some precipitation is found in the 2nd PDF component.
845 0 : precip_frac_1(j,k) = one
846 : precip_frac_2(j,k) = ( precip_frac(j,k) - mixt_frac(j,k) ) &
847 0 : / ( one - mixt_frac(j,k) )
848 :
849 : if ( precip_frac_2(j,k) > one &
850 0 : .and. abs(precip_frac(j,k) - one) < abs(precip_frac(j,k) + one) / 2 * eps ) then
851 :
852 : ! Set precip_frac_2 = 1.
853 0 : precip_frac_2(j,k) = one
854 :
855 0 : elseif ( precip_frac_2(j,k) < precip_frac_tol(j) ) then
856 :
857 : ! Since precipitation is found in the 2nd PDF component, it
858 : ! must have a value of at least precip_frac_tol.
859 0 : precip_frac_2(j,k) = precip_frac_tol(j)
860 :
861 : ! Recalculate precip_frac_1
862 : precip_frac_1(j,k) &
863 : = ( precip_frac(j,k) &
864 : - ( one - mixt_frac(j,k) ) * precip_frac_2(j,k) ) &
865 0 : / mixt_frac(j,k)
866 :
867 : ! Double check precip_frac_1
868 0 : if ( precip_frac_1(j,k) > one ) then
869 :
870 0 : precip_frac_1(j,k) = one
871 :
872 : precip_frac_2(j,k) = ( precip_frac(j,k) - mixt_frac(j,k) ) &
873 0 : / ( one - mixt_frac(j,k) )
874 :
875 0 : elseif ( precip_frac_1(j,k) < precip_frac_tol(j) ) then
876 :
877 0 : precip_frac_1(j,k) = precip_frac_tol(j)
878 :
879 : ! fp = a*fp1+(1-a)*fp2 solving for fp2
880 : precip_frac_2(j,k) = precip_frac_1(j,k) * &
881 : ( ( ( precip_frac(j,k) / precip_frac_1(j,k)) &
882 0 : - mixt_frac(j,k) ) / ( one - mixt_frac(j,k) ) )
883 :
884 : end if
885 :
886 : end if ! precip_frac_2(k) < precip_frac_tol
887 :
888 : endif ! precip_frac(k) <= mixt_frac(k)
889 :
890 : else ! all( hydromet(k,:) < hydromet_tol(:) )
891 : ! There aren't any hydrometeors found at the grid level.
892 0 : precip_frac_1(j,k) = zero
893 0 : precip_frac_2(j,k) = zero
894 : end if
895 :
896 : end do
897 : end do
898 :
899 0 : elseif ( abs(upsilon_precip_frac_rat - zero) < &
900 : abs(upsilon_precip_frac_rat + zero) / 2 * eps ) then
901 :
902 0 : do k = 1, nz, 1
903 0 : do j = 1, ngrdcol
904 :
905 0 : if ( any( hydromet(j,k,:) >= hydromet_tol(:) ) ) then
906 :
907 0 : if ( precip_frac(j,k) <= ( one - mixt_frac(j,k) ) ) then
908 :
909 : ! All the precipitation is found in the 2nd PDF component.
910 0 : precip_frac_1(j,k) = zero
911 0 : precip_frac_2(j,k) = precip_frac(j,k) / ( one - mixt_frac(j,k) )
912 :
913 : else ! precip_frac(k) > ( 1 - mixt_frac(k) )
914 :
915 : ! Some precipitation is found in the 1st PDF component.
916 : precip_frac_1(j,k) = ( precip_frac(j,k) - ( one - mixt_frac(j,k) ) ) &
917 0 : / mixt_frac(j,k)
918 0 : precip_frac_2(j,k) = one
919 :
920 : if ( precip_frac_1(j,k) > one &
921 0 : .and. abs(precip_frac(j,k) - one) < abs(precip_frac(j,k) + one) / 2 * eps ) then
922 :
923 : ! Set precip_frac_1 = 1.
924 0 : precip_frac_1(j,k) = one
925 :
926 0 : elseif ( precip_frac_1(j,k) < precip_frac_tol(j) ) then
927 :
928 : ! Since precipitation is found in the 1st PDF component, it
929 : ! must have a value of at least precip_frac_tol.
930 0 : precip_frac_1(j,k) = precip_frac_tol(j)
931 :
932 : ! Recalculate precip_frac_2
933 : precip_frac_2(j,k) = ( precip_frac(j,k) &
934 : - mixt_frac(j,k) * precip_frac_1(j,k) ) &
935 0 : / ( one - mixt_frac(j,k) )
936 :
937 : ! Double check precip_frac_2
938 0 : if ( precip_frac_2(j,k) > one ) then
939 :
940 0 : precip_frac_2(j,k) = one
941 :
942 : precip_frac_1(j,k) = ( ( precip_frac(j,k) - one ) + mixt_frac(j,k) ) &
943 0 : / mixt_frac(j,k)
944 :
945 0 : elseif ( precip_frac_2(j,k) < precip_frac_tol(j) ) then
946 :
947 0 : precip_frac_2(j,k) = precip_frac_tol(j)
948 :
949 : ! fp = a*fp1+(1-a)*fp2 solving for fp1
950 : precip_frac_1(j,k) = ( precip_frac(j,k) - precip_frac_2(j,k) ) / mixt_frac(j,k) &
951 0 : + precip_frac_2(j,k)
952 :
953 : endif
954 :
955 : end if ! precip_frac_1(k) < precip_frac_tol
956 :
957 : endif ! precip_frac(k) <= ( 1 - mixt_frac(k) )
958 :
959 : else ! all( hydromet(k,:) < hydromet_tol(:) )
960 : ! There aren't any hydrometeors found at the grid level.
961 0 : precip_frac_1(j,k) = zero
962 0 : precip_frac_2(j,k) = zero
963 : end if
964 :
965 : end do
966 : end do
967 :
968 : else ! 0 < upsilon_precip_frac_rat < 1
969 :
970 0 : do k = 1, nz, 1
971 0 : do j = 1, ngrdcol
972 :
973 0 : if ( any( hydromet(j,k,:) >= hydromet_tol(:) ) ) then
974 : ! Precipitation is found in both PDF components. Each component
975 : ! must have a precipitation fraction that is at least
976 : ! precip_frac_tol and that does not exceed 1.
977 :
978 : ! Calculate precipitation fraction in the 1st PDF component.
979 : precip_frac_1(j,k) &
980 0 : = upsilon_precip_frac_rat * precip_frac(j,k) / mixt_frac(j,k)
981 :
982 : ! Special cases for precip_frac_1
983 0 : if ( precip_frac_1(j,k) > one ) then
984 0 : precip_frac_1(j,k) = one
985 0 : elseif ( precip_frac_1(j,k) < precip_frac_tol(j) ) then
986 0 : precip_frac_1(j,k) = precip_frac_tol(j)
987 : endif
988 :
989 : ! Calculate precipitation fraction in the 2nd PDF component.
990 : precip_frac_2(j,k) = ( precip_frac(j,k) &
991 : - mixt_frac(j,k) * precip_frac_1(j,k) ) &
992 0 : / ( one - mixt_frac(j,k) )
993 :
994 : ! Special case for precip_frac_2
995 0 : if ( precip_frac_2(j,k) > one ) then
996 :
997 : ! Set precip_frac_2 to 1.
998 0 : precip_frac_2(j,k) = one
999 :
1000 : ! Recalculate precipitation fraction in the 1st PDF component.
1001 : precip_frac_1(j,k) &
1002 0 : = ( precip_frac(j,k) - ( one - mixt_frac(j,k) ) ) / mixt_frac(j,k)
1003 :
1004 : ! Double check precip_frac_1
1005 0 : if ( precip_frac_1(j,k) > one ) then
1006 :
1007 0 : precip_frac_1(j,k) = one
1008 :
1009 : precip_frac_2(j,k) = ( precip_frac(j,k) - mixt_frac(j,k) ) &
1010 0 : / ( one - mixt_frac(j,k) )
1011 :
1012 0 : elseif ( precip_frac_1(j,k) < precip_frac_tol(j) ) then
1013 :
1014 0 : precip_frac_1(j,k) = precip_frac_tol(j)
1015 :
1016 : ! fp = a*fp1+(1-a)*fp2 solving for fp2
1017 : precip_frac_2(j,k) = precip_frac_1(j,k) &
1018 : * ( ( ( precip_frac(j,k) / precip_frac_1(j,k)) &
1019 0 : - mixt_frac(j,k) ) / ( one - mixt_frac(j,k) ) )
1020 :
1021 : endif
1022 :
1023 0 : elseif ( precip_frac_2(j,k) < precip_frac_tol(j) ) then
1024 :
1025 : ! Set precip_frac_2 to precip_frac_tol.
1026 0 : precip_frac_2(j,k) = precip_frac_tol(j)
1027 :
1028 : ! Recalculate precipitation fraction in the 1st PDF component.
1029 : precip_frac_1(j,k) = ( precip_frac(j,k) &
1030 : - ( one - mixt_frac(j,k) ) * precip_frac_2(j,k) ) &
1031 0 : / mixt_frac(j,k)
1032 :
1033 : ! Double check precip_frac_1
1034 0 : if ( precip_frac_1(j,k) > one ) then
1035 :
1036 0 : precip_frac_1(j,k) = one
1037 :
1038 : precip_frac_2(j,k) = ( precip_frac(j,k) - mixt_frac(j,k) ) &
1039 0 : / ( one - mixt_frac(j,k) )
1040 :
1041 0 : elseif ( precip_frac_1(j,k) < precip_frac_tol(j) ) then
1042 :
1043 0 : precip_frac_1(j,k) = precip_frac_tol(j)
1044 :
1045 : ! fp = a*fp1+(1-a)*fp2 solving for fp2
1046 : precip_frac_2(j,k) = precip_frac_1(j,k) &
1047 : * ( ( ( precip_frac(j,k) / precip_frac_1(j,k)) &
1048 0 : - mixt_frac(j,k) ) / ( one - mixt_frac(j,k) ) )
1049 :
1050 : end if
1051 :
1052 : end if ! Special cases for precip_frac_2
1053 :
1054 : else ! all( hydromet(k,:) < hydromet_tol(:) )
1055 : ! There aren't any hydrometeors found at the grid level.
1056 0 : precip_frac_1(j,k) = zero
1057 0 : precip_frac_2(j,k) = zero
1058 : end if
1059 :
1060 : end do
1061 : end do
1062 :
1063 : end if ! upsilon_precip_frac_rat
1064 :
1065 0 : return
1066 :
1067 : end subroutine component_precip_frac_specify
1068 :
1069 : !=============================================================================
1070 0 : subroutine precip_frac_assert_check( nz, hydromet_dim, hydromet_tol, &
1071 0 : hydromet, mixt_frac, precip_frac, &
1072 0 : precip_frac_1, precip_frac_2, &
1073 : precip_frac_tol )
1074 :
1075 : ! Description:
1076 : ! Assertion check for the precipitation fraction code.
1077 :
1078 : ! References:
1079 : !-----------------------------------------------------------------------
1080 :
1081 : use constants_clubb, only: &
1082 : one, & ! Constant(s)
1083 : zero, &
1084 : fstderr, &
1085 : eps
1086 :
1087 : use clubb_precision, only: &
1088 : core_rknd ! Variable(s)
1089 :
1090 : use error_code, only: &
1091 : err_code, & ! Error Indicator
1092 : clubb_fatal_error ! Constant
1093 :
1094 : implicit none
1095 :
1096 : ! Input Variables
1097 : integer, intent(in) :: &
1098 : nz, & ! Number of model vertical grid levels
1099 : hydromet_dim ! Number of hydrometeor species
1100 :
1101 : real( kind = core_rknd ), dimension(hydromet_dim), intent(in) :: &
1102 : hydromet_tol ! Tolerance values for all hydrometeors [units vary]
1103 :
1104 : real( kind = core_rknd ), dimension(nz,hydromet_dim), intent(in) :: &
1105 : hydromet ! Mean of hydrometeor, hm (overall) [units vary]
1106 :
1107 : real( kind = core_rknd ), dimension(nz), intent(in) :: &
1108 : mixt_frac, & ! Mixture fraction [-]
1109 : precip_frac, & ! Precipitation fraction (overall) [-]
1110 : precip_frac_1, & ! Precipitation fraction (1st PDF component) [-]
1111 : precip_frac_2 ! Precipitation fraction (2nd PDF component) [-]
1112 :
1113 : real( kind = core_rknd ), intent(in) :: &
1114 : precip_frac_tol ! Minimum precip. frac. when hydromet. are present [-]
1115 :
1116 : ! Local Variables
1117 : integer :: k ! Loop index
1118 :
1119 :
1120 : ! Loop over all vertical levels.
1121 0 : do k = 1, nz, 1
1122 :
1123 0 : if ( any( hydromet(k,:) >= hydromet_tol(:) ) ) then
1124 :
1125 : ! Overall precipitation fraction cannot be less than precip_frac_tol
1126 : ! when a hydrometeor is present at a grid level.
1127 0 : if ( precip_frac(k) < precip_frac_tol ) then
1128 : write(fstderr,*) "precip_frac < precip_frac_tol when " &
1129 0 : // "a hydrometeor is present"
1130 0 : write(fstderr,*) "level = ", k
1131 0 : write(fstderr,*) "precip_frac = ", precip_frac(k), &
1132 0 : "precip_frac_tol = ", precip_frac_tol
1133 0 : err_code = clubb_fatal_error
1134 0 : return
1135 : endif
1136 :
1137 : ! Overall precipitation fraction cannot exceed 1.
1138 0 : if ( precip_frac(k) > one ) then
1139 0 : write(fstderr,*) "precip_frac > 1"
1140 0 : write(fstderr,*) "level = ", k
1141 0 : write(fstderr,*) "precip_frac = ", precip_frac(k)
1142 0 : err_code = clubb_fatal_error
1143 0 : return
1144 : endif
1145 :
1146 : ! Precipitation fraction in the 1st PDF component is allowed to be 0
1147 : ! when all the precipitation is found in the 2nd PDF component.
1148 : ! Otherwise, it cannot be less than precip_frac_tol when a hydrometeor
1149 : ! is present at a grid level. In other words, it cannot have a value
1150 : ! that is greater than 0 but less than precip_frac_tol
1151 : if ( precip_frac_1(k) > zero &
1152 0 : .and. precip_frac_1(k) < precip_frac_tol ) then
1153 0 : write(fstderr,*) "0 < precip_frac_1 < precip_frac_tol"
1154 0 : write(fstderr,*) "level = ", k
1155 0 : write(fstderr,*) "precip_frac_1 = ", precip_frac_1(k), &
1156 0 : "precip_frac_tol = ", precip_frac_tol
1157 0 : err_code = clubb_fatal_error
1158 0 : return
1159 : endif
1160 :
1161 : ! Precipitation fraction in the 1st PDF component cannot exceed 1.
1162 0 : if ( precip_frac_1(k) > one ) then
1163 0 : write(fstderr,*) "precip_frac_1 > 1"
1164 0 : write(fstderr,*) "level = ", k
1165 0 : write(fstderr,*) "precip_frac_1 = ", precip_frac_1(k)
1166 0 : err_code = clubb_fatal_error
1167 0 : return
1168 : endif
1169 :
1170 : ! Precipiation fraction in the 1st PDF component cannot be negative.
1171 0 : if ( precip_frac_1(k) < zero ) then
1172 0 : write(fstderr,*) "precip_frac_1 < 0"
1173 0 : write(fstderr,*) "level = ", k
1174 0 : write(fstderr,*) "precip_frac_1 = ", precip_frac_1(k)
1175 0 : err_code = clubb_fatal_error
1176 0 : return
1177 : endif
1178 :
1179 : ! Precipitation fraction in the 2nd PDF component is allowed to be 0
1180 : ! when all the precipitation is found in the 1st PDF component.
1181 : ! Otherwise, it cannot be less than precip_frac_tol when a hydrometeor
1182 : ! is present at a grid level. In other words, it cannot have a value
1183 : ! that is greater than 0 but less than precip_frac_tol
1184 : if ( precip_frac_2(k) > zero &
1185 0 : .and. precip_frac_2(k) < precip_frac_tol ) then
1186 0 : write(fstderr,*) "0 < precip_frac_2 < precip_frac_tol"
1187 0 : write(fstderr,*) "level = ", k
1188 0 : write(fstderr,*) "precip_frac_2 = ", precip_frac_2(k), &
1189 0 : "precip_frac_tol = ", precip_frac_tol
1190 0 : err_code = clubb_fatal_error
1191 0 : return
1192 : endif
1193 :
1194 : ! Precipitation fraction in the 2nd PDF component cannot exceed 1.
1195 0 : if ( precip_frac_2(k) > one ) then
1196 0 : write(fstderr,*) "precip_frac_2 > 1"
1197 0 : write(fstderr,*) "level = ", k
1198 0 : write(fstderr,*) "precip_frac_2 = ", precip_frac_2(k)
1199 0 : err_code = clubb_fatal_error
1200 0 : return
1201 : endif
1202 :
1203 : ! Precipiation fraction in the 2nd PDF component cannot be negative.
1204 0 : if ( precip_frac_2(k) < zero ) then
1205 0 : write(fstderr,*) "precip_frac_2 < 0"
1206 0 : write(fstderr,*) "level = ", k
1207 0 : write(fstderr,*) "precip_frac_2 = ", precip_frac_2(k)
1208 0 : err_code = clubb_fatal_error
1209 0 : return
1210 : endif
1211 :
1212 : else ! all( hydromet(k,:) < hydromet_tol(:) )
1213 :
1214 : ! Overall precipitation fraction must be 0 when no hydrometeors are
1215 : ! found at a grid level.
1216 0 : if ( abs(precip_frac(k)) > eps) then
1217 0 : write(fstderr,*) "precip_frac /= 0 when no hydrometeors are found"
1218 0 : write(fstderr,*) "level = ", k
1219 0 : write(fstderr,*) "precip_frac = ", precip_frac(k)
1220 0 : err_code = clubb_fatal_error
1221 0 : return
1222 : endif
1223 :
1224 : ! Precipitation fraction in the 1st PDF component must be 0 when no
1225 : ! hydrometeors are found at a grid level.
1226 0 : if ( abs(precip_frac_1(k)) > eps) then
1227 : write(fstderr,*) "precip_frac_1 /= 0 when no hydrometeors " &
1228 0 : // "are found"
1229 0 : write(fstderr,*) "level = ", k
1230 0 : write(fstderr,*) "precip_frac_1 = ", precip_frac_1(k)
1231 0 : err_code = clubb_fatal_error
1232 0 : return
1233 : endif
1234 :
1235 : ! Precipitation fraction in the 2nd PDF component must be 0 when no
1236 : ! hydrometeors are found at a grid level.
1237 0 : if ( abs(precip_frac_2(k)) > eps) then
1238 : write(fstderr,*) "precip_frac_2 /= 0 when no hydrometeors " &
1239 0 : // "are found"
1240 0 : write(fstderr,*) "level = ", k
1241 0 : write(fstderr,*) "precip_frac_2 = ", precip_frac_2(k)
1242 0 : err_code = clubb_fatal_error
1243 0 : return
1244 : endif
1245 :
1246 : endif ! any( hydromet(k,:) >= hydromet_tol(:) )
1247 :
1248 : ! The precipitation fraction equation is:
1249 : !
1250 : ! precip_frac
1251 : ! = mixt_frac * precip_frac_1 + ( 1 - mixt_frac ) * precip_frac_2;
1252 : !
1253 : ! which means that:
1254 : !
1255 : ! precip_frac
1256 : ! - ( mixt_frac * precip_frac_1 + ( 1 - mixt_frac ) * precip_frac_2 )
1257 : ! = 0.
1258 : !
1259 : ! Check that this is true with numerical round off.
1260 0 : if ( ( precip_frac(k) &
1261 : - ( mixt_frac(k) * precip_frac_1(k) &
1262 : + ( one - mixt_frac(k) ) * precip_frac_2(k) ) ) &
1263 0 : > ( epsilon( precip_frac(k) ) * precip_frac(k) ) ) then
1264 : write(fstderr,*) "mixt_frac * precip_frac_1 " &
1265 : // "+ ( 1 - mixt_frac ) * precip_frac_2 " &
1266 0 : // "/= precip_frac within numerical roundoff"
1267 0 : write(fstderr,*) "level = ", k
1268 : write(fstderr,*) "mixt_frac * precip_frac_1 " &
1269 0 : // "+ ( 1 - mixt_frac ) * precip_frac_2 = ", &
1270 0 : mixt_frac(k) * precip_frac_1(k) &
1271 0 : + ( one - mixt_frac(k) ) * precip_frac_2(k)
1272 0 : write(fstderr,*) "precip_frac = ", precip_frac(k)
1273 0 : err_code = clubb_fatal_error
1274 0 : return
1275 : endif
1276 :
1277 : enddo ! k = 1, nz, 1
1278 :
1279 :
1280 : return
1281 :
1282 : end subroutine precip_frac_assert_check
1283 :
1284 : !===============================================================================
1285 :
1286 : end module precipitation_fraction
|