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