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

          Line data    Source code
       1             : ! $Id$
       2             : !===============================================================================
       3             : module new_pdf_main
       4             : 
       5             :   ! Description:
       6             :   ! The portion of CLUBB's multivariate, two-component PDF that is the
       7             :   ! trivariate, two-component normal PDF of vertical velocity (w), total water
       8             :   ! mixing ratio (rt), and liquid water potential temperature (thl).
       9             : 
      10             :   ! References:
      11             :   !-------------------------------------------------------------------------
      12             : 
      13             :   implicit none
      14             : 
      15             :   public :: new_pdf_driver    ! Procedure(s)
      16             : 
      17             :   private :: calc_responder_var,     & ! Procedure(s)
      18             :              calc_F_x_zeta_x_setter, &
      19             :              calc_F_x_responder
      20             : 
      21             :   private
      22             : 
      23             :   contains
      24             : 
      25             :   !=============================================================================
      26           0 :   subroutine new_pdf_driver( nz, ngrdcol, wm, rtm, thlm, wp2, rtp2, thlp2, Skw,  & ! In
      27           0 :                              wprtp, wpthlp, rtpthlp,                    & ! In
      28           0 :                              clubb_params,                              & ! In
      29           0 :                              Skrt, Skthl,                               & ! I/O
      30           0 :                              mu_w_1, mu_w_2,                            & ! Out
      31           0 :                              mu_rt_1, mu_rt_2,                          & ! Out
      32           0 :                              mu_thl_1, mu_thl_2,                        & ! Out
      33           0 :                              sigma_w_1_sqd, sigma_w_2_sqd,              & ! Out
      34           0 :                              sigma_rt_1_sqd, sigma_rt_2_sqd,            & ! Out
      35           0 :                              sigma_thl_1_sqd, sigma_thl_2_sqd,          & ! Out
      36           0 :                              mixt_frac,                                 & ! Out
      37             :                              pdf_implicit_coefs_terms,                  & ! Out
      38           0 :                              F_w, F_rt, F_thl, min_F_w, max_F_w,        & ! Out
      39           0 :                              min_F_rt, max_F_rt, min_F_thl, max_F_thl )   ! Out
      40             :                              
      41             : 
      42             :     ! Description:
      43             :     ! Selects which variable is used to set the mixture fraction for the PDF
      44             :     ! ("the setter") and which variables are handled after that mixture fraction
      45             :     ! has been set ("the responders").  Traditionally, w has been used to set
      46             :     ! the PDF.
      47             : 
      48             :     ! References:
      49             :     !-----------------------------------------------------------------------
      50             : 
      51             :     use constants_clubb, only: &
      52             :         four,                & ! Variable(s)
      53             :         two,                 &
      54             :         one,                 &
      55             :         zero,                &
      56             :         rt_tol,              &
      57             :         thl_tol,             &
      58             :         max_mag_correlation
      59             : 
      60             :     use new_pdf, only: &
      61             :         calc_setter_var_params,     & ! Procedure(s)
      62             :         calc_coef_wp4_implicit,     &
      63             :         calc_coef_wpxp2_implicit,   &
      64             :         calc_coefs_wp2xp_semiimpl,  &
      65             :         calc_coefs_wpxpyp_semiimpl
      66             : 
      67             :     use pdf_parameter_module, only: &
      68             :         implicit_coefs_terms    ! Variable Type
      69             : 
      70             :     use model_flags, only: &
      71             :         l_explicit_turbulent_adv_wp3,  & ! Variable(s)
      72             :         l_explicit_turbulent_adv_wpxp, &
      73             :         l_explicit_turbulent_adv_xpyp
      74             : 
      75             :     use clubb_precision, only: &
      76             :         core_rknd    ! Variable(s)
      77             : 
      78             :      use parameter_indices, only: &
      79             :         nparams,                       & ! Variable(s)
      80             :         islope_coef_spread_DG_means_w, &
      81             :         ipdf_component_stdev_factor_w, &
      82             :         icoef_spread_DG_means_rt,      &
      83             :         icoef_spread_DG_means_thl
      84             : 
      85             :     implicit none
      86             : 
      87             :     integer, intent(in) :: &
      88             :       nz, &
      89             :       ngrdcol
      90             : 
      91             :     ! Input Variables
      92             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
      93             :       wm,      & ! Mean of w (overall)                 [m/s]
      94             :       rtm,     & ! Mean of rt (overall)                [kg/kg]
      95             :       thlm,    & ! Mean of thl (overall)               [K]
      96             :       wp2,     & ! Variance of w (overall)             [m^2/s^2]
      97             :       rtp2,    & ! Variance of rt (overall)            [kg^2/kg^2]
      98             :       thlp2,   & ! Variance of thl (overall)           [K^2]
      99             :       Skw,     & ! Skewness of w (overall)             [-]
     100             :       wprtp,   & ! Covariance of w and rt (overall)    [(m/s)kg/kg]
     101             :       wpthlp,  & ! Covariance of w and thl (overall)   [(m/s)K]
     102             :       rtpthlp    ! Covariance of rt and thl (overall)  [(kg/kg)K]
     103             : 
     104             :     real( kind = core_rknd ), dimension(ngrdcol,nparams), intent(in) :: &
     105             :       clubb_params    ! Array of CLUBB's tunable parameters    [units vary]
     106             : 
     107             :     ! Input/Output Variables
     108             :     ! These variables are input/output because their values may be clipped.
     109             :     ! Otherwise, as long as it is not necessary to clip them, their values
     110             :     ! will stay the same.
     111             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(inout) :: &
     112             :       Skrt,  & ! Skewness of rt (overall)            [-]
     113             :       Skthl    ! Skewness of thl (overall)           [-]
     114             : 
     115             :     ! Output Variables
     116             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
     117             :       mu_w_1,          & ! Mean of w (1st PDF component)        [m/s]
     118             :       mu_w_2,          & ! Mean of w (2nd PDF component)        [m/s]
     119             :       mu_rt_1,         & ! Mean of rt (1st PDF component)       [kg/kg]
     120             :       mu_rt_2,         & ! Mean of rt (2nd PDF component)       [kg/kg]
     121             :       mu_thl_1,        & ! Mean of thl (1st PDF component)      [K]
     122             :       mu_thl_2,        & ! Mean of thl (2nd PDF component)      [K]
     123             :       sigma_w_1_sqd,   & ! Variance of w (1st PDF component)    [m^2/s^2]
     124             :       sigma_w_2_sqd,   & ! Variance of w (2nd PDF component)    [m^2/s^2]
     125             :       sigma_rt_1_sqd,  & ! Variance of rt (1st PDF component)   [kg^2/kg^2]
     126             :       sigma_rt_2_sqd,  & ! Variance of rt (2nd PDF component)   [kg^2/kg^2]
     127             :       sigma_thl_1_sqd, & ! Variance of thl (1st PDF component)  [K^2]
     128             :       sigma_thl_2_sqd, & ! Variance of thl (2nd PDF component)  [K^2]
     129             :       mixt_frac          ! Mixture fraction                     [-]
     130             : 
     131             :     type(implicit_coefs_terms), intent(inout) :: &
     132             :       pdf_implicit_coefs_terms    ! Implicit coefs / explicit terms [units vary]
     133             : 
     134             :     ! Output only for recording statistics.
     135             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
     136             :       F_w,   & ! Parameter for the spread of the PDF component means of w    [-]
     137             :       F_rt,  & ! Parameter for the spread of the PDF component means of rt   [-]
     138             :       F_thl    ! Parameter for the spread of the PDF component means of thl  [-]
     139             : 
     140             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
     141             :       min_F_w,   & ! Minimum allowable value of parameter F_w      [-]
     142             :       max_F_w,   & ! Maximum allowable value of parameter F_w      [-]
     143             :       min_F_rt,  & ! Minimum allowable value of parameter F_rt     [-]
     144             :       max_F_rt,  & ! Maximum allowable value of parameter F_rt     [-]
     145             :       min_F_thl, & ! Minimum allowable value of parameter F_thl    [-]
     146             :       max_F_thl    ! Maximum allowable value of parameter F_thl    [-]
     147             : 
     148             :     ! Local Variables
     149             :     real( kind = core_rknd ), dimension(nz) :: &
     150           0 :       sigma_w_1,   & ! Standard deviation of w (1st PDF component)      [m/s]
     151           0 :       sigma_w_2,   & ! Standard deviation of w (2nd PDF component)      [m/s]
     152           0 :       sgn_wprtp,   & ! Sign of the covariance of w and rt (overall)     [-]
     153           0 :       sgn_wpthlp,  & ! Sign of the covariance of w and thl (overall)    [-]
     154           0 :       sgn_wp2        ! Sign of the variance of w (overall); always pos. [-]
     155             : 
     156             :     real( kind = core_rknd ), dimension(nz) :: &
     157           0 :       coef_sigma_w_1_sqd,   & ! sigma_w_1^2 = coef_sigma_w_1_sqd * <w'^2>    [-]
     158           0 :       coef_sigma_w_2_sqd,   & ! sigma_w_2^2 = coef_sigma_w_2_sqd * <w'^2>    [-]
     159           0 :       coef_sigma_rt_1_sqd,  & ! sigma_rt_1^2 = coef_sigma_rt_1_sqd * <rt'^2> [-]
     160           0 :       coef_sigma_rt_2_sqd,  & ! sigma_rt_2^2 = coef_sigma_rt_2_sqd * <rt'^2> [-]
     161           0 :       coef_sigma_thl_1_sqd, & ! sigma_thl_1^2=coef_sigma_thl_1_sqd*<thl'^2>  [-]
     162           0 :       coef_sigma_thl_2_sqd    ! sigma_thl_2^2=coef_sigma_thl_2_sqd*<thl'^2>  [-]
     163             : 
     164             :     real( kind = core_rknd ), dimension(nz) :: &
     165           0 :       max_Skx2_pos_Skx_sgn_wpxp, & ! Maximum Skx^2 when Skx*sgn(<w'x'>) >= 0 [-]
     166           0 :       max_Skx2_neg_Skx_sgn_wpxp    ! Maximum Skx^2 when Skx*sgn(<w'x'>) < 0  [-]
     167             : 
     168             :     real( kind = core_rknd ), dimension(nz) :: &
     169           0 :       zeta_w    ! Parameter for the PDF component variances of w           [-]
     170             : 
     171             :     real ( kind = core_rknd ) :: &
     172             :       lambda_w    ! Param. that increases or decreases Skw dependence  [-]
     173             : 
     174             :     real ( kind = core_rknd ), dimension(nz) :: &
     175           0 :       exp_factor_rt,   & ! Factor of the form 1 - exp{} that reduces F_rt   [-]
     176           0 :       exp_factor_thl,  & ! Don't reduce F_thl by exp_factor_thl        [-]
     177           0 :       adj_corr_rt_thl    ! Adjusted (overall) correlation of rt and theta-l [-]
     178             : 
     179             :     real ( kind = core_rknd ), dimension(nz) :: &
     180           0 :       coef_wp4_implicit,     & ! <w'^4> = coef_wp4_implicit * <w'^2>^2       [-]
     181           0 :       coef_wprtp2_implicit,  & ! <w'rt'^2> = coef_wprtp2_implicit*<rt'^2>  [m/s]
     182           0 :       coef_wpthlp2_implicit    ! <w'thl'^2>=coef_wpthlp2_implicit*<thl'^2> [m/s]
     183             : 
     184             :     ! <w'^2 rt'> = coef_wp2rtp_implicit * <w'rt'> + term_wp2rtp_explicit
     185             :     real ( kind = core_rknd ), dimension(nz) :: &
     186           0 :       coef_wp2rtp_implicit, & ! Coefficient that is multiplied by <w'rt'>  [m/s]
     187           0 :       term_wp2rtp_explicit    ! Term that is on the RHS          [m^2/s^2 kg/kg]
     188             : 
     189             :     ! <w'^2 thl'> = coef_wp2thlp_implicit * <w'thl'> + term_wp2thlp_explicit
     190             :     real ( kind = core_rknd ), dimension(nz) :: &
     191           0 :       coef_wp2thlp_implicit, & ! Coef. that is multiplied by <w'thl'>      [m/s]
     192           0 :       term_wp2thlp_explicit    ! Term that is on the RHS             [m^2/s^2 K]
     193             : 
     194             :     ! <w'rt'thl'> = coef_wprtpthlp_implicit*<rt'thl'> + term_wprtpthlp_explicit
     195             :     real ( kind = core_rknd ), dimension(nz) :: &
     196           0 :       coef_wprtpthlp_implicit, & ! Coef. that is multiplied by <rt'thl'>   [m/s]
     197           0 :       term_wprtpthlp_explicit    ! Term that is on the RHS         [m/s(kg/kg)K]
     198             :       
     199             :     integer :: &
     200             :       i
     201             :       
     202             :     ! ------------------------- Begin Code -------------------------
     203             :       
     204           0 :     do  i = 1, ngrdcol
     205             : 
     206             :       ! Calculate sgn( <w'rt'> ).
     207           0 :       where ( wprtp(i,:) >= zero )
     208             :          sgn_wprtp = one
     209             :       elsewhere ! wprtp < 0
     210             :          sgn_wprtp = -one
     211             :       endwhere ! wprtp >= 0
     212             : 
     213             :       ! Calculate sgn( <w'thl'> ).
     214           0 :       where ( wpthlp(i,:) >= zero )
     215             :          sgn_wpthlp = one
     216             :       elsewhere ! wpthlp < 0
     217             :          sgn_wpthlp = -one
     218             :       endwhere ! wpthlp >= 0
     219             : 
     220             :       ! Sign of the variance of w (overall), which is always positive.
     221           0 :       sgn_wp2 = one
     222             : 
     223           0 :       lambda_w = 0.5_core_rknd
     224             : 
     225             :       ! Calculate the adjusted (overall) correlation of rt and theta-l, and the
     226             :       ! value of exp_factor_rt.
     227           0 :       where ( rtp2(i,:) >= rt_tol**2 .and. thlp2(i,:) >= thl_tol**2 )
     228             :          adj_corr_rt_thl = rtpthlp(i,:) / sqrt( rtp2(i,:) * thlp2(i,:) ) * sgn_wprtp * sgn_wpthlp
     229             :          adj_corr_rt_thl = min( max( adj_corr_rt_thl, -max_mag_correlation ), &
     230             :                                 max_mag_correlation )
     231             :          exp_factor_rt = one &
     232             :                          - exp( -0.2_core_rknd * ( adj_corr_rt_thl + one )**5 )
     233             :       elsewhere ! <rt'^2> < rt_tol^2 or <thl'^2> < thl_tol^2
     234             :          adj_corr_rt_thl = zero  ! adj_corr_rt_thl is undefined in this scenario.
     235             :          exp_factor_rt = one     ! Set exp_factor_rt to 1.
     236             :       endwhere ! <rt'^2> >= rt_tol^2 and <thl'^2> >= thl_tol^2
     237             : 
     238             :       ! The value of F_thl is not reduced by exp_factor_thl.
     239           0 :       exp_factor_thl = one
     240             : 
     241             : 
     242             :       ! Vertical velocity, w, will always be the setter variable.
     243             :       call calc_F_x_zeta_x_setter( nz, Skw(i,:),                      & ! In
     244             :                                    clubb_params(i,islope_coef_spread_DG_means_w), & ! In
     245             :                                    clubb_params(i,ipdf_component_stdev_factor_w), & ! In
     246             :                                    lambda_w,                     & ! In
     247             :                                    F_w(i,:), zeta_w,                  & ! Out
     248           0 :                                    min_F_w(i,:), max_F_w(i,:)              ) ! Out
     249             : 
     250             :       ! Calculate the PDF parameters, including mixture fraction, for the
     251             :       ! setter variable, w.
     252             :       call calc_setter_var_params( nz, wm(i,:), wp2(i,:), Skw(i,:), sgn_wp2, & ! In
     253             :                                    F_w(i,:), zeta_w,               & ! In
     254             :                                    mu_w_1(i,:), mu_w_2(i,:), sigma_w_1, & ! Out
     255             :                                    sigma_w_2, mixt_frac(i,:),      & ! Out
     256             :                                    coef_sigma_w_1_sqd,        & ! Out
     257           0 :                                    coef_sigma_w_2_sqd         ) ! Out
     258             : 
     259           0 :       sigma_w_1_sqd(i,:) = sigma_w_1**2
     260           0 :       sigma_w_2_sqd(i,:) = sigma_w_2**2
     261             : 
     262             :       ! Calculate the upper limit on the magnitude of skewness for responding
     263             :       ! variables.
     264             :       max_Skx2_pos_Skx_sgn_wpxp = four * ( one - mixt_frac(i,:) )**2 &
     265           0 :                                   / ( mixt_frac(i,:) * ( two - mixt_frac(i,:) ) )
     266             : 
     267           0 :       max_Skx2_neg_Skx_sgn_wpxp = four * mixt_frac(i,:)**2 / ( one - mixt_frac(i,:)**2 )
     268             : 
     269             :       ! Calculate the PDF parameters for responder variable rt.
     270             :       call calc_responder_var( nz, rtm(i,:), rtp2(i,:), sgn_wprtp, mixt_frac(i,:), & ! In
     271             :                                clubb_params(i,icoef_spread_DG_means_rt),         & ! In
     272             :                                exp_factor_rt,                   & ! In
     273             :                                max_Skx2_pos_Skx_sgn_wpxp,       & ! In
     274             :                                max_Skx2_neg_Skx_sgn_wpxp,       & ! In
     275             :                                Skrt(i,:),                            & ! In/Out
     276             :                                mu_rt_1(i,:), mu_rt_2(i,:),                & ! Out
     277             :                                sigma_rt_1_sqd(i,:), sigma_rt_2_sqd(i,:),  & ! Out
     278             :                                coef_sigma_rt_1_sqd,             & ! Out
     279             :                                coef_sigma_rt_2_sqd,             & ! Out
     280           0 :                                F_rt(i,:), min_F_rt(i,:), max_F_rt(i,:)         ) ! Out
     281             : 
     282             :       ! Calculate the PDF parameters for responder variable thl.
     283             :       call calc_responder_var( nz, thlm(i,:), thlp2(i,:), sgn_wpthlp, mixt_frac(i,:), & ! In
     284             :                                clubb_params(i,icoef_spread_DG_means_thl),           & ! In
     285             :                                exp_factor_thl,                     & ! In
     286             :                                max_Skx2_pos_Skx_sgn_wpxp,          & ! In
     287             :                                max_Skx2_neg_Skx_sgn_wpxp,          & ! In
     288             :                                Skthl(i,:),                              & ! In/Out
     289             :                                mu_thl_1(i,:), mu_thl_2(i,:),                 & ! Out
     290             :                                sigma_thl_1_sqd(i,:), sigma_thl_2_sqd(i,:),   & ! Out
     291             :                                coef_sigma_thl_1_sqd,               & ! Out
     292             :                                coef_sigma_thl_2_sqd,               & ! Out
     293           0 :                                F_thl(i,:), min_F_thl(i,:), max_F_thl(i,:)         ) ! Out
     294             : 
     295             : 
     296             :       if ( .not. l_explicit_turbulent_adv_wp3 ) then
     297             : 
     298             :          ! Turbulent advection of <w'^3> is being handled implicitly.
     299             : 
     300             :          ! <w'^4> = coef_wp4_implicit * <w'^2>^2.
     301             :          coef_wp4_implicit &
     302             :          = calc_coef_wp4_implicit( nz, mixt_frac(i,:), F_w(i,:), &
     303             :                                    coef_sigma_w_1_sqd, &
     304           0 :                                    coef_sigma_w_2_sqd )
     305             : 
     306             :       else ! l_explicit_turbulent_adv_wp3
     307             : 
     308             :          ! Turbulent advection of <w'^3> is being handled explicitly.
     309             :          coef_wp4_implicit = zero
     310             : 
     311             :       endif ! .not. l_explicit_turbulent_adv_wp3
     312             : 
     313             :       if ( .not. l_explicit_turbulent_adv_xpyp ) then
     314             : 
     315             :          ! Turbulent advection of <rt'^2> and <thl'^2> is being handled
     316             :          ! implicitly.  Turbulent advection of <rt'thl'> is being handled
     317             :          ! semi-implicitly.
     318             : 
     319             :          ! <w'rt'^2> = coef_wprtp2_implicit * <rt'^2>
     320             :          coef_wprtp2_implicit &
     321             :          = calc_coef_wpxp2_implicit( nz, wp2(i,:), rtp2(i,:), wprtp(i,:), sgn_wprtp, &
     322             :                                      mixt_frac(i,:), F_w(i,:), F_rt(i,:), &
     323             :                                      coef_sigma_w_1_sqd, &
     324             :                                      coef_sigma_w_2_sqd, &
     325             :                                      coef_sigma_rt_1_sqd, &
     326           0 :                                      coef_sigma_rt_2_sqd  )
     327             : 
     328             :          ! <w'thl'^2> = coef_wpthlp2_implicit * <thl'^2>
     329             :          coef_wpthlp2_implicit &
     330             :          = calc_coef_wpxp2_implicit( nz, wp2(i,:), thlp2(i,:), wpthlp(i,:), sgn_wpthlp, &
     331             :                                      mixt_frac(i,:), F_w(i,:), F_thl(i,:), &
     332             :                                      coef_sigma_w_1_sqd, &
     333             :                                      coef_sigma_w_2_sqd, &
     334             :                                      coef_sigma_thl_1_sqd, &
     335           0 :                                      coef_sigma_thl_2_sqd  )
     336             : 
     337             :          ! <w'rt'thl'> = coef_wprtpthlp_implicit * <rt'thl'>
     338             :          !               + term_wprtpthlp_explicit
     339             :          call calc_coefs_wpxpyp_semiimpl( nz, wp2(i,:), rtp2(i,:), thlp2(i,:), wprtp(i,:),       & ! In
     340             :                                           wpthlp(i,:), sgn_wprtp, sgn_wpthlp, & ! In
     341             :                                           mixt_frac(i,:), F_w(i,:), F_rt(i,:), F_thl(i,:),   & ! In
     342             :                                           coef_sigma_w_1_sqd  ,          & ! In
     343             :                                           coef_sigma_w_2_sqd,            & ! In
     344             :                                           coef_sigma_rt_1_sqd,           & ! In
     345             :                                           coef_sigma_rt_2_sqd,           & ! In
     346             :                                           coef_sigma_thl_1_sqd,          & ! In
     347             :                                           coef_sigma_thl_2_sqd,          & ! In
     348             :                                           coef_wprtpthlp_implicit,       & ! Out
     349           0 :                                           term_wprtpthlp_explicit        ) ! Out
     350             : 
     351             :       else ! l_explicit_turbulent_adv_xpyp
     352             : 
     353             :          ! Turbulent advection of <rt'^2>, <thl'^2>, and <rt'thl'> is being
     354             :          ! handled explicitly.
     355             :          coef_wprtp2_implicit = zero
     356             :          coef_wpthlp2_implicit = zero
     357             :          coef_wprtpthlp_implicit = zero
     358             :          term_wprtpthlp_explicit = zero
     359             : 
     360             :       endif ! .not. l_explicit_turbulent_adv_xpyp
     361             : 
     362             :       if ( .not. l_explicit_turbulent_adv_wpxp ) then
     363             : 
     364             :          ! Turbulent advection of <w'rt'> and <w'thl'> is being handled
     365             :          ! semi-implicitly.
     366             : 
     367             :          ! <w'^2 rt'> = coef_wp2rtp_implicit * <w'rt'> + term_wp2rtp_explicit
     368             :          call calc_coefs_wp2xp_semiimpl( nz, wp2(i,:), rtp2(i,:), sgn_wprtp, & ! In
     369             :                                          mixt_frac(i,:), F_w(i,:), F_rt(i,:), & ! In
     370             :                                          coef_sigma_w_1_sqd,   & ! In
     371             :                                          coef_sigma_w_2_sqd,   & ! In
     372             :                                          coef_sigma_rt_1_sqd,  & ! In
     373             :                                          coef_sigma_rt_2_sqd,  & ! In
     374             :                                          coef_wp2rtp_implicit, & ! Out
     375           0 :                                          term_wp2rtp_explicit  ) ! Out
     376             : 
     377             :          ! <w'^2 thl'> = coef_wp2thlp_implicit * <w'thl'> + term_wp2thlp_explicit
     378             :          call calc_coefs_wp2xp_semiimpl( nz, wp2(i,:), thlp2(i,:), sgn_wpthlp, & ! In
     379             :                                          mixt_frac(i,:), F_w(i,:), F_thl(i,:),  & ! In
     380             :                                          coef_sigma_w_1_sqd,     & ! In
     381             :                                          coef_sigma_w_2_sqd,     & ! In
     382             :                                          coef_sigma_thl_1_sqd,   & ! In
     383             :                                          coef_sigma_thl_2_sqd,   & ! In
     384             :                                          coef_wp2thlp_implicit,  & ! Out
     385           0 :                                          term_wp2thlp_explicit   ) ! Out
     386             : 
     387             :       else ! l_explicit_turbulent_adv_wpxp
     388             : 
     389             :          ! Turbulent advection of <w'rt'> and <w'thl'> is being handled
     390             :          ! explicitly.
     391             :          coef_wp2rtp_implicit = zero
     392             :          term_wp2rtp_explicit = zero
     393             :          coef_wp2thlp_implicit = zero
     394             :          term_wp2thlp_explicit = zero
     395             : 
     396             :       endif ! .not. l_explicit_turbulent_adv_wpxp
     397             : 
     398             :       ! Pack the implicit coefficients and explicit terms into a single type
     399             :       ! variable for output.
     400           0 :       pdf_implicit_coefs_terms%coef_wp4_implicit(i,:) = coef_wp4_implicit
     401           0 :       pdf_implicit_coefs_terms%coef_wprtp2_implicit(i,:) = coef_wprtp2_implicit
     402           0 :       pdf_implicit_coefs_terms%coef_wpthlp2_implicit(i,:) = coef_wpthlp2_implicit
     403           0 :       pdf_implicit_coefs_terms%coef_wp2rtp_implicit(i,:) = coef_wp2rtp_implicit
     404           0 :       pdf_implicit_coefs_terms%term_wp2rtp_explicit(i,:) = term_wp2rtp_explicit
     405           0 :       pdf_implicit_coefs_terms%coef_wp2thlp_implicit(i,:) = coef_wp2thlp_implicit
     406           0 :       pdf_implicit_coefs_terms%term_wp2thlp_explicit(i,:) = term_wp2thlp_explicit
     407           0 :       pdf_implicit_coefs_terms%coef_wprtpthlp_implicit(i,:) = coef_wprtpthlp_implicit
     408           0 :       pdf_implicit_coefs_terms%term_wprtpthlp_explicit(i,:) = term_wprtpthlp_explicit
     409             : 
     410             :     end do
     411             : 
     412           0 :     return
     413             : 
     414             :   end subroutine new_pdf_driver
     415             : 
     416             :   !=============================================================================
     417           0 :   subroutine calc_responder_var( nz, xm, xp2, sgn_wpxp, mixt_frac, & ! In
     418             :                                  coef_spread_DG_means_x,       & ! In
     419           0 :                                  exp_factor_x,                 & ! In
     420           0 :                                  max_Skx2_pos_Skx_sgn_wpxp,    & ! In
     421           0 :                                  max_Skx2_neg_Skx_sgn_wpxp,    & ! In
     422           0 :                                  Skx,                          & ! In/Out
     423           0 :                                  mu_x_1, mu_x_2,               & ! Out
     424           0 :                                  sigma_x_1_sqd, sigma_x_2_sqd, & ! Out
     425           0 :                                  coef_sigma_x_1_sqd,           & ! Out
     426           0 :                                  coef_sigma_x_2_sqd,           & ! Out
     427           0 :                                  F_x, min_F_x, max_F_x         ) ! Out
     428             : 
     429             :     ! Description:
     430             :     ! This is the sub-driver for a responder variable.  The upper limits of the
     431             :     ! magnitude of Skx are calculated, and Skx is clipped when its magnitude
     432             :     ! exceeds the upper limits.  The limits of the F_x parameter are calculated,
     433             :     ! and the value of F_x is set within those limits.  Then, the PDF parameters
     434             :     ! for responder variable x are calculated.
     435             : 
     436             :     ! References:
     437             :     !-----------------------------------------------------------------------
     438             : 
     439             :     use constants_clubb, only: &
     440             :         zero    ! Variable(s)
     441             : 
     442             :     use new_pdf, only: &
     443             :         calc_limits_F_x_responder, & ! Procedure(s)
     444             :         calc_responder_params
     445             : 
     446             :     use clubb_precision, only: &
     447             :         core_rknd    ! Variable(s)
     448             : 
     449             :     implicit none
     450             : 
     451             :     integer, intent(in) :: &
     452             :       nz
     453             : 
     454             :     ! Input Variables
     455             :     real( kind = core_rknd ), dimension(nz), intent(in) :: &
     456             :       xm,           & ! Mean of x (overall)                       [units vary]
     457             :       xp2,          & ! Variance of x (overall)               [(units vary)^2]
     458             :       sgn_wpxp,     & ! Sign of the covariance of w and x                  [-]
     459             :       mixt_frac,    & ! Mixture fraction                                   [-]
     460             :       exp_factor_x    ! Factor of the form 1 - exp{}; reduces F_x          [-]
     461             : 
     462             :     real( kind = core_rknd ), intent(in) :: &
     463             :       coef_spread_DG_means_x    ! Coef.: spread betw. PDF comp. means of x   [-]
     464             : 
     465             :     real( kind = core_rknd ), dimension(nz), intent(in) :: &
     466             :       max_Skx2_pos_Skx_sgn_wpxp, & ! Maximum Skx^2 when Skx*sgn(<w'x'>) >= 0 [-]
     467             :       max_Skx2_neg_Skx_sgn_wpxp    ! Maximum Skx^2 when Skx*sgn(<w'x'>) < 0  [-]
     468             : 
     469             :     ! Input/Output Variable
     470             :     real( kind = core_rknd ), dimension(nz), intent(inout) :: &
     471             :       Skx    ! Skewness of x (overall)              [-]
     472             : 
     473             :     ! Output Variables
     474             :     real( kind = core_rknd ), dimension(nz), intent(out) :: &
     475             :       mu_x_1,        & ! Mean of x (1st PDF component)        [units vary]
     476             :       mu_x_2,        & ! Mean of x (2nd PDF component)        [units vary]
     477             :       sigma_x_1_sqd, & ! Variance of x (1st PDF component)    [(units vary)^2]
     478             :       sigma_x_2_sqd    ! Variance of x (2nd PDF component)    [(units vary)^2]
     479             : 
     480             :     real( kind = core_rknd ), dimension(nz), intent(out) :: &
     481             :       coef_sigma_x_1_sqd, & ! sigma_x_1^2 = coef_sigma_x_1_sqd * <x'^2>    [-]
     482             :       coef_sigma_x_2_sqd    ! sigma_x_2^2 = coef_sigma_x_2_sqd * <x'^2>    [-]
     483             : 
     484             :     ! Output only for recording statistics.
     485             :     real( kind = core_rknd ), dimension(nz), intent(out) :: &
     486             :       F_x,     & ! Param. for the spread betw. the PDF component means of x  [-]
     487             :       min_F_x, & ! Minimum allowable value of parameter F_x                  [-]
     488             :       max_F_x    ! Maximum allowable value of parameter F_x                  [-]
     489             : 
     490             : 
     491             :     ! Calculate the upper limit of the magnitude of Skx.
     492           0 :     where ( Skx * sgn_wpxp >= zero )
     493             :        where ( Skx**2 >= max_Skx2_pos_Skx_sgn_wpxp )
     494             :           where ( Skx >= zero )
     495             :              Skx = sqrt( 0.99_core_rknd * max_Skx2_pos_Skx_sgn_wpxp )
     496             :           elsewhere
     497             :              Skx = -sqrt( 0.99_core_rknd * max_Skx2_pos_Skx_sgn_wpxp )
     498             :           endwhere
     499             :        endwhere ! Skx^2 >= max_Skx2_pos_Skx_sgn_wpxp
     500             :     elsewhere ! Skx * sgn( <w'x'> ) < 0
     501             :        where ( Skx**2 >= max_Skx2_neg_Skx_sgn_wpxp )
     502             :           where ( Skx >= zero )
     503             :              Skx = sqrt( 0.99_core_rknd * max_Skx2_neg_Skx_sgn_wpxp )
     504             :           elsewhere
     505             :              Skx = -sqrt( 0.99_core_rknd * max_Skx2_neg_Skx_sgn_wpxp )
     506             :           endwhere
     507             :        endwhere ! Skx^2 >= max_Skx2_neg_Skx_sgn_wpxp
     508             :     endwhere ! Skx * sgn( <w'x'> ) >= 0
     509             : 
     510             :     call calc_limits_F_x_responder( nz, mixt_frac, Skx, sgn_wpxp,  & ! In
     511             :                                     max_Skx2_pos_Skx_sgn_wpxp, & ! In
     512             :                                     max_Skx2_neg_Skx_sgn_wpxp, & ! In
     513           0 :                                     min_F_x, max_F_x )           ! Out
     514             : 
     515             :     ! F_x must have a value between min_F_x and max_F_x.
     516             :     F_x = calc_F_x_responder( nz, coef_spread_DG_means_x, exp_factor_x, &
     517           0 :                               min_F_x, max_F_x )
     518             : 
     519             :     call calc_responder_params( nz, xm, xp2, Skx, sgn_wpxp,       & ! In
     520             :                                 F_x, mixt_frac,               & ! In
     521             :                                 mu_x_1, mu_x_2,               & ! Out
     522             :                                 sigma_x_1_sqd, sigma_x_2_sqd, & ! Out
     523             :                                 coef_sigma_x_1_sqd,           & ! Out
     524           0 :                                 coef_sigma_x_2_sqd            ) ! Out
     525             : 
     526             : 
     527           0 :     return
     528             : 
     529             :   end subroutine calc_responder_var
     530             : 
     531             :   !=============================================================================
     532           0 :   subroutine calc_F_x_zeta_x_setter( nz, Skx,                          & ! In
     533             :                                      slope_coef_spread_DG_means_x, & ! In
     534             :                                      pdf_component_stdev_factor_x, & ! In
     535             :                                      lambda,                       & ! In
     536           0 :                                      F_x, zeta_x,                  & ! Out
     537           0 :                                      min_F_x, max_F_x              ) ! Out
     538             : 
     539             :     ! Description:
     540             :     ! Calculates the values of F_x and zeta_x for the setter variable (which is
     541             :     ! the variable that sets the mixture fraction).
     542             :     !
     543             :     ! The value of F_x is calculated between 0 (min_F_x) and 1 (max_F_x).  The
     544             :     ! equation is:
     545             :     !
     546             :     ! F_x = max_F_x + ( min_F_x - max_F_x )
     547             :     !                 * exp{ -|Skx|^lambda / slope_coef_spread_DG_means_x };
     548             :     !
     549             :     ! which reduces to:
     550             :     !
     551             :     ! F_x = 1 - exp{ -|Skx|^lambda / slope_coef_spread_DG_means_x };
     552             :     !
     553             :     ! where lambda > 0 and slope_coef_spread_DG_means_x > 0.  As |Skx| goes
     554             :     ! toward 0, the value of F_x goes toward 0, and as |Skx| becomes large, the
     555             :     ! value of F_x goes toward 1.  When slope_coef_spread_DG_means_x is small,
     556             :     ! the value of F_x tends toward 1, and when slope_coef_spread_DG_means_x is
     557             :     ! large, the value of F_x tends toward 0.  When lambda is small, the value
     558             :     ! of F_x is less dependent on Skx, and when lambda is large, the value of
     559             :     ! F_x is more dependent on Skx.
     560             :     !
     561             :     ! Mathematically, this equation will always produce a value of F_x that
     562             :     ! falls between min_F_x and max_F_x.  However, in order to prevent a value
     563             :     ! of F_x from being calculated outside the bounds of min_F_x and max_F_x
     564             :     ! owing to numerical underflow or loss of precision, this equation can be
     565             :     ! rewritten as:
     566             :     !
     567             :     ! F_x
     568             :     ! = min_F_x * exp{ -|Skx|^lambda / slope_coef_spread_DG_means_x }
     569             :     !   + max_F_x * ( 1 - exp{ -|Skx|^lambda / slope_coef_spread_DG_means_x } ).
     570             :     !
     571             :     ! The value of zeta_x used to adjust the PDF component standard devations:
     572             :     !
     573             :     ! 1 + zeta_x = ( mixt_frac * sigma_x_1^2 )
     574             :     !              / ( ( 1 - mixt_frac ) * sigma_x_2^2 );
     575             :     !
     576             :     ! where zeta_x > -1.  The sign of zeta_x is used to easily determine if
     577             :     ! mixt_frac * sigma_x_1^2 is greater than ( 1 - mixt_frac ) * sigma_x_2^2
     578             :     ! (when zeta_x is positive), mixt_frac * sigma_x_1^2 is less than
     579             :     ! ( 1 - mixt_frac ) * sigma_x_2^2 (when zeta_x is negative), or
     580             :     ! mixt_frac * sigma_x_1^2 is equal to ( 1 - mixt_frac ) * sigma_x_2^2 (when
     581             :     ! zeta_x is 0).
     582             :     !
     583             :     ! In order to allow for a tunable parameter that is the pure ratio of
     584             :     ! mixt_frac * sigma_x_1^2 to ( 1 - mixt_frac ) * sigma_x_2^2, zeta_x is
     585             :     ! related to the parameter pdf_component_stdev_factor_x, where:
     586             :     !
     587             :     ! 1 + zeta_x = pdf_component_stdev_factor_x.
     588             : 
     589             :     ! References:
     590             :     !-----------------------------------------------------------------------
     591             : 
     592             :     use constants_clubb, only: &
     593             :         one,  & ! Variable(s)
     594             :         zero
     595             : 
     596             :     use clubb_precision, only: &
     597             :         core_rknd    ! Variable(s)
     598             : 
     599             :     implicit none
     600             : 
     601             :     integer, intent(in) :: &
     602             :       nz
     603             : 
     604             :     ! Input Variables
     605             :     real( kind = core_rknd ), dimension(nz), intent(in) :: &
     606             :       Skx    ! Skewness of x (overall)              [-]
     607             : 
     608             :     real( kind = core_rknd ), intent(in) :: &
     609             :       slope_coef_spread_DG_means_x, & ! Slope coef: spread PDF comp. means x [-]
     610             :       pdf_component_stdev_factor_x, & ! Param.: PDF comp. standard devs.; x  [-]
     611             :       lambda                          ! Param. for Skx dependence            [-]
     612             : 
     613             :     ! Output Variables
     614             :     real( kind = core_rknd ), dimension(nz), intent(out) :: &
     615             :       F_x,     & ! Parameter for the spread of the PDF component means of x  [-]
     616             :       zeta_x,  & ! Parameter for the PDF component variances of x            [-]
     617             :       min_F_x, & ! Minimum allowable value of parameter F_x                  [-]
     618             :       max_F_x    ! Maximum allowable value of parameter F_x                  [-]
     619             : 
     620             :     ! Local Variable
     621             :     real( kind = core_rknd ), dimension(nz) :: &
     622           0 :       exp_Skx_interp_factor    ! Function to interp. between min. and max.   [-]
     623             : 
     624             : 
     625             :     ! Set min_F_x to 0 and max_F_x to 1 for the setter variable.
     626           0 :     where ( abs( Skx ) > zero )
     627             :        min_F_x = 1.0e-3_core_rknd
     628             :     elsewhere
     629             :        min_F_x = zero
     630             :     endwhere
     631           0 :     max_F_x = one
     632             : 
     633             :     ! F_x must have a value between min_F_x and max_F_x.
     634             :     exp_Skx_interp_factor &
     635           0 :     = exp( -abs(Skx)**lambda / slope_coef_spread_DG_means_x )
     636             : 
     637             :     F_x = min_F_x * exp_Skx_interp_factor &
     638           0 :           + max_F_x * ( one - exp_Skx_interp_factor )
     639             : 
     640             :     ! The value of zeta_x must be greater than -1.
     641           0 :     zeta_x = pdf_component_stdev_factor_x - one
     642             : 
     643             : 
     644           0 :     return
     645             : 
     646             :   end subroutine calc_F_x_zeta_x_setter
     647             : 
     648             :   !=============================================================================
     649           0 :   function calc_F_x_responder( nz, coef_spread_DG_means_x, exp_factor_x, &
     650           0 :                                min_F_x, max_F_x ) &
     651           0 :   result( F_x )
     652             : 
     653             :     ! Description:
     654             :     ! Calculates the value of F_x as a tunable function between min_F_x and
     655             :     ! max_F_x.
     656             :     !
     657             :     ! The value of F_x is calculated between min_F_x and max_F_x.  The equation
     658             :     ! is:
     659             :     !
     660             :     ! F_x = min_F_x + ( max_F_x - min_F_x )
     661             :     !                 * coef_spread_DG_means_x * exp_factor_x;
     662             :     !
     663             :     ! where 0 <= coef_spread_DG_means_x <= 1.  As coef_spread_DG_means_x
     664             :     ! goes toward 0, the value of F_x goes toward min_F_x, and as
     665             :     ! coef_spread_DG_means_x goes toward 1, the value of F_x goes toward
     666             :     ! max_F_x.  The exp_factor_x is a factor of the form 1 - exp{ }.  The range
     667             :     ! of values of exp_factor_x is 0 <= exp_factor_x <= 1.  Here, exp_factor_x
     668             :     ! is used to reduce the value of F_x under special conditions.
     669             :     !
     670             :     ! The equation for exp_factor_x depends on which responder variable is being
     671             :     ! solved for.  For rt, the F_rt equation is:
     672             :     !
     673             :     ! F_rt = min_F_rt + ( max_F_rt - min_F_rt )
     674             :     !                   * coef_spread_DG_means_rt * exp_factor_rt;
     675             :     !
     676             :     ! where exp_factor_rt is given by:
     677             :     !
     678             :     ! exp_factor_rt = 1 - exp{ -0.2 * ( adj_corr_rt_thl + 1 )^5 };
     679             :     !
     680             :     ! and where the adjusted (overall) correlation of rt and theta-l, denoted
     681             :     ! adj_corr_rt_thl, is given by:
     682             :     !
     683             :     ! adj_corr_rt_thl = <rt'thl'> / ( sqrt( <rt'^2> ) * sqrt( <thl'^2> ) )
     684             :     !                   * sgn( <w'rt'> ) * sgn( <w'thl'> ).
     685             :     !
     686             :     ! The range of values of the adjusted (overall) correlation of rt and
     687             :     ! theta-l is -1 <= adj_corr_rt_thl <= 1.  The adjusted (overall) correlation
     688             :     ! of rt and theta-l has the same absolute magnitude as the (overall)
     689             :     ! correlation of rt and theta-l, but the sign of the adjusted correlation is
     690             :     ! dependent on the agreement between the signs of <w'rt'> and <w'thl'>.
     691             :     ! When all three covariances are in sign agreement, which is when:
     692             :     !
     693             :     ! sgn( <w'rt'> ) * sgn( <w'thl'> ) * sgn( <rt'thl'> ) = 1,
     694             :     !
     695             :     ! the value of adj_corr_rt_thl is positive.  However, when the three
     696             :     ! covariances aren't consistent in sign, which is when:
     697             :     !
     698             :     ! sgn( <w'rt'> ) * sgn( <w'thl'> ) * sgn( <rt'thl'> ) = -1,
     699             :     !
     700             :     ! the value of adj_corr_rt_thl is negative.
     701             :     !
     702             :     ! For theta-l, the F_thl equation is:
     703             :     !
     704             :     ! F_thl = min_F_thl + ( max_F_thl - min_F_thl ) * coef_spread_DG_means_thl;
     705             :     !
     706             :     ! where exp_factor_thl is set to 1 (1 - exp{-inf} = 1) because reducing
     707             :     ! F_thl through the use of exp_factor_thl is not desired for theta-l.
     708             :     !
     709             :     ! The direction of the two PDF component means of x (mu_x_1 and mu_x_2) with
     710             :     ! respect to each other and the overall mean of x (<x>) is given by
     711             :     ! sgn( <w'x'> ).  When sgn( <w'x'> ) = 1, ( mu_x_1 - <x> ) >= 0 and
     712             :     ! ( mu_x_2 - <x> ) <= 0.  When sgn( <w'x'> ) = -1, ( mu_x_1 - <x> ) <= 0 and
     713             :     ! ( mu_x_2 - <x> ) >= 0.  This helps to promote realizability of the PDF
     714             :     ! component correlations (corr_w_x_1 and corr_w_x_2).
     715             :     !
     716             :     ! The realizability of the PDF component correlations corr_w_rt_1,
     717             :     ! corr_w_rt_2, corr_w_thl_1, and corr_w_thl_2 is promoted by using
     718             :     ! sgn( <w'rt'> ) and sgn( <w'thl'> ), respectively to choose the direction
     719             :     ! of the PDF component means of rt and theta-l.  However, when the three
     720             :     ! covariances (<w'rt'>, <w'thl'>, and <rt'thl'>) aren't consistent in sign
     721             :     ! (as shown above), it becomes difficult to keep the PDF component
     722             :     ! correlations corr_rt_thl_1 and corr_rt_thl_2 in the realizable range.
     723             :     ! However, in this situation, the above equation for exp_factor_rt brings
     724             :     ! F_rt closer to min_F_rt, which promotes realizability for corr_rt_thl_1
     725             :     ! and corr_rt_thl_2.
     726             :     !
     727             :     ! Mathematically, this equation will always produce a value of F_x that
     728             :     ! falls between min_F_x and max_F_x.  However, in order to prevent a value
     729             :     ! of F_x from being calculated outside the bounds of min_F_x and max_F_x
     730             :     ! owing to numerical underflow or loss of precision, this equation can be
     731             :     ! rewritten as:
     732             :     !
     733             :     ! F_x = min_F_x * ( 1 - coef_spread_DG_means_x * exp_factor_x )
     734             :     !       + max_F_x * coef_spread_DG_means_x * exp_factor_x.
     735             : 
     736             :     ! References:
     737             :     !-----------------------------------------------------------------------
     738             : 
     739             :     use constants_clubb, only: &
     740             :         one    ! Variable(s)
     741             : 
     742             :     use clubb_precision, only: &
     743             :         core_rknd    ! Variable(s)
     744             : 
     745             :     implicit none
     746             : 
     747             :     integer, intent(in) :: &
     748             :       nz
     749             : 
     750             :     ! Input Variables
     751             :     real( kind = core_rknd ), intent(in) :: &
     752             :       coef_spread_DG_means_x    ! Coef.: spread betw. PDF comp. means of x   [-]
     753             : 
     754             :     real( kind = core_rknd ), dimension(nz), intent(in) :: &
     755             :       exp_factor_x, & ! Factor of the form 1 - exp{}; reduces F_x  [-]
     756             :       min_F_x,      & ! Minimum allowable value of parameter F_x   [-]
     757             :       max_F_x         ! Maximum allowable value of parameter F_x   [-]
     758             : 
     759             :     ! Return Variable
     760             :     real( kind = core_rknd ), dimension(nz) :: &
     761             :       F_x    ! Parameter for the spread between the PDF component means of x [-]
     762             : 
     763             : 
     764             :     ! F_x must have a value between min_F_x and max_F_x.
     765             :     F_x = min_F_x * ( one - coef_spread_DG_means_x * exp_factor_x ) &
     766           0 :           + max_F_x * coef_spread_DG_means_x * exp_factor_x
     767             : 
     768             : 
     769           0 :     return
     770             : 
     771           0 :   end function calc_F_x_responder
     772             : 
     773             :   !=============================================================================
     774             : 
     775             : end module new_pdf_main

Generated by: LCOV version 1.14