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