Line data Source code
1 : ! $Id$
2 : !===============================================================================
3 : module new_pdf
4 :
5 : ! Description:
6 : ! The portion of CLUBB's multivariate, two-component PDF that is the
7 : ! trivariate, two-component normal PDF of vertical velocity (w), total water
8 : ! mixing ratio (rt), and liquid water potential temperature (thl).
9 :
10 : ! References:
11 : ! Griffin and Larson (2018)
12 : !-------------------------------------------------------------------------
13 :
14 : implicit none
15 :
16 : public :: calc_mixture_fraction, & ! Procedure(s)
17 : calc_setter_var_params, &
18 : calc_responder_params, &
19 : calc_limits_F_x_responder, &
20 : calc_coef_wp4_implicit, &
21 : calc_coef_wpxp2_implicit, &
22 : calc_coefs_wp2xp_semiimpl, &
23 : calc_coefs_wpxpyp_semiimpl
24 :
25 : private :: sort_roots ! Procedure(s)
26 :
27 : private
28 :
29 : contains
30 :
31 : !=============================================================================
32 : !
33 : ! DESCRIPTION OF THE METHOD FOR THE VARIABLE THAT SETS THE MIXTURE FRACTION
34 : ! =========================================================================
35 : !
36 : ! Many times, w has been used as the variable that sets the mixture fraction
37 : ! for the PDF. There are five PDF parameters that need to be calculated,
38 : ! which are mu_w_1 (the mean of w is the 1st PDF component), mu_w_2 (the mean
39 : ! of w in the 2nd PDF component), sigma_w_1 (the standard deviation of w in
40 : ! the 1st PDF component), sigma_w_2 (the standard deviation of w in the 2nd
41 : ! PDF component), and mixt_frac (the mixture fraction, which is the weight of
42 : ! the 1st PDF component). In order to solve for these five parameters, five
43 : ! equations are needed. These five equations are the equations for <w>,
44 : ! <w'^2>, and <w'^3> as found by integrating over the PDF. Additionally, two
45 : ! more equations, which involve tunable parameters F_w and zeta_w, and which
46 : ! are used to help control the spread of the PDF component means and the size
47 : ! of the PDF component standard deviations compared to each other,
48 : ! respectively, are used in this equation set. The five equations are:
49 : !
50 : ! <w> = mixt_frac * mu_w_1 + ( 1 - mixt_frac ) * mu_w_2;
51 : !
52 : ! <w'^2> = mixt_frac * ( ( mu_w_1 - <w> )^2 + sigma_w_1^2 )
53 : ! + ( 1 - mixt_frac ) * ( ( mu_w_2 - <w> )^2 + sigma_w_2^2 );
54 : !
55 : ! <w'^3> = mixt_frac * ( mu_w_1 - <w> )
56 : ! * ( ( mu_w_1 - <w> )^2 + 3 * sigma_w_1^2 )
57 : ! + ( 1 - mixt_frac ) * ( mu_w_2 - <w> )
58 : ! * ( ( mu_w_2 - <w> )^2 + 3 * sigma_w_2^2 );
59 : !
60 : ! mu_w_1 - <w> = sqrt(F_w) * ( sqrt( 1 - mixt_frac ) / sqrt( mixt_frac ) )
61 : ! * sqrt( <w'^2> );
62 : !
63 : ! where 0 <= F_w <= 1; and
64 : !
65 : ! 1 + zeta_w = ( mixt_frac * sigma_w_1^2 )
66 : ! / ( ( 1 - mixt_frac ) * sigma_w_2^2 );
67 : !
68 : ! where zeta_w > -1.
69 : !
70 : ! Following convention for w, mu_w_1 is defined to be greater than or equal to
71 : ! mu_w_2 (and is also greater than or equal to <w>, while mu_w_2 is less than
72 : ! or equal to <w>). This is relationship is found in the mu_w_1 - <w>
73 : ! equation above.
74 : !
75 : ! The resulting equations for the five PDF parameters are:
76 : !
77 : ! mixt_frac
78 : ! = ( 4 * F_w^3
79 : ! + 18 * F_w * ( zeta_w + 1 ) * ( 1 - F_w ) / ( zeta_w + 2 )
80 : ! + 6 * F_w^2 * ( 1 - F_w ) / ( zeta_w + 2 )
81 : ! + Skw^2
82 : ! - Skw * sqrt( 4 * F_w^3
83 : ! + 12 * F_w^2 * ( 1 - F_w )
84 : ! + 36 * F_w * ( zeta_w + 1 ) * ( 1 - F_w )^2
85 : ! / ( zeta_w + 2 )^2
86 : ! + Skw^2 ) )
87 : ! / ( 2 * F_w * ( F_w - 3 )^2 + 2 * Skw^2 );
88 : !
89 : ! mu_w_1 = <w> + sqrt( F_w * ( ( 1 - mixt_frac ) / mixt_frac ) * <w'^2> );
90 : !
91 : ! mu_w_2 = <w> - ( mixt_frac / ( 1 - mixt_frac ) ) * ( mu_w_1 - <w> );
92 : !
93 : ! sigma_w_1
94 : ! = sqrt( ( ( zeta_w + 1 ) * ( 1 - F_w ) )
95 : ! / ( ( zeta_w + 2 ) * mixt_frac ) * <w'^2> ); and
96 : !
97 : ! sigma_w_2
98 : ! = sqrt( ( mixt_frac * sigma_w_1^2 )
99 : ! / ( ( 1 - mixt_frac ) * ( 1 + zeta_w ) ) );
100 : !
101 : ! where Skw is the skewness of w, and Skw = <w'^3> / <w'^2>^(3/2).
102 : !
103 : ! This method works for all values of F_w (where 0 <= F_w <= 1) and zeta_w
104 : ! (where zeta_w > -1).
105 : !
106 : !
107 : ! Generalized equations for any variable, x, that sets the mixture fraction:
108 : !
109 : ! A slight alteration is made to the above equations in order to have any
110 : ! variable, x, set the mixture fraction. The same five PDF parameters need to
111 : ! be calculated for the setting variable, which are mu_x_1 (the mean of x in
112 : ! the 1st PDF component), mu_x_2 (the mean of x in the 2nd PDF component),
113 : ! sigma_x_1 (the standard deviation of x in the 1st PDF component), sigma_x_2
114 : ! (the standard deviation of x in the 2nd PDF component), and mixt_frac (the
115 : ! mixture fraction). Again, five equations are needed, and they are the
116 : ! equations for <x>, <x'^2>, and <x'^3> as found by integrating over the PDF,
117 : ! as well as the equations that involve tunable parameters F_x and zeta_x.
118 : ! However, the equation for F_x is multiplied by a new variable,
119 : ! sgn( <w'x'> ), where <w'x'> is the covariance of w and x, and sgn( <w'x'> )
120 : ! is given by:
121 : !
122 : ! sgn( <w'x'> ) = | 1, when <w'x'> >= 0;
123 : ! | -1, when <w'x'> < 0.
124 : !
125 : ! The five equations are:
126 : !
127 : ! <x> = mixt_frac * mu_x_1 + ( 1 - mixt_frac ) * mu_x_2;
128 : !
129 : ! <x'^2> = mixt_frac * ( ( mu_x_1 - <x> )^2 + sigma_x_1^2 )
130 : ! + ( 1 - mixt_frac ) * ( ( mu_x_2 - <x> )^2 + sigma_x_2^2 );
131 : !
132 : ! <x'^3> = mixt_frac * ( mu_x_1 - <x> )
133 : ! * ( ( mu_x_1 - <x> )^2 + 3 * sigma_x_1^2 )
134 : ! + ( 1 - mixt_frac ) * ( mu_x_2 - <x> )
135 : ! * ( ( mu_x_2 - <x> )^2 + 3 * sigma_x_2^2 );
136 : !
137 : ! mu_x_1 - <x> = sqrt(F_x) * ( sqrt( 1 - mixt_frac ) / sqrt( mixt_frac ) )
138 : ! * sqrt( <x'^2> ) * sgn( <w'x'> );
139 : !
140 : ! where 0 <= F_x <= 1; and
141 : !
142 : ! 1 + zeta_x = ( mixt_frac * sigma_x_1^2 )
143 : ! / ( ( 1 - mixt_frac ) * sigma_x_2^2 );
144 : !
145 : ! where zeta_x > -1.
146 : !
147 : ! The only equations that are altered are the equation for mu_x_1 and the
148 : ! equation for mixt_frac, which now both contain a sgn( <w'x'> ). The mu_x_2
149 : ! equation is not altered, but the sign of mu_x_2 - <x> will be the opposite
150 : ! of the sign of mu_x_1 - <x>. The resulting equations for the five PDF
151 : ! parameters are:
152 : !
153 : ! mixt_frac
154 : ! = ( 4 * F_x^3
155 : ! + 18 * F_x * ( zeta_x + 1 ) * ( 1 - F_x ) / ( zeta_x + 2 )
156 : ! + 6 * F_x^2 * ( 1 - F_x ) / ( zeta_x + 2 )
157 : ! + Skx^2
158 : ! - Skx * sgn( <w'x'> )
159 : ! * sqrt( 4 * F_x^3
160 : ! + 12 * F_x^2 * ( 1 - F_x )
161 : ! + 36 * F_x * ( zeta_x + 1 ) * ( 1 - F_x )^2
162 : ! / ( zeta_x + 2 )^2
163 : ! + Skx^2 ) )
164 : ! / ( 2 * F_x * ( F_x - 3 )^2 + 2 * Skx^2 );
165 : !
166 : ! mu_x_1 = <x> + sqrt( F_x * ( ( 1 - mixt_frac ) / mixt_frac ) * <x'^2> )
167 : ! * sgn( <w'x'> );
168 : !
169 : ! mu_x_2 = <x> - ( mixt_frac / ( 1 - mixt_frac ) ) * ( mu_x_1 - <x> );
170 : !
171 : ! sigma_x_1
172 : ! = sqrt( ( ( zeta_x + 1 ) * ( 1 - F_x ) )
173 : ! / ( ( zeta_x + 2 ) * mixt_frac ) * <x'^2> ); and
174 : !
175 : ! sigma_x_2
176 : ! = sqrt( ( mixt_frac * sigma_x_1^2 )
177 : ! / ( ( 1 - mixt_frac ) * ( 1 + zeta_x ) ) );
178 : !
179 : ! where Skx is the skewness of x, and Skx = <x'^3> / <x'^2>^(3/2).
180 : !
181 : ! This method works for all values of F_x (where 0 <= F_x <= 1) and zeta_x
182 : ! (where zeta_x > -1).
183 : !
184 : ! When the generalized form is solved for w (x = w), sgn( <w'^2> ) = 1, and
185 : ! the equations are unaltered from the equations listed above for w.
186 : !
187 : !
188 : ! Special case:
189 : !
190 : ! When Skx = 0 and F_x = 0, the equation for mixt_frac is undefined. The
191 : ! equation for mixture fraction in this scenario can be derived by using the
192 : ! above equation for mixture fraction and then setting Skx = 0. The resulting
193 : ! equation becomes:
194 : !
195 : ! mixt_frac
196 : ! = ( 4 * F_x^3
197 : ! + 18 * F_x * ( zeta_x + 1 ) * ( 1 - F_x ) / ( zeta_x + 2 )
198 : ! + 6 * F_x^2 * ( 1 - F_x ) / ( zeta_x + 2 ) )
199 : ! / ( 2 * F_x * ( F_x - 3 )^2 ).
200 : !
201 : ! All of the terms in the numerator and denominator contain a F_x, so this
202 : ! equation can be rewritten as:
203 : !
204 : ! mixt_frac
205 : ! = ( 4 * F_x^2
206 : ! + 18 * ( zeta_x + 1 ) * ( 1 - F_x ) / ( zeta_x + 2 )
207 : ! + 6 * F_x * ( 1 - F_x ) / ( zeta_x + 2 ) )
208 : ! / ( 2 * ( F_x - 3 )^2 ).
209 : !
210 : ! Now setting F_x = 0, the equation becomes:
211 : !
212 : ! mixt_frac = ( 18 * ( zeta_x + 1 ) / ( zeta_x + 2 ) ) / 18;
213 : !
214 : ! which can be rewritten as:
215 : !
216 : ! mixt_frac = ( zeta_x + 1 ) / ( zeta_x + 2 ).
217 : !
218 : ! When F_x = 0, Skx must have a value of 0 in order for the PDF to function
219 : ! correctly. When F_x = 0, mu_x_1 = mu_x_2. When the two PDF component means
220 : ! are equal to each other (and to the overall mean, <x>), the only value of
221 : ! Skx that can be represented is a value of 0. In the equation for mixture
222 : ! fraction, when F_x is set to 0, but | Skx | > 0, mixt_frac will either have
223 : ! a value of 0 or 1, depending on whether Skx is positive or negative,
224 : ! respectively.
225 : !
226 : ! The value of F_x should be set as a function of Skx. The value F_x should
227 : ! go toward 0 as | Skx | (or Skx^2) goes toward 0. The value of F_x should
228 : ! go toward 1 as | Skx | (or Skx^2) goes to infinity.
229 : !
230 : !
231 : ! Tunable parameters:
232 : !
233 : ! 1) F_x: This parameter controls the spread of the PDF component means. The
234 : ! range of this parameter is 0 <= F_x <= 1. When F_x = 0, the two
235 : ! PDF component means (mu_x_1 and mu_x_2) are equal to each other
236 : ! (and Skx must equal 0). All of the variance of x is accounted for
237 : ! by the PDF component standard deviations (sigma_x_1 and sigma_x_2).
238 : ! When F_x = 1, mu_x_1 and mu_x_2 are spread as far apart as they can
239 : ! be. Both PDF component standard deviations (sigma_x_1 and
240 : ! sigma_x_2) are equal to 0, and all of the variance of x is
241 : ! accounted for by the spread of the PDF component means.
242 : !
243 : ! When sigma_x_1 = sigma_x_2 = 0, the equation for <x'^2> becomes:
244 : !
245 : ! <x'^2> = mixt_frac * ( mu_x_1 - <x> )^2
246 : ! + ( 1 - mixt_frac ) * ( mu_x_2 - <x> )^2.
247 : !
248 : ! Substituting the equation for <x> into the above equation for
249 : ! mu_x_2 - <x>, the above equation becomes:
250 : !
251 : ! <x'^2> = ( mixt_frac / ( 1 - mixt_frac ) ) * ( mu_x_1 - <x> )^2;
252 : !
253 : ! which can be rewritten as:
254 : !
255 : ! ( mu_x_1 - <x> )^2 = ( ( 1 - mixt_frac ) / mixt_frac ) * <x'^2>.
256 : !
257 : ! Taking the square root of the above equation:
258 : !
259 : ! mu_x_1 - <x> = +/- ( sqrt( 1 - mixt_frac ) / sqrt(mixt_frac) )
260 : ! * sqrt( <x'^2> ).
261 : !
262 : ! This equation can be compared to the equation for mu_x_1 - <x> in
263 : ! the set of 5 equations, which is:
264 : !
265 : ! mu_x_1 - <x>
266 : ! = sqrt(F_x) * ( sqrt( 1 - mixt_frac ) / sqrt( mixt_frac ) )
267 : ! * sqrt( <x'^2> ) * sgn( <w'x'> ).
268 : !
269 : ! The above equations give another example of the meaning of F_x.
270 : ! The value of sqrt(F_x) is ratio of mu_x_1 - <x> to its maximum
271 : ! value (or minimum value, depending on sgn( <w'x'> )), which is:
272 : !
273 : ! sqrt( ( ( 1 - mixt_frac ) / mixt_frac ) * <x'^2> ) * sgn( <w'x'> ).
274 : !
275 : !
276 : ! 2) zeta_x: This parameter controls the size of the PDF component standard
277 : ! deviations compared to each other. The equation for zeta_x is:
278 : !
279 : ! 1 + zeta_x = ( mixt_frac * sigma_x_1^2 )
280 : ! / ( ( 1 - mixt_frac ) * sigma_x_2^2 ).
281 : !
282 : ! When zeta_x > 0, mixt_frac * sigma_x_1^2 increases at the
283 : ! expense of ( 1 - mixt_frac ) * sigma_x_2^2, which decreases in
284 : ! this variance-preserving equation set. When zeta_x = 0, then
285 : ! mixt_frac * sigma_x_1^2 = ( 1 - mixt_frac ) * sigma_x_2^2.
286 : ! When -1 < zeta_x < 0, ( 1 - mixt_frac ) * sigma_x_2^2 increases
287 : ! at the expense of mixt_frac * sigma_x_1^2, which decreases. As
288 : ! a result, greater values of zeta_x cause the 1st PDF component
289 : ! to become broader while the 2nd PDF component becomes narrower,
290 : ! and smaller values of zeta_x cause the 1st PDF component to
291 : ! become narrower while the 2nd PDF component becomes broader.
292 : !
293 : ! Symmetry
294 : !
295 : ! When zeta_x = 0, the PDF is always symmetric. In other words,
296 : ! the PDF at any positive value of Skx (for example, Skx = 2.5)
297 : ! will look like a mirror-image (reflection across the y-axis)
298 : ! of the PDF at a negative value of Skx of the same magnitude (in
299 : ! this example, Skx = -2.5). However, when zeta_x /= 0, the PDF
300 : ! loses this quality and is not symmetric.
301 : !
302 : ! When symmetry is desired at values of zeta_x besides zeta_x = 0,
303 : ! the solution is to turn zeta_x into a function of Skx. A basic
304 : ! example of a zeta_x skewness equation that produces a symmetric
305 : ! PDF for values of zeta_x other than 0 is:
306 : !
307 : ! zeta_x = | zeta_x_in, when Skx >= 0;
308 : ! | ( 1 / ( 1 + zeta_x_in ) ) - 1, when Skx < 0.
309 : !
310 : !
311 : ! Notes:
312 : !
313 : ! When F_x = 0 (which can only happen when Skx = 0), mu_x_1 = mu_x_2, and
314 : ! mixt_frac = ( zeta_x + 1 ) / ( zeta_x + 2 ). When these equations are
315 : ! substituted into the equations for sigma_x_1 and sigma_x_2, the result is
316 : ! sigma_x_1 = sigma_x_2 = sqrt( <x'^2> ). This means that the distribution
317 : ! becomes a single Gaussian when F_x = 0 (and Skx = 0). This happens
318 : ! regardless of the value of zeta_x.
319 : !
320 : ! The equations for the PDF component means and standard deviations can also
321 : ! be written as:
322 : !
323 : ! mu_x_1 = <x> + sqrt( F_x * ( ( 1 - mixt_frac ) / mixt_frac ) * <x'^2> )
324 : ! * sgn( <w'x'> );
325 : !
326 : ! mu_x_2 = <x> - sqrt( F_x * ( mixt_frac / ( 1 - mixt_frac ) ) * <x'^2> )
327 : ! * sgn( <w'x'> );
328 : !
329 : ! sigma_x_1 = sqrt( coef_sigma_x_1_sqd * <x'^2> ); and
330 : !
331 : ! sigma_x_2 = sqrt( coef_sigma_x_2_sqd * <x'^2> ); where
332 : !
333 : ! coef_sigma_x_1_sqd = ( ( zeta_x + 1 ) * ( 1 - F_x ) )
334 : ! / ( ( zeta_x + 2 ) * mixt_frac ); and
335 : !
336 : ! coef_sigma_x_2_sqd = ( 1 - F_x ) / ( ( zeta_x + 2 ) * ( 1 - mixt_frac ) ).
337 : !
338 : ! The above equations can be substituted into an equation for a variable that
339 : ! has been derived by integrating over the PDF. Many variables like this are
340 : ! used in parts of the predictive equation set. These substitutions allow
341 : ! some terms to solved implicitly or semi-implicitly in the predictive
342 : ! equations.
343 : !
344 : !
345 : ! Brian Griffin; September 2017.
346 : !
347 : !=============================================================================
348 0 : function calc_mixture_fraction( nz, Skx, F_x, zeta_x, sgn_wpxp ) &
349 0 : result( mixt_frac )
350 :
351 : ! Description:
352 : ! Calculates mixture fraction.
353 :
354 : ! References:
355 : ! Griffin and Larson (2018)
356 : !-----------------------------------------------------------------------
357 :
358 : use constants_clubb, only: &
359 : thirty_six, & ! Constant(s)
360 : eighteen, &
361 : twelve, &
362 : six, &
363 : four, &
364 : three, &
365 : two, &
366 : one, &
367 : zero, &
368 : fstderr
369 :
370 : use clubb_precision, only: &
371 : core_rknd ! Variable(s)
372 :
373 : implicit none
374 :
375 : integer, intent(in) :: &
376 : nz
377 :
378 : ! Input Variables
379 : real( kind = core_rknd ), dimension(nz), intent(in) :: &
380 : Skx, & ! Skewness of x [-]
381 : F_x, & ! Parameter for the spread of the PDF component means of x [-]
382 : zeta_x, & ! Parameter for the PDF component variances of x [-]
383 : sgn_wpxp ! Sign of the covariance of w and x [-]
384 :
385 : ! Return Variable
386 : real( kind = core_rknd ), dimension(nz) :: &
387 : mixt_frac ! Mixture fraction [-]
388 :
389 : ! Local Variable
390 : ! Flag that turns off when conditions aren't right for calculating mixt_frac
391 : logical, dimension(nz) :: &
392 0 : l_calc_mixt_frac
393 :
394 :
395 : ! Initialize l_calc_mixt_frac
396 0 : l_calc_mixt_frac = .true.
397 :
398 : ! Calculate mixture fraction, which is the weight of the 1st PDF component.
399 : ! The 2nd PDF component has a weight of 1 - mixt_frac.
400 0 : where ( F_x > zero )
401 :
402 : mixt_frac &
403 : = ( four * F_x**3 &
404 : + eighteen * F_x &
405 : * ( zeta_x + one ) * ( one - F_x ) / ( zeta_x + two ) &
406 : + six * F_x**2 * ( one - F_x ) / ( zeta_x + two ) &
407 : + Skx**2 &
408 : - Skx * sgn_wpxp * sqrt( four * F_x**3 &
409 : + twelve * F_x**2 * ( one - F_x ) &
410 : + thirty_six * F_x &
411 : * ( zeta_x + one ) * ( one - F_x )**2 &
412 : / ( zeta_x + two )**2 &
413 : + Skx**2 ) ) &
414 : / ( two * F_x * ( F_x - three )**2 + two * Skx**2 )
415 :
416 : elsewhere ! F_x = 0
417 :
418 : where ( abs( Skx ) > zero )
419 :
420 : l_calc_mixt_frac = .false.
421 :
422 : elsewhere ! Skx = 0
423 :
424 : mixt_frac = ( zeta_x + one ) / ( zeta_x + two )
425 :
426 : endwhere ! | Skx | > 0
427 :
428 : endwhere ! F_x > 0
429 :
430 :
431 0 : if ( any( .not. l_calc_mixt_frac ) ) then
432 0 : write(fstderr,*) "Mixture fraction cannot be calculated."
433 : write(fstderr,*) "The value of F_x must be greater than 0 when " &
434 0 : // "| Skx | > 0."
435 0 : error stop
436 : endif ! any( .not. l_valid_mixt_frac )
437 :
438 :
439 0 : return
440 :
441 0 : end function calc_mixture_fraction
442 :
443 : !=============================================================================
444 0 : subroutine calc_setter_var_params( nz, xm, xp2, Skx, sgn_wpxp, & ! In
445 0 : F_x, zeta_x, & ! In
446 0 : mu_x_1, mu_x_2, sigma_x_1, & ! Out
447 0 : sigma_x_2, mixt_frac, & ! Out
448 0 : coef_sigma_x_1_sqd, & ! Out
449 0 : coef_sigma_x_2_sqd ) ! Out
450 :
451 : ! Description:
452 : ! Calculates the PDF component means, the PDF component standard deviations,
453 : ! and the mixture fraction for the variable that sets the PDF.
454 :
455 : ! References:
456 : ! Griffin and Larson (2018)
457 : !-----------------------------------------------------------------------
458 :
459 : use constants_clubb, only: &
460 : two, & ! Variable(s)
461 : one
462 :
463 : use clubb_precision, only: &
464 : core_rknd ! Variable(s)
465 :
466 : implicit none
467 :
468 : integer, intent(in) :: &
469 : nz
470 :
471 : ! Input Variables
472 : real( kind = core_rknd ), dimension(nz), intent(in) :: &
473 : xm, & ! Mean of x (overall) [units vary]
474 : xp2, & ! Variance of x (overall) [(units vary)^2]
475 : Skx, & ! Skewness of x [-]
476 : sgn_wpxp, & ! Sign of the covariance of w and x (overall) [-]
477 : F_x, & ! Parameter for the spread of the PDF component means of x [-]
478 : zeta_x ! Parameter for the PDF component variances of x [-]
479 :
480 : ! Output Variables
481 : real( kind = core_rknd ), dimension(nz), intent(out) :: &
482 : mu_x_1, & ! Mean of x (1st PDF component) [units vary]
483 : mu_x_2, & ! Mean of x (2nd PDF component) [units vary]
484 : sigma_x_1, & ! Standard deviation of x (1st PDF component) [units vary]
485 : sigma_x_2, & ! Standard deviation of x (2nd PDF component) [units vary]
486 : mixt_frac ! Mixture fraction [-]
487 :
488 : real( kind = core_rknd ), dimension(nz), intent(out) :: &
489 : coef_sigma_x_1_sqd, & ! sigma_x_1^2 = coef_sigma_x_1_sqd * <x'^2> [-]
490 : coef_sigma_x_2_sqd ! sigma_x_2^2 = coef_sigma_x_2_sqd * <x'^2> [-]
491 :
492 :
493 : ! Calculate the mixture fraction.
494 0 : mixt_frac = calc_mixture_fraction( nz, Skx, F_x, zeta_x, sgn_wpxp )
495 :
496 : ! Calculate the mean of x in the 1st PDF component.
497 : mu_x_1 = xm + sqrt( F_x * ( ( one - mixt_frac ) / mixt_frac ) * xp2 ) &
498 0 : * sgn_wpxp
499 :
500 : ! Calculate the mean of x in the 2nd PDF component.
501 0 : mu_x_2 = xm - ( mixt_frac / ( one - mixt_frac ) ) * ( mu_x_1 - xm )
502 :
503 : ! Calculate the standard deviation of x in the 1st PDF component.
504 : ! sigma_x_1 = sqrt( ( ( zeta_x + 1 ) * ( 1 - F_x ) )
505 : ! / ( ( zeta_x + 2 ) * mixt_frac ) * <x'^2> )
506 : coef_sigma_x_1_sqd = ( ( zeta_x + one ) * ( one - F_x ) ) &
507 0 : / ( ( zeta_x + two ) * mixt_frac )
508 :
509 0 : sigma_x_1 = sqrt( coef_sigma_x_1_sqd * xp2 )
510 :
511 : ! Calculate the standard deviation of x in the 2nd PDF component.
512 : ! sigma_x_2 = sqrt( ( mixt_frac * sigma_x_1^2 )
513 : ! / ( ( 1 - mixt_frac ) * ( 1 + zeta_x ) ) )
514 : ! = sqrt( ( 1 - F_x )
515 : ! / ( ( zeta_x + 2 ) * ( 1 - mixt_frac ) ) * <x'^2> )
516 : coef_sigma_x_2_sqd = ( one - F_x ) &
517 0 : / ( ( zeta_x + two ) * ( one - mixt_frac ) )
518 :
519 0 : sigma_x_2 = sqrt( coef_sigma_x_2_sqd * xp2 )
520 :
521 :
522 0 : return
523 :
524 : end subroutine calc_setter_var_params
525 :
526 : !=============================================================================
527 : !
528 : ! DESCRIPTION OF THE METHOD FOR EACH RESPONDING VARIABLE
529 : ! ======================================================
530 : !
531 : ! In order to find equations for the four PDF parameters for each responding
532 : ! variable, which are mu_x_1, mu_x_2, sigma_x_1, and sigma_x_2 (where x stands
533 : ! for a responding variable here), four equations are needed. These four
534 : ! equations are the equations for <x>, <x'^2>, and <x'^3> as found by
535 : ! integrating over the PDF. Additionally, one more equation, which involves
536 : ! a tunable parameter F_x, and which is used to help control the spread of the
537 : ! PDF component means, is used in this equation set. The four equations are:
538 : !
539 : ! <x> = mixt_frac * mu_x_1 + ( 1 - mixt_frac ) * mu_x_2;
540 : !
541 : ! <x'^2> = mixt_frac * ( ( mu_x_1 - <x> )^2 + sigma_x_1^2 )
542 : ! + ( 1 - mixt_frac ) * ( ( mu_x_2 - <x> )^2 + sigma_x_2^2 );
543 : !
544 : ! <x'^3> = mixt_frac * ( mu_x_1 - <x> )
545 : ! * ( ( mu_x_1 - <x> )^2 + 3 * sigma_x_1^2 )
546 : ! + ( 1 - mixt_frac ) * ( mu_x_2 - <x> )
547 : ! * ( ( mu_x_2 - <x> )^2 + 3 * sigma_x_2^2 ); and
548 : !
549 : ! mu_x_1 - <x> = sqrt(F_x) * ( sqrt( 1 - mixt_frac ) / sqrt( mixt_frac ) )
550 : ! * sqrt( <x'^2> ) * sgn( <w'x'> );
551 : !
552 : ! where 0 <= F_x <= 1, and where sgn( <w'x'> ) is given by:
553 : !
554 : ! sgn( <w'x'> ) = | 1, when <w'x'> >= 0;
555 : ! | -1, when <w'x'> < 0.
556 : !
557 : ! The resulting equations for the four PDF parameters are:
558 : !
559 : ! mu_x_1 = <x> + sqrt( F_x * ( ( 1 - mixt_frac ) / mixt_frac ) * <x'^2> )
560 : ! * sgn( <w'x'> );
561 : !
562 : ! mu_x_2 = <x> - ( mixt_frac / ( 1 - mixt_frac ) ) * ( mu_x_1 - <x> );
563 : !
564 : ! sigma_x_1^2
565 : ! = ( ( sqrt( mixt_frac * ( 1 - mixt_frac ) ) * Skx * sgn( <w'x'> )
566 : ! - ( 1 + mixt_frac ) * F_x^1.5 + 3 * mixt_frac * sqrt( F_x ) )
567 : ! / ( 3 * mixt_frac * sqrt( F_x ) ) )
568 : ! * <x'^2>; and
569 : !
570 : ! sigma_x_2^2 = ( ( 1 - F_x ) / ( 1 - mixt_frac ) ) * <x'^2>
571 : ! - ( mixt_frac / ( 1 - mixt_frac ) ) * sigma_x_1^2;
572 : !
573 : ! where Skx is the skewness of x, and Skx = <x'^3> / <x'^2>^(3/2).
574 : !
575 : !
576 : ! Special case:
577 : !
578 : ! When Skx = 0 and F_x = 0, the equations for sigma_x_1^2 and sigma_x_2^2 are
579 : ! both undefined. The equations for sigma_x_1^2 and sigma_x_2^2 in this
580 : ! scenario can be derived by using the above equations for sigma_x_1^2 and
581 : ! sigma_x_2^2 and then setting Skx = 0. The resulting equation for
582 : ! sigma_x_1^2 becomes:
583 : !
584 : ! sigma_x_1^2
585 : ! = ( ( - ( 1 + mixt_frac ) * F_x^1.5 + 3 * mixt_frac * sqrt( F_x ) )
586 : ! / ( 3 * mixt_frac * sqrt( F_x ) ) )
587 : ! * <x'^2>.
588 : !
589 : ! All of the terms in the numerator and denominator contain a sqrt( F_x ),
590 : ! so this equation can be rewritten as:
591 : !
592 : ! sigma_x_1^2
593 : ! = ( ( - ( 1 + mixt_frac ) * F_x + 3 * mixt_frac ) / ( 3 * mixt_frac ) )
594 : ! * <x'^2>.
595 : !
596 : ! Now setting F_x = 0, the equation becomes:
597 : !
598 : ! sigma_x_1^2 = ( ( 3 * mixt_frac ) / ( 3 * mixt_frac ) ) * <x'^2>;
599 : !
600 : ! which can be rewritten as:
601 : !
602 : ! sigma_x_1^2 = <x'^2>.
603 : !
604 : ! Substituting the equation for sigma_x_1^2 into the equation for sigma_x_2^2,
605 : ! and also setting F_x = 0, the equation for sigma_x_2^2 becomes:
606 : !
607 : ! sigma_x_2^2 = ( 1 / ( 1 - mixt_frac ) ) * <x'^2>
608 : ! - ( mixt_frac / ( 1 - mixt_frac ) ) * <x'^2>;
609 : !
610 : ! which can be rewritten as:
611 : !
612 : ! sigma_x_2^2
613 : ! = ( ( 1 / ( 1 - mixt_frac ) ) - ( mixt_frac / ( 1 - mixt_frac ) ) )
614 : ! * <x'^2>.
615 : !
616 : ! This equation becomes:
617 : !
618 : ! sigma_x_2^2 = ( ( 1 - mixt_frac ) / ( 1 - mixt_frac ) ) * <x'^2>;
619 : !
620 : ! which can be rewritten as:
621 : !
622 : ! sigma_x_2^2 = <x'^2>.
623 : !
624 : ! When F_x = 0, Skx must have a value of 0 in order for the PDF to function
625 : ! correctly. When F_x = 0, mu_x_1 = mu_x_2. When the two PDF component means
626 : ! are equal to each other (and to the overall mean, <x>), the only value of
627 : ! Skx that can be represented is a value of 0. The equations that place
628 : ! limits on F_x for a responding variable (below) calculate the minimum
629 : ! allowable value of F_x to be greater than 0 when | Skx | > 0.
630 : !
631 : ! The value of F_x should be set as a function of Skx. The value F_x should
632 : ! go toward 0 as | Skx | (or Skx^2) goes toward 0. The value of F_x should
633 : ! go toward 1 as | Skx | (or Skx^2) goes to infinity. However, the value of
634 : ! F_x must also be between the minimum and maximum allowable values of F_x for
635 : ! a responding variable (below).
636 : !
637 : !
638 : ! Tunable parameter:
639 : !
640 : ! F_x: This parameter controls the spread of the PDF component means. The
641 : ! range of this parameter is 0 <= F_x <= 1. When F_x = 0, the two PDF
642 : ! component means (mu_x_1 and mu_x_2) are equal to each other (and Skx
643 : ! must equal 0). All of the variance of x is accounted for by the PDF
644 : ! component standard deviations (sigma_x_1 and sigma_x_2). When
645 : ! F_x = 1, mu_x_1 and mu_x_2 are spread as far apart as they can be.
646 : ! Both PDF component standard deviations (sigma_x_1 and sigma_x_2) are
647 : ! equal to 0, and all of the variance of x is accounted for by the
648 : ! spread of the PDF component means.
649 : !
650 : !
651 : ! Limits on F_x:
652 : !
653 : ! Since the PDF parameters for this variable need to work with the mixture
654 : ! fraction that has been provided by the setting variable, the method does
655 : ! not work for all values of F_x and Skx. However, the limits of Skx and F_x
656 : ! can always be calculated. The limits are based on keeping the values of
657 : ! sigma_x_1 and sigma_x_2 greater than or equal to 0. The equation for
658 : ! keeping the value of sigma_x_1 greater than or equal to 0 is:
659 : !
660 : ! - ( 1 + mixt_frac ) * sqrt( F_x )^3 + 3 * mixt_frac * sqrt( F_x )
661 : ! + sqrt( mixt_frac * ( 1 - mixt_frac ) ) * Skx * sgn( <w'x'> ) >= 0.
662 : !
663 : ! The roots of sqrt( F_x ) can be solved by an equation of the form:
664 : !
665 : ! A * sqrt( F_x )^3 + B * sqrt( F_x )^2 + C * sqrt( F_x ) + D = 0;
666 : !
667 : ! where:
668 : !
669 : ! A = - ( 1 + mixt_frac );
670 : ! B = 0;
671 : ! C = 3 * mixt_frac; and
672 : ! D = sqrt( mixt_frac * ( 1 - mixt_frac ) ) * Skx * sgn( <w'x'> ).
673 : !
674 : ! The equation for keeping the value of sigma_x_2 greater than or equal to 0
675 : ! is:
676 : !
677 : ! - ( 2 - mixt_frac ) * sqrt( F_x )^3 + 3 * ( 1 - mixt_frac ) * sqrt( F_x )
678 : ! - sqrt( mixt_frac * ( 1 - mixt_frac ) ) * Skx * sgn( <w'x'> ) >= 0.
679 : !
680 : ! The roots of sqrt( F_x ) can be solved by an equation of the form:
681 : !
682 : ! A * sqrt( F_x )^3 + B * sqrt( F_x )^2 + C * sqrt( F_x ) + D = 0;
683 : !
684 : ! where:
685 : !
686 : ! A = - ( 2 - mixt_frac );
687 : ! B = 0;
688 : ! C = 3 * ( 1 - mixt_frac ); and
689 : ! D = - sqrt( mixt_frac * ( 1 - mixt_frac ) ) * Skx * sgn( <w'x'> ).
690 : !
691 : ! After careful analysis of the above equations, the following properties
692 : ! emerge:
693 : !
694 : ! When Skx * sgn( <w'x'> ) >= 0,
695 : ! Skx^2 < 4 * ( 1 - mixt_frac )^2 / ( mixt_frac * ( 2 - mixt_frac ) )
696 : ! is required; and
697 : ! when Skx * sgn( <w'x'> ) < 0,
698 : ! Skx^2 < 4 * mixt_frac^2 / ( 1 - mixt_frac^2 ) is required.
699 : !
700 : ! Whenever Skx^2 exceeds these limits, Skx must be limited (preserving its
701 : ! sign) in order to have any value of F_x that will work in the equation set.
702 : !
703 : ! When Skx is found to be within the above limits (or after it has been
704 : ! limited to fall within its limits), the range of valid values of F_x can be
705 : ! found according to the following:
706 : !
707 : ! When Skx * sgn( <w'x'> ) >= 0:
708 : !
709 : ! When 4 * mixt_frac^2 / ( 1 - mixt_frac^2 ) < Skx^2
710 : ! < 4 * ( 1 - mixt_frac )^2 / ( mixt_frac * ( 2 - mixt_frac ) ):
711 : !
712 : ! Minimum sqrt( F_x ): 2nd root (middle-valued root; also smallest
713 : ! positive) of the second equation (sigma_x_2
714 : ! based).
715 : !
716 : ! Maximum sqrt( F_x ): Minimum of the largest root of the second
717 : ! equation (sigma_x_2 based) and the only* root
718 : ! of the first equation (sigma_x_1 based).
719 : !
720 : ! When Skx^2 <= 4 * mixt_frac^2 / ( 1 - mixt_frac^2 ):
721 : !
722 : ! Minimum sqrt( F_x ): 2nd root (middle-valued root; also smallest
723 : ! positive) of the second equation (sigma_x_2
724 : ! based).
725 : !
726 : ! Maximum sqrt( F_x ): Minimum of the largest root of the second
727 : ! equation (sigma_x_2 based) and the largest
728 : ! root of the first equation (sigma_x_1 based).
729 : !
730 : ! When Skx * sgn( <w'x'> ) < 0:
731 : !
732 : ! When 4 * ( 1 - mixt_frac )^2 / ( mixt_frac * ( 2 - mixt_frac ) )
733 : ! < Skx^2 < 4 * mixt_frac^2 / ( 1 - mixt_frac^2 ):
734 : !
735 : ! Minimum sqrt( F_x ): 2nd root (middle-valued root; also smallest
736 : ! positive) of the first equation (sigma_x_1
737 : ! based).
738 : !
739 : ! Maximum sqrt( F_x ): Minimum of the largest root of the first
740 : ! equation (sigma_x_1 based) and the only* root
741 : ! of the second equation (sigma_x_2 based).
742 : !
743 : ! When Skx^2
744 : ! <= 4 * ( 1 - mixt_frac )^2 / ( mixt_frac * ( 2 - mixt_frac ) ):
745 : !
746 : ! Minimum sqrt( F_x ): 2nd root (middle-valued root; also smallest
747 : ! positive) of the first equation (sigma_x_1
748 : ! based).
749 : !
750 : ! Maximum sqrt( F_x ): Minimum of the largest root of the first
751 : ! equation (sigma_x_1 based) and the largest
752 : ! root of the second equation (sigma_x_2
753 : ! based).
754 : !
755 : ! Here, "only* root" means the the only root that isn't a complex root.
756 : !
757 : ! The value of sqrt( F_x ) is also limited with a minimum of 0 and a maximum
758 : ! of 1. The minimum and maximum allowable values of F_x are found by taking
759 : ! the square of the minimum and maximum allowable values of sqrt( F_x ),
760 : ! respectively.
761 : !
762 : !
763 : ! Notes:
764 : !
765 : ! When F_x = 0 (which can only happen when Skx = 0), mu_x_1 = mu_x_2, and
766 : ! sigma_x_1 = sigma_x_2 = sqrt( <x'^2> ). This means that the distribution
767 : ! becomes a single Gaussian when F_x = 0 (and Skx = 0).
768 : !
769 : ! The equations for the PDF component means and standard deviations can also
770 : ! be written as:
771 : !
772 : ! mu_x_1 = <x> + sqrt( F_x * ( ( 1 - mixt_frac ) / mixt_frac ) * <x'^2> )
773 : ! * sgn( <w'x'> );
774 : !
775 : ! mu_x_2 = <x> - sqrt( F_x * ( mixt_frac / ( 1 - mixt_frac ) ) * <x'^2> )
776 : ! * sgn( <w'x'> );
777 : !
778 : ! sigma_x_1 = sqrt( coef_sigma_x_1_sqd * <x'^2> ); and
779 : !
780 : ! sigma_x_2 = sqrt( coef_sigma_x_2_sqd * <x'^2> ); where
781 : !
782 : ! coef_sigma_x_1_sqd
783 : ! = ( sqrt( mixt_frac * ( 1 - mixt_frac ) ) * Skx * sgn( <w'x'> )
784 : ! - ( 1 + mixt_frac ) * F_x^1.5 + 3 * mixt_frac * sqrt( F_x ) )
785 : ! / ( 3 * mixt_frac * sqrt( F_x ) )
786 : ! = sqrt( mixt_frac * ( 1 - mixt_frac ) ) * Skx * sgn( <w'x'> )
787 : ! / ( 3 * mixt_frac * sqrt( F_x ) )
788 : ! - ( 1 + mixt_frac ) * F_x / ( 3 * mixt_frac )
789 : ! + 1; and
790 : !
791 : ! coef_sigma_x_2_sqd
792 : ! = ( 1 - F_x ) / ( 1 - mixt_frac )
793 : ! - mixt_frac / ( 1 - mixt_frac )
794 : ! * ( sqrt( mixt_frac * ( 1 - mixt_frac ) ) * Skx * sgn( <w'x'> )
795 : ! / ( 3 * mixt_frac * sqrt( F_x ) )
796 : ! - ( 1 + mixt_frac ) * F_x / ( 3 * mixt_frac )
797 : ! + 1 )
798 : ! = ( ( 1 - F_x ) - mixt_frac * coef_sigma_x_1_sqd ) / ( 1 - mixt_frac ).
799 : !
800 : ! The above equations can be substituted into an equation for a variable that
801 : ! has been derived by integrating over the PDF. Many variables like this are
802 : ! used in parts of the predictive equation set. These substitutions allow
803 : ! some terms to solved implicitly or semi-implicitly in the predictive
804 : ! equations.
805 : !
806 : !
807 : ! Brian Griffin; September 2017.
808 : !
809 : !=============================================================================
810 0 : subroutine calc_responder_params( nz, xm, xp2, Skx, sgn_wpxp, & ! In
811 0 : F_x, mixt_frac, & ! In
812 0 : mu_x_1, mu_x_2, & ! Out
813 0 : sigma_x_1_sqd, sigma_x_2_sqd, & ! Out
814 0 : coef_sigma_x_1_sqd, & ! Out
815 0 : coef_sigma_x_2_sqd ) ! Out
816 :
817 : ! Description:
818 : ! Calculates the PDF component means and the PDF component standard
819 : ! deviations for a responding variable (a variable that is not used to set
820 : ! the mixture fraction).
821 :
822 : ! References:
823 : ! Griffin and Larson (2018)
824 : !-----------------------------------------------------------------------
825 :
826 : use constants_clubb, only: &
827 : three, & ! Variable(s)
828 : one, &
829 : zero
830 :
831 : use clubb_precision, only: &
832 : core_rknd ! Variable(s)
833 :
834 : implicit none
835 :
836 : integer, intent(in) :: &
837 : nz
838 :
839 : ! Input Variables
840 : real( kind = core_rknd ), dimension(nz), intent(in) :: &
841 : xm, & ! Mean of x (overall) [units vary]
842 : xp2, & ! Variance of x (overall) [(units vary)^2]
843 : Skx, & ! Skewness of x [-]
844 : sgn_wpxp, & ! Sign of the covariance of w and x (overall) [-]
845 : F_x, & ! Parameter for the spread of the PDF component means of x [-]
846 : mixt_frac ! Mixture fraction [-]
847 :
848 : ! Output Variables
849 : real( kind = core_rknd ), dimension(nz), intent(out) :: &
850 : mu_x_1, & ! Mean of x (1st PDF component) [units vary]
851 : mu_x_2, & ! Mean of x (2nd PDF component) [units vary]
852 : sigma_x_1_sqd, & ! Variance of x (1st PDF component) [(units vary)^2]
853 : sigma_x_2_sqd ! Variance of x (2nd PDF component) [(units vary)^2]
854 :
855 : real( kind = core_rknd ), dimension(nz), intent(out) :: &
856 : coef_sigma_x_1_sqd, & ! sigma_x_1^2 = coef_sigma_x_1_sqd * <x'^2> [-]
857 : coef_sigma_x_2_sqd ! sigma_x_2^2 = coef_sigma_x_2_sqd * <x'^2> [-]
858 :
859 :
860 0 : where ( F_x > zero )
861 :
862 : ! Calculate the mean of x in the 1st PDF component.
863 : mu_x_1 = xm + sqrt( F_x * ( ( one - mixt_frac ) / mixt_frac ) * xp2 ) &
864 : * sgn_wpxp
865 :
866 : ! Calculate the mean of x in the 2nd PDF component.
867 : mu_x_2 = xm - ( mixt_frac / ( one - mixt_frac ) ) * ( mu_x_1 - xm )
868 :
869 : ! Calculate the variance of x in the 1st PDF component.
870 : ! sigma_x_1^2
871 : ! = ( ( sqrt( mixt_frac * ( 1 - mixt_frac ) ) * Skx * sgn( <w'x'> )
872 : ! - ( 1 + mixt_frac ) * F_x^1.5 + 3 * mixt_frac * sqrt( F_x ) )
873 : ! / ( 3 * mixt_frac * sqrt( F_x ) ) ) * <x'^2>
874 : ! = ( sqrt( mixt_frac * ( 1 - mixt_frac ) ) * Skx * sgn( <w'x'> )
875 : ! / ( 3 * mixt_frac * sqrt( F_x ) )
876 : ! - ( 1 + mixt_frac ) * F_x / ( 3 * mixt_frac )
877 : ! + 1 ) * <x'^2>
878 : coef_sigma_x_1_sqd &
879 : = sqrt( mixt_frac * ( one - mixt_frac ) ) * Skx * sgn_wpxp &
880 : / ( three * mixt_frac * sqrt( F_x ) ) &
881 : - ( one + mixt_frac ) * F_x / ( three * mixt_frac ) &
882 : + one
883 :
884 : sigma_x_1_sqd = coef_sigma_x_1_sqd * xp2
885 :
886 : ! Calculate the variance of x in the 2nd PDF component.
887 : ! sigma_x_2^2
888 : ! = ( ( 1 - F_x ) / ( 1 - mixt_frac )
889 : ! - mixt_frac / ( 1 - mixt_frac )
890 : ! * ( sqrt( mixt_frac * ( 1 - mixt_frac ) ) * Skx * sgn( <w'x'> )
891 : ! / ( 3 * mixt_frac * sqrt( F_x ) )
892 : ! - ( 1 + mixt_frac ) * F_x / ( 3 * mixt_frac )
893 : ! + 1 ) ) * <x'^2>
894 : ! = ( ( ( 1 - F_x ) - mixt_frac * coef_sigma_x_1_sqd )
895 : ! / ( 1 - mixt_frac ) ) * <x'^2>
896 : coef_sigma_x_2_sqd &
897 : = ( ( one - F_x ) - mixt_frac * coef_sigma_x_1_sqd ) &
898 : / ( one - mixt_frac )
899 :
900 : sigma_x_2_sqd = coef_sigma_x_2_sqd * xp2
901 :
902 : elsewhere ! F_x = 0
903 :
904 : ! When F_x has a value of 0, the PDF becomes a single Gaussian. This
905 : ! only works when Skx = 0. However, when Skx /= 0, the value of min_F_x
906 : ! is greater than 0, preventing a problem where F_x = 0 but | Skx | > 0.
907 : mu_x_1 = xm
908 : mu_x_2 = xm
909 : sigma_x_1_sqd = xp2
910 : sigma_x_2_sqd = xp2
911 : coef_sigma_x_1_sqd = one
912 : coef_sigma_x_2_sqd = one
913 :
914 : endwhere ! F_x > 0
915 :
916 :
917 0 : return
918 :
919 : end subroutine calc_responder_params
920 :
921 : !=============================================================================
922 0 : subroutine calc_limits_F_x_responder( nz, mixt_frac, Skx, sgn_wpxp, & ! In
923 0 : max_Skx2_pos_Skx_sgn_wpxp, & ! In
924 0 : max_Skx2_neg_Skx_sgn_wpxp, & ! In
925 0 : min_F_x, max_F_x ) ! Out
926 :
927 : ! Description:
928 : ! Calculates the minimum and maximum allowable values for F_x for a
929 : ! responding variable.
930 :
931 : ! References:
932 : !-----------------------------------------------------------------------
933 :
934 : use constants_clubb, only: &
935 : three, & ! Variable(s)
936 : two, &
937 : one, &
938 : zero
939 :
940 : use calc_roots, only: &
941 : cubic_solve ! Procedure(s)
942 :
943 : use clubb_precision, only: &
944 : core_rknd ! Variable(s)
945 :
946 : implicit none
947 :
948 : integer, intent(in) :: &
949 : nz
950 :
951 : ! Input Variables
952 : real( kind = core_rknd ), dimension(nz), intent(in) :: &
953 : mixt_frac, & ! Mixture fraction [-]
954 : Skx, & ! Skewness of x [-]
955 : sgn_wpxp ! Sign of covariance of w and x [-]
956 :
957 : real( kind = core_rknd ), dimension(nz), intent(in) :: &
958 : max_Skx2_pos_Skx_sgn_wpxp, & ! Maximum Skx^2 when Skx*sgn(<w'x'>) >= 0 [-]
959 : max_Skx2_neg_Skx_sgn_wpxp ! Maximum Skx^2 when Skx*sgn(<w'x'>) < 0 [-]
960 :
961 : ! Output Variables
962 : real( kind = core_rknd ), dimension(nz), intent(out) :: &
963 : min_F_x, & ! Minimum allowable value of F_x [-]
964 : max_F_x ! Maximum allowable value of F_x [-]
965 :
966 : ! Local Variables
967 : real( kind = core_rknd ), dimension(nz) :: &
968 0 : coef_A_1, & ! Coef. A in Ax^3 + Bx^2 + Cx + D = 0 (1st PDF comp. lim.) [-]
969 0 : coef_B_1, & ! Coef. B in Ax^3 + Bx^2 + Cx + D = 0 (1st PDF comp. lim.) [-]
970 0 : coef_C_1, & ! Coef. C in Ax^3 + Bx^2 + Cx + D = 0 (1st PDF comp. lim.) [-]
971 0 : coef_D_1, & ! Coef. D in Ax^3 + Bx^2 + Cx + D = 0 (1st PDF comp. lim.) [-]
972 0 : coef_A_2, & ! Coef. A in Ax^3 + Bx^2 + Cx + D = 0 (2nd PDF comp. lim.) [-]
973 0 : coef_B_2, & ! Coef. B in Ax^3 + Bx^2 + Cx + D = 0 (2nd PDF comp. lim.) [-]
974 0 : coef_C_2, & ! Coef. C in Ax^3 + Bx^2 + Cx + D = 0 (2nd PDF comp. lim.) [-]
975 0 : coef_D_2 ! Coef. D in Ax^3 + Bx^2 + Cx + D = 0 (2nd PDF comp. lim.) [-]
976 :
977 : complex( kind = core_rknd ), dimension(nz,3) :: &
978 0 : sqrt_F_x_roots_1, & ! Roots of sqrt(F_x) for the sigma_x_1 term [-]
979 0 : sqrt_F_x_roots_2 ! Roots of sqrt(F_x) for the sigma_x_2 term [-]
980 :
981 : real( kind = core_rknd ), dimension(nz,3) :: &
982 0 : sqrt_F_x_roots_1_sorted, & ! Sorted roots of sqrt(F_x): sigma_x_1 term [-]
983 0 : sqrt_F_x_roots_2_sorted ! Sorted roots of sqrt(F_x): sigma_x_2 term [-]
984 :
985 : real( kind = core_rknd ), dimension(nz) :: &
986 0 : min_sqrt_F_x, & ! Minimum allowable value of sqrt(F_x) [-]
987 0 : max_sqrt_F_x ! Maximum allowable value of sqrt(F_x) [-]
988 :
989 :
990 : ! Set up the coefficients in the equation for the limit of sqrt(F_x) based
991 : ! on the 1st PDF component standard deviation (sigma_x_1) being greater than
992 : ! or equal to 0. This equation has the form:
993 : ! A * sqrt(F_x)^3 + B * sqrt(F_x)^2 + C * sqrt(F_x) + D = 0.
994 0 : coef_A_1 = -( one + mixt_frac )
995 0 : coef_B_1 = zero
996 0 : coef_C_1 = three * mixt_frac
997 0 : coef_D_1 = sqrt( mixt_frac * ( one - mixt_frac ) ) * Skx * sgn_wpxp
998 :
999 : ! Solve for the roots (values of sqrt(F_x)) that satisfy the above equation.
1000 : sqrt_F_x_roots_1 &
1001 0 : = cubic_solve( nz, coef_A_1, coef_B_1, coef_C_1, coef_D_1 )
1002 :
1003 : ! Sort the values of the roots (values of sqrt(F_x)) from smallest to
1004 : ! largest. Ignore any complex component of the roots. The code below that
1005 : ! uses sqrt_F_x_roots_1_sorted already factors the appropriate roots to use
1006 : ! into account.
1007 : sqrt_F_x_roots_1_sorted &
1008 0 : = sort_roots( nz, real( sqrt_F_x_roots_1, kind = core_rknd ) )
1009 :
1010 : ! Set up the coefficients in the equation for the limit of sqrt(F_x) based
1011 : ! on the 2nd PDF component standard deviation (sigma_x_2) being greater than
1012 : ! or equal to 0. This equation has the form:
1013 : ! A * sqrt(F_x)^3 + B * sqrt(F_x)^2 + C * sqrt(F_x) + D = 0.
1014 0 : coef_A_2 = -( two - mixt_frac )
1015 0 : coef_B_2 = zero
1016 0 : coef_C_2 = three * ( one - mixt_frac )
1017 0 : coef_D_2 = -sqrt( mixt_frac * ( one - mixt_frac ) ) * Skx * sgn_wpxp
1018 :
1019 : ! Solve for the roots (values of sqrt(F_x)) that satisfy the above equation.
1020 : sqrt_F_x_roots_2 &
1021 0 : = cubic_solve( nz, coef_A_2, coef_B_2, coef_C_2, coef_D_2 )
1022 :
1023 : ! Sort the values of the roots (values of sqrt(F_x)) from smallest to
1024 : ! largest. Ignore any complex component of the roots. The code below that
1025 : ! uses sqrt_F_x_roots_2_sorted already factors the appropriate roots to use
1026 : ! into account.
1027 : sqrt_F_x_roots_2_sorted &
1028 0 : = sort_roots( nz, real( sqrt_F_x_roots_2, kind = core_rknd ) )
1029 :
1030 :
1031 : ! Find the minimum and maximum allowable values of sqrt(F_x) based on Skx
1032 : ! and sgn( <w'x'> ).
1033 0 : where ( Skx * sgn_wpxp >= zero )
1034 :
1035 : where ( Skx**2 > max_Skx2_neg_Skx_sgn_wpxp )
1036 :
1037 : min_sqrt_F_x = sqrt_F_x_roots_2_sorted(:,2)
1038 : max_sqrt_F_x = min( real( sqrt_F_x_roots_1(:,1), kind = core_rknd ), &
1039 : sqrt_F_x_roots_2_sorted(:,3) )
1040 :
1041 : elsewhere ! Skx^2 <= max_Skx2_neg_Skx_sgn_wpxp
1042 :
1043 : min_sqrt_F_x = sqrt_F_x_roots_2_sorted(:,2)
1044 : max_sqrt_F_x = min( sqrt_F_x_roots_1_sorted(:,3), &
1045 : sqrt_F_x_roots_2_sorted(:,3) )
1046 :
1047 : endwhere ! Skx**2 > max_Skx2_neg_Skx_sgn_wpxp
1048 :
1049 : elsewhere ! Skx * sgn( <w'x'> ) < 0
1050 :
1051 : where ( Skx**2 > max_Skx2_pos_Skx_sgn_wpxp )
1052 :
1053 : min_sqrt_F_x = sqrt_F_x_roots_1_sorted(:,2)
1054 : max_sqrt_F_x = min( real( sqrt_F_x_roots_2(:,1), kind = core_rknd ), &
1055 : sqrt_F_x_roots_1_sorted(:,3) )
1056 :
1057 : elsewhere ! Skx^2 <= max_Skx2_pos_Skx_sgn_wpxp
1058 :
1059 : min_sqrt_F_x = sqrt_F_x_roots_1_sorted(:,2)
1060 : max_sqrt_F_x = min( sqrt_F_x_roots_1_sorted(:,3), &
1061 : sqrt_F_x_roots_2_sorted(:,3) )
1062 :
1063 : endwhere ! Skx**2 > max_Skx2_pos_Skx_sgn_wpxp
1064 :
1065 : endwhere ! Skx * sgn( <w'x'> ) >= 0
1066 :
1067 :
1068 : ! The minimum and maximum are also limited by 0 and 1, respectively.
1069 0 : min_sqrt_F_x = max( min_sqrt_F_x, zero )
1070 0 : max_sqrt_F_x = min( max_sqrt_F_x, one )
1071 :
1072 : ! The minimum and maximum allowable values for F_x are the squares of the
1073 : ! minimum and maximum allowable values for sqrt(F_x).
1074 0 : min_F_x = min_sqrt_F_x**2
1075 0 : max_F_x = max_sqrt_F_x**2
1076 :
1077 :
1078 0 : return
1079 :
1080 : end subroutine calc_limits_F_x_responder
1081 :
1082 : !=============================================================================
1083 0 : function sort_roots( nz, roots ) &
1084 0 : result ( roots_sorted )
1085 :
1086 : ! Description:
1087 : ! Sorts roots from smallest to largest.
1088 :
1089 : ! References:
1090 : !-----------------------------------------------------------------------
1091 :
1092 : use clubb_precision, only: &
1093 : core_rknd ! Variable(s)
1094 :
1095 : implicit none
1096 :
1097 : integer, intent(in) :: &
1098 : nz
1099 :
1100 : ! Input Variable
1101 : real( kind = core_rknd ), dimension(nz,3), intent(in) :: &
1102 : roots ! Roots [-]
1103 :
1104 : ! Return Variable
1105 : real( kind = core_rknd ), dimension(nz,3) :: &
1106 : roots_sorted ! Roots sorted from smallest to largest [-]
1107 :
1108 :
1109 0 : where ( roots(:,1) <= roots(:,2) .and. roots(:,1) <= roots(:,3) )
1110 :
1111 : ! The value of roots(1) is the smallest root.
1112 : roots_sorted(:,1) = roots(:,1)
1113 :
1114 : where ( roots(:,2) <= roots(:,3) )
1115 :
1116 : ! The value of roots(2) is the middle-valued root and the value of
1117 : ! roots(3) is the largest root.
1118 : roots_sorted(:,2) = roots(:,2)
1119 : roots_sorted(:,3) = roots(:,3)
1120 :
1121 : elsewhere ! roots(3) < roots(2)
1122 :
1123 : ! The value of roots(3) is the middle-valued root and the value of
1124 : ! roots(2) is the largest root.
1125 : roots_sorted(:,2) = roots(:,3)
1126 : roots_sorted(:,3) = roots(:,2)
1127 :
1128 : endwhere ! roots(2) <= roots(3)
1129 :
1130 : elsewhere ( roots(:,2) < roots(:,1) .and. roots(:,2) <= roots(:,3) )
1131 :
1132 : ! The value of roots(2) is the smallest root.
1133 : roots_sorted(:,1) = roots(:,2)
1134 :
1135 : where ( roots(:,1) <= roots(:,3) )
1136 :
1137 : ! The value of roots(1) is the middle-valued root and the value of
1138 : ! roots(3) is the largest root.
1139 : roots_sorted(:,2) = roots(:,1)
1140 : roots_sorted(:,3) = roots(:,3)
1141 :
1142 : elsewhere ! roots(3) < roots(1)
1143 :
1144 : ! The value of roots(3) is the middle-valued root and the value of
1145 : ! roots(1) is the largest root.
1146 : roots_sorted(:,2) = roots(:,3)
1147 : roots_sorted(:,3) = roots(:,1)
1148 :
1149 : endwhere ! roots(1) <= roots(3)
1150 :
1151 : elsewhere ! roots(3) < roots(1) .and. roots(3) < roots(2)
1152 :
1153 : ! The value of roots(3) is the smallest root.
1154 : roots_sorted(:,1) = roots(:,3)
1155 :
1156 : where ( roots(:,1) <= roots(:,2) )
1157 :
1158 : ! The value of roots(1) is the middle-valued root and the value of
1159 : ! roots(2) is the largest root.
1160 : roots_sorted(:,2) = roots(:,1)
1161 : roots_sorted(:,3) = roots(:,2)
1162 :
1163 : elsewhere ! roots(2) < roots(1)
1164 :
1165 : ! The value of roots(2) is the middle-valued root and the value of
1166 : ! roots(1) is the largest root.
1167 : roots_sorted(:,2) = roots(:,2)
1168 : roots_sorted(:,3) = roots(:,1)
1169 :
1170 : endwhere ! roots(1) <= roots(2)
1171 :
1172 : endwhere ! roots(1) <= roots(2) .and. roots(1) <= roots(3)
1173 :
1174 :
1175 0 : return
1176 :
1177 0 : end function sort_roots
1178 :
1179 : !=============================================================================
1180 0 : function calc_coef_wp4_implicit( nz, mixt_frac, F_w, &
1181 0 : coef_sigma_w_1_sqd, &
1182 0 : coef_sigma_w_2_sqd ) &
1183 0 : result( coef_wp4_implicit )
1184 :
1185 : ! Description:
1186 : ! The predictive equation for <w'^3> contains a turbulent advection term of
1187 : ! the form:
1188 : !
1189 : ! - ( 1 / rho_ds ) * d ( rho_ds * <w'^4> ) / dz;
1190 : !
1191 : ! where z is height, rho_ds is the dry, base-state density, and <w'^4> is
1192 : ! calculated by integrating over the PDF. The equation for <w'^4> is:
1193 : !
1194 : ! <w'^4> = mixt_frac * ( 3 * sigma_w_1^4
1195 : ! + 6 * ( mu_w_1 - <w> )^2 * sigma_w_1^2
1196 : ! + ( mu_w_1 - <w> )^4 )
1197 : ! + ( 1 - mixt_frac ) * ( 3 * sigma_w_2^4
1198 : ! + 6 * ( mu_w_2 - <w> )^2 * sigma_w_2^2
1199 : ! + ( mu_w_2 - <w> )^4 ).
1200 : !
1201 : ! The following substitutions are made into the above equation:
1202 : !
1203 : ! mu_w_1 - <w> = sqrt(F_w) * sqrt( ( 1 - mixt_frac ) / mixt_frac )
1204 : ! * sqrt( <w'^2> );
1205 : !
1206 : ! mu_w_2 - <w> = - sqrt(F_w) * sqrt( mixt_frac / ( 1 - mixt_frac ) )
1207 : ! * sqrt( <w'^2> );
1208 : !
1209 : ! sigma_w_1 = sqrt( coef_sigma_w_1_sqd * <w'^2> ); and
1210 : !
1211 : ! sigma_w_2 = sqrt( coef_sigma_w_2_sqd * <w'^2> ).
1212 : !
1213 : ! When w is the setting variable, coef_sigma_w_1_sqd and coef_sigma_w_2_sqd
1214 : ! are given by:
1215 : !
1216 : ! coef_sigma_w_1_sqd = ( ( zeta_w + 1 ) * ( 1 - F_w ) )
1217 : ! / ( ( zeta_w + 2 ) * mixt_frac ); and
1218 : !
1219 : ! coef_sigma_w_2_sqd = ( 1 - F_w ) / ( ( zeta_w + 2 ) * ( 1 - mixt_frac ) ).
1220 : !
1221 : ! When w is a responding variable, coef_sigma_w_1_sqd and coef_sigma_w_2_sqd
1222 : ! are given by:
1223 : !
1224 : ! coef_sigma_w_1_sqd = sqrt( mixt_frac * ( 1 - mixt_frac ) ) * Skw
1225 : ! / ( 3 * mixt_frac * sqrt( F_w ) )
1226 : ! - ( 1 + mixt_frac ) * F_w / ( 3 * mixt_frac )
1227 : ! + 1; and
1228 : !
1229 : ! coef_sigma_w_2_sqd = ( ( 1 - F_w ) - mixt_frac * coef_sigma_w_1_sqd )
1230 : ! / ( 1 - mixt_frac ).
1231 : !
1232 : ! The equation for <w'4> becomes:
1233 : !
1234 : ! <w'^4> = ( 3 * mixt_frac * coef_sigma_w_1_sqd^2
1235 : ! + 6 * F_w * ( 1 - mixt_frac ) * coef_sigma_w_1_sqd
1236 : ! + F_w^2 * ( 1 - mixt_frac )^2 / mixt_frac
1237 : ! + 3 * ( 1 - mixt_frac ) * coef_sigma_w_2_sqd^2
1238 : ! + 6 * F_w * mixt_frac * coef_sigma_w_2_sqd
1239 : ! + F_w^2 * mixt_frac^2 / ( 1 - mixt_frac ) ) * <w'^2>^2.
1240 : !
1241 : ! This equation is of the form:
1242 : !
1243 : ! <w'^4> = coef_wp4_implicit * <w'^2>^2;
1244 : !
1245 : ! where:
1246 : !
1247 : ! coef_wp4_implicit = 3 * mixt_frac * coef_sigma_w_1_sqd^2
1248 : ! + 6 * F_w * ( 1 - mixt_frac ) * coef_sigma_w_1_sqd
1249 : ! + F_w^2 * ( 1 - mixt_frac )^2 / mixt_frac
1250 : ! + 3 * ( 1 - mixt_frac ) * coef_sigma_w_2_sqd^2
1251 : ! + 6 * F_w * mixt_frac * coef_sigma_w_2_sqd
1252 : ! + F_w^2 * mixt_frac^2 / ( 1 - mixt_frac ).
1253 : !
1254 : ! While the <w'^4> term is found in the <w'^3> predictive equation and not
1255 : ! the <w'^2> predictive equation, the <w'^3> and <w'^2> predictive equations
1256 : ! are solved together. This allows the term containing <w'^4> to be solved
1257 : ! implicitly.
1258 :
1259 : ! References:
1260 : !-----------------------------------------------------------------------
1261 :
1262 : use constants_clubb, only: &
1263 : six, & ! Variable(s)
1264 : three, &
1265 : one
1266 :
1267 : use clubb_precision, only: &
1268 : core_rknd ! Procedure(s)
1269 :
1270 : implicit none
1271 :
1272 : integer, intent(in) :: &
1273 : nz
1274 :
1275 : ! Input Variables
1276 : real ( kind = core_rknd ), dimension(nz), intent(in) :: &
1277 : mixt_frac, & ! Mixture fraction [-]
1278 : F_w, & ! Parameter: spread of the PDF comp. means of w [-]
1279 : coef_sigma_w_1_sqd, & ! sigma_w_1^2 = coef_sigma_w_1_sqd * <w'^2> [-]
1280 : coef_sigma_w_2_sqd ! sigma_w_2^2 = coef_sigma_w_2_sqd * <w'^2> [-]
1281 :
1282 : ! Return Variable
1283 : real ( kind = core_rknd ), dimension(nz) :: &
1284 : coef_wp4_implicit ! Coef.: <w'^4> = coef_wp4_implicit * <w'^2>^2 [-]
1285 :
1286 :
1287 : ! Calculate coef_wp4_implicit.
1288 : coef_wp4_implicit = three * mixt_frac * coef_sigma_w_1_sqd**2 &
1289 : + six * F_w * ( one - mixt_frac ) * coef_sigma_w_1_sqd &
1290 : + F_w**2 * ( one - mixt_frac )**2 / mixt_frac &
1291 : + three * ( one - mixt_frac ) * coef_sigma_w_2_sqd**2 &
1292 : + six * F_w * mixt_frac * coef_sigma_w_2_sqd &
1293 0 : + F_w**2 * mixt_frac**2 / ( one - mixt_frac )
1294 :
1295 :
1296 0 : return
1297 :
1298 0 : end function calc_coef_wp4_implicit
1299 :
1300 : !=============================================================================
1301 0 : function calc_coef_wpxp2_implicit( nz, wp2, xp2, wpxp, sgn_wpxp, &
1302 0 : mixt_frac, F_w, F_x, &
1303 0 : coef_sigma_w_1_sqd, &
1304 0 : coef_sigma_w_2_sqd, &
1305 0 : coef_sigma_x_1_sqd, &
1306 0 : coef_sigma_x_2_sqd ) &
1307 0 : result( coef_wpxp2_implicit )
1308 :
1309 : ! Description:
1310 : ! The predictive equation for <x'^2> contains a turbulent advection term of
1311 : ! the form:
1312 : !
1313 : ! - ( 1 / rho_ds ) * d ( rho_ds * <w'x'^2> ) / dz;
1314 : !
1315 : ! where z is height, rho_ds is the dry, base-state density, and <w'x'^2> is
1316 : ! calculated by integrating over the PDF. The equation for <w'x'^2> is:
1317 : !
1318 : ! <w'x'^2>
1319 : ! = mixt_frac * ( ( mu_w_1 - <w> ) * ( ( mu_x_1 - <x> )^2 + sigma_x_1^2 )
1320 : ! + 2 * corr_w_x_1 * sigma_w_1 * sigma_x_1
1321 : ! * ( mu_x_1 - <x> ) )
1322 : ! + ( 1 - mixt_frac ) * ( ( mu_w_2 - <w> )
1323 : ! * ( ( mu_x_2 - <x> )^2 + sigma_x_2^2 )
1324 : ! + 2 * corr_w_x_2 * sigma_w_2 * sigma_x_2
1325 : ! * ( mu_x_2 - <x> ) ).
1326 : !
1327 : ! The following substitutions are made into the above equation:
1328 : !
1329 : ! mu_w_1 - <w> = sqrt(F_w) * sqrt( ( 1 - mixt_frac ) / mixt_frac )
1330 : ! * sqrt( <w'^2> );
1331 : !
1332 : ! mu_w_2 - <w> = - sqrt(F_w) * sqrt( mixt_frac / ( 1 - mixt_frac ) )
1333 : ! * sqrt( <w'^2> );
1334 : !
1335 : ! mu_x_1 - <x> = sqrt(F_x) * sqrt( ( 1 - mixt_frac ) / mixt_frac )
1336 : ! * sqrt( <x'^2> ) * sgn( <w'x'> );
1337 : !
1338 : ! mu_x_2 - <x> = - sqrt(F_x) * sqrt( mixt_frac / ( 1 - mixt_frac ) )
1339 : ! * sqrt( <x'^2> ) * sgn( <w'x'> );
1340 : !
1341 : ! sigma_w_1 = sqrt( coef_sigma_w_1_sqd * <w'^2> );
1342 : !
1343 : ! sigma_w_2 = sqrt( coef_sigma_w_2_sqd * <w'^2> );
1344 : !
1345 : ! sigma_x_1 = sqrt( coef_sigma_x_1_sqd * <x'^2> ); and
1346 : !
1347 : ! sigma_x_2 = sqrt( coef_sigma_x_2_sqd * <x'^2> ).
1348 : !
1349 : ! Either w can be the setting variable and x can be a responding variable,
1350 : ! x can be the setting variable and w can be a responding variable, or both
1351 : ! w and x can be responding variables.
1352 : !
1353 : ! When w is the setting variable, coef_sigma_w_1_sqd and coef_sigma_w_2_sqd
1354 : ! are given by:
1355 : !
1356 : ! coef_sigma_w_1_sqd = ( ( zeta_w + 1 ) * ( 1 - F_w ) )
1357 : ! / ( ( zeta_w + 2 ) * mixt_frac ); and
1358 : !
1359 : ! coef_sigma_w_2_sqd = ( 1 - F_w ) / ( ( zeta_w + 2 ) * ( 1 - mixt_frac ) ).
1360 : !
1361 : ! When w is a responding variable, coef_sigma_w_1_sqd and coef_sigma_w_2_sqd
1362 : ! are given by:
1363 : !
1364 : ! coef_sigma_w_1_sqd = sqrt( mixt_frac * ( 1 - mixt_frac ) ) * Skw
1365 : ! / ( 3 * mixt_frac * sqrt( F_w ) )
1366 : ! - ( 1 + mixt_frac ) * F_w / ( 3 * mixt_frac )
1367 : ! + 1; and
1368 : !
1369 : ! coef_sigma_w_2_sqd = ( ( 1 - F_w ) - mixt_frac * coef_sigma_w_1_sqd )
1370 : ! / ( 1 - mixt_frac ).
1371 : !
1372 : ! When x is the setting variable, coef_sigma_x_1_sqd and coef_sigma_x_2_sqd
1373 : ! are given by:
1374 : !
1375 : ! coef_sigma_x_1_sqd = ( ( zeta_x + 1 ) * ( 1 - F_x ) )
1376 : ! / ( ( zeta_x + 2 ) * mixt_frac ); and
1377 : !
1378 : ! coef_sigma_x_2_sqd = ( 1 - F_x ) / ( ( zeta_x + 2 ) * ( 1 - mixt_frac ) ).
1379 : !
1380 : ! When x is a responding variable, coef_sigma_x_1_sqd and coef_sigma_x_2_sqd
1381 : ! are given by:
1382 : !
1383 : ! coef_sigma_x_1_sqd = sqrt( mixt_frac * ( 1 - mixt_frac ) )
1384 : ! * Skx * sgn( <w'x'> )
1385 : ! / ( 3 * mixt_frac * sqrt( F_x ) )
1386 : ! - ( 1 + mixt_frac ) * F_x / ( 3 * mixt_frac )
1387 : ! + 1; and
1388 : !
1389 : ! coef_sigma_x_2_sqd = ( ( 1 - F_x ) - mixt_frac * coef_sigma_x_1_sqd )
1390 : ! / ( 1 - mixt_frac ).
1391 : !
1392 : ! Additionally:
1393 : !
1394 : ! corr_w_x_1 = corr_w_x_2
1395 : ! = ( <w'x'> - mixt_frac * ( mu_w_1 - <w> ) * ( mu_x_1 - <x> )
1396 : ! - ( 1 - mixt_frac ) * ( mu_w_2 - <w> ) * ( mu_x_2 - <x> ) )
1397 : ! / ( mixt_frac * sigma_w_1 * sigma_x_1
1398 : ! + ( 1 - mixt_frac ) * sigma_w_2 * sigma_x_2 );
1399 : !
1400 : ! where -1 <= corr_w_x_1 = corr_w_x_2 <= 1. This equation can be rewritten
1401 : ! as:
1402 : !
1403 : ! corr_w_x_1 = corr_w_x_2
1404 : ! = ( <w'x'>
1405 : ! - sqrt( F_w ) * sqrt( F_x ) * sgn( <w'x'> )
1406 : ! * sqrt( <w'^2> ) * sqrt( <x'^2 > ) )
1407 : ! / ( ( mixt_frac * sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd )
1408 : ! + ( 1 - mixt_frac )
1409 : ! * sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) )
1410 : ! * sqrt( <w'^2> ) * sqrt( <x'^2> ) ).
1411 : !
1412 : ! The equation for <w'x'^2> becomes:
1413 : !
1414 : ! <w'x'^2>
1415 : ! = sqrt( mixt_frac * ( 1 - mixt_frac ) ) * sqrt( <w'^2> )
1416 : ! * ( sqrt( F_w ) * F_x
1417 : ! * ( ( 1 - mixt_frac ) / mixt_frac - mixt_frac / ( 1 - mixt_frac ) )
1418 : ! + sqrt( F_w ) * ( coef_sigma_x_1_sqd - coef_sigma_x_2_sqd )
1419 : ! + ( 2 * sqrt( F_x ) * sgn( <w'x'> ) * <w'x'>
1420 : ! * ( sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd )
1421 : ! - sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) ) )
1422 : ! / ( ( mixt_frac * sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd )
1423 : ! + ( 1 - mixt_frac )
1424 : ! * sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) )
1425 : ! * sqrt( <w'^2> ) * sqrt( <x'^2> ) )
1426 : ! - ( 2 * sqrt( F_w ) * F_x
1427 : ! * ( sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd )
1428 : ! - sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) ) )
1429 : ! / ( mixt_frac * sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd )
1430 : ! + ( 1 - mixt_frac )
1431 : ! * sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) ) )
1432 : ! * <x'^2>
1433 : !
1434 : ! This equation is of the form:
1435 : !
1436 : ! <w'x'^2> = coef_wpxp2_implicit * <x'^2>;
1437 : !
1438 : ! where:
1439 : !
1440 : ! coef_wpxp2_implicit
1441 : ! = sqrt( mixt_frac * ( 1 - mixt_frac ) ) * sqrt( <w'^2> )
1442 : ! * ( sqrt( F_w ) * F_x
1443 : ! * ( ( 1 - mixt_frac ) / mixt_frac - mixt_frac / ( 1 - mixt_frac ) )
1444 : ! + sqrt( F_w ) * ( coef_sigma_x_1_sqd - coef_sigma_x_2_sqd )
1445 : ! + ( 2 * sqrt( F_x ) * sgn( <w'x'> ) * <w'x'>
1446 : ! * ( sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd )
1447 : ! - sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) ) )
1448 : ! / ( ( mixt_frac * sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd )
1449 : ! + ( 1 - mixt_frac )
1450 : ! * sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) )
1451 : ! * sqrt( <w'^2> ) * sqrt( <x'^2> ) )
1452 : ! - ( 2 * sqrt( F_w ) * F_x
1453 : ! * ( sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd )
1454 : ! - sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) ) )
1455 : ! / ( mixt_frac * sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd )
1456 : ! + ( 1 - mixt_frac )
1457 : ! * sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) ) ).
1458 : !
1459 : ! In the special case that coef_sigma_w_1_sqd * coef_sigma_x_1_sqd = 0 and
1460 : ! coef_sigma_w_2_sqd * coef_sigma_x_2_sqd = 0, the above equation is
1461 : ! undefined. However, the equation for this special case can be derived by
1462 : ! taking the original equation for <w'x'^2> and setting both
1463 : ! sigma_w_1 * sigma_x_1 = 0 and sigma_w_2 * sigma_x_2 = 0. The equation
1464 : ! becomes:
1465 : !
1466 : ! <w'x'^2>
1467 : ! = mixt_frac * ( ( mu_w_1 - <w> ) * ( ( mu_x_1 - <x> )^2 + sigma_x_1^2 )
1468 : ! + ( 1 - mixt_frac ) * ( ( mu_w_2 - <w> )
1469 : ! * ( ( mu_x_2 - <x> )^2 + sigma_x_2^2 );
1470 : !
1471 : ! and making the same substitutions as before, it can be rewritten as:
1472 : !
1473 : ! <w'x'^2>
1474 : ! = sqrt( mixt_frac * ( 1 - mixt_frac ) ) * sqrt( <w'^2> ) * sqrt( F_w )
1475 : ! * ( F_x * ( ( 1 - mixt_frac ) / mixt_frac
1476 : ! - mixt_frac / ( 1 - mixt_frac ) )
1477 : ! + ( coef_sigma_x_1_sqd - coef_sigma_x_2_sqd ) ) * <x'^2>.
1478 : !
1479 : ! The coefficient in this special case is:
1480 : !
1481 : ! coef_wpxp2_implicit
1482 : ! = sqrt( mixt_frac * ( 1 - mixt_frac ) ) * sqrt( <w'^2> ) * sqrt( F_w )
1483 : ! * ( F_x * ( ( 1 - mixt_frac ) / mixt_frac
1484 : ! - mixt_frac / ( 1 - mixt_frac ) )
1485 : ! + ( coef_sigma_x_1_sqd - coef_sigma_x_2_sqd ) ).
1486 :
1487 : ! References:
1488 : !-----------------------------------------------------------------------
1489 :
1490 : use constants_clubb, only: &
1491 : two, & ! Variable(s)
1492 : one, &
1493 : zero
1494 :
1495 : use clubb_precision, only: &
1496 : core_rknd ! Procedure(s)
1497 :
1498 : implicit none
1499 :
1500 : integer, intent(in) :: &
1501 : nz
1502 :
1503 : ! Input Variables
1504 : real ( kind = core_rknd ), dimension(nz), intent(in) :: &
1505 : wp2, & ! Variance of w (overall) [m^2/s^2]
1506 : xp2, & ! Variance of x (overall) [(units vary)^2]
1507 : wpxp, & ! Covariance of w and x [m/s (units vary)]
1508 : sgn_wpxp, & ! Sign of the covariance of w and x [-]
1509 : mixt_frac, & ! Mixture fraction [-]
1510 : F_w, & ! Parameter: spread of the PDF comp. means of w [-]
1511 : F_x, & ! Parameter: spread of the PDF comp. means of x [-]
1512 : coef_sigma_w_1_sqd, & ! sigma_w_1^2 = coef_sigma_w_1_sqd * <w'^2> [-]
1513 : coef_sigma_w_2_sqd, & ! sigma_w_2^2 = coef_sigma_w_2_sqd * <w'^2> [-]
1514 : coef_sigma_x_1_sqd, & ! sigma_x_1^2 = coef_sigma_x_1_sqd * <x'^2> [-]
1515 : coef_sigma_x_2_sqd ! sigma_x_2^2 = coef_sigma_x_2_sqd * <x'^2> [-]
1516 :
1517 : ! Return Variable
1518 : real ( kind = core_rknd ), dimension(nz) :: &
1519 : coef_wpxp2_implicit ! Coef.: <w'x'^2> = coef_wpxp2_implicit * <x'^2> [m/s]
1520 :
1521 : ! Local Variable
1522 : real ( kind = core_rknd ), dimension(nz) :: &
1523 0 : coefs_factor ! Factor involving coef_sigma_... coefficients [-]
1524 :
1525 :
1526 : ! Calculate coef_wpxp2_implicit.
1527 : where ( ( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd > zero &
1528 : .or. coef_sigma_w_2_sqd * coef_sigma_x_2_sqd > zero ) &
1529 0 : .and. ( wp2 * xp2 > zero ) )
1530 :
1531 : coefs_factor &
1532 : = ( sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd ) &
1533 : - sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) ) &
1534 : / ( mixt_frac * sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd ) &
1535 : + ( one - mixt_frac ) &
1536 : * sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) )
1537 :
1538 : coef_wpxp2_implicit &
1539 : = sqrt( mixt_frac * ( one - mixt_frac ) ) * sqrt( wp2 ) &
1540 : * ( sqrt( F_w ) * F_x &
1541 : * ( ( one - mixt_frac ) / mixt_frac &
1542 : - mixt_frac / ( one - mixt_frac ) ) &
1543 : + sqrt( F_w ) * ( coef_sigma_x_1_sqd - coef_sigma_x_2_sqd ) &
1544 : + two * sqrt( F_x ) * coefs_factor * sgn_wpxp * wpxp &
1545 : / ( sqrt( wp2 ) * sqrt( xp2 ) ) &
1546 : - two * sqrt( F_w ) * F_x * coefs_factor )
1547 :
1548 : elsewhere ! ( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd = 0
1549 : ! and coef_sigma_w_2_sqd * coef_sigma_x_2_sqd = 0 )
1550 : ! or wp2 * xp2 = 0
1551 :
1552 : coef_wpxp2_implicit &
1553 : = sqrt( mixt_frac * ( one - mixt_frac ) ) * sqrt( wp2 ) * sqrt( F_w ) &
1554 : * ( F_x * ( ( one - mixt_frac ) / mixt_frac &
1555 : - mixt_frac / ( one - mixt_frac ) ) &
1556 : + ( coef_sigma_x_1_sqd - coef_sigma_x_2_sqd ) )
1557 :
1558 : endwhere
1559 :
1560 :
1561 0 : return
1562 :
1563 0 : end function calc_coef_wpxp2_implicit
1564 :
1565 : !=============================================================================
1566 0 : subroutine calc_coefs_wp2xp_semiimpl( nz, wp2, xp2, sgn_wpxp, & ! In
1567 0 : mixt_frac, F_w, F_x, & ! In
1568 0 : coef_sigma_w_1_sqd, & ! In
1569 0 : coef_sigma_w_2_sqd, & ! In
1570 0 : coef_sigma_x_1_sqd, & ! In
1571 0 : coef_sigma_x_2_sqd, & ! In
1572 0 : coef_wp2xp_implicit, & ! Out
1573 0 : term_wp2xp_explicit ) ! Out
1574 :
1575 : ! Description:
1576 : ! The predictive equation for <w'x'> contains a turbulent advection term of
1577 : ! the form:
1578 : !
1579 : ! - ( 1 / rho_ds ) * d ( rho_ds * <w'^2 x'> ) / dz;
1580 : !
1581 : ! where z is height, rho_ds is the dry, base-state density, and <w'^2 x'> is
1582 : ! calculated by integrating over the PDF. The equation for <w'^2 x'> is:
1583 : !
1584 : ! <w'^2 x'>
1585 : ! = mixt_frac * ( ( mu_x_1 - <x> ) * ( ( mu_w_1 - <w> )^2 + sigma_w_1^2 )
1586 : ! + 2 * corr_w_x_1 * sigma_w_1 * sigma_x_1
1587 : ! * ( mu_w_1 - <w> ) )
1588 : ! + ( 1 - mixt_frac ) * ( ( mu_x_2 - <x> )
1589 : ! * ( ( mu_w_2 - <w> )^2 + sigma_w_2^2 )
1590 : ! + 2 * corr_w_x_2 * sigma_w_2 * sigma_x_2
1591 : ! * ( mu_w_2 - <w> ) ).
1592 : !
1593 : ! The following substitutions are made into the above equation:
1594 : !
1595 : ! mu_w_1 - <w> = sqrt(F_w) * sqrt( ( 1 - mixt_frac ) / mixt_frac )
1596 : ! * sqrt( <w'^2> );
1597 : !
1598 : ! mu_w_2 - <w> = - sqrt(F_w) * sqrt( mixt_frac / ( 1 - mixt_frac ) )
1599 : ! * sqrt( <w'^2> );
1600 : !
1601 : ! mu_x_1 - <x> = sqrt(F_x) * sqrt( ( 1 - mixt_frac ) / mixt_frac )
1602 : ! * sqrt( <x'^2> ) * sgn( <w'x'> );
1603 : !
1604 : ! mu_x_2 - <x> = - sqrt(F_x) * sqrt( mixt_frac / ( 1 - mixt_frac ) )
1605 : ! * sqrt( <x'^2> ) * sgn( <w'x'> );
1606 : !
1607 : ! sigma_w_1 = sqrt( coef_sigma_w_1_sqd * <w'^2> );
1608 : !
1609 : ! sigma_w_2 = sqrt( coef_sigma_w_2_sqd * <w'^2> );
1610 : !
1611 : ! sigma_x_1 = sqrt( coef_sigma_x_1_sqd * <x'^2> ); and
1612 : !
1613 : ! sigma_x_2 = sqrt( coef_sigma_x_2_sqd * <x'^2> ).
1614 : !
1615 : ! Either w can be the setting variable and x can be a responding variable,
1616 : ! x can be the setting variable and w can be a responding variable, or both
1617 : ! w and x can be responding variables.
1618 : !
1619 : ! When w is the setting variable, coef_sigma_w_1_sqd and coef_sigma_w_2_sqd
1620 : ! are given by:
1621 : !
1622 : ! coef_sigma_w_1_sqd = ( ( zeta_w + 1 ) * ( 1 - F_w ) )
1623 : ! / ( ( zeta_w + 2 ) * mixt_frac ); and
1624 : !
1625 : ! coef_sigma_w_2_sqd = ( 1 - F_w ) / ( ( zeta_w + 2 ) * ( 1 - mixt_frac ) ).
1626 : !
1627 : ! When w is a responding variable, coef_sigma_w_1_sqd and coef_sigma_w_2_sqd
1628 : ! are given by:
1629 : !
1630 : ! coef_sigma_w_1_sqd = sqrt( mixt_frac * ( 1 - mixt_frac ) ) * Skw
1631 : ! / ( 3 * mixt_frac * sqrt( F_w ) )
1632 : ! - ( 1 + mixt_frac ) * F_w / ( 3 * mixt_frac )
1633 : ! + 1; and
1634 : !
1635 : ! coef_sigma_w_2_sqd = ( ( 1 - F_w ) - mixt_frac * coef_sigma_w_1_sqd )
1636 : ! / ( 1 - mixt_frac ).
1637 : !
1638 : ! When x is the setting variable, coef_sigma_x_1_sqd and coef_sigma_x_2_sqd
1639 : ! are given by:
1640 : !
1641 : ! coef_sigma_x_1_sqd = ( ( zeta_x + 1 ) * ( 1 - F_x ) )
1642 : ! / ( ( zeta_x + 2 ) * mixt_frac ); and
1643 : !
1644 : ! coef_sigma_x_2_sqd = ( 1 - F_x ) / ( ( zeta_x + 2 ) * ( 1 - mixt_frac ) ).
1645 : !
1646 : ! When x is a responding variable, coef_sigma_x_1_sqd and coef_sigma_x_2_sqd
1647 : ! are given by:
1648 : !
1649 : ! coef_sigma_x_1_sqd = sqrt( mixt_frac * ( 1 - mixt_frac ) )
1650 : ! * Skx * sgn( <w'x'> )
1651 : ! / ( 3 * mixt_frac * sqrt( F_x ) )
1652 : ! - ( 1 + mixt_frac ) * F_x / ( 3 * mixt_frac )
1653 : ! + 1; and
1654 : !
1655 : ! coef_sigma_x_2_sqd = ( ( 1 - F_x ) - mixt_frac * coef_sigma_x_1_sqd )
1656 : ! / ( 1 - mixt_frac ).
1657 : !
1658 : ! Additionally:
1659 : !
1660 : ! corr_w_x_1 = corr_w_x_2
1661 : ! = ( <w'x'> - mixt_frac * ( mu_w_1 - <w> ) * ( mu_x_1 - <x> )
1662 : ! - ( 1 - mixt_frac ) * ( mu_w_2 - <w> ) * ( mu_x_2 - <x> ) )
1663 : ! / ( mixt_frac * sigma_w_1 * sigma_x_1
1664 : ! + ( 1 - mixt_frac ) * sigma_w_2 * sigma_x_2 );
1665 : !
1666 : ! where -1 <= corr_w_x_1 = corr_w_x_2 <= 1. This equation can be rewritten
1667 : ! as:
1668 : !
1669 : ! corr_w_x_1 = corr_w_x_2
1670 : ! = ( <w'x'>
1671 : ! - sqrt( F_w ) * sqrt( F_x ) * sgn( <w'x'> )
1672 : ! * sqrt( <w'^2> ) * sqrt( <x'^2 > ) )
1673 : ! / ( ( mixt_frac * sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd )
1674 : ! + ( 1 - mixt_frac )
1675 : ! * sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) )
1676 : ! * sqrt( <w'^2> ) * sqrt( <x'^2> ) ).
1677 : !
1678 : ! The equation for <w'^2 x'> becomes:
1679 : !
1680 : ! <w'^2 x'>
1681 : ! = sqrt( mixt_frac * ( 1 - mixt_frac ) ) * sqrt( F_x )
1682 : ! * sqrt( <x'^2> ) * <w'^2> * sgn( <w'x'> )
1683 : ! * ( F_w * ( ( 1 - mixt_frac ) / mixt_frac
1684 : ! - mixt_frac / ( 1 - mixt_frac ) )
1685 : ! + ( coef_sigma_w_1_sqd - coef_sigma_w_2_sqd )
1686 : ! - ( 2 * F_w * ( sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd )
1687 : ! - sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) )
1688 : ! / ( mixt_frac * sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd )
1689 : ! + ( 1 - mixt_frac )
1690 : ! * sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) ) ) )
1691 : ! + ( sqrt( mixt_frac * ( 1 - mixt_frac ) )
1692 : ! * 2 * sqrt( F_w ) * sqrt( <w'^2> )
1693 : ! * ( sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd )
1694 : ! - sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) )
1695 : ! / ( mixt_frac * sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd )
1696 : ! + ( 1 - mixt_frac )
1697 : ! * sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) ) )
1698 : ! * <w'x'>
1699 : !
1700 : ! This equation is of the form:
1701 : !
1702 : ! <w'^2 x'> = coef_wp2xp_implicit * <w'x'> + term_wp2xp_explicit;
1703 : !
1704 : ! where:
1705 : !
1706 : ! coef_wp2xp_implicit
1707 : ! = sqrt( mixt_frac * ( 1 - mixt_frac ) ) * 2 * sqrt( F_w ) * sqrt( <w'^2> )
1708 : ! * ( sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd )
1709 : ! - sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) )
1710 : ! / ( mixt_frac * sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd )
1711 : ! + ( 1 - mixt_frac )
1712 : ! * sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) ); and
1713 : !
1714 : ! term_wp2xp_explicit
1715 : ! = sqrt( mixt_frac * ( 1 - mixt_frac ) ) * sqrt( F_x )
1716 : ! * sqrt( <x'^2> ) * <w'^2> * sgn( <w'x'> )
1717 : ! * ( F_w * ( ( 1 - mixt_frac ) / mixt_frac
1718 : ! - mixt_frac / ( 1 - mixt_frac ) )
1719 : ! + ( coef_sigma_w_1_sqd - coef_sigma_w_2_sqd )
1720 : ! - ( 2 * F_w * ( sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd )
1721 : ! - sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) )
1722 : ! / ( mixt_frac * sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd )
1723 : ! + ( 1 - mixt_frac )
1724 : ! * sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) ) ) ).
1725 : !
1726 : ! In the special case that coef_sigma_w_1_sqd * coef_sigma_x_1_sqd = 0 and
1727 : ! coef_sigma_w_2_sqd * coef_sigma_x_2_sqd = 0, the above equation is
1728 : ! undefined. However, the equation for this special case can be derived by
1729 : ! taking the original equation for <w'^2 x'> and setting both
1730 : ! sigma_w_1 * sigma_x_1 = 0 and sigma_w_2 * sigma_x_2 = 0. The equation
1731 : ! becomes:
1732 : !
1733 : ! <w'^2 x'>
1734 : ! = mixt_frac * ( ( mu_x_1 - <x> ) * ( ( mu_w_1 - <w> )^2 + sigma_w_1^2 )
1735 : ! + ( 1 - mixt_frac ) * ( ( mu_x_2 - <x> )
1736 : ! * ( ( mu_w_2 - <w> )^2 + sigma_w_2^2 );
1737 : !
1738 : ! and making the same substitutions as before, it can be rewritten as:
1739 : !
1740 : ! <w'^2 x'>
1741 : ! = sqrt( mixt_frac * ( 1 - mixt_frac ) )
1742 : ! * sqrt( F_x ) * sqrt( <x'^2> ) * <w'^2> * sgn( <w'x'> )
1743 : ! * ( F_w * ( ( 1 - mixt_frac ) / mixt_frac
1744 : ! - mixt_frac / ( 1 - mixt_frac ) )
1745 : ! + ( coef_sigma_w_1_sqd - coef_sigma_w_2_sqd ) ).
1746 : !
1747 : ! Likewise, the equation for <w'x'> in this special case becomes:
1748 : !
1749 : ! <w'x'> = mixt_frac * ( mu_w_1 - <w> ) * ( mu_x_1 - <x> )
1750 : ! + ( 1 - mixt_frac ) * ( mu_w_2 - <w> ) * ( mu_x_2 - <x> );
1751 : !
1752 : ! and making the same substitutions as before, it can be rewritten as:
1753 : !
1754 : ! <w'x'> = sqrt( F_w ) * sqrt( F_x )
1755 : ! * sqrt( <w'^2> ) * sqrt( <x'^2> ) * sgn( <w'x'> ).
1756 : !
1757 : ! The equation for <w'^2 x'> can be rewritten as:
1758 : !
1759 : ! <w'^2 x'>
1760 : ! = sqrt( mixt_frac * ( 1 - mixt_frac ) ) * sqrt( F_w ) * sqrt( <w'^2> )
1761 : ! * ( ( 1 - mixt_frac ) / mixt_frac - mixt_frac / ( 1 - mixt_frac ) )
1762 : ! * <w'x'>
1763 : ! + sqrt( mixt_frac * ( 1 - mixt_frac ) )
1764 : ! * sqrt( F_x ) * sqrt( <x'^2> ) * <w'^2> * sgn( <w'x'> )
1765 : ! * ( coef_sigma_w_1_sqd - coef_sigma_w_2_sqd ).
1766 : !
1767 : ! The coefficients in this special case are:
1768 : !
1769 : ! coef_wp2xp_implicit
1770 : ! = sqrt( mixt_frac * ( 1 - mixt_frac ) ) * sqrt( F_w ) * sqrt( <w'^2> )
1771 : ! * ( ( 1 - mixt_frac ) / mixt_frac - mixt_frac / ( 1 - mixt_frac ) ); and
1772 : !
1773 : ! term_wp2xp_explicit
1774 : ! = sqrt( mixt_frac * ( 1 - mixt_frac ) )
1775 : ! * sqrt( F_x ) * sqrt( <x'^2> ) * <w'^2> * sgn( <w'x'> )
1776 : ! * ( coef_sigma_w_1_sqd - coef_sigma_w_2_sqd ).
1777 :
1778 : ! References:
1779 : !-----------------------------------------------------------------------
1780 :
1781 : use constants_clubb, only: &
1782 : two, & ! Variable(s)
1783 : one, &
1784 : zero
1785 :
1786 : use clubb_precision, only: &
1787 : core_rknd ! Procedure(s)
1788 :
1789 : implicit none
1790 :
1791 : integer, intent(in) :: &
1792 : nz
1793 :
1794 : ! Input Variables
1795 : real ( kind = core_rknd ), dimension(nz), intent(in) :: &
1796 : wp2, & ! Variance of w (overall) [m^2/s^2]
1797 : xp2, & ! Variance of x (overall) [(units vary)^2]
1798 : sgn_wpxp, & ! Sign of the covariance of w and x [-]
1799 : mixt_frac, & ! Mixture fraction [-]
1800 : F_w, & ! Parameter: spread of the PDF comp. means of w [-]
1801 : F_x, & ! Parameter: spread of the PDF comp. means of x [-]
1802 : coef_sigma_w_1_sqd, & ! sigma_w_1^2 = coef_sigma_w_1_sqd * <w'^2> [-]
1803 : coef_sigma_w_2_sqd, & ! sigma_w_2^2 = coef_sigma_w_2_sqd * <w'^2> [-]
1804 : coef_sigma_x_1_sqd, & ! sigma_x_1^2 = coef_sigma_x_1_sqd * <x'^2> [-]
1805 : coef_sigma_x_2_sqd ! sigma_x_2^2 = coef_sigma_x_2_sqd * <x'^2> [-]
1806 :
1807 : ! Output Variables
1808 : ! Coefs.: <w'^2 x'> = coef_wp2xp_implicit * <w'x'> + term_wp2xp_explicit
1809 : real ( kind = core_rknd ), dimension(nz), intent(out) :: &
1810 : coef_wp2xp_implicit, & ! Coefficient that is multiplied by <w'x'> [m/s]
1811 : term_wp2xp_explicit ! Term that is on the RHS [m^2/s^2 (units vary)]
1812 :
1813 : ! Local Variable
1814 : real ( kind = core_rknd ), dimension(nz) :: &
1815 0 : coefs_factor ! Factor involving coef_sigma_... coefficients [-]
1816 :
1817 :
1818 : ! Calculate coef_wp2xp_implicit and term_wp2xp_explicit.
1819 : where ( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd > zero &
1820 0 : .or. coef_sigma_w_2_sqd * coef_sigma_x_2_sqd > zero )
1821 :
1822 : coefs_factor &
1823 : = ( sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd ) &
1824 : - sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) ) &
1825 : / ( mixt_frac * sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd ) &
1826 : + ( one - mixt_frac ) &
1827 : * sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) )
1828 :
1829 : coef_wp2xp_implicit &
1830 : = sqrt( mixt_frac * ( one - mixt_frac ) ) &
1831 : * two * sqrt( F_w ) * sqrt( wp2 ) * coefs_factor
1832 :
1833 : term_wp2xp_explicit &
1834 : = sqrt( mixt_frac * ( one - mixt_frac ) ) &
1835 : * sqrt( F_x ) * sqrt( xp2 ) * wp2 * sgn_wpxp &
1836 : * ( F_w * ( ( one - mixt_frac ) / mixt_frac &
1837 : - mixt_frac / ( one - mixt_frac ) ) &
1838 : + ( coef_sigma_w_1_sqd - coef_sigma_w_2_sqd ) &
1839 : - two * F_w * coefs_factor )
1840 :
1841 : elsewhere ! coef_sigma_w_1_sqd * coef_sigma_x_1_sqd = 0
1842 : ! and coef_sigma_w_2_sqd * coef_sigma_x_2_sqd = 0
1843 :
1844 : coef_wp2xp_implicit &
1845 : = sqrt( mixt_frac * ( one - mixt_frac ) ) &
1846 : * sqrt( F_w ) * sqrt( wp2 ) &
1847 : * ( ( one - mixt_frac ) / mixt_frac &
1848 : - mixt_frac / ( one - mixt_frac ) )
1849 :
1850 : term_wp2xp_explicit &
1851 : = sqrt( mixt_frac * ( one - mixt_frac ) ) &
1852 : * sqrt( F_x ) * sqrt( xp2 ) * wp2 * sgn_wpxp &
1853 : * ( coef_sigma_w_1_sqd - coef_sigma_w_2_sqd )
1854 :
1855 : endwhere
1856 :
1857 :
1858 0 : return
1859 :
1860 : end subroutine calc_coefs_wp2xp_semiimpl
1861 :
1862 : !=============================================================================
1863 0 : subroutine calc_coefs_wpxpyp_semiimpl( nz, wp2, xp2, yp2, wpxp, & ! In
1864 0 : wpyp, sgn_wpxp, sgn_wpyp, & ! In
1865 0 : mixt_frac, F_w, F_x, F_y, & ! In
1866 0 : coef_sigma_w_1_sqd, & ! In
1867 0 : coef_sigma_w_2_sqd, & ! In
1868 0 : coef_sigma_x_1_sqd, & ! In
1869 0 : coef_sigma_x_2_sqd, & ! In
1870 0 : coef_sigma_y_1_sqd, & ! In
1871 0 : coef_sigma_y_2_sqd, & ! In
1872 0 : coef_wpxpyp_implicit, & ! Out
1873 0 : term_wpxpyp_explicit ) ! Out
1874 :
1875 : ! Description:
1876 : ! The predictive equation for <w'x'y'> contains a turbulent advection term
1877 : ! of the form:
1878 : !
1879 : ! - ( 1 / rho_ds ) * d ( rho_ds * <w'x'y'> ) / dz;
1880 : !
1881 : ! where z is height, rho_ds is the dry, base-state density, and <w'x'y'> is
1882 : ! calculated by integrating over the PDF. The equation for <w'x'y'> is:
1883 : !
1884 : ! <w'x'y'>
1885 : ! = mixt_frac
1886 : ! * ( ( mu_w_1 - <w> ) * ( mu_x_1 - <x> ) * ( mu_y_1 - <y> )
1887 : ! + corr_x_y_1 * sigma_x_1 * sigma_y_1 * ( mu_w_1 - <w> )
1888 : ! + corr_w_y_1 * sigma_w_1 * sigma_y_1 * ( mu_x_1 - <x> )
1889 : ! + corr_w_x_1 * sigma_w_1 * sigma_x_1 * ( mu_y_1 - <y> ) )
1890 : ! + ( 1 - mixt_frac )
1891 : ! * ( ( mu_w_2 - <w> ) * ( mu_x_2 - <x> ) * ( mu_y_2 - <y> )
1892 : ! + corr_x_y_2 * sigma_x_2 * sigma_y_2 * ( mu_w_2 - <w> )
1893 : ! + corr_w_y_2 * sigma_w_2 * sigma_y_2 * ( mu_x_2 - <x> )
1894 : ! + corr_w_x_2 * sigma_w_2 * sigma_x_2 * ( mu_y_2 - <y> ) ).
1895 : !
1896 : ! The following substitutions are made into the above equation:
1897 : !
1898 : ! mu_w_1 - <w> = sqrt(F_w) * sqrt( ( 1 - mixt_frac ) / mixt_frac )
1899 : ! * sqrt( <w'^2> );
1900 : !
1901 : ! mu_w_2 - <w> = - sqrt(F_w) * sqrt( mixt_frac / ( 1 - mixt_frac ) )
1902 : ! * sqrt( <w'^2> );
1903 : !
1904 : ! mu_x_1 - <x> = sqrt(F_x) * sqrt( ( 1 - mixt_frac ) / mixt_frac )
1905 : ! * sqrt( <x'^2> ) * sgn( <w'x'> );
1906 : !
1907 : ! mu_x_2 - <x> = - sqrt(F_x) * sqrt( mixt_frac / ( 1 - mixt_frac ) )
1908 : ! * sqrt( <x'^2> ) * sgn( <w'x'> );
1909 : !
1910 : ! mu_y_1 - <y> = sqrt(F_y) * sqrt( ( 1 - mixt_frac ) / mixt_frac )
1911 : ! * sqrt( <y'^2> ) * sgn( <w'y'> );
1912 : !
1913 : ! mu_y_2 - <y> = - sqrt(F_y) * sqrt( mixt_frac / ( 1 - mixt_frac ) )
1914 : ! * sqrt( <y'^2> ) * sgn( <w'y'> );
1915 : !
1916 : ! sigma_w_1 = sqrt( coef_sigma_w_1_sqd * <w'^2> );
1917 : !
1918 : ! sigma_w_2 = sqrt( coef_sigma_w_2_sqd * <w'^2> );
1919 : !
1920 : ! sigma_x_1 = sqrt( coef_sigma_x_1_sqd * <x'^2> );
1921 : !
1922 : ! sigma_x_2 = sqrt( coef_sigma_x_2_sqd * <x'^2> );
1923 : !
1924 : ! sigma_y_1 = sqrt( coef_sigma_y_1_sqd * <y'^2> ); and
1925 : !
1926 : ! sigma_y_2 = sqrt( coef_sigma_y_2_sqd * <y'^2> ).
1927 : !
1928 : ! Either w can be the setting variable and both x and y can be responding
1929 : ! variables, x can be the setting variable and both w and y can be
1930 : ! responding variables, y can be the setting variable and both w and x can
1931 : ! be responding variables, or all of w, x, and y can be responding
1932 : ! variables.
1933 : !
1934 : ! When w is the setting variable, coef_sigma_w_1_sqd and coef_sigma_w_2_sqd
1935 : ! are given by:
1936 : !
1937 : ! coef_sigma_w_1_sqd = ( ( zeta_w + 1 ) * ( 1 - F_w ) )
1938 : ! / ( ( zeta_w + 2 ) * mixt_frac ); and
1939 : !
1940 : ! coef_sigma_w_2_sqd = ( 1 - F_w ) / ( ( zeta_w + 2 ) * ( 1 - mixt_frac ) ).
1941 : !
1942 : ! When w is a responding variable, coef_sigma_w_1_sqd and coef_sigma_w_2_sqd
1943 : ! are given by:
1944 : !
1945 : ! coef_sigma_w_1_sqd = sqrt( mixt_frac * ( 1 - mixt_frac ) ) * Skw
1946 : ! / ( 3 * mixt_frac * sqrt( F_w ) )
1947 : ! - ( 1 + mixt_frac ) * F_w / ( 3 * mixt_frac )
1948 : ! + 1; and
1949 : !
1950 : ! coef_sigma_w_2_sqd = ( ( 1 - F_w ) - mixt_frac * coef_sigma_w_1_sqd )
1951 : ! / ( 1 - mixt_frac ).
1952 : !
1953 : ! When x is the setting variable, coef_sigma_x_1_sqd and coef_sigma_x_2_sqd
1954 : ! are given by:
1955 : !
1956 : ! coef_sigma_x_1_sqd = ( ( zeta_x + 1 ) * ( 1 - F_x ) )
1957 : ! / ( ( zeta_x + 2 ) * mixt_frac ); and
1958 : !
1959 : ! coef_sigma_x_2_sqd = ( 1 - F_x ) / ( ( zeta_x + 2 ) * ( 1 - mixt_frac ) ).
1960 : !
1961 : ! When x is a responding variable, coef_sigma_x_1_sqd and coef_sigma_x_2_sqd
1962 : ! are given by:
1963 : !
1964 : ! coef_sigma_x_1_sqd = sqrt( mixt_frac * ( 1 - mixt_frac ) )
1965 : ! * Skx * sgn( <w'x'> )
1966 : ! / ( 3 * mixt_frac * sqrt( F_x ) )
1967 : ! - ( 1 + mixt_frac ) * F_x / ( 3 * mixt_frac )
1968 : ! + 1; and
1969 : !
1970 : ! coef_sigma_x_2_sqd = ( ( 1 - F_x ) - mixt_frac * coef_sigma_x_1_sqd )
1971 : ! / ( 1 - mixt_frac ).
1972 : !
1973 : ! When y is the setting variable, coef_sigma_y_1_sqd and coef_sigma_y_2_sqd
1974 : ! are given by:
1975 : !
1976 : ! coef_sigma_y_1_sqd = ( ( zeta_y + 1 ) * ( 1 - F_y ) )
1977 : ! / ( ( zeta_y + 2 ) * mixt_frac ); and
1978 : !
1979 : ! coef_sigma_y_2_sqd = ( 1 - F_y ) / ( ( zeta_y + 2 ) * ( 1 - mixt_frac ) ).
1980 : !
1981 : ! When y is a responding variable, coef_sigma_y_1_sqd and coef_sigma_y_2_sqd
1982 : ! are given by:
1983 : !
1984 : ! coef_sigma_y_1_sqd = sqrt( mixt_frac * ( 1 - mixt_frac ) )
1985 : ! * Sky * sgn( <w'y'> )
1986 : ! / ( 3 * mixt_frac * sqrt( F_y ) )
1987 : ! - ( 1 + mixt_frac ) * F_y / ( 3 * mixt_frac )
1988 : ! + 1; and
1989 : !
1990 : ! coef_sigma_y_2_sqd = ( ( 1 - F_y ) - mixt_frac * coef_sigma_y_1_sqd )
1991 : ! / ( 1 - mixt_frac ).
1992 : !
1993 : ! Additionally:
1994 : !
1995 : ! corr_w_x_1 = corr_w_x_2
1996 : ! = ( <w'x'> - mixt_frac * ( mu_w_1 - <w> ) * ( mu_x_1 - <x> )
1997 : ! - ( 1 - mixt_frac ) * ( mu_w_2 - <w> ) * ( mu_x_2 - <x> ) )
1998 : ! / ( mixt_frac * sigma_w_1 * sigma_x_1
1999 : ! + ( 1 - mixt_frac ) * sigma_w_2 * sigma_x_2 );
2000 : !
2001 : ! where -1 <= corr_w_x_1 = corr_w_x_2 <= 1. This equation can be rewritten
2002 : ! as:
2003 : !
2004 : ! corr_w_x_1 = corr_w_x_2
2005 : ! = ( <w'x'>
2006 : ! - sqrt( F_w ) * sqrt( F_x ) * sgn( <w'x'> )
2007 : ! * sqrt( <w'^2> ) * sqrt( <x'^2 > ) )
2008 : ! / ( ( mixt_frac * sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd )
2009 : ! + ( 1 - mixt_frac )
2010 : ! * sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) )
2011 : ! * sqrt( <w'^2> ) * sqrt( <x'^2> ) ).
2012 : !
2013 : ! Likewise:
2014 : !
2015 : ! corr_w_y_1 = corr_w_y_2
2016 : ! = ( <w'y'>
2017 : ! - sqrt( F_w ) * sqrt( F_y ) * sgn( <w'y'> )
2018 : ! * sqrt( <w'^2> ) * sqrt( <y'^2 > ) )
2019 : ! / ( ( mixt_frac * sqrt( coef_sigma_w_1_sqd * coef_sigma_y_1_sqd )
2020 : ! + ( 1 - mixt_frac )
2021 : ! * sqrt( coef_sigma_w_2_sqd * coef_sigma_y_2_sqd ) )
2022 : ! * sqrt( <w'^2> ) * sqrt( <y'^2> ) ); and
2023 : !
2024 : ! corr_x_y_1 = corr_x_y_2
2025 : ! = ( <x'y'>
2026 : ! - sqrt( F_x ) * sqrt( F_y ) * sgn( <w'x'> ) * sgn( <w'y'> )
2027 : ! * sqrt( <x'^2> ) * sqrt( <y'^2 > ) )
2028 : ! / ( ( mixt_frac * sqrt( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd )
2029 : ! + ( 1 - mixt_frac )
2030 : ! * sqrt( coef_sigma_x_2_sqd * coef_sigma_y_2_sqd ) )
2031 : ! * sqrt( <x'^2> ) * sqrt( <y'^2> ) ).
2032 : !
2033 : ! The equation for <w'x'y'> becomes:
2034 : !
2035 : ! <w'x'y'>
2036 : ! = sqrt( mixt_frac * ( 1 - mixt_frac ) ) * sqrt( F_w ) * sqrt( <w'^2> )
2037 : ! * ( sqrt( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd )
2038 : ! - sqrt( coef_sigma_x_2_sqd * coef_sigma_y_2_sqd ) )
2039 : ! / ( mixt_frac * sqrt( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd )
2040 : ! + ( 1 - mixt_frac )
2041 : ! * sqrt( coef_sigma_x_2_sqd * coef_sigma_y_2_sqd ) )
2042 : ! * <x'y'>
2043 : ! + sqrt( mixt_frac * ( 1 - mixt_frac ) ) * sqrt( F_w ) * sqrt( <w'^2> )
2044 : ! * sqrt( F_x ) * sqrt( <x'^2> ) * sgn( <w'x'> )
2045 : ! * sqrt( F_y ) * sqrt( <y'^2> ) * sgn( <w'y'> )
2046 : ! * ( ( 1 - mixt_frac ) / mixt_frac - mixt_frac / ( 1 - mixt_frac )
2047 : ! - ( sqrt( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd )
2048 : ! - sqrt( coef_sigma_x_2_sqd * coef_sigma_y_2_sqd ) )
2049 : ! / ( mixt_frac * sqrt( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd )
2050 : ! + ( 1 - mixt_frac )
2051 : ! * sqrt( coef_sigma_x_2_sqd * coef_sigma_y_2_sqd ) )
2052 : ! - ( sqrt( coef_sigma_w_1_sqd * coef_sigma_y_1_sqd )
2053 : ! - sqrt( coef_sigma_w_2_sqd * coef_sigma_y_2_sqd ) )
2054 : ! / ( mixt_frac * sqrt( coef_sigma_w_1_sqd * coef_sigma_y_1_sqd )
2055 : ! + ( 1 - mixt_frac )
2056 : ! * sqrt( coef_sigma_w_2_sqd * coef_sigma_y_2_sqd ) )
2057 : ! - ( sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd )
2058 : ! - sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) )
2059 : ! / ( mixt_frac * sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd )
2060 : ! + ( 1 - mixt_frac )
2061 : ! * sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) ) )
2062 : ! + sqrt( mixt_frac * ( 1 - mixt_frac ) )
2063 : ! * sqrt( F_x ) * sqrt( <x'^2> ) * sgn( <w'x'> )
2064 : ! * ( sqrt( coef_sigma_w_1_sqd * coef_sigma_y_1_sqd )
2065 : ! - sqrt( coef_sigma_w_2_sqd * coef_sigma_y_2_sqd ) )
2066 : ! / ( mixt_frac * sqrt( coef_sigma_w_1_sqd * coef_sigma_y_1_sqd )
2067 : ! + ( 1 - mixt_frac )
2068 : ! * sqrt( coef_sigma_w_2_sqd * coef_sigma_y_2_sqd ) )
2069 : ! * <w'y'>
2070 : ! + sqrt( mixt_frac * ( 1 - mixt_frac ) )
2071 : ! * sqrt( F_y ) * sqrt( <y'^2> ) * sgn( <w'y'> )
2072 : ! * ( sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd )
2073 : ! - sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) )
2074 : ! / ( mixt_frac * sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd )
2075 : ! + ( 1 - mixt_frac )
2076 : ! * sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) )
2077 : ! * <w'x'>.
2078 : !
2079 : ! This equation is of the form:
2080 : !
2081 : ! <w'^2 x'> = coef_wp2xp_implicit * <w'x'> + term_wp2xp_explicit;
2082 : !
2083 : ! where:
2084 : !
2085 : ! coef_wpxpyp_implicit
2086 : ! = sqrt( mixt_frac * ( 1 - mixt_frac ) ) * sqrt( F_w ) * sqrt( <w'^2> )
2087 : ! * ( sqrt( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd )
2088 : ! - sqrt( coef_sigma_x_2_sqd * coef_sigma_y_2_sqd ) )
2089 : ! / ( mixt_frac * sqrt( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd )
2090 : ! + ( 1 - mixt_frac )
2091 : ! * sqrt( coef_sigma_x_2_sqd * coef_sigma_y_2_sqd ) ); and
2092 : !
2093 : ! term_wpxpyp_explicit
2094 : ! = sqrt( mixt_frac * ( 1 - mixt_frac ) ) * sqrt( F_w ) * sqrt( <w'^2> )
2095 : ! * sqrt( F_x ) * sqrt( <x'^2> ) * sgn( <w'x'> )
2096 : ! * sqrt( F_y ) * sqrt( <y'^2> ) * sgn( <w'y'> )
2097 : ! * ( ( 1 - mixt_frac ) / mixt_frac - mixt_frac / ( 1 - mixt_frac )
2098 : ! - ( sqrt( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd )
2099 : ! - sqrt( coef_sigma_x_2_sqd * coef_sigma_y_2_sqd ) )
2100 : ! / ( mixt_frac * sqrt( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd )
2101 : ! + ( 1 - mixt_frac )
2102 : ! * sqrt( coef_sigma_x_2_sqd * coef_sigma_y_2_sqd ) )
2103 : ! - ( sqrt( coef_sigma_w_1_sqd * coef_sigma_y_1_sqd )
2104 : ! - sqrt( coef_sigma_w_2_sqd * coef_sigma_y_2_sqd ) )
2105 : ! / ( mixt_frac * sqrt( coef_sigma_w_1_sqd * coef_sigma_y_1_sqd )
2106 : ! + ( 1 - mixt_frac )
2107 : ! * sqrt( coef_sigma_w_2_sqd * coef_sigma_y_2_sqd ) )
2108 : ! - ( sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd )
2109 : ! - sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) )
2110 : ! / ( mixt_frac * sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd )
2111 : ! + ( 1 - mixt_frac )
2112 : ! * sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) ) )
2113 : ! + sqrt( mixt_frac * ( 1 - mixt_frac ) )
2114 : ! * sqrt( F_x ) * sqrt( <x'^2> ) * sgn( <w'x'> )
2115 : ! * ( sqrt( coef_sigma_w_1_sqd * coef_sigma_y_1_sqd )
2116 : ! - sqrt( coef_sigma_w_2_sqd * coef_sigma_y_2_sqd ) )
2117 : ! / ( mixt_frac * sqrt( coef_sigma_w_1_sqd * coef_sigma_y_1_sqd )
2118 : ! + ( 1 - mixt_frac )
2119 : ! * sqrt( coef_sigma_w_2_sqd * coef_sigma_y_2_sqd ) )
2120 : ! * <w'y'>
2121 : ! + sqrt( mixt_frac * ( 1 - mixt_frac ) )
2122 : ! * sqrt( F_y ) * sqrt( <y'^2> ) * sgn( <w'y'> )
2123 : ! * ( sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd )
2124 : ! - sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) )
2125 : ! / ( mixt_frac * sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd )
2126 : ! + ( 1 - mixt_frac )
2127 : ! * sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) )
2128 : ! * <w'x'>.
2129 : !
2130 : ! There are also special cases for the above equations.
2131 :
2132 : ! References:
2133 : !-----------------------------------------------------------------------
2134 :
2135 : use constants_clubb, only: &
2136 : one, & ! Variable(s)
2137 : zero
2138 :
2139 : use clubb_precision, only: &
2140 : core_rknd ! Procedure(s)
2141 :
2142 : implicit none
2143 :
2144 : integer, intent(in) :: &
2145 : nz
2146 :
2147 : ! Input Variables
2148 : real ( kind = core_rknd ), dimension(nz), intent(in) :: &
2149 : wp2, & ! Variance of w (overall) [m^2/s^2]
2150 : xp2, & ! Variance of x (overall) [(x units)^2]
2151 : yp2, & ! Variance of y (overall) [(y units)^2]
2152 : wpxp, & ! Covariance of w and x (overall) [m/s(x units)]
2153 : wpyp, & ! Covariance of w and y (overall) [m/s(y units)]
2154 : sgn_wpxp, & ! Sign of the covariance of w and x [-]
2155 : sgn_wpyp, & ! Sign of the covariance of w and y [-]
2156 : mixt_frac, & ! Mixture fraction [-]
2157 : F_w, & ! Parameter: spread of the PDF comp. means of w [-]
2158 : F_x, & ! Parameter: spread of the PDF comp. means of x [-]
2159 : F_y, & ! Parameter: spread of the PDF comp. means of y [-]
2160 : coef_sigma_w_1_sqd, & ! sigma_w_1^2 = coef_sigma_w_1_sqd * <w'^2> [-]
2161 : coef_sigma_w_2_sqd, & ! sigma_w_2^2 = coef_sigma_w_2_sqd * <w'^2> [-]
2162 : coef_sigma_x_1_sqd, & ! sigma_x_1^2 = coef_sigma_x_1_sqd * <x'^2> [-]
2163 : coef_sigma_x_2_sqd, & ! sigma_x_2^2 = coef_sigma_x_2_sqd * <x'^2> [-]
2164 : coef_sigma_y_1_sqd, & ! sigma_y_1^2 = coef_sigma_y_1_sqd * <y'^2> [-]
2165 : coef_sigma_y_2_sqd ! sigma_y_2^2 = coef_sigma_y_2_sqd * <y'^2> [-]
2166 :
2167 : ! Output Variables
2168 : ! Coefs.: <w'x'y'> = coef_wpxpyp_implicit * <x'y'> + term_wpxpyp_explicit
2169 : real ( kind = core_rknd ), dimension(nz), intent(out) :: &
2170 : coef_wpxpyp_implicit, & ! Coefficient that is multiplied by <x'y'> [m/s]
2171 : term_wpxpyp_explicit ! Term that is on the RHS [m/s(x units)(y units)]
2172 :
2173 : ! Local Variables
2174 : real ( kind = core_rknd ), dimension(nz) :: &
2175 0 : coefs_factor_wx, & ! Factor involving coef_sigma_... w and x coefs [-]
2176 0 : coefs_factor_wy, & ! Factor involving coef_sigma_... w and y coefs [-]
2177 0 : coefs_factor_xy ! Factor involving coef_sigma_... x and y coefs [-]
2178 :
2179 :
2180 : ! Calculate coefs_factor_wx.
2181 : where ( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd > zero &
2182 0 : .or. coef_sigma_w_2_sqd * coef_sigma_x_2_sqd > zero )
2183 :
2184 : ! coefs_factor_wx
2185 : ! = ( sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd )
2186 : ! - sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) )
2187 : ! / ( mixt_frac * sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd )
2188 : ! + ( 1 - mixt_frac )
2189 : ! * sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) )
2190 : coefs_factor_wx &
2191 : = ( sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd ) &
2192 : - sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) ) &
2193 : / ( mixt_frac * sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd ) &
2194 : + ( one - mixt_frac ) &
2195 : * sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) )
2196 :
2197 : elsewhere ! coef_sigma_w_1_sqd * coef_sigma_x_1_sqd = 0
2198 : ! and coef_sigma_w_2_sqd * coef_sigma_x_2_sqd = 0
2199 :
2200 : ! When coef_sigma_w_1_sqd * coef_sigma_x_1_sqd = 0 and
2201 : ! coef_sigma_w_2_sqd * coef_sigma_x_2_sqd = 0, the value of
2202 : ! coefs_factor_wx is undefined. However, setting coefs_factor_wx to a
2203 : ! value of 0 in this scenario allows for the use of general form
2204 : ! equations below for coef_wpxpyp_implicit and term_wpxpyp_explicit.
2205 : coefs_factor_wx = zero
2206 :
2207 : endwhere
2208 :
2209 : ! Calculate coefs_factor_wy.
2210 : where ( coef_sigma_w_1_sqd * coef_sigma_y_1_sqd > zero &
2211 0 : .or. coef_sigma_w_2_sqd * coef_sigma_y_2_sqd > zero )
2212 :
2213 : ! coefs_factor_wy
2214 : ! = ( sqrt( coef_sigma_w_1_sqd * coef_sigma_y_1_sqd )
2215 : ! - sqrt( coef_sigma_w_2_sqd * coef_sigma_y_2_sqd ) )
2216 : ! / ( mixt_frac * sqrt( coef_sigma_w_1_sqd * coef_sigma_y_1_sqd )
2217 : ! + ( 1 - mixt_frac )
2218 : ! * sqrt( coef_sigma_w_2_sqd * coef_sigma_y_2_sqd ) )
2219 : coefs_factor_wy &
2220 : = ( sqrt( coef_sigma_w_1_sqd * coef_sigma_y_1_sqd ) &
2221 : - sqrt( coef_sigma_w_2_sqd * coef_sigma_y_2_sqd ) ) &
2222 : / ( mixt_frac * sqrt( coef_sigma_w_1_sqd * coef_sigma_y_1_sqd ) &
2223 : + ( one - mixt_frac ) &
2224 : * sqrt( coef_sigma_w_2_sqd * coef_sigma_y_2_sqd ) )
2225 :
2226 : elsewhere ! coef_sigma_w_1_sqd * coef_sigma_y_1_sqd = 0
2227 : ! and coef_sigma_w_2_sqd * coef_sigma_y_2_sqd = 0
2228 :
2229 : ! When coef_sigma_w_1_sqd * coef_sigma_y_1_sqd = 0 and
2230 : ! coef_sigma_w_2_sqd * coef_sigma_y_2_sqd = 0, the value of
2231 : ! coefs_factor_wy is undefined. However, setting coefs_factor_wy to a
2232 : ! value of 0 in this scenario allows for the use of general form
2233 : ! equations below for coef_wpxpyp_implicit and term_wpxpyp_explicit.
2234 : coefs_factor_wy = zero
2235 :
2236 : endwhere
2237 :
2238 : ! Calculate coefs_factor_xy.
2239 : where ( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd > zero &
2240 0 : .or. coef_sigma_x_2_sqd * coef_sigma_y_2_sqd > zero )
2241 :
2242 : ! coefs_factor_xy
2243 : ! = ( sqrt( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd )
2244 : ! - sqrt( coef_sigma_x_2_sqd * coef_sigma_y_2_sqd ) )
2245 : ! / ( mixt_frac * sqrt( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd )
2246 : ! + ( 1 - mixt_frac )
2247 : ! * sqrt( coef_sigma_x_2_sqd * coef_sigma_y_2_sqd ) )
2248 : coefs_factor_xy &
2249 : = ( sqrt( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd ) &
2250 : - sqrt( coef_sigma_x_2_sqd * coef_sigma_y_2_sqd ) ) &
2251 : / ( mixt_frac * sqrt( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd ) &
2252 : + ( one - mixt_frac ) &
2253 : * sqrt( coef_sigma_x_2_sqd * coef_sigma_y_2_sqd ) )
2254 :
2255 : elsewhere ! coef_sigma_x_1_sqd * coef_sigma_y_1_sqd = 0
2256 : ! and coef_sigma_x_2_sqd * coef_sigma_y_2_sqd = 0
2257 :
2258 : ! When coef_sigma_x_1_sqd * coef_sigma_y_1_sqd = 0 and
2259 : ! coef_sigma_x_2_sqd * coef_sigma_y_2_sqd = 0, the value of
2260 : ! coefs_factor_xy is undefined. However, setting coefs_factor_xy to a
2261 : ! value of 0 in this scenario allows for the use of general form
2262 : ! equations below for coef_wpxpyp_implicit and term_wpxpyp_explicit.
2263 : coefs_factor_xy = zero
2264 :
2265 : endwhere
2266 :
2267 :
2268 : ! Calculate coef_wpxpyp_implicit and term_wpxpyp_explicit.
2269 : where ( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd > zero &
2270 0 : .or. coef_sigma_x_2_sqd * coef_sigma_y_2_sqd > zero )
2271 :
2272 : coef_wpxpyp_implicit &
2273 : = sqrt( mixt_frac * ( one - mixt_frac ) ) &
2274 : * sqrt( F_w ) * sqrt( wp2 ) * coefs_factor_xy
2275 :
2276 : term_wpxpyp_explicit &
2277 : = sqrt( mixt_frac * ( one - mixt_frac ) ) * sqrt( F_w ) * sqrt( wp2 ) &
2278 : * sqrt( F_x ) * sqrt( xp2 ) * sgn_wpxp &
2279 : * sqrt( F_y ) * sqrt( yp2 ) * sgn_wpyp &
2280 : * ( ( one - mixt_frac ) / mixt_frac - mixt_frac / ( one - mixt_frac ) &
2281 : - coefs_factor_xy - coefs_factor_wy - coefs_factor_wx ) &
2282 : + sqrt( mixt_frac * ( one - mixt_frac ) ) &
2283 : * sqrt( F_x ) * sqrt( xp2 ) * sgn_wpxp * coefs_factor_wy * wpyp &
2284 : + sqrt( mixt_frac * ( one - mixt_frac ) ) &
2285 : * sqrt( F_y ) * sqrt( yp2 ) * sgn_wpyp * coefs_factor_wx * wpxp
2286 :
2287 : elsewhere ! coef_sigma_x_1_sqd * coef_sigma_y_1_sqd = 0
2288 : ! and coef_sigma_x_2_sqd * coef_sigma_y_2_sqd = 0
2289 :
2290 : coef_wpxpyp_implicit &
2291 : = sqrt( mixt_frac * ( one - mixt_frac ) ) * sqrt( F_w ) * sqrt( wp2 ) &
2292 : * ( ( one - mixt_frac ) / mixt_frac - mixt_frac / ( one - mixt_frac ) &
2293 : - coefs_factor_wy - coefs_factor_wx )
2294 :
2295 : term_wpxpyp_explicit &
2296 : = sqrt( mixt_frac * ( one - mixt_frac ) ) &
2297 : * sqrt( F_x ) * sqrt( xp2 ) * sgn_wpxp * coefs_factor_wy * wpyp &
2298 : + sqrt( mixt_frac * ( one - mixt_frac ) ) &
2299 : * sqrt( F_y ) * sqrt( yp2 ) * sgn_wpyp * coefs_factor_wx * wpxp
2300 :
2301 : endwhere
2302 :
2303 :
2304 0 : return
2305 :
2306 : end subroutine calc_coefs_wpxpyp_semiimpl
2307 :
2308 : !=============================================================================
2309 :
2310 : end module new_pdf
|
