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