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