LCOV - code coverage report
Current view: top level - physics/clubb/src/CLUBB_core - new_tsdadg_pdf.F90 (source / functions) Hit Total Coverage
Test: coverage.info Lines: 0 121 0.0 %
Date: 2024-12-17 17:57:11 Functions: 0 4 0.0 %

          Line data    Source code
       1             : ! $Id$
       2             : !===============================================================================
       3             : module new_tsdadg_pdf
       4             : 
       5             :   ! Description:
       6             :   ! The new trivariate, skewness-dependent, analytic double Gaussian (TSDADG)
       7             :   ! PDF.
       8             : 
       9             :   implicit none
      10             : 
      11             :   public :: tsdadg_pdf_driver,      & ! Procedure(s)
      12             :             calc_setter_parameters, &
      13             :             calc_L_x_Skx_fnc
      14             : 
      15             :   private :: calc_respnder_parameters    ! Procedure(s)
      16             : 
      17             :   private  ! default scope
      18             : 
      19             :   contains
      20             : 
      21             :   !=============================================================================
      22           0 :   subroutine tsdadg_pdf_driver( nz, wm, rtm, thlm, wp2, rtp2, thlp2,    & ! In
      23           0 :                                 Skw, Skrt, Skthl, wprtp, wpthlp,    & ! In
      24           0 :                                 mu_w_1, mu_w_2,                     & ! Out
      25           0 :                                 mu_rt_1, mu_rt_2,                   & ! Out
      26           0 :                                 mu_thl_1, mu_thl_2,                 & ! Out
      27           0 :                                 sigma_w_1_sqd, sigma_w_2_sqd,       & ! Out
      28           0 :                                 sigma_rt_1_sqd, sigma_rt_2_sqd,     & ! Out
      29           0 :                                 sigma_thl_1_sqd, sigma_thl_2_sqd,   & ! Out
      30           0 :                                 mixt_frac )                           ! Out
      31             : 
      32             : 
      33             :     ! Description:
      34             :     ! Selects which variable is used to set the mixture fraction for the PDF
      35             :     ! ("the setter") and which variables are handled after that mixture fraction
      36             :     ! has been set ("the responders").  Traditionally, w has been used to set
      37             :     ! the PDF.  However, here, the variable with the greatest magnitude of
      38             :     ! skewness is used to set the PDF.
      39             : 
      40             :     ! References:
      41             :     !-----------------------------------------------------------------------
      42             : 
      43             :     use constants_clubb, only: &
      44             :         one,     & ! Variable(s)
      45             :         zero,    &
      46             :         fstderr
      47             : 
      48             :     use clubb_precision, only: &
      49             :         core_rknd    ! Variable(s)
      50             : 
      51             :     implicit none
      52             : 
      53             :     integer, intent(in) :: &
      54             :       nz
      55             : 
      56             :     ! Input Variables
      57             :     real( kind = core_rknd ), dimension(nz), intent(in) :: &
      58             :       wm,     & ! Mean of w (overall)                [m/s]
      59             :       rtm,    & ! Mean of rt (overall)               [kg/kg]
      60             :       thlm,   & ! Mean of thl (overall)              [K]
      61             :       wp2,    & ! Variance of w (overall)            [m^2/s^2]
      62             :       rtp2,   & ! Variance of rt (overall)           [kg^2/kg^2]
      63             :       thlp2,  & ! Variance of thl (overall)          [K^2]
      64             :       Skw,    & ! Skewness of w (overall)            [-]
      65             :       Skrt,   & ! Skewness of rt (overall)           [-]
      66             :       Skthl,  & ! Skewness of thl (overall)          [-]
      67             :       wprtp,  & ! Covariance of w and rt (overall)   [(m/s)kg/kg]
      68             :       wpthlp    ! Covariance of w and thl (overall)  [(m/s)K]
      69             : 
      70             :     ! Output Variables
      71             :     real( kind = core_rknd ), dimension(nz), intent(out) :: &
      72             :       mu_w_1,          & ! Mean of w (1st PDF component)        [m/s]
      73             :       mu_w_2,          & ! Mean of w (2nd PDF component)        [m/s]
      74             :       mu_rt_1,         & ! Mean of rt (1st PDF component)       [kg/kg]
      75             :       mu_rt_2,         & ! Mean of rt (2nd PDF component)       [kg/kg]
      76             :       mu_thl_1,        & ! Mean of thl (1st PDF component)      [K]
      77             :       mu_thl_2,        & ! Mean of thl (2nd PDF component)      [K]
      78             :       sigma_w_1_sqd,   & ! Variance of w (1st PDF component)    [m^2/s^2]
      79             :       sigma_w_2_sqd,   & ! Variance of w (2nd PDF component)    [m^2/s^2]
      80             :       sigma_rt_1_sqd,  & ! Variance of rt (1st PDF component)   [kg^2/kg^2]
      81             :       sigma_rt_2_sqd,  & ! Variance of rt (2nd PDF component)   [kg^2/kg^2]
      82             :       sigma_thl_1_sqd, & ! Variance of thl (1st PDF component)  [K^2]
      83             :       sigma_thl_2_sqd, & ! Variance of thl (2nd PDF component)  [K^2]
      84             :       mixt_frac          ! Mixture fraction                     [-]
      85             : 
      86             :     ! Local Variables
      87             :     real( kind = core_rknd ), dimension(nz) :: &
      88           0 :       big_L_w_1,   & ! Parameter for the 1st PDF comp. mean of w            [-]
      89           0 :       big_L_w_2,   & ! Parameter for the 2nd PDF comp. mean of w (setter)   [-]
      90           0 :       big_L_rt_1,  & ! Parameter for the 1st PDF comp. mean of rt           [-]
      91           0 :       big_L_rt_2,  & ! Parameter for the 2nd PDF comp. mean of rt (setter)  [-]
      92           0 :       big_L_thl_1, & ! Parameter for the 1st PDF comp. mean of thl          [-]
      93           0 :       big_L_thl_2    ! Parameter for the 2nd PDF comp. mean of thl (setter) [-]
      94             : 
      95             :     real( kind = core_rknd ) :: &
      96             :       small_l_w_1,   & ! Param. for the 1st PDF comp. mean of w            [-]
      97             :       small_l_w_2,   & ! Param. for the 2nd PDF comp. mean of w (setter)   [-]
      98             :       small_l_rt_1,  & ! Param. for the 1st PDF comp. mean of rt           [-]
      99             :       small_l_rt_2,  & ! Param. for the 2nd PDF comp. mean of rt (setter)  [-]
     100             :       small_l_thl_1, & ! Param. for the 1st PDF comp. mean of thl          [-]
     101             :       small_l_thl_2    ! Param. for the 2nd PDF comp. mean of thl (setter) [-]
     102             : 
     103             :     real( kind = core_rknd ), dimension(nz) :: &
     104           0 :       sgn_wprtp,  & ! Sign of the covariance of w and rt (overall)         [-]
     105           0 :       sgn_wpthlp, & ! Sign of the covariance of w and thl (overall)        [-]
     106           0 :       sgn_wp2       ! Sign of the variance of w (overall); always positive [-]
     107             : 
     108             :     real( kind = core_rknd ), dimension(nz) :: &
     109           0 :       coef_sigma_w_1_sqd,   & ! sigma_w_1^2 = coef_sigma_w_1_sqd * <w'^2>    [-]
     110           0 :       coef_sigma_w_2_sqd,   & ! sigma_w_2^2 = coef_sigma_w_2_sqd * <w'^2>    [-]
     111           0 :       coef_sigma_rt_1_sqd,  & ! sigma_rt_1^2 = coef_sigma_rt_1_sqd * <rt'^2> [-]
     112           0 :       coef_sigma_rt_2_sqd,  & ! sigma_rt_2^2 = coef_sigma_rt_2_sqd * <rt'^2> [-]
     113           0 :       coef_sigma_thl_1_sqd, & ! sigma_thl_1^2=coef_sigma_thl_1_sqd*<thl'^2>  [-]
     114           0 :       coef_sigma_thl_2_sqd    ! sigma_thl_2^2=coef_sigma_thl_2_sqd*<thl'^2>  [-]
     115             : 
     116             :     integer :: k    ! Vertical level loop index
     117             : 
     118             : 
     119             :     ! Calculate sgn( <w'rt'> ).
     120           0 :     where ( wprtp >= zero )
     121             :        sgn_wprtp = one
     122             :     elsewhere ! wprtp < 0
     123             :        sgn_wprtp = -one
     124             :     endwhere ! wprtp >= 0
     125             : 
     126             :     ! Calculate sgn( <w'thl'> ).
     127           0 :     where ( wpthlp >= zero )
     128             :        sgn_wpthlp = one
     129             :     elsewhere ! wpthlp < 0
     130             :        sgn_wpthlp = -one
     131             :     endwhere ! wpthlp >= 0
     132             : 
     133             :     ! The sign of the variance of w is always positive.
     134           0 :     sgn_wp2 = one
     135             : 
     136           0 :     small_l_w_1 = 0.75_core_rknd
     137           0 :     small_l_w_2 = 0.5_core_rknd
     138           0 :     small_l_rt_1 = 0.75_core_rknd
     139           0 :     small_l_rt_2 = 0.5_core_rknd
     140           0 :     small_l_thl_1 = 0.75_core_rknd
     141           0 :     small_l_thl_2 = 0.5_core_rknd
     142             : 
     143           0 :     do k = 1, nz, 1
     144             : 
     145           0 :        call calc_L_x_Skx_fnc( Skw(k), sgn_wp2(k),        & ! In
     146             :                               small_l_w_1, small_l_w_2,  & ! In
     147           0 :                               big_L_w_1(k), big_L_w_2(k) ) ! Out
     148             : 
     149             :        call calc_L_x_Skx_fnc( Skrt(k), sgn_wprtp(k),       & ! In
     150             :                               small_l_rt_1, small_l_rt_2,  & ! In
     151             :                               big_L_rt_1(k), big_L_rt_2(k) ) ! Out
     152             : 
     153             :        call calc_L_x_Skx_fnc( Skthl(k), sgn_wpthlp(k),       & ! In
     154             :                               small_l_thl_1, small_l_thl_2,  & ! In
     155             :                               big_L_thl_1(k), big_L_thl_2(k) ) ! Out
     156             : 
     157             : 
     158             :        ! The variable with the greatest magnitude of skewness will be the setter
     159             :        ! variable and the other variables will be responder variables.
     160             :        if ( abs( Skw(k) ) >= abs( Skrt(k) ) &
     161           0 :             .and. abs( Skw(k) ) >= abs( Skthl(k) ) ) then
     162             : 
     163             :           ! The variable w has the greatest magnitude of skewness.
     164             : 
     165             :           call calc_setter_parameters( wm(k), wp2(k),                  & ! In
     166             :                                        Skw(k), sgn_wp2(k),             & ! In
     167             :                                        big_L_w_1(k), big_L_w_2(k),     & ! In
     168             :                                        mu_w_1(k), mu_w_2(k),           & ! Out
     169             :                                        sigma_w_1_sqd(k),               & ! Out
     170             :                                        sigma_w_2_sqd(k), mixt_frac(k), & ! Out
     171             :                                        coef_sigma_w_1_sqd(k),          & ! Out
     172             :                                        coef_sigma_w_2_sqd(k)           ) ! Out
     173             : 
     174             :           call calc_respnder_parameters( rtm(k), rtp2(k),             & ! In
     175             :                                          Skrt(k), sgn_wprtp(k),       & ! In
     176             :                                          mixt_frac(k), big_L_rt_1(k), & ! In
     177             :                                          mu_rt_1(k), mu_rt_2(k),      & ! Out
     178             :                                          sigma_rt_1_sqd(k),           & ! Out
     179             :                                          sigma_rt_2_sqd(k),           & ! Out
     180             :                                          coef_sigma_rt_1_sqd(k),      & ! Out
     181             :                                          coef_sigma_rt_2_sqd(k)       ) ! Out
     182             : 
     183             :           call calc_respnder_parameters( thlm(k), thlp2(k),            & ! In
     184             :                                          Skthl(k), sgn_wpthlp(k),      & ! In
     185             :                                          mixt_frac(k), big_L_thl_1(k), & ! In
     186             :                                          mu_thl_1(k), mu_thl_2(k),     & ! Out
     187             :                                          sigma_thl_1_sqd(k),           & ! Out
     188             :                                          sigma_thl_2_sqd(k),           & ! Out
     189             :                                          coef_sigma_thl_1_sqd(k),      & ! Out
     190             :                                          coef_sigma_thl_2_sqd(k)       ) ! Out
     191             : 
     192             : 
     193             :        elseif ( abs( Skrt(k) ) > abs( Skw(k) ) &
     194           0 :                 .and. abs( Skrt(k) ) >= abs( Skthl(k) ) ) then
     195             : 
     196             :           ! The variable rt has the greatest magnitude of skewness.
     197             : 
     198             :           call calc_setter_parameters( rtm(k), rtp2(k),                 & ! In
     199             :                                        Skrt(k), sgn_wprtp(k),           & ! In
     200             :                                        big_L_rt_1(k), big_L_rt_2(k),    & ! In
     201             :                                        mu_rt_1(k), mu_rt_2(k),          & ! Out
     202             :                                        sigma_rt_1_sqd(k),               & ! Out
     203             :                                        sigma_rt_2_sqd(k), mixt_frac(k), & ! Out
     204             :                                        coef_sigma_rt_1_sqd(k),          & ! Out
     205             :                                        coef_sigma_rt_2_sqd(k)           ) ! Out
     206             : 
     207             :           call calc_respnder_parameters( wm(k), wp2(k),              & ! In
     208             :                                          Skw(k), sgn_wp2(k),         & ! In
     209             :                                          mixt_frac(k), big_L_w_1(k), & ! In
     210             :                                          mu_w_1(k), mu_w_2(k),       & ! Out
     211             :                                          sigma_w_1_sqd(k),           & ! Out
     212             :                                          sigma_w_2_sqd(k),           & ! Out
     213             :                                          coef_sigma_w_1_sqd(k),      & ! Out
     214             :                                          coef_sigma_w_2_sqd(k)       ) ! Out
     215             : 
     216             :           call calc_respnder_parameters( thlm(k), thlp2(k),            & ! In
     217             :                                          Skthl(k), sgn_wpthlp(k),      & ! In
     218             :                                          mixt_frac(k), big_L_thl_1(k), & ! In
     219             :                                          mu_thl_1(k), mu_thl_2(k),     & ! Out
     220             :                                          sigma_thl_1_sqd(k),           & ! Out
     221             :                                          sigma_thl_2_sqd(k),           & ! Out
     222             :                                          coef_sigma_thl_1_sqd(k),      & ! Out
     223             :                                          coef_sigma_thl_2_sqd(k)       ) ! Out
     224             : 
     225             : 
     226             :        else ! abs( Skthl ) > abs( Skw ) .and. abs( Skthl ) > abs( Skrt )
     227             : 
     228             :           ! The variable thl has the greatest magnitude of skewness.
     229             : 
     230             :           call calc_setter_parameters( thlm(k), thlp2(k),                & ! In
     231             :                                        Skthl(k), sgn_wpthlp(k),          & ! In
     232             :                                        big_L_thl_1(k), big_L_thl_2(k),   & ! In
     233             :                                        mu_thl_1(k), mu_thl_2(k),         & ! Out
     234             :                                        sigma_thl_1_sqd(k),               & ! Out
     235             :                                        sigma_thl_2_sqd(k), mixt_frac(k), & ! Out
     236             :                                        coef_sigma_thl_1_sqd(k),          & ! Out
     237             :                                        coef_sigma_thl_2_sqd(k)           ) ! Out
     238             : 
     239             :           call calc_respnder_parameters( wm(k), wp2(k),              & ! In
     240             :                                          Skw(k), sgn_wp2(k),         & ! In
     241             :                                          mixt_frac(k), big_L_w_1(k), & ! In
     242             :                                          mu_w_1(k), mu_w_2(k),       & ! Out
     243             :                                          sigma_w_1_sqd(k),           & ! Out
     244             :                                          sigma_w_2_sqd(k),           & ! Out
     245             :                                          coef_sigma_w_1_sqd(k),      & ! Out
     246             :                                          coef_sigma_w_2_sqd(k)       ) ! Out
     247             : 
     248             :           call calc_respnder_parameters( rtm(k), rtp2(k),             & ! In
     249             :                                          Skrt(k), sgn_wprtp(k),       & ! In
     250             :                                          mixt_frac(k), big_L_rt_1(k), & ! In
     251             :                                          mu_rt_1(k), mu_rt_2(k),      & ! Out
     252             :                                          sigma_rt_1_sqd(k),           & ! Out
     253             :                                          sigma_rt_2_sqd(k),           & ! Out
     254             :                                          coef_sigma_rt_1_sqd(k),      & ! Out
     255             :                                          coef_sigma_rt_2_sqd(k)       ) ! Out
     256             : 
     257             : 
     258             :        endif ! Find variable with the greatest magnitude of skewness.
     259             : 
     260             : 
     261           0 :        if ( sigma_w_1_sqd(k) < zero ) then
     262             :           write(fstderr,*) "WARNING:  New TSDADG PDF.  The variance of w in " &
     263             :                            // "the 1st PDF component is negative and is " &
     264           0 :                            // "being clipped to 0."
     265           0 :           write(fstderr,*) "sigma_w_1^2 (before clipping) = ", sigma_w_1_sqd(k)
     266           0 :           sigma_w_1_sqd(k) = zero
     267           0 :           coef_sigma_w_1_sqd(k) = zero
     268             :        endif ! sigma_w_1_sqd < 0
     269             : 
     270           0 :        if ( sigma_w_2_sqd(k) < zero ) then
     271             :           write(fstderr,*) "WARNING:  New TSDADG PDF.  The variance of w in " &
     272             :                            // "the 2nd PDF component is negative and is " &
     273           0 :                            // "being clipped to 0."
     274           0 :           write(fstderr,*) "sigma_w_2^2 (before clipping) = ", sigma_w_2_sqd(k)
     275           0 :           sigma_w_2_sqd(k) = zero
     276           0 :           coef_sigma_w_2_sqd(k) = zero
     277             :        endif ! sigma_w_2_sqd < 0
     278             : 
     279           0 :        if ( sigma_rt_1_sqd(k) < zero ) then
     280             :           write(fstderr,*) "WARNING:  New TSDADG PDF.  The variance of rt in " &
     281             :                            // "the 1st PDF component is negative and is " &
     282           0 :                            // "being clipped to 0."
     283           0 :           write(fstderr,*) "sigma_rt_1^2 (before clipping) = ", &
     284           0 :                            sigma_rt_1_sqd(k)
     285           0 :           sigma_rt_1_sqd(k) = zero
     286           0 :           coef_sigma_rt_1_sqd(k) = zero
     287             :        endif ! sigma_rt_1_sqd < 0
     288             : 
     289           0 :        if ( sigma_rt_2_sqd(k) < zero ) then
     290             :           write(fstderr,*) "WARNING:  New TSDADG PDF.  The variance of rt in " &
     291             :                            // "the 2nd PDF component is negative and is " &
     292           0 :                            // "being clipped to 0."
     293           0 :           write(fstderr,*) "sigma_rt_2^2 (before clipping) = ", &
     294           0 :                            sigma_rt_2_sqd(k)
     295           0 :           sigma_rt_2_sqd(k) = zero
     296           0 :           coef_sigma_rt_2_sqd(k) = zero
     297             :        endif ! sigma_rt_2_sqd < 0
     298             : 
     299           0 :        if ( sigma_thl_1_sqd(k) < zero ) then
     300             :           write(fstderr,*) "WARNING:  New TSDADG PDF.  The variance of thl " &
     301             :                            // "in the 1st PDF component is negative and is " &
     302           0 :                            // "being clipped to 0."
     303           0 :           write(fstderr,*) "sigma_thl_1^2 (before clipping) = ", &
     304           0 :                            sigma_thl_1_sqd(k)
     305           0 :           sigma_thl_1_sqd(k) = zero
     306           0 :           coef_sigma_thl_1_sqd(k) = zero
     307             :        endif ! sigma_thl_1_sqd < 0
     308             : 
     309           0 :        if ( sigma_thl_2_sqd(k) < zero ) then
     310             :           write(fstderr,*) "WARNING:  New TSDADG PDF.  The variance of thl " &
     311             :                            // "in the 2nd PDF component is negative and is " &
     312           0 :                            // "being clipped to 0."
     313           0 :           write(fstderr,*) "sigma_thl_2^2 (before clipping) = ", &
     314           0 :                            sigma_thl_2_sqd(k)
     315           0 :           sigma_thl_2_sqd(k) = zero
     316           0 :           coef_sigma_thl_2_sqd(k) = zero
     317             :        endif ! sigma_thl_2_sqd < 0
     318             : 
     319             :     enddo ! k = 1, nz, 1
     320             : 
     321             : 
     322           0 :     return
     323             : 
     324             :   end subroutine tsdadg_pdf_driver
     325             : 
     326             :   !=============================================================================
     327             :   !
     328             :   ! DESCRIPTION OF THE METHOD FOR THE VARIABLE THAT SETS THE MIXTURE FRACTION
     329             :   ! =========================================================================
     330             :   !
     331             :   ! There are five PDF parameters that need to be calculated for the setting
     332             :   ! variable, which are mu_x_1 (the mean of x in the 1st PDF component), mu_x_2
     333             :   ! (the mean of x in the 2nd PDF component), sigma_x_1 (the standard deviation
     334             :   ! of x in the 1st PDF component), sigma_x_2 (the standard deviation of x in
     335             :   ! the 2nd PDF component), and mixt_frac (the mixture fraction).  In order to
     336             :   ! solve for these five parameters, five equations are needed.  These five
     337             :   ! equations are the equations for <x>, <x'^2>, and <x'^3> as found by
     338             :   ! integrating over the PDF.  Additionally, two more equations, which involve
     339             :   ! tunable parameters L_x_1 and L_x_2, and which are used to help control the
     340             :   ! spread of the PDF component means in the 1st PDF component and the 2nd PDF
     341             :   ! component, respectively, are used in this equation set.  The five equations
     342             :   ! are:
     343             :   !
     344             :   ! <x> = mixt_frac * mu_x_1 + ( 1 - mixt_frac ) * mu_x_2;
     345             :   !
     346             :   ! <x'^2> = mixt_frac * ( ( mu_x_1 - <x> )^2 + sigma_x_1^2 )
     347             :   !          + ( 1 - mixt_frac ) * ( ( mu_x_2 - <x> )^2 + sigma_x_2^2 );
     348             :   !
     349             :   ! <x'^3> = mixt_frac * ( mu_x_1 - <x> )
     350             :   !                    * ( ( mu_x_1 - <x> )^2 + 3 * sigma_x_1^2 )
     351             :   !          + ( 1 - mixt_frac ) * ( mu_x_2 - <x> )
     352             :   !                              * ( ( mu_x_2 - <x> )^2 + 3 * sigma_x_2^2 );
     353             :   !
     354             :   ! mu_x_1 - <x> = L_x_1
     355             :   !                * sqrt( ( 1 + Skx * sgn( <w'x'> ) / sqrt( 4 + Skx^2 ) )
     356             :   !                        / ( 1 - Skx * sgn( <w'x'> ) / sqrt( 4 + Skx^2 ) ) )
     357             :   !                * sqrt( <x'^2> ) * sgn( <w'x'> ); and
     358             :   !
     359             :   ! mu_x_2 - <x> = -L_x_2
     360             :   !                 * sqrt( ( 1 - Skx * sgn( <w'x'> ) / sqrt( 4 + Skx^2 ) )
     361             :   !                         / ( 1 + Skx * sgn( <w'x'> ) / sqrt( 4 + Skx^2 ) ) )
     362             :   !                 * sqrt( <x'^2> ) * sgn( <w'x'> );
     363             :   !
     364             :   ! where 0 <= L_x_1 <= 1, 0 <= L_x_2 <= 1, Skx is the skewness of x, such that
     365             :   ! Skx = <x'^3> / <x'^2>^(3/2), and sgn( <w'x'> ) is the sign of <w'x'>, such
     366             :   ! that:
     367             :   !
     368             :   ! sgn( <w'x'> ) = |  1, when <w'x'> >= 0;
     369             :   !                 | -1, when <w'x'> < 0.
     370             :   !
     371             :   ! The resulting equations for the five PDF parameters are:
     372             :   !
     373             :   ! mu_x_1 = <x> + L_x_1
     374             :   !                * sqrt( ( 1 + Skx * sgn( <w'x'> ) / sqrt( 4 + Skx^2 ) )
     375             :   !                        / ( 1 - Skx * sgn( <w'x'> ) / sqrt( 4 + Skx^2 ) ) )
     376             :   !                * sqrt( <x'^2> ) * sgn( <w'x'> );
     377             :   !
     378             :   ! mu_x_2 = <x> - L_x_2
     379             :   !                * sqrt( ( 1 - Skx * sgn( <w'x'> ) / sqrt( 4 + Skx^2 ) )
     380             :   !                        / ( 1 + Skx * sgn( <w'x'> ) / sqrt( 4 + Skx^2 ) ) )
     381             :   !                * sqrt( <x'^2> ) * sgn( <w'x'> );
     382             :   !
     383             :   ! mixt_frac = 1 / ( 1 + abs( mu_x_1_nrmlized / mu_x_2_nrmlized ) );
     384             :   !
     385             :   ! sigma_x_1 = sqrt( ( 1 - mixt_frac * mu_x_1_nrmlized^2
     386             :   !                     - ( 1 - mixt_frac ) * mu_x_2_nrmlized^2
     387             :   !                     + ( 1 - mixt_frac )
     388             :   !                       * ( Skx / ( 3 * mixt_frac * mu_x_1_nrmlized )
     389             :   !                           - mu_x_1_nrmlized^2 / 3
     390             :   !                           + mu_x_2_nrmlized^2 / 3 ) )
     391             :   !                   * <x'^2> ); and
     392             :   !
     393             :   ! sigma_x_2 = sqrt( ( 1 - mixt_frac * mu_x_1_nrmlized^2
     394             :   !                     - ( 1 - mixt_frac ) * mu_x_2_nrmlized^2
     395             :   !                     - mixt_frac
     396             :   !                       * ( Skx / ( 3 * mixt_frac * mu_x_1_nrmlized )
     397             :   !                           - mu_x_1_nrmlized^2 / 3
     398             :   !                           + mu_x_2_nrmlized^2 / 3 ) )
     399             :   !                   * <x'^2> ); where
     400             :   !
     401             :   ! mu_x_1_nrmlized = ( mu_x_1 - <x> ) / sqrt( <x'^2> ); and
     402             :   !
     403             :   ! mu_x_2_nrmlized = ( mu_x_2 - <x> ) / sqrt( <x'^2> ).
     404             :   !
     405             :   !
     406             :   ! Notes:
     407             :   !
     408             :   ! This method does NOT work for all values of L_x_1 and L_x_2 (where
     409             :   ! 0 <= L_x_1 <= 1 and 0 <= L_x_2 <= 1).  Only a subregion of this parameter
     410             :   ! space produces valid results.
     411             :   !
     412             :   ! When both L_x_1 = 0 and L_x_2 = 0, mu_x_1 = mu_x_2 = <x> (which can only
     413             :   ! happen when Skx = 0).  In this scenario, the above equations for mixt_frac,
     414             :   ! sigma_x_1, and sigma_x_2 are all undefined.  In this special case, the
     415             :   ! distribution reduces to a single Gaussian, so the following values are set:
     416             :   ! mixt_frac = 1/2, sigma_x_1 = sqrt( <x'^2> ), and sigma_x_2 = sqrt( <x'^2> ).
     417             :   !
     418             :   !
     419             :   ! Tunable parameters:
     420             :   !
     421             :   ! The parameter L_x_1 controls the 1st PDF component mean while L_x_2 controls
     422             :   ! the 2nd PDF component mean.  The equations involving the tunable parameters
     423             :   ! L_x_1 and L_x_2 (the mu_x_1 and mu_x_2 equations) are based on the values of
     424             :   ! mu_x_1 and mu_x_2 when sigma_x_1 = sigma_x_2 = 0.  In this scenario, the
     425             :   ! equation for mixture fraction reduces to:
     426             :   !
     427             :   ! mixt_frac = (1/2) * ( 1 +/- Skx / sqrt( 4 + Skx^2 ) ).
     428             :   !
     429             :   ! The +/- is dependent on the sign of ( mu_x_1 - <x> ) vs. ( mu_x_2 - <x> ).
     430             :   ! This is dependent on sgn( <w'x'> ), and the mixture fraction equation is
     431             :   ! written as:
     432             :   !
     433             :   ! mixt_frac = (1/2) * ( 1 - Skx * sgn( <w'x'> ) / sqrt( 4 + Skx^2 ) ).
     434             :   !
     435             :   ! Meanwhile, the equation for 1 - mixt_frac is:
     436             :   !
     437             :   ! 1 - mixt_frac = (1/2) * ( 1 + Skx * sgn( <w'x'> ) / sqrt( 4 + Skx^2 ) ).
     438             :   !
     439             :   ! When sigma_x_1 = sigma_x_2 = 0, the equations for mu_x_1 and mu_x_2 are:
     440             :   !
     441             :   ! mu_x_1 = <x> + sqrt( ( 1 - mixt_frac ) / mixt_frac ) * sqrt( <x'^2> )
     442             :   !                * sgn( <w'x'> ); and
     443             :   !
     444             :   ! mu_x_2 = <x> - sqrt( mixt_frac / ( 1 - mixt_frac ) ) * sqrt( <x'^2> )
     445             :   !                * sgn( <w'x'> ).
     446             :   !
     447             :   ! Substituting the equations for mixt_frac and 1 - mixt_frac into the
     448             :   ! equations for mu_x_1 and mu_x_2 (when sigma_x_1 = sigma_x_2 = 0), the
     449             :   ! equations for mu_x_1 and mu_x_2 become:
     450             :   !
     451             :   ! mu_x_1 = <x> + sqrt( ( 1 + Skx * sgn( <w'x'> ) / sqrt( 4 + Skx^2 ) )
     452             :   !                      / ( 1 - Skx * sgn( <w'x'> ) / sqrt( 4 + Skx^2 ) ) )
     453             :   !                * sqrt( <x'^2> ) * sgn( <w'x'> ); and
     454             :   !
     455             :   ! mu_x_2 = <x> - sqrt( ( 1 - Skx * sgn( <w'x'> ) / sqrt( 4 + Skx^2 ) )
     456             :   !                      / ( 1 + Skx * sgn( <w'x'> ) / sqrt( 4 + Skx^2 ) ) )
     457             :   !                * sqrt( <x'^2> ) * sgn( <w'x'> ).
     458             :   !
     459             :   ! These equations represent the maximum deviation of mu_x_1 and mu_x_2 from
     460             :   ! the overall mean, <x>.  The range of parameters of L_x_i is 0 <= L_x_i <= 1.
     461             :   ! When L_x_1 = L_x_2 = 0, the value of mu_x_1 = mu_x_2 = <x> (and the
     462             :   ! distribution becomes a single Gaussian).  When L_x_i = 1, the value of
     463             :   ! mu_x_i - <x> is at its maximum magnitude.
     464             :   !
     465             :   ! The values of L_x_1 and L_x_2 are also calculated by skewness functions.
     466             :   ! Those functions are:
     467             :   !
     468             :   ! L_x_1 = l_x_1 * abs( Skx ) / sqrt( 4 + Skx^2 ); and
     469             :   ! L_x_2 = l_x_2 * abs( Skx ) / sqrt( 4 + Skx^2 );
     470             :   !
     471             :   ! where both l_x_1 and l_x_2 are tunable parameters.
     472             :   !
     473             :   ! As previously stated, this method does not work for all combinations of
     474             :   ! L_x_1 and L_x_2, but rather only for a subregion of parameter space.  This
     475             :   ! applies to l_x_1 and l_x_2, as well.  The conditions on l_x_1 and l_x_2 are:
     476             :   !
     477             :   ! 2/3 < l_x_1 < 1 and 0 < l_x_2 < 1; when Skx * sgn( <w'x'> ) >= 0; and
     478             :   ! 0 < l_x_1 < 1 and 2/3 < l_x_2 < 1; when Skx * sgn( <w'x'> ) < 0.
     479             :   !
     480             :   ! The condition that l_x_1 > 2/3 prevents a negative PDF component variance
     481             :   ! when Skx = 0.
     482             :   !
     483             :   !
     484             :   ! Equations for PDF component standard deviations:
     485             :   !
     486             :   ! The equations for the PDF component standard deviations can also be written
     487             :   ! as:
     488             :   !
     489             :   ! sigma_x_1 = sqrt( coef_sigma_x_1_sqd * <x'^2> ); and
     490             :   !
     491             :   ! sigma_x_2 = sqrt( coef_sigma_x_2_sqd * <x'^2> ); where
     492             :   !
     493             :   ! coef_sigma_x_1_sqd = 1 - mixt_frac * mu_x_1_nrmlized^2
     494             :   !                      - ( 1 - mixt_frac ) * mu_x_2_nrmlized^2
     495             :   !                      + ( 1 - mixt_frac )
     496             :   !                        * ( Skx / ( 3 * mixt_frac * mu_x_1_nrmlized )
     497             :   !                            - mu_x_1_nrmlized^2 / 3
     498             :   !                            + mu_x_2_nrmlized^2 / 3 ); and
     499             :   !
     500             :   ! coef_sigma_x_2_sqd = 1 - mixt_frac * mu_x_1_nrmlized^2
     501             :   !                      - ( 1 - mixt_frac ) * mu_x_2_nrmlized^2
     502             :   !                      - mixt_frac
     503             :   !                        * ( Skx / ( 3 * mixt_frac * mu_x_1_nrmlized )
     504             :   !                            - mu_x_1_nrmlized^2 / 3
     505             :   !                            + mu_x_2_nrmlized^2 / 3 ).
     506             :   !
     507             :   ! The above equations can be substituted into an equation for a variable that
     508             :   ! has been derived by integrating over the PDF.  Many variables like this are
     509             :   ! used in parts of the predictive equation set.  These substitutions allow
     510             :   ! some terms to solved implicitly or semi-implicitly in the predictive
     511             :   ! equations.
     512             :   !
     513             :   !
     514             :   !=============================================================================
     515           0 :   subroutine calc_setter_parameters( xm, xp2, Skx, sgn_wpxp,        & ! In
     516             :                                      big_L_x_1, big_L_x_2,          & ! In
     517             :                                      mu_x_1, mu_x_2, sigma_x_1_sqd, & ! Out
     518             :                                      sigma_x_2_sqd, mixt_frac,      & ! Out
     519             :                                      coef_sigma_x_1_sqd,            & ! Out
     520             :                                      coef_sigma_x_2_sqd             ) ! Out
     521             : 
     522             :     ! Description:
     523             :     ! Calculates the PDF component means, the PDF component standard deviations,
     524             :     ! and the mixture fraction for the variable that sets the PDF.
     525             : 
     526             :     ! References:
     527             :     !-----------------------------------------------------------------------
     528             : 
     529             :     use constants_clubb, only: &
     530             :         four,     & ! Variable(s)
     531             :         three,    &
     532             :         one,      &
     533             :         one_half, &
     534             :         zero,     &
     535             :         eps
     536             : 
     537             :     use clubb_precision, only: &
     538             :         core_rknd    ! Variable(s)
     539             : 
     540             :     implicit none
     541             : 
     542             :     ! Input Variables
     543             :     real( kind = core_rknd ), intent(in) :: &
     544             :       xm,        & ! Mean of x (overall)                            [units vary]
     545             :       xp2,       & ! Variance of x (overall)                    [(units vary)^2]
     546             :       Skx,       & ! Skewness of x                                           [-]
     547             :       sgn_wpxp,  & ! Sign of the covariance of w and x (overall)             [-]
     548             :       big_L_x_1, & ! Parameter for the spread of the 1st PDF comp. mean of x [-]
     549             :       big_L_x_2    ! Parameter for the spread of the 2nd PDF comp. mean of x [-]
     550             : 
     551             :     ! Output Variables
     552             :     real( kind = core_rknd ), intent(out) :: &
     553             :       mu_x_1,        & ! Mean of x (1st PDF component)              [units vary]
     554             :       mu_x_2,        & ! Mean of x (2nd PDF component)              [units vary]
     555             :       sigma_x_1_sqd, & ! Variance of x (1st PDF component)      [(units vary)^2]
     556             :       sigma_x_2_sqd, & ! Variance of x (2nd PDF component)      [(units vary)^2]
     557             :       mixt_frac        ! Mixture fraction                b                   [-]
     558             : 
     559             :     real( kind = core_rknd ), intent(out) :: &
     560             :       coef_sigma_x_1_sqd, & ! sigma_x_1^2 = coef_sigma_x_1_sqd * <x'^2>      [-]
     561             :       coef_sigma_x_2_sqd    ! sigma_x_2^2 = coef_sigma_x_2_sqd * <x'^2>      [-]
     562             : 
     563             :     ! Local Variables
     564             :     real( kind = core_rknd ) :: &
     565             :       mu_x_1_nrmlized, & ! Normalized mean of x (1st PDF component)          [-]
     566             :       mu_x_2_nrmlized    ! Normalized mean of x (2nd PDF component)          [-]
     567             : 
     568             :     real( kind = core_rknd ) :: &
     569             :       factor_plus,               &
     570             :       factor_minus,              &
     571             :       sqrt_factor_plus_ov_minus, &
     572             :       sqrt_factor_minus_ov_plus, &
     573             :       mu_x_1_nrmlized_thresh
     574             : 
     575             : 
     576             :     ! Calculate the factors in the PDF component mean equations.
     577           0 :     factor_plus = one + Skx * sgn_wpxp / sqrt( four + Skx**2 )
     578             : 
     579           0 :     factor_minus = one - Skx * sgn_wpxp / sqrt( four + Skx**2 )
     580             : 
     581           0 :     sqrt_factor_plus_ov_minus = sqrt( factor_plus / factor_minus )
     582           0 :     sqrt_factor_minus_ov_plus = sqrt( factor_minus / factor_plus )
     583             : 
     584             :     ! Calculate the normalized mean of x in the 1st PDF component.
     585           0 :     mu_x_1_nrmlized = big_L_x_1 * sqrt_factor_plus_ov_minus * sgn_wpxp
     586             : 
     587             :     ! Calculate the normalized mean of x in the 2nd PDF component.
     588           0 :     mu_x_2_nrmlized = -big_L_x_2 * sqrt_factor_minus_ov_plus * sgn_wpxp
     589             : 
     590             :     ! Calculate the mean of x in the 1st PDF component.
     591           0 :     mu_x_1 = xm + mu_x_1_nrmlized * sqrt( xp2 )
     592             : 
     593             :     ! Calculate the mean of x in the 2nd PDF component.
     594           0 :     mu_x_2 = xm + mu_x_2_nrmlized * sqrt( xp2 )
     595             : 
     596             :     ! Calculate the mixture fraction.
     597             :     if ( abs( mu_x_1_nrmlized ) >= eps &
     598           0 :          .and. abs( mu_x_2_nrmlized ) >= eps ) then
     599           0 :        mixt_frac = one / ( one + abs( mu_x_1_nrmlized / mu_x_2_nrmlized ) )
     600             :     elseif ( abs( mu_x_1_nrmlized ) >= eps &
     601           0 :              .and. abs( mu_x_2_nrmlized ) < eps ) then
     602           0 :        mixt_frac = one / ( one + abs( mu_x_1_nrmlized / eps ) )
     603             :     elseif ( abs( mu_x_1_nrmlized ) < eps &
     604           0 :              .and. abs( mu_x_2_nrmlized ) >= eps ) then
     605           0 :        mixt_frac = one / ( one + abs( eps / mu_x_2_nrmlized ) )
     606             :     else ! abs( mu_x_1_nrmlized ) < eps and abs( mu_x_2_nrmlized ) < eps
     607           0 :        mixt_frac = one_half
     608             :     endif
     609             : 
     610             :     ! Use a minimum magnitude value of mu_x_1_nrmlized in the denominator of a
     611             :     ! term in the PDF component variance equations in order to prevent a
     612             :     ! divide-by-zero error.
     613           0 :     if ( mu_x_1_nrmlized >= zero ) then
     614           0 :        mu_x_1_nrmlized_thresh = max( mu_x_1_nrmlized, eps )
     615             :     else ! mu_x_1_nrmlized < 0
     616           0 :        mu_x_1_nrmlized_thresh = min( mu_x_1_nrmlized, -eps )
     617             :     endif ! mu_x_1_nrmlized >= 0
     618             : 
     619             :     ! Calculate the variance of x in the 1st PDF component.
     620             :     coef_sigma_x_1_sqd &
     621             :     = one - mixt_frac * mu_x_1_nrmlized**2 &
     622             :       - ( one - mixt_frac ) * mu_x_2_nrmlized**2 &
     623             :       + ( one - mixt_frac ) &
     624             :         * ( Skx / ( three * mixt_frac * mu_x_1_nrmlized_thresh ) &
     625           0 :             - mu_x_1_nrmlized**2 / three + mu_x_2_nrmlized**2 / three )
     626             : 
     627           0 :     sigma_x_1_sqd = coef_sigma_x_1_sqd * xp2
     628             : 
     629             :     ! Calculate the variance of x in the 2nd PDF component.
     630             :     coef_sigma_x_2_sqd & 
     631             :     = one - mixt_frac * mu_x_1_nrmlized**2 &
     632             :       - ( one - mixt_frac ) * mu_x_2_nrmlized**2 &
     633             :       - mixt_frac &
     634             :         * ( Skx / ( three * mixt_frac * mu_x_1_nrmlized_thresh ) &
     635           0 :             - mu_x_1_nrmlized**2 / three + mu_x_2_nrmlized**2 / three )
     636             : 
     637           0 :     sigma_x_2_sqd = coef_sigma_x_2_sqd * xp2
     638             : 
     639             : 
     640           0 :     return
     641             : 
     642             :   end subroutine calc_setter_parameters
     643             : 
     644             :   !=============================================================================
     645           0 :   subroutine calc_L_x_Skx_fnc( Skx, sgn_wpxp,            & ! In
     646             :                                small_l_x_1, small_l_x_2, & ! In
     647             :                                big_L_x_1, big_L_x_2      ) ! Out
     648             : 
     649             :     ! Description:
     650             :     ! Calculates the values of big_L_x_1 and big_L_x_2 as functions of Skx.
     651             : 
     652             :     ! References:
     653             :     !-----------------------------------------------------------------------
     654             : 
     655             :     use constants_clubb, only: &
     656             :         four, & ! Variable(s)
     657             :         zero
     658             : 
     659             :     use clubb_precision, only: &
     660             :         core_rknd    ! Variable(s)
     661             : 
     662             :     implicit none
     663             : 
     664             :     ! Input Variables
     665             :     real( kind = core_rknd ), intent(in) :: &
     666             :       Skx,         & ! Skewness of x (overall)                               [-]
     667             :       sgn_wpxp,    & ! Sign of the covariance of w and x (overall)           [-]
     668             :       small_l_x_1, & ! Param. for the spread of the 1st PDF comp. mean of x  [-]
     669             :       small_l_x_2    ! Param. for the spread of the 2nd PDF comp. mean of x  [-]
     670             : 
     671             :     ! Output Variable
     672             :     real( kind = core_rknd ), intent(out) :: &
     673             :       big_L_x_1, & ! Parameter for the spread of the 1st PDF comp. mean of x [-]
     674             :       big_L_x_2    ! Parameter for the spread of the 2nd PDF comp. mean of x [-]
     675             : 
     676             :     ! Local Variable
     677             :     real( kind = core_rknd ) :: &
     678             :       Skx_fnc_factor
     679             : 
     680             : 
     681             :     ! The values of L_x_1 and L_x_2 are calculated by skewness functions.
     682             :     ! Those functions are:
     683             :     !
     684             :     ! L_x_1 = l_x_1 * abs( Skx ) / sqrt( 4 + Skx^2 ); and
     685             :     ! L_x_2 = l_x_2 * abs( Skx ) / sqrt( 4 + Skx^2 ).
     686             :     !
     687             :     ! The conditions on l_x_1 and l_x_2 are:
     688             :     !
     689             :     ! 2/3 < l_x_1 < 1 and 0 < l_x_2 < 1; when Skx * sgn( <w'x'> ) >= 0; and
     690             :     ! 0 < l_x_1 < 1 and 2/3 < l_x_2 < 1; when Skx * sgn( <w'x'> ) < 0.
     691             :     !
     692             :     ! For simplicity, this can also be accomplished by setting 2/3 < l_x_1 < 1
     693             :     ! and 0 < l_x_2 < 1, and then using the following equations.
     694             :     !
     695             :     ! When Skx * sgn( <w'x'> ) >= 0:
     696             :     ! L_x_1 = l_x_1 * abs( Skx ) / sqrt( 4 + Skx^2 ); and
     697             :     ! L_x_2 = l_x_2 * abs( Skx ) / sqrt( 4 + Skx^2 );
     698             :     !
     699             :     ! otherwise, when Skx * sgn( <w'x'> ) < 0, switch l_x_1 and l_x_2:
     700             :     ! L_x_1 = l_x_2 * abs( Skx ) / sqrt( 4 + Skx^2 ); and
     701             :     ! L_x_2 = l_x_1 * abs( Skx ) / sqrt( 4 + Skx^2 ).
     702           0 :     Skx_fnc_factor = abs( Skx ) / sqrt( four + Skx**2 )
     703             : 
     704           0 :     if ( Skx * sgn_wpxp >= zero ) then
     705           0 :        big_L_x_1 = small_l_x_1 * Skx_fnc_factor
     706           0 :        big_L_x_2 = small_l_x_2 * Skx_fnc_factor
     707             :     else ! Skx * sgn_wpxp < 0
     708           0 :        big_L_x_1 = small_l_x_2 * Skx_fnc_factor
     709           0 :        big_L_x_2 = small_l_x_1 * Skx_fnc_factor
     710             :     endif
     711             : 
     712             : 
     713           0 :     return
     714             : 
     715             :   end subroutine calc_L_x_Skx_fnc
     716             : 
     717             :   !=============================================================================
     718             :   !
     719             :   ! DESCRIPTION OF THE METHOD FOR EACH RESPONDING VARIABLE
     720             :   ! ======================================================
     721             :   !
     722             :   ! In order to find equations for the four PDF parameters for each responding
     723             :   ! variable, which are mu_x_1, mu_x_2, sigma_x_1, and sigma_x_2 (where x stands
     724             :   ! for a responding variable here), four equations are needed.  These four
     725             :   ! equations are the equations for <x>, <x'^2>, and <x'^3> as found by
     726             :   ! integrating over the PDF.  Additionally, one more equation, which involves
     727             :   ! tunable parameter L_x_1, and which is used to help control the mean of the
     728             :   ! 1st PDF component, is used in this equation set.  The four equations are:
     729             :   !
     730             :   ! <x> = mixt_frac * mu_x_1 + ( 1 - mixt_frac ) * mu_x_2;
     731             :   !
     732             :   ! <x'^2> = mixt_frac * ( ( mu_x_1 - <x> )^2 + sigma_x_1^2 )
     733             :   !          + ( 1 - mixt_frac ) * ( ( mu_x_2 - <x> )^2 + sigma_x_2^2 );
     734             :   !
     735             :   ! <x'^3> = mixt_frac * ( mu_x_1 - <x> )
     736             :   !                    * ( ( mu_x_1 - <x> )^2 + 3 * sigma_x_1^2 )
     737             :   !          + ( 1 - mixt_frac ) * ( mu_x_2 - <x> )
     738             :   !                              * ( ( mu_x_2 - <x> )^2 + 3 * sigma_x_2^2 ); and
     739             :   !
     740             :   ! mu_x_1 - <x> = L_x_1
     741             :   !                * sqrt( ( 1 + Skx * sgn( <w'x'> ) / sqrt( 4 + Skx^2 ) )
     742             :   !                        / ( 1 - Skx * sgn( <w'x'> ) / sqrt( 4 + Skx^2 ) ) )
     743             :   !                * sqrt( <x'^2> ) * sgn( <w'x'> );
     744             :   !
     745             :   ! where 0 <= L_x_1 <= 1, Skx is the skewness of x, such that
     746             :   ! Skx = <x'^3> / <x'^2>^(3/2), and sgn( <w'x'> ) is given by:
     747             :   !
     748             :   ! sgn( <w'x'> ) = |  1, when <w'x'> >= 0;
     749             :   !                 | -1, when <w'x'> < 0.
     750             :   !
     751             :   ! The resulting equations for the four PDF parameters are:
     752             :   !
     753             :   ! mu_x_1 = <x> + L_x_1
     754             :   !                * sqrt( ( 1 + Skx * sgn( <w'x'> ) / sqrt( 4 + Skx^2 ) )
     755             :   !                        / ( 1 - Skx * sgn( <w'x'> ) / sqrt( 4 + Skx^2 ) ) )
     756             :   !                * sqrt( <x'^2> ) * sgn( <w'x'> );
     757             :   !
     758             :   ! mu_x_2 = <x> - ( mixt_frac / ( 1 - mixt_frac ) ) * ( mu_x_1 - <x> );
     759             :   !
     760             :   ! sigma_x_1 = sqrt( ( 1 - mixt_frac * mu_x_1_nrmlized^2
     761             :   !                     - ( 1 - mixt_frac ) * mu_x_2_nrmlized^2
     762             :   !                     + ( 1 - mixt_frac )
     763             :   !                       * ( Skx / ( 3 * mixt_frac * mu_x_1_nrmlized )
     764             :   !                           - mu_x_1_nrmlized^2 / 3
     765             :   !                           + mu_x_2_nrmlized^2 / 3 ) )
     766             :   !                   * <x'^2> ); and
     767             :   !
     768             :   ! sigma_x_2 = sqrt( ( 1 - mixt_frac * mu_x_1_nrmlized^2
     769             :   !                     - ( 1 - mixt_frac ) * mu_x_2_nrmlized^2
     770             :   !                     - mixt_frac
     771             :   !                       * ( Skx / ( 3 * mixt_frac * mu_x_1_nrmlized )
     772             :   !                           - mu_x_1_nrmlized^2 / 3
     773             :   !                           + mu_x_2_nrmlized^2 / 3 ) )
     774             :   !                   * <x'^2> ); where
     775             :   !
     776             :   ! mu_x_1_nrmlized = ( mu_x_1 - <x> ) / sqrt( <x'^2> ); and
     777             :   !
     778             :   ! mu_x_2_nrmlized = ( mu_x_2 - <x> ) / sqrt( <x'^2> ).
     779             :   !
     780             :   !
     781             :   ! Notes:
     782             :   !
     783             :   ! When L_x_1 = 0, mu_x_1 = mu_x_2 = <x>, sigma_x_1^2 = sigma_x_2^2 = <x'^2>,
     784             :   ! and the distribution reduces to a single Gaussian.
     785             :   !
     786             :   !
     787             :   ! Equations for PDF component standard deviations:
     788             :   !
     789             :   ! The equations for the PDF component standard deviations can also be written
     790             :   ! as:
     791             :   !
     792             :   ! sigma_x_1 = sqrt( coef_sigma_x_1_sqd * <x'^2> ); and
     793             :   !
     794             :   ! sigma_x_2 = sqrt( coef_sigma_x_2_sqd * <x'^2> ); where
     795             :   !
     796             :   ! coef_sigma_x_1_sqd = 1 - mixt_frac * mu_x_1_nrmlized^2
     797             :   !                      - ( 1 - mixt_frac ) * mu_x_2_nrmlized^2
     798             :   !                      + ( 1 - mixt_frac )
     799             :   !                        * ( Skx / ( 3 * mixt_frac * mu_x_1_nrmlized )
     800             :   !                            - mu_x_1_nrmlized^2 / 3
     801             :   !                            + mu_x_2_nrmlized^2 / 3 ); and
     802             :   !
     803             :   ! coef_sigma_x_2_sqd = 1 - mixt_frac * mu_x_1_nrmlized^2
     804             :   !                      - ( 1 - mixt_frac ) * mu_x_2_nrmlized^2
     805             :   !                      - mixt_frac
     806             :   !                        * ( Skx / ( 3 * mixt_frac * mu_x_1_nrmlized )
     807             :   !                            - mu_x_1_nrmlized^2 / 3
     808             :   !                            + mu_x_2_nrmlized^2 / 3 ).
     809             :   !
     810             :   ! The above equations can be substituted into an equation for a variable that
     811             :   ! has been derived by integrating over the PDF.  Many variables like this are
     812             :   ! used in parts of the predictive equation set.  These substitutions allow
     813             :   ! some terms to solved implicitly or semi-implicitly in the predictive
     814             :   ! equations.
     815             :   !
     816             :   !
     817             :   !=============================================================================
     818           0 :   subroutine calc_respnder_parameters( xm, xp2, Skx, sgn_wpxp,       & ! In
     819             :                                        mixt_frac, big_L_x_1,         & ! In
     820             :                                        mu_x_1, mu_x_2,               & ! Out
     821             :                                        sigma_x_1_sqd, sigma_x_2_sqd, & ! Out
     822             :                                        coef_sigma_x_1_sqd,           & ! Out
     823             :                                        coef_sigma_x_2_sqd            ) ! Out
     824             : 
     825             :     ! Description:
     826             :     ! Calculates the PDF component means, the PDF component standard deviations,
     827             :     ! and the mixture fraction for the variable that sets the PDF.
     828             : 
     829             :     ! References:
     830             :     !-----------------------------------------------------------------------
     831             : 
     832             :     use constants_clubb, only: &
     833             :         four,  & ! Variable(s)
     834             :         three, &
     835             :         one,   &
     836             :         zero,  &
     837             :         eps
     838             : 
     839             :     use clubb_precision, only: &
     840             :         core_rknd    ! Variable(s)
     841             : 
     842             :     implicit none
     843             : 
     844             :     ! Input Variables
     845             :     real( kind = core_rknd ), intent(in) :: &
     846             :       xm,        & ! Mean of x (overall)                            [units vary]
     847             :       xp2,       & ! Variance of x (overall)                    [(units vary)^2]
     848             :       Skx,       & ! Skewness of x                                           [-]
     849             :       sgn_wpxp,  & ! Sign of the covariance of w and x (overall)             [-]
     850             :       mixt_frac, & ! Mixture fraction                                        [-]
     851             :       big_L_x_1    ! Parameter for the spread of the 1st PDF comp. mean of x [-]
     852             : 
     853             :     ! Output Variables
     854             :     real( kind = core_rknd ), intent(out) :: &
     855             :       mu_x_1,        & ! Mean of x (1st PDF component)              [units vary]
     856             :       mu_x_2,        & ! Mean of x (2nd PDF component)              [units vary]
     857             :       sigma_x_1_sqd, & ! Variance of x (1st PDF component)      [(units vary)^2]
     858             :       sigma_x_2_sqd    ! Variance of x (2nd PDF component)      [(units vary)^2]
     859             : 
     860             :     real( kind = core_rknd ), intent(out) :: &
     861             :       coef_sigma_x_1_sqd, & ! sigma_x_1^2 = coef_sigma_x_1_sqd * <x'^2>      [-]
     862             :       coef_sigma_x_2_sqd    ! sigma_x_2^2 = coef_sigma_x_2_sqd * <x'^2>      [-]
     863             : 
     864             :     ! Local Variables
     865             :     real( kind = core_rknd ) :: &
     866             :       mu_x_1_nrmlized, & ! Normalized mean of x (1st PDF component)          [-]
     867             :       mu_x_2_nrmlized    ! Normalized mean of x (2nd PDF component)          [-]
     868             : 
     869             :     real( kind = core_rknd ) :: &
     870             :       factor_plus,               &
     871             :       factor_minus,              &
     872             :       sqrt_factor_plus_ov_minus, &
     873             :       mu_x_1_nrmlized_thresh
     874             : 
     875             : 
     876             :     ! Calculate the factors in the PDF component mean equations.
     877           0 :     factor_plus = one + Skx * sgn_wpxp / sqrt( four + Skx**2 )
     878             : 
     879           0 :     factor_minus = one - Skx * sgn_wpxp / sqrt( four + Skx**2 )
     880             : 
     881           0 :     sqrt_factor_plus_ov_minus = sqrt( factor_plus / factor_minus )
     882             : 
     883             :     ! Calculate the normalized mean of x in the 1st PDF component.
     884           0 :     mu_x_1_nrmlized = big_L_x_1 * sqrt_factor_plus_ov_minus * sgn_wpxp
     885             : 
     886             :     ! Calculate the normalized mean of x in the 2nd PDF component.
     887           0 :     mu_x_2_nrmlized = - ( mixt_frac / ( one - mixt_frac ) ) * mu_x_1_nrmlized
     888             : 
     889             :     ! Calculate the mean of x in the 1st PDF component.
     890           0 :     mu_x_1 = xm + mu_x_1_nrmlized * sqrt( xp2 )
     891             : 
     892             :     ! Calculate the mean of x in the 2nd PDF component.
     893           0 :     mu_x_2 = xm + mu_x_2_nrmlized * sqrt( xp2 )
     894             : 
     895             :     ! Use a minimum magnitude value of mu_x_1_nrmlized in the denominator of a
     896             :     ! term in the PDF component variance equations in order to prevent a
     897             :     ! divide-by-zero error.
     898           0 :     if ( mu_x_1_nrmlized >= zero ) then
     899           0 :        mu_x_1_nrmlized_thresh = max( mu_x_1_nrmlized, eps )
     900             :     else ! mu_x_1_nrmlized < 0
     901           0 :        mu_x_1_nrmlized_thresh = min( mu_x_1_nrmlized, -eps )
     902             :     endif ! mu_x_1_nrmlized >= 0
     903             : 
     904             :     ! Calculate the variance of x in the 1st PDF component.
     905             :     coef_sigma_x_1_sqd &
     906             :     = one - mixt_frac * mu_x_1_nrmlized**2 &
     907             :       - ( one - mixt_frac ) * mu_x_2_nrmlized**2 &
     908             :       + ( one - mixt_frac ) &
     909             :         * ( Skx / ( three * mixt_frac * mu_x_1_nrmlized_thresh ) &
     910           0 :             - mu_x_1_nrmlized**2 / three + mu_x_2_nrmlized**2 / three )
     911             : 
     912           0 :     sigma_x_1_sqd = coef_sigma_x_1_sqd * xp2
     913             : 
     914             :     ! Calculate the variance of x in the 2nd PDF component.
     915             :     coef_sigma_x_2_sqd & 
     916             :     = one - mixt_frac * mu_x_1_nrmlized**2 &
     917             :       - ( one - mixt_frac ) * mu_x_2_nrmlized**2 &
     918             :       - mixt_frac &
     919             :         * ( Skx / ( three * mixt_frac * mu_x_1_nrmlized_thresh ) &
     920           0 :             - mu_x_1_nrmlized**2 / three + mu_x_2_nrmlized**2 / three )
     921             : 
     922           0 :     sigma_x_2_sqd = coef_sigma_x_2_sqd * xp2
     923             : 
     924             : 
     925           0 :     return
     926             : 
     927             :   end subroutine calc_respnder_parameters
     928             : 
     929             :   !=============================================================================
     930             : 
     931             : end module new_tsdadg_pdf

Generated by: LCOV version 1.14