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