Line data Source code
1 : !---------------------------------------------------------------------------
2 : ! $Id$
3 : !===============================================================================
4 : module pdf_closure_module
5 :
6 : ! Options for the two component normal (double Gaussian) PDF type to use for
7 : ! the w, rt, and theta-l (or w, chi, and eta) portion of CLUBB's multivariate,
8 : ! two-component PDF.
9 : use model_flags, only: &
10 : iiPDF_ADG1, & ! ADG1 PDF
11 : iiPDF_ADG2, & ! ADG2 PDF
12 : iiPDF_3D_Luhar, & ! 3D Luhar PDF
13 : iiPDF_new, & ! new PDF
14 : iiPDF_TSDADG, & ! new TSDADG PDF
15 : iiPDF_LY93, & ! Lewellen and Yoh (1993)
16 : iiPDF_new_hybrid ! new hybrid PDF
17 :
18 : implicit none
19 :
20 : public :: pdf_closure, &
21 : calc_wp4_pdf, &
22 : calc_wp2xp_pdf, &
23 : calc_wpxp2_pdf, &
24 : calc_wpxpyp_pdf, &
25 : calc_w_up_in_cloud
26 :
27 : private ! Set Default Scope
28 :
29 : contains
30 : !------------------------------------------------------------------------
31 :
32 : !#######################################################################
33 : !#######################################################################
34 : ! If you change the argument list of pdf_closure you also have to
35 : ! change the calls to this function in the host models CAM, WRF, SAM
36 : ! and GFDL.
37 : !#######################################################################
38 : !#######################################################################
39 17870112 : subroutine pdf_closure( nz, ngrdcol, sclr_dim, sclr_tol, &
40 17870112 : hydromet_dim, p_in_Pa, exner, thv_ds, &
41 17870112 : wm, wp2, wp3, &
42 17870112 : Skw, Skthl_in, Skrt_in, Sku_in, Skv_in, &
43 17870112 : rtm, rtp2, wprtp, &
44 17870112 : thlm, thlp2, wpthlp, &
45 17870112 : um, up2, upwp, &
46 17870112 : vm, vp2, vpwp, &
47 17870112 : rtpthlp, &
48 17870112 : sclrm, wpsclrp, sclrp2, &
49 17870112 : sclrprtp, sclrpthlp, Sksclr_in, &
50 17870112 : gamma_Skw_fnc, &
51 : #ifdef GFDL
52 : RH_crit, do_liquid_only_in_clubb, & ! h1g, 2010-06-15
53 : #endif
54 17870112 : wphydrometp, wp2hmp, &
55 17870112 : rtphmp, thlphmp, &
56 17870112 : clubb_params, &
57 : saturation_formula, &
58 : stats_metadata, &
59 : iiPDF_type, &
60 17870112 : l_mix_rat_hm, &
61 17870112 : sigma_sqd_w, &
62 : pdf_params, pdf_implicit_coefs_terms, &
63 17870112 : wpup2, wpvp2, &
64 17870112 : wp2up2, wp2vp2, wp4, &
65 17870112 : wprtp2, wp2rtp, &
66 17870112 : wpthlp2, wp2thlp, wprtpthlp, &
67 17870112 : cloud_frac, ice_supersat_frac, &
68 17870112 : rcm, wpthvp, wp2thvp, rtpthvp, &
69 17870112 : thlpthvp, wprcp, wp2rcp, rtprcp, &
70 17870112 : thlprcp, rcp2, &
71 17870112 : uprcp, vprcp, &
72 17870112 : w_up_in_cloud, w_down_in_cloud, &
73 17870112 : cloudy_updraft_frac, cloudy_downdraft_frac, &
74 17870112 : F_w, F_rt, F_thl, &
75 17870112 : min_F_w, max_F_w, &
76 17870112 : min_F_rt, max_F_rt, &
77 17870112 : min_F_thl, max_F_thl, &
78 17870112 : wpsclrprtp, wpsclrp2, sclrpthvp, &
79 17870112 : wpsclrpthlp, sclrprcp, wp2sclrp, &
80 17870112 : rc_coef )
81 :
82 :
83 : ! Description:
84 : ! Subroutine that computes pdf parameters analytically.
85 : !
86 : ! Based of the original formulation, but with some tweaks
87 : ! to remove some of the less realistic assumptions and
88 : ! improve transport terms.
89 :
90 : ! Corrected version that should remove inconsistency
91 :
92 : ! References:
93 : ! The shape of CLUBB's PDF is given by the expression in
94 : ! https://arxiv.org/pdf/1711.03675v1.pdf#nameddest=url:clubb_pdf
95 :
96 : ! Eqn. 29, 30, 31, 32 & 33 on p. 3547 of
97 : ! ``A PDF-Based Model for Boundary Layer Clouds. Part I:
98 : ! Method and Model Description'' Golaz, et al. (2002)
99 : ! JAS, Vol. 59, pp. 3540--3551.
100 : !----------------------------------------------------------------------
101 :
102 : use grid_class, only: &
103 : grid ! Type
104 :
105 : use constants_clubb, only: & ! Constants
106 : three, & ! 3
107 : one, & ! 1
108 : one_half, & ! 1/2
109 : zero, & ! 0
110 : Cp, & ! Dry air specific heat at constant p [J/kg/K]
111 : Lv, & ! Latent heat of vaporization [J/kg]
112 : ep1, & ! (1.0-ep)/ep; ep1 = 0.61 [-]
113 : ep2, & ! 1.0/ep; ep2 = 1.61 [-]
114 : rt_tol, & ! Tolerance for r_t [kg/kg]
115 : thl_tol, & ! Tolerance for th_l [K]
116 : fstderr, &
117 : zero_threshold, &
118 : eps, &
119 : w_tol
120 :
121 : use parameters_model, only: &
122 : mixt_frac_max_mag ! Variable(s)
123 :
124 : use parameter_indices, only: &
125 : nparams, & ! Variable(s)
126 : ibeta, &
127 : iSkw_denom_coef, &
128 : ipdf_component_stdev_factor_w, &
129 : islope_coef_spread_DG_means_w
130 :
131 : use pdf_parameter_module, only: &
132 : pdf_parameter, & ! Variable Type
133 : implicit_coefs_terms
134 :
135 : use new_pdf_main, only: &
136 : new_pdf_driver ! Procedure(s)
137 :
138 : use new_hybrid_pdf_main, only: &
139 : new_hybrid_pdf_driver ! Procedure(s)
140 :
141 : use adg1_adg2_3d_luhar_pdf, only: &
142 : ADG1_pdf_driver, & ! Procedure(s)
143 : ADG2_pdf_driver, &
144 : Luhar_3D_pdf_driver
145 :
146 : use new_tsdadg_pdf, only: &
147 : tsdadg_pdf_driver ! Procedure(s)
148 :
149 : use LY93_pdf, only: &
150 : LY93_driver ! Procedure(s)
151 :
152 : use pdf_utilities, only: &
153 : calc_comp_corrs_binormal, & ! Procedure(s)
154 : calc_corr_chi_x, &
155 : calc_corr_eta_x
156 :
157 : use model_flags, only: &
158 : l_explicit_turbulent_adv_xpyp ! Variable(s)
159 :
160 : use numerical_check, only: &
161 : pdf_closure_check ! Procedure(s)
162 :
163 : use saturation, only: &
164 : sat_mixrat_liq, & ! Procedure(s)
165 : sat_mixrat_ice
166 :
167 : use clubb_precision, only: &
168 : core_rknd ! Variable(s)
169 :
170 : use error_code, only: &
171 : clubb_at_least_debug_level, & ! Procedure
172 : err_code, & ! Error Indicator
173 : clubb_fatal_error ! Constant
174 :
175 : use stats_variables, only: &
176 : stats_metadata_type
177 :
178 : implicit none
179 :
180 : !----------------------------- Input Variables -----------------------------
181 : integer, intent(in) :: &
182 : nz, & ! Number of vertical levels
183 : ngrdcol, & ! Number of grid columns
184 : hydromet_dim, & ! Number of hydrometeor species
185 : sclr_dim ! Number of passive scalars
186 :
187 : real( kind = core_rknd ), intent(in), dimension(sclr_dim) :: &
188 : sclr_tol ! Threshold(s) on the passive scalars [units vary]
189 :
190 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
191 : p_in_Pa, & ! Pressure [Pa]
192 : exner, & ! Exner function [-]
193 : thv_ds, & ! Dry, base-state theta_v (ref. th_l here) [K]
194 : wm, & ! mean w-wind component (vertical velocity) [m/s]
195 : wp2, & ! w'^2 [m^2/s^2]
196 : wp3, & ! w'^3 [m^3/s^3]
197 : Skw, & ! Skewness of w [-]
198 : Skthl_in, & ! Skewness of thl [-]
199 : Skrt_in, & ! Skewness of rt [-]
200 : Sku_in, & ! Skewness of u [-]
201 : Skv_in, & ! Skewness of v [-]
202 : rtm, & ! Mean total water mixing ratio [kg/kg]
203 : rtp2, & ! r_t'^2 [(kg/kg)^2]
204 : wprtp, & ! w'r_t' [(kg/kg)(m/s)]
205 : thlm, & ! Mean liquid water potential temperature [K]
206 : thlp2, & ! th_l'^2 [K^2]
207 : wpthlp, & ! w'th_l' [K(m/s)]
208 : rtpthlp ! r_t'th_l' [K(kg/kg)]
209 :
210 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
211 : um, & ! Grid-mean eastward wind [m/s]
212 : up2, & ! u'^2 [(m/s)^2]
213 : upwp, & ! u'w' [(m/s)^2]
214 : vm, & ! Grid-mean northward wind [m/s]
215 : vp2, & ! v'^2 [(m/s)^2]
216 : vpwp ! v'w' [(m/s)^2]
217 :
218 : real( kind = core_rknd ), dimension(ngrdcol,nz, sclr_dim), intent(in) :: &
219 : sclrm, & ! Mean passive scalar [units vary]
220 : wpsclrp, & ! w' sclr' [units vary]
221 : sclrp2, & ! sclr'^2 [units vary]
222 : sclrprtp, & ! sclr' r_t' [units vary]
223 : sclrpthlp, & ! sclr' th_l' [units vary]
224 : Sksclr_in ! Skewness of sclr [-]
225 :
226 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
227 : gamma_Skw_fnc ! Gamma as a function of skewness [-]
228 :
229 : #ifdef GFDL
230 : ! critial relative humidity for nucleation
231 : real( kind = core_rknd ), dimension(ngrdcol, nz, min(1,sclr_dim), 2 ), intent(in) :: & ! h1g, 2010-06-15
232 : RH_crit ! critical relative humidity for droplet and ice nucleation
233 : ! ---> h1g, 2012-06-14
234 : logical, intent(in) :: do_liquid_only_in_clubb
235 : ! <--- h1g, 2012-06-14
236 : #endif
237 :
238 : real( kind = core_rknd ), dimension(ngrdcol,nz,hydromet_dim), intent(in) :: &
239 : wphydrometp, & ! Covariance of w and a hydrometeor [(m/s) <hm units>]
240 : wp2hmp, & ! Third-order moment: < w'^2 hm' > [(m/s)^2 <hm units>]
241 : rtphmp, & ! Covariance of rt and a hydrometeor [(kg/kg) <hm units>]
242 : thlphmp ! Covariance of thl and a hydrometeor [K <hm units>]
243 :
244 : real( kind = core_rknd ), dimension(ngrdcol,nparams), intent(in) :: &
245 : clubb_params ! Array of CLUBB's tunable parameters [units vary]
246 :
247 : integer, intent(in) :: &
248 : iiPDF_type ! Selected option for the two-component normal (double
249 : ! Gaussian) PDF type to use for the w, rt, and theta-l (or
250 : ! w, chi, and eta) portion of CLUBB's multivariate,
251 : ! two-component PDF.
252 :
253 : logical, dimension(hydromet_dim), intent(in) :: &
254 : l_mix_rat_hm ! if true, then the quantity is a hydrometeor mixing ratio
255 :
256 : integer, intent(in) :: &
257 : saturation_formula ! Integer that stores the saturation formula to be used
258 :
259 : type (stats_metadata_type), intent(in) :: &
260 : stats_metadata
261 :
262 : !----------------------------- InOut Variables -----------------------------
263 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(inout) :: &
264 : ! If iiPDF_type == iiPDF_ADG2, this gets overwritten. Therefore,
265 : ! intent(inout). Otherwise it should be intent(in)
266 : sigma_sqd_w ! Width of individual w plumes [-]
267 :
268 : type(pdf_parameter), intent(inout) :: &
269 : pdf_params ! pdf paramters [units vary]
270 :
271 : type(implicit_coefs_terms), intent(inout) :: &
272 : pdf_implicit_coefs_terms ! Implicit coefs / explicit terms [units vary]
273 :
274 : !----------------------------- Output Variables -----------------------------
275 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
276 : wpup2, & ! w'u'^2 [m^3/s^3]
277 : wpvp2, & ! w'v'^2 [m^3/s^3]
278 : wp2up2, & ! w'^2u'^2 [m^2/s^4]
279 : wp2vp2, & ! w'^2v'^2 [m^2/s^4]
280 : wp4, & ! w'^4 [m^4/s^4]
281 : wprtp2, & ! w' r_t' [(m kg)/(s kg)]
282 : wp2rtp, & ! w'^2 r_t' [(m^2 kg)/(s^2 kg)]
283 : wpthlp2, & ! w' th_l'^2 [(m K^2)/s]
284 : wp2thlp, & ! w'^2 th_l' [(m^2 K)/s^2]
285 : cloud_frac, & ! Cloud fraction [-]
286 : ice_supersat_frac, & ! Ice cloud fracion [-]
287 : rcm, & ! Mean liquid water [kg/kg]
288 : wpthvp, & ! Buoyancy flux [(K m)/s]
289 : wp2thvp, & ! w'^2 th_v' [(m^2 K)/s^2]
290 : rtpthvp, & ! r_t' th_v' [(kg K)/kg]
291 : thlpthvp, & ! th_l' th_v' [K^2]
292 : wprcp, & ! w' r_c' [(m kg)/(s kg)]
293 : wp2rcp, & ! w'^2 r_c' [(m^2 kg)/(s^2 kg)]
294 : rtprcp, & ! r_t' r_c' [(kg^2)/(kg^2)]
295 : thlprcp, & ! th_l' r_c' [(K kg)/kg]
296 : rcp2, & ! r_c'^2 [(kg^2)/(kg^2)]
297 : wprtpthlp, & ! w' r_t' th_l' [(m kg K)/(s kg)]
298 : w_up_in_cloud, & ! cloudy updraft vel [m/s]
299 : w_down_in_cloud, & ! cloudy downdraft vel [m/s]
300 : cloudy_updraft_frac, & ! cloudy updraft fraction [-]
301 : cloudy_downdraft_frac ! cloudy downdraft fraction [-]
302 :
303 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
304 : uprcp, & ! u' r_c' [(m kg)/(s kg)]
305 : vprcp ! v' r_c' [(m kg)/(s kg)]
306 :
307 : ! Parameters output only for recording statistics (new PDF).
308 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
309 : F_w, & ! Parameter for the spread of the PDF component means of w [-]
310 : F_rt, & ! Parameter for the spread of the PDF component means of rt [-]
311 : F_thl ! Parameter for the spread of the PDF component means of thl [-]
312 :
313 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
314 : min_F_w, & ! Minimum allowable value of parameter F_w [-]
315 : max_F_w, & ! Maximum allowable value of parameter F_w [-]
316 : min_F_rt, & ! Minimum allowable value of parameter F_rt [-]
317 : max_F_rt, & ! Maximum allowable value of parameter F_rt [-]
318 : min_F_thl, & ! Minimum allowable value of parameter F_thl [-]
319 : max_F_thl ! Maximum allowable value of parameter F_thl [-]
320 :
321 : ! Output (passive scalar variables)
322 : real( kind = core_rknd ), intent(out), dimension(ngrdcol,nz,sclr_dim) :: &
323 : sclrpthvp, &
324 : sclrprcp, &
325 : wpsclrp2, &
326 : wpsclrprtp, &
327 : wpsclrpthlp, &
328 : wp2sclrp
329 :
330 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
331 : rc_coef ! Coefficient on X'r_c' in X'th_v' equation [K/(kg/kg)]
332 :
333 : !----------------------------- Local Variables -----------------------------
334 :
335 : ! Variables that are stored in derived data type pdf_params.
336 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
337 35740224 : u_1, & ! Mean of eastward wind (1st PDF component) [m/s]
338 35740224 : u_2, & ! Mean of eastward wind (2nd PDF component) [m/s]
339 35740224 : varnce_u_1, & ! Variance of u (1st PDF component) [m^2/s^2]
340 35740224 : varnce_u_2, & ! Variance of u (2nd PDF component) [m^2/s^2]
341 35740224 : v_1, & ! Mean of northward wind (1st PDF component) [m/s]
342 35740224 : v_2, & ! Mean of northward wind (2nd PDF component) [m/s]
343 35740224 : varnce_v_1, & ! Variance of v (1st PDF component) [m^2/s^2]
344 35740224 : varnce_v_2, & ! Variance of v (2nd PDF component) [m^2/s^2]
345 35740224 : alpha_u, & ! Factor relating to normalized variance for u [-]
346 35740224 : alpha_v ! Factor relating to normalized variance for v [-]
347 :
348 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
349 35740224 : corr_u_w_1, & ! Correlation of u and w (1st PDF component) [-]
350 35740224 : corr_u_w_2, & ! Correlation of u and w (2nd PDF component) [-]
351 35740224 : corr_v_w_1, & ! Correlation of v and w (1st PDF component) [-]
352 35740224 : corr_v_w_2 ! Correlation of v and w (2nd PDF component) [-]
353 :
354 : ! Note: alpha coefficients = 0.5 * ( 1 - correlations^2 ).
355 : ! These are used to calculate the scalar widths
356 : ! varnce_thl_1, varnce_thl_2, varnce_rt_1, and varnce_rt_2 as in
357 : ! Eq. (34) of Larson and Golaz (2005)
358 :
359 : ! Passive scalar local variables
360 :
361 : real( kind = core_rknd ), dimension(ngrdcol,nz,sclr_dim) :: &
362 35740224 : sclr1, sclr2, &
363 35740224 : varnce_sclr1, varnce_sclr2, &
364 35740224 : alpha_sclr, &
365 35740224 : corr_sclr_thl_1, corr_sclr_thl_2, &
366 35740224 : corr_sclr_rt_1, corr_sclr_rt_2, &
367 35740224 : corr_w_sclr_1, corr_w_sclr_2
368 :
369 : logical :: &
370 : l_scalar_calc, & ! True if sclr_dim > 0
371 : l_calc_ice_supersat_frac ! True if we should calculate ice_supersat_frac
372 :
373 : ! Quantities needed to predict higher order moments
374 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
375 35740224 : tl1, tl2
376 :
377 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
378 35740224 : sqrt_wp2, & ! Square root of wp2 [m/s]
379 35740224 : Skthl, & ! Skewness of thl [-]
380 35740224 : Skrt, & ! Skewness of rt [-]
381 35740224 : Sku, & ! Skewness of u [-]
382 35740224 : Skv ! Skewness of v [-]
383 :
384 : real( kind = core_rknd ), dimension(ngrdcol,nz,sclr_dim) :: &
385 35740224 : Sksclr ! Skewness of rt [-]
386 :
387 : ! Thermodynamic quantity
388 :
389 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
390 35740224 : wprcp_contrib_comp_1, & ! <w'rc'> contrib. (1st PDF comp.) [m/s(kg/kg)]
391 35740224 : wprcp_contrib_comp_2, & ! <w'rc'> contrib. (2nd PDF comp.) [m/s(kg/kg)]
392 35740224 : wp2rcp_contrib_comp_1, & ! <w'^2rc'> contrib. (1st comp) [m^2/s^2(kg/kg)]
393 35740224 : wp2rcp_contrib_comp_2, & ! <w'^2rc'> contrib. (2nd comp) [m^2/s^2(kg/kg)]
394 35740224 : rtprcp_contrib_comp_1, & ! <rt'rc'> contrib. (1st PDF comp.) [kg^2/kg^2]
395 35740224 : rtprcp_contrib_comp_2, & ! <rt'rc'> contrib. (2nd PDF comp.) [kg^2/kg^2]
396 35740224 : thlprcp_contrib_comp_1, & ! <thl'rc'> contrib. (1st PDF comp.) [K(kg/kg)]
397 35740224 : thlprcp_contrib_comp_2, & ! <thl'rc'> contrib. (2nd PDF comp.) [K(kg/kg)]
398 35740224 : uprcp_contrib_comp_1, & ! <u'rc'> contrib. (1st PDF comp.) [m/s(kg/kg)]
399 35740224 : uprcp_contrib_comp_2, & ! <u'rc'> contrib. (2nd PDF comp.) [m/s(kg/kg)]
400 35740224 : vprcp_contrib_comp_1, & ! <v'rc'> contrib. (1st PDF comp.) [m/s(kg/kg)]
401 35740224 : vprcp_contrib_comp_2 ! <v'rc'> contrib. (2nd PDF comp.) [m/s(kg/kg)]
402 :
403 : ! variables for computing ice cloud fraction
404 : real( kind = core_rknd), dimension(ngrdcol,nz) :: &
405 35740224 : rc_1_ice, rc_2_ice
406 :
407 : ! To test pdf parameters
408 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
409 35740224 : wm_clubb_pdf, &
410 35740224 : rtm_clubb_pdf, &
411 35740224 : thlm_clubb_pdf, &
412 35740224 : wp2_clubb_pdf, &
413 35740224 : rtp2_clubb_pdf, &
414 35740224 : thlp2_clubb_pdf, &
415 35740224 : wp3_clubb_pdf, &
416 35740224 : rtp3_clubb_pdf, &
417 35740224 : thlp3_clubb_pdf, &
418 35740224 : Skw_clubb_pdf, &
419 35740224 : Skrt_clubb_pdf, &
420 35740224 : Skthl_clubb_pdf, &
421 35740224 : rsatl_1, &
422 35740224 : rsatl_2
423 :
424 : real( kind = core_rknd ) :: &
425 : Skw_denom_coef ! CLUBB tunable parameter beta
426 :
427 : logical, parameter :: &
428 : l_liq_ice_loading_test = .false. ! Temp. flag liq./ice water loading test
429 :
430 : integer :: k, i, j, hm_idx ! Indices
431 :
432 : #ifdef GFDL
433 : real ( kind = core_rknd ), parameter :: t1_combined = 273.16, &
434 : t2_combined = 268.16, &
435 : t3_combined = 238.16
436 : #endif
437 :
438 : !----------------------------- Begin Code -----------------------------
439 :
440 : !$acc enter data create( u_1, u_2, varnce_u_1, varnce_u_2, v_1, v_2, &
441 : !$acc varnce_v_1, varnce_v_2, alpha_u, alpha_v, &
442 : !$acc corr_u_w_1, corr_u_w_2, corr_v_w_1, corr_v_w_2, &
443 : !$acc tl1, tl2, sqrt_wp2, Skthl, &
444 : !$acc Skrt, Sku, Skv, wprcp_contrib_comp_1, wprcp_contrib_comp_2, &
445 : !$acc wp2rcp_contrib_comp_1, wp2rcp_contrib_comp_2, &
446 : !$acc rtprcp_contrib_comp_1, rtprcp_contrib_comp_2, &
447 : !$acc thlprcp_contrib_comp_1, thlprcp_contrib_comp_2, &
448 : !$acc uprcp_contrib_comp_1, uprcp_contrib_comp_2, &
449 : !$acc vprcp_contrib_comp_1, vprcp_contrib_comp_2, &
450 : !$acc rc_1_ice, rc_2_ice, rsatl_1, rsatl_2 )
451 :
452 : !$acc enter data if( sclr_dim > 0 ) &
453 : !$acc create( sclr1, sclr2, varnce_sclr1, varnce_sclr2, &
454 : !$acc alpha_sclr, corr_sclr_thl_1, corr_sclr_thl_2, &
455 : !$acc corr_sclr_rt_1, corr_sclr_rt_2, corr_w_sclr_1, &
456 : !$acc corr_w_sclr_2, Sksclr )
457 :
458 : ! Check whether the passive scalars are present.
459 17870112 : if ( sclr_dim > 0 ) then
460 0 : l_scalar_calc = .true.
461 : else
462 17870112 : l_scalar_calc = .false.
463 : end if
464 :
465 : ! Initialize to default values to prevent a runtime error
466 17870112 : if ( ( iiPDF_type /= iiPDF_ADG1 ) .and. ( iiPDF_type /= iiPDF_ADG2 ) ) then
467 :
468 0 : do k = 1, nz
469 0 : do i = 1, ngrdcol
470 0 : pdf_params%alpha_thl(i,k) = one_half
471 0 : pdf_params%alpha_rt(i,k) = one_half
472 : end do
473 : end do
474 :
475 : ! This allows for skewness to be clipped locally without passing the updated
476 : ! value back out.
477 0 : do k = 1, nz
478 0 : do i = 1, ngrdcol
479 0 : Skrt(i,k) = Skrt_in(i,k)
480 0 : Skthl(i,k) = Skthl_in(i,k)
481 0 : Sku(i,k) = Sku_in(i,k)
482 0 : Skv(i,k) = Skv_in(i,k)
483 : end do
484 : end do
485 :
486 0 : do j = 1, sclr_dim
487 0 : do k = 1, nz
488 0 : do i = 1, ngrdcol
489 :
490 0 : Sksclr(i,k,j) = Sksclr_in(i,k,j)
491 :
492 0 : if ( l_scalar_calc ) then
493 0 : alpha_sclr(i,k,j) = one_half
494 : end if
495 :
496 : end do
497 : end do
498 : end do
499 :
500 : end if
501 :
502 : ! Initialize to 0 to prevent a runtime error
503 17870112 : if ( iiPDF_type /= iiPDF_new .and. iiPDF_type /= iiPDF_new_hybrid ) then
504 : ! Stats only variables, setting to zero
505 1536829632 : do k = 1, nz
506 25380961632 : do i = 1, ngrdcol
507 23844132000 : F_w(i,k) = zero
508 23844132000 : F_rt(i,k) = zero
509 23844132000 : F_thl(i,k) = zero
510 23844132000 : min_F_w(i,k) = zero
511 23844132000 : max_F_w(i,k) = zero
512 23844132000 : min_F_rt(i,k) = zero
513 23844132000 : max_F_rt(i,k) = zero
514 23844132000 : min_F_thl(i,k) = zero
515 25363091520 : max_F_thl(i,k) = zero
516 : end do
517 : end do
518 : end if
519 :
520 : ! To avoid recomputing
521 : !$acc parallel loop gang vector collapse(2) default(present)
522 1536829632 : do k = 1, nz
523 25380961632 : do i = 1, ngrdcol
524 25363091520 : sqrt_wp2(i,k) = sqrt( wp2(i,k) )
525 : end do
526 : end do
527 : !$acc end parallel loop
528 :
529 : ! Select the PDF closure method for the two-component PDF used by CLUBB for
530 : ! w, rt, theta-l, and passive scalar variables.
531 : ! Calculate the mixture fraction for the multivariate PDF, as well as both
532 : ! PDF component means and both PDF component variances for each of w, rt,
533 : ! theta-l, and passive scalar variables.
534 : if ( iiPDF_type == iiPDF_ADG1 ) then ! use ADG1
535 :
536 : call ADG1_pdf_driver( nz, ngrdcol, sclr_dim, sclr_tol, & ! In
537 : wm, rtm, thlm, um, vm, & ! In
538 : wp2, rtp2, thlp2, up2, vp2, & ! In
539 : Skw, wprtp, wpthlp, upwp, vpwp, sqrt_wp2, & ! In
540 : sigma_sqd_w, clubb_params(:,ibeta), mixt_frac_max_mag, & ! In
541 : sclrm, sclrp2, wpsclrp, l_scalar_calc, & ! In
542 : pdf_params%w_1, pdf_params%w_2, & ! Out
543 : pdf_params%rt_1, pdf_params%rt_2, & ! Out
544 : pdf_params%thl_1, pdf_params%thl_2, & ! Out
545 : u_1, u_2, v_1, v_2, & ! Out
546 : pdf_params%varnce_w_1, pdf_params%varnce_w_2, & ! Out
547 : pdf_params%varnce_rt_1, pdf_params%varnce_rt_2, & ! Out
548 : pdf_params%varnce_thl_1, pdf_params%varnce_thl_2, & ! Out
549 : varnce_u_1, varnce_u_2, & ! Out
550 : varnce_v_1, varnce_v_2, & ! Out
551 : pdf_params%mixt_frac, & ! Out
552 : pdf_params%alpha_rt, pdf_params%alpha_thl, & ! Out
553 : alpha_u, alpha_v, & ! Out
554 : sclr1, sclr2, varnce_sclr1, & ! Out
555 17870112 : varnce_sclr2, alpha_sclr ) ! Out
556 :
557 : elseif ( iiPDF_type == iiPDF_ADG2 ) then ! use ADG2
558 :
559 : call ADG2_pdf_driver( nz, ngrdcol, sclr_dim, sclr_tol, & ! In
560 : wm, rtm, thlm, wp2, rtp2, thlp2, & ! In
561 : Skw, wprtp, wpthlp, sqrt_wp2, clubb_params(:,ibeta), & ! In
562 : sclrm, sclrp2, wpsclrp, l_scalar_calc, & ! In
563 : pdf_params%w_1, pdf_params%w_2, & ! Out
564 : pdf_params%rt_1, pdf_params%rt_2, & ! Out
565 : pdf_params%thl_1, pdf_params%thl_2, & ! Out
566 : pdf_params%varnce_w_1, pdf_params%varnce_w_2, & ! Out
567 : pdf_params%varnce_rt_1, pdf_params%varnce_rt_2, & ! Out
568 : pdf_params%varnce_thl_1, pdf_params%varnce_thl_2, & ! Out
569 : pdf_params%mixt_frac, & ! Out
570 : pdf_params%alpha_rt, pdf_params%alpha_thl, & ! Out
571 : sigma_sqd_w, sclr1, sclr2, & ! Out
572 0 : varnce_sclr1, varnce_sclr2, alpha_sclr ) ! Out
573 :
574 : elseif ( iiPDF_type == iiPDF_3D_Luhar ) then ! use 3D Luhar
575 0 : do i = 1, ngrdcol
576 : call Luhar_3D_pdf_driver( nz, &
577 0 : wm(i,:), rtm(i,:), thlm(i,:), wp2(i,:), rtp2(i,:), thlp2(i,:), & ! In
578 0 : Skw(i,:), Skrt(i,:), Skthl(i,:), wprtp(i,:), wpthlp(i,:), & ! In
579 0 : pdf_params%w_1(i,:), pdf_params%w_2(i,:), & ! Out
580 0 : pdf_params%rt_1(i,:), pdf_params%rt_2(i,:), & ! Out
581 0 : pdf_params%thl_1(i,:), pdf_params%thl_2(i,:), & ! Out
582 0 : pdf_params%varnce_w_1(i,:), pdf_params%varnce_w_2(i,:), & ! Out
583 0 : pdf_params%varnce_rt_1(i,:), pdf_params%varnce_rt_2(i,:), & ! Out
584 0 : pdf_params%varnce_thl_1(i,:), pdf_params%varnce_thl_2(i,:), & ! Out
585 0 : pdf_params%mixt_frac(i,:) ) ! Out
586 : end do
587 : elseif ( iiPDF_type == iiPDF_new ) then ! use new PDF
588 : call new_pdf_driver( nz, ngrdcol, wm, rtm, thlm, wp2, rtp2, thlp2, Skw, & ! In
589 : wprtp, wpthlp, rtpthlp, & ! In
590 : clubb_params, & ! In
591 : Skrt, Skthl, & ! In/Out
592 : pdf_params%w_1, pdf_params%w_2, & ! Out
593 : pdf_params%rt_1, pdf_params%rt_2, & ! Out
594 : pdf_params%thl_1, pdf_params%thl_2, & ! Out
595 : pdf_params%varnce_w_1, pdf_params%varnce_w_2, & ! Out
596 : pdf_params%varnce_rt_1, pdf_params%varnce_rt_2, & ! Out
597 : pdf_params%varnce_thl_1, pdf_params%varnce_thl_2, & ! Out
598 : pdf_params%mixt_frac, & ! Out
599 : pdf_implicit_coefs_terms, & ! Out
600 : F_w, F_rt, F_thl, min_F_w, max_F_w, & ! Out
601 0 : min_F_rt, max_F_rt, min_F_thl, max_F_thl ) ! Out
602 : elseif ( iiPDF_type == iiPDF_TSDADG ) then
603 0 : do i = 1, ngrdcol
604 : call tsdadg_pdf_driver( nz, &
605 0 : wm(i,:), rtm(i,:), thlm(i,:), wp2(i,:), rtp2(i,:), thlp2(i,:), & ! In
606 0 : Skw(i,:), Skrt(i,:), Skthl(i,:), wprtp(i,:), wpthlp(i,:), & ! In
607 0 : pdf_params%w_1(i,:), pdf_params%w_2(i,:), & ! Out
608 0 : pdf_params%rt_1(i,:), pdf_params%rt_2(i,:), & ! Out
609 0 : pdf_params%thl_1(i,:), pdf_params%thl_2(i,:), & ! Out
610 0 : pdf_params%varnce_w_1(i,:), pdf_params%varnce_w_2(i,:), & ! Out
611 0 : pdf_params%varnce_rt_1(i,:), pdf_params%varnce_rt_2(i,:), & ! Out
612 0 : pdf_params%varnce_thl_1(i,:), pdf_params%varnce_thl_2(i,:), & ! Out
613 0 : pdf_params%mixt_frac(i,:) ) ! Out
614 : end do
615 : elseif ( iiPDF_type == iiPDF_LY93 ) then ! use LY93
616 0 : do i = 1, ngrdcol
617 0 : call LY93_driver( nz, wm(i,:), rtm(i,:), thlm(i,:), wp2(i,:), rtp2(i,:), & ! In
618 0 : thlp2(i,:), Skw(i,:), Skrt(i,:), Skthl(i,:), & ! In
619 0 : pdf_params%w_1(i,:), pdf_params%w_2(i,:), & ! Out
620 0 : pdf_params%rt_1(i,:), pdf_params%rt_2(i,:), & ! Out
621 0 : pdf_params%thl_1(i,:), pdf_params%thl_2(i,:), & ! Out
622 0 : pdf_params%varnce_w_1(i,:), pdf_params%varnce_w_2(i,:), & ! Out
623 0 : pdf_params%varnce_rt_1(i,:), pdf_params%varnce_rt_2(i,:), & ! Out
624 0 : pdf_params%varnce_thl_1(i,:), pdf_params%varnce_thl_2(i,:), & ! Out
625 0 : pdf_params%mixt_frac(i,:) ) ! Out
626 : end do
627 : elseif ( iiPDF_type == iiPDF_new_hybrid ) then ! use new hybrid PDF
628 : call new_hybrid_pdf_driver( nz, ngrdcol, sclr_dim, & ! In
629 : wm, rtm, thlm, um, vm, & ! In
630 : wp2, rtp2, thlp2, up2, vp2, & ! In
631 : Skw, wprtp, wpthlp, upwp, vpwp, & ! In
632 : sclrm, sclrp2, wpsclrp, & ! In
633 : gamma_Skw_fnc, & ! In
634 : clubb_params(:,islope_coef_spread_DG_means_w), & ! In
635 : clubb_params(:,ipdf_component_stdev_factor_w), & ! In
636 : Skrt, Skthl, Sku, Skv, Sksclr, & ! I/O
637 : pdf_params%w_1, pdf_params%w_2, & ! Out
638 : pdf_params%rt_1, pdf_params%rt_2, & ! Out
639 : pdf_params%thl_1, pdf_params%thl_2, & ! Out
640 : u_1, u_2, v_1, v_2, & ! Out
641 : pdf_params%varnce_w_1, & ! Out
642 : pdf_params%varnce_w_2, & ! Out
643 : pdf_params%varnce_rt_1, & ! Out
644 : pdf_params%varnce_rt_2, & ! Out
645 : pdf_params%varnce_thl_1, & ! Out
646 : pdf_params%varnce_thl_2, & ! Out
647 : varnce_u_1, varnce_u_2, & ! Out
648 : varnce_v_1, varnce_v_2, & ! Out
649 : sclr1, sclr2, & ! Out
650 : varnce_sclr1, varnce_sclr2, & ! Out
651 : pdf_params%mixt_frac, & ! Out
652 : pdf_implicit_coefs_terms, & ! Out
653 0 : F_w, min_F_w, max_F_w ) ! Out
654 :
655 : ! The calculation of skewness of rt, thl, u, v, and scalars is hard-wired
656 : ! for use with the ADG1 code, which contains the variable sigma_sqd_w.
657 : ! In order to use an equivalent expression for these skewnesses using the
658 : ! new hybrid PDF (without doing more recoding), set the value of
659 : ! sigma_sqd_w to 1 - F_w.
660 0 : do k = 1, nz
661 0 : do i = 1, ngrdcol
662 0 : sigma_sqd_w(i,k) = one - F_w(i,k)
663 : end do
664 : end do
665 :
666 : end if ! iiPDF_type
667 :
668 : ! Calculate the PDF component correlations of rt and thl.
669 : call calc_comp_corrs_binormal( nz, ngrdcol, & ! In
670 : rtpthlp, rtm, thlm, & ! In
671 : pdf_params%rt_1, pdf_params%rt_2, & ! In
672 : pdf_params%thl_1, pdf_params%thl_2, & ! In
673 : pdf_params%varnce_rt_1, pdf_params%varnce_rt_2, & ! In
674 : pdf_params%varnce_thl_1, pdf_params%varnce_thl_2, & ! In
675 : pdf_params%mixt_frac, & ! In
676 17870112 : pdf_params%corr_rt_thl_1, pdf_params%corr_rt_thl_2 ) ! Out
677 :
678 : if ( iiPDF_type == iiPDF_ADG1 .or. iiPDF_type == iiPDF_ADG2 &
679 17870112 : .or. iiPDF_type == iiPDF_new_hybrid ) then
680 :
681 : ! These PDF types define corr_w_rt_1, corr_w_rt_2, corr_w_thl_1, and
682 : ! corr_w_thl_2 to all have a value of 0, so skip the calculation.
683 : ! The values of corr_u_w_1, corr_u_w_2, corr_v_w_1, and corr_v_w_2 are
684 : ! all defined to be 0, as well.
685 : !$acc parallel loop gang vector collapse(2) default(present)
686 1536829632 : do k = 1, nz
687 25380961632 : do i = 1, ngrdcol
688 23844132000 : pdf_params%corr_w_rt_1(i,k) = zero
689 23844132000 : pdf_params%corr_w_rt_2(i,k) = zero
690 23844132000 : pdf_params%corr_w_thl_1(i,k) = zero
691 23844132000 : pdf_params%corr_w_thl_2(i,k) = zero
692 23844132000 : corr_u_w_1(i,k) = zero
693 23844132000 : corr_u_w_2(i,k) = zero
694 23844132000 : corr_v_w_1(i,k) = zero
695 25363091520 : corr_v_w_2(i,k) = zero
696 : end do
697 : end do
698 : !$acc end parallel loop
699 :
700 : else
701 :
702 : ! Calculate the PDF component correlations of w and rt.
703 : call calc_comp_corrs_binormal( nz, ngrdcol, & ! In
704 : wprtp, wm, rtm, & ! In
705 : pdf_params%w_1, pdf_params%w_2, & ! In
706 : pdf_params%rt_1, pdf_params%rt_2, & ! In
707 : pdf_params%varnce_w_1, pdf_params%varnce_w_2, & ! In
708 : pdf_params%varnce_rt_1, pdf_params%varnce_rt_2, & ! In
709 : pdf_params%mixt_frac, & ! In
710 0 : pdf_params%corr_w_rt_1, pdf_params%corr_w_rt_2 ) ! Out
711 :
712 : ! Calculate the PDF component correlations of w and thl.
713 : call calc_comp_corrs_binormal( nz, ngrdcol, & ! In
714 : wpthlp, wm, thlm, & ! In
715 : pdf_params%w_1, pdf_params%w_2, & ! In
716 : pdf_params%thl_1, pdf_params%thl_2, & ! In
717 : pdf_params%varnce_w_1, pdf_params%varnce_w_2, & ! In
718 : pdf_params%varnce_thl_1, pdf_params%varnce_thl_2, & ! In
719 : pdf_params%mixt_frac, & ! In
720 0 : pdf_params%corr_w_thl_1, pdf_params%corr_w_thl_2 ) ! Out
721 : end if
722 :
723 17870112 : if ( l_scalar_calc ) then
724 :
725 : ! Calculate the PDF component correlations of a passive scalar and thl.
726 0 : do j = 1, sclr_dim
727 : call calc_comp_corrs_binormal( nz, ngrdcol, & ! In
728 : sclrpthlp(:,:,j), sclrm(:,:,j), thlm, & ! In
729 : sclr1(:,:,j), sclr2(:,:,j), & ! In
730 : pdf_params%thl_1, pdf_params%thl_2, & ! In
731 : varnce_sclr1(:,:,j), varnce_sclr2(:,:,j), & ! In
732 : pdf_params%varnce_thl_1, pdf_params%varnce_thl_2, & ! In
733 : pdf_params%mixt_frac, & ! In
734 0 : corr_sclr_thl_1(:,:,j), corr_sclr_thl_2(:,:,j) ) ! Out
735 : end do
736 :
737 : ! Calculate the PDF component correlations of a passive scalar and rt.
738 0 : do j = 1, sclr_dim
739 : call calc_comp_corrs_binormal( nz, ngrdcol, & ! In
740 : sclrprtp(:,:,j), sclrm(:,:,j), rtm, & ! In
741 : sclr1(:,:,j), sclr2(:,:,j), & ! In
742 : pdf_params%rt_1, pdf_params%rt_2, & ! In
743 : varnce_sclr1(:,:,j), varnce_sclr2(:,:,j), & ! In
744 : pdf_params%varnce_rt_1, pdf_params%varnce_rt_2, & ! In
745 : pdf_params%mixt_frac, & ! In
746 0 : corr_sclr_rt_1(:,:,j), corr_sclr_rt_2(:,:,j) ) ! Out
747 : end do
748 :
749 : if ( iiPDF_type == iiPDF_ADG1 .or. iiPDF_type == iiPDF_ADG2 &
750 0 : .or. iiPDF_type == iiPDF_new_hybrid ) then
751 :
752 : ! These PDF types define all PDF component correlations involving w
753 : ! to have a value of 0, so skip the calculation.
754 : !$acc parallel loop gang vector collapse(2) default(present)
755 0 : do j = 1, sclr_dim
756 0 : do k = 1, nz
757 0 : do i = 1, ngrdcol
758 0 : corr_w_sclr_1(i,k,j) = zero
759 0 : corr_w_sclr_2(i,k,j) = zero
760 : end do
761 : end do
762 : end do
763 : !$acc end parallel loop
764 :
765 : else
766 :
767 : ! Calculate the PDF component correlations of w and a passive scalar.
768 0 : do j = 1, sclr_dim
769 : call calc_comp_corrs_binormal( nz, ngrdcol, & ! In
770 : wpsclrp(:,:,j), wm, sclrm(:,:,j), & ! In
771 : pdf_params%w_1, pdf_params%w_2, & ! In
772 : sclr1(:,:,j), sclr2(:,:,j), & ! In
773 : pdf_params%varnce_w_1, pdf_params%varnce_w_2, & ! In
774 : varnce_sclr1(:,:,j), varnce_sclr2(:,:,j), & ! In
775 : pdf_params%mixt_frac, & ! In
776 0 : corr_w_sclr_1(:,:,j), corr_w_sclr_2(:,:,j) ) ! Out
777 :
778 : end do
779 :
780 : end if
781 :
782 : end if
783 :
784 :
785 : ! Compute higher order moments (these are interactive)
786 : call calc_wp2xp_pdf( nz, ngrdcol, &
787 : wm, rtm, pdf_params%w_1, pdf_params%w_2, &
788 : pdf_params%rt_1, pdf_params%rt_2, &
789 : pdf_params%varnce_w_1, pdf_params%varnce_w_2, &
790 : pdf_params%varnce_rt_1, pdf_params%varnce_rt_2, &
791 : pdf_params%corr_w_rt_1, pdf_params%corr_w_rt_2, &
792 : pdf_params%mixt_frac, &
793 17870112 : wp2rtp )
794 :
795 : call calc_wp2xp_pdf( nz, ngrdcol, &
796 : wm, thlm, pdf_params%w_1, pdf_params%w_2, &
797 : pdf_params%thl_1, pdf_params%thl_2, &
798 : pdf_params%varnce_w_1, pdf_params%varnce_w_2, &
799 : pdf_params%varnce_thl_1, pdf_params%varnce_thl_2, &
800 : pdf_params%corr_w_thl_1, pdf_params%corr_w_thl_2, &
801 : pdf_params%mixt_frac, &
802 17870112 : wp2thlp )
803 :
804 : ! Compute higher order moments (these may be interactive)
805 : call calc_wpxp2_pdf( nz, ngrdcol, &
806 : wm, um, pdf_params%w_1, pdf_params%w_2, &
807 : u_1, u_2, &
808 : pdf_params%varnce_w_1, pdf_params%varnce_w_2, &
809 : varnce_u_1, varnce_u_2, &
810 : corr_u_w_1, corr_u_w_2, &
811 : pdf_params%mixt_frac, &
812 17870112 : wpup2 )
813 :
814 : call calc_wpxp2_pdf( nz, ngrdcol, &
815 : wm, vm, pdf_params%w_1, pdf_params%w_2, &
816 : v_1, v_2, &
817 : pdf_params%varnce_w_1, pdf_params%varnce_w_2, &
818 : varnce_v_1, varnce_v_2, &
819 : corr_v_w_1, corr_v_w_2, &
820 : pdf_params%mixt_frac, &
821 17870112 : wpvp2 )
822 :
823 : call calc_wp2xp2_pdf( nz, ngrdcol, &
824 : wm, um, pdf_params%w_1, pdf_params%w_2, &
825 : u_1, u_2, &
826 : pdf_params%varnce_w_1, pdf_params%varnce_w_2, &
827 : varnce_u_1, varnce_u_2, &
828 : corr_u_w_1, corr_u_w_2, &
829 : pdf_params%mixt_frac, &
830 17870112 : wp2up2 )
831 :
832 : call calc_wp2xp2_pdf( nz, ngrdcol, &
833 : wm, vm, pdf_params%w_1, pdf_params%w_2, &
834 : v_1, v_2, &
835 : pdf_params%varnce_w_1, pdf_params%varnce_w_2, &
836 : varnce_v_1, varnce_v_2, &
837 : corr_v_w_1, corr_v_w_2, &
838 : pdf_params%mixt_frac, &
839 17870112 : wp2vp2 )
840 :
841 : call calc_wp4_pdf( nz, ngrdcol, &
842 : wm, pdf_params%w_1, pdf_params%w_2, &
843 : pdf_params%varnce_w_1, pdf_params%varnce_w_2, &
844 : pdf_params%mixt_frac, &
845 17870112 : wp4 )
846 :
847 17870112 : if ( l_explicit_turbulent_adv_xpyp .or. stats_metadata%iwprtp2 > 0 ) then
848 : call calc_wpxp2_pdf( nz, ngrdcol, &
849 : wm, rtm, pdf_params%w_1, pdf_params%w_2, &
850 : pdf_params%rt_1, pdf_params%rt_2, &
851 : pdf_params%varnce_w_1, pdf_params%varnce_w_2, &
852 : pdf_params%varnce_rt_1, pdf_params%varnce_rt_2, &
853 : pdf_params%corr_w_rt_1, pdf_params%corr_w_rt_2, &
854 : pdf_params%mixt_frac, &
855 0 : wprtp2 )
856 : end if
857 :
858 17870112 : if ( l_explicit_turbulent_adv_xpyp .or. stats_metadata%iwpthlp2 > 0 ) then
859 : call calc_wpxp2_pdf( nz, ngrdcol, &
860 : wm, thlm, pdf_params%w_1, pdf_params%w_2, &
861 : pdf_params%thl_1, pdf_params%thl_2, &
862 : pdf_params%varnce_w_1, pdf_params%varnce_w_2, &
863 : pdf_params%varnce_thl_1, pdf_params%varnce_thl_2, &
864 : pdf_params%corr_w_thl_1, pdf_params%corr_w_thl_2, &
865 : pdf_params%mixt_frac, &
866 0 : wpthlp2 )
867 : end if
868 :
869 17870112 : if ( l_explicit_turbulent_adv_xpyp .or. stats_metadata%iwprtpthlp > 0 ) then
870 :
871 : call calc_wpxpyp_pdf( nz, ngrdcol, &
872 : wm, rtm, thlm, pdf_params%w_1, pdf_params%w_2, &
873 : pdf_params%rt_1, pdf_params%rt_2, &
874 : pdf_params%thl_1, pdf_params%thl_2, &
875 : pdf_params%varnce_w_1, pdf_params%varnce_w_2, &
876 : pdf_params%varnce_rt_1, pdf_params%varnce_rt_2, &
877 : pdf_params%varnce_thl_1, pdf_params%varnce_thl_2, &
878 : pdf_params%corr_w_rt_1, pdf_params%corr_w_rt_2, &
879 : pdf_params%corr_w_thl_1, pdf_params%corr_w_thl_2, &
880 : pdf_params%corr_rt_thl_1, pdf_params%corr_rt_thl_2, &
881 : pdf_params%mixt_frac, &
882 0 : wprtpthlp )
883 : end if
884 :
885 :
886 : ! Scalar Addition to higher order moments
887 17870112 : if ( l_scalar_calc ) then
888 :
889 0 : do j = 1, sclr_dim
890 : call calc_wp2xp_pdf( nz, ngrdcol, &
891 : wm, sclrm(:,:,j), pdf_params%w_1, pdf_params%w_2, &
892 : sclr1(:,:,j), sclr2(:,:,j), &
893 : pdf_params%varnce_w_1, pdf_params%varnce_w_2, &
894 : varnce_sclr1(:,:,j), varnce_sclr2(:,:,j), &
895 : corr_w_sclr_1(:,:,j), corr_w_sclr_2(:,:,j), &
896 : pdf_params%mixt_frac, &
897 0 : wp2sclrp(:,:,j) )
898 : end do
899 :
900 0 : do j = 1, sclr_dim
901 : call calc_wpxp2_pdf( nz, ngrdcol, &
902 : wm, sclrm(:,:,j), pdf_params%w_1, pdf_params%w_2, &
903 : sclr1(:,:,j), sclr2(:,:,j), &
904 : pdf_params%varnce_w_1, pdf_params%varnce_w_2, &
905 : varnce_sclr1(:,:,j), varnce_sclr2(:,:,j), &
906 : corr_w_sclr_1(:,:,j), corr_w_sclr_2(:,:,j), &
907 : pdf_params%mixt_frac, &
908 0 : wpsclrp2(:,:,j) )
909 : end do
910 :
911 0 : do j = 1, sclr_dim
912 : call calc_wpxpyp_pdf( nz, ngrdcol, &
913 : wm, sclrm(:,:,j), rtm, pdf_params%w_1, pdf_params%w_2, &
914 : sclr1(:,:,j), sclr2(:,:,j), &
915 : pdf_params%rt_1, pdf_params%rt_2, &
916 : pdf_params%varnce_w_1, pdf_params%varnce_w_2, &
917 : varnce_sclr1(:,:,j), varnce_sclr2(:,:,j), &
918 : pdf_params%varnce_rt_1, pdf_params%varnce_rt_2, &
919 : corr_w_sclr_1(:,:,j), corr_w_sclr_2(:,:,j), &
920 : pdf_params%corr_w_rt_1, pdf_params%corr_w_rt_2, &
921 : corr_sclr_rt_1(:,:,j), corr_sclr_rt_2(:,:,j), &
922 : pdf_params%mixt_frac, &
923 0 : wpsclrprtp(:,:,j) )
924 : end do
925 :
926 0 : do j = 1, sclr_dim
927 : call calc_wpxpyp_pdf( nz, ngrdcol, &
928 : wm, sclrm(:,:,j), thlm, pdf_params%w_1, pdf_params%w_2, &
929 : sclr1(:,:,j), sclr2(:,:,j), &
930 : pdf_params%thl_1, pdf_params%thl_2, &
931 : pdf_params%varnce_w_1, pdf_params%varnce_w_2, &
932 : varnce_sclr1(:,:,j), varnce_sclr2(:,:,j), &
933 : pdf_params%varnce_thl_1, pdf_params%varnce_thl_2, &
934 : corr_w_sclr_1(:,:,j), corr_w_sclr_2(:,:,j), &
935 : pdf_params%corr_w_thl_1, pdf_params%corr_w_thl_2, &
936 : corr_sclr_thl_1(:,:,j), corr_sclr_thl_2(:,:,j), &
937 : pdf_params%mixt_frac, &
938 0 : wpsclrpthlp(:,:,j) )
939 : end do
940 :
941 : end if
942 :
943 : ! Compute higher order moments that include theta_v.
944 :
945 : ! First compute some preliminary quantities.
946 : ! "1" denotes first Gaussian; "2" denotes 2nd Gaussian
947 : ! liq water temp (Sommeria & Deardorff 1977 (SD), eqn. 3)
948 : !$acc parallel loop gang vector collapse(2) default(present)
949 1536829632 : do k = 1, nz
950 25380961632 : do i = 1, ngrdcol
951 23844132000 : tl1(i,k) = pdf_params%thl_1(i,k)*exner(i,k)
952 25363091520 : tl2(i,k) = pdf_params%thl_2(i,k)*exner(i,k)
953 : end do
954 : end do
955 : !$acc end parallel loop
956 :
957 : #ifdef GFDL
958 : if ( sclr_dim > 0 .and. (.not. do_liquid_only_in_clubb) ) then ! h1g, 2010-06-16 begin mod
959 :
960 : do i = 1, ngrdcol
961 : where ( tl1(i,:) > t1_combined )
962 : pdf_params%rsatl_1(i,:) = sat_mixrat_liq( nz, p_in_Pa(i,:), tl1(i,:) )
963 : elsewhere ( tl1(i,:) > t2_combined )
964 : pdf_params%rsatl_1(i,:) = sat_mixrat_liq( nz, p_in_Pa(i,:), tl1(i,:) ) &
965 : * (tl1(i,:) - t2_combined)/(t1_combined - t2_combined) &
966 : + sat_mixrat_ice( nz, p_in_Pa(i,:), tl1(i,:) ) &
967 : * (t1_combined - tl1(i,:))/(t1_combined - t2_combined)
968 : elsewhere ( tl1(i,:) > t3_combined )
969 : pdf_params%rsatl_1(i,:) = sat_mixrat_ice( nz, p_in_Pa(i,:), tl1(i,:) ) &
970 : + sat_mixrat_ice( nz, p_in_Pa(i,:), tl1(i,:) ) * (RH_crit(i, :, 1, 1) -one ) &
971 : * ( t2_combined -tl1(i,:))/(t2_combined - t3_combined)
972 : elsewhere
973 : pdf_params%rsatl_1(i,:) = sat_mixrat_ice( nz, p_in_Pa(i,:), tl1(i,:) ) * RH_crit(i, :, 1, 1)
974 : endwhere
975 :
976 : where ( tl2(i,:) > t1_combined )
977 : pdf_params%rsatl_2(i,:) = sat_mixrat_liq( nz, p_in_Pa(i,:), tl2(i,:) )
978 : elsewhere ( tl2(i,:) > t2_combined )
979 : pdf_params%rsatl_2(i,:) = sat_mixrat_liq( nz, p_in_Pa(i,:), tl2(i,:) ) &
980 : * (tl2(i,:) - t2_combined)/(t1_combined - t2_combined) &
981 : + sat_mixrat_ice( nz, p_in_Pa(i,:), tl2(i,:) ) &
982 : * (t1_combined - tl2(i,:))/(t1_combined - t2_combined)
983 : elsewhere ( tl2(i,:) > t3_combined )
984 : pdf_params%rsatl_2(i,:) = sat_mixrat_ice( nz, p_in_Pa(i,:), tl2(i,:) ) &
985 : + sat_mixrat_ice( nz, p_in_Pa(i,:), tl2(i,:) )* (RH_crit(i, :, 1, 2) -one) &
986 : * ( t2_combined -tl2(i,:))/(t2_combined - t3_combined)
987 : elsewhere
988 : pdf_params%rsatl_2(i,:) = sat_mixrat_ice( nz, p_in_Pa(i,:), tl2(i,:) ) * RH_crit(i, :, 1, 2)
989 : endwhere
990 :
991 : end do
992 :
993 : else ! sclr_dim <= 0 or do_liquid_only_in_clubb = .T.
994 :
995 : pdf_params%rsatl_1 = sat_mixrat_liq( nz, ngrdcol, p_in_Pa, tl1, saturation_formula )
996 : pdf_params%rsatl_2 = sat_mixrat_liq( nz, ngrdcol, p_in_Pa, tl2, saturation_formula )
997 :
998 : end if !sclr_dim > 0
999 :
1000 : ! Determine whether to compute ice_supersat_frac. We do not compute
1001 : ! ice_supersat_frac for GFDL (unless do_liquid_only_in_clubb is true),
1002 : ! because liquid and ice are both fed into rtm, ruining the calculation.
1003 : if (do_liquid_only_in_clubb) then
1004 : l_calc_ice_supersat_frac = .true.
1005 : else
1006 : l_calc_ice_supersat_frac = .false.
1007 : end if
1008 :
1009 : #else
1010 17870112 : rsatl_1 = sat_mixrat_liq( nz, ngrdcol, p_in_Pa, tl1, saturation_formula )
1011 17870112 : rsatl_2 = sat_mixrat_liq( nz, ngrdcol, p_in_Pa, tl2, saturation_formula ) ! h1g, 2010-06-16 end mod
1012 :
1013 : !$acc parallel loop gang vector collapse(2) default(present)
1014 1536829632 : do k = 1, nz
1015 25380961632 : do i = 1, ngrdcol
1016 23844132000 : pdf_params%rsatl_1(i,k) = rsatl_1(i,k)
1017 25363091520 : pdf_params%rsatl_2(i,k) = rsatl_2(i,k)
1018 : end do
1019 : end do
1020 : !$acc end parallel loop
1021 :
1022 17870112 : l_calc_ice_supersat_frac = .true.
1023 : #endif
1024 :
1025 : call transform_pdf_chi_eta_component( nz, ngrdcol, &
1026 : tl1, pdf_params%rsatl_1, pdf_params%rt_1, exner, & ! In
1027 : pdf_params%varnce_thl_1, pdf_params%varnce_rt_1, & ! In
1028 : pdf_params%corr_rt_thl_1, pdf_params%chi_1, & ! In
1029 : pdf_params%crt_1, pdf_params%cthl_1, & ! Out
1030 : pdf_params%stdev_chi_1, pdf_params%stdev_eta_1, & ! Out
1031 : pdf_params%covar_chi_eta_1, & ! Out
1032 17870112 : pdf_params%corr_chi_eta_1 ) ! Out
1033 :
1034 :
1035 : ! Calculate cloud fraction component for pdf 1
1036 : call calc_liquid_cloud_frac_component( nz, ngrdcol, &
1037 : pdf_params%chi_1, pdf_params%stdev_chi_1, & ! In
1038 17870112 : pdf_params%cloud_frac_1, pdf_params%rc_1 ) ! Out
1039 :
1040 : ! Calc ice_supersat_frac
1041 : if ( l_calc_ice_supersat_frac ) then
1042 :
1043 : call calc_ice_cloud_frac_component( nz, ngrdcol, &
1044 : pdf_params%chi_1, pdf_params%stdev_chi_1, &
1045 : pdf_params%rc_1, pdf_params%cloud_frac_1, &
1046 : p_in_Pa, tl1, &
1047 : pdf_params%rsatl_1, pdf_params%crt_1, &
1048 : saturation_formula, &
1049 17870112 : pdf_params%ice_supersat_frac_1, rc_1_ice )
1050 : end if
1051 :
1052 : call transform_pdf_chi_eta_component( nz, ngrdcol, &
1053 : tl2, pdf_params%rsatl_2, pdf_params%rt_2, exner, & ! In
1054 : pdf_params%varnce_thl_2, pdf_params%varnce_rt_2, & ! In
1055 : pdf_params%corr_rt_thl_2, pdf_params%chi_2, & ! In
1056 : pdf_params%crt_2, pdf_params%cthl_2, & ! Out
1057 : pdf_params%stdev_chi_2, pdf_params%stdev_eta_2, & ! Out
1058 : pdf_params%covar_chi_eta_2, & ! Out
1059 17870112 : pdf_params%corr_chi_eta_2 ) ! Out
1060 :
1061 :
1062 : ! Calculate cloud fraction component for pdf 2
1063 : call calc_liquid_cloud_frac_component( nz, ngrdcol, &
1064 : pdf_params%chi_2, pdf_params%stdev_chi_2, & ! In
1065 17870112 : pdf_params%cloud_frac_2, pdf_params%rc_2 ) ! Out
1066 :
1067 : ! Calc ice_supersat_frac
1068 : if ( l_calc_ice_supersat_frac ) then
1069 :
1070 : call calc_ice_cloud_frac_component( nz, ngrdcol, &
1071 : pdf_params%chi_2, pdf_params%stdev_chi_2, &
1072 : pdf_params%rc_2, pdf_params%cloud_frac_2, &
1073 : p_in_Pa, tl2, &
1074 : pdf_params%rsatl_2, pdf_params%crt_2, &
1075 : saturation_formula, &
1076 17870112 : pdf_params%ice_supersat_frac_2, rc_2_ice )
1077 :
1078 : ! Compute ice cloud fraction, ice_supersat_frac
1079 : !$acc parallel loop gang vector collapse(2) default(present)
1080 1536829632 : do k = 1, nz
1081 25380961632 : do i = 1, ngrdcol
1082 47688264000 : ice_supersat_frac(i,k) = pdf_params%mixt_frac(i,k) &
1083 0 : * pdf_params%ice_supersat_frac_1(i,k) &
1084 : + ( one - pdf_params%mixt_frac(i,k) ) &
1085 73051355520 : * pdf_params%ice_supersat_frac_2(i,k)
1086 : end do
1087 : end do
1088 : !$acc end parallel loop
1089 :
1090 : else
1091 :
1092 : ! ice_supersat_frac will be garbage if computed as above
1093 : !$acc parallel loop gang vector collapse(2) default(present)
1094 : do k = 1, nz
1095 : do i = 1, ngrdcol
1096 : ice_supersat_frac(i,k) = 0.0_core_rknd
1097 : end do
1098 : end do
1099 : !$acc end parallel loop
1100 :
1101 : if (clubb_at_least_debug_level( 1 )) then
1102 : write(fstderr,*) "Warning: ice_supersat_frac has garbage values if &
1103 : & do_liquid_only_in_clubb = .false."
1104 : end if
1105 :
1106 : end if ! l_calc_ice_supersat_frac
1107 :
1108 :
1109 : ! Compute cloud fraction and mean cloud water mixing ratio.
1110 : ! Reference:
1111 : ! https://arxiv.org/pdf/1711.03675v1.pdf#nameddest=url:anl_int_cloud_terms
1112 : !$acc parallel loop gang vector collapse(2) default(present)
1113 1536829632 : do k = 1, nz
1114 25380961632 : do i = 1, ngrdcol
1115 47688264000 : cloud_frac(i,k) = pdf_params%mixt_frac(i,k) * pdf_params%cloud_frac_1(i,k) &
1116 71532396000 : + ( one - pdf_params%mixt_frac(i,k) ) * pdf_params%cloud_frac_2(i,k)
1117 0 : rcm(i,k) = pdf_params%mixt_frac(i,k) * pdf_params%rc_1(i,k) + ( one - pdf_params%mixt_frac(i,k) ) &
1118 23844132000 : * pdf_params%rc_2(i,k)
1119 25363091520 : rcm(i,k) = max( zero_threshold, rcm(i,k) )
1120 : end do
1121 : end do
1122 : !$acc end parallel loop
1123 :
1124 : if ( iiPDF_type == iiPDF_ADG1 .or. iiPDF_type == iiPDF_ADG2 &
1125 17870112 : .or. iiPDF_type == iiPDF_new_hybrid ) then
1126 :
1127 : ! corr_w_rt and corr_w_thl are zero for these pdf types so
1128 : ! corr_w_chi and corr_w_eta are zero as well
1129 : !$acc parallel loop gang vector collapse(2) default(present)
1130 1536829632 : do k = 1, nz
1131 25380961632 : do i = 1, ngrdcol
1132 23844132000 : pdf_params%corr_w_chi_1(i,k) = zero
1133 23844132000 : pdf_params%corr_w_chi_2(i,k) = zero
1134 23844132000 : pdf_params%corr_w_eta_1(i,k) = zero
1135 25363091520 : pdf_params%corr_w_eta_2(i,k) = zero
1136 : end do
1137 : end do
1138 : !$acc end parallel loop
1139 :
1140 : else
1141 :
1142 : ! Correlation of w and chi for each component.
1143 : pdf_params%corr_w_chi_1 &
1144 0 : = calc_corr_chi_x( pdf_params%crt_1, pdf_params%cthl_1, &
1145 0 : sqrt(pdf_params%varnce_rt_1), sqrt(pdf_params%varnce_thl_1), &
1146 0 : pdf_params%stdev_chi_1, &
1147 0 : pdf_params%corr_w_rt_1, pdf_params%corr_w_thl_1 )
1148 :
1149 : pdf_params%corr_w_chi_2 &
1150 0 : = calc_corr_chi_x( pdf_params%crt_2, pdf_params%cthl_2, &
1151 0 : sqrt(pdf_params%varnce_rt_2), sqrt(pdf_params%varnce_thl_2), &
1152 0 : pdf_params%stdev_chi_2, pdf_params%corr_w_rt_2, &
1153 0 : pdf_params%corr_w_thl_2 )
1154 :
1155 : ! Correlation of w and eta for each component.
1156 : pdf_params%corr_w_eta_1 &
1157 0 : = calc_corr_eta_x( pdf_params%crt_1, pdf_params%cthl_1, &
1158 0 : sqrt(pdf_params%varnce_rt_1), sqrt(pdf_params%varnce_thl_1), &
1159 0 : pdf_params%stdev_eta_1, pdf_params%corr_w_rt_1, &
1160 0 : pdf_params%corr_w_thl_1 )
1161 :
1162 : pdf_params%corr_w_eta_2 &
1163 0 : = calc_corr_eta_x( pdf_params%crt_2, pdf_params%cthl_2, &
1164 0 : sqrt(pdf_params%varnce_rt_2), sqrt(pdf_params%varnce_thl_2), &
1165 0 : pdf_params%stdev_eta_2, pdf_params%corr_w_rt_2, &
1166 0 : pdf_params%corr_w_thl_2 )
1167 :
1168 : end if
1169 :
1170 :
1171 : ! Compute moments that depend on theta_v
1172 : !
1173 : ! The moments that depend on th_v' are calculated based on an approximated
1174 : ! and linearized form of the theta_v equation:
1175 : !
1176 : ! theta_v = theta_l + { (R_v/R_d) - 1 } * thv_ds * r_t
1177 : ! + [ {L_v/(C_p*exner)} - (R_v/R_d) * thv_ds ] * r_c;
1178 : !
1179 : ! and therefore:
1180 : !
1181 : ! th_v' = th_l' + { (R_v/R_d) - 1 } * thv_ds * r_t'
1182 : ! + [ {L_v/(C_p*exner)} - (R_v/R_d) * thv_ds ] * r_c';
1183 : !
1184 : ! where thv_ds is used as a reference value to approximate theta_l.
1185 : !
1186 : ! Reference:
1187 : ! https://arxiv.org/pdf/1711.03675v1.pdf#nameddest=url:anl_int_buoy_terms
1188 :
1189 : ! Calculate the contributions to <w'rc'>, <w'^2 rc'>, <rt'rc'>, and
1190 : ! <thl'rc'> from the 1st PDF component.
1191 : call calc_xprcp_component( nz, ngrdcol, & ! In
1192 : wm, rtm, thlm, um, vm, rcm, & ! In
1193 : pdf_params%w_1, pdf_params%rt_1, & ! In
1194 : pdf_params%thl_1, u_1, v_1, & ! In
1195 : pdf_params%varnce_w_1, pdf_params%chi_1, & ! In
1196 : pdf_params%stdev_chi_1, pdf_params%stdev_eta_1, & ! In
1197 : pdf_params%corr_w_chi_1, pdf_params%corr_chi_eta_1, & ! In
1198 : pdf_params%crt_1, pdf_params%cthl_1, & ! In
1199 : pdf_params%rc_1, pdf_params%cloud_frac_1, iiPDF_type, & ! In
1200 : wprcp_contrib_comp_1, wp2rcp_contrib_comp_1, & ! Out
1201 : rtprcp_contrib_comp_1, thlprcp_contrib_comp_1, & ! Out
1202 17870112 : uprcp_contrib_comp_1, vprcp_contrib_comp_1 ) ! Out
1203 :
1204 : call calc_xprcp_component( nz, ngrdcol, & ! In
1205 : wm, rtm, thlm, um, vm, rcm, & ! In
1206 : pdf_params%w_2, pdf_params%rt_2, & ! In
1207 : pdf_params%thl_2, u_2, v_2, & ! In
1208 : pdf_params%varnce_w_2, pdf_params%chi_2, & ! In
1209 : pdf_params%stdev_chi_2, pdf_params%stdev_eta_2, & ! In
1210 : pdf_params%corr_w_chi_2, pdf_params%corr_chi_eta_2, & ! In
1211 : pdf_params%crt_2, pdf_params%cthl_2, & ! In
1212 : pdf_params%rc_2, pdf_params%cloud_frac_2, iiPDF_type, & ! In
1213 : wprcp_contrib_comp_2, wp2rcp_contrib_comp_2, & ! Out
1214 : rtprcp_contrib_comp_2, thlprcp_contrib_comp_2, & ! Out
1215 17870112 : uprcp_contrib_comp_2, vprcp_contrib_comp_2 ) ! Out
1216 :
1217 :
1218 : ! Calculate rc_coef, which is the coefficient on <x'rc'> in the <x'thv'> equation.
1219 : !$acc parallel loop gang vector collapse(2) default(present)
1220 1536829632 : do k = 1, nz
1221 25380961632 : do i = 1, ngrdcol
1222 :
1223 23844132000 : rc_coef(i,k) = Lv / ( exner(i,k) * Cp ) - ep2 * thv_ds(i,k)
1224 :
1225 : ! Calculate <w'rc'>, <w'^2 rc'>, <rt'rc'>, and <thl'rc'>.
1226 0 : wprcp(i,k) = pdf_params%mixt_frac(i,k) * wprcp_contrib_comp_1(i,k) &
1227 23844132000 : + ( one - pdf_params%mixt_frac(i,k) ) * wprcp_contrib_comp_2(i,k)
1228 :
1229 0 : wp2rcp(i,k) = pdf_params%mixt_frac(i,k) * wp2rcp_contrib_comp_1(i,k) &
1230 23844132000 : + ( one - pdf_params%mixt_frac(i,k) ) * wp2rcp_contrib_comp_2(i,k)
1231 :
1232 0 : rtprcp(i,k) = pdf_params%mixt_frac(i,k) * rtprcp_contrib_comp_1(i,k) &
1233 23844132000 : + ( one - pdf_params%mixt_frac(i,k) ) * rtprcp_contrib_comp_2(i,k)
1234 :
1235 0 : thlprcp(i,k) = pdf_params%mixt_frac(i,k) * thlprcp_contrib_comp_1(i,k) &
1236 23844132000 : + ( one - pdf_params%mixt_frac(i,k) ) * thlprcp_contrib_comp_2(i,k)
1237 :
1238 0 : uprcp(i,k) = pdf_params%mixt_frac(i,k) * uprcp_contrib_comp_1(i,k) &
1239 23844132000 : + ( one - pdf_params%mixt_frac(i,k) ) * uprcp_contrib_comp_2(i,k)
1240 :
1241 0 : vprcp(i,k) = pdf_params%mixt_frac(i,k) * vprcp_contrib_comp_1(i,k) &
1242 25363091520 : + ( one - pdf_params%mixt_frac(i,k) ) * vprcp_contrib_comp_2(i,k)
1243 : end do
1244 : end do
1245 : !$acc end parallel loop
1246 :
1247 : ! Calculate <w'thv'>, <w'^2 thv'>, <rt'thv'>, and <thl'thv'>.
1248 : !$acc parallel loop gang vector collapse(2) default(present)
1249 1536829632 : do k = 1, nz
1250 25380961632 : do i = 1, ngrdcol
1251 23844132000 : wpthvp(i,k) = wpthlp(i,k) + ep1 * thv_ds(i,k) * wprtp(i,k) + rc_coef(i,k) * wprcp(i,k)
1252 23844132000 : wp2thvp(i,k) = wp2thlp(i,k) + ep1 * thv_ds(i,k) * wp2rtp(i,k) + rc_coef(i,k) * wp2rcp(i,k)
1253 23844132000 : rtpthvp(i,k) = rtpthlp(i,k) + ep1 * thv_ds(i,k) * rtp2(i,k) + rc_coef(i,k) * rtprcp(i,k)
1254 25363091520 : thlpthvp(i,k)= thlp2(i,k) + ep1 * thv_ds(i,k) * rtpthlp(i,k) + rc_coef(i,k) * thlprcp(i,k)
1255 : end do
1256 : end do
1257 : !$acc end parallel loop
1258 :
1259 : ! Add the precipitation loading term in the <x'thv'> equation.
1260 : if ( l_liq_ice_loading_test ) then
1261 :
1262 : do hm_idx = 1, hydromet_dim, 1
1263 :
1264 : if ( l_mix_rat_hm(hm_idx) ) then
1265 : !$acc parallel loop gang vector collapse(2) default(present)
1266 : do k = 1, nz
1267 : do i = 1, ngrdcol
1268 : wp2thvp(i,k) = wp2thvp(i,k) - thv_ds(i,k) * wp2hmp(i,k,hm_idx)
1269 : wpthvp(i,k) = wpthvp(i,k) - thv_ds(i,k) * wphydrometp(i,k,hm_idx)
1270 : thlpthvp(i,k) = thlpthvp(i,k) - thv_ds(i,k) * thlphmp(i,k,hm_idx)
1271 : rtpthvp(i,k) = rtpthvp(i,k) - thv_ds(i,k) * rtphmp(i,k,hm_idx)
1272 : end do
1273 : end do
1274 : !$acc end parallel loop
1275 : end if
1276 :
1277 : end do
1278 :
1279 : end if
1280 :
1281 : ! Account for subplume correlation of scalar, theta_v.
1282 : ! See Eqs. A13, A8 from Larson et al. (2002) ``Small-scale...''
1283 : ! where the ``scalar'' in this paper is w.
1284 17870112 : if ( l_scalar_calc ) then
1285 :
1286 : !$acc parallel loop gang vector collapse(3) default(present)
1287 0 : do j = 1, sclr_dim
1288 0 : do k = 1, nz
1289 0 : do i = 1, ngrdcol
1290 :
1291 0 : sclrprcp(i,k,j) &
1292 0 : = pdf_params%mixt_frac(i,k) * ( ( sclr1(i,k,j) - sclrm(i,k,j) ) * pdf_params%rc_1(i,k) ) &
1293 : + ( one - pdf_params%mixt_frac(i,k) ) * ( ( sclr2(i,k,j) - sclrm(i,k,j) ) &
1294 0 : * pdf_params%rc_2(i,k) ) &
1295 0 : + pdf_params%mixt_frac(i,k) * corr_sclr_rt_1(i,k,j) * pdf_params%crt_1(i,k) &
1296 0 : * sqrt( varnce_sclr1(i,k,j) * pdf_params%varnce_rt_1(i,k) ) &
1297 0 : * pdf_params%cloud_frac_1(i,k) &
1298 0 : + ( one - pdf_params%mixt_frac(i,k) ) * corr_sclr_rt_2(i,k,j) * pdf_params%crt_2(i,k) &
1299 0 : * sqrt( varnce_sclr2(i,k,j) * pdf_params%varnce_rt_2(i,k) ) &
1300 0 : * pdf_params%cloud_frac_2(i,k) &
1301 0 : - pdf_params%mixt_frac(i,k) * corr_sclr_thl_1(i,k,j) * pdf_params%cthl_1(i,k) &
1302 0 : * sqrt( varnce_sclr1(i,k,j) * pdf_params%varnce_thl_1(i,k) ) &
1303 : * pdf_params%cloud_frac_1(i,k) &
1304 0 : - ( one - pdf_params%mixt_frac(i,k) ) * corr_sclr_thl_2(i,k,j) * pdf_params%cthl_2(i,k) &
1305 0 : * sqrt( varnce_sclr2(i,k,j) * pdf_params%varnce_thl_2(i,k) ) &
1306 0 : * pdf_params%cloud_frac_2(i,k)
1307 :
1308 : sclrpthvp(i,k,j) = sclrpthlp(i,k,j) + ep1*thv_ds(i,k)*sclrprtp(i,k,j) &
1309 0 : + rc_coef(i,k)*sclrprcp(i,k,j)
1310 :
1311 : end do
1312 : end do
1313 : end do ! i=1, sclr_dim
1314 : !$acc end parallel loop
1315 :
1316 : end if ! l_scalar_calc
1317 :
1318 :
1319 : ! Compute variance of liquid water mixing ratio.
1320 : ! This is not needed for closure. Statistical Analysis only.
1321 :
1322 : #ifndef CLUBB_CAM
1323 : ! if CLUBB is used in CAM we want this variable computed no matter what
1324 : if ( stats_metadata%ircp2 > 0 ) then
1325 : #endif
1326 : !$acc parallel loop gang vector collapse(2) default(present)
1327 1536829632 : do k = 1,nz
1328 25380961632 : do i = 1, ngrdcol
1329 47688264000 : rcp2(i,k) = pdf_params%mixt_frac(i,k) &
1330 0 : * ( pdf_params%chi_1(i,k)*pdf_params%rc_1(i,k) &
1331 0 : + pdf_params%cloud_frac_1(i,k)*pdf_params%stdev_chi_1(i,k)**2 ) &
1332 : + ( one-pdf_params%mixt_frac(i,k) ) &
1333 0 : * ( pdf_params%chi_2(i,k)*pdf_params%rc_2(i,k) &
1334 0 : + pdf_params%cloud_frac_2(i,k)*pdf_params%stdev_chi_2(i,k)**2 ) &
1335 71532396000 : - rcm(i,k)**2
1336 25363091520 : rcp2(i,k) = max( zero_threshold, rcp2(i,k) )
1337 :
1338 : end do
1339 : end do
1340 : !$acc end parallel loop
1341 : #ifndef CLUBB_CAM
1342 : ! if CLUBB is used in CAM we want this variable computed no matter what
1343 : end if
1344 : #endif
1345 :
1346 : if ( ( iiPDF_type == iiPDF_ADG1 .or. iiPDF_type == iiPDF_ADG2 &
1347 : .or. iiPDF_type == iiPDF_new_hybrid ) &
1348 17870112 : .and. ( stats_metadata%iw_up_in_cloud > 0 .or. stats_metadata%iw_down_in_cloud > 0 ) ) then
1349 :
1350 : call calc_w_up_in_cloud( nz, ngrdcol, & ! In
1351 : pdf_params%mixt_frac, & ! In
1352 : pdf_params%cloud_frac_1, pdf_params%cloud_frac_2, & ! In
1353 : pdf_params%w_1, pdf_params%w_2, & ! In
1354 : pdf_params%varnce_w_1, pdf_params%varnce_w_2, & ! In
1355 : w_up_in_cloud, w_down_in_cloud, & ! Out
1356 0 : cloudy_updraft_frac, cloudy_downdraft_frac ) ! Out
1357 :
1358 : else
1359 : !$acc parallel loop gang vector collapse(2) default(present)
1360 1536829632 : do k = 1,nz
1361 25380961632 : do i = 1, ngrdcol
1362 23844132000 : w_up_in_cloud(i,k) = zero
1363 25363091520 : w_down_in_cloud(i,k) = zero
1364 : end do
1365 : end do
1366 : end if
1367 :
1368 : #ifdef TUNER
1369 :
1370 : !$acc update host( pdf_params, pdf_params%thl_1, pdf_params%thl_2 )
1371 :
1372 : ! Check the first levels (and first gridcolumn) for reasonable temperatures
1373 : ! greater than 190K and less than 1000K
1374 : ! This is necessary because for certain parameter sets we can get floating point errors
1375 : do i=1, min( 10, size(pdf_params%thl_1(1,:)) )
1376 : if ( pdf_params%thl_1(1,i) < 190. ) then
1377 : write(fstderr,*) "Fatal error: pdf_params%thl_1 =", pdf_params%thl_1(1,i), &
1378 : " < 190K at first grid column and grid level i = ", i
1379 : err_code = clubb_fatal_error
1380 : return
1381 : end if
1382 : if ( pdf_params%thl_2(1,i) < 190. ) then
1383 : write(fstderr,*) "Fatal error: pdf_params%thl_2 =", pdf_params%thl_2(1,i), &
1384 : " < 190K at first grid column and grid level i = ", i
1385 : err_code = clubb_fatal_error
1386 : return
1387 : end if
1388 : if ( pdf_params%thl_1(1,i) > 1000. ) then
1389 : write(fstderr,*) "Fatal error: pdf_params%thl_1 =", pdf_params%thl_1(1,i), &
1390 : " > 1000K at first grid column and grid level i = ", i
1391 : err_code = clubb_fatal_error
1392 : return
1393 : end if
1394 : if ( pdf_params%thl_2(1,i) > 1000. ) then
1395 : write(fstderr,*) "Fatal error: pdf_params%thl_2 =", pdf_params%thl_2(1,i), &
1396 : " > 1000K at first grid column and grid level i = ", i
1397 : err_code = clubb_fatal_error
1398 : return
1399 : end if
1400 : end do
1401 : #endif /*TUNER*/
1402 :
1403 17870112 : if ( clubb_at_least_debug_level( 2 ) ) then
1404 :
1405 : !$acc update host( wp4, wprtp2, wp2rtp, wpthlp2, wp2thlp, cloud_frac, &
1406 : !$acc rcm, wpthvp, wp2thvp, rtpthvp, thlpthvp, wprcp, wp2rcp, &
1407 : !$acc rtprcp, thlprcp, rcp2, wprtpthlp, &
1408 : !$acc pdf_params%w_1, pdf_params%w_2, &
1409 : !$acc pdf_params%varnce_w_1, pdf_params%varnce_w_2, &
1410 : !$acc pdf_params%rt_1, pdf_params%rt_2, &
1411 : !$acc pdf_params%varnce_rt_1, pdf_params%varnce_rt_2, &
1412 : !$acc pdf_params%thl_1, pdf_params%thl_2, &
1413 : !$acc pdf_params%varnce_thl_1, pdf_params%varnce_thl_2, &
1414 : !$acc pdf_params%corr_w_rt_1, pdf_params%corr_w_rt_2, &
1415 : !$acc pdf_params%corr_w_thl_1, pdf_params%corr_w_thl_2, &
1416 : !$acc pdf_params%corr_rt_thl_1, pdf_params%corr_rt_thl_2,&
1417 : !$acc pdf_params%alpha_thl, pdf_params%alpha_rt, &
1418 : !$acc pdf_params%crt_1, pdf_params%crt_2, pdf_params%cthl_1, &
1419 : !$acc pdf_params%cthl_2, pdf_params%chi_1, &
1420 : !$acc pdf_params%chi_2, pdf_params%stdev_chi_1, &
1421 : !$acc pdf_params%stdev_chi_2, pdf_params%stdev_eta_1, &
1422 : !$acc pdf_params%stdev_eta_2, pdf_params%covar_chi_eta_1, &
1423 : !$acc pdf_params%covar_chi_eta_2, pdf_params%corr_w_chi_1, &
1424 : !$acc pdf_params%corr_w_chi_2, pdf_params%corr_w_eta_1, &
1425 : !$acc pdf_params%corr_w_eta_2, pdf_params%corr_chi_eta_1, &
1426 : !$acc pdf_params%corr_chi_eta_2, pdf_params%rsatl_1, &
1427 : !$acc pdf_params%rsatl_2, pdf_params%rc_1, pdf_params%rc_2, &
1428 : !$acc pdf_params%cloud_frac_1, pdf_params%cloud_frac_2, &
1429 : !$acc pdf_params%mixt_frac, pdf_params%ice_supersat_frac_1, &
1430 : !$acc pdf_params%ice_supersat_frac_2 )
1431 :
1432 : !$acc update host( sclrpthvp, sclrprcp, wpsclrp2, wpsclrprtp, wpsclrpthlp, wp2sclrp ) &
1433 : !$acc if ( sclr_dim > 0 )
1434 :
1435 0 : do i = 1, ngrdcol
1436 :
1437 : call pdf_closure_check( &
1438 : nz, sclr_dim, &
1439 0 : wp4(i,:), wprtp2(i,:), wp2rtp(i,:), wpthlp2(i,:), & ! intent(in)
1440 0 : wp2thlp(i,:), cloud_frac(i,:), rcm(i,:), wpthvp(i,:), wp2thvp(i,:), & ! intent(in)
1441 0 : rtpthvp(i,:), thlpthvp(i,:), wprcp(i,:), wp2rcp(i,:), & ! intent(in)
1442 0 : rtprcp(i,:), thlprcp(i,:), rcp2(i,:), wprtpthlp(i,:), & ! intent(in)
1443 0 : pdf_params%crt_1(i,:), pdf_params%crt_2(i,:), & ! intent(in)
1444 0 : pdf_params%cthl_1(i,:), pdf_params%cthl_2(i,:), & ! intent(in)
1445 : pdf_params, & ! intent(in)
1446 0 : sclrpthvp(i,:,:), sclrprcp(i,:,:), wpsclrp2(i,:,:), & ! intent(in)
1447 0 : wpsclrprtp(i,:,:), wpsclrpthlp(i,:,:), wp2sclrp(i,:,:), & ! intent(in)
1448 0 : stats_metadata ) ! intent(in)
1449 : end do
1450 : end if
1451 :
1452 : ! Error Reporting
1453 : ! Joshua Fasching February 2008
1454 17870112 : if ( clubb_at_least_debug_level( 2 ) ) then
1455 0 : if ( err_code == clubb_fatal_error ) then
1456 :
1457 : !$acc update host( p_in_Pa, exner, thv_ds, wm, wp2, wp3, sigma_sqd_w, &
1458 : !$acc rtm, rtp2, wprtp, thlm, thlp2, wpthlp, rtpthlp, ice_supersat_frac )
1459 :
1460 : !$acc update host( sclrm, wpsclrp, sclrp2, sclrprtp, sclrpthlp ) if ( sclr_dim > 0 )
1461 :
1462 0 : write(fstderr,*) "Error in pdf_closure_new"
1463 :
1464 0 : write(fstderr,*) "Intent(in)"
1465 :
1466 0 : write(fstderr,*) "p_in_Pa = ", p_in_Pa
1467 0 : write(fstderr,*) "exner = ", exner
1468 0 : write(fstderr,*) "thv_ds = ", thv_ds
1469 0 : write(fstderr,*) "wm = ", wm
1470 0 : write(fstderr,*) "wp2 = ", wp2
1471 0 : write(fstderr,*) "wp3 = ", wp3
1472 0 : write(fstderr,*) "sigma_sqd_w = ", sigma_sqd_w
1473 0 : write(fstderr,*) "rtm = ", rtm
1474 0 : write(fstderr,*) "rtp2 = ", rtp2
1475 0 : write(fstderr,*) "wprtp = ", wprtp
1476 0 : write(fstderr,*) "thlm = ", thlm
1477 0 : write(fstderr,*) "thlp2 = ", thlp2
1478 0 : write(fstderr,*) "wpthlp = ", wpthlp
1479 0 : write(fstderr,*) "rtpthlp = ", rtpthlp
1480 :
1481 0 : if ( sclr_dim > 0 ) then
1482 0 : write(fstderr,*) "sclrm = ", sclrm
1483 0 : write(fstderr,*) "wpsclrp = ", wpsclrp
1484 0 : write(fstderr,*) "sclrp2 = ", sclrp2
1485 0 : write(fstderr,*) "sclrprtp = ", sclrprtp
1486 0 : write(fstderr,*) "sclrpthlp = ", sclrpthlp
1487 : end if
1488 :
1489 0 : write(fstderr,*) "Intent(out)"
1490 :
1491 0 : write(fstderr,*) "wp4 = ", wp4
1492 0 : if ( l_explicit_turbulent_adv_xpyp .or. stats_metadata%iwprtp2 > 0 ) then
1493 0 : write(fstderr,*) "wprtp2 = ", wprtp2
1494 : end if
1495 0 : write(fstderr,*) "wp2rtp = ", wp2rtp
1496 0 : if ( l_explicit_turbulent_adv_xpyp .or. stats_metadata%iwpthlp2 > 0 ) then
1497 0 : write(fstderr,*) "wpthlp2 = ", wpthlp2
1498 : end if
1499 0 : write(fstderr,*) "cloud_frac = ", cloud_frac
1500 0 : write(fstderr,*) "ice_supersat_frac = ", ice_supersat_frac
1501 0 : write(fstderr,*) "rcm = ", rcm
1502 0 : write(fstderr,*) "wpthvp = ", wpthvp
1503 0 : write(fstderr,*) "wp2thvp = ", wp2thvp
1504 0 : write(fstderr,*) "rtpthvp = ", rtpthvp
1505 0 : write(fstderr,*) "thlpthvp = ", thlpthvp
1506 0 : write(fstderr,*) "wprcp = ", wprcp
1507 0 : write(fstderr,*) "wp2rcp = ", wp2rcp
1508 0 : write(fstderr,*) "rtprcp = ", rtprcp
1509 0 : write(fstderr,*) "thlprcp = ", thlprcp
1510 : #ifndef CLUBB_CAM
1511 : ! if CLUBB is used in CAM we want this variable computed no matter what
1512 : if ( stats_metadata%ircp2 > 0 ) then
1513 : #endif
1514 0 : write(fstderr,*) "rcp2 = ", rcp2
1515 : #ifndef CLUBB_CAM
1516 : end if
1517 : #endif
1518 0 : if ( l_explicit_turbulent_adv_xpyp .or. stats_metadata%iwprtpthlp > 0 ) then
1519 0 : write(fstderr,*) "wprtpthlp = ", wprtpthlp
1520 : end if
1521 0 : write(fstderr,*) "rcp2 = ", rcp2
1522 0 : write(fstderr,*) "wprtpthlp = ", wprtpthlp
1523 0 : write(fstderr,*) "pdf_params%w_1 = ", pdf_params%w_1
1524 0 : write(fstderr,*) "pdf_params%w_2 = ", pdf_params%w_2
1525 0 : write(fstderr,*) "pdf_params%varnce_w_1 = ", pdf_params%varnce_w_1
1526 0 : write(fstderr,*) "pdf_params%varnce_w_2 = ", pdf_params%varnce_w_2
1527 0 : write(fstderr,*) "pdf_params%rt_1 = ", pdf_params%rt_1
1528 0 : write(fstderr,*) "pdf_params%rt_2 = ", pdf_params%rt_2
1529 0 : write(fstderr,*) "pdf_params%varnce_rt_1 = ", pdf_params%varnce_rt_1
1530 0 : write(fstderr,*) "pdf_params%varnce_rt_2 = ", pdf_params%varnce_rt_2
1531 0 : write(fstderr,*) "pdf_params%thl_1 = ", pdf_params%thl_1
1532 0 : write(fstderr,*) "pdf_params%thl_2 = ", pdf_params%thl_2
1533 0 : write(fstderr,*) "pdf_params%varnce_thl_1 = ", pdf_params%varnce_thl_1
1534 0 : write(fstderr,*) "pdf_params%varnce_thl_2 = ", pdf_params%varnce_thl_2
1535 0 : write(fstderr,*) "pdf_params%corr_w_rt_1 = ", pdf_params%corr_w_rt_1
1536 0 : write(fstderr,*) "pdf_params%corr_w_rt_2 = ", pdf_params%corr_w_rt_2
1537 0 : write(fstderr,*) "pdf_params%corr_w_thl_1 = ", pdf_params%corr_w_thl_1
1538 0 : write(fstderr,*) "pdf_params%corr_w_thl_2 = ", pdf_params%corr_w_thl_2
1539 0 : write(fstderr,*) "pdf_params%corr_rt_thl_1 = ", pdf_params%corr_rt_thl_1
1540 0 : write(fstderr,*) "pdf_params%corr_rt_thl_2 = ", pdf_params%corr_rt_thl_2
1541 0 : write(fstderr,*) "pdf_params%alpha_thl = ", pdf_params%alpha_thl
1542 0 : write(fstderr,*) "pdf_params%alpha_rt = ", pdf_params%alpha_rt
1543 0 : write(fstderr,*) "pdf_params%crt_1 = ", pdf_params%crt_1
1544 0 : write(fstderr,*) "pdf_params%crt_2 = ", pdf_params%crt_2
1545 0 : write(fstderr,*) "pdf_params%cthl_1 = ", pdf_params%cthl_1
1546 0 : write(fstderr,*) "pdf_params%cthl_2 = ", pdf_params%cthl_2
1547 0 : write(fstderr,*) "pdf_params%chi_1 = ", pdf_params%chi_1
1548 0 : write(fstderr,*) "pdf_params%chi_2 = ", pdf_params%chi_2
1549 0 : write(fstderr,*) "pdf_params%stdev_chi_1 = ", pdf_params%stdev_chi_1
1550 0 : write(fstderr,*) "pdf_params%stdev_chi_2 = ", pdf_params%stdev_chi_2
1551 0 : write(fstderr,*) "pdf_params%stdev_eta_1 = ", pdf_params%stdev_eta_1
1552 0 : write(fstderr,*) "pdf_params%stdev_eta_2 = ", pdf_params%stdev_eta_2
1553 0 : write(fstderr,*) "pdf_params%covar_chi_eta_1 = ", &
1554 0 : pdf_params%covar_chi_eta_1
1555 0 : write(fstderr,*) "pdf_params%covar_chi_eta_2 = ", &
1556 0 : pdf_params%covar_chi_eta_2
1557 0 : write(fstderr,*) "pdf_params%corr_w_chi_1 = ", pdf_params%corr_w_chi_1
1558 0 : write(fstderr,*) "pdf_params%corr_w_chi_2 = ", pdf_params%corr_w_chi_2
1559 0 : write(fstderr,*) "pdf_params%corr_w_eta_1 = ", pdf_params%corr_w_eta_1
1560 0 : write(fstderr,*) "pdf_params%corr_w_eta_2 = ", pdf_params%corr_w_eta_2
1561 0 : write(fstderr,*) "pdf_params%corr_chi_eta_1 = ", &
1562 0 : pdf_params%corr_chi_eta_1
1563 0 : write(fstderr,*) "pdf_params%corr_chi_eta_2 = ", &
1564 0 : pdf_params%corr_chi_eta_2
1565 0 : write(fstderr,*) "pdf_params%rsatl_1 = ", pdf_params%rsatl_1
1566 0 : write(fstderr,*) "pdf_params%rsatl_2 = ", pdf_params%rsatl_2
1567 0 : write(fstderr,*) "pdf_params%rc_1 = ", pdf_params%rc_1
1568 0 : write(fstderr,*) "pdf_params%rc_2 = ", pdf_params%rc_2
1569 0 : write(fstderr,*) "pdf_params%cloud_frac_1 = ", pdf_params%cloud_frac_1
1570 0 : write(fstderr,*) "pdf_params%cloud_frac_2 = ", pdf_params%cloud_frac_2
1571 0 : write(fstderr,*) "pdf_params%mixt_frac = ", pdf_params%mixt_frac
1572 0 : write(fstderr,*) "pdf_params%ice_supersat_frac_1 = ", &
1573 0 : pdf_params%ice_supersat_frac_1
1574 0 : write(fstderr,*) "pdf_params%ice_supersat_frac_2 = ", &
1575 0 : pdf_params%ice_supersat_frac_2
1576 :
1577 0 : if ( sclr_dim > 0 )then
1578 0 : write(fstderr,*) "sclrpthvp = ", sclrpthvp
1579 0 : write(fstderr,*) "sclrprcp = ", sclrprcp
1580 0 : write(fstderr,*) "wpsclrp2 = ", wpsclrp2
1581 0 : write(fstderr,*) "wpsclrprtp = ", wpsclrprtp
1582 0 : write(fstderr,*) "wpsclrpthlp = ", wpsclrpthlp
1583 0 : write(fstderr,*) "wp2sclrp = ", wp2sclrp
1584 : end if
1585 :
1586 0 : return
1587 :
1588 : end if ! Fatal error
1589 :
1590 0 : do i = 1, ngrdcol
1591 :
1592 : ! Error check pdf parameters and moments to ensure consistency
1593 0 : if ( iiPDF_type == iiPDF_3D_Luhar ) then
1594 :
1595 : ! Means
1596 0 : wm_clubb_pdf(i,:) = pdf_params%mixt_frac(i,:) * pdf_params%w_1(i,:) &
1597 0 : + ( one - pdf_params%mixt_frac(i,:) ) * pdf_params%w_2(i,:)
1598 :
1599 0 : do k = 1, nz, 1
1600 0 : if ( abs( ( wm_clubb_pdf(i,k) - wm(i,k) ) &
1601 0 : / max( wm(i,k), eps ) ) > .05_core_rknd ) then
1602 0 : write(fstderr,*) "wm error at thlm = ", thlm(i,k), &
1603 : ( ( wm_clubb_pdf(i,k) - wm(i,k) ) &
1604 0 : / max( wm(i,k), eps ) )
1605 : end if
1606 : end do ! k = 1, nz, 1
1607 :
1608 0 : rtm_clubb_pdf(i,:) = pdf_params%mixt_frac(i,:) * pdf_params%rt_1(i,:) &
1609 0 : + ( one - pdf_params%mixt_frac(i,:) ) * pdf_params%rt_2(i,:)
1610 :
1611 0 : do k = 1, nz, 1
1612 0 : if ( abs( ( rtm_clubb_pdf(i,k) - rtm(i,k) ) &
1613 0 : / max( rtm(i,k), eps ) ) > .05_core_rknd ) then
1614 0 : write(fstderr,*) "rtm error at thlm = ", thlm(i,k), &
1615 : ( ( rtm_clubb_pdf(i,k) - rtm(i,k) ) &
1616 0 : / max( rtm(i,k), eps ) )
1617 : end if
1618 : end do ! k = 1, nz, 1
1619 :
1620 0 : thlm_clubb_pdf(i,:) = pdf_params%mixt_frac(i,:) * pdf_params%thl_1(i,:) &
1621 0 : + ( one - pdf_params%mixt_frac(i,:) ) * pdf_params%thl_2(i,:)
1622 :
1623 0 : do k = 1, nz, 1
1624 0 : if ( abs( ( thlm_clubb_pdf(i,k) - thlm(i,k) ) / thlm(i,k) ) &
1625 0 : > .05_core_rknd ) then
1626 0 : write(fstderr,*) "thlm error at thlm = ", thlm(i,k), &
1627 0 : ( ( thlm_clubb_pdf(i,k) - thlm(i,k) ) / thlm(i,k) )
1628 : end if
1629 : end do ! k = 1, nz, 1
1630 :
1631 : ! Variances
1632 0 : wp2_clubb_pdf(i,:) = pdf_params%mixt_frac(i,:) &
1633 0 : * ( ( pdf_params%w_1(i,:) - wm(i,:) )**2 + pdf_params%varnce_w_1(i,:) ) &
1634 : + ( one - pdf_params%mixt_frac(i,:) ) &
1635 0 : * ( ( pdf_params%w_2(i,:) - wm(i,:) )**2 + pdf_params%varnce_w_2(i,:) )
1636 :
1637 0 : do k = 1, nz, 1
1638 0 : if ( wp2(i,k) > w_tol**2 ) then
1639 0 : if ( abs( ( wp2_clubb_pdf(i,k) - wp2(i,k) ) / wp2(i,k) ) &
1640 : > .05_core_rknd ) then
1641 0 : write(fstderr,*) "wp2 error at thlm = ", thlm(i,k), &
1642 0 : ( ( wp2_clubb_pdf(i,k) - wp2(i,k) ) / wp2(i,k) )
1643 : end if
1644 : end if
1645 : end do ! k = 1, nz, 1
1646 :
1647 : rtp2_clubb_pdf(i,:) &
1648 0 : = pdf_params%mixt_frac(i,:) &
1649 0 : * ( ( pdf_params%rt_1(i,:) - rtm(i,:) )**2 + pdf_params%varnce_rt_1(i,:) ) &
1650 : + ( one - pdf_params%mixt_frac(i,:) ) &
1651 0 : * ( ( pdf_params%rt_2(i,:) - rtm(i,:) )**2 + pdf_params%varnce_rt_2(i,:) )
1652 :
1653 0 : do k = 1, nz, 1
1654 0 : if ( rtp2(i,k) > rt_tol**2 ) then
1655 0 : if ( abs( ( rtp2_clubb_pdf(i,k) - rtp2(i,k) ) / rtp2(i,k) ) &
1656 : > .05_core_rknd ) then
1657 0 : write(fstderr,*) "rtp2 error at thlm = ", thlm(i,k), &
1658 0 : "Error = ", ( ( rtp2_clubb_pdf(i,k) - rtp2(i,k) ) / rtp2(i,k) )
1659 : end if
1660 : end if
1661 : end do ! k = 1, nz, 1
1662 :
1663 : thlp2_clubb_pdf(i,:) &
1664 0 : = pdf_params%mixt_frac(i,:) &
1665 0 : * ( ( pdf_params%thl_1(i,:) - thlm(i,:) )**2 + pdf_params%varnce_thl_1(i,:) ) &
1666 : + ( one - pdf_params%mixt_frac(i,:) ) &
1667 0 : * ( ( pdf_params%thl_2(i,:) - thlm(i,:) )**2 + pdf_params%varnce_thl_2(i,:) )
1668 :
1669 0 : do k = 1, nz, 1
1670 0 : if( thlp2(i,k) > thl_tol**2 ) then
1671 0 : if ( abs( ( thlp2_clubb_pdf(i,k) - thlp2(i,k) ) / thlp2(i,k) ) &
1672 : > .05_core_rknd ) then
1673 0 : write(fstderr,*) "thlp2 error at thlm = ", thlm(i,k), &
1674 0 : "Error = ", ( ( thlp2_clubb_pdf(i,k) - thlp2(i,k) ) / thlp2(i,k) )
1675 : end if
1676 : end if
1677 : end do ! k = 1, nz, 1
1678 :
1679 : ! Third order moments
1680 : wp3_clubb_pdf(i,:) &
1681 0 : = pdf_params%mixt_frac(i,:) * ( pdf_params%w_1(i,:) - wm(i,:) ) &
1682 0 : * ( ( pdf_params%w_1(i,:) - wm(i,:) )**2 + three * pdf_params%varnce_w_1(i,:) ) &
1683 0 : + ( one - pdf_params%mixt_frac(i,:) ) * ( pdf_params%w_2(i,:) - wm(i,:) ) &
1684 0 : * ( ( pdf_params%w_2(i,:) - wm(i,:) )**2 + three * pdf_params%varnce_w_2(i,:) )
1685 :
1686 : rtp3_clubb_pdf(i,:) &
1687 0 : = pdf_params%mixt_frac(i,:) * ( pdf_params%rt_1(i,:) - rtm(i,:) ) &
1688 0 : * ( ( pdf_params%rt_1(i,:) - rtm(i,:) )**2 + three * pdf_params%varnce_rt_1(i,:) ) &
1689 0 : + ( one - pdf_params%mixt_frac(i,:) ) * ( pdf_params%rt_2(i,:) - rtm(i,:) ) &
1690 0 : * ( ( pdf_params%rt_2(i,:) - rtm(i,:) )**2 + three * pdf_params%varnce_rt_2(i,:) )
1691 :
1692 : thlp3_clubb_pdf(i,:) &
1693 0 : = pdf_params%mixt_frac(i,:) * ( pdf_params%thl_1(i,:) - thlm(i,:) ) &
1694 0 : * ( ( pdf_params%thl_1(i,:) - thlm(i,:) )**2 + three * pdf_params%varnce_thl_1(i,:) ) &
1695 0 : + ( one - pdf_params%mixt_frac(i,:) ) * ( pdf_params%thl_2(i,:) - thlm(i,:) ) &
1696 0 : * ( ( pdf_params%thl_2(i,:) - thlm(i,:) )**2 + three * pdf_params%varnce_thl_2(i,:) )
1697 :
1698 : ! Skewness
1699 0 : Skw_denom_coef = clubb_params(i,iSkw_denom_coef)
1700 :
1701 : Skw_clubb_pdf(i,:) &
1702 : = wp3_clubb_pdf(i,:) &
1703 0 : / ( wp2_clubb_pdf(i,:) + Skw_denom_coef * w_tol**2 )**1.5_core_rknd
1704 :
1705 0 : do k = 1, nz, 1
1706 0 : if ( Skw(i,k) > .05_core_rknd ) then
1707 0 : if( abs( ( Skw_clubb_pdf(i,k) - Skw(i,k) ) / Skw(i,k) ) &
1708 : > .25_core_rknd ) then
1709 0 : write(fstderr,*) "Skw error at thlm = ", thlm(i,k), &
1710 0 : "Error = ", ( ( Skw_clubb_pdf(i,k) - Skw(i,k) ) / Skw(i,k) ), &
1711 0 : Skw_clubb_pdf(i,k), Skw(i,k)
1712 : end if
1713 : end if
1714 : end do ! k = 1, nz, 1
1715 :
1716 : Skrt_clubb_pdf(i,:) &
1717 : = rtp3_clubb_pdf(i,:) &
1718 0 : / ( rtp2_clubb_pdf(i,:) + Skw_denom_coef * rt_tol**2 )**1.5_core_rknd
1719 :
1720 0 : do k = 1, nz, 1
1721 0 : if ( Skrt(i,k) > .05_core_rknd ) then
1722 0 : if( abs( ( Skrt_clubb_pdf(i,k) - Skrt(i,k) ) / Skrt(i,k) ) &
1723 : > .25_core_rknd ) then
1724 0 : write(fstderr,*) "Skrt error at thlm = ", thlm(i,k), &
1725 0 : "Error = ", ( ( Skrt_clubb_pdf(i,k) - Skrt(i,k) ) / Skrt(i,k) ), &
1726 0 : Skrt_clubb_pdf(i,k), Skrt(i,k)
1727 : end if
1728 : end if
1729 : end do ! k = 1, nz, 1
1730 :
1731 : Skthl_clubb_pdf(i,:) &
1732 : = thlp3_clubb_pdf(i,:) &
1733 0 : / ( thlp2_clubb_pdf(i,:) + Skw_denom_coef * thl_tol**2 )**1.5_core_rknd
1734 :
1735 0 : do k = 1, nz, 1
1736 0 : if ( Skthl(i,k) > .05_core_rknd ) then
1737 0 : if ( abs( ( Skthl_clubb_pdf(i,k) - Skthl(i,k) ) / Skthl(i,k) ) &
1738 : > .25_core_rknd ) then
1739 0 : write(fstderr,*) "Skthl error at thlm = ", thlm(i,k), &
1740 0 : "Error = ", ( ( Skthl_clubb_pdf(i,k) - Skthl(i,k) ) / Skthl(i,k) ), &
1741 0 : Skthl_clubb_pdf(i,k), Skthl(i,k)
1742 : end if
1743 : end if
1744 : end do ! k = 1, nz, 1
1745 :
1746 : end if ! iiPDF_type == iiPDF_3D_Luhar
1747 :
1748 : end do
1749 :
1750 : end if ! clubb_at_least_debug_level
1751 :
1752 : !$acc exit data delete( u_1, u_2, varnce_u_1, varnce_u_2, v_1, v_2, &
1753 : !$acc varnce_v_1, varnce_v_2, alpha_u, alpha_v, &
1754 : !$acc corr_u_w_1, corr_u_w_2, corr_v_w_1, corr_v_w_2, &
1755 : !$acc tl1, tl2, sqrt_wp2, Skthl, &
1756 : !$acc Skrt, Sku, Skv, wprcp_contrib_comp_1, wprcp_contrib_comp_2, &
1757 : !$acc wp2rcp_contrib_comp_1, wp2rcp_contrib_comp_2, &
1758 : !$acc rtprcp_contrib_comp_1, rtprcp_contrib_comp_2, &
1759 : !$acc thlprcp_contrib_comp_1, thlprcp_contrib_comp_2, &
1760 : !$acc uprcp_contrib_comp_1, uprcp_contrib_comp_2, &
1761 : !$acc vprcp_contrib_comp_1, vprcp_contrib_comp_2, &
1762 : !$acc rc_1_ice, rc_2_ice, rsatl_1, rsatl_2 )
1763 :
1764 : !$acc exit data if( sclr_dim > 0 ) &
1765 : !$acc delete( sclr1, sclr2, varnce_sclr1, varnce_sclr2, &
1766 : !$acc alpha_sclr, corr_sclr_thl_1, corr_sclr_thl_2, &
1767 : !$acc corr_sclr_rt_1, corr_sclr_rt_2, corr_w_sclr_1, &
1768 : !$acc corr_w_sclr_2, Sksclr )
1769 :
1770 : return
1771 :
1772 : end subroutine pdf_closure
1773 :
1774 : !===============================================================================================
1775 35740224 : subroutine transform_pdf_chi_eta_component( nz, ngrdcol, &
1776 35740224 : tl, rsatl, rt, exner, & ! In
1777 35740224 : varnce_thl, varnce_rt, & ! In
1778 35740224 : corr_rt_thl, chi, & ! In
1779 35740224 : crt, cthl, & ! Out
1780 35740224 : stdev_chi, stdev_eta, & ! Out
1781 35740224 : covar_chi_eta, & ! Out
1782 35740224 : corr_chi_eta ) ! Out
1783 :
1784 : use clubb_precision, only: &
1785 : core_rknd ! Variable(s)
1786 :
1787 : use constants_clubb, only: &
1788 : zero, one, two, &
1789 : ep, Lv, Rd, Cp, &
1790 : chi_tol, &
1791 : eta_tol, &
1792 : max_mag_correlation
1793 :
1794 : implicit none
1795 :
1796 : integer, intent(in) :: &
1797 : ngrdcol, & ! Number of grid columns
1798 : nz ! Number of vertical level
1799 :
1800 : ! ----------- Input Variables -----------
1801 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
1802 : tl, &
1803 : rsatl, &
1804 : rt, &
1805 : varnce_thl, &
1806 : varnce_rt, &
1807 : corr_rt_thl, &
1808 : exner
1809 :
1810 : ! ----------- Output Variables -----------
1811 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
1812 : chi, & ! s from Lewellen and Yoh 1993 (LY) eqn. 1
1813 : crt, & ! Coefficients for s'
1814 : cthl, & ! Coefficients for s'
1815 : stdev_chi, & ! Standard deviation of chi for each component.
1816 : stdev_eta, & ! Standard deviation of eta for each component.
1817 : covar_chi_eta, & ! Covariance of chi and eta for each component.
1818 : corr_chi_eta ! Correlation of chi and eta for each component.
1819 :
1820 : ! ----------- Local Variables -----------
1821 : real( kind = core_rknd ) :: &
1822 : varnce_rt_term, &
1823 : corr_rt_thl_term, &
1824 : varnce_thl_term, &
1825 : varnce_chi, &
1826 : varnce_eta, &
1827 : beta, &
1828 : invrs_beta_rsatl_p1
1829 :
1830 : real( kind = core_rknd ), parameter :: &
1831 : chi_tol_sqd = chi_tol**2, &
1832 : eta_tol_sqd = eta_tol**2, &
1833 : Cp_on_Lv = Cp / Lv
1834 :
1835 : ! Loop variable
1836 : integer :: k, i
1837 :
1838 : ! ----------- Begin Code -----------
1839 :
1840 : !$acc parallel loop gang vector collapse(2) default(present)
1841 3073659264 : do k = 1, nz
1842 50761923264 : do i = 1, ngrdcol
1843 :
1844 : ! SD's beta (eqn. 8)
1845 47688264000 : beta = ep * Lv**2 / ( Rd * Cp * tl(i,k)**2 )
1846 :
1847 47688264000 : invrs_beta_rsatl_p1 = one / ( one + beta * rsatl(i,k) )
1848 :
1849 : ! s from Lewellen and Yoh 1993 (LY) eqn. 1
1850 47688264000 : chi(i,k) = ( rt(i,k) - rsatl(i,k) ) * invrs_beta_rsatl_p1
1851 :
1852 : ! For each normal distribution in the sum of two normal distributions,
1853 : ! s' = crt * rt' + cthl * thl';
1854 : ! therefore, x's' = crt * x'rt' + cthl * x'thl'.
1855 : ! Larson et al. May, 2001.
1856 47688264000 : crt(i,k) = invrs_beta_rsatl_p1
1857 : cthl(i,k) = ( one + beta * rt(i,k) ) * invrs_beta_rsatl_p1**2 &
1858 50726183040 : * Cp_on_Lv * beta * rsatl(i,k) * exner(i,k)
1859 :
1860 : end do
1861 : end do
1862 : !$acc end parallel loop
1863 :
1864 : ! Calculate covariance, correlation, and standard deviation of
1865 : ! chi and eta for each component
1866 : ! Include subplume correlation of qt, thl
1867 : !$acc parallel loop gang vector collapse(2) default(present)
1868 3073659264 : do k = 1, nz
1869 50761923264 : do i = 1, ngrdcol
1870 :
1871 47688264000 : varnce_rt_term = crt(i,k)**2 * varnce_rt(i,k)
1872 47688264000 : varnce_thl_term = cthl(i,k)**2 * varnce_thl(i,k)
1873 :
1874 47688264000 : covar_chi_eta(i,k) = varnce_rt_term - varnce_thl_term
1875 :
1876 : corr_rt_thl_term = two * corr_rt_thl(i,k) * crt(i,k) * cthl(i,k) &
1877 47688264000 : * sqrt( varnce_rt(i,k) * varnce_thl(i,k) )
1878 :
1879 47688264000 : varnce_chi = varnce_rt_term - corr_rt_thl_term + varnce_thl_term
1880 47688264000 : varnce_eta = varnce_rt_term + corr_rt_thl_term + varnce_thl_term
1881 :
1882 : ! We need to introduce a threshold value for the variance of chi and eta
1883 50726183040 : if ( varnce_chi < chi_tol_sqd .or. varnce_eta < eta_tol_sqd ) then
1884 :
1885 1164135890 : if ( varnce_chi < chi_tol_sqd ) then
1886 1163815327 : stdev_chi(i,k) = zero ! Treat chi as a delta function
1887 : else
1888 320563 : stdev_chi(i,k) = sqrt( varnce_chi )
1889 : end if
1890 :
1891 1164135890 : if ( varnce_eta < eta_tol_sqd ) then
1892 1163965456 : stdev_eta(i,k) = zero ! Treat eta as a delta function
1893 : else
1894 170434 : stdev_eta(i,k) = sqrt( varnce_eta )
1895 : end if
1896 :
1897 1164135890 : corr_chi_eta(i,k) = zero
1898 :
1899 : else
1900 :
1901 46524128110 : stdev_chi(i,k) = sqrt( varnce_chi )
1902 46524128110 : stdev_eta(i,k) = sqrt( varnce_eta )
1903 :
1904 46524128110 : corr_chi_eta(i,k) = covar_chi_eta(i,k) / ( stdev_chi(i,k) * stdev_eta(i,k) )
1905 : corr_chi_eta(i,k) = min( max_mag_correlation, &
1906 46524128110 : max( -max_mag_correlation, corr_chi_eta(i,k) ) )
1907 :
1908 : end if
1909 :
1910 : end do
1911 : end do
1912 : !$acc end parallel loop
1913 :
1914 35740224 : end subroutine transform_pdf_chi_eta_component
1915 :
1916 : !=============================================================================
1917 17870112 : subroutine calc_wp4_pdf( nz, ngrdcol, &
1918 17870112 : wm, w_1, w_2, &
1919 17870112 : varnce_w_1, varnce_w_2, &
1920 17870112 : mixt_frac, &
1921 17870112 : wp4 )
1922 :
1923 : ! Description:
1924 : ! Calculates <w'^4> by integrating over the PDF of w. The integral is:
1925 : !
1926 : ! <w'^4> = INT(-inf:inf) ( w - <w> )^4 P(w) dw;
1927 : !
1928 : ! where <w> is the overall mean of w and P(w) is a two-component normal
1929 : ! distribution of w. The integrated equation is:
1930 : !
1931 : ! <w'^4> = mixt_frac * ( 3 * sigma_w_1^4
1932 : ! + 6 * ( mu_w_1 - <w> )^2 * sigma_w_1^2
1933 : ! + ( mu_w_1 - <w> )^4 )
1934 : ! + ( 1 - mixt_frac ) * ( 3 * sigma_w_2^4
1935 : ! + 6 * ( mu_w_2 - <w> )^2 * sigma_w_2^2
1936 : ! + ( mu_w_2 - <w> )^4 );
1937 : !
1938 : ! where mu_w_1 is the mean of w in the 1st PDF component, mu_w_2 is the mean
1939 : ! of w in the 2nd PDF component, sigma_w_1 is the standard deviation of w in
1940 : ! the 1st PDF component, sigma_w_2 is the standard deviation of w in the 2nd
1941 : ! PDF component, and mixt_frac is the mixture fraction, which is the weight
1942 : ! of the 1st PDF component.
1943 :
1944 : ! References:
1945 : !-----------------------------------------------------------------------
1946 :
1947 : use constants_clubb, only: &
1948 : six, & ! Variable(s)
1949 : three, &
1950 : one
1951 :
1952 : use clubb_precision, only: &
1953 : core_rknd ! Variable(s)
1954 :
1955 : implicit none
1956 :
1957 : integer, intent(in) :: &
1958 : ngrdcol, & ! Number of grid columns
1959 : nz ! Number of vertical level
1960 :
1961 : ! Input Variables
1962 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
1963 : wm, & ! Mean of w (overall) [m/s]
1964 : w_1, & ! Mean of w (1st PDF component) [m/s]
1965 : w_2, & ! Mean of w (2nd PDF component) [m/s]
1966 : varnce_w_1, & ! Variance of w (1st PDF component) [m^2/s^2]
1967 : varnce_w_2, & ! Variance of w (2nd PDF component) [m^2/s^2]
1968 : mixt_frac ! Mixture fraction [-]
1969 :
1970 : ! Output Variable
1971 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
1972 : wp4 ! <w'^4> [m^4/s^4]
1973 :
1974 : ! Local Variables
1975 : integer :: i, k
1976 :
1977 : !$acc parallel loop gang vector collapse(2) default(present)
1978 1536829632 : do k = 1, nz
1979 25380961632 : do i = 1, ngrdcol
1980 :
1981 : ! Calculate <w'^4> by integrating over the PDF.
1982 47688264000 : wp4(i,k) = mixt_frac(i,k) * ( three * varnce_w_1(i,k)**2 &
1983 : + six * ( ( w_1(i,k) - wm(i,k) )**2 ) * varnce_w_1(i,k) &
1984 : + ( w_1(i,k) - wm(i,k) )**4 ) &
1985 : + ( one - mixt_frac(i,k) ) * ( three * varnce_w_2(i,k)**2 &
1986 : + six * ( (w_2(i,k) - wm(i,k) )**2 )*varnce_w_2(i,k) &
1987 73051355520 : + ( w_2(i,k) - wm(i,k) )**4 )
1988 : end do
1989 : end do
1990 : !$acc end parallel loop
1991 :
1992 17870112 : return
1993 :
1994 : end subroutine calc_wp4_pdf
1995 :
1996 : !=============================================================================
1997 35740224 : subroutine calc_wp2xp2_pdf( nz, ngrdcol, &
1998 35740224 : wm, xm, w_1, &
1999 35740224 : w_2, x_1, x_2, &
2000 35740224 : varnce_w_1, varnce_w_2, &
2001 35740224 : varnce_x_1, varnce_x_2, &
2002 35740224 : corr_w_x_1, corr_w_x_2, &
2003 35740224 : mixt_frac, &
2004 35740224 : wp2xp2 )
2005 :
2006 : ! Description:
2007 : ! Calculates <w'^2x'^2> by integrating over the PDF of w and x. The
2008 : ! integral
2009 : ! is:
2010 : !
2011 : ! <w'^2x'^2>
2012 : ! = INT(-inf:inf) INT(-inf:inf) ( w - <w> )^2 ( x - <x> )^2 P(w,x) dx dw;
2013 : !
2014 : ! where <w> is the overall mean of w, <x> is the overall mean of x, and
2015 : ! P(w,x) is a two-component bivariate normal distribution of w and x. The
2016 : ! integrated equation is:
2017 : !
2018 : ! <w'^2x'^2>
2019 : ! = mixt_frac
2020 : ! * ( ( mu_w_1 - <w> )**2 * ( ( mu_x_1 - <x> )**2 + sigma_x_1^2 )
2021 : ! + four * corr_w_x_1 * sigma_w_1 * sigma_x_1 * ( mu_x_1 - <x> ) * (
2022 : ! mu_w_1 - <w> )
2023 : ! + ( ( mu_x_1 - <x> )**2 + ( 1 + 2*corr_w_x_1**2 ) * sigma_x_1^2 ) *
2024 : ! sigma_w_1^2 )
2025 : ! + ( one - mixt_frac )
2026 : ! * ( ( mu_w_2 - <w> )**2 * ( ( mu_x_2 - <x> )**2 + sigma_x_2^2 )
2027 : ! + four * corr_w_x_2 * sigma_w_2 * sigma_x_2 * ( mu_x_2 - <x> ) * (
2028 : ! mu_w_2 - <w> )
2029 : ! + ( ( mu_x_2 - <x> )**2 + ( 1 + 2*corr_w_x_2**2 ) * sigma_x_2^2 ) *
2030 : ! sigma_w_2^2 )
2031 : !
2032 : ! where mu_w_1 is the mean of w in the 1st PDF component, mu_w_2 is the mean
2033 : ! of w in the 2nd PDF component, mu_x_1 is the mean of x in the 1st PDF
2034 : ! component, mu_x_2 is the mean of x in the 2nd PDF component, sigma_w_1 is
2035 : ! the standard deviation of w in the 1st PDF component, sigma_w_2 is the
2036 : ! standard deviation of w in the 2nd PDF component, sigma_x_1 is the
2037 : ! standard deviation of x in the 1st PDF component, sigma_x_2 is the
2038 : ! standard deviation of x in the 2nd PDF component, corr_w_x_1 is the
2039 : ! correlation of w and x in the 1st PDF component, corr_w_x_2 is the
2040 : ! correlation of w and x in the 2nd PDF component, and mixt_frac is the
2041 : ! mixture fraction, which is the weight of the 1st PDF component.
2042 :
2043 : ! References:
2044 : !-----------------------------------------------------------------------
2045 :
2046 : use constants_clubb, only: &
2047 : one, & ! Variable(s)
2048 : four
2049 :
2050 : use clubb_precision, only: &
2051 : core_rknd ! Variable(s)
2052 :
2053 : implicit none
2054 :
2055 : integer, intent(in) :: &
2056 : ngrdcol, & ! Number of grid columns
2057 : nz ! Number of vertical level
2058 :
2059 : ! Input Variables
2060 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
2061 : wm, & ! Mean of w (overall) [m/s]
2062 : xm, & ! Mean of x (overall) [units vary]
2063 : w_1, & ! Mean of w (1st PDF component) [m/s]
2064 : w_2, & ! Mean of w (2nd PDF component) [m/s]
2065 : x_1, & ! Mean of x (1st PDF component) [units vary]
2066 : x_2, & ! Mean of x (2nd PDF component) [units vary]
2067 : varnce_w_1, & ! Variance of w (1st PDF component) [m^2/s^2]
2068 : varnce_w_2, & ! Variance of w (2nd PDF component) [m^2/s^2]
2069 : varnce_x_1, & ! Variance of x (1st PDF component) [(units vary)^2]
2070 : varnce_x_2, & ! Variance of x (2nd PDF component) [(units vary)^2]
2071 : corr_w_x_1, & ! Correlation of w and x (1st PDF comp.) [-]
2072 : corr_w_x_2, & ! Correlation of w and x (2nd PDF comp.) [-]
2073 : mixt_frac ! Mixture fraction [-]
2074 :
2075 : ! Output Variable
2076 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
2077 : wp2xp2 ! <w'^2x'^2> [m^2/s^2 (units vary)^2]
2078 :
2079 : ! Local Variable
2080 : integer :: i, k
2081 :
2082 : !$acc parallel loop gang vector collapse(2) default(present)
2083 3073659264 : do k = 1, nz
2084 50761923264 : do i = 1, ngrdcol
2085 :
2086 : ! Calculate <w'x'^2> by integrating over the PDF.
2087 95376528000 : wp2xp2(i,k) = mixt_frac(i,k) &
2088 : * ( ( w_1(i,k) - wm(i,k) )**2 * ( ( x_1(i,k) - xm(i,k) )**2 + varnce_x_1(i,k) ) &
2089 : + four * corr_w_x_1(i,k) * sqrt( varnce_w_1(i,k) * varnce_x_1(i,k) ) &
2090 : * ( x_1(i,k) - xm(i,k) ) * ( w_1(i,k) - wm(i,k) ) &
2091 : + ( ( x_1(i,k) - xm(i,k) )**2 &
2092 : + ( 1 + 2*corr_w_x_1(i,k)**2 ) * varnce_x_1(i,k) ) * varnce_w_1(i,k) ) &
2093 : + ( one - mixt_frac(i,k) ) &
2094 : * ( ( w_2(i,k) - wm(i,k) )**2 * ( ( x_2(i,k) - xm(i,k) )**2 + varnce_x_2(i,k) ) &
2095 : + four * corr_w_x_2(i,k) * sqrt( varnce_w_2(i,k) * varnce_x_2(i,k) ) &
2096 : * ( x_2(i,k) - xm(i,k) ) * ( w_2(i,k) - wm(i,k) ) &
2097 : + ( ( x_2(i,k) - xm(i,k) )**2 &
2098 >14610*10^7 : + ( 1 + 2*corr_w_x_2(i,k)**2 ) * varnce_x_2(i,k) ) * varnce_w_2(i,k) )
2099 : end do
2100 : end do
2101 : !$acc end parallel loop
2102 :
2103 35740224 : return
2104 :
2105 : end subroutine calc_wp2xp2_pdf
2106 :
2107 : !=============================================================================
2108 35740224 : subroutine calc_wp2xp_pdf( nz, ngrdcol, &
2109 35740224 : wm, xm, w_1, w_2, &
2110 35740224 : x_1, x_2, &
2111 35740224 : varnce_w_1, varnce_w_2, &
2112 35740224 : varnce_x_1, varnce_x_2, &
2113 35740224 : corr_w_x_1, corr_w_x_2, &
2114 35740224 : mixt_frac, &
2115 35740224 : wp2xp )
2116 :
2117 : ! Description:
2118 : ! Calculates <w'^2 x'> by integrating over the PDF of w and x. The integral
2119 : ! is:
2120 : !
2121 : ! <w'^2 x'>
2122 : ! = INT(-inf:inf) INT(-inf:inf) ( w - <w> )^2 ( x - <x> ) P(w,x) dx dw;
2123 : !
2124 : ! where <w> is the overall mean of w, <x> is the overall mean of x, and
2125 : ! P(w,x) is a two-component bivariate normal distribution of w and x. The
2126 : ! integrated equation is:
2127 : !
2128 : ! <w'^2 x'>
2129 : ! = mixt_frac * ( ( mu_x_1 - <x> ) * ( ( mu_w_1 - <w> )^2 + sigma_w_1^2 )
2130 : ! + 2 * corr_w_x_1 * sigma_w_1 * sigma_x_1
2131 : ! * ( mu_w_1 - <w> ) )
2132 : ! + ( 1 - mixt_frac ) * ( ( mu_x_2 - <x> )
2133 : ! * ( ( mu_w_2 - <w> )^2 + sigma_w_2^2 )
2134 : ! + 2 * corr_w_x_2 * sigma_w_2 * sigma_x_2
2135 : ! * ( mu_w_2 - <w> ) );
2136 : !
2137 : ! where mu_w_1 is the mean of w in the 1st PDF component, mu_w_2 is the mean
2138 : ! of w in the 2nd PDF component, mu_x_1 is the mean of x in the 1st PDF
2139 : ! component, mu_x_2 is the mean of x in the 2nd PDF component, sigma_w_1 is
2140 : ! the standard deviation of w in the 1st PDF component, sigma_w_2 is the
2141 : ! standard deviation of w in the 2nd PDF component, sigma_x_1 is the
2142 : ! standard deviation of x in the 1st PDF component, sigma_x_2 is the
2143 : ! standard deviation of x in the 2nd PDF component, corr_w_x_1 is the
2144 : ! correlation of w and x in the 1st PDF component, corr_w_x_2 is the
2145 : ! correlation of w and x in the 2nd PDF component, and mixt_frac is the
2146 : ! mixture fraction, which is the weight of the 1st PDF component.
2147 :
2148 : ! References:
2149 : !-----------------------------------------------------------------------
2150 :
2151 : use constants_clubb, only: &
2152 : two, & ! Variable(s)
2153 : one
2154 :
2155 : use clubb_precision, only: &
2156 : core_rknd ! Variable(s)
2157 :
2158 : implicit none
2159 :
2160 : integer, intent(in) :: &
2161 : ngrdcol, & ! Number of grid columns
2162 : nz ! Number of vertical level
2163 :
2164 : ! Input Variables
2165 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
2166 : wm, & ! Mean of w (overall) [m/s]
2167 : xm, & ! Mean of x (overall) [units vary]
2168 : w_1, & ! Mean of w (1st PDF component) [m/s]
2169 : w_2, & ! Mean of w (2nd PDF component) [m/s]
2170 : x_1, & ! Mean of x (1st PDF component) [units vary]
2171 : x_2, & ! Mean of x (2nd PDF component) [units vary]
2172 : varnce_w_1, & ! Variance of w (1st PDF component) [m^2/s^2]
2173 : varnce_w_2, & ! Variance of w (2nd PDF component) [m^2/s^2]
2174 : varnce_x_1, & ! Variance of x (1st PDF component) [(units vary)^2]
2175 : varnce_x_2, & ! Variance of x (2nd PDF component) [(units vary)^2]
2176 : corr_w_x_1, & ! Correlation of w and x (1st PDF comp.) [-]
2177 : corr_w_x_2, & ! Correlation of w and x (2nd PDF comp.) [-]
2178 : mixt_frac ! Mixture fraction [-]
2179 :
2180 : ! Output Variable
2181 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
2182 : wp2xp ! <w'^2 x'> [m^2/s^2 (units vary)]
2183 :
2184 : ! Local Variables
2185 : integer :: i, k
2186 :
2187 :
2188 : ! Calculate <w'^2 x'> by integrating over the PDF.
2189 : !$acc parallel loop gang vector collapse(2) default(present)
2190 3073659264 : do k = 1, nz
2191 50761923264 : do i = 1, ngrdcol
2192 :
2193 95376528000 : wp2xp(i,k) = mixt_frac(i,k) &
2194 : * ( ( ( w_1(i,k) - wm(i,k) )**2 + varnce_w_1(i,k) ) * ( x_1(i,k) - xm(i,k) ) &
2195 : + two * corr_w_x_1(i,k) * sqrt( varnce_w_1(i,k) * varnce_x_1(i,k) ) &
2196 : * ( w_1(i,k) - wm(i,k) ) ) &
2197 : + ( one - mixt_frac(i,k) ) &
2198 : * ( ( ( w_2(i,k) - wm(i,k) )**2 + varnce_w_2(i,k) ) * ( x_2(i,k) - xm(i,k) ) &
2199 : + two * corr_w_x_2(i,k) * sqrt( varnce_w_2(i,k) * varnce_x_2(i,k) ) &
2200 >14610*10^7 : * ( w_2(i,k) - wm(i,k) ) )
2201 : end do
2202 : end do
2203 : !$acc end parallel loop
2204 :
2205 35740224 : return
2206 :
2207 : end subroutine calc_wp2xp_pdf
2208 :
2209 : !=============================================================================
2210 35740224 : subroutine calc_wpxp2_pdf( nz, ngrdcol, &
2211 35740224 : wm, xm, w_1, &
2212 35740224 : w_2, x_1, x_2, &
2213 35740224 : varnce_w_1, varnce_w_2, &
2214 35740224 : varnce_x_1, varnce_x_2, &
2215 35740224 : corr_w_x_1, corr_w_x_2, &
2216 35740224 : mixt_frac, &
2217 35740224 : wpxp2 )
2218 :
2219 : ! Description:
2220 : ! Calculates <w'x'^2> by integrating over the PDF of w and x. The integral
2221 : ! is:
2222 : !
2223 : ! <w'x'^2>
2224 : ! = INT(-inf:inf) INT(-inf:inf) ( w - <w> ) ( x - <x> )^2 P(w,x) dx dw;
2225 : !
2226 : ! where <w> is the overall mean of w, <x> is the overall mean of x, and
2227 : ! P(w,x) is a two-component bivariate normal distribution of w and x. The
2228 : ! integrated equation is:
2229 : !
2230 : ! <w'x'^2>
2231 : ! = mixt_frac * ( ( mu_w_1 - <w> ) * ( ( mu_x_1 - <x> )^2 + sigma_x_1^2 )
2232 : ! + 2 * corr_w_x_1 * sigma_w_1 * sigma_x_1
2233 : ! * ( mu_x_1 - <x> ) )
2234 : ! + ( 1 - mixt_frac ) * ( ( mu_w_2 - <w> )
2235 : ! * ( ( mu_x_2 - <x> )^2 + sigma_x_2^2 )
2236 : ! + 2 * corr_w_x_2 * sigma_w_2 * sigma_x_2
2237 : ! * ( mu_x_2 - <x> ) );
2238 : !
2239 : ! where mu_w_1 is the mean of w in the 1st PDF component, mu_w_2 is the mean
2240 : ! of w in the 2nd PDF component, mu_x_1 is the mean of x in the 1st PDF
2241 : ! component, mu_x_2 is the mean of x in the 2nd PDF component, sigma_w_1 is
2242 : ! the standard deviation of w in the 1st PDF component, sigma_w_2 is the
2243 : ! standard deviation of w in the 2nd PDF component, sigma_x_1 is the
2244 : ! standard deviation of x in the 1st PDF component, sigma_x_2 is the
2245 : ! standard deviation of x in the 2nd PDF component, corr_w_x_1 is the
2246 : ! correlation of w and x in the 1st PDF component, corr_w_x_2 is the
2247 : ! correlation of w and x in the 2nd PDF component, and mixt_frac is the
2248 : ! mixture fraction, which is the weight of the 1st PDF component.
2249 :
2250 : ! References:
2251 : !-----------------------------------------------------------------------
2252 :
2253 : use constants_clubb, only: &
2254 : two, & ! Variable(s)
2255 : one
2256 :
2257 : use clubb_precision, only: &
2258 : core_rknd ! Variable(s)
2259 :
2260 : implicit none
2261 :
2262 : integer, intent(in) :: &
2263 : ngrdcol, & ! Number of grid columns
2264 : nz ! Number of vertical level
2265 :
2266 : ! Input Variables
2267 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
2268 : wm, & ! Mean of w (overall) [m/s]
2269 : xm, & ! Mean of x (overall) [units vary]
2270 : w_1, & ! Mean of w (1st PDF component) [m/s]
2271 : w_2, & ! Mean of w (2nd PDF component) [m/s]
2272 : x_1, & ! Mean of x (1st PDF component) [units vary]
2273 : x_2, & ! Mean of x (2nd PDF component) [units vary]
2274 : varnce_w_1, & ! Variance of w (1st PDF component) [m^2/s^2]
2275 : varnce_w_2, & ! Variance of w (2nd PDF component) [m^2/s^2]
2276 : varnce_x_1, & ! Variance of x (1st PDF component) [(units vary)^2]
2277 : varnce_x_2, & ! Variance of x (2nd PDF component) [(units vary)^2]
2278 : corr_w_x_1, & ! Correlation of w and x (1st PDF comp.) [-]
2279 : corr_w_x_2, & ! Correlation of w and x (2nd PDF comp.) [-]
2280 : mixt_frac ! Mixture fraction [-]
2281 :
2282 : ! Return Variable
2283 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
2284 : wpxp2 ! <w'x'^2> [m/s (units vary)^2]
2285 :
2286 : ! Local Variables
2287 : integer :: i, k
2288 :
2289 : !$acc parallel loop gang vector collapse(2) default(present)
2290 3073659264 : do k = 1, nz
2291 50761923264 : do i = 1, ngrdcol
2292 :
2293 : ! Calculate <w'x'^2> by integrating over the PDF.
2294 95376528000 : wpxp2(i,k) = mixt_frac(i,k) &
2295 : * ( ( w_1(i,k) - wm(i,k) ) * ( ( x_1(i,k) - xm(i,k) )**2 + varnce_x_1(i,k) ) &
2296 : + two * corr_w_x_1(i,k) * sqrt( varnce_w_1(i,k) * varnce_x_1(i,k) ) &
2297 : * ( x_1(i,k) - xm(i,k) ) ) &
2298 : + ( one - mixt_frac(i,k) ) &
2299 : * ( ( w_2(i,k) - wm(i,k) ) * ( ( x_2(i,k) - xm(i,k) )**2 + varnce_x_2(i,k) ) &
2300 : + two * corr_w_x_2(i,k) * sqrt( varnce_w_2(i,k) * varnce_x_2(i,k) ) &
2301 >14610*10^7 : * ( x_2(i,k) - xm(i,k) ) )
2302 : end do
2303 : end do
2304 : !$acc end parallel loop
2305 :
2306 35740224 : return
2307 :
2308 : end subroutine calc_wpxp2_pdf
2309 :
2310 : !=============================================================================
2311 0 : subroutine calc_wpxpyp_pdf( nz, ngrdcol, &
2312 0 : wm, xm, ym, w_1, w_2, &
2313 0 : x_1, x_2, &
2314 0 : y_1, y_2, &
2315 0 : varnce_w_1, varnce_w_2, &
2316 0 : varnce_x_1, varnce_x_2, &
2317 0 : varnce_y_1, varnce_y_2, &
2318 0 : corr_w_x_1, corr_w_x_2, &
2319 0 : corr_w_y_1, corr_w_y_2, &
2320 0 : corr_x_y_1, corr_x_y_2, &
2321 0 : mixt_frac, &
2322 0 : wpxpyp )
2323 :
2324 : ! Description:
2325 : ! Calculates <w'x'y'> by integrating over the PDF of w, x, and y. The
2326 : ! integral is:
2327 : !
2328 : ! <w'x'y'>
2329 : ! = INT(-inf:inf) INT(-inf:inf) INT(-inf:inf)
2330 : ! ( w - <w> ) ( x - <x> ) ( y - <y> ) P(w,x,y) dy dx dw;
2331 : !
2332 : ! where <w> is the overall mean of w, <x> is the overall mean of x, <y> is
2333 : ! the overall mean of y, and P(w,x,y) is a two-component trivariate normal
2334 : ! distribution of w, x, and y. The integrated equation is:
2335 : !
2336 : ! <w'x'y'>
2337 : ! = mixt_frac
2338 : ! * ( ( mu_w_1 - <w> ) * ( mu_x_1 - <x> ) * ( mu_y_1 - <y> )
2339 : ! + corr_x_y_1 * sigma_x_1 * sigma_y_1 * ( mu_w_1 - <w> )
2340 : ! + corr_w_y_1 * sigma_w_1 * sigma_y_1 * ( mu_x_1 - <x> )
2341 : ! + corr_w_x_1 * sigma_w_1 * sigma_x_1 * ( mu_y_1 - <y> ) )
2342 : ! + ( 1 - mixt_frac )
2343 : ! * ( ( mu_w_2 - <w> ) * ( mu_x_2 - <x> ) * ( mu_y_2 - <y> )
2344 : ! + corr_x_y_2 * sigma_x_2 * sigma_y_2 * ( mu_w_2 - <w> )
2345 : ! + corr_w_y_2 * sigma_w_2 * sigma_y_2 * ( mu_x_2 - <x> )
2346 : ! + corr_w_x_2 * sigma_w_2 * sigma_x_2 * ( mu_y_2 - <y> ) );
2347 : !
2348 : ! where mu_w_1 is the mean of w in the 1st PDF component, mu_w_2 is the mean
2349 : ! of w in the 2nd PDF component, mu_x_1 is the mean of x in the 1st PDF
2350 : ! component, mu_x_2 is the mean of x in the 2nd PDF component, mu_y_1 is the
2351 : ! mean of y in the 1st PDF component, mu_y_2 is the mean of y in the 2nd PDF
2352 : ! component, sigma_w_1 is the standard deviation of w in the 1st PDF
2353 : ! component, sigma_w_2 is the standard deviation of w in the 2nd PDF
2354 : ! component, sigma_x_1 is the standard deviation of x in the 1st PDF
2355 : ! component, sigma_x_2 is the standard deviation of x in the 2nd PDF
2356 : ! component, sigma_y_1 is the standard deviation of y in the 1st PDF
2357 : ! component, sigma_y_2 is the standard deviation of y in the 2nd PDF
2358 : ! component, corr_w_x_1 is the correlation of w and x in the 1st PDF
2359 : ! component, corr_w_x_2 is the correlation of w and x in the 2nd PDF
2360 : ! component, corr_w_y_1 is the correlation of w and y in the 1st PDF
2361 : ! component, corr_w_y_2 is the correlation of w and y in the 2nd PDF
2362 : ! component, corr_x_y_1 is the correlation of x and y in the 1st PDF
2363 : ! component, corr_x_y_2 is the correlation of x and y in the 2nd PDF
2364 : ! component, and mixt_frac is the mixture fraction, which is the weight of
2365 : ! the 1st PDF component.
2366 :
2367 : ! References:
2368 : !-----------------------------------------------------------------------
2369 :
2370 : use grid_class, only: &
2371 : grid ! Type
2372 :
2373 : use constants_clubb, only: &
2374 : one ! Variable(s)
2375 :
2376 : use clubb_precision, only: &
2377 : core_rknd ! Variable(s)
2378 :
2379 : implicit none
2380 :
2381 : integer, intent(in) :: &
2382 : ngrdcol, & ! Number of grid columns
2383 : nz ! Number of vertical level
2384 :
2385 : ! Input Variables
2386 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
2387 : wm, & ! Mean of w (overall) [m/s]
2388 : xm, & ! Mean of x (overall) [x units]
2389 : ym, & ! Mean of y (overall) [y units]
2390 : w_1, & ! Mean of w (1st PDF component) [m/s]
2391 : w_2, & ! Mean of w (2nd PDF component) [m/s]
2392 : x_1, & ! Mean of x (1st PDF component) [x units]
2393 : x_2, & ! Mean of x (2nd PDF component) [x units]
2394 : y_1, & ! Mean of y (1st PDF component) [y units]
2395 : y_2, & ! Mean of y (2nd PDF component) [y units]
2396 : varnce_w_1, & ! Variance of w (1st PDF component) [m^2/s^2]
2397 : varnce_w_2, & ! Variance of w (2nd PDF component) [m^2/s^2]
2398 : varnce_x_1, & ! Variance of x (1st PDF component) [(x units)^2]
2399 : varnce_x_2, & ! Variance of x (2nd PDF component) [(x units)^2]
2400 : varnce_y_1, & ! Variance of y (1st PDF component) [(y units)^2]
2401 : varnce_y_2, & ! Variance of y (2nd PDF component) [(y units)^2]
2402 : corr_w_x_1, & ! Correlation of w and x (1st PDF component) [-]
2403 : corr_w_x_2, & ! Correlation of w and x (2nd PDF component) [-]
2404 : corr_w_y_1, & ! Correlation of w and y (1st PDF component) [-]
2405 : corr_w_y_2, & ! Correlation of w and y (2nd PDF component) [-]
2406 : corr_x_y_1, & ! Correlation of x and y (1st PDF component) [-]
2407 : corr_x_y_2, & ! Correlation of x and y (2nd PDF component) [-]
2408 : mixt_frac ! Mixture fraction [-]
2409 :
2410 : ! Output Variable
2411 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
2412 : wpxpyp ! <w'x'y'> [m/s (units vary)]
2413 :
2414 : ! Local Variables
2415 : integer :: i, k
2416 :
2417 :
2418 : ! Calculate <w'x'y'> by integrating over the PDF.
2419 : !$acc parallel loop gang vector collapse(2) default(present)
2420 0 : do k = 1, nz
2421 0 : do i = 1, ngrdcol
2422 0 : wpxpyp(i,k) &
2423 : = mixt_frac(i,k) &
2424 : * ( ( w_1(i,k) - wm(i,k) ) * ( x_1(i,k) - xm(i,k) ) * ( y_1(i,k) - ym(i,k) ) &
2425 : + corr_x_y_1(i,k)*sqrt( varnce_x_1(i,k)*varnce_y_1(i,k) )*( w_1(i,k)-wm(i,k) ) &
2426 : + corr_w_y_1(i,k)*sqrt( varnce_w_1(i,k)*varnce_y_1(i,k) )*( x_1(i,k)-xm(i,k) ) &
2427 : + corr_w_x_1(i,k)*sqrt( varnce_w_1(i,k)*varnce_x_1(i,k) )*( y_1(i,k)-ym(i,k) ) ) &
2428 : + ( one - mixt_frac(i,k) ) &
2429 : * ( ( w_2(i,k) - wm(i,k) )*( x_2(i,k) - xm(i,k) ) * ( y_2(i,k) - ym(i,k) ) &
2430 : + corr_x_y_2(i,k)*sqrt( varnce_x_2(i,k)*varnce_y_2(i,k) )*( w_2(i,k)-wm(i,k) ) &
2431 : + corr_w_y_2(i,k)*sqrt( varnce_w_2(i,k)*varnce_y_2(i,k) )*( x_2(i,k)-xm(i,k) ) &
2432 0 : + corr_w_x_2(i,k)*sqrt( varnce_w_2(i,k)*varnce_x_2(i,k) )*( y_2(i,k)-ym(i,k) ) )
2433 : end do
2434 : end do
2435 : !$acc end parallel loop
2436 :
2437 0 : return
2438 :
2439 : end subroutine calc_wpxpyp_pdf
2440 :
2441 : !=============================================================================
2442 35740224 : subroutine calc_liquid_cloud_frac_component( nz, ngrdcol, &
2443 35740224 : mean_chi, stdev_chi, &
2444 35740224 : cloud_frac, rc )
2445 : ! Description:
2446 : ! Calculates the PDF component cloud water mixing ratio, rc_i, and cloud
2447 : ! fraction, cloud_frac_i, for the ith PDF component.
2448 : !
2449 : ! The equation for cloud water mixing ratio, rc, at any point is:
2450 : !
2451 : ! rc = chi * H(chi);
2452 : !
2453 : ! and the equation for cloud fraction at a point, fc, is:
2454 : !
2455 : ! fc = H(chi);
2456 : !
2457 : ! where where extended liquid water mixing ratio, chi, is equal to cloud
2458 : ! water mixing ratio, rc, when positive. When the atmosphere is saturated
2459 : ! at this point, cloud water is found, and rc = chi, while fc = 1.
2460 : ! Otherwise, clear air is found at this point, and rc = fc = 0.
2461 : !
2462 : ! The mean of rc and fc is calculated by integrating over the PDF, such
2463 : ! that:
2464 : !
2465 : ! <rc> = INT(-inf:inf) chi * H(chi) * P(chi) dchi; and
2466 : !
2467 : ! cloud_frac = <fc> = INT(-inf:inf) H(chi) * P(chi) dchi.
2468 : !
2469 : ! This can be rewritten as:
2470 : !
2471 : ! <rc> = INT(0:inf) chi * P(chi) dchi; and
2472 : !
2473 : ! cloud_frac = <fc> = INT(0:inf) P(chi) dchi;
2474 : !
2475 : ! and further rewritten as:
2476 : !
2477 : ! <rc> = SUM(i=1,N) mixt_frac_i INT(0:inf) chi * P_i(chi) dchi; and
2478 : !
2479 : ! cloud_frac = SUM(i=1,N) mixt_frac_i INT(0:inf) P_i(chi) dchi;
2480 : !
2481 : ! where N is the number of PDF components. The equation for mean rc in the
2482 : ! ith PDF component is:
2483 : !
2484 : ! rc_i = INT(0:inf) chi * P_i(chi) dchi;
2485 : !
2486 : ! and the equation for cloud fraction in the ith PDF component is:
2487 : !
2488 : ! cloud_frac_i = INT(0:inf) P_i(chi) dchi.
2489 : !
2490 : ! The component values are related to the overall values by:
2491 : !
2492 : ! <rc> = SUM(i=1,N) mixt_frac_i * rc_i; and
2493 : !
2494 : ! cloud_frac = SUM(i=1,N) mixt_frac_i * cloud_frac_i.
2495 :
2496 : ! References:
2497 : !----------------------------------------------------------------------
2498 :
2499 : use constants_clubb, only: &
2500 : chi_tol, & ! Tolerance for pdf parameter chi [kg/kg]
2501 : sqrt_2pi, & ! sqrt(2*pi)
2502 : sqrt_2, & ! sqrt(2)
2503 : one, & ! 1
2504 : one_half, & ! 1/2
2505 : zero, & ! 0
2506 : max_num_stdevs, &
2507 : eps
2508 :
2509 : use clubb_precision, only: &
2510 : core_rknd ! Precision
2511 :
2512 : implicit none
2513 :
2514 : integer, intent(in) :: &
2515 : ngrdcol, & ! Number of grid columns
2516 : nz ! Number of vertical level
2517 :
2518 : !----------- Input Variables -----------
2519 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
2520 : mean_chi, & ! Mean of chi (old s) (ith PDF component) [kg/kg]
2521 : stdev_chi ! Standard deviation of chi (ith PDF component) [kg/kg]
2522 :
2523 : !----------- Output Variables -----------
2524 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
2525 : cloud_frac, & ! Cloud fraction (ith PDF component) [-]
2526 : rc ! Mean cloud water mixing ratio (ith PDF comp.) [kg/kg]
2527 :
2528 : !----------- Local Variables -----------
2529 : real( kind = core_rknd), parameter :: &
2530 : invrs_sqrt_2 = one / sqrt_2, &
2531 : invrs_sqrt_2pi = one / sqrt_2pi
2532 :
2533 : real( kind = core_rknd ) :: &
2534 : zeta
2535 :
2536 : integer :: k, i ! Vertical loop index
2537 :
2538 : !----------- Begin Code -----------
2539 : !$acc parallel loop gang vector collapse(2) default(present)
2540 3073659264 : do k = 1, nz
2541 50761923264 : do i = 1, ngrdcol
2542 :
2543 95376528000 : if ( ( abs( mean_chi(i,k) ) <= eps .and. stdev_chi(i,k) <= chi_tol ) &
2544 >14610*10^7 : .or. ( mean_chi(i,k) < - max_num_stdevs * stdev_chi(i,k) ) ) then
2545 :
2546 : ! The mean of chi is at saturation and does not vary in the ith PDF component
2547 44856697767 : cloud_frac(i,k) = zero
2548 44856697767 : rc(i,k) = zero
2549 :
2550 2831566233 : elseif ( mean_chi(i,k) > max_num_stdevs * stdev_chi(i,k) ) then
2551 :
2552 : ! The mean of chi is multiple standard deviations above the saturation point.
2553 : ! Thus, all cloud in the ith PDF component.
2554 297535127 : cloud_frac(i,k) = one
2555 297535127 : rc(i,k) = mean_chi(i,k)
2556 :
2557 : else
2558 :
2559 : ! The mean of chi is within max_num_stdevs of the saturation point.
2560 : ! Thus, layer is partly cloudy, requires calculation.
2561 :
2562 2534031106 : zeta = mean_chi(i,k) / stdev_chi(i,k)
2563 :
2564 2534031106 : cloud_frac(i,k) = one_half * ( one + erf( zeta * invrs_sqrt_2 ) )
2565 :
2566 : rc(i,k) = mean_chi(i,k) * cloud_frac(i,k) &
2567 2534031106 : + stdev_chi(i,k) * exp( - one_half * zeta**2 ) * invrs_sqrt_2pi
2568 :
2569 : end if
2570 :
2571 : end do
2572 : end do
2573 : !$acc end parallel loop
2574 :
2575 35740224 : return
2576 :
2577 : end subroutine calc_liquid_cloud_frac_component
2578 :
2579 : !=============================================================================
2580 35740224 : subroutine calc_ice_cloud_frac_component( nz, ngrdcol, &
2581 35740224 : mean_chi, stdev_chi, &
2582 35740224 : rc_in, cloud_frac, &
2583 35740224 : p_in_Pa, tl, &
2584 35740224 : rsatl, crt, &
2585 : saturation_formula, &
2586 35740224 : ice_supersat_frac, rc )
2587 : ! Description:
2588 : ! A version of the cloud fraction calculation modified to work
2589 : ! for layers that are potentially below freezing. If there are
2590 : ! no below freezing levels, the ice_supersat_frac calculation is
2591 : ! the same as cloud_frac.
2592 : !
2593 : ! For the below freezing levels, the saturation point will be
2594 : ! non-zero, thus we need to calculate chi_at_ice_sat.
2595 : !
2596 : ! The description of the equations are located in the description
2597 : ! of calc_liquid_cloud_frac_component.
2598 : !----------------------------------------------------------------------
2599 :
2600 : use constants_clubb, only: &
2601 : chi_tol, & ! Tolerance for pdf parameter chi [kg/kg]
2602 : T_freeze_K, & ! Freezing point of water [K]
2603 : sqrt_2pi, & ! sqrt(2*pi)
2604 : sqrt_2, & ! sqrt(2)
2605 : one, & ! 1
2606 : one_half, & ! 1/2
2607 : zero, & ! 0
2608 : max_num_stdevs, &
2609 : eps
2610 :
2611 : use clubb_precision, only: &
2612 : core_rknd ! Precision
2613 :
2614 : use saturation, only: &
2615 : sat_mixrat_ice
2616 :
2617 : implicit none
2618 :
2619 : ! ---------------------- Input Variables ----------------------
2620 : integer, intent(in) :: &
2621 : ngrdcol, & ! Number of grid columns
2622 : nz ! Number of vertical level
2623 :
2624 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
2625 : mean_chi, & ! Mean of chi (old s) (ith PDF component) [kg/kg]
2626 : stdev_chi, & ! Standard deviation of chi (ith PDF component) [kg/kg]
2627 : rc_in, & ! Mean cloud water mixing ratio (ith PDF comp.) [kg/kg]
2628 : cloud_frac, & ! Cloud fraction [-]
2629 : p_in_Pa, & ! Pressure [Pa]
2630 : rsatl, & ! Saturation mixing ratio of liquid [kg/kg]
2631 : crt, & ! r_t coef. in chi/eta eqns. [-]
2632 : tl ! Quantities needed to predict higher order moments
2633 : ! tl = thl*exner
2634 :
2635 : integer, intent(in) :: &
2636 : saturation_formula ! Integer that stores the saturation formula to be used
2637 :
2638 : ! ---------------------- Output Variables ----------------------
2639 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
2640 : ice_supersat_frac, & ! Ice supersaturation fraction [-]
2641 : rc ! Mean cloud ice mixing ratio (ith PDF comp.) [kg/kg]
2642 :
2643 : ! ---------------------- Local Variables----------------------
2644 : real( kind = core_rknd), parameter :: &
2645 : invrs_sqrt_2 = one / sqrt_2, &
2646 : invrs_sqrt_2pi = one / sqrt_2pi
2647 :
2648 : real( kind = core_rknd ) :: &
2649 : zeta, &
2650 : chi_at_ice_sat
2651 :
2652 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
2653 71480448 : rsat_ice
2654 :
2655 : integer :: k, i ! Loop indices
2656 :
2657 : logical :: &
2658 : l_any_below_freezing
2659 :
2660 : ! ---------------------- Begin Code ----------------------
2661 :
2662 35740224 : l_any_below_freezing = .false.
2663 :
2664 : ! If a grid boxes is above freezing, then the calculation is the
2665 : ! same as the cloud_frac calculation
2666 : !$acc parallel loop gang vector collapse(2) default(present) &
2667 : !$acc reduction(.or.:l_any_below_freezing)
2668 3073659264 : do k = 1, nz
2669 50761923264 : do i = 1, ngrdcol
2670 50726183040 : if ( tl(i,k) > T_freeze_K ) then
2671 9318427162 : ice_supersat_frac(i,k) = cloud_frac(i,k)
2672 9318427162 : rc(i,k) = rc_in(i,k)
2673 : else
2674 : l_any_below_freezing = .true.
2675 : end if
2676 : end do
2677 : end do
2678 : !$acc end parallel loop
2679 :
2680 : ! If all grid boxes are above freezing, then the calculation is the
2681 : ! same as the cloud_frac calculation
2682 35740224 : if ( .not. l_any_below_freezing ) then
2683 : return
2684 : end if
2685 :
2686 : !$acc data create( rsat_ice )
2687 :
2688 : ! Calculate the saturation mixing ratio of ice
2689 35740224 : rsat_ice = sat_mixrat_ice( nz, ngrdcol, p_in_Pa, tl, saturation_formula )
2690 :
2691 : !$acc parallel loop gang vector collapse(2) default(present)
2692 3073659264 : do k = 1, nz
2693 50761923264 : do i = 1, ngrdcol
2694 :
2695 50726183040 : if ( tl(i,k) <= T_freeze_K ) then
2696 :
2697 : ! Temperature is freezing, we must compute chi_at_ice_sat and
2698 : ! calculate the new cloud_frac_component
2699 38369836838 : chi_at_ice_sat = crt(i,k) * ( rsat_ice(i,k) - rsatl(i,k) )
2700 :
2701 : if ( ( abs( mean_chi(i,k)-chi_at_ice_sat ) <= eps .and. stdev_chi(i,k) <= chi_tol ) &
2702 38369836838 : .or. ( mean_chi(i,k)-chi_at_ice_sat < - max_num_stdevs * stdev_chi(i,k) ) ) then
2703 :
2704 : ! The mean of chi is at saturation and does not vary in the ith PDF component
2705 35394308879 : ice_supersat_frac(i,k) = zero
2706 35394308879 : rc(i,k) = zero
2707 :
2708 2975527959 : elseif ( mean_chi(i,k)-chi_at_ice_sat > max_num_stdevs * stdev_chi(i,k) ) then
2709 :
2710 : ! The mean of chi is multiple standard deviations above the saturation point.
2711 : ! Thus, all cloud in the ith PDF component.
2712 1695417368 : ice_supersat_frac(i,k) = one
2713 1695417368 : rc(i,k) = mean_chi(i,k)-chi_at_ice_sat
2714 :
2715 : else
2716 :
2717 : ! The mean of chi is within max_num_stdevs of the saturation point.
2718 : ! Thus, layer is partly cloudy, requires calculation.
2719 :
2720 1280110591 : zeta = (mean_chi(i,k)-chi_at_ice_sat) / stdev_chi(i,k)
2721 :
2722 1280110591 : ice_supersat_frac(i,k) = one_half * ( one + erf( zeta * invrs_sqrt_2 ) )
2723 :
2724 : rc(i,k) = (mean_chi(i,k)-chi_at_ice_sat) * ice_supersat_frac(i,k) &
2725 1280110591 : + stdev_chi(i,k) * exp( - one_half * zeta**2 ) * invrs_sqrt_2pi
2726 :
2727 : end if
2728 :
2729 : end if
2730 :
2731 : end do
2732 : end do
2733 : !$acc end parallel loop
2734 :
2735 : !$acc end data
2736 :
2737 : return
2738 :
2739 : end subroutine calc_ice_cloud_frac_component
2740 :
2741 : !=============================================================================
2742 35740224 : subroutine calc_xprcp_component( nz, ngrdcol, & ! In
2743 35740224 : wm, rtm, thlm, um, vm, rcm, & ! In
2744 35740224 : w_i, rt_i, & ! In
2745 35740224 : thl_i, u_i, v_i, & ! In
2746 35740224 : varnce_w_i, chi_i, & ! In
2747 35740224 : stdev_chi_i, stdev_eta_i, & ! In
2748 35740224 : corr_w_chi_i, corr_chi_eta_i, & ! In
2749 : ! corr_u_w_i, corr_v_w_i, & ! In
2750 35740224 : crt_i, cthl_i, & ! In
2751 35740224 : rc_i, cloud_frac_i, iiPDF_type, & ! In
2752 35740224 : wprcp_contrib_comp_i, wp2rcp_contrib_comp_i, & ! Out
2753 35740224 : rtprcp_contrib_comp_i, thlprcp_contrib_comp_i, & ! Out
2754 35740224 : uprcp_contrib_comp_i, vprcp_contrib_comp_i ) ! Out
2755 :
2756 : ! Description:
2757 : ! Calculates the contribution to <w'rc'>, <w'^2 rc'>, <rt'rc'>, and
2758 : ! <thl'rc'> from the ith PDF component.
2759 : !
2760 : !
2761 : ! <w'rc'>
2762 : ! -------
2763 : !
2764 : ! The value of <w'rc'> is calculated by integrating over the PDF:
2765 : !
2766 : ! <w'rc'>
2767 : ! = INT(-inf:inf) INT(-inf:inf) INT(-inf:inf)
2768 : ! ( w - <w> ) ( rc - <rc> ) P(w,rt,thl) dthl drt dw;
2769 : !
2770 : ! where <w> is the overall mean of w, <rc> is the overall mean of rc, and
2771 : ! P(w,rt,thl) is a two-component trivariate normal distribution of w, rt,
2772 : ! and thl. This equation is rewritten as:
2773 : !
2774 : ! <w'rc'>
2775 : ! = mixt_frac
2776 : ! * INT(-inf:inf) INT(-inf:inf) INT(-inf:inf)
2777 : ! ( w - <w> ) ( rc - <rc> ) P_1(w,rt,thl) dthl drt dw
2778 : ! + ( 1 - mixt_frac )
2779 : ! * INT(-inf:inf) INT(-inf:inf) INT(-inf:inf)
2780 : ! ( w - <w> ) ( rc - <rc> ) P_2(w,rt,thl) dthl drt dw;
2781 : !
2782 : ! where mixt_frac is the mixture fraction, which is the weight of the 1st
2783 : ! PDF component, and where P_1(w,rt,thl) and P_2(w,rt,thl) are the equations
2784 : ! for the trivariate normal PDF of w, rt, and thl in the 1st and 2nd PDF
2785 : ! components, respectively. The contribution from the ith PDF component is:
2786 : !
2787 : ! INT(-inf:inf) INT(-inf:inf) INT(-inf:inf)
2788 : ! ( w - <w> ) ( rc - <rc> ) P_i(w,rt,thl) dthl drt dw;
2789 : !
2790 : ! where P_i(w,rt,thl) is the trivariate normal PDF of w, rt, and thl in the
2791 : ! ith PDF component. The PDF undergoes a PDF transformation in each PDF
2792 : ! component, which is a change of variables and a translation, stretching,
2793 : ! and rotation of the axes. The PDF becomes a trivariate normal PDF that is
2794 : ! written in terms of w, chi, and eta coordinates. Cloud water mixing
2795 : ! ratio, rc, is written in terms of extended liquid water mixing ratio, chi,
2796 : ! such that:
2797 : !
2798 : ! rc = chi H(chi);
2799 : !
2800 : ! where H(chi) is the Heaviside step function. The contribution from the
2801 : ! ith PDF component to <w'rc'> can be written as:
2802 : !
2803 : ! INT(-inf:inf) INT(-inf:inf)
2804 : ! ( w - <w> ) ( chi H(chi) - <rc> ) P_i(w,chi) dchi dw;
2805 : !
2806 : ! where P_i(w,chi) is the bivariate normal PDF of w and chi in the ith PDF
2807 : ! component. The solved equation for the <w'rc'> contribution from the ith
2808 : ! PDF component (wprcp_contrib_comp_i) is:
2809 : !
2810 : ! wprcp_contrib_comp_i
2811 : ! = INT(-inf:inf) INT(-inf:inf)
2812 : ! ( w - <w> ) ( chi H(chi) - <rc> ) P_i(w,chi) dchi dw
2813 : ! = ( mu_w_i - <w> )
2814 : ! * ( mu_chi_i * 1/2 * ( 1 + erf( mu_chi_i / ( sqrt(2) * sigma_chi_i ) ) )
2815 : ! + 1/sqrt(2*pi) * sigma_chi_i
2816 : ! * exp{ - mu_chi_i^2 / ( 2 * sigma_chi_i^2 ) } - <rc> )
2817 : ! + corr_w_chi_i * sigma_w_i * sigma_chi_i
2818 : ! * 1/2 * ( 1 + erf( mu_chi_i / ( sqrt(2) * sigma_chi_i ) ) );
2819 : !
2820 : ! where mu_w_i is the mean of w in the ith PDF component, mu_chi_i is the
2821 : ! mean of chi in the ith PDF component, sigma_w_i is the standard deviation
2822 : ! of w in the ith PDF component, sigma_chi_i is the standard deviation of
2823 : ! chi in the ith PDF component, and corr_w_chi_i is the correlation of w and
2824 : ! chi in the ith PDF component.
2825 : !
2826 : ! Special case: sigma_chi_i = 0.
2827 : !
2828 : ! In the special case that sigma_chi_i = 0, chi, as well as rc, are constant
2829 : ! in the ith PDF component. The equation becomes:
2830 : !
2831 : ! wprcp_contrib_comp_i
2832 : ! = | ( mu_w_i - <w> ) * ( mu_chi_i - <rc> ); when mu_chi_i > 0;
2833 : ! | ( mu_w_i - <w> ) * ( -<rc> ); when mu_chi_i <= 0.
2834 : !
2835 : !
2836 : ! <w'^2 rc'>
2837 : ! ----------
2838 : !
2839 : ! The value of <w'^2 rc'> is calculated by integrating over the PDF:
2840 : !
2841 : ! <w'^2 rc'>
2842 : ! = INT(-inf:inf) INT(-inf:inf) INT(-inf:inf)
2843 : ! ( w - <w> )^2 ( rc - <rc> ) P(w,rt,thl) dthl drt dw.
2844 : !
2845 : ! This equation is rewritten as:
2846 : !
2847 : ! <w'^2 rc'>
2848 : ! = mixt_frac
2849 : ! * INT(-inf:inf) INT(-inf:inf) INT(-inf:inf)
2850 : ! ( w - <w> )^2 ( rc - <rc> ) P_1(w,rt,thl) dthl drt dw
2851 : ! + ( 1 - mixt_frac )
2852 : ! * INT(-inf:inf) INT(-inf:inf) INT(-inf:inf)
2853 : ! ( w - <w> )^2 ( rc - <rc> ) P_2(w,rt,thl) dthl drt dw.
2854 : !
2855 : ! The contribution from the ith PDF component is:
2856 : !
2857 : ! INT(-inf:inf) INT(-inf:inf) INT(-inf:inf)
2858 : ! ( w - <w> )^2 ( rc - <rc> ) P_i(w,rt,thl) dthl drt dw.
2859 : !
2860 : ! The PDF undergoes a PDF transformation in each PDF component, and becomes
2861 : ! a trivariate normal PDF that is written in terms of w, chi, and eta
2862 : ! coordinates. The contribution from the ith PDF component to <w'^2 rc'>
2863 : ! can be written as:
2864 : !
2865 : ! INT(-inf:inf) INT(-inf:inf)
2866 : ! ( w - <w> )^2 ( chi H(chi) - <rc> ) P_i(w,chi) dchi dw.
2867 : !
2868 : ! The solved equation for the <w'^2 rc'> contribution from the ith PDF
2869 : ! component (wp2rcp_contrib_comp_i) is:
2870 : !
2871 : ! wp2rcp_contrib_comp_i
2872 : ! = INT(-inf:inf) INT(-inf:inf)
2873 : ! ( w - <w> )^2 ( chi H(chi) - <rc> ) P_i(w,chi) dchi dw
2874 : ! = ( ( mu_w_i - <w> )^2 + sigma_w_i^2 )
2875 : ! * ( mu_chi_i * 1/2 * ( 1 + erf( mu_chi_i / ( sqrt(2) * sigma_chi_i ) ) )
2876 : ! + 1/sqrt(2*pi) * sigma_chi_i
2877 : ! * exp{ - mu_chi_i^2 / ( 2 * sigma_chi_i^2 ) } - <rc> )
2878 : ! + ( mu_w_i - <w> ) * corr_w_chi_i * sigma_w_i * sigma_chi_i
2879 : ! * ( 1 + erf( mu_chi_i / ( sqrt(2) * sigma_chi_i ) ) )
2880 : ! + 1/sqrt(2*pi) * corr_w_chi_i^2 * sigma_w_i^2 * sigma_chi_i
2881 : ! * exp{ - mu_chi_i^2 / ( 2 * sigma_chi_i^2 ) }.
2882 : !
2883 : ! Special case: sigma_chi_i = 0.
2884 : !
2885 : ! In the special case that sigma_chi_i = 0, chi, as well as rc, are constant
2886 : ! in the ith PDF component. The equation becomes:
2887 : !
2888 : ! wp2rcp_contrib_comp_i
2889 : ! = | ( ( mu_w_i - <w> )^2 + sigma_w_i^2 ) * ( mu_chi_i - <rc> );
2890 : ! | when mu_chi_i > 0;
2891 : ! | ( ( mu_w_i - <w> )^2 + sigma_w_i^2 ) * ( -<rc> );
2892 : ! | when mu_chi_i <= 0.
2893 : !
2894 : !
2895 : ! <rt'rc'>
2896 : ! --------
2897 : !
2898 : ! The value of <rt'rc'> is calculated by integrating over the PDF:
2899 : !
2900 : ! <rt'rc'>
2901 : ! = INT(-inf:inf) INT(-inf:inf)
2902 : ! ( rt - <rt> ) ( rc - <rc> ) P(rt,thl) dthl drt;
2903 : !
2904 : ! where <rt> is the overall mean of rt, and where P(rt,thl) is a
2905 : ! two-component bivariate normal distribution of rt and thl. This equation
2906 : ! is rewritten as:
2907 : !
2908 : ! <rt'rc'>
2909 : ! = mixt_frac
2910 : ! * INT(-inf:inf) INT(-inf:inf)
2911 : ! ( rt - <rt> ) ( rc - <rc> ) P_1(rt,thl) dthl drt
2912 : ! + ( 1 - mixt_frac )
2913 : ! * INT(-inf:inf) INT(-inf:inf)
2914 : ! ( rt - <rt> ) ( rc - <rc> ) P_2(rt,thl) dthl drt;
2915 : !
2916 : ! where P_1(rt,thl) and P_2(rt,thl) are the equations for the bivariate
2917 : ! normal PDF of rt and thl in the 1st and 2nd PDF components, respectively.
2918 : ! The contribution from the ith PDF component is:
2919 : !
2920 : ! INT(-inf:inf) INT(-inf:inf)
2921 : ! ( rt - <rt> ) ( rc - <rc> ) P_i(rt,thl) dthl drt;
2922 : !
2923 : ! where P_i(rt,thl) is the bivariate normal PDF of rt and thl in the ith PDF
2924 : ! component. The PDF undergoes a PDF transformation in each PDF component,
2925 : ! and becomes a bivariate normal PDF that is written in terms of chi and
2926 : ! eta coordinates. Total water mixing ratio, rt, is rewritten in terms of
2927 : ! chi and eta by:
2928 : !
2929 : ! rt = mu_rt_i
2930 : ! + ( ( eta - mu_eta_i ) + ( chi - mu_chi_i ) ) / ( 2 * crt_i );
2931 : !
2932 : ! where mu_rt_i is the mean of rt in the ith PDF component, mu_eta_i is the
2933 : ! mean of eta in the ith PDF component, and crt_i is a coefficient on rt in
2934 : ! the chi/eta transformation equations. The contribution from the ith PDF
2935 : ! component to <rt'rc'> can be written as:
2936 : !
2937 : ! INT(-inf:inf) INT(-inf:inf)
2938 : ! ( mu_rt_i - <rt> + ( eta - mu_eta_i ) / ( 2 * crt_i )
2939 : ! + ( chi - mu_chi_i ) / ( 2 * crt_i ) )
2940 : ! * ( chi H(chi) - <rc> ) P_i(chi,eta) deta dchi;
2941 : !
2942 : ! where P_i(chi,eta) is the bivariate normal PDF of chi and eta in the ith
2943 : ! PDF component. The solved equation for the <rt'rc'> contribution from the
2944 : ! ith PDF component (rtprcp_contrib_comp_i) is:
2945 : !
2946 : ! rtprcp_contrib_comp_i
2947 : ! = INT(-inf:inf) INT(-inf:inf)
2948 : ! ( mu_rt_i - <rt> + ( eta - mu_eta_i ) / ( 2 * crt_i )
2949 : ! + ( chi - mu_chi_i ) / ( 2 * crt_i ) )
2950 : ! * ( chi H(chi) - <rc> ) P_i(chi,eta) deta dchi
2951 : ! = ( mu_rt_i - <rt> )
2952 : ! * ( mu_chi_i * 1/2 * ( 1 + erf( mu_chi_i / ( sqrt(2) * sigma_chi_i ) ) )
2953 : ! + 1/sqrt(2*pi) * sigma_chi_i
2954 : ! * exp{ - mu_chi_i^2 / ( 2 * sigma_chi_i^2 ) } - <rc> )
2955 : ! + ( corr_chi_eta_i * sigma_eta_i + sigma_chi_i ) / ( 2 * crt_i )
2956 : ! * sigma_chi_i
2957 : ! * 1/2 * ( 1 + erf( mu_chi_i / ( sqrt(2) * sigma_chi_i ) ) );
2958 : !
2959 : ! where sigma_eta_i is the standard deviation of eta in the ith PDF
2960 : ! component and corr_chi_eta_i is the correlation of chi and eta in the ith
2961 : ! PDF component.
2962 : !
2963 : ! Special case: sigma_chi_i = 0.
2964 : !
2965 : ! In the special case that sigma_chi_i = 0, chi, as well as rc, are constant
2966 : ! in the ith PDF component. The equation becomes:
2967 : !
2968 : ! rtprcp_contrib_comp_i
2969 : ! = | ( mu_rt_i - <rt> ) * ( mu_chi_i - <rc> ); when mu_chi_i > 0;
2970 : ! | ( mu_rt_i - <rt> ) * ( -<rc> ); when mu_chi_i <= 0.
2971 : !
2972 : !
2973 : ! <thl'rc'>
2974 : ! ---------
2975 : !
2976 : ! The value of <thl'rc'> is calculated by integrating over the PDF:
2977 : !
2978 : ! <thl'rc'>
2979 : ! = INT(-inf:inf) INT(-inf:inf)
2980 : ! ( thl - <thl> ) ( rc - <rc> ) P(rt,thl) dthl drt;
2981 : !
2982 : ! where <thl> is the overall mean of thl. This equation is rewritten as:
2983 : !
2984 : ! <thl'rc'>
2985 : ! = mixt_frac
2986 : ! * INT(-inf:inf) INT(-inf:inf)
2987 : ! ( thl - <thl> ) ( rc - <rc> ) P_1(rt,thl) dthl drt
2988 : ! + ( 1 - mixt_frac )
2989 : ! * INT(-inf:inf) INT(-inf:inf)
2990 : ! ( thl - <thl> ) ( rc - <rc> ) P_2(rt,thl) dthl drt.
2991 : !
2992 : ! The contribution from the ith PDF component is:
2993 : !
2994 : ! INT(-inf:inf) INT(-inf:inf)
2995 : ! ( thl - <thl> ) ( rc - <rc> ) P_i(rt,thl) dthl drt.
2996 : !
2997 : ! The PDF undergoes a PDF transformation in each PDF component, and becomes
2998 : ! a bivariate normal PDF that is written in terms of chi and eta
2999 : ! coordinates. Liquid water potential temperature, thl, is rewritten in
3000 : ! terms of chi and eta by:
3001 : !
3002 : ! thl = mu_thl_i
3003 : ! + ( ( eta - mu_eta_i ) - ( chi - mu_chi_i ) ) / ( 2 * cthl_i );
3004 : !
3005 : ! where mu_thl_i is the mean of thl in the ith PDF component and cthl_i is a
3006 : ! coefficient on thl in the chi/eta transformation equations. The
3007 : ! contribution from the ith PDF component to <thl'rc'> can be written as:
3008 : !
3009 : ! INT(-inf:inf) INT(-inf:inf)
3010 : ! ( mu_thl_i - <thl> + ( eta - mu_eta_i ) / ( 2 * cthl_i )
3011 : ! - ( chi - mu_chi_i ) / ( 2 * cthl_i ) )
3012 : ! * ( chi H(chi) - <rc> ) P_i(chi,eta) deta dchi.
3013 : !
3014 : ! The solved equation for the <thl'rc'> contribution from the ith PDF
3015 : ! component (thlprcp_contrib_comp_i) is:
3016 : !
3017 : ! thlprcp_contrib_comp_i
3018 : ! = INT(-inf:inf) INT(-inf:inf)
3019 : ! ( mu_thl_i - <thl> + ( eta - mu_eta_i ) / ( 2 * cthl_i )
3020 : ! - ( chi - mu_chi_i ) / ( 2 * cthl_i ) )
3021 : ! * ( chi H(chi) - <rc> ) P_i(chi,eta) deta dchi
3022 : ! = ( mu_thl_i - <thl> )
3023 : ! * ( mu_chi_i * 1/2 * ( 1 + erf( mu_chi_i / ( sqrt(2) * sigma_chi_i ) ) )
3024 : ! + 1/sqrt(2*pi) * sigma_chi_i
3025 : ! * exp{ - mu_chi_i^2 / ( 2 * sigma_chi_i^2 ) } - <rc> )
3026 : ! + ( corr_chi_eta_i * sigma_eta_i - sigma_chi_i ) / ( 2 * cthl_i )
3027 : ! * sigma_chi_i
3028 : ! * 1/2 * ( 1 + erf( mu_chi_i / ( sqrt(2) * sigma_chi_i ) ) ).
3029 : !
3030 : ! Special case: sigma_chi_i = 0.
3031 : !
3032 : ! In the special case that sigma_chi_i = 0, chi, as well as rc, are constant
3033 : ! in the ith PDF component. The equation becomes:
3034 : !
3035 : ! thlprcp_contrib_comp_i
3036 : ! = | ( mu_thl_i - <thl> ) * ( mu_chi_i - <rc> ); when mu_chi_i > 0;
3037 : ! | ( mu_thl_i - <thl> ) * ( -<rc> ); when mu_chi_i <= 0.
3038 : !
3039 : !
3040 : ! Use equations for PDF component cloud fraction cloud water mixing ratio
3041 : ! -----------------------------------------------------------------------
3042 : !
3043 : ! The equation for cloud fraction in the ith PDF component, fc_i, is:
3044 : !
3045 : ! fc_i = 1/2 * ( 1 + erf( mu_chi_i / ( sqrt(2) * sigma_chi_i ) ) ).
3046 : !
3047 : ! In the special case that sigma_chi_i = 0, the equation becomes:
3048 : !
3049 : ! fc_i = | 1; when mu_chi_i > 0;
3050 : ! | 0; when mu_chi_i <= 0.
3051 : !
3052 : ! The equation for mean cloud water mixing ratio in the ith PDF component,
3053 : ! rc_i, is:
3054 : !
3055 : ! rc_i
3056 : ! = mu_chi_i * 1/2 * ( 1 + erf( mu_chi_i / ( sqrt(2) * sigma_chi_i ) ) )
3057 : ! + 1/sqrt(2*pi) * sigma_chi_i
3058 : ! * exp{ - mu_chi_i^2 / ( 2 * sigma_chi_i^2 ) }
3059 : ! = mu_chi_i * fc_i
3060 : ! + 1/sqrt(2*pi) * sigma_chi_i
3061 : ! * exp{ - mu_chi_i^2 / ( 2 * sigma_chi_i^2 ) }.
3062 : !
3063 : ! In the special case that sigma_chi_i = 0, the equation becomes:
3064 : !
3065 : ! rc_i = | mu_chi_i; when mu_chi_i > 0;
3066 : ! | 0; when mu_chi_i <= 0.
3067 : !
3068 : ! The above equations can be substituted into the equations for
3069 : ! wprcp_contrib_comp_i, wp2rcp_contrib_comp_i, rtprcp_contrib_comp_i, and
3070 : ! thlprcp_contrib_comp_i. The new equations are:
3071 : !
3072 : ! wprcp_contrib_comp_i
3073 : ! = ( mu_w_i - <w> ) * ( rc_i - <rc> )
3074 : ! + corr_w_chi_i * sigma_w_i * sigma_chi_i * fc_i;
3075 : !
3076 : ! wp2rcp_contrib_comp_i
3077 : ! = ( ( mu_w_i - <w> )^2 + sigma_w_i^2 ) * ( rc_i - <rc> )
3078 : ! + 2 * ( mu_w_i - <w> ) * corr_w_chi_i * sigma_w_i * sigma_chi_i * fc_i
3079 : ! + 1/sqrt(2*pi) * corr_w_chi_i^2 * sigma_w_i^2 * sigma_chi_i
3080 : ! * exp{ - mu_chi_i^2 / ( 2 * sigma_chi_i^2 ) };
3081 : !
3082 : ! rtprcp_contrib_comp_i
3083 : ! = ( mu_rt_i - <rt> ) * ( rc_i - <rc> )
3084 : ! + ( corr_chi_eta_i * sigma_eta_i + sigma_chi_i ) / ( 2 * crt_i )
3085 : ! * sigma_chi_i * fc_i; and
3086 : !
3087 : ! thlprcp_contrib_comp_i
3088 : ! = ( mu_thl_i - <thl> ) * ( rc_i - <rc> )
3089 : ! + ( corr_chi_eta_i * sigma_eta_i - sigma_chi_i ) / ( 2 * cthl_i )
3090 : ! * sigma_chi_i * fc_i.
3091 : !
3092 : ! While the above equations reduce to their listed versions in the special
3093 : ! case that sigma_chi_i = 0, those versions are faster to calculate. When
3094 : ! mu_chi_i > 0, they are:
3095 : !
3096 : ! wprcp_contrib_comp_i = ( mu_w_i - <w> ) * ( mu_chi_i - <rc> );
3097 : ! wp2rcp_contrib_comp_i
3098 : ! = ( ( mu_w_i - <w> )^2 + sigma_w_i^2 ) * ( mu_chi_i - <rc> );
3099 : ! rtprcp_contrib_comp_i = ( mu_rt_i - <rt> ) * ( mu_chi_i - <rc> ); and
3100 : ! thlprcp_contrib_comp_i = ( mu_thl_i - <thl> ) * ( mu_chi_i - <rc> );
3101 : !
3102 : ! and when mu_chi_i <= 0, they are:
3103 : !
3104 : ! wprcp_contrib_comp_i = - ( mu_w_i - <w> ) * <rc>;
3105 : ! wp2rcp_contrib_comp_i = - ( ( mu_w_i - <w> )^2 + sigma_w_i^2 ) * <rc>;
3106 : ! rtprcp_contrib_comp_i = - ( mu_rt_i - <rt> ) * <rc>; and
3107 : ! thlprcp_contrib_comp_i = - ( mu_thl_i - <thl> ) * <rc>.
3108 :
3109 : ! References:
3110 : !-----------------------------------------------------------------------
3111 :
3112 : use grid_class, only: &
3113 : grid ! Type
3114 :
3115 : use constants_clubb, only: &
3116 : sqrt_2pi, & ! Variable(s)
3117 : two, &
3118 : zero, &
3119 : chi_tol
3120 :
3121 : use clubb_precision, only: &
3122 : core_rknd ! Variable(s)
3123 :
3124 : implicit none
3125 :
3126 : integer, intent(in) :: &
3127 : ngrdcol, & ! Number of grid columns
3128 : nz ! Number of vertical level
3129 :
3130 : ! Input Variables
3131 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
3132 : wm, & ! Mean of w (overall) [m/s]
3133 : rtm, & ! Mean of rt (overall) [kg/kg]
3134 : thlm, & ! Mean of thl (overall) [K]
3135 : um, & ! Mean of eastward wind (overall) [m/s]
3136 : vm, & ! Mean of northward wind (overall) [m/s]
3137 : rcm, & ! Mean of rc (overall) [kg/kg]
3138 : w_i, & ! Mean of w (ith PDF component) [m/s]
3139 : rt_i, & ! Mean of rt (ith PDF component) [kg/kg]
3140 : thl_i, & ! Mean of thl (ith PDF component) [K]
3141 : u_i, & ! Mean of eastward wind (ith PDF component) [m/s]
3142 : v_i, & ! Mean of northward wind (ith PDF component) [m/s]
3143 : varnce_w_i, & ! Variance of w (ith PDF component) [m^2/s^2]
3144 : chi_i, & ! Mean of chi (ith PDF component) [kg/kg]
3145 : stdev_chi_i, & ! Standard deviation of chi (ith PDF comp.) [kg/kg]
3146 : stdev_eta_i, & ! Standard deviation of eta (ith PDF comp.) [kg/kg]
3147 : corr_w_chi_i, & ! Correlation of w and chi (ith PDF component) [-]
3148 : corr_chi_eta_i, & ! Correlation of chi and eta (ith PDF comp.) [-]
3149 : ! corr_u_w_i, & ! Correlation of u and w (ith PDF component) [-]
3150 : ! corr_v_w_i, & ! Correlation of v and w (ith PDF component) [-]
3151 : crt_i, & ! Coef. on rt in chi/eta eqns. (ith PDF comp.) [-]
3152 : cthl_i, & ! Coef. on thl: chi/eta eqns. (ith PDF comp.) [kg/kg/K]
3153 : rc_i, & ! Mean of rc (ith PDF component) [kg/kg]
3154 : cloud_frac_i ! Cloud fraction (ith PDF component) [-]
3155 :
3156 : integer, intent(in) :: &
3157 : iiPDF_type ! Selected option for the two-component normal (double
3158 : ! Gaussian) PDF type to use for the w, rt, and theta-l (or
3159 : ! w, chi, and eta) portion of CLUBB's multivariate,
3160 : ! two-component PDF.
3161 :
3162 : ! Output Variables
3163 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
3164 : wprcp_contrib_comp_i, & ! <w'rc'> contrib. (ith PDF comp.) [m/s(kg/kg)]
3165 : wp2rcp_contrib_comp_i, & ! <w'^2rc'> contrib. (ith comp) [m^2/s^2(kg/kg)]
3166 : rtprcp_contrib_comp_i, & ! <rt'rc'> contrib. (ith PDF comp.) [kg^2/kg^2]
3167 : thlprcp_contrib_comp_i, & ! <thl'rc'> contrib. (ith PDF comp.) [K(kg/kg)]
3168 : uprcp_contrib_comp_i, & ! <u'rc'> contrib. (ith PDF comp.) [m/s(kg/kg)]
3169 : vprcp_contrib_comp_i ! <v'rc'> contrib. (ith PDF comp.) [m/s(kg/kg)]
3170 :
3171 : ! Local Variables
3172 : integer :: i, k
3173 :
3174 : ! ---------------------- Begin Code ------------------
3175 :
3176 : ! Changing these conditionals may result in inconsistencies with the conditional
3177 : ! statements located in calc_cloud_frac_component
3178 : !$acc parallel loop gang vector collapse(2) default(present)
3179 3073659264 : do k = 1, nz
3180 50761923264 : do i = 1, ngrdcol
3181 :
3182 47688264000 : wprcp_contrib_comp_i(i,k) = ( w_i(i,k) - wm(i,k) ) * ( rc_i(i,k) - rcm(i,k) )
3183 :
3184 : wp2rcp_contrib_comp_i(i,k) = ( ( w_i(i,k) - wm(i,k) )**2 + varnce_w_i(i,k) ) &
3185 47688264000 : * ( rc_i(i,k) - rcm(i,k) )
3186 :
3187 : rtprcp_contrib_comp_i(i,k) = ( rt_i(i,k) - rtm(i,k) ) * ( rc_i(i,k) - rcm(i,k) ) &
3188 : + ( corr_chi_eta_i(i,k) * stdev_eta_i(i,k) + stdev_chi_i(i,k) ) &
3189 47688264000 : / ( two * crt_i(i,k) ) * stdev_chi_i(i,k) * cloud_frac_i(i,k)
3190 :
3191 : thlprcp_contrib_comp_i(i,k) = ( thl_i(i,k) - thlm(i,k) ) * ( rc_i(i,k) - rcm(i,k) ) &
3192 : + ( corr_chi_eta_i(i,k) * stdev_eta_i(i,k) - stdev_chi_i(i,k) ) &
3193 47688264000 : / ( two * cthl_i(i,k) ) * stdev_chi_i(i,k) * cloud_frac_i(i,k)
3194 :
3195 47688264000 : uprcp_contrib_comp_i(i,k) = ( u_i(i,k) - um(i,k) ) * ( rc_i(i,k) - rcm(i,k) )
3196 :
3197 50726183040 : vprcp_contrib_comp_i(i,k) = ( v_i(i,k) - vm(i,k) ) * ( rc_i(i,k) - rcm(i,k) )
3198 :
3199 : end do
3200 : end do
3201 : !$acc end parallel loop
3202 :
3203 : ! If iiPDF_type isn't iiPDF_ADG1, iiPDF_ADG2, or iiPDF_new_hybrid, so
3204 : ! corr_w_chi_i /= 0 (and perhaps corr_u_w_i /= 0).
3205 35740224 : if ( .not. ( iiPDF_type == iiPDF_ADG1 .or. iiPDF_type == iiPDF_ADG2 &
3206 : .or. iiPDF_type == iiPDF_new_hybrid ) ) then
3207 :
3208 : ! Chi varies significantly in the ith PDF component (stdev_chi > chi_tol)
3209 : ! and there is some cloud (0 < cloud_frac <= 1)
3210 0 : do k = 1, nz
3211 0 : do i = 1, ngrdcol
3212 0 : if ( stdev_chi_i(i,k) > chi_tol .and. cloud_frac_i(i,k) > zero ) then
3213 :
3214 : wprcp_contrib_comp_i(i,k) = wprcp_contrib_comp_i(i,k) &
3215 : + corr_w_chi_i(i,k) * sqrt( varnce_w_i(i,k) ) &
3216 0 : * stdev_chi_i(i,k) * cloud_frac_i(i,k)
3217 :
3218 : wp2rcp_contrib_comp_i(i,k) = wp2rcp_contrib_comp_i(i,k) &
3219 : + two * ( w_i(i,k) - wm(i,k) ) * corr_w_chi_i(i,k) &
3220 : * sqrt( varnce_w_i(i,k) ) * stdev_chi_i(i,k) &
3221 : * cloud_frac_i(i,k) &
3222 : + corr_w_chi_i(i,k)**2 * varnce_w_i(i,k) &
3223 : * stdev_chi_i(i,k) &
3224 : * exp( -chi_i(i,k)**2 / ( two*stdev_chi_i(i,k)**2 ) ) &
3225 0 : / sqrt_2pi
3226 :
3227 : ! In principle, uprcp_contrib_comp_i might depend on corr_u_w_i here.
3228 : end if
3229 : end do
3230 : end do
3231 : end if
3232 :
3233 35740224 : return
3234 :
3235 : end subroutine calc_xprcp_component
3236 :
3237 : !=============================================================================
3238 0 : subroutine calc_w_up_in_cloud( nz, ngrdcol, &
3239 0 : mixt_frac, cloud_frac_1, cloud_frac_2, &
3240 0 : w_1, w_2, varnce_w_1, varnce_w_2, &
3241 0 : w_up_in_cloud, w_down_in_cloud, &
3242 0 : cloudy_updraft_frac, cloudy_downdraft_frac )
3243 :
3244 : ! Description:
3245 : ! Subroutine that computes the mean cloudy updraft (and also calculates
3246 : ! the mean cloudy downdraft).
3247 : !
3248 : ! In order to activate aerosol, we'd like to feed the activation scheme
3249 : ! a vertical velocity that's representative of cloudy updrafts. For skewed
3250 : ! layers, like cumulus layers, this might be an improvement over the square
3251 : ! root of wp2 that's currently used. At the same time, it would be simpler
3252 : ! and less expensive than feeding SILHS samples into the aerosol code
3253 : ! (see larson-group/e3sm#19 and larson-group/e3sm#26).
3254 : !
3255 : ! The formulas are only valid for certain PDFs in CLUBB (ADG1, ADG2,
3256 : ! new hybrid), hence we omit calculation if another PDF type is used.
3257 : !
3258 : ! References: https://www.overleaf.com/project/614a136d47846639af22ae34
3259 : !----------------------------------------------------------------------
3260 :
3261 : use constants_clubb, only: &
3262 : sqrt_2pi, & ! sqrt(2*pi)
3263 : sqrt_2, & ! sqrt(2)
3264 : one, & ! 1
3265 : one_half, & ! 1/2
3266 : zero, & ! 0
3267 : max_num_stdevs, &
3268 : eps
3269 :
3270 : use clubb_precision, only: &
3271 : core_rknd ! Precision
3272 :
3273 : implicit none
3274 :
3275 : integer, intent(in) :: &
3276 : ngrdcol, & ! Number of grid columns
3277 : nz ! Number of vertical level
3278 :
3279 : !----------- Input Variables -----------
3280 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
3281 : mixt_frac, & ! mixture fraction [-]
3282 : cloud_frac_1, & ! cloud fraction (1st PDF component) [-]
3283 : cloud_frac_2, & ! cloud fraction (2nd PDF component) [-]
3284 : w_1, & ! upward velocity (1st PDF component) [m/s]
3285 : w_2, & ! upward velocity (2nd PDF component) [m/s]
3286 : varnce_w_1, & ! standard deviation of w (1st PDF component) [m^2/s^2]
3287 : varnce_w_2 ! standard deviation of w (2nd PDF component) [m^2/s^2]
3288 :
3289 : !----------- Output Variables -----------
3290 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
3291 : w_up_in_cloud, & ! mean cloudy updraft speed [m/s]
3292 : w_down_in_cloud, & ! mean cloudy downdraft speed [m/s]
3293 : cloudy_updraft_frac, & ! cloudy updraft fraction [-]
3294 : cloudy_downdraft_frac ! cloudy downdraft fraction [-]
3295 :
3296 : !----------- Local Variables -----------
3297 : real( kind = core_rknd ) :: &
3298 : w_up_1, w_up_2, & ! integral of w and Heaviside fnc, where w > 0
3299 : w_down_1, w_down_2, & ! integral of w and Heaviside fnc, where w < 0
3300 : stdev_w_1, stdev_w_2, & ! Standard deviation of w
3301 : ratio_w_1, & ! mu_w_1 / ( sqrt(2) * sigma_w_1 )
3302 : ratio_w_2, & ! mu_w_2 / ( sqrt(2) * sigma_w_2 )
3303 : erf_ratio_w_1, & ! erf( ratio_w_1 )
3304 : erf_ratio_w_2, & ! erf( ratio_w_2 )
3305 : exp_neg_ratio_w_1_sqd, & ! exp( -ratio_w_1^2 )
3306 : exp_neg_ratio_w_2_sqd, & ! exp( -ratio_w_2^2 )
3307 : updraft_frac_1, & ! Fraction of 1st PDF comp. where w > 0
3308 : updraft_frac_2, & ! Fraction of 2nd PDF comp. where w > 0
3309 : downdraft_frac_1, & ! Fraction of 1st PDF comp. where w < 0
3310 : downdraft_frac_2 ! Fraction of 2nd PDF comp. where w < 0
3311 :
3312 : integer :: i, k
3313 :
3314 : !$acc parallel loop gang vector collapse(2) default(present)
3315 0 : do k = 1, nz
3316 0 : do i = 1, ngrdcol
3317 :
3318 0 : stdev_w_1 = sqrt(varnce_w_1(i,k))
3319 0 : stdev_w_2 = sqrt(varnce_w_2(i,k))
3320 :
3321 : ! Calculate quantities in the 1st PDF component.
3322 0 : if ( w_1(i,k) > max_num_stdevs * stdev_w_1 ) then
3323 :
3324 : ! The mean of w in the 1st PDF component is more than
3325 : ! max_num_stdevs standard deviations above 0.
3326 : ! The entire 1st PDF component is found in an updraft (w > 0).
3327 : w_up_1 = w_1(i,k)
3328 : updraft_frac_1 = one
3329 : w_down_1 = zero
3330 : downdraft_frac_1 = zero
3331 :
3332 0 : elseif ( w_1(i,k) < - max_num_stdevs * stdev_w_1 ) then
3333 :
3334 : ! The mean of w in the 1st PDF component is more than
3335 : ! max_num_stdevs standard deviations below 0.
3336 : ! The entire 1st PDF component is found in a downdraft (w < 0).
3337 : w_up_1 = zero
3338 : updraft_frac_1 = zero
3339 : w_down_1 = w_1(i,k)
3340 : downdraft_frac_1 = one
3341 :
3342 : else
3343 :
3344 : ! The 1st PDF component contains both updraft and downdraft.
3345 0 : ratio_w_1 = w_1(i,k) / ( sqrt_2 * max(eps, stdev_w_1) )
3346 0 : erf_ratio_w_1 = erf( ratio_w_1 )
3347 0 : exp_neg_ratio_w_1_sqd = exp( -ratio_w_1**2 )
3348 :
3349 : w_up_1 &
3350 : = one_half * w_1(i,k) * ( one + erf_ratio_w_1 ) &
3351 0 : + ( stdev_w_1 / sqrt_2pi ) * exp_neg_ratio_w_1_sqd
3352 :
3353 0 : updraft_frac_1 = one_half * ( one + erf_ratio_w_1 )
3354 :
3355 : w_down_1 &
3356 : = one_half * w_1(i,k) * ( one - erf_ratio_w_1 ) &
3357 0 : - ( stdev_w_1 / sqrt_2pi ) * exp_neg_ratio_w_1_sqd
3358 :
3359 : !downdraft_frac_1 = one_half * ( one - erf_ratio_w_1 )
3360 0 : downdraft_frac_1 = one - updraft_frac_1
3361 :
3362 : endif
3363 :
3364 : ! Calculate quantities in the 2nd PDF component.
3365 0 : if ( w_2(i,k) > max_num_stdevs * stdev_w_2 ) then
3366 :
3367 : ! The mean of w in the 2nd PDF component is more than
3368 : ! max_num_stdevs standard deviations above 0.
3369 : ! The entire 2nd PDF component is found in an updraft (w > 0).
3370 : w_up_2 = w_2(i,k)
3371 : updraft_frac_2 = one
3372 : w_down_2 = zero
3373 : downdraft_frac_2 = zero
3374 :
3375 0 : elseif ( w_2(i,k) < - max_num_stdevs * stdev_w_2 ) then
3376 :
3377 : ! The mean of w in the 2nd PDF component is more than
3378 : ! max_num_stdevs standard deviations below 0.
3379 : ! The entire 2nd PDF component is found in a downdraft (w < 0).
3380 : w_up_2 = zero
3381 : updraft_frac_2 = zero
3382 : w_down_2 = w_2(i,k)
3383 : downdraft_frac_2 = one
3384 :
3385 : else
3386 :
3387 : ! The 2nd PDF component contains both updraft and downdraft.
3388 0 : ratio_w_2 = w_2(i,k) / ( sqrt_2 * max(eps, stdev_w_2) )
3389 0 : erf_ratio_w_2 = erf( ratio_w_2 )
3390 0 : exp_neg_ratio_w_2_sqd = exp( -ratio_w_2**2 )
3391 :
3392 : w_up_2 &
3393 : = one_half * w_2(i,k) * ( one + erf_ratio_w_2 ) &
3394 0 : + ( stdev_w_2 / sqrt_2pi ) * exp_neg_ratio_w_2_sqd
3395 :
3396 0 : updraft_frac_2 = one_half * ( one + erf_ratio_w_2 )
3397 :
3398 : w_down_2 &
3399 : = one_half * w_2(i,k) * ( one - erf_ratio_w_2 ) &
3400 0 : - ( stdev_w_2 / sqrt_2pi ) * exp_neg_ratio_w_2_sqd
3401 :
3402 : !downdraft_frac_2 = one_half * ( one - erf_ratio_w_2 )
3403 0 : downdraft_frac_2 = one - updraft_frac_2
3404 :
3405 : endif
3406 :
3407 : ! Calculate the total cloudy updraft fraction.
3408 : cloudy_updraft_frac(i,k) &
3409 : = mixt_frac(i,k) * cloud_frac_1(i,k) * updraft_frac_1 &
3410 0 : + ( one - mixt_frac(i,k) ) * cloud_frac_2(i,k) * updraft_frac_2
3411 :
3412 : ! Calculate the total cloudy downdraft fraction.
3413 : cloudy_downdraft_frac(i,k) &
3414 : = mixt_frac(i,k) * cloud_frac_1(i,k) * downdraft_frac_1 &
3415 0 : + ( one - mixt_frac(i,k) ) * cloud_frac_2(i,k) * downdraft_frac_2
3416 :
3417 : ! Calculate the mean vertical velocity found in a cloudy updraft.
3418 : w_up_in_cloud(i,k) &
3419 : = ( mixt_frac(i,k) * cloud_frac_1(i,k) * w_up_1 &
3420 : + ( one - mixt_frac(i,k) ) * cloud_frac_2(i,k) * w_up_2 ) &
3421 0 : / max( eps, cloudy_updraft_frac(i,k) )
3422 :
3423 : ! Calculate the mean vertical velocity found in a cloudy downdraft.
3424 : w_down_in_cloud(i,k) &
3425 : = ( mixt_frac(i,k) * cloud_frac_1(i,k) * w_down_1 &
3426 : + ( one - mixt_frac(i,k) ) * cloud_frac_2(i,k) * w_down_2 ) &
3427 0 : / max( eps, cloudy_downdraft_frac(i,k) )
3428 :
3429 : end do
3430 : end do
3431 : !$acc end parallel loop
3432 :
3433 0 : return
3434 :
3435 : end subroutine calc_w_up_in_cloud
3436 :
3437 : !=============================================================================
3438 : function interp_var_array( n_points, nz, k, z_vals, var )
3439 :
3440 : ! Description:
3441 : ! Interpolates a variable to an array of values about a given level
3442 :
3443 : ! References
3444 : !-----------------------------------------------------------------------
3445 :
3446 : use clubb_precision, only: &
3447 : core_rknd ! Constant
3448 :
3449 : implicit none
3450 :
3451 : ! Input Variables
3452 : integer, intent(in) :: &
3453 : n_points, & ! Number of points to interpolate to (must be odd and >= 3)
3454 : nz, & ! Total number of vertical levels
3455 : k ! Center of interpolation array
3456 :
3457 : real( kind = core_rknd ), dimension(nz), intent(in) :: &
3458 : z_vals, & ! Height at each vertical level [m]
3459 : var ! Variable values on grid [units vary]
3460 :
3461 : ! Output Variables
3462 : real( kind = core_rknd ), dimension(n_points) :: &
3463 : interp_var_array ! Interpolated values of variable [units vary]
3464 :
3465 : ! Local Variables
3466 : real( kind = core_rknd ) :: &
3467 : dz ! Distance between vertical levels
3468 :
3469 : real( kind = core_rknd ) :: &
3470 : z_val ! Height at some sub-grid level
3471 :
3472 : integer :: &
3473 : i, & ! Loop iterator
3474 :
3475 : subgrid_lev_count ! Number of refined grid points located between
3476 : ! two defined grid levels
3477 :
3478 : !-----------------------------------------------------------------------
3479 :
3480 : !----- Begin Code -----
3481 :
3482 : ! Place a point at each of k-1, k, and k+1.
3483 : interp_var_array(1) = var_value_integer_height( nz, k-1, z_vals, var )
3484 : interp_var_array((n_points+1)/2) = var_value_integer_height( nz, k, z_vals, var )
3485 : interp_var_array(n_points) = var_value_integer_height( nz, k+1, z_vals, var )
3486 :
3487 : subgrid_lev_count = (n_points - 3) / 2
3488 :
3489 : ! Lower half
3490 : if ( k == 1 ) then
3491 : dz = (z_vals(2) - z_vals(1)) / real( subgrid_lev_count+1, kind=core_rknd )
3492 : else
3493 : dz = (z_vals(k) - z_vals(k-1)) / real( subgrid_lev_count+1, kind=core_rknd )
3494 : end if
3495 : do i=1, subgrid_lev_count
3496 : z_val = z_vals(k) - real( i, kind=core_rknd ) * dz
3497 : interp_var_array(1+i) &
3498 : = var_subgrid_interp( nz, k, z_vals, var, z_val, l_below=.true. )
3499 : end do
3500 :
3501 : ! Upper half
3502 : if ( k == nz ) then
3503 : dz = ( z_vals(nz) - z_vals(nz-1) ) / real( subgrid_lev_count+1, kind=core_rknd )
3504 : else
3505 : dz = ( z_vals(k+1) - z_vals(k) ) / real( subgrid_lev_count+1, kind=core_rknd )
3506 : end if
3507 : do i=1, (n_points-3)/2
3508 : z_val = z_vals(k) + real( i, kind=core_rknd ) * dz
3509 : interp_var_array((n_points+1)/2+i) &
3510 : = var_subgrid_interp( nz, k, z_vals, var, z_val, l_below=.false. )
3511 : end do
3512 :
3513 : return
3514 : end function interp_var_array
3515 :
3516 : !=============================================================================
3517 : function var_value_integer_height( nz, k, z_vals, var_grid_value ) result( var_value )
3518 :
3519 : ! Description
3520 : ! Returns the value of a variable at an integer height between 0 and
3521 : ! nz+1 inclusive, using extrapolation when k==0 or k==nz+1
3522 :
3523 : ! References
3524 : !-----------------------------------------------------------------------
3525 :
3526 : use clubb_precision, only: &
3527 : core_rknd ! Constant
3528 :
3529 : use interpolation, only: &
3530 : mono_cubic_interp ! Procedure
3531 :
3532 : implicit none
3533 :
3534 : ! Input Variables
3535 : integer, intent(in) :: &
3536 : nz, & ! Total number of vertical levels
3537 : k ! Level to resolve variable value
3538 :
3539 : real( kind = core_rknd ), dimension(nz), intent(in) :: &
3540 : z_vals, & ! Height at each vertical level [m]
3541 : var_grid_value ! Value of variable at each grid level [units vary]
3542 :
3543 : ! Output Variables
3544 : real( kind = core_rknd ) :: &
3545 : var_value ! Value of variable at height level [units vary]
3546 :
3547 : ! Local Variables
3548 : integer :: km1, k00, kp1, kp2
3549 : !-----------------------------------------------------------------------
3550 :
3551 : !----- Begin Code -----
3552 :
3553 : if ( k >= 1 .and. k <= nz ) then
3554 : ! This is the simple case. No extrapolation necessary.
3555 : var_value = var_grid_value(k)
3556 : else if ( k == 0 ) then
3557 : ! Extrapolate below the lower boundary
3558 : km1 = nz
3559 : k00 = 1
3560 : kp1 = 2
3561 : kp2 = 3
3562 : var_value = mono_cubic_interp( z_vals(1)-(z_vals(2)-z_vals(1)), &
3563 : km1, k00, kp1, kp2, &
3564 : z_vals(km1), z_vals(k00), z_vals(kp1), z_vals(kp2), &
3565 : var_grid_value(km1), var_grid_value(k00), &
3566 : var_grid_value(kp1), var_grid_value(kp2) )
3567 : else if ( k == nz+1 ) then
3568 : ! Extrapolate above the upper boundary
3569 : km1 = nz
3570 : k00 = nz-1
3571 : kp1 = nz
3572 : kp2 = nz
3573 : var_value = mono_cubic_interp( z_vals(nz)+(z_vals(nz)-z_vals(nz-1)), &
3574 : km1, k00, kp1, kp2, &
3575 : z_vals(km1), z_vals(k00), z_vals(kp1), z_vals(kp2), &
3576 : var_grid_value(km1), var_grid_value(k00), &
3577 : var_grid_value(kp1), var_grid_value(kp2) )
3578 : else
3579 : ! Invalid height requested
3580 : var_value = -999._core_rknd
3581 : end if ! k > 1 .and. k < nz
3582 : return
3583 : end function var_value_integer_height
3584 :
3585 : !=============================================================================
3586 : function var_subgrid_interp( nz, k, z_vals, var, z_interp, l_below ) result( var_value )
3587 :
3588 : ! Description
3589 : ! Interpolates (or extrapolates) a variable to a value between grid
3590 : ! levels
3591 :
3592 : ! References
3593 : !-----------------------------------------------------------------------
3594 :
3595 : use clubb_precision, only: &
3596 : core_rknd ! Constant
3597 :
3598 : use interpolation, only: &
3599 : mono_cubic_interp ! Procedure
3600 :
3601 : implicit none
3602 :
3603 : ! Input Variables
3604 : integer, intent(in) :: &
3605 : nz, & ! Number of vertical levels
3606 : k ! Grid level near interpolation target
3607 :
3608 : real( kind = core_rknd ), dimension(nz), intent(in) :: &
3609 : z_vals, & ! Height at each grid level [m]
3610 : var ! Variable values at grid levels [units vary]
3611 :
3612 : real( kind = core_rknd ), intent(in) :: &
3613 : z_interp ! Interpolation target height [m]
3614 :
3615 : logical, intent(in) :: &
3616 : l_below ! True if z_interp < z_vals(k), false otherwise
3617 :
3618 : ! Output Variable
3619 : real( kind = core_rknd ) :: &
3620 : var_value ! Interpolated value of variable [units vary]
3621 :
3622 : ! Local Variables
3623 : integer :: km1, k00, kp1, kp2 ! Parameters for call to mono_cubic_interp
3624 : !----------------------------------------------------------------------
3625 :
3626 : !----- Begin Code -----
3627 : if ( l_below ) then
3628 :
3629 : if ( k == 1 ) then ! Extrapolation
3630 : km1 = nz
3631 : k00 = 1
3632 : kp1 = 2
3633 : kp2 = 3
3634 : else if ( k == 2 ) then
3635 : km1 = 1
3636 : k00 = 1
3637 : kp1 = 2
3638 : kp2 = 3
3639 : else if ( k == nz ) then
3640 : km1 = nz-2
3641 : k00 = nz-1
3642 : kp1 = nz
3643 : kp2 = nz
3644 : else
3645 : km1 = k-2
3646 : k00 = k-1
3647 : kp1 = k
3648 : kp2 = k+1
3649 : end if ! k == 1
3650 :
3651 : else ! .not. l_below
3652 :
3653 : if ( k == 1 ) then
3654 : km1 = 1
3655 : k00 = 1
3656 : kp1 = 2
3657 : kp2 = 3
3658 : else if ( k == nz-1 ) then
3659 : km1 = nz-2
3660 : k00 = nz-1
3661 : kp1 = nz
3662 : kp2 = nz
3663 : else if ( k == nz ) then ! Extrapolation
3664 : km1 = nz
3665 : k00 = nz-1
3666 : kp1 = nz
3667 : kp2 = nz
3668 : else
3669 : km1 = k-1
3670 : k00 = k
3671 : kp1 = k+1
3672 : kp2 = k+2
3673 : end if ! k == 1
3674 :
3675 : end if ! l_below
3676 :
3677 : ! Now perform the interpolation
3678 : var_value = mono_cubic_interp( z_interp, km1, k00, kp1, kp2, &
3679 : z_vals(km1), z_vals(k00), z_vals(kp1), z_vals(kp2), &
3680 : var(km1), var(k00), var(kp1), var(kp2) )
3681 :
3682 : return
3683 :
3684 : end function var_subgrid_interp
3685 :
3686 : !=============================================================================
3687 :
3688 : end module pdf_closure_module
|