LCOV - code coverage report
Current view: top level - physics/clubb/src/CLUBB_core - mixing_length.F90 (source / functions) Hit Total Coverage
Test: coverage.info Lines: 176 402 43.8 %
Date: 2024-12-17 17:57:11 Functions: 2 3 66.7 %

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

Generated by: LCOV version 1.14