Line data Source code
1 : ! $Id$
2 : !===============================================================================
3 : module new_tsdadg_pdf
4 :
5 : ! Description:
6 : ! The new trivariate, skewness-dependent, analytic double Gaussian (TSDADG)
7 : ! PDF.
8 :
9 : implicit none
10 :
11 : public :: tsdadg_pdf_driver, & ! Procedure(s)
12 : calc_setter_parameters, &
13 : calc_L_x_Skx_fnc
14 :
15 : private :: calc_respnder_parameters ! Procedure(s)
16 :
17 : private ! default scope
18 :
19 : contains
20 :
21 : !=============================================================================
22 0 : subroutine tsdadg_pdf_driver( nz, wm, rtm, thlm, wp2, rtp2, thlp2, & ! In
23 0 : Skw, Skrt, Skthl, wprtp, wpthlp, & ! In
24 0 : mu_w_1, mu_w_2, & ! Out
25 0 : mu_rt_1, mu_rt_2, & ! Out
26 0 : mu_thl_1, mu_thl_2, & ! Out
27 0 : sigma_w_1_sqd, sigma_w_2_sqd, & ! Out
28 0 : sigma_rt_1_sqd, sigma_rt_2_sqd, & ! Out
29 0 : sigma_thl_1_sqd, sigma_thl_2_sqd, & ! Out
30 0 : mixt_frac ) ! Out
31 :
32 :
33 : ! Description:
34 : ! Selects which variable is used to set the mixture fraction for the PDF
35 : ! ("the setter") and which variables are handled after that mixture fraction
36 : ! has been set ("the responders"). Traditionally, w has been used to set
37 : ! the PDF. However, here, the variable with the greatest magnitude of
38 : ! skewness is used to set the PDF.
39 :
40 : ! References:
41 : !-----------------------------------------------------------------------
42 :
43 : use constants_clubb, only: &
44 : one, & ! Variable(s)
45 : zero, &
46 : fstderr
47 :
48 : use clubb_precision, only: &
49 : core_rknd ! Variable(s)
50 :
51 : implicit none
52 :
53 : integer, intent(in) :: &
54 : nz
55 :
56 : ! Input Variables
57 : real( kind = core_rknd ), dimension(nz), intent(in) :: &
58 : wm, & ! Mean of w (overall) [m/s]
59 : rtm, & ! Mean of rt (overall) [kg/kg]
60 : thlm, & ! Mean of thl (overall) [K]
61 : wp2, & ! Variance of w (overall) [m^2/s^2]
62 : rtp2, & ! Variance of rt (overall) [kg^2/kg^2]
63 : thlp2, & ! Variance of thl (overall) [K^2]
64 : Skw, & ! Skewness of w (overall) [-]
65 : Skrt, & ! Skewness of rt (overall) [-]
66 : Skthl, & ! Skewness of thl (overall) [-]
67 : wprtp, & ! Covariance of w and rt (overall) [(m/s)kg/kg]
68 : wpthlp ! Covariance of w and thl (overall) [(m/s)K]
69 :
70 : ! Output Variables
71 : real( kind = core_rknd ), dimension(nz), intent(out) :: &
72 : mu_w_1, & ! Mean of w (1st PDF component) [m/s]
73 : mu_w_2, & ! Mean of w (2nd PDF component) [m/s]
74 : mu_rt_1, & ! Mean of rt (1st PDF component) [kg/kg]
75 : mu_rt_2, & ! Mean of rt (2nd PDF component) [kg/kg]
76 : mu_thl_1, & ! Mean of thl (1st PDF component) [K]
77 : mu_thl_2, & ! Mean of thl (2nd PDF component) [K]
78 : sigma_w_1_sqd, & ! Variance of w (1st PDF component) [m^2/s^2]
79 : sigma_w_2_sqd, & ! Variance of w (2nd PDF component) [m^2/s^2]
80 : sigma_rt_1_sqd, & ! Variance of rt (1st PDF component) [kg^2/kg^2]
81 : sigma_rt_2_sqd, & ! Variance of rt (2nd PDF component) [kg^2/kg^2]
82 : sigma_thl_1_sqd, & ! Variance of thl (1st PDF component) [K^2]
83 : sigma_thl_2_sqd, & ! Variance of thl (2nd PDF component) [K^2]
84 : mixt_frac ! Mixture fraction [-]
85 :
86 : ! Local Variables
87 : real( kind = core_rknd ), dimension(nz) :: &
88 0 : big_L_w_1, & ! Parameter for the 1st PDF comp. mean of w [-]
89 0 : big_L_w_2, & ! Parameter for the 2nd PDF comp. mean of w (setter) [-]
90 0 : big_L_rt_1, & ! Parameter for the 1st PDF comp. mean of rt [-]
91 0 : big_L_rt_2, & ! Parameter for the 2nd PDF comp. mean of rt (setter) [-]
92 0 : big_L_thl_1, & ! Parameter for the 1st PDF comp. mean of thl [-]
93 0 : big_L_thl_2 ! Parameter for the 2nd PDF comp. mean of thl (setter) [-]
94 :
95 : real( kind = core_rknd ) :: &
96 : small_l_w_1, & ! Param. for the 1st PDF comp. mean of w [-]
97 : small_l_w_2, & ! Param. for the 2nd PDF comp. mean of w (setter) [-]
98 : small_l_rt_1, & ! Param. for the 1st PDF comp. mean of rt [-]
99 : small_l_rt_2, & ! Param. for the 2nd PDF comp. mean of rt (setter) [-]
100 : small_l_thl_1, & ! Param. for the 1st PDF comp. mean of thl [-]
101 : small_l_thl_2 ! Param. for the 2nd PDF comp. mean of thl (setter) [-]
102 :
103 : real( kind = core_rknd ), dimension(nz) :: &
104 0 : sgn_wprtp, & ! Sign of the covariance of w and rt (overall) [-]
105 0 : sgn_wpthlp, & ! Sign of the covariance of w and thl (overall) [-]
106 0 : sgn_wp2 ! Sign of the variance of w (overall); always positive [-]
107 :
108 : real( kind = core_rknd ), dimension(nz) :: &
109 0 : coef_sigma_w_1_sqd, & ! sigma_w_1^2 = coef_sigma_w_1_sqd * <w'^2> [-]
110 0 : coef_sigma_w_2_sqd, & ! sigma_w_2^2 = coef_sigma_w_2_sqd * <w'^2> [-]
111 0 : coef_sigma_rt_1_sqd, & ! sigma_rt_1^2 = coef_sigma_rt_1_sqd * <rt'^2> [-]
112 0 : coef_sigma_rt_2_sqd, & ! sigma_rt_2^2 = coef_sigma_rt_2_sqd * <rt'^2> [-]
113 0 : coef_sigma_thl_1_sqd, & ! sigma_thl_1^2=coef_sigma_thl_1_sqd*<thl'^2> [-]
114 0 : coef_sigma_thl_2_sqd ! sigma_thl_2^2=coef_sigma_thl_2_sqd*<thl'^2> [-]
115 :
116 : integer :: k ! Vertical level loop index
117 :
118 :
119 : ! Calculate sgn( <w'rt'> ).
120 0 : where ( wprtp >= zero )
121 : sgn_wprtp = one
122 : elsewhere ! wprtp < 0
123 : sgn_wprtp = -one
124 : endwhere ! wprtp >= 0
125 :
126 : ! Calculate sgn( <w'thl'> ).
127 0 : where ( wpthlp >= zero )
128 : sgn_wpthlp = one
129 : elsewhere ! wpthlp < 0
130 : sgn_wpthlp = -one
131 : endwhere ! wpthlp >= 0
132 :
133 : ! The sign of the variance of w is always positive.
134 0 : sgn_wp2 = one
135 :
136 0 : small_l_w_1 = 0.75_core_rknd
137 0 : small_l_w_2 = 0.5_core_rknd
138 0 : small_l_rt_1 = 0.75_core_rknd
139 0 : small_l_rt_2 = 0.5_core_rknd
140 0 : small_l_thl_1 = 0.75_core_rknd
141 0 : small_l_thl_2 = 0.5_core_rknd
142 :
143 0 : do k = 1, nz, 1
144 :
145 0 : call calc_L_x_Skx_fnc( Skw(k), sgn_wp2(k), & ! In
146 : small_l_w_1, small_l_w_2, & ! In
147 0 : big_L_w_1(k), big_L_w_2(k) ) ! Out
148 :
149 : call calc_L_x_Skx_fnc( Skrt(k), sgn_wprtp(k), & ! In
150 : small_l_rt_1, small_l_rt_2, & ! In
151 : big_L_rt_1(k), big_L_rt_2(k) ) ! Out
152 :
153 : call calc_L_x_Skx_fnc( Skthl(k), sgn_wpthlp(k), & ! In
154 : small_l_thl_1, small_l_thl_2, & ! In
155 : big_L_thl_1(k), big_L_thl_2(k) ) ! Out
156 :
157 :
158 : ! The variable with the greatest magnitude of skewness will be the setter
159 : ! variable and the other variables will be responder variables.
160 : if ( abs( Skw(k) ) >= abs( Skrt(k) ) &
161 0 : .and. abs( Skw(k) ) >= abs( Skthl(k) ) ) then
162 :
163 : ! The variable w has the greatest magnitude of skewness.
164 :
165 : call calc_setter_parameters( wm(k), wp2(k), & ! In
166 : Skw(k), sgn_wp2(k), & ! In
167 : big_L_w_1(k), big_L_w_2(k), & ! In
168 : mu_w_1(k), mu_w_2(k), & ! Out
169 : sigma_w_1_sqd(k), & ! Out
170 : sigma_w_2_sqd(k), mixt_frac(k), & ! Out
171 : coef_sigma_w_1_sqd(k), & ! Out
172 : coef_sigma_w_2_sqd(k) ) ! Out
173 :
174 : call calc_respnder_parameters( rtm(k), rtp2(k), & ! In
175 : Skrt(k), sgn_wprtp(k), & ! In
176 : mixt_frac(k), big_L_rt_1(k), & ! In
177 : mu_rt_1(k), mu_rt_2(k), & ! Out
178 : sigma_rt_1_sqd(k), & ! Out
179 : sigma_rt_2_sqd(k), & ! Out
180 : coef_sigma_rt_1_sqd(k), & ! Out
181 : coef_sigma_rt_2_sqd(k) ) ! Out
182 :
183 : call calc_respnder_parameters( thlm(k), thlp2(k), & ! In
184 : Skthl(k), sgn_wpthlp(k), & ! In
185 : mixt_frac(k), big_L_thl_1(k), & ! In
186 : mu_thl_1(k), mu_thl_2(k), & ! Out
187 : sigma_thl_1_sqd(k), & ! Out
188 : sigma_thl_2_sqd(k), & ! Out
189 : coef_sigma_thl_1_sqd(k), & ! Out
190 : coef_sigma_thl_2_sqd(k) ) ! Out
191 :
192 :
193 : elseif ( abs( Skrt(k) ) > abs( Skw(k) ) &
194 0 : .and. abs( Skrt(k) ) >= abs( Skthl(k) ) ) then
195 :
196 : ! The variable rt has the greatest magnitude of skewness.
197 :
198 : call calc_setter_parameters( rtm(k), rtp2(k), & ! In
199 : Skrt(k), sgn_wprtp(k), & ! In
200 : big_L_rt_1(k), big_L_rt_2(k), & ! In
201 : mu_rt_1(k), mu_rt_2(k), & ! Out
202 : sigma_rt_1_sqd(k), & ! Out
203 : sigma_rt_2_sqd(k), mixt_frac(k), & ! Out
204 : coef_sigma_rt_1_sqd(k), & ! Out
205 : coef_sigma_rt_2_sqd(k) ) ! Out
206 :
207 : call calc_respnder_parameters( wm(k), wp2(k), & ! In
208 : Skw(k), sgn_wp2(k), & ! In
209 : mixt_frac(k), big_L_w_1(k), & ! In
210 : mu_w_1(k), mu_w_2(k), & ! Out
211 : sigma_w_1_sqd(k), & ! Out
212 : sigma_w_2_sqd(k), & ! Out
213 : coef_sigma_w_1_sqd(k), & ! Out
214 : coef_sigma_w_2_sqd(k) ) ! Out
215 :
216 : call calc_respnder_parameters( thlm(k), thlp2(k), & ! In
217 : Skthl(k), sgn_wpthlp(k), & ! In
218 : mixt_frac(k), big_L_thl_1(k), & ! In
219 : mu_thl_1(k), mu_thl_2(k), & ! Out
220 : sigma_thl_1_sqd(k), & ! Out
221 : sigma_thl_2_sqd(k), & ! Out
222 : coef_sigma_thl_1_sqd(k), & ! Out
223 : coef_sigma_thl_2_sqd(k) ) ! Out
224 :
225 :
226 : else ! abs( Skthl ) > abs( Skw ) .and. abs( Skthl ) > abs( Skrt )
227 :
228 : ! The variable thl has the greatest magnitude of skewness.
229 :
230 : call calc_setter_parameters( thlm(k), thlp2(k), & ! In
231 : Skthl(k), sgn_wpthlp(k), & ! In
232 : big_L_thl_1(k), big_L_thl_2(k), & ! In
233 : mu_thl_1(k), mu_thl_2(k), & ! Out
234 : sigma_thl_1_sqd(k), & ! Out
235 : sigma_thl_2_sqd(k), mixt_frac(k), & ! Out
236 : coef_sigma_thl_1_sqd(k), & ! Out
237 : coef_sigma_thl_2_sqd(k) ) ! Out
238 :
239 : call calc_respnder_parameters( wm(k), wp2(k), & ! In
240 : Skw(k), sgn_wp2(k), & ! In
241 : mixt_frac(k), big_L_w_1(k), & ! In
242 : mu_w_1(k), mu_w_2(k), & ! Out
243 : sigma_w_1_sqd(k), & ! Out
244 : sigma_w_2_sqd(k), & ! Out
245 : coef_sigma_w_1_sqd(k), & ! Out
246 : coef_sigma_w_2_sqd(k) ) ! Out
247 :
248 : call calc_respnder_parameters( rtm(k), rtp2(k), & ! In
249 : Skrt(k), sgn_wprtp(k), & ! In
250 : mixt_frac(k), big_L_rt_1(k), & ! In
251 : mu_rt_1(k), mu_rt_2(k), & ! Out
252 : sigma_rt_1_sqd(k), & ! Out
253 : sigma_rt_2_sqd(k), & ! Out
254 : coef_sigma_rt_1_sqd(k), & ! Out
255 : coef_sigma_rt_2_sqd(k) ) ! Out
256 :
257 :
258 : endif ! Find variable with the greatest magnitude of skewness.
259 :
260 :
261 0 : if ( sigma_w_1_sqd(k) < zero ) then
262 : write(fstderr,*) "WARNING: New TSDADG PDF. The variance of w in " &
263 : // "the 1st PDF component is negative and is " &
264 0 : // "being clipped to 0."
265 0 : write(fstderr,*) "sigma_w_1^2 (before clipping) = ", sigma_w_1_sqd(k)
266 0 : sigma_w_1_sqd(k) = zero
267 0 : coef_sigma_w_1_sqd(k) = zero
268 : endif ! sigma_w_1_sqd < 0
269 :
270 0 : if ( sigma_w_2_sqd(k) < zero ) then
271 : write(fstderr,*) "WARNING: New TSDADG PDF. The variance of w in " &
272 : // "the 2nd PDF component is negative and is " &
273 0 : // "being clipped to 0."
274 0 : write(fstderr,*) "sigma_w_2^2 (before clipping) = ", sigma_w_2_sqd(k)
275 0 : sigma_w_2_sqd(k) = zero
276 0 : coef_sigma_w_2_sqd(k) = zero
277 : endif ! sigma_w_2_sqd < 0
278 :
279 0 : if ( sigma_rt_1_sqd(k) < zero ) then
280 : write(fstderr,*) "WARNING: New TSDADG PDF. The variance of rt in " &
281 : // "the 1st PDF component is negative and is " &
282 0 : // "being clipped to 0."
283 0 : write(fstderr,*) "sigma_rt_1^2 (before clipping) = ", &
284 0 : sigma_rt_1_sqd(k)
285 0 : sigma_rt_1_sqd(k) = zero
286 0 : coef_sigma_rt_1_sqd(k) = zero
287 : endif ! sigma_rt_1_sqd < 0
288 :
289 0 : if ( sigma_rt_2_sqd(k) < zero ) then
290 : write(fstderr,*) "WARNING: New TSDADG PDF. The variance of rt in " &
291 : // "the 2nd PDF component is negative and is " &
292 0 : // "being clipped to 0."
293 0 : write(fstderr,*) "sigma_rt_2^2 (before clipping) = ", &
294 0 : sigma_rt_2_sqd(k)
295 0 : sigma_rt_2_sqd(k) = zero
296 0 : coef_sigma_rt_2_sqd(k) = zero
297 : endif ! sigma_rt_2_sqd < 0
298 :
299 0 : if ( sigma_thl_1_sqd(k) < zero ) then
300 : write(fstderr,*) "WARNING: New TSDADG PDF. The variance of thl " &
301 : // "in the 1st PDF component is negative and is " &
302 0 : // "being clipped to 0."
303 0 : write(fstderr,*) "sigma_thl_1^2 (before clipping) = ", &
304 0 : sigma_thl_1_sqd(k)
305 0 : sigma_thl_1_sqd(k) = zero
306 0 : coef_sigma_thl_1_sqd(k) = zero
307 : endif ! sigma_thl_1_sqd < 0
308 :
309 0 : if ( sigma_thl_2_sqd(k) < zero ) then
310 : write(fstderr,*) "WARNING: New TSDADG PDF. The variance of thl " &
311 : // "in the 2nd PDF component is negative and is " &
312 0 : // "being clipped to 0."
313 0 : write(fstderr,*) "sigma_thl_2^2 (before clipping) = ", &
314 0 : sigma_thl_2_sqd(k)
315 0 : sigma_thl_2_sqd(k) = zero
316 0 : coef_sigma_thl_2_sqd(k) = zero
317 : endif ! sigma_thl_2_sqd < 0
318 :
319 : enddo ! k = 1, nz, 1
320 :
321 :
322 0 : return
323 :
324 : end subroutine tsdadg_pdf_driver
325 :
326 : !=============================================================================
327 : !
328 : ! DESCRIPTION OF THE METHOD FOR THE VARIABLE THAT SETS THE MIXTURE FRACTION
329 : ! =========================================================================
330 : !
331 : ! There are five PDF parameters that need to be calculated for the setting
332 : ! variable, which are mu_x_1 (the mean of x in the 1st PDF component), mu_x_2
333 : ! (the mean of x in the 2nd PDF component), sigma_x_1 (the standard deviation
334 : ! of x in the 1st PDF component), sigma_x_2 (the standard deviation of x in
335 : ! the 2nd PDF component), and mixt_frac (the mixture fraction). In order to
336 : ! solve for these five parameters, five equations are needed. These five
337 : ! equations are the equations for <x>, <x'^2>, and <x'^3> as found by
338 : ! integrating over the PDF. Additionally, two more equations, which involve
339 : ! tunable parameters L_x_1 and L_x_2, and which are used to help control the
340 : ! spread of the PDF component means in the 1st PDF component and the 2nd PDF
341 : ! component, respectively, are used in this equation set. The five equations
342 : ! are:
343 : !
344 : ! <x> = mixt_frac * mu_x_1 + ( 1 - mixt_frac ) * mu_x_2;
345 : !
346 : ! <x'^2> = mixt_frac * ( ( mu_x_1 - <x> )^2 + sigma_x_1^2 )
347 : ! + ( 1 - mixt_frac ) * ( ( mu_x_2 - <x> )^2 + sigma_x_2^2 );
348 : !
349 : ! <x'^3> = mixt_frac * ( mu_x_1 - <x> )
350 : ! * ( ( mu_x_1 - <x> )^2 + 3 * sigma_x_1^2 )
351 : ! + ( 1 - mixt_frac ) * ( mu_x_2 - <x> )
352 : ! * ( ( mu_x_2 - <x> )^2 + 3 * sigma_x_2^2 );
353 : !
354 : ! mu_x_1 - <x> = L_x_1
355 : ! * sqrt( ( 1 + Skx * sgn( <w'x'> ) / sqrt( 4 + Skx^2 ) )
356 : ! / ( 1 - Skx * sgn( <w'x'> ) / sqrt( 4 + Skx^2 ) ) )
357 : ! * sqrt( <x'^2> ) * sgn( <w'x'> ); and
358 : !
359 : ! mu_x_2 - <x> = -L_x_2
360 : ! * sqrt( ( 1 - Skx * sgn( <w'x'> ) / sqrt( 4 + Skx^2 ) )
361 : ! / ( 1 + Skx * sgn( <w'x'> ) / sqrt( 4 + Skx^2 ) ) )
362 : ! * sqrt( <x'^2> ) * sgn( <w'x'> );
363 : !
364 : ! where 0 <= L_x_1 <= 1, 0 <= L_x_2 <= 1, Skx is the skewness of x, such that
365 : ! Skx = <x'^3> / <x'^2>^(3/2), and sgn( <w'x'> ) is the sign of <w'x'>, such
366 : ! that:
367 : !
368 : ! sgn( <w'x'> ) = | 1, when <w'x'> >= 0;
369 : ! | -1, when <w'x'> < 0.
370 : !
371 : ! The resulting equations for the five PDF parameters are:
372 : !
373 : ! mu_x_1 = <x> + L_x_1
374 : ! * sqrt( ( 1 + Skx * sgn( <w'x'> ) / sqrt( 4 + Skx^2 ) )
375 : ! / ( 1 - Skx * sgn( <w'x'> ) / sqrt( 4 + Skx^2 ) ) )
376 : ! * sqrt( <x'^2> ) * sgn( <w'x'> );
377 : !
378 : ! mu_x_2 = <x> - L_x_2
379 : ! * sqrt( ( 1 - Skx * sgn( <w'x'> ) / sqrt( 4 + Skx^2 ) )
380 : ! / ( 1 + Skx * sgn( <w'x'> ) / sqrt( 4 + Skx^2 ) ) )
381 : ! * sqrt( <x'^2> ) * sgn( <w'x'> );
382 : !
383 : ! mixt_frac = 1 / ( 1 + abs( mu_x_1_nrmlized / mu_x_2_nrmlized ) );
384 : !
385 : ! sigma_x_1 = sqrt( ( 1 - mixt_frac * mu_x_1_nrmlized^2
386 : ! - ( 1 - mixt_frac ) * mu_x_2_nrmlized^2
387 : ! + ( 1 - mixt_frac )
388 : ! * ( Skx / ( 3 * mixt_frac * mu_x_1_nrmlized )
389 : ! - mu_x_1_nrmlized^2 / 3
390 : ! + mu_x_2_nrmlized^2 / 3 ) )
391 : ! * <x'^2> ); and
392 : !
393 : ! sigma_x_2 = sqrt( ( 1 - mixt_frac * mu_x_1_nrmlized^2
394 : ! - ( 1 - mixt_frac ) * mu_x_2_nrmlized^2
395 : ! - mixt_frac
396 : ! * ( Skx / ( 3 * mixt_frac * mu_x_1_nrmlized )
397 : ! - mu_x_1_nrmlized^2 / 3
398 : ! + mu_x_2_nrmlized^2 / 3 ) )
399 : ! * <x'^2> ); where
400 : !
401 : ! mu_x_1_nrmlized = ( mu_x_1 - <x> ) / sqrt( <x'^2> ); and
402 : !
403 : ! mu_x_2_nrmlized = ( mu_x_2 - <x> ) / sqrt( <x'^2> ).
404 : !
405 : !
406 : ! Notes:
407 : !
408 : ! This method does NOT work for all values of L_x_1 and L_x_2 (where
409 : ! 0 <= L_x_1 <= 1 and 0 <= L_x_2 <= 1). Only a subregion of this parameter
410 : ! space produces valid results.
411 : !
412 : ! When both L_x_1 = 0 and L_x_2 = 0, mu_x_1 = mu_x_2 = <x> (which can only
413 : ! happen when Skx = 0). In this scenario, the above equations for mixt_frac,
414 : ! sigma_x_1, and sigma_x_2 are all undefined. In this special case, the
415 : ! distribution reduces to a single Gaussian, so the following values are set:
416 : ! mixt_frac = 1/2, sigma_x_1 = sqrt( <x'^2> ), and sigma_x_2 = sqrt( <x'^2> ).
417 : !
418 : !
419 : ! Tunable parameters:
420 : !
421 : ! The parameter L_x_1 controls the 1st PDF component mean while L_x_2 controls
422 : ! the 2nd PDF component mean. The equations involving the tunable parameters
423 : ! L_x_1 and L_x_2 (the mu_x_1 and mu_x_2 equations) are based on the values of
424 : ! mu_x_1 and mu_x_2 when sigma_x_1 = sigma_x_2 = 0. In this scenario, the
425 : ! equation for mixture fraction reduces to:
426 : !
427 : ! mixt_frac = (1/2) * ( 1 +/- Skx / sqrt( 4 + Skx^2 ) ).
428 : !
429 : ! The +/- is dependent on the sign of ( mu_x_1 - <x> ) vs. ( mu_x_2 - <x> ).
430 : ! This is dependent on sgn( <w'x'> ), and the mixture fraction equation is
431 : ! written as:
432 : !
433 : ! mixt_frac = (1/2) * ( 1 - Skx * sgn( <w'x'> ) / sqrt( 4 + Skx^2 ) ).
434 : !
435 : ! Meanwhile, the equation for 1 - mixt_frac is:
436 : !
437 : ! 1 - mixt_frac = (1/2) * ( 1 + Skx * sgn( <w'x'> ) / sqrt( 4 + Skx^2 ) ).
438 : !
439 : ! When sigma_x_1 = sigma_x_2 = 0, the equations for mu_x_1 and mu_x_2 are:
440 : !
441 : ! mu_x_1 = <x> + sqrt( ( 1 - mixt_frac ) / mixt_frac ) * sqrt( <x'^2> )
442 : ! * sgn( <w'x'> ); and
443 : !
444 : ! mu_x_2 = <x> - sqrt( mixt_frac / ( 1 - mixt_frac ) ) * sqrt( <x'^2> )
445 : ! * sgn( <w'x'> ).
446 : !
447 : ! Substituting the equations for mixt_frac and 1 - mixt_frac into the
448 : ! equations for mu_x_1 and mu_x_2 (when sigma_x_1 = sigma_x_2 = 0), the
449 : ! equations for mu_x_1 and mu_x_2 become:
450 : !
451 : ! mu_x_1 = <x> + sqrt( ( 1 + Skx * sgn( <w'x'> ) / sqrt( 4 + Skx^2 ) )
452 : ! / ( 1 - Skx * sgn( <w'x'> ) / sqrt( 4 + Skx^2 ) ) )
453 : ! * sqrt( <x'^2> ) * sgn( <w'x'> ); and
454 : !
455 : ! mu_x_2 = <x> - sqrt( ( 1 - Skx * sgn( <w'x'> ) / sqrt( 4 + Skx^2 ) )
456 : ! / ( 1 + Skx * sgn( <w'x'> ) / sqrt( 4 + Skx^2 ) ) )
457 : ! * sqrt( <x'^2> ) * sgn( <w'x'> ).
458 : !
459 : ! These equations represent the maximum deviation of mu_x_1 and mu_x_2 from
460 : ! the overall mean, <x>. The range of parameters of L_x_i is 0 <= L_x_i <= 1.
461 : ! When L_x_1 = L_x_2 = 0, the value of mu_x_1 = mu_x_2 = <x> (and the
462 : ! distribution becomes a single Gaussian). When L_x_i = 1, the value of
463 : ! mu_x_i - <x> is at its maximum magnitude.
464 : !
465 : ! The values of L_x_1 and L_x_2 are also calculated by skewness functions.
466 : ! Those functions are:
467 : !
468 : ! L_x_1 = l_x_1 * abs( Skx ) / sqrt( 4 + Skx^2 ); and
469 : ! L_x_2 = l_x_2 * abs( Skx ) / sqrt( 4 + Skx^2 );
470 : !
471 : ! where both l_x_1 and l_x_2 are tunable parameters.
472 : !
473 : ! As previously stated, this method does not work for all combinations of
474 : ! L_x_1 and L_x_2, but rather only for a subregion of parameter space. This
475 : ! applies to l_x_1 and l_x_2, as well. The conditions on l_x_1 and l_x_2 are:
476 : !
477 : ! 2/3 < l_x_1 < 1 and 0 < l_x_2 < 1; when Skx * sgn( <w'x'> ) >= 0; and
478 : ! 0 < l_x_1 < 1 and 2/3 < l_x_2 < 1; when Skx * sgn( <w'x'> ) < 0.
479 : !
480 : ! The condition that l_x_1 > 2/3 prevents a negative PDF component variance
481 : ! when Skx = 0.
482 : !
483 : !
484 : ! Equations for PDF component standard deviations:
485 : !
486 : ! The equations for the PDF component standard deviations can also be written
487 : ! as:
488 : !
489 : ! sigma_x_1 = sqrt( coef_sigma_x_1_sqd * <x'^2> ); and
490 : !
491 : ! sigma_x_2 = sqrt( coef_sigma_x_2_sqd * <x'^2> ); where
492 : !
493 : ! coef_sigma_x_1_sqd = 1 - mixt_frac * mu_x_1_nrmlized^2
494 : ! - ( 1 - mixt_frac ) * mu_x_2_nrmlized^2
495 : ! + ( 1 - mixt_frac )
496 : ! * ( Skx / ( 3 * mixt_frac * mu_x_1_nrmlized )
497 : ! - mu_x_1_nrmlized^2 / 3
498 : ! + mu_x_2_nrmlized^2 / 3 ); and
499 : !
500 : ! coef_sigma_x_2_sqd = 1 - mixt_frac * mu_x_1_nrmlized^2
501 : ! - ( 1 - mixt_frac ) * mu_x_2_nrmlized^2
502 : ! - mixt_frac
503 : ! * ( Skx / ( 3 * mixt_frac * mu_x_1_nrmlized )
504 : ! - mu_x_1_nrmlized^2 / 3
505 : ! + mu_x_2_nrmlized^2 / 3 ).
506 : !
507 : ! The above equations can be substituted into an equation for a variable that
508 : ! has been derived by integrating over the PDF. Many variables like this are
509 : ! used in parts of the predictive equation set. These substitutions allow
510 : ! some terms to solved implicitly or semi-implicitly in the predictive
511 : ! equations.
512 : !
513 : !
514 : !=============================================================================
515 0 : subroutine calc_setter_parameters( xm, xp2, Skx, sgn_wpxp, & ! In
516 : big_L_x_1, big_L_x_2, & ! In
517 : mu_x_1, mu_x_2, sigma_x_1_sqd, & ! Out
518 : sigma_x_2_sqd, mixt_frac, & ! Out
519 : coef_sigma_x_1_sqd, & ! Out
520 : coef_sigma_x_2_sqd ) ! Out
521 :
522 : ! Description:
523 : ! Calculates the PDF component means, the PDF component standard deviations,
524 : ! and the mixture fraction for the variable that sets the PDF.
525 :
526 : ! References:
527 : !-----------------------------------------------------------------------
528 :
529 : use constants_clubb, only: &
530 : four, & ! Variable(s)
531 : three, &
532 : one, &
533 : one_half, &
534 : zero, &
535 : eps
536 :
537 : use clubb_precision, only: &
538 : core_rknd ! Variable(s)
539 :
540 : implicit none
541 :
542 : ! Input Variables
543 : real( kind = core_rknd ), intent(in) :: &
544 : xm, & ! Mean of x (overall) [units vary]
545 : xp2, & ! Variance of x (overall) [(units vary)^2]
546 : Skx, & ! Skewness of x [-]
547 : sgn_wpxp, & ! Sign of the covariance of w and x (overall) [-]
548 : big_L_x_1, & ! Parameter for the spread of the 1st PDF comp. mean of x [-]
549 : big_L_x_2 ! Parameter for the spread of the 2nd PDF comp. mean of x [-]
550 :
551 : ! Output Variables
552 : real( kind = core_rknd ), intent(out) :: &
553 : mu_x_1, & ! Mean of x (1st PDF component) [units vary]
554 : mu_x_2, & ! Mean of x (2nd PDF component) [units vary]
555 : sigma_x_1_sqd, & ! Variance of x (1st PDF component) [(units vary)^2]
556 : sigma_x_2_sqd, & ! Variance of x (2nd PDF component) [(units vary)^2]
557 : mixt_frac ! Mixture fraction b [-]
558 :
559 : real( kind = core_rknd ), intent(out) :: &
560 : coef_sigma_x_1_sqd, & ! sigma_x_1^2 = coef_sigma_x_1_sqd * <x'^2> [-]
561 : coef_sigma_x_2_sqd ! sigma_x_2^2 = coef_sigma_x_2_sqd * <x'^2> [-]
562 :
563 : ! Local Variables
564 : real( kind = core_rknd ) :: &
565 : mu_x_1_nrmlized, & ! Normalized mean of x (1st PDF component) [-]
566 : mu_x_2_nrmlized ! Normalized mean of x (2nd PDF component) [-]
567 :
568 : real( kind = core_rknd ) :: &
569 : factor_plus, &
570 : factor_minus, &
571 : sqrt_factor_plus_ov_minus, &
572 : sqrt_factor_minus_ov_plus, &
573 : mu_x_1_nrmlized_thresh
574 :
575 :
576 : ! Calculate the factors in the PDF component mean equations.
577 0 : factor_plus = one + Skx * sgn_wpxp / sqrt( four + Skx**2 )
578 :
579 0 : factor_minus = one - Skx * sgn_wpxp / sqrt( four + Skx**2 )
580 :
581 0 : sqrt_factor_plus_ov_minus = sqrt( factor_plus / factor_minus )
582 0 : sqrt_factor_minus_ov_plus = sqrt( factor_minus / factor_plus )
583 :
584 : ! Calculate the normalized mean of x in the 1st PDF component.
585 0 : mu_x_1_nrmlized = big_L_x_1 * sqrt_factor_plus_ov_minus * sgn_wpxp
586 :
587 : ! Calculate the normalized mean of x in the 2nd PDF component.
588 0 : mu_x_2_nrmlized = -big_L_x_2 * sqrt_factor_minus_ov_plus * sgn_wpxp
589 :
590 : ! Calculate the mean of x in the 1st PDF component.
591 0 : mu_x_1 = xm + mu_x_1_nrmlized * sqrt( xp2 )
592 :
593 : ! Calculate the mean of x in the 2nd PDF component.
594 0 : mu_x_2 = xm + mu_x_2_nrmlized * sqrt( xp2 )
595 :
596 : ! Calculate the mixture fraction.
597 : if ( abs( mu_x_1_nrmlized ) >= eps &
598 0 : .and. abs( mu_x_2_nrmlized ) >= eps ) then
599 0 : mixt_frac = one / ( one + abs( mu_x_1_nrmlized / mu_x_2_nrmlized ) )
600 : elseif ( abs( mu_x_1_nrmlized ) >= eps &
601 0 : .and. abs( mu_x_2_nrmlized ) < eps ) then
602 0 : mixt_frac = one / ( one + abs( mu_x_1_nrmlized / eps ) )
603 : elseif ( abs( mu_x_1_nrmlized ) < eps &
604 0 : .and. abs( mu_x_2_nrmlized ) >= eps ) then
605 0 : mixt_frac = one / ( one + abs( eps / mu_x_2_nrmlized ) )
606 : else ! abs( mu_x_1_nrmlized ) < eps and abs( mu_x_2_nrmlized ) < eps
607 0 : mixt_frac = one_half
608 : endif
609 :
610 : ! Use a minimum magnitude value of mu_x_1_nrmlized in the denominator of a
611 : ! term in the PDF component variance equations in order to prevent a
612 : ! divide-by-zero error.
613 0 : if ( mu_x_1_nrmlized >= zero ) then
614 0 : mu_x_1_nrmlized_thresh = max( mu_x_1_nrmlized, eps )
615 : else ! mu_x_1_nrmlized < 0
616 0 : mu_x_1_nrmlized_thresh = min( mu_x_1_nrmlized, -eps )
617 : endif ! mu_x_1_nrmlized >= 0
618 :
619 : ! Calculate the variance of x in the 1st PDF component.
620 : coef_sigma_x_1_sqd &
621 : = one - mixt_frac * mu_x_1_nrmlized**2 &
622 : - ( one - mixt_frac ) * mu_x_2_nrmlized**2 &
623 : + ( one - mixt_frac ) &
624 : * ( Skx / ( three * mixt_frac * mu_x_1_nrmlized_thresh ) &
625 0 : - mu_x_1_nrmlized**2 / three + mu_x_2_nrmlized**2 / three )
626 :
627 0 : sigma_x_1_sqd = coef_sigma_x_1_sqd * xp2
628 :
629 : ! Calculate the variance of x in the 2nd PDF component.
630 : coef_sigma_x_2_sqd &
631 : = one - mixt_frac * mu_x_1_nrmlized**2 &
632 : - ( one - mixt_frac ) * mu_x_2_nrmlized**2 &
633 : - mixt_frac &
634 : * ( Skx / ( three * mixt_frac * mu_x_1_nrmlized_thresh ) &
635 0 : - mu_x_1_nrmlized**2 / three + mu_x_2_nrmlized**2 / three )
636 :
637 0 : sigma_x_2_sqd = coef_sigma_x_2_sqd * xp2
638 :
639 :
640 0 : return
641 :
642 : end subroutine calc_setter_parameters
643 :
644 : !=============================================================================
645 0 : subroutine calc_L_x_Skx_fnc( Skx, sgn_wpxp, & ! In
646 : small_l_x_1, small_l_x_2, & ! In
647 : big_L_x_1, big_L_x_2 ) ! Out
648 :
649 : ! Description:
650 : ! Calculates the values of big_L_x_1 and big_L_x_2 as functions of Skx.
651 :
652 : ! References:
653 : !-----------------------------------------------------------------------
654 :
655 : use constants_clubb, only: &
656 : four, & ! Variable(s)
657 : zero
658 :
659 : use clubb_precision, only: &
660 : core_rknd ! Variable(s)
661 :
662 : implicit none
663 :
664 : ! Input Variables
665 : real( kind = core_rknd ), intent(in) :: &
666 : Skx, & ! Skewness of x (overall) [-]
667 : sgn_wpxp, & ! Sign of the covariance of w and x (overall) [-]
668 : small_l_x_1, & ! Param. for the spread of the 1st PDF comp. mean of x [-]
669 : small_l_x_2 ! Param. for the spread of the 2nd PDF comp. mean of x [-]
670 :
671 : ! Output Variable
672 : real( kind = core_rknd ), intent(out) :: &
673 : big_L_x_1, & ! Parameter for the spread of the 1st PDF comp. mean of x [-]
674 : big_L_x_2 ! Parameter for the spread of the 2nd PDF comp. mean of x [-]
675 :
676 : ! Local Variable
677 : real( kind = core_rknd ) :: &
678 : Skx_fnc_factor
679 :
680 :
681 : ! The values of L_x_1 and L_x_2 are calculated by skewness functions.
682 : ! Those functions are:
683 : !
684 : ! L_x_1 = l_x_1 * abs( Skx ) / sqrt( 4 + Skx^2 ); and
685 : ! L_x_2 = l_x_2 * abs( Skx ) / sqrt( 4 + Skx^2 ).
686 : !
687 : ! The conditions on l_x_1 and l_x_2 are:
688 : !
689 : ! 2/3 < l_x_1 < 1 and 0 < l_x_2 < 1; when Skx * sgn( <w'x'> ) >= 0; and
690 : ! 0 < l_x_1 < 1 and 2/3 < l_x_2 < 1; when Skx * sgn( <w'x'> ) < 0.
691 : !
692 : ! For simplicity, this can also be accomplished by setting 2/3 < l_x_1 < 1
693 : ! and 0 < l_x_2 < 1, and then using the following equations.
694 : !
695 : ! When Skx * sgn( <w'x'> ) >= 0:
696 : ! L_x_1 = l_x_1 * abs( Skx ) / sqrt( 4 + Skx^2 ); and
697 : ! L_x_2 = l_x_2 * abs( Skx ) / sqrt( 4 + Skx^2 );
698 : !
699 : ! otherwise, when Skx * sgn( <w'x'> ) < 0, switch l_x_1 and l_x_2:
700 : ! L_x_1 = l_x_2 * abs( Skx ) / sqrt( 4 + Skx^2 ); and
701 : ! L_x_2 = l_x_1 * abs( Skx ) / sqrt( 4 + Skx^2 ).
702 0 : Skx_fnc_factor = abs( Skx ) / sqrt( four + Skx**2 )
703 :
704 0 : if ( Skx * sgn_wpxp >= zero ) then
705 0 : big_L_x_1 = small_l_x_1 * Skx_fnc_factor
706 0 : big_L_x_2 = small_l_x_2 * Skx_fnc_factor
707 : else ! Skx * sgn_wpxp < 0
708 0 : big_L_x_1 = small_l_x_2 * Skx_fnc_factor
709 0 : big_L_x_2 = small_l_x_1 * Skx_fnc_factor
710 : endif
711 :
712 :
713 0 : return
714 :
715 : end subroutine calc_L_x_Skx_fnc
716 :
717 : !=============================================================================
718 : !
719 : ! DESCRIPTION OF THE METHOD FOR EACH RESPONDING VARIABLE
720 : ! ======================================================
721 : !
722 : ! In order to find equations for the four PDF parameters for each responding
723 : ! variable, which are mu_x_1, mu_x_2, sigma_x_1, and sigma_x_2 (where x stands
724 : ! for a responding variable here), four equations are needed. These four
725 : ! equations are the equations for <x>, <x'^2>, and <x'^3> as found by
726 : ! integrating over the PDF. Additionally, one more equation, which involves
727 : ! tunable parameter L_x_1, and which is used to help control the mean of the
728 : ! 1st PDF component, is used in this equation set. The four equations are:
729 : !
730 : ! <x> = mixt_frac * mu_x_1 + ( 1 - mixt_frac ) * mu_x_2;
731 : !
732 : ! <x'^2> = mixt_frac * ( ( mu_x_1 - <x> )^2 + sigma_x_1^2 )
733 : ! + ( 1 - mixt_frac ) * ( ( mu_x_2 - <x> )^2 + sigma_x_2^2 );
734 : !
735 : ! <x'^3> = mixt_frac * ( mu_x_1 - <x> )
736 : ! * ( ( mu_x_1 - <x> )^2 + 3 * sigma_x_1^2 )
737 : ! + ( 1 - mixt_frac ) * ( mu_x_2 - <x> )
738 : ! * ( ( mu_x_2 - <x> )^2 + 3 * sigma_x_2^2 ); and
739 : !
740 : ! mu_x_1 - <x> = L_x_1
741 : ! * sqrt( ( 1 + Skx * sgn( <w'x'> ) / sqrt( 4 + Skx^2 ) )
742 : ! / ( 1 - Skx * sgn( <w'x'> ) / sqrt( 4 + Skx^2 ) ) )
743 : ! * sqrt( <x'^2> ) * sgn( <w'x'> );
744 : !
745 : ! where 0 <= L_x_1 <= 1, Skx is the skewness of x, such that
746 : ! Skx = <x'^3> / <x'^2>^(3/2), and sgn( <w'x'> ) is given by:
747 : !
748 : ! sgn( <w'x'> ) = | 1, when <w'x'> >= 0;
749 : ! | -1, when <w'x'> < 0.
750 : !
751 : ! The resulting equations for the four PDF parameters are:
752 : !
753 : ! mu_x_1 = <x> + L_x_1
754 : ! * sqrt( ( 1 + Skx * sgn( <w'x'> ) / sqrt( 4 + Skx^2 ) )
755 : ! / ( 1 - Skx * sgn( <w'x'> ) / sqrt( 4 + Skx^2 ) ) )
756 : ! * sqrt( <x'^2> ) * sgn( <w'x'> );
757 : !
758 : ! mu_x_2 = <x> - ( mixt_frac / ( 1 - mixt_frac ) ) * ( mu_x_1 - <x> );
759 : !
760 : ! sigma_x_1 = sqrt( ( 1 - mixt_frac * mu_x_1_nrmlized^2
761 : ! - ( 1 - mixt_frac ) * mu_x_2_nrmlized^2
762 : ! + ( 1 - mixt_frac )
763 : ! * ( Skx / ( 3 * mixt_frac * mu_x_1_nrmlized )
764 : ! - mu_x_1_nrmlized^2 / 3
765 : ! + mu_x_2_nrmlized^2 / 3 ) )
766 : ! * <x'^2> ); and
767 : !
768 : ! sigma_x_2 = sqrt( ( 1 - mixt_frac * mu_x_1_nrmlized^2
769 : ! - ( 1 - mixt_frac ) * mu_x_2_nrmlized^2
770 : ! - mixt_frac
771 : ! * ( Skx / ( 3 * mixt_frac * mu_x_1_nrmlized )
772 : ! - mu_x_1_nrmlized^2 / 3
773 : ! + mu_x_2_nrmlized^2 / 3 ) )
774 : ! * <x'^2> ); where
775 : !
776 : ! mu_x_1_nrmlized = ( mu_x_1 - <x> ) / sqrt( <x'^2> ); and
777 : !
778 : ! mu_x_2_nrmlized = ( mu_x_2 - <x> ) / sqrt( <x'^2> ).
779 : !
780 : !
781 : ! Notes:
782 : !
783 : ! When L_x_1 = 0, mu_x_1 = mu_x_2 = <x>, sigma_x_1^2 = sigma_x_2^2 = <x'^2>,
784 : ! and the distribution reduces to a single Gaussian.
785 : !
786 : !
787 : ! Equations for PDF component standard deviations:
788 : !
789 : ! The equations for the PDF component standard deviations can also be written
790 : ! as:
791 : !
792 : ! sigma_x_1 = sqrt( coef_sigma_x_1_sqd * <x'^2> ); and
793 : !
794 : ! sigma_x_2 = sqrt( coef_sigma_x_2_sqd * <x'^2> ); where
795 : !
796 : ! coef_sigma_x_1_sqd = 1 - mixt_frac * mu_x_1_nrmlized^2
797 : ! - ( 1 - mixt_frac ) * mu_x_2_nrmlized^2
798 : ! + ( 1 - mixt_frac )
799 : ! * ( Skx / ( 3 * mixt_frac * mu_x_1_nrmlized )
800 : ! - mu_x_1_nrmlized^2 / 3
801 : ! + mu_x_2_nrmlized^2 / 3 ); and
802 : !
803 : ! coef_sigma_x_2_sqd = 1 - mixt_frac * mu_x_1_nrmlized^2
804 : ! - ( 1 - mixt_frac ) * mu_x_2_nrmlized^2
805 : ! - mixt_frac
806 : ! * ( Skx / ( 3 * mixt_frac * mu_x_1_nrmlized )
807 : ! - mu_x_1_nrmlized^2 / 3
808 : ! + mu_x_2_nrmlized^2 / 3 ).
809 : !
810 : ! The above equations can be substituted into an equation for a variable that
811 : ! has been derived by integrating over the PDF. Many variables like this are
812 : ! used in parts of the predictive equation set. These substitutions allow
813 : ! some terms to solved implicitly or semi-implicitly in the predictive
814 : ! equations.
815 : !
816 : !
817 : !=============================================================================
818 0 : subroutine calc_respnder_parameters( xm, xp2, Skx, sgn_wpxp, & ! In
819 : mixt_frac, big_L_x_1, & ! In
820 : mu_x_1, mu_x_2, & ! Out
821 : sigma_x_1_sqd, sigma_x_2_sqd, & ! Out
822 : coef_sigma_x_1_sqd, & ! Out
823 : coef_sigma_x_2_sqd ) ! Out
824 :
825 : ! Description:
826 : ! Calculates the PDF component means, the PDF component standard deviations,
827 : ! and the mixture fraction for the variable that sets the PDF.
828 :
829 : ! References:
830 : !-----------------------------------------------------------------------
831 :
832 : use constants_clubb, only: &
833 : four, & ! Variable(s)
834 : three, &
835 : one, &
836 : zero, &
837 : eps
838 :
839 : use clubb_precision, only: &
840 : core_rknd ! Variable(s)
841 :
842 : implicit none
843 :
844 : ! Input Variables
845 : real( kind = core_rknd ), intent(in) :: &
846 : xm, & ! Mean of x (overall) [units vary]
847 : xp2, & ! Variance of x (overall) [(units vary)^2]
848 : Skx, & ! Skewness of x [-]
849 : sgn_wpxp, & ! Sign of the covariance of w and x (overall) [-]
850 : mixt_frac, & ! Mixture fraction [-]
851 : big_L_x_1 ! Parameter for the spread of the 1st PDF comp. mean of x [-]
852 :
853 : ! Output Variables
854 : real( kind = core_rknd ), intent(out) :: &
855 : mu_x_1, & ! Mean of x (1st PDF component) [units vary]
856 : mu_x_2, & ! Mean of x (2nd PDF component) [units vary]
857 : sigma_x_1_sqd, & ! Variance of x (1st PDF component) [(units vary)^2]
858 : sigma_x_2_sqd ! Variance of x (2nd PDF component) [(units vary)^2]
859 :
860 : real( kind = core_rknd ), intent(out) :: &
861 : coef_sigma_x_1_sqd, & ! sigma_x_1^2 = coef_sigma_x_1_sqd * <x'^2> [-]
862 : coef_sigma_x_2_sqd ! sigma_x_2^2 = coef_sigma_x_2_sqd * <x'^2> [-]
863 :
864 : ! Local Variables
865 : real( kind = core_rknd ) :: &
866 : mu_x_1_nrmlized, & ! Normalized mean of x (1st PDF component) [-]
867 : mu_x_2_nrmlized ! Normalized mean of x (2nd PDF component) [-]
868 :
869 : real( kind = core_rknd ) :: &
870 : factor_plus, &
871 : factor_minus, &
872 : sqrt_factor_plus_ov_minus, &
873 : mu_x_1_nrmlized_thresh
874 :
875 :
876 : ! Calculate the factors in the PDF component mean equations.
877 0 : factor_plus = one + Skx * sgn_wpxp / sqrt( four + Skx**2 )
878 :
879 0 : factor_minus = one - Skx * sgn_wpxp / sqrt( four + Skx**2 )
880 :
881 0 : sqrt_factor_plus_ov_minus = sqrt( factor_plus / factor_minus )
882 :
883 : ! Calculate the normalized mean of x in the 1st PDF component.
884 0 : mu_x_1_nrmlized = big_L_x_1 * sqrt_factor_plus_ov_minus * sgn_wpxp
885 :
886 : ! Calculate the normalized mean of x in the 2nd PDF component.
887 0 : mu_x_2_nrmlized = - ( mixt_frac / ( one - mixt_frac ) ) * mu_x_1_nrmlized
888 :
889 : ! Calculate the mean of x in the 1st PDF component.
890 0 : mu_x_1 = xm + mu_x_1_nrmlized * sqrt( xp2 )
891 :
892 : ! Calculate the mean of x in the 2nd PDF component.
893 0 : mu_x_2 = xm + mu_x_2_nrmlized * sqrt( xp2 )
894 :
895 : ! Use a minimum magnitude value of mu_x_1_nrmlized in the denominator of a
896 : ! term in the PDF component variance equations in order to prevent a
897 : ! divide-by-zero error.
898 0 : if ( mu_x_1_nrmlized >= zero ) then
899 0 : mu_x_1_nrmlized_thresh = max( mu_x_1_nrmlized, eps )
900 : else ! mu_x_1_nrmlized < 0
901 0 : mu_x_1_nrmlized_thresh = min( mu_x_1_nrmlized, -eps )
902 : endif ! mu_x_1_nrmlized >= 0
903 :
904 : ! Calculate the variance of x in the 1st PDF component.
905 : coef_sigma_x_1_sqd &
906 : = one - mixt_frac * mu_x_1_nrmlized**2 &
907 : - ( one - mixt_frac ) * mu_x_2_nrmlized**2 &
908 : + ( one - mixt_frac ) &
909 : * ( Skx / ( three * mixt_frac * mu_x_1_nrmlized_thresh ) &
910 0 : - mu_x_1_nrmlized**2 / three + mu_x_2_nrmlized**2 / three )
911 :
912 0 : sigma_x_1_sqd = coef_sigma_x_1_sqd * xp2
913 :
914 : ! Calculate the variance of x in the 2nd PDF component.
915 : coef_sigma_x_2_sqd &
916 : = one - mixt_frac * mu_x_1_nrmlized**2 &
917 : - ( one - mixt_frac ) * mu_x_2_nrmlized**2 &
918 : - mixt_frac &
919 : * ( Skx / ( three * mixt_frac * mu_x_1_nrmlized_thresh ) &
920 0 : - mu_x_1_nrmlized**2 / three + mu_x_2_nrmlized**2 / three )
921 :
922 0 : sigma_x_2_sqd = coef_sigma_x_2_sqd * xp2
923 :
924 :
925 0 : return
926 :
927 : end subroutine calc_respnder_parameters
928 :
929 : !=============================================================================
930 :
931 : end module new_tsdadg_pdf
|