Line data Source code
1 : !-------------------------------------------------------------------------
2 : ! $Id$
3 : !===============================================================================
4 : module pdf_utilities
5 :
6 : implicit none
7 :
8 : private ! Set default scope to private
9 :
10 : public :: mean_L2N, &
11 : mean_L2N_dp, &
12 : stdev_L2N, &
13 : stdev_L2N_dp, &
14 : corr_NL2NN, &
15 : corr_NL2NN_dp, &
16 : corr_NN2NL, &
17 : corr_LL2NN, &
18 : corr_LL2NN_dp, &
19 : corr_NN2LL, &
20 : compute_mean_binormal, &
21 : compute_variance_binormal, &
22 : calc_comp_corrs_binormal, &
23 : calc_corr_chi_x, &
24 : calc_corr_eta_x, &
25 : calc_corr_rt_x, &
26 : calc_corr_thl_x, &
27 : calc_xp2
28 :
29 : contains
30 :
31 : !=============================================================================
32 0 : elemental function mean_L2N( mu_x, sigma2_on_mu2 ) &
33 : result( mu_x_n )
34 :
35 : ! Description:
36 : ! For a lognormally-distributed variable x, this function finds the mean of
37 : ! ln x (mu_x_n) for the ith component of the PDF, given the mean of x (mu_x)
38 : ! and the variance of x (sigma_sqd_x) for the ith component of the PDF. The
39 : ! value ln x is distributed normally when x is distributed lognormally.
40 :
41 : ! References:
42 : ! Garvey, P. R., 2000: Probability methods for cost uncertainty analysis.
43 : ! Marcel Dekker, 401 pp.
44 : ! -- App. B.
45 : !-----------------------------------------------------------------------
46 :
47 : use constants_clubb, only: &
48 : one ! Constant(s)
49 :
50 : use clubb_precision, only: &
51 : core_rknd ! Variable(s)
52 :
53 : implicit none
54 :
55 : ! Input Variables
56 : real( kind = core_rknd ), intent(in) :: &
57 : mu_x, & ! Mean of x (ith PDF component) [-]
58 : sigma2_on_mu2 ! Ratio: sigma_x^2 / mu_x^2 (ith PDF component) [-]
59 :
60 : ! Return Variable
61 : real( kind = core_rknd ) :: &
62 : mu_x_n ! Mean of ln x (ith PDF component) [-]
63 :
64 :
65 : ! Find the mean of ln x for the ith component of the PDF.
66 : ! The max( mu_x / sqrt( 1 + sigma_x^2 / mu_x^2 ), tiny( mu_x ) ) statement
67 : ! is used to prevent taking ln 0, which will produce a result of -infinity.
68 : ! This would happen when mu_x is 0. However, this code should not be
69 : ! entered when mu_x has a value of 0.
70 0 : mu_x_n = log( max( mu_x / sqrt( one + sigma2_on_mu2 ), tiny( mu_x ) ) )
71 :
72 :
73 : return
74 :
75 : end function mean_L2N
76 :
77 : !=============================================================================
78 0 : elemental function mean_L2N_dp( mu_x, sigma2_on_mu2 ) &
79 : result( mu_x_n )
80 :
81 : ! Description:
82 : ! For a lognormally-distributed variable x, this function finds the mean of
83 : ! ln x (mu_x_n) for the ith component of the PDF, given the mean of x (mu_x)
84 : ! and the variance of x (sigma_sqd_x) for the ith component of the PDF. The
85 : ! value ln x is distributed normally when x is distributed lognormally.
86 : ! This function uses double precision variables.
87 :
88 : ! References:
89 : ! Garvey, P. R., 2000: Probability methods for cost uncertainty analysis.
90 : ! Marcel Dekker, 401 pp.
91 : ! -- App. B.
92 : !-----------------------------------------------------------------------
93 :
94 : use constants_clubb, only: &
95 : one_dp ! Constant(s)
96 :
97 : use clubb_precision, only: &
98 : dp ! double precision
99 :
100 : implicit none
101 :
102 : ! Input Variables
103 : real( kind = dp ), intent(in) :: &
104 : mu_x, & ! Mean of x (ith PDF component) [-]
105 : sigma2_on_mu2 ! Ratio: sigma_x^2 / mu_x^2 (ith PDF component) [-]
106 :
107 : ! Return Variable
108 : real( kind = dp ) :: &
109 : mu_x_n ! Mean of ln x (ith PDF component) [-]
110 :
111 :
112 : ! Find the mean of ln x for the ith component of the PDF.
113 0 : mu_x_n = log( mu_x / sqrt( one_dp + sigma2_on_mu2 ) )
114 :
115 :
116 : return
117 :
118 : end function mean_L2N_dp
119 :
120 : !=============================================================================
121 0 : elemental function stdev_L2N( sigma2_on_mu2 ) &
122 : result( sigma_x_n )
123 :
124 : ! Description:
125 : ! For a lognormally-distributed variable x, this function finds the standard
126 : ! deviation of ln x (sigma_x_n) for the ith component of the PDF, given the
127 : ! mean of x (mu_x) and the variance of x (sigma_sqd_x) for the ith component
128 : ! of the PDF. The value ln x is distributed normally when x is distributed
129 : ! lognormally.
130 :
131 : ! References:
132 : ! Garvey, P. R., 2000: Probability methods for cost uncertainty analysis.
133 : ! Marcel Dekker, 401 pp.
134 : ! -- App. B.
135 : !-----------------------------------------------------------------------
136 :
137 : use constants_clubb, only: &
138 : one ! Constant(s)
139 :
140 : use clubb_precision, only: &
141 : core_rknd ! Variable(s)
142 :
143 : implicit none
144 :
145 : ! Input Variables
146 : real( kind = core_rknd ), intent(in) :: &
147 : sigma2_on_mu2 ! Ratio: sigma_x^2 / mu_x^2 (ith PDF component) [-]
148 :
149 : ! Return Variable
150 : real( kind = core_rknd ) :: &
151 : sigma_x_n ! Standard deviation of ln x (ith PDF component) [-]
152 :
153 :
154 : ! Find the standard deviation of ln x for the ith component of the PDF.
155 0 : sigma_x_n = sqrt( log( one + sigma2_on_mu2 ) )
156 :
157 :
158 : return
159 :
160 : end function stdev_L2N
161 :
162 : !=============================================================================
163 0 : elemental function stdev_L2N_dp( sigma2_on_mu2 ) &
164 : result( sigma_x_n )
165 :
166 : ! Description:
167 : ! For a lognormally-distributed variable x, this function finds the standard
168 : ! deviation of ln x (sigma_x_n) for the ith component of the PDF, given the
169 : ! mean of x (mu_x) and the variance of x (sigma_sqd_x) for the ith component
170 : ! of the PDF. The value ln x is distributed normally when x is distributed
171 : ! lognormally.
172 : ! This function uses double precision variables.
173 :
174 : ! References:
175 : ! Garvey, P. R., 2000: Probability methods for cost uncertainty analysis.
176 : ! Marcel Dekker, 401 pp.
177 : ! -- App. B.
178 : !-----------------------------------------------------------------------
179 :
180 : use constants_clubb, only: &
181 : one_dp ! Constant(s)
182 :
183 : use clubb_precision, only: &
184 : dp ! double precision
185 :
186 : implicit none
187 :
188 : ! Input Variables
189 : real( kind = dp ), intent(in) :: &
190 : sigma2_on_mu2 ! Ratio: sigma_x^2 / mu_x^2 (ith PDF component) [-]
191 :
192 : ! Return Variable
193 : real( kind = dp ) :: &
194 : sigma_x_n ! Standard deviation of ln x (ith PDF component) [-]
195 :
196 :
197 : ! Find the standard deviation of ln x for the ith component of the PDF.
198 0 : sigma_x_n = sqrt( log( one_dp + sigma2_on_mu2 ) )
199 :
200 :
201 : return
202 :
203 : end function stdev_L2N_dp
204 :
205 : !=============================================================================
206 0 : elemental function corr_NL2NN( corr_x_y, sigma_y_n, y_sigma2_on_mu2 ) &
207 : result( corr_x_y_n )
208 :
209 : ! Description:
210 : ! For a normally-distributed variable x and a lognormally-distributed
211 : ! variable y, this function finds the correlation of x and ln y (corr_x_y_n)
212 : ! for the ith component of the PDF, given the correlation of x and y
213 : ! (corr_x_y) and the standard deviation of ln y (sigma_y_n) for the ith
214 : ! component of the PDF. The value ln y is distributed normally when y is
215 : ! distributed lognormally.
216 :
217 : ! References:
218 : ! Garvey, P. R., 2000: Probability methods for cost uncertainty analysis.
219 : ! Marcel Dekker, 401 pp.
220 : ! -- Eq. B-1.
221 : !-----------------------------------------------------------------------
222 :
223 : use constants_clubb, only: &
224 : max_mag_correlation, & ! Constant(s)
225 : zero
226 :
227 : use clubb_precision, only: &
228 : core_rknd ! Variable(s)
229 :
230 : implicit none
231 :
232 : ! Input Variables
233 : real( kind = core_rknd ), intent(in) :: &
234 : corr_x_y, & ! Correlation of x and y (ith PDF component) [-]
235 : sigma_y_n, & ! Standard deviation of ln y (ith PDF component) [-]
236 : y_sigma2_on_mu2 ! Ratio: sigma_y^2 / mu_y^2 (ith PDF component) [-]
237 :
238 : ! Return Variable
239 : real( kind = core_rknd ) :: &
240 : corr_x_y_n ! Correlation of x and ln y (ith PDF component) [-]
241 :
242 :
243 : ! Find the correlation of x and ln y for the ith component of the PDF.
244 : ! When sigma_y = 0 and mu_y > 0, y_sigma2_on_mu2 = 0. This results in
245 : ! sigma_y_n = 0. The resulting corr_x_y_n is undefined. However, the
246 : ! divide-by-zero problem needs to be addressed in the code.
247 0 : if ( sigma_y_n > zero ) then
248 0 : corr_x_y_n = corr_x_y * sqrt( y_sigma2_on_mu2 ) / sigma_y_n
249 : else ! sigma_y_n = 0
250 : ! The value of sqrt( y_sigma2_on_mu2 ) / sigma_y_n can be rewritten as:
251 : ! sqrt( y_sigma2_on_mu2 ) / sqrt( ln( 1 + y_sigma2_on_mu2 ) ).
252 : ! This can be further rewritten as:
253 : ! sqrt( y_sigma2_on_mu2 / ln( 1 + y_sigma2_on_mu2 ) ),
254 : ! which has a limit of 1 as y_sigma2_on_mu2 approaches 0 from the right.
255 : ! When sigma_y_n = 0, the value of corr_x_y_n is undefined, so set it
256 : ! to corr_x_y.
257 0 : corr_x_y_n = corr_x_y
258 : endif ! sigma_y_n > 0
259 :
260 : ! Clip the magnitude of the correlation of x and ln y in the ith PDF
261 : ! component, just in case the correlation (ith PDF component) of x and y and
262 : ! the standard deviation (ith PDF component) of ln y are inconsistent,
263 : ! resulting in an unrealizable value for corr_x_y_n.
264 0 : if ( corr_x_y_n > max_mag_correlation ) then
265 : corr_x_y_n = max_mag_correlation
266 0 : elseif ( corr_x_y_n < -max_mag_correlation ) then
267 0 : corr_x_y_n = -max_mag_correlation
268 : endif
269 :
270 :
271 : return
272 :
273 : end function corr_NL2NN
274 :
275 : !=============================================================================
276 0 : elemental function corr_NL2NN_dp( corr_x_y, sigma_y_n, y_sigma2_on_mu2 ) &
277 : result( corr_x_y_n )
278 :
279 : ! Description:
280 : ! For a normally-distributed variable x and a lognormally-distributed
281 : ! variable y, this function finds the correlation of x and ln y (corr_x_y_n)
282 : ! for the ith component of the PDF, given the correlation of x and y
283 : ! (corr_x_y) and the standard deviation of ln y (sigma_y_n) for the ith
284 : ! component of the PDF. The value ln y is distributed normally when y is
285 : ! distributed lognormally.
286 : ! This function uses double precision variables.
287 :
288 : ! References:
289 : ! Garvey, P. R., 2000: Probability methods for cost uncertainty analysis.
290 : ! Marcel Dekker, 401 pp.
291 : ! -- Eq. B-1.
292 : !-----------------------------------------------------------------------
293 :
294 : use constants_clubb, only: &
295 : max_mag_correlation, & ! Constant(s)
296 : zero_dp
297 :
298 : use clubb_precision, only: &
299 : dp ! double precision
300 :
301 : implicit none
302 :
303 : ! Input Variables
304 : real( kind = dp ), intent(in) :: &
305 : corr_x_y, & ! Correlation of x and y (ith PDF component) [-]
306 : sigma_y_n, & ! Standard deviation of ln y (ith PDF component) [-]
307 : y_sigma2_on_mu2 ! Ratio: sigma_y^2 / mu_y^2 (ith PDF component) [-]
308 :
309 : ! Return Variable
310 : real( kind = dp ) :: &
311 : corr_x_y_n ! Correlation of x and ln y (ith PDF component) [-]
312 :
313 :
314 : ! Find the correlation of x and ln y for the ith component of the PDF.
315 : ! When sigma_y = 0 and mu_y > 0, y_sigma2_on_mu2 = 0. This results in
316 : ! sigma_y_n = 0. The resulting corr_x_y_n is undefined. However, the
317 : ! divide-by-zero problem needs to be addressed in the code.
318 0 : if ( sigma_y_n > zero_dp ) then
319 0 : corr_x_y_n = corr_x_y * sqrt( y_sigma2_on_mu2 ) / sigma_y_n
320 : else ! sigma_y_n = 0
321 : ! The value of sqrt( y_sigma2_on_mu2 ) / sigma_y_n can be rewritten as:
322 : ! sqrt( y_sigma2_on_mu2 ) / sqrt( ln( 1 + y_sigma2_on_mu2 ) ).
323 : ! This can be further rewritten as:
324 : ! sqrt( y_sigma2_on_mu2 / ln( 1 + y_sigma2_on_mu2 ) ),
325 : ! which has a limit of 1 as y_sigma2_on_mu2 approaches 0 from the right.
326 : ! When sigma_y_n = 0, the value of corr_x_y_n is undefined, so set it
327 : ! to corr_x_y.
328 0 : corr_x_y_n = corr_x_y
329 : endif ! sigma_y_n > 0
330 :
331 : ! Clip the magnitude of the correlation of x and ln y in the ith PDF
332 : ! component, just in case the correlation (ith PDF component) of x and y and
333 : ! the standard deviation (ith PDF component) of ln y are inconsistent,
334 : ! resulting in an unrealizable value for corr_x_y_n.
335 0 : if ( corr_x_y_n > real( max_mag_correlation, kind = dp ) ) then
336 : corr_x_y_n = real( max_mag_correlation, kind = dp )
337 0 : elseif ( corr_x_y_n < -real( max_mag_correlation, kind = dp ) ) then
338 0 : corr_x_y_n = -real( max_mag_correlation, kind = dp )
339 : endif
340 :
341 :
342 : return
343 :
344 : end function corr_NL2NN_dp
345 :
346 : !=============================================================================
347 0 : subroutine corr_NN2NL( nz, ngrdcol, &
348 0 : corr_x_y_n, sigma_y_n, y_sigma2_on_mu2, &
349 0 : corr_x_y )
350 :
351 : ! Description:
352 : ! For a normally-distributed variable x and a lognormally-distributed
353 : ! variable y, this function finds the correlation of x and y (corr_x_y) for
354 : ! the ith component of the PDF, given the correlation of x and ln y
355 : ! (corr_x_y_n) and the standard deviation of ln y (sigma_y_n) for the ith
356 : ! component of the PDF. The value ln y is distributed normally when y is
357 : ! distributed lognormally.
358 :
359 : ! References:
360 : ! Garvey, P. R., 2000: Probability methods for cost uncertainty analysis.
361 : ! Marcel Dekker, 401 pp.
362 : ! -- Eq. B-1.
363 : !-----------------------------------------------------------------------
364 :
365 : use constants_clubb, only: &
366 : max_mag_correlation, & ! Constant(s)
367 : zero
368 :
369 : use clubb_precision, only: &
370 : core_rknd ! Variable(s)
371 :
372 : implicit none
373 :
374 : ! Input Variables
375 : integer, intent(in) :: &
376 : nz, & ! Number of vertical levels
377 : ngrdcol ! Number of grid columns
378 :
379 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
380 : corr_x_y_n, & ! Correlation of x and ln y (ith PDF component) [-]
381 : sigma_y_n, & ! Standard deviation of ln y (ith PDF component) [-]
382 : y_sigma2_on_mu2 ! Ratio: sigma_y^2 / mu_y^2 (ith PDF component) [-]
383 :
384 : ! Output Variable
385 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
386 : corr_x_y ! Correlation of x and y (ith PDF component) [-]
387 :
388 : ! Local variables
389 : integer :: &
390 : j, k ! Loop iterators
391 :
392 :
393 : ! Find the correlation of x and y for the ith component of the PDF.
394 : ! When sigma_y = 0 and mu_y > 0, y_sigma2_on_mu2 = 0. This results in
395 : ! sigma_y_n = 0. The resulting corr_x_y and corr_x_y_n are undefined.
396 : ! However, the divide-by-zero problem needs to be addressed in the code.
397 0 : do k = 1, nz
398 0 : do j = 1, ngrdcol
399 :
400 0 : if ( sigma_y_n(j,k) > zero ) then
401 : ! Use the maximum of y_sigma2_on_mu2 and tiny( y_sigma2_on_mu2 ) instead
402 : ! of just y_sigma2_on_mu2. The value of y_sigma2_on_mu2 must already be
403 : ! greater than 0 in order for this block of code to be entered (when
404 : ! y_sigma2_on_mu2 = 0, sigma_y_n = 0, and this block of code is not
405 : ! entered).
406 : corr_x_y(j,k) = corr_x_y_n(j,k) * sigma_y_n(j,k) &
407 0 : / sqrt( max( y_sigma2_on_mu2(j,k), tiny( y_sigma2_on_mu2(j,k) ) ) )
408 : else ! sigma_y_n = 0
409 : ! The value of sigma_y_n / sqrt( y_sigma2_on_mu2 ) can be rewritten as:
410 : ! sqrt( ln( 1 + y_sigma2_on_mu2 ) ) / sqrt( y_sigma2_on_mu2 ).
411 : ! This can be further rewritten as:
412 : ! sqrt( ln( 1 + y_sigma2_on_mu2 ) / y_sigma2_on_mu2 ),
413 : ! which has a limit of 1 as y_sigma2_on_mu2 approaches 0 from the right.
414 : ! When sigma_y_n = 0, the value of corr_x_y is undefined, so set it
415 : ! to corr_x_y_n.
416 0 : corr_x_y(j,k) = corr_x_y_n(j,k)
417 : endif ! sigma_y_n > 0
418 :
419 : ! Clip the magnitude of the correlation of x and y in the ith PDF component,
420 : ! just in case the correlation (ith PDF component) of x and ln y and the
421 : ! standard deviation (ith PDF component) of ln y are inconsistent, resulting
422 : ! in an unrealizable value for corr_x_y.
423 0 : if ( corr_x_y(j,k) > max_mag_correlation ) then
424 0 : corr_x_y(j,k) = max_mag_correlation
425 0 : elseif ( corr_x_y(j,k) < -max_mag_correlation ) then
426 0 : corr_x_y(j,k) = -max_mag_correlation
427 : endif
428 :
429 : end do
430 : end do
431 :
432 0 : return
433 :
434 : end subroutine corr_NN2NL
435 :
436 : !=============================================================================
437 0 : elemental function corr_LL2NN( corr_x_y, sigma_x_n, sigma_y_n, &
438 : x_sigma2_on_mu2, y_sigma2_on_mu2 ) &
439 : result( corr_x_y_n )
440 :
441 : ! Description:
442 : ! For lognormally-distributed variables x and y, this function finds the
443 : ! correlation of ln x and ln y (corr_x_y_n) for the ith component of the
444 : ! PDF, given the correlation of x and y (corr_x_y), the standard deviation
445 : ! of ln x (sigma_x_n), and the standard deviation of ln y (sigma_y_n) for
446 : ! the ith component of the PDF. The value of ln x (or ln y) is distributed
447 : ! normally when x (or y) is distributed lognormally.
448 :
449 : ! References:
450 : ! Garvey, P. R., 2000: Probability methods for cost uncertainty analysis.
451 : ! Marcel Dekker, 401 pp.
452 : ! -- Eq. C-3.
453 : !-----------------------------------------------------------------------
454 :
455 : use constants_clubb, only: &
456 : one, & ! Constant(s)
457 : zero, &
458 : max_mag_correlation
459 :
460 : use clubb_precision, only: &
461 : core_rknd ! Variable(s)
462 :
463 : implicit none
464 :
465 : ! Input Variables
466 : real( kind = core_rknd ), intent(in) :: &
467 : corr_x_y, & ! Correlation of x and y (ith PDF component) [-]
468 : sigma_x_n, & ! Standard deviation of ln x (ith PDF component) [-]
469 : sigma_y_n, & ! Standard deviation of ln y (ith PDF component) [-]
470 : x_sigma2_on_mu2, & ! Ratio: sigma_x^2 / mu_x^2 (ith PDF component) [-]
471 : y_sigma2_on_mu2 ! Ratio: sigma_y^2 / mu_y^2 (ith PDF component) [-]
472 :
473 : ! Return Variable
474 : real( kind = core_rknd ) :: &
475 : corr_x_y_n ! Correlation of ln x and ln y (ith PDF component) [-]
476 :
477 : ! Local Variable
478 : real( kind = core_rknd ) :: &
479 : log_arg ! Input into the ln function [-]
480 :
481 :
482 : ! Find the correlation of ln x and ln y for the ith component of the PDF.
483 : ! When sigma_x = 0 and mu_x > 0, x_sigma2_on_mu2 = 0. This results in
484 : ! sigma_x_n = 0. The resulting corr_x_y_n is undefined. The same holds
485 : ! true when sigma_y = 0 and mu_y > 0. However, the divide-by-zero problem
486 : ! needs to be addressed in the code.
487 0 : if ( sigma_x_n > zero .and. sigma_y_n > zero ) then
488 : ! corr_x_y_n = log( one + corr_x_y * sqrt( exp( sigma_x_n**2 ) - one ) &
489 : ! * sqrt( exp( sigma_y_n**2 ) - one ) ) &
490 : ! / ( sigma_x_n * sigma_y_n )
491 0 : log_arg = one + corr_x_y * sqrt( x_sigma2_on_mu2 * y_sigma2_on_mu2 )
492 0 : corr_x_y_n = log( log_arg ) / ( sigma_x_n * sigma_y_n )
493 : else ! sigma_x_n = 0 or sigma_y_n = 0
494 : ! The value of corr_x_y_n is undefined, so set it to corr_x_y.
495 0 : corr_x_y_n = corr_x_y
496 : endif ! sigma_x_n > 0 and sigma_y_n > 0
497 :
498 : ! Clip the magnitude of the correlation of ln x and ln y in the ith PDF
499 : ! component, just in case the correlation (ith PDF component) of x and y,
500 : ! the standard deviation (ith PDF component) of ln x, and the standard
501 : ! deviation (ith PDF component) of ln y are inconsistent, resulting in an
502 : ! unrealizable value for corr_x_y_n.
503 0 : if ( corr_x_y_n > max_mag_correlation ) then
504 : corr_x_y_n = max_mag_correlation
505 0 : elseif ( corr_x_y_n < -max_mag_correlation ) then
506 0 : corr_x_y_n = -max_mag_correlation
507 : endif
508 :
509 :
510 : return
511 :
512 : end function corr_LL2NN
513 :
514 : !=============================================================================
515 0 : elemental function corr_LL2NN_dp( corr_x_y, sigma_x_n, sigma_y_n, &
516 : x_sigma2_on_mu2, y_sigma2_on_mu2 ) &
517 : result( corr_x_y_n )
518 :
519 : ! Description:
520 : ! For lognormally-distributed variables x and y, this function finds the
521 : ! correlation of ln x and ln y (corr_x_y_n) for the ith component of the
522 : ! PDF, given the correlation of x and y (corr_x_y), the standard deviation
523 : ! of ln x (sigma_x_n), and the standard deviation of ln y (sigma_y_n) for
524 : ! the ith component of the PDF. The value of ln x (or ln y) is distributed
525 : ! normally when x (or y) is distributed lognormally.
526 : ! This function uses double precision variables.
527 :
528 : ! References:
529 : ! Garvey, P. R., 2000: Probability methods for cost uncertainty analysis.
530 : ! Marcel Dekker, 401 pp.
531 : ! -- Eq. C-3.
532 : !-----------------------------------------------------------------------
533 :
534 : use constants_clubb, only: &
535 : one_dp, & ! Constant(s)
536 : zero_dp, &
537 : max_mag_correlation
538 :
539 : use clubb_precision, only: &
540 : dp ! double precision
541 :
542 : implicit none
543 :
544 : ! Input Variables
545 : real( kind = dp ), intent(in) :: &
546 : corr_x_y, & ! Correlation of x and y (ith PDF component) [-]
547 : sigma_x_n, & ! Standard deviation of ln x (ith PDF component) [-]
548 : sigma_y_n, & ! Standard deviation of ln y (ith PDF component) [-]
549 : x_sigma2_on_mu2, & ! Ratio: sigma_x^2 / mu_x^2 (ith PDF component) [-]
550 : y_sigma2_on_mu2 ! Ratio: sigma_y^2 / mu_y^2 (ith PDF component) [-]
551 :
552 :
553 : ! Return Variable
554 : real( kind = dp ) :: &
555 : corr_x_y_n ! Correlation of ln x and ln y (ith PDF component) [-]
556 :
557 :
558 : ! Find the correlation of ln x and ln y for the ith component of the PDF.
559 : ! When sigma_x = 0 and mu_x > 0, x_sigma2_on_mu2 = 0. This results in
560 : ! sigma_x_n = 0. The resulting corr_x_y_n is undefined. The same holds
561 : ! true when sigma_y = 0 and mu_y > 0. However, the divide-by-zero problem
562 : ! needs to be addressed in the code.
563 0 : if ( sigma_x_n > zero_dp .and. sigma_y_n > zero_dp ) then
564 : corr_x_y_n &
565 : = log( one_dp + corr_x_y * sqrt( x_sigma2_on_mu2 * y_sigma2_on_mu2 ) ) &
566 0 : / ( sigma_x_n * sigma_y_n )
567 : else ! sigma_x_n = 0 or sigma_y_n = 0
568 : ! The value of corr_x_y_n is undefined, so set it to corr_x_y.
569 0 : corr_x_y_n = corr_x_y
570 : endif ! sigma_x_n > 0 and sigma_y_n > 0
571 :
572 : ! Clip the magnitude of the correlation of ln x and ln y in the ith PDF
573 : ! component, just in case the correlation (ith PDF component) of x and y,
574 : ! the standard deviation (ith PDF component) of ln x, and the standard
575 : ! deviation (ith PDF component) of ln y are inconsistent, resulting in an
576 : ! unrealizable value for corr_x_y_n.
577 0 : if ( corr_x_y_n > real( max_mag_correlation, kind = dp ) ) then
578 : corr_x_y_n = real( max_mag_correlation, kind = dp )
579 0 : elseif ( corr_x_y_n < -real( max_mag_correlation, kind = dp ) ) then
580 0 : corr_x_y_n = -real( max_mag_correlation, kind = dp )
581 : endif
582 :
583 :
584 : return
585 :
586 : end function corr_LL2NN_dp
587 :
588 : !=============================================================================
589 0 : subroutine corr_NN2LL( nz, ngrdcol, &
590 0 : corr_x_y_n, sigma_x_n, sigma_y_n, &
591 0 : x_sigma2_on_mu2, y_sigma2_on_mu2, &
592 0 : corr_x_y )
593 :
594 : ! Description:
595 : ! For lognormally-distributed variables x and y, this function finds the
596 : ! correlation of x and y (corr_x_y) for the ith component of the PDF, given
597 : ! the correlation of ln x and ln y (corr_x_y_n), the standard deviation of
598 : ! ln x (sigma_x_n), and the standard deviation of ln y (sigma_y_n) for
599 : ! the ith component of the PDF. The value of ln x (or ln y) is distributed
600 : ! normally when x (or y) is distributed lognormally.
601 :
602 : ! References:
603 : ! Garvey, P. R., 2000: Probability methods for cost uncertainty analysis.
604 : ! Marcel Dekker, 401 pp.
605 : ! -- Eq. C-3.
606 : !-----------------------------------------------------------------------
607 :
608 : use constants_clubb, only: &
609 : one, & ! Constant(s)
610 : zero, &
611 : max_mag_correlation
612 :
613 : use clubb_precision, only: &
614 : core_rknd ! Variable(s)
615 :
616 : implicit none
617 :
618 : ! Input Variables
619 : integer, intent(in) :: &
620 : nz, & ! Number of vertical levels
621 : ngrdcol ! Number of grid columns
622 :
623 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
624 : corr_x_y_n, & ! Correlation of ln x and ln y (ith PDF component) [-]
625 : sigma_x_n, & ! Standard deviation of ln x (ith PDF component) [-]
626 : sigma_y_n, & ! Standard deviation of ln y (ith PDF component) [-]
627 : x_sigma2_on_mu2, & ! Ratio: sigma_x^2 / mu_x^2 (ith PDF component) [-]
628 : y_sigma2_on_mu2 ! Ratio: sigma_y^2 / mu_y^2 (ith PDF component) [-]
629 :
630 : ! Output Variable
631 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
632 : corr_x_y ! Correlation of x and y (ith PDF component) [-]
633 :
634 : ! Local variables
635 : integer :: &
636 : j, k ! Loop iterators
637 :
638 : ! Find the correlation of x and y for the ith component of the PDF.
639 : ! When sigma_x = 0 and mu_x > 0, x_sigma2_on_mu2 = 0. This results in
640 : ! sigma_x_n = 0. The resulting corr_x_y and corr_x_y_n are undefined. The
641 : ! same holds true when sigma_y = 0 and mu_y > 0. However, the
642 : ! divide-by-zero problem needs to be addressed in the code.
643 0 : do k = 1, nz
644 0 : do j = 1, ngrdcol
645 :
646 0 : if ( sigma_x_n(j,k) > zero .and. sigma_y_n(j,k) > zero ) then
647 : ! corr_x_y = ( exp( sigma_x_n * sigma_y_n * corr_x_y_n ) - one ) &
648 : ! / ( sqrt( exp( sigma_x_n**2 ) - one ) &
649 : ! * sqrt( exp( sigma_y_n**2 ) - one ) )
650 : corr_x_y(j,k) = ( exp( sigma_x_n(j,k) * sigma_y_n(j,k) * corr_x_y_n(j,k) ) - one ) &
651 0 : / sqrt( x_sigma2_on_mu2(j,k) * y_sigma2_on_mu2(j,k) )
652 : else ! sigma_x_n = 0 or sigma_y_n = 0
653 : ! The value of corr_x_y is undefined, so set it to corr_x_y_n.
654 0 : corr_x_y(j,k) = corr_x_y_n(j,k)
655 : endif ! sigma_x_n > 0 and sigma_y_n > 0
656 :
657 : ! Clip the magnitude of the correlation of x and y in the ith PDF component,
658 : ! just in case the correlation (ith PDF component) of ln x and ln y, the
659 : ! standard deviation (ith PDF component) of ln x, and the standard deviation
660 : ! (ith PDF component) of ln y are inconsistent, resulting in an unrealizable
661 : ! value for corr_x_y.
662 0 : if ( corr_x_y(j,k) > max_mag_correlation ) then
663 0 : corr_x_y(j,k) = max_mag_correlation
664 0 : elseif ( corr_x_y(j,k) < -max_mag_correlation ) then
665 0 : corr_x_y(j,k) = -max_mag_correlation
666 : endif
667 :
668 : end do
669 : end do
670 :
671 0 : return
672 :
673 : end subroutine corr_NN2LL
674 :
675 : !=============================================================================
676 0 : elemental function compute_mean_binormal( mu_x_1, mu_x_2, mixt_frac ) &
677 : result( xm )
678 :
679 : ! Description:
680 : ! Computes the overall grid-box mean of a binormal distribution from the
681 : ! mean of each component
682 :
683 : ! References:
684 : ! None
685 : !-----------------------------------------------------------------------
686 :
687 : use clubb_precision, only: &
688 : core_rknd ! Constant
689 :
690 : use constants_clubb, only: &
691 : one ! Constant
692 :
693 : implicit none
694 :
695 : ! Input Variables
696 : real( kind = core_rknd ), intent(in) :: &
697 : mu_x_1, & ! First PDF component mean of 'x' [?]
698 : mu_x_2, & ! Second PDF component mean of 'x' [?]
699 : mixt_frac ! Weight of the first PDF component [-]
700 :
701 : ! Output Variables
702 : real( kind = core_rknd ) :: &
703 : xm ! Mean of 'x' (overall) [?]
704 :
705 : !-----------------------------------------------------------------------
706 :
707 : !----- Begin Code -----
708 0 : xm = mixt_frac * mu_x_1 + ( one - mixt_frac ) * mu_x_2
709 :
710 :
711 : return
712 :
713 : end function compute_mean_binormal
714 :
715 : !=============================================================================
716 0 : elemental function compute_variance_binormal( xm, mu_x_1, mu_x_2, &
717 : stdev_x_1, stdev_x_2, &
718 : mixt_frac ) &
719 : result( xp2 )
720 :
721 : ! Description:
722 : ! Computes the overall grid-box variance of a binormal distribution from the
723 : ! variance of each component.
724 :
725 : ! References:
726 : ! None
727 : !-----------------------------------------------------------------------
728 :
729 : use clubb_precision, only: &
730 : core_rknd ! Constant
731 :
732 : use constants_clubb, only: &
733 : one ! Constant
734 :
735 : implicit none
736 :
737 : ! Input Variables
738 : real( kind = core_rknd ), intent(in) :: &
739 : xm, & ! Overall mean of 'x' [?]
740 : mu_x_1, & ! First PDF component mean of 'x' [?]
741 : mu_x_2, & ! Second PDF component mean of 'x' [?]
742 : stdev_x_1, & ! Standard deviation of 'x' in the first PDF component [?]
743 : stdev_x_2, & ! Standard deviation of 'x' in the second PDF component [?]
744 : mixt_frac ! Weight of the first PDF component [-]
745 :
746 : ! Output Variables
747 : real( kind = core_rknd ) :: &
748 : xp2 ! Variance of 'x' (overall) [?^2]
749 :
750 : !-----------------------------------------------------------------------
751 :
752 : !----- Begin Code -----
753 : xp2 = mixt_frac * ( ( mu_x_1 - xm )**2 + stdev_x_1**2 ) &
754 0 : + ( one - mixt_frac ) * ( ( mu_x_2 - xm )**2 + stdev_x_2**2 )
755 :
756 :
757 : return
758 :
759 : end function compute_variance_binormal
760 :
761 : !=============================================================================
762 17870112 : subroutine calc_comp_corrs_binormal( nz, ngrdcol, & ! In
763 17870112 : xpyp, xm, ym, & ! In
764 17870112 : mu_x_1, mu_x_2, & ! In
765 17870112 : mu_y_1, mu_y_2, & ! In
766 17870112 : sigma_x_1_sqd, & ! In
767 17870112 : sigma_x_2_sqd, & ! In
768 17870112 : sigma_y_1_sqd, & ! In
769 17870112 : sigma_y_2_sqd, & ! In
770 17870112 : mixt_frac, & ! In
771 17870112 : corr_x_y_1, & ! Out
772 17870112 : corr_x_y_2 ) ! Out
773 :
774 : ! Description:
775 : ! Calculates the PDF component correlations of variables x and y, where
776 : ! x and y are both distributed as two-component normals (or binormals).
777 : ! The PDF component correlations are set equal to each other.
778 : !
779 : ! The overall covariance of x and y, <x'y'>, can be expressed in terms of
780 : ! PDF parameters by integrating over the PDF:
781 : !
782 : ! <x'y'> = INT(-inf:inf) INT(-inf:inf) ( x - <x> ) ( y - <y> ) P(x,y) dy dx;
783 : !
784 : ! where <x> is the overall mean of x, <y> is the overall mean of y, and
785 : ! P(x,y) is the equation for the two-component normal PDF of x and y.
786 : !
787 : ! The integral is evaluated, and the equation for <x'y'> is:
788 : !
789 : ! <x'y'> = mixt_frac * ( ( mu_x_1 - <x> ) * ( mu_y_1 - <y> )
790 : ! + corr_x_y_1 * sigma_x_1 * sigma_y_1 )
791 : ! + ( 1 - mixt_frac ) * ( ( mu_x_2 - <x> ) * ( mu_y_2 - <y> )
792 : ! + corr_x_y_2 * sigma_x_2 * sigma_y_2 );
793 : !
794 : ! where mu_x_1 is the mean of x in the 1st PDF component, mu_x_2 is the mean
795 : ! of x in the 2nd PDF component, mu_y_1 is the mean of y in the 1st PDF
796 : ! component, mu_y_2 is the mean of y in the 2nd PDF component, sigma_x_1 is
797 : ! the standard deviation of x in the 1st PDF component, sigma_x_2 is the
798 : ! standard deviation of x in the 2nd PDF component, sigma_y_1 is the
799 : ! standard deviation of y in the 1st PDF component, sigma_y_2 is the
800 : ! standard deviation of y in the 2nd PDF component, corr_x_y_1 is the
801 : ! correlation of x and y in the 1st PDF component, corr_x_y_2 is the
802 : ! correlation of x and y in the 2nd PDF component, and mixt_frac is the
803 : ! mixture fraction (weight of the 1st PDF component).
804 : !
805 : ! This equation can be rewritten as:
806 : !
807 : ! <x'y'> = mixt_frac * ( mu_x_1 - <x> ) * ( mu_y_1 - <y> )
808 : ! + mixt_frac * corr_x_y_1 * sigma_x_1 * sigma_y_1
809 : ! + ( 1 - mixt_frac ) * ( mu_x_2 - <x> ) * ( mu_y_2 - <y> )
810 : ! + ( 1 - mixt_frac ) * corr_x_y_2 * sigma_x_2 * sigma_y_2.
811 : !
812 : ! Setting the two PDF component correlations equal to each other
813 : ! (corr_x_y_1 = corr_x_y_2), the equation can be solved for the PDF
814 : ! component correlations:
815 : !
816 : ! corr_x_y_1 = corr_x_y_2
817 : ! = ( <x'y'> - mixt_frac * ( mu_x_1 - <x> ) * ( mu_y_1 - <y> )
818 : ! - ( 1 - mixt_frac ) * ( mu_x_2 - <x> ) * ( mu_y_2 - <y> ) )
819 : ! / ( mixt_frac * sigma_x_1 * sigma_y_1
820 : ! + ( 1 - mixt_frac ) * sigma_x_2 * sigma_y_2 );
821 : !
822 : ! where -1 <= corr_x_y_1 = corr_x_y_2 <= 1.
823 : !
824 : ! When sigma_x_1 * sigma_y_1 = 0 and sigma_x_2 * sigma_y_2 = 0, at least one
825 : ! of x or y are constant within each PDF component, and both PDF component
826 : ! correlations are undefined.
827 :
828 : ! References:
829 : !-----------------------------------------------------------------------
830 :
831 : use constants_clubb, only: &
832 : max_mag_correlation, & ! Variable(s)
833 : one, &
834 : zero
835 :
836 : use clubb_precision, only: &
837 : core_rknd ! Variable(s)
838 :
839 : implicit none
840 :
841 : ! ---------------- Input Variables ----------------
842 : integer, intent(in) :: &
843 : nz, & ! Number of vertical levels
844 : ngrdcol ! Number of grid columns
845 :
846 : real ( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
847 : xpyp, & ! Covariance of x and y (overall) [(x units)(y units)]
848 : xm, & ! Mean of x (overall) [x units]
849 : ym, & ! Mean of y (overall) [y units]
850 : mu_x_1, & ! Mean of x (1st PDF component) [x units]
851 : mu_x_2, & ! Mean of x (2nd PDF component) [x units]
852 : mu_y_1, & ! Mean of y (1st PDF component) [y units]
853 : mu_y_2, & ! Mean of y (2nd PDF component) [y units]
854 : sigma_x_1_sqd, & ! Variance of x (1st PDF component) [(x units)^2]
855 : sigma_x_2_sqd, & ! Variance of x (2nd PDF component) [(x units)^2]
856 : sigma_y_1_sqd, & ! Variance of y (1st PDF component) [(y units)^2]
857 : sigma_y_2_sqd, & ! Variance of y (2nd PDF component) [(y units)^2]
858 : mixt_frac ! Mixture fraction [-]
859 :
860 : ! ---------------- Output Variables ----------------
861 : real ( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
862 : corr_x_y_1, & ! Correlation of x and y (1st PDF component) [-]
863 : corr_x_y_2 ! Correlation of x and y (2nd PDF component) [-]
864 :
865 : ! ---------------- Local Variables ----------------
866 : integer :: i, k
867 :
868 : !$acc parallel loop gang vector collapse(2) default(present)
869 1536829632 : do k = 1, nz
870 25380961632 : do i = 1, ngrdcol
871 47688264000 : if ( sigma_x_1_sqd(i,k) * sigma_y_1_sqd(i,k) > zero &
872 71532396000 : .or. sigma_x_2_sqd(i,k) * sigma_y_2_sqd(i,k) > zero ) then
873 :
874 : ! Calculate corr_x_y_1 (which also equals corr_x_y_2).
875 : corr_x_y_1(i,k) &
876 : = ( xpyp(i,k) &
877 : - mixt_frac(i,k) * ( mu_x_1(i,k) - xm(i,k) ) * ( mu_y_1(i,k) - ym(i,k) ) &
878 : - ( one - mixt_frac(i,k) ) * ( mu_x_2(i,k) - xm(i,k) ) * ( mu_y_2(i,k) - ym(i,k))) &
879 : / ( mixt_frac(i,k) * sqrt( sigma_x_1_sqd(i,k) * sigma_y_1_sqd(i,k) ) &
880 23701892776 : + ( one - mixt_frac(i,k) ) * sqrt( sigma_x_2_sqd(i,k) * sigma_y_2_sqd(i,k) ) )
881 :
882 : ! The correlation must fall within the bounds of
883 : ! -max_mag_correlation < corr_x_y_1 (= corr_x_y_2) < max_mag_correlation
884 : corr_x_y_1(i,k) = max( -max_mag_correlation, &
885 23701892776 : min( max_mag_correlation, corr_x_y_1(i,k) ) )
886 :
887 : else ! sigma_x_1^2 * sigma_y_1^2 = 0 and sigma_x_2^2 * sigma_y_2^2 = 0.
888 :
889 : ! The correlation is undefined (output as 0).
890 142239224 : corr_x_y_1(i,k) = zero
891 :
892 : endif
893 : ! Set corr_x_y_2 equal to corr_x_y_1.
894 25363091520 : corr_x_y_2(i,k) = corr_x_y_1(i,k)
895 : end do
896 : end do
897 : !$acc end parallel loop
898 :
899 17870112 : return
900 :
901 : end subroutine calc_comp_corrs_binormal
902 :
903 : !=============================================================================
904 0 : elemental function calc_corr_chi_x( crt_i, cthl_i, &
905 : sigma_rt_i, sigma_thl_i, &
906 : sigma_chi_i, &
907 : corr_rt_x_i, corr_thl_x_i ) &
908 : result( corr_chi_x_i )
909 :
910 : ! Description:
911 : ! This function calculates the correlation of extended liquid water mixing
912 : ! ratio, chi (old s), and a generic variable x, within the ith component of
913 : ! the PDF. The variable chi can be split into mean and turbulent
914 : ! components, such that:
915 : !
916 : ! chi = <chi> + chi';
917 : !
918 : ! where < > denotes a mean field an ' denotes a turbulent component.
919 : !
920 : ! The linearized equation for chi' is given in Larson et al. (2001), where
921 : ! within the ith component of the PDF:
922 : !
923 : ! chi_(i)' = Coef_rt(i) * r_t(i)' - Coef_thl(i) * th_l(i)'.
924 : !
925 : ! The equation for chi' can be multiplied by x'. The equation becomes:
926 : !
927 : ! chi'x'_(i) = Coef_rt(i) * r_t'x'_(i) - Coef_thl(i) * th_l'x'_(i).
928 : !
929 : ! Averaging both sides, the covariance <chi'x'> is given by the equation:
930 : !
931 : ! <chi'x'_(i)> = Coef_rt(i) * <r_t'x'_(i)> - Coef_thl(i) * <th_l'x'_(i)>.
932 : !
933 : ! This equation can be rewritten as:
934 : !
935 : ! sigma_chi(i) * sigma_x(i) * corr_chi_x(i)
936 : ! = Coef_rt(i) * sigma_rt(i) * sigma_x(i) * corr_rt_x(i)
937 : ! - Coef_thl(i) * sigma_thl(i) * sigma_x(i) * corr_thl_x(i).
938 : !
939 : ! This equation can be solved for corr_chi_x(i):
940 : !
941 : ! corr_chi_x(i)
942 : ! = Coef_rt(i) * ( sigma_rt(i) / sigma_chi(i) ) * corr_rt_x(i)
943 : ! - Coef_thl(i) * ( sigma_thl(i) / sigma_chi(i) ) * corr_thl_x(i).
944 : !
945 : ! The correlation of chi and x within the ith component of the PDF is
946 : ! calculated.
947 :
948 : ! References:
949 : ! Eq. (13) and Eq. (14) of Larson, V. E., R. Wood, P. R. Field, J.-C. Golaz,
950 : ! T. H. Vonder Haar, W. R. Cotton, 2001: Systematic Biases in the
951 : ! Microphysics and Thermodynamics of Numerical Models That Ignore
952 : ! Subgrid-Scale Variability. J. Atmos. Sci., 58, 1117--1128,
953 : ! doi:https://doi.org/10.1175/1520-0469(2001)058%3C1117:SBITMA%3E2.0.CO;2.
954 : !
955 : ! Eq. (A29) of Griffin, B. M., 2016: Improving the Subgrid-Scale
956 : ! Representation of Hydrometeors and Microphysical Feedback Effects Using a
957 : ! Multivariate PDF. Doctoral dissertation, University of
958 : ! Wisconsin -- Milwaukee, Milwaukee, WI, Paper 1144, 165 pp., URL
959 : ! http://dc.uwm.edu/cgi/viewcontent.cgi?article=2149&context=etd.
960 : !-----------------------------------------------------------------------
961 :
962 : use constants_clubb, only: &
963 : zero, & ! Constant(s)
964 : chi_tol, &
965 : max_mag_correlation
966 :
967 : use clubb_precision, only: &
968 : core_rknd ! Variable(s)
969 :
970 : implicit none
971 :
972 : ! Input Variables
973 : real( kind = core_rknd ), intent(in) :: &
974 : crt_i, & ! Coefficient of r_t for chi (old s) (ith PDF comp.) [-]
975 : cthl_i, & ! Coefficient of th_l for chi (ith PDF comp.) [(kg/kg)/K]
976 : sigma_rt_i, & ! Standard deviation of r_t (ith PDF component) [kg/kg]
977 : sigma_thl_i, & ! Standard deviation of th_l (ith PDF component) [K]
978 : sigma_chi_i, & ! Standard deviation of chi (ith PDF component) [kg/kg]
979 : corr_rt_x_i, & ! Correlation of r_t and x (ith PDF component) [-]
980 : corr_thl_x_i ! Correlation of th_l and x (ith PDF component) [-]
981 :
982 : ! Return Variable
983 : real( kind = core_rknd ) :: &
984 : corr_chi_x_i ! Correlation of chi and x (ith PDF component) [-]
985 :
986 :
987 : ! Calculate the correlation of chi and x in the ith PDF component.
988 0 : if ( sigma_chi_i > chi_tol ) then
989 :
990 : corr_chi_x_i = crt_i * ( sigma_rt_i / sigma_chi_i ) * corr_rt_x_i &
991 0 : - cthl_i * ( sigma_thl_i / sigma_chi_i ) * corr_thl_x_i
992 :
993 : else ! sigma_chi_i = 0
994 :
995 : ! The standard deviation of chi in the ith PDF component is 0. This
996 : ! means that chi is constant within the ith PDF component, and the ith
997 : ! PDF component covariance of chi and x is also 0. The correlation of
998 : ! chi and x is undefined in the ith PDF component, so a value of 0 will
999 : ! be used.
1000 : corr_chi_x_i = zero
1001 :
1002 : endif
1003 :
1004 : ! Clip the magnitude of the correlation of chi and x in the ith PDF
1005 : ! component, just in case the correlations and standard deviations used in
1006 : ! calculating it are inconsistent, resulting in an unrealizable value for
1007 : ! corr_chi_x_i.
1008 0 : if ( corr_chi_x_i > max_mag_correlation ) then
1009 : corr_chi_x_i = max_mag_correlation
1010 0 : elseif ( corr_chi_x_i < -max_mag_correlation ) then
1011 0 : corr_chi_x_i = -max_mag_correlation
1012 : endif
1013 :
1014 :
1015 : return
1016 :
1017 : end function calc_corr_chi_x
1018 :
1019 : !=============================================================================
1020 0 : elemental function calc_corr_eta_x( crt_i, cthl_i, &
1021 : sigma_rt_i, sigma_thl_i, &
1022 : sigma_eta_i, corr_rt_x_i, &
1023 : corr_thl_x_i ) &
1024 : result( corr_eta_x_i )
1025 :
1026 : ! Description:
1027 : ! This function calculates the correlation of the variable that is
1028 : ! orthogonal to extended liquid water mixing ratio in a PDF transformation,
1029 : ! eta (old t), and a generic variable x, within the ith component of
1030 : ! the PDF.
1031 :
1032 : ! References:
1033 : ! Eq. (A30) of Griffin, B. M., 2016: Improving the Subgrid-Scale
1034 : ! Representation of Hydrometeors and Microphysical Feedback Effects Using a
1035 : ! Multivariate PDF. Doctoral dissertation, University of
1036 : ! Wisconsin -- Milwaukee, Milwaukee, WI, Paper 1144, 165 pp., URL
1037 : ! http://dc.uwm.edu/cgi/viewcontent.cgi?article=2149&context=etd.
1038 : !-----------------------------------------------------------------------
1039 :
1040 : use constants_clubb, only: &
1041 : zero, & ! Constant(s)
1042 : eta_tol, &
1043 : max_mag_correlation
1044 :
1045 : use clubb_precision, only: &
1046 : core_rknd ! Variable(s)
1047 :
1048 : implicit none
1049 :
1050 : ! Input Variables
1051 : real( kind = core_rknd ), intent(in) :: &
1052 : crt_i, & ! Coefficient of r_t for chi (old s) (ith PDF comp.) [-]
1053 : cthl_i, & ! Coefficient of th_l for chi (ith PDF comp.) [(kg/kg)/K]
1054 : sigma_rt_i, & ! Standard deviation of r_t (ith PDF component) [kg/kg]
1055 : sigma_thl_i, & ! Standard deviation of th_l (ith PDF component) [K]
1056 : sigma_eta_i, & ! Standard deviation of eta (ith PDF component) [kg/kg]
1057 : corr_rt_x_i, & ! Correlation of r_t and x (ith PDF component) [-]
1058 : corr_thl_x_i ! Correlation of th_l and x (ith PDF component) [-]
1059 :
1060 : ! Return Variable
1061 : real( kind = core_rknd ) &
1062 : corr_eta_x_i ! Correlation of eta and x (ith PDF component) [-]
1063 :
1064 :
1065 : ! Calculate the correlation of eta and x in the ith PDF component.
1066 0 : if ( sigma_eta_i > eta_tol ) then
1067 :
1068 : corr_eta_x_i = crt_i * ( sigma_rt_i / sigma_eta_i ) * corr_rt_x_i &
1069 0 : + cthl_i * ( sigma_thl_i / sigma_eta_i ) * corr_thl_x_i
1070 :
1071 : else ! sigma_eta_i = 0
1072 :
1073 : ! The standard deviation of eta in the ith PDF component is 0. This
1074 : ! means that eta is constant within the ith PDF component, and the ith
1075 : ! PDF component covariance of eta and x is also 0. The correlation of
1076 : ! eta and x is undefined in the ith PDF component, so a value of 0 will
1077 : ! be used.
1078 : corr_eta_x_i = zero
1079 :
1080 : endif
1081 :
1082 : ! Clip the magnitude of the correlation of eta and x in the ith PDF
1083 : ! component, just in case the correlations and standard deviations used in
1084 : ! calculating it are inconsistent, resulting in an unrealizable value for
1085 : ! corr_eta_x_i.
1086 0 : if ( corr_eta_x_i > max_mag_correlation ) then
1087 : corr_eta_x_i = max_mag_correlation
1088 0 : elseif ( corr_eta_x_i < -max_mag_correlation ) then
1089 0 : corr_eta_x_i = -max_mag_correlation
1090 : endif
1091 :
1092 :
1093 : return
1094 :
1095 : end function calc_corr_eta_x
1096 :
1097 : !=============================================================================
1098 0 : elemental function calc_corr_rt_x( crt_i, sigma_rt_i, sigma_chi_i, &
1099 : sigma_eta_i, corr_chi_x_i, corr_eta_x_i ) &
1100 : result( corr_rt_x_i )
1101 :
1102 : ! Description:
1103 : ! This function calculates the correlation of rt and x based on the
1104 : ! correlation of chi and x and the correlation of eta and x.
1105 :
1106 : ! References:
1107 : !-----------------------------------------------------------------------
1108 :
1109 : use constants_clubb, only: &
1110 : two, & ! Constant(s)
1111 : zero, &
1112 : rt_tol, &
1113 : max_mag_correlation
1114 :
1115 : use clubb_precision, only: &
1116 : core_rknd ! Variable(s)
1117 :
1118 : implicit none
1119 :
1120 : ! Input Variables
1121 : real( kind = core_rknd ), intent(in) :: &
1122 : crt_i, & ! Coef. of r_t in chi/eta eqns. (ith PDF component) [-]
1123 : sigma_rt_i, & ! Standard deviation of r_t (ith PDF component) [kg/kg]
1124 : sigma_chi_i, & ! Standard deviation of chi (ith PDF component) [kg/kg]
1125 : sigma_eta_i, & ! Standard deviation of eta (ith PDF component) [kg/kg]
1126 : corr_chi_x_i, & ! Correlation of chi and x (ith PDF component) [-]
1127 : corr_eta_x_i ! Correlation of eta and x (ith PDF component) [-]
1128 :
1129 : ! Return Variable
1130 : real( kind = core_rknd ) :: &
1131 : corr_rt_x_i ! Correlation of rt and x (ith PDF component) [-]
1132 :
1133 :
1134 : ! Calculate the correlation of rt and x in the ith PDF component.
1135 0 : if ( sigma_rt_i > rt_tol ) then
1136 :
1137 : corr_rt_x_i = ( sigma_eta_i * corr_eta_x_i &
1138 : + sigma_chi_i * corr_chi_x_i ) &
1139 0 : / ( two * crt_i * sigma_rt_i )
1140 :
1141 : else ! sigma_rt_i = 0
1142 :
1143 : ! The standard deviation of rt in the ith PDF component is 0. This means
1144 : ! that rt is constant within the ith PDF component, and the ith PDF
1145 : ! component covariance of rt and x is also 0. The correlation of rt and
1146 : ! x is undefined in the ith PDF component, so a value of 0 will be used.
1147 : corr_rt_x_i = zero
1148 :
1149 : endif
1150 :
1151 : ! Clip the magnitude of the correlation of rt and x in the ith PDF
1152 : ! component, just in case the correlations and standard deviations used in
1153 : ! calculating it are inconsistent, resulting in an unrealizable value for
1154 : ! corr_rt_x_i.
1155 0 : if ( corr_rt_x_i > max_mag_correlation ) then
1156 : corr_rt_x_i = max_mag_correlation
1157 0 : elseif ( corr_rt_x_i < -max_mag_correlation ) then
1158 0 : corr_rt_x_i = -max_mag_correlation
1159 : endif
1160 :
1161 :
1162 : return
1163 :
1164 : end function calc_corr_rt_x
1165 :
1166 : !=============================================================================
1167 0 : elemental function calc_corr_thl_x( cthl_i, sigma_thl_i, sigma_chi_i, &
1168 : sigma_eta_i, corr_chi_x_i, &
1169 : corr_eta_x_i ) &
1170 : result( corr_thl_x_i )
1171 :
1172 : ! Description:
1173 : ! This function calculates the correlation of thl and x based on the
1174 : ! correlation of chi and x and the correlation of eta and x.
1175 :
1176 : ! References:
1177 : !-----------------------------------------------------------------------
1178 :
1179 : use constants_clubb, only: &
1180 : two, & ! Constant(s)
1181 : zero, &
1182 : thl_tol, &
1183 : max_mag_correlation
1184 :
1185 : use clubb_precision, only: &
1186 : core_rknd ! Variable(s)
1187 :
1188 : implicit none
1189 :
1190 : ! Input Variables
1191 : real( kind = core_rknd ), intent(in) :: &
1192 : cthl_i, & ! Coef. of thl: chi/eta eqns. (ith PDF comp.) [(kg/kg)/K]
1193 : sigma_thl_i, & ! Standard deviation of thl (ith PDF component) [K]
1194 : sigma_chi_i, & ! Standard deviation of chi (ith PDF component) [kg/kg]
1195 : sigma_eta_i, & ! Standard deviation of eta (ith PDF component) [kg/kg]
1196 : corr_chi_x_i, & ! Correlation of chi and x (ith PDF component) [-]
1197 : corr_eta_x_i ! Correlation of eta and x (ith PDF component) [-]
1198 :
1199 : ! Return Variable
1200 : real( kind = core_rknd ) :: &
1201 : corr_thl_x_i ! Correlation of thl and x (ith PDF component) [-]
1202 :
1203 :
1204 : ! Calculate the correlation of thl and x in the ith PDF component.
1205 0 : if ( sigma_thl_i > thl_tol ) then
1206 :
1207 : corr_thl_x_i = ( sigma_eta_i * corr_eta_x_i &
1208 : - sigma_chi_i * corr_chi_x_i ) &
1209 0 : / ( two * cthl_i * sigma_thl_i )
1210 :
1211 : else ! sigma_thl_i = 0
1212 :
1213 : ! The standard deviation of thl in the ith PDF component is 0. This
1214 : ! means that thl is constant within the ith PDF component, and the ith
1215 : ! PDF component covariance of thl and x is also 0. The correlation of
1216 : ! thl and x is undefined in the ith PDF component, so a value of 0 will
1217 : ! be used.
1218 : corr_thl_x_i = zero
1219 :
1220 : endif
1221 :
1222 : ! Clip the magnitude of the correlation of thl and x in the ith PDF
1223 : ! component, just in case the correlations and standard deviations used in
1224 : ! calculating it are inconsistent, resulting in an unrealizable value for
1225 : ! corr_thl_x_i.
1226 0 : if ( corr_thl_x_i > max_mag_correlation ) then
1227 : corr_thl_x_i = max_mag_correlation
1228 0 : elseif ( corr_thl_x_i < -max_mag_correlation ) then
1229 0 : corr_thl_x_i = -max_mag_correlation
1230 : endif
1231 :
1232 :
1233 : return
1234 :
1235 : end function calc_corr_thl_x
1236 :
1237 : !=============================================================================
1238 0 : elemental function calc_xp2( mu_x_1, mu_x_2, &
1239 : mu_x_1_n, mu_x_2_n, &
1240 : sigma_x_1, sigma_x_2, &
1241 : sigma_x_1_n, sigma_x_2_n, &
1242 : mixt_frac, x_frac_1, x_frac_2, &
1243 : x_mean ) &
1244 : result( xp2 )
1245 :
1246 : ! Description:
1247 : ! Calculates the overall variance of x, <x'^2>, where the distribution of x
1248 : ! is a combination of a lognormal distribution and/or 0 in each PDF
1249 : ! component. The fraction of each component where x is lognormally
1250 : ! distributed (amd greater than 0) is x_frac_i (x_frac_1 and x_frac_2 for
1251 : ! PDF components 1 and 2, respectively). The fraction of each component
1252 : ! where x has a value of 0 is ( 1 - x_frac_i ). This function should be
1253 : ! called to calculate the total variance for x when <x'^2> is not provided
1254 : ! by a predictive (or other) equation.
1255 : !
1256 : ! This function is used to calculate the overall variance for rain water
1257 : ! mixing ratio, <r_r'^2>, and the overall variance for rain drop
1258 : ! concentration, <N_r'^2>. The ratio of variance to mean-value-squared is
1259 : ! specified for the in-precip values of r_r and N_r within each PDF
1260 : ! component, allowing for the calculation of sigma_rr_i and sigma_Nr_i,
1261 : ! as well as sigma_rr_i_n and sigma_Nr_i_n.
1262 :
1263 : ! References:
1264 : !-----------------------------------------------------------------------
1265 :
1266 : use constants_clubb, only: &
1267 : two, & ! Constant(s)
1268 : one, &
1269 : zero, &
1270 : eps
1271 :
1272 : use clubb_precision, only: &
1273 : core_rknd ! Variable(s)
1274 :
1275 : implicit none
1276 :
1277 : ! Input Variables
1278 : real( kind = core_rknd ), intent(in) :: &
1279 : mu_x_1, & ! Mean of x (1st PDF comp.) in x_frac [-]
1280 : mu_x_2, & ! Mean of x (2nd PDF comp.) in x_frac [-]
1281 : mu_x_1_n, & ! Mean of ln x (1st PDF comp.) in x_frac [-]
1282 : mu_x_2_n, & ! Mean of ln x (2nd PDF comp.) in x_frac [-]
1283 : sigma_x_1, & ! Standard deviation of x (1st PDF comp.) in x_frac [-]
1284 : sigma_x_2, & ! Standard deviation of x (2nd PDF comp.) in x_frac [-]
1285 : sigma_x_1_n, & ! Standard deviation of ln x (1st PDF comp.) in x_frac [-]
1286 : sigma_x_2_n, & ! Standard deviation of ln x (2nd PDF comp.) in x_frac [-]
1287 : mixt_frac, & ! Mixture fraction [-]
1288 : x_frac_1, & ! Fraction: x distributed lognormally (1st PDF comp.) [-]
1289 : x_frac_2, & ! Fraction: x distributed lognormally (2nd PDF comp.) [-]
1290 : x_mean ! Overall mean value of x [-]
1291 :
1292 : ! Return Variable
1293 : real( kind = core_rknd ) :: &
1294 : xp2 ! Overall variance of x [-]
1295 :
1296 :
1297 : ! Calculate overall variance of x, <x'^2>.
1298 0 : if ( abs(sigma_x_1) < eps .and. abs(sigma_x_2) < eps ) then
1299 :
1300 : ! The value of x is constant within both PDF components.
1301 : xp2 = ( mixt_frac * x_frac_1 * mu_x_1**2 &
1302 : + ( one - mixt_frac ) * x_frac_2 * mu_x_2**2 &
1303 : ) &
1304 0 : - x_mean**2
1305 :
1306 :
1307 0 : elseif ( abs(sigma_x_1) < eps ) then
1308 :
1309 : ! The value of x is constant within the 1st PDF component.
1310 : xp2 = ( mixt_frac * x_frac_1 * mu_x_1**2 &
1311 : + ( one - mixt_frac ) * x_frac_2 &
1312 : * exp( two * mu_x_2_n + two * sigma_x_2_n**2 ) &
1313 : ) &
1314 0 : - x_mean**2
1315 :
1316 :
1317 0 : elseif ( abs(sigma_x_2) < eps ) then
1318 :
1319 : ! The value of x is constant within the 2nd PDF component.
1320 : xp2 = ( mixt_frac * x_frac_1 &
1321 : * exp( two * mu_x_1_n + two * sigma_x_1_n**2 ) &
1322 : + ( one - mixt_frac ) * x_frac_2 * mu_x_2**2 &
1323 : ) &
1324 0 : - x_mean**2
1325 :
1326 :
1327 : else ! sigma_x_1 and sigma_x_2 > 0
1328 :
1329 : ! The value of x varies within both PDF component.
1330 : xp2 = ( mixt_frac * x_frac_1 &
1331 : * exp( two * mu_x_1_n + two * sigma_x_1_n**2 ) &
1332 : + ( one - mixt_frac ) * x_frac_2 &
1333 : * exp( two * mu_x_2_n + two * sigma_x_2_n**2 ) &
1334 : ) &
1335 0 : - x_mean**2
1336 :
1337 :
1338 : endif
1339 :
1340 :
1341 : ! As a check, prevent negative values for hydrometeor variances due to
1342 : ! numerical loss of precision error.
1343 0 : if ( xp2 < zero ) then
1344 0 : xp2 = zero
1345 : endif
1346 :
1347 :
1348 : return
1349 :
1350 : end function calc_xp2
1351 :
1352 : !===============================================================================
1353 :
1354 : end module pdf_utilities
|