LCOV - code coverage report
Current view: | top level - physics/cam - subcol_SILHS.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 0 | 485 | 0.0 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 0 | 11 | 0.0 %

              Line data    Source code
           1             : module subcol_SILHS
           2             :    !---------------------------------------------------------------------------
           3             :    ! Purpose:
           4             :    !
           5             :    ! Implement a subcolumn scheme based on the Subgrid Importance Latin Hypercube
           6             :    ! Sampling (SILHS) functionality of the CLUBB moist turbulence parameterization.
           7             :    !
           8             :    !---------------------------------------------------------------------------
           9             :
          10             :    use shr_kind_mod,     only: r8=>shr_kind_r8, r4=>shr_kind_r4, i4=>shr_kind_i4
          11             :    use physics_types,    only: physics_state, physics_tend, physics_ptend
          12             :    use ppgrid,           only: pcols, psubcols, pver, pverp, begchunk, endchunk
          13             :    use constituents,     only: pcnst, cnst_get_ind
          14             :    use cam_abortutils,   only: endrun
          15             :    use cam_logfile,      only: iulog
          16             :    use cam_history,      only: addfld, add_default, outfld, horiz_only
          17             :    use ref_pres,         only: top_lev => trop_cloud_top_lev
          18             : #ifdef CLUBB_SGS
          19             : #ifdef SILHS
          20             :   use clubb_intr, only: &
          21             :         clubb_config_flags, &
          22             :         clubb_params_single_col, &
          23             :         stats_metadata, &
          24             :         stats_zt, stats_zm, stats_sfc, &
          25             :         pdf_params_chnk, &
          26             :         hm_metadata, &
          27             :         hydromet_dim, &
          28             :         pdf_dim
          29             :
          30             :    use clubb_api_module, only: &
          31             :         hmp2_ip_on_hmm2_ip_slope_type, &
          32             :         hmp2_ip_on_hmm2_ip_intrcpt_type, &
          33             :         precipitation_fractions, &
          34             :         stats, &
          35             :         core_rknd
          36             :
          37             :    use silhs_api_module, only: &
          38             :         silhs_config_flags_type
          39             : #endif
          40             : #endif
          41             :    use physconst,     only: cpair, gravit, latvap, latice, rair, rga, cappa
          42             :
          43             :    implicit none
          44             :    private
          45             :    save
          46             :
          47             :    public :: subcol_register_SILHS  !
          48             :    public :: subcol_init_SILHS      ! Initialize
          49             :    public :: subcol_gen_SILHS       ! Generate subcolumn fields by calling SILHS
          50             :    public :: subcol_readnl_SILHS    ! SILHS namelist reader
          51             :    public :: subcol_ptend_avg_SILHS
          52             :    public :: subcol_SILHS_var_covar_driver
          53             :    public :: subcol_SILHS_fill_holes_conserv
          54             :    public :: subcol_SILHS_hydromet_conc_tend_lim
          55             :    public :: init_state_subcol
          56             :    private :: fill_holes_sedimentation
          57             :    private :: fill_holes_same_phase_vert
          58             : #ifdef SILHS
          59             :    ! Calc subcol mean ! Calc subcol variance
          60             :    private :: meansc
          61             :    private :: stdsc
          62             :
          63             :    type (stats), target :: stats_lh_zt,   &
          64             :                            stats_lh_sfc
          65             :    !$omp threadprivate(stats_lh_zt, stats_lh_sfc)
          66             :
          67             :   real( kind = core_rknd ), dimension(:,:), allocatable :: &
          68             :       corr_array_n_cloud, &
          69             :       corr_array_n_below
          70             :
          71             : #endif
          72             :
          73             :    !-----
          74             :    ! Private module vars
          75             :    !-----
          76             :
          77             :    ! constituent indicies
          78             :    integer :: &
          79             :       ixq      = 0, &
          80             :       ixcldliq = 0, &
          81             :       ixnumliq = 0, &
          82             :       ixcldice = 0, &
          83             :       ixnumice = 0, &
          84             :       ixrain   = 0, &
          85             :       ixnumrain= 0, &
          86             :       ixsnow   = 0, &
          87             :       ixnumsnow= 0
          88             :
          89             :    ! Pbuf indicies
          90             :    integer :: thlm_idx, rtm_idx, ice_supersat_idx, &
          91             :               alst_idx, cld_idx, qrain_idx, qsnow_idx, &
          92             :               nrain_idx, nsnow_idx, ztodt_idx, tke_idx, kvh_idx, &
          93             :               prec_pcw_idx, snow_pcw_idx, prec_str_idx, snow_str_idx, &
          94             :               qcsedten_idx, qrsedten_idx, qisedten_idx, qssedten_idx, &
          95             :               vtrmc_idx, umr_idx, vtrmi_idx, ums_idx, qcsevap_idx, qisevap_idx
          96             :
          97             :    logical :: subcol_SILHS_weight  ! if set, sets up weights for averaging subcolumns for SILHS
          98             :    integer :: subcol_SILHS_numsubcol ! number of subcolumns for this run
          99             :    logical :: docldfracscaling = .false. ! Weight tendencies by cloud fraction
         100             :
         101             :    character(len=256) :: subcol_SILHS_corr_file_path
         102             :    character(len=16)  :: subcol_SILHS_corr_file_name
         103             :
         104             :    logical :: subcol_SILHS_q_to_micro, &
         105             :               subcol_SILHS_n_to_micro, &
         106             :               subcol_SILHS_use_clear_col, &
         107             :               subcol_SILHS_meanice, &
         108             :               subcol_SILHS_constrainmn
         109             :
         110             :    logical :: subcol_SILHS_var_covar_src
         111             :
         112             :    real(r8) :: subcol_SILHS_ncnp2_on_ncnm2
         113             :
         114             :    ! There may or may not be a better place to put this.
         115             :    real(r8), parameter :: p0_clubb = 100000._r8
         116             :
         117             :
         118             : !   real(r8) :: subcol_SILHS_c6rt, subcol_SILHS_c7, subcol_SILHS_c8, subcol_SILHS_c11, &
         119             : !               subcol_SILHS_c11b, subcol_SILHS_gamma_coef, &
         120             : !               subcol_SILHS_mult_coef, subcol_SILHS_mu
         121             :
         122             :    real(r8) :: ztodt  ! model timestep
         123             : #ifdef CLUBB_SGS
         124             : #ifdef SILHS
         125             :     type(hmp2_ip_on_hmm2_ip_slope_type) :: subcol_SILHS_hmp2_ip_on_hmm2_ip_slope
         126             :     type(hmp2_ip_on_hmm2_ip_intrcpt_type) :: subcol_SILHS_hmp2_ip_on_hmm2_ip_intrcpt
         127             :
         128             :     type(silhs_config_flags_type) :: silhs_config_flags
         129             : #endif
         130             : #endif
         131             :
         132             : contains
         133             :
         134           0 :    subroutine subcol_register_SILHS()
         135             :
         136             :       !--------------------------------
         137             :       ! Register fields needed by SILHS in the physics buffer
         138             :       ! Currently, most fields needed by SILHS but calculated by CLUBB are registered
         139             :       ! by clubb in clubb_intr.F90.
         140             :       !
         141             :       !--------------------------------
         142             : #ifdef CLUBB_SGS
         143             : #ifdef SILHS
         144             : #endif
         145             : #endif
         146           0 :    end subroutine subcol_register_SILHS
         147             :
         148             :
         149           0 :    subroutine subcol_readnl_SILHS(nlfile)
         150             : #ifdef CLUBB_SGS
         151             : #ifdef SILHS
         152             :       use namelist_utils,   only: find_group_name
         153             :       use spmd_utils,       only: masterproc, masterprocid, mpicom
         154             :       use spmd_utils,       only: mpi_integer, mpi_logical, mpi_character, mpir8, iam
         155             :       use clubb_api_module, only: core_rknd
         156             :       use silhs_api_module, only: set_default_silhs_config_flags_api, &
         157             :                                   initialize_silhs_config_flags_type_api, &
         158             :                                   print_silhs_config_flags_api
         159             : #endif
         160             : #endif
         161             :       character(len=*), intent(in) :: nlfile  ! filepath for file containing namelist input
         162             :
         163             :       ! Local variables
         164             :       integer :: unitn, ierr
         165             : #ifdef CLUBB_SGS
         166             : #ifdef SILHS
         167             :
         168             :       integer :: &
         169             :           cluster_allocation_strategy
         170             :
         171             :       logical :: &
         172             :           subcol_silhs_l_lh_importance_sampling, &
         173             :           subcol_silhs_l_Lscale_vert_avg, &
         174             :           subcol_silhs_l_lh_straight_mc, &
         175             :           subcol_silhs_l_lh_clustered_sampling, &
         176             :           subcol_silhs_l_rcm_in_cloud_k_lh_start, &
         177             :           subcol_silhs_l_random_k_lh_start, &
         178             :           subcol_silhs_l_max_overlap_in_cloud, &
         179             :           subcol_silhs_l_lh_instant_var_covar_src, &
         180             :           subcol_silhs_l_lh_limit_weights, &
         181             :           subcol_silhs_l_lh_var_frac, &
         182             :           subcol_silhs_l_lh_normalize_weights
         183             :
         184             :       namelist /subcol_SILHS_nl/ subcol_SILHS_weight, &
         185             :                                  subcol_SILHS_numsubcol, &
         186             :                                  subcol_SILHS_corr_file_path, &
         187             :                                  subcol_SILHS_corr_file_name, &
         188             :                                  subcol_SILHS_q_to_micro, &
         189             :                                  subcol_SILHS_n_to_micro, &
         190             :                                  subcol_SILHS_ncnp2_on_ncnm2, &
         191             :                                  subcol_SILHS_hmp2_ip_on_hmm2_ip_slope, &
         192             :                                  subcol_SILHS_hmp2_ip_on_hmm2_ip_intrcpt, &
         193             :                                  subcol_SILHS_meanice, &
         194             :                                  subcol_SILHS_use_clear_col, &
         195             :                                  subcol_SILHS_constrainmn, &
         196             :                                  subcol_SILHS_var_covar_src
         197             : !                                 subcol_SILHS_c6rt, subcol_SILHS_c7, &
         198             : !                                 subcol_SILHS_c8, subcol_SILHS_c11, subcol_SILHS_c11b, &
         199             : !                                 subcol_SILHS_gamma_coef, subcol_SILHS_mult_coef, subcol_SILHS_mu
         200             :
         201             :      namelist /silhs_config_flags_nl/ subcol_silhs_l_lh_importance_sampling, &
         202             :                                       subcol_silhs_l_Lscale_vert_avg, &
         203             :                                       subcol_silhs_l_lh_straight_mc, &
         204             :                                       subcol_silhs_l_lh_clustered_sampling, &
         205             :                                       subcol_silhs_l_rcm_in_cloud_k_lh_start, &
         206             :                                       subcol_silhs_l_random_k_lh_start, &
         207             :                                       subcol_silhs_l_max_overlap_in_cloud, &
         208             :                                       subcol_silhs_l_lh_instant_var_covar_src, &
         209             :                                       subcol_silhs_l_lh_limit_weights, &
         210             :                                       subcol_silhs_l_lh_var_frac, &
         211             :                                       subcol_silhs_l_lh_normalize_weights
         212             :
         213             :       !-----------------------------------------------------------------------------
         214             :       ! Set defaults
         215             :
         216             :       ! Eric Raut changed a default.
         217             :       subcol_SILHS_hmp2_ip_on_hmm2_ip_slope%Ni = 0.0_core_rknd
         218             :       subcol_SILHS_hmp2_ip_on_hmm2_ip_intrcpt%Ni = 0.5_core_rknd
         219             :
         220             :       if (masterproc) then
         221             :          open( newunit=unitn, file=trim(nlfile), status='old' )
         222             :          call find_group_name(unitn, 'subcol_SILHS_nl', status=ierr)
         223             :          if (ierr == 0) then
         224             :             read(unitn, subcol_SILHS_nl, iostat=ierr)
         225             :             if (ierr /= 0) then
         226             :                call endrun('subcol_readnl_SILHS: ERROR reading namelist')
         227             :             end if
         228             :          end if
         229             :          close(unitn)
         230             :       end if
         231             :
         232             :       ! Set default silhs_config_flags entires
         233             :       call set_default_silhs_config_flags_api( cluster_allocation_strategy, &
         234             :                                                subcol_silhs_l_lh_importance_sampling, &
         235             :                                                subcol_silhs_l_Lscale_vert_avg, &
         236             :                                                subcol_silhs_l_lh_straight_mc, &
         237             :                                                subcol_silhs_l_lh_clustered_sampling, &
         238             :                                                subcol_silhs_l_rcm_in_cloud_k_lh_start, &
         239             :                                                subcol_silhs_l_random_k_lh_start, &
         240             :                                                subcol_silhs_l_max_overlap_in_cloud, &
         241             :                                                subcol_silhs_l_lh_instant_var_covar_src, &
         242             :                                                subcol_silhs_l_lh_limit_weights, &
         243             :                                                subcol_silhs_l_lh_var_frac, &
         244             :                                                subcol_silhs_l_lh_normalize_weights )
         245             :
         246             :       ! Get silhs_config_flags entries from namelist
         247             :       if (masterproc) then
         248             :         open( newunit=unitn, file=trim(nlfile), status='old' )
         249             :         call find_group_name(unitn, 'silhs_config_flags_nl', status=ierr)
         250             :         if (ierr == 0) then
         251             :           read(unitn, silhs_config_flags_nl, iostat=ierr)
         252             :           if (ierr /= 0) then
         253             :             call endrun('silhs_config_flags_nl: ERROR reading namelist')
         254             :           end if
         255             :         end if
         256             :         close(unitn)
         257             :       end if
         258             :
         259             :       ! Save silhs_config_flags entries into module variable silhs_config_flags
         260             :       call initialize_silhs_config_flags_type_api( cluster_allocation_strategy, &
         261             :                                                    subcol_silhs_l_lh_importance_sampling, &
         262             :                                                    subcol_silhs_l_Lscale_vert_avg, &
         263             :                                                    subcol_silhs_l_lh_straight_mc, &
         264             :                                                    subcol_silhs_l_lh_clustered_sampling, &
         265             :                                                    subcol_silhs_l_rcm_in_cloud_k_lh_start, &
         266             :                                                    subcol_silhs_l_random_k_lh_start, &
         267             :                                                    subcol_silhs_l_max_overlap_in_cloud, &
         268             :                                                    subcol_silhs_l_lh_instant_var_covar_src, &
         269             :                                                    subcol_silhs_l_lh_limit_weights, &
         270             :                                                    subcol_silhs_l_lh_var_frac, &
         271             :                                                    subcol_silhs_l_lh_normalize_weights, &
         272             :                                                    silhs_config_flags )
         273             :
         274             :       ! Print the SILHS configurable flags
         275             :       call print_silhs_config_flags_api( iulog, silhs_config_flags ) ! Intent(in)
         276             :
         277             : #ifdef SPMD
         278             :       ! Broadcast namelist variables
         279             :       call mpi_bcast(subcol_SILHS_weight,    1, mpi_logical, masterprocid, mpicom, ierr)
         280             :       call mpi_bcast(subcol_SILHS_numsubcol, 1, mpi_integer, masterprocid, mpicom, ierr)
         281             :       call mpi_bcast(subcol_SILHS_corr_file_path, len(subcol_SILHS_corr_file_path), &
         282             :                      mpi_character, masterprocid, mpicom, ierr)
         283             :       call mpi_bcast(subcol_SILHS_corr_file_name, len(subcol_SILHS_corr_file_name), &
         284             :                      mpi_character, masterprocid, mpicom, ierr)
         285             :       call mpi_bcast(subcol_SILHS_use_clear_col, 1, mpi_logical, masterprocid, mpicom, ierr)
         286             :       call mpi_bcast(subcol_SILHS_constrainmn, 1, mpi_logical, masterprocid, mpicom, ierr)
         287             :       call mpi_bcast(subcol_SILHS_meanice, 1, mpi_logical, masterprocid, mpicom, ierr)
         288             :       call mpi_bcast(subcol_SILHS_q_to_micro, 1, mpi_logical, masterprocid, mpicom, ierr)
         289             :       call mpi_bcast(subcol_SILHS_n_to_micro, 1, mpi_logical, masterprocid, mpicom, ierr)
         290             :       call mpi_bcast(subcol_SILHS_var_covar_src,1,mpi_logical,masterprocid, mpicom, ierr)
         291             :       call mpi_bcast(subcol_SILHS_ncnp2_on_ncnm2, 1, mpir8, masterprocid, mpicom, ierr)
         292             :       call mpi_bcast(subcol_SILHS_hmp2_ip_on_hmm2_ip_slope%rr, 1, mpir8, masterprocid, mpicom, ierr)
         293             :       call mpi_bcast(subcol_SILHS_hmp2_ip_on_hmm2_ip_slope%Nr, 1, mpir8, masterprocid, mpicom, ierr)
         294             :       call mpi_bcast(subcol_SILHS_hmp2_ip_on_hmm2_ip_slope%ri, 1, mpir8, masterprocid, mpicom, ierr)
         295             :       call mpi_bcast(subcol_SILHS_hmp2_ip_on_hmm2_ip_slope%Ni, 1, mpir8, masterprocid, mpicom, ierr)
         296             :       call mpi_bcast(subcol_SILHS_hmp2_ip_on_hmm2_ip_slope%rs, 1, mpir8, masterprocid, mpicom, ierr)
         297             :       call mpi_bcast(subcol_SILHS_hmp2_ip_on_hmm2_ip_slope%Ns, 1, mpir8, masterprocid, mpicom, ierr)
         298             :       call mpi_bcast(subcol_SILHS_hmp2_ip_on_hmm2_ip_intrcpt%rr, 1, mpir8, masterprocid, mpicom, ierr)
         299             :       call mpi_bcast(subcol_SILHS_hmp2_ip_on_hmm2_ip_intrcpt%Nr, 1, mpir8, masterprocid, mpicom, ierr)
         300             :       call mpi_bcast(subcol_SILHS_hmp2_ip_on_hmm2_ip_intrcpt%ri, 1, mpir8, masterprocid, mpicom, ierr)
         301             :       call mpi_bcast(subcol_SILHS_hmp2_ip_on_hmm2_ip_intrcpt%Ni, 1, mpir8, masterprocid, mpicom, ierr)
         302             :       call mpi_bcast(subcol_SILHS_hmp2_ip_on_hmm2_ip_intrcpt%rs, 1, mpir8, masterprocid, mpicom, ierr)
         303             :       call mpi_bcast(subcol_SILHS_hmp2_ip_on_hmm2_ip_intrcpt%Ns, 1, mpir8, masterprocid, mpicom, ierr)
         304             : !      call mpi_bcast(subcol_SILHS_c6rt, 1, mpir8, masterprocid, mpicom, ierr)
         305             : !      call mpi_bcast(subcol_SILHS_c7, 1, mpir8, masterprocid, mpicom, ierr)
         306             : !      call mpi_bcast(subcol_SILHS_c8, 1, mpir8, masterprocid, mpicom, ierr)
         307             : !      call mpi_bcast(subcol_SILHS_c11, 1, mpir8, masterprocid, mpicom, ierr)
         308             : !      call mpi_bcast(subcol_SILHS_c11b, 1, mpir8, masterprocid, mpicom, ierr)
         309             : !      call mpi_bcast(subcol_SILHS_gamma_coef, 1, mpir8, masterprocid, mpicom, ierr)
         310             : !      call mpi_bcast(subcol_SILHS_mult_coef, 1, mpir8, masterprocid, mpicom, ierr)
         311             : !      call mpi_bcast(subcol_SILHS_mu, 1, mpir8, masterprocid, mpicom, ierr)
         312             :       call mpi_bcast(silhs_config_flags%l_lh_importance_sampling, 1, mpi_logical, masterprocid, mpicom, ierr)
         313             :       call mpi_bcast(silhs_config_flags%l_Lscale_vert_avg, 1, mpi_logical, masterprocid, mpicom, ierr)
         314             :       call mpi_bcast(silhs_config_flags%l_lh_straight_mc, 1, mpi_logical, masterprocid, mpicom, ierr)
         315             :       call mpi_bcast(silhs_config_flags%l_lh_clustered_sampling, 1, mpi_logical, masterprocid, mpicom, ierr)
         316             :       call mpi_bcast(silhs_config_flags%l_rcm_in_cloud_k_lh_start, 1, mpi_logical, masterprocid, mpicom, ierr)
         317             :       call mpi_bcast(silhs_config_flags%l_random_k_lh_start, 1, mpi_logical, masterprocid, mpicom, ierr)
         318             :       call mpi_bcast(silhs_config_flags%l_max_overlap_in_cloud, 1, mpi_logical, masterprocid, mpicom, ierr)
         319             :       call mpi_bcast(silhs_config_flags%l_lh_instant_var_covar_src, 1, mpi_logical, masterprocid, mpicom, ierr)
         320             :       call mpi_bcast(silhs_config_flags%l_lh_limit_weights, 1, mpi_logical, masterprocid, mpicom, ierr)
         321             :       call mpi_bcast(silhs_config_flags%l_lh_var_frac, 1, mpi_logical, masterprocid, mpicom, ierr)
         322             :       call mpi_bcast(silhs_config_flags%l_lh_normalize_weights, 1, mpi_logical, masterprocid, mpicom, ierr)
         323             :
         324             : ! SPMD
         325             : #endif
         326             : ! SILHS
         327             : #endif
         328             : ! CLUBB_SGS
         329             : #endif
         330           0 :    end subroutine subcol_readnl_SILHS
         331             :
         332             :
         333           0 :    subroutine subcol_init_SILHS(pbuf2d)
         334             :
         335             :       !--------------------------------
         336             :       ! Read in parameters and initialize SILHS PDF fields.
         337             :       ! Set up indexes into Pbuf fields.
         338             :       ! Register history outputs.
         339             :       !--------------------------------
         340             :
         341             :       use physics_buffer,          only: physics_buffer_desc, pbuf_get_field, &
         342             :                                          dtype_r8, pbuf_get_index
         343             : #ifdef CLUBB_SGS
         344             : #ifdef SILHS
         345             :       use clubb_api_module,        only: core_rknd, &
         346             :                                          setup_corr_varnce_array_api, &
         347             :                                          init_pdf_hydromet_arrays_api, &
         348             :                                          set_clubb_debug_level_api
         349             :
         350             : #endif
         351             : #endif
         352             :
         353             :       type(physics_buffer_desc), pointer :: pbuf2d(:,:)
         354             :
         355             : #ifdef CLUBB_SGS
         356             : #ifdef SILHS
         357             :
         358             :       integer :: iunit = 501 ! Default value, will get iunit from CAM
         359             :       !character(len=*), parameter :: default_corr_case = "arm_97"
         360             :       character(len=*), parameter :: &
         361             :             cloud_file_ext  = "_corr_array_cloud.in", & ! File extensions for corr files
         362             :             below_file_ext  = "_corr_array_below.in"
         363             :       character(len=256) :: corr_file_path_cloud, corr_file_path_below
         364             :
         365             :       ! To set up CLUBB hydromet indices
         366             :       integer :: &
         367             :           iirr,         & ! Hydrometeor array index for rain water mixing ratio, rr
         368             :           iirs,         & ! Hydrometeor array index for snow mixing ratio, rs
         369             :           iiri,         & ! Hydrometeor array index for ice mixing ratio, ri
         370             :           iirg,         & ! Hydrometeor array index for graupel mixing ratio, rg
         371             :           iiNr,         & ! Hydrometeor array index for rain drop concentration, Nr
         372             :           iiNs,         & ! Hydrometeor array index for snow concentration, Ns
         373             :           iiNi,         & ! Hydrometeor array index for ice concentration, Ni
         374             :           iiNg            ! Hydrometeor array index for graupel concentration, Ng
         375             :
         376             :       integer :: l, ierr=0  ! Loop variable, error check
         377             :
         378             :       ! Set CLUBB's debug level
         379             :       ! This is called in module clubb_intr; no need to do it here.
         380             : !      call set_clubb_debug_level_api( 0 )
         381             :
         382             :       !-------------------------------
         383             :       ! CLUBB-SILHS Parameters (global module variables)
         384             :       !-------------------------------
         385             :       clubb_config_flags%l_fix_w_chi_eta_correlations = .true.
         386             :       clubb_config_flags%l_diagnose_correlations = .false.
         387             :       clubb_config_flags%l_calc_w_corr = .false.
         388             : !      l_prescribed_avg_deltaz = .false.
         389             :       clubb_config_flags%l_use_cloud_cover = .false.
         390             :       clubb_config_flags%l_const_Nc_in_cloud = .true.
         391             :
         392             :       ! Values from the namelist
         393             :       docldfracscaling = subcol_SILHS_use_clear_col
         394             :
         395             :       ! Namelist "tuning" or set correlations
         396             :       ! KTC Todo: Move these to a tuning "in" file or into the namelist
         397             :       ! JHTODO: we might want these on CLUBB's API and ultimatively on a namelist for tuning
         398             : !      C6rt = subcol_SILHS_c6rt
         399             : !      C7 = subcol_SILHS_c7                                      ! to all ice clouds
         400             : !      C8 = subcol_SILHS_c8
         401             : !      C11 = subcol_SILHS_c11
         402             : !      C11b = subcol_SILHS_c11b
         403             : !      gamma_coef = subcol_SILHS_gamma_coef
         404             : !      mult_coef = subcol_SILHS_mult_coef
         405             : !      mu = subcol_SILHS_mu
         406             :
         407             :       !call set_clubb_debug_level( 0 )  !#KTCtodo: Add a namelist variable to set debug level
         408             :
         409             :       ! Get constituent indices
         410             :       call cnst_get_ind('Q', ixq)
         411             :       call cnst_get_ind('CLDLIQ', ixcldliq)
         412             :       call cnst_get_ind('NUMLIQ', ixnumliq)
         413             :       call cnst_get_ind('CLDICE', ixcldice)
         414             :       call cnst_get_ind('NUMICE', ixnumice)
         415             :       call cnst_get_ind('RAINQM', ixrain, abort=.false.)
         416             :       call cnst_get_ind('NUMRAI', ixnumrain, abort=.false.)
         417             :       call cnst_get_ind('SNOWQM', ixsnow, abort=.false.)
         418             :       call cnst_get_ind('NUMSNO', ixnumsnow, abort=.false.)
         419             :
         420             :       ! Get physics buffer indexes
         421             :       thlm_idx = pbuf_get_index('THLM')
         422             :       rtm_idx = pbuf_get_index('RTM')
         423             :       cld_idx = pbuf_get_index('CLD')
         424             :       alst_idx = pbuf_get_index('ALST')  ! SILHS expects clubb's cloud_frac liq stratus fraction
         425             :       ztodt_idx = pbuf_get_index('ZTODT')
         426             :       ice_supersat_idx = pbuf_get_index('ISS_FRAC')
         427             :       tke_idx = pbuf_get_index('tke')
         428             :       kvh_idx = pbuf_get_index('kvh')
         429             :       prec_pcw_idx = pbuf_get_index('PREC_PCW')
         430             :       snow_pcw_idx = pbuf_get_index('SNOW_PCW')
         431             :       prec_str_idx = pbuf_get_index('PREC_STR')
         432             :       snow_str_idx = pbuf_get_index('SNOW_STR')
         433             :       qcsedten_idx = pbuf_get_index('QCSEDTEN')
         434             :       qrsedten_idx = pbuf_get_index('QRSEDTEN')
         435             :       qisedten_idx = pbuf_get_index('QISEDTEN')
         436             :       qssedten_idx = pbuf_get_index('QSSEDTEN')
         437             :       vtrmc_idx = pbuf_get_index('VTRMC')
         438             :       umr_idx = pbuf_get_index('UMR')
         439             :       vtrmi_idx = pbuf_get_index('VTRMI')
         440             :       ums_idx = pbuf_get_index('UMS')
         441             :       qcsevap_idx = pbuf_get_index('QCSEVAP')
         442             :       qisevap_idx = pbuf_get_index('QISEVAP')
         443             :       qrain_idx = pbuf_get_index('QRAIN')
         444             :       qsnow_idx = pbuf_get_index('QSNOW')
         445             :       nrain_idx = pbuf_get_index('NRAIN')
         446             :       nsnow_idx = pbuf_get_index('NSNOW')
         447             :
         448             :       !-------------------------------
         449             :       ! Set up SILHS hydrometeors #KTCtodo: move microphys specification to config time,
         450             :       !        Steve wants to set up a microphysics query so I can ask the microphysics
         451             :       !        scheme which hydrometeors to use. For the future.
         452             :       !-------------------------------
         453             :       iirr = 1
         454             :       iirs = 3
         455             :       iiri = 5
         456             :       iirg = -1
         457             :
         458             :       iiNr = 2
         459             :       iiNs = 4
         460             :       iiNi = 6
         461             :       iiNg = -1
         462             :
         463             :       hydromet_dim = 6
         464             :
         465             :       ! Set up pdf indices, hydromet indicies, hydromet arrays, and hydromet variance ratios
         466             :       call init_pdf_hydromet_arrays_api( 1.0_core_rknd, 1.0_core_rknd, hydromet_dim,  & ! intent(in)
         467             :                                          iirr, iiNr, iiri, iiNi,                      & ! intent(in)
         468             :                                          iirs, iiNs, iirg, iiNg,                      & ! intent(in)
         469             :                                          subcol_SILHS_ncnp2_on_ncnm2,                 & ! intent(in)
         470             :                                          hm_metadata, pdf_dim,                        & ! intent(out)
         471             :                                          subcol_SILHS_hmp2_ip_on_hmm2_ip_slope,       & ! optional(in)
         472             :                                          subcol_SILHS_hmp2_ip_on_hmm2_ip_intrcpt )      ! optional(in)
         473             :
         474             :       !-------------------------------
         475             :       ! Set up hydrometeors and correlation arrays for SILHS
         476             :       !-------------------------------
         477             :       allocate( corr_array_n_cloud(pdf_dim,pdf_dim), corr_array_n_below(pdf_dim,pdf_dim), stat=ierr)
         478             :       if( ierr /= 0 ) call endrun(' subcol_init_SILHS: failed to allocate corr_array fields ')
         479             :
         480             :       corr_file_path_cloud = trim( subcol_SILHS_corr_file_path )//trim( subcol_SILHS_corr_file_name )//cloud_file_ext
         481             :       corr_file_path_below = trim( subcol_SILHS_corr_file_path )//trim( subcol_SILHS_corr_file_name )//below_file_ext
         482             :
         483             :       call setup_corr_varnce_array_api( corr_file_path_cloud, corr_file_path_below, &
         484             :                                         pdf_dim, hm_metadata, newunit(iunit), &
         485             :                                         clubb_config_flags%l_fix_w_chi_eta_correlations,               & ! In
         486             :                                         corr_array_n_cloud, corr_array_n_below )
         487             :
         488             :       !-------------------------------
         489             :       ! Register output fields from SILHS
         490             :       !-------------------------------
         491             :       call addfld('SILHS_NCLD_SCOL', (/'psubcols', 'ilev    '/), 'I', 'm^-3', &
         492             :            'Subcolumn Cloud Number Concentration', flag_xyfill=.true., fill_value=1.e30_r8)
         493             :       call addfld('SILHS_NRAIN_SCOL', (/'psubcols', 'ilev    '/), 'I', 'm^-3', &
         494             :            'Subcolumn Number Concentration of Rain from SILHS', flag_xyfill=.true., fill_value=1.e30_r8)
         495             :       call addfld('SILHS_OMEGA_SCOL', (/'psubcols', 'ilev    '/), 'I', 'Pa/s', &
         496             :            'Subcolumn vertical pressure velocity', flag_xyfill=.true., fill_value=1.e30_r8)
         497             :       call addfld('SILHS_RCM_SCOL', (/'psubcols', 'ilev    '/), 'I', 'kg/kg', &
         498             :            'Subcolumn Cloud Liquid Water from SILHS', flag_xyfill=.true., fill_value=1.e30_r8)
         499             :       call addfld('SILHS_RICLD_SCOL', (/'psubcols', 'ilev    '/), 'I', 'kg/kg', &
         500             :            'Subcolumn Cloud Ice Water from SILHS', flag_xyfill=.true., fill_value=1.e30_r8)
         501             :       call addfld('SILHS_NICLD_SCOL', (/'psubcols', 'ilev    '/), 'I', 'kg/kg', &
         502             :            'Subcolumn Cloud Ice Number Conc from SILHS', flag_xyfill=.true., fill_value=1.e30_r8)
         503             :       call addfld('SILHS_RRAIN_SCOL', (/'psubcols', 'ilev    '/), 'I', 'kg/kg', &
         504             :            'Subcolumn Precipitating Liquid Water from SILHS', flag_xyfill=.true., fill_value=1.e30_r8)
         505             :       call addfld('SILHS_RT_SCOL', (/'psubcols', 'ilev    '/), 'I', 'kg/kg ', &
         506             :            'Subcolumn Total Water from SILHS', flag_xyfill=.true., fill_value=1.e30_r8)
         507             :       call addfld('SILHS_THLM_SCOL', (/'psubcols', 'ilev    '/), 'I', 'K', &
         508             :            'Subcolumn liquid water pot temperature', flag_xyfill=.true., fill_value=1.e30_r8)
         509             :       call addfld('SILHS_WEIGHT_SCOL', (/'psubcols'/), 'I', 'frac', &
         510             :            'Weights for each subcolumn', flag_xyfill=.true., fill_value=1.e30_r8)
         511             :       call addfld('SILHS_WM_SCOL', (/'psubcols', 'ilev    '/), 'I', 'm/s', &
         512             :            'Subcolumn vertical velocity from SILHS', flag_xyfill=.true., fill_value=1.e30_r8)
         513             :
         514             :       call addfld('NR_IN_LH', (/ 'lev' /), 'I', 'm^-3', &
         515             :                   'Num Rain Conc as input to SILHS')
         516             :      call addfld('SILHS_RTM', (/ 'ilev' /), 'I', 'kg/kg', &
         517             :                   'Input total water mixing ratio')
         518             :      call addfld('SILHS_THLM', (/ 'ilev' /), 'I', 'K', &
         519             :                   'Input liquid water potential temperature')
         520             :      call addfld('SILHS_QC_IN', (/ 'lev' /), 'I', 'kg/kg', &
         521             :                   'Input cloud water mixing ratio')
         522             :      call addfld('SILHS_QI_IN', (/ 'lev' /), 'I', 'kg/kg', &
         523             :                   'Input cloud ice mixing ratio')
         524             :      call addfld('SILHS_NC_IN', (/ 'lev' /), 'I', '#/kg', &
         525             :                   'Input cloud water number concentration')
         526             :      call addfld('SILHS_NI_IN', (/ 'lev' /), 'I', '#/kg', &
         527             :                   'Input cloud ice number concentration')
         528             :      call addfld('AKM_CLUBB', (/ 'ilev' /), 'I', '(kg/kg)/s', &
         529             :                   'Exact Kessler autoconversion')
         530             :      call addfld('AKM_LH_CLUBB', (/ 'ilev' /), 'I', '(kg/kg)/s', &
         531             :                   'Monte Carlo estimate of Kessler autoconversion')
         532             :      call addfld('INVS_EXNER', (/ 'lev' /), 'I', 'none', &
         533             :                   'inverse EXNER function from state in subcol_SILHS')
         534             :      call addfld('SILHS_ZTODT', horiz_only, 'I', 's', &
         535             :                   'Length of Physics timestep (for debugging)')
         536             :      if ( subcol_SILHS_constrainmn ) then
         537             :         call addfld('SILHS_MSC_CLDICE', (/ 'lev' /), 'A', 'kg/kg', &
         538             :                     'Mean Cloud Ice across subcolumns')
         539             :         call addfld('SILHS_STDSC_CLDICE', (/ 'lev' /), 'A', 'kg/kg', &
         540             :                     'Standard deviation of Ice across subcolumns')
         541             :         if ( ixsnow > 0 ) then
         542             :            call addfld('SILHS_MSC_CLDLIQ', (/ 'lev' /), 'A', 'kg/kg', &
         543             :                        'Mean Cloud Liquid across subcolumns')
         544             :            call addfld('SILHS_STDSC_CLDLIQ', (/ 'lev' /), 'A', 'kg/kg', &
         545             :                        'Standard deviation of Liquid across subcolumns')
         546             :            call addfld('SILHS_MSC_Q', (/ 'lev' /), 'A', 'kg/kg', &
         547             :                        'Mean water vapor across subcolumns')
         548             :            call addfld('SILHS_STDSC_Q', (/ 'lev' /), 'A', 'kg/kg', &
         549             :                        'Standard deviation of water vapor across subcolumns')
         550             :         endif ! ixsnow > 0
         551             :      endif ! subcol_SILHS_constrainmn
         552             :      call addfld('SILHS_EFF_CLDFRAC', (/ 'lev' /), 'A', 'frac', &
         553             :                   'Calculated cloud fraction from subcolumn liq or ice')
         554             :
         555             :      call addfld('SILHS_CLUBB_PRECIP_FRAC', (/ 'lev' /), 'A', 'frac', &
         556             :                   'Precipitation fraction from CLUBB (set_up_pdf_params_incl_hydromet)')
         557             :      call addfld('SILHS_CLUBB_ICE_SS_FRAC', (/ 'lev' /), 'A', 'frac', &
         558             :                   'Ice supersaturation fraction from CLUBB')
         559             :
         560             :      call addfld ('QVHFTEN', (/ 'lev' /), 'A', 'kg/kg/s', 'Water vapor mixing ratio tendency from hole filling')
         561             :      call addfld ('QCHFTEN', (/ 'lev' /), 'A', 'kg/kg/s', 'Cloud water mixing ratio tendency from hole filling')
         562             :      call addfld ('QRHFTEN', (/ 'lev' /), 'A', 'kg/kg/s', 'Rain water mixing ratio tendency from hole filling')
         563             :      call addfld ('QIHFTEN', (/ 'lev' /), 'A', 'kg/kg/s', 'Cloud ice mixing ratio tendency from hole filling')
         564             :      call addfld ('QSHFTEN', (/ 'lev' /), 'A', 'kg/kg/s', 'Snow mixing ratio tendency from hole filling')
         565             :      call addfld ('THFTEN', (/ 'lev' /), 'A', 'K/s', 'Temperature tendency from hole filling')
         566             :
         567             : #endif
         568             : #endif
         569           0 :    end subroutine subcol_init_SILHS
         570             : !==============================================================!
         571           0 :    subroutine init_state_subcol(state, tend, state_sc, tend_sc)
         572             :
         573           0 :      use ppgrid,                 only : pver, pverp, pcols
         574             :
         575             :      use subcol_utils,           only : subcol_set_subcols
         576             :
         577             :      implicit none
         578             :
         579             :      type(physics_state), intent(inout) :: state
         580             :      type(physics_tend),  intent(inout) :: tend
         581             :      type(physics_state), intent(inout) :: state_sc        ! sub-column state
         582             :      type(physics_tend),  intent(inout) :: tend_sc         ! sub-column tend
         583             :
         584             :      integer, dimension(pcols) :: numsubcol_arr ! To set up the state struct
         585             :
         586           0 :      numsubcol_arr(:) = 0  ! Start over each chunk
         587           0 :      numsubcol_arr(:state%ngrdcol) = subcol_SILHS_numsubcol ! Only set for valid grid columns
         588           0 :      call subcol_set_subcols(state, tend, numsubcol_arr, state_sc, tend_sc)
         589             :
         590           0 :    end subroutine init_state_subcol
         591             : !==================================================================!
         592           0 :    subroutine subcol_gen_SILHS(state, tend, state_sc, tend_sc, pbuf)
         593             :       !-------------------------------
         594             :       ! This is where the subcolumns are created, and the call to
         595             :       !      generate_silhs_sample_mod_api
         596             :       !    goes out. Variables needed to make this call are pulled from the
         597             :       !    pbuf, from module data, and calculated based on the CAM state.
         598             :       !-------------------------------
         599             :
         600           0 :       use physics_buffer,         only : physics_buffer_desc, pbuf_get_index, &
         601             :                                          pbuf_get_field
         602             :       use time_manager,           only : get_nstep
         603             :       use subcol_utils,           only : subcol_set_subcols, subcol_set_weight
         604             :       use subcol_pack_mod,        only : subcol_pack
         605             :       use phys_control,           only : phys_getopts
         606             :       use spmd_utils,             only : masterproc
         607             :       use shr_const_mod,          only : SHR_CONST_PI, SHR_CONST_RHOFW
         608             :
         609             : #ifdef CLUBB_SGS
         610             : #ifdef SILHS
         611             :       use clubb_api_module,       only : setup_pdf_parameters_api, &
         612             :
         613             :                                          zm2zt_api, setup_grid_heights_api, &
         614             :
         615             :                                          core_rknd, &
         616             :
         617             :                                          w_tol_sqd, zero_threshold, &
         618             :                                          em_min, cloud_frac_min, & ! rc_tol, &
         619             :
         620             :                                          genrand_intg, genrand_init_api, &
         621             :
         622             :                                          nparams, ic_K, &
         623             :                                          read_parameters_api, &
         624             :                                          Cp, Lv, &
         625             :                                          grid, setup_grid_api, &
         626             :                                          init_precip_fracs_api
         627             :
         628             :       use silhs_api_module, only :       generate_silhs_sample_api, & ! Ncn_to_Nc, &
         629             :                                          clip_transform_silhs_output_api, &
         630             :                                          est_kessler_microphys_api, &
         631             :                                          vert_decorr_coef
         632             :
         633             : #endif
         634             : #endif
         635             :
         636             :       ! CAM data structures
         637             :       type(physics_state), intent(inout) :: state
         638             :       type(physics_tend),  intent(inout) :: tend
         639             :       type(physics_state), intent(inout) :: state_sc        ! sub-column state
         640             :       type(physics_tend),  intent(inout) :: tend_sc         ! sub-column tend
         641             :       type(physics_buffer_desc), pointer :: pbuf(:)
         642             :
         643             : #ifdef CLUBB_SGS
         644             : #ifdef SILHS
         645             :       !----------------
         646             :       ! Local variables
         647             :       !----------------
         648             :       logical, parameter :: &
         649             :                  l_implemented = .true.   ! Implemented in a host model
         650             :       logical, parameter :: rx_Nc = .false. ! Use NC calculated based on grid mean effective radius
         651             :       integer, parameter :: &
         652             :                  grid_type = 3            ! The 3 grid centered on momentum levels
         653             :       real(r8), parameter :: cldmin = 0.001_r8 ! To use when cld frac = 0.0 to be consistant with micro_mg
         654             :       real(r8), parameter :: min_num_conc = 1.0e-12_r8
         655             :       real(r8), parameter :: qsmall = 1.0e-18_r8  ! Microphysics cut-off for cloud
         656             :
         657             :       integer :: i, j, k, ngrdcol, ncol, lchnk, stncol
         658             :       real(r8) :: sfc_elevation(state%ngrdcol)  ! Surface elevation
         659             :
         660             :       real(r8), dimension(state%ngrdcol,pverp-top_lev+1) :: zt_g, zi_g ! Thermo & Momentum grids for clubb
         661             :
         662             :       real(r8), dimension(pverp) :: scfrac     ! cloud fraction based on sc distributions
         663             :       real(r8) :: msc, std, maxcldfrac, maxsccldfrac
         664             :       real(r8) :: scale = 1.0_r8
         665             :
         666             :       real(r8) :: c_K ! CLUBB parameter c_K (for eddy diffusivity)
         667             :
         668             :       integer( kind = genrand_intg ) :: &
         669             :         lh_seed    ! Seed used in random number generator that will be different
         670             :                    ! for each column, yet reproducible for a restart run
         671             :
         672             :       !----------------
         673             :       ! Required for set_up_pdf_params_incl_hydromet
         674             :       !----------------
         675             :       real(r8), dimension(state%ngrdcol,pverp-top_lev+1) :: cld_frac_in  ! Cloud fraction
         676             :
         677             :       real(r8), dimension(state%ngrdcol, pverp-top_lev+1, pdf_dim, pdf_dim) :: &
         678             :                           corr_array_1, corr_array_2  ! Correlation matrix for pdf components
         679             :
         680             :       real(r8), dimension(state%ngrdcol, pverp-top_lev+1, pdf_dim) :: &
         681             :                           mu_x_1, mu_x_2, &    ! Mean array for PDF components
         682             :                           sigma_x_1, sigma_x_2 ! Std dev arr for PDF components
         683             :
         684             :       real(r8), dimension(state%ngrdcol, pverp-top_lev+1, pdf_dim, pdf_dim) :: &
         685             :                           corr_cholesky_mtx_1, corr_cholesky_mtx_2  ! Transposed corr cholesky mtx
         686             :
         687             :       real(r8), dimension(state%ngrdcol, pverp-top_lev+1) :: Nc_in_cloud
         688             :       real(r8), dimension(state%ngrdcol, pverp-top_lev+1) :: ice_supersat_frac_in
         689             :       real(r8), dimension(state%ngrdcol, pverp-top_lev+1, hydromet_dim) :: hydrometp2
         690             :
         691             :
         692             :       !----------------
         693             :       ! Input to generate_silhs_sample
         694             :       !----------------
         695             :       integer :: iter                            ! CLUBB iteration
         696             :       integer :: num_subcols                     ! Number of subcolumns
         697             :       integer, parameter :: sequence_length = 1  ! Number of timesteps btn subcol calls
         698             :
         699             :       real(r8), dimension(state%ngrdcol,pverp-top_lev+1) :: rho_ds_zt    ! Dry static density (kg/m^3) on thermo levs
         700             :       real(r8), dimension(state%ngrdcol,pver)  :: dz_g         ! thickness of layer
         701             :       real(r8), dimension(state%ngrdcol,pverp-top_lev+1) :: delta_zm     ! Difference in u wind altitudes
         702             :
         703             :       real(r8), dimension(state%ngrdcol,pverp-top_lev+1,hydromet_dim) :: hydromet  ! Hydrometeor species
         704             :       real(r8), dimension(state%ngrdcol,pverp-top_lev+1,hydromet_dim) :: wphydrometp  ! Hydrometeor flux
         705             :       real(r8), dimension(state%ngrdcol,pverp-top_lev+1)              :: Ncm ! Mean cloud droplet concentration, <N_c>
         706             :
         707             :       real(r8), dimension(state%ngrdcol,pverp-top_lev+1) :: tke       ! TKE
         708             :       real(r8), dimension(state%ngrdcol,pverp-top_lev+1) :: khzm      ! Eddy diffusivity coef
         709             :       real(r8), dimension(state%ngrdcol,pverp-top_lev+1) :: Lscale_zm ! CLUBB's length scale on momentum (zm) levels
         710             :       real(r8), dimension(state%ngrdcol,pverp-top_lev+1) :: Lscale    ! CLUBB's length scale
         711             :
         712             :       logical, parameter :: &
         713             :          l_calc_weights_all_levs = .false. ! .false. if all time steps use the same
         714             :                                           !   weights at all vertical grid levels
         715             :       logical :: &
         716             :         l_calc_weights_all_levs_itime, & ! .true. if we calculate sample weights separately at all
         717             :                                          !    grid levels at the current time step
         718             :         l_rad_itime                      ! .true. if we calculate radiation at the current time step
         719             :
         720             :       !---------------
         721             :       !Output from generate_silhs_sample
         722             :       !--------------
         723             :       real(r8), dimension(state%ngrdcol,subcol_SILHS_numsubcol,pverp-top_lev+1,pdf_dim) :: X_nl_all_levs ! Sample transformed to normal-lognormal
         724             :       real(r8), dimension(state%ngrdcol,subcol_SILHS_numsubcol,pverp-top_lev+1)   :: lh_sample_point_weights ! Subcolumn weights
         725             :       integer, dimension(state%ngrdcol,subcol_SILHS_numsubcol,pverp-top_lev+1)    :: X_mixt_comp_all_levs ! Which Mixture Component
         726             :
         727             :       real(r8), dimension(state%ngrdcol,pverp-top_lev+1, subcol_SILHS_numsubcol) :: &
         728             :                   rc_all_points, & ! Calculate RCM from LH output
         729             :                   rain_all_pts,  & ! Calculate Rain from LH output
         730             :                   nrain_all_pts, & ! Calculate Rain Conc from LH
         731             :                   snow_all_pts,  & ! Calculate Snow from LH output
         732             :                   nsnow_all_pts, & ! Calculate Snow Conc from LH
         733             :                   w_all_points,  & ! Calculate W from LH output
         734             :                   ice_all_pts,   & ! Calculate Cld Ice from LH output
         735             :                   nice_all_pts,  & ! Calculate Num cld ice from LH
         736             :                   nclw_all_pts     ! Calculate Num cld wat from LH
         737             :
         738             :       !----------------
         739             :       ! Output from clip_transform_silhs_output_api
         740             :       !----------------
         741             :       real( kind = core_rknd ), dimension(state%ngrdcol,subcol_SILHS_numsubcol,pverp-top_lev+1) :: &
         742             :         lh_rt_clipped,  & ! rt generated from silhs sample points
         743             :         lh_thl_clipped, & ! thl generated from silhs sample points
         744             :         lh_rc_clipped,  & ! rc generated from silhs sample points
         745             :         lh_rv_clipped,  & ! rv generated from silhs sample points
         746             :         lh_Nc_clipped     ! Nc generated from silhs sample points
         747             :
         748             :       logical, parameter :: &
         749             :         l_use_Ncn_to_Nc = .true.  ! Whether to call Ncn_to_Nc (.true.) or not (.false.);
         750             :                                   ! Ncn_to_Nc might cause problems with the MG microphysics
         751             :                                   ! since the changes made here (Nc-tendency) are not fed into
         752             :                                   ! the microphysics
         753             :
         754             :
         755             :       !----------------
         756             :       ! Output to history
         757             :       !----------------
         758             :       ! V. Larson note: These variables are on the zt (full) levels: why do they
         759             :       ! have dimension pverp?  The pverp level corresponds to the CLUBB
         760             :       ! below-ground level.
         761             :       ! The variables in this paragraph are oriented like CAM variables (k=1 is
         762             :       ! the model top).
         763             :       ! They are flipped versions of CLUBB variables, for the entire chunk.
         764             :       real(r8), dimension(pcols*psubcols, pverp) :: RT_lh_out
         765             :       real(r8), dimension(pcols*psubcols, pverp) :: THL_lh_out
         766             :       real(r8), dimension(pcols*psubcols, pverp) :: OMEGA_lh_out
         767             :       real(r8), dimension(pcols*psubcols, pverp) :: WM_lh_out
         768             :       real(r8), dimension(pcols*psubcols, pverp) :: RVM_lh_out
         769             :       real(r8), dimension(pcols*psubcols, pverp) :: RCM_lh_out
         770             :       real(r8), dimension(pcols*psubcols, pverp) :: NCLW_lh_out
         771             :       real(r8), dimension(pcols*psubcols, pverp) :: ICE_lh_out
         772             :       real(r8), dimension(pcols*psubcols, pverp) :: NICE_lh_out
         773             :       real(r8), dimension(pcols*psubcols, pverp) :: RAIN_lh_out
         774             :       real(r8), dimension(pcols*psubcols, pverp) :: NRAIN_lh_out
         775             :       real(r8), dimension(pcols*psubcols, pverp) :: SNOW_lh_out
         776             :       real(r8), dimension(pcols*psubcols, pverp) :: NSNOW_lh_out
         777             :
         778             :       real(r8), dimension(state_sc%psetcols) :: weights ! Subcol weights
         779             :
         780             :       real(r8), dimension(pcols, pver) :: meansc_ice
         781             :       real(r8), dimension(pcols, pver) :: stdsc_ice
         782             :
         783             :       real(r8), dimension(pcols, pver) :: meansc_liq
         784             :       real(r8), dimension(pcols, pver) :: stdsc_liq
         785             :
         786             :       real(r8), dimension(pcols, pver) :: meansc_vap
         787             :       real(r8), dimension(pcols, pver) :: stdsc_vap
         788             :       real(r8), dimension(pcols, pver) :: grmn_eff_rad
         789             :       real(r8), dimension(pcols, pver) :: eff_cldfrac
         790             :       real(r8), dimension(pcols, pver) :: precip_frac_out
         791             :
         792             :       real(r8) :: tmp_mean, diff_mean, rcubed
         793             :
         794             :       !----------------
         795             :       ! Output from Est_Kessler_microphys
         796             :       !----------------
         797             :       real(r8), dimension(state%ngrdcol,pverp-top_lev+1) :: lh_Akm     ! Monte Carlo estimate of Kessler Autoconversion
         798             :       real(r8), dimension(state%ngrdcol,pverp-top_lev+1) :: AKm        ! Exact Kessler autoconversion
         799             :       real(r8), dimension(state%ngrdcol,pverp-top_lev+1) :: AKstd      ! Exact Stdev of gba Kessler
         800             :       real(r8), dimension(state%ngrdcol,pverp-top_lev+1) :: AKstd_cld  ! Exact w/in cloud stdev of gba Kessler
         801             :       real(r8), dimension(state%ngrdcol,pverp-top_lev+1) :: AKm_rcm    ! Exact local gba Kessler auto based on rcm
         802             :       real(r8), dimension(state%ngrdcol,pverp-top_lev+1) :: AKm_rcc    ! Exact local gba Kessler based on w/in cloud rc
         803             :       real(r8), dimension(state%ngrdcol,pverp-top_lev+1) :: lh_rcm_avg ! LH estimate of grid box avg liquid water
         804             :       real(r8), dimension(pcols,pverp) :: lh_AKm_out, AKm_out
         805             :
         806             :       !----------------
         807             :       ! Needed to update State
         808             :       !----------------
         809             :       real(r8), dimension(pver)  :: Temp_prof  ! Subcolumn LWPT converted to Abs Temp
         810             :       real(r8), dimension(pver)  :: SE_prof    ! Static Energy calculated from Abs Temp
         811             :       real(r8), dimension(pver)  :: No_cloud = 0.0_r8     ! Clear air condensate profile
         812             :       real(r8), dimension(pcols, pver)  :: invs_exner  ! inverse exner sent to conversion codw
         813             :                                                        ! pcols for output to history
         814             :       real(r8) :: eff_rad_coef = 1.0_r8/(4.0_r8/3.0_r8*SHR_CONST_RHOFW*SHR_CONST_PI)
         815             :
         816             :       !----------------
         817             :       ! Pointers
         818             :       !----------------
         819             :       real(r8), pointer, dimension(:) :: ztodt_ptr
         820             :       real(r8), pointer, dimension(:,:) :: thlm      ! Mean temperature
         821             :       real(r8), pointer, dimension(:,:) :: ice_supersat_frac ! ice cloud fraction
         822             :       real(r8), pointer, dimension(:,:) :: rtm       ! mean moisture mixing ratio
         823             :       real(r8), pointer, dimension(:,:) :: cld       ! CAM cloud fraction
         824             :       real(r8), pointer, dimension(:,:) :: alst      ! CLUBB liq cloud fraction
         825             :       real(r8), pointer, dimension(:,:) :: qrain     ! micro_mg rain from previous step
         826             :       real(r8), pointer, dimension(:,:) :: qsnow
         827             :       real(r8), pointer, dimension(:,:) :: nrain     ! micro_mg rain num conc
         828             :       real(r8), pointer, dimension(:,:) :: nsnow
         829             :
         830             :       real(r8), pointer, dimension(:,:) :: tke_in    ! TKE
         831             :       real(r8), pointer, dimension(:,:) :: khzm_in   ! Eddy diffusivity coef
         832             :
         833             :       logical, parameter :: l_est_kessler_microphys = .false.
         834             :       logical, parameter :: l_outfld_subcol         = .false.
         835             :
         836             :       type(grid) :: gr
         837             :
         838             :       type(precipitation_fractions) :: precip_fracs
         839             :
         840             :       ! Used as shortcuts to avoid typing hm_metadata%iiPDF_xx
         841             :       integer :: &
         842             :         iiPDF_chi, iiPDF_rr, iiPDF_w, iiPDF_Nr, &
         843             :         iiPDF_ri, iiPDF_Ni, iiPDF_Ncn, iiPDF_rs, iiPDF_Ns, &
         844             :         iirr, iiNr, iirs, iiri, &
         845             :         iirg, iiNs, iiNi, iiNg
         846             :
         847             :       !------------------------------------------------
         848             :       !                     Begin Code
         849             :       !------------------------------------------------
         850             :
         851             : #ifdef SILHS_OPENACC
         852             :       if ( l_est_kessler_microphys ) then
         853             :         call endrun('subcol_gen error: compilation with OpenACC requires l_est_kessler_microphys = .false.')
         854             :       end if
         855             :
         856             :       if ( subcol_SILHS_constrainmn ) then
         857             :         call endrun('subcol_gen error: compilation with OpenACC requires subcol_SILHS_constrainmn = .false.')
         858             :       end if
         859             :
         860             :       if ( subcol_SILHS_weight ) then
         861             :         call endrun('subcol_gen error: Importance sampling is not enabled for SILHS when using OpenACC. Set subcol_SILHS_weight to false.')
         862             :       end if
         863             : #endif
         864             :
         865             :       if (.not. allocated(state_sc%lat)) then
         866             :          call endrun('subcol_gen error: state_sc must be allocated before calling subcol_gen')
         867             :       end if
         868             :
         869             :       if( rx_Nc ) then
         870             :          call endrun('subcol_gen_SILHS: rx_Nc not enabled')
         871             :       endif
         872             :
         873             :       if (subcol_SILHS_meanice) then
         874             :         call endrun('subcol_gen_SILHS: subcol_SILHS_meanice = T not currently available')
         875             :       end if
         876             :
         877             :       ! Determine num of columns and which chunk we're working on and what timestep
         878             :       ngrdcol = state%ngrdcol
         879             :       ncol = state%ncol
         880             :       lchnk = state%lchnk
         881             :       iter = get_nstep() ! #KTCtodo: The model iteration is passed into SILHS without taking
         882             :                          !           substepping into account. I may need to change this in
         883             :                          !           the future. Also, why does SILHS need an iter, but CLUBB
         884             :                          !           does not?
         885             :                          ! #ERDBG:   The model iteration number is not used in SILHS unless
         886             :                          !           sequence_length > 1, but nobody runs with that option.
         887             :
         888             :       ! Copy hm_metadata indices to shortcuts
         889             :       iiPDF_chi = hm_metadata%iiPDF_chi
         890             :       iiPDF_Ncn = hm_metadata%iiPDF_Ncn
         891             :       iiPDF_rr  = hm_metadata%iiPDF_rr
         892             :       iiPDF_w   = hm_metadata%iiPDF_w
         893             :       iiPDF_Nr  = hm_metadata%iiPDF_Nr
         894             :       iiPDF_ri  = hm_metadata%iiPDF_ri
         895             :       iiPDF_Ni  = hm_metadata%iiPDF_Ni
         896             :       iiPDF_rs  = hm_metadata%iiPDF_rs
         897             :       iiPDF_Ns  = hm_metadata%iiPDF_Ns
         898             :       iirr      = hm_metadata%iirr
         899             :       iiNr      = hm_metadata%iiNr
         900             :       iirs      = hm_metadata%iirs
         901             :       iiri      = hm_metadata%iiri
         902             :       iirg      = hm_metadata%iirg
         903             :       iiNs      = hm_metadata%iiNs
         904             :       iiNi      = hm_metadata%iiNi
         905             :       iiNg      = hm_metadata%iiNg
         906             :
         907             :       !----------------
         908             :       ! Establish associations between pointers and physics buffer fields
         909             :       !----------------
         910             :       call pbuf_get_field(pbuf, thlm_idx, thlm)
         911             :       call pbuf_get_field(pbuf, ztodt_idx, ztodt_ptr)
         912             :       call pbuf_get_field(pbuf, ice_supersat_idx, ice_supersat_frac)
         913             :       call pbuf_get_field(pbuf, rtm_idx, rtm)
         914             :       call pbuf_get_field(pbuf, alst_idx, alst)
         915             :       call pbuf_get_field(pbuf, cld_idx, cld)
         916             :       call pbuf_get_field(pbuf, qrain_idx, qrain)
         917             :       call pbuf_get_field(pbuf, qsnow_idx, qsnow)
         918             :       call pbuf_get_field(pbuf, nrain_idx, nrain)
         919             :       call pbuf_get_field(pbuf, nsnow_idx, nsnow)
         920             :       call pbuf_get_field(pbuf, tke_idx, tke_in)
         921             :       call pbuf_get_field(pbuf, kvh_idx, khzm_in)
         922             :
         923             :       ! Pull c_K from clubb parameters.
         924             :       c_K = clubb_params_single_col(ic_K)
         925             :
         926             :       !----------------
         927             :       ! Copy state and populate numbers and values of sub-columns
         928             :       !----------------
         929             :       ztodt = ztodt_ptr(1)
         930             :       num_subcols = subcol_SILHS_numsubcol
         931             :
         932             :       ! The number of vertical grid levels used in CLUBB is pverp, which is originally
         933             :       ! set in the call to setup_clubb_core_api from subroutine clubb_ini_cam.  This
         934             :       ! is stored in CLUBB in the object gr%nz.  This isn't changed in CLUBB.
         935             :       ! However, when SILHS is used, SILHS only uses pverp - top_lev + 1 vertical grid
         936             :       ! levels and also uses the gr%nz object.  The value of gr%nz needs to be reset
         937             :       ! for SILHS here and then set again for CLUBB in subroutine clubb_tend_cam.
         938             :       gr%nz = pverp - top_lev + 1
         939             :
         940             :       ! Calculate sample weights separately at all grid levels when
         941             :       ! radiation is not called
         942             :       l_calc_weights_all_levs_itime = .false. ! subcol_utils cannot compute weighted avgs
         943             :                                               !   when the weights vary with height.
         944             :                                               !   Don't set to true until this is fixed!!
         945             :
         946             :
         947             :       ! Setup the CLUBB vertical grid object. This must be done for each
         948             :       ! column as the z-distance between hybrid pressure levels can
         949             :       ! change easily.
         950             :       ! Define the CLUBB momentum grid (in height, units of m)
         951             :       do k = 1, pverp-top_lev+1
         952             :         do i = 1, ngrdcol
         953             :           zi_g(i,k) = state%zi(i,pverp-k+1)-state%zi(i,pverp)
         954             :         end do
         955             :       end do
         956             :
         957             :
         958             :       ! Define the CLUBB thermodynamic grid (in units of m)
         959             :       do k = 1, pver-top_lev+1
         960             :         do i = 1, ngrdcol
         961             :           zt_g(i,k+1) = state%zm(i,pver-k+1)-state%zi(i,pverp)
         962             :
         963             :           ! Thermodynamic ghost point is below surface
         964             :           zt_g(i,1) = -1._r8*zt_g(i,2)
         965             :         end do
         966             :       end do
         967             :
         968             :       do i=1, ncol
         969             :         !  Set the elevation of the surface
         970             :         sfc_elevation(i) = state%zi(i,pver+1)
         971             :       end do
         972             :
         973             :       !  Heights need to be set at each timestep.
         974             :       call setup_grid_api( pverp+1-top_lev, ncol, sfc_elevation(1:ncol), l_implemented,  & ! intent(in)
         975             :                            grid_type, zi_g(1:ncol,2), zi_g(1:ncol,1), zi_g(1:ncol,pverp+1-top_lev),   & ! intent(in)
         976             :                            zi_g(1:ncol,:), zt_g(1:ncol,:),                              & ! intent(in)
         977             :                            gr )
         978             :
         979             :       ! Calculate the distance between grid levels on the host model grid,
         980             :       ! using host model grid indices.
         981             :       do k = top_lev, pver
         982             :         do i = 1, ngrdcol
         983             :           dz_g(i,k) = state%zi(i,k)-state%zi(i,k+1)
         984             :         end do
         985             :       end do
         986             :
         987             :       ! Inverse delta_zm is required for the 3-level L-scale averaging
         988             :       do k = 1, pver-top_lev+1
         989             :         do i = 1, ngrdcol
         990             :           delta_zm(i,k+1) = state%zi(i,pverp-k)-state%zi(i,pverp-k+1)
         991             :
         992             :           ! Handle CLUBB sub-sfc ghost point as done in clubb grid_class.F90
         993             :           delta_zm(i,1) = delta_zm(i,2)
         994             :         end do
         995             :       end do
         996             :
         997             :       ! Compute dry static density on CLUBB vertical grid
         998             :       do k = 1, pver-top_lev+1
         999             :         do i = 1, ngrdcol
        1000             :           rho_ds_zt(i,k+1) = (rga)*state%pdel(i,pverp-k)/dz_g(i,pverp-k)
        1001             :
        1002             :           ! CLUBB ghost point under the surface
        1003             :           rho_ds_zt(i,1) = rho_ds_zt(i,2)
        1004             :         end do
        1005             :       end do
        1006             :
        1007             :       ! Set up hydromet array, flipped from CAM vert grid to CLUBB
        1008             :       if ( iirr > 0 ) then
        1009             :         ! If ixrain and family are greater than zero, then MG2 is
        1010             :         ! being used, and rain and snow are part of state. Otherwise,
        1011             :         ! diagnostic rain and snow from MG1 are used in hydromet.
        1012             :         if (ixrain > 0) then
        1013             :           do k = 1, pver-top_lev+1
        1014             :             do i = 1, ngrdcol
        1015             :               hydromet(i,k+1,iirr) = state%q(i,pver-k+1,ixrain)
        1016             :             end do
        1017             :           end do
        1018             :         else
        1019             :           do k = 1, pver-top_lev+1
        1020             :             do i = 1, ngrdcol
        1021             :               hydromet(i,k+1,iirr) = qrain(i,pver-k+1)
        1022             :             end do
        1023             :           end do
        1024             :         endif
        1025             :       endif
        1026             :
        1027             :       if ( iiNr > 0 ) then
        1028             :         if (ixnumrain > 0) then
        1029             :           do k = 1, pver-top_lev+1
        1030             :             do i = 1, ngrdcol
        1031             :               hydromet(i,k+1,iiNr) = state%q(i,pver-k+1,ixnumrain)
        1032             :             end do
        1033             :           end do
        1034             :         else
        1035             :           do k = 1, pver-top_lev+1
        1036             :             do i = 1, ngrdcol
        1037             :               hydromet(i,k+1,iiNr) = nrain(i,pver-k+1)
        1038             :             end do
        1039             :           end do
        1040             :         endif
        1041             :       endif
        1042             :
        1043             :       if ( iirs > 0 ) then
        1044             :         if (ixsnow > 0) then
        1045             :           do k = 1, pver-top_lev+1
        1046             :             do i = 1, ngrdcol
        1047             :               hydromet(i,k+1,iirs) = state%q(i,pver-k+1,ixsnow)
        1048             :             end do
        1049             :           end do
        1050             :         else
        1051             :           do k = 1, pver-top_lev+1
        1052             :             do i = 1, ngrdcol
        1053             :               hydromet(i,k+1,iirs) = qsnow(i,pver-k+1)
        1054             :             end do
        1055             :           end do
        1056             :         endif
        1057             :       endif
        1058             :
        1059             :       if ( iiNs > 0 ) then
        1060             :         if (ixnumsnow > 0) then
        1061             :           do k = 1, pver-top_lev+1
        1062             :             do i = 1, ngrdcol
        1063             :               hydromet(i,k+1,iiNs) = state%q(i,pver-k+1,ixnumsnow)
        1064             :             end do
        1065             :           end do
        1066             :         else
        1067             :           do k = 1, pver-top_lev+1
        1068             :             do i = 1, ngrdcol
        1069             :               hydromet(i,k+1,iiNs) = nsnow(i,pver-k+1)
        1070             :             end do
        1071             :           end do
        1072             :         endif
        1073             :       endif
        1074             :
        1075             :       if ( iiri > 0 ) then
        1076             :         do k = 1, pver-top_lev+1
        1077             :           do i = 1, ngrdcol
        1078             :             hydromet(i,k+1,iiri) = state%q(i,pver-k+1,ixcldice)
        1079             :           end do
        1080             :         end do
        1081             :       endif
        1082             :
        1083             :       if ( iiNi > 0 ) then
        1084             :         do k = 1, pver-top_lev+1
        1085             :           do i = 1, ngrdcol
        1086             :             hydromet(i,k+1,iiNi) = state%q(i,pver-k+1,ixnumice)
        1087             :           end do
        1088             :         end do
        1089             :       endif
        1090             :
        1091             :       do k = 1, hydromet_dim ! ghost point below the surface
        1092             :         do i = 1, ngrdcol
        1093             :           hydromet(i,1,k) = hydromet(i,2,k)
        1094             :         end do
        1095             :       end do
        1096             :
        1097             :       do k = 1, pver-top_lev+1
        1098             :         do i = 1, ngrdcol
        1099             :           Ncm(i,k+1) = state%q(i,pver-k+1,ixnumliq)
        1100             :           Ncm(i,1) = Ncm(i,2)
        1101             :        end do
        1102             :       end do
        1103             :
        1104             :       ! Convert from CAM vertical grid to CLUBB
        1105             :       do k = 1, pverp-top_lev+1
        1106             :         do i = 1, ngrdcol
        1107             :           ice_supersat_frac_in(i,k) = ice_supersat_frac(i,pverp-k+1)
        1108             :         end do
        1109             :       end do
        1110             :
        1111             :
        1112             :       do k = 1, pver-top_lev+1
        1113             :         do i = 1, ngrdcol
        1114             :           cld_frac_in(i,k+1) = alst(i,pver-k+1)
        1115             :           cld_frac_in(i,1) = cld_frac_in(i,2) ! Ghost pt below surface
        1116             :         end do
        1117             :       end do
        1118             :
        1119             :       ! Calculate a clubb-specific exner function
        1120             :       ! (This is grid mean, as pressure levels do not change in
        1121             :       !  the subcolumn state)
        1122             :       do k = 1, pver-top_lev+1
        1123             :         do i = 1, ngrdcol
        1124             :           invs_exner(i,k) = ((state%pmid(i,k)/p0_clubb)**(cappa))
        1125             :         end do
        1126             :       end do
        1127             :
        1128             :       ! Call setup_pdf_parameters to get the CLUBB PDF ready for SILHS
        1129             :       ! Compute Num concentration of cloud nuclei
        1130             :       do k = 1, pverp-top_lev+1
        1131             :         do i = 1, ngrdcol
        1132             :           Nc_in_cloud(i,k) = Ncm(i,k) / max( cld_frac_in(i,k), cloud_frac_min )
        1133             :         end do
        1134             :       end do
        1135             :
        1136             :       ! The variable wphydrometp is only used when l_calc_w_corr is enabled.
        1137             :       ! The l_calc_w_corr flag is turned off by default, so wphydrometp will
        1138             :       ! simply be set to 0 to simplify matters.
        1139             :       wphydrometp = 0.0_r8
        1140             :
        1141             :       do k = 1, pverp-top_lev+1
        1142             :         do i = 1, ngrdcol
        1143             :           khzm(i,k) = khzm_in(i,pverp-k+1)
        1144             :         end do
        1145             :       end do
        1146             :
        1147             :       ! Allocate 2D arrays in precip_fracs for all grid columns and vertical levels
        1148             :       call init_precip_fracs_api( pverp-top_lev+1, ngrdcol, &
        1149             :                                   precip_fracs )
        1150             :
        1151             :       call setup_pdf_parameters_api( gr, pverp-top_lev+1, ngrdcol, pdf_dim, hydromet_dim, ztodt,  & ! In
        1152             :                                      Nc_in_cloud, cld_frac_in, khzm, &                              ! In
        1153             :                                      ice_supersat_frac_in, hydromet, wphydrometp, &                 ! In
        1154             :                                      corr_array_n_cloud, corr_array_n_below, &                      ! In
        1155             :                                      hm_metadata, &                                                 ! In
        1156             :                                      pdf_params_chnk(lchnk), &                                      ! In
        1157             :                                      clubb_params_single_col, &                                     ! In
        1158             :                                      clubb_config_flags%iiPDF_type, &                               ! In
        1159             :                                      clubb_config_flags%l_use_precip_frac, &                        ! In
        1160             :                                      clubb_config_flags%l_predict_upwp_vpwp, &                      ! In
        1161             :                                      clubb_config_flags%l_diagnose_correlations, &                  ! In
        1162             :                                      clubb_config_flags%l_calc_w_corr, &                            ! In
        1163             :                                      clubb_config_flags%l_const_Nc_in_cloud, &                      ! In
        1164             :                                      clubb_config_flags%l_fix_w_chi_eta_correlations, &             ! In
        1165             :                                      stats_metadata, &                                              ! In
        1166             :                                      stats_zt, stats_zm, stats_sfc, &                               ! In
        1167             :                                      hydrometp2, &                                                  ! Inout
        1168             :                                      mu_x_1, mu_x_2, &                                              ! Out
        1169             :                                      sigma_x_1, sigma_x_2, &                                        ! Out
        1170             :                                      corr_array_1, corr_array_2, &                                  ! Out
        1171             :                                      corr_cholesky_mtx_1, corr_cholesky_mtx_2, &                    ! Out
        1172             :                                      precip_fracs )                                                 ! Inout
        1173             :
        1174             :       ! In order for Lscale to be used properly, it needs to be passed out of
        1175             :       ! advance_clubb_core, saved to the pbuf, and then pulled out of the
        1176             :       ! pbuf for use here.  The profile of Lscale is passed into subroutine
        1177             :       ! generate_silhs_sample_api for use in calculating the vertical
        1178             :       ! correlation coefficient.  Rather than output Lscale directly, its
        1179             :       ! value can be calculated from other fields that are already output to
        1180             :       ! pbuf.  The equation relating Lscale to eddy diffusivity is:
        1181             :       !
        1182             :       ! Kh = c_K * Lscale * sqrt( TKE ).
        1183             :       !
        1184             :       ! Both Kh and TKE are written to the pbuf, and c_K is easily extracted
        1185             :       ! from CLUBB's tunable parameters.  The equation for Lscale is:
        1186             :       !
        1187             :       ! Lscale = Kh / ( c_K * sqrt( TKE ) ).
        1188             :       !
        1189             :       ! Since Kh and TKE are output on momentum (interface) grid levels, the
        1190             :       ! resulting calculation of Lscale is also found on momentum levels.  It
        1191             :       ! needs to be interpolated back to thermodynamic (midpoint) grid levels
        1192             :       ! for further use.
        1193             :       do k = 1, pverp-top_lev+1
        1194             :         do i = 1, ngrdcol
        1195             :           tke(i,k) = tke_in(i,pverp-k+1)
        1196             :         end do
        1197             :       end do
        1198             :
        1199             :       do k = 1, pverp-top_lev+1
        1200             :         do i = 1, ngrdcol
        1201             :           Lscale_zm(i,k) = khzm(i,k) / ( c_K * sqrt( max( tke(i,k), em_min ) ) )
        1202             :         end do
        1203             :       end do
        1204             :
        1205             :       do i = 1, ngrdcol
        1206             :         Lscale(i,1) = Lscale_zm(i,1) + ( Lscale_zm(i,2) - Lscale_zm(i,1) ) &
        1207             :                                      * ( zt_g(i,1) - zi_g(i,1) ) / ( zi_g(i,2) - zi_g(i,1) )
        1208             :       end do
        1209             :
        1210             :       do k = 2, pverp-top_lev+1
        1211             :         do i = 1, ngrdcol
        1212             :           Lscale(i,k) = Lscale_zm(i,k-1) + ( Lscale_zm(i,k) - Lscale_zm(i,k-1) ) &
        1213             :                                          * ( zt_g(i,k) - zi_g(i,k-1) ) / ( zi_g(i,k) - zi_g(i,k-1) )
        1214             :         end do
        1215             :       end do
        1216             :
        1217             :       do k = 2, pverp-top_lev+1
        1218             :         do i = 1, ngrdcol
        1219             :           Lscale(i,:) = max( Lscale(i,:), 0.01_r8 )
        1220             :         end do
        1221             :       end do
        1222             :
        1223             :       !$acc data create( X_mixt_comp_all_levs, X_nl_all_levs, lh_rc_clipped, lh_Nc_clipped, &
        1224             :       !$acc&             lh_sample_point_weights, lh_rt_clipped, lh_rt_clipped, &
        1225             :       !$acc&             lh_rv_clipped, lh_thl_clipped, THL_lh_out, &
        1226             :       !$acc&             RT_lh_out, RCM_lh_out, NCLW_lh_out, ICE_lh_out, &
        1227             :       !$acc&             NICE_lh_out, RVM_lh_out, THL_lh_out, RAIN_lh_out, &
        1228             :       !$acc&             NRAIN_lh_out, SNOW_lh_out, NSNOW_lh_out, WM_lh_out, &
        1229             :       !$acc&             OMEGA_lh_out ) &
        1230             :       !$acc&     copyin( state, state%zm, state%phis, rho_ds_zt, invs_exner ) &
        1231             :       !$acc&     copyout( state%t, state%s, state%omega, state_sc%q )
        1232             :       !$acc& async(1)
        1233             :
        1234             :       ! Set the seed to the random number generator based on a quantity that
        1235             :       ! will be reproducible for restarts.
        1236             :       lh_seed = int( 1.0e4_r8 * rtm(1,pver), kind = genrand_intg )
        1237             :
        1238             :       ! Let's generate some subcolumns!!!!!
        1239             :       call generate_silhs_sample_api( &
        1240             :                     iter, pdf_dim, num_subcols, sequence_length, pverp-top_lev+1, ngrdcol, & ! In
        1241             :                     l_calc_weights_all_levs_itime, &                      ! In
        1242             :                     pdf_params_chnk(lchnk), delta_zm, Lscale, &           ! In
        1243             :                     lh_seed, hm_metadata, &                               ! In
        1244             :                     rho_ds_zt, &                                          ! In
        1245             :                     mu_x_1, mu_x_2, sigma_x_1, sigma_x_2, &               ! In
        1246             :                     corr_cholesky_mtx_1, corr_cholesky_mtx_2, &           ! In
        1247             :                     precip_fracs, silhs_config_flags, &                   ! In
        1248             :                     vert_decorr_coef, &                                   ! In
        1249             :                     stats_metadata, &                                     ! In
        1250             :                     stats_lh_zt, stats_lh_sfc, &                          ! InOut
        1251             :                     X_nl_all_levs, X_mixt_comp_all_levs, &                ! Out
        1252             :                     lh_sample_point_weights)                              ! Out
        1253             :
        1254             :       ! Extract clipped variables from subcolumns
        1255             :       call clip_transform_silhs_output_api( gr, pverp-top_lev+1, ngrdcol, num_subcols, & ! In
        1256             :                                             pdf_dim, hydromet_dim, hm_metadata,        & ! In
        1257             :                                             X_mixt_comp_all_levs,                      & ! In
        1258             :                                             X_nl_all_levs,                             & ! In
        1259             :                                             pdf_params_chnk(lchnk),                    & ! In
        1260             :                                             l_use_Ncn_to_Nc,                           & ! In
        1261             :                                             lh_rt_clipped, lh_thl_clipped,             & ! Out
        1262             :                                             lh_rc_clipped, lh_rv_clipped,              & ! Out
        1263             :                                             lh_Nc_clipped )                              ! Out
        1264             :       !$acc wait
        1265             :
        1266             :       if ( l_est_kessler_microphys ) then
        1267             :         call endrun('subcol_SILHS: l_est_kessler_microphys = T is not currently supported')
        1268             :       end if
        1269             :
        1270             :       !-------------------------------------------------------------------------
        1271             :       !            Convert from CLUBB vertical grid to CAM grid
        1272             :       !------------------------------------------------------------------------
        1273             :       ! This kernel is executed in stream 1:
        1274             :       !$acc parallel loop collapse(3) default(present) async(1)
        1275             :       do k = top_lev, pverp
        1276             :         do j = 1, num_subcols
        1277             :           do i = 1, ngrdcol
        1278             :             RT_lh_out(   num_subcols*(i-1)+j,k ) = lh_rt_clipped(i,j,pverp-k+1)
        1279             :             RCM_lh_out(  num_subcols*(i-1)+j,k ) = lh_rc_clipped(i,j,pverp-k+1)
        1280             :             NCLW_lh_out( num_subcols*(i-1)+j,k ) = lh_Nc_clipped(i,j,pverp-k+1)
        1281             :             RVM_lh_out(  num_subcols*(i-1)+j,k ) = lh_rv_clipped(i,j,pverp-k+1)
        1282             :             THL_lh_out(  num_subcols*(i-1)+j,k ) = lh_thl_clipped(i,j,pverp-k+1)
        1283             :           end do
        1284             :         end do
        1285             :       end do
        1286             :
        1287             :       ! This kernel is executed in stream 2:
        1288             :       !$acc parallel loop collapse(3) default(present) async(2)
        1289             :       do k = top_lev, pverp
        1290             :         do j = 1, num_subcols
        1291             :           do i = 1, ngrdcol
        1292             :             ICE_lh_out(   num_subcols*(i-1)+j,k ) = X_nl_all_levs(i,j,pverp-k+1,iiPDF_ri)
        1293             :             NICE_lh_out(  num_subcols*(i-1)+j,k ) = X_nl_all_levs(i,j,pverp-k+1,iiPDF_Ni)
        1294             :             RAIN_lh_out(  num_subcols*(i-1)+j,k ) = X_nl_all_levs(i,j,pverp-k+1,iiPDF_rr)
        1295             :             NRAIN_lh_out( num_subcols*(i-1)+j,k ) = X_nl_all_levs(i,j,pverp-k+1,iiPDF_Nr)
        1296             :             SNOW_lh_out(  num_subcols*(i-1)+j,k ) = X_nl_all_levs(i,j,pverp-k+1,iiPDF_rs)
        1297             :             NSNOW_lh_out( num_subcols*(i-1)+j,k ) = X_nl_all_levs(i,j,pverp-k+1,iiPDF_Ns)
        1298             :             WM_lh_out(    num_subcols*(i-1)+j,k ) = X_nl_all_levs(i,j,pverp-k+1,iiPDF_w)
        1299             :           end do
        1300             :         end do
        1301             :       end do
        1302             :
        1303             :       ! This kernel is executed in stream 2 because WM_lh_out comes from stream 2:
        1304             :       !$acc parallel loop collapse(3) default(present) async(2)
        1305             :       do k = top_lev, pverp
        1306             :         do j = 1, num_subcols
        1307             :           do i = 1, ngrdcol
        1308             :             OMEGA_lh_out( num_subcols*(i-1)+j,k ) = -1._r8 * WM_lh_out(num_subcols*(i-1)+j,k) &
        1309             :                                                            * rho_ds_zt(i,pverp-k+1) * gravit
        1310             :           end do
        1311             :         end do
        1312             :       end do
        1313             :
        1314             :       if ( l_est_kessler_microphys ) then
        1315             :         do k = top_lev, pverp
        1316             :           do j = 1, num_subcols
        1317             :             do i = 1, ngrdcol
        1318             :               AKm_out(i,k) = AKm(i,pverp-k+1)
        1319             :               lh_AKm_out(i,k) = lh_AKm(i,pverp-k+1)
        1320             :             end do
        1321             :           end do
        1322             :         end do
        1323             :       end if
        1324             :
        1325             :       ! Pack up weights
        1326             :       ! Using grid level 2 always won't work if weights vary with height.
        1327             :       call subcol_pack(lchnk, lh_sample_point_weights(:,:,2), weights )
        1328             :       call subcol_set_weight(lchnk, weights)
        1329             :
        1330             :       ! Constrain the sample distribution of cloud water and ice to the same mean
        1331             :       ! as the grid to prevent negative condensate errors
        1332             :       if(subcol_SILHS_constrainmn) then
        1333             :
        1334             :         do i = 1, ngrdcol
        1335             :
        1336             :           stncol = num_subcols*(i-1)
        1337             :
        1338             :           call subcol_constrainmn( num_subcols, ICE_lh_out(stncol+1:stncol+num_subcols,:), &
        1339             :                                    weights(stncol+1:stncol+num_subcols), &
        1340             :                                    state%q(i,:,ixcldice), meansc_ice(i,:), stdsc_ice(i,:) )
        1341             :           if ( ixrain > 0 ) &
        1342             :           call subcol_constrainmn( num_subcols, RAIN_lh_out(stncol+1:stncol+num_subcols,:), &
        1343             :                                    weights(stncol+1:stncol+num_subcols), &
        1344             :                                    state%q(i,:,ixrain) )
        1345             :           if ( ixsnow > 0 ) &
        1346             :           call subcol_constrainmn( num_subcols, SNOW_lh_out(stncol+1:stncol+num_subcols,:), &
        1347             :                                    weights(stncol+1:stncol+num_subcols), &
        1348             :                                    state%q(i,:,ixsnow) )
        1349             :           call subcol_constrainmn( num_subcols, RCM_lh_out(stncol+1:stncol+num_subcols,:), &
        1350             :                                    weights(stncol+1:stncol+num_subcols), &
        1351             :                                    state%q(i,:,ixcldliq), meansc_liq(i,:), stdsc_liq(i,:) )
        1352             :           call subcol_constrainmn( num_subcols, RVM_lh_out(stncol+1:stncol+num_subcols,:), &
        1353             :                                    weights(stncol+1:stncol+num_subcols), &
        1354             :                                    state%q(i,:,ixq), meansc_vap(i,:), stdsc_vap(i,:) )
        1355             :           call subcol_constrainmn( num_subcols, NICE_lh_out(stncol+1:stncol+num_subcols,:), &
        1356             :                                    weights(stncol+1:stncol+num_subcols), &
        1357             :                                    state%q(i,:,ixnumice) )
        1358             :           if ( ixnumrain > 0 ) &
        1359             :           call subcol_constrainmn( num_subcols, NRAIN_lh_out(stncol+1:stncol+num_subcols,:), &
        1360             :                                    weights(stncol+1:stncol+num_subcols), &
        1361             :                                    state%q(i,:,ixnumrain) )
        1362             :           if ( ixnumsnow > 0 ) &
        1363             :           call subcol_constrainmn( num_subcols, NSNOW_lh_out(stncol+1:stncol+num_subcols,:), &
        1364             :                                    weights(stncol+1:stncol+num_subcols), &
        1365             :                                    state%q(i,:,ixnumsnow) )
        1366             :           call subcol_constrainmn( num_subcols, NCLW_lh_out(stncol+1:stncol+num_subcols,:), &
        1367             :                                    weights(stncol+1:stncol+num_subcols), &
        1368             :                                    state%q(i,:,ixnumliq) )
        1369             :           do k = top_lev, pver
        1370             :             ! Look for exceptionally large values of condensate
        1371             :             if(ANY(ICE_lh_out(stncol+1:stncol+num_subcols,k) .gt. 0.01_r8)) then
        1372             :               ! Clip the large values
        1373             :               where(ICE_lh_out(stncol+1:stncol+num_subcols,k) .gt. 0.01_r8)
        1374             :                  ICE_lh_out(stncol+1:stncol+num_subcols,k) = 0.01_r8
        1375             :                  NICE_lh_out(stncol+1:stncol+num_subcols,k) = 1.5e+7_r8
        1376             :               end where
        1377             :               ! Recalculate the weighted subcolumn mean
        1378             :               tmp_mean = meansc( ICE_lh_out( stncol+1:stncol+num_subcols, k ), &
        1379             :                                      weights(stncol+1:stncol+num_subcols), &
        1380             :                                      real(num_subcols,r8) )
        1381             :               ! Calculate the difference between the weighted mean and grid mean
        1382             :               diff_mean = state%q(i,k,ixcldice)-tmp_mean
        1383             :               ! Add the difference to each subcolumn
        1384             :               ICE_lh_out(stncol+1:stncol+num_subcols,k) = &
        1385             :                  ICE_lh_out(stncol+1:stncol+num_subcols,k)+diff_mean
        1386             :               ! Recalculate the weight subcolumn mean for ice num conc
        1387             :               tmp_mean = meansc( NICE_lh_out( stncol+1:stncol+num_subcols, k ), &
        1388             :                                      weights(stncol+1:stncol+num_subcols), &
        1389             :                                      real(num_subcols,r8) )
        1390             :               ! Calculate the difference between the weighted mean and grid mean
        1391             :               diff_mean = state%q(i,k,ixnumice)-tmp_mean
        1392             :               ! Add the difference to each subcolumn
        1393             :               if(diff_mean.gt.0.0_r8) then
        1394             :                  NICE_lh_out(stncol+1:stncol+num_subcols,k) = &
        1395             :                      NICE_lh_out(stncol+1:stncol+num_subcols,k)+diff_mean
        1396             :               else ! just use the grid mean in each subcolumn
        1397             :                  NICE_lh_out(stncol+1:stncol+num_subcols,k) = &
        1398             :                      state%q(i,k,ixnumice)
        1399             :               end if
        1400             :               ! Test adjusted means for debugging
        1401             :               tmp_mean = meansc( ICE_lh_out( stncol+1:stncol+num_subcols, k ), &
        1402             :                                      weights(stncol+1:stncol+num_subcols), &
        1403             :                                      real(num_subcols,r8) )
        1404             :               diff_mean = state%q(i,k,ixcldice)-tmp_mean
        1405             :               tmp_mean = meansc( NICE_lh_out( stncol+1:stncol+num_subcols, k ), &
        1406             :                                      weights(stncol+1:stncol+num_subcols), &
        1407             :                                      real(num_subcols,r8) )
        1408             :               diff_mean = state%q(i,k,ixnumice)-tmp_mean
        1409             :             endif
        1410             :           end do ! k = top_lev, pver
        1411             :         end do
        1412             :       endif ! subcol_silhs_constrainm
        1413             :
        1414             :
        1415             :       !---------------------------------------------------
        1416             :       !            Updating state variables
        1417             :       !---------------------------------------------------
        1418             :       ! Code to update the state variables for interactive runs
        1419             :       ! This kernel is executed in stream 3, but waits for stream 1
        1420             :       ! because THL_lh_out and RCM_lh_out come from stream 1:
        1421             :       !$acc parallel loop collapse(3) default(present) wait(1) async(3)
        1422             :       do k = 1, pver-top_lev+1
        1423             :         do j = 1, num_subcols
        1424             :           do i = 1, ngrdcol
        1425             :
        1426             :             state_sc%t(num_subcols*(i-1)+j,k) = THL_lh_out(num_subcols*(i-1)+j,k) * invs_exner(i,k) &
        1427             :                                                 + Lv * RCM_lh_out(num_subcols*(i-1)+j,k) / Cp
        1428             :
        1429             :             state_sc%s(num_subcols*(i-1)+j,k) = cpair * state_sc%t(num_subcols*(i-1)+j,k) &
        1430             :                                                 + gravit * state%zm(i,k) + state%phis(i)
        1431             :           end do
        1432             :         end do
        1433             :       end do
        1434             :
        1435             :       ! This kernel is executed in stream 4, but waits for stream 1 and 2
        1436             :       ! because RVM_lh_out is from stream 1 and OMEGA_lh_out is from stream 2:
        1437             :       !$acc parallel loop collapse(3) default(present) wait(1,2) async(4)
        1438             :       do k = 1, pver-top_lev+1
        1439             :         do j = 1, num_subcols
        1440             :           do i = 1, ngrdcol
        1441             :             ! Vertical Velocity is not part of the energy conservation checks, but
        1442             :             ! we need to be careful here, because the SILHS output VV is noisy.
        1443             :             state_sc%omega(num_subcols*(i-1)+j,k) = OMEGA_lh_out(num_subcols*(i-1)+j,k)
        1444             :             state_sc%q(num_subcols*(i-1)+j,k,ixq) = RVM_lh_out(num_subcols*(i-1)+j,k)
        1445             :           end do
        1446             :         end do
        1447             :       end do
        1448             :
        1449             :
        1450             :       if (subcol_SILHS_q_to_micro) then ! Send SILHS predicted constituents to microp
        1451             :
        1452             :         ! This kernel is executed in stream 5, but waits for stream 1 and 2
        1453             :         ! because RCM_lh_out is from stream 1 and ICE_lh_out is from stream 2:
        1454             :         !$acc parallel loop collapse(3) default(present) wait(1,2) async(5)
        1455             :         do k = 1, pver-top_lev+1
        1456             :           do j = 1, num_subcols
        1457             :             do i = 1, ngrdcol
        1458             :               state_sc%q(num_subcols*(i-1)+j,k,ixcldliq) = RCM_lh_out(num_subcols*(i-1)+j,k)
        1459             :               state_sc%q(num_subcols*(i-1)+j,k,ixcldice) = ICE_lh_out(num_subcols*(i-1)+j,k)
        1460             :             end do
        1461             :           end do
        1462             :         end do
        1463             :
        1464             :         if (ixrain > 0) then
        1465             :           ! This kernel is executed in stream 6, but waits for stream 2
        1466             :           ! because RAIN_lh_out is from stream 2:
        1467             :           !$acc parallel loop collapse(3) default(present) wait(2) async(6)
        1468             :           do k = 1, pver-top_lev+1
        1469             :             do j = 1, num_subcols
        1470             :               do i = 1, ngrdcol
        1471             :                 state_sc%q(num_subcols*(i-1)+j,k,ixrain) = RAIN_lh_out(num_subcols*(i-1)+j,k)
        1472             :               end do
        1473             :             end do
        1474             :           end do
        1475             :         end if
        1476             :
        1477             :         if (ixsnow > 0) then
        1478             :           ! This kernel is executed in stream 7, but waits for stream 2
        1479             :           ! because SNOW_lh_out is from stream 2:
        1480             :           !$acc parallel loop collapse(3) default(present) wait(2) async(7)
        1481             :           do k = 1, pver-top_lev+1
        1482             :             do j = 1, num_subcols
        1483             :               do i = 1, ngrdcol
        1484             :                 state_sc%q(num_subcols*(i-1)+j,k,ixsnow) = SNOW_lh_out(num_subcols*(i-1)+j,k)
        1485             :               end do
        1486             :             end do
        1487             :           end do
        1488             :         end if
        1489             :
        1490             :       else
        1491             :
        1492             :         do k = 1, pver-top_lev+1
        1493             :           do j = 1, num_subcols
        1494             :             do i = 1, ngrdcol
        1495             :               state_sc%q(num_subcols*(i-1)+j,k,ixcldliq) = state%q(i,k,ixcldliq)
        1496             :               state_sc%q(num_subcols*(i-1)+j,k,ixcldice) = state%q(i,k,ixcldice)
        1497             :               if (ixrain > 0) then
        1498             :                  state_sc%q(num_subcols*(i-1)+j,k,ixrain) = state%q(i,k,ixrain)
        1499             :               end if
        1500             :               if (ixsnow > 0) then
        1501             :                  state_sc%q(num_subcols*(i-1)+j,k,ixsnow) = state%q(i,k,ixsnow)
        1502             :               end if
        1503             :             end do
        1504             :           end do
        1505             :         end do
        1506             :
        1507             :       endif
        1508             :
        1509             :       if (subcol_SILHS_n_to_micro) then ! Send SILHS predicted number conc to microp
        1510             :
        1511             :         ! This kernel is executed in stream 8, but waits for stream 1 and 2
        1512             :         ! because NCLW_lh_out is from stream 1 and NICE_lh_out is from stream 2:
        1513             :         !$acc parallel loop collapse(3) default(present) wait(1,2) async(8)
        1514             :         do k = 1, pver-top_lev+1
        1515             :           do j = 1, num_subcols
        1516             :             do i = 1, ngrdcol
        1517             :               state_sc%q(num_subcols*(i-1)+j,k,ixnumice) = NICE_lh_out(num_subcols*(i-1)+j,k)
        1518             :               state_sc%q(num_subcols*(i-1)+j,k,ixnumliq) = NCLW_lh_out(num_subcols*(i-1)+j,k)
        1519             :             end do
        1520             :           end do
        1521             :         end do
        1522             :
        1523             :         if (ixnumrain > 0) then
        1524             :           ! This kernel is executed in stream 9, but waits for stream 2
        1525             :           ! because NRAIN_lh_out is from stream 2:
        1526             :           !$acc parallel loop collapse(3) default(present) wait(2) async(9)
        1527             :           do k = 1, pver-top_lev+1
        1528             :             do j = 1, num_subcols
        1529             :               do i = 1, ngrdcol
        1530             :                 state_sc%q(num_subcols*(i-1)+j,k,ixnumrain) = NRAIN_lh_out(num_subcols*(i-1)+j,k)
        1531             :               end do
        1532             :             end do
        1533             :           end do
        1534             :         end if
        1535             :
        1536             :         if (ixnumsnow > 0) then
        1537             :           ! This kernel is executed in stream 10, but waits for stream 2
        1538             :           ! because NSNOW_lh_out is from stream 2:
        1539             :           !$acc parallel loop collapse(3) default(present) wait(2) async(10)
        1540             :           do k = 1, pver-top_lev+1
        1541             :             do j = 1, num_subcols
        1542             :               do i = 1, ngrdcol
        1543             :                 state_sc%q(num_subcols*(i-1)+j,k,ixnumsnow) = NSNOW_lh_out(num_subcols*(i-1)+j,k)
        1544             :               end do
        1545             :             end do
        1546             :           end do
        1547             :         end if
        1548             :
        1549             :       else
        1550             :
        1551             :         do k = 1, pver-top_lev+1
        1552             :           do j = 1, num_subcols
        1553             :             do i = 1, ngrdcol
        1554             :               state_sc%q(num_subcols*(i-1)+j,k,ixnumliq) = state%q(i,k,ixnumliq)
        1555             :               state_sc%q(num_subcols*(i-1)+j,k,ixnumice) = state%q(i,k,ixnumice)
        1556             :               if (ixnumrain > 0) then
        1557             :                  state_sc%q(num_subcols*(i-1)+j,k,ixnumrain) = state%q(i,k,ixnumrain)
        1558             :               end if
        1559             :               if (ixnumsnow > 0) then
        1560             :                  state_sc%q(num_subcols*(i-1)+j,k,ixnumsnow) = state%q(i,k,ixnumsnow)
        1561             :               end if
        1562             :             end do
        1563             :           end do
        1564             :         end do
        1565             :
        1566             :       endif
        1567             :
        1568             :       ! This kernel is executed in stream 8, because state_sc%q(:,:,ixnumliq) and
        1569             :       !  state_sc%q(:,:,ixnumice) are from stream 8
        1570             :       !$acc parallel loop collapse(3) default(present) async(8)
        1571             :       do k = 1, pver-top_lev+1
        1572             :         do j = 1, num_subcols
        1573             :           do i = 1, ngrdcol
        1574             :             ! Change liq and ice (and rain and snow) num conc zeros to min values (1e-12)
        1575             :             if (state_sc%q(num_subcols*(i-1)+j,k,ixnumliq) .lt. min_num_conc) then
        1576             :                 state_sc%q(num_subcols*(i-1)+j,k,ixnumliq) = min_num_conc
        1577             :             end if
        1578             :
        1579             :             if (state_sc%q(num_subcols*(i-1)+j,k,ixnumice) .lt. min_num_conc) then
        1580             :                 state_sc%q(num_subcols*(i-1)+j,k,ixnumice) = min_num_conc
        1581             :             end if
        1582             :           end do
        1583             :         end do
        1584             :       end do
        1585             :
        1586             :       if (ixnumrain > 0) then
        1587             :         ! This kernel is executed in stream 9, because state_sc%q(:,:,ixnumrain) is
        1588             :         ! from stream 9
        1589             :         !$acc parallel loop collapse(3) default(present) async(9)
        1590             :         do k = 1, pver-top_lev+1
        1591             :           do j = 1, num_subcols
        1592             :             do i = 1, ngrdcol
        1593             :                if(state_sc%q(num_subcols*(i-1)+j,k,ixnumrain) .lt. min_num_conc) then
        1594             :                   state_sc%q(num_subcols*(i-1)+j,k,ixnumrain) = min_num_conc
        1595             :                end if
        1596             :             end do
        1597             :           end do
        1598             :         end do
        1599             :       endif
        1600             :
        1601             :       if (ixnumsnow > 0) then
        1602             :         ! This kernel is executed in stream 10, because state_sc%q(:,:,ixnumsnow) is
        1603             :         ! from stream 10
        1604             :         !$acc parallel loop collapse(3) default(present) async(10)
        1605             :         do k = 1, pver-top_lev+1
        1606             :           do j = 1, num_subcols
        1607             :             do i = 1, ngrdcol
        1608             :               if(state_sc%q(num_subcols*(i-1)+j,k,ixnumsnow) .lt. min_num_conc) then
        1609             :                  state_sc%q(num_subcols*(i-1)+j,k,ixnumsnow) = min_num_conc
        1610             :               end if
        1611             :             end do
        1612             :           end do
        1613             :         end do
        1614             :       endif
        1615             :
        1616             :       if ( l_outfld_subcol ) then
        1617             :
        1618             :         do k = 1, pver-top_lev+1
        1619             :           do i = 1, ngrdcol
        1620             :             do j = 1, num_subcols
        1621             :
        1622             :               ! Calc effective cloud fraction for testing
        1623             :               if ( ( lh_rc_clipped(i,j,pverp-k+1) .gt. qsmall ) &
        1624             :                      .or. ( X_nl_all_levs(i,j,pverp-k+1,iiPDF_ri) .gt. qsmall ) ) then
        1625             :                  eff_cldfrac(i,k) = eff_cldfrac(i,k) + lh_sample_point_weights(i,j,pverp-k+1)
        1626             :               else
        1627             :                 eff_cldfrac(i,k) = 0.0_r8
        1628             :               endif
        1629             :
        1630             :              end do
        1631             :
        1632             :              eff_cldfrac(i,k) = eff_cldfrac(i,k)/real(num_subcols, kind=r8)
        1633             :
        1634             :           end do
        1635             :         end do
        1636             :
        1637             :         ! Pack precip_frac for output
        1638             :         do k = 2, pverp-top_lev+1
        1639             :           do i = 1, ngrdcol
        1640             :             precip_frac_out(i,pver-k+2) = precip_fracs%precip_frac(i,k)
        1641             :           end do
        1642             :         end do
        1643             :
        1644             :         call outfld( 'SILHS_THLM_SCOL', THL_lh_out, pcols*psubcols, lchnk )
        1645             :         call outfld( 'SILHS_RT_SCOL', RT_lh_out, pcols*psubcols, lchnk )
        1646             :         call outfld( 'SILHS_OMEGA_SCOL', OMEGA_lh_out, pcols*psubcols, lchnk )
        1647             :         call outfld( 'SILHS_WM_SCOL', WM_lh_out, pcols*psubcols, lchnk )
        1648             :         call outfld( 'SILHS_RCM_SCOL', RCM_lh_out, pcols*psubcols, lchnk )
        1649             :         call outfld( 'SILHS_RICLD_SCOL', ICE_lh_out, pcols*psubcols, lchnk )
        1650             :         call outfld( 'SILHS_NICLD_SCOL', NICE_lh_out, pcols*psubcols, lchnk )
        1651             :         call outfld( 'SILHS_NCLD_SCOL', NCLW_lh_out, pcols*psubcols, lchnk )
        1652             :         call outfld( 'SILHS_RRAIN_SCOL', RAIN_lh_out, pcols*psubcols, lchnk )
        1653             :         call outfld( 'SILHS_NRAIN_SCOL', NRAIN_lh_out, pcols*psubcols, lchnk )
        1654             :         call outfld( 'SILHS_WEIGHT_SCOL', weights, pcols*psubcols, lchnk )
        1655             :         call outfld( 'NR_IN_LH', nrain, pcols, lchnk )
        1656             :         call outfld( 'SILHS_RTM', rtm, pcols, lchnk )
        1657             :         call outfld( 'SILHS_THLM', thlm, pcols, lchnk )
        1658             :         call outfld( 'SILHS_QC_IN', state%q(:,:,ixcldliq), pcols, lchnk )
        1659             :         call outfld( 'SILHS_QI_IN', state%q(:,:,ixcldice), pcols, lchnk )
        1660             :         call outfld( 'SILHS_NC_IN', state%q(:,:,ixnumliq), pcols, lchnk )
        1661             :         call outfld( 'SILHS_NI_IN', state%q(:,:,ixnumice), pcols, lchnk )
        1662             :         if ( l_est_kessler_microphys ) then
        1663             :           call outfld( 'AKM_CLUBB', AKm_out, pcols, lchnk )
        1664             :           call outfld( 'AKM_LH_CLUBB', lh_AKm_out, pcols, lchnk )
        1665             :         end if
        1666             :         call outfld( 'INVS_EXNER', invs_exner, pcols, lchnk )
        1667             :         call outfld( 'SILHS_ZTODT', ztodt_ptr, pcols, lchnk )
        1668             :         if ( subcol_SILHS_constrainmn ) then
        1669             :           call outfld( 'SILHS_MSC_CLDICE', meansc_ice, pcols, lchnk )
        1670             :           call outfld( 'SILHS_STDSC_CLDICE', stdsc_ice, pcols, lchnk )
        1671             :           if ( ixsnow > 0 ) then
        1672             :             call outfld( 'SILHS_MSC_CLDLIQ', meansc_liq, pcols, lchnk )
        1673             :             call outfld( 'SILHS_STDSC_CLDLIQ', stdsc_liq, pcols, lchnk )
        1674             :             call outfld( 'SILHS_MSC_Q', meansc_vap, pcols, lchnk )
        1675             :             call outfld( 'SILHS_STDSC_Q', stdsc_vap, pcols, lchnk )
        1676             :           endif ! ixsnow > 0
        1677             :         endif ! subcol_SILHS_constrainmn
        1678             :         call outfld( 'SILHS_EFF_CLDFRAC', eff_cldfrac, pcols, lchnk )
        1679             :         call outfld( 'SILHS_CLUBB_PRECIP_FRAC', precip_frac_out, pcols, lchnk )
        1680             :         call outfld( 'SILHS_CLUBB_ICE_SS_FRAC', ice_supersat_frac, pcols, lchnk )
        1681             :       end if
        1682             :
        1683             :       !$acc end data
        1684             :       !$acc wait
        1685             :
        1686             : #endif
        1687             : #endif
        1688           0 :    end subroutine subcol_gen_SILHS
        1689             :
        1690           0 :    subroutine subcol_ptend_avg_SILHS(ptend_sc, ngrdcol, lchnk, ptend)
        1691           0 :       use physics_buffer,   only: physics_buffer_desc
        1692             :       use subcol_utils,     only: subcol_ptend_get_firstsubcol, subcol_ptend_avg_shr, &
        1693             :                                   subcol_get_weight, subcol_get_filter, &
        1694             :                                   is_filter_set, is_weight_set
        1695             :
        1696             :       !-----------------------------------
        1697             :       ! Average the subcolumns dimension (pcols*psubcols) to the grid dimension (pcols)
        1698             :       !-----------------------------------
        1699             :
        1700             :       type(physics_ptend), intent(in)             :: ptend_sc        ! intent in
        1701             :       integer,  intent(in)                        :: ngrdcol       ! # grid cols
        1702             :       integer,  intent(in)                        :: lchnk         ! chunk index
        1703             :       type(physics_ptend), intent(inout)          :: ptend
        1704             :       ! Because we can't get a state passed in here, we might have to use values from the
        1705             :       ! subcolumn generation. This would make any conservation checks invalid if this
        1706             :       ! function is called after another parameterization... hmm.
        1707             :
        1708           0 :        call subcol_ptend_avg_shr(ptend_sc, ngrdcol, lchnk, ptend, is_filter_set(), is_weight_set())
        1709             :
        1710           0 :    end subroutine subcol_ptend_avg_SILHS
        1711             :
        1712           0 :    subroutine subcol_SILHS_var_covar_driver &
        1713             :               ( ztodt, state_sc, ptend_sc, &
        1714             :                 pbuf )
        1715             :
        1716             :      ! This subroutine calculates microphysical effects on five variances and
        1717             :      ! covariances: rtp2, thlp2, wprtp, wpthlp, and rtpthlp.
        1718             :      !
        1719             :      ! This code is experimental!!
        1720             :
        1721           0 :      use physics_buffer,          only: physics_buffer_desc, pbuf_get_index, pbuf_get_field
        1722             : #ifdef CLUBB_SGS
        1723             : #ifdef SILHS
        1724             :      use subcol_utils,            only: subcol_get_weight
        1725             :      use subcol_pack_mod,         only: subcol_unpack, subcol_get_nsubcol
        1726             :      use clubb_api_module,        only: T_in_K2thlm_api, &
        1727             :                                         init_pdf_params_api, &
        1728             :                                         copy_multi_pdf_params_to_single,&
        1729             :                                         pdf_parameter
        1730             :      use silhs_api_module,        only: lh_microphys_var_covar_driver_api
        1731             : #endif
        1732             : #endif
        1733             :
        1734             :      implicit none
        1735             :
        1736             :      ! Parameters
        1737             :      !  This fill value is set to catch errors; it should not be read.
        1738             :      real(r8), parameter                   :: fillvalue = -999._r8
        1739             :
        1740             :      ! Input Variables
        1741             :      real(r8), intent(in)                  :: ztodt        ! model time increment
        1742             :      type(physics_state), intent(in)       :: state_sc     ! state for sub-columns
        1743             :      type(physics_ptend), intent(in)       :: ptend_sc     ! ptend for sub-columns
        1744             :
        1745             :      ! Pointers
        1746             :      type(physics_buffer_desc), pointer    :: pbuf(:)
        1747             :
        1748             : #ifdef CLUBB_SGS
        1749             : #ifdef SILHS
        1750             :      ! Local Variables
        1751             :      integer :: lchnk, ngrdcol, igrdcol, isubcol, ns, k
        1752             :      integer, dimension(pcols) :: nsubcol
        1753             :      real(r8), dimension(pcols*psubcols)       :: weights_packed
        1754             :      real(r8), dimension(pcols,psubcols)       :: weights
        1755             :      real(r8), dimension(pcols,psubcols,pverp) :: rc_all, rv_all, rt_all, w_all, thl_all
        1756             :      real(r8), dimension(pcols,psubcols,pver ) :: s_all, t_all, zm_all, omega_all, pmid_all
        1757             :      real(r8), dimension(pcols,psubcols)       :: phis_all
        1758             :      real(r8), dimension(pcols,psubcols,pver ) :: stend, ttend
        1759             :      real(r8), dimension(pcols,psubcols,pverp) :: thltend, qctend, qvtend
        1760             :
        1761             :      real(r8), dimension(pcols,psubcols,pver)  :: dz_g, pdel_all, rho
        1762             :      real(r8), dimension(pcols,psubcols,pverp) :: zi_all
        1763             :
        1764             :      real(r8), dimension(pcols,psubcols,pver ) :: exner
        1765             :
        1766             :      ! Inputs to lh_microphys_var_covar_driver
        1767             :      real(r8), dimension(pcols,psubcols,pverp) :: rt_all_clubb, thl_all_clubb, w_all_clubb, &
        1768             :                                                   qctend_clubb, qvtend_clubb, thltend_clubb
        1769             :      real(r8), dimension(pcols,psubcols,pverp-top_lev+1) :: height_depndt_weights
        1770             :
        1771             :      ! Outputs from lh_microphys_var_covar_driver
        1772             :      real(r8), dimension(:,:), pointer :: rtp2_mc_zt, thlp2_mc_zt, wprtp_mc_zt, &
        1773             :                                           wpthlp_mc_zt, rtpthlp_mc_zt
        1774             :
        1775             :      ! pbuf indices
        1776             :      integer :: &
        1777             :        rtp2_mc_zt_idx,    &
        1778             :        thlp2_mc_zt_idx,   &
        1779             :        wprtp_mc_zt_idx,   &
        1780             :        wpthlp_mc_zt_idx,  &
        1781             :        rtpthlp_mc_zt_idx
        1782             :
        1783             :      type(pdf_parameter) :: pdf_params_single_col
        1784             :
        1785             :      !----- Begin Code -----
        1786             :
        1787             :      call init_pdf_params_api( pverp+1-top_lev, 1, pdf_params_single_col )
        1788             :
        1789             :      ! Don't do anything if this option isn't enabled.
        1790             :      if ( .not. subcol_SILHS_var_covar_src ) return
        1791             :
        1792             :      lchnk = state_sc%lchnk
        1793             :      ngrdcol  = state_sc%ngrdcol
        1794             :
        1795             :      ! Obtain indices
        1796             :      rtp2_mc_zt_idx = pbuf_get_index('rtp2_mc_zt')
        1797             :      thlp2_mc_zt_idx = pbuf_get_index('thlp2_mc_zt')
        1798             :      wprtp_mc_zt_idx = pbuf_get_index('wprtp_mc_zt')
        1799             :      wpthlp_mc_zt_idx = pbuf_get_index('wpthlp_mc_zt')
        1800             :      rtpthlp_mc_zt_idx = pbuf_get_index('rtpthlp_mc_zt')
        1801             :
        1802             :      ! Obtain pbuf fields for output
        1803             :      call pbuf_get_field(pbuf, rtp2_mc_zt_idx, rtp2_mc_zt)
        1804             :      call pbuf_get_field(pbuf, thlp2_mc_zt_idx, thlp2_mc_zt)
        1805             :      call pbuf_get_field(pbuf, wprtp_mc_zt_idx, wprtp_mc_zt)
        1806             :      call pbuf_get_field(pbuf, wpthlp_mc_zt_idx, wpthlp_mc_zt)
        1807             :      call pbuf_get_field(pbuf, rtpthlp_mc_zt_idx, rtpthlp_mc_zt)
        1808             :
        1809             :      ! Unpack needed tendencies from subcolumn ptends
        1810             :      call subcol_unpack(lchnk, ptend_sc%s(:,:), stend, fillvalue)
        1811             :      call subcol_unpack(lchnk, ptend_sc%q(:,:,ixcldliq), qctend(:,:,1:pver), fillvalue)
        1812             :      call subcol_unpack(lchnk, ptend_sc%q(:,:,ixq), qvtend(:,:,1:pver), fillvalue)
        1813             :
        1814             :      ! Unpack sample point values from subcolumn states
        1815             :      call subcol_unpack(lchnk, state_sc%q(:,:,ixcldliq), rc_all(:,:,1:pver), fillvalue)
        1816             :      call subcol_unpack(lchnk, state_sc%q(:,:,ixq),      rv_all(:,:,1:pver), fillvalue)
        1817             :      call subcol_unpack(lchnk, state_sc%omega(:,:),      omega_all (:,:,:),  fillvalue)
        1818             :      call subcol_unpack(lchnk, state_sc%s(:,:),          s_all,              fillvalue)
        1819             :      call subcol_unpack(lchnk, state_sc%zm,              zm_all,             fillvalue)
        1820             :      call subcol_unpack(lchnk, state_sc%phis,            phis_all,           fillvalue)
        1821             :      call subcol_unpack(lchnk, state_sc%zi,              zi_all,             fillvalue)
        1822             :      call subcol_unpack(lchnk, state_sc%pdel,            pdel_all,           fillvalue)
        1823             :      call subcol_unpack(lchnk, state_sc%pmid,            pmid_all,           fillvalue)
        1824             :
        1825             :      ! Initialize fields to fillvalue.
        1826             :      rt_all  = fillvalue
        1827             :      thl_all = fillvalue
        1828             :      w_all   = fillvalue
        1829             :      qctend  = fillvalue
        1830             :      qvtend  = fillvalue
        1831             :      thltend = fillvalue
        1832             :
        1833             :      ! How many subcolumns in each column?
        1834             :      call subcol_get_nsubcol(lchnk, nsubcol)
        1835             :
        1836             :      do igrdcol = 1, ngrdcol
        1837             :         do isubcol = 1, nsubcol(igrdcol)
        1838             :
        1839             :            rt_all(igrdcol,isubcol,top_lev:pver) = rc_all(igrdcol,isubcol,top_lev:pver) &
        1840             :                                                   + rv_all(igrdcol,isubcol,top_lev:pver)
        1841             :
        1842             :            ! Compute dry static density on CLUBB vertical grid
        1843             :            do k = top_lev, pver
        1844             :               dz_g(igrdcol,isubcol,k) = zi_all(igrdcol,isubcol,k) - zi_all(igrdcol,isubcol,k+1) ! thickness
        1845             :               rho(igrdcol,isubcol,k) = (rga)*pdel_all(igrdcol,isubcol,k)/dz_g(igrdcol,isubcol,k)
        1846             :            end do
        1847             :
        1848             :            ! Compute w from omega
        1849             :            w_all(igrdcol,isubcol,top_lev:pver) = -omega_all(igrdcol,isubcol,top_lev:pver) &
        1850             :                                                   / ( rho(igrdcol,isubcol,top_lev:pver) * gravit )
        1851             :
        1852             :            ! Convert stend and s_all to ttend and t_all
        1853             :            !  Note 1: With subcolumns, cpair is truly a constant (I think).
        1854             :            !  Note 2: For tendencies, the extra terns zm and phis should
        1855             :            !          not be included in the calculation.
        1856             :            ttend(igrdcol,isubcol,top_lev:pver) = stend(igrdcol,isubcol,top_lev:pver) / cpair
        1857             :
        1858             :            do k = top_lev, pver
        1859             :               t_all(igrdcol,isubcol,k) = ( s_all(igrdcol,isubcol,k) &
        1860             :                                            - gravit * zm_all(igrdcol,isubcol,k) &
        1861             :                                            - phis_all(igrdcol,isubcol) ) / cpair
        1862             :            end do ! k = 1, pver
        1863             :
        1864             :            ! This formula is taken from earlier in this file.
        1865             :            exner(igrdcol,isubcol,top_lev:pver) &
        1866             :            = ( pmid_all(igrdcol,isubcol,top_lev:pver) / p0_clubb )**(rair/cpair)
        1867             :
        1868             :            ! Note: all tendencies or all means should be used in the call to
        1869             :            !       T_in_K2thlm_api (with the exception of exner)
        1870             :            do k = top_lev, pver
        1871             :               thltend(igrdcol,isubcol,k) &
        1872             :               = T_in_K2thlm_api( ttend(igrdcol,isubcol,k), exner(igrdcol,isubcol,k), &
        1873             :                                  qctend(igrdcol,isubcol,k) )
        1874             :               thl_all(igrdcol,isubcol,k) &
        1875             :               = T_in_K2thlm_api( t_all(igrdcol,isubcol,k), exner(igrdcol,isubcol,k), &
        1876             :                                  rc_all(igrdcol,isubcol,k) )
        1877             :            end do ! k = 1, pver
        1878             :
        1879             :            ! Add ghost points
        1880             :            rt_all (igrdcol,isubcol,pverp) = rt_all (igrdcol,isubcol,pver)
        1881             :            thl_all(igrdcol,isubcol,pverp) = thl_all(igrdcol,isubcol,pver)
        1882             :            w_all  (igrdcol,isubcol,pverp) = w_all  (igrdcol,isubcol,pver)
        1883             :            qctend (igrdcol,isubcol,pverp) = qctend (igrdcol,isubcol,pver)
        1884             :            qvtend (igrdcol,isubcol,pverp) = qvtend (igrdcol,isubcol,pver)
        1885             :            thltend(igrdcol,isubcol,pverp) = thltend(igrdcol,isubcol,pver)
        1886             :
        1887             :            ! Flip inputs to CLUBB's grid. Note the dimension ordering change.
        1888             :            rt_all_clubb(igrdcol,isubcol,1:pverp) = clubb_flip_grid( rt_all(igrdcol,isubcol,1:pverp) )
        1889             :            thl_all_clubb(igrdcol,isubcol,1:pverp) = clubb_flip_grid( thl_all(igrdcol,isubcol,1:pverp) )
        1890             :            w_all_clubb(igrdcol,isubcol,1:pverp) = clubb_flip_grid( w_all(igrdcol,isubcol,1:pverp) )
        1891             :            qctend_clubb(igrdcol,isubcol,1:pverp) = clubb_flip_grid( qctend(igrdcol,isubcol,1:pverp) )
        1892             :            qvtend_clubb(igrdcol,isubcol,1:pverp) = clubb_flip_grid( qvtend(igrdcol,isubcol,1:pverp) )
        1893             :            thltend_clubb(igrdcol,isubcol,1:pverp) = clubb_flip_grid( thltend(igrdcol,isubcol,1:pverp) )
        1894             :
        1895             :         end do ! isubcol = 1, nsubcol(igrdcol)
        1896             :      end do ! igrdcol = 1, ngrdcol
        1897             :
        1898             :      ! Obtain weights
        1899             :      call subcol_get_weight(lchnk, weights_packed)
        1900             :      call subcol_unpack(lchnk, weights_packed, weights, fillvalue)
        1901             :
        1902             :      ! Call lh_microphys_var_covar_driver for each column
        1903             :      do igrdcol=1, ngrdcol
        1904             :        ns = nsubcol(igrdcol)
        1905             :
        1906             :        ! This code assumes that the weights are height independent.
        1907             :        ! It will have to change once the weights vary with altitude!
        1908             :        ! I'm not sure whether the grid will need to be flipped.
        1909             :        do k = 1, pverp-top_lev+1
        1910             :           height_depndt_weights(igrdcol,1:ns,k) = weights(igrdcol,1:ns)
        1911             :        end do
        1912             :
        1913             :        ! Copy the igrdcol column from the multicolumn pdf_params_chnk to the single column
        1914             :        ! version of pdf_params_single_col since lh_microphys_var_covar_driver_api only
        1915             :        ! works over 1 column currently
        1916             :        call copy_multi_pdf_params_to_single( pdf_params_chnk(lchnk), igrdcol, &
        1917             :                                              pdf_params_single_col )
        1918             :
        1919             :        ! Make the call!!!!!
        1920             :        call lh_microphys_var_covar_driver_api &
        1921             :             ( pverp-top_lev+1, ns, ztodt, height_depndt_weights(igrdcol,1:ns,1:pverp-top_lev+1), &
        1922             :               pdf_params_single_col, &
        1923             :               rt_all_clubb(igrdcol,1:ns,1:pverp-top_lev+1), thl_all_clubb(igrdcol,1:ns,1:pverp-top_lev+1), &
        1924             :               w_all_clubb(igrdcol,1:ns,1:pverp-top_lev+1), qctend_clubb(igrdcol,1:ns,1:pverp-top_lev+1), &
        1925             :               qvtend_clubb(igrdcol,1:ns,1:pverp-top_lev+1), thltend_clubb(igrdcol,1:ns,1:pverp-top_lev+1), &
        1926             :               silhs_config_flags%l_lh_instant_var_covar_src, &
        1927             :               rtp2_mc_zt(igrdcol,1:pverp-top_lev+1), thlp2_mc_zt(igrdcol,1:pverp-top_lev+1), &
        1928             :               wprtp_mc_zt(igrdcol,1:pverp-top_lev+1), wpthlp_mc_zt(igrdcol,1:pverp-top_lev+1), &
        1929             :               rtpthlp_mc_zt(igrdcol,1:pverp-top_lev+1) )
        1930             :
        1931             :        ! The *_mc_zt microphysics tendencies are passed out of SILHS and back
        1932             :        ! to CLUBB without being used at all in the rest of the host model code.
        1933             :        ! The arrays aren't flipped for the *_mc_zt microphysics tendencies, and
        1934             :        ! they don't need to be.
        1935             :
        1936             :        ! CLUBB used pverp vertical levels, but SILHS only uses
        1937             :        ! pverp - top_lev + 1 vertical levels.
        1938             :        ! Fill the upper levels with 0s when necessary.
        1939             :        if ( pverp > pverp-top_lev+1 ) then
        1940             :           rtp2_mc_zt(igrdcol,pverp-top_lev+2:pverp) = 0.0_r8
        1941             :           thlp2_mc_zt(igrdcol,pverp-top_lev+2:pverp) = 0.0_r8
        1942             :           wprtp_mc_zt(igrdcol,pverp-top_lev+2:pverp) = 0.0_r8
        1943             :           wpthlp_mc_zt(igrdcol,pverp-top_lev+2:pverp) = 0.0_r8
        1944             :           rtpthlp_mc_zt(igrdcol,pverp-top_lev+2:pverp) = 0.0_r8
        1945             :        endif ! pverp > pverp-top_lev+1
        1946             :
        1947             :      end do ! igrdcol = 1, ngrdcol
        1948             : #endif
        1949             : #endif
        1950             :
        1951           0 :      return
        1952           0 :    end subroutine subcol_SILHS_var_covar_driver
        1953             : #ifdef SILHS
        1954             :    real(r8) function meansc(arr_in, w_in, ns) result(val)
        1955             :       real(r8), intent(in) :: ns                         ! Length of Array
        1956             :       real(r8), dimension(int(ns)), intent(in) :: arr_in      ! Input array
        1957             :       real(r8), dimension(int(ns)), intent(in) :: w_in        ! Weights
        1958             :       real(r8) :: acc  ! accumulator
        1959             :       integer :: i
        1960             :       acc = 0
        1961             :       val = 0
        1962             :       do i=1,ns
        1963             :          acc = acc + arr_in(i)*w_in(i)
        1964             :       end do
        1965             :       val = acc/ns
        1966             :    end function
        1967             :
        1968             :    real(r8) function stdsc(arr_in, w_in, mn_in, ns) result(val)
        1969             :       real(r8), intent(in) :: ns  ! Number of elements (subcolumns)
        1970             :       real(r8), dimension(int(ns)), intent(in) :: arr_in, w_in  !Input array and weights
        1971             :       real(r8), intent(in) :: mn_in   ! The mean of arr_in
        1972             :       real(r8) :: accvar, var
        1973             :       integer :: i
        1974             :       accvar = 0
        1975             :       do i=1,ns
        1976             :          accvar = accvar + ((arr_in(i)-mn_in)**2)*w_in(i)
        1977             :       end do
        1978             :       var = accvar/ns
        1979             :       val = sqrt(var)
        1980             :    end function
        1981             :
        1982             :    subroutine THL_profile(nz, ABST_prof, ex_prof, rcm_prof, THL_prof)
        1983             :
        1984             :       use clubb_api_module,              only : T_in_K2thlm_api
        1985             :
        1986             :       integer,                 intent(in)  :: nz         ! Num vert levels
        1987             :       real(r8), dimension(nz), intent(in)  :: ABST_prof  ! Abs Temp prof
        1988             :       real(r8), dimension(nz), intent(in)  :: ex_prof    ! Profile of Exner func
        1989             :       real(r8), dimension(nz), intent(in)  :: rcm_prof   ! Profile of Cld Wat MR
        1990             :       real(r8), dimension(nz), intent(out) :: THL_prof  ! LWPT prof
        1991             :       integer :: i
        1992             :
        1993             :       do i=1,nz
        1994             :          THL_prof(i) = T_in_K2thlm_api(ABST_prof(i), ex_prof(i), rcm_prof(i))
        1995             :       end do
        1996             :
        1997             :    end subroutine
        1998             :
        1999             :    subroutine subcol_constrainmn( num_subcols, samples, weights, grid_mean, mean_sc, std_sc )
        2000             :
        2001             :       ! Input/Output Variables
        2002             :       integer, intent(in) :: num_subcols
        2003             :       real(r8), dimension(num_subcols, pverp), intent(inout) :: samples
        2004             :       real(r8), dimension(num_subcols), intent(in) :: weights
        2005             :       real(r8), dimension(pverp), intent(in) :: grid_mean
        2006             :       real(r8), dimension(pver), intent(out), optional :: mean_sc, std_sc
        2007             :
        2008             :       ! Local Variables
        2009             :       real(r8) :: meansc_loc, adj_rat
        2010             :       integer :: k
        2011             :    !------------------------------------------------------------------
        2012             :       !----- Begin Code -----
        2013             :       do k=1, pver
        2014             :          meansc_loc = meansc( samples(:,k), weights(:), real(num_subcols, r8) )
        2015             :
        2016             :          if (present(mean_sc)) &
        2017             :             mean_sc(k) = meansc_loc
        2018             :          if (present(std_sc)) &
        2019             :             std_sc(k) = stdsc( samples(:,k), weights(:), meansc_loc, &
        2020             :                                real(num_subcols, r8) )
        2021             :
        2022             :          if ( meansc_loc > 0.0_r8 ) then
        2023             :             adj_rat = grid_mean(k)/meansc_loc
        2024             :          else
        2025             :             ! If the mean is zero, then zero out all subcolumns to avoid
        2026             :             ! negative samples
        2027             :             adj_rat = 0.0_r8
        2028             :          end if
        2029             :          samples(:,k) = samples(:,k) * adj_rat
        2030             :       end do
        2031             :    end subroutine subcol_constrainmn
        2032             :
        2033             :    ! =============================================================================== !
        2034             :    !                                                                                 !
        2035             :    ! =============================================================================== !
        2036             :    function clubb_flip_grid ( profile ) result( profile_flipped )
        2037             :
        2038             :      ! Description:
        2039             :      !   Swaps the elements in profile so they are in reverse order. CAM and
        2040             :      !   CLUBB's grids are flipped with respect to each other.
        2041             :      !
        2042             :      !   Usage:
        2043             :      !     clubb_var = clubb_flip_grid( cam_var )
        2044             :      !     cam_var   = clubb_flip_grid( clubb_var )
        2045             :
        2046             :      implicit none
        2047             :
        2048             :      ! Input Variable
        2049             :      real(r8), dimension(pverp), intent(in) :: profile
        2050             :
        2051             :      ! Output Variable
        2052             :      real(r8), dimension(pverp) :: profile_flipped
        2053             :
        2054             :      ! Local Variable
        2055             :      integer :: k
        2056             :
        2057             :      do k=1, pverp
        2058             :        profile_flipped(k) = profile(pverp-k+1)
        2059             :      end do ! k=1, pverp
        2060             :
        2061             :      return
        2062             :    end function clubb_flip_grid
        2063             :    ! =============================================================================== !
        2064             :    !                                                                                 !
        2065             :    ! =============================================================================== !
        2066             : #endif
        2067             :    !============================================================================
        2068           0 :    subroutine subcol_SILHS_fill_holes_conserv( state, dt, ptend, pbuf )
        2069             :
        2070             :      ! The William F. Buckley Jr. Conservative Hole Filler.
        2071             :
        2072             :      ! Description:
        2073             :      ! Stops holes from forming in a hydrometeor mixing ratio by reducing the
        2074             :      ! microphysics tendency of that hydrometeor mixing ratio which would
        2075             :      ! otherwise cause that hydrometeor mixing ratio to have a negative value
        2076             :      ! once the microphysics tendency is applied.  This code is used to prevent
        2077             :      ! holes in water mass, not number concentration.
        2078             :      !
        2079             :      ! This subroutine is called after microphysics has completed and after
        2080             :      ! microphysics fields from subcolumns have been averaged back to grid
        2081             :      ! columns, but before the grid-column microphysics tendencies have been
        2082             :      ! applied in physics_update.  This code is meant for use with the SILHS
        2083             :      ! subcolumn approach.  This code needs to be applied to grid columns, not
        2084             :      ! subcolumns.
        2085             :      !
        2086             :      ! This code adjusts the tendencies (ptend) before they are used to update
        2087             :      ! the grid mean fields (state variables).
        2088             :      !
        2089             :      ! The column-integrated total water needs to be conserved during
        2090             :      ! microphysics.  The conserved amount includes the amount of water that
        2091             :      ! precipitated to the ground from sedimentation during microphysics.
        2092             :      ! The conservation equation for each grid column is:
        2093             :      !
        2094             :      ! SUM(k=top_lev:pver) ( rv_start(k) + rc_start(k) + rr_start(k)
        2095             :      !                       + ri_start(k) + rs_start(k) ) * pdel(k) / g
        2096             :      ! = SUM(k=top_lev:pver) ( rv(k) + rc(k) + rr(k) + ri(k) + rs(k) )
        2097             :      !                       * pdel(k) / g
        2098             :      !   + prect * dt * 1000;
        2099             :      !
        2100             :      ! where rv_start, rc_start, rr_start, ri_start, and rs_start are water
        2101             :      ! vapor, cloud water, rain water, cloud ice, and snow mixing ratios before
        2102             :      ! microphysics is called; rv, rc, rr, ri, and rs are water vapor, cloud
        2103             :      ! water, rain water, cloud ice, and snow mixing ratios after being updated
        2104             :      ! by microphysics; pdel is the pressure difference between vertical levels,
        2105             :      ! g is gravity, and prect * dt * 1000 is the total amount of water (from
        2106             :      ! all precipitating hydrometeors) that sedimented to the ground during
        2107             :      ! microphysics (dt is the timestep used for microphysics).  The units of
        2108             :      ! column-integrated total water are kg (water) / m^2.
        2109             :      !
        2110             :      ! All the updated hydrometeor fields are related to the hydrometeor fields
        2111             :      ! at the start by:
        2112             :      !
        2113             :      ! rv(k) = rv_start(k) + rv_tend(k) * dt;
        2114             :      ! rc(k) = rc_start(k) + rc_tend(k) * dt;
        2115             :      ! rr(k) = rr_start(k) + rr_tend(k) * dt;
        2116             :      ! ri(k) = ri_start(k) + ri_tend(k) * dt; and
        2117             :      ! rs(k) = rs_start(k) + rs_tend(k) * dt;
        2118             :      !
        2119             :      ! where rv_tend, rc_tend, rr_tend, ri_tend, and rs_tend are water vapor,
        2120             :      ! cloud water, rain water, cloud ice, and snow mixing ratio tendencies
        2121             :      ! from microphysics, which includes the sum of microphysics process rates
        2122             :      ! and sedimentation.  When these equations are applied to the equation
        2123             :      ! for column-integrated total water, that equation becomes:
        2124             :      !
        2125             :      ! SUM(k=top_lev:pver) ( rv_tend(k) + rc_tend(k) + rr_tend(k)
        2126             :      !                       + ri_tend(k) + rs_tend(k) ) * dt * pdel(k) / g
        2127             :      ! + prect * dt * 1000 = 0.
        2128             :      !
        2129             :      ! As stated above, the hydrometeor tendencies are the sum of tendencies
        2130             :      ! from microphysics process rates and tendencies from sedimentation:
        2131             :      !
        2132             :      ! rv_tend(k) = rv_mc_tend(k);
        2133             :      ! rc_tend(k) = rc_mc_tend(k) + rc_sed_tend(k);
        2134             :      ! rr_tend(k) = rr_mc_tend(k) + rr_sed_tend(k);
        2135             :      ! ri_tend(k) = ri_mc_tend(k) + ri_sed_tend(k); and
        2136             :      ! rs_tend(k) = rs_mc_tend(k) + rs_sed_tend(k);
        2137             :      !
        2138             :      ! where rv_mc_tend, rc_mc_tend, rr_mc_tend, ri_mc_tend, and rs_mc_tend are
        2139             :      ! the tendencies of water vapor, cloud water, rain water, cloud ice, and
        2140             :      ! snow from microphysics process rates, and rc_sed_tend, rr_sed_tend,
        2141             :      ! ri_sed_tend, and rs_sed_tend are the tendencies of cloud water,
        2142             :      ! rain water, cloud ice, and snow from sedimentation.  When these equations
        2143             :      ! are applied to the equation for column-integrated total water, that
        2144             :      ! equation becomes:
        2145             :      !
        2146             :      ! SUM(k=top_lev:pver) ( rv_mc_tend(k) + rc_mc_tend(k) + rr_mc_tend(k)
        2147             :      !                       + ri_mc_tend(k) + rs_mc_tend(k) )
        2148             :      !                     * dt * pdel(k) / g
        2149             :      ! + SUM(k=top_lev:pver) ( rc_sed_tend(k) + rr_sed_tend(k) + ri_sed_tend(k)
        2150             :      !                         + rs_sed_tend(k) ) * dt * pdel(k) / g
        2151             :      ! + prect * dt * 1000 = 0.
        2152             :      !
        2153             :      ! At any vertical level, the tendencies from microphysics process rates
        2154             :      ! (mc_tend variables) must balance:
        2155             :      !
        2156             :      ! rv_mc_tend(k) + rc_mc_tend(k) + rr_mc_tend(k)
        2157             :      ! + ri_mc_tend(k) + rs_mc_tend(k) = 0; for all k from top_lev to pver.
        2158             :      !
        2159             :      ! The column-integrated total water equation can be applied to
        2160             :      ! sedimentation:
        2161             :      !
        2162             :      ! SUM(k=top_lev:pver) ( rc_sed_tend(k) + rr_sed_tend(k) + ri_sed_tend(k)
        2163             :      !                       + rs_sed_tend(k) ) * dt * pdel(k) / g
        2164             :      ! + prect * dt * 1000 = 0.
        2165             :      !
        2166             :      ! The total precipitation rate, prect, can be split into liquid
        2167             :      ! precipitation rate, precl, and frozen precipitation rate, preci:
        2168             :      !
        2169             :      ! prect = precl + preci.
        2170             :      !
        2171             :      ! The microphysics code outputs prect and preci, so precl can be calculated
        2172             :      ! by precl = prect - preci.  The column-integrated total water equation can
        2173             :      ! be split into:
        2174             :      !
        2175             :      ! SUM(k=top_lev:pver) ( rc_sed_tend(k) + rr_sed_tend(k) )
        2176             :      !                     * dt * pdel(k) / g
        2177             :      ! + precl * dt * 1000 = 0; and
        2178             :      !
        2179             :      ! SUM(k=top_lev:pver) ( ri_sed_tend(k) + rs_sed_tend(k) )
        2180             :      !                     * dt * pdel(k) / g
        2181             :      ! + preci * dt * 1000 = 0.
        2182             :      !
        2183             :      ! Overall, the conservation methods used in this subroutine are:
        2184             :      !
        2185             :      ! 1) When adjusting the tendencies from microphysics process rates,
        2186             :      !    conserve:
        2187             :      !
        2188             :      !    rv_mc_tend(k) + rc_mc_tend(k) + rr_mc_tend(k)
        2189             :      !    + ri_mc_tend(k) + rs_mc_tend(k) = 0; for all k from top_lev to pver.
        2190             :      !
        2191             :      ! 2) When adjusting the tendencies from microphysics process rates, adjust
        2192             :      !    dry static energy appropriately.  The change in dry static energy
        2193             :      !    is necessary because of phase changes.  This "puts back" the extra dry
        2194             :      !    static energy that was "taken out" when an excessive phase-changing
        2195             :      !    process rate was produced by microphysics.
        2196             :      !
        2197             :      ! 3) When adjusting the hydrometeor tendency from sedimentation of a
        2198             :      !    liquid hydrometeor (cloud water or rain water), conserve:
        2199             :      !
        2200             :      !    SUM(k=top_lev:pver) ( rc_sed_tend(k) + rr_sed_tend(k) )
        2201             :      !                        * dt * pdel(k) / g
        2202             :      !    + precl * dt * 1000 = 0.
        2203             :      !
        2204             :      ! 4) When adjusting the hydrometeor tendency from sedimentation of a
        2205             :      !    frozen hydrometeor (cloud ice or snow), conserve:
        2206             :      !
        2207             :      !    SUM(k=top_lev:pver) ( ri_sed_tend(k) + rs_sed_tend(k) )
        2208             :      !                        * dt * pdel(k) / g
        2209             :      !    + preci * dt * 1000 = 0.
        2210             :      !
        2211             :      ! The conservative hole filler works as follows.  The total microphysics
        2212             :      ! tendency for each hydrometeor is provided in ptend.  This is the sum of
        2213             :      ! the microphysics process rate tendency and sedimentation tendency for
        2214             :      ! each hydrometeor.  The sedimentation tendency is provided in pbuf.  The
        2215             :      ! sedimentation tendency is subtracted off the total microphysics tendency
        2216             :      ! to produce the microphysics process rate tendency for each hydrometeor.
        2217             :      ! The microphysics process rate tendency is adjusted when necessary so that
        2218             :      ! holes in the hydrometeor are not produced by microphysics process rates.
        2219             :      ! When a hydrometeor's negative microphysics process rate tendency needs to
        2220             :      ! be made smaller in magnitude to avoid a hole, all hydrometeor tendencies
        2221             :      ! that are positive at that grid level are also decreased proportionately
        2222             :      ! to maintain a balance.  Dry static energy tendency is also adjusted
        2223             :      ! appropriately when necessary.  After this, the vertical integral of each
        2224             :      ! hydrometeor species is greater than or equal to 0.
        2225             :      !
        2226             :      ! The sedimentation tendency is then added back onto the new microphysics
        2227             :      ! process rate tendency to produce a new total microphysics tendency for
        2228             :      ! each hydrometeor.  Since the sedimentation tendency was based on the old
        2229             :      ! value of hydrometeor, before the hole-filling adjustment, it is possible
        2230             :      ! that the new total microphysics tendency may produce holes.  When this
        2231             :      ! happens, sedimentation hole filling fills holes in the vertical profile
        2232             :      ! of each hydrometeor.  Holes are filled using mass from other vertical
        2233             :      ! levels for the same hydrometeor (or from a same-phase hydrometeor when
        2234             :      ! necessary).  Since the vertical integral of sedimentation tendency
        2235             :      ! (including surface precipitation rate) is 0, the vertical integral of the
        2236             :      ! hydrometeor must be greater than or equal to 0, which means that all
        2237             :      ! holes can be filled.  The result is that all holes in any hydrometeor
        2238             :      ! mixing ratio are filled completely and conservatively.  The value of
        2239             :      ! ptend is updated appropriately so that it can be applied later in
        2240             :      ! physics_update.
        2241             :
        2242             :      !----------------------------------------------------------------------
        2243             :
        2244           0 :      use physics_buffer, only: &
        2245             :          physics_buffer_desc, &
        2246             :          pbuf_get_field
        2247             :
        2248             :      use ppgrid, only: &
        2249             :          pcols
        2250             :
        2251             :      use constituents, only: &
        2252             :          qmin
        2253             :
        2254             :      use ref_pres, only: &
        2255             :          top_lev => trop_cloud_top_lev
        2256             :
        2257             :      implicit none
        2258             :
        2259             :      ! Input Variables
        2260             :      type(physics_state), intent(in) :: state     ! Physics state variables
        2261             :      real(r8), intent(in) :: dt                   ! Time step duration
        2262             :
        2263             :      ! Input/Output Variables
        2264             :      type(physics_ptend),  intent(inout) :: ptend  ! Parameterization tendencies
        2265             :      type(physics_buffer_desc), pointer :: pbuf(:) ! Physics buffer
        2266             :
        2267             :      ! Local Variables
        2268             :      real(r8), dimension(pcols,pver) :: &
        2269             :        rv_start, & ! Water vapor mixing ratio at start of microphysics  [kg/kg]
        2270             :        rc_start, & ! Cloud water mixing ratio at start of microphysics  [kg/kg]
        2271             :        rr_start, & ! Rain water mixing ratio at start of microphysics   [kg/kg]
        2272             :        ri_start, & ! Cloud ice mixing ratio at start of microphysics    [kg/kg]
        2273             :        rs_start    ! Snow mixing ratio at start of microphysics         [kg/kg]
        2274             :
        2275             :      real(r8), dimension(pcols,pver) :: &
        2276             :        rv_tend, & ! Water vapor mixing ratio tendency  [kg/kg/s]
        2277             :        rc_tend, & ! Cloud water mixing ratio tendency  [kg/kg/s]
        2278             :        rr_tend, & ! Rain water mixing ratio tendency   [kg/kg/s]
        2279             :        ri_tend, & ! Cloud ice mixing ratio tendency    [kg/kg/s]
        2280             :        rs_tend, & ! Snow mixing ratio tendency         [kg/kg/s]
        2281             :        stend      ! Dry static energy tendency         [J/kg/s]
        2282             :
        2283             :      real(r8), dimension(:), pointer :: &
        2284           0 :        prect,    & ! Total microphysics precipitation rate (surface)      [m/s]
        2285           0 :        preci,    & ! Ice-phase microphysics precipitation rate (surface)  [m/s]
        2286           0 :        prec_str, & ! Total surface precipitation rate from stratoform     [m/s]
        2287           0 :        snow_str    ! Snow surface precipitation rate from stratoform      [m/s]
        2288             :
        2289             :      real(r8), dimension(:,:), pointer :: &
        2290           0 :        rc_sed_tend, & ! Mean cloud water sedimentation tendency        [kg/kg/s]
        2291           0 :        rr_sed_tend, & ! Mean rain water sedimentation tendency         [kg/kg/s]
        2292           0 :        ri_sed_tend, & ! Mean cloud ice sedimentation tendency          [kg/kg/s]
        2293           0 :        rs_sed_tend, & ! Mean snow sedimentation tendency               [kg/kg/s]
        2294           0 :        vtrmc,       & ! Mean cloud water sedimentation velocity        [m/s]
        2295           0 :        umr,         & ! Mean rain water sedimentation velocity         [m/s]
        2296           0 :        vtrmi,       & ! Mean cloud ice sedimentation velocity          [m/s]
        2297           0 :        ums,         & ! Mean snow sedimentation velocity               [m/s]
        2298           0 :        rc_sed_evap, & ! Mean evap of cloud water during sedimentation  [kg/kg/s]
        2299           0 :        ri_sed_subl    ! Mean subl of cloud ice during sedimentation    [kg/kg/s]
        2300             :
        2301             :      real(r8), dimension(pcols,pver) :: &
        2302             :        rv_mc_tend, & ! Water vapor mixing ratio microphysics tendency  [kg/kg/s]
        2303             :        rc_mc_tend, & ! Cloud water mixing ratio microphysics tendency  [kg/kg/s]
        2304             :        rr_mc_tend, & ! Rain water mixing ratio microphysics tendency   [kg/kg/s]
        2305             :        ri_mc_tend, & ! Cloud ice mixing ratio microphysics tendency    [kg/kg/s]
        2306             :        rs_mc_tend    ! Snow mixing ratio microphysics tendency         [kg/kg/s]
        2307             :
        2308             :      real(r8) :: &
        2309             :        rv_curr, & ! Current water vapor mixing ratio    [kg/kg]
        2310             :        rc_curr, & ! Current cloud water mixing ratio    [kg/kg]
        2311             :        rr_curr, & ! Current rain water mixing ratio     [kg/kg]
        2312             :        ri_curr, & ! Current cloud ice mixing ratio      [kg/kg]
        2313             :        rs_curr    ! Current snow mixing ratio           [kg/kg]
        2314             :
        2315             :      logical :: &
        2316             :        l_pos_rv_mc_tend, & ! Flag for positive water vapor mixing ratio mc tend.
        2317             :        l_pos_rc_mc_tend, & ! Flag for positive cloud water mixing ratio mc tend.
        2318             :        l_pos_rr_mc_tend, & ! Flag for positive rain water mixing ratio mc tend.
        2319             :        l_pos_ri_mc_tend, & ! Flag for positive cloud ice mixing ratio mc tend.
        2320             :        l_pos_rs_mc_tend    ! Flag for positive snow mixing ratio mc tend.
        2321             :
        2322             :      real(r8) :: &
        2323             :        mc_tend_max_mag,     & ! Max. allowable mag. of (neg.) mc tend [kg/kg/s]
        2324             :        mc_tend_correction,  & ! Amnt. correction necessary to mc tend [kg/kg/s]
        2325             :        total_mc_positive,   & ! Total of all positive mc tendencies   [kg/kg/s]
        2326             :        mc_correction_ratio    ! Ratio: mc_tend_correction/total_mc_positive [-]
        2327             :
        2328             :      real(r8), dimension(pcols) :: &
        2329             :        precl    ! Liquid-phase precipitation rate (surface)        [m/s]
        2330             :
        2331             :      ! Budgeting terms for hole filling.
        2332             :      ! These variables are for use in stats output.
        2333             :      real(r8), dimension(pcols,pver) :: &
        2334             :        rv_hf_tend, & ! Water vapor mixing ratio hole-filling tendency  [kg/kg/s]
        2335             :        rc_hf_tend, & ! Cloud water mixing ratio hole-filling tendency  [kg/kg/s]
        2336             :        rr_hf_tend, & ! Rain water mixing ratio hole-filling tendency   [kg/kg/s]
        2337             :        ri_hf_tend, & ! Cloud ice mixing ratio hole-filling tendency    [kg/kg/s]
        2338             :        rs_hf_tend, & ! Snow mixing ratio hole-filling tendency         [kg/kg/s]
        2339             :        s_hf_tend     ! Dry static energy hole-filling tendency         [J/kg/s]
        2340             :
        2341             :      integer :: ncol  ! Number of grid columns
        2342             :
        2343             :      integer :: icol, k  ! Loop indices
        2344             :
        2345             :      ! Flag to perform hole filling after the original sedimentation tendency
        2346             :      ! is added back on to the new microphysics process tendency.  This calls
        2347             :      ! the sedimentation hole filler.
        2348             :      logical, parameter :: &
        2349             :        l_sed_hole_fill = .true.
        2350             :
        2351             :      logical, parameter :: &
        2352             :        l_check_conservation = .true. ! Flag to perform water conservation check
        2353             :
        2354             :      ! Vertically-integrated grand total water (rv+rc+rr+ri+rs)    [kg/m^2]
        2355             :      real(r8), dimension(pcols) :: &
        2356             :        grand_total_water_column_start,  & ! Column integral at start
        2357             :        grand_total_water_column_finish    ! Column integral at finish
        2358             :
        2359             :      ! Vertically-integrated total water energy    [J/m^2]
        2360             :      real(r8), dimension(pcols) :: &
        2361             :        total_energy_column_start,  & ! Column integral at start
        2362             :        total_energy_column_finish    ! Column integral at finish
        2363             :
        2364             :      real(r8), dimension(pcols) :: &
        2365             :        tot_water_rel_err,  & ! Relative error: vert-integrated grand total water
        2366             :        tot_energy_rel_err    ! Relative error: vert-integrated total energy
        2367             :
        2368             :      real(r8), parameter :: &
        2369             :        err_thresh = 1.0e-14_r8  ! Threshold of relative error
        2370             :
        2371             :
        2372             :      ! Get the number of grid columns.
        2373           0 :      ncol = state%ncol
        2374             :
        2375             :      ! Get fields from the pbuf.
        2376           0 :      call pbuf_get_field(pbuf, prec_pcw_idx, prect)
        2377           0 :      call pbuf_get_field(pbuf, snow_pcw_idx, preci)
        2378           0 :      call pbuf_get_field(pbuf, prec_str_idx, prec_str)
        2379           0 :      call pbuf_get_field(pbuf, snow_str_idx, snow_str)
        2380           0 :      call pbuf_get_field(pbuf, qcsedten_idx, rc_sed_tend)
        2381           0 :      call pbuf_get_field(pbuf, qrsedten_idx, rr_sed_tend)
        2382           0 :      call pbuf_get_field(pbuf, qisedten_idx, ri_sed_tend)
        2383           0 :      call pbuf_get_field(pbuf, qssedten_idx, rs_sed_tend)
        2384           0 :      call pbuf_get_field(pbuf, vtrmc_idx, vtrmc)
        2385           0 :      call pbuf_get_field(pbuf, umr_idx, umr)
        2386           0 :      call pbuf_get_field(pbuf, vtrmi_idx, vtrmi)
        2387           0 :      call pbuf_get_field(pbuf, ums_idx, ums)
        2388           0 :      call pbuf_get_field(pbuf, qcsevap_idx, rc_sed_evap)
        2389           0 :      call pbuf_get_field(pbuf, qisevap_idx, ri_sed_subl)
        2390             :
        2391             :      ! Calculate liquid precipitation rate (precl) from the total precipitation
        2392             :      ! rate (prect) and the frozen preciptation rate (preci).  This should never
        2393             :      ! be negative, but just to be safe, threshold at 0.
        2394           0 :      precl(:ncol) = max( prect(:ncol) - preci(:ncol), 0.0_r8 )
        2395             :
        2396             :      ! Perform total water and total energy conservation checks.
        2397             :      if ( l_check_conservation ) then
        2398             :
        2399             :         ! Calculate total water in each column.
        2400             :         ! This calculation is the vertically-integrated grand total water (where
        2401             :         ! grand total water is the sum of water vapor, cloud water, rain water,
        2402             :         ! cloud ice, and snow, as well as the amount of water that precipitated
        2403             :         ! to the surface) in each grid column after microphysics, but at the
        2404             :         ! start of hole filling.
        2405           0 :         do icol = 1, ncol
        2406           0 :            grand_total_water_column_start(icol) = 0.0_r8
        2407           0 :            do k = top_lev, pver
        2408             :               grand_total_water_column_start(icol) &
        2409             :               = grand_total_water_column_start(icol) &
        2410           0 :                 + ( state%q(icol,k,1) + ptend%q(icol,k,1) * dt &
        2411           0 :                     + state%q(icol,k,ixcldliq) &
        2412           0 :                     + ptend%q(icol,k,ixcldliq) * dt &
        2413           0 :                     + state%q(icol,k,ixcldice) &
        2414           0 :                     + ptend%q(icol,k,ixcldice) * dt ) &
        2415           0 :                   * state%pdel(icol,k) * rga
        2416           0 :               if ( ixrain > 0 ) then
        2417             :                  grand_total_water_column_start(icol) &
        2418             :                  = grand_total_water_column_start(icol) &
        2419           0 :                    + ( state%q(icol,k,ixrain) + ptend%q(icol,k,ixrain) * dt ) &
        2420           0 :                      * state%pdel(icol,k) * rga
        2421             :               endif
        2422           0 :               if ( ixsnow > 0 ) then
        2423             :                  grand_total_water_column_start(icol) &
        2424             :                  = grand_total_water_column_start(icol) &
        2425           0 :                    + ( state%q(icol,k,ixsnow) + ptend%q(icol,k,ixsnow) * dt ) &
        2426           0 :                      * state%pdel(icol,k) * rga
        2427             :               endif
        2428             :            end do ! k = top_lev, pver
        2429             :            grand_total_water_column_start(icol) &
        2430             :            = grand_total_water_column_start(icol) &
        2431           0 :              + prect(icol) * dt * 1000.0_r8
        2432             :         end do ! icol = 1, ncol
        2433             :
        2434             :         ! Calculate total energy in each column.
        2435             :         ! This calculation is the vertically-integrated total energy in each
        2436             :         ! grid column after microphysics, but at the start of hole filling.
        2437             :         ! Since the microphysics and hole filling code does not directly change
        2438             :         ! kinetic energy, 0.5 * ( u^2 + v^2 ), it can be skipped as part of the
        2439             :         ! energy conservation check.
        2440           0 :         do icol = 1, ncol
        2441           0 :            total_energy_column_start(icol) = 0.0_r8
        2442           0 :            do k = top_lev, pver
        2443             :               total_energy_column_start(icol) &
        2444             :               = total_energy_column_start(icol) &
        2445           0 :                 + ( state%s(icol,k) + ptend%s(icol,k) * dt &
        2446             :                     + ( latvap + latice ) &
        2447           0 :                       * ( state%q(icol,k,1) + ptend%q(icol,k,1) * dt ) &
        2448           0 :                     + latice * ( state%q(icol,k,ixcldliq) &
        2449           0 :                                  + ptend%q(icol,k,ixcldliq) * dt ) ) &
        2450           0 :                   * state%pdel(icol,k) * rga
        2451           0 :               if ( ixrain > 0 ) then
        2452             :                  total_energy_column_start(icol) &
        2453             :                  = total_energy_column_start(icol) &
        2454           0 :                    + latice * ( state%q(icol,k,ixrain) &
        2455           0 :                                 + ptend%q(icol,k,ixrain) * dt ) &
        2456           0 :                      * state%pdel(icol,k) * rga
        2457             :               endif
        2458             :            end do ! k = top_lev, pver
        2459             :            total_energy_column_start(icol) &
        2460             :            = total_energy_column_start(icol) &
        2461           0 :              + latice * precl(icol) * dt * 1000.0_r8
        2462             :         end do ! icol = 1, ncol
        2463             :
        2464             :      endif ! l_check_conservation
        2465             :
        2466             :      ! The fields within state haven't been updated yet, since this is before
        2467             :      ! the call to physics_update.
        2468           0 :      rv_start = state%q(:,:,1)
        2469           0 :      rc_start = state%q(:,:,ixcldliq)
        2470           0 :      if ( ixrain > 0 ) then
        2471           0 :         rr_start = state%q(:,:,ixrain)
        2472             :      endif
        2473           0 :      ri_start = state%q(:,:,ixcldice)
        2474           0 :      if ( ixsnow > 0 ) then
        2475           0 :         rs_start = state%q(:,:,ixsnow)
        2476             :      endif
        2477             :
        2478             :      ! Unpack the current total tendencies for hydrometeor mixing ratio fields.
        2479           0 :      rv_tend = ptend%q(:,:,1)
        2480           0 :      rc_tend = ptend%q(:,:,ixcldliq)
        2481           0 :      if ( ixrain > 0 ) then
        2482           0 :         rr_tend = ptend%q(:,:,ixrain)
        2483             :      endif
        2484           0 :      ri_tend = ptend%q(:,:,ixcldice)
        2485           0 :      if ( ixsnow > 0 ) then
        2486           0 :         rs_tend = ptend%q(:,:,ixsnow)
        2487             :      endif
        2488             :
        2489             :      ! Unpack the current tendency for dry static energy.
        2490           0 :      stend = ptend%s
        2491             :
        2492             :      ! The total hydrometeor tendencies are the sum of microphysics process
        2493             :      ! rates and sedimentation rates.  Calculate the microphysics process
        2494             :      ! tendencies by subtracting the sedimentation tendencies from the overall
        2495             :      ! tendencies.
        2496             :      ! The sedimentation tendencies for cloud water (rc_sed_tend) and cloud ice
        2497             :      ! (ri_sed_tend) include the evaporation of cloud water during sedimentation
        2498             :      ! and the sublimation of cloud ice during sedimentation, respectively.  The
        2499             :      ! true sedimentation of cloud water is the sum of rc_sed_tend and
        2500             :      ! rc_sed_evap, and the true sedimentation of cloud ice is the sum of
        2501             :      ! ri_sed_tend and ri_sed_subl.  Subtract off only the true sedimentation
        2502             :      ! rates, as evaporation and sublimation need to be included in the
        2503             :      ! microphysics process rates.
        2504           0 :      rv_mc_tend(:ncol,:) = rv_tend(:ncol,:)
        2505           0 :      rc_mc_tend(:ncol,:) = rc_tend(:ncol,:) - ( rc_sed_tend(:ncol,:) + rc_sed_evap(:ncol,:) )
        2506           0 :      if ( ixrain > 0 ) then
        2507           0 :         rr_mc_tend(:ncol,:) = rr_tend(:ncol,:) - rr_sed_tend(:ncol,:)
        2508             :      endif
        2509           0 :      ri_mc_tend(:ncol,:) = ri_tend(:ncol,:) - ( ri_sed_tend(:ncol,:) + ri_sed_subl(:ncol,:) )
        2510           0 :      if ( ixsnow > 0 ) then
        2511           0 :         rs_mc_tend(:ncol,:) = rs_tend(:ncol,:) - rs_sed_tend(:ncol,:)
        2512             :      endif
        2513             :
        2514             :      ! This section adjusts microphysics process rate tendencies so that the
        2515             :      ! resulting values of all hydrometeor mixing ratios are greater than or
        2516             :      ! equal to qmin after this section is complete.  Once sedimentation is
        2517             :      ! added back on after this section, some of the hydrometeor mixing ratios
        2518             :      ! may become less than qmin again.
        2519             :      !
        2520             :      ! This section, which again is concerned only with adjusting microphysics
        2521             :      ! process rates, makes use of the following two principles:
        2522             :      !
        2523             :      ! 1) When adjusting the tendencies from microphysics process rates,
        2524             :      !    conserve:
        2525             :      !
        2526             :      !    rv_mc_tend(k) + rc_mc_tend(k) + rr_mc_tend(k)
        2527             :      !    + ri_mc_tend(k) + rs_mc_tend(k) = 0; for all k from top_lev to pver.
        2528             :      !
        2529             :      ! 2) When adjusting the tendencies from microphysics process rates, adjust
        2530             :      !    dry static energy appropriately.  The change in dry static energy
        2531             :      !    is necessary because of phase changes.  This "puts back" the extra dry
        2532             :      !    static energy that was "taken out" when an excessive phase-changing
        2533             :      !    process rate was produced by microphysics.
        2534             :
        2535             :      ! Loop over all columns, performing any tendency adjustments one column
        2536             :      ! at a time.
        2537           0 :      do icol = 1, ncol
        2538             :
        2539             :         ! Loop over all vertical levels, performing any microphysics process
        2540             :         ! tendency adjustments one level at a time.
        2541           0 :         do k = top_lev, pver
        2542             :
        2543             :            ! Find which hydrometeors have positive microphysics process
        2544             :            ! tendencies at this level.
        2545           0 :            if ( rv_mc_tend(icol,k) >= 0.0_r8 ) then
        2546             :               l_pos_rv_mc_tend = .true.
        2547             :            else
        2548           0 :               l_pos_rv_mc_tend = .false.
        2549             :            endif
        2550           0 :            if ( rc_mc_tend(icol,k) >= 0.0_r8 ) then
        2551             :               l_pos_rc_mc_tend = .true.
        2552             :            else
        2553           0 :               l_pos_rc_mc_tend = .false.
        2554             :            endif
        2555           0 :            if ( ixrain > 0 ) then
        2556           0 :               if ( rr_mc_tend(icol,k) >= 0.0_r8 ) then
        2557             :                  l_pos_rr_mc_tend = .true.
        2558             :               else
        2559           0 :                  l_pos_rr_mc_tend = .false.
        2560             :               endif
        2561             :            endif
        2562           0 :            if ( ri_mc_tend(icol,k) >= 0.0_r8 ) then
        2563             :               l_pos_ri_mc_tend = .true.
        2564             :            else
        2565           0 :               l_pos_ri_mc_tend = .false.
        2566             :            endif
        2567           0 :            if ( ixsnow > 0 ) then
        2568           0 :               if ( rs_mc_tend(icol,k) >= 0.0_r8 ) then
        2569             :                  l_pos_rs_mc_tend = .true.
        2570             :               else
        2571           0 :                  l_pos_rs_mc_tend = .false.
        2572             :               endif
        2573             :            endif
        2574             :
        2575             :            !!! Check for holes in water vapor mixing ratio
        2576           0 :            if ( .not. l_pos_rv_mc_tend ) then
        2577             :
        2578             :               ! Calculate the water vapor mixing ratio as it would be with the
        2579             :               ! current microphysics process tendency.
        2580           0 :               rv_curr = rv_start(icol,k) + rv_mc_tend(icol,k) * dt
        2581             :
        2582           0 :               if ( rv_curr < qmin(1) ) then
        2583             :
        2584             :                  ! Microphysics processes are causing a hole in water vapor
        2585             :                  ! mixing ratio.
        2586             :
        2587             :                  ! Calculate the maximum allowable magnitude of (negative) water
        2588             :                  ! vapor microphysics process tendency.
        2589           0 :                  mc_tend_max_mag = ( qmin(1) - rv_start(icol,k) ) / dt
        2590             :
        2591             :                  ! Calculate the amount of the correction that needs to be made
        2592             :                  ! to the water vapor mixing ratio microphysics process
        2593             :                  ! tendency.  This number is positive.
        2594           0 :                  mc_tend_correction = mc_tend_max_mag - rv_mc_tend(icol,k)
        2595             :
        2596             :                  ! Calculate the total amount of positive microphysics process
        2597             :                  ! tendencies for all hydrometeor mixing ratios.
        2598           0 :                  total_mc_positive = 0.0_r8
        2599           0 :                  if ( l_pos_rc_mc_tend ) then
        2600           0 :                     total_mc_positive = total_mc_positive + rc_mc_tend(icol,k)
        2601             :                  endif
        2602           0 :                  if ( ixrain > 0 .and. l_pos_rr_mc_tend ) then
        2603           0 :                     total_mc_positive = total_mc_positive + rr_mc_tend(icol,k)
        2604             :                  endif
        2605           0 :                  if ( l_pos_ri_mc_tend ) then
        2606           0 :                     total_mc_positive = total_mc_positive + ri_mc_tend(icol,k)
        2607             :                  endif
        2608           0 :                  if ( ixsnow > 0 .and. l_pos_rs_mc_tend ) then
        2609           0 :                     total_mc_positive = total_mc_positive + rs_mc_tend(icol,k)
        2610             :                  endif
        2611             :
        2612             :                  ! Calculate the correction ratio.
        2613             :                  ! In principle, this should never be greater than 1 outside of
        2614             :                  ! numerical round-off errors.  This is limited at 1 to be safe.
        2615             :                  mc_correction_ratio &
        2616             :                  = min( mc_tend_correction &
        2617           0 :                         / max( total_mc_positive, 1.0e-30_r8 ), 1.0_r8 )
        2618             :
        2619             :                  ! Adjust (decrease) the tendencies of all positive hydrometeor
        2620             :                  ! mixing ratio tendencies to balance the adjustment (increase)
        2621             :                  ! to the excessively negative water vapor mixing ratio.
        2622             :                  ! Transfer dry static energy appropriately (in response to the
        2623             :                  ! excessive depletion of water vapor).
        2624           0 :                  if ( l_pos_rc_mc_tend ) then
        2625             :                     ! Changing cloud water to water vapor cools and reduces
        2626             :                     ! dry static energy.
        2627             :                     stend(icol,k) &
        2628             :                     = stend(icol,k) &
        2629           0 :                       - latvap * mc_correction_ratio * rc_mc_tend(icol,k)
        2630             :                     ! Update cloud water mixing ratio microphysics tendency.
        2631             :                     rc_mc_tend(icol,k) &
        2632           0 :                     = rc_mc_tend(icol,k) * ( 1.0_r8 - mc_correction_ratio )
        2633             :                  endif
        2634           0 :                  if ( ixrain > 0 .and. l_pos_rr_mc_tend ) then
        2635             :                     ! Changing rain water to water vapor cools and reduces
        2636             :                     ! dry static energy.
        2637             :                     stend(icol,k) &
        2638             :                     = stend(icol,k) &
        2639           0 :                       - latvap * mc_correction_ratio * rr_mc_tend(icol,k)
        2640             :                     ! Update rain water mixing ratio microphysics tendency.
        2641             :                     rr_mc_tend(icol,k) &
        2642           0 :                     = rr_mc_tend(icol,k) * ( 1.0_r8 - mc_correction_ratio )
        2643             :                  endif
        2644           0 :                  if ( l_pos_ri_mc_tend ) then
        2645             :                     ! Changing cloud ice to water vapor cools and reduces
        2646             :                     ! dry static energy.
        2647             :                     stend(icol,k) &
        2648             :                     = stend(icol,k) &
        2649             :                       - ( latvap + latice ) &
        2650           0 :                         * mc_correction_ratio * ri_mc_tend(icol,k)
        2651             :                     ! Update cloud ice mixing ratio microphysics tendency.
        2652             :                     ri_mc_tend(icol,k) &
        2653           0 :                     = ri_mc_tend(icol,k) * ( 1.0_r8 - mc_correction_ratio )
        2654             :                  endif
        2655           0 :                  if ( ixsnow > 0 .and. l_pos_rs_mc_tend ) then
        2656             :                     ! Changing snow to water vapor cools and reduces dry
        2657             :                     ! static energy.
        2658             :                     stend(icol,k) &
        2659             :                     = stend(icol,k) &
        2660             :                       - ( latvap + latice ) &
        2661           0 :                         * mc_correction_ratio * rs_mc_tend(icol,k)
        2662             :                     ! Update snow mixing ratio microphysics tendency.
        2663             :                     rs_mc_tend(icol,k) &
        2664           0 :                     = rs_mc_tend(icol,k) * ( 1.0_r8 - mc_correction_ratio )
        2665             :                  endif
        2666             :
        2667             :                  ! Calculate the new water vapor mixing ratio microphysics
        2668             :                  ! process tendency.  This should be equal to the maximum
        2669             :                  ! magnitude (negative) amount allowed, mc_tend_max_mag.
        2670             :                  rv_mc_tend(icol,k) &
        2671             :                  = rv_mc_tend(icol,k) &
        2672           0 :                    + mc_correction_ratio * total_mc_positive
        2673             :
        2674             :               endif ! rv_curr < qmin(1)
        2675             :
        2676             :            endif ! .not. l_pos_rv_mc_tend
        2677             :
        2678             :            !!! Check for holes in cloud water mixing ratio
        2679           0 :            if ( .not. l_pos_rc_mc_tend ) then
        2680             :
        2681             :               ! Calculate the cloud water mixing ratio as it would be with the
        2682             :               ! current microphysics process tendency.
        2683           0 :               rc_curr = rc_start(icol,k) + rc_mc_tend(icol,k) * dt
        2684             :
        2685           0 :               if ( rc_curr < qmin(ixcldliq) ) then
        2686             :
        2687             :                  ! Microphysics processes are causing a hole in cloud water
        2688             :                  ! mixing ratio.
        2689             :
        2690             :                  ! Calculate the maximum allowable magnitude of (negative) cloud
        2691             :                  ! water microphysics process tendency.
        2692           0 :                  mc_tend_max_mag = ( qmin(ixcldliq) - rc_start(icol,k) ) / dt
        2693             :
        2694             :                  ! Calculate the amount of the correction that needs to be made
        2695             :                  ! to the cloud water mixing ratio microphysics process
        2696             :                  ! tendency.  This number is positive.
        2697           0 :                  mc_tend_correction = mc_tend_max_mag - rc_mc_tend(icol,k)
        2698             :
        2699             :                  ! Calculate the total amount of positive microphysics process
        2700             :                  ! tendencies for all hydrometeor mixing ratios.
        2701           0 :                  total_mc_positive = 0.0_r8
        2702           0 :                  if ( l_pos_rv_mc_tend ) then
        2703           0 :                     total_mc_positive = total_mc_positive + rv_mc_tend(icol,k)
        2704             :                  endif
        2705           0 :                  if ( ixrain > 0 .and. l_pos_rr_mc_tend ) then
        2706           0 :                     total_mc_positive = total_mc_positive + rr_mc_tend(icol,k)
        2707             :                  endif
        2708           0 :                  if ( l_pos_ri_mc_tend ) then
        2709           0 :                     total_mc_positive = total_mc_positive + ri_mc_tend(icol,k)
        2710             :                  endif
        2711           0 :                  if ( ixsnow > 0 .and. l_pos_rs_mc_tend ) then
        2712           0 :                     total_mc_positive = total_mc_positive + rs_mc_tend(icol,k)
        2713             :                  endif
        2714             :
        2715             :                  ! Calculate the correction ratio.
        2716             :                  ! In principle, this should never be greater than 1 outside of
        2717             :                  ! numerical round-off errors.  This is limited at 1 to be safe.
        2718             :                  mc_correction_ratio &
        2719             :                  = min( mc_tend_correction &
        2720           0 :                         / max( total_mc_positive, 1.0e-30_r8 ), 1.0_r8 )
        2721             :
        2722             :                  ! Adjust (decrease) the tendencies of all positive hydrometeor
        2723             :                  ! mixing ratio tendencies to balance the adjustment (increase)
        2724             :                  ! to the excessively negative cloud water mixing ratio.
        2725             :                  ! Transfer dry static energy appropriately (in response to the
        2726             :                  ! excessive depletion of cloud water).
        2727           0 :                  if ( l_pos_rv_mc_tend ) then
        2728             :                     ! Changing water vapor to cloud water heats and increases
        2729             :                     ! dry static energy.
        2730             :                     stend(icol,k) &
        2731             :                     = stend(icol,k) &
        2732           0 :                       + latvap * mc_correction_ratio * rv_mc_tend(icol,k)
        2733             :                     ! Update water vapor mixing ratio microphysics tendency.
        2734             :                     rv_mc_tend(icol,k) &
        2735           0 :                     = rv_mc_tend(icol,k) * ( 1.0_r8 - mc_correction_ratio )
        2736             :                  endif
        2737           0 :                  if ( ixrain > 0 .and. l_pos_rr_mc_tend ) then
        2738             :                     ! Changing rain water to cloud water does not change
        2739             :                     ! dry static energy.
        2740             :                     ! Update rain water mixing ratio microphysics tendency.
        2741             :                     rr_mc_tend(icol,k) &
        2742           0 :                     = rr_mc_tend(icol,k) * ( 1.0_r8 - mc_correction_ratio )
        2743             :                  endif
        2744           0 :                  if ( l_pos_ri_mc_tend ) then
        2745             :                     ! Changing cloud ice to cloud water cools and reduces
        2746             :                     ! dry static energy.
        2747             :                     stend(icol,k) &
        2748             :                     = stend(icol,k) &
        2749           0 :                       - latice * mc_correction_ratio * ri_mc_tend(icol,k)
        2750             :                     ! Update cloud ice mixing ratio microphysics tendency.
        2751             :                     ri_mc_tend(icol,k) &
        2752           0 :                     = ri_mc_tend(icol,k) * ( 1.0_r8 - mc_correction_ratio )
        2753             :                  endif
        2754           0 :                  if ( ixsnow > 0 .and. l_pos_rs_mc_tend ) then
        2755             :                     ! Changing snow to cloud water cools and reduces dry
        2756             :                     ! static energy.
        2757             :                     stend(icol,k) &
        2758             :                     = stend(icol,k) &
        2759           0 :                       - latice * mc_correction_ratio * rs_mc_tend(icol,k)
        2760             :                     ! Update snow mixing ratio microphysics tendency.
        2761             :                     rs_mc_tend(icol,k) &
        2762           0 :                     = rs_mc_tend(icol,k) * ( 1.0_r8 - mc_correction_ratio )
        2763             :                  endif
        2764             :
        2765             :                  ! Calculate the new cloud water mixing ratio microphysics
        2766             :                  ! process tendency.  This should be equal to the maximum
        2767             :                  ! magnitude (negative) amount allowed, mc_tend_max_mag.
        2768             :                  rc_mc_tend(icol,k) &
        2769             :                  = rc_mc_tend(icol,k) &
        2770           0 :                    + mc_correction_ratio * total_mc_positive
        2771             :
        2772             :               endif ! rc_curr < qmin(ixcldliq)
        2773             :
        2774             :            endif ! .not. l_pos_rc_mc_tend
        2775             :
        2776             :            !!! Check for holes in rain water mixing ratio
        2777           0 :            if ( ixrain > 0 .and. ( .not. l_pos_rr_mc_tend ) ) then
        2778             :
        2779             :               ! Calculate the rain water mixing ratio as it would be with the
        2780             :               ! current microphysics process tendency.
        2781           0 :               rr_curr = rr_start(icol,k) + rr_mc_tend(icol,k) * dt
        2782             :
        2783           0 :               if ( rr_curr < qmin(ixrain) ) then
        2784             :
        2785             :                  ! Microphysics processes are causing a hole in rain water
        2786             :                  ! mixing ratio.
        2787             :
        2788             :                  ! Calculate the maximum allowable magnitude of (negative) rain
        2789             :                  ! water microphysics process tendency.
        2790           0 :                  mc_tend_max_mag = ( qmin(ixrain) - rr_start(icol,k) ) / dt
        2791             :
        2792             :                  ! Calculate the amount of the correction that needs to be made
        2793             :                  ! to the rain water mixing ratio microphysics process
        2794             :                  ! tendency.  This number is positive.
        2795           0 :                  mc_tend_correction = mc_tend_max_mag - rr_mc_tend(icol,k)
        2796             :
        2797             :                  ! Calculate the total amount of positive microphysics process
        2798             :                  ! tendencies for all hydrometeor mixing ratios.
        2799           0 :                  total_mc_positive = 0.0_r8
        2800           0 :                  if ( l_pos_rv_mc_tend ) then
        2801           0 :                     total_mc_positive = total_mc_positive + rv_mc_tend(icol,k)
        2802             :                  endif
        2803           0 :                  if ( l_pos_rc_mc_tend ) then
        2804           0 :                     total_mc_positive = total_mc_positive + rc_mc_tend(icol,k)
        2805             :                  endif
        2806           0 :                  if ( l_pos_ri_mc_tend ) then
        2807           0 :                     total_mc_positive = total_mc_positive + ri_mc_tend(icol,k)
        2808             :                  endif
        2809           0 :                  if ( ixsnow > 0 .and. l_pos_rs_mc_tend ) then
        2810           0 :                     total_mc_positive = total_mc_positive + rs_mc_tend(icol,k)
        2811             :                  endif
        2812             :
        2813             :                  ! Calculate the correction ratio.
        2814             :                  ! In principle, this should never be greater than 1 outside of
        2815             :                  ! numerical round-off errors.  This is limited at 1 to be safe.
        2816             :                  mc_correction_ratio &
        2817             :                  = min( mc_tend_correction &
        2818           0 :                         / max( total_mc_positive, 1.0e-30_r8 ), 1.0_r8 )
        2819             :
        2820             :                  ! Adjust (decrease) the tendencies of all positive hydrometeor
        2821             :                  ! mixing ratio tendencies to balance the adjustment (increase)
        2822             :                  ! to the excessively negative rain water mixing ratio.
        2823             :                  ! Transfer dry static energy appropriately (in response to the
        2824             :                  ! excessive depletion of rain water).
        2825           0 :                  if ( l_pos_rv_mc_tend ) then
        2826             :                     ! Changing water vapor to rain water heats and increases
        2827             :                     ! dry static energy.
        2828             :                     stend(icol,k) &
        2829             :                     = stend(icol,k) &
        2830           0 :                       + latvap * mc_correction_ratio * rv_mc_tend(icol,k)
        2831             :                     ! Update water vapor mixing ratio microphysics tendency.
        2832             :                     rv_mc_tend(icol,k) &
        2833           0 :                     = rv_mc_tend(icol,k) * ( 1.0_r8 - mc_correction_ratio )
        2834             :                  endif
        2835           0 :                  if ( l_pos_rc_mc_tend ) then
        2836             :                     ! Changing cloud water to rain water does not change
        2837             :                     ! dry static energy.
        2838             :                     ! Update cloud water mixing ratio microphysics tendency.
        2839             :                     rc_mc_tend(icol,k) &
        2840           0 :                     = rc_mc_tend(icol,k) * ( 1.0_r8 - mc_correction_ratio )
        2841             :                  endif
        2842           0 :                  if ( l_pos_ri_mc_tend ) then
        2843             :                     ! Changing cloud ice to rain water cools and reduces
        2844             :                     ! dry static energy.
        2845             :                     stend(icol,k) &
        2846             :                     = stend(icol,k) &
        2847           0 :                       - latice * mc_correction_ratio * ri_mc_tend(icol,k)
        2848             :                     ! Update cloud ice mixing ratio microphysics tendency.
        2849             :                     ri_mc_tend(icol,k) &
        2850           0 :                     = ri_mc_tend(icol,k) * ( 1.0_r8 - mc_correction_ratio )
        2851             :                  endif
        2852           0 :                  if ( ixsnow > 0 .and. l_pos_rs_mc_tend ) then
        2853             :                     ! Changing snow to rain water cools and reduces dry
        2854             :                     ! static energy.
        2855             :                     stend(icol,k) &
        2856             :                     = stend(icol,k) &
        2857           0 :                       - latice * mc_correction_ratio * rs_mc_tend(icol,k)
        2858             :                     ! Update snow mixing ratio microphysics tendency.
        2859             :                     rs_mc_tend(icol,k) &
        2860           0 :                     = rs_mc_tend(icol,k) * ( 1.0_r8 - mc_correction_ratio )
        2861             :                  endif
        2862             :
        2863             :                  ! Calculate the new rain water mixing ratio microphysics
        2864             :                  ! process tendency.  This should be equal to the maximum
        2865             :                  ! magnitude (negative) amount allowed, mc_tend_max_mag.
        2866             :                  rr_mc_tend(icol,k) &
        2867             :                  = rr_mc_tend(icol,k) &
        2868           0 :                    + mc_correction_ratio * total_mc_positive
        2869             :
        2870             :               endif ! rr_curr < qmin(ixrain)
        2871             :
        2872             :            endif ! ixrain > 0 .and. ( .not. l_pos_rr_mc_tend )
        2873             :
        2874             :            !!! Check for holes in cloud ice mixing ratio
        2875           0 :            if ( .not. l_pos_ri_mc_tend ) then
        2876             :
        2877             :               ! Calculate the cloud ice mixing ratio as it would be with the
        2878             :               ! current microphysics process tendency.
        2879           0 :               ri_curr = ri_start(icol,k) + ri_mc_tend(icol,k) * dt
        2880             :
        2881           0 :               if ( ri_curr < qmin(ixcldice) ) then
        2882             :
        2883             :                  ! Microphysics processes are causing a hole in cloud ice
        2884             :                  ! mixing ratio.
        2885             :
        2886             :                  ! Calculate the maximum allowable magnitude of (negative) cloud
        2887             :                  ! ice microphysics process tendency.
        2888           0 :                  mc_tend_max_mag = ( qmin(ixcldice) - ri_start(icol,k) ) / dt
        2889             :
        2890             :                  ! Calculate the amount of the correction that needs to be made
        2891             :                  ! to the cloud ice mixing ratio microphysics process
        2892             :                  ! tendency.  This number is positive.
        2893           0 :                  mc_tend_correction = mc_tend_max_mag - ri_mc_tend(icol,k)
        2894             :
        2895             :                  ! Calculate the total amount of positive microphysics process
        2896             :                  ! tendencies for all hydrometeor mixing ratios.
        2897           0 :                  total_mc_positive = 0.0_r8
        2898           0 :                  if ( l_pos_rv_mc_tend ) then
        2899           0 :                     total_mc_positive = total_mc_positive + rv_mc_tend(icol,k)
        2900             :                  endif
        2901           0 :                  if ( l_pos_rc_mc_tend ) then
        2902           0 :                     total_mc_positive = total_mc_positive + rc_mc_tend(icol,k)
        2903             :                  endif
        2904           0 :                  if ( ixrain > 0 .and. l_pos_rr_mc_tend ) then
        2905           0 :                     total_mc_positive = total_mc_positive + rr_mc_tend(icol,k)
        2906             :                  endif
        2907           0 :                  if ( ixsnow > 0 .and. l_pos_rs_mc_tend ) then
        2908           0 :                     total_mc_positive = total_mc_positive + rs_mc_tend(icol,k)
        2909             :                  endif
        2910             :
        2911             :                  ! Calculate the correction ratio.
        2912             :                  ! In principle, this should never be greater than 1 outside of
        2913             :                  ! numerical round-off errors.  This is limited at 1 to be safe.
        2914             :                  mc_correction_ratio &
        2915             :                  = min( mc_tend_correction &
        2916           0 :                         / max( total_mc_positive, 1.0e-30_r8 ), 1.0_r8 )
        2917             :
        2918             :                  ! Adjust (decrease) the tendencies of all positive hydrometeor
        2919             :                  ! mixing ratio tendencies to balance the adjustment (increase)
        2920             :                  ! to the excessively negative cloud ice mixing ratio.
        2921             :                  ! Transfer dry static energy appropriately (in response to the
        2922             :                  ! excessive depletion of cloud ice).
        2923           0 :                  if ( l_pos_rv_mc_tend ) then
        2924             :                     ! Changing water vapor to cloud ice heats and increases
        2925             :                     ! dry static energy.
        2926             :                     stend(icol,k) &
        2927             :                     = stend(icol,k) &
        2928             :                       + ( latvap + latice ) &
        2929           0 :                         * mc_correction_ratio * rv_mc_tend(icol,k)
        2930             :                     ! Update water vapor mixing ratio microphysics tendency.
        2931             :                     rv_mc_tend(icol,k) &
        2932           0 :                     = rv_mc_tend(icol,k) * ( 1.0_r8 - mc_correction_ratio )
        2933             :                  endif
        2934           0 :                  if ( l_pos_rc_mc_tend ) then
        2935             :                     ! Changing cloud water to cloud ice heats and increases
        2936             :                     ! dry static energy.
        2937             :                     stend(icol,k) &
        2938             :                     = stend(icol,k) &
        2939           0 :                       + latice * mc_correction_ratio * rc_mc_tend(icol,k)
        2940             :                     ! Update cloud water mixing ratio microphysics tendency.
        2941             :                     rc_mc_tend(icol,k) &
        2942           0 :                     = rc_mc_tend(icol,k) * ( 1.0_r8 - mc_correction_ratio )
        2943             :                  endif
        2944           0 :                  if ( ixrain > 0 .and. l_pos_rr_mc_tend ) then
        2945             :                     ! Changing rain water to cloud ice heats and increases
        2946             :                     ! dry static energy.
        2947             :                     stend(icol,k) &
        2948             :                     = stend(icol,k) &
        2949           0 :                       + latice * mc_correction_ratio * rr_mc_tend(icol,k)
        2950             :                     ! Update rain water mixing ratio microphysics tendency.
        2951             :                     rr_mc_tend(icol,k) &
        2952           0 :                     = rr_mc_tend(icol,k) * ( 1.0_r8 - mc_correction_ratio )
        2953             :                  endif
        2954           0 :                  if ( ixsnow > 0 .and. l_pos_rs_mc_tend ) then
        2955             :                     ! Changing snow to cloud ice does not change dry static
        2956             :                     ! energy.
        2957             :                     ! Update snow mixing ratio microphysics tendency.
        2958             :                     rs_mc_tend(icol,k) &
        2959           0 :                     = rs_mc_tend(icol,k) * ( 1.0_r8 - mc_correction_ratio )
        2960             :                  endif
        2961             :
        2962             :                  ! Calculate the new cloud ice mixing ratio microphysics
        2963             :                  ! process tendency.  This should be equal to the maximum
        2964             :                  ! magnitude (negative) amount allowed, mc_tend_max_mag.
        2965             :                  ri_mc_tend(icol,k) &
        2966             :                  = ri_mc_tend(icol,k) &
        2967           0 :                    + mc_correction_ratio * total_mc_positive
        2968             :
        2969             :               endif ! ri_curr < qmin(ixcldice)
        2970             :
        2971             :            endif ! .not. l_pos_ri_mc_tend
        2972             :
        2973             :            !!! Check for holes in snow mixing ratio
        2974           0 :            if ( ixsnow > 0 .and. ( .not. l_pos_rs_mc_tend ) ) then
        2975             :
        2976             :               ! Calculate the snow mixing ratio as it would be with the
        2977             :               ! current microphysics process tendency.
        2978           0 :               rs_curr = rs_start(icol,k) + rs_mc_tend(icol,k) * dt
        2979             :
        2980           0 :               if ( rs_curr < qmin(ixsnow) ) then
        2981             :
        2982             :                  ! Microphysics processes are causing a hole in snow mixing
        2983             :                  ! ratio.
        2984             :
        2985             :                  ! Calculate the maximum allowable magnitude of (negative) snow
        2986             :                  ! microphysics process tendency.
        2987           0 :                  mc_tend_max_mag = ( qmin(ixsnow) - rs_start(icol,k) ) / dt
        2988             :
        2989             :                  ! Calculate the amount of the correction that needs to be made
        2990             :                  ! to the snow mixing ratio microphysics process tendency.
        2991             :                  ! This number is positive.
        2992           0 :                  mc_tend_correction = mc_tend_max_mag - rs_mc_tend(icol,k)
        2993             :
        2994             :                  ! Calculate the total amount of positive microphysics process
        2995             :                  ! tendencies for all hydrometeor mixing ratios.
        2996           0 :                  total_mc_positive = 0.0_r8
        2997           0 :                  if ( l_pos_rv_mc_tend ) then
        2998           0 :                     total_mc_positive = total_mc_positive + rv_mc_tend(icol,k)
        2999             :                  endif
        3000           0 :                  if ( l_pos_rc_mc_tend ) then
        3001           0 :                     total_mc_positive = total_mc_positive + rc_mc_tend(icol,k)
        3002             :                  endif
        3003           0 :                  if ( ixrain > 0 .and. l_pos_rr_mc_tend ) then
        3004           0 :                     total_mc_positive = total_mc_positive + rr_mc_tend(icol,k)
        3005             :                  endif
        3006           0 :                  if ( l_pos_ri_mc_tend ) then
        3007           0 :                     total_mc_positive = total_mc_positive + ri_mc_tend(icol,k)
        3008             :                  endif
        3009             :
        3010             :                  ! Calculate the correction ratio.
        3011             :                  ! In principle, this should never be greater than 1 outside of
        3012             :                  ! numerical round-off errors.  This is limited at 1 to be safe.
        3013             :                  mc_correction_ratio &
        3014             :                  = min( mc_tend_correction &
        3015           0 :                         / max( total_mc_positive, 1.0e-30_r8 ), 1.0_r8 )
        3016             :
        3017             :                  ! Adjust (decrease) the tendencies of all positive hydrometeor
        3018             :                  ! mixing ratio tendencies to balance the adjustment (increase)
        3019             :                  ! to the excessively negative snow mixing ratio.
        3020             :                  ! Transfer dry static energy appropriately (in response to the
        3021             :                  ! excessive depletion of snow).
        3022           0 :                  if ( l_pos_rv_mc_tend ) then
        3023             :                     ! Changing water vapor to snow heats and increases dry
        3024             :                     ! static energy.
        3025             :                     stend(icol,k) &
        3026             :                     = stend(icol,k) &
        3027             :                       + ( latvap + latice ) &
        3028           0 :                         * mc_correction_ratio * rv_mc_tend(icol,k)
        3029             :                     ! Update water vapor mixing ratio microphysics tendency.
        3030             :                     rv_mc_tend(icol,k) &
        3031           0 :                     = rv_mc_tend(icol,k) * ( 1.0_r8 - mc_correction_ratio )
        3032             :                  endif
        3033           0 :                  if ( l_pos_rc_mc_tend ) then
        3034             :                     ! Changing cloud water to snow heats and increases dry
        3035             :                     ! static energy.
        3036             :                     stend(icol,k) &
        3037             :                     = stend(icol,k) &
        3038           0 :                       + latice * mc_correction_ratio * rc_mc_tend(icol,k)
        3039             :                     ! Update cloud water mixing ratio microphysics tendency.
        3040             :                     rc_mc_tend(icol,k) &
        3041           0 :                     = rc_mc_tend(icol,k) * ( 1.0_r8 - mc_correction_ratio )
        3042             :                  endif
        3043           0 :                  if ( ixrain > 0 .and. l_pos_rr_mc_tend ) then
        3044             :                     ! Changing rain water to snow heats and increases dry
        3045             :                     ! static energy.
        3046             :                     stend(icol,k) &
        3047             :                     = stend(icol,k) &
        3048           0 :                       + latice * mc_correction_ratio * rr_mc_tend(icol,k)
        3049             :                     ! Update rain water mixing ratio microphysics tendency.
        3050             :                     rr_mc_tend(icol,k) &
        3051           0 :                     = rr_mc_tend(icol,k) * ( 1.0_r8 - mc_correction_ratio )
        3052             :                  endif
        3053           0 :                  if ( l_pos_ri_mc_tend ) then
        3054             :                     ! Changing cloud ice to snow does not change dry static
        3055             :                     ! energy.
        3056             :                     ! Update cloud ice mixing ratio microphysics tendency.
        3057             :                     ri_mc_tend(icol,k) &
        3058           0 :                     = ri_mc_tend(icol,k) * ( 1.0_r8 - mc_correction_ratio )
        3059             :                  endif
        3060             :
        3061             :                  ! Calculate the new snow mixing ratio microphysics process
        3062             :                  ! tendency.  This should be equal to the maximum magnitude
        3063             :                  ! (negative) amount allowed, mc_tend_max_mag.
        3064             :                  rs_mc_tend(icol,k) &
        3065             :                  = rs_mc_tend(icol,k) &
        3066           0 :                    + mc_correction_ratio * total_mc_positive
        3067             :
        3068             :               endif ! rs_curr < qmin(ixsnow)
        3069             :
        3070             :            endif ! ixsnow > 0 .and. ( .not. l_pos_rs_mc_tend )
        3071             :
        3072             :         end do ! k = top_lev, pver
        3073             :
        3074             :      end do ! icol = 1, ncol
        3075             :
        3076             :      ! Calculate the new overall tendencies by adding the sedimentation
        3077             :      ! tendencies back onto the new microphysics process tendencies.
        3078             :      ! For cloud water and cloud ice, the sedimentation tendencies that are
        3079             :      ! added back on are the true sedimentation tendencies.  For cloud water,
        3080             :      ! this is the sum of rc_sed_tend and rc_sed_evap, and for cloud ice, this
        3081             :      ! is the sum of ri_sed_tend and ri_sed_subl.
        3082           0 :      rv_tend(:ncol,:) = rv_mc_tend(:ncol,:)
        3083           0 :      rc_tend(:ncol,:) = rc_mc_tend(:ncol,:) + ( rc_sed_tend(:ncol,:) + rc_sed_evap(:ncol,:) )
        3084           0 :      if ( ixrain > 0 ) then
        3085           0 :         rr_tend(:ncol,:) = rr_mc_tend(:ncol,:) + rr_sed_tend(:ncol,:)
        3086             :      endif
        3087           0 :      ri_tend(:ncol,:) = ri_mc_tend(:ncol,:) + ( ri_sed_tend(:ncol,:) + ri_sed_subl(:ncol,:) )
        3088           0 :      if ( ixsnow > 0 ) then
        3089           0 :         rs_tend(:ncol,:) = rs_mc_tend(:ncol,:) + rs_sed_tend(:ncol,:)
        3090             :      endif
        3091             :
        3092             :      ! Now that the original sedimentation tendency has been added to the
        3093             :      ! new microphysics process tendency, the new total microphysics tendency
        3094             :      ! can still cause holes to form.  After the microphysics process rates were
        3095             :      ! adjusted, the values of the hydrometeor fields were greater than or equal
        3096             :      ! to 0 at all grid levels, which means their vertical integrals were also
        3097             :      ! greater than or equal to 0.  Sedimentation by itself has a vertical
        3098             :      ! integral of 0 (including the amount that sedimented to the surface).
        3099             :      ! This means that after the microphysics process rates have been adjusted
        3100             :      ! and sedimentation has been added back on, the resulting hydrometeor
        3101             :      ! fields all still have vertical integrals that are greater than or equal
        3102             :      ! to 0.  Holes that develop at any particular grid level can be filled.
        3103             :      ! These holes can be filled conservatively using the sedimentation hole
        3104             :      ! filler.
        3105             :      if ( l_sed_hole_fill ) then
        3106             :
        3107             :         ! This section makes use of the following principle:
        3108             :         !
        3109             :         ! 3) When adjusting the hydrometeor tendency from sedimentation of a
        3110             :         !    liquid hydrometeor (cloud water or rain water), conserve:
        3111             :         !
        3112             :         !    SUM(k=top_lev:pver) ( rc_sed_tend(k) + rr_sed_tend(k) )
        3113             :         !                        * dt * pdel(k) / g
        3114             :         !    + precl * dt * 1000 = 0.
        3115             :
        3116             :         ! Call the sedimentation hole filler for rain water mixing ratio.
        3117             :         ! This can update rr_tend and precl.
        3118           0 :         if ( ixrain > 0 ) then
        3119             :            call fill_holes_sedimentation( dt, ncol, rr_start, state%pdel, &
        3120           0 :                                           umr, state%zi, qmin(ixrain), &
        3121           0 :                                           rr_tend, precl )
        3122             :         endif ! ixrain > 0
        3123             :
        3124             :         ! Call the sedimentation hole filler for cloud water mixing ratio.
        3125             :         ! This can update rc_tend and precl.
        3126             :         call fill_holes_sedimentation( dt, ncol, rc_start, state%pdel, &
        3127           0 :                                        vtrmc, state%zi, qmin(ixcldliq), &
        3128           0 :                                        rc_tend, precl )
        3129             :
        3130             :         ! Occasionally, a situation can occur where filling a hole in rain can
        3131             :         ! deplete all the surface liquid-phase precipitation (precl), resulting
        3132             :         ! in not enough water mass in the vertical profile of cloud water to
        3133             :         ! fill a hole in cloud water.  When this happens, there must be liquid
        3134             :         ! water found in the vertical profile of rain, so pull the water from
        3135             :         ! rain to fill any remaining holes in cloud water.
        3136           0 :         if ( ixrain > 0 ) then
        3137             :            call fill_holes_same_phase_vert( dt, ncol, rc_start, rr_start, &
        3138           0 :                                             state%pdel, qmin(ixcldliq), &
        3139           0 :                                             qmin(ixrain), &
        3140           0 :                                             rc_tend, rr_tend )
        3141             :         endif ! ixrain > 0
        3142             :
        3143             :         ! This section makes use of the following principle:
        3144             :         !
        3145             :         ! 4) When adjusting the hydrometeor tendency from sedimentation of a
        3146             :         !    frozen hydrometeor (cloud ice or snow), conserve:
        3147             :         !
        3148             :         !    SUM(k=top_lev:pver) ( ri_sed_tend(k) + rs_sed_tend(k) )
        3149             :         !                        * dt * pdel(k) / g
        3150             :         !    + preci * dt * 1000 = 0.
        3151             :
        3152             :         ! Call the sedimentation hole filler for snow mixing ratio.
        3153             :         ! This can update rs_tend and preci.
        3154           0 :         if ( ixsnow > 0 ) then
        3155             :            call fill_holes_sedimentation( dt, ncol, rs_start, state%pdel, &
        3156           0 :                                           ums, state%zi, qmin(ixsnow), &
        3157           0 :                                           rs_tend, preci )
        3158             :         endif ! ixsnow > 0
        3159             :
        3160             :         ! Call the sedimentation hole filler for cloud ice mixing ratio.
        3161             :         ! This can update ri_tend and preci.
        3162             :         call fill_holes_sedimentation( dt, ncol, ri_start, state%pdel, &
        3163           0 :                                        vtrmi, state%zi, qmin(ixcldice), &
        3164           0 :                                        ri_tend, preci )
        3165             :
        3166             :         ! Occasionally, a situation can occur where filling a hole in snow can
        3167             :         ! deplete all the surface ice-phase precipitation (preci), resulting
        3168             :         ! in not enough water mass in the vertical profile of cloud ice to
        3169             :         ! fill a hole in cloud ice.  When this happens, there must be ice-phase
        3170             :         ! water found in the vertical profile of snow, so pull the water from
        3171             :         ! snow to fill any remaining holes in cloud ice.
        3172           0 :         if ( ixsnow > 0 ) then
        3173             :            call fill_holes_same_phase_vert( dt, ncol, ri_start, rs_start, &
        3174           0 :                                             state%pdel, qmin(ixcldice), &
        3175           0 :                                             qmin(ixsnow), &
        3176           0 :                                             ri_tend, rs_tend )
        3177             :         endif  ! ixsnow > 0
        3178             :
        3179             :         ! Update the total precipitation rate (prect) from the updated liquid
        3180             :         ! precipitation rate (precl) and the updated frozen preciptation rate
        3181             :         ! (preci).
        3182           0 :         prect(:ncol) = precl(:ncol) + preci(:ncol)
        3183             :
        3184             :         ! The MG code sets prec_str equal to prect (prec_pcw) and snow_str equal
        3185             :         ! to preci (snow_pcw).  The prec_str and snow_str variables are used
        3186             :         ! in the calculations for energy and water conservation.  Since prect
        3187             :         ! and preci are adjusted here, when necessary, prec_str and snow_str
        3188             :         ! also need to be adjusted.
        3189           0 :         prec_str(:ncol) = prect(:ncol)
        3190           0 :         snow_str(:ncol) = preci(:ncol)
        3191             :
        3192             :      endif ! l_sed_hole_fill
        3193             :
        3194             :      ! The updated total microphysics tendencies after hole filling have not
        3195             :      ! been used to update ptend yet, so record the budget terms for hole
        3196             :      ! filling first.
        3197           0 :      rv_hf_tend = rv_tend - ptend%q(:,:,1)
        3198           0 :      rc_hf_tend = rc_tend - ptend%q(:,:,ixcldliq)
        3199           0 :      if ( ixrain > 0 ) then
        3200           0 :         rr_hf_tend = rr_tend - ptend%q(:,:,ixrain)
        3201             :      endif ! ixrain > 0
        3202           0 :      ri_hf_tend = ri_tend - ptend%q(:,:,ixcldice)
        3203           0 :      if ( ixsnow > 0 ) then
        3204           0 :         rs_hf_tend = rs_tend - ptend%q(:,:,ixsnow)
        3205             :      endif ! ixsnow > 0
        3206             :
        3207             :      ! The updated dry static energy tendency after hole filling has not been
        3208             :      ! used to update ptend yet, so record the budget term for hole filling
        3209             :      ! first.
        3210           0 :      s_hf_tend = stend - ptend%s
        3211             :
        3212             :      ! Pack the current total tendencies for hydrometeor mixing ratio fields.
        3213           0 :      ptend%q(:,:,1) = rv_tend
        3214           0 :      ptend%q(:,:,ixcldliq) = rc_tend
        3215           0 :      if ( ixrain > 0 ) then
        3216           0 :         ptend%q(:,:,ixrain) = rr_tend
        3217             :      endif
        3218           0 :      ptend%q(:,:,ixcldice) = ri_tend
        3219           0 :      if ( ixsnow > 0 ) then
        3220           0 :         ptend%q(:,:,ixsnow) = rs_tend
        3221             :      endif
        3222             :
        3223             :      ! Pack the current tendency for dry static energy.
        3224           0 :      ptend%s = stend
        3225             :
        3226             :      ! Output stats for hole filling tendencies.
        3227           0 :      call outfld( 'QVHFTEN', rv_hf_tend, pcols, state%lchnk )
        3228           0 :      call outfld( 'QCHFTEN', rc_hf_tend, pcols, state%lchnk )
        3229           0 :      call outfld( 'QRHFTEN', rr_hf_tend, pcols, state%lchnk )
        3230           0 :      call outfld( 'QIHFTEN', ri_hf_tend, pcols, state%lchnk )
        3231           0 :      call outfld( 'QSHFTEN', rs_hf_tend, pcols, state%lchnk )
        3232           0 :      call outfld( 'THFTEN', s_hf_tend / cpair, pcols, state%lchnk )
        3233             :
        3234             :      ! Perform total water and total energy conservation checks.
        3235             :      if ( l_check_conservation ) then
        3236             :
        3237             :         ! Calculate total water in each grid column.
        3238             :         ! This calculation is the vertically-integrated grand total water
        3239             :         ! in each grid column updated for all microphysics and hole filling.
        3240             :         ! This includes the amount that precipitated to the surface.
        3241           0 :         do icol = 1, ncol
        3242           0 :            grand_total_water_column_finish(icol) = 0.0_r8
        3243           0 :            do k = top_lev, pver
        3244             :               grand_total_water_column_finish(icol) &
        3245             :               = grand_total_water_column_finish(icol) &
        3246           0 :                 + ( state%q(icol,k,1) &
        3247           0 :                     + ptend%q(icol,k,1) * dt &
        3248           0 :                     + state%q(icol,k,ixcldliq) &
        3249           0 :                     + ptend%q(icol,k,ixcldliq) * dt &
        3250           0 :                     + state%q(icol,k,ixcldice) &
        3251           0 :                     + ptend%q(icol,k,ixcldice) * dt ) &
        3252           0 :                   * state%pdel(icol,k) * rga
        3253           0 :               if ( ixrain > 0 ) then
        3254             :                  grand_total_water_column_finish(icol) &
        3255             :                  = grand_total_water_column_finish(icol) &
        3256           0 :                    + ( state%q(icol,k,ixrain) + ptend%q(icol,k,ixrain) * dt ) &
        3257           0 :                      * state%pdel(icol,k) * rga
        3258             :               endif
        3259           0 :               if ( ixsnow > 0 ) then
        3260             :                  grand_total_water_column_finish(icol) &
        3261             :                  = grand_total_water_column_finish(icol) &
        3262           0 :                    + ( state%q(icol,k,ixsnow) + ptend%q(icol,k,ixsnow) * dt ) &
        3263           0 :                      * state%pdel(icol,k) * rga
        3264             :               endif
        3265             :            end do ! k = top_lev, pver
        3266             :            grand_total_water_column_finish(icol) &
        3267             :            = grand_total_water_column_finish(icol) &
        3268           0 :              + prect(icol) * dt * 1000.0_r8
        3269             :         end do ! icol = 1, ncol
        3270             :
        3271             :         ! Calculate total energy in each column.
        3272             :         ! This calculation is the vertically-integrated total energy in each
        3273             :         ! grid column updated for all microphysics and hole filling.  This
        3274             :         ! includes the amount that precipitated to the surface.  Since, the
        3275             :         ! microphysics code does not directly change kinetic energy,
        3276             :         ! 0.5 * ( u^2 + v^2 ), it can be skipped as part of the energy
        3277             :         ! conservation check.
        3278           0 :         do icol = 1, ncol
        3279           0 :            total_energy_column_finish(icol) = 0.0_r8
        3280           0 :            do k = top_lev, pver
        3281             :               total_energy_column_finish(icol) &
        3282             :               = total_energy_column_finish(icol) &
        3283           0 :                 + ( state%s(icol,k) + ptend%s(icol,k) * dt &
        3284             :                     + ( latvap + latice ) &
        3285           0 :                       * ( state%q(icol,k,1) + ptend%q(icol,k,1) * dt ) &
        3286           0 :                     + latice * ( state%q(icol,k,ixcldliq) &
        3287           0 :                                  + ptend%q(icol,k,ixcldliq) * dt ) ) &
        3288           0 :                   * state%pdel(icol,k) * rga
        3289           0 :               if ( ixrain > 0 ) then
        3290             :                  total_energy_column_finish(icol) &
        3291             :                  = total_energy_column_finish(icol) &
        3292           0 :                    + latice * ( state%q(icol,k,ixrain) &
        3293           0 :                                 + ptend%q(icol,k,ixrain) * dt ) &
        3294           0 :                      * state%pdel(icol,k) * rga
        3295             :               endif
        3296             :            end do ! k = top_lev, pver
        3297             :            total_energy_column_finish(icol) &
        3298             :            = total_energy_column_finish(icol) &
        3299           0 :              + latice * precl(icol) * dt * 1000.0_r8
        3300             :         end do ! icol = 1, ncol
        3301             :
        3302             :         ! Calculate the total relative error in each grid column.
        3303           0 :         do icol = 1, ncol
        3304             :
        3305           0 :            tot_water_rel_err(icol) &
        3306             :            = abs( ( grand_total_water_column_finish(icol) &
        3307             :                     - grand_total_water_column_start(icol) ) ) &
        3308             :              / min( grand_total_water_column_finish(icol), &
        3309           0 :                     grand_total_water_column_start(icol) )
        3310             :
        3311             :            tot_energy_rel_err(icol) &
        3312             :            = abs( ( total_energy_column_finish(icol) &
        3313             :                     - total_energy_column_start(icol) ) ) &
        3314             :              / min( total_energy_column_finish(icol), &
        3315           0 :                     total_energy_column_start(icol) )
        3316             :
        3317             :         end do ! icol = 1, ncol
        3318             :
        3319             :         ! Print an error message if any total water relative error is found to
        3320             :         ! be greater than the threshold.
        3321           0 :         if ( any( tot_water_rel_err(:ncol) >= err_thresh ) ) then
        3322           0 :            write(iulog,*) "Water conservation error reported in hole filling"
        3323           0 :            do icol = 1, ncol
        3324           0 :               if ( tot_water_rel_err(icol) >= err_thresh ) then
        3325           0 :                  write(iulog,*) "Column = ", icol, &
        3326           0 :                           "Relative error = ", tot_water_rel_err(icol), &
        3327           0 :                           "Column-integrated grand total water at start = ", &
        3328           0 :                           grand_total_water_column_start(icol), &
        3329           0 :                           "Column-integrated grand total water at finish = ", &
        3330           0 :                           grand_total_water_column_finish(icol)
        3331             :               endif ! tot_water_rel_err(icol) >= err_thresh
        3332             :            end do ! icol = 1, ncol
        3333             :         endif ! any( tot_water_rel_err >= err_thresh )
        3334             :
        3335             :         ! Print an error message if any total energy relative error is found to
        3336             :         ! be greater than the threshold.
        3337           0 :         if ( any( tot_energy_rel_err(:ncol) >= err_thresh ) ) then
        3338           0 :            write(iulog,*) "Energy conservation error reported in hole filling"
        3339           0 :            do icol = 1, ncol
        3340           0 :               if ( tot_energy_rel_err(icol) >= err_thresh ) then
        3341           0 :                  write(iulog,*) "Column = ", icol, &
        3342           0 :                           "Relative error = ", tot_energy_rel_err(icol), &
        3343           0 :                           "Column-integrated total energy at start = ", &
        3344           0 :                           total_energy_column_start(icol), &
        3345           0 :                           "Column-integrated total energy at finish = ", &
        3346           0 :                           total_energy_column_finish(icol)
        3347             :               endif ! tot_energy_rel_err(icol) >= err_thresh
        3348             :            end do ! icol = 1, ncol
        3349             :         endif ! any( tot_energy_rel_err >= err_thresh )
        3350             :
        3351             :      endif ! l_check_conservation
        3352             :
        3353             :
        3354           0 :      return
        3355             :
        3356           0 :    end subroutine subcol_SILHS_fill_holes_conserv
        3357             :
        3358             :    !============================================================================
        3359           0 :    subroutine fill_holes_sedimentation( dt, ncol, hm_start, pdel, &
        3360             :                                         fallspeed_m_per_s, zi, qmin_hm, &
        3361             :                                         hm_tend, prec )
        3362             :
        3363             :      ! Description:
        3364             :      ! After hydrometeor tendencies from microphysics processes were adjusted
        3365             :      ! so that holes don't form in a hydrometeor field from microphysics
        3366             :      ! processes, the sedimentation tendency was added back on to produce an
        3367             :      ! updated total microphysics tendency.  The first-order "up-gravity"
        3368             :      ! sedimentation method that was originally used is positive definite.
        3369             :      ! However, since the microphysics process tendencies were altered so that
        3370             :      ! holes don't form, it is possible that adding the old sedimentation
        3371             :      ! tendencies back onto the new microphysics process tendencies could
        3372             :      ! produce new total microphysics tendencies that cause holes to form.
        3373             :      !
        3374             :      ! In this subroutine, holes in a hydrometeor field are checked for after
        3375             :      ! the updated microphysics tendency is applied.  If any are found, they are
        3376             :      ! filled from positive hydrometeor mass found at grid levels below where
        3377             :      ! the hole is found.  The levels that are used to fill are within range
        3378             :      ! based on fallspeed of the hydrometeor.  If the level that contains the
        3379             :      ! hole is within proximity to the surface, then the water that sedimented
        3380             :      ! to the surface can be used to fill the hole, as well.
        3381             :      !
        3382             :      ! If there isn't enough total hydrometeor mass within the fall range to
        3383             :      ! fill the hole, then positive hydrometeor mass from levels below the fall
        3384             :      ! range is to be added to the total available mass to fill the hole.  Mass
        3385             :      ! is added one level at a time until enough mass is found to fill the hole
        3386             :      ! or until the surface is reached and the surface precipitation is added to
        3387             :      ! the total available fill mass.
        3388             :      !
        3389             :      ! After this, if there still isn't enough available mass to fill the hole,
        3390             :      ! then positive hydrometeor mass is added from all levels above the hole to
        3391             :      ! the total mass that is available to fill the hole.
        3392             :
        3393             :      !----------------------------------------------------------------------
        3394             :
        3395           0 :      use ppgrid, only: &
        3396             :          pcols
        3397             :
        3398             :      use ref_pres, only: &
        3399             :          top_lev => trop_cloud_top_lev
        3400             :
        3401             :      implicit none
        3402             :
        3403             :      ! Input Variables
        3404             :      real(r8), intent(in) :: dt                   ! Time step duration
        3405             :
        3406             :      integer, intent(in) :: ncol                  ! Number of grid columns
        3407             :
        3408             :      real(r8), dimension(pcols,pver), intent(in) :: &
        3409             :        hm_start, & ! Hydrometeor mixing ratio at start of microphysics  [kg/kg]
        3410             :        pdel        ! Pressure difference between grid levels            [Pa]
        3411             :
        3412             :      real(r8), dimension(pcols,pver), intent(in) :: &
        3413             :        fallspeed_m_per_s    ! Hydrometeor mixing ratio fall speed     [m/s]
        3414             :
        3415             :      real(r8), dimension(pcols,pverp), intent(in) :: &
        3416             :        zi    ! Height of momentum (interface) grid levels    [m]
        3417             :
        3418             :      real(r8), intent(in) :: &
        3419             :        qmin_hm    ! Minimum threshold value of hydrometeor mixing ratio  [kg/kg]
        3420             :
        3421             :      ! Input/Output Variables
        3422             :      real(r8), dimension(pcols,pver), intent(inout) :: &
        3423             :        hm_tend    ! Hydrometeor mixing ratio tendency  [kg/kg/s]
        3424             :
        3425             :      real(r8), dimension(pcols), intent(inout) :: &
        3426             :        prec    ! Precipitation rate (surface)        [m/s]
        3427             :
        3428             :      ! Local Variables
        3429             :      real(r8), dimension(pver) :: &
        3430             :        hm_update, & ! Hydrometeor mixing ratio; start of sed. hole fill  [kg/kg]
        3431             :        hm_curr      ! Current value of hydrometeor mixing ratio          [kg/kg]
        3432             :
        3433             :      real(r8) :: &
        3434             :        total_hole,          & ! Total mass of hole in hydrometeor       [kg/m^2]
        3435             :        total_fill_mass,     & ! Total mass available to fill hole       [kg/m^2]
        3436             :        hole_fillmass_ratio, & ! Ratio: total_hole / total_fill_mass     [-]
        3437             :        fallspeed_Pa_per_s,  & ! Hydrometeor mixing ratio fall speed     [Pa/s]
        3438             :        total_fall_Pa,       & ! Pressure "distance" hydrometeor fell    [Pa]
        3439             :        sum_pdel               ! Sum of pdel over levels                 [Pa]
        3440             :
        3441             :      logical, dimension(pver) :: &
        3442             :        l_pos_hm  ! Flag for a hydrometeor having a positive (>= qmin_hm) value
        3443             :
        3444             :      ! Flag for whether surface precipitation mass needs to be included in
        3445             :      ! the total_fill_mass for hole filling.
        3446             :      logical :: l_reached_surface
        3447             :
        3448             :      ! Flag for whether hydrometeor mass from levels above the hole needs to be
        3449             :      ! included in the total_fill_mass for hole filling.
        3450             :      logical :: l_fill_from_above
        3451             :
        3452             :      integer :: icol  ! Grid column index
        3453             :
        3454             :      integer :: k, idx  ! Vertical grid level indices
        3455             :
        3456             :      ! Index of the lowest vertical grid level that needs to be included in the
        3457             :      ! total_fill_mass for hole filling.
        3458             :      integer :: lowest_level_idx
        3459             :
        3460             :
        3461             :      ! Loop over all columns, performing any adjustments one column at a time.
        3462           0 :      do icol = 1, ncol
        3463             :
        3464             :         ! Calculate the updated value of the hydrometeor field based on the
        3465             :         ! updated microphysics tendency.  Since the original sedimentation
        3466             :         ! tendency has been added to the updated microphysics process tendency
        3467             :         ! to produce the updated total microphysics tendency (hm_tend), the
        3468             :         ! updated value of the hydrometeor field (hm_update) could be negative.
        3469           0 :         hm_update = hm_start(icol,:) + hm_tend(icol,:) * dt
        3470           0 :         hm_curr = hm_update
        3471             :
        3472             :         ! Check for any holes in the vertical profile
        3473           0 :         if ( any( hm_curr(top_lev:pver) < qmin_hm ) ) then
        3474             :
        3475             :            ! At least one hole is found in this hydrometeor species in this
        3476             :            ! grid column.  The holes must be filled conservatively.
        3477             :
        3478             :            ! Check which levels have values of the hydrometeor that are at or
        3479             :            ! above the minimum threshold value.
        3480           0 :            do k = top_lev, pver
        3481           0 :               if ( hm_curr(k) >= qmin_hm ) then
        3482           0 :                 l_pos_hm(k) = .true.
        3483             :               else ! hm_curr < qmin_hm
        3484           0 :                 l_pos_hm(k) = .false.
        3485             :               endif ! hm_curr >= qmin_hm
        3486             :            end do ! k = top_lev, pver
        3487             :
        3488           0 :            do k = pver, top_lev, -1
        3489             :
        3490           0 :               if ( .not. l_pos_hm(k) ) then
        3491             :
        3492             :                  ! A hole is found in the hydrometeor at this grid level.
        3493             :
        3494             :                  ! Calculate the total hydrometeor mass of the hole that needs
        3495             :                  ! to be filled.
        3496             :                  ! The value of the hydrometeor mixing ratio is negative, but
        3497             :                  ! the value of total_hole is positive.
        3498           0 :                  total_hole = ( qmin_hm - hm_curr(k) ) * pdel(icol,k) * rga
        3499             :
        3500             :                  ! Calculate the total hydrometeor mass available from below
        3501             :                  ! to fill the hole.
        3502           0 :                  if ( k == pver ) then
        3503             :
        3504             :                     ! A hole is found at the lowermost level.
        3505             :                     ! The only place the hydrometeor could have sedimented
        3506             :                     ! to is the surface, so fill from only the surface.
        3507           0 :                     l_reached_surface = .true.
        3508             :
        3509             :                     ! Calculate the available amount of hydrometeor mass to
        3510             :                     ! fill the hole.
        3511           0 :                     total_fill_mass = prec(icol) * dt * 1000.0_r8
        3512             :
        3513             :                  else ! top_lev <= k < pver
        3514             :
        3515             :                     ! Calculate the hydrometeor fallspeed in Pa/s.
        3516             :                     ! In MG2, the equation for this is given by:
        3517             :                     !
        3518             :                     ! fallspeed([Pa/s]) = g * rho * fallspeed([m/s]).
        3519             :                     !
        3520             :                     ! The value of rho is typically calculated from the
        3521             :                     ! hydrostatic approximation:
        3522             :                     !
        3523             :                     ! rho = - ( 1 / g ) * dp/dz.
        3524             :                     !
        3525             :                     ! The equation for fallspeed in Pa/s becomes:
        3526             :                     !
        3527             :                     ! fallspeed([Pa/s]) = - dp/dz * fallspeed([m/s]).
        3528             :                     fallspeed_Pa_per_s &
        3529             :                     = fallspeed_m_per_s(icol,k) &
        3530           0 :                       * pdel(icol,k) / ( zi(icol,k) - zi(icol,k+1) )
        3531             :
        3532             :                     ! Calculate the fall "distance" in Pa.
        3533           0 :                     total_fall_Pa = fallspeed_Pa_per_s * dt
        3534             :
        3535             :                     ! Find the index of the vertical level that the hydrometeor
        3536             :                     ! sedimented to in one timestep.  It must sediment at least
        3537             :                     ! one level.
        3538           0 :                     sum_pdel = 0.0_r8
        3539           0 :                     idx = k + 1
        3540           0 :                     do
        3541             :                        ! Update the total pressure difference between the
        3542             :                        ! level of origin and the current level.
        3543           0 :                        sum_pdel = sum_pdel + pdel(icol,idx)
        3544           0 :                        if ( sum_pdel >= total_fall_Pa ) then
        3545             :                           ! The total pressure difference between the level of
        3546             :                           ! origin and the current level exceeds the total
        3547             :                           ! hydrometeor fall "distance" (in Pa).
        3548             :                           lowest_level_idx = idx
        3549             :                           l_reached_surface = .false.
        3550             :                           exit
        3551             :                        else ! sum_pdel < total_fall_Pa
        3552             :                           ! The total hydrometeor fall "distance" (in Pa)
        3553             :                           ! exceeds the total pressure difference between the
        3554             :                           ! level of origin and the current level.
        3555           0 :                           if ( idx == pver ) then
        3556             :                              ! The lowest level of the model has been reached.
        3557             :                              ! The hydrometeor sedimented to the surface.
        3558             :                              lowest_level_idx = pver
        3559             :                              l_reached_surface = .true.
        3560             :                              exit
        3561             :                           else ! idx < pver
        3562             :                              ! Increment idx and keep going.
        3563           0 :                              idx = idx + 1
        3564             :                           endif ! idx == pver
        3565             :                        endif ! sum_pdel >= total_fall_Pa
        3566             :                     end do
        3567             :
        3568             :                     ! Calculate the available amount of hydrometeor mass to
        3569             :                     ! fill the hole.
        3570             :                     total_fill_mass = 0.0_r8
        3571             :                     if ( l_reached_surface ) then
        3572             :                        ! The hydrometeor sedimented to the surface, so
        3573             :                        ! automatically loop down to pver and include the
        3574             :                        ! surface mass.
        3575           0 :                        do idx = k+1, pver, 1
        3576           0 :                           if ( l_pos_hm(idx) ) then
        3577             :                              total_fill_mass &
        3578             :                              = total_fill_mass &
        3579             :                                + ( hm_curr(idx) - qmin_hm ) &
        3580           0 :                                  * pdel(icol,idx) * rga
        3581             :                           endif ! l_pos_hm(idx)
        3582             :                        end do ! idx = k+1, pver, 1
        3583             :                        ! Contribution to total fill mass from the surface.
        3584             :                        total_fill_mass &
        3585           0 :                        = total_fill_mass + prec(icol) * dt * 1000.0_r8
        3586             :                     else ! .not. l_reached_surface
        3587             :                        ! The hydrometeor sedimented to lowest_level_idx.
        3588             :                        idx = k + 1
        3589             :                        do
        3590           0 :                           if ( l_pos_hm(idx) ) then
        3591             :                              total_fill_mass &
        3592             :                              = total_fill_mass &
        3593             :                                + ( hm_curr(idx) - qmin_hm ) &
        3594           0 :                                  * pdel(icol,idx) * rga
        3595             :                           endif ! l_pos_hm(idx)
        3596           0 :                           if ( idx >= lowest_level_idx ) then
        3597             :                              ! Check if enough mass has been gathered in
        3598             :                              ! total_fill_mass to fill the hole.
        3599           0 :                              if ( total_fill_mass >= total_hole ) then
        3600             :                                 ! There has been enough total_fill_mass
        3601             :                                 ! gathered to completely fill the hole.
        3602             :                                 lowest_level_idx = idx
        3603             :                                 exit
        3604             :                              else ! total_fill_mass < total_hole
        3605             :                                 ! Even though lowest_level_idx has been reached,
        3606             :                                 ! more total_fill_mass needs to be added in
        3607             :                                 ! order to completely fill the hole, so keep
        3608             :                                 ! going.
        3609           0 :                                 if ( idx == pver ) then
        3610             :                                    ! The lowest vertical level has already been
        3611             :                                    ! reached, so go to the surface.
        3612           0 :                                    lowest_level_idx = pver
        3613           0 :                                    l_reached_surface = .true.
        3614             :                                    ! Contribution to total fill mass from the
        3615             :                                    ! surface.
        3616             :                                    total_fill_mass &
        3617             :                                    = total_fill_mass &
        3618           0 :                                      + prec(icol) * dt * 1000.0_r8
        3619           0 :                                    exit
        3620             :                                 else ! idx < pver
        3621             :                                    ! Haven't reached pver yet, so increment
        3622             :                                    ! and keep going.
        3623           0 :                                    idx = idx + 1
        3624             :                                 endif ! idx == pver
        3625             :                              endif ! total_fill_mass >= total_hole
        3626             :                           else ! idx < lowest_level_idx
        3627             :                              ! Haven't reached lowest_level_idx yet, so
        3628             :                              ! increment and keep going.
        3629           0 :                              idx = idx + 1
        3630             :                           endif ! idx >= lowest_level_idx
        3631             :                        end do
        3632             :                     endif ! l_reached_surface
        3633             :
        3634             :                  endif ! k == pver
        3635             :
        3636             :                  ! If mass has been added all the way down to the surface and
        3637             :                  ! there's still not enough mass to fill the hole, then fill the
        3638             :                  ! hole pulling mass from above.
        3639           0 :                  if ( total_fill_mass >= total_hole ) then
        3640             :                     l_fill_from_above = .false.
        3641             :                  else ! total_fill_mass < total_hole
        3642           0 :                     l_fill_from_above = .true.
        3643           0 :                     do idx = top_lev, k-1, 1
        3644           0 :                        if ( l_pos_hm(idx) ) then
        3645             :                           total_fill_mass &
        3646             :                           = total_fill_mass &
        3647             :                             + ( hm_curr(idx) - qmin_hm ) &
        3648           0 :                               * pdel(icol,idx) * rga
        3649             :                        endif ! l_pos_hm(idx)
        3650             :                     end do ! idx = top_lev, k-1, 1
        3651             :                  endif ! total_fill_mass >= total_hole
        3652             :
        3653             :                  ! Calculate the ratio of total hole to total fill mass.  This
        3654             :                  ! should not exceed 1 except as a result of numerical round-off
        3655             :                  ! errors.  Use thresholding to be safe.
        3656             :                  hole_fillmass_ratio &
        3657             :                  = min( total_hole / max( total_fill_mass, 1.0e-30_r8 ), &
        3658           0 :                         1.0_r8 )
        3659             :
        3660           0 :                  if ( k < pver ) then
        3661             :                     ! Modify (reduce) the amount of the hydrometeor at levels
        3662             :                     ! that were used to fill the hole.
        3663           0 :                     do idx = k+1, lowest_level_idx
        3664           0 :                        if ( l_pos_hm(idx) ) then
        3665             :                           ! Since pdel at a grid level does not change and
        3666             :                           ! gravit is constant, the only variable that needs to
        3667             :                           ! be modified proportionately is hm_curr.
        3668             :                           hm_curr(idx) &
        3669             :                           = qmin_hm &
        3670             :                             + ( hm_curr(idx) - qmin_hm ) &
        3671           0 :                               * ( 1.0_r8 - hole_fillmass_ratio )
        3672             :                        endif ! l_pos_hm(idx)
        3673             :                     end do ! idx = k+1, lowest_level_idx
        3674             :                  endif ! k < pver
        3675             :
        3676           0 :                  if ( l_reached_surface ) then
        3677             :                     ! Modify (reduce) the amount of surface precipitation in
        3678             :                     ! order to fill the hole.  Since dt and 1000 are constants,
        3679             :                     ! the only variable that needs to be modified
        3680             :                     ! proportionately is prec.
        3681           0 :                     prec(icol) = prec(icol) * ( 1.0_r8 - hole_fillmass_ratio )
        3682             :                  endif ! l_reached_surface
        3683             :
        3684           0 :                  if ( l_fill_from_above ) then
        3685             :                     ! Modify (reduce) the amount of the hydrometeor at levels
        3686             :                     ! that were used to fill the hole.
        3687           0 :                     do idx = top_lev, k-1
        3688           0 :                        if ( l_pos_hm(idx) ) then
        3689             :                           ! Since pdel at a grid level does not change and
        3690             :                           ! gravit is constant, the only variable that needs to
        3691             :                           ! be modified proportionately is hm_curr.
        3692             :                           hm_curr(idx) &
        3693             :                           = qmin_hm &
        3694             :                             + ( hm_curr(idx) - qmin_hm ) &
        3695           0 :                               * ( 1.0_r8 - hole_fillmass_ratio )
        3696             :                        endif ! l_pos_hm(idx)
        3697             :                     end do ! idx = top_lev, k-1
        3698             :                  endif ! l_fill_from_above
        3699             :
        3700             :                  ! Update the value of the hydrometeor at the level where the
        3701             :                  ! hole was found.  Mathematically, as long as the available
        3702             :                  ! mass was able to fill the entire hole, the new value of the
        3703             :                  ! hydrometeor mixing ratio (hm_curr) should be qmin_hm.
        3704             :                  hm_curr(k) &
        3705             :                  = hm_curr(k) &
        3706             :                    + hole_fillmass_ratio * total_fill_mass &
        3707           0 :                      * gravit / pdel(icol,k)
        3708             :
        3709             :               endif ! .not. l_pos_hm(k)
        3710             :
        3711             :            end do ! k = pver, top_lev, -1
        3712             :
        3713             :         endif ! any( hm_curr(top_lev:pver) < qmin_hm )
        3714             :
        3715             :         ! Update the value of total microphysics tendency after hole filling.
        3716           0 :         hm_tend(icol,:) = hm_tend(icol,:) + ( hm_curr - hm_update ) / dt
        3717             :
        3718             :      end do ! icol = 1, ncol
        3719             :
        3720             :
        3721           0 :      return
        3722             :
        3723           0 :    end subroutine fill_holes_sedimentation
        3724             :
        3725             :    !============================================================================
        3726           0 :    subroutine fill_holes_same_phase_vert( dt, ncol, hm_start, hm_start_filler, &
        3727             :                                           pdel, qmin_hm, qmin_hm_filler, &
        3728             :                                           hm_tend, hm_tend_filler )
        3729             :
        3730             :      ! Description:
        3731             :      ! Fills remaining holes in a hydrometeor with mass from the the vertical
        3732             :      ! profile of another hydrometeor of the same phase.  Remaining holes in
        3733             :      ! cloud water are filled with rain water and remaining holes in snow are
        3734             :      ! filled with cloud ice.
        3735             :      !
        3736             :      ! This subroutine, combined with subroutine fill_holes_sedimentation, fill
        3737             :      ! holes making use of the following principles:
        3738             :      !
        3739             :      ! 3) When adjusting the hydrometeor tendency from sedimentation of a
        3740             :      !    liquid hydrometeor (cloud water or rain water), conserve:
        3741             :      !
        3742             :      !    SUM(k=top_lev:pver) ( rc_sed_tend(k) + rr_sed_tend(k) )
        3743             :      !                        * dt * pdel(k) / g
        3744             :      !    + precl * dt * 1000 = 0.
        3745             :      !
        3746             :      ! 4) When adjusting the hydrometeor tendency from sedimentation of a
        3747             :      !    frozen hydrometeor (cloud ice or snow), conserve:
        3748             :      !
        3749             :      !    SUM(k=top_lev:pver) ( ri_sed_tend(k) + rs_sed_tend(k) )
        3750             :      !                        * dt * pdel(k) / g
        3751             :      !    + preci * dt * 1000 = 0.
        3752             :      !
        3753             :      ! These two equations (one for liquid-phase hydrometeors and one for
        3754             :      ! ice-phase hydrometeors) could be further split into one equation for
        3755             :      ! each hydrometeor if there was prec output for each hydrometeor.  However,
        3756             :      ! there's only prec output for ice-phase precipitation rate and total
        3757             :      ! precipitation rate (liquid preciptation rate is total rate minus
        3758             :      ! ice-phase rate).
        3759             :      !
        3760             :      ! Since only liquid-phase precipitation rate (precl) and ice-phase
        3761             :      ! precipitation rate (preci) are available, and there are two hydrometeors
        3762             :      ! in each category, one hydrometeor from each category must fill before
        3763             :      ! the other hydrometeor from its category and get priority access to precl
        3764             :      ! or preci.  Since a vast majority of liquid precipitation comes from rain
        3765             :      ! rather than sedimenting cloud water, rain is filled before cloud water
        3766             :      ! and gets priority access to precl.  Likewise, since a vast majority of
        3767             :      ! frozen precipitation comes from snow rather than sedimenting cloud ice,
        3768             :      ! snow is filled before cloud ice and gets priority access to preci.
        3769             :      !
        3770             :      ! The order of sedimentation hole filling is as follows.  First, a level
        3771             :      ! with a hole in it is identified.  The fall distance for the hydrometeor
        3772             :      ! that originated at a level is calculated.  Total mass to fill the hole is
        3773             :      ! calculated from all levels within the fall range that have positive
        3774             :      ! values of the hydrometeor.  The amount that precipitated to the surface
        3775             :      ! is also included if the hydrometeor fell that far.  If that isn't enough
        3776             :      ! mass to fill the hole, then levels that are lower in the profile are
        3777             :      ! included (if the hydrometeor has a positive value) until enough mass is
        3778             :      ! found to fill the hole or until the surface is reached.  If there isn't
        3779             :      ! enough mass found in all levels below the hole, including the amount that
        3780             :      ! precipitated to the ground, to fill the hole, then the hydrometeor mass
        3781             :      ! from all levels above the hole (again, where a positive value of the
        3782             :      ! hydrometeor is found) are included in the total available mass to fill
        3783             :      ! the hole.
        3784             :      !
        3785             :      ! Occasionally, a situation can occur where both hydrometeors in a category
        3786             :      ! contributed to surface precipitation rate, and filling a hole in rain
        3787             :      ! (or snow) can deplete all the surface precl (or preci), resulting in not
        3788             :      ! enough water mass in the vertical profile (including the surface) of
        3789             :      ! cloud water (or cloud ice) to fill a hole in cloud water (or cloud ice).
        3790             :      ! When this happens, there must still be liquid water (or frozen water)
        3791             :      ! found in the vertical profile of rain (or snow), so pull the water from
        3792             :      ! rain (or snow) to fill any remaining holes in cloud water (or cloud ice).
        3793             :
        3794             :      !----------------------------------------------------------------------
        3795             :
        3796           0 :      use ppgrid, only: &
        3797             :          pcols
        3798             :
        3799             :      use ref_pres, only: &
        3800             :          top_lev => trop_cloud_top_lev
        3801             :
        3802             :      implicit none
        3803             :
        3804             :      ! Input Variables
        3805             :      real(r8), intent(in) :: dt                   ! Time step duration
        3806             :
        3807             :      integer, intent(in) :: ncol                  ! Number of grid columns
        3808             :
        3809             :      real(r8), dimension(pcols,pver), intent(in) :: &
        3810             :        hm_start,        & ! Hydrometeor mixing ratio (microphys start)   [kg/kg]
        3811             :        hm_start_filler, & ! Filler hydromet mix ratio (microphys start)  [kg/kg]
        3812             :        pdel               ! Pressure difference between grid levels      [Pa]
        3813             :
        3814             :      real(r8), intent(in) :: &
        3815             :        qmin_hm,        & ! Minimum threshold hydrometeor mixing ratio  [kg/kg]
        3816             :        qmin_hm_filler    ! Min threshold filler hydromet mixing ratio  [kg/kg]
        3817             :
        3818             :      ! Input/Output Variables
        3819             :      real(r8), dimension(pcols,pver), intent(inout) :: &
        3820             :        hm_tend,        & ! Hydrometeor mixing ratio tendency         [kg/kg/s]
        3821             :        hm_tend_filler    ! Filler hydrometeor mixing ratio tendency  [kg/kg/s]
        3822             :
        3823             :      ! Local Variables
        3824             :      real(r8), dimension(pver) :: &
        3825             :        hm_update,        & ! Hydrometeor mixing ratio; start           [kg/kg]
        3826             :        hm_update_filler, & ! Filler Hydrometeor mixing ratio; start    [kg/kg]
        3827             :        hm_curr,          & ! Current hydrometeor mixing ratio          [kg/kg]
        3828             :        hm_curr_filler      ! Current filler hydrometeor mixing ratio   [kg/kg]
        3829             :
        3830             :      real(r8) :: &
        3831             :        total_hole,          & ! Total mass of hole in hydrometeor       [kg/m^2]
        3832             :        total_fill_mass,     & ! Total mass available to fill hole       [kg/m^2]
        3833             :        hole_fillmass_ratio    ! Ratio: total_hole / total_fill_mass     [-]
        3834             :
        3835             :      logical, dimension(pver) :: &
        3836             :        l_pos_hm,        & ! Flag: hydrometeor has positive (>= qmin_hm) value
        3837             :        l_pos_hm_filler    ! Flag: filler hydrometeor has positive value
        3838             :
        3839             :      integer :: icol  ! Grid column index
        3840             :
        3841             :      integer :: k, idx  ! Vertical grid level indices
        3842             :
        3843             :
        3844             :      ! Loop over all columns, performing any adjustments one column at a time.
        3845           0 :      do icol = 1, ncol
        3846             :
        3847             :         ! Calculate the updated value of the hydrometeor field based on the
        3848             :         ! updated microphysics tendency.
        3849           0 :         hm_update = hm_start(icol,:) + hm_tend(icol,:) * dt
        3850           0 :         hm_curr = hm_update
        3851             :
        3852             :         ! Calculate the updated value of the filler hydrometeor field based on
        3853             :         ! the updated microphysics tendency.
        3854           0 :         hm_update_filler = hm_start_filler(icol,:) + hm_tend_filler(icol,:) * dt
        3855           0 :         hm_curr_filler = hm_update_filler
        3856             :
        3857             :         ! Check for any holes in the vertical profile
        3858           0 :         if ( any( hm_curr(top_lev:pver) < qmin_hm ) ) then
        3859             :
        3860             :            ! At least one hole is found in this hydrometeor species in this
        3861             :            ! grid column.  The holes must be filled conservatively.
        3862             :
        3863             :            ! Check which levels have values of the hydrometeor that are at or
        3864             :            ! above the minimum threshold value.
        3865           0 :            do k = top_lev, pver
        3866             :               ! Check for the hydrometeor that might need to be filled.
        3867           0 :               if ( hm_curr(k) >= qmin_hm ) then
        3868           0 :                 l_pos_hm(k) = .true.
        3869             :               else ! hm_curr < qmin_hm
        3870           0 :                 l_pos_hm(k) = .false.
        3871             :               endif ! hm_curr >= qmin_hm
        3872             :               ! Check for the filler hydrometeor, as some levels might have
        3873             :               ! numerical round-off level, small negative values.
        3874           0 :               if ( hm_curr_filler(k) >= qmin_hm_filler ) then
        3875           0 :                 l_pos_hm_filler(k) = .true.
        3876             :               else ! hm_curr_filler < qmin_hm_filler
        3877           0 :                 l_pos_hm_filler(k) = .false.
        3878             :               endif ! hm_curr_filler >= qmin_hm_filler
        3879             :            end do ! k = top_lev, pver
        3880             :
        3881           0 :            do k = top_lev, pver
        3882             :
        3883           0 :               if ( .not. l_pos_hm(k) ) then
        3884             :
        3885             :                  ! A hole is found in the hydrometeor at this grid level.
        3886             :
        3887             :                  ! Calculate the total hydrometeor mass of the hole that needs
        3888             :                  ! to be filled.
        3889             :                  ! The value of the hydrometeor mixing ratio is negative, but
        3890             :                  ! the value of total_hole is positive.
        3891           0 :                  total_hole = ( qmin_hm - hm_curr(k) ) * pdel(icol,k) * rga
        3892             :
        3893             :                  ! Calculate the total hydrometeor mass available from the
        3894             :                  ! filler hydrometeor to fill the hole.
        3895           0 :                  total_fill_mass = 0.0_r8
        3896           0 :                  do idx = top_lev, pver, 1
        3897           0 :                     if ( l_pos_hm_filler(idx) ) then
        3898             :                         total_fill_mass &
        3899             :                         = total_fill_mass &
        3900             :                           + ( hm_curr_filler(idx) - qmin_hm_filler ) &
        3901           0 :                             * pdel(icol,idx) * rga
        3902             :                      endif ! l_pos_hm_filler(idx)
        3903             :                  end do ! idx = top_lev, pver, 1
        3904             :
        3905             :                  ! Calculate the ratio of total hole to total fill mass.  This
        3906             :                  ! should not exceed 1 except as a result of numerical round-off
        3907             :                  ! errors.  Use thresholding to be safe.
        3908             :                  hole_fillmass_ratio &
        3909             :                  = min( total_hole / max( total_fill_mass, 1.0e-30_r8 ), &
        3910           0 :                         1.0_r8 )
        3911             :
        3912             :                  ! Modify (reduce) the amount of the filler hydrometeor.
        3913           0 :                  do idx = top_lev, pver
        3914           0 :                     if ( l_pos_hm_filler(idx) ) then
        3915             :                        ! Since pdel at a grid level does not change and gravit
        3916             :                        ! is constant, the only variable that needs to be
        3917             :                        ! modified proportionately is hm_curr_filler.
        3918             :                        hm_curr_filler(idx) &
        3919             :                        = qmin_hm_filler &
        3920             :                          + ( hm_curr_filler(idx) - qmin_hm_filler ) &
        3921           0 :                            * ( 1.0_r8 - hole_fillmass_ratio )
        3922             :                     endif ! l_pos_hm_filler(idx)
        3923             :                  end do ! idx = top_lev, pver
        3924             :
        3925             :                  ! Update the value of the hydrometeor at the level where the
        3926             :                  ! hole was found.  Mathematically, as long as the available
        3927             :                  ! mass was able to fill the entire hole, the new value of the
        3928             :                  ! hydrometeor mixing ratio (hm_curr) should be qmin_hm.
        3929             :                  hm_curr(k) &
        3930             :                  = hm_curr(k) &
        3931             :                    + hole_fillmass_ratio * total_fill_mass &
        3932           0 :                      * gravit / pdel(icol,k)
        3933             :
        3934             :               endif ! .not. l_pos_hm(k)
        3935             :
        3936             :            end do ! k = top_lev, pver
        3937             :
        3938             :         endif ! any( hm_curr(top_lev:pver) < qmin_hm )
        3939             :
        3940             :         ! Update the value of total microphysics tendency after hole filling.
        3941           0 :         hm_tend(icol,:) = hm_tend(icol,:) + ( hm_curr - hm_update ) / dt
        3942             :
        3943             :         ! Update the value of total microphysics tendency after hole filling for
        3944             :         ! the filler hydrometeor.
        3945             :         hm_tend_filler(icol,:) &
        3946           0 :         = hm_tend_filler(icol,:) + ( hm_curr_filler - hm_update_filler ) / dt
        3947             :
        3948             :      end do ! icol = 1, ncol
        3949             :
        3950             :
        3951           0 :      return
        3952             :
        3953           0 :    end subroutine fill_holes_same_phase_vert
        3954             :
        3955             :    !============================================================================
        3956           0 :    subroutine subcol_SILHS_hydromet_conc_tend_lim( state, dt, ptend )
        3957             :
        3958             :      ! Description:
        3959             :      ! Limits the values of mean hydrometeor concentrations so that the mean
        3960             :      ! drop size for the hydrometeor type remains reasonable and does not become
        3961             :      ! too large.
        3962             :
        3963             :      !----------------------------------------------------------------------
        3964             :
        3965           0 :      use shr_const_mod, only: &
        3966             :          shr_const_pi, &
        3967             :          shr_const_rhofw
        3968             :
        3969             :      use constituents, only: &
        3970             :          qmin
        3971             :
        3972             :      use ref_pres, only: &
        3973             :          top_lev => trop_cloud_top_lev
        3974             :
        3975             :      implicit none
        3976             :
        3977             :      ! Input Variables
        3978             :      type(physics_state), intent(in) :: state     ! Physics state variables
        3979             :      real(r8), intent(in) :: dt                   ! Time step duration
        3980             :
        3981             :      ! Input/Output Variable
        3982             :      type(physics_ptend),  intent(inout) :: ptend  ! Parameterization tendencies
        3983             :
        3984             :      ! Local Variables
        3985             :      real( r8 ) :: &
        3986             :        rcm_update, & ! New value of mean cloud water mixing ratio    [kg/kg]
        3987             :        rrm_update, & ! New value of mean rain water mixing ratio     [kg/kg]
        3988             :        rim_update, & ! New value of mean ice mixing ratio            [kg/kg]
        3989             :        rsm_update    ! New value of mean snow mixing ratio           [kg/kg]
        3990             :
        3991             :      real( r8 ) :: &
        3992             :        Nc_tend_min, & ! Minimum value of cloud droplet conc. tendency [num/kg/s]
        3993             :        Nr_tend_min, & ! Minimum value of rain drop conc. tendency     [num/kg/s]
        3994             :        Ni_tend_min, & ! Minimum value of ice conc. tendency           [num/kg/s]
        3995             :        Ns_tend_min    ! Minimum value of snow conc. tendency          [num/kg/s]
        3996             :
        3997             :      real( r8 ), parameter :: &
        3998             :        four_thirds = 4.0_r8/3.0_r8,    & ! 4/3
        3999             :        rho_ice     = 917.0_r8,         & ! Density of ice            [kg/m^3]
        4000             :        rho_lw      = shr_const_rhofw,  & ! Density of liquid water   [kg/m^3]
        4001             :        pi          = shr_const_pi        ! Pi
        4002             :
        4003             :      real( r8 ), parameter :: &
        4004             :        mvr_cloud_max = 1.6E-5_r8, & ! Max. avg. mean vol. rad. cloud   [m]
        4005             :        mvr_rain_max  = 5.0E-3_r8, & ! Max. avg. mean vol. rad. rain    [m]
        4006             :        mvr_ice_max   = 1.3E-4_r8, & ! Max. avg. mean vol. rad. ice     [m]
        4007             :        mvr_snow_max  = 1.0E-2_r8    ! Max. avg. mean vol. rad. snow    [m]
        4008             :
        4009             :      ! Calculate the coefficient for the minimum mean cloud droplet
        4010             :      ! concentration, where <Nc>|_min = Ncm_min_coef * <rc> and has units of
        4011             :      ! 1/kg.
        4012             :      real( r8 ), parameter :: &
        4013             :        Ncm_min_coef = 1.0_r8 / ( four_thirds * pi * rho_lw * mvr_cloud_max**3 )
        4014             :
        4015             :      ! Calculate the coefficient for the minimum mean rain drop concentration,
        4016             :      ! where <Nr>|_min = Nrm_min_coef * <rr> and has units of 1/kg.
        4017             :      real( r8 ), parameter :: &
        4018             :        Nrm_min_coef = 1.0_r8 / ( four_thirds * pi * rho_lw * mvr_rain_max**3 )
        4019             :
        4020             :      ! Calculate the coefficient for the minimum mean ice crystal concentration,
        4021             :      ! where <Ni>|_min = Nim_min_coef * <ri> and has units of 1/kg.
        4022             :      real( r8 ), parameter :: &
        4023             :        Nim_min_coef = 1.0_r8 / ( four_thirds * pi * rho_ice * mvr_ice_max**3 )
        4024             :
        4025             :      ! Calculate the coefficient for the minimum mean snow flake concentration,
        4026             :      ! where <Ns>|_min = Nsm_min_coef * <rs> and has units of 1/kg.
        4027             :      real( r8 ), parameter :: &
        4028             :        Nsm_min_coef = 1.0_r8 / ( four_thirds * pi * rho_ice * mvr_snow_max**3 )
        4029             :
        4030             :      integer :: ncol  ! Number of grid columns
        4031             :
        4032             :      integer :: icol    ! Column loop index
        4033             :
        4034             :      integer :: k    ! Vertical level loop index
        4035             :
        4036             :
        4037             :      ! Get the number of grid columns.
        4038           0 :      ncol = state%ncol
        4039             :
        4040             :      ! Loop over all grid columns.
        4041           0 :      do icol = 1, ncol
        4042             :
        4043             :         ! Loop over all vertical levels from top_lev to pver.
        4044           0 :         do k = top_lev, pver
        4045             :
        4046             :            ! Cloud droplet concentration
        4047           0 :            if ( ixcldliq > 0 .and. ixnumliq > 0 ) then
        4048             :
        4049             :               ! Calculate the value of cloud water mixing ratio after the
        4050             :               ! update.
        4051             :               rcm_update &
        4052           0 :               = max( state%q(icol,k,ixcldliq) + ptend%q(icol,k,ixcldliq) * dt, &
        4053           0 :                      qmin(ixcldliq) )
        4054             :
        4055             :               ! Calculate the limiting cloud droplet concentration tendency so
        4056             :               ! that cloud maintains a reasonable (not too big) mean volume
        4057             :               ! radius.
        4058             :               Nc_tend_min &
        4059           0 :               = ( Ncm_min_coef * rcm_update - state%q(icol,k,ixnumliq) ) / dt
        4060             :
        4061             :               ! The cloud droplet concentration tendency needs to be the greater
        4062             :               ! of the current Nc_tend and Nc_tend_min.
        4063           0 :               ptend%q(icol,k,ixnumliq) &
        4064           0 :               = max( ptend%q(icol,k,ixnumliq), Nc_tend_min )
        4065             :
        4066             :            endif ! ixcldliq > 0 .and. ixnumliq > 0
        4067             :
        4068             :            ! Rain drop concentration
        4069           0 :            if ( ixrain > 0 .and. ixnumrain > 0 ) then
        4070             :
        4071             :               ! Calculate the value of rain water mixing ratio after the update.
        4072             :               rrm_update &
        4073           0 :               = max( state%q(icol,k,ixrain) + ptend%q(icol,k,ixrain) * dt, &
        4074           0 :                      qmin(ixrain) )
        4075             :
        4076             :               ! Calculate the limiting rain drop concentration tendency so that
        4077             :               ! rain maintains a reasonable (not too big) mean volume radius.
        4078             :               Nr_tend_min &
        4079           0 :               = ( Nrm_min_coef * rrm_update - state%q(icol,k,ixnumrain) ) / dt
        4080             :
        4081             :               ! The rain drop concentration tendency needs to be the greater of
        4082             :               ! the current Nr_tend and Nr_tend_min.
        4083           0 :               ptend%q(icol,k,ixnumrain) &
        4084           0 :               = max( ptend%q(icol,k,ixnumrain), Nr_tend_min )
        4085             :
        4086             :            endif ! ixrain > 0 .and. ixnumrain > 0
        4087             :
        4088             :            ! Ice crystal concentration
        4089           0 :            if ( ixcldice > 0 .and. ixnumice > 0 ) then
        4090             :
        4091             :               ! Calculate the value of ice mixing ratio after the update.
        4092             :               rim_update &
        4093           0 :               = max( state%q(icol,k,ixcldice) + ptend%q(icol,k,ixcldice) * dt, &
        4094           0 :                      qmin(ixcldice) )
        4095             :
        4096             :               ! Calculate the limiting ice crystal concentration tendency so
        4097             :               ! that ice maintains a reasonable (not too big) mean volume
        4098             :               ! radius.
        4099             :               Ni_tend_min &
        4100           0 :               = ( Nim_min_coef * rim_update - state%q(icol,k,ixnumice) ) / dt
        4101             :
        4102             :               ! The ice crystal concentration tendency needs to be the greater
        4103             :               ! of the current Ni_tend and Ni_tend_min.
        4104           0 :               ptend%q(icol,k,ixnumice) &
        4105           0 :               = max( ptend%q(icol,k,ixnumice), Ni_tend_min )
        4106             :
        4107             :            endif ! ixcldice > 0 .and. ixnumice > 0
        4108             :
        4109             :            ! Snow flake concentration
        4110           0 :            if ( ixsnow > 0 .and. ixnumsnow > 0 ) then
        4111             :
        4112             :               ! Calculate the value of snow mixing ratio after the update.
        4113             :               rsm_update &
        4114           0 :               = max( state%q(icol,k,ixsnow) + ptend%q(icol,k,ixsnow) * dt, &
        4115           0 :                      qmin(ixsnow) )
        4116             :
        4117             :               ! Calculate the limiting snow flake concentration tendency so that
        4118             :               ! snow maintains a reasonable (not too big) mean volume radius.
        4119             :               Ns_tend_min &
        4120           0 :               = ( Nsm_min_coef * rsm_update - state%q(icol,k,ixnumsnow) ) / dt
        4121             :
        4122             :               ! The snow flake concentration tendency needs to be the greater of
        4123             :               ! the current Ns_tend and Ns_tend_min.
        4124           0 :               ptend%q(icol,k,ixnumsnow) &
        4125           0 :               = max( ptend%q(icol,k,ixnumsnow), Ns_tend_min )
        4126             :
        4127             :            endif ! ixsnow > 0 .and. ixnumsnow > 0
        4128             :
        4129             :         end do ! k = top_lev, pver
        4130             :
        4131             :      end do ! icol = 1, ncol
        4132             :
        4133             :
        4134           0 :      return
        4135             :
        4136           0 :    end subroutine subcol_SILHS_hydromet_conc_tend_lim
        4137             :
        4138             :    !============================================================================
        4139             :
        4140             :    ! Getunit and Freeunit are depreciated in Fortran going forward, so this is a
        4141             :    ! small function to get an unused stream identifier to send to setup_corr_varnce_array_api
        4142             :    ! or any other silhs/clubb functions that require a unit number argument
        4143             :    ! This comes directly from the Fortran wiki
        4144             :    integer function newunit(unit)
        4145             :      integer, intent(out), optional :: unit
        4146             :
        4147             :      integer, parameter :: LUN_MIN=10, LUN_MAX=1000
        4148             :      logical :: opened
        4149             :      integer :: lun
        4150             :
        4151             :      newunit=-1
        4152             :      do lun=LUN_MIN,LUN_MAX
        4153             :         inquire(unit=lun,opened=opened)
        4154             :         if (.not. opened) then
        4155             :            newunit=lun
        4156             :            exit
        4157             :         end if
        4158             :      end do
        4159             :      if (present(unit)) unit=newunit
        4160           0 :    end function newunit
        4161             :
        4162             : end module subcol_SILHS

Generated by: LCOV version 1.14