Line data Source code
1 : !-------------------------------------------------------------------------
2 : ! $Id$
3 : !===============================================================================
4 : module setup_clubb_pdf_params
5 :
6 : implicit none
7 :
8 : private
9 :
10 : public :: setup_pdf_parameters, &
11 : compute_mean_stdev, &
12 : calc_comp_mu_sigma_hm, &
13 : norm_transform_mean_stdev, &
14 : comp_corr_norm, &
15 : denorm_transform_corr
16 :
17 : private :: component_corr_w_x, &
18 : component_corr_chi_eta, &
19 : component_corr_w_hm_n_ip, &
20 : component_corr_x_hm_n_ip, &
21 : component_corr_hmx_hmy_n_ip, &
22 : component_corr_eta_hm_n_ip, &
23 : calc_corr_w_hm_n, &
24 : pdf_param_hm_stats, &
25 : pdf_param_ln_hm_stats, &
26 : pack_hydromet_pdf_params, &
27 : compute_rtp2_from_chi
28 :
29 : ! The component correlations of w and r_t and the component correlations of
30 : ! w and theta_l must both be 0 when using the ADG1 PDF. In other words, w
31 : ! and r_t (theta_l) have overall covariance w'r_t' (w'theta_l'), but the
32 : ! individual component covariance and correlation are defined to be 0.
33 : ! Since the component covariances (or correlations) of w and chi (old s) and
34 : ! of w and eta (old t) are based on the covariances (or correlations) of w
35 : ! and r_t and of w and theta_l, the individual component correlation and
36 : ! covariance of w and chi, as well of as w and eta, are defined to be 0.
37 : !
38 : ! The PDF component correlation of w and x, calculated as part of the PDF
39 : ! parameters, comes out to be 0 (within numerical roundoff) when the ADG1
40 : ! PDF is used. However, when l_fix_w_chi_eta_correlations is enabled, the
41 : ! component correlations of PDF variables are specified, and the values
42 : ! that were calculated as part of the PDF parameters are ignored in the
43 : ! correlation array. When this happens, and when the ADG1 PDF is selected,
44 : ! the l_follow_ADG1_PDF_standards flag can be used to force the component
45 : ! correlations of w and x to have a value of 0, following the ADG1 standard,
46 : ! rather than whatever value might be specified in the correlation array.
47 : !
48 : ! Note: the ADG2 PDF follows the same standards as the ADG1 PDF.
49 : logical, parameter :: &
50 : l_follow_ADG1_PDF_standards = .true.
51 :
52 : contains
53 :
54 : !=============================================================================
55 0 : subroutine setup_pdf_parameters( gr, nz, ngrdcol, pdf_dim, hydromet_dim, dt, & ! Intent(in)
56 0 : Nc_in_cloud, cloud_frac, Kh_zm, & ! Intent(in)
57 0 : ice_supersat_frac, hydromet, wphydrometp, & ! Intent(in)
58 0 : corr_array_n_cloud, corr_array_n_below, & ! Intent(in)
59 : hm_metadata, & ! Intent(in)
60 : pdf_params, & ! Intent(in)
61 : clubb_params, & ! Intent(in)
62 : iiPDF_type, & ! Intent(in)
63 : l_use_precip_frac, & ! Intent(in)
64 : l_predict_upwp_vpwp, & ! Intent(in)
65 : l_diagnose_correlations, & ! Intent(in)
66 : l_calc_w_corr, & ! Intent(in)
67 : l_const_Nc_in_cloud, & ! Intent(in)
68 : l_fix_w_chi_eta_correlations, & ! Intent(in)
69 : stats_metadata, & ! Intent(in)
70 0 : stats_zt, stats_zm, stats_sfc, & ! intent(inout)
71 0 : hydrometp2, & ! Intent(out)
72 0 : mu_x_1_n, mu_x_2_n, & ! Intent(out)
73 0 : sigma_x_1_n, sigma_x_2_n, & ! Intent(out)
74 0 : corr_array_1_n, corr_array_2_n, & ! Intent(out)
75 0 : corr_cholesky_mtx_1, corr_cholesky_mtx_2, & ! Intent(out)
76 : precip_fracs, & ! Intent(out)
77 0 : hydromet_pdf_params ) ! Intent(out)
78 :
79 : ! Description:
80 :
81 : ! References:
82 : !-----------------------------------------------------------------------
83 :
84 : use grid_class, only: &
85 : grid, & ! Type
86 : zm2zt, & ! Procedure(s)
87 : zt2zm
88 :
89 : use constants_clubb, only: &
90 : one, & ! Constant(s)
91 : zero, &
92 : rc_tol, &
93 : Ncn_tol, &
94 : cloud_frac_min, &
95 : fstderr, &
96 : zero_threshold
97 :
98 : use pdf_parameter_module, only: &
99 : pdf_parameter ! Variable(s)
100 :
101 : use hydromet_pdf_parameter_module, only: &
102 : hydromet_pdf_parameter, & ! Type
103 : precipitation_fractions, &
104 : init_hydromet_pdf_params ! Procedure
105 :
106 : use precipitation_fraction, only: &
107 : precip_fraction
108 :
109 : use Nc_Ncn_eqns, only: &
110 : Nc_in_cloud_to_Ncnm ! Procedure(s)
111 :
112 : use parameter_indices, only: &
113 : nparams, & ! Variable(s)
114 : ic_K_hm, &
115 : iomicron, &
116 : izeta_vrnce_rat
117 :
118 : use pdf_utilities, only: &
119 : calc_xp2, & ! Procedure(s)
120 : compute_mean_binormal, &
121 : compute_variance_binormal, &
122 : stdev_L2N, &
123 : corr_NN2NL
124 :
125 : use clip_explicit, only: &
126 : clip_covar_level, & ! Procedure(s)
127 : clip_wphydrometp ! Variables(s)
128 :
129 : use clubb_precision, only: &
130 : core_rknd ! Variable(s)
131 :
132 : use matrix_operations, only: &
133 : Cholesky_factor ! Procedure(s)
134 :
135 : use stats_type_utilities, only: &
136 : ! stat_update_var, & ! Procedure(s)
137 : stat_update_var_pt
138 :
139 : use diagnose_correlations_module, only: &
140 : diagnose_correlations, & ! Procedure(s)
141 : calc_cholesky_corr_mtx_approx
142 :
143 : use corr_varnce_module, only: &
144 : hm_metadata_type, & ! Type
145 : assert_corr_symmetric ! Procedure(s)
146 :
147 : use error_code, only: &
148 : clubb_at_least_debug_level, & ! Procedure
149 : err_code, & ! Error Indicator
150 : clubb_fatal_error ! Constant
151 :
152 : use stats_type, only: stats ! Type
153 :
154 : use advance_helper_module, only : &
155 : calc_xpwp ! Procedure(s)
156 :
157 : use stats_variables, only: &
158 : stats_metadata_type
159 :
160 : implicit none
161 :
162 : !------------------------ Input Variables ------------------------
163 : integer, intent(in) :: &
164 : nz, & ! Number of model vertical grid levels
165 : pdf_dim, & ! Number of variables in the correlation array
166 : ngrdcol, & ! Number of grid columns
167 : hydromet_dim ! Number of hydrometeor species
168 :
169 : type (grid), target, intent(in) :: gr
170 :
171 : real( kind = core_rknd ), intent(in) :: &
172 : dt ! Model timestep [s]
173 :
174 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
175 : Nc_in_cloud, & ! Mean (in-cloud) cloud droplet conc. [num/kg]
176 : cloud_frac, & ! Cloud fraction [-]
177 : Kh_zm, & ! Eddy diffusivity coef. on momentum levels [m^2/s]
178 : ice_supersat_frac ! Ice supersaturation fraction [-]
179 :
180 : real( kind = core_rknd ), dimension(ngrdcol,nz,hydromet_dim), intent(in) :: &
181 : hydromet, & ! Mean of hydrometeor, hm (overall) (t-levs.) [units]
182 : wphydrometp ! Covariance < w'h_m' > (momentum levels) [(m/s)units]
183 :
184 : real( kind = core_rknd ), dimension(pdf_dim,pdf_dim), intent(in) :: &
185 : corr_array_n_cloud, & ! Prescribed normal space corr. array in cloud [-]
186 : corr_array_n_below ! Prescribed normal space corr. array below cl. [-]
187 :
188 : type (hm_metadata_type), intent(in) :: &
189 : hm_metadata
190 :
191 : type(pdf_parameter), intent(in) :: &
192 : pdf_params ! PDF parameters [units vary]
193 :
194 : real( kind = core_rknd ), dimension(nparams), intent(in) :: &
195 : clubb_params ! Array of CLUBB's tunable parameters [units vary]
196 :
197 : integer, intent(in) :: &
198 : iiPDF_type ! Selected option for the two-component normal (double
199 : ! Gaussian) PDF type to use for the w, rt, and theta-l (or
200 : ! w, chi, and eta) portion of CLUBB's multivariate,
201 : ! two-component PDF.
202 :
203 : logical, intent(in) :: &
204 : l_use_precip_frac, & ! Flag to use precipitation fraction in KK microphysics. The
205 : ! precipitation fraction is automatically set to 1 when this
206 : ! flag is turned off.
207 : l_predict_upwp_vpwp, & ! Flag to predict <u'w'> and <v'w'> along with <u> and <v>
208 : ! alongside the advancement of <rt>, <w'rt'>, <thl>,
209 : ! <wpthlp>, <sclr>, and <w'sclr'> in subroutine
210 : ! advance_xm_wpxp. Otherwise, <u'w'> and <v'w'> are still
211 : ! approximated by eddy diffusivity when <u> and <v> are
212 : ! advanced in subroutine advance_windm_edsclrm.
213 : l_diagnose_correlations, & ! Diagnose correlations instead of using fixed ones
214 : l_calc_w_corr, & ! Calculate the correlations between w and the hydrometeors
215 : l_const_Nc_in_cloud, & ! Use a constant cloud droplet conc. within cloud (K&K)
216 : l_fix_w_chi_eta_correlations ! Use a fixed correlation for s and t Mellor(chi/eta)
217 :
218 : type (stats_metadata_type), intent(in) :: &
219 : stats_metadata
220 :
221 : !------------------------ Input/Output Variables ------------------------
222 : type (stats), target, dimension(ngrdcol), intent(inout) :: &
223 : stats_zt, &
224 : stats_zm, &
225 : stats_sfc
226 :
227 : !------------------------ Output Variables ------------------------
228 : real( kind = core_rknd ), dimension(ngrdcol,nz,hydromet_dim), intent(out) :: &
229 : hydrometp2 ! Variance of a hydrometeor (overall) (m-levs.) [units^2]
230 :
231 : real( kind = core_rknd ), dimension(ngrdcol,nz,pdf_dim), intent(out) :: &
232 : mu_x_1_n, & ! Mean array (normal space): PDF vars. (comp. 1) [un. vary]
233 : mu_x_2_n, & ! Mean array (normal space): PDF vars. (comp. 2) [un. vary]
234 : sigma_x_1_n, & ! Std. dev. array (normal space): PDF vars (comp. 1) [u.v.]
235 : sigma_x_2_n ! Std. dev. array (normal space): PDF vars (comp. 2) [u.v.]
236 :
237 : real( kind = core_rknd ), dimension(ngrdcol,nz,pdf_dim,pdf_dim), intent(out) :: &
238 : corr_array_1_n, & ! Corr. array (normal space) of PDF vars. (comp. 1) [-]
239 : corr_array_2_n ! Corr. array (normal space) of PDF vars. (comp. 2) [-]
240 :
241 : real( kind = core_rknd ), dimension(ngrdcol,nz,pdf_dim,pdf_dim), intent(out) :: &
242 : corr_cholesky_mtx_1, & ! Transposed corr. cholesky matrix, 1st comp. [-]
243 : corr_cholesky_mtx_2 ! Transposed corr. cholesky matrix, 2nd comp. [-]
244 :
245 : ! This is only an output, but it contains allocated arrays, so we need to treat it as inout
246 : type(precipitation_fractions), intent(inout) :: &
247 : precip_fracs ! Precipitation fractions [-]
248 :
249 : type(hydromet_pdf_parameter), dimension(ngrdcol,nz), optional, intent(out) :: &
250 : hydromet_pdf_params ! Hydrometeor PDF parameters [units vary]
251 :
252 : !------------------------ Local Variables ------------------------
253 :
254 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
255 0 : Kh_zm_c_K_hm ! Eddy diffusivity coef. on momentum levels [m^2/s]
256 :
257 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
258 0 : sigma_w_1, & ! Standard deviation of w (1st PDF component) [m/s]
259 0 : sigma_w_2 ! Standard deviation of w (2nd PDF component) [m/s]
260 :
261 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
262 : ice_supersat_frac_1, & ! Ice supersaturation fraction (1st PDF comp.) [-]
263 : ice_supersat_frac_2 ! Ice supersaturation fraction (2nd PDF comp.) [-]
264 :
265 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
266 0 : Ncnm ! Mean cloud nuclei concentration, < N_cn > [num/kg]
267 :
268 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
269 0 : rcm_pdf ! Liquid water mixing ratio [kg/kg]
270 :
271 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
272 0 : wpNcnp_zm, & ! Covariance of N_cn and w (momentum levs.) [(m/s)(num/kg)]
273 0 : wpNcnp_zt ! Covariance of N_cn and w on thermo. levels [(m/s)(num/kg)]
274 :
275 : real( kind = core_rknd ), dimension(ngrdcol,nz,hydromet_dim) :: &
276 0 : hm_1, & ! Mean of a precip. hydrometeor (1st PDF component) [units vary]
277 0 : hm_2 ! Mean of a precip. hydrometeor (2nd PDF component) [units vary]
278 :
279 : real( kind = core_rknd ), dimension(ngrdcol,nz,hydromet_dim) :: &
280 0 : hydrometp2_zt, & ! Variance of a hydrometeor (overall); t-lev [units^2]
281 0 : wphydrometp_zt ! Covariance of w and hm interp. to t-levs. [(m/s)units]
282 :
283 : real( kind = core_rknd ), dimension(ngrdcol,nz,pdf_dim,pdf_dim) :: &
284 0 : corr_array_1, & ! Correlation array of PDF vars. (comp. 1) [-]
285 0 : corr_array_2 ! Correlation array of PDF vars. (comp. 2) [-]
286 :
287 : real( kind = core_rknd ), dimension(ngrdcol,nz,pdf_dim) :: &
288 0 : mu_x_1, & ! Mean array of PDF vars. (1st PDF component) [units vary]
289 0 : mu_x_2, & ! Mean array of PDF vars. (2nd PDF component) [units vary]
290 0 : sigma_x_1, & ! Standard deviation array of PDF vars (comp. 1) [units vary]
291 0 : sigma_x_2 ! Standard deviation array of PDF vars (comp. 2) [units vary]
292 :
293 : real( kind = core_rknd ), dimension(pdf_dim) :: &
294 0 : corr_array_scaling
295 :
296 : real( kind = core_rknd ), dimension(ngrdcol,nz,pdf_dim) :: &
297 0 : sigma2_on_mu2_ip_1, & ! Ratio array sigma_hm_1^2/mu_hm_1^2 [-]
298 0 : sigma2_on_mu2_ip_2 ! Ratio array sigma_hm_2^2/mu_hm_2^2 [-]
299 :
300 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
301 0 : const_Ncnp2_on_Ncnm2, & ! Prescribed ratio: <Ncn'^2> to <Ncn>^2 [-]
302 0 : stdev_const_Ncnp2_on_Ncnm2, & ! Standard deviation of const_Ncnp2_on_Ncnm2
303 0 : const_corr_chi_Ncn, & ! Prescribed corr. of chi (old s) & Ncn [-]
304 0 : const_corr_chi_Ncn_n_cloud ! Prescribed corr. of chi & Ncn in cloud [-]
305 :
306 : real( kind = core_rknd ), dimension(ngrdcol) :: &
307 0 : precip_frac_tol ! Min. precip. frac. when hydromet. present [-]
308 :
309 : real( kind = core_rknd ), dimension(ngrdcol,nz,hydromet_dim) :: &
310 0 : wphydrometp_chnge ! Change in wphydrometp_zt: covar. clip. [(m/s)units]
311 :
312 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
313 0 : wm_zt, & ! Mean vertical velocity, <w>, on thermo. levels [m/s]
314 0 : wp2_zt ! Variance of w, <w'^2> (interp. to t-levs.) [m^2/s^2]
315 :
316 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
317 0 : rtp2_zt_from_chi, &
318 0 : rtp2_zm_from_chi
319 :
320 : real( kind = core_rknd ) :: &
321 : omicron, & ! Relative width parameter, omicron = R / Rmax [-]
322 : zeta_vrnce_rat ! Width parameter for sigma_hm_1^2 / mu_hm_1^2 [-]
323 :
324 : logical :: l_corr_array_scaling
325 :
326 : ! Flags used for covariance clipping of <w'hm'>.
327 : logical, parameter :: &
328 : l_first_clip_ts = .true., & ! First instance of clipping in a timestep.
329 : l_last_clip_ts = .true. ! Last instance of clipping in a timestep.
330 :
331 : character(len=10) :: &
332 : hydromet_name ! Name of a hydrometeor
333 :
334 : integer :: k, i, j ! Loop indices
335 :
336 : ! ---- Begin Code ----
337 :
338 : ! Assertion check
339 : ! Check that all hydrometeors are positive otherwise exit the program
340 0 : if ( clubb_at_least_debug_level( 0 ) ) then
341 0 : if ( any(hydromet < zero_threshold) ) then
342 0 : do j = 1, hydromet_dim
343 0 : hydromet_name = hm_metadata%hydromet_list(j)
344 0 : do k = 1, nz
345 0 : do i = 1, ngrdcol
346 0 : if ( hydromet(i,k,j) < zero_threshold ) then
347 : ! Write error message
348 0 : write(fstderr,*) " at beginning of setup_pdf_parameters: ", &
349 0 : trim( hydromet_name )//" = ", &
350 0 : hydromet(i,k,j), " < ", zero_threshold, &
351 0 : " at k = ", k
352 0 : err_code = clubb_fatal_error
353 0 : return
354 : end if ! hydromet(k,i) < 0
355 : end do ! k = 1, nz
356 : end do ! i = 1, hydromet_dim
357 : end do
358 : end if
359 : end if !clubb_at_least_debug_level( 0 )
360 :
361 : ! Setup some of the PDF parameters
362 0 : do k = 1, nz
363 0 : do i = 1, ngrdcol
364 0 : sigma_w_1(i,k) = sqrt( pdf_params%varnce_w_1(i,k) )
365 0 : sigma_w_2(i,k) = sqrt( pdf_params%varnce_w_2(i,k) )
366 : end do
367 : end do
368 :
369 : ! Compute rcm_pdf for use within SILHS
370 0 : do i = 1, ngrdcol
371 0 : rcm_pdf(i,:) = compute_mean_binormal( pdf_params%rc_1(i,:), pdf_params%rc_2(i,:), &
372 0 : pdf_params%mixt_frac(i,:) )
373 : end do
374 :
375 : ! Note on hydrometeor PDF shape:
376 : ! To use a single lognormal over the entire grid level, turn off the
377 : ! l_use_precip_frac flag and set omicron to 1 and zeta_vrnce_rat to 0 in
378 : ! tunable_parameters.in.
379 : ! To use a single delta-lognormal (single lognormal in-precip.), enable the
380 : ! l_use_precip_frac flag and set omicron to 1 and zeta_vrnce_rat to 0 in
381 : ! tunable_parameters.in.
382 : ! Otherwise, with l_use_precip_frac enabled and omicron and zeta_vrnce_rat
383 : ! values that are not 1 and 0, respectively, the PDF shape is a double
384 : ! delta-lognormal (two independent lognormals in-precip.).
385 :
386 : ! Calculate precipitation fraction.
387 0 : if ( l_use_precip_frac ) then
388 :
389 : call precip_fraction( nz, ngrdcol, hydromet_dim, & ! In
390 : hydromet(:,:,:), cloud_frac(:,:), pdf_params%cloud_frac_1(:,:), & ! In
391 : hm_metadata%l_mix_rat_hm, hm_metadata%l_frozen_hm, & ! In
392 : hm_metadata%hydromet_tol, & ! In
393 : pdf_params%cloud_frac_2(:,:), ice_supersat_frac(:,:), & ! In
394 : pdf_params%ice_supersat_frac_1(:,:), pdf_params%ice_supersat_frac_2(:,:), & ! In
395 : pdf_params%mixt_frac(:,:), clubb_params, & ! In
396 : stats_metadata, & ! In
397 : stats_sfc(:), & ! Inout
398 : precip_fracs%precip_frac(:,:), & ! Out
399 : precip_fracs%precip_frac_1(:,:), & ! Out
400 : precip_fracs%precip_frac_2(:,:), & ! Out
401 0 : precip_frac_tol(:) ) ! Out
402 :
403 0 : if ( err_code == clubb_fatal_error ) then
404 0 : write(fstderr,*) " in setup_pdf_parameters after calling precip_fraction"
405 0 : return
406 : end if
407 :
408 : else
409 :
410 0 : precip_fracs%precip_frac(:,:) = one
411 0 : precip_fracs%precip_frac_1(:,:) = one
412 0 : precip_fracs%precip_frac_2(:,:) = one
413 0 : precip_frac_tol(:) = cloud_frac_min
414 :
415 : end if
416 :
417 : ! Calculate <N_cn> from Nc_in_cloud, whether Nc_in_cloud is predicted or
418 : ! based on a prescribed value, and whether the value is constant or varying
419 : ! over the grid level.
420 0 : if ( .not. l_const_Nc_in_cloud ) then
421 : ! Ncn varies at each vertical level.
422 0 : do k = 1, nz
423 0 : do i = 1, ngrdcol
424 0 : const_Ncnp2_on_Ncnm2(i,k) = hm_metadata%Ncnp2_on_Ncnm2
425 : end do
426 : end do
427 :
428 0 : stdev_const_Ncnp2_on_Ncnm2(:,:) = stdev_L2N( const_Ncnp2_on_Ncnm2(:,:) )
429 :
430 : else ! l_const_Nc_in_cloud
431 : ! Ncn is constant at each vertical level.
432 0 : const_Ncnp2_on_Ncnm2(:,:) = zero
433 0 : stdev_const_Ncnp2_on_Ncnm2(:,:) = zero
434 : end if
435 :
436 0 : const_corr_chi_Ncn_n_cloud(:,:) = corr_array_n_cloud( hm_metadata%iiPDF_Ncn,&
437 0 : hm_metadata%iiPDF_chi )
438 :
439 : call corr_NN2NL( nz, ngrdcol, & ! intent(in)
440 : const_corr_chi_Ncn_n_cloud(:,:), & ! intent(in)
441 : stdev_const_Ncnp2_on_Ncnm2(:,:), const_Ncnp2_on_Ncnm2(:,:), & ! intent(in)
442 0 : const_corr_chi_Ncn(:,:) ) ! intent(out)
443 :
444 : Ncnm(:,:) = Nc_in_cloud_to_Ncnm( &
445 0 : pdf_params%chi_1(:,:), pdf_params%chi_2(:,:), pdf_params%stdev_chi_1(:,:), &
446 0 : pdf_params%stdev_chi_2(:,:), pdf_params%mixt_frac(:,:), Nc_in_cloud(:,:), &
447 0 : pdf_params%cloud_frac_1(:,:), pdf_params%cloud_frac_2(:,:), &
448 0 : const_Ncnp2_on_Ncnm2(:,:), const_corr_chi_Ncn(:,:) )
449 :
450 : ! Boundary Condition.
451 : ! At thermodynamic level k = 1, which is below the model lower boundary, the
452 : ! value of Ncnm does not matter.
453 0 : Ncnm(:,1) = Nc_in_cloud(:,1)
454 :
455 :
456 : ! Calculate the overall variance of a precipitating hydrometeor (hm),
457 : !<hm'^2>.
458 0 : do j = 1, hydromet_dim
459 0 : do k = 1, nz, 1
460 0 : do i = 1, ngrdcol
461 :
462 0 : if ( hydromet(i,k,j) >= hm_metadata%hydromet_tol(j) ) then
463 : ! There is some of the hydrometeor species found at level k.
464 : ! Calculate the variance (overall) of the hydrometeor.
465 : hydrometp2_zt(i,k,j) &
466 0 : = ( ( hm_metadata%hmp2_ip_on_hmm2_ip(j) + one ) / precip_fracs%precip_frac(i,k) - one ) &
467 0 : * hydromet(i,k,j)**2
468 : else
469 0 : hydrometp2_zt(i,k,j) = zero
470 : end if
471 :
472 : end do ! k = 1, nz, 1
473 : end do
474 : end do
475 :
476 : ! Interpolate the overall variance of a hydrometeor, <hm'^2>, to its home on
477 : ! momentum grid levels.
478 0 : do i = 1, hydromet_dim, 1
479 0 : hydrometp2(:,:,i) = zt2zm( nz, ngrdcol, gr, hydrometp2_zt(:,:,i), zero_threshold )
480 0 : hydrometp2(:,nz,i) = zero
481 : end do
482 :
483 :
484 : ! Calculate correlations involving w and Ncn by first calculating the
485 : ! overall covariance of w and Ncn using the down-gradient approximation.
486 0 : if ( l_calc_w_corr ) then
487 :
488 : ! Recalculate wm_zt and wp2_zt. Mean vertical velocity may not be easy to
489 : ! pass into this subroutine from a host model, and wp2_zt needs to have a
490 : ! value consistent with the value it had when the PDF parameters involving w
491 : ! were originally set in subroutine pdf_closure. The variable wp2 has since
492 : ! been advanced, resulting a new wp2_zt. However, the value of wp2 here
493 : ! needs to be consistent with wp2 at the time the PDF parameters were
494 : ! calculated.
495 :
496 : ! Calculate the overall mean of vertical velocity, w, on thermodynamic
497 : ! levels.
498 0 : wm_zt(:,:) = compute_mean_binormal( pdf_params%w_1(:,:), pdf_params%w_2(:,:), &
499 0 : pdf_params%mixt_frac(:,:) )
500 :
501 : ! Calculate the overall variance of vertical velocity on thermodynamic
502 : ! levels.
503 : wp2_zt(:,:) = compute_variance_binormal( wm_zt(:,:), &
504 0 : pdf_params%w_1(:,:), pdf_params%w_2(:,:), &
505 : sigma_w_1(:,:), sigma_w_2(:,:), &
506 0 : pdf_params%mixt_frac(:,:) )
507 :
508 : ! Interpolate the covariances (overall) of w and precipitating
509 : ! hydrometeors to thermodynamic grid levels.
510 0 : do i = 1, hydromet_dim
511 0 : wphydrometp_zt(:,:,i) = zm2zt( nz, ngrdcol, gr, wphydrometp(:,:,i) )
512 : end do
513 :
514 0 : do j = 1, hydromet_dim
515 0 : do k = 1, nz
516 0 : do i = 1, ngrdcol
517 0 : if ( hydromet(i,k,j) < hm_metadata%hydromet_tol(j) ) then
518 0 : wphydrometp_zt(i,k,j) = zero
519 : end if
520 : end do
521 : end do
522 : end do
523 :
524 : ! When the mean value of a precipitating hydrometeor is below tolerance
525 : ! value, it is considered to have a value of 0, and the precipitating
526 : ! hydrometeor does not vary over the grid level. Any covariances
527 : ! involving that precipitating hydrometeor also have values of 0 at that
528 : ! grid level.
529 0 : do j = 1, hydromet_dim
530 0 : do k = 1, nz, 1
531 0 : do i = 1, ngrdcol
532 :
533 : ! Clip the value of covariance <w'hm'> on thermodynamic levels.
534 : call clip_covar_level( clip_wphydrometp, k, l_first_clip_ts, & ! In
535 0 : l_last_clip_ts, dt, wp2_zt(i,k), & ! In
536 0 : hydrometp2_zt(i,k,j), & ! In
537 : l_predict_upwp_vpwp, & ! In
538 : stats_metadata, & ! In
539 : stats_zm(i), & ! Inout
540 : wphydrometp_zt(i,k,j), & ! Inout
541 0 : wphydrometp_chnge(i,k,j) ) ! Out
542 :
543 : end do ! k = 1, nz, 1
544 : end do ! i = 1, hydromet_dim, 1
545 : end do
546 :
547 0 : Kh_zm_c_K_hm(:,:) = -clubb_params(ic_K_hm) * Kh_zm(:,:)
548 :
549 : call calc_xpwp( nz, ngrdcol, gr, &
550 : Kh_zm_c_K_hm, Ncnm, &
551 0 : wpNcnp_zm )
552 :
553 : ! Boundary conditions; We are assuming zero flux at the top.
554 0 : wpNcnp_zm(:,nz) = zero
555 :
556 : ! Interpolate the covariances to thermodynamic grid levels.
557 0 : wpNcnp_zt(:,:) = zm2zt( nz, ngrdcol, gr, wpNcnp_zm(:,:) )
558 :
559 : ! When the mean value of Ncn is below tolerance value, it is considered
560 : ! to have a value of 0, and Ncn does not vary over the grid level. Any
561 : ! covariance involving Ncn also has a value of 0 at that grid level.
562 0 : do k = 1, nz, 1
563 0 : do i = 1, ngrdcol
564 0 : if ( Ncnm(i,k) <= Ncn_tol ) then
565 0 : wpNcnp_zt(i,k) = zero
566 : end if
567 : end do ! k = 1, nz, 1
568 : end do
569 :
570 : end if ! l_calc_w_corr
571 :
572 : ! Unpack CLUBB parameters
573 0 : omicron = clubb_params(iomicron)
574 0 : zeta_vrnce_rat = clubb_params(izeta_vrnce_rat)
575 :
576 : !!! Calculate the means and standard deviations involving PDF variables
577 : !!! -- w, chi, eta, N_cn, and any precipitating hydrometeors (hm in-precip)
578 : !!! -- for each PDF component.
579 : call compute_mean_stdev( nz, ngrdcol, hydromet_dim, &
580 : hydromet(:,:,:), hydrometp2_zt(:,:,:), & ! Intent(in)
581 : hm_metadata, & ! Intent(in)
582 : Ncnm(:,:), pdf_params%mixt_frac(:,:), & ! Intent(in)
583 : precip_fracs%precip_frac(:,:), & ! Intent(in)
584 : precip_fracs%precip_frac_1(:,:), & ! Intent(in)
585 : precip_fracs%precip_frac_2(:,:), & ! Intent(in)
586 : precip_frac_tol(:), & ! Intent(in)
587 : pdf_params%w_1(:,:), & ! Intent(in)
588 : pdf_params%w_2(:,:), & ! Intent(in)
589 : sigma_w_1(:,:), & ! Intent(in)
590 : sigma_w_1(:,:), & ! Intent(in)
591 : pdf_params%chi_1(:,:), & ! Intent(in)
592 : pdf_params%chi_2(:,:), & ! Intent(in)
593 : pdf_params%stdev_chi_1(:,:), & ! Intent(in)
594 : pdf_params%stdev_chi_2(:,:), & ! Intent(in)
595 : pdf_params%stdev_eta_1(:,:), & ! Intent(in)
596 : pdf_params%stdev_eta_2(:,:), & ! Intent(in)
597 : pdf_params%thl_1(:,:), & ! Intent(in)
598 : pdf_params%thl_2(:,:), & ! Intent(in)
599 : pdf_dim, & ! Intent(in)
600 : omicron, zeta_vrnce_rat, & ! Intent(in)
601 : l_const_Nc_in_cloud, & ! Intent(in)
602 : mu_x_1(:,:,:), mu_x_2(:,:,:), & ! Intent(out)
603 : sigma_x_1(:,:,:), sigma_x_2(:,:,:), & ! Intent(out)
604 : hm_1(:,:,:), hm_2(:,:,:), & ! Intent(out)
605 : sigma2_on_mu2_ip_1(:,:,:), & ! Intent(out)
606 0 : sigma2_on_mu2_ip_2(:,:,:) ) ! Intent(out)
607 :
608 : !!! Transform the component means and standard deviations involving
609 : !!! precipitating hydrometeors (hm in-precip) and N_cn -- ln hm and
610 : !!! ln N_cn -- to normal space for each PDF component.
611 : call norm_transform_mean_stdev( nz, ngrdcol, & ! Intent(in)
612 : hydromet_dim, hm_metadata, & ! Intent(in)
613 : hm_1(:,:,:), hm_2(:,:,:), & ! Intent(in)
614 : Ncnm(:,:), pdf_dim, & ! Intent(in)
615 : mu_x_1(:,:,:), mu_x_2(:,:,:), & ! Intent(in)
616 : sigma_x_1(:,:,:), sigma_x_2(:,:,:), & ! Intent(in)
617 : sigma2_on_mu2_ip_1(:,:,:), & ! Intent(in)
618 : sigma2_on_mu2_ip_2(:,:,:), & ! Intent(in)
619 : l_const_Nc_in_cloud, & ! Intent(in)
620 : mu_x_1_n(:,:,:), mu_x_2_n(:,:,:), & ! Intent(out)
621 0 : sigma_x_1_n(:,:,:), sigma_x_2_n(:,:,:) ) ! Intent(out)
622 :
623 :
624 : !!! Calculate the normal space correlations.
625 : !!! The normal space correlations are the same as the true correlations
626 : !!! except when at least one of the variables involved is a precipitating
627 : !!! hydrometeor or Ncn. In these cases, the normal space correlation
628 : !!! involves the natural logarithm of the precipitating hydrometeors, ln hm
629 : !!! (for example, ln r_r and ln N_r), and ln N_cn for each PDF component.
630 0 : if ( l_diagnose_correlations ) then
631 :
632 0 : do k = 2, nz, 1
633 0 : do i = 1, ngrdcol
634 :
635 0 : if ( rcm_pdf(i,k) > rc_tol ) then
636 :
637 : call diagnose_correlations( pdf_dim, hm_metadata%iiPDF_w, & ! In
638 : corr_array_n_cloud, & ! In
639 : l_calc_w_corr, & ! In
640 0 : corr_array_1_n(i,k,:,:) ) ! Out
641 :
642 : call diagnose_correlations( pdf_dim, hm_metadata%iiPDF_w, & ! In
643 : corr_array_n_cloud, & ! In
644 : l_calc_w_corr, & ! In
645 0 : corr_array_2_n(i,k,:,:) ) ! Out
646 :
647 : else
648 :
649 : call diagnose_correlations( pdf_dim, hm_metadata%iiPDF_w, & ! In
650 : corr_array_n_below, & ! In
651 : l_calc_w_corr, & ! In
652 0 : corr_array_1_n(i,k,:,:) ) ! Out
653 :
654 : call diagnose_correlations( pdf_dim, hm_metadata%iiPDF_w, & ! In
655 : corr_array_n_below, & ! In
656 : l_calc_w_corr, & ! In
657 0 : corr_array_2_n(i,k,:,:) ) ! Out
658 :
659 : end if
660 :
661 : end do
662 : end do
663 :
664 0 : do k = 2, nz, 1
665 0 : do i = 1, ngrdcol
666 : call calc_cholesky_corr_mtx_approx ( &
667 : pdf_dim, hm_metadata%iiPDF_w, & ! intent(in)
668 0 : corr_array_1_n(i,k,:,:), & ! intent(in)
669 0 : corr_cholesky_mtx_1(i,k,:,:), corr_array_1_n(i,k,:,:) ) ! intent(out)
670 : end do
671 : end do
672 :
673 0 : do k = 2, nz, 1
674 0 : do i = 1, ngrdcol
675 : call calc_cholesky_corr_mtx_approx( &
676 : pdf_dim, hm_metadata%iiPDF_w, & ! intent(in)
677 0 : corr_array_2_n(i,k,:,:), & ! intent(in)
678 0 : corr_cholesky_mtx_2(i,k,:,:), corr_array_2_n(i,k,:,:) ) ! intent(out)
679 : end do
680 : end do
681 :
682 : else ! if .not. l_diagnose_correlations
683 :
684 : if ( l_fix_w_chi_eta_correlations &
685 0 : .and. .not. l_calc_w_corr ) then
686 :
687 : ! When the flags are set this way, the correlation matrices do not vary with any vertical
688 : ! values, and instead are determined entirely by prescribed values. This results in there
689 : ! being only two unique correlation matrices, one for when the grid box is in cloud
690 : ! and one for when it is not. So instead of setting up correlation matrices for all
691 : ! grid boxes then calculating their Cholesky decomps, we can simply set up two correlation
692 : ! matrices, one for in cloud and one for out cloud, calculate the corresponding
693 : ! Cholesky decompositions, then use the value of rc at each grid box to determine whether
694 : ! we assign the in cloud or out of cloud matrices to that grid box.
695 : call calc_corr_norm_and_cholesky_factor( nz, ngrdcol, iiPDF_type, & ! intent(in)
696 : pdf_dim, hm_metadata, & ! intent(in)
697 : pdf_params%rc_1, pdf_params%rc_2, & ! intent(in)
698 : corr_array_n_cloud, corr_array_n_below, & ! intent(out)
699 : corr_array_1_n, corr_array_2_n, & ! intent(out)
700 0 : corr_cholesky_mtx_1, corr_cholesky_mtx_2 ) ! out
701 :
702 : else
703 :
704 : ! The correlation matrices can vary with vertical values. So we need to set the
705 : ! correlation matrices up for each grid box, then find the Cholesky decomp for each
706 : ! grid box individually. This is very computationally expensive.
707 :
708 : call comp_corr_norm( nz, pdf_dim, ngrdcol, hydromet_dim, hm_metadata, & ! In
709 : wm_zt(:,:), pdf_params%rc_1(:,:), pdf_params%rc_2(:,:), & ! In
710 : pdf_params%mixt_frac(:,:), & ! In
711 : precip_fracs%precip_frac_1(:,:), & ! In
712 : precip_fracs%precip_frac_2(:,:), & ! In
713 : wpNcnp_zt(:,:), wphydrometp_zt(:,:,:), & ! In
714 : mu_x_1(:,:,:), mu_x_2(:,:,:), & ! In
715 : sigma_x_1(:,:,:), sigma_x_2(:,:,:), & ! In
716 : sigma_x_1_n(:,:,:), sigma_x_2_n(:,:,:), & ! In
717 : corr_array_n_cloud, corr_array_n_below, & ! In
718 : pdf_params, & ! In
719 : iiPDF_type, & ! In
720 : l_calc_w_corr, & ! In
721 : l_fix_w_chi_eta_correlations, & ! In
722 0 : corr_array_1_n(:,:,:,:), corr_array_2_n(:,:,:,:) ) ! Out
723 :
724 : ! Compute choleksy factorization for the correlation matrix of 1st PDF component
725 0 : do k = 2, nz, 1
726 0 : do i = 1, ngrdcol
727 0 : call Cholesky_factor( pdf_dim, corr_array_1_n(i,k,:,:), & ! In
728 0 : corr_array_scaling, corr_cholesky_mtx_1(i,k,:,:), & ! Out
729 0 : l_corr_array_scaling ) ! Out
730 : end do
731 : end do
732 :
733 : ! Compute choleksy factorization for the correlation matrix of 2nd PDF component
734 0 : do k = 2, nz, 1
735 0 : do i = 1, ngrdcol
736 0 : call Cholesky_factor( pdf_dim, corr_array_2_n(i,k,:,:), & ! In
737 0 : corr_array_scaling, corr_cholesky_mtx_2(i,k,:,:), & ! Out
738 0 : l_corr_array_scaling ) ! Out
739 : end do
740 : end do
741 :
742 : end if
743 :
744 : end if ! l_diagnose_correlations
745 :
746 : ! hydromet_pdf_params is optional, so if it is not present we simply skip the steps
747 : ! to compute it.
748 0 : if ( present(hydromet_pdf_params) ) then
749 :
750 : !!! Calculate the true correlations for each PDF component.
751 : call denorm_transform_corr( nz, ngrdcol, pdf_dim, hm_metadata, & ! In
752 : sigma_x_1_n(:,:,:), sigma_x_2_n(:,:,:), & ! In
753 : sigma2_on_mu2_ip_1(:,:,:), sigma2_on_mu2_ip_2(:,:,:), & ! In
754 : corr_array_1_n(:,:,:,:), & ! In
755 : corr_array_2_n(:,:,:,:), & ! In
756 0 : corr_array_1(:,:,:,:), corr_array_2(:,:,:,:) ) ! Out
757 :
758 : !!! Pack the PDF parameters
759 : call pack_hydromet_pdf_params( nz, ngrdcol, & ! In
760 : hydromet_dim, hm_metadata, & ! In
761 : hm_1(:,:,:), hm_2(:,:,:), pdf_dim, & ! In
762 : mu_x_1(:,:,:), mu_x_2(:,:,:), & ! In
763 : sigma_x_1(:,:,:), sigma_x_2(:,:,:), & ! In
764 : corr_array_1(:,:,:,:), corr_array_2(:,:,:,:), & ! In
765 0 : hydromet_pdf_params(:,:) ) ! Out
766 :
767 0 : do i = 1, ngrdcol
768 0 : call init_hydromet_pdf_params( hydromet_pdf_params(i,1) ) ! intent(out)
769 : end do
770 :
771 : end if
772 :
773 0 : if ( stats_metadata%l_stats_samp ) then
774 :
775 0 : do j = 1, hydromet_dim
776 0 : if ( stats_metadata%ihmp2_zt(j) > 0 ) then
777 : ! Variance (overall) of the hydrometeor, <hm'^2>.
778 : ! call stat_update_var( stats_metadata%ihmp2_zt(i), hydrometp2_zt(:,i), stats_zt )
779 : ! Switch back to using stat_update_var once the code is generalized
780 : ! to pass in the number of vertical levels.
781 0 : do k = 1, nz, 1
782 0 : do i = 1, ngrdcol
783 0 : call stat_update_var_pt( stats_metadata%ihmp2_zt(j), k, hydrometp2_zt(i,k,j), & ! intent(in)
784 0 : stats_zt(i) ) ! intent(inout)
785 : end do ! k = 1, nz, 1
786 : end do
787 : end if
788 : end do
789 :
790 0 : do i = 1, ngrdcol
791 : call pdf_param_hm_stats( nz, pdf_dim, hydromet_dim, hm_metadata, & ! intent(in)
792 0 : hm_1(i,:,:), hm_2(i,:,:), & ! intent(in)
793 0 : mu_x_1(i,:,:), mu_x_2(i,:,:), & ! intent(in)
794 0 : sigma_x_1(i,:,:), sigma_x_2(i,:,:), & ! intent(in)
795 0 : corr_array_1(i,:,:,:), corr_array_2(i,:,:,:), & ! intent(in)
796 : stats_metadata, & ! intent(in)
797 0 : stats_zt(i) ) ! intent(inout)
798 : end do
799 :
800 : !!! Statistics for normal space PDF parameters involving hydrometeors.
801 0 : do i = 1, ngrdcol
802 : call pdf_param_ln_hm_stats( nz, pdf_dim, hm_metadata, & ! intent(in)
803 0 : mu_x_1_n(i,:,:), mu_x_2_n(i,:,:), & ! intent(in)
804 0 : sigma_x_1_n(i,:,:), sigma_x_2_n(i,:,:), & ! intent(in)
805 0 : corr_array_1_n(i,:,:,:), & ! intent(in)
806 0 : corr_array_2_n(i,:,:,:), & ! intent(in)
807 : stats_metadata, & ! intent(in)
808 0 : stats_zt(i) ) ! intent(inout)
809 : end do
810 :
811 0 : if ( stats_metadata%irtp2_from_chi > 0 ) then
812 :
813 0 : do i = 1, ngrdcol
814 : rtp2_zt_from_chi(i,:) &
815 0 : = compute_rtp2_from_chi( pdf_params%stdev_chi_1(i,:), pdf_params%stdev_chi_2(i,:), &
816 0 : pdf_params%stdev_eta_1(i,:), pdf_params%stdev_eta_2(i,:), &
817 0 : pdf_params%rt_1(i,:), pdf_params%rt_2(j,:), &
818 0 : pdf_params%crt_1(i,:), pdf_params%crt_2(i,:), &
819 0 : pdf_params%mixt_frac(i,:), &
820 0 : corr_array_1_n(i,:,hm_metadata%iiPDF_chi,hm_metadata%iiPDF_eta), &
821 0 : corr_array_2_n(i,:,hm_metadata%iiPDF_chi,hm_metadata%iiPDF_eta) )
822 : end do
823 :
824 0 : rtp2_zm_from_chi = zt2zm( nz, ngrdcol, gr, rtp2_zt_from_chi, zero_threshold )
825 :
826 : ! Switch back to using stat_update_var once the code is generalized
827 : ! to pass in the number of vertical levels.
828 : ! call stat_update_var( irtp2_from_chi, zt2zm( rtp2_zt_from_chi ), &
829 : ! stats_zm )
830 0 : do k = 1, nz, 1
831 0 : do i = 1, ngrdcol
832 0 : call stat_update_var_pt( stats_metadata%irtp2_from_chi, k, rtp2_zm_from_chi(i,k), & !in
833 0 : stats_zm(i) ) ! intent(inout)
834 : end do
835 : end do ! k = 1, nz, 1
836 :
837 : end if
838 :
839 0 : do i = 1, ngrdcol
840 : ! Switch back to using stat_update_var once the code is generalized
841 : ! to pass in the number of vertical levels.
842 0 : if ( stats_metadata%iprecip_frac > 0 ) then
843 : ! Overall precipitation fraction.
844 : ! call stat_update_var( stats_metadata%iprecip_frac, precip_frac, stats_zt )
845 0 : do k = 1, nz, 1
846 0 : call stat_update_var_pt( stats_metadata%iprecip_frac, k, precip_fracs%precip_frac(i,k), & ! intent(in)
847 0 : stats_zt(i) ) ! intent(inout)
848 : end do ! k = 1, nz, 1
849 : end if
850 :
851 0 : if ( stats_metadata%iprecip_frac_1 > 0 ) then
852 : ! Precipitation fraction in PDF component 1.
853 : ! call stat_update_var( stats_metadata%iprecip_frac_1, precip_frac_1, stats_zt )
854 0 : do k = 1, nz, 1
855 0 : call stat_update_var_pt( stats_metadata%iprecip_frac_1, k, precip_fracs%precip_frac_1(i,k), & ! In
856 0 : stats_zt(i) ) ! intent(inout)
857 : end do ! k = 1, nz, 1
858 : end if
859 :
860 0 : if ( stats_metadata%iprecip_frac_2 > 0 ) then
861 : ! Precipitation fraction in PDF component 2.
862 : ! call stat_update_var( stats_metadata%iprecip_frac_2, precip_frac_2, stats_zt )
863 0 : do k = 1, nz, 1
864 0 : call stat_update_var_pt( stats_metadata%iprecip_frac_2, k, precip_fracs%precip_frac_2(i,k), & ! In
865 0 : stats_zt(i) ) ! intent(inout)
866 : end do ! k = 1, nz, 1
867 : end if
868 :
869 0 : if ( stats_metadata%iNcnm > 0 ) then
870 : ! Mean simplified cloud nuclei concentration (overall).
871 : ! call stat_update_var( stats_metadata%iNcnm, Ncnm, stats_zt )
872 0 : do k = 1, nz, 1
873 0 : call stat_update_var_pt( stats_metadata%iNcnm, k, Ncnm(i,k), & ! intent(in)
874 0 : stats_zt(i) ) ! intent(inout)
875 : end do ! k = 1, nz, 1
876 : end if
877 : end do
878 :
879 : end if
880 :
881 : ! Boundary conditions for the output variables at k=1.
882 0 : mu_x_1_n(:,1,:) = zero
883 0 : mu_x_2_n(:,1,:) = zero
884 0 : sigma_x_1_n(:,1,:) = zero
885 0 : sigma_x_2_n(:,1,:) = zero
886 0 : corr_array_1_n(:,1,:,:) = zero
887 0 : corr_array_2_n(:,1,:,:) = zero
888 0 : corr_cholesky_mtx_1(:,1,:,:) = zero
889 0 : corr_cholesky_mtx_2(:,1,:,:) = zero
890 :
891 0 : precip_fracs%precip_frac(:,1) = zero
892 0 : precip_fracs%precip_frac_1(:,1) = zero
893 0 : precip_fracs%precip_frac_2(:,1) = zero
894 :
895 0 : if ( clubb_at_least_debug_level( 2 ) ) then
896 0 : do k = 2, nz
897 0 : do i = 1, ngrdcol
898 0 : call assert_corr_symmetric( corr_array_1_n(i,k,:,:), pdf_dim) ! intent(in)
899 0 : call assert_corr_symmetric( corr_array_2_n(i,k,:,:), pdf_dim) ! intent(in)
900 : end do
901 : end do
902 : end if
903 :
904 : return
905 :
906 0 : end subroutine setup_pdf_parameters
907 :
908 : !=============================================================================
909 0 : subroutine compute_mean_stdev( nz, ngrdcol, hydromet_dim, & ! Intent(in)
910 0 : hydromet, hydrometp2_zt, & ! Intent(in)
911 : hm_metadata, & ! Intent(in)
912 0 : Ncnm, mixt_frac, & ! Intent(in)
913 0 : precip_frac, & ! Intent(in)
914 0 : precip_frac_1, & ! Intent(in)
915 0 : precip_frac_2, & ! Intent(in)
916 0 : precip_frac_tol, & ! Intent(in)
917 0 : w_1, w_2, & ! Intent(in)
918 0 : stdev_w_1, stdev_w_2, & ! Intent(in)
919 0 : chi_1, chi_2, & ! Intent(in)
920 0 : stdev_chi_1, stdev_chi_2, & ! Intent(in)
921 0 : stdev_eta_1, stdev_eta_2, & ! Intent(in)
922 0 : thl_1, thl_2, & ! Intent(in)
923 : pdf_dim, & ! Intent(in)
924 : omicron, zeta_vrnce_rat, & ! Intent(in)
925 : l_const_Nc_in_cloud, & ! Intent(in)
926 0 : mu_x_1, mu_x_2, & ! Intent(out)
927 0 : sigma_x_1, sigma_x_2, & ! Intent(out)
928 0 : hm_1, hm_2, & ! Intent(out)
929 0 : sigma_hm_1_sqd_on_mu_hm_1_sqd, & ! Intent(out)
930 0 : sigma_hm_2_sqd_on_mu_hm_2_sqd ) ! Intent(out)
931 :
932 : ! Description:
933 : ! Calculates the means and standard deviations (for each PDF component) of
934 : ! chi, eta, w, Ncn, and the precipitating hydrometeors. For the
935 : ! precipitating hydrometeors, the component means and standard deviations
936 : ! are in-precip.
937 :
938 : ! References:
939 : !-----------------------------------------------------------------------
940 :
941 : use constants_clubb, only: &
942 : zero ! Constant(s)
943 :
944 : use corr_varnce_module, only: &
945 : hm_metadata_type
946 :
947 : use index_mapping, only: &
948 : pdf2hydromet_idx ! Procedure(s)
949 :
950 : use clubb_precision, only: &
951 : core_rknd ! Variable(s)
952 :
953 : implicit none
954 :
955 : !----------------------- Input Variables -----------------------
956 : integer, intent(in) :: &
957 : nz, & ! Number of model vertical grid levels
958 : pdf_dim, & ! Number of PDF variables
959 : ngrdcol, & ! Number of grid columns
960 : hydromet_dim ! Number of hydrometeor species
961 :
962 : real( kind = core_rknd ), dimension(ngrdcol,nz,hydromet_dim), intent(in) :: &
963 : hydromet, & ! Mean of a hydrometeor (overall) [hm units]
964 : hydrometp2_zt ! Variance of a hydrometeor (overall) [(hm units)^2]
965 :
966 : type (hm_metadata_type), intent(in) :: &
967 : hm_metadata
968 :
969 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
970 : Ncnm, & ! Mean simplified cloud nuclei concentration [num/kg]
971 : mixt_frac, & ! Mixture fraction [-]
972 : precip_frac, & ! Precipitation fraction (overall) [-]
973 : precip_frac_1, & ! Precipitation fraction (1st PDF component) [-]
974 : precip_frac_2 ! Precipitation fraction (2nd PDF component) [-]
975 :
976 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
977 : w_1, & ! Mean of w (1st PDF component) [m/s]
978 : w_2, & ! Mean of w (2nd PDF component) [m/s]
979 : stdev_w_1, & ! Variance of w (1st PDF component) [m^2/s^2]
980 : stdev_w_2, & ! Variance of w (2nd PDF component) [m^2/s^2]
981 : chi_1, & ! Mean of chi (1st PDF component) [kg/kg]
982 : chi_2, & ! Mean of chi (2nd PDF component) [kg/kg]
983 : stdev_chi_1, & ! Standard deviation of chi (1st PDF component) [kg/kg]
984 : stdev_chi_2, & ! Standard deviation of chi (2nd PDF component) [kg/kg]
985 : stdev_eta_1, & ! Standard deviation of eta (1st PDF component) [kg/kg]
986 : stdev_eta_2, & ! Standard deviation of eta (2nd PDF component) [kg/kg]
987 : thl_1, & ! Mean of thl (1st PDF component) [K]
988 : thl_2 ! Mean of thl (2nd PDF component) [K]
989 :
990 : real( kind = core_rknd ), dimension(ngrdcol), intent(in) :: &
991 : precip_frac_tol ! Minimum precip. frac. when hydromet. are present [-]
992 :
993 : real( kind = core_rknd ), intent(in) :: &
994 : omicron, & ! Relative width parameter, omicron = R / Rmax [-]
995 : zeta_vrnce_rat ! Width parameter for sigma_hm_1^2 / mu_hm_1^2 [-]
996 :
997 : logical, intent(in) :: &
998 : l_const_Nc_in_cloud ! Use a constant cloud droplet conc. within cloud (K&K)
999 :
1000 : !----------------------- Output Variables -----------------------
1001 : ! Note: This code assumes to be these arrays in the same order as the
1002 : ! correlation arrays, etc., which is determined by the iiPDF indices.
1003 : ! The order should be as follows: chi, eta, w, Ncn, <precip. hydrometeors>
1004 : ! (indices increasing from left to right).
1005 : real( kind = core_rknd ), dimension(ngrdcol,nz,pdf_dim), intent(out) :: &
1006 : mu_x_1, & ! Mean array of PDF vars. (1st PDF component) [units vary]
1007 : mu_x_2, & ! Mean array of PDF vars. (2nd PDF component) [units vary]
1008 : sigma_x_1, & ! Standard deviation array of PDF vars (comp. 1) [units vary]
1009 : sigma_x_2 ! Standard deviation array of PDF vars (comp. 2) [units vary]
1010 :
1011 : real( kind = core_rknd ), dimension(ngrdcol,nz,hydromet_dim), intent(out) :: &
1012 : hm_1, & ! Mean of a precip. hydrometeor (1st PDF component) [units vary]
1013 : hm_2 ! Mean of a precip. hydrometeor (2nd PDF component) [units vary]
1014 :
1015 : real( kind = core_rknd ), dimension(ngrdcol,nz,pdf_dim), intent(out) :: &
1016 : sigma_hm_1_sqd_on_mu_hm_1_sqd, & ! Ratio sigma_hm_1^2 / mu_hm_1^2 [-]
1017 : sigma_hm_2_sqd_on_mu_hm_2_sqd ! Ratio sigma_hm_2^2 / mu_hm_2^2 [-]
1018 :
1019 : !----------------------- Local Variables -----------------------
1020 : integer :: ivar ! Loop iterator
1021 :
1022 : integer :: hm_idx ! Hydrometeor array index.
1023 :
1024 : !----------------------- Begin Code -----------------------
1025 :
1026 : !!! Initialize output variables.
1027 0 : hm_1(:,:,:) = zero
1028 0 : hm_2(:,:,:) = zero
1029 0 : sigma_hm_1_sqd_on_mu_hm_1_sqd(:,:,:) = zero
1030 0 : sigma_hm_2_sqd_on_mu_hm_2_sqd(:,:,:) = zero
1031 :
1032 :
1033 : !!! Enter the PDF parameters.
1034 :
1035 : !!! Vertical velocity, w.
1036 :
1037 : ! Mean of vertical velocity, w, in PDF component 1.
1038 0 : mu_x_1(:,:,hm_metadata%iiPDF_w) = w_1(:,:)
1039 :
1040 : ! Mean of vertical velocity, w, in PDF component 2.
1041 0 : mu_x_2(:,:,hm_metadata%iiPDF_w) = w_2(:,:)
1042 :
1043 : ! Standard deviation of vertical velocity, w, in PDF component 1.
1044 0 : sigma_x_1(:,:,hm_metadata%iiPDF_w) = stdev_w_1(:,:)
1045 :
1046 : ! Standard deviation of vertical velocity, w, in PDF component 2.
1047 0 : sigma_x_2(:,:,hm_metadata%iiPDF_w) = stdev_w_2(:,:)
1048 :
1049 :
1050 : !!! Extended liquid water mixing ratio, chi.
1051 :
1052 : ! Mean of extended liquid water mixing ratio, chi (old s),
1053 : ! in PDF component 1.
1054 0 : mu_x_1(:,:,hm_metadata%iiPDF_chi) = chi_1(:,:)
1055 :
1056 : ! Mean of extended liquid water mixing ratio, chi (old s),
1057 : ! in PDF component 2.
1058 0 : mu_x_2(:,:,hm_metadata%iiPDF_chi) = chi_2(:,:)
1059 :
1060 : ! Standard deviation of extended liquid water mixing ratio, chi (old s),
1061 : ! in PDF component 1.
1062 0 : sigma_x_1(:,:,hm_metadata%iiPDF_chi) = stdev_chi_1(:,:)
1063 :
1064 : ! Standard deviation of extended liquid water mixing ratio, chi (old s),
1065 : ! in PDF component 2.
1066 0 : sigma_x_2(:,:,hm_metadata%iiPDF_chi) = stdev_chi_2(:,:)
1067 :
1068 :
1069 : !!! Coordinate orthogonal to chi, eta.
1070 :
1071 : ! Mean of eta (old t) in PDF component 1.
1072 : ! Set the component mean values of eta to 0.
1073 : ! The component mean values of eta are not important. They can be set to
1074 : ! anything. They cancel out in the model code. However, the best thing to
1075 : ! do is to set them to 0 and avoid any kind of numerical error.
1076 0 : mu_x_1(:,:,hm_metadata%iiPDF_eta) = zero
1077 :
1078 : ! Mean of eta (old t) in PDF component 2.
1079 : ! Set the component mean values of eta to 0.
1080 : ! The component mean values of eta are not important. They can be set to
1081 : ! anything. They cancel out in the model code. However, the best thing to
1082 : ! do is to set them to 0 and avoid any kind of numerical error.
1083 0 : mu_x_2(:,:,hm_metadata%iiPDF_eta) = zero
1084 :
1085 : ! Standard deviation of eta (old t) in PDF component 1.
1086 0 : sigma_x_1(:,:,hm_metadata%iiPDF_eta) = stdev_eta_1(:,:)
1087 :
1088 : ! Standard deviation of eta (old t) in PDF component 2.
1089 0 : sigma_x_2(:,:,hm_metadata%iiPDF_eta) = stdev_eta_2(:,:)
1090 :
1091 :
1092 : !!! Simplified cloud nuclei concentration, Ncn.
1093 :
1094 : ! Mean of simplified cloud nuclei concentration, Ncn, in PDF component 1.
1095 0 : mu_x_1(:,:,hm_metadata%iiPDF_Ncn) = Ncnm(:,:)
1096 :
1097 : ! Mean of simplified cloud nuclei concentration, Ncn, in PDF component 2.
1098 0 : mu_x_2(:,:,hm_metadata%iiPDF_Ncn) = Ncnm(:,:)
1099 :
1100 : ! Standard deviation of simplified cloud nuclei concentration, Ncn,
1101 : ! in PDF component 1.
1102 0 : if ( .not. l_const_Nc_in_cloud ) then
1103 :
1104 : ! Ncn varies in both PDF components.
1105 0 : sigma_x_1(:,:,hm_metadata%iiPDF_Ncn) = sqrt( hm_metadata%Ncnp2_on_Ncnm2 ) * Ncnm(:,:)
1106 :
1107 0 : sigma_x_2(:,:,hm_metadata%iiPDF_Ncn) = sqrt( hm_metadata%Ncnp2_on_Ncnm2 ) * Ncnm(:,:)
1108 :
1109 : ! Ncn is not an official hydrometeor. However, both the
1110 : ! sigma_hm_1_sqd_on_mu_hm_1_sqd and sigma_hm_2_sqd_on_mu_hm_2_sqd arrays
1111 : ! have size pdf_dim, and both sigma_Ncn_1^2/mu_Ncn_1^2 and
1112 : ! sigma_Ncn_2^2/mu_Ncn_2^2 need to be output as part of these arrays.
1113 0 : sigma_hm_1_sqd_on_mu_hm_1_sqd(:,:,hm_metadata%iiPDF_Ncn) = hm_metadata%Ncnp2_on_Ncnm2
1114 0 : sigma_hm_2_sqd_on_mu_hm_2_sqd(:,:,hm_metadata%iiPDF_Ncn) = hm_metadata%Ncnp2_on_Ncnm2
1115 :
1116 : else ! l_const_Nc_in_cloud
1117 :
1118 : ! Ncn is constant in both PDF components.
1119 0 : sigma_x_1(:,:,hm_metadata%iiPDF_Ncn) = zero
1120 :
1121 0 : sigma_x_2(:,:,hm_metadata%iiPDF_Ncn) = zero
1122 :
1123 : ! Ncn is not an official hydrometeor. However, both the
1124 : ! sigma_hm_1_sqd_on_mu_hm_1_sqd and sigma_hm_2_sqd_on_mu_hm_2_sqd arrays
1125 : ! have size pdf_dim, and both sigma_Ncn_1^2/mu_Ncn_1^2 and
1126 : ! sigma_Ncn_2^2/mu_Ncn_2^2 need to be output as part of these arrays.
1127 0 : sigma_hm_1_sqd_on_mu_hm_1_sqd(:,:,hm_metadata%iiPDF_Ncn) = zero
1128 0 : sigma_hm_2_sqd_on_mu_hm_2_sqd(:,:,hm_metadata%iiPDF_Ncn) = zero
1129 :
1130 : end if ! .not. l_const_Nc_in_cloud
1131 :
1132 :
1133 : !!! Precipitating hydrometeor species.
1134 0 : do ivar = hm_metadata%iiPDF_Ncn+1, pdf_dim, 1
1135 :
1136 0 : hm_idx = pdf2hydromet_idx(ivar,hm_metadata)
1137 :
1138 : call calc_comp_mu_sigma_hm( nz, ngrdcol, & ! In
1139 : hydromet(:,:,hm_idx), hydrometp2_zt(:,:,hm_idx), & ! In
1140 0 : hm_metadata%hmp2_ip_on_hmm2_ip(hm_idx), & ! In
1141 : mixt_frac(:,:), precip_frac(:,:), & ! In
1142 : precip_frac_1(:,:), precip_frac_2(:,:), & ! In
1143 0 : hm_metadata%hydromet_tol(hm_idx), & ! In
1144 : precip_frac_tol(:), & ! In
1145 : thl_1(:,:), thl_2(:,:), & ! In
1146 : omicron, zeta_vrnce_rat, & ! In
1147 : mu_x_1(:,:,ivar), mu_x_2(:,:,ivar), & ! Out
1148 : sigma_x_1(:,:,ivar), sigma_x_2(:,:,ivar), & ! Out
1149 : hm_1(:,:,hm_idx), hm_2(:,:,hm_idx), & ! Out
1150 : sigma_hm_1_sqd_on_mu_hm_1_sqd(:,:,ivar), & ! Out
1151 0 : sigma_hm_2_sqd_on_mu_hm_2_sqd(:,:,ivar) ) ! Out
1152 :
1153 : end do ! ivar = iiPDF_Ncn+1, pdf_dim, 1
1154 :
1155 :
1156 0 : return
1157 :
1158 : end subroutine compute_mean_stdev
1159 :
1160 : !=============================================================================
1161 0 : subroutine calc_corr_norm_and_cholesky_factor( nz, ngrdcol, iiPDF_type, &
1162 : pdf_dim, hm_metadata, &
1163 0 : rc_1, rc_2, &
1164 0 : corr_array_n_cloud, corr_array_n_below, &
1165 0 : corr_array_1_n, corr_array_2_n, &
1166 0 : corr_cholesky_mtx_1, corr_cholesky_mtx_2 )
1167 :
1168 : ! Description: This subroutine computes the correlation arrays and correlation
1169 : ! Cholesky matrices of PDF vars for both components. Here, we assume that
1170 : ! there are only two unique correlation arrays, which allows us to compute
1171 : ! these two unique arrays and their corresponding Cholesky decompositions,
1172 : ! then use rc to determine which one to assign to each grid column and
1173 : ! vertical level. If the correlation arrays vary based on vertically varying
1174 : ! values, then this subroutine is not appropriate.
1175 : !
1176 : ! References:
1177 : ! https://github.com/larson-group/cam/issues/129#issuecomment-816205563
1178 : !-----------------------------------------------------------------------
1179 :
1180 : use constants_clubb, only: &
1181 : rc_tol, &
1182 : zero
1183 :
1184 : use clubb_precision, only: &
1185 : core_rknd ! Variable(s)
1186 :
1187 : use corr_varnce_module, only: &
1188 : hm_metadata_type
1189 :
1190 : use model_flags, only: &
1191 : iiPDF_ADG1, & ! Variable(s)
1192 : iiPDF_ADG2, &
1193 : iiPDF_new_hybrid
1194 :
1195 : use pdf_parameter_module, only: &
1196 : pdf_parameter ! Variable(s)
1197 :
1198 : use matrix_operations, only: &
1199 : Cholesky_factor ! Procedure(s)
1200 :
1201 : implicit none
1202 :
1203 : ! Input Variables
1204 : integer, intent(in) :: &
1205 : nz, & ! Number of vertical levels
1206 : pdf_dim, & ! Number of variables in the corr/mean/stdev arrays
1207 : ngrdcol ! Number of grid columns
1208 :
1209 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
1210 : rc_1, & ! Mean of r_c (1st PDF component) [kg/kg]
1211 : rc_2 ! Mean of r_c (2nd PDF component) [kg/kg]
1212 :
1213 : real( kind = core_rknd ), dimension(pdf_dim,pdf_dim), intent(in) :: &
1214 : corr_array_n_cloud, & ! Prescribed correlation array in cloud [-]
1215 : corr_array_n_below ! Prescribed correlation array below cloud [-]
1216 :
1217 : type (hm_metadata_type), intent(in) :: &
1218 : hm_metadata
1219 :
1220 : integer, intent(in) :: &
1221 : iiPDF_type ! Selected option for the two-component normal (double
1222 : ! Gaussian) PDF type to use for the w, rt, and theta-l (or
1223 : ! w, chi, and eta) portion of CLUBB's multivariate,
1224 : ! two-component PDF.
1225 :
1226 : ! Output Variables
1227 : real( kind = core_rknd ), dimension(ngrdcol,nz,pdf_dim,pdf_dim), intent(out) :: &
1228 : corr_array_1_n, & ! Corr. array (normal space) of PDF vars. (comp. 1) [-]
1229 : corr_array_2_n, & ! Corr. array (normal space) of PDF vars. (comp. 2) [-]
1230 : corr_cholesky_mtx_1, & ! Transposed corr. cholesky matrix, 1st comp. [-]
1231 : corr_cholesky_mtx_2 ! Transposed corr. cholesky matrix, 2nd comp. [-]
1232 :
1233 : ! Local Variables
1234 : real( kind = core_rknd ), dimension(pdf_dim,pdf_dim) :: &
1235 0 : corr_array_cloud, & ! General in cloud corr. matrix
1236 0 : corr_array_below, & ! General out of cloud corr. matrix
1237 0 : corr_cholesky_mtx_cloud, & ! General in cloud Cholesky matrix
1238 0 : corr_cholesky_mtx_below ! General out of cloud Cholesky matrix
1239 :
1240 : logical :: &
1241 : l_corr_array_scaling ! Dummy variable that we need for calling Cholesky_factor
1242 :
1243 : real( kind = core_rknd ), dimension(pdf_dim) :: &
1244 0 : corr_array_scaling ! Dummy variable that we need for calling Cholesky_factor
1245 :
1246 : integer :: ivar, jvar, i, k ! Indices
1247 :
1248 : !-------------------- Begin Code --------------------
1249 :
1250 : ! Initialize correlation arrays with prescribed values
1251 0 : corr_array_cloud(:,:) = corr_array_n_cloud(:,:)
1252 0 : corr_array_below(:,:) = corr_array_n_below(:,:)
1253 :
1254 : ! The ADG1 PDF fixes the correlation of w and rt and the correlation of
1255 : ! w and theta_l to be 0, which means the correlation of w and chi and the
1256 : ! correlation of w and eta must also be 0.
1257 : if ( ( iiPDF_type == iiPDF_ADG1 .or. iiPDF_type == iiPDF_ADG2 &
1258 : .or. iiPDF_type == iiPDF_new_hybrid ) &
1259 0 : .and. l_follow_ADG1_PDF_standards ) then
1260 :
1261 0 : corr_array_cloud(hm_metadata%iiPDF_w,hm_metadata%iiPDF_chi) = zero
1262 0 : corr_array_below(hm_metadata%iiPDF_w,hm_metadata%iiPDF_chi) = zero
1263 :
1264 0 : corr_array_cloud(hm_metadata%iiPDF_w,hm_metadata%iiPDF_eta) = zero
1265 0 : corr_array_below(hm_metadata%iiPDF_w,hm_metadata%iiPDF_eta) = zero
1266 :
1267 : end if
1268 :
1269 : ! Ncn is an inherently in-cloud property, so replace out of cloud correlation values
1270 : ! with in cloud ones.
1271 0 : corr_array_below(hm_metadata%iiPDF_Ncn,hm_metadata%iiPDF_chi) &
1272 0 : = corr_array_cloud(hm_metadata%iiPDF_Ncn,hm_metadata%iiPDF_chi)
1273 :
1274 0 : corr_array_below(hm_metadata%iiPDF_Ncn,hm_metadata%iiPDF_eta) &
1275 0 : = corr_array_cloud(hm_metadata%iiPDF_Ncn,hm_metadata%iiPDF_eta)
1276 :
1277 :
1278 : ! Estimates the correlation of the natural logarithm of a
1279 : ! hydrometeor species and eta using the correlation of chi and eta and the
1280 : ! correlation of chi and the natural logarithm of the hydrometeor. This
1281 : ! facilitates the Cholesky decomposability of the correlation array that will
1282 : ! inevitably be decomposed for SILHS purposes. Without this estimation, we
1283 : ! have found that the resulting correlation matrix cannot be decomposed.
1284 0 : do jvar = hm_metadata%iiPDF_Ncn+1, pdf_dim
1285 :
1286 0 : corr_array_cloud(jvar,hm_metadata%iiPDF_eta) &
1287 : = corr_array_cloud(hm_metadata%iiPDF_eta,hm_metadata%iiPDF_chi) &
1288 0 : * corr_array_cloud(jvar,hm_metadata%iiPDF_chi)
1289 :
1290 : corr_array_below(jvar,hm_metadata%iiPDF_eta) &
1291 : = corr_array_below(hm_metadata%iiPDF_eta,hm_metadata%iiPDF_chi) &
1292 0 : * corr_array_below(jvar,hm_metadata%iiPDF_chi)
1293 :
1294 : end do
1295 :
1296 : ! Calc in cloud Cholesky
1297 : call Cholesky_factor( pdf_dim, corr_array_cloud(:,:), & ! In
1298 : corr_array_scaling, corr_cholesky_mtx_cloud(:,:), & ! Out
1299 0 : l_corr_array_scaling ) ! Out
1300 :
1301 : ! Calc out of cloud Cholesky
1302 : call Cholesky_factor( pdf_dim, corr_array_below(:,:), & ! In
1303 : corr_array_scaling, corr_cholesky_mtx_below(:,:), & ! Out
1304 0 : l_corr_array_scaling ) ! Out
1305 :
1306 : ! Use rc_1 to determine which correlation and Cholesky matrices to assign to 1st PDF
1307 : ! Correlation matrices are symmetric, so we copy ij values to ji indices as well
1308 : ! Cholesky matrices are lower triangular, so we only copy those values
1309 0 : do ivar = 1, pdf_dim
1310 0 : do jvar = ivar, pdf_dim
1311 0 : do k = 1, nz
1312 0 : do i = 1, ngrdcol
1313 :
1314 0 : if ( rc_1(i,k) > rc_tol ) then
1315 : ! Assign in cloud matrices to 1st PDF component
1316 0 : corr_array_1_n(i,k,jvar,ivar) = corr_array_cloud(jvar,ivar)
1317 0 : corr_array_1_n(i,k,ivar,jvar) = corr_array_cloud(jvar,ivar)
1318 0 : corr_cholesky_mtx_1(i,k,jvar,ivar) = corr_cholesky_mtx_cloud(jvar,ivar)
1319 : else
1320 : ! Assign out of cloud matrices to 1st PDF component
1321 0 : corr_array_1_n(i,k,jvar,ivar) = corr_array_below(jvar,ivar)
1322 0 : corr_array_1_n(i,k,ivar,jvar) = corr_array_below(jvar,ivar)
1323 0 : corr_cholesky_mtx_1(i,k,jvar,ivar) = corr_cholesky_mtx_below(jvar,ivar)
1324 : end if
1325 :
1326 : end do
1327 : end do
1328 : end do
1329 : end do
1330 :
1331 : ! Use rc_1 to determine which correlation and Cholesky matrices to assign to 2nd PDF
1332 : ! Correlation matrices are symmetric, so we copy ij values to ji indices as well
1333 : ! Cholesky matrices are lower triangular, so we only copy those values
1334 0 : do ivar = 1, pdf_dim
1335 0 : do jvar = ivar, pdf_dim
1336 0 : do k = 1, nz
1337 0 : do i = 1, ngrdcol
1338 :
1339 0 : if ( rc_2(i,k) > rc_tol ) then
1340 : ! Assign in cloud matrices to 2nd PDF component
1341 0 : corr_array_2_n(i,k,jvar,ivar) = corr_array_cloud(jvar,ivar)
1342 0 : corr_array_2_n(i,k,ivar,jvar) = corr_array_cloud(jvar,ivar)
1343 0 : corr_cholesky_mtx_2(i,k,jvar,ivar) = corr_cholesky_mtx_cloud(jvar,ivar)
1344 : else
1345 : ! Assign out of cloud matrices to 2nd PDF component
1346 0 : corr_array_2_n(i,k,jvar,ivar) = corr_array_below(jvar,ivar)
1347 0 : corr_array_2_n(i,k,ivar,jvar) = corr_array_below(jvar,ivar)
1348 0 : corr_cholesky_mtx_2(i,k,jvar,ivar) = corr_cholesky_mtx_below(jvar,ivar)
1349 : end if
1350 :
1351 : end do
1352 : end do
1353 : end do
1354 : end do
1355 :
1356 : ! Set upper triangular parts of Cholesky matrices to zero
1357 0 : do ivar = 1, pdf_dim
1358 0 : do jvar = 1, ivar-1
1359 0 : corr_cholesky_mtx_1(:,:,jvar,ivar) = zero
1360 0 : corr_cholesky_mtx_2(:,:,jvar,ivar) = zero
1361 : end do
1362 : end do
1363 :
1364 0 : end subroutine calc_corr_norm_and_cholesky_factor
1365 :
1366 : !=============================================================================
1367 0 : subroutine comp_corr_norm( nz, pdf_dim, ngrdcol, hydromet_dim, hm_metadata, &
1368 0 : wm_zt, rc_1, rc_2, &
1369 0 : mixt_frac, &
1370 0 : precip_frac_1, &
1371 0 : precip_frac_2, &
1372 0 : wpNcnp_zt, wphydrometp_zt, &
1373 0 : mu_x_1, mu_x_2, sigma_x_1, sigma_x_2, &
1374 0 : sigma_x_1_n, sigma_x_2_n, &
1375 0 : corr_array_n_cloud, corr_array_n_below, &
1376 : pdf_params, &
1377 : iiPDF_type, &
1378 : l_calc_w_corr, &
1379 : l_fix_w_chi_eta_correlations, &
1380 0 : corr_array_1_n, corr_array_2_n )
1381 :
1382 : ! Description:
1383 :
1384 : ! References:
1385 : !-----------------------------------------------------------------------
1386 :
1387 : use constants_clubb, only: &
1388 : Ncn_tol, &
1389 : one, &
1390 : zero
1391 :
1392 : use diagnose_correlations_module, only: &
1393 : calc_mean, & ! Procedure(s)
1394 : calc_w_corr
1395 :
1396 : use index_mapping, only: &
1397 : pdf2hydromet_idx ! Procedure(s)
1398 :
1399 : use clubb_precision, only: &
1400 : core_rknd ! Variable(s)
1401 :
1402 : use corr_varnce_module, only: &
1403 : hm_metadata_type
1404 :
1405 : use pdf_parameter_module, only: &
1406 : pdf_parameter ! Variable(s)
1407 :
1408 : use matrix_operations, only: &
1409 : mirror_lower_triangular_matrix
1410 :
1411 : implicit none
1412 :
1413 : !------------------------ Input Variables ------------------------
1414 : integer, intent(in) :: &
1415 : nz, & ! Number of vertical levels
1416 : pdf_dim, & ! Number of variables in the corr/mean/stdev arrays
1417 : ngrdcol, & ! Number of grid columns
1418 : hydromet_dim ! Number of hydrometeor species
1419 :
1420 : type (hm_metadata_type), intent(in) :: &
1421 : hm_metadata
1422 :
1423 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
1424 : wm_zt, & ! Mean vertical velocity, <w>, on thermo. levels [m/s]
1425 : rc_1, & ! Mean of r_c (1st PDF component) [kg/kg]
1426 : rc_2, & ! Mean of r_c (2nd PDF component) [kg/kg]
1427 : mixt_frac, & ! Mixture fraction [-]
1428 : precip_frac_1, & ! Precipitation fraction (1st PDF component) [-]
1429 : precip_frac_2, & ! Precipitation fraction (2nd PDF component) [-]
1430 : wpNcnp_zt ! Covariance of w and N_cn on t-levs. [(m/s) num/kg]
1431 :
1432 : real( kind = core_rknd ), dimension(ngrdcol,nz,hydromet_dim), intent(in) :: &
1433 : wphydrometp_zt ! Covariance of w and hm interp. to t-levs. [(m/s)u.v.]
1434 :
1435 : real( kind = core_rknd ), dimension(ngrdcol,nz,pdf_dim), intent(in) :: &
1436 : mu_x_1, & ! Mean of x array (1st PDF component) [units vary]
1437 : mu_x_2, & ! Mean of x array (2nd PDF component) [units vary]
1438 : sigma_x_1, & ! Standard deviation of x array (1st PDF comp.) [un. vary]
1439 : sigma_x_2, & ! Standard deviation of x array (2nd PDF comp.) [un. vary]
1440 : sigma_x_1_n, & ! Std. dev. array (normal space): PDF vars (comp. 1) [u.v.]
1441 : sigma_x_2_n ! Std. dev. array (normal space): PDF vars (comp. 2) [u.v.]
1442 :
1443 : real( kind = core_rknd ), dimension(pdf_dim,pdf_dim), &
1444 : intent(in) :: &
1445 : corr_array_n_cloud, & ! Prescribed correlation array in cloud [-]
1446 : corr_array_n_below ! Prescribed correlation array below cloud [-]
1447 :
1448 : type(pdf_parameter), intent(in) :: &
1449 : pdf_params ! PDF parameters [units vary]
1450 :
1451 : integer, intent(in) :: &
1452 : iiPDF_type ! Selected option for the two-component normal (double
1453 : ! Gaussian) PDF type to use for the w, rt, and theta-l (or
1454 : ! w, chi, and eta) portion of CLUBB's multivariate,
1455 : ! two-component PDF.
1456 :
1457 : logical, intent(in) :: &
1458 : l_calc_w_corr, & ! Calculate the correlations between w and the hydrometeors
1459 : l_fix_w_chi_eta_correlations ! Use a fixed correlation for s and t Mellor(chi/eta)
1460 :
1461 : !------------------------ Output Variables ------------------------
1462 : real( kind = core_rknd ), dimension(ngrdcol,nz,pdf_dim,pdf_dim), intent(out) :: &
1463 : corr_array_1_n, & ! Corr. array (normal space) of PDF vars. (comp. 1) [-]
1464 : corr_array_2_n ! Corr. array (normal space) of PDF vars. (comp. 2) [-]
1465 :
1466 : !------------------------ Local Variables ------------------------
1467 : real( kind = core_rknd ), dimension(ngrdcol,nz,pdf_dim) :: &
1468 0 : corr_w_hm_1_n, & ! Correlation of w and ln hm (1st PDF component) ip [-]
1469 0 : corr_w_hm_2_n ! Correlation of w and ln hm (2nd PDF component) ip [-]
1470 :
1471 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
1472 0 : corr_w_Ncn_1_n, & ! Correlation of w and ln Ncn (1st PDF component) [-]
1473 0 : corr_w_Ncn_2_n ! Correlation of w and ln Ncn (2nd PDF component) [-]
1474 :
1475 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
1476 0 : ones_vector, & ! Vector of 1s
1477 0 : Ncn_tol_in, & ! Tolerance value for Ncn
1478 0 : hydromet_tol_in ! Tolerance value for hydromet
1479 :
1480 : logical :: &
1481 : l_limit_corr_chi_eta ! Flag to limit the correlation of chi and eta [-]
1482 :
1483 : integer :: ivar, jvar, hm_idx, i, k ! Indices
1484 :
1485 : ! Just to avoid typing hm_metadata%iiPDF_x everywhere
1486 : integer :: &
1487 : iiPDF_chi, &
1488 : iiPDF_eta, &
1489 : iiPDF_w, &
1490 : iiPDF_Ncn
1491 :
1492 : !------------------------ Begin Code ------------------------
1493 :
1494 0 : iiPDF_chi = hm_metadata%iiPDF_chi
1495 0 : iiPDF_eta = hm_metadata%iiPDF_eta
1496 0 : iiPDF_w = hm_metadata%iiPDF_w
1497 0 : iiPDF_Ncn = hm_metadata%iiPDF_Ncn
1498 :
1499 :
1500 : !!! Normal space correlations
1501 :
1502 : ! Initialize corr_w_hm_1_n and corr_w_hm_2_n arrays to 0.
1503 0 : corr_w_hm_1_n(:,:,:) = zero
1504 0 : corr_w_hm_2_n(:,:,:) = zero
1505 :
1506 :
1507 : ! Set ones_vector to a vector of 1s.
1508 0 : ones_vector = one
1509 :
1510 : ! Calculate normal space correlations involving w by first calculating total
1511 : ! covariances involving w (<w'Ncn'>, etc.) using the down-gradient
1512 : ! approximation.
1513 0 : if ( l_calc_w_corr ) then
1514 :
1515 0 : Ncn_tol_in = Ncn_tol
1516 :
1517 : ! Calculate the correlation of w and ln Ncn in each PDF component.
1518 : ! The subroutine calc_corr_w_hm_n can be used to do this as long as a
1519 : ! value of 1 is sent in for precip_frac_1 and precip_frac_2.
1520 0 : jvar = iiPDF_Ncn
1521 0 : do i = 1, ngrdcol
1522 0 : call calc_corr_w_hm_n( wm_zt(i,:), wpNcnp_zt(i,:), &
1523 0 : mu_x_1(i,:,iiPDF_w), mu_x_2(i,:,iiPDF_w), &
1524 0 : mu_x_1(i,:,jvar), mu_x_2(i,:,jvar), &
1525 : sigma_x_1(i,:,iiPDF_w), sigma_x_2(i,:,iiPDF_w), &
1526 : sigma_x_1(i,:,jvar), sigma_x_2(i,:,jvar), &
1527 : sigma_x_1_n(i,:,jvar), sigma_x_2_n(i,:,jvar), &
1528 : mixt_frac(i,:), ones_vector(i,:), &
1529 : ones_vector(i,:), Ncn_tol_in(i,:), &
1530 0 : corr_w_Ncn_1_n(i,:), corr_w_Ncn_2_n(i,:) )
1531 : end do
1532 :
1533 : ! Calculate the correlation of w and the natural logarithm of the
1534 : ! hydrometeor for each PDF component and each hydrometeor type.
1535 0 : do jvar = iiPDF_Ncn+1, pdf_dim
1536 :
1537 0 : hm_idx = pdf2hydromet_idx(jvar,hm_metadata)
1538 :
1539 0 : hydromet_tol_in(:,:) = hm_metadata%hydromet_tol(hm_idx)
1540 :
1541 0 : do i = 1, ngrdcol
1542 0 : call calc_corr_w_hm_n( wm_zt(i,:), wphydrometp_zt(i,:,hm_idx), & ! intent(in)
1543 0 : mu_x_1(i,:,iiPDF_w), mu_x_2(i,:,iiPDF_w), & ! intent(in)
1544 0 : mu_x_1(i,:,jvar), mu_x_2(i,:,jvar), & ! intent(in)
1545 : sigma_x_1(i,:,iiPDF_w), sigma_x_2(i,:,iiPDF_w), & ! intent(in)
1546 : sigma_x_1(i,:,jvar), sigma_x_2(i,:,jvar), & ! intent(in)
1547 : sigma_x_1_n(i,:,jvar), sigma_x_2_n(i,:,jvar), & ! intent(in)
1548 : mixt_frac(i,:), precip_frac_1(i,:), & ! intent(in)
1549 : precip_frac_2(i,:), hydromet_tol_in(i,:), & ! intent(in)
1550 0 : corr_w_hm_1_n(i,:,jvar), corr_w_hm_2_n(i,:,jvar) ) ! intent(out)
1551 : end do
1552 :
1553 : end do ! jvar = iiPDF_Ncn+1, pdf_dim
1554 :
1555 : end if ! l_calc_w_corr
1556 :
1557 : ! In order to decompose the normal space correlation matrix,
1558 : ! we must not have a perfect correlation of chi and
1559 : ! eta. Thus, we impose a limitation.
1560 0 : l_limit_corr_chi_eta = .true.
1561 :
1562 :
1563 : ! Initialize the normal space correlation arrays
1564 0 : corr_array_1_n(:,:,:,:) = zero
1565 0 : corr_array_2_n(:,:,:,:) = zero
1566 :
1567 : !!! The corr_arrays are assumed to be lower triangular matrices
1568 : ! Set diagonal elements to 1
1569 0 : do ivar=1, pdf_dim
1570 0 : do k = 1, nz
1571 0 : do i = 1, ngrdcol
1572 0 : corr_array_1_n(i,k,ivar,ivar) = one
1573 0 : corr_array_2_n(i,k,ivar,ivar) = one
1574 : end do
1575 : end do
1576 : end do
1577 :
1578 :
1579 : !!! This code assumes the following order in the prescribed correlation
1580 : !!! arrays (iiPDF indices):
1581 : !!! chi, eta, w, Ncn, <hydrometeors> (indices increasing from left to right)
1582 0 : if ( l_fix_w_chi_eta_correlations ) then
1583 : ! Correlation of chi (old s) and eta (old t)
1584 : call component_corr_chi_eta( nz, ngrdcol, & ! intent(in)
1585 : rc_1(:,:), & ! intent(in)
1586 : rc_2(:,:), & ! intent(in)
1587 0 : corr_array_n_cloud(iiPDF_eta,iiPDF_chi), & ! intent(in)
1588 : corr_array_n_below(iiPDF_eta,iiPDF_chi), & ! intent(in)
1589 : l_limit_corr_chi_eta, & ! intent(in)
1590 : corr_array_1_n(:,:,iiPDF_eta,iiPDF_chi), & ! intent(out)
1591 0 : corr_array_2_n(:,:,iiPDF_eta,iiPDF_chi) ) ! intent(out)
1592 : else
1593 :
1594 : ! Preferred, more accurate version.
1595 0 : do i = 1, ngrdcol
1596 0 : corr_array_1_n(i,:,iiPDF_eta,iiPDF_chi) = pdf_params%corr_chi_eta_1(i,:)
1597 0 : corr_array_2_n(i,:,iiPDF_eta,iiPDF_chi) = pdf_params%corr_chi_eta_2(i,:)
1598 : end do
1599 :
1600 : end if
1601 :
1602 0 : if ( l_fix_w_chi_eta_correlations ) then
1603 : ! Correlation of chi (old s) and w
1604 : call component_corr_w_x( nz, ngrdcol, & ! intent(in)
1605 : rc_1(:,:), & ! intent(in)
1606 : rc_2(:,:), & ! intent(in)
1607 0 : corr_array_n_cloud(iiPDF_w,iiPDF_chi), & ! intent(in)
1608 : corr_array_n_below(iiPDF_w,iiPDF_chi), & ! intent(in)
1609 : iiPDF_type, & ! intent(in)
1610 : corr_array_1_n(:,:,iiPDF_w,iiPDF_chi), & ! intent(out)
1611 0 : corr_array_2_n(:,:,iiPDF_w,iiPDF_chi) ) ! intent(out)
1612 : else
1613 :
1614 : ! Preferred, more accurate version.
1615 0 : do i = 1, ngrdcol
1616 0 : corr_array_1_n(i,:,iiPDF_w,iiPDF_chi) = pdf_params%corr_w_chi_1(i,:)
1617 0 : corr_array_2_n(i,:,iiPDF_w,iiPDF_chi) = pdf_params%corr_w_chi_2(i,:)
1618 : end do
1619 :
1620 : end if
1621 :
1622 :
1623 : ! Correlation of chi (old s) and ln Ncn, corr_array_n_cloud used twice because
1624 : ! Ncn is an inherently in-cloud property.
1625 : call component_corr_x_hm_n_ip( nz, ngrdcol, & ! intent(in)
1626 : rc_1(:,:), & ! intent(in)
1627 : rc_2(:,:), & ! intent(in)
1628 0 : corr_array_n_cloud(iiPDF_Ncn,iiPDF_chi), & ! intent(in)
1629 : corr_array_n_cloud(iiPDF_Ncn,iiPDF_chi), & ! intent(in)
1630 : corr_array_1_n(:,:,iiPDF_Ncn,iiPDF_chi), & ! intent(out)
1631 0 : corr_array_2_n(:,:,iiPDF_Ncn,iiPDF_chi) ) ! intent(out)
1632 :
1633 : ! Correlation of chi (old s) and the natural logarithm of the hydrometeors
1634 0 : do jvar = iiPDF_Ncn+1, pdf_dim
1635 : call component_corr_x_hm_n_ip( nz, ngrdcol, & ! intent(in)
1636 : rc_1(:,:), & ! intent(in)
1637 : rc_2(:,:), & ! intent(in)
1638 0 : corr_array_n_cloud(jvar,iiPDF_chi), & ! intent(in)
1639 : corr_array_n_below(jvar,iiPDF_chi), & ! intent(in)
1640 : corr_array_1_n(:,:,jvar,iiPDF_chi), & ! intent(out)
1641 0 : corr_array_2_n(:,:,jvar,iiPDF_chi) ) ! intent(out)
1642 : end do
1643 :
1644 : ! Correlation of eta (old t) and w
1645 0 : if ( l_fix_w_chi_eta_correlations ) then
1646 : ! Correlation of chi (old s) and w
1647 : call component_corr_w_x( nz, ngrdcol, & ! intent(in)
1648 : rc_1(:,:), & ! intent(in)
1649 : rc_2(:,:), & ! intent(in)
1650 0 : corr_array_n_cloud(iiPDF_w,iiPDF_eta), & ! intent(in)
1651 : corr_array_n_below(iiPDF_w,iiPDF_eta), & ! intent(in)
1652 : iiPDF_type, & ! intent(in)
1653 : corr_array_1_n(:,:,iiPDF_w,iiPDF_eta), & ! intent(out)
1654 0 : corr_array_2_n(:,:,iiPDF_w,iiPDF_eta) ) ! intent(out)
1655 : else
1656 :
1657 : ! Preferred, more accurate version.
1658 0 : do i = 1, ngrdcol
1659 0 : corr_array_1_n(i,:,iiPDF_w,iiPDF_chi) = pdf_params%corr_w_chi_1(i,:)
1660 0 : corr_array_2_n(i,:,iiPDF_w,iiPDF_chi) = pdf_params%corr_w_chi_2(i,:)
1661 : end do
1662 :
1663 : end if
1664 :
1665 : ! Correlation of eta (old t) and ln Ncn, corr_array_n_cloud used twice because
1666 : ! Ncn is an inherently in-cloud property.
1667 : call component_corr_x_hm_n_ip( nz, ngrdcol, & ! intent(in)
1668 : rc_1(:,:), & ! intent(in)
1669 : rc_2(:,:), & ! intent(in)
1670 0 : corr_array_n_cloud(iiPDF_Ncn,iiPDF_eta), & ! intent(in)
1671 : corr_array_n_cloud(iiPDF_Ncn,iiPDF_eta), & ! intent(in)
1672 : corr_array_1_n(:,:,iiPDF_Ncn,iiPDF_eta), & ! intent(out)
1673 0 : corr_array_2_n(:,:,iiPDF_Ncn,iiPDF_eta) ) ! intent(out)
1674 :
1675 :
1676 : ! Correlation of eta (old t) and the natural logarithm of the hydrometeors
1677 0 : do jvar = iiPDF_Ncn+1, pdf_dim
1678 : call component_corr_eta_hm_n_ip( nz, ngrdcol, & ! intent(in)
1679 : corr_array_1_n(:,:,iiPDF_eta,iiPDF_chi), & ! intent(in)
1680 : corr_array_1_n(:,:,jvar,iiPDF_chi), & ! intent(in)
1681 : corr_array_2_n(:,:,iiPDF_eta,iiPDF_chi), & ! intent(in)
1682 : corr_array_2_n(:,:,jvar,iiPDF_chi), & ! intent(in)
1683 : corr_array_1_n(:,:,jvar,iiPDF_eta), & ! intent(out)
1684 0 : corr_array_2_n(:,:,jvar,iiPDF_eta) ) ! intent(out)
1685 : end do
1686 :
1687 :
1688 : ! Correlation of w and ln Ncn
1689 : call component_corr_w_hm_n_ip( nz, ngrdcol, & ! intent(in)
1690 : corr_w_Ncn_1_n(:,:), rc_1(:,:), & ! intent(in)
1691 : corr_w_Ncn_2_n(:,:), rc_2(:,:), & ! intent(in)
1692 0 : corr_array_n_cloud(iiPDF_Ncn,iiPDF_w), & ! intent(in)
1693 : corr_array_n_below(iiPDF_Ncn,iiPDF_w), & ! intent(in)
1694 : l_calc_w_corr, & ! intent(in)
1695 : corr_array_1_n(:,:,iiPDF_Ncn,iiPDF_w), & ! intent(out)
1696 0 : corr_array_2_n(:,:,iiPDF_Ncn,iiPDF_w) ) ! intent(out)
1697 :
1698 : ! Correlation of w and the natural logarithm of the hydrometeors
1699 0 : do jvar = iiPDF_Ncn+1, pdf_dim
1700 : call component_corr_w_hm_n_ip( nz, ngrdcol, & ! intent(in)
1701 : corr_w_hm_1_n(:,:,jvar), rc_1(:,:), & ! intent(in)
1702 : corr_w_hm_2_n(:,:,jvar), rc_2(:,:), & ! intent(in)
1703 0 : corr_array_n_cloud(jvar,iiPDF_w), & ! intent(in)
1704 : corr_array_n_below(jvar,iiPDF_w), & ! intent(in)
1705 : l_calc_w_corr, & ! intent(in)
1706 : corr_array_1_n(:,:,jvar,iiPDF_w), & ! intent(out)
1707 0 : corr_array_2_n(:,:,jvar,iiPDF_w) ) ! intent(out)
1708 : end do
1709 :
1710 : ! Correlation of ln Ncn and the natural logarithm of the hydrometeors
1711 0 : do jvar = iiPDF_Ncn+1, pdf_dim
1712 :
1713 : call component_corr_hmx_hmy_n_ip( nz, ngrdcol, & ! intent(in)
1714 : rc_1(:,:), & ! intent(in)
1715 : rc_2(:,:), & ! intent(in)
1716 0 : corr_array_n_cloud(jvar,iiPDF_Ncn), & ! intent(in)
1717 : corr_array_n_below(jvar,iiPDF_Ncn), & ! intent(in)
1718 : corr_array_1_n(:,:,jvar,iiPDF_Ncn), & ! intent(out)
1719 0 : corr_array_2_n(:,:,jvar,iiPDF_Ncn) ) ! intent(out)
1720 : end do
1721 :
1722 : ! Correlation of the natural logarithm of two hydrometeors
1723 0 : do ivar = iiPDF_Ncn+1, pdf_dim-1
1724 0 : do jvar = ivar+1, pdf_dim
1725 :
1726 :
1727 : call component_corr_hmx_hmy_n_ip( nz, ngrdcol, & ! intent(in)
1728 : rc_1(:,:), & ! intent(in)
1729 : rc_2(:,:), & ! intent(in)
1730 0 : corr_array_n_cloud(jvar,ivar), & ! intent(in)
1731 : corr_array_n_below(jvar,ivar), & ! intent(in)
1732 : corr_array_1_n(:,:,jvar,ivar), & ! intent(out)
1733 0 : corr_array_2_n(:,:,jvar,ivar) ) ! intent(out)
1734 :
1735 : end do ! jvar
1736 : end do ! ivar
1737 :
1738 : ! For ease of use later in the code, we make the correlation arrays
1739 : ! symmetrical
1740 0 : do k = 2, nz, 1
1741 0 : do i = 1, ngrdcol
1742 0 : call mirror_lower_triangular_matrix( pdf_dim, corr_array_1_n(i,k,:,:) )
1743 0 : call mirror_lower_triangular_matrix( pdf_dim, corr_array_2_n(i,k,:,:) )
1744 : end do
1745 : end do
1746 :
1747 0 : return
1748 :
1749 : end subroutine comp_corr_norm
1750 :
1751 : !=============================================================================
1752 0 : subroutine calc_comp_mu_sigma_hm( nz, ngrdcol, & ! In
1753 0 : hmm, hmp2, & ! In
1754 : hmp2_ip_on_hmm2_ip, & ! In
1755 0 : mixt_frac, precip_frac, & ! In
1756 0 : precip_frac_1, precip_frac_2, & ! In
1757 0 : hm_tol, precip_frac_tol, & ! In
1758 0 : mu_thl_1, mu_thl_2, & ! In
1759 : omicron, zeta_vrnce_rat_in, & ! In
1760 0 : mu_hm_1, mu_hm_2, & ! Out
1761 0 : sigma_hm_1, sigma_hm_2, & ! Out
1762 0 : hm_1, hm_2, & ! Out
1763 0 : sigma_hm_1_sqd_on_mu_hm_1_sqd, & ! Out
1764 0 : sigma_hm_2_sqd_on_mu_hm_2_sqd ) ! Out
1765 :
1766 :
1767 : ! Description:
1768 : ! Calculates the in-precipitation mean and in-precipitation standard
1769 : ! deviation in both PDF components for a precipitating hydrometeor.
1770 : !
1771 : ! When precipitation is found in both PDF components (precip_frac_1 > 0 and
1772 : ! precip_frac_2 > 0), the method that solves for in-precip. mean and
1773 : ! in-precip. standard deviation in each PDF component, preserving overall
1774 : ! mean and overall variance, is used. When precipitation fraction is found
1775 : ! in one PDF component but not the other one (precip_frac_1 > 0 and
1776 : ! precip_frac_2 = 0, or precip_frac_1 = 0 and precip_frac_2 > 0), the
1777 : ! calculation of component in-precip. mean and in-precip. standard deviation
1778 : ! is simple. When precipitation is not found in either component
1779 : ! (precip_frac_1 = 0 and precip_frac_2 = 0), there isn't any precipitation
1780 : ! found overall (at that grid level).
1781 : !
1782 : !
1783 : ! DESCRIPTION OF THE METHOD THAT SOLVES FOR TWO IN-PRECIPITATION COMPONENTS
1784 : ! =========================================================================
1785 : !
1786 : ! OVERVIEW
1787 : !
1788 : ! The goal is to calculate the in-precip. mean of the hydrometeor field in
1789 : ! each PDF component (mu_hm_1 and mu_hm_2) in a scenario when there is
1790 : ! precipitation found in both PDF components. The fields provided are the
1791 : ! overall mean of the hydrometeor, <hm>, the overall variance of the
1792 : ! hydrometeor, <hm’^2>, the mixture fraction, a, the overall precipitation
1793 : ! fraction, f_p, and the precipitation fraction in each PDF component
1794 : ! (f_p_1 and f_p_2).
1795 : !
1796 : ! The PDF equation for <hm> is:
1797 : !
1798 : ! <hm> = a * f_p_1 * mu_hm_1 + ( 1- a ) * f_p_2 * mu_hm_2.
1799 : !
1800 : ! Likewise, the PDF equation for <hm’^2> is:
1801 : !
1802 : ! <hm’^2> = a * f_p_1 * ( mu_hm_1^2 + sigma_hm_1^2 )
1803 : ! + ( 1 - a ) * f_p_2 * ( mu_hm_2^2 + sigma_hm_2^2 )
1804 : ! - <hm>^2;
1805 : !
1806 : ! where sigma_hm_1 and sigma_hm_2 are the in-precip. standard deviations of
1807 : ! the hydrometeor field in each PDF component. This can be rewritten as:
1808 : !
1809 : ! <hm’^2>
1810 : ! = a * f_p_1 * ( 1 + sigma_hm_1^2 / mu_hm_1^2 ) * mu_hm_1^2
1811 : ! + ( 1 - a ) * f_p_2 * ( 1 + sigma_hm_2^2 / mu_hm_2^2 ) * mu_hm_2^2
1812 : ! - <hm>^2.
1813 : !
1814 : ! The ratio of sigma_hm_2^2 to mu_hm_2^2 is denoted R:
1815 : !
1816 : ! R = sigma_hm_2^2 / mu_hm_2^2.
1817 : !
1818 : ! In order to allow sigma_hm_1^2 / mu_hm_1^2 to have a different ratio, the
1819 : ! parameter zeta is introduced, such that:
1820 : !
1821 : ! R * ( 1 + zeta ) = sigma_hm_1^2 / mu_hm_1^2;
1822 : !
1823 : ! where zeta > -1. When -1 < zeta < 0, the ratio sigma_hm_2^2 / mu_hm_2^2
1824 : ! grows at the expense of sigma_hm_1^2 / mu_hm_1^2, which narrows. When
1825 : ! zeta = 0, the ratio sigma_hm_1^2 / mu_hm_1^2 is the same as
1826 : ! sigma_hm_2^2 / mu_hm_2^2. When zeta > 0, sigma_hm_1^2 / mu_hm_1^2 grows
1827 : ! at the expense of sigma_hm_2^2 / mu_hm_2^2, which narrows. The component
1828 : ! variances are written as:
1829 : !
1830 : ! sigma_hm_1^2 = R * ( 1 + zeta ) * mu_hm_1^2; and
1831 : ! sigma_hm_2^2 = R * mu_hm_2^2,
1832 : !
1833 : ! and the component standard deviations are simply:
1834 : !
1835 : ! sigma_hm_1 = sqrt( R * ( 1 + zeta ) ) * mu_hm_1; and
1836 : ! sigma_hm_2 = sqrt( R ) * mu_hm_2.
1837 : !
1838 : ! The equation for <hm’^2> can be rewritten as:
1839 : !
1840 : ! <hm’^2> = a * f_p_1 * ( 1 + R * ( 1 + zeta ) ) * mu_hm_1^2
1841 : ! + ( 1 - a ) * f_p_2 * ( 1 + R ) * mu_hm_2^2
1842 : ! - <hm>^2.
1843 : !
1844 : !
1845 : ! HYDROMETEOR IN-PRECIP. VARIANCE:
1846 : ! THE SPREAD OF THE MEANS VS. THE STANDARD DEVIATIONS
1847 : !
1848 : ! Part I: Minimum and Maximum Values for R
1849 : !
1850 : ! The in-precip. variance of the hydrometeor is accounted for through a
1851 : ! combination of the variance of each PDF component and the spread between
1852 : ! the means of each PDF component. At one extreme, the standard deviation
1853 : ! of each component could be set to 0 and the in-prccip. variance could be
1854 : ! accounted for by spreading the PDF component (in-precip.) means far apart.
1855 : ! The value of R in this scenario would be its minimum possible value, which
1856 : ! is 0. At the other extreme, the means of each component could be set
1857 : ! equal to each other and the in-precip. variance could be accounted for
1858 : ! entirely by the PDF component (in-precip.) standard deviations. The value
1859 : ! of R in this scenario would be its maximum possible value, which is Rmax.
1860 : !
1861 : ! In order to calculate the value of Rmax, use the equation set but set
1862 : ! mu_hm_1 = mu_hm_2 and R = Rmax. When this happens:
1863 : !
1864 : ! <hm> = ( a * f_p_1 + ( 1- a ) * f_p_2 ) * mu_hm_i;
1865 : !
1866 : ! and since f_p = a * f_p_1 + ( 1 - a ) * f_p_2:
1867 : !
1868 : ! mu_hm_i = <hm> / f_p = <hm|_ip>;
1869 : !
1870 : ! where <hm|_ip> is the in-precip. mean of the hydrometeor. The equation
1871 : ! for hydrometeor variance in this scenario becomes:
1872 : !
1873 : ! <hm’^2> = <hm|_ip>^2 * ( a * f_p_1 * ( 1 + Rmax * ( 1 + zeta ) )
1874 : ! + ( 1 - a ) * f_p_2 * ( 1 + Rmax ) )
1875 : ! - <hm>^2.
1876 : !
1877 : ! The general equation for the in-precip. variance of a hydrometeor,
1878 : ! <hm|_ip’^2>, is given by:
1879 : !
1880 : ! <hm|_ip’^2> = ( <hm’^2> + <hm>^2 - f_p * <hm|_ip>^2 ) / f_p;
1881 : !
1882 : ! which can be rewritten as:
1883 : !
1884 : ! <hm’^2> + <hm>^2 = f_p * ( <hm|_ip’^2> + <hm|_ip>^2 ).
1885 : !
1886 : ! When the above equation is substituted into the modified PDF equation for
1887 : ! <hm’^2>, Rmax is solved for and the equation is:
1888 : !
1889 : ! Rmax = ( f_p / ( a * f_p_1 * ( 1 + zeta ) + ( 1 - a ) * f_p_2 ) )
1890 : ! * ( <hm|_ip’^2> / <hm|_ip>^2 ).
1891 : !
1892 : ! Here, in the scenario that zeta = 0, both PDF components have the same
1893 : ! mean and same variance, which reduces the in-precip. distribution to an
1894 : ! assumed single lognormal, and the above equation reduces to:
1895 : !
1896 : ! Rmax = <hm|_ip’^2> / <hm|_ip>^2;
1897 : !
1898 : ! which is what is expected in that case.
1899 : !
1900 : !
1901 : ! Part II: Enter omicron
1902 : !
1903 : ! A parameter is used to prescribe the ratio of R to its maximum value,
1904 : ! Rmax. The prescribed parameter is called omicron, where:
1905 : !
1906 : ! R = omicron * Rmax;
1907 : !
1908 : ! where 0 <= omicron <= 1. When omicron = 0, the standard deviation of each
1909 : ! PDF component is 0, and mu_hm_1 is spread as far away from mu_hm_2 as it
1910 : ! needs to be to account for the in-precip. variance. When omicron = 1,
1911 : ! mu_hm_1 is equal to mu_hm_2, and the standard deviations of the PDF
1912 : ! components account for all of the in-precip. variance (and when zeta = 0,
1913 : ! the PDF shape is a single lognormal in-precip.). At intermediate values
1914 : ! of omicron, the means of each PDF component are somewhat spread and each
1915 : ! PDF component has some width. The modified parameters are listed below.
1916 : !
1917 : ! The ratio of sigma_hm_2^2 to mu_hm_2^2 is:
1918 : !
1919 : ! sigma_hm_2^2 / mu_hm_2^2 = omicron * Rmax;
1920 : !
1921 : ! and the ratio of sigma_hm_1^2 / mu_hm_1^2 is:
1922 : !
1923 : ! sigma_hm_1^2 / mu_hm_1^2 = omicron * Rmax * ( 1 + zeta ).
1924 : !
1925 : ! The component variances are written as:
1926 : !
1927 : ! sigma_hm_1^2 = omicron * Rmax * ( 1 + zeta ) * mu_hm_1^2; and
1928 : ! sigma_hm_2^2 = omicron * Rmax * mu_hm_2^2,
1929 : !
1930 : ! and the component standard deviations are simply:
1931 : !
1932 : ! sigma_hm_1 = sqrt( omicron * Rmax * ( 1 + zeta ) ) * mu_hm_1; and
1933 : ! sigma_hm_2 = sqrt( omicron * Rmax ) * mu_hm_2.
1934 : !
1935 : ! The equation set becomes:
1936 : !
1937 : ! [1] <hm> = a * f_p_1 * mu_hm_1 + ( 1- a ) * f_p_2 * mu_hm_2; and
1938 : !
1939 : ! [2] <hm’^2>
1940 : ! = a * f_p_1 * ( 1 + omicron * Rmax * ( 1 + zeta ) ) * mu_hm_1^2
1941 : ! + ( 1 - a ) * f_p_2 * ( 1 + omicron * Rmax ) * mu_hm_2^2
1942 : ! - <hm>^2.
1943 : !
1944 : !
1945 : ! SOLVING THE EQUATION SET FOR MU_HM_1 AND MU_HM_2.
1946 : !
1947 : ! The above system of two equations can be solved for mu_hm_1 and mu_hm_2.
1948 : ! All other quantities in the equation set are known quantities. The
1949 : ! equation for <hm> is rewritten to isolate mu_hm_2:
1950 : !
1951 : ! mu_hm_2 = ( <hm> - a * f_p_1 * mu_hm_1 ) / ( ( 1 - a ) * f_p_2 ).
1952 : !
1953 : ! The above equation is substituted into the equation for <hm’^2>. The
1954 : ! equation for <hm’^2> is rewritten, resulting in:
1955 : !
1956 : ! [ a * f_p_1 * ( 1 + omicron * Rmax * ( 1 + zeta ) )
1957 : ! + a^2 * f_p_1^2 * ( 1 + omicron * Rmax ) / ( ( 1 - a ) * f_p_2 ) ]
1958 : ! * mu_hm_1^2
1959 : ! + [ - 2 * <hm> * a * f_p_1 * ( 1 + omicron * Rmax )
1960 : ! / ( ( 1 - a ) * f_p_2 ) ] * mu_hm_1
1961 : ! + [ - ( <hm’^2>
1962 : ! + ( 1 - ( 1 + omicron * Rmax ) / ( ( 1 - a ) * f_p_2 ) )
1963 : ! * <hm>^2 ) ]
1964 : ! = 0.
1965 : !
1966 : ! This equation is of the form:
1967 : !
1968 : ! A * mu_hm_1^2 + B * mu_hm_1 + C = 0;
1969 : !
1970 : ! so the solution for mu_hm_1 is:
1971 : !
1972 : ! mu_hm_1 = ( -B +/- sqrt( B^2 - 4*A*C ) ) / (2*A);
1973 : !
1974 : ! where:
1975 : !
1976 : ! A = a * f_p_1 * ( 1 + omicron * Rmax * ( 1 + zeta ) )
1977 : ! + a^2 * f_p_1^2 * ( 1 + omicron * Rmax ) / ( ( 1 - a ) * f_p_2 );
1978 : !
1979 : ! B = - 2 * <hm> * a * f_p_1 * ( 1 + omicron * Rmax )
1980 : ! / ( ( 1 - a ) * f_p_2 );
1981 : !
1982 : ! and
1983 : !
1984 : ! C = - ( <hm’^2>
1985 : ! + ( 1 - ( 1 + omicron * Rmax ) / ( ( 1 - a ) * f_p_2 ) )
1986 : ! * <hm>^2 ).
1987 : !
1988 : ! The signs of the coefficients:
1989 : !
1990 : ! 1) coefficient A is always positive,
1991 : ! 2) coefficient B is always negative (this means that -B is always
1992 : ! positive), and
1993 : ! 3) coefficient C can be positive, negative, or zero.
1994 : !
1995 : ! Since ( 1 - ( 1 + omicron * Rmax ) / ( ( 1 - a ) * f_p_2 ) ) * <hm>^2 is
1996 : ! always negative and <hm’^2> is always positive, the sign of coefficient C
1997 : ! depends on which term is greater in magnitude.
1998 : !
1999 : ! When <hm’^2> is greater, the sign of coefficient C is negative. This
2000 : ! means that -4*A*C is positive, which in turn means that
2001 : ! sqrt( B^2 - 4*A*C ) is greater in magnitude than -B. If the subtraction
2002 : ! option of the +/- were to be chosen, the value of mu_hm_1 would be
2003 : ! negative in this scenerio. So the natural thing to do would be to always
2004 : ! choose the addition option. However, this method requires that mu_hm_1
2005 : ! equals mu_hm_2 when omicron = 1. When zeta >= 0, this happens when the
2006 : ! addition option is chosen, but not when the subtraction option is chosen.
2007 : ! However, when zeta < 0, this happens when the subtraction option is
2008 : ! chosen, but not when the addition option is chosen. So, the equation for
2009 : ! mu_hm_1 becomes:
2010 : !
2011 : ! mu_hm_1 = ( -B + sqrt( B^2 - 4*A*C ) ) / (2*A); when zeta >= 0; and
2012 : ! mu_hm_1 = ( -B - sqrt( B^2 - 4*A*C ) ) / (2*A); when zeta < 0.
2013 : !
2014 : ! Once this is set, of course:
2015 : !
2016 : ! mu_hm_2 = ( <hm> - a * f_p_1 * mu_hm_1 ) / ( ( 1 - a ) * f_p_2 ).
2017 : !
2018 : ! The system has been solved and the in-precip. PDF component means have
2019 : ! been found!
2020 : !
2021 : !
2022 : ! NOTES
2023 : !
2024 : ! Note 1:
2025 : !
2026 : ! The term B^2 - 4*A*C has been analyzed, and mathematically:
2027 : !
2028 : ! B^2 - 4*A*C >= 0
2029 : !
2030 : ! always holds true. Additionally, the minimum value:
2031 : !
2032 : ! B^2 - 4*A*C = 0,
2033 : !
2034 : ! can only occur when omicron = 1 and zeta = 0 (or alternatively to
2035 : ! zeta = 0, Rmax = 0, but this only occurs when <hm|_ip'^2> / <hm|_ip>^2 has
2036 : ! a value of 0).
2037 : !
2038 : ! Numerically, when omicron = 1 and zeta = 0, B^2 - 4*A*C can produce very
2039 : ! small (on the order of epsilon) negative values. This is due to numerical
2040 : ! round off error. When this happens, the erroneous small, negative value
2041 : ! of B^2 - 4*A*C is simply reset to the value it's supposed to have, which
2042 : ! is 0.
2043 : !
2044 : !
2045 : ! Note 2:
2046 : !
2047 : ! As the value of <hm|_ip'^2> / <hm|_ip>^2 increases and as the value of
2048 : ! omicron decreases (narrowing the in-precip standard deviations and
2049 : ! increasing the spread between the in-precip means), a situtation arises
2050 : ! where the value of one of the component means will become negative. This
2051 : ! is because there is a limit to the amount of in-precip variance that can
2052 : ! be represented by this kind of distribution. In order to prevent
2053 : ! out-of-bounds values of mu_hm_1 or mu_hm_2, lower limits will be
2054 : ! declared, called mu_hm_1_min and mu_hm_2_min. The value of the
2055 : ! hydrometeor in-precip. component mean will be limited from going any
2056 : ! smaller (or negative) at this value. From there, the value of the other
2057 : ! hydrometeor in-precip. component mean is easy to calculate. Then, both
2058 : ! values will be entered into the calculation of hydrometeor variance, which
2059 : ! will be rewritten to solve for R. Then, both the hydrometeor mean and
2060 : ! hydrometeor variance will be preserved with a valid distribution.
2061 : !
2062 : ! In this emergency scenario, the value of R is:
2063 : !
2064 : ! R = ( <hm'^2> + <hm>^2 - a * f_p_1 * mu_hm_1^2
2065 : ! - ( 1 - a ) * f_p_2 * mu_hm_2^2 )
2066 : ! / ( a * f_p_1 * ( 1 + zeta ) * mu_hm_1^2
2067 : ! + ( 1 - a ) * f_p_2 * mu_hm_2^2 ).
2068 : !
2069 : ! The minimum values of the in-precip. component means are bounded by:
2070 : !
2071 : ! mu_hm_1_min >= hm_tol / f_p_1; and
2072 : ! mu_hm_2_min >= hm_tol / f_p_2.
2073 : !
2074 : ! These are set this way because hm_1 ( = mu_hm_1 * f_p_1 ) and
2075 : ! hm_2 ( = mu_hm_2 * f_p_2 ) need to have values of at least hm_tol when
2076 : ! precipitation is found in both PDF components.
2077 : !
2078 : ! However, an in-precip. component mean value of hm_tol / f_p_1 or
2079 : ! hm_tol / f_p_2 often produces a distribution where one component centers
2080 : ! around values that are too small to be a good match with data taken from
2081 : ! Large Eddy Simulations (LES). It is desirable to increase the minimum
2082 : ! threshold of mu_hm_1 and mu_hm_2.
2083 : !
2084 : ! As the minimum threshold increases, the value of the in-precip. component
2085 : ! mean that is from the component that is not being set to the minimum
2086 : ! threshold decreases. If the minimum threshold were to be boosted as high
2087 : ! as <hm> / f_p (in most cases, <hm> / f_p >> hm_tol / f_p_i), both
2088 : ! components would have a value of <hm> / f_p. The minimum threshold should
2089 : ! not be set this high.
2090 : !
2091 : ! Additionally, the minimum threshold for one in-precip. component mean
2092 : ! cannot be set so high as to drive the other in-precip. component mean
2093 : ! below hm_tol / f_p_i. (This doesn't come into play unless <hm> is close
2094 : ! to hm_tol.) The upper limit for the in-precip. mean values are:
2095 : !
2096 : ! mu_hm_1|_(upper. lim.) = ( <hm> - ( 1 - a ) * f_p_2 * ( hm_tol / f_p_2 ) )
2097 : ! / ( a * f_p_1 ); and
2098 : !
2099 : ! mu_hm_2|_(upper. lim.) = ( <hm> - a * f_p_1 * ( hm_tol / f_p_1 ) )
2100 : ! / ( ( 1 - a ) * f_p_2 );
2101 : !
2102 : ! which reduces to:
2103 : !
2104 : ! mu_hm_1|_(upper. lim.) = ( <hm> - ( 1 - a ) * hm_tol ) / ( a * f_p_1 );
2105 : ! and
2106 : ! mu_hm_2|_(upper. lim.) = ( <hm> - a * hm_tol ) / ( ( 1 - a ) * f_p_2 ).
2107 : !
2108 : ! An appropriate minimum value for mu_hm_1 can be set by:
2109 : !
2110 : ! mu_hm_1_min = | min( hm_tol / f_p_1
2111 : ! | + mu_hm_min_coef * ( <hm> / f_p - hm_tol / f_p_1 ),
2112 : ! | ( <hm> - ( 1 - a ) * hm_tol ) / ( a * f_p_1 ) );
2113 : ! | where <hm> / f_p > hm_tol / f_p_1;
2114 : ! | hm_tol / f_p_1;
2115 : ! | where <hm> / f_p <= hm_tol / f_p_1;
2116 : !
2117 : ! and similarly for mu_hm_2:
2118 : !
2119 : ! mu_hm_2_min = | min( hm_tol / f_p_2
2120 : ! | + mu_hm_min_coef * ( <hm> / f_p - hm_tol / f_p_2 ),
2121 : ! | ( <hm> - a * hm_tol ) / ( ( 1 - a ) * f_p_2 ) );
2122 : ! | where <hm> / f_p > hm_tol / f_p_2;
2123 : ! | hm_tol / f_p_2;
2124 : ! | where <hm> / f_p <= hm_tol / f_p_2;
2125 : !
2126 : ! where mu_hm_min_coef is a coefficient that has a value
2127 : ! 0 <= mu_hm_min_coef < 1. When the value of mu_hm_min_coef is 0,
2128 : ! mu_hm_1_min reverts to hm_tol / f_p_1 and mu_hm_2_min reverts to
2129 : ! hm_tol / f_p_2. An appropriate value for mu_hm_min_coef should be small,
2130 : ! such as 0.01 - 0.05.
2131 : !
2132 : !
2133 : ! Note 3:
2134 : !
2135 : ! When the value of zeta >= 0, the value of mu_hm_1 tends to be larger than
2136 : ! the value of mu_hm_2. Likewise when the value of zeta < 0, the value of
2137 : ! mu_hm_2 tends to be larger than the value of mu_hm_1. Since most cloud
2138 : ! water and cloud fraction tends to be found in PDF component 1, it is
2139 : ! advantageous to have the larger in-precip. component mean of the
2140 : ! hydrometeor also found in PDF component 1. The recommended value of zeta
2141 : ! is a value greater than or equal to 0.
2142 : !
2143 : !
2144 : ! Update:
2145 : !
2146 : ! In order to better represent the increase in <th_l'^2> near the ground
2147 : ! from the evaporation of rain, the code is modified to tend towards a
2148 : ! negative correlation of th_l and hm. When mu_thl_1 <= mu_thl_2,
2149 : ! mu_hm_1 >= mu_hm_2; otherwise, when mu_thl_1 > mu_thl_2,
2150 : ! mu_hm_1 <= mu_hm_2.
2151 : !
2152 : ! In the original derivation, mu_hm_1 >= mu_hm_2 when zeta >= 0, and
2153 : ! mu_hm_1 < mu_hm_2 when zeta < 0, where zeta is a tunable or adjustable
2154 : ! parameter. In order to allow the relationship of mu_hm_1 to mu_hm_2 to
2155 : ! depend on the relationship of mu_thl_1 to mu_thl_2, the value of zeta must
2156 : ! also depend on the relationship of mu_thl_1 to mu_thl_2.
2157 : !
2158 : ! When mu_thl_1 <= mu_thl_2:
2159 : !
2160 : ! The relationship of mu_hm_1 to mu_hm_2 is mu_hm_1 >= mu_hm_2, so
2161 : ! zeta >= 0. The tunable value of zeta is referred to as zeta_in. When
2162 : ! zeta_in is already greater than 0 (meaning sigma_hm_1^2 / mu_hm_1^2 is
2163 : ! greater than sigma_hm_2^2 / mu_hm_2^2), zeta is simply set to zeta_in. In
2164 : ! other words, the component with the greater mean value of the hydrometeor
2165 : ! also has the greater value of component variance to the square of the
2166 : ! component mean. However, when zeta_in is less than 0, zeta must be
2167 : ! adjusted to be greater than 0. The following equation is used to set the
2168 : ! value of zeta:
2169 : !
2170 : ! | zeta_in; when zeta_in >= 0
2171 : ! zeta = |
2172 : ! | ( 1 / ( 1 + zeta_in ) ) - 1; when zeta_in < 0.
2173 : !
2174 : ! Previously, when zeta_in < 0, mu_hm_1 < mu_hm_2, and
2175 : ! sigma_hm_1^2 / mu_hm_1^2 < sigma_hm_2^2 / mu_hm_2^2. Now that zeta has to
2176 : ! be greater than 0, sigma_hm_1^2 / mu_hm_1^2 > sigma_hm_2^2 / mu_hm_2^2.
2177 : ! The ratio of the greater variance-over-mean-squared to the smaller
2178 : ! variance-over-mean-squared remains the same when using the equation for
2179 : ! zeta listed above.
2180 : !
2181 : ! When mu_thl_1 > mu_thl_2:
2182 : !
2183 : ! The relationship of mu_hm_1 to mu_hm_2 is mu_hm_1 <= mu_hm_2, so
2184 : ! zeta <= 0. When zeta_in is already less than 0 (meaning
2185 : ! sigma_hm_1^2 / mu_hm_1^2 is less than sigma_hm_2^2 / mu_hm_2^2), zeta is
2186 : ! simply set to zeta_in. In other words, the component with the greater
2187 : ! mean value of the hydrometeor also has the greater value of component
2188 : ! variance to the square of the component mean. However, when zeta_in is
2189 : ! greater than 0, zeta must be adjusted to be less than 0. The following
2190 : ! equation is used to set the value of zeta:
2191 : !
2192 : ! | zeta_in; when zeta_in <= 0
2193 : ! zeta = |
2194 : ! | ( 1 / ( 1 + zeta_in ) ) - 1; when zeta_in > 0.
2195 : !
2196 : ! Previously, when zeta_in > 0, mu_hm_1 > mu_hm_2, and
2197 : ! sigma_hm_1^2 / mu_hm_1^2 > sigma_hm_2^2 / mu_hm_2^2. Now that zeta has to
2198 : ! be less than 0, sigma_hm_1^2 / mu_hm_1^2 < sigma_hm_2^2 / mu_hm_2^2.
2199 : ! The ratio of the greater variance-over-mean-squared to the smaller
2200 : ! variance-over-mean-squared remains the same when using the equation for
2201 : ! zeta listed above.
2202 : !
2203 : !
2204 : ! Brian Griffin; February 2015.
2205 :
2206 : ! References:
2207 : !-----------------------------------------------------------------------
2208 :
2209 : use constants_clubb, only: &
2210 : four, & ! Constant(s)
2211 : two, &
2212 : one, &
2213 : zero
2214 :
2215 : use clubb_precision, only: &
2216 : core_rknd ! Variable(s)
2217 :
2218 : implicit none
2219 :
2220 : ! Input Variables
2221 : integer, intent(in) :: &
2222 : nz, & ! Number of model vertical grid levels
2223 : ngrdcol ! Number of grid columns
2224 :
2225 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
2226 : hmm, & ! Hydrometeor mean (overall), <hm> [hm units]
2227 : hmp2, & ! Hydrometeor variance (overall), <hm'^2> [hm units^2]
2228 : mixt_frac, & ! Mixture fraction [-]
2229 : precip_frac, & ! Precipitation fraction (overall) [-]
2230 : precip_frac_1, & ! Precipitation fraction (1st PDF component) [-]
2231 : precip_frac_2, & ! Precipitation fraction (2nd PDF component) [-]
2232 : mu_thl_1, & ! Mean of th_l (1st PDF component) [K]
2233 : mu_thl_2 ! Mean of th_l (2nd PDF component) [K]
2234 :
2235 : real( kind = core_rknd ), intent(in) :: &
2236 : hmp2_ip_on_hmm2_ip, & ! Ratio <hm|_ip'^2> / <hm|_ip>^2 [-]
2237 : hm_tol ! Tolerance value of hydrometeor [hm units]
2238 :
2239 : real( kind = core_rknd ), dimension(ngrdcol), intent(in) :: &
2240 : precip_frac_tol ! Min. precip. frac. when hydromet. are present [-]
2241 :
2242 : real( kind = core_rknd ), intent(in) :: &
2243 : omicron, & ! Relative width parameter, omicron = R / Rmax [-]
2244 : zeta_vrnce_rat_in ! Width parameter for sigma_hm_1^2 / mu_hm_1^2 [-]
2245 :
2246 : ! Output Variables
2247 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
2248 : mu_hm_1, & ! Mean of hm (1st PDF component) in-precip (ip) [hm un]
2249 : mu_hm_2, & ! Mean of hm (2nd PDF component) ip [hm un]
2250 : sigma_hm_1, & ! Standard deviation of hm (1st PDF component) ip [hm un]
2251 : sigma_hm_2, & ! Standard deviation of hm (2nd PDF component) ip [hm un]
2252 : hm_1, & ! Mean of hm (1st PDF component) [hm un]
2253 : hm_2 ! Mean of hm (2nd PDF component) [hm un]
2254 :
2255 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
2256 : sigma_hm_1_sqd_on_mu_hm_1_sqd, & ! Ratio sigma_hm_1**2 / mu_hm_1**2 [-]
2257 : sigma_hm_2_sqd_on_mu_hm_2_sqd ! Ratio sigma_hm_2**2 / mu_hm_2**2 [-]
2258 :
2259 : ! Local Variables
2260 : real( kind = core_rknd ) :: &
2261 : Rmax, & ! Maximum possible value of ratio R [-]
2262 : coef_A, & ! Coefficient A in A*mu_hm_1^2 + B*mu_hm_1 + C = 0 [-]
2263 : coef_B, & ! Coefficient B in A*mu_hm_1^2 + B*mu_hm_1 + C = 0 [hm un]
2264 : coef_C, & ! Coefficient C in A*mu_hm_1^2 + B*mu_hm_1 + C = 0 [hm un^2]
2265 : Bsqd_m_4AC ! Value B^2 - 4*A*C in quadratic eqn. for mu_hm_1 [hm un^2]
2266 :
2267 : real( kind = core_rknd ) :: &
2268 : zeta_vrnce_rat ! Width parameter for sigma_hm_1^2 / mu_hm_1^2 [-]
2269 :
2270 : real( kind = core_rknd ) :: &
2271 : mu_hm_1_min, & ! Minimum value of mu_hm_1 (precip. in both comps.) [hm un]
2272 : mu_hm_2_min ! Minimum value of mu_hm_2 (precip. in both comps.) [hm un]
2273 :
2274 : real( kind = core_rknd ), parameter :: &
2275 : mu_hm_min_coef = 0.01_core_rknd ! Coef. for mu_hm_1_min and mu_hm_2_min
2276 :
2277 : integer :: i, k ! Loop iterators
2278 :
2279 0 : do k = 1, nz
2280 0 : do i = 1, ngrdcol
2281 :
2282 0 : if ( hmm(i,k) >= hm_tol &
2283 : .and. precip_frac_1(i,k) >= precip_frac_tol(i) &
2284 0 : .and. precip_frac_2(i,k) >= precip_frac_tol(i) ) then
2285 :
2286 : ! Adjust the value of zeta based on the relationship of mu_thl_1 to
2287 : ! mu_thl_2.
2288 0 : if ( mu_thl_1(i,k) <= mu_thl_2(i,k) ) then
2289 0 : if ( zeta_vrnce_rat_in >= zero ) then
2290 : zeta_vrnce_rat = zeta_vrnce_rat_in
2291 : else ! zeta_vrnce_rat_in < 0
2292 0 : zeta_vrnce_rat = ( one / ( one + zeta_vrnce_rat_in ) ) - one
2293 : end if ! zeta_vrnce_rat_in >= 0
2294 : else ! mu_thl_1 > mu_thl_2
2295 0 : if ( zeta_vrnce_rat_in <= zero ) then
2296 : zeta_vrnce_rat = zeta_vrnce_rat_in
2297 : else ! zeta_vrnce_rat_in > 0
2298 0 : zeta_vrnce_rat = ( one / ( one + zeta_vrnce_rat_in ) ) - one
2299 : end if ! zeta_vrnce_rat_in <= 0
2300 : end if ! mu_thl_1 <= mu_thl_2
2301 :
2302 : ! Calculate the value of Rmax.
2303 : ! Rmax = ( f_p / ( a * f_p_1 * ( 1 + zeta ) + ( 1 - a ) * f_p_2 ) )
2304 : ! * ( <hm|_ip’^2> / <hm|_ip>^2 ).
2305 : ! The parameter zeta is written in the code as zeta_vrnce_rat.
2306 : Rmax = ( precip_frac(i,k) &
2307 : / ( mixt_frac(i,k) * precip_frac_1(i,k) * ( one + zeta_vrnce_rat ) &
2308 : + ( one - mixt_frac(i,k) ) * precip_frac_2(i,k) ) ) &
2309 0 : * hmp2_ip_on_hmm2_ip
2310 :
2311 : ! Calculate the value of coefficient A.
2312 : ! A = a * f_p_1 * ( 1 + omicron * Rmax * ( 1 + zeta ) )
2313 : ! + a^2 * f_p_1^2 * ( 1 + omicron * Rmax ) / ( ( 1 - a ) * f_p_2 ).
2314 : coef_A = mixt_frac(i,k) * precip_frac_1(i,k) &
2315 : * ( one + omicron * Rmax * ( one + zeta_vrnce_rat ) ) &
2316 : + mixt_frac(i,k)**2 * precip_frac_1(i,k)**2 &
2317 : * ( one + omicron * Rmax ) &
2318 0 : / ( ( one - mixt_frac(i,k) ) * precip_frac_2(i,k) )
2319 :
2320 : ! Calculate the value of coefficient B.
2321 : ! B = - 2 * <hm> * a * f_p_1 * ( 1 + omicron * Rmax )
2322 : ! / ( ( 1 - a ) * f_p_2 ).
2323 : coef_B = -two * hmm(i,k) * mixt_frac(i,k) * precip_frac_1(i,k) &
2324 : * ( one + omicron * Rmax ) &
2325 0 : / ( ( one - mixt_frac(i,k) ) * precip_frac_2(i,k) )
2326 :
2327 : ! Calculate the value of coefficient C.
2328 : ! C = - ( <hm’^2>
2329 : ! + ( 1 - ( 1 + omicron * Rmax ) / ( ( 1 - a ) * f_p_2 ) )
2330 : ! * <hm>^2 ).
2331 : coef_C = - ( hmp2(i,k) + ( one &
2332 : - ( one + omicron * Rmax ) &
2333 : / ( ( one - mixt_frac(i,k) ) * precip_frac_2(i,k) ) &
2334 0 : ) * hmm(i,k)**2 )
2335 :
2336 : ! Calculate value of B^2 - 4*A*C.
2337 0 : Bsqd_m_4AC = coef_B**2 - four * coef_A * coef_C
2338 :
2339 : ! Mathematically, the value of B^2 - 4*A*C cannot be less than 0.
2340 : ! Numerically, this can happen when numerical round off error causes an
2341 : ! epsilon-sized negative value. When this happens, reset the value of
2342 : ! B^2 - 4*A*C to 0.
2343 0 : if ( Bsqd_m_4AC < zero ) then
2344 0 : Bsqd_m_4AC = zero
2345 : end if
2346 :
2347 : ! Calculate the mean (in-precip.) of the hydrometeor in the 1st PDF
2348 : ! component.
2349 0 : if ( mu_thl_1(i,k) <= mu_thl_2(i,k) ) then
2350 0 : mu_hm_1(i,k) = ( -coef_B + sqrt( Bsqd_m_4AC ) ) / ( two * coef_A )
2351 : else ! mu_thl_1 > mu_thl_2
2352 0 : mu_hm_1(i,k) = ( -coef_B - sqrt( Bsqd_m_4AC ) ) / ( two * coef_A )
2353 : end if ! mu_thl_1 <= mu_thl_2
2354 :
2355 : ! Calculate the mean (in-precip.) of the hydrometeor in the 2nd PDF
2356 : ! component.
2357 : mu_hm_2(i,k) = ( hmm(i,k) - mixt_frac(i,k) * precip_frac_1(i,k) * mu_hm_1(i,k) ) &
2358 0 : / ( ( one - mixt_frac(i,k) ) * precip_frac_2(i,k) )
2359 :
2360 : ! Calculate the value of the ratio R (which is sigma_hm_2^2 / mu_hm_2^2),
2361 : ! where R = omicron * Rmax. The name of the variable used for R is
2362 : ! sigma_hm_2_sqd_on_mu_hm_2_sqd.
2363 0 : sigma_hm_2_sqd_on_mu_hm_2_sqd(i,k) = omicron * Rmax
2364 :
2365 : ! Calculate minimum allowable values for mu_hm_1 and mu_hm_2.
2366 0 : if ( hmm(i,k) / precip_frac(i,k) > hm_tol / precip_frac_1(i,k) ) then
2367 : mu_hm_1_min &
2368 : = min( hm_tol / precip_frac_1(i,k) &
2369 : + mu_hm_min_coef * ( hmm(i,k) / precip_frac(i,k) &
2370 : - hm_tol / precip_frac_1(i,k) ), &
2371 : ( hmm(i,k) - ( one - mixt_frac(i,k) ) * hm_tol ) &
2372 0 : / ( mixt_frac(i,k) * precip_frac_1(i,k) ) )
2373 : else ! hmm / precip_frac <= hm_tol / precip_frac_1
2374 : mu_hm_1_min = hm_tol / precip_frac_1(i,k)
2375 : end if
2376 0 : if ( hmm(i,k) / precip_frac(i,k) > hm_tol / precip_frac_2(i,k) ) then
2377 : mu_hm_2_min &
2378 : = min( hm_tol / precip_frac_2(i,k) &
2379 : + mu_hm_min_coef * ( hmm(i,k) / precip_frac(i,k) &
2380 : - hm_tol / precip_frac_2(i,k) ), &
2381 : ( hmm(i,k) - mixt_frac(i,k) * hm_tol ) &
2382 0 : / ( ( one - mixt_frac(i,k) ) * precip_frac_2(i,k) ) )
2383 : else ! hmm / precip_frac <= hm_tol / precip_frac_2
2384 : mu_hm_2_min = hm_tol / precip_frac_2(i,k)
2385 : end if
2386 :
2387 : ! Handle the "emergency" situation when the specified value of omicron is
2388 : ! too small for the value of <hm|_ip'^2> / <hm|_ip>^2, resulting in a
2389 : ! component mean that is too small (below tolerance value) or negative.
2390 0 : if ( mu_hm_1(i,k) < mu_hm_1_min ) then
2391 :
2392 : ! Set the value of mu_hm_1 to the threshold positive value.
2393 0 : mu_hm_1(i,k) = mu_hm_1_min
2394 :
2395 : ! Recalculate the mean (in-precip.) of the hydrometeor in the 2nd PDF
2396 : ! component.
2397 : mu_hm_2(i,k) = ( hmm(i,k) - mixt_frac(i,k) * precip_frac_1(i,k) * mu_hm_1(i,k) ) &
2398 0 : / ( ( one - mixt_frac(i,k) ) * precip_frac_2(i,k) )
2399 :
2400 : ! Recalculate the value of R ( sigma_hm_2^2 / mu_hm_2^2 ) in this
2401 : ! scenario.
2402 : ! R = ( <hm'^2> + <hm>^2 - a * f_p_1 * mu_hm_1^2
2403 : ! - ( 1 - a ) * f_p_2 * mu_hm_2^2 )
2404 : ! / ( a * f_p_1 * ( 1 + zeta ) * mu_hm_1^2
2405 : ! + ( 1 - a ) * f_p_2 * mu_hm_2^2 ).
2406 : sigma_hm_2_sqd_on_mu_hm_2_sqd(i,k) &
2407 : = ( hmp2(i,k) + hmm(i,k)**2 - mixt_frac(i,k) * precip_frac_1(i,k) * mu_hm_1(i,k)**2 &
2408 : - ( one - mixt_frac(i,k) ) * precip_frac_2(i,k) * mu_hm_2(i,k)**2 ) &
2409 : / ( mixt_frac(i,k) * precip_frac_1(i,k) &
2410 : * ( one + zeta_vrnce_rat ) * mu_hm_1(i,k)**2 &
2411 0 : + ( one - mixt_frac(i,k) ) * precip_frac_2(i,k) * mu_hm_2(i,k)**2 )
2412 :
2413 : ! Mathematically, this ratio can never be less than 0. In case
2414 : ! numerical round off error produces a negative value in extreme
2415 : ! cases, reset the value of R to 0.
2416 0 : if ( sigma_hm_2_sqd_on_mu_hm_2_sqd(i,k) < zero ) then
2417 0 : sigma_hm_2_sqd_on_mu_hm_2_sqd(i,k) = zero
2418 : end if
2419 :
2420 0 : elseif ( mu_hm_2(i,k) < mu_hm_2_min ) then
2421 :
2422 : ! Set the value of mu_hm_2 to the threshold positive value.
2423 0 : mu_hm_2(i,k) = mu_hm_2_min
2424 :
2425 : ! Recalculate the mean (in-precip.) of the hydrometeor in the 1st PDF
2426 : ! component.
2427 : mu_hm_1(i,k) = ( hmm(i,k) - ( one - mixt_frac(i,k) ) &
2428 : * precip_frac_2(i,k) * mu_hm_2(i,k) ) &
2429 0 : / ( mixt_frac(i,k) * precip_frac_1(i,k) )
2430 :
2431 : ! Recalculate the value of R ( sigma_hm_2^2 / mu_hm_2^2 ) in this
2432 : ! scenario.
2433 : ! R = ( <hm'^2> + <hm>^2 - a * f_p_1 * mu_hm_1^2
2434 : ! - ( 1 - a ) * f_p_2 * mu_hm_2^2 )
2435 : ! / ( a * f_p_1 * ( 1 + zeta ) * mu_hm_1^2
2436 : ! + ( 1 - a ) * f_p_2 * mu_hm_2^2 ).
2437 : sigma_hm_2_sqd_on_mu_hm_2_sqd(i,k) &
2438 : = ( hmp2(i,k) + hmm(i,k)**2 - mixt_frac(i,k) * precip_frac_1(i,k) * mu_hm_1(i,k)**2 &
2439 : - ( one - mixt_frac(i,k) ) * precip_frac_2(i,k) * mu_hm_2(i,k)**2 ) &
2440 : / ( mixt_frac(i,k) * precip_frac_1(i,k) &
2441 : * ( one + zeta_vrnce_rat ) * mu_hm_1(i,k)**2 &
2442 0 : + ( one - mixt_frac(i,k) ) * precip_frac_2(i,k) * mu_hm_2(i,k)**2 )
2443 :
2444 : ! Mathematically, this ratio can never be less than 0. In case
2445 : ! numerical round off error produces a negative value in extreme
2446 : ! cases, reset the value of R to 0.
2447 0 : if ( sigma_hm_2_sqd_on_mu_hm_2_sqd(i,k) < zero ) then
2448 0 : sigma_hm_2_sqd_on_mu_hm_2_sqd(i,k) = zero
2449 : end if
2450 :
2451 : end if
2452 :
2453 : ! Calculate the standard deviation (in-precip.) of the hydrometeor in the
2454 : ! 1st PDF component.
2455 : sigma_hm_1(i,k) = sqrt( sigma_hm_2_sqd_on_mu_hm_2_sqd(i,k) &
2456 : * ( one + zeta_vrnce_rat ) ) &
2457 0 : * mu_hm_1(i,k)
2458 :
2459 : ! Calculate the standard deviation (in-precip.) of the hydrometeor in the
2460 : ! 2nd PDF component.
2461 0 : sigma_hm_2(i,k) = sqrt( sigma_hm_2_sqd_on_mu_hm_2_sqd(i,k) ) * mu_hm_2(i,k)
2462 :
2463 : ! Calculate the mean of the hydrometeor in the 1st PDF component.
2464 0 : hm_1(i,k) = max( mu_hm_1(i,k) * precip_frac_1(i,k), hm_tol )
2465 :
2466 : ! Calculate the mean of the hydrometeor in the 1st PDF component.
2467 0 : hm_2(i,k) = max( mu_hm_2(i,k) * precip_frac_2(i,k), hm_tol )
2468 :
2469 : ! Calculate the ratio of sigma_hm_1^2 / mu_hm_1^2.
2470 0 : sigma_hm_1_sqd_on_mu_hm_1_sqd(i,k) = sigma_hm_1(i,k)**2 / mu_hm_1(i,k)**2
2471 :
2472 : ! The value of R, sigma_hm_2_sqd_on_mu_hm_2_sqd, has already been
2473 : ! calculated.
2474 :
2475 :
2476 0 : elseif ( hmm(i,k) >= hm_tol .and. precip_frac_1(i,k) >= precip_frac_tol(i) ) then
2477 :
2478 : ! Precipitation is found in the 1st PDF component, but not in the 2nd
2479 : ! PDF component (precip_frac_2 = 0).
2480 0 : mu_hm_1(i,k) = hmm(i,k) / ( mixt_frac(i,k) * precip_frac_1(i,k) )
2481 0 : mu_hm_2(i,k) = zero
2482 :
2483 : sigma_hm_1(i,k) = sqrt( max( ( hmp2(i,k) + hmm(i,k)**2 &
2484 : - mixt_frac(i,k) * precip_frac_1(i,k) * mu_hm_1(i,k)**2 ) &
2485 : / ( mixt_frac(i,k) * precip_frac_1(i,k) ), &
2486 0 : zero ) )
2487 0 : sigma_hm_2(i,k) = zero
2488 :
2489 0 : hm_1(i,k) = mu_hm_1(i,k) * precip_frac_1(i,k)
2490 0 : hm_2(i,k) = zero
2491 :
2492 0 : sigma_hm_1_sqd_on_mu_hm_1_sqd(i,k) = sigma_hm_1(i,k)**2 / mu_hm_1(i,k)**2
2493 : ! The ratio sigma_hm_2^2 / mu_hm_2^2 is undefined.
2494 0 : sigma_hm_2_sqd_on_mu_hm_2_sqd(i,k) = zero
2495 :
2496 :
2497 0 : elseif ( hmm(i,k) >= hm_tol .and. precip_frac_2(i,k) >= precip_frac_tol(i) ) then
2498 :
2499 : ! Precipitation is found in the 2nd PDF component, but not in the 1st
2500 : ! PDF component (precip_frac_1 = 0).
2501 0 : mu_hm_1(i,k) = zero
2502 0 : mu_hm_2(i,k) = hmm(i,k) / ( ( one - mixt_frac(i,k) ) * precip_frac_2(i,k) )
2503 :
2504 0 : sigma_hm_1(i,k) = zero
2505 : sigma_hm_2(i,k) &
2506 : = sqrt( max( ( hmp2(i,k) + hmm(i,k)**2 &
2507 : - ( one - mixt_frac(i,k) ) * precip_frac_2(i,k) * mu_hm_2(i,k)**2 ) &
2508 : / ( ( one - mixt_frac(i,k) ) * precip_frac_2(i,k) ), &
2509 0 : zero ) )
2510 :
2511 0 : hm_1(i,k) = zero
2512 0 : hm_2(i,k) = mu_hm_2(i,k) * precip_frac_2(i,k)
2513 :
2514 : ! The ratio sigma_hm_1^2 / mu_hm_1^2 is undefined.
2515 0 : sigma_hm_1_sqd_on_mu_hm_1_sqd(i,k) = zero
2516 0 : sigma_hm_2_sqd_on_mu_hm_2_sqd(i,k) = sigma_hm_2(i,k)**2 / mu_hm_2(i,k)**2
2517 :
2518 :
2519 : else ! hm < hm_tol or ( precip_frac_1 = 0 and precip_frac_2 = 0 ).
2520 :
2521 : ! Precipitation is not found in either PDF component.
2522 0 : mu_hm_1(i,k) = zero
2523 0 : mu_hm_2(i,k) = zero
2524 :
2525 0 : sigma_hm_1(i,k) = zero
2526 0 : sigma_hm_2(i,k) = zero
2527 :
2528 0 : hm_1(i,k) = zero
2529 0 : hm_2(i,k) = zero
2530 :
2531 : ! The ratio sigma_hm_1^2 / mu_hm_1^2 is undefined.
2532 0 : sigma_hm_1_sqd_on_mu_hm_1_sqd(i,k) = zero
2533 : ! The ratio sigma_hm_2^2 / mu_hm_2^2 is undefined.
2534 0 : sigma_hm_2_sqd_on_mu_hm_2_sqd(i,k) = zero
2535 :
2536 :
2537 : end if ! hmm >= hm_tol and precip_frac_1 >= precip_frac_tol
2538 : ! and precip_frac_2 >= precip_frac_tol
2539 : end do
2540 : end do
2541 :
2542 0 : return
2543 :
2544 : end subroutine calc_comp_mu_sigma_hm
2545 :
2546 : !=============================================================================
2547 0 : subroutine component_corr_w_x( nz, ngrdcol, &
2548 0 : rc_1, &
2549 0 : rc_2, &
2550 : corr_w_x_NN_cloud, corr_w_x_NN_below, &
2551 : iiPDF_type, &
2552 0 : corr_w_x_1, &
2553 0 : corr_w_x_2 )
2554 :
2555 : ! Description:
2556 : ! Calculates the correlation of w and x within the ith PDF component.
2557 : ! Here, x is a variable with a normally distributed individual marginal PDF,
2558 : ! such as chi or eta.
2559 :
2560 : ! References:
2561 : !-----------------------------------------------------------------------
2562 :
2563 : use constants_clubb, only: &
2564 : zero, &
2565 : rc_tol
2566 :
2567 : use model_flags, only: &
2568 : iiPDF_ADG1, & ! Variable(s)
2569 : iiPDF_ADG2, &
2570 : iiPDF_new_hybrid
2571 :
2572 : use clubb_precision, only: &
2573 : core_rknd ! Variable(s)
2574 :
2575 : implicit none
2576 :
2577 : ! Input Variables
2578 : integer, intent(in) :: &
2579 : nz, & ! Number of model vertical grid levels
2580 : ngrdcol ! Number of grid columns
2581 :
2582 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
2583 : rc_1, & ! Mean cloud water mixing ratio (1st PDF comp.) [kg/kg]
2584 : rc_2 ! Mean cloud water mixing ratio (2nd PDF comp.) [kg/kg]
2585 :
2586 : real( kind = core_rknd ), intent(in) :: &
2587 : corr_w_x_NN_cloud, & ! Corr. of w and x (ith PDF comp.); cloudy levs [-]
2588 : corr_w_x_NN_below ! Corr. of w and x (ith PDF comp.); clear levs [-]
2589 :
2590 : integer, intent(in) :: &
2591 : iiPDF_type ! Selected option for the two-component normal (double
2592 : ! Gaussian) PDF type to use for the w, rt, and theta-l (or
2593 : ! w, chi, and eta) portion of CLUBB's multivariate,
2594 : ! two-component PDF.
2595 :
2596 : ! Output Variable
2597 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
2598 : corr_w_x_1, & ! Correlation of w and x (1st PDF component) [-]
2599 : corr_w_x_2 ! Correlation of w and x (2nd PDF component) [-]
2600 :
2601 : ! Local Variables
2602 : integer :: i, k
2603 :
2604 : ! Correlation of w and x in the ith PDF component.
2605 :
2606 : ! The PDF variables chi and eta result from a transformation of the PDF
2607 : ! involving r_t and theta_l. The correlation of w and x (whether x is chi
2608 : ! or eta) depends on the correlation of w and r_t, the correlation of w and
2609 : ! theta_l, as well as the variances of r_t and theta_l, and other factors.
2610 : ! The correlation of w and x is subject to change at every vertical level
2611 : ! and model time step, and is calculated as part of the CLUBB PDF
2612 : ! parameters.
2613 :
2614 : ! The ADG1 PDF fixes the correlation of w and rt and the correlation of
2615 : ! w and theta_l to be 0, which means the correlation of w and chi and the
2616 : ! correlation of w and eta must also be 0.
2617 : if ( ( iiPDF_type == iiPDF_ADG1 .or. iiPDF_type == iiPDF_ADG2 &
2618 : .or. iiPDF_type == iiPDF_new_hybrid ) &
2619 0 : .and. l_follow_ADG1_PDF_standards ) then
2620 :
2621 0 : corr_w_x_1(:,:) = zero
2622 0 : corr_w_x_2(:,:) = zero
2623 :
2624 : else ! use prescribed paramter values
2625 :
2626 0 : do k = 1, nz
2627 0 : do i = 1, ngrdcol
2628 :
2629 0 : if ( rc_1(i,k) > rc_tol ) then
2630 0 : corr_w_x_1(i,k) = corr_w_x_NN_cloud
2631 : else
2632 0 : corr_w_x_1(i,k) = corr_w_x_NN_below
2633 : end if
2634 :
2635 0 : if ( rc_2(i,k) > rc_tol ) then
2636 0 : corr_w_x_2(i,k) = corr_w_x_NN_cloud
2637 : else
2638 0 : corr_w_x_2(i,k) = corr_w_x_NN_below
2639 : end if
2640 :
2641 : end do
2642 : end do
2643 :
2644 : end if ! iiPDF_type
2645 :
2646 0 : return
2647 :
2648 : end subroutine component_corr_w_x
2649 :
2650 : !=============================================================================
2651 0 : subroutine component_corr_chi_eta( nz, ngrdcol, &
2652 0 : rc_1, &
2653 0 : rc_2, &
2654 : corr_chi_eta_NN_cloud, &
2655 : corr_chi_eta_NN_below, &
2656 : l_limit_corr_chi_eta, &
2657 0 : corr_chi_eta_1, &
2658 0 : corr_chi_eta_2 )
2659 :
2660 : ! Description:
2661 : ! Calculates the correlation of chi (old s) and eta (old t) within the
2662 : ! ith PDF component.
2663 :
2664 : ! References:
2665 : !-----------------------------------------------------------------------
2666 :
2667 : use constants_clubb, only: &
2668 : rc_tol, &
2669 : max_mag_correlation
2670 :
2671 : use clubb_precision, only: &
2672 : core_rknd ! Constant
2673 :
2674 : implicit none
2675 :
2676 : ! Input Variables
2677 : integer, intent(in) :: &
2678 : nz, & ! Number of model vertical grid levels
2679 : ngrdcol ! Number of grid columns
2680 :
2681 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
2682 : rc_1, & ! Mean cloud water mix. rat. (1st PDF comp.) [kg/kg]&
2683 : rc_2 ! Mean cloud water mix. rat. (2nd PDF comp.) [kg/kg]
2684 :
2685 : real( kind = core_rknd ), intent(in) :: &
2686 : corr_chi_eta_NN_cloud, & ! Corr. of chi & eta (ith PDF comp.); cloudy [-]
2687 : corr_chi_eta_NN_below ! Corr. of chi & eta (ith PDF comp.); clear [-]
2688 :
2689 : logical, intent(in) :: &
2690 : l_limit_corr_chi_eta ! We must limit the correlation of chi and eta if
2691 : ! we are to take the Cholesky decomposition of the
2692 : ! resulting correlation matrix. This is because a
2693 : ! perfect correlation of chi and eta was found to
2694 : ! be unrealizable.
2695 :
2696 : ! Output Variable
2697 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
2698 : corr_chi_eta_1, & ! Correlation of chi and eta (1st PDF component) [-]
2699 : corr_chi_eta_2 ! Correlation of chi and eta (2nd PDF component) [-]
2700 :
2701 : ! Local Variables
2702 : integer :: i, k
2703 :
2704 : ! Correlation of chi (old s) and eta (old t) in the ith PDF component.
2705 :
2706 : ! The PDF variables chi and eta result from a transformation of the PDF
2707 : ! involving r_t and theta_l. The correlation of chi and eta depends on the
2708 : ! correlation of r_t and theta_l, as well as the variances of r_t and
2709 : ! theta_l, and other factors. The correlation of chi and eta is subject to
2710 : ! change at every vertical level and model time step, and is calculated as
2711 : ! part of the CLUBB PDF parameters.
2712 :
2713 : ! WARNING: this code is inconsistent with the rest of CLUBB's PDF. This
2714 : ! code is necessary because SILHS is lazy and wussy, and only
2715 : ! wants to declare correlation arrays at the start of the model
2716 : ! run, rather than updating them throughout the model run.
2717 :
2718 :
2719 0 : do k = 1, nz
2720 0 : do i = 1, ngrdcol
2721 :
2722 0 : if ( rc_1(i,k) > rc_tol ) then
2723 0 : corr_chi_eta_1(i,k) = corr_chi_eta_NN_cloud
2724 : else
2725 0 : corr_chi_eta_1(i,k) = corr_chi_eta_NN_below
2726 : end if
2727 : end do
2728 : end do
2729 :
2730 0 : do k = 1, nz
2731 0 : do i = 1, ngrdcol
2732 0 : if ( rc_2(i,k) > rc_tol ) then
2733 0 : corr_chi_eta_2(i,k) = corr_chi_eta_NN_cloud
2734 : else
2735 0 : corr_chi_eta_2(i,k) = corr_chi_eta_NN_below
2736 : end if
2737 : end do
2738 : end do
2739 :
2740 :
2741 :
2742 : ! We cannot have a perfect correlation of chi (old s) and eta (old t) if we
2743 : ! plan to decompose this matrix and we don't want the Cholesky_factor code
2744 : ! to throw a fit.
2745 0 : if ( l_limit_corr_chi_eta ) then
2746 :
2747 0 : do k = 1, nz
2748 0 : do i = 1, ngrdcol
2749 0 : corr_chi_eta_1(i,k) = max( min( corr_chi_eta_1(i,k), max_mag_correlation ), &
2750 0 : -max_mag_correlation )
2751 : end do
2752 : end do
2753 :
2754 0 : do k = 1, nz
2755 0 : do i = 1, ngrdcol
2756 0 : corr_chi_eta_2(i,k) = max( min( corr_chi_eta_2(i,k), max_mag_correlation ), &
2757 0 : -max_mag_correlation )
2758 : end do
2759 : end do
2760 :
2761 : end if
2762 :
2763 0 : return
2764 :
2765 : end subroutine component_corr_chi_eta
2766 :
2767 : !=============================================================================
2768 0 : subroutine component_corr_w_hm_n_ip( nz, ngrdcol, &
2769 0 : corr_w_hm_1_n_in, rc_1, &
2770 0 : corr_w_hm_2_n_in, rc_2, &
2771 : corr_w_hm_n_NL_cloud, &
2772 : corr_w_hm_n_NL_below, &
2773 : l_calc_w_corr, &
2774 0 : corr_w_hm_1_n, &
2775 0 : corr_w_hm_2_n )
2776 :
2777 : ! Description:
2778 : ! Calculates the in-precip correlation of w and the natural logarithm of a
2779 : ! hydrometeor species within the ith PDF component.
2780 :
2781 : ! References:
2782 : !-----------------------------------------------------------------------
2783 :
2784 : use constants_clubb, only: &
2785 : rc_tol
2786 :
2787 : use clubb_precision, only: &
2788 : core_rknd ! Variable(s)
2789 :
2790 : implicit none
2791 :
2792 : ! Input Variables
2793 : integer, intent(in) :: &
2794 : nz, & ! Number of model vertical grid levels
2795 : ngrdcol ! Number of grid columns
2796 :
2797 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
2798 : corr_w_hm_1_n_in, & ! Correlation of w and ln hm (1st PDF comp.) ip [-]
2799 : corr_w_hm_2_n_in, & ! Correlation of w and ln hm (2nd PDF comp.) ip [-]
2800 : rc_1, & ! Mean cloud water mix. ratio (1st PDF comp.) [kg/kg]
2801 : rc_2 ! Mean cloud water mix. ratio (2nd PDF comp.) [kg/kg]
2802 :
2803 : real( kind = core_rknd ), intent(in) :: &
2804 : corr_w_hm_n_NL_cloud, & ! Corr. of w & ln hm (ith PDF comp.) ip; cloud [-]
2805 : corr_w_hm_n_NL_below ! Corr. of w & ln hm (ith PDF comp.) ip; clear [-]
2806 :
2807 : logical, intent(in) :: &
2808 : l_calc_w_corr ! Calculate the correlations between w and the hydrometeors
2809 :
2810 : ! Output Variables
2811 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
2812 : corr_w_hm_1_n, & ! Correlation of w and ln hm (1st PDF component) ip [-]
2813 : corr_w_hm_2_n ! Correlation of w and ln hm (2nd PDF component) ip [-]
2814 :
2815 : ! Local Variables
2816 : integer :: i, k
2817 :
2818 :
2819 : ! Correlation (in-precip) of w and the natural logarithm of the hydrometeor
2820 : ! in the ith PDF component.
2821 0 : if ( l_calc_w_corr ) then
2822 :
2823 0 : corr_w_hm_1_n(:,:) = corr_w_hm_1_n_in(:,:)
2824 0 : corr_w_hm_2_n(:,:) = corr_w_hm_2_n_in(:,:)
2825 :
2826 :
2827 : else
2828 :
2829 0 : do k = 1, nz
2830 0 : do i = 1, ngrdcol
2831 :
2832 0 : if ( rc_1(i,k) > rc_tol ) then
2833 0 : corr_w_hm_1_n(i,k) = corr_w_hm_n_NL_cloud
2834 : else
2835 0 : corr_w_hm_1_n(i,k) = corr_w_hm_n_NL_below
2836 : end if
2837 :
2838 0 : if ( rc_2(i,k) > rc_tol ) then
2839 0 : corr_w_hm_2_n(i,k) = corr_w_hm_n_NL_cloud
2840 : else
2841 0 : corr_w_hm_2_n(i,k) = corr_w_hm_n_NL_below
2842 : end if
2843 : end do
2844 : end do
2845 :
2846 : end if ! l_calc_w_corr
2847 :
2848 0 : return
2849 :
2850 : end subroutine component_corr_w_hm_n_ip
2851 :
2852 : !=============================================================================
2853 0 : subroutine component_corr_x_hm_n_ip( nz, ngrdcol, &
2854 0 : rc_1, &
2855 0 : rc_2, &
2856 : corr_x_hm_n_NL_cloud, &
2857 : corr_x_hm_n_NL_below, &
2858 0 : corr_x_hm_1_n, &
2859 0 : corr_x_hm_2_n )
2860 :
2861 : ! Description:
2862 : ! Calculates the in-precip correlation of x and a hydrometeor species
2863 : ! within the ith PDF component. Here, x is a variable with a normally
2864 : ! distributed individual marginal PDF, such as chi or eta.
2865 :
2866 : ! References:
2867 : !-----------------------------------------------------------------------
2868 :
2869 : use constants_clubb, only: &
2870 : rc_tol
2871 :
2872 : use clubb_precision, only: &
2873 : core_rknd ! Variable(s)
2874 :
2875 : implicit none
2876 :
2877 : ! Input Variables
2878 : integer, intent(in) :: &
2879 : nz, & ! Number of model vertical grid levels
2880 : ngrdcol ! Number of grid columns
2881 :
2882 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
2883 : rc_1, & ! Mean cloud water mixing ratio (1st PDF comp.) [kg/kg]
2884 : rc_2 ! Mean cloud water mixing ratio (2nd PDF comp.) [kg/kg]
2885 :
2886 :
2887 : real( kind = core_rknd ), intent(in) :: &
2888 : corr_x_hm_n_NL_cloud, & ! Corr. of x and ln hm (ith PDF comp.) ip [-]
2889 : corr_x_hm_n_NL_below ! Corr. of x and ln hm (ith PDF comp.) ip [-]
2890 :
2891 : ! Output Variables
2892 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
2893 : corr_x_hm_1_n, & ! Correlation of x and ln hm (1st PDF component) ip [-]
2894 : corr_x_hm_2_n ! Correlation of x and ln hm (2nd PDF component) ip [-]
2895 :
2896 : ! Local Variables
2897 : integer :: i, k
2898 :
2899 0 : do k = 1, nz
2900 0 : do i = 1, ngrdcol
2901 :
2902 0 : if ( rc_1(i,k) > rc_tol ) then
2903 0 : corr_x_hm_1_n(i,k) = corr_x_hm_n_NL_cloud
2904 : else
2905 0 : corr_x_hm_1_n(i,k) = corr_x_hm_n_NL_below
2906 : end if
2907 :
2908 0 : if ( rc_2(i,k) > rc_tol ) then
2909 0 : corr_x_hm_2_n(i,k) = corr_x_hm_n_NL_cloud
2910 : else
2911 0 : corr_x_hm_2_n(i,k) = corr_x_hm_n_NL_below
2912 : end if
2913 :
2914 : end do
2915 : end do
2916 :
2917 0 : return
2918 :
2919 : end subroutine component_corr_x_hm_n_ip
2920 :
2921 : !=============================================================================
2922 0 : subroutine component_corr_hmx_hmy_n_ip( nz, ngrdcol, &
2923 0 : rc_1, &
2924 0 : rc_2, &
2925 : corr_hmx_hmy_n_LL_cloud, &
2926 : corr_hmx_hmy_n_LL_below, &
2927 0 : corr_hmx_hmy_1_n, &
2928 0 : corr_hmx_hmy_2_n )
2929 :
2930 : ! Description:
2931 : ! Calculates the in-precip correlation of the natural logarithms of
2932 : ! hydrometeor x and hydrometeor y within the ith PDF component.
2933 :
2934 : ! References:
2935 : !-----------------------------------------------------------------------
2936 :
2937 : use constants_clubb, only: &
2938 : rc_tol
2939 :
2940 : use clubb_precision, only: &
2941 : core_rknd ! Variable(s)
2942 :
2943 : implicit none
2944 :
2945 : ! Input Variables
2946 : integer, intent(in) :: &
2947 : nz, & ! Number of model vertical grid levels
2948 : ngrdcol ! Number of grid columns
2949 :
2950 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
2951 : rc_1, & ! Mean cloud water mixing ratio (ith PDF comp.) [kg/kg]
2952 : rc_2 ! Mean cloud water mixing ratio (ith PDF comp.) [kg/kg]
2953 :
2954 : real( kind = core_rknd ), intent(in) :: &
2955 : corr_hmx_hmy_n_LL_cloud, & ! Corr.: ln hmx & ln hmy (ith PDF comp.) ip [-]
2956 : corr_hmx_hmy_n_LL_below ! Corr.: ln hmx & ln hmy (ith PDF comp.) ip [-]
2957 :
2958 : ! Output Variable
2959 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
2960 : corr_hmx_hmy_1_n, & ! Corr. of ln hmx & ln hmy (ith PDF comp.) ip [-]
2961 : corr_hmx_hmy_2_n ! Corr. of ln hmx & ln hmy (ith PDF comp.) ip [-]
2962 :
2963 : ! Local Variable
2964 : integer :: i, k
2965 :
2966 :
2967 : ! Correlation (in-precip) of the natural logarithms of hydrometeor x and
2968 : ! hydrometeor y in the ith PDF component.
2969 0 : do k = 1, nz
2970 0 : do i = 1, ngrdcol
2971 :
2972 0 : if ( rc_1(i,k) > rc_tol ) then
2973 0 : corr_hmx_hmy_1_n(i,k) = corr_hmx_hmy_n_LL_cloud
2974 : else
2975 0 : corr_hmx_hmy_1_n(i,k) = corr_hmx_hmy_n_LL_below
2976 : end if
2977 :
2978 0 : if ( rc_2(i,k) > rc_tol ) then
2979 0 : corr_hmx_hmy_2_n(i,k) = corr_hmx_hmy_n_LL_cloud
2980 : else
2981 0 : corr_hmx_hmy_2_n(i,k) = corr_hmx_hmy_n_LL_below
2982 : end if
2983 :
2984 : end do
2985 : end do
2986 :
2987 0 : return
2988 :
2989 : end subroutine component_corr_hmx_hmy_n_ip
2990 :
2991 : !=============================================================================
2992 0 : subroutine component_corr_eta_hm_n_ip( nz, ngrdcol, &
2993 0 : corr_chi_eta_1, &
2994 0 : corr_chi_hm_n_1, &
2995 0 : corr_chi_eta_2, &
2996 0 : corr_chi_hm_n_2, &
2997 0 : corr_eta_hm_n_1, &
2998 0 : corr_eta_hm_n_2)
2999 : ! Description:
3000 : ! Estimates the correlation of eta and the natural logarithm of a
3001 : ! hydrometeor species using the correlation of chi and eta and the
3002 : ! correlation of chi and the natural logarithm of the hydrometeor. This
3003 : ! facilities the Cholesky decomposability of the correlation array that will
3004 : ! inevitably be decomposed for SILHS purposes. Without this estimation, we
3005 : ! have found that the resulting correlation matrix cannot be decomposed.
3006 :
3007 : ! References:
3008 : !-----------------------------------------------------------------------
3009 :
3010 : use clubb_precision, only: &
3011 : core_rknd ! Constant
3012 :
3013 : implicit none
3014 :
3015 : ! Input Variables
3016 : integer, intent(in) :: &
3017 : nz, & ! Number of model vertical grid levels
3018 : ngrdcol ! Number of grid columns
3019 :
3020 :
3021 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
3022 : corr_chi_eta_1, & ! Component correlation of chi and eta [-]
3023 : corr_chi_eta_2, & ! Component correlation of chi and eta [-]
3024 : corr_chi_hm_n_1, & ! Component correlation of chi and ln hm [-]
3025 : corr_chi_hm_n_2 ! Component correlation of chi and ln hm [-]
3026 :
3027 : ! Output Variables
3028 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
3029 : corr_eta_hm_n_1, & ! Component correlation of eta and ln hm [-]
3030 : corr_eta_hm_n_2 ! Component correlation of eta and ln hm [-]
3031 :
3032 :
3033 0 : corr_eta_hm_n_1 = corr_chi_eta_1 * corr_chi_hm_n_1
3034 0 : corr_eta_hm_n_2 = corr_chi_eta_2 * corr_chi_hm_n_2
3035 :
3036 :
3037 0 : return
3038 :
3039 : end subroutine component_corr_eta_hm_n_ip
3040 :
3041 : !=============================================================================
3042 0 : subroutine norm_transform_mean_stdev( nz, ngrdcol, &
3043 : hydromet_dim, hm_metadata, &
3044 0 : hm_1, hm_2, &
3045 0 : Ncnm, pdf_dim, &
3046 0 : mu_x_1, mu_x_2, &
3047 0 : sigma_x_1, sigma_x_2, &
3048 0 : sigma2_on_mu2_ip_1, &
3049 0 : sigma2_on_mu2_ip_2, &
3050 : l_const_Nc_in_cloud, &
3051 0 : mu_x_1_n, mu_x_2_n, &
3052 0 : sigma_x_1_n, sigma_x_2_n )
3053 :
3054 : ! Description:
3055 : ! Transforms the means and the standard deviations of PDF variables that
3056 : ! have assumed lognormal distributions -- which are precipitating
3057 : ! hydrometeors (in precipitation) and N_cn -- to normal space for each PDF
3058 : ! component.
3059 :
3060 : ! References:
3061 : !-----------------------------------------------------------------------
3062 :
3063 : use constants_clubb, only: &
3064 : Ncn_tol, & ! Constant(s)
3065 : zero
3066 :
3067 : use pdf_utilities, only: &
3068 : mean_L2N, & ! Procedure(s)
3069 : stdev_L2N
3070 :
3071 : use index_mapping, only: &
3072 : pdf2hydromet_idx ! Procedure(s)
3073 :
3074 : use corr_varnce_module, only: &
3075 : hm_metadata_type
3076 :
3077 : use clubb_precision, only: &
3078 : core_rknd ! Variable(s)
3079 :
3080 : implicit none
3081 :
3082 : !--------------------------- Input Variables ---------------------------
3083 : integer, intent(in) :: &
3084 : nz, & ! Number of model vertical grid levels
3085 : pdf_dim, & ! Number of variables in CLUBB's PDF
3086 : ngrdcol, & ! Number of grid columns
3087 : hydromet_dim ! Number of hydrometeor species
3088 :
3089 : type (hm_metadata_type), intent(in) :: &
3090 : hm_metadata
3091 :
3092 : real( kind = core_rknd ), dimension(ngrdcol,nz,hydromet_dim), intent(in) :: &
3093 : hm_1, & ! Mean of a precip. hydrometeor (1st PDF component) [units vary]
3094 : hm_2 ! Mean of a precip. hydrometeor (2nd PDF component) [units vary]
3095 :
3096 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
3097 : Ncnm ! Mean cloud nuclei concentration, < N_cn > [num/kg]
3098 :
3099 : real( kind = core_rknd ), dimension(ngrdcol,nz,pdf_dim), intent(in) :: &
3100 : mu_x_1, & ! Mean array of PDF vars. (1st PDF component) [units vary]
3101 : mu_x_2, & ! Mean array of PDF vars. (2nd PDF component) [units vary]
3102 : sigma_x_1, & ! Standard deviation array of PDF vars (comp. 1) [units vary]
3103 : sigma_x_2 ! Standard deviation array of PDF vars (comp. 2) [units vary]
3104 :
3105 : real( kind = core_rknd ), dimension(ngrdcol,nz,pdf_dim), intent(in) :: &
3106 : sigma2_on_mu2_ip_1, & ! Prescribed ratio array: sigma_hm_1^2/mu_hm_1^2 [-]
3107 : sigma2_on_mu2_ip_2 ! Prescribed ratio array: sigma_hm_2^2/mu_hm_2^2 [-]
3108 :
3109 : logical, intent(in) :: &
3110 : l_const_Nc_in_cloud ! Use a constant cloud droplet conc. within cloud (K&K)
3111 :
3112 : !--------------------------- Output Variables ---------------------------
3113 : real( kind = core_rknd ), dimension(ngrdcol,nz,pdf_dim), intent(out) :: &
3114 : mu_x_1_n, & ! Mean array (normal space): PDF vars. (comp. 1) [un. vary]
3115 : mu_x_2_n, & ! Mean array (normal space): PDF vars. (comp. 2) [un. vary]
3116 : sigma_x_1_n, & ! Std. dev. array (normal space): PDF vars (comp. 1) [u.v.]
3117 : sigma_x_2_n ! Std. dev. array (normal space): PDF vars (comp. 2) [u.v.]
3118 :
3119 : !--------------------------- Local Variable ---------------------------
3120 : integer :: ivar, hm_idx, k, i ! Indices
3121 :
3122 : ! Just to avoid typing hm_metadata%iiPDF_x everywhere
3123 : integer :: &
3124 : iiPDF_w, &
3125 : iiPDF_chi, &
3126 : iiPDF_eta, &
3127 : iiPDF_Ncn
3128 :
3129 : !--------------------------- Begin Code ---------------------------
3130 :
3131 0 : iiPDF_chi = hm_metadata%iiPDF_chi
3132 0 : iiPDF_eta = hm_metadata%iiPDF_eta
3133 0 : iiPDF_w = hm_metadata%iiPDF_w
3134 0 : iiPDF_Ncn = hm_metadata%iiPDF_Ncn
3135 :
3136 : ! The means and standard deviations in each PDF component of w, chi (old s),
3137 : ! and eta (old t) do not need to be transformed to normal space, since w,
3138 : ! chi, and eta already follow assumed normal distributions in each PDF
3139 : ! component. The normal space means and standard deviations are the same as
3140 : ! the actual means and standard deviations.
3141 0 : do k = 1, nz
3142 0 : mu_x_1_n(:,k,iiPDF_chi) = mu_x_1(:,k,iiPDF_chi)
3143 0 : mu_x_1_n(:,k,iiPDF_eta) = mu_x_1(:,k,iiPDF_eta)
3144 0 : mu_x_1_n(:,k,iiPDF_w) = mu_x_1(:,k,iiPDF_w)
3145 : end do
3146 :
3147 0 : do k = 1, nz
3148 0 : mu_x_2_n(:,k,iiPDF_chi) = mu_x_2(:,k,iiPDF_chi)
3149 0 : mu_x_2_n(:,k,iiPDF_eta) = mu_x_2(:,k,iiPDF_eta)
3150 0 : mu_x_2_n(:,k,iiPDF_w) = mu_x_2(:,k,iiPDF_w)
3151 : end do
3152 :
3153 0 : do k = 1, nz
3154 0 : sigma_x_1_n(:,k,iiPDF_chi) = sigma_x_1(:,k,iiPDF_chi)
3155 0 : sigma_x_1_n(:,k,iiPDF_eta) = sigma_x_1(:,k,iiPDF_eta)
3156 0 : sigma_x_1_n(:,k,iiPDF_w) = sigma_x_1(:,k,iiPDF_w)
3157 : end do
3158 :
3159 0 : do k = 1, nz
3160 0 : sigma_x_2_n(:,k,iiPDF_chi) = sigma_x_2(:,k,iiPDF_chi)
3161 0 : sigma_x_2_n(:,k,iiPDF_eta) = sigma_x_2(:,k,iiPDF_eta)
3162 0 : sigma_x_2_n(:,k,iiPDF_w) = sigma_x_2(:,k,iiPDF_w)
3163 : end do
3164 :
3165 : !!! Transform the mean and standard deviation to normal space in each PDF
3166 : !!! component for variables that have an assumed lognormal distribution,
3167 : !!! given the mean and standard deviation in each PDF component for those
3168 : !!! variables. A precipitating hydrometeor has an assumed lognormal
3169 : !!! distribution in precipitation in each PDF component. Simplified cloud
3170 : !!! nuclei concentration, N_cn, has an assumed lognormal distribution in
3171 : !!! each PDF component, and furthermore, mu_Ncn_1 = mu_Ncn_2 and
3172 : !!! sigma_Ncn_1 = sigma_Ncn_2, so N_cn has an assumed single lognormal
3173 : !!! distribution over the entire domain.
3174 :
3175 : ! Normal space mean of simplified cloud nuclei concentration, N_cn,
3176 : ! in PDF component 1.
3177 0 : do k = 1, nz
3178 0 : do i = 1, ngrdcol
3179 :
3180 0 : if ( Ncnm(i,k) >= Ncn_tol ) then
3181 :
3182 0 : mu_x_1_n(i,k,iiPDF_Ncn) = mean_L2N( mu_x_1(i,k,iiPDF_Ncn), &
3183 0 : sigma2_on_mu2_ip_1(i,k,iiPDF_Ncn) )
3184 :
3185 : else
3186 :
3187 : ! Mean simplified cloud nuclei concentration in PDF component 1 is less
3188 : ! than the tolerance amount. It is considered to have a value of 0.
3189 : ! There are not any cloud nuclei or cloud at this grid level. The value
3190 : ! of mu_Ncn_1_n should be -inf. It will be set to -huge for purposes of
3191 : ! assigning it a value.
3192 0 : mu_x_1_n(i,k,iiPDF_Ncn) = -huge( mu_x_1(i,k,iiPDF_Ncn) )
3193 :
3194 : end if
3195 :
3196 : end do
3197 : end do
3198 :
3199 : ! Normal space mean of simplified cloud nuclei concentration, N_cn,
3200 : ! in PDF component 2.
3201 0 : do k = 1, nz
3202 0 : do i = 1, ngrdcol
3203 :
3204 0 : if ( Ncnm(i,k) >= Ncn_tol ) then
3205 :
3206 0 : mu_x_2_n(i,k,iiPDF_Ncn) = mean_L2N( mu_x_2(i,k,iiPDF_Ncn), &
3207 0 : sigma2_on_mu2_ip_1(i,k,iiPDF_Ncn) )
3208 :
3209 : else
3210 :
3211 : ! Mean simplified cloud nuclei concentration in PDF component 1 is less
3212 : ! than the tolerance amount. It is considered to have a value of 0.
3213 : ! There are not any cloud nuclei or cloud at this grid level. The value
3214 : ! of mu_Ncn_1_n should be -inf. It will be set to -huge for purposes of
3215 : ! assigning it a value.
3216 0 : mu_x_2_n(i,k,iiPDF_Ncn) = -huge( mu_x_2(i,k,iiPDF_Ncn) )
3217 :
3218 : end if
3219 :
3220 : end do
3221 : end do
3222 :
3223 : ! Normal space standard deviation of simplified cloud nuclei concentration,
3224 : ! N_cn, in PDF components 1 and 2.
3225 0 : if ( l_const_Nc_in_cloud ) then
3226 : ! Ncn does not vary in the grid box.
3227 0 : sigma_x_1_n(:,:,iiPDF_Ncn) = zero
3228 0 : sigma_x_2_n(:,:,iiPDF_Ncn) = zero
3229 : else
3230 : ! Ncn (perhaps) varies in the grid box.
3231 0 : sigma_x_1_n(:,:,iiPDF_Ncn) = stdev_L2N( sigma2_on_mu2_ip_1(:,:,iiPDF_Ncn) )
3232 0 : sigma_x_2_n(:,:,iiPDF_Ncn) = stdev_L2N( sigma2_on_mu2_ip_2(:,:,iiPDF_Ncn) )
3233 : end if
3234 :
3235 : ! Normal space precipitating hydrometeor means and standard deviations.
3236 0 : do ivar = iiPDF_Ncn+1, pdf_dim, 1
3237 :
3238 0 : hm_idx = pdf2hydromet_idx(ivar,hm_metadata)
3239 :
3240 : ! Normal space mean of a precipitating hydrometeor, hm, in PDF
3241 : ! component 1.
3242 0 : do k = 1, nz
3243 0 : do i = 1, ngrdcol
3244 :
3245 0 : if ( hm_1(i,k,hm_idx) >= hm_metadata%hydromet_tol(hm_idx) ) then
3246 :
3247 0 : mu_x_1_n(i,k,ivar) = mean_L2N( mu_x_1(i,k,ivar), &
3248 0 : sigma2_on_mu2_ip_1(i,k,ivar) )
3249 :
3250 : else
3251 :
3252 : ! The mean of a precipitating hydrometeor in PDF component 1 is less
3253 : ! than its tolerance amount. It is considered to have a value of 0.
3254 : ! There is not any of this precipitating hydrometeor in the 1st PDF
3255 : ! component at this grid level. The in-precip mean of this
3256 : ! precipitating hydrometeor (1st PDF component) is also 0. The value
3257 : ! of mu_hm_1_n should be -inf. It will be set to -huge for purposes
3258 : ! of assigning it a value.
3259 0 : mu_x_1_n(i,k,ivar) = -huge( mu_x_1(i,k,ivar) )
3260 :
3261 : end if
3262 :
3263 : end do
3264 : end do
3265 :
3266 : ! Normal space standard deviation of a precipitating hydrometeor, hm, in
3267 : ! PDF component 1.
3268 0 : sigma_x_1_n(:,:,ivar) = stdev_L2N( sigma2_on_mu2_ip_1(:,:,ivar) )
3269 :
3270 : ! Normal space mean of a precipitating hydrometeor, hm, in PDF
3271 : ! component 2.
3272 0 : do k = 1, nz
3273 0 : do i = 1, ngrdcol
3274 :
3275 0 : if ( hm_2(i,k,hm_idx) >= hm_metadata%hydromet_tol(hm_idx) ) then
3276 :
3277 : mu_x_2_n(i,k,ivar) = mean_L2N( mu_x_2(i,k,ivar), &
3278 0 : sigma2_on_mu2_ip_2(i,k,ivar) )
3279 :
3280 : else
3281 :
3282 : ! The mean of a precipitating hydrometeor in PDF component 2 is less
3283 : ! than its tolerance amount. It is considered to have a value of 0.
3284 : ! There is not any of this precipitating hydrometeor in the 2nd PDF
3285 : ! component at this grid level. The in-precip mean of this
3286 : ! precipitating hydrometeor (2nd PDF component) is also 0. The value
3287 : ! of mu_hm_2_n should be -inf. It will be set to -huge for purposes
3288 : ! of assigning it a value.
3289 0 : mu_x_2_n(i,k,ivar) = -huge( mu_x_2(i,k,ivar) )
3290 :
3291 : end if
3292 :
3293 : end do
3294 : end do
3295 :
3296 : ! Normal space standard deviation of a precipitating hydrometeor, hm, in
3297 : ! PDF component 2.
3298 0 : sigma_x_2_n(:,:,ivar) = stdev_L2N( sigma2_on_mu2_ip_2(:,:,ivar) )
3299 :
3300 : end do ! ivar = iiPDF_Ncn+1, pdf_dim, 1
3301 :
3302 :
3303 0 : return
3304 :
3305 : end subroutine norm_transform_mean_stdev
3306 :
3307 : !=============================================================================
3308 0 : subroutine denorm_transform_corr( nz, ngrdcol, pdf_dim, hm_metadata, &
3309 0 : sigma_x_1_n, sigma_x_2_n, &
3310 0 : sigma2_on_mu2_ip_1, sigma2_on_mu2_ip_2, &
3311 0 : corr_array_1_n, &
3312 0 : corr_array_2_n, &
3313 0 : corr_array_1, corr_array_2 )
3314 :
3315 : ! Description:
3316 : ! Calculates the true or "real-space" correlations between PDF variables,
3317 : ! where at least one of the variables that is part of a correlation has an
3318 : ! assumed lognormal distribution -- which are the precipitating hydrometeors
3319 : ! (in precipitation) and N_cn.
3320 :
3321 : ! References:
3322 : !-----------------------------------------------------------------------
3323 :
3324 : use constants_clubb, only : &
3325 : one
3326 :
3327 : use pdf_utilities, only: &
3328 : corr_NN2NL, & ! Procedure(s)
3329 : corr_NN2LL
3330 :
3331 : use corr_varnce_module, only: &
3332 : hm_metadata_type
3333 :
3334 : use clubb_precision, only: &
3335 : core_rknd ! Variable(s)
3336 :
3337 : implicit none
3338 :
3339 : !--------------------------- Input Variables ---------------------------
3340 : integer, intent(in) :: &
3341 : nz, & ! Number of vertical levels
3342 : pdf_dim, & ! Number of variables
3343 : ngrdcol ! Number of grid columns
3344 :
3345 : type (hm_metadata_type), intent(in) :: &
3346 : hm_metadata
3347 :
3348 : real( kind = core_rknd ), dimension(ngrdcol,nz,pdf_dim), intent(in) :: &
3349 : sigma_x_1_n, & ! Std. dev. array (normal space): PDF vars (comp. 1) [u.v.]
3350 : sigma_x_2_n ! Std. dev. array (normal space): PDF vars (comp. 2) [u.v.]
3351 :
3352 : real( kind = core_rknd ), dimension(ngrdcol,nz,pdf_dim), intent(in) :: &
3353 : sigma2_on_mu2_ip_1, & ! Ratio array sigma_hm_1^2/mu_hm_1^2 [-]
3354 : sigma2_on_mu2_ip_2 ! Ratio array sigma_hm_2^2/mu_hm_2^2 [-]
3355 :
3356 : real( kind = core_rknd ), dimension(ngrdcol,nz,pdf_dim,pdf_dim), &
3357 : intent(in) :: &
3358 : corr_array_1_n, & ! Corr. array (normal space) of PDF vars. (comp. 1) [-]
3359 : corr_array_2_n ! Corr. array (normal space) of PDF vars. (comp. 2) [-]
3360 :
3361 : !--------------------------- Output Variables ---------------------------
3362 : real( kind = core_rknd ), dimension(ngrdcol,nz,pdf_dim,pdf_dim), &
3363 : intent(out) :: &
3364 : corr_array_1, & ! Correlation array of PDF vars. (comp. 1) [-]
3365 : corr_array_2 ! Correlation array of PDF vars. (comp. 2) [-]
3366 :
3367 : !--------------------------- Local Variables ---------------------------
3368 : integer :: ivar, jvar ! Loop indices
3369 :
3370 : ! Just to avoid typing hm_metadata%iiPDF_x everywhere
3371 : integer :: &
3372 : iiPDF_chi, &
3373 : iiPDF_eta, &
3374 : iiPDF_w, &
3375 : iiPDF_Ncn
3376 :
3377 : !------------------------ Begin Code ------------------------
3378 :
3379 0 : iiPDF_chi = hm_metadata%iiPDF_chi
3380 0 : iiPDF_eta = hm_metadata%iiPDF_eta
3381 0 : iiPDF_w = hm_metadata%iiPDF_w
3382 0 : iiPDF_Ncn = hm_metadata%iiPDF_Ncn
3383 :
3384 : ! Initialize diagonal elements to one
3385 0 : do ivar = 1, pdf_dim
3386 0 : corr_array_1(:,:,ivar,ivar) = one
3387 0 : corr_array_2(:,:,ivar,ivar) = one
3388 : end do
3389 :
3390 : ! The correlations in each PDF component between two of w, chi (old s), and
3391 : ! eta (old t) do not need to be transformed to standard space, since w, chi,
3392 : ! and eta follow assumed normal distributions in each PDF component. The
3393 : ! normal space correlations between any two of these variables are the same
3394 : ! as the actual correlations.
3395 0 : corr_array_1(:,:,iiPDF_eta,iiPDF_chi) = corr_array_1_n(:,:,iiPDF_eta,iiPDF_chi)
3396 0 : corr_array_1(:,:,iiPDF_w,iiPDF_chi) = corr_array_1_n(:,:,iiPDF_w,iiPDF_chi)
3397 0 : corr_array_1(:,:,iiPDF_w,iiPDF_eta) = corr_array_1_n(:,:,iiPDF_w,iiPDF_eta)
3398 :
3399 0 : corr_array_2(:,:,iiPDF_eta,iiPDF_chi) = corr_array_2_n(:,:,iiPDF_eta,iiPDF_chi)
3400 0 : corr_array_2(:,:,iiPDF_w,iiPDF_chi) = corr_array_2_n(:,:,iiPDF_w,iiPDF_chi)
3401 0 : corr_array_2(:,:,iiPDF_w,iiPDF_eta) = corr_array_2_n(:,:,iiPDF_w,iiPDF_eta)
3402 :
3403 : !!! Calculate the true correlation of variables that have an assumed normal
3404 : !!! distribution and variables that have an assumed lognormal distribution
3405 : !!! for the ith PDF component, given their normal space correlation and the
3406 : !!! normal space standard deviation of the variable with the assumed
3407 : !!! lognormal distribution.
3408 :
3409 : ! Transform the correlations between chi/eta/w and N_cn to standard space.
3410 :
3411 : ! Transform the correlation of chi (old s) and N_cn to standard space in PDF component 1.
3412 : call corr_NN2NL( nz, ngrdcol, & ! intent(in)
3413 : corr_array_1_n(:,:,iiPDF_Ncn,iiPDF_chi), & ! intent(in)
3414 : sigma_x_1_n(:,:,iiPDF_Ncn), sigma2_on_mu2_ip_1(:,:,iiPDF_Ncn), & ! intent(in)
3415 0 : corr_array_1(:,:,iiPDF_Ncn,iiPDF_chi) ) ! intent(out)
3416 :
3417 : ! Transform the correlation of eta (old t) and N_cn to standard space in PDF component 1.
3418 : call corr_NN2NL( nz, ngrdcol, & ! intent(in)
3419 : corr_array_1_n(:,:,iiPDF_Ncn,iiPDF_eta), & ! intent(in)
3420 : sigma_x_1_n(:,:,iiPDF_Ncn), sigma2_on_mu2_ip_1(:,:,iiPDF_Ncn), & ! intent(in)
3421 0 : corr_array_1(:,:,iiPDF_Ncn,iiPDF_eta) )! intent(out)
3422 :
3423 : ! Transform the correlation of w and N_cn to standard space in PDF component 1.
3424 : call corr_NN2NL( nz, ngrdcol, & ! intent(in)
3425 : corr_array_1_n(:,:,iiPDF_Ncn,iiPDF_w), & ! intent(in)
3426 : sigma_x_1_n(:,:,iiPDF_Ncn), sigma2_on_mu2_ip_1(:,:,iiPDF_Ncn), & ! intent(in)
3427 0 : corr_array_1(:,:,iiPDF_Ncn,iiPDF_w) ) ! intent(out)
3428 :
3429 : ! Transform the correlation of chi (old s) and N_cn to standard space in PDF component 2.
3430 : call corr_NN2NL( nz, ngrdcol, & ! intent(in)
3431 : corr_array_2_n(:,:,iiPDF_Ncn,iiPDF_chi), & ! intent(in)
3432 : sigma_x_2_n(:,:,iiPDF_Ncn), sigma2_on_mu2_ip_1(:,:,iiPDF_Ncn), & ! intent(in)
3433 0 : corr_array_2(:,:,iiPDF_Ncn,iiPDF_chi) ) ! intent(out)
3434 :
3435 : ! Transform the correlation of eta (old t) and N_cn to standard space in PDF component 2.
3436 : call corr_NN2NL( nz, ngrdcol, & ! intent(in)
3437 : corr_array_2_n(:,:,iiPDF_Ncn,iiPDF_eta), & ! intent(in)
3438 : sigma_x_2_n(:,:,iiPDF_Ncn), sigma2_on_mu2_ip_1(:,:,iiPDF_Ncn), & ! intent(in)
3439 0 : corr_array_2(:,:,iiPDF_Ncn,iiPDF_eta) ) ! intent(out)
3440 :
3441 : ! Transform the correlation of w and N_cn to standard space in PDF component 2.
3442 : call corr_NN2NL( nz, ngrdcol, & ! intent(in)
3443 : corr_array_2_n(:,:,iiPDF_Ncn,iiPDF_w), & ! intent(in)
3444 : sigma_x_2_n(:,:,iiPDF_Ncn), sigma2_on_mu2_ip_1(:,:,iiPDF_Ncn), & ! intent(in)
3445 0 : corr_array_2(:,:,iiPDF_Ncn,iiPDF_w) ) ! intent(out)
3446 :
3447 : ! Transform the correlations (in-precip) between chi/eta/w and the
3448 : ! precipitating hydrometeors to standard space.
3449 0 : do ivar = iiPDF_chi, iiPDF_w
3450 0 : do jvar = iiPDF_Ncn+1, pdf_dim
3451 :
3452 : ! Transform the correlation (in-precip) between w, chi, or eta and a
3453 : ! precipitating hydrometeor, hm, to standard space in PDF component 1.
3454 : call corr_NN2NL( nz, ngrdcol, & ! intent(in)
3455 : corr_array_1_n(:,:,jvar,ivar), & ! intent(in)
3456 : sigma_x_1_n(:,:,jvar), sigma2_on_mu2_ip_1(:,:,jvar), & ! intent(in)
3457 0 : corr_array_1(:,:,jvar,ivar) ) ! intent(out)
3458 :
3459 : ! Transform the correlation (in-precip) between w, chi, or eta and a
3460 : ! precipitating hydrometeor, hm, to standard space in PDF component 2.
3461 : call corr_NN2NL( nz, ngrdcol, & ! intent(in)
3462 : corr_array_2_n(:,:,jvar,ivar), & ! intent(in)
3463 : sigma_x_2_n(:,:,jvar), sigma2_on_mu2_ip_2(:,:,jvar), & ! intent(in)
3464 0 : corr_array_2(:,:,jvar,ivar) ) ! intent(out)
3465 :
3466 : end do ! jvar = iiPDF_Ncn+1, pdf_dim
3467 : end do ! ivar = iiPDF_chi, iiPDF_w
3468 :
3469 :
3470 : !!! Calculate the true correlation of two variables that both have an
3471 : !!! assumed lognormal distribution for the ith PDF component, given their
3472 : !!! normal space correlation and both of their normal space standard
3473 : !!! deviations.
3474 :
3475 : ! Transform the correlations (in-precip) between N_cn and the precipitating
3476 : ! hydrometeors to standard space.
3477 0 : do jvar = iiPDF_Ncn+1, pdf_dim
3478 :
3479 : ! Transform the correlation (in-precip) between N_cn and a precipitating
3480 : ! hydrometeor, hm, to standard space in PDF component 1.
3481 : call corr_NN2LL( nz, ngrdcol, & ! intent(in)
3482 : corr_array_1_n(:,:,jvar,iiPDF_Ncn), & ! intent(in)
3483 : sigma_x_1_n(:,:,iiPDF_Ncn), sigma_x_1_n(:,:,jvar), & ! intent(in)
3484 : sigma2_on_mu2_ip_1(:,:,iiPDF_Ncn), sigma2_on_mu2_ip_1(:,:,jvar), & ! in
3485 0 : corr_array_1(:,:,jvar,iiPDF_Ncn)) ! intent(out)
3486 :
3487 : ! Transform the correlation (in-precip) between N_cn and a precipitating
3488 : ! hydrometeor, hm, to standard space in PDF component 2.
3489 : call corr_NN2LL( nz, ngrdcol, & ! intent(in)
3490 : corr_array_2_n(:,:,jvar,iiPDF_Ncn), & ! intent(in)
3491 : sigma_x_2_n(:,:,iiPDF_Ncn), sigma_x_2_n(:,:,jvar), & ! intent(in)
3492 : sigma2_on_mu2_ip_1(:,:,iiPDF_Ncn), sigma2_on_mu2_ip_2(:,:,jvar), & ! in
3493 0 : corr_array_2(:,:,jvar,iiPDF_Ncn)) ! intent(out)
3494 :
3495 : end do ! jvar = ivar+1, pdf_dim
3496 :
3497 : ! Transform the correlations (in-precip) between two precipitating
3498 : ! hydrometeors to standard space.
3499 0 : do ivar = iiPDF_Ncn+1, pdf_dim-1
3500 0 : do jvar = ivar+1, pdf_dim
3501 :
3502 : ! Transform the correlation (in-precip) between two precipitating
3503 : ! hydrometeors (for example, r_r and N_r) to standard space in PDF
3504 : ! component 1.
3505 : call corr_NN2LL( nz, ngrdcol, & ! intent(in)
3506 : corr_array_1_n(:,:,jvar,ivar), & ! intent(in)
3507 : sigma_x_1_n(:,:,ivar), sigma_x_1_n(:,:,jvar), & ! intent(in)
3508 : sigma2_on_mu2_ip_1(:,:,ivar), sigma2_on_mu2_ip_1(:,:,jvar), & ! in
3509 0 : corr_array_1(:,:,jvar,ivar)) ! intent(out)
3510 :
3511 : ! Transform the correlation (in-precip) between two precipitating
3512 : ! hydrometeors (for example, r_r and N_r) to standard space in PDF
3513 : ! component 2.
3514 : call corr_NN2LL( nz, ngrdcol, & ! intent(in)
3515 : corr_array_2_n(:,:,jvar,ivar), & ! intent(in)
3516 : sigma_x_2_n(:,:,ivar), sigma_x_2_n(:,:,jvar), & ! intent(in)
3517 : sigma2_on_mu2_ip_2(:,:,ivar), sigma2_on_mu2_ip_2(:,:,jvar), & ! in
3518 0 : corr_array_2(:,:,jvar,ivar)) ! intent(out)
3519 :
3520 : end do ! jvar = ivar+1, pdf_dim
3521 : end do ! ivar = iiPDF_Ncn+1, pdf_dim-1
3522 :
3523 0 : return
3524 :
3525 : end subroutine denorm_transform_corr
3526 :
3527 : !=============================================================================
3528 0 : elemental subroutine calc_corr_w_hm_n( wm, wphydrometp, &
3529 : mu_w_1, mu_w_2, &
3530 : mu_hm_1, mu_hm_2, &
3531 : sigma_w_1, sigma_w_2, &
3532 : sigma_hm_1, sigma_hm_2, &
3533 : sigma_hm_1_n, sigma_hm_2_n, &
3534 : mixt_frac, precip_frac_1, &
3535 : precip_frac_2, hm_tol, &
3536 : corr_w_hm_1_n, corr_w_hm_2_n )
3537 :
3538 : ! Description:
3539 : ! Calculates the PDF component correlation (in-precip) between vertical
3540 : ! velocity, w, and the natural logarithm of a hydrometeor, ln hm. The
3541 : ! overall covariance of w and hm, <w'hm'> can be written in terms of the PDF
3542 : ! parameters. When both w and hm vary in both PDF components, the equation
3543 : ! is written as:
3544 : !
3545 : ! <w'hm'> = mixt_frac * precip_frac_1
3546 : ! * ( mu_w_1 - <w>
3547 : ! + corr_w_hm_1_n * sigma_w_1 * sigma_hm_1_n ) * mu_hm_1
3548 : ! + ( 1 - mixt_frac ) * precip_frac_2
3549 : ! * ( mu_w_2 - <w>
3550 : ! + corr_w_hm_2_n * sigma_w_2 * sigma_hm_2_n ) * mu_hm_2.
3551 : !
3552 : ! The overall covariance is provided, so the component correlation is solved
3553 : ! by setting corr_w_hm_1_n = corr_w_hm_2_n ( = corr_w_hm_n ). The equation
3554 : ! is:
3555 : !
3556 : ! corr_w_hm_n
3557 : ! = ( <w'hm'>
3558 : ! - mixt_frac * precip_frac_1 * ( mu_w_1 - <w> ) * mu_hm_1
3559 : ! - ( 1 - mixt_frac ) * precip_frac_2 * ( mu_w_2 - <w> ) * mu_hm_2 )
3560 : ! / ( mixt_frac * precip_frac_1 * sigma_w_1 * sigma_hm_1_n * mu_hm_1
3561 : ! + ( 1 - mixt_frac ) * precip_frac_2
3562 : ! * sigma_w_2 * sigma_hm_2_n * mu_hm_2 );
3563 : !
3564 : ! again, where corr_w_hm_1_n = corr_w_hm_2_n = corr_w_hm_n. When either w
3565 : ! or hm is constant in one PDF component, but both w and hm vary in the
3566 : ! other PDF component, the equation for <w'hm'> is written as:
3567 : !
3568 : ! <w'hm'> = mixt_frac * precip_frac_1
3569 : ! * ( mu_w_1 - <w>
3570 : ! + corr_w_hm_1_n * sigma_w_1 * sigma_hm_1_n ) * mu_hm_1
3571 : ! + ( 1 - mixt_frac ) * precip_frac_2
3572 : ! * ( mu_w_2 - <w> ) * mu_hm_2.
3573 : !
3574 : ! In the above equation, either w or hm (or both) is (are) constant in PDF
3575 : ! component 2, but both w and hm vary in PDF component 1. When both w and
3576 : ! hm vary in PDF component 2, but at least one of w or hm is constant in PDF
3577 : ! component 1, the equation is similar. The above equation can be rewritten
3578 : ! to solve for corr_w_hm_1_n, such that:
3579 : !
3580 : ! corr_w_hm_1_n
3581 : ! = ( <w'hm'>
3582 : ! - mixt_frac * precip_frac_1 * ( mu_w_1 - <w> ) * mu_hm_1
3583 : ! - ( 1 - mixt_frac ) * precip_frac_2 * ( mu_w_2 - <w> ) * mu_hm_2 )
3584 : ! / ( mixt_frac * precip_frac_1 * sigma_w_1 * sigma_hm_1_n * mu_hm_1 ).
3585 : !
3586 : ! Since either w or hm is constant in PDF component 2, corr_w_hm_2_n is
3587 : ! undefined. When both w and hm vary in PDF component 2, but at least one
3588 : ! of w or hm is constant in PDF component 1, the equation is similar, but
3589 : ! is in terms of corr_w_hm_2_n, while corr_w_hm_1_n is undefined. When
3590 : ! either w or hm is constant in both PDF components, the equation for
3591 : ! <w'hm'> is:
3592 : !
3593 : ! <w'hm'> = mixt_frac * precip_frac_1
3594 : ! * ( mu_w_1 - <w> ) * mu_hm_1
3595 : ! + ( 1 - mixt_frac ) * precip_frac_2
3596 : ! * ( mu_w_2 - <w> ) * mu_hm_2.
3597 : !
3598 : ! When this is the case, both corr_w_hm_1_n and corr_w_hm_2_n are undefined.
3599 :
3600 : ! References:
3601 : !-----------------------------------------------------------------------
3602 :
3603 : use constants_clubb, only: &
3604 : one, & ! Constant(s)
3605 : zero, &
3606 : max_mag_correlation, &
3607 : w_tol
3608 :
3609 : use clubb_precision, only: &
3610 : core_rknd ! Variable(s)
3611 :
3612 : implicit none
3613 :
3614 : ! Input Variables
3615 : real( kind = core_rknd ), intent(in) :: &
3616 : wm, & ! Mean vertical velocity (overall), <w> [m/s]
3617 : wphydrometp, & ! Covariance of w and hm (overall), <w'hm'> [m/s(hm un)]
3618 : mu_w_1, & ! Mean of w (1st PDF component) [m/s]
3619 : mu_w_2, & ! Mean of w (2nd PDF component) [m/s]
3620 : mu_hm_1, & ! Mean of hm (1st PDF component) in-precip (ip) [hm un]
3621 : mu_hm_2, & ! Mean of hm (2nd PDF component) ip [hm un]
3622 : sigma_w_1, & ! Standard deviation of w (1st PDF component) [m/s]
3623 : sigma_w_2, & ! Standard deviation of w (2nd PDF component) [m/s]
3624 : sigma_hm_1, & ! Standard deviation of hm (1st PDF component) ip [hm un]
3625 : sigma_hm_2, & ! Standard deviation of hm (2nd PDF component) ip [hm un]
3626 : sigma_hm_1_n, & ! Standard deviation of ln hm (1st PDF component) ip [-]
3627 : sigma_hm_2_n, & ! Standard deviation of ln hm (2nd PDF component) ip [-]
3628 : mixt_frac, & ! Mixture fraction [-]
3629 : precip_frac_1, & ! Precipitation fraction (1st PDF component) [-]
3630 : precip_frac_2 ! Precipitation fraction (2nd PDF component) [-]
3631 :
3632 : real( kind = core_rknd ), intent(in) :: &
3633 : hm_tol ! Hydrometeor tolerance value [hm un]
3634 :
3635 : ! Output Variables
3636 : real( kind = core_rknd ), intent(out) :: &
3637 : corr_w_hm_1_n, & ! Correlation of w and ln hm (1st PDF component) ip [-]
3638 : corr_w_hm_2_n ! Correlation of w and ln hm (2nd PDF component) ip [-]
3639 :
3640 : ! Local Variables
3641 : real( kind = core_rknd ) :: &
3642 : corr_w_hm_n ! Correlation of w and ln hm (both PDF components) ip [-]
3643 :
3644 :
3645 : ! Calculate the PDF component correlation of vertical velocity, w, and the
3646 : ! natural logarithm of a hydrometeor, ln hm, in precipitation.
3647 : if ( sigma_w_1 > w_tol .and. sigma_hm_1 > hm_tol .and. &
3648 0 : sigma_w_2 > w_tol .and. sigma_hm_2 > hm_tol ) then
3649 :
3650 : ! Both w and hm vary in both PDF components.
3651 : ! Calculate corr_w_hm_n (where corr_w_hm_1_n = corr_w_hm_2_n
3652 : ! = corr_w_hm_n).
3653 : corr_w_hm_n &
3654 : = ( wphydrometp &
3655 : - mixt_frac * precip_frac_1 * ( mu_w_1 - wm ) * mu_hm_1 &
3656 : - ( one - mixt_frac ) * precip_frac_2 * ( mu_w_2 - wm ) * mu_hm_2 ) &
3657 : / ( mixt_frac * precip_frac_1 * sigma_w_1 * sigma_hm_1_n * mu_hm_1 &
3658 : + ( one - mixt_frac ) * precip_frac_2 &
3659 0 : * sigma_w_2 * sigma_hm_2_n * mu_hm_2 )
3660 :
3661 : ! Check that the PDF component correlations have reasonable values.
3662 0 : if ( corr_w_hm_n > max_mag_correlation ) then
3663 : corr_w_hm_n = max_mag_correlation
3664 0 : elseif ( corr_w_hm_n < -max_mag_correlation ) then
3665 0 : corr_w_hm_n = -max_mag_correlation
3666 : end if
3667 :
3668 : ! The PDF component correlations between w and ln hm (in-precip) are
3669 : ! equal.
3670 0 : corr_w_hm_1_n = corr_w_hm_n
3671 0 : corr_w_hm_2_n = corr_w_hm_n
3672 :
3673 :
3674 0 : elseif ( sigma_w_1 > w_tol .and. sigma_hm_1 > hm_tol ) then
3675 :
3676 : ! Both w and hm vary in PDF component 1, but at least one of w and hm is
3677 : ! constant in PDF component 2.
3678 : ! Calculate the PDF component 1 correlation of w and ln hm (in-precip).
3679 : corr_w_hm_1_n &
3680 : = ( wphydrometp &
3681 : - mixt_frac * precip_frac_1 * ( mu_w_1 - wm ) * mu_hm_1 &
3682 : - ( one - mixt_frac ) * precip_frac_2 * ( mu_w_2 - wm ) * mu_hm_2 ) &
3683 0 : / ( mixt_frac * precip_frac_1 * sigma_w_1 * sigma_hm_1_n * mu_hm_1 )
3684 :
3685 : ! Check that the PDF component 1 correlation has a reasonable value.
3686 0 : if ( corr_w_hm_1_n > max_mag_correlation ) then
3687 0 : corr_w_hm_1_n = max_mag_correlation
3688 0 : elseif ( corr_w_hm_1_n < -max_mag_correlation ) then
3689 0 : corr_w_hm_1_n = -max_mag_correlation
3690 : end if
3691 :
3692 : ! The PDF component 2 correlation is undefined.
3693 0 : corr_w_hm_2_n = zero
3694 :
3695 :
3696 0 : elseif ( sigma_w_2 > w_tol .and. sigma_hm_2 > hm_tol ) then
3697 :
3698 : ! Both w and hm vary in PDF component 2, but at least one of w and hm is
3699 : ! constant in PDF component 1.
3700 : ! Calculate the PDF component 2 correlation of w and ln hm (in-precip).
3701 : corr_w_hm_2_n &
3702 : = ( wphydrometp &
3703 : - mixt_frac * precip_frac_1 * ( mu_w_1 - wm ) * mu_hm_1 &
3704 : - ( one - mixt_frac ) * precip_frac_2 * ( mu_w_2 - wm ) * mu_hm_2 ) &
3705 : / ( ( one - mixt_frac ) * precip_frac_2 &
3706 0 : * sigma_w_2 * sigma_hm_2_n * mu_hm_2 )
3707 :
3708 : ! Check that the PDF component 2 correlation has a reasonable value.
3709 0 : if ( corr_w_hm_2_n > max_mag_correlation ) then
3710 0 : corr_w_hm_2_n = max_mag_correlation
3711 0 : elseif ( corr_w_hm_2_n < -max_mag_correlation ) then
3712 0 : corr_w_hm_2_n = -max_mag_correlation
3713 : end if
3714 :
3715 : ! The PDF component 1 correlation is undefined.
3716 0 : corr_w_hm_1_n = zero
3717 :
3718 :
3719 : else ! sigma_w_1 * sigma_hm_1 = 0 .and. sigma_w_2 * sigma_hm_2 = 0.
3720 :
3721 : ! At least one of w and hm is constant in both PDF components.
3722 :
3723 : ! The PDF component 1 and component 2 correlations are both undefined.
3724 0 : corr_w_hm_1_n = zero
3725 0 : corr_w_hm_2_n = zero
3726 :
3727 :
3728 : end if
3729 :
3730 :
3731 0 : return
3732 :
3733 : end subroutine calc_corr_w_hm_n
3734 :
3735 : !=============================================================================
3736 0 : subroutine pdf_param_hm_stats( nz, pdf_dim, hydromet_dim, hm_metadata, &
3737 0 : hm_1, hm_2, &
3738 0 : mu_x_1, mu_x_2, &
3739 0 : sigma_x_1, sigma_x_2, &
3740 0 : corr_array_1, corr_array_2, &
3741 : stats_metadata, &
3742 : stats_zt )
3743 :
3744 : ! Description:
3745 : ! Record statistics for standard PDF parameters involving hydrometeors.
3746 :
3747 : ! References:
3748 : !-----------------------------------------------------------------------
3749 :
3750 : use index_mapping, only: &
3751 : pdf2hydromet_idx ! Procedure(s)
3752 :
3753 : use corr_varnce_module, only: &
3754 : hm_metadata_type
3755 :
3756 : use clubb_precision, only: &
3757 : core_rknd ! Variable(s)
3758 :
3759 : use stats_type_utilities, only: &
3760 : stat_update_var_pt ! Procedure(s)
3761 :
3762 : use stats_variables, only: &
3763 : stats_metadata_type
3764 :
3765 : use stats_type, only: stats ! Type
3766 :
3767 : implicit none
3768 :
3769 : !--------------------------- Input Variables ---------------------------
3770 : integer, intent(in) :: &
3771 : nz, & ! Number of vertical levels
3772 : pdf_dim, & ! Number of variables in the correlation array
3773 : hydromet_dim ! Number of hydrometeor species
3774 :
3775 : type (hm_metadata_type), intent(in) :: &
3776 : hm_metadata
3777 :
3778 : real( kind = core_rknd ), dimension(nz,hydromet_dim), intent(in) :: &
3779 : hm_1, & ! Mean of a precip. hydrometeor (1st PDF component) [units vary]
3780 : hm_2 ! Mean of a precip. hydrometeor (2nd PDF component) [units vary]
3781 :
3782 : real( kind = core_rknd ), dimension(nz,pdf_dim), intent(in) :: &
3783 : mu_x_1, & ! Mean array of PDF vars. (1st PDF component) [units vary]
3784 : mu_x_2, & ! Mean array of PDF vars. (2nd PDF component) [units vary]
3785 : sigma_x_1, & ! Standard deviation array of PDF vars (comp. 1) [units vary]
3786 : sigma_x_2 ! Standard deviation array of PDF vars (comp. 2) [units vary]
3787 :
3788 : real( kind = core_rknd ), dimension(nz,pdf_dim,pdf_dim), intent(in) :: &
3789 : corr_array_1, & ! Correlation array of PDF vars. (comp. 1) [-]
3790 : corr_array_2 ! Correlation array of PDF vars. (comp. 2) [-]
3791 :
3792 : type (stats_metadata_type), intent(in) :: &
3793 : stats_metadata
3794 :
3795 : !--------------------------- InOut Variable ---------------------------
3796 : type (stats), target, intent(inout) :: &
3797 : stats_zt
3798 :
3799 : !--------------------------- Local Variable ---------------------------
3800 : integer :: ivar, jvar, hm_idx, hm_idx_ivar, hm_idx_jvar, k ! Indices
3801 :
3802 : ! Just to avoid typing hm_metadata%iiPDF_x everywhere
3803 : integer :: &
3804 : iiPDF_chi, &
3805 : iiPDF_eta, &
3806 : iiPDF_w, &
3807 : iiPDF_Ncn
3808 :
3809 : !------------------------ Begin Code ------------------------
3810 :
3811 0 : iiPDF_chi = hm_metadata%iiPDF_chi
3812 0 : iiPDF_eta = hm_metadata%iiPDF_eta
3813 0 : iiPDF_w = hm_metadata%iiPDF_w
3814 0 : iiPDF_Ncn = hm_metadata%iiPDF_Ncn
3815 :
3816 : !!! Output the statistics for hydrometeor PDF parameters.
3817 :
3818 : ! Statistics
3819 0 : if ( stats_metadata%l_stats_samp ) then
3820 :
3821 : ! Switch back to using stat_update_var once the code is generalized
3822 : ! to pass in the number of vertical levels.
3823 0 : do ivar = 1, hydromet_dim, 1
3824 :
3825 0 : if ( stats_metadata%ihm_1(ivar) > 0 ) then
3826 : ! Mean of the precipitating hydrometeor in PDF component 1.
3827 : ! call stat_update_var( stats_metadata%ihm_1(ivar), hm_1(:,ivar), stats_zt )
3828 0 : do k = 1, nz, 1
3829 0 : call stat_update_var_pt( stats_metadata%ihm_1(ivar), k, hm_1(k,ivar), & ! intent(in)
3830 0 : stats_zt ) ! intent(inout)
3831 : end do ! k = 1, nz, 1
3832 : end if
3833 :
3834 0 : if ( stats_metadata%ihm_2(ivar) > 0 ) then
3835 : ! Mean of the precipitating hydrometeor in PDF component 2.
3836 : ! call stat_update_var( stats_metadata%ihm_2(ivar), hm_2(:,ivar), stats_zt )
3837 0 : do k = 1, nz, 1
3838 0 : call stat_update_var_pt( stats_metadata%ihm_2(ivar), k, hm_2(k,ivar), & ! intent(in)
3839 0 : stats_zt ) ! intent(inout)
3840 : end do ! k = 1, nz, 1
3841 : end if
3842 :
3843 : end do ! ivar = 1, hydromet_dim, 1
3844 :
3845 0 : do ivar = iiPDF_Ncn+1, pdf_dim, 1
3846 :
3847 0 : hm_idx = pdf2hydromet_idx(ivar,hm_metadata)
3848 :
3849 : ! Mean of the precipitating hydrometeor (in-precip)
3850 : ! in PDF component 1.
3851 0 : if ( stats_metadata%imu_hm_1(hm_idx) > 0 ) then
3852 : ! call stat_update_var( stats_metadata%imu_hm_1(hm_idx), mu_x_1(ivar,:), stats_zt )
3853 0 : do k = 1, nz, 1
3854 0 : call stat_update_var_pt( stats_metadata%imu_hm_1(hm_idx), k, mu_x_1(k,ivar), & ! intent(in)
3855 0 : stats_zt ) ! intent(inout)
3856 : end do ! k = 1, nz, 1
3857 : end if
3858 :
3859 : ! Mean of the precipitating hydrometeor (in-precip)
3860 : ! in PDF component 2.
3861 0 : if ( stats_metadata%imu_hm_2(hm_idx) > 0 ) then
3862 : ! call stat_update_var( stats_metadata%imu_hm_2(hm_idx), mu_x_2(ivar,:), stats_zt )
3863 0 : do k = 1, nz, 1
3864 0 : call stat_update_var_pt( stats_metadata%imu_hm_2(hm_idx), k, mu_x_2(k,ivar), & ! intent(in)
3865 0 : stats_zt ) ! intent(inout)
3866 : end do ! k = 1, nz, 1
3867 : end if
3868 :
3869 : end do ! ivar = iiPDF_Ncn+1, pdf_dim, 1
3870 :
3871 : ! Mean of cloud nuclei concentration in PDF component 1.
3872 0 : if ( stats_metadata%imu_Ncn_1 > 0 ) then
3873 : ! call stat_update_var( stats_metadata%imu_Ncn_1, mu_x_1(iiPDF_Ncn,:), stats_zt )
3874 0 : do k = 1, nz, 1
3875 0 : call stat_update_var_pt( stats_metadata%imu_Ncn_1, k, mu_x_1(k,iiPDF_Ncn), & ! intent(in)
3876 0 : stats_zt ) ! intent(inout)
3877 : end do ! k = 1, nz, 1
3878 : end if
3879 :
3880 : ! Mean of cloud nuclei concentration in PDF component 2.
3881 0 : if ( stats_metadata%imu_Ncn_2 > 0 ) then
3882 : ! call stat_update_var( stats_metadata%imu_Ncn_2, mu_x_2(iiPDF_Ncn,:), stats_zt )
3883 0 : do k = 1, nz, 1
3884 0 : call stat_update_var_pt( stats_metadata%imu_Ncn_2, k, mu_x_2(k,iiPDF_Ncn), & ! intent(in)
3885 0 : stats_zt ) ! intent(inout)
3886 : end do ! k = 1, nz, 1
3887 : end if
3888 :
3889 0 : do ivar = iiPDF_Ncn+1, pdf_dim, 1
3890 :
3891 0 : hm_idx = pdf2hydromet_idx(ivar,hm_metadata)
3892 :
3893 : ! Standard deviation of the precipitating hydrometeor (in-precip)
3894 : ! in PDF component 1.
3895 0 : if ( stats_metadata%isigma_hm_1(hm_idx) > 0 ) then
3896 : ! call stat_update_var( stats_metadata%isigma_hm_1(hm_idx), sigma_x_1(ivar,:), &
3897 : ! stats_zt )
3898 0 : do k = 1, nz, 1
3899 0 : call stat_update_var_pt( stats_metadata%isigma_hm_1(hm_idx), k, sigma_x_1(k,ivar), & ! intent(in)
3900 0 : stats_zt ) ! intent(inout)
3901 : end do ! k = 1, nz, 1
3902 : end if
3903 :
3904 : ! Standard deviation of the precipitating hydrometeor (in-precip)
3905 : ! in PDF component 2.
3906 0 : if ( stats_metadata%isigma_hm_2(hm_idx) > 0 ) then
3907 : ! call stat_update_var( stats_metadata%isigma_hm_2(hm_idx), sigma_x_2(ivar,:), &
3908 : ! stats_zt )
3909 0 : do k = 1, nz, 1
3910 0 : call stat_update_var_pt( stats_metadata%isigma_hm_2(hm_idx), k, sigma_x_2(k,ivar), & ! intent(in)
3911 0 : stats_zt ) ! intent(inout)
3912 : end do ! k = 1, nz, 1
3913 : end if
3914 :
3915 : end do ! ivar = iiPDF_Ncn+1, pdf_dim, 1
3916 :
3917 : ! Standard deviation of cloud nuclei concentration in PDF component 1.
3918 0 : if ( stats_metadata%isigma_Ncn_1 > 0 ) then
3919 : ! call stat_update_var( stats_metadata%isigma_Ncn_1, sigma_x_1(iiPDF_Ncn,:), stats_zt )
3920 0 : do k = 1, nz, 1
3921 0 : call stat_update_var_pt( stats_metadata%isigma_Ncn_1, k, sigma_x_1(k,iiPDF_Ncn), & ! intent(in)
3922 0 : stats_zt ) ! intent(inout)
3923 : end do ! k = 1, nz, 1
3924 : end if
3925 :
3926 : ! Standard deviation of cloud nuclei concentration in PDF component 2.
3927 0 : if ( stats_metadata%isigma_Ncn_2 > 0 ) then
3928 : ! call stat_update_var( stats_metadata%isigma_Ncn_2, sigma_x_2(iiPDF_Ncn,:), stats_zt )
3929 0 : do k = 1, nz, 1
3930 0 : call stat_update_var_pt( stats_metadata%isigma_Ncn_2, k, sigma_x_2(k,iiPDF_Ncn), & ! intent(in)
3931 0 : stats_zt ) ! intent(inout)
3932 : end do ! k = 1, nz, 1
3933 : end if
3934 :
3935 : ! Correlation of w and chi (old s) in PDF component 1 found in the
3936 : ! correlation array.
3937 : ! The true correlation of w and chi in each PDF component is solved for
3938 : ! by an equation and is part of CLUBB's PDF parameters. However, there
3939 : ! is an option in CLUBB, l_fix_w_chi_eta_correlations, that sets the
3940 : ! component correlation of w and chi to a constant, prescribed value
3941 : ! because of SILHS. The correlation of w and chi in PDF component 1
3942 : ! that is calculated by an equation is stored in stats as "corr_w_chi_1".
3943 : ! Here, "corr_w_chi_1_ca" outputs whatever value is found in the
3944 : ! correlation array, whether or not it matches "corr_w_chi_1".
3945 0 : if ( stats_metadata%icorr_w_chi_1_ca > 0 ) then
3946 : ! call stat_update_var( stats_metadata%icorr_w_chi_1_ca, &
3947 : ! corr_array_1(iiPDF_w,iiPDF_chi,:), stats_zt )
3948 0 : do k = 1, nz, 1
3949 : call stat_update_var_pt( stats_metadata%icorr_w_chi_1_ca, k, &
3950 0 : corr_array_1(k,iiPDF_w,iiPDF_chi), & ! intent(in)
3951 0 : stats_zt ) ! intent(inout)
3952 : end do ! k = 1, nz, 1
3953 : end if
3954 :
3955 : ! Correlation of w and chi (old s) in PDF component 2 found in the
3956 : ! correlation array.
3957 : ! The true correlation of w and chi in each PDF component is solved for
3958 : ! by an equation and is part of CLUBB's PDF parameters. However, there
3959 : ! is an option in CLUBB, l_fix_w_chi_eta_correlations, that sets the
3960 : ! component correlation of w and chi to a constant, prescribed value
3961 : ! because of SILHS. The correlation of w and chi in PDF component 2
3962 : ! that is calculated by an equation is stored in stats as "corr_w_chi_2".
3963 : ! Here, "corr_w_chi_2_ca" outputs whatever value is found in the
3964 : ! correlation array, whether or not it matches "corr_w_chi_2".
3965 0 : if ( stats_metadata%icorr_w_chi_2_ca > 0 ) then
3966 : ! call stat_update_var( stats_metadata%icorr_w_chi_2_ca, &
3967 : ! corr_array_2(iiPDF_w,iiPDF_chi,:), stats_zt )
3968 0 : do k = 1, nz, 1
3969 : call stat_update_var_pt( stats_metadata%icorr_w_chi_2_ca, k, & ! intent(in)
3970 0 : corr_array_2(k,iiPDF_w,iiPDF_chi), & ! intent(in)
3971 0 : stats_zt ) ! intent(inout)
3972 : end do ! k = 1, nz, 1
3973 : end if
3974 :
3975 : ! Correlation of w and eta (old t) in PDF component 1 found in the
3976 : ! correlation array.
3977 : ! The true correlation of w and eta in each PDF component is solved for
3978 : ! by an equation and is part of CLUBB's PDF parameters. However, there
3979 : ! is an option in CLUBB, l_fix_w_chi_eta_correlations, that sets the
3980 : ! component correlation of w and eta to a constant, prescribed value
3981 : ! because of SILHS. The correlation of w and eta in PDF component 1
3982 : ! that is calculated by an equation is stored in stats as "corr_w_eta_1".
3983 : ! Here, "corr_w_eta_1_ca" outputs whatever value is found in the
3984 : ! correlation array, whether or not it matches "corr_w_eta_1".
3985 0 : if ( stats_metadata%icorr_w_eta_1_ca > 0 ) then
3986 : ! call stat_update_var( stats_metadata%icorr_w_eta_1_ca, &
3987 : ! corr_array_1(iiPDF_w,iiPDF_eta,:), stats_zt )
3988 0 : do k = 1, nz, 1
3989 : call stat_update_var_pt( stats_metadata%icorr_w_eta_1_ca, k, & ! intent(in)
3990 0 : corr_array_1(k,iiPDF_w,iiPDF_eta), & ! intent(in)
3991 0 : stats_zt ) ! intent(inout)
3992 : end do ! k = 1, nz, 1
3993 : end if
3994 :
3995 : ! Correlation of w and eta (old t) in PDF component 2 found in the
3996 : ! correlation array.
3997 : ! The true correlation of w and eta in each PDF component is solved for
3998 : ! by an equation and is part of CLUBB's PDF parameters. However, there
3999 : ! is an option in CLUBB, l_fix_w_chi_eta_correlations, that sets the
4000 : ! component correlation of w and eta to a constant, prescribed value
4001 : ! because of SILHS. The correlation of w and eta in PDF component 2
4002 : ! that is calculated by an equation is stored in stats as "corr_w_eta_2".
4003 : ! Here, "corr_w_eta_2_ca" outputs whatever value is found in the
4004 : ! correlation array, whether or not it matches "corr_w_eta_2".
4005 0 : if ( stats_metadata%icorr_w_eta_2_ca > 0 ) then
4006 : ! call stat_update_var( stats_metadata%icorr_w_eta_2_ca, &
4007 : ! corr_array_2(iiPDF_w,iiPDF_eta,:), stats_zt )
4008 0 : do k = 1, nz, 1
4009 : call stat_update_var_pt( stats_metadata%icorr_w_eta_2_ca, k, & ! intent(in)
4010 0 : corr_array_2(k,iiPDF_w,iiPDF_eta), & ! intent(in)
4011 0 : stats_zt ) ! intent(inout)
4012 : end do ! k = 1, nz, 1
4013 : end if
4014 :
4015 0 : do ivar = iiPDF_Ncn+1, pdf_dim, 1
4016 :
4017 0 : hm_idx = pdf2hydromet_idx(ivar,hm_metadata)
4018 :
4019 : ! Correlation (in-precip) of w and the precipitating hydrometeor
4020 : ! in PDF component 1.
4021 0 : if ( stats_metadata%icorr_w_hm_1(hm_idx) > 0 ) then
4022 : ! call stat_update_var( stats_metadata%icorr_w_hm_1(hm_idx), &
4023 : ! corr_array_1(ivar,iiPDF_w,:), stats_zt )
4024 0 : do k = 1, nz, 1
4025 : call stat_update_var_pt( stats_metadata%icorr_w_hm_1(hm_idx), k, & ! intent(in)
4026 0 : corr_array_1(k,ivar,iiPDF_w), & ! intent(in)
4027 0 : stats_zt ) ! intent(inout)
4028 : end do ! k = 1, nz, 1
4029 : end if
4030 :
4031 : ! Correlation (in-precip) of w and the precipitating hydrometeor
4032 : ! in PDF component 2.
4033 0 : if ( stats_metadata%icorr_w_hm_2(hm_idx) > 0 ) then
4034 : ! call stat_update_var( stats_metadata%icorr_w_hm_2(hm_idx), &
4035 : ! corr_array_2(ivar,iiPDF_w,:), stats_zt )
4036 0 : do k = 1, nz, 1
4037 : call stat_update_var_pt( stats_metadata%icorr_w_hm_2(hm_idx), k, & ! intent(in)
4038 0 : corr_array_2(k,ivar,iiPDF_w), & ! intent(in)
4039 0 : stats_zt ) ! intent(inout)
4040 : end do ! k = 1, nz, 1
4041 : end if
4042 :
4043 : end do ! ivar = iiPDF_Ncn+1, pdf_dim, 1
4044 :
4045 : ! Correlation of w and N_cn in PDF component 1.
4046 0 : if ( stats_metadata%icorr_w_Ncn_1 > 0 ) then
4047 : ! call stat_update_var( stats_metadata%icorr_w_Ncn_1, &
4048 : ! corr_array_1(iiPDF_Ncn,iiPDF_w,:), stats_zt )
4049 0 : do k = 1, nz, 1
4050 : call stat_update_var_pt( stats_metadata%icorr_w_Ncn_1, k, & ! intent(in)
4051 0 : corr_array_1(k,iiPDF_Ncn,iiPDF_w), & ! intent(in)
4052 0 : stats_zt ) ! intent(inout)
4053 : end do ! k = 1, nz, 1
4054 : end if
4055 :
4056 : ! Correlation of w and N_cn in PDF component 2.
4057 0 : if ( stats_metadata%icorr_w_Ncn_2 > 0 ) then
4058 : ! call stat_update_var( stats_metadata%icorr_w_Ncn_2, &
4059 : ! corr_array_2(iiPDF_Ncn,iiPDF_w,:), stats_zt )
4060 0 : do k = 1, nz, 1
4061 : call stat_update_var_pt( stats_metadata%icorr_w_Ncn_2, k, & ! intent(in)
4062 0 : corr_array_2(k,iiPDF_Ncn,iiPDF_w), & ! intent(in)
4063 0 : stats_zt ) ! intent(inout)
4064 : end do ! k = 1, nz, 1
4065 : end if
4066 :
4067 : ! Correlation of chi (old s) and eta (old t) in PDF component 1 found in
4068 : ! the correlation array.
4069 : ! The true correlation of chi and eta in each PDF component is solved for
4070 : ! by an equation and is part of CLUBB's PDF parameters. However, there
4071 : ! is an option in CLUBB, l_fix_w_chi_eta_correlations, that sets the
4072 : ! component correlation of chi and eta to a constant, prescribed value
4073 : ! because of SILHS. The correlation of chi and eta in PDF component 1
4074 : ! that is calculated by an equation is stored in stats as
4075 : ! "corr_chi_eta_1". Here, "corr_chi_eta_1_ca" outputs whatever value is
4076 : ! found in the correlation array, whether or not it matches
4077 : ! "corr_chi_eta_1".
4078 0 : if ( stats_metadata%icorr_chi_eta_1_ca > 0 ) then
4079 : ! call stat_update_var( stats_metadata%icorr_chi_eta_1_ca, &
4080 : ! corr_array_1(iiPDF_eta,iiPDF_chi,:), stats_zt )
4081 0 : do k = 1, nz, 1
4082 : call stat_update_var_pt( stats_metadata%icorr_chi_eta_1_ca, k, & ! intent(in)
4083 0 : corr_array_1(k,iiPDF_eta,iiPDF_chi), & ! intent(in)
4084 0 : stats_zt ) ! intent(inout)
4085 : end do ! k = 1, nz, 1
4086 : end if
4087 :
4088 : ! Correlation of chi (old s) and eta (old t) in PDF component 2 found in
4089 : ! the correlation array.
4090 : ! The true correlation of chi and eta in each PDF component is solved for
4091 : ! by an equation and is part of CLUBB's PDF parameters. However, there
4092 : ! is an option in CLUBB, l_fix_w_chi_eta_correlations, that sets the
4093 : ! component correlation of chi and eta to a constant, prescribed value
4094 : ! because of SILHS. The correlation of chi and eta in PDF component 2
4095 : ! that is calculated by an equation is stored in stats as
4096 : ! "corr_chi_eta_2". Here, "corr_chi_eta_2_ca" outputs whatever value is
4097 : ! found in the correlation array, whether or not it matches
4098 : ! "corr_chi_eta_2".
4099 0 : if ( stats_metadata%icorr_chi_eta_2_ca > 0 ) then
4100 : ! call stat_update_var( stats_metadata%icorr_chi_eta_2_ca, &
4101 : ! corr_array_2(iiPDF_eta,iiPDF_chi,:), stats_zt )
4102 0 : do k = 1, nz, 1
4103 : call stat_update_var_pt( stats_metadata%icorr_chi_eta_2_ca, k, & ! intent(in)
4104 0 : corr_array_2(k,iiPDF_eta,iiPDF_chi), & ! intent(in)
4105 0 : stats_zt ) ! intent(inout)
4106 : end do ! k = 1, nz, 1
4107 : end if
4108 :
4109 0 : do ivar = iiPDF_Ncn+1, pdf_dim, 1
4110 :
4111 0 : hm_idx = pdf2hydromet_idx(ivar,hm_metadata)
4112 :
4113 : ! Correlation (in-precip) of chi (old s) and the precipitating
4114 : ! hydrometeor in PDF component 1.
4115 0 : if ( stats_metadata%icorr_chi_hm_1(hm_idx) > 0 ) then
4116 : ! call stat_update_var( stats_metadata%icorr_chi_hm_1(hm_idx), &
4117 : ! corr_array_1(ivar,iiPDF_chi,:), stats_zt )
4118 0 : do k = 1, nz, 1
4119 : call stat_update_var_pt( stats_metadata%icorr_chi_hm_1(hm_idx), k, & ! intent(in)
4120 0 : corr_array_1(k,ivar,iiPDF_chi), & ! intent(in)
4121 0 : stats_zt ) ! intent(inout)
4122 : end do ! k = 1, nz, 1
4123 : end if
4124 :
4125 : ! Correlation (in-precip) of chi (old s) and the precipitating
4126 : ! hydrometeor in PDF component 2.
4127 0 : if ( stats_metadata%icorr_chi_hm_2(hm_idx) > 0 ) then
4128 : ! call stat_update_var( stats_metadata%icorr_chi_hm_2(hm_idx), &
4129 : ! corr_array_2(ivar,iiPDF_chi,:), stats_zt )
4130 0 : do k = 1, nz, 1
4131 : call stat_update_var_pt( stats_metadata%icorr_chi_hm_2(hm_idx), k, & ! intent(in)
4132 0 : corr_array_2(k,ivar,iiPDF_chi), & ! intent(in)
4133 0 : stats_zt ) ! intent(inout)
4134 : end do ! k = 1, nz, 1
4135 : end if
4136 :
4137 : end do ! ivar = iiPDF_Ncn+1, pdf_dim, 1
4138 :
4139 : ! Correlation of chi (old s) and N_cn in PDF component 1.
4140 0 : if ( stats_metadata%icorr_chi_Ncn_1 > 0 ) then
4141 : ! call stat_update_var( stats_metadata%icorr_chi_Ncn_1, &
4142 : ! corr_array_1(iiPDF_Ncn,iiPDF_chi,:), stats_zt )
4143 0 : do k = 1, nz, 1
4144 : call stat_update_var_pt( stats_metadata%icorr_chi_Ncn_1, k, & ! intent(in)
4145 0 : corr_array_1(k,iiPDF_Ncn,iiPDF_chi), & ! intent(inout)
4146 0 : stats_zt ) ! intent(in)
4147 : end do ! k = 1, nz, 1
4148 : end if
4149 :
4150 : ! Correlation of chi (old s) and N_cn in PDF component 2.
4151 0 : if ( stats_metadata%icorr_chi_Ncn_2 > 0 ) then
4152 : ! call stat_update_var( stats_metadata%icorr_chi_Ncn_2, &
4153 : ! corr_array_2(iiPDF_Ncn,iiPDF_chi,:), stats_zt )
4154 0 : do k = 1, nz, 1
4155 : call stat_update_var_pt( stats_metadata%icorr_chi_Ncn_2, k, & ! intent(in)
4156 0 : corr_array_2(k,iiPDF_Ncn,iiPDF_chi), & ! intent(in)
4157 0 : stats_zt ) ! intent(inout)
4158 : end do ! k = 1, nz, 1
4159 : end if
4160 :
4161 0 : do ivar = iiPDF_Ncn+1, pdf_dim, 1
4162 :
4163 0 : hm_idx = pdf2hydromet_idx(ivar,hm_metadata)
4164 :
4165 : ! Correlation (in-precip) of eta (old t) and the precipitating
4166 : ! hydrometeor in PDF component 1.
4167 0 : if ( stats_metadata%icorr_eta_hm_1(hm_idx) > 0 ) then
4168 : ! call stat_update_var( stats_metadata%icorr_eta_hm_1(hm_idx), &
4169 : ! corr_array_1(ivar,iiPDF_eta,:), stats_zt )
4170 0 : do k = 1, nz, 1
4171 : call stat_update_var_pt( stats_metadata%icorr_eta_hm_1(hm_idx), k, & ! intent(in)
4172 0 : corr_array_1(k,ivar,iiPDF_eta), & ! intent(in)
4173 0 : stats_zt ) ! intent(inout)
4174 : end do ! k = 1, nz, 1
4175 : end if
4176 :
4177 : ! Correlation (in-precip) of eta (old t) and the precipitating
4178 : ! hydrometeor in PDF component 2.
4179 0 : if ( stats_metadata%icorr_eta_hm_2(hm_idx) > 0 ) then
4180 : ! call stat_update_var( stats_metadata%icorr_eta_hm_2(hm_idx), &
4181 : ! corr_array_2(ivar,iiPDF_eta,:), stats_zt )
4182 0 : do k = 1, nz, 1
4183 : call stat_update_var_pt( stats_metadata%icorr_eta_hm_2(hm_idx), k, & ! intent(in)
4184 0 : corr_array_2(k,ivar,iiPDF_eta), & ! intent(in)
4185 0 : stats_zt ) ! intent(inout)
4186 : end do ! k = 1, nz, 1
4187 : end if
4188 :
4189 : end do ! ivar = iiPDF_Ncn+1, pdf_dim, 1
4190 :
4191 : ! Correlation of eta (old t) and N_cn in PDF component 1.
4192 0 : if ( stats_metadata%icorr_eta_Ncn_1 > 0 ) then
4193 : ! call stat_update_var( stats_metadata%icorr_eta_Ncn_1, &
4194 : ! corr_array_1(iiPDF_Ncn,iiPDF_eta,:), stats_zt )
4195 0 : do k = 1, nz, 1
4196 : call stat_update_var_pt( stats_metadata%icorr_eta_Ncn_1, k, & ! intent(in)
4197 0 : corr_array_1(k,iiPDF_Ncn,iiPDF_eta), & ! intent(in)
4198 0 : stats_zt ) ! intent(inout)
4199 : end do ! k = 1, nz, 1
4200 : end if
4201 :
4202 : ! Correlation of eta (old t) and N_cn in PDF component 2.
4203 0 : if ( stats_metadata%icorr_eta_Ncn_2 > 0 ) then
4204 : ! call stat_update_var( stats_metadata%icorr_eta_Ncn_2, &
4205 : ! corr_array_2(iiPDF_Ncn,iiPDF_eta,:), stats_zt )
4206 0 : do k = 1, nz, 1
4207 : call stat_update_var_pt( stats_metadata%icorr_eta_Ncn_2, k, & ! intent(in)
4208 0 : corr_array_2(k,iiPDF_Ncn,iiPDF_eta), & ! intent(in)
4209 0 : stats_zt ) ! intent(inout)
4210 : end do ! k = 1, nz, 1
4211 : end if
4212 :
4213 0 : do ivar = iiPDF_Ncn+1, pdf_dim, 1
4214 :
4215 0 : hm_idx = pdf2hydromet_idx(ivar,hm_metadata)
4216 :
4217 : ! Correlation (in-precip) of N_cn and the precipitating
4218 : ! hydrometeor in PDF component 1.
4219 0 : if ( stats_metadata%icorr_Ncn_hm_1(hm_idx) > 0 ) then
4220 : ! call stat_update_var( stats_metadata%icorr_Ncn_hm_1(hm_idx), &
4221 : ! corr_array_1(ivar,iiPDF_Ncn,:), stats_zt )
4222 0 : do k = 1, nz, 1
4223 : call stat_update_var_pt( stats_metadata%icorr_Ncn_hm_1(hm_idx), k, & ! intent(in)
4224 0 : corr_array_1(k,ivar,iiPDF_Ncn), & ! intent(in)
4225 0 : stats_zt ) ! intent(inout)
4226 : end do ! k = 1, nz, 1
4227 : end if
4228 :
4229 : ! Correlation (in-precip) of N_cn and the precipitating
4230 : ! hydrometeor in PDF component 2.
4231 0 : if ( stats_metadata%icorr_Ncn_hm_2(hm_idx) > 0 ) then
4232 : ! call stat_update_var( stats_metadata%icorr_Ncn_hm_2(hm_idx), &
4233 : ! corr_array_2(ivar,iiPDF_Ncn,:), stats_zt )
4234 0 : do k = 1, nz, 1
4235 : call stat_update_var_pt( stats_metadata%icorr_Ncn_hm_2(hm_idx), k, & ! intent(in)
4236 0 : corr_array_2(k,ivar,iiPDF_Ncn), & ! intent(in)
4237 0 : stats_zt ) ! intent(inout)
4238 : end do ! k = 1, nz, 1
4239 : end if
4240 :
4241 : end do ! ivar = iiPDF_Ncn+1, pdf_dim, 1
4242 :
4243 0 : do ivar = iiPDF_Ncn+1, pdf_dim, 1
4244 :
4245 0 : hm_idx_ivar = pdf2hydromet_idx(ivar,hm_metadata)
4246 :
4247 0 : do jvar = ivar+1, pdf_dim, 1
4248 :
4249 0 : hm_idx_jvar = pdf2hydromet_idx(jvar,hm_metadata)
4250 :
4251 : ! Correlation (in-precip) of two different hydrometeors (hmx and
4252 : ! hmy) in PDF component 1.
4253 0 : if ( stats_metadata%icorr_hmx_hmy_1(hm_idx_jvar,hm_idx_ivar) > 0 ) then
4254 : ! call stat_update_var( stats_metadata%icorr_hmx_hmy_1(hm_idx_jvar,hm_idx_ivar), &
4255 : ! corr_array_1(jvar,ivar,:), stats_zt )
4256 0 : do k = 1, nz, 1
4257 : call stat_update_var_pt( stats_metadata%icorr_hmx_hmy_1(hm_idx_jvar,hm_idx_ivar), k, & ! in
4258 0 : corr_array_1(k,jvar,ivar), & ! intent(in)
4259 0 : stats_zt ) ! intent(inout)
4260 : end do ! k = 1, nz, 1
4261 : end if
4262 :
4263 : ! Correlation (in-precip) of two different hydrometeors (hmx and
4264 : ! hmy) in PDF component 2.
4265 0 : if ( stats_metadata%icorr_hmx_hmy_2(hm_idx_jvar,hm_idx_ivar) > 0 ) then
4266 : ! call stat_update_var( stats_metadata%icorr_hmx_hmy_2(hm_idx_jvar,hm_idx_ivar), &
4267 : ! corr_array_2(jvar,ivar,:), stats_zt )
4268 0 : do k = 1, nz, 1
4269 : call stat_update_var_pt( stats_metadata%icorr_hmx_hmy_2(hm_idx_jvar,hm_idx_ivar), k, & ! in
4270 0 : corr_array_2(k,jvar,ivar), & ! intent(in)
4271 0 : stats_zt ) ! intent(inout)
4272 : end do ! k = 1, nz, 1
4273 : end if
4274 :
4275 : end do ! jvar = ivar+1, pdf_dim, 1
4276 : end do ! ivar = iiPDF_Ncn+1, pdf_dim, 1
4277 :
4278 : end if ! stats_metadata%l_stats_samp
4279 :
4280 :
4281 0 : return
4282 :
4283 : end subroutine pdf_param_hm_stats
4284 :
4285 : !=============================================================================
4286 0 : subroutine pdf_param_ln_hm_stats( nz, pdf_dim, hm_metadata, &
4287 0 : mu_x_1_n, mu_x_2_n, &
4288 0 : sigma_x_1_n, sigma_x_2_n, &
4289 0 : corr_array_1_n, &
4290 0 : corr_array_2_n, &
4291 : stats_metadata, &
4292 : stats_zt )
4293 :
4294 : ! Description:
4295 : ! Record statistics for normal space PDF parameters involving hydrometeors.
4296 :
4297 : ! References:
4298 : !-----------------------------------------------------------------------
4299 :
4300 : use index_mapping, only: &
4301 : pdf2hydromet_idx ! Procedure(s)
4302 :
4303 : use corr_varnce_module, only: &
4304 : hm_metadata_type
4305 :
4306 : use clubb_precision, only: &
4307 : core_rknd ! Variable(s)
4308 :
4309 : use stats_type_utilities, only: &
4310 : stat_update_var_pt ! Procedure(s)
4311 :
4312 : use stats_type, only: &
4313 : stats ! Type
4314 :
4315 : use stats_variables, only: &
4316 : stats_metadata_type
4317 :
4318 : implicit none
4319 :
4320 : !------------------------ Input Variables ------------------------
4321 : integer, intent(in) :: &
4322 : nz, & ! Number of vertical levels
4323 : pdf_dim ! Number of variables in the correlation array
4324 :
4325 : type (hm_metadata_type), intent(in) :: &
4326 : hm_metadata
4327 :
4328 : real( kind = core_rknd ), dimension(nz,pdf_dim), intent(in) :: &
4329 : mu_x_1_n, & ! Mean array (normal space): PDF vars. (comp. 1) [un. vary]
4330 : mu_x_2_n, & ! Mean array (normal space): PDF vars. (comp. 2) [un. vary]
4331 : sigma_x_1_n, & ! Std. dev. array (normal space): PDF vars (comp. 1) [u.v.]
4332 : sigma_x_2_n ! Std. dev. array (normal space): PDF vars (comp. 2) [u.v.]
4333 :
4334 : real( kind = core_rknd ), dimension(nz,pdf_dim,pdf_dim), &
4335 : intent(in) :: &
4336 : corr_array_1_n, & ! Corr. array (normal space) of PDF vars. (comp. 1) [-]
4337 : corr_array_2_n ! Corr. array (normal space) of PDF vars. (comp. 2) [-]
4338 :
4339 : type (stats_metadata_type), intent(in) :: &
4340 : stats_metadata
4341 :
4342 : !------------------------ InOut Variables ------------------------
4343 : type (stats), target, intent(inout) :: &
4344 : stats_zt
4345 :
4346 : !------------------------ Local Variables ------------------------
4347 : real( kind = core_rknd ), dimension(nz) :: &
4348 0 : mu_hm_1_n, & ! Mean of ln hm (1st PDF component) [units vary]
4349 0 : mu_hm_2_n, & ! Mean of ln hm (2nd PDF component) [units vary]
4350 0 : mu_Ncn_1_n, & ! Mean of ln Ncn (1st PDF component) [ln(num/kg)]
4351 0 : mu_Ncn_2_n ! Mean of ln Ncn (2nd PDF component) [ln(num/kg)]
4352 :
4353 : integer :: ivar, jvar, hm_idx, hm_idx_ivar, hm_idx_jvar, k ! Indices
4354 :
4355 : ! Just to avoid typing hm_metadata%iiPDF_x everywhere
4356 : integer :: &
4357 : iiPDF_chi, &
4358 : iiPDF_eta, &
4359 : iiPDF_w, &
4360 : iiPDF_Ncn
4361 :
4362 : !------------------------ Begin Code ------------------------
4363 :
4364 0 : iiPDF_chi = hm_metadata%iiPDF_chi
4365 0 : iiPDF_eta = hm_metadata%iiPDF_eta
4366 0 : iiPDF_w = hm_metadata%iiPDF_w
4367 0 : iiPDF_Ncn = hm_metadata%iiPDF_Ncn
4368 :
4369 : !!! Output the statistics for normal space hydrometeor PDF parameters.
4370 :
4371 : ! Statistics
4372 0 : if ( stats_metadata%l_stats_samp ) then
4373 :
4374 : ! Switch back to using stat_update_var once the code is generalized
4375 : ! to pass in the number of vertical levels.
4376 0 : do ivar = iiPDF_Ncn+1, pdf_dim, 1
4377 :
4378 0 : hm_idx = pdf2hydromet_idx(ivar,hm_metadata)
4379 :
4380 : ! Mean (in-precip) of ln hm in PDF component 1.
4381 0 : if ( stats_metadata%imu_hm_1_n(hm_idx) > 0 ) then
4382 0 : where ( mu_x_1_n(:,ivar) > real( -huge( 0.0 ), kind = core_rknd ) )
4383 0 : mu_hm_1_n = mu_x_1_n(:,ivar)
4384 : elsewhere
4385 : ! When hm_1 is 0 (or below tolerance value), mu_hm_1_n is -inf,
4386 : ! and is set to -huge for the default CLUBB kind. Some
4387 : ! compilers have issues outputting to stats files (in single
4388 : ! precision) when the default CLUBB kind is in double precision.
4389 : ! Set to -huge for single precision.
4390 : mu_hm_1_n = real( -huge( 0.0 ), kind = core_rknd )
4391 : endwhere
4392 : ! call stat_update_var( stats_metadata%imu_hm_1_n(hm_idx), mu_hm_1_n, stats_zt )
4393 0 : do k = 1, nz, 1
4394 0 : call stat_update_var_pt( stats_metadata%imu_hm_1_n(hm_idx), k, mu_hm_1_n(k), & ! intent(in)
4395 0 : stats_zt ) ! intent(inout)
4396 : end do ! k = 1, nz, 1
4397 : end if
4398 :
4399 : ! Mean (in-precip) of ln hm in PDF component 2.
4400 0 : if ( stats_metadata%imu_hm_2_n(hm_idx) > 0 ) then
4401 0 : where ( mu_x_2_n(:,ivar) > real( -huge( 0.0 ), kind = core_rknd ) )
4402 0 : mu_hm_2_n = mu_x_2_n(:,ivar)
4403 : elsewhere
4404 : ! When hm_2 is 0 (or below tolerance value), mu_hm_2_n is -inf,
4405 : ! and is set to -huge for the default CLUBB kind. Some
4406 : ! compilers have issues outputting to stats files (in single
4407 : ! precision) when the default CLUBB kind is in double precision.
4408 : ! Set to -huge for single precision.
4409 : mu_hm_2_n = real( -huge( 0.0 ), kind = core_rknd )
4410 : endwhere
4411 : ! call stat_update_var( stats_metadata%imu_hm_2_n(hm_idx), mu_hm_2_n, stats_zt )
4412 0 : do k = 1, nz, 1
4413 0 : call stat_update_var_pt( stats_metadata%imu_hm_2_n(hm_idx), k, mu_hm_2_n(k), & ! intent(in)
4414 0 : stats_zt ) ! intent(inout)
4415 : end do ! k = 1, nz, 1
4416 : end if
4417 :
4418 : end do ! ivar = iiPDF_Ncn+1, pdf_dim, 1
4419 :
4420 : ! Mean of ln N_cn in PDF component 1.
4421 0 : if ( stats_metadata%imu_Ncn_1_n > 0 ) then
4422 : where ( mu_x_1_n(:,iiPDF_Ncn) &
4423 0 : > real( -huge( 0.0 ), kind = core_rknd ) )
4424 0 : mu_Ncn_1_n = mu_x_1_n(:,iiPDF_Ncn)
4425 : elsewhere
4426 : ! When Ncnm is 0 (or below tolerance value), mu_Ncn_1_n is -inf,
4427 : ! and is set to -huge for the default CLUBB kind. Some compilers
4428 : ! have issues outputting to stats files (in single precision) when
4429 : ! the default CLUBB kind is in double precision.
4430 : ! Set to -huge for single precision.
4431 : mu_Ncn_1_n = real( -huge( 0.0 ), kind = core_rknd )
4432 : endwhere
4433 : ! call stat_update_var( stats_metadata%imu_Ncn_1_n, mu_Ncn_1_n, stats_zt )
4434 0 : do k = 1, nz, 1
4435 0 : call stat_update_var_pt( stats_metadata%imu_Ncn_1_n, k, mu_Ncn_1_n(k), & ! intent(in)
4436 0 : stats_zt ) ! intent(inout)
4437 : end do ! k = 1, nz, 1
4438 : end if
4439 :
4440 : ! Mean of ln N_cn in PDF component 2.
4441 0 : if ( stats_metadata%imu_Ncn_2_n > 0 ) then
4442 : where ( mu_x_2_n(:,iiPDF_Ncn) &
4443 0 : > real( -huge( 0.0 ), kind = core_rknd ) )
4444 0 : mu_Ncn_2_n = mu_x_2_n(:,iiPDF_Ncn)
4445 : elsewhere
4446 : ! When Ncnm is 0 (or below tolerance value), mu_Ncn_2_n is -inf,
4447 : ! and is set to -huge for the default CLUBB kind. Some compilers
4448 : ! have issues outputting to stats files (in single precision) when
4449 : ! the default CLUBB kind is in double precision.
4450 : ! Set to -huge for single precision.
4451 : mu_Ncn_2_n = real( -huge( 0.0 ), kind = core_rknd )
4452 : endwhere
4453 : ! call stat_update_var( stats_metadata%imu_Ncn_2_n, mu_Ncn_2_n, stats_zt )
4454 0 : do k = 1, nz, 1
4455 0 : call stat_update_var_pt( stats_metadata%imu_Ncn_2_n, k, mu_Ncn_2_n(k), & ! intent(in)
4456 0 : stats_zt ) ! intent(inout)
4457 : end do ! k = 1, nz, 1
4458 : end if
4459 :
4460 0 : do ivar = iiPDF_Ncn+1, pdf_dim, 1
4461 :
4462 0 : hm_idx = pdf2hydromet_idx(ivar,hm_metadata)
4463 :
4464 : ! Standard deviation (in-precip) of ln hm in PDF component 1.
4465 0 : if ( stats_metadata%isigma_hm_1_n(hm_idx) > 0 ) then
4466 : ! call stat_update_var( stats_metadata%isigma_hm_1_n(hm_idx), sigma_x_1_n(ivar,:), &
4467 : ! stats_zt )
4468 0 : do k = 1, nz, 1
4469 : call stat_update_var_pt( stats_metadata%isigma_hm_1_n(hm_idx), k, & ! intent(in)
4470 0 : sigma_x_1_n(k,ivar), & ! intent(in)
4471 0 : stats_zt ) ! intent(inout)
4472 : end do ! k = 1, nz, 1
4473 : end if
4474 :
4475 : ! Standard deviation (in-precip) of ln hm in PDF component 2.
4476 0 : if ( stats_metadata%isigma_hm_2_n(hm_idx) > 0 ) then
4477 : ! call stat_update_var( stats_metadata%isigma_hm_2_n(hm_idx), sigma_x_2_n(ivar,:), &
4478 : ! stats_zt )
4479 0 : do k = 1, nz, 1
4480 : call stat_update_var_pt( stats_metadata%isigma_hm_2_n(hm_idx), k, & ! intent(in)
4481 0 : sigma_x_2_n(k,ivar), & ! intent(in)
4482 0 : stats_zt ) ! intent(inout)
4483 : end do ! k = 1, nz, 1
4484 : end if
4485 :
4486 : end do ! ivar = iiPDF_Ncn+1, pdf_dim, 1
4487 :
4488 : ! Standard deviation of ln N_cn in PDF component 1.
4489 0 : if ( stats_metadata%isigma_Ncn_1_n > 0 ) then
4490 : ! call stat_update_var( stats_metadata%isigma_Ncn_1_n, sigma_x_1_n(iiPDF_Ncn,:), &
4491 : ! stats_zt )
4492 0 : do k = 1, nz, 1
4493 : call stat_update_var_pt( stats_metadata%isigma_Ncn_1_n, k, & ! intent(in)
4494 0 : sigma_x_1_n(k,iiPDF_Ncn), & ! intent(in)
4495 0 : stats_zt ) ! intent(inout)
4496 : end do ! k = 1, nz, 1
4497 : end if
4498 :
4499 : ! Standard deviation of ln N_cn in PDF component 2.
4500 0 : if ( stats_metadata%isigma_Ncn_2_n > 0 ) then
4501 : ! call stat_update_var( stats_metadata%isigma_Ncn_2_n, sigma_x_2_n(iiPDF_Ncn,:), &
4502 : ! stats_zt )
4503 0 : do k = 1, nz, 1
4504 : call stat_update_var_pt( stats_metadata%isigma_Ncn_2_n, k, & ! intent(in)
4505 0 : sigma_x_2_n(k,iiPDF_Ncn), & ! intent(in)
4506 0 : stats_zt ) ! intent(inout)
4507 : end do ! k = 1, nz, 1
4508 : end if
4509 :
4510 0 : do ivar = iiPDF_Ncn+1, pdf_dim, 1
4511 :
4512 0 : hm_idx = pdf2hydromet_idx(ivar,hm_metadata)
4513 :
4514 : ! Correlation (in-precip) of w and ln hm in PDF component 1.
4515 0 : if ( stats_metadata%icorr_w_hm_1_n(hm_idx) > 0 ) then
4516 : ! call stat_update_var( stats_metadata%icorr_w_hm_1_n(hm_idx), &
4517 : ! corr_array_1_n(ivar,iiPDF_w,:), stats_zt )
4518 0 : do k = 1, nz, 1
4519 : call stat_update_var_pt( stats_metadata%icorr_w_hm_1_n(hm_idx), k, & ! intent(in)
4520 0 : corr_array_1_n(k,ivar,iiPDF_w), & ! intent(in)
4521 0 : stats_zt ) ! intent(inout)
4522 : end do ! k = 1, nz, 1
4523 : end if
4524 :
4525 : ! Correlation (in-precip) of w and ln hm in PDF component 2.
4526 0 : if ( stats_metadata%icorr_w_hm_2_n(hm_idx) > 0 ) then
4527 : ! call stat_update_var( stats_metadata%icorr_w_hm_2_n(hm_idx), &
4528 : ! corr_array_2_n(ivar,iiPDF_w,:), stats_zt )
4529 0 : do k = 1, nz, 1
4530 : call stat_update_var_pt( stats_metadata%icorr_w_hm_2_n(hm_idx), k, & ! intent(in)
4531 0 : corr_array_2_n(k,ivar,iiPDF_w), & ! intent(in)
4532 0 : stats_zt ) ! intent(inout)
4533 : end do ! k = 1, nz, 1
4534 : end if
4535 :
4536 : end do ! ivar = iiPDF_Ncn+1, pdf_dim, 1
4537 :
4538 : ! Correlation of w and ln N_cn in PDF component 1.
4539 0 : if ( stats_metadata%icorr_w_Ncn_1_n > 0 ) then
4540 : ! call stat_update_var( stats_metadata%icorr_w_Ncn_1_n, &
4541 : ! corr_array_1_n(iiPDF_Ncn,iiPDF_w,:), stats_zt )
4542 0 : do k = 1, nz, 1
4543 : call stat_update_var_pt( stats_metadata%icorr_w_Ncn_1_n, k, & ! intent(in)
4544 0 : corr_array_1_n(k,iiPDF_Ncn,iiPDF_w), & ! intent(in)
4545 0 : stats_zt ) ! intent(inout)
4546 : end do ! k = 1, nz, 1
4547 : end if
4548 :
4549 : ! Correlation of w and ln N_cn in PDF component 2.
4550 0 : if ( stats_metadata%icorr_w_Ncn_2_n > 0 ) then
4551 : ! call stat_update_var( stats_metadata%icorr_w_Ncn_2_n, &
4552 : ! corr_array_2_n(iiPDF_Ncn,iiPDF_w,:), stats_zt )
4553 0 : do k = 1, nz, 1
4554 : call stat_update_var_pt( stats_metadata%icorr_w_Ncn_2_n, k, & ! intent(in)
4555 0 : corr_array_2_n(k,iiPDF_Ncn,iiPDF_w), & ! intent(in)
4556 0 : stats_zt ) ! intent(inout)
4557 : end do ! k = 1, nz, 1
4558 : end if
4559 :
4560 0 : do ivar = iiPDF_Ncn+1, pdf_dim, 1
4561 :
4562 0 : hm_idx = pdf2hydromet_idx(ivar,hm_metadata)
4563 :
4564 : ! Correlation (in-precip) of chi (old s) and ln hm in PDF component 1.
4565 0 : if ( stats_metadata%icorr_chi_hm_1_n(hm_idx) > 0 ) then
4566 : ! call stat_update_var( stats_metadata%icorr_chi_hm_1_n(hm_idx), &
4567 : ! corr_array_1_n(ivar,iiPDF_chi,:), stats_zt )
4568 0 : do k = 1, nz, 1
4569 : call stat_update_var_pt( stats_metadata%icorr_chi_hm_1_n(hm_idx), k, & ! intent(in)
4570 0 : corr_array_1_n(k,ivar,iiPDF_chi), & ! intent(in)
4571 0 : stats_zt ) ! intent(inout)
4572 : end do ! k = 1, nz, 1
4573 : end if
4574 :
4575 : ! Correlation (in-precip) of chi( old s) and ln hm in PDF component 2.
4576 0 : if ( stats_metadata%icorr_chi_hm_2_n(hm_idx) > 0 ) then
4577 : ! call stat_update_var( stats_metadata%icorr_chi_hm_2_n(hm_idx), &
4578 : ! corr_array_2_n(ivar,iiPDF_chi,:), stats_zt )
4579 0 : do k = 1, nz, 1
4580 : call stat_update_var_pt( stats_metadata%icorr_chi_hm_2_n(hm_idx), k, & ! intent(in)
4581 0 : corr_array_2_n(k,ivar,iiPDF_chi), & ! intent(in)
4582 0 : stats_zt ) ! intent(inout)
4583 : end do ! k = 1, nz, 1
4584 : end if
4585 :
4586 : end do ! ivar = iiPDF_Ncn+1, pdf_dim, 1
4587 :
4588 : ! Correlation of chi (old s) and ln N_cn in PDF component 1.
4589 0 : if ( stats_metadata%icorr_chi_Ncn_1_n > 0 ) then
4590 : ! call stat_update_var( stats_metadata%icorr_chi_Ncn_1_n, &
4591 : ! corr_array_1_n(iiPDF_Ncn,iiPDF_chi,:), &
4592 : ! stats_zt )
4593 0 : do k = 1, nz, 1
4594 : call stat_update_var_pt( stats_metadata%icorr_chi_Ncn_1_n, k, & ! intent(in)
4595 0 : corr_array_1_n(k,iiPDF_Ncn,iiPDF_chi), & ! intent(in)
4596 0 : stats_zt ) ! intent(inout)
4597 : end do ! k = 1, nz, 1
4598 : end if
4599 :
4600 : ! Correlation of chi(old s) and ln N_cn in PDF component 2.
4601 0 : if ( stats_metadata%icorr_chi_Ncn_2_n > 0 ) then
4602 : ! call stat_update_var( stats_metadata%icorr_chi_Ncn_2_n, &
4603 : ! corr_array_2_n(iiPDF_Ncn,iiPDF_chi,:), &
4604 : ! stats_zt )
4605 0 : do k = 1, nz, 1
4606 : call stat_update_var_pt( stats_metadata%icorr_chi_Ncn_2_n, k, & ! intent(in)
4607 0 : corr_array_2_n(k,iiPDF_Ncn,iiPDF_chi), & ! intent(in)
4608 0 : stats_zt ) ! intent(inout)
4609 : end do ! k = 1, nz, 1
4610 : end if
4611 :
4612 0 : do ivar = iiPDF_Ncn+1, pdf_dim, 1
4613 :
4614 0 : hm_idx = pdf2hydromet_idx(ivar,hm_metadata)
4615 :
4616 : ! Correlation (in-precip) of eta (old t) and ln hm in PDF component 1.
4617 0 : if ( stats_metadata%icorr_eta_hm_1_n(hm_idx) > 0 ) then
4618 : ! call stat_update_var( stats_metadata%icorr_eta_hm_1_n(hm_idx), &
4619 : ! corr_array_1_n(ivar,iiPDF_eta,:), stats_zt )
4620 0 : do k = 1, nz, 1
4621 : call stat_update_var_pt( stats_metadata%icorr_eta_hm_1_n(hm_idx), k, & ! intent(in)
4622 0 : corr_array_1_n(k,ivar,iiPDF_eta), & ! intent(in)
4623 0 : stats_zt ) ! intent(inout)
4624 : end do ! k = 1, nz, 1
4625 : end if
4626 :
4627 : ! Correlation (in-precip) of eta (old t) and ln hm in PDF component 2.
4628 0 : if ( stats_metadata%icorr_eta_hm_2_n(hm_idx) > 0 ) then
4629 : ! call stat_update_var( stats_metadata%icorr_eta_hm_2_n(hm_idx), &
4630 : ! corr_array_2_n(ivar,iiPDF_eta,:), stats_zt )
4631 0 : do k = 1, nz, 1
4632 : call stat_update_var_pt( stats_metadata%icorr_eta_hm_2_n(hm_idx), k, & ! intent(in)
4633 0 : corr_array_2_n(k,ivar,iiPDF_eta), & ! intent(in)
4634 0 : stats_zt ) ! intent(inout)
4635 : end do ! k = 1, nz, 1
4636 : end if
4637 :
4638 : end do ! ivar = iiPDF_Ncn+1, pdf_dim, 1
4639 :
4640 : ! Correlation of eta (old t) and ln N_cn in PDF component 1.
4641 0 : if ( stats_metadata%icorr_eta_Ncn_1_n > 0 ) then
4642 : ! call stat_update_var( stats_metadata%icorr_eta_Ncn_1_n, &
4643 : ! corr_array_1_n(iiPDF_Ncn,iiPDF_eta,:), &
4644 : ! stats_zt )
4645 0 : do k = 1, nz, 1
4646 : call stat_update_var_pt( stats_metadata%icorr_eta_Ncn_1_n, k, & ! intent(in)
4647 0 : corr_array_1_n(k,iiPDF_Ncn,iiPDF_eta), & ! intent(in)
4648 0 : stats_zt ) ! intent(inout)
4649 : end do ! k = 1, nz, 1
4650 : end if
4651 :
4652 : ! Correlation of eta (old t) and ln N_cn in PDF component 2.
4653 0 : if ( stats_metadata%icorr_eta_Ncn_2_n > 0 ) then
4654 : ! call stat_update_var( stats_metadata%icorr_eta_Ncn_2_n, &
4655 : ! corr_array_2_n(iiPDF_Ncn,iiPDF_eta,:), &
4656 : ! stats_zt )
4657 0 : do k = 1, nz, 1
4658 : call stat_update_var_pt( stats_metadata%icorr_eta_Ncn_2_n, k, & ! intent(in)
4659 0 : corr_array_2_n(k,iiPDF_Ncn,iiPDF_eta), & ! intent(in)
4660 0 : stats_zt ) ! intent(inout)
4661 : end do ! k = 1, nz, 1
4662 : end if
4663 :
4664 0 : do ivar = iiPDF_Ncn+1, pdf_dim, 1
4665 :
4666 0 : hm_idx = pdf2hydromet_idx(ivar,hm_metadata)
4667 :
4668 : ! Correlation (in-precip) of ln N_cn and ln hm in PDF
4669 : ! component 1.
4670 0 : if ( stats_metadata%icorr_Ncn_hm_1_n(hm_idx) > 0 ) then
4671 : ! call stat_update_var( stats_metadata%icorr_Ncn_hm_1_n(hm_idx), &
4672 : ! corr_array_1_n(ivar,iiPDF_Ncn,:), stats_zt )
4673 0 : do k = 1, nz, 1
4674 : call stat_update_var_pt( stats_metadata%icorr_Ncn_hm_1_n(hm_idx), k, & ! intent(in)
4675 0 : corr_array_1_n(k,ivar,iiPDF_Ncn), & ! intent(in)
4676 0 : stats_zt ) ! intent(inout)
4677 : end do ! k = 1, nz, 1
4678 : end if
4679 :
4680 : ! Correlation (in-precip) of ln N_cn and ln hm in PDF
4681 : ! component 2.
4682 0 : if ( stats_metadata%icorr_Ncn_hm_2_n(hm_idx) > 0 ) then
4683 : ! call stat_update_var( stats_metadata%icorr_Ncn_hm_2_n(hm_idx), &
4684 : ! corr_array_2_n(ivar,iiPDF_Ncn,:), stats_zt )
4685 0 : do k = 1, nz, 1
4686 : call stat_update_var_pt( stats_metadata%icorr_Ncn_hm_2_n(hm_idx), k, & ! intent(in)
4687 0 : corr_array_2_n(k,ivar,iiPDF_Ncn), & ! intent(in)
4688 0 : stats_zt )! intent(inout)
4689 : end do ! k = 1, nz, 1
4690 : end if
4691 :
4692 : end do ! ivar = iiPDF_Ncn+1, pdf_dim, 1
4693 :
4694 0 : do ivar = iiPDF_Ncn+1, pdf_dim, 1
4695 :
4696 0 : hm_idx_ivar = pdf2hydromet_idx(ivar,hm_metadata)
4697 :
4698 0 : do jvar = ivar+1, pdf_dim, 1
4699 :
4700 0 : hm_idx_jvar= pdf2hydromet_idx(jvar,hm_metadata)
4701 :
4702 : ! Correlation (in-precip) of ln hmx and ln hmy (two different
4703 : ! hydrometeors) in PDF component 1.
4704 0 : if ( stats_metadata%icorr_hmx_hmy_1_n(hm_idx_jvar,hm_idx_ivar) > 0 ) then
4705 : ! call stat_update_var( &
4706 : ! stats_metadata%icorr_hmx_hmy_1_n(hm_idx_jvar,hm_idx_ivar), &
4707 : ! corr_array_1_n(jvar,ivar,:), stats_zt )
4708 0 : do k = 1, nz, 1
4709 : call stat_update_var_pt( &
4710 : stats_metadata%icorr_hmx_hmy_1_n(hm_idx_jvar,hm_idx_ivar), k, & ! intent(in)
4711 0 : corr_array_1_n(k,jvar,ivar), & ! intent(in)
4712 0 : stats_zt ) ! intent(inout)
4713 : end do ! k = 1, nz, 1
4714 : end if
4715 :
4716 : ! Correlation (in-precip) of ln hmx and ln hmy (two different
4717 : ! hydrometeors) in PDF component 2.
4718 0 : if ( stats_metadata%icorr_hmx_hmy_2_n(hm_idx_jvar,hm_idx_ivar) > 0 ) then
4719 : ! call stat_update_var( &
4720 : ! stats_metadata%icorr_hmx_hmy_2_n(hm_idx_jvar,hm_idx_ivar), &
4721 : ! corr_array_2_n(jvar,ivar,:), stats_zt )
4722 0 : do k = 1, nz, 1
4723 : call stat_update_var_pt( &
4724 : stats_metadata%icorr_hmx_hmy_2_n(hm_idx_jvar,hm_idx_ivar), k, & ! intent(in)
4725 0 : corr_array_2_n(k,jvar,ivar), & ! intent(in)
4726 0 : stats_zt ) ! intent(inout)
4727 : end do ! k = 1, nz, 1
4728 : end if
4729 :
4730 : end do ! jvar = ivar+1, pdf_dim, 1
4731 : end do ! ivar = iiPDF_Ncn+1, pdf_dim, 1
4732 :
4733 : end if ! stats_metadata%l_stats_samp
4734 :
4735 :
4736 0 : return
4737 :
4738 : end subroutine pdf_param_ln_hm_stats
4739 :
4740 : !=============================================================================
4741 0 : subroutine pack_hydromet_pdf_params( nz, ngrdcol, & ! In
4742 : hydromet_dim, hm_metadata, & ! In
4743 0 : hm_1, hm_2, pdf_dim, mu_x_1, & ! In
4744 0 : mu_x_2, sigma_x_1, sigma_x_2, & ! In
4745 0 : corr_array_1, corr_array_2, & ! In
4746 0 : hydromet_pdf_params ) ! Out
4747 :
4748 : ! Description:
4749 : ! Pack the standard means and variances involving hydrometeors, as well as a
4750 : ! few other variables, into the structure hydromet_pdf_params.
4751 :
4752 : ! References:
4753 : !-----------------------------------------------------------------------
4754 :
4755 : use constants_clubb, only: &
4756 : one ! Constant(s)
4757 :
4758 : use hydromet_pdf_parameter_module, only: &
4759 : hydromet_pdf_parameter ! Variable(s)
4760 :
4761 : use index_mapping, only: &
4762 : hydromet2pdf_idx ! Procedure(s)
4763 :
4764 : use corr_varnce_module, only: &
4765 : hm_metadata_type
4766 :
4767 : use clubb_precision, only: &
4768 : core_rknd ! Variable(s)
4769 :
4770 : implicit none
4771 :
4772 : !-------------------------- Input Variables --------------------------
4773 : integer, intent(in) :: &
4774 : nz, & ! Number of vertical grid levels
4775 : ngrdcol, & ! Number of grid columns
4776 : hydromet_dim ! Number of hydrometeor species
4777 :
4778 : real( kind = core_rknd ), dimension(ngrdcol,nz,hydromet_dim), intent(in) :: &
4779 : hm_1, & ! Mean of a precip. hydrometeor (1st PDF component) [units vary]
4780 : hm_2 ! Mean of a precip. hydrometeor (2nd PDF component) [units vary]
4781 :
4782 : integer, intent(in) :: &
4783 : pdf_dim ! Number of variables in the mean/stdev arrays
4784 :
4785 : type (hm_metadata_type), intent(in) :: &
4786 : hm_metadata
4787 :
4788 : real( kind = core_rknd ), dimension(ngrdcol,nz,pdf_dim), intent(in) :: &
4789 : mu_x_1, & ! Mean array of PDF vars. (1st PDF component) [units vary]
4790 : mu_x_2, & ! Mean array of PDF vars. (2nd PDF component) [units vary]
4791 : sigma_x_1, & ! Standard deviation array of PDF vars (comp. 1) [units vary]
4792 : sigma_x_2 ! Standard deviation array of PDF vars (comp. 2) [units vary]
4793 :
4794 : real( kind = core_rknd ), dimension(ngrdcol,nz,pdf_dim,pdf_dim), intent(in) :: &
4795 : corr_array_1, & ! Correlation array of PDF vars. (comp. 1) [-]
4796 : corr_array_2 ! Correlation array of PDF vars. (comp. 2) [-]
4797 :
4798 : !-------------------------- Output Variable --------------------------
4799 : type(hydromet_pdf_parameter), dimension(ngrdcol,nz), intent(out) :: &
4800 : hydromet_pdf_params ! Hydrometeor PDF parameters [units vary]
4801 :
4802 : !-------------------------- Local Variables --------------------------
4803 : integer :: ivar, jvar, pdf_idx ! Indices
4804 :
4805 : ! Just to avoid typing hm_metadata%iiPDF_x everywhere
4806 : integer :: &
4807 : iiPDF_chi, &
4808 : iiPDF_eta, &
4809 : iiPDF_w, &
4810 : iiPDF_Ncn
4811 :
4812 : !------------------------ Begin Code ------------------------
4813 :
4814 0 : iiPDF_chi = hm_metadata%iiPDF_chi
4815 0 : iiPDF_eta = hm_metadata%iiPDF_eta
4816 0 : iiPDF_w = hm_metadata%iiPDF_w
4817 0 : iiPDF_Ncn = hm_metadata%iiPDF_Ncn
4818 :
4819 :
4820 : ! Pack remaining means and standard deviations into hydromet_pdf_params.
4821 0 : do ivar = 1, hydromet_dim, 1
4822 :
4823 0 : pdf_idx = hydromet2pdf_idx(ivar,hm_metadata)
4824 :
4825 : ! Mean of a hydrometeor (overall) in the 1st PDF component.
4826 0 : hydromet_pdf_params(:,:)%hm_1(ivar) = hm_1(:,:,ivar)
4827 : ! Mean of a hydrometeor (overall) in the 2nd PDF component.
4828 0 : hydromet_pdf_params(:,:)%hm_2(ivar) = hm_2(:,:,ivar)
4829 :
4830 : ! Mean of a hydrometeor (in-precip) in the 1st PDF component.
4831 0 : hydromet_pdf_params(:,:)%mu_hm_1(ivar) = mu_x_1(:,:,pdf_idx)
4832 : ! Mean of a hydrometeor (in-precip) in the 2nd PDF component.
4833 0 : hydromet_pdf_params(:,:)%mu_hm_2(ivar) = mu_x_2(:,:,pdf_idx)
4834 :
4835 : ! Standard deviation of a hydrometeor (in-precip) in the
4836 : ! 1st PDF component.
4837 0 : hydromet_pdf_params(:,:)%sigma_hm_1(ivar) = sigma_x_1(:,:,pdf_idx)
4838 : ! Standard deviation of a hydrometeor (in-precip) in the
4839 : ! 2nd PDF component.
4840 0 : hydromet_pdf_params(:,:)%sigma_hm_2(ivar) = sigma_x_2(:,:,pdf_idx)
4841 :
4842 : ! Correlation (in-precip) of w and a hydrometeor in the 1st PDF
4843 : ! component.
4844 0 : hydromet_pdf_params(:,:)%corr_w_hm_1(ivar) = corr_array_1(:,:,pdf_idx,iiPDF_w)
4845 :
4846 : ! Correlation (in-precip) of w and a hydrometeor in the 2nd PDF
4847 : ! component.
4848 0 : hydromet_pdf_params(:,:)%corr_w_hm_2(ivar) = corr_array_2(:,:,pdf_idx,iiPDF_w)
4849 :
4850 : ! Correlation (in-precip) of chi and a hydrometeor in the 1st PDF
4851 : ! component.
4852 0 : hydromet_pdf_params(:,:)%corr_chi_hm_1(ivar) &
4853 0 : = corr_array_1(:,:,pdf_idx,iiPDF_chi)
4854 :
4855 : ! Correlation (in-precip) of chi and a hydrometeor in the 2nd PDF
4856 : ! component.
4857 0 : hydromet_pdf_params(:,:)%corr_chi_hm_2(ivar) &
4858 0 : = corr_array_2(:,:,pdf_idx,iiPDF_chi)
4859 :
4860 : ! Correlation (in-precip) of eta and a hydrometeor in the 1st PDF
4861 : ! component.
4862 0 : hydromet_pdf_params(:,:)%corr_eta_hm_1(ivar) &
4863 0 : = corr_array_1(:,:,pdf_idx,iiPDF_eta)
4864 :
4865 : ! Correlation (in-precip) of eta and a hydrometeor in the 2nd PDF
4866 : ! component.
4867 0 : hydromet_pdf_params(:,:)%corr_eta_hm_2(ivar) &
4868 0 : = corr_array_2(:,:,pdf_idx,iiPDF_eta)
4869 :
4870 : ! Correlation (in-precip) of two hydrometeors, hmx and hmy, in the 1st
4871 : ! PDF component.
4872 0 : hydromet_pdf_params(:,:)%corr_hmx_hmy_1(ivar,ivar) = one
4873 :
4874 0 : do jvar = ivar+1, hydromet_dim, 1
4875 :
4876 0 : hydromet_pdf_params(:,:)%corr_hmx_hmy_1(jvar,ivar) &
4877 0 : = corr_array_1(:,:,hydromet2pdf_idx(jvar,hm_metadata),pdf_idx)
4878 :
4879 0 : hydromet_pdf_params(:,:)%corr_hmx_hmy_1(ivar,jvar) &
4880 0 : = hydromet_pdf_params(:,:)%corr_hmx_hmy_1(jvar,ivar)
4881 :
4882 : end do ! jvar = ivar+1, hydromet_dim, 1
4883 :
4884 : ! Correlation (in-precip) of two hydrometeors, hmx and hmy, in the 2nd
4885 : ! PDF component.
4886 0 : hydromet_pdf_params(:,:)%corr_hmx_hmy_2(ivar,ivar) = one
4887 :
4888 0 : do jvar = ivar+1, hydromet_dim, 1
4889 :
4890 0 : hydromet_pdf_params(:,:)%corr_hmx_hmy_2(jvar,ivar) &
4891 0 : = corr_array_2(:,:,hydromet2pdf_idx(jvar,hm_metadata),pdf_idx)
4892 :
4893 0 : hydromet_pdf_params(:,:)%corr_hmx_hmy_2(ivar,jvar) &
4894 0 : = hydromet_pdf_params(:,:)%corr_hmx_hmy_2(jvar,ivar)
4895 :
4896 : end do ! jvar = ivar+1, hydromet_dim, 1
4897 :
4898 : end do ! ivar = 1, hydromet_dim, 1
4899 :
4900 : ! Mean of Ncn (overall) in the 1st PDF component.
4901 0 : hydromet_pdf_params(:,:)%mu_Ncn_1 = mu_x_1(:,:,iiPDF_Ncn)
4902 : ! Mean of Ncn (overall) in the 2nd PDF component.
4903 0 : hydromet_pdf_params(:,:)%mu_Ncn_2 = mu_x_2(:,:,iiPDF_Ncn)
4904 :
4905 : ! Standard deviation of Ncn (overall) in the 1st PDF component.
4906 0 : hydromet_pdf_params(:,:)%sigma_Ncn_1 = sigma_x_1(:,:,iiPDF_Ncn)
4907 : ! Standard deviation of Ncn (overall) in the 2nd PDF component.
4908 0 : hydromet_pdf_params(:,:)%sigma_Ncn_2 = sigma_x_2(:,:,iiPDF_Ncn)
4909 :
4910 :
4911 0 : return
4912 :
4913 : end subroutine pack_hydromet_pdf_params
4914 :
4915 : !=============================================================================
4916 0 : elemental function compute_rtp2_from_chi( sigma_chi_1, sigma_chi_2, &
4917 : sigma_eta_1, sigma_eta_2, &
4918 : rt_1, rt_2, &
4919 : crt_1, crt_2, &
4920 : mixt_frac, &
4921 : corr_chi_eta_1, corr_chi_eta_2 ) &
4922 : result( rtp2_zt_from_chi )
4923 :
4924 : ! Description:
4925 : ! Compute the variance of rt from the distribution of chi and eta. The
4926 : ! resulting variance will be consistent with CLUBB's extended PDF
4927 : ! involving chi and eta, including if l_fix_w_chi_eta_correlations = .true..
4928 :
4929 : ! References:
4930 : ! None
4931 : !-----------------------------------------------------------------------
4932 :
4933 : use clubb_precision, only: &
4934 : core_rknd ! Constant
4935 :
4936 : use pdf_utilities, only: &
4937 : compute_variance_binormal ! Procedure
4938 :
4939 : use constants_clubb, only: &
4940 : one_half, & ! Constant(s)
4941 : one, &
4942 : two
4943 :
4944 : implicit none
4945 :
4946 : ! Input Variables
4947 : real( kind = core_rknd ), intent(in) :: &
4948 : sigma_chi_1, & ! Standard deviation of chi (1st PDF comp.) [kg/kg]
4949 : sigma_chi_2, & ! Standard deviation of chi (2nd PDF comp.) [kg/kg]
4950 : sigma_eta_1, & ! Standard deviation of eta (1st PDF comp.) [kg/kg]
4951 : sigma_eta_2, & ! Standard deviation of eta (2nd PDF comp.) [kg/kg]
4952 : crt_1, & ! Coef. of r_t in chi/eta eqns. (1st comp.) [-]
4953 : crt_2, & ! Coef. of r_t in chi/eta eqns. (2nd comp.) [-]
4954 : rt_1, & ! Mean of rt (1st PDF component) [kg/kg]
4955 : rt_2, & ! Mean of rt (2nd PDF component) [kg/kg]
4956 : mixt_frac ! Weight of 1st gaussian PDF component [-]
4957 :
4958 : real( kind = core_rknd ), intent(in) :: &
4959 : corr_chi_eta_1, & ! Correlation of chi and eta in 1st PDF component [-]
4960 : corr_chi_eta_2 ! Correlation of chi and eta in 2nd PDF component [-]
4961 :
4962 : ! Output Variable
4963 : real( kind = core_rknd ) :: &
4964 : rtp2_zt_from_chi ! Grid-box variance of rtp2 on thermo. levels [kg/kg]
4965 :
4966 : ! Local Variables
4967 : real( kind = core_rknd ) :: &
4968 : varnce_rt_1_zt_from_chi, varnce_rt_2_zt_from_chi
4969 :
4970 : real( kind = core_rknd ) :: &
4971 : rtm, & ! Mean of rt (overall) [kg/kg]
4972 : sigma_rt_1_from_chi, & ! Standard deviation of rt (1st PDF comp.) [kg/kg]
4973 : sigma_rt_2_from_chi ! Standard deviation of rt (2nd PDF comp.) [kg/kg]
4974 :
4975 :
4976 : !-----------------------------------------------------------------------
4977 :
4978 : !----- Begin Code -----
4979 :
4980 : varnce_rt_1_zt_from_chi &
4981 : = ( corr_chi_eta_1 * sigma_chi_1 * sigma_eta_1 &
4982 : + one_half * sigma_chi_1**2 + one_half * sigma_eta_1**2 ) &
4983 0 : / ( two * crt_1**2 )
4984 :
4985 : varnce_rt_2_zt_from_chi &
4986 : = ( corr_chi_eta_2 * sigma_chi_2 * sigma_eta_2 &
4987 : + one_half * sigma_chi_2**2 + one_half * sigma_eta_2**2 ) &
4988 0 : / ( two * crt_2**2 )
4989 :
4990 0 : rtm = mixt_frac*rt_1 + (one-mixt_frac)*rt_2
4991 :
4992 0 : sigma_rt_1_from_chi = sqrt( varnce_rt_1_zt_from_chi )
4993 0 : sigma_rt_2_from_chi = sqrt( varnce_rt_2_zt_from_chi )
4994 :
4995 : rtp2_zt_from_chi &
4996 : = compute_variance_binormal( rtm, rt_1, rt_2, sigma_rt_1_from_chi, &
4997 0 : sigma_rt_2_from_chi, mixt_frac )
4998 :
4999 :
5000 : return
5001 :
5002 : end function compute_rtp2_from_chi
5003 :
5004 : !===============================================================================
5005 :
5006 : end module setup_clubb_pdf_params
|