Line data Source code
1 : ! $Id$
2 : !===============================================================================
3 : module LY93_pdf
4 :
5 : ! Description:
6 : ! The multivariate, two-component PDF of Lewellen and Yoh (1993).
7 :
8 : ! References:
9 : ! Lewellen, W. S. and Yoh, S., 1993. Binormal Model of Ensemble Partial
10 : ! Cloudiness. J. Atmos. Sci., 50, 9, 1228--1237.
11 : !-------------------------------------------------------------------------
12 :
13 : implicit none
14 :
15 : public :: LY93_driver, & ! Procedure(s)
16 : calc_mixt_frac_LY93, &
17 : calc_params_LY93
18 :
19 : private ! default scope
20 :
21 : contains
22 :
23 : !=============================================================================
24 0 : subroutine LY93_driver( nz, wm, rtm, thlm, wp2, rtp2, & ! In
25 0 : thlp2, Skw, Skrt, Skthl, & ! In
26 0 : mu_w_1, mu_w_2, & ! Out
27 0 : mu_rt_1, mu_rt_2, & ! Out
28 0 : mu_thl_1, mu_thl_2, & ! Out
29 0 : sigma_w_1_sqd, sigma_w_2_sqd, & ! Out
30 0 : sigma_rt_1_sqd, sigma_rt_2_sqd, & ! Out
31 0 : sigma_thl_1_sqd, sigma_thl_2_sqd, & ! Out
32 0 : mixt_frac ) ! Out
33 :
34 : ! Description:
35 : ! Calculates the mixture fraction and the PDF component means and PDF
36 : ! component variances of w, rt, and theta-l following Lewellen and Yoh.
37 :
38 : ! References:
39 : ! Lewellen, W. S. and Yoh, S., 1993. Binormal Model of Ensemble Partial
40 : ! Cloudiness. J. Atmos. Sci., 50, 9, 1228--1237.
41 : !-----------------------------------------------------------------------
42 :
43 : use clubb_precision, only: &
44 : core_rknd ! Variable(s)
45 :
46 : implicit none
47 :
48 : integer, intent(in) :: &
49 : nz
50 :
51 : ! Input Variables
52 : real( kind = core_rknd), dimension(nz), intent(in) :: &
53 : wm, & ! Mean of w (overall) [m/s]
54 : wp2, & ! Variance of w (overall) [m^2/s^2]
55 : Skw, & ! Skewness of w (overall) [-]
56 : rtm, & ! Mean of rt (overall) [kg/kg]
57 : rtp2, & ! Variance of rt (overall) [kg^2/kg^2]
58 : Skrt, & ! Skewness of rt (overall) [-]
59 : thlm, & ! Mean of thl (overall) [K]
60 : thlp2, & ! Variance of thl (overall) [K^2]
61 : Skthl ! Skewness of thl (overall) [-]
62 :
63 : ! Output Variables
64 : real( kind = core_rknd), dimension(nz), intent(out) :: &
65 : mu_w_1, & ! Mean of w (1st PDF component) [m/s]
66 : mu_w_2, & ! Mean of w (2nd PDF component) [m/s]
67 : mu_rt_1, & ! Mean of rt (1st PDF component) [kg/kg]
68 : mu_rt_2, & ! Mean of rt (2nd PDF component) [kg/kg]
69 : mu_thl_1, & ! Mean of thl (1st PDF component) [K]
70 : mu_thl_2, & ! Mean of thl (2nd PDF component) [K]
71 : sigma_w_1_sqd, & ! Variance of w (1st PDF component) [m^2/s^2]
72 : sigma_w_2_sqd, & ! Variance of w (2nd PDF component) [m^2/s^2]
73 : sigma_rt_1_sqd, & ! Variance of rt (1st PDF component) [m^2/s^2]
74 : sigma_rt_2_sqd, & ! Variance of rt (2nd PDF component) [m^2/s^2]
75 : sigma_thl_1_sqd, & ! Variance of thl (1st PDF component) [m^2/s^2]
76 : sigma_thl_2_sqd, & ! Variance of thl (2nd PDF component) [m^2/s^2]
77 : mixt_frac ! Mixture fraction [-]
78 :
79 : ! Local Variables
80 : real( kind = core_rknd), dimension(nz) :: &
81 0 : Sk_max ! Maximum of magnitudes of skewness [-]
82 :
83 :
84 : ! Find the maximum of the magnitudes of skewness.
85 0 : Sk_max = max( abs( Skw ), abs( Skrt ), abs( Skthl ) )
86 :
87 : ! Calculate mixture fraction.
88 0 : mixt_frac = calc_mixt_frac_LY93( nz, Sk_max )
89 :
90 : ! Calculate the PDF parameters for w.
91 : call calc_params_LY93( nz, wm, wp2, Skw, mixt_frac, & ! In
92 : mu_w_1, mu_w_2, & ! Out
93 0 : sigma_w_1_sqd, sigma_w_2_sqd ) ! Out
94 :
95 : ! Calculate the PDF parameters for rt.
96 : call calc_params_LY93( nz, rtm, rtp2, Skrt, mixt_frac, & ! In
97 : mu_rt_1, mu_rt_2, & ! Out
98 0 : sigma_rt_1_sqd, sigma_rt_2_sqd ) ! Out
99 :
100 : ! Calculate the PDF parameters for thl.
101 : call calc_params_LY93( nz, thlm, thlp2, Skthl, mixt_frac, & ! In
102 : mu_thl_1, mu_thl_2, & ! Out
103 0 : sigma_thl_1_sqd, sigma_thl_2_sqd ) ! Out
104 :
105 :
106 0 : return
107 :
108 : end subroutine LY93_driver
109 :
110 : !=============================================================================
111 0 : function calc_mixt_frac_LY93( nz, Sk_max ) &
112 0 : result( mixt_frac )
113 :
114 : ! Description:
115 : ! Calculates mixture fraction iteratively according to Lewellen and Yoh.
116 :
117 : ! References:
118 : ! Eq. (21) of Lewellen, W. S. and Yoh, S., 1993. Binormal Model of Ensemble
119 : ! Partial Cloudiness. J. Atmos. Sci., 50, 9, 1228--1237.
120 : !-----------------------------------------------------------------------
121 :
122 : use constants_clubb, only: &
123 : one, & ! Constant(s)
124 : three_fourths, &
125 : one_half, &
126 : zero
127 :
128 : use clubb_precision, only: &
129 : core_rknd ! Variable(s)
130 :
131 : implicit none
132 :
133 : integer, intent(in) :: &
134 : nz
135 :
136 : ! Input Variable
137 : real( kind = core_rknd), dimension(nz), intent(in) :: &
138 : Sk_max ! Maximum of magnitudes of skewness [-]
139 :
140 : ! Return Variable
141 : real( kind = core_rknd), dimension(nz) :: &
142 : mixt_frac ! Mixture fraction [-]
143 :
144 : ! Local Variables
145 : real( kind = core_rknd) :: &
146 : mixt_frac_low, & ! Low value of mixture frac. in iterative solver [-]
147 : mixt_frac_high, & ! High value of mixture frac.in iterative solver [-]
148 : expr_equal_zero ! Expr. mixt_frac^6 - Sk_max * ( 1 - mixt_frac ) [-]
149 :
150 : ! Tolerance for mixture fraction in solver [-]
151 : real( kind = core_rknd) :: &
152 : LY_mixt_frac_tol = 1.0e-4_core_rknd
153 :
154 : integer :: k ! Vertical level index
155 :
156 :
157 0 : do k = 1, nz, 1
158 :
159 0 : if ( Sk_max(k) > 0.84_core_rknd ) then
160 :
161 : mixt_frac_low = one_half
162 : mixt_frac_high = one
163 :
164 : do ! solve iteratively for mixture fraction
165 :
166 0 : mixt_frac(k) = one_half * ( mixt_frac_low + mixt_frac_high )
167 :
168 : expr_equal_zero &
169 0 : = mixt_frac(k)**6 - Sk_max(k)**2 * ( one - mixt_frac(k) )
170 :
171 0 : if ( abs( expr_equal_zero ) < LY_mixt_frac_tol ) then
172 : ! Mixture fraction has been solved for within the specificed
173 : ! tolerance.
174 : exit
175 : else
176 0 : if ( expr_equal_zero > zero ) then
177 : mixt_frac_high = mixt_frac(k)
178 : else ! expr_equal_zero < 0
179 0 : mixt_frac_low = mixt_frac(k)
180 : endif
181 : endif
182 :
183 : enddo ! solve iteratively for mixture fraction
184 :
185 : else ! Sk_max <= 0.84
186 :
187 0 : mixt_frac(k) = three_fourths
188 :
189 : endif
190 :
191 : enddo ! k = 1, nz, 1
192 :
193 :
194 0 : return
195 :
196 0 : end function calc_mixt_frac_LY93
197 :
198 : !=============================================================================
199 0 : subroutine calc_params_LY93( nz, xm, xp2, Skx, mixt_frac, & ! In
200 0 : mu_x_1, mu_x_2, & ! Out
201 0 : sigma_x_1_sqd, sigma_x_2_sqd ) ! Out
202 :
203 : ! Description:
204 : ! Calculates the PDF component means and PDF component variances for
205 : ! variable x according to Lewellen and Yoh.
206 :
207 : ! References:
208 : ! Eq. (14), Eq. (15), Eq. (16), Eq. (17), and Eq. (18) of
209 : ! Lewellen, W. S. and Yoh, S., 1993. Binormal Model of Ensemble Partial
210 : ! Cloudiness. J. Atmos. Sci., 50, 9, 1228--1237.
211 : !-----------------------------------------------------------------------
212 :
213 : use constants_clubb, only: &
214 : three, & ! Constant(s)
215 : one, &
216 : one_third, &
217 : zero
218 :
219 : use clubb_precision, only: &
220 : core_rknd ! Variable(s)
221 :
222 : implicit none
223 :
224 : integer, intent(in) :: &
225 : nz
226 :
227 : ! Input Variables
228 : real( kind = core_rknd), dimension(nz), intent(in) :: &
229 : xm, & ! Mean of x (overall) [units vary]
230 : xp2, & ! Variance of x (overall) [(units vary)^2]
231 : Skx, & ! Skewness of x (overall) [-]
232 : mixt_frac ! Mixture fraction [-]
233 :
234 : ! Output Variables
235 : real( kind = core_rknd), dimension(nz), intent(out) :: &
236 : mu_x_1, & ! Mean of x (1st PDF component) [units vary]
237 : mu_x_2, & ! Mean of x (2nd PDF component) [units vary]
238 : sigma_x_1_sqd, & ! Variance of x (1st PDF component) [(units vary)^2]
239 : sigma_x_2_sqd ! Variance of x (2nd PDF component) [(units vary)^2]
240 :
241 : ! Local Variables
242 : real( kind = core_rknd), dimension(nz) :: &
243 0 : sgn_Skx, & ! Sign of Skx [-]
244 0 : B_x ! Spread of the PDF component means function [units vary]
245 :
246 :
247 : ! Find the sign of Skx
248 0 : where ( Skx >= zero )
249 : sgn_Skx = one
250 : elsewhere ! Skx < 0
251 : sgn_Skx = -one
252 : endwhere
253 :
254 : ! Calculate B_x, the LY function for the spread of the PDF component means.
255 : B_x = sgn_Skx * sqrt( xp2 ) &
256 0 : * ( abs( Skx ) / ( one - mixt_frac ) )**one_third
257 :
258 : ! Calculate the mean of x in the 1st PDF component.
259 0 : mu_x_1 = xm - B_x * ( one - mixt_frac )
260 :
261 : ! Calculate the mean of x in the 2nd PDF component.
262 0 : mu_x_2 = xm + B_x * mixt_frac
263 :
264 : ! Calculate the variance of x in the 1st PDF component.
265 : sigma_x_1_sqd = xp2 - B_x**2 * ( one - mixt_frac ) &
266 : * ( one + mixt_frac + mixt_frac**2 ) &
267 0 : / ( three * mixt_frac )
268 :
269 : ! Calculate the variance of x in the 2nd PDF component.
270 0 : sigma_x_2_sqd = xp2 + B_x**2 * ( one - mixt_frac )**2 / three
271 :
272 :
273 0 : return
274 :
275 : end subroutine calc_params_LY93
276 :
277 : !=============================================================================
278 :
279 : end module LY93_pdf
|