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