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

              Line data    Source code
           1             : ! $Id$
           2             : !===============================================================================
           3             : module new_hybrid_pdf
           4             :
           5             :   ! Description:
           6             :   ! The portion of CLUBB's multivariate, two-component PDF that is the
           7             :   ! multivariate, two-component normal PDF of vertical velocity (w), total water
           8             :   ! mixing ratio (rt), liquid water potential temperature (thl), and optionally,
           9             :   ! the west-east horizontal wind component (u), the south-north horizontal wind
          10             :   ! component (v), and passive scalars (sclr).
          11             :
          12             :   ! References:
          13             :   ! Griffin and Larson (2020)
          14             :   !-------------------------------------------------------------------------
          15             :
          16             :   implicit none
          17             :
          18             :   public :: calculate_mixture_fraction,  & ! Procedure(s)
          19             :             calculate_w_params,          &
          20             :             calculate_responder_params,  &
          21             :             calculate_coef_wp4_implicit, &
          22             :             calc_coef_wp2xp_implicit,    &
          23             :             calc_coefs_wpxp2_semiimpl,   &
          24             :             calc_coefs_wpxpyp_semiimpl
          25             :
          26             :   private
          27             :
          28             :   contains
          29             :
          30             :   !=============================================================================
          31             :   !
          32             :   ! DESCRIPTION OF THE METHOD FOR THE VARIABLE THAT SETS THE MIXTURE FRACTION
          33             :   ! =========================================================================
          34             :   !
          35             :   ! The variable that sets the mixture fraction for the PDF is w.  There are
          36             :   ! five PDF parameters that need to be calculated, which are mu_w_1 (the mean
          37             :   ! of w is the 1st PDF component), mu_w_2 (the mean of w in the 2nd PDF
          38             :   ! component), sigma_w_1 (the standard deviation of w in the 1st PDF
          39             :   ! component), sigma_w_2 (the standard deviation of w in the 2nd PDF
          40             :   ! component), and mixt_frac (the mixture fraction, which is the weight of the
          41             :   ! 1st PDF component).  In order to solve for these five parameters, five
          42             :   ! equations are needed.  These five equations are the equations for <w>,
          43             :   ! <w'^2>, and <w'^3> as found by integrating over the PDF.  Additionally, two
          44             :   ! more equations, which involve tunable parameters F_w and zeta_w, and which
          45             :   ! are used to help control the spread of the PDF component means and the size
          46             :   ! of the PDF component standard deviations compared to each other,
          47             :   ! respectively, are used in this equation set.  The five equations are:
          48             :   !
          49             :   ! <w> = mixt_frac * mu_w_1 + ( 1 - mixt_frac ) * mu_w_2;
          50             :   !
          51             :   ! <w'^2> = mixt_frac * ( ( mu_w_1 - <w> )^2 + sigma_w_1^2 )
          52             :   !          + ( 1 - mixt_frac ) * ( ( mu_w_2 - <w> )^2 + sigma_w_2^2 );
          53             :   !
          54             :   ! <w'^3> = mixt_frac * ( mu_w_1 - <w> )
          55             :   !                    * ( ( mu_w_1 - <w> )^2 + 3 * sigma_w_1^2 )
          56             :   !          + ( 1 - mixt_frac ) * ( mu_w_2 - <w> )
          57             :   !                              * ( ( mu_w_2 - <w> )^2 + 3 * sigma_w_2^2 );
          58             :   !
          59             :   ! mu_w_1 - <w> = sqrt(F_w) * ( sqrt( 1 - mixt_frac ) / sqrt( mixt_frac ) )
          60             :   !                * sqrt( <w'^2> );
          61             :   !
          62             :   ! where 0 <= F_w <= 1; and
          63             :   !
          64             :   ! 1 + zeta_w = ( mixt_frac * sigma_w_1^2 )
          65             :   !              / ( ( 1 - mixt_frac ) * sigma_w_2^2 );
          66             :   !
          67             :   ! where zeta_w > -1.
          68             :   !
          69             :   ! Following convention for w, mu_w_1 is defined to be greater than or equal to
          70             :   ! mu_w_2 (and is also greater than or equal to <w>, while mu_w_2 is less than
          71             :   ! or equal to <w>).  This is relationship is found in the mu_w_1 - <w>
          72             :   ! equation above.
          73             :   !
          74             :   ! The resulting equations for the five PDF parameters are:
          75             :   !
          76             :   ! mixt_frac
          77             :   ! = ( 4 * F_w^3
          78             :   !     + 18 * F_w * ( zeta_w + 1 ) * ( 1 - F_w ) / ( zeta_w + 2 )
          79             :   !     + 6 * F_w^2 * ( 1 - F_w ) / ( zeta_w + 2 )
          80             :   !     + Skw^2
          81             :   !     - Skw * sqrt( 4 * F_w^3
          82             :   !                   + 12 * F_w^2 * ( 1 - F_w )
          83             :   !                   + 36 * F_w * ( zeta_w + 1 ) * ( 1 - F_w )^2
          84             :   !                     / ( zeta_w + 2 )^2
          85             :   !                   + Skw^2 ) )
          86             :   !   / ( 2 * F_w * ( F_w - 3 )^2 + 2 * Skw^2 );
          87             :   !
          88             :   ! mu_w_1 = <w> + sqrt( F_w * ( ( 1 - mixt_frac ) / mixt_frac ) * <w'^2> );
          89             :   !
          90             :   ! mu_w_2 = <w> - ( mixt_frac / ( 1 - mixt_frac ) ) * ( mu_w_1 - <w> );
          91             :   !
          92             :   ! sigma_w_1 = sqrt( ( ( zeta_w + 1 ) * ( 1 - F_w ) )
          93             :   !                   / ( ( zeta_w + 2 ) * mixt_frac ) * <w'^2> ); and
          94             :   !
          95             :   ! sigma_w_2 = sqrt( ( mixt_frac * sigma_w_1^2 )
          96             :   !                   / ( ( 1 - mixt_frac ) * ( 1 + zeta_w ) ) );
          97             :   !
          98             :   ! where Skw is the skewness of w, and Skw = <w'^3> / <w'^2>^(3/2).
          99             :   !
         100             :   ! This method works for all values of F_w (where 0 <= F_w <= 1) and zeta_w
         101             :   ! (where zeta_w > -1).
         102             :   !
         103             :   !
         104             :   ! Special case:
         105             :   !
         106             :   ! When Skw = 0 and F_w = 0, the equation for mixt_frac is undefined.  The
         107             :   ! equation for mixture fraction in this scenario can be derived by using the
         108             :   ! above equation for mixture fraction and then setting Skw = 0.  The resulting
         109             :   ! equation becomes:
         110             :   !
         111             :   ! mixt_frac
         112             :   ! = ( 4 * F_w^3
         113             :   !     + 18 * F_w * ( zeta_w + 1 ) * ( 1 - F_w ) / ( zeta_w + 2 )
         114             :   !     + 6 * F_w^2 * ( 1 - F_w ) / ( zeta_w + 2 ) )
         115             :   !   / ( 2 * F_w * ( F_w - 3 )^2 ).
         116             :   !
         117             :   ! All of the terms in the numerator and denominator contain a F_w, so this
         118             :   ! equation can be rewritten as:
         119             :   !
         120             :   ! mixt_frac
         121             :   ! = ( 4 * F_w^2
         122             :   !     + 18 * ( zeta_w + 1 ) * ( 1 - F_w ) / ( zeta_w + 2 )
         123             :   !     + 6 * F_w * ( 1 - F_w ) / ( zeta_w + 2 ) )
         124             :   !   / ( 2 * ( F_w - 3 )^2 ).
         125             :   !
         126             :   ! Now setting F_w = 0, the equation becomes:
         127             :   !
         128             :   ! mixt_frac = ( 18 * ( zeta_w + 1 ) / ( zeta_w + 2 ) ) / 18;
         129             :   !
         130             :   ! which can be rewritten as:
         131             :   !
         132             :   ! mixt_frac = ( zeta_w + 1 ) / ( zeta_w + 2 ).
         133             :   !
         134             :   ! When F_w = 0, Skw must have a value of 0 in order for the PDF to function
         135             :   ! correctly.  When F_w = 0, mu_w_1 = mu_w_2.  When the two PDF component means
         136             :   ! are equal to each other (and to the overall mean, <w>), the only value of
         137             :   ! Skw that can be represented is a value of 0.  In the equation for mixture
         138             :   ! fraction, when F_w is set to 0, but | Skw | > 0, mixt_frac will either have
         139             :   ! a value of 0 or 1, depending on whether Skw is positive or negative,
         140             :   ! respectively.
         141             :   !
         142             :   ! The value of F_w should be set as a function of Skw.  The value F_w should
         143             :   ! go toward 0 as | Skw | (or Skw^2) goes toward 0.  The value of F_w should
         144             :   ! go toward 1 as | Skw | (or Skw^2) goes to infinity.
         145             :   !
         146             :   !
         147             :   ! Tunable parameters:
         148             :   !
         149             :   ! 1) F_w:  This parameter controls the spread of the PDF component means.  The
         150             :   !          range of this parameter is 0 <= F_w <= 1.  When F_w = 0, the two
         151             :   !          PDF component means (mu_w_1 and mu_w_2) are equal to each other
         152             :   !          (and Skw must equal 0).  All of the variance of w is accounted for
         153             :   !          by the PDF component standard deviations (sigma_w_1 and sigma_w_2).
         154             :   !          When F_w = 1, mu_w_1 and mu_w_2 are spread as far apart as they can
         155             :   !          be.  Both PDF component standard deviations (sigma_w_1 and
         156             :   !          sigma_w_2) are equal to 0, and all of the variance of w is
         157             :   !          accounted for by the spread of the PDF component means.
         158             :   !
         159             :   !          When sigma_w_1 = sigma_w_2 = 0, the equation for <w'^2> becomes:
         160             :   !
         161             :   !          <w'^2> = mixt_frac * ( mu_w_1 - <w> )^2
         162             :   !                   + ( 1 - mixt_frac ) * ( mu_w_2 - <w> )^2.
         163             :   !
         164             :   !          Substituting the equation for <w> into the above equation for
         165             :   !          mu_w_2 - <w>, the above equation becomes:
         166             :   !
         167             :   !          <w'^2> = ( mixt_frac / ( 1 - mixt_frac ) ) * ( mu_w_1 - <w> )^2;
         168             :   !
         169             :   !          which can be rewritten as:
         170             :   !
         171             :   !          ( mu_w_1 - <w> )^2 = ( ( 1 - mixt_frac ) / mixt_frac ) * <w'^2>.
         172             :   !
         173             :   !          Taking the square root of the above equation:
         174             :   !
         175             :   !          mu_w_1 - <w> = +/- ( sqrt( 1 - mixt_frac ) / sqrt(mixt_frac) )
         176             :   !                             * sqrt( <w'^2> ).
         177             :   !
         178             :   !          This equation can be compared to the equation for mu_w_1 - <w> in
         179             :   !          the set of 5 equations, which is:
         180             :   !
         181             :   !          mu_w_1 - <w>
         182             :   !          = sqrt(F_w) * ( sqrt( 1 - mixt_frac ) / sqrt( mixt_frac ) )
         183             :   !            * sqrt( <w'^2> ).
         184             :   !
         185             :   !          The above equations give another example of the meaning of F_w.
         186             :   !          The value of sqrt(F_w) is ratio of mu_w_1 - <w> to its maximum
         187             :   !          value, which is:
         188             :   !
         189             :   !          sqrt( ( ( 1 - mixt_frac ) / mixt_frac ) * <w'^2> ).
         190             :   !
         191             :   !
         192             :   ! 2) zeta_w:  This parameter controls the size of the PDF component standard
         193             :   !             deviations compared to each other.  The equation for zeta_w is:
         194             :   !
         195             :   !             1 + zeta_w = ( mixt_frac * sigma_w_1^2 )
         196             :   !                          / ( ( 1 - mixt_frac ) * sigma_w_2^2 ).
         197             :   !
         198             :   !             When zeta_w > 0, mixt_frac * sigma_w_1^2 increases at the
         199             :   !             expense of ( 1 - mixt_frac ) * sigma_w_2^2, which decreases in
         200             :   !             this variance-preserving equation set.  When zeta_w = 0, then
         201             :   !             mixt_frac * sigma_w_1^2 = ( 1 - mixt_frac ) * sigma_w_2^2.
         202             :   !             When -1 < zeta_w < 0, ( 1 - mixt_frac ) * sigma_w_2^2 increases
         203             :   !             at the expense of mixt_frac * sigma_w_1^2, which decreases.  As
         204             :   !             a result, greater values of zeta_w cause the 1st PDF component
         205             :   !             to become broader while the 2nd PDF component becomes narrower,
         206             :   !             and smaller values of zeta_w cause the 1st PDF component to
         207             :   !             become narrower while the 2nd PDF component becomes broader.
         208             :   !
         209             :   !             Symmetry
         210             :   !
         211             :   !             When zeta_w = 0, the PDF is always symmetric.  In other words,
         212             :   !             the PDF at any positive value of Skw (for example, Skw = 2.5)
         213             :   !             will look like a mirror-image (reflection across the y-axis)
         214             :   !             of the PDF at a negative value of Skw of the same magnitude (in
         215             :   !             this example, Skw = -2.5).  However, when zeta_w /= 0, the PDF
         216             :   !             loses this quality and is not symmetric.
         217             :   !
         218             :   !             When symmetry is desired at values of zeta_w besides zeta_w = 0,
         219             :   !             the solution is to turn zeta_w into a function of Skw.  A basic
         220             :   !             example of a zeta_w skewness equation that produces a symmetric
         221             :   !             PDF for values of zeta_w other than 0 is:
         222             :   !
         223             :   !             zeta_w = | zeta_w_in,                      when Skw >= 0;
         224             :   !                      | ( 1 / ( 1 + zeta_w_in ) ) - 1,  when Skw < 0.
         225             :   !
         226             :   !
         227             :   ! Notes:
         228             :   !
         229             :   ! When F_w = 0 (which can only happen when Skw = 0), mu_w_1 = mu_w_2, and
         230             :   ! mixt_frac = ( zeta_w + 1 ) / ( zeta_w + 2 ).  When these equations are
         231             :   ! substituted into the equations for sigma_w_1 and sigma_w_2, the result is
         232             :   ! sigma_w_1 = sigma_w_2 = sqrt( <w'^2> ).  This means that the distribution
         233             :   ! becomes a single Gaussian when F_w = 0 (and Skw = 0).  This happens
         234             :   ! regardless of the value of zeta_w.
         235             :   !
         236             :   ! The equations for the PDF component means and standard deviations can also
         237             :   ! be written as:
         238             :   !
         239             :   ! mu_w_1 = <w> + sqrt( F_w * ( ( 1 - mixt_frac ) / mixt_frac ) * <w'^2> );
         240             :   !
         241             :   ! mu_w_2 = <w> - sqrt( F_w * ( mixt_frac / ( 1 - mixt_frac ) ) * <w'^2> );
         242             :   !
         243             :   ! sigma_w_1 = sqrt( coef_sigma_w_1_sqd * <w'^2> ); and
         244             :   !
         245             :   ! sigma_w_2 = sqrt( coef_sigma_w_2_sqd * <w'^2> ); where
         246             :   !
         247             :   ! coef_sigma_w_1_sqd = ( ( zeta_w + 1 ) * ( 1 - F_w ) )
         248             :   !                      / ( ( zeta_w + 2 ) * mixt_frac ); and
         249             :   !
         250             :   ! coef_sigma_w_2_sqd = ( 1 - F_w ) / ( ( zeta_w + 2 ) * ( 1 - mixt_frac ) ).
         251             :   !
         252             :   ! The above equations can be substituted into an equation for a variable that
         253             :   ! has been derived by integrating over the PDF.  Many variables like this are
         254             :   ! used in parts of the predictive equation set.  These substitutions allow
         255             :   ! some terms to solved implicitly or semi-implicitly in the predictive
         256             :   ! equations.
         257             :   !
         258             :   !
         259             :   ! Brian Griffin; September 2017.
         260             :   !
         261             :   !=============================================================================
         262           0 :   elemental function calculate_mixture_fraction( Skw, F_w, zeta_w ) &
         263             :   result( mixt_frac )
         264             :
         265             :     ! Description:
         266             :     ! Calculates mixture fraction.
         267             :
         268             :     ! References:
         269             :     ! Griffin and Larson (2020)
         270             :     !-----------------------------------------------------------------------
         271             :
         272             :     use constants_clubb, only: &
         273             :         thirty_six, & ! Constant(s)
         274             :         eighteen,   &
         275             :         twelve,     &
         276             :         six,        &
         277             :         four,       &
         278             :         three,      &
         279             :         two,        &
         280             :         one,        &
         281             :         zero
         282             :
         283             :     use clubb_precision, only: &
         284             :         core_rknd    ! Variable(s)
         285             :
         286             :     implicit none
         287             :
         288             :     ! Input Variables
         289             :     real( kind = core_rknd ), intent(in) :: &
         290             :       Skw,    & ! Skewness of w                                            [-]
         291             :       F_w,    & ! Parameter for the spread of the PDF component means of w [-]
         292             :       zeta_w    ! Parameter for the PDF component variances of w           [-]
         293             :
         294             :     ! Return Variable
         295             :     real( kind = core_rknd ) :: &
         296             :       mixt_frac    ! Mixture fraction    [-]
         297             :
         298             :
         299             :     ! Calculate mixture fraction, which is the weight of the 1st PDF component.
         300             :     ! The 2nd PDF component has a weight of 1 - mixt_frac.
         301           0 :     if ( F_w > zero ) then
         302             :
         303             :        mixt_frac &
         304             :        = ( four * F_w**3 &
         305             :            + eighteen * F_w &
         306             :              * ( zeta_w + one ) * ( one - F_w ) / ( zeta_w + two ) &
         307             :            + six * F_w**2 * ( one - F_w ) / ( zeta_w + two ) &
         308             :            + Skw**2 &
         309             :            - Skw * sqrt( four * F_w**3 &
         310             :                          + twelve * F_w**2 * ( one - F_w ) &
         311             :                          + thirty_six * F_w &
         312             :                            * ( zeta_w + one ) * ( one - F_w )**2 &
         313             :                            / ( zeta_w + two )**2 &
         314             :                          + Skw**2 ) ) &
         315           0 :          / ( two * F_w * ( F_w - three )**2 + two * Skw**2 )
         316             :
         317             :     else ! F_w = 0
         318             :
         319           0 :        if ( abs( Skw ) > zero ) then
         320             :
         321             :           ! When F_w = 0, | Skw | must have a value of 0.  In a scenario where
         322             :           ! F_w = 0 and | Skw | > 0, the mixture fraction (and the rest of the
         323             :           ! PDF parameters) can't be calculated.  Since mixture fraction can
         324             :           ! only have values 0 < mixt_frac < 1, set mixt_frac to -1 in this
         325             :           ! scenario.
         326             :           mixt_frac = -one
         327             :
         328             :        else ! Skw = 0
         329             :
         330           0 :           mixt_frac = ( zeta_w + one ) / ( zeta_w + two )
         331             :
         332             :        endif ! | Skw | > 0
         333             :
         334             :     endif ! F_w > 0
         335             :
         336             :
         337             :     return
         338             :
         339             :   end function calculate_mixture_fraction
         340             :
         341             :   !=============================================================================
         342           0 :   elemental subroutine calculate_w_params( wm, wp2, Skw, F_w, zeta_w, & ! In
         343             :                                            mu_w_1, mu_w_2, sigma_w_1, & ! Out
         344             :                                            sigma_w_2, mixt_frac,      & ! Out
         345             :                                            coef_sigma_w_1_sqd,        & ! Out
         346             :                                            coef_sigma_w_2_sqd         ) ! Out
         347             :
         348             :     ! Description:
         349             :     ! Calculates the PDF component means, the PDF component standard deviations,
         350             :     ! and the mixture fraction for the variable that sets the PDF.
         351             :
         352             :     ! References:
         353             :     ! Griffin and Larson (2020)
         354             :     !-----------------------------------------------------------------------
         355             :
         356             :     use constants_clubb, only: &
         357             :         two,  & ! Variable(s)
         358             :         one,  &
         359             :         zero
         360             :
         361             :     use clubb_precision, only: &
         362             :         core_rknd    ! Variable(s)
         363             :
         364             :     implicit none
         365             :
         366             :     ! Input Variables
         367             :     real( kind = core_rknd ), intent(in) :: &
         368             :       wm,     & ! Mean of w (overall)                                     [m/s]
         369             :       wp2,    & ! Variance of w (overall)                             [m^2/s^2]
         370             :       Skw,    & ! Skewness of w                                             [-]
         371             :       F_w,    & ! Parameter for the spread of the PDF component means of w  [-]
         372             :       zeta_w    ! Parameter for the PDF component variances of w            [-]
         373             :
         374             :     ! Output Variables
         375             :     real( kind = core_rknd ), intent(out) :: &
         376             :       mu_w_1,    & ! Mean of w (1st PDF component)                  [m/s]
         377             :       mu_w_2,    & ! Mean of w (2nd PDF component)                  [m/s]
         378             :       sigma_w_1, & ! Standard deviation of w (1st PDF component)    [m/s]
         379             :       sigma_w_2, & ! Standard deviation of w (2nd PDF component)    [m/s]
         380             :       mixt_frac    ! Mixture fraction                               [-]
         381             :
         382             :     real( kind = core_rknd ), intent(out) :: &
         383             :       coef_sigma_w_1_sqd, & ! sigma_w_1^2 = coef_sigma_w_1_sqd * <w'^2>     [-]
         384             :       coef_sigma_w_2_sqd    ! sigma_w_2^2 = coef_sigma_w_2_sqd * <w'^2>     [-]
         385             :
         386             :
         387             :     ! Calculate the mixture fraction.
         388           0 :     mixt_frac = calculate_mixture_fraction( Skw, F_w, zeta_w )
         389             :
         390           0 :     if ( mixt_frac > zero .and. mixt_frac < one ) then
         391             :
         392             :        ! Calculate the mean of w in the 1st PDF component.
         393           0 :        mu_w_1 = wm + sqrt( F_w * ( ( one - mixt_frac ) / mixt_frac ) * wp2 )
         394             :
         395             :        ! Calculate the mean of w in the 2nd PDF component.
         396           0 :        mu_w_2 = wm - ( mixt_frac / ( one - mixt_frac ) ) * ( mu_w_1 - wm )
         397             :
         398             :        ! Calculate the standard deviation of w in the 1st PDF component.
         399             :        ! sigma_w_1 = sqrt( ( ( zeta_w + 1 ) * ( 1 - F_w ) )
         400             :        !                   / ( ( zeta_w + 2 ) * mixt_frac ) * <w'^2> )
         401             :        coef_sigma_w_1_sqd = ( ( zeta_w + one ) * ( one - F_w ) ) &
         402           0 :                             / ( ( zeta_w + two ) * mixt_frac )
         403             :
         404           0 :        sigma_w_1 = sqrt( coef_sigma_w_1_sqd * wp2 )
         405             :
         406             :        ! Calculate the standard deviation of w in the 2nd PDF component.
         407             :        ! sigma_w_2 = sqrt( ( mixt_frac * sigma_w_1^2 )
         408             :        !                   / ( ( 1 - mixt_frac ) * ( 1 + zeta_w ) ) )
         409             :        !           = sqrt( ( 1 - F_w )
         410             :        !                   / ( ( zeta_w + 2 ) * ( 1 - mixt_frac ) ) * <w'^2> )
         411             :        coef_sigma_w_2_sqd = ( one - F_w ) &
         412           0 :                             / ( ( zeta_w + two ) * ( one - mixt_frac ) )
         413             :
         414           0 :        sigma_w_2 = sqrt( coef_sigma_w_2_sqd * wp2 )
         415             :
         416             :     else ! mixt_frac <= 0 or mixt_frac >= 1
         417             :
         418             :        ! The mixture fraction produced is invalid.  This should only happen in
         419             :        ! the scenario where F_w = 0 and | Skw | > 0, where the value of
         420             :        ! mixt_frac has been set to -1.  Set all output variables to 0 in this
         421             :        ! scenario.  Since F_w is a function of skewness, the mixture fraction
         422             :        ! and the PDF should always be valid, and this section of code shouldn't
         423             :        ! be entered.
         424           0 :        mu_w_1 = zero
         425           0 :        mu_w_2 = zero
         426           0 :        sigma_w_1 = zero
         427           0 :        sigma_w_2 = zero
         428           0 :        coef_sigma_w_1_sqd = zero
         429           0 :        coef_sigma_w_2_sqd = zero
         430             :
         431             :     endif ! 0 < mixt_frac < 1
         432             :
         433             :
         434           0 :     return
         435             :
         436             :   end subroutine calculate_w_params
         437             :
         438             :   !=============================================================================
         439             :   !
         440             :   ! DESCRIPTION OF THE METHOD FOR EACH RESPONDING VARIABLE
         441             :   ! ======================================================
         442             :   !
         443             :   ! In order to find equations for the four PDF parameters for each responding
         444             :   ! variable, which are mu_x_1, mu_x_2, sigma_x_1, and sigma_x_2 (where x stands
         445             :   ! for a responding variable here), four equations are needed.  These four
         446             :   ! equations are the equations for <x>, <x'^2>, <x'^3>, and <w'x'> as found by
         447             :   ! integrating over the PDF.  The four equations are:
         448             :   !
         449             :   ! <x> = mixt_frac * mu_x_1 + ( 1 - mixt_frac ) * mu_x_2;
         450             :   !
         451             :   ! <x'^2> = mixt_frac * ( ( mu_x_1 - <x> )^2 + sigma_x_1^2 )
         452             :   !          + ( 1 - mixt_frac ) * ( ( mu_x_2 - <x> )^2 + sigma_x_2^2 );
         453             :   !
         454             :   ! <x'^3> = mixt_frac * ( mu_x_1 - <x> )
         455             :   !                    * ( ( mu_x_1 - <x> )^2 + 3 * sigma_x_1^2 )
         456             :   !          + ( 1 - mixt_frac ) * ( mu_x_2 - <x> )
         457             :   !                              * ( ( mu_x_2 - <x> )^2 + 3 * sigma_x_2^2 ); and
         458             :   !
         459             :   ! <w'x'> = mixt_frac * ( mu_w_1 - <w> ) * ( mu_x_1 - <x> )
         460             :   !          + ( 1 - mixt_frac ) * ( mu_w_2 - <w> ) * ( mu_x_2 - <x> );
         461             :   !
         462             :   ! where the correlations that are normally found in the <w'x'> equation,
         463             :   ! corr_w_x_1 and corr_w_x_2, have both been set to 0.
         464             :   !
         465             :   ! The equations for mu_w_1 - <w> and mu_w_2 - <w> are:
         466             :   !
         467             :   ! mu_w_1 - <w> = sqrt( F_w * ( ( 1 - mixt_frac ) / mixt_frac ) * <w'^2> );
         468             :   !
         469             :   ! mu_w_2 - <w> = - sqrt( F_w * ( mixt_frac / ( 1 - mixt_frac ) ) * <w'^2> );
         470             :   !
         471             :   ! The resulting equations for the four PDF parameters are:
         472             :   !
         473             :   ! mu_x_1 = <x> + sqrt( ( 1 - mixt_frac ) / mixt_frac )
         474             :   !                * <w'x'> / sqrt( F_w * <w'^2> );
         475             :   !
         476             :   ! mu_x_2 = <x> - sqrt( mixt_frac / ( 1 - mixt_frac ) )
         477             :   !                * <w'x'> / sqrt( F_w * <w'^2> );
         478             :   !
         479             :   ! sigma_x_1^2 = ( 1 + sqrt( ( 1 - mixt_frac ) / mixt_frac )
         480             :   !                     * Skx * sqrt( F_w * <w'^2> * <x'^2> ) / ( 3 * <w'x'> )
         481             :   !                 - ( ( 1 + mixt_frac ) / mixt_frac )
         482             :   !                   * <w'x'>^2 / ( 3 * F_w * <w'^2> * <x'^2> ) )
         483             :   !               * <x'^2>; and
         484             :   !
         485             :   ! sigma_x_2^2 = ( 1 - sqrt( mixt_frac / ( 1 - mixt_frac ) )
         486             :   !                     * Skx * sqrt( F_w * <w'^2> * <x'^2> ) / ( 3 * <w'x'> )
         487             :   !                 + ( ( mixt_frac - 2 ) / ( 1 - mixt_frac ) )
         488             :   !                   * <w'x'>^2 / ( 3 * F_w * <w'^2> * <x'^2> ) )
         489             :   !               * <x'^2>;
         490             :   !
         491             :   ! where Skx is the skewness of x, and Skx = <x'^3> / <x'^2>^(3/2).
         492             :   !
         493             :   !
         494             :   ! Limits on F_w:
         495             :   !
         496             :   ! The only limits placed on the value of F_w from the w equation set itself
         497             :   ! are 0 <= F_w <= 1.  However, use of the above equation set for responder
         498             :   ! variable x forces an additional limit to be placed on the value of F_w.
         499             :   ! That additional limit restricts the range of F_w to:
         500             :   !
         501             :   ! <w'x'>^2 / ( <w'^2> * <x'^2> ) <= F_w <= 1.
         502             :   !
         503             :   ! Furthermore, when there is more than one responder variable, F_w is limited
         504             :   ! by the most restrictive cases, such that:
         505             :   !
         506             :   ! max( <w'x'>^2 / ( <w'^2> * <x'^2> ), for all variables x ) <= F_w <= 1.
         507             :   !
         508             :   !
         509             :   ! Limits on Skx:
         510             :   !
         511             :   ! Since the PDF parameters for this variable need to work with the mixture
         512             :   ! fraction that has been provided by the setting variable, which is w, the
         513             :   ! method does not work for all values of Skx.  However, the limits of Skx can
         514             :   ! always be calculated.  The limits on Skw given by:
         515             :   !
         516             :   ! when <w'x'> > 0:
         517             :   !
         518             :   ! ( 1 + mixt_frac ) / sqrt( mixt_frac * ( 1 - mixt_frac ) )
         519             :   ! * <w'x'>^3 / ( F_w * <w'^2> * <x'^2> )^(3/2)
         520             :   ! - sqrt( mixt_frac / ( 1 - mixt_frac ) )
         521             :   !   * 3 * <w'x'> / sqrt( F_w * <w'^2> * <x'^2> )
         522             :   ! <= Skx <=
         523             :   ! ( mixt_frac - 2 ) / sqrt( mixt_frac * ( 1 - mixt_frac ) )
         524             :   ! * <w'x'>^3 / ( F_w * <w'^2> * <x'^2> )^(3/2)
         525             :   ! + sqrt( ( 1 - mixt_frac ) / mixt_frac )
         526             :   !   * 3 * <w'x'> / sqrt( F_w * <w'^2> * <x'^2> );
         527             :   !
         528             :   ! when <w'x'> < 0:
         529             :   !
         530             :   ! ( mixt_frac - 2 ) / sqrt( mixt_frac * ( 1 - mixt_frac ) )
         531             :   ! * <w'x'>^3 / ( F_w * <w'^2> * <x'^2> )^(3/2)
         532             :   ! + sqrt( ( 1 - mixt_frac ) / mixt_frac )
         533             :   !   * 3 * <w'x'> / sqrt( F_w * <w'^2> * <x'^2> )
         534             :   ! <= Skx <=
         535             :   ! ( 1 + mixt_frac ) / sqrt( mixt_frac * ( 1 - mixt_frac ) )
         536             :   ! * <w'x'>^3 / ( F_w * <w'^2> * <x'^2> )^(3/2)
         537             :   ! - sqrt( mixt_frac / ( 1 - mixt_frac ) )
         538             :   !   * 3 * <w'x'> / sqrt( F_w * <w'^2> * <x'^2> );
         539             :   !
         540             :   ! and when <w'x'> = 0, Skx = 0.
         541             :   !
         542             :   !
         543             :   ! Special cases:
         544             :   !
         545             :   ! When <w'x'> = 0, mu_x_1 = mu_x_2 = <x>, and the value of Skx must be 0.
         546             :   ! Since both <w'x'> = 0 and Skx = 0, the equations for sigma_x_1^2 and
         547             :   ! sigma_x_2^2 are both undefined.  In this situation, the equations for the
         548             :   ! PDF parameters of x are:
         549             :   !
         550             :   ! mu_x_1 = mu_x_2 = <x>; and
         551             :   ! sigma_x_1^2 = sigma_x_2^2 = <x'^2>.
         552             :   !
         553             :   ! The value of F_w is allowed to be 0 only when <w'x'> = 0 (for all variables
         554             :   ! x).  When <w'^2> = 0 and/or <x'^2> = 0, this means that <w'x'> = 0, as well.
         555             :   ! In all these situations, the equation set for the situation when <w'x'> = 0
         556             :   ! is used.  This means that the distribution becomes a single Gaussian when
         557             :   ! <w'x'> = 0 (and Skx = 0).
         558             :   !
         559             :   !
         560             :   ! Notes:
         561             :   !
         562             :   ! The equations for the PDF component means and standard deviations can also
         563             :   ! be written as:
         564             :   !
         565             :   ! mu_x_1 = <x> + sqrt( ( 1 - mixt_frac ) / mixt_frac )
         566             :   !                * <w'x'> / sqrt( F_w * <w'^2> );
         567             :   !
         568             :   ! mu_x_2 = <x> - sqrt( mixt_frac / ( 1 - mixt_frac ) )
         569             :   !                * <w'x'> / sqrt( F_w * <w'^2> );
         570             :   !
         571             :   ! sigma_x_1 = sqrt( coef_sigma_x_1_sqd * <x'^2> ); and
         572             :   !
         573             :   ! sigma_x_2 = sqrt( coef_sigma_x_2_sqd * <x'^2> ); where
         574             :   !
         575             :   ! coef_sigma_x_1_sqd
         576             :   ! = 1 + sqrt( ( 1 - mixt_frac ) / mixt_frac )
         577             :   !       * Skx * sqrt( F_w * <w'^2> * <x'^2> ) / ( 3 * <w'x'> )
         578             :   !   - ( ( 1 + mixt_frac ) / mixt_frac )
         579             :   !     * <w'x'>^2 / ( 3 * F_w * <w'^2> * <x'^2> ); and
         580             :   !
         581             :   ! coef_sigma_x_2_sqd
         582             :   ! = 1 - sqrt( mixt_frac / ( 1 - mixt_frac ) )
         583             :   !       * Skx * sqrt( F_w * <w'^2> * <x'^2> ) / ( 3 * <w'x'> )
         584             :   !   + ( ( mixt_frac - 2 ) / ( 1 - mixt_frac ) )
         585             :   !     * <w'x'>^2 / ( 3 * F_w * <w'^2> * <x'^2> ).
         586             :   !
         587             :   ! The above equations can be substituted into an equation for a variable that
         588             :   ! has been derived by integrating over the PDF.  Many variables like this are
         589             :   ! used in parts of the predictive equation set.  These substitutions allow
         590             :   ! some terms to solved implicitly or semi-implicitly in the predictive
         591             :   ! equations.
         592             :   !
         593             :   !
         594             :   ! Brian Griffin; September 2019.
         595             :   !
         596             :   !=============================================================================
         597           0 :   elemental subroutine calculate_responder_params( xm, xp2, Skx, wpxp,  & ! In
         598             :                                                    wp2, F_w, mixt_frac, & ! In
         599             :                                                    mu_x_1, mu_x_2,      & ! Out
         600             :                                                    sigma_x_1_sqd,       & ! Out
         601             :                                                    sigma_x_2_sqd,       & ! Out
         602             :                                                    coef_sigma_x_1_sqd,  & ! Out
         603             :                                                    coef_sigma_x_2_sqd   ) ! Out
         604             :
         605             :     ! Description:
         606             :     ! Calculates the PDF component means and the PDF component standard
         607             :     ! deviations for a responding variable (a variable that is not used to set
         608             :     ! the mixture fraction).
         609             :
         610             :     ! References:
         611             :     ! Griffin and Larson (2020)
         612             :     !-----------------------------------------------------------------------
         613             :
         614             :     use constants_clubb, only: &
         615             :         three, & ! Variable(s)
         616             :         two,   &
         617             :         one,   &
         618             :         zero
         619             :
         620             :     use clubb_precision, only: &
         621             :         core_rknd    ! Variable(s)
         622             :
         623             :     implicit none
         624             :
         625             :     ! Input Variables
         626             :     real( kind = core_rknd ), intent(in) :: &
         627             :       xm,       & ! Mean of x (overall)                             [units vary]
         628             :       xp2,      & ! Variance of x (overall)                     [(units vary)^2]
         629             :       Skx,      & ! Skewness of x                                            [-]
         630             :       wpxp,     & ! Covariance of w and x (overall)            [m/s(units vary)]
         631             :       wp2,      & ! Variance of w (overall)                            [m^2/s^2]
         632             :       F_w,      & ! Parameter for the spread of the PDF component means of w [-]
         633             :       mixt_frac   ! Mixture fraction                                         [-]
         634             :
         635             :     ! Output Variables
         636             :     real( kind = core_rknd ), intent(out) :: &
         637             :       mu_x_1,        & ! Mean of x (1st PDF component)              [units vary]
         638             :       mu_x_2,        & ! Mean of x (2nd PDF component)              [units vary]
         639             :       sigma_x_1_sqd, & ! Variance of x (1st PDF component)      [(units vary)^2]
         640             :       sigma_x_2_sqd    ! Variance of x (2nd PDF component)      [(units vary)^2]
         641             :
         642             :     real( kind = core_rknd ), intent(out) :: &
         643             :       coef_sigma_x_1_sqd, & ! sigma_x_1^2 = coef_sigma_x_1_sqd * <x'^2>      [-]
         644             :       coef_sigma_x_2_sqd    ! sigma_x_2^2 = coef_sigma_x_2_sqd * <x'^2>      [-]
         645             :
         646             :
         647           0 :     if ( abs( wpxp ) > zero ) then
         648             :
         649             :        ! Note:  when |<w'x'>| > 0, F_w, <w'^2>, and <x'^2> must all have values
         650             :        !        greater than 0.
         651             :
         652             :        ! Calculate the mean of x in the 1st PDF component.
         653             :        mu_x_1 = xm + sqrt( ( one - mixt_frac ) / mixt_frac ) &
         654           0 :                      * wpxp / sqrt( F_w * wp2 )
         655             :
         656             :        ! Calculate the mean of x in the 2nd PDF component.
         657             :        mu_x_2 = xm - sqrt( mixt_frac / ( one - mixt_frac ) ) &
         658           0 :                      * wpxp / sqrt( F_w * wp2 )
         659             :
         660             :        ! Calculate the variance of x in the 1st PDF component.
         661             :        ! sigma_x_1^2
         662             :        ! = ( 1 + sqrt( ( 1 - mixt_frac ) / mixt_frac )
         663             :        !         * Skx * sqrt( F_w * <w'^2> * <x'^2> ) / ( 3 * <w'x'> )
         664             :        !     - ( ( 1 + mixt_frac ) / mixt_frac )
         665             :        !       * <w'x'>^2 / ( 3 * F_w * <w'^2> * <x'^2> ) )
         666             :        !   * <x'^2>
         667             :        coef_sigma_x_1_sqd &
         668             :        = one + sqrt( ( one - mixt_frac ) / mixt_frac ) &
         669             :                * Skx * sqrt( F_w * wp2 * xp2 ) / ( three * wpxp ) &
         670             :          - ( ( one + mixt_frac ) / mixt_frac ) &
         671           0 :            * wpxp**2 / ( three * F_w * wp2 * xp2 )
         672             :
         673             :        ! Mathematically, the value of coef_sigma_x_1_sqd cannot be less than 0.
         674             :        ! Numerically, this can happen when numerical round off error causes an
         675             :        ! epsilon-sized negative value.  When this happens, reset the value of
         676             :        ! coef_sigma_x_1_sqd to 0.
         677           0 :        coef_sigma_x_1_sqd = max( coef_sigma_x_1_sqd, zero )
         678             :
         679           0 :        sigma_x_1_sqd = coef_sigma_x_1_sqd * xp2
         680             :
         681             :        ! Calculate the variance of x in the 2nd PDF component.
         682             :        ! sigma_x_2^2
         683             :        ! = ( 1 - sqrt( mixt_frac / ( 1 - mixt_frac ) )
         684             :        !         * Skx * sqrt( F_w * <w'^2> * <x'^2> ) / ( 3 * <w'x'> )
         685             :        !     + ( ( mixt_frac - 2 ) / ( 1 - mixt_frac ) )
         686             :        !       * <w'x'>^2 / ( 3 * F_w * <w'^2> * <x'^2> ) )
         687             :        !   * <x'^2>
         688             :        coef_sigma_x_2_sqd &
         689             :        = one - sqrt( mixt_frac / ( one - mixt_frac ) ) &
         690             :                * Skx * sqrt( F_w * wp2 * xp2 ) / ( three * wpxp ) &
         691             :          + ( ( mixt_frac - two ) / ( one - mixt_frac ) ) &
         692           0 :            * wpxp**2 / ( three * F_w * wp2 * xp2 )
         693             :
         694             :        ! Mathematically, the value of coef_sigma_x_2_sqd cannot be less than 0.
         695             :        ! Numerically, this can happen when numerical round off error causes an
         696             :        ! epsilon-sized negative value.  When this happens, reset the value of
         697             :        ! coef_sigma_x_2_sqd to 0.
         698           0 :        coef_sigma_x_2_sqd = max( coef_sigma_x_2_sqd, zero )
         699             :
         700           0 :        sigma_x_2_sqd = coef_sigma_x_2_sqd * xp2
         701             :
         702             :     else ! | <w'x'> | = 0
         703             :
         704             :        ! When <w'x'> has a value of 0, the PDF becomes a single Gaussian.  This
         705             :        ! only works when Skx = 0.  However, when Skx /= 0, the value of min_F_x
         706             :        ! is greater than 0, preventing a problem where F_x = 0 but | Skx | > 0.
         707           0 :        mu_x_1 = xm
         708           0 :        mu_x_2 = xm
         709           0 :        sigma_x_1_sqd = xp2
         710           0 :        sigma_x_2_sqd = xp2
         711           0 :        coef_sigma_x_1_sqd = one
         712           0 :        coef_sigma_x_2_sqd = one
         713             :
         714             :     endif ! | <w'x'> | > 0
         715             :
         716             :
         717           0 :     return
         718             :
         719             :   end subroutine calculate_responder_params
         720             :
         721             :   !=============================================================================
         722           0 :   elemental function calculate_coef_wp4_implicit( mixt_frac, F_w, &
         723             :                                                   coef_sigma_w_1_sqd, &
         724             :                                                   coef_sigma_w_2_sqd ) &
         725             :   result( coef_wp4_implicit )
         726             :
         727             :     ! Description:
         728             :     ! The predictive equation for <w'^3> contains a turbulent advection term of
         729             :     ! the form:
         730             :     !
         731             :     ! - ( 1 / rho_ds ) * d ( rho_ds * <w'^4> ) / dz;
         732             :     !
         733             :     ! where z is height, rho_ds is the dry, base-state density, and <w'^4> is
         734             :     ! calculated by integrating over the PDF.  The equation for <w'^4> is:
         735             :     !
         736             :     ! <w'^4> = mixt_frac * ( 3 * sigma_w_1^4
         737             :     !                        + 6 * ( mu_w_1 - <w> )^2 * sigma_w_1^2
         738             :     !                        + ( mu_w_1 - <w> )^4 )
         739             :     !          + ( 1 - mixt_frac ) * ( 3 * sigma_w_2^4
         740             :     !                                  + 6 * ( mu_w_2 - <w> )^2 * sigma_w_2^2
         741             :     !                                  + ( mu_w_2 - <w> )^4 ).
         742             :     !
         743             :     ! The following substitutions are made into the above equation:
         744             :     !
         745             :     ! mu_w_1 - <w> = sqrt(F_w) * sqrt( ( 1 - mixt_frac ) / mixt_frac )
         746             :     !                * sqrt( <w'^2> );
         747             :     !
         748             :     ! mu_w_2 - <w> = - sqrt(F_w) * sqrt( mixt_frac / ( 1 - mixt_frac ) )
         749             :     !                  * sqrt( <w'^2> );
         750             :     !
         751             :     ! sigma_w_1 = sqrt( coef_sigma_w_1_sqd * <w'^2> ); and
         752             :     !
         753             :     ! sigma_w_2 = sqrt( coef_sigma_w_2_sqd * <w'^2> ).
         754             :     !
         755             :     ! The equations for coef_sigma_w_1_sqd and coef_sigma_w_2_sqd are:
         756             :     !
         757             :     ! coef_sigma_w_1_sqd = ( ( zeta_w + 1 ) * ( 1 - F_w ) )
         758             :     !                      / ( ( zeta_w + 2 ) * mixt_frac ); and
         759             :     !
         760             :     ! coef_sigma_w_2_sqd = ( 1 - F_w ) / ( ( zeta_w + 2 ) * ( 1 - mixt_frac ) ).
         761             :     !
         762             :     ! The equation for <w'4> becomes:
         763             :     !
         764             :     ! <w'^4> = ( 3 * mixt_frac * coef_sigma_w_1_sqd^2
         765             :     !            + 6 * F_w * ( 1 - mixt_frac ) * coef_sigma_w_1_sqd
         766             :     !            + F_w^2 * ( 1 - mixt_frac )^2 / mixt_frac
         767             :     !            + 3 * ( 1 - mixt_frac ) * coef_sigma_w_2_sqd^2
         768             :     !            + 6 * F_w * mixt_frac * coef_sigma_w_2_sqd
         769             :     !            + F_w^2 * mixt_frac^2 / ( 1 - mixt_frac ) ) * <w'^2>^2.
         770             :     !
         771             :     ! This equation is of the form:
         772             :     !
         773             :     ! <w'^4> = coef_wp4_implicit * <w'^2>^2;
         774             :     !
         775             :     ! where:
         776             :     !
         777             :     ! coef_wp4_implicit = 3 * mixt_frac * coef_sigma_w_1_sqd^2
         778             :     !                     + 6 * F_w * ( 1 - mixt_frac ) * coef_sigma_w_1_sqd
         779             :     !                     + F_w^2 * ( 1 - mixt_frac )^2 / mixt_frac
         780             :     !                     + 3 * ( 1 - mixt_frac ) * coef_sigma_w_2_sqd^2
         781             :     !                     + 6 * F_w * mixt_frac * coef_sigma_w_2_sqd
         782             :     !                     + F_w^2 * mixt_frac^2 / ( 1 - mixt_frac ).
         783             :     !
         784             :     ! While the <w'^4> term is found in the <w'^3> predictive equation and not
         785             :     ! the <w'^2> predictive equation, the <w'^3> and <w'^2> predictive equations
         786             :     ! are solved together.  This allows the term containing <w'^4> to be solved
         787             :     ! implicitly.
         788             :
         789             :     ! References:
         790             :     !-----------------------------------------------------------------------
         791             :
         792             :     use constants_clubb, only: &
         793             :         six,   & ! Variable(s)
         794             :         three, &
         795             :         one
         796             :
         797             :     use clubb_precision, only: &
         798             :         core_rknd    ! Procedure(s)
         799             :
         800             :     implicit none
         801             :
         802             :     ! Input Variables
         803             :     real ( kind = core_rknd ), intent(in) :: &
         804             :       mixt_frac,          & ! Mixture fraction                               [-]
         805             :       F_w,                & ! Parameter: spread of the PDF comp. means of w  [-]
         806             :       coef_sigma_w_1_sqd, & ! sigma_w_1^2 = coef_sigma_w_1_sqd * <w'^2>      [-]
         807             :       coef_sigma_w_2_sqd    ! sigma_w_2^2 = coef_sigma_w_2_sqd * <w'^2>      [-]
         808             :
         809             :     ! Return Variable
         810             :     real ( kind = core_rknd ) :: &
         811             :       coef_wp4_implicit    ! Coef.: <w'^4> = coef_wp4_implicit * <w'^2>^2    [-]
         812             :
         813             :
         814             :     ! Calculate coef_wp4_implicit.
         815             :     coef_wp4_implicit = three * mixt_frac * coef_sigma_w_1_sqd**2 &
         816             :                         + six * F_w * ( one - mixt_frac ) * coef_sigma_w_1_sqd &
         817             :                         + F_w**2 * ( one - mixt_frac )**2 / mixt_frac &
         818             :                         + three * ( one - mixt_frac ) * coef_sigma_w_2_sqd**2 &
         819             :                         + six * F_w * mixt_frac * coef_sigma_w_2_sqd &
         820           0 :                         + F_w**2 * mixt_frac**2 / ( one - mixt_frac )
         821             :
         822             :
         823             :     return
         824             :
         825             :   end function calculate_coef_wp4_implicit
         826             :
         827             :   !=============================================================================
         828           0 :   elemental function calc_coef_wp2xp_implicit( wp2, mixt_frac, F_w, &
         829             :                                                coef_sigma_w_1_sqd,  &
         830             :                                                coef_sigma_w_2_sqd   ) &
         831             :   result( coef_wp2xp_implicit )
         832             :
         833             :     ! Description:
         834             :     ! The predictive equation for <w'x'> contains a turbulent advection term of
         835             :     ! the form:
         836             :     !
         837             :     ! - ( 1 / rho_ds ) * d ( rho_ds * <w'^2 x'> ) / dz;
         838             :     !
         839             :     ! where z is height, rho_ds is the dry, base-state density, and <w'^2 x'> is
         840             :     ! calculated by integrating over the PDF.  The equation for <w'^2 x'> is:
         841             :     !
         842             :     ! <w'^2 x'>
         843             :     ! = mixt_frac * ( ( mu_x_1 - <x> ) * ( ( mu_w_1 - <w> )^2 + sigma_w_1^2 )
         844             :     !                 + 2 * corr_w_x_1 * sigma_w_1 * sigma_x_1
         845             :     !                   * ( mu_w_1 - <w> ) )
         846             :     !   + ( 1 - mixt_frac ) * ( ( mu_x_2 - <x> )
         847             :     !                           * ( ( mu_w_2 - <w> )^2 + sigma_w_2^2 )
         848             :     !                           + 2 * corr_w_x_2 * sigma_w_2 * sigma_x_2
         849             :     !                             * ( mu_w_2 - <w> ) ).
         850             :     !
         851             :     ! The following substitutions are made into the above equation:
         852             :     !
         853             :     ! mu_w_1 - <w> = sqrt(F_w) * sqrt( ( 1 - mixt_frac ) / mixt_frac )
         854             :     !                * sqrt( <w'^2> );
         855             :     !
         856             :     ! mu_w_2 - <w> = - sqrt(F_w) * sqrt( mixt_frac / ( 1 - mixt_frac ) )
         857             :     !                  * sqrt( <w'^2> );
         858             :     !
         859             :     ! mu_x_1 - <x> = sqrt( ( 1 - mixt_frac ) / mixt_frac )
         860             :     !                * <w'x'> / sqrt( F_w * <w'^2> );
         861             :     !
         862             :     ! mu_x_2 - <x> = - sqrt( mixt_frac / ( 1 - mixt_frac ) )
         863             :     !                  * <w'x'> / sqrt( F_w * <w'^2> );
         864             :     !
         865             :     ! sigma_w_1 = sqrt( coef_sigma_w_1_sqd * <w'^2> );
         866             :     !
         867             :     ! sigma_w_2 = sqrt( coef_sigma_w_2_sqd * <w'^2> );
         868             :     !
         869             :     ! sigma_x_1 = sqrt( coef_sigma_x_1_sqd * <x'^2> ); and
         870             :     !
         871             :     ! sigma_x_2 = sqrt( coef_sigma_x_2_sqd * <x'^2> ).
         872             :     !
         873             :     ! The equations for coef_sigma_w_1_sqd and coef_sigma_w_2_sqd are:
         874             :     !
         875             :     ! coef_sigma_w_1_sqd = ( ( zeta_w + 1 ) * ( 1 - F_w ) )
         876             :     !                      / ( ( zeta_w + 2 ) * mixt_frac ); and
         877             :     !
         878             :     ! coef_sigma_w_2_sqd = ( 1 - F_w ) / ( ( zeta_w + 2 ) * ( 1 - mixt_frac ) ).
         879             :     !
         880             :     ! The equations for coef_sigma_x_1_sqd and coef_sigma_x_2_sqd are:
         881             :     !
         882             :     ! coef_sigma_x_1_sqd
         883             :     ! = 1 + sqrt( ( 1 - mixt_frac ) / mixt_frac )
         884             :     !       * Skx * sqrt( F_w * <w'^2> * <x'^2> ) / ( 3 * <w'x'> )
         885             :     !   - ( ( 1 + mixt_frac ) / mixt_frac )
         886             :     !     * <w'x'>^2 / ( 3 * F_w * <w'^2> * <x'^2> ); and
         887             :     !
         888             :     ! coef_sigma_x_2_sqd
         889             :     ! = 1 - sqrt( mixt_frac / ( 1 - mixt_frac ) )
         890             :     !       * Skx * sqrt( F_w * <w'^2> * <x'^2> ) / ( 3 * <w'x'> )
         891             :     !   + ( ( mixt_frac - 2 ) / ( 1 - mixt_frac ) )
         892             :     !     * <w'x'>^2 / ( 3 * F_w * <w'^2> * <x'^2> ).
         893             :     !
         894             :     ! Additionally, corr_w_x_1 = corr_w_x_2 = 0.
         895             :     !
         896             :     ! The equation for <w'^2 x'> becomes:
         897             :     !
         898             :     ! <w'^2 x'> = sqrt( mixt_frac * ( 1 - mixt_frac ) )
         899             :     !             * ( F_w * ( ( 1 - mixt_frac ) / mixt_frac
         900             :     !                         - mixt_frac / ( 1 - mixt_frac ) )
         901             :     !                 + coef_sigma_w_1_sqd - coef_sigma_w_2_sqd )
         902             :     !             * sqrt( <w'^2> / F_w ) * <w'x'>.
         903             :     !
         904             :     ! This equation is of the form:
         905             :     !
         906             :     ! <w'^2 x'> = coef_wp2xp_implicit * <w'x'>;
         907             :     !
         908             :     ! where:
         909             :     !
         910             :     ! coef_wp2xp_implicit = sqrt( mixt_frac * ( 1 - mixt_frac ) )
         911             :     !                       * ( F_w * ( ( 1 - mixt_frac ) / mixt_frac
         912             :     !                                   - mixt_frac / ( 1 - mixt_frac ) )
         913             :     !                           + coef_sigma_w_1_sqd - coef_sigma_w_2_sqd )
         914             :     !                       * sqrt( <w'^2> / F_w ).
         915             :     !
         916             :     ! In the special case that F_w = 0, <w'x'> must have a value of 0, in which
         917             :     ! case mu_x_1 - <x> = mu_x_2 - <x> = 0, and <w'^2 x'> = 0.  The value of
         918             :     ! coef_wp2xp_implicit when F_w = 0 and <w'x'> = 0 can be calculated since
         919             :     ! Skw must also have a value of 0 when F_w = 0.  When F_w = 0 and Skw = 0,
         920             :     ! mixt_frac = ( zeta_w + 1 ) / ( zeta_w + 2 ).  When this happens,
         921             :     ! coef_sigma_w_1_sqd - coef_sigma_w_2_sqd = 0.  The equation for
         922             :     ! coef_wp2xp_implicit becomes:
         923             :     !
         924             :     ! coef_wp2xp_implicit = sqrt( mixt_frac * ( 1 - mixt_frac ) )
         925             :     !                       * ( ( 1 - mixt_frac ) / mixt_frac
         926             :     !                           - mixt_frac / ( 1 - mixt_frac ) )
         927             :     !                       * sqrt( F_w * <w'^2> );
         928             :     !
         929             :     ! and since F_w = 0, coef_wp2xp_implicit = 0.
         930             :
         931             :     ! References:
         932             :     !-----------------------------------------------------------------------
         933             :
         934             :     use constants_clubb, only: &
         935             :         one,  & ! Variable(s)
         936             :         zero
         937             :
         938             :     use clubb_precision, only: &
         939             :         core_rknd    ! Procedure(s)
         940             :
         941             :     implicit none
         942             :
         943             :     ! Input Variables
         944             :     real ( kind = core_rknd ), intent(in) :: &
         945             :       wp2,                & ! Variance of w (overall)                  [m^2/s^2]
         946             :       mixt_frac,          & ! Mixture fraction                               [-]
         947             :       F_w,                & ! Parameter: spread of the PDF comp. means of w  [-]
         948             :       coef_sigma_w_1_sqd, & ! sigma_w_1^2 = coef_sigma_w_1_sqd * <w'^2>      [-]
         949             :       coef_sigma_w_2_sqd    ! sigma_w_2^2 = coef_sigma_w_2_sqd * <w'^2>      [-]
         950             :
         951             :     ! Return Variable
         952             :     ! Coefficient: <w'^2 x'> = coef_wp2xp_implicit * <w'x'>
         953             :     real ( kind = core_rknd ) :: &
         954             :       coef_wp2xp_implicit    ! Coefficient that is multiplied by <w'x'>    [m/s]
         955             :
         956             :
         957             :     ! Calculate coef_wp2xp_implicit.
         958           0 :     if ( F_w > 0 ) then
         959             :
         960             :        coef_wp2xp_implicit &
         961             :        = sqrt( mixt_frac * ( one - mixt_frac ) ) &
         962             :          * ( F_w * ( ( one - mixt_frac ) / mixt_frac &
         963             :                      - mixt_frac / ( one - mixt_frac ) ) &
         964             :              + coef_sigma_w_1_sqd - coef_sigma_w_2_sqd ) &
         965           0 :          * sqrt( wp2 / F_w )
         966             :
         967             :     else ! F_w = 0
         968             :
         969             :        coef_wp2xp_implicit = zero
         970             :
         971             :     endif
         972             :
         973             :
         974             :     return
         975             :
         976             :   end function calc_coef_wp2xp_implicit
         977             :
         978             :   !=============================================================================
         979           0 :   elemental subroutine calc_coefs_wpxp2_semiimpl( wp2, wpxp,           & ! In
         980             :                                                   mixt_frac, F_w,      & ! In
         981             :                                                   coef_sigma_x_1_sqd,  & ! In
         982             :                                                   coef_sigma_x_2_sqd,  & ! In
         983             :                                                   coef_wpxp2_implicit, & ! Out
         984             :                                                   term_wpxp2_explicit )  ! Out
         985             :
         986             :     ! Description:
         987             :     ! The predictive equation for <x'^2> contains a turbulent advection term of
         988             :     ! the form:
         989             :     !
         990             :     ! - ( 1 / rho_ds ) * d ( rho_ds * <w'x'^2> ) / dz;
         991             :     !
         992             :     ! where z is height, rho_ds is the dry, base-state density, and <w'x'^2> is
         993             :     ! calculated by integrating over the PDF.  The equation for <w'x'^2> is:
         994             :     !
         995             :     ! <w'x'^2>
         996             :     ! = mixt_frac * ( ( mu_w_1 - <w> ) * ( ( mu_x_1 - <x> )^2 + sigma_x_1^2 )
         997             :     !                 + 2 * corr_w_x_1 * sigma_w_1 * sigma_x_1
         998             :     !                   * ( mu_x_1 - <x> ) )
         999             :     !   + ( 1 - mixt_frac ) * ( ( mu_w_2 - <w> )
        1000             :     !                           * ( ( mu_x_2 - <x> )^2 + sigma_x_2^2 )
        1001             :     !                           + 2 * corr_w_x_2 * sigma_w_2 * sigma_x_2
        1002             :     !                             * ( mu_x_2 - <x> ) ).
        1003             :     !
        1004             :     ! The following substitutions are made into the above equation:
        1005             :     !
        1006             :     ! mu_w_1 - <w> = sqrt(F_w) * sqrt( ( 1 - mixt_frac ) / mixt_frac )
        1007             :     !                * sqrt( <w'^2> );
        1008             :     !
        1009             :     ! mu_w_2 - <w> = - sqrt(F_w) * sqrt( mixt_frac / ( 1 - mixt_frac ) )
        1010             :     !                  * sqrt( <w'^2> );
        1011             :     !
        1012             :     ! mu_x_1 - <x> = sqrt( ( 1 - mixt_frac ) / mixt_frac )
        1013             :     !                * <w'x'> / sqrt( F_w * <w'^2> );
        1014             :     !
        1015             :     ! mu_x_2 - <x> = - sqrt( mixt_frac / ( 1 - mixt_frac ) )
        1016             :     !                  * <w'x'> / sqrt( F_w * <w'^2> );
        1017             :     !
        1018             :     ! sigma_w_1 = sqrt( coef_sigma_w_1_sqd * <w'^2> );
        1019             :     !
        1020             :     ! sigma_w_2 = sqrt( coef_sigma_w_2_sqd * <w'^2> );
        1021             :     !
        1022             :     ! sigma_x_1 = sqrt( coef_sigma_x_1_sqd * <x'^2> ); and
        1023             :     !
        1024             :     ! sigma_x_2 = sqrt( coef_sigma_x_2_sqd * <x'^2> ).
        1025             :     !
        1026             :     ! The equations for coef_sigma_w_1_sqd and coef_sigma_w_2_sqd are:
        1027             :     !
        1028             :     ! coef_sigma_w_1_sqd = ( ( zeta_w + 1 ) * ( 1 - F_w ) )
        1029             :     !                      / ( ( zeta_w + 2 ) * mixt_frac ); and
        1030             :     !
        1031             :     ! coef_sigma_w_2_sqd = ( 1 - F_w ) / ( ( zeta_w + 2 ) * ( 1 - mixt_frac ) ).
        1032             :     !
        1033             :     ! The equations for coef_sigma_x_1_sqd and coef_sigma_x_2_sqd are:
        1034             :     !
        1035             :     ! coef_sigma_x_1_sqd
        1036             :     ! = 1 + sqrt( ( 1 - mixt_frac ) / mixt_frac )
        1037             :     !       * Skx * sqrt( F_w * <w'^2> * <x'^2> ) / ( 3 * <w'x'> )
        1038             :     !   - ( ( 1 + mixt_frac ) / mixt_frac )
        1039             :     !     * <w'x'>^2 / ( 3 * F_w * <w'^2> * <x'^2> ); and
        1040             :     !
        1041             :     ! coef_sigma_x_2_sqd
        1042             :     ! = 1 - sqrt( mixt_frac / ( 1 - mixt_frac ) )
        1043             :     !       * Skx * sqrt( F_w * <w'^2> * <x'^2> ) / ( 3 * <w'x'> )
        1044             :     !   + ( ( mixt_frac - 2 ) / ( 1 - mixt_frac ) )
        1045             :     !     * <w'x'>^2 / ( 3 * F_w * <w'^2> * <x'^2> ).
        1046             :     !
        1047             :     ! Additionally, corr_w_x_1 = corr_w_x_2 = 0.
        1048             :     !
        1049             :     ! The equation for <w'x'^2> becomes:
        1050             :     !
        1051             :     ! <w'x'^2>
        1052             :     ! = sqrt( mixt_frac * ( 1 - mixt_frac ) ) * sqrt( F_w * <w'^2> )
        1053             :     !   * ( coef_sigma_x_1_sqd - coef_sigma_x_2_sqd ) * <x'^2>
        1054             :     !   + sqrt( mixt_frac * ( 1 - mixt_frac ) )
        1055             :     !     * <w'x'>^2 / sqrt( F_w * <w'^2> )
        1056             :     !     * ( ( 1 - mixt_frac ) / mixt_frac  - mixt_frac / ( 1 - mixt_frac ) ).
        1057             :     !
        1058             :     ! This equation is of the form:
        1059             :     !
        1060             :     ! <w'x'^2> = coef_wpxp2_implicit * <x'^2> + term_wpxp2_explicit;
        1061             :     !
        1062             :     ! where:
        1063             :     !
        1064             :     ! coef_wpxp2_implicit
        1065             :     ! = sqrt( mixt_frac * ( 1 - mixt_frac ) ) * sqrt( F_w * <w'^2> )
        1066             :     !   * ( coef_sigma_x_1_sqd - coef_sigma_x_2_sqd ); and
        1067             :     !
        1068             :     ! term_wpxp2_explicit
        1069             :     ! = sqrt( mixt_frac * ( 1 - mixt_frac ) ) * <w'x'>^2 / sqrt( F_w * <w'^2> )
        1070             :     !   * ( ( 1 - mixt_frac ) / mixt_frac  - mixt_frac / ( 1 - mixt_frac ) ).
        1071             :     !
        1072             :     ! In the special case that F_w = 0, mu_w_1 - <w> = mu_w_2 - <w> = 0, and
        1073             :     ! <w'x'^2> = 0.  The value of coef_wp2xp_implicit = 0 when F_w = 0.
        1074             :     ! Likewise, the value of <w'x'> = 0 when F_w = 0, which makes
        1075             :     ! term_wp2xp_explicit undefined in this scenario.  However, since both
        1076             :     ! <w'x'^2> = 0 and coef_wp2xp_implicit = 0, term_wp2xp_explicit = 0 in this
        1077             :     ! situation.
        1078             :
        1079             :     ! References:
        1080             :     !-----------------------------------------------------------------------
        1081             :
        1082             :     use constants_clubb, only: &
        1083             :         one,  & ! Variable(s)
        1084             :         zero
        1085             :
        1086             :     use clubb_precision, only: &
        1087             :         core_rknd    ! Procedure(s)
        1088             :
        1089             :     implicit none
        1090             :
        1091             :     ! Input Variables
        1092             :     real ( kind = core_rknd ), intent(in) :: &
        1093             :       wp2,                & ! Variance of w (overall)                  [m^2/s^2]
        1094             :       wpxp,               & ! Covariance of w and x           [m/s (units vary)]
        1095             :       mixt_frac,          & ! Mixture fraction                               [-]
        1096             :       F_w,                & ! Parameter: spread of the PDF comp. means of w  [-]
        1097             :       coef_sigma_x_1_sqd, & ! sigma_x_1^2 = coef_sigma_x_1_sqd * <x'^2>      [-]
        1098             :       coef_sigma_x_2_sqd    ! sigma_x_2^2 = coef_sigma_x_2_sqd * <x'^2>      [-]
        1099             :
        1100             :     ! Output Variables
        1101             :     ! Coefs.: <w'x'^2> = coef_wpxp2_implicit * <x'^2> + term_wpxp2_explicit
        1102             :     real ( kind = core_rknd ), intent(out) :: &
        1103             :       coef_wpxp2_implicit, & ! Coefficient that is multiplied by <x'^2>  [m/s]
        1104             :       term_wpxp2_explicit    ! Term that is on the RHS    [m/s (units vary)^2]
        1105             :
        1106             :
        1107             :     ! Calculate coef_wpxp2_implicit and term_wpxp2_explicit.
        1108           0 :     if ( F_w > 0 .and. wp2 > 0 ) then
        1109             :
        1110             :        coef_wpxp2_implicit &
        1111             :        = sqrt( mixt_frac * ( one - mixt_frac ) ) * sqrt( F_w * wp2 ) &
        1112           0 :          * ( coef_sigma_x_1_sqd - coef_sigma_x_2_sqd )
        1113             :
        1114             :        term_wpxp2_explicit &
        1115             :        = sqrt( mixt_frac * ( one - mixt_frac ) ) * wpxp**2 / sqrt( F_w * wp2 ) &
        1116           0 :          * ( ( one - mixt_frac ) / mixt_frac - mixt_frac / ( one - mixt_frac ) )
        1117             :
        1118             :     else ! F_w = 0 or wp2 = 0
        1119             :
        1120           0 :        coef_wpxp2_implicit = zero
        1121           0 :        term_wpxp2_explicit = zero
        1122             :
        1123             :     endif ! F_w > 0 and wp2 > 0
        1124             :
        1125             :
        1126           0 :     return
        1127             :
        1128             :   end subroutine calc_coefs_wpxp2_semiimpl
        1129             :
        1130             :   !=============================================================================
        1131           0 :   elemental subroutine calc_coefs_wpxpyp_semiimpl( wp2, wpxp, wpyp,      & ! In
        1132             :                                                    mixt_frac, F_w,       & ! In
        1133             :                                                    coef_sigma_x_1_sqd,   & ! In
        1134             :                                                    coef_sigma_x_2_sqd,   & ! In
        1135             :                                                    coef_sigma_y_1_sqd,   & ! In
        1136             :                                                    coef_sigma_y_2_sqd,   & ! In
        1137             :                                                    coef_wpxpyp_implicit, & ! Out
        1138             :                                                    term_wpxpyp_explicit  ) ! Out
        1139             :
        1140             :     ! Description:
        1141             :     ! The predictive equation for <w'x'y'> contains a turbulent advection term
        1142             :     ! of the form:
        1143             :     !
        1144             :     ! - ( 1 / rho_ds ) * d ( rho_ds * <w'x'y'> ) / dz;
        1145             :     !
        1146             :     ! where z is height, rho_ds is the dry, base-state density, and <w'x'y'> is
        1147             :     ! calculated by integrating over the PDF.  The equation for <w'x'y'> is:
        1148             :     !
        1149             :     ! <w'x'y'>
        1150             :     ! = mixt_frac
        1151             :     !   * ( ( mu_w_1 - <w> ) * ( mu_x_1 - <x> ) * ( mu_y_1 - <y> )
        1152             :     !       + corr_x_y_1 * sigma_x_1 * sigma_y_1 * ( mu_w_1 - <w> )
        1153             :     !       + corr_w_y_1 * sigma_w_1 * sigma_y_1 * ( mu_x_1 - <x> )
        1154             :     !       + corr_w_x_1 * sigma_w_1 * sigma_x_1 * ( mu_y_1 - <y> ) )
        1155             :     !   + ( 1 - mixt_frac )
        1156             :     !     * ( ( mu_w_2 - <w> ) * ( mu_x_2 - <x> ) * ( mu_y_2 - <y> )
        1157             :     !         + corr_x_y_2 * sigma_x_2 * sigma_y_2 * ( mu_w_2 - <w> )
        1158             :     !         + corr_w_y_2 * sigma_w_2 * sigma_y_2 * ( mu_x_2 - <x> )
        1159             :     !         + corr_w_x_2 * sigma_w_2 * sigma_x_2 * ( mu_y_2 - <y> ) ).
        1160             :     !
        1161             :     ! The following substitutions are made into the above equation:
        1162             :     !
        1163             :     ! mu_w_1 - <w> = sqrt(F_w) * sqrt( ( 1 - mixt_frac ) / mixt_frac )
        1164             :     !                * sqrt( <w'^2> );
        1165             :     !
        1166             :     ! mu_w_2 - <w> = - sqrt(F_w) * sqrt( mixt_frac / ( 1 - mixt_frac ) )
        1167             :     !                  * sqrt( <w'^2> );
        1168             :     !
        1169             :     ! mu_x_1 - <x> = sqrt( ( 1 - mixt_frac ) / mixt_frac )
        1170             :     !                * <w'x'> / sqrt( F_w * <w'^2> );
        1171             :     !
        1172             :     ! mu_x_2 - <x> = - sqrt( mixt_frac / ( 1 - mixt_frac ) )
        1173             :     !                  * <w'x'> / sqrt( F_w * <w'^2> );
        1174             :     !
        1175             :     ! mu_y_1 - <y> = sqrt( ( 1 - mixt_frac ) / mixt_frac )
        1176             :     !                * <w'y'> / sqrt( F_w * <w'^2> );
        1177             :     !
        1178             :     ! mu_y_2 - <y> = - sqrt( mixt_frac / ( 1 - mixt_frac ) )
        1179             :     !                  * <w'y'> / sqrt( F_w * <w'^2> );
        1180             :     !
        1181             :     ! sigma_w_1 = sqrt( coef_sigma_w_1_sqd * <w'^2> );
        1182             :     !
        1183             :     ! sigma_w_2 = sqrt( coef_sigma_w_2_sqd * <w'^2> );
        1184             :     !
        1185             :     ! sigma_x_1 = sqrt( coef_sigma_x_1_sqd * <x'^2> );
        1186             :     !
        1187             :     ! sigma_x_2 = sqrt( coef_sigma_x_2_sqd * <x'^2> );
        1188             :     !
        1189             :     ! sigma_y_1 = sqrt( coef_sigma_y_1_sqd * <y'^2> ); and
        1190             :     !
        1191             :     ! sigma_y_2 = sqrt( coef_sigma_y_2_sqd * <y'^2> ).
        1192             :     !
        1193             :     ! The equations for coef_sigma_w_1_sqd and coef_sigma_w_2_sqd are:
        1194             :     !
        1195             :     ! coef_sigma_w_1_sqd = ( ( zeta_w + 1 ) * ( 1 - F_w ) )
        1196             :     !                      / ( ( zeta_w + 2 ) * mixt_frac ); and
        1197             :     !
        1198             :     ! coef_sigma_w_2_sqd = ( 1 - F_w ) / ( ( zeta_w + 2 ) * ( 1 - mixt_frac ) ).
        1199             :     !
        1200             :     ! The equations for coef_sigma_x_1_sqd and coef_sigma_x_2_sqd are:
        1201             :     !
        1202             :     ! coef_sigma_x_1_sqd
        1203             :     ! = 1 + sqrt( ( 1 - mixt_frac ) / mixt_frac )
        1204             :     !       * Skx * sqrt( F_w * <w'^2> * <x'^2> ) / ( 3 * <w'x'> )
        1205             :     !   - ( ( 1 + mixt_frac ) / mixt_frac )
        1206             :     !     * <w'x'>^2 / ( 3 * F_w * <w'^2> * <x'^2> ); and
        1207             :     !
        1208             :     ! coef_sigma_x_2_sqd
        1209             :     ! = 1 - sqrt( mixt_frac / ( 1 - mixt_frac ) )
        1210             :     !       * Skx * sqrt( F_w * <w'^2> * <x'^2> ) / ( 3 * <w'x'> )
        1211             :     !   + ( ( mixt_frac - 2 ) / ( 1 - mixt_frac ) )
        1212             :     !     * <w'x'>^2 / ( 3 * F_w * <w'^2> * <x'^2> ).
        1213             :     !
        1214             :     ! The equations for coef_sigma_y_1_sqd and coef_sigma_y_2_sqd are:
        1215             :     !
        1216             :     ! coef_sigma_y_1_sqd
        1217             :     ! = 1 + sqrt( ( 1 - mixt_frac ) / mixt_frac )
        1218             :     !       * Sky * sqrt( F_w * <w'^2> * <y'^2> ) / ( 3 * <w'y'> )
        1219             :     !   - ( ( 1 + mixt_frac ) / mixt_frac )
        1220             :     !     * <w'y'>^2 / ( 3 * F_w * <w'^2> * <y'^2> ); and
        1221             :     !
        1222             :     ! coef_sigma_y_2_sqd
        1223             :     ! = 1 - sqrt( mixt_frac / ( 1 - mixt_frac ) )
        1224             :     !       * Sky * sqrt( F_w * <w'^2> * <y'^2> ) / ( 3 * <w'y'> )
        1225             :     !   + ( ( mixt_frac - 2 ) / ( 1 - mixt_frac ) )
        1226             :     !     * <w'y'>^2 / ( 3 * F_w * <w'^2> * <y'^2> ).
        1227             :     !
        1228             :     ! Additionally, corr_w_x_1 = corr_w_x_2 = corr_w_y_1 = corr_w_y_2 = 0; and:
        1229             :     !
        1230             :     ! corr_x_y_1 = corr_x_y_2
        1231             :     ! = ( <x'y'> - mixt_frac * ( mu_x_1 - <x> ) * ( mu_y_1 - <y> )
        1232             :     !            - ( 1 - mixt_frac ) * ( mu_x_2 - <x> ) * ( mu_y_2 - <y> ) )
        1233             :     !   / ( mixt_frac * sigma_x_1 * sigma_y_1
        1234             :     !       + ( 1 - mixt_frac ) * sigma_x_2 * sigma_y_2 );
        1235             :     !
        1236             :     ! where -1 <= corr_x_y_1 = corr_x_y_2 <= 1.  This equation can be rewritten
        1237             :     ! as:
        1238             :     !
        1239             :     ! corr_x_y_1 = corr_x_y_2
        1240             :     ! = ( <x'y'>
        1241             :     !     - sqrt( F_x ) * sqrt( F_y ) * sgn( <w'x'> ) * sgn( <w'y'> )
        1242             :     !       * sqrt( <x'^2> ) * sqrt( <y'^2 > ) )
        1243             :     !   / ( ( mixt_frac * sqrt( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd )
        1244             :     !         + ( 1 - mixt_frac )
        1245             :     !           * sqrt( coef_sigma_x_2_sqd * coef_sigma_y_2_sqd ) )
        1246             :     !       * sqrt( <x'^2> ) * sqrt( <y'^2> ) ).
        1247             :     !
        1248             :     ! The equation for <w'x'y'> becomes:
        1249             :     !
        1250             :     ! <w'x'y'>
        1251             :     ! = sqrt( mixt_frac * ( 1 - mixt_frac ) ) * sqrt( F_w * <w'^2> )
        1252             :     !   * ( sqrt( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd )
        1253             :     !             - sqrt( coef_sigma_x_2_sqd * coef_sigma_y_2_sqd ) )
        1254             :     !     / ( mixt_frac * sqrt( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd )
        1255             :     !         + ( 1 - mixt_frac )
        1256             :     !           * sqrt( coef_sigma_x_2_sqd * coef_sigma_y_2_sqd ) )
        1257             :     !   * <x'y'>
        1258             :     !   + sqrt( mixt_frac * ( 1 - mixt_frac ) )
        1259             :     !     * <w'x'> * <w'y'> / sqrt( F_w * <w'^2> )
        1260             :     !     * ( ( 1 - mixt_frac ) / mixt_frac - mixt_frac / ( 1 - mixt_frac )
        1261             :     !         - ( sqrt( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd )
        1262             :     !             - sqrt( coef_sigma_x_2_sqd * coef_sigma_y_2_sqd ) )
        1263             :     !           / ( mixt_frac * sqrt( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd )
        1264             :     !               + ( 1 - mixt_frac )
        1265             :     !                 * sqrt( coef_sigma_x_2_sqd * coef_sigma_y_2_sqd ) ) )
        1266             :     !
        1267             :     ! This equation is of the form:
        1268             :     !
        1269             :     ! <w'x'y'> = coef_wpxpyp_implicit * <x'y'> + term_wpxpyp_explicit;
        1270             :     !
        1271             :     ! where:
        1272             :     !
        1273             :     ! coef_wpxpyp_implicit
        1274             :     ! = sqrt( mixt_frac * ( 1 - mixt_frac ) ) * sqrt( F_w * <w'^2> )
        1275             :     !   * ( sqrt( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd )
        1276             :     !             - sqrt( coef_sigma_x_2_sqd * coef_sigma_y_2_sqd ) )
        1277             :     !     / ( mixt_frac * sqrt( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd )
        1278             :     !         + ( 1 - mixt_frac )
        1279             :     !           * sqrt( coef_sigma_x_2_sqd * coef_sigma_y_2_sqd ) ); and
        1280             :     !
        1281             :     ! term_wpxpyp_explicit
        1282             :     ! = sqrt( mixt_frac * ( 1 - mixt_frac ) )
        1283             :     !   * <w'x'> * <w'y'> / sqrt( F_w * <w'^2> )
        1284             :     !   * ( ( 1 - mixt_frac ) / mixt_frac - mixt_frac / ( 1 - mixt_frac )
        1285             :     !       - ( sqrt( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd )
        1286             :     !           - sqrt( coef_sigma_x_2_sqd * coef_sigma_y_2_sqd ) )
        1287             :     !         / ( mixt_frac * sqrt( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd )
        1288             :     !             + ( 1 - mixt_frac )
        1289             :     !               * sqrt( coef_sigma_x_2_sqd * coef_sigma_y_2_sqd ) ) ).
        1290             :     !
        1291             :     ! There are also special cases for the above equations.  In the scenario
        1292             :     ! that sigma_x_1 * sigma_y_1 = sigma_x_2 * sigma_y_2 = 0, and equations for
        1293             :     ! coef_wpxpyp_implicit and term_wpxpyp_explicit become:
        1294             :     !
        1295             :     ! coef_wpxpyp_implicit
        1296             :     ! = sqrt( mixt_frac * ( 1 - mixt_frac ) ) * sqrt( F_w * <w'^2> )
        1297             :     !   * ( ( 1 - mixt_frac ) / mixt_frac - mixt_frac / ( 1 - mixt_frac ) ); and
        1298             :     !
        1299             :     ! term_wpxpyp_explicit = 0.
        1300             :     !
        1301             :     ! In the scenario where F_w = 0, mu_w_1 - <w> = mu_w_2 - <w> = 0, and
        1302             :     ! <w'x'> = <w'y'> = 0, which means mu_x_1 - <x> = mu_x_2 - <x> = 0, as well
        1303             :     ! as mu_y_1 - <y> = mu_y_2 - <y> = 0, and as a result, <w'x'y'> = 0.  When
        1304             :     ! F_w = 0, coef_wpxpyp_implicit = 0.  Since <w'x'y'> = 0 and
        1305             :     ! coef_wpxpyp_implicit = 0 when F_w = 0, term_wpxpyp_explicit = 0.
        1306             :
        1307             :     ! References:
        1308             :     !-----------------------------------------------------------------------
        1309             :
        1310             :     use constants_clubb, only: &
        1311             :         one,  & ! Variable(s)
        1312             :         zero
        1313             :
        1314             :     use clubb_precision, only: &
        1315             :         core_rknd    ! Procedure(s)
        1316             :
        1317             :     implicit none
        1318             :
        1319             :     ! Input Variables
        1320             :     real ( kind = core_rknd ), intent(in) :: &
        1321             :       wp2,                & ! Variance of w (overall)                  [m^2/s^2]
        1322             :       wpxp,               & ! Covariance of w and x (overall)     [m/s(x units)]
        1323             :       wpyp,               & ! Covariance of w and y (overall)     [m/s(y units)]
        1324             :       mixt_frac,          & ! Mixture fraction                               [-]
        1325             :       F_w,                & ! Parameter: spread of the PDF comp. means of w  [-]
        1326             :       coef_sigma_x_1_sqd, & ! sigma_x_1^2 = coef_sigma_x_1_sqd * <x'^2>      [-]
        1327             :       coef_sigma_x_2_sqd, & ! sigma_x_2^2 = coef_sigma_x_2_sqd * <x'^2>      [-]
        1328             :       coef_sigma_y_1_sqd, & ! sigma_y_1^2 = coef_sigma_y_1_sqd * <y'^2>      [-]
        1329             :       coef_sigma_y_2_sqd    ! sigma_y_2^2 = coef_sigma_y_2_sqd * <y'^2>      [-]
        1330             :
        1331             :     ! Output Variables
        1332             :     ! Coefs.: <w'x'y'> = coef_wpxpyp_implicit * <x'y'> + term_wpxpyp_explicit
        1333             :     real ( kind = core_rknd ), intent(out) :: &
        1334             :       coef_wpxpyp_implicit, & ! Coefficient that is multiplied by <x'y'>   [m/s]
        1335             :       term_wpxpyp_explicit    ! Term that is on the RHS  [m/s(x units)(y units)]
        1336             :
        1337             :     ! Local Variables
        1338             :     real ( kind = core_rknd ) :: &
        1339             :       coefs_factor_xy    ! Factor involving coef_sigma_... x and y coefs     [-]
        1340             :
        1341             :
        1342             :     ! Calculate coef_wpxpyp_implicit and term_wpxpyp_explicit.
        1343             :     if ( ( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd > zero &
        1344             :            .or. coef_sigma_x_2_sqd * coef_sigma_y_2_sqd > zero ) &
        1345           0 :           .and. F_w > zero .and. wp2 > zero ) then
        1346             :
        1347             :        ! coefs_factor_xy
        1348             :        ! = ( sqrt( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd )
        1349             :        !     - sqrt( coef_sigma_x_2_sqd * coef_sigma_y_2_sqd ) )
        1350             :        !   / ( mixt_frac * sqrt( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd )
        1351             :        !       + ( 1 - mixt_frac )
        1352             :        !         * sqrt( coef_sigma_x_2_sqd * coef_sigma_y_2_sqd ) )
        1353             :        coefs_factor_xy &
        1354             :        = ( sqrt( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd ) &
        1355             :            - sqrt( coef_sigma_x_2_sqd * coef_sigma_y_2_sqd ) ) &
        1356             :          / ( mixt_frac * sqrt( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd ) &
        1357             :              + ( one - mixt_frac ) &
        1358           0 :                * sqrt( coef_sigma_x_2_sqd * coef_sigma_y_2_sqd ) )
        1359             :
        1360             :        coef_wpxpyp_implicit &
        1361             :        = sqrt( mixt_frac * ( 1 - mixt_frac ) ) &
        1362           0 :          * sqrt( F_w * wp2 ) * coefs_factor_xy
        1363             :
        1364             :        term_wpxpyp_explicit &
        1365             :        = sqrt( mixt_frac * ( one - mixt_frac ) ) &
        1366             :          * wpxp * wpyp / sqrt( F_w * wp2 ) &
        1367             :          * ( ( one - mixt_frac ) / mixt_frac &
        1368             :              - mixt_frac / ( one - mixt_frac ) &
        1369           0 :              - coefs_factor_xy )
        1370             :
        1371             :     else ! ( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd = 0
        1372             :          !   and coef_sigma_x_2_sqd * coef_sigma_y_2_sqd = 0 )
        1373             :          ! or F_w = 0 or wp2 = 0
        1374             :
        1375           0 :        if ( F_w > 0 ) then
        1376             :
        1377             :           coef_wpxpyp_implicit &
        1378             :           = sqrt( mixt_frac * ( one - mixt_frac ) ) * sqrt( F_w * wp2 ) &
        1379             :             * ( ( one - mixt_frac ) / mixt_frac &
        1380           0 :                 - mixt_frac / ( one - mixt_frac ) )
        1381             :
        1382             :        else ! F_w = 0
        1383             :
        1384           0 :           coef_wpxpyp_implicit = zero
        1385             :
        1386             :        endif
        1387             :
        1388           0 :        term_wpxpyp_explicit = zero
        1389             :
        1390             :     endif
        1391             :
        1392             :
        1393           0 :     return
        1394             :
        1395             :   end subroutine calc_coefs_wpxpyp_semiimpl
        1396             :
        1397             :   !=============================================================================
        1398             :
        1399             : end module new_hybrid_pdf

Generated by: LCOV version 1.14