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 22:39:59 |  | 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