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

              Line data    Source code
           1             : ! $Id$
           2             : !===============================================================================
           3             : module new_hybrid_pdf_main
           4             :
           5             :   ! Description:
           6             :   ! The portion of CLUBB's multivariate, two-component PDF that is the
           7             :   ! multivariate, two-component normal PDF of vertical velocity (w), total water
           8             :   ! mixing ratio (rt), liquid water potential temperature (thl), and optionally,
           9             :   ! the west-east horizontal wind component (u), the south-north horizontal wind
          10             :   ! component (v), and passive scalars (sclr).
          11             :
          12             :   ! References:
          13             :   !-------------------------------------------------------------------------
          14             :
          15             :   implicit none
          16             :
          17             :   public :: new_hybrid_pdf_driver    ! Procedure(s)
          18             :
          19             :   private :: calc_responder_driver, & ! Procedure(s)
          20             :              calc_F_w_zeta_w
          21             :
          22             :   private
          23             :
          24             :   contains
          25             :
          26             :   !=============================================================================
          27           0 :   subroutine new_hybrid_pdf_driver( nz, ngrdcol, sclr_dim,              & ! In
          28           0 :                                     wm, rtm, thlm, um, vm,              & ! In
          29           0 :                                     wp2, rtp2, thlp2, up2, vp2,         & ! In
          30           0 :                                     Skw, wprtp, wpthlp, upwp, vpwp,     & ! In
          31           0 :                                     sclrm, sclrp2, wpsclrp,             & ! In
          32           0 :                                     gamma_Skw_fnc,                      & ! In
          33           0 :                                     slope_coef_spread_DG_means_w,       & ! In
          34           0 :                                     pdf_component_stdev_factor_w,       & ! In
          35           0 :                                     Skrt, Skthl, Sku, Skv, Sksclr,      & ! I/O
          36           0 :                                     mu_w_1, mu_w_2,                     & ! Out
          37           0 :                                     mu_rt_1, mu_rt_2,                   & ! Out
          38           0 :                                     mu_thl_1, mu_thl_2,                 & ! Out
          39           0 :                                     mu_u_1, mu_u_2, mu_v_1, mu_v_2,     & ! Out
          40           0 :                                     sigma_w_1_sqd, sigma_w_2_sqd,       & ! Out
          41           0 :                                     sigma_rt_1_sqd, sigma_rt_2_sqd,     & ! Out
          42           0 :                                     sigma_thl_1_sqd, sigma_thl_2_sqd,   & ! Out
          43           0 :                                     sigma_u_1_sqd, sigma_u_2_sqd,       & ! Out
          44           0 :                                     sigma_v_1_sqd, sigma_v_2_sqd,       & ! Out
          45           0 :                                     mu_sclr_1, mu_sclr_2,               & ! Out
          46           0 :                                     sigma_sclr_1_sqd, sigma_sclr_2_sqd, & ! Out
          47           0 :                                     mixt_frac,                          & ! Out
          48             :                                     pdf_implicit_coefs_terms,           & ! Out
          49           0 :                                     F_w, min_F_w, max_F_w               ) ! Out
          50             :
          51             :
          52             :     ! Description:
          53             :     ! Calculate the PDF parameters for w (including mixture fraction), rt,
          54             :     ! theta-l, and optionally, u, v, and passive scalar variables.
          55             :
          56             :     ! References:
          57             :     !-----------------------------------------------------------------------
          58             :
          59             :     use constants_clubb, only: &
          60             :         zero,     & ! Constant(s)
          61             :         fstderr
          62             :
          63             :     use new_hybrid_pdf, only: &
          64             :         calculate_w_params,          & ! Procedure(s)
          65             :         calculate_coef_wp4_implicit, &
          66             :         calc_coef_wp2xp_implicit,    &
          67             :         calc_coefs_wpxp2_semiimpl,   &
          68             :         calc_coefs_wpxpyp_semiimpl
          69             :
          70             :     use pdf_parameter_module, only: &
          71             :         implicit_coefs_terms    ! Variable Type
          72             :
          73             :     use model_flags, only: &
          74             :         l_explicit_turbulent_adv_wp3,  & ! Variable(s)
          75             :         l_explicit_turbulent_adv_wpxp, &
          76             :         l_explicit_turbulent_adv_xpyp
          77             :
          78             :     use clubb_precision, only: &
          79             :         core_rknd    ! Variable(s)
          80             :
          81             :     implicit none
          82             :
          83             :     integer, intent(in) :: &
          84             :       nz, &
          85             :       ngrdcol, &
          86             :       sclr_dim
          87             :
          88             :     ! Input Variables
          89             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
          90             :       wm,      & ! Mean of w (overall)                 [m/s]
          91             :       rtm,     & ! Mean of rt (overall)                [kg/kg]
          92             :       thlm,    & ! Mean of thl (overall)               [K]
          93             :       um,      & ! Mean of u (overall)                 [m/s]
          94             :       vm,      & ! Mean of v (overall)                 [m/s]
          95             :       wp2,     & ! Variance of w (overall)             [m^2/s^2]
          96             :       rtp2,    & ! Variance of rt (overall)            [kg^2/kg^2]
          97             :       thlp2,   & ! Variance of thl (overall)           [K^2]
          98             :       up2,     & ! Variance of u (overall)             [m^2/s^2]
          99             :       vp2,     & ! Variance of v (overall)             [m^2/s^2]
         100             :       Skw,     & ! Skewness of w (overall)             [-]
         101             :       wprtp,   & ! Covariance of w and rt (overall)    [(m/s)kg/kg]
         102             :       wpthlp,  & ! Covariance of w and thl (overall)   [(m/s)K]
         103             :       upwp,    & ! Covariance of u and w (overall)     [m^2/s^2]
         104             :       vpwp       ! Covariance of v and w (overall)     [m^2/s^2]
         105             :
         106             :     real( kind = core_rknd ), dimension(ngrdcol,nz,sclr_dim), intent(in) :: &
         107             :       sclrm,   & ! Mean of sclr (overall)                [units vary]
         108             :       sclrp2,  & ! Variance of sclr (overall)            [(units vary)^2]
         109             :       wpsclrp    ! Covariance of w and sclr (overall)    [(m/s)(units vary)]
         110             :
         111             :     ! Tunable parameter gamma.
         112             :     ! When gamma goes to 0, the standard deviations of w in each PDF component
         113             :     ! become small, and the spread between the two PDF component means of w
         114             :     ! becomes large.  F_w goes to min_F_w.
         115             :     ! When gamma goes to 1, the standard deviations of w in each PDF component
         116             :     ! become large, and the spread between the two PDF component means of w
         117             :     ! becomes small.  F_w goes to max_F_w.
         118             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
         119             :       gamma_Skw_fnc    ! Value of parameter gamma from tunable Skw function  [-]
         120             :
         121             :     real( kind = core_rknd ), dimension(ngrdcol), intent(in) :: &
         122             :       ! Slope coefficient for the spread between the PDF component means of w.
         123             :       slope_coef_spread_DG_means_w, &
         124             :       ! Parameter to adjust the PDF component standard deviations of w.
         125             :       pdf_component_stdev_factor_w
         126             :
         127             :     ! Input/Output Variables
         128             :     ! These variables are input/output because their values may be clipped.
         129             :     ! Otherwise, as long as it is not necessary to clip them, their values
         130             :     ! will stay the same.
         131             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(inout) :: &
         132             :       Skrt,  & ! Skewness of rt (overall)            [-]
         133             :       Skthl, & ! Skewness of thl (overall)           [-]
         134             :       Sku,   & ! Skewness of u (overall)             [-]
         135             :       Skv      ! Skewness of v (overall)             [-]
         136             :
         137             :     real( kind = core_rknd ), dimension(ngrdcol,nz,sclr_dim), intent(inout) :: &
         138             :       Sksclr    ! Skewness of sclr (overall)         [-]
         139             :
         140             :     ! Output Variables
         141             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
         142             :       mu_w_1,          & ! Mean of w (1st PDF component)        [m/s]
         143             :       mu_w_2,          & ! Mean of w (2nd PDF component)        [m/s]
         144             :       mu_rt_1,         & ! Mean of rt (1st PDF component)       [kg/kg]
         145             :       mu_rt_2,         & ! Mean of rt (2nd PDF component)       [kg/kg]
         146             :       mu_thl_1,        & ! Mean of thl (1st PDF component)      [K]
         147             :       mu_thl_2,        & ! Mean of thl (2nd PDF component)      [K]
         148             :       mu_u_1,          & ! Mean of u (1st PDF component)        [m/s]
         149             :       mu_u_2,          & ! Mean of u (2nd PDF component)        [m/s]
         150             :       mu_v_1,          & ! Mean of v (1st PDF component)        [m/s]
         151             :       mu_v_2,          & ! Mean of v (2nd PDF component)        [m/s]
         152             :       sigma_w_1_sqd,   & ! Variance of w (1st PDF component)    [m^2/s^2]
         153             :       sigma_w_2_sqd,   & ! Variance of w (2nd PDF component)    [m^2/s^2]
         154             :       sigma_rt_1_sqd,  & ! Variance of rt (1st PDF component)   [kg^2/kg^2]
         155             :       sigma_rt_2_sqd,  & ! Variance of rt (2nd PDF component)   [kg^2/kg^2]
         156             :       sigma_thl_1_sqd, & ! Variance of thl (1st PDF component)  [K^2]
         157             :       sigma_thl_2_sqd, & ! Variance of thl (2nd PDF component)  [K^2]
         158             :       sigma_u_1_sqd,   & ! Variance of u (1st PDF component)    [m^2/s^2]
         159             :       sigma_u_2_sqd,   & ! Variance of u (2nd PDF component)    [m^2/s^2]
         160             :       sigma_v_1_sqd,   & ! Variance of v (1st PDF component)    [m^2/s^2]
         161             :       sigma_v_2_sqd,   & ! Variance of v (2nd PDF component)    [m^2/s^2]
         162             :       mixt_frac          ! Mixture fraction                     [-]
         163             :
         164             :     real( kind = core_rknd ), dimension(ngrdcol,nz,sclr_dim), intent(out) :: &
         165             :       mu_sclr_1,        & ! Mean of sclr (1st PDF component)      [units vary]
         166             :       mu_sclr_2,        & ! Mean of sclr (2nd PDF component)      [units vary]
         167             :       sigma_sclr_1_sqd, & ! Variance of sclr (1st PDF component)  [(un. vary)^2]
         168             :       sigma_sclr_2_sqd    ! Variance of sclr (2nd PDF component)  [(un. vary)^2]
         169             :
         170             :     type(implicit_coefs_terms), intent(inout) :: &
         171             :       pdf_implicit_coefs_terms    ! Implicit coefs / explicit terms [units vary]
         172             :
         173             :     ! Output only for recording statistics.
         174             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
         175             :       F_w,     & ! Parameter for the spread of the PDF component means of w  [-]
         176             :       min_F_w, & ! Minimum allowable value of parameter F_w                  [-]
         177             :       max_F_w    ! Maximum allowable value of parameter F_w                  [-]
         178             :
         179             :     ! Local Variables
         180             :     real( kind = core_rknd ), dimension(nz) :: &
         181           0 :       sigma_w_1, & ! Standard deviation of w (1st PDF component)      [m/s]
         182           0 :       sigma_w_2    ! Standard deviation of w (2nd PDF component)      [m/s]
         183             :
         184             :     real( kind = core_rknd ), dimension(nz) :: &
         185           0 :       coef_sigma_w_1_sqd,   & ! sigma_w_1^2 = coef_sigma_w_1_sqd * <w'^2>    [-]
         186           0 :       coef_sigma_w_2_sqd,   & ! sigma_w_2^2 = coef_sigma_w_2_sqd * <w'^2>    [-]
         187           0 :       coef_sigma_rt_1_sqd,  & ! sigma_rt_1^2 = coef_sigma_rt_1_sqd * <rt'^2> [-]
         188           0 :       coef_sigma_rt_2_sqd,  & ! sigma_rt_2^2 = coef_sigma_rt_2_sqd * <rt'^2> [-]
         189           0 :       coef_sigma_thl_1_sqd, & ! sigma_thl_1^2=coef_sigma_thl_1_sqd*<thl'^2>  [-]
         190           0 :       coef_sigma_thl_2_sqd, & ! sigma_thl_2^2=coef_sigma_thl_2_sqd*<thl'^2>  [-]
         191           0 :       coef_sigma_u_1_sqd,   & ! sigma_u_1^2 = coef_sigma_u_1_sqd * <u'^2>    [-]
         192           0 :       coef_sigma_u_2_sqd,   & ! sigma_u_2^2 = coef_sigma_u_2_sqd * <u'^2>    [-]
         193           0 :       coef_sigma_v_1_sqd,   & ! sigma_v_1^2 = coef_sigma_v_1_sqd * <v'^2>    [-]
         194           0 :       coef_sigma_v_2_sqd      ! sigma_v_2^2 = coef_sigma_v_2_sqd * <v'^2>    [-]
         195             :
         196             :     ! sigma_sclr_1^2 = coef_sigma_sclr_1_sqd * <sclr'^2>
         197             :     ! sigma_sclr_2^2 = coef_sigma_sclr_2_sqd * <sclr'^2>
         198             :     real( kind = core_rknd ), dimension(nz,sclr_dim) :: &
         199           0 :       coef_sigma_sclr_1_sqd, & ! Coefficient that is multiplied by <sclr'^2> [-]
         200           0 :       coef_sigma_sclr_2_sqd    ! Coefficient that is multiplied by <sclr'^2> [-]
         201             :
         202             :     real( kind = core_rknd ), dimension(nz) :: &
         203           0 :       zeta_w    ! Parameter for the PDF component variances of w           [-]
         204             :
         205             :     real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
         206             :       ! Slope coefficient for the spread between the PDF component means of w.
         207           0 :       slope_coef_spread_DG_means_w_in, &
         208             :       ! Parameter to adjust the PDF component standard deviations of w.
         209           0 :       pdf_component_stdev_factor_w_in
         210             :
         211             :     real( kind = core_rknd ), dimension(nz) :: &
         212           0 :       coef_wp4_implicit,     & ! <w'^4> = coef_wp4_implicit * <w'^2>^2     [-]
         213           0 :       coef_wp2rtp_implicit,  & ! <w'^2 rt'> = coef_wp2rtp_implicit*<w'rt'> [m/s]
         214           0 :       coef_wp2thlp_implicit, & ! <w'^2 thl'>=coef_wp2thlp_implicit*<w'thl'>[m/s]
         215           0 :       coef_wp2up_implicit,   & ! <w'^2 u'> = coef_wp2up_implicit * <u'w'>  [m/s]
         216           0 :       coef_wp2vp_implicit      ! <w'^2 v'> = coef_wp2vp_implicit * <v'w'>  [m/s]
         217             :
         218             :     ! <w'^2 sclr'> = coef_wp2sclrp_implicit * <w'sclr'>
         219             :     real( kind = core_rknd ), dimension(nz,sclr_dim) :: &
         220           0 :       coef_wp2sclrp_implicit    ! Coef. that is multiplied by <w'sclr'>    [m/s]
         221             :
         222             :     ! <w'rt'^2> = coef_wprtp2_implicit * <rt'^2> + term_wprtp2_explicit
         223             :     real( kind = core_rknd ), dimension(nz) :: &
         224           0 :       coef_wprtp2_implicit, & ! Coefficient that is multiplied by <rt'^2>  [m/s]
         225           0 :       term_wprtp2_explicit    ! Term that is on the RHS          [m/s kg^2/kg^2]
         226             :
         227             :     ! <w'thl'^2> = coef_wpthlp2_implicit * <thl'^2> + term_wpthlp2_explicit
         228             :     real( kind = core_rknd ), dimension(nz) :: &
         229           0 :       coef_wpthlp2_implicit, & ! Coef. that is multiplied by <thl'^2>      [m/s]
         230           0 :       term_wpthlp2_explicit    ! Term that is on the RHS               [m/s K^2]
         231             :
         232             :     ! <w'rt'thl'> = coef_wprtpthlp_implicit*<rt'thl'> + term_wprtpthlp_explicit
         233             :     real( kind = core_rknd ), dimension(nz) :: &
         234           0 :       coef_wprtpthlp_implicit, & ! Coef. that is multiplied by <rt'thl'>   [m/s]
         235           0 :       term_wprtpthlp_explicit    ! Term that is on the RHS         [m/s(kg/kg)K]
         236             :
         237             :     ! <w'u'^2> = coef_wpup2_implicit * <u'^2> + term_wpup2_explicit
         238             :     real( kind = core_rknd ), dimension(nz) :: &
         239           0 :       coef_wpup2_implicit, & ! Coefficient that is multiplied by <u'^2>    [m/s]
         240           0 :       term_wpup2_explicit    ! Term that is on the RHS                 [m^3/s^3]
         241             :
         242             :     ! <w'v'^2> = coef_wpvp2_implicit * <v'^2> + term_wpvp2_explicit
         243             :     real( kind = core_rknd ), dimension(nz) :: &
         244           0 :       coef_wpvp2_implicit, & ! Coefficient that is multiplied by <v'^2>    [m/s]
         245           0 :       term_wpvp2_explicit    ! Term that is on the RHS                 [m^3/s^3]
         246             :
         247             :     ! <w'sclr'^2> = coef_wpsclrp2_implicit * <sclr'^2> + term_wpsclrp2_explicit
         248             :     real( kind = core_rknd ), dimension(nz,sclr_dim) :: &
         249           0 :       coef_wpsclrp2_implicit, & ! Coef. that is multiplied by <sclr'^2>    [m/s]
         250           0 :       term_wpsclrp2_explicit    ! Term that is on the RHS    [m/s(units vary)^2]
         251             :
         252             :     ! <w'rt'sclr'> = coef_wprtpsclrp_implicit * <sclr'rt'>
         253             :     !                + term_wprtpsclrp_explicit
         254             :     real( kind = core_rknd ), dimension(nz,sclr_dim) :: &
         255           0 :       coef_wprtpsclrp_implicit, & ! Coef. that is multiplied by <sclr'rt'> [m/s]
         256           0 :       term_wprtpsclrp_explicit    ! Term that is on the RHS [m/s(kg/kg)(un. v.)]
         257             :
         258             :     ! <w'thl'sclr'> = coef_wpthlpsclrp_implicit * <sclr'thl'>
         259             :     !                 + term_wpthlpsclrp_explicit
         260             :     real( kind = core_rknd ), dimension(nz,sclr_dim) :: &
         261           0 :       coef_wpthlpsclrp_implicit, & ! Coef. that is mult. by <sclr'thl'>    [m/s]
         262           0 :       term_wpthlpsclrp_explicit    ! Term that is on the RHS  [(m/s)K(un. vary)]
         263             :
         264             :     ! Variables to interface with code for the jth scalar
         265             :     real( kind = core_rknd ), dimension(nz) :: &
         266           0 :       sclrjm,   & ! Mean of sclr j (overall)                [units vary]
         267           0 :       sclrjp2,  & ! Variance of sclr j (overall)            [(units vary)^2]
         268           0 :       wpsclrjp, & ! Covariance of w and sclr j (overall)    [(m/s)(units vary)]
         269           0 :       Sksclrj     ! Skewness of rt (overall)                [-]
         270             :
         271             :     real( kind = core_rknd ), dimension(nz) :: &
         272           0 :       mu_sclrj_1,        & ! Mean of sclr j (1st PDF component) [units vary]
         273           0 :       mu_sclrj_2,        & ! Mean of sclr j (2nd PDF component) [units vary]
         274           0 :       sigma_sclrj_1_sqd, & ! Variance of sclr j (1st PDF comp.) [(units vary)^2]
         275           0 :       sigma_sclrj_2_sqd    ! Variance of sclr j (2nd PDF comp.) [(units vary)^2]
         276             :
         277             :     ! sigma_sclrj_1^2 = coef_sigma_sclrj_1_sqd * <sclrj'^2>
         278             :     ! sigma_sclrj_2^2 = coef_sigma_sclrj_2_sqd * <sclrj'^2>
         279             :     real( kind = core_rknd ), dimension(nz) :: &
         280           0 :       coef_sigma_sclrj_1_sqd, & ! Coef. that is multiplied by <sclrj'^2>    [-]
         281           0 :       coef_sigma_sclrj_2_sqd    ! Coef. that is multiplied by <sclrj'^2>    [-]
         282             :
         283             :     ! <w'sclrj'^2> = coef_wpsclrjp2_implicit * <sclrj'^2>
         284             :     !                + term_wpsclrjp2_explicit
         285             :     real( kind = core_rknd ), dimension(nz) :: &
         286           0 :       coef_wpsclrjp2_implicit, & ! Coef. that is multiplied by <sclrj'^2>  [m/s]
         287           0 :       term_wpsclrjp2_explicit    ! Term that is on the RHS   [m/s(units vary)^2]
         288             :
         289             :     ! <w'rt'sclrj'> = coef_wprtpsclrjp_implicit * <sclrj'rt'>
         290             :     !                 + term_wprtpsclrjp_explicit
         291             :     real( kind = core_rknd ), dimension(nz) :: &
         292           0 :       coef_wprtpsclrjp_implicit, & ! Coef. that is mult. by <sclr'rt'>     [m/s]
         293           0 :       term_wprtpsclrjp_explicit    ! Term that is on the RHS [m/s(kg/kg)(un.v.)]
         294             :
         295             :     ! <w'thl'sclrj'> = coef_wpthlpsclrjp_implicit * <sclrj'thl'>
         296             :     !                  + term_wpthlpsclrjp_explicit
         297             :     real( kind = core_rknd ), dimension(nz) :: &
         298           0 :       coef_wpthlpsclrjp_implicit, & ! Coef. that is mult. by <sclrj'thl'>  [m/s]
         299           0 :       term_wpthlpsclrjp_explicit    ! Term that is on the RHS [(m/s)K(un. vary)]
         300             :
         301             :     real( kind = core_rknd ), dimension(nz) :: &
         302           0 :       max_corr_w_sclr_sqd    ! Max value of wpsclrp^2 / ( wp2 * sclrp2 )     [-]
         303             :
         304             :     real( kind = core_rknd ), dimension(nz) :: &
         305           0 :       zeros    ! Vector of 0s (size nz)    [-]
         306             :
         307             :     real( kind = core_rknd ), dimension(nz,sclr_dim) :: &
         308           0 :       zero_array    ! Array of 0s (size nz x sclr_dim)    [-]
         309             :
         310             :     integer :: i, k, j  ! Loop indices
         311             :
         312           0 :     do i = 1, ngrdcol
         313             :
         314             :       ! Calculate the maximum value of the square of the correlation of w and a
         315             :       ! scalar when scalars are used.
         316             :       ! Initialize max_corr_w_sclr_sqd to 0.  It needs to retain this value even
         317             :       ! when sclr_dim = 0.
         318           0 :       max_corr_w_sclr_sqd = zero
         319           0 :       if ( sclr_dim > 0 ) then
         320           0 :          do k = 1, nz, 1
         321           0 :             do j = 1, sclr_dim, 1
         322           0 :                if ( wp2(i,k) * sclrp2(i,k,j) > zero ) then
         323             :                   max_corr_w_sclr_sqd(k) = max( wpsclrp(i,k,j)**2 &
         324             :                                                 / ( wp2(i,k) * sclrp2(i,k,j) ), &
         325           0 :                                                 max_corr_w_sclr_sqd(k) )
         326             :                endif ! wp2(k) * sclrp2(k,j) > 0
         327             :             enddo ! j = 1, sclr_dim, 1
         328             :          enddo ! k = 1, nz, 1
         329             :       endif ! sclr_dim > 0
         330             :
         331           0 :       slope_coef_spread_DG_means_w_in(i,:) = slope_coef_spread_DG_means_w(i)
         332           0 :       pdf_component_stdev_factor_w_in(i,:) = pdf_component_stdev_factor_w(i)
         333             :
         334             :       ! Calculate the values of PDF tunable parameters F_w and zeta_w.
         335             :       ! Vertical velocity, w, will always be the setter variable.
         336             :       call calc_F_w_zeta_w( Skw(i,:), wprtp(i,:), wpthlp(i,:), upwp(i,:), vpwp(i,:),  & ! In
         337             :                             wp2(i,:), rtp2(i,:), thlp2(i,:), up2(i,:), vp2(i,:),      & ! In
         338             :                             gamma_Skw_fnc(i,:),                   & ! In
         339             :                             slope_coef_spread_DG_means_w_in(i,:), & ! In
         340             :                             pdf_component_stdev_factor_w_in(i,:), & ! In
         341             :                             max_corr_w_sclr_sqd,             & ! In
         342           0 :                             F_w(i,:), zeta_w, min_F_w(i,:), max_F_w(i,:)    ) ! Out
         343             :
         344             :       ! Calculate the PDF parameters, including mixture fraction, for the
         345             :       ! setter variable, w.
         346             :       call calculate_w_params( wm(i,:), wp2(i,:), Skw(i,:), F_w(i,:), zeta_w, & ! In
         347             :                                mu_w_1(i,:), mu_w_2(i,:), sigma_w_1, & ! Out
         348             :                                sigma_w_2, mixt_frac(i,:),      & ! Out
         349             :                                coef_sigma_w_1_sqd,        & ! Out
         350           0 :                                coef_sigma_w_2_sqd         ) ! Out
         351             :
         352           0 :       sigma_w_1_sqd(i,:) = sigma_w_1**2
         353           0 :       sigma_w_2_sqd(i,:) = sigma_w_2**2
         354             :
         355           0 :       if ( any( mixt_frac(i,:) < zero ) ) then
         356           0 :          write(fstderr,*) "Mixture fraction cannot be calculated."
         357             :          write(fstderr,*) "The value of F_w must be greater than 0 when " &
         358           0 :                           // "| Skw | > 0."
         359           0 :          error stop
         360             :       endif
         361             :
         362             :       ! Calculate the PDF parameters for responder variable rt.
         363             :       call calc_responder_driver( rtm(i,:), rtp2(i,:), wprtp(i,:), wp2(i,:), & ! In
         364             :                                   mixt_frac(i,:), F_w(i,:),        & ! In
         365             :                                   Skrt(i,:),                  & ! In/Out
         366             :                                   mu_rt_1(i,:), mu_rt_2(i,:),      & ! Out
         367             :                                   sigma_rt_1_sqd(i,:),        & ! Out
         368             :                                   sigma_rt_2_sqd(i,:),        & ! Out
         369             :                                   coef_sigma_rt_1_sqd,   & ! Out
         370           0 :                                   coef_sigma_rt_2_sqd    ) ! Out
         371             :
         372             :       ! Calculate the PDF parameters for responder variable thl.
         373             :       call calc_responder_driver( thlm(i,:), thlp2(i,:), wpthlp(i,:), wp2(i,:), & ! In
         374             :                                   mixt_frac(i,:), F_w(i,:),           & ! In
         375             :                                   Skthl(i,:),                    & ! In/Out
         376             :                                   mu_thl_1(i,:), mu_thl_2(i,:),       & ! Out
         377             :                                   sigma_thl_1_sqd(i,:),          & ! Out
         378             :                                   sigma_thl_2_sqd(i,:),          & ! Out
         379             :                                   coef_sigma_thl_1_sqd,     & ! Out
         380           0 :                                   coef_sigma_thl_2_sqd      ) ! Out
         381             :
         382             :       ! Calculate the PDF parameters for responder variable u.
         383             :       call calc_responder_driver( um(i,:), up2(i,:), upwp(i,:), wp2(i,:), & ! In
         384             :                                   mixt_frac(i,:), F_w(i,:),     & ! In
         385             :                                   Sku(i,:),                & ! In/Out
         386             :                                   mu_u_1(i,:), mu_u_2(i,:),     & ! Out
         387             :                                   sigma_u_1_sqd(i,:),      & ! Out
         388             :                                   sigma_u_2_sqd(i,:),      & ! Out
         389             :                                   coef_sigma_u_1_sqd, & ! Out
         390           0 :                                   coef_sigma_u_2_sqd  ) ! Out
         391             :
         392             :       ! Calculate the PDF parameters for responder variable v.
         393             :       call calc_responder_driver( vm(i,:), vp2(i,:), vpwp(i,:), wp2(i,:), & ! In
         394             :                                   mixt_frac(i,:), F_w(i,:),     & ! In
         395             :                                   Skv(i,:),                & ! In/Out
         396             :                                   mu_v_1(i,:), mu_v_2(i,:),     & ! Out
         397             :                                   sigma_v_1_sqd(i,:),      & ! Out
         398             :                                   sigma_v_2_sqd(i,:),      & ! Out
         399             :                                   coef_sigma_v_1_sqd, & ! Out
         400           0 :                                   coef_sigma_v_2_sqd  ) ! Out
         401             :
         402             :       ! Calculate the PDF parameters for responder variable sclr.
         403           0 :       if ( sclr_dim > 0 ) then
         404             :
         405           0 :          do j = 1, sclr_dim, 1
         406             :
         407           0 :             do k = 1, nz, 1
         408           0 :                sclrjm(k) = sclrm(i,k,j)
         409           0 :                sclrjp2(k) = sclrp2(i,k,j)
         410           0 :                wpsclrjp(k) = wpsclrp(i,k,j)
         411           0 :                Sksclrj(k) = Sksclr(i,k,j)
         412             :             enddo ! k = 1, nz, 1
         413             :
         414             :             call calc_responder_driver( sclrjm, sclrjp2, wpsclrjp, wp2(i,:), & ! In
         415             :                                         mixt_frac(i,:), F_w(i,:),                 & ! In
         416             :                                         Sksclrj,                        & ! In/Out
         417             :                                         mu_sclrj_1, mu_sclrj_2,         & ! Out
         418             :                                         sigma_sclrj_1_sqd,              & ! Out
         419             :                                         sigma_sclrj_2_sqd,              & ! Out
         420             :                                         coef_sigma_sclrj_1_sqd,         & ! Out
         421           0 :                                         coef_sigma_sclrj_2_sqd          ) ! Out
         422             :
         423           0 :             do k = 1, nz, 1
         424           0 :                Sksclr(i,k,j) = Sksclrj(k)
         425           0 :                mu_sclr_1(i,k,j) = mu_sclrj_1(k)
         426           0 :                mu_sclr_2(i,k,j) = mu_sclrj_2(k)
         427           0 :                sigma_sclr_1_sqd(i,k,j) = sigma_sclrj_1_sqd(k)
         428           0 :                sigma_sclr_2_sqd(i,k,j) = sigma_sclrj_2_sqd(k)
         429           0 :                coef_sigma_sclr_1_sqd(k,j) = coef_sigma_sclrj_1_sqd(k)
         430           0 :                coef_sigma_sclr_2_sqd(k,j) = coef_sigma_sclrj_2_sqd(k)
         431             :             enddo ! k = 1, nz, 1
         432             :
         433             :          enddo ! j = 1, sclr_dim, 1
         434             :
         435             :       endif ! sclr_dim > 0
         436             :
         437             :
         438             :       if ( .not. l_explicit_turbulent_adv_wp3 ) then
         439             :
         440             :          ! Turbulent advection of <w'^3> is being handled implicitly.
         441             :
         442             :          ! <w'^4> = coef_wp4_implicit * <w'^2>^2.
         443             :          coef_wp4_implicit &
         444             :          = calculate_coef_wp4_implicit( mixt_frac(i,:), F_w(i,:), &
         445             :                                         coef_sigma_w_1_sqd, &
         446           0 :                                         coef_sigma_w_2_sqd )
         447             :
         448             :       else ! l_explicit_turbulent_adv_wp3
         449             :
         450             :          ! Turbulent advection of <w'^3> is being handled explicitly.
         451             :          coef_wp4_implicit = zero
         452             :
         453             :       endif ! .not. l_explicit_turbulent_adv_wp3
         454             :
         455             :       if ( .not. l_explicit_turbulent_adv_wpxp ) then
         456             :
         457             :          ! Turbulent advection of <w'rt'>, <w'thl'>, <u'w'>, <v'w'>, and <w'sclr'>
         458             :          ! are all being handled implicitly.
         459             :
         460             :          ! <w'^2 rt'> = coef_wp2rtp_implicit * <w'rt'>
         461             :          coef_wp2rtp_implicit = calc_coef_wp2xp_implicit( wp2(i,:), mixt_frac(i,:), F_w(i,:), &
         462             :                                                           coef_sigma_w_1_sqd,  &
         463           0 :                                                           coef_sigma_w_2_sqd   )
         464             :
         465             :          ! <w'^2 thl'> = coef_wp2thlp_implicit * <w'thl'>;
         466             :          ! <w'^2 u'> = coef_wp2up_implicit * <u'w'>; and
         467             :          ! <w'^2 v'> = coef_wp2vp_implicit * <v'w'>;
         468             :          ! where each coef_wp2xp_implicit is the same as coef_wp2rtp_implicit.
         469           0 :          coef_wp2thlp_implicit = coef_wp2rtp_implicit
         470           0 :          coef_wp2up_implicit = coef_wp2rtp_implicit
         471           0 :          coef_wp2vp_implicit = coef_wp2rtp_implicit
         472             :
         473             :          ! <w'^2 sclr'> = coef_wp2sclrp_implicit * <w'sclr'>;
         474             :          ! where each coef_wp2xp_implicit is the same as coef_wp2rtp_implicit.
         475           0 :          if ( sclr_dim > 0 ) then
         476           0 :             do k = 1, nz, 1
         477           0 :                do j = 1, sclr_dim, 1
         478           0 :                   coef_wp2sclrp_implicit(k,j) = coef_wp2rtp_implicit(k)
         479             :                enddo ! j = 1, sclr_dim, 1
         480             :             enddo ! k = 1, nz, 1
         481             :          endif ! sclr_dim > 0
         482             :
         483             :       else ! l_explicit_turbulent_adv_wpxp
         484             :
         485             :          ! Turbulent advection of <w'rt'>, <w'thl'>, <u'w'>, <v'w'>, and <w'sclr'>
         486             :          ! are all being handled explicitly.
         487             :          coef_wp2rtp_implicit = zero
         488             :          coef_wp2thlp_implicit = zero
         489             :          coef_wp2up_implicit = zero
         490             :          coef_wp2vp_implicit = zero
         491             :          if ( sclr_dim > 0 ) then
         492             :             coef_wp2sclrp_implicit = zero
         493             :          endif ! sclr_dim > 0
         494             :
         495             :       endif ! .not. l_explicit_turbulent_adv_wpxp
         496             :
         497             :       if ( .not. l_explicit_turbulent_adv_xpyp ) then
         498             :
         499             :          ! Turbulent advection of <rt'^2>, <thl'^2>, <rt'thl'>, <u'^2>, <v'^2>,
         500             :          ! <sclr'^2>, <sclr'rt'>, and <sclr'thl'> are all being handled
         501             :          ! semi-implicitly.
         502             :
         503             :          ! <w'rt'^2> = coef_wprtp2_implicit * <rt'^2> + term_wprtp2_explicit
         504             :          call calc_coefs_wpxp2_semiimpl( wp2(i,:), wprtp(i,:),           & ! In
         505             :                                          mixt_frac(i,:), F_w(i,:),       & ! In
         506             :                                          coef_sigma_rt_1_sqd,  & ! In
         507             :                                          coef_sigma_rt_2_sqd,  & ! In
         508             :                                          coef_wprtp2_implicit, & ! Out
         509           0 :                                          term_wprtp2_explicit )  ! Out
         510             :
         511             :          ! <w'thl'^2> = coef_wpthlp2_implicit * <thl'^2> + term_wprtp2_explicit
         512             :          call calc_coefs_wpxp2_semiimpl( wp2(i,:), wpthlp(i,:),           & ! In
         513             :                                          mixt_frac(i,:), F_w(i,:),        & ! In
         514             :                                          coef_sigma_thl_1_sqd,  & ! In
         515             :                                          coef_sigma_thl_2_sqd,  & ! In
         516             :                                          coef_wpthlp2_implicit, & ! Out
         517           0 :                                          term_wpthlp2_explicit )  ! Out
         518             :
         519             :          ! <w'rt'thl'> = coef_wprtpthlp_implicit * <rt'thl'>
         520             :          !               + term_wprtpthlp_explicit
         521             :          call calc_coefs_wpxpyp_semiimpl( wp2(i,:), wprtp(i,:), wpthlp(i,:),      & ! In
         522             :                                           mixt_frac(i,:), F_w(i,:),          & ! In
         523             :                                           coef_sigma_rt_1_sqd,     & ! In
         524             :                                           coef_sigma_rt_2_sqd,     & ! In
         525             :                                           coef_sigma_thl_1_sqd,    & ! In
         526             :                                           coef_sigma_thl_2_sqd,    & ! In
         527             :                                           coef_wprtpthlp_implicit, & ! Out
         528           0 :                                           term_wprtpthlp_explicit  ) ! Out
         529             :
         530             :          ! <w'u'^2> = coef_wpup2_implicit * <u'^2> + term_wpup2_explicit
         531             :          call calc_coefs_wpxp2_semiimpl( wp2(i,:), upwp(i,:),           & ! In
         532             :                                          mixt_frac(i,:), F_w(i,:),      & ! In
         533             :                                          coef_sigma_u_1_sqd,  & ! In
         534             :                                          coef_sigma_u_2_sqd,  & ! In
         535             :                                          coef_wpup2_implicit, & ! Out
         536           0 :                                          term_wpup2_explicit )  ! Out
         537             :
         538             :          ! <w'v'^2> = coef_wpvp2_implicit * <v'^2> + term_wpvp2_explicit
         539             :          call calc_coefs_wpxp2_semiimpl( wp2(i,:), vpwp(i,:),           & ! In
         540             :                                          mixt_frac(i,:), F_w(i,:),      & ! In
         541             :                                          coef_sigma_v_1_sqd,  & ! In
         542             :                                          coef_sigma_v_2_sqd,  & ! In
         543             :                                          coef_wpvp2_implicit, & ! Out
         544           0 :                                          term_wpvp2_explicit )  ! Out
         545             :
         546           0 :          if ( sclr_dim > 0 ) then
         547             :
         548           0 :            do j = 1, sclr_dim, 1
         549             :
         550           0 :               do k = 1, nz, 1
         551           0 :                  wpsclrjp(k) = wpsclrp(i,k,j)
         552           0 :                  coef_sigma_sclrj_1_sqd(k) = coef_sigma_sclr_1_sqd(k,j)
         553           0 :                  coef_sigma_sclrj_2_sqd(k) = coef_sigma_sclr_2_sqd(k,j)
         554             :               enddo ! k = 1, nz, 1
         555             :
         556             :               ! <w'sclr'^2> = coef_wpsclrp2_implicit * <sclr'^2>
         557             :               !               + term_wpsclrp2_explicit
         558             :               call calc_coefs_wpxp2_semiimpl( wp2(i,:), wpsclrjp,           & ! In
         559             :                                               mixt_frac(i,:), F_w(i,:),          & ! In
         560             :                                               coef_sigma_sclrj_1_sqd,  & ! In
         561             :                                               coef_sigma_sclrj_2_sqd,  & ! In
         562             :                                               coef_wpsclrjp2_implicit, & ! Out
         563           0 :                                               term_wpsclrjp2_explicit )  ! Out
         564             :
         565             :               ! <w'rt'sclr'> = coef_wprtpsclrp_implicit * <sclr'rt'>
         566             :               !                + term_wprtpsclrp_explicit
         567             :               call calc_coefs_wpxpyp_semiimpl( wp2(i,:), wprtp(i,:), wpsclrjp,      & ! In
         568             :                                                mixt_frac(i,:), F_w(i,:),            & ! In
         569             :                                                coef_sigma_rt_1_sqd,       & ! In
         570             :                                                coef_sigma_rt_2_sqd,       & ! In
         571             :                                                coef_sigma_sclrj_1_sqd,    & ! In
         572             :                                                coef_sigma_sclrj_2_sqd,    & ! In
         573             :                                                coef_wprtpsclrjp_implicit, & ! Out
         574           0 :                                                term_wprtpsclrjp_explicit  ) ! Out
         575             :
         576             :               ! <w'thl'sclr'> = coef_wpthlpsclrp_implicit * <sclr'thl'>
         577             :               !                 + term_wpthlpsclrp_explicit
         578             :               call calc_coefs_wpxpyp_semiimpl( wp2(i,:), wpthlp(i,:), wpsclrjp,      & ! In
         579             :                                                mixt_frac(i,:), F_w(i,:),             & ! In
         580             :                                                coef_sigma_thl_1_sqd,       & ! In
         581             :                                                coef_sigma_thl_2_sqd,       & ! In
         582             :                                                coef_sigma_sclrj_1_sqd,     & ! In
         583             :                                                coef_sigma_sclrj_2_sqd,     & ! In
         584             :                                                coef_wpthlpsclrjp_implicit, & ! Out
         585           0 :                                                term_wpthlpsclrjp_explicit  ) ! Out
         586             :
         587           0 :               do k = 1, nz, 1
         588           0 :                  coef_wpsclrp2_implicit(k,j) = coef_wpsclrjp2_implicit(k)
         589           0 :                  term_wpsclrp2_explicit(k,j) = term_wpsclrjp2_explicit(k)
         590           0 :                  coef_wprtpsclrp_implicit(k,j) = coef_wprtpsclrjp_implicit(k)
         591           0 :                  term_wprtpsclrp_explicit(k,j) = term_wprtpsclrjp_explicit(k)
         592           0 :                  coef_wpthlpsclrp_implicit(k,j) = coef_wpthlpsclrjp_implicit(k)
         593           0 :                  term_wpthlpsclrp_explicit(k,j) = term_wpthlpsclrjp_explicit(k)
         594             :               enddo ! k = 1, nz, 1
         595             :
         596             :            enddo ! j = 1, sclr_dim, 1
         597             :
         598             :          endif ! sclr_dim > 0
         599             :
         600             :       else ! l_explicit_turbulent_adv_xpyp
         601             :
         602             :          ! Turbulent advection of <rt'^2>, <thl'^2>, <rt'thl'>, <u'^2>, <v'^2>,
         603             :          ! <sclr'^2>, <sclr'rt'>, and <sclr'thl'> are being handled explicitly.
         604             :          coef_wprtp2_implicit = zero
         605             :          term_wprtp2_explicit = zero
         606             :          coef_wpthlp2_implicit = zero
         607             :          term_wpthlp2_explicit = zero
         608             :          coef_wprtpthlp_implicit = zero
         609             :          term_wprtpthlp_explicit = zero
         610             :          coef_wpup2_implicit = zero
         611             :          term_wpup2_explicit = zero
         612             :          coef_wpvp2_implicit = zero
         613             :          term_wpvp2_explicit = zero
         614             :          if ( sclr_dim > 0 ) then
         615             :             coef_wpsclrp2_implicit = zero
         616             :             term_wpsclrp2_explicit = zero
         617             :             coef_wprtpsclrp_implicit = zero
         618             :             term_wprtpsclrp_explicit = zero
         619             :             coef_wpthlpsclrp_implicit = zero
         620             :             term_wpthlpsclrp_explicit = zero
         621             :          endif ! sclr_dim > 0
         622             :
         623             :       endif ! .not. l_explicit_turbulent_adv_xpyp
         624             :
         625             :       ! Set up a vector of 0s and an array of 0s to help write results back to
         626             :       ! type variable pdf_implicit_coefs_terms.
         627           0 :       zeros = zero
         628           0 :       zero_array = zero
         629             :
         630             :       ! Pack the implicit coefficients and explicit terms into a single type
         631             :       ! variable for output.
         632           0 :       pdf_implicit_coefs_terms%coef_wp4_implicit(i,:) = coef_wp4_implicit
         633           0 :       pdf_implicit_coefs_terms%coef_wp2rtp_implicit(i,:) = coef_wp2rtp_implicit
         634           0 :       pdf_implicit_coefs_terms%term_wp2rtp_explicit(i,:) = zeros
         635           0 :       pdf_implicit_coefs_terms%coef_wp2thlp_implicit(i,:) = coef_wp2thlp_implicit
         636           0 :       pdf_implicit_coefs_terms%term_wp2thlp_explicit(i,:) = zeros
         637           0 :       pdf_implicit_coefs_terms%coef_wp2up_implicit(i,:) = coef_wp2up_implicit
         638           0 :       pdf_implicit_coefs_terms%term_wp2up_explicit(i,:) = zeros
         639           0 :       pdf_implicit_coefs_terms%coef_wp2vp_implicit(i,:) = coef_wp2vp_implicit
         640           0 :       pdf_implicit_coefs_terms%term_wp2vp_explicit(i,:) = zeros
         641           0 :       pdf_implicit_coefs_terms%coef_wprtp2_implicit(i,:) = coef_wprtp2_implicit
         642           0 :       pdf_implicit_coefs_terms%term_wprtp2_explicit(i,:) = term_wprtp2_explicit
         643           0 :       pdf_implicit_coefs_terms%coef_wpthlp2_implicit(i,:) = coef_wpthlp2_implicit
         644           0 :       pdf_implicit_coefs_terms%term_wpthlp2_explicit(i,:) = term_wpthlp2_explicit
         645           0 :       pdf_implicit_coefs_terms%coef_wprtpthlp_implicit(i,:) = coef_wprtpthlp_implicit
         646           0 :       pdf_implicit_coefs_terms%term_wprtpthlp_explicit(i,:) = term_wprtpthlp_explicit
         647           0 :       pdf_implicit_coefs_terms%coef_wpup2_implicit(i,:) = coef_wpup2_implicit
         648           0 :       pdf_implicit_coefs_terms%term_wpup2_explicit(i,:) = term_wpup2_explicit
         649           0 :       pdf_implicit_coefs_terms%coef_wpvp2_implicit(i,:) = coef_wpvp2_implicit
         650           0 :       pdf_implicit_coefs_terms%term_wpvp2_explicit(i,:) = term_wpvp2_explicit
         651           0 :       if ( sclr_dim > 0 ) then
         652           0 :          pdf_implicit_coefs_terms%coef_wp2sclrp_implicit(i,:,:) = coef_wp2sclrp_implicit
         653           0 :          pdf_implicit_coefs_terms%term_wp2sclrp_explicit(i,:,:) = zero_array
         654           0 :          pdf_implicit_coefs_terms%coef_wpsclrp2_implicit(i,:,:) = coef_wpsclrp2_implicit
         655           0 :          pdf_implicit_coefs_terms%term_wpsclrp2_explicit(i,:,:) = term_wpsclrp2_explicit
         656           0 :          pdf_implicit_coefs_terms%coef_wprtpsclrp_implicit(i,:,:) &
         657           0 :          = coef_wprtpsclrp_implicit
         658           0 :          pdf_implicit_coefs_terms%term_wprtpsclrp_explicit(i,:,:) &
         659           0 :          = term_wprtpsclrp_explicit
         660           0 :          pdf_implicit_coefs_terms%coef_wpthlpsclrp_implicit(i,:,:) &
         661           0 :          = coef_wpthlpsclrp_implicit
         662           0 :          pdf_implicit_coefs_terms%term_wpthlpsclrp_explicit(i,:,:) &
         663           0 :          = term_wpthlpsclrp_explicit
         664             :       endif ! sclr_dim > 0
         665             :
         666             :     end do
         667             :
         668           0 :     return
         669             :
         670             :   end subroutine new_hybrid_pdf_driver
         671             :
         672             :   !=============================================================================
         673           0 :   elemental subroutine calc_responder_driver( xm, xp2, wpxp, wp2, & ! In
         674             :                                               mixt_frac, F_w,     & ! In
         675             :                                               Skx,                & ! In/Out
         676             :                                               mu_x_1, mu_x_2,     & ! Out
         677             :                                               sigma_x_1_sqd,      & ! Out
         678             :                                               sigma_x_2_sqd,      & ! Out
         679             :                                               coef_sigma_x_1_sqd, & ! Out
         680             :                                               coef_sigma_x_2_sqd  ) ! Out
         681             :
         682             :     ! Description:
         683             :     ! This is the sub-driver for a responder variable.  The limits of the range
         684             :     ! of values of Skx that the PDF is able to represent are calculated, and
         685             :     ! Skx is clipped when it exceeds the upper or lower limit.  Then, the PDF
         686             :     ! parameters for responder variable x are calculated.
         687             :
         688             :     ! References:
         689             :     !-----------------------------------------------------------------------
         690             :
         691             :     use constants_clubb, only: &
         692             :         three,        & ! Constant(s)
         693             :         two,          &
         694             :         three_halves, &
         695             :         one,          &
         696             :         zero
         697             :
         698             :     use new_hybrid_pdf, only: &
         699             :         calculate_responder_params    ! Procedure(s)
         700             :
         701             :     use clubb_precision, only: &
         702             :         core_rknd    ! Variable(s)
         703             :
         704             :     implicit none
         705             :
         706             :     ! Input Variables
         707             :     real( kind = core_rknd ), intent(in) :: &
         708             :       xm,        & ! Mean of x (overall)                          [units vary]
         709             :       xp2,       & ! Variance of x (overall)                  [(units vary)^2]
         710             :       wpxp,      & ! Covariance of w and x (overall)        [m/s (units vary)]
         711             :       wp2,       & ! Variance of w (overall)                         [m^2/s^2]
         712             :       mixt_frac, & ! Mixture fraction                                      [-]
         713             :       F_w          ! Parameter for the spread of the PDF comp. means of w  [-]
         714             :
         715             :     ! Input/Output Variable
         716             :     real( kind = core_rknd ), intent(inout) :: &
         717             :       Skx    ! Skewness of x (overall)              [-]
         718             :
         719             :     ! Output Variables
         720             :     real( kind = core_rknd ), intent(out) :: &
         721             :       mu_x_1,        & ! Mean of x (1st PDF component)        [units vary]
         722             :       mu_x_2,        & ! Mean of x (2nd PDF component)        [units vary]
         723             :       sigma_x_1_sqd, & ! Variance of x (1st PDF component)    [(units vary)^2]
         724             :       sigma_x_2_sqd    ! Variance of x (2nd PDF component)    [(units vary)^2]
         725             :
         726             :     real( kind = core_rknd ), intent(out) :: &
         727             :       coef_sigma_x_1_sqd, & ! sigma_x_1^2 = coef_sigma_x_1_sqd * <x'^2>    [-]
         728             :       coef_sigma_x_2_sqd    ! sigma_x_2^2 = coef_sigma_x_2_sqd * <x'^2>    [-]
         729             :
         730             :     ! Local Variables
         731             :     real( kind = core_rknd ) :: &
         732             :       corr_w_x, & ! Correlation (overall) of w and x                  [-]
         733             :       min_Skx,  & ! Minimum Skx that can be represented by the PDF    [-]
         734             :       max_Skx     ! Maximum Skx that can be represented by the PDF    [-]
         735             :
         736             :
         737             :     ! Calculate the limits of Skx.
         738           0 :     if ( F_w > zero ) then
         739             :
         740             :        ! Calculate the overall correlation of w and x.
         741           0 :        if ( wp2 * xp2 > zero ) then
         742             :
         743             :           ! The overall correlation of w and x is:
         744             :           !
         745             :           ! corr_w_x = <w'x'> / sqrt( <w'^2> * <x'^2> ).
         746           0 :           corr_w_x = wpxp / sqrt( wp2 * xp2 )
         747             :
         748             :        else ! <w'^2> * <x'^2> = 0
         749             :
         750             :           ! When <w'^2> = 0 or <x'^2> = 0, <w'x'> = 0.  The correlation of w
         751             :           ! and x is undefined, however, since <w'x'> = 0, Skx = 0.  Setting
         752             :           ! corr_w_x = 0 in this scenario will set max_Skx = min_Skx = 0.
         753             :           corr_w_x = zero
         754             :
         755             :        endif
         756             :
         757           0 :        if ( wpxp >= zero ) then
         758             :
         759             :           min_Skx = ( one + mixt_frac ) &
         760             :                     / sqrt( mixt_frac * ( one - mixt_frac ) ) &
         761             :                     * corr_w_x**3 / F_w**three_halves &
         762             :                     - sqrt( mixt_frac / ( one - mixt_frac ) ) &
         763           0 :                       * three * corr_w_x / sqrt( F_w )
         764             :
         765             :           max_Skx = ( mixt_frac - two ) &
         766             :                     / sqrt( mixt_frac * ( one - mixt_frac ) ) &
         767             :                     * corr_w_x**3 / F_w**three_halves &
         768             :                     + sqrt( ( one - mixt_frac ) / mixt_frac ) &
         769           0 :                       * three * corr_w_x / sqrt( F_w )
         770             :
         771             :        else ! <w'x'> < 0
         772             :
         773             :           min_Skx = ( mixt_frac - two ) &
         774             :                     / sqrt( mixt_frac * ( one - mixt_frac ) ) &
         775             :                     * corr_w_x**3 / F_w**three_halves &
         776             :                     + sqrt( ( one - mixt_frac ) / mixt_frac ) &
         777           0 :                       * three * corr_w_x / sqrt( F_w )
         778             :
         779             :           max_Skx = ( one + mixt_frac ) &
         780             :                     / sqrt( mixt_frac * ( one - mixt_frac ) ) &
         781             :                     * corr_w_x**3 / F_w**three_halves &
         782             :                     - sqrt( mixt_frac / ( one - mixt_frac ) ) &
         783           0 :                       * three * corr_w_x / sqrt( F_w )
         784             :
         785             :        endif ! <w'x'> >= 0
         786             :
         787             :     else ! F_w > 0
         788             :
         789             :        ! When F_w = 0, <w'x'> = 0, and Skx = 0.
         790             :        min_Skx = zero
         791             :        max_Skx = zero
         792             :
         793             :     endif ! F_w > 0
         794             :
         795             :     ! Limit Skx so that min_Skx <= Skx <= max_Skx.
         796           0 :     if ( Skx > max_Skx ) then
         797           0 :        Skx = max_Skx
         798           0 :     elseif ( Skx < min_Skx ) then
         799           0 :        Skx = min_Skx
         800             :     endif ! Skx >= max_Skx
         801             :
         802             :     ! Calculate the PDF parameters for responder variable x.
         803             :     call calculate_responder_params( xm, xp2, Skx, wpxp,  & ! In
         804             :                                      wp2, F_w, mixt_frac, & ! In
         805             :                                      mu_x_1, mu_x_2,      & ! Out
         806             :                                      sigma_x_1_sqd,       & ! Out
         807             :                                      sigma_x_2_sqd,       & ! Out
         808             :                                      coef_sigma_x_1_sqd,  & ! Out
         809           0 :                                      coef_sigma_x_2_sqd   ) ! Out
         810             :
         811             :
         812           0 :     return
         813             :
         814             :   end subroutine calc_responder_driver
         815             :
         816             :   !=============================================================================
         817           0 :   elemental subroutine calc_F_w_zeta_w( Skw, wprtp, wpthlp, upwp, vpwp, & ! In
         818             :                                         wp2, rtp2, thlp2, up2, vp2,     & ! In
         819             :                                         gamma_Skw_fnc,                  & ! In
         820             :                                         slope_coef_spread_DG_means_w,   & ! In
         821             :                                         pdf_component_stdev_factor_w,   & ! In
         822             :                                         max_corr_w_sclr_sqd,            & ! In
         823             :                                         F_w, zeta_w, min_F_w, max_F_w   ) ! Out
         824             :
         825             :     ! Description:
         826             :     ! Calculates the values of F_w and zeta_w for w, which is the setter
         827             :     ! variable (which is the variable that sets the mixture fraction).
         828             :     !
         829             :     ! Based purely on the PDF of w and not considering restrictions caused by
         830             :     ! other variables in the multivariate PDF, the value of F_w is calculated
         831             :     ! between 0 (min_F_w) and 1 (max_F_w).  However, other variables in the PDF
         832             :     ! place restrictions on the minimum value of F_w.  The range of acceptible
         833             :     ! values for F_w is given by:
         834             :     !
         835             :     ! max( <w'x'>^2 / ( <w'^2> * <x'^2> ), for all variables x ) <= F_w <= 1.
         836             :     !
         837             :     ! Additionally, the value of F_w should still increase with an increasing
         838             :     ! magnitude of Skw.  The value of F_w can be 0 only when Skw = 0 (as well as
         839             :     ! all <w'x'> = 0, as well).  The F_w skewness equation used is:
         840             :     !
         841             :     ! F_w = max_F_w + ( min_F_w - max_F_w )
         842             :     !                 * exp{ -|Skw|^lambda / slope_coef_spread_DG_means_w };
         843             :     !
         844             :     ! where lambda > 0 and slope_coef_spread_DG_means_w > 0.  As |Skw| goes
         845             :     ! toward 0, the value of F_w goes toward min_F_w, and as |Skw| becomes
         846             :     ! large, the value of F_w goes toward max_F_w.  When the value of
         847             :     ! slope_coef_spread_DG_means_w is small, the value of F_w tends toward
         848             :     ! max_F_w, and when the value of slope_coef_spread_DG_means_w is large, the
         849             :     ! value of F_w tends toward min_F_w.  When lambda is small, the value of F_w
         850             :     ! is less dependent on Skw, and when lambda is large, the value of F_w is
         851             :     ! more dependent on Skw.
         852             :     !
         853             :     ! Mathematically, this equation will always produce a value of F_w that
         854             :     ! falls between min_F_w and max_F_w.  However, in order to prevent a value
         855             :     ! of F_w from being calculated outside the bounds of min_F_w and max_F_w
         856             :     ! owing to numerical underflow or loss of precision, this equation can be
         857             :     ! rewritten as:
         858             :     !
         859             :     ! F_w
         860             :     ! = min_F_w * exp{ -|Skw|^lambda / slope_coef_spread_DG_means_w }
         861             :     !   + max_F_w * ( 1 - exp{ -|Skw|^lambda / slope_coef_spread_DG_means_w } ).
         862             :     !
         863             :     ! The value of zeta_w used to adjust the PDF component standard devations:
         864             :     !
         865             :     ! 1 + zeta_w = ( mixt_frac * sigma_w_1^2 )
         866             :     !              / ( ( 1 - mixt_frac ) * sigma_w_2^2 );
         867             :     !
         868             :     ! where zeta_w > -1.  The sign of zeta_w is used to easily determine if
         869             :     ! mixt_frac * sigma_w_1^2 is greater than ( 1 - mixt_frac ) * sigma_w_2^2
         870             :     ! (when zeta_w is positive), mixt_frac * sigma_w_1^2 is less than
         871             :     ! ( 1 - mixt_frac ) * sigma_w_2^2 (when zeta_w is negative), or
         872             :     ! mixt_frac * sigma_w_1^2 is equal to ( 1 - mixt_frac ) * sigma_w_2^2 (when
         873             :     ! zeta_w is 0).
         874             :     !
         875             :     ! In order to allow for a tunable parameter that is the pure ratio of
         876             :     ! mixt_frac * sigma_w_1^2 to ( 1 - mixt_frac ) * sigma_w_2^2, zeta_w is
         877             :     ! related to the parameter pdf_component_stdev_factor_w, where:
         878             :     !
         879             :     ! 1 + zeta_w = pdf_component_stdev_factor_w.
         880             :
         881             :     ! References:
         882             :     !-----------------------------------------------------------------------
         883             :
         884             :     use constants_clubb, only: &
         885             :         one,                      & ! Variable(s)
         886             :         zero,                     &
         887             :         max_mag_correlation_flux
         888             :
         889             :     use clubb_precision, only: &
         890             :         core_rknd    ! Variable(s)
         891             :
         892             :     implicit none
         893             :
         894             :     ! Input Variables
         895             :     real( kind = core_rknd ), intent(in) :: &
         896             :       Skw,    & ! Skewness of w (overall)              [-]
         897             :       wprtp,  & ! Covariance (overall) of w and rt     [m/s (kg/kg)]
         898             :       wpthlp, & ! Covariance (overall) of w and thl    [m/s K]
         899             :       upwp,   & ! Covariance (overall) of u and w      [m^2/s^2]
         900             :       vpwp,   & ! Covariance (overall) of u and v      [m^2/s^2]
         901             :       wp2,    & ! Variance (overall) of w              [m^2/s^2]
         902             :       rtp2,   & ! Variance (overall) of rt             [(kg/kg)^2]
         903             :       thlp2,  & ! Variance (overall) of thl            [K^2]
         904             :       up2,    & ! Variance (overall) of u              [m^2/s^2]
         905             :       vp2       ! Variance (overall) of v              [m^2/s^2]
         906             :
         907             :     ! Tunable parameter gamma.
         908             :     ! When gamma goes to 0, the standard deviations of w in each PDF component
         909             :     ! become small, and the spread between the two PDF component means of w
         910             :     ! becomes large.  F_w goes to min_F_w.
         911             :     ! When gamma goes to 1, the standard deviations of w in each PDF component
         912             :     ! become large, and the spread between the two PDF component means of w
         913             :     ! becomes small.  F_w goes to max_F_w.
         914             :     real( kind = core_rknd ), intent(in) :: &
         915             :       gamma_Skw_fnc    ! Value of parameter gamma from tunable Skw function  [-]
         916             :
         917             :     real( kind = core_rknd ), intent(in) :: &
         918             :       ! Slope coefficient for the spread between the PDF component means of w.
         919             :       slope_coef_spread_DG_means_w, &
         920             :       ! Parameter to adjust the PDF component standard deviations of w.
         921             :       pdf_component_stdev_factor_w
         922             :
         923             :     real( kind = core_rknd ), intent(in) :: &
         924             :       max_corr_w_sclr_sqd    ! Max value of wpsclrp^2 / ( wp2 * sclrp2 )     [-]
         925             :
         926             :     ! Output Variables
         927             :     real( kind = core_rknd ), intent(out) :: &
         928             :       F_w,     & ! Parameter for the spread of the PDF component means of w  [-]
         929             :       zeta_w,  & ! Parameter for the PDF component variances of w            [-]
         930             :       min_F_w, & ! Minimum allowable value of parameter F_w                  [-]
         931             :       max_F_w    ! Maximum allowable value of parameter F_w                  [-]
         932             :
         933             :     ! Local Variables
         934             :     real( kind = core_rknd ) :: &
         935             :       corr_w_rt_sqd,    & ! Value of wprtp^2 / ( wp2 * rtp2 )              [-]
         936             :       corr_w_thl_sqd,   & ! Value of wpthlp^2 / ( wp2 * thlp2 )            [-]
         937             :       corr_u_w_sqd,     & ! Value of upwp^2 / ( up2 * wp2 )                [-]
         938             :       corr_v_w_sqd,     & ! Value of vpwp^2 / ( vp2 * wp2 )                [-]
         939             :       max_corr_w_x_sqd    ! Max. val. of wpxp^2/(wp2*xp2) for all vars. x  [-]
         940             :
         941             :     real( kind = core_rknd ) :: &
         942             :       exp_Skw_interp_factor, & ! Function to interp. between min. and max.   [-]
         943             :       zeta_w_star
         944             :
         945             :     real( kind = core_rknd ), parameter :: &
         946             :       lambda = 0.5_core_rknd    ! Parameter for Skw dependence    [-]
         947             :
         948             :
         949             :     ! Find the maximum value of <w'x'>^2 / ( <w'^2> * <x'^2> ) for all
         950             :     ! variables x that are Double Gaussian PDF responder variables.  This
         951             :     ! includes rt and theta-l.  When l_predict_upwp_vpwp is enabled, u and v are
         952             :     ! also calculated as part of the PDF, and they are included as well.
         953             :     ! Additionally, when sclr_dim > 0, passive scalars (sclr) are also included.
         954             :
         955             :     ! Calculate the squared value of the correlation of w and rt.
         956           0 :     if ( wp2 * rtp2 > zero ) then
         957           0 :        corr_w_rt_sqd = wprtp**2 / ( wp2 * rtp2 )
         958             :     else
         959             :        corr_w_rt_sqd = zero
         960             :     endif
         961             :
         962             :     ! Calculate the squared value of the correlation of w and thl.
         963           0 :     if ( wp2 * thlp2 > zero ) then
         964           0 :        corr_w_thl_sqd = wpthlp**2 / ( wp2 * thlp2 )
         965             :     else
         966             :        corr_w_thl_sqd = zero
         967             :     endif
         968             :
         969             :     ! Calculate the squared value of the correlation of u and w.
         970           0 :     if ( up2 * wp2 > zero ) then
         971           0 :        corr_u_w_sqd = upwp**2 / ( up2 * wp2 )
         972             :     else
         973             :        corr_u_w_sqd = zero
         974             :     endif
         975             :
         976             :     ! Calculate the squared value of the correlation of v and w.
         977           0 :     if ( vp2 * wp2 > zero ) then
         978           0 :        corr_v_w_sqd = vpwp**2 / ( vp2 * wp2 )
         979             :     else
         980             :        corr_v_w_sqd = zero
         981             :     endif
         982             :
         983             :     ! Calculate the maximum value of corr_w_rt^2, corr_w_thl^2, corr_u_w^2, and
         984             :     ! corr_v_w^2 in the calculation of max(corr_w_x^2).  Include the maximum
         985             :     ! value of corr_w_sclr^2 (for all sclrs) when scalars are found.
         986             :     max_corr_w_x_sqd = max( corr_w_rt_sqd, corr_w_thl_sqd, &
         987           0 :                             corr_u_w_sqd, corr_v_w_sqd, max_corr_w_sclr_sqd )
         988             :
         989           0 :     max_corr_w_x_sqd = min( max_corr_w_x_sqd, max_mag_correlation_flux**2 )
         990             :
         991             :     ! Calculate min_F_w and set max_F_w to 1.
         992           0 :     if ( abs( Skw ) > zero ) then
         993           0 :        min_F_w = max( max_corr_w_x_sqd, 1.0e-3_core_rknd )
         994             :     else
         995           0 :        min_F_w = max_corr_w_x_sqd
         996             :     endif
         997           0 :     max_F_w = one
         998             :
         999             :     ! F_w must have a value between min_F_w and max_F_w.
        1000             :     exp_Skw_interp_factor &
        1001             :     = exp( -abs(Skw)**lambda / slope_coef_spread_DG_means_w )
        1002             :
        1003             : !    F_w = min_F_w * exp_Skw_interp_factor &
        1004             : !          + max_F_w * ( one - exp_Skw_interp_factor )
        1005             :
        1006             :     ! For now, use a formulation similar to what is used for ADG1.
        1007             :     ! Tunable parameter gamma.
        1008             :     ! When gamma goes to 0, the standard deviations of w in each PDF component
        1009             :     ! become small, and the spread between the two PDF component means of w
        1010             :     ! becomes large.  F_w goes to min_F_w.
        1011             :     ! When gamma goes to 1, the standard deviations of w in each PDF component
        1012             :     ! become large, and the spread between the two PDF component means of w
        1013             :     ! becomes small.  F_w goes to max_F_w.
        1014           0 :     F_w = max_F_w - gamma_Skw_fnc * ( max_F_w - min_F_w )
        1015             :
        1016             :     ! The value of zeta_w must be greater than -1.
        1017           0 :     zeta_w_star = pdf_component_stdev_factor_w - one
        1018             :
        1019             :     ! Make the PDF of w symmetric.  In other words, the PDF at a value of
        1020             :     ! positive skewness will look like a mirror image of the PDF at the
        1021             :     ! opposite value of negative skewness.
        1022           0 :     if ( Skw >= zero ) then
        1023           0 :        zeta_w = zeta_w_star
        1024             :     else ! Skw < 0
        1025           0 :        zeta_w = - zeta_w_star / ( zeta_w_star + one )
        1026             :     endif
        1027             :
        1028             :
        1029           0 :     return
        1030             :
        1031             :   end subroutine calc_F_w_zeta_w
        1032             :
        1033             :   !=============================================================================
        1034             :
        1035             : end module new_hybrid_pdf_main

Generated by: LCOV version 1.14