Line data Source code
1 : ! $Id$
2 : !===============================================================================
3 : module new_hybrid_pdf_main
4 :
5 : ! Description:
6 : ! The portion of CLUBB's multivariate, two-component PDF that is the
7 : ! multivariate, two-component normal PDF of vertical velocity (w), total water
8 : ! mixing ratio (rt), liquid water potential temperature (thl), and optionally,
9 : ! the west-east horizontal wind component (u), the south-north horizontal wind
10 : ! component (v), and passive scalars (sclr).
11 :
12 : ! References:
13 : !-------------------------------------------------------------------------
14 :
15 : implicit none
16 :
17 : public :: new_hybrid_pdf_driver ! Procedure(s)
18 :
19 : private :: calc_responder_driver, & ! Procedure(s)
20 : calc_F_w_zeta_w
21 :
22 : private
23 :
24 : contains
25 :
26 : !=============================================================================
27 0 : subroutine new_hybrid_pdf_driver( nz, ngrdcol, sclr_dim, & ! In
28 0 : wm, rtm, thlm, um, vm, & ! In
29 0 : wp2, rtp2, thlp2, up2, vp2, & ! In
30 0 : Skw, wprtp, wpthlp, upwp, vpwp, & ! In
31 0 : sclrm, sclrp2, wpsclrp, & ! In
32 0 : gamma_Skw_fnc, & ! In
33 0 : slope_coef_spread_DG_means_w, & ! In
34 0 : pdf_component_stdev_factor_w, & ! In
35 0 : Skrt, Skthl, Sku, Skv, Sksclr, & ! I/O
36 0 : mu_w_1, mu_w_2, & ! Out
37 0 : mu_rt_1, mu_rt_2, & ! Out
38 0 : mu_thl_1, mu_thl_2, & ! Out
39 0 : mu_u_1, mu_u_2, mu_v_1, mu_v_2, & ! Out
40 0 : sigma_w_1_sqd, sigma_w_2_sqd, & ! Out
41 0 : sigma_rt_1_sqd, sigma_rt_2_sqd, & ! Out
42 0 : sigma_thl_1_sqd, sigma_thl_2_sqd, & ! Out
43 0 : sigma_u_1_sqd, sigma_u_2_sqd, & ! Out
44 0 : sigma_v_1_sqd, sigma_v_2_sqd, & ! Out
45 0 : mu_sclr_1, mu_sclr_2, & ! Out
46 0 : sigma_sclr_1_sqd, sigma_sclr_2_sqd, & ! Out
47 0 : mixt_frac, & ! Out
48 : pdf_implicit_coefs_terms, & ! Out
49 0 : F_w, min_F_w, max_F_w ) ! Out
50 :
51 :
52 : ! Description:
53 : ! Calculate the PDF parameters for w (including mixture fraction), rt,
54 : ! theta-l, and optionally, u, v, and passive scalar variables.
55 :
56 : ! References:
57 : !-----------------------------------------------------------------------
58 :
59 : use constants_clubb, only: &
60 : zero, & ! Constant(s)
61 : fstderr
62 :
63 : use new_hybrid_pdf, only: &
64 : calculate_w_params, & ! Procedure(s)
65 : calculate_coef_wp4_implicit, &
66 : calc_coef_wp2xp_implicit, &
67 : calc_coefs_wpxp2_semiimpl, &
68 : calc_coefs_wpxpyp_semiimpl
69 :
70 : use pdf_parameter_module, only: &
71 : implicit_coefs_terms ! Variable Type
72 :
73 : use model_flags, only: &
74 : l_explicit_turbulent_adv_wp3, & ! Variable(s)
75 : l_explicit_turbulent_adv_wpxp, &
76 : l_explicit_turbulent_adv_xpyp
77 :
78 : use clubb_precision, only: &
79 : core_rknd ! Variable(s)
80 :
81 : implicit none
82 :
83 : integer, intent(in) :: &
84 : nz, &
85 : ngrdcol, &
86 : sclr_dim
87 :
88 : ! Input Variables
89 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
90 : wm, & ! Mean of w (overall) [m/s]
91 : rtm, & ! Mean of rt (overall) [kg/kg]
92 : thlm, & ! Mean of thl (overall) [K]
93 : um, & ! Mean of u (overall) [m/s]
94 : vm, & ! Mean of v (overall) [m/s]
95 : wp2, & ! Variance of w (overall) [m^2/s^2]
96 : rtp2, & ! Variance of rt (overall) [kg^2/kg^2]
97 : thlp2, & ! Variance of thl (overall) [K^2]
98 : up2, & ! Variance of u (overall) [m^2/s^2]
99 : vp2, & ! Variance of v (overall) [m^2/s^2]
100 : Skw, & ! Skewness of w (overall) [-]
101 : wprtp, & ! Covariance of w and rt (overall) [(m/s)kg/kg]
102 : wpthlp, & ! Covariance of w and thl (overall) [(m/s)K]
103 : upwp, & ! Covariance of u and w (overall) [m^2/s^2]
104 : vpwp ! Covariance of v and w (overall) [m^2/s^2]
105 :
106 : real( kind = core_rknd ), dimension(ngrdcol,nz,sclr_dim), intent(in) :: &
107 : sclrm, & ! Mean of sclr (overall) [units vary]
108 : sclrp2, & ! Variance of sclr (overall) [(units vary)^2]
109 : wpsclrp ! Covariance of w and sclr (overall) [(m/s)(units vary)]
110 :
111 : ! Tunable parameter gamma.
112 : ! When gamma goes to 0, the standard deviations of w in each PDF component
113 : ! become small, and the spread between the two PDF component means of w
114 : ! becomes large. F_w goes to min_F_w.
115 : ! When gamma goes to 1, the standard deviations of w in each PDF component
116 : ! become large, and the spread between the two PDF component means of w
117 : ! becomes small. F_w goes to max_F_w.
118 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
119 : gamma_Skw_fnc ! Value of parameter gamma from tunable Skw function [-]
120 :
121 : real( kind = core_rknd ), dimension(ngrdcol), intent(in) :: &
122 : ! Slope coefficient for the spread between the PDF component means of w.
123 : slope_coef_spread_DG_means_w, &
124 : ! Parameter to adjust the PDF component standard deviations of w.
125 : pdf_component_stdev_factor_w
126 :
127 : ! Input/Output Variables
128 : ! These variables are input/output because their values may be clipped.
129 : ! Otherwise, as long as it is not necessary to clip them, their values
130 : ! will stay the same.
131 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(inout) :: &
132 : Skrt, & ! Skewness of rt (overall) [-]
133 : Skthl, & ! Skewness of thl (overall) [-]
134 : Sku, & ! Skewness of u (overall) [-]
135 : Skv ! Skewness of v (overall) [-]
136 :
137 : real( kind = core_rknd ), dimension(ngrdcol,nz,sclr_dim), intent(inout) :: &
138 : Sksclr ! Skewness of sclr (overall) [-]
139 :
140 : ! Output Variables
141 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
142 : mu_w_1, & ! Mean of w (1st PDF component) [m/s]
143 : mu_w_2, & ! Mean of w (2nd PDF component) [m/s]
144 : mu_rt_1, & ! Mean of rt (1st PDF component) [kg/kg]
145 : mu_rt_2, & ! Mean of rt (2nd PDF component) [kg/kg]
146 : mu_thl_1, & ! Mean of thl (1st PDF component) [K]
147 : mu_thl_2, & ! Mean of thl (2nd PDF component) [K]
148 : mu_u_1, & ! Mean of u (1st PDF component) [m/s]
149 : mu_u_2, & ! Mean of u (2nd PDF component) [m/s]
150 : mu_v_1, & ! Mean of v (1st PDF component) [m/s]
151 : mu_v_2, & ! Mean of v (2nd PDF component) [m/s]
152 : sigma_w_1_sqd, & ! Variance of w (1st PDF component) [m^2/s^2]
153 : sigma_w_2_sqd, & ! Variance of w (2nd PDF component) [m^2/s^2]
154 : sigma_rt_1_sqd, & ! Variance of rt (1st PDF component) [kg^2/kg^2]
155 : sigma_rt_2_sqd, & ! Variance of rt (2nd PDF component) [kg^2/kg^2]
156 : sigma_thl_1_sqd, & ! Variance of thl (1st PDF component) [K^2]
157 : sigma_thl_2_sqd, & ! Variance of thl (2nd PDF component) [K^2]
158 : sigma_u_1_sqd, & ! Variance of u (1st PDF component) [m^2/s^2]
159 : sigma_u_2_sqd, & ! Variance of u (2nd PDF component) [m^2/s^2]
160 : sigma_v_1_sqd, & ! Variance of v (1st PDF component) [m^2/s^2]
161 : sigma_v_2_sqd, & ! Variance of v (2nd PDF component) [m^2/s^2]
162 : mixt_frac ! Mixture fraction [-]
163 :
164 : real( kind = core_rknd ), dimension(ngrdcol,nz,sclr_dim), intent(out) :: &
165 : mu_sclr_1, & ! Mean of sclr (1st PDF component) [units vary]
166 : mu_sclr_2, & ! Mean of sclr (2nd PDF component) [units vary]
167 : sigma_sclr_1_sqd, & ! Variance of sclr (1st PDF component) [(un. vary)^2]
168 : sigma_sclr_2_sqd ! Variance of sclr (2nd PDF component) [(un. vary)^2]
169 :
170 : type(implicit_coefs_terms), intent(inout) :: &
171 : pdf_implicit_coefs_terms ! Implicit coefs / explicit terms [units vary]
172 :
173 : ! Output only for recording statistics.
174 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
175 : F_w, & ! Parameter for the spread of the PDF component means of w [-]
176 : min_F_w, & ! Minimum allowable value of parameter F_w [-]
177 : max_F_w ! Maximum allowable value of parameter F_w [-]
178 :
179 : ! Local Variables
180 : real( kind = core_rknd ), dimension(nz) :: &
181 0 : sigma_w_1, & ! Standard deviation of w (1st PDF component) [m/s]
182 0 : sigma_w_2 ! Standard deviation of w (2nd PDF component) [m/s]
183 :
184 : real( kind = core_rknd ), dimension(nz) :: &
185 0 : coef_sigma_w_1_sqd, & ! sigma_w_1^2 = coef_sigma_w_1_sqd * <w'^2> [-]
186 0 : coef_sigma_w_2_sqd, & ! sigma_w_2^2 = coef_sigma_w_2_sqd * <w'^2> [-]
187 0 : coef_sigma_rt_1_sqd, & ! sigma_rt_1^2 = coef_sigma_rt_1_sqd * <rt'^2> [-]
188 0 : coef_sigma_rt_2_sqd, & ! sigma_rt_2^2 = coef_sigma_rt_2_sqd * <rt'^2> [-]
189 0 : coef_sigma_thl_1_sqd, & ! sigma_thl_1^2=coef_sigma_thl_1_sqd*<thl'^2> [-]
190 0 : coef_sigma_thl_2_sqd, & ! sigma_thl_2^2=coef_sigma_thl_2_sqd*<thl'^2> [-]
191 0 : coef_sigma_u_1_sqd, & ! sigma_u_1^2 = coef_sigma_u_1_sqd * <u'^2> [-]
192 0 : coef_sigma_u_2_sqd, & ! sigma_u_2^2 = coef_sigma_u_2_sqd * <u'^2> [-]
193 0 : coef_sigma_v_1_sqd, & ! sigma_v_1^2 = coef_sigma_v_1_sqd * <v'^2> [-]
194 0 : coef_sigma_v_2_sqd ! sigma_v_2^2 = coef_sigma_v_2_sqd * <v'^2> [-]
195 :
196 : ! sigma_sclr_1^2 = coef_sigma_sclr_1_sqd * <sclr'^2>
197 : ! sigma_sclr_2^2 = coef_sigma_sclr_2_sqd * <sclr'^2>
198 : real( kind = core_rknd ), dimension(nz,sclr_dim) :: &
199 0 : coef_sigma_sclr_1_sqd, & ! Coefficient that is multiplied by <sclr'^2> [-]
200 0 : coef_sigma_sclr_2_sqd ! Coefficient that is multiplied by <sclr'^2> [-]
201 :
202 : real( kind = core_rknd ), dimension(nz) :: &
203 0 : zeta_w ! Parameter for the PDF component variances of w [-]
204 :
205 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
206 : ! Slope coefficient for the spread between the PDF component means of w.
207 0 : slope_coef_spread_DG_means_w_in, &
208 : ! Parameter to adjust the PDF component standard deviations of w.
209 0 : pdf_component_stdev_factor_w_in
210 :
211 : real( kind = core_rknd ), dimension(nz) :: &
212 0 : coef_wp4_implicit, & ! <w'^4> = coef_wp4_implicit * <w'^2>^2 [-]
213 0 : coef_wp2rtp_implicit, & ! <w'^2 rt'> = coef_wp2rtp_implicit*<w'rt'> [m/s]
214 0 : coef_wp2thlp_implicit, & ! <w'^2 thl'>=coef_wp2thlp_implicit*<w'thl'>[m/s]
215 0 : coef_wp2up_implicit, & ! <w'^2 u'> = coef_wp2up_implicit * <u'w'> [m/s]
216 0 : coef_wp2vp_implicit ! <w'^2 v'> = coef_wp2vp_implicit * <v'w'> [m/s]
217 :
218 : ! <w'^2 sclr'> = coef_wp2sclrp_implicit * <w'sclr'>
219 : real( kind = core_rknd ), dimension(nz,sclr_dim) :: &
220 0 : coef_wp2sclrp_implicit ! Coef. that is multiplied by <w'sclr'> [m/s]
221 :
222 : ! <w'rt'^2> = coef_wprtp2_implicit * <rt'^2> + term_wprtp2_explicit
223 : real( kind = core_rknd ), dimension(nz) :: &
224 0 : coef_wprtp2_implicit, & ! Coefficient that is multiplied by <rt'^2> [m/s]
225 0 : term_wprtp2_explicit ! Term that is on the RHS [m/s kg^2/kg^2]
226 :
227 : ! <w'thl'^2> = coef_wpthlp2_implicit * <thl'^2> + term_wpthlp2_explicit
228 : real( kind = core_rknd ), dimension(nz) :: &
229 0 : coef_wpthlp2_implicit, & ! Coef. that is multiplied by <thl'^2> [m/s]
230 0 : term_wpthlp2_explicit ! Term that is on the RHS [m/s K^2]
231 :
232 : ! <w'rt'thl'> = coef_wprtpthlp_implicit*<rt'thl'> + term_wprtpthlp_explicit
233 : real( kind = core_rknd ), dimension(nz) :: &
234 0 : coef_wprtpthlp_implicit, & ! Coef. that is multiplied by <rt'thl'> [m/s]
235 0 : term_wprtpthlp_explicit ! Term that is on the RHS [m/s(kg/kg)K]
236 :
237 : ! <w'u'^2> = coef_wpup2_implicit * <u'^2> + term_wpup2_explicit
238 : real( kind = core_rknd ), dimension(nz) :: &
239 0 : coef_wpup2_implicit, & ! Coefficient that is multiplied by <u'^2> [m/s]
240 0 : term_wpup2_explicit ! Term that is on the RHS [m^3/s^3]
241 :
242 : ! <w'v'^2> = coef_wpvp2_implicit * <v'^2> + term_wpvp2_explicit
243 : real( kind = core_rknd ), dimension(nz) :: &
244 0 : coef_wpvp2_implicit, & ! Coefficient that is multiplied by <v'^2> [m/s]
245 0 : term_wpvp2_explicit ! Term that is on the RHS [m^3/s^3]
246 :
247 : ! <w'sclr'^2> = coef_wpsclrp2_implicit * <sclr'^2> + term_wpsclrp2_explicit
248 : real( kind = core_rknd ), dimension(nz,sclr_dim) :: &
249 0 : coef_wpsclrp2_implicit, & ! Coef. that is multiplied by <sclr'^2> [m/s]
250 0 : term_wpsclrp2_explicit ! Term that is on the RHS [m/s(units vary)^2]
251 :
252 : ! <w'rt'sclr'> = coef_wprtpsclrp_implicit * <sclr'rt'>
253 : ! + term_wprtpsclrp_explicit
254 : real( kind = core_rknd ), dimension(nz,sclr_dim) :: &
255 0 : coef_wprtpsclrp_implicit, & ! Coef. that is multiplied by <sclr'rt'> [m/s]
256 0 : term_wprtpsclrp_explicit ! Term that is on the RHS [m/s(kg/kg)(un. v.)]
257 :
258 : ! <w'thl'sclr'> = coef_wpthlpsclrp_implicit * <sclr'thl'>
259 : ! + term_wpthlpsclrp_explicit
260 : real( kind = core_rknd ), dimension(nz,sclr_dim) :: &
261 0 : coef_wpthlpsclrp_implicit, & ! Coef. that is mult. by <sclr'thl'> [m/s]
262 0 : term_wpthlpsclrp_explicit ! Term that is on the RHS [(m/s)K(un. vary)]
263 :
264 : ! Variables to interface with code for the jth scalar
265 : real( kind = core_rknd ), dimension(nz) :: &
266 0 : sclrjm, & ! Mean of sclr j (overall) [units vary]
267 0 : sclrjp2, & ! Variance of sclr j (overall) [(units vary)^2]
268 0 : wpsclrjp, & ! Covariance of w and sclr j (overall) [(m/s)(units vary)]
269 0 : Sksclrj ! Skewness of rt (overall) [-]
270 :
271 : real( kind = core_rknd ), dimension(nz) :: &
272 0 : mu_sclrj_1, & ! Mean of sclr j (1st PDF component) [units vary]
273 0 : mu_sclrj_2, & ! Mean of sclr j (2nd PDF component) [units vary]
274 0 : sigma_sclrj_1_sqd, & ! Variance of sclr j (1st PDF comp.) [(units vary)^2]
275 0 : sigma_sclrj_2_sqd ! Variance of sclr j (2nd PDF comp.) [(units vary)^2]
276 :
277 : ! sigma_sclrj_1^2 = coef_sigma_sclrj_1_sqd * <sclrj'^2>
278 : ! sigma_sclrj_2^2 = coef_sigma_sclrj_2_sqd * <sclrj'^2>
279 : real( kind = core_rknd ), dimension(nz) :: &
280 0 : coef_sigma_sclrj_1_sqd, & ! Coef. that is multiplied by <sclrj'^2> [-]
281 0 : coef_sigma_sclrj_2_sqd ! Coef. that is multiplied by <sclrj'^2> [-]
282 :
283 : ! <w'sclrj'^2> = coef_wpsclrjp2_implicit * <sclrj'^2>
284 : ! + term_wpsclrjp2_explicit
285 : real( kind = core_rknd ), dimension(nz) :: &
286 0 : coef_wpsclrjp2_implicit, & ! Coef. that is multiplied by <sclrj'^2> [m/s]
287 0 : term_wpsclrjp2_explicit ! Term that is on the RHS [m/s(units vary)^2]
288 :
289 : ! <w'rt'sclrj'> = coef_wprtpsclrjp_implicit * <sclrj'rt'>
290 : ! + term_wprtpsclrjp_explicit
291 : real( kind = core_rknd ), dimension(nz) :: &
292 0 : coef_wprtpsclrjp_implicit, & ! Coef. that is mult. by <sclr'rt'> [m/s]
293 0 : term_wprtpsclrjp_explicit ! Term that is on the RHS [m/s(kg/kg)(un.v.)]
294 :
295 : ! <w'thl'sclrj'> = coef_wpthlpsclrjp_implicit * <sclrj'thl'>
296 : ! + term_wpthlpsclrjp_explicit
297 : real( kind = core_rknd ), dimension(nz) :: &
298 0 : coef_wpthlpsclrjp_implicit, & ! Coef. that is mult. by <sclrj'thl'> [m/s]
299 0 : term_wpthlpsclrjp_explicit ! Term that is on the RHS [(m/s)K(un. vary)]
300 :
301 : real( kind = core_rknd ), dimension(nz) :: &
302 0 : max_corr_w_sclr_sqd ! Max value of wpsclrp^2 / ( wp2 * sclrp2 ) [-]
303 :
304 : real( kind = core_rknd ), dimension(nz) :: &
305 0 : zeros ! Vector of 0s (size nz) [-]
306 :
307 : real( kind = core_rknd ), dimension(nz,sclr_dim) :: &
308 0 : zero_array ! Array of 0s (size nz x sclr_dim) [-]
309 :
310 : integer :: i, k, j ! Loop indices
311 :
312 0 : do i = 1, ngrdcol
313 :
314 : ! Calculate the maximum value of the square of the correlation of w and a
315 : ! scalar when scalars are used.
316 : ! Initialize max_corr_w_sclr_sqd to 0. It needs to retain this value even
317 : ! when sclr_dim = 0.
318 0 : max_corr_w_sclr_sqd = zero
319 0 : if ( sclr_dim > 0 ) then
320 0 : do k = 1, nz, 1
321 0 : do j = 1, sclr_dim, 1
322 0 : if ( wp2(i,k) * sclrp2(i,k,j) > zero ) then
323 : max_corr_w_sclr_sqd(k) = max( wpsclrp(i,k,j)**2 &
324 : / ( wp2(i,k) * sclrp2(i,k,j) ), &
325 0 : max_corr_w_sclr_sqd(k) )
326 : endif ! wp2(k) * sclrp2(k,j) > 0
327 : enddo ! j = 1, sclr_dim, 1
328 : enddo ! k = 1, nz, 1
329 : endif ! sclr_dim > 0
330 :
331 0 : slope_coef_spread_DG_means_w_in(i,:) = slope_coef_spread_DG_means_w(i)
332 0 : pdf_component_stdev_factor_w_in(i,:) = pdf_component_stdev_factor_w(i)
333 :
334 : ! Calculate the values of PDF tunable parameters F_w and zeta_w.
335 : ! Vertical velocity, w, will always be the setter variable.
336 : call calc_F_w_zeta_w( Skw(i,:), wprtp(i,:), wpthlp(i,:), upwp(i,:), vpwp(i,:), & ! In
337 : wp2(i,:), rtp2(i,:), thlp2(i,:), up2(i,:), vp2(i,:), & ! In
338 : gamma_Skw_fnc(i,:), & ! In
339 : slope_coef_spread_DG_means_w_in(i,:), & ! In
340 : pdf_component_stdev_factor_w_in(i,:), & ! In
341 : max_corr_w_sclr_sqd, & ! In
342 0 : F_w(i,:), zeta_w, min_F_w(i,:), max_F_w(i,:) ) ! Out
343 :
344 : ! Calculate the PDF parameters, including mixture fraction, for the
345 : ! setter variable, w.
346 : call calculate_w_params( wm(i,:), wp2(i,:), Skw(i,:), F_w(i,:), zeta_w, & ! In
347 : mu_w_1(i,:), mu_w_2(i,:), sigma_w_1, & ! Out
348 : sigma_w_2, mixt_frac(i,:), & ! Out
349 : coef_sigma_w_1_sqd, & ! Out
350 0 : coef_sigma_w_2_sqd ) ! Out
351 :
352 0 : sigma_w_1_sqd(i,:) = sigma_w_1**2
353 0 : sigma_w_2_sqd(i,:) = sigma_w_2**2
354 :
355 0 : if ( any( mixt_frac(i,:) < zero ) ) then
356 0 : write(fstderr,*) "Mixture fraction cannot be calculated."
357 : write(fstderr,*) "The value of F_w must be greater than 0 when " &
358 0 : // "| Skw | > 0."
359 0 : error stop
360 : endif
361 :
362 : ! Calculate the PDF parameters for responder variable rt.
363 : call calc_responder_driver( rtm(i,:), rtp2(i,:), wprtp(i,:), wp2(i,:), & ! In
364 : mixt_frac(i,:), F_w(i,:), & ! In
365 : Skrt(i,:), & ! In/Out
366 : mu_rt_1(i,:), mu_rt_2(i,:), & ! Out
367 : sigma_rt_1_sqd(i,:), & ! Out
368 : sigma_rt_2_sqd(i,:), & ! Out
369 : coef_sigma_rt_1_sqd, & ! Out
370 0 : coef_sigma_rt_2_sqd ) ! Out
371 :
372 : ! Calculate the PDF parameters for responder variable thl.
373 : call calc_responder_driver( thlm(i,:), thlp2(i,:), wpthlp(i,:), wp2(i,:), & ! In
374 : mixt_frac(i,:), F_w(i,:), & ! In
375 : Skthl(i,:), & ! In/Out
376 : mu_thl_1(i,:), mu_thl_2(i,:), & ! Out
377 : sigma_thl_1_sqd(i,:), & ! Out
378 : sigma_thl_2_sqd(i,:), & ! Out
379 : coef_sigma_thl_1_sqd, & ! Out
380 0 : coef_sigma_thl_2_sqd ) ! Out
381 :
382 : ! Calculate the PDF parameters for responder variable u.
383 : call calc_responder_driver( um(i,:), up2(i,:), upwp(i,:), wp2(i,:), & ! In
384 : mixt_frac(i,:), F_w(i,:), & ! In
385 : Sku(i,:), & ! In/Out
386 : mu_u_1(i,:), mu_u_2(i,:), & ! Out
387 : sigma_u_1_sqd(i,:), & ! Out
388 : sigma_u_2_sqd(i,:), & ! Out
389 : coef_sigma_u_1_sqd, & ! Out
390 0 : coef_sigma_u_2_sqd ) ! Out
391 :
392 : ! Calculate the PDF parameters for responder variable v.
393 : call calc_responder_driver( vm(i,:), vp2(i,:), vpwp(i,:), wp2(i,:), & ! In
394 : mixt_frac(i,:), F_w(i,:), & ! In
395 : Skv(i,:), & ! In/Out
396 : mu_v_1(i,:), mu_v_2(i,:), & ! Out
397 : sigma_v_1_sqd(i,:), & ! Out
398 : sigma_v_2_sqd(i,:), & ! Out
399 : coef_sigma_v_1_sqd, & ! Out
400 0 : coef_sigma_v_2_sqd ) ! Out
401 :
402 : ! Calculate the PDF parameters for responder variable sclr.
403 0 : if ( sclr_dim > 0 ) then
404 :
405 0 : do j = 1, sclr_dim, 1
406 :
407 0 : do k = 1, nz, 1
408 0 : sclrjm(k) = sclrm(i,k,j)
409 0 : sclrjp2(k) = sclrp2(i,k,j)
410 0 : wpsclrjp(k) = wpsclrp(i,k,j)
411 0 : Sksclrj(k) = Sksclr(i,k,j)
412 : enddo ! k = 1, nz, 1
413 :
414 : call calc_responder_driver( sclrjm, sclrjp2, wpsclrjp, wp2(i,:), & ! In
415 : mixt_frac(i,:), F_w(i,:), & ! In
416 : Sksclrj, & ! In/Out
417 : mu_sclrj_1, mu_sclrj_2, & ! Out
418 : sigma_sclrj_1_sqd, & ! Out
419 : sigma_sclrj_2_sqd, & ! Out
420 : coef_sigma_sclrj_1_sqd, & ! Out
421 0 : coef_sigma_sclrj_2_sqd ) ! Out
422 :
423 0 : do k = 1, nz, 1
424 0 : Sksclr(i,k,j) = Sksclrj(k)
425 0 : mu_sclr_1(i,k,j) = mu_sclrj_1(k)
426 0 : mu_sclr_2(i,k,j) = mu_sclrj_2(k)
427 0 : sigma_sclr_1_sqd(i,k,j) = sigma_sclrj_1_sqd(k)
428 0 : sigma_sclr_2_sqd(i,k,j) = sigma_sclrj_2_sqd(k)
429 0 : coef_sigma_sclr_1_sqd(k,j) = coef_sigma_sclrj_1_sqd(k)
430 0 : coef_sigma_sclr_2_sqd(k,j) = coef_sigma_sclrj_2_sqd(k)
431 : enddo ! k = 1, nz, 1
432 :
433 : enddo ! j = 1, sclr_dim, 1
434 :
435 : endif ! sclr_dim > 0
436 :
437 :
438 : if ( .not. l_explicit_turbulent_adv_wp3 ) then
439 :
440 : ! Turbulent advection of <w'^3> is being handled implicitly.
441 :
442 : ! <w'^4> = coef_wp4_implicit * <w'^2>^2.
443 : coef_wp4_implicit &
444 : = calculate_coef_wp4_implicit( mixt_frac(i,:), F_w(i,:), &
445 : coef_sigma_w_1_sqd, &
446 0 : coef_sigma_w_2_sqd )
447 :
448 : else ! l_explicit_turbulent_adv_wp3
449 :
450 : ! Turbulent advection of <w'^3> is being handled explicitly.
451 : coef_wp4_implicit = zero
452 :
453 : endif ! .not. l_explicit_turbulent_adv_wp3
454 :
455 : if ( .not. l_explicit_turbulent_adv_wpxp ) then
456 :
457 : ! Turbulent advection of <w'rt'>, <w'thl'>, <u'w'>, <v'w'>, and <w'sclr'>
458 : ! are all being handled implicitly.
459 :
460 : ! <w'^2 rt'> = coef_wp2rtp_implicit * <w'rt'>
461 : coef_wp2rtp_implicit = calc_coef_wp2xp_implicit( wp2(i,:), mixt_frac(i,:), F_w(i,:), &
462 : coef_sigma_w_1_sqd, &
463 0 : coef_sigma_w_2_sqd )
464 :
465 : ! <w'^2 thl'> = coef_wp2thlp_implicit * <w'thl'>;
466 : ! <w'^2 u'> = coef_wp2up_implicit * <u'w'>; and
467 : ! <w'^2 v'> = coef_wp2vp_implicit * <v'w'>;
468 : ! where each coef_wp2xp_implicit is the same as coef_wp2rtp_implicit.
469 0 : coef_wp2thlp_implicit = coef_wp2rtp_implicit
470 0 : coef_wp2up_implicit = coef_wp2rtp_implicit
471 0 : coef_wp2vp_implicit = coef_wp2rtp_implicit
472 :
473 : ! <w'^2 sclr'> = coef_wp2sclrp_implicit * <w'sclr'>;
474 : ! where each coef_wp2xp_implicit is the same as coef_wp2rtp_implicit.
475 0 : if ( sclr_dim > 0 ) then
476 0 : do k = 1, nz, 1
477 0 : do j = 1, sclr_dim, 1
478 0 : coef_wp2sclrp_implicit(k,j) = coef_wp2rtp_implicit(k)
479 : enddo ! j = 1, sclr_dim, 1
480 : enddo ! k = 1, nz, 1
481 : endif ! sclr_dim > 0
482 :
483 : else ! l_explicit_turbulent_adv_wpxp
484 :
485 : ! Turbulent advection of <w'rt'>, <w'thl'>, <u'w'>, <v'w'>, and <w'sclr'>
486 : ! are all being handled explicitly.
487 : coef_wp2rtp_implicit = zero
488 : coef_wp2thlp_implicit = zero
489 : coef_wp2up_implicit = zero
490 : coef_wp2vp_implicit = zero
491 : if ( sclr_dim > 0 ) then
492 : coef_wp2sclrp_implicit = zero
493 : endif ! sclr_dim > 0
494 :
495 : endif ! .not. l_explicit_turbulent_adv_wpxp
496 :
497 : if ( .not. l_explicit_turbulent_adv_xpyp ) then
498 :
499 : ! Turbulent advection of <rt'^2>, <thl'^2>, <rt'thl'>, <u'^2>, <v'^2>,
500 : ! <sclr'^2>, <sclr'rt'>, and <sclr'thl'> are all being handled
501 : ! semi-implicitly.
502 :
503 : ! <w'rt'^2> = coef_wprtp2_implicit * <rt'^2> + term_wprtp2_explicit
504 : call calc_coefs_wpxp2_semiimpl( wp2(i,:), wprtp(i,:), & ! In
505 : mixt_frac(i,:), F_w(i,:), & ! In
506 : coef_sigma_rt_1_sqd, & ! In
507 : coef_sigma_rt_2_sqd, & ! In
508 : coef_wprtp2_implicit, & ! Out
509 0 : term_wprtp2_explicit ) ! Out
510 :
511 : ! <w'thl'^2> = coef_wpthlp2_implicit * <thl'^2> + term_wprtp2_explicit
512 : call calc_coefs_wpxp2_semiimpl( wp2(i,:), wpthlp(i,:), & ! In
513 : mixt_frac(i,:), F_w(i,:), & ! In
514 : coef_sigma_thl_1_sqd, & ! In
515 : coef_sigma_thl_2_sqd, & ! In
516 : coef_wpthlp2_implicit, & ! Out
517 0 : term_wpthlp2_explicit ) ! Out
518 :
519 : ! <w'rt'thl'> = coef_wprtpthlp_implicit * <rt'thl'>
520 : ! + term_wprtpthlp_explicit
521 : call calc_coefs_wpxpyp_semiimpl( wp2(i,:), wprtp(i,:), wpthlp(i,:), & ! In
522 : mixt_frac(i,:), F_w(i,:), & ! In
523 : coef_sigma_rt_1_sqd, & ! In
524 : coef_sigma_rt_2_sqd, & ! In
525 : coef_sigma_thl_1_sqd, & ! In
526 : coef_sigma_thl_2_sqd, & ! In
527 : coef_wprtpthlp_implicit, & ! Out
528 0 : term_wprtpthlp_explicit ) ! Out
529 :
530 : ! <w'u'^2> = coef_wpup2_implicit * <u'^2> + term_wpup2_explicit
531 : call calc_coefs_wpxp2_semiimpl( wp2(i,:), upwp(i,:), & ! In
532 : mixt_frac(i,:), F_w(i,:), & ! In
533 : coef_sigma_u_1_sqd, & ! In
534 : coef_sigma_u_2_sqd, & ! In
535 : coef_wpup2_implicit, & ! Out
536 0 : term_wpup2_explicit ) ! Out
537 :
538 : ! <w'v'^2> = coef_wpvp2_implicit * <v'^2> + term_wpvp2_explicit
539 : call calc_coefs_wpxp2_semiimpl( wp2(i,:), vpwp(i,:), & ! In
540 : mixt_frac(i,:), F_w(i,:), & ! In
541 : coef_sigma_v_1_sqd, & ! In
542 : coef_sigma_v_2_sqd, & ! In
543 : coef_wpvp2_implicit, & ! Out
544 0 : term_wpvp2_explicit ) ! Out
545 :
546 0 : if ( sclr_dim > 0 ) then
547 :
548 0 : do j = 1, sclr_dim, 1
549 :
550 0 : do k = 1, nz, 1
551 0 : wpsclrjp(k) = wpsclrp(i,k,j)
552 0 : coef_sigma_sclrj_1_sqd(k) = coef_sigma_sclr_1_sqd(k,j)
553 0 : coef_sigma_sclrj_2_sqd(k) = coef_sigma_sclr_2_sqd(k,j)
554 : enddo ! k = 1, nz, 1
555 :
556 : ! <w'sclr'^2> = coef_wpsclrp2_implicit * <sclr'^2>
557 : ! + term_wpsclrp2_explicit
558 : call calc_coefs_wpxp2_semiimpl( wp2(i,:), wpsclrjp, & ! In
559 : mixt_frac(i,:), F_w(i,:), & ! In
560 : coef_sigma_sclrj_1_sqd, & ! In
561 : coef_sigma_sclrj_2_sqd, & ! In
562 : coef_wpsclrjp2_implicit, & ! Out
563 0 : term_wpsclrjp2_explicit ) ! Out
564 :
565 : ! <w'rt'sclr'> = coef_wprtpsclrp_implicit * <sclr'rt'>
566 : ! + term_wprtpsclrp_explicit
567 : call calc_coefs_wpxpyp_semiimpl( wp2(i,:), wprtp(i,:), wpsclrjp, & ! In
568 : mixt_frac(i,:), F_w(i,:), & ! In
569 : coef_sigma_rt_1_sqd, & ! In
570 : coef_sigma_rt_2_sqd, & ! In
571 : coef_sigma_sclrj_1_sqd, & ! In
572 : coef_sigma_sclrj_2_sqd, & ! In
573 : coef_wprtpsclrjp_implicit, & ! Out
574 0 : term_wprtpsclrjp_explicit ) ! Out
575 :
576 : ! <w'thl'sclr'> = coef_wpthlpsclrp_implicit * <sclr'thl'>
577 : ! + term_wpthlpsclrp_explicit
578 : call calc_coefs_wpxpyp_semiimpl( wp2(i,:), wpthlp(i,:), wpsclrjp, & ! In
579 : mixt_frac(i,:), F_w(i,:), & ! In
580 : coef_sigma_thl_1_sqd, & ! In
581 : coef_sigma_thl_2_sqd, & ! In
582 : coef_sigma_sclrj_1_sqd, & ! In
583 : coef_sigma_sclrj_2_sqd, & ! In
584 : coef_wpthlpsclrjp_implicit, & ! Out
585 0 : term_wpthlpsclrjp_explicit ) ! Out
586 :
587 0 : do k = 1, nz, 1
588 0 : coef_wpsclrp2_implicit(k,j) = coef_wpsclrjp2_implicit(k)
589 0 : term_wpsclrp2_explicit(k,j) = term_wpsclrjp2_explicit(k)
590 0 : coef_wprtpsclrp_implicit(k,j) = coef_wprtpsclrjp_implicit(k)
591 0 : term_wprtpsclrp_explicit(k,j) = term_wprtpsclrjp_explicit(k)
592 0 : coef_wpthlpsclrp_implicit(k,j) = coef_wpthlpsclrjp_implicit(k)
593 0 : term_wpthlpsclrp_explicit(k,j) = term_wpthlpsclrjp_explicit(k)
594 : enddo ! k = 1, nz, 1
595 :
596 : enddo ! j = 1, sclr_dim, 1
597 :
598 : endif ! sclr_dim > 0
599 :
600 : else ! l_explicit_turbulent_adv_xpyp
601 :
602 : ! Turbulent advection of <rt'^2>, <thl'^2>, <rt'thl'>, <u'^2>, <v'^2>,
603 : ! <sclr'^2>, <sclr'rt'>, and <sclr'thl'> are being handled explicitly.
604 : coef_wprtp2_implicit = zero
605 : term_wprtp2_explicit = zero
606 : coef_wpthlp2_implicit = zero
607 : term_wpthlp2_explicit = zero
608 : coef_wprtpthlp_implicit = zero
609 : term_wprtpthlp_explicit = zero
610 : coef_wpup2_implicit = zero
611 : term_wpup2_explicit = zero
612 : coef_wpvp2_implicit = zero
613 : term_wpvp2_explicit = zero
614 : if ( sclr_dim > 0 ) then
615 : coef_wpsclrp2_implicit = zero
616 : term_wpsclrp2_explicit = zero
617 : coef_wprtpsclrp_implicit = zero
618 : term_wprtpsclrp_explicit = zero
619 : coef_wpthlpsclrp_implicit = zero
620 : term_wpthlpsclrp_explicit = zero
621 : endif ! sclr_dim > 0
622 :
623 : endif ! .not. l_explicit_turbulent_adv_xpyp
624 :
625 : ! Set up a vector of 0s and an array of 0s to help write results back to
626 : ! type variable pdf_implicit_coefs_terms.
627 0 : zeros = zero
628 0 : zero_array = zero
629 :
630 : ! Pack the implicit coefficients and explicit terms into a single type
631 : ! variable for output.
632 0 : pdf_implicit_coefs_terms%coef_wp4_implicit(i,:) = coef_wp4_implicit
633 0 : pdf_implicit_coefs_terms%coef_wp2rtp_implicit(i,:) = coef_wp2rtp_implicit
634 0 : pdf_implicit_coefs_terms%term_wp2rtp_explicit(i,:) = zeros
635 0 : pdf_implicit_coefs_terms%coef_wp2thlp_implicit(i,:) = coef_wp2thlp_implicit
636 0 : pdf_implicit_coefs_terms%term_wp2thlp_explicit(i,:) = zeros
637 0 : pdf_implicit_coefs_terms%coef_wp2up_implicit(i,:) = coef_wp2up_implicit
638 0 : pdf_implicit_coefs_terms%term_wp2up_explicit(i,:) = zeros
639 0 : pdf_implicit_coefs_terms%coef_wp2vp_implicit(i,:) = coef_wp2vp_implicit
640 0 : pdf_implicit_coefs_terms%term_wp2vp_explicit(i,:) = zeros
641 0 : pdf_implicit_coefs_terms%coef_wprtp2_implicit(i,:) = coef_wprtp2_implicit
642 0 : pdf_implicit_coefs_terms%term_wprtp2_explicit(i,:) = term_wprtp2_explicit
643 0 : pdf_implicit_coefs_terms%coef_wpthlp2_implicit(i,:) = coef_wpthlp2_implicit
644 0 : pdf_implicit_coefs_terms%term_wpthlp2_explicit(i,:) = term_wpthlp2_explicit
645 0 : pdf_implicit_coefs_terms%coef_wprtpthlp_implicit(i,:) = coef_wprtpthlp_implicit
646 0 : pdf_implicit_coefs_terms%term_wprtpthlp_explicit(i,:) = term_wprtpthlp_explicit
647 0 : pdf_implicit_coefs_terms%coef_wpup2_implicit(i,:) = coef_wpup2_implicit
648 0 : pdf_implicit_coefs_terms%term_wpup2_explicit(i,:) = term_wpup2_explicit
649 0 : pdf_implicit_coefs_terms%coef_wpvp2_implicit(i,:) = coef_wpvp2_implicit
650 0 : pdf_implicit_coefs_terms%term_wpvp2_explicit(i,:) = term_wpvp2_explicit
651 0 : if ( sclr_dim > 0 ) then
652 0 : pdf_implicit_coefs_terms%coef_wp2sclrp_implicit(i,:,:) = coef_wp2sclrp_implicit
653 0 : pdf_implicit_coefs_terms%term_wp2sclrp_explicit(i,:,:) = zero_array
654 0 : pdf_implicit_coefs_terms%coef_wpsclrp2_implicit(i,:,:) = coef_wpsclrp2_implicit
655 0 : pdf_implicit_coefs_terms%term_wpsclrp2_explicit(i,:,:) = term_wpsclrp2_explicit
656 0 : pdf_implicit_coefs_terms%coef_wprtpsclrp_implicit(i,:,:) &
657 0 : = coef_wprtpsclrp_implicit
658 0 : pdf_implicit_coefs_terms%term_wprtpsclrp_explicit(i,:,:) &
659 0 : = term_wprtpsclrp_explicit
660 0 : pdf_implicit_coefs_terms%coef_wpthlpsclrp_implicit(i,:,:) &
661 0 : = coef_wpthlpsclrp_implicit
662 0 : pdf_implicit_coefs_terms%term_wpthlpsclrp_explicit(i,:,:) &
663 0 : = term_wpthlpsclrp_explicit
664 : endif ! sclr_dim > 0
665 :
666 : end do
667 :
668 0 : return
669 :
670 : end subroutine new_hybrid_pdf_driver
671 :
672 : !=============================================================================
673 0 : elemental subroutine calc_responder_driver( xm, xp2, wpxp, wp2, & ! In
674 : mixt_frac, F_w, & ! In
675 : Skx, & ! In/Out
676 : mu_x_1, mu_x_2, & ! Out
677 : sigma_x_1_sqd, & ! Out
678 : sigma_x_2_sqd, & ! Out
679 : coef_sigma_x_1_sqd, & ! Out
680 : coef_sigma_x_2_sqd ) ! Out
681 :
682 : ! Description:
683 : ! This is the sub-driver for a responder variable. The limits of the range
684 : ! of values of Skx that the PDF is able to represent are calculated, and
685 : ! Skx is clipped when it exceeds the upper or lower limit. Then, the PDF
686 : ! parameters for responder variable x are calculated.
687 :
688 : ! References:
689 : !-----------------------------------------------------------------------
690 :
691 : use constants_clubb, only: &
692 : three, & ! Constant(s)
693 : two, &
694 : three_halves, &
695 : one, &
696 : zero
697 :
698 : use new_hybrid_pdf, only: &
699 : calculate_responder_params ! Procedure(s)
700 :
701 : use clubb_precision, only: &
702 : core_rknd ! Variable(s)
703 :
704 : implicit none
705 :
706 : ! Input Variables
707 : real( kind = core_rknd ), intent(in) :: &
708 : xm, & ! Mean of x (overall) [units vary]
709 : xp2, & ! Variance of x (overall) [(units vary)^2]
710 : wpxp, & ! Covariance of w and x (overall) [m/s (units vary)]
711 : wp2, & ! Variance of w (overall) [m^2/s^2]
712 : mixt_frac, & ! Mixture fraction [-]
713 : F_w ! Parameter for the spread of the PDF comp. means of w [-]
714 :
715 : ! Input/Output Variable
716 : real( kind = core_rknd ), intent(inout) :: &
717 : Skx ! Skewness of x (overall) [-]
718 :
719 : ! Output Variables
720 : real( kind = core_rknd ), intent(out) :: &
721 : mu_x_1, & ! Mean of x (1st PDF component) [units vary]
722 : mu_x_2, & ! Mean of x (2nd PDF component) [units vary]
723 : sigma_x_1_sqd, & ! Variance of x (1st PDF component) [(units vary)^2]
724 : sigma_x_2_sqd ! Variance of x (2nd PDF component) [(units vary)^2]
725 :
726 : real( kind = core_rknd ), intent(out) :: &
727 : coef_sigma_x_1_sqd, & ! sigma_x_1^2 = coef_sigma_x_1_sqd * <x'^2> [-]
728 : coef_sigma_x_2_sqd ! sigma_x_2^2 = coef_sigma_x_2_sqd * <x'^2> [-]
729 :
730 : ! Local Variables
731 : real( kind = core_rknd ) :: &
732 : corr_w_x, & ! Correlation (overall) of w and x [-]
733 : min_Skx, & ! Minimum Skx that can be represented by the PDF [-]
734 : max_Skx ! Maximum Skx that can be represented by the PDF [-]
735 :
736 :
737 : ! Calculate the limits of Skx.
738 0 : if ( F_w > zero ) then
739 :
740 : ! Calculate the overall correlation of w and x.
741 0 : if ( wp2 * xp2 > zero ) then
742 :
743 : ! The overall correlation of w and x is:
744 : !
745 : ! corr_w_x = <w'x'> / sqrt( <w'^2> * <x'^2> ).
746 0 : corr_w_x = wpxp / sqrt( wp2 * xp2 )
747 :
748 : else ! <w'^2> * <x'^2> = 0
749 :
750 : ! When <w'^2> = 0 or <x'^2> = 0, <w'x'> = 0. The correlation of w
751 : ! and x is undefined, however, since <w'x'> = 0, Skx = 0. Setting
752 : ! corr_w_x = 0 in this scenario will set max_Skx = min_Skx = 0.
753 : corr_w_x = zero
754 :
755 : endif
756 :
757 0 : if ( wpxp >= zero ) then
758 :
759 : min_Skx = ( one + mixt_frac ) &
760 : / sqrt( mixt_frac * ( one - mixt_frac ) ) &
761 : * corr_w_x**3 / F_w**three_halves &
762 : - sqrt( mixt_frac / ( one - mixt_frac ) ) &
763 0 : * three * corr_w_x / sqrt( F_w )
764 :
765 : max_Skx = ( mixt_frac - two ) &
766 : / sqrt( mixt_frac * ( one - mixt_frac ) ) &
767 : * corr_w_x**3 / F_w**three_halves &
768 : + sqrt( ( one - mixt_frac ) / mixt_frac ) &
769 0 : * three * corr_w_x / sqrt( F_w )
770 :
771 : else ! <w'x'> < 0
772 :
773 : min_Skx = ( mixt_frac - two ) &
774 : / sqrt( mixt_frac * ( one - mixt_frac ) ) &
775 : * corr_w_x**3 / F_w**three_halves &
776 : + sqrt( ( one - mixt_frac ) / mixt_frac ) &
777 0 : * three * corr_w_x / sqrt( F_w )
778 :
779 : max_Skx = ( one + mixt_frac ) &
780 : / sqrt( mixt_frac * ( one - mixt_frac ) ) &
781 : * corr_w_x**3 / F_w**three_halves &
782 : - sqrt( mixt_frac / ( one - mixt_frac ) ) &
783 0 : * three * corr_w_x / sqrt( F_w )
784 :
785 : endif ! <w'x'> >= 0
786 :
787 : else ! F_w > 0
788 :
789 : ! When F_w = 0, <w'x'> = 0, and Skx = 0.
790 : min_Skx = zero
791 : max_Skx = zero
792 :
793 : endif ! F_w > 0
794 :
795 : ! Limit Skx so that min_Skx <= Skx <= max_Skx.
796 0 : if ( Skx > max_Skx ) then
797 0 : Skx = max_Skx
798 0 : elseif ( Skx < min_Skx ) then
799 0 : Skx = min_Skx
800 : endif ! Skx >= max_Skx
801 :
802 : ! Calculate the PDF parameters for responder variable x.
803 : call calculate_responder_params( xm, xp2, Skx, wpxp, & ! In
804 : wp2, F_w, mixt_frac, & ! In
805 : mu_x_1, mu_x_2, & ! Out
806 : sigma_x_1_sqd, & ! Out
807 : sigma_x_2_sqd, & ! Out
808 : coef_sigma_x_1_sqd, & ! Out
809 0 : coef_sigma_x_2_sqd ) ! Out
810 :
811 :
812 0 : return
813 :
814 : end subroutine calc_responder_driver
815 :
816 : !=============================================================================
817 0 : elemental subroutine calc_F_w_zeta_w( Skw, wprtp, wpthlp, upwp, vpwp, & ! In
818 : wp2, rtp2, thlp2, up2, vp2, & ! In
819 : gamma_Skw_fnc, & ! In
820 : slope_coef_spread_DG_means_w, & ! In
821 : pdf_component_stdev_factor_w, & ! In
822 : max_corr_w_sclr_sqd, & ! In
823 : F_w, zeta_w, min_F_w, max_F_w ) ! Out
824 :
825 : ! Description:
826 : ! Calculates the values of F_w and zeta_w for w, which is the setter
827 : ! variable (which is the variable that sets the mixture fraction).
828 : !
829 : ! Based purely on the PDF of w and not considering restrictions caused by
830 : ! other variables in the multivariate PDF, the value of F_w is calculated
831 : ! between 0 (min_F_w) and 1 (max_F_w). However, other variables in the PDF
832 : ! place restrictions on the minimum value of F_w. The range of acceptible
833 : ! values for F_w is given by:
834 : !
835 : ! max( <w'x'>^2 / ( <w'^2> * <x'^2> ), for all variables x ) <= F_w <= 1.
836 : !
837 : ! Additionally, the value of F_w should still increase with an increasing
838 : ! magnitude of Skw. The value of F_w can be 0 only when Skw = 0 (as well as
839 : ! all <w'x'> = 0, as well). The F_w skewness equation used is:
840 : !
841 : ! F_w = max_F_w + ( min_F_w - max_F_w )
842 : ! * exp{ -|Skw|^lambda / slope_coef_spread_DG_means_w };
843 : !
844 : ! where lambda > 0 and slope_coef_spread_DG_means_w > 0. As |Skw| goes
845 : ! toward 0, the value of F_w goes toward min_F_w, and as |Skw| becomes
846 : ! large, the value of F_w goes toward max_F_w. When the value of
847 : ! slope_coef_spread_DG_means_w is small, the value of F_w tends toward
848 : ! max_F_w, and when the value of slope_coef_spread_DG_means_w is large, the
849 : ! value of F_w tends toward min_F_w. When lambda is small, the value of F_w
850 : ! is less dependent on Skw, and when lambda is large, the value of F_w is
851 : ! more dependent on Skw.
852 : !
853 : ! Mathematically, this equation will always produce a value of F_w that
854 : ! falls between min_F_w and max_F_w. However, in order to prevent a value
855 : ! of F_w from being calculated outside the bounds of min_F_w and max_F_w
856 : ! owing to numerical underflow or loss of precision, this equation can be
857 : ! rewritten as:
858 : !
859 : ! F_w
860 : ! = min_F_w * exp{ -|Skw|^lambda / slope_coef_spread_DG_means_w }
861 : ! + max_F_w * ( 1 - exp{ -|Skw|^lambda / slope_coef_spread_DG_means_w } ).
862 : !
863 : ! The value of zeta_w used to adjust the PDF component standard devations:
864 : !
865 : ! 1 + zeta_w = ( mixt_frac * sigma_w_1^2 )
866 : ! / ( ( 1 - mixt_frac ) * sigma_w_2^2 );
867 : !
868 : ! where zeta_w > -1. The sign of zeta_w is used to easily determine if
869 : ! mixt_frac * sigma_w_1^2 is greater than ( 1 - mixt_frac ) * sigma_w_2^2
870 : ! (when zeta_w is positive), mixt_frac * sigma_w_1^2 is less than
871 : ! ( 1 - mixt_frac ) * sigma_w_2^2 (when zeta_w is negative), or
872 : ! mixt_frac * sigma_w_1^2 is equal to ( 1 - mixt_frac ) * sigma_w_2^2 (when
873 : ! zeta_w is 0).
874 : !
875 : ! In order to allow for a tunable parameter that is the pure ratio of
876 : ! mixt_frac * sigma_w_1^2 to ( 1 - mixt_frac ) * sigma_w_2^2, zeta_w is
877 : ! related to the parameter pdf_component_stdev_factor_w, where:
878 : !
879 : ! 1 + zeta_w = pdf_component_stdev_factor_w.
880 :
881 : ! References:
882 : !-----------------------------------------------------------------------
883 :
884 : use constants_clubb, only: &
885 : one, & ! Variable(s)
886 : zero, &
887 : max_mag_correlation_flux
888 :
889 : use clubb_precision, only: &
890 : core_rknd ! Variable(s)
891 :
892 : implicit none
893 :
894 : ! Input Variables
895 : real( kind = core_rknd ), intent(in) :: &
896 : Skw, & ! Skewness of w (overall) [-]
897 : wprtp, & ! Covariance (overall) of w and rt [m/s (kg/kg)]
898 : wpthlp, & ! Covariance (overall) of w and thl [m/s K]
899 : upwp, & ! Covariance (overall) of u and w [m^2/s^2]
900 : vpwp, & ! Covariance (overall) of u and v [m^2/s^2]
901 : wp2, & ! Variance (overall) of w [m^2/s^2]
902 : rtp2, & ! Variance (overall) of rt [(kg/kg)^2]
903 : thlp2, & ! Variance (overall) of thl [K^2]
904 : up2, & ! Variance (overall) of u [m^2/s^2]
905 : vp2 ! Variance (overall) of v [m^2/s^2]
906 :
907 : ! Tunable parameter gamma.
908 : ! When gamma goes to 0, the standard deviations of w in each PDF component
909 : ! become small, and the spread between the two PDF component means of w
910 : ! becomes large. F_w goes to min_F_w.
911 : ! When gamma goes to 1, the standard deviations of w in each PDF component
912 : ! become large, and the spread between the two PDF component means of w
913 : ! becomes small. F_w goes to max_F_w.
914 : real( kind = core_rknd ), intent(in) :: &
915 : gamma_Skw_fnc ! Value of parameter gamma from tunable Skw function [-]
916 :
917 : real( kind = core_rknd ), intent(in) :: &
918 : ! Slope coefficient for the spread between the PDF component means of w.
919 : slope_coef_spread_DG_means_w, &
920 : ! Parameter to adjust the PDF component standard deviations of w.
921 : pdf_component_stdev_factor_w
922 :
923 : real( kind = core_rknd ), intent(in) :: &
924 : max_corr_w_sclr_sqd ! Max value of wpsclrp^2 / ( wp2 * sclrp2 ) [-]
925 :
926 : ! Output Variables
927 : real( kind = core_rknd ), intent(out) :: &
928 : F_w, & ! Parameter for the spread of the PDF component means of w [-]
929 : zeta_w, & ! Parameter for the PDF component variances of w [-]
930 : min_F_w, & ! Minimum allowable value of parameter F_w [-]
931 : max_F_w ! Maximum allowable value of parameter F_w [-]
932 :
933 : ! Local Variables
934 : real( kind = core_rknd ) :: &
935 : corr_w_rt_sqd, & ! Value of wprtp^2 / ( wp2 * rtp2 ) [-]
936 : corr_w_thl_sqd, & ! Value of wpthlp^2 / ( wp2 * thlp2 ) [-]
937 : corr_u_w_sqd, & ! Value of upwp^2 / ( up2 * wp2 ) [-]
938 : corr_v_w_sqd, & ! Value of vpwp^2 / ( vp2 * wp2 ) [-]
939 : max_corr_w_x_sqd ! Max. val. of wpxp^2/(wp2*xp2) for all vars. x [-]
940 :
941 : real( kind = core_rknd ) :: &
942 : exp_Skw_interp_factor, & ! Function to interp. between min. and max. [-]
943 : zeta_w_star
944 :
945 : real( kind = core_rknd ), parameter :: &
946 : lambda = 0.5_core_rknd ! Parameter for Skw dependence [-]
947 :
948 :
949 : ! Find the maximum value of <w'x'>^2 / ( <w'^2> * <x'^2> ) for all
950 : ! variables x that are Double Gaussian PDF responder variables. This
951 : ! includes rt and theta-l. When l_predict_upwp_vpwp is enabled, u and v are
952 : ! also calculated as part of the PDF, and they are included as well.
953 : ! Additionally, when sclr_dim > 0, passive scalars (sclr) are also included.
954 :
955 : ! Calculate the squared value of the correlation of w and rt.
956 0 : if ( wp2 * rtp2 > zero ) then
957 0 : corr_w_rt_sqd = wprtp**2 / ( wp2 * rtp2 )
958 : else
959 : corr_w_rt_sqd = zero
960 : endif
961 :
962 : ! Calculate the squared value of the correlation of w and thl.
963 0 : if ( wp2 * thlp2 > zero ) then
964 0 : corr_w_thl_sqd = wpthlp**2 / ( wp2 * thlp2 )
965 : else
966 : corr_w_thl_sqd = zero
967 : endif
968 :
969 : ! Calculate the squared value of the correlation of u and w.
970 0 : if ( up2 * wp2 > zero ) then
971 0 : corr_u_w_sqd = upwp**2 / ( up2 * wp2 )
972 : else
973 : corr_u_w_sqd = zero
974 : endif
975 :
976 : ! Calculate the squared value of the correlation of v and w.
977 0 : if ( vp2 * wp2 > zero ) then
978 0 : corr_v_w_sqd = vpwp**2 / ( vp2 * wp2 )
979 : else
980 : corr_v_w_sqd = zero
981 : endif
982 :
983 : ! Calculate the maximum value of corr_w_rt^2, corr_w_thl^2, corr_u_w^2, and
984 : ! corr_v_w^2 in the calculation of max(corr_w_x^2). Include the maximum
985 : ! value of corr_w_sclr^2 (for all sclrs) when scalars are found.
986 : max_corr_w_x_sqd = max( corr_w_rt_sqd, corr_w_thl_sqd, &
987 0 : corr_u_w_sqd, corr_v_w_sqd, max_corr_w_sclr_sqd )
988 :
989 0 : max_corr_w_x_sqd = min( max_corr_w_x_sqd, max_mag_correlation_flux**2 )
990 :
991 : ! Calculate min_F_w and set max_F_w to 1.
992 0 : if ( abs( Skw ) > zero ) then
993 0 : min_F_w = max( max_corr_w_x_sqd, 1.0e-3_core_rknd )
994 : else
995 0 : min_F_w = max_corr_w_x_sqd
996 : endif
997 0 : max_F_w = one
998 :
999 : ! F_w must have a value between min_F_w and max_F_w.
1000 : exp_Skw_interp_factor &
1001 : = exp( -abs(Skw)**lambda / slope_coef_spread_DG_means_w )
1002 :
1003 : ! F_w = min_F_w * exp_Skw_interp_factor &
1004 : ! + max_F_w * ( one - exp_Skw_interp_factor )
1005 :
1006 : ! For now, use a formulation similar to what is used for ADG1.
1007 : ! Tunable parameter gamma.
1008 : ! When gamma goes to 0, the standard deviations of w in each PDF component
1009 : ! become small, and the spread between the two PDF component means of w
1010 : ! becomes large. F_w goes to min_F_w.
1011 : ! When gamma goes to 1, the standard deviations of w in each PDF component
1012 : ! become large, and the spread between the two PDF component means of w
1013 : ! becomes small. F_w goes to max_F_w.
1014 0 : F_w = max_F_w - gamma_Skw_fnc * ( max_F_w - min_F_w )
1015 :
1016 : ! The value of zeta_w must be greater than -1.
1017 0 : zeta_w_star = pdf_component_stdev_factor_w - one
1018 :
1019 : ! Make the PDF of w symmetric. In other words, the PDF at a value of
1020 : ! positive skewness will look like a mirror image of the PDF at the
1021 : ! opposite value of negative skewness.
1022 0 : if ( Skw >= zero ) then
1023 0 : zeta_w = zeta_w_star
1024 : else ! Skw < 0
1025 0 : zeta_w = - zeta_w_star / ( zeta_w_star + one )
1026 : endif
1027 :
1028 :
1029 0 : return
1030 :
1031 : end subroutine calc_F_w_zeta_w
1032 :
1033 : !=============================================================================
1034 :
1035 : end module new_hybrid_pdf_main
|