Line data Source code
1 : !-----------------------------------------------------------------------
2 : ! $Id: mixing_length.F90 8664 2018-05-10 20:21:35Z huebler@uwm.edu $
3 : !===============================================================================
4 : module mixing_length
5 :
6 : implicit none
7 :
8 : private ! Default Scope
9 :
10 : public :: calc_Lscale_directly, &
11 : diagnose_Lscale_from_tau
12 :
13 : contains
14 :
15 : !=============================================================================
16 8935056 : subroutine compute_mixing_length( nz, ngrdcol, gr, thvm, thlm, &
17 8935056 : rtm, em, Lscale_max, p_in_Pa, &
18 8935056 : exner, thv_ds, mu, lmin, &
19 : saturation_formula, &
20 : l_implemented, &
21 : stats_metadata, &
22 8935056 : Lscale, Lscale_up, Lscale_down )
23 :
24 : ! Description:
25 : ! Larson's 5th moist, nonlocal length scale
26 : !
27 : ! References:
28 : ! Section 3b ( /Eddy length formulation/ ) of
29 : ! ``A PDF-Based Model for Boundary Layer Clouds. Part I:
30 : ! Method and Model Description'' Golaz, et al. (2002)
31 : ! JAS, Vol. 59, pp. 3540--3551.
32 : !
33 : ! Notes:
34 : !
35 : ! The equation for the rate of change of theta_l and r_t of the parcel with
36 : ! respect to height, due to entrainment, is:
37 : !
38 : ! d(thl_par)/dz = - mu * ( thl_parcel - thl_environment );
39 : !
40 : ! d(rt_par)/dz = - mu * ( rt_parcel - rt_environment );
41 : !
42 : ! where mu is the entrainment rate,
43 : ! such that:
44 : !
45 : ! mu = (1/m)*(dm/dz);
46 : !
47 : ! where m is the mass of the parcel. The value of mu is set to be a
48 : ! constant.
49 : !
50 : ! The differential equations are solved for given the boundary condition
51 : ! and given the fact that the value of thl_environment and rt_environment
52 : ! are treated as changing linearly for a parcel of air from one grid level
53 : ! to the next.
54 : !
55 : ! For the special case where entrainment rate, mu, is set to 0,
56 : ! thl_parcel and rt_parcel remain constant
57 : !
58 : !
59 : ! The equation for Lscale_up is:
60 : !
61 : ! INT(z_i:z_i+Lscale_up) g * ( thv_par - thvm ) / thvm dz = -em(z_i);
62 : !
63 : ! and for Lscale_down
64 : !
65 : ! INT(z_i-Lscale_down:z_i) g * ( thv_par - thvm ) / thvm dz = em(z_i);
66 : !
67 : ! where thv_par is theta_v of the parcel, thvm is the mean
68 : ! environmental value of theta_v, z_i is the altitude that the parcel
69 : ! started from, and em is the mean value of TKE at
70 : ! altitude z_i (which gives the parcel its initial boost)
71 : !
72 : ! The increment of CAPE (convective air potential energy) for any two
73 : ! successive vertical levels is:
74 : !
75 : ! Upwards:
76 : ! CAPE_incr = INT(z_0:z_1) g * ( thv_par - thvm ) / thvm dz
77 : !
78 : ! Downwards:
79 : ! CAPE_incr = INT(z_(-1):z_0) g * ( thv_par - thvm ) / thvm dz
80 : !
81 : ! Thus, the derivative of CAPE with respect to height is:
82 : !
83 : ! dCAPE/dz = g * ( thv_par - thvm ) / thvm.
84 : !
85 : ! A purely trapezoidal rule is used between levels, and is considered
86 : ! to vary linearly at all altitudes. Thus, dCAPE/dz is considered to be
87 : ! of the form: A * (z-zo) + dCAPE/dz|_(z_0),
88 : ! where A = ( dCAPE/dz|_(z_1) - dCAPE/dz|_(z_0) ) / ( z_1 - z_0 )
89 : !
90 : ! The integral is evaluated to find the CAPE increment between two
91 : ! successive vertical levels. The result either adds to or depletes
92 : ! from the total amount of energy that keeps the parcel ascending/descending.
93 : !
94 : !
95 : ! IMPORTANT NOTE:
96 : ! This subroutine has been optimized by adding precalculations, rearranging
97 : ! equations to avoid divides, and modifying the algorithm entirely.
98 : ! -Gunther Huebler
99 : !
100 : ! The algorithm previously used looped over every grid level, following a
101 : ! a parcel up from its initial grid level to its max. The very nature of this
102 : ! algorithm is an N^2
103 : !--------------------------------------------------------------------------------
104 :
105 : ! mu = (1/M) dM/dz > 0. mu=0 for no entrainment.
106 : ! Siebesma recommends mu=2e-3, although most schemes use mu=1e-4
107 : ! When mu was fixed, we used the value mu = 6.e-4
108 :
109 : use constants_clubb, only: & ! Variable(s)
110 : Cp, & ! Dry air specific heat at constant pressure [J/kg/K]
111 : Rd, & ! Dry air gas constant [J/kg/K]
112 : ep, & ! Rd / Rv [-]
113 : ep1, & ! (1-ep)/ep [-]
114 : ep2, & ! 1/ep [-]
115 : Lv, & ! Latent heat of vaporiztion [J/kg/K]
116 : grav, & ! Gravitational acceleration [m/s^2]
117 : fstderr, &
118 : zero_threshold, &
119 : eps, &
120 : one_half, &
121 : one, &
122 : two, &
123 : zero
124 :
125 : use grid_class, only: &
126 : grid, & ! Type
127 : zm2zt ! Procedure(s)
128 :
129 : use numerical_check, only: &
130 : length_check ! Procedure(s)
131 :
132 : use clubb_precision, only: &
133 : core_rknd ! Variable(s)
134 :
135 : use error_code, only: &
136 : clubb_at_least_debug_level, & ! Procedure
137 : err_code, & ! Error Indicator
138 : clubb_fatal_error ! Constant
139 :
140 : use saturation, only: &
141 : sat_mixrat_liq ! Procedure(s)
142 :
143 : use stats_variables, only: &
144 : stats_metadata_type
145 :
146 : implicit none
147 :
148 : ! Constant Parameters
149 : real( kind = core_rknd ), parameter :: &
150 : zlmin = 0.1_core_rknd, & ! Minimum value for Lscale [m]
151 : Lscale_sfclyr_depth = 500._core_rknd ! [m]
152 :
153 : !--------------------------------- Input Variables ---------------------------------
154 : integer, intent(in) :: &
155 : nz, &
156 : ngrdcol
157 :
158 : type (grid), target, intent(in) :: &
159 : gr
160 :
161 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
162 : thvm, & ! Virtual potential temp. on themodynamic level [K]
163 : thlm, & ! Liquid potential temp. on themodynamic level [K]
164 : rtm, & ! Total water mixing ratio on themodynamic level [kg/kg]
165 : em, & ! em = 3/2 * w'^2; on momentum level [m^2/s^2]
166 : exner, & ! Exner function on thermodynamic level [-]
167 : p_in_Pa, & ! Pressure on thermodynamic level [Pa]
168 : thv_ds ! Dry, base-state theta_v on thermodynamic level [K]
169 : ! Note: thv_ds used as a reference theta_l here
170 :
171 : real( kind = core_rknd ), dimension(ngrdcol), intent(in) :: &
172 : Lscale_max ! Maximum allowable value for Lscale [m]
173 :
174 : real( kind = core_rknd ), dimension(ngrdcol), intent(in) :: &
175 : mu ! mu Fractional extrainment rate per unit altitude [1/m]
176 :
177 : real( kind = core_rknd ), intent(in) :: &
178 : lmin ! CLUBB tunable parameter lmin
179 :
180 : integer, intent(in) :: &
181 : saturation_formula ! Integer that stores the saturation formula to be used
182 :
183 : logical, intent(in) :: &
184 : l_implemented ! Flag for CLUBB being implemented in a larger model
185 :
186 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
187 : Lscale, & ! Mixing length [m]
188 : Lscale_up, & ! Mixing length up [m]
189 : Lscale_down ! Mixing length down [m]
190 :
191 : type (stats_metadata_type), intent(in) :: &
192 : stats_metadata
193 :
194 : !--------------------------------- Local Variables ---------------------------------
195 :
196 : integer :: i, j, k, start_index
197 :
198 : real( kind = core_rknd ) :: tke, CAPE_incr
199 :
200 : real( kind = core_rknd ) :: dCAPE_dz_j, dCAPE_dz_j_minus_1, dCAPE_dz_j_plus_1
201 :
202 : ! Temporary 2D arrays to store calculations to speed runtime
203 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
204 17870112 : exp_mu_dzm, &
205 17870112 : invrs_dzm_on_mu, &
206 17870112 : grav_on_thvm, &
207 17870112 : Lv_coef, &
208 17870112 : entrain_coef, &
209 17870112 : thl_par_j_precalc, &
210 17870112 : rt_par_j_precalc, &
211 17870112 : tl_par_1, &
212 17870112 : rt_par_1, &
213 17870112 : rsatl_par_1, &
214 17870112 : thl_par_1, &
215 17870112 : dCAPE_dz_1, &
216 17870112 : s_par_1, &
217 17870112 : rc_par_1, &
218 17870112 : CAPE_incr_1, &
219 17870112 : thv_par_1
220 :
221 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
222 17870112 : tke_i
223 :
224 : ! Minimum value for Lscale that will taper off with height
225 : real( kind = core_rknd ) :: lminh
226 :
227 : ! Parcel quantities at grid level j
228 : real( kind = core_rknd ) :: thl_par_j, rt_par_j, rc_par_j, thv_par_j
229 :
230 : ! Used in latent heating calculation
231 : real( kind = core_rknd ) :: tl_par_j, rsatl_par_j, s_par_j
232 :
233 : ! Variables to make L nonlocal
234 : real( kind = core_rknd ) :: Lscale_up_max_alt, Lscale_down_min_alt
235 :
236 : ! Variables used to precalculate values
237 : real( kind = core_rknd ) :: &
238 : Lv2_coef, &
239 : tl_par_j_sqd, &
240 : invrs_dCAPE_diff, &
241 : invrs_Lscale_sfclyr_depth
242 :
243 : ! --------------------------------- Begin Code ---------------------------------
244 :
245 : !$acc enter data create( exp_mu_dzm, invrs_dzm_on_mu, grav_on_thvm, Lv_coef, &
246 : !$acc entrain_coef, thl_par_j_precalc, rt_par_j_precalc, &
247 : !$acc tl_par_1, rt_par_1, rsatl_par_1, thl_par_1, dCAPE_dz_1, &
248 : !$acc s_par_1, rc_par_1, CAPE_incr_1, thv_par_1, tke_i )
249 :
250 : !$acc parallel loop gang vector default(present)
251 149194656 : do i = 1, ngrdcol
252 149194656 : if ( abs(mu(i)) < eps ) then
253 0 : err_code = clubb_fatal_error
254 0 : print *, "mu = ", mu(i)
255 : end if
256 : end do
257 : !$acc end parallel loop
258 :
259 8935056 : if ( err_code == clubb_fatal_error ) then
260 0 : write(fstderr,*) "Entrainment rate mu cannot be 0"
261 0 : error stop "Fatal error in subroutine compute_mixing_length"
262 : end if
263 :
264 : ! Calculate initial turbulent kinetic energy for each grid level
265 8935056 : tke_i = zm2zt( nz, ngrdcol, gr, em )
266 :
267 : ! Initialize arrays and precalculate values for computational efficiency
268 : !$acc parallel loop gang vector collapse(2) default(present)
269 149194656 : do i = 1, ngrdcol
270 12071260656 : do k = 1, nz
271 :
272 : ! Initialize up and down arrays
273 11922066000 : Lscale_up(i,k) = zlmin
274 11922066000 : Lscale_down(i,k) = zlmin
275 :
276 : ! Precalculate values to avoid unnecessary calculations later
277 11922066000 : exp_mu_dzm(i,k) = exp( -mu(i) * gr%dzm(i,k) )
278 11922066000 : invrs_dzm_on_mu(i,k) = ( gr%invrs_dzm(i,k) ) / mu(i)
279 11922066000 : grav_on_thvm(i,k) = grav / thvm(i,k)
280 11922066000 : Lv_coef(i,k) = Lv / ( exner(i,k) * cp ) - ep2 * thv_ds(i,k)
281 12062325600 : entrain_coef(i,k) = ( one - exp_mu_dzm(i,k) ) * invrs_dzm_on_mu(i,k)
282 :
283 : end do
284 : end do
285 : !$acc end parallel loop
286 :
287 : !$acc parallel loop gang vector default(present)
288 149194656 : do i = 1, ngrdcol
289 :
290 : ! Avoid uninitialized memory (these values are not used in Lscale)
291 140259600 : Lscale_up(i,1) = zero
292 149194656 : Lscale_down(i,1) = zero
293 : end do
294 : !$acc end parallel loop
295 :
296 : ! Precalculations of single values to avoid unnecessary calculations later
297 : Lv2_coef = ep * Lv**2 / ( Rd * cp )
298 : invrs_Lscale_sfclyr_depth = one / Lscale_sfclyr_depth
299 :
300 :
301 : ! ---------------- Upwards Length Scale Calculation ----------------
302 :
303 : ! Precalculate values for upward Lscale, these are useful only if a parcel can rise
304 : ! more than one level. They are used in the equations that calculate thl and rt
305 : ! recursively for a parcel as it ascends
306 :
307 : !$acc parallel loop gang vector collapse(2) default(present)
308 149194656 : do i = 1, ngrdcol
309 11790741456 : do j = 2, nz-1
310 :
311 34924640400 : thl_par_j_precalc(i,j) = thlm(i,j) - thlm(i,j-1) * exp_mu_dzm(i,j-1) &
312 34924640400 : - ( thlm(i,j) - thlm(i,j-1) ) * entrain_coef(i,j-1)
313 :
314 : rt_par_j_precalc(i,j) = rtm(i,j) - rtm(i,j-1) * exp_mu_dzm(i,j-1) &
315 11781806400 : - ( rtm(i,j) - rtm(i,j-1) ) * entrain_coef(i,j-1)
316 : end do
317 : end do
318 : !$acc end parallel loop
319 :
320 : ! Calculate the initial change in TKE for each level. This is done for computational
321 : ! efficiency, it helps because there will be at least one calculation for each grid level,
322 : ! meaning the first one can be done for every grid level and therefore the calculations can
323 : ! be vectorized, clubb:ticket:834. After the initial calculation however, it is uncertain
324 : ! how many more iterations should be done for each individual grid level, and calculating
325 : ! one change in TKE for each level until all are exhausted will result in many unnessary
326 : ! and expensive calculations.
327 :
328 : ! Calculate initial thl, tl, and rt for parcels at each grid level
329 : !$acc parallel loop gang vector collapse(2) default(present)
330 149194656 : do i = 1, ngrdcol
331 11790741456 : do j = 3, nz
332 :
333 11641546800 : thl_par_1(i,j) = thlm(i,j) - ( thlm(i,j) - thlm(i,j-1) ) * entrain_coef(i,j-1)
334 :
335 11641546800 : tl_par_1(i,j) = thl_par_1(i,j) * exner(i,j)
336 :
337 11781806400 : rt_par_1(i,j) = rtm(i,j) - ( rtm(i,j) - rtm(i,j-1) ) * entrain_coef(i,j-1)
338 :
339 : end do
340 : end do
341 : !$acc end parallel loop
342 :
343 :
344 : ! Caclculate initial rsatl for parcels at each grid level
345 :
346 : ! The entire pressure and temperature arrays are passed as
347 : ! argument and the sub-arrays are choosen using
348 : ! start_index. This workaround is used to solve
349 : ! subarray issues with OpenACC.
350 : ! rsatl_par_1(i,3:) = sat_mixrat_liq_acc( nz-2, ngrdcol, p_in_Pa(i,3:), tl_par_1(i,3:) )
351 : ! since subarray 3:, the start_index is 3 and it is an optional argument
352 8935056 : start_index = 3
353 8935056 : rsatl_par_1 = sat_mixrat_liq( nz, ngrdcol, p_in_Pa, tl_par_1, saturation_formula, start_index )
354 :
355 : ! Calculate initial dCAPE_dz and CAPE_incr for parcels at each grid level
356 : !$acc parallel loop gang vector default(present)
357 149194656 : do i = 1, ngrdcol
358 11781806400 : do j = 3, nz
359 :
360 11641546800 : tl_par_j_sqd = tl_par_1(i,j)**2
361 :
362 : ! s from Lewellen and Yoh 1993 (LY) eqn. 1
363 : ! s = ( rt - rsatl ) / ( 1 + beta * rsatl )
364 : ! and SD's beta (eqn. 8),
365 : ! beta = ep * ( Lv / ( Rd * tl ) ) * ( Lv / ( cp * tl ) )
366 : !
367 : ! Simplified by multiplying top and bottom by tl^2 to avoid a divide and precalculating
368 : ! ep * Lv**2 / ( Rd * cp )
369 : s_par_1(i,j) = ( rt_par_1(i,j) - rsatl_par_1(i,j) ) * tl_par_j_sqd &
370 11641546800 : / ( tl_par_j_sqd + Lv2_coef * rsatl_par_1(i,j) )
371 :
372 11641546800 : rc_par_1(i,j) = max( s_par_1(i,j), zero_threshold )
373 :
374 : ! theta_v of entraining parcel at grid level j
375 11641546800 : thv_par_1(i,j) = thl_par_1(i,j) + ep1 * thv_ds(i,j) * rt_par_1(i,j) + Lv_coef(i,j) * rc_par_1(i,j)
376 :
377 :
378 : ! dCAPE/dz = g * ( thv_par - thvm ) / thvm.
379 11641546800 : dCAPE_dz_1(i,j) = grav_on_thvm(i,j) * ( thv_par_1(i,j) - thvm(i,j) )
380 :
381 : ! CAPE_incr = INT(z_0:z_1) g * ( thv_par - thvm ) / thvm dz
382 : ! Trapezoidal estimate between grid levels, dCAPE at z_0 = 0 for this initial calculation
383 11781806400 : CAPE_incr_1(i,j) = one_half * dCAPE_dz_1(i,j) * gr%dzm(i,j-1)
384 :
385 : end do
386 :
387 :
388 : ! Calculate Lscale_up for each grid level. If the TKE from a parcel has not been completely
389 : ! exhausted by the initial change then continue the exhaustion calculations here for a single
390 : ! grid level at a time until the TKE is exhausted.
391 :
392 140259600 : Lscale_up_max_alt = zero ! Set initial max value for Lscale_up to 0
393 11650481856 : do k = 2, nz-2
394 :
395 : ! If the initial turbulent kinetic energy (tke) has not been exhausted for this grid level
396 11501287200 : if ( tke_i(i,k) + CAPE_incr_1(i,k+1) > zero ) then
397 :
398 : ! Calculate new TKE for parcel
399 915834623 : tke = tke_i(i,k) + CAPE_incr_1(i,k+1)
400 :
401 : ! Set j to 2 levels above current Lscale_up level, this is because we've already
402 : ! determined that the parcel can rise at least 1 full level
403 915834623 : j = k + 2
404 :
405 : ! Set initial thl, rt, and dCAPE_dz to the values found by the intial calculations
406 915834623 : thl_par_j = thl_par_1(i,k+1)
407 915834623 : rt_par_j = rt_par_1(i,k+1)
408 915834623 : dCAPE_dz_j_minus_1 = dCAPE_dz_1(i,k+1)
409 :
410 :
411 : ! Continue change in TKE calculations until it is exhausted or the max grid
412 : ! level has been reached. j is the next grid level above the level that can
413 : ! be reached for a parcel starting at level k. If TKE is exhausted in this loop
414 : ! that means the parcel starting at k cannot reach level j, but has reached j-1
415 4657532814 : do while ( j < nz )
416 :
417 : ! thl, rt of parcel are conserved except for entrainment
418 : !
419 : ! The values of thl_env and rt_env are treated as changing linearly for a parcel
420 : ! of air ascending from level j-1 to level j
421 :
422 : ! theta_l of the parcel starting at grid level k, and currenly
423 : ! at grid level j
424 : !
425 : ! d(thl_par)/dz = - mu * ( thl_par - thl_env )
426 4657532814 : thl_par_j = thl_par_j_precalc(i,j) + thl_par_j * exp_mu_dzm(i,j-1)
427 :
428 :
429 : ! r_t of the parcel starting at grid level k, and currenly
430 : ! at grid level j
431 : !
432 : ! d(rt_par)/dz = - mu * ( rt_par - rt_env )
433 4657532814 : rt_par_j = rt_par_j_precalc(i,j) + rt_par_j * exp_mu_dzm(i,j-1)
434 :
435 :
436 : ! Include effects of latent heating on Lscale_up 6/12/00
437 : ! Use thermodynamic formula of Bougeault 1981 JAS Vol. 38, 2416
438 : ! Probably should use properties of bump 1 in Gaussian, not mean!!!
439 :
440 4657532814 : tl_par_j = thl_par_j*exner(i,j)
441 :
442 4657532814 : rsatl_par_j = sat_mixrat_liq( p_in_Pa(i,j), tl_par_j, saturation_formula )
443 :
444 4657532814 : tl_par_j_sqd = tl_par_j**2
445 :
446 : ! s from Lewellen and Yoh 1993 (LY) eqn. 1
447 : ! s = ( rt - rsatl ) / ( 1 + beta * rsatl )
448 : ! and SD's beta (eqn. 8),
449 : ! beta = ep * ( Lv / ( Rd * tl ) ) * ( Lv / ( cp * tl ) )
450 : !
451 : ! Simplified by multiplying top and bottom by tl^2 to avoid a
452 : ! divide and precalculating ep * Lv**2 / ( Rd * cp )
453 : s_par_j = ( rt_par_j - rsatl_par_j ) * tl_par_j_sqd &
454 4657532814 : / ( tl_par_j_sqd + Lv2_coef * rsatl_par_j )
455 :
456 4657532814 : rc_par_j = max( s_par_j, zero_threshold )
457 :
458 : ! theta_v of entraining parcel at grid level j
459 : thv_par_j = thl_par_j + ep1 * thv_ds(i,j) * rt_par_j &
460 4657532814 : + Lv_coef(i,j) * rc_par_j
461 :
462 : ! dCAPE/dz = g * ( thv_par - thvm ) / thvm.
463 4657532814 : dCAPE_dz_j = grav_on_thvm(i,j) * ( thv_par_j - thvm(i,j) )
464 :
465 : ! CAPE_incr = INT(z_0:z_1) g * ( thv_par - thvm ) / thvm dz
466 : ! Trapezoidal estimate between grid levels j and j-1
467 4657532814 : CAPE_incr = one_half * ( dCAPE_dz_j + dCAPE_dz_j_minus_1 ) * gr%dzm(i,j-1)
468 :
469 : ! Exit loop early if tke has been exhaused between level j and j+1
470 4657532814 : if ( tke + CAPE_incr <= zero ) then
471 : exit
472 : end if
473 :
474 : ! Save previous dCAPE value for next cycle
475 3741698191 : dCAPE_dz_j_minus_1 = dCAPE_dz_j
476 :
477 : ! Caclulate new TKE and increment j
478 3741698191 : tke = tke + CAPE_incr
479 4657532814 : j = j + 1
480 :
481 : enddo
482 :
483 :
484 : ! Add full grid level thickness for each grid level that was passed without the TKE
485 : ! being exhausted, difference between starting level (k) and last level passed (j-1)
486 915834623 : Lscale_up(i,k) = Lscale_up(i,k) + gr%zt(i,j-1) - gr%zt(i,k)
487 :
488 :
489 915834623 : if ( j < nz ) then
490 :
491 : ! Loop terminated early, meaning TKE was completely exhaused at grid level j.
492 : ! Add the thickness z - z_0 (where z_0 < z <= z_1) to Lscale_up.
493 :
494 915834623 : if ( abs( dCAPE_dz_j - dCAPE_dz_j_minus_1 ) * 2 <= &
495 : abs( dCAPE_dz_j + dCAPE_dz_j_minus_1 ) * eps ) then
496 :
497 : ! Special case where dCAPE/dz|_(z_1) - dCAPE/dz|_(z_0) = 0
498 : ! Find the remaining distance z - z_0 that it takes to
499 : ! exhaust the remaining TKE
500 :
501 0 : Lscale_up(i,k) = Lscale_up(i,k) + ( - tke / dCAPE_dz_j )
502 :
503 : else
504 :
505 : ! Case used for most scenarios where dCAPE/dz|_(z_1) /= dCAPE/dz|_(z_0)
506 : ! Find the remaining distance z - z_0 that it takes to exhaust the
507 : ! remaining TKE (tke_i), using the quadratic formula (only the
508 : ! negative (-) root works in this scenario).
509 915834623 : invrs_dCAPE_diff = one / ( dCAPE_dz_j - dCAPE_dz_j_minus_1 )
510 :
511 : Lscale_up(i,k) = Lscale_up(i,k) &
512 0 : - dCAPE_dz_j_minus_1 * invrs_dCAPE_diff * gr%dzm(i,j-1) &
513 : - sqrt( dCAPE_dz_j_minus_1**2 &
514 0 : - two * tke * gr%invrs_dzm(i,j-1) &
515 : * ( dCAPE_dz_j - dCAPE_dz_j_minus_1 ) ) &
516 915834623 : * invrs_dCAPE_diff * gr%dzm(i,j-1)
517 : endif
518 :
519 : end if
520 :
521 : else ! TKE for parcel at level (k) was exhaused before one full grid level
522 :
523 : ! Find the remaining distance z - z_0 that it takes to exhaust the
524 : ! remaining TKE (tke_i), using the quadratic formula. Simplified
525 : ! since dCAPE_dz_j_minus_1 = 0.0
526 : Lscale_up(i,k) = Lscale_up(i,k) - sqrt( - two * tke_i(i,k) &
527 0 : * gr%dzm(i,k) * dCAPE_dz_1(i,k+1) ) &
528 10585452577 : / dCAPE_dz_1(i,k+1)
529 : endif
530 :
531 :
532 : ! If a parcel at a previous grid level can rise past the parcel at this grid level
533 : ! then this one should also be able to rise up to that height. This feature insures
534 : ! that the profile of Lscale_up will be smooth, thus reducing numerical instability.
535 11641546800 : if ( gr%zt(i,k) + Lscale_up(i,k) < Lscale_up_max_alt ) then
536 :
537 : ! A lower starting parcel can ascend higher than this one, set height to the max
538 : ! that any lower starting parcel can ascend to
539 1191995387 : Lscale_up(i,k) = Lscale_up_max_alt - gr%zt(i,k)
540 : else
541 :
542 : ! This parcel can ascend higher than any below it, save final height
543 : Lscale_up_max_alt = Lscale_up(i,k) + gr%zt(i,k)
544 : end if
545 :
546 :
547 : end do
548 : end do
549 : !$acc end parallel loop
550 :
551 : ! ---------------- Downwards Length Scale Calculation ----------------
552 :
553 : ! Precalculate values for downward Lscale, these are useful only if a parcel can descend
554 : ! more than one level. They are used in the equations that calculate thl and rt
555 : ! recursively for a parcel as it descends
556 : !$acc parallel loop gang vector collapse(2) default(present)
557 149194656 : do i = 1, ngrdcol
558 11790741456 : do j = 2, nz-1
559 :
560 34924640400 : thl_par_j_precalc(i,j) = thlm(i,j) - thlm(i,j+1) * exp_mu_dzm(i,j) &
561 34924640400 : - ( thlm(i,j) - thlm(i,j+1) ) * entrain_coef(i,j)
562 :
563 : rt_par_j_precalc(i,j) = rtm(i,j) - rtm(i,j+1) * exp_mu_dzm(i,j) &
564 11781806400 : - ( rtm(i,j) - rtm(i,j+1) ) * entrain_coef(i,j)
565 : end do
566 : end do
567 : !$acc end parallel loop
568 :
569 : ! Calculate the initial change in TKE for each level. This is done for computational
570 : ! efficiency, it helps because there will be at least one calculation for each grid level,
571 : ! meaning the first one can be done for every grid level and therefore the calculations can
572 : ! be vectorized, clubb:ticket:834. After the initial calculation however, it is uncertain
573 : ! how many more iterations should be done for each individual grid level, and calculating
574 : ! one change in TKE for each level until all are exhausted will result in many unnessary
575 : ! and expensive calculations.
576 :
577 : ! Calculate initial thl, tl, and rt for parcels at each grid level
578 : !$acc parallel loop gang vector collapse(2) default(present)
579 149194656 : do i = 1, ngrdcol
580 11790741456 : do j = 2, nz-1
581 :
582 11641546800 : thl_par_1(i,j) = thlm(i,j) - ( thlm(i,j) - thlm(i,j+1) ) * entrain_coef(i,j)
583 :
584 11641546800 : tl_par_1(i,j) = thl_par_1(i,j) * exner(i,j)
585 :
586 11781806400 : rt_par_1(i,j) = rtm(i,j) - ( rtm(i,j) - rtm(i,j+1) ) * entrain_coef(i,j)
587 :
588 : end do
589 : end do
590 : !$acc end parallel loop
591 :
592 : ! Caclculate initial rsatl for parcels at each grid level, this function is elemental
593 :
594 : ! The entire pressure and temperature arrays are passed as
595 : ! argument and the sub-arrays are choosen using
596 : ! start_index. This workaround is used to solve
597 : ! subarray issues with OpenACC.
598 : ! rsatl_par_1(i,2:) = sat_mixrat_liq_acc( nz-1, p_in_Pa(i,2:), tl_par_1(i,2:) )
599 : ! since subarray 2:, the start_index is 2 and it is an optional argument
600 8935056 : start_index = 2
601 8935056 : rsatl_par_1 = sat_mixrat_liq( nz, ngrdcol, p_in_Pa, tl_par_1, saturation_formula, start_index )
602 :
603 : ! Calculate initial dCAPE_dz and CAPE_incr for parcels at each grid level
604 : !$acc parallel loop gang vector default(present)
605 149194656 : do i = 1, ngrdcol
606 11781806400 : do j = 2, nz-1
607 :
608 11641546800 : tl_par_j_sqd = tl_par_1(i,j)**2
609 :
610 : ! s from Lewellen and Yoh 1993 (LY) eqn. 1
611 : ! s = ( rt - rsatl ) / ( 1 + beta * rsatl )
612 : ! and SD's beta (eqn. 8),
613 : ! beta = ep * ( Lv / ( Rd * tl ) ) * ( Lv / ( cp * tl ) )
614 : !
615 : ! Simplified by multiplying top and bottom by tl^2 to avoid a divide and precalculating
616 : ! ep * Lv**2 / ( Rd * cp )
617 : s_par_1(i,j) = ( rt_par_1(i,j) - rsatl_par_1(i,j) ) * tl_par_j_sqd &
618 11641546800 : / ( tl_par_j_sqd + Lv2_coef * rsatl_par_1(i,j) )
619 :
620 11641546800 : rc_par_1(i,j) = max( s_par_1(i,j), zero_threshold )
621 :
622 : ! theta_v of entraining parcel at grid level j
623 11641546800 : thv_par_1(i,j) = thl_par_1(i,j) + ep1 * thv_ds(i,j) * rt_par_1(i,j) + Lv_coef(i,j) * rc_par_1(i,j)
624 :
625 : ! dCAPE/dz = g * ( thv_par - thvm ) / thvm.
626 11641546800 : dCAPE_dz_1(i,j) = grav_on_thvm(i,j) * ( thv_par_1(i,j) - thvm(i,j) )
627 :
628 : ! CAPE_incr = INT(z_0:z_1) g * ( thv_par - thvm ) / thvm dz
629 : ! Trapezoidal estimate between grid levels, dCAPE at z_0 = 0 for this initial calculation
630 11781806400 : CAPE_incr_1(i,j) = one_half * dCAPE_dz_1(i,j) * gr%dzm(i,j)
631 :
632 : end do
633 :
634 :
635 : ! Calculate Lscale_down for each grid level. If the TKE from a parcel has not been completely
636 : ! exhausted by the initial change then continue the exhaustion calculations here for a single
637 : ! grid level at a time until the TKE is exhausted.
638 :
639 140259600 : Lscale_down_min_alt = gr%zt(i,nz) ! Set initial min value for Lscale_down to max zt
640 11790741456 : do k = nz, 3, -1
641 :
642 : ! If the initial turbulent kinetic energy (tke) has not been exhausted for this grid level
643 11641546800 : if ( tke_i(i,k) - CAPE_incr_1(i,k-1) > zero ) then
644 :
645 : ! Calculate new TKE for parcel
646 741946634 : tke = tke_i(i,k) - CAPE_incr_1(i,k-1)
647 :
648 : ! Set j to 2 levels below current Lscale_down level, this is because we've already
649 : ! determined that the parcel can descend at least 1 full level
650 741946634 : j = k - 2
651 :
652 : ! Set initial thl, rt, and dCAPE_dz to the values found by the intial calculations
653 741946634 : thl_par_j = thl_par_1(i,k-1)
654 741946634 : rt_par_j = rt_par_1(i,k-1)
655 741946634 : dCAPE_dz_j_plus_1 = dCAPE_dz_1(i,k-1)
656 :
657 :
658 : ! Continue change in TKE calculations until it is exhausted or the min grid
659 : ! level has been reached. j is the next grid level below the level that can
660 : ! be reached for a parcel starting at level k. If TKE is exhausted in this loop
661 : ! that means the parcel starting at k cannot sink to level j, but can sink to j+1
662 2215833479 : do while ( j >= 2 )
663 :
664 : ! thl, rt of parcel are conserved except for entrainment
665 : !
666 : ! The values of thl_env and rt_env are treated as changing linearly for a parcel
667 : ! of air descending from level j to level j-1
668 :
669 : ! theta_l of the parcel starting at grid level k, and currenly
670 : ! at grid level j
671 : !
672 : ! d(thl_par)/dz = - mu * ( thl_par - thl_env )
673 1711112817 : thl_par_j = thl_par_j_precalc(i,j) + thl_par_j * exp_mu_dzm(i,j)
674 :
675 :
676 : ! r_t of the parcel starting at grid level k, and currenly
677 : ! at grid level j
678 : !
679 : ! d(rt_par)/dz = - mu * ( rt_par - rt_env )
680 1711112817 : rt_par_j = rt_par_j_precalc(i,j) + rt_par_j * exp_mu_dzm(i,j)
681 :
682 :
683 : ! Include effects of latent heating on Lscale_up 6/12/00
684 : ! Use thermodynamic formula of Bougeault 1981 JAS Vol. 38, 2416
685 : ! Probably should use properties of bump 1 in Gaussian, not mean!!!
686 :
687 1711112817 : tl_par_j = thl_par_j*exner(i,j)
688 :
689 1711112817 : rsatl_par_j = sat_mixrat_liq( p_in_Pa(i,j), tl_par_j, saturation_formula )
690 :
691 1711112817 : tl_par_j_sqd = tl_par_j**2
692 :
693 : ! s from Lewellen and Yoh 1993 (LY) eqn. 1
694 : ! s = ( rt - rsatl ) / ( 1 + beta * rsatl )
695 : ! and SD's beta (eqn. 8),
696 : ! beta = ep * ( Lv / ( Rd * tl ) ) * ( Lv / ( cp * tl ) )
697 : !
698 : ! Simplified by multiplying top and bottom by tl^2 to avoid a
699 : ! divide and precalculating ep * Lv**2 / ( Rd * cp )
700 : s_par_j = (rt_par_j - rsatl_par_j) * tl_par_j_sqd &
701 1711112817 : / ( tl_par_j_sqd + Lv2_coef * rsatl_par_j )
702 :
703 1711112817 : rc_par_j = max( s_par_j, zero_threshold )
704 :
705 : ! theta_v of entraining parcel at grid level j
706 1711112817 : thv_par_j = thl_par_j + ep1 * thv_ds(i,j) * rt_par_j + Lv_coef(i,j) * rc_par_j
707 :
708 : ! dCAPE/dz = g * ( thv_par - thvm ) / thvm.
709 1711112817 : dCAPE_dz_j = grav_on_thvm(i,j) * ( thv_par_j - thvm(i,j) )
710 :
711 : ! CAPE_incr = INT(z_0:z_1) g * ( thv_par - thvm ) / thvm dz
712 : ! Trapezoidal estimate between grid levels j+1 and j
713 1711112817 : CAPE_incr = one_half * ( dCAPE_dz_j + dCAPE_dz_j_plus_1 ) * gr%dzm(i,j)
714 :
715 : ! Exit loop early if tke has been exhaused between level j+1 and j
716 1711112817 : if ( tke - CAPE_incr <= zero ) then
717 : exit
718 : endif
719 :
720 : ! Save previous dCAPE value for next cycle
721 1473886845 : dCAPE_dz_j_plus_1 = dCAPE_dz_j
722 :
723 : ! Caclulate new TKE and increment j
724 1473886845 : tke = tke - CAPE_incr
725 1711112817 : j = j - 1
726 :
727 : enddo
728 :
729 : ! Add full grid level thickness for each grid level that was passed without the TKE
730 : ! being exhausted, difference between starting level (k) and last level passed (j+1)
731 741946634 : Lscale_down(i,k) = Lscale_down(i,k) + gr%zt(i,k) - gr%zt(i,j+1)
732 :
733 :
734 741946634 : if ( j >= 2 ) then
735 :
736 : ! Loop terminated early, meaning TKE was completely exhaused at grid level j.
737 : ! Add the thickness z - z_0 (where z_0 < z <= z_1) to Lscale_up.
738 :
739 237225972 : if ( abs( dCAPE_dz_j - dCAPE_dz_j_plus_1 ) * 2 <= &
740 : abs( dCAPE_dz_j + dCAPE_dz_j_plus_1 ) * eps ) then
741 :
742 : ! Special case where dCAPE/dz|_(z_(-1)) - dCAPE/dz|_(z_0) = 0
743 : ! Find the remaining distance z_0 - z that it takes to
744 : ! exhaust the remaining TKE
745 :
746 0 : Lscale_down(i,k) = Lscale_down(i,k) + ( tke / dCAPE_dz_j )
747 :
748 : else
749 :
750 : ! Case used for most scenarios where dCAPE/dz|_(z_(-1)) /= dCAPE/dz|_(z_0)
751 : ! Find the remaining distance z_0 - z that it takes to exhaust the
752 : ! remaining TKE (tke_i), using the quadratic formula (only the
753 : ! negative (-) root works in this scenario) -- however, the
754 : ! negative (-) root is divided by another negative (-) factor,
755 : ! which results in an overall plus (+) sign in front of the
756 : ! square root term in the equation below).
757 237225972 : invrs_dCAPE_diff = one / ( dCAPE_dz_j - dCAPE_dz_j_plus_1 )
758 :
759 : Lscale_down(i,k) = Lscale_down(i,k) &
760 0 : - dCAPE_dz_j_plus_1 * invrs_dCAPE_diff * gr%dzm(i,j) &
761 : + sqrt( dCAPE_dz_j_plus_1**2 &
762 0 : + two * tke * gr%invrs_dzm(i,j) &
763 : * ( dCAPE_dz_j - dCAPE_dz_j_plus_1 ) ) &
764 237225972 : * invrs_dCAPE_diff * gr%dzm(i,j)
765 : endif
766 :
767 : end if
768 :
769 : else ! TKE for parcel at level (k) was exhaused before one full grid level
770 :
771 : ! Find the remaining distance z_0 - z that it takes to exhaust the
772 : ! remaining TKE (tke_i), using the quadratic formula. Simplified
773 : ! since dCAPE_dz_j_plus_1 = 0.0
774 10899600166 : Lscale_down(i,k) = Lscale_down(i,k) + sqrt( two * tke_i(i,k) &
775 0 : * gr%dzm(i,k-1) * dCAPE_dz_1(i,k-1) ) &
776 10899600166 : / dCAPE_dz_1(i,k-1)
777 : endif
778 :
779 : ! If a parcel at a previous grid level can descend past the parcel at this grid level
780 : ! then this one should also be able to descend down to that height. This feature insures
781 : ! that the profile of Lscale_down will be smooth, thus reducing numerical instability.
782 11781806400 : if ( gr%zt(i,k) - Lscale_down(i,k) > Lscale_down_min_alt ) then
783 235041601 : Lscale_down(i,k) = gr%zt(i,k) - Lscale_down_min_alt
784 : else
785 : Lscale_down_min_alt = gr%zt(i,k) - Lscale_down(i,k)
786 : end if
787 :
788 : end do
789 : end do
790 : !$acc end parallel loop
791 :
792 : ! ---------------- Final Lscale Calculation ----------------
793 :
794 : !$acc parallel loop gang vector default(present)
795 149194656 : do i = 1, ngrdcol
796 11922066000 : do k = 2, nz, 1
797 :
798 : ! Make lminh a linear function starting at value lmin at the bottom
799 : ! and going to zero at 500 meters in altitude.
800 11781806400 : if( l_implemented ) then
801 :
802 : ! Within a host model, increase mixing length in 500 m layer above *ground*
803 0 : lminh = max( zero_threshold, Lscale_sfclyr_depth - ( gr%zt(i,k) - gr%zm(i,1) ) ) &
804 11781806400 : * lmin * invrs_Lscale_sfclyr_depth
805 : else
806 :
807 : ! In standalone mode, increase mixing length in 500 m layer above *mean sea level*
808 0 : lminh = max( zero_threshold, Lscale_sfclyr_depth - gr%zt(i,k) ) &
809 0 : * lmin * invrs_Lscale_sfclyr_depth
810 : end if
811 :
812 11781806400 : Lscale_up(i,k) = max( lminh, Lscale_up(i,k) )
813 11781806400 : Lscale_down(i,k) = max( lminh, Lscale_down(i,k) )
814 :
815 : ! When L is large, turbulence is weakly damped
816 : ! When L is small, turbulence is strongly damped
817 : ! Use a geometric mean to determine final Lscale so that L tends to become small
818 : ! if either Lscale_up or Lscale_down becomes small.
819 11922066000 : Lscale(i,k) = sqrt( Lscale_up(i,k) * Lscale_down(i,k) )
820 :
821 : enddo
822 :
823 : ! Set the value of Lscale at the upper and lower boundaries.
824 140259600 : Lscale(i,1) = Lscale(i,2)
825 140259600 : Lscale(i,nz) = Lscale(i,nz-1)
826 :
827 : ! Vince Larson limited Lscale to allow host
828 : ! model to take over deep convection. 13 Feb 2008.
829 12071260656 : Lscale(i,:) = min( Lscale(i,:), Lscale_max(i) )
830 :
831 : end do
832 : !$acc end parallel loop
833 :
834 : ! Ensure that no Lscale values are NaN
835 8935056 : if ( clubb_at_least_debug_level( 1 ) ) then
836 :
837 : !$acc update host( Lscale, Lscale_up, Lscale_down, &
838 : !$acc thvm, thlm, rtm, em, exner, p_in_Pa, thv_ds )
839 :
840 0 : do i = 1, ngrdcol
841 0 : call length_check( nz, Lscale(i,:), Lscale_up(i,:), Lscale_down(i,:) ) ! intent(in)
842 : end do
843 :
844 0 : if ( err_code == clubb_fatal_error ) then
845 :
846 0 : write(fstderr,*) "Errors in compute_mixing_length subroutine"
847 :
848 0 : write(fstderr,*) "Intent(in)"
849 :
850 0 : write(fstderr,*) "thvm = ", thvm
851 0 : write(fstderr,*) "thlm = ", thlm
852 0 : write(fstderr,*) "rtm = ", rtm
853 0 : write(fstderr,*) "em = ", em
854 0 : write(fstderr,*) "exner = ", exner
855 0 : write(fstderr,*) "p_in_Pa = ", p_in_Pa
856 0 : write(fstderr,*) "thv_ds = ", thv_ds
857 :
858 0 : write(fstderr,*) "Intent(out)"
859 :
860 0 : write(fstderr,*) "Lscale = ", Lscale
861 0 : write(fstderr,*) "Lscale_up = ", Lscale_up
862 0 : write(fstderr,*) "Lscale_down = ", Lscale_down
863 :
864 : endif ! Fatal error
865 :
866 : end if
867 :
868 : !$acc exit data delete( exp_mu_dzm, invrs_dzm_on_mu, grav_on_thvm, Lv_coef, &
869 : !$acc entrain_coef, thl_par_j_precalc, rt_par_j_precalc, &
870 : !$acc tl_par_1, rt_par_1, rsatl_par_1, thl_par_1, dCAPE_dz_1, &
871 : !$acc s_par_1, rc_par_1, CAPE_incr_1, thv_par_1, tke_i )
872 :
873 8935056 : return
874 :
875 : end subroutine compute_mixing_length
876 :
877 : !===============================================================================
878 8935056 : subroutine calc_Lscale_directly ( ngrdcol, nz, gr, &
879 8935056 : l_implemented, p_in_Pa, exner, rtm, &
880 8935056 : thlm, thvm, newmu, rtp2_zt, thlp2_zt, rtpthlp_zt, &
881 8935056 : pdf_params, em, thv_ds_zt, Lscale_max, lmin, &
882 8935056 : clubb_params, &
883 : saturation_formula, &
884 : l_Lscale_plume_centered, &
885 : stats_metadata, &
886 8935056 : stats_zt, &
887 8935056 : Lscale, Lscale_up, Lscale_down)
888 :
889 : use constants_clubb, only: &
890 : thl_tol, &
891 : rt_tol, &
892 : one_half, &
893 : one_third, &
894 : one, &
895 : three, &
896 : unused_var
897 :
898 : use parameter_indices, only: &
899 : nparams, &
900 : iLscale_mu_coef, &
901 : iLscale_pert_coef
902 :
903 : use grid_class, only: &
904 : grid ! Type
905 :
906 : use clubb_precision, only: &
907 : core_rknd
908 :
909 : use stats_variables, only: &
910 : stats_metadata_type
911 :
912 : use pdf_parameter_module, only: &
913 : pdf_parameter
914 :
915 : use stats_type_utilities, only: &
916 : stat_update_var
917 :
918 : use error_code, only: &
919 : clubb_at_least_debug_level, & ! Procedure
920 : err_code, & ! Error Indicator
921 : clubb_fatal_error ! Constant
922 :
923 : use constants_clubb, only: &
924 : fstderr ! Variable(s)
925 :
926 : use stats_type, only: &
927 : stats ! Type
928 :
929 : implicit none
930 :
931 : !--------------------------------- Input Variables ---------------------------------
932 : integer, intent(in) :: &
933 : nz, &
934 : ngrdcol
935 :
936 : type (grid), target, intent(in) :: &
937 : gr
938 :
939 : logical, intent(in) :: &
940 : l_implemented ! True if CLUBB is being run within a large-scale hostmodel,
941 : ! rather than a standalone single-column model.
942 :
943 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
944 : rtp2_zt, &
945 : thlp2_zt, &
946 : rtpthlp_zt, &
947 : thlm, &
948 : thvm, &
949 : rtm, &
950 : em, &
951 : p_in_Pa, & ! Air pressure (thermodynamic levels) [Pa]
952 : exner, &
953 : thv_ds_zt
954 :
955 : real( kind = core_rknd ), dimension(ngrdcol), intent(in) :: &
956 : newmu, &
957 : Lscale_max
958 :
959 : real( kind = core_rknd ), intent(in) :: &
960 : lmin
961 :
962 : type (pdf_parameter), intent(in) :: &
963 : pdf_params ! PDF Parameters [units vary]
964 :
965 : real( kind = core_rknd ), dimension(ngrdcol,nparams), intent(in) :: &
966 : clubb_params ! Array of CLUBB's tunable parameters [units vary]
967 :
968 : integer, intent(in) :: &
969 : saturation_formula ! Integer that stores the saturation formula to be used
970 :
971 : logical, intent(in) :: &
972 : l_Lscale_plume_centered ! Alternate that uses the PDF to compute the perturbed values
973 :
974 : type (stats_metadata_type), intent(in) :: &
975 : stats_metadata
976 :
977 : !--------------------------------- InOut Variables ---------------------------------
978 : type (stats), target, dimension(ngrdcol), intent(inout) :: &
979 : stats_zt
980 :
981 : !--------------------------------- Output Variables ---------------------------------
982 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
983 : Lscale, & ! Mixing length [m]
984 : Lscale_up, & ! Mixing length up [m]
985 : Lscale_down ! Mixing length down [m]
986 :
987 : !--------------------------------- Local Variables ---------------------------------
988 : integer :: k, i
989 :
990 : logical, parameter :: &
991 : l_avg_Lscale = .false. ! Lscale is calculated in subroutine compute_mixing_length
992 : ! if l_avg_Lscale is true, compute_mixing_length is called two additional
993 : ! times with
994 : ! perturbed values of rtm and thlm. An average value of Lscale
995 : ! from the three calls to compute_mixing_length is then calculated.
996 : ! This reduces temporal noise in RICO, BOMEX, LBA, and other cases.
997 :
998 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
999 8935056 : sign_rtpthlp_zt, & ! Sign of the covariance rtpthlp [-]
1000 17870112 : Lscale_pert_1, Lscale_pert_2, & ! For avg. calculation of Lscale [m]
1001 8935056 : thlm_pert_1, thlm_pert_2, & ! For avg. calculation of Lscale [K]
1002 8935056 : rtm_pert_1, rtm_pert_2, & ! For avg. calculation of Lscale [kg/kg]
1003 8935056 : thlm_pert_pos_rt, thlm_pert_neg_rt, & ! For avg. calculation of Lscale [K]
1004 8935056 : rtm_pert_pos_rt, rtm_pert_neg_rt ! For avg. calculation of Lscale [kg/kg]
1005 :
1006 : real( kind = core_rknd ), dimension(ngrdcol) :: &
1007 8935056 : mu_pert_1, mu_pert_2, &
1008 8935056 : mu_pert_pos_rt, mu_pert_neg_rt ! For l_Lscale_plume_centered
1009 :
1010 : real( kind = core_rknd ) :: &
1011 : Lscale_pert_coef
1012 :
1013 : !Lscale_weight Uncomment this if you need to use this vairable at some
1014 : !point.
1015 :
1016 : !--------------------------------- Begin Code ---------------------------------
1017 :
1018 : !$acc enter data create( sign_rtpthlp_zt, Lscale_pert_1, Lscale_pert_2, &
1019 : !$acc thlm_pert_1, thlm_pert_2, rtm_pert_1, rtm_pert_2, &
1020 : !$acc thlm_pert_pos_rt, thlm_pert_neg_rt, rtm_pert_pos_rt, &
1021 : !$acc rtm_pert_neg_rt, &
1022 : !$acc mu_pert_1, mu_pert_2, mu_pert_pos_rt, mu_pert_neg_rt )
1023 :
1024 8935056 : if ( clubb_at_least_debug_level( 0 ) ) then
1025 :
1026 8935056 : if ( l_Lscale_plume_centered .and. .not. l_avg_Lscale ) then
1027 0 : write(fstderr,*) "l_Lscale_plume_centered requires l_avg_Lscale"
1028 0 : write(fstderr,*) "Fatal error in advance_clubb_core"
1029 0 : err_code = clubb_fatal_error
1030 0 : return
1031 : end if
1032 :
1033 : end if
1034 :
1035 : if ( l_avg_Lscale .and. .not. l_Lscale_plume_centered ) then
1036 :
1037 : ! Call compute length two additional times with perturbed values
1038 : ! of rtm and thlm so that an average value of Lscale may be calculated.
1039 :
1040 : !$acc parallel loop gang vector collapse(2) default(present)
1041 : do k = 1, nz, 1
1042 : do i = 1, ngrdcol
1043 : sign_rtpthlp_zt(i,k) = sign( one, rtpthlp_zt(i,k) )
1044 : end do
1045 : end do
1046 : !$acc end parallel loop
1047 :
1048 : !$acc parallel loop gang vector collapse(2) default(present)
1049 : do k = 1, nz, 1
1050 : do i = 1, ngrdcol
1051 : rtm_pert_1(i,k) = rtm(i,k) + clubb_params(i,iLscale_pert_coef) &
1052 : * sqrt( max( rtp2_zt(i,k), rt_tol**2 ) )
1053 : end do
1054 : end do
1055 : !$acc end parallel loop
1056 :
1057 : !$acc parallel loop gang vector collapse(2) default(present)
1058 : do k = 1, nz, 1
1059 : do i = 1, ngrdcol
1060 : thlm_pert_1(i,k) = thlm(i,k) + sign_rtpthlp_zt(i,k) * clubb_params(i,iLscale_pert_coef) &
1061 : * sqrt( max( thlp2_zt(i,k), thl_tol**2 ) )
1062 : end do
1063 : end do
1064 : !$acc end parallel loop
1065 :
1066 : !$acc parallel loop gang vector default(present)
1067 : do i = 1, ngrdcol
1068 : mu_pert_1(i) = newmu(i) / clubb_params(i,iLscale_mu_coef)
1069 : end do
1070 : !$acc end parallel loop
1071 :
1072 : call compute_mixing_length( nz, ngrdcol, gr, thvm, thlm_pert_1, & ! In
1073 : rtm_pert_1, em, Lscale_max, p_in_Pa, & ! In
1074 : exner, thv_ds_zt, mu_pert_1, lmin, & ! In
1075 : saturation_formula, & ! In
1076 : l_implemented, & ! In
1077 : stats_metadata, & ! In
1078 : Lscale_pert_1, Lscale_up, Lscale_down ) ! Out
1079 :
1080 :
1081 : !$acc parallel loop gang vector collapse(2) default(present)
1082 : do k = 1, nz, 1
1083 : do i = 1, ngrdcol
1084 : rtm_pert_2(i,k) = rtm(i,k) - clubb_params(i,iLscale_pert_coef) &
1085 : * sqrt( max( rtp2_zt(i,k), rt_tol**2 ) )
1086 : end do
1087 : end do
1088 : !$acc end parallel loop
1089 :
1090 : !$acc parallel loop gang vector collapse(2) default(present)
1091 : do k = 1, nz, 1
1092 : do i = 1, ngrdcol
1093 : thlm_pert_2(i,k) = thlm(i,k) - sign_rtpthlp_zt(i,k) * clubb_params(i,iLscale_pert_coef) &
1094 : * sqrt( max( thlp2_zt(i,k), thl_tol**2 ) )
1095 : end do
1096 : end do
1097 : !$acc end parallel loop
1098 :
1099 : !$acc parallel loop gang vector default(present)
1100 : do i = 1, ngrdcol
1101 : mu_pert_2(i) = newmu(i) * clubb_params(i,iLscale_mu_coef)
1102 : end do
1103 : !$acc end parallel loop
1104 :
1105 : call compute_mixing_length( nz, ngrdcol, gr, thvm, thlm_pert_2, & ! In
1106 : rtm_pert_2, em, Lscale_max, p_in_Pa, & ! In
1107 : exner, thv_ds_zt, mu_pert_2, lmin, & ! In
1108 : saturation_formula, & ! In
1109 : l_implemented, & ! In
1110 : stats_metadata, & ! In
1111 : Lscale_pert_2, Lscale_up, Lscale_down ) ! Out
1112 :
1113 : else if ( l_avg_Lscale .and. l_Lscale_plume_centered ) then
1114 :
1115 : ! Take the values of thl and rt based one 1st or 2nd plume
1116 :
1117 : !$acc parallel loop gang vector collapse(2) default(present)
1118 : do k = 1, nz
1119 : do i = 1, ngrdcol
1120 : sign_rtpthlp_zt(i,k) = sign( one, rtpthlp_zt(i,k) )
1121 : end do
1122 : end do
1123 : !$acc end parallel loop
1124 :
1125 : !$acc parallel loop gang vector collapse(2) default(present)
1126 : do k = 1, nz
1127 : do i = 1, ngrdcol
1128 :
1129 : Lscale_pert_coef = clubb_params(i,iLscale_pert_coef)
1130 :
1131 : if ( pdf_params%rt_1(i,k) > pdf_params%rt_2(i,k) ) then
1132 :
1133 : rtm_pert_pos_rt(i,k) = pdf_params%rt_1(i,k) &
1134 : + Lscale_pert_coef * sqrt( max( pdf_params%varnce_rt_1(i,k), rt_tol**2 ) )
1135 :
1136 : thlm_pert_pos_rt(i,k) = pdf_params%thl_1(i,k) + ( sign_rtpthlp_zt(i,k) * Lscale_pert_coef &
1137 : * sqrt( max( pdf_params%varnce_thl_1(i,k), thl_tol**2 ) ) )
1138 :
1139 : thlm_pert_neg_rt(i,k) = pdf_params%thl_2(i,k) - ( sign_rtpthlp_zt(i,k) * Lscale_pert_coef &
1140 : * sqrt( max( pdf_params%varnce_thl_2(i,k), thl_tol**2 ) ) )
1141 :
1142 : rtm_pert_neg_rt(i,k) = pdf_params%rt_2(i,k) &
1143 : - Lscale_pert_coef * sqrt( max( pdf_params%varnce_rt_2(i,k), rt_tol**2 ) )
1144 :
1145 : !Lscale_weight = pdf_params%mixt_frac(i,k)
1146 :
1147 : else
1148 :
1149 : rtm_pert_pos_rt(i,k) = pdf_params%rt_2(i,k) &
1150 : + Lscale_pert_coef * sqrt( max( pdf_params%varnce_rt_2(i,k), rt_tol**2 ) )
1151 :
1152 : thlm_pert_pos_rt(i,k) = pdf_params%thl_2(i,k) + ( sign_rtpthlp_zt(i,k) * Lscale_pert_coef &
1153 : * sqrt( max( pdf_params%varnce_thl_2(i,k), thl_tol**2 ) ) )
1154 :
1155 : thlm_pert_neg_rt(i,k) = pdf_params%thl_1(i,k) - ( sign_rtpthlp_zt(i,k) * Lscale_pert_coef &
1156 : * sqrt( max( pdf_params%varnce_thl_1(i,k), thl_tol**2 ) ) )
1157 :
1158 : rtm_pert_neg_rt(i,k) = pdf_params%rt_1(i,k) &
1159 : - Lscale_pert_coef * sqrt( max( pdf_params%varnce_rt_1(i,k), rt_tol**2 ) )
1160 :
1161 : !Lscale_weight = 1.0_core_rknd - pdf_params%mixt_frac(i,k)
1162 :
1163 : end if
1164 :
1165 : end do
1166 : end do
1167 : !$acc end parallel loop
1168 :
1169 : !$acc parallel loop gang vector default(present)
1170 : do i = 1, ngrdcol
1171 : mu_pert_pos_rt(i) = newmu(i) / clubb_params(i,iLscale_mu_coef)
1172 : mu_pert_neg_rt(i) = newmu(i) * clubb_params(i,iLscale_mu_coef)
1173 : end do
1174 : !$acc end parallel loop
1175 :
1176 : ! Call length with perturbed values of thl and rt
1177 : call compute_mixing_length( nz, ngrdcol, gr, thvm, thlm_pert_pos_rt, & ! In
1178 : rtm_pert_pos_rt, em, Lscale_max, p_in_Pa, & ! In
1179 : exner, thv_ds_zt, mu_pert_pos_rt, lmin, & ! In
1180 : saturation_formula, & ! In
1181 : l_implemented, & ! In
1182 : stats_metadata, & ! In
1183 : Lscale_pert_1, Lscale_up, Lscale_down ) ! Out
1184 :
1185 : call compute_mixing_length( nz, ngrdcol, gr, thvm, thlm_pert_neg_rt, & ! In
1186 : rtm_pert_neg_rt, em, Lscale_max, p_in_Pa, & ! In
1187 : exner, thv_ds_zt, mu_pert_neg_rt, lmin, & ! In
1188 : saturation_formula, & ! In
1189 : l_implemented, & ! In
1190 : stats_metadata, & ! In
1191 : Lscale_pert_2, Lscale_up, Lscale_down ) ! Out
1192 : else
1193 : !$acc parallel loop gang vector collapse(2) default(present)
1194 768414816 : do k = 1, nz
1195 12690480816 : do i = 1, ngrdcol
1196 11922066000 : Lscale_pert_1(i,k) = unused_var ! Undefined
1197 12681545760 : Lscale_pert_2(i,k) = unused_var ! Undefined
1198 : end do
1199 : end do
1200 : !$acc end parallel loop
1201 : end if ! l_avg_Lscale
1202 :
1203 8935056 : if ( stats_metadata%l_stats_samp ) then
1204 : !$acc update host( Lscale_pert_1, Lscale_pert_2 )
1205 0 : do i = 1, ngrdcol
1206 0 : call stat_update_var( stats_metadata%iLscale_pert_1, Lscale_pert_1(i,:), & ! intent(in)
1207 0 : stats_zt(i) ) ! intent(inout)
1208 : call stat_update_var( stats_metadata%iLscale_pert_2, Lscale_pert_2(i,:), & ! intent(in)
1209 0 : stats_zt(i) ) ! intent(inout)
1210 : end do
1211 : end if ! stats_metadata%l_stats_samp
1212 :
1213 :
1214 : ! ********** NOTE: **********
1215 : ! This call to compute_mixing_length must be last. Otherwise, the values
1216 : ! of
1217 : ! Lscale_up and Lscale_down in stats will be based on perturbation length
1218 : ! scales
1219 : ! rather than the mean length scale.
1220 :
1221 : ! Diagnose CLUBB's turbulent mixing length scale.
1222 : call compute_mixing_length( nz, ngrdcol, gr, thvm, thlm, & ! In
1223 : rtm, em, Lscale_max, p_in_Pa, & ! In
1224 : exner, thv_ds_zt, newmu, lmin, & ! In
1225 : saturation_formula, & ! In
1226 : l_implemented, & ! In
1227 : stats_metadata, & ! In
1228 8935056 : Lscale, Lscale_up, Lscale_down ) ! Out
1229 :
1230 : if ( l_avg_Lscale ) then
1231 : if ( l_Lscale_plume_centered ) then
1232 : ! Weighted average of mean, pert_1, & pert_2
1233 : ! Lscale = 0.5_core_rknd * ( Lscale + Lscale_weight*Lscale_pert_1 &
1234 : ! + (1.0_core_rknd-Lscale_weight)*Lscale_pert_2
1235 : ! )
1236 : ! Weighted average of just the perturbed values
1237 : ! Lscale = Lscale_weight*Lscale_pert_1 +
1238 : ! (1.0_core_rknd-Lscale_weight)*Lscale_pert_2
1239 :
1240 : ! Un-weighted average of just the perturbed values
1241 : !$acc parallel loop gang vector collapse(2) default(present)
1242 : do k = 1, nz
1243 : do i = 1, ngrdcol
1244 : Lscale(i,k) = one_half *( Lscale_pert_1(i,k) + Lscale_pert_2(i,k) )
1245 : end do
1246 : end do
1247 : !$acc end parallel loop
1248 : else
1249 : !$acc parallel loop gang vector collapse(2) default(present)
1250 : do k = 1, nz
1251 : do i = 1, ngrdcol
1252 : Lscale(i,k) = one_third * ( Lscale(i,k) + Lscale_pert_1(i,k) + Lscale_pert_2(i,k) )
1253 : end do
1254 : end do
1255 : !$acc end parallel loop
1256 : end if
1257 : end if
1258 :
1259 : !$acc exit data delete( sign_rtpthlp_zt, Lscale_pert_1, Lscale_pert_2, &
1260 : !$acc thlm_pert_1, thlm_pert_2, rtm_pert_1, rtm_pert_2, &
1261 : !$acc thlm_pert_pos_rt, thlm_pert_neg_rt, rtm_pert_pos_rt, &
1262 : !$acc rtm_pert_neg_rt, &
1263 : !$acc mu_pert_1, mu_pert_2, mu_pert_pos_rt, mu_pert_neg_rt )
1264 :
1265 8935056 : return
1266 :
1267 : end subroutine calc_Lscale_directly
1268 :
1269 :
1270 :
1271 : !===============================================================================
1272 :
1273 0 : subroutine diagnose_Lscale_from_tau( nz, ngrdcol, gr, &
1274 0 : upwp_sfc, vpwp_sfc, um, vm, & !intent in
1275 0 : exner, p_in_Pa, & !intent in
1276 0 : rtm, thlm, thvm, & !intent in
1277 0 : rcm, ice_supersat_frac, &! intent in
1278 0 : em, sqrt_em_zt, & ! intent in
1279 : ufmin, tau_const, & ! intent in
1280 0 : sfc_elevation, Lscale_max, & ! intent in
1281 0 : clubb_params, & ! intent in
1282 : stats_metadata, & ! intent in
1283 : saturation_formula, & ! intent in
1284 : l_e3sm_config, & ! intent in
1285 : l_brunt_vaisala_freq_moist, & !intent in
1286 : l_use_thvm_in_bv_freq, &! intent in
1287 : l_smooth_Heaviside_tau_wpxp, & ! intent in
1288 : l_modify_limiters_for_cnvg_test, & ! intent in
1289 0 : stats_zm, & ! intent inout
1290 0 : brunt_vaisala_freq_sqd, brunt_vaisala_freq_sqd_mixed, & ! intent out
1291 0 : brunt_vaisala_freq_sqd_dry, brunt_vaisala_freq_sqd_moist, & ! intent out
1292 0 : Ri_zm, & ! intent out
1293 0 : invrs_tau_zt, invrs_tau_zm, & ! intent out
1294 0 : invrs_tau_sfc, invrs_tau_no_N2_zm, invrs_tau_bkgnd, & ! intent out
1295 0 : invrs_tau_shear, invrs_tau_N2_iso, & ! intent out
1296 0 : invrs_tau_wp2_zm, invrs_tau_xp2_zm, & ! intent out
1297 0 : invrs_tau_wp3_zm, invrs_tau_wp3_zt, invrs_tau_wpxp_zm, & ! intent out
1298 0 : tau_max_zm, tau_max_zt, tau_zm, tau_zt, & !intent out
1299 0 : Lscale, Lscale_up, Lscale_down)! intent out
1300 : ! Description:
1301 : ! Diagnose inverse damping time scales (invrs_tau_...) and turbulent mixing length (Lscale)
1302 : ! References:
1303 : ! Guo et al.(2021, JAMES)
1304 : !--------------------------------------------------------------------------------------------------
1305 :
1306 : use advance_helper_module, only: &
1307 : calc_brunt_vaisala_freq_sqd, &
1308 : smooth_heaviside_peskin, &
1309 : smooth_min, smooth_max
1310 :
1311 : use stats_type_utilities, only: &
1312 : stat_update_var ! Procedure
1313 :
1314 : use constants_clubb, only: &
1315 : one_fourth, &
1316 : one_half, &
1317 : vonk, &
1318 : zero, &
1319 : one, &
1320 : two, &
1321 : em_min, &
1322 : zero_threshold, &
1323 : eps
1324 :
1325 : use grid_class, only: &
1326 : grid, & ! Type
1327 : zt2zm, &
1328 : zm2zt, &
1329 : zm2zt2zm, &
1330 : zt2zm2zt, &
1331 : ddzt
1332 :
1333 : use clubb_precision, only: &
1334 : core_rknd
1335 :
1336 : use parameter_indices, only: &
1337 : nparams, & ! Variable(s)
1338 : iC_invrs_tau_bkgnd, &
1339 : iC_invrs_tau_shear, &
1340 : iC_invrs_tau_sfc, &
1341 : iC_invrs_tau_N2, &
1342 : iC_invrs_tau_N2_wp2 , &
1343 : iC_invrs_tau_N2_wpxp, &
1344 : iC_invrs_tau_N2_xp2, &
1345 : iC_invrs_tau_wpxp_N2_thresh, &
1346 : iC_invrs_tau_N2_clear_wp3, &
1347 : iC_invrs_tau_wpxp_Ri, &
1348 : ialtitude_threshold, &
1349 : ibv_efold, &
1350 : iwpxp_Ri_exp, &
1351 : iz_displace
1352 :
1353 : use error_code, only: &
1354 : err_code, &
1355 : clubb_fatal_error, &
1356 : clubb_at_least_debug_level
1357 :
1358 : use stats_variables, only: &
1359 : stats_metadata_type
1360 :
1361 : use stats_type, only: stats ! Type
1362 :
1363 : implicit none
1364 :
1365 : !--------------------------------- Input Variables ---------------------------------
1366 : integer, intent(in) :: &
1367 : nz, &
1368 : ngrdcol
1369 :
1370 : type (grid), target, intent(in) :: &
1371 : gr
1372 :
1373 : real( kind = core_rknd ), dimension(ngrdcol), intent(in) :: &
1374 : upwp_sfc, &
1375 : vpwp_sfc
1376 :
1377 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
1378 : um, &
1379 : vm, &
1380 : exner, &
1381 : p_in_Pa, &
1382 : rtm, &
1383 : thlm, &
1384 : thvm, &
1385 : rcm, &
1386 : ice_supersat_frac, &
1387 : em, &
1388 : sqrt_em_zt
1389 :
1390 : real(kind = core_rknd), intent(in) :: &
1391 : ufmin, &
1392 : tau_const
1393 :
1394 : real(kind = core_rknd), dimension(ngrdcol), intent(in) :: &
1395 : sfc_elevation, &
1396 : Lscale_max
1397 :
1398 : real( kind = core_rknd ), dimension(ngrdcol,nparams), intent(in) :: &
1399 : clubb_params ! Array of CLUBB's tunable parameters [units vary]
1400 :
1401 : type (stats_metadata_type), intent(in) :: &
1402 : stats_metadata
1403 :
1404 : integer, intent(in) :: &
1405 : saturation_formula ! Integer that stores the saturation formula to be used
1406 :
1407 : logical, intent(in) :: &
1408 : l_e3sm_config, &
1409 : l_brunt_vaisala_freq_moist, & ! Use a different formula for the Brunt-Vaisala frequency in
1410 : ! saturated atmospheres (from Durran and Klemp, 1982)
1411 : l_use_thvm_in_bv_freq, & ! Use thvm in the calculation of Brunt-Vaisala frequency
1412 : l_smooth_Heaviside_tau_wpxp ! Use the smoothed Heaviside 'Peskin' function
1413 : ! to compute invrs_tau_wpxp_zm
1414 :
1415 : ! Flag to activate modifications on limiters for convergence test
1416 : ! (smoothed max and min for Cx_fnc_Richardson in advance_helper_module.F90)
1417 : ! (remove the clippings on brunt_vaisala_freq_sqd_smth in mixing_length.F90)
1418 : ! (reduce threshold on limiters for Ri_zm in mixing_length.F90)
1419 : logical, intent(in) :: &
1420 : l_modify_limiters_for_cnvg_test
1421 :
1422 : !--------------------------- Input/Output Variables ---------------------------
1423 : type (stats), target, intent(inout), dimension(ngrdcol) :: &
1424 : stats_zm
1425 :
1426 : !--------------------------------- Output Variables ---------------------------------
1427 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
1428 : brunt_vaisala_freq_sqd, &
1429 : brunt_vaisala_freq_sqd_mixed, &
1430 : brunt_vaisala_freq_sqd_dry, &
1431 : brunt_vaisala_freq_sqd_moist, &
1432 : Ri_zm, &
1433 : invrs_tau_zt, &
1434 : invrs_tau_zm, &
1435 : invrs_tau_sfc, &
1436 : invrs_tau_no_N2_zm, &
1437 : invrs_tau_bkgnd, &
1438 : invrs_tau_shear, &
1439 : invrs_tau_N2_iso, &
1440 : invrs_tau_wp2_zm, &
1441 : invrs_tau_xp2_zm, &
1442 : invrs_tau_wp3_zm, &
1443 : invrs_tau_wp3_zt, &
1444 : invrs_tau_wpxp_zm, &
1445 : tau_max_zm, &
1446 : tau_max_zt, &
1447 : tau_zm, &
1448 : tau_zt, &
1449 : Lscale, &
1450 : Lscale_up, &
1451 : Lscale_down
1452 :
1453 : !--------------------------------- Local Variables ---------------------------------
1454 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
1455 0 : brunt_freq_pos, &
1456 0 : brunt_vaisala_freq_sqd_smth, & ! smoothed Buoyancy frequency squared, N^2 [s^-2]
1457 0 : brunt_freq_out_cloud, &
1458 0 : smooth_thlm, &
1459 0 : bvf_thresh, & ! temporatory array
1460 0 : H_invrs_tau_wpxp_N2 ! Heaviside function for clippings of invrs_tau_wpxp_N2
1461 :
1462 : real( kind = core_rknd ), dimension(ngrdcol) :: &
1463 0 : ustar
1464 :
1465 : real( kind = core_rknd ) :: &
1466 : C_invrs_tau_N2, &
1467 : C_invrs_tau_N2_wp2
1468 :
1469 : real( kind = core_rknd ), parameter :: &
1470 : min_max_smth_mag = 1.0e-9_core_rknd, & ! "base" smoothing magnitude before scaling
1471 : ! for the respective data structure. See
1472 : ! https://github.com/larson-group/clubb/issues/965#issuecomment-1119816722
1473 : ! for a plot on how output behaves with varying min_max_smth_mag
1474 : heaviside_smth_range = 1.0e-0_core_rknd ! range where Heaviside function is smoothed
1475 :
1476 : logical, parameter :: l_smooth_min_max = .false. ! whether to apply smooth min and max functions
1477 :
1478 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
1479 0 : ddzt_um, &
1480 0 : ddzt_vm, &
1481 0 : ddzt_umvm_sqd, &
1482 0 : ddzt_umvm_sqd_clipped, &
1483 0 : norm_ddzt_umvm, &
1484 0 : smooth_norm_ddzt_umvm, &
1485 0 : brunt_vaisala_freq_clipped, &
1486 0 : ice_supersat_frac_zm, &
1487 0 : invrs_tau_shear_smooth, &
1488 0 : Ri_zm_clipped, &
1489 0 : Ri_zm_smooth, &
1490 0 : em_clipped, &
1491 0 : tau_zm_unclipped, &
1492 0 : tau_zt_unclipped, &
1493 0 : tmp_calc, &
1494 0 : tmp_calc_max, &
1495 0 : tmp_calc_min_max
1496 :
1497 : real( kind = core_rknd ), dimension(ngrdcol) :: &
1498 0 : tmp_calc_ngrdcol
1499 :
1500 : integer :: i, k
1501 :
1502 : !--------------------------------- Begin Code ---------------------------------
1503 :
1504 : !$acc enter data create( brunt_freq_pos, brunt_vaisala_freq_sqd_smth, brunt_freq_out_cloud, &
1505 : !$acc smooth_thlm, bvf_thresh, H_invrs_tau_wpxp_N2, ustar, &
1506 : !$acc ddzt_um, ddzt_vm, norm_ddzt_umvm, smooth_norm_ddzt_umvm, &
1507 : !$acc brunt_vaisala_freq_clipped, &
1508 : !$acc ice_supersat_frac_zm, invrs_tau_shear_smooth, &
1509 : !$acc ddzt_umvm_sqd, tmp_calc_ngrdcol )
1510 :
1511 : !$acc enter data if( l_smooth_min_max .or. l_modify_limiters_for_cnvg_test ) &
1512 : !$acc create( Ri_zm_clipped, ddzt_umvm_sqd_clipped, &
1513 : !$acc tau_zm_unclipped, tau_zt_unclipped, Ri_zm_smooth, em_clipped, &
1514 : !$acc tmp_calc, tmp_calc_max, tmp_calc_min_max )
1515 :
1516 : !$acc parallel loop gang vector default(present)
1517 0 : do i = 1, ngrdcol
1518 0 : if ( gr%zm(i,1) - sfc_elevation(i) + clubb_params(i,iz_displace) < eps ) then
1519 0 : err_code = clubb_fatal_error
1520 : end if
1521 : end do
1522 : !$acc end parallel loop
1523 :
1524 0 : if ( clubb_at_least_debug_level(0) ) then
1525 0 : if ( err_code == clubb_fatal_error ) then
1526 0 : error stop "Lowest zm grid level is below ground in CLUBB."
1527 : end if
1528 : end if
1529 :
1530 : ! Smooth thlm by interpolating to zm then back to zt
1531 0 : smooth_thlm = zt2zm2zt( nz, ngrdcol, gr, thlm )
1532 :
1533 : call calc_brunt_vaisala_freq_sqd( nz, ngrdcol, gr, smooth_thlm, & ! intent(in)
1534 : exner, rtm, rcm, p_in_Pa, thvm, & ! intent(in)
1535 : ice_supersat_frac, & ! intent(in)
1536 : saturation_formula, & ! intent(in)
1537 : l_brunt_vaisala_freq_moist, & ! intent(in)
1538 : l_use_thvm_in_bv_freq, & ! intent(in)
1539 : clubb_params(:,ibv_efold), & ! intent(in)
1540 : brunt_vaisala_freq_sqd, & ! intent(out)
1541 : brunt_vaisala_freq_sqd_mixed,& ! intent(out)
1542 : brunt_vaisala_freq_sqd_dry, & ! intent(out)
1543 0 : brunt_vaisala_freq_sqd_moist ) ! intent(out)
1544 :
1545 : !$acc parallel loop gang vector default(present)
1546 0 : do i = 1, ngrdcol
1547 0 : tmp_calc_ngrdcol(i) = ( upwp_sfc(i)**2 + vpwp_sfc(i)**2 )**one_fourth
1548 : end do
1549 : !$acc end parallel loop
1550 :
1551 : if ( l_smooth_min_max ) then
1552 :
1553 : ! tmp_calc_ngrdcol used here because openacc has an issue with
1554 : ! a variable being both input and output of a function
1555 : ustar = smooth_max( ngrdcol, tmp_calc_ngrdcol, ufmin, min_max_smth_mag )
1556 :
1557 : else
1558 :
1559 : !$acc parallel loop gang vector default(present)
1560 0 : do i = 1, ngrdcol
1561 0 : ustar(i) = max( tmp_calc_ngrdcol(i), ufmin )
1562 : end do
1563 : !$acc end parallel loop
1564 :
1565 : end if
1566 :
1567 : !$acc parallel loop gang vector collapse(2) default(present)
1568 0 : do k = 1, nz
1569 0 : do i = 1, ngrdcol
1570 0 : invrs_tau_bkgnd(i,k) = clubb_params(i,iC_invrs_tau_bkgnd) / tau_const
1571 : end do
1572 : end do
1573 : !$acc end parallel loop
1574 :
1575 0 : ddzt_um = ddzt( nz, ngrdcol, gr, um )
1576 0 : ddzt_vm = ddzt( nz, ngrdcol, gr, vm )
1577 :
1578 : !$acc parallel loop gang vector collapse(2) default(present)
1579 0 : do k = 1, nz
1580 0 : do i = 1, ngrdcol
1581 0 : ddzt_umvm_sqd(i,k) = ddzt_um(i,k)**2 + ddzt_vm(i,k)**2
1582 0 : norm_ddzt_umvm(i,k) = sqrt( ddzt_umvm_sqd(i,k) )
1583 : end do
1584 : end do
1585 : !$acc end parallel loop
1586 :
1587 0 : smooth_norm_ddzt_umvm = zm2zt2zm( nz, ngrdcol, gr, norm_ddzt_umvm )
1588 :
1589 : !$acc parallel loop gang vector collapse(2) default(present)
1590 0 : do k = 1, nz
1591 0 : do i = 1, ngrdcol
1592 0 : invrs_tau_shear_smooth(i,k) = clubb_params(i,iC_invrs_tau_shear) * smooth_norm_ddzt_umvm(i,k)
1593 : end do
1594 : end do
1595 : !$acc end parallel loop
1596 :
1597 : ! Enforce that invrs_tau_shear is positive
1598 : invrs_tau_shear = smooth_max( nz, ngrdcol, invrs_tau_shear_smooth, &
1599 0 : zero_threshold, min_max_smth_mag )
1600 :
1601 : !$acc parallel loop gang vector collapse(2) default(present)
1602 0 : do k = 1, nz
1603 0 : do i = 1, ngrdcol
1604 0 : invrs_tau_sfc(i,k) = clubb_params(i,iC_invrs_tau_sfc) &
1605 : * ( ustar(i) / vonk ) &
1606 0 : / ( gr%zm(i,k) - sfc_elevation(i) + clubb_params(i,iz_displace) )
1607 : !C_invrs_tau_sfc * ( wp2 / vonk /ustar ) / ( gr%zm(1,:) -sfc_elevation + z_displace )
1608 : end do
1609 : end do
1610 : !$acc end parallel loop
1611 :
1612 : !$acc parallel loop gang vector collapse(2) default(present)
1613 0 : do k = 1, nz
1614 0 : do i = 1, ngrdcol
1615 0 : invrs_tau_no_N2_zm(i,k) = invrs_tau_bkgnd(i,k) + invrs_tau_sfc(i,k) + invrs_tau_shear(i,k)
1616 : end do
1617 : end do
1618 : !$acc end parallel loop
1619 :
1620 : !The min function below smooths the slope discontinuity in brunt freq
1621 : ! and thereby allows tau to remain large in Sc layers in which thlm may
1622 : ! be slightly stably stratified.
1623 0 : if ( l_modify_limiters_for_cnvg_test ) then
1624 :
1625 : !Remove the limiters to improve the solution convergence
1626 0 : brunt_vaisala_freq_sqd_smth = zm2zt2zm( nz,ngrdcol,gr, brunt_vaisala_freq_sqd_mixed )
1627 :
1628 : else ! default method
1629 :
1630 : if ( l_smooth_min_max ) then
1631 :
1632 : !$acc parallel loop gang vector collapse(2) default(present)
1633 : do k = 1, nz
1634 : do i = 1, ngrdcol
1635 : tmp_calc(i,k) = 1.e8_core_rknd * abs(brunt_vaisala_freq_sqd_mixed(i,k))**3
1636 : end do
1637 : end do
1638 : !$acc end parallel loop
1639 :
1640 : brunt_vaisala_freq_clipped = smooth_min( nz, ngrdcol, &
1641 : brunt_vaisala_freq_sqd_mixed, &
1642 : tmp_calc, &
1643 : 1.0e-4_core_rknd * min_max_smth_mag)
1644 :
1645 : brunt_vaisala_freq_sqd_smth = zm2zt2zm( nz, ngrdcol, gr, brunt_vaisala_freq_clipped )
1646 :
1647 : else
1648 :
1649 : !$acc parallel loop gang vector collapse(2) default(present)
1650 0 : do k = 1, nz
1651 0 : do i = 1, ngrdcol
1652 0 : brunt_vaisala_freq_clipped(i,k) = min( brunt_vaisala_freq_sqd_mixed(i,k), &
1653 0 : 1.e8_core_rknd * abs(brunt_vaisala_freq_sqd_mixed(i,k))**3)
1654 : end do
1655 : end do
1656 : !$acc end parallel loop
1657 :
1658 0 : brunt_vaisala_freq_sqd_smth = zm2zt2zm( nz, ngrdcol, gr, brunt_vaisala_freq_clipped )
1659 :
1660 : end if
1661 :
1662 : end if
1663 :
1664 0 : if ( stats_metadata%l_stats_samp ) then
1665 :
1666 : !$acc update host( brunt_vaisala_freq_sqd_smth )
1667 :
1668 0 : do i = 1, ngrdcol
1669 0 : call stat_update_var(stats_metadata%ibrunt_vaisala_freq_sqd_smth, brunt_vaisala_freq_sqd_smth(i,:), & ! intent(in)
1670 0 : stats_zm(i)) ! intent(inout)
1671 : end do
1672 : end if
1673 :
1674 0 : if ( l_modify_limiters_for_cnvg_test ) then
1675 :
1676 : !$acc parallel loop gang vector collapse(2) default(present)
1677 0 : do k = 1, nz
1678 0 : do i = 1, ngrdcol
1679 0 : Ri_zm_clipped(i,k) = max( 0.0_core_rknd, brunt_vaisala_freq_sqd_smth(i,k) ) &
1680 0 : / max( ddzt_umvm_sqd(i,k), 1.0e-12_core_rknd )
1681 : end do
1682 : end do
1683 : !$acc end parallel loop
1684 :
1685 0 : Ri_zm = zm2zt2zm( nz, ngrdcol, gr, Ri_zm_clipped )
1686 :
1687 : else ! default method
1688 :
1689 : if ( l_smooth_min_max ) then
1690 :
1691 : brunt_vaisala_freq_clipped = smooth_max( nz, ngrdcol, 1.0e-7_core_rknd, brunt_vaisala_freq_sqd_smth, &
1692 : 1.0e-4_core_rknd * min_max_smth_mag )
1693 :
1694 : ddzt_umvm_sqd_clipped = smooth_max( nz, ngrdcol, ddzt_umvm_sqd, 1.0e-7_core_rknd, &
1695 : 1.0e-6_core_rknd * min_max_smth_mag )
1696 :
1697 : !$acc parallel loop gang vector collapse(2) default(present)
1698 : do k = 1, nz
1699 : do i = 1, ngrdcol
1700 : Ri_zm(i,k) = brunt_vaisala_freq_clipped(i,k) / ddzt_umvm_sqd_clipped(i,k)
1701 : end do
1702 : end do
1703 : !$acc end parallel loop
1704 :
1705 : else
1706 :
1707 : !$acc parallel loop gang vector collapse(2) default(present)
1708 0 : do k = 1, nz
1709 0 : do i = 1, ngrdcol
1710 0 : Ri_zm(i,k) = max( 1.0e-7_core_rknd, brunt_vaisala_freq_sqd_smth(i,k) ) &
1711 0 : / max( ddzt_umvm_sqd(i,k), 1.0e-7_core_rknd )
1712 : end do
1713 : end do
1714 : !$acc end parallel loop
1715 :
1716 : end if
1717 :
1718 : end if
1719 :
1720 : if ( l_smooth_min_max ) then
1721 :
1722 : brunt_vaisala_freq_clipped = smooth_max( nz, ngrdcol, zero_threshold, &
1723 : brunt_vaisala_freq_sqd_smth, &
1724 : 1.0e-4_core_rknd * min_max_smth_mag )
1725 :
1726 : !$acc parallel loop gang vector collapse(2) default(present)
1727 : do k = 1, nz
1728 : do i = 1, ngrdcol
1729 : brunt_freq_pos(i,k) = sqrt( brunt_vaisala_freq_clipped(i,k) )
1730 : end do
1731 : end do
1732 : !$acc end parallel loop
1733 :
1734 : else
1735 :
1736 : !$acc parallel loop gang vector collapse(2) default(present)
1737 0 : do k = 1, nz
1738 0 : do i = 1, ngrdcol
1739 0 : brunt_freq_pos(i,k) = sqrt( max( zero_threshold, brunt_vaisala_freq_sqd_smth(i,k) ) )
1740 : end do
1741 : end do
1742 : !$acc end parallel loop
1743 :
1744 : end if
1745 :
1746 0 : ice_supersat_frac_zm = zt2zm( nz, ngrdcol, gr, ice_supersat_frac, zero_threshold )
1747 :
1748 : if ( l_smooth_min_max ) then
1749 :
1750 : ! roll this back as well once checks have passed
1751 : !$acc parallel loop gang vector collapse(2) default(present)
1752 : do k = 1, nz
1753 : do i = 1, ngrdcol
1754 : tmp_calc(i,k) = one - ice_supersat_frac_zm(i,k) / 0.001_core_rknd
1755 : end do
1756 : end do
1757 : !$acc end parallel loop
1758 :
1759 : tmp_calc_max = smooth_max( nz, ngrdcol, zero_threshold, tmp_calc, &
1760 : min_max_smth_mag)
1761 :
1762 : tmp_calc_min_max = smooth_min( nz, ngrdcol, one, tmp_calc_max, &
1763 : min_max_smth_mag )
1764 :
1765 : !$acc parallel loop gang vector collapse(2) default(present)
1766 : do k = 1, nz
1767 : do i = 1, ngrdcol
1768 : brunt_freq_out_cloud(i,k) = brunt_freq_pos(i,k) * tmp_calc_min_max(i,k)
1769 : end do
1770 : end do
1771 : !$acc end parallel loop
1772 :
1773 : else
1774 :
1775 : !$acc parallel loop gang vector collapse(2) default(present)
1776 0 : do k = 1, nz
1777 0 : do i = 1, ngrdcol
1778 0 : brunt_freq_out_cloud(i,k) &
1779 : = brunt_freq_pos(i,k) &
1780 : * min(one, max(zero_threshold, &
1781 0 : one - ( ( ice_supersat_frac_zm(i,k) / 0.001_core_rknd) )))
1782 : end do
1783 : end do
1784 : !$acc end parallel loop
1785 :
1786 : end if
1787 :
1788 : !$acc parallel loop gang vector collapse(2) default(present)
1789 0 : do k = 1, nz
1790 0 : do i = 1, ngrdcol
1791 0 : if ( gr%zm(i,k) < clubb_params(i,ialtitude_threshold) ) then
1792 0 : brunt_freq_out_cloud(i,k) = zero
1793 : end if
1794 : end do
1795 : end do
1796 : !$acc end parallel loop
1797 :
1798 : ! Write both bv extra terms to invrs_taus to disk
1799 0 : if ( stats_metadata%l_stats_samp ) then
1800 :
1801 : !$acc update host( brunt_freq_pos, brunt_freq_out_cloud )
1802 :
1803 0 : do i = 1, ngrdcol
1804 0 : call stat_update_var(stats_metadata%ibrunt_freq_pos, brunt_freq_pos(i,:), & ! intent(in)
1805 0 : stats_zm(i)) ! intent(inout)
1806 : call stat_update_var(stats_metadata%ibrunt_freq_out_cloud, brunt_freq_out_cloud(i,:), & ! intent(in)
1807 0 : stats_zm(i)) ! intent(inout)
1808 : end do
1809 : end if
1810 :
1811 : ! This time scale is used optionally for the return-to-isotropy term. It
1812 : ! omits invrs_tau_sfc based on the rationale that the isotropization
1813 : ! rate shouldn't be enhanced near the ground.
1814 : !$acc parallel loop gang vector collapse(2) default(present)
1815 0 : do k = 1, nz
1816 0 : do i = 1, ngrdcol
1817 :
1818 0 : C_invrs_tau_N2 = clubb_params(i,iC_invrs_tau_N2)
1819 0 : C_invrs_tau_N2_wp2 = clubb_params(i,iC_invrs_tau_N2_wp2)
1820 :
1821 0 : invrs_tau_N2_iso(i,k) = invrs_tau_bkgnd(i,k) + invrs_tau_shear(i,k) &
1822 0 : + C_invrs_tau_N2_wp2 * brunt_freq_pos(i,k)
1823 :
1824 : invrs_tau_wp2_zm(i,k) = invrs_tau_no_N2_zm(i,k) + &
1825 : C_invrs_tau_N2 * brunt_freq_pos(i,k) + &
1826 0 : C_invrs_tau_N2_wp2 * brunt_freq_out_cloud(i,k)
1827 :
1828 0 : invrs_tau_zm(i,k) = invrs_tau_no_N2_zm(i,k) + C_invrs_tau_N2 * brunt_freq_pos(i,k)
1829 : end do
1830 : end do
1831 : !$acc end parallel loop
1832 :
1833 :
1834 0 : if ( l_e3sm_config ) then
1835 :
1836 : !$acc parallel loop gang vector collapse(2) default(present)
1837 0 : do k = 1, nz
1838 0 : do i = 1, ngrdcol
1839 0 : invrs_tau_zm(i,k) = one_half * invrs_tau_zm(i,k)
1840 : end do
1841 : end do
1842 : !$acc end parallel loop
1843 :
1844 : !$acc parallel loop gang vector collapse(2) default(present)
1845 0 : do k = 1, nz
1846 0 : do i = 1, ngrdcol
1847 0 : invrs_tau_xp2_zm(i,k) = invrs_tau_no_N2_zm(i,k) &
1848 : + clubb_params(i,iC_invrs_tau_N2_xp2) * brunt_freq_pos(i,k) & ! 0
1849 : + clubb_params(i,iC_invrs_tau_sfc) * two &
1850 : * sqrt(em(i,k)) &
1851 0 : / ( gr%zm(i,k) - sfc_elevation(i) + clubb_params(i,iz_displace) ) ! small
1852 : end do
1853 : end do
1854 : !$acc end parallel loop
1855 :
1856 : if ( l_smooth_min_max ) then
1857 :
1858 : brunt_vaisala_freq_clipped = smooth_max( nz, ngrdcol, 1.0e-7_core_rknd, &
1859 : brunt_vaisala_freq_sqd_smth, &
1860 : 1.0e-4_core_rknd * min_max_smth_mag )
1861 :
1862 : !$acc parallel loop gang vector collapse(2) default(present)
1863 : do k = 1, nz
1864 : do i = 1, ngrdcol
1865 : tmp_calc(i,k) = sqrt( ddzt_umvm_sqd(i,k) / brunt_vaisala_freq_clipped(i,k) )
1866 : end do
1867 : end do
1868 : !$acc end parallel loop
1869 :
1870 : tmp_calc_max = smooth_max( nz, ngrdcol, tmp_calc, &
1871 : 0.3_core_rknd, 0.3_core_rknd * min_max_smth_mag )
1872 :
1873 : tmp_calc_min_max = smooth_min( nz, ngrdcol, tmp_calc_max, &
1874 : one, min_max_smth_mag )
1875 :
1876 : !$acc parallel loop gang vector collapse(2) default(present)
1877 : do k = 1, nz
1878 : do i = 1, ngrdcol
1879 : invrs_tau_xp2_zm(i,k) = tmp_calc_min_max(i,k) * invrs_tau_xp2_zm(i,k)
1880 : end do
1881 : end do
1882 : !$acc end parallel loop
1883 :
1884 : else
1885 :
1886 : !$acc parallel loop gang vector collapse(2) default(present)
1887 0 : do k = 1, nz
1888 0 : do i = 1, ngrdcol
1889 0 : invrs_tau_xp2_zm(i,k) &
1890 : = min( max( sqrt( ddzt_umvm_sqd(i,k) &
1891 : / max( 1.0e-7_core_rknd, brunt_vaisala_freq_sqd_smth(i,k) ) ), &
1892 : 0.3_core_rknd ), one ) &
1893 0 : * invrs_tau_xp2_zm(i,k)
1894 : end do
1895 : end do
1896 : !$acc end parallel loop
1897 :
1898 : end if
1899 :
1900 : !$acc parallel loop gang vector collapse(2) default(present)
1901 0 : do k = 1, nz
1902 0 : do i = 1, ngrdcol
1903 0 : invrs_tau_wpxp_zm(i,k) = two * invrs_tau_zm(i,k) &
1904 0 : + clubb_params(i,iC_invrs_tau_N2_wpxp) * brunt_freq_out_cloud(i,k)
1905 : end do
1906 : end do
1907 : !$acc end parallel loop
1908 :
1909 : else ! l_e3sm_config = false
1910 :
1911 : !$acc parallel loop gang vector collapse(2) default(present)
1912 0 : do k = 1, nz
1913 0 : do i = 1, ngrdcol
1914 0 : invrs_tau_xp2_zm(i,k) = invrs_tau_no_N2_zm(i,k) + &
1915 : clubb_params(i,iC_invrs_tau_N2) * brunt_freq_pos(i,k) + &
1916 0 : clubb_params(i,iC_invrs_tau_N2_xp2) * brunt_freq_out_cloud(i,k)
1917 : end do
1918 : end do
1919 : !$acc end parallel loop
1920 :
1921 0 : ice_supersat_frac_zm = zt2zm( nz, ngrdcol, gr, ice_supersat_frac, zero_threshold )
1922 :
1923 : ! !$acc parallel loop gang vector collapse(2) default(present)
1924 : ! do k = 1, nz
1925 : ! do i = 1, ngrdcol
1926 : ! if ( ice_supersat_frac_zm(i,k) <= 0.01_core_rknd &
1927 : ! .and. invrs_tau_xp2_zm(i,k) >= 0.003_core_rknd ) then
1928 : ! invrs_tau_xp2_zm(i,k) = 0.003_core_rknd
1929 : ! end if
1930 : ! end do
1931 : ! end do
1932 : ! !$acc end parallel loop
1933 :
1934 : !$acc parallel loop gang vector collapse(2) default(present)
1935 0 : do k = 1, nz
1936 0 : do i = 1, ngrdcol
1937 0 : invrs_tau_wpxp_zm(i,k) = invrs_tau_no_N2_zm(i,k) + &
1938 : clubb_params(i,iC_invrs_tau_N2) * brunt_freq_pos(i,k) + &
1939 0 : clubb_params(i,iC_invrs_tau_N2_wpxp) * brunt_freq_out_cloud(i,k)
1940 : end do
1941 : end do
1942 : !$acc end parallel loop
1943 :
1944 : end if ! l_e3sm_config
1945 :
1946 0 : if ( l_smooth_Heaviside_tau_wpxp ) then
1947 :
1948 : !$acc parallel loop gang vector collapse(2) default(present)
1949 0 : do k = 1, nz
1950 0 : do i = 1, ngrdcol
1951 0 : bvf_thresh(i,k) = brunt_vaisala_freq_sqd_smth(i,k) &
1952 0 : / clubb_params(i,iC_invrs_tau_wpxp_N2_thresh) - one
1953 : end do
1954 : end do
1955 : !$acc end parallel loop
1956 :
1957 0 : H_invrs_tau_wpxp_N2 = smooth_heaviside_peskin( nz, ngrdcol, bvf_thresh, heaviside_smth_range )
1958 :
1959 : else ! l_smooth_Heaviside_tau_wpxp = .false.
1960 :
1961 : !$acc parallel loop gang vector collapse(2) default(present)
1962 0 : do k = 1, nz
1963 0 : do i = 1, ngrdcol
1964 0 : if ( brunt_vaisala_freq_sqd_smth(i,k) > clubb_params(i,iC_invrs_tau_wpxp_N2_thresh) ) then
1965 0 : H_invrs_tau_wpxp_N2(i,k) = one
1966 : else
1967 0 : H_invrs_tau_wpxp_N2(i,k) = zero
1968 : end if
1969 : end do
1970 : end do
1971 : !$acc end parallel loop
1972 :
1973 : end if ! l_smooth_Heaviside_tau_wpxp
1974 :
1975 : if ( l_smooth_min_max ) then
1976 :
1977 : Ri_zm_smooth = smooth_max( nz, ngrdcol, Ri_zm, zero, &
1978 : 12.0_core_rknd * min_max_smth_mag )
1979 :
1980 : !$acc parallel loop gang vector collapse(2) default(present)
1981 : do k = 1, nz
1982 : do i = 1, ngrdcol
1983 : tmp_calc(i,k) = clubb_params(i,iC_invrs_tau_wpxp_Ri) * Ri_zm_smooth(i,k)**clubb_params(i,iwpxp_Ri_exp)
1984 : end do
1985 : end do
1986 : !$acc end parallel loop
1987 :
1988 : Ri_zm_smooth = smooth_min( nz, ngrdcol, tmp_calc, &
1989 : 12.0_core_rknd, 12.0_core_rknd * min_max_smth_mag )
1990 :
1991 : !$acc parallel loop gang vector collapse(2) default(present)
1992 : do k = 1, nz
1993 : do i = 1, ngrdcol
1994 :
1995 : if ( gr%zm(i,k) > clubb_params(i,ialtitude_threshold) ) then
1996 : invrs_tau_wpxp_zm(i,k) = invrs_tau_wpxp_zm(i,k) &
1997 : * ( one + H_invrs_tau_wpxp_N2(i,k) &
1998 : * Ri_zm_smooth(i,k) )
1999 :
2000 : end if
2001 : end do
2002 : end do
2003 : !$acc end parallel loop
2004 :
2005 : else ! l_smooth_min_max
2006 :
2007 : !$acc parallel loop gang vector collapse(2) default(present)
2008 0 : do k = 1, nz
2009 0 : do i = 1, ngrdcol
2010 0 : if ( gr%zm(i,k) > clubb_params(i,ialtitude_threshold) ) then
2011 0 : invrs_tau_wpxp_zm(i,k) = invrs_tau_wpxp_zm(i,k) &
2012 : * ( one + H_invrs_tau_wpxp_N2(i,k) &
2013 : * min( clubb_params(i,iC_invrs_tau_wpxp_Ri) &
2014 0 : * max( Ri_zm(i,k), zero)**clubb_params(i,iwpxp_Ri_exp), 12.0_core_rknd ) )
2015 : end if
2016 : end do
2017 : end do
2018 : !$acc end parallel loop
2019 :
2020 : end if
2021 :
2022 : !$acc parallel loop gang vector collapse(2) default(present)
2023 0 : do k = 1, nz
2024 0 : do i = 1, ngrdcol
2025 0 : invrs_tau_wp3_zm(i,k) = invrs_tau_wp2_zm(i,k) &
2026 0 : + clubb_params(i,iC_invrs_tau_N2_clear_wp3) * brunt_freq_out_cloud(i,k)
2027 : end do
2028 : end do
2029 : !$acc end parallel loop
2030 :
2031 : ! Calculate the maximum allowable value of time-scale tau,
2032 : ! which depends of the value of Lscale_max.
2033 : if ( l_smooth_min_max ) then
2034 :
2035 : em_clipped = smooth_max( nz, ngrdcol, em, em_min, min_max_smth_mag )
2036 :
2037 : !$acc parallel loop gang vector collapse(2) default(present)
2038 : do k = 1, nz
2039 : do i = 1, ngrdcol
2040 : tau_max_zt(i,k) = Lscale_max(i) / sqrt_em_zt(i,k)
2041 : tau_max_zm(i,k) = Lscale_max(i) / sqrt( em_clipped(i,k) )
2042 : end do
2043 : end do
2044 : !$acc end parallel loop
2045 :
2046 : else
2047 :
2048 : !$acc parallel loop gang vector collapse(2) default(present)
2049 0 : do k = 1, nz
2050 0 : do i = 1, ngrdcol
2051 0 : tau_max_zt(i,k) = Lscale_max(i) / sqrt_em_zt(i,k)
2052 0 : tau_max_zm(i,k) = Lscale_max(i) / sqrt( max( em(i,k), em_min ) )
2053 : end do
2054 : end do
2055 : !$acc end parallel loop
2056 :
2057 : end if
2058 :
2059 : if ( l_smooth_min_max ) then
2060 :
2061 : !$acc parallel loop gang vector collapse(2) default(present)
2062 : do k = 1, nz
2063 : do i = 1, ngrdcol
2064 : tau_zm_unclipped(i,k) = one / invrs_tau_zm(i,k)
2065 : end do
2066 : end do
2067 : !$acc end parallel loop
2068 :
2069 : tau_zm = smooth_min( nz, ngrdcol, tau_zm_unclipped, &
2070 : tau_max_zm, 1.0e3_core_rknd * min_max_smth_mag )
2071 :
2072 : tau_zt_unclipped = zm2zt( nz, ngrdcol, gr, tau_zm )
2073 :
2074 : tau_zt = smooth_min( nz, ngrdcol, tau_zt_unclipped, tau_max_zt, 1.0e3_core_rknd * min_max_smth_mag )
2075 :
2076 : else
2077 :
2078 : !$acc parallel loop gang vector collapse(2) default(present)
2079 0 : do k = 1, nz
2080 0 : do i = 1, ngrdcol
2081 0 : tau_zm(i,k) = min( one / invrs_tau_zm(i,k), tau_max_zm(i,k) )
2082 : end do
2083 : end do
2084 : !$acc end parallel loop
2085 :
2086 0 : tau_zt = zm2zt( nz, ngrdcol, gr, tau_zm )
2087 :
2088 : !$acc parallel loop gang vector collapse(2) default(present)
2089 0 : do k = 1, nz
2090 0 : do i = 1, ngrdcol
2091 0 : tau_zt(i,k) = min( tau_zt(i,k), tau_max_zt(i,k) )
2092 : end do
2093 : end do
2094 : !$acc end parallel loop
2095 :
2096 : end if
2097 :
2098 0 : invrs_tau_zt = zm2zt( nz, ngrdcol, gr, invrs_tau_zm )
2099 0 : invrs_tau_wp3_zt = zm2zt( nz, ngrdcol, gr, invrs_tau_wp3_zm )
2100 :
2101 : !$acc parallel loop gang vector collapse(2) default(present)
2102 0 : do k = 1, nz
2103 0 : do i = 1, ngrdcol
2104 :
2105 0 : Lscale(i,k) = tau_zt(i,k) * sqrt_em_zt(i,k)
2106 :
2107 : ! Lscale_up and Lscale_down aren't calculated with this option.
2108 : ! They are set to 0 for stats output.
2109 0 : Lscale_up(i,k) = zero
2110 0 : Lscale_down(i,k) = zero
2111 :
2112 : end do
2113 : end do
2114 : !$acc end parallel loop
2115 :
2116 : !$acc exit data delete( brunt_freq_pos, brunt_vaisala_freq_sqd_smth, brunt_freq_out_cloud, &
2117 : !$acc smooth_thlm, bvf_thresh, H_invrs_tau_wpxp_N2, ustar, &
2118 : !$acc ddzt_um, ddzt_vm, norm_ddzt_umvm, smooth_norm_ddzt_umvm, &
2119 : !$acc brunt_vaisala_freq_clipped, &
2120 : !$acc ice_supersat_frac_zm, invrs_tau_shear_smooth, &
2121 : !$acc ddzt_umvm_sqd, tmp_calc_ngrdcol )
2122 :
2123 : !$acc exit data if( l_smooth_min_max .or. l_modify_limiters_for_cnvg_test ) &
2124 : !$acc delete( Ri_zm_clipped, ddzt_umvm_sqd_clipped, &
2125 : !$acc tau_zm_unclipped, tau_zt_unclipped, Ri_zm_smooth, em_clipped, &
2126 : !$acc tmp_calc, tmp_calc_max, tmp_calc_min_max )
2127 :
2128 0 : return
2129 :
2130 : end subroutine diagnose_Lscale_from_tau
2131 :
2132 : end module mixing_length
|