Line data Source code
1 : !-----------------------------------------------------------------------
2 : ! $Id$
3 : !===============================================================================
4 : module mono_flux_limiter
5 :
6 : implicit none
7 :
8 : private ! Default Scope
9 :
10 : public :: monotonic_turbulent_flux_limit, &
11 : calc_turb_adv_range
12 :
13 : private :: mfl_xm_lhs, &
14 : mfl_xm_rhs, &
15 : mfl_xm_solve, &
16 : mean_vert_vel_up_down
17 :
18 : ! Private named constants to avoid string comparisons
19 : ! NOTE: These values must match the values for xm_wpxp_thlm
20 : ! and xm_wpxp_rtm given in advance_xm_wpxp_module!
21 : integer, parameter, private :: &
22 : mono_flux_thlm = 1, & ! Named constant for thlm mono_flux calls
23 : mono_flux_rtm = 2, & ! Named constant for rtm mono_flux calls
24 : mono_flux_um = 4, & ! Named constant for um mono_flux calls
25 : mono_flux_vm = 5 ! Named constant for vm mono_flux calls
26 :
27 : integer, parameter :: &
28 : ndiags3 = 3
29 :
30 : contains
31 :
32 : !=============================================================================
33 35740224 : subroutine monotonic_turbulent_flux_limit( nz, ngrdcol, gr, solve_type, dt, xm_old, &
34 35740224 : xp2, wm_zt, xm_forcing, &
35 35740224 : rho_ds_zm, rho_ds_zt, &
36 35740224 : invrs_rho_ds_zm, invrs_rho_ds_zt, &
37 : xp2_threshold, xm_tol, l_implemented, &
38 35740224 : low_lev_effect, high_lev_effect, &
39 : tridiag_solve_method, &
40 : l_upwind_xm_ma, &
41 : l_mono_flux_lim_spikefix, &
42 : stats_metadata, &
43 35740224 : stats_zt, stats_zm, &
44 35740224 : xm, wpxp )
45 :
46 : ! Description:
47 : ! Limits the value of w'x' and corrects the value of xm when the xm turbulent
48 : ! advection term is not monotonic. A monotonic turbulent advection scheme
49 : ! will not create new extrema for variable x, based only on turbulent
50 : ! advection (not considering mean advection and xm forcings).
51 : !
52 : ! Montonic turbulent advection
53 : ! ----------------------------
54 : !
55 : ! A monotonic turbulent advection scheme does not allow new extrema for
56 : ! variable x to be created (by means of turbulent advection). In a
57 : ! monotonic turbulent advection scheme, when only the effects of turbulent
58 : ! advection are considered (neglecting forcings and mean advection), the
59 : ! value of variable x at a given point should not increase above the
60 : ! greatest value of variable x at nearby points, nor decrease below the
61 : ! smallest value of variable x at nearby points. Nearby points are points
62 : ! that are close enough to the given point so that the value of variable x
63 : ! at the given point is effected by the values of variable x at the nearby
64 : ! points by means of transfer by turbulent winds during a time step. Again,
65 : ! a monotonic scheme insures that advection only transfers around values of
66 : ! variable x and does not create new extrema for variable x. A monotonic
67 : ! turbulent advection scheme is useful because the turbulent advection term
68 : ! (w'x') may go numerically unstable, resulting in large instabilities in
69 : ! the mean field (xm). A monotonic turbulent advection scheme will limit
70 : ! the change in xm, and also in w'x'.
71 : !
72 : ! The following example illustrates the concept of monotonic turbulent
73 : ! advection. Three successive vertical grid levels are shown (k-1, k, and
74 : ! k+1). Three point values of theta-l are listed at every vertical grid
75 : ! level. All three vertical levels have a mean theta-l (thlm) of 288.0 K.
76 : ! A circulation is occuring (in the direction of the arrows) in the vertical
77 : ! (w wind component) and in the horizontal (u and/or v wind components),
78 : ! such that the mean value of vertical velocity (wmm) is 0, but there is a
79 : ! turbulent component such that w'^2 > 0.
80 : !
81 : ! level = k+1 || --- 287.0 K --- 288.0 K --- 289.0 K --- || thlm = 288.0
82 : ! || / \--------------------->| ||
83 : ! || | | || wmm = 0; wp2 > 0
84 : ! || |<---------------------\ / ||
85 : ! level = k || --- 288.0 K --- 288.0 K --- 288.0 K --- || thlm = 288.0
86 : ! || |<---------------------/ \ ||
87 : ! || | | || wmm = 0; wp2 > 0
88 : ! || \ /--------------------->| ||
89 : ! level = k-1 || --- 287.5 K --- 288.0 K --- 288.5 K --- || thlm = 288.0
90 : !
91 : ! Neglecting any contributions from thlm forcings (effects of radiation,
92 : ! microphysics, large-scale horizontal advection, etc.), the values of
93 : ! theta-l as shown will be altered by only turbulent advection. As a side
94 : ! note, the contribution of mean advection will be 0 since wmm = 0. The
95 : ! diagram shows that the value of theta-l at the point on the right at level
96 : ! k will increase. However, the values of theta-l at the other two points
97 : ! at level k will remain the same. Thus, the value of thlm at level k will
98 : ! become greater than 288.0 K. In the same manner, the values of thlm at
99 : ! the other two vertical levels (k-1 and k+1) will become smaller than
100 : ! 288.0 K. However, the monotonic turbulent advection scheme insures that
101 : ! any theta-l point value cannot become smaller than the smallest theta-l
102 : ! point value (287.0 K) or larger than the largest theta-l point value
103 : ! (289.0 K). Since all theta-l point values must fall between 287.0 K and
104 : ! 289.0 K, the level averages of theta-l (thlm) must fall between 287.0 K
105 : ! and 289.0 K. Thus, any values of the turbulent flux, w'th_l', that would
106 : ! cause thlm to rise above 289.0 K or fall below 287.0 K, not considering
107 : ! the effect of other terms on thlm (such as forcings), are faulty and need
108 : ! to be limited appropriately. The values of thlm also need to be corrected
109 : ! appropriately.
110 : !
111 : ! Formula for the limitation of w'x' and xm
112 : ! -----------------------------------------
113 : !
114 : ! The equation for change in the mean field, xm, over time is:
115 : !
116 : ! d(xm)/dt = -w*d(xm)/dz - (1/rho_ds) * d( rho_ds * w'x' )/dz + xm_forcing;
117 : !
118 : ! where w*d(xm)/dz is the mean advection component,
119 : ! (1/rho_ds) * d( rho_ds * w'x' )/dz is the turbulent advection component,
120 : ! and xm_forcing is the xm forcing component. The d(xm)/dt time tendency
121 : ! component is discretized as:
122 : !
123 : ! xm(k,<t+1>)/dt = xm(k,<t>)/dt - w*d(xm)/dz
124 : ! - (1/rho_ds) * d( rho_ds * w'x' )/dz + xm_forcing.
125 : !
126 : ! The value of xm after it has been advanced to timestep (t+1) must be in an
127 : ! appropriate range based on the values of xm at timestep (t), the amount of
128 : ! xm forcings applied over the ensuing time step, and the amount of mean
129 : ! advection applied over the ensuing time step. This is exactly the same
130 : ! thing as saying that the value of xm(k,<t+1>), with the contribution of
131 : ! turbulent advection included, must fall into a certain range based on the
132 : ! value of xm(k,<t+1>) without the contribution of the turbulent advection
133 : ! component over the last time step. The following inequality is used to
134 : ! limit the value of xm(k,<t+1>):
135 : !
136 : ! MIN{ xm(k-1,<t>) + dt*xm_forcing(k-1) - dt*wm_zt(k-1)*d(xm)/dz|_(k-1)
137 : ! - x_max_dev_low(k-1,<t>),
138 : ! xm(k,<t>) + dt*xm_forcing(k) - dt*wm_zt(k)*d(xm)/dz|_(k)
139 : ! - x_max_dev_low(k,<t>),
140 : ! xm(k+1,<t>) + dt*xm_forcing(k+1) - dt*wm_zt(k+1)*d(xm)/dz|_(k+1)
141 : ! - x_max_dev_low(k+1,<t>) }
142 : ! <= xm(k,<t+1>) <=
143 : ! MAX{ xm(k-1,<t>) + dt*xm_forcing(k-1) - dt*wm_zt(k-1)*d(xm)/dz|_(k-1)
144 : ! + x_max_dev_high(k-1,<t>),
145 : ! xm(k,<t>) + dt*xm_forcing(k) - dt*wm_zt(k)*d(xm)/dz|_(k)
146 : ! + x_max_dev_high(k,<t>),
147 : ! xm(k+1,<t>) + dt*xm_forcing(k+1) - dt*wm_zt(k+1)*d(xm)/dz|_(k+1)
148 : ! + x_max_dev_high(k+1,<t>) };
149 : !
150 : ! where x_max_dev_low is the absolute value of the deviation from the mean
151 : ! of the smallest point value of variable x at the given vertical level and
152 : ! timestep; and where x_max_dev_high is the deviation from the mean of the
153 : ! largest point value of variable x at the given vertical level and
154 : ! timestep. For example, at vertical level (k+1) and timestep (t):
155 : !
156 : ! x_max_dev_low(k+1,<t>) = | MIN( x(k+1,<t>) ) - xm(k+1,<t>) |;
157 : ! x_max_dev_high(k+1,<t>) = MAX( x(k+1,<t>) ) - xm(k+1,<t>).
158 : !
159 : ! The inequality shown above only takes into account values from the central
160 : ! level, one-level-below the central level, and one-level-above the central
161 : ! level. This is the minimal amount of vertical levels that can have their
162 : ! values taken into consideration. Any vertical level that can have it's
163 : ! properties advect to the given level during the course of a single time
164 : ! step can be taken into consideration. However, only three levels will be
165 : ! considered in this example for the sake of simplicity.
166 : !
167 : ! The inequality will be written in more simple terms:
168 : !
169 : ! xm_lower_lim_allowable(k) <= xm(k,<t+1>) <= xm_upper_lim_allowable(k).
170 : !
171 : ! The inequality can now be related to the turbulent flux, w'x'(k,<t+1>),
172 : ! through a substitution that is made for xm(k,<t+1>), such that:
173 : !
174 : ! xm(k,<t+1>) = xm(k,<t>) + dt*xm_forcing(k) - dt*wm_zt(k)*d(xm)/dz|_(k)
175 : ! - dt * (1/rho_ds) * d( rho_ds * w'x' )/dz|_(k).
176 : !
177 : ! The inequality becomes:
178 : !
179 : ! xm_lower_lim_allowable(k)
180 : ! <=
181 : ! xm(k,<t>) + dt*xm_forcing(k) - dt*wm_zt(k)*d(xm)/dz|_(k)
182 : ! - dt * (1/rho_ds) * d( rho_ds * w'x' )/dz|_(k)
183 : ! <=
184 : ! xm_upper_lim_allowable(k).
185 : !
186 : ! The inequality is rearranged, and the turbulent advection term,
187 : ! d(w'x')/dz, is discretized:
188 : !
189 : ! xm_lower_lim_allowable(k)
190 : ! - [ xm(k,<t>) + dt*xm_forcing(k) - dt*wm_zt(k)*d(xm)/dz|_(k) ]
191 : ! <=
192 : ! - dt * (1/rho_ds_zt(k))
193 : ! * invrs_dzt(k)
194 : ! * [ rho_ds_zm(k) * w'x'(k,<t+1>)
195 : ! - rho_ds_zm(k-1) * w'x'(k-1,<t+1>) ]
196 : ! <=
197 : ! xm_upper_lim_allowable(k)
198 : ! - [ xm(k,<t>) + dt*xm_forcing(k) - dt*wm_zt(k)*d(xm)/dz|_(k) ];
199 : !
200 : ! where invrs_dzt(k) = 1 / ( zm(k) - zm(k-1) ).
201 : !
202 : ! Multiplying the inequality by -rho_ds_zt(k)/(dz*invrs_dzt(k)):
203 : !
204 : ! rho_ds_zt(k)/(dz*invrs_dzt(k))
205 : ! * [ xm(k,<t>) + dt*xm_forcing(k) - dt*wm_zt(k)*d(xm)/dz|_(k)
206 : ! - xm_lower_lim_allowable(k) ]
207 : ! >=
208 : ! rho_ds_zm(k) * w'x'(k,<t+1>) - rho_ds_zm(k-1) * w'x'(k-1,<t+1>)
209 : ! >=
210 : ! rho_ds_zt(k)/(dz*invrs_dzt(k))
211 : ! * [ xm(k,<t>) + dt*xm_forcing(k) - dt*wm_zt(k)*d(xm)/dz|_(k)
212 : ! - xm_upper_lim_allowable(k) ].
213 : !
214 : ! Note: The inequality symbols have been flipped due to multiplication
215 : ! involving a (-) sign.
216 : !
217 : ! Adding rho_ds_zm(k-1) * w'x'(k-1,<t+1>) to the inequality:
218 : !
219 : ! rho_ds_zt(k)/(dz*invrs_dzt(k))
220 : ! * [ xm(k,<t>) + dt*xm_forcing(k) - dt*wm_zt(k)*d(xm)/dz|_(k)
221 : ! - xm_lower_lim_allowable(k) ]
222 : ! + rho_ds_zm(k-1) * w'x'(k-1,<t+1>)
223 : ! >= rho_ds_zm(k) * w'x'(k,<t+1>) >=
224 : ! rho_ds_zt(k)/(dz*invrs_dzt(k))
225 : ! * [ xm(k,<t>) + dt*xm_forcing(k) - dt*wm_zt(k)*d(xm)/dz|_(k)
226 : ! - xm_upper_lim_allowable(k) ]
227 : ! + rho_ds_zm(k-1) * w'x'(k-1,<t+1>).
228 : !
229 : ! The inequality is then rearranged to be based around w'x'(k,<t+1>):
230 : !
231 : ! (1/rho_ds_zm(k))
232 : ! * [ rho_ds_zt(k)/(dt*invrs_dzt(k))
233 : ! * { xm(k,<t>) + dt*xm_forcing(k) - dt*wm_zt(k)*d(xm)/dz|_(k)
234 : ! - xm_lower_lim_allowable(k) }
235 : ! + rho_ds_zm(k-1) * w'x'(k-1,<t+1>) ]
236 : ! >= w'x'(k,<t+1>) >=
237 : ! (1/rho_ds_zm(k))
238 : ! * [ rho_ds_zt(k)/(dt*invrs_dzt(k))
239 : ! * { xm(k,<t>) + dt*xm_forcing(k) - dt*wm_zt(k)*d(xm)/dz|_(k)
240 : ! - xm_upper_lim_allowable(k) }
241 : ! + rho_ds_zm(k-1) * w'x'(k-1,<t+1>) ].
242 : !
243 : ! The values of w'x' are found on the momentum levels, while the values of
244 : ! xm are found on the thermodynamic levels. Additionally, the values of
245 : ! rho_ds_zm are found on the momentum levels, and the values of rho_ds_zt
246 : ! are found on the thermodynamic levels. The inequality is applied to
247 : ! w'x'(k,<t+1>) from vertical levels 2 through the second-highest level
248 : ! (gr%nz-1). The value of w'x' at level 1 is a set surface (or lowest
249 : ! level) flux. The value of w'x' at the highest level is also a set value,
250 : ! and therefore is not altered.
251 : !
252 : ! Approximating maximum and minimum values of x at any given vertical level
253 : ! -------------------------------------------------------------------------
254 : !
255 : ! The CLUBB code provides means, variances, and covariances for certain
256 : ! variables at all vertical levels. However, there is no way to find the
257 : ! maximum or minimum point value of any variable on any vertical level.
258 : ! Without that information, x_max_dev_low and x_max_dev_high can't be found,
259 : ! and the inequality above is useless. However, there is a way to
260 : ! approximate the maximum and minimum point values at any given vertical
261 : ! level. The maximum and minimum point values can be approximated through
262 : ! the use of the variance, x'^2.
263 : !
264 : ! Just as the mean value of x, which is xm, and the turbulent flux of x,
265 : ! which is w'x', are known, so is the variance of x, which is x'^2. The
266 : ! standard deviation of x is the square root of the variance of x. The
267 : ! distribution of x along the horizontal plane (at vertical level k) is
268 : ! approximated to be the sum of two normal (or Gaussian) distributions.
269 : ! Most of the values in a normal distribution are found within 2 standard
270 : ! deviations from the mean. Thus, the maximum point value of x along the
271 : ! horizontal plance at any vertical level can be approximated as:
272 : ! xm + 2*sqrt(x'^2). Likewise, the minimum value of x along the horizontal
273 : ! plane at any vertical level can be approximated as: xm - 2*sqrt(x'^2).
274 : !
275 : ! The values of x'^2 are found on the momentum levels. The values of xm
276 : ! are found on the thermodynamic levels. Thus, the values of x'^2 are
277 : ! interpolated to the thermodynamic levels in order to find the maximum
278 : ! and minimum point values of variable x.
279 : !
280 : ! The one downfall of this method is that instabilities can arise in the
281 : ! model where unphysically large values of x'^2 are produced. Thus, this
282 : ! allows for an unphysically large deviation of xm from its values at the
283 : ! previous time step due to turbulent advection. Thus, for purposes of
284 : ! determining the maximum and minimum point values of x, a upper limit
285 : ! is placed on x'^2, in order to limit the standard deviation of x. This
286 : ! limit is only applied in this subroutine, and is not applied to x'^2
287 : ! elsewhere in the model code.
288 :
289 : ! References:
290 : !-----------------------------------------------------------------------
291 :
292 : use grid_class, only: &
293 : grid, & ! Type
294 : zm2zt ! Procedure(s)
295 :
296 : use constants_clubb, only: &
297 : zero_threshold, &
298 : eps, &
299 : fstderr
300 :
301 : use error_code, only: &
302 : clubb_at_least_debug_level, & ! Procedure
303 : err_code, & ! Error Indicator
304 : clubb_fatal_error ! Constant
305 :
306 : use clubb_precision, only: &
307 : core_rknd ! Variable(s)
308 :
309 : use advance_helper_module, only: &
310 : vertical_integral ! Procedure(s)
311 :
312 : use stats_type_utilities, only: &
313 : stat_begin_update, & ! Procedure(s)
314 : stat_end_update, &
315 : stat_update_var
316 :
317 : use stats_variables, only: &
318 : stats_metadata_type
319 :
320 : use stats_type, only: stats ! Type
321 :
322 : implicit none
323 :
324 : ! Constant Parameters
325 :
326 : ! Flag for using a semi-implicit, tridiagonal method to solve for xm(t+1)
327 : ! when xm(t+1) needs to be changed.
328 : logical, parameter :: l_mfl_xm_imp_adj = .true.
329 :
330 : !----------------------- Input Variables -----------------------
331 : integer, intent(in) :: &
332 : nz, &
333 : ngrdcol
334 :
335 : type (grid), target, intent(in) :: gr
336 :
337 : integer, intent(in) :: &
338 : solve_type ! Variables being solved for.
339 :
340 : real( kind = core_rknd ), intent(in) :: &
341 : dt ! Model timestep length [s]
342 :
343 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
344 : xm_old, & ! xm at previous time step (thermo. levs.) [units vary]
345 : xp2, & ! x'^2 (momentum levels) [units vary]
346 : wm_zt, & ! w wind component on thermodynamic levels [m/s]
347 : xm_forcing, & ! xm forcings (thermodynamic levels) [units vary]
348 : rho_ds_zm, & ! Dry, static density on momentum levels [kg/m^3]
349 : rho_ds_zt, & ! Dry, static density on thermo. levels [kg/m^3]
350 : invrs_rho_ds_zm, & ! Inv. dry, static density @ moment. levs. [m^3/kg]
351 : invrs_rho_ds_zt ! Inv. dry, static density @ thermo. levs. [m^3/kg]
352 :
353 : real( kind = core_rknd ), intent(in) :: &
354 : xp2_threshold, & ! Lower limit of x'^2 [units vary]
355 : xm_tol ! Lower limit of maxdev [units vary]
356 :
357 : logical, intent(in) :: &
358 : l_implemented ! Flag for CLUBB being implemented in a larger model.
359 :
360 : integer, dimension(ngrdcol,nz), intent(in) :: &
361 : low_lev_effect, & ! Index of lowest level that has an effect (for lev. k)
362 : high_lev_effect ! Index of highest level that has an effect (for lev. k)
363 :
364 : integer, intent(in) :: &
365 : tridiag_solve_method ! Specifier for method to solve tridiagonal systems
366 :
367 : logical, intent(in) :: &
368 : l_upwind_xm_ma, & ! This flag determines whether we want to use an upwind differencing
369 : ! approximation rather than a centered differencing for turbulent or
370 : ! mean advection terms. It affects rtm, thlm, sclrm, um and vm.
371 : l_mono_flux_lim_spikefix ! Flag to implement monotonic flux limiter code that
372 : ! eliminates spurious drying tendencies at model top
373 :
374 : type (stats_metadata_type), intent(in) :: &
375 : stats_metadata
376 :
377 : !----------------------- Input/Output Variables -----------------------
378 : type (stats), target, dimension(ngrdcol), intent(inout) :: &
379 : stats_zt, &
380 : stats_zm
381 :
382 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(inout) :: &
383 : xm, & ! xm at current time step (thermodynamic levels) [units vary]
384 : wpxp ! w'x' (momentum levels) [units vary]
385 :
386 : !----------------------- Local Variables -----------------------
387 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
388 71480448 : xp2_zt, & ! x'^2 interpolated to thermodynamic levels [units vary]
389 71480448 : xm_enter_mfl, & ! xm as it enters the MFL [units vary]
390 71480448 : xm_without_ta, & ! Value of xm without turb. adv. contrib. [units vary]
391 71480448 : wpxp_net_adjust ! Net amount of adjustment needed on w'x' [units vary]
392 :
393 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
394 71480448 : min_x_allowable_lev, & ! Smallest usuable value of x at lev k [units vary]
395 71480448 : max_x_allowable_lev, & ! Largest usuable value of x at lev k [units vary]
396 71480448 : min_x_allowable, & ! Smallest usuable x within k +/- num_levs [units vary]
397 71480448 : max_x_allowable, & ! Largest usuable x within k +/- num_levs [units vary]
398 71480448 : wpxp_mfl_max, & ! Upper limit on w'x'(k) [units vary]
399 71480448 : wpxp_mfl_min ! Lower limit on w'x'(k) [units vary]
400 :
401 : real( kind = core_rknd ) :: &
402 : max_xp2, & ! Maximum allowable x'^2 [units vary]
403 : max_dev, & ! Determines approximate upper/lower limit of x [units vary]
404 : m_adv_term, & ! Contribution of mean advection to d(xm)/dt [units vary]
405 : xm_density_weighted, & ! Density weighted xm at domain top [units vary]
406 : xm_adj_coef, & ! Coeffecient to eliminate spikes at domain top [units vary]
407 : xm_vert_integral, & ! Vertical integral of xm [units_vary]
408 : dxm_dt_mfl_adjust, & ! Rate of change of adjustment to xm [units vary]
409 : dz ! zm grid spacing at top of domain [m]
410 :
411 : real( kind = core_rknd ), dimension(ndiags3,ngrdcol,nz) :: &
412 71480448 : lhs_mfl_xm ! Left hand side of tridiagonal matrix
413 :
414 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
415 71480448 : rhs_mfl_xm ! Right hand side of tridiagonal matrix equation
416 :
417 : integer :: &
418 : k, km1, i, j ! Array indices
419 :
420 : ! integer, parameter :: &
421 : ! num_levs = 10 ! Number of levels above and below level k to look for
422 : ! ! maxima and minima of variable x.
423 :
424 : integer :: &
425 : low_lev, & ! Lowest level (from level k) to look for x minima and maxima
426 : high_lev ! Highest level (from level k) to look for x minima and maxima
427 :
428 : integer :: &
429 : iwpxp_mfl, &
430 : ixm_mfl
431 :
432 : logical, dimension(ngrdcol) :: &
433 71480448 : l_adjustment_needed ! Indicates if we need an adjustment for a column
434 :
435 : logical:: &
436 : l_any_adjustment_needed
437 :
438 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
439 71480448 : xm_mfl
440 :
441 : real( kind = core_rknd ) :: &
442 : max_tmp, &
443 : min_tmp
444 :
445 : real( kind = core_rknd ), dimension(nz) :: &
446 71480448 : tmp_in
447 :
448 : !---------------------------- Begin Code ----------------------------
449 :
450 : !$acc enter data create( xp2_zt, xm_enter_mfl, xm_without_ta, wpxp_net_adjust, &
451 : !$acc min_x_allowable_lev, max_x_allowable_lev, min_x_allowable, &
452 : !$acc max_x_allowable, wpxp_mfl_max, wpxp_mfl_min, lhs_mfl_xm, &
453 : !$acc rhs_mfl_xm, l_adjustment_needed, xm_mfl )
454 :
455 44675280 : select case( solve_type )
456 : case ( mono_flux_rtm ) ! rtm/wprtp
457 8935056 : iwpxp_mfl = stats_metadata%iwprtp_mfl
458 8935056 : ixm_mfl = stats_metadata%irtm_mfl
459 8935056 : max_xp2 = 5.0e-6_core_rknd
460 : case ( mono_flux_thlm ) ! thlm/wpthlp
461 8935056 : iwpxp_mfl = stats_metadata%iwpthlp_mfl
462 8935056 : ixm_mfl = stats_metadata%ithlm_mfl
463 8935056 : max_xp2 = 5.0_core_rknd
464 : case ( mono_flux_um ) ! um/upwp
465 8935056 : iwpxp_mfl = stats_metadata%iupwp_mfl
466 8935056 : ixm_mfl = stats_metadata%ium_mfl
467 8935056 : max_xp2 = 10.0_core_rknd
468 : case ( mono_flux_vm ) ! vm/vpwp
469 8935056 : iwpxp_mfl = stats_metadata%ivpwp_mfl
470 8935056 : ixm_mfl = stats_metadata%ivm_mfl
471 8935056 : max_xp2 = 10.0_core_rknd
472 : case default ! passive scalars are involved
473 0 : iwpxp_mfl = 0
474 0 : ixm_mfl = 0
475 35740224 : max_xp2 = 5.0_core_rknd
476 : end select
477 :
478 :
479 35740224 : if ( stats_metadata%l_stats_samp ) then
480 : !$acc update host( wpxp, xm )
481 0 : do i = 1, ngrdcol
482 0 : call stat_begin_update( nz, iwpxp_mfl, wpxp(i,:) / dt, & ! intent(in)
483 0 : stats_zm(i) ) ! intent(inout)
484 0 : tmp_in(1) = 0.0_core_rknd
485 0 : tmp_in(2:nz) = xm(i,2:nz)
486 : call stat_begin_update( nz, ixm_mfl, tmp_in / dt, & ! intent(in)
487 0 : stats_zt(i) ) ! intent(inout)
488 : end do
489 : endif
490 35740224 : if ( stats_metadata%l_stats_samp .and. solve_type == mono_flux_thlm ) then
491 : !$acc update host( xm, xm_old, wpxp )
492 0 : do i = 1, ngrdcol
493 0 : tmp_in(1) = 0.0_core_rknd
494 0 : tmp_in(2:nz) = xm(i,2:nz)
495 : call stat_update_var( stats_metadata%ithlm_enter_mfl, tmp_in, & ! intent(in)
496 0 : stats_zt(i) ) ! intent(inout)
497 0 : tmp_in(1) = 0.0_core_rknd
498 0 : tmp_in(2:nz) = xm_old(i,2:nz)
499 : call stat_update_var( stats_metadata%ithlm_old, tmp_in, & ! intent(in)
500 : stats_zt(i) ) ! intent(inout)
501 : call stat_update_var( stats_metadata%iwpthlp_enter_mfl, wpxp(i,:), & ! intent(in)
502 0 : stats_zm(i) ) ! intent(inout)
503 : end do
504 35740224 : elseif ( stats_metadata%l_stats_samp .and. solve_type == mono_flux_rtm ) then
505 : !$acc update host( xm, xm_old, wpxp )
506 0 : do i = 1, ngrdcol
507 0 : tmp_in(1) = 0.0_core_rknd
508 0 : tmp_in(2:nz) = xm(i,2:nz)
509 : call stat_update_var( stats_metadata%irtm_enter_mfl, tmp_in, & ! intent(in)
510 0 : stats_zt(i) ) ! intent(inout)
511 0 : tmp_in(1) = 0.0_core_rknd
512 0 : tmp_in(2:nz) = xm_old(i,2:nz)
513 : call stat_update_var( stats_metadata%irtm_old, tmp_in, & ! intent(in)
514 : stats_zt(i) ) ! intent(inout)
515 : call stat_update_var( stats_metadata%iwprtp_enter_mfl, wpxp(i,:), & ! intent(in)
516 0 : stats_zm(i) ) ! intent(inout)
517 : end do
518 : endif
519 :
520 :
521 : !$acc parallel loop gang vector collapse(2) default(present)
522 3073659264 : do k = 1, nz
523 50761923264 : do i = 1, ngrdcol
524 : ! Initialize arrays.
525 47688264000 : wpxp_net_adjust(i,k) = 0.0_core_rknd
526 :
527 : ! Store the value of xm as it enters the mfl
528 50726183040 : xm_enter_mfl(i,k) = xm(i,k)
529 : end do
530 : end do
531 : !$acc end parallel loop
532 :
533 : ! Interpolate x'^2 to thermodynamic levels.
534 35740224 : xp2_zt(:,:) = zm2zt( nz, ngrdcol, gr, xp2(:,:) )
535 :
536 : ! Place an upper limit on xp2_zt.
537 : ! For purposes of this subroutine, an upper limit has been placed on the
538 : ! variance, x'^2. This does not effect the value of x'^2 anywhere else in
539 : ! the model code. The upper limit is a reasonable upper limit. This is
540 : ! done to prevent unphysically large standard deviations caused by numerical
541 : ! instabilities in the x'^2 profile.
542 : !$acc parallel loop gang vector collapse(2) default(present)
543 3073659264 : do k = 1, nz
544 50761923264 : do i = 1, ngrdcol
545 50726183040 : xp2_zt(i,k) = min( max( xp2_zt(i,k), xp2_threshold ), max_xp2 )
546 : end do
547 : end do
548 : !$acc end parallel loop
549 :
550 : ! Find the maximum and minimum usuable values of variable x at each
551 : ! vertical level. Start from level 2, which is the first level above
552 : ! the ground (or above the model surface). This computation needs to be
553 : ! performed for all vertical levels above the ground (or model surface).
554 : !$acc parallel loop gang vector collapse(2) default(present)
555 3037919040 : do k = 2, nz, 1
556 50165144640 : do i = 1, ngrdcol
557 :
558 47127225600 : km1 = max( k-1, 1 )
559 : !kp1 = min( k+1, gr%nz )
560 :
561 : ! Most values are found within +/- 2 standard deviations from the mean.
562 : ! Use +/- 2 standard deviations from the mean as the maximum/minimum
563 : ! values.
564 : ! max_dev = 2.0_core_rknd*stnd_dev_x
565 :
566 : ! Set a minimum on max_dev
567 47127225600 : max_dev = max(2.0_core_rknd * sqrt( xp2_zt(i,k) ), xm_tol)
568 :
569 : ! Calculate the contribution of the mean advection term:
570 : ! m_adv_term = -wm_zt(k)*d(xm)/dz|_(k).
571 : ! Note: mean advection is not applied to xm at level gr%nz.
572 : !if ( .not. l_implemented .and. k < gr%nz ) then
573 : ! tmp(1:3) = term_ma_zt_lhs( wm_zt(k), gr%invrs_dzt(1,k), k )
574 : ! m_adv_term = - tmp(1) * xm(kp1) &
575 : ! - tmp(2) * xm(k) &
576 : ! - tmp(3) * xm(km1)
577 : !else
578 : ! m_adv_term = 0.0_core_rknd
579 : !endif
580 :
581 : ! Shut off to avoid using new, possibly corrupt mean advection term
582 47127225600 : m_adv_term = 0.0_core_rknd
583 :
584 : ! Find the value of xm without the contribution from the turbulent
585 : ! advection term.
586 : ! Note: the contribution of xm_forcing at level gr%nz should be 0.
587 : xm_without_ta(i,k) = xm_old(i,k) + dt*xm_forcing(i,k) &
588 47127225600 : + dt*m_adv_term
589 :
590 : ! Find the minimum usuable value of variable x at each vertical level.
591 47127225600 : if ( solve_type /= mono_flux_um .and. solve_type /= mono_flux_vm ) then
592 :
593 : ! Since variable x must be one of theta_l, r_t, or a scalar, all of
594 : ! which are positive definite quantities, the value must be >= 0.
595 : min_x_allowable_lev(i,k) &
596 23563612800 : = max( xm_without_ta(i,k) - max_dev, zero_threshold )
597 :
598 : else ! solve_type == mono_flux_um .or. solve_type == mono_flux_vm
599 :
600 : ! Variable x must be one of u or v.
601 23563612800 : min_x_allowable_lev(i,k) = xm_without_ta(i,k) - max_dev
602 :
603 : endif ! solve_type /= mono_flux_um .and. solve_type /= mono_flux_vm
604 :
605 : ! Find the maximum usuable value of variable x at each vertical level.
606 50129404416 : max_x_allowable_lev(i,k) = xm_without_ta(i,k) + max_dev
607 : end do
608 : end do
609 : !$acc end parallel loop
610 :
611 : ! Boundary condition on xm_without_ta
612 : !$acc parallel loop gang vector default(present)
613 596778624 : do i = 1, ngrdcol
614 561038400 : xm_without_ta(i,1) = xm(i,1)
615 561038400 : min_x_allowable_lev(i,1) = min_x_allowable_lev(i,2)
616 596778624 : max_x_allowable_lev(i,1) = max_x_allowable_lev(i,2)
617 : end do
618 : !$acc end parallel loop
619 :
620 : ! Find the maximum and minimum usuable values of x that can effect the value
621 : ! of x at level k. Then, find the upper and lower limits of w'x'. Reset
622 : ! the value of w'x' if it is outside of those limits, and store the amount
623 : ! of adjustment that was needed to w'x'.
624 : ! The values of w'x' at level 1 and at level gr%nz are set values and
625 : ! are not altered.
626 :
627 : ! Find the smallest and largest value of all relevant levels for variable x.
628 : !$acc parallel loop gang vector collapse(2) default(present)
629 3002178816 : do k = 2, nz-1
630 49568366016 : do i = 1, ngrdcol
631 :
632 46566187200 : low_lev = max( low_lev_effect(i,k), 2)
633 46566187200 : high_lev = min( high_lev_effect(i,k), nz )
634 :
635 46566187200 : min_tmp = min_x_allowable_lev(i,low_lev)
636 46566187200 : max_tmp = max_x_allowable_lev(i,low_lev)
637 :
638 57233821848 : do j = low_lev+1, high_lev
639 10667634648 : min_tmp = min( min_tmp, min_x_allowable_lev(i,j) )
640 57233821848 : max_tmp = max( max_tmp, max_x_allowable_lev(i,j) )
641 : end do
642 :
643 46566187200 : min_x_allowable(i,k) = min_tmp
644 49532625792 : max_x_allowable(i,k) = max_tmp
645 :
646 : end do
647 : end do
648 : !$acc end parallel loop
649 :
650 : !$acc parallel loop gang vector default(present)
651 596778624 : do i = 1, ngrdcol
652 47162965824 : do k = 2, nz-1
653 :
654 : ! Find the upper limit for w'x' for a monotonic turbulent flux.
655 : ! The following "if" statement ensures there are no "spikes" at the top of the column,
656 : ! which can cause unphysical rtm and thlm tendencies over the height of the column.
657 : ! The fix essentially turns off the monotonic flux limiter for these special cases,
658 : ! but tests show that it still performs well otherwise and runs stably.
659 : if ( l_mono_flux_lim_spikefix .and. solve_type == mono_flux_rtm &
660 93132374400 : .and. abs( wpxp(i,k-1) ) > 1 / ( dt * gr%invrs_dzt(i,k) ) &
661 46566187200 : * ( xm_without_ta(i,k) - min_x_allowable(i,k) ) &
662 >18626*10^7 : .and. wpxp(i,k-1) < 0.0_core_rknd ) then
663 532105 : wpxp_mfl_max(i,k) = 0.0_core_rknd
664 : else
665 : wpxp_mfl_max(i,k) &
666 : = invrs_rho_ds_zm(i,k) &
667 : * ( ( rho_ds_zt(i,k) / (dt*gr%invrs_dzt(i,k)) ) &
668 : * ( xm_without_ta(i,k) - min_x_allowable(i,k) ) &
669 46565655095 : + rho_ds_zm(i,k-1) * wpxp(i,k-1) )
670 : endif
671 :
672 : ! Find the lower limit for w'x' for a monotonic turbulent flux.
673 : wpxp_mfl_min(i,k) &
674 : = invrs_rho_ds_zm(i,k) &
675 0 : * ( ( rho_ds_zt(i,k) / (dt*gr%invrs_dzt(i,k)) ) &
676 : * ( xm_without_ta(i,k) - max_x_allowable(i,k) ) &
677 46566187200 : + rho_ds_zm(i,k-1) * wpxp(i,k-1) )
678 :
679 47127225600 : if ( wpxp(i,k) > wpxp_mfl_max(i,k) ) then
680 :
681 : ! This block of print statements can be uncommented for debugging.
682 : !print *, "k = ", k
683 : !print *, "wpxp too large (mfl)"
684 : !print *, "xm(t) = ", xm_old(k)
685 : !print *, "xm(t+1) entering mfl = ", xm(k)
686 : !print *, "xm(t+1) without ta = ", xm_without_ta(k)
687 : !print *, "max x allowable = ", max_x_allowable(k)
688 : !print *, "min x allowable = ", min_x_allowable(k)
689 : !print *, "1/rho_ds_zm(k) = ", invrs_rho_ds_zm(k)
690 : !print *, "rho_ds_zt(k) = ", rho_ds_zt(k)
691 : !print *, "rho_ds_zt(k)*(delta_zt/dt) = ", &
692 : ! real( rho_ds_zt(k) / (dt*gr%invrs_dzt(1,k)) )
693 : !print *, "xm without ta - min x allow = ", &
694 : ! xm_without_ta(k) - min_x_allowable(k)
695 : !print *, "rho_ds_zm(km1) = ", rho_ds_zm(km1)
696 : !print *, "wpxp(km1) = ", wpxp(km1)
697 : !print *, "rho_ds_zm(km1) * wpxp(km1) = ", rho_ds_zm(km1) * wpxp(km1)
698 : !print *, "wpxp upper lim = ", wpxp_mfl_max(k)
699 : !print *, "wpxp before adjustment = ", wpxp(k)
700 :
701 : ! Determine the net amount of adjustment needed for w'x'.
702 388202 : wpxp_net_adjust(i,k) = wpxp_mfl_max(i,k) - wpxp(i,k)
703 :
704 : ! Reset the value of w'x' to the upper limit allowed by the
705 : ! monotonic flux limiter.
706 388202 : wpxp(i,k) = wpxp_mfl_max(i,k)
707 :
708 46565798998 : elseif ( wpxp(i,k) < wpxp_mfl_min(i,k) ) then
709 :
710 : ! This block of print statements can be uncommented for debugging.
711 : !print *, "k = ", k
712 : !print *, "wpxp too small (mfl)"
713 : !print *, "xm(t) = ", xm_old(k)
714 : !print *, "xm(t+1) entering mfl = ", xm(k)
715 : !print *, "xm(t+1) without ta = ", xm_without_ta(k)
716 : !print *, "max x allowable = ", max_x_allowable(k)
717 : !print *, "min x allowable = ", min_x_allowable(k)
718 : !print *, "1/rho_ds_zm(k) = ", invrs_rho_ds_zm(k)
719 : !print *, "rho_ds_zt(k) = ", rho_ds_zt(k)
720 : !print *, "rho_ds_zt(k)*(delta_zt/dt) = ", &
721 : ! real( rho_ds_zt(k) / (dt*gr%invrs_dzt(1,k)) )
722 : !print *, "xm without ta - max x allow = ", &
723 : ! xm_without_ta(k) - max_x_allowable(k)
724 : !print *, "rho_ds_zm(km1) = ", rho_ds_zm(km1)
725 : !print *, "wpxp(km1) = ", wpxp(km1)
726 : !print *, "rho_ds_zm(km1) * wpxp(km1) = ", rho_ds_zm(km1) * wpxp(km1)
727 : !print *, "wpxp lower lim = ", wpxp_mfl_min(k)
728 : !print *, "wpxp before adjustment = ", wpxp(k)
729 :
730 : ! Determine the net amount of adjustment needed for w'x'.
731 840623 : wpxp_net_adjust(i,k) = wpxp_mfl_min(i,k) - wpxp(i,k)
732 :
733 : ! Reset the value of w'x' to the lower limit allowed by the
734 : ! monotonic flux limiter.
735 840623 : wpxp(i,k) = wpxp_mfl_min(i,k)
736 :
737 : ! This block of code can be uncommented for debugging.
738 : !else
739 : !
740 : ! ! wpxp(k) is okay.
741 : ! if ( wpxp_net_adjust(km1) /= 0.0_core_rknd ) then
742 : ! print *, "k = ", k
743 : ! print *, "wpxp is in an acceptable range (mfl)"
744 : ! print *, "xm(t) = ", xm_old(k)
745 : ! print *, "xm(t+1) entering mfl = ", xm(k)
746 : ! print *, "xm(t+1) without ta = ", xm_without_ta(k)
747 : ! print *, "max x allowable = ", max_x_allowable(k)
748 : ! print *, "min x allowable = ", min_x_allowable(k)
749 : ! print *, "1/rho_ds_zm(k) = ", invrs_rho_ds_zm(k)
750 : ! print *, "rho_ds_zt(k) = ", rho_ds_zt(k)
751 : ! print *, "rho_ds_zt(k)*(delta_zt/dt) = ", &
752 : ! real( rho_ds_zt(k) / (dt*gr%invrs_dzt(1,k)) )
753 : ! print *, "xm without ta - min x allow = ", &
754 : ! xm_without_ta(k) - min_x_allowable(k)
755 : ! print *, "xm without ta - max x allow = ", &
756 : ! xm_without_ta(k) - max_x_allowable(k)
757 : ! print *, "rho_ds_zm(km1) = ", rho_ds_zm(km1)
758 : ! print *, "wpxp(km1) = ", wpxp(km1)
759 : ! print *, "rho_ds_zm(km1) * wpxp(km1) = ", &
760 : ! rho_ds_zm(km1) * wpxp(km1)
761 : ! print *, "wpxp upper lim = ", wpxp_mfl_max(k)
762 : ! print *, "wpxp lower lim = ", wpxp_mfl_min(k)
763 : ! print *, "wpxp (stays the same) = ", wpxp(k)
764 : ! endif
765 : !
766 : endif
767 : end do
768 : end do
769 : !$acc end parallel loop
770 :
771 : ! Boundary conditions
772 : !$acc parallel loop gang vector default(present)
773 596778624 : do i = 1, ngrdcol
774 561038400 : min_x_allowable(i,1) = 0._core_rknd
775 561038400 : max_x_allowable(i,1) = 0._core_rknd
776 :
777 561038400 : min_x_allowable(i,nz) = 0._core_rknd
778 561038400 : max_x_allowable(i,nz) = 0._core_rknd
779 :
780 561038400 : wpxp_mfl_min(i,1) = 0._core_rknd
781 561038400 : wpxp_mfl_max(i,1) = 0._core_rknd
782 :
783 561038400 : wpxp_mfl_min(i,nz) = 0._core_rknd
784 596778624 : wpxp_mfl_max(i,nz) = 0._core_rknd
785 : end do
786 : !$acc end parallel loop
787 :
788 35740224 : if ( stats_metadata%l_stats_samp .and. solve_type == mono_flux_thlm ) then
789 : !$acc update host( xm_without_ta, min_x_allowable, wpxp_mfl_min, &
790 : !$acc wpxp_mfl_max, max_x_allowable )
791 0 : do i = 1, ngrdcol
792 0 : tmp_in(1) = 0.0_core_rknd
793 0 : tmp_in(2:nz) = xm_without_ta(i,2:nz)
794 : call stat_update_var( stats_metadata%ithlm_without_ta, tmp_in, & ! intent(in)
795 0 : stats_zt(i) ) ! intent(inout)
796 0 : tmp_in(1) = 0.0_core_rknd
797 0 : tmp_in(2:nz) = min_x_allowable(i,2:nz)
798 : call stat_update_var( stats_metadata%ithlm_mfl_min, tmp_in, & ! intent(in)
799 : stats_zt(i) ) ! intent(inout)
800 0 : tmp_in(1) = 0.0_core_rknd
801 0 : tmp_in(2:nz) = max_x_allowable(i,2:nz)
802 : call stat_update_var( stats_metadata%ithlm_mfl_max, tmp_in, & ! intent(in)
803 : stats_zt(i) ) ! intent(inout)
804 : call stat_update_var( stats_metadata%iwpthlp_mfl_min, wpxp_mfl_min(i,:), & ! intent(in)
805 0 : stats_zm(i) ) ! intent(inout)
806 : call stat_update_var( stats_metadata%iwpthlp_mfl_max, wpxp_mfl_max(i,:), & ! intent(in)
807 0 : stats_zm(i) ) ! intent(inout)
808 : end do
809 35740224 : elseif ( stats_metadata%l_stats_samp .and. solve_type == mono_flux_rtm ) then
810 : !$acc update host( xm_without_ta, min_x_allowable, max_x_allowable, &
811 : !$acc wpxp_mfl_min, wpxp_mfl_max )
812 0 : do i = 1, ngrdcol
813 0 : tmp_in(1) = 0.0_core_rknd
814 0 : tmp_in(2:nz) = xm_without_ta(i,2:nz)
815 : call stat_update_var( stats_metadata%irtm_without_ta, tmp_in, & ! intent(in)
816 0 : stats_zt(i) ) ! intent(inout)
817 0 : tmp_in(1) = 0.0_core_rknd
818 0 : tmp_in(2:nz) = min_x_allowable(i,2:nz)
819 : call stat_update_var( stats_metadata%irtm_mfl_min, tmp_in, & ! intent(in)
820 : stats_zt(i) ) ! intent(inout)
821 0 : tmp_in(1) = 0.0_core_rknd
822 0 : tmp_in(2:nz) = max_x_allowable(i,2:nz)
823 : call stat_update_var( stats_metadata%irtm_mfl_max, tmp_in, & ! intent(in)
824 : stats_zt(i) ) ! intent(inout)
825 : call stat_update_var( stats_metadata%iwprtp_mfl_min, wpxp_mfl_min(i,:), & ! intent(in)
826 0 : stats_zm(i) ) ! intent(inout)
827 : call stat_update_var( stats_metadata%iwprtp_mfl_max, wpxp_mfl_max(i,:), & ! intent(in)
828 0 : stats_zm(i) ) ! intent(inout)
829 : end do
830 : endif
831 :
832 596778624 : l_any_adjustment_needed = .false.
833 :
834 : !$acc parallel loop gang vector default(present)
835 596778624 : do i = 1, ngrdcol
836 596778624 : l_adjustment_needed(i) = .false.
837 : end do
838 : !$acc end parallel loop
839 :
840 : !$acc parallel loop gang vector collapse(2) default(present) &
841 : !$acc reduction(.or.:l_any_adjustment_needed)
842 596778624 : do i = 1, ngrdcol
843 48285042624 : do k = 1, nz
844 48249302400 : if ( abs(wpxp_net_adjust(i,k)) > eps ) then
845 1228823 : l_adjustment_needed(i) = .true.
846 1228823 : l_any_adjustment_needed = .true.
847 : end if
848 : end do
849 : end do
850 : !$acc end parallel loop
851 :
852 35740224 : if ( l_any_adjustment_needed ) then
853 :
854 : ! Reset the value of xm to compensate for the change to w'x'.
855 :
856 : if ( l_mfl_xm_imp_adj ) then
857 :
858 : ! A tridiagonal matrix is used to semi-implicitly re-solve for the
859 : ! values of xm at timestep index (t+1).
860 :
861 : ! Set up the left-hand side of the tridiagonal matrix equation.
862 : call mfl_xm_lhs( nz, ngrdcol, dt, gr%weights_zt2zm, & ! intent(in)
863 : gr%invrs_dzt, gr%invrs_dzm, & ! intent(in)
864 : wm_zt, l_implemented, l_upwind_xm_ma, & ! intent(in)
865 405213 : lhs_mfl_xm ) ! intent(out)
866 :
867 : ! Set up the right-hand side of tridiagonal matrix equation.
868 : call mfl_xm_rhs( nz, ngrdcol, dt, xm_old, wpxp, xm_forcing, & ! intent(in)
869 : gr%invrs_dzt, rho_ds_zm, invrs_rho_ds_zt, & ! intent(in)
870 405213 : rhs_mfl_xm ) ! intent(out)
871 :
872 : ! Solve the tridiagonal matrix equation.
873 : call mfl_xm_solve( nz, ngrdcol, solve_type, tridiag_solve_method, & ! intent(in)
874 : lhs_mfl_xm, rhs_mfl_xm, & ! intent(inout)
875 405213 : xm_mfl ) ! intent(out)
876 :
877 : ! If an adjustment is for a column
878 : !$acc parallel loop gang vector collapse(2) default(present)
879 34848318 : do k = 1, nz
880 576011358 : do i = 1, ngrdcol
881 575606145 : if ( l_adjustment_needed(i) ) then
882 87588505 : xm(i,k) = xm_mfl(i,k)
883 : end if
884 : end do
885 : end do
886 : !$acc end parallel loop
887 :
888 : ! Check for errors
889 405213 : if ( clubb_at_least_debug_level( 0 ) ) then
890 405213 : if ( err_code == clubb_fatal_error ) return
891 : end if
892 :
893 : else ! l_mfl_xm_imp_adj = .false.
894 :
895 : ! An explicit adjustment is made to the values of xm at timestep
896 : ! index (t+1), which is based upon the array of the amounts of w'x'
897 : ! adjustments.
898 :
899 : !$acc parallel loop gang vector collapse(2) default(present)
900 : do k = 2, nz, 1
901 : do i = 1, ngrdcol
902 :
903 : if ( l_adjustment_needed(i) ) then
904 :
905 : ! The rate of change of the adjustment to xm due to the monotonic
906 : ! flux limiter.
907 : dxm_dt_mfl_adjust = - invrs_rho_ds_zt(i,k) * gr%invrs_dzt(i,k) &
908 : * ( rho_ds_zm(i,k) * wpxp_net_adjust(i,k) &
909 : - rho_ds_zm(i,k-1) * wpxp_net_adjust(i,k-1) )
910 :
911 : ! The net change to xm due to the monotonic flux limiter is the
912 : ! rate of change multiplied by the time step length. Add the
913 : ! product to xm to find the new xm resulting from the monotonic
914 : ! flux limiter.
915 : xm(i,k) = xm(i,k) + dxm_dt_mfl_adjust * dt
916 : end if
917 :
918 : end do
919 : end do
920 : !$acc end parallel loop
921 :
922 : ! Boundary condition on xm
923 : !$acc parallel loop gang vector default(present)
924 : do i = 1, ngrdcol
925 : xm(i,1) = xm(i,2)
926 : end do
927 : !$acc end parallel loop
928 :
929 : endif ! l_mfl_xm_imp_adj
930 :
931 : ! This code can be uncommented for debugging.
932 : !do k = 1, gr%nz, 1
933 : ! print *, "k = ", k, "xm(t) = ", xm_old(k), "new xm(t+1) = ", xm(k)
934 : !enddo
935 :
936 : !Ensure there are no spikes at the top of the domain
937 : !$acc parallel loop gang vector default(present)
938 6771837 : do i = 1, ngrdcol
939 :
940 6771837 : if (abs( xm(i,nz) - xm_enter_mfl(i,nz) ) > 10._core_rknd * xm_tol) then
941 159 : dz = gr%zm(i,nz) - gr%zm(i,nz - 1)
942 :
943 : xm_density_weighted = rho_ds_zt(i,nz) &
944 : * (xm(i,nz) - xm_enter_mfl(i,nz)) &
945 159 : * dz
946 :
947 13356 : xm_vert_integral = sum(rho_ds_zt(i,2:nz-1) * xm(i,2:nz-1) * gr%dzt(i,2:nz-1) )
948 :
949 : !Check to ensure the vertical integral is not zero to avoid a divide
950 : !by zero error
951 159 : if ( abs(xm_vert_integral) < eps ) then
952 : #ifndef CLUBB_GPU
953 0 : write(fstderr,*) "Vertical integral of xm is zero;", &
954 0 : "mfl will remove spike at top of domain,", &
955 0 : "but it will not conserve xm."
956 : #endif
957 : !Remove the spike at the top of the domain
958 0 : xm(i,nz) = xm_enter_mfl(i,nz)
959 : else
960 159 : xm_adj_coef = xm_density_weighted / xm_vert_integral
961 :
962 : !xm_adj_coef can not be smaller than -1
963 159 : if (xm_adj_coef < -0.99_core_rknd) then
964 : #ifndef CLUBB_GPU
965 : write(fstderr,*) "xm_adj_coef in mfl less than -0.99, " &
966 0 : // "mx_adj_coef set to -0.99"
967 : #endif
968 0 : xm_adj_coef = -0.99_core_rknd
969 : endif
970 :
971 : !Apply the adjustment
972 13674 : xm(i,:) = xm(i,:) * (1._core_rknd + xm_adj_coef)
973 :
974 : !Remove the spike at the top of the domain
975 159 : xm(i,nz) = xm_enter_mfl(i,nz)
976 :
977 : !This code can be uncommented to ensure conservation
978 : !if (abs(sum(rho_ds_zt(2:gr%nz) * xm(2:gr%nz) / gr%invrs_dzt(2:gr%nz)) - &
979 : ! sum(rho_ds_zt(2:gr%nz) * xm_enter_mfl(2:gr%nz) / gr%invrs_dzt(2:gr%nz)))&
980 : ! > (1000 * xm_tol)) then
981 : ! write(fstderr,*) "NON-CONSERVATION in MFL", trim( solve_type ), &
982 : ! abs(sum(rho_ds_zt(2:gr%nz) * xm(2:gr%nz) / gr%invrs_dzt(2:gr%nz)) - &
983 : ! sum(rho_ds_zt(2:gr%nz) * xm_enter_mfl(2:gr%nz) / &
984 : ! gr%invrs_dzt(2:gr%nz)))
985 : !
986 : ! write(fstderr,*) "XM_ENTER_MFL=", xm_enter_mfl
987 : ! write(fstderr,*) "XM_AFTER_SPIKE_REMOVAL", xm
988 : ! write(fstderr,*) "XM_TOL", xm_tol
989 : ! write(fstderr,*) "XM_ADJ_COEF", xm_adj_coef
990 : !endif
991 :
992 : endif ! xm_vert_integral < eps
993 : endif ! spike at domain top
994 : end do
995 : !$acc end parallel loop
996 :
997 : end if
998 :
999 35740224 : if ( stats_metadata%l_stats_samp ) then
1000 : !$acc update host( wpxp, xm )
1001 0 : do i = 1, ngrdcol
1002 :
1003 0 : call stat_end_update( nz, iwpxp_mfl, wpxp(i,:) / dt, & ! intent(in)
1004 0 : stats_zm(i) ) ! intent(inout)
1005 :
1006 0 : tmp_in(1) = 0.0_core_rknd
1007 0 : tmp_in(2:nz) = xm(i,2:nz)
1008 : call stat_end_update( nz, ixm_mfl, tmp_in / dt, & ! intent(in)
1009 0 : stats_zt(i) ) ! intent(inout)
1010 :
1011 0 : if ( solve_type == mono_flux_thlm ) then
1012 0 : tmp_in(1) = 0.0_core_rknd
1013 0 : tmp_in(2:nz) = xm(i,2:nz)
1014 : call stat_update_var( stats_metadata%ithlm_exit_mfl, tmp_in, & ! intent(in)
1015 : stats_zt(i) ) ! intent(inout)
1016 : call stat_update_var( stats_metadata%iwpthlp_exit_mfl, wpxp(i,:), & ! intent(in)
1017 0 : stats_zm(i) ) ! intent(inout)
1018 0 : elseif ( solve_type == mono_flux_rtm ) then
1019 0 : tmp_in(1) = 0.0_core_rknd
1020 0 : tmp_in(2:nz) = xm(i,2:nz)
1021 : call stat_update_var( stats_metadata%irtm_exit_mfl, tmp_in, & ! intent(in)
1022 : stats_zt(i) ) ! intent(inout)
1023 : call stat_update_var( stats_metadata%iwprtp_exit_mfl, wpxp(i,:), & ! intent(in)
1024 0 : stats_zm(i) ) ! intent(inout)
1025 : endif
1026 : end do
1027 : endif
1028 :
1029 : !$acc exit data delete( xp2_zt, xm_enter_mfl, xm_without_ta, wpxp_net_adjust, &
1030 : !$acc min_x_allowable_lev, max_x_allowable_lev, min_x_allowable, &
1031 : !$acc max_x_allowable, wpxp_mfl_max, wpxp_mfl_min, lhs_mfl_xm, &
1032 : !$acc rhs_mfl_xm, l_adjustment_needed, xm_mfl )
1033 :
1034 : return
1035 :
1036 : end subroutine monotonic_turbulent_flux_limit
1037 :
1038 : !=============================================================================
1039 405213 : subroutine mfl_xm_lhs( nz, ngrdcol, dt, weights_zt2zm, &
1040 405213 : invrs_dzt, invrs_dzm, &
1041 405213 : wm_zt, l_implemented, l_upwind_xm_ma, &
1042 405213 : lhs )
1043 :
1044 : ! Description:
1045 : ! This subroutine is part of the process of re-solving for xm at timestep
1046 : ! index (t+1). This is done because the original solving process produced
1047 : ! values outside of what is deemed acceptable by the monotonic flux limiter.
1048 : ! Unlike the original formulation for advancing xm one timestep, which
1049 : ! combines w'x' and xm in a band-diagonal solver, this formulation uses a
1050 : ! tridiagonal solver to solve for only the value of xm(t+1), for w'x'(t+1)
1051 : ! is known.
1052 : !
1053 : ! Subroutine mfl_xm_lhs sets up the left-hand side of the matrix equation.
1054 :
1055 : use grid_class, only: &
1056 : grid ! Type
1057 :
1058 : use mean_adv, only: &
1059 : term_ma_zt_lhs ! Procedure(s)
1060 :
1061 : use clubb_precision, only: &
1062 : core_rknd ! Variable(s)
1063 :
1064 : implicit none
1065 :
1066 : ! Constant parameters
1067 : integer, parameter :: &
1068 : t_above = 1, & ! Index for upper thermodynamic level grid weight.
1069 : t_below = 2 ! Index for lower thermodynamic level grid weight.
1070 :
1071 : integer, parameter :: &
1072 : k_tdiag = 2 ! Thermodynamic main diagonal index.
1073 :
1074 : !---------------------------- Input Variables ----------------------------
1075 : integer, intent(in) :: &
1076 : nz, &
1077 : ngrdcol
1078 :
1079 : real( kind = core_rknd ), intent(in) :: &
1080 : dt ! Model timestep length [s]
1081 :
1082 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
1083 : wm_zt, & ! w wind component on thermodynamic levels [m/s]
1084 : invrs_dzt, &
1085 : invrs_dzm
1086 :
1087 : real( kind = core_rknd ), dimension(ngrdcol,nz,t_above:t_below), intent(in) :: &
1088 : weights_zt2zm
1089 :
1090 : logical, intent(in) :: &
1091 : l_implemented ! Flag for CLUBB being implemented in a larger model.
1092 :
1093 : logical, intent(in) :: &
1094 : l_upwind_xm_ma ! This flag determines whether we want to use an upwind differencing
1095 : ! approximation rather than a centered differencing for turbulent or
1096 : ! mean advection terms. It affects rtm, thlm, sclrm, um and vm.
1097 :
1098 : !---------------------------- Output Variables ----------------------------
1099 : real( kind = core_rknd ), dimension(ndiags3,ngrdcol,nz), intent(out) :: &
1100 : lhs ! Left hand side of tridiagonal matrix
1101 :
1102 : !---------------------------- Local Variables ----------------------------
1103 : integer :: i, k, b ! Array index
1104 :
1105 : !---------------------------- Begin Code ----------------------------
1106 :
1107 : ! The xm loop runs between k = 2 and k = nz. The value of xm at
1108 : ! level k = 1, which is below the model surface, is simply set equal to the
1109 : ! value of xm at level k = 2 after the solve has been completed.
1110 :
1111 : ! Setup LHS of the tridiagonal system
1112 :
1113 : ! LHS xm mean advection (ma) term.
1114 405213 : if ( .not. l_implemented ) then
1115 :
1116 : call term_ma_zt_lhs( nz, ngrdcol, wm_zt, weights_zt2zm, & ! intent(in)
1117 : invrs_dzt, invrs_dzm, & ! intent(in)
1118 : l_upwind_xm_ma, & ! intent(in)
1119 0 : lhs ) ! intent(out)
1120 : else
1121 : !$acc parallel loop gang vector collapse(3) default(present)
1122 34848318 : do k = 1, nz
1123 576011358 : do i = 1, ngrdcol
1124 2199095265 : do b = 1, ndiags3
1125 2164652160 : lhs(b,i,k) = 0.0_core_rknd
1126 : end do
1127 : end do
1128 : end do
1129 : !$acc end parallel loop
1130 : endif
1131 :
1132 : !$acc parallel loop gang vector collapse(2) default(present)
1133 34443105 : do k = 2, nz, 1
1134 569239521 : do i = 1, ngrdcol
1135 : ! LHS xm time tendency.
1136 568834308 : lhs(k_tdiag,i,k) = lhs(k_tdiag,i,k) + 1.0_core_rknd / dt
1137 : end do
1138 : end do ! xm loop: 2..nz
1139 : !$acc end parallel loop
1140 :
1141 : ! Boundary conditions.
1142 :
1143 : ! Lower boundary
1144 : !$acc parallel loop gang vector collapse(2) default(present)
1145 34848318 : do k = 1, nz
1146 576011358 : do i = 1, ngrdcol
1147 2164652160 : lhs(:,i,1) = 0.0_core_rknd
1148 575606145 : lhs(k_tdiag,i,1) = 1.0_core_rknd
1149 : end do
1150 : end do
1151 : !$acc end parallel loop
1152 :
1153 405213 : return
1154 :
1155 : end subroutine mfl_xm_lhs
1156 :
1157 : !=============================================================================
1158 405213 : subroutine mfl_xm_rhs( nz, ngrdcol, dt, xm_old, wpxp, xm_forcing, &
1159 405213 : invrs_dzt, rho_ds_zm, invrs_rho_ds_zt, &
1160 405213 : rhs )
1161 :
1162 : ! Description:
1163 : ! This subroutine is part of the process of re-solving for xm at timestep
1164 : ! index (t+1). This is done because the original solving process produced
1165 : ! values outside of what is deemed acceptable by the monotonic flux limiter.
1166 : ! Unlike the original formulation for advancing xm one timestep, which
1167 : ! combines w'x' and xm in a band-diagonal solver, this formulation uses a
1168 : ! tridiagonal solver to solve for only the value of xm(t+1), for w'x'(t+1)
1169 : ! is known.
1170 : !
1171 : ! Subroutine mfl_xm_rhs sets up the right-hand side of the matrix equation.
1172 :
1173 : use clubb_precision, only: &
1174 : core_rknd ! Variable(s)
1175 :
1176 : implicit none
1177 :
1178 : !---------------------------- Input Variables ----------------------------
1179 : integer, intent(in) :: &
1180 : nz, &
1181 : ngrdcol
1182 :
1183 : real( kind = core_rknd ), intent(in) :: &
1184 : dt ! Model timestep length [s]
1185 :
1186 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
1187 : invrs_dzt, & ! The inverse spacing between momentum grid levels;
1188 : ! centered over thermodynamic grid levels.
1189 : xm_old, & ! xm; timestep (t) (thermodynamic levels) [units vary]
1190 : wpxp, & ! w'x'; timestep (t+1); limited (m-levs.) [units vary]
1191 : xm_forcing, & ! xm forcings (thermodynamic levels) [units vary]
1192 : rho_ds_zm, & ! Dry, static density on momentum levels [kg/m^3]
1193 : invrs_rho_ds_zt ! Inv. dry, static density @ thermo. levs. [m^3/kg]
1194 :
1195 : !---------------------------- Output Variable ----------------------------
1196 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
1197 : rhs ! Right hand side of tridiagonal matrix equation
1198 :
1199 : !---------------------------- Local Variables ----------------------------
1200 : integer :: i, k ! Array indices
1201 :
1202 : !---------------------------- Begin Code ----------------------------
1203 :
1204 : ! The xm loop runs between k = 2 and k = gr%nz. The value of xm at
1205 : ! level k = 1, which is below the model surface, is simply set equal to the
1206 : ! value of xm at level k = 2 after the solve has been completed.
1207 :
1208 : !$acc parallel loop gang vector collapse(2) default(present)
1209 34443105 : do k = 2, nz, 1
1210 569239521 : do i = 1, ngrdcol
1211 :
1212 : ! RHS xm time tendency.
1213 534796416 : rhs(i,k) = xm_old(i,k) / dt
1214 :
1215 : ! RHS xm turbulent advection (ta) term.
1216 : ! Note: Normally, the turbulent advection (ta) term is treated
1217 : ! implicitly when advancing xm one timestep, as both xm and w'x'
1218 : ! are advanced together from timestep index (t) to timestep
1219 : ! index (t+1). However, in this case, both xm and w'x' have
1220 : ! already been advanced one timestep. However, w'x'(t+1) has been
1221 : ! limited after the fact, and therefore it's values at timestep
1222 : ! index (t+1) are known. Thus, in re-solving for xm(t+1), the
1223 : ! derivative of w'x'(t+1) can be placed on the right-hand side of
1224 : ! the d(xm)/dt equation.
1225 : rhs(i,k) = rhs(i,k) &
1226 : - invrs_rho_ds_zt(i,k) * invrs_dzt(i,k) &
1227 534796416 : * ( rho_ds_zm(i,k) * wpxp(i,k) - rho_ds_zm(i,k-1) * wpxp(i,k-1) )
1228 :
1229 : ! RHS xm forcings.
1230 : ! Note: xm forcings include the effects of microphysics,
1231 : ! cloud water sedimentation, radiation, and any
1232 : ! imposed forcings on xm.
1233 568834308 : rhs(i,k) = rhs(i,k) + xm_forcing(i,k)
1234 :
1235 : end do
1236 : end do ! xm loop: 2..gr%nz
1237 : !$acc end parallel loop
1238 :
1239 : ! Boundary conditions
1240 :
1241 : ! Lower Boundary
1242 : ! The value of xm at the lower boundary will remain the same. However, the
1243 : ! value of xm at the lower boundary gets overwritten after the matrix is
1244 : ! solved for the next timestep, such that xm(1) = xm(2).
1245 : !$acc parallel loop gang vector default(present)
1246 6771837 : do i = 1, ngrdcol
1247 6771837 : rhs(i,1) = xm_old(i,1)
1248 : end do
1249 : !$acc end parallel loop
1250 :
1251 405213 : return
1252 :
1253 : end subroutine mfl_xm_rhs
1254 :
1255 : !=============================================================================
1256 405213 : subroutine mfl_xm_solve( nz, ngrdcol, solve_type, tridiag_solve_method, &
1257 405213 : lhs, rhs, &
1258 405213 : xm )
1259 :
1260 : ! Description:
1261 : ! This subroutine is part of the process of re-solving for xm at timestep
1262 : ! index (t+1). This is done because the original solving process produced
1263 : ! values outside of what is deemed acceptable by the monotonic flux limiter.
1264 : ! Unlike the original formulation for advancing xm one timestep, which
1265 : ! combines w'x' and xm in a band-diagonal solver, this formulation uses a
1266 : ! tridiagonal solver to solve for only the value of xm(t+1), for w'x'(t+1)
1267 : ! is known.
1268 : !
1269 : ! Subroutine mfl_xm_solve solves the tridiagonal matrix equation for xm at
1270 : ! timestep index (t+1).
1271 :
1272 : use matrix_solver_wrapper, only: &
1273 : tridiag_solve ! Procedure(s)
1274 :
1275 : use clubb_precision, only: &
1276 : core_rknd
1277 :
1278 : use error_code, only: &
1279 : clubb_at_least_debug_level, & ! Procedure
1280 : err_code, & ! Error Indicator
1281 : clubb_fatal_error ! Constant
1282 :
1283 : implicit none
1284 :
1285 : ! Constant parameters
1286 : integer, parameter :: &
1287 : kp1_tdiag = 1, & ! Thermodynamic superdiagonal index.
1288 : k_tdiag = 2, & ! Thermodynamic main diagonal index.
1289 : km1_tdiag = 3 ! Thermodynamic subdiagonal index.
1290 :
1291 : !---------------------------- Input Variables ----------------------------
1292 : integer, intent(in) :: &
1293 : nz, &
1294 : ngrdcol
1295 :
1296 : integer, intent(in) :: &
1297 : solve_type ! Variables being solved for.
1298 :
1299 : integer, intent(in) :: &
1300 : tridiag_solve_method ! Specifier for method to solve tridiagonal systems
1301 :
1302 : !---------------------------- InOut Variables ----------------------------
1303 : real( kind = core_rknd ), dimension(ndiags3,ngrdcol,nz), intent(inout) :: &
1304 : lhs ! Left hand side of tridiagonal matrix
1305 :
1306 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(inout) :: &
1307 : rhs ! Right hand side of tridiagonal matrix equation
1308 :
1309 : !---------------------------- Output Variables ----------------------------
1310 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
1311 : xm ! Value of variable being solved for at timestep (t+1) [units vary]
1312 :
1313 : !---------------------------- Local Variable ----------------------------
1314 : character(len=10) :: &
1315 : solve_type_str ! solve_type as a string for debug output purposes
1316 :
1317 : integer :: i
1318 :
1319 : !---------------------------- Begin Code ----------------------------
1320 :
1321 430063 : select case( solve_type )
1322 : case ( mono_flux_rtm )
1323 24850 : solve_type_str = "rtm"
1324 : case ( mono_flux_thlm )
1325 203659 : solve_type_str = "thlm"
1326 : case default
1327 405213 : solve_type_str = "scalars"
1328 : end select
1329 :
1330 : ! Solve for xm at timestep index (t+1) using the tridiagonal solver.
1331 : call tridiag_solve( solve_type_str, tridiag_solve_method, & ! Intent(in)
1332 : ngrdcol, nz, & ! Intent(in)
1333 : lhs, rhs, & ! Intent(inout)
1334 405213 : xm ) ! Intent(out)
1335 :
1336 : ! Check for errors
1337 405213 : if ( clubb_at_least_debug_level( 0 ) ) then
1338 405213 : if ( err_code == clubb_fatal_error ) then
1339 : return
1340 : end if
1341 : end if
1342 :
1343 : ! Boundary condition on xm
1344 : !$acc parallel loop gang vector default(present)
1345 6771837 : do i = 1, ngrdcol
1346 6771837 : xm(i,1) = xm(i,2)
1347 : end do
1348 : !$acc end parallel loop
1349 :
1350 : return
1351 : end subroutine mfl_xm_solve
1352 :
1353 : !=============================================================================
1354 8935056 : subroutine calc_turb_adv_range( nz, ngrdcol, gr, dt, &
1355 8935056 : w_1_zm, w_2_zm, varnce_w_1_zm, varnce_w_2_zm, &
1356 8935056 : mixt_frac_zm, &
1357 : stats_metadata, &
1358 8935056 : stats_zm, &
1359 8935056 : low_lev_effect, high_lev_effect )
1360 :
1361 : ! Description:
1362 : ! Calculates the lowermost and uppermost thermodynamic grid levels that can
1363 : ! effect the base (or central) thermodynamic level through the effects of
1364 : ! turbulent advection over the course of one time step. This is used as
1365 : ! part of the monotonic turbulent advection scheme.
1366 : !
1367 : ! One method is to use the vertical velocity at each level to determine the
1368 : ! amount of time that it takes to travel across that particular grid level.
1369 : ! The method is to keep on advancing one grid level until either (a) the
1370 : ! total sum of time taken reaches or exceeds the model time step length,
1371 : ! (b) the top or bottom of the model is reached, or (c) a level is reached
1372 : ! where the vertical velocity component (with turbulence included) is
1373 : ! oriented completely opposite of the direction of travel towards the base
1374 : ! (or central) thermodynamic level. An example of situation (c) would be,
1375 : ! while starting from a higher altitude and searching downward for all
1376 : ! upward vertical velocity components, encountering a strong downdraft
1377 : ! where the vertical velocity at every single point is oriented downward.
1378 : ! Such a situation would occur when the mean vertical velocity (wm_zm)
1379 : ! exceeds any turbulent component (w') that would be oriented upwards.
1380 : !
1381 : ! Another method is to simply set the thickness (in meters) of the layer
1382 : ! that turbulent advection is allowed to act over, for purposes of the
1383 : ! monotonic turbulent advection scheme. The lowermost and uppermost
1384 : ! grid level that can effect the base (or central) thermodynamic level
1385 : ! is computed based on the thickness and altitude of each level.
1386 :
1387 : ! References:
1388 : !-----------------------------------------------------------------------
1389 :
1390 : use grid_class, only: &
1391 : grid ! Type
1392 :
1393 : use clubb_precision, only: &
1394 : core_rknd ! Variable(s)
1395 :
1396 : use stats_type, only: &
1397 : stats ! Type
1398 :
1399 : use stats_variables, only: &
1400 : stats_metadata_type
1401 :
1402 : implicit none
1403 :
1404 : ! Constant parameters
1405 : logical, parameter :: &
1406 : l_constant_thickness = .false. ! Toggle constant or variable thickness.
1407 :
1408 : real( kind = core_rknd ), parameter :: &
1409 : const_thick = 150.0_core_rknd ! Constant thickness value [m]
1410 :
1411 : !------------------------- Input Variables -------------------------
1412 : integer, intent(in) :: &
1413 : nz, &
1414 : ngrdcol
1415 :
1416 : type (grid), target, intent(in) :: gr
1417 :
1418 : real( kind = core_rknd ), intent(in) :: &
1419 : dt ! Model timestep length [s]
1420 :
1421 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
1422 : w_1_zm, & ! Mean w (1st PDF component) [m/s]
1423 : w_2_zm, & ! Mean w (2nd PDF component) [m/s]
1424 : varnce_w_1_zm, & ! Variance of w (1st PDF component) [m^2/s^2]
1425 : varnce_w_2_zm, & ! Variance of w (2nd PDF component) [m^2/s^2]
1426 : mixt_frac_zm ! Weight of 1st PDF component (Sk_w dependent) [-]
1427 :
1428 : type (stats_metadata_type), intent(in) :: &
1429 : stats_metadata
1430 :
1431 : !------------------------- Inout Variables -------------------------
1432 : type (stats), target, dimension(ngrdcol), intent(inout) :: &
1433 : stats_zm
1434 :
1435 : !------------------------- Output Variables -------------------------
1436 : integer, dimension(ngrdcol,nz), intent(out) :: &
1437 : low_lev_effect, & ! Index of lowest level that has an effect (for lev. k)
1438 : high_lev_effect ! Index of highest level that has an effect (for lev. k)
1439 :
1440 : !------------------------- Local Variables -------------------------
1441 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
1442 17870112 : vert_vel_up, & ! Average upwards vertical velocity component [m/s]
1443 17870112 : vert_vel_down, & ! Average downwards vertical velocity component [m/s]
1444 17870112 : w_min ! Minimum velocity to affect adjacent levels [m/s]
1445 :
1446 : real(kind = core_rknd ) :: &
1447 : dt_one_grid_lev, & ! Amount of time to travel one grid box [s]
1448 : dt_all_grid_levs, & ! Running count of amount of time taken to travel [s]
1449 : invrs_dt ! Inverse of timestep, used to reduce divides [1/s]
1450 :
1451 : integer :: k, i, j
1452 :
1453 : !------------------------- Begin Code -------------------------
1454 :
1455 : !$acc enter data create( vert_vel_up, vert_vel_down, w_min )
1456 :
1457 : if ( l_constant_thickness ) then ! thickness is a constant value.
1458 :
1459 : ! The value of w'x' may only be altered between levels 3 and gr%nz-2.
1460 : do k = 3, nz-2, 1
1461 : do i = 1, ngrdcol
1462 :
1463 : ! Compute the number of levels that effect the central thermodynamic
1464 : ! level through upwards motion (traveling from lower levels to reach
1465 : ! the central thermodynamic level).
1466 :
1467 : ! Start with the index of the thermodynamic level immediately below
1468 : ! the central thermodynamic level.
1469 : j = k - 1
1470 :
1471 : do ! loop downwards until answer is found.
1472 :
1473 : if ( gr%zt(i,k) - gr%zt(i,j) >= const_thick ) then
1474 :
1475 : ! Stop, the current grid level is the lowest level that can
1476 : ! be considered.
1477 : low_lev_effect(i,k) = j
1478 :
1479 : exit
1480 :
1481 : else
1482 :
1483 : ! Thermodynamic level 1 cannot be considered because it is
1484 : ! located below the surface or below the bottom of the model.
1485 : ! The lowest level that can be considered is thermodynamic
1486 : ! level 2.
1487 : if ( j == 2 ) then
1488 :
1489 : ! The current level (level 2) is the lowest level that can
1490 : ! be considered.
1491 : low_lev_effect(i,k) = j
1492 :
1493 : exit
1494 :
1495 : else
1496 :
1497 : ! Increment to the next vertical level down.
1498 : j = j - 1
1499 :
1500 : end if
1501 :
1502 : end if
1503 :
1504 : end do ! downwards loop
1505 :
1506 : end do
1507 : end do ! k = 3, gr%nz-2
1508 :
1509 : ! Compute the number of levels that effect the central thermodynamic
1510 : ! level through downwards motion (traveling from higher levels to
1511 : ! reach the central thermodynamic level).
1512 :
1513 : do k = 3, nz-2, 1
1514 : do i = 1, ngrdcol
1515 :
1516 : ! Start with the index of the thermodynamic level immediately above
1517 : ! the central thermodynamic level.
1518 : j = k + 1
1519 :
1520 : do ! loop upwards until answer is found.
1521 :
1522 : if ( gr%zt(i,j) - gr%zt(i,k) >= const_thick ) then
1523 :
1524 : ! Stop, the current grid level is the highest level that can
1525 : ! be considered.
1526 : high_lev_effect(i,k) = j
1527 :
1528 : exit
1529 :
1530 : else
1531 :
1532 : ! The highest level that can be considered is thermodynamic
1533 : ! level gr%nz.
1534 : if ( j == nz ) then
1535 :
1536 : ! The current level (level gr%nz) is the highest level
1537 : ! that can be considered.
1538 : high_lev_effect(i,k) = j
1539 :
1540 : exit
1541 :
1542 : else
1543 :
1544 : ! Increment to the next vertical level up.
1545 : j = j + 1
1546 :
1547 : end if
1548 :
1549 : end if
1550 :
1551 : end do ! upwards loop
1552 :
1553 : end do
1554 : end do ! k = 3, gr%nz-2
1555 :
1556 : else ! thickness based on vertical velocity and time step length.
1557 :
1558 8935056 : invrs_dt = 1.0_core_rknd / dt
1559 :
1560 : !$acc parallel loop gang vector collapse(2) default(present)
1561 768414816 : do k = 1, nz
1562 12690480816 : do i = 1, ngrdcol
1563 12681545760 : w_min(i,k) = gr%dzm(i,k) * invrs_dt
1564 : end do
1565 : end do
1566 : !$acc end parallel loop
1567 :
1568 : ! Find the average upwards vertical velocity and the average downwards
1569 : ! vertical velocity.
1570 : ! Note: A level that has all vertical wind moving downwards will have a
1571 : ! vert_vel_up value that is 0, and vice versa.
1572 : call mean_vert_vel_up_down( nz, ngrdcol, &
1573 : w_1_zm, w_2_zm, varnce_w_1_zm, varnce_w_2_zm, & ! In
1574 : mixt_frac_zm, 0.0_core_rknd, w_min, & ! In
1575 : stats_metadata, & ! In
1576 : stats_zm, & ! intent(inout)
1577 8935056 : vert_vel_down, vert_vel_up ) ! intent(out)
1578 :
1579 : ! The value of w'x' may only be altered between levels 3 and gr%nz-2.
1580 : !$acc parallel loop gang vector collapse(2) default(present)
1581 732674592 : do k = 3, nz-2, 1
1582 12093702192 : do i = 1, ngrdcol
1583 :
1584 : ! Compute the number of levels that effect the central thermodynamic
1585 : ! level through upwards motion (traveling from lower levels to reach
1586 : ! the central thermodynamic level).
1587 :
1588 : ! Initialize the overall delta t counter to 0.
1589 11361027600 : dt_all_grid_levs = 0.0_core_rknd
1590 :
1591 : ! Start with the index of the thermodynamic level immediately below
1592 : ! the central thermodynamic level.
1593 12294598700 : do j = k-1, 2, -1
1594 :
1595 11501599853 : low_lev_effect(i,k) = j
1596 :
1597 : ! Continue if there is some component of upwards vertical velocity.
1598 12831868236 : if ( vert_vel_up(i,j) > 0.0_core_rknd ) then
1599 :
1600 : ! Compute the amount of time it takes to travel one grid level
1601 : ! upwards: delta_t = delta_z / vert_vel_up.
1602 1470840636 : dt_one_grid_lev = gr%dzm(i,j) / vert_vel_up(i,j)
1603 :
1604 :
1605 : ! Total time elapsed for crossing all grid levels that have been
1606 : ! passed, thus far.
1607 1470840636 : dt_all_grid_levs = dt_all_grid_levs + dt_one_grid_lev
1608 :
1609 : ! Stop if has taken more than one model time step (overall) to
1610 : ! travel the entire extent of the current vertical grid level.
1611 1470840636 : if ( dt_all_grid_levs >= dt ) then
1612 :
1613 : ! The current level is the lowest level that can be
1614 : ! considered.
1615 : exit
1616 :
1617 : endif
1618 :
1619 : ! Stop if there isn't a component of upwards vertical velocity.
1620 : else
1621 :
1622 : ! The current level cannot be considered. The lowest level that
1623 : ! can be considered is one-level-above the current level.
1624 10030759217 : low_lev_effect(i,k) = j + 1
1625 :
1626 10030759217 : exit
1627 :
1628 : endif
1629 :
1630 : enddo ! downwards loop
1631 :
1632 : end do
1633 : enddo ! k = 3, gr%nz-2
1634 : !$acc end parallel loop
1635 :
1636 :
1637 : ! Compute the number of levels that effect the central thermodynamic
1638 : ! level through downwards motion (traveling from higher levels to
1639 : ! reach the central thermodynamic level).
1640 :
1641 : !$acc parallel loop gang vector collapse(2) default(present)
1642 732674592 : do k = 3, nz-2, 1
1643 12093702192 : do i = 1, ngrdcol
1644 :
1645 : ! Initialize the overall delta t counter to 0.
1646 11361027600 : dt_all_grid_levs = 0.0_core_rknd
1647 :
1648 : ! Start with the index of the thermodynamic level immediately above
1649 : ! the central thermodynamic level.
1650 12219259704 : do j = k+1, nz
1651 :
1652 11495378304 : high_lev_effect(i,k) = j
1653 :
1654 : ! Continue if there is some component of downwards vertical velocity.
1655 12416836026 : if ( vert_vel_down(i,j-1) < 0.0_core_rknd ) then
1656 :
1657 : ! Compute the amount of time it takes to travel one grid level
1658 : ! downwards: delta_t = - delta_z / vert_vel_down.
1659 : ! Note: There is a (-) sign in front of delta_z because the
1660 : ! distance traveled is downwards. Since vert_vel_down
1661 : ! has a negative value, dt_one_grid_lev will be a
1662 : ! positive value.
1663 1055808426 : dt_one_grid_lev = -gr%dzm(i,j-1) / vert_vel_down(i,j-1)
1664 :
1665 : ! Total time elapsed for crossing all grid levels that have been
1666 : ! passed, thus far.
1667 1055808426 : dt_all_grid_levs = dt_all_grid_levs + dt_one_grid_lev
1668 :
1669 : ! Stop if has taken more than one model time step (overall) to
1670 : ! travel the entire extent of the current vertical grid level.
1671 1055808426 : if ( dt_all_grid_levs >= dt ) then
1672 :
1673 : ! The current level is the highest level that can be
1674 : ! considered.
1675 : exit
1676 :
1677 : endif
1678 :
1679 : ! Stop if there isn't a component of downwards vertical velocity.
1680 : else
1681 :
1682 : ! The current level cannot be considered. The highest level
1683 : ! that can be considered is one-level-below the current level.
1684 10439569878 : high_lev_effect(i,k) = j - 1
1685 :
1686 10439569878 : exit
1687 :
1688 : end if
1689 :
1690 : end do ! upwards loop
1691 :
1692 : end do
1693 : enddo ! k = 3, gr%nz-2
1694 : !$acc end parallel loop
1695 :
1696 : end if ! l_constant_thickness
1697 :
1698 :
1699 : ! Information for levels 1, 2, gr%nz-1, and gr%nz is not needed.
1700 : ! However, set the values at these levels for purposes of not having odd
1701 : ! values in the arrays.
1702 : !$acc parallel loop gang vector default(present)
1703 149194656 : do i = 1, ngrdcol
1704 140259600 : low_lev_effect(i,1) = 1
1705 140259600 : high_lev_effect(i,1) = 1
1706 140259600 : low_lev_effect(i,2) = 2
1707 140259600 : high_lev_effect(i,2) = 2
1708 140259600 : low_lev_effect(i,nz-1) = nz-1
1709 140259600 : high_lev_effect(i,nz-1) = nz
1710 140259600 : low_lev_effect(i,nz) = nz
1711 149194656 : high_lev_effect(i,nz) = nz
1712 : end do
1713 : !$acc end parallel loop
1714 :
1715 : !$acc exit data delete( vert_vel_up, vert_vel_down, w_min )
1716 :
1717 8935056 : return
1718 :
1719 : end subroutine calc_turb_adv_range
1720 :
1721 : !=============================================================================
1722 8935056 : subroutine mean_vert_vel_up_down( nz, ngrdcol, &
1723 8935056 : w_1_zm, w_2_zm, varnce_w_1_zm, varnce_w_2_zm, &
1724 8935056 : mixt_frac_zm, w_ref, w_min, &
1725 : stats_metadata, &
1726 8935056 : stats_zm, &
1727 8935056 : mean_w_down, mean_w_up )
1728 :
1729 : ! Description
1730 : ! The values of vertical velocity, along a horizontal plane at any given
1731 : ! vertical level, are not allowed by CLUBB to be uniform. In other words,
1732 : ! there must be some variance in vertical velocity. This subroutine
1733 : ! calculates the mean of all values of vertical velocity, at any given
1734 : ! vertical level, that are greater than a certain reference velocity. This
1735 : ! subroutine also calculates the mean of all values of vertical velocity, at
1736 : ! any given vertical level, that are less than a certain reference velocity.
1737 : ! The reference velocity is usually 0 m/s, in which case this subroutine
1738 : ! calculates the average positive (upward) velocity and the average negative
1739 : ! (downward) velocity. However, the reference velocity may be other values,
1740 : ! such as wm_zm, which is the overall mean vertical velocity. If the
1741 : ! reference velocity is wm_zm, this subroutine calculates the average of all
1742 : ! values of w that are on the positive ("upward") side of the mean and the
1743 : ! average of all values of w that are on the negative ("downward") side of
1744 : ! the mean. These mean positive and negative vertical velocities are useful
1745 : ! in determining how long, on average, it takes a parcel of air, being
1746 : ! driven by subgrid updrafts or downdrafts, to traverse the length of the
1747 : ! vertical grid level.
1748 : !
1749 : ! Method
1750 : ! ------
1751 : !
1752 : ! The CLUBB model uses a joint PDF of vertical velocity, liquid water
1753 : ! potential temperature, and total water mixing ratio to determine subgrid
1754 : ! variability.
1755 : !
1756 : ! The values of vertical velocity, w, along an undefined horizontal plane
1757 : ! at any vertical level, are considered to approximately follow a
1758 : ! distribution that is a mixture of two normal (or Gaussian) distributions.
1759 : ! The values of w that are a part of the 1st normal distribution are
1760 : ! referred to as w_1, and the values of w that are part of the 2nd normal
1761 : ! distribution are referred to as w_2. Note that these distributions
1762 : ! overlap, and there are many values of w that are found in both w_1 and w_2.
1763 : !
1764 : ! The probability density function (PDF) for w, P(w), is:
1765 : !
1766 : ! P(w) = mixt_frac*P(w_1) + (1-mixt_frac)*P(w_2);
1767 : !
1768 : ! where "mixt_frac" is the weight of the 1st normal distribution, and P(w_1) and
1769 : ! P(w_2) are the equations for the 1st and 2nd normal distributions,
1770 : ! respectively:
1771 : !
1772 : ! P(w_1) = 1 / ( sigma_w_1 * sqrt(2*PI) )
1773 : ! * EXP[ -(w_1-mu_w_1)^2 / (2*sigma_w_1^2) ]; and
1774 : !
1775 : ! P(w_2) = 1 / ( sigma_w_2 * sqrt(2*PI) )
1776 : ! * EXP[ -(w_2-mu_w_2)^2 / (2*sigma_w_2^2) ].
1777 : !
1778 : ! The mean of the 1st normal distribution is mu_w_1, and the standard
1779 : ! deviation of the 1st normal distribution is sigma_w_1. The mean of the
1780 : ! 2nd normal distribution is mu_w_2, and the standard deviation of the 2nd
1781 : ! normal distribution is sigma_w_2.
1782 : !
1783 : ! The average value of w, distributed according to the probability
1784 : ! distribution, between limits alpha and beta, is:
1785 : !
1786 : ! <w|_(alpha:beta)> = INT(alpha:beta) w P(w) dw.
1787 : !
1788 : ! The average value of w over a certain domain is used to determine the
1789 : ! average positive and negative (as compared to the reference velocity)
1790 : ! values of w at any vertical level.
1791 : !
1792 : ! Average Negative Vertical Velocity
1793 : ! ----------------------------------
1794 : !
1795 : ! The average of all values of w in the distribution that are below the
1796 : ! reference velocity, w|_ref, is the mean value of w over the domain
1797 : ! -inf <= w <= w|_ref, such that:
1798 : !
1799 : ! <w|_(-inf:w|_ref)> = INT(-inf:w|_ref) w P(w) dw.
1800 : ! = mixt_frac * INT(-inf:w|_ref) w_1 P(w_1) dw_1
1801 : ! + (1-mixt_frac) * INT(-inf:w|_ref) w_2 P(w_2) dw_2.
1802 : !
1803 : ! For each normal distribution in the mixture of normal distribution, i
1804 : ! (where "i" can be 1 or 2):
1805 : !
1806 : ! INT(-inf:w|_ref) wi P(wi) dwi =
1807 : ! - ( sigma_wi / sqrt(2*PI) ) * EXP[ -(w|_ref-mu_wi)^2 / (2*sigma_wi^2) ]
1808 : ! + mu_wi * (1/2)*[ 1 + erf( (w|_ref-mu_wi) / (sqrt(2)*sigma_wi) ) ];
1809 : !
1810 : ! where mu_wi is the mean of w for the ith normal distribution, sigma_wi is
1811 : ! the standard deviations of w for the ith normal distribution, and erf( )
1812 : ! is the error function.
1813 : !
1814 : ! The mean of all values of w <= w|_ref is:
1815 : !
1816 : ! <w|_(-inf:w|_ref)> =
1817 : ! mixt_frac * { - ( sigma_w_1 / sqrt(2*PI) )
1818 : ! * EXP[ -(w|_ref-mu_w_1)^2 / (2*sigma_w_1^2) ]
1819 : ! + mu_w_1 * (1/2)
1820 : ! *[1 + erf( (w|_ref-mu_w_1) / (sqrt(2)*sigma_w_1) )] }
1821 : ! + (1-mixt_frac) * { - ( sigma_w_2 / sqrt(2*PI) )
1822 : ! * EXP[ -(w|_ref-mu_w_2)^2 / (2*sigma_w_2^2) ]
1823 : ! + mu_w_2 * (1/2)
1824 : ! *[1 + erf( (w|_ref-mu_w_2) / (sqrt(2)*sigma_w_2) )] }.
1825 : !
1826 : ! Average Positive Vertical Velocity
1827 : ! ----------------------------------
1828 : !
1829 : ! The average of all values of w in the distribution that are above the
1830 : ! reference velocity, w|_ref, is the mean value of w over the domain
1831 : ! w|_ref <= w <= inf, such that:
1832 : !
1833 : ! <w|_(w|_ref:inf)> = INT(w|_ref:inf) w P(w) dw.
1834 : ! = mixt_frac * INT(w|_ref:inf) w_1 P(w_1) dw_1
1835 : ! + (1-mixt_frac) * INT(w|_ref:inf) w_2 P(w_2) dw_2.
1836 : !
1837 : ! For each normal distribution in the mixture of normal distribution, i
1838 : ! (where "i" can be 1 or 2):
1839 : !
1840 : ! INT(w|_ref:inf) wi P(wi) dwi =
1841 : ! ( sigma_wi / sqrt(2*PI) ) * EXP[ -(w|_ref-mu_wi)^2 / (2*sigma_wi^2) ]
1842 : ! + mu_wi * (1/2)*[ 1 - erf( (w|_ref-mu_wi) / (sqrt(2)*sigma_wi) ) ];
1843 : !
1844 : ! where mu_wi is the mean of w for the ith normal distribution, sigma_wi is
1845 : ! the standard deviations of w for the ith normal distribution, and erf( )
1846 : ! is the error function.
1847 : !
1848 : ! The mean of all values of w >= w|_ref is:
1849 : !
1850 : ! <w|_(w|_ref:inf)> =
1851 : ! mixt_frac * { ( sigma_w_1 / sqrt(2*PI) )
1852 : ! * EXP[ -(w|_ref-mu_w_1)^2 / (2*sigma_w_1^2) ]
1853 : ! + mu_w_1 * (1/2)
1854 : ! *[1 - erf( (w|_ref-mu_w_1) / (sqrt(2)*sigma_w_1) )] }
1855 : ! + (1-mixt_frac) * { ( sigma_w_2 / sqrt(2*PI) )
1856 : ! * EXP[ -(w|_ref-mu_w_2)^2 / (2*sigma_w_2^2) ]
1857 : ! + mu_w_2 * (1/2)
1858 : ! *[1 - erf( (w|_ref-mu_w_2) / (sqrt(2)*sigma_w_2) )] }.
1859 : !
1860 : ! Special Limitations:
1861 : ! --------------------
1862 : !
1863 : ! A normal distribution has a domain from -inf to inf. However, the mixture
1864 : ! of normal distributions is an approximation of the distribution of values
1865 : ! of w along a horizontal plane at any given vertical level. Vertical
1866 : ! velocity, w, has absolute minimum and maximum values (that cannot be
1867 : ! predicted by the PDF). The absolute maximum and minimum for each normal
1868 : ! distribution is most likely found within 2 or 3 standard deviations of the
1869 : ! mean for the relevant normal distribution. In other words, for each
1870 : ! normal distribution in the mixture of normal distributions, all the values
1871 : ! of w are found within 2 or 3 standard deviations on both sides of the
1872 : ! mean. Therefore, if one (or both) of the normal distributions has a mean
1873 : ! that is more than 3 standard deviations away from the reference velocity,
1874 : ! then that entire w distribution is found on ONE side of the reference
1875 : ! velocity.
1876 : !
1877 : ! Therefore:
1878 : !
1879 : ! a) where mu_wi + 3*sigma_wi <= w|_ref:
1880 : !
1881 : ! The entire ith normal distribution of w is on the negative side of
1882 : ! w|_ref; and
1883 : !
1884 : ! INT(-inf:w|_ref) wi P(wi) dwi = mu_wi; and
1885 : ! INT(inf:w|_ref) wi P(wi) dwi = 0.
1886 : !
1887 : ! b) where mu_wi - 3*sigma_wi >= w|_ref:
1888 : !
1889 : ! The entire ith normal distribution of w is on the positive side of
1890 : ! w|_ref; and
1891 : !
1892 : ! INT(-inf:w|_ref) wi P(wi) dwi = 0; and
1893 : ! INT(inf:w|_ref) wi P(wi) dwi = mu_wi.
1894 : !
1895 : ! Notes: A value of 3 standard deviations above and below the mean of the
1896 : ! ith normal distribution was chosen for the approximate maximum and
1897 : ! minimum values of the ith normal distribution because 99.7% of
1898 : ! values in a normal distribution are found within 3 standard
1899 : ! deviations from the mean (compared to 95.4% for 2 standard
1900 : ! deviations). The value of 3 standard deviations provides for a
1901 : ! reasonable estimate of the absolute maximum and minimum of w, while
1902 : ! covering a great majority of the normal distribution.
1903 : !
1904 : ! In addition to approximating the up and down components of w
1905 : ! by checking if the pdfs are greater than 3 standard deviations
1906 : ! from the mean, there is now a case to approximate when w is
1907 : ! too small in general. The input array, w_min, contains the
1908 : ! minimum values of vertical velocity that would be required
1909 : ! at a given grid level for that grid box to be able to affect
1910 : ! the adjacent levels. If the magnitude of w at a given level
1911 : ! is less than 3 standard deviations below w_min for that level,
1912 : ! then there is no significant portion of the air from that grid
1913 : ! box that is capable of interacting with the next level, and
1914 : ! the upward and downward components for that pdf are set to 0.
1915 : !
1916 : ! References:
1917 : !-----------------------------------------------------------------------
1918 :
1919 : use grid_class, only: &
1920 : grid ! Type
1921 :
1922 : use stats_type_utilities, only: &
1923 : stat_update_var ! Procedure(s)
1924 :
1925 : use stats_variables, only: &
1926 : stats_metadata_type
1927 :
1928 : use clubb_precision, only: &
1929 : core_rknd ! Variable(s)
1930 :
1931 : use stats_type, only: stats ! Type
1932 :
1933 : implicit none
1934 :
1935 : !------------------------- Input Variables -------------------------
1936 : integer, intent(in) :: &
1937 : nz, &
1938 : ngrdcol
1939 :
1940 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
1941 : w_1_zm, & ! Mean w (1st PDF component) [m/s]
1942 : w_2_zm, & ! Mean w (2nd PDF component) [m/s]
1943 : varnce_w_1_zm, & ! Variance of w (1st PDF component) [m^2/s^2]
1944 : varnce_w_2_zm, & ! Variance of w (2nd PDF component) [m^2/s^2]
1945 : mixt_frac_zm, & ! Weight of 1st PDF component (Sk_w dependent) [-]
1946 : w_min ! Minimum velocity to affect adjacent level [m/s]
1947 :
1948 : real( kind = core_rknd ), intent(in) :: &
1949 : w_ref ! Reference velocity, w|_ref (normally = 0) [m/s]
1950 :
1951 : type (stats_metadata_type), intent(in) :: &
1952 : stats_metadata
1953 :
1954 : !------------------------- Inout Variables -------------------------
1955 : type (stats), target, dimension(ngrdcol), intent(inout) :: &
1956 : stats_zm
1957 :
1958 : !------------------------- Output Variables -------------------------
1959 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
1960 : mean_w_down, & ! Overall mean w (<= w|_ref) [m/s]
1961 : mean_w_up ! Overall mean w (>= w|_ref) [m/s]
1962 :
1963 : !------------------------- Local Variables -------------------------
1964 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
1965 17870112 : mean_w_down_1st, & ! Mean w (<= w|_ref) from 1st normal distribution [m/s]
1966 17870112 : mean_w_down_2nd, & ! Mean w (<= w|_ref) from 2nd normal distribution [m/s]
1967 17870112 : mean_w_up_1st, & ! Mean w (>= w|_ref) from 1st normal distribution [m/s]
1968 17870112 : mean_w_up_2nd ! Mean w (>= w|_ref) from 2nd normal distribution [m/s]
1969 :
1970 : integer :: i, k
1971 :
1972 : !------------------------- Begin Code -------------------------
1973 :
1974 : !$acc enter data create( mean_w_down_1st, mean_w_down_2nd, mean_w_up_1st, mean_w_up_2nd )
1975 :
1976 : call calc_mean_w_up_down_component( nz, ngrdcol, & ! intent(in)
1977 : w_1_zm, varnce_w_1_zm, & ! intent(in)
1978 : w_ref, w_min, & ! intent(in)
1979 8935056 : mean_w_down_1st, mean_w_up_1st ) ! intent(out)
1980 :
1981 : call calc_mean_w_up_down_component( nz, ngrdcol, & ! intent(in)
1982 : w_2_zm, varnce_w_2_zm, & ! intent(in)
1983 : w_ref, w_min, & ! intent(in)
1984 8935056 : mean_w_down_2nd, mean_w_up_2nd ) ! intent(out)
1985 :
1986 : ! Overall mean of downwards w.
1987 : !$acc parallel loop gang vector collapse(2) default(present)
1988 768414816 : do k = 1, nz
1989 12690480816 : do i = 1, ngrdcol
1990 23844132000 : mean_w_down(i,k) = mixt_frac_zm(i,k) * mean_w_down_1st(i,k) &
1991 36525677760 : + ( 1.0_core_rknd - mixt_frac_zm(i,k) ) * mean_w_down_2nd(i,k)
1992 : end do
1993 : end do
1994 : !$acc end parallel loop
1995 :
1996 : ! Overall mean of upwards w.
1997 : !$acc parallel loop gang vector collapse(2) default(present)
1998 768414816 : do k = 1, nz
1999 12690480816 : do i = 1, ngrdcol
2000 23844132000 : mean_w_up(i,k) = mixt_frac_zm(i,k) * mean_w_up_1st(i,k) &
2001 36525677760 : + ( 1.0_core_rknd - mixt_frac_zm(i,k) ) * mean_w_up_2nd(i,k)
2002 : end do
2003 : end do
2004 : !$acc end parallel loop
2005 :
2006 8935056 : if ( stats_metadata%l_stats_samp ) then
2007 : !$acc update host( mean_w_up, mean_w_down )
2008 0 : do i = 1, ngrdcol
2009 0 : call stat_update_var( stats_metadata%imean_w_up, mean_w_up(i,:), & ! intent(in)
2010 0 : stats_zm(i) ) ! intent(inout)
2011 :
2012 : call stat_update_var( stats_metadata%imean_w_down, mean_w_down(i,:), & ! intent(in)
2013 0 : stats_zm(i) ) ! intent(inout)
2014 : end do
2015 : end if ! stats_metadata%l_stats_samp
2016 :
2017 : !$acc exit data delete( mean_w_down_1st, mean_w_down_2nd, mean_w_up_1st, mean_w_up_2nd )
2018 :
2019 8935056 : return
2020 :
2021 : end subroutine mean_vert_vel_up_down
2022 :
2023 : !=============================================================================
2024 17870112 : subroutine calc_mean_w_up_down_component( nz, ngrdcol, &
2025 17870112 : w_i_zm, varnce_w_i, &
2026 17870112 : w_ref, w_min, &
2027 17870112 : mean_w_down_i, mean_w_up_i )
2028 :
2029 : ! Description: This procedure is used to split the PDF component of
2030 : ! vertical velocity into upward and downward components.
2031 : !
2032 : ! The method used is described in the description of
2033 : ! mean_vert_vel_up_down, which calls this function.
2034 : !
2035 : ! Notes: The calculation has been updated to optionally use intel's
2036 : ! mkl_vml functions to allow vectorized calculations. Not all
2037 : ! grid levels require expensive calculations though, so the
2038 : ! strategy is as follows
2039 : ! 1. Keep track of which levels do need the calculation
2040 : ! 2. Store those the relavent values from those levels in
2041 : ! a contigous array
2042 : ! 3. Perform vectorized calculation on contiguous arrays
2043 : ! using mkl_vml functions
2044 : ! 4. Unpack results from contiguous array into output array
2045 : ! Enabling this faster version requires compilation with MKL, by
2046 : ! using -DMKL as a compiler flag
2047 : !
2048 : !-----------------------------------------------------------------------
2049 :
2050 : use grid_class, only: &
2051 : grid ! Type
2052 :
2053 : use constants_clubb, only: &
2054 : sqrt_2pi, & ! Constant(s)
2055 : sqrt_2, &
2056 : one
2057 : #ifdef MKL
2058 : use constants_clubb, only: &
2059 : one_half ! Constant(s)
2060 : #endif
2061 :
2062 : use clubb_precision, only: &
2063 : core_rknd ! Variable(s)
2064 :
2065 : implicit none
2066 :
2067 : integer, intent(in) :: &
2068 : nz, &
2069 : ngrdcol
2070 :
2071 : !------------------------- Input Variables -------------------------
2072 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
2073 : w_i_zm, & ! Mean of w [m/s]
2074 : varnce_w_i, & ! Variance of w [m^2/s^2]
2075 : w_min ! Minimum velocity required to affect adjacent level [m/s]
2076 :
2077 : real( kind = core_rknd ), intent(in) :: &
2078 : w_ref ! Reference velocity, w|_ref (normally = 0) [m/s]
2079 :
2080 : !------------------------- Output Variables -------------------------
2081 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
2082 : mean_w_down_i, & ! Mean w (<= w|_ref) from normal distribution [m/s]
2083 : mean_w_up_i ! Mean w (>= w|_ref) from normal distribution [m/s]
2084 :
2085 : !------------------------- Local Variables -------------------------
2086 : real( kind = core_rknd ) :: &
2087 : erf_cache, & ! erf/cdfnorm values
2088 : exp_cache ! exp() values
2089 :
2090 : real( kind = core_rknd ) :: &
2091 : sigma_w_i, & ! Variance of w (1st PDF component) [m^2/s^2]
2092 : invrs_sqrt_2pi ! The inverse of sqrt(2*pi), calculated to save divide operations
2093 :
2094 : integer :: i, k ! Vertical loop index
2095 :
2096 : !------------------------- Begin Code -------------------------
2097 :
2098 17870112 : invrs_sqrt_2pi = one / sqrt_2pi
2099 :
2100 : ! Loop over momentum levels from 2 to nz-1. Levels 1 and nz
2101 : ! are not needed.
2102 : !$acc parallel loop gang vector collapse(2) default(present)
2103 1501089408 : do k = 2, nz-1
2104 24784183008 : do i = 1, ngrdcol
2105 :
2106 : ! Standard deviation of w for the normal distribution.
2107 23283093600 : sigma_w_i = sqrt( varnce_w_i(i,k) )
2108 :
2109 24766312896 : if( abs( w_i_zm(i,k) ) + 3.0_core_rknd*sigma_w_i <= w_min(i,k) ) then
2110 :
2111 : ! The entire normal is too weak to affect adjacent grid levels
2112 : ! w is considered to be 0 in both up and down directions
2113 20908946477 : mean_w_down_i(i,k) = 0.0_core_rknd
2114 20908946477 : mean_w_up_i(i,k) = 0.0_core_rknd
2115 :
2116 2374147123 : elseif ( w_i_zm(i,k) + 3._core_rknd*sigma_w_i <= w_ref ) then
2117 :
2118 : ! The entire normal is on the negative side of w|_ref.
2119 9626801 : mean_w_down_i(i,k) = w_i_zm(i,k)
2120 9626801 : mean_w_up_i(i,k) = 0.0_core_rknd
2121 :
2122 2364520322 : elseif ( w_i_zm(i,k) - 3._core_rknd*sigma_w_i >= w_ref ) then
2123 :
2124 : ! The entire normal is on the positive side of w|_ref.
2125 363315792 : mean_w_down_i(i,k) = 0.0_core_rknd
2126 363315792 : mean_w_up_i(i,k) = w_i_zm(i,k)
2127 :
2128 : else
2129 :
2130 : ! The normal has significant values on both sides of w_ref.
2131 2001204530 : exp_cache = exp( -(w_ref-w_i_zm(i,k))**2 / (2.0_core_rknd*sigma_w_i**2) )
2132 :
2133 2001204530 : erf_cache = erf( (w_ref-w_i_zm(i,k)) / (sqrt_2*sigma_w_i ) )
2134 :
2135 : mean_w_down_i(i,k) = - sigma_w_i * invrs_sqrt_2pi * exp_cache &
2136 2001204530 : + w_i_zm(i,k) * 0.5_core_rknd*( 1.0_core_rknd + erf_cache)
2137 :
2138 : mean_w_up_i(i,k) = + sigma_w_i * invrs_sqrt_2pi * exp_cache &
2139 2001204530 : + w_i_zm(i,k) * 0.5_core_rknd*( 1.0_core_rknd - erf_cache)
2140 :
2141 : end if
2142 :
2143 : end do
2144 : end do ! k = 2, gr%nz
2145 : !$acc end parallel loop
2146 :
2147 : ! Upper and lower levels are not used, set to 0 to besafe and avoid NaN problems
2148 : !$acc parallel loop gang vector default(present)
2149 298389312 : do i = 1, ngrdcol
2150 280519200 : mean_w_down_i(i,1) = 0.0_core_rknd
2151 280519200 : mean_w_up_i(i,1) = 0.0_core_rknd
2152 :
2153 280519200 : mean_w_down_i(i,nz) = 0.0_core_rknd
2154 298389312 : mean_w_up_i(i,nz) = 0.0_core_rknd
2155 : end do
2156 : !$acc end parallel loop
2157 :
2158 17870112 : return
2159 :
2160 : end subroutine calc_mean_w_up_down_component
2161 :
2162 : !===============================================================================
2163 :
2164 : end module mono_flux_limiter
|