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