Line data Source code
1 :
2 : module advance_helper_module
3 :
4 : ! Description:
5 : ! This module contains helper methods for the advance_* modules.
6 : !------------------------------------------------------------------------
7 :
8 : implicit none
9 :
10 : public :: &
11 : set_boundary_conditions_lhs, &
12 : set_boundary_conditions_rhs, &
13 : calc_stability_correction, &
14 : calc_brunt_vaisala_freq_sqd, &
15 : compute_Cx_fnc_Richardson, &
16 : wp2_term_splat_lhs, &
17 : wp3_term_splat_lhs, &
18 : smooth_min, smooth_max, &
19 : smooth_heaviside_peskin, &
20 : calc_xpwp, &
21 : vertical_avg, &
22 : vertical_integral, &
23 : Lscale_width_vert_avg
24 :
25 : interface calc_xpwp
26 : module procedure calc_xpwp_1D
27 : module procedure calc_xpwp_2D
28 : end interface
29 :
30 : private ! Set Default Scope
31 :
32 : !===============================================================================
33 : interface smooth_min
34 :
35 : ! These functions smooth the output of the min function
36 : ! by introducing a varyingly steep path between the two input variables.
37 : ! The degree to which smoothing is applied depends on the value of 'smth_coef'.
38 : ! If 'smth_coef' goes toward 0, the output of the min function will be
39 : ! 0.5 * ((a+b) - abs(a-b))
40 : ! If a > b, then this comes out to be b. Likewise, if a < b, abs(a-b)=b-a so we get a.
41 : ! Increasing the smoothing coefficient will lead to a greater degree of smoothing
42 : ! in the smooth min and max functions. Generally, the coefficient should roughly scale
43 : ! with the magnitude of data in the data structure that is to be smoothed, in order to
44 : ! obtain a sensible degree of smoothing (not too much, not too little).
45 :
46 : module procedure smooth_min_sclr_idx
47 : module procedure smooth_min_array_scalar
48 : module procedure smooth_min_arrays
49 : module procedure smooth_min_scalars
50 :
51 : end interface
52 :
53 : !===============================================================================
54 : interface smooth_max
55 :
56 : ! These functions smooth the output of the max functions
57 : ! by introducing a varyingly steep path between the two input variables.
58 : ! The degree to which smoothing is applied depends on the value of 'smth_coef'.
59 : ! If 'smth_coef' goes toward 0, the output of the max function will be
60 : ! 0.5 * ((a+b) + abs(a-b))
61 : ! If a > b, then this comes out to be a. Likewise, if a < b, abs(a-b)=b-a so we get b.
62 : ! Increasing the smoothing coefficient will lead to a greater degree of smoothing
63 : ! in the smooth min and max functions. Generally, the coefficient should roughly scale
64 : ! with the magnitude of data in the data structure that is to be smoothed, in order to
65 : ! obtain a sensible degree of smoothing (not too much, not too little).
66 :
67 : module procedure smooth_max_sclr_idx
68 : module procedure smooth_max_array_scalar
69 : module procedure smooth_max_arrays
70 : module procedure smooth_max_array_1D_scalar
71 : module procedure smooth_max_scalars
72 :
73 : end interface
74 :
75 : !===============================================================================
76 : contains
77 :
78 : !---------------------------------------------------------------------------
79 0 : subroutine set_boundary_conditions_lhs( diag_index, low_bound, high_bound, &
80 0 : lhs, &
81 : diag_index2, low_bound2, high_bound2 )
82 :
83 : ! Description:
84 : ! Sets the boundary conditions for a left-hand side LAPACK matrix.
85 : !
86 : ! References:
87 : ! none
88 : !---------------------------------------------------------------------------
89 :
90 : use clubb_precision, only: &
91 : core_rknd ! Variable(s)
92 :
93 : use constants_clubb, only: &
94 : one, zero
95 :
96 : implicit none
97 :
98 : ! Exernal
99 : intrinsic :: present
100 :
101 : ! Input Variables
102 : integer, intent(in) :: &
103 : diag_index, low_bound, high_bound ! boundary indexes for the first variable
104 :
105 : ! Input / Output Variables
106 : real( kind = core_rknd ), dimension(:,:), intent(inout) :: &
107 : lhs ! left hand side of the LAPACK matrix equation
108 :
109 : ! Optional Input Variables
110 : integer, intent(in), optional :: &
111 : diag_index2, low_bound2, high_bound2 ! boundary indexes for the second variable
112 :
113 : ! --------------------- BEGIN CODE ----------------------
114 :
115 0 : if ( ( present( low_bound2 ) .or. present( high_bound2 ) ) .and. &
116 : ( .not. present( diag_index2 ) ) ) then
117 :
118 0 : error stop "Boundary index provided without diag_index."
119 :
120 : end if
121 :
122 : ! Set the lower boundaries for the first variable
123 0 : lhs(:,low_bound) = zero
124 0 : lhs(diag_index,low_bound) = one
125 :
126 : ! Set the upper boundaries for the first variable
127 0 : lhs(:,high_bound) = zero
128 0 : lhs(diag_index,high_bound) = one
129 :
130 : ! Set the lower boundaries for the second variable, if it is provided
131 0 : if ( present( low_bound2 ) ) then
132 :
133 0 : lhs(:,low_bound2) = zero
134 0 : lhs(diag_index2,low_bound2) = one
135 :
136 : end if
137 :
138 : ! Set the upper boundaries for the second variable, if it is provided
139 0 : if ( present( high_bound2 ) ) then
140 :
141 0 : lhs(:,high_bound2) = zero
142 0 : lhs(diag_index2,high_bound2) = one
143 :
144 : end if
145 :
146 0 : return
147 : end subroutine set_boundary_conditions_lhs
148 :
149 : !--------------------------------------------------------------------------
150 0 : subroutine set_boundary_conditions_rhs( &
151 : low_value, low_bound, high_value, high_bound, &
152 0 : rhs, &
153 : low_value2, low_bound2, high_value2, high_bound2 )
154 :
155 : ! Description:
156 : ! Sets the boundary conditions for a right-hand side LAPACK vector.
157 : !
158 : ! References:
159 : ! none
160 : !---------------------------------------------------------------------------
161 :
162 : use clubb_precision, only: &
163 : core_rknd ! Variable(s)
164 :
165 : implicit none
166 :
167 : ! Exernal
168 : intrinsic :: present
169 :
170 : ! Input Variables
171 :
172 : ! The values for the first variable
173 : real( kind = core_rknd ), intent(in) :: low_value, high_value
174 :
175 : ! The bounds for the first variable
176 : integer, intent(in) :: low_bound, high_bound
177 :
178 : ! Input / Output Variables
179 :
180 : ! The right-hand side vector
181 : real( kind = core_rknd ), dimension(:), intent(inout) :: rhs
182 :
183 : ! Optional Input Variables
184 :
185 : ! The values for the second variable
186 : real( kind = core_rknd ), intent(in), optional :: low_value2, high_value2
187 :
188 : ! The bounds for the second variable
189 : integer, intent(in), optional :: low_bound2, high_bound2
190 :
191 :
192 : ! -------------------- BEGIN CODE ------------------------
193 :
194 : ! Stop execution if a boundary was provided without a value
195 0 : if ( (present( low_bound2 ) .and. (.not. present( low_value2 ))) .or. &
196 : (present( high_bound2 ) .and. (.not. present( high_value2 ))) ) then
197 :
198 0 : error stop "Boundary condition provided without value."
199 :
200 : end if
201 :
202 : ! Set the lower and upper bounds for the first variable
203 0 : rhs(low_bound) = low_value
204 0 : rhs(high_bound) = high_value
205 :
206 : ! If a lower bound was given for the second variable, set it
207 0 : if ( present( low_bound2 ) ) then
208 0 : rhs(low_bound2) = low_value2
209 : end if
210 :
211 : ! If an upper bound was given for the second variable, set it
212 0 : if ( present( high_bound2 ) ) then
213 0 : rhs(high_bound2) = high_value2
214 : end if
215 :
216 0 : return
217 : end subroutine set_boundary_conditions_rhs
218 :
219 : !===============================================================================
220 8935056 : subroutine calc_stability_correction( nz, ngrdcol, gr, &
221 8935056 : thlm, Lscale_zm, em, &
222 8935056 : exner, rtm, rcm, &
223 8935056 : p_in_Pa, thvm, ice_supersat_frac, &
224 8935056 : lambda0_stability_coef, &
225 8935056 : bv_efold, &
226 : saturation_formula, &
227 : l_brunt_vaisala_freq_moist, &
228 : l_use_thvm_in_bv_freq, &
229 8935056 : stability_correction )
230 : !
231 : ! Description:
232 : ! Stability Factor
233 : !
234 : ! References:
235 : !
236 : !--------------------------------------------------------------------
237 :
238 : use constants_clubb, only: &
239 : zero, one, three ! Constant(s)
240 :
241 : use grid_class, only: &
242 : grid ! Type
243 :
244 : use clubb_precision, only: &
245 : core_rknd ! Variable(s)
246 :
247 : implicit none
248 :
249 : ! ---------------- Input Variables ----------------
250 : integer, intent(in) :: &
251 : nz, &
252 : ngrdcol
253 :
254 : type (grid), target, intent(in) :: gr
255 :
256 : real( kind = core_rknd ), intent(in), dimension(ngrdcol,nz) :: &
257 : Lscale_zm, & ! Turbulent mixing length interp to m levs [m]
258 : em, & ! Turbulent Kinetic Energy (TKE) [m^2/s^2]
259 : thlm, & ! th_l (thermo. levels) [K]
260 : exner, & ! Exner function [-]
261 : rtm, & ! total water mixing ratio, r_t [kg/kg]
262 : rcm, & ! cloud water mixing ratio, r_c [kg/kg]
263 : p_in_Pa, & ! Air pressure [Pa]
264 : thvm, & ! Virtual potential temperature [K]
265 : ice_supersat_frac
266 :
267 : real( kind = core_rknd ), dimension(ngrdcol), intent(in) :: &
268 : lambda0_stability_coef, & ! CLUBB tunable parameter lambda0_stability_coef
269 : bv_efold ! Control parameter for inverse e-folding of
270 : ! cloud fraction in the mixed Brunt Vaisala frequency
271 :
272 : integer, intent(in) :: &
273 : saturation_formula ! Integer that stores the saturation formula to be used
274 :
275 : logical, intent(in) :: &
276 : l_brunt_vaisala_freq_moist, & ! Use a different formula for the Brunt-Vaisala frequency in
277 : ! saturated atmospheres (from Durran and Klemp, 1982)
278 : l_use_thvm_in_bv_freq ! Use thvm in the calculation of Brunt-Vaisala frequency
279 :
280 : ! ---------------- Output Variables ----------------
281 : real( kind = core_rknd ), intent(out), dimension(ngrdcol,nz) :: &
282 : stability_correction
283 :
284 : ! ---------------- Local Variables ----------------
285 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
286 17870112 : brunt_vaisala_freq_sqd, & ! []
287 17870112 : brunt_vaisala_freq_sqd_mixed, &
288 17870112 : brunt_vaisala_freq_sqd_dry, & ! []
289 17870112 : brunt_vaisala_freq_sqd_moist, &
290 17870112 : lambda0_stability
291 :
292 : integer :: i, k
293 :
294 : !------------ Begin Code --------------
295 :
296 : !$acc enter data create( brunt_vaisala_freq_sqd, brunt_vaisala_freq_sqd_mixed, &
297 : !$acc brunt_vaisala_freq_sqd_moist, brunt_vaisala_freq_sqd_dry, &
298 : !$acc lambda0_stability )
299 :
300 : call calc_brunt_vaisala_freq_sqd( nz, ngrdcol, gr, thlm, & ! intent(in)
301 : exner, rtm, rcm, p_in_Pa, thvm, & ! intent(in)
302 : ice_supersat_frac, & ! intent(in)
303 : saturation_formula, & ! intent(in)
304 : l_brunt_vaisala_freq_moist, & ! intent(in)
305 : l_use_thvm_in_bv_freq, & ! intent(in)
306 : bv_efold, & ! intent(in)
307 : brunt_vaisala_freq_sqd, & ! intent(out)
308 : brunt_vaisala_freq_sqd_mixed,& ! intent(out)
309 : brunt_vaisala_freq_sqd_dry, & ! intent(out)
310 8935056 : brunt_vaisala_freq_sqd_moist ) ! intent(out)
311 :
312 : !$acc parallel loop gang vector collapse(2) default(present)
313 768414816 : do k = 1, nz
314 12690480816 : do i = 1, ngrdcol
315 12681545760 : if ( brunt_vaisala_freq_sqd(i,k) > zero ) then
316 11650400594 : lambda0_stability(i,k) = lambda0_stability_coef(i)
317 : else
318 271665406 : lambda0_stability(i,k) = zero
319 : end if
320 : end do
321 : end do
322 : !$acc end parallel loop
323 :
324 : !$acc parallel loop gang vector collapse(2) default(present)
325 768414816 : do k = 1, nz
326 12690480816 : do i = 1, ngrdcol
327 23844132000 : stability_correction(i,k) = one + min( lambda0_stability(i,k) * brunt_vaisala_freq_sqd(i,k) &
328 36525677760 : * Lscale_zm(i,k)**2 / em(i,k), three )
329 : end do
330 : end do
331 : !$acc end parallel loop
332 :
333 : !$acc exit data delete( brunt_vaisala_freq_sqd, brunt_vaisala_freq_sqd_mixed, &
334 : !$acc brunt_vaisala_freq_sqd_moist, brunt_vaisala_freq_sqd_dry, &
335 : !$acc lambda0_stability )
336 :
337 8935056 : return
338 :
339 : end subroutine calc_stability_correction
340 :
341 : !===============================================================================
342 17870112 : subroutine calc_brunt_vaisala_freq_sqd( nz, ngrdcol, gr, thlm, &
343 17870112 : exner, rtm, rcm, p_in_Pa, thvm, &
344 17870112 : ice_supersat_frac, &
345 : saturation_formula, &
346 : l_brunt_vaisala_freq_moist, &
347 : l_use_thvm_in_bv_freq, &
348 17870112 : bv_efold, &
349 17870112 : brunt_vaisala_freq_sqd, &
350 17870112 : brunt_vaisala_freq_sqd_mixed,&
351 17870112 : brunt_vaisala_freq_sqd_dry, &
352 17870112 : brunt_vaisala_freq_sqd_moist )
353 :
354 : ! Description:
355 : ! Calculate the Brunt-Vaisala frequency squared, N^2.
356 :
357 : ! References:
358 : ! ?
359 : !-----------------------------------------------------------------------
360 :
361 : use clubb_precision, only: &
362 : core_rknd ! Konstant
363 :
364 : use constants_clubb, only: &
365 : grav, & ! Constant(s)
366 : Lv, &
367 : Cp, &
368 : Rd, &
369 : ep, &
370 : one, &
371 : one_half, &
372 : zero_threshold
373 :
374 : use parameters_model, only: &
375 : T0 ! Variable!
376 :
377 : use grid_class, only: &
378 : grid, & ! Type
379 : ddzt, & ! Procedure(s)
380 : zt2zm
381 :
382 : use T_in_K_module, only: &
383 : thlm2T_in_K ! Procedure
384 :
385 : use saturation, only: &
386 : sat_mixrat_liq ! Procedure
387 :
388 : implicit none
389 :
390 : !---------------------------- Input Variables ----------------------------
391 : integer, intent(in) :: &
392 : nz, &
393 : ngrdcol
394 :
395 : type (grid), target, intent(in) :: gr
396 :
397 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
398 : thlm, & ! th_l (thermo. levels) [K]
399 : exner, & ! Exner function [-]
400 : rtm, & ! total water mixing ratio, r_t [kg/kg]
401 : rcm, & ! cloud water mixing ratio, r_c [kg/kg]
402 : p_in_Pa, & ! Air pressure [Pa]
403 : thvm, & ! Virtual potential temperature [K]
404 : ice_supersat_frac
405 :
406 : integer, intent(in) :: &
407 : saturation_formula ! Integer that stores the saturation formula to be used
408 :
409 : logical, intent(in) :: &
410 : l_brunt_vaisala_freq_moist, & ! Use a different formula for the Brunt-Vaisala frequency in
411 : ! saturated atmospheres (from Durran and Klemp, 1982)
412 : l_use_thvm_in_bv_freq ! Use thvm in the calculation of Brunt-Vaisala frequency
413 :
414 : real( kind = core_rknd ), dimension(ngrdcol), intent(in) :: &
415 : bv_efold ! Control parameter for inverse e-folding of
416 : ! cloud fraction in the mixed Brunt Vaisala frequency
417 :
418 : !---------------------------- Output Variables ----------------------------
419 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
420 : brunt_vaisala_freq_sqd, & ! Brunt-Vaisala frequency squared, N^2 [1/s^2]
421 : brunt_vaisala_freq_sqd_mixed, &
422 : brunt_vaisala_freq_sqd_dry,&
423 : brunt_vaisala_freq_sqd_moist
424 :
425 : !---------------------------- Local Variables ----------------------------
426 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
427 35740224 : T_in_K, T_in_K_zm, rsat, rsat_zm, thm, thm_zm, ddzt_thlm, &
428 35740224 : ddzt_thm, ddzt_rsat, ddzt_rtm, thvm_zm, ddzt_thvm, ice_supersat_frac_zm
429 :
430 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
431 35740224 : stat_dry, stat_liq, ddzt_stat_liq, ddzt_stat_liq_zm, &
432 35740224 : stat_dry_virtual, stat_dry_virtual_zm, ddzt_rtm_zm
433 :
434 : integer :: i, k
435 :
436 : !---------------------------- Begin Code ----------------------------
437 :
438 : !$acc data copyin( gr, gr%zt, &
439 : !$acc thlm, exner, rtm, rcm, p_in_Pa, thvm, ice_supersat_frac, bv_efold ) &
440 : !$acc copyout( brunt_vaisala_freq_sqd, brunt_vaisala_freq_sqd_mixed, &
441 : !$acc brunt_vaisala_freq_sqd_dry, brunt_vaisala_freq_sqd_moist ) &
442 : !$acc create( T_in_K, T_in_K_zm, rsat, rsat_zm, thm, thm_zm, ddzt_thlm, &
443 : !$acc ddzt_thm, ddzt_rsat, ddzt_rtm, thvm_zm, ddzt_thvm, stat_dry, &
444 : !$acc stat_liq, ddzt_stat_liq, ddzt_stat_liq_zm, stat_dry_virtual, &
445 : !$acc stat_dry_virtual_zm, ddzt_rtm_zm, ice_supersat_frac_zm )
446 :
447 17870112 : ddzt_thlm = ddzt( nz, ngrdcol, gr, thlm )
448 17870112 : thvm_zm = zt2zm( nz, ngrdcol, gr, thvm, zero_threshold )
449 17870112 : ddzt_thvm = ddzt( nz, ngrdcol, gr, thvm )
450 :
451 : ! Dry Brunt-Vaisala frequency
452 17870112 : if ( l_use_thvm_in_bv_freq ) then
453 :
454 : !$acc parallel loop gang vector collapse(2) default(present)
455 0 : do k = 1, nz
456 0 : do i = 1, ngrdcol
457 0 : brunt_vaisala_freq_sqd(i,k) = ( grav / thvm_zm(i,k) ) * ddzt_thvm(i,k)
458 : end do
459 : end do
460 : !$acc end parallel loop
461 :
462 : else
463 :
464 : !$acc parallel loop gang vector collapse(2) default(present)
465 1536829632 : do k = 1, nz
466 25380961632 : do i = 1, ngrdcol
467 25363091520 : brunt_vaisala_freq_sqd(i,k) = ( grav / T0 ) * ddzt_thlm(i,k)
468 : end do
469 : end do
470 : !$acc end parallel loop
471 :
472 : end if
473 :
474 17870112 : T_in_K = thlm2T_in_K( nz, ngrdcol, thlm, exner, rcm )
475 17870112 : T_in_K_zm = zt2zm( nz, ngrdcol, gr, T_in_K, zero_threshold )
476 17870112 : rsat = sat_mixrat_liq( nz, ngrdcol, p_in_Pa, T_in_K, saturation_formula )
477 17870112 : rsat_zm = zt2zm( nz, ngrdcol, gr, rsat, zero_threshold )
478 17870112 : ddzt_rsat = ddzt( nz, ngrdcol, gr, rsat )
479 :
480 : !$acc parallel loop gang vector collapse(2) default(present)
481 1536829632 : do k = 1, nz
482 25380961632 : do i = 1, ngrdcol
483 25363091520 : thm(i,k) = thlm(i,k) + Lv/(Cp*exner(i,k)) * rcm(i,k)
484 : end do
485 : end do
486 : !$acc end parallel loop
487 :
488 17870112 : thm_zm = zt2zm( nz, ngrdcol, gr, thm, zero_threshold )
489 17870112 : ddzt_thm = ddzt( nz, ngrdcol, gr, thm )
490 17870112 : ddzt_rtm = ddzt( nz, ngrdcol, gr, rtm )
491 :
492 : !$acc parallel loop gang vector collapse(2) default(present)
493 1536829632 : do k = 1, nz
494 25380961632 : do i = 1, ngrdcol
495 23844132000 : stat_dry(i,k) = Cp * T_in_K(i,k) + grav * gr%zt(i,k)
496 25363091520 : stat_liq(i,k) = stat_dry(i,k) - Lv * rcm(i,k)
497 : end do
498 : end do
499 : !$acc end parallel loop
500 :
501 17870112 : ddzt_stat_liq = ddzt( nz, ngrdcol, gr, stat_liq )
502 17870112 : ddzt_stat_liq_zm = zt2zm( nz, ngrdcol, gr, ddzt_stat_liq)
503 :
504 : !$acc parallel loop gang vector collapse(2) default(present)
505 1536829632 : do k = 1, nz
506 25380961632 : do i = 1, ngrdcol
507 47688264000 : stat_dry_virtual(i,k) = stat_dry(i,k) + Cp * T_in_K(i,k) &
508 73051355520 : *( 0.608 * ( rtm(i,k) - rcm(i,k) )- rcm(i,k) )
509 : end do
510 : end do
511 : !$acc end parallel loop
512 :
513 17870112 : stat_dry_virtual_zm = zt2zm( nz, ngrdcol, gr, stat_dry_virtual, zero_threshold )
514 17870112 : ddzt_rtm_zm = zt2zm( nz, ngrdcol, gr, ddzt_rtm )
515 :
516 : !$acc parallel loop gang vector collapse(2) default(present)
517 1536829632 : do k = 1, nz
518 25380961632 : do i = 1, ngrdcol
519 25363091520 : brunt_vaisala_freq_sqd_dry(i,k) = ( grav / thm_zm(i,k) )* ddzt_thm(i,k)
520 : end do
521 : end do
522 : !$acc end parallel loop
523 :
524 : !$acc parallel loop gang vector collapse(2) default(present)
525 1536829632 : do k = 1, nz
526 25380961632 : do i = 1, ngrdcol
527 : ! In-cloud Brunt-Vaisala frequency. This is Eq. (36) of Durran and
528 : ! Klemp (1982)
529 47688264000 : brunt_vaisala_freq_sqd_moist(i,k) = &
530 : grav * ( ((one + Lv*rsat_zm(i,k) / (Rd*T_in_K_zm(i,k))) / &
531 : (one + ep*(Lv**2)*rsat_zm(i,k)/(Cp*Rd*T_in_K_zm(i,k)**2))) * &
532 : ( (one/thm_zm(i,k) * ddzt_thm(i,k)) + (Lv/(Cp*T_in_K_zm(i,k)))*ddzt_rsat(i,k)) - &
533 73051355520 : ddzt_rtm(i,k) )
534 : end do
535 : end do ! k=1, gr%nz
536 : !$acc end parallel loop
537 :
538 17870112 : ice_supersat_frac_zm = zt2zm( nz, ngrdcol, gr, ice_supersat_frac, zero_threshold )
539 :
540 : !$acc parallel loop gang vector collapse(2) default(present)
541 1536829632 : do k = 1, nz
542 25380961632 : do i = 1, ngrdcol
543 47688264000 : brunt_vaisala_freq_sqd_mixed(i,k) = &
544 : brunt_vaisala_freq_sqd_moist(i,k) + &
545 : exp( - bv_efold(i) * ice_supersat_frac_zm(i,k) ) * &
546 73051355520 : ( brunt_vaisala_freq_sqd_dry(i,k) - brunt_vaisala_freq_sqd_moist(i,k) )
547 : end do
548 : end do
549 : !$acc end parallel loop
550 :
551 17870112 : if ( l_brunt_vaisala_freq_moist ) then
552 :
553 0 : brunt_vaisala_freq_sqd = brunt_vaisala_freq_sqd_moist
554 :
555 : end if
556 :
557 : !$acc end data
558 :
559 17870112 : return
560 :
561 : end subroutine calc_brunt_vaisala_freq_sqd
562 :
563 : !===============================================================================
564 0 : subroutine compute_Cx_fnc_Richardson( nz, ngrdcol, gr, &
565 0 : thlm, um, vm, em, Lscale, exner, rtm, &
566 0 : rcm, p_in_Pa, thvm, rho_ds_zm, &
567 0 : ice_supersat_frac, &
568 0 : clubb_params, &
569 : saturation_formula, &
570 : l_brunt_vaisala_freq_moist, &
571 : l_use_thvm_in_bv_freq, &
572 : l_use_shear_Richardson, &
573 : l_modify_limiters_for_cnvg_test, &
574 : stats_metadata, &
575 0 : stats_zm, &
576 0 : Cx_fnc_Richardson )
577 :
578 : ! Description:
579 : ! Compute Cx as a function of the Richardson number
580 :
581 : ! References:
582 : ! cam:ticket:59
583 : !-----------------------------------------------------------------------
584 :
585 : use clubb_precision, only: &
586 : core_rknd ! Konstant
587 :
588 : use grid_class, only: &
589 : grid, & ! Type
590 : ddzt, & ! Procedure(s)
591 : zt2zm, &
592 : zm2zt2zm
593 :
594 : use constants_clubb, only: &
595 : one, &
596 : zero, &
597 : zero_threshold
598 :
599 : use interpolation, only: &
600 : linear_interp_factor ! Procedure
601 :
602 : use parameter_indices, only: &
603 : nparams, & ! Variable(s)
604 : iCx_min, &
605 : iCx_max, &
606 : iRichardson_num_min, &
607 : iRichardson_num_max, &
608 : ibv_efold
609 :
610 : use stats_variables, only: &
611 : stats_metadata_type
612 :
613 : use stats_type_utilities, only: &
614 : stat_update_var ! Procedure
615 :
616 : use stats_type, only: stats ! Type
617 :
618 : implicit none
619 :
620 : !------------------------------ Constant Parameters ------------------------------
621 : real( kind = core_rknd ), parameter :: &
622 : Richardson_num_divisor_threshold = 1.0e-6_core_rknd, &
623 : Cx_fnc_Richardson_below_ground_value = one
624 :
625 : logical, parameter :: &
626 : l_Cx_fnc_Richardson_vert_avg = .false. ! Vertically average Cx_fnc_Richardson over a
627 : ! distance of Lscale
628 :
629 : real( kind = core_rknd ), parameter :: &
630 : min_max_smth_mag = 1.0e-9_core_rknd ! "base" smoothing magnitude before scaling
631 : ! for the respective data structure. See
632 : ! https://github.com/larson-group/clubb/issues/965#issuecomment-1119816722
633 : ! for a plot on how output behaves with varying min_max_smth_mag
634 :
635 : !------------------------------ Input Variables ------------------------------
636 : integer, intent(in) :: &
637 : nz, &
638 : ngrdcol
639 :
640 : type (grid), target, intent(in) :: gr
641 :
642 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
643 : thlm, & ! th_l (liquid water potential temperature) [K]
644 : um, & ! u mean wind component (thermodynamic levels) [m/s]
645 : vm, & ! v mean wind component (thermodynamic levels) [m/s]
646 : em, & ! Turbulent Kinetic Energy (TKE) [m^2/s^2]
647 : Lscale, & ! Turbulent mixing length [m]
648 : exner, & ! Exner function [-]
649 : rtm, & ! total water mixing ratio, r_t [kg/kg]
650 : rcm, & ! cloud water mixing ratio, r_c [kg/kg]
651 : p_in_Pa, & ! Air pressure [Pa]
652 : thvm, & ! Virtual potential temperature [K]
653 : rho_ds_zm, & ! Dry static density on momentum levels [kg/m^3]
654 : ice_supersat_frac ! ice cloud fraction
655 :
656 : real( kind = core_rknd ), dimension(ngrdcol,nparams), intent(in) :: &
657 : clubb_params ! Array of CLUBB's tunable parameters [units vary]
658 :
659 : integer, intent(in) :: &
660 : saturation_formula ! Integer that stores the saturation formula to be used
661 :
662 : logical, intent(in) :: &
663 : l_brunt_vaisala_freq_moist, & ! Use a different formula for the Brunt-Vaisala frequency in
664 : ! saturated atmospheres (from Durran and Klemp, 1982)
665 : l_use_thvm_in_bv_freq, & ! Use thvm in the calculation of Brunt-Vaisala frequency
666 : l_use_shear_Richardson ! Use shear in the calculation of Richardson number
667 :
668 : ! Flag to activate modifications on limiters for convergence test
669 : ! (smoothed max and min for Cx_fnc_Richardson in advance_helper_module.F90)
670 : ! (remove the clippings on brunt_vaisala_freq_sqd_smth in mixing_length.F90)
671 : ! (reduce threshold on limiters for sqrt_Ri_zm in mixing_length.F90)
672 : logical, intent(in) :: &
673 : l_modify_limiters_for_cnvg_test
674 :
675 : type (stats_metadata_type), intent(in) :: &
676 : stats_metadata
677 :
678 : !------------------------------ InOut Variable ------------------------------
679 : type (stats), target, dimension(ngrdcol), intent(inout) :: &
680 : stats_zm
681 :
682 : !------------------------------ Output Variable ------------------------------
683 : real( kind = core_rknd), dimension(ngrdcol,nz), intent(out) :: &
684 : Cx_fnc_Richardson
685 :
686 : !------------------------------ Local Variables ------------------------------
687 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
688 0 : brunt_vaisala_freq_sqd, &
689 0 : brunt_vaisala_freq_sqd_mixed,&
690 0 : brunt_vaisala_freq_sqd_dry, &
691 0 : brunt_vaisala_freq_sqd_moist, &
692 0 : fnc_Richardson, &
693 0 : fnc_Richardson_clipped, &
694 0 : fnc_Richardson_smooth, &
695 0 : Ri_zm, &
696 0 : ddzt_um, &
697 0 : ddzt_vm, &
698 0 : shear_sqd, &
699 0 : Lscale_zm, &
700 0 : Cx_fnc_interp, &
701 0 : Cx_fnc_Richardson_avg
702 :
703 : real ( kind = core_rknd ) :: &
704 : invrs_min_max_diff, &
705 : invrs_num_div_thresh
706 :
707 : real( kind = core_rknd ) :: &
708 : Richardson_num_max, & ! CLUBB tunable parameter Richardson_num_max
709 : Richardson_num_min, & ! CLUBB tunable parameter Richardson_num_min
710 : Cx_max, & ! CLUBB tunable parameter max of Cx_fnc_Richardson
711 : Cx_min ! CLUBB tunable parameter min of Cx_fnc_Richardson
712 :
713 : integer :: smth_type = 1
714 :
715 : integer :: i, k
716 :
717 : !------------------------------ Begin Code ------------------------------
718 :
719 : !$acc enter data create( brunt_vaisala_freq_sqd, brunt_vaisala_freq_sqd_mixed, &
720 : !$acc brunt_vaisala_freq_sqd_dry, brunt_vaisala_freq_sqd_moist, &
721 : !$acc Cx_fnc_interp, &
722 : !$acc Ri_zm, ddzt_um, ddzt_vm, shear_sqd, Lscale_zm, &
723 : !$acc Cx_fnc_Richardson_avg, fnc_Richardson, &
724 : !$acc fnc_Richardson_clipped, fnc_Richardson_smooth )
725 :
726 : call calc_brunt_vaisala_freq_sqd( nz, ngrdcol, gr, thlm, & ! intent(in)
727 : exner, rtm, rcm, p_in_Pa, thvm, & ! intent(in)
728 : ice_supersat_frac, & ! intent(in)
729 : saturation_formula, & ! intent(in)
730 : l_brunt_vaisala_freq_moist, & ! intent(in)
731 : l_use_thvm_in_bv_freq, & ! intent(in)
732 : clubb_params(:,ibv_efold), & ! intent(in)
733 : brunt_vaisala_freq_sqd, & ! intent(out)
734 : brunt_vaisala_freq_sqd_mixed,& ! intent(out)
735 : brunt_vaisala_freq_sqd_dry, & ! intent(out)
736 0 : brunt_vaisala_freq_sqd_moist ) ! intent(out)
737 :
738 0 : invrs_num_div_thresh = one / Richardson_num_divisor_threshold
739 :
740 0 : Lscale_zm = zt2zm( nz, ngrdcol, gr, Lscale, zero_threshold )
741 :
742 : ! Calculate shear_sqd
743 0 : ddzt_um = ddzt( nz, ngrdcol, gr, um )
744 0 : ddzt_vm = ddzt( nz, ngrdcol, gr, vm )
745 :
746 : !$acc parallel loop gang vector collapse(2) default(present)
747 0 : do k = 1, nz
748 0 : do i = 1, ngrdcol
749 0 : shear_sqd(i,k) = ddzt_um(i,k)**2 + ddzt_vm(i,k)**2
750 : end do
751 : end do
752 : !$acc end parallel loop
753 :
754 0 : if ( stats_metadata%l_stats_samp ) then
755 : !$acc update host(shear_sqd)
756 0 : do i = 1, ngrdcol
757 0 : call stat_update_var( stats_metadata%ishear_sqd, shear_sqd(i,:), & ! intent(in)
758 0 : stats_zm(i) ) ! intent(inout)
759 : end do
760 : end if
761 :
762 0 : if ( l_use_shear_Richardson ) then
763 :
764 : !$acc parallel loop gang vector collapse(2) default(present)
765 0 : do k = 1, nz
766 0 : do i = 1, ngrdcol
767 0 : Ri_zm(i,k) = max( 1.0e-7_core_rknd, brunt_vaisala_freq_sqd_mixed(i,k) ) &
768 0 : / max( shear_sqd(i,k), 1.0e-7_core_rknd )
769 : end do
770 : end do
771 : !$acc end parallel loop
772 :
773 : else
774 : !$acc parallel loop gang vector collapse(2) default(present)
775 0 : do k = 1, nz
776 0 : do i = 1, ngrdcol
777 0 : Ri_zm(i,k) = brunt_vaisala_freq_sqd(i,k) * invrs_num_div_thresh
778 : end do
779 : end do
780 : !$acc end parallel loop
781 : end if
782 :
783 : ! Cx_fnc_Richardson is interpolated based on the value of Richardson_num
784 : ! The min function ensures that Cx does not exceed Cx_max, regardless of the
785 : ! value of Richardson_num_max.
786 0 : if ( l_modify_limiters_for_cnvg_test ) then
787 :
788 : !$acc parallel loop gang vector collapse(2) default(present)
789 0 : do k = 1, nz
790 0 : do i = 1, ngrdcol
791 :
792 0 : Richardson_num_max = clubb_params(i,iRichardson_num_max)
793 0 : Richardson_num_min = clubb_params(i,iRichardson_num_min)
794 :
795 0 : invrs_min_max_diff = one / ( Richardson_num_max - Richardson_num_min )
796 :
797 0 : fnc_Richardson(i,k) = ( Ri_zm(i,k) - clubb_params(i,iRichardson_num_min) ) * invrs_min_max_diff
798 : end do
799 : end do
800 :
801 : fnc_Richardson_clipped = smooth_min( nz, ngrdcol, one, &
802 : fnc_Richardson, &
803 0 : min_max_smth_mag )
804 :
805 : fnc_Richardson_smooth = smooth_max( nz, ngrdcol, zero, &
806 : fnc_Richardson_clipped, &
807 0 : min_max_smth_mag )
808 :
809 : ! use smoothed max amd min to achive smoothed profile and avoid discontinuities
810 : !$acc parallel loop gang vector collapse(2) default(present)
811 0 : do k = 1, nz
812 0 : do i = 1, ngrdcol
813 :
814 0 : Cx_max = clubb_params(i,iCx_max)
815 0 : Cx_min = clubb_params(i,iCx_min)
816 :
817 0 : Cx_fnc_interp(i,k) = fnc_Richardson_smooth(i,k) * ( Cx_max - Cx_min ) + Cx_min
818 : end do
819 : end do
820 :
821 0 : Cx_fnc_Richardson = zm2zt2zm( nz, ngrdcol, gr, Cx_fnc_interp )
822 :
823 : else ! default method
824 :
825 : !$acc parallel loop gang vector collapse(2) default(present)
826 0 : do k = 1, nz
827 0 : do i = 1, ngrdcol
828 :
829 0 : Richardson_num_max = clubb_params(i,iRichardson_num_max)
830 0 : Richardson_num_min = clubb_params(i,iRichardson_num_min)
831 0 : Cx_max = clubb_params(i,iCx_max)
832 0 : Cx_min = clubb_params(i,iCx_min)
833 :
834 0 : invrs_min_max_diff = one / ( Richardson_num_max - Richardson_num_min )
835 :
836 0 : Cx_fnc_Richardson(i,k) = ( max(min(Richardson_num_max, Ri_zm(i,k)), Richardson_num_min) &
837 : - Richardson_num_min ) &
838 0 : * invrs_min_max_diff * ( Cx_max - Cx_min ) + Cx_min
839 : end do
840 : end do
841 : !$acc end parallel loop
842 :
843 : end if
844 :
845 : if ( l_Cx_fnc_Richardson_vert_avg ) then
846 : Cx_fnc_Richardson = Lscale_width_vert_avg( nz, ngrdcol, gr, smth_type, &
847 : Cx_fnc_Richardson, Lscale_zm, rho_ds_zm, &
848 : Cx_fnc_Richardson_below_ground_value )
849 :
850 : !$acc parallel loop gang vector collapse(2) default(present)
851 : do k = 1, nz
852 : do i = 1, ngrdcol
853 : Cx_fnc_Richardson(i,k) = Cx_fnc_Richardson_avg(i,k)
854 : end do
855 : end do
856 : !$acc end parallel loop
857 : end if
858 :
859 : ! On some compilers, roundoff error can result in Cx_fnc_Richardson being
860 : ! slightly outside the range [0,1]. Thus, it is clipped here.
861 : !$acc parallel loop gang vector collapse(2) default(present)
862 0 : do k = 1, nz
863 0 : do i = 1, ngrdcol
864 0 : Cx_fnc_Richardson(i,k) = max( zero, min( one, Cx_fnc_Richardson(i,k) ) )
865 : end do
866 : end do
867 : !$acc end parallel loop
868 :
869 : !$acc exit data delete( brunt_vaisala_freq_sqd, brunt_vaisala_freq_sqd_mixed, &
870 : !$acc brunt_vaisala_freq_sqd_dry, brunt_vaisala_freq_sqd_moist, &
871 : !$acc Cx_fnc_interp, Ri_zm, &
872 : !$acc ddzt_um, ddzt_vm, shear_sqd, Lscale_zm, &
873 : !$acc Cx_fnc_Richardson_avg, fnc_Richardson, &
874 : !$acc fnc_Richardson_clipped, fnc_Richardson_smooth )
875 :
876 0 : return
877 :
878 : end subroutine compute_Cx_fnc_Richardson
879 : !----------------------------------------------------------------------
880 :
881 : !----------------------------------------------------------------------
882 8935056 : function Lscale_width_vert_avg( nz, ngrdcol, gr, smth_type, &
883 8935056 : var_profile, Lscale_zm, rho_ds_zm, &
884 : var_below_ground_value )&
885 35740224 : result (Lscale_width_vert_avg_output)
886 :
887 : ! Description:
888 : ! Averages a profile with a running mean of width Lscale_zm
889 :
890 : ! References:
891 : ! cam:ticket:59
892 :
893 : use clubb_precision, only: &
894 : core_rknd ! Precision
895 :
896 : use grid_class, only: &
897 : grid ! Type
898 :
899 : use constants_clubb, only: &
900 : zero
901 :
902 : implicit none
903 :
904 : !-------------------------- Input Variables --------------------------
905 : integer, intent(in) :: &
906 : nz, &
907 : ngrdcol, &
908 : smth_type
909 :
910 : type (grid), target, intent(in) :: gr
911 :
912 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
913 : var_profile, & ! Profile on momentum levels
914 : Lscale_zm, & ! Lscale on momentum levels
915 : rho_ds_zm ! Dry static energy on momentum levels!
916 :
917 : real( kind = core_rknd ), intent(in) :: &
918 : var_below_ground_value ! Value to use below ground
919 :
920 : !-------------------------- Result Variable --------------------------
921 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
922 : Lscale_width_vert_avg_output ! Vertically averaged profile (on momentum levels)
923 :
924 : !-------------------------- Local Variables --------------------------
925 : integer :: &
926 : k, i, & ! Loop variable
927 : k_avg_lower, &
928 : k_avg_upper, &
929 : k_avg
930 :
931 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
932 17870112 : one_half_avg_width, &
933 17870112 : numer_terms, &
934 8935056 : denom_terms
935 :
936 : integer :: &
937 : n_below_ground_levels
938 :
939 : real( kind = core_rknd ) :: &
940 : numer_integral, & ! Integral in the numerator (see description)
941 : denom_integral ! Integral in the denominator (see description)
942 :
943 : !-------------------------- Begin Code --------------------------
944 :
945 : !$acc enter data create( one_half_avg_width, numer_terms, denom_terms )
946 :
947 8935056 : if ( smth_type == 1 ) then
948 : !$acc parallel loop gang vector collapse(2) default(present)
949 0 : do k = 1, nz
950 0 : do i = 1, ngrdcol
951 0 : one_half_avg_width(i,k) = max( Lscale_zm(i,k), 500.0_core_rknd )
952 : end do
953 : end do
954 8935056 : else if (smth_type == 2 ) then
955 : !$acc parallel loop gang vector collapse(2) default(present)
956 768414816 : do k = 1, nz
957 12690480816 : do i = 1, ngrdcol
958 12681545760 : one_half_avg_width(i,k) = 60.0_core_rknd
959 : end do
960 : end do
961 : endif
962 :
963 : ! Pre calculate numerator and denominator terms
964 : !$acc parallel loop gang vector collapse(2) default(present)
965 768414816 : do k = 1, nz
966 12690480816 : do i = 1, ngrdcol
967 11922066000 : numer_terms(i,k) = rho_ds_zm(i,k) * gr%dzm(i,k) * var_profile(i,k)
968 12681545760 : denom_terms(i,k) = rho_ds_zm(i,k) * gr%dzm(i,k)
969 : end do
970 : end do
971 :
972 : ! For every grid level
973 : !$acc parallel loop gang vector collapse(2) default(present)
974 768414816 : do k = 1, nz
975 12690480816 : do i = 1, ngrdcol
976 :
977 : !-----------------------------------------------------------------------
978 : ! Hunt down all vertical levels with one_half_avg_width(k) of gr%zm(k).
979 : !
980 : ! Note: Outdated explanation of version that improves CPU performance
981 : ! below. Reworked due to it requiring a k dependency. Now we
982 : ! begin looking for k_avg_upper and k_avg_lower starting at
983 : ! the kth level.
984 : !
985 : ! Outdated but potentially useful note:
986 : ! k_avg_upper and k_avg_lower can be saved each loop iteration, this
987 : ! reduces iterations beacuse the kth values are likely to be within
988 : ! one or two grid levels of the k-1th values. Less searching is required
989 : ! by starting the search at the previous values and incrementing or
990 : ! decrement as needed.
991 : !-----------------------------------------------------------------------
992 :
993 : ! Determine if k_avg_upper needs to increment
994 12081472528 : do k_avg_upper = k, nz-1
995 12081472528 : if ( gr%zm(i,k_avg_upper+1) - gr%zm(i,k) > one_half_avg_width(i,k) ) then
996 : exit
997 : end if
998 : end do
999 :
1000 : ! Determine if k_avg_lower needs to decrement
1001 12081472528 : do k_avg_lower = k, 2, -1
1002 12081472528 : if ( gr%zm(i,k) - gr%zm(i,k_avg_lower-1) > one_half_avg_width(i,k) ) then
1003 : exit
1004 : end if
1005 : end do
1006 :
1007 : ! Compute the number of levels below ground to include.
1008 11922066000 : if ( k_avg_lower > 1 ) then
1009 :
1010 : ! k=1, the lowest "real" level, is not included in the average, so no
1011 : ! below-ground levels should be included.
1012 24162945056 : n_below_ground_levels = 0
1013 :
1014 : numer_integral = zero
1015 : denom_integral = zero
1016 :
1017 : else
1018 :
1019 : ! The number of below-ground levels included is equal to the distance
1020 : ! below the lowest level spanned by one_half_avg_width(k)
1021 : ! divided by the distance between vertical levels below ground; the
1022 : ! latter is assumed to be the same as the distance between the first and
1023 : ! second vertical levels.
1024 841557600 : n_below_ground_levels = int( ( one_half_avg_width(i,k)-(gr%zm(i,k)-gr%zm(i,1)) ) / &
1025 1122076800 : ( gr%zm(i,2)-gr%zm(i,1) ) )
1026 :
1027 280519200 : numer_integral = n_below_ground_levels * denom_terms(i,1) * var_below_ground_value
1028 : denom_integral = n_below_ground_levels * denom_terms(i,1)
1029 :
1030 : end if
1031 :
1032 : ! Add numerator and denominator terms for all above-ground levels
1033 24162945056 : do k_avg = k_avg_lower, k_avg_upper
1034 :
1035 12240879056 : numer_integral = numer_integral + numer_terms(i,k_avg)
1036 24162945056 : denom_integral = denom_integral + denom_terms(i,k_avg)
1037 :
1038 : end do
1039 :
1040 12681545760 : Lscale_width_vert_avg_output(i,k) = numer_integral / denom_integral
1041 :
1042 : end do
1043 : end do
1044 :
1045 : !$acc exit data delete( one_half_avg_width, numer_terms, denom_terms )
1046 :
1047 8935056 : return
1048 :
1049 8935056 : end function Lscale_width_vert_avg
1050 :
1051 : !=============================================================================
1052 8935056 : subroutine wp2_term_splat_lhs( nz, ngrdcol, gr, C_wp2_splat, &
1053 8935056 : brunt_vaisala_freq_sqd_splat, &
1054 8935056 : lhs_splat_wp2 )
1055 :
1056 : ! Description
1057 : ! DESCRIBE TERM
1058 :
1059 : ! References:
1060 : !-----------------------------------------------------------------------
1061 :
1062 : use grid_class, only: &
1063 : grid, & ! Type
1064 : zm2zt2zm
1065 :
1066 : use constants_clubb, only: &
1067 : zero
1068 :
1069 : use clubb_precision, only: &
1070 : core_rknd ! Variable(s)
1071 :
1072 : implicit none
1073 :
1074 : ! --------------------- Input Variables ---------------------
1075 : integer, intent(in) :: &
1076 : nz, &
1077 : ngrdcol
1078 :
1079 : type (grid), target, intent(in) :: &
1080 : gr
1081 :
1082 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
1083 : brunt_vaisala_freq_sqd_splat ! Inverse time-scale tau at momentum levels [1/s^2]
1084 :
1085 : real( kind = core_rknd ), dimension(ngrdcol), intent(in) :: &
1086 : C_wp2_splat ! Model parameter C_wp2_splat [ -]
1087 :
1088 : ! --------------------- Output Variable ---------------------
1089 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
1090 : lhs_splat_wp2 ! LHS coefficient of wp2 splatting term [1/s]
1091 :
1092 : ! --------------------- Local Variables ---------------------
1093 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
1094 17870112 : brunt_vaisala_freq_splat_clipped, &
1095 17870112 : brunt_vaisala_freq_splat_smooth
1096 :
1097 : integer :: i, k
1098 :
1099 : !----------------------------- Begin Code -----------------------------
1100 :
1101 : !$acc enter data create( brunt_vaisala_freq_splat_clipped, brunt_vaisala_freq_splat_smooth )
1102 :
1103 : !$acc parallel loop gang vector collapse(2) default(present)
1104 768414816 : do k = 1, nz
1105 12690480816 : do i = 1, ngrdcol
1106 23844132000 : brunt_vaisala_freq_splat_clipped(i,k) &
1107 36525677760 : = sqrt( max( zero, brunt_vaisala_freq_sqd_splat(i,k) ) )
1108 : end do
1109 : end do
1110 : !$acc end parallel loop
1111 :
1112 : brunt_vaisala_freq_splat_smooth = zm2zt2zm( nz, ngrdcol, gr, &
1113 8935056 : brunt_vaisala_freq_splat_clipped )
1114 :
1115 : !$acc parallel loop gang vector collapse(2) default(present)
1116 768414816 : do k = 1, nz
1117 12690480816 : do i = 1, ngrdcol
1118 12681545760 : lhs_splat_wp2(i,k) = + C_wp2_splat(i) * brunt_vaisala_freq_splat_smooth(i,k)
1119 : end do
1120 : end do
1121 : !$acc end parallel loop
1122 :
1123 : !$acc exit data delete( brunt_vaisala_freq_splat_clipped, brunt_vaisala_freq_splat_smooth )
1124 :
1125 8935056 : return
1126 :
1127 : end subroutine wp2_term_splat_lhs
1128 :
1129 : !=============================================================================
1130 8935056 : subroutine wp3_term_splat_lhs( nz, ngrdcol, gr, C_wp2_splat, &
1131 8935056 : brunt_vaisala_freq_sqd_splat, &
1132 8935056 : lhs_splat_wp3 )
1133 :
1134 : ! Description
1135 : ! DESCRIBE TERM
1136 :
1137 : ! References:
1138 : !-----------------------------------------------------------------------
1139 :
1140 : use grid_class, only: &
1141 : grid, & ! Type
1142 : zm2zt
1143 :
1144 : use constants_clubb, only: &
1145 : zero, &
1146 : one_half, &
1147 : three
1148 :
1149 : use clubb_precision, only: &
1150 : core_rknd ! Variable(s)
1151 :
1152 : implicit none
1153 :
1154 : ! --------------------- Input Variables ---------------------
1155 : integer, intent(in) :: &
1156 : nz, &
1157 : ngrdcol
1158 :
1159 : type (grid), target, intent(in) :: &
1160 : gr
1161 :
1162 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
1163 : brunt_vaisala_freq_sqd_splat ! Inverse time-scale tau at momentum levels [1/s^2]
1164 :
1165 : real( kind = core_rknd ), dimension(ngrdcol), intent(in) :: &
1166 : C_wp2_splat ! Model parameter C_wp2_splat [-]
1167 :
1168 : ! --------------------- Output Variable ---------------------
1169 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
1170 : lhs_splat_wp3 ! LHS coefficient of wp3 splatting term [1/s]
1171 :
1172 : ! --------------------- Local Variables ---------------------
1173 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
1174 17870112 : brunt_vaisala_freq_splat_clipped, &
1175 17870112 : brunt_vaisala_freq_splat_smooth
1176 :
1177 : integer :: i, k
1178 :
1179 : !----------------------------- Begin Code -----------------------------
1180 :
1181 : !$acc enter data create( brunt_vaisala_freq_splat_clipped, brunt_vaisala_freq_splat_smooth )
1182 :
1183 : !$acc parallel loop gang vector collapse(2) default(present)
1184 768414816 : do k = 1, nz
1185 12690480816 : do i = 1, ngrdcol
1186 23844132000 : brunt_vaisala_freq_splat_clipped(i,k) &
1187 36525677760 : = sqrt( max( zero, brunt_vaisala_freq_sqd_splat(i,k) ) )
1188 : end do
1189 : end do
1190 : !$acc end parallel loop
1191 :
1192 : brunt_vaisala_freq_splat_smooth = zm2zt( nz, ngrdcol, gr, &
1193 8935056 : brunt_vaisala_freq_splat_clipped )
1194 :
1195 : !$acc parallel loop gang vector collapse(2) default(present)
1196 768414816 : do k = 1, nz
1197 12690480816 : do i = 1, ngrdcol
1198 23844132000 : lhs_splat_wp3(i,k) = + one_half * three * C_wp2_splat(i) &
1199 36525677760 : * brunt_vaisala_freq_splat_smooth(i,k)
1200 : end do
1201 : end do
1202 : !$acc end parallel loop
1203 :
1204 : !$acc exit data delete( brunt_vaisala_freq_splat_clipped, brunt_vaisala_freq_splat_smooth )
1205 :
1206 8935056 : return
1207 :
1208 : end subroutine wp3_term_splat_lhs
1209 :
1210 : !===============================================================================
1211 0 : function smooth_min_sclr_idx( nz, ngrdcol, input_var1, input_var2, smth_coef ) &
1212 0 : result( output_var )
1213 :
1214 : ! Description:
1215 : ! Computes a smoothed version of the min function, using one scalar and
1216 : ! one 1d array as inputs. For more details, see the interface in this file.
1217 :
1218 : ! References:
1219 : ! See clubb:ticket:894, updated version: 965
1220 : !----------------------------------------------------------------------
1221 :
1222 : use clubb_precision, only: &
1223 : core_rknd ! Constant(s)
1224 :
1225 : use constants_clubb, only: &
1226 : one_half
1227 :
1228 : implicit none
1229 :
1230 : integer, intent(in) :: &
1231 : nz, &
1232 : ngrdcol
1233 :
1234 : !----------------------------- Input Variables -----------------------------
1235 : real ( kind = core_rknd ), intent(in) :: &
1236 : input_var1, & ! Units vary
1237 : smth_coef ! "intensity" of the smoothing. Should be of a similar magnitude to
1238 : ! that of the data structures input_var1 and input_var2
1239 :
1240 : real ( kind = core_rknd ), dimension(ngrdcol, nz), intent(in) :: &
1241 : input_var2 ! Units vary
1242 :
1243 : !----------------------------- Output Variables -----------------------------
1244 : real( kind = core_rknd ), dimension(ngrdcol, nz) :: &
1245 : output_var ! Same unit as input_var1 and input_var2
1246 :
1247 : !----------------------------- Local Variables -----------------------------
1248 : integer :: i, k
1249 :
1250 : !----------------------------- Begin Code -----------------------------
1251 :
1252 : !$acc data copyin( input_var2 ) &
1253 : !$acc copyout( output_var )
1254 :
1255 : !$acc parallel loop gang vector collapse(2) default(present)
1256 0 : do k = 1, nz
1257 0 : do i = 1, ngrdcol
1258 0 : output_var(i,k) = one_half * ( (input_var1+input_var2(i,k)) - &
1259 0 : sqrt((input_var1-input_var2(i,k))**2 + smth_coef**2) )
1260 : end do
1261 : end do
1262 : !$acc end parallel loop
1263 :
1264 : !$acc end data
1265 :
1266 0 : return
1267 :
1268 0 : end function smooth_min_sclr_idx
1269 :
1270 : !===============================================================================
1271 0 : function smooth_min_array_scalar( nz, ngrdcol, input_var1, input_var2, smth_coef ) &
1272 0 : result( output_var )
1273 :
1274 : ! Description:
1275 : ! Computes a smoothed version of the min function, using one scalar and
1276 : ! one 1d array as inputs. For more details, see the interface in this file.
1277 :
1278 : ! References:
1279 : ! See clubb:ticket:894, updated version: 965
1280 : !----------------------------------------------------------------------
1281 :
1282 : use clubb_precision, only: &
1283 : core_rknd ! Constant(s)
1284 :
1285 : use constants_clubb, only: &
1286 : one_half
1287 :
1288 : implicit none
1289 :
1290 : !----------------------------- Input Variables -----------------------------
1291 : integer, intent(in) :: &
1292 : nz, &
1293 : ngrdcol
1294 :
1295 : real ( kind = core_rknd ), dimension(ngrdcol, nz), intent(in) :: &
1296 : input_var1 ! Units vary
1297 :
1298 : real ( kind = core_rknd ), intent(in) :: &
1299 : input_var2, & ! Units vary
1300 : smth_coef ! "intensity" of the smoothing. Should be of a similar magnitude to
1301 : ! that of the data structures input_var1 and input_var2
1302 :
1303 : !----------------------------- Output Variables -----------------------------
1304 : real( kind = core_rknd ), dimension(ngrdcol, nz) :: &
1305 : output_var ! Same unit as input_var1 and input_var2
1306 :
1307 : !----------------------------- Local Variables -----------------------------
1308 : integer :: i, k
1309 :
1310 : !----------------------------- Begin Code -----------------------------
1311 :
1312 : !$acc data copyin( input_var1 ) &
1313 : !$acc copyout( output_var )
1314 :
1315 : !$acc parallel loop gang vector collapse(2) default(present)
1316 0 : do k = 1, nz
1317 0 : do i = 1, ngrdcol
1318 0 : output_var(i,k) = one_half * ( (input_var1(i,k)+input_var2) - &
1319 0 : sqrt((input_var1(i,k)-input_var2)**2 + smth_coef**2) )
1320 : end do
1321 : end do
1322 : !$acc end parallel loop
1323 :
1324 : !$acc end data
1325 :
1326 0 : return
1327 :
1328 0 : end function smooth_min_array_scalar
1329 :
1330 : !===============================================================================
1331 0 : function smooth_min_arrays( nz, ngrdcol, input_var1, input_var2, smth_coef ) &
1332 0 : result( output_var )
1333 :
1334 : ! Description:
1335 : ! Computes a smoothed version of the min function, using two 1d arrays as inputs.
1336 : ! For more details, see the interface in this file.
1337 :
1338 : ! References:
1339 : ! See clubb:ticket:894, updated version: 965
1340 : !----------------------------------------------------------------------
1341 :
1342 : use clubb_precision, only: &
1343 : core_rknd ! Constant(s)
1344 :
1345 : use constants_clubb, only: &
1346 : one_half
1347 :
1348 : implicit none
1349 :
1350 : !----------------------------- Input Variables-----------------------------
1351 : integer, intent(in) :: &
1352 : nz, &
1353 : ngrdcol
1354 :
1355 : real ( kind = core_rknd ), dimension(ngrdcol, nz), intent(in) :: &
1356 : input_var1, & ! Units vary
1357 : input_var2 ! Units vary
1358 :
1359 : real ( kind = core_rknd ), intent(in) :: &
1360 : smth_coef ! "intensity" of the smoothing. Should be of a similar magnitude to
1361 : ! that of the data structures input_var1 and input_var2
1362 :
1363 : !----------------------------- Output Variables -----------------------------
1364 : real( kind = core_rknd ), dimension(ngrdcol, nz) :: &
1365 : output_var ! Same unit as input_var1 and input_var2
1366 :
1367 : !----------------------------- Local Variables -----------------------------
1368 : integer :: i, k
1369 :
1370 : !----------------------------- Begin Code -----------------------------
1371 :
1372 : !$acc data copyin( input_var1, input_var2 ) &
1373 : !$acc copyout( output_var )
1374 :
1375 : !$acc parallel loop gang vector collapse(2) default(present)
1376 0 : do k = 1, nz
1377 0 : do i = 1, ngrdcol
1378 0 : output_var(i,k) = one_half * ( (input_var1(i,k)+input_var2(i,k)) - &
1379 0 : sqrt((input_var1(i,k)-input_var2(i,k))**2 + smth_coef**2) )
1380 : end do
1381 : end do
1382 : !$acc end parallel loop
1383 :
1384 : !$acc end data
1385 :
1386 0 : return
1387 :
1388 0 : end function smooth_min_arrays
1389 :
1390 : !===============================================================================
1391 0 : function smooth_min_scalars( input_var1, input_var2, smth_coef ) &
1392 : result( output_var )
1393 :
1394 : ! Description:
1395 : ! Computes a smoothed version of the min function, using two scalars as inputs.
1396 : ! For more details, see the interface in this file.
1397 :
1398 : ! References:
1399 : ! See clubb:ticket: 965
1400 : !----------------------------------------------------------------------
1401 :
1402 : use clubb_precision, only: &
1403 : core_rknd ! Constant(s)
1404 :
1405 : use constants_clubb, only: &
1406 : one_half
1407 :
1408 : implicit none
1409 :
1410 : ! Input Variables
1411 : real ( kind = core_rknd ), intent(in) :: &
1412 : input_var1, & ! Units vary
1413 : input_var2, & ! Units vary
1414 : smth_coef ! "intensity" of the smoothing. Should be of a similar magnitude to
1415 : ! that of the data structures input_var1 and input_var2
1416 :
1417 : ! Output Variables
1418 : real( kind = core_rknd ) :: &
1419 : output_var ! Same unit as input_var1 and input_var2
1420 :
1421 : !----------------------------------------------------------------------
1422 :
1423 : output_var = one_half * ( (input_var1+input_var2) - &
1424 0 : sqrt((input_var1-input_var2)**2 + smth_coef**2) )
1425 :
1426 : return
1427 : end function smooth_min_scalars
1428 :
1429 : !===============================================================================
1430 0 : function smooth_max_sclr_idx( nz, ngrdcol, input_var1, input_var2, smth_coef ) &
1431 0 : result( output_var )
1432 :
1433 : ! Description:
1434 : ! Computes a smoothed version of the max function, using one scalar and
1435 : ! one 1d array as inputs. For more details, see the interface in this file.
1436 :
1437 : ! References:
1438 : ! See clubb:ticket:894, updated version: 965
1439 : !----------------------------------------------------------------------
1440 :
1441 : use clubb_precision, only: &
1442 : core_rknd ! Constant(s)
1443 :
1444 : use constants_clubb, only: &
1445 : one_half
1446 :
1447 : implicit none
1448 :
1449 : !----------------------------- Input Variables -----------------------------
1450 : integer, intent(in) :: &
1451 : nz, &
1452 : ngrdcol
1453 :
1454 : real ( kind = core_rknd ), intent(in) :: &
1455 : input_var1, & ! Units vary
1456 : smth_coef ! "intensity" of the smoothing. Should be of a similar magnitude to
1457 : ! that of the data structures input_var1 and input_var2
1458 :
1459 : real ( kind = core_rknd ), dimension(ngrdcol, nz), intent(in) :: &
1460 : input_var2 ! Units vary
1461 :
1462 : !----------------------------- Output Variables -----------------------------
1463 : real( kind = core_rknd ), dimension(ngrdcol, nz) :: &
1464 : output_var ! Same unit as input_var1 and input_var2
1465 :
1466 : !----------------------------- Local Variables -----------------------------
1467 : integer :: i, k
1468 :
1469 : !----------------------------- Begin Code -----------------------------
1470 :
1471 : !$acc data copyin( input_var2 ) &
1472 : !$acc copyout( output_var )
1473 :
1474 : !$acc parallel loop gang vector collapse(2) default(present)
1475 0 : do k = 1, nz
1476 0 : do i = 1, ngrdcol
1477 0 : output_var(i,k) = one_half * ( (input_var1+input_var2(i,k)) + &
1478 0 : sqrt((input_var1-input_var2(i,k))**2 + smth_coef**2) )
1479 : end do
1480 : end do
1481 : !$acc end parallel loop
1482 :
1483 : !$acc end data
1484 :
1485 0 : return
1486 :
1487 0 : end function smooth_max_sclr_idx
1488 :
1489 : !===============================================================================
1490 0 : function smooth_max_array_scalar( nz, ngrdcol, input_var1, input_var2, smth_coef ) &
1491 0 : result( output_var )
1492 :
1493 : ! Description:
1494 : ! Computes a smoothed version of the max function, using one scalar and
1495 : ! one 1d array as inputs. For more details, see the interface in this file.
1496 :
1497 : ! References:
1498 : ! See clubb:ticket:894, updated version: 965
1499 : !----------------------------------------------------------------------
1500 :
1501 : use clubb_precision, only: &
1502 : core_rknd ! Constant(s)
1503 :
1504 : use constants_clubb, only: &
1505 : one_half
1506 :
1507 : implicit none
1508 :
1509 : !----------------------------- Input Variables -----------------------------
1510 : integer, intent(in) :: &
1511 : nz, &
1512 : ngrdcol
1513 :
1514 : real ( kind = core_rknd ), dimension(ngrdcol, nz), intent(in) :: &
1515 : input_var1 ! Units vary
1516 :
1517 : real ( kind = core_rknd ), intent(in) :: &
1518 : input_var2, & ! Units vary
1519 : smth_coef ! "intensity" of the smoothing. Should be of a similar magnitude to
1520 : ! that of the data structures input_var1 and input_var2
1521 :
1522 : !----------------------------- Output Variables -----------------------------
1523 : real( kind = core_rknd ), dimension(ngrdcol, nz) :: &
1524 : output_var ! Same unit as input_var1 and input_var2
1525 :
1526 : !----------------------------- Local Variables -----------------------------
1527 : integer :: i, k
1528 :
1529 : !----------------------------- Begin Code -----------------------------
1530 :
1531 : !$acc data copyin( input_var1 ) &
1532 : !$acc copyout( output_var )
1533 :
1534 : !$acc parallel loop gang vector collapse(2) default(present)
1535 0 : do k = 1, nz
1536 0 : do i = 1, ngrdcol
1537 0 : output_var(i,k) = one_half * ( ( input_var1(i,k) + input_var2 ) + &
1538 0 : sqrt(( input_var1(i,k) - input_var2 )**2 + smth_coef**2) )
1539 : end do
1540 : end do
1541 : !$acc end parallel loop
1542 :
1543 : !$acc end data
1544 :
1545 0 : return
1546 :
1547 0 : end function smooth_max_array_scalar
1548 :
1549 : !===============================================================================
1550 0 : function smooth_max_array_1D_scalar( ngrdcol, input_var1, input_var2, smth_coef ) &
1551 0 : result( output_var )
1552 :
1553 : ! Description:
1554 : ! Computes a smoothed version of the max function, using one scalar and
1555 : ! one 1d array as inputs. For more details, see the interface in this file.
1556 :
1557 : ! References:
1558 : ! See clubb:ticket:894, updated version: 965
1559 : !----------------------------------------------------------------------
1560 :
1561 : use clubb_precision, only: &
1562 : core_rknd ! Constant(s)
1563 :
1564 : use constants_clubb, only: &
1565 : one_half
1566 :
1567 : implicit none
1568 :
1569 : !----------------------------- Input Variables -----------------------------
1570 : integer, intent(in) :: &
1571 : ngrdcol
1572 :
1573 : real ( kind = core_rknd ), dimension(ngrdcol), intent(in) :: &
1574 : input_var1 ! Units vary
1575 :
1576 : real ( kind = core_rknd ), intent(in) :: &
1577 : input_var2, & ! Units vary
1578 : smth_coef ! "intensity" of the smoothing. Should be of a similar magnitude to
1579 : ! that of the data structures input_var1 and input_var2
1580 :
1581 : !----------------------------- Output Variables -----------------------------
1582 : real( kind = core_rknd ), dimension(ngrdcol) :: &
1583 : output_var ! Same unit as input_var1 and input_var2
1584 :
1585 : !----------------------------- Local Variables -----------------------------
1586 : integer :: i
1587 :
1588 : !----------------------------- Begin Code -----------------------------
1589 :
1590 : !$acc data copyin( input_var1 ) &
1591 : !$acc copyout( output_var )
1592 :
1593 : !$acc parallel loop gang vector default(present)
1594 0 : do i = 1, ngrdcol
1595 0 : output_var(i) = one_half * ( ( input_var1(i) + input_var2 ) + &
1596 0 : sqrt(( input_var1(i) - input_var2 )**2 + smth_coef**2) )
1597 : end do
1598 : !$acc end parallel loop
1599 :
1600 : !$acc end data
1601 :
1602 0 : return
1603 :
1604 0 : end function smooth_max_array_1D_scalar
1605 :
1606 : !===============================================================================
1607 0 : function smooth_max_arrays( nz, ngrdcol, input_var1, input_var2, smth_coef ) &
1608 0 : result( output_var )
1609 :
1610 : ! Description:
1611 : ! Computes a smoothed version of the max function, using two 1d arrays as inputs.
1612 : ! For more details, see the interface in this file.
1613 :
1614 : ! References:
1615 : ! See clubb:ticket:894, updated version: 965
1616 : !----------------------------------------------------------------------
1617 :
1618 : use clubb_precision, only: &
1619 : core_rknd ! Constant(s)
1620 :
1621 : use constants_clubb, only: &
1622 : one_half
1623 :
1624 : implicit none
1625 :
1626 : !----------------------------- Input Variables -----------------------------
1627 : integer, intent(in) :: &
1628 : nz, &
1629 : ngrdcol
1630 :
1631 : real ( kind = core_rknd ), dimension(ngrdcol, nz), intent(in) :: &
1632 : input_var1, & ! Units vary
1633 : input_var2 ! Units vary
1634 :
1635 : real( kind = core_rknd ), intent(in) :: &
1636 : smth_coef ! "intensity" of the smoothing. Should be of a similar magnitude to
1637 : ! that of the data structures input_var1 and input_var2
1638 :
1639 : !----------------------------- Output Variables -----------------------------
1640 : real( kind = core_rknd ), dimension(ngrdcol, nz) :: &
1641 : output_var ! Same unit as input_var1 and input_var2
1642 :
1643 : !----------------------------- Local Variables -----------------------------
1644 : integer :: i, k
1645 :
1646 : !----------------------------- Begin Code -----------------------------
1647 :
1648 : !$acc data copyin( input_var1, input_var2 ) &
1649 : !$acc copyout( output_var )
1650 :
1651 : !$acc parallel loop gang vector collapse(2) default(present)
1652 0 : do k = 1, nz
1653 0 : do i = 1, ngrdcol
1654 0 : output_var(i,k) = one_half * ( (input_var1(i,k)+input_var2(i,k)) + &
1655 0 : sqrt((input_var1(i,k)-input_var2(i,k))**2 + smth_coef**2) )
1656 : end do
1657 : end do
1658 : !$acc end parallel loop
1659 :
1660 : !$acc end data
1661 :
1662 0 : return
1663 :
1664 0 : end function smooth_max_arrays
1665 :
1666 : !===============================================================================
1667 0 : function smooth_max_scalars( input_var1, input_var2, smth_coef ) &
1668 : result( output_var )
1669 :
1670 : ! Description:
1671 : ! Computes a smoothed version of the max function, using two scalars as inputs.
1672 : ! For more details, see the interface in this file.
1673 :
1674 : ! References:
1675 : ! See clubb:ticket: 965
1676 : !----------------------------------------------------------------------
1677 :
1678 : use clubb_precision, only: &
1679 : core_rknd ! Constant(s)
1680 :
1681 : use constants_clubb, only: &
1682 : one_half
1683 :
1684 : implicit none
1685 :
1686 : !----------------------------- Input Variables -----------------------------
1687 : real ( kind = core_rknd ), intent(in) :: &
1688 : input_var1, & ! Units vary
1689 : input_var2, & ! Units vary
1690 : smth_coef ! "intensity" of the smoothing. Should be of a similar magnitude to
1691 : ! that of the data structures input_var1 and input_var2
1692 :
1693 : !----------------------------- Output Variables -----------------------------
1694 : real( kind = core_rknd ) :: &
1695 : output_var ! Same unit as input_var1 and input_var2
1696 :
1697 : !----------------------------- Local Variables -----------------------------
1698 : integer :: i, k
1699 :
1700 : !----------------------------- Begin Code -----------------------------
1701 :
1702 : output_var = one_half * ( (input_var1+input_var2) + &
1703 0 : sqrt((input_var1-input_var2)**2 + smth_coef**2) )
1704 : return
1705 :
1706 : end function smooth_max_scalars
1707 :
1708 0 : function smooth_heaviside_peskin( nz, ngrdcol, input, smth_range ) &
1709 0 : result( smth_output )
1710 :
1711 : ! Description:
1712 : ! Computes a smoothed heaviside function as in
1713 : ! [Lin, Lee et al., 2005, A level set characteristic Galerkin
1714 : ! finite element method for free surface flows], equation (2)
1715 :
1716 : ! References:
1717 : ! See clubb:ticket:965
1718 : !----------------------------------------------------------------------
1719 :
1720 : use clubb_precision, only: &
1721 : core_rknd ! Constant(s)
1722 :
1723 : use constants_clubb, only: &
1724 : pi, invrs_pi, one, one_half, zero
1725 :
1726 : implicit none
1727 :
1728 : !------------------------- Input Variables -------------------------
1729 : integer, intent(in) :: &
1730 : nz, &
1731 : ngrdcol
1732 :
1733 : real ( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
1734 : input ! Units vary
1735 :
1736 : real ( kind = core_rknd ), intent(in) :: &
1737 : smth_range ! Smooth Heaviside function on [-smth_range, smth_range]
1738 :
1739 : !------------------------- Output Variables -------------------------
1740 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
1741 : smth_output ! Same units as input
1742 :
1743 : !------------------------- Local Variables -------------------------
1744 : real ( kind = core_rknd ) :: &
1745 : input_over_smth_range ! input divided by smth_range
1746 :
1747 : integer :: i, k
1748 :
1749 : !------------------------- Begin Code -------------------------
1750 :
1751 : !$acc data copyin( input ) &
1752 : !$acc copyout( smth_output )
1753 :
1754 : !$acc parallel loop gang vector collapse(2) default(present)
1755 0 : do k = 1, nz
1756 0 : do i = 1, ngrdcol
1757 :
1758 0 : if ( input(i,k) < -smth_range ) then
1759 0 : smth_output(i,k) = zero
1760 0 : elseif ( input(i,k) > smth_range ) then
1761 0 : smth_output(i,k) = one
1762 : else
1763 : ! Note that this case will only ever be reached if smth_range != 0,
1764 : ! so this division is fine and should not cause any issues
1765 0 : input_over_smth_range = input(i,k) / smth_range
1766 : smth_output(i,k) = one_half &
1767 : * (one + input_over_smth_range &
1768 0 : + invrs_pi * sin(pi * input_over_smth_range))
1769 : end if
1770 :
1771 : end do
1772 : end do
1773 : !$acc end parallel loop
1774 :
1775 : !$acc end data
1776 :
1777 0 : return
1778 :
1779 0 : end function smooth_heaviside_peskin
1780 :
1781 : !===============================================================================
1782 0 : subroutine calc_xpwp_1D( gr, Km_zm, xm, &
1783 0 : xpwp )
1784 :
1785 : ! Description:
1786 : ! Compute x'w' from x<k>, x<k+1>, Kh and invrs_dzm
1787 :
1788 : ! References:
1789 : ! None
1790 : !-----------------------------------------------------------------------
1791 :
1792 : use clubb_precision, only: &
1793 : core_rknd ! Variable(s)
1794 :
1795 : use grid_class, only: &
1796 : grid
1797 :
1798 : implicit none
1799 :
1800 : ! ----------------------- Input variables -----------------------
1801 : type (grid), target, intent(in) :: gr
1802 :
1803 : real( kind = core_rknd ), dimension(gr%nz), intent(in) :: &
1804 : Km_zm, & ! Eddy diff. (k momentum level) [m^2/s]
1805 : xm ! x (k thermo level) [units vary]
1806 :
1807 : ! ----------------------- Output variable -----------------------
1808 : real( kind = core_rknd ), dimension(gr%nz), intent(out) :: &
1809 : xpwp ! x'w' [(units vary)(m/s)]
1810 :
1811 : integer :: k
1812 :
1813 : ! ----------------------- Begin Code -----------------------
1814 :
1815 : ! Solve for x'w' at all intermediate model levels.
1816 0 : do k = 1, gr%nz-1
1817 0 : xpwp(k) = Km_zm(k) * gr%invrs_dzm(1,k) * ( xm(k+1) - xm(k) )
1818 : end do
1819 :
1820 0 : return
1821 : end subroutine calc_xpwp_1D
1822 :
1823 : !===============================================================================
1824 232311456 : subroutine calc_xpwp_2D( nz, ngrdcol, gr, &
1825 232311456 : Km_zm, xm, &
1826 232311456 : xpwp )
1827 :
1828 : ! Description:
1829 : ! Compute x'w' from x<k>, x<k+1>, Kh and invrs_dzm
1830 :
1831 : ! References:
1832 : ! None
1833 : !-----------------------------------------------------------------------
1834 :
1835 : use clubb_precision, only: &
1836 : core_rknd ! Variable(s)
1837 :
1838 : use grid_class, only: &
1839 : grid
1840 :
1841 : implicit none
1842 :
1843 : ! ----------------------- Input variables -----------------------
1844 : integer, intent(in) :: &
1845 : nz, &
1846 : ngrdcol
1847 :
1848 : type (grid), target, intent(in) :: gr
1849 :
1850 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
1851 : Km_zm, & ! Eddy diff. (k momentum level) [m^2/s]
1852 : xm ! x (k thermo level) [units vary]
1853 :
1854 : ! ----------------------- Output variable -----------------------
1855 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
1856 : xpwp ! x'w' [(units vary)(m/s)]
1857 :
1858 : integer :: i, k
1859 :
1860 : ! ----------------------- Begin Code -----------------------
1861 :
1862 : !$acc data copyin( gr, gr%invrs_dzm, Km_zm, xm ) &
1863 : !$acc copyout( xpwp )
1864 :
1865 : ! Solve for x'w' at all intermediate model levels.
1866 : !$acc parallel loop gang vector collapse(2) default(present)
1867 19746473760 : do k = 1, nz-1
1868 >32607*10^7 : do i = 1, ngrdcol
1869 >32584*10^7 : xpwp(i,k) = Km_zm(i,k) * gr%invrs_dzm(i,k) * ( xm(i,k+1) - xm(i,k) )
1870 : end do
1871 : end do
1872 : !$acc end parallel loop
1873 :
1874 : !$acc end data
1875 :
1876 232311456 : return
1877 :
1878 : end subroutine calc_xpwp_2D
1879 :
1880 : !=============================================================================
1881 0 : function vertical_avg( total_idx, rho_ds, field, dz )
1882 :
1883 : ! Description:
1884 : ! Computes the density-weighted vertical average of a field.
1885 : !
1886 : ! The average value of a function, f, over a set domain, [a,b], is
1887 : ! calculated by the equation:
1888 : !
1889 : ! f_avg = ( INT(a:b) f*g ) / ( INT(a:b) g );
1890 : !
1891 : ! as long as f is continous and g is nonnegative and integrable. Therefore,
1892 : ! the density-weighted (by dry, static, base-static density) vertical
1893 : ! average value of any model field, x, is calculated by the equation:
1894 : !
1895 : ! x_avg|_z = ( INT(z_bot:z_top) x rho_ds dz )
1896 : ! / ( INT(z_bot:z_top) rho_ds dz );
1897 : !
1898 : ! where z_bot is the bottom of the vertical domain, and z_top is the top of
1899 : ! the vertical domain.
1900 : !
1901 : ! This calculation is done slightly differently depending on whether x is a
1902 : ! thermodynamic-level or a momentum-level variable.
1903 : !
1904 : ! Thermodynamic-level computation:
1905 :
1906 : !
1907 : ! For numerical purposes, INT(z_bot:z_top) x rho_ds dz, which is the
1908 : ! numerator integral, is calculated as:
1909 : !
1910 : ! SUM(k_bot:k_top) x(k) rho_ds(k) delta_z(k);
1911 : !
1912 : ! where k is the index of the given thermodynamic level, x and rho_ds are
1913 : ! both thermodynamic-level variables, and delta_z(k) = zm(k) - zm(k-1). The
1914 : ! indices k_bot and k_top are the indices of the respective lower and upper
1915 : ! thermodynamic levels involved in the integration.
1916 : !
1917 : ! Likewise, INT(z_bot:z_top) rho_ds dz, which is the denominator integral,
1918 : ! is calculated as:
1919 : !
1920 : ! SUM(k_bot:k_top) rho_ds(k) delta_z(k).
1921 : !
1922 : ! The first (k=1) thermodynamic level is below ground (or below the
1923 : ! official lower boundary at the first momentum level), so it should not
1924 : ! count in a vertical average, whether that vertical average is used for
1925 : ! the hole-filling scheme or for statistical purposes. Begin no lower
1926 : ! than level k=2, which is the first thermodynamic level above ground (or
1927 : ! above the model lower boundary).
1928 : !
1929 : ! For cases where hole-filling over the entire (global) vertical domain
1930 : ! is desired, or where statistics over the entire (global) vertical
1931 : ! domain are desired, the lower (thermodynamic-level) index of k = 2 and
1932 : ! the upper (thermodynamic-level) index of k = gr%nz, means that the
1933 : ! overall vertical domain will be gr%zm(1,gr%nz) - gr%zm(1,1).
1934 : !
1935 : !
1936 : ! Momentum-level computation:
1937 : !
1938 : ! For numerical purposes, INT(z_bot:z_top) x rho_ds dz, which is the
1939 : ! numerator integral, is calculated as:
1940 : !
1941 : ! SUM(k_bot:k_top) x(k) rho_ds(k) delta_z(k);
1942 : !
1943 : ! where k is the index of the given momentum level, x and rho_ds are both
1944 : ! momentum-level variables, and delta_z(k) = zt(k+1) - zt(k). The indices
1945 : ! k_bot and k_top are the indices of the respective lower and upper momentum
1946 : ! levels involved in the integration.
1947 : !
1948 : ! Likewise, INT(z_bot:z_top) rho_ds dz, which is the denominator integral,
1949 : ! is calculated as:
1950 : !
1951 : ! SUM(k_bot:k_top) rho_ds(k) delta_z(k).
1952 : !
1953 : ! The first (k=1) momentum level is right at ground level (or right at
1954 : ! the official lower boundary). The momentum level variables that call
1955 : ! the hole-filling scheme have set values at the surface (or lower
1956 : ! boundary), and those set values should not be changed. Therefore, the
1957 : ! vertical average (for purposes of hole-filling) should not include the
1958 : ! surface level (or lower boundary level). For hole-filling purposes,
1959 : ! begin no lower than level k=2, which is the second momentum level above
1960 : ! ground (or above the model lower boundary). Likewise, the value at the
1961 : ! model upper boundary (k=gr%nz) is also set for momentum level
1962 : ! variables. That value should also not be changed.
1963 : !
1964 : ! However, this function is also used to keep track (for statistical
1965 : ! purposes) of the vertical average of certain variables. In that case,
1966 : ! the vertical average needs to be taken over the entire vertical domain
1967 : ! (level 1 to level gr%nz).
1968 : !
1969 : !
1970 : ! In both the thermodynamic-level computation and the momentum-level
1971 : ! computation, the numerator integral is divided by the denominator integral
1972 : ! in order to find the average value (over the vertical domain) of x.
1973 :
1974 : ! References:
1975 : ! None
1976 : !-----------------------------------------------------------------------
1977 :
1978 : use clubb_precision, only: &
1979 : core_rknd ! Variable(s)
1980 :
1981 : implicit none
1982 :
1983 : ! Input variables
1984 : integer, intent(in) :: &
1985 : total_idx ! The total numer of indices within the range of averaging
1986 :
1987 : real( kind = core_rknd ), dimension(total_idx), intent(in) :: &
1988 : rho_ds, & ! Dry, static density on either thermodynamic or momentum levels [kg/m^3]
1989 : field, & ! The field (e.g. wp2) to be vertically averaged [Units vary]
1990 : dz ! Reciprocal of thermodynamic or momentum level thickness [1/m]
1991 : ! depending on whether we're on zt or zm grid.
1992 : ! Note: The rho_ds and field points need to be arranged from
1993 : ! lowest to highest in altitude, with rho_ds(1) and
1994 : ! field(1) actually their respective values at level k = 1.
1995 :
1996 : ! Output variable
1997 : real( kind = core_rknd ) :: &
1998 : vertical_avg ! Vertical average of field [Units of field]
1999 :
2000 : ! Local variables
2001 : real( kind = core_rknd ) :: &
2002 : numer_integral, & ! Integral in the numerator (see description)
2003 : denom_integral ! Integral in the denominator (see description)
2004 :
2005 :
2006 : integer :: k
2007 :
2008 : !-----------------------------------------------------------------------
2009 :
2010 : ! Initialize variable
2011 0 : numer_integral = 0.0_core_rknd
2012 0 : denom_integral = 0.0_core_rknd
2013 :
2014 : ! Compute the numerator and denominator integral.
2015 : ! Multiply rho_ds at level k by the level thickness
2016 : ! at level k. Then, sum over all vertical levels.
2017 0 : do k=1, total_idx
2018 :
2019 0 : numer_integral = numer_integral + rho_ds(k) * dz(k) * field(k)
2020 0 : denom_integral = denom_integral + rho_ds(k) * dz(k)
2021 :
2022 : end do
2023 :
2024 : ! Find the vertical average of 'field'.
2025 0 : vertical_avg = numer_integral / denom_integral
2026 : !vertical_avg = sum( rho_ds(:) * dz(:) * field(:) ) / sum( rho_ds(:) * dz(:) )
2027 :
2028 : return
2029 : end function vertical_avg
2030 :
2031 : !=============================================================================
2032 0 : function vertical_integral( total_idx, rho_ds, &
2033 0 : field, dz )
2034 :
2035 : ! Description:
2036 : ! Computes the vertical integral. rho_ds, field, and dz must all be
2037 : ! of size total_idx and should all start at the same index.
2038 : !
2039 :
2040 : ! References:
2041 : ! None
2042 : !-----------------------------------------------------------------------
2043 :
2044 : use clubb_precision, only: &
2045 : core_rknd ! Variable(s)
2046 :
2047 : implicit none
2048 :
2049 : ! Input variables
2050 : integer, intent(in) :: &
2051 : total_idx ! The total numer of indices within the range of averaging
2052 :
2053 : real( kind = core_rknd ), dimension(total_idx), intent(in) :: &
2054 : rho_ds, & ! Dry, static density [kg/m^3]
2055 : field, & ! The field to be vertically averaged [Units vary]
2056 : dz ! Level thickness [1/m]
2057 : ! Note: The rho_ds and field points need to be arranged from
2058 : ! lowest to highest in altitude, with rho_ds(1) and
2059 : ! field(1) actually their respective values at level k = k_start.
2060 :
2061 : ! Local variables
2062 : real( kind = core_rknd ) :: &
2063 : vertical_integral ! Integral in the numerator (see description)
2064 :
2065 : !-----------------------------------------------------------------------
2066 :
2067 : ! Assertion checks: that k_start <= gr%nz - 1
2068 : ! that k_end >= 2
2069 : ! that k_start <= k_end
2070 :
2071 :
2072 : ! Initializing vertical_integral to avoid a compiler warning.
2073 0 : vertical_integral = 0.0_core_rknd
2074 :
2075 : ! Compute the integral.
2076 : ! Multiply the field at level k by rho_ds at level k and by
2077 : ! the level thickness at level k. Then, sum over all vertical levels.
2078 : ! Note: The values of the field and rho_ds are passed into this function
2079 : ! so that field(1) and rho_ds(1) are actually the field and rho_ds
2080 : ! at level k_start.
2081 0 : vertical_integral = sum( field * rho_ds * dz )
2082 :
2083 : !print *, vertical_integral
2084 :
2085 : return
2086 : end function vertical_integral
2087 :
2088 :
2089 : end module advance_helper_module
|