Line data Source code
1 : !-------------------------------------------------------------------------------
2 : ! $Id$
3 : !===============================================================================
4 : module clip_explicit
5 :
6 : implicit none
7 :
8 : private
9 :
10 : public :: clip_covars_denom, &
11 : clip_covar, &
12 : clip_covar_level, &
13 : clip_variance, &
14 : clip_skewness, &
15 : clip_skewness_core
16 :
17 : ! Named constants to avoid string comparisons
18 : integer, parameter, public :: &
19 : clip_rtp2 = 1, & ! Named constant for rtp2 clipping
20 : clip_thlp2 = 2, & ! Named constant for thlp2 clipping
21 : clip_rtpthlp = 3, & ! Named constant for rtpthlp clipping
22 : clip_up2 = 5, & ! Named constant for up2 clipping
23 : clip_vp2 = 6, & ! Named constant for vp2 clipping
24 : ! clip_scalar = 7, & ! Named constant for scalar clipping
25 : clip_wprtp = 8, & ! Named constant for wprtp clipping
26 : clip_wpthlp = 9, & ! Named constant for wpthlp clipping
27 : clip_upwp = 10, & ! Named constant for upwp clipping
28 : clip_vpwp = 11, & ! Named constant for vpwp clipping
29 : clip_wp2 = 12, & ! Named constant for wp2 clipping
30 : clip_wpsclrp = 13, & ! Named constant for wp scalar clipping
31 : clip_sclrp2 = 14, & ! Named constant for sclrp2 clipping
32 : clip_sclrprtp = 15, & ! Named constant for sclrprtp clipping
33 : clip_sclrpthlp = 16, & ! Named constant for sclrpthlp clipping
34 : clip_wphydrometp = 17 ! Named constant for wphydrometp clipping
35 :
36 : contains
37 :
38 : !=============================================================================
39 17870112 : subroutine clip_covars_denom( nz, ngrdcol, sclr_dim, gr, dt, &
40 17870112 : rtp2, thlp2, up2, vp2, wp2, &
41 17870112 : sclrp2, wprtp_cl_num, wpthlp_cl_num, &
42 : wpsclrp_cl_num, upwp_cl_num, vpwp_cl_num, &
43 : l_predict_upwp_vpwp, &
44 : l_tke_aniso, &
45 : l_linearize_pbl_winds, &
46 : stats_metadata, &
47 17870112 : stats_zm, &
48 17870112 : wprtp, wpthlp, upwp, vpwp, wpsclrp, &
49 17870112 : upwp_pert, vpwp_pert )
50 :
51 : ! Description:
52 : ! Some of the covariances found in the CLUBB model code need to be clipped
53 : ! multiple times during each timestep to ensure that the correlation between
54 : ! the two relevant variables stays between -1 and 1 at all times during the
55 : ! model run. The covariances that need to be clipped multiple times are
56 : ! w'r_t', w'th_l', w'sclr', u'w', and v'w'. One of the times that each one
57 : ! of these covariances is clipped is immediately after each one is set.
58 : ! However, each covariance still needs to be clipped two more times during
59 : ! each timestep (once after advance_xp2_xpyp is called and once after
60 : ! advance_wp2_wp3 is called). This subroutine handles the times that the
61 : ! covariances are clipped away from the time that they are set. In other
62 : ! words, this subroutine clips the covariances after the denominator terms
63 : ! in the relevant correlation equation have been altered, ensuring that
64 : ! all correlations will remain between -1 and 1 at all times.
65 :
66 : ! References:
67 : ! None
68 : !-----------------------------------------------------------------------
69 :
70 : use grid_class, only: &
71 : grid ! Type
72 :
73 : use clubb_precision, only: &
74 : core_rknd ! Variable(s)
75 :
76 : use stats_type, only: &
77 : stats ! Type
78 :
79 : use stats_variables, only: &
80 : stats_metadata_type
81 :
82 : implicit none
83 :
84 : ! --------------------- Input Variables ---------------------
85 : integer, intent(in) :: &
86 : nz, &
87 : ngrdcol, &
88 : sclr_dim
89 :
90 : type (grid), target, intent(in) :: gr
91 :
92 : real( kind = core_rknd ), intent(in) :: &
93 : dt ! Timestep [s]
94 :
95 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
96 : rtp2, & ! r_t'^2 [(kg/kg)^2]
97 : thlp2, & ! theta_l'^2 [K^2]
98 : up2, & ! u'^2 [m^2/s^2]
99 : vp2, & ! v'^2 [m^2/s^2]
100 : wp2 ! w'^2 [m^2/s^2]
101 :
102 : real( kind = core_rknd ), dimension(ngrdcol,nz,sclr_dim), intent(in) :: &
103 : sclrp2 ! sclr'^2 [{units vary}^2]
104 :
105 : integer, intent(in) :: &
106 : wprtp_cl_num, &
107 : wpthlp_cl_num, &
108 : wpsclrp_cl_num, &
109 : upwp_cl_num, &
110 : vpwp_cl_num
111 :
112 : logical, intent(in) :: &
113 : l_predict_upwp_vpwp, & ! Flag to predict <u'w'> and <v'w'> along with <u> and <v> alongside
114 : ! the advancement of <rt>, <w'rt'>, <thl>, <wpthlp>, <sclr>, and
115 : ! <w'sclr'> in subroutine advance_xm_wpxp. Otherwise, <u'w'> and
116 : ! <v'w'> are still approximated by eddy diffusivity when <u> and <v>
117 : ! are advanced in subroutine advance_windm_edsclrm.
118 : l_tke_aniso, & ! For anisotropic turbulent kinetic energy, i.e. TKE = 1/2
119 : ! (u'^2 + v'^2 + w'^2)
120 : l_linearize_pbl_winds ! Flag (used by E3SM) to linearize PBL winds
121 :
122 : type (stats_metadata_type), intent(in) :: &
123 : stats_metadata
124 :
125 : ! --------------------- Input/Output Variables ---------------------
126 : type (stats), target, dimension(ngrdcol), intent(inout) :: &
127 : stats_zm
128 :
129 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(inout) :: &
130 : wprtp, & ! w'r_t' [(kg/kg) m/s]
131 : wpthlp, & ! w'theta_l' [K m/s]
132 : upwp, & ! u'w' [m^2/s^2]
133 : vpwp ! v'w' [m^2/s^2]
134 :
135 : real( kind = core_rknd ), dimension(ngrdcol,nz,sclr_dim), intent(inout) :: &
136 : wpsclrp ! w'sclr' [units m/s]
137 :
138 : ! Variables used to track perturbed version of winds.
139 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(inout) :: &
140 : upwp_pert, & ! perturbed <u'w'> [m^2/s^2]
141 : vpwp_pert ! perturbed <v'w'> [m^2/s^2]
142 :
143 : ! --------------------- Local Variables ---------------------
144 : logical :: &
145 : l_first_clip_ts, & ! First instance of clipping in a timestep.
146 : l_last_clip_ts ! Last instance of clipping in a timestep.
147 :
148 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
149 35740224 : wprtp_chnge, & ! Net change in w'r_t' due to clipping [(kg/kg) m/s]
150 35740224 : wpthlp_chnge, & ! Net change in w'th_l' due to clipping [K m/s]
151 35740224 : upwp_chnge, & ! Net change in u'w' due to clipping [m^2/s^2]
152 35740224 : vpwp_chnge ! Net change in v'w' due to clipping [m^2/s^2]
153 :
154 : real( kind = core_rknd ), dimension(ngrdcol,nz,sclr_dim) :: &
155 35740224 : wpsclrp_chnge ! Net change in w'sclr' due to clipping [{units vary}]
156 :
157 : integer :: sclr, i ! scalar array index.
158 :
159 : ! --------------------- Begin Code ---------------------
160 :
161 : !$acc enter data create( wprtp_chnge, wpthlp_chnge, upwp_chnge, vpwp_chnge )
162 : !$acc enter data if( sclr_dim > 0 ) create( wpsclrp_chnge )
163 :
164 : !!! Clipping for w'r_t'
165 : !
166 : ! Clipping w'r_t' at each vertical level, based on the
167 : ! correlation of w and r_t at each vertical level, such that:
168 : ! corr_(w,r_t) = w'r_t' / [ sqrt(w'^2) * sqrt(r_t'^2) ];
169 : ! -1 <= corr_(w,r_t) <= 1.
170 : !
171 : ! Since w'^2, r_t'^2, and w'r_t' are each advanced in different
172 : ! subroutines from each other in advance_clubb_core, clipping for w'r_t'
173 : ! is done three times during each timestep (once after each variable has
174 : ! been updated).
175 : !
176 : ! This subroutine handles the first and third instances of
177 : ! w'r_t' clipping.
178 : ! The first instance of w'r_t' clipping takes place after
179 : ! r_t'^2 is updated in advance_xp2_xpyp.
180 : ! The third instance of w'r_t' clipping takes place after
181 : ! w'^2 is updated in advance_wp2_wp3.
182 :
183 : ! Used within subroutine clip_covar.
184 17870112 : if ( wprtp_cl_num == 1 ) then
185 0 : l_first_clip_ts = .true.
186 0 : l_last_clip_ts = .false.
187 17870112 : elseif ( wprtp_cl_num == 2 ) then
188 8935056 : l_first_clip_ts = .false.
189 8935056 : l_last_clip_ts = .false.
190 8935056 : elseif ( wprtp_cl_num == 3 ) then
191 8935056 : l_first_clip_ts = .false.
192 8935056 : l_last_clip_ts = .true.
193 : endif
194 :
195 : ! Clip w'r_t'
196 : call clip_covar( nz, ngrdcol, gr, clip_wprtp, l_first_clip_ts, & ! intent(in)
197 : l_last_clip_ts, dt, wp2, rtp2, & ! intent(in)
198 : l_predict_upwp_vpwp, & ! intent(in)
199 : stats_metadata, & ! intent(in)
200 : stats_zm, & ! intent(inout)
201 17870112 : wprtp, wprtp_chnge ) ! intent(inout)
202 :
203 : !!! Clipping for w'th_l'
204 : !
205 : ! Clipping w'th_l' at each vertical level, based on the
206 : ! correlation of w and th_l at each vertical level, such that:
207 : ! corr_(w,th_l) = w'th_l' / [ sqrt(w'^2) * sqrt(th_l'^2) ];
208 : ! -1 <= corr_(w,th_l) <= 1.
209 : !
210 : ! Since w'^2, th_l'^2, and w'th_l' are each advanced in different
211 : ! subroutines from each other in advance_clubb_core, clipping for w'th_l'
212 : ! is done three times during each timestep (once after each variable has
213 : ! been updated).
214 : !
215 : ! This subroutine handles the first and third instances of
216 : ! w'th_l' clipping.
217 : ! The first instance of w'th_l' clipping takes place after
218 : ! th_l'^2 is updated in advance_xp2_xpyp.
219 : ! The third instance of w'th_l' clipping takes place after
220 : ! w'^2 is updated in advance_wp2_wp3.
221 :
222 : ! Used within subroutine clip_covar.
223 17870112 : if ( wpthlp_cl_num == 1 ) then
224 0 : l_first_clip_ts = .true.
225 0 : l_last_clip_ts = .false.
226 17870112 : elseif ( wpthlp_cl_num == 2 ) then
227 8935056 : l_first_clip_ts = .false.
228 8935056 : l_last_clip_ts = .false.
229 8935056 : elseif ( wpthlp_cl_num == 3 ) then
230 8935056 : l_first_clip_ts = .false.
231 8935056 : l_last_clip_ts = .true.
232 : endif
233 :
234 : ! Clip w'th_l'
235 : call clip_covar( nz, ngrdcol, gr, clip_wpthlp, l_first_clip_ts, & ! intent(in)
236 : l_last_clip_ts, dt, wp2, thlp2, & ! intent(in)
237 : l_predict_upwp_vpwp, & ! intent(in)
238 : stats_metadata, & ! intent(in)
239 : stats_zm, & ! intent(inout)
240 17870112 : wpthlp, wpthlp_chnge ) ! intent(inout)
241 :
242 : !!! Clipping for w'sclr'
243 : !
244 : ! Clipping w'sclr' at each vertical level, based on the
245 : ! correlation of w and sclr at each vertical level, such that:
246 : ! corr_(w,sclr) = w'sclr' / [ sqrt(w'^2) * sqrt(sclr'^2) ];
247 : ! -1 <= corr_(w,sclr) <= 1.
248 : !
249 : ! Since w'^2, sclr'^2, and w'sclr' are each advanced in different
250 : ! subroutines from each other in advance_clubb_core, clipping for w'sclr'
251 : ! is done three times during each timestep (once after each variable has
252 : ! been updated).
253 : !
254 : ! This subroutine handles the first and third instances of
255 : ! w'sclr' clipping.
256 : ! The first instance of w'sclr' clipping takes place after
257 : ! sclr'^2 is updated in advance_xp2_xpyp.
258 : ! The third instance of w'sclr' clipping takes place after
259 : ! w'^2 is updated in advance_wp2_wp3.
260 :
261 : ! Used within subroutine clip_covar.
262 17870112 : if ( wpsclrp_cl_num == 1 ) then
263 0 : l_first_clip_ts = .true.
264 0 : l_last_clip_ts = .false.
265 17870112 : elseif ( wpsclrp_cl_num == 2 ) then
266 8935056 : l_first_clip_ts = .false.
267 8935056 : l_last_clip_ts = .false.
268 8935056 : elseif ( wpsclrp_cl_num == 3 ) then
269 8935056 : l_first_clip_ts = .false.
270 8935056 : l_last_clip_ts = .true.
271 : endif
272 :
273 : ! Clip w'sclr'
274 17870112 : do sclr = 1, sclr_dim
275 : call clip_covar( nz, ngrdcol, gr, clip_wpsclrp, l_first_clip_ts, & ! intent(in)
276 : l_last_clip_ts, dt, wp2(:,:), sclrp2(:,:,sclr), & ! intent(in)
277 : l_predict_upwp_vpwp, & ! intent(in)
278 : stats_metadata, & ! intent(in)
279 : stats_zm, & ! intent(inout)
280 17870112 : wpsclrp(:,:,sclr), wpsclrp_chnge(:,:,sclr) ) ! intent(inout)
281 : enddo
282 :
283 :
284 : !!! Clipping for u'w'
285 : !
286 : ! Clipping u'w' at each vertical level, based on the
287 : ! correlation of u and w at each vertical level, such that:
288 : ! corr_(u,w) = u'w' / [ sqrt(u'^2) * sqrt(w'^2) ];
289 : ! -1 <= corr_(u,w) <= 1.
290 : !
291 : ! Since w'^2, u'^2, and u'w' are each advanced in different
292 : ! subroutines from each other in advance_clubb_core, clipping for u'w'
293 : ! is done three times during each timestep (once after each variable has
294 : ! been updated).
295 : !
296 : ! This subroutine handles the first and second instances of
297 : ! u'w' clipping.
298 : ! The first instance of u'w' clipping takes place after
299 : ! u'^2 is updated in advance_xp2_xpyp.
300 : ! The second instance of u'w' clipping takes place after
301 : ! w'^2 is updated in advance_wp2_wp3.
302 :
303 : ! Used within subroutine clip_covar.
304 17870112 : if ( upwp_cl_num == 1 ) then
305 0 : l_first_clip_ts = .true.
306 0 : l_last_clip_ts = .false.
307 17870112 : elseif ( upwp_cl_num == 2 ) then
308 8935056 : l_first_clip_ts = .false.
309 8935056 : l_last_clip_ts = .false.
310 8935056 : elseif ( upwp_cl_num == 3 ) then
311 8935056 : l_first_clip_ts = .false.
312 8935056 : l_last_clip_ts = .true.
313 : endif
314 :
315 : ! Clip u'w'
316 17870112 : if ( l_tke_aniso ) then
317 : call clip_covar( nz, ngrdcol, gr, clip_upwp, l_first_clip_ts, & ! intent(in)
318 : l_last_clip_ts, dt, wp2, up2, & ! intent(in)
319 : l_predict_upwp_vpwp, & ! intent(in)
320 : stats_metadata, & ! intent(in)
321 : stats_zm, & ! intent(inout)
322 17870112 : upwp, upwp_chnge ) ! intent(inout)
323 :
324 17870112 : if ( l_linearize_pbl_winds ) then
325 : call clip_covar( nz, ngrdcol, gr, clip_upwp, l_first_clip_ts, & ! intent(in)
326 : l_last_clip_ts, dt, wp2, up2, & ! intent(in)
327 : l_predict_upwp_vpwp, & ! intent(in)
328 : stats_metadata, & ! intent(in)
329 : stats_zm, & ! intent(inout)
330 0 : upwp_pert, upwp_chnge ) ! intent(inout)
331 : endif ! l_linearize_pbl_winds
332 : else
333 : ! In this case, up2 = wp2, and the variable `up2' does not interact
334 : call clip_covar( nz, ngrdcol, gr, clip_upwp, l_first_clip_ts, & ! intent(in)
335 : l_last_clip_ts, dt, wp2, wp2, & ! intent(in)
336 : l_predict_upwp_vpwp, & ! intent(in)
337 : stats_metadata, & ! intent(in)
338 : stats_zm, & ! intent(inout)
339 0 : upwp, upwp_chnge ) ! intent(inout)
340 :
341 0 : if ( l_linearize_pbl_winds ) then
342 : call clip_covar( nz, ngrdcol, gr, clip_upwp, l_first_clip_ts, & ! intent(in)
343 : l_last_clip_ts, dt, wp2, wp2, & ! intent(in)
344 : l_predict_upwp_vpwp, & ! intent(in)
345 : stats_metadata, & ! intent(in)
346 : stats_zm, & ! intent(inout)
347 0 : upwp_pert, upwp_chnge ) ! intent(inout)
348 : endif ! l_linearize_pbl_winds
349 : end if
350 :
351 :
352 :
353 : !!! Clipping for v'w'
354 : !
355 : ! Clipping v'w' at each vertical level, based on the
356 : ! correlation of v and w at each vertical level, such that:
357 : ! corr_(v,w) = v'w' / [ sqrt(v'^2) * sqrt(w'^2) ];
358 : ! -1 <= corr_(v,w) <= 1.
359 : !
360 : ! Since w'^2, v'^2, and v'w' are each advanced in different
361 : ! subroutines from each other in advance_clubb_core, clipping for v'w'
362 : ! is done three times during each timestep (once after each variable has
363 : ! been updated).
364 : !
365 : ! This subroutine handles the first and second instances of
366 : ! v'w' clipping.
367 : ! The first instance of v'w' clipping takes place after
368 : ! v'^2 is updated in advance_xp2_xpyp.
369 : ! The second instance of v'w' clipping takes place after
370 : ! w'^2 is updated in advance_wp2_wp3.
371 :
372 : ! Used within subroutine clip_covar.
373 17870112 : if ( vpwp_cl_num == 1 ) then
374 0 : l_first_clip_ts = .true.
375 0 : l_last_clip_ts = .false.
376 17870112 : elseif ( vpwp_cl_num == 2 ) then
377 8935056 : l_first_clip_ts = .false.
378 8935056 : l_last_clip_ts = .false.
379 8935056 : elseif ( vpwp_cl_num == 3 ) then
380 8935056 : l_first_clip_ts = .false.
381 8935056 : l_last_clip_ts = .true.
382 : endif
383 :
384 17870112 : if ( l_tke_aniso ) then
385 : call clip_covar( nz, ngrdcol, gr, clip_vpwp, l_first_clip_ts, & ! intent(in)
386 : l_last_clip_ts, dt, wp2, vp2, & ! intent(in)
387 : l_predict_upwp_vpwp, & ! intent(in)
388 : stats_metadata, & ! intent(in)
389 : stats_zm, & ! intent(inout)
390 17870112 : vpwp, vpwp_chnge ) ! intent(inout)
391 :
392 17870112 : if ( l_linearize_pbl_winds ) then
393 : call clip_covar( nz, ngrdcol, gr, clip_vpwp, l_first_clip_ts, & ! intent(in)
394 : l_last_clip_ts, dt, wp2, vp2, & ! intent(in)
395 : l_predict_upwp_vpwp, & ! intent(in)
396 : stats_metadata, & ! intent(in)
397 : stats_zm, & ! intent(inout)
398 0 : vpwp_pert, vpwp_chnge ) ! intent(inout)
399 : endif ! l_linearize_pbl_winds
400 : else
401 : ! In this case, vp2 = wp2, and the variable `vp2' does not interact
402 : call clip_covar( nz, ngrdcol, gr, clip_vpwp, l_first_clip_ts, & ! intent(in)
403 : l_last_clip_ts, dt, wp2, wp2, & ! intent(in)
404 : l_predict_upwp_vpwp, & ! intent(in)
405 : stats_metadata, & ! intent(in)
406 : stats_zm, & ! intent(inout)
407 0 : vpwp, vpwp_chnge ) ! intent(inout)
408 :
409 0 : if ( l_linearize_pbl_winds ) then
410 : call clip_covar( nz, ngrdcol, gr, clip_vpwp, l_first_clip_ts, & ! intent(in)
411 : l_last_clip_ts, dt, wp2, wp2, & ! intent(in)
412 : l_predict_upwp_vpwp, & ! intent(in)
413 : stats_metadata, & ! intent(in) stats_metadata, & ! intent(in)
414 : stats_zm, & ! intent(inout)
415 0 : vpwp_pert, vpwp_chnge ) ! intent(inout)
416 : endif ! l_linearize_pbl_winds
417 : end if
418 :
419 : !$acc exit data delete( wprtp_chnge, wpthlp_chnge, upwp_chnge, vpwp_chnge )
420 : !$acc exit data if( sclr_dim > 0 ) delete( wpsclrp_chnge )
421 :
422 17870112 : return
423 : end subroutine clip_covars_denom
424 :
425 : !=============================================================================
426 116155728 : subroutine clip_covar( nz, ngrdcol, gr, solve_type, l_first_clip_ts, &
427 116155728 : l_last_clip_ts, dt, xp2, yp2, &
428 : l_predict_upwp_vpwp, &
429 : stats_metadata, &
430 116155728 : stats_zm, &
431 116155728 : xpyp, xpyp_chnge )
432 :
433 : ! Description:
434 : ! Clipping the value of covariance x'y' based on the correlation between x
435 : ! and y.
436 : !
437 : ! The correlation between variables x and y is:
438 : !
439 : ! corr_(x,y) = x'y' / [ sqrt(x'^2) * sqrt(y'^2) ];
440 : !
441 : ! where x'^2 is the variance of x, y'^2 is the variance of y, and x'y' is
442 : ! the covariance of x and y.
443 : !
444 : ! The correlation of two variables must always have a value between -1
445 : ! and 1, such that:
446 : !
447 : ! -1 <= corr_(x,y) <= 1.
448 : !
449 : ! Therefore, there is an upper limit on x'y', such that:
450 : !
451 : ! x'y' <= [ sqrt(x'^2) * sqrt(y'^2) ];
452 : !
453 : ! and a lower limit on x'y', such that:
454 : !
455 : ! x'y' >= -[ sqrt(x'^2) * sqrt(y'^2) ].
456 : !
457 : ! The values of x'y', x'^2, and y'^2 are all found on momentum levels.
458 : !
459 : ! The value of x'y' may need to be clipped whenever x'y', x'^2, or y'^2 is
460 : ! updated.
461 : !
462 : ! The following covariances are found in the code:
463 : !
464 : ! w'r_t', w'th_l', w'sclr', (computed in advance_xm_wpxp);
465 : ! r_t'th_l', sclr'r_t', sclr'th_l', (computed in advance_xp2_xpyp);
466 : ! u'w', v'w', w'edsclr' (computed in advance_windm_edsclrm);
467 : ! and w'hm' (computed in setup_pdf_parameters).
468 :
469 : ! References:
470 : ! None
471 : !-----------------------------------------------------------------------
472 :
473 : use grid_class, only: &
474 : grid ! Type
475 :
476 : use constants_clubb, only: &
477 : max_mag_correlation, & ! Constant(s)
478 : max_mag_correlation_flux
479 :
480 : use clubb_precision, only: &
481 : core_rknd ! Variable(s)
482 :
483 : use stats_type_utilities, only: &
484 : stat_begin_update, & ! Procedure(s)
485 : stat_modify, &
486 : stat_end_update
487 :
488 : use stats_variables, only: &
489 : stats_metadata_type
490 :
491 : use stats_type, only: stats ! Type
492 :
493 : implicit none
494 :
495 : ! -------------------------- Input Variables --------------------------
496 : integer, intent(in) :: &
497 : nz, &
498 : ngrdcol
499 :
500 : type (grid), target, intent(in) :: gr
501 :
502 : integer, intent(in) :: &
503 : solve_type ! Variable being solved; used for STATS.
504 :
505 : logical, intent(in) :: &
506 : l_first_clip_ts, & ! First instance of clipping in a timestep.
507 : l_last_clip_ts ! Last instance of clipping in a timestep.
508 :
509 : real( kind = core_rknd ), intent(in) :: &
510 : dt ! Model timestep; used here for STATS [s]
511 :
512 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
513 : xp2, & ! Variance of x, x'^2 (momentum levels) [{x units}^2]
514 : yp2 ! Variance of y, y'^2 (momentum levels) [{y units}^2]
515 :
516 : logical, intent(in) :: &
517 : l_predict_upwp_vpwp ! Flag to predict <u'w'> and <v'w'> along with <u> and <v> alongside the
518 : ! advancement of <rt>, <w'rt'>, <thl>, <wpthlp>, <sclr>, and <w'sclr'> in
519 : ! subroutine advance_xm_wpxp. Otherwise, <u'w'> and <v'w'> are still
520 : ! approximated by eddy diffusivity when <u> and <v> are advanced in
521 : ! subroutine advance_windm_edsclrm.
522 :
523 : type (stats_metadata_type), intent(in) :: &
524 : stats_metadata
525 :
526 : ! -------------------------- InOut Variables --------------------------
527 : type (stats), target, dimension(ngrdcol), intent(inout) :: &
528 : stats_zm
529 :
530 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(inout) :: &
531 : xpyp ! Covariance of x and y, x'y' (momentum levels) [{x units}*{y units}]
532 :
533 : !-------------------------- Output Variable --------------------------
534 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
535 : xpyp_chnge ! Net change in x'y' due to clipping [{x units}*{y units}]
536 :
537 :
538 : ! -------------------------- Local Variables --------------------------
539 : real( kind = core_rknd ) :: &
540 : max_mag_corr, & ! Maximum magnitude of a correlation allowed
541 : xpyp_bound
542 :
543 : integer :: i, k ! Array index
544 :
545 : integer :: &
546 : ixpyp_cl
547 :
548 : ! -------------------------- Begin Code --------------------------
549 :
550 142960896 : select case ( solve_type )
551 : case ( clip_wprtp ) ! wprtp clipping budget term
552 26805168 : ixpyp_cl = stats_metadata%iwprtp_cl
553 : case ( clip_wpthlp ) ! wpthlp clipping budget term
554 26805168 : ixpyp_cl = stats_metadata%iwpthlp_cl
555 : case ( clip_rtpthlp ) ! rtpthlp clipping budget term
556 8935056 : ixpyp_cl = stats_metadata%irtpthlp_cl
557 : case ( clip_upwp ) ! upwp clipping budget term
558 26805168 : if ( l_predict_upwp_vpwp ) then
559 26805168 : ixpyp_cl = stats_metadata%iupwp_cl
560 : else
561 0 : ixpyp_cl = 0
562 : endif ! l_predict_upwp_vpwp
563 : case ( clip_vpwp ) ! vpwp clipping budget term
564 26805168 : if ( l_predict_upwp_vpwp ) then
565 26805168 : ixpyp_cl = stats_metadata%ivpwp_cl
566 : else
567 0 : ixpyp_cl = 0
568 : endif ! l_predict_upwp_vpwp
569 : case default ! scalars (or upwp/vpwp) are involved
570 116155728 : ixpyp_cl = 0
571 : end select
572 :
573 :
574 116155728 : if ( stats_metadata%l_stats_samp ) then
575 :
576 : !$acc update host( xpyp )
577 :
578 0 : if ( l_first_clip_ts ) then
579 0 : do i = 1, ngrdcol
580 0 : call stat_begin_update( nz, ixpyp_cl, xpyp(i,:) / dt, & ! intent(in)
581 0 : stats_zm(i) ) ! intent(inout)
582 : end do
583 : else
584 0 : do i = 1, ngrdcol
585 0 : call stat_modify( nz, ixpyp_cl, -xpyp(i,:) / dt, & ! intent(in)
586 0 : stats_zm(i) ) ! intent(inout)
587 : end do
588 : endif
589 : endif
590 :
591 : ! When clipping for wprtp or wpthlp, use the special value for
592 : ! max_mag_correlation_flux. For all other correlations, use
593 : ! max_mag_correlation.
594 : if ( ( solve_type == clip_wprtp ) .or. ( solve_type == clip_wpthlp ) ) then
595 : max_mag_corr = max_mag_correlation_flux
596 : else ! All other covariances
597 : max_mag_corr = max_mag_correlation
598 : endif ! solve_type
599 :
600 : ! The value of x'y' at the surface (or lower boundary) is a set value that
601 : ! is either specified or determined elsewhere in a surface subroutine. It
602 : ! is ensured elsewhere that the correlation between x and y at the surface
603 : ! (or lower boundary) is between -1 and 1. Thus, the covariance clipping
604 : ! code does not need to be invoked at the lower boundary. Likewise, the
605 : ! value of x'y' is set at the upper boundary, so the covariance clipping
606 : ! code does not need to be invoked at the upper boundary.
607 : ! Note that if clipping were applied at the lower boundary, momentum will
608 : ! not be conserved, therefore it should never be added.
609 : !$acc parallel loop gang vector collapse(2) default(present)
610 9757081152 : do k = 2, nz-1
611 >16109*10^7 : do i = 1, ngrdcol
612 >15134*10^7 : xpyp_bound = max_mag_corr * sqrt( xp2(i,k) * yp2(i,k) )
613 :
614 : ! Clipping for xpyp at an upper limit corresponding with a correlation
615 : ! between x and y of max_mag_corr.
616 >16098*10^7 : if ( xpyp(i,k) > xpyp_bound ) then
617 :
618 34070823 : xpyp_chnge(i,k) = xpyp_bound - xpyp(i,k)
619 34070823 : xpyp(i,k) = xpyp_bound
620 :
621 : ! Clipping for xpyp at a lower limit corresponding with a correlation
622 : ! between x and y of -max_mag_corr.
623 >15130*10^7 : else if ( xpyp(i,k) < -xpyp_bound ) then
624 :
625 30885521 : xpyp_chnge(i,k) = -xpyp_bound - xpyp(i,k)
626 30885521 : xpyp(i,k) = -xpyp_bound
627 :
628 : else
629 :
630 >15127*10^7 : xpyp_chnge(i,k) = 0.0_core_rknd
631 :
632 : end if
633 : end do
634 : end do
635 : !$acc end parallel loop
636 :
637 : ! Since there is no covariance clipping at the upper or lower boundaries,
638 : ! the change in x'y' due to covariance clipping at those levels is 0.
639 : !$acc parallel loop gang vector default(present)
640 1939530528 : do i = 1, ngrdcol
641 1823374800 : xpyp_chnge(i,1) = 0.0_core_rknd
642 1939530528 : xpyp_chnge(i,nz) = 0.0_core_rknd
643 : end do
644 : !$acc end parallel loop
645 :
646 116155728 : if ( stats_metadata%l_stats_samp ) then
647 :
648 : !$acc update host( xpyp )
649 :
650 0 : if ( l_last_clip_ts ) then
651 0 : do i = 1, ngrdcol
652 0 : call stat_end_update( nz, ixpyp_cl, xpyp(i,:) / dt, & ! intent(in)
653 0 : stats_zm(i) ) ! intent(inout)
654 : end do
655 : else
656 0 : do i = 1, ngrdcol
657 0 : call stat_modify( nz, ixpyp_cl, xpyp(i,:) / dt, & ! intent(in)
658 0 : stats_zm(i) ) ! intent(inout)
659 : end do
660 : endif
661 : endif
662 :
663 116155728 : return
664 :
665 : end subroutine clip_covar
666 :
667 : !=============================================================================
668 0 : subroutine clip_covar_level( solve_type, level, l_first_clip_ts, &
669 : l_last_clip_ts, dt, xp2, yp2, &
670 : l_predict_upwp_vpwp, &
671 : stats_metadata, &
672 : stats_zm, &
673 : xpyp, xpyp_chnge )
674 :
675 : ! Description:
676 : ! Clipping the value of covariance x'y' based on the correlation between x
677 : ! and y. This is all done at a single vertical level.
678 : !
679 : ! The correlation between variables x and y is:
680 : !
681 : ! corr_(x,y) = x'y' / [ sqrt(x'^2) * sqrt(y'^2) ];
682 : !
683 : ! where x'^2 is the variance of x, y'^2 is the variance of y, and x'y' is
684 : ! the covariance of x and y.
685 : !
686 : ! The correlation of two variables must always have a value between -1
687 : ! and 1, such that:
688 : !
689 : ! -1 <= corr_(x,y) <= 1.
690 : !
691 : ! Therefore, there is an upper limit on x'y', such that:
692 : !
693 : ! x'y' <= [ sqrt(x'^2) * sqrt(y'^2) ];
694 : !
695 : ! and a lower limit on x'y', such that:
696 : !
697 : ! x'y' >= -[ sqrt(x'^2) * sqrt(y'^2) ].
698 : !
699 : ! The values of x'y', x'^2, and y'^2 are all found on momentum levels.
700 : !
701 : ! The value of x'y' may need to be clipped whenever x'y', x'^2, or y'^2 is
702 : ! updated.
703 : !
704 : ! The following covariances are found in the code:
705 : !
706 : ! w'r_t', w'th_l', w'sclr', (computed in advance_xm_wpxp);
707 : ! r_t'th_l', sclr'r_t', sclr'th_l', (computed in advance_xp2_xpyp);
708 : ! u'w', v'w', w'edsclr' (computed in advance_windm_edsclrm);
709 : ! and w'hm' (computed in setup_pdf_parameters).
710 :
711 : ! References:
712 : ! None
713 : !-----------------------------------------------------------------------
714 :
715 : use constants_clubb, only: &
716 : max_mag_correlation, & ! Constant(s)
717 : max_mag_correlation_flux, &
718 : zero
719 :
720 : use clubb_precision, only: &
721 : core_rknd ! Variable(s)
722 :
723 : use stats_type_utilities, only: &
724 : stat_begin_update_pt, & ! Procedure(s)
725 : stat_modify_pt, &
726 : stat_end_update_pt
727 :
728 : use stats_variables, only: &
729 : stats_metadata_type
730 :
731 : use stats_type, only: stats ! Type
732 :
733 : implicit none
734 :
735 : type (stats), target, intent(inout) :: &
736 : stats_zm
737 :
738 : !------------------------- Input Variables -------------------------
739 : integer, intent(in) :: &
740 : solve_type, & ! Variable being solved; used for STATS
741 : level ! Vertical level index
742 :
743 : logical, intent(in) :: &
744 : l_first_clip_ts, & ! First instance of clipping in a timestep.
745 : l_last_clip_ts ! Last instance of clipping in a timestep.
746 :
747 : real( kind = core_rknd ), intent(in) :: &
748 : dt ! Model timestep; used here for STATS [s]
749 :
750 : real( kind = core_rknd ), intent(in) :: &
751 : xp2, & ! Variance of x, <x'^2> [{x units}^2]
752 : yp2 ! Variance of y, <y'^2> [{y units}^2]
753 :
754 : logical, intent(in) :: &
755 : l_predict_upwp_vpwp ! Flag to predict <u'w'> and <v'w'> along with <u> and <v> alongside the
756 : ! advancement of <rt>, <w'rt'>, <thl>, <wpthlp>, <sclr>, and <w'sclr'> in
757 : ! subroutine advance_xm_wpxp. Otherwise, <u'w'> and <v'w'> are still
758 : ! approximated by eddy diffusivity when <u> and <v> are advanced in
759 : ! subroutine advance_windm_edsclrm.
760 :
761 : type (stats_metadata_type), intent(in) :: &
762 : stats_metadata
763 :
764 : !------------------------- InOut Variable -------------------------
765 : real( kind = core_rknd ), intent(inout) :: &
766 : xpyp ! Covariance of x and y, <x'y'> [{x units}*{y units}]
767 :
768 : !------------------------- Output Variable -------------------------
769 : real( kind = core_rknd ), intent(out) :: &
770 : xpyp_chnge ! Net change in <x'y'> due to clipping [{x units}*{y units}]
771 :
772 :
773 : !------------------------- Local Variables -------------------------
774 : real( kind = core_rknd ) :: &
775 : max_mag_corr ! Maximum magnitude of a correlation allowed
776 :
777 : integer :: &
778 : ixpyp_cl ! Statistics index
779 :
780 : !------------------------- Begin Code -------------------------
781 :
782 0 : select case ( solve_type )
783 : case ( clip_wprtp ) ! wprtp clipping budget term
784 0 : ixpyp_cl = stats_metadata%iwprtp_cl
785 : case ( clip_wpthlp ) ! wpthlp clipping budget term
786 0 : ixpyp_cl = stats_metadata%iwpthlp_cl
787 : case ( clip_rtpthlp ) ! rtpthlp clipping budget term
788 0 : ixpyp_cl = stats_metadata%irtpthlp_cl
789 : case ( clip_upwp ) ! upwp clipping budget term
790 0 : if ( l_predict_upwp_vpwp ) then
791 0 : ixpyp_cl = stats_metadata%iupwp_cl
792 : else
793 0 : ixpyp_cl = 0
794 : endif ! l_predict_upwp_vpwp
795 : case ( clip_vpwp ) ! vpwp clipping budget term
796 0 : if ( l_predict_upwp_vpwp ) then
797 0 : ixpyp_cl = stats_metadata%ivpwp_cl
798 : else
799 0 : ixpyp_cl = 0
800 : endif ! l_predict_upwp_vpwp
801 : case default ! scalars (or upwp/vpwp) are involved
802 0 : ixpyp_cl = 0
803 : end select
804 :
805 :
806 0 : if ( stats_metadata%l_stats_samp ) then
807 0 : if ( l_first_clip_ts ) then
808 : call stat_begin_update_pt( ixpyp_cl, level, xpyp / dt, & ! intent(in)
809 0 : stats_zm ) ! intent(inout)
810 : else
811 : call stat_modify_pt( ixpyp_cl, level, -xpyp / dt, & ! intent(in)
812 0 : stats_zm ) ! intent(inout)
813 : endif
814 : endif
815 :
816 : ! When clipping for wprtp or wpthlp, use the special value for
817 : ! max_mag_correlation_flux. For all other correlations, use
818 : ! max_mag_correlation.
819 : if ( ( solve_type == clip_wprtp ) .or. ( solve_type == clip_wpthlp ) ) then
820 : max_mag_corr = max_mag_correlation_flux
821 : else ! All other covariances
822 : max_mag_corr = max_mag_correlation
823 : endif ! solve_type
824 :
825 : ! The value of x'y' at the surface (or lower boundary) is a set value that
826 : ! is either specified or determined elsewhere in a surface subroutine. It
827 : ! is ensured elsewhere that the correlation between x and y at the surface
828 : ! (or lower boundary) is between -1 and 1. Thus, the covariance clipping
829 : ! code does not need to be invoked at the lower boundary. Likewise, the
830 : ! value of x'y' is set at the upper boundary, so the covariance clipping
831 : ! code does not need to be invoked at the upper boundary.
832 : ! Note that if clipping were applied at the lower boundary, momentum will
833 : ! not be conserved, therefore it should never be added.
834 :
835 : ! Clipping for xpyp at an upper limit corresponding with a correlation
836 : ! between x and y of max_mag_corr.
837 0 : if ( xpyp > max_mag_corr * sqrt( xp2 * yp2 ) ) then
838 :
839 0 : xpyp_chnge = max_mag_corr * sqrt( xp2 * yp2 ) - xpyp
840 :
841 0 : xpyp = max_mag_corr * sqrt( xp2 * yp2 )
842 :
843 : ! Clipping for xpyp at a lower limit corresponding with a correlation
844 : ! between x and y of -max_mag_corr.
845 0 : elseif ( xpyp < -max_mag_corr * sqrt( xp2 * yp2 ) ) then
846 :
847 0 : xpyp_chnge = -max_mag_corr * sqrt( xp2 * yp2 ) - xpyp
848 :
849 0 : xpyp = -max_mag_corr * sqrt( xp2 * yp2 )
850 :
851 : else
852 :
853 0 : xpyp_chnge = zero
854 :
855 : endif
856 :
857 0 : if ( stats_metadata%l_stats_samp ) then
858 0 : if ( l_last_clip_ts ) then
859 : call stat_end_update_pt( ixpyp_cl, level, xpyp / dt, & ! intent(in)
860 0 : stats_zm ) ! intent(inout)
861 : else
862 : call stat_modify_pt( ixpyp_cl, level, xpyp / dt, & ! intent(in)
863 0 : stats_zm ) ! intent(inout)
864 : endif
865 : endif
866 :
867 :
868 0 : return
869 : end subroutine clip_covar_level
870 :
871 : !=============================================================================
872 44675280 : subroutine clip_variance( nz, ngrdcol, gr, solve_type, dt, threshold, &
873 : stats_metadata, &
874 44675280 : stats_zm, &
875 44675280 : xp2 )
876 :
877 : ! Description:
878 : ! Clipping the value of variance x'^2 based on a minimum threshold value.
879 : ! The threshold value must be greater than or equal to 0.
880 : !
881 : ! The values of x'^2 are found on the momentum levels.
882 : !
883 : ! The following variances are found in the code:
884 : !
885 : ! r_t'^2, th_l'^2, u'^2, v'^2, sclr'^2, (computed in advance_xp2_xpyp);
886 : ! w'^2 (computed in advance_wp2_wp3).
887 :
888 : ! References:
889 : ! None
890 : !-----------------------------------------------------------------------
891 :
892 : use grid_class, only: &
893 : grid ! Type
894 :
895 : use clubb_precision, only: &
896 : core_rknd ! Variable(s)
897 :
898 : use stats_type_utilities, only: &
899 : stat_begin_update, & ! Procedure(s)
900 : stat_end_update
901 :
902 : use stats_variables, only: &
903 : stats_metadata_type
904 :
905 : use stats_type, only: stats ! Type
906 :
907 : implicit none
908 :
909 : ! -------------------- Input Variables --------------------
910 : integer, intent(in) :: &
911 : nz, &
912 : ngrdcol
913 :
914 : type (grid), target, intent(in) :: gr
915 :
916 : integer, intent(in) :: &
917 : solve_type ! Variable being solved; used for STATS.
918 :
919 : real( kind = core_rknd ), intent(in) :: &
920 : dt ! Model timestep; used here for STATS [s]
921 :
922 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
923 : threshold ! Minimum value of x'^2 [{x units}^2]
924 :
925 : type (stats_metadata_type), intent(in) :: &
926 : stats_metadata
927 :
928 : ! -------------------- InOut Variables --------------------
929 : type (stats), target, dimension(ngrdcol), intent(inout) :: &
930 : stats_zm
931 :
932 : ! -------------------- Output Variable --------------------
933 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(inout) :: &
934 : xp2 ! Variance of x, x'^2 (momentum levels) [{x units}^2]
935 :
936 : ! -------------------- Local Variables --------------------
937 : integer :: i, k ! Array index
938 :
939 : integer :: &
940 : ixp2_cl
941 :
942 : ! -------------------- Begin Code --------------------
943 :
944 : !$acc data copyin( threshold ) &
945 : !$acc copy( xp2 )
946 :
947 53610336 : select case ( solve_type )
948 : case ( clip_wp2 ) ! wp2 clipping budget term
949 8935056 : ixp2_cl = stats_metadata%iwp2_cl
950 : case ( clip_rtp2 ) ! rtp2 clipping budget term
951 8935056 : ixp2_cl = stats_metadata%irtp2_cl
952 : case ( clip_thlp2 ) ! thlp2 clipping budget term
953 8935056 : ixp2_cl = stats_metadata%ithlp2_cl
954 : case ( clip_up2 ) ! up2 clipping budget term
955 8935056 : ixp2_cl = stats_metadata%iup2_cl
956 : case ( clip_vp2 ) ! vp2 clipping budget term
957 8935056 : ixp2_cl = stats_metadata%ivp2_cl
958 : case default ! scalars are involved
959 44675280 : ixp2_cl = 0
960 : end select
961 :
962 :
963 44675280 : if ( stats_metadata%l_stats_samp ) then
964 : !$acc update host( xp2 )
965 0 : do i = 1, ngrdcol
966 0 : call stat_begin_update( nz, ixp2_cl, xp2(i,:) / dt, & ! intent(in)
967 0 : stats_zm(i) ) ! intent(inout)
968 : end do
969 : end if
970 :
971 : ! Limit the value of x'^2 at threshold.
972 : ! The value of x'^2 at the surface (or lower boundary) is a set value that
973 : ! is determined elsewhere in a surface subroutine. Thus, the variance
974 : ! clipping code does not need to be invoked at the lower boundary.
975 : ! Likewise, the value of x'^2 is set at the upper boundary, so the variance
976 : ! clipping code does not need to be invoked at the upper boundary.
977 : !
978 : ! charlass on 09/11/2013: I changed the clipping so that also the surface
979 : ! level is clipped. I did this because we discovered that there are slightly
980 : ! negative values in thlp2(1) and rtp2(1) when running quarter_ss case with
981 : ! WRF-CLUBB (see wrf:ticket:51#comment:33)
982 : !$acc parallel loop gang vector collapse(2) default(present)
983 3797398800 : do k = 1, nz-1, 1
984 62706430800 : do i = 1, ngrdcol
985 62661755520 : if ( xp2(i,k) < threshold(i,k) ) then
986 57712676 : xp2(i,k) = threshold(i,k)
987 : end if
988 : end do
989 : end do
990 : !$acc end parallel loop
991 :
992 44675280 : if ( stats_metadata%l_stats_samp ) then
993 : !$acc update host( xp2 )
994 0 : do i = 1, ngrdcol
995 0 : call stat_end_update( nz, ixp2_cl, xp2(i,:) / dt, & ! intent(in)
996 0 : stats_zm(i) ) ! intent(inout)
997 : end do
998 : end if
999 :
1000 : !$acc end data
1001 :
1002 44675280 : return
1003 :
1004 : end subroutine clip_variance
1005 :
1006 : !=============================================================================
1007 8935056 : subroutine clip_skewness( nz, ngrdcol, gr, dt, sfc_elevation, & ! intent(in)
1008 8935056 : Skw_max_mag, wp2_zt, & ! intent(in)
1009 : l_use_wp3_lim_with_smth_Heaviside, & ! intent(in)
1010 : stats_metadata, & ! intent(in)
1011 8935056 : stats_zt, & ! intent(inout)
1012 8935056 : wp3 ) ! intent(out)
1013 :
1014 : ! Description:
1015 : ! Clipping the value of w'^3 based on the skewness of w, Sk_w.
1016 : !
1017 : ! Aditionally, to prevent possible crashes due to wp3 growing too large,
1018 : ! abs(wp3) will be clipped to 100.
1019 : !
1020 : ! The skewness of w is:
1021 : !
1022 : ! Sk_w = w'^3 / (w'^2)^(3/2).
1023 : !
1024 : ! The value of Sk_w is limited to a range between an upper limit and a lower
1025 : ! limit. The values of the limits depend on whether the level altitude is
1026 : ! within 100 meters of the surface.
1027 : !
1028 : ! For altitudes less than or equal to 100 meters above ground level (AGL):
1029 : !
1030 : ! -0.2_core_rknd*sqrt(2) <= Sk_w <= 0.2_core_rknd*sqrt(2);
1031 : !
1032 : ! while for all altitudes greater than 100 meters AGL:
1033 : !
1034 : ! -4.5_core_rknd <= Sk_w <= 4.5_core_rknd.
1035 : !
1036 : ! Therefore, there is an upper limit on w'^3, such that:
1037 : !
1038 : ! w'^3 <= threshold_magnitude * (w'^2)^(3/2);
1039 : !
1040 : ! and a lower limit on w'^3, such that:
1041 : !
1042 : ! w'^3 >= -threshold_magnitude * (w'^2)^(3/2).
1043 : !
1044 : ! The values of w'^3 are found on the thermodynamic levels, while the values
1045 : ! of w'^2 are found on the momentum levels. Therefore, the values of w'^2
1046 : ! are interpolated to the thermodynamic levels before being used to
1047 : ! calculate the upper and lower limits for w'^3.
1048 :
1049 : ! References:
1050 : ! None
1051 : !-----------------------------------------------------------------------
1052 :
1053 : use grid_class, only: &
1054 : grid ! Type
1055 :
1056 : use clubb_precision, only: &
1057 : core_rknd ! Variable(s)
1058 :
1059 : use stats_type_utilities, only: &
1060 : stat_begin_update, & ! Procedure(s)
1061 : stat_end_update
1062 :
1063 : use stats_variables, only: &
1064 : stats_metadata_type
1065 :
1066 : use stats_type, only: stats ! Type
1067 :
1068 : implicit none
1069 :
1070 : ! ----------------------- Input Variables -----------------------
1071 : integer, intent(in) :: &
1072 : nz, &
1073 : ngrdcol
1074 :
1075 : type (grid), target, intent(in) :: gr
1076 :
1077 : real( kind = core_rknd ), intent(in) :: &
1078 : dt ! Model timestep; used here for STATS [s]
1079 :
1080 : real( kind = core_rknd ), dimension(ngrdcol), intent(in) :: &
1081 : sfc_elevation, & ! Elevation of ground level [m AMSL]
1082 : Skw_max_mag ! Maximum allowable magnitude of Skewness [-]
1083 :
1084 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
1085 : wp2_zt ! w'^2 interpolated to thermodyamic levels [m^2/s^2]
1086 :
1087 : ! Flag to activate modifications on wp3 limiters for convergence test
1088 : ! (use smooth Heaviside 'Preskin' function in the calculation of
1089 : ! clip_skewness for wp3)
1090 : logical, intent(in):: &
1091 : l_use_wp3_lim_with_smth_Heaviside
1092 :
1093 : type (stats_metadata_type), intent(in) :: &
1094 : stats_metadata
1095 :
1096 : ! ----------------------- Input/Output Variables -----------------------
1097 : type (stats), target, dimension(ngrdcol), intent(inout) :: &
1098 : stats_zt
1099 :
1100 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(inout) :: &
1101 : wp3 ! w'^3 (thermodynamic levels) [m^3/s^3]
1102 :
1103 : ! ----------------------- Local Variables -----------------------
1104 : integer :: i
1105 :
1106 : ! ----------------------- Begin Code -----------------------
1107 :
1108 : !$acc data copyin( gr, gr%zt, &
1109 : !$acc sfc_elevation, wp2_zt ) &
1110 : !$acc copy( wp3 )
1111 :
1112 8935056 : if ( stats_metadata%l_stats_samp ) then
1113 :
1114 : !$acc update host( wp3 )
1115 :
1116 0 : do i = 1, ngrdcol
1117 0 : call stat_begin_update( nz, stats_metadata%iwp3_cl, wp3(i,:) / dt, & ! intent(in)
1118 0 : stats_zt(i) ) ! intent(inout)
1119 : end do
1120 : end if
1121 :
1122 : call clip_skewness_core( nz, ngrdcol, gr, sfc_elevation, & ! intent(in)
1123 : Skw_max_mag, wp2_zt, & ! intent(in)
1124 : l_use_wp3_lim_with_smth_Heaviside, & ! intent(in)
1125 8935056 : wp3 ) ! intent(inout)
1126 :
1127 8935056 : if ( stats_metadata%l_stats_samp ) then
1128 :
1129 : !$acc update host( wp3 )
1130 :
1131 0 : do i = 1, ngrdcol
1132 0 : call stat_end_update( nz, stats_metadata%iwp3_cl, wp3(i,:) / dt, & ! intent(in)
1133 0 : stats_zt(i) ) ! intent(inout)
1134 : end do
1135 : end if
1136 :
1137 : !$acc end data
1138 :
1139 8935056 : return
1140 :
1141 : end subroutine clip_skewness
1142 :
1143 : !=============================================================================
1144 8935056 : subroutine clip_skewness_core( nz, ngrdcol, gr, sfc_elevation, &
1145 8935056 : Skw_max_mag, wp2_zt, &
1146 : l_use_wp3_lim_with_smth_Heaviside, &
1147 8935056 : wp3 )
1148 :
1149 : use grid_class, only: &
1150 : grid ! Type
1151 :
1152 : use clubb_precision, only: &
1153 : core_rknd ! Variable(s)
1154 :
1155 : use advance_helper_module, only: &
1156 : smooth_heaviside_peskin
1157 :
1158 : implicit none
1159 :
1160 : !----------------------- Input Variables -----------------------
1161 : integer, intent(in) :: &
1162 : nz, &
1163 : ngrdcol
1164 :
1165 : type (grid), target, intent(in) :: gr
1166 :
1167 : real( kind = core_rknd ), dimension(ngrdcol), intent(in) :: &
1168 : sfc_elevation, & ! Elevation of ground level [m AMSL]
1169 : Skw_max_mag ! Maximum allowable magnitude of Skewness [-]
1170 :
1171 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
1172 : wp2_zt ! w'^2 interpolated to thermodyamic levels [m^2/s^2]
1173 :
1174 : ! Flag to activate modifications on wp3 limiters for convergence test
1175 : ! (use smooth Heaviside 'Preskin' function in the calculation of clip_skewness for wp3)
1176 : logical, intent(in):: &
1177 : l_use_wp3_lim_with_smth_Heaviside
1178 :
1179 : !----------------------- Input/Output Variables -----------------------
1180 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(inout) :: &
1181 : wp3 ! w'^3 (thermodynamic levels) [m^3/s^3]
1182 :
1183 : !----------------------- Local Variables -----------------------
1184 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
1185 17870112 : wp2_zt_cubed, & ! Variance of vertical velocity cubed (w^2_{zt}^3) [m^6/s^6]
1186 17870112 : wp3_lim_sqd ! Keeps absolute value of Sk_w from becoming > limit [m^6/s^6]
1187 :
1188 : integer :: i, k ! Vertical array index.
1189 :
1190 : real( kind = core_rknd ), parameter :: &
1191 : wp3_max = 100._core_rknd ! Threshold for wp3 [m^3/s^3]
1192 :
1193 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
1194 17870112 : zagl_thresh, & ! temporatory array
1195 17870112 : H_zagl ! Heaviside function for clippings of wp3_lim_sqd
1196 :
1197 : !----------------------- Begin Code-----------------------
1198 :
1199 : !$acc enter data create( wp2_zt_cubed, wp3_lim_sqd, zagl_thresh, H_zagl )
1200 :
1201 : ! Compute the upper and lower limits of w'^3 at every level,
1202 : ! based on the skewness of w, Sk_w, such that:
1203 : ! Sk_w = w'^3 / (w'^2)^(3/2);
1204 : ! -4.5 <= Sk_w <= 4.5;
1205 : ! or, if the level altitude is within 100 meters of the surface,
1206 : ! -0.2*sqrt(2) <= Sk_w <= 0.2*sqrt(2).
1207 :
1208 : ! The normal magnitude limit of skewness of w in the CLUBB code is 4.5.
1209 : ! However, according to Andre et al. (1976b & 1978), wp3 should not exceed
1210 : ! [2*(wp2^3)]^(1/2) at any level. However, this term should be multiplied
1211 : ! by 0.2 close to the surface to include surface effects. There already is
1212 : ! a wp3 clipping term in place for all other altitudes, but this term will
1213 : ! be included for the surface layer only. Therefore, the lowest level wp3
1214 : ! should not exceed 0.2 * sqrt(2) * wp2^(3/2). Brian Griffin. 12/18/05.
1215 :
1216 : ! To lower compute time, we squared both sides of the equation and compute
1217 : ! wp2^3 only once. -dschanen 9 Oct 2008
1218 : !$acc parallel loop gang vector collapse(2) default(present)
1219 768414816 : do k = 1, nz
1220 12690480816 : do i = 1, ngrdcol
1221 12681545760 : wp2_zt_cubed(i,k) = wp2_zt(i,k)**3
1222 : end do
1223 : end do
1224 : !$acc end parallel loop
1225 :
1226 8935056 : if ( l_use_wp3_lim_with_smth_Heaviside ) then
1227 :
1228 : !implement a smoothed Heaviside function to avoid discontinuities
1229 : !$acc parallel loop gang vector collapse(2) default(present)
1230 0 : do k = 1, nz
1231 0 : do i = 1, ngrdcol
1232 0 : zagl_thresh(i,k) = ( gr%zt(i,k) - sfc_elevation(i) ) / 100.0_core_rknd
1233 0 : zagl_thresh(i,k) = zagl_thresh(i,k) - 1.0_core_rknd
1234 : end do
1235 : end do
1236 : !$acc end parallel loop
1237 :
1238 0 : H_zagl(:,:) = smooth_heaviside_peskin(nz, ngrdcol, zagl_thresh(:,:), 0.6_core_rknd)
1239 :
1240 : !$acc parallel loop gang vector collapse(2) default(present)
1241 0 : do k = 1, nz
1242 0 : do i = 1, ngrdcol
1243 0 : wp3_lim_sqd(i,k) = wp2_zt_cubed(i,k) &
1244 : * ( H_zagl(i,k) * Skw_max_mag(i)**2 &
1245 : + (1.0_core_rknd - H_zagl(i,k)) &
1246 0 : * 0.0021_core_rknd *Skw_max_mag(i)**2 )
1247 : end do
1248 : end do
1249 : !$acc end parallel loop
1250 :
1251 : else ! default method
1252 :
1253 : !$acc parallel loop gang vector collapse(2) default(present)
1254 768414816 : do k = 1, nz
1255 12690480816 : do i = 1, ngrdcol
1256 12681545760 : if ( gr%zt(i,k) - sfc_elevation(i) <= 100.0_core_rknd ) then ! Clip for 100 m. AGL.
1257 : !wp3_upper_lim(k) = 0.2_core_rknd * sqrt_2 * wp2_zt(k)**(3.0_core_rknd/2.0_core_rknd)
1258 : !wp3_lower_lim(k) = -0.2_core_rknd * sqrt_2 * wp2_zt(k)**(3.0_core_rknd/2.0_core_rknd)
1259 420778800 : wp3_lim_sqd(i,k) = 0.0021_core_rknd * Skw_max_mag(i)**2 * wp2_zt_cubed(i,k)
1260 : else ! Clip skewness consistently with a.
1261 : !wp3_upper_lim(k) = 4.5_core_rknd * wp2_zt(k)**(3.0_core_rknd/2.0_core_rknd)
1262 : !wp3_lower_lim(k) = -4.5_core_rknd * wp2_zt(k)**(3.0_core_rknd/2.0_core_rknd)
1263 11501287200 : wp3_lim_sqd(i,k) = Skw_max_mag(i)**2 * wp2_zt_cubed(i,k) ! Skw_max_mag = 4.5_core_rknd^2
1264 : endif
1265 : end do
1266 : end do
1267 : !$acc end parallel loop
1268 :
1269 : end if
1270 :
1271 : ! Clipping for w'^3 at an upper and lower limit corresponding with
1272 : ! the appropriate value of Sk_w.
1273 : !$acc parallel loop gang vector collapse(2) default(present)
1274 768414816 : do k = 1, nz
1275 12690480816 : do i = 1, ngrdcol
1276 : ! Set the magnitude to the wp3 limit and apply the sign of the current wp3
1277 12681545760 : if ( wp3(i,k)**2 > wp3_lim_sqd(i,k) ) then
1278 201653678 : wp3(i,k) = sign( sqrt( wp3_lim_sqd(i,k) ), wp3(i,k) )
1279 : end if
1280 : end do
1281 : end do
1282 : !$acc end parallel loop
1283 :
1284 : ! Clipping abs(wp3) to 100. This keeps wp3 from growing too large in some
1285 : ! deep convective cases, which helps prevent these cases from blowing up.
1286 : !$acc parallel loop gang vector collapse(2) default(present)
1287 768414816 : do k = 1, nz
1288 12690480816 : do i = 1, ngrdcol
1289 12681545760 : if ( abs(wp3(i,k)) > wp3_max ) then
1290 4 : wp3(i,k) = sign( wp3_max, wp3(i,k) ) ! Known magic number
1291 : end if
1292 : end do
1293 : end do
1294 : !$acc end parallel loop
1295 :
1296 : !$acc exit data delete( wp2_zt_cubed, wp3_lim_sqd, zagl_thresh, H_zagl )
1297 :
1298 8935056 : end subroutine clip_skewness_core
1299 :
1300 : !===============================================================================
1301 :
1302 : end module clip_explicit
|