LCOV - code coverage report
Current view: | top level - physics/clubb/src/CLUBB_core - advance_helper_module.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 147 | 336 | 43.8 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 6 | 22 | 27.3 %

              Line data    Source code
           1             :
           2             : module advance_helper_module
           3             :
           4             : ! Description:
           5             : !   This module contains helper methods for the advance_* modules.
           6             : !------------------------------------------------------------------------
           7             :
           8             :   implicit none
           9             :
          10             :   public :: &
          11             :     set_boundary_conditions_lhs, &
          12             :     set_boundary_conditions_rhs, &
          13             :     calc_stability_correction,   &
          14             :     calc_brunt_vaisala_freq_sqd, &
          15             :     compute_Cx_fnc_Richardson, &
          16             :     wp2_term_splat_lhs, &
          17             :     wp3_term_splat_lhs, &
          18             :     smooth_min, smooth_max, &
          19             :     smooth_heaviside_peskin, &
          20             :     calc_xpwp, &
          21             :     vertical_avg, &
          22             :     vertical_integral, &
          23             :     Lscale_width_vert_avg
          24             :
          25             :   interface calc_xpwp
          26             :     module procedure calc_xpwp_1D
          27             :     module procedure calc_xpwp_2D
          28             :   end interface
          29             :
          30             :   private ! Set Default Scope
          31             :
          32             : !===============================================================================
          33             :   interface smooth_min
          34             :
          35             :     ! These functions smooth the output of the min function
          36             :     ! by introducing a varyingly steep path between the two input variables.
          37             :     ! The degree to which smoothing is applied depends on the value of 'smth_coef'.
          38             :     ! If 'smth_coef' goes toward 0, the output of the min function will be
          39             :     !        0.5 * ((a+b) - abs(a-b))
          40             :     ! If a > b, then this comes out to be b. Likewise, if a < b, abs(a-b)=b-a so we get a.
          41             :     ! Increasing the smoothing coefficient will lead to a greater degree of smoothing
          42             :     ! in the smooth min and max functions. Generally, the coefficient should roughly scale
          43             :     ! with the magnitude of data in the data structure that is to be smoothed, in order to
          44             :     ! obtain a sensible degree of smoothing (not too much, not too little).
          45             :
          46             :     module procedure smooth_min_sclr_idx
          47             :     module procedure smooth_min_array_scalar
          48             :     module procedure smooth_min_arrays
          49             :     module procedure smooth_min_scalars
          50             :
          51             :   end interface
          52             :
          53             : !===============================================================================
          54             :   interface smooth_max
          55             :
          56             :     ! These functions smooth the output of the max functions
          57             :     ! by introducing a varyingly steep path between the two input variables.
          58             :     ! The degree to which smoothing is applied depends on the value of 'smth_coef'.
          59             :     ! If 'smth_coef' goes toward 0, the output of the max function will be
          60             :     !        0.5 * ((a+b) + abs(a-b))
          61             :     ! If a > b, then this comes out to be a. Likewise, if a < b, abs(a-b)=b-a so we get b.
          62             :     ! Increasing the smoothing coefficient will lead to a greater degree of smoothing
          63             :     ! in the smooth min and max functions. Generally, the coefficient should roughly scale
          64             :     ! with the magnitude of data in the data structure that is to be smoothed, in order to
          65             :     ! obtain a sensible degree of smoothing (not too much, not too little).
          66             :
          67             :     module procedure smooth_max_sclr_idx
          68             :     module procedure smooth_max_array_scalar
          69             :     module procedure smooth_max_arrays
          70             :     module procedure smooth_max_array_1D_scalar
          71             :     module procedure smooth_max_scalars
          72             :
          73             :   end interface
          74             :
          75             : !===============================================================================
          76             :   contains
          77             :
          78             :   !---------------------------------------------------------------------------
          79           0 :   subroutine set_boundary_conditions_lhs( diag_index, low_bound, high_bound, &
          80           0 :                                           lhs, &
          81             :                                           diag_index2, low_bound2, high_bound2 )
          82             :
          83             :   ! Description:
          84             :   !   Sets the boundary conditions for a left-hand side LAPACK matrix.
          85             :   !
          86             :   ! References:
          87             :   !   none
          88             :   !---------------------------------------------------------------------------
          89             :
          90             :     use clubb_precision, only: &
          91             :         core_rknd ! Variable(s)
          92             :
          93             :     use constants_clubb, only: &
          94             :         one, zero
          95             :
          96             :     implicit none
          97             :
          98             :     ! Exernal
          99             :     intrinsic :: present
         100             :
         101             :     ! Input Variables
         102             :     integer, intent(in) :: &
         103             :       diag_index, low_bound, high_bound ! boundary indexes for the first variable
         104             :
         105             :     ! Input / Output Variables
         106             :     real( kind = core_rknd ), dimension(:,:), intent(inout) :: &
         107             :       lhs ! left hand side of the LAPACK matrix equation
         108             :
         109             :     ! Optional Input Variables
         110             :     integer, intent(in), optional :: &
         111             :       diag_index2, low_bound2, high_bound2 ! boundary indexes for the second variable
         112             :
         113             :     ! --------------------- BEGIN CODE ----------------------
         114             :
         115           0 :     if ( ( present( low_bound2 ) .or. present( high_bound2 ) ) .and. &
         116             :          ( .not. present( diag_index2 ) ) ) then
         117             :
         118           0 :       error stop "Boundary index provided without diag_index."
         119             :
         120             :     end if
         121             :
         122             :     ! Set the lower boundaries for the first variable
         123           0 :     lhs(:,low_bound) = zero
         124           0 :     lhs(diag_index,low_bound) = one
         125             :
         126             :     ! Set the upper boundaries for the first variable
         127           0 :     lhs(:,high_bound) = zero
         128           0 :     lhs(diag_index,high_bound) = one
         129             :
         130             :     ! Set the lower boundaries for the second variable, if it is provided
         131           0 :     if ( present( low_bound2 ) ) then
         132             :
         133           0 :       lhs(:,low_bound2) = zero
         134           0 :       lhs(diag_index2,low_bound2) = one
         135             :
         136             :     end if
         137             :
         138             :     ! Set the upper boundaries for the second variable, if it is provided
         139           0 :     if ( present( high_bound2 ) ) then
         140             :
         141           0 :       lhs(:,high_bound2) = zero
         142           0 :       lhs(diag_index2,high_bound2) = one
         143             :
         144             :     end if
         145             :
         146           0 :     return
         147             :   end subroutine set_boundary_conditions_lhs
         148             :
         149             :   !--------------------------------------------------------------------------
         150           0 :   subroutine set_boundary_conditions_rhs( &
         151             :                low_value, low_bound, high_value, high_bound, &
         152           0 :                rhs, &
         153             :                low_value2, low_bound2, high_value2, high_bound2 )
         154             :
         155             :   ! Description:
         156             :   !   Sets the boundary conditions for a right-hand side LAPACK vector.
         157             :   !
         158             :   ! References:
         159             :   !   none
         160             :   !---------------------------------------------------------------------------
         161             :
         162             :     use clubb_precision, only: &
         163             :         core_rknd ! Variable(s)
         164             :
         165             :     implicit none
         166             :
         167             :     ! Exernal
         168             :     intrinsic :: present
         169             :
         170             :     ! Input Variables
         171             :
         172             :     ! The values for the first variable
         173             :     real( kind = core_rknd ), intent(in) :: low_value, high_value
         174             :
         175             :     ! The bounds for the first variable
         176             :     integer, intent(in) :: low_bound, high_bound
         177             :
         178             :     ! Input / Output Variables
         179             :
         180             :     ! The right-hand side vector
         181             :     real( kind = core_rknd ), dimension(:), intent(inout) :: rhs
         182             :
         183             :     ! Optional Input Variables
         184             :
         185             :     ! The values for the second variable
         186             :     real( kind = core_rknd ), intent(in), optional :: low_value2, high_value2
         187             :
         188             :     ! The bounds for the second variable
         189             :     integer, intent(in), optional :: low_bound2, high_bound2
         190             :
         191             :
         192             :     ! -------------------- BEGIN CODE ------------------------
         193             :
         194             :     ! Stop execution if a boundary was provided without a value
         195           0 :     if ( (present( low_bound2 ) .and. (.not. present( low_value2 ))) .or. &
         196             :          (present( high_bound2 ) .and. (.not. present( high_value2 ))) ) then
         197             :
         198           0 :       error stop "Boundary condition provided without value."
         199             :
         200             :     end if
         201             :
         202             :     ! Set the lower and upper bounds for the first variable
         203           0 :     rhs(low_bound) = low_value
         204           0 :     rhs(high_bound) = high_value
         205             :
         206             :     ! If a lower bound was given for the second variable, set it
         207           0 :     if ( present( low_bound2 ) ) then
         208           0 :       rhs(low_bound2) = low_value2
         209             :     end if
         210             :
         211             :     ! If an upper bound was given for the second variable, set it
         212           0 :     if ( present( high_bound2 ) ) then
         213           0 :       rhs(high_bound2) = high_value2
         214             :     end if
         215             :
         216           0 :     return
         217             :   end subroutine set_boundary_conditions_rhs
         218             :
         219             :   !===============================================================================
         220     8935056 :   subroutine calc_stability_correction( nz, ngrdcol, gr, &
         221     8935056 :                                         thlm, Lscale_zm, em, &
         222     8935056 :                                         exner, rtm, rcm, &
         223     8935056 :                                         p_in_Pa, thvm, ice_supersat_frac, &
         224     8935056 :                                         lambda0_stability_coef, &
         225     8935056 :                                         bv_efold, &
         226             :                                         saturation_formula, &
         227             :                                         l_brunt_vaisala_freq_moist, &
         228             :                                         l_use_thvm_in_bv_freq, &
         229     8935056 :                                         stability_correction )
         230             :   !
         231             :   ! Description:
         232             :   !   Stability Factor
         233             :   !
         234             :   ! References:
         235             :   !
         236             :   !--------------------------------------------------------------------
         237             :
         238             :     use constants_clubb, only: &
         239             :         zero, one, three    ! Constant(s)
         240             :
         241             :     use grid_class, only:  &
         242             :         grid    ! Type
         243             :
         244             :     use clubb_precision, only:  &
         245             :         core_rknd ! Variable(s)
         246             :
         247             :     implicit none
         248             :
         249             :     ! ---------------- Input Variables ----------------
         250             :     integer, intent(in) :: &
         251             :       nz, &
         252             :       ngrdcol
         253             :
         254             :     type (grid), target, intent(in) :: gr
         255             :
         256             :     real( kind = core_rknd ), intent(in), dimension(ngrdcol,nz) :: &
         257             :       Lscale_zm,       & ! Turbulent mixing length interp to m levs  [m]
         258             :       em,              & ! Turbulent Kinetic Energy (TKE)            [m^2/s^2]
         259             :       thlm,            & ! th_l (thermo. levels)                     [K]
         260             :       exner,           & ! Exner function                            [-]
         261             :       rtm,             & ! total water mixing ratio, r_t             [kg/kg]
         262             :       rcm,             & ! cloud water mixing ratio, r_c             [kg/kg]
         263             :       p_in_Pa,         & ! Air pressure                              [Pa]
         264             :       thvm,            & ! Virtual potential temperature             [K]
         265             :       ice_supersat_frac
         266             :
         267             :     real( kind = core_rknd ), dimension(ngrdcol), intent(in) :: &
         268             :       lambda0_stability_coef, &     ! CLUBB tunable parameter lambda0_stability_coef
         269             :       bv_efold                      ! Control parameter for inverse e-folding of
         270             :                                     ! cloud fraction in the mixed Brunt Vaisala frequency
         271             :
         272             :     integer, intent(in) :: &
         273             :       saturation_formula ! Integer that stores the saturation formula to be used
         274             :
         275             :     logical, intent(in) :: &
         276             :       l_brunt_vaisala_freq_moist, & ! Use a different formula for the Brunt-Vaisala frequency in
         277             :                                     ! saturated atmospheres (from Durran and Klemp, 1982)
         278             :       l_use_thvm_in_bv_freq         ! Use thvm in the calculation of Brunt-Vaisala frequency
         279             :
         280             :     ! ---------------- Output Variables ----------------
         281             :     real( kind = core_rknd ), intent(out), dimension(ngrdcol,nz) :: &
         282             :       stability_correction
         283             :
         284             :     ! ---------------- Local Variables ----------------
         285             :     real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
         286    17870112 :       brunt_vaisala_freq_sqd, & !  []
         287    17870112 :       brunt_vaisala_freq_sqd_mixed, &
         288    17870112 :       brunt_vaisala_freq_sqd_dry, & !  []
         289    17870112 :       brunt_vaisala_freq_sqd_moist, &
         290    17870112 :       lambda0_stability
         291             :
         292             :     integer :: i, k
         293             :
         294             :     !------------ Begin Code --------------
         295             :
         296             :     !$acc enter data create( brunt_vaisala_freq_sqd, brunt_vaisala_freq_sqd_mixed, &
         297             :     !$acc                    brunt_vaisala_freq_sqd_moist, brunt_vaisala_freq_sqd_dry, &
         298             :     !$acc                    lambda0_stability )
         299             :
         300             :     call calc_brunt_vaisala_freq_sqd( nz, ngrdcol, gr, thlm, &          ! intent(in)
         301             :                                       exner, rtm, rcm, p_in_Pa, thvm, & ! intent(in)
         302             :                                       ice_supersat_frac, &              ! intent(in)
         303             :                                       saturation_formula, &             ! intent(in)
         304             :                                       l_brunt_vaisala_freq_moist, &     ! intent(in)
         305             :                                       l_use_thvm_in_bv_freq, &          ! intent(in)
         306             :                                       bv_efold, &                       ! intent(in)
         307             :                                       brunt_vaisala_freq_sqd, &         ! intent(out)
         308             :                                       brunt_vaisala_freq_sqd_mixed,&    ! intent(out)
         309             :                                       brunt_vaisala_freq_sqd_dry, &     ! intent(out)
         310     8935056 :                                       brunt_vaisala_freq_sqd_moist )    ! intent(out)
         311             :
         312             :     !$acc parallel loop gang vector collapse(2) default(present)
         313   768414816 :     do k = 1, nz
         314 12690480816 :       do i = 1, ngrdcol
         315 12681545760 :         if ( brunt_vaisala_freq_sqd(i,k) > zero  ) then
         316 11650400594 :           lambda0_stability(i,k) = lambda0_stability_coef(i)
         317             :         else
         318   271665406 :           lambda0_stability(i,k) = zero
         319             :         end if
         320             :       end do
         321             :     end do
         322             :     !$acc end parallel loop
         323             :
         324             :     !$acc parallel loop gang vector collapse(2) default(present)
         325   768414816 :     do k = 1, nz
         326 12690480816 :       do i = 1, ngrdcol
         327 23844132000 :         stability_correction(i,k) = one + min( lambda0_stability(i,k) * brunt_vaisala_freq_sqd(i,k) &
         328 36525677760 :                                                * Lscale_zm(i,k)**2 / em(i,k), three )
         329             :       end do
         330             :     end do
         331             :     !$acc end parallel loop
         332             :
         333             :     !$acc exit data delete( brunt_vaisala_freq_sqd, brunt_vaisala_freq_sqd_mixed, &
         334             :     !$acc                   brunt_vaisala_freq_sqd_moist, brunt_vaisala_freq_sqd_dry, &
         335             :     !$acc                   lambda0_stability )
         336             :
         337     8935056 :     return
         338             :
         339             :   end subroutine calc_stability_correction
         340             :
         341             :   !===============================================================================
         342    17870112 :   subroutine calc_brunt_vaisala_freq_sqd(  nz, ngrdcol, gr, thlm, &
         343    17870112 :                                            exner, rtm, rcm, p_in_Pa, thvm, &
         344    17870112 :                                            ice_supersat_frac, &
         345             :                                            saturation_formula, &
         346             :                                            l_brunt_vaisala_freq_moist, &
         347             :                                            l_use_thvm_in_bv_freq, &
         348    17870112 :                                            bv_efold, &
         349    17870112 :                                            brunt_vaisala_freq_sqd, &
         350    17870112 :                                            brunt_vaisala_freq_sqd_mixed,&
         351    17870112 :                                            brunt_vaisala_freq_sqd_dry, &
         352    17870112 :                                            brunt_vaisala_freq_sqd_moist )
         353             :
         354             :   ! Description:
         355             :   !   Calculate the Brunt-Vaisala frequency squared, N^2.
         356             :
         357             :   ! References:
         358             :   !   ?
         359             :   !-----------------------------------------------------------------------
         360             :
         361             :     use clubb_precision, only: &
         362             :         core_rknd ! Konstant
         363             :
         364             :     use constants_clubb, only: &
         365             :         grav, & ! Constant(s)
         366             :         Lv, &
         367             :         Cp, &
         368             :         Rd, &
         369             :         ep, &
         370             :         one, &
         371             :         one_half, &
         372             :         zero_threshold
         373             :
         374             :     use parameters_model, only: &
         375             :         T0 ! Variable!
         376             :
         377             :     use grid_class, only: &
         378             :         grid, & ! Type
         379             :         ddzt,   &  ! Procedure(s)
         380             :         zt2zm
         381             :
         382             :     use T_in_K_module, only: &
         383             :         thlm2T_in_K ! Procedure
         384             :
         385             :     use saturation, only: &
         386             :         sat_mixrat_liq ! Procedure
         387             :
         388             :     implicit none
         389             :
         390             :     !---------------------------- Input Variables ----------------------------
         391             :     integer, intent(in) :: &
         392             :       nz, &
         393             :       ngrdcol
         394             :
         395             :     type (grid), target, intent(in) :: gr
         396             :
         397             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
         398             :       thlm,    &  ! th_l (thermo. levels)              [K]
         399             :       exner,   &  ! Exner function                     [-]
         400             :       rtm,     &  ! total water mixing ratio, r_t      [kg/kg]
         401             :       rcm,     &  ! cloud water mixing ratio, r_c      [kg/kg]
         402             :       p_in_Pa, &  ! Air pressure                       [Pa]
         403             :       thvm,    &  ! Virtual potential temperature      [K]
         404             :       ice_supersat_frac
         405             :
         406             :     integer, intent(in) :: &
         407             :       saturation_formula ! Integer that stores the saturation formula to be used
         408             :
         409             :     logical, intent(in) :: &
         410             :       l_brunt_vaisala_freq_moist, & ! Use a different formula for the Brunt-Vaisala frequency in
         411             :                                     ! saturated atmospheres (from Durran and Klemp, 1982)
         412             :       l_use_thvm_in_bv_freq         ! Use thvm in the calculation of Brunt-Vaisala frequency
         413             :
         414             :     real( kind = core_rknd ), dimension(ngrdcol), intent(in) :: &
         415             :       bv_efold                      ! Control parameter for inverse e-folding of
         416             :                                     ! cloud fraction in the mixed Brunt Vaisala frequency
         417             :
         418             :     !---------------------------- Output Variables ----------------------------
         419             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
         420             :       brunt_vaisala_freq_sqd, & ! Brunt-Vaisala frequency squared, N^2 [1/s^2]
         421             :       brunt_vaisala_freq_sqd_mixed, &
         422             :       brunt_vaisala_freq_sqd_dry,&
         423             :       brunt_vaisala_freq_sqd_moist
         424             :
         425             :     !---------------------------- Local Variables ----------------------------
         426             :     real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
         427    35740224 :       T_in_K, T_in_K_zm, rsat, rsat_zm, thm, thm_zm, ddzt_thlm, &
         428    35740224 :       ddzt_thm, ddzt_rsat, ddzt_rtm, thvm_zm, ddzt_thvm, ice_supersat_frac_zm
         429             :
         430             :     real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
         431    35740224 :       stat_dry, stat_liq, ddzt_stat_liq, ddzt_stat_liq_zm, &
         432    35740224 :       stat_dry_virtual, stat_dry_virtual_zm,  ddzt_rtm_zm
         433             :
         434             :     integer :: i, k
         435             :
         436             :     !---------------------------- Begin Code ----------------------------
         437             :
         438             :     !$acc data copyin( gr, gr%zt, &
         439             :     !$acc              thlm, exner, rtm, rcm, p_in_Pa, thvm, ice_supersat_frac, bv_efold ) &
         440             :     !$acc      copyout( brunt_vaisala_freq_sqd, brunt_vaisala_freq_sqd_mixed, &
         441             :     !$acc               brunt_vaisala_freq_sqd_dry, brunt_vaisala_freq_sqd_moist ) &
         442             :     !$acc       create( T_in_K, T_in_K_zm, rsat, rsat_zm, thm, thm_zm, ddzt_thlm, &
         443             :     !$acc               ddzt_thm, ddzt_rsat, ddzt_rtm, thvm_zm, ddzt_thvm, stat_dry, &
         444             :     !$acc               stat_liq, ddzt_stat_liq, ddzt_stat_liq_zm, stat_dry_virtual, &
         445             :     !$acc               stat_dry_virtual_zm, ddzt_rtm_zm, ice_supersat_frac_zm )
         446             :
         447    17870112 :     ddzt_thlm = ddzt( nz, ngrdcol, gr, thlm )
         448    17870112 :     thvm_zm   = zt2zm( nz, ngrdcol, gr, thvm, zero_threshold )
         449    17870112 :     ddzt_thvm = ddzt( nz, ngrdcol, gr, thvm )
         450             :
         451             :     ! Dry Brunt-Vaisala frequency
         452    17870112 :     if ( l_use_thvm_in_bv_freq ) then
         453             :
         454             :       !$acc parallel loop gang vector collapse(2) default(present)
         455           0 :       do k = 1, nz
         456           0 :         do i = 1, ngrdcol
         457           0 :           brunt_vaisala_freq_sqd(i,k) = ( grav / thvm_zm(i,k) ) * ddzt_thvm(i,k)
         458             :         end do
         459             :       end do
         460             :       !$acc end parallel loop
         461             :
         462             :     else
         463             :
         464             :       !$acc parallel loop gang vector collapse(2) default(present)
         465  1536829632 :       do k = 1, nz
         466 25380961632 :         do i = 1, ngrdcol
         467 25363091520 :           brunt_vaisala_freq_sqd(i,k) = ( grav / T0 ) * ddzt_thlm(i,k)
         468             :         end do
         469             :       end do
         470             :       !$acc end parallel loop
         471             :
         472             :     end if
         473             :
         474    17870112 :     T_in_K    = thlm2T_in_K( nz, ngrdcol, thlm, exner, rcm )
         475    17870112 :     T_in_K_zm = zt2zm( nz, ngrdcol, gr, T_in_K, zero_threshold )
         476    17870112 :     rsat      = sat_mixrat_liq( nz, ngrdcol, p_in_Pa, T_in_K, saturation_formula )
         477    17870112 :     rsat_zm   = zt2zm( nz, ngrdcol, gr, rsat, zero_threshold )
         478    17870112 :     ddzt_rsat = ddzt( nz, ngrdcol, gr, rsat )
         479             :
         480             :     !$acc parallel loop gang vector collapse(2) default(present)
         481  1536829632 :     do k = 1, nz
         482 25380961632 :       do i = 1, ngrdcol
         483 25363091520 :         thm(i,k) = thlm(i,k) + Lv/(Cp*exner(i,k)) * rcm(i,k)
         484             :       end do
         485             :     end do
         486             :     !$acc end parallel loop
         487             :
         488    17870112 :     thm_zm   = zt2zm( nz, ngrdcol, gr, thm, zero_threshold )
         489    17870112 :     ddzt_thm = ddzt( nz, ngrdcol, gr, thm )
         490    17870112 :     ddzt_rtm = ddzt( nz, ngrdcol, gr, rtm )
         491             :
         492             :     !$acc parallel loop gang vector collapse(2) default(present)
         493  1536829632 :     do k = 1, nz
         494 25380961632 :       do i = 1, ngrdcol
         495 23844132000 :         stat_dry(i,k)  =  Cp * T_in_K(i,k) + grav * gr%zt(i,k)
         496 25363091520 :         stat_liq(i,k)  =  stat_dry(i,k) - Lv * rcm(i,k)
         497             :       end do
         498             :     end do
         499             :     !$acc end parallel loop
         500             :
         501    17870112 :     ddzt_stat_liq    = ddzt( nz, ngrdcol, gr, stat_liq )
         502    17870112 :     ddzt_stat_liq_zm = zt2zm( nz, ngrdcol, gr, ddzt_stat_liq)
         503             :
         504             :     !$acc parallel loop gang vector collapse(2) default(present)
         505  1536829632 :     do k = 1, nz
         506 25380961632 :       do i = 1, ngrdcol
         507 47688264000 :         stat_dry_virtual(i,k) = stat_dry(i,k) + Cp * T_in_K(i,k) &
         508 73051355520 :                                                 *( 0.608 * ( rtm(i,k) - rcm(i,k) )- rcm(i,k) )
         509             :       end do
         510             :     end do
         511             :     !$acc end parallel loop
         512             :
         513    17870112 :     stat_dry_virtual_zm = zt2zm( nz, ngrdcol, gr, stat_dry_virtual, zero_threshold )
         514    17870112 :     ddzt_rtm_zm         = zt2zm( nz, ngrdcol, gr, ddzt_rtm )
         515             :
         516             :     !$acc parallel loop gang vector collapse(2) default(present)
         517  1536829632 :     do k = 1, nz
         518 25380961632 :       do i = 1, ngrdcol
         519 25363091520 :         brunt_vaisala_freq_sqd_dry(i,k) = ( grav / thm_zm(i,k) )* ddzt_thm(i,k)
         520             :       end do
         521             :     end do
         522             :     !$acc end parallel loop
         523             :
         524             :     !$acc parallel loop gang vector collapse(2) default(present)
         525  1536829632 :     do k = 1, nz
         526 25380961632 :       do i = 1, ngrdcol
         527             :         ! In-cloud Brunt-Vaisala frequency. This is Eq. (36) of Durran and
         528             :         ! Klemp (1982)
         529 47688264000 :         brunt_vaisala_freq_sqd_moist(i,k) = &
         530             :           grav * ( ((one + Lv*rsat_zm(i,k) / (Rd*T_in_K_zm(i,k))) / &
         531             :           (one + ep*(Lv**2)*rsat_zm(i,k)/(Cp*Rd*T_in_K_zm(i,k)**2))) * &
         532             :           ( (one/thm_zm(i,k) * ddzt_thm(i,k)) + (Lv/(Cp*T_in_K_zm(i,k)))*ddzt_rsat(i,k)) - &
         533 73051355520 :           ddzt_rtm(i,k) )
         534             :       end do
         535             :     end do ! k=1, gr%nz
         536             :     !$acc end parallel loop
         537             :
         538    17870112 :     ice_supersat_frac_zm = zt2zm( nz, ngrdcol, gr, ice_supersat_frac, zero_threshold )
         539             :
         540             :     !$acc parallel loop gang vector collapse(2) default(present)
         541  1536829632 :     do k = 1, nz
         542 25380961632 :       do i = 1, ngrdcol
         543 47688264000 :          brunt_vaisala_freq_sqd_mixed(i,k) = &
         544             :              brunt_vaisala_freq_sqd_moist(i,k) + &
         545             :                  exp( - bv_efold(i) * ice_supersat_frac_zm(i,k) ) * &
         546 73051355520 :                  ( brunt_vaisala_freq_sqd_dry(i,k) - brunt_vaisala_freq_sqd_moist(i,k) )
         547             :       end do
         548             :     end do
         549             :     !$acc end parallel loop
         550             :
         551    17870112 :     if ( l_brunt_vaisala_freq_moist ) then
         552             :
         553           0 :       brunt_vaisala_freq_sqd = brunt_vaisala_freq_sqd_moist
         554             :
         555             :     end if
         556             :
         557             :     !$acc end data
         558             :
         559    17870112 :     return
         560             :
         561             :   end subroutine calc_brunt_vaisala_freq_sqd
         562             :
         563             : !===============================================================================
         564           0 :   subroutine compute_Cx_fnc_Richardson( nz, ngrdcol, gr, &
         565           0 :                                         thlm, um, vm, em, Lscale, exner, rtm, &
         566           0 :                                         rcm, p_in_Pa, thvm, rho_ds_zm, &
         567           0 :                                         ice_supersat_frac, &
         568           0 :                                         clubb_params, &
         569             :                                         saturation_formula, &
         570             :                                         l_brunt_vaisala_freq_moist, &
         571             :                                         l_use_thvm_in_bv_freq, &
         572             :                                         l_use_shear_Richardson, &
         573             :                                         l_modify_limiters_for_cnvg_test, &
         574             :                                         stats_metadata, &
         575           0 :                                         stats_zm, &
         576           0 :                                         Cx_fnc_Richardson )
         577             :
         578             :   ! Description:
         579             :   !   Compute Cx as a function of the Richardson number
         580             :
         581             :   ! References:
         582             :   !   cam:ticket:59
         583             :   !-----------------------------------------------------------------------
         584             :
         585             :     use clubb_precision, only: &
         586             :         core_rknd  ! Konstant
         587             :
         588             :     use grid_class, only: &
         589             :         grid, & ! Type
         590             :         ddzt, & ! Procedure(s)
         591             :         zt2zm, &
         592             :         zm2zt2zm
         593             :
         594             :     use constants_clubb, only: &
         595             :         one, &
         596             :         zero, &
         597             :         zero_threshold
         598             :
         599             :     use interpolation, only: &
         600             :         linear_interp_factor ! Procedure
         601             :
         602             :     use parameter_indices, only: &
         603             :         nparams,             & ! Variable(s)
         604             :         iCx_min,             &
         605             :         iCx_max,             &
         606             :         iRichardson_num_min, &
         607             :         iRichardson_num_max, &
         608             :         ibv_efold
         609             :
         610             :     use stats_variables, only: &
         611             :         stats_metadata_type
         612             :
         613             :     use stats_type_utilities, only: &
         614             :         stat_update_var      ! Procedure
         615             :
         616             :     use stats_type, only: stats ! Type
         617             :
         618             :     implicit none
         619             :
         620             :     !------------------------------ Constant Parameters ------------------------------
         621             :     real( kind = core_rknd ), parameter :: &
         622             :       Richardson_num_divisor_threshold = 1.0e-6_core_rknd, &
         623             :       Cx_fnc_Richardson_below_ground_value = one
         624             :
         625             :     logical, parameter :: &
         626             :       l_Cx_fnc_Richardson_vert_avg = .false.    ! Vertically average Cx_fnc_Richardson over a
         627             :                                                 !  distance of Lscale
         628             :
         629             :     real( kind = core_rknd ), parameter :: &
         630             :       min_max_smth_mag = 1.0e-9_core_rknd ! "base" smoothing magnitude before scaling
         631             :                                           ! for the respective data structure. See
         632             :                                           ! https://github.com/larson-group/clubb/issues/965#issuecomment-1119816722
         633             :                                           ! for a plot on how output behaves with varying min_max_smth_mag
         634             :
         635             :     !------------------------------ Input Variables ------------------------------
         636             :     integer, intent(in) :: &
         637             :       nz, &
         638             :       ngrdcol
         639             :
         640             :     type (grid), target, intent(in) :: gr
         641             :
         642             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
         643             :       thlm,      & ! th_l (liquid water potential temperature)      [K]
         644             :       um,        & ! u mean wind component (thermodynamic levels)   [m/s]
         645             :       vm,        & ! v mean wind component (thermodynamic levels)   [m/s]
         646             :       em,        & ! Turbulent Kinetic Energy (TKE)                 [m^2/s^2]
         647             :       Lscale,    & ! Turbulent mixing length                        [m]
         648             :       exner,     & ! Exner function                                 [-]
         649             :       rtm,       & ! total water mixing ratio, r_t                  [kg/kg]
         650             :       rcm,       & ! cloud water mixing ratio, r_c                  [kg/kg]
         651             :       p_in_Pa,   & ! Air pressure                                   [Pa]
         652             :       thvm,      & ! Virtual potential temperature                  [K]
         653             :       rho_ds_zm, &  ! Dry static density on momentum levels          [kg/m^3]
         654             :       ice_supersat_frac  ! ice cloud fraction
         655             :
         656             :     real( kind = core_rknd ), dimension(ngrdcol,nparams), intent(in) :: &
         657             :       clubb_params    ! Array of CLUBB's tunable parameters    [units vary]
         658             :
         659             :     integer, intent(in) :: &
         660             :       saturation_formula ! Integer that stores the saturation formula to be used
         661             :
         662             :     logical, intent(in) :: &
         663             :       l_brunt_vaisala_freq_moist, & ! Use a different formula for the Brunt-Vaisala frequency in
         664             :                                     ! saturated atmospheres (from Durran and Klemp, 1982)
         665             :       l_use_thvm_in_bv_freq,      & ! Use thvm in the calculation of Brunt-Vaisala frequency
         666             :       l_use_shear_Richardson        ! Use shear in the calculation of Richardson number
         667             :
         668             :     ! Flag to activate modifications on limiters for convergence test
         669             :     ! (smoothed max and min for Cx_fnc_Richardson in advance_helper_module.F90)
         670             :     ! (remove the clippings on brunt_vaisala_freq_sqd_smth in mixing_length.F90)
         671             :     ! (reduce threshold on limiters for sqrt_Ri_zm in mixing_length.F90)
         672             :     logical, intent(in) :: &
         673             :       l_modify_limiters_for_cnvg_test
         674             :
         675             :     type (stats_metadata_type), intent(in) :: &
         676             :       stats_metadata
         677             :
         678             :     !------------------------------ InOut Variable ------------------------------
         679             :     type (stats), target, dimension(ngrdcol), intent(inout) :: &
         680             :       stats_zm
         681             :
         682             :     !------------------------------ Output Variable ------------------------------
         683             :     real( kind = core_rknd), dimension(ngrdcol,nz), intent(out) :: &
         684             :       Cx_fnc_Richardson
         685             :
         686             :     !------------------------------ Local Variables ------------------------------
         687             :     real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
         688           0 :       brunt_vaisala_freq_sqd, &
         689           0 :       brunt_vaisala_freq_sqd_mixed,&
         690           0 :       brunt_vaisala_freq_sqd_dry, &
         691           0 :       brunt_vaisala_freq_sqd_moist, &
         692           0 :       fnc_Richardson, &
         693           0 :       fnc_Richardson_clipped, &
         694           0 :       fnc_Richardson_smooth, &
         695           0 :       Ri_zm, &
         696           0 :       ddzt_um, &
         697           0 :       ddzt_vm, &
         698           0 :       shear_sqd, &
         699           0 :       Lscale_zm, &
         700           0 :       Cx_fnc_interp, &
         701           0 :       Cx_fnc_Richardson_avg
         702             :
         703             :     real ( kind = core_rknd ) :: &
         704             :       invrs_min_max_diff, &
         705             :       invrs_num_div_thresh
         706             :
         707             :     real( kind = core_rknd ) :: &
         708             :       Richardson_num_max, & ! CLUBB tunable parameter Richardson_num_max
         709             :       Richardson_num_min, & ! CLUBB tunable parameter Richardson_num_min
         710             :       Cx_max,             & ! CLUBB tunable parameter max of Cx_fnc_Richardson
         711             :       Cx_min                ! CLUBB tunable parameter min of Cx_fnc_Richardson
         712             :
         713             :     integer :: smth_type = 1
         714             :
         715             :     integer :: i, k
         716             :
         717             :     !------------------------------ Begin Code ------------------------------
         718             :
         719             :     !$acc enter data create( brunt_vaisala_freq_sqd, brunt_vaisala_freq_sqd_mixed, &
         720             :     !$acc                    brunt_vaisala_freq_sqd_dry, brunt_vaisala_freq_sqd_moist, &
         721             :     !$acc                    Cx_fnc_interp, &
         722             :     !$acc                    Ri_zm, ddzt_um, ddzt_vm, shear_sqd, Lscale_zm, &
         723             :     !$acc                    Cx_fnc_Richardson_avg, fnc_Richardson, &
         724             :     !$acc                    fnc_Richardson_clipped, fnc_Richardson_smooth )
         725             :
         726             :     call calc_brunt_vaisala_freq_sqd( nz, ngrdcol, gr, thlm, &          ! intent(in)
         727             :                                       exner, rtm, rcm, p_in_Pa, thvm, & ! intent(in)
         728             :                                       ice_supersat_frac, &              ! intent(in)
         729             :                                       saturation_formula, &             ! intent(in)
         730             :                                       l_brunt_vaisala_freq_moist, &     ! intent(in)
         731             :                                       l_use_thvm_in_bv_freq, &          ! intent(in)
         732             :                                       clubb_params(:,ibv_efold), &       ! intent(in)
         733             :                                       brunt_vaisala_freq_sqd, &         ! intent(out)
         734             :                                       brunt_vaisala_freq_sqd_mixed,&    ! intent(out)
         735             :                                       brunt_vaisala_freq_sqd_dry, &     ! intent(out)
         736           0 :                                       brunt_vaisala_freq_sqd_moist )    ! intent(out)
         737             :
         738           0 :     invrs_num_div_thresh = one / Richardson_num_divisor_threshold
         739             :
         740           0 :     Lscale_zm = zt2zm( nz, ngrdcol, gr, Lscale, zero_threshold )
         741             :
         742             :     ! Calculate shear_sqd
         743           0 :     ddzt_um = ddzt( nz, ngrdcol, gr, um )
         744           0 :     ddzt_vm = ddzt( nz, ngrdcol, gr, vm )
         745             :
         746             :     !$acc parallel loop gang vector collapse(2) default(present)
         747           0 :     do k = 1, nz
         748           0 :       do i = 1, ngrdcol
         749           0 :         shear_sqd(i,k) = ddzt_um(i,k)**2 + ddzt_vm(i,k)**2
         750             :       end do
         751             :     end do
         752             :     !$acc end parallel loop
         753             :
         754           0 :     if ( stats_metadata%l_stats_samp ) then
         755             :       !$acc update host(shear_sqd)
         756           0 :       do i = 1, ngrdcol
         757           0 :         call stat_update_var( stats_metadata%ishear_sqd, shear_sqd(i,:), & ! intent(in)
         758           0 :                               stats_zm(i) )               ! intent(inout)
         759             :       end do
         760             :     end if
         761             :
         762           0 :     if ( l_use_shear_Richardson ) then
         763             :
         764             :       !$acc parallel loop gang vector collapse(2) default(present)
         765           0 :       do k = 1, nz
         766           0 :         do i = 1, ngrdcol
         767           0 :           Ri_zm(i,k) = max( 1.0e-7_core_rknd, brunt_vaisala_freq_sqd_mixed(i,k) ) &
         768           0 :                        / max( shear_sqd(i,k), 1.0e-7_core_rknd )
         769             :         end do
         770             :       end do
         771             :       !$acc end parallel loop
         772             :
         773             :     else
         774             :       !$acc parallel loop gang vector collapse(2) default(present)
         775           0 :       do k = 1, nz
         776           0 :         do i = 1, ngrdcol
         777           0 :           Ri_zm(i,k) = brunt_vaisala_freq_sqd(i,k) * invrs_num_div_thresh
         778             :         end do
         779             :       end do
         780             :       !$acc end parallel loop
         781             :     end if
         782             :
         783             :     ! Cx_fnc_Richardson is interpolated based on the value of Richardson_num
         784             :     ! The min function ensures that Cx does not exceed Cx_max, regardless of the
         785             :     !     value of Richardson_num_max.
         786           0 :     if ( l_modify_limiters_for_cnvg_test ) then
         787             :
         788             :       !$acc parallel loop gang vector collapse(2) default(present)
         789           0 :       do k = 1, nz
         790           0 :         do i = 1, ngrdcol
         791             :
         792           0 :           Richardson_num_max = clubb_params(i,iRichardson_num_max)
         793           0 :           Richardson_num_min = clubb_params(i,iRichardson_num_min)
         794             :
         795           0 :           invrs_min_max_diff = one / ( Richardson_num_max - Richardson_num_min )
         796             :
         797           0 :           fnc_Richardson(i,k) = ( Ri_zm(i,k) - clubb_params(i,iRichardson_num_min) ) * invrs_min_max_diff
         798             :         end do
         799             :       end do
         800             :
         801             :       fnc_Richardson_clipped = smooth_min( nz, ngrdcol, one, &
         802             :                                            fnc_Richardson, &
         803           0 :                                            min_max_smth_mag )
         804             :
         805             :       fnc_Richardson_smooth = smooth_max( nz, ngrdcol, zero, &
         806             :                                           fnc_Richardson_clipped, &
         807           0 :                                           min_max_smth_mag )
         808             :
         809             :       ! use smoothed max amd min to achive smoothed profile and avoid discontinuities
         810             :       !$acc parallel loop gang vector collapse(2) default(present)
         811           0 :       do k = 1, nz
         812           0 :         do i = 1, ngrdcol
         813             :
         814           0 :           Cx_max = clubb_params(i,iCx_max)
         815           0 :           Cx_min = clubb_params(i,iCx_min)
         816             :
         817           0 :           Cx_fnc_interp(i,k) = fnc_Richardson_smooth(i,k) * ( Cx_max - Cx_min ) + Cx_min
         818             :         end do
         819             :       end do
         820             :
         821           0 :       Cx_fnc_Richardson = zm2zt2zm( nz, ngrdcol, gr, Cx_fnc_interp )
         822             :
         823             :     else ! default method
         824             :
         825             :       !$acc parallel loop gang vector collapse(2) default(present)
         826           0 :       do k = 1, nz
         827           0 :         do i = 1, ngrdcol
         828             :
         829           0 :           Richardson_num_max = clubb_params(i,iRichardson_num_max)
         830           0 :           Richardson_num_min = clubb_params(i,iRichardson_num_min)
         831           0 :           Cx_max = clubb_params(i,iCx_max)
         832           0 :           Cx_min = clubb_params(i,iCx_min)
         833             :
         834           0 :           invrs_min_max_diff = one / ( Richardson_num_max - Richardson_num_min )
         835             :
         836           0 :           Cx_fnc_Richardson(i,k) = ( max(min(Richardson_num_max, Ri_zm(i,k)), Richardson_num_min) &
         837             :                                      - Richardson_num_min )  &
         838           0 :                                    * invrs_min_max_diff * ( Cx_max - Cx_min ) + Cx_min
         839             :         end do
         840             :       end do
         841             :       !$acc end parallel loop
         842             :
         843             :     end if
         844             :
         845             :     if ( l_Cx_fnc_Richardson_vert_avg ) then
         846             :       Cx_fnc_Richardson = Lscale_width_vert_avg( nz, ngrdcol, gr, smth_type, &
         847             :                                                  Cx_fnc_Richardson, Lscale_zm, rho_ds_zm, &
         848             :                                                  Cx_fnc_Richardson_below_ground_value )
         849             :
         850             :       !$acc parallel loop gang vector collapse(2) default(present)
         851             :       do k = 1, nz
         852             :         do i = 1, ngrdcol
         853             :           Cx_fnc_Richardson(i,k) = Cx_fnc_Richardson_avg(i,k)
         854             :         end do
         855             :       end do
         856             :       !$acc end parallel loop
         857             :     end if
         858             :
         859             :     ! On some compilers, roundoff error can result in Cx_fnc_Richardson being
         860             :     ! slightly outside the range [0,1]. Thus, it is clipped here.
         861             :     !$acc parallel loop gang vector collapse(2) default(present)
         862           0 :     do k = 1, nz
         863           0 :       do i = 1, ngrdcol
         864           0 :         Cx_fnc_Richardson(i,k) = max( zero, min( one, Cx_fnc_Richardson(i,k) ) )
         865             :       end do
         866             :     end do
         867             :     !$acc end parallel loop
         868             :
         869             :     !$acc exit data delete( brunt_vaisala_freq_sqd, brunt_vaisala_freq_sqd_mixed, &
         870             :     !$acc                   brunt_vaisala_freq_sqd_dry, brunt_vaisala_freq_sqd_moist, &
         871             :     !$acc                   Cx_fnc_interp, Ri_zm, &
         872             :     !$acc                   ddzt_um, ddzt_vm, shear_sqd, Lscale_zm, &
         873             :     !$acc                   Cx_fnc_Richardson_avg, fnc_Richardson, &
         874             :     !$acc                   fnc_Richardson_clipped, fnc_Richardson_smooth )
         875             :
         876           0 :     return
         877             :
         878             :   end subroutine compute_Cx_fnc_Richardson
         879             :   !----------------------------------------------------------------------
         880             :
         881             :   !----------------------------------------------------------------------
         882     8935056 :   function Lscale_width_vert_avg( nz, ngrdcol, gr, smth_type, &
         883     8935056 :                                   var_profile, Lscale_zm, rho_ds_zm, &
         884             :                                   var_below_ground_value )&
         885    35740224 :   result (Lscale_width_vert_avg_output)
         886             :
         887             :   ! Description:
         888             :   !   Averages a profile with a running mean of width Lscale_zm
         889             :
         890             :   ! References:
         891             :   !   cam:ticket:59
         892             :
         893             :     use clubb_precision, only: &
         894             :         core_rknd ! Precision
         895             :
         896             :     use grid_class, only: &
         897             :         grid ! Type
         898             :
         899             :     use constants_clubb, only: &
         900             :         zero
         901             :
         902             :     implicit none
         903             :
         904             :     !-------------------------- Input Variables --------------------------
         905             :     integer, intent(in) :: &
         906             :       nz, &
         907             :       ngrdcol, &
         908             :       smth_type
         909             :
         910             :     type (grid), target, intent(in) :: gr
         911             :
         912             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
         913             :       var_profile, &      ! Profile on momentum levels
         914             :       Lscale_zm, &        ! Lscale on momentum levels
         915             :       rho_ds_zm           ! Dry static energy on momentum levels!
         916             :
         917             :     real( kind = core_rknd ), intent(in) :: &
         918             :       var_below_ground_value ! Value to use below ground
         919             :
         920             :     !-------------------------- Result Variable --------------------------
         921             :     real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
         922             :       Lscale_width_vert_avg_output ! Vertically averaged profile (on momentum levels)
         923             :
         924             :     !-------------------------- Local Variables --------------------------
         925             :     integer :: &
         926             :         k, i,        & ! Loop variable
         927             :         k_avg_lower, &
         928             :         k_avg_upper, &
         929             :         k_avg
         930             :
         931             :     real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
         932    17870112 :       one_half_avg_width, &
         933    17870112 :       numer_terms, &
         934     8935056 :       denom_terms
         935             :
         936             :     integer :: &
         937             :       n_below_ground_levels
         938             :
         939             :     real( kind = core_rknd ) :: &
         940             :       numer_integral, & ! Integral in the numerator (see description)
         941             :       denom_integral    ! Integral in the denominator (see description)
         942             :
         943             :     !-------------------------- Begin Code --------------------------
         944             :
         945             :     !$acc enter data create( one_half_avg_width, numer_terms, denom_terms )
         946             :
         947     8935056 :     if ( smth_type == 1 ) then
         948             :       !$acc parallel loop gang vector collapse(2) default(present)
         949           0 :       do k = 1, nz
         950           0 :         do i = 1, ngrdcol
         951           0 :           one_half_avg_width(i,k) = max( Lscale_zm(i,k), 500.0_core_rknd )
         952             :         end do
         953             :       end do
         954     8935056 :     else if (smth_type == 2 ) then
         955             :       !$acc parallel loop gang vector collapse(2) default(present)
         956   768414816 :       do k = 1, nz
         957 12690480816 :         do i = 1, ngrdcol
         958 12681545760 :           one_half_avg_width(i,k) = 60.0_core_rknd
         959             :         end do
         960             :       end do
         961             :     endif
         962             :
         963             :     ! Pre calculate numerator and denominator terms
         964             :     !$acc parallel loop gang vector collapse(2) default(present)
         965   768414816 :     do k = 1, nz
         966 12690480816 :       do i = 1, ngrdcol
         967 11922066000 :         numer_terms(i,k) = rho_ds_zm(i,k) * gr%dzm(i,k) * var_profile(i,k)
         968 12681545760 :         denom_terms(i,k) = rho_ds_zm(i,k) * gr%dzm(i,k)
         969             :       end do
         970             :     end do
         971             :
         972             :     ! For every grid level
         973             :     !$acc parallel loop gang vector collapse(2) default(present)
         974   768414816 :     do k = 1, nz
         975 12690480816 :       do i = 1, ngrdcol
         976             :
         977             :         !-----------------------------------------------------------------------
         978             :         ! Hunt down all vertical levels with one_half_avg_width(k) of gr%zm(k).
         979             :         !
         980             :         ! Note: Outdated explanation of version that improves CPU performance
         981             :         !       below. Reworked due to it requiring a k dependency. Now we
         982             :         !       begin looking for k_avg_upper and k_avg_lower starting at
         983             :         !       the kth level.
         984             :         !
         985             :         ! Outdated but potentially useful note:
         986             :         !     k_avg_upper and k_avg_lower can be saved each loop iteration, this
         987             :         !     reduces iterations beacuse the kth values are likely to be within
         988             :         !     one or two grid levels of the k-1th values. Less searching is required
         989             :         !     by starting the search at the previous values and incrementing or
         990             :         !     decrement as needed.
         991             :         !-----------------------------------------------------------------------
         992             :
         993             :         ! Determine if k_avg_upper needs to increment
         994 12081472528 :         do k_avg_upper = k, nz-1
         995 12081472528 :           if ( gr%zm(i,k_avg_upper+1) - gr%zm(i,k) > one_half_avg_width(i,k) ) then
         996             :             exit
         997             :           end if
         998             :         end do
         999             :
        1000             :         ! Determine if k_avg_lower needs to decrement
        1001 12081472528 :         do k_avg_lower = k, 2, -1
        1002 12081472528 :           if ( gr%zm(i,k) - gr%zm(i,k_avg_lower-1) > one_half_avg_width(i,k) ) then
        1003             :             exit
        1004             :           end if
        1005             :         end do
        1006             :
        1007             :         ! Compute the number of levels below ground to include.
        1008 11922066000 :         if ( k_avg_lower > 1 ) then
        1009             :
        1010             :           ! k=1, the lowest "real" level, is not included in the average, so no
        1011             :           ! below-ground levels should be included.
        1012 24162945056 :           n_below_ground_levels = 0
        1013             :
        1014             :           numer_integral = zero
        1015             :           denom_integral = zero
        1016             :
        1017             :         else
        1018             :
        1019             :           ! The number of below-ground levels included is equal to the distance
        1020             :           ! below the lowest level spanned by one_half_avg_width(k)
        1021             :           ! divided by the distance between vertical levels below ground; the
        1022             :           ! latter is assumed to be the same as the distance between the first and
        1023             :           ! second vertical levels.
        1024   841557600 :           n_below_ground_levels = int( ( one_half_avg_width(i,k)-(gr%zm(i,k)-gr%zm(i,1)) ) / &
        1025  1122076800 :                                       ( gr%zm(i,2)-gr%zm(i,1) ) )
        1026             :
        1027   280519200 :           numer_integral = n_below_ground_levels * denom_terms(i,1) * var_below_ground_value
        1028             :           denom_integral = n_below_ground_levels * denom_terms(i,1)
        1029             :
        1030             :         end if
        1031             :
        1032             :         ! Add numerator and denominator terms for all above-ground levels
        1033 24162945056 :         do k_avg = k_avg_lower, k_avg_upper
        1034             :
        1035 12240879056 :           numer_integral = numer_integral + numer_terms(i,k_avg)
        1036 24162945056 :           denom_integral = denom_integral + denom_terms(i,k_avg)
        1037             :
        1038             :         end do
        1039             :
        1040 12681545760 :         Lscale_width_vert_avg_output(i,k) = numer_integral / denom_integral
        1041             :
        1042             :       end do
        1043             :     end do
        1044             :
        1045             :     !$acc exit data delete( one_half_avg_width, numer_terms, denom_terms )
        1046             :
        1047     8935056 :     return
        1048             :
        1049     8935056 :   end function Lscale_width_vert_avg
        1050             :
        1051             :  !=============================================================================
        1052     8935056 :   subroutine wp2_term_splat_lhs( nz, ngrdcol, gr, C_wp2_splat, &
        1053     8935056 :                                  brunt_vaisala_freq_sqd_splat, &
        1054     8935056 :                                  lhs_splat_wp2 )
        1055             :
        1056             :     ! Description
        1057             :     ! DESCRIBE TERM
        1058             :
        1059             :     ! References:
        1060             :     !-----------------------------------------------------------------------
        1061             :
        1062             :     use grid_class, only:  &
        1063             :         grid, & ! Type
        1064             :         zm2zt2zm
        1065             :
        1066             :     use constants_clubb, only: &
        1067             :         zero
        1068             :
        1069             :     use clubb_precision, only: &
        1070             :         core_rknd    ! Variable(s)
        1071             :
        1072             :     implicit none
        1073             :
        1074             :     ! --------------------- Input Variables ---------------------
        1075             :     integer, intent(in) :: &
        1076             :       nz, &
        1077             :       ngrdcol
        1078             :
        1079             :     type (grid), target, intent(in) :: &
        1080             :       gr
        1081             :
        1082             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
        1083             :       brunt_vaisala_freq_sqd_splat  ! Inverse time-scale tau at momentum levels  [1/s^2]
        1084             :
        1085             :     real( kind = core_rknd ), dimension(ngrdcol), intent(in) :: &
        1086             :       C_wp2_splat    ! Model parameter C_wp2_splat             [ -]
        1087             :
        1088             :     ! --------------------- Output Variable ---------------------
        1089             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
        1090             :       lhs_splat_wp2    ! LHS coefficient of wp2 splatting term  [1/s]
        1091             :
        1092             :     ! --------------------- Local Variables ---------------------
        1093             :     real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
        1094    17870112 :       brunt_vaisala_freq_splat_clipped, &
        1095    17870112 :       brunt_vaisala_freq_splat_smooth
        1096             :
        1097             :     integer :: i, k
        1098             :
        1099             :     !----------------------------- Begin Code -----------------------------
        1100             :
        1101             :     !$acc enter data create( brunt_vaisala_freq_splat_clipped, brunt_vaisala_freq_splat_smooth )
        1102             :
        1103             :     !$acc parallel loop gang vector collapse(2) default(present)
        1104   768414816 :     do k = 1, nz
        1105 12690480816 :       do i = 1, ngrdcol
        1106 23844132000 :         brunt_vaisala_freq_splat_clipped(i,k) &
        1107 36525677760 :                 = sqrt( max( zero, brunt_vaisala_freq_sqd_splat(i,k) ) )
        1108             :       end do
        1109             :     end do
        1110             :     !$acc end parallel loop
        1111             :
        1112             :     brunt_vaisala_freq_splat_smooth = zm2zt2zm( nz, ngrdcol, gr, &
        1113     8935056 :                                                 brunt_vaisala_freq_splat_clipped )
        1114             :
        1115             :     !$acc parallel loop gang vector collapse(2) default(present)
        1116   768414816 :     do k = 1, nz
        1117 12690480816 :       do i = 1, ngrdcol
        1118 12681545760 :         lhs_splat_wp2(i,k) = + C_wp2_splat(i) * brunt_vaisala_freq_splat_smooth(i,k)
        1119             :       end do
        1120             :     end do
        1121             :     !$acc end parallel loop
        1122             :
        1123             :     !$acc exit data delete( brunt_vaisala_freq_splat_clipped, brunt_vaisala_freq_splat_smooth )
        1124             :
        1125     8935056 :     return
        1126             :
        1127             :   end subroutine wp2_term_splat_lhs
        1128             :
        1129             :  !=============================================================================
        1130     8935056 :   subroutine wp3_term_splat_lhs( nz, ngrdcol, gr, C_wp2_splat, &
        1131     8935056 :                                  brunt_vaisala_freq_sqd_splat, &
        1132     8935056 :                                  lhs_splat_wp3 )
        1133             :
        1134             :     ! Description
        1135             :     ! DESCRIBE TERM
        1136             :
        1137             :     ! References:
        1138             :     !-----------------------------------------------------------------------
        1139             :
        1140             :     use grid_class, only:  &
        1141             :         grid, & ! Type
        1142             :         zm2zt
        1143             :
        1144             :     use constants_clubb, only: &
        1145             :         zero, &
        1146             :         one_half, &
        1147             :         three
        1148             :
        1149             :     use clubb_precision, only: &
        1150             :         core_rknd    ! Variable(s)
        1151             :
        1152             :     implicit none
        1153             :
        1154             :     ! --------------------- Input Variables ---------------------
        1155             :     integer, intent(in) :: &
        1156             :       nz, &
        1157             :       ngrdcol
        1158             :
        1159             :     type (grid), target, intent(in) :: &
        1160             :       gr
        1161             :
        1162             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
        1163             :       brunt_vaisala_freq_sqd_splat  ! Inverse time-scale tau at momentum levels  [1/s^2]
        1164             :
        1165             :     real( kind = core_rknd ), dimension(ngrdcol), intent(in) :: &
        1166             :       C_wp2_splat    ! Model parameter C_wp2_splat              [-]
        1167             :
        1168             :     ! --------------------- Output Variable ---------------------
        1169             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
        1170             :       lhs_splat_wp3    ! LHS coefficient of wp3 splatting term [1/s]
        1171             :
        1172             :     ! --------------------- Local Variables ---------------------
        1173             :     real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
        1174    17870112 :       brunt_vaisala_freq_splat_clipped, &
        1175    17870112 :       brunt_vaisala_freq_splat_smooth
        1176             :
        1177             :     integer :: i, k
        1178             :
        1179             :     !----------------------------- Begin Code -----------------------------
        1180             :
        1181             :     !$acc enter data create( brunt_vaisala_freq_splat_clipped, brunt_vaisala_freq_splat_smooth )
        1182             :
        1183             :     !$acc parallel loop gang vector collapse(2) default(present)
        1184   768414816 :     do k = 1, nz
        1185 12690480816 :       do i = 1, ngrdcol
        1186 23844132000 :         brunt_vaisala_freq_splat_clipped(i,k) &
        1187 36525677760 :                 = sqrt( max( zero, brunt_vaisala_freq_sqd_splat(i,k) ) )
        1188             :       end do
        1189             :     end do
        1190             :     !$acc end parallel loop
        1191             :
        1192             :     brunt_vaisala_freq_splat_smooth = zm2zt( nz, ngrdcol, gr, &
        1193     8935056 :                                              brunt_vaisala_freq_splat_clipped )
        1194             :
        1195             :     !$acc parallel loop gang vector collapse(2) default(present)
        1196   768414816 :     do k = 1, nz
        1197 12690480816 :       do i = 1, ngrdcol
        1198 23844132000 :         lhs_splat_wp3(i,k) = + one_half * three * C_wp2_splat(i) &
        1199 36525677760 :                                * brunt_vaisala_freq_splat_smooth(i,k)
        1200             :       end do
        1201             :     end do
        1202             :     !$acc end parallel loop
        1203             :
        1204             :     !$acc exit data delete( brunt_vaisala_freq_splat_clipped, brunt_vaisala_freq_splat_smooth )
        1205             :
        1206     8935056 :     return
        1207             :
        1208             :   end subroutine wp3_term_splat_lhs
        1209             :
        1210             : !===============================================================================
        1211           0 :   function smooth_min_sclr_idx( nz, ngrdcol, input_var1, input_var2, smth_coef ) &
        1212           0 :   result( output_var )
        1213             :
        1214             :   ! Description:
        1215             :   !   Computes a smoothed version of the min function, using one scalar and
        1216             :   !   one 1d array as inputs. For more details, see the interface in this file.
        1217             :
        1218             :   ! References:
        1219             :   !   See clubb:ticket:894, updated version: 965
        1220             :   !----------------------------------------------------------------------
        1221             :
        1222             :     use clubb_precision, only: &
        1223             :         core_rknd                     ! Constant(s)
        1224             :
        1225             :     use constants_clubb, only: &
        1226             :         one_half
        1227             :
        1228             :     implicit none
        1229             :
        1230             :     integer, intent(in) :: &
        1231             :       nz, &
        1232             :       ngrdcol
        1233             :
        1234             :     !----------------------------- Input Variables -----------------------------
        1235             :     real ( kind = core_rknd ), intent(in) :: &
        1236             :       input_var1, &       ! Units vary
        1237             :       smth_coef           ! "intensity" of the smoothing. Should be of a similar magnitude to
        1238             :                           ! that of the data structures input_var1 and input_var2
        1239             :
        1240             :     real ( kind = core_rknd ), dimension(ngrdcol, nz), intent(in) :: &
        1241             :       input_var2          ! Units vary
        1242             :
        1243             :     !----------------------------- Output Variables -----------------------------
        1244             :     real( kind = core_rknd ), dimension(ngrdcol, nz) :: &
        1245             :       output_var          ! Same unit as input_var1 and input_var2
        1246             :
        1247             :     !----------------------------- Local Variables -----------------------------
        1248             :     integer :: i, k
        1249             :
        1250             :     !----------------------------- Begin Code -----------------------------
        1251             :
        1252             :     !$acc data copyin( input_var2 ) &
        1253             :     !$acc     copyout( output_var )
        1254             :
        1255             :     !$acc parallel loop gang vector collapse(2) default(present)
        1256           0 :     do k = 1, nz
        1257           0 :       do i = 1, ngrdcol
        1258           0 :         output_var(i,k) = one_half * ( (input_var1+input_var2(i,k)) - &
        1259           0 :                                   sqrt((input_var1-input_var2(i,k))**2 + smth_coef**2) )
        1260             :       end do
        1261             :     end do
        1262             :     !$acc end parallel loop
        1263             :
        1264             :     !$acc end data
        1265             :
        1266           0 :     return
        1267             :
        1268           0 :   end function smooth_min_sclr_idx
        1269             :
        1270             : !===============================================================================
        1271           0 :   function smooth_min_array_scalar( nz, ngrdcol, input_var1, input_var2, smth_coef ) &
        1272           0 :   result( output_var )
        1273             :
        1274             :   ! Description:
        1275             :   !   Computes a smoothed version of the min function, using one scalar and
        1276             :   !   one 1d array as inputs. For more details, see the interface in this file.
        1277             :
        1278             :   ! References:
        1279             :   !   See clubb:ticket:894, updated version: 965
        1280             :   !----------------------------------------------------------------------
        1281             :
        1282             :     use clubb_precision, only: &
        1283             :         core_rknd                     ! Constant(s)
        1284             :
        1285             :     use constants_clubb, only: &
        1286             :         one_half
        1287             :
        1288             :     implicit none
        1289             :
        1290             :     !----------------------------- Input Variables -----------------------------
        1291             :     integer, intent(in) :: &
        1292             :       nz, &
        1293             :       ngrdcol
        1294             :
        1295             :     real ( kind = core_rknd ), dimension(ngrdcol, nz), intent(in) :: &
        1296             :       input_var1          ! Units vary
        1297             :
        1298             :     real ( kind = core_rknd ), intent(in) :: &
        1299             :       input_var2, &       ! Units vary
        1300             :       smth_coef           ! "intensity" of the smoothing. Should be of a similar magnitude to
        1301             :                           ! that of the data structures input_var1 and input_var2
        1302             :
        1303             :     !----------------------------- Output Variables -----------------------------
        1304             :     real( kind = core_rknd ), dimension(ngrdcol, nz) :: &
        1305             :       output_var          ! Same unit as input_var1 and input_var2
        1306             :
        1307             :     !----------------------------- Local Variables -----------------------------
        1308             :     integer :: i, k
        1309             :
        1310             :     !----------------------------- Begin Code -----------------------------
        1311             :
        1312             :     !$acc data copyin( input_var1 ) &
        1313             :     !$acc     copyout( output_var )
        1314             :
        1315             :     !$acc parallel loop gang vector collapse(2) default(present)
        1316           0 :     do k = 1, nz
        1317           0 :       do i = 1, ngrdcol
        1318           0 :         output_var(i,k) = one_half * ( (input_var1(i,k)+input_var2) - &
        1319           0 :                                   sqrt((input_var1(i,k)-input_var2)**2 + smth_coef**2) )
        1320             :       end do
        1321             :     end do
        1322             :     !$acc end parallel loop
        1323             :
        1324             :     !$acc end data
        1325             :
        1326           0 :     return
        1327             :
        1328           0 :   end function smooth_min_array_scalar
        1329             :
        1330             : !===============================================================================
        1331           0 :   function smooth_min_arrays( nz, ngrdcol, input_var1, input_var2, smth_coef ) &
        1332           0 :   result( output_var )
        1333             :
        1334             :   ! Description:
        1335             :   !   Computes a smoothed version of the min function, using two 1d arrays as inputs.
        1336             :   !   For more details, see the interface in this file.
        1337             :
        1338             :   ! References:
        1339             :   !   See clubb:ticket:894, updated version: 965
        1340             :   !----------------------------------------------------------------------
        1341             :
        1342             :     use clubb_precision, only: &
        1343             :         core_rknd                     ! Constant(s)
        1344             :
        1345             :     use constants_clubb, only: &
        1346             :         one_half
        1347             :
        1348             :     implicit none
        1349             :
        1350             :     !----------------------------- Input Variables-----------------------------
        1351             :     integer, intent(in) :: &
        1352             :       nz, &
        1353             :       ngrdcol
        1354             :
        1355             :     real ( kind = core_rknd ), dimension(ngrdcol, nz), intent(in) :: &
        1356             :       input_var1, &       ! Units vary
        1357             :       input_var2          ! Units vary
        1358             :
        1359             :     real ( kind = core_rknd ), intent(in) :: &
        1360             :       smth_coef           ! "intensity" of the smoothing. Should be of a similar magnitude to
        1361             :                           ! that of the data structures input_var1 and input_var2
        1362             :
        1363             :     !----------------------------- Output Variables -----------------------------
        1364             :     real( kind = core_rknd ), dimension(ngrdcol, nz) :: &
        1365             :       output_var          ! Same unit as input_var1 and input_var2
        1366             :
        1367             :     !----------------------------- Local Variables -----------------------------
        1368             :     integer :: i, k
        1369             :
        1370             :     !----------------------------- Begin Code -----------------------------
        1371             :
        1372             :     !$acc data copyin( input_var1, input_var2 ) &
        1373             :     !$acc     copyout( output_var )
        1374             :
        1375             :     !$acc parallel loop gang vector collapse(2) default(present)
        1376           0 :     do k = 1, nz
        1377           0 :       do i = 1, ngrdcol
        1378           0 :         output_var(i,k) = one_half * ( (input_var1(i,k)+input_var2(i,k)) - &
        1379           0 :                                   sqrt((input_var1(i,k)-input_var2(i,k))**2 + smth_coef**2) )
        1380             :       end do
        1381             :     end do
        1382             :     !$acc end parallel loop
        1383             :
        1384             :     !$acc end data
        1385             :
        1386           0 :     return
        1387             :
        1388           0 :   end function smooth_min_arrays
        1389             :
        1390             : !===============================================================================
        1391           0 :   function smooth_min_scalars( input_var1, input_var2, smth_coef ) &
        1392             :   result( output_var )
        1393             :
        1394             :   ! Description:
        1395             :   !   Computes a smoothed version of the min function, using two scalars as inputs.
        1396             :   !   For more details, see the interface in this file.
        1397             :
        1398             :   ! References:
        1399             :   !   See clubb:ticket: 965
        1400             :   !----------------------------------------------------------------------
        1401             :
        1402             :     use clubb_precision, only: &
        1403             :         core_rknd                     ! Constant(s)
        1404             :
        1405             :     use constants_clubb, only: &
        1406             :         one_half
        1407             :
        1408             :     implicit none
        1409             :
        1410             :   ! Input Variables
        1411             :     real ( kind = core_rknd ), intent(in) :: &
        1412             :       input_var1, &       ! Units vary
        1413             :       input_var2, &       ! Units vary
        1414             :       smth_coef           ! "intensity" of the smoothing. Should be of a similar magnitude to
        1415             :                           ! that of the data structures input_var1 and input_var2
        1416             :
        1417             :   ! Output Variables
        1418             :     real( kind = core_rknd ) :: &
        1419             :       output_var          ! Same unit as input_var1 and input_var2
        1420             :
        1421             :   !----------------------------------------------------------------------
        1422             :
        1423             :     output_var = one_half * ( (input_var1+input_var2) - &
        1424           0 :                               sqrt((input_var1-input_var2)**2 + smth_coef**2) )
        1425             :
        1426             :     return
        1427             :   end function smooth_min_scalars
        1428             :
        1429             : !===============================================================================
        1430           0 :   function smooth_max_sclr_idx( nz, ngrdcol, input_var1, input_var2, smth_coef ) &
        1431           0 :   result( output_var )
        1432             :
        1433             :   ! Description:
        1434             :   !   Computes a smoothed version of the max function, using one scalar and
        1435             :   !   one 1d array as inputs. For more details, see the interface in this file.
        1436             :
        1437             :   ! References:
        1438             :   !   See clubb:ticket:894, updated version: 965
        1439             :   !----------------------------------------------------------------------
        1440             :
        1441             :     use clubb_precision, only: &
        1442             :         core_rknd                     ! Constant(s)
        1443             :
        1444             :     use constants_clubb, only: &
        1445             :         one_half
        1446             :
        1447             :     implicit none
        1448             :
        1449             :     !----------------------------- Input Variables -----------------------------
        1450             :     integer, intent(in) :: &
        1451             :       nz, &
        1452             :       ngrdcol
        1453             :
        1454             :     real ( kind = core_rknd ), intent(in) :: &
        1455             :       input_var1, &       ! Units vary
        1456             :       smth_coef           ! "intensity" of the smoothing. Should be of a similar magnitude to
        1457             :                           ! that of the data structures input_var1 and input_var2
        1458             :
        1459             :     real ( kind = core_rknd ), dimension(ngrdcol, nz), intent(in) :: &
        1460             :       input_var2          ! Units vary
        1461             :
        1462             :     !----------------------------- Output Variables -----------------------------
        1463             :     real( kind = core_rknd ), dimension(ngrdcol, nz) :: &
        1464             :       output_var          ! Same unit as input_var1 and input_var2
        1465             :
        1466             :     !----------------------------- Local Variables -----------------------------
        1467             :     integer :: i, k
        1468             :
        1469             :     !----------------------------- Begin Code -----------------------------
        1470             :
        1471             :     !$acc data copyin( input_var2 ) &
        1472             :     !$acc     copyout( output_var )
        1473             :
        1474             :     !$acc parallel loop gang vector collapse(2) default(present)
        1475           0 :     do k = 1, nz
        1476           0 :       do i = 1, ngrdcol
        1477           0 :         output_var(i,k) = one_half * ( (input_var1+input_var2(i,k)) + &
        1478           0 :                                   sqrt((input_var1-input_var2(i,k))**2 + smth_coef**2) )
        1479             :       end do
        1480             :     end do
        1481             :     !$acc end parallel loop
        1482             :
        1483             :     !$acc end data
        1484             :
        1485           0 :     return
        1486             :
        1487           0 :   end function smooth_max_sclr_idx
        1488             :
        1489             : !===============================================================================
        1490           0 :   function smooth_max_array_scalar( nz, ngrdcol, input_var1, input_var2, smth_coef ) &
        1491           0 :   result( output_var )
        1492             :
        1493             :   ! Description:
        1494             :   !   Computes a smoothed version of the max function, using one scalar and
        1495             :   !   one 1d array as inputs. For more details, see the interface in this file.
        1496             :
        1497             :   ! References:
        1498             :   !   See clubb:ticket:894, updated version: 965
        1499             :   !----------------------------------------------------------------------
        1500             :
        1501             :     use clubb_precision, only: &
        1502             :         core_rknd                     ! Constant(s)
        1503             :
        1504             :     use constants_clubb, only: &
        1505             :         one_half
        1506             :
        1507             :     implicit none
        1508             :
        1509             :     !----------------------------- Input Variables -----------------------------
        1510             :     integer, intent(in) :: &
        1511             :       nz, &
        1512             :       ngrdcol
        1513             :
        1514             :     real ( kind = core_rknd ), dimension(ngrdcol, nz), intent(in) :: &
        1515             :       input_var1          ! Units vary
        1516             :
        1517             :     real ( kind = core_rknd ), intent(in) :: &
        1518             :       input_var2, &       ! Units vary
        1519             :       smth_coef           ! "intensity" of the smoothing. Should be of a similar magnitude to
        1520             :                           ! that of the data structures input_var1 and input_var2
        1521             :
        1522             :     !----------------------------- Output Variables -----------------------------
        1523             :     real( kind = core_rknd ), dimension(ngrdcol, nz) :: &
        1524             :       output_var          ! Same unit as input_var1 and input_var2
        1525             :
        1526             :     !----------------------------- Local Variables -----------------------------
        1527             :     integer :: i, k
        1528             :
        1529             :     !----------------------------- Begin Code -----------------------------
        1530             :
        1531             :     !$acc data copyin( input_var1 ) &
        1532             :     !$acc     copyout( output_var )
        1533             :
        1534             :     !$acc parallel loop gang vector collapse(2) default(present)
        1535           0 :     do k = 1, nz
        1536           0 :       do i = 1, ngrdcol
        1537           0 :         output_var(i,k) = one_half * ( ( input_var1(i,k) + input_var2 ) + &
        1538           0 :                                   sqrt(( input_var1(i,k) - input_var2 )**2 + smth_coef**2) )
        1539             :       end do
        1540             :     end do
        1541             :     !$acc end parallel loop
        1542             :
        1543             :     !$acc end data
        1544             :
        1545           0 :     return
        1546             :
        1547           0 :   end function smooth_max_array_scalar
        1548             :
        1549             : !===============================================================================
        1550           0 :   function smooth_max_array_1D_scalar( ngrdcol, input_var1, input_var2, smth_coef ) &
        1551           0 :   result( output_var )
        1552             :
        1553             :   ! Description:
        1554             :   !   Computes a smoothed version of the max function, using one scalar and
        1555             :   !   one 1d array as inputs. For more details, see the interface in this file.
        1556             :
        1557             :   ! References:
        1558             :   !   See clubb:ticket:894, updated version: 965
        1559             :   !----------------------------------------------------------------------
        1560             :
        1561             :     use clubb_precision, only: &
        1562             :         core_rknd                     ! Constant(s)
        1563             :
        1564             :     use constants_clubb, only: &
        1565             :         one_half
        1566             :
        1567             :     implicit none
        1568             :
        1569             :     !----------------------------- Input Variables -----------------------------
        1570             :     integer, intent(in) :: &
        1571             :       ngrdcol
        1572             :
        1573             :     real ( kind = core_rknd ), dimension(ngrdcol), intent(in) :: &
        1574             :       input_var1          ! Units vary
        1575             :
        1576             :     real ( kind = core_rknd ), intent(in) :: &
        1577             :       input_var2, &       ! Units vary
        1578             :       smth_coef           ! "intensity" of the smoothing. Should be of a similar magnitude to
        1579             :                           ! that of the data structures input_var1 and input_var2
        1580             :
        1581             :     !----------------------------- Output Variables -----------------------------
        1582             :     real( kind = core_rknd ), dimension(ngrdcol) :: &
        1583             :       output_var          ! Same unit as input_var1 and input_var2
        1584             :
        1585             :     !----------------------------- Local Variables -----------------------------
        1586             :     integer :: i
        1587             :
        1588             :     !----------------------------- Begin Code -----------------------------
        1589             :
        1590             :     !$acc data copyin( input_var1 ) &
        1591             :     !$acc     copyout( output_var )
        1592             :
        1593             :     !$acc parallel loop gang vector default(present)
        1594           0 :     do i = 1, ngrdcol
        1595           0 :       output_var(i) = one_half * ( ( input_var1(i) + input_var2 ) + &
        1596           0 :                                 sqrt(( input_var1(i) - input_var2 )**2 + smth_coef**2) )
        1597             :     end do
        1598             :     !$acc end parallel loop
        1599             :
        1600             :     !$acc end data
        1601             :
        1602           0 :     return
        1603             :
        1604           0 :   end function smooth_max_array_1D_scalar
        1605             :
        1606             : !===============================================================================
        1607           0 :   function smooth_max_arrays( nz, ngrdcol, input_var1, input_var2, smth_coef ) &
        1608           0 :   result( output_var )
        1609             :
        1610             :   ! Description:
        1611             :   !   Computes a smoothed version of the max function, using two 1d arrays as inputs.
        1612             :   !   For more details, see the interface in this file.
        1613             :
        1614             :   ! References:
        1615             :   !   See clubb:ticket:894, updated version: 965
        1616             :   !----------------------------------------------------------------------
        1617             :
        1618             :     use clubb_precision, only: &
        1619             :         core_rknd                     ! Constant(s)
        1620             :
        1621             :     use constants_clubb, only: &
        1622             :         one_half
        1623             :
        1624             :     implicit none
        1625             :
        1626             :     !----------------------------- Input Variables -----------------------------
        1627             :     integer, intent(in) :: &
        1628             :       nz, &
        1629             :       ngrdcol
        1630             :
        1631             :     real ( kind = core_rknd ), dimension(ngrdcol, nz), intent(in) :: &
        1632             :       input_var1, &       ! Units vary
        1633             :       input_var2          ! Units vary
        1634             :
        1635             :     real( kind = core_rknd ), intent(in) :: &
        1636             :       smth_coef           ! "intensity" of the smoothing. Should be of a similar magnitude to
        1637             :                           ! that of the data structures input_var1 and input_var2
        1638             :
        1639             :     !----------------------------- Output Variables -----------------------------
        1640             :     real( kind = core_rknd ), dimension(ngrdcol, nz) :: &
        1641             :       output_var          ! Same unit as input_var1 and input_var2
        1642             :
        1643             :     !----------------------------- Local Variables -----------------------------
        1644             :     integer :: i, k
        1645             :
        1646             :     !----------------------------- Begin Code -----------------------------
        1647             :
        1648             :     !$acc data copyin( input_var1, input_var2 ) &
        1649             :     !$acc     copyout( output_var )
        1650             :
        1651             :     !$acc parallel loop gang vector collapse(2) default(present)
        1652           0 :     do k = 1, nz
        1653           0 :       do i = 1, ngrdcol
        1654           0 :         output_var(i,k) = one_half * ( (input_var1(i,k)+input_var2(i,k)) + &
        1655           0 :                                   sqrt((input_var1(i,k)-input_var2(i,k))**2 + smth_coef**2) )
        1656             :       end do
        1657             :     end do
        1658             :     !$acc end parallel loop
        1659             :
        1660             :     !$acc end data
        1661             :
        1662           0 :     return
        1663             :
        1664           0 :   end function smooth_max_arrays
        1665             :
        1666             : !===============================================================================
        1667           0 :   function smooth_max_scalars( input_var1, input_var2, smth_coef ) &
        1668             :   result( output_var )
        1669             :
        1670             :   ! Description:
        1671             :   !   Computes a smoothed version of the max function, using two scalars as inputs.
        1672             :   !   For more details, see the interface in this file.
        1673             :
        1674             :   ! References:
        1675             :   !   See clubb:ticket: 965
        1676             :   !----------------------------------------------------------------------
        1677             :
        1678             :     use clubb_precision, only: &
        1679             :         core_rknd                     ! Constant(s)
        1680             :
        1681             :     use constants_clubb, only: &
        1682             :         one_half
        1683             :
        1684             :     implicit none
        1685             :
        1686             :     !----------------------------- Input Variables -----------------------------
        1687             :     real ( kind = core_rknd ), intent(in) :: &
        1688             :       input_var1, &       ! Units vary
        1689             :       input_var2, &       ! Units vary
        1690             :       smth_coef           ! "intensity" of the smoothing. Should be of a similar magnitude to
        1691             :                           ! that of the data structures input_var1 and input_var2
        1692             :
        1693             :     !----------------------------- Output Variables -----------------------------
        1694             :     real( kind = core_rknd ) :: &
        1695             :       output_var          ! Same unit as input_var1 and input_var2
        1696             :
        1697             :     !----------------------------- Local Variables -----------------------------
        1698             :     integer :: i, k
        1699             :
        1700             :     !----------------------------- Begin Code -----------------------------
        1701             :
        1702             :     output_var = one_half * ( (input_var1+input_var2) + &
        1703           0 :                               sqrt((input_var1-input_var2)**2 + smth_coef**2) )
        1704             :     return
        1705             :
        1706             :   end function smooth_max_scalars
        1707             :
        1708           0 :   function smooth_heaviside_peskin( nz, ngrdcol, input, smth_range ) &
        1709           0 :     result( smth_output )
        1710             :
        1711             :   ! Description:
        1712             :   !   Computes a smoothed heaviside function as in
        1713             :   !       [Lin, Lee et al., 2005, A level set characteristic Galerkin
        1714             :   !       finite element method for free surface flows], equation (2)
        1715             :
        1716             :   ! References:
        1717             :   !   See clubb:ticket:965
        1718             :   !----------------------------------------------------------------------
        1719             :
        1720             :     use clubb_precision, only: &
        1721             :         core_rknd                     ! Constant(s)
        1722             :
        1723             :     use constants_clubb, only: &
        1724             :         pi, invrs_pi, one, one_half, zero
        1725             :
        1726             :     implicit none
        1727             :
        1728             :     !------------------------- Input Variables -------------------------
        1729             :     integer, intent(in) :: &
        1730             :       nz, &
        1731             :       ngrdcol
        1732             :
        1733             :     real ( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
        1734             :       input    ! Units vary
        1735             :
        1736             :     real ( kind = core_rknd ), intent(in) :: &
        1737             :       smth_range  ! Smooth Heaviside function on [-smth_range, smth_range]
        1738             :
        1739             :     !------------------------- Output Variables -------------------------
        1740             :     real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
        1741             :       smth_output    ! Same units as input
        1742             :
        1743             :     !------------------------- Local Variables -------------------------
        1744             :     real ( kind = core_rknd ) :: &
        1745             :       input_over_smth_range  ! input divided by smth_range
        1746             :
        1747             :     integer :: i, k
        1748             :
        1749             :     !------------------------- Begin Code -------------------------
        1750             :
        1751             :     !$acc data copyin( input ) &
        1752             :     !$acc     copyout( smth_output )
        1753             :
        1754             :     !$acc parallel loop gang vector collapse(2) default(present)
        1755           0 :     do k = 1, nz
        1756           0 :       do i = 1, ngrdcol
        1757             :
        1758           0 :         if ( input(i,k) < -smth_range ) then
        1759           0 :           smth_output(i,k) = zero
        1760           0 :         elseif ( input(i,k) > smth_range ) then
        1761           0 :            smth_output(i,k) = one
        1762             :         else
        1763             :           ! Note that this case will only ever be reached if smth_range != 0,
        1764             :           ! so this division is fine and should not cause any issues
        1765           0 :           input_over_smth_range = input(i,k) / smth_range
        1766             :           smth_output(i,k) = one_half &
        1767             :                              * (one + input_over_smth_range &
        1768           0 :                                + invrs_pi * sin(pi * input_over_smth_range))
        1769             :         end if
        1770             :
        1771             :       end do
        1772             :     end do
        1773             :     !$acc end parallel loop
        1774             :
        1775             :     !$acc end data
        1776             :
        1777           0 :     return
        1778             :
        1779           0 :   end function smooth_heaviside_peskin
        1780             :
        1781             :   !===============================================================================
        1782           0 :   subroutine calc_xpwp_1D( gr, Km_zm, xm, &
        1783           0 :                            xpwp )
        1784             :
        1785             :     ! Description:
        1786             :     ! Compute x'w' from x<k>, x<k+1>, Kh and invrs_dzm
        1787             :
        1788             :     ! References:
        1789             :     ! None
        1790             :     !-----------------------------------------------------------------------
        1791             :
        1792             :     use clubb_precision, only: &
        1793             :         core_rknd ! Variable(s)
        1794             :
        1795             :     use grid_class, only: &
        1796             :       grid
        1797             :
        1798             :     implicit none
        1799             :
        1800             :     ! ----------------------- Input variables -----------------------
        1801             :     type (grid), target, intent(in) :: gr
        1802             :
        1803             :     real( kind = core_rknd ), dimension(gr%nz), intent(in) :: &
        1804             :       Km_zm,     & ! Eddy diff. (k momentum level)                 [m^2/s]
        1805             :       xm           ! x (k thermo level)                            [units vary]
        1806             :
        1807             :     ! ----------------------- Output variable -----------------------
        1808             :     real( kind = core_rknd ), dimension(gr%nz), intent(out) :: &
        1809             :       xpwp ! x'w'   [(units vary)(m/s)]
        1810             :
        1811             :     integer :: k
        1812             :
        1813             :     ! ----------------------- Begin Code -----------------------
        1814             :
        1815             :     ! Solve for x'w' at all intermediate model levels.
        1816           0 :     do k = 1, gr%nz-1
        1817           0 :       xpwp(k) = Km_zm(k) * gr%invrs_dzm(1,k) * ( xm(k+1) - xm(k) )
        1818             :     end do
        1819             :
        1820           0 :     return
        1821             :   end subroutine calc_xpwp_1D
        1822             :
        1823             :   !===============================================================================
        1824   232311456 :   subroutine calc_xpwp_2D( nz, ngrdcol, gr, &
        1825   232311456 :                         Km_zm, xm, &
        1826   232311456 :                         xpwp )
        1827             :
        1828             :     ! Description:
        1829             :     ! Compute x'w' from x<k>, x<k+1>, Kh and invrs_dzm
        1830             :
        1831             :     ! References:
        1832             :     ! None
        1833             :     !-----------------------------------------------------------------------
        1834             :
        1835             :     use clubb_precision, only: &
        1836             :         core_rknd ! Variable(s)
        1837             :
        1838             :     use grid_class, only: &
        1839             :       grid
        1840             :
        1841             :     implicit none
        1842             :
        1843             :     ! ----------------------- Input variables -----------------------
        1844             :     integer, intent(in) :: &
        1845             :       nz, &
        1846             :       ngrdcol
        1847             :
        1848             :     type (grid), target, intent(in) :: gr
        1849             :
        1850             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
        1851             :       Km_zm,     & ! Eddy diff. (k momentum level)                 [m^2/s]
        1852             :       xm           ! x (k thermo level)                            [units vary]
        1853             :
        1854             :     ! ----------------------- Output variable -----------------------
        1855             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
        1856             :       xpwp ! x'w'   [(units vary)(m/s)]
        1857             :
        1858             :     integer :: i, k
        1859             :
        1860             :     ! ----------------------- Begin Code -----------------------
        1861             :
        1862             :     !$acc data copyin( gr, gr%invrs_dzm, Km_zm, xm ) &
        1863             :     !$acc     copyout( xpwp )
        1864             :
        1865             :     ! Solve for x'w' at all intermediate model levels.
        1866             :     !$acc parallel loop gang vector collapse(2) default(present)
        1867 19746473760 :     do k = 1, nz-1
        1868 >32607*10^7 :       do i = 1, ngrdcol
        1869 >32584*10^7 :         xpwp(i,k) = Km_zm(i,k) * gr%invrs_dzm(i,k) * ( xm(i,k+1) - xm(i,k) )
        1870             :       end do
        1871             :     end do
        1872             :     !$acc end parallel loop
        1873             :
        1874             :     !$acc end data
        1875             :
        1876   232311456 :     return
        1877             :
        1878             :   end subroutine calc_xpwp_2D
        1879             :
        1880             :   !=============================================================================
        1881           0 :   function vertical_avg( total_idx, rho_ds, field, dz )
        1882             :
        1883             :     ! Description:
        1884             :     ! Computes the density-weighted vertical average of a field.
        1885             :     !
        1886             :     ! The average value of a function, f, over a set domain, [a,b], is
        1887             :     ! calculated by the equation:
        1888             :     !
        1889             :     ! f_avg = ( INT(a:b) f*g ) / ( INT(a:b) g );
        1890             :     !
        1891             :     ! as long as f is continous and g is nonnegative and integrable.  Therefore,
        1892             :     ! the density-weighted (by dry, static, base-static density) vertical
        1893             :     ! average value of any model field, x, is calculated by the equation:
        1894             :     !
        1895             :     ! x_avg|_z = ( INT(z_bot:z_top) x rho_ds dz )
        1896             :     !            / ( INT(z_bot:z_top) rho_ds dz );
        1897             :     !
        1898             :     ! where z_bot is the bottom of the vertical domain, and z_top is the top of
        1899             :     ! the vertical domain.
        1900             :     !
        1901             :     ! This calculation is done slightly differently depending on whether x is a
        1902             :     ! thermodynamic-level or a momentum-level variable.
        1903             :     !
        1904             :     ! Thermodynamic-level computation:
        1905             :
        1906             :     !
        1907             :     ! For numerical purposes, INT(z_bot:z_top) x rho_ds dz, which is the
        1908             :     ! numerator integral, is calculated as:
        1909             :     !
        1910             :     ! SUM(k_bot:k_top) x(k) rho_ds(k) delta_z(k);
        1911             :     !
        1912             :     ! where k is the index of the given thermodynamic level, x and rho_ds are
        1913             :     ! both thermodynamic-level variables, and delta_z(k) = zm(k) - zm(k-1).  The
        1914             :     ! indices k_bot and k_top are the indices of the respective lower and upper
        1915             :     ! thermodynamic levels involved in the integration.
        1916             :     !
        1917             :     ! Likewise, INT(z_bot:z_top) rho_ds dz, which is the denominator integral,
        1918             :     ! is calculated as:
        1919             :     !
        1920             :     ! SUM(k_bot:k_top) rho_ds(k) delta_z(k).
        1921             :     !
        1922             :     ! The first (k=1) thermodynamic level is below ground (or below the
        1923             :     ! official lower boundary at the first momentum level), so it should not
        1924             :     ! count in a vertical average, whether that vertical average is used for
        1925             :     ! the hole-filling scheme or for statistical purposes. Begin no lower
        1926             :     ! than level k=2, which is the first thermodynamic level above ground (or
        1927             :     ! above the model lower boundary).
        1928             :     !
        1929             :     ! For cases where hole-filling over the entire (global) vertical domain
        1930             :     ! is desired, or where statistics over the entire (global) vertical
        1931             :     ! domain are desired, the lower (thermodynamic-level) index of k = 2 and
        1932             :     ! the upper (thermodynamic-level) index of k = gr%nz, means that the
        1933             :     ! overall vertical domain will be gr%zm(1,gr%nz) - gr%zm(1,1).
        1934             :     !
        1935             :     !
        1936             :     ! Momentum-level computation:
        1937             :     !
        1938             :     ! For numerical purposes, INT(z_bot:z_top) x rho_ds dz, which is the
        1939             :     ! numerator integral, is calculated as:
        1940             :     !
        1941             :     ! SUM(k_bot:k_top) x(k) rho_ds(k) delta_z(k);
        1942             :     !
        1943             :     ! where k is the index of the given momentum level, x and rho_ds are both
        1944             :     ! momentum-level variables, and delta_z(k) = zt(k+1) - zt(k).  The indices
        1945             :     ! k_bot and k_top are the indices of the respective lower and upper momentum
        1946             :     ! levels involved in the integration.
        1947             :     !
        1948             :     ! Likewise, INT(z_bot:z_top) rho_ds dz, which is the denominator integral,
        1949             :     ! is calculated as:
        1950             :     !
        1951             :     ! SUM(k_bot:k_top) rho_ds(k) delta_z(k).
        1952             :     !
        1953             :     ! The first (k=1) momentum level is right at ground level (or right at
        1954             :     ! the official lower boundary).  The momentum level variables that call
        1955             :     ! the hole-filling scheme have set values at the surface (or lower
        1956             :     ! boundary), and those set values should not be changed.  Therefore, the
        1957             :     ! vertical average (for purposes of hole-filling) should not include the
        1958             :     ! surface level (or lower boundary level).  For hole-filling purposes,
        1959             :     ! begin no lower than level k=2, which is the second momentum level above
        1960             :     ! ground (or above the model lower boundary).  Likewise, the value at the
        1961             :     ! model upper boundary (k=gr%nz) is also set for momentum level
        1962             :     ! variables.  That value should also not be changed.
        1963             :     !
        1964             :     ! However, this function is also used to keep track (for statistical
        1965             :     ! purposes) of the vertical average of certain variables.  In that case,
        1966             :     ! the vertical average needs to be taken over the entire vertical domain
        1967             :     ! (level 1 to level gr%nz).
        1968             :     !
        1969             :     !
        1970             :     ! In both the thermodynamic-level computation and the momentum-level
        1971             :     ! computation, the numerator integral is divided by the denominator integral
        1972             :     ! in order to find the average value (over the vertical domain) of x.
        1973             :
        1974             :     ! References:
        1975             :     ! None
        1976             :     !-----------------------------------------------------------------------
        1977             :
        1978             :     use clubb_precision, only: &
        1979             :         core_rknd ! Variable(s)
        1980             :
        1981             :     implicit none
        1982             :
        1983             :     ! Input variables
        1984             :     integer, intent(in) :: &
        1985             :       total_idx ! The total numer of indices within the range of averaging
        1986             :
        1987             :     real( kind = core_rknd ), dimension(total_idx), intent(in) ::  &
        1988             :       rho_ds, & ! Dry, static density on either thermodynamic or momentum levels    [kg/m^3]
        1989             :       field,  & ! The field (e.g. wp2) to be vertically averaged                    [Units vary]
        1990             :       dz  ! Reciprocal of thermodynamic or momentum level thickness           [1/m]
        1991             :                 ! depending on whether we're on zt or zm grid.
        1992             :     ! Note:  The rho_ds and field points need to be arranged from
        1993             :     !        lowest to highest in altitude, with rho_ds(1) and
        1994             :     !        field(1) actually their respective values at level k = 1.
        1995             :
        1996             :     ! Output variable
        1997             :     real( kind = core_rknd ) :: &
        1998             :       vertical_avg  ! Vertical average of field    [Units of field]
        1999             :
        2000             :     ! Local variables
        2001             :     real( kind = core_rknd ) :: &
        2002             :       numer_integral, & ! Integral in the numerator (see description)
        2003             :       denom_integral    ! Integral in the denominator (see description)
        2004             :
        2005             :
        2006             :     integer :: k
        2007             :
        2008             :     !-----------------------------------------------------------------------
        2009             :
        2010             :     ! Initialize variable
        2011           0 :     numer_integral = 0.0_core_rknd
        2012           0 :     denom_integral = 0.0_core_rknd
        2013             :
        2014             :     ! Compute the numerator and denominator integral.
        2015             :     ! Multiply rho_ds at level k by the level thickness
        2016             :     ! at level k.  Then, sum over all vertical levels.
        2017           0 :     do k=1, total_idx
        2018             :
        2019           0 :         numer_integral = numer_integral + rho_ds(k) * dz(k) * field(k)
        2020           0 :         denom_integral = denom_integral + rho_ds(k) * dz(k)
        2021             :
        2022             :     end do
        2023             :
        2024             :     ! Find the vertical average of 'field'.
        2025           0 :     vertical_avg = numer_integral / denom_integral
        2026             :     !vertical_avg = sum( rho_ds(:) * dz(:) * field(:) ) / sum( rho_ds(:) * dz(:) )
        2027             :
        2028             :     return
        2029             :   end function vertical_avg
        2030             :
        2031             :   !=============================================================================
        2032           0 :   function vertical_integral( total_idx, rho_ds, &
        2033           0 :                                        field, dz )
        2034             :
        2035             :     ! Description:
        2036             :     ! Computes the vertical integral. rho_ds, field, and dz must all be
        2037             :     ! of size total_idx and should all start at the same index.
        2038             :     !
        2039             :
        2040             :     ! References:
        2041             :     ! None
        2042             :     !-----------------------------------------------------------------------
        2043             :
        2044             :     use clubb_precision, only: &
        2045             :         core_rknd ! Variable(s)
        2046             :
        2047             :     implicit none
        2048             :
        2049             :     ! Input variables
        2050             :     integer, intent(in) :: &
        2051             :       total_idx  ! The total numer of indices within the range of averaging
        2052             :
        2053             :     real( kind = core_rknd ), dimension(total_idx), intent(in) ::  &
        2054             :       rho_ds,  & ! Dry, static density                   [kg/m^3]
        2055             :       field,   & ! The field to be vertically averaged   [Units vary]
        2056             :       dz         ! Level thickness                       [1/m]
        2057             :     ! Note:  The rho_ds and field points need to be arranged from
        2058             :     !        lowest to highest in altitude, with rho_ds(1) and
        2059             :     !        field(1) actually their respective values at level k = k_start.
        2060             :
        2061             :     ! Local variables
        2062             :     real( kind = core_rknd ) :: &
        2063             :       vertical_integral ! Integral in the numerator (see description)
        2064             :
        2065             :     !-----------------------------------------------------------------------
        2066             :
        2067             :     !  Assertion checks: that k_start <= gr%nz - 1
        2068             :     !                    that k_end   >= 2
        2069             :     !                    that k_start <= k_end
        2070             :
        2071             :
        2072             :     ! Initializing vertical_integral to avoid a compiler warning.
        2073           0 :     vertical_integral = 0.0_core_rknd
        2074             :
        2075             :     ! Compute the integral.
        2076             :     ! Multiply the field at level k by rho_ds at level k and by
        2077             :     ! the level thickness at level k.  Then, sum over all vertical levels.
        2078             :     ! Note:  The values of the field and rho_ds are passed into this function
        2079             :     !        so that field(1) and rho_ds(1) are actually the field and rho_ds
        2080             :     !        at level k_start.
        2081           0 :     vertical_integral = sum( field * rho_ds * dz )
        2082             :
        2083             :     !print *, vertical_integral
        2084             :
        2085             :     return
        2086             :   end function vertical_integral
        2087             :
        2088             :
        2089             : end module advance_helper_module

Generated by: LCOV version 1.14