LCOV - code coverage report
Current view: | top level - physics/clubb/src/CLUBB_core - new_pdf.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 0 | 122 | 0.0 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 0 | 9 | 0.0 %

              Line data    Source code
           1             : ! $Id$
           2             : !===============================================================================
           3             : module new_pdf
           4             :
           5             :   ! Description:
           6             :   ! The portion of CLUBB's multivariate, two-component PDF that is the
           7             :   ! trivariate, two-component normal PDF of vertical velocity (w), total water
           8             :   ! mixing ratio (rt), and liquid water potential temperature (thl).
           9             :
          10             :   ! References:
          11             :   ! Griffin and Larson (2018)
          12             :   !-------------------------------------------------------------------------
          13             :
          14             :   implicit none
          15             :
          16             :   public :: calc_mixture_fraction,      & ! Procedure(s)
          17             :             calc_setter_var_params,     &
          18             :             calc_responder_params,      &
          19             :             calc_limits_F_x_responder,  &
          20             :             calc_coef_wp4_implicit,     &
          21             :             calc_coef_wpxp2_implicit,   &
          22             :             calc_coefs_wp2xp_semiimpl,  &
          23             :             calc_coefs_wpxpyp_semiimpl
          24             :
          25             :   private :: sort_roots    ! Procedure(s)
          26             :
          27             :   private
          28             :
          29             :   contains
          30             :
          31             :   !=============================================================================
          32             :   !
          33             :   ! DESCRIPTION OF THE METHOD FOR THE VARIABLE THAT SETS THE MIXTURE FRACTION
          34             :   ! =========================================================================
          35             :   !
          36             :   ! Many times, w has been used as the variable that sets the mixture fraction
          37             :   ! for the PDF.  There are five PDF parameters that need to be calculated,
          38             :   ! which are mu_w_1 (the mean of w is the 1st PDF component), mu_w_2 (the mean
          39             :   ! of w in the 2nd PDF component), sigma_w_1 (the standard deviation of w in
          40             :   ! the 1st PDF component), sigma_w_2 (the standard deviation of w in the 2nd
          41             :   ! PDF component), and mixt_frac (the mixture fraction, which is the weight of
          42             :   ! the 1st PDF component).  In order to solve for these five parameters, five
          43             :   ! equations are needed.  These five equations are the equations for <w>,
          44             :   ! <w'^2>, and <w'^3> as found by integrating over the PDF.  Additionally, two
          45             :   ! more equations, which involve tunable parameters F_w and zeta_w, and which
          46             :   ! are used to help control the spread of the PDF component means and the size
          47             :   ! of the PDF component standard deviations compared to each other,
          48             :   ! respectively, are used in this equation set.  The five equations are:
          49             :   !
          50             :   ! <w> = mixt_frac * mu_w_1 + ( 1 - mixt_frac ) * mu_w_2;
          51             :   !
          52             :   ! <w'^2> = mixt_frac * ( ( mu_w_1 - <w> )^2 + sigma_w_1^2 )
          53             :   !          + ( 1 - mixt_frac ) * ( ( mu_w_2 - <w> )^2 + sigma_w_2^2 );
          54             :   !
          55             :   ! <w'^3> = mixt_frac * ( mu_w_1 - <w> )
          56             :   !                    * ( ( mu_w_1 - <w> )^2 + 3 * sigma_w_1^2 )
          57             :   !          + ( 1 - mixt_frac ) * ( mu_w_2 - <w> )
          58             :   !                              * ( ( mu_w_2 - <w> )^2 + 3 * sigma_w_2^2 );
          59             :   !
          60             :   ! mu_w_1 - <w> = sqrt(F_w) * ( sqrt( 1 - mixt_frac ) / sqrt( mixt_frac ) )
          61             :   !                * sqrt( <w'^2> );
          62             :   !
          63             :   ! where 0 <= F_w <= 1; and
          64             :   !
          65             :   ! 1 + zeta_w = ( mixt_frac * sigma_w_1^2 )
          66             :   !              / ( ( 1 - mixt_frac ) * sigma_w_2^2 );
          67             :   !
          68             :   ! where zeta_w > -1.
          69             :   !
          70             :   ! Following convention for w, mu_w_1 is defined to be greater than or equal to
          71             :   ! mu_w_2 (and is also greater than or equal to <w>, while mu_w_2 is less than
          72             :   ! or equal to <w>).  This is relationship is found in the mu_w_1 - <w>
          73             :   ! equation above.
          74             :   !
          75             :   ! The resulting equations for the five PDF parameters are:
          76             :   !
          77             :   ! mixt_frac
          78             :   ! = ( 4 * F_w^3
          79             :   !     + 18 * F_w * ( zeta_w + 1 ) * ( 1 - F_w ) / ( zeta_w + 2 )
          80             :   !     + 6 * F_w^2 * ( 1 - F_w ) / ( zeta_w + 2 )
          81             :   !     + Skw^2
          82             :   !     - Skw * sqrt( 4 * F_w^3
          83             :   !                   + 12 * F_w^2 * ( 1 - F_w )
          84             :   !                   + 36 * F_w * ( zeta_w + 1 ) * ( 1 - F_w )^2
          85             :   !                     / ( zeta_w + 2 )^2
          86             :   !                   + Skw^2 ) )
          87             :   !   / ( 2 * F_w * ( F_w - 3 )^2 + 2 * Skw^2 );
          88             :   !
          89             :   ! mu_w_1 = <w> + sqrt( F_w * ( ( 1 - mixt_frac ) / mixt_frac ) * <w'^2> );
          90             :   !
          91             :   ! mu_w_2 = <w> - ( mixt_frac / ( 1 - mixt_frac ) ) * ( mu_w_1 - <w> );
          92             :   !
          93             :   ! sigma_w_1
          94             :   ! = sqrt( ( ( zeta_w + 1 ) * ( 1 - F_w ) )
          95             :   !         / ( ( zeta_w + 2 ) * mixt_frac ) * <w'^2> ); and
          96             :   !
          97             :   ! sigma_w_2
          98             :   ! = sqrt( ( mixt_frac * sigma_w_1^2 )
          99             :   !         / ( ( 1 - mixt_frac ) * ( 1 + zeta_w ) ) );
         100             :   !
         101             :   ! where Skw is the skewness of w, and Skw = <w'^3> / <w'^2>^(3/2).
         102             :   !
         103             :   ! This method works for all values of F_w (where 0 <= F_w <= 1) and zeta_w
         104             :   ! (where zeta_w > -1).
         105             :   !
         106             :   !
         107             :   ! Generalized equations for any variable, x, that sets the mixture fraction:
         108             :   !
         109             :   ! A slight alteration is made to the above equations in order to have any
         110             :   ! variable, x, set the mixture fraction.  The same five PDF parameters need to
         111             :   ! be calculated for the setting variable, which are mu_x_1 (the mean of x in
         112             :   ! the 1st PDF component), mu_x_2 (the mean of x in the 2nd PDF component),
         113             :   ! sigma_x_1 (the standard deviation of x in the 1st PDF component), sigma_x_2
         114             :   ! (the standard deviation of x in the 2nd PDF component), and mixt_frac (the
         115             :   ! mixture fraction).  Again, five equations are needed, and they are the
         116             :   ! equations for <x>, <x'^2>, and <x'^3> as found by integrating over the PDF,
         117             :   ! as well as the equations that involve tunable parameters F_x and zeta_x.
         118             :   ! However, the equation for F_x is multiplied by a new variable,
         119             :   ! sgn( <w'x'> ), where <w'x'> is the covariance of w and x, and sgn( <w'x'> )
         120             :   ! is given by:
         121             :   !
         122             :   ! sgn( <w'x'> ) = |  1, when <w'x'> >= 0;
         123             :   !                 | -1, when <w'x'> < 0.
         124             :   !
         125             :   ! The five equations are:
         126             :   !
         127             :   ! <x> = mixt_frac * mu_x_1 + ( 1 - mixt_frac ) * mu_x_2;
         128             :   !
         129             :   ! <x'^2> = mixt_frac * ( ( mu_x_1 - <x> )^2 + sigma_x_1^2 )
         130             :   !          + ( 1 - mixt_frac ) * ( ( mu_x_2 - <x> )^2 + sigma_x_2^2 );
         131             :   !
         132             :   ! <x'^3> = mixt_frac * ( mu_x_1 - <x> )
         133             :   !                    * ( ( mu_x_1 - <x> )^2 + 3 * sigma_x_1^2 )
         134             :   !          + ( 1 - mixt_frac ) * ( mu_x_2 - <x> )
         135             :   !                              * ( ( mu_x_2 - <x> )^2 + 3 * sigma_x_2^2 );
         136             :   !
         137             :   ! mu_x_1 - <x> = sqrt(F_x) * ( sqrt( 1 - mixt_frac ) / sqrt( mixt_frac ) )
         138             :   !                * sqrt( <x'^2> ) * sgn( <w'x'> );
         139             :   !
         140             :   ! where 0 <= F_x <= 1; and
         141             :   !
         142             :   ! 1 + zeta_x = ( mixt_frac * sigma_x_1^2 )
         143             :   !              / ( ( 1 - mixt_frac ) * sigma_x_2^2 );
         144             :   !
         145             :   ! where zeta_x > -1.
         146             :   !
         147             :   ! The only equations that are altered are the equation for mu_x_1 and the
         148             :   ! equation for mixt_frac, which now both contain a sgn( <w'x'> ).  The mu_x_2
         149             :   ! equation is not altered, but the sign of mu_x_2 - <x> will be the opposite
         150             :   ! of the sign of mu_x_1 - <x>.  The resulting equations for the five PDF
         151             :   ! parameters are:
         152             :   !
         153             :   ! mixt_frac
         154             :   ! = ( 4 * F_x^3
         155             :   !     + 18 * F_x * ( zeta_x + 1 ) * ( 1 - F_x ) / ( zeta_x + 2 )
         156             :   !     + 6 * F_x^2 * ( 1 - F_x ) / ( zeta_x + 2 )
         157             :   !     + Skx^2
         158             :   !     - Skx * sgn( <w'x'> )
         159             :   !           * sqrt( 4 * F_x^3
         160             :   !                   + 12 * F_x^2 * ( 1 - F_x )
         161             :   !                   + 36 * F_x * ( zeta_x + 1 ) * ( 1 - F_x )^2
         162             :   !                     / ( zeta_x + 2 )^2
         163             :   !                   + Skx^2 ) )
         164             :   !   / ( 2 * F_x * ( F_x - 3 )^2 + 2 * Skx^2 );
         165             :   !
         166             :   ! mu_x_1 = <x> + sqrt( F_x * ( ( 1 - mixt_frac ) / mixt_frac ) * <x'^2> )
         167             :   !                * sgn( <w'x'> );
         168             :   !
         169             :   ! mu_x_2 = <x> - ( mixt_frac / ( 1 - mixt_frac ) ) * ( mu_x_1 - <x> );
         170             :   !
         171             :   ! sigma_x_1
         172             :   ! = sqrt( ( ( zeta_x + 1 ) * ( 1 - F_x ) )
         173             :   !         / ( ( zeta_x + 2 ) * mixt_frac ) * <x'^2> ); and
         174             :   !
         175             :   ! sigma_x_2
         176             :   ! = sqrt( ( mixt_frac * sigma_x_1^2 )
         177             :   !         / ( ( 1 - mixt_frac ) * ( 1 + zeta_x ) ) );
         178             :   !
         179             :   ! where Skx is the skewness of x, and Skx = <x'^3> / <x'^2>^(3/2).
         180             :   !
         181             :   ! This method works for all values of F_x (where 0 <= F_x <= 1) and zeta_x
         182             :   ! (where zeta_x > -1).
         183             :   !
         184             :   ! When the generalized form is solved for w (x = w), sgn( <w'^2> ) = 1, and
         185             :   ! the equations are unaltered from the equations listed above for w.
         186             :   !
         187             :   !
         188             :   ! Special case:
         189             :   !
         190             :   ! When Skx = 0 and F_x = 0, the equation for mixt_frac is undefined.  The
         191             :   ! equation for mixture fraction in this scenario can be derived by using the
         192             :   ! above equation for mixture fraction and then setting Skx = 0.  The resulting
         193             :   ! equation becomes:
         194             :   !
         195             :   ! mixt_frac
         196             :   ! = ( 4 * F_x^3
         197             :   !     + 18 * F_x * ( zeta_x + 1 ) * ( 1 - F_x ) / ( zeta_x + 2 )
         198             :   !     + 6 * F_x^2 * ( 1 - F_x ) / ( zeta_x + 2 ) )
         199             :   !   / ( 2 * F_x * ( F_x - 3 )^2 ).
         200             :   !
         201             :   ! All of the terms in the numerator and denominator contain a F_x, so this
         202             :   ! equation can be rewritten as:
         203             :   !
         204             :   ! mixt_frac
         205             :   ! = ( 4 * F_x^2
         206             :   !     + 18 * ( zeta_x + 1 ) * ( 1 - F_x ) / ( zeta_x + 2 )
         207             :   !     + 6 * F_x * ( 1 - F_x ) / ( zeta_x + 2 ) )
         208             :   !   / ( 2 * ( F_x - 3 )^2 ).
         209             :   !
         210             :   ! Now setting F_x = 0, the equation becomes:
         211             :   !
         212             :   ! mixt_frac = ( 18 * ( zeta_x + 1 ) / ( zeta_x + 2 ) ) / 18;
         213             :   !
         214             :   ! which can be rewritten as:
         215             :   !
         216             :   ! mixt_frac = ( zeta_x + 1 ) / ( zeta_x + 2 ).
         217             :   !
         218             :   ! When F_x = 0, Skx must have a value of 0 in order for the PDF to function
         219             :   ! correctly.  When F_x = 0, mu_x_1 = mu_x_2.  When the two PDF component means
         220             :   ! are equal to each other (and to the overall mean, <x>), the only value of
         221             :   ! Skx that can be represented is a value of 0.  In the equation for mixture
         222             :   ! fraction, when F_x is set to 0, but | Skx | > 0, mixt_frac will either have
         223             :   ! a value of 0 or 1, depending on whether Skx is positive or negative,
         224             :   ! respectively.
         225             :   !
         226             :   ! The value of F_x should be set as a function of Skx.  The value F_x should
         227             :   ! go toward 0 as | Skx | (or Skx^2) goes toward 0.  The value of F_x should
         228             :   ! go toward 1 as | Skx | (or Skx^2) goes to infinity.
         229             :   !
         230             :   !
         231             :   ! Tunable parameters:
         232             :   !
         233             :   ! 1) F_x:  This parameter controls the spread of the PDF component means.  The
         234             :   !          range of this parameter is 0 <= F_x <= 1.  When F_x = 0, the two
         235             :   !          PDF component means (mu_x_1 and mu_x_2) are equal to each other
         236             :   !          (and Skx must equal 0).  All of the variance of x is accounted for
         237             :   !          by the PDF component standard deviations (sigma_x_1 and sigma_x_2).
         238             :   !          When F_x = 1, mu_x_1 and mu_x_2 are spread as far apart as they can
         239             :   !          be.  Both PDF component standard deviations (sigma_x_1 and
         240             :   !          sigma_x_2) are equal to 0, and all of the variance of x is
         241             :   !          accounted for by the spread of the PDF component means.
         242             :   !
         243             :   !          When sigma_x_1 = sigma_x_2 = 0, the equation for <x'^2> becomes:
         244             :   !
         245             :   !          <x'^2> = mixt_frac * ( mu_x_1 - <x> )^2
         246             :   !                   + ( 1 - mixt_frac ) * ( mu_x_2 - <x> )^2.
         247             :   !
         248             :   !          Substituting the equation for <x> into the above equation for
         249             :   !          mu_x_2 - <x>, the above equation becomes:
         250             :   !
         251             :   !          <x'^2> = ( mixt_frac / ( 1 - mixt_frac ) ) * ( mu_x_1 - <x> )^2;
         252             :   !
         253             :   !          which can be rewritten as:
         254             :   !
         255             :   !          ( mu_x_1 - <x> )^2 = ( ( 1 - mixt_frac ) / mixt_frac ) * <x'^2>.
         256             :   !
         257             :   !          Taking the square root of the above equation:
         258             :   !
         259             :   !          mu_x_1 - <x> = +/- ( sqrt( 1 - mixt_frac ) / sqrt(mixt_frac) )
         260             :   !                             * sqrt( <x'^2> ).
         261             :   !
         262             :   !          This equation can be compared to the equation for mu_x_1 - <x> in
         263             :   !          the set of 5 equations, which is:
         264             :   !
         265             :   !          mu_x_1 - <x>
         266             :   !          = sqrt(F_x) * ( sqrt( 1 - mixt_frac ) / sqrt( mixt_frac ) )
         267             :   !            * sqrt( <x'^2> ) * sgn( <w'x'> ).
         268             :   !
         269             :   !          The above equations give another example of the meaning of F_x.
         270             :   !          The value of sqrt(F_x) is ratio of mu_x_1 - <x> to its maximum
         271             :   !          value (or minimum value, depending on sgn( <w'x'> )), which is:
         272             :   !
         273             :   !          sqrt( ( ( 1 - mixt_frac ) / mixt_frac ) * <x'^2> ) * sgn( <w'x'> ).
         274             :   !
         275             :   !
         276             :   ! 2) zeta_x:  This parameter controls the size of the PDF component standard
         277             :   !             deviations compared to each other.  The equation for zeta_x is:
         278             :   !
         279             :   !             1 + zeta_x = ( mixt_frac * sigma_x_1^2 )
         280             :   !                          / ( ( 1 - mixt_frac ) * sigma_x_2^2 ).
         281             :   !
         282             :   !             When zeta_x > 0, mixt_frac * sigma_x_1^2 increases at the
         283             :   !             expense of ( 1 - mixt_frac ) * sigma_x_2^2, which decreases in
         284             :   !             this variance-preserving equation set.  When zeta_x = 0, then
         285             :   !             mixt_frac * sigma_x_1^2 = ( 1 - mixt_frac ) * sigma_x_2^2.
         286             :   !             When -1 < zeta_x < 0, ( 1 - mixt_frac ) * sigma_x_2^2 increases
         287             :   !             at the expense of mixt_frac * sigma_x_1^2, which decreases.  As
         288             :   !             a result, greater values of zeta_x cause the 1st PDF component
         289             :   !             to become broader while the 2nd PDF component becomes narrower,
         290             :   !             and smaller values of zeta_x cause the 1st PDF component to
         291             :   !             become narrower while the 2nd PDF component becomes broader.
         292             :   !
         293             :   !             Symmetry
         294             :   !
         295             :   !             When zeta_x = 0, the PDF is always symmetric.  In other words,
         296             :   !             the PDF at any positive value of Skx (for example, Skx = 2.5)
         297             :   !             will look like a mirror-image (reflection across the y-axis)
         298             :   !             of the PDF at a negative value of Skx of the same magnitude (in
         299             :   !             this example, Skx = -2.5).  However, when zeta_x /= 0, the PDF
         300             :   !             loses this quality and is not symmetric.
         301             :   !
         302             :   !             When symmetry is desired at values of zeta_x besides zeta_x = 0,
         303             :   !             the solution is to turn zeta_x into a function of Skx.  A basic
         304             :   !             example of a zeta_x skewness equation that produces a symmetric
         305             :   !             PDF for values of zeta_x other than 0 is:
         306             :   !
         307             :   !             zeta_x = | zeta_x_in,                      when Skx >= 0;
         308             :   !                      | ( 1 / ( 1 + zeta_x_in ) ) - 1,  when Skx < 0.
         309             :   !
         310             :   !
         311             :   ! Notes:
         312             :   !
         313             :   ! When F_x = 0 (which can only happen when Skx = 0), mu_x_1 = mu_x_2, and
         314             :   ! mixt_frac = ( zeta_x + 1 ) / ( zeta_x + 2 ).  When these equations are
         315             :   ! substituted into the equations for sigma_x_1 and sigma_x_2, the result is
         316             :   ! sigma_x_1 = sigma_x_2 = sqrt( <x'^2> ).  This means that the distribution
         317             :   ! becomes a single Gaussian when F_x = 0 (and Skx = 0).  This happens
         318             :   ! regardless of the value of zeta_x.
         319             :   !
         320             :   ! The equations for the PDF component means and standard deviations can also
         321             :   ! be written as:
         322             :   !
         323             :   ! mu_x_1 = <x> + sqrt( F_x * ( ( 1 - mixt_frac ) / mixt_frac ) * <x'^2> )
         324             :   !                * sgn( <w'x'> );
         325             :   !
         326             :   ! mu_x_2 = <x> - sqrt( F_x * ( mixt_frac / ( 1 - mixt_frac ) ) * <x'^2> )
         327             :   !                * sgn( <w'x'> );
         328             :   !
         329             :   ! sigma_x_1 = sqrt( coef_sigma_x_1_sqd * <x'^2> ); and
         330             :   !
         331             :   ! sigma_x_2 = sqrt( coef_sigma_x_2_sqd * <x'^2> ); where
         332             :   !
         333             :   ! coef_sigma_x_1_sqd = ( ( zeta_x + 1 ) * ( 1 - F_x ) )
         334             :   !                      / ( ( zeta_x + 2 ) * mixt_frac ); and
         335             :   !
         336             :   ! coef_sigma_x_2_sqd = ( 1 - F_x ) / ( ( zeta_x + 2 ) * ( 1 - mixt_frac ) ).
         337             :   !
         338             :   ! The above equations can be substituted into an equation for a variable that
         339             :   ! has been derived by integrating over the PDF.  Many variables like this are
         340             :   ! used in parts of the predictive equation set.  These substitutions allow
         341             :   ! some terms to solved implicitly or semi-implicitly in the predictive
         342             :   ! equations.
         343             :   !
         344             :   !
         345             :   ! Brian Griffin; September 2017.
         346             :   !
         347             :   !=============================================================================
         348           0 :   function calc_mixture_fraction( nz, Skx, F_x, zeta_x, sgn_wpxp ) &
         349           0 :   result( mixt_frac )
         350             :
         351             :     ! Description:
         352             :     ! Calculates mixture fraction.
         353             :
         354             :     ! References:
         355             :     ! Griffin and Larson (2018)
         356             :     !-----------------------------------------------------------------------
         357             :
         358             :     use constants_clubb, only: &
         359             :         thirty_six, & ! Constant(s)
         360             :         eighteen,   &
         361             :         twelve,     &
         362             :         six,        &
         363             :         four,       &
         364             :         three,      &
         365             :         two,        &
         366             :         one,        &
         367             :         zero,       &
         368             :         fstderr
         369             :
         370             :     use clubb_precision, only: &
         371             :         core_rknd    ! Variable(s)
         372             :
         373             :     implicit none
         374             :
         375             :     integer, intent(in) :: &
         376             :       nz
         377             :
         378             :     ! Input Variables
         379             :     real( kind = core_rknd ), dimension(nz), intent(in) :: &
         380             :       Skx,      & ! Skewness of x                                            [-]
         381             :       F_x,      & ! Parameter for the spread of the PDF component means of x [-]
         382             :       zeta_x,   & ! Parameter for the PDF component variances of x           [-]
         383             :       sgn_wpxp    ! Sign of the covariance of w and x                        [-]
         384             :
         385             :     ! Return Variable
         386             :     real( kind = core_rknd ), dimension(nz) :: &
         387             :       mixt_frac    ! Mixture fraction    [-]
         388             :
         389             :     ! Local Variable
         390             :     ! Flag that turns off when conditions aren't right for calculating mixt_frac
         391             :     logical, dimension(nz) :: &
         392           0 :       l_calc_mixt_frac
         393             :
         394             :
         395             :     ! Initialize l_calc_mixt_frac
         396           0 :     l_calc_mixt_frac = .true.
         397             :
         398             :     ! Calculate mixture fraction, which is the weight of the 1st PDF component.
         399             :     ! The 2nd PDF component has a weight of 1 - mixt_frac.
         400           0 :     where ( F_x > zero )
         401             :
         402             :        mixt_frac &
         403             :        = ( four * F_x**3 &
         404             :            + eighteen * F_x &
         405             :              * ( zeta_x + one ) * ( one - F_x ) / ( zeta_x + two ) &
         406             :            + six * F_x**2 * ( one - F_x ) / ( zeta_x + two ) &
         407             :            + Skx**2 &
         408             :            - Skx * sgn_wpxp * sqrt( four * F_x**3 &
         409             :                                     + twelve * F_x**2 * ( one - F_x ) &
         410             :                                     + thirty_six * F_x &
         411             :                                       * ( zeta_x + one ) * ( one - F_x )**2 &
         412             :                                       / ( zeta_x + two )**2 &
         413             :                                     + Skx**2 ) ) &
         414             :          / ( two * F_x * ( F_x - three )**2 + two * Skx**2 )
         415             :
         416             :     elsewhere ! F_x = 0
         417             :
         418             :        where ( abs( Skx ) > zero )
         419             :
         420             :           l_calc_mixt_frac = .false.
         421             :
         422             :        elsewhere ! Skx = 0
         423             :
         424             :           mixt_frac = ( zeta_x + one ) / ( zeta_x + two )
         425             :
         426             :        endwhere ! | Skx | > 0
         427             :
         428             :     endwhere ! F_x > 0
         429             :
         430             :
         431           0 :     if ( any( .not. l_calc_mixt_frac ) ) then
         432           0 :        write(fstderr,*) "Mixture fraction cannot be calculated."
         433             :        write(fstderr,*) "The value of F_x must be greater than 0 when " &
         434           0 :                         // "| Skx | > 0."
         435           0 :        error stop
         436             :     endif ! any( .not. l_valid_mixt_frac )
         437             :
         438             :
         439           0 :     return
         440             :
         441           0 :   end function calc_mixture_fraction
         442             :
         443             :   !=============================================================================
         444           0 :   subroutine calc_setter_var_params( nz, xm, xp2, Skx, sgn_wpxp,    & ! In
         445           0 :                                      F_x, zeta_x,               & ! In
         446           0 :                                      mu_x_1, mu_x_2, sigma_x_1, & ! Out
         447           0 :                                      sigma_x_2, mixt_frac,      & ! Out
         448           0 :                                      coef_sigma_x_1_sqd,        & ! Out
         449           0 :                                      coef_sigma_x_2_sqd         ) ! Out
         450             :
         451             :     ! Description:
         452             :     ! Calculates the PDF component means, the PDF component standard deviations,
         453             :     ! and the mixture fraction for the variable that sets the PDF.
         454             :
         455             :     ! References:
         456             :     ! Griffin and Larson (2018)
         457             :     !-----------------------------------------------------------------------
         458             :
         459             :     use constants_clubb, only: &
         460             :         two, & ! Variable(s)
         461             :         one
         462             :
         463             :     use clubb_precision, only: &
         464             :         core_rknd    ! Variable(s)
         465             :
         466             :     implicit none
         467             :
         468             :     integer, intent(in) :: &
         469             :       nz
         470             :
         471             :     ! Input Variables
         472             :     real( kind = core_rknd ), dimension(nz), intent(in) :: &
         473             :       xm,       & ! Mean of x (overall)                             [units vary]
         474             :       xp2,      & ! Variance of x (overall)                     [(units vary)^2]
         475             :       Skx,      & ! Skewness of x                                            [-]
         476             :       sgn_wpxp, & ! Sign of the covariance of w and x (overall)              [-]
         477             :       F_x,      & ! Parameter for the spread of the PDF component means of x [-]
         478             :       zeta_x      ! Parameter for the PDF component variances of x           [-]
         479             :
         480             :     ! Output Variables
         481             :     real( kind = core_rknd ), dimension(nz), intent(out) :: &
         482             :       mu_x_1,    & ! Mean of x (1st PDF component)                  [units vary]
         483             :       mu_x_2,    & ! Mean of x (2nd PDF component)                  [units vary]
         484             :       sigma_x_1, & ! Standard deviation of x (1st PDF component)    [units vary]
         485             :       sigma_x_2, & ! Standard deviation of x (2nd PDF component)    [units vary]
         486             :       mixt_frac    ! Mixture fraction                                        [-]
         487             :
         488             :     real( kind = core_rknd ), dimension(nz), intent(out) :: &
         489             :       coef_sigma_x_1_sqd, & ! sigma_x_1^2 = coef_sigma_x_1_sqd * <x'^2>      [-]
         490             :       coef_sigma_x_2_sqd    ! sigma_x_2^2 = coef_sigma_x_2_sqd * <x'^2>      [-]
         491             :
         492             :
         493             :     ! Calculate the mixture fraction.
         494           0 :     mixt_frac = calc_mixture_fraction( nz, Skx, F_x, zeta_x, sgn_wpxp )
         495             :
         496             :     ! Calculate the mean of x in the 1st PDF component.
         497             :     mu_x_1 = xm + sqrt( F_x * ( ( one - mixt_frac ) / mixt_frac ) * xp2 ) &
         498           0 :                   * sgn_wpxp
         499             :
         500             :     ! Calculate the mean of x in the 2nd PDF component.
         501           0 :     mu_x_2 = xm - ( mixt_frac / ( one - mixt_frac ) ) * ( mu_x_1 - xm )
         502             :
         503             :     ! Calculate the standard deviation of x in the 1st PDF component.
         504             :     ! sigma_x_1 = sqrt( ( ( zeta_x + 1 ) * ( 1 - F_x ) )
         505             :     !                   / ( ( zeta_x + 2 ) * mixt_frac ) * <x'^2> )
         506             :     coef_sigma_x_1_sqd = ( ( zeta_x + one ) * ( one - F_x ) ) &
         507           0 :                          / ( ( zeta_x + two ) * mixt_frac )
         508             :
         509           0 :     sigma_x_1 = sqrt( coef_sigma_x_1_sqd * xp2 )
         510             :
         511             :     ! Calculate the standard deviation of x in the 2nd PDF component.
         512             :     ! sigma_x_2 = sqrt( ( mixt_frac * sigma_x_1^2 )
         513             :     !                   / ( ( 1 - mixt_frac ) * ( 1 + zeta_x ) ) )
         514             :     !           = sqrt( ( 1 - F_x )
         515             :     !                   / ( ( zeta_x + 2 ) * ( 1 - mixt_frac ) ) * <x'^2> )
         516             :     coef_sigma_x_2_sqd = ( one - F_x ) &
         517           0 :                          / ( ( zeta_x + two ) * ( one - mixt_frac ) )
         518             :
         519           0 :     sigma_x_2 = sqrt( coef_sigma_x_2_sqd * xp2 )
         520             :
         521             :
         522           0 :     return
         523             :
         524             :   end subroutine calc_setter_var_params
         525             :
         526             :   !=============================================================================
         527             :   !
         528             :   ! DESCRIPTION OF THE METHOD FOR EACH RESPONDING VARIABLE
         529             :   ! ======================================================
         530             :   !
         531             :   ! In order to find equations for the four PDF parameters for each responding
         532             :   ! variable, which are mu_x_1, mu_x_2, sigma_x_1, and sigma_x_2 (where x stands
         533             :   ! for a responding variable here), four equations are needed.  These four
         534             :   ! equations are the equations for <x>, <x'^2>, and <x'^3> as found by
         535             :   ! integrating over the PDF.  Additionally, one more equation, which involves
         536             :   ! a tunable parameter F_x, and which is used to help control the spread of the
         537             :   ! PDF component means, is used in this equation set.  The four equations are:
         538             :   !
         539             :   ! <x> = mixt_frac * mu_x_1 + ( 1 - mixt_frac ) * mu_x_2;
         540             :   !
         541             :   ! <x'^2> = mixt_frac * ( ( mu_x_1 - <x> )^2 + sigma_x_1^2 )
         542             :   !          + ( 1 - mixt_frac ) * ( ( mu_x_2 - <x> )^2 + sigma_x_2^2 );
         543             :   !
         544             :   ! <x'^3> = mixt_frac * ( mu_x_1 - <x> )
         545             :   !                    * ( ( mu_x_1 - <x> )^2 + 3 * sigma_x_1^2 )
         546             :   !          + ( 1 - mixt_frac ) * ( mu_x_2 - <x> )
         547             :   !                              * ( ( mu_x_2 - <x> )^2 + 3 * sigma_x_2^2 ); and
         548             :   !
         549             :   ! mu_x_1 - <x> = sqrt(F_x) * ( sqrt( 1 - mixt_frac ) / sqrt( mixt_frac ) )
         550             :   !                * sqrt( <x'^2> ) * sgn( <w'x'> );
         551             :   !
         552             :   ! where 0 <= F_x <= 1, and where sgn( <w'x'> ) is given by:
         553             :   !
         554             :   ! sgn( <w'x'> ) = |  1, when <w'x'> >= 0;
         555             :   !                 | -1, when <w'x'> < 0.
         556             :   !
         557             :   ! The resulting equations for the four PDF parameters are:
         558             :   !
         559             :   ! mu_x_1 = <x> + sqrt( F_x * ( ( 1 - mixt_frac ) / mixt_frac ) * <x'^2> )
         560             :   !                * sgn( <w'x'> );
         561             :   !
         562             :   ! mu_x_2 = <x> - ( mixt_frac / ( 1 - mixt_frac ) ) * ( mu_x_1 - <x> );
         563             :   !
         564             :   ! sigma_x_1^2
         565             :   ! = ( ( sqrt( mixt_frac * ( 1 - mixt_frac ) ) * Skx * sgn( <w'x'> )
         566             :   !       - ( 1 + mixt_frac ) * F_x^1.5 + 3 * mixt_frac * sqrt( F_x ) )
         567             :   !     / ( 3 * mixt_frac * sqrt( F_x ) ) )
         568             :   !   * <x'^2>; and
         569             :   !
         570             :   ! sigma_x_2^2 = ( ( 1 - F_x ) / ( 1 - mixt_frac ) ) * <x'^2>
         571             :   !               - ( mixt_frac / ( 1 - mixt_frac ) ) * sigma_x_1^2;
         572             :   !
         573             :   ! where Skx is the skewness of x, and Skx = <x'^3> / <x'^2>^(3/2).
         574             :   !
         575             :   !
         576             :   ! Special case:
         577             :   !
         578             :   ! When Skx = 0 and F_x = 0, the equations for sigma_x_1^2 and sigma_x_2^2 are
         579             :   ! both undefined.  The equations for sigma_x_1^2 and sigma_x_2^2 in this
         580             :   ! scenario can be derived by using the above equations for sigma_x_1^2 and
         581             :   ! sigma_x_2^2 and then setting Skx = 0.  The resulting equation for
         582             :   ! sigma_x_1^2 becomes:
         583             :   !
         584             :   ! sigma_x_1^2
         585             :   ! = ( ( - ( 1 + mixt_frac ) * F_x^1.5 + 3 * mixt_frac * sqrt( F_x ) )
         586             :   !     / ( 3 * mixt_frac * sqrt( F_x ) ) )
         587             :   !   * <x'^2>.
         588             :   !
         589             :   ! All of the terms in the numerator and denominator contain a sqrt( F_x ),
         590             :   ! so this equation can be rewritten as:
         591             :   !
         592             :   ! sigma_x_1^2
         593             :   ! = ( ( - ( 1 + mixt_frac ) * F_x + 3 * mixt_frac ) / ( 3 * mixt_frac ) )
         594             :   !   * <x'^2>.
         595             :   !
         596             :   ! Now setting F_x = 0, the equation becomes:
         597             :   !
         598             :   ! sigma_x_1^2 = ( ( 3 * mixt_frac ) / ( 3 * mixt_frac ) ) * <x'^2>;
         599             :   !
         600             :   ! which can be rewritten as:
         601             :   !
         602             :   ! sigma_x_1^2 = <x'^2>.
         603             :   !
         604             :   ! Substituting the equation for sigma_x_1^2 into the equation for sigma_x_2^2,
         605             :   ! and also setting F_x = 0, the equation for sigma_x_2^2 becomes:
         606             :   !
         607             :   ! sigma_x_2^2 = ( 1 / ( 1 - mixt_frac ) ) * <x'^2>
         608             :   !               - ( mixt_frac / ( 1 - mixt_frac ) ) * <x'^2>;
         609             :   !
         610             :   ! which can be rewritten as:
         611             :   !
         612             :   ! sigma_x_2^2
         613             :   ! = ( ( 1 / ( 1 - mixt_frac ) ) - ( mixt_frac / ( 1 - mixt_frac ) ) )
         614             :   !   * <x'^2>.
         615             :   !
         616             :   ! This equation becomes:
         617             :   !
         618             :   ! sigma_x_2^2 = ( ( 1 - mixt_frac ) / ( 1 - mixt_frac ) ) * <x'^2>;
         619             :   !
         620             :   ! which can be rewritten as:
         621             :   !
         622             :   ! sigma_x_2^2 = <x'^2>.
         623             :   !
         624             :   ! When F_x = 0, Skx must have a value of 0 in order for the PDF to function
         625             :   ! correctly.  When F_x = 0, mu_x_1 = mu_x_2.  When the two PDF component means
         626             :   ! are equal to each other (and to the overall mean, <x>), the only value of
         627             :   ! Skx that can be represented is a value of 0.  The equations that place
         628             :   ! limits on F_x for a responding variable (below) calculate the minimum
         629             :   ! allowable value of F_x to be greater than 0 when | Skx | > 0.
         630             :   !
         631             :   ! The value of F_x should be set as a function of Skx.  The value F_x should
         632             :   ! go toward 0 as | Skx | (or Skx^2) goes toward 0.  The value of F_x should
         633             :   ! go toward 1 as | Skx | (or Skx^2) goes to infinity.  However, the value of
         634             :   ! F_x must also be between the minimum and maximum allowable values of F_x for
         635             :   ! a responding variable (below).
         636             :   !
         637             :   !
         638             :   ! Tunable parameter:
         639             :   !
         640             :   ! F_x:  This parameter controls the spread of the PDF component means.  The
         641             :   !       range of this parameter is 0 <= F_x <= 1.  When F_x = 0, the two PDF
         642             :   !       component means (mu_x_1 and mu_x_2) are equal to each other (and Skx
         643             :   !       must equal 0).  All of the variance of x is accounted for by the PDF
         644             :   !       component standard deviations (sigma_x_1 and sigma_x_2).  When
         645             :   !       F_x = 1, mu_x_1 and mu_x_2 are spread as far apart as they can be.
         646             :   !       Both PDF component standard deviations (sigma_x_1 and sigma_x_2) are
         647             :   !       equal to 0, and all of the variance of x is accounted for by the
         648             :   !       spread of the PDF component means.
         649             :   !
         650             :   !
         651             :   ! Limits on F_x:
         652             :   !
         653             :   ! Since the PDF parameters for this variable need to work with the mixture
         654             :   ! fraction that has been provided by the setting variable, the method does
         655             :   ! not work for all values of F_x and Skx.  However, the limits of Skx and F_x
         656             :   ! can always be calculated.  The limits are based on keeping the values of
         657             :   ! sigma_x_1 and sigma_x_2 greater than or equal to 0.  The equation for
         658             :   ! keeping the value of sigma_x_1 greater than or equal to 0 is:
         659             :   !
         660             :   ! - ( 1 + mixt_frac ) * sqrt( F_x )^3 + 3 * mixt_frac * sqrt( F_x )
         661             :   ! + sqrt( mixt_frac * ( 1 - mixt_frac ) ) * Skx * sgn( <w'x'> ) >= 0.
         662             :   !
         663             :   ! The roots of sqrt( F_x ) can be solved by an equation of the form:
         664             :   !
         665             :   ! A * sqrt( F_x )^3 + B * sqrt( F_x )^2 + C * sqrt( F_x ) + D = 0;
         666             :   !
         667             :   ! where:
         668             :   !
         669             :   ! A = - ( 1 + mixt_frac );
         670             :   ! B = 0;
         671             :   ! C = 3 * mixt_frac; and
         672             :   ! D = sqrt( mixt_frac * ( 1 - mixt_frac ) ) * Skx * sgn( <w'x'> ).
         673             :   !
         674             :   ! The equation for keeping the value of sigma_x_2 greater than or equal to 0
         675             :   ! is:
         676             :   !
         677             :   ! - ( 2 - mixt_frac ) * sqrt( F_x )^3 + 3 * ( 1 - mixt_frac ) * sqrt( F_x )
         678             :   ! - sqrt( mixt_frac * ( 1 - mixt_frac ) ) * Skx * sgn( <w'x'> ) >= 0.
         679             :   !
         680             :   ! The roots of sqrt( F_x ) can be solved by an equation of the form:
         681             :   !
         682             :   ! A * sqrt( F_x )^3 + B * sqrt( F_x )^2 + C * sqrt( F_x ) + D = 0;
         683             :   !
         684             :   ! where:
         685             :   !
         686             :   ! A = - ( 2 - mixt_frac );
         687             :   ! B = 0;
         688             :   ! C = 3 * ( 1 - mixt_frac ); and
         689             :   ! D = - sqrt( mixt_frac * ( 1 - mixt_frac ) ) * Skx * sgn( <w'x'> ).
         690             :   !
         691             :   ! After careful analysis of the above equations, the following properties
         692             :   ! emerge:
         693             :   !
         694             :   ! When Skx * sgn( <w'x'> ) >= 0,
         695             :   !    Skx^2 < 4 * ( 1 - mixt_frac )^2 / ( mixt_frac * ( 2 - mixt_frac ) )
         696             :   !    is required; and
         697             :   ! when Skx * sgn( <w'x'> ) < 0,
         698             :   !    Skx^2 < 4 * mixt_frac^2 / ( 1 - mixt_frac^2 ) is required.
         699             :   !
         700             :   ! Whenever Skx^2 exceeds these limits, Skx must be limited (preserving its
         701             :   ! sign) in order to have any value of F_x that will work in the equation set.
         702             :   !
         703             :   ! When Skx is found to be within the above limits (or after it has been
         704             :   ! limited to fall within its limits), the range of valid values of F_x can be
         705             :   ! found according to the following:
         706             :   !
         707             :   ! When Skx * sgn( <w'x'> ) >= 0:
         708             :   !
         709             :   !     When  4 * mixt_frac^2 / ( 1 - mixt_frac^2 )  <  Skx^2
         710             :   !           <  4 * ( 1 - mixt_frac )^2 / ( mixt_frac * ( 2 - mixt_frac ) ):
         711             :   !
         712             :   !          Minimum sqrt( F_x ):  2nd root (middle-valued root; also smallest
         713             :   !                                positive) of the second equation (sigma_x_2
         714             :   !                                based).
         715             :   !
         716             :   !          Maximum sqrt( F_x ):  Minimum of the largest root of the second
         717             :   !                                equation (sigma_x_2 based) and the only* root
         718             :   !                                of the first equation (sigma_x_1 based).
         719             :   !
         720             :   !     When  Skx^2  <=  4 * mixt_frac^2 / ( 1 - mixt_frac^2 ):
         721             :   !
         722             :   !          Minimum sqrt( F_x ):  2nd root (middle-valued root; also smallest
         723             :   !                                positive) of the second equation (sigma_x_2
         724             :   !                                based).
         725             :   !
         726             :   !          Maximum sqrt( F_x ):  Minimum of the largest root of the second
         727             :   !                                equation (sigma_x_2 based) and the largest
         728             :   !                                root of the first equation (sigma_x_1 based).
         729             :   !
         730             :   ! When Skx * sgn( <w'x'> ) < 0:
         731             :   !
         732             :   !     When  4 * ( 1 - mixt_frac )^2 / ( mixt_frac * ( 2 - mixt_frac ) )
         733             :   !           <  Skx^2  <  4 * mixt_frac^2 / ( 1 - mixt_frac^2 ):
         734             :   !
         735             :   !          Minimum sqrt( F_x ):  2nd root (middle-valued root; also smallest
         736             :   !                                positive) of the first equation (sigma_x_1
         737             :   !                                based).
         738             :   !
         739             :   !          Maximum sqrt( F_x ):  Minimum of the largest root of the first
         740             :   !                                equation (sigma_x_1 based) and the only* root
         741             :   !                                of the second equation (sigma_x_2 based).
         742             :   !
         743             :   !     When  Skx^2
         744             :   !           <=  4 * ( 1 - mixt_frac )^2 / ( mixt_frac * ( 2 - mixt_frac ) ):
         745             :   !
         746             :   !          Minimum sqrt( F_x ):  2nd root (middle-valued root; also smallest
         747             :   !                                positive) of the first equation (sigma_x_1
         748             :   !                                based).
         749             :   !
         750             :   !          Maximum sqrt( F_x ):  Minimum of the largest root of the first
         751             :   !                                equation (sigma_x_1 based) and the largest
         752             :   !                                root of the second equation (sigma_x_2
         753             :   !                                based).
         754             :   !
         755             :   ! Here, "only* root" means the the only root that isn't a complex root.
         756             :   !
         757             :   ! The value of sqrt( F_x ) is also limited with a minimum of 0 and a maximum
         758             :   ! of 1.  The minimum and maximum allowable values of F_x are found by taking
         759             :   ! the square of the minimum and maximum allowable values of sqrt( F_x ),
         760             :   ! respectively.
         761             :   !
         762             :   !
         763             :   ! Notes:
         764             :   !
         765             :   ! When F_x = 0 (which can only happen when Skx = 0), mu_x_1 = mu_x_2, and
         766             :   ! sigma_x_1 = sigma_x_2 = sqrt( <x'^2> ).  This means that the distribution
         767             :   ! becomes a single Gaussian when F_x = 0 (and Skx = 0).
         768             :   !
         769             :   ! The equations for the PDF component means and standard deviations can also
         770             :   ! be written as:
         771             :   !
         772             :   ! mu_x_1 = <x> + sqrt( F_x * ( ( 1 - mixt_frac ) / mixt_frac ) * <x'^2> )
         773             :   !                * sgn( <w'x'> );
         774             :   !
         775             :   ! mu_x_2 = <x> - sqrt( F_x * ( mixt_frac / ( 1 - mixt_frac ) ) * <x'^2> )
         776             :   !                * sgn( <w'x'> );
         777             :   !
         778             :   ! sigma_x_1 = sqrt( coef_sigma_x_1_sqd * <x'^2> ); and
         779             :   !
         780             :   ! sigma_x_2 = sqrt( coef_sigma_x_2_sqd * <x'^2> ); where
         781             :   !
         782             :   ! coef_sigma_x_1_sqd
         783             :   ! = ( sqrt( mixt_frac * ( 1 - mixt_frac ) ) * Skx * sgn( <w'x'> )
         784             :   !     - ( 1 + mixt_frac ) * F_x^1.5 + 3 * mixt_frac * sqrt( F_x ) )
         785             :   !   / ( 3 * mixt_frac * sqrt( F_x ) )
         786             :   ! = sqrt( mixt_frac * ( 1 - mixt_frac ) ) * Skx * sgn( <w'x'> )
         787             :   !   / ( 3 * mixt_frac * sqrt( F_x ) )
         788             :   !   - ( 1 + mixt_frac ) * F_x / ( 3 * mixt_frac )
         789             :   !   + 1; and
         790             :   !
         791             :   ! coef_sigma_x_2_sqd
         792             :   ! = ( 1 - F_x ) / ( 1 - mixt_frac )
         793             :   !   - mixt_frac / ( 1 - mixt_frac )
         794             :   !     * ( sqrt( mixt_frac * ( 1 - mixt_frac ) ) * Skx * sgn( <w'x'> )
         795             :   !         / ( 3 * mixt_frac * sqrt( F_x ) )
         796             :   !         - ( 1 + mixt_frac ) * F_x / ( 3 * mixt_frac )
         797             :   !         + 1 )
         798             :   ! = ( ( 1 - F_x ) - mixt_frac * coef_sigma_x_1_sqd ) / ( 1 - mixt_frac ).
         799             :   !
         800             :   ! The above equations can be substituted into an equation for a variable that
         801             :   ! has been derived by integrating over the PDF.  Many variables like this are
         802             :   ! used in parts of the predictive equation set.  These substitutions allow
         803             :   ! some terms to solved implicitly or semi-implicitly in the predictive
         804             :   ! equations.
         805             :   !
         806             :   !
         807             :   ! Brian Griffin; September 2017.
         808             :   !
         809             :   !=============================================================================
         810           0 :   subroutine calc_responder_params( nz, xm, xp2, Skx, sgn_wpxp,       & ! In
         811           0 :                                     F_x, mixt_frac,               & ! In
         812           0 :                                     mu_x_1, mu_x_2,               & ! Out
         813           0 :                                     sigma_x_1_sqd, sigma_x_2_sqd, & ! Out
         814           0 :                                     coef_sigma_x_1_sqd,           & ! Out
         815           0 :                                     coef_sigma_x_2_sqd            ) ! Out
         816             :
         817             :     ! Description:
         818             :     ! Calculates the PDF component means and the PDF component standard
         819             :     ! deviations for a responding variable (a variable that is not used to set
         820             :     ! the mixture fraction).
         821             :
         822             :     ! References:
         823             :     ! Griffin and Larson (2018)
         824             :     !-----------------------------------------------------------------------
         825             :
         826             :     use constants_clubb, only: &
         827             :         three, & ! Variable(s)
         828             :         one,   &
         829             :         zero
         830             :
         831             :     use clubb_precision, only: &
         832             :         core_rknd    ! Variable(s)
         833             :
         834             :     implicit none
         835             :
         836             :     integer, intent(in) :: &
         837             :       nz
         838             :
         839             :     ! Input Variables
         840             :     real( kind = core_rknd ), dimension(nz), intent(in) :: &
         841             :       xm,       & ! Mean of x (overall)                             [units vary]
         842             :       xp2,      & ! Variance of x (overall)                     [(units vary)^2]
         843             :       Skx,      & ! Skewness of x                                            [-]
         844             :       sgn_wpxp, & ! Sign of the covariance of w and x (overall)              [-]
         845             :       F_x,      & ! Parameter for the spread of the PDF component means of x [-]
         846             :       mixt_frac   ! Mixture fraction                                         [-]
         847             :
         848             :     ! Output Variables
         849             :     real( kind = core_rknd ), dimension(nz), intent(out) :: &
         850             :       mu_x_1,        & ! Mean of x (1st PDF component)              [units vary]
         851             :       mu_x_2,        & ! Mean of x (2nd PDF component)              [units vary]
         852             :       sigma_x_1_sqd, & ! Variance of x (1st PDF component)      [(units vary)^2]
         853             :       sigma_x_2_sqd    ! Variance of x (2nd PDF component)      [(units vary)^2]
         854             :
         855             :     real( kind = core_rknd ), dimension(nz), intent(out) :: &
         856             :       coef_sigma_x_1_sqd, & ! sigma_x_1^2 = coef_sigma_x_1_sqd * <x'^2>      [-]
         857             :       coef_sigma_x_2_sqd    ! sigma_x_2^2 = coef_sigma_x_2_sqd * <x'^2>      [-]
         858             :
         859             :
         860           0 :     where ( F_x > zero )
         861             :
         862             :        ! Calculate the mean of x in the 1st PDF component.
         863             :        mu_x_1 = xm + sqrt( F_x * ( ( one - mixt_frac ) / mixt_frac ) * xp2 ) &
         864             :                      * sgn_wpxp
         865             :
         866             :        ! Calculate the mean of x in the 2nd PDF component.
         867             :        mu_x_2 = xm - ( mixt_frac / ( one - mixt_frac ) ) * ( mu_x_1 - xm )
         868             :
         869             :        ! Calculate the variance of x in the 1st PDF component.
         870             :        ! sigma_x_1^2
         871             :        ! = ( ( sqrt( mixt_frac * ( 1 - mixt_frac ) ) * Skx * sgn( <w'x'> )
         872             :        !       - ( 1 + mixt_frac ) * F_x^1.5 + 3 * mixt_frac * sqrt( F_x ) )
         873             :        !     / ( 3 * mixt_frac * sqrt( F_x ) ) ) * <x'^2>
         874             :        ! = ( sqrt( mixt_frac * ( 1 - mixt_frac ) ) * Skx * sgn( <w'x'> )
         875             :        !     / ( 3 * mixt_frac * sqrt( F_x ) )
         876             :        !     - ( 1 + mixt_frac ) * F_x / ( 3 * mixt_frac )
         877             :        !     + 1 ) * <x'^2>
         878             :        coef_sigma_x_1_sqd &
         879             :        = sqrt( mixt_frac * ( one - mixt_frac ) ) * Skx * sgn_wpxp &
         880             :          / ( three * mixt_frac * sqrt( F_x ) ) &
         881             :          - ( one + mixt_frac ) * F_x / ( three * mixt_frac ) &
         882             :          + one
         883             :
         884             :        sigma_x_1_sqd = coef_sigma_x_1_sqd * xp2
         885             :
         886             :        ! Calculate the variance of x in the 2nd PDF component.
         887             :        ! sigma_x_2^2
         888             :        ! = ( ( 1 - F_x ) / ( 1 - mixt_frac )
         889             :        !     - mixt_frac / ( 1 - mixt_frac )
         890             :        !       * ( sqrt( mixt_frac * ( 1 - mixt_frac ) ) * Skx * sgn( <w'x'> )
         891             :        !           / ( 3 * mixt_frac * sqrt( F_x ) )
         892             :        !           - ( 1 + mixt_frac ) * F_x / ( 3 * mixt_frac )
         893             :        !           + 1 ) ) * <x'^2>
         894             :        ! = ( ( ( 1 - F_x ) - mixt_frac * coef_sigma_x_1_sqd )
         895             :        !     / ( 1 - mixt_frac ) ) * <x'^2>
         896             :        coef_sigma_x_2_sqd &
         897             :        = ( ( one - F_x ) - mixt_frac * coef_sigma_x_1_sqd ) &
         898             :          / ( one - mixt_frac )
         899             :
         900             :        sigma_x_2_sqd = coef_sigma_x_2_sqd * xp2
         901             :
         902             :     elsewhere ! F_x = 0
         903             :
         904             :        ! When F_x has a value of 0, the PDF becomes a single Gaussian.  This
         905             :        ! only works when Skx = 0.  However, when Skx /= 0, the value of min_F_x
         906             :        ! is greater than 0, preventing a problem where F_x = 0 but | Skx | > 0.
         907             :        mu_x_1 = xm
         908             :        mu_x_2 = xm
         909             :        sigma_x_1_sqd = xp2
         910             :        sigma_x_2_sqd = xp2
         911             :        coef_sigma_x_1_sqd = one
         912             :        coef_sigma_x_2_sqd = one
         913             :
         914             :     endwhere ! F_x > 0
         915             :
         916             :
         917           0 :     return
         918             :
         919             :   end subroutine calc_responder_params
         920             :
         921             :   !=============================================================================
         922           0 :   subroutine calc_limits_F_x_responder( nz, mixt_frac, Skx, sgn_wpxp,  & ! In
         923           0 :                                         max_Skx2_pos_Skx_sgn_wpxp, & ! In
         924           0 :                                         max_Skx2_neg_Skx_sgn_wpxp, & ! In
         925           0 :                                         min_F_x, max_F_x )           ! Out
         926             :
         927             :     ! Description:
         928             :     ! Calculates the minimum and maximum allowable values for F_x for a
         929             :     ! responding variable.
         930             :
         931             :     ! References:
         932             :     !-----------------------------------------------------------------------
         933             :
         934             :     use constants_clubb, only: &
         935             :         three, & ! Variable(s)
         936             :         two,   &
         937             :         one,   &
         938             :         zero
         939             :
         940             :     use calc_roots, only: &
         941             :         cubic_solve    ! Procedure(s)
         942             :
         943             :     use clubb_precision, only: &
         944             :         core_rknd    ! Variable(s)
         945             :
         946             :     implicit none
         947             :
         948             :     integer, intent(in) :: &
         949             :       nz
         950             :
         951             :     ! Input Variables
         952             :     real( kind = core_rknd ), dimension(nz), intent(in) :: &
         953             :       mixt_frac, & ! Mixture fraction                 [-]
         954             :       Skx,       & ! Skewness of x                    [-]
         955             :       sgn_wpxp     ! Sign of covariance of w and x    [-]
         956             :
         957             :     real( kind = core_rknd ), dimension(nz), intent(in) :: &
         958             :       max_Skx2_pos_Skx_sgn_wpxp, & ! Maximum Skx^2 when Skx*sgn(<w'x'>) >= 0 [-]
         959             :       max_Skx2_neg_Skx_sgn_wpxp    ! Maximum Skx^2 when Skx*sgn(<w'x'>) < 0  [-]
         960             :
         961             :     ! Output Variables
         962             :     real( kind = core_rknd ), dimension(nz), intent(out) :: &
         963             :       min_F_x, & ! Minimum allowable value of F_x    [-]
         964             :       max_F_x    ! Maximum allowable value of F_x    [-]
         965             :
         966             :     ! Local Variables
         967             :     real( kind = core_rknd ), dimension(nz) :: &
         968           0 :       coef_A_1, & ! Coef. A in Ax^3 + Bx^2 + Cx + D = 0 (1st PDF comp. lim.) [-]
         969           0 :       coef_B_1, & ! Coef. B in Ax^3 + Bx^2 + Cx + D = 0 (1st PDF comp. lim.) [-]
         970           0 :       coef_C_1, & ! Coef. C in Ax^3 + Bx^2 + Cx + D = 0 (1st PDF comp. lim.) [-]
         971           0 :       coef_D_1, & ! Coef. D in Ax^3 + Bx^2 + Cx + D = 0 (1st PDF comp. lim.) [-]
         972           0 :       coef_A_2, & ! Coef. A in Ax^3 + Bx^2 + Cx + D = 0 (2nd PDF comp. lim.) [-]
         973           0 :       coef_B_2, & ! Coef. B in Ax^3 + Bx^2 + Cx + D = 0 (2nd PDF comp. lim.) [-]
         974           0 :       coef_C_2, & ! Coef. C in Ax^3 + Bx^2 + Cx + D = 0 (2nd PDF comp. lim.) [-]
         975           0 :       coef_D_2    ! Coef. D in Ax^3 + Bx^2 + Cx + D = 0 (2nd PDF comp. lim.) [-]
         976             :
         977             :     complex( kind = core_rknd ), dimension(nz,3) :: &
         978           0 :       sqrt_F_x_roots_1, & ! Roots of sqrt(F_x) for the sigma_x_1 term    [-]
         979           0 :       sqrt_F_x_roots_2    ! Roots of sqrt(F_x) for the sigma_x_2 term    [-]
         980             :
         981             :     real( kind = core_rknd ), dimension(nz,3) :: &
         982           0 :       sqrt_F_x_roots_1_sorted, & ! Sorted roots of sqrt(F_x): sigma_x_1 term [-]
         983           0 :       sqrt_F_x_roots_2_sorted    ! Sorted roots of sqrt(F_x): sigma_x_2 term [-]
         984             :
         985             :     real( kind = core_rknd ), dimension(nz) :: &
         986           0 :       min_sqrt_F_x, & ! Minimum allowable value of sqrt(F_x)    [-]
         987           0 :       max_sqrt_F_x    ! Maximum allowable value of sqrt(F_x)    [-]
         988             :
         989             :
         990             :     ! Set up the coefficients in the equation for the limit of sqrt(F_x) based
         991             :     ! on the 1st PDF component standard deviation (sigma_x_1) being greater than
         992             :     ! or equal to 0.  This equation has the form:
         993             :     ! A * sqrt(F_x)^3 + B * sqrt(F_x)^2 + C * sqrt(F_x) + D = 0.
         994           0 :     coef_A_1 = -( one + mixt_frac )
         995           0 :     coef_B_1 = zero
         996           0 :     coef_C_1 = three * mixt_frac
         997           0 :     coef_D_1 = sqrt( mixt_frac * ( one - mixt_frac ) ) * Skx * sgn_wpxp
         998             :
         999             :     ! Solve for the roots (values of sqrt(F_x)) that satisfy the above equation.
        1000             :     sqrt_F_x_roots_1 &
        1001           0 :     = cubic_solve( nz, coef_A_1, coef_B_1, coef_C_1, coef_D_1 )
        1002             :
        1003             :     ! Sort the values of the roots (values of sqrt(F_x)) from smallest to
        1004             :     ! largest.  Ignore any complex component of the roots.  The code below that
        1005             :     ! uses sqrt_F_x_roots_1_sorted already factors the appropriate roots to use
        1006             :     ! into account.
        1007             :     sqrt_F_x_roots_1_sorted &
        1008           0 :     = sort_roots( nz, real( sqrt_F_x_roots_1, kind = core_rknd ) )
        1009             :
        1010             :     ! Set up the coefficients in the equation for the limit of sqrt(F_x) based
        1011             :     ! on the 2nd PDF component standard deviation (sigma_x_2) being greater than
        1012             :     ! or equal to 0.  This equation has the form:
        1013             :     ! A * sqrt(F_x)^3 + B * sqrt(F_x)^2 + C * sqrt(F_x) + D = 0.
        1014           0 :     coef_A_2 = -( two - mixt_frac )
        1015           0 :     coef_B_2 = zero
        1016           0 :     coef_C_2 = three * ( one - mixt_frac )
        1017           0 :     coef_D_2 = -sqrt( mixt_frac * ( one - mixt_frac ) ) * Skx * sgn_wpxp
        1018             :
        1019             :     ! Solve for the roots (values of sqrt(F_x)) that satisfy the above equation.
        1020             :     sqrt_F_x_roots_2 &
        1021           0 :     = cubic_solve( nz, coef_A_2, coef_B_2, coef_C_2, coef_D_2 )
        1022             :
        1023             :     ! Sort the values of the roots (values of sqrt(F_x)) from smallest to
        1024             :     ! largest.  Ignore any complex component of the roots.  The code below that
        1025             :     ! uses sqrt_F_x_roots_2_sorted already factors the appropriate roots to use
        1026             :     ! into account.
        1027             :     sqrt_F_x_roots_2_sorted &
        1028           0 :     = sort_roots( nz, real( sqrt_F_x_roots_2, kind = core_rknd ) )
        1029             :
        1030             :
        1031             :     ! Find the minimum and maximum allowable values of sqrt(F_x) based on Skx
        1032             :     ! and sgn( <w'x'> ).
        1033           0 :     where ( Skx * sgn_wpxp >= zero )
        1034             :
        1035             :        where ( Skx**2 > max_Skx2_neg_Skx_sgn_wpxp )
        1036             :
        1037             :           min_sqrt_F_x = sqrt_F_x_roots_2_sorted(:,2)
        1038             :           max_sqrt_F_x = min( real( sqrt_F_x_roots_1(:,1), kind = core_rknd ), &
        1039             :                               sqrt_F_x_roots_2_sorted(:,3) )
        1040             :
        1041             :        elsewhere ! Skx^2 <= max_Skx2_neg_Skx_sgn_wpxp
        1042             :
        1043             :           min_sqrt_F_x = sqrt_F_x_roots_2_sorted(:,2)
        1044             :           max_sqrt_F_x = min( sqrt_F_x_roots_1_sorted(:,3), &
        1045             :                               sqrt_F_x_roots_2_sorted(:,3) )
        1046             :
        1047             :        endwhere ! Skx**2 > max_Skx2_neg_Skx_sgn_wpxp
        1048             :
        1049             :     elsewhere ! Skx * sgn( <w'x'> ) < 0
        1050             :
        1051             :        where ( Skx**2 > max_Skx2_pos_Skx_sgn_wpxp )
        1052             :
        1053             :           min_sqrt_F_x = sqrt_F_x_roots_1_sorted(:,2)
        1054             :           max_sqrt_F_x = min( real( sqrt_F_x_roots_2(:,1), kind = core_rknd ), &
        1055             :                               sqrt_F_x_roots_1_sorted(:,3) )
        1056             :
        1057             :        elsewhere ! Skx^2 <= max_Skx2_pos_Skx_sgn_wpxp
        1058             :
        1059             :           min_sqrt_F_x = sqrt_F_x_roots_1_sorted(:,2)
        1060             :           max_sqrt_F_x = min( sqrt_F_x_roots_1_sorted(:,3), &
        1061             :                               sqrt_F_x_roots_2_sorted(:,3) )
        1062             :
        1063             :        endwhere ! Skx**2 > max_Skx2_pos_Skx_sgn_wpxp
        1064             :
        1065             :     endwhere ! Skx * sgn( <w'x'> ) >= 0
        1066             :
        1067             :
        1068             :     ! The minimum and maximum are also limited by 0 and 1, respectively.
        1069           0 :     min_sqrt_F_x = max( min_sqrt_F_x, zero )
        1070           0 :     max_sqrt_F_x = min( max_sqrt_F_x, one )
        1071             :
        1072             :     ! The minimum and maximum allowable values for F_x are the squares of the
        1073             :     ! minimum and maximum allowable values for sqrt(F_x).
        1074           0 :     min_F_x = min_sqrt_F_x**2
        1075           0 :     max_F_x = max_sqrt_F_x**2
        1076             :
        1077             :
        1078           0 :     return
        1079             :
        1080             :   end subroutine calc_limits_F_x_responder
        1081             :
        1082             :   !=============================================================================
        1083           0 :   function sort_roots( nz, roots ) &
        1084           0 :   result ( roots_sorted )
        1085             :
        1086             :     ! Description:
        1087             :     ! Sorts roots from smallest to largest.
        1088             :
        1089             :     ! References:
        1090             :     !-----------------------------------------------------------------------
        1091             :
        1092             :     use clubb_precision, only: &
        1093             :         core_rknd    ! Variable(s)
        1094             :
        1095             :     implicit none
        1096             :
        1097             :     integer, intent(in) :: &
        1098             :       nz
        1099             :
        1100             :     ! Input Variable
        1101             :     real( kind = core_rknd ), dimension(nz,3), intent(in) :: &
        1102             :       roots    ! Roots    [-]
        1103             :
        1104             :     ! Return Variable
        1105             :     real( kind = core_rknd ), dimension(nz,3) :: &
        1106             :       roots_sorted    ! Roots sorted from smallest to largest    [-]
        1107             :
        1108             :
        1109           0 :     where ( roots(:,1) <= roots(:,2) .and. roots(:,1) <= roots(:,3) )
        1110             :
        1111             :        ! The value of roots(1) is the smallest root.
        1112             :        roots_sorted(:,1) = roots(:,1)
        1113             :
        1114             :        where ( roots(:,2) <= roots(:,3) )
        1115             :
        1116             :           ! The value of roots(2) is the middle-valued root and the value of
        1117             :           ! roots(3) is the largest root.
        1118             :           roots_sorted(:,2) = roots(:,2)
        1119             :           roots_sorted(:,3) = roots(:,3)
        1120             :
        1121             :        elsewhere ! roots(3) < roots(2)
        1122             :
        1123             :           ! The value of roots(3) is the middle-valued root and the value of
        1124             :           ! roots(2) is the largest root.
        1125             :           roots_sorted(:,2) = roots(:,3)
        1126             :           roots_sorted(:,3) = roots(:,2)
        1127             :
        1128             :        endwhere ! roots(2) <= roots(3)
        1129             :
        1130             :     elsewhere ( roots(:,2) < roots(:,1) .and. roots(:,2) <= roots(:,3) )
        1131             :
        1132             :        ! The value of roots(2) is the smallest root.
        1133             :        roots_sorted(:,1) = roots(:,2)
        1134             :
        1135             :        where ( roots(:,1) <= roots(:,3) )
        1136             :
        1137             :           ! The value of roots(1) is the middle-valued root and the value of
        1138             :           ! roots(3) is the largest root.
        1139             :           roots_sorted(:,2) = roots(:,1)
        1140             :           roots_sorted(:,3) = roots(:,3)
        1141             :
        1142             :        elsewhere ! roots(3) < roots(1)
        1143             :
        1144             :           ! The value of roots(3) is the middle-valued root and the value of
        1145             :           ! roots(1) is the largest root.
        1146             :           roots_sorted(:,2) = roots(:,3)
        1147             :           roots_sorted(:,3) = roots(:,1)
        1148             :
        1149             :        endwhere ! roots(1) <= roots(3)
        1150             :
        1151             :     elsewhere ! roots(3) < roots(1) .and. roots(3) < roots(2)
        1152             :
        1153             :        ! The value of roots(3) is the smallest root.
        1154             :        roots_sorted(:,1) = roots(:,3)
        1155             :
        1156             :        where ( roots(:,1) <= roots(:,2) )
        1157             :
        1158             :           ! The value of roots(1) is the middle-valued root and the value of
        1159             :           ! roots(2) is the largest root.
        1160             :           roots_sorted(:,2) = roots(:,1)
        1161             :           roots_sorted(:,3) = roots(:,2)
        1162             :
        1163             :        elsewhere ! roots(2) < roots(1)
        1164             :
        1165             :           ! The value of roots(2) is the middle-valued root and the value of
        1166             :           ! roots(1) is the largest root.
        1167             :           roots_sorted(:,2) = roots(:,2)
        1168             :           roots_sorted(:,3) = roots(:,1)
        1169             :
        1170             :        endwhere ! roots(1) <= roots(2)
        1171             :
        1172             :     endwhere ! roots(1) <= roots(2) .and. roots(1) <= roots(3)
        1173             :
        1174             :
        1175           0 :     return
        1176             :
        1177           0 :   end function sort_roots
        1178             :
        1179             :   !=============================================================================
        1180           0 :   function calc_coef_wp4_implicit( nz, mixt_frac, F_w, &
        1181           0 :                                    coef_sigma_w_1_sqd, &
        1182           0 :                                    coef_sigma_w_2_sqd ) &
        1183           0 :   result( coef_wp4_implicit )
        1184             :
        1185             :     ! Description:
        1186             :     ! The predictive equation for <w'^3> contains a turbulent advection term of
        1187             :     ! the form:
        1188             :     !
        1189             :     ! - ( 1 / rho_ds ) * d ( rho_ds * <w'^4> ) / dz;
        1190             :     !
        1191             :     ! where z is height, rho_ds is the dry, base-state density, and <w'^4> is
        1192             :     ! calculated by integrating over the PDF.  The equation for <w'^4> is:
        1193             :     !
        1194             :     ! <w'^4> = mixt_frac * ( 3 * sigma_w_1^4
        1195             :     !                        + 6 * ( mu_w_1 - <w> )^2 * sigma_w_1^2
        1196             :     !                        + ( mu_w_1 - <w> )^4 )
        1197             :     !          + ( 1 - mixt_frac ) * ( 3 * sigma_w_2^4
        1198             :     !                                  + 6 * ( mu_w_2 - <w> )^2 * sigma_w_2^2
        1199             :     !                                  + ( mu_w_2 - <w> )^4 ).
        1200             :     !
        1201             :     ! The following substitutions are made into the above equation:
        1202             :     !
        1203             :     ! mu_w_1 - <w> = sqrt(F_w) * sqrt( ( 1 - mixt_frac ) / mixt_frac )
        1204             :     !                * sqrt( <w'^2> );
        1205             :     !
        1206             :     ! mu_w_2 - <w> = - sqrt(F_w) * sqrt( mixt_frac / ( 1 - mixt_frac ) )
        1207             :     !                  * sqrt( <w'^2> );
        1208             :     !
        1209             :     ! sigma_w_1 = sqrt( coef_sigma_w_1_sqd * <w'^2> ); and
        1210             :     !
        1211             :     ! sigma_w_2 = sqrt( coef_sigma_w_2_sqd * <w'^2> ).
        1212             :     !
        1213             :     ! When w is the setting variable, coef_sigma_w_1_sqd and coef_sigma_w_2_sqd
        1214             :     ! are given by:
        1215             :     !
        1216             :     ! coef_sigma_w_1_sqd = ( ( zeta_w + 1 ) * ( 1 - F_w ) )
        1217             :     !                      / ( ( zeta_w + 2 ) * mixt_frac ); and
        1218             :     !
        1219             :     ! coef_sigma_w_2_sqd = ( 1 - F_w ) / ( ( zeta_w + 2 ) * ( 1 - mixt_frac ) ).
        1220             :     !
        1221             :     ! When w is a responding variable, coef_sigma_w_1_sqd and coef_sigma_w_2_sqd
        1222             :     ! are given by:
        1223             :     !
        1224             :     ! coef_sigma_w_1_sqd = sqrt( mixt_frac * ( 1 - mixt_frac ) ) * Skw
        1225             :     !                      / ( 3 * mixt_frac * sqrt( F_w ) )
        1226             :     !                      - ( 1 + mixt_frac ) * F_w / ( 3 * mixt_frac )
        1227             :     !                      + 1; and
        1228             :     !
        1229             :     ! coef_sigma_w_2_sqd = ( ( 1 - F_w ) - mixt_frac * coef_sigma_w_1_sqd )
        1230             :     !                      / ( 1 - mixt_frac ).
        1231             :     !
        1232             :     ! The equation for <w'4> becomes:
        1233             :     !
        1234             :     ! <w'^4> = ( 3 * mixt_frac * coef_sigma_w_1_sqd^2
        1235             :     !            + 6 * F_w * ( 1 - mixt_frac ) * coef_sigma_w_1_sqd
        1236             :     !            + F_w^2 * ( 1 - mixt_frac )^2 / mixt_frac
        1237             :     !            + 3 * ( 1 - mixt_frac ) * coef_sigma_w_2_sqd^2
        1238             :     !            + 6 * F_w * mixt_frac * coef_sigma_w_2_sqd
        1239             :     !            + F_w^2 * mixt_frac^2 / ( 1 - mixt_frac ) ) * <w'^2>^2.
        1240             :     !
        1241             :     ! This equation is of the form:
        1242             :     !
        1243             :     ! <w'^4> = coef_wp4_implicit * <w'^2>^2;
        1244             :     !
        1245             :     ! where:
        1246             :     !
        1247             :     ! coef_wp4_implicit = 3 * mixt_frac * coef_sigma_w_1_sqd^2
        1248             :     !                     + 6 * F_w * ( 1 - mixt_frac ) * coef_sigma_w_1_sqd
        1249             :     !                     + F_w^2 * ( 1 - mixt_frac )^2 / mixt_frac
        1250             :     !                     + 3 * ( 1 - mixt_frac ) * coef_sigma_w_2_sqd^2
        1251             :     !                     + 6 * F_w * mixt_frac * coef_sigma_w_2_sqd
        1252             :     !                     + F_w^2 * mixt_frac^2 / ( 1 - mixt_frac ).
        1253             :     !
        1254             :     ! While the <w'^4> term is found in the <w'^3> predictive equation and not
        1255             :     ! the <w'^2> predictive equation, the <w'^3> and <w'^2> predictive equations
        1256             :     ! are solved together.  This allows the term containing <w'^4> to be solved
        1257             :     ! implicitly.
        1258             :
        1259             :     ! References:
        1260             :     !-----------------------------------------------------------------------
        1261             :
        1262             :     use constants_clubb, only: &
        1263             :         six,   & ! Variable(s)
        1264             :         three, &
        1265             :         one
        1266             :
        1267             :     use clubb_precision, only: &
        1268             :         core_rknd    ! Procedure(s)
        1269             :
        1270             :     implicit none
        1271             :
        1272             :     integer, intent(in) :: &
        1273             :       nz
        1274             :
        1275             :     ! Input Variables
        1276             :     real ( kind = core_rknd ), dimension(nz), intent(in) :: &
        1277             :       mixt_frac,          & ! Mixture fraction                               [-]
        1278             :       F_w,                & ! Parameter: spread of the PDF comp. means of w  [-]
        1279             :       coef_sigma_w_1_sqd, & ! sigma_w_1^2 = coef_sigma_w_1_sqd * <w'^2>      [-]
        1280             :       coef_sigma_w_2_sqd    ! sigma_w_2^2 = coef_sigma_w_2_sqd * <w'^2>      [-]
        1281             :
        1282             :     ! Return Variable
        1283             :     real ( kind = core_rknd ), dimension(nz) :: &
        1284             :       coef_wp4_implicit    ! Coef.: <w'^4> = coef_wp4_implicit * <w'^2>^2    [-]
        1285             :
        1286             :
        1287             :     ! Calculate coef_wp4_implicit.
        1288             :     coef_wp4_implicit = three * mixt_frac * coef_sigma_w_1_sqd**2 &
        1289             :                         + six * F_w * ( one - mixt_frac ) * coef_sigma_w_1_sqd &
        1290             :                         + F_w**2 * ( one - mixt_frac )**2 / mixt_frac &
        1291             :                         + three * ( one - mixt_frac ) * coef_sigma_w_2_sqd**2 &
        1292             :                         + six * F_w * mixt_frac * coef_sigma_w_2_sqd &
        1293           0 :                         + F_w**2 * mixt_frac**2 / ( one - mixt_frac )
        1294             :
        1295             :
        1296           0 :     return
        1297             :
        1298           0 :   end function calc_coef_wp4_implicit
        1299             :
        1300             :   !=============================================================================
        1301           0 :   function calc_coef_wpxp2_implicit( nz, wp2, xp2, wpxp, sgn_wpxp, &
        1302           0 :                                      mixt_frac, F_w, F_x, &
        1303           0 :                                      coef_sigma_w_1_sqd, &
        1304           0 :                                      coef_sigma_w_2_sqd, &
        1305           0 :                                      coef_sigma_x_1_sqd, &
        1306           0 :                                      coef_sigma_x_2_sqd  ) &
        1307           0 :   result( coef_wpxp2_implicit )
        1308             :
        1309             :     ! Description:
        1310             :     ! The predictive equation for <x'^2> contains a turbulent advection term of
        1311             :     ! the form:
        1312             :     !
        1313             :     ! - ( 1 / rho_ds ) * d ( rho_ds * <w'x'^2> ) / dz;
        1314             :     !
        1315             :     ! where z is height, rho_ds is the dry, base-state density, and <w'x'^2> is
        1316             :     ! calculated by integrating over the PDF.  The equation for <w'x'^2> is:
        1317             :     !
        1318             :     ! <w'x'^2>
        1319             :     ! = mixt_frac * ( ( mu_w_1 - <w> ) * ( ( mu_x_1 - <x> )^2 + sigma_x_1^2 )
        1320             :     !                 + 2 * corr_w_x_1 * sigma_w_1 * sigma_x_1
        1321             :     !                   * ( mu_x_1 - <x> ) )
        1322             :     !   + ( 1 - mixt_frac ) * ( ( mu_w_2 - <w> )
        1323             :     !                           * ( ( mu_x_2 - <x> )^2 + sigma_x_2^2 )
        1324             :     !                           + 2 * corr_w_x_2 * sigma_w_2 * sigma_x_2
        1325             :     !                             * ( mu_x_2 - <x> ) ).
        1326             :     !
        1327             :     ! The following substitutions are made into the above equation:
        1328             :     !
        1329             :     ! mu_w_1 - <w> = sqrt(F_w) * sqrt( ( 1 - mixt_frac ) / mixt_frac )
        1330             :     !                * sqrt( <w'^2> );
        1331             :     !
        1332             :     ! mu_w_2 - <w> = - sqrt(F_w) * sqrt( mixt_frac / ( 1 - mixt_frac ) )
        1333             :     !                  * sqrt( <w'^2> );
        1334             :     !
        1335             :     ! mu_x_1 - <x> = sqrt(F_x) * sqrt( ( 1 - mixt_frac ) / mixt_frac )
        1336             :     !                * sqrt( <x'^2> ) * sgn( <w'x'> );
        1337             :     !
        1338             :     ! mu_x_2 - <x> = - sqrt(F_x) * sqrt( mixt_frac / ( 1 - mixt_frac ) )
        1339             :     !                  * sqrt( <x'^2> ) * sgn( <w'x'> );
        1340             :     !
        1341             :     ! sigma_w_1 = sqrt( coef_sigma_w_1_sqd * <w'^2> );
        1342             :     !
        1343             :     ! sigma_w_2 = sqrt( coef_sigma_w_2_sqd * <w'^2> );
        1344             :     !
        1345             :     ! sigma_x_1 = sqrt( coef_sigma_x_1_sqd * <x'^2> ); and
        1346             :     !
        1347             :     ! sigma_x_2 = sqrt( coef_sigma_x_2_sqd * <x'^2> ).
        1348             :     !
        1349             :     ! Either w can be the setting variable and x can be a responding variable,
        1350             :     ! x can be the setting variable and w can be a responding variable, or both
        1351             :     ! w and x can be responding variables.
        1352             :     !
        1353             :     ! When w is the setting variable, coef_sigma_w_1_sqd and coef_sigma_w_2_sqd
        1354             :     ! are given by:
        1355             :     !
        1356             :     ! coef_sigma_w_1_sqd = ( ( zeta_w + 1 ) * ( 1 - F_w ) )
        1357             :     !                      / ( ( zeta_w + 2 ) * mixt_frac ); and
        1358             :     !
        1359             :     ! coef_sigma_w_2_sqd = ( 1 - F_w ) / ( ( zeta_w + 2 ) * ( 1 - mixt_frac ) ).
        1360             :     !
        1361             :     ! When w is a responding variable, coef_sigma_w_1_sqd and coef_sigma_w_2_sqd
        1362             :     ! are given by:
        1363             :     !
        1364             :     ! coef_sigma_w_1_sqd = sqrt( mixt_frac * ( 1 - mixt_frac ) ) * Skw
        1365             :     !                      / ( 3 * mixt_frac * sqrt( F_w ) )
        1366             :     !                      - ( 1 + mixt_frac ) * F_w / ( 3 * mixt_frac )
        1367             :     !                      + 1; and
        1368             :     !
        1369             :     ! coef_sigma_w_2_sqd = ( ( 1 - F_w ) - mixt_frac * coef_sigma_w_1_sqd )
        1370             :     !                      / ( 1 - mixt_frac ).
        1371             :     !
        1372             :     ! When x is the setting variable, coef_sigma_x_1_sqd and coef_sigma_x_2_sqd
        1373             :     ! are given by:
        1374             :     !
        1375             :     ! coef_sigma_x_1_sqd = ( ( zeta_x + 1 ) * ( 1 - F_x ) )
        1376             :     !                      / ( ( zeta_x + 2 ) * mixt_frac ); and
        1377             :     !
        1378             :     ! coef_sigma_x_2_sqd = ( 1 - F_x ) / ( ( zeta_x + 2 ) * ( 1 - mixt_frac ) ).
        1379             :     !
        1380             :     ! When x is a responding variable, coef_sigma_x_1_sqd and coef_sigma_x_2_sqd
        1381             :     ! are given by:
        1382             :     !
        1383             :     ! coef_sigma_x_1_sqd = sqrt( mixt_frac * ( 1 - mixt_frac ) )
        1384             :     !                      * Skx * sgn( <w'x'> )
        1385             :     !                      / ( 3 * mixt_frac * sqrt( F_x ) )
        1386             :     !                      - ( 1 + mixt_frac ) * F_x / ( 3 * mixt_frac )
        1387             :     !                      + 1; and
        1388             :     !
        1389             :     ! coef_sigma_x_2_sqd = ( ( 1 - F_x ) - mixt_frac * coef_sigma_x_1_sqd )
        1390             :     !                      / ( 1 - mixt_frac ).
        1391             :     !
        1392             :     ! Additionally:
        1393             :     !
        1394             :     ! corr_w_x_1 = corr_w_x_2
        1395             :     ! = ( <w'x'> - mixt_frac * ( mu_w_1 - <w> ) * ( mu_x_1 - <x> )
        1396             :     !            - ( 1 - mixt_frac ) * ( mu_w_2 - <w> ) * ( mu_x_2 - <x> ) )
        1397             :     !   / ( mixt_frac * sigma_w_1 * sigma_x_1
        1398             :     !       + ( 1 - mixt_frac ) * sigma_w_2 * sigma_x_2 );
        1399             :     !
        1400             :     ! where -1 <= corr_w_x_1 = corr_w_x_2 <= 1.  This equation can be rewritten
        1401             :     ! as:
        1402             :     !
        1403             :     ! corr_w_x_1 = corr_w_x_2
        1404             :     ! = ( <w'x'>
        1405             :     !     - sqrt( F_w ) * sqrt( F_x ) * sgn( <w'x'> )
        1406             :     !       * sqrt( <w'^2> ) * sqrt( <x'^2 > ) )
        1407             :     !   / ( ( mixt_frac * sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd )
        1408             :     !         + ( 1 - mixt_frac )
        1409             :     !           * sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) )
        1410             :     !       * sqrt( <w'^2> ) * sqrt( <x'^2> ) ).
        1411             :     !
        1412             :     ! The equation for <w'x'^2> becomes:
        1413             :     !
        1414             :     ! <w'x'^2>
        1415             :     ! = sqrt( mixt_frac * ( 1 - mixt_frac ) ) * sqrt( <w'^2> )
        1416             :     !   * ( sqrt( F_w ) * F_x
        1417             :     !       * ( ( 1 - mixt_frac ) / mixt_frac - mixt_frac / ( 1 - mixt_frac ) )
        1418             :     !       + sqrt( F_w ) * ( coef_sigma_x_1_sqd - coef_sigma_x_2_sqd )
        1419             :     !       + ( 2 * sqrt( F_x ) * sgn( <w'x'> ) * <w'x'>
        1420             :     !           * ( sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd )
        1421             :     !               - sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) ) )
        1422             :     !         / ( ( mixt_frac * sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd )
        1423             :     !               + ( 1 - mixt_frac )
        1424             :     !                 * sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) )
        1425             :     !             * sqrt( <w'^2> ) * sqrt( <x'^2> ) )
        1426             :     !       - ( 2 * sqrt( F_w ) * F_x
        1427             :     !           * ( sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd )
        1428             :     !               - sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) ) )
        1429             :     !         / ( mixt_frac * sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd )
        1430             :     !             + ( 1 - mixt_frac )
        1431             :     !               * sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) ) )
        1432             :     !   * <x'^2>
        1433             :     !
        1434             :     ! This equation is of the form:
        1435             :     !
        1436             :     ! <w'x'^2> = coef_wpxp2_implicit * <x'^2>;
        1437             :     !
        1438             :     ! where:
        1439             :     !
        1440             :     ! coef_wpxp2_implicit
        1441             :     ! = sqrt( mixt_frac * ( 1 - mixt_frac ) ) * sqrt( <w'^2> )
        1442             :     !   * ( sqrt( F_w ) * F_x
        1443             :     !       * ( ( 1 - mixt_frac ) / mixt_frac - mixt_frac / ( 1 - mixt_frac ) )
        1444             :     !       + sqrt( F_w ) * ( coef_sigma_x_1_sqd - coef_sigma_x_2_sqd )
        1445             :     !       + ( 2 * sqrt( F_x ) * sgn( <w'x'> ) * <w'x'>
        1446             :     !           * ( sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd )
        1447             :     !               - sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) ) )
        1448             :     !         / ( ( mixt_frac * sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd )
        1449             :     !               + ( 1 - mixt_frac )
        1450             :     !                 * sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) )
        1451             :     !             * sqrt( <w'^2> ) * sqrt( <x'^2> ) )
        1452             :     !       - ( 2 * sqrt( F_w ) * F_x
        1453             :     !           * ( sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd )
        1454             :     !               - sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) ) )
        1455             :     !         / ( mixt_frac * sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd )
        1456             :     !             + ( 1 - mixt_frac )
        1457             :     !               * sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) ) ).
        1458             :     !
        1459             :     ! In the special case that coef_sigma_w_1_sqd * coef_sigma_x_1_sqd = 0 and
        1460             :     ! coef_sigma_w_2_sqd * coef_sigma_x_2_sqd = 0, the above equation is
        1461             :     ! undefined.  However, the equation for this special case can be derived by
        1462             :     ! taking the original equation for <w'x'^2> and setting both
        1463             :     ! sigma_w_1 * sigma_x_1 = 0 and sigma_w_2 * sigma_x_2 = 0.  The equation
        1464             :     ! becomes:
        1465             :     !
        1466             :     ! <w'x'^2>
        1467             :     ! = mixt_frac * ( ( mu_w_1 - <w> ) * ( ( mu_x_1 - <x> )^2 + sigma_x_1^2 )
        1468             :     !   + ( 1 - mixt_frac ) * ( ( mu_w_2 - <w> )
        1469             :     !                           * ( ( mu_x_2 - <x> )^2 + sigma_x_2^2 );
        1470             :     !
        1471             :     ! and making the same substitutions as before, it can be rewritten as:
        1472             :     !
        1473             :     ! <w'x'^2>
        1474             :     ! = sqrt( mixt_frac * ( 1 - mixt_frac ) ) * sqrt( <w'^2> ) * sqrt( F_w )
        1475             :     !   * ( F_x * ( ( 1 - mixt_frac ) / mixt_frac
        1476             :     !               - mixt_frac / ( 1 - mixt_frac ) )
        1477             :     !       + ( coef_sigma_x_1_sqd - coef_sigma_x_2_sqd ) ) * <x'^2>.
        1478             :     !
        1479             :     ! The coefficient in this special case is:
        1480             :     !
        1481             :     ! coef_wpxp2_implicit
        1482             :     ! = sqrt( mixt_frac * ( 1 - mixt_frac ) ) * sqrt( <w'^2> ) * sqrt( F_w )
        1483             :     !   * ( F_x * ( ( 1 - mixt_frac ) / mixt_frac
        1484             :     !               - mixt_frac / ( 1 - mixt_frac ) )
        1485             :     !       + ( coef_sigma_x_1_sqd - coef_sigma_x_2_sqd ) ).
        1486             :
        1487             :     ! References:
        1488             :     !-----------------------------------------------------------------------
        1489             :
        1490             :     use constants_clubb, only: &
        1491             :         two,  & ! Variable(s)
        1492             :         one,  &
        1493             :         zero
        1494             :
        1495             :     use clubb_precision, only: &
        1496             :         core_rknd    ! Procedure(s)
        1497             :
        1498             :     implicit none
        1499             :
        1500             :     integer, intent(in) :: &
        1501             :       nz
        1502             :
        1503             :     ! Input Variables
        1504             :     real ( kind = core_rknd ), dimension(nz), intent(in) :: &
        1505             :       wp2,                & ! Variance of w (overall)                  [m^2/s^2]
        1506             :       xp2,                & ! Variance of x (overall)           [(units vary)^2]
        1507             :       wpxp,               & ! Covariance of w and x           [m/s (units vary)]
        1508             :       sgn_wpxp,           & ! Sign of the covariance of w and x              [-]
        1509             :       mixt_frac,          & ! Mixture fraction                               [-]
        1510             :       F_w,                & ! Parameter: spread of the PDF comp. means of w  [-]
        1511             :       F_x,                & ! Parameter: spread of the PDF comp. means of x  [-]
        1512             :       coef_sigma_w_1_sqd, & ! sigma_w_1^2 = coef_sigma_w_1_sqd * <w'^2>      [-]
        1513             :       coef_sigma_w_2_sqd, & ! sigma_w_2^2 = coef_sigma_w_2_sqd * <w'^2>      [-]
        1514             :       coef_sigma_x_1_sqd, & ! sigma_x_1^2 = coef_sigma_x_1_sqd * <x'^2>      [-]
        1515             :       coef_sigma_x_2_sqd    ! sigma_x_2^2 = coef_sigma_x_2_sqd * <x'^2>      [-]
        1516             :
        1517             :     ! Return Variable
        1518             :     real ( kind = core_rknd ), dimension(nz) :: &
        1519             :       coef_wpxp2_implicit ! Coef.: <w'x'^2> = coef_wpxp2_implicit * <x'^2> [m/s]
        1520             :
        1521             :     ! Local Variable
        1522             :     real ( kind = core_rknd ), dimension(nz) :: &
        1523           0 :       coefs_factor    ! Factor involving coef_sigma_... coefficients         [-]
        1524             :
        1525             :
        1526             :     ! Calculate coef_wpxp2_implicit.
        1527             :     where ( ( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd > zero &
        1528             :               .or. coef_sigma_w_2_sqd * coef_sigma_x_2_sqd > zero ) &
        1529           0 :             .and. ( wp2 * xp2 > zero ) )
        1530             :
        1531             :        coefs_factor &
        1532             :        = ( sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd ) &
        1533             :            - sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) ) &
        1534             :          / ( mixt_frac * sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd ) &
        1535             :              + ( one - mixt_frac ) &
        1536             :                * sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) )
        1537             :
        1538             :        coef_wpxp2_implicit &
        1539             :        = sqrt( mixt_frac * ( one - mixt_frac ) ) * sqrt( wp2 ) &
        1540             :          * ( sqrt( F_w ) * F_x &
        1541             :              * ( ( one - mixt_frac ) / mixt_frac &
        1542             :                  - mixt_frac / ( one - mixt_frac ) ) &
        1543             :              + sqrt( F_w ) * ( coef_sigma_x_1_sqd - coef_sigma_x_2_sqd ) &
        1544             :              + two * sqrt( F_x ) * coefs_factor * sgn_wpxp * wpxp &
        1545             :                / ( sqrt( wp2 ) * sqrt( xp2 ) ) &
        1546             :              - two * sqrt( F_w ) * F_x * coefs_factor )
        1547             :
        1548             :     elsewhere ! ( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd = 0
        1549             :               !   and coef_sigma_w_2_sqd * coef_sigma_x_2_sqd = 0 )
        1550             :               ! or wp2 * xp2 = 0
        1551             :
        1552             :        coef_wpxp2_implicit &
        1553             :        = sqrt( mixt_frac * ( one - mixt_frac ) ) * sqrt( wp2 ) * sqrt( F_w ) &
        1554             :          * ( F_x * ( ( one - mixt_frac ) / mixt_frac &
        1555             :                      - mixt_frac / ( one - mixt_frac ) ) &
        1556             :              + ( coef_sigma_x_1_sqd - coef_sigma_x_2_sqd ) )
        1557             :
        1558             :     endwhere
        1559             :
        1560             :
        1561           0 :     return
        1562             :
        1563           0 :   end function calc_coef_wpxp2_implicit
        1564             :
        1565             :   !=============================================================================
        1566           0 :   subroutine calc_coefs_wp2xp_semiimpl( nz, wp2, xp2, sgn_wpxp,  & ! In
        1567           0 :                                         mixt_frac, F_w, F_x, & ! In
        1568           0 :                                         coef_sigma_w_1_sqd,  & ! In
        1569           0 :                                         coef_sigma_w_2_sqd,  & ! In
        1570           0 :                                         coef_sigma_x_1_sqd,  & ! In
        1571           0 :                                         coef_sigma_x_2_sqd,  & ! In
        1572           0 :                                         coef_wp2xp_implicit, & ! Out
        1573           0 :                                         term_wp2xp_explicit  ) ! Out
        1574             :
        1575             :     ! Description:
        1576             :     ! The predictive equation for <w'x'> contains a turbulent advection term of
        1577             :     ! the form:
        1578             :     !
        1579             :     ! - ( 1 / rho_ds ) * d ( rho_ds * <w'^2 x'> ) / dz;
        1580             :     !
        1581             :     ! where z is height, rho_ds is the dry, base-state density, and <w'^2 x'> is
        1582             :     ! calculated by integrating over the PDF.  The equation for <w'^2 x'> is:
        1583             :     !
        1584             :     ! <w'^2 x'>
        1585             :     ! = mixt_frac * ( ( mu_x_1 - <x> ) * ( ( mu_w_1 - <w> )^2 + sigma_w_1^2 )
        1586             :     !                 + 2 * corr_w_x_1 * sigma_w_1 * sigma_x_1
        1587             :     !                   * ( mu_w_1 - <w> ) )
        1588             :     !   + ( 1 - mixt_frac ) * ( ( mu_x_2 - <x> )
        1589             :     !                           * ( ( mu_w_2 - <w> )^2 + sigma_w_2^2 )
        1590             :     !                           + 2 * corr_w_x_2 * sigma_w_2 * sigma_x_2
        1591             :     !                             * ( mu_w_2 - <w> ) ).
        1592             :     !
        1593             :     ! The following substitutions are made into the above equation:
        1594             :     !
        1595             :     ! mu_w_1 - <w> = sqrt(F_w) * sqrt( ( 1 - mixt_frac ) / mixt_frac )
        1596             :     !                * sqrt( <w'^2> );
        1597             :     !
        1598             :     ! mu_w_2 - <w> = - sqrt(F_w) * sqrt( mixt_frac / ( 1 - mixt_frac ) )
        1599             :     !                  * sqrt( <w'^2> );
        1600             :     !
        1601             :     ! mu_x_1 - <x> = sqrt(F_x) * sqrt( ( 1 - mixt_frac ) / mixt_frac )
        1602             :     !                * sqrt( <x'^2> ) * sgn( <w'x'> );
        1603             :     !
        1604             :     ! mu_x_2 - <x> = - sqrt(F_x) * sqrt( mixt_frac / ( 1 - mixt_frac ) )
        1605             :     !                  * sqrt( <x'^2> ) * sgn( <w'x'> );
        1606             :     !
        1607             :     ! sigma_w_1 = sqrt( coef_sigma_w_1_sqd * <w'^2> );
        1608             :     !
        1609             :     ! sigma_w_2 = sqrt( coef_sigma_w_2_sqd * <w'^2> );
        1610             :     !
        1611             :     ! sigma_x_1 = sqrt( coef_sigma_x_1_sqd * <x'^2> ); and
        1612             :     !
        1613             :     ! sigma_x_2 = sqrt( coef_sigma_x_2_sqd * <x'^2> ).
        1614             :     !
        1615             :     ! Either w can be the setting variable and x can be a responding variable,
        1616             :     ! x can be the setting variable and w can be a responding variable, or both
        1617             :     ! w and x can be responding variables.
        1618             :     !
        1619             :     ! When w is the setting variable, coef_sigma_w_1_sqd and coef_sigma_w_2_sqd
        1620             :     ! are given by:
        1621             :     !
        1622             :     ! coef_sigma_w_1_sqd = ( ( zeta_w + 1 ) * ( 1 - F_w ) )
        1623             :     !                      / ( ( zeta_w + 2 ) * mixt_frac ); and
        1624             :     !
        1625             :     ! coef_sigma_w_2_sqd = ( 1 - F_w ) / ( ( zeta_w + 2 ) * ( 1 - mixt_frac ) ).
        1626             :     !
        1627             :     ! When w is a responding variable, coef_sigma_w_1_sqd and coef_sigma_w_2_sqd
        1628             :     ! are given by:
        1629             :     !
        1630             :     ! coef_sigma_w_1_sqd = sqrt( mixt_frac * ( 1 - mixt_frac ) ) * Skw
        1631             :     !                      / ( 3 * mixt_frac * sqrt( F_w ) )
        1632             :     !                      - ( 1 + mixt_frac ) * F_w / ( 3 * mixt_frac )
        1633             :     !                      + 1; and
        1634             :     !
        1635             :     ! coef_sigma_w_2_sqd = ( ( 1 - F_w ) - mixt_frac * coef_sigma_w_1_sqd )
        1636             :     !                      / ( 1 - mixt_frac ).
        1637             :     !
        1638             :     ! When x is the setting variable, coef_sigma_x_1_sqd and coef_sigma_x_2_sqd
        1639             :     ! are given by:
        1640             :     !
        1641             :     ! coef_sigma_x_1_sqd = ( ( zeta_x + 1 ) * ( 1 - F_x ) )
        1642             :     !                      / ( ( zeta_x + 2 ) * mixt_frac ); and
        1643             :     !
        1644             :     ! coef_sigma_x_2_sqd = ( 1 - F_x ) / ( ( zeta_x + 2 ) * ( 1 - mixt_frac ) ).
        1645             :     !
        1646             :     ! When x is a responding variable, coef_sigma_x_1_sqd and coef_sigma_x_2_sqd
        1647             :     ! are given by:
        1648             :     !
        1649             :     ! coef_sigma_x_1_sqd = sqrt( mixt_frac * ( 1 - mixt_frac ) )
        1650             :     !                      * Skx * sgn( <w'x'> )
        1651             :     !                      / ( 3 * mixt_frac * sqrt( F_x ) )
        1652             :     !                      - ( 1 + mixt_frac ) * F_x / ( 3 * mixt_frac )
        1653             :     !                      + 1; and
        1654             :     !
        1655             :     ! coef_sigma_x_2_sqd = ( ( 1 - F_x ) - mixt_frac * coef_sigma_x_1_sqd )
        1656             :     !                      / ( 1 - mixt_frac ).
        1657             :     !
        1658             :     ! Additionally:
        1659             :     !
        1660             :     ! corr_w_x_1 = corr_w_x_2
        1661             :     ! = ( <w'x'> - mixt_frac * ( mu_w_1 - <w> ) * ( mu_x_1 - <x> )
        1662             :     !            - ( 1 - mixt_frac ) * ( mu_w_2 - <w> ) * ( mu_x_2 - <x> ) )
        1663             :     !   / ( mixt_frac * sigma_w_1 * sigma_x_1
        1664             :     !       + ( 1 - mixt_frac ) * sigma_w_2 * sigma_x_2 );
        1665             :     !
        1666             :     ! where -1 <= corr_w_x_1 = corr_w_x_2 <= 1.  This equation can be rewritten
        1667             :     ! as:
        1668             :     !
        1669             :     ! corr_w_x_1 = corr_w_x_2
        1670             :     ! = ( <w'x'>
        1671             :     !     - sqrt( F_w ) * sqrt( F_x ) * sgn( <w'x'> )
        1672             :     !       * sqrt( <w'^2> ) * sqrt( <x'^2 > ) )
        1673             :     !   / ( ( mixt_frac * sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd )
        1674             :     !         + ( 1 - mixt_frac )
        1675             :     !           * sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) )
        1676             :     !       * sqrt( <w'^2> ) * sqrt( <x'^2> ) ).
        1677             :     !
        1678             :     ! The equation for <w'^2 x'> becomes:
        1679             :     !
        1680             :     ! <w'^2 x'>
        1681             :     ! = sqrt( mixt_frac * ( 1 - mixt_frac ) ) * sqrt( F_x )
        1682             :     !   * sqrt( <x'^2> ) * <w'^2> * sgn( <w'x'> )
        1683             :     !   * ( F_w * ( ( 1 - mixt_frac ) / mixt_frac
        1684             :     !               - mixt_frac / ( 1 - mixt_frac ) )
        1685             :     !       + ( coef_sigma_w_1_sqd - coef_sigma_w_2_sqd )
        1686             :     !       - ( 2 * F_w * ( sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd )
        1687             :     !                       - sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) )
        1688             :     !         / ( mixt_frac * sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd )
        1689             :     !             + ( 1 - mixt_frac )
        1690             :     !               * sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) ) ) )
        1691             :     !   + ( sqrt( mixt_frac * ( 1 - mixt_frac ) )
        1692             :     !       * 2 * sqrt( F_w ) * sqrt( <w'^2> )
        1693             :     !       * ( sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd )
        1694             :     !           - sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) )
        1695             :     !         / ( mixt_frac * sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd )
        1696             :     !             + ( 1 - mixt_frac )
        1697             :     !               * sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) ) )
        1698             :     !     * <w'x'>
        1699             :     !
        1700             :     ! This equation is of the form:
        1701             :     !
        1702             :     ! <w'^2 x'> = coef_wp2xp_implicit * <w'x'> + term_wp2xp_explicit;
        1703             :     !
        1704             :     ! where:
        1705             :     !
        1706             :     ! coef_wp2xp_implicit
        1707             :     ! = sqrt( mixt_frac * ( 1 - mixt_frac ) ) * 2 * sqrt( F_w ) * sqrt( <w'^2> )
        1708             :     !   * ( sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd )
        1709             :     !       - sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) )
        1710             :     !     / ( mixt_frac * sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd )
        1711             :     !         + ( 1 - mixt_frac )
        1712             :     !           * sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) ); and
        1713             :     !
        1714             :     ! term_wp2xp_explicit
        1715             :     ! = sqrt( mixt_frac * ( 1 - mixt_frac ) ) * sqrt( F_x )
        1716             :     !   * sqrt( <x'^2> ) * <w'^2> * sgn( <w'x'> )
        1717             :     !   * ( F_w * ( ( 1 - mixt_frac ) / mixt_frac
        1718             :     !               - mixt_frac / ( 1 - mixt_frac ) )
        1719             :     !       + ( coef_sigma_w_1_sqd - coef_sigma_w_2_sqd )
        1720             :     !       - ( 2 * F_w * ( sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd )
        1721             :     !                       - sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) )
        1722             :     !         / ( mixt_frac * sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd )
        1723             :     !             + ( 1 - mixt_frac )
        1724             :     !               * sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) ) ) ).
        1725             :     !
        1726             :     ! In the special case that coef_sigma_w_1_sqd * coef_sigma_x_1_sqd = 0 and
        1727             :     ! coef_sigma_w_2_sqd * coef_sigma_x_2_sqd = 0, the above equation is
        1728             :     ! undefined.  However, the equation for this special case can be derived by
        1729             :     ! taking the original equation for <w'^2 x'> and setting both
        1730             :     ! sigma_w_1 * sigma_x_1 = 0 and sigma_w_2 * sigma_x_2 = 0.  The equation
        1731             :     ! becomes:
        1732             :     !
        1733             :     ! <w'^2 x'>
        1734             :     ! = mixt_frac * ( ( mu_x_1 - <x> ) * ( ( mu_w_1 - <w> )^2 + sigma_w_1^2 )
        1735             :     !   + ( 1 - mixt_frac ) * ( ( mu_x_2 - <x> )
        1736             :     !                           * ( ( mu_w_2 - <w> )^2 + sigma_w_2^2 );
        1737             :     !
        1738             :     ! and making the same substitutions as before, it can be rewritten as:
        1739             :     !
        1740             :     ! <w'^2 x'>
        1741             :     ! = sqrt( mixt_frac * ( 1 - mixt_frac ) )
        1742             :     !   * sqrt( F_x ) * sqrt( <x'^2> ) * <w'^2> * sgn( <w'x'> )
        1743             :     !   * ( F_w * ( ( 1 - mixt_frac ) / mixt_frac
        1744             :     !               - mixt_frac / ( 1 - mixt_frac ) )
        1745             :     !       + ( coef_sigma_w_1_sqd - coef_sigma_w_2_sqd ) ).
        1746             :     !
        1747             :     ! Likewise, the equation for <w'x'> in this special case becomes:
        1748             :     !
        1749             :     ! <w'x'> = mixt_frac * ( mu_w_1 - <w> ) * ( mu_x_1 - <x> )
        1750             :     !          + ( 1 - mixt_frac ) * ( mu_w_2 - <w> ) * ( mu_x_2 - <x> );
        1751             :     !
        1752             :     ! and making the same substitutions as before, it can be rewritten as:
        1753             :     !
        1754             :     ! <w'x'> = sqrt( F_w ) * sqrt( F_x )
        1755             :     !          * sqrt( <w'^2> ) * sqrt( <x'^2> ) * sgn( <w'x'> ).
        1756             :     !
        1757             :     ! The equation for <w'^2 x'> can be rewritten as:
        1758             :     !
        1759             :     ! <w'^2 x'>
        1760             :     ! = sqrt( mixt_frac * ( 1 - mixt_frac ) ) * sqrt( F_w ) * sqrt( <w'^2> )
        1761             :     !   * ( ( 1 - mixt_frac ) / mixt_frac - mixt_frac / ( 1 - mixt_frac ) )
        1762             :     !   * <w'x'>
        1763             :     !   + sqrt( mixt_frac * ( 1 - mixt_frac ) )
        1764             :     !     * sqrt( F_x ) * sqrt( <x'^2> ) * <w'^2> * sgn( <w'x'> )
        1765             :     !     * ( coef_sigma_w_1_sqd - coef_sigma_w_2_sqd ).
        1766             :     !
        1767             :     ! The coefficients in this special case are:
        1768             :     !
        1769             :     ! coef_wp2xp_implicit
        1770             :     ! = sqrt( mixt_frac * ( 1 - mixt_frac ) ) * sqrt( F_w ) * sqrt( <w'^2> )
        1771             :     !   * ( ( 1 - mixt_frac ) / mixt_frac - mixt_frac / ( 1 - mixt_frac ) ); and
        1772             :     !
        1773             :     ! term_wp2xp_explicit
        1774             :     ! = sqrt( mixt_frac * ( 1 - mixt_frac ) )
        1775             :     !   * sqrt( F_x ) * sqrt( <x'^2> ) * <w'^2> * sgn( <w'x'> )
        1776             :     !   * ( coef_sigma_w_1_sqd - coef_sigma_w_2_sqd ).
        1777             :
        1778             :     ! References:
        1779             :     !-----------------------------------------------------------------------
        1780             :
        1781             :     use constants_clubb, only: &
        1782             :         two,  & ! Variable(s)
        1783             :         one,  &
        1784             :         zero
        1785             :
        1786             :     use clubb_precision, only: &
        1787             :         core_rknd    ! Procedure(s)
        1788             :
        1789             :     implicit none
        1790             :
        1791             :     integer, intent(in) :: &
        1792             :       nz
        1793             :
        1794             :     ! Input Variables
        1795             :     real ( kind = core_rknd ), dimension(nz), intent(in) :: &
        1796             :       wp2,                & ! Variance of w (overall)                  [m^2/s^2]
        1797             :       xp2,                & ! Variance of x (overall)           [(units vary)^2]
        1798             :       sgn_wpxp,           & ! Sign of the covariance of w and x              [-]
        1799             :       mixt_frac,          & ! Mixture fraction                               [-]
        1800             :       F_w,                & ! Parameter: spread of the PDF comp. means of w  [-]
        1801             :       F_x,                & ! Parameter: spread of the PDF comp. means of x  [-]
        1802             :       coef_sigma_w_1_sqd, & ! sigma_w_1^2 = coef_sigma_w_1_sqd * <w'^2>      [-]
        1803             :       coef_sigma_w_2_sqd, & ! sigma_w_2^2 = coef_sigma_w_2_sqd * <w'^2>      [-]
        1804             :       coef_sigma_x_1_sqd, & ! sigma_x_1^2 = coef_sigma_x_1_sqd * <x'^2>      [-]
        1805             :       coef_sigma_x_2_sqd    ! sigma_x_2^2 = coef_sigma_x_2_sqd * <x'^2>      [-]
        1806             :
        1807             :     ! Output Variables
        1808             :     ! Coefs.: <w'^2 x'> = coef_wp2xp_implicit * <w'x'> + term_wp2xp_explicit
        1809             :     real ( kind = core_rknd ), dimension(nz), intent(out) :: &
        1810             :       coef_wp2xp_implicit, & ! Coefficient that is multiplied by <w'x'>    [m/s]
        1811             :       term_wp2xp_explicit    ! Term that is on the RHS    [m^2/s^2 (units vary)]
        1812             :
        1813             :     ! Local Variable
        1814             :     real ( kind = core_rknd ), dimension(nz) :: &
        1815           0 :       coefs_factor    ! Factor involving coef_sigma_... coefficients         [-]
        1816             :
        1817             :
        1818             :     ! Calculate coef_wp2xp_implicit and term_wp2xp_explicit.
        1819             :     where ( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd > zero &
        1820           0 :             .or. coef_sigma_w_2_sqd * coef_sigma_x_2_sqd > zero )
        1821             :
        1822             :        coefs_factor &
        1823             :        = ( sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd ) &
        1824             :            - sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) ) &
        1825             :          / ( mixt_frac * sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd ) &
        1826             :              + ( one - mixt_frac ) &
        1827             :                * sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) )
        1828             :
        1829             :        coef_wp2xp_implicit &
        1830             :        = sqrt( mixt_frac * ( one - mixt_frac ) ) &
        1831             :          * two * sqrt( F_w ) * sqrt( wp2 ) * coefs_factor
        1832             :
        1833             :        term_wp2xp_explicit &
        1834             :        = sqrt( mixt_frac * ( one - mixt_frac ) ) &
        1835             :          * sqrt( F_x ) * sqrt( xp2 ) * wp2 * sgn_wpxp &
        1836             :          * ( F_w * ( ( one - mixt_frac ) / mixt_frac &
        1837             :                      - mixt_frac / ( one - mixt_frac ) ) &
        1838             :              + ( coef_sigma_w_1_sqd - coef_sigma_w_2_sqd ) &
        1839             :              - two * F_w * coefs_factor )
        1840             :
        1841             :     elsewhere ! coef_sigma_w_1_sqd * coef_sigma_x_1_sqd = 0
        1842             :               ! and coef_sigma_w_2_sqd * coef_sigma_x_2_sqd = 0
        1843             :
        1844             :        coef_wp2xp_implicit &
        1845             :        = sqrt( mixt_frac * ( one - mixt_frac ) ) &
        1846             :          * sqrt( F_w ) * sqrt( wp2 ) &
        1847             :          * ( ( one - mixt_frac ) / mixt_frac &
        1848             :              - mixt_frac / ( one - mixt_frac ) )
        1849             :
        1850             :        term_wp2xp_explicit &
        1851             :        = sqrt( mixt_frac * ( one - mixt_frac ) ) &
        1852             :          * sqrt( F_x ) * sqrt( xp2 ) * wp2 * sgn_wpxp &
        1853             :          * ( coef_sigma_w_1_sqd - coef_sigma_w_2_sqd )
        1854             :
        1855             :     endwhere
        1856             :
        1857             :
        1858           0 :     return
        1859             :
        1860             :   end subroutine calc_coefs_wp2xp_semiimpl
        1861             :
        1862             :   !=============================================================================
        1863           0 :   subroutine calc_coefs_wpxpyp_semiimpl( nz, wp2, xp2, yp2, wpxp,      & ! In
        1864           0 :                                          wpyp, sgn_wpxp, sgn_wpyp, & ! In
        1865           0 :                                          mixt_frac, F_w, F_x, F_y, & ! In
        1866           0 :                                          coef_sigma_w_1_sqd,       & ! In
        1867           0 :                                          coef_sigma_w_2_sqd,       & ! In
        1868           0 :                                          coef_sigma_x_1_sqd,       & ! In
        1869           0 :                                          coef_sigma_x_2_sqd,       & ! In
        1870           0 :                                          coef_sigma_y_1_sqd,       & ! In
        1871           0 :                                          coef_sigma_y_2_sqd,       & ! In
        1872           0 :                                          coef_wpxpyp_implicit,     & ! Out
        1873           0 :                                          term_wpxpyp_explicit      ) ! Out
        1874             :
        1875             :     ! Description:
        1876             :     ! The predictive equation for <w'x'y'> contains a turbulent advection term
        1877             :     ! of the form:
        1878             :     !
        1879             :     ! - ( 1 / rho_ds ) * d ( rho_ds * <w'x'y'> ) / dz;
        1880             :     !
        1881             :     ! where z is height, rho_ds is the dry, base-state density, and <w'x'y'> is
        1882             :     ! calculated by integrating over the PDF.  The equation for <w'x'y'> is:
        1883             :     !
        1884             :     ! <w'x'y'>
        1885             :     ! = mixt_frac
        1886             :     !   * ( ( mu_w_1 - <w> ) * ( mu_x_1 - <x> ) * ( mu_y_1 - <y> )
        1887             :     !       + corr_x_y_1 * sigma_x_1 * sigma_y_1 * ( mu_w_1 - <w> )
        1888             :     !       + corr_w_y_1 * sigma_w_1 * sigma_y_1 * ( mu_x_1 - <x> )
        1889             :     !       + corr_w_x_1 * sigma_w_1 * sigma_x_1 * ( mu_y_1 - <y> ) )
        1890             :     !   + ( 1 - mixt_frac )
        1891             :     !     * ( ( mu_w_2 - <w> ) * ( mu_x_2 - <x> ) * ( mu_y_2 - <y> )
        1892             :     !         + corr_x_y_2 * sigma_x_2 * sigma_y_2 * ( mu_w_2 - <w> )
        1893             :     !         + corr_w_y_2 * sigma_w_2 * sigma_y_2 * ( mu_x_2 - <x> )
        1894             :     !         + corr_w_x_2 * sigma_w_2 * sigma_x_2 * ( mu_y_2 - <y> ) ).
        1895             :     !
        1896             :     ! The following substitutions are made into the above equation:
        1897             :     !
        1898             :     ! mu_w_1 - <w> = sqrt(F_w) * sqrt( ( 1 - mixt_frac ) / mixt_frac )
        1899             :     !                * sqrt( <w'^2> );
        1900             :     !
        1901             :     ! mu_w_2 - <w> = - sqrt(F_w) * sqrt( mixt_frac / ( 1 - mixt_frac ) )
        1902             :     !                  * sqrt( <w'^2> );
        1903             :     !
        1904             :     ! mu_x_1 - <x> = sqrt(F_x) * sqrt( ( 1 - mixt_frac ) / mixt_frac )
        1905             :     !                * sqrt( <x'^2> ) * sgn( <w'x'> );
        1906             :     !
        1907             :     ! mu_x_2 - <x> = - sqrt(F_x) * sqrt( mixt_frac / ( 1 - mixt_frac ) )
        1908             :     !                  * sqrt( <x'^2> ) * sgn( <w'x'> );
        1909             :     !
        1910             :     ! mu_y_1 - <y> = sqrt(F_y) * sqrt( ( 1 - mixt_frac ) / mixt_frac )
        1911             :     !                * sqrt( <y'^2> ) * sgn( <w'y'> );
        1912             :     !
        1913             :     ! mu_y_2 - <y> = - sqrt(F_y) * sqrt( mixt_frac / ( 1 - mixt_frac ) )
        1914             :     !                  * sqrt( <y'^2> ) * sgn( <w'y'> );
        1915             :     !
        1916             :     ! sigma_w_1 = sqrt( coef_sigma_w_1_sqd * <w'^2> );
        1917             :     !
        1918             :     ! sigma_w_2 = sqrt( coef_sigma_w_2_sqd * <w'^2> );
        1919             :     !
        1920             :     ! sigma_x_1 = sqrt( coef_sigma_x_1_sqd * <x'^2> );
        1921             :     !
        1922             :     ! sigma_x_2 = sqrt( coef_sigma_x_2_sqd * <x'^2> );
        1923             :     !
        1924             :     ! sigma_y_1 = sqrt( coef_sigma_y_1_sqd * <y'^2> ); and
        1925             :     !
        1926             :     ! sigma_y_2 = sqrt( coef_sigma_y_2_sqd * <y'^2> ).
        1927             :     !
        1928             :     ! Either w can be the setting variable and both x and y can be responding
        1929             :     ! variables, x can be the setting variable and both w and y can be
        1930             :     ! responding variables, y can be the setting variable and both w and x can
        1931             :     ! be responding variables, or all of w, x, and y can be responding
        1932             :     ! variables.
        1933             :     !
        1934             :     ! When w is the setting variable, coef_sigma_w_1_sqd and coef_sigma_w_2_sqd
        1935             :     ! are given by:
        1936             :     !
        1937             :     ! coef_sigma_w_1_sqd = ( ( zeta_w + 1 ) * ( 1 - F_w ) )
        1938             :     !                      / ( ( zeta_w + 2 ) * mixt_frac ); and
        1939             :     !
        1940             :     ! coef_sigma_w_2_sqd = ( 1 - F_w ) / ( ( zeta_w + 2 ) * ( 1 - mixt_frac ) ).
        1941             :     !
        1942             :     ! When w is a responding variable, coef_sigma_w_1_sqd and coef_sigma_w_2_sqd
        1943             :     ! are given by:
        1944             :     !
        1945             :     ! coef_sigma_w_1_sqd = sqrt( mixt_frac * ( 1 - mixt_frac ) ) * Skw
        1946             :     !                      / ( 3 * mixt_frac * sqrt( F_w ) )
        1947             :     !                      - ( 1 + mixt_frac ) * F_w / ( 3 * mixt_frac )
        1948             :     !                      + 1; and
        1949             :     !
        1950             :     ! coef_sigma_w_2_sqd = ( ( 1 - F_w ) - mixt_frac * coef_sigma_w_1_sqd )
        1951             :     !                      / ( 1 - mixt_frac ).
        1952             :     !
        1953             :     ! When x is the setting variable, coef_sigma_x_1_sqd and coef_sigma_x_2_sqd
        1954             :     ! are given by:
        1955             :     !
        1956             :     ! coef_sigma_x_1_sqd = ( ( zeta_x + 1 ) * ( 1 - F_x ) )
        1957             :     !                      / ( ( zeta_x + 2 ) * mixt_frac ); and
        1958             :     !
        1959             :     ! coef_sigma_x_2_sqd = ( 1 - F_x ) / ( ( zeta_x + 2 ) * ( 1 - mixt_frac ) ).
        1960             :     !
        1961             :     ! When x is a responding variable, coef_sigma_x_1_sqd and coef_sigma_x_2_sqd
        1962             :     ! are given by:
        1963             :     !
        1964             :     ! coef_sigma_x_1_sqd = sqrt( mixt_frac * ( 1 - mixt_frac ) )
        1965             :     !                      * Skx * sgn( <w'x'> )
        1966             :     !                      / ( 3 * mixt_frac * sqrt( F_x ) )
        1967             :     !                      - ( 1 + mixt_frac ) * F_x / ( 3 * mixt_frac )
        1968             :     !                      + 1; and
        1969             :     !
        1970             :     ! coef_sigma_x_2_sqd = ( ( 1 - F_x ) - mixt_frac * coef_sigma_x_1_sqd )
        1971             :     !                      / ( 1 - mixt_frac ).
        1972             :     !
        1973             :     ! When y is the setting variable, coef_sigma_y_1_sqd and coef_sigma_y_2_sqd
        1974             :     ! are given by:
        1975             :     !
        1976             :     ! coef_sigma_y_1_sqd = ( ( zeta_y + 1 ) * ( 1 - F_y ) )
        1977             :     !                      / ( ( zeta_y + 2 ) * mixt_frac ); and
        1978             :     !
        1979             :     ! coef_sigma_y_2_sqd = ( 1 - F_y ) / ( ( zeta_y + 2 ) * ( 1 - mixt_frac ) ).
        1980             :     !
        1981             :     ! When y is a responding variable, coef_sigma_y_1_sqd and coef_sigma_y_2_sqd
        1982             :     ! are given by:
        1983             :     !
        1984             :     ! coef_sigma_y_1_sqd = sqrt( mixt_frac * ( 1 - mixt_frac ) )
        1985             :     !                      * Sky * sgn( <w'y'> )
        1986             :     !                      / ( 3 * mixt_frac * sqrt( F_y ) )
        1987             :     !                      - ( 1 + mixt_frac ) * F_y / ( 3 * mixt_frac )
        1988             :     !                      + 1; and
        1989             :     !
        1990             :     ! coef_sigma_y_2_sqd = ( ( 1 - F_y ) - mixt_frac * coef_sigma_y_1_sqd )
        1991             :     !                      / ( 1 - mixt_frac ).
        1992             :     !
        1993             :     ! Additionally:
        1994             :     !
        1995             :     ! corr_w_x_1 = corr_w_x_2
        1996             :     ! = ( <w'x'> - mixt_frac * ( mu_w_1 - <w> ) * ( mu_x_1 - <x> )
        1997             :     !            - ( 1 - mixt_frac ) * ( mu_w_2 - <w> ) * ( mu_x_2 - <x> ) )
        1998             :     !   / ( mixt_frac * sigma_w_1 * sigma_x_1
        1999             :     !       + ( 1 - mixt_frac ) * sigma_w_2 * sigma_x_2 );
        2000             :     !
        2001             :     ! where -1 <= corr_w_x_1 = corr_w_x_2 <= 1.  This equation can be rewritten
        2002             :     ! as:
        2003             :     !
        2004             :     ! corr_w_x_1 = corr_w_x_2
        2005             :     ! = ( <w'x'>
        2006             :     !     - sqrt( F_w ) * sqrt( F_x ) * sgn( <w'x'> )
        2007             :     !       * sqrt( <w'^2> ) * sqrt( <x'^2 > ) )
        2008             :     !   / ( ( mixt_frac * sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd )
        2009             :     !         + ( 1 - mixt_frac )
        2010             :     !           * sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) )
        2011             :     !       * sqrt( <w'^2> ) * sqrt( <x'^2> ) ).
        2012             :     !
        2013             :     ! Likewise:
        2014             :     !
        2015             :     ! corr_w_y_1 = corr_w_y_2
        2016             :     ! = ( <w'y'>
        2017             :     !     - sqrt( F_w ) * sqrt( F_y ) * sgn( <w'y'> )
        2018             :     !       * sqrt( <w'^2> ) * sqrt( <y'^2 > ) )
        2019             :     !   / ( ( mixt_frac * sqrt( coef_sigma_w_1_sqd * coef_sigma_y_1_sqd )
        2020             :     !         + ( 1 - mixt_frac )
        2021             :     !           * sqrt( coef_sigma_w_2_sqd * coef_sigma_y_2_sqd ) )
        2022             :     !       * sqrt( <w'^2> ) * sqrt( <y'^2> ) ); and
        2023             :     !
        2024             :     ! corr_x_y_1 = corr_x_y_2
        2025             :     ! = ( <x'y'>
        2026             :     !     - sqrt( F_x ) * sqrt( F_y ) * sgn( <w'x'> ) * sgn( <w'y'> )
        2027             :     !       * sqrt( <x'^2> ) * sqrt( <y'^2 > ) )
        2028             :     !   / ( ( mixt_frac * sqrt( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd )
        2029             :     !         + ( 1 - mixt_frac )
        2030             :     !           * sqrt( coef_sigma_x_2_sqd * coef_sigma_y_2_sqd ) )
        2031             :     !       * sqrt( <x'^2> ) * sqrt( <y'^2> ) ).
        2032             :     !
        2033             :     ! The equation for <w'x'y'> becomes:
        2034             :     !
        2035             :     ! <w'x'y'>
        2036             :     ! = sqrt( mixt_frac * ( 1 - mixt_frac ) ) * sqrt( F_w ) * sqrt( <w'^2> )
        2037             :     !   * ( sqrt( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd )
        2038             :     !             - sqrt( coef_sigma_x_2_sqd * coef_sigma_y_2_sqd ) )
        2039             :     !     / ( mixt_frac * sqrt( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd )
        2040             :     !         + ( 1 - mixt_frac )
        2041             :     !           * sqrt( coef_sigma_x_2_sqd * coef_sigma_y_2_sqd ) )
        2042             :     !   * <x'y'>
        2043             :     !   + sqrt( mixt_frac * ( 1 - mixt_frac ) ) * sqrt( F_w ) * sqrt( <w'^2> )
        2044             :     !     * sqrt( F_x ) * sqrt( <x'^2> ) * sgn( <w'x'> )
        2045             :     !     * sqrt( F_y ) * sqrt( <y'^2> ) * sgn( <w'y'> )
        2046             :     !     * ( ( 1 - mixt_frac ) / mixt_frac - mixt_frac / ( 1 - mixt_frac )
        2047             :     !         - ( sqrt( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd )
        2048             :     !             - sqrt( coef_sigma_x_2_sqd * coef_sigma_y_2_sqd ) )
        2049             :     !           / ( mixt_frac * sqrt( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd )
        2050             :     !               + ( 1 - mixt_frac )
        2051             :     !                 * sqrt( coef_sigma_x_2_sqd * coef_sigma_y_2_sqd ) )
        2052             :     !         - ( sqrt( coef_sigma_w_1_sqd * coef_sigma_y_1_sqd )
        2053             :     !             - sqrt( coef_sigma_w_2_sqd * coef_sigma_y_2_sqd ) )
        2054             :     !           / ( mixt_frac * sqrt( coef_sigma_w_1_sqd * coef_sigma_y_1_sqd )
        2055             :     !               + ( 1 - mixt_frac )
        2056             :     !                 * sqrt( coef_sigma_w_2_sqd * coef_sigma_y_2_sqd ) )
        2057             :     !         - ( sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd )
        2058             :     !             - sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) )
        2059             :     !           / ( mixt_frac * sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd )
        2060             :     !               + ( 1 - mixt_frac )
        2061             :     !                 * sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) ) )
        2062             :     !   + sqrt( mixt_frac * ( 1 - mixt_frac ) )
        2063             :     !     * sqrt( F_x ) * sqrt( <x'^2> ) * sgn( <w'x'> )
        2064             :     !     * ( sqrt( coef_sigma_w_1_sqd * coef_sigma_y_1_sqd )
        2065             :     !         - sqrt( coef_sigma_w_2_sqd * coef_sigma_y_2_sqd ) )
        2066             :     !       / ( mixt_frac * sqrt( coef_sigma_w_1_sqd * coef_sigma_y_1_sqd )
        2067             :     !           + ( 1 - mixt_frac )
        2068             :     !             * sqrt( coef_sigma_w_2_sqd * coef_sigma_y_2_sqd ) )
        2069             :     !     * <w'y'>
        2070             :     !   + sqrt( mixt_frac * ( 1 - mixt_frac ) )
        2071             :     !     * sqrt( F_y ) * sqrt( <y'^2> ) * sgn( <w'y'> )
        2072             :     !     * ( sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd )
        2073             :     !         - sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) )
        2074             :     !       / ( mixt_frac * sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd )
        2075             :     !           + ( 1 - mixt_frac )
        2076             :     !             * sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) )
        2077             :     !     * <w'x'>.
        2078             :     !
        2079             :     ! This equation is of the form:
        2080             :     !
        2081             :     ! <w'^2 x'> = coef_wp2xp_implicit * <w'x'> + term_wp2xp_explicit;
        2082             :     !
        2083             :     ! where:
        2084             :     !
        2085             :     ! coef_wpxpyp_implicit
        2086             :     ! = sqrt( mixt_frac * ( 1 - mixt_frac ) ) * sqrt( F_w ) * sqrt( <w'^2> )
        2087             :     !   * ( sqrt( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd )
        2088             :     !             - sqrt( coef_sigma_x_2_sqd * coef_sigma_y_2_sqd ) )
        2089             :     !     / ( mixt_frac * sqrt( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd )
        2090             :     !         + ( 1 - mixt_frac )
        2091             :     !           * sqrt( coef_sigma_x_2_sqd * coef_sigma_y_2_sqd ) ); and
        2092             :     !
        2093             :     ! term_wpxpyp_explicit
        2094             :     ! = sqrt( mixt_frac * ( 1 - mixt_frac ) ) * sqrt( F_w ) * sqrt( <w'^2> )
        2095             :     !   * sqrt( F_x ) * sqrt( <x'^2> ) * sgn( <w'x'> )
        2096             :     !   * sqrt( F_y ) * sqrt( <y'^2> ) * sgn( <w'y'> )
        2097             :     !   * ( ( 1 - mixt_frac ) / mixt_frac - mixt_frac / ( 1 - mixt_frac )
        2098             :     !       - ( sqrt( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd )
        2099             :     !           - sqrt( coef_sigma_x_2_sqd * coef_sigma_y_2_sqd ) )
        2100             :     !         / ( mixt_frac * sqrt( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd )
        2101             :     !             + ( 1 - mixt_frac )
        2102             :     !               * sqrt( coef_sigma_x_2_sqd * coef_sigma_y_2_sqd ) )
        2103             :     !       - ( sqrt( coef_sigma_w_1_sqd * coef_sigma_y_1_sqd )
        2104             :     !           - sqrt( coef_sigma_w_2_sqd * coef_sigma_y_2_sqd ) )
        2105             :     !         / ( mixt_frac * sqrt( coef_sigma_w_1_sqd * coef_sigma_y_1_sqd )
        2106             :     !             + ( 1 - mixt_frac )
        2107             :     !               * sqrt( coef_sigma_w_2_sqd * coef_sigma_y_2_sqd ) )
        2108             :     !       - ( sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd )
        2109             :     !           - sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) )
        2110             :     !         / ( mixt_frac * sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd )
        2111             :     !             + ( 1 - mixt_frac )
        2112             :     !               * sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) ) )
        2113             :     !   + sqrt( mixt_frac * ( 1 - mixt_frac ) )
        2114             :     !     * sqrt( F_x ) * sqrt( <x'^2> ) * sgn( <w'x'> )
        2115             :     !     * ( sqrt( coef_sigma_w_1_sqd * coef_sigma_y_1_sqd )
        2116             :     !         - sqrt( coef_sigma_w_2_sqd * coef_sigma_y_2_sqd ) )
        2117             :     !       / ( mixt_frac * sqrt( coef_sigma_w_1_sqd * coef_sigma_y_1_sqd )
        2118             :     !           + ( 1 - mixt_frac )
        2119             :     !             * sqrt( coef_sigma_w_2_sqd * coef_sigma_y_2_sqd ) )
        2120             :     !     * <w'y'>
        2121             :     !   + sqrt( mixt_frac * ( 1 - mixt_frac ) )
        2122             :     !     * sqrt( F_y ) * sqrt( <y'^2> ) * sgn( <w'y'> )
        2123             :     !     * ( sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd )
        2124             :     !         - sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) )
        2125             :     !       / ( mixt_frac * sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd )
        2126             :     !           + ( 1 - mixt_frac )
        2127             :     !             * sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) )
        2128             :     !     * <w'x'>.
        2129             :     !
        2130             :     ! There are also special cases for the above equations.
        2131             :
        2132             :     ! References:
        2133             :     !-----------------------------------------------------------------------
        2134             :
        2135             :     use constants_clubb, only: &
        2136             :         one,  & ! Variable(s)
        2137             :         zero
        2138             :
        2139             :     use clubb_precision, only: &
        2140             :         core_rknd    ! Procedure(s)
        2141             :
        2142             :     implicit none
        2143             :
        2144             :     integer, intent(in) :: &
        2145             :       nz
        2146             :
        2147             :     ! Input Variables
        2148             :     real ( kind = core_rknd ), dimension(nz), intent(in) :: &
        2149             :       wp2,                & ! Variance of w (overall)                  [m^2/s^2]
        2150             :       xp2,                & ! Variance of x (overall)              [(x units)^2]
        2151             :       yp2,                & ! Variance of y (overall)              [(y units)^2]
        2152             :       wpxp,               & ! Covariance of w and x (overall)     [m/s(x units)]
        2153             :       wpyp,               & ! Covariance of w and y (overall)     [m/s(y units)]
        2154             :       sgn_wpxp,           & ! Sign of the covariance of w and x              [-]
        2155             :       sgn_wpyp,           & ! Sign of the covariance of w and y              [-]
        2156             :       mixt_frac,          & ! Mixture fraction                               [-]
        2157             :       F_w,                & ! Parameter: spread of the PDF comp. means of w  [-]
        2158             :       F_x,                & ! Parameter: spread of the PDF comp. means of x  [-]
        2159             :       F_y,                & ! Parameter: spread of the PDF comp. means of y  [-]
        2160             :       coef_sigma_w_1_sqd, & ! sigma_w_1^2 = coef_sigma_w_1_sqd * <w'^2>      [-]
        2161             :       coef_sigma_w_2_sqd, & ! sigma_w_2^2 = coef_sigma_w_2_sqd * <w'^2>      [-]
        2162             :       coef_sigma_x_1_sqd, & ! sigma_x_1^2 = coef_sigma_x_1_sqd * <x'^2>      [-]
        2163             :       coef_sigma_x_2_sqd, & ! sigma_x_2^2 = coef_sigma_x_2_sqd * <x'^2>      [-]
        2164             :       coef_sigma_y_1_sqd, & ! sigma_y_1^2 = coef_sigma_y_1_sqd * <y'^2>      [-]
        2165             :       coef_sigma_y_2_sqd    ! sigma_y_2^2 = coef_sigma_y_2_sqd * <y'^2>      [-]
        2166             :
        2167             :     ! Output Variables
        2168             :     ! Coefs.: <w'x'y'> = coef_wpxpyp_implicit * <x'y'> + term_wpxpyp_explicit
        2169             :     real ( kind = core_rknd ), dimension(nz), intent(out) :: &
        2170             :       coef_wpxpyp_implicit, & ! Coefficient that is multiplied by <x'y'>   [m/s]
        2171             :       term_wpxpyp_explicit    ! Term that is on the RHS  [m/s(x units)(y units)]
        2172             :
        2173             :     ! Local Variables
        2174             :     real ( kind = core_rknd ), dimension(nz) :: &
        2175           0 :       coefs_factor_wx, & ! Factor involving coef_sigma_... w and x coefs     [-]
        2176           0 :       coefs_factor_wy, & ! Factor involving coef_sigma_... w and y coefs     [-]
        2177           0 :       coefs_factor_xy    ! Factor involving coef_sigma_... x and y coefs     [-]
        2178             :
        2179             :
        2180             :     ! Calculate coefs_factor_wx.
        2181             :     where ( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd > zero &
        2182           0 :             .or. coef_sigma_w_2_sqd * coef_sigma_x_2_sqd > zero )
        2183             :
        2184             :        ! coefs_factor_wx
        2185             :        ! = ( sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd )
        2186             :        !     - sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) )
        2187             :        !   / ( mixt_frac * sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd )
        2188             :        !       + ( 1 - mixt_frac )
        2189             :        !         * sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) )
        2190             :        coefs_factor_wx &
        2191             :        = ( sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd ) &
        2192             :            - sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) ) &
        2193             :          / ( mixt_frac * sqrt( coef_sigma_w_1_sqd * coef_sigma_x_1_sqd ) &
        2194             :              + ( one - mixt_frac ) &
        2195             :                * sqrt( coef_sigma_w_2_sqd * coef_sigma_x_2_sqd ) )
        2196             :
        2197             :     elsewhere ! coef_sigma_w_1_sqd * coef_sigma_x_1_sqd = 0
        2198             :               ! and coef_sigma_w_2_sqd * coef_sigma_x_2_sqd = 0
        2199             :
        2200             :        ! When coef_sigma_w_1_sqd * coef_sigma_x_1_sqd = 0 and
        2201             :        ! coef_sigma_w_2_sqd * coef_sigma_x_2_sqd = 0, the value of
        2202             :        ! coefs_factor_wx is undefined.  However, setting coefs_factor_wx to a
        2203             :        ! value of 0 in this scenario allows for the use of general form
        2204             :        ! equations below for coef_wpxpyp_implicit and term_wpxpyp_explicit.
        2205             :        coefs_factor_wx = zero
        2206             :
        2207             :     endwhere
        2208             :
        2209             :     ! Calculate coefs_factor_wy.
        2210             :     where ( coef_sigma_w_1_sqd * coef_sigma_y_1_sqd > zero &
        2211           0 :             .or. coef_sigma_w_2_sqd * coef_sigma_y_2_sqd > zero )
        2212             :
        2213             :        ! coefs_factor_wy
        2214             :        ! = ( sqrt( coef_sigma_w_1_sqd * coef_sigma_y_1_sqd )
        2215             :        !     - sqrt( coef_sigma_w_2_sqd * coef_sigma_y_2_sqd ) )
        2216             :        !   / ( mixt_frac * sqrt( coef_sigma_w_1_sqd * coef_sigma_y_1_sqd )
        2217             :        !       + ( 1 - mixt_frac )
        2218             :        !         * sqrt( coef_sigma_w_2_sqd * coef_sigma_y_2_sqd ) )
        2219             :        coefs_factor_wy &
        2220             :        = ( sqrt( coef_sigma_w_1_sqd * coef_sigma_y_1_sqd ) &
        2221             :            - sqrt( coef_sigma_w_2_sqd * coef_sigma_y_2_sqd ) ) &
        2222             :          / ( mixt_frac * sqrt( coef_sigma_w_1_sqd * coef_sigma_y_1_sqd ) &
        2223             :              + ( one - mixt_frac ) &
        2224             :                * sqrt( coef_sigma_w_2_sqd * coef_sigma_y_2_sqd ) )
        2225             :
        2226             :     elsewhere ! coef_sigma_w_1_sqd * coef_sigma_y_1_sqd = 0
        2227             :               ! and coef_sigma_w_2_sqd * coef_sigma_y_2_sqd = 0
        2228             :
        2229             :        ! When coef_sigma_w_1_sqd * coef_sigma_y_1_sqd = 0 and
        2230             :        ! coef_sigma_w_2_sqd * coef_sigma_y_2_sqd = 0, the value of
        2231             :        ! coefs_factor_wy is undefined.  However, setting coefs_factor_wy to a
        2232             :        ! value of 0 in this scenario allows for the use of general form
        2233             :        ! equations below for coef_wpxpyp_implicit and term_wpxpyp_explicit.
        2234             :        coefs_factor_wy = zero
        2235             :
        2236             :     endwhere
        2237             :
        2238             :     ! Calculate coefs_factor_xy.
        2239             :     where ( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd > zero &
        2240           0 :             .or. coef_sigma_x_2_sqd * coef_sigma_y_2_sqd > zero )
        2241             :
        2242             :        ! coefs_factor_xy
        2243             :        ! = ( sqrt( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd )
        2244             :        !     - sqrt( coef_sigma_x_2_sqd * coef_sigma_y_2_sqd ) )
        2245             :        !   / ( mixt_frac * sqrt( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd )
        2246             :        !       + ( 1 - mixt_frac )
        2247             :        !         * sqrt( coef_sigma_x_2_sqd * coef_sigma_y_2_sqd ) )
        2248             :        coefs_factor_xy &
        2249             :        = ( sqrt( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd ) &
        2250             :            - sqrt( coef_sigma_x_2_sqd * coef_sigma_y_2_sqd ) ) &
        2251             :          / ( mixt_frac * sqrt( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd ) &
        2252             :              + ( one - mixt_frac ) &
        2253             :                * sqrt( coef_sigma_x_2_sqd * coef_sigma_y_2_sqd ) )
        2254             :
        2255             :     elsewhere ! coef_sigma_x_1_sqd * coef_sigma_y_1_sqd = 0
        2256             :               ! and coef_sigma_x_2_sqd * coef_sigma_y_2_sqd = 0
        2257             :
        2258             :        ! When coef_sigma_x_1_sqd * coef_sigma_y_1_sqd = 0 and
        2259             :        ! coef_sigma_x_2_sqd * coef_sigma_y_2_sqd = 0, the value of
        2260             :        ! coefs_factor_xy is undefined.  However, setting coefs_factor_xy to a
        2261             :        ! value of 0 in this scenario allows for the use of general form
        2262             :        ! equations below for coef_wpxpyp_implicit and term_wpxpyp_explicit.
        2263             :        coefs_factor_xy = zero
        2264             :
        2265             :     endwhere
        2266             :
        2267             :
        2268             :     ! Calculate coef_wpxpyp_implicit and term_wpxpyp_explicit.
        2269             :     where ( coef_sigma_x_1_sqd * coef_sigma_y_1_sqd > zero &
        2270           0 :             .or. coef_sigma_x_2_sqd * coef_sigma_y_2_sqd > zero )
        2271             :
        2272             :        coef_wpxpyp_implicit &
        2273             :        = sqrt( mixt_frac * ( one - mixt_frac ) ) &
        2274             :          * sqrt( F_w ) * sqrt( wp2 ) * coefs_factor_xy
        2275             :
        2276             :        term_wpxpyp_explicit &
        2277             :        = sqrt( mixt_frac * ( one - mixt_frac ) ) * sqrt( F_w ) * sqrt( wp2 ) &
        2278             :          * sqrt( F_x ) * sqrt( xp2 ) * sgn_wpxp &
        2279             :          * sqrt( F_y ) * sqrt( yp2 ) * sgn_wpyp &
        2280             :          * ( ( one - mixt_frac ) / mixt_frac - mixt_frac / ( one - mixt_frac ) &
        2281             :              - coefs_factor_xy - coefs_factor_wy - coefs_factor_wx ) &
        2282             :          + sqrt( mixt_frac * ( one - mixt_frac ) ) &
        2283             :            * sqrt( F_x ) * sqrt( xp2 ) * sgn_wpxp * coefs_factor_wy * wpyp &
        2284             :          + sqrt( mixt_frac * ( one - mixt_frac ) ) &
        2285             :            * sqrt( F_y ) * sqrt( yp2 ) * sgn_wpyp * coefs_factor_wx * wpxp
        2286             :
        2287             :     elsewhere ! coef_sigma_x_1_sqd * coef_sigma_y_1_sqd = 0
        2288             :               ! and coef_sigma_x_2_sqd * coef_sigma_y_2_sqd = 0
        2289             :
        2290             :        coef_wpxpyp_implicit &
        2291             :        = sqrt( mixt_frac * ( one - mixt_frac ) ) * sqrt( F_w ) * sqrt( wp2 ) &
        2292             :          * ( ( one - mixt_frac ) / mixt_frac - mixt_frac / ( one - mixt_frac ) &
        2293             :              - coefs_factor_wy - coefs_factor_wx )
        2294             :
        2295             :        term_wpxpyp_explicit &
        2296             :        = sqrt( mixt_frac * ( one - mixt_frac ) ) &
        2297             :          * sqrt( F_x ) * sqrt( xp2 ) * sgn_wpxp * coefs_factor_wy * wpyp &
        2298             :          + sqrt( mixt_frac * ( one - mixt_frac ) ) &
        2299             :            * sqrt( F_y ) * sqrt( yp2 ) * sgn_wpyp * coefs_factor_wx * wpxp
        2300             :
        2301             :     endwhere
        2302             :
        2303             :
        2304           0 :     return
        2305             :
        2306             :   end subroutine calc_coefs_wpxpyp_semiimpl
        2307             :
        2308             :   !=============================================================================
        2309             :
        2310             : end module new_pdf

Generated by: LCOV version 1.14