LCOV - code coverage report
Current view: | top level - physics/clubb/src/CLUBB_core - mono_flux_limiter.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 259 | 327 | 79.2 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 7 | 7 | 100.0 %

              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

Generated by: LCOV version 1.14