Line data Source code
1 : ! $Id$
2 : !===============================================================================
3 : module adg1_adg2_3d_luhar_pdf
4 :
5 : ! Description:
6 : ! Contains code to calculate the PDF parameters using the Analytic Double
7 : ! Gaussian 1 (ADG1) PDF closure, the Analytic Double Gaussian 2 (ADG2) PDF
8 : ! closure, or the 3D Luhar PDF closure. The ADG1 PDF is the original PDF used
9 : ! by CLUBB in interactive runs.
10 :
11 : ! References:
12 : !-------------------------------------------------------------------------
13 :
14 : implicit none
15 :
16 : public :: ADG1_pdf_driver, & ! Procedure(s)
17 : ADG2_pdf_driver, &
18 : Luhar_3D_pdf_driver, &
19 : ADG1_w_closure, &
20 : calc_Luhar_params, &
21 : close_Luhar_pdf
22 :
23 : private :: ADG1_ADG2_responder_params, & ! Procedure(s)
24 : backsolve_Luhar_params, &
25 : max_cubic_root
26 :
27 : private
28 :
29 : contains
30 :
31 : !=============================================================================
32 17870112 : subroutine ADG1_pdf_driver( nz, ngrdcol, sclr_dim, sclr_tol, & ! In
33 17870112 : wm, rtm, thlm, um, vm, & ! In
34 17870112 : wp2, rtp2, thlp2, up2, vp2, & ! In
35 17870112 : Skw, wprtp, wpthlp, upwp, vpwp, sqrt_wp2, & ! In
36 17870112 : sigma_sqd_w, beta, mixt_frac_max_mag, & ! In
37 17870112 : sclrm, sclrp2, wpsclrp, l_scalar_calc, & ! In
38 17870112 : w_1, w_2, & ! Out
39 17870112 : rt_1, rt_2, & ! Out
40 17870112 : thl_1, thl_2, & ! Out
41 17870112 : u_1, u_2, v_1, v_2, & ! Out
42 17870112 : varnce_w_1, varnce_w_2, & ! Out
43 17870112 : varnce_rt_1, varnce_rt_2, & ! Out
44 17870112 : varnce_thl_1, varnce_thl_2, & ! Out
45 17870112 : varnce_u_1, varnce_u_2, & ! Out
46 17870112 : varnce_v_1, varnce_v_2, & ! Out
47 17870112 : mixt_frac, & ! Out
48 17870112 : alpha_rt, alpha_thl, & ! Out
49 17870112 : alpha_u, alpha_v, & ! Out
50 17870112 : sclr_1, sclr_2, varnce_sclr_1, & ! Out
51 17870112 : varnce_sclr_2, alpha_sclr ) ! Out
52 :
53 : ! Calculates the PDF parameters using the Analytic Double Gaussian 1 (ADG1)
54 : ! PDF closure.
55 :
56 : ! References:
57 : !-----------------------------------------------------------------------
58 :
59 : use constants_clubb, only: &
60 : rt_tol, & ! Constant(s)
61 : thl_tol
62 :
63 : use clubb_precision, only: &
64 : core_rknd ! Variable(s)
65 :
66 : implicit none
67 :
68 : integer, intent(in) :: &
69 : nz, & ! Number of vertical levels
70 : ngrdcol, & ! Number of grid columns
71 : sclr_dim ! Number of passive scalars
72 :
73 : real( kind = core_rknd ), intent(in), dimension(sclr_dim) :: &
74 : sclr_tol ! Threshold(s) on the passive scalars [units vary]
75 :
76 : ! Input Variables
77 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
78 : wm, & ! Mean of w-wind comp. (vert. vel.) [m/s]
79 : rtm, & ! Mean of total water mixing ratio [kg/kg]
80 : thlm, & ! Mean of liquid water potential temp. [K]
81 : um, & ! Mean of eastward wind [m/s]
82 : vm, & ! Mean of northward wind [m/s]
83 : wp2, & ! Variance of w (overall) [m^2/s^2]
84 : rtp2, & ! Variance of r_t (overall) [(kg/kg)^2]
85 : thlp2, & ! Variance of th_l (overall) [K^2]
86 : up2, & ! Variance of eastward wind (overall) [m^2/s^2]
87 : vp2, & ! Variance of northward wind (overall) [m^2/s^2]
88 : Skw, & ! Skewness of w [-]
89 : wprtp, & ! Covariance of w and r_t [(kg/kg)(m/s)]
90 : wpthlp, & ! Covariance of w and th_l [K(m/s)]
91 : upwp, & ! Covariance of u and w [m^2/s^2]
92 : vpwp, & ! Covariance of v and w [m^2/s^2]
93 : sqrt_wp2, & ! Square root of variance of w [m/s]
94 : sigma_sqd_w ! Width of individual w plumes [-]
95 :
96 : real( kind = core_rknd ), dimension(ngrdcol), intent(in) :: &
97 : beta ! CLUBB tunable parameter beta [-]
98 :
99 : real( kind = core_rknd ), intent(in) :: &
100 : mixt_frac_max_mag ! Maximum allowable mag. of mixt_frac [-]
101 :
102 : real( kind = core_rknd ), dimension(ngrdcol, nz, sclr_dim), intent(in) :: &
103 : sclrm, & ! Mean of passive scalar (overall) [units vary]
104 : sclrp2, & ! Variance of passive scalar (overall) [(units vary)^2]
105 : wpsclrp ! Covariance of w and passive scalar [m/s (units vary)]
106 :
107 : logical, intent(in) :: &
108 : l_scalar_calc ! Flag to perform calculations for passive scalars
109 :
110 : ! Output Variables
111 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
112 : w_1, & ! Mean of w (1st PDF component) [m/s]
113 : w_2, & ! Mean of w (2nd PDF component) [m/s]
114 : rt_1, & ! Mean of r_t (1st PDF component) [kg/kg]
115 : rt_2, & ! Mean of r_t (2nd PDF component) [kg/kg]
116 : thl_1, & ! Mean of th_l (1st PDF component) [K]
117 : thl_2, & ! Mean of th_l (2nd PDF component) [K]
118 : u_1, & ! Mean of u (eastward wind) (1st PDF component) [m/s]
119 : u_2, & ! Mean of u (eastward wind) (2nd PDF component) [m/s]
120 : v_1, & ! Mean of v (northward wind) (1st PDF component) [m/s]
121 : v_2, & ! Mean of v (northward wind) (2nd PDF component) [m/s]
122 : varnce_w_1, & ! Variance of w (1st PDF component) [m^2/s^2]
123 : varnce_w_2, & ! Variance of w (2nd PDF component) [m^2/s^2]
124 : varnce_rt_1, & ! Variance of r_t (1st PDF component) [kg^2/kg^2]
125 : varnce_rt_2, & ! Variance of r_t (2nd PDF component) [kg^2/kg^2]
126 : varnce_thl_1, & ! Variance of th_l (1st PDF component) [K^2]
127 : varnce_thl_2, & ! Variance of th_l (2nd PDF component) [K^2]
128 : varnce_u_1, & ! Variance of u wind (1st PDF component) [m^2/s^2]
129 : varnce_u_2, & ! Variance of u wind (2nd PDF component) [m^2/s^2]
130 : varnce_v_1, & ! Variance of v wind (1st PDF component) [m^2/s^2]
131 : varnce_v_2, & ! Variance of v wind (2nd PDF component) [m^2/s^2]
132 : mixt_frac, & ! Mixture fraction (weight of 1st PDF component) [-]
133 : alpha_thl, & ! Factor relating to normalized variance for th_l [-]
134 : alpha_rt, & ! Factor relating to normalized variance for r_t [-]
135 : alpha_u, & ! Factor relating to normalized variance for u wind [-]
136 : alpha_v ! Factor relating to normalized variance for v wind [-]
137 :
138 : real( kind = core_rknd ), dimension(ngrdcol, nz, sclr_dim), intent(out) :: &
139 : sclr_1, & ! Mean of passive scalar (1st PDF component) [units vary]
140 : sclr_2, & ! Mean of passive scalar (2nd PDF component) [units vary]
141 : varnce_sclr_1, & ! Variance of pass. sclr (1st PDF comp.) [(units vary)^2]
142 : varnce_sclr_2, & ! Variance of pass. sclr (2nd PDF comp.) [(units vary)^2]
143 : alpha_sclr ! Factor relating to normalized variance for sclr [-]
144 :
145 : ! Local Variables
146 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
147 35740224 : w_1_n, & ! Normalized mean of w (1st PDF component) [-]
148 35740224 : w_2_n ! Normalized mean of w (2nd PDF component) [-]
149 :
150 : integer :: j ! Loop index
151 :
152 : !$acc enter data create( w_1_n, w_2_n )
153 :
154 : ! Calculate the mixture fraction and the PDF component means and variances
155 : ! of w.
156 : call ADG1_w_closure( nz, ngrdcol, wm, wp2, Skw, sigma_sqd_w, & ! In
157 : sqrt_wp2, mixt_frac_max_mag, & ! In
158 : w_1, w_2, w_1_n, w_2_n, & ! Out
159 17870112 : varnce_w_1, varnce_w_2, mixt_frac ) ! Out
160 :
161 : ! Calculate the PDF component means and variances of rt.
162 : call ADG1_ADG2_responder_params( nz, ngrdcol, &
163 : rtm, rtp2, wp2, sqrt_wp2, & ! In
164 : wprtp, w_1_n, w_2_n, mixt_frac, & ! In
165 : sigma_sqd_w, beta, rt_tol, & ! In
166 : rt_1, rt_2, varnce_rt_1, & ! Out
167 17870112 : varnce_rt_2, alpha_rt ) ! Out
168 :
169 : ! Calculate the PDF component means and variances of thl.
170 : call ADG1_ADG2_responder_params( nz, ngrdcol, &
171 : thlm, thlp2, wp2, sqrt_wp2, & ! In
172 : wpthlp, w_1_n, w_2_n, mixt_frac, & ! In
173 : sigma_sqd_w, beta, thl_tol, & ! In
174 : thl_1, thl_2, varnce_thl_1, & ! Out
175 17870112 : varnce_thl_2, alpha_thl ) ! Out
176 :
177 : ! Calculate the PDF component means and variances of u wind.
178 : call ADG1_ADG2_responder_params( nz, ngrdcol, &
179 : um, up2, wp2, sqrt_wp2, & ! In
180 : upwp, w_1_n, w_2_n, mixt_frac, & ! In
181 : sigma_sqd_w, beta, thl_tol, & ! In
182 : u_1, u_2, varnce_u_1, & ! Out
183 17870112 : varnce_u_2, alpha_u ) ! Out
184 :
185 : ! Calculate the PDF component means and variances of v wind.
186 : call ADG1_ADG2_responder_params( nz, ngrdcol, &
187 : vm, vp2, wp2, sqrt_wp2, & ! In
188 : vpwp, w_1_n, w_2_n, mixt_frac, & ! In
189 : sigma_sqd_w, beta, thl_tol, & ! In
190 : v_1, v_2, varnce_v_1, & ! Out
191 17870112 : varnce_v_2, alpha_v ) ! Out
192 :
193 : ! Calculate the PDF component means and variances of passive scalars.
194 17870112 : if ( l_scalar_calc ) then
195 0 : do j = 1, sclr_dim
196 : call ADG1_ADG2_responder_params( nz, ngrdcol, &
197 : sclrm(:,:,j), sclrp2(:,:,j), & ! In
198 : wp2, sqrt_wp2, wpsclrp(:,:,j), & ! In
199 : w_1_n, w_2_n, mixt_frac, & ! In
200 : sigma_sqd_w, beta, & ! In
201 0 : sclr_tol(j), & ! In
202 : sclr_1(:,:,j), sclr_2(:,:,j), & ! Out
203 : varnce_sclr_1(:,:,j), & ! Out
204 : varnce_sclr_2(:,:,j), & ! Out
205 0 : alpha_sclr(:,:,j) ) ! Out
206 : enddo ! i=1, sclr_dim
207 : endif ! l_scalar_calc
208 :
209 : !$acc exit data delete( w_1_n, w_2_n )
210 :
211 17870112 : return
212 :
213 : end subroutine ADG1_pdf_driver
214 :
215 : !=============================================================================
216 0 : subroutine ADG2_pdf_driver( nz, ngrdcol, sclr_dim, sclr_tol, & ! In
217 0 : wm, rtm, thlm, wp2, rtp2, thlp2, & ! In
218 0 : Skw, wprtp, wpthlp, sqrt_wp2, beta, & ! In
219 0 : sclrm, sclrp2, wpsclrp, l_scalar_calc, & ! In
220 0 : w_1, w_2, & ! Out
221 0 : rt_1, rt_2, & ! Out
222 0 : thl_1, thl_2, & ! Out
223 0 : varnce_w_1, varnce_w_2, & ! Out
224 0 : varnce_rt_1, varnce_rt_2, & ! Out
225 0 : varnce_thl_1, varnce_thl_2, & ! Out
226 0 : mixt_frac, & ! Out
227 0 : alpha_rt, alpha_thl, & ! Out
228 0 : sigma_sqd_w, sclr_1, sclr_2, & ! Out
229 0 : varnce_sclr_1, varnce_sclr_2, alpha_sclr ) ! Out
230 :
231 : ! Calculates the PDF parameters using the Analytic Double Gaussian 2 (ADG2)
232 : ! PDF closure.
233 :
234 : ! References:
235 : !-----------------------------------------------------------------------
236 :
237 : use grid_class, only: &
238 : grid ! Type
239 :
240 : use constants_clubb, only: &
241 : one, & ! Constant(s)
242 : w_tol_sqd, &
243 : rt_tol, &
244 : thl_tol
245 :
246 : use clubb_precision, only: &
247 : core_rknd ! Variable(s)
248 :
249 : implicit none
250 :
251 : integer, intent(in) :: &
252 : nz, & ! Number of vertical levels
253 : ngrdcol, & ! Number of grid columns
254 : sclr_dim ! Number of passive scalars
255 :
256 : real( kind = core_rknd ), intent(in), dimension(sclr_dim) :: &
257 : sclr_tol ! Threshold(s) on the passive scalars [units vary]
258 :
259 : ! Input Variables
260 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
261 : wm, & ! Mean of w-wind component (vertical velocity) [m/s]
262 : rtm, & ! Mean of total water mixing ratio [kg/kg]
263 : thlm, & ! Mean of liquid water potential temperature [K]
264 : wp2, & ! Variance of w (overall) [m^2/s^2]
265 : rtp2, & ! Variance of r_t (overall) [(kg/kg)^2]
266 : thlp2, & ! Variance of th_l (overall) [K^2]
267 : Skw, & ! Skewness of w [-]
268 : wprtp, & ! Covariance of w and r_t [(kg/kg)(m/s)]
269 : wpthlp, & ! Covariance of w and th_l [K(m/s)]
270 : sqrt_wp2 ! Square root of variance of w [m/s]
271 :
272 : real ( kind = core_rknd ), dimension(ngrdcol), intent(in) :: &
273 : beta ! CLUBB tunable parameter beta [-]
274 :
275 : real( kind = core_rknd ), dimension(ngrdcol,nz, sclr_dim), intent(in) :: &
276 : sclrm, & ! Mean of passive scalar (overall) [units vary]
277 : sclrp2, & ! Variance of passive scalar (overall) [(units vary)^2]
278 : wpsclrp ! Covariance of w and passive scalar [m/s (units vary)]
279 :
280 : logical, intent(in) :: &
281 : l_scalar_calc ! Flag to perform calculations for passive scalars
282 :
283 : ! Output Variables
284 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
285 : w_1, & ! Mean of w (1st PDF component) [m/s]
286 : w_2, & ! Mean of w (2nd PDF component) [m/s]
287 : rt_1, & ! Mean of r_t (1st PDF component) [kg/kg]
288 : rt_2, & ! Mean of r_t (2nd PDF component) [kg/kg]
289 : thl_1, & ! Mean of th_l (1st PDF component) [K]
290 : thl_2, & ! Mean of th_l (2nd PDF component) [K]
291 : varnce_w_1, & ! Variance of w (1st PDF component) [m^2/s^2]
292 : varnce_w_2, & ! Variance of w (2nd PDF component) [m^2/s^2]
293 : varnce_rt_1, & ! Variance of r_t (1st PDF component) [kg^2/kg^2]
294 : varnce_rt_2, & ! Variance of r_t (2nd PDF component) [kg^2/kg^2]
295 : varnce_thl_1, & ! Variance of th_l (1st PDF component) [K^2]
296 : varnce_thl_2, & ! Variance of th_l (2nd PDF component) [K^2]
297 : mixt_frac, & ! Mixture fraction (weight of 1st PDF component) [-]
298 : alpha_thl, & ! Factor relating to normalized variance for th_l [-]
299 : alpha_rt, & ! Factor relating to normalized variance for r_t [-]
300 : sigma_sqd_w ! Width of individual w plumes (ADG1 parameter) [-]
301 :
302 : real( kind = core_rknd ), dimension(ngrdcol,nz, sclr_dim), intent(out) :: &
303 : sclr_1, & ! Mean of passive scalar (1st PDF component) [units vary]
304 : sclr_2, & ! Mean of passive scalar (2nd PDF component) [units vary]
305 : varnce_sclr_1, & ! Variance of pass. sclr (1st PDF comp.) [(units vary)^2]
306 : varnce_sclr_2, & ! Variance of pass. sclr (2nd PDF comp.) [(units vary)^2]
307 : alpha_sclr ! Factor relating to normalized variance for sclr [-]
308 :
309 : ! Local Variables
310 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
311 0 : w_1_n, & ! Normalized mean of w (1st PDF component) [-]
312 0 : w_2_n, & ! Normalized mean of w (2nd PDF component) [-]
313 0 : small_m_w, & ! Luhar's m (tunable parameter) [-]
314 0 : big_m_w, & ! Luhar's M (a function of Luhar's m) [-]
315 0 : sigma_sqd_w_1, & ! Normalized width parameter of w (1st PDF component) [-]
316 0 : sigma_sqd_w_2 ! Normalized width parameter of w (2nd PDF component) [-]
317 :
318 : integer :: j, i ! Loop index
319 :
320 :
321 : ! Calculate the mixture fraction and the PDF component means and variances
322 : ! of w.
323 :
324 : ! Reproduce ADG2_w_closure using separate functions
325 0 : do i = 1, ngrdcol
326 :
327 0 : call calc_Luhar_params( nz, Skw(i,:), wp2(i,:), & ! intent(in)
328 0 : wp2(i,:), w_tol_sqd, & ! intent(in)
329 0 : mixt_frac(i,:), big_m_w(i,:), small_m_w(i,:) ) ! intent(out)
330 :
331 0 : call close_Luhar_pdf( nz, wm(i,:), wp2(i,:), mixt_frac(i,:), & ! intent(in)
332 0 : small_m_w(i,:), wp2(i,:), w_tol_sqd, & ! intent(in)
333 0 : sigma_sqd_w_1(i,:), sigma_sqd_w_2(i,:), & ! intent(out)
334 0 : varnce_w_1(i,:), varnce_w_2(i,:), & ! intent(out)
335 0 : w_1_n(i,:), w_2_n(i,:), w_1(i,:), w_2(i,:) ) ! intent(out)
336 :
337 : ! Overwrite sigma_sqd_w for consistency with ADG1
338 0 : sigma_sqd_w(i,:) = min( one / ( one + small_m_w(i,:)**2 ), 0.99_core_rknd )
339 :
340 : end do
341 :
342 : ! Calculate the PDF component means and variances of rt.
343 : call ADG1_ADG2_responder_params( nz, ngrdcol, & ! In
344 : rtm, rtp2, wp2, sqrt_wp2, & ! In
345 : wprtp, w_1_n, w_2_n, mixt_frac, & ! In
346 : sigma_sqd_w, beta, rt_tol, & ! In
347 : rt_1, rt_2, varnce_rt_1, & ! Out
348 0 : varnce_rt_2, alpha_rt ) ! Out
349 :
350 : ! Calculate the PDF component means and variances of thl.
351 : call ADG1_ADG2_responder_params( nz, ngrdcol, & ! In
352 : thlm, thlp2, wp2, sqrt_wp2, & ! In
353 : wpthlp, w_1_n, w_2_n, mixt_frac, & ! In
354 : sigma_sqd_w, beta, thl_tol, & ! In
355 : thl_1, thl_2, varnce_thl_1, & ! Out
356 0 : varnce_thl_2, alpha_thl ) ! Out
357 :
358 : ! Calculate the PDF component means and variances of passive scalars.
359 0 : if ( l_scalar_calc ) then
360 0 : do j = 1, sclr_dim
361 : call ADG1_ADG2_responder_params( nz, ngrdcol, & ! In
362 : sclrm(:,:,j), sclrp2(:,:,j), & ! In
363 : wp2, sqrt_wp2, wpsclrp(:,:,j), & ! In
364 : w_1_n, w_2_n, mixt_frac, & ! In
365 : sigma_sqd_w, beta, & ! In
366 0 : sclr_tol(j), & ! In
367 : sclr_1(:,:,j), sclr_2(:,:,j), & ! Out
368 : varnce_sclr_1(:,:,j), & ! Out
369 : varnce_sclr_2(:,:,j), & ! Out
370 0 : alpha_sclr(:,:,j) ) ! Out
371 : enddo ! i=1, sclr_dim
372 : endif ! l_scalar_calc
373 :
374 0 : return
375 :
376 : end subroutine ADG2_pdf_driver
377 :
378 : !=============================================================================
379 0 : subroutine Luhar_3D_pdf_driver( nz, wm, rtm, thlm, wp2, rtp2, thlp2, & ! In
380 0 : Skw, Skrt, Skthl, wprtp, wpthlp, & ! In
381 0 : w_1, w_2, & ! Out
382 0 : rt_1, rt_2, & ! Out
383 0 : thl_1, thl_2, & ! Out
384 0 : varnce_w_1, varnce_w_2, & ! Out
385 0 : varnce_rt_1, varnce_rt_2, & ! Out
386 0 : varnce_thl_1, varnce_thl_2, & ! Out
387 0 : mixt_frac ) ! Out
388 :
389 : ! Calculates the PDF parameters using the 3D Luhar PDF closure.
390 :
391 : ! References:
392 : !-----------------------------------------------------------------------
393 :
394 : use constants_clubb, only: &
395 : w_tol_sqd, & ! Constant(s)
396 : rt_tol, &
397 : thl_tol
398 :
399 : use clubb_precision, only: &
400 : core_rknd ! Variable(s)
401 :
402 : implicit none
403 :
404 : integer, intent(in) :: &
405 : nz
406 :
407 : ! Input Variables
408 : real( kind = core_rknd ), dimension(nz), intent(in) :: &
409 : wm, & ! Mean of w-wind component (vertical velocity) [m/s]
410 : rtm, & ! Mean of total water mixing ratio [kg/kg]
411 : thlm, & ! Mean of liquid water potential temperature [K]
412 : wp2, & ! Variance of w (overall) [m^2/s^2]
413 : rtp2, & ! Variance of r_t (overall) [(kg/kg)^2]
414 : thlp2, & ! Variance of th_l (overall) [K^2]
415 : Skw, & ! Skewness of w [-]
416 : Skrt, & ! Skewness of rt [-]
417 : Skthl, & ! Skewness of th_l [-]
418 : wprtp, & ! Covariance of w and r_t [(kg/kg)(m/s)]
419 : wpthlp ! Covariance of w and th_l [K(m/s)]
420 :
421 : ! Output Variables
422 : real( kind = core_rknd ), dimension(nz), intent(out) :: &
423 : w_1, & ! Mean of w (1st PDF component) [m/s]
424 : w_2, & ! Mean of w (2nd PDF component) [m/s]
425 : rt_1, & ! Mean of r_t (1st PDF component) [kg/kg]
426 : rt_2, & ! Mean of r_t (2nd PDF component) [kg/kg]
427 : thl_1, & ! Mean of th_l (1st PDF component) [K]
428 : thl_2, & ! Mean of th_l (2nd PDF component) [K]
429 : varnce_w_1, & ! Variance of w (1st PDF component) [m^2/s^2]
430 : varnce_w_2, & ! Variance of w (2nd PDF component) [m^2/s^2]
431 : varnce_rt_1, & ! Variance of r_t (1st PDF component) [kg^2/kg^2]
432 : varnce_rt_2, & ! Variance of r_t (2nd PDF component) [kg^2/kg^2]
433 : varnce_thl_1, & ! Variance of th_l (1st PDF component) [K^2]
434 : varnce_thl_2, & ! Variance of th_l (2nd PDF component) [K^2]
435 : mixt_frac ! Mixture fraction (weight of 1st PDF component) [-]
436 :
437 : real( kind = core_rknd ), dimension(nz) :: &
438 0 : w_1_n, & ! Normalized mean of w (1st PDF component) [-]
439 0 : w_2_n, & ! Normalized mean of w (2nd PDF component) [-]
440 0 : rt_1_n, & ! Normalized mean of rt (1st PDF component) [-]
441 0 : rt_2_n, & ! Normalized mean of rt (2nd PDF component) [-]
442 0 : thl_1_n, & ! Normalized mean of thl (1st PDF component) [-]
443 0 : thl_2_n, & ! Normalized mean of thl (2nd PDF component) [-]
444 0 : small_m_w, & ! Luhar's m (tunable parameter) for w [-]
445 0 : big_m_w, & ! Luhar's M (a function of Luhar's m) for w [-]
446 0 : small_m_rt, & ! Luhar's m (tunable parameter) for rt [-]
447 0 : big_m_rt, & ! Luhar's M (a function of Luhar's m) for rt [-]
448 0 : small_m_thl, & ! Luhar's m (tunable parameter) for thl [-]
449 0 : big_m_thl, & ! Luhar's M (a function of Luhar's m) for thl [-]
450 0 : sigma_sqd_w_1, & ! Normalized width parameter of w (1st PDF comp.) [-]
451 0 : sigma_sqd_w_2, & ! Normalized width parameter of w (2nd PDF comp.) [-]
452 0 : sigma_sqd_rt_1, & ! Normalized width parameter of rt (1st PDF comp.) [-]
453 0 : sigma_sqd_rt_2, & ! Normalized width parameter of rt (2nd PDF comp.) [-]
454 0 : sigma_sqd_thl_1, & ! Normalized width parameter of thl (1st PDF comp.) [-]
455 0 : sigma_sqd_thl_2 ! Normalized width parameter of thl (2nd PDF comp.) [-]
456 :
457 : integer :: k ! Vertical loop index
458 :
459 :
460 0 : do k = 1, nz, 1
461 :
462 0 : if ( ( abs( Skw(k) ) >= abs( Skthl(k) ) ) &
463 0 : .and. ( abs( Skw(k) ) >= abs( Skrt(k) ) ) ) then
464 :
465 : ! w has the greatest magnitude of skewness.
466 :
467 : ! Solve for the w PDF
468 : call calc_Luhar_params( 1, Skw(k), wp2(k), & ! In
469 : wp2(k), w_tol_sqd, & ! In
470 0 : mixt_frac(k), big_m_w(k), small_m_w(k) ) ! Out
471 :
472 : call close_Luhar_pdf( 1, wm(k), wp2(k), mixt_frac(k), & ! In
473 : small_m_w(k), wp2(k), w_tol_sqd, & ! In
474 : sigma_sqd_w_1(k), sigma_sqd_w_2(k), & ! Out
475 : varnce_w_1(k), varnce_w_2(k), & ! Out
476 0 : w_1_n(k), w_2_n(k), w_1(k), w_2(k) ) ! Out
477 :
478 : ! Solve for the thl PDF
479 : call backsolve_Luhar_params( Skw(k), Skthl(k), & ! In
480 : big_m_w(k), mixt_frac(k), & ! In
481 : big_m_thl(k), small_m_thl(k) ) ! Out
482 :
483 : call close_Luhar_pdf( 1, thlm(k), thlp2(k), mixt_frac(k), & ! In
484 : small_m_thl(k), wpthlp(k), thl_tol**2, & ! In
485 : sigma_sqd_thl_1(k), sigma_sqd_thl_2(k), & ! Out
486 : varnce_thl_1(k), varnce_thl_2(k), & ! Out
487 : thl_1_n(k), thl_2_n(k), & ! Out
488 0 : thl_1(k), thl_2(k) ) ! Out
489 :
490 : ! Solve for the rt PDF
491 : call backsolve_Luhar_params( Skw(k), Skrt(k), & ! In
492 : big_m_w(k), mixt_frac(k), & ! In
493 : big_m_rt(k), small_m_rt(k) ) ! Out
494 :
495 : call close_Luhar_pdf( 1, rtm(k), rtp2(k), mixt_frac(k), & ! In
496 : small_m_rt(k), wprtp(k), rt_tol**2, & ! In
497 : sigma_sqd_rt_1(k), sigma_sqd_rt_2(k), & ! Out
498 : varnce_rt_1(k), varnce_rt_2(k), & ! Out
499 0 : rt_1_n(k), rt_2_n(k), rt_1(k), rt_2(k) ) ! Out
500 :
501 : elseif ( ( abs( Skthl(k) ) > abs( Skw(k) ) ) &
502 0 : .and. ( abs( Skthl(k) ) >= abs( Skrt(k) ) ) ) then
503 :
504 : ! theta-l has the greatest magnitude of skewness.
505 :
506 : ! Solve for the thl PDF
507 : call calc_Luhar_params( 1, Skthl(k), wpthlp(k), & ! In
508 : thlp2(k), thl_tol**2, & ! In
509 : mixt_frac(k), big_m_thl(k), & ! Out
510 0 : small_m_thl(k) ) ! Out
511 :
512 : ! Solve for the thl PDF
513 : call close_Luhar_pdf( 1, thlm(k), thlp2(k), mixt_frac(k), & ! In
514 : small_m_thl(k), wpthlp(k), thl_tol**2, & ! In
515 : sigma_sqd_thl_1(k), sigma_sqd_thl_2(k), & ! Out
516 : varnce_thl_1(k), varnce_thl_2(k), & ! Out
517 : thl_1_n(k), thl_2_n(k), & ! Out
518 0 : thl_1(k), thl_2(k) ) ! Out
519 :
520 : ! Solve for the w PDF
521 : call backsolve_Luhar_params( Skthl(k), Skw(k), & ! In
522 : big_m_thl(k), mixt_frac(k), & ! In
523 : big_m_w(k), small_m_w(k) ) ! Out
524 :
525 : call close_Luhar_pdf( 1, wm(k), wp2(k), mixt_frac(k), & ! In
526 : small_m_w(k), wp2(k), w_tol_sqd, & ! In
527 : sigma_sqd_w_1(k), sigma_sqd_w_2(k), & ! Out
528 : varnce_w_1(k), varnce_w_2(k), & ! Out
529 0 : w_1_n(k), w_2_n(k), w_1(k), w_2(k) ) ! Out
530 :
531 : ! Solve for the rt PDF
532 : call backsolve_Luhar_params( Skthl(k), Skrt(k), & ! In
533 : big_m_thl(k), mixt_frac(k), & ! In
534 : big_m_rt(k), small_m_rt(k) ) ! Out
535 :
536 : call close_Luhar_pdf( 1, rtm(k), rtp2(k), mixt_frac(k), & ! In
537 : small_m_rt(k), wprtp(k), rt_tol**2, & ! In
538 : sigma_sqd_rt_1(k), sigma_sqd_rt_2(k), & ! Out
539 : varnce_rt_1(k), varnce_rt_2(k), & ! Out
540 0 : rt_1_n(k), rt_2_n(k), rt_1(k), rt_2(k) ) ! Out
541 :
542 : else
543 :
544 : ! rt has the greatest magnitude of skewness.
545 :
546 : ! Solve for the rt PDF
547 : call calc_Luhar_params( 1, Skrt(k), wprtp(k), & ! In
548 : rtp2(k), rt_tol**2, & ! In
549 : mixt_frac(k), big_m_rt(k), & ! Out
550 0 : small_m_rt(k) ) ! Out
551 :
552 : ! Solve for the rt PDF
553 : call close_Luhar_pdf( 1, rtm(k), rtp2(k), mixt_frac(k), & ! In
554 : small_m_rt(k), wprtp(k), rt_tol**2, & ! In
555 : sigma_sqd_rt_1(k), sigma_sqd_rt_2(k), & ! Out
556 : varnce_rt_1(k), varnce_rt_2(k), & ! Out
557 0 : rt_1_n(k), rt_2_n(k), rt_1(k), rt_2(k) ) ! Out
558 :
559 : ! Solve for the w PDF
560 : call backsolve_Luhar_params( Skrt(k), Skw(k), & ! In
561 : big_m_rt(k), mixt_frac(k), & ! In
562 : big_m_w(k), small_m_w(k) ) ! Out
563 :
564 : call close_Luhar_pdf( 1, wm(k), wp2(k), mixt_frac(k), & ! In
565 : small_m_w(k), wp2(k), w_tol_sqd, & ! In
566 : sigma_sqd_w_1(k), sigma_sqd_w_2(k), & ! Out
567 : varnce_w_1(k), varnce_w_2(k), & ! Out
568 0 : w_1_n(k), w_2_n(k), w_1(k), w_2(k) ) ! OUt
569 :
570 : ! Solve for the thl PDF
571 : call backsolve_Luhar_params( Skrt(k), Skthl(k), & ! In
572 : big_m_rt(k), mixt_frac(k), & ! In
573 : big_m_thl(k), small_m_thl(k) ) ! Out
574 :
575 : call close_Luhar_pdf( 1, thlm(k), thlp2(k), mixt_frac(k), & ! In
576 : small_m_thl(k), wpthlp(k), thl_tol**2, & ! In
577 : sigma_sqd_thl_1(k), sigma_sqd_thl_2(k), & ! Out
578 : varnce_thl_1(k), varnce_thl_2(k), & ! Out
579 : thl_1_n(k), thl_2_n(k), & ! Out
580 0 : thl_1(k), thl_2(k) ) ! Out
581 :
582 : endif
583 :
584 : enddo ! k = 1, nz, 1
585 :
586 :
587 0 : return
588 :
589 : end subroutine Luhar_3D_pdf_driver
590 :
591 : !======================================================================
592 17870112 : subroutine ADG1_w_closure( nz, ngrdcol, wm, wp2, Skw, sigma_sqd_w, & !In
593 17870112 : sqrt_wp2, mixt_frac_max_mag, & !In
594 17870112 : w_1, w_2, w_1_n, w_2_n, & !Out
595 17870112 : varnce_w_1, varnce_w_2, mixt_frac ) !Out
596 :
597 :
598 : ! Description:
599 : ! Calculates the mixture fraction, the PDF component means of w, and the PDF
600 : ! component variances of w for the Analytic Double Gaussian 1 (ADG1)
601 : ! closure. It sets the widths of both w Gaussians to be the same, and
602 : ! furthermore, relates them to the overall variance of w (<w'^2>) by a
603 : ! parameter, sigma_sqd_w. The equation is:
604 : !
605 : ! sigma_w_1^2 = sigma_w_2^2 = sigma_sqd_w * <w'^2>;
606 : !
607 : ! where sigma_w_1^2 is the variance of w in the 1st PDF component,
608 : ! sigma_w_2^2 is the variance of w in the 2nd PDF component, and parameter
609 : ! sigma_sqd_w must have a value within the range 0 <= sigma_sqd_w <= 1.
610 : !
611 : ! References:
612 : ! Golaz, J-C., V. E. Larson, and W. R. Cotton, 2002a: A PDF-based model for
613 : ! boundary layer clouds. Part I: Method and model description. J. Atmos.
614 : ! Sci., 59, 3540–3551.
615 : !
616 : ! Vincent E. Larson and Jean-Christophe Golaz, 2005: Using Probability
617 : ! Density Functions to Derive Consistent Closure Relationships among
618 : ! Higher-Order Moments. Mon. Wea. Rev., 133, 1023–1042.
619 : !-----------------------------------------------------------------------
620 :
621 : use constants_clubb, only: &
622 : four, & ! Constant(s)
623 : one, &
624 : one_half, &
625 : zero, &
626 : w_tol_sqd
627 :
628 : use clubb_precision, only: &
629 : core_rknd ! Precision
630 :
631 : implicit none
632 :
633 : integer, intent(in) :: &
634 : ngrdcol, & ! Number of grid columns
635 : nz ! Number of vertical level
636 :
637 : ! Input Variables
638 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
639 : wm, & ! Mean of w (overall) [m/s]
640 : wp2, & ! Variance of w (overall) [m^2/s^2]
641 : Skw, & ! Skewness of w [-]
642 : sigma_sqd_w, & ! Widths of each w Gaussian [-]
643 : sqrt_wp2 ! Square root of the variance of w [m/s]
644 :
645 : real( kind = core_rknd ), intent(in) :: &
646 : mixt_frac_max_mag ! Maximum allowable mag. of mixt_frac [-]
647 :
648 : ! Output Variables
649 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
650 : w_1, & ! Mean of w (1st PDF component) [m/s]
651 : w_2, & ! Mean of w (2nd PDF component) [m/s]
652 : w_1_n, & ! Normalized mean of w (1st PDF component) [-]
653 : w_2_n, & ! Normalized mean of w (2nd PDF component) [-]
654 : varnce_w_1, & ! Variance of w (1st PDF component) [m^2/s^2]
655 : varnce_w_2, & ! Variance of w (2nd PDF component) [m^2/s^2]
656 : mixt_frac ! Mixture fraction [-]
657 :
658 : ! Local Variable
659 : integer :: k, i
660 :
661 : !----- Begin Code -----
662 :
663 : !$acc parallel loop gang vector collapse(2) default(present)
664 1536829632 : do k = 1, nz
665 25380961632 : do i = 1, ngrdcol
666 :
667 25363091520 : if ( wp2(i,k) > w_tol_sqd ) then
668 :
669 : ! Width (standard deviation) parameters are non-zero
670 :
671 : ! The variable "mixt_frac" is the weight of the 1st PDF component. The
672 : ! weight of the 2nd PDF component is "1-mixt_frac". If there isn't any
673 : ! skewness of w (Sk_w = 0 because w'^3 = 0), mixt_frac = 0.5, and both
674 : ! PDF components are equally weighted. If there is positive skewness of
675 : ! w (Sk_w > 0 because w'^3 > 0), 0 < mixt_frac < 0.5, and the 2nd PDF
676 : ! component has greater weight than does the 1st PDF component. If there
677 : ! is negative skewness of w (Sk_w < 0 because w'^3 < 0),
678 : ! 0.5 < mixt_frac < 1, and the 1st PDF component has greater weight than
679 : ! does the 2nd PDF component.
680 5897779378 : if ( abs( Skw(i,k) ) <= 1.0e-5_core_rknd ) then
681 442655339 : mixt_frac(i,k) = one_half
682 : else
683 : mixt_frac(i,k) &
684 : = one_half &
685 5455124039 : * ( one - Skw(i,k) / sqrt( four * ( one - sigma_sqd_w(i,k) )**3 + Skw(i,k)**2 ) )
686 : endif
687 :
688 : ! Clip mixt_frac, and 1 - mixt_frac, to avoid dividing by a small number.
689 : ! Formula for mixt_frac_max_mag =
690 : ! 1 - ( 1/2 * ( 1 - Skw_max
691 : ! / sqrt( 4*( 1 - sigma_sqd_w )^3 + Skw_max^2 ) ) ),
692 : ! where sigma_sqd_w is fixed at 0.4 for this calculation.
693 : mixt_frac(i,k) = min( max( mixt_frac(i,k), one - mixt_frac_max_mag ), &
694 5897779378 : mixt_frac_max_mag )
695 :
696 : ! The normalized mean of w for Gaussian "plume" 1 is w_1_n. It's value
697 : ! will always be greater than 0. As an example, a value of 1.0 would
698 : ! indicate that the actual mean of w for Gaussian "plume" 1 is found 1.0
699 : ! standard deviation above the overall mean for w.
700 : w_1_n(i,k) &
701 5897779378 : = sqrt( ( ( one - mixt_frac(i,k) ) / mixt_frac(i,k) ) * ( one - sigma_sqd_w(i,k) ) )
702 : ! The normalized mean of w for Gaussian "plume" 2 is w_2_n. It's value
703 : ! will always be less than 0. As an example, a value of -0.5 would
704 : ! indicate that the actual mean of w for Gaussian "plume" 2 is found 0.5
705 : ! standard deviations below the overall mean for w.
706 : w_2_n(i,k) &
707 5897779378 : = -sqrt( ( mixt_frac(i,k) / ( one - mixt_frac(i,k) ) ) * ( one - sigma_sqd_w(i,k) ) )
708 : ! The mean of w for Gaussian "plume" 1 is w_1.
709 5897779378 : w_1(i,k) = wm(i,k) + sqrt_wp2(i,k) * w_1_n(i,k)
710 : ! The mean of w for Gaussian "plume" 2 is w_2.
711 5897779378 : w_2(i,k) = wm(i,k) + sqrt_wp2(i,k) * w_2_n(i,k)
712 :
713 : ! The variance of w for Gaussian "plume" 1 for varnce_w_1.
714 5897779378 : varnce_w_1(i,k) = sigma_sqd_w(i,k) * wp2(i,k)
715 : ! The variance of w for Gaussian "plume" 2 for varnce_w_2.
716 : ! The variance in both Gaussian "plumes" is defined to be the same.
717 5897779378 : varnce_w_2(i,k) = sigma_sqd_w(i,k) * wp2(i,k)
718 :
719 : else
720 :
721 : ! Vertical velocity doesn't vary.
722 17946352622 : mixt_frac(i,k) = one_half
723 17946352622 : w_1_n(i,k) = sqrt( one - sigma_sqd_w(i,k) )
724 17946352622 : w_2_n(i,k) = -sqrt( one - sigma_sqd_w(i,k) )
725 17946352622 : w_1(i,k) = wm(i,k)
726 17946352622 : w_2(i,k) = wm(i,k)
727 17946352622 : varnce_w_1(i,k) = zero
728 17946352622 : varnce_w_2(i,k) = zero
729 :
730 : endif ! Widths non-zero
731 : end do
732 : end do
733 : !$acc end parallel loop
734 17870112 : return
735 :
736 :
737 :
738 : end subroutine ADG1_w_closure
739 :
740 :
741 : !=============================================================================
742 0 : subroutine calc_Luhar_params( nz, Skx, wpxp, & ! In
743 0 : xp2, x_tol_sqd, & ! In
744 0 : mixt_frac, big_m, small_m ) ! Out
745 :
746 : ! Description:
747 : ! For the Luhar closure, this subroutine takes Skx (and Skw) as input and
748 : ! outputs the mixture fraction, big_m, and small_m. This code was written
749 : ! using the equations and nomenclature of Larson et al. (2002) Appendix
750 : ! section e.
751 : !
752 : ! The relationship between skewness of x (Skx), mixture fraction (a), and
753 : ! Luhar's small m (m) is given by:
754 : !
755 : ! Skx^2 = ( m^2 * ( m^2 + 3 )^2 / ( m^2 + 1 )^3 )
756 : ! * ( 1 - 2*a )^2 / ( a * ( 1 - a ) ).
757 : !
758 : ! Luhar's large M (M) is used to more easily express the factor involving
759 : ! the m's:
760 : !
761 : ! M = ( m^2 + 1 )^3 / ( m^2 * ( m^2 + 3 )^2 ).
762 : !
763 : ! The equation involving skewness of x becomes:
764 : !
765 : ! Skx^2 = ( 1 / M ) * ( 1 - 2*a )^2 / ( a * ( 1 - a ) );
766 : !
767 : ! or:
768 : !
769 : ! M * Skx^2 = ( 1 - 2*a )^2 / ( a * ( 1 - a ) ).
770 : !
771 : ! This equation can be rewritten as:
772 : !
773 : ! ( a * ( 1 - a ) ) * M * Skx^2 = ( 1 - 2*a )^2;
774 : !
775 : ! as well as:
776 : !
777 : ! ( a - a^2 ) * M * Skx^2 = 1 - 4*a + 4*a^2;
778 : !
779 : ! and eventually as:
780 : !
781 : ! ( 4 + M * Skx^2 ) * a^2 - ( 4 + M * Skx^2 ) * a + 1 = 0.
782 : !
783 : ! Solving the quadratic equation for a:
784 : !
785 : ! a = (1/2) * ( 1 +- Skx * sqrt( 1 / ( 4/M + Skx^2 ) ) ).
786 : !
787 : ! Since by definition, mu_w_1 >= mu_w_2, a < 0.5 when Skw > 0, the equation
788 : ! for mixture fraction is:
789 : !
790 : ! a = (1/2) * ( 1 - Skx * sqrt( 1 / ( 4/M + Skx^2 ) ) ).
791 : !
792 : ! For 3-D Luhar, the variable (w, rt, or theta-l) with the greatest
793 : ! magnitude of skewness is used to calculate mixture fraction. Since it is
794 : ! desirable to still have a < 0.5 when Skw > 0 and a > 0.5 when Skw < 0, the
795 : ! sign function is used. The value of Skx is replaced by:
796 : !
797 : ! Skx|_adj = sgn( <w'x'> ) * Skx;
798 : !
799 : ! where
800 : !
801 : ! sgn( <w'x'> ) = | 1 when <w'x'> >= 0
802 : ! | -1 when <w'x'> < 0.
803 : !
804 : ! Since Skx|_adj^2 = ( sgn( <w'x'> ) * Skx )^2 = ( sgn( <w'x'> ) )^2 * Skx^2
805 : ! = Skx^2, the equation for mixture fraction is:
806 : !
807 : ! a = (1/2) * ( 1 - sgn( <w'x'> ) * Skx * sqrt( 1 / ( 4/M + Skx^2 ) ) ).
808 : !
809 : ! When using the ADG2 closure or when using the 3-D Luhar closure when the
810 : ! variable with the greatest magnitude of skewness is w, Skw = Skx and
811 : ! sgn( <w'^2> ) is always equal to 1, reducing the equation to its previous
812 : ! form.
813 :
814 : ! References:
815 : ! Vincent E. Larson, Jean-Christophe Golaz, and William R. Cotton, 2002:
816 : ! Small-Scale and Mesoscale Variability in Cloudy Boundary Layers: Joint
817 : ! Probability Density Functions. J. Atmos. Sci., 59, 3519–3539.
818 : !-----------------------------------------------------------------------
819 :
820 : use constants_clubb, only: &
821 : four, & ! Constant(s)
822 : three, &
823 : one, &
824 : two_thirds, &
825 : one_half, &
826 : one_third, &
827 : zero
828 :
829 : use clubb_precision, only: &
830 : core_rknd ! Precision
831 :
832 : implicit none
833 :
834 : ! Input Variables
835 : integer, intent(in) :: &
836 : nz ! Number of vertical levels
837 :
838 : real( kind = core_rknd ), dimension(nz), intent(in) :: &
839 : Skx, & ! Skewness of x ( <x'^3> / <x'2>^(3/2) ) [-]
840 : wpxp, & ! Covariance of w and x [m/s (x units)]
841 : xp2 ! Variance of x [(x units)^2]
842 :
843 : real( kind = core_rknd ), intent(in) :: &
844 : x_tol_sqd ! Tolerance value of x squared [(x units)^2]
845 :
846 : ! Output Variables
847 : real( kind = core_rknd ), dimension(nz), intent(out) :: &
848 : mixt_frac, & ! Mixture fraction [-]
849 : big_m, & ! Luhar's M [-]
850 : small_m ! Luhar's m [-]
851 :
852 : ! Local Variables
853 : real( kind = core_rknd ), dimension(nz) :: &
854 0 : small_m_sqd, & ! Luhar's m^2 [-]
855 0 : sgn_wpxp ! Sgn(<w'x'>); 1 when <w'x'> >= 0 or -1 when <w'x'> < 0 [-]
856 :
857 :
858 0 : where ( xp2 > x_tol_sqd )
859 :
860 : ! Width (standard deviation) parameters are non-zero
861 :
862 : ! Calculate Luhar's m (small m).
863 : ! If Skx is very small, then small_m will tend to zero which risks
864 : ! divide-by-zero. To ameliorate this problem, we enforce abs( x_1_n )
865 : ! and abs( x_2_n ) > 0.05.
866 : ! Note: Luhar's small_m (m) is the only tunable parameter in the Luhar
867 : ! closure, so this equation can be changed. However, the value
868 : ! of m should go toward 0 as Skx goes toward 0 so that the double
869 : ! Gaussian reduces to a single Gaussian when the distribution is
870 : ! unskewed.
871 : small_m = max( two_thirds * abs( Skx )**one_third, 0.05_core_rknd )
872 :
873 : ! Calculate m^2.
874 : small_m_sqd = small_m**2
875 :
876 : ! Calculate Luhar's M (big M).
877 : big_m = ( one + small_m_sqd )**3 &
878 : / ( ( three + small_m_sqd )**2 * small_m_sqd )
879 :
880 : ! Calculate sgn( <w'x'> ).
881 : where ( wpxp >= zero )
882 : sgn_wpxp = one
883 : elsewhere ! <w'x'> < 0
884 : sgn_wpxp = -one
885 : endwhere ! <w'x'> >= 0
886 :
887 : ! Calculate mixture fraction.
888 : mixt_frac = one_half &
889 : * ( one - sgn_wpxp * Skx &
890 : * sqrt( one / ( ( four / big_m ) + Skx**2 ) ) )
891 :
892 : elsewhere
893 :
894 : ! Variable x doesn't vary.
895 : mixt_frac = one_half
896 :
897 : ! For output purposes, set small_m and big_m to 0.
898 : small_m = zero
899 : big_m = zero
900 :
901 : endwhere ! Widths non-zero
902 :
903 :
904 0 : return
905 :
906 : end subroutine calc_Luhar_params
907 :
908 : !=============================================================================
909 0 : subroutine close_Luhar_pdf( nz, xm, xp2, mixt_frac, & ! In
910 0 : small_m, wpxp, x_tol_sqd, & ! In
911 0 : sigma_sqd_x_1, sigma_sqd_x_2, & ! Out
912 0 : varnce_x_1, varnce_x_2, & ! Out
913 0 : x_1_n, x_2_n, x_1, x_2 ) ! Out
914 :
915 : ! Description:
916 : ! For the Luhar closure, this subroutine takes Skx, xm, xp2, and mixt_frac,
917 : ! big_m, and small_m (calculated in calc_Luhar_params) as input and outputs
918 : ! the PDF component means and variances of a variable x in the joint-PDF
919 : ! according to Luhar et al. (1996). This code was written using the
920 : ! equations and nomenclature of Larson et al. (2002) Appendix section e.
921 :
922 : ! References:
923 : ! Vincent E. Larson, Jean-Christophe Golaz, and William R. Cotton, 2002:
924 : ! Small-Scale and Mesoscale Variability in Cloudy Boundary Layers: Joint
925 : ! Probability Density Functions. J. Atmos. Sci., 59, 3519–3539.
926 : !-----------------------------------------------------------------------
927 :
928 : use constants_clubb, only: &
929 : one, & ! Constant(s)
930 : zero
931 :
932 : use clubb_precision, only: &
933 : core_rknd ! Precision
934 :
935 : implicit none
936 :
937 : ! Input Variables
938 : integer, intent(in) :: &
939 : nz ! Number of vertical levels
940 :
941 : real( kind = core_rknd ), dimension(nz), intent(in) :: &
942 : xm, & ! Mean (overall) of x, <x> [(x units)]
943 : xp2, & ! Variance (overall) of x, <x'^2> [(x units)^2]
944 : mixt_frac, & ! Mixture fraction [-]
945 : small_m, & ! Luhar's small m [-]
946 : wpxp ! Covariance of w and x, <w'x'> [m/s (x units)]
947 :
948 : real( kind = core_rknd ), intent(in) :: &
949 : x_tol_sqd ! Tolerance value of x squared [(x units)^2]
950 :
951 : ! Output Variables
952 : real( kind = core_rknd ), dimension(nz), intent(out) :: &
953 : sigma_sqd_x_1, & ! Normalized width parameter of x (1st PDF component) [-]
954 : sigma_sqd_x_2, & ! Normalized width parameter of x (1st PDF component) [-]
955 : varnce_x_1, & ! Variance of x (1st PDF component) [(x units)^2]
956 : varnce_x_2, & ! Variance of x (2nd PDF component) [(x units)^2]
957 : x_1_n, & ! Normalized mean of x (1st PDF component) [-]
958 : x_2_n, & ! Normalized mean of x (2nd PDF component) [-]
959 : x_1, & ! Mean of x (1st PDF component) [(x units)]
960 : x_2 ! Mean of x (2nd PDF component) [(x units)]
961 :
962 : ! Local Variables
963 : real( kind = core_rknd), dimension(nz) :: &
964 0 : sqrt_xp2, & ! Square root of the variance of x [(x units)]
965 0 : sgn_wpxp ! Sgn( <w'x'> ); 1 when <w'x'> >= 0 or -1 when <w'x'> < 0 [-]
966 :
967 :
968 0 : where ( xp2 > x_tol_sqd )
969 :
970 : ! Width (standard deviation) parameters are non-zero
971 :
972 : ! Calculate sgn( <w'x'> ).
973 : where ( wpxp >= zero )
974 : sgn_wpxp = one
975 : elsewhere ! <w'x'> < 0
976 : sgn_wpxp = -one
977 : endwhere ! <w'x'> >= 0
978 :
979 : ! Calculate the square root of the overall variance of x.
980 : sqrt_xp2 = sqrt( xp2 )
981 :
982 : ! Normalized width parameter of x in the 1st PDF component.
983 : sigma_sqd_x_1 &
984 : = ( one - mixt_frac ) / ( mixt_frac * ( one + small_m**2 ) )
985 :
986 : ! The variance of x in the 1st PDF component.
987 : varnce_x_1 = sigma_sqd_x_1 * xp2
988 :
989 : ! Normalized width parameter of x in the 2nd PDF component.
990 : sigma_sqd_x_2 &
991 : = mixt_frac / ( ( one - mixt_frac ) * ( one + small_m**2 ) )
992 :
993 : ! The variance of x in the 2nd PDF component.
994 : varnce_x_2 = sigma_sqd_x_2 * xp2
995 :
996 : ! Normalized mean of x in the 1st PDF component.
997 : x_1_n = sgn_wpxp * small_m * sqrt( sigma_sqd_x_1 )
998 :
999 : ! Normalized mean of x in the 2nd PDF component.
1000 : x_2_n = -sgn_wpxp * small_m * sqrt( sigma_sqd_x_2 )
1001 :
1002 : ! The mean of x in the 1st PDF component.
1003 : x_1 = xm + sqrt_xp2 * x_1_n
1004 :
1005 : ! The mean of x in the 2nd PDF component.
1006 : x_2 = xm + sqrt_xp2 * x_2_n
1007 :
1008 : elsewhere
1009 :
1010 : ! Variable x doesn't vary.
1011 : sigma_sqd_x_1 = ( one / ( one + small_m**2 ) )
1012 : sigma_sqd_x_2 = ( one / ( one + small_m**2 ) )
1013 : varnce_x_1 = zero
1014 : varnce_x_2 = zero
1015 : x_1_n = sgn_wpxp * small_m * sqrt( sigma_sqd_x_1 )
1016 : x_2_n = -sgn_wpxp * small_m * sqrt( sigma_sqd_x_2 )
1017 : x_1 = xm
1018 : x_2 = xm
1019 :
1020 : endwhere ! Widths non-zero
1021 :
1022 :
1023 0 : return
1024 :
1025 : end subroutine close_Luhar_pdf
1026 :
1027 : !=============================================================================
1028 71480448 : subroutine ADG1_ADG2_responder_params( nz, ngrdcol, & ! In
1029 71480448 : xm, xp2, wp2, sqrt_wp2, & ! In
1030 71480448 : wpxp, w_1_n, w_2_n, mixt_frac, & ! In
1031 71480448 : sigma_sqd_w, beta, x_tol, & ! In
1032 71480448 : x_1, x_2, varnce_x_1, & ! Out
1033 71480448 : varnce_x_2, alpha_x ) ! Out
1034 :
1035 : ! Description:
1036 : ! Calculates the PDF component means and PDF component variances for thl,
1037 : ! rt, or sclr when either the ADG1 PDF or ADG2 PDF are used.
1038 : !
1039 : ! The normalized variance for thl, rt, and sclr for "plume" 1 is:
1040 : !
1041 : ! { 1 - [1/(1-sigma_sqd_w)]*[ <w'x'>^2 / (<w'^2> * <x'^2>) ] / mixt_frac }
1042 : ! * { (1/3)*beta + mixt_frac*( 1 - (2/3)*beta ) };
1043 : !
1044 : ! where "x" stands for thl, rt, or sclr; "mixt_frac" is the weight of
1045 : ! Gaussian "plume" 1, and 0 <= beta <= 3.
1046 : !
1047 : ! The factor { (1/3)*beta + mixt_frac*( 1 - (2/3)*beta ) } does not depend
1048 : ! on which varable "x" stands for. The factor is multiplied by 2 and
1049 : ! defined as width_factor_1.
1050 : !
1051 : ! The factor
1052 : ! { 1 - [1/(1-sigma_sqd_w)]*[ (w'x')^2 / (w'^2 * x'^2) ] / mixt_frac }
1053 : ! depends on which variable "x" stands for. It is multiplied by 1/2 and
1054 : ! defined as alpha_x, where "x" stands for thl, rt, or sclr.
1055 :
1056 : ! Vince Larson added a dimensionless factor so that the
1057 : ! width of plumes in theta_l, rt can vary.
1058 : ! beta is a constant defined in CLUBB's tunable parameters.
1059 : ! Set 0<beta<3.
1060 : ! beta=1.5_core_rknd recovers Chris Golaz' simplified formula.
1061 : ! 3 Nov 2003
1062 :
1063 : ! References:
1064 :
1065 : use constants_clubb, only: &
1066 : two, & ! Constant(s)
1067 : one, &
1068 : two_thirds, &
1069 : one_half, &
1070 : zero_threshold, &
1071 : w_tol_sqd
1072 :
1073 : use clubb_precision, only: &
1074 : core_rknd ! Variable(s)
1075 :
1076 : implicit none
1077 :
1078 : integer, intent(in) :: &
1079 : ngrdcol, & ! Number of grid columns
1080 : nz ! Number of vertical level
1081 :
1082 : ! Input Variables
1083 : real ( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
1084 : xm, & ! Mean of x (overall) [units vary]
1085 : xp2, & ! Variance of x (overall) [(units vary)^2]
1086 : wp2, & ! Variance of w (overall) [m^2/s^2]
1087 : sqrt_wp2, & ! Square root of the variance of w [m/s]
1088 : wpxp, & ! Covariance of w and x [m/s (units vary)]
1089 : w_1_n, & ! Normalized mean of w (1st PDF comp.) [-]
1090 : w_2_n, & ! Normalized mean of w (2nd PDF comp.) [-]
1091 : mixt_frac, & ! Mixture fraction [-]
1092 : sigma_sqd_w ! Width of individual w plumes [-]
1093 :
1094 : real ( kind = core_rknd ), dimension(ngrdcol), intent(in) :: &
1095 : beta ! CLUBB tunable parameter beta [-]
1096 :
1097 : real ( kind = core_rknd ), intent(in) :: &
1098 : x_tol ! Tolerance value for x [units vary]
1099 :
1100 : ! Output Variables
1101 : real ( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
1102 : x_1, & ! Mean of x (1st PDF component) [units vary]
1103 : x_2, & ! Mean of x (2nd PDF component) [units vary]
1104 : varnce_x_1, & ! Variance of x (1st PDF component) [(units vary)^2]
1105 : varnce_x_2, & ! Variance of x (2nd PDF component) [(units vary)^2]
1106 : alpha_x ! Factor relating to normalized variance for x [-]
1107 :
1108 : ! Local Variables
1109 : ! variables for a generalization of Chris Golaz' closure
1110 : ! varies width of plumes in theta_l, rt
1111 : real ( kind = core_rknd ) :: &
1112 : width_factor_1 ! Width factor relating to PDF component 1 [-]
1113 :
1114 : integer :: &
1115 : k, i ! Vertical loop index
1116 :
1117 : !----- Begin Code -----
1118 : !$acc parallel loop gang vector collapse(2) default(present)
1119 6147318528 : do k = 1, nz
1120 >10152*10^7 : do i = 1, ngrdcol
1121 :
1122 >10145*10^7 : if ( wp2(i,k) <= w_tol_sqd .or. xp2(i,k) <= x_tol**2 ) then
1123 :
1124 : ! Vertical velocity doesn't vary. Variable x is treated as a single
1125 : ! Gaussian in this situation.
1126 72247479667 : x_1(i,k) = xm(i,k)
1127 72247479667 : x_2(i,k) = xm(i,k)
1128 72247479667 : varnce_x_1(i,k) = xp2(i,k)
1129 72247479667 : varnce_x_2(i,k) = xp2(i,k)
1130 72247479667 : alpha_x(i,k) = one_half
1131 :
1132 : else
1133 :
1134 23129048333 : x_1(i,k) = xm(i,k) - wpxp(i,k) / ( sqrt_wp2(i,k) * w_2_n(i,k) )
1135 23129048333 : x_2(i,k) = xm(i,k) - wpxp(i,k) / ( sqrt_wp2(i,k) * w_1_n(i,k) )
1136 :
1137 : alpha_x(i,k) = one_half &
1138 : * ( one - wpxp(i,k) * wpxp(i,k) &
1139 23129048333 : / ( ( one - sigma_sqd_w(i,k) ) * wp2(i,k) * xp2(i,k) ) )
1140 :
1141 23129048333 : alpha_x(i,k) = max( min( alpha_x(i,k), one ), zero_threshold )
1142 :
1143 : width_factor_1 = two_thirds * beta(i) &
1144 23129048333 : + two * mixt_frac(i,k) * ( one - two_thirds * beta(i) )
1145 :
1146 : ! Vince Larson multiplied original expressions by width_factor_1,2
1147 : ! to generalize scalar skewnesses. 05 Nov 03
1148 23129048333 : varnce_x_1(i,k) = width_factor_1 * xp2(i,k) * alpha_x(i,k) / mixt_frac(i,k)
1149 : varnce_x_2(i,k) = ( two - width_factor_1 ) * xp2(i,k) &
1150 23129048333 : * alpha_x(i,k) / ( one - mixt_frac(i,k) )
1151 :
1152 : end if
1153 :
1154 : end do
1155 : end do
1156 : !$acc end parallel loop
1157 :
1158 71480448 : return
1159 :
1160 : end subroutine ADG1_ADG2_responder_params
1161 :
1162 : !=============================================================================
1163 0 : subroutine backsolve_Luhar_params( Sk_max, Skx, & ! In
1164 : big_m_max, mixt_frac, & ! In
1165 : big_m_x, small_m_x ) ! Out
1166 :
1167 : ! Description:
1168 : ! This subroutine calculates Luhar's big_m and small_m for the variate 'x'
1169 : ! consistent with the mixture fraction of the variate with the largest
1170 : ! skewness.
1171 : !
1172 : ! The relationship between skewness of x (Skx), mixture fraction (a), and
1173 : ! Luhar's small m (m) is given by:
1174 : !
1175 : ! Skx^2 = ( m^2 * ( m^2 + 3 )^2 / ( m^2 + 1 )^3 )
1176 : ! * ( 1 - 2*a )^2 / ( a * ( 1 - a ) ).
1177 : !
1178 : ! Moving the factor involving mixture fraction to the right-hand side:
1179 : !
1180 : ! ( ( a * ( 1 - a ) ) / ( 1 - 2*a )^2 ) * Skx^2
1181 : ! = m^2 * ( m^2 + 3 )^2 / ( m^2 + 1 )^3.
1182 : !
1183 : ! This can be rewritten as:
1184 : !
1185 : ! ( ( a * ( 1 - a ) ) / ( 1 - 2*a )^2 ) * Skx^2
1186 : ! = ( m^6 + 6*m^4 + 9*m^2 ) / ( m^6 + 3*m^4 + 3*m^2 + 1 ).
1187 : !
1188 : ! Setting alpha = ( ( a * ( 1 - a ) ) / ( 1 - 2*a )^2 ) * Skx^2, the
1189 : ! equation can be rewritten as:
1190 : !
1191 : ! ( m^6 + 3*m^4 + 3*m^2 + 1 ) * alpha = m^6 + 6*m^4 + 9*m^2.
1192 : !
1193 : ! This can be rearranged and rewritten as:
1194 : !
1195 : ! ( alpha - 1 ) * m^6 + ( 3 * alpha - 6 ) * m^4
1196 : ! + ( 3 * alpha - 9 ) * m^2 + alpha = 0.
1197 : !
1198 : ! This can be rewritten again as:
1199 : !
1200 : ! ( alpha - 1 ) * (m^2)^3 + ( 3 * alpha - 6 ) * (m^2)^2
1201 : ! + ( 3 * alpha - 9 ) * (m^2) + alpha = 0.
1202 : !
1203 : ! The goal is to solve for m^2, and then take the square root of m^2 to
1204 : ! solve for m. This can be accomplished by using the cubic formula (with
1205 : ! the l_use_cubic_backsolve option), or else by a quadratic approximation.
1206 :
1207 : ! References:
1208 : !-----------------------------------------------------------------------
1209 :
1210 : use constants_clubb, only: &
1211 : nine, & ! Constant(s)
1212 : six, &
1213 : five, &
1214 : four, &
1215 : three, &
1216 : two, &
1217 : one, &
1218 : one_half, &
1219 : zero, &
1220 : eps
1221 :
1222 : use clubb_precision, only: &
1223 : core_rknd ! Precision
1224 :
1225 : implicit none
1226 :
1227 : ! Input Variables
1228 : real( kind = core_rknd ), intent(in) :: &
1229 : Sk_max, & ! Maximum skewness [-]
1230 : Skx, & ! Skewness of the variate solving small_m and big_m for [-]
1231 : big_m_max, & ! Luhar's big_m of the variate with maximum skewness [-]
1232 : mixt_frac ! Mixture fraction [-]
1233 :
1234 : ! Output Variables
1235 : real( kind = core_rknd ), intent(out) :: &
1236 : big_m_x, & ! Luhar's big_m for the variate being solved for [-]
1237 : small_m_x ! Luhar's small_m for the variate being solved for [-]
1238 :
1239 : ! Local Variables
1240 : real( kind = core_rknd ) :: &
1241 : alpha, & ! 1 / big_m_x
1242 : a, & ! For readability, quadratic equation
1243 : b, &
1244 : c, &
1245 : alpha_upr, &
1246 : alpha_low, &
1247 : discrim
1248 :
1249 : real( kind = core_rknd ), dimension(1) :: &
1250 : a_coef_in, & ! Coefficient a (of x^3) in a*x^3 + b*x^2 + c^x + d = 0 [-]
1251 : b_coef_in, & ! Coefficient b (of x^2) in a*x^3 + b*x^2 + c^x + d = 0 [-]
1252 : c_coef_in, & ! Coefficient c (of x) in a*x^3 + b*x^2 + c^x + d = 0 [-]
1253 : d_coef_in ! Coefficient d in a*x^3 + b*x^2 + c^x + d = 0 [-]
1254 :
1255 : ! Flag to backsolve for m^2 using cubic formula
1256 : logical, parameter :: &
1257 : l_use_cubic_backsolve = .true.
1258 :
1259 :
1260 : if ( l_use_cubic_backsolve ) then
1261 :
1262 0 : if ( abs( mixt_frac - one_half ) < 0.001_core_rknd ) then
1263 :
1264 : ! When mixture fraction = 0.5 (based on the variable with the largest
1265 : ! magnitude of skewness), all variables must have a skewness of 0.
1266 : ! Set m to the minimum threshold of 0.05.
1267 0 : small_m_x = 0.05_core_rknd
1268 :
1269 : ! Calculate the corresponding value of big_m_x.
1270 : big_m_x = ( one + small_m_x**2 )**3 &
1271 0 : / ( ( three + small_m_x**2 )**2 * small_m_x**2 )
1272 :
1273 0 : elseif ( abs(Skx) < eps ) then
1274 :
1275 : ! Mixture fraction /= 0.5 because the variable with the largest
1276 : ! magnitude of skewness has a skewness /= 0. However, variable x has
1277 : ! a skewness of 0. In order to reproduce the correct skewness for
1278 : ! variable x, set m to 0 (regardless of minimum thresholds used in
1279 : ! other parts of the code).
1280 0 : small_m_x = zero
1281 :
1282 : ! The value of big_m_x should be inf. Set it to huge. This is not
1283 : ! used in any calculation, anyway.
1284 0 : big_m_x = huge( big_m_x )
1285 :
1286 : else ! mixt_frac /= 0.5 and Skx /= 0
1287 :
1288 : ! Backsolve for m, given mixt_frac and Skx.
1289 :
1290 : ! alpha = 1/M is given by:
1291 : ! [ mixt_frac * ( 1 - mixt_frac ) / ( 1 - 2 * mixt_frac )^2 ] * Skx^2.
1292 : alpha = ( mixt_frac * ( one - mixt_frac ) &
1293 0 : / ( one - two * mixt_frac )**2 ) * Skx**2
1294 :
1295 : ! Calculate big_m_x.
1296 0 : big_m_x = one / alpha
1297 :
1298 : ! Solve the cubic equation for m^2:
1299 : ! ( alpha - 1 ) * (m^2)^3 + ( 3 * alpha - 6 ) * (m^2)^2
1300 : ! + ( 3 * alpha - 9 ) * (m^2) + alpha = 0.
1301 : ! The largest root is preferred.
1302 0 : a_coef_in(1) = alpha - one
1303 0 : b_coef_in(1) = three * alpha - six
1304 0 : c_coef_in(1) = three * alpha - nine
1305 0 : d_coef_in(1) = alpha
1306 : small_m_x &
1307 : = sqrt( max( max_cubic_root( a_coef_in, b_coef_in, &
1308 : c_coef_in, d_coef_in ), &
1309 0 : 0.05_core_rknd**2 ) )
1310 :
1311 : endif
1312 :
1313 : else ! original formualation
1314 :
1315 : alpha = ( Skx**2 / (max(Sk_max**2, eps) * big_m_max) ) ! 1 / big_m_x
1316 :
1317 : ! This limit keeps the discriminant >= 0
1318 : alpha_upr = two * sqrt( 13.0_core_rknd ) - five
1319 :
1320 : alpha_low = eps
1321 :
1322 : ! For this approximation, alpha must be less than 2*sqrt(13) - 5 to get a
1323 : ! real ans.
1324 : alpha = min( alpha, alpha_upr )
1325 :
1326 : ! For testing, eliminate possibility of divide by zero
1327 : alpha = max( alpha, alpha_low )
1328 :
1329 : ! Use a piece-wise approximation
1330 : if ( alpha < one ) then
1331 :
1332 : a = max( three * alpha - six, eps) ! Prevent divide by zero
1333 : b = three * alpha - nine
1334 : c = alpha
1335 :
1336 : discrim = b**2 - four * a * c
1337 : small_m_x = sqrt( ( - b - sqrt( discrim ) ) / ( two * a ) )
1338 :
1339 : else ! alpha >= 1
1340 :
1341 : ! For this approximation, alpha must be less than 2*sqrt(13) - 5
1342 : ! to get a real ans.
1343 : alpha = min( alpha, two )
1344 :
1345 : a = max( six * alpha - nine, eps) ! Prevent divide by zero
1346 : b = -six
1347 : c = two * alpha - one
1348 :
1349 : discrim = b**2 - four * a * c
1350 : small_m_x = sqrt( ( - b - sqrt( discrim ) ) / ( two * a ) )
1351 :
1352 : endif ! alpha < 1
1353 :
1354 : ! Clip consistently with subroutine calc_Luhar_params
1355 : small_m_x = max( 5e-2_core_rknd, small_m_x)
1356 :
1357 : big_m_x = one / alpha
1358 :
1359 : endif ! l_use_cubic_backsolve
1360 :
1361 :
1362 0 : return
1363 :
1364 : end subroutine backsolve_Luhar_params
1365 :
1366 : !=============================================================================
1367 0 : function max_cubic_root( a_coef, b_coef, c_coef, d_coef ) &
1368 : result( max_root )
1369 :
1370 : ! Description:
1371 : ! Calculates the largest root that results from solving a cubic equation of
1372 : ! the form a*x^3 + b*x^2 + c*x + d = 0.
1373 : !
1374 : ! This is done to backsolve for m^2 for the 3-D Luhar closure, given the
1375 : ! values of mixt_frac and Skx.
1376 :
1377 : ! References:
1378 : !-----------------------------------------------------------------------
1379 :
1380 : use constants_clubb, only: &
1381 : eps ! Constant(s)
1382 :
1383 : use calc_roots, only: &
1384 : cubic_solve, & ! Procedure(s)
1385 : quadratic_solve
1386 :
1387 : use clubb_precision, only: &
1388 : core_rknd ! Variable(s)
1389 :
1390 : implicit none
1391 :
1392 : ! Input Variables
1393 : real( kind = core_rknd ), dimension(1), intent(in) :: &
1394 : a_coef, & ! Coefficient a (of x^3) in a*x^3 + b*x^2 + c^x + d = 0 [-]
1395 : b_coef, & ! Coefficient b (of x^2) in a*x^3 + b*x^2 + c^x + d = 0 [-]
1396 : c_coef, & ! Coefficient c (of x) in a*x^3 + b*x^2 + c^x + d = 0 [-]
1397 : d_coef ! Coefficient d in a*x^3 + b*x^2 + c^x + d = 0 [-]
1398 :
1399 : ! Return Variable
1400 : real( kind = core_rknd ) :: &
1401 : max_root ! Maximum root that solves the cubic equation [-]
1402 :
1403 : ! Local Variables
1404 : complex( kind = core_rknd ), dimension(1,3) :: &
1405 : cubic_roots ! Roots of x that satisfy a*x^3 + b*x^2 + c*x + d = 0 [-]
1406 :
1407 : complex( kind = core_rknd ), dimension(1,2) :: &
1408 : quadratic_roots ! Roots of x that satisfy b*x^2 + c*x + d = 0 [-]
1409 :
1410 : real( kind = core_rknd ) :: &
1411 : a_coef_thresh, & ! Minimum threshold of |a| to use cubic solver [-]
1412 : b_coef_thresh ! Minimum threshold of |b| to use quadratic solver [-]
1413 :
1414 :
1415 : ! Calculate a minimum threshold for |a| to call this a cubic equation.
1416 : a_coef_thresh = 0.001_core_rknd &
1417 0 : * max( abs(b_coef(1)), abs(c_coef(1)), abs(d_coef(1)) )
1418 :
1419 : ! Calculate a minimum threshold for |b| to call this a quadratic equation.
1420 : ! This only matters when |a| <= a_coef_thresh.
1421 0 : b_coef_thresh = 0.001_core_rknd * max( abs(c_coef(1)), abs(d_coef(1)) )
1422 :
1423 0 : if ( abs( a_coef(1) ) > a_coef_thresh ) then
1424 :
1425 : ! The equation is a cubic equation.
1426 0 : cubic_roots = cubic_solve( 1, a_coef, b_coef, c_coef, d_coef )
1427 :
1428 0 : if ( abs(aimag( cubic_roots(1,2) )) < eps .and. &
1429 : abs(aimag( cubic_roots(1,3) )) < eps ) then
1430 :
1431 : ! Find the maximum root of the three roots.
1432 : max_root = max( real( cubic_roots(1,1), kind = core_rknd ), &
1433 : real( cubic_roots(1,2), kind = core_rknd ), &
1434 0 : real( cubic_roots(1,3), kind = core_rknd ) )
1435 :
1436 : else ! cubic_roots(2) and cubic_roots(3) are complex.
1437 :
1438 0 : max_root = real( cubic_roots(1,1), kind = core_rknd )
1439 :
1440 : endif
1441 :
1442 0 : elseif ( abs( b_coef(1) ) > b_coef_thresh ) then
1443 :
1444 : ! The equation is a quadratic equation, since a = 0, but b /= 0.
1445 : ! This should very rarely occur for 3-D Luhar. When it does, the result
1446 : ! will always be two real-valued roots.
1447 0 : quadratic_roots = quadratic_solve( 1, b_coef, c_coef, d_coef )
1448 :
1449 : ! Find the maximum root of the two roots.
1450 : max_root = max( real( quadratic_roots(1,1), kind = core_rknd ), &
1451 0 : real( quadratic_roots(1,2), kind = core_rknd ) )
1452 :
1453 : else ! |a| = 0 and |b| = 0
1454 :
1455 : ! The equation is a linear equation.
1456 : ! This won't happen for 3-D Luhar.
1457 0 : max_root = - d_coef(1) / c_coef(1)
1458 :
1459 : endif ! |a| > 0
1460 :
1461 :
1462 : return
1463 :
1464 : end function max_cubic_root
1465 :
1466 : !=============================================================================
1467 :
1468 : end module adg1_adg2_3d_luhar_pdf
|