LCOV - code coverage report
Current view: | top level - physics/clubb/src/CLUBB_core - adg1_adg2_3d_luhar_pdf.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 76 | 203 | 37.4 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 3 | 9 | 33.3 %

              Line data    Source code
           1             : ! $Id$
           2             : !===============================================================================
           3             : module adg1_adg2_3d_luhar_pdf
           4             :
           5             :   ! Description:
           6             :   ! Contains code to calculate the PDF parameters using the Analytic Double
           7             :   ! Gaussian 1 (ADG1) PDF closure, the Analytic Double Gaussian 2 (ADG2) PDF
           8             :   ! closure, or the 3D Luhar PDF closure.  The ADG1 PDF is the original PDF used
           9             :   ! by CLUBB in interactive runs.
          10             :
          11             :   ! References:
          12             :   !-------------------------------------------------------------------------
          13             :
          14             :   implicit none
          15             :
          16             :   public :: ADG1_pdf_driver,     & ! Procedure(s)
          17             :             ADG2_pdf_driver,     &
          18             :             Luhar_3D_pdf_driver, &
          19             :             ADG1_w_closure,      &
          20             :             calc_Luhar_params,   &
          21             :             close_Luhar_pdf
          22             :
          23             :   private :: ADG1_ADG2_responder_params, & ! Procedure(s)
          24             :              backsolve_Luhar_params,     &
          25             :              max_cubic_root
          26             :
          27             :   private
          28             :
          29             :   contains
          30             :
          31             :   !=============================================================================
          32    17870112 :   subroutine ADG1_pdf_driver( nz, ngrdcol, sclr_dim, sclr_tol,          & ! In
          33    17870112 :                               wm, rtm, thlm, um, vm,                    & ! In
          34    17870112 :                               wp2, rtp2, thlp2, up2, vp2,               & ! In
          35    17870112 :                               Skw, wprtp, wpthlp, upwp, vpwp, sqrt_wp2, & ! In
          36    17870112 :                               sigma_sqd_w, beta, mixt_frac_max_mag,     & ! In
          37    17870112 :                               sclrm, sclrp2, wpsclrp, l_scalar_calc,    & ! In
          38    17870112 :                               w_1, w_2,                                 & ! Out
          39    17870112 :                               rt_1, rt_2,                               & ! Out
          40    17870112 :                               thl_1, thl_2,                             & ! Out
          41    17870112 :                               u_1, u_2, v_1, v_2,                       & ! Out
          42    17870112 :                               varnce_w_1, varnce_w_2,                   & ! Out
          43    17870112 :                               varnce_rt_1, varnce_rt_2,                 & ! Out
          44    17870112 :                               varnce_thl_1, varnce_thl_2,               & ! Out
          45    17870112 :                               varnce_u_1, varnce_u_2,                   & ! Out
          46    17870112 :                               varnce_v_1, varnce_v_2,                   & ! Out
          47    17870112 :                               mixt_frac,                                & ! Out
          48    17870112 :                               alpha_rt, alpha_thl,                      & ! Out
          49    17870112 :                               alpha_u, alpha_v,                         & ! Out
          50    17870112 :                               sclr_1, sclr_2, varnce_sclr_1,            & ! Out
          51    17870112 :                               varnce_sclr_2, alpha_sclr )                 ! Out
          52             :
          53             :     ! Calculates the PDF parameters using the Analytic Double Gaussian 1 (ADG1)
          54             :     ! PDF closure.
          55             :
          56             :     ! References:
          57             :     !-----------------------------------------------------------------------
          58             :
          59             :     use constants_clubb, only: &
          60             :         rt_tol,  & ! Constant(s)
          61             :         thl_tol
          62             :
          63             :     use clubb_precision, only: &
          64             :         core_rknd    ! Variable(s)
          65             :
          66             :     implicit none
          67             :
          68             :     integer, intent(in) :: &
          69             :       nz,           & ! Number of vertical levels
          70             :       ngrdcol,      & ! Number of grid columns
          71             :       sclr_dim        ! Number of passive scalars
          72             :
          73             :     real( kind = core_rknd ), intent(in), dimension(sclr_dim) :: &
          74             :       sclr_tol          ! Threshold(s) on the passive scalars  [units vary]
          75             :
          76             :     ! Input Variables
          77             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) ::  &
          78             :       wm,          & ! Mean of w-wind comp. (vert. vel.)    [m/s]
          79             :       rtm,         & ! Mean of total water mixing ratio     [kg/kg]
          80             :       thlm,        & ! Mean of liquid water potential temp. [K]
          81             :       um,          & ! Mean of eastward wind                [m/s]
          82             :       vm,          & ! Mean of northward wind               [m/s]
          83             :       wp2,         & ! Variance of w (overall)              [m^2/s^2]
          84             :       rtp2,        & ! Variance of r_t (overall)            [(kg/kg)^2]
          85             :       thlp2,       & ! Variance of th_l (overall)           [K^2]
          86             :       up2,         & ! Variance of eastward wind (overall)  [m^2/s^2]
          87             :       vp2,         & ! Variance of northward wind (overall) [m^2/s^2]
          88             :       Skw,         & ! Skewness of w                        [-]
          89             :       wprtp,       & ! Covariance of w and r_t              [(kg/kg)(m/s)]
          90             :       wpthlp,      & ! Covariance of w and th_l             [K(m/s)]
          91             :       upwp,        & ! Covariance of u and w                [m^2/s^2]
          92             :       vpwp,        & ! Covariance of v and w                [m^2/s^2]
          93             :       sqrt_wp2,    & ! Square root of variance of w         [m/s]
          94             :       sigma_sqd_w    ! Width of individual w plumes         [-]
          95             :
          96             :     real( kind = core_rknd ), dimension(ngrdcol), intent(in) ::  &
          97             :       beta                    ! CLUBB tunable parameter beta         [-]
          98             :
          99             :     real( kind = core_rknd ), intent(in) ::  &
         100             :       mixt_frac_max_mag       ! Maximum allowable mag. of mixt_frac  [-]
         101             :
         102             :     real( kind = core_rknd ), dimension(ngrdcol, nz, sclr_dim), intent(in) ::  &
         103             :       sclrm,   & ! Mean of passive scalar (overall)           [units vary]
         104             :       sclrp2,  & ! Variance of passive scalar (overall)       [(units vary)^2]
         105             :       wpsclrp    ! Covariance of w and passive scalar         [m/s (units vary)]
         106             :
         107             :     logical, intent(in) :: &
         108             :       l_scalar_calc    ! Flag to perform calculations for passive scalars
         109             :
         110             :     ! Output Variables
         111             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) ::  &
         112             :       w_1,          & ! Mean of w (1st PDF component)                      [m/s]
         113             :       w_2,          & ! Mean of w (2nd PDF component)                      [m/s]
         114             :       rt_1,         & ! Mean of r_t (1st PDF component)                  [kg/kg]
         115             :       rt_2,         & ! Mean of r_t (2nd PDF component)                  [kg/kg]
         116             :       thl_1,        & ! Mean of th_l (1st PDF component)                     [K]
         117             :       thl_2,        & ! Mean of th_l (2nd PDF component)                     [K]
         118             :       u_1,          & ! Mean of u (eastward wind) (1st PDF component)      [m/s]
         119             :       u_2,          & ! Mean of u (eastward wind) (2nd PDF component)      [m/s]
         120             :       v_1,          & ! Mean of v (northward wind) (1st PDF component)     [m/s]
         121             :       v_2,          & ! Mean of v (northward wind) (2nd PDF component)     [m/s]
         122             :       varnce_w_1,   & ! Variance of w (1st PDF component)              [m^2/s^2]
         123             :       varnce_w_2,   & ! Variance of w (2nd PDF component)              [m^2/s^2]
         124             :       varnce_rt_1,  & ! Variance of r_t (1st PDF component)          [kg^2/kg^2]
         125             :       varnce_rt_2,  & ! Variance of r_t (2nd PDF component)          [kg^2/kg^2]
         126             :       varnce_thl_1, & ! Variance of th_l (1st PDF component)               [K^2]
         127             :       varnce_thl_2, & ! Variance of th_l (2nd PDF component)               [K^2]
         128             :       varnce_u_1,   & ! Variance of u wind (1st PDF component)         [m^2/s^2]
         129             :       varnce_u_2,   & ! Variance of u wind (2nd PDF component)         [m^2/s^2]
         130             :       varnce_v_1,   & ! Variance of v wind (1st PDF component)         [m^2/s^2]
         131             :       varnce_v_2,   & ! Variance of v wind (2nd PDF component)         [m^2/s^2]
         132             :       mixt_frac,    & ! Mixture fraction (weight of 1st PDF component)       [-]
         133             :       alpha_thl,    & ! Factor relating to normalized variance for th_l      [-]
         134             :       alpha_rt,     & ! Factor relating to normalized variance for r_t       [-]
         135             :       alpha_u,      & ! Factor relating to normalized variance for u wind    [-]
         136             :       alpha_v         ! Factor relating to normalized variance for v wind    [-]
         137             :
         138             :     real( kind = core_rknd ), dimension(ngrdcol, nz, sclr_dim), intent(out) ::  &
         139             :       sclr_1,        & ! Mean of passive scalar (1st PDF component) [units vary]
         140             :       sclr_2,        & ! Mean of passive scalar (2nd PDF component) [units vary]
         141             :       varnce_sclr_1, & ! Variance of pass. sclr (1st PDF comp.) [(units vary)^2]
         142             :       varnce_sclr_2, & ! Variance of pass. sclr (2nd PDF comp.) [(units vary)^2]
         143             :       alpha_sclr       ! Factor relating to normalized variance for sclr     [-]
         144             :
         145             :     ! Local Variables
         146             :     real( kind = core_rknd ), dimension(ngrdcol,nz) ::  &
         147    35740224 :       w_1_n, & ! Normalized mean of w (1st PDF component)     [-]
         148    35740224 :       w_2_n    ! Normalized mean of w (2nd PDF component)     [-]
         149             :
         150             :     integer :: j  ! Loop index
         151             :
         152             :     !$acc enter data create( w_1_n, w_2_n )
         153             :
         154             :     ! Calculate the mixture fraction and the PDF component means and variances
         155             :     ! of w.
         156             :     call ADG1_w_closure( nz, ngrdcol, wm, wp2, Skw, sigma_sqd_w, &        ! In
         157             :                          sqrt_wp2, mixt_frac_max_mag, &      ! In
         158             :                          w_1, w_2, w_1_n, w_2_n, &           ! Out
         159    17870112 :                          varnce_w_1, varnce_w_2, mixt_frac ) ! Out
         160             :
         161             :     ! Calculate the PDF component means and variances of rt.
         162             :     call ADG1_ADG2_responder_params( nz, ngrdcol, &
         163             :                                      rtm, rtp2, wp2, sqrt_wp2, &   ! In
         164             :                                      wprtp, w_1_n, w_2_n, mixt_frac, & ! In
         165             :                                      sigma_sqd_w, beta, rt_tol, &      ! In
         166             :                                      rt_1, rt_2, varnce_rt_1, &        ! Out
         167    17870112 :                                      varnce_rt_2, alpha_rt )           ! Out
         168             :
         169             :     ! Calculate the PDF component means and variances of thl.
         170             :     call ADG1_ADG2_responder_params( nz, ngrdcol, &
         171             :                                      thlm, thlp2, wp2, sqrt_wp2, &  ! In
         172             :                                      wpthlp, w_1_n, w_2_n, mixt_frac, & ! In
         173             :                                      sigma_sqd_w, beta, thl_tol, &      ! In
         174             :                                      thl_1, thl_2, varnce_thl_1, &      ! Out
         175    17870112 :                                      varnce_thl_2, alpha_thl )          ! Out
         176             :
         177             :     ! Calculate the PDF component means and variances of u wind.
         178             :     call ADG1_ADG2_responder_params( nz, ngrdcol, &
         179             :                                      um, up2, wp2, sqrt_wp2, &    ! In
         180             :                                      upwp, w_1_n, w_2_n, mixt_frac, & ! In
         181             :                                      sigma_sqd_w, beta, thl_tol, &    ! In
         182             :                                      u_1, u_2, varnce_u_1, &          ! Out
         183    17870112 :                                      varnce_u_2, alpha_u )            ! Out
         184             :
         185             :     ! Calculate the PDF component means and variances of v wind.
         186             :     call ADG1_ADG2_responder_params( nz, ngrdcol, &
         187             :                                      vm, vp2, wp2, sqrt_wp2, &    ! In
         188             :                                      vpwp, w_1_n, w_2_n, mixt_frac, & ! In
         189             :                                      sigma_sqd_w, beta, thl_tol, &    ! In
         190             :                                      v_1, v_2, varnce_v_1, &          ! Out
         191    17870112 :                                      varnce_v_2, alpha_v )            ! Out
         192             :
         193             :     ! Calculate the PDF component means and variances of passive scalars.
         194    17870112 :     if ( l_scalar_calc ) then
         195           0 :        do j = 1, sclr_dim
         196             :           call ADG1_ADG2_responder_params( nz, ngrdcol, &
         197             :                                            sclrm(:,:,j), sclrp2(:,:,j), & ! In
         198             :                                            wp2, sqrt_wp2, wpsclrp(:,:,j), & ! In
         199             :                                            w_1_n, w_2_n, mixt_frac,     & ! In
         200             :                                            sigma_sqd_w, beta,           & ! In
         201           0 :                                            sclr_tol(j),                 & ! In
         202             :                                            sclr_1(:,:,j), sclr_2(:,:,j),    & ! Out
         203             :                                            varnce_sclr_1(:,:,j),          & ! Out
         204             :                                            varnce_sclr_2(:,:,j),          & ! Out
         205           0 :                                            alpha_sclr(:,:,j) )              ! Out
         206             :        enddo ! i=1, sclr_dim
         207             :     endif ! l_scalar_calc
         208             :
         209             :     !$acc exit data delete( w_1_n, w_2_n )
         210             :
         211    17870112 :     return
         212             :
         213             :   end subroutine ADG1_pdf_driver
         214             :
         215             :   !=============================================================================
         216           0 :   subroutine ADG2_pdf_driver( nz, ngrdcol, sclr_dim, sclr_tol,          & ! In
         217           0 :                               wm, rtm, thlm, wp2, rtp2, thlp2,          & ! In
         218           0 :                               Skw, wprtp, wpthlp, sqrt_wp2, beta,       & ! In
         219           0 :                               sclrm, sclrp2, wpsclrp, l_scalar_calc,    & ! In
         220           0 :                               w_1, w_2,                                 & ! Out
         221           0 :                               rt_1, rt_2,                               & ! Out
         222           0 :                               thl_1, thl_2,                             & ! Out
         223           0 :                               varnce_w_1, varnce_w_2,                   & ! Out
         224           0 :                               varnce_rt_1, varnce_rt_2,                 & ! Out
         225           0 :                               varnce_thl_1, varnce_thl_2,               & ! Out
         226           0 :                               mixt_frac,                                & ! Out
         227           0 :                               alpha_rt, alpha_thl,                      & ! Out
         228           0 :                               sigma_sqd_w, sclr_1, sclr_2,              & ! Out
         229           0 :                               varnce_sclr_1, varnce_sclr_2, alpha_sclr )  ! Out
         230             :
         231             :     ! Calculates the PDF parameters using the Analytic Double Gaussian 2 (ADG2)
         232             :     ! PDF closure.
         233             :
         234             :     ! References:
         235             :     !-----------------------------------------------------------------------
         236             :
         237             :     use grid_class, only: &
         238             :         grid ! Type
         239             :
         240             :     use constants_clubb, only: &
         241             :         one,       & ! Constant(s)
         242             :         w_tol_sqd, &
         243             :         rt_tol,    &
         244             :         thl_tol
         245             :
         246             :     use clubb_precision, only: &
         247             :         core_rknd    ! Variable(s)
         248             :
         249             :     implicit none
         250             :
         251             :     integer, intent(in) :: &
         252             :       nz,           & ! Number of vertical levels
         253             :       ngrdcol,      & ! Number of grid columns
         254             :       sclr_dim        ! Number of passive scalars
         255             :
         256             :     real( kind = core_rknd ), intent(in), dimension(sclr_dim) :: &
         257             :       sclr_tol          ! Threshold(s) on the passive scalars  [units vary]
         258             :
         259             :     ! Input Variables
         260             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) ::  &
         261             :       wm,       & ! Mean of w-wind component (vertical velocity)  [m/s]
         262             :       rtm,      & ! Mean of total water mixing ratio              [kg/kg]
         263             :       thlm,     & ! Mean of liquid water potential temperature    [K]
         264             :       wp2,      & ! Variance of w (overall)                       [m^2/s^2]
         265             :       rtp2,     & ! Variance of r_t (overall)                     [(kg/kg)^2]
         266             :       thlp2,    & ! Variance of th_l (overall)                    [K^2]
         267             :       Skw,      & ! Skewness of w                                 [-]
         268             :       wprtp,    & ! Covariance of w and r_t                       [(kg/kg)(m/s)]
         269             :       wpthlp,   & ! Covariance of w and th_l                      [K(m/s)]
         270             :       sqrt_wp2    ! Square root of variance of w                  [m/s]
         271             :
         272             :     real ( kind = core_rknd ), dimension(ngrdcol), intent(in) :: &
         273             :       beta       ! CLUBB tunable parameter beta               [-]
         274             :
         275             :     real( kind = core_rknd ), dimension(ngrdcol,nz, sclr_dim), intent(in) ::  &
         276             :       sclrm,   & ! Mean of passive scalar (overall)           [units vary]
         277             :       sclrp2,  & ! Variance of passive scalar (overall)       [(units vary)^2]
         278             :       wpsclrp    ! Covariance of w and passive scalar         [m/s (units vary)]
         279             :
         280             :     logical, intent(in) :: &
         281             :       l_scalar_calc    ! Flag to perform calculations for passive scalars
         282             :
         283             :     ! Output Variables
         284             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) ::  &
         285             :       w_1,          & ! Mean of w (1st PDF component)                      [m/s]
         286             :       w_2,          & ! Mean of w (2nd PDF component)                      [m/s]
         287             :       rt_1,         & ! Mean of r_t (1st PDF component)                  [kg/kg]
         288             :       rt_2,         & ! Mean of r_t (2nd PDF component)                  [kg/kg]
         289             :       thl_1,        & ! Mean of th_l (1st PDF component)                     [K]
         290             :       thl_2,        & ! Mean of th_l (2nd PDF component)                     [K]
         291             :       varnce_w_1,   & ! Variance of w (1st PDF component)              [m^2/s^2]
         292             :       varnce_w_2,   & ! Variance of w (2nd PDF component)              [m^2/s^2]
         293             :       varnce_rt_1,  & ! Variance of r_t (1st PDF component)          [kg^2/kg^2]
         294             :       varnce_rt_2,  & ! Variance of r_t (2nd PDF component)          [kg^2/kg^2]
         295             :       varnce_thl_1, & ! Variance of th_l (1st PDF component)               [K^2]
         296             :       varnce_thl_2, & ! Variance of th_l (2nd PDF component)               [K^2]
         297             :       mixt_frac,    & ! Mixture fraction (weight of 1st PDF component)       [-]
         298             :       alpha_thl,    & ! Factor relating to normalized variance for th_l      [-]
         299             :       alpha_rt,     & ! Factor relating to normalized variance for r_t       [-]
         300             :       sigma_sqd_w     ! Width of individual w plumes (ADG1 parameter)        [-]
         301             :
         302             :     real( kind = core_rknd ), dimension(ngrdcol,nz, sclr_dim), intent(out) ::  &
         303             :       sclr_1,        & ! Mean of passive scalar (1st PDF component) [units vary]
         304             :       sclr_2,        & ! Mean of passive scalar (2nd PDF component) [units vary]
         305             :       varnce_sclr_1, & ! Variance of pass. sclr (1st PDF comp.) [(units vary)^2]
         306             :       varnce_sclr_2, & ! Variance of pass. sclr (2nd PDF comp.) [(units vary)^2]
         307             :       alpha_sclr       ! Factor relating to normalized variance for sclr     [-]
         308             :
         309             :     ! Local Variables
         310             :     real( kind = core_rknd ), dimension(ngrdcol,nz) ::  &
         311           0 :       w_1_n,         & ! Normalized mean of w (1st PDF component)            [-]
         312           0 :       w_2_n,         & ! Normalized mean of w (2nd PDF component)            [-]
         313           0 :       small_m_w,     & ! Luhar's m (tunable parameter)                       [-]
         314           0 :       big_m_w,       & ! Luhar's M (a function of Luhar's m)                 [-]
         315           0 :       sigma_sqd_w_1, & ! Normalized width parameter of w (1st PDF component) [-]
         316           0 :       sigma_sqd_w_2    ! Normalized width parameter of w (2nd PDF component) [-]
         317             :
         318             :     integer :: j, i  ! Loop index
         319             :
         320             :
         321             :     ! Calculate the mixture fraction and the PDF component means and variances
         322             :     ! of w.
         323             :
         324             :     ! Reproduce ADG2_w_closure using separate functions
         325           0 :     do i = 1, ngrdcol
         326             :
         327           0 :       call calc_Luhar_params( nz, Skw(i,:), wp2(i,:), &                      ! intent(in)
         328           0 :                               wp2(i,:), w_tol_sqd, &                         ! intent(in)
         329           0 :                               mixt_frac(i,:), big_m_w(i,:), small_m_w(i,:) ) ! intent(out)
         330             :
         331           0 :       call close_Luhar_pdf( nz, wm(i,:), wp2(i,:), mixt_frac(i,:),       & ! intent(in)
         332           0 :                             small_m_w(i,:), wp2(i,:), w_tol_sqd,         & ! intent(in)
         333           0 :                             sigma_sqd_w_1(i,:), sigma_sqd_w_2(i,:),      & ! intent(out)
         334           0 :                             varnce_w_1(i,:), varnce_w_2(i,:),            & ! intent(out)
         335           0 :                             w_1_n(i,:), w_2_n(i,:), w_1(i,:), w_2(i,:) )   ! intent(out)
         336             :
         337             :         ! Overwrite sigma_sqd_w for consistency with ADG1
         338           0 :         sigma_sqd_w(i,:) = min( one / ( one + small_m_w(i,:)**2 ), 0.99_core_rknd )
         339             :
         340             :     end do
         341             :
         342             :     ! Calculate the PDF component means and variances of rt.
         343             :     call ADG1_ADG2_responder_params( nz, ngrdcol, &                    ! In
         344             :                                      rtm, rtp2, wp2, sqrt_wp2, &       ! In
         345             :                                      wprtp, w_1_n, w_2_n, mixt_frac, & ! In
         346             :                                      sigma_sqd_w, beta, rt_tol, &      ! In
         347             :                                      rt_1, rt_2, varnce_rt_1, &        ! Out
         348           0 :                                      varnce_rt_2, alpha_rt )           ! Out
         349             :
         350             :     ! Calculate the PDF component means and variances of thl.
         351             :     call ADG1_ADG2_responder_params( nz, ngrdcol, &                     ! In
         352             :                                      thlm, thlp2, wp2, sqrt_wp2, &      ! In
         353             :                                      wpthlp, w_1_n, w_2_n, mixt_frac, & ! In
         354             :                                      sigma_sqd_w, beta, thl_tol, &      ! In
         355             :                                      thl_1, thl_2, varnce_thl_1, &      ! Out
         356           0 :                                      varnce_thl_2, alpha_thl )          ! Out
         357             :
         358             :     ! Calculate the PDF component means and variances of passive scalars.
         359           0 :     if ( l_scalar_calc ) then
         360           0 :        do j = 1, sclr_dim
         361             :           call ADG1_ADG2_responder_params( nz, ngrdcol,                   & ! In
         362             :                                            sclrm(:,:,j), sclrp2(:,:,j),   & ! In
         363             :                                            wp2, sqrt_wp2, wpsclrp(:,:,j), & ! In
         364             :                                            w_1_n, w_2_n, mixt_frac,       & ! In
         365             :                                            sigma_sqd_w, beta,             & ! In
         366           0 :                                            sclr_tol(j),                   & ! In
         367             :                                            sclr_1(:,:,j), sclr_2(:,:,j),  & ! Out
         368             :                                            varnce_sclr_1(:,:,j),          & ! Out
         369             :                                            varnce_sclr_2(:,:,j),          & ! Out
         370           0 :                                            alpha_sclr(:,:,j) )              ! Out
         371             :        enddo ! i=1, sclr_dim
         372             :     endif ! l_scalar_calc
         373             :
         374           0 :     return
         375             :
         376             :   end subroutine ADG2_pdf_driver
         377             :
         378             :   !=============================================================================
         379           0 :   subroutine Luhar_3D_pdf_driver( nz, wm, rtm, thlm, wp2, rtp2, thlp2,  & ! In
         380           0 :                                   Skw, Skrt, Skthl, wprtp, wpthlp,      & ! In
         381           0 :                                   w_1, w_2,                             & ! Out
         382           0 :                                   rt_1, rt_2,                           & ! Out
         383           0 :                                   thl_1, thl_2,                         & ! Out
         384           0 :                                   varnce_w_1, varnce_w_2,               & ! Out
         385           0 :                                   varnce_rt_1,  varnce_rt_2,            & ! Out
         386           0 :                                   varnce_thl_1, varnce_thl_2,           & ! Out
         387           0 :                                   mixt_frac )                             ! Out
         388             :
         389             :     ! Calculates the PDF parameters using the 3D Luhar PDF closure.
         390             :
         391             :     ! References:
         392             :     !-----------------------------------------------------------------------
         393             :
         394             :     use constants_clubb, only: &
         395             :         w_tol_sqd, & ! Constant(s)
         396             :         rt_tol,    &
         397             :         thl_tol
         398             :
         399             :     use clubb_precision, only: &
         400             :         core_rknd    ! Variable(s)
         401             :
         402             :     implicit none
         403             :
         404             :     integer, intent(in) :: &
         405             :       nz
         406             :
         407             :     ! Input Variables
         408             :     real( kind = core_rknd ), dimension(nz), intent(in) ::  &
         409             :       wm,     & ! Mean of w-wind component (vertical velocity)  [m/s]
         410             :       rtm,    & ! Mean of total water mixing ratio              [kg/kg]
         411             :       thlm,   & ! Mean of liquid water potential temperature    [K]
         412             :       wp2,    & ! Variance of w (overall)                       [m^2/s^2]
         413             :       rtp2,   & ! Variance of r_t (overall)                     [(kg/kg)^2]
         414             :       thlp2,  & ! Variance of th_l (overall)                    [K^2]
         415             :       Skw,    & ! Skewness of w                                 [-]
         416             :       Skrt,   & ! Skewness of rt                                [-]
         417             :       Skthl,  & ! Skewness of th_l                              [-]
         418             :       wprtp,  & ! Covariance of w and r_t                       [(kg/kg)(m/s)]
         419             :       wpthlp    ! Covariance of w and th_l                      [K(m/s)]
         420             :
         421             :     ! Output Variables
         422             :     real( kind = core_rknd ), dimension(nz), intent(out) ::  &
         423             :       w_1,          & ! Mean of w (1st PDF component)                      [m/s]
         424             :       w_2,          & ! Mean of w (2nd PDF component)                      [m/s]
         425             :       rt_1,         & ! Mean of r_t (1st PDF component)                  [kg/kg]
         426             :       rt_2,         & ! Mean of r_t (2nd PDF component)                  [kg/kg]
         427             :       thl_1,        & ! Mean of th_l (1st PDF component)                     [K]
         428             :       thl_2,        & ! Mean of th_l (2nd PDF component)                     [K]
         429             :       varnce_w_1,   & ! Variance of w (1st PDF component)              [m^2/s^2]
         430             :       varnce_w_2,   & ! Variance of w (2nd PDF component)              [m^2/s^2]
         431             :       varnce_rt_1,  & ! Variance of r_t (1st PDF component)          [kg^2/kg^2]
         432             :       varnce_rt_2,  & ! Variance of r_t (2nd PDF component)          [kg^2/kg^2]
         433             :       varnce_thl_1, & ! Variance of th_l (1st PDF component)               [K^2]
         434             :       varnce_thl_2, & ! Variance of th_l (2nd PDF component)               [K^2]
         435             :       mixt_frac       ! Mixture fraction (weight of 1st PDF component)       [-]
         436             :
         437             :     real( kind = core_rknd ), dimension(nz) ::  &
         438           0 :       w_1_n,           & ! Normalized mean of w (1st PDF component)          [-]
         439           0 :       w_2_n,           & ! Normalized mean of w (2nd PDF component)          [-]
         440           0 :       rt_1_n,          & ! Normalized mean of rt (1st PDF component)         [-]
         441           0 :       rt_2_n,          & ! Normalized mean of rt (2nd PDF component)         [-]
         442           0 :       thl_1_n,         & ! Normalized mean of thl (1st PDF component)        [-]
         443           0 :       thl_2_n,         & ! Normalized mean of thl (2nd PDF component)        [-]
         444           0 :       small_m_w,       & ! Luhar's m (tunable parameter) for w               [-]
         445           0 :       big_m_w,         & ! Luhar's M (a function of Luhar's m) for w         [-]
         446           0 :       small_m_rt,      & ! Luhar's m (tunable parameter) for rt              [-]
         447           0 :       big_m_rt,        & ! Luhar's M (a function of Luhar's m) for rt        [-]
         448           0 :       small_m_thl,     & ! Luhar's m (tunable parameter) for thl             [-]
         449           0 :       big_m_thl,       & ! Luhar's M (a function of Luhar's m) for thl       [-]
         450           0 :       sigma_sqd_w_1,   & ! Normalized width parameter of w (1st PDF comp.)   [-]
         451           0 :       sigma_sqd_w_2,   & ! Normalized width parameter of w (2nd PDF comp.)   [-]
         452           0 :       sigma_sqd_rt_1,  & ! Normalized width parameter of rt (1st PDF comp.)  [-]
         453           0 :       sigma_sqd_rt_2,  & ! Normalized width parameter of rt (2nd PDF comp.)  [-]
         454           0 :       sigma_sqd_thl_1, & ! Normalized width parameter of thl (1st PDF comp.) [-]
         455           0 :       sigma_sqd_thl_2    ! Normalized width parameter of thl (2nd PDF comp.) [-]
         456             :
         457             :     integer :: k    ! Vertical loop index
         458             :
         459             :
         460           0 :     do k = 1, nz, 1
         461             :
         462           0 :        if ( ( abs( Skw(k) ) >= abs( Skthl(k) ) ) &
         463           0 :             .and. ( abs( Skw(k) ) >= abs( Skrt(k) ) ) ) then
         464             :
         465             :           ! w has the greatest magnitude of skewness.
         466             :
         467             :           ! Solve for the w PDF
         468             :           call calc_Luhar_params( 1, Skw(k), wp2(k), &                     ! In
         469             :                                   wp2(k), w_tol_sqd, &                     ! In
         470           0 :                                   mixt_frac(k), big_m_w(k), small_m_w(k) ) ! Out
         471             :
         472             :           call close_Luhar_pdf( 1, wm(k), wp2(k), mixt_frac(k),     & ! In
         473             :                                 small_m_w(k), wp2(k), w_tol_sqd,    & ! In
         474             :                                 sigma_sqd_w_1(k), sigma_sqd_w_2(k), & ! Out
         475             :                                 varnce_w_1(k), varnce_w_2(k),       & ! Out
         476           0 :                                 w_1_n(k), w_2_n(k), w_1(k), w_2(k)  ) ! Out
         477             :
         478             :           ! Solve for the thl PDF
         479             :           call backsolve_Luhar_params( Skw(k), Skthl(k),            & ! In
         480             :                                        big_m_w(k), mixt_frac(k),    & ! In
         481             :                                        big_m_thl(k), small_m_thl(k) ) ! Out
         482             :
         483             :           call close_Luhar_pdf( 1, thlm(k), thlp2(k), mixt_frac(k),     & ! In
         484             :                                 small_m_thl(k), wpthlp(k), thl_tol**2,  & ! In
         485             :                                 sigma_sqd_thl_1(k), sigma_sqd_thl_2(k), & ! Out
         486             :                                 varnce_thl_1(k), varnce_thl_2(k),       & ! Out
         487             :                                 thl_1_n(k), thl_2_n(k),                 & ! Out
         488           0 :                                 thl_1(k), thl_2(k)                      ) ! Out
         489             :
         490             :           ! Solve for the rt PDF
         491             :           call backsolve_Luhar_params( Skw(k), Skrt(k),           & ! In
         492             :                                        big_m_w(k), mixt_frac(k),  & ! In
         493             :                                        big_m_rt(k), small_m_rt(k) ) ! Out
         494             :
         495             :           call close_Luhar_pdf( 1, rtm(k), rtp2(k), mixt_frac(k),      & ! In
         496             :                                 small_m_rt(k), wprtp(k), rt_tol**2,    & ! In
         497             :                                 sigma_sqd_rt_1(k), sigma_sqd_rt_2(k),  & ! Out
         498             :                                 varnce_rt_1(k), varnce_rt_2(k),        & ! Out
         499           0 :                                 rt_1_n(k), rt_2_n(k), rt_1(k), rt_2(k) ) ! Out
         500             :
         501             :        elseif ( ( abs( Skthl(k) ) > abs( Skw(k) ) ) &
         502           0 :                   .and. ( abs( Skthl(k) ) >= abs( Skrt(k) ) ) ) then
         503             :
         504             :           ! theta-l has the greatest magnitude of skewness.
         505             :
         506             :           ! Solve for the thl PDF
         507             :           call calc_Luhar_params( 1, Skthl(k), wpthlp(k),     & ! In
         508             :                                   thlp2(k), thl_tol**2,       & ! In
         509             :                                   mixt_frac(k), big_m_thl(k), & ! Out
         510           0 :                                   small_m_thl(k)              ) ! Out
         511             :
         512             :           ! Solve for the thl PDF
         513             :           call close_Luhar_pdf( 1, thlm(k), thlp2(k), mixt_frac(k),     & ! In
         514             :                                 small_m_thl(k), wpthlp(k), thl_tol**2,  & ! In
         515             :                                 sigma_sqd_thl_1(k), sigma_sqd_thl_2(k), & ! Out
         516             :                                 varnce_thl_1(k), varnce_thl_2(k),       & ! Out
         517             :                                 thl_1_n(k), thl_2_n(k),                 & ! Out
         518           0 :                                 thl_1(k), thl_2(k)                      ) ! Out
         519             :
         520             :           ! Solve for the w PDF
         521             :           call backsolve_Luhar_params( Skthl(k), Skw(k),           & ! In
         522             :                                        big_m_thl(k), mixt_frac(k), & ! In
         523             :                                        big_m_w(k), small_m_w(k) )    ! Out
         524             :
         525             :           call close_Luhar_pdf( 1, wm(k), wp2(k), mixt_frac(k),     & ! In
         526             :                                 small_m_w(k), wp2(k), w_tol_sqd,    & ! In
         527             :                                 sigma_sqd_w_1(k), sigma_sqd_w_2(k), & ! Out
         528             :                                 varnce_w_1(k), varnce_w_2(k),       & ! Out
         529           0 :                                 w_1_n(k), w_2_n(k), w_1(k), w_2(k)  ) ! Out
         530             :
         531             :           ! Solve for the rt PDF
         532             :           call backsolve_Luhar_params( Skthl(k), Skrt(k),           & ! In
         533             :                                        big_m_thl(k), mixt_frac(k),  & ! In
         534             :                                        big_m_rt(k), small_m_rt(k) )   ! Out
         535             :
         536             :           call close_Luhar_pdf( 1, rtm(k), rtp2(k), mixt_frac(k),      & ! In
         537             :                                 small_m_rt(k), wprtp(k), rt_tol**2,    & ! In
         538             :                                 sigma_sqd_rt_1(k), sigma_sqd_rt_2(k),  & ! Out
         539             :                                 varnce_rt_1(k), varnce_rt_2(k),        & ! Out
         540           0 :                                 rt_1_n(k), rt_2_n(k), rt_1(k), rt_2(k) ) ! Out
         541             :
         542             :        else
         543             :
         544             :           ! rt has the greatest magnitude of skewness.
         545             :
         546             :           ! Solve for the rt PDF
         547             :           call calc_Luhar_params( 1, Skrt(k), wprtp(k),      & ! In
         548             :                                   rtp2(k), rt_tol**2,        & ! In
         549             :                                   mixt_frac(k), big_m_rt(k), & ! Out
         550           0 :                                   small_m_rt(k)              ) ! Out
         551             :
         552             :           ! Solve for the rt PDF
         553             :           call close_Luhar_pdf( 1, rtm(k), rtp2(k), mixt_frac(k),      & ! In
         554             :                                 small_m_rt(k), wprtp(k), rt_tol**2,    & ! In
         555             :                                 sigma_sqd_rt_1(k), sigma_sqd_rt_2(k),  & ! Out
         556             :                                 varnce_rt_1(k), varnce_rt_2(k),        & ! Out
         557           0 :                                 rt_1_n(k), rt_2_n(k), rt_1(k), rt_2(k) ) ! Out
         558             :
         559             :           ! Solve for the w PDF
         560             :           call backsolve_Luhar_params( Skrt(k), Skw(k),           & ! In
         561             :                                        big_m_rt(k), mixt_frac(k), & ! In
         562             :                                        big_m_w(k), small_m_w(k) )   ! Out
         563             :
         564             :           call close_Luhar_pdf( 1, wm(k), wp2(k), mixt_frac(k),     & ! In
         565             :                                 small_m_w(k), wp2(k), w_tol_sqd,    & ! In
         566             :                                 sigma_sqd_w_1(k), sigma_sqd_w_2(k), & ! Out
         567             :                                 varnce_w_1(k), varnce_w_2(k),       & ! Out
         568           0 :                                 w_1_n(k), w_2_n(k), w_1(k), w_2(k)  ) ! OUt
         569             :
         570             :           ! Solve for the thl PDF
         571             :           call backsolve_Luhar_params( Skrt(k), Skthl(k),           & ! In
         572             :                                        big_m_rt(k), mixt_frac(k),   & ! In
         573             :                                        big_m_thl(k), small_m_thl(k) ) ! Out
         574             :
         575             :           call close_Luhar_pdf( 1, thlm(k), thlp2(k), mixt_frac(k),     & ! In
         576             :                                 small_m_thl(k), wpthlp(k), thl_tol**2,  & ! In
         577             :                                 sigma_sqd_thl_1(k), sigma_sqd_thl_2(k), & ! Out
         578             :                                 varnce_thl_1(k), varnce_thl_2(k),       & ! Out
         579             :                                 thl_1_n(k), thl_2_n(k),                 & ! Out
         580           0 :                                 thl_1(k), thl_2(k)                      ) ! Out
         581             :
         582             :        endif
         583             :
         584             :     enddo ! k = 1, nz, 1
         585             :
         586             :
         587           0 :     return
         588             :
         589             :   end subroutine Luhar_3D_pdf_driver
         590             :
         591             :   !======================================================================
         592    17870112 :   subroutine ADG1_w_closure( nz, ngrdcol, wm, wp2, Skw, sigma_sqd_w, & !In
         593    17870112 :                              sqrt_wp2, mixt_frac_max_mag, &            !In
         594    17870112 :                              w_1, w_2, w_1_n, w_2_n, &                 !Out
         595    17870112 :                              varnce_w_1, varnce_w_2, mixt_frac )       !Out
         596             :
         597             :
         598             :     ! Description:
         599             :     ! Calculates the mixture fraction, the PDF component means of w, and the PDF
         600             :     ! component variances of w for the Analytic Double Gaussian 1 (ADG1)
         601             :     ! closure.  It sets the widths of both w Gaussians to be the same, and
         602             :     ! furthermore, relates them to the overall variance of w (<w'^2>) by a
         603             :     ! parameter, sigma_sqd_w.  The equation is:
         604             :     !
         605             :     ! sigma_w_1^2 = sigma_w_2^2 = sigma_sqd_w * <w'^2>;
         606             :     !
         607             :     ! where sigma_w_1^2 is the variance of w in the 1st PDF component,
         608             :     ! sigma_w_2^2 is the variance of w in the 2nd PDF component, and parameter
         609             :     ! sigma_sqd_w must have a value within the range 0 <= sigma_sqd_w <= 1.
         610             :     !
         611             :     ! References:
         612             :     ! Golaz, J-C., V. E. Larson, and W. R. Cotton, 2002a: A PDF-based model for
         613             :     ! boundary layer clouds. Part I: Method and model description. J. Atmos.
         614             :     ! Sci., 59, 3540–3551.
         615             :     !
         616             :     ! Vincent E. Larson and Jean-Christophe Golaz, 2005: Using Probability
         617             :     ! Density Functions to Derive Consistent Closure Relationships among
         618             :     ! Higher-Order Moments. Mon. Wea. Rev., 133, 1023–1042.
         619             :     !-----------------------------------------------------------------------
         620             :
         621             :     use constants_clubb, only: &
         622             :         four,      & ! Constant(s)
         623             :         one,       &
         624             :         one_half,  &
         625             :         zero,      &
         626             :         w_tol_sqd
         627             :
         628             :     use clubb_precision, only: &
         629             :         core_rknd     ! Precision
         630             :
         631             :     implicit none
         632             :
         633             :     integer, intent(in) :: &
         634             :       ngrdcol,  & ! Number of grid columns
         635             :       nz          ! Number of vertical level
         636             :
         637             :     ! Input Variables
         638             :     real( kind = core_rknd ), dimension(ngrdcol,nz),  intent(in) :: &
         639             :       wm,                & ! Mean of w (overall)                   [m/s]
         640             :       wp2,               & ! Variance of w (overall)               [m^2/s^2]
         641             :       Skw,               & ! Skewness of w                         [-]
         642             :       sigma_sqd_w,       & ! Widths of each w Gaussian             [-]
         643             :       sqrt_wp2             ! Square root of the variance of w      [m/s]
         644             :
         645             :     real( kind = core_rknd ),  intent(in) :: &
         646             :       mixt_frac_max_mag    ! Maximum allowable mag. of mixt_frac   [-]
         647             :
         648             :     ! Output Variables
         649             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
         650             :       w_1,        & ! Mean of w (1st PDF component)                [m/s]
         651             :       w_2,        & ! Mean of w (2nd PDF component)                [m/s]
         652             :       w_1_n,      & ! Normalized mean of w (1st PDF component)     [-]
         653             :       w_2_n,      & ! Normalized mean of w (2nd PDF component)     [-]
         654             :       varnce_w_1, & ! Variance of w (1st PDF component)            [m^2/s^2]
         655             :       varnce_w_2, & ! Variance of w (2nd PDF component)            [m^2/s^2]
         656             :       mixt_frac     ! Mixture fraction                             [-]
         657             :
         658             :     ! Local Variable
         659             :     integer :: k, i
         660             :
         661             :     !----- Begin Code -----
         662             :
         663             :     !$acc parallel loop gang vector collapse(2) default(present)
         664  1536829632 :     do k = 1, nz
         665 25380961632 :        do i = 1, ngrdcol
         666             :
         667 25363091520 :           if ( wp2(i,k) > w_tol_sqd ) then
         668             :
         669             :              ! Width (standard deviation) parameters are non-zero
         670             :
         671             :              ! The variable "mixt_frac" is the weight of the 1st PDF component.  The
         672             :              ! weight of the 2nd PDF component is "1-mixt_frac".  If there isn't any
         673             :              ! skewness of w (Sk_w = 0 because w'^3 = 0), mixt_frac = 0.5, and both
         674             :              ! PDF components are equally weighted.  If there is positive skewness of
         675             :              ! w (Sk_w > 0 because w'^3 > 0), 0 < mixt_frac < 0.5, and the 2nd PDF
         676             :              ! component has greater weight than does the 1st PDF component.  If there
         677             :              ! is negative skewness of w (Sk_w < 0 because w'^3 < 0),
         678             :              ! 0.5 < mixt_frac < 1, and the 1st PDF component has greater weight than
         679             :              ! does the 2nd PDF component.
         680  5993468723 :              if ( abs( Skw(i,k) ) <= 1.0e-5_core_rknd ) then
         681   442112203 :                 mixt_frac(i,k) = one_half
         682             :              else
         683             :                 mixt_frac(i,k) &
         684             :                  = one_half &
         685  5551356520 :                      * ( one - Skw(i,k) / sqrt( four * ( one - sigma_sqd_w(i,k) )**3 + Skw(i,k)**2 ) )
         686             :              endif
         687             :
         688             :              ! Clip mixt_frac, and 1 - mixt_frac, to avoid dividing by a small number.
         689             :              ! Formula for mixt_frac_max_mag =
         690             :              ! 1 - ( 1/2 * ( 1 - Skw_max
         691             :              !                   / sqrt( 4*( 1 - sigma_sqd_w )^3 + Skw_max^2 ) ) ),
         692             :              ! where sigma_sqd_w is fixed at 0.4 for this calculation.
         693             :              mixt_frac(i,k) = min( max( mixt_frac(i,k), one - mixt_frac_max_mag ), &
         694  5993468723 :                                    mixt_frac_max_mag )
         695             :
         696             :              ! The normalized mean of w for Gaussian "plume" 1 is w_1_n.  It's value
         697             :              ! will always be greater than 0.  As an example, a value of 1.0 would
         698             :              ! indicate that the actual mean of w for Gaussian "plume" 1 is found 1.0
         699             :              ! standard deviation above the overall mean for w.
         700             :              w_1_n(i,k) &
         701  5993468723 :               = sqrt( ( ( one - mixt_frac(i,k) ) / mixt_frac(i,k) ) * ( one - sigma_sqd_w(i,k) ) )
         702             :              ! The normalized mean of w for Gaussian "plume" 2 is w_2_n.  It's value
         703             :              ! will always be less than 0.  As an example, a value of -0.5 would
         704             :              ! indicate that the actual mean of w for Gaussian "plume" 2 is found 0.5
         705             :              ! standard deviations below the overall mean for w.
         706             :              w_2_n(i,k) &
         707  5993468723 :               = -sqrt( ( mixt_frac(i,k) / ( one - mixt_frac(i,k) ) ) * ( one - sigma_sqd_w(i,k) ) )
         708             :              ! The mean of w for Gaussian "plume" 1 is w_1.
         709  5993468723 :              w_1(i,k) = wm(i,k) + sqrt_wp2(i,k) * w_1_n(i,k)
         710             :              ! The mean of w for Gaussian "plume" 2 is w_2.
         711  5993468723 :              w_2(i,k) = wm(i,k) + sqrt_wp2(i,k) * w_2_n(i,k)
         712             :
         713             :              ! The variance of w for Gaussian "plume" 1 for varnce_w_1.
         714  5993468723 :              varnce_w_1(i,k) = sigma_sqd_w(i,k) * wp2(i,k)
         715             :              ! The variance of w for Gaussian "plume" 2 for varnce_w_2.
         716             :              ! The variance in both Gaussian "plumes" is defined to be the same.
         717  5993468723 :              varnce_w_2(i,k) = sigma_sqd_w(i,k) * wp2(i,k)
         718             :
         719             :           else
         720             :
         721             :              ! Vertical velocity doesn't vary.
         722 17850663277 :              mixt_frac(i,k)  = one_half
         723 17850663277 :              w_1_n(i,k)      = sqrt( one - sigma_sqd_w(i,k) )
         724 17850663277 :              w_2_n(i,k)      = -sqrt( one - sigma_sqd_w(i,k) )
         725 17850663277 :              w_1(i,k)        = wm(i,k)
         726 17850663277 :              w_2(i,k)        = wm(i,k)
         727 17850663277 :              varnce_w_1(i,k) = zero
         728 17850663277 :              varnce_w_2(i,k) = zero
         729             :
         730             :           endif  ! Widths non-zero
         731             :        end do
         732             :     end do
         733             :     !$acc end parallel loop
         734    17870112 :     return
         735             :
         736             :
         737             :
         738             :   end subroutine ADG1_w_closure
         739             :
         740             :
         741             :   !=============================================================================
         742           0 :   subroutine calc_Luhar_params( nz, Skx, wpxp,            & ! In
         743           0 :                                 xp2, x_tol_sqd,           & ! In
         744           0 :                                 mixt_frac, big_m, small_m ) ! Out
         745             :
         746             :     ! Description:
         747             :     ! For the Luhar closure, this subroutine takes Skx (and Skw) as input and
         748             :     ! outputs the mixture fraction, big_m, and small_m. This code was written
         749             :     ! using the equations and nomenclature of Larson et al. (2002) Appendix
         750             :     ! section e.
         751             :     !
         752             :     ! The relationship between skewness of x (Skx), mixture fraction (a), and
         753             :     ! Luhar's small m (m) is given by:
         754             :     !
         755             :     ! Skx^2 = ( m^2 * ( m^2 + 3 )^2 / ( m^2 + 1 )^3 )
         756             :     !         * ( 1 - 2*a )^2 / ( a * ( 1 - a ) ).
         757             :     !
         758             :     ! Luhar's large M (M) is used to more easily express the factor involving
         759             :     ! the m's:
         760             :     !
         761             :     ! M = ( m^2 + 1 )^3 / ( m^2 * ( m^2 + 3 )^2 ).
         762             :     !
         763             :     ! The equation involving skewness of x becomes:
         764             :     !
         765             :     ! Skx^2 = ( 1 / M ) * ( 1 - 2*a )^2 / ( a * ( 1 - a ) );
         766             :     !
         767             :     ! or:
         768             :     !
         769             :     ! M * Skx^2 = ( 1 - 2*a )^2 / ( a * ( 1 - a ) ).
         770             :     !
         771             :     ! This equation can be rewritten as:
         772             :     !
         773             :     ! ( a * ( 1 - a ) ) * M * Skx^2 = ( 1 - 2*a )^2;
         774             :     !
         775             :     ! as well as:
         776             :     !
         777             :     ! ( a - a^2 ) * M * Skx^2 = 1 - 4*a + 4*a^2;
         778             :     !
         779             :     ! and eventually as:
         780             :     !
         781             :     ! ( 4 + M * Skx^2 ) * a^2 - ( 4 + M * Skx^2 ) * a + 1 = 0.
         782             :     !
         783             :     ! Solving the quadratic equation for a:
         784             :     !
         785             :     ! a = (1/2) * ( 1 +- Skx * sqrt( 1 / ( 4/M + Skx^2 ) ) ).
         786             :     !
         787             :     ! Since by definition, mu_w_1 >= mu_w_2, a < 0.5 when Skw > 0, the equation
         788             :     ! for mixture fraction is:
         789             :     !
         790             :     ! a = (1/2) * ( 1 - Skx * sqrt( 1 / ( 4/M + Skx^2 ) ) ).
         791             :     !
         792             :     ! For 3-D Luhar, the variable (w, rt, or theta-l) with the greatest
         793             :     ! magnitude of skewness is used to calculate mixture fraction.  Since it is
         794             :     ! desirable to still have a < 0.5 when Skw > 0 and a > 0.5 when Skw < 0, the
         795             :     ! sign function is used.  The value of Skx is replaced by:
         796             :     !
         797             :     ! Skx|_adj = sgn( <w'x'> ) * Skx;
         798             :     !
         799             :     ! where
         800             :     !
         801             :     ! sgn( <w'x'> ) = | 1 when <w'x'> >= 0
         802             :     !                 | -1 when <w'x'> < 0.
         803             :     !
         804             :     ! Since Skx|_adj^2 = ( sgn( <w'x'> ) * Skx )^2 = ( sgn( <w'x'> ) )^2 * Skx^2
         805             :     ! = Skx^2, the equation for mixture fraction is:
         806             :     !
         807             :     ! a = (1/2) * ( 1 - sgn( <w'x'> ) * Skx * sqrt( 1 / ( 4/M + Skx^2 ) ) ).
         808             :     !
         809             :     ! When using the ADG2 closure or when using the 3-D Luhar closure when the
         810             :     ! variable with the greatest magnitude of skewness is w, Skw = Skx and
         811             :     ! sgn( <w'^2> ) is always equal to 1, reducing the equation to its previous
         812             :     ! form.
         813             :
         814             :     ! References:
         815             :     ! Vincent E. Larson, Jean-Christophe Golaz, and William R. Cotton, 2002:
         816             :     ! Small-Scale and Mesoscale Variability in Cloudy Boundary Layers: Joint
         817             :     ! Probability Density Functions. J. Atmos. Sci., 59, 3519–3539.
         818             :     !-----------------------------------------------------------------------
         819             :
         820             :     use constants_clubb, only: &
         821             :         four,       & ! Constant(s)
         822             :         three,      &
         823             :         one,        &
         824             :         two_thirds, &
         825             :         one_half,   &
         826             :         one_third,  &
         827             :         zero
         828             :
         829             :     use clubb_precision, only: &
         830             :         core_rknd     ! Precision
         831             :
         832             :     implicit none
         833             :
         834             :     ! Input Variables
         835             :     integer, intent(in) :: &
         836             :       nz    ! Number of vertical levels
         837             :
         838             :     real( kind = core_rknd ), dimension(nz), intent(in) :: &
         839             :       Skx,  & ! Skewness of x ( <x'^3> / <x'2>^(3/2) )      [-]
         840             :       wpxp, & ! Covariance of w and x                       [m/s (x units)]
         841             :       xp2     ! Variance of x                               [(x units)^2]
         842             :
         843             :     real( kind = core_rknd ), intent(in) :: &
         844             :       x_tol_sqd    ! Tolerance value of x squared           [(x units)^2]
         845             :
         846             :     ! Output Variables
         847             :     real( kind = core_rknd ), dimension(nz), intent(out) :: &
         848             :       mixt_frac, & ! Mixture fraction                       [-]
         849             :       big_m,     & ! Luhar's M                              [-]
         850             :       small_m      ! Luhar's m                              [-]
         851             :
         852             :     ! Local Variables
         853             :     real( kind = core_rknd ), dimension(nz) :: &
         854           0 :       small_m_sqd, & ! Luhar's m^2                                           [-]
         855           0 :       sgn_wpxp       ! Sgn(<w'x'>); 1 when <w'x'> >= 0 or -1 when <w'x'> < 0 [-]
         856             :
         857             :
         858           0 :     where ( xp2 > x_tol_sqd )
         859             :
         860             :        ! Width (standard deviation) parameters are non-zero
         861             :
         862             :        ! Calculate Luhar's m (small m).
         863             :        ! If Skx is very small, then small_m will tend to zero which risks
         864             :        ! divide-by-zero.  To ameliorate this problem, we enforce abs( x_1_n )
         865             :        ! and abs( x_2_n ) > 0.05.
         866             :        ! Note:  Luhar's small_m (m) is the only tunable parameter in the Luhar
         867             :        !        closure, so this equation can be changed.  However, the value
         868             :        !        of m should go toward 0 as Skx goes toward 0 so that the double
         869             :        !        Gaussian reduces to a single Gaussian when the distribution is
         870             :        !        unskewed.
         871             :        small_m = max( two_thirds * abs( Skx )**one_third, 0.05_core_rknd )
         872             :
         873             :        ! Calculate m^2.
         874             :        small_m_sqd = small_m**2
         875             :
         876             :        ! Calculate Luhar's M (big M).
         877             :        big_m = ( one + small_m_sqd )**3 &
         878             :                / ( ( three + small_m_sqd )**2 * small_m_sqd )
         879             :
         880             :        ! Calculate sgn( <w'x'> ).
         881             :        where ( wpxp >= zero )
         882             :           sgn_wpxp = one
         883             :        elsewhere ! <w'x'> < 0
         884             :           sgn_wpxp = -one
         885             :        endwhere ! <w'x'> >= 0
         886             :
         887             :        ! Calculate mixture fraction.
         888             :        mixt_frac = one_half &
         889             :                    * ( one - sgn_wpxp * Skx &
         890             :                              * sqrt( one / ( ( four / big_m ) + Skx**2 ) ) )
         891             :
         892             :     elsewhere
         893             :
         894             :        ! Variable x doesn't vary.
         895             :        mixt_frac = one_half
         896             :
         897             :        ! For output purposes, set small_m and big_m to 0.
         898             :        small_m = zero
         899             :        big_m   = zero
         900             :
         901             :     endwhere  ! Widths non-zero
         902             :
         903             :
         904           0 :     return
         905             :
         906             :   end subroutine calc_Luhar_params
         907             :
         908             :   !=============================================================================
         909           0 :   subroutine close_Luhar_pdf( nz, xm, xp2, mixt_frac,       & ! In
         910           0 :                               small_m, wpxp, x_tol_sqd,     & ! In
         911           0 :                               sigma_sqd_x_1, sigma_sqd_x_2, & ! Out
         912           0 :                               varnce_x_1, varnce_x_2,       & ! Out
         913           0 :                               x_1_n, x_2_n, x_1, x_2 )        ! Out
         914             :
         915             :     ! Description:
         916             :     ! For the Luhar closure, this subroutine takes Skx, xm, xp2, and mixt_frac,
         917             :     ! big_m, and small_m (calculated in calc_Luhar_params) as input and outputs
         918             :     ! the PDF component means and variances of a variable x in the joint-PDF
         919             :     ! according to Luhar et al. (1996).  This code was written using the
         920             :     ! equations and nomenclature of Larson et al. (2002) Appendix section e.
         921             :
         922             :     ! References:
         923             :     !    Vincent E. Larson, Jean-Christophe Golaz, and William R. Cotton, 2002:
         924             :     !    Small-Scale and Mesoscale Variability in Cloudy Boundary Layers: Joint
         925             :     !    Probability Density Functions. J. Atmos. Sci., 59, 3519–3539.
         926             :     !-----------------------------------------------------------------------
         927             :
         928             :     use constants_clubb, only: &
         929             :         one,  & ! Constant(s)
         930             :         zero
         931             :
         932             :     use clubb_precision, only: &
         933             :         core_rknd     ! Precision
         934             :
         935             :     implicit none
         936             :
         937             :     ! Input Variables
         938             :     integer, intent(in) :: &
         939             :       nz    ! Number of vertical levels
         940             :
         941             :     real( kind = core_rknd ), dimension(nz), intent(in) :: &
         942             :       xm,        & ! Mean (overall) of x, <x>                    [(x units)]
         943             :       xp2,       & ! Variance (overall) of x, <x'^2>             [(x units)^2]
         944             :       mixt_frac, & ! Mixture fraction                            [-]
         945             :       small_m,   & ! Luhar's small m                             [-]
         946             :       wpxp         ! Covariance of w and x, <w'x'>               [m/s (x units)]
         947             :
         948             :     real( kind = core_rknd ), intent(in) :: &
         949             :       x_tol_sqd    ! Tolerance value of x squared                [(x units)^2]
         950             :
         951             :     ! Output Variables
         952             :     real( kind = core_rknd ), dimension(nz), intent(out) :: &
         953             :       sigma_sqd_x_1, & ! Normalized width parameter of x (1st PDF component) [-]
         954             :       sigma_sqd_x_2, & ! Normalized width parameter of x (1st PDF component) [-]
         955             :       varnce_x_1,    & ! Variance of x (1st PDF component)         [(x units)^2]
         956             :       varnce_x_2,    & ! Variance of x (2nd PDF component)         [(x units)^2]
         957             :       x_1_n,         & ! Normalized mean of x (1st PDF component)            [-]
         958             :       x_2_n,         & ! Normalized mean of x (2nd PDF component)            [-]
         959             :       x_1,           & ! Mean of x (1st PDF component)               [(x units)]
         960             :       x_2              ! Mean of x (2nd PDF component)               [(x units)]
         961             :
         962             :     ! Local Variables
         963             :     real( kind = core_rknd), dimension(nz) :: &
         964           0 :       sqrt_xp2, & ! Square root of the variance of x                 [(x units)]
         965           0 :       sgn_wpxp    ! Sgn( <w'x'> ); 1 when <w'x'> >= 0 or -1 when <w'x'> < 0  [-]
         966             :
         967             :
         968           0 :     where ( xp2 > x_tol_sqd )
         969             :
         970             :        ! Width (standard deviation) parameters are non-zero
         971             :
         972             :        ! Calculate sgn( <w'x'> ).
         973             :        where ( wpxp >= zero )
         974             :           sgn_wpxp = one
         975             :        elsewhere ! <w'x'> < 0
         976             :           sgn_wpxp = -one
         977             :        endwhere ! <w'x'> >= 0
         978             :
         979             :        ! Calculate the square root of the overall variance of x.
         980             :        sqrt_xp2 = sqrt( xp2 )
         981             :
         982             :        ! Normalized width parameter of x in the 1st PDF component.
         983             :        sigma_sqd_x_1 &
         984             :        = ( one - mixt_frac ) / ( mixt_frac * ( one + small_m**2 ) )
         985             :
         986             :        ! The variance of x in the 1st PDF component.
         987             :        varnce_x_1 = sigma_sqd_x_1 * xp2
         988             :
         989             :        ! Normalized width parameter of x in the 2nd PDF component.
         990             :        sigma_sqd_x_2 &
         991             :        = mixt_frac / ( ( one - mixt_frac ) * ( one + small_m**2 ) )
         992             :
         993             :        ! The variance of x in the 2nd PDF component.
         994             :        varnce_x_2 = sigma_sqd_x_2 * xp2
         995             :
         996             :        ! Normalized mean of x in the 1st PDF component.
         997             :        x_1_n = sgn_wpxp * small_m * sqrt( sigma_sqd_x_1 )
         998             :
         999             :        ! Normalized mean of x in the 2nd PDF component.
        1000             :        x_2_n = -sgn_wpxp * small_m * sqrt( sigma_sqd_x_2 )
        1001             :
        1002             :        ! The mean of x in the 1st PDF component.
        1003             :        x_1 = xm + sqrt_xp2 * x_1_n
        1004             :
        1005             :        ! The mean of x in the 2nd PDF component.
        1006             :        x_2 = xm + sqrt_xp2 * x_2_n
        1007             :
        1008             :     elsewhere
        1009             :
        1010             :        ! Variable x doesn't vary.
        1011             :        sigma_sqd_x_1 = ( one / ( one + small_m**2 ) )
        1012             :        sigma_sqd_x_2 = ( one / ( one + small_m**2 ) )
        1013             :        varnce_x_1    = zero
        1014             :        varnce_x_2    = zero
        1015             :        x_1_n         = sgn_wpxp * small_m * sqrt( sigma_sqd_x_1 )
        1016             :        x_2_n         = -sgn_wpxp * small_m * sqrt( sigma_sqd_x_2 )
        1017             :        x_1           = xm
        1018             :        x_2           = xm
        1019             :
        1020             :     endwhere  ! Widths non-zero
        1021             :
        1022             :
        1023           0 :     return
        1024             :
        1025             :   end subroutine close_Luhar_pdf
        1026             :
        1027             :   !=============================================================================
        1028    71480448 :   subroutine ADG1_ADG2_responder_params( nz, ngrdcol, &                   ! In
        1029    71480448 :                                          xm, xp2, wp2, sqrt_wp2, &        ! In
        1030    71480448 :                                          wpxp, w_1_n, w_2_n, mixt_frac, & ! In
        1031    71480448 :                                          sigma_sqd_w, beta, x_tol, &      ! In
        1032    71480448 :                                          x_1, x_2, varnce_x_1, &          ! Out
        1033    71480448 :                                          varnce_x_2, alpha_x )            ! Out
        1034             :
        1035             :     ! Description:
        1036             :     ! Calculates the PDF component means and PDF component variances for thl,
        1037             :     ! rt, or sclr when either the ADG1 PDF or ADG2 PDF are used.
        1038             :     !
        1039             :     ! The normalized variance for thl, rt, and sclr for "plume" 1 is:
        1040             :     !
        1041             :     ! { 1 - [1/(1-sigma_sqd_w)]*[ <w'x'>^2 / (<w'^2> * <x'^2>) ] / mixt_frac }
        1042             :     ! * { (1/3)*beta + mixt_frac*( 1 - (2/3)*beta ) };
        1043             :     !
        1044             :     ! where "x" stands for thl, rt, or sclr; "mixt_frac" is the weight of
        1045             :     ! Gaussian "plume" 1, and 0 <= beta <= 3.
        1046             :     !
        1047             :     ! The factor { (1/3)*beta + mixt_frac*( 1 - (2/3)*beta ) } does not depend
        1048             :     ! on which varable "x" stands for.  The factor is multiplied by 2 and
        1049             :     ! defined as width_factor_1.
        1050             :     !
        1051             :     ! The factor
        1052             :     ! { 1 - [1/(1-sigma_sqd_w)]*[ (w'x')^2 / (w'^2 * x'^2) ] / mixt_frac }
        1053             :     ! depends on which variable "x" stands for.  It is multiplied by 1/2 and
        1054             :     ! defined as alpha_x, where "x" stands for thl, rt, or sclr.
        1055             :
        1056             :     ! Vince Larson added a dimensionless factor so that the
        1057             :     ! width of plumes in theta_l, rt can vary.
        1058             :     ! beta is a constant defined in CLUBB's tunable parameters.
        1059             :     !   Set 0<beta<3.
        1060             :     ! beta=1.5_core_rknd recovers Chris Golaz' simplified formula.
        1061             :     ! 3 Nov 2003
        1062             :
        1063             :     ! References:
        1064             :
        1065             :     use constants_clubb, only: &
        1066             :         two,            & ! Constant(s)
        1067             :         one,            &
        1068             :         two_thirds,     &
        1069             :         one_half,       &
        1070             :         zero_threshold, &
        1071             :         w_tol_sqd
        1072             :
        1073             :     use clubb_precision, only: &
        1074             :         core_rknd    ! Variable(s)
        1075             :
        1076             :     implicit none
        1077             :
        1078             :     integer, intent(in) :: &
        1079             :       ngrdcol,  & ! Number of grid columns
        1080             :       nz          ! Number of vertical level
        1081             :
        1082             :     ! Input Variables
        1083             :     real ( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
        1084             :       xm,          & ! Mean of x (overall)                    [units vary]
        1085             :       xp2,         & ! Variance of x (overall)                [(units vary)^2]
        1086             :       wp2,         & ! Variance of w (overall)                [m^2/s^2]
        1087             :       sqrt_wp2,    & ! Square root of the variance of w       [m/s]
        1088             :       wpxp,        & ! Covariance of w and x                  [m/s (units vary)]
        1089             :       w_1_n,       & ! Normalized mean of w (1st PDF comp.)   [-]
        1090             :       w_2_n,       & ! Normalized mean of w (2nd PDF comp.)   [-]
        1091             :       mixt_frac,   & ! Mixture fraction                       [-]
        1092             :       sigma_sqd_w    ! Width of individual w plumes           [-]
        1093             :
        1094             :     real ( kind = core_rknd ), dimension(ngrdcol), intent(in) :: &
        1095             :       beta           ! CLUBB tunable parameter beta           [-]
        1096             :
        1097             :     real ( kind = core_rknd ), intent(in) :: &
        1098             :       x_tol          ! Tolerance value for x                  [units vary]
        1099             :
        1100             :     ! Output Variables
        1101             :     real ( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
        1102             :       x_1,        & ! Mean of x (1st PDF component)             [units vary]
        1103             :       x_2,        & ! Mean of x (2nd PDF component)             [units vary]
        1104             :       varnce_x_1, & ! Variance of x (1st PDF component)         [(units vary)^2]
        1105             :       varnce_x_2, & ! Variance of x (2nd PDF component)         [(units vary)^2]
        1106             :       alpha_x       ! Factor relating to normalized variance for x           [-]
        1107             :
        1108             :     ! Local Variables
        1109             :     ! variables for a generalization of Chris Golaz' closure
        1110             :     ! varies width of plumes in theta_l, rt
        1111             :     real ( kind = core_rknd ) :: &
        1112             :       width_factor_1 ! Width factor relating to PDF component 1    [-]
        1113             :
        1114             :     integer :: &
        1115             :       k, i     ! Vertical loop index
        1116             :
        1117             :     !----- Begin Code -----
        1118             :     !$acc parallel loop gang vector collapse(2) default(present)
        1119  6147318528 :     do k = 1, nz
        1120 >10152*10^7 :       do i = 1, ngrdcol
        1121             :
        1122 >10145*10^7 :         if ( wp2(i,k) <= w_tol_sqd .or. xp2(i,k) <= x_tol**2 ) then
        1123             :
        1124             :             ! Vertical velocity doesn't vary.  Variable x is treated as a single
        1125             :             ! Gaussian in this situation.
        1126 71844499189 :             x_1(i,k)        = xm(i,k)
        1127 71844499189 :             x_2(i,k)        = xm(i,k)
        1128 71844499189 :             varnce_x_1(i,k) = xp2(i,k)
        1129 71844499189 :             varnce_x_2(i,k) = xp2(i,k)
        1130 71844499189 :             alpha_x(i,k)    = one_half
        1131             :
        1132             :         else
        1133             :
        1134 23532028811 :             x_1(i,k) = xm(i,k) - wpxp(i,k) / ( sqrt_wp2(i,k) * w_2_n(i,k) )
        1135 23532028811 :             x_2(i,k) = xm(i,k) - wpxp(i,k) / ( sqrt_wp2(i,k) * w_1_n(i,k) )
        1136             :
        1137             :             alpha_x(i,k) = one_half &
        1138             :                          * ( one - wpxp(i,k) * wpxp(i,k) &
        1139 23532028811 :                                    / ( ( one - sigma_sqd_w(i,k) ) * wp2(i,k) * xp2(i,k) ) )
        1140             :
        1141 23532028811 :             alpha_x(i,k) = max( min( alpha_x(i,k), one ), zero_threshold )
        1142             :
        1143             :             width_factor_1 = two_thirds * beta(i) &
        1144 23532028811 :                              + two * mixt_frac(i,k) * ( one - two_thirds * beta(i) )
        1145             :
        1146             :             ! Vince Larson multiplied original expressions by width_factor_1,2
        1147             :             !   to generalize scalar skewnesses.  05 Nov 03
        1148 23532028811 :             varnce_x_1(i,k) = width_factor_1 * xp2(i,k) * alpha_x(i,k) / mixt_frac(i,k)
        1149             :             varnce_x_2(i,k) = ( two - width_factor_1 ) * xp2(i,k) &
        1150 23532028811 :                               * alpha_x(i,k) / ( one - mixt_frac(i,k) )
        1151             :
        1152             :         end if
        1153             :
        1154             :       end do
        1155             :     end do
        1156             :     !$acc end parallel loop
        1157             :
        1158    71480448 :     return
        1159             :
        1160             :   end subroutine ADG1_ADG2_responder_params
        1161             :
        1162             :   !=============================================================================
        1163           0 :   subroutine backsolve_Luhar_params( Sk_max, Skx,          & ! In
        1164             :                                      big_m_max, mixt_frac, & ! In
        1165             :                                      big_m_x, small_m_x    ) ! Out
        1166             :
        1167             :     ! Description:
        1168             :     ! This subroutine calculates Luhar's big_m and small_m for the variate 'x'
        1169             :     ! consistent with the mixture fraction of the variate with the largest
        1170             :     ! skewness.
        1171             :     !
        1172             :     ! The relationship between skewness of x (Skx), mixture fraction (a), and
        1173             :     ! Luhar's small m (m) is given by:
        1174             :     !
        1175             :     ! Skx^2 = ( m^2 * ( m^2 + 3 )^2 / ( m^2 + 1 )^3 )
        1176             :     !         * ( 1 - 2*a )^2 / ( a * ( 1 - a ) ).
        1177             :     !
        1178             :     ! Moving the factor involving mixture fraction to the right-hand side:
        1179             :     !
        1180             :     ! ( ( a * ( 1 - a ) ) / ( 1 - 2*a )^2 ) * Skx^2
        1181             :     ! = m^2 * ( m^2 + 3 )^2 / ( m^2 + 1 )^3.
        1182             :     !
        1183             :     ! This can be rewritten as:
        1184             :     !
        1185             :     ! ( ( a * ( 1 - a ) ) / ( 1 - 2*a )^2 ) * Skx^2
        1186             :     ! = ( m^6 + 6*m^4 + 9*m^2 ) / ( m^6 + 3*m^4 + 3*m^2 + 1 ).
        1187             :     !
        1188             :     ! Setting alpha = ( ( a * ( 1 - a ) ) / ( 1 - 2*a )^2 ) * Skx^2, the
        1189             :     ! equation can be rewritten as:
        1190             :     !
        1191             :     ! ( m^6 + 3*m^4 + 3*m^2 + 1 ) * alpha = m^6 + 6*m^4 + 9*m^2.
        1192             :     !
        1193             :     ! This can be rearranged and rewritten as:
        1194             :     !
        1195             :     ! ( alpha - 1 ) * m^6 + ( 3 * alpha - 6 ) * m^4
        1196             :     ! + ( 3 * alpha - 9 ) * m^2 + alpha = 0.
        1197             :     !
        1198             :     ! This can be rewritten again as:
        1199             :     !
        1200             :     ! ( alpha - 1 ) * (m^2)^3 + ( 3 * alpha - 6 ) * (m^2)^2
        1201             :     ! + ( 3 * alpha - 9 ) * (m^2) + alpha = 0.
        1202             :     !
        1203             :     ! The goal is to solve for m^2, and then take the square root of m^2 to
        1204             :     ! solve for m.  This can be accomplished by using the cubic formula (with
        1205             :     ! the l_use_cubic_backsolve option), or else by a quadratic approximation.
        1206             :
        1207             :     ! References:
        1208             :     !-----------------------------------------------------------------------
        1209             :
        1210             :     use constants_clubb, only: &
        1211             :         nine,     & ! Constant(s)
        1212             :         six,      &
        1213             :         five,     &
        1214             :         four,     &
        1215             :         three,    &
        1216             :         two,      &
        1217             :         one,      &
        1218             :         one_half, &
        1219             :         zero,     &
        1220             :         eps
        1221             :
        1222             :     use clubb_precision, only: &
        1223             :         core_rknd     ! Precision
        1224             :
        1225             :     implicit none
        1226             :
        1227             :     ! Input Variables
        1228             :     real( kind = core_rknd ), intent(in) :: &
        1229             :       Sk_max,    & ! Maximum skewness                                      [-]
        1230             :       Skx,       & ! Skewness of the variate solving small_m and big_m for [-]
        1231             :       big_m_max, & ! Luhar's big_m of the variate with maximum skewness    [-]
        1232             :       mixt_frac    ! Mixture fraction                                      [-]
        1233             :
        1234             :     ! Output Variables
        1235             :     real( kind = core_rknd ), intent(out) :: &
        1236             :       big_m_x,   & ! Luhar's big_m for the variate being solved for    [-]
        1237             :       small_m_x    ! Luhar's small_m for the variate being solved for  [-]
        1238             :
        1239             :     ! Local Variables
        1240             :     real( kind = core_rknd ) :: &
        1241             :       alpha,     & ! 1 / big_m_x
        1242             :       a,         & ! For readability, quadratic equation
        1243             :       b,         &
        1244             :       c,         &
        1245             :       alpha_upr, &
        1246             :       alpha_low, &
        1247             :       discrim
        1248             :
        1249             :     real( kind = core_rknd ), dimension(1) :: &
        1250             :       a_coef_in, & ! Coefficient a (of x^3) in a*x^3 + b*x^2 + c^x + d = 0   [-]
        1251             :       b_coef_in, & ! Coefficient b (of x^2) in a*x^3 + b*x^2 + c^x + d = 0   [-]
        1252             :       c_coef_in, & ! Coefficient c (of x) in a*x^3 + b*x^2 + c^x + d = 0     [-]
        1253             :       d_coef_in    ! Coefficient d in a*x^3 + b*x^2 + c^x + d = 0            [-]
        1254             :
        1255             :     ! Flag to backsolve for m^2 using cubic formula
        1256             :     logical, parameter :: &
        1257             :       l_use_cubic_backsolve = .true.
        1258             :
        1259             :
        1260             :     if ( l_use_cubic_backsolve ) then
        1261             :
        1262           0 :        if ( abs( mixt_frac - one_half ) < 0.001_core_rknd ) then
        1263             :
        1264             :           ! When mixture fraction = 0.5 (based on the variable with the largest
        1265             :           ! magnitude of skewness), all variables must have a skewness of 0.
        1266             :           ! Set m to the minimum threshold of 0.05.
        1267           0 :           small_m_x = 0.05_core_rknd
        1268             :
        1269             :           ! Calculate the corresponding value of big_m_x.
        1270             :           big_m_x = ( one + small_m_x**2 )**3 &
        1271           0 :                     / ( ( three + small_m_x**2 )**2 * small_m_x**2 )
        1272             :
        1273           0 :        elseif ( abs(Skx) < eps ) then
        1274             :
        1275             :           ! Mixture fraction /= 0.5 because the variable with the largest
        1276             :           ! magnitude of skewness has a skewness /= 0.  However, variable x has
        1277             :           ! a skewness of 0.  In order to reproduce the correct skewness for
        1278             :           ! variable x, set m to 0 (regardless of minimum thresholds used in
        1279             :           ! other parts of the code).
        1280           0 :           small_m_x = zero
        1281             :
        1282             :           ! The value of big_m_x should be inf.  Set it to huge.  This is not
        1283             :           ! used in any calculation, anyway.
        1284           0 :           big_m_x = huge( big_m_x )
        1285             :
        1286             :        else  ! mixt_frac /= 0.5 and Skx /= 0
        1287             :
        1288             :           ! Backsolve for m, given mixt_frac and Skx.
        1289             :
        1290             :           ! alpha = 1/M is given by:
        1291             :           ! [ mixt_frac * ( 1 - mixt_frac ) / ( 1 - 2 * mixt_frac )^2 ] * Skx^2.
        1292             :           alpha = ( mixt_frac * ( one - mixt_frac ) &
        1293           0 :                     / ( one - two * mixt_frac )**2 ) * Skx**2
        1294             :
        1295             :           ! Calculate big_m_x.
        1296           0 :           big_m_x = one / alpha
        1297             :
        1298             :           ! Solve the cubic equation for m^2:
        1299             :           ! ( alpha - 1 ) * (m^2)^3 + ( 3 * alpha - 6 ) * (m^2)^2
        1300             :           ! + ( 3 * alpha - 9 ) * (m^2) + alpha = 0.
        1301             :           ! The largest root is preferred.
        1302           0 :           a_coef_in(1) = alpha - one
        1303           0 :           b_coef_in(1) = three * alpha - six
        1304           0 :           c_coef_in(1) = three * alpha - nine
        1305           0 :           d_coef_in(1) = alpha
        1306             :           small_m_x &
        1307             :           = sqrt( max( max_cubic_root( a_coef_in, b_coef_in, &
        1308             :                                        c_coef_in, d_coef_in ), &
        1309           0 :                        0.05_core_rknd**2 ) )
        1310             :
        1311             :        endif
        1312             :
        1313             :     else ! original formualation
        1314             :
        1315             :        alpha = ( Skx**2 / (max(Sk_max**2, eps) * big_m_max) )  ! 1 / big_m_x
        1316             :
        1317             :        ! This limit keeps the discriminant >= 0
        1318             :        alpha_upr = two * sqrt( 13.0_core_rknd ) - five
        1319             :
        1320             :        alpha_low = eps
        1321             :
        1322             :        ! For this approximation, alpha must be less than 2*sqrt(13) - 5 to get a
        1323             :        ! real ans.
        1324             :        alpha = min( alpha, alpha_upr )
        1325             :
        1326             :        ! For testing, eliminate possibility of divide by zero
        1327             :        alpha = max( alpha, alpha_low )
        1328             :
        1329             :        ! Use a piece-wise approximation
        1330             :        if ( alpha < one ) then
        1331             :
        1332             :           a = max( three * alpha - six, eps) ! Prevent divide by zero
        1333             :           b = three * alpha - nine
        1334             :           c = alpha
        1335             :
        1336             :           discrim = b**2 - four * a * c
        1337             :           small_m_x = sqrt( ( - b - sqrt( discrim ) ) / ( two * a ) )
        1338             :
        1339             :        else ! alpha >= 1
        1340             :
        1341             :           ! For this approximation, alpha must be less than 2*sqrt(13) - 5
        1342             :           ! to get a real ans.
        1343             :           alpha = min( alpha, two )
        1344             :
        1345             :           a = max( six * alpha - nine, eps) ! Prevent divide by zero
        1346             :           b = -six
        1347             :           c = two * alpha - one
        1348             :
        1349             :           discrim = b**2 - four * a * c
        1350             :           small_m_x = sqrt( ( - b - sqrt( discrim ) ) / ( two * a ) )
        1351             :
        1352             :        endif ! alpha < 1
        1353             :
        1354             :        ! Clip consistently with subroutine calc_Luhar_params
        1355             :        small_m_x = max( 5e-2_core_rknd, small_m_x)
        1356             :
        1357             :        big_m_x = one / alpha
        1358             :
        1359             :     endif ! l_use_cubic_backsolve
        1360             :
        1361             :
        1362           0 :     return
        1363             :
        1364             :   end subroutine backsolve_Luhar_params
        1365             :
        1366             :   !=============================================================================
        1367           0 :   function max_cubic_root( a_coef, b_coef, c_coef, d_coef ) &
        1368             :   result( max_root )
        1369             :
        1370             :     ! Description:
        1371             :     ! Calculates the largest root that results from solving a cubic equation of
        1372             :     ! the form a*x^3 + b*x^2 + c*x + d = 0.
        1373             :     !
        1374             :     ! This is done to backsolve for m^2 for the 3-D Luhar closure, given the
        1375             :     ! values of mixt_frac and Skx.
        1376             :
        1377             :     ! References:
        1378             :     !-----------------------------------------------------------------------
        1379             :
        1380             :     use constants_clubb, only: &
        1381             :         eps                ! Constant(s)
        1382             :
        1383             :     use calc_roots, only: &
        1384             :         cubic_solve,     & ! Procedure(s)
        1385             :         quadratic_solve
        1386             :
        1387             :     use clubb_precision, only: &
        1388             :         core_rknd ! Variable(s)
        1389             :
        1390             :     implicit none
        1391             :
        1392             :     ! Input Variables
        1393             :     real( kind = core_rknd ), dimension(1), intent(in) :: &
        1394             :       a_coef, & ! Coefficient a (of x^3) in a*x^3 + b*x^2 + c^x + d = 0    [-]
        1395             :       b_coef, & ! Coefficient b (of x^2) in a*x^3 + b*x^2 + c^x + d = 0    [-]
        1396             :       c_coef, & ! Coefficient c (of x) in a*x^3 + b*x^2 + c^x + d = 0      [-]
        1397             :       d_coef    ! Coefficient d in a*x^3 + b*x^2 + c^x + d = 0             [-]
        1398             :
        1399             :     ! Return Variable
        1400             :     real( kind = core_rknd ) :: &
        1401             :       max_root    ! Maximum root that solves the cubic equation            [-]
        1402             :
        1403             :     ! Local Variables
        1404             :     complex( kind = core_rknd ), dimension(1,3) :: &
        1405             :       cubic_roots    ! Roots of x that satisfy a*x^3 + b*x^2 + c*x + d = 0 [-]
        1406             :
        1407             :     complex( kind = core_rknd ), dimension(1,2) :: &
        1408             :       quadratic_roots    ! Roots of x that satisfy b*x^2 + c*x + d = 0     [-]
        1409             :
        1410             :     real( kind = core_rknd ) :: &
        1411             :       a_coef_thresh, & ! Minimum threshold of |a| to use cubic solver      [-]
        1412             :       b_coef_thresh    ! Minimum threshold of |b| to use quadratic solver  [-]
        1413             :
        1414             :
        1415             :     ! Calculate a minimum threshold for |a| to call this a cubic equation.
        1416             :     a_coef_thresh = 0.001_core_rknd &
        1417           0 :                     * max( abs(b_coef(1)), abs(c_coef(1)), abs(d_coef(1)) )
        1418             :
        1419             :     ! Calculate a minimum threshold for |b| to call this a quadratic equation.
        1420             :     ! This only matters when |a| <= a_coef_thresh.
        1421           0 :     b_coef_thresh = 0.001_core_rknd * max( abs(c_coef(1)), abs(d_coef(1)) )
        1422             :
        1423           0 :     if ( abs( a_coef(1) ) > a_coef_thresh ) then
        1424             :
        1425             :        ! The equation is a cubic equation.
        1426           0 :        cubic_roots = cubic_solve( 1, a_coef, b_coef, c_coef, d_coef )
        1427             :
        1428           0 :        if ( abs(aimag( cubic_roots(1,2) )) < eps .and.  &
        1429             :             abs(aimag( cubic_roots(1,3) )) < eps ) then
        1430             :
        1431             :           ! Find the maximum root of the three roots.
        1432             :           max_root = max( real( cubic_roots(1,1), kind = core_rknd ), &
        1433             :                           real( cubic_roots(1,2), kind = core_rknd ), &
        1434           0 :                           real( cubic_roots(1,3), kind = core_rknd ) )
        1435             :
        1436             :        else  ! cubic_roots(2) and cubic_roots(3) are complex.
        1437             :
        1438           0 :           max_root = real( cubic_roots(1,1), kind = core_rknd )
        1439             :
        1440             :        endif
        1441             :
        1442           0 :     elseif ( abs( b_coef(1) ) > b_coef_thresh ) then
        1443             :
        1444             :        ! The equation is a quadratic equation, since a = 0, but b /= 0.
        1445             :        ! This should very rarely occur for 3-D Luhar.  When it does, the result
        1446             :        ! will always be two real-valued roots.
        1447           0 :        quadratic_roots = quadratic_solve( 1, b_coef, c_coef, d_coef )
        1448             :
        1449             :        ! Find the maximum root of the two roots.
        1450             :        max_root = max( real( quadratic_roots(1,1), kind = core_rknd ), &
        1451           0 :                        real( quadratic_roots(1,2), kind = core_rknd ) )
        1452             :
        1453             :     else ! |a| = 0 and |b| = 0
        1454             :
        1455             :        ! The equation is a linear equation.
        1456             :        ! This won't happen for 3-D Luhar.
        1457           0 :        max_root = - d_coef(1) / c_coef(1)
        1458             :
        1459             :     endif ! |a| > 0
        1460             :
        1461             :
        1462             :     return
        1463             :
        1464             :   end function max_cubic_root
        1465             :
        1466             :   !=============================================================================
        1467             :
        1468             : end module adg1_adg2_3d_luhar_pdf

Generated by: LCOV version 1.14