Line data Source code
1 : ! $Id$
2 : !===============================================================================
3 : module new_pdf_main
4 :
5 : ! Description:
6 : ! The portion of CLUBB's multivariate, two-component PDF that is the
7 : ! trivariate, two-component normal PDF of vertical velocity (w), total water
8 : ! mixing ratio (rt), and liquid water potential temperature (thl).
9 :
10 : ! References:
11 : !-------------------------------------------------------------------------
12 :
13 : implicit none
14 :
15 : public :: new_pdf_driver ! Procedure(s)
16 :
17 : private :: calc_responder_var, & ! Procedure(s)
18 : calc_F_x_zeta_x_setter, &
19 : calc_F_x_responder
20 :
21 : private
22 :
23 : contains
24 :
25 : !=============================================================================
26 0 : subroutine new_pdf_driver( nz, ngrdcol, wm, rtm, thlm, wp2, rtp2, thlp2, Skw, & ! In
27 0 : wprtp, wpthlp, rtpthlp, & ! In
28 0 : clubb_params, & ! In
29 0 : Skrt, Skthl, & ! I/O
30 0 : mu_w_1, mu_w_2, & ! Out
31 0 : mu_rt_1, mu_rt_2, & ! Out
32 0 : mu_thl_1, mu_thl_2, & ! Out
33 0 : sigma_w_1_sqd, sigma_w_2_sqd, & ! Out
34 0 : sigma_rt_1_sqd, sigma_rt_2_sqd, & ! Out
35 0 : sigma_thl_1_sqd, sigma_thl_2_sqd, & ! Out
36 0 : mixt_frac, & ! Out
37 : pdf_implicit_coefs_terms, & ! Out
38 0 : F_w, F_rt, F_thl, min_F_w, max_F_w, & ! Out
39 0 : min_F_rt, max_F_rt, min_F_thl, max_F_thl ) ! Out
40 :
41 :
42 : ! Description:
43 : ! Selects which variable is used to set the mixture fraction for the PDF
44 : ! ("the setter") and which variables are handled after that mixture fraction
45 : ! has been set ("the responders"). Traditionally, w has been used to set
46 : ! the PDF.
47 :
48 : ! References:
49 : !-----------------------------------------------------------------------
50 :
51 : use constants_clubb, only: &
52 : four, & ! Variable(s)
53 : two, &
54 : one, &
55 : zero, &
56 : rt_tol, &
57 : thl_tol, &
58 : max_mag_correlation
59 :
60 : use new_pdf, only: &
61 : calc_setter_var_params, & ! Procedure(s)
62 : calc_coef_wp4_implicit, &
63 : calc_coef_wpxp2_implicit, &
64 : calc_coefs_wp2xp_semiimpl, &
65 : calc_coefs_wpxpyp_semiimpl
66 :
67 : use pdf_parameter_module, only: &
68 : implicit_coefs_terms ! Variable Type
69 :
70 : use model_flags, only: &
71 : l_explicit_turbulent_adv_wp3, & ! Variable(s)
72 : l_explicit_turbulent_adv_wpxp, &
73 : l_explicit_turbulent_adv_xpyp
74 :
75 : use clubb_precision, only: &
76 : core_rknd ! Variable(s)
77 :
78 : use parameter_indices, only: &
79 : nparams, & ! Variable(s)
80 : islope_coef_spread_DG_means_w, &
81 : ipdf_component_stdev_factor_w, &
82 : icoef_spread_DG_means_rt, &
83 : icoef_spread_DG_means_thl
84 :
85 : implicit none
86 :
87 : integer, intent(in) :: &
88 : nz, &
89 : ngrdcol
90 :
91 : ! Input Variables
92 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
93 : wm, & ! Mean of w (overall) [m/s]
94 : rtm, & ! Mean of rt (overall) [kg/kg]
95 : thlm, & ! Mean of thl (overall) [K]
96 : wp2, & ! Variance of w (overall) [m^2/s^2]
97 : rtp2, & ! Variance of rt (overall) [kg^2/kg^2]
98 : thlp2, & ! Variance of thl (overall) [K^2]
99 : Skw, & ! Skewness of w (overall) [-]
100 : wprtp, & ! Covariance of w and rt (overall) [(m/s)kg/kg]
101 : wpthlp, & ! Covariance of w and thl (overall) [(m/s)K]
102 : rtpthlp ! Covariance of rt and thl (overall) [(kg/kg)K]
103 :
104 : real( kind = core_rknd ), dimension(ngrdcol,nparams), intent(in) :: &
105 : clubb_params ! Array of CLUBB's tunable parameters [units vary]
106 :
107 : ! Input/Output Variables
108 : ! These variables are input/output because their values may be clipped.
109 : ! Otherwise, as long as it is not necessary to clip them, their values
110 : ! will stay the same.
111 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(inout) :: &
112 : Skrt, & ! Skewness of rt (overall) [-]
113 : Skthl ! Skewness of thl (overall) [-]
114 :
115 : ! Output Variables
116 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
117 : mu_w_1, & ! Mean of w (1st PDF component) [m/s]
118 : mu_w_2, & ! Mean of w (2nd PDF component) [m/s]
119 : mu_rt_1, & ! Mean of rt (1st PDF component) [kg/kg]
120 : mu_rt_2, & ! Mean of rt (2nd PDF component) [kg/kg]
121 : mu_thl_1, & ! Mean of thl (1st PDF component) [K]
122 : mu_thl_2, & ! Mean of thl (2nd PDF component) [K]
123 : sigma_w_1_sqd, & ! Variance of w (1st PDF component) [m^2/s^2]
124 : sigma_w_2_sqd, & ! Variance of w (2nd PDF component) [m^2/s^2]
125 : sigma_rt_1_sqd, & ! Variance of rt (1st PDF component) [kg^2/kg^2]
126 : sigma_rt_2_sqd, & ! Variance of rt (2nd PDF component) [kg^2/kg^2]
127 : sigma_thl_1_sqd, & ! Variance of thl (1st PDF component) [K^2]
128 : sigma_thl_2_sqd, & ! Variance of thl (2nd PDF component) [K^2]
129 : mixt_frac ! Mixture fraction [-]
130 :
131 : type(implicit_coefs_terms), intent(inout) :: &
132 : pdf_implicit_coefs_terms ! Implicit coefs / explicit terms [units vary]
133 :
134 : ! Output only for recording statistics.
135 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
136 : F_w, & ! Parameter for the spread of the PDF component means of w [-]
137 : F_rt, & ! Parameter for the spread of the PDF component means of rt [-]
138 : F_thl ! Parameter for the spread of the PDF component means of thl [-]
139 :
140 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
141 : min_F_w, & ! Minimum allowable value of parameter F_w [-]
142 : max_F_w, & ! Maximum allowable value of parameter F_w [-]
143 : min_F_rt, & ! Minimum allowable value of parameter F_rt [-]
144 : max_F_rt, & ! Maximum allowable value of parameter F_rt [-]
145 : min_F_thl, & ! Minimum allowable value of parameter F_thl [-]
146 : max_F_thl ! Maximum allowable value of parameter F_thl [-]
147 :
148 : ! Local Variables
149 : real( kind = core_rknd ), dimension(nz) :: &
150 0 : sigma_w_1, & ! Standard deviation of w (1st PDF component) [m/s]
151 0 : sigma_w_2, & ! Standard deviation of w (2nd PDF component) [m/s]
152 0 : sgn_wprtp, & ! Sign of the covariance of w and rt (overall) [-]
153 0 : sgn_wpthlp, & ! Sign of the covariance of w and thl (overall) [-]
154 0 : sgn_wp2 ! Sign of the variance of w (overall); always pos. [-]
155 :
156 : real( kind = core_rknd ), dimension(nz) :: &
157 0 : coef_sigma_w_1_sqd, & ! sigma_w_1^2 = coef_sigma_w_1_sqd * <w'^2> [-]
158 0 : coef_sigma_w_2_sqd, & ! sigma_w_2^2 = coef_sigma_w_2_sqd * <w'^2> [-]
159 0 : coef_sigma_rt_1_sqd, & ! sigma_rt_1^2 = coef_sigma_rt_1_sqd * <rt'^2> [-]
160 0 : coef_sigma_rt_2_sqd, & ! sigma_rt_2^2 = coef_sigma_rt_2_sqd * <rt'^2> [-]
161 0 : coef_sigma_thl_1_sqd, & ! sigma_thl_1^2=coef_sigma_thl_1_sqd*<thl'^2> [-]
162 0 : coef_sigma_thl_2_sqd ! sigma_thl_2^2=coef_sigma_thl_2_sqd*<thl'^2> [-]
163 :
164 : real( kind = core_rknd ), dimension(nz) :: &
165 0 : max_Skx2_pos_Skx_sgn_wpxp, & ! Maximum Skx^2 when Skx*sgn(<w'x'>) >= 0 [-]
166 0 : max_Skx2_neg_Skx_sgn_wpxp ! Maximum Skx^2 when Skx*sgn(<w'x'>) < 0 [-]
167 :
168 : real( kind = core_rknd ), dimension(nz) :: &
169 0 : zeta_w ! Parameter for the PDF component variances of w [-]
170 :
171 : real ( kind = core_rknd ) :: &
172 : lambda_w ! Param. that increases or decreases Skw dependence [-]
173 :
174 : real ( kind = core_rknd ), dimension(nz) :: &
175 0 : exp_factor_rt, & ! Factor of the form 1 - exp{} that reduces F_rt [-]
176 0 : exp_factor_thl, & ! Don't reduce F_thl by exp_factor_thl [-]
177 0 : adj_corr_rt_thl ! Adjusted (overall) correlation of rt and theta-l [-]
178 :
179 : real ( kind = core_rknd ), dimension(nz) :: &
180 0 : coef_wp4_implicit, & ! <w'^4> = coef_wp4_implicit * <w'^2>^2 [-]
181 0 : coef_wprtp2_implicit, & ! <w'rt'^2> = coef_wprtp2_implicit*<rt'^2> [m/s]
182 0 : coef_wpthlp2_implicit ! <w'thl'^2>=coef_wpthlp2_implicit*<thl'^2> [m/s]
183 :
184 : ! <w'^2 rt'> = coef_wp2rtp_implicit * <w'rt'> + term_wp2rtp_explicit
185 : real ( kind = core_rknd ), dimension(nz) :: &
186 0 : coef_wp2rtp_implicit, & ! Coefficient that is multiplied by <w'rt'> [m/s]
187 0 : term_wp2rtp_explicit ! Term that is on the RHS [m^2/s^2 kg/kg]
188 :
189 : ! <w'^2 thl'> = coef_wp2thlp_implicit * <w'thl'> + term_wp2thlp_explicit
190 : real ( kind = core_rknd ), dimension(nz) :: &
191 0 : coef_wp2thlp_implicit, & ! Coef. that is multiplied by <w'thl'> [m/s]
192 0 : term_wp2thlp_explicit ! Term that is on the RHS [m^2/s^2 K]
193 :
194 : ! <w'rt'thl'> = coef_wprtpthlp_implicit*<rt'thl'> + term_wprtpthlp_explicit
195 : real ( kind = core_rknd ), dimension(nz) :: &
196 0 : coef_wprtpthlp_implicit, & ! Coef. that is multiplied by <rt'thl'> [m/s]
197 0 : term_wprtpthlp_explicit ! Term that is on the RHS [m/s(kg/kg)K]
198 :
199 : integer :: &
200 : i
201 :
202 : ! ------------------------- Begin Code -------------------------
203 :
204 0 : do i = 1, ngrdcol
205 :
206 : ! Calculate sgn( <w'rt'> ).
207 0 : where ( wprtp(i,:) >= zero )
208 : sgn_wprtp = one
209 : elsewhere ! wprtp < 0
210 : sgn_wprtp = -one
211 : endwhere ! wprtp >= 0
212 :
213 : ! Calculate sgn( <w'thl'> ).
214 0 : where ( wpthlp(i,:) >= zero )
215 : sgn_wpthlp = one
216 : elsewhere ! wpthlp < 0
217 : sgn_wpthlp = -one
218 : endwhere ! wpthlp >= 0
219 :
220 : ! Sign of the variance of w (overall), which is always positive.
221 0 : sgn_wp2 = one
222 :
223 0 : lambda_w = 0.5_core_rknd
224 :
225 : ! Calculate the adjusted (overall) correlation of rt and theta-l, and the
226 : ! value of exp_factor_rt.
227 0 : where ( rtp2(i,:) >= rt_tol**2 .and. thlp2(i,:) >= thl_tol**2 )
228 : adj_corr_rt_thl = rtpthlp(i,:) / sqrt( rtp2(i,:) * thlp2(i,:) ) * sgn_wprtp * sgn_wpthlp
229 : adj_corr_rt_thl = min( max( adj_corr_rt_thl, -max_mag_correlation ), &
230 : max_mag_correlation )
231 : exp_factor_rt = one &
232 : - exp( -0.2_core_rknd * ( adj_corr_rt_thl + one )**5 )
233 : elsewhere ! <rt'^2> < rt_tol^2 or <thl'^2> < thl_tol^2
234 : adj_corr_rt_thl = zero ! adj_corr_rt_thl is undefined in this scenario.
235 : exp_factor_rt = one ! Set exp_factor_rt to 1.
236 : endwhere ! <rt'^2> >= rt_tol^2 and <thl'^2> >= thl_tol^2
237 :
238 : ! The value of F_thl is not reduced by exp_factor_thl.
239 0 : exp_factor_thl = one
240 :
241 :
242 : ! Vertical velocity, w, will always be the setter variable.
243 : call calc_F_x_zeta_x_setter( nz, Skw(i,:), & ! In
244 : clubb_params(i,islope_coef_spread_DG_means_w), & ! In
245 : clubb_params(i,ipdf_component_stdev_factor_w), & ! In
246 : lambda_w, & ! In
247 : F_w(i,:), zeta_w, & ! Out
248 0 : min_F_w(i,:), max_F_w(i,:) ) ! Out
249 :
250 : ! Calculate the PDF parameters, including mixture fraction, for the
251 : ! setter variable, w.
252 : call calc_setter_var_params( nz, wm(i,:), wp2(i,:), Skw(i,:), sgn_wp2, & ! In
253 : F_w(i,:), zeta_w, & ! In
254 : mu_w_1(i,:), mu_w_2(i,:), sigma_w_1, & ! Out
255 : sigma_w_2, mixt_frac(i,:), & ! Out
256 : coef_sigma_w_1_sqd, & ! Out
257 0 : coef_sigma_w_2_sqd ) ! Out
258 :
259 0 : sigma_w_1_sqd(i,:) = sigma_w_1**2
260 0 : sigma_w_2_sqd(i,:) = sigma_w_2**2
261 :
262 : ! Calculate the upper limit on the magnitude of skewness for responding
263 : ! variables.
264 : max_Skx2_pos_Skx_sgn_wpxp = four * ( one - mixt_frac(i,:) )**2 &
265 0 : / ( mixt_frac(i,:) * ( two - mixt_frac(i,:) ) )
266 :
267 0 : max_Skx2_neg_Skx_sgn_wpxp = four * mixt_frac(i,:)**2 / ( one - mixt_frac(i,:)**2 )
268 :
269 : ! Calculate the PDF parameters for responder variable rt.
270 : call calc_responder_var( nz, rtm(i,:), rtp2(i,:), sgn_wprtp, mixt_frac(i,:), & ! In
271 : clubb_params(i,icoef_spread_DG_means_rt), & ! In
272 : exp_factor_rt, & ! In
273 : max_Skx2_pos_Skx_sgn_wpxp, & ! In
274 : max_Skx2_neg_Skx_sgn_wpxp, & ! In
275 : Skrt(i,:), & ! In/Out
276 : mu_rt_1(i,:), mu_rt_2(i,:), & ! Out
277 : sigma_rt_1_sqd(i,:), sigma_rt_2_sqd(i,:), & ! Out
278 : coef_sigma_rt_1_sqd, & ! Out
279 : coef_sigma_rt_2_sqd, & ! Out
280 0 : F_rt(i,:), min_F_rt(i,:), max_F_rt(i,:) ) ! Out
281 :
282 : ! Calculate the PDF parameters for responder variable thl.
283 : call calc_responder_var( nz, thlm(i,:), thlp2(i,:), sgn_wpthlp, mixt_frac(i,:), & ! In
284 : clubb_params(i,icoef_spread_DG_means_thl), & ! In
285 : exp_factor_thl, & ! In
286 : max_Skx2_pos_Skx_sgn_wpxp, & ! In
287 : max_Skx2_neg_Skx_sgn_wpxp, & ! In
288 : Skthl(i,:), & ! In/Out
289 : mu_thl_1(i,:), mu_thl_2(i,:), & ! Out
290 : sigma_thl_1_sqd(i,:), sigma_thl_2_sqd(i,:), & ! Out
291 : coef_sigma_thl_1_sqd, & ! Out
292 : coef_sigma_thl_2_sqd, & ! Out
293 0 : F_thl(i,:), min_F_thl(i,:), max_F_thl(i,:) ) ! Out
294 :
295 :
296 : if ( .not. l_explicit_turbulent_adv_wp3 ) then
297 :
298 : ! Turbulent advection of <w'^3> is being handled implicitly.
299 :
300 : ! <w'^4> = coef_wp4_implicit * <w'^2>^2.
301 : coef_wp4_implicit &
302 : = calc_coef_wp4_implicit( nz, mixt_frac(i,:), F_w(i,:), &
303 : coef_sigma_w_1_sqd, &
304 0 : coef_sigma_w_2_sqd )
305 :
306 : else ! l_explicit_turbulent_adv_wp3
307 :
308 : ! Turbulent advection of <w'^3> is being handled explicitly.
309 : coef_wp4_implicit = zero
310 :
311 : endif ! .not. l_explicit_turbulent_adv_wp3
312 :
313 : if ( .not. l_explicit_turbulent_adv_xpyp ) then
314 :
315 : ! Turbulent advection of <rt'^2> and <thl'^2> is being handled
316 : ! implicitly. Turbulent advection of <rt'thl'> is being handled
317 : ! semi-implicitly.
318 :
319 : ! <w'rt'^2> = coef_wprtp2_implicit * <rt'^2>
320 : coef_wprtp2_implicit &
321 : = calc_coef_wpxp2_implicit( nz, wp2(i,:), rtp2(i,:), wprtp(i,:), sgn_wprtp, &
322 : mixt_frac(i,:), F_w(i,:), F_rt(i,:), &
323 : coef_sigma_w_1_sqd, &
324 : coef_sigma_w_2_sqd, &
325 : coef_sigma_rt_1_sqd, &
326 0 : coef_sigma_rt_2_sqd )
327 :
328 : ! <w'thl'^2> = coef_wpthlp2_implicit * <thl'^2>
329 : coef_wpthlp2_implicit &
330 : = calc_coef_wpxp2_implicit( nz, wp2(i,:), thlp2(i,:), wpthlp(i,:), sgn_wpthlp, &
331 : mixt_frac(i,:), F_w(i,:), F_thl(i,:), &
332 : coef_sigma_w_1_sqd, &
333 : coef_sigma_w_2_sqd, &
334 : coef_sigma_thl_1_sqd, &
335 0 : coef_sigma_thl_2_sqd )
336 :
337 : ! <w'rt'thl'> = coef_wprtpthlp_implicit * <rt'thl'>
338 : ! + term_wprtpthlp_explicit
339 : call calc_coefs_wpxpyp_semiimpl( nz, wp2(i,:), rtp2(i,:), thlp2(i,:), wprtp(i,:), & ! In
340 : wpthlp(i,:), sgn_wprtp, sgn_wpthlp, & ! In
341 : mixt_frac(i,:), F_w(i,:), F_rt(i,:), F_thl(i,:), & ! In
342 : coef_sigma_w_1_sqd , & ! In
343 : coef_sigma_w_2_sqd, & ! In
344 : coef_sigma_rt_1_sqd, & ! In
345 : coef_sigma_rt_2_sqd, & ! In
346 : coef_sigma_thl_1_sqd, & ! In
347 : coef_sigma_thl_2_sqd, & ! In
348 : coef_wprtpthlp_implicit, & ! Out
349 0 : term_wprtpthlp_explicit ) ! Out
350 :
351 : else ! l_explicit_turbulent_adv_xpyp
352 :
353 : ! Turbulent advection of <rt'^2>, <thl'^2>, and <rt'thl'> is being
354 : ! handled explicitly.
355 : coef_wprtp2_implicit = zero
356 : coef_wpthlp2_implicit = zero
357 : coef_wprtpthlp_implicit = zero
358 : term_wprtpthlp_explicit = zero
359 :
360 : endif ! .not. l_explicit_turbulent_adv_xpyp
361 :
362 : if ( .not. l_explicit_turbulent_adv_wpxp ) then
363 :
364 : ! Turbulent advection of <w'rt'> and <w'thl'> is being handled
365 : ! semi-implicitly.
366 :
367 : ! <w'^2 rt'> = coef_wp2rtp_implicit * <w'rt'> + term_wp2rtp_explicit
368 : call calc_coefs_wp2xp_semiimpl( nz, wp2(i,:), rtp2(i,:), sgn_wprtp, & ! In
369 : mixt_frac(i,:), F_w(i,:), F_rt(i,:), & ! In
370 : coef_sigma_w_1_sqd, & ! In
371 : coef_sigma_w_2_sqd, & ! In
372 : coef_sigma_rt_1_sqd, & ! In
373 : coef_sigma_rt_2_sqd, & ! In
374 : coef_wp2rtp_implicit, & ! Out
375 0 : term_wp2rtp_explicit ) ! Out
376 :
377 : ! <w'^2 thl'> = coef_wp2thlp_implicit * <w'thl'> + term_wp2thlp_explicit
378 : call calc_coefs_wp2xp_semiimpl( nz, wp2(i,:), thlp2(i,:), sgn_wpthlp, & ! In
379 : mixt_frac(i,:), F_w(i,:), F_thl(i,:), & ! In
380 : coef_sigma_w_1_sqd, & ! In
381 : coef_sigma_w_2_sqd, & ! In
382 : coef_sigma_thl_1_sqd, & ! In
383 : coef_sigma_thl_2_sqd, & ! In
384 : coef_wp2thlp_implicit, & ! Out
385 0 : term_wp2thlp_explicit ) ! Out
386 :
387 : else ! l_explicit_turbulent_adv_wpxp
388 :
389 : ! Turbulent advection of <w'rt'> and <w'thl'> is being handled
390 : ! explicitly.
391 : coef_wp2rtp_implicit = zero
392 : term_wp2rtp_explicit = zero
393 : coef_wp2thlp_implicit = zero
394 : term_wp2thlp_explicit = zero
395 :
396 : endif ! .not. l_explicit_turbulent_adv_wpxp
397 :
398 : ! Pack the implicit coefficients and explicit terms into a single type
399 : ! variable for output.
400 0 : pdf_implicit_coefs_terms%coef_wp4_implicit(i,:) = coef_wp4_implicit
401 0 : pdf_implicit_coefs_terms%coef_wprtp2_implicit(i,:) = coef_wprtp2_implicit
402 0 : pdf_implicit_coefs_terms%coef_wpthlp2_implicit(i,:) = coef_wpthlp2_implicit
403 0 : pdf_implicit_coefs_terms%coef_wp2rtp_implicit(i,:) = coef_wp2rtp_implicit
404 0 : pdf_implicit_coefs_terms%term_wp2rtp_explicit(i,:) = term_wp2rtp_explicit
405 0 : pdf_implicit_coefs_terms%coef_wp2thlp_implicit(i,:) = coef_wp2thlp_implicit
406 0 : pdf_implicit_coefs_terms%term_wp2thlp_explicit(i,:) = term_wp2thlp_explicit
407 0 : pdf_implicit_coefs_terms%coef_wprtpthlp_implicit(i,:) = coef_wprtpthlp_implicit
408 0 : pdf_implicit_coefs_terms%term_wprtpthlp_explicit(i,:) = term_wprtpthlp_explicit
409 :
410 : end do
411 :
412 0 : return
413 :
414 : end subroutine new_pdf_driver
415 :
416 : !=============================================================================
417 0 : subroutine calc_responder_var( nz, xm, xp2, sgn_wpxp, mixt_frac, & ! In
418 : coef_spread_DG_means_x, & ! In
419 0 : exp_factor_x, & ! In
420 0 : max_Skx2_pos_Skx_sgn_wpxp, & ! In
421 0 : max_Skx2_neg_Skx_sgn_wpxp, & ! In
422 0 : Skx, & ! In/Out
423 0 : mu_x_1, mu_x_2, & ! Out
424 0 : sigma_x_1_sqd, sigma_x_2_sqd, & ! Out
425 0 : coef_sigma_x_1_sqd, & ! Out
426 0 : coef_sigma_x_2_sqd, & ! Out
427 0 : F_x, min_F_x, max_F_x ) ! Out
428 :
429 : ! Description:
430 : ! This is the sub-driver for a responder variable. The upper limits of the
431 : ! magnitude of Skx are calculated, and Skx is clipped when its magnitude
432 : ! exceeds the upper limits. The limits of the F_x parameter are calculated,
433 : ! and the value of F_x is set within those limits. Then, the PDF parameters
434 : ! for responder variable x are calculated.
435 :
436 : ! References:
437 : !-----------------------------------------------------------------------
438 :
439 : use constants_clubb, only: &
440 : zero ! Variable(s)
441 :
442 : use new_pdf, only: &
443 : calc_limits_F_x_responder, & ! Procedure(s)
444 : calc_responder_params
445 :
446 : use clubb_precision, only: &
447 : core_rknd ! Variable(s)
448 :
449 : implicit none
450 :
451 : integer, intent(in) :: &
452 : nz
453 :
454 : ! Input Variables
455 : real( kind = core_rknd ), dimension(nz), intent(in) :: &
456 : xm, & ! Mean of x (overall) [units vary]
457 : xp2, & ! Variance of x (overall) [(units vary)^2]
458 : sgn_wpxp, & ! Sign of the covariance of w and x [-]
459 : mixt_frac, & ! Mixture fraction [-]
460 : exp_factor_x ! Factor of the form 1 - exp{}; reduces F_x [-]
461 :
462 : real( kind = core_rknd ), intent(in) :: &
463 : coef_spread_DG_means_x ! Coef.: spread betw. PDF comp. means of x [-]
464 :
465 : real( kind = core_rknd ), dimension(nz), intent(in) :: &
466 : max_Skx2_pos_Skx_sgn_wpxp, & ! Maximum Skx^2 when Skx*sgn(<w'x'>) >= 0 [-]
467 : max_Skx2_neg_Skx_sgn_wpxp ! Maximum Skx^2 when Skx*sgn(<w'x'>) < 0 [-]
468 :
469 : ! Input/Output Variable
470 : real( kind = core_rknd ), dimension(nz), intent(inout) :: &
471 : Skx ! Skewness of x (overall) [-]
472 :
473 : ! Output Variables
474 : real( kind = core_rknd ), dimension(nz), intent(out) :: &
475 : mu_x_1, & ! Mean of x (1st PDF component) [units vary]
476 : mu_x_2, & ! Mean of x (2nd PDF component) [units vary]
477 : sigma_x_1_sqd, & ! Variance of x (1st PDF component) [(units vary)^2]
478 : sigma_x_2_sqd ! Variance of x (2nd PDF component) [(units vary)^2]
479 :
480 : real( kind = core_rknd ), dimension(nz), intent(out) :: &
481 : coef_sigma_x_1_sqd, & ! sigma_x_1^2 = coef_sigma_x_1_sqd * <x'^2> [-]
482 : coef_sigma_x_2_sqd ! sigma_x_2^2 = coef_sigma_x_2_sqd * <x'^2> [-]
483 :
484 : ! Output only for recording statistics.
485 : real( kind = core_rknd ), dimension(nz), intent(out) :: &
486 : F_x, & ! Param. for the spread betw. the PDF component means of x [-]
487 : min_F_x, & ! Minimum allowable value of parameter F_x [-]
488 : max_F_x ! Maximum allowable value of parameter F_x [-]
489 :
490 :
491 : ! Calculate the upper limit of the magnitude of Skx.
492 0 : where ( Skx * sgn_wpxp >= zero )
493 : where ( Skx**2 >= max_Skx2_pos_Skx_sgn_wpxp )
494 : where ( Skx >= zero )
495 : Skx = sqrt( 0.99_core_rknd * max_Skx2_pos_Skx_sgn_wpxp )
496 : elsewhere
497 : Skx = -sqrt( 0.99_core_rknd * max_Skx2_pos_Skx_sgn_wpxp )
498 : endwhere
499 : endwhere ! Skx^2 >= max_Skx2_pos_Skx_sgn_wpxp
500 : elsewhere ! Skx * sgn( <w'x'> ) < 0
501 : where ( Skx**2 >= max_Skx2_neg_Skx_sgn_wpxp )
502 : where ( Skx >= zero )
503 : Skx = sqrt( 0.99_core_rknd * max_Skx2_neg_Skx_sgn_wpxp )
504 : elsewhere
505 : Skx = -sqrt( 0.99_core_rknd * max_Skx2_neg_Skx_sgn_wpxp )
506 : endwhere
507 : endwhere ! Skx^2 >= max_Skx2_neg_Skx_sgn_wpxp
508 : endwhere ! Skx * sgn( <w'x'> ) >= 0
509 :
510 : call calc_limits_F_x_responder( nz, mixt_frac, Skx, sgn_wpxp, & ! In
511 : max_Skx2_pos_Skx_sgn_wpxp, & ! In
512 : max_Skx2_neg_Skx_sgn_wpxp, & ! In
513 0 : min_F_x, max_F_x ) ! Out
514 :
515 : ! F_x must have a value between min_F_x and max_F_x.
516 : F_x = calc_F_x_responder( nz, coef_spread_DG_means_x, exp_factor_x, &
517 0 : min_F_x, max_F_x )
518 :
519 : call calc_responder_params( nz, xm, xp2, Skx, sgn_wpxp, & ! In
520 : F_x, mixt_frac, & ! In
521 : mu_x_1, mu_x_2, & ! Out
522 : sigma_x_1_sqd, sigma_x_2_sqd, & ! Out
523 : coef_sigma_x_1_sqd, & ! Out
524 0 : coef_sigma_x_2_sqd ) ! Out
525 :
526 :
527 0 : return
528 :
529 : end subroutine calc_responder_var
530 :
531 : !=============================================================================
532 0 : subroutine calc_F_x_zeta_x_setter( nz, Skx, & ! In
533 : slope_coef_spread_DG_means_x, & ! In
534 : pdf_component_stdev_factor_x, & ! In
535 : lambda, & ! In
536 0 : F_x, zeta_x, & ! Out
537 0 : min_F_x, max_F_x ) ! Out
538 :
539 : ! Description:
540 : ! Calculates the values of F_x and zeta_x for the setter variable (which is
541 : ! the variable that sets the mixture fraction).
542 : !
543 : ! The value of F_x is calculated between 0 (min_F_x) and 1 (max_F_x). The
544 : ! equation is:
545 : !
546 : ! F_x = max_F_x + ( min_F_x - max_F_x )
547 : ! * exp{ -|Skx|^lambda / slope_coef_spread_DG_means_x };
548 : !
549 : ! which reduces to:
550 : !
551 : ! F_x = 1 - exp{ -|Skx|^lambda / slope_coef_spread_DG_means_x };
552 : !
553 : ! where lambda > 0 and slope_coef_spread_DG_means_x > 0. As |Skx| goes
554 : ! toward 0, the value of F_x goes toward 0, and as |Skx| becomes large, the
555 : ! value of F_x goes toward 1. When slope_coef_spread_DG_means_x is small,
556 : ! the value of F_x tends toward 1, and when slope_coef_spread_DG_means_x is
557 : ! large, the value of F_x tends toward 0. When lambda is small, the value
558 : ! of F_x is less dependent on Skx, and when lambda is large, the value of
559 : ! F_x is more dependent on Skx.
560 : !
561 : ! Mathematically, this equation will always produce a value of F_x that
562 : ! falls between min_F_x and max_F_x. However, in order to prevent a value
563 : ! of F_x from being calculated outside the bounds of min_F_x and max_F_x
564 : ! owing to numerical underflow or loss of precision, this equation can be
565 : ! rewritten as:
566 : !
567 : ! F_x
568 : ! = min_F_x * exp{ -|Skx|^lambda / slope_coef_spread_DG_means_x }
569 : ! + max_F_x * ( 1 - exp{ -|Skx|^lambda / slope_coef_spread_DG_means_x } ).
570 : !
571 : ! The value of zeta_x used to adjust the PDF component standard devations:
572 : !
573 : ! 1 + zeta_x = ( mixt_frac * sigma_x_1^2 )
574 : ! / ( ( 1 - mixt_frac ) * sigma_x_2^2 );
575 : !
576 : ! where zeta_x > -1. The sign of zeta_x is used to easily determine if
577 : ! mixt_frac * sigma_x_1^2 is greater than ( 1 - mixt_frac ) * sigma_x_2^2
578 : ! (when zeta_x is positive), mixt_frac * sigma_x_1^2 is less than
579 : ! ( 1 - mixt_frac ) * sigma_x_2^2 (when zeta_x is negative), or
580 : ! mixt_frac * sigma_x_1^2 is equal to ( 1 - mixt_frac ) * sigma_x_2^2 (when
581 : ! zeta_x is 0).
582 : !
583 : ! In order to allow for a tunable parameter that is the pure ratio of
584 : ! mixt_frac * sigma_x_1^2 to ( 1 - mixt_frac ) * sigma_x_2^2, zeta_x is
585 : ! related to the parameter pdf_component_stdev_factor_x, where:
586 : !
587 : ! 1 + zeta_x = pdf_component_stdev_factor_x.
588 :
589 : ! References:
590 : !-----------------------------------------------------------------------
591 :
592 : use constants_clubb, only: &
593 : one, & ! Variable(s)
594 : zero
595 :
596 : use clubb_precision, only: &
597 : core_rknd ! Variable(s)
598 :
599 : implicit none
600 :
601 : integer, intent(in) :: &
602 : nz
603 :
604 : ! Input Variables
605 : real( kind = core_rknd ), dimension(nz), intent(in) :: &
606 : Skx ! Skewness of x (overall) [-]
607 :
608 : real( kind = core_rknd ), intent(in) :: &
609 : slope_coef_spread_DG_means_x, & ! Slope coef: spread PDF comp. means x [-]
610 : pdf_component_stdev_factor_x, & ! Param.: PDF comp. standard devs.; x [-]
611 : lambda ! Param. for Skx dependence [-]
612 :
613 : ! Output Variables
614 : real( kind = core_rknd ), dimension(nz), intent(out) :: &
615 : F_x, & ! Parameter for the spread of the PDF component means of x [-]
616 : zeta_x, & ! Parameter for the PDF component variances of x [-]
617 : min_F_x, & ! Minimum allowable value of parameter F_x [-]
618 : max_F_x ! Maximum allowable value of parameter F_x [-]
619 :
620 : ! Local Variable
621 : real( kind = core_rknd ), dimension(nz) :: &
622 0 : exp_Skx_interp_factor ! Function to interp. between min. and max. [-]
623 :
624 :
625 : ! Set min_F_x to 0 and max_F_x to 1 for the setter variable.
626 0 : where ( abs( Skx ) > zero )
627 : min_F_x = 1.0e-3_core_rknd
628 : elsewhere
629 : min_F_x = zero
630 : endwhere
631 0 : max_F_x = one
632 :
633 : ! F_x must have a value between min_F_x and max_F_x.
634 : exp_Skx_interp_factor &
635 0 : = exp( -abs(Skx)**lambda / slope_coef_spread_DG_means_x )
636 :
637 : F_x = min_F_x * exp_Skx_interp_factor &
638 0 : + max_F_x * ( one - exp_Skx_interp_factor )
639 :
640 : ! The value of zeta_x must be greater than -1.
641 0 : zeta_x = pdf_component_stdev_factor_x - one
642 :
643 :
644 0 : return
645 :
646 : end subroutine calc_F_x_zeta_x_setter
647 :
648 : !=============================================================================
649 0 : function calc_F_x_responder( nz, coef_spread_DG_means_x, exp_factor_x, &
650 0 : min_F_x, max_F_x ) &
651 0 : result( F_x )
652 :
653 : ! Description:
654 : ! Calculates the value of F_x as a tunable function between min_F_x and
655 : ! max_F_x.
656 : !
657 : ! The value of F_x is calculated between min_F_x and max_F_x. The equation
658 : ! is:
659 : !
660 : ! F_x = min_F_x + ( max_F_x - min_F_x )
661 : ! * coef_spread_DG_means_x * exp_factor_x;
662 : !
663 : ! where 0 <= coef_spread_DG_means_x <= 1. As coef_spread_DG_means_x
664 : ! goes toward 0, the value of F_x goes toward min_F_x, and as
665 : ! coef_spread_DG_means_x goes toward 1, the value of F_x goes toward
666 : ! max_F_x. The exp_factor_x is a factor of the form 1 - exp{ }. The range
667 : ! of values of exp_factor_x is 0 <= exp_factor_x <= 1. Here, exp_factor_x
668 : ! is used to reduce the value of F_x under special conditions.
669 : !
670 : ! The equation for exp_factor_x depends on which responder variable is being
671 : ! solved for. For rt, the F_rt equation is:
672 : !
673 : ! F_rt = min_F_rt + ( max_F_rt - min_F_rt )
674 : ! * coef_spread_DG_means_rt * exp_factor_rt;
675 : !
676 : ! where exp_factor_rt is given by:
677 : !
678 : ! exp_factor_rt = 1 - exp{ -0.2 * ( adj_corr_rt_thl + 1 )^5 };
679 : !
680 : ! and where the adjusted (overall) correlation of rt and theta-l, denoted
681 : ! adj_corr_rt_thl, is given by:
682 : !
683 : ! adj_corr_rt_thl = <rt'thl'> / ( sqrt( <rt'^2> ) * sqrt( <thl'^2> ) )
684 : ! * sgn( <w'rt'> ) * sgn( <w'thl'> ).
685 : !
686 : ! The range of values of the adjusted (overall) correlation of rt and
687 : ! theta-l is -1 <= adj_corr_rt_thl <= 1. The adjusted (overall) correlation
688 : ! of rt and theta-l has the same absolute magnitude as the (overall)
689 : ! correlation of rt and theta-l, but the sign of the adjusted correlation is
690 : ! dependent on the agreement between the signs of <w'rt'> and <w'thl'>.
691 : ! When all three covariances are in sign agreement, which is when:
692 : !
693 : ! sgn( <w'rt'> ) * sgn( <w'thl'> ) * sgn( <rt'thl'> ) = 1,
694 : !
695 : ! the value of adj_corr_rt_thl is positive. However, when the three
696 : ! covariances aren't consistent in sign, which is when:
697 : !
698 : ! sgn( <w'rt'> ) * sgn( <w'thl'> ) * sgn( <rt'thl'> ) = -1,
699 : !
700 : ! the value of adj_corr_rt_thl is negative.
701 : !
702 : ! For theta-l, the F_thl equation is:
703 : !
704 : ! F_thl = min_F_thl + ( max_F_thl - min_F_thl ) * coef_spread_DG_means_thl;
705 : !
706 : ! where exp_factor_thl is set to 1 (1 - exp{-inf} = 1) because reducing
707 : ! F_thl through the use of exp_factor_thl is not desired for theta-l.
708 : !
709 : ! The direction of the two PDF component means of x (mu_x_1 and mu_x_2) with
710 : ! respect to each other and the overall mean of x (<x>) is given by
711 : ! sgn( <w'x'> ). When sgn( <w'x'> ) = 1, ( mu_x_1 - <x> ) >= 0 and
712 : ! ( mu_x_2 - <x> ) <= 0. When sgn( <w'x'> ) = -1, ( mu_x_1 - <x> ) <= 0 and
713 : ! ( mu_x_2 - <x> ) >= 0. This helps to promote realizability of the PDF
714 : ! component correlations (corr_w_x_1 and corr_w_x_2).
715 : !
716 : ! The realizability of the PDF component correlations corr_w_rt_1,
717 : ! corr_w_rt_2, corr_w_thl_1, and corr_w_thl_2 is promoted by using
718 : ! sgn( <w'rt'> ) and sgn( <w'thl'> ), respectively to choose the direction
719 : ! of the PDF component means of rt and theta-l. However, when the three
720 : ! covariances (<w'rt'>, <w'thl'>, and <rt'thl'>) aren't consistent in sign
721 : ! (as shown above), it becomes difficult to keep the PDF component
722 : ! correlations corr_rt_thl_1 and corr_rt_thl_2 in the realizable range.
723 : ! However, in this situation, the above equation for exp_factor_rt brings
724 : ! F_rt closer to min_F_rt, which promotes realizability for corr_rt_thl_1
725 : ! and corr_rt_thl_2.
726 : !
727 : ! Mathematically, this equation will always produce a value of F_x that
728 : ! falls between min_F_x and max_F_x. However, in order to prevent a value
729 : ! of F_x from being calculated outside the bounds of min_F_x and max_F_x
730 : ! owing to numerical underflow or loss of precision, this equation can be
731 : ! rewritten as:
732 : !
733 : ! F_x = min_F_x * ( 1 - coef_spread_DG_means_x * exp_factor_x )
734 : ! + max_F_x * coef_spread_DG_means_x * exp_factor_x.
735 :
736 : ! References:
737 : !-----------------------------------------------------------------------
738 :
739 : use constants_clubb, only: &
740 : one ! Variable(s)
741 :
742 : use clubb_precision, only: &
743 : core_rknd ! Variable(s)
744 :
745 : implicit none
746 :
747 : integer, intent(in) :: &
748 : nz
749 :
750 : ! Input Variables
751 : real( kind = core_rknd ), intent(in) :: &
752 : coef_spread_DG_means_x ! Coef.: spread betw. PDF comp. means of x [-]
753 :
754 : real( kind = core_rknd ), dimension(nz), intent(in) :: &
755 : exp_factor_x, & ! Factor of the form 1 - exp{}; reduces F_x [-]
756 : min_F_x, & ! Minimum allowable value of parameter F_x [-]
757 : max_F_x ! Maximum allowable value of parameter F_x [-]
758 :
759 : ! Return Variable
760 : real( kind = core_rknd ), dimension(nz) :: &
761 : F_x ! Parameter for the spread between the PDF component means of x [-]
762 :
763 :
764 : ! F_x must have a value between min_F_x and max_F_x.
765 : F_x = min_F_x * ( one - coef_spread_DG_means_x * exp_factor_x ) &
766 0 : + max_F_x * coef_spread_DG_means_x * exp_factor_x
767 :
768 :
769 0 : return
770 :
771 0 : end function calc_F_x_responder
772 :
773 : !=============================================================================
774 :
775 : end module new_pdf_main
|