Line data Source code
1 : ! $Id$
2 : !===============================================================================
3 : module new_hybrid_pdf
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 : ! Griffin and Larson (2020)
14 : !-------------------------------------------------------------------------
15 :
16 : implicit none
17 :
18 : public :: calculate_mixture_fraction, & ! Procedure(s)
19 : calculate_w_params, &
20 : calculate_responder_params, &
21 : calculate_coef_wp4_implicit, &
22 : calc_coef_wp2xp_implicit, &
23 : calc_coefs_wpxp2_semiimpl, &
24 : calc_coefs_wpxpyp_semiimpl
25 :
26 : private
27 :
28 : contains
29 :
30 : !=============================================================================
31 : !
32 : ! DESCRIPTION OF THE METHOD FOR THE VARIABLE THAT SETS THE MIXTURE FRACTION
33 : ! =========================================================================
34 : !
35 : ! The variable that sets the mixture fraction for the PDF is w. There are
36 : ! five PDF parameters that need to be calculated, which are mu_w_1 (the mean
37 : ! of w is the 1st PDF component), mu_w_2 (the mean of w in the 2nd PDF
38 : ! component), sigma_w_1 (the standard deviation of w in the 1st PDF
39 : ! component), sigma_w_2 (the standard deviation of w in the 2nd PDF
40 : ! component), and mixt_frac (the mixture fraction, which is the weight of the
41 : ! 1st PDF component). In order to solve for these five parameters, five
42 : ! equations are needed. These five equations are the equations for <w>,
43 : ! <w'^2>, and <w'^3> as found by integrating over the PDF. Additionally, two
44 : ! more equations, which involve tunable parameters F_w and zeta_w, and which
45 : ! are used to help control the spread of the PDF component means and the size
46 : ! of the PDF component standard deviations compared to each other,
47 : ! respectively, are used in this equation set. The five equations are:
48 : !
49 : ! <w> = mixt_frac * mu_w_1 + ( 1 - mixt_frac ) * mu_w_2;
50 : !
51 : ! <w'^2> = mixt_frac * ( ( mu_w_1 - <w> )^2 + sigma_w_1^2 )
52 : ! + ( 1 - mixt_frac ) * ( ( mu_w_2 - <w> )^2 + sigma_w_2^2 );
53 : !
54 : ! <w'^3> = mixt_frac * ( mu_w_1 - <w> )
55 : ! * ( ( mu_w_1 - <w> )^2 + 3 * sigma_w_1^2 )
56 : ! + ( 1 - mixt_frac ) * ( mu_w_2 - <w> )
57 : ! * ( ( mu_w_2 - <w> )^2 + 3 * sigma_w_2^2 );
58 : !
59 : ! mu_w_1 - <w> = sqrt(F_w) * ( sqrt( 1 - mixt_frac ) / sqrt( mixt_frac ) )
60 : ! * sqrt( <w'^2> );
61 : !
62 : ! where 0 <= F_w <= 1; and
63 : !
64 : ! 1 + zeta_w = ( mixt_frac * sigma_w_1^2 )
65 : ! / ( ( 1 - mixt_frac ) * sigma_w_2^2 );
66 : !
67 : ! where zeta_w > -1.
68 : !
69 : ! Following convention for w, mu_w_1 is defined to be greater than or equal to
70 : ! mu_w_2 (and is also greater than or equal to <w>, while mu_w_2 is less than
71 : ! or equal to <w>). This is relationship is found in the mu_w_1 - <w>
72 : ! equation above.
73 : !
74 : ! The resulting equations for the five PDF parameters are:
75 : !
76 : ! mixt_frac
77 : ! = ( 4 * F_w^3
78 : ! + 18 * F_w * ( zeta_w + 1 ) * ( 1 - F_w ) / ( zeta_w + 2 )
79 : ! + 6 * F_w^2 * ( 1 - F_w ) / ( zeta_w + 2 )
80 : ! + Skw^2
81 : ! - Skw * sqrt( 4 * F_w^3
82 : ! + 12 * F_w^2 * ( 1 - F_w )
83 : ! + 36 * F_w * ( zeta_w + 1 ) * ( 1 - F_w )^2
84 : ! / ( zeta_w + 2 )^2
85 : ! + Skw^2 ) )
86 : ! / ( 2 * F_w * ( F_w - 3 )^2 + 2 * Skw^2 );
87 : !
88 : ! mu_w_1 = <w> + sqrt( F_w * ( ( 1 - mixt_frac ) / mixt_frac ) * <w'^2> );
89 : !
90 : ! mu_w_2 = <w> - ( mixt_frac / ( 1 - mixt_frac ) ) * ( mu_w_1 - <w> );
91 : !
92 : ! sigma_w_1 = sqrt( ( ( zeta_w + 1 ) * ( 1 - F_w ) )
93 : ! / ( ( zeta_w + 2 ) * mixt_frac ) * <w'^2> ); and
94 : !
95 : ! sigma_w_2 = sqrt( ( mixt_frac * sigma_w_1^2 )
96 : ! / ( ( 1 - mixt_frac ) * ( 1 + zeta_w ) ) );
97 : !
98 : ! where Skw is the skewness of w, and Skw = <w'^3> / <w'^2>^(3/2).
99 : !
100 : ! This method works for all values of F_w (where 0 <= F_w <= 1) and zeta_w
101 : ! (where zeta_w > -1).
102 : !
103 : !
104 : ! Special case:
105 : !
106 : ! When Skw = 0 and F_w = 0, the equation for mixt_frac is undefined. The
107 : ! equation for mixture fraction in this scenario can be derived by using the
108 : ! above equation for mixture fraction and then setting Skw = 0. The resulting
109 : ! equation becomes:
110 : !
111 : ! mixt_frac
112 : ! = ( 4 * F_w^3
113 : ! + 18 * F_w * ( zeta_w + 1 ) * ( 1 - F_w ) / ( zeta_w + 2 )
114 : ! + 6 * F_w^2 * ( 1 - F_w ) / ( zeta_w + 2 ) )
115 : ! / ( 2 * F_w * ( F_w - 3 )^2 ).
116 : !
117 : ! All of the terms in the numerator and denominator contain a F_w, so this
118 : ! equation can be rewritten as:
119 : !
120 : ! mixt_frac
121 : ! = ( 4 * F_w^2
122 : ! + 18 * ( zeta_w + 1 ) * ( 1 - F_w ) / ( zeta_w + 2 )
123 : ! + 6 * F_w * ( 1 - F_w ) / ( zeta_w + 2 ) )
124 : ! / ( 2 * ( F_w - 3 )^2 ).
125 : !
126 : ! Now setting F_w = 0, the equation becomes:
127 : !
128 : ! mixt_frac = ( 18 * ( zeta_w + 1 ) / ( zeta_w + 2 ) ) / 18;
129 : !
130 : ! which can be rewritten as:
131 : !
132 : ! mixt_frac = ( zeta_w + 1 ) / ( zeta_w + 2 ).
133 : !
134 : ! When F_w = 0, Skw must have a value of 0 in order for the PDF to function
135 : ! correctly. When F_w = 0, mu_w_1 = mu_w_2. When the two PDF component means
136 : ! are equal to each other (and to the overall mean, <w>), the only value of
137 : ! Skw that can be represented is a value of 0. In the equation for mixture
138 : ! fraction, when F_w is set to 0, but | Skw | > 0, mixt_frac will either have
139 : ! a value of 0 or 1, depending on whether Skw is positive or negative,
140 : ! respectively.
141 : !
142 : ! The value of F_w should be set as a function of Skw. The value F_w should
143 : ! go toward 0 as | Skw | (or Skw^2) goes toward 0. The value of F_w should
144 : ! go toward 1 as | Skw | (or Skw^2) goes to infinity.
145 : !
146 : !
147 : ! Tunable parameters:
148 : !
149 : ! 1) F_w: This parameter controls the spread of the PDF component means. The
150 : ! range of this parameter is 0 <= F_w <= 1. When F_w = 0, the two
151 : ! PDF component means (mu_w_1 and mu_w_2) are equal to each other
152 : ! (and Skw must equal 0). All of the variance of w is accounted for
153 : ! by the PDF component standard deviations (sigma_w_1 and sigma_w_2).
154 : ! When F_w = 1, mu_w_1 and mu_w_2 are spread as far apart as they can
155 : ! be. Both PDF component standard deviations (sigma_w_1 and
156 : ! sigma_w_2) are equal to 0, and all of the variance of w is
157 : ! accounted for by the spread of the PDF component means.
158 : !
159 : ! When sigma_w_1 = sigma_w_2 = 0, the equation for <w'^2> becomes:
160 : !
161 : ! <w'^2> = mixt_frac * ( mu_w_1 - <w> )^2
162 : ! + ( 1 - mixt_frac ) * ( mu_w_2 - <w> )^2.
163 : !
164 : ! Substituting the equation for <w> into the above equation for
165 : ! mu_w_2 - <w>, the above equation becomes:
166 : !
167 : ! <w'^2> = ( mixt_frac / ( 1 - mixt_frac ) ) * ( mu_w_1 - <w> )^2;
168 : !
169 : ! which can be rewritten as:
170 : !
171 : ! ( mu_w_1 - <w> )^2 = ( ( 1 - mixt_frac ) / mixt_frac ) * <w'^2>.
172 : !
173 : ! Taking the square root of the above equation:
174 : !
175 : ! mu_w_1 - <w> = +/- ( sqrt( 1 - mixt_frac ) / sqrt(mixt_frac) )
176 : ! * sqrt( <w'^2> ).
177 : !
178 : ! This equation can be compared to the equation for mu_w_1 - <w> in
179 : ! the set of 5 equations, which is:
180 : !
181 : ! mu_w_1 - <w>
182 : ! = sqrt(F_w) * ( sqrt( 1 - mixt_frac ) / sqrt( mixt_frac ) )
183 : ! * sqrt( <w'^2> ).
184 : !
185 : ! The above equations give another example of the meaning of F_w.
186 : ! The value of sqrt(F_w) is ratio of mu_w_1 - <w> to its maximum
187 : ! value, which is:
188 : !
189 : ! sqrt( ( ( 1 - mixt_frac ) / mixt_frac ) * <w'^2> ).
190 : !
191 : !
192 : ! 2) zeta_w: This parameter controls the size of the PDF component standard
193 : ! deviations compared to each other. The equation for zeta_w is:
194 : !
195 : ! 1 + zeta_w = ( mixt_frac * sigma_w_1^2 )
196 : ! / ( ( 1 - mixt_frac ) * sigma_w_2^2 ).
197 : !
198 : ! When zeta_w > 0, mixt_frac * sigma_w_1^2 increases at the
199 : ! expense of ( 1 - mixt_frac ) * sigma_w_2^2, which decreases in
200 : ! this variance-preserving equation set. When zeta_w = 0, then
201 : ! mixt_frac * sigma_w_1^2 = ( 1 - mixt_frac ) * sigma_w_2^2.
202 : ! When -1 < zeta_w < 0, ( 1 - mixt_frac ) * sigma_w_2^2 increases
203 : ! at the expense of mixt_frac * sigma_w_1^2, which decreases. As
204 : ! a result, greater values of zeta_w cause the 1st PDF component
205 : ! to become broader while the 2nd PDF component becomes narrower,
206 : ! and smaller values of zeta_w cause the 1st PDF component to
207 : ! become narrower while the 2nd PDF component becomes broader.
208 : !
209 : ! Symmetry
210 : !
211 : ! When zeta_w = 0, the PDF is always symmetric. In other words,
212 : ! the PDF at any positive value of Skw (for example, Skw = 2.5)
213 : ! will look like a mirror-image (reflection across the y-axis)
214 : ! of the PDF at a negative value of Skw of the same magnitude (in
215 : ! this example, Skw = -2.5). However, when zeta_w /= 0, the PDF
216 : ! loses this quality and is not symmetric.
217 : !
218 : ! When symmetry is desired at values of zeta_w besides zeta_w = 0,
219 : ! the solution is to turn zeta_w into a function of Skw. A basic
220 : ! example of a zeta_w skewness equation that produces a symmetric
221 : ! PDF for values of zeta_w other than 0 is:
222 : !
223 : ! zeta_w = | zeta_w_in, when Skw >= 0;
224 : ! | ( 1 / ( 1 + zeta_w_in ) ) - 1, when Skw < 0.
225 : !
226 : !
227 : ! Notes:
228 : !
229 : ! When F_w = 0 (which can only happen when Skw = 0), mu_w_1 = mu_w_2, and
230 : ! mixt_frac = ( zeta_w + 1 ) / ( zeta_w + 2 ). When these equations are
231 : ! substituted into the equations for sigma_w_1 and sigma_w_2, the result is
232 : ! sigma_w_1 = sigma_w_2 = sqrt( <w'^2> ). This means that the distribution
233 : ! becomes a single Gaussian when F_w = 0 (and Skw = 0). This happens
234 : ! regardless of the value of zeta_w.
235 : !
236 : ! The equations for the PDF component means and standard deviations can also
237 : ! be written as:
238 : !
239 : ! mu_w_1 = <w> + sqrt( F_w * ( ( 1 - mixt_frac ) / mixt_frac ) * <w'^2> );
240 : !
241 : ! mu_w_2 = <w> - sqrt( F_w * ( mixt_frac / ( 1 - mixt_frac ) ) * <w'^2> );
242 : !
243 : ! sigma_w_1 = sqrt( coef_sigma_w_1_sqd * <w'^2> ); and
244 : !
245 : ! sigma_w_2 = sqrt( coef_sigma_w_2_sqd * <w'^2> ); where
246 : !
247 : ! coef_sigma_w_1_sqd = ( ( zeta_w + 1 ) * ( 1 - F_w ) )
248 : ! / ( ( zeta_w + 2 ) * mixt_frac ); and
249 : !
250 : ! coef_sigma_w_2_sqd = ( 1 - F_w ) / ( ( zeta_w + 2 ) * ( 1 - mixt_frac ) ).
251 : !
252 : ! The above equations can be substituted into an equation for a variable that
253 : ! has been derived by integrating over the PDF. Many variables like this are
254 : ! used in parts of the predictive equation set. These substitutions allow
255 : ! some terms to solved implicitly or semi-implicitly in the predictive
256 : ! equations.
257 : !
258 : !
259 : ! Brian Griffin; September 2017.
260 : !
261 : !=============================================================================
262 0 : elemental function calculate_mixture_fraction( Skw, F_w, zeta_w ) &
263 : result( mixt_frac )
264 :
265 : ! Description:
266 : ! Calculates mixture fraction.
267 :
268 : ! References:
269 : ! Griffin and Larson (2020)
270 : !-----------------------------------------------------------------------
271 :
272 : use constants_clubb, only: &
273 : thirty_six, & ! Constant(s)
274 : eighteen, &
275 : twelve, &
276 : six, &
277 : four, &
278 : three, &
279 : two, &
280 : one, &
281 : zero
282 :
283 : use clubb_precision, only: &
284 : core_rknd ! Variable(s)
285 :
286 : implicit none
287 :
288 : ! Input Variables
289 : real( kind = core_rknd ), intent(in) :: &
290 : Skw, & ! Skewness of w [-]
291 : F_w, & ! Parameter for the spread of the PDF component means of w [-]
292 : zeta_w ! Parameter for the PDF component variances of w [-]
293 :
294 : ! Return Variable
295 : real( kind = core_rknd ) :: &
296 : mixt_frac ! Mixture fraction [-]
297 :
298 :
299 : ! Calculate mixture fraction, which is the weight of the 1st PDF component.
300 : ! The 2nd PDF component has a weight of 1 - mixt_frac.
301 0 : if ( F_w > zero ) then
302 :
303 : mixt_frac &
304 : = ( four * F_w**3 &
305 : + eighteen * F_w &
306 : * ( zeta_w + one ) * ( one - F_w ) / ( zeta_w + two ) &
307 : + six * F_w**2 * ( one - F_w ) / ( zeta_w + two ) &
308 : + Skw**2 &
309 : - Skw * sqrt( four * F_w**3 &
310 : + twelve * F_w**2 * ( one - F_w ) &
311 : + thirty_six * F_w &
312 : * ( zeta_w + one ) * ( one - F_w )**2 &
313 : / ( zeta_w + two )**2 &
314 : + Skw**2 ) ) &
315 0 : / ( two * F_w * ( F_w - three )**2 + two * Skw**2 )
316 :
317 : else ! F_w = 0
318 :
319 0 : if ( abs( Skw ) > zero ) then
320 :
321 : ! When F_w = 0, | Skw | must have a value of 0. In a scenario where
322 : ! F_w = 0 and | Skw | > 0, the mixture fraction (and the rest of the
323 : ! PDF parameters) can't be calculated. Since mixture fraction can
324 : ! only have values 0 < mixt_frac < 1, set mixt_frac to -1 in this
325 : ! scenario.
326 : mixt_frac = -one
327 :
328 : else ! Skw = 0
329 :
330 0 : mixt_frac = ( zeta_w + one ) / ( zeta_w + two )
331 :
332 : endif ! | Skw | > 0
333 :
334 : endif ! F_w > 0
335 :
336 :
337 : return
338 :
339 : end function calculate_mixture_fraction
340 :
341 : !=============================================================================
342 0 : elemental subroutine calculate_w_params( wm, wp2, Skw, F_w, zeta_w, & ! In
343 : mu_w_1, mu_w_2, sigma_w_1, & ! Out
344 : sigma_w_2, mixt_frac, & ! Out
345 : coef_sigma_w_1_sqd, & ! Out
346 : coef_sigma_w_2_sqd ) ! Out
347 :
348 : ! Description:
349 : ! Calculates the PDF component means, the PDF component standard deviations,
350 : ! and the mixture fraction for the variable that sets the PDF.
351 :
352 : ! References:
353 : ! Griffin and Larson (2020)
354 : !-----------------------------------------------------------------------
355 :
356 : use constants_clubb, only: &
357 : two, & ! Variable(s)
358 : one, &
359 : zero
360 :
361 : use clubb_precision, only: &
362 : core_rknd ! Variable(s)
363 :
364 : implicit none
365 :
366 : ! Input Variables
367 : real( kind = core_rknd ), intent(in) :: &
368 : wm, & ! Mean of w (overall) [m/s]
369 : wp2, & ! Variance of w (overall) [m^2/s^2]
370 : Skw, & ! Skewness of w [-]
371 : F_w, & ! Parameter for the spread of the PDF component means of w [-]
372 : zeta_w ! Parameter for the PDF component variances of w [-]
373 :
374 : ! Output Variables
375 : real( kind = core_rknd ), intent(out) :: &
376 : mu_w_1, & ! Mean of w (1st PDF component) [m/s]
377 : mu_w_2, & ! Mean of w (2nd PDF component) [m/s]
378 : sigma_w_1, & ! Standard deviation of w (1st PDF component) [m/s]
379 : sigma_w_2, & ! Standard deviation of w (2nd PDF component) [m/s]
380 : mixt_frac ! Mixture fraction [-]
381 :
382 : real( kind = core_rknd ), intent(out) :: &
383 : coef_sigma_w_1_sqd, & ! sigma_w_1^2 = coef_sigma_w_1_sqd * <w'^2> [-]
384 : coef_sigma_w_2_sqd ! sigma_w_2^2 = coef_sigma_w_2_sqd * <w'^2> [-]
385 :
386 :
387 : ! Calculate the mixture fraction.
388 0 : mixt_frac = calculate_mixture_fraction( Skw, F_w, zeta_w )
389 :
390 0 : if ( mixt_frac > zero .and. mixt_frac < one ) then
391 :
392 : ! Calculate the mean of w in the 1st PDF component.
393 0 : mu_w_1 = wm + sqrt( F_w * ( ( one - mixt_frac ) / mixt_frac ) * wp2 )
394 :
395 : ! Calculate the mean of w in the 2nd PDF component.
396 0 : mu_w_2 = wm - ( mixt_frac / ( one - mixt_frac ) ) * ( mu_w_1 - wm )
397 :
398 : ! Calculate the standard deviation of w in the 1st PDF component.
399 : ! sigma_w_1 = sqrt( ( ( zeta_w + 1 ) * ( 1 - F_w ) )
400 : ! / ( ( zeta_w + 2 ) * mixt_frac ) * <w'^2> )
401 : coef_sigma_w_1_sqd = ( ( zeta_w + one ) * ( one - F_w ) ) &
402 0 : / ( ( zeta_w + two ) * mixt_frac )
403 :
404 0 : sigma_w_1 = sqrt( coef_sigma_w_1_sqd * wp2 )
405 :
406 : ! Calculate the standard deviation of w in the 2nd PDF component.
407 : ! sigma_w_2 = sqrt( ( mixt_frac * sigma_w_1^2 )
408 : ! / ( ( 1 - mixt_frac ) * ( 1 + zeta_w ) ) )
409 : ! = sqrt( ( 1 - F_w )
410 : ! / ( ( zeta_w + 2 ) * ( 1 - mixt_frac ) ) * <w'^2> )
411 : coef_sigma_w_2_sqd = ( one - F_w ) &
412 0 : / ( ( zeta_w + two ) * ( one - mixt_frac ) )
413 :
414 0 : sigma_w_2 = sqrt( coef_sigma_w_2_sqd * wp2 )
415 :
416 : else ! mixt_frac <= 0 or mixt_frac >= 1
417 :
418 : ! The mixture fraction produced is invalid. This should only happen in
419 : ! the scenario where F_w = 0 and | Skw | > 0, where the value of
420 : ! mixt_frac has been set to -1. Set all output variables to 0 in this
421 : ! scenario. Since F_w is a function of skewness, the mixture fraction
422 : ! and the PDF should always be valid, and this section of code shouldn't
423 : ! be entered.
424 0 : mu_w_1 = zero
425 0 : mu_w_2 = zero
426 0 : sigma_w_1 = zero
427 0 : sigma_w_2 = zero
428 0 : coef_sigma_w_1_sqd = zero
429 0 : coef_sigma_w_2_sqd = zero
430 :
431 : endif ! 0 < mixt_frac < 1
432 :
433 :
434 0 : return
435 :
436 : end subroutine calculate_w_params
437 :
438 : !=============================================================================
439 : !
440 : ! DESCRIPTION OF THE METHOD FOR EACH RESPONDING VARIABLE
441 : ! ======================================================
442 : !
443 : ! In order to find equations for the four PDF parameters for each responding
444 : ! variable, which are mu_x_1, mu_x_2, sigma_x_1, and sigma_x_2 (where x stands
445 : ! for a responding variable here), four equations are needed. These four
446 : ! equations are the equations for <x>, <x'^2>, <x'^3>, and <w'x'> as found by
447 : ! integrating over the PDF. The four equations are:
448 : !
449 : ! <x> = mixt_frac * mu_x_1 + ( 1 - mixt_frac ) * mu_x_2;
450 : !
451 : ! <x'^2> = mixt_frac * ( ( mu_x_1 - <x> )^2 + sigma_x_1^2 )
452 : ! + ( 1 - mixt_frac ) * ( ( mu_x_2 - <x> )^2 + sigma_x_2^2 );
453 : !
454 : ! <x'^3> = mixt_frac * ( mu_x_1 - <x> )
455 : ! * ( ( mu_x_1 - <x> )^2 + 3 * sigma_x_1^2 )
456 : ! + ( 1 - mixt_frac ) * ( mu_x_2 - <x> )
457 : ! * ( ( mu_x_2 - <x> )^2 + 3 * sigma_x_2^2 ); and
458 : !
459 : ! <w'x'> = mixt_frac * ( mu_w_1 - <w> ) * ( mu_x_1 - <x> )
460 : ! + ( 1 - mixt_frac ) * ( mu_w_2 - <w> ) * ( mu_x_2 - <x> );
461 : !
462 : ! where the correlations that are normally found in the <w'x'> equation,
463 : ! corr_w_x_1 and corr_w_x_2, have both been set to 0.
464 : !
465 : ! The equations for mu_w_1 - <w> and mu_w_2 - <w> are:
466 : !
467 : ! mu_w_1 - <w> = sqrt( F_w * ( ( 1 - mixt_frac ) / mixt_frac ) * <w'^2> );
468 : !
469 : ! mu_w_2 - <w> = - sqrt( F_w * ( mixt_frac / ( 1 - mixt_frac ) ) * <w'^2> );
470 : !
471 : ! The resulting equations for the four PDF parameters are:
472 : !
473 : ! mu_x_1 = <x> + sqrt( ( 1 - mixt_frac ) / mixt_frac )
474 : ! * <w'x'> / sqrt( F_w * <w'^2> );
475 : !
476 : ! mu_x_2 = <x> - sqrt( mixt_frac / ( 1 - mixt_frac ) )
477 : ! * <w'x'> / sqrt( F_w * <w'^2> );
478 : !
479 : ! sigma_x_1^2 = ( 1 + sqrt( ( 1 - mixt_frac ) / mixt_frac )
480 : ! * Skx * sqrt( F_w * <w'^2> * <x'^2> ) / ( 3 * <w'x'> )
481 : ! - ( ( 1 + mixt_frac ) / mixt_frac )
482 : ! * <w'x'>^2 / ( 3 * F_w * <w'^2> * <x'^2> ) )
483 : ! * <x'^2>; and
484 : !
485 : ! sigma_x_2^2 = ( 1 - sqrt( mixt_frac / ( 1 - mixt_frac ) )
486 : ! * Skx * sqrt( F_w * <w'^2> * <x'^2> ) / ( 3 * <w'x'> )
487 : ! + ( ( mixt_frac - 2 ) / ( 1 - mixt_frac ) )
488 : ! * <w'x'>^2 / ( 3 * F_w * <w'^2> * <x'^2> ) )
489 : ! * <x'^2>;
490 : !
491 : ! where Skx is the skewness of x, and Skx = <x'^3> / <x'^2>^(3/2).
492 : !
493 : !
494 : ! Limits on F_w:
495 : !
496 : ! The only limits placed on the value of F_w from the w equation set itself
497 : ! are 0 <= F_w <= 1. However, use of the above equation set for responder
498 : ! variable x forces an additional limit to be placed on the value of F_w.
499 : ! That additional limit restricts the range of F_w to:
500 : !
501 : ! <w'x'>^2 / ( <w'^2> * <x'^2> ) <= F_w <= 1.
502 : !
503 : ! Furthermore, when there is more than one responder variable, F_w is limited
504 : ! by the most restrictive cases, such that:
505 : !
506 : ! max( <w'x'>^2 / ( <w'^2> * <x'^2> ), for all variables x ) <= F_w <= 1.
507 : !
508 : !
509 : ! Limits on Skx:
510 : !
511 : ! Since the PDF parameters for this variable need to work with the mixture
512 : ! fraction that has been provided by the setting variable, which is w, the
513 : ! method does not work for all values of Skx. However, the limits of Skx can
514 : ! always be calculated. The limits on Skw given by:
515 : !
516 : ! when <w'x'> > 0:
517 : !
518 : ! ( 1 + mixt_frac ) / sqrt( mixt_frac * ( 1 - mixt_frac ) )
519 : ! * <w'x'>^3 / ( F_w * <w'^2> * <x'^2> )^(3/2)
520 : ! - sqrt( mixt_frac / ( 1 - mixt_frac ) )
521 : ! * 3 * <w'x'> / sqrt( F_w * <w'^2> * <x'^2> )
522 : ! <= Skx <=
523 : ! ( mixt_frac - 2 ) / sqrt( mixt_frac * ( 1 - mixt_frac ) )
524 : ! * <w'x'>^3 / ( F_w * <w'^2> * <x'^2> )^(3/2)
525 : ! + sqrt( ( 1 - mixt_frac ) / mixt_frac )
526 : ! * 3 * <w'x'> / sqrt( F_w * <w'^2> * <x'^2> );
527 : !
528 : ! when <w'x'> < 0:
529 : !
530 : ! ( mixt_frac - 2 ) / sqrt( mixt_frac * ( 1 - mixt_frac ) )
531 : ! * <w'x'>^3 / ( F_w * <w'^2> * <x'^2> )^(3/2)
532 : ! + sqrt( ( 1 - mixt_frac ) / mixt_frac )
533 : ! * 3 * <w'x'> / sqrt( F_w * <w'^2> * <x'^2> )
534 : ! <= Skx <=
535 : ! ( 1 + mixt_frac ) / sqrt( mixt_frac * ( 1 - mixt_frac ) )
536 : ! * <w'x'>^3 / ( F_w * <w'^2> * <x'^2> )^(3/2)
537 : ! - sqrt( mixt_frac / ( 1 - mixt_frac ) )
538 : ! * 3 * <w'x'> / sqrt( F_w * <w'^2> * <x'^2> );
539 : !
540 : ! and when <w'x'> = 0, Skx = 0.
541 : !
542 : !
543 : ! Special cases:
544 : !
545 : ! When <w'x'> = 0, mu_x_1 = mu_x_2 = <x>, and the value of Skx must be 0.
546 : ! Since both <w'x'> = 0 and Skx = 0, the equations for sigma_x_1^2 and
547 : ! sigma_x_2^2 are both undefined. In this situation, the equations for the
548 : ! PDF parameters of x are:
549 : !
550 : ! mu_x_1 = mu_x_2 = <x>; and
551 : ! sigma_x_1^2 = sigma_x_2^2 = <x'^2>.
552 : !
553 : ! The value of F_w is allowed to be 0 only when <w'x'> = 0 (for all variables
554 : ! x). When <w'^2> = 0 and/or <x'^2> = 0, this means that <w'x'> = 0, as well.
555 : ! In all these situations, the equation set for the situation when <w'x'> = 0
556 : ! is used. This means that the distribution becomes a single Gaussian when
557 : ! <w'x'> = 0 (and Skx = 0).
558 : !
559 : !
560 : ! Notes:
561 : !
562 : ! The equations for the PDF component means and standard deviations can also
563 : ! be written as:
564 : !
565 : ! mu_x_1 = <x> + sqrt( ( 1 - mixt_frac ) / mixt_frac )
566 : ! * <w'x'> / sqrt( F_w * <w'^2> );
567 : !
568 : ! mu_x_2 = <x> - sqrt( mixt_frac / ( 1 - mixt_frac ) )
569 : ! * <w'x'> / sqrt( F_w * <w'^2> );
570 : !
571 : ! sigma_x_1 = sqrt( coef_sigma_x_1_sqd * <x'^2> ); and
572 : !
573 : ! sigma_x_2 = sqrt( coef_sigma_x_2_sqd * <x'^2> ); where
574 : !
575 : ! coef_sigma_x_1_sqd
576 : ! = 1 + sqrt( ( 1 - mixt_frac ) / mixt_frac )
577 : ! * Skx * sqrt( F_w * <w'^2> * <x'^2> ) / ( 3 * <w'x'> )
578 : ! - ( ( 1 + mixt_frac ) / mixt_frac )
579 : ! * <w'x'>^2 / ( 3 * F_w * <w'^2> * <x'^2> ); and
580 : !
581 : ! coef_sigma_x_2_sqd
582 : ! = 1 - sqrt( mixt_frac / ( 1 - mixt_frac ) )
583 : ! * Skx * sqrt( F_w * <w'^2> * <x'^2> ) / ( 3 * <w'x'> )
584 : ! + ( ( mixt_frac - 2 ) / ( 1 - mixt_frac ) )
585 : ! * <w'x'>^2 / ( 3 * F_w * <w'^2> * <x'^2> ).
586 : !
587 : ! The above equations can be substituted into an equation for a variable that
588 : ! has been derived by integrating over the PDF. Many variables like this are
589 : ! used in parts of the predictive equation set. These substitutions allow
590 : ! some terms to solved implicitly or semi-implicitly in the predictive
591 : ! equations.
592 : !
593 : !
594 : ! Brian Griffin; September 2019.
595 : !
596 : !=============================================================================
597 0 : elemental subroutine calculate_responder_params( xm, xp2, Skx, wpxp, & ! In
598 : wp2, F_w, mixt_frac, & ! In
599 : mu_x_1, mu_x_2, & ! Out
600 : sigma_x_1_sqd, & ! Out
601 : sigma_x_2_sqd, & ! Out
602 : coef_sigma_x_1_sqd, & ! Out
603 : coef_sigma_x_2_sqd ) ! Out
604 :
605 : ! Description:
606 : ! Calculates the PDF component means and the PDF component standard
607 : ! deviations for a responding variable (a variable that is not used to set
608 : ! the mixture fraction).
609 :
610 : ! References:
611 : ! Griffin and Larson (2020)
612 : !-----------------------------------------------------------------------
613 :
614 : use constants_clubb, only: &
615 : three, & ! Variable(s)
616 : two, &
617 : one, &
618 : zero
619 :
620 : use clubb_precision, only: &
621 : core_rknd ! Variable(s)
622 :
623 : implicit none
624 :
625 : ! Input Variables
626 : real( kind = core_rknd ), intent(in) :: &
627 : xm, & ! Mean of x (overall) [units vary]
628 : xp2, & ! Variance of x (overall) [(units vary)^2]
629 : Skx, & ! Skewness of x [-]
630 : wpxp, & ! Covariance of w and x (overall) [m/s(units vary)]
631 : wp2, & ! Variance of w (overall) [m^2/s^2]
632 : F_w, & ! Parameter for the spread of the PDF component means of w [-]
633 : mixt_frac ! Mixture fraction [-]
634 :
635 : ! Output Variables
636 : real( kind = core_rknd ), intent(out) :: &
637 : mu_x_1, & ! Mean of x (1st PDF component) [units vary]
638 : mu_x_2, & ! Mean of x (2nd PDF component) [units vary]
639 : sigma_x_1_sqd, & ! Variance of x (1st PDF component) [(units vary)^2]
640 : sigma_x_2_sqd ! Variance of x (2nd PDF component) [(units vary)^2]
641 :
642 : real( kind = core_rknd ), intent(out) :: &
643 : coef_sigma_x_1_sqd, & ! sigma_x_1^2 = coef_sigma_x_1_sqd * <x'^2> [-]
644 : coef_sigma_x_2_sqd ! sigma_x_2^2 = coef_sigma_x_2_sqd * <x'^2> [-]
645 :
646 :
647 0 : if ( abs( wpxp ) > zero ) then
648 :
649 : ! Note: when |<w'x'>| > 0, F_w, <w'^2>, and <x'^2> must all have values
650 : ! greater than 0.
651 :
652 : ! Calculate the mean of x in the 1st PDF component.
653 : mu_x_1 = xm + sqrt( ( one - mixt_frac ) / mixt_frac ) &
654 0 : * wpxp / sqrt( F_w * wp2 )
655 :
656 : ! Calculate the mean of x in the 2nd PDF component.
657 : mu_x_2 = xm - sqrt( mixt_frac / ( one - mixt_frac ) ) &
658 0 : * wpxp / sqrt( F_w * wp2 )
659 :
660 : ! Calculate the variance of x in the 1st PDF component.
661 : ! sigma_x_1^2
662 : ! = ( 1 + sqrt( ( 1 - mixt_frac ) / mixt_frac )
663 : ! * Skx * sqrt( F_w * <w'^2> * <x'^2> ) / ( 3 * <w'x'> )
664 : ! - ( ( 1 + mixt_frac ) / mixt_frac )
665 : ! * <w'x'>^2 / ( 3 * F_w * <w'^2> * <x'^2> ) )
666 : ! * <x'^2>
667 : coef_sigma_x_1_sqd &
668 : = one + sqrt( ( one - mixt_frac ) / mixt_frac ) &
669 : * Skx * sqrt( F_w * wp2 * xp2 ) / ( three * wpxp ) &
670 : - ( ( one + mixt_frac ) / mixt_frac ) &
671 0 : * wpxp**2 / ( three * F_w * wp2 * xp2 )
672 :
673 : ! Mathematically, the value of coef_sigma_x_1_sqd cannot be less than 0.
674 : ! Numerically, this can happen when numerical round off error causes an
675 : ! epsilon-sized negative value. When this happens, reset the value of
676 : ! coef_sigma_x_1_sqd to 0.
677 0 : coef_sigma_x_1_sqd = max( coef_sigma_x_1_sqd, zero )
678 :
679 0 : sigma_x_1_sqd = coef_sigma_x_1_sqd * xp2
680 :
681 : ! Calculate the variance of x in the 2nd PDF component.
682 : ! sigma_x_2^2
683 : ! = ( 1 - sqrt( mixt_frac / ( 1 - mixt_frac ) )
684 : ! * Skx * sqrt( F_w * <w'^2> * <x'^2> ) / ( 3 * <w'x'> )
685 : ! + ( ( mixt_frac - 2 ) / ( 1 - mixt_frac ) )
686 : ! * <w'x'>^2 / ( 3 * F_w * <w'^2> * <x'^2> ) )
687 : ! * <x'^2>
688 : coef_sigma_x_2_sqd &
689 : = one - sqrt( mixt_frac / ( one - mixt_frac ) ) &
690 : * Skx * sqrt( F_w * wp2 * xp2 ) / ( three * wpxp ) &
691 : + ( ( mixt_frac - two ) / ( one - mixt_frac ) ) &
692 0 : * wpxp**2 / ( three * F_w * wp2 * xp2 )
693 :
694 : ! Mathematically, the value of coef_sigma_x_2_sqd cannot be less than 0.
695 : ! Numerically, this can happen when numerical round off error causes an
696 : ! epsilon-sized negative value. When this happens, reset the value of
697 : ! coef_sigma_x_2_sqd to 0.
698 0 : coef_sigma_x_2_sqd = max( coef_sigma_x_2_sqd, zero )
699 :
700 0 : sigma_x_2_sqd = coef_sigma_x_2_sqd * xp2
701 :
702 : else ! | <w'x'> | = 0
703 :
704 : ! When <w'x'> has a value of 0, the PDF becomes a single Gaussian. This
705 : ! only works when Skx = 0. However, when Skx /= 0, the value of min_F_x
706 : ! is greater than 0, preventing a problem where F_x = 0 but | Skx | > 0.
707 0 : mu_x_1 = xm
708 0 : mu_x_2 = xm
709 0 : sigma_x_1_sqd = xp2
710 0 : sigma_x_2_sqd = xp2
711 0 : coef_sigma_x_1_sqd = one
712 0 : coef_sigma_x_2_sqd = one
713 :
714 : endif ! | <w'x'> | > 0
715 :
716 :
717 0 : return
718 :
719 : end subroutine calculate_responder_params
720 :
721 : !=============================================================================
722 0 : elemental function calculate_coef_wp4_implicit( mixt_frac, F_w, &
723 : coef_sigma_w_1_sqd, &
724 : coef_sigma_w_2_sqd ) &
725 : result( coef_wp4_implicit )
726 :
727 : ! Description:
728 : ! The predictive equation for <w'^3> contains a turbulent advection term of
729 : ! the form:
730 : !
731 : ! - ( 1 / rho_ds ) * d ( rho_ds * <w'^4> ) / dz;
732 : !
733 : ! where z is height, rho_ds is the dry, base-state density, and <w'^4> is
734 : ! calculated by integrating over the PDF. The equation for <w'^4> is:
735 : !
736 : ! <w'^4> = mixt_frac * ( 3 * sigma_w_1^4
737 : ! + 6 * ( mu_w_1 - <w> )^2 * sigma_w_1^2
738 : ! + ( mu_w_1 - <w> )^4 )
739 : ! + ( 1 - mixt_frac ) * ( 3 * sigma_w_2^4
740 : ! + 6 * ( mu_w_2 - <w> )^2 * sigma_w_2^2
741 : ! + ( mu_w_2 - <w> )^4 ).
742 : !
743 : ! The following substitutions are made into the above equation:
744 : !
745 : ! mu_w_1 - <w> = sqrt(F_w) * sqrt( ( 1 - mixt_frac ) / mixt_frac )
746 : ! * sqrt( <w'^2> );
747 : !
748 : ! mu_w_2 - <w> = - sqrt(F_w) * sqrt( mixt_frac / ( 1 - mixt_frac ) )
749 : ! * sqrt( <w'^2> );
750 : !
751 : ! sigma_w_1 = sqrt( coef_sigma_w_1_sqd * <w'^2> ); and
752 : !
753 : ! sigma_w_2 = sqrt( coef_sigma_w_2_sqd * <w'^2> ).
754 : !
755 : ! The equations for coef_sigma_w_1_sqd and coef_sigma_w_2_sqd are:
756 : !
757 : ! coef_sigma_w_1_sqd = ( ( zeta_w + 1 ) * ( 1 - F_w ) )
758 : ! / ( ( zeta_w + 2 ) * mixt_frac ); and
759 : !
760 : ! coef_sigma_w_2_sqd = ( 1 - F_w ) / ( ( zeta_w + 2 ) * ( 1 - mixt_frac ) ).
761 : !
762 : ! The equation for <w'4> becomes:
763 : !
764 : ! <w'^4> = ( 3 * mixt_frac * coef_sigma_w_1_sqd^2
765 : ! + 6 * F_w * ( 1 - mixt_frac ) * coef_sigma_w_1_sqd
766 : ! + F_w^2 * ( 1 - mixt_frac )^2 / mixt_frac
767 : ! + 3 * ( 1 - mixt_frac ) * coef_sigma_w_2_sqd^2
768 : ! + 6 * F_w * mixt_frac * coef_sigma_w_2_sqd
769 : ! + F_w^2 * mixt_frac^2 / ( 1 - mixt_frac ) ) * <w'^2>^2.
770 : !
771 : ! This equation is of the form:
772 : !
773 : ! <w'^4> = coef_wp4_implicit * <w'^2>^2;
774 : !
775 : ! where:
776 : !
777 : ! coef_wp4_implicit = 3 * mixt_frac * coef_sigma_w_1_sqd^2
778 : ! + 6 * F_w * ( 1 - mixt_frac ) * coef_sigma_w_1_sqd
779 : ! + F_w^2 * ( 1 - mixt_frac )^2 / mixt_frac
780 : ! + 3 * ( 1 - mixt_frac ) * coef_sigma_w_2_sqd^2
781 : ! + 6 * F_w * mixt_frac * coef_sigma_w_2_sqd
782 : ! + F_w^2 * mixt_frac^2 / ( 1 - mixt_frac ).
783 : !
784 : ! While the <w'^4> term is found in the <w'^3> predictive equation and not
785 : ! the <w'^2> predictive equation, the <w'^3> and <w'^2> predictive equations
786 : ! are solved together. This allows the term containing <w'^4> to be solved
787 : ! implicitly.
788 :
789 : ! References:
790 : !-----------------------------------------------------------------------
791 :
792 : use constants_clubb, only: &
793 : six, & ! Variable(s)
794 : three, &
795 : one
796 :
797 : use clubb_precision, only: &
798 : core_rknd ! Procedure(s)
799 :
800 : implicit none
801 :
802 : ! Input Variables
803 : real ( kind = core_rknd ), intent(in) :: &
804 : mixt_frac, & ! Mixture fraction [-]
805 : F_w, & ! Parameter: spread of the PDF comp. means of w [-]
806 : coef_sigma_w_1_sqd, & ! sigma_w_1^2 = coef_sigma_w_1_sqd * <w'^2> [-]
807 : coef_sigma_w_2_sqd ! sigma_w_2^2 = coef_sigma_w_2_sqd * <w'^2> [-]
808 :
809 : ! Return Variable
810 : real ( kind = core_rknd ) :: &
811 : coef_wp4_implicit ! Coef.: <w'^4> = coef_wp4_implicit * <w'^2>^2 [-]
812 :
813 :
814 : ! Calculate coef_wp4_implicit.
815 : coef_wp4_implicit = three * mixt_frac * coef_sigma_w_1_sqd**2 &
816 : + six * F_w * ( one - mixt_frac ) * coef_sigma_w_1_sqd &
817 : + F_w**2 * ( one - mixt_frac )**2 / mixt_frac &
818 : + three * ( one - mixt_frac ) * coef_sigma_w_2_sqd**2 &
819 : + six * F_w * mixt_frac * coef_sigma_w_2_sqd &
820 0 : + F_w**2 * mixt_frac**2 / ( one - mixt_frac )
821 :
822 :
823 : return
824 :
825 : end function calculate_coef_wp4_implicit
826 :
827 : !=============================================================================
828 0 : elemental function calc_coef_wp2xp_implicit( wp2, mixt_frac, F_w, &
829 : coef_sigma_w_1_sqd, &
830 : coef_sigma_w_2_sqd ) &
831 : result( coef_wp2xp_implicit )
832 :
833 : ! Description:
834 : ! The predictive equation for <w'x'> contains a turbulent advection term of
835 : ! the form:
836 : !
837 : ! - ( 1 / rho_ds ) * d ( rho_ds * <w'^2 x'> ) / dz;
838 : !
839 : ! where z is height, rho_ds is the dry, base-state density, and <w'^2 x'> is
840 : ! calculated by integrating over the PDF. The equation for <w'^2 x'> is:
841 : !
842 : ! <w'^2 x'>
843 : ! = mixt_frac * ( ( mu_x_1 - <x> ) * ( ( mu_w_1 - <w> )^2 + sigma_w_1^2 )
844 : ! + 2 * corr_w_x_1 * sigma_w_1 * sigma_x_1
845 : ! * ( mu_w_1 - <w> ) )
846 : ! + ( 1 - mixt_frac ) * ( ( mu_x_2 - <x> )
847 : ! * ( ( mu_w_2 - <w> )^2 + sigma_w_2^2 )
848 : ! + 2 * corr_w_x_2 * sigma_w_2 * sigma_x_2
849 : ! * ( mu_w_2 - <w> ) ).
850 : !
851 : ! The following substitutions are made into the above equation:
852 : !
853 : ! mu_w_1 - <w> = sqrt(F_w) * sqrt( ( 1 - mixt_frac ) / mixt_frac )
854 : ! * sqrt( <w'^2> );
855 : !
856 : ! mu_w_2 - <w> = - sqrt(F_w) * sqrt( mixt_frac / ( 1 - mixt_frac ) )
857 : ! * sqrt( <w'^2> );
858 : !
859 : ! mu_x_1 - <x> = sqrt( ( 1 - mixt_frac ) / mixt_frac )
860 : ! * <w'x'> / sqrt( F_w * <w'^2> );
861 : !
862 : ! mu_x_2 - <x> = - sqrt( mixt_frac / ( 1 - mixt_frac ) )
863 : ! * <w'x'> / sqrt( F_w * <w'^2> );
864 : !
865 : ! sigma_w_1 = sqrt( coef_sigma_w_1_sqd * <w'^2> );
866 : !
867 : ! sigma_w_2 = sqrt( coef_sigma_w_2_sqd * <w'^2> );
868 : !
869 : ! sigma_x_1 = sqrt( coef_sigma_x_1_sqd * <x'^2> ); and
870 : !
871 : ! sigma_x_2 = sqrt( coef_sigma_x_2_sqd * <x'^2> ).
872 : !
873 : ! The equations for coef_sigma_w_1_sqd and coef_sigma_w_2_sqd are:
874 : !
875 : ! coef_sigma_w_1_sqd = ( ( zeta_w + 1 ) * ( 1 - F_w ) )
876 : ! / ( ( zeta_w + 2 ) * mixt_frac ); and
877 : !
878 : ! coef_sigma_w_2_sqd = ( 1 - F_w ) / ( ( zeta_w + 2 ) * ( 1 - mixt_frac ) ).
879 : !
880 : ! The equations for coef_sigma_x_1_sqd and coef_sigma_x_2_sqd are:
881 : !
882 : ! coef_sigma_x_1_sqd
883 : ! = 1 + sqrt( ( 1 - mixt_frac ) / mixt_frac )
884 : ! * Skx * sqrt( F_w * <w'^2> * <x'^2> ) / ( 3 * <w'x'> )
885 : ! - ( ( 1 + mixt_frac ) / mixt_frac )
886 : ! * <w'x'>^2 / ( 3 * F_w * <w'^2> * <x'^2> ); and
887 : !
888 : ! coef_sigma_x_2_sqd
889 : ! = 1 - sqrt( mixt_frac / ( 1 - mixt_frac ) )
890 : ! * Skx * sqrt( F_w * <w'^2> * <x'^2> ) / ( 3 * <w'x'> )
891 : ! + ( ( mixt_frac - 2 ) / ( 1 - mixt_frac ) )
892 : ! * <w'x'>^2 / ( 3 * F_w * <w'^2> * <x'^2> ).
893 : !
894 : ! Additionally, corr_w_x_1 = corr_w_x_2 = 0.
895 : !
896 : ! The equation for <w'^2 x'> becomes:
897 : !
898 : ! <w'^2 x'> = sqrt( mixt_frac * ( 1 - mixt_frac ) )
899 : ! * ( F_w * ( ( 1 - mixt_frac ) / mixt_frac
900 : ! - mixt_frac / ( 1 - mixt_frac ) )
901 : ! + coef_sigma_w_1_sqd - coef_sigma_w_2_sqd )
902 : ! * sqrt( <w'^2> / F_w ) * <w'x'>.
903 : !
904 : ! This equation is of the form:
905 : !
906 : ! <w'^2 x'> = coef_wp2xp_implicit * <w'x'>;
907 : !
908 : ! where:
909 : !
910 : ! coef_wp2xp_implicit = sqrt( mixt_frac * ( 1 - mixt_frac ) )
911 : ! * ( F_w * ( ( 1 - mixt_frac ) / mixt_frac
912 : ! - mixt_frac / ( 1 - mixt_frac ) )
913 : ! + coef_sigma_w_1_sqd - coef_sigma_w_2_sqd )
914 : ! * sqrt( <w'^2> / F_w ).
915 : !
916 : ! In the special case that F_w = 0, <w'x'> must have a value of 0, in which
917 : ! case mu_x_1 - <x> = mu_x_2 - <x> = 0, and <w'^2 x'> = 0. The value of
918 : ! coef_wp2xp_implicit when F_w = 0 and <w'x'> = 0 can be calculated since
919 : ! Skw must also have a value of 0 when F_w = 0. When F_w = 0 and Skw = 0,
920 : ! mixt_frac = ( zeta_w + 1 ) / ( zeta_w + 2 ). When this happens,
921 : ! coef_sigma_w_1_sqd - coef_sigma_w_2_sqd = 0. The equation for
922 : ! coef_wp2xp_implicit becomes:
923 : !
924 : ! coef_wp2xp_implicit = sqrt( mixt_frac * ( 1 - mixt_frac ) )
925 : ! * ( ( 1 - mixt_frac ) / mixt_frac
926 : ! - mixt_frac / ( 1 - mixt_frac ) )
927 : ! * sqrt( F_w * <w'^2> );
928 : !
929 : ! and since F_w = 0, coef_wp2xp_implicit = 0.
930 :
931 : ! References:
932 : !-----------------------------------------------------------------------
933 :
934 : use constants_clubb, only: &
935 : one, & ! Variable(s)
936 : zero
937 :
938 : use clubb_precision, only: &
939 : core_rknd ! Procedure(s)
940 :
941 : implicit none
942 :
943 : ! Input Variables
944 : real ( kind = core_rknd ), intent(in) :: &
945 : wp2, & ! Variance of w (overall) [m^2/s^2]
946 : mixt_frac, & ! Mixture fraction [-]
947 : F_w, & ! Parameter: spread of the PDF comp. means of w [-]
948 : coef_sigma_w_1_sqd, & ! sigma_w_1^2 = coef_sigma_w_1_sqd * <w'^2> [-]
949 : coef_sigma_w_2_sqd ! sigma_w_2^2 = coef_sigma_w_2_sqd * <w'^2> [-]
950 :
951 : ! Return Variable
952 : ! Coefficient: <w'^2 x'> = coef_wp2xp_implicit * <w'x'>
953 : real ( kind = core_rknd ) :: &
954 : coef_wp2xp_implicit ! Coefficient that is multiplied by <w'x'> [m/s]
955 :
956 :
957 : ! Calculate coef_wp2xp_implicit.
958 0 : if ( F_w > 0 ) then
959 :
960 : coef_wp2xp_implicit &
961 : = sqrt( mixt_frac * ( one - mixt_frac ) ) &
962 : * ( F_w * ( ( one - mixt_frac ) / mixt_frac &
963 : - mixt_frac / ( one - mixt_frac ) ) &
964 : + coef_sigma_w_1_sqd - coef_sigma_w_2_sqd ) &
965 0 : * sqrt( wp2 / F_w )
966 :
967 : else ! F_w = 0
968 :
969 : coef_wp2xp_implicit = zero
970 :
971 : endif
972 :
973 :
974 : return
975 :
976 : end function calc_coef_wp2xp_implicit
977 :
978 : !=============================================================================
979 0 : elemental subroutine calc_coefs_wpxp2_semiimpl( wp2, wpxp, & ! In
980 : mixt_frac, F_w, & ! In
981 : coef_sigma_x_1_sqd, & ! In
982 : coef_sigma_x_2_sqd, & ! In
983 : coef_wpxp2_implicit, & ! Out
984 : term_wpxp2_explicit ) ! Out
985 :
986 : ! Description:
987 : ! The predictive equation for <x'^2> contains a turbulent advection term of
988 : ! the form:
989 : !
990 : ! - ( 1 / rho_ds ) * d ( rho_ds * <w'x'^2> ) / dz;
991 : !
992 : ! where z is height, rho_ds is the dry, base-state density, and <w'x'^2> is
993 : ! calculated by integrating over the PDF. The equation for <w'x'^2> is:
994 : !
995 : ! <w'x'^2>
996 : ! = mixt_frac * ( ( mu_w_1 - <w> ) * ( ( mu_x_1 - <x> )^2 + sigma_x_1^2 )
997 : ! + 2 * corr_w_x_1 * sigma_w_1 * sigma_x_1
998 : ! * ( mu_x_1 - <x> ) )
999 : ! + ( 1 - mixt_frac ) * ( ( mu_w_2 - <w> )
1000 : ! * ( ( mu_x_2 - <x> )^2 + sigma_x_2^2 )
1001 : ! + 2 * corr_w_x_2 * sigma_w_2 * sigma_x_2
1002 : ! * ( mu_x_2 - <x> ) ).
1003 : !
1004 : ! The following substitutions are made into the above equation:
1005 : !
1006 : ! mu_w_1 - <w> = sqrt(F_w) * sqrt( ( 1 - mixt_frac ) / mixt_frac )
1007 : ! * sqrt( <w'^2> );
1008 : !
1009 : ! mu_w_2 - <w> = - sqrt(F_w) * sqrt( mixt_frac / ( 1 - mixt_frac ) )
1010 : ! * sqrt( <w'^2> );
1011 : !
1012 : ! mu_x_1 - <x> = sqrt( ( 1 - mixt_frac ) / mixt_frac )
1013 : ! * <w'x'> / sqrt( F_w * <w'^2> );
1014 : !
1015 : ! mu_x_2 - <x> = - sqrt( mixt_frac / ( 1 - mixt_frac ) )
1016 : ! * <w'x'> / sqrt( F_w * <w'^2> );
1017 : !
1018 : ! sigma_w_1 = sqrt( coef_sigma_w_1_sqd * <w'^2> );
1019 : !
1020 : ! sigma_w_2 = sqrt( coef_sigma_w_2_sqd * <w'^2> );
1021 : !
1022 : ! sigma_x_1 = sqrt( coef_sigma_x_1_sqd * <x'^2> ); and
1023 : !
1024 : ! sigma_x_2 = sqrt( coef_sigma_x_2_sqd * <x'^2> ).
1025 : !
1026 : ! The equations for coef_sigma_w_1_sqd and coef_sigma_w_2_sqd are:
1027 : !
1028 : ! coef_sigma_w_1_sqd = ( ( zeta_w + 1 ) * ( 1 - F_w ) )
1029 : ! / ( ( zeta_w + 2 ) * mixt_frac ); and
1030 : !
1031 : ! coef_sigma_w_2_sqd = ( 1 - F_w ) / ( ( zeta_w + 2 ) * ( 1 - mixt_frac ) ).
1032 : !
1033 : ! The equations for coef_sigma_x_1_sqd and coef_sigma_x_2_sqd are:
1034 : !
1035 : ! coef_sigma_x_1_sqd
1036 : ! = 1 + sqrt( ( 1 - mixt_frac ) / mixt_frac )
1037 : ! * Skx * sqrt( F_w * <w'^2> * <x'^2> ) / ( 3 * <w'x'> )
1038 : ! - ( ( 1 + mixt_frac ) / mixt_frac )
1039 : ! * <w'x'>^2 / ( 3 * F_w * <w'^2> * <x'^2> ); and
1040 : !
1041 : ! coef_sigma_x_2_sqd
1042 : ! = 1 - sqrt( mixt_frac / ( 1 - mixt_frac ) )
1043 : ! * Skx * sqrt( F_w * <w'^2> * <x'^2> ) / ( 3 * <w'x'> )
1044 : ! + ( ( mixt_frac - 2 ) / ( 1 - mixt_frac ) )
1045 : ! * <w'x'>^2 / ( 3 * F_w * <w'^2> * <x'^2> ).
1046 : !
1047 : ! Additionally, corr_w_x_1 = corr_w_x_2 = 0.
1048 : !
1049 : ! The equation for <w'x'^2> becomes:
1050 : !
1051 : ! <w'x'^2>
1052 : ! = sqrt( mixt_frac * ( 1 - mixt_frac ) ) * sqrt( F_w * <w'^2> )
1053 : ! * ( coef_sigma_x_1_sqd - coef_sigma_x_2_sqd ) * <x'^2>
1054 : ! + sqrt( mixt_frac * ( 1 - mixt_frac ) )
1055 : ! * <w'x'>^2 / sqrt( F_w * <w'^2> )
1056 : ! * ( ( 1 - mixt_frac ) / mixt_frac - mixt_frac / ( 1 - mixt_frac ) ).
1057 : !
1058 : ! This equation is of the form:
1059 : !
1060 : ! <w'x'^2> = coef_wpxp2_implicit * <x'^2> + term_wpxp2_explicit;
1061 : !
1062 : ! where:
1063 : !
1064 : ! coef_wpxp2_implicit
1065 : ! = sqrt( mixt_frac * ( 1 - mixt_frac ) ) * sqrt( F_w * <w'^2> )
1066 : ! * ( coef_sigma_x_1_sqd - coef_sigma_x_2_sqd ); and
1067 : !
1068 : ! term_wpxp2_explicit
1069 : ! = sqrt( mixt_frac * ( 1 - mixt_frac ) ) * <w'x'>^2 / sqrt( F_w * <w'^2> )
1070 : ! * ( ( 1 - mixt_frac ) / mixt_frac - mixt_frac / ( 1 - mixt_frac ) ).
1071 : !
1072 : ! In the special case that F_w = 0, mu_w_1 - <w> = mu_w_2 - <w> = 0, and
1073 : ! <w'x'^2> = 0. The value of coef_wp2xp_implicit = 0 when F_w = 0.
1074 : ! Likewise, the value of <w'x'> = 0 when F_w = 0, which makes
1075 : ! term_wp2xp_explicit undefined in this scenario. However, since both
1076 : ! <w'x'^2> = 0 and coef_wp2xp_implicit = 0, term_wp2xp_explicit = 0 in this
1077 : ! situation.
1078 :
1079 : ! References:
1080 : !-----------------------------------------------------------------------
1081 :
1082 : use constants_clubb, only: &
1083 : one, & ! Variable(s)
1084 : zero
1085 :
1086 : use clubb_precision, only: &
1087 : core_rknd ! Procedure(s)
1088 :
1089 : implicit none
1090 :
1091 : ! Input Variables
1092 : real ( kind = core_rknd ), intent(in) :: &
1093 : wp2, & ! Variance of w (overall) [m^2/s^2]
1094 : wpxp, & ! Covariance of w and x [m/s (units vary)]
1095 : mixt_frac, & ! Mixture fraction [-]
1096 : F_w, & ! Parameter: spread of the PDF comp. means of w [-]
1097 : coef_sigma_x_1_sqd, & ! sigma_x_1^2 = coef_sigma_x_1_sqd * <x'^2> [-]
1098 : coef_sigma_x_2_sqd ! sigma_x_2^2 = coef_sigma_x_2_sqd * <x'^2> [-]
1099 :
1100 : ! Output Variables
1101 : ! Coefs.: <w'x'^2> = coef_wpxp2_implicit * <x'^2> + term_wpxp2_explicit
1102 : real ( kind = core_rknd ), intent(out) :: &
1103 : coef_wpxp2_implicit, & ! Coefficient that is multiplied by <x'^2> [m/s]
1104 : term_wpxp2_explicit ! Term that is on the RHS [m/s (units vary)^2]
1105 :
1106 :
1107 : ! Calculate coef_wpxp2_implicit and term_wpxp2_explicit.
1108 0 : if ( F_w > 0 .and. wp2 > 0 ) then
1109 :
1110 : coef_wpxp2_implicit &
1111 : = sqrt( mixt_frac * ( one - mixt_frac ) ) * sqrt( F_w * wp2 ) &
1112 0 : * ( coef_sigma_x_1_sqd - coef_sigma_x_2_sqd )
1113 :
1114 : term_wpxp2_explicit &
1115 : = sqrt( mixt_frac * ( one - mixt_frac ) ) * wpxp**2 / sqrt( F_w * wp2 ) &
1116 0 : * ( ( one - mixt_frac ) / mixt_frac - mixt_frac / ( one - mixt_frac ) )
1117 :
1118 : else ! F_w = 0 or wp2 = 0
1119 :
1120 0 : coef_wpxp2_implicit = zero
1121 0 : term_wpxp2_explicit = zero
1122 :
1123 : endif ! F_w > 0 and wp2 > 0
1124 :
1125 :
1126 0 : return
1127 :
1128 : end subroutine calc_coefs_wpxp2_semiimpl
1129 :
1130 : !=============================================================================
1131 0 : elemental subroutine calc_coefs_wpxpyp_semiimpl( wp2, wpxp, wpyp, & ! In
1132 : mixt_frac, F_w, & ! In
1133 : coef_sigma_x_1_sqd, & ! In
1134 : coef_sigma_x_2_sqd, & ! In
1135 : coef_sigma_y_1_sqd, & ! In
1136 : coef_sigma_y_2_sqd, & ! In
1137 : coef_wpxpyp_implicit, & ! Out
1138 : term_wpxpyp_explicit ) ! Out
1139 :
1140 : ! Description:
1141 : ! The predictive equation for <w'x'y'> contains a turbulent advection term
1142 : ! of the form:
1143 : !
1144 : ! - ( 1 / rho_ds ) * d ( rho_ds * <w'x'y'> ) / dz;
1145 : !
1146 : ! where z is height, rho_ds is the dry, base-state density, and <w'x'y'> is
1147 : ! calculated by integrating over the PDF. The equation for <w'x'y'> is:
1148 : !
1149 : ! <w'x'y'>
1150 : ! = mixt_frac
1151 : ! * ( ( mu_w_1 - <w> ) * ( mu_x_1 - <x> ) * ( mu_y_1 - <y> )
1152 : ! + corr_x_y_1 * sigma_x_1 * sigma_y_1 * ( mu_w_1 - <w> )
1153 : ! + corr_w_y_1 * sigma_w_1 * sigma_y_1 * ( mu_x_1 - <x> )
1154 : ! + corr_w_x_1 * sigma_w_1 * sigma_x_1 * ( mu_y_1 - <y> ) )
1155 : ! + ( 1 - mixt_frac )
1156 : ! * ( ( mu_w_2 - <w> ) * ( mu_x_2 - <x> ) * ( mu_y_2 - <y> )
1157 : ! + corr_x_y_2 * sigma_x_2 * sigma_y_2 * ( mu_w_2 - <w> )
1158 : ! + corr_w_y_2 * sigma_w_2 * sigma_y_2 * ( mu_x_2 - <x> )
1159 : ! + corr_w_x_2 * sigma_w_2 * sigma_x_2 * ( mu_y_2 - <y> ) ).
1160 : !
1161 : ! The following substitutions are made into the above equation:
1162 : !
1163 : ! mu_w_1 - <w> = sqrt(F_w) * sqrt( ( 1 - mixt_frac ) / mixt_frac )
1164 : ! * sqrt( <w'^2> );
1165 : !
1166 : ! mu_w_2 - <w> = - sqrt(F_w) * sqrt( mixt_frac / ( 1 - mixt_frac ) )
1167 : ! * sqrt( <w'^2> );
1168 : !
1169 : ! mu_x_1 - <x> = sqrt( ( 1 - mixt_frac ) / mixt_frac )
1170 : ! * <w'x'> / sqrt( F_w * <w'^2> );
1171 : !
1172 : ! mu_x_2 - <x> = - sqrt( mixt_frac / ( 1 - mixt_frac ) )
1173 : ! * <w'x'> / sqrt( F_w * <w'^2> );
1174 : !
1175 : ! mu_y_1 - <y> = sqrt( ( 1 - mixt_frac ) / mixt_frac )
1176 : ! * <w'y'> / sqrt( F_w * <w'^2> );
1177 : !
1178 : ! mu_y_2 - <y> = - sqrt( mixt_frac / ( 1 - mixt_frac ) )
1179 : ! * <w'y'> / sqrt( F_w * <w'^2> );
1180 : !
1181 : ! sigma_w_1 = sqrt( coef_sigma_w_1_sqd * <w'^2> );
1182 : !
1183 : ! sigma_w_2 = sqrt( coef_sigma_w_2_sqd * <w'^2> );
1184 : !
1185 : ! sigma_x_1 = sqrt( coef_sigma_x_1_sqd * <x'^2> );
1186 : !
1187 : ! sigma_x_2 = sqrt( coef_sigma_x_2_sqd * <x'^2> );
1188 : !
1189 : ! sigma_y_1 = sqrt( coef_sigma_y_1_sqd * <y'^2> ); and
1190 : !
1191 : ! sigma_y_2 = sqrt( coef_sigma_y_2_sqd * <y'^2> ).
1192 : !
1193 : ! The equations for coef_sigma_w_1_sqd and coef_sigma_w_2_sqd are:
1194 : !
1195 : ! coef_sigma_w_1_sqd = ( ( zeta_w + 1 ) * ( 1 - F_w ) )
1196 : ! / ( ( zeta_w + 2 ) * mixt_frac ); and
1197 : !
1198 : ! coef_sigma_w_2_sqd = ( 1 - F_w ) / ( ( zeta_w + 2 ) * ( 1 - mixt_frac ) ).
1199 : !
1200 : ! The equations for coef_sigma_x_1_sqd and coef_sigma_x_2_sqd are:
1201 : !
1202 : ! coef_sigma_x_1_sqd
1203 : ! = 1 + sqrt( ( 1 - mixt_frac ) / mixt_frac )
1204 : ! * Skx * sqrt( F_w * <w'^2> * <x'^2> ) / ( 3 * <w'x'> )
1205 : ! - ( ( 1 + mixt_frac ) / mixt_frac )
1206 : ! * <w'x'>^2 / ( 3 * F_w * <w'^2> * <x'^2> ); and
1207 : !
1208 : ! coef_sigma_x_2_sqd
1209 : ! = 1 - sqrt( mixt_frac / ( 1 - mixt_frac ) )
1210 : ! * Skx * sqrt( F_w * <w'^2> * <x'^2> ) / ( 3 * <w'x'> )
1211 : ! + ( ( mixt_frac - 2 ) / ( 1 - mixt_frac ) )
1212 : ! * <w'x'>^2 / ( 3 * F_w * <w'^2> * <x'^2> ).
1213 : !
1214 : ! The equations for coef_sigma_y_1_sqd and coef_sigma_y_2_sqd are:
1215 : !
1216 : ! coef_sigma_y_1_sqd
1217 : ! = 1 + sqrt( ( 1 - mixt_frac ) / mixt_frac )
1218 : ! * Sky * sqrt( F_w * <w'^2> * <y'^2> ) / ( 3 * <w'y'> )
1219 : ! - ( ( 1 + mixt_frac ) / mixt_frac )
1220 : ! * <w'y'>^2 / ( 3 * F_w * <w'^2> * <y'^2> ); and
1221 : !
1222 : ! coef_sigma_y_2_sqd
1223 : ! = 1 - sqrt( mixt_frac / ( 1 - mixt_frac ) )
1224 : ! * Sky * sqrt( F_w * <w'^2> * <y'^2> ) / ( 3 * <w'y'> )
1225 : ! + ( ( mixt_frac - 2 ) / ( 1 - mixt_frac ) )
1226 : ! * <w'y'>^2 / ( 3 * F_w * <w'^2> * <y'^2> ).
1227 : !
1228 : ! Additionally, corr_w_x_1 = corr_w_x_2 = corr_w_y_1 = corr_w_y_2 = 0; and:
1229 : !
1230 : ! corr_x_y_1 = corr_x_y_2
1231 : ! = ( <x'y'> - mixt_frac * ( mu_x_1 - <x> ) * ( mu_y_1 - <y> )
1232 : ! - ( 1 - mixt_frac ) * ( mu_x_2 - <x> ) * ( mu_y_2 - <y> ) )
1233 : ! / ( mixt_frac * sigma_x_1 * sigma_y_1
1234 : ! + ( 1 - mixt_frac ) * sigma_x_2 * sigma_y_2 );
1235 : !
1236 : ! where -1 <= corr_x_y_1 = corr_x_y_2 <= 1. This equation can be rewritten
1237 : ! as:
1238 : !
1239 : ! corr_x_y_1 = corr_x_y_2
1240 : ! = ( <x'y'>
1241 : ! - sqrt( F_x ) * sqrt( F_y ) * sgn( <w'x'> ) * sgn( <w'y'> )
1242 : ! * sqrt( <x'^2> ) * sqrt( <y'^2 > ) )
1243 : ! / ( ( mixt_frac * sqrt( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd )
1244 : ! + ( 1 - mixt_frac )
1245 : ! * sqrt( coef_sigma_x_2_sqd * coef_sigma_y_2_sqd ) )
1246 : ! * sqrt( <x'^2> ) * sqrt( <y'^2> ) ).
1247 : !
1248 : ! The equation for <w'x'y'> becomes:
1249 : !
1250 : ! <w'x'y'>
1251 : ! = sqrt( mixt_frac * ( 1 - mixt_frac ) ) * sqrt( F_w * <w'^2> )
1252 : ! * ( sqrt( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd )
1253 : ! - sqrt( coef_sigma_x_2_sqd * coef_sigma_y_2_sqd ) )
1254 : ! / ( mixt_frac * sqrt( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd )
1255 : ! + ( 1 - mixt_frac )
1256 : ! * sqrt( coef_sigma_x_2_sqd * coef_sigma_y_2_sqd ) )
1257 : ! * <x'y'>
1258 : ! + sqrt( mixt_frac * ( 1 - mixt_frac ) )
1259 : ! * <w'x'> * <w'y'> / sqrt( F_w * <w'^2> )
1260 : ! * ( ( 1 - mixt_frac ) / mixt_frac - mixt_frac / ( 1 - mixt_frac )
1261 : ! - ( sqrt( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd )
1262 : ! - sqrt( coef_sigma_x_2_sqd * coef_sigma_y_2_sqd ) )
1263 : ! / ( mixt_frac * sqrt( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd )
1264 : ! + ( 1 - mixt_frac )
1265 : ! * sqrt( coef_sigma_x_2_sqd * coef_sigma_y_2_sqd ) ) )
1266 : !
1267 : ! This equation is of the form:
1268 : !
1269 : ! <w'x'y'> = coef_wpxpyp_implicit * <x'y'> + term_wpxpyp_explicit;
1270 : !
1271 : ! where:
1272 : !
1273 : ! coef_wpxpyp_implicit
1274 : ! = sqrt( mixt_frac * ( 1 - mixt_frac ) ) * sqrt( F_w * <w'^2> )
1275 : ! * ( sqrt( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd )
1276 : ! - sqrt( coef_sigma_x_2_sqd * coef_sigma_y_2_sqd ) )
1277 : ! / ( mixt_frac * sqrt( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd )
1278 : ! + ( 1 - mixt_frac )
1279 : ! * sqrt( coef_sigma_x_2_sqd * coef_sigma_y_2_sqd ) ); and
1280 : !
1281 : ! term_wpxpyp_explicit
1282 : ! = sqrt( mixt_frac * ( 1 - mixt_frac ) )
1283 : ! * <w'x'> * <w'y'> / sqrt( F_w * <w'^2> )
1284 : ! * ( ( 1 - mixt_frac ) / mixt_frac - mixt_frac / ( 1 - mixt_frac )
1285 : ! - ( sqrt( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd )
1286 : ! - sqrt( coef_sigma_x_2_sqd * coef_sigma_y_2_sqd ) )
1287 : ! / ( mixt_frac * sqrt( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd )
1288 : ! + ( 1 - mixt_frac )
1289 : ! * sqrt( coef_sigma_x_2_sqd * coef_sigma_y_2_sqd ) ) ).
1290 : !
1291 : ! There are also special cases for the above equations. In the scenario
1292 : ! that sigma_x_1 * sigma_y_1 = sigma_x_2 * sigma_y_2 = 0, and equations for
1293 : ! coef_wpxpyp_implicit and term_wpxpyp_explicit become:
1294 : !
1295 : ! coef_wpxpyp_implicit
1296 : ! = sqrt( mixt_frac * ( 1 - mixt_frac ) ) * sqrt( F_w * <w'^2> )
1297 : ! * ( ( 1 - mixt_frac ) / mixt_frac - mixt_frac / ( 1 - mixt_frac ) ); and
1298 : !
1299 : ! term_wpxpyp_explicit = 0.
1300 : !
1301 : ! In the scenario where F_w = 0, mu_w_1 - <w> = mu_w_2 - <w> = 0, and
1302 : ! <w'x'> = <w'y'> = 0, which means mu_x_1 - <x> = mu_x_2 - <x> = 0, as well
1303 : ! as mu_y_1 - <y> = mu_y_2 - <y> = 0, and as a result, <w'x'y'> = 0. When
1304 : ! F_w = 0, coef_wpxpyp_implicit = 0. Since <w'x'y'> = 0 and
1305 : ! coef_wpxpyp_implicit = 0 when F_w = 0, term_wpxpyp_explicit = 0.
1306 :
1307 : ! References:
1308 : !-----------------------------------------------------------------------
1309 :
1310 : use constants_clubb, only: &
1311 : one, & ! Variable(s)
1312 : zero
1313 :
1314 : use clubb_precision, only: &
1315 : core_rknd ! Procedure(s)
1316 :
1317 : implicit none
1318 :
1319 : ! Input Variables
1320 : real ( kind = core_rknd ), intent(in) :: &
1321 : wp2, & ! Variance of w (overall) [m^2/s^2]
1322 : wpxp, & ! Covariance of w and x (overall) [m/s(x units)]
1323 : wpyp, & ! Covariance of w and y (overall) [m/s(y units)]
1324 : mixt_frac, & ! Mixture fraction [-]
1325 : F_w, & ! Parameter: spread of the PDF comp. means of w [-]
1326 : coef_sigma_x_1_sqd, & ! sigma_x_1^2 = coef_sigma_x_1_sqd * <x'^2> [-]
1327 : coef_sigma_x_2_sqd, & ! sigma_x_2^2 = coef_sigma_x_2_sqd * <x'^2> [-]
1328 : coef_sigma_y_1_sqd, & ! sigma_y_1^2 = coef_sigma_y_1_sqd * <y'^2> [-]
1329 : coef_sigma_y_2_sqd ! sigma_y_2^2 = coef_sigma_y_2_sqd * <y'^2> [-]
1330 :
1331 : ! Output Variables
1332 : ! Coefs.: <w'x'y'> = coef_wpxpyp_implicit * <x'y'> + term_wpxpyp_explicit
1333 : real ( kind = core_rknd ), intent(out) :: &
1334 : coef_wpxpyp_implicit, & ! Coefficient that is multiplied by <x'y'> [m/s]
1335 : term_wpxpyp_explicit ! Term that is on the RHS [m/s(x units)(y units)]
1336 :
1337 : ! Local Variables
1338 : real ( kind = core_rknd ) :: &
1339 : coefs_factor_xy ! Factor involving coef_sigma_... x and y coefs [-]
1340 :
1341 :
1342 : ! Calculate coef_wpxpyp_implicit and term_wpxpyp_explicit.
1343 : if ( ( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd > zero &
1344 : .or. coef_sigma_x_2_sqd * coef_sigma_y_2_sqd > zero ) &
1345 0 : .and. F_w > zero .and. wp2 > zero ) then
1346 :
1347 : ! coefs_factor_xy
1348 : ! = ( sqrt( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd )
1349 : ! - sqrt( coef_sigma_x_2_sqd * coef_sigma_y_2_sqd ) )
1350 : ! / ( mixt_frac * sqrt( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd )
1351 : ! + ( 1 - mixt_frac )
1352 : ! * sqrt( coef_sigma_x_2_sqd * coef_sigma_y_2_sqd ) )
1353 : coefs_factor_xy &
1354 : = ( sqrt( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd ) &
1355 : - sqrt( coef_sigma_x_2_sqd * coef_sigma_y_2_sqd ) ) &
1356 : / ( mixt_frac * sqrt( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd ) &
1357 : + ( one - mixt_frac ) &
1358 0 : * sqrt( coef_sigma_x_2_sqd * coef_sigma_y_2_sqd ) )
1359 :
1360 : coef_wpxpyp_implicit &
1361 : = sqrt( mixt_frac * ( 1 - mixt_frac ) ) &
1362 0 : * sqrt( F_w * wp2 ) * coefs_factor_xy
1363 :
1364 : term_wpxpyp_explicit &
1365 : = sqrt( mixt_frac * ( one - mixt_frac ) ) &
1366 : * wpxp * wpyp / sqrt( F_w * wp2 ) &
1367 : * ( ( one - mixt_frac ) / mixt_frac &
1368 : - mixt_frac / ( one - mixt_frac ) &
1369 0 : - coefs_factor_xy )
1370 :
1371 : else ! ( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd = 0
1372 : ! and coef_sigma_x_2_sqd * coef_sigma_y_2_sqd = 0 )
1373 : ! or F_w = 0 or wp2 = 0
1374 :
1375 0 : if ( F_w > 0 ) then
1376 :
1377 : coef_wpxpyp_implicit &
1378 : = sqrt( mixt_frac * ( one - mixt_frac ) ) * sqrt( F_w * wp2 ) &
1379 : * ( ( one - mixt_frac ) / mixt_frac &
1380 0 : - mixt_frac / ( one - mixt_frac ) )
1381 :
1382 : else ! F_w = 0
1383 :
1384 0 : coef_wpxpyp_implicit = zero
1385 :
1386 : endif
1387 :
1388 0 : term_wpxpyp_explicit = zero
1389 :
1390 : endif
1391 :
1392 :
1393 0 : return
1394 :
1395 : end subroutine calc_coefs_wpxpyp_semiimpl
1396 :
1397 : !=============================================================================
1398 :
1399 : end module new_hybrid_pdf
|
