LCOV - code coverage report
Current view: | top level - physics/clubb/src/CLUBB_core - pdf_closure_module.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 347 | 724 | 47.9 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 9 | 11 | 81.8 %

              Line data    Source code
           1             : !---------------------------------------------------------------------------
           2             : ! $Id$
           3             : !===============================================================================
           4             : module pdf_closure_module
           5             :
           6             :   ! Options for the two component normal (double Gaussian) PDF type to use for
           7             :   ! the w, rt, and theta-l (or w, chi, and eta) portion of CLUBB's multivariate,
           8             :   ! two-component PDF.
           9             :   use model_flags, only: &
          10             :       iiPDF_ADG1,       & ! ADG1 PDF
          11             :       iiPDF_ADG2,       & ! ADG2 PDF
          12             :       iiPDF_3D_Luhar,   & ! 3D Luhar PDF
          13             :       iiPDF_new,        & ! new PDF
          14             :       iiPDF_TSDADG,     & ! new TSDADG PDF
          15             :       iiPDF_LY93,       & ! Lewellen and Yoh (1993)
          16             :       iiPDF_new_hybrid    ! new hybrid PDF
          17             :
          18             :   implicit none
          19             :
          20             :   public :: pdf_closure, &
          21             :             calc_wp4_pdf, &
          22             :             calc_wp2xp_pdf, &
          23             :             calc_wpxp2_pdf, &
          24             :             calc_wpxpyp_pdf, &
          25             :             calc_w_up_in_cloud
          26             :
          27             :   private ! Set Default Scope
          28             :
          29             :   contains
          30             : !------------------------------------------------------------------------
          31             :
          32             :   !#######################################################################
          33             :   !#######################################################################
          34             :   ! If you change the argument list of pdf_closure you also have to
          35             :   ! change the calls to this function in the host models CAM, WRF, SAM
          36             :   ! and GFDL.
          37             :   !#######################################################################
          38             :   !#######################################################################
          39    17870112 :   subroutine pdf_closure( nz, ngrdcol, sclr_dim, sclr_tol,            &
          40    17870112 :                           hydromet_dim, p_in_Pa, exner, thv_ds,       &
          41    17870112 :                           wm, wp2, wp3,                               &
          42    17870112 :                           Skw, Skthl_in, Skrt_in, Sku_in, Skv_in,     &
          43    17870112 :                           rtm, rtp2, wprtp,                           &
          44    17870112 :                           thlm, thlp2, wpthlp,                        &
          45    17870112 :                           um, up2, upwp,                              &
          46    17870112 :                           vm, vp2, vpwp,                              &
          47    17870112 :                           rtpthlp,                                    &
          48    17870112 :                           sclrm, wpsclrp, sclrp2,                     &
          49    17870112 :                           sclrprtp, sclrpthlp, Sksclr_in,             &
          50    17870112 :                           gamma_Skw_fnc,                              &
          51             : #ifdef GFDL
          52             :                           RH_crit, do_liquid_only_in_clubb,           & ! h1g, 2010-06-15
          53             : #endif
          54    17870112 :                           wphydrometp, wp2hmp,                        &
          55    17870112 :                           rtphmp, thlphmp,                            &
          56    17870112 :                           clubb_params,                               &
          57             :                           saturation_formula,                         &
          58             :                           stats_metadata,                             &
          59             :                           iiPDF_type,                                 &
          60    17870112 :                           l_mix_rat_hm,                               &
          61    17870112 :                           sigma_sqd_w,                                &
          62             :                           pdf_params, pdf_implicit_coefs_terms,       &
          63    17870112 :                           wpup2, wpvp2,                               &
          64    17870112 :                           wp2up2, wp2vp2, wp4,                        &
          65    17870112 :                           wprtp2, wp2rtp,                             &
          66    17870112 :                           wpthlp2, wp2thlp, wprtpthlp,                &
          67    17870112 :                           cloud_frac, ice_supersat_frac,              &
          68    17870112 :                           rcm, wpthvp, wp2thvp, rtpthvp,              &
          69    17870112 :                           thlpthvp, wprcp, wp2rcp, rtprcp,            &
          70    17870112 :                           thlprcp, rcp2,                              &
          71    17870112 :                           uprcp, vprcp,                               &
          72    17870112 :                           w_up_in_cloud, w_down_in_cloud,             &
          73    17870112 :                           cloudy_updraft_frac, cloudy_downdraft_frac, &
          74    17870112 :                           F_w, F_rt, F_thl,                           &
          75    17870112 :                           min_F_w, max_F_w,                           &
          76    17870112 :                           min_F_rt, max_F_rt,                         &
          77    17870112 :                           min_F_thl, max_F_thl,                       &
          78    17870112 :                           wpsclrprtp, wpsclrp2, sclrpthvp,            &
          79    17870112 :                           wpsclrpthlp, sclrprcp, wp2sclrp,            &
          80    17870112 :                           rc_coef                                     )
          81             :
          82             :
          83             :     ! Description:
          84             :     ! Subroutine that computes pdf parameters analytically.
          85             :     !
          86             :     ! Based of the original formulation, but with some tweaks
          87             :     ! to remove some of the less realistic assumptions and
          88             :     ! improve transport terms.
          89             :
          90             :     !   Corrected version that should remove inconsistency
          91             :
          92             :     ! References:
          93             :     !   The shape of CLUBB's PDF is given by the expression in
          94             :     !   https://arxiv.org/pdf/1711.03675v1.pdf#nameddest=url:clubb_pdf
          95             :
          96             :     !   Eqn. 29, 30, 31, 32 & 33  on p. 3547 of
          97             :     !   ``A PDF-Based Model for Boundary Layer Clouds. Part I:
          98             :     !   Method and Model Description'' Golaz, et al. (2002)
          99             :     !   JAS, Vol. 59, pp. 3540--3551.
         100             :     !----------------------------------------------------------------------
         101             :
         102             :     use grid_class, only: &
         103             :         grid ! Type
         104             :
         105             :     use constants_clubb, only: &  ! Constants
         106             :         three,          & ! 3
         107             :         one,            & ! 1
         108             :         one_half,       & ! 1/2
         109             :         zero,           & ! 0
         110             :         Cp,             & ! Dry air specific heat at constant p [J/kg/K]
         111             :         Lv,             & ! Latent heat of vaporization         [J/kg]
         112             :         ep1,            & ! (1.0-ep)/ep; ep1 = 0.61             [-]
         113             :         ep2,            & ! 1.0/ep;      ep2 = 1.61             [-]
         114             :         rt_tol,         & ! Tolerance for r_t                   [kg/kg]
         115             :         thl_tol,        & ! Tolerance for th_l                  [K]
         116             :         fstderr,        &
         117             :         zero_threshold, &
         118             :         eps, &
         119             :         w_tol
         120             :
         121             :     use parameters_model, only: &
         122             :         mixt_frac_max_mag  ! Variable(s)
         123             :
         124             :     use parameter_indices, only: &
         125             :         nparams,                       & ! Variable(s)
         126             :         ibeta,                         &
         127             :         iSkw_denom_coef,               &
         128             :         ipdf_component_stdev_factor_w, &
         129             :         islope_coef_spread_DG_means_w
         130             :
         131             :     use pdf_parameter_module, only:  &
         132             :         pdf_parameter,        & ! Variable Type
         133             :         implicit_coefs_terms
         134             :
         135             :     use new_pdf_main, only: &
         136             :         new_pdf_driver    ! Procedure(s)
         137             :
         138             :     use new_hybrid_pdf_main, only: &
         139             :         new_hybrid_pdf_driver    ! Procedure(s)
         140             :
         141             :     use adg1_adg2_3d_luhar_pdf, only: &
         142             :         ADG1_pdf_driver,     & ! Procedure(s)
         143             :         ADG2_pdf_driver,     &
         144             :         Luhar_3D_pdf_driver
         145             :
         146             :     use new_tsdadg_pdf, only: &
         147             :         tsdadg_pdf_driver    ! Procedure(s)
         148             :
         149             :     use LY93_pdf, only: &
         150             :         LY93_driver    ! Procedure(s)
         151             :
         152             :     use pdf_utilities, only: &
         153             :         calc_comp_corrs_binormal, & ! Procedure(s)
         154             :         calc_corr_chi_x,          &
         155             :         calc_corr_eta_x
         156             :
         157             :     use model_flags, only: &
         158             :         l_explicit_turbulent_adv_xpyp ! Variable(s)
         159             :
         160             :     use numerical_check, only:  &
         161             :         pdf_closure_check ! Procedure(s)
         162             :
         163             :     use saturation, only:  &
         164             :         sat_mixrat_liq, & ! Procedure(s)
         165             :         sat_mixrat_ice
         166             :
         167             :     use clubb_precision, only: &
         168             :         core_rknd ! Variable(s)
         169             :
         170             :     use error_code, only: &
         171             :         clubb_at_least_debug_level,  & ! Procedure
         172             :         err_code,                    & ! Error Indicator
         173             :         clubb_fatal_error              ! Constant
         174             :
         175             :     use stats_variables, only: &
         176             :         stats_metadata_type
         177             :
         178             :     implicit none
         179             :
         180             :     !----------------------------- Input Variables -----------------------------
         181             :     integer, intent(in) :: &
         182             :       nz,           & ! Number of vertical levels
         183             :       ngrdcol,      & ! Number of grid columns
         184             :       hydromet_dim, & ! Number of hydrometeor species
         185             :       sclr_dim        ! Number of passive scalars
         186             :
         187             :     real( kind = core_rknd ), intent(in), dimension(sclr_dim) :: &
         188             :       sclr_tol          ! Threshold(s) on the passive scalars  [units vary]
         189             :
         190             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) ::  &
         191             :       p_in_Pa,     & ! Pressure                                   [Pa]
         192             :       exner,       & ! Exner function                             [-]
         193             :       thv_ds,      & ! Dry, base-state theta_v (ref. th_l here)   [K]
         194             :       wm,          & ! mean w-wind component (vertical velocity)  [m/s]
         195             :       wp2,         & ! w'^2                                       [m^2/s^2]
         196             :       wp3,         & ! w'^3                                       [m^3/s^3]
         197             :       Skw,         & ! Skewness of w                              [-]
         198             :       Skthl_in,    & ! Skewness of thl                            [-]
         199             :       Skrt_in,     & ! Skewness of rt                             [-]
         200             :       Sku_in,      & ! Skewness of u                              [-]
         201             :       Skv_in,      & ! Skewness of v                              [-]
         202             :       rtm,         & ! Mean total water mixing ratio              [kg/kg]
         203             :       rtp2,        & ! r_t'^2                                     [(kg/kg)^2]
         204             :       wprtp,       & ! w'r_t'                                     [(kg/kg)(m/s)]
         205             :       thlm,        & ! Mean liquid water potential temperature    [K]
         206             :       thlp2,       & ! th_l'^2                                    [K^2]
         207             :       wpthlp,      & ! w'th_l'                                    [K(m/s)]
         208             :       rtpthlp        ! r_t'th_l'                                  [K(kg/kg)]
         209             :
         210             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) ::  &
         211             :       um,          & ! Grid-mean eastward wind     [m/s]
         212             :       up2,         & ! u'^2                        [(m/s)^2]
         213             :       upwp,        & ! u'w'                        [(m/s)^2]
         214             :       vm,          & ! Grid-mean northward wind    [m/s]
         215             :       vp2,         & ! v'^2                        [(m/s)^2]
         216             :       vpwp           ! v'w'                        [(m/s)^2]
         217             :
         218             :     real( kind = core_rknd ), dimension(ngrdcol,nz, sclr_dim), intent(in) ::  &
         219             :       sclrm,       & ! Mean passive scalar        [units vary]
         220             :       wpsclrp,     & ! w' sclr'                   [units vary]
         221             :       sclrp2,      & ! sclr'^2                    [units vary]
         222             :       sclrprtp,    & ! sclr' r_t'                 [units vary]
         223             :       sclrpthlp,   & ! sclr' th_l'                [units vary]
         224             :       Sksclr_in      ! Skewness of sclr           [-]
         225             :
         226             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) ::  &
         227             :       gamma_Skw_fnc    ! Gamma as a function of skewness            [-]
         228             :
         229             : #ifdef  GFDL
         230             :     ! critial relative humidity for nucleation
         231             :     real( kind = core_rknd ), dimension(ngrdcol, nz, min(1,sclr_dim), 2 ), intent(in) ::  & ! h1g, 2010-06-15
         232             :        RH_crit     ! critical relative humidity for droplet and ice nucleation
         233             : ! ---> h1g, 2012-06-14
         234             :     logical, intent(in)                 ::  do_liquid_only_in_clubb
         235             : ! <--- h1g, 2012-06-14
         236             : #endif
         237             :
         238             :     real( kind = core_rknd ), dimension(ngrdcol,nz,hydromet_dim), intent(in) :: &
         239             :       wphydrometp, & ! Covariance of w and a hydrometeor    [(m/s) <hm units>]
         240             :       wp2hmp,      & ! Third-order moment:  < w'^2 hm' >    [(m/s)^2 <hm units>]
         241             :       rtphmp,      & ! Covariance of rt and a hydrometeor   [(kg/kg) <hm units>]
         242             :       thlphmp        ! Covariance of thl and a hydrometeor  [K <hm units>]
         243             :
         244             :     real( kind = core_rknd ), dimension(ngrdcol,nparams), intent(in) :: &
         245             :       clubb_params    ! Array of CLUBB's tunable parameters    [units vary]
         246             :
         247             :     integer, intent(in) :: &
         248             :       iiPDF_type    ! Selected option for the two-component normal (double
         249             :                     ! Gaussian) PDF type to use for the w, rt, and theta-l (or
         250             :                     ! w, chi, and eta) portion of CLUBB's multivariate,
         251             :                     ! two-component PDF.
         252             :
         253             :     logical, dimension(hydromet_dim), intent(in) :: &
         254             :       l_mix_rat_hm   ! if true, then the quantity is a hydrometeor mixing ratio
         255             :
         256             :     integer, intent(in) :: &
         257             :       saturation_formula ! Integer that stores the saturation formula to be used
         258             :
         259             :     type (stats_metadata_type), intent(in) :: &
         260             :       stats_metadata
         261             :
         262             :     !----------------------------- InOut Variables -----------------------------
         263             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(inout) :: &
         264             :       ! If iiPDF_type == iiPDF_ADG2, this gets overwritten. Therefore,
         265             :       ! intent(inout). Otherwise it should be intent(in)
         266             :       sigma_sqd_w   ! Width of individual w plumes               [-]
         267             :
         268             :     type(pdf_parameter), intent(inout) :: &
         269             :       pdf_params     ! pdf paramters         [units vary]
         270             :
         271             :     type(implicit_coefs_terms), intent(inout) :: &
         272             :       pdf_implicit_coefs_terms    ! Implicit coefs / explicit terms [units vary]
         273             :
         274             :     !----------------------------- Output Variables -----------------------------
         275             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) ::  &
         276             :       wpup2,                 & ! w'u'^2                     [m^3/s^3]
         277             :       wpvp2,                 & ! w'v'^2                     [m^3/s^3]
         278             :       wp2up2,                & ! w'^2u'^2                   [m^2/s^4]
         279             :       wp2vp2,                & ! w'^2v'^2                   [m^2/s^4]
         280             :       wp4,                   & ! w'^4                       [m^4/s^4]
         281             :       wprtp2,                & ! w' r_t'                    [(m kg)/(s kg)]
         282             :       wp2rtp,                & ! w'^2 r_t'                  [(m^2 kg)/(s^2 kg)]
         283             :       wpthlp2,               & ! w' th_l'^2                 [(m K^2)/s]
         284             :       wp2thlp,               & ! w'^2 th_l'                 [(m^2 K)/s^2]
         285             :       cloud_frac,            & ! Cloud fraction             [-]
         286             :       ice_supersat_frac,     & ! Ice cloud fracion          [-]
         287             :       rcm,                   & ! Mean liquid water          [kg/kg]
         288             :       wpthvp,                & ! Buoyancy flux              [(K m)/s]
         289             :       wp2thvp,               & ! w'^2 th_v'                 [(m^2 K)/s^2]
         290             :       rtpthvp,               & ! r_t' th_v'                 [(kg K)/kg]
         291             :       thlpthvp,              & ! th_l' th_v'                [K^2]
         292             :       wprcp,                 & ! w' r_c'                    [(m kg)/(s kg)]
         293             :       wp2rcp,                & ! w'^2 r_c'                  [(m^2 kg)/(s^2 kg)]
         294             :       rtprcp,                & ! r_t' r_c'                  [(kg^2)/(kg^2)]
         295             :       thlprcp,               & ! th_l' r_c'                 [(K kg)/kg]
         296             :       rcp2,                  & ! r_c'^2                     [(kg^2)/(kg^2)]
         297             :       wprtpthlp,             & ! w' r_t' th_l'              [(m kg K)/(s kg)]
         298             :       w_up_in_cloud,         & ! cloudy updraft vel         [m/s]
         299             :       w_down_in_cloud,       & ! cloudy downdraft vel       [m/s]
         300             :       cloudy_updraft_frac,   & ! cloudy updraft fraction    [-]
         301             :       cloudy_downdraft_frac    ! cloudy downdraft fraction  [-]
         302             :
         303             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) ::  &
         304             :       uprcp,              & ! u' r_c'               [(m kg)/(s kg)]
         305             :       vprcp                 ! v' r_c'               [(m kg)/(s kg)]
         306             :
         307             :     ! Parameters output only for recording statistics (new PDF).
         308             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
         309             :       F_w,   & ! Parameter for the spread of the PDF component means of w    [-]
         310             :       F_rt,  & ! Parameter for the spread of the PDF component means of rt   [-]
         311             :       F_thl    ! Parameter for the spread of the PDF component means of thl  [-]
         312             :
         313             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
         314             :       min_F_w,   & ! Minimum allowable value of parameter F_w      [-]
         315             :       max_F_w,   & ! Maximum allowable value of parameter F_w      [-]
         316             :       min_F_rt,  & ! Minimum allowable value of parameter F_rt     [-]
         317             :       max_F_rt,  & ! Maximum allowable value of parameter F_rt     [-]
         318             :       min_F_thl, & ! Minimum allowable value of parameter F_thl    [-]
         319             :       max_F_thl    ! Maximum allowable value of parameter F_thl    [-]
         320             :
         321             :     ! Output (passive scalar variables)
         322             :     real( kind = core_rknd ), intent(out), dimension(ngrdcol,nz,sclr_dim) ::  &
         323             :       sclrpthvp, &
         324             :       sclrprcp, &
         325             :       wpsclrp2, &
         326             :       wpsclrprtp, &
         327             :       wpsclrpthlp, &
         328             :       wp2sclrp
         329             :
         330             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
         331             :       rc_coef    ! Coefficient on X'r_c' in X'th_v' equation    [K/(kg/kg)]
         332             :
         333             :     !----------------------------- Local Variables -----------------------------
         334             :
         335             :     ! Variables that are stored in derived data type pdf_params.
         336             :     real( kind = core_rknd ), dimension(ngrdcol,nz) ::  &
         337    35740224 :       u_1,           & ! Mean of eastward wind (1st PDF component)         [m/s]
         338    35740224 :       u_2,           & ! Mean of eastward wind (2nd PDF component)         [m/s]
         339    35740224 :       varnce_u_1,    & ! Variance of u (1st PDF component)             [m^2/s^2]
         340    35740224 :       varnce_u_2,    & ! Variance of u (2nd PDF component)             [m^2/s^2]
         341    35740224 :       v_1,           & ! Mean of northward wind (1st PDF component)        [m/s]
         342    35740224 :       v_2,           & ! Mean of northward wind (2nd PDF component)        [m/s]
         343    35740224 :       varnce_v_1,    & ! Variance of v (1st PDF component)             [m^2/s^2]
         344    35740224 :       varnce_v_2,    & ! Variance of v (2nd PDF component)             [m^2/s^2]
         345    35740224 :       alpha_u,       & ! Factor relating to normalized variance for u        [-]
         346    35740224 :       alpha_v          ! Factor relating to normalized variance for v        [-]
         347             :
         348             :     real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
         349    35740224 :       corr_u_w_1,      & ! Correlation of u and w   (1st PDF component)      [-]
         350    35740224 :       corr_u_w_2,      & ! Correlation of u and w   (2nd PDF component)      [-]
         351    35740224 :       corr_v_w_1,      & ! Correlation of v and w   (1st PDF component)      [-]
         352    35740224 :       corr_v_w_2         ! Correlation of v and w   (2nd PDF component)      [-]
         353             :
         354             :     ! Note:  alpha coefficients = 0.5 * ( 1 - correlations^2 ).
         355             :     !        These are used to calculate the scalar widths
         356             :     !        varnce_thl_1, varnce_thl_2, varnce_rt_1, and varnce_rt_2 as in
         357             :     !        Eq. (34) of Larson and Golaz (2005)
         358             :
         359             :     ! Passive scalar local variables
         360             :
         361             :     real( kind = core_rknd ), dimension(ngrdcol,nz,sclr_dim) ::  &
         362    35740224 :       sclr1, sclr2,  &
         363    35740224 :       varnce_sclr1, varnce_sclr2, &
         364    35740224 :       alpha_sclr,  &
         365    35740224 :       corr_sclr_thl_1, corr_sclr_thl_2, &
         366    35740224 :       corr_sclr_rt_1, corr_sclr_rt_2, &
         367    35740224 :       corr_w_sclr_1, corr_w_sclr_2
         368             :
         369             :     logical :: &
         370             :       l_scalar_calc, &  ! True if sclr_dim > 0
         371             :       l_calc_ice_supersat_frac ! True if we should calculate ice_supersat_frac
         372             :
         373             :     ! Quantities needed to predict higher order moments
         374             :     real( kind = core_rknd ), dimension(ngrdcol,nz) ::  &
         375    35740224 :       tl1, tl2
         376             :
         377             :     real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
         378    35740224 :       sqrt_wp2, & ! Square root of wp2          [m/s]
         379    35740224 :       Skthl,    & ! Skewness of thl             [-]
         380    35740224 :       Skrt,     & ! Skewness of rt              [-]
         381    35740224 :       Sku,      & ! Skewness of u               [-]
         382    35740224 :       Skv         ! Skewness of v               [-]
         383             :
         384             :     real( kind = core_rknd ), dimension(ngrdcol,nz,sclr_dim) :: &
         385    35740224 :       Sksclr      ! Skewness of rt              [-]
         386             :
         387             :     ! Thermodynamic quantity
         388             :
         389             :     real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
         390    35740224 :       wprcp_contrib_comp_1,   & ! <w'rc'> contrib. (1st PDF comp.)  [m/s(kg/kg)]
         391    35740224 :       wprcp_contrib_comp_2,   & ! <w'rc'> contrib. (2nd PDF comp.)  [m/s(kg/kg)]
         392    35740224 :       wp2rcp_contrib_comp_1,  & ! <w'^2rc'> contrib. (1st comp) [m^2/s^2(kg/kg)]
         393    35740224 :       wp2rcp_contrib_comp_2,  & ! <w'^2rc'> contrib. (2nd comp) [m^2/s^2(kg/kg)]
         394    35740224 :       rtprcp_contrib_comp_1,  & ! <rt'rc'> contrib. (1st PDF comp.)  [kg^2/kg^2]
         395    35740224 :       rtprcp_contrib_comp_2,  & ! <rt'rc'> contrib. (2nd PDF comp.)  [kg^2/kg^2]
         396    35740224 :       thlprcp_contrib_comp_1, & ! <thl'rc'> contrib. (1st PDF comp.)  [K(kg/kg)]
         397    35740224 :       thlprcp_contrib_comp_2, & ! <thl'rc'> contrib. (2nd PDF comp.)  [K(kg/kg)]
         398    35740224 :       uprcp_contrib_comp_1,   & ! <u'rc'> contrib. (1st PDF comp.)  [m/s(kg/kg)]
         399    35740224 :       uprcp_contrib_comp_2,   & ! <u'rc'> contrib. (2nd PDF comp.)  [m/s(kg/kg)]
         400    35740224 :       vprcp_contrib_comp_1,   & ! <v'rc'> contrib. (1st PDF comp.)  [m/s(kg/kg)]
         401    35740224 :       vprcp_contrib_comp_2      ! <v'rc'> contrib. (2nd PDF comp.)  [m/s(kg/kg)]
         402             :
         403             :     ! variables for computing ice cloud fraction
         404             :     real( kind = core_rknd), dimension(ngrdcol,nz) :: &
         405    35740224 :       rc_1_ice, rc_2_ice
         406             :
         407             :     ! To test pdf parameters
         408             :     real( kind = core_rknd ), dimension(ngrdcol,nz) ::  &
         409    35740224 :       wm_clubb_pdf,    &
         410    35740224 :       rtm_clubb_pdf,   &
         411    35740224 :       thlm_clubb_pdf,  &
         412    35740224 :       wp2_clubb_pdf,   &
         413    35740224 :       rtp2_clubb_pdf,  &
         414    35740224 :       thlp2_clubb_pdf, &
         415    35740224 :       wp3_clubb_pdf,   &
         416    35740224 :       rtp3_clubb_pdf,  &
         417    35740224 :       thlp3_clubb_pdf, &
         418    35740224 :       Skw_clubb_pdf,   &
         419    35740224 :       Skrt_clubb_pdf,  &
         420    35740224 :       Skthl_clubb_pdf, &
         421    35740224 :       rsatl_1, &
         422    35740224 :       rsatl_2
         423             :
         424             :     real( kind = core_rknd ) :: &
         425             :       Skw_denom_coef     ! CLUBB tunable parameter beta
         426             :
         427             :     logical, parameter :: &
         428             :       l_liq_ice_loading_test = .false. ! Temp. flag liq./ice water loading test
         429             :
         430             :     integer :: k, i, j, hm_idx   ! Indices
         431             :
         432             : #ifdef GFDL
         433             :     real ( kind = core_rknd ), parameter :: t1_combined = 273.16, &
         434             :                                             t2_combined = 268.16, &
         435             :                                             t3_combined = 238.16
         436             : #endif
         437             :
         438             :     !----------------------------- Begin Code -----------------------------
         439             :
         440             :     !$acc enter data create( u_1, u_2, varnce_u_1, varnce_u_2, v_1, v_2, &
         441             :     !$acc                 varnce_v_1, varnce_v_2, alpha_u, alpha_v, &
         442             :     !$acc                 corr_u_w_1, corr_u_w_2, corr_v_w_1, corr_v_w_2, &
         443             :     !$acc                 tl1, tl2, sqrt_wp2, Skthl, &
         444             :     !$acc                 Skrt, Sku, Skv, wprcp_contrib_comp_1, wprcp_contrib_comp_2, &
         445             :     !$acc                 wp2rcp_contrib_comp_1, wp2rcp_contrib_comp_2, &
         446             :     !$acc                 rtprcp_contrib_comp_1, rtprcp_contrib_comp_2, &
         447             :     !$acc                 thlprcp_contrib_comp_1, thlprcp_contrib_comp_2, &
         448             :     !$acc                 uprcp_contrib_comp_1, uprcp_contrib_comp_2, &
         449             :     !$acc                 vprcp_contrib_comp_1, vprcp_contrib_comp_2, &
         450             :     !$acc                 rc_1_ice, rc_2_ice, rsatl_1, rsatl_2 )
         451             :
         452             :     !$acc enter data if( sclr_dim > 0 ) &
         453             :     !$acc            create( sclr1, sclr2, varnce_sclr1, varnce_sclr2, &
         454             :     !$acc                    alpha_sclr, corr_sclr_thl_1, corr_sclr_thl_2, &
         455             :     !$acc                    corr_sclr_rt_1, corr_sclr_rt_2, corr_w_sclr_1, &
         456             :     !$acc                    corr_w_sclr_2, Sksclr )
         457             :
         458             :     ! Check whether the passive scalars are present.
         459    17870112 :     if ( sclr_dim > 0 ) then
         460           0 :       l_scalar_calc = .true.
         461             :     else
         462    17870112 :       l_scalar_calc = .false.
         463             :     end if
         464             :
         465             :     ! Initialize to default values to prevent a runtime error
         466    17870112 :     if ( ( iiPDF_type /= iiPDF_ADG1 ) .and. ( iiPDF_type /= iiPDF_ADG2 ) ) then
         467             :
         468           0 :       do k = 1, nz
         469           0 :         do i = 1, ngrdcol
         470           0 :           pdf_params%alpha_thl(i,k) = one_half
         471           0 :           pdf_params%alpha_rt(i,k) = one_half
         472             :         end do
         473             :       end do
         474             :
         475             :       ! This allows for skewness to be clipped locally without passing the updated
         476             :       ! value back out.
         477           0 :       do k = 1, nz
         478           0 :         do i = 1, ngrdcol
         479           0 :           Skrt(i,k) = Skrt_in(i,k)
         480           0 :           Skthl(i,k) = Skthl_in(i,k)
         481           0 :           Sku(i,k) = Sku_in(i,k)
         482           0 :           Skv(i,k) = Skv_in(i,k)
         483             :         end do
         484             :       end do
         485             :
         486           0 :       do j = 1, sclr_dim
         487           0 :         do k = 1, nz
         488           0 :           do i = 1, ngrdcol
         489             :
         490           0 :             Sksclr(i,k,j) = Sksclr_in(i,k,j)
         491             :
         492           0 :             if ( l_scalar_calc ) then
         493           0 :                 alpha_sclr(i,k,j) = one_half
         494             :             end if
         495             :
         496             :           end do
         497             :         end do
         498             :       end do
         499             :
         500             :     end if
         501             :
         502             :     ! Initialize to 0 to prevent a runtime error
         503    17870112 :     if ( iiPDF_type /= iiPDF_new .and. iiPDF_type /= iiPDF_new_hybrid ) then
         504             :       ! Stats only variables, setting to zero
         505  1536829632 :       do k = 1, nz
         506 25380961632 :         do i = 1, ngrdcol
         507 23844132000 :           F_w(i,k) = zero
         508 23844132000 :           F_rt(i,k) = zero
         509 23844132000 :           F_thl(i,k) = zero
         510 23844132000 :           min_F_w(i,k) = zero
         511 23844132000 :           max_F_w(i,k) = zero
         512 23844132000 :           min_F_rt(i,k) = zero
         513 23844132000 :           max_F_rt(i,k) = zero
         514 23844132000 :           min_F_thl(i,k) = zero
         515 25363091520 :           max_F_thl(i,k) = zero
         516             :         end do
         517             :       end do
         518             :     end if
         519             :
         520             :     ! To avoid recomputing
         521             :     !$acc parallel loop gang vector collapse(2) default(present)
         522  1536829632 :     do k = 1, nz
         523 25380961632 :       do i = 1, ngrdcol
         524 25363091520 :         sqrt_wp2(i,k) = sqrt( wp2(i,k) )
         525             :       end do
         526             :     end do
         527             :     !$acc end parallel loop
         528             :
         529             :     ! Select the PDF closure method for the two-component PDF used by CLUBB for
         530             :     ! w, rt, theta-l, and passive scalar variables.
         531             :     ! Calculate the mixture fraction for the multivariate PDF, as well as both
         532             :     ! PDF component means and both PDF component variances for each of w, rt,
         533             :     ! theta-l, and passive scalar variables.
         534             :     if ( iiPDF_type == iiPDF_ADG1 ) then ! use ADG1
         535             :
         536             :       call ADG1_pdf_driver( nz, ngrdcol, sclr_dim, sclr_tol,                        & ! In
         537             :                             wm, rtm, thlm, um, vm,                                  & ! In
         538             :                             wp2, rtp2, thlp2, up2, vp2,                             & ! In
         539             :                             Skw, wprtp, wpthlp, upwp, vpwp, sqrt_wp2,               & ! In
         540             :                             sigma_sqd_w, clubb_params(:,ibeta), mixt_frac_max_mag,  & ! In
         541             :                             sclrm, sclrp2, wpsclrp, l_scalar_calc,                  & ! In
         542             :                             pdf_params%w_1, pdf_params%w_2,                         & ! Out
         543             :                             pdf_params%rt_1, pdf_params%rt_2,                       & ! Out
         544             :                             pdf_params%thl_1, pdf_params%thl_2,                     & ! Out
         545             :                             u_1, u_2, v_1, v_2,                                     & ! Out
         546             :                             pdf_params%varnce_w_1, pdf_params%varnce_w_2,           & ! Out
         547             :                             pdf_params%varnce_rt_1, pdf_params%varnce_rt_2,         & ! Out
         548             :                             pdf_params%varnce_thl_1, pdf_params%varnce_thl_2,       & ! Out
         549             :                             varnce_u_1, varnce_u_2,                                 & ! Out
         550             :                             varnce_v_1, varnce_v_2,                                 & ! Out
         551             :                             pdf_params%mixt_frac,                                   & ! Out
         552             :                             pdf_params%alpha_rt, pdf_params%alpha_thl,              & ! Out
         553             :                             alpha_u, alpha_v,                                       & ! Out
         554             :                             sclr1, sclr2, varnce_sclr1,                             & ! Out
         555    17870112 :                             varnce_sclr2, alpha_sclr )                                ! Out
         556             :
         557             :     elseif ( iiPDF_type == iiPDF_ADG2 ) then ! use ADG2
         558             :
         559             :       call ADG2_pdf_driver( nz, ngrdcol, sclr_dim, sclr_tol,                      & ! In
         560             :                             wm, rtm, thlm, wp2, rtp2, thlp2,                      & ! In
         561             :                             Skw, wprtp, wpthlp, sqrt_wp2, clubb_params(:,ibeta),  & ! In
         562             :                             sclrm, sclrp2, wpsclrp, l_scalar_calc,                & ! In
         563             :                             pdf_params%w_1, pdf_params%w_2,                       & ! Out
         564             :                             pdf_params%rt_1, pdf_params%rt_2,                     & ! Out
         565             :                             pdf_params%thl_1, pdf_params%thl_2,                   & ! Out
         566             :                             pdf_params%varnce_w_1, pdf_params%varnce_w_2,         & ! Out
         567             :                             pdf_params%varnce_rt_1, pdf_params%varnce_rt_2,       & ! Out
         568             :                             pdf_params%varnce_thl_1, pdf_params%varnce_thl_2,     & ! Out
         569             :                             pdf_params%mixt_frac,                                 & ! Out
         570             :                             pdf_params%alpha_rt, pdf_params%alpha_thl,            & ! Out
         571             :                             sigma_sqd_w, sclr1, sclr2,                            & ! Out
         572           0 :                             varnce_sclr1, varnce_sclr2, alpha_sclr )                ! Out
         573             :
         574             :     elseif ( iiPDF_type == iiPDF_3D_Luhar ) then ! use 3D Luhar
         575           0 :       do i = 1, ngrdcol
         576             :         call Luhar_3D_pdf_driver( nz, &
         577           0 :                            wm(i,:), rtm(i,:), thlm(i,:), wp2(i,:), rtp2(i,:), thlp2(i,:),                             & ! In
         578           0 :                            Skw(i,:), Skrt(i,:), Skthl(i,:), wprtp(i,:), wpthlp(i,:),                             & ! In
         579           0 :                            pdf_params%w_1(i,:), pdf_params%w_2(i,:),                    & ! Out
         580           0 :                            pdf_params%rt_1(i,:), pdf_params%rt_2(i,:),                  & ! Out
         581           0 :                            pdf_params%thl_1(i,:), pdf_params%thl_2(i,:),                & ! Out
         582           0 :                            pdf_params%varnce_w_1(i,:), pdf_params%varnce_w_2(i,:),      & ! Out
         583           0 :                            pdf_params%varnce_rt_1(i,:), pdf_params%varnce_rt_2(i,:),    & ! Out
         584           0 :                            pdf_params%varnce_thl_1(i,:), pdf_params%varnce_thl_2(i,:),  & ! Out
         585           0 :                            pdf_params%mixt_frac(i,:) )                                    ! Out
         586             :       end do
         587             :     elseif ( iiPDF_type == iiPDF_new ) then ! use new PDF
         588             :       call new_pdf_driver( nz, ngrdcol, wm, rtm, thlm, wp2, rtp2, thlp2, Skw, & ! In
         589             :                            wprtp, wpthlp, rtpthlp,                            & ! In
         590             :                            clubb_params,                                      & ! In
         591             :                            Skrt, Skthl,                                       & ! In/Out
         592             :                            pdf_params%w_1, pdf_params%w_2,                    & ! Out
         593             :                            pdf_params%rt_1, pdf_params%rt_2,                  & ! Out
         594             :                            pdf_params%thl_1, pdf_params%thl_2,                & ! Out
         595             :                            pdf_params%varnce_w_1, pdf_params%varnce_w_2,      & ! Out
         596             :                            pdf_params%varnce_rt_1, pdf_params%varnce_rt_2,    & ! Out
         597             :                            pdf_params%varnce_thl_1, pdf_params%varnce_thl_2,  & ! Out
         598             :                            pdf_params%mixt_frac,                              & ! Out
         599             :                            pdf_implicit_coefs_terms,                          & ! Out
         600             :                            F_w, F_rt, F_thl, min_F_w, max_F_w,                & ! Out
         601           0 :                            min_F_rt, max_F_rt, min_F_thl, max_F_thl )           ! Out
         602             :     elseif ( iiPDF_type == iiPDF_TSDADG ) then
         603           0 :       do i = 1, ngrdcol
         604             :         call tsdadg_pdf_driver( nz, &
         605           0 :                           wm(i,:), rtm(i,:), thlm(i,:), wp2(i,:), rtp2(i,:), thlp2(i,:),                                & ! In
         606           0 :                           Skw(i,:), Skrt(i,:), Skthl(i,:), wprtp(i,:), wpthlp(i,:),                                & ! In
         607           0 :                           pdf_params%w_1(i,:), pdf_params%w_2(i,:),                       & ! Out
         608           0 :                           pdf_params%rt_1(i,:), pdf_params%rt_2(i,:),                     & ! Out
         609           0 :                           pdf_params%thl_1(i,:), pdf_params%thl_2(i,:),                   & ! Out
         610           0 :                           pdf_params%varnce_w_1(i,:), pdf_params%varnce_w_2(i,:),         & ! Out
         611           0 :                           pdf_params%varnce_rt_1(i,:), pdf_params%varnce_rt_2(i,:),       & ! Out
         612           0 :                           pdf_params%varnce_thl_1(i,:), pdf_params%varnce_thl_2(i,:),     & ! Out
         613           0 :                           pdf_params%mixt_frac(i,:) )                                       ! Out
         614             :       end do
         615             :     elseif ( iiPDF_type == iiPDF_LY93 ) then ! use LY93
         616           0 :       do i = 1, ngrdcol
         617           0 :         call LY93_driver( nz, wm(i,:), rtm(i,:), thlm(i,:), wp2(i,:), rtp2(i,:),                          & ! In
         618           0 :                           thlp2(i,:), Skw(i,:), Skrt(i,:), Skthl(i,:),                           & ! In
         619           0 :                           pdf_params%w_1(i,:), pdf_params%w_2(i,:),                    & ! Out
         620           0 :                           pdf_params%rt_1(i,:), pdf_params%rt_2(i,:),                  & ! Out
         621           0 :                           pdf_params%thl_1(i,:), pdf_params%thl_2(i,:),                & ! Out
         622           0 :                           pdf_params%varnce_w_1(i,:), pdf_params%varnce_w_2(i,:),      & ! Out
         623           0 :                           pdf_params%varnce_rt_1(i,:), pdf_params%varnce_rt_2(i,:),    & ! Out
         624           0 :                           pdf_params%varnce_thl_1(i,:), pdf_params%varnce_thl_2(i,:),  & ! Out
         625           0 :                           pdf_params%mixt_frac(i,:) )                               ! Out
         626             :       end do
         627             :     elseif ( iiPDF_type == iiPDF_new_hybrid ) then ! use new hybrid PDF
         628             :       call new_hybrid_pdf_driver( nz, ngrdcol, sclr_dim,                          & ! In
         629             :                                   wm, rtm, thlm, um, vm,                          & ! In
         630             :                                   wp2, rtp2, thlp2, up2, vp2,                     & ! In
         631             :                                   Skw, wprtp, wpthlp, upwp, vpwp,                 & ! In
         632             :                                   sclrm, sclrp2, wpsclrp,                         & ! In
         633             :                                   gamma_Skw_fnc,                                  & ! In
         634             :                                   clubb_params(:,islope_coef_spread_DG_means_w),  & ! In
         635             :                                   clubb_params(:,ipdf_component_stdev_factor_w),  & ! In
         636             :                                   Skrt, Skthl, Sku, Skv, Sksclr,                  & ! I/O
         637             :                                   pdf_params%w_1, pdf_params%w_2,                 & ! Out
         638             :                                   pdf_params%rt_1, pdf_params%rt_2,               & ! Out
         639             :                                   pdf_params%thl_1, pdf_params%thl_2,             & ! Out
         640             :                                   u_1, u_2, v_1, v_2,                             & ! Out
         641             :                                   pdf_params%varnce_w_1,                          & ! Out
         642             :                                   pdf_params%varnce_w_2,                          & ! Out
         643             :                                   pdf_params%varnce_rt_1,                         & ! Out
         644             :                                   pdf_params%varnce_rt_2,                         & ! Out
         645             :                                   pdf_params%varnce_thl_1,                        & ! Out
         646             :                                   pdf_params%varnce_thl_2,                        & ! Out
         647             :                                   varnce_u_1, varnce_u_2,                         & ! Out
         648             :                                   varnce_v_1, varnce_v_2,                         & ! Out
         649             :                                   sclr1, sclr2,                                   & ! Out
         650             :                                   varnce_sclr1, varnce_sclr2,                     & ! Out
         651             :                                   pdf_params%mixt_frac,                           & ! Out
         652             :                                   pdf_implicit_coefs_terms,                       & ! Out
         653           0 :                                   F_w, min_F_w, max_F_w )                           ! Out
         654             :
         655             :       ! The calculation of skewness of rt, thl, u, v, and scalars is hard-wired
         656             :       ! for use with the ADG1 code, which contains the variable sigma_sqd_w.
         657             :       ! In order to use an equivalent expression for these skewnesses using the
         658             :       ! new hybrid PDF (without doing more recoding), set the value of
         659             :       ! sigma_sqd_w to 1 - F_w.
         660           0 :       do k = 1, nz
         661           0 :         do i = 1, ngrdcol
         662           0 :           sigma_sqd_w(i,k) = one - F_w(i,k)
         663             :         end do
         664             :       end do
         665             :
         666             :     end if ! iiPDF_type
         667             :
         668             :     ! Calculate the PDF component correlations of rt and thl.
         669             :     call calc_comp_corrs_binormal( nz, ngrdcol,                                           & ! In
         670             :                                    rtpthlp, rtm, thlm,                                    & ! In
         671             :                                    pdf_params%rt_1, pdf_params%rt_2,                      & ! In
         672             :                                    pdf_params%thl_1, pdf_params%thl_2,                    & ! In
         673             :                                    pdf_params%varnce_rt_1, pdf_params%varnce_rt_2,        & ! In
         674             :                                    pdf_params%varnce_thl_1, pdf_params%varnce_thl_2,      & ! In
         675             :                                    pdf_params%mixt_frac,                                  & ! In
         676    17870112 :                                    pdf_params%corr_rt_thl_1, pdf_params%corr_rt_thl_2 )     ! Out
         677             :
         678             :     if ( iiPDF_type == iiPDF_ADG1 .or. iiPDF_type == iiPDF_ADG2 &
         679    17870112 :          .or. iiPDF_type == iiPDF_new_hybrid ) then
         680             :
         681             :       ! These PDF types define corr_w_rt_1, corr_w_rt_2, corr_w_thl_1, and
         682             :       ! corr_w_thl_2 to all have a value of 0, so skip the calculation.
         683             :       ! The values of corr_u_w_1, corr_u_w_2, corr_v_w_1, and corr_v_w_2 are
         684             :       ! all defined to be 0, as well.
         685             :       !$acc parallel loop gang vector collapse(2) default(present)
         686  1536829632 :       do k = 1, nz
         687 25380961632 :         do i = 1, ngrdcol
         688 23844132000 :           pdf_params%corr_w_rt_1(i,k)  = zero
         689 23844132000 :           pdf_params%corr_w_rt_2(i,k)  = zero
         690 23844132000 :           pdf_params%corr_w_thl_1(i,k) = zero
         691 23844132000 :           pdf_params%corr_w_thl_2(i,k) = zero
         692 23844132000 :           corr_u_w_1(i,k)   = zero
         693 23844132000 :           corr_u_w_2(i,k)   = zero
         694 23844132000 :           corr_v_w_1(i,k)   = zero
         695 25363091520 :           corr_v_w_2(i,k)   = zero
         696             :         end do
         697             :       end do
         698             :       !$acc end parallel loop
         699             :
         700             :     else
         701             :
         702             :       ! Calculate the PDF component correlations of w and rt.
         703             :       call calc_comp_corrs_binormal( nz, ngrdcol,                                      & ! In
         704             :                                      wprtp, wm, rtm,                                   & ! In
         705             :                                      pdf_params%w_1, pdf_params%w_2,                   & ! In
         706             :                                      pdf_params%rt_1, pdf_params%rt_2,                 & ! In
         707             :                                      pdf_params%varnce_w_1, pdf_params%varnce_w_2,     & ! In
         708             :                                      pdf_params%varnce_rt_1, pdf_params%varnce_rt_2,   & ! In
         709             :                                      pdf_params%mixt_frac,                             & ! In
         710           0 :                                      pdf_params%corr_w_rt_1, pdf_params%corr_w_rt_2 )    ! Out
         711             :
         712             :       ! Calculate the PDF component correlations of w and thl.
         713             :       call calc_comp_corrs_binormal( nz, ngrdcol,                                          & ! In
         714             :                                      wpthlp, wm, thlm,                                     & ! In
         715             :                                      pdf_params%w_1, pdf_params%w_2,                       & ! In
         716             :                                      pdf_params%thl_1, pdf_params%thl_2,                   & ! In
         717             :                                      pdf_params%varnce_w_1, pdf_params%varnce_w_2,         & ! In
         718             :                                      pdf_params%varnce_thl_1, pdf_params%varnce_thl_2,     & ! In
         719             :                                      pdf_params%mixt_frac,                                 & ! In
         720           0 :                                      pdf_params%corr_w_thl_1, pdf_params%corr_w_thl_2 )      ! Out
         721             :     end if
         722             :
         723    17870112 :     if ( l_scalar_calc ) then
         724             :
         725             :       ! Calculate the PDF component correlations of a passive scalar and thl.
         726           0 :       do j = 1, sclr_dim
         727             :         call calc_comp_corrs_binormal( nz, ngrdcol,                                       & ! In
         728             :                                        sclrpthlp(:,:,j), sclrm(:,:,j), thlm,              & ! In
         729             :                                        sclr1(:,:,j), sclr2(:,:,j),                        & ! In
         730             :                                        pdf_params%thl_1, pdf_params%thl_2,                & ! In
         731             :                                        varnce_sclr1(:,:,j), varnce_sclr2(:,:,j),          & ! In
         732             :                                        pdf_params%varnce_thl_1, pdf_params%varnce_thl_2,  & ! In
         733             :                                        pdf_params%mixt_frac,                              & ! In
         734           0 :                                        corr_sclr_thl_1(:,:,j), corr_sclr_thl_2(:,:,j) )     ! Out
         735             :       end do
         736             :
         737             :       ! Calculate the PDF component correlations of a passive scalar and rt.
         738           0 :       do j = 1, sclr_dim
         739             :         call calc_comp_corrs_binormal( nz, ngrdcol,                                       & ! In
         740             :                                        sclrprtp(:,:,j), sclrm(:,:,j), rtm,                & ! In
         741             :                                        sclr1(:,:,j), sclr2(:,:,j),                        & ! In
         742             :                                        pdf_params%rt_1, pdf_params%rt_2,                  & ! In
         743             :                                        varnce_sclr1(:,:,j), varnce_sclr2(:,:,j),          & ! In
         744             :                                        pdf_params%varnce_rt_1, pdf_params%varnce_rt_2,    & ! In
         745             :                                        pdf_params%mixt_frac,                              & ! In
         746           0 :                                        corr_sclr_rt_1(:,:,j), corr_sclr_rt_2(:,:,j) )       ! Out
         747             :       end do
         748             :
         749             :       if ( iiPDF_type == iiPDF_ADG1 .or. iiPDF_type == iiPDF_ADG2 &
         750           0 :            .or. iiPDF_type == iiPDF_new_hybrid ) then
         751             :
         752             :         ! These PDF types define all PDF component correlations involving w
         753             :         ! to have a value of 0, so skip the calculation.
         754             :         !$acc parallel loop gang vector collapse(2) default(present)
         755           0 :         do j = 1, sclr_dim
         756           0 :           do k = 1, nz
         757           0 :             do i = 1, ngrdcol
         758           0 :               corr_w_sclr_1(i,k,j) = zero
         759           0 :               corr_w_sclr_2(i,k,j) = zero
         760             :             end do
         761             :           end do
         762             :         end do
         763             :         !$acc end parallel loop
         764             :
         765             :       else
         766             :
         767             :         ! Calculate the PDF component correlations of w and a passive scalar.
         768           0 :         do j = 1, sclr_dim
         769             :           call calc_comp_corrs_binormal( nz, ngrdcol,                                   & ! In
         770             :                                          wpsclrp(:,:,j), wm, sclrm(:,:,j),              & ! In
         771             :                                          pdf_params%w_1, pdf_params%w_2,                & ! In
         772             :                                          sclr1(:,:,j), sclr2(:,:,j),                    & ! In
         773             :                                          pdf_params%varnce_w_1, pdf_params%varnce_w_2,  & ! In
         774             :                                          varnce_sclr1(:,:,j), varnce_sclr2(:,:,j),      & ! In
         775             :                                          pdf_params%mixt_frac,                          & ! In
         776           0 :                                          corr_w_sclr_1(:,:,j), corr_w_sclr_2(:,:,j) )     ! Out
         777             :
         778             :         end do
         779             :
         780             :       end if
         781             :
         782             :     end if
         783             :
         784             :
         785             :     ! Compute higher order moments (these are interactive)
         786             :     call calc_wp2xp_pdf( nz, ngrdcol,                                       &
         787             :                          wm, rtm, pdf_params%w_1, pdf_params%w_2,           &
         788             :                          pdf_params%rt_1, pdf_params%rt_2,                  &
         789             :                          pdf_params%varnce_w_1, pdf_params%varnce_w_2,      &
         790             :                          pdf_params%varnce_rt_1, pdf_params%varnce_rt_2,    &
         791             :                          pdf_params%corr_w_rt_1, pdf_params%corr_w_rt_2,    &
         792             :                          pdf_params%mixt_frac,                              &
         793    17870112 :                          wp2rtp )
         794             :
         795             :     call calc_wp2xp_pdf( nz, ngrdcol,                                          &
         796             :                          wm, thlm, pdf_params%w_1, pdf_params%w_2,             &
         797             :                          pdf_params%thl_1, pdf_params%thl_2,                   &
         798             :                          pdf_params%varnce_w_1, pdf_params%varnce_w_2,         &
         799             :                          pdf_params%varnce_thl_1, pdf_params%varnce_thl_2,     &
         800             :                          pdf_params%corr_w_thl_1, pdf_params%corr_w_thl_2,     &
         801             :                          pdf_params%mixt_frac,                                 &
         802    17870112 :                          wp2thlp )
         803             :
         804             :     ! Compute higher order moments (these may be interactive)
         805             :     call calc_wpxp2_pdf( nz, ngrdcol, &
         806             :                          wm, um, pdf_params%w_1, pdf_params%w_2, &
         807             :                          u_1, u_2, &
         808             :                          pdf_params%varnce_w_1, pdf_params%varnce_w_2,   &
         809             :                          varnce_u_1, varnce_u_2, &
         810             :                          corr_u_w_1, corr_u_w_2, &
         811             :                          pdf_params%mixt_frac, &
         812    17870112 :                          wpup2 )
         813             :
         814             :     call calc_wpxp2_pdf( nz, ngrdcol, &
         815             :                          wm, vm, pdf_params%w_1, pdf_params%w_2, &
         816             :                          v_1, v_2, &
         817             :                          pdf_params%varnce_w_1, pdf_params%varnce_w_2,   &
         818             :                          varnce_v_1, varnce_v_2, &
         819             :                          corr_v_w_1, corr_v_w_2, &
         820             :                          pdf_params%mixt_frac, &
         821    17870112 :                          wpvp2 )
         822             :
         823             :     call calc_wp2xp2_pdf( nz, ngrdcol, &
         824             :                           wm, um, pdf_params%w_1, pdf_params%w_2, &
         825             :                           u_1, u_2, &
         826             :                           pdf_params%varnce_w_1, pdf_params%varnce_w_2, &
         827             :                           varnce_u_1, varnce_u_2, &
         828             :                           corr_u_w_1, corr_u_w_2, &
         829             :                           pdf_params%mixt_frac, &
         830    17870112 :                           wp2up2 )
         831             :
         832             :     call calc_wp2xp2_pdf( nz, ngrdcol, &
         833             :                           wm, vm, pdf_params%w_1, pdf_params%w_2, &
         834             :                           v_1, v_2, &
         835             :                           pdf_params%varnce_w_1, pdf_params%varnce_w_2, &
         836             :                           varnce_v_1, varnce_v_2, &
         837             :                           corr_v_w_1, corr_v_w_2, &
         838             :                           pdf_params%mixt_frac, &
         839    17870112 :                           wp2vp2 )
         840             :
         841             :     call calc_wp4_pdf( nz, ngrdcol, &
         842             :                        wm, pdf_params%w_1, pdf_params%w_2, &
         843             :                        pdf_params%varnce_w_1, pdf_params%varnce_w_2,    &
         844             :                        pdf_params%mixt_frac, &
         845    17870112 :                        wp4 )
         846             :
         847    17870112 :     if ( l_explicit_turbulent_adv_xpyp .or. stats_metadata%iwprtp2 > 0 ) then
         848             :       call calc_wpxp2_pdf( nz, ngrdcol, &
         849             :                            wm, rtm, pdf_params%w_1, pdf_params%w_2,        &
         850             :                            pdf_params%rt_1, pdf_params%rt_2,               &
         851             :                            pdf_params%varnce_w_1, pdf_params%varnce_w_2,   &
         852             :                            pdf_params%varnce_rt_1, pdf_params%varnce_rt_2, &
         853             :                            pdf_params%corr_w_rt_1, pdf_params%corr_w_rt_2, &
         854             :                            pdf_params%mixt_frac, &
         855           0 :                            wprtp2 )
         856             :     end if
         857             :
         858    17870112 :     if ( l_explicit_turbulent_adv_xpyp .or. stats_metadata%iwpthlp2 > 0 ) then
         859             :       call calc_wpxp2_pdf( nz, ngrdcol, &
         860             :                            wm, thlm, pdf_params%w_1, pdf_params%w_2,          &
         861             :                            pdf_params%thl_1, pdf_params%thl_2,                &
         862             :                            pdf_params%varnce_w_1, pdf_params%varnce_w_2,      &
         863             :                            pdf_params%varnce_thl_1, pdf_params%varnce_thl_2,  &
         864             :                            pdf_params%corr_w_thl_1, pdf_params%corr_w_thl_2,  &
         865             :                            pdf_params%mixt_frac, &
         866           0 :                            wpthlp2 )
         867             :     end if
         868             :
         869    17870112 :     if ( l_explicit_turbulent_adv_xpyp .or. stats_metadata%iwprtpthlp > 0 ) then
         870             :
         871             :       call calc_wpxpyp_pdf( nz, ngrdcol, &
         872             :                             wm, rtm, thlm, pdf_params%w_1, pdf_params%w_2,      &
         873             :                             pdf_params%rt_1, pdf_params%rt_2,                   &
         874             :                             pdf_params%thl_1, pdf_params%thl_2,                 &
         875             :                             pdf_params%varnce_w_1, pdf_params%varnce_w_2,       &
         876             :                             pdf_params%varnce_rt_1, pdf_params%varnce_rt_2,     &
         877             :                             pdf_params%varnce_thl_1, pdf_params%varnce_thl_2,   &
         878             :                             pdf_params%corr_w_rt_1, pdf_params%corr_w_rt_2,     &
         879             :                             pdf_params%corr_w_thl_1, pdf_params%corr_w_thl_2,   &
         880             :                             pdf_params%corr_rt_thl_1, pdf_params%corr_rt_thl_2, &
         881             :                             pdf_params%mixt_frac, &
         882           0 :                             wprtpthlp )
         883             :     end if
         884             :
         885             :
         886             :     ! Scalar Addition to higher order moments
         887    17870112 :     if ( l_scalar_calc ) then
         888             :
         889           0 :       do j = 1, sclr_dim
         890             :         call calc_wp2xp_pdf( nz, ngrdcol,                                       &
         891             :                              wm, sclrm(:,:,j), pdf_params%w_1, pdf_params%w_2,  &
         892             :                              sclr1(:,:,j), sclr2(:,:,j),                        &
         893             :                              pdf_params%varnce_w_1, pdf_params%varnce_w_2,      &
         894             :                              varnce_sclr1(:,:,j), varnce_sclr2(:,:,j),          &
         895             :                              corr_w_sclr_1(:,:,j), corr_w_sclr_2(:,:,j),        &
         896             :                              pdf_params%mixt_frac,                              &
         897           0 :                              wp2sclrp(:,:,j) )
         898             :       end do
         899             :
         900           0 :       do j = 1, sclr_dim
         901             :         call calc_wpxp2_pdf( nz, ngrdcol, &
         902             :                              wm, sclrm(:,:,j), pdf_params%w_1, pdf_params%w_2, &
         903             :                              sclr1(:,:,j), sclr2(:,:,j),                         &
         904             :                              pdf_params%varnce_w_1, pdf_params%varnce_w_2,   &
         905             :                              varnce_sclr1(:,:,j), varnce_sclr2(:,:,j),           &
         906             :                              corr_w_sclr_1(:,:,j), corr_w_sclr_2(:,:,j),         &
         907             :                              pdf_params%mixt_frac, &
         908           0 :                              wpsclrp2(:,:,j) )
         909             :       end do
         910             :
         911           0 :       do j = 1, sclr_dim
         912             :         call calc_wpxpyp_pdf( nz, ngrdcol, &
         913             :                               wm, sclrm(:,:,j), rtm, pdf_params%w_1, pdf_params%w_2,  &
         914             :                               sclr1(:,:,j), sclr2(:,:,j),                             &
         915             :                               pdf_params%rt_1, pdf_params%rt_2,                       &
         916             :                               pdf_params%varnce_w_1, pdf_params%varnce_w_2,           &
         917             :                               varnce_sclr1(:,:,j), varnce_sclr2(:,:,j),               &
         918             :                               pdf_params%varnce_rt_1, pdf_params%varnce_rt_2,         &
         919             :                               corr_w_sclr_1(:,:,j), corr_w_sclr_2(:,:,j),             &
         920             :                               pdf_params%corr_w_rt_1, pdf_params%corr_w_rt_2,         &
         921             :                               corr_sclr_rt_1(:,:,j), corr_sclr_rt_2(:,:,j),           &
         922             :                               pdf_params%mixt_frac, &
         923           0 :                               wpsclrprtp(:,:,j) )
         924             :       end do
         925             :
         926           0 :       do j = 1, sclr_dim
         927             :         call calc_wpxpyp_pdf( nz, ngrdcol, &
         928             :                               wm, sclrm(:,:,j), thlm, pdf_params%w_1, pdf_params%w_2,   &
         929             :                               sclr1(:,:,j), sclr2(:,:,j),                               &
         930             :                               pdf_params%thl_1, pdf_params%thl_2,                       &
         931             :                               pdf_params%varnce_w_1, pdf_params%varnce_w_2,             &
         932             :                               varnce_sclr1(:,:,j), varnce_sclr2(:,:,j),                 &
         933             :                               pdf_params%varnce_thl_1, pdf_params%varnce_thl_2,         &
         934             :                               corr_w_sclr_1(:,:,j), corr_w_sclr_2(:,:,j),               &
         935             :                               pdf_params%corr_w_thl_1, pdf_params%corr_w_thl_2,         &
         936             :                               corr_sclr_thl_1(:,:,j), corr_sclr_thl_2(:,:,j),           &
         937             :                               pdf_params%mixt_frac, &
         938           0 :                               wpsclrpthlp(:,:,j) )
         939             :       end do
         940             :
         941             :     end if
         942             :
         943             :     ! Compute higher order moments that include theta_v.
         944             :
         945             :     ! First compute some preliminary quantities.
         946             :     ! "1" denotes first Gaussian; "2" denotes 2nd Gaussian
         947             :     ! liq water temp (Sommeria & Deardorff 1977 (SD), eqn. 3)
         948             :     !$acc parallel loop gang vector collapse(2) default(present)
         949  1536829632 :     do k = 1, nz
         950 25380961632 :       do i = 1, ngrdcol
         951 23844132000 :         tl1(i,k)  = pdf_params%thl_1(i,k)*exner(i,k)
         952 25363091520 :         tl2(i,k)  = pdf_params%thl_2(i,k)*exner(i,k)
         953             :       end do
         954             :     end do
         955             :     !$acc end parallel loop
         956             :
         957             : #ifdef GFDL
         958             :     if ( sclr_dim > 0  .and.  (.not. do_liquid_only_in_clubb) ) then ! h1g, 2010-06-16 begin mod
         959             :
         960             :       do i = 1, ngrdcol
         961             :         where ( tl1(i,:) > t1_combined )
         962             :           pdf_params%rsatl_1(i,:) = sat_mixrat_liq( nz, p_in_Pa(i,:), tl1(i,:) )
         963             :         elsewhere ( tl1(i,:) > t2_combined )
         964             :           pdf_params%rsatl_1(i,:) = sat_mixrat_liq( nz, p_in_Pa(i,:), tl1(i,:) ) &
         965             :                     * (tl1(i,:) - t2_combined)/(t1_combined - t2_combined) &
         966             :                     + sat_mixrat_ice( nz, p_in_Pa(i,:), tl1(i,:) ) &
         967             :                       * (t1_combined - tl1(i,:))/(t1_combined - t2_combined)
         968             :         elsewhere ( tl1(i,:) > t3_combined )
         969             :           pdf_params%rsatl_1(i,:) = sat_mixrat_ice( nz, p_in_Pa(i,:), tl1(i,:) ) &
         970             :                     + sat_mixrat_ice( nz, p_in_Pa(i,:), tl1(i,:) ) * (RH_crit(i, :, 1, 1) -one ) &
         971             :                       * ( t2_combined -tl1(i,:))/(t2_combined - t3_combined)
         972             :         elsewhere
         973             :           pdf_params%rsatl_1(i,:) = sat_mixrat_ice( nz, p_in_Pa(i,:), tl1(i,:) ) * RH_crit(i, :, 1, 1)
         974             :         endwhere
         975             :
         976             :         where ( tl2(i,:) > t1_combined )
         977             :           pdf_params%rsatl_2(i,:) = sat_mixrat_liq( nz, p_in_Pa(i,:), tl2(i,:) )
         978             :         elsewhere ( tl2(i,:) > t2_combined )
         979             :           pdf_params%rsatl_2(i,:) = sat_mixrat_liq( nz, p_in_Pa(i,:), tl2(i,:) ) &
         980             :                     * (tl2(i,:) - t2_combined)/(t1_combined - t2_combined) &
         981             :                     + sat_mixrat_ice( nz, p_in_Pa(i,:), tl2(i,:) ) &
         982             :                       * (t1_combined - tl2(i,:))/(t1_combined - t2_combined)
         983             :         elsewhere ( tl2(i,:) > t3_combined )
         984             :           pdf_params%rsatl_2(i,:) = sat_mixrat_ice( nz, p_in_Pa(i,:), tl2(i,:) ) &
         985             :                     + sat_mixrat_ice( nz, p_in_Pa(i,:), tl2(i,:) )* (RH_crit(i, :, 1, 2) -one) &
         986             :                       * ( t2_combined -tl2(i,:))/(t2_combined - t3_combined)
         987             :         elsewhere
         988             :           pdf_params%rsatl_2(i,:) = sat_mixrat_ice( nz, p_in_Pa(i,:), tl2(i,:) ) * RH_crit(i, :, 1, 2)
         989             :         endwhere
         990             :
         991             :       end do
         992             :
         993             :     else ! sclr_dim <= 0  or  do_liquid_only_in_clubb = .T.
         994             :
         995             :       pdf_params%rsatl_1 = sat_mixrat_liq( nz, ngrdcol, p_in_Pa, tl1, saturation_formula )
         996             :       pdf_params%rsatl_2 = sat_mixrat_liq( nz, ngrdcol, p_in_Pa, tl2, saturation_formula )
         997             :
         998             :     end if !sclr_dim > 0
         999             :
        1000             :     ! Determine whether to compute ice_supersat_frac. We do not compute
        1001             :     ! ice_supersat_frac for GFDL (unless do_liquid_only_in_clubb is true),
        1002             :     ! because liquid and ice are both fed into rtm, ruining the calculation.
        1003             :     if (do_liquid_only_in_clubb) then
        1004             :       l_calc_ice_supersat_frac = .true.
        1005             :     else
        1006             :       l_calc_ice_supersat_frac = .false.
        1007             :     end if
        1008             :
        1009             : #else
        1010    17870112 :     rsatl_1 = sat_mixrat_liq( nz, ngrdcol, p_in_Pa, tl1, saturation_formula  )
        1011    17870112 :     rsatl_2 = sat_mixrat_liq( nz, ngrdcol, p_in_Pa, tl2, saturation_formula  ) ! h1g, 2010-06-16 end mod
        1012             :
        1013             :     !$acc parallel loop gang vector collapse(2) default(present)
        1014  1536829632 :     do k = 1, nz
        1015 25380961632 :       do i = 1, ngrdcol
        1016 23844132000 :         pdf_params%rsatl_1(i,k) = rsatl_1(i,k)
        1017 25363091520 :         pdf_params%rsatl_2(i,k) = rsatl_2(i,k)
        1018             :       end do
        1019             :     end do
        1020             :     !$acc end parallel loop
        1021             :
        1022    17870112 :     l_calc_ice_supersat_frac = .true.
        1023             : #endif
        1024             :
        1025             :     call transform_pdf_chi_eta_component( nz, ngrdcol, &
        1026             :                                           tl1, pdf_params%rsatl_1, pdf_params%rt_1, exner,  & ! In
        1027             :                                           pdf_params%varnce_thl_1, pdf_params%varnce_rt_1,  & ! In
        1028             :                                           pdf_params%corr_rt_thl_1, pdf_params%chi_1,       & ! In
        1029             :                                           pdf_params%crt_1, pdf_params%cthl_1,              & ! Out
        1030             :                                           pdf_params%stdev_chi_1, pdf_params%stdev_eta_1,   & ! Out
        1031             :                                           pdf_params%covar_chi_eta_1,                       & ! Out
        1032    17870112 :                                           pdf_params%corr_chi_eta_1 )                         ! Out
        1033             :
        1034             :
        1035             :     ! Calculate cloud fraction component for pdf 1
        1036             :     call calc_liquid_cloud_frac_component( nz, ngrdcol, &
        1037             :                                            pdf_params%chi_1, pdf_params%stdev_chi_1,    & ! In
        1038    17870112 :                                            pdf_params%cloud_frac_1, pdf_params%rc_1 )     ! Out
        1039             :
        1040             :     ! Calc ice_supersat_frac
        1041             :     if ( l_calc_ice_supersat_frac ) then
        1042             :
        1043             :       call calc_ice_cloud_frac_component( nz, ngrdcol, &
        1044             :                                           pdf_params%chi_1, pdf_params%stdev_chi_1, &
        1045             :                                           pdf_params%rc_1, pdf_params%cloud_frac_1, &
        1046             :                                           p_in_Pa, tl1, &
        1047             :                                           pdf_params%rsatl_1, pdf_params%crt_1, &
        1048             :                                           saturation_formula, &
        1049    17870112 :                                           pdf_params%ice_supersat_frac_1, rc_1_ice )
        1050             :     end if
        1051             :
        1052             :     call transform_pdf_chi_eta_component( nz, ngrdcol, &
        1053             :                                           tl2, pdf_params%rsatl_2, pdf_params%rt_2, exner,  & ! In
        1054             :                                           pdf_params%varnce_thl_2, pdf_params%varnce_rt_2,  & ! In
        1055             :                                           pdf_params%corr_rt_thl_2, pdf_params%chi_2,       & ! In
        1056             :                                           pdf_params%crt_2, pdf_params%cthl_2,              & ! Out
        1057             :                                           pdf_params%stdev_chi_2, pdf_params%stdev_eta_2,   & ! Out
        1058             :                                           pdf_params%covar_chi_eta_2,                       & ! Out
        1059    17870112 :                                           pdf_params%corr_chi_eta_2 )                         ! Out
        1060             :
        1061             :
        1062             :     ! Calculate cloud fraction component for pdf 2
        1063             :     call calc_liquid_cloud_frac_component( nz, ngrdcol, &
        1064             :                                            pdf_params%chi_2, pdf_params%stdev_chi_2,    & ! In
        1065    17870112 :                                            pdf_params%cloud_frac_2, pdf_params%rc_2 )     ! Out
        1066             :
        1067             :     ! Calc ice_supersat_frac
        1068             :     if ( l_calc_ice_supersat_frac ) then
        1069             :
        1070             :       call calc_ice_cloud_frac_component( nz, ngrdcol, &
        1071             :                                           pdf_params%chi_2, pdf_params%stdev_chi_2, &
        1072             :                                           pdf_params%rc_2, pdf_params%cloud_frac_2, &
        1073             :                                           p_in_Pa, tl2, &
        1074             :                                           pdf_params%rsatl_2, pdf_params%crt_2, &
        1075             :                                           saturation_formula, &
        1076    17870112 :                                           pdf_params%ice_supersat_frac_2, rc_2_ice )
        1077             :
        1078             :       ! Compute ice cloud fraction, ice_supersat_frac
        1079             :       !$acc parallel loop gang vector collapse(2) default(present)
        1080  1536829632 :       do k = 1, nz
        1081 25380961632 :         do i = 1, ngrdcol
        1082 47688264000 :           ice_supersat_frac(i,k) = pdf_params%mixt_frac(i,k) &
        1083           0 :                                    * pdf_params%ice_supersat_frac_1(i,k) &
        1084             :                                    + ( one - pdf_params%mixt_frac(i,k) ) &
        1085 73051355520 :                                      * pdf_params%ice_supersat_frac_2(i,k)
        1086             :         end do
        1087             :       end do
        1088             :       !$acc end parallel loop
        1089             :
        1090             :     else
        1091             :
        1092             :       ! ice_supersat_frac will be garbage if computed as above
        1093             :       !$acc parallel loop gang vector collapse(2) default(present)
        1094             :       do k = 1, nz
        1095             :         do i = 1, ngrdcol
        1096             :           ice_supersat_frac(i,k) = 0.0_core_rknd
        1097             :         end do
        1098             :       end do
        1099             :       !$acc end parallel loop
        1100             :
        1101             :       if (clubb_at_least_debug_level( 1 )) then
        1102             :           write(fstderr,*) "Warning: ice_supersat_frac has garbage values if &
        1103             :                           & do_liquid_only_in_clubb = .false."
        1104             :       end if
        1105             :
        1106             :     end if ! l_calc_ice_supersat_frac
        1107             :
        1108             :
        1109             :     ! Compute cloud fraction and mean cloud water mixing ratio.
        1110             :     ! Reference:
        1111             :     ! https://arxiv.org/pdf/1711.03675v1.pdf#nameddest=url:anl_int_cloud_terms
        1112             :     !$acc parallel loop gang vector collapse(2) default(present)
        1113  1536829632 :     do k = 1, nz
        1114 25380961632 :       do i = 1, ngrdcol
        1115 47688264000 :         cloud_frac(i,k) = pdf_params%mixt_frac(i,k) * pdf_params%cloud_frac_1(i,k) &
        1116 71532396000 :                      + ( one - pdf_params%mixt_frac(i,k) ) * pdf_params%cloud_frac_2(i,k)
        1117           0 :         rcm(i,k) = pdf_params%mixt_frac(i,k) * pdf_params%rc_1(i,k) + ( one - pdf_params%mixt_frac(i,k) ) &
        1118 23844132000 :                                                                  * pdf_params%rc_2(i,k)
        1119 25363091520 :         rcm(i,k) = max( zero_threshold, rcm(i,k) )
        1120             :       end do
        1121             :     end do
        1122             :     !$acc end parallel loop
        1123             :
        1124             :     if ( iiPDF_type == iiPDF_ADG1 .or. iiPDF_type == iiPDF_ADG2 &
        1125    17870112 :          .or. iiPDF_type == iiPDF_new_hybrid ) then
        1126             :
        1127             :       ! corr_w_rt and corr_w_thl are zero for these pdf types so
        1128             :       ! corr_w_chi and corr_w_eta are zero as well
        1129             :       !$acc parallel loop gang vector collapse(2) default(present)
        1130  1536829632 :       do k = 1, nz
        1131 25380961632 :         do i = 1, ngrdcol
        1132 23844132000 :           pdf_params%corr_w_chi_1(i,k) = zero
        1133 23844132000 :           pdf_params%corr_w_chi_2(i,k) = zero
        1134 23844132000 :           pdf_params%corr_w_eta_1(i,k) = zero
        1135 25363091520 :           pdf_params%corr_w_eta_2(i,k) = zero
        1136             :         end do
        1137             :       end do
        1138             :       !$acc end parallel loop
        1139             :
        1140             :     else
        1141             :
        1142             :       ! Correlation of w and chi for each component.
        1143             :       pdf_params%corr_w_chi_1 &
        1144           0 :       = calc_corr_chi_x( pdf_params%crt_1, pdf_params%cthl_1, &
        1145           0 :                          sqrt(pdf_params%varnce_rt_1), sqrt(pdf_params%varnce_thl_1), &
        1146           0 :                          pdf_params%stdev_chi_1, &
        1147           0 :                          pdf_params%corr_w_rt_1, pdf_params%corr_w_thl_1 )
        1148             :
        1149             :       pdf_params%corr_w_chi_2 &
        1150           0 :       = calc_corr_chi_x( pdf_params%crt_2, pdf_params%cthl_2, &
        1151           0 :                          sqrt(pdf_params%varnce_rt_2), sqrt(pdf_params%varnce_thl_2), &
        1152           0 :                          pdf_params%stdev_chi_2, pdf_params%corr_w_rt_2, &
        1153           0 :                          pdf_params%corr_w_thl_2 )
        1154             :
        1155             :       ! Correlation of w and eta for each component.
        1156             :       pdf_params%corr_w_eta_1 &
        1157           0 :       = calc_corr_eta_x( pdf_params%crt_1, pdf_params%cthl_1, &
        1158           0 :                          sqrt(pdf_params%varnce_rt_1), sqrt(pdf_params%varnce_thl_1), &
        1159           0 :                          pdf_params%stdev_eta_1, pdf_params%corr_w_rt_1, &
        1160           0 :                          pdf_params%corr_w_thl_1 )
        1161             :
        1162             :       pdf_params%corr_w_eta_2 &
        1163           0 :       = calc_corr_eta_x( pdf_params%crt_2, pdf_params%cthl_2, &
        1164           0 :                          sqrt(pdf_params%varnce_rt_2), sqrt(pdf_params%varnce_thl_2), &
        1165           0 :                          pdf_params%stdev_eta_2, pdf_params%corr_w_rt_2, &
        1166           0 :                          pdf_params%corr_w_thl_2 )
        1167             :
        1168             :     end if
        1169             :
        1170             :
        1171             :     ! Compute moments that depend on theta_v
        1172             :     !
        1173             :     ! The moments that depend on th_v' are calculated based on an approximated
        1174             :     ! and linearized form of the theta_v equation:
        1175             :     !
        1176             :     ! theta_v = theta_l + { (R_v/R_d) - 1 } * thv_ds * r_t
        1177             :     !                   + [ {L_v/(C_p*exner)} - (R_v/R_d) * thv_ds ] * r_c;
        1178             :     !
        1179             :     ! and therefore:
        1180             :     !
        1181             :     ! th_v' = th_l' + { (R_v/R_d) - 1 } * thv_ds * r_t'
        1182             :     !               + [ {L_v/(C_p*exner)} - (R_v/R_d) * thv_ds ] * r_c';
        1183             :     !
        1184             :     ! where thv_ds is used as a reference value to approximate theta_l.
        1185             :     !
        1186             :     ! Reference:
        1187             :     ! https://arxiv.org/pdf/1711.03675v1.pdf#nameddest=url:anl_int_buoy_terms
        1188             :
        1189             :     ! Calculate the contributions to <w'rc'>, <w'^2 rc'>, <rt'rc'>, and
        1190             :     ! <thl'rc'> from the 1st PDF component.
        1191             :     call calc_xprcp_component( nz, ngrdcol,                                          & ! In
        1192             :                                wm, rtm, thlm, um, vm, rcm,                           & ! In
        1193             :                                pdf_params%w_1, pdf_params%rt_1,                      & ! In
        1194             :                                pdf_params%thl_1, u_1, v_1,                           & ! In
        1195             :                                pdf_params%varnce_w_1, pdf_params%chi_1,              & ! In
        1196             :                                pdf_params%stdev_chi_1, pdf_params%stdev_eta_1,       & ! In
        1197             :                                pdf_params%corr_w_chi_1, pdf_params%corr_chi_eta_1,   & ! In
        1198             :                                pdf_params%crt_1, pdf_params%cthl_1,                  & ! In
        1199             :                                pdf_params%rc_1, pdf_params%cloud_frac_1, iiPDF_type, & ! In
        1200             :                                wprcp_contrib_comp_1, wp2rcp_contrib_comp_1,          & ! Out
        1201             :                                rtprcp_contrib_comp_1, thlprcp_contrib_comp_1,        & ! Out
        1202    17870112 :                                uprcp_contrib_comp_1, vprcp_contrib_comp_1 )            ! Out
        1203             :
        1204             :     call calc_xprcp_component( nz, ngrdcol,                                          & ! In
        1205             :                                wm, rtm, thlm, um, vm, rcm,                           & ! In
        1206             :                                pdf_params%w_2, pdf_params%rt_2,                      & ! In
        1207             :                                pdf_params%thl_2, u_2, v_2,                           & ! In
        1208             :                                pdf_params%varnce_w_2, pdf_params%chi_2,              & ! In
        1209             :                                pdf_params%stdev_chi_2, pdf_params%stdev_eta_2,       & ! In
        1210             :                                pdf_params%corr_w_chi_2, pdf_params%corr_chi_eta_2,   & ! In
        1211             :                                pdf_params%crt_2, pdf_params%cthl_2,                  & ! In
        1212             :                                pdf_params%rc_2, pdf_params%cloud_frac_2, iiPDF_type, & ! In
        1213             :                                wprcp_contrib_comp_2, wp2rcp_contrib_comp_2,          & ! Out
        1214             :                                rtprcp_contrib_comp_2, thlprcp_contrib_comp_2,        & ! Out
        1215    17870112 :                                uprcp_contrib_comp_2, vprcp_contrib_comp_2 )            ! Out
        1216             :
        1217             :
        1218             :     ! Calculate rc_coef, which is the coefficient on <x'rc'> in the <x'thv'> equation.
        1219             :     !$acc parallel loop gang vector collapse(2) default(present)
        1220  1536829632 :     do k = 1, nz
        1221 25380961632 :       do i = 1, ngrdcol
        1222             :
        1223 23844132000 :         rc_coef(i,k) = Lv / ( exner(i,k) * Cp ) - ep2 * thv_ds(i,k)
        1224             :
        1225             :         ! Calculate <w'rc'>, <w'^2 rc'>, <rt'rc'>, and <thl'rc'>.
        1226           0 :         wprcp(i,k) = pdf_params%mixt_frac(i,k) * wprcp_contrib_comp_1(i,k) &
        1227 23844132000 :                      + ( one - pdf_params%mixt_frac(i,k) ) * wprcp_contrib_comp_2(i,k)
        1228             :
        1229           0 :         wp2rcp(i,k) = pdf_params%mixt_frac(i,k) * wp2rcp_contrib_comp_1(i,k) &
        1230 23844132000 :                       + ( one - pdf_params%mixt_frac(i,k) ) * wp2rcp_contrib_comp_2(i,k)
        1231             :
        1232           0 :         rtprcp(i,k) = pdf_params%mixt_frac(i,k) * rtprcp_contrib_comp_1(i,k) &
        1233 23844132000 :                       + ( one - pdf_params%mixt_frac(i,k) ) * rtprcp_contrib_comp_2(i,k)
        1234             :
        1235           0 :         thlprcp(i,k) = pdf_params%mixt_frac(i,k) * thlprcp_contrib_comp_1(i,k) &
        1236 23844132000 :                        + ( one - pdf_params%mixt_frac(i,k) ) * thlprcp_contrib_comp_2(i,k)
        1237             :
        1238           0 :         uprcp(i,k) = pdf_params%mixt_frac(i,k) * uprcp_contrib_comp_1(i,k) &
        1239 23844132000 :                      + ( one - pdf_params%mixt_frac(i,k) ) * uprcp_contrib_comp_2(i,k)
        1240             :
        1241           0 :         vprcp(i,k) = pdf_params%mixt_frac(i,k) * vprcp_contrib_comp_1(i,k) &
        1242 25363091520 :                      + ( one - pdf_params%mixt_frac(i,k) ) * vprcp_contrib_comp_2(i,k)
        1243             :       end do
        1244             :     end do
        1245             :     !$acc end parallel loop
        1246             :
        1247             :     ! Calculate <w'thv'>, <w'^2 thv'>, <rt'thv'>, and <thl'thv'>.
        1248             :     !$acc parallel loop gang vector collapse(2) default(present)
        1249  1536829632 :     do k = 1, nz
        1250 25380961632 :       do i = 1, ngrdcol
        1251 23844132000 :         wpthvp(i,k)  = wpthlp(i,k)  + ep1 * thv_ds(i,k) * wprtp(i,k)   + rc_coef(i,k) * wprcp(i,k)
        1252 23844132000 :         wp2thvp(i,k) = wp2thlp(i,k) + ep1 * thv_ds(i,k) * wp2rtp(i,k)  + rc_coef(i,k) * wp2rcp(i,k)
        1253 23844132000 :         rtpthvp(i,k) = rtpthlp(i,k) + ep1 * thv_ds(i,k) * rtp2(i,k)    + rc_coef(i,k) * rtprcp(i,k)
        1254 25363091520 :         thlpthvp(i,k)= thlp2(i,k)   + ep1 * thv_ds(i,k) * rtpthlp(i,k) + rc_coef(i,k) * thlprcp(i,k)
        1255             :       end do
        1256             :     end do
        1257             :     !$acc end parallel loop
        1258             :
        1259             :     ! Add the precipitation loading term in the <x'thv'> equation.
        1260             :     if ( l_liq_ice_loading_test ) then
        1261             :
        1262             :        do hm_idx = 1, hydromet_dim, 1
        1263             :
        1264             :           if ( l_mix_rat_hm(hm_idx) ) then
        1265             :             !$acc parallel loop gang vector collapse(2) default(present)
        1266             :             do k = 1, nz
        1267             :               do i = 1, ngrdcol
        1268             :                 wp2thvp(i,k)  = wp2thvp(i,k)  - thv_ds(i,k) * wp2hmp(i,k,hm_idx)
        1269             :                 wpthvp(i,k)   = wpthvp(i,k)   - thv_ds(i,k) * wphydrometp(i,k,hm_idx)
        1270             :                 thlpthvp(i,k) = thlpthvp(i,k) - thv_ds(i,k) * thlphmp(i,k,hm_idx)
        1271             :                 rtpthvp(i,k)  = rtpthvp(i,k)  - thv_ds(i,k) * rtphmp(i,k,hm_idx)
        1272             :               end do
        1273             :             end do
        1274             :             !$acc end parallel loop
        1275             :           end if
        1276             :
        1277             :        end do
        1278             :
        1279             :     end if
        1280             :
        1281             :     ! Account for subplume correlation of scalar, theta_v.
        1282             :     ! See Eqs. A13, A8 from Larson et al. (2002) ``Small-scale...''
        1283             :     !  where the ``scalar'' in this paper is w.
        1284    17870112 :     if ( l_scalar_calc ) then
        1285             :
        1286             :       !$acc parallel loop gang vector collapse(3) default(present)
        1287           0 :       do j = 1, sclr_dim
        1288           0 :         do k = 1, nz
        1289           0 :           do i = 1, ngrdcol
        1290             :
        1291           0 :             sclrprcp(i,k,j) &
        1292           0 :             = pdf_params%mixt_frac(i,k) * ( ( sclr1(i,k,j) - sclrm(i,k,j) ) * pdf_params%rc_1(i,k) ) &
        1293             :               + ( one - pdf_params%mixt_frac(i,k) ) * ( ( sclr2(i,k,j) - sclrm(i,k,j) ) &
        1294           0 :                                                         * pdf_params%rc_2(i,k) ) &
        1295           0 :               + pdf_params%mixt_frac(i,k) * corr_sclr_rt_1(i,k,j) * pdf_params%crt_1(i,k) &
        1296           0 :                 * sqrt( varnce_sclr1(i,k,j) * pdf_params%varnce_rt_1(i,k) ) &
        1297           0 :                 * pdf_params%cloud_frac_1(i,k) &
        1298           0 :               + ( one - pdf_params%mixt_frac(i,k) ) * corr_sclr_rt_2(i,k,j) * pdf_params%crt_2(i,k) &
        1299           0 :                 * sqrt( varnce_sclr2(i,k,j) * pdf_params%varnce_rt_2(i,k) ) &
        1300           0 :                 * pdf_params%cloud_frac_2(i,k) &
        1301           0 :               - pdf_params%mixt_frac(i,k) * corr_sclr_thl_1(i,k,j) * pdf_params%cthl_1(i,k) &
        1302           0 :                 * sqrt( varnce_sclr1(i,k,j) * pdf_params%varnce_thl_1(i,k) ) &
        1303             :                 * pdf_params%cloud_frac_1(i,k) &
        1304           0 :               - ( one - pdf_params%mixt_frac(i,k) ) * corr_sclr_thl_2(i,k,j) * pdf_params%cthl_2(i,k) &
        1305           0 :                 * sqrt( varnce_sclr2(i,k,j) * pdf_params%varnce_thl_2(i,k) ) &
        1306           0 :                 * pdf_params%cloud_frac_2(i,k)
        1307             :
        1308             :             sclrpthvp(i,k,j) = sclrpthlp(i,k,j) + ep1*thv_ds(i,k)*sclrprtp(i,k,j) &
        1309           0 :                              + rc_coef(i,k)*sclrprcp(i,k,j)
        1310             :
        1311             :           end do
        1312             :         end do
        1313             :       end do ! i=1, sclr_dim
        1314             :       !$acc end parallel loop
        1315             :
        1316             :     end if ! l_scalar_calc
        1317             :
        1318             :
        1319             :     ! Compute variance of liquid water mixing ratio.
        1320             :     ! This is not needed for closure.  Statistical Analysis only.
        1321             :
        1322             : #ifndef CLUBB_CAM
        1323             :       !  if CLUBB is used in CAM we want this variable computed no matter what
        1324             :       if ( stats_metadata%ircp2 > 0 ) then
        1325             : #endif
        1326             :     !$acc parallel loop gang vector collapse(2) default(present)
        1327  1536829632 :     do k = 1,nz
        1328 25380961632 :       do i = 1, ngrdcol
        1329 47688264000 :         rcp2(i,k) = pdf_params%mixt_frac(i,k) &
        1330           0 :                     * ( pdf_params%chi_1(i,k)*pdf_params%rc_1(i,k) &
        1331           0 :                         + pdf_params%cloud_frac_1(i,k)*pdf_params%stdev_chi_1(i,k)**2 ) &
        1332             :                     + ( one-pdf_params%mixt_frac(i,k) ) &
        1333           0 :                       * ( pdf_params%chi_2(i,k)*pdf_params%rc_2(i,k) &
        1334           0 :                           + pdf_params%cloud_frac_2(i,k)*pdf_params%stdev_chi_2(i,k)**2 ) &
        1335 71532396000 :                     - rcm(i,k)**2
        1336 25363091520 :         rcp2(i,k) = max( zero_threshold, rcp2(i,k) )
        1337             :
        1338             :       end do
        1339             :     end do
        1340             :     !$acc end parallel loop
        1341             : #ifndef CLUBB_CAM
        1342             :       !  if CLUBB is used in CAM we want this variable computed no matter what
        1343             :       end if
        1344             : #endif
        1345             :
        1346             :     if ( ( iiPDF_type == iiPDF_ADG1 .or. iiPDF_type == iiPDF_ADG2 &
        1347             :            .or. iiPDF_type == iiPDF_new_hybrid ) &
        1348    17870112 :          .and. ( stats_metadata%iw_up_in_cloud > 0 .or. stats_metadata%iw_down_in_cloud > 0 ) ) then
        1349             :
        1350             :       call calc_w_up_in_cloud( nz, ngrdcol, &                                      ! In
        1351             :                                pdf_params%mixt_frac, &                             ! In
        1352             :                                pdf_params%cloud_frac_1, pdf_params%cloud_frac_2, & ! In
        1353             :                                pdf_params%w_1, pdf_params%w_2, &                   ! In
        1354             :                                pdf_params%varnce_w_1, pdf_params%varnce_w_2, &     ! In
        1355             :                                w_up_in_cloud, w_down_in_cloud, &                   ! Out
        1356           0 :                                cloudy_updraft_frac, cloudy_downdraft_frac )        ! Out
        1357             :
        1358             :     else
        1359             :       !$acc parallel loop gang vector collapse(2) default(present)
        1360  1536829632 :       do k = 1,nz
        1361 25380961632 :         do i = 1, ngrdcol
        1362 23844132000 :           w_up_in_cloud(i,k) = zero
        1363 25363091520 :           w_down_in_cloud(i,k) = zero
        1364             :         end do
        1365             :       end do
        1366             :     end if
        1367             :
        1368             : #ifdef TUNER
        1369             :
        1370             :     !$acc update host( pdf_params, pdf_params%thl_1, pdf_params%thl_2 )
        1371             :
        1372             :     ! Check the first levels (and first gridcolumn) for reasonable temperatures
        1373             :     ! greater than 190K and less than 1000K
        1374             :     ! This is necessary because for certain parameter sets we can get floating point errors
        1375             :     do i=1, min( 10, size(pdf_params%thl_1(1,:)) )
        1376             :         if ( pdf_params%thl_1(1,i) < 190. ) then
        1377             :             write(fstderr,*) "Fatal error: pdf_params%thl_1 =", pdf_params%thl_1(1,i), &
        1378             :                              " < 190K at first grid column and grid level i = ", i
        1379             :             err_code = clubb_fatal_error
        1380             :             return
        1381             :         end if
        1382             :         if ( pdf_params%thl_2(1,i) < 190. ) then
        1383             :             write(fstderr,*) "Fatal error: pdf_params%thl_2 =", pdf_params%thl_2(1,i), &
        1384             :                              " < 190K at first grid column and grid level i = ", i
        1385             :             err_code = clubb_fatal_error
        1386             :             return
        1387             :         end if
        1388             :         if ( pdf_params%thl_1(1,i) > 1000. ) then
        1389             :             write(fstderr,*) "Fatal error: pdf_params%thl_1 =", pdf_params%thl_1(1,i), &
        1390             :                              " > 1000K at first grid column and grid level i = ", i
        1391             :             err_code = clubb_fatal_error
        1392             :             return
        1393             :         end if
        1394             :         if ( pdf_params%thl_2(1,i) > 1000. ) then
        1395             :             write(fstderr,*) "Fatal error: pdf_params%thl_2 =", pdf_params%thl_2(1,i), &
        1396             :                              " > 1000K at first grid column and grid level i = ", i
        1397             :             err_code = clubb_fatal_error
        1398             :             return
        1399             :         end if
        1400             :     end do
        1401             : #endif /*TUNER*/
        1402             :
        1403    17870112 :     if ( clubb_at_least_debug_level( 2 ) ) then
        1404             :
        1405             :       !$acc update host( wp4, wprtp2, wp2rtp, wpthlp2, wp2thlp, cloud_frac, &
        1406             :       !$acc              rcm, wpthvp, wp2thvp, rtpthvp, thlpthvp, wprcp, wp2rcp, &
        1407             :       !$acc              rtprcp, thlprcp, rcp2, wprtpthlp, &
        1408             :       !$acc              pdf_params%w_1, pdf_params%w_2, &
        1409             :       !$acc              pdf_params%varnce_w_1, pdf_params%varnce_w_2, &
        1410             :       !$acc              pdf_params%rt_1, pdf_params%rt_2, &
        1411             :       !$acc              pdf_params%varnce_rt_1, pdf_params%varnce_rt_2,  &
        1412             :       !$acc              pdf_params%thl_1, pdf_params%thl_2, &
        1413             :       !$acc              pdf_params%varnce_thl_1, pdf_params%varnce_thl_2, &
        1414             :       !$acc              pdf_params%corr_w_rt_1, pdf_params%corr_w_rt_2,  &
        1415             :       !$acc              pdf_params%corr_w_thl_1, pdf_params%corr_w_thl_2, &
        1416             :       !$acc              pdf_params%corr_rt_thl_1, pdf_params%corr_rt_thl_2,&
        1417             :       !$acc              pdf_params%alpha_thl, pdf_params%alpha_rt, &
        1418             :       !$acc              pdf_params%crt_1, pdf_params%crt_2, pdf_params%cthl_1, &
        1419             :       !$acc              pdf_params%cthl_2, pdf_params%chi_1, &
        1420             :       !$acc              pdf_params%chi_2, pdf_params%stdev_chi_1, &
        1421             :       !$acc              pdf_params%stdev_chi_2, pdf_params%stdev_eta_1, &
        1422             :       !$acc              pdf_params%stdev_eta_2, pdf_params%covar_chi_eta_1, &
        1423             :       !$acc              pdf_params%covar_chi_eta_2, pdf_params%corr_w_chi_1, &
        1424             :       !$acc              pdf_params%corr_w_chi_2, pdf_params%corr_w_eta_1, &
        1425             :       !$acc              pdf_params%corr_w_eta_2, pdf_params%corr_chi_eta_1, &
        1426             :       !$acc              pdf_params%corr_chi_eta_2, pdf_params%rsatl_1, &
        1427             :       !$acc              pdf_params%rsatl_2, pdf_params%rc_1, pdf_params%rc_2, &
        1428             :       !$acc              pdf_params%cloud_frac_1, pdf_params%cloud_frac_2,  &
        1429             :       !$acc              pdf_params%mixt_frac, pdf_params%ice_supersat_frac_1, &
        1430             :       !$acc              pdf_params%ice_supersat_frac_2 )
        1431             :
        1432             :       !$acc update host( sclrpthvp, sclrprcp, wpsclrp2, wpsclrprtp, wpsclrpthlp, wp2sclrp ) &
        1433             :       !$acc if ( sclr_dim > 0 )
        1434             :
        1435           0 :       do i = 1, ngrdcol
        1436             :
        1437             :         call pdf_closure_check( &
        1438             :                nz, sclr_dim, &
        1439           0 :                wp4(i,:), wprtp2(i,:), wp2rtp(i,:), wpthlp2(i,:), & ! intent(in)
        1440           0 :                wp2thlp(i,:), cloud_frac(i,:), rcm(i,:), wpthvp(i,:), wp2thvp(i,:), &  ! intent(in)
        1441           0 :                rtpthvp(i,:), thlpthvp(i,:), wprcp(i,:), wp2rcp(i,:), & ! intent(in)
        1442           0 :                rtprcp(i,:), thlprcp(i,:), rcp2(i,:), wprtpthlp(i,:), & ! intent(in)
        1443           0 :                pdf_params%crt_1(i,:), pdf_params%crt_2(i,:), & ! intent(in)
        1444           0 :                pdf_params%cthl_1(i,:), pdf_params%cthl_2(i,:), & ! intent(in)
        1445             :                pdf_params, & ! intent(in)
        1446           0 :                sclrpthvp(i,:,:), sclrprcp(i,:,:), wpsclrp2(i,:,:), & ! intent(in)
        1447           0 :                wpsclrprtp(i,:,:), wpsclrpthlp(i,:,:), wp2sclrp(i,:,:), & ! intent(in)
        1448           0 :                stats_metadata ) ! intent(in)
        1449             :       end do
        1450             :     end if
        1451             :
        1452             :     ! Error Reporting
        1453             :     ! Joshua Fasching February 2008
        1454    17870112 :     if ( clubb_at_least_debug_level( 2 ) ) then
        1455           0 :       if ( err_code == clubb_fatal_error ) then
        1456             :
        1457             :         !$acc update host( p_in_Pa, exner, thv_ds, wm, wp2, wp3, sigma_sqd_w, &
        1458             :         !$acc              rtm, rtp2, wprtp, thlm, thlp2, wpthlp, rtpthlp, ice_supersat_frac )
        1459             :
        1460             :         !$acc update host( sclrm, wpsclrp, sclrp2, sclrprtp, sclrpthlp ) if ( sclr_dim > 0 )
        1461             :
        1462           0 :         write(fstderr,*) "Error in pdf_closure_new"
        1463             :
        1464           0 :         write(fstderr,*) "Intent(in)"
        1465             :
        1466           0 :         write(fstderr,*) "p_in_Pa = ", p_in_Pa
        1467           0 :         write(fstderr,*) "exner = ", exner
        1468           0 :         write(fstderr,*) "thv_ds = ", thv_ds
        1469           0 :         write(fstderr,*) "wm = ", wm
        1470           0 :         write(fstderr,*) "wp2 = ", wp2
        1471           0 :         write(fstderr,*) "wp3 = ", wp3
        1472           0 :         write(fstderr,*) "sigma_sqd_w = ", sigma_sqd_w
        1473           0 :         write(fstderr,*) "rtm = ", rtm
        1474           0 :         write(fstderr,*) "rtp2 = ", rtp2
        1475           0 :         write(fstderr,*) "wprtp = ", wprtp
        1476           0 :         write(fstderr,*) "thlm = ", thlm
        1477           0 :         write(fstderr,*) "thlp2 = ", thlp2
        1478           0 :         write(fstderr,*) "wpthlp = ", wpthlp
        1479           0 :         write(fstderr,*) "rtpthlp = ", rtpthlp
        1480             :
        1481           0 :         if ( sclr_dim > 0 ) then
        1482           0 :           write(fstderr,*) "sclrm = ", sclrm
        1483           0 :           write(fstderr,*) "wpsclrp = ", wpsclrp
        1484           0 :           write(fstderr,*) "sclrp2 = ", sclrp2
        1485           0 :           write(fstderr,*) "sclrprtp = ", sclrprtp
        1486           0 :           write(fstderr,*) "sclrpthlp = ", sclrpthlp
        1487             :         end if
        1488             :
        1489           0 :         write(fstderr,*) "Intent(out)"
        1490             :
        1491           0 :         write(fstderr,*) "wp4 = ", wp4
        1492           0 :         if ( l_explicit_turbulent_adv_xpyp .or. stats_metadata%iwprtp2 > 0 ) then
        1493           0 :           write(fstderr,*) "wprtp2 = ", wprtp2
        1494             :         end if
        1495           0 :         write(fstderr,*) "wp2rtp = ", wp2rtp
        1496           0 :         if ( l_explicit_turbulent_adv_xpyp .or. stats_metadata%iwpthlp2 > 0 ) then
        1497           0 :           write(fstderr,*) "wpthlp2 = ", wpthlp2
        1498             :         end if
        1499           0 :         write(fstderr,*) "cloud_frac = ", cloud_frac
        1500           0 :         write(fstderr,*) "ice_supersat_frac = ", ice_supersat_frac
        1501           0 :         write(fstderr,*) "rcm = ", rcm
        1502           0 :         write(fstderr,*) "wpthvp = ", wpthvp
        1503           0 :         write(fstderr,*) "wp2thvp = ", wp2thvp
        1504           0 :         write(fstderr,*) "rtpthvp = ", rtpthvp
        1505           0 :         write(fstderr,*) "thlpthvp = ", thlpthvp
        1506           0 :         write(fstderr,*) "wprcp = ", wprcp
        1507           0 :         write(fstderr,*) "wp2rcp = ", wp2rcp
        1508           0 :         write(fstderr,*) "rtprcp = ", rtprcp
        1509           0 :         write(fstderr,*) "thlprcp = ", thlprcp
        1510             : #ifndef CLUBB_CAM
        1511             :         !  if CLUBB is used in CAM we want this variable computed no matter what
        1512             :         if ( stats_metadata%ircp2 > 0 ) then
        1513             : #endif
        1514           0 :           write(fstderr,*) "rcp2 = ", rcp2
        1515             : #ifndef CLUBB_CAM
        1516             :         end if
        1517             : #endif
        1518           0 :         if ( l_explicit_turbulent_adv_xpyp .or. stats_metadata%iwprtpthlp > 0 ) then
        1519           0 :           write(fstderr,*) "wprtpthlp = ", wprtpthlp
        1520             :         end if
        1521           0 :         write(fstderr,*) "rcp2 = ", rcp2
        1522           0 :         write(fstderr,*) "wprtpthlp = ", wprtpthlp
        1523           0 :         write(fstderr,*) "pdf_params%w_1 = ", pdf_params%w_1
        1524           0 :         write(fstderr,*) "pdf_params%w_2 = ", pdf_params%w_2
        1525           0 :         write(fstderr,*) "pdf_params%varnce_w_1 = ", pdf_params%varnce_w_1
        1526           0 :         write(fstderr,*) "pdf_params%varnce_w_2 = ", pdf_params%varnce_w_2
        1527           0 :         write(fstderr,*) "pdf_params%rt_1 = ", pdf_params%rt_1
        1528           0 :         write(fstderr,*) "pdf_params%rt_2 = ", pdf_params%rt_2
        1529           0 :         write(fstderr,*) "pdf_params%varnce_rt_1 = ", pdf_params%varnce_rt_1
        1530           0 :         write(fstderr,*) "pdf_params%varnce_rt_2 = ", pdf_params%varnce_rt_2
        1531           0 :         write(fstderr,*) "pdf_params%thl_1 = ", pdf_params%thl_1
        1532           0 :         write(fstderr,*) "pdf_params%thl_2 = ", pdf_params%thl_2
        1533           0 :         write(fstderr,*) "pdf_params%varnce_thl_1 = ", pdf_params%varnce_thl_1
        1534           0 :         write(fstderr,*) "pdf_params%varnce_thl_2 = ", pdf_params%varnce_thl_2
        1535           0 :         write(fstderr,*) "pdf_params%corr_w_rt_1 = ", pdf_params%corr_w_rt_1
        1536           0 :         write(fstderr,*) "pdf_params%corr_w_rt_2 = ", pdf_params%corr_w_rt_2
        1537           0 :         write(fstderr,*) "pdf_params%corr_w_thl_1 = ", pdf_params%corr_w_thl_1
        1538           0 :         write(fstderr,*) "pdf_params%corr_w_thl_2 = ", pdf_params%corr_w_thl_2
        1539           0 :         write(fstderr,*) "pdf_params%corr_rt_thl_1 = ", pdf_params%corr_rt_thl_1
        1540           0 :         write(fstderr,*) "pdf_params%corr_rt_thl_2 = ", pdf_params%corr_rt_thl_2
        1541           0 :         write(fstderr,*) "pdf_params%alpha_thl = ", pdf_params%alpha_thl
        1542           0 :         write(fstderr,*) "pdf_params%alpha_rt = ", pdf_params%alpha_rt
        1543           0 :         write(fstderr,*) "pdf_params%crt_1 = ", pdf_params%crt_1
        1544           0 :         write(fstderr,*) "pdf_params%crt_2 = ", pdf_params%crt_2
        1545           0 :         write(fstderr,*) "pdf_params%cthl_1 = ", pdf_params%cthl_1
        1546           0 :         write(fstderr,*) "pdf_params%cthl_2 = ", pdf_params%cthl_2
        1547           0 :         write(fstderr,*) "pdf_params%chi_1 = ", pdf_params%chi_1
        1548           0 :         write(fstderr,*) "pdf_params%chi_2 = ", pdf_params%chi_2
        1549           0 :         write(fstderr,*) "pdf_params%stdev_chi_1 = ", pdf_params%stdev_chi_1
        1550           0 :         write(fstderr,*) "pdf_params%stdev_chi_2 = ", pdf_params%stdev_chi_2
        1551           0 :         write(fstderr,*) "pdf_params%stdev_eta_1 = ", pdf_params%stdev_eta_1
        1552           0 :         write(fstderr,*) "pdf_params%stdev_eta_2 = ", pdf_params%stdev_eta_2
        1553           0 :         write(fstderr,*) "pdf_params%covar_chi_eta_1 = ", &
        1554           0 :                          pdf_params%covar_chi_eta_1
        1555           0 :         write(fstderr,*) "pdf_params%covar_chi_eta_2 = ", &
        1556           0 :                          pdf_params%covar_chi_eta_2
        1557           0 :         write(fstderr,*) "pdf_params%corr_w_chi_1 = ", pdf_params%corr_w_chi_1
        1558           0 :         write(fstderr,*) "pdf_params%corr_w_chi_2 = ", pdf_params%corr_w_chi_2
        1559           0 :         write(fstderr,*) "pdf_params%corr_w_eta_1 = ", pdf_params%corr_w_eta_1
        1560           0 :         write(fstderr,*) "pdf_params%corr_w_eta_2 = ", pdf_params%corr_w_eta_2
        1561           0 :         write(fstderr,*) "pdf_params%corr_chi_eta_1 = ", &
        1562           0 :                          pdf_params%corr_chi_eta_1
        1563           0 :         write(fstderr,*) "pdf_params%corr_chi_eta_2 = ", &
        1564           0 :                          pdf_params%corr_chi_eta_2
        1565           0 :         write(fstderr,*) "pdf_params%rsatl_1 = ", pdf_params%rsatl_1
        1566           0 :         write(fstderr,*) "pdf_params%rsatl_2 = ", pdf_params%rsatl_2
        1567           0 :         write(fstderr,*) "pdf_params%rc_1 = ", pdf_params%rc_1
        1568           0 :         write(fstderr,*) "pdf_params%rc_2 = ", pdf_params%rc_2
        1569           0 :         write(fstderr,*) "pdf_params%cloud_frac_1 = ", pdf_params%cloud_frac_1
        1570           0 :         write(fstderr,*) "pdf_params%cloud_frac_2 = ", pdf_params%cloud_frac_2
        1571           0 :         write(fstderr,*) "pdf_params%mixt_frac = ", pdf_params%mixt_frac
        1572           0 :         write(fstderr,*) "pdf_params%ice_supersat_frac_1 = ", &
        1573           0 :                          pdf_params%ice_supersat_frac_1
        1574           0 :         write(fstderr,*) "pdf_params%ice_supersat_frac_2 = ", &
        1575           0 :                          pdf_params%ice_supersat_frac_2
        1576             :
        1577           0 :         if ( sclr_dim > 0 )then
        1578           0 :           write(fstderr,*) "sclrpthvp = ", sclrpthvp
        1579           0 :           write(fstderr,*) "sclrprcp = ", sclrprcp
        1580           0 :           write(fstderr,*) "wpsclrp2 = ", wpsclrp2
        1581           0 :           write(fstderr,*) "wpsclrprtp = ", wpsclrprtp
        1582           0 :           write(fstderr,*) "wpsclrpthlp = ", wpsclrpthlp
        1583           0 :           write(fstderr,*) "wp2sclrp = ", wp2sclrp
        1584             :         end if
        1585             :
        1586           0 :         return
        1587             :
        1588             :       end if ! Fatal error
        1589             :
        1590           0 :       do i = 1, ngrdcol
        1591             :
        1592             :         ! Error check pdf parameters and moments to ensure consistency
        1593           0 :         if ( iiPDF_type == iiPDF_3D_Luhar ) then
        1594             :
        1595             :           ! Means
        1596           0 :           wm_clubb_pdf(i,:) = pdf_params%mixt_frac(i,:) * pdf_params%w_1(i,:) &
        1597           0 :                          + ( one - pdf_params%mixt_frac(i,:) ) * pdf_params%w_2(i,:)
        1598             :
        1599           0 :           do k = 1, nz, 1
        1600           0 :              if ( abs( ( wm_clubb_pdf(i,k) - wm(i,k) ) &
        1601           0 :                        / max( wm(i,k), eps ) ) > .05_core_rknd ) then
        1602           0 :                 write(fstderr,*) "wm error at thlm = ", thlm(i,k), &
        1603             :                                  ( ( wm_clubb_pdf(i,k) - wm(i,k) ) &
        1604           0 :                                    / max( wm(i,k), eps ) )
        1605             :              end if
        1606             :           end do ! k = 1, nz, 1
        1607             :
        1608           0 :           rtm_clubb_pdf(i,:) = pdf_params%mixt_frac(i,:) * pdf_params%rt_1(i,:) &
        1609           0 :                           + ( one - pdf_params%mixt_frac(i,:) ) * pdf_params%rt_2(i,:)
        1610             :
        1611           0 :           do k = 1, nz, 1
        1612           0 :              if ( abs( ( rtm_clubb_pdf(i,k) - rtm(i,k) ) &
        1613           0 :                        / max( rtm(i,k), eps ) ) > .05_core_rknd ) then
        1614           0 :                 write(fstderr,*) "rtm error at thlm = ", thlm(i,k), &
        1615             :                                  ( ( rtm_clubb_pdf(i,k) - rtm(i,k) ) &
        1616           0 :                                    / max( rtm(i,k), eps ) )
        1617             :              end if
        1618             :           end do ! k = 1, nz, 1
        1619             :
        1620           0 :           thlm_clubb_pdf(i,:) = pdf_params%mixt_frac(i,:) * pdf_params%thl_1(i,:) &
        1621           0 :                            + ( one - pdf_params%mixt_frac(i,:) ) * pdf_params%thl_2(i,:)
        1622             :
        1623           0 :           do k = 1, nz, 1
        1624           0 :              if ( abs( ( thlm_clubb_pdf(i,k) - thlm(i,k) ) / thlm(i,k) ) &
        1625           0 :                   > .05_core_rknd ) then
        1626           0 :                 write(fstderr,*) "thlm error at thlm = ", thlm(i,k), &
        1627           0 :                                  ( ( thlm_clubb_pdf(i,k) - thlm(i,k) ) / thlm(i,k) )
        1628             :              end if
        1629             :           end do ! k = 1, nz, 1
        1630             :
        1631             :           ! Variances
        1632           0 :           wp2_clubb_pdf(i,:) = pdf_params%mixt_frac(i,:) &
        1633           0 :                           * ( ( pdf_params%w_1(i,:) - wm(i,:) )**2 + pdf_params%varnce_w_1(i,:) ) &
        1634             :                           + ( one - pdf_params%mixt_frac(i,:) ) &
        1635           0 :                             * ( ( pdf_params%w_2(i,:) - wm(i,:) )**2 + pdf_params%varnce_w_2(i,:) )
        1636             :
        1637           0 :           do k = 1, nz, 1
        1638           0 :              if ( wp2(i,k) > w_tol**2 ) then
        1639           0 :                 if ( abs( ( wp2_clubb_pdf(i,k) - wp2(i,k) ) / wp2(i,k) ) &
        1640             :                      > .05_core_rknd ) then
        1641           0 :                    write(fstderr,*) "wp2 error at thlm = ", thlm(i,k), &
        1642           0 :                                     ( ( wp2_clubb_pdf(i,k) - wp2(i,k) ) / wp2(i,k) )
        1643             :                 end if
        1644             :              end if
        1645             :           end do ! k = 1, nz, 1
        1646             :
        1647             :           rtp2_clubb_pdf(i,:) &
        1648           0 :           = pdf_params%mixt_frac(i,:) &
        1649           0 :             * ( ( pdf_params%rt_1(i,:) - rtm(i,:) )**2 + pdf_params%varnce_rt_1(i,:) ) &
        1650             :             + ( one - pdf_params%mixt_frac(i,:) ) &
        1651           0 :               * ( ( pdf_params%rt_2(i,:) - rtm(i,:) )**2 + pdf_params%varnce_rt_2(i,:) )
        1652             :
        1653           0 :           do k = 1, nz, 1
        1654           0 :              if ( rtp2(i,k) > rt_tol**2 ) then
        1655           0 :                 if ( abs( ( rtp2_clubb_pdf(i,k) - rtp2(i,k) ) / rtp2(i,k) ) &
        1656             :                      > .05_core_rknd ) then
        1657           0 :                    write(fstderr,*) "rtp2 error at thlm = ", thlm(i,k), &
        1658           0 :                    "Error = ", ( ( rtp2_clubb_pdf(i,k) - rtp2(i,k) ) / rtp2(i,k) )
        1659             :                 end if
        1660             :              end if
        1661             :           end do ! k = 1, nz, 1
        1662             :
        1663             :           thlp2_clubb_pdf(i,:) &
        1664           0 :           = pdf_params%mixt_frac(i,:) &
        1665           0 :             * ( ( pdf_params%thl_1(i,:) - thlm(i,:) )**2 + pdf_params%varnce_thl_1(i,:) ) &
        1666             :             + ( one - pdf_params%mixt_frac(i,:) ) &
        1667           0 :               * ( ( pdf_params%thl_2(i,:) - thlm(i,:) )**2 + pdf_params%varnce_thl_2(i,:) )
        1668             :
        1669           0 :           do k = 1, nz, 1
        1670           0 :              if( thlp2(i,k) > thl_tol**2 ) then
        1671           0 :                 if ( abs( ( thlp2_clubb_pdf(i,k) - thlp2(i,k) ) / thlp2(i,k) ) &
        1672             :                      > .05_core_rknd ) then
        1673           0 :                    write(fstderr,*) "thlp2 error at thlm = ", thlm(i,k), &
        1674           0 :                    "Error = ", ( ( thlp2_clubb_pdf(i,k) - thlp2(i,k) ) / thlp2(i,k) )
        1675             :                 end if
        1676             :              end if
        1677             :           end do ! k = 1, nz, 1
        1678             :
        1679             :           ! Third order moments
        1680             :           wp3_clubb_pdf(i,:) &
        1681           0 :           = pdf_params%mixt_frac(i,:) * ( pdf_params%w_1(i,:) - wm(i,:) ) &
        1682           0 :             * ( ( pdf_params%w_1(i,:) - wm(i,:) )**2 + three * pdf_params%varnce_w_1(i,:) ) &
        1683           0 :             + ( one - pdf_params%mixt_frac(i,:) ) * ( pdf_params%w_2(i,:) - wm(i,:) ) &
        1684           0 :               * ( ( pdf_params%w_2(i,:) - wm(i,:) )**2 + three * pdf_params%varnce_w_2(i,:) )
        1685             :
        1686             :           rtp3_clubb_pdf(i,:) &
        1687           0 :           = pdf_params%mixt_frac(i,:) * ( pdf_params%rt_1(i,:) - rtm(i,:) ) &
        1688           0 :             * ( ( pdf_params%rt_1(i,:) - rtm(i,:) )**2 + three * pdf_params%varnce_rt_1(i,:) ) &
        1689           0 :             + ( one - pdf_params%mixt_frac(i,:) ) * ( pdf_params%rt_2(i,:) - rtm(i,:) ) &
        1690           0 :               * ( ( pdf_params%rt_2(i,:) - rtm(i,:) )**2 + three * pdf_params%varnce_rt_2(i,:) )
        1691             :
        1692             :           thlp3_clubb_pdf(i,:) &
        1693           0 :           = pdf_params%mixt_frac(i,:) * ( pdf_params%thl_1(i,:) - thlm(i,:) ) &
        1694           0 :             * ( ( pdf_params%thl_1(i,:) - thlm(i,:) )**2 + three * pdf_params%varnce_thl_1(i,:) ) &
        1695           0 :             + ( one - pdf_params%mixt_frac(i,:) ) * ( pdf_params%thl_2(i,:) - thlm(i,:) ) &
        1696           0 :               * ( ( pdf_params%thl_2(i,:) - thlm(i,:) )**2 + three * pdf_params%varnce_thl_2(i,:) )
        1697             :
        1698             :           ! Skewness
        1699           0 :           Skw_denom_coef = clubb_params(i,iSkw_denom_coef)
        1700             :
        1701             :           Skw_clubb_pdf(i,:) &
        1702             :           = wp3_clubb_pdf(i,:) &
        1703           0 :             / ( wp2_clubb_pdf(i,:) + Skw_denom_coef * w_tol**2 )**1.5_core_rknd
        1704             :
        1705           0 :           do k = 1, nz, 1
        1706           0 :              if ( Skw(i,k) > .05_core_rknd ) then
        1707           0 :                 if( abs( ( Skw_clubb_pdf(i,k) - Skw(i,k) ) / Skw(i,k) ) &
        1708             :                     > .25_core_rknd ) then
        1709           0 :                    write(fstderr,*) "Skw error at thlm = ", thlm(i,k), &
        1710           0 :                    "Error = ", ( ( Skw_clubb_pdf(i,k) - Skw(i,k) ) / Skw(i,k) ), &
        1711           0 :                    Skw_clubb_pdf(i,k), Skw(i,k)
        1712             :                 end if
        1713             :              end if
        1714             :           end do ! k = 1, nz, 1
        1715             :
        1716             :           Skrt_clubb_pdf(i,:) &
        1717             :           = rtp3_clubb_pdf(i,:) &
        1718           0 :             / ( rtp2_clubb_pdf(i,:) + Skw_denom_coef * rt_tol**2 )**1.5_core_rknd
        1719             :
        1720           0 :           do k = 1, nz, 1
        1721           0 :              if ( Skrt(i,k) > .05_core_rknd ) then
        1722           0 :                 if( abs( ( Skrt_clubb_pdf(i,k) - Skrt(i,k) ) / Skrt(i,k) ) &
        1723             :                     > .25_core_rknd ) then
        1724           0 :                    write(fstderr,*) "Skrt error at thlm = ", thlm(i,k), &
        1725           0 :                    "Error = ", ( ( Skrt_clubb_pdf(i,k) - Skrt(i,k) ) / Skrt(i,k) ), &
        1726           0 :                    Skrt_clubb_pdf(i,k), Skrt(i,k)
        1727             :                 end if
        1728             :              end if
        1729             :           end do ! k = 1, nz, 1
        1730             :
        1731             :           Skthl_clubb_pdf(i,:) &
        1732             :           = thlp3_clubb_pdf(i,:) &
        1733           0 :             / ( thlp2_clubb_pdf(i,:) + Skw_denom_coef * thl_tol**2 )**1.5_core_rknd
        1734             :
        1735           0 :           do k = 1, nz, 1
        1736           0 :              if ( Skthl(i,k) > .05_core_rknd ) then
        1737           0 :                 if ( abs( ( Skthl_clubb_pdf(i,k) - Skthl(i,k) ) / Skthl(i,k) ) &
        1738             :                      > .25_core_rknd ) then
        1739           0 :                    write(fstderr,*) "Skthl error at thlm = ", thlm(i,k), &
        1740           0 :                    "Error = ", ( ( Skthl_clubb_pdf(i,k) - Skthl(i,k) ) / Skthl(i,k) ), &
        1741           0 :                    Skthl_clubb_pdf(i,k), Skthl(i,k)
        1742             :                 end if
        1743             :              end if
        1744             :           end do ! k = 1, nz, 1
        1745             :
        1746             :         end if ! iiPDF_type == iiPDF_3D_Luhar
        1747             :
        1748             :       end do
        1749             :
        1750             :     end if ! clubb_at_least_debug_level
        1751             :
        1752             :     !$acc exit data delete( u_1, u_2, varnce_u_1, varnce_u_2, v_1, v_2, &
        1753             :     !$acc                   varnce_v_1, varnce_v_2, alpha_u, alpha_v, &
        1754             :     !$acc                   corr_u_w_1, corr_u_w_2, corr_v_w_1, corr_v_w_2, &
        1755             :     !$acc                   tl1, tl2, sqrt_wp2, Skthl, &
        1756             :     !$acc                   Skrt, Sku, Skv, wprcp_contrib_comp_1, wprcp_contrib_comp_2, &
        1757             :     !$acc                   wp2rcp_contrib_comp_1, wp2rcp_contrib_comp_2, &
        1758             :     !$acc                   rtprcp_contrib_comp_1, rtprcp_contrib_comp_2, &
        1759             :     !$acc                   thlprcp_contrib_comp_1, thlprcp_contrib_comp_2, &
        1760             :     !$acc                   uprcp_contrib_comp_1, uprcp_contrib_comp_2, &
        1761             :     !$acc                   vprcp_contrib_comp_1, vprcp_contrib_comp_2, &
        1762             :     !$acc                   rc_1_ice, rc_2_ice, rsatl_1, rsatl_2 )
        1763             :
        1764             :     !$acc exit data if( sclr_dim > 0 ) &
        1765             :     !$acc           delete( sclr1, sclr2, varnce_sclr1, varnce_sclr2, &
        1766             :     !$acc                   alpha_sclr, corr_sclr_thl_1, corr_sclr_thl_2, &
        1767             :     !$acc                   corr_sclr_rt_1, corr_sclr_rt_2, corr_w_sclr_1, &
        1768             :     !$acc                   corr_w_sclr_2, Sksclr )
        1769             :
        1770             :     return
        1771             :
        1772             :   end subroutine pdf_closure
        1773             :
        1774             :   !===============================================================================================
        1775    35740224 :   subroutine transform_pdf_chi_eta_component( nz, ngrdcol, &
        1776    35740224 :                                               tl, rsatl, rt, exner,     & ! In
        1777    35740224 :                                               varnce_thl, varnce_rt,    & ! In
        1778    35740224 :                                               corr_rt_thl, chi,         & ! In
        1779    35740224 :                                               crt, cthl,                & ! Out
        1780    35740224 :                                               stdev_chi, stdev_eta,     & ! Out
        1781    35740224 :                                               covar_chi_eta,            & ! Out
        1782    35740224 :                                               corr_chi_eta )              ! Out
        1783             :
        1784             :     use clubb_precision, only: &
        1785             :         core_rknd ! Variable(s)
        1786             :
        1787             :     use constants_clubb, only: &
        1788             :         zero, one, two, &
        1789             :         ep, Lv, Rd, Cp, &
        1790             :         chi_tol, &
        1791             :         eta_tol, &
        1792             :         max_mag_correlation
        1793             :
        1794             :     implicit none
        1795             :
        1796             :     integer, intent(in) :: &
        1797             :       ngrdcol,  & ! Number of grid columns
        1798             :       nz          ! Number of vertical level
        1799             :
        1800             :     ! ----------- Input Variables -----------
        1801             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) ::  &
        1802             :       tl, &
        1803             :       rsatl, &
        1804             :       rt, &
        1805             :       varnce_thl, &
        1806             :       varnce_rt, &
        1807             :       corr_rt_thl, &
        1808             :       exner
        1809             :
        1810             :     ! ----------- Output Variables -----------
        1811             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) ::  &
        1812             :       chi, &            ! s from Lewellen and Yoh 1993 (LY) eqn. 1
        1813             :       crt, &            ! Coefficients for s'
        1814             :       cthl, &           ! Coefficients for s'
        1815             :       stdev_chi, &      ! Standard deviation of chi for each component.
        1816             :       stdev_eta, &      ! Standard deviation of eta for each component.
        1817             :       covar_chi_eta, &  ! Covariance of chi and eta for each component.
        1818             :       corr_chi_eta      ! Correlation of chi and eta for each component.
        1819             :
        1820             :     ! ----------- Local Variables -----------
        1821             :     real( kind = core_rknd ) :: &
        1822             :       varnce_rt_term, &
        1823             :       corr_rt_thl_term, &
        1824             :       varnce_thl_term, &
        1825             :       varnce_chi, &
        1826             :       varnce_eta, &
        1827             :       beta, &
        1828             :       invrs_beta_rsatl_p1
        1829             :
        1830             :     real( kind = core_rknd ), parameter :: &
        1831             :       chi_tol_sqd = chi_tol**2, &
        1832             :       eta_tol_sqd = eta_tol**2, &
        1833             :       Cp_on_Lv = Cp / Lv
        1834             :
        1835             :     ! Loop variable
        1836             :     integer :: k, i
        1837             :
        1838             :     ! ----------- Begin Code -----------
        1839             :
        1840             :     !$acc parallel loop gang vector collapse(2) default(present)
        1841  3073659264 :     do k = 1, nz
        1842 50761923264 :       do i = 1, ngrdcol
        1843             :
        1844             :         ! SD's beta (eqn. 8)
        1845 47688264000 :         beta = ep * Lv**2 / ( Rd * Cp * tl(i,k)**2 )
        1846             :
        1847 47688264000 :         invrs_beta_rsatl_p1 = one / ( one + beta * rsatl(i,k) )
        1848             :
        1849             :         ! s from Lewellen and Yoh 1993 (LY) eqn. 1
        1850 47688264000 :         chi(i,k) = ( rt(i,k) - rsatl(i,k) ) * invrs_beta_rsatl_p1
        1851             :
        1852             :         ! For each normal distribution in the sum of two normal distributions,
        1853             :         ! s' = crt * rt'  +  cthl * thl';
        1854             :         ! therefore, x's' = crt * x'rt'  +  cthl * x'thl'.
        1855             :         ! Larson et al. May, 2001.
        1856 47688264000 :         crt(i,k)  = invrs_beta_rsatl_p1
        1857             :         cthl(i,k) = ( one + beta * rt(i,k) ) * invrs_beta_rsatl_p1**2 &
        1858 50726183040 :                     * Cp_on_Lv * beta * rsatl(i,k) * exner(i,k)
        1859             :
        1860             :       end do
        1861             :     end do
        1862             :     !$acc end parallel loop
        1863             :
        1864             :     ! Calculate covariance, correlation, and standard deviation of
        1865             :     ! chi and eta for each component
        1866             :     ! Include subplume correlation of qt, thl
        1867             :     !$acc parallel loop gang vector collapse(2) default(present)
        1868  3073659264 :     do k = 1, nz
        1869 50761923264 :       do i = 1, ngrdcol
        1870             :
        1871 47688264000 :         varnce_rt_term = crt(i,k)**2 * varnce_rt(i,k)
        1872 47688264000 :         varnce_thl_term = cthl(i,k)**2 * varnce_thl(i,k)
        1873             :
        1874 47688264000 :         covar_chi_eta(i,k) = varnce_rt_term - varnce_thl_term
        1875             :
        1876             :         corr_rt_thl_term = two * corr_rt_thl(i,k) * crt(i,k) * cthl(i,k) &
        1877 47688264000 :                            * sqrt( varnce_rt(i,k) * varnce_thl(i,k) )
        1878             :
        1879 47688264000 :         varnce_chi = varnce_rt_term - corr_rt_thl_term + varnce_thl_term
        1880 47688264000 :         varnce_eta = varnce_rt_term + corr_rt_thl_term + varnce_thl_term
        1881             :
        1882             :         ! We need to introduce a threshold value for the variance of chi and eta
        1883 50726183040 :         if ( varnce_chi < chi_tol_sqd .or. varnce_eta < eta_tol_sqd ) then
        1884             :
        1885  1164135890 :             if ( varnce_chi < chi_tol_sqd ) then
        1886  1163815327 :                 stdev_chi(i,k) = zero  ! Treat chi as a delta function
        1887             :             else
        1888      320563 :                 stdev_chi(i,k) = sqrt( varnce_chi )
        1889             :             end if
        1890             :
        1891  1164135890 :             if ( varnce_eta < eta_tol_sqd ) then
        1892  1163965456 :                 stdev_eta(i,k) = zero  ! Treat eta as a delta function
        1893             :             else
        1894      170434 :                 stdev_eta(i,k) = sqrt( varnce_eta )
        1895             :             end if
        1896             :
        1897  1164135890 :             corr_chi_eta(i,k) = zero
        1898             :
        1899             :         else
        1900             :
        1901 46524128110 :             stdev_chi(i,k) = sqrt( varnce_chi )
        1902 46524128110 :             stdev_eta(i,k) = sqrt( varnce_eta )
        1903             :
        1904 46524128110 :             corr_chi_eta(i,k) = covar_chi_eta(i,k) / ( stdev_chi(i,k) * stdev_eta(i,k) )
        1905             :             corr_chi_eta(i,k) = min( max_mag_correlation, &
        1906 46524128110 :                                    max( -max_mag_correlation, corr_chi_eta(i,k) ) )
        1907             :
        1908             :         end if
        1909             :
        1910             :       end do
        1911             :     end do
        1912             :     !$acc end parallel loop
        1913             :
        1914    35740224 :   end subroutine transform_pdf_chi_eta_component
        1915             :
        1916             :   !=============================================================================
        1917    17870112 :   subroutine calc_wp4_pdf( nz, ngrdcol, &
        1918    17870112 :                            wm, w_1, w_2, &
        1919    17870112 :                            varnce_w_1, varnce_w_2,    &
        1920    17870112 :                            mixt_frac, &
        1921    17870112 :                            wp4 )
        1922             :
        1923             :     ! Description:
        1924             :     ! Calculates <w'^4> by integrating over the PDF of w.  The integral is:
        1925             :     !
        1926             :     ! <w'^4> = INT(-inf:inf) ( w - <w> )^4 P(w) dw;
        1927             :     !
        1928             :     ! where <w> is the overall mean of w and P(w) is a two-component normal
        1929             :     ! distribution of w.  The integrated equation is:
        1930             :     !
        1931             :     ! <w'^4> = mixt_frac * ( 3 * sigma_w_1^4
        1932             :     !                        + 6 * ( mu_w_1 - <w> )^2 * sigma_w_1^2
        1933             :     !                        + ( mu_w_1 - <w> )^4 )
        1934             :     !          + ( 1 - mixt_frac ) * ( 3 * sigma_w_2^4
        1935             :     !                                  + 6 * ( mu_w_2 - <w> )^2 * sigma_w_2^2
        1936             :     !                                  + ( mu_w_2 - <w> )^4 );
        1937             :     !
        1938             :     ! where mu_w_1 is the mean of w in the 1st PDF component, mu_w_2 is the mean
        1939             :     ! of w in the 2nd PDF component, sigma_w_1 is the standard deviation of w in
        1940             :     ! the 1st PDF component, sigma_w_2 is the standard deviation of w in the 2nd
        1941             :     ! PDF component, and mixt_frac is the mixture fraction, which is the weight
        1942             :     ! of the 1st PDF component.
        1943             :
        1944             :     ! References:
        1945             :     !-----------------------------------------------------------------------
        1946             :
        1947             :     use constants_clubb, only: &
        1948             :         six,   & ! Variable(s)
        1949             :         three, &
        1950             :         one
        1951             :
        1952             :     use clubb_precision, only: &
        1953             :         core_rknd    ! Variable(s)
        1954             :
        1955             :     implicit none
        1956             :
        1957             :     integer, intent(in) :: &
        1958             :       ngrdcol,  & ! Number of grid columns
        1959             :       nz          ! Number of vertical level
        1960             :
        1961             :     ! Input Variables
        1962             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) ::  &
        1963             :       wm,         & ! Mean of w (overall)                           [m/s]
        1964             :       w_1,        & ! Mean of w (1st PDF component)                 [m/s]
        1965             :       w_2,        & ! Mean of w (2nd PDF component)                 [m/s]
        1966             :       varnce_w_1, & ! Variance of w (1st PDF component)             [m^2/s^2]
        1967             :       varnce_w_2, & ! Variance of w (2nd PDF component)             [m^2/s^2]
        1968             :       mixt_frac     ! Mixture fraction                              [-]
        1969             :
        1970             :     ! Output Variable
        1971             :     real( kind = core_rknd ), dimension(ngrdcol,nz) ::  &
        1972             :       wp4    ! <w'^4>                   [m^4/s^4]
        1973             :
        1974             :     ! Local Variables
        1975             :     integer :: i, k
        1976             :
        1977             :     !$acc parallel loop gang vector collapse(2) default(present)
        1978  1536829632 :     do k = 1, nz
        1979 25380961632 :       do i = 1, ngrdcol
        1980             :
        1981             :         ! Calculate <w'^4> by integrating over the PDF.
        1982 47688264000 :         wp4(i,k) = mixt_frac(i,k) * ( three * varnce_w_1(i,k)**2 &
        1983             :                             + six * ( ( w_1(i,k) - wm(i,k) )**2 ) * varnce_w_1(i,k) &
        1984             :                             + ( w_1(i,k) - wm(i,k) )**4 ) &
        1985             :                    + ( one - mixt_frac(i,k) ) * ( three * varnce_w_2(i,k)**2 &
        1986             :                                         + six * ( (w_2(i,k) - wm(i,k) )**2 )*varnce_w_2(i,k) &
        1987 73051355520 :                                         + ( w_2(i,k) - wm(i,k) )**4 )
        1988             :       end do
        1989             :     end do
        1990             :     !$acc end parallel loop
        1991             :
        1992    17870112 :     return
        1993             :
        1994             :   end subroutine calc_wp4_pdf
        1995             :
        1996             :   !=============================================================================
        1997    35740224 :   subroutine calc_wp2xp2_pdf( nz, ngrdcol,             &
        1998    35740224 :                               wm, xm, w_1,             &
        1999    35740224 :                               w_2, x_1, x_2,           &
        2000    35740224 :                               varnce_w_1, varnce_w_2,  &
        2001    35740224 :                               varnce_x_1, varnce_x_2,  &
        2002    35740224 :                               corr_w_x_1, corr_w_x_2,  &
        2003    35740224 :                               mixt_frac, &
        2004    35740224 :                               wp2xp2 )
        2005             :
        2006             :     ! Description:
        2007             :     ! Calculates <w'^2x'^2> by integrating over the PDF of w and x.  The
        2008             :     ! integral
        2009             :     ! is:
        2010             :     !
        2011             :     ! <w'^2x'^2>
        2012             :     ! = INT(-inf:inf) INT(-inf:inf) ( w - <w> )^2 ( x - <x> )^2 P(w,x) dx dw;
        2013             :     !
        2014             :     ! where <w> is the overall mean of w, <x> is the overall mean of x, and
        2015             :     ! P(w,x) is a two-component bivariate normal distribution of w and x.  The
        2016             :     ! integrated equation is:
        2017             :     !
        2018             :     ! <w'^2x'^2>
        2019             :     !   = mixt_frac
        2020             :     !      * ( ( mu_w_1 - <w> )**2 * ( ( mu_x_1 - <x> )**2 + sigma_x_1^2 )
        2021             :     !      + four * corr_w_x_1 * sigma_w_1 * sigma_x_1 * ( mu_x_1 - <x> ) * (
        2022             :     !      mu_w_1 - <w> )
        2023             :     !      + ( ( mu_x_1 - <x> )**2 + ( 1 + 2*corr_w_x_1**2 ) * sigma_x_1^2 ) *
        2024             :     !      sigma_w_1^2 )
        2025             :     !    + ( one - mixt_frac )
        2026             :     !      * ( ( mu_w_2 - <w> )**2 * ( ( mu_x_2 - <x> )**2 + sigma_x_2^2 )
        2027             :     !      + four * corr_w_x_2 * sigma_w_2 * sigma_x_2 * ( mu_x_2 - <x> ) * (
        2028             :     !      mu_w_2 - <w> )
        2029             :     !      + ( ( mu_x_2 - <x> )**2 + ( 1 + 2*corr_w_x_2**2 ) * sigma_x_2^2 ) *
        2030             :     !      sigma_w_2^2 )
        2031             :     !
        2032             :     ! where mu_w_1 is the mean of w in the 1st PDF component, mu_w_2 is the mean
        2033             :     ! of w in the 2nd PDF component, mu_x_1 is the mean of x in the 1st PDF
        2034             :     ! component, mu_x_2 is the mean of x in the 2nd PDF component, sigma_w_1 is
        2035             :     ! the standard deviation of w in the 1st PDF component, sigma_w_2 is the
        2036             :     ! standard deviation of w in the 2nd PDF component, sigma_x_1 is the
        2037             :     ! standard deviation of x in the 1st PDF component, sigma_x_2 is the
        2038             :     ! standard deviation of x in the 2nd PDF component, corr_w_x_1 is the
        2039             :     ! correlation of w and x in the 1st PDF component, corr_w_x_2 is the
        2040             :     ! correlation of w and x in the 2nd PDF component, and mixt_frac is the
        2041             :     ! mixture fraction, which is the weight of the 1st PDF component.
        2042             :
        2043             :     ! References:
        2044             :     !-----------------------------------------------------------------------
        2045             :
        2046             :     use constants_clubb, only: &
        2047             :         one,   & ! Variable(s)
        2048             :         four
        2049             :
        2050             :     use clubb_precision, only: &
        2051             :         core_rknd    ! Variable(s)
        2052             :
        2053             :     implicit none
        2054             :
        2055             :     integer, intent(in) :: &
        2056             :       ngrdcol,  & ! Number of grid columns
        2057             :       nz          ! Number of vertical level
        2058             :
        2059             :     ! Input Variables
        2060             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) ::  &
        2061             :       wm,         & ! Mean of w (overall)                       [m/s]
        2062             :       xm,         & ! Mean of x (overall)                       [units vary]
        2063             :       w_1,        & ! Mean of w (1st PDF component)             [m/s]
        2064             :       w_2,        & ! Mean of w (2nd PDF component)             [m/s]
        2065             :       x_1,        & ! Mean of x (1st PDF component)             [units vary]
        2066             :       x_2,        & ! Mean of x (2nd PDF component)             [units vary]
        2067             :       varnce_w_1, & ! Variance of w (1st PDF component)         [m^2/s^2]
        2068             :       varnce_w_2, & ! Variance of w (2nd PDF component)         [m^2/s^2]
        2069             :       varnce_x_1, & ! Variance of x (1st PDF component)         [(units vary)^2]
        2070             :       varnce_x_2, & ! Variance of x (2nd PDF component)         [(units vary)^2]
        2071             :       corr_w_x_1, & ! Correlation of w and x (1st PDF comp.)    [-]
        2072             :       corr_w_x_2, & ! Correlation of w and x (2nd PDF comp.)    [-]
        2073             :       mixt_frac     ! Mixture fraction                          [-]
        2074             :
        2075             :     ! Output Variable
        2076             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) ::  &
        2077             :       wp2xp2        ! <w'^2x'^2>                   [m^2/s^2 (units vary)^2]
        2078             :
        2079             :     ! Local Variable
        2080             :     integer :: i, k
        2081             :
        2082             :     !$acc parallel loop gang vector collapse(2) default(present)
        2083  3073659264 :     do k = 1, nz
        2084 50761923264 :       do i = 1, ngrdcol
        2085             :
        2086             :         ! Calculate <w'x'^2> by integrating over the PDF.
        2087 95376528000 :         wp2xp2(i,k) = mixt_frac(i,k) &
        2088             :                 * ( ( w_1(i,k) - wm(i,k) )**2 * ( ( x_1(i,k) - xm(i,k) )**2 + varnce_x_1(i,k) ) &
        2089             :                 + four * corr_w_x_1(i,k) * sqrt( varnce_w_1(i,k) * varnce_x_1(i,k) ) &
        2090             :                                          * ( x_1(i,k) - xm(i,k) ) * ( w_1(i,k) - wm(i,k) ) &
        2091             :                 + ( ( x_1(i,k) - xm(i,k) )**2 &
        2092             :                 + ( 1 + 2*corr_w_x_1(i,k)**2 ) * varnce_x_1(i,k) ) * varnce_w_1(i,k) ) &
        2093             :                 + ( one - mixt_frac(i,k) ) &
        2094             :                 * ( ( w_2(i,k) - wm(i,k) )**2 * ( ( x_2(i,k) - xm(i,k) )**2 + varnce_x_2(i,k) ) &
        2095             :                 + four * corr_w_x_2(i,k) * sqrt( varnce_w_2(i,k) * varnce_x_2(i,k) ) &
        2096             :                                          * ( x_2(i,k) - xm(i,k) ) * ( w_2(i,k) - wm(i,k) ) &
        2097             :                 + ( ( x_2(i,k) - xm(i,k) )**2 &
        2098 >14610*10^7 :                 + ( 1 + 2*corr_w_x_2(i,k)**2 ) * varnce_x_2(i,k) ) * varnce_w_2(i,k) )
        2099             :       end do
        2100             :     end do
        2101             :     !$acc end parallel loop
        2102             :
        2103    35740224 :     return
        2104             :
        2105             :   end subroutine calc_wp2xp2_pdf
        2106             :
        2107             :   !=============================================================================
        2108    35740224 :   subroutine calc_wp2xp_pdf( nz, ngrdcol,             &
        2109    35740224 :                              wm, xm, w_1, w_2,        &
        2110    35740224 :                              x_1, x_2,                &
        2111    35740224 :                              varnce_w_1, varnce_w_2,  &
        2112    35740224 :                              varnce_x_1, varnce_x_2,  &
        2113    35740224 :                              corr_w_x_1, corr_w_x_2,  &
        2114    35740224 :                              mixt_frac, &
        2115    35740224 :                              wp2xp )
        2116             :
        2117             :     ! Description:
        2118             :     ! Calculates <w'^2 x'> by integrating over the PDF of w and x.  The integral
        2119             :     ! is:
        2120             :     !
        2121             :     ! <w'^2 x'>
        2122             :     ! = INT(-inf:inf) INT(-inf:inf) ( w - <w> )^2 ( x - <x> ) P(w,x) dx dw;
        2123             :     !
        2124             :     ! where <w> is the overall mean of w, <x> is the overall mean of x, and
        2125             :     ! P(w,x) is a two-component bivariate normal distribution of w and x.  The
        2126             :     ! integrated equation is:
        2127             :     !
        2128             :     ! <w'^2 x'>
        2129             :     ! = mixt_frac * ( ( mu_x_1 - <x> ) * ( ( mu_w_1 - <w> )^2 + sigma_w_1^2 )
        2130             :     !                 + 2 * corr_w_x_1 * sigma_w_1 * sigma_x_1
        2131             :     !                   * ( mu_w_1 - <w> ) )
        2132             :     !   + ( 1 - mixt_frac ) * ( ( mu_x_2 - <x> )
        2133             :     !                           * ( ( mu_w_2 - <w> )^2 + sigma_w_2^2 )
        2134             :     !                           + 2 * corr_w_x_2 * sigma_w_2 * sigma_x_2
        2135             :     !                             * ( mu_w_2 - <w> ) );
        2136             :     !
        2137             :     ! where mu_w_1 is the mean of w in the 1st PDF component, mu_w_2 is the mean
        2138             :     ! of w in the 2nd PDF component, mu_x_1 is the mean of x in the 1st PDF
        2139             :     ! component, mu_x_2 is the mean of x in the 2nd PDF component, sigma_w_1 is
        2140             :     ! the standard deviation of w in the 1st PDF component, sigma_w_2 is the
        2141             :     ! standard deviation of w in the 2nd PDF component, sigma_x_1 is the
        2142             :     ! standard deviation of x in the 1st PDF component, sigma_x_2 is the
        2143             :     ! standard deviation of x in the 2nd PDF component, corr_w_x_1 is the
        2144             :     ! correlation of w and x in the 1st PDF component, corr_w_x_2 is the
        2145             :     ! correlation of w and x in the 2nd PDF component, and mixt_frac is the
        2146             :     ! mixture fraction, which is the weight of the 1st PDF component.
        2147             :
        2148             :     ! References:
        2149             :     !-----------------------------------------------------------------------
        2150             :
        2151             :     use constants_clubb, only: &
        2152             :         two,   & ! Variable(s)
        2153             :         one
        2154             :
        2155             :     use clubb_precision, only: &
        2156             :         core_rknd    ! Variable(s)
        2157             :
        2158             :     implicit none
        2159             :
        2160             :     integer, intent(in) :: &
        2161             :       ngrdcol,  & ! Number of grid columns
        2162             :       nz          ! Number of vertical level
        2163             :
        2164             :     ! Input Variables
        2165             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) ::  &
        2166             :       wm,         & ! Mean of w (overall)                       [m/s]
        2167             :       xm,         & ! Mean of x (overall)                       [units vary]
        2168             :       w_1,        & ! Mean of w (1st PDF component)             [m/s]
        2169             :       w_2,        & ! Mean of w (2nd PDF component)             [m/s]
        2170             :       x_1,        & ! Mean of x (1st PDF component)             [units vary]
        2171             :       x_2,        & ! Mean of x (2nd PDF component)             [units vary]
        2172             :       varnce_w_1, & ! Variance of w (1st PDF component)         [m^2/s^2]
        2173             :       varnce_w_2, & ! Variance of w (2nd PDF component)         [m^2/s^2]
        2174             :       varnce_x_1, & ! Variance of x (1st PDF component)         [(units vary)^2]
        2175             :       varnce_x_2, & ! Variance of x (2nd PDF component)         [(units vary)^2]
        2176             :       corr_w_x_1, & ! Correlation of w and x (1st PDF comp.)    [-]
        2177             :       corr_w_x_2, & ! Correlation of w and x (2nd PDF comp.)    [-]
        2178             :       mixt_frac     ! Mixture fraction                          [-]
        2179             :
        2180             :     ! Output Variable
        2181             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) ::  &
        2182             :       wp2xp    ! <w'^2 x'>                   [m^2/s^2 (units vary)]
        2183             :
        2184             :     ! Local Variables
        2185             :     integer :: i, k
        2186             :
        2187             :
        2188             :     ! Calculate <w'^2 x'> by integrating over the PDF.
        2189             :     !$acc parallel loop gang vector collapse(2) default(present)
        2190  3073659264 :     do k = 1, nz
        2191 50761923264 :       do i = 1, ngrdcol
        2192             :
        2193 95376528000 :         wp2xp(i,k)  = mixt_frac(i,k) &
        2194             :                    * ( ( ( w_1(i,k) - wm(i,k) )**2 + varnce_w_1(i,k) ) * ( x_1(i,k) - xm(i,k) ) &
        2195             :                        + two * corr_w_x_1(i,k) * sqrt( varnce_w_1(i,k) * varnce_x_1(i,k) ) &
        2196             :                          * ( w_1(i,k) - wm(i,k) ) ) &
        2197             :                    + ( one - mixt_frac(i,k) ) &
        2198             :                      * ( ( ( w_2(i,k) - wm(i,k) )**2 + varnce_w_2(i,k) ) * ( x_2(i,k) - xm(i,k) ) &
        2199             :                          + two * corr_w_x_2(i,k) * sqrt( varnce_w_2(i,k) * varnce_x_2(i,k) ) &
        2200 >14610*10^7 :                            * ( w_2(i,k) - wm(i,k) ) )
        2201             :       end do
        2202             :     end do
        2203             :     !$acc end parallel loop
        2204             :
        2205    35740224 :     return
        2206             :
        2207             :   end subroutine calc_wp2xp_pdf
        2208             :
        2209             :   !=============================================================================
        2210    35740224 :   subroutine calc_wpxp2_pdf( nz, ngrdcol,             &
        2211    35740224 :                              wm, xm, w_1,             &
        2212    35740224 :                              w_2, x_1, x_2,           &
        2213    35740224 :                              varnce_w_1, varnce_w_2,  &
        2214    35740224 :                              varnce_x_1, varnce_x_2,  &
        2215    35740224 :                              corr_w_x_1, corr_w_x_2,  &
        2216    35740224 :                              mixt_frac, &
        2217    35740224 :                              wpxp2 )
        2218             :
        2219             :     ! Description:
        2220             :     ! Calculates <w'x'^2> by integrating over the PDF of w and x.  The integral
        2221             :     ! is:
        2222             :     !
        2223             :     ! <w'x'^2>
        2224             :     ! = INT(-inf:inf) INT(-inf:inf) ( w - <w> ) ( x - <x> )^2 P(w,x) dx dw;
        2225             :     !
        2226             :     ! where <w> is the overall mean of w, <x> is the overall mean of x, and
        2227             :     ! P(w,x) is a two-component bivariate normal distribution of w and x.  The
        2228             :     ! integrated equation is:
        2229             :     !
        2230             :     ! <w'x'^2>
        2231             :     ! = mixt_frac * ( ( mu_w_1 - <w> ) * ( ( mu_x_1 - <x> )^2 + sigma_x_1^2 )
        2232             :     !                 + 2 * corr_w_x_1 * sigma_w_1 * sigma_x_1
        2233             :     !                   * ( mu_x_1 - <x> ) )
        2234             :     !   + ( 1 - mixt_frac ) * ( ( mu_w_2 - <w> )
        2235             :     !                           * ( ( mu_x_2 - <x> )^2 + sigma_x_2^2 )
        2236             :     !                           + 2 * corr_w_x_2 * sigma_w_2 * sigma_x_2
        2237             :     !                             * ( mu_x_2 - <x> ) );
        2238             :     !
        2239             :     ! where mu_w_1 is the mean of w in the 1st PDF component, mu_w_2 is the mean
        2240             :     ! of w in the 2nd PDF component, mu_x_1 is the mean of x in the 1st PDF
        2241             :     ! component, mu_x_2 is the mean of x in the 2nd PDF component, sigma_w_1 is
        2242             :     ! the standard deviation of w in the 1st PDF component, sigma_w_2 is the
        2243             :     ! standard deviation of w in the 2nd PDF component, sigma_x_1 is the
        2244             :     ! standard deviation of x in the 1st PDF component, sigma_x_2 is the
        2245             :     ! standard deviation of x in the 2nd PDF component, corr_w_x_1 is the
        2246             :     ! correlation of w and x in the 1st PDF component, corr_w_x_2 is the
        2247             :     ! correlation of w and x in the 2nd PDF component, and mixt_frac is the
        2248             :     ! mixture fraction, which is the weight of the 1st PDF component.
        2249             :
        2250             :     ! References:
        2251             :     !-----------------------------------------------------------------------
        2252             :
        2253             :     use constants_clubb, only: &
        2254             :         two,   & ! Variable(s)
        2255             :         one
        2256             :
        2257             :     use clubb_precision, only: &
        2258             :         core_rknd    ! Variable(s)
        2259             :
        2260             :     implicit none
        2261             :
        2262             :     integer, intent(in) :: &
        2263             :       ngrdcol,  & ! Number of grid columns
        2264             :       nz          ! Number of vertical level
        2265             :
        2266             :     ! Input Variables
        2267             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) ::  &
        2268             :       wm,         & ! Mean of w (overall)                       [m/s]
        2269             :       xm,         & ! Mean of x (overall)                       [units vary]
        2270             :       w_1,        & ! Mean of w (1st PDF component)             [m/s]
        2271             :       w_2,        & ! Mean of w (2nd PDF component)             [m/s]
        2272             :       x_1,        & ! Mean of x (1st PDF component)             [units vary]
        2273             :       x_2,        & ! Mean of x (2nd PDF component)             [units vary]
        2274             :       varnce_w_1, & ! Variance of w (1st PDF component)         [m^2/s^2]
        2275             :       varnce_w_2, & ! Variance of w (2nd PDF component)         [m^2/s^2]
        2276             :       varnce_x_1, & ! Variance of x (1st PDF component)         [(units vary)^2]
        2277             :       varnce_x_2, & ! Variance of x (2nd PDF component)         [(units vary)^2]
        2278             :       corr_w_x_1, & ! Correlation of w and x (1st PDF comp.)    [-]
        2279             :       corr_w_x_2, & ! Correlation of w and x (2nd PDF comp.)    [-]
        2280             :       mixt_frac     ! Mixture fraction                          [-]
        2281             :
        2282             :     ! Return Variable
        2283             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) ::  &
        2284             :       wpxp2    ! <w'x'^2>                   [m/s (units vary)^2]
        2285             :
        2286             :     ! Local Variables
        2287             :     integer :: i, k
        2288             :
        2289             :     !$acc parallel loop gang vector collapse(2) default(present)
        2290  3073659264 :     do k = 1, nz
        2291 50761923264 :       do i = 1, ngrdcol
        2292             :
        2293             :         ! Calculate <w'x'^2> by integrating over the PDF.
        2294 95376528000 :         wpxp2(i,k) = mixt_frac(i,k) &
        2295             :                 * ( ( w_1(i,k) - wm(i,k) ) * ( ( x_1(i,k) - xm(i,k) )**2 + varnce_x_1(i,k) ) &
        2296             :                     + two * corr_w_x_1(i,k) * sqrt( varnce_w_1(i,k) * varnce_x_1(i,k) ) &
        2297             :                       * ( x_1(i,k) - xm(i,k) ) ) &
        2298             :                 + ( one - mixt_frac(i,k) ) &
        2299             :                   * ( ( w_2(i,k) - wm(i,k) ) * ( ( x_2(i,k) - xm(i,k) )**2 + varnce_x_2(i,k) ) &
        2300             :                       + two * corr_w_x_2(i,k) * sqrt( varnce_w_2(i,k) * varnce_x_2(i,k) ) &
        2301 >14610*10^7 :                         * ( x_2(i,k) - xm(i,k) ) )
        2302             :       end do
        2303             :     end do
        2304             :     !$acc end parallel loop
        2305             :
        2306    35740224 :     return
        2307             :
        2308             :   end subroutine calc_wpxp2_pdf
        2309             :
        2310             :   !=============================================================================
        2311           0 :   subroutine calc_wpxpyp_pdf( nz, ngrdcol, &
        2312           0 :                               wm, xm, ym, w_1, w_2,   &
        2313           0 :                               x_1, x_2,               &
        2314           0 :                               y_1, y_2,               &
        2315           0 :                               varnce_w_1, varnce_w_2, &
        2316           0 :                               varnce_x_1, varnce_x_2, &
        2317           0 :                               varnce_y_1, varnce_y_2, &
        2318           0 :                               corr_w_x_1, corr_w_x_2, &
        2319           0 :                               corr_w_y_1, corr_w_y_2, &
        2320           0 :                               corr_x_y_1, corr_x_y_2, &
        2321           0 :                               mixt_frac, &
        2322           0 :                               wpxpyp )
        2323             :
        2324             :     ! Description:
        2325             :     ! Calculates <w'x'y'> by integrating over the PDF of w, x, and y.  The
        2326             :     ! integral is:
        2327             :     !
        2328             :     ! <w'x'y'>
        2329             :     ! = INT(-inf:inf) INT(-inf:inf) INT(-inf:inf)
        2330             :     !   ( w - <w> ) ( x - <x> ) ( y - <y> ) P(w,x,y) dy dx dw;
        2331             :     !
        2332             :     ! where <w> is the overall mean of w, <x> is the overall mean of x, <y> is
        2333             :     ! the overall mean of y, and P(w,x,y) is a two-component trivariate normal
        2334             :     ! distribution of w, x, and y.  The integrated equation is:
        2335             :     !
        2336             :     ! <w'x'y'>
        2337             :     ! = mixt_frac
        2338             :     !   * ( ( mu_w_1 - <w> ) * ( mu_x_1 - <x> ) * ( mu_y_1 - <y> )
        2339             :     !       + corr_x_y_1 * sigma_x_1 * sigma_y_1 * ( mu_w_1 - <w> )
        2340             :     !       + corr_w_y_1 * sigma_w_1 * sigma_y_1 * ( mu_x_1 - <x> )
        2341             :     !       + corr_w_x_1 * sigma_w_1 * sigma_x_1 * ( mu_y_1 - <y> ) )
        2342             :     !   + ( 1 - mixt_frac )
        2343             :     !     * ( ( mu_w_2 - <w> ) * ( mu_x_2 - <x> ) * ( mu_y_2 - <y> )
        2344             :     !         + corr_x_y_2 * sigma_x_2 * sigma_y_2 * ( mu_w_2 - <w> )
        2345             :     !         + corr_w_y_2 * sigma_w_2 * sigma_y_2 * ( mu_x_2 - <x> )
        2346             :     !         + corr_w_x_2 * sigma_w_2 * sigma_x_2 * ( mu_y_2 - <y> ) );
        2347             :     !
        2348             :     ! where mu_w_1 is the mean of w in the 1st PDF component, mu_w_2 is the mean
        2349             :     ! of w in the 2nd PDF component, mu_x_1 is the mean of x in the 1st PDF
        2350             :     ! component, mu_x_2 is the mean of x in the 2nd PDF component, mu_y_1 is the
        2351             :     ! mean of y in the 1st PDF component, mu_y_2 is the mean of y in the 2nd PDF
        2352             :     ! component, sigma_w_1 is the standard deviation of w in the 1st PDF
        2353             :     ! component, sigma_w_2 is the standard deviation of w in the 2nd PDF
        2354             :     ! component, sigma_x_1 is the standard deviation of x in the 1st PDF
        2355             :     ! component, sigma_x_2 is the standard deviation of x in the 2nd PDF
        2356             :     ! component, sigma_y_1 is the standard deviation of y in the 1st PDF
        2357             :     ! component, sigma_y_2 is the standard deviation of y in the 2nd PDF
        2358             :     ! component, corr_w_x_1 is the correlation of w and x in the 1st PDF
        2359             :     ! component, corr_w_x_2 is the correlation of w and x in the 2nd PDF
        2360             :     ! component, corr_w_y_1 is the correlation of w and y in the 1st PDF
        2361             :     ! component, corr_w_y_2 is the correlation of w and y in the 2nd PDF
        2362             :     ! component, corr_x_y_1 is the correlation of x and y in the 1st PDF
        2363             :     ! component, corr_x_y_2 is the correlation of x and y in the 2nd PDF
        2364             :     ! component, and mixt_frac is the mixture fraction, which is the weight of
        2365             :     ! the 1st PDF component.
        2366             :
        2367             :     ! References:
        2368             :     !-----------------------------------------------------------------------
        2369             :
        2370             :     use grid_class, only: &
        2371             :         grid ! Type
        2372             :
        2373             :     use constants_clubb, only: &
        2374             :         one    ! Variable(s)
        2375             :
        2376             :     use clubb_precision, only: &
        2377             :         core_rknd    ! Variable(s)
        2378             :
        2379             :     implicit none
        2380             :
        2381             :     integer, intent(in) :: &
        2382             :       ngrdcol,  & ! Number of grid columns
        2383             :       nz          ! Number of vertical level
        2384             :
        2385             :     ! Input Variables
        2386             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) ::  &
        2387             :       wm,         & ! Mean of w (overall)                          [m/s]
        2388             :       xm,         & ! Mean of x (overall)                          [x units]
        2389             :       ym,         & ! Mean of y (overall)                          [y units]
        2390             :       w_1,        & ! Mean of w (1st PDF component)                [m/s]
        2391             :       w_2,        & ! Mean of w (2nd PDF component)                [m/s]
        2392             :       x_1,        & ! Mean of x (1st PDF component)                [x units]
        2393             :       x_2,        & ! Mean of x (2nd PDF component)                [x units]
        2394             :       y_1,        & ! Mean of y (1st PDF component)                [y units]
        2395             :       y_2,        & ! Mean of y (2nd PDF component)                [y units]
        2396             :       varnce_w_1, & ! Variance of w (1st PDF component)            [m^2/s^2]
        2397             :       varnce_w_2, & ! Variance of w (2nd PDF component)            [m^2/s^2]
        2398             :       varnce_x_1, & ! Variance of x (1st PDF component)            [(x units)^2]
        2399             :       varnce_x_2, & ! Variance of x (2nd PDF component)            [(x units)^2]
        2400             :       varnce_y_1, & ! Variance of y (1st PDF component)            [(y units)^2]
        2401             :       varnce_y_2, & ! Variance of y (2nd PDF component)            [(y units)^2]
        2402             :       corr_w_x_1, & ! Correlation of w and x (1st PDF component)   [-]
        2403             :       corr_w_x_2, & ! Correlation of w and x (2nd PDF component)   [-]
        2404             :       corr_w_y_1, & ! Correlation of w and y (1st PDF component)   [-]
        2405             :       corr_w_y_2, & ! Correlation of w and y (2nd PDF component)   [-]
        2406             :       corr_x_y_1, & ! Correlation of x and y (1st PDF component)   [-]
        2407             :       corr_x_y_2, & ! Correlation of x and y (2nd PDF component)   [-]
        2408             :       mixt_frac     ! Mixture fraction                             [-]
        2409             :
        2410             :     ! Output Variable
        2411             :     real( kind = core_rknd ), dimension(ngrdcol,nz) ::  &
        2412             :       wpxpyp    ! <w'x'y'>                   [m/s (units vary)]
        2413             :
        2414             :     ! Local Variables
        2415             :     integer :: i, k
        2416             :
        2417             :
        2418             :     ! Calculate <w'x'y'> by integrating over the PDF.
        2419             :     !$acc parallel loop gang vector collapse(2) default(present)
        2420           0 :     do k = 1, nz
        2421           0 :       do i = 1, ngrdcol
        2422           0 :         wpxpyp(i,k) &
        2423             :         = mixt_frac(i,k) &
        2424             :           * ( ( w_1(i,k) - wm(i,k) ) * ( x_1(i,k) - xm(i,k) ) * ( y_1(i,k) - ym(i,k) ) &
        2425             :               + corr_x_y_1(i,k)*sqrt( varnce_x_1(i,k)*varnce_y_1(i,k) )*( w_1(i,k)-wm(i,k) ) &
        2426             :               + corr_w_y_1(i,k)*sqrt( varnce_w_1(i,k)*varnce_y_1(i,k) )*( x_1(i,k)-xm(i,k) ) &
        2427             :               + corr_w_x_1(i,k)*sqrt( varnce_w_1(i,k)*varnce_x_1(i,k) )*( y_1(i,k)-ym(i,k) ) ) &
        2428             :           + ( one - mixt_frac(i,k) ) &
        2429             :             * ( ( w_2(i,k) - wm(i,k) )*( x_2(i,k) - xm(i,k) ) * ( y_2(i,k) - ym(i,k) ) &
        2430             :                 + corr_x_y_2(i,k)*sqrt( varnce_x_2(i,k)*varnce_y_2(i,k) )*( w_2(i,k)-wm(i,k) ) &
        2431             :                 + corr_w_y_2(i,k)*sqrt( varnce_w_2(i,k)*varnce_y_2(i,k) )*( x_2(i,k)-xm(i,k) ) &
        2432           0 :                 + corr_w_x_2(i,k)*sqrt( varnce_w_2(i,k)*varnce_x_2(i,k) )*( y_2(i,k)-ym(i,k) ) )
        2433             :       end do
        2434             :     end do
        2435             :     !$acc end parallel loop
        2436             :
        2437           0 :     return
        2438             :
        2439             :   end subroutine calc_wpxpyp_pdf
        2440             :
        2441             :   !=============================================================================
        2442    35740224 :   subroutine calc_liquid_cloud_frac_component( nz, ngrdcol, &
        2443    35740224 :                                                mean_chi, stdev_chi, &
        2444    35740224 :                                                cloud_frac, rc )
        2445             :     ! Description:
        2446             :     ! Calculates the PDF component cloud water mixing ratio, rc_i, and cloud
        2447             :     ! fraction, cloud_frac_i, for the ith PDF component.
        2448             :     !
        2449             :     ! The equation for cloud water mixing ratio, rc, at any point is:
        2450             :     !
        2451             :     ! rc = chi * H(chi);
        2452             :     !
        2453             :     ! and the equation for cloud fraction at a point, fc, is:
        2454             :     !
        2455             :     ! fc = H(chi);
        2456             :     !
        2457             :     ! where where extended liquid water mixing ratio, chi, is equal to cloud
        2458             :     ! water mixing ratio, rc, when positive.  When the atmosphere is saturated
        2459             :     ! at this point, cloud water is found, and rc = chi, while fc = 1.
        2460             :     ! Otherwise, clear air is found at this point, and rc = fc = 0.
        2461             :     !
        2462             :     ! The mean of rc and fc is calculated by integrating over the PDF, such
        2463             :     ! that:
        2464             :     !
        2465             :     ! <rc> = INT(-inf:inf) chi * H(chi) * P(chi) dchi; and
        2466             :     !
        2467             :     ! cloud_frac = <fc> = INT(-inf:inf) H(chi) * P(chi) dchi.
        2468             :     !
        2469             :     ! This can be rewritten as:
        2470             :     !
        2471             :     ! <rc> = INT(0:inf) chi * P(chi) dchi; and
        2472             :     !
        2473             :     ! cloud_frac = <fc> = INT(0:inf) P(chi) dchi;
        2474             :     !
        2475             :     ! and further rewritten as:
        2476             :     !
        2477             :     ! <rc> = SUM(i=1,N) mixt_frac_i INT(0:inf) chi * P_i(chi) dchi; and
        2478             :     !
        2479             :     ! cloud_frac = SUM(i=1,N) mixt_frac_i INT(0:inf) P_i(chi) dchi;
        2480             :     !
        2481             :     ! where N is the number of PDF components.  The equation for mean rc in the
        2482             :     ! ith PDF component is:
        2483             :     !
        2484             :     ! rc_i = INT(0:inf) chi * P_i(chi) dchi;
        2485             :     !
        2486             :     ! and the equation for cloud fraction in the ith PDF component is:
        2487             :     !
        2488             :     ! cloud_frac_i = INT(0:inf) P_i(chi) dchi.
        2489             :     !
        2490             :     ! The component values are related to the overall values by:
        2491             :     !
        2492             :     ! <rc> = SUM(i=1,N) mixt_frac_i * rc_i; and
        2493             :     !
        2494             :     ! cloud_frac = SUM(i=1,N) mixt_frac_i * cloud_frac_i.
        2495             :
        2496             :     ! References:
        2497             :     !----------------------------------------------------------------------
        2498             :
        2499             :     use constants_clubb, only: &
        2500             :         chi_tol,        & ! Tolerance for pdf parameter chi       [kg/kg]
        2501             :         sqrt_2pi,       & ! sqrt(2*pi)
        2502             :         sqrt_2,         & ! sqrt(2)
        2503             :         one,            & ! 1
        2504             :         one_half,       & ! 1/2
        2505             :         zero,           & ! 0
        2506             :         max_num_stdevs, &
        2507             :         eps
        2508             :
        2509             :     use clubb_precision, only: &
        2510             :         core_rknd     ! Precision
        2511             :
        2512             :     implicit none
        2513             :
        2514             :     integer, intent(in) :: &
        2515             :       ngrdcol,  & ! Number of grid columns
        2516             :       nz          ! Number of vertical level
        2517             :
        2518             :     !----------- Input Variables -----------
        2519             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
        2520             :       mean_chi,  & ! Mean of chi (old s) (ith PDF component)           [kg/kg]
        2521             :       stdev_chi    ! Standard deviation of chi (ith PDF component)     [kg/kg]
        2522             :
        2523             :     !----------- Output Variables -----------
        2524             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
        2525             :       cloud_frac, & ! Cloud fraction (ith PDF component)               [-]
        2526             :       rc            ! Mean cloud water mixing ratio (ith PDF comp.)    [kg/kg]
        2527             :
        2528             :     !----------- Local Variables -----------
        2529             :     real( kind = core_rknd), parameter :: &
        2530             :       invrs_sqrt_2 = one / sqrt_2, &
        2531             :       invrs_sqrt_2pi = one / sqrt_2pi
        2532             :
        2533             :     real( kind = core_rknd ) :: &
        2534             :       zeta
        2535             :
        2536             :     integer :: k, i    ! Vertical loop index
        2537             :
        2538             :     !----------- Begin Code -----------
        2539             :     !$acc parallel loop gang vector collapse(2) default(present)
        2540  3073659264 :     do k = 1, nz
        2541 50761923264 :       do i = 1, ngrdcol
        2542             :
        2543 95376528000 :         if ( ( abs( mean_chi(i,k) ) <= eps .and. stdev_chi(i,k) <= chi_tol ) &
        2544 >14610*10^7 :                .or. ( mean_chi(i,k) < - max_num_stdevs * stdev_chi(i,k) ) ) then
        2545             :
        2546             :             ! The mean of chi is at saturation and does not vary in the ith PDF component
        2547 44856697767 :             cloud_frac(i,k) = zero
        2548 44856697767 :             rc(i,k)         = zero
        2549             :
        2550  2831566233 :         elseif ( mean_chi(i,k) > max_num_stdevs * stdev_chi(i,k) ) then
        2551             :
        2552             :             ! The mean of chi is multiple standard deviations above the saturation point.
        2553             :             ! Thus, all cloud in the ith PDF component.
        2554   297535127 :             cloud_frac(i,k) = one
        2555   297535127 :             rc(i,k)         = mean_chi(i,k)
        2556             :
        2557             :         else
        2558             :
        2559             :             ! The mean of chi is within max_num_stdevs of the saturation point.
        2560             :             ! Thus, layer is partly cloudy, requires calculation.
        2561             :
        2562  2534031106 :             zeta = mean_chi(i,k) / stdev_chi(i,k)
        2563             :
        2564  2534031106 :             cloud_frac(i,k) = one_half * ( one + erf( zeta * invrs_sqrt_2 )  )
        2565             :
        2566             :             rc(i,k) = mean_chi(i,k) * cloud_frac(i,k) &
        2567  2534031106 :                       + stdev_chi(i,k) * exp( - one_half * zeta**2 ) * invrs_sqrt_2pi
        2568             :
        2569             :         end if
        2570             :
        2571             :       end do
        2572             :     end do
        2573             :     !$acc end parallel loop
        2574             :
        2575    35740224 :     return
        2576             :
        2577             :   end subroutine calc_liquid_cloud_frac_component
        2578             :
        2579             :   !=============================================================================
        2580    35740224 :   subroutine calc_ice_cloud_frac_component( nz, ngrdcol, &
        2581    35740224 :                                             mean_chi, stdev_chi, &
        2582    35740224 :                                             rc_in, cloud_frac, &
        2583    35740224 :                                             p_in_Pa, tl, &
        2584    35740224 :                                             rsatl, crt, &
        2585             :                                             saturation_formula, &
        2586    35740224 :                                             ice_supersat_frac, rc )
        2587             :   ! Description:
        2588             :   !   A version of the cloud fraction calculation modified to work
        2589             :   !   for layers that are potentially below freezing. If there are
        2590             :   !   no below freezing levels, the ice_supersat_frac calculation is
        2591             :   !   the same as cloud_frac.
        2592             :   !
        2593             :   !   For the below freezing levels, the saturation point will be
        2594             :   !   non-zero, thus we need to calculate chi_at_ice_sat.
        2595             :   !
        2596             :   !   The description of the equations are located in the description
        2597             :   !   of calc_liquid_cloud_frac_component.
        2598             :   !----------------------------------------------------------------------
        2599             :
        2600             :     use constants_clubb, only: &
        2601             :         chi_tol,        & ! Tolerance for pdf parameter chi       [kg/kg]
        2602             :         T_freeze_K,     & ! Freezing point of water             [K]
        2603             :         sqrt_2pi,       & ! sqrt(2*pi)
        2604             :         sqrt_2,         & ! sqrt(2)
        2605             :         one,            & ! 1
        2606             :         one_half,       & ! 1/2
        2607             :         zero,           & ! 0
        2608             :         max_num_stdevs, &
        2609             :         eps
        2610             :
        2611             :     use clubb_precision, only: &
        2612             :         core_rknd     ! Precision
        2613             :
        2614             :     use saturation, only:  &
        2615             :         sat_mixrat_ice
        2616             :
        2617             :     implicit none
        2618             :
        2619             :     ! ---------------------- Input Variables ----------------------
        2620             :     integer, intent(in) :: &
        2621             :       ngrdcol,  & ! Number of grid columns
        2622             :       nz          ! Number of vertical level
        2623             :
        2624             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
        2625             :       mean_chi,   & ! Mean of chi (old s) (ith PDF component)           [kg/kg]
        2626             :       stdev_chi,  & ! Standard deviation of chi (ith PDF component)     [kg/kg]
        2627             :       rc_in,      & ! Mean cloud water mixing ratio (ith PDF comp.)     [kg/kg]
        2628             :       cloud_frac, & ! Cloud fraction                                    [-]
        2629             :       p_in_Pa,    & ! Pressure                                          [Pa]
        2630             :       rsatl,      & ! Saturation mixing ratio of liquid                 [kg/kg]
        2631             :       crt,        & ! r_t coef. in chi/eta eqns.                        [-]
        2632             :       tl            ! Quantities needed to predict higher order moments
        2633             :                     ! tl = thl*exner
        2634             :
        2635             :     integer, intent(in) :: &
        2636             :       saturation_formula ! Integer that stores the saturation formula to be used
        2637             :
        2638             :     ! ---------------------- Output Variables ----------------------
        2639             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
        2640             :       ice_supersat_frac,  & ! Ice supersaturation fraction                [-]
        2641             :       rc                    ! Mean cloud ice mixing ratio (ith PDF comp.) [kg/kg]
        2642             :
        2643             :     ! ---------------------- Local Variables----------------------
        2644             :     real( kind = core_rknd), parameter :: &
        2645             :       invrs_sqrt_2 = one / sqrt_2, &
        2646             :       invrs_sqrt_2pi = one / sqrt_2pi
        2647             :
        2648             :     real( kind = core_rknd ) :: &
        2649             :       zeta, &
        2650             :       chi_at_ice_sat
        2651             :
        2652             :     real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
        2653    71480448 :       rsat_ice
        2654             :
        2655             :     integer :: k, i    ! Loop indices
        2656             :
        2657             :     logical :: &
        2658             :       l_any_below_freezing
        2659             :
        2660             :     ! ---------------------- Begin Code ----------------------
        2661             :
        2662    35740224 :     l_any_below_freezing = .false.
        2663             :
        2664             :     ! If a grid boxes is above freezing, then the calculation is the
        2665             :     ! same as the cloud_frac calculation
        2666             :     !$acc parallel loop gang vector collapse(2) default(present) &
        2667             :     !$acc          reduction(.or.:l_any_below_freezing)
        2668  3073659264 :     do k = 1, nz
        2669 50761923264 :       do i = 1, ngrdcol
        2670 50726183040 :         if ( tl(i,k) > T_freeze_K ) then
        2671  9318427162 :           ice_supersat_frac(i,k) = cloud_frac(i,k)
        2672  9318427162 :           rc(i,k)                = rc_in(i,k)
        2673             :         else
        2674             :           l_any_below_freezing = .true.
        2675             :         end if
        2676             :       end do
        2677             :     end do
        2678             :     !$acc end parallel loop
        2679             :
        2680             :     ! If all grid boxes are above freezing, then the calculation is the
        2681             :     ! same as the cloud_frac calculation
        2682    35740224 :     if ( .not. l_any_below_freezing ) then
        2683             :       return
        2684             :     end if
        2685             :
        2686             :     !$acc data create( rsat_ice )
        2687             :
        2688             :     ! Calculate the saturation mixing ratio of ice
        2689    35740224 :     rsat_ice = sat_mixrat_ice( nz, ngrdcol, p_in_Pa, tl, saturation_formula )
        2690             :
        2691             :     !$acc parallel loop gang vector collapse(2) default(present)
        2692  3073659264 :     do k = 1, nz
        2693 50761923264 :       do i = 1, ngrdcol
        2694             :
        2695 50726183040 :         if ( tl(i,k) <= T_freeze_K ) then
        2696             :
        2697             :           ! Temperature is freezing, we must compute chi_at_ice_sat and
        2698             :           ! calculate the new cloud_frac_component
        2699 38369836838 :           chi_at_ice_sat = crt(i,k) * ( rsat_ice(i,k) - rsatl(i,k) )
        2700             :
        2701             :           if ( ( abs( mean_chi(i,k)-chi_at_ice_sat ) <= eps .and. stdev_chi(i,k) <= chi_tol ) &
        2702 38369836838 :            .or. ( mean_chi(i,k)-chi_at_ice_sat < - max_num_stdevs * stdev_chi(i,k) ) ) then
        2703             :
        2704             :             ! The mean of chi is at saturation and does not vary in the ith PDF component
        2705 35394308879 :             ice_supersat_frac(i,k) = zero
        2706 35394308879 :             rc(i,k)         = zero
        2707             :
        2708  2975527959 :           elseif ( mean_chi(i,k)-chi_at_ice_sat > max_num_stdevs * stdev_chi(i,k) ) then
        2709             :
        2710             :             ! The mean of chi is multiple standard deviations above the saturation point.
        2711             :             ! Thus, all cloud in the ith PDF component.
        2712  1695417368 :             ice_supersat_frac(i,k) = one
        2713  1695417368 :             rc(i,k)         = mean_chi(i,k)-chi_at_ice_sat
        2714             :
        2715             :           else
        2716             :
        2717             :             ! The mean of chi is within max_num_stdevs of the saturation point.
        2718             :             ! Thus, layer is partly cloudy, requires calculation.
        2719             :
        2720  1280110591 :             zeta = (mean_chi(i,k)-chi_at_ice_sat) / stdev_chi(i,k)
        2721             :
        2722  1280110591 :             ice_supersat_frac(i,k) = one_half * ( one + erf( zeta * invrs_sqrt_2 )  )
        2723             :
        2724             :             rc(i,k) = (mean_chi(i,k)-chi_at_ice_sat) * ice_supersat_frac(i,k) &
        2725  1280110591 :                       + stdev_chi(i,k) * exp( - one_half * zeta**2 ) * invrs_sqrt_2pi
        2726             :
        2727             :           end if
        2728             :
        2729             :         end if
        2730             :
        2731             :       end do
        2732             :     end do
        2733             :     !$acc end parallel loop
        2734             :
        2735             :     !$acc end data
        2736             :
        2737             :     return
        2738             :
        2739             :   end subroutine calc_ice_cloud_frac_component
        2740             :
        2741             :   !=============================================================================
        2742    35740224 :   subroutine calc_xprcp_component( nz, ngrdcol,                                     & ! In
        2743    35740224 :                                    wm, rtm, thlm, um, vm, rcm,                      & ! In
        2744    35740224 :                                    w_i, rt_i,                                       & ! In
        2745    35740224 :                                    thl_i, u_i, v_i,                                 & ! In
        2746    35740224 :                                    varnce_w_i, chi_i,                               & ! In
        2747    35740224 :                                    stdev_chi_i, stdev_eta_i,                        & ! In
        2748    35740224 :                                    corr_w_chi_i, corr_chi_eta_i,                    & ! In
        2749             : !                                  corr_u_w_i, corr_v_w_i,                          & ! In
        2750    35740224 :                                    crt_i, cthl_i,                                   & ! In
        2751    35740224 :                                    rc_i, cloud_frac_i, iiPDF_type,                  & ! In
        2752    35740224 :                                    wprcp_contrib_comp_i, wp2rcp_contrib_comp_i,     & ! Out
        2753    35740224 :                                    rtprcp_contrib_comp_i, thlprcp_contrib_comp_i,   & ! Out
        2754    35740224 :                                    uprcp_contrib_comp_i, vprcp_contrib_comp_i )       ! Out
        2755             :
        2756             :     ! Description:
        2757             :     ! Calculates the contribution to <w'rc'>, <w'^2 rc'>, <rt'rc'>, and
        2758             :     ! <thl'rc'> from the ith PDF component.
        2759             :     !
        2760             :     !
        2761             :     ! <w'rc'>
        2762             :     ! -------
        2763             :     !
        2764             :     ! The value of <w'rc'> is calculated by integrating over the PDF:
        2765             :     !
        2766             :     ! <w'rc'>
        2767             :     ! = INT(-inf:inf) INT(-inf:inf) INT(-inf:inf)
        2768             :     !   ( w - <w> ) ( rc - <rc> ) P(w,rt,thl) dthl drt dw;
        2769             :     !
        2770             :     ! where <w> is the overall mean of w, <rc> is the overall mean of rc, and
        2771             :     ! P(w,rt,thl) is a two-component trivariate normal distribution of w, rt,
        2772             :     ! and thl.  This equation is rewritten as:
        2773             :     !
        2774             :     ! <w'rc'>
        2775             :     ! = mixt_frac
        2776             :     !   * INT(-inf:inf) INT(-inf:inf) INT(-inf:inf)
        2777             :     !     ( w - <w> ) ( rc - <rc> ) P_1(w,rt,thl) dthl drt dw
        2778             :     !   + ( 1 - mixt_frac )
        2779             :     !     * INT(-inf:inf) INT(-inf:inf) INT(-inf:inf)
        2780             :     !       ( w - <w> ) ( rc - <rc> ) P_2(w,rt,thl) dthl drt dw;
        2781             :     !
        2782             :     ! where mixt_frac is the mixture fraction, which is the weight of the 1st
        2783             :     ! PDF component, and where P_1(w,rt,thl) and P_2(w,rt,thl) are the equations
        2784             :     ! for the trivariate normal PDF of w, rt, and thl in the 1st and 2nd PDF
        2785             :     ! components, respectively.  The contribution from the ith PDF component is:
        2786             :     !
        2787             :     ! INT(-inf:inf) INT(-inf:inf) INT(-inf:inf)
        2788             :     ! ( w - <w> ) ( rc - <rc> ) P_i(w,rt,thl) dthl drt dw;
        2789             :     !
        2790             :     ! where P_i(w,rt,thl) is the trivariate normal PDF of w, rt, and thl in the
        2791             :     ! ith PDF component.  The PDF undergoes a PDF transformation in each PDF
        2792             :     ! component, which is a change of variables and a translation, stretching,
        2793             :     ! and rotation of the axes.  The PDF becomes a trivariate normal PDF that is
        2794             :     ! written in terms of w, chi, and eta coordinates.  Cloud water mixing
        2795             :     ! ratio, rc, is written in terms of extended liquid water mixing ratio, chi,
        2796             :     ! such that:
        2797             :     !
        2798             :     ! rc = chi H(chi);
        2799             :     !
        2800             :     ! where H(chi) is the Heaviside step function.  The contribution from the
        2801             :     ! ith PDF component to <w'rc'> can be written as:
        2802             :     !
        2803             :     ! INT(-inf:inf) INT(-inf:inf)
        2804             :     ! ( w - <w> ) ( chi H(chi) - <rc> ) P_i(w,chi) dchi dw;
        2805             :     !
        2806             :     ! where P_i(w,chi) is the bivariate normal PDF of w and chi in the ith PDF
        2807             :     ! component.  The solved equation for the <w'rc'> contribution from the ith
        2808             :     ! PDF component (wprcp_contrib_comp_i) is:
        2809             :     !
        2810             :     ! wprcp_contrib_comp_i
        2811             :     ! = INT(-inf:inf) INT(-inf:inf)
        2812             :     !   ( w - <w> ) ( chi H(chi) - <rc> ) P_i(w,chi) dchi dw
        2813             :     ! = ( mu_w_i - <w> )
        2814             :     !   * ( mu_chi_i * 1/2 * ( 1 + erf( mu_chi_i / ( sqrt(2) * sigma_chi_i ) ) )
        2815             :     !       + 1/sqrt(2*pi) * sigma_chi_i
        2816             :     !         * exp{ - mu_chi_i^2 / ( 2 * sigma_chi_i^2 ) } - <rc> )
        2817             :     !   + corr_w_chi_i * sigma_w_i * sigma_chi_i
        2818             :     !     * 1/2 * ( 1 + erf( mu_chi_i / ( sqrt(2) * sigma_chi_i ) ) );
        2819             :     !
        2820             :     ! where mu_w_i is the mean of w in the ith PDF component, mu_chi_i is the
        2821             :     ! mean of chi in the ith PDF component, sigma_w_i is the standard deviation
        2822             :     ! of w in the ith PDF component, sigma_chi_i is the standard deviation of
        2823             :     ! chi in the ith PDF component, and corr_w_chi_i is the correlation of w and
        2824             :     ! chi in the ith PDF component.
        2825             :     !
        2826             :     ! Special case:  sigma_chi_i = 0.
        2827             :     !
        2828             :     ! In the special case that sigma_chi_i = 0, chi, as well as rc, are constant
        2829             :     ! in the ith PDF component.  The equation becomes:
        2830             :     !
        2831             :     ! wprcp_contrib_comp_i
        2832             :     ! = | ( mu_w_i - <w> ) * ( mu_chi_i - <rc> ); when mu_chi_i > 0;
        2833             :     !   | ( mu_w_i - <w> ) * ( -<rc> ); when mu_chi_i <= 0.
        2834             :     !
        2835             :     !
        2836             :     ! <w'^2 rc'>
        2837             :     ! ----------
        2838             :     !
        2839             :     ! The value of <w'^2 rc'> is calculated by integrating over the PDF:
        2840             :     !
        2841             :     ! <w'^2 rc'>
        2842             :     ! = INT(-inf:inf) INT(-inf:inf) INT(-inf:inf)
        2843             :     !   ( w - <w> )^2 ( rc - <rc> ) P(w,rt,thl) dthl drt dw.
        2844             :     !
        2845             :     ! This equation is rewritten as:
        2846             :     !
        2847             :     ! <w'^2 rc'>
        2848             :     ! = mixt_frac
        2849             :     !   * INT(-inf:inf) INT(-inf:inf) INT(-inf:inf)
        2850             :     !     ( w - <w> )^2 ( rc - <rc> ) P_1(w,rt,thl) dthl drt dw
        2851             :     !   + ( 1 - mixt_frac )
        2852             :     !     * INT(-inf:inf) INT(-inf:inf) INT(-inf:inf)
        2853             :     !       ( w - <w> )^2 ( rc - <rc> ) P_2(w,rt,thl) dthl drt dw.
        2854             :     !
        2855             :     ! The contribution from the ith PDF component is:
        2856             :     !
        2857             :     ! INT(-inf:inf) INT(-inf:inf) INT(-inf:inf)
        2858             :     ! ( w - <w> )^2 ( rc - <rc> ) P_i(w,rt,thl) dthl drt dw.
        2859             :     !
        2860             :     ! The PDF undergoes a PDF transformation in each PDF component, and becomes
        2861             :     ! a trivariate normal PDF that is written in terms of w, chi, and eta
        2862             :     ! coordinates.  The contribution from the ith PDF component to <w'^2 rc'>
        2863             :     ! can be written as:
        2864             :     !
        2865             :     ! INT(-inf:inf) INT(-inf:inf)
        2866             :     ! ( w - <w> )^2 ( chi H(chi) - <rc> ) P_i(w,chi) dchi dw.
        2867             :     !
        2868             :     ! The solved equation for the <w'^2 rc'> contribution from the ith PDF
        2869             :     ! component (wp2rcp_contrib_comp_i) is:
        2870             :     !
        2871             :     ! wp2rcp_contrib_comp_i
        2872             :     ! = INT(-inf:inf) INT(-inf:inf)
        2873             :     !   ( w - <w> )^2 ( chi H(chi) - <rc> ) P_i(w,chi) dchi dw
        2874             :     ! = ( ( mu_w_i - <w> )^2 + sigma_w_i^2 )
        2875             :     !   * ( mu_chi_i * 1/2 * ( 1 + erf( mu_chi_i / ( sqrt(2) * sigma_chi_i ) ) )
        2876             :     !       + 1/sqrt(2*pi) * sigma_chi_i
        2877             :     !         * exp{ - mu_chi_i^2 / ( 2 * sigma_chi_i^2 ) } - <rc> )
        2878             :     !   + ( mu_w_i - <w> ) * corr_w_chi_i * sigma_w_i * sigma_chi_i
        2879             :     !     * ( 1 + erf( mu_chi_i / ( sqrt(2) * sigma_chi_i ) ) )
        2880             :     !   + 1/sqrt(2*pi) * corr_w_chi_i^2 * sigma_w_i^2 * sigma_chi_i
        2881             :     !     * exp{ - mu_chi_i^2 / ( 2 * sigma_chi_i^2 ) }.
        2882             :     !
        2883             :     ! Special case:  sigma_chi_i = 0.
        2884             :     !
        2885             :     ! In the special case that sigma_chi_i = 0, chi, as well as rc, are constant
        2886             :     ! in the ith PDF component.  The equation becomes:
        2887             :     !
        2888             :     ! wp2rcp_contrib_comp_i
        2889             :     ! = | ( ( mu_w_i - <w> )^2 + sigma_w_i^2 ) * ( mu_chi_i - <rc> );
        2890             :     !   |     when mu_chi_i > 0;
        2891             :     !   | ( ( mu_w_i - <w> )^2 + sigma_w_i^2 ) * ( -<rc> );
        2892             :     !   |     when mu_chi_i <= 0.
        2893             :     !
        2894             :     !
        2895             :     ! <rt'rc'>
        2896             :     ! --------
        2897             :     !
        2898             :     ! The value of <rt'rc'> is calculated by integrating over the PDF:
        2899             :     !
        2900             :     ! <rt'rc'>
        2901             :     ! = INT(-inf:inf) INT(-inf:inf)
        2902             :     !   ( rt - <rt> ) ( rc - <rc> ) P(rt,thl) dthl drt;
        2903             :     !
        2904             :     ! where <rt> is the overall mean of rt, and where P(rt,thl) is a
        2905             :     ! two-component bivariate normal distribution of rt and thl.  This equation
        2906             :     ! is rewritten as:
        2907             :     !
        2908             :     ! <rt'rc'>
        2909             :     ! = mixt_frac
        2910             :     !   * INT(-inf:inf) INT(-inf:inf)
        2911             :     !     ( rt - <rt> ) ( rc - <rc> ) P_1(rt,thl) dthl drt
        2912             :     !   + ( 1 - mixt_frac )
        2913             :     !     * INT(-inf:inf) INT(-inf:inf)
        2914             :     !       ( rt - <rt> ) ( rc - <rc> ) P_2(rt,thl) dthl drt;
        2915             :     !
        2916             :     ! where P_1(rt,thl) and P_2(rt,thl) are the equations for the bivariate
        2917             :     ! normal PDF of rt and thl in the 1st and 2nd PDF components, respectively.
        2918             :     ! The contribution from the ith PDF component is:
        2919             :     !
        2920             :     ! INT(-inf:inf) INT(-inf:inf)
        2921             :     ! ( rt - <rt> ) ( rc - <rc> ) P_i(rt,thl) dthl drt;
        2922             :     !
        2923             :     ! where P_i(rt,thl) is the bivariate normal PDF of rt and thl in the ith PDF
        2924             :     ! component.  The PDF undergoes a PDF transformation in each PDF component,
        2925             :     ! and becomes a bivariate normal PDF that is written in terms of chi and
        2926             :     ! eta coordinates.  Total water mixing ratio, rt, is rewritten in terms of
        2927             :     ! chi and eta by:
        2928             :     !
        2929             :     ! rt = mu_rt_i
        2930             :     !      + ( ( eta - mu_eta_i ) + ( chi - mu_chi_i ) ) / ( 2 * crt_i );
        2931             :     !
        2932             :     ! where mu_rt_i is the mean of rt in the ith PDF component, mu_eta_i is the
        2933             :     ! mean of eta in the ith PDF component, and crt_i is a coefficient on rt in
        2934             :     ! the chi/eta transformation equations.  The contribution from the ith PDF
        2935             :     ! component to <rt'rc'> can be written as:
        2936             :     !
        2937             :     ! INT(-inf:inf) INT(-inf:inf)
        2938             :     ! ( mu_rt_i - <rt> + ( eta - mu_eta_i ) / ( 2 * crt_i )
        2939             :     !   + ( chi - mu_chi_i ) / ( 2 * crt_i ) )
        2940             :     ! * ( chi H(chi) - <rc> ) P_i(chi,eta) deta dchi;
        2941             :     !
        2942             :     ! where P_i(chi,eta) is the bivariate normal PDF of chi and eta in the ith
        2943             :     ! PDF component.  The solved equation for the <rt'rc'> contribution from the
        2944             :     ! ith PDF component (rtprcp_contrib_comp_i) is:
        2945             :     !
        2946             :     ! rtprcp_contrib_comp_i
        2947             :     ! = INT(-inf:inf) INT(-inf:inf)
        2948             :     !   ( mu_rt_i - <rt> + ( eta - mu_eta_i ) / ( 2 * crt_i )
        2949             :     !     + ( chi - mu_chi_i ) / ( 2 * crt_i ) )
        2950             :     !   * ( chi H(chi) - <rc> ) P_i(chi,eta) deta dchi
        2951             :     ! = ( mu_rt_i - <rt> )
        2952             :     !   * ( mu_chi_i * 1/2 * ( 1 + erf( mu_chi_i / ( sqrt(2) * sigma_chi_i ) ) )
        2953             :     !       + 1/sqrt(2*pi) * sigma_chi_i
        2954             :     !         * exp{ - mu_chi_i^2 / ( 2 * sigma_chi_i^2 ) } - <rc> )
        2955             :     !   + ( corr_chi_eta_i * sigma_eta_i + sigma_chi_i ) / ( 2 * crt_i )
        2956             :     !     * sigma_chi_i
        2957             :     !     * 1/2 * ( 1 + erf( mu_chi_i / ( sqrt(2) * sigma_chi_i ) ) );
        2958             :     !
        2959             :     ! where sigma_eta_i is the standard deviation of eta in the ith PDF
        2960             :     ! component and corr_chi_eta_i is the correlation of chi and eta in the ith
        2961             :     ! PDF component.
        2962             :     !
        2963             :     ! Special case:  sigma_chi_i = 0.
        2964             :     !
        2965             :     ! In the special case that sigma_chi_i = 0, chi, as well as rc, are constant
        2966             :     ! in the ith PDF component.  The equation becomes:
        2967             :     !
        2968             :     ! rtprcp_contrib_comp_i
        2969             :     ! = | ( mu_rt_i - <rt> ) * ( mu_chi_i - <rc> ); when mu_chi_i > 0;
        2970             :     !   | ( mu_rt_i - <rt> ) * ( -<rc> ); when mu_chi_i <= 0.
        2971             :     !
        2972             :     !
        2973             :     ! <thl'rc'>
        2974             :     ! ---------
        2975             :     !
        2976             :     ! The value of <thl'rc'> is calculated by integrating over the PDF:
        2977             :     !
        2978             :     ! <thl'rc'>
        2979             :     ! = INT(-inf:inf) INT(-inf:inf)
        2980             :     !   ( thl - <thl> ) ( rc - <rc> ) P(rt,thl) dthl drt;
        2981             :     !
        2982             :     ! where <thl> is the overall mean of thl.  This equation is rewritten as:
        2983             :     !
        2984             :     ! <thl'rc'>
        2985             :     ! = mixt_frac
        2986             :     !   * INT(-inf:inf) INT(-inf:inf)
        2987             :     !     ( thl - <thl> ) ( rc - <rc> ) P_1(rt,thl) dthl drt
        2988             :     !   + ( 1 - mixt_frac )
        2989             :     !     * INT(-inf:inf) INT(-inf:inf)
        2990             :     !       ( thl - <thl> ) ( rc - <rc> ) P_2(rt,thl) dthl drt.
        2991             :     !
        2992             :     ! The contribution from the ith PDF component is:
        2993             :     !
        2994             :     ! INT(-inf:inf) INT(-inf:inf)
        2995             :     ! ( thl - <thl> ) ( rc - <rc> ) P_i(rt,thl) dthl drt.
        2996             :     !
        2997             :     ! The PDF undergoes a PDF transformation in each PDF component, and becomes
        2998             :     ! a bivariate normal PDF that is written in terms of chi and eta
        2999             :     ! coordinates.  Liquid water potential temperature, thl, is rewritten in
        3000             :     ! terms of chi and eta by:
        3001             :     !
        3002             :     ! thl = mu_thl_i
        3003             :     !       + ( ( eta - mu_eta_i ) - ( chi - mu_chi_i ) ) / ( 2 * cthl_i );
        3004             :     !
        3005             :     ! where mu_thl_i is the mean of thl in the ith PDF component and cthl_i is a
        3006             :     ! coefficient on thl in the chi/eta transformation equations.  The
        3007             :     ! contribution from the ith PDF component to <thl'rc'> can be written as:
        3008             :     !
        3009             :     ! INT(-inf:inf) INT(-inf:inf)
        3010             :     ! ( mu_thl_i - <thl> + ( eta - mu_eta_i ) / ( 2 * cthl_i )
        3011             :     !   - ( chi - mu_chi_i ) / ( 2 * cthl_i ) )
        3012             :     ! * ( chi H(chi) - <rc> ) P_i(chi,eta) deta dchi.
        3013             :     !
        3014             :     ! The solved equation for the <thl'rc'> contribution from the ith PDF
        3015             :     ! component (thlprcp_contrib_comp_i) is:
        3016             :     !
        3017             :     ! thlprcp_contrib_comp_i
        3018             :     ! = INT(-inf:inf) INT(-inf:inf)
        3019             :     !   ( mu_thl_i - <thl> + ( eta - mu_eta_i ) / ( 2 * cthl_i )
        3020             :     !     - ( chi - mu_chi_i ) / ( 2 * cthl_i ) )
        3021             :     !   * ( chi H(chi) - <rc> ) P_i(chi,eta) deta dchi
        3022             :     ! = ( mu_thl_i - <thl> )
        3023             :     !   * ( mu_chi_i * 1/2 * ( 1 + erf( mu_chi_i / ( sqrt(2) * sigma_chi_i ) ) )
        3024             :     !       + 1/sqrt(2*pi) * sigma_chi_i
        3025             :     !         * exp{ - mu_chi_i^2 / ( 2 * sigma_chi_i^2 ) } - <rc> )
        3026             :     !   + ( corr_chi_eta_i * sigma_eta_i - sigma_chi_i ) / ( 2 * cthl_i )
        3027             :     !     * sigma_chi_i
        3028             :     !     * 1/2 * ( 1 + erf( mu_chi_i / ( sqrt(2) * sigma_chi_i ) ) ).
        3029             :     !
        3030             :     ! Special case:  sigma_chi_i = 0.
        3031             :     !
        3032             :     ! In the special case that sigma_chi_i = 0, chi, as well as rc, are constant
        3033             :     ! in the ith PDF component.  The equation becomes:
        3034             :     !
        3035             :     ! thlprcp_contrib_comp_i
        3036             :     ! = | ( mu_thl_i - <thl> ) * ( mu_chi_i - <rc> ); when mu_chi_i > 0;
        3037             :     !   | ( mu_thl_i - <thl> ) * ( -<rc> ); when mu_chi_i <= 0.
        3038             :     !
        3039             :     !
        3040             :     ! Use equations for PDF component cloud fraction cloud water mixing ratio
        3041             :     ! -----------------------------------------------------------------------
        3042             :     !
        3043             :     ! The equation for cloud fraction in the ith PDF component, fc_i, is:
        3044             :     !
        3045             :     ! fc_i = 1/2 * ( 1 + erf( mu_chi_i / ( sqrt(2) * sigma_chi_i ) ) ).
        3046             :     !
        3047             :     ! In the special case that sigma_chi_i = 0, the equation becomes:
        3048             :     !
        3049             :     ! fc_i = | 1; when mu_chi_i > 0;
        3050             :     !        | 0; when mu_chi_i <= 0.
        3051             :     !
        3052             :     ! The equation for mean cloud water mixing ratio in the ith PDF component,
        3053             :     ! rc_i, is:
        3054             :     !
        3055             :     ! rc_i
        3056             :     ! = mu_chi_i * 1/2 * ( 1 + erf( mu_chi_i / ( sqrt(2) * sigma_chi_i ) ) )
        3057             :     !   + 1/sqrt(2*pi) * sigma_chi_i
        3058             :     !     * exp{ - mu_chi_i^2 / ( 2 * sigma_chi_i^2 ) }
        3059             :     ! = mu_chi_i * fc_i
        3060             :     !   + 1/sqrt(2*pi) * sigma_chi_i
        3061             :     !     * exp{ - mu_chi_i^2 / ( 2 * sigma_chi_i^2 ) }.
        3062             :     !
        3063             :     ! In the special case that sigma_chi_i = 0, the equation becomes:
        3064             :     !
        3065             :     ! rc_i = | mu_chi_i; when mu_chi_i > 0;
        3066             :     !        | 0; when mu_chi_i <= 0.
        3067             :     !
        3068             :     ! The above equations can be substituted into the equations for
        3069             :     ! wprcp_contrib_comp_i, wp2rcp_contrib_comp_i, rtprcp_contrib_comp_i, and
        3070             :     ! thlprcp_contrib_comp_i.  The new equations are:
        3071             :     !
        3072             :     ! wprcp_contrib_comp_i
        3073             :     ! = ( mu_w_i - <w> ) * ( rc_i - <rc> )
        3074             :     !   + corr_w_chi_i * sigma_w_i * sigma_chi_i * fc_i;
        3075             :     !
        3076             :     ! wp2rcp_contrib_comp_i
        3077             :     ! = ( ( mu_w_i - <w> )^2 + sigma_w_i^2 ) * ( rc_i - <rc> )
        3078             :     !   + 2 * ( mu_w_i - <w> ) * corr_w_chi_i * sigma_w_i * sigma_chi_i * fc_i
        3079             :     !   + 1/sqrt(2*pi) * corr_w_chi_i^2 * sigma_w_i^2 * sigma_chi_i
        3080             :     !     * exp{ - mu_chi_i^2 / ( 2 * sigma_chi_i^2 ) };
        3081             :     !
        3082             :     ! rtprcp_contrib_comp_i
        3083             :     ! = ( mu_rt_i - <rt> ) * ( rc_i - <rc> )
        3084             :     !   + ( corr_chi_eta_i * sigma_eta_i + sigma_chi_i ) / ( 2 * crt_i )
        3085             :     !     * sigma_chi_i * fc_i; and
        3086             :     !
        3087             :     ! thlprcp_contrib_comp_i
        3088             :     ! = ( mu_thl_i - <thl> ) * ( rc_i - <rc> )
        3089             :     !   + ( corr_chi_eta_i * sigma_eta_i - sigma_chi_i ) / ( 2 * cthl_i )
        3090             :     !     * sigma_chi_i * fc_i.
        3091             :     !
        3092             :     ! While the above equations reduce to their listed versions in the special
        3093             :     ! case that sigma_chi_i = 0, those versions are faster to calculate.  When
        3094             :     ! mu_chi_i > 0, they are:
        3095             :     !
        3096             :     ! wprcp_contrib_comp_i = ( mu_w_i - <w> ) * ( mu_chi_i - <rc> );
        3097             :     ! wp2rcp_contrib_comp_i
        3098             :     ! = ( ( mu_w_i - <w> )^2 + sigma_w_i^2 ) * ( mu_chi_i - <rc> );
        3099             :     ! rtprcp_contrib_comp_i = ( mu_rt_i - <rt> ) * ( mu_chi_i - <rc> ); and
        3100             :     ! thlprcp_contrib_comp_i = ( mu_thl_i - <thl> ) * ( mu_chi_i - <rc> );
        3101             :     !
        3102             :     ! and when mu_chi_i <= 0, they are:
        3103             :     !
        3104             :     ! wprcp_contrib_comp_i = - ( mu_w_i - <w> ) * <rc>;
        3105             :     ! wp2rcp_contrib_comp_i = - ( ( mu_w_i - <w> )^2 + sigma_w_i^2 ) * <rc>;
        3106             :     ! rtprcp_contrib_comp_i = - ( mu_rt_i - <rt> ) * <rc>; and
        3107             :     ! thlprcp_contrib_comp_i = - ( mu_thl_i - <thl> ) * <rc>.
        3108             :
        3109             :     ! References:
        3110             :     !-----------------------------------------------------------------------
        3111             :
        3112             :     use grid_class, only: &
        3113             :         grid ! Type
        3114             :
        3115             :     use constants_clubb, only: &
        3116             :         sqrt_2pi,       & ! Variable(s)
        3117             :         two,            &
        3118             :         zero,           &
        3119             :         chi_tol
        3120             :
        3121             :     use clubb_precision, only: &
        3122             :         core_rknd    ! Variable(s)
        3123             :
        3124             :     implicit none
        3125             :
        3126             :     integer, intent(in) :: &
        3127             :       ngrdcol,  & ! Number of grid columns
        3128             :       nz          ! Number of vertical level
        3129             :
        3130             :     ! Input Variables
        3131             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
        3132             :       wm,             & ! Mean of w (overall)                          [m/s]
        3133             :       rtm,            & ! Mean of rt (overall)                         [kg/kg]
        3134             :       thlm,           & ! Mean of thl (overall)                        [K]
        3135             :       um,             & ! Mean of eastward wind (overall)              [m/s]
        3136             :       vm,             & ! Mean of northward wind (overall)             [m/s]
        3137             :       rcm,            & ! Mean of rc (overall)                         [kg/kg]
        3138             :       w_i,            & ! Mean of w (ith PDF component)                [m/s]
        3139             :       rt_i,           & ! Mean of rt (ith PDF component)               [kg/kg]
        3140             :       thl_i,          & ! Mean of thl (ith PDF component)              [K]
        3141             :       u_i,            & ! Mean of eastward wind (ith PDF component)    [m/s]
        3142             :       v_i,            & ! Mean of northward wind (ith PDF component)   [m/s]
        3143             :       varnce_w_i,     & ! Variance of w (ith PDF component)            [m^2/s^2]
        3144             :       chi_i,          & ! Mean of chi (ith PDF component)              [kg/kg]
        3145             :       stdev_chi_i,    & ! Standard deviation of chi (ith PDF comp.)    [kg/kg]
        3146             :       stdev_eta_i,    & ! Standard deviation of eta (ith PDF comp.)    [kg/kg]
        3147             :       corr_w_chi_i,   & ! Correlation of w and chi (ith PDF component) [-]
        3148             :       corr_chi_eta_i, & ! Correlation of chi and eta (ith PDF comp.)   [-]
        3149             : !     corr_u_w_i,     & ! Correlation of u and w (ith PDF component)   [-]
        3150             : !     corr_v_w_i,     & ! Correlation of v and w (ith PDF component)   [-]
        3151             :       crt_i,          & ! Coef. on rt in chi/eta eqns. (ith PDF comp.) [-]
        3152             :       cthl_i,         & ! Coef. on thl: chi/eta eqns. (ith PDF comp.)  [kg/kg/K]
        3153             :       rc_i,           & ! Mean of rc (ith PDF component)               [kg/kg]
        3154             :       cloud_frac_i      ! Cloud fraction (ith PDF component)           [-]
        3155             :
        3156             :     integer, intent(in) :: &
        3157             :       iiPDF_type    ! Selected option for the two-component normal (double
        3158             :                     ! Gaussian) PDF type to use for the w, rt, and theta-l (or
        3159             :                     ! w, chi, and eta) portion of CLUBB's multivariate,
        3160             :                     ! two-component PDF.
        3161             :
        3162             :     ! Output Variables
        3163             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
        3164             :       wprcp_contrib_comp_i,   & ! <w'rc'> contrib. (ith PDF comp.)  [m/s(kg/kg)]
        3165             :       wp2rcp_contrib_comp_i,  & ! <w'^2rc'> contrib. (ith comp) [m^2/s^2(kg/kg)]
        3166             :       rtprcp_contrib_comp_i,  & ! <rt'rc'> contrib. (ith PDF comp.)  [kg^2/kg^2]
        3167             :       thlprcp_contrib_comp_i, & ! <thl'rc'> contrib. (ith PDF comp.)  [K(kg/kg)]
        3168             :       uprcp_contrib_comp_i,   & ! <u'rc'> contrib. (ith PDF comp.)  [m/s(kg/kg)]
        3169             :       vprcp_contrib_comp_i      ! <v'rc'> contrib. (ith PDF comp.)  [m/s(kg/kg)]
        3170             :
        3171             :     ! Local Variables
        3172             :     integer :: i, k
        3173             :
        3174             :     ! ---------------------- Begin Code ------------------
        3175             :
        3176             :     ! Changing these conditionals may result in inconsistencies with the conditional
        3177             :     ! statements located in calc_cloud_frac_component
        3178             :     !$acc parallel loop gang vector collapse(2) default(present)
        3179  3073659264 :     do k = 1, nz
        3180 50761923264 :       do i = 1, ngrdcol
        3181             :
        3182 47688264000 :         wprcp_contrib_comp_i(i,k) = ( w_i(i,k) - wm(i,k) ) * ( rc_i(i,k) - rcm(i,k) )
        3183             :
        3184             :         wp2rcp_contrib_comp_i(i,k) = ( ( w_i(i,k) - wm(i,k) )**2 + varnce_w_i(i,k) ) &
        3185 47688264000 :                                      * ( rc_i(i,k) - rcm(i,k) )
        3186             :
        3187             :         rtprcp_contrib_comp_i(i,k) = ( rt_i(i,k) - rtm(i,k) ) * ( rc_i(i,k) - rcm(i,k) ) &
        3188             :                                 + ( corr_chi_eta_i(i,k) * stdev_eta_i(i,k) + stdev_chi_i(i,k) ) &
        3189 47688264000 :                                   / ( two * crt_i(i,k) ) * stdev_chi_i(i,k) * cloud_frac_i(i,k)
        3190             :
        3191             :         thlprcp_contrib_comp_i(i,k) = ( thl_i(i,k) - thlm(i,k) ) * ( rc_i(i,k) - rcm(i,k) ) &
        3192             :                                  + ( corr_chi_eta_i(i,k) * stdev_eta_i(i,k) - stdev_chi_i(i,k) ) &
        3193 47688264000 :                                    / ( two * cthl_i(i,k) ) * stdev_chi_i(i,k) * cloud_frac_i(i,k)
        3194             :
        3195 47688264000 :         uprcp_contrib_comp_i(i,k) = ( u_i(i,k) - um(i,k) ) * ( rc_i(i,k) - rcm(i,k) )
        3196             :
        3197 50726183040 :         vprcp_contrib_comp_i(i,k) = ( v_i(i,k) - vm(i,k) ) * ( rc_i(i,k) - rcm(i,k) )
        3198             :
        3199             :       end do
        3200             :     end do
        3201             :     !$acc end parallel loop
        3202             :
        3203             :     ! If iiPDF_type isn't iiPDF_ADG1, iiPDF_ADG2, or iiPDF_new_hybrid, so
        3204             :     ! corr_w_chi_i /= 0 (and perhaps corr_u_w_i /= 0).
        3205    35740224 :     if ( .not. ( iiPDF_type == iiPDF_ADG1 .or. iiPDF_type == iiPDF_ADG2 &
        3206             :                  .or. iiPDF_type == iiPDF_new_hybrid ) ) then
        3207             :
        3208             :         ! Chi varies significantly in the ith PDF component (stdev_chi > chi_tol)
        3209             :         ! and there is some cloud (0 < cloud_frac <= 1)
        3210           0 :         do k = 1, nz
        3211           0 :           do i = 1, ngrdcol
        3212           0 :             if ( stdev_chi_i(i,k) > chi_tol .and. cloud_frac_i(i,k) > zero ) then
        3213             :
        3214             :               wprcp_contrib_comp_i(i,k) = wprcp_contrib_comp_i(i,k) &
        3215             :                                           + corr_w_chi_i(i,k) * sqrt( varnce_w_i(i,k) ) &
        3216           0 :                                             * stdev_chi_i(i,k) * cloud_frac_i(i,k)
        3217             :
        3218             :               wp2rcp_contrib_comp_i(i,k) = wp2rcp_contrib_comp_i(i,k) &
        3219             :                                            + two * ( w_i(i,k) - wm(i,k) ) * corr_w_chi_i(i,k) &
        3220             :                                              * sqrt( varnce_w_i(i,k) ) * stdev_chi_i(i,k) &
        3221             :                                              * cloud_frac_i(i,k) &
        3222             :                                            + corr_w_chi_i(i,k)**2 * varnce_w_i(i,k) &
        3223             :                                              * stdev_chi_i(i,k) &
        3224             :                                              * exp( -chi_i(i,k)**2 / ( two*stdev_chi_i(i,k)**2 ) ) &
        3225           0 :                                                / sqrt_2pi
        3226             :
        3227             :             ! In principle, uprcp_contrib_comp_i might depend on corr_u_w_i here.
        3228             :           end if
        3229             :         end do
        3230             :       end do
        3231             :     end if
        3232             :
        3233    35740224 :     return
        3234             :
        3235             :   end subroutine calc_xprcp_component
        3236             :
        3237             :   !=============================================================================
        3238           0 :   subroutine calc_w_up_in_cloud( nz, ngrdcol, &
        3239           0 :                                  mixt_frac, cloud_frac_1, cloud_frac_2, &
        3240           0 :                                  w_1, w_2, varnce_w_1, varnce_w_2, &
        3241           0 :                                  w_up_in_cloud, w_down_in_cloud, &
        3242           0 :                                  cloudy_updraft_frac, cloudy_downdraft_frac )
        3243             :
        3244             :     ! Description:
        3245             :     ! Subroutine that computes the mean cloudy updraft (and also calculates
        3246             :     ! the mean cloudy downdraft).
        3247             :     !
        3248             :     ! In order to activate aerosol, we'd like to feed the activation scheme
        3249             :     ! a vertical velocity that's representative of cloudy updrafts. For skewed
        3250             :     ! layers, like cumulus layers, this might be an improvement over the square
        3251             :     ! root of wp2 that's currently used. At the same time, it would be simpler
        3252             :     ! and less expensive than feeding SILHS samples into the aerosol code
        3253             :     ! (see larson-group/e3sm#19 and larson-group/e3sm#26).
        3254             :     !
        3255             :     ! The formulas are only valid for certain PDFs in CLUBB (ADG1, ADG2,
        3256             :     ! new hybrid), hence we omit calculation if another PDF type is used.
        3257             :     !
        3258             :     ! References: https://www.overleaf.com/project/614a136d47846639af22ae34
        3259             :     !----------------------------------------------------------------------
        3260             :
        3261             :     use constants_clubb, only: &
        3262             :         sqrt_2pi,       & ! sqrt(2*pi)
        3263             :         sqrt_2,         & ! sqrt(2)
        3264             :         one,            & ! 1
        3265             :         one_half,       & ! 1/2
        3266             :         zero,           & ! 0
        3267             :         max_num_stdevs, &
        3268             :         eps
        3269             :
        3270             :     use clubb_precision, only: &
        3271             :         core_rknd     ! Precision
        3272             :
        3273             :     implicit none
        3274             :
        3275             :     integer, intent(in) :: &
        3276             :       ngrdcol,  & ! Number of grid columns
        3277             :       nz          ! Number of vertical level
        3278             :
        3279             :     !----------- Input Variables -----------
        3280             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
        3281             :       mixt_frac, &      ! mixture fraction                             [-]
        3282             :       cloud_frac_1, &   ! cloud fraction (1st PDF component)           [-]
        3283             :       cloud_frac_2, &   ! cloud fraction (2nd PDF component)           [-]
        3284             :       w_1, &            ! upward velocity (1st PDF component)          [m/s]
        3285             :       w_2, &            ! upward velocity (2nd PDF component)          [m/s]
        3286             :       varnce_w_1, &     ! standard deviation of w (1st PDF component)  [m^2/s^2]
        3287             :       varnce_w_2        ! standard deviation of w (2nd PDF component)  [m^2/s^2]
        3288             :
        3289             :     !----------- Output Variables -----------
        3290             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
        3291             :       w_up_in_cloud,         & ! mean cloudy updraft speed             [m/s]
        3292             :       w_down_in_cloud,       & ! mean cloudy downdraft speed           [m/s]
        3293             :       cloudy_updraft_frac,   & ! cloudy updraft fraction               [-]
        3294             :       cloudy_downdraft_frac    ! cloudy downdraft fraction             [-]
        3295             :
        3296             :     !----------- Local Variables -----------
        3297             :     real( kind = core_rknd ) :: &
        3298             :       w_up_1, w_up_2, &        ! integral of w and Heaviside fnc, where w > 0
        3299             :       w_down_1, w_down_2, &    ! integral of w and Heaviside fnc, where w < 0
        3300             :       stdev_w_1, stdev_w_2, &  ! Standard deviation of w
        3301             :       ratio_w_1, &             ! mu_w_1 / ( sqrt(2) * sigma_w_1 )
        3302             :       ratio_w_2, &             ! mu_w_2 / ( sqrt(2) * sigma_w_2 )
        3303             :       erf_ratio_w_1, &         ! erf( ratio_w_1 )
        3304             :       erf_ratio_w_2, &         ! erf( ratio_w_2 )
        3305             :       exp_neg_ratio_w_1_sqd, & ! exp( -ratio_w_1^2 )
        3306             :       exp_neg_ratio_w_2_sqd, & ! exp( -ratio_w_2^2 )
        3307             :       updraft_frac_1, &        ! Fraction of 1st PDF comp. where w > 0
        3308             :       updraft_frac_2, &        ! Fraction of 2nd PDF comp. where w > 0
        3309             :       downdraft_frac_1, &      ! Fraction of 1st PDF comp. where w < 0
        3310             :       downdraft_frac_2         ! Fraction of 2nd PDF comp. where w < 0
        3311             :
        3312             :     integer :: i, k
        3313             :
        3314             :     !$acc parallel loop gang vector collapse(2) default(present)
        3315           0 :     do k = 1, nz
        3316           0 :       do i = 1, ngrdcol
        3317             :
        3318           0 :         stdev_w_1 = sqrt(varnce_w_1(i,k))
        3319           0 :         stdev_w_2 = sqrt(varnce_w_2(i,k))
        3320             :
        3321             :         ! Calculate quantities in the 1st PDF component.
        3322           0 :         if ( w_1(i,k) > max_num_stdevs * stdev_w_1 ) then
        3323             :
        3324             :            ! The mean of w in the 1st PDF component is more than
        3325             :            ! max_num_stdevs standard deviations above 0.
        3326             :            ! The entire 1st PDF component is found in an updraft (w > 0).
        3327             :            w_up_1 = w_1(i,k)
        3328             :            updraft_frac_1 = one
        3329             :            w_down_1 = zero
        3330             :            downdraft_frac_1 = zero
        3331             :
        3332           0 :         elseif ( w_1(i,k) < - max_num_stdevs * stdev_w_1 ) then
        3333             :
        3334             :            ! The mean of w in the 1st PDF component is more than
        3335             :            ! max_num_stdevs standard deviations below 0.
        3336             :            ! The entire 1st PDF component is found in a downdraft (w < 0).
        3337             :            w_up_1 = zero
        3338             :            updraft_frac_1 = zero
        3339             :            w_down_1 = w_1(i,k)
        3340             :            downdraft_frac_1 = one
        3341             :
        3342             :         else
        3343             :
        3344             :            ! The 1st PDF component contains both updraft and downdraft.
        3345           0 :            ratio_w_1 = w_1(i,k) / ( sqrt_2 * max(eps, stdev_w_1) )
        3346           0 :            erf_ratio_w_1 = erf( ratio_w_1 )
        3347           0 :            exp_neg_ratio_w_1_sqd = exp( -ratio_w_1**2 )
        3348             :
        3349             :            w_up_1 &
        3350             :            = one_half * w_1(i,k) * ( one + erf_ratio_w_1 ) &
        3351           0 :              + ( stdev_w_1 / sqrt_2pi ) * exp_neg_ratio_w_1_sqd
        3352             :
        3353           0 :            updraft_frac_1 = one_half * ( one + erf_ratio_w_1 )
        3354             :
        3355             :            w_down_1 &
        3356             :            = one_half * w_1(i,k) * ( one - erf_ratio_w_1 ) &
        3357           0 :              - ( stdev_w_1 / sqrt_2pi ) * exp_neg_ratio_w_1_sqd
        3358             :
        3359             :            !downdraft_frac_1 = one_half * ( one - erf_ratio_w_1 )
        3360           0 :            downdraft_frac_1 = one - updraft_frac_1
        3361             :
        3362             :         endif
        3363             :
        3364             :         ! Calculate quantities in the 2nd PDF component.
        3365           0 :         if ( w_2(i,k) > max_num_stdevs * stdev_w_2 ) then
        3366             :
        3367             :            ! The mean of w in the 2nd PDF component is more than
        3368             :            ! max_num_stdevs standard deviations above 0.
        3369             :            ! The entire 2nd PDF component is found in an updraft (w > 0).
        3370             :            w_up_2 = w_2(i,k)
        3371             :            updraft_frac_2 = one
        3372             :            w_down_2 = zero
        3373             :            downdraft_frac_2 = zero
        3374             :
        3375           0 :         elseif ( w_2(i,k) < - max_num_stdevs * stdev_w_2 ) then
        3376             :
        3377             :            ! The mean of w in the 2nd PDF component is more than
        3378             :            ! max_num_stdevs standard deviations below 0.
        3379             :            ! The entire 2nd PDF component is found in a downdraft (w < 0).
        3380             :            w_up_2 = zero
        3381             :            updraft_frac_2 = zero
        3382             :            w_down_2 = w_2(i,k)
        3383             :            downdraft_frac_2 = one
        3384             :
        3385             :         else
        3386             :
        3387             :            ! The 2nd PDF component contains both updraft and downdraft.
        3388           0 :            ratio_w_2 = w_2(i,k) / ( sqrt_2 * max(eps, stdev_w_2) )
        3389           0 :            erf_ratio_w_2 = erf( ratio_w_2 )
        3390           0 :            exp_neg_ratio_w_2_sqd = exp( -ratio_w_2**2 )
        3391             :
        3392             :            w_up_2 &
        3393             :            = one_half * w_2(i,k) * ( one + erf_ratio_w_2 ) &
        3394           0 :              + ( stdev_w_2 / sqrt_2pi ) * exp_neg_ratio_w_2_sqd
        3395             :
        3396           0 :            updraft_frac_2 = one_half * ( one + erf_ratio_w_2 )
        3397             :
        3398             :            w_down_2 &
        3399             :            = one_half * w_2(i,k) * ( one - erf_ratio_w_2 ) &
        3400           0 :              - ( stdev_w_2 / sqrt_2pi ) * exp_neg_ratio_w_2_sqd
        3401             :
        3402             :            !downdraft_frac_2 = one_half * ( one - erf_ratio_w_2 )
        3403           0 :            downdraft_frac_2 = one - updraft_frac_2
        3404             :
        3405             :         endif
        3406             :
        3407             :         ! Calculate the total cloudy updraft fraction.
        3408             :         cloudy_updraft_frac(i,k) &
        3409             :         = mixt_frac(i,k) * cloud_frac_1(i,k) * updraft_frac_1 &
        3410           0 :           + ( one - mixt_frac(i,k) ) * cloud_frac_2(i,k) * updraft_frac_2
        3411             :
        3412             :         ! Calculate the total cloudy downdraft fraction.
        3413             :         cloudy_downdraft_frac(i,k) &
        3414             :         = mixt_frac(i,k) * cloud_frac_1(i,k) * downdraft_frac_1 &
        3415           0 :           + ( one - mixt_frac(i,k) ) * cloud_frac_2(i,k) * downdraft_frac_2
        3416             :
        3417             :         ! Calculate the mean vertical velocity found in a cloudy updraft.
        3418             :         w_up_in_cloud(i,k) &
        3419             :         = ( mixt_frac(i,k) * cloud_frac_1(i,k) * w_up_1 &
        3420             :             + ( one - mixt_frac(i,k) ) * cloud_frac_2(i,k) * w_up_2 ) &
        3421           0 :           / max( eps, cloudy_updraft_frac(i,k) )
        3422             :
        3423             :         ! Calculate the mean vertical velocity found in a cloudy downdraft.
        3424             :         w_down_in_cloud(i,k) &
        3425             :         = ( mixt_frac(i,k) * cloud_frac_1(i,k) * w_down_1 &
        3426             :             + ( one - mixt_frac(i,k) ) * cloud_frac_2(i,k) * w_down_2 ) &
        3427           0 :           / max( eps, cloudy_downdraft_frac(i,k) )
        3428             :
        3429             :       end do
        3430             :     end do
        3431             :     !$acc end parallel loop
        3432             :
        3433           0 :     return
        3434             :
        3435             :   end subroutine calc_w_up_in_cloud
        3436             :
        3437             :   !=============================================================================
        3438             :   function interp_var_array( n_points, nz, k, z_vals, var )
        3439             :
        3440             :   ! Description:
        3441             :   !   Interpolates a variable to an array of values about a given level
        3442             :
        3443             :   ! References
        3444             :   !-----------------------------------------------------------------------
        3445             :
        3446             :     use clubb_precision, only: &
        3447             :         core_rknd           ! Constant
        3448             :
        3449             :     implicit none
        3450             :
        3451             :   ! Input Variables
        3452             :     integer, intent(in) :: &
        3453             :       n_points, & ! Number of points to interpolate to (must be odd and >= 3)
        3454             :       nz,       & ! Total number of vertical levels
        3455             :       k           ! Center of interpolation array
        3456             :
        3457             :     real( kind = core_rknd ), dimension(nz), intent(in) :: &
        3458             :       z_vals, &         ! Height at each vertical level           [m]
        3459             :       var               ! Variable values on grid                 [units vary]
        3460             :
        3461             :   ! Output Variables
        3462             :     real( kind = core_rknd ), dimension(n_points) :: &
        3463             :       interp_var_array  ! Interpolated values of variable         [units vary]
        3464             :
        3465             :   ! Local Variables
        3466             :     real( kind = core_rknd ) :: &
        3467             :       dz    ! Distance between vertical levels
        3468             :
        3469             :     real( kind = core_rknd ) :: &
        3470             :       z_val             ! Height at some sub-grid level
        3471             :
        3472             :     integer :: &
        3473             :       i, &                      ! Loop iterator
        3474             :
        3475             :       subgrid_lev_count         ! Number of refined grid points located between
        3476             :                               ! two defined grid levels
        3477             :
        3478             :   !-----------------------------------------------------------------------
        3479             :
        3480             :     !----- Begin Code -----
        3481             :
        3482             :     ! Place a point at each of k-1, k, and k+1.
        3483             :     interp_var_array(1) = var_value_integer_height( nz, k-1, z_vals, var )
        3484             :     interp_var_array((n_points+1)/2) = var_value_integer_height( nz, k, z_vals, var )
        3485             :     interp_var_array(n_points) = var_value_integer_height( nz, k+1, z_vals, var )
        3486             :
        3487             :     subgrid_lev_count = (n_points - 3) / 2
        3488             :
        3489             :     ! Lower half
        3490             :     if ( k == 1 ) then
        3491             :       dz = (z_vals(2) - z_vals(1)) / real( subgrid_lev_count+1, kind=core_rknd )
        3492             :     else
        3493             :       dz = (z_vals(k) - z_vals(k-1)) / real( subgrid_lev_count+1, kind=core_rknd )
        3494             :     end if
        3495             :     do i=1, subgrid_lev_count
        3496             :       z_val = z_vals(k) - real( i, kind=core_rknd ) * dz
        3497             :       interp_var_array(1+i) &
        3498             :       = var_subgrid_interp( nz, k, z_vals, var, z_val, l_below=.true. )
        3499             :     end do
        3500             :
        3501             :     ! Upper half
        3502             :     if ( k == nz ) then
        3503             :       dz = ( z_vals(nz) - z_vals(nz-1) ) / real( subgrid_lev_count+1, kind=core_rknd )
        3504             :     else
        3505             :       dz = ( z_vals(k+1) - z_vals(k) ) / real( subgrid_lev_count+1, kind=core_rknd )
        3506             :     end if
        3507             :     do i=1, (n_points-3)/2
        3508             :       z_val = z_vals(k) + real( i, kind=core_rknd ) * dz
        3509             :       interp_var_array((n_points+1)/2+i) &
        3510             :       = var_subgrid_interp( nz, k, z_vals, var, z_val, l_below=.false. )
        3511             :     end do
        3512             :
        3513             :     return
        3514             :   end function interp_var_array
        3515             :
        3516             :   !=============================================================================
        3517             :   function var_value_integer_height( nz, k, z_vals, var_grid_value ) result( var_value )
        3518             :
        3519             :   ! Description
        3520             :   !   Returns the value of a variable at an integer height between 0 and
        3521             :   !   nz+1 inclusive, using extrapolation when k==0 or k==nz+1
        3522             :
        3523             :   ! References
        3524             :   !-----------------------------------------------------------------------
        3525             :
        3526             :     use clubb_precision, only: &
        3527             :         core_rknd       ! Constant
        3528             :
        3529             :     use interpolation, only: &
        3530             :         mono_cubic_interp  ! Procedure
        3531             :
        3532             :     implicit none
        3533             :
        3534             :     ! Input Variables
        3535             :     integer, intent(in) :: &
        3536             :       nz,    & ! Total number of vertical levels
        3537             :       k        ! Level to resolve variable value
        3538             :
        3539             :     real( kind = core_rknd ), dimension(nz), intent(in) :: &
        3540             :       z_vals,            & ! Height at each vertical level                  [m]
        3541             :       var_grid_value       ! Value of variable at each grid level           [units vary]
        3542             :
        3543             :     ! Output Variables
        3544             :     real( kind = core_rknd ) :: &
        3545             :       var_value            ! Value of variable at height level              [units vary]
        3546             :
        3547             :     ! Local Variables
        3548             :     integer :: km1, k00, kp1, kp2
        3549             :   !-----------------------------------------------------------------------
        3550             :
        3551             :     !----- Begin Code -----
        3552             :
        3553             :     if ( k >= 1 .and. k <= nz ) then
        3554             :       ! This is the simple case. No extrapolation necessary.
        3555             :       var_value = var_grid_value(k)
        3556             :     else if ( k == 0 ) then
        3557             :       ! Extrapolate below the lower boundary
        3558             :       km1 = nz
        3559             :       k00 = 1
        3560             :       kp1 = 2
        3561             :       kp2 = 3
        3562             :       var_value = mono_cubic_interp( z_vals(1)-(z_vals(2)-z_vals(1)), &
        3563             :                                      km1, k00, kp1, kp2, &
        3564             :                                      z_vals(km1), z_vals(k00), z_vals(kp1), z_vals(kp2), &
        3565             :                                      var_grid_value(km1), var_grid_value(k00), &
        3566             :                                      var_grid_value(kp1), var_grid_value(kp2) )
        3567             :     else if ( k == nz+1 ) then
        3568             :       ! Extrapolate above the upper boundary
        3569             :       km1 = nz
        3570             :       k00 = nz-1
        3571             :       kp1 = nz
        3572             :       kp2 = nz
        3573             :       var_value = mono_cubic_interp( z_vals(nz)+(z_vals(nz)-z_vals(nz-1)), &
        3574             :                                      km1, k00, kp1, kp2, &
        3575             :                                      z_vals(km1), z_vals(k00), z_vals(kp1), z_vals(kp2), &
        3576             :                                      var_grid_value(km1), var_grid_value(k00), &
        3577             :                                      var_grid_value(kp1), var_grid_value(kp2) )
        3578             :     else
        3579             :       ! Invalid height requested
        3580             :       var_value = -999._core_rknd
        3581             :     end if ! k > 1 .and. k < nz
        3582             :     return
        3583             :   end function var_value_integer_height
        3584             :
        3585             :   !=============================================================================
        3586             :   function var_subgrid_interp( nz, k, z_vals, var, z_interp, l_below ) result( var_value )
        3587             :
        3588             :   ! Description
        3589             :   !   Interpolates (or extrapolates) a variable to a value between grid
        3590             :   !   levels
        3591             :
        3592             :   ! References
        3593             :   !-----------------------------------------------------------------------
        3594             :
        3595             :     use clubb_precision, only: &
        3596             :         core_rknd       ! Constant
        3597             :
        3598             :     use interpolation, only: &
        3599             :         mono_cubic_interp   ! Procedure
        3600             :
        3601             :     implicit none
        3602             :
        3603             :     ! Input Variables
        3604             :     integer, intent(in) :: &
        3605             :       nz, &         ! Number of vertical levels
        3606             :       k             ! Grid level near interpolation target
        3607             :
        3608             :     real( kind = core_rknd ), dimension(nz), intent(in) :: &
        3609             :       z_vals, &     ! Height at each grid level          [m]
        3610             :       var           ! Variable values at grid levels     [units vary]
        3611             :
        3612             :     real( kind = core_rknd ), intent(in) :: &
        3613             :       z_interp      ! Interpolation target height        [m]
        3614             :
        3615             :     logical, intent(in) :: &
        3616             :       l_below       ! True if z_interp < z_vals(k), false otherwise
        3617             :
        3618             :     ! Output Variable
        3619             :     real( kind = core_rknd ) :: &
        3620             :       var_value     ! Interpolated value of variable     [units vary]
        3621             :
        3622             :     ! Local Variables
        3623             :     integer :: km1, k00, kp1, kp2 ! Parameters for call to mono_cubic_interp
        3624             :   !----------------------------------------------------------------------
        3625             :
        3626             :     !----- Begin Code -----
        3627             :     if ( l_below ) then
        3628             :
        3629             :       if ( k == 1 ) then ! Extrapolation
        3630             :         km1 = nz
        3631             :         k00 = 1
        3632             :         kp1 = 2
        3633             :         kp2 = 3
        3634             :       else if ( k == 2 ) then
        3635             :         km1 = 1
        3636             :         k00 = 1
        3637             :         kp1 = 2
        3638             :         kp2 = 3
        3639             :       else if ( k == nz ) then
        3640             :         km1 = nz-2
        3641             :         k00 = nz-1
        3642             :         kp1 = nz
        3643             :         kp2 = nz
        3644             :       else
        3645             :         km1 = k-2
        3646             :         k00 = k-1
        3647             :         kp1 = k
        3648             :         kp2 = k+1
        3649             :       end if ! k == 1
        3650             :
        3651             :     else ! .not. l_below
        3652             :
        3653             :       if ( k == 1 ) then
        3654             :         km1 = 1
        3655             :         k00 = 1
        3656             :         kp1 = 2
        3657             :         kp2 = 3
        3658             :       else if ( k == nz-1 ) then
        3659             :         km1 = nz-2
        3660             :         k00 = nz-1
        3661             :         kp1 = nz
        3662             :         kp2 = nz
        3663             :       else if ( k == nz ) then ! Extrapolation
        3664             :         km1 = nz
        3665             :         k00 = nz-1
        3666             :         kp1 = nz
        3667             :         kp2 = nz
        3668             :       else
        3669             :         km1 = k-1
        3670             :         k00 = k
        3671             :         kp1 = k+1
        3672             :         kp2 = k+2
        3673             :       end if ! k == 1
        3674             :
        3675             :     end if ! l_below
        3676             :
        3677             :     ! Now perform the interpolation
        3678             :     var_value = mono_cubic_interp( z_interp, km1, k00, kp1, kp2, &
        3679             :                                    z_vals(km1), z_vals(k00), z_vals(kp1), z_vals(kp2), &
        3680             :                                    var(km1), var(k00), var(kp1), var(kp2) )
        3681             :
        3682             :     return
        3683             :
        3684             :   end function var_subgrid_interp
        3685             :
        3686             :   !=============================================================================
        3687             :
        3688             : end module pdf_closure_module

Generated by: LCOV version 1.14