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 17:57:11 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  5897779378 :              if ( abs( Skw(i,k) ) <= 1.0e-5_core_rknd ) then
     681   442655339 :                 mixt_frac(i,k) = one_half
     682             :              else
     683             :                 mixt_frac(i,k) &
     684             :                  = one_half &
     685  5455124039 :                      * ( 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  5897779378 :                                    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  5897779378 :               = 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  5897779378 :               = -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  5897779378 :              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  5897779378 :              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  5897779378 :              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  5897779378 :              varnce_w_2(i,k) = sigma_sqd_w(i,k) * wp2(i,k)
     718             : 
     719             :           else
     720             : 
     721             :              ! Vertical velocity doesn't vary.
     722 17946352622 :              mixt_frac(i,k)  = one_half
     723 17946352622 :              w_1_n(i,k)      = sqrt( one - sigma_sqd_w(i,k) )
     724 17946352622 :              w_2_n(i,k)      = -sqrt( one - sigma_sqd_w(i,k) )
     725 17946352622 :              w_1(i,k)        = wm(i,k)
     726 17946352622 :              w_2(i,k)        = wm(i,k)
     727 17946352622 :              varnce_w_1(i,k) = zero
     728 17946352622 :              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 72247479667 :             x_1(i,k)        = xm(i,k)
    1127 72247479667 :             x_2(i,k)        = xm(i,k)
    1128 72247479667 :             varnce_x_1(i,k) = xp2(i,k)
    1129 72247479667 :             varnce_x_2(i,k) = xp2(i,k)
    1130 72247479667 :             alpha_x(i,k)    = one_half
    1131             : 
    1132             :         else
    1133             : 
    1134 23129048333 :             x_1(i,k) = xm(i,k) - wpxp(i,k) / ( sqrt_wp2(i,k) * w_2_n(i,k) )
    1135 23129048333 :             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 23129048333 :                                    / ( ( one - sigma_sqd_w(i,k) ) * wp2(i,k) * xp2(i,k) ) )
    1140             : 
    1141 23129048333 :             alpha_x(i,k) = max( min( alpha_x(i,k), one ), zero_threshold )
    1142             : 
    1143             :             width_factor_1 = two_thirds * beta(i) &
    1144 23129048333 :                              + 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 23129048333 :             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 23129048333 :                               * 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