LCOV - code coverage report
Current view: | top level - physics/pumas-frozen - micro_pumas_v1.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 0 | 1661 | 0.0 %
Date: | 2025-01-13 21:54:50 |  | Functions: | 0 | 5 | 0.0 %

              Line data    Source code
           1             : module micro_pumas_v1
           2             : !---------------------------------------------------------------------------------
           3             : ! Purpose:
           4             : !   MG microphysics version 3.0 - Update of MG microphysics with
           5             : !                                 prognostic hail OR graupel.
           6             : !
           7             : ! Author: Andrew Gettelman, Hugh Morrison
           8             : !
           9             : ! Version 3 history: Sep 2016: development begun for hail, graupel
          10             : !
          11             : ! Version 2 history: Sep 2011: Development begun.
          12             : !                    Feb 2013: Added of prognostic precipitation.
          13             : !                    Aug 2015: Published and released version
          14             : ! Contributions from:  Sean Santos, Peter Caldwell, Xiaohong Liu and Steve Ghan
          15             : !
          16             : ! invoked in CAM by specifying -microphys=mg3
          17             : !
          18             : ! References:
          19             : !
          20             : !           Gettelman, A. and H. Morrison, Advanced Two-Moment Microphysics for Global Models.
          21             : !
          22             : !           Part I: Off line tests and comparisons with other schemes.
          23             : !
          24             : !           J. Climate, 28, 1268-1287. doi: 10.1175/JCLI-D-14-00102.1, 2015.
          25             : !
          26             : !
          27             : !
          28             : !           Gettelman, A., H. Morrison, S. Santos, P. Bogenschutz and P. H. Caldwell
          29             : !
          30             : !           Advanced Two-Moment Microphysics for Global Models.
          31             : !
          32             : !           Part II: Global model solutions and Aerosol-Cloud Interactions.
          33             : !
          34             : !           J. Climate, 28, 1288-1307. doi:10.1175/JCLI-D-14-00103.1 , 2015.
          35             : !
          36             : ! for questions contact Hugh Morrison, Andrew Gettelman
          37             : ! e-mail: morrison@ucar.edu, andrew@ucar.edu
          38             : !---------------------------------------------------------------------------------
          39             : !
          40             : ! NOTE: Modified to allow other microphysics packages (e.g. CARMA) to do ice
          41             : ! microphysics in cooperation with the MG liquid microphysics. This is
          42             : ! controlled by the do_cldice variable.
          43             : !
          44             : ! If do_cldice is false, then MG microphysics should not update CLDICE or
          45             : ! NUMICE; it is assumed that the other microphysics scheme will have updated
          46             : ! CLDICE and NUMICE. The other microphysics should handle the following
          47             : ! processes that would have been done by MG:
          48             : !   - Detrainment (liquid and ice)
          49             : !   - Homogeneous ice nucleation
          50             : !   - Heterogeneous ice nucleation
          51             : !   - Bergeron process
          52             : !   - Melting of ice
          53             : !   - Freezing of cloud drops
          54             : !   - Autoconversion (ice -> snow)
          55             : !   - Growth/Sublimation of ice
          56             : !   - Sedimentation of ice
          57             : !
          58             : ! This option has not been updated since the introduction of prognostic
          59             : ! precipitation, and probably should be adjusted to cover snow as well.
          60             : !
          61             : !---------------------------------------------------------------------------------
          62             : ! Version 3.O based on micro_mg2_0.F90 and WRF3.8.1 module_mp_morr_two_moment.F
          63             : !---------------------------------------------------------------------------------
          64             : ! Based on micro_mg (restructuring of former cldwat2m_micro)
          65             : ! Author: Andrew Gettelman, Hugh Morrison.
          66             : ! Contributions from: Xiaohong Liu and Steve Ghan
          67             : ! December 2005-May 2010
          68             : ! Description in: Morrison and Gettelman, 2008. J. Climate (MG2008)
          69             : !                 Gettelman et al., 2010 J. Geophys. Res. - Atmospheres (G2010)
          70             : ! for questions contact Hugh Morrison, Andrew Gettelman
          71             : ! e-mail: morrison@ucar.edu, andrew@ucar.edu
          72             : !---------------------------------------------------------------------------------
          73             : ! Code comments added by HM, 093011
          74             : ! General code structure:
          75             : !
          76             : ! Code is divided into two main subroutines:
          77             : !   subroutine micro_pumas_init --> initializes microphysics routine, should be called
          78             : !                                  once at start of simulation
          79             : !   subroutine micro_pumas_tend --> main microphysics routine to be called each time step
          80             : !                                this also calls several smaller subroutines to calculate
          81             : !                                microphysical processes and other utilities
          82             : !
          83             : ! List of external functions:
          84             : !   qsat_water --> for calculating saturation vapor pressure with respect to liquid water
          85             : !   qsat_ice --> for calculating saturation vapor pressure with respect to ice
          86             : !   gamma   --> standard mathematical gamma function
          87             : ! .........................................................................
          88             : ! List of inputs through use statement in fortran90:
          89             : ! Variable Name                      Description                Units
          90             : ! .........................................................................
          91             : ! gravit          acceleration due to gravity                    m s-2
          92             : ! rair            dry air gas constant for air                  J kg-1 K-1
          93             : ! tmelt           temperature of melting point for water          K
          94             : ! cpair           specific heat at constant pressure for dry air J kg-1 K-1
          95             : ! rh2o            gas constant for water vapor                  J kg-1 K-1
          96             : ! latvap          latent heat of vaporization                   J kg-1
          97             : ! latice          latent heat of fusion                         J kg-1
          98             : ! qsat_water      external function for calculating liquid water
          99             : !                 saturation vapor pressure/humidity              -
         100             : ! qsat_ice        external function for calculating ice
         101             : !                 saturation vapor pressure/humidity              pa
         102             : ! rhmini          relative humidity threshold parameter for
         103             : !                 nucleating ice                                  -
         104             : ! .........................................................................
         105             : ! NOTE: List of all inputs/outputs passed through the call/subroutine statement
         106             : !       for micro_pumas_tend is given below at the start of subroutine micro_pumas_tend.
         107             : !---------------------------------------------------------------------------------
         108             :
         109             : ! Procedures required:
         110             : ! 1) An implementation of the gamma function (if not intrinsic).
         111             : ! 2) saturation vapor pressure and specific humidity over water
         112             : ! 3) svp over ice
         113             :
         114             : #ifndef HAVE_GAMMA_INTRINSICS
         115             : use shr_spfn_mod, only: gamma => shr_spfn_gamma
         116             : #endif
         117             :
         118             : use wv_sat_methods, only: &
         119             :      qsat_water => wv_sat_qsat_water_vect, &
         120             :      qsat_ice => wv_sat_qsat_ice_vect
         121             :
         122             : ! Parameters from the utilities module.
         123             : use micro_pumas_utils, only: &
         124             :      r8, &
         125             :      pi, &
         126             :      omsm, &
         127             :      qsmall, &
         128             :      mincld, &
         129             :      rhosn, &
         130             :      rhoi, &
         131             :      rhow, &
         132             :      rhows, &
         133             :      ac, bc, &
         134             :      ai, bi, &
         135             :      aj, bj, &
         136             :      ar, br, &
         137             :      as, bs, &
         138             :      ag, bg, &
         139             :      ah, bh, &
         140             :      rhog,rhoh, &
         141             :      mi0, &
         142             :      rising_factorial
         143             :
         144             : implicit none
         145             : private
         146             : save
         147             :
         148             : public :: &
         149             :      micro_pumas_init, &
         150             :      micro_pumas_get_cols, &
         151             :      micro_pumas_tend
         152             :
         153             : ! Switches for specification rather than prediction of droplet and crystal number
         154             : ! note: number will be adjusted as needed to keep mean size within bounds,
         155             : ! even when specified droplet or ice number is used
         156             : !
         157             : ! If constant cloud ice number is set (nicons = .true.),
         158             : ! then all microphysical processes except mass transfer due to ice nucleation
         159             : ! (mnuccd) are based on the fixed cloud ice number. Calculation of
         160             : ! mnuccd follows from the prognosed ice crystal number ni.
         161             :
         162             : logical :: nccons ! nccons = .true. to specify constant cloud droplet number
         163             : logical :: nicons ! nicons = .true. to specify constant cloud ice number
         164             : logical :: ngcons ! ngcons = .true. to specify constant graupel number
         165             : logical :: nrcons ! constant rain number
         166             : logical :: nscons ! constant snow number
         167             :
         168             : ! specified ice and droplet number concentrations
         169             : ! note: these are local in-cloud values, not grid-mean
         170             : real(r8) :: ncnst  ! droplet num concentration when nccons=.true. (m-3)
         171             : real(r8) :: ninst  ! ice num concentration when nicons=.true. (m-3)
         172             : real(r8) :: ngnst   ! graupel num concentration when ngcons=.true. (m-3)
         173             : real(r8) :: nrnst
         174             : real(r8) :: nsnst
         175             :
         176             : ! IFS Switches....
         177             : ! Switch to turn off evaporation of sedimenting condensate
         178             : ! Found to interact badly in some models with diagnostic cloud fraction
         179             : logical :: evap_sed_off
         180             :
         181             : ! Remove RH conditional from ice nucleation
         182             : logical :: icenuc_rh_off
         183             :
         184             : ! Internally: Meyers Ice Nucleation
         185             : logical :: icenuc_use_meyers
         186             :
         187             : ! Scale evaporation as IFS does (*0.3)
         188             : logical :: evap_scl_ifs
         189             :
         190             : ! Evap RH threhold following ifs
         191             : logical :: evap_rhthrsh_ifs
         192             :
         193             : ! Rain freezing at 0C following ifs
         194             :
         195             : logical :: rainfreeze_ifs
         196             :
         197             : ! Snow sedimentation = 1 m/s
         198             :
         199             : logical :: ifs_sed
         200             :
         201             : ! Precipitation fall speed, prevent zero velocity if precip above
         202             :
         203             : logical :: precip_fall_corr
         204             :
         205             : !--ag
         206             :
         207             : !=========================================================
         208             : ! Private module parameters
         209             : !=========================================================
         210             :
         211             : !Range of cloudsat reflectivities (dBz) for analytic simulator
         212             : real(r8), parameter :: csmin = -30._r8
         213             : real(r8), parameter :: csmax = 26._r8
         214             : real(r8), parameter :: mindbz = -99._r8
         215             : real(r8), parameter :: minrefl = 1.26e-10_r8    ! minrefl = 10._r8**(mindbz/10._r8)
         216             :
         217             : integer, parameter  :: MG_PRECIP_FRAC_INCLOUD = 101
         218             : integer, parameter  :: MG_PRECIP_FRAC_OVERLAP = 102
         219             :
         220             : ! autoconversion size threshold for cloud ice to snow (m)
         221             : real(r8) :: dcs
         222             :
         223             : ! minimum mass of new crystal due to freezing of cloud droplets done
         224             : ! externally (kg)
         225             : real(r8), parameter :: mi0l_min = 4._r8/3._r8*pi*rhow*(4.e-6_r8)**3
         226             :
         227             : ! Ice number sublimation parameter. Assume some decrease in ice number with sublimation if non-zero. Else, no decrease in number with sublimation.
         228             : real(r8), parameter :: sublim_factor =0.0_r8      !number sublimation factor.
         229             :
         230             : integer, parameter :: VLENS = 128  ! vector length of a GPU compute kernel
         231             :
         232             : !=========================================================
         233             : ! Constants set in initialization
         234             : !=========================================================
         235             :
         236             : ! Set using arguments to micro_pumas_init
         237             : real(r8) :: g           ! gravity
         238             : real(r8) :: r           ! dry air gas constant
         239             : real(r8) :: rv          ! water vapor gas constant
         240             : real(r8) :: cpp         ! specific heat of dry air
         241             : real(r8) :: tmelt       ! freezing point of water (K)
         242             :
         243             : ! latent heats of:
         244             : real(r8) :: xxlv        ! vaporization
         245             : real(r8) :: xlf         ! freezing
         246             : real(r8) :: xxls        ! sublimation
         247             :
         248             : real(r8) :: rhmini      ! Minimum rh for ice cloud fraction > 0.
         249             :
         250             : ! flags
         251             : logical :: microp_uniform
         252             : logical :: do_cldice
         253             : logical :: use_hetfrz_classnuc
         254             : logical :: do_hail
         255             : logical :: do_graupel
         256             :
         257             : real(r8) :: rhosu       ! typical 850mn air density
         258             :
         259             : real(r8) :: icenuct     ! ice nucleation temperature: currently -5 degrees C
         260             :
         261             : real(r8) :: snowmelt    ! what temp to melt all snow: currently 2 degrees C
         262             : real(r8) :: rainfrze    ! what temp to freeze all rain: currently -5 degrees C
         263             :
         264             : ! additional constants to help speed up code
         265             : real(r8) :: gamma_br_plus1
         266             : real(r8) :: gamma_br_plus4
         267             : real(r8) :: gamma_bs_plus1
         268             : real(r8) :: gamma_bs_plus4
         269             : real(r8) :: gamma_bi_plus1
         270             : real(r8) :: gamma_bi_plus4
         271             : real(r8) :: gamma_bj_plus1
         272             : real(r8) :: gamma_bj_plus4
         273             : real(r8) :: gamma_bg_plus1
         274             : real(r8) :: gamma_bg_plus4
         275             : real(r8) :: xxlv_squared
         276             : real(r8) :: xxls_squared
         277             :
         278             : character(len=16)  :: micro_mg_precip_frac_method  ! type of precipitation fraction method
         279             : real(r8)           :: micro_mg_berg_eff_factor     ! berg efficiency factor
         280             :
         281             : real(r8)           :: micro_mg_accre_enhan_fact     ! accretion enhancment factor
         282             : real(r8)           :: micro_mg_autocon_fact     ! autoconversion prefactor
         283             : real(r8)           :: micro_mg_autocon_nd_exp     ! autoconversion Nd exponent factor
         284             : real(r8)           :: micro_mg_autocon_lwp_exp  !autoconversion LWP exponent
         285             : real(r8)           :: micro_mg_homog_size ! size of freezing homogeneous ice
         286             : real(r8)           :: micro_mg_vtrmi_factor
         287             : real(r8)           :: micro_mg_effi_factor
         288             : real(r8)           :: micro_mg_iaccr_factor
         289             : real(r8)           :: micro_mg_max_nicons
         290             :
         291             : logical  :: remove_supersat ! If true, remove supersaturation after sedimentation loop
         292             : logical  :: do_sb_physics ! do SB 2001 autoconversion or accretion physics
         293             :
         294             : !$acc declare create (nccons,nicons,ngcons,nrcons,nscons,ncnst,ninst,ngnst,   &
         295             : !$acc                 nrnst,nsnst,evap_sed_off,icenuc_rh_off,evap_scl_ifs,    &
         296             : !$acc                 icenuc_use_meyers,evap_rhthrsh_ifs,rainfreeze_ifs,      &
         297             : !$acc                 ifs_sed,precip_fall_corr,dcs,                           &
         298             : !$acc                 g,r,rv,cpp,tmelt,xxlv,xlf,xxls,rhmini,microp_uniform,   &
         299             : !$acc                 do_cldice,use_hetfrz_classnuc,do_hail,do_graupel,rhosu, &
         300             : !$acc                 icenuct,snowmelt,rainfrze,xxlv_squared,xxls_squared,    &
         301             : !$acc                 gamma_br_plus1,gamma_br_plus4,gamma_bs_plus1,           &
         302             : !$acc                 gamma_bs_plus4,gamma_bi_plus1,gamma_bi_plus4,           &
         303             : !$acc                 gamma_bj_plus1,gamma_bj_plus4,gamma_bg_plus1,           &
         304             : !$acc                 gamma_bg_plus4,micro_mg_berg_eff_factor,                &
         305             : !$acc                 remove_supersat,do_sb_physics)
         306             :
         307             : !===============================================================================
         308             : contains
         309             : !===============================================================================
         310             :
         311           0 : subroutine micro_pumas_init( &
         312             :      kind, gravit, rair, rh2o, cpair,    &
         313             :      tmelt_in, latvap, latice,           &
         314             :      rhmini_in, micro_mg_dcs,            &
         315             :      micro_mg_do_hail_in,micro_mg_do_graupel_in, &
         316             :      microp_uniform_in, do_cldice_in, use_hetfrz_classnuc_in, &
         317           0 :      micro_mg_precip_frac_method_in, micro_mg_berg_eff_factor_in, &
         318             :      micro_mg_accre_enhan_fact_in, micro_mg_autocon_fact_in, &
         319             :      micro_mg_autocon_nd_exp_in, micro_mg_autocon_lwp_exp_in, micro_mg_homog_size_in, &
         320             :      micro_mg_vtrmi_factor_in, micro_mg_effi_factor_in,  micro_mg_iaccr_factor_in,&
         321             :      micro_mg_max_nicons_in, &
         322             :      remove_supersat_in, do_sb_physics_in, &
         323             :      micro_mg_evap_sed_off_in, micro_mg_icenuc_rh_off_in, micro_mg_icenuc_use_meyers_in, &
         324             :      micro_mg_evap_scl_ifs_in, micro_mg_evap_rhthrsh_ifs_in, &
         325             :      micro_mg_rainfreeze_ifs_in,  micro_mg_ifs_sed_in, micro_mg_precip_fall_corr, &
         326             :      nccons_in, nicons_in, ncnst_in, ninst_in, ngcons_in, ngnst_in, &
         327             :      nrcons_in, nrnst_in, nscons_in, nsnst_in, &
         328           0 :      errstring)
         329             :
         330             :   use micro_pumas_utils, only: micro_pumas_utils_init
         331             :
         332             :   !-----------------------------------------------------------------------
         333             :   !
         334             :   ! Purpose:
         335             :   ! initialize constants for MG microphysics
         336             :   !
         337             :   ! Author: Andrew Gettelman Dec 2005
         338             :   !
         339             :   !-----------------------------------------------------------------------
         340             :
         341             :   integer,  intent(in)  :: kind         ! Kind used for reals
         342             :   real(r8), intent(in)  :: gravit
         343             :   real(r8), intent(in)  :: rair
         344             :   real(r8), intent(in)  :: rh2o
         345             :   real(r8), intent(in)  :: cpair
         346             :   real(r8), intent(in)  :: tmelt_in     ! Freezing point of water (K)
         347             :   real(r8), intent(in)  :: latvap
         348             :   real(r8), intent(in)  :: latice
         349             :   real(r8), intent(in)  :: rhmini_in    ! Minimum rh for ice cloud fraction > 0.
         350             :   real(r8), intent(in)  :: micro_mg_dcs
         351             :
         352             : !MG3 dense precipitating ice. Note, only 1 can be true, or both false.
         353             :   logical,  intent(in)  :: micro_mg_do_graupel_in    ! .true. = configure with graupel
         354             :                                                    ! .false. = no graupel (hail possible)
         355             :   logical,  intent(in)  :: micro_mg_do_hail_in    ! .true. = configure with hail
         356             :                                                    ! .false. = no hail (graupel possible)
         357             :   logical,  intent(in)  :: microp_uniform_in    ! .true. = configure uniform for sub-columns
         358             :                                             ! .false. = use w/o sub-columns (standard)
         359             :   logical,  intent(in)  :: do_cldice_in     ! .true. = do all processes (standard)
         360             :                                             ! .false. = skip all processes affecting
         361             :                                             !           cloud ice
         362             :   logical,  intent(in)  :: use_hetfrz_classnuc_in ! use heterogeneous freezing
         363             :
         364             :   character(len=16),intent(in)  :: micro_mg_precip_frac_method_in  ! type of precipitation fraction method
         365             :   real(r8),         intent(in)  :: micro_mg_berg_eff_factor_in     ! berg efficiency factor
         366             :   real(r8),         intent(in)  :: micro_mg_accre_enhan_fact_in     !accretion enhancment factor
         367             :   real(r8),         intent(in) ::  micro_mg_autocon_fact_in    !autconversion prefactor
         368             :   real(r8),         intent(in) ::  micro_mg_autocon_nd_exp_in !autconversion exponent factor
         369             :   real(r8),         intent(in) ::  micro_mg_autocon_lwp_exp_in    !autconversion exponent factor
         370             :   real(r8),         intent(in) ::  micro_mg_homog_size_in  ! size of homoegenous freezing ice
         371             :   real(r8),         intent(in)  :: micro_mg_vtrmi_factor_in    !factor for ice fall velocity
         372             :   real(r8),         intent(in)  :: micro_mg_effi_factor_in    !factor for ice effective radius
         373             :   real(r8),         intent(in)  :: micro_mg_iaccr_factor_in  ! ice accretion factor
         374             :   real(r8),         intent(in)  :: micro_mg_max_nicons_in ! maximum number ice crystal allowed
         375             :
         376             :   logical,  intent(in)  ::  remove_supersat_in ! If true, remove supersaturation after sedimentation loop
         377             :   logical,  intent(in)  ::  do_sb_physics_in ! do SB autoconversion and accretion physics
         378             :
         379             : ! IFS-like Switches
         380             :
         381             :   logical, intent(in) :: micro_mg_evap_sed_off_in ! Turn off evaporation/sublimation based on cloud fraction for sedimenting condensate
         382             :
         383             :   logical, intent(in) :: micro_mg_icenuc_rh_off_in ! Remove RH conditional from ice nucleation
         384             :   logical, intent(in) :: micro_mg_icenuc_use_meyers_in ! Internally: Meyers Ice Nucleation
         385             :   logical, intent(in) :: micro_mg_evap_scl_ifs_in ! Scale evaporation as IFS does (*0.3)
         386             :   logical, intent(in) :: micro_mg_evap_rhthrsh_ifs_in ! Evap RH threhold following ifs
         387             :   logical, intent(in) :: micro_mg_rainfreeze_ifs_in ! Rain freezing temp following ifs
         388             :   logical, intent(in) :: micro_mg_ifs_sed_in ! snow sedimentation = 1m/s following ifs
         389             :   logical, intent(in) :: micro_mg_precip_fall_corr ! ensure rain fall speed non-zero if rain above in column
         390             :
         391             :   logical, intent(in)   :: nccons_in
         392             :   logical, intent(in)   :: nicons_in
         393             :   real(r8), intent(in)  :: ncnst_in
         394             :   real(r8), intent(in)  :: ninst_in
         395             :
         396             :   logical, intent(in)   :: ngcons_in
         397             :   real(r8), intent(in)  :: ngnst_in
         398             :   logical, intent(in)   :: nrcons_in
         399             :   real(r8), intent(in)  :: nrnst_in
         400             :   logical, intent(in)   :: nscons_in
         401             :   real(r8), intent(in)  :: nsnst_in
         402             :
         403             :   character(128), intent(out) :: errstring    ! Output status (non-blank for error return)
         404             :
         405             :   !-----------------------------------------------------------------------
         406             :
         407           0 :   dcs = micro_mg_dcs
         408             :
         409             :   ! Initialize subordinate utilities module.
         410             :   call micro_pumas_utils_init(kind, rair, rh2o, cpair, tmelt_in, latvap, latice, &
         411           0 :        dcs, errstring)
         412             :
         413           0 :   if (trim(errstring) /= "") return
         414             :
         415             :   ! declarations for MG code (transforms variable names)
         416             :
         417           0 :   g= gravit                 ! gravity
         418           0 :   r= rair                   ! dry air gas constant: note units(phys_constants are in J/K/kmol)
         419           0 :   rv= rh2o                  ! water vapor gas constant
         420           0 :   cpp = cpair               ! specific heat of dry air
         421           0 :   tmelt = tmelt_in
         422           0 :   rhmini = rhmini_in
         423           0 :   micro_mg_precip_frac_method = micro_mg_precip_frac_method_in
         424           0 :   micro_mg_berg_eff_factor    = micro_mg_berg_eff_factor_in
         425           0 :   micro_mg_accre_enhan_fact   =  micro_mg_accre_enhan_fact_in
         426           0 :   micro_mg_autocon_fact  = micro_mg_autocon_fact_in
         427           0 :   micro_mg_autocon_nd_exp = micro_mg_autocon_nd_exp_in
         428           0 :   micro_mg_autocon_lwp_exp = micro_mg_autocon_lwp_exp_in
         429           0 :   micro_mg_homog_size   = micro_mg_homog_size_in
         430           0 :   micro_mg_vtrmi_factor = micro_mg_vtrmi_factor_in
         431           0 :   micro_mg_effi_factor = micro_mg_effi_factor_in
         432           0 :   micro_mg_iaccr_factor = micro_mg_iaccr_factor_in
         433           0 :   micro_mg_max_nicons = micro_mg_max_nicons_in
         434           0 :   remove_supersat          = remove_supersat_in
         435           0 :   do_sb_physics               = do_sb_physics_in
         436             :
         437           0 :   nccons = nccons_in
         438           0 :   nicons = nicons_in
         439           0 :   ncnst  = ncnst_in
         440           0 :   ninst  = ninst_in
         441           0 :   ngcons = ngcons_in
         442           0 :   ngnst  = ngnst_in
         443           0 :   nscons = nscons_in
         444           0 :   nsnst  = nsnst_in
         445           0 :   nrcons = nrcons_in
         446           0 :   nrnst  = nrnst_in
         447             :
         448             :   ! latent heats
         449             :
         450           0 :   xxlv = latvap         ! latent heat vaporization
         451           0 :   xlf  = latice         ! latent heat freezing
         452           0 :   xxls = xxlv + xlf     ! latent heat of sublimation
         453             :
         454             :   ! flags
         455           0 :   microp_uniform = microp_uniform_in
         456           0 :   do_cldice  = do_cldice_in
         457           0 :   use_hetfrz_classnuc = use_hetfrz_classnuc_in
         458           0 :   do_hail = micro_mg_do_hail_in
         459           0 :   do_graupel = micro_mg_do_graupel_in
         460           0 :   evap_sed_off = micro_mg_evap_sed_off_in
         461           0 :   icenuc_rh_off = micro_mg_icenuc_rh_off_in
         462           0 :   icenuc_use_meyers = micro_mg_icenuc_use_meyers_in
         463           0 :   evap_scl_ifs = micro_mg_evap_scl_ifs_in
         464           0 :   evap_rhthrsh_ifs = micro_mg_evap_rhthrsh_ifs_in
         465           0 :   rainfreeze_ifs = micro_mg_rainfreeze_ifs_in
         466           0 :   ifs_sed = micro_mg_ifs_sed_in
         467           0 :   precip_fall_corr = micro_mg_precip_fall_corr
         468             :   ! typical air density at 850 mb
         469             :
         470           0 :   rhosu = 85000._r8/(rair * tmelt)
         471             :
         472             :   ! Maximum temperature at which snow is allowed to exist
         473           0 :   snowmelt = tmelt + 2._r8
         474             :   ! Minimum temperature at which rain is allowed to exist
         475           0 :    if (rainfreeze_ifs) then
         476           0 :       rainfrze = tmelt
         477             :    else
         478           0 :       rainfrze = tmelt - 40._r8
         479             :    end if
         480             :
         481             :
         482             :   ! Ice nucleation temperature
         483           0 :   icenuct  = tmelt - 5._r8
         484             :
         485             :   ! Define constants to help speed up code (this limits calls to gamma function)
         486           0 :   gamma_br_plus1=gamma(1._r8+br)
         487           0 :   gamma_br_plus4=gamma(4._r8+br)
         488           0 :   gamma_bs_plus1=gamma(1._r8+bs)
         489           0 :   gamma_bs_plus4=gamma(4._r8+bs)
         490           0 :   gamma_bi_plus1=gamma(1._r8+bi)
         491           0 :   gamma_bi_plus4=gamma(4._r8+bi)
         492           0 :   gamma_bj_plus1=gamma(1._r8+bj)
         493           0 :   gamma_bj_plus4=gamma(4._r8+bj)
         494           0 :   gamma_bg_plus1=gamma(1._r8)
         495           0 :   gamma_bg_plus4=gamma(4._r8)
         496           0 :   if (do_hail) then
         497           0 :      gamma_bg_plus1 = gamma(1._r8+bh)
         498           0 :      gamma_bg_plus4 = gamma(4._r8+bh)
         499             :   end if
         500           0 :   if (do_graupel) then
         501           0 :      gamma_bg_plus1 = gamma(1._r8+bg)
         502           0 :      gamma_bg_plus4 = gamma(4._r8+bg)
         503             :   end if
         504             :
         505           0 :   xxlv_squared=xxlv**2
         506           0 :   xxls_squared=xxls**2
         507             :
         508             :   !$acc update device (nccons,nicons,ngcons,nrcons,nscons,ncnst,ninst,ngnst,   &
         509             :   !$acc                nrnst,nsnst,evap_sed_off,icenuc_rh_off,evap_scl_ifs,    &
         510             :   !$acc                icenuc_use_meyers,evap_rhthrsh_ifs,rainfreeze_ifs,      &
         511             :   !$acc                ifs_sed,precip_fall_corr,dcs,                           &
         512             :   !$acc                g,r,rv,cpp,tmelt,xxlv,xlf,xxls,rhmini,microp_uniform,   &
         513             :   !$acc                do_cldice,use_hetfrz_classnuc,do_hail,do_graupel,rhosu, &
         514             :   !$acc                icenuct,snowmelt,rainfrze,xxlv_squared,xxls_squared,    &
         515             :   !$acc                gamma_br_plus1,gamma_br_plus4,gamma_bs_plus1,           &
         516             :   !$acc                gamma_bs_plus4,gamma_bi_plus1,gamma_bi_plus4,           &
         517             :   !$acc                gamma_bj_plus1,gamma_bj_plus4,gamma_bg_plus1,           &
         518             :   !$acc                gamma_bg_plus4,micro_mg_berg_eff_factor,                &
         519             :   !$acc                remove_supersat,do_sb_physics)
         520             :
         521             : end subroutine micro_pumas_init
         522             :
         523             : !===============================================================================
         524             : !microphysics routine for each timestep goes here...
         525             :
         526           0 : subroutine micro_pumas_tend ( &
         527             :      mgncol,             nlev,               deltatin,           &
         528           0 :      t,                            q,                            &
         529           0 :      qcn,                          qin,                          &
         530           0 :      ncn,                          nin,                          &
         531           0 :      qrn,                          qsn,                          &
         532           0 :      nrn,                          nsn,                          &
         533           0 :      qgr,                          ngr,                          &
         534           0 :      relvar,                       accre_enhan,                  &
         535           0 :      p,                            pdel,                         &
         536           0 :      cldn,    liqcldf,        icecldf,       qsatfac,            &
         537           0 :      qcsinksum_rate1ord,                                         &
         538           0 :      naai,                         npccn,                        &
         539           0 :      rndst,                        nacon,                        &
         540           0 :      tlat,                         qvlat,                        &
         541           0 :      qctend,                       qitend,                       &
         542           0 :      nctend,                       nitend,                       &
         543           0 :      qrtend,                       qstend,                       &
         544           0 :      nrtend,                       nstend,                       &
         545           0 :      qgtend,                       ngtend,                       &
         546           0 :      effc,               effc_fn,            effi,               &
         547           0 :      sadice,                       sadsnow,                      &
         548           0 :      prect,                        preci,                        &
         549           0 :      nevapr,                       evapsnow,                     &
         550           0 :      am_evp_st,                                                  &
         551           0 :      prain,                        prodsnow,                     &
         552           0 :      cmeout,                       deffi,                        &
         553           0 :      pgamrad,                      lamcrad,                      &
         554           0 :      qsout,                        dsout,                        &
         555           0 :      qgout,     ngout,             dgout,                        &
         556           0 :      lflx,               iflx,                                   &
         557           0 :      gflx,                                                       &
         558           0 :      rflx,               sflx,               qrout,              &
         559           0 :      reff_rain,          reff_snow,          reff_grau,          &
         560           0 :      qcsevap,            qisevap,            qvres,              &
         561           0 :      cmeitot,            vtrmc,              vtrmi,              &
         562           0 :      umr,                          ums,                          &
         563           0 :      umg,                qgsedten,                               &
         564           0 :      qcsedten,                     qisedten,                     &
         565           0 :      qrsedten,                     qssedten,                     &
         566           0 :      pratot,                       prctot,                       &
         567           0 :      mnuccctot,          mnuccttot,          msacwitot,          &
         568           0 :      psacwstot,          bergstot,           bergtot,            &
         569           0 :      melttot,    meltstot,   meltgtot,       homotot,            &
         570           0 :      qcrestot,           prcitot,            praitot,            &
         571           0 :      qirestot,           mnuccrtot,    mnudeptot,       mnuccritot, pracstot,           &
         572           0 :      meltsdttot,         frzrdttot,          mnuccdtot,          &
         573           0 :      pracgtot,           psacwgtot,          pgsacwtot,          &
         574           0 :      pgracstot,          prdgtot,           &
         575           0 :      qmultgtot,          qmultrgtot,         psacrtot,           &
         576           0 :      npracgtot,          nscngtot,           ngracstot,          &
         577           0 :      nmultgtot,          nmultrgtot,         npsacwgtot,         &
         578           0 :      nrout,                        nsout,                        &
         579           0 :      refl,               arefl,              areflz,             &
         580           0 :      frefl,              csrfl,              acsrfl,             &
         581           0 :      fcsrfl,                       rercld,                       &
         582           0 :      ncai,                         ncal,                         &
         583           0 :      qrout2,                       qsout2,                       &
         584           0 :      nrout2,                       nsout2,                       &
         585           0 :      drout2,                       dsout2,                       &
         586           0 :      qgout2,        ngout2,        dgout2,    freqg,                   &
         587           0 :      freqs,                        freqr,                        &
         588           0 :      nfice,                        qcrat,                        &
         589           0 :      errstring, & ! Below arguments are "optional" (pass null pointers to omit).
         590           0 :      tnd_qsnow,          tnd_nsnow,          re_ice,             &
         591           0 :      prer_evap,                                                      &
         592           0 :      frzimm,             frzcnt,             frzdep)
         593             :
         594             :   ! Constituent properties.
         595             :   use micro_pumas_utils, only: &
         596             :        mg_liq_props, &
         597             :        mg_ice_props, &
         598             :        mg_rain_props, &
         599             :        mg_graupel_props, &
         600             :        mg_hail_props, &
         601             :        mg_snow_props
         602             :
         603             :   ! Size calculation functions.
         604             :   use micro_pumas_utils, only: &
         605             :        size_dist_param_liq, &
         606             :        size_dist_param_basic, &
         607             :        avg_diameter, &
         608             :        avg_diameter_vec
         609             :
         610             :   ! Microphysical processes.
         611             :   use micro_pumas_utils, only: &
         612             :        ice_deposition_sublimation, &
         613             :        sb2001v2_liq_autoconversion,&
         614             :        sb2001v2_accre_cld_water_rain,&
         615             :        kk2000_liq_autoconversion, &
         616             :        ice_autoconversion, &
         617             :        immersion_freezing, &
         618             :        contact_freezing, &
         619             :        snow_self_aggregation, &
         620             :        accrete_cloud_water_snow, &
         621             :        secondary_ice_production, &
         622             :        accrete_rain_snow, &
         623             :        heterogeneous_rain_freezing, &
         624             :        accrete_cloud_water_rain, &
         625             :        self_collection_rain, &
         626             :        accrete_cloud_ice_snow, &
         627             :        evaporate_sublimate_precip, &
         628             :        bergeron_process_snow, &
         629             :        graupel_collecting_snow, &
         630             :        graupel_collecting_rain, &
         631             :        graupel_collecting_cld_water, &
         632             :        graupel_riming_liquid_snow, &
         633             :        graupel_rain_riming_snow, &
         634             :        graupel_rime_splintering, &
         635             :        evaporate_sublimate_precip_graupel
         636             :
         637             :   !Authors: Hugh Morrison, Andrew Gettelman, NCAR, Peter Caldwell, LLNL
         638             :   ! e-mail: morrison@ucar.edu, andrew@ucar.edu
         639             :
         640             :   ! input arguments
         641             :   integer,  intent(in) :: mgncol         ! number of microphysics columns
         642             :   integer,  intent(in) :: nlev           ! number of layers
         643             :   real(r8), intent(in) :: deltatin       ! time step (s)
         644             :   real(r8), intent(in) :: t(mgncol,nlev) ! input temperature (K)
         645             :   real(r8), intent(in) :: q(mgncol,nlev) ! input h20 vapor mixing ratio (kg/kg)
         646             :
         647             :   ! note: all input cloud variables are grid-averaged
         648             :   real(r8), intent(in) :: qcn(mgncol,nlev)       ! cloud water mixing ratio (kg/kg)
         649             :   real(r8), intent(in) :: qin(mgncol,nlev)       ! cloud ice mixing ratio (kg/kg)
         650             :   real(r8), intent(in) :: ncn(mgncol,nlev)       ! cloud water number conc (1/kg)
         651             :   real(r8), intent(in) :: nin(mgncol,nlev)       ! cloud ice number conc (1/kg)
         652             :
         653             :   real(r8), intent(in) :: qrn(mgncol,nlev)       ! rain mixing ratio (kg/kg)
         654             :   real(r8), intent(in) :: qsn(mgncol,nlev)       ! snow mixing ratio (kg/kg)
         655             :   real(r8), intent(in) :: nrn(mgncol,nlev)       ! rain number conc (1/kg)
         656             :   real(r8), intent(in) :: nsn(mgncol,nlev)       ! snow number conc (1/kg)
         657             :   real(r8), intent(in) :: qgr(mgncol,nlev)       ! graupel/hail mixing ratio (kg/kg)
         658             :   real(r8), intent(in) :: ngr(mgncol,nlev)       ! graupel/hail number conc (1/kg)
         659             :
         660             :   real(r8), intent(in) :: relvar(mgncol,nlev)      ! cloud water relative variance (-)
         661             :   real(r8), intent(in) :: accre_enhan(mgncol,nlev) ! optional accretion
         662             :                                              ! enhancement factor (-)
         663             :
         664             :   real(r8), intent(in) :: p(mgncol,nlev)        ! air pressure (pa)
         665             :   real(r8), intent(in) :: pdel(mgncol,nlev)     ! pressure difference across level (pa)
         666             :
         667             :   real(r8), intent(in) :: cldn(mgncol,nlev)      ! cloud fraction (no units)
         668             :   real(r8), intent(in) :: liqcldf(mgncol,nlev)   ! liquid cloud fraction (no units)
         669             :   real(r8), intent(in) :: icecldf(mgncol,nlev)   ! ice cloud fraction (no units)
         670             :   real(r8), intent(in) :: qsatfac(mgncol,nlev)   ! subgrid cloud water saturation scaling factor (no units)
         671             :
         672             :   ! used for scavenging
         673             :   ! Inputs for aerosol activation
         674             :   real(r8), intent(in) :: naai(mgncol,nlev)     ! ice nucleation number (from microp_aero_ts) (1/kg)
         675             :   real(r8), intent(in) :: npccn(mgncol,nlev)   ! ccn activated number tendency (from microp_aero_ts) (1/kg*s)
         676             :
         677             :   ! Note that for these variables, the dust bin is assumed to be the last index.
         678             :   ! (For example, in CAM, the last dimension is always size 4.)
         679             :   real(r8), intent(in) :: rndst(:,:,:)  ! radius of each dust bin, for contact freezing (from microp_aero_ts) (m)
         680             :   real(r8), intent(in) :: nacon(:,:,:) ! number in each dust bin, for contact freezing  (from microp_aero_ts) (1/m^3)
         681             :
         682             :   ! output arguments
         683             :
         684             :   real(r8), intent(out) :: qcsinksum_rate1ord(mgncol,nlev) ! 1st order rate for
         685             :   ! direct cw to precip conversion
         686             :   real(r8), intent(out) :: tlat(mgncol,nlev)         ! latent heating rate       (W/kg)
         687             :   real(r8), intent(out) :: qvlat(mgncol,nlev)        ! microphysical tendency qv (1/s)
         688             :   real(r8), intent(out) :: qctend(mgncol,nlev)       ! microphysical tendency qc (1/s)
         689             :   real(r8), intent(out) :: qitend(mgncol,nlev)       ! microphysical tendency qi (1/s)
         690             :   real(r8), intent(out) :: nctend(mgncol,nlev)       ! microphysical tendency nc (1/(kg*s))
         691             :   real(r8), intent(out) :: nitend(mgncol,nlev)       ! microphysical tendency ni (1/(kg*s))
         692             :
         693             :   real(r8), intent(out) :: qrtend(mgncol,nlev)       ! microphysical tendency qr (1/s)
         694             :   real(r8), intent(out) :: qstend(mgncol,nlev)       ! microphysical tendency qs (1/s)
         695             :   real(r8), intent(out) :: nrtend(mgncol,nlev)       ! microphysical tendency nr (1/(kg*s))
         696             :   real(r8), intent(out) :: nstend(mgncol,nlev)       ! microphysical tendency ns (1/(kg*s))
         697             :   real(r8), intent(out) :: qgtend(mgncol,nlev)       ! microphysical tendency qg (1/s)
         698             :   real(r8), intent(out) :: ngtend(mgncol,nlev)       ! microphysical tendency ng (1/(kg*s))
         699             :
         700             :   real(r8), intent(out) :: effc(mgncol,nlev)         ! droplet effective radius (micron)
         701             :   real(r8), intent(out) :: effc_fn(mgncol,nlev)      ! droplet effective radius, assuming nc = 1.e8 kg-1
         702             :   real(r8), intent(out) :: effi(mgncol,nlev)         ! cloud ice effective radius (micron)
         703             :   real(r8), intent(out) :: sadice(mgncol,nlev)       ! cloud ice surface area density (cm2/cm3)
         704             :   real(r8), intent(out) :: sadsnow(mgncol,nlev)      ! cloud snow surface area density (cm2/cm3)
         705             :   real(r8), intent(out) :: prect(mgncol)             ! surface precip rate (m/s)
         706             :   real(r8), intent(out) :: preci(mgncol)             ! cloud ice/snow precip rate (m/s)
         707             :   real(r8), intent(out) :: nevapr(mgncol,nlev)       ! evaporation rate of rain + snow (1/s)
         708             :   real(r8), intent(out) :: evapsnow(mgncol,nlev)     ! sublimation rate of snow (1/s)
         709             :   real(r8), intent(out) :: am_evp_st(mgncol,nlev)    ! stratiform evaporation area (frac)
         710             :   real(r8), intent(out) :: prain(mgncol,nlev)        ! production of rain + snow (1/s)
         711             :   real(r8), intent(out) :: prodsnow(mgncol,nlev)     ! production of snow (1/s)
         712             :   real(r8), intent(out) :: cmeout(mgncol,nlev)       ! evap/sub of cloud (1/s)
         713             :   real(r8), intent(out) :: deffi(mgncol,nlev)        ! ice effective diameter for optics (radiation) (micron)
         714             :   real(r8), intent(out) :: pgamrad(mgncol,nlev)      ! ice gamma parameter for optics (radiation) (no units)
         715             :   real(r8), intent(out) :: lamcrad(mgncol,nlev)      ! slope of droplet distribution for optics (radiation) (1/m)
         716             :   real(r8), intent(out) :: qsout(mgncol,nlev)        ! snow mixing ratio (kg/kg)
         717             :   real(r8), intent(out) :: dsout(mgncol,nlev)        ! snow diameter (m)
         718             :   real(r8), intent(out) :: lflx(mgncol,nlev+1)       ! grid-box average liquid condensate flux (kg m^-2 s^-1)
         719             :   real(r8), intent(out) :: iflx(mgncol,nlev+1)       ! grid-box average ice condensate flux (kg m^-2 s^-1)
         720             :   real(r8), intent(out) :: rflx(mgncol,nlev+1)       ! grid-box average rain flux (kg m^-2 s^-1)
         721             :   real(r8), intent(out) :: sflx(mgncol,nlev+1)       ! grid-box average snow flux (kg m^-2 s^-1)
         722             :   real(r8), intent(out) :: gflx(mgncol,nlev+1)       ! grid-box average graupel/hail flux (kg m^-2 s^-1)
         723             :
         724             :   real(r8), intent(out) :: qrout(mgncol,nlev)        ! grid-box average rain mixing ratio (kg/kg)
         725             :   real(r8), intent(out) :: reff_rain(mgncol,nlev)    ! rain effective radius (micron)
         726             :   real(r8), intent(out) :: reff_snow(mgncol,nlev)    ! snow effective radius (micron)
         727             :   real(r8), intent(out) :: reff_grau(mgncol,nlev)    ! graupel effective radius (micron)
         728             :   real(r8), intent(out) :: qcsevap(mgncol,nlev)      ! cloud water evaporation due to sedimentation (1/s)
         729             :   real(r8), intent(out) :: qisevap(mgncol,nlev)      ! cloud ice sublimation due to sublimation (1/s)
         730             :   real(r8), intent(out) :: qvres(mgncol,nlev)        ! residual condensation term to ensure RH < 100% (1/s)
         731             :   real(r8), intent(out) :: cmeitot(mgncol,nlev)      ! grid-mean cloud ice sub/dep (1/s)
         732             :   real(r8), intent(out) :: vtrmc(mgncol,nlev)        ! mass-weighted cloud water fallspeed (m/s)
         733             :   real(r8), intent(out) :: vtrmi(mgncol,nlev)        ! mass-weighted cloud ice fallspeed (m/s)
         734             :   real(r8), intent(out) :: umr(mgncol,nlev)          ! mass weighted rain fallspeed (m/s)
         735             :   real(r8), intent(out) :: ums(mgncol,nlev)          ! mass weighted snow fallspeed (m/s)
         736             :   real(r8), intent(out) :: umg(mgncol,nlev)          ! mass weighted graupel/hail fallspeed (m/s)
         737             :   real(r8), intent(out) :: qgsedten(mgncol,nlev)     ! qg sedimentation tendency (1/s)
         738             :   real(r8), intent(out) :: qcsedten(mgncol,nlev)     ! qc sedimentation tendency (1/s)
         739             :   real(r8), intent(out) :: qisedten(mgncol,nlev)     ! qi sedimentation tendency (1/s)
         740             :   real(r8), intent(out) :: qrsedten(mgncol,nlev)     ! qr sedimentation tendency (1/s)
         741             :   real(r8), intent(out) :: qssedten(mgncol,nlev)     ! qs sedimentation tendency (1/s)
         742             :
         743             :   ! microphysical process rates for output (mixing ratio tendencies) (all have units of 1/s)
         744             :   real(r8), intent(out) :: pratot(mgncol,nlev)          ! accretion of cloud by rain
         745             :   real(r8), intent(out) :: prctot(mgncol,nlev)          ! autoconversion of cloud to rain
         746             :   real(r8), intent(out) :: mnuccctot(mgncol,nlev)       ! mixing ratio tend due to immersion freezing
         747             :   real(r8), intent(out) :: mnuccttot(mgncol,nlev)       ! mixing ratio tend due to contact freezing
         748             :   real(r8), intent(out) :: msacwitot(mgncol,nlev)       ! mixing ratio tend due to H-M splintering
         749             :   real(r8), intent(out) :: psacwstot(mgncol,nlev)       ! collection of cloud water by snow
         750             :   real(r8), intent(out) :: bergstot(mgncol,nlev)        ! bergeron process on snow
         751             :   real(r8), intent(out) :: bergtot(mgncol,nlev)         ! bergeron process on cloud ice
         752             :   real(r8), intent(out) :: melttot(mgncol,nlev)         ! melting of cloud ice
         753             :   real(r8), intent(out) :: meltstot(mgncol,nlev)         ! melting of snow
         754             :   real(r8), intent(out) :: meltgtot(mgncol,nlev)         ! melting of graupel
         755             :   real(r8), intent(out) :: mnudeptot(mgncol,nlev)       ! deposition nucleation to ice
         756             :   real(r8), intent(out) :: homotot(mgncol,nlev)         ! homogeneous freezing cloud water
         757             :   real(r8), intent(out) :: qcrestot(mgncol,nlev)        ! residual cloud condensation due to removal of excess supersat
         758             :   real(r8), intent(out) :: prcitot(mgncol,nlev)         ! autoconversion of cloud ice to snow
         759             :   real(r8), intent(out) :: praitot(mgncol,nlev)         ! accretion of cloud ice by snow
         760             :   real(r8), intent(out) :: qirestot(mgncol,nlev)        ! residual ice deposition due to removal of excess supersat
         761             :   real(r8), intent(out) :: mnuccrtot(mgncol,nlev)       ! mixing ratio tendency due to heterogeneous freezing of rain to snow (1/s)
         762             :   real(r8), intent(out) :: mnuccritot(mgncol,nlev)      ! mixing ratio tendency due to heterogeneous freezing of rain to ice (1/s)
         763             :   real(r8), intent(out) :: pracstot(mgncol,nlev)        ! mixing ratio tendency due to accretion of rain by snow (1/s)
         764             :   real(r8), intent(out) :: meltsdttot(mgncol,nlev)      ! latent heating rate due to melting of snow  (W/kg)
         765             :   real(r8), intent(out) :: frzrdttot(mgncol,nlev)       ! latent heating rate due to homogeneous freezing of rain (W/kg)
         766             :   real(r8), intent(out) :: mnuccdtot(mgncol,nlev)       ! mass tendency from ice nucleation
         767             :   real(r8), intent(out) :: pracgtot(mgncol,nlev)        ! change in q collection rain by graupel  (precipf)
         768             :   real(r8), intent(out) :: psacwgtot(mgncol,nlev)       ! change in q collection droplets by graupel (lcldm)
         769             :   real(r8), intent(out) :: pgsacwtot(mgncol,nlev)       ! conversion q to graupel due to collection droplets by snow  (lcldm)
         770             :   real(r8), intent(out) :: pgracstot(mgncol,nlev)       ! conversion q to graupel due to collection rain by snow (precipf)
         771             :   real(r8), intent(out) :: prdgtot(mgncol,nlev)         ! dep of graupel (precipf)
         772             :   real(r8), intent(out) :: qmultgtot(mgncol,nlev)       ! change q due to ice mult droplets/graupel  (lcldm)
         773             :   real(r8), intent(out) :: qmultrgtot(mgncol,nlev)      ! change q due to ice mult rain/graupel (precipf)
         774             :   real(r8), intent(out) :: psacrtot(mgncol,nlev)        ! conversion due to coll of snow by rain (precipf)
         775             :   real(r8), intent(out) :: npracgtot(mgncol,nlev)       ! change n collection rain by graupel  (precipf)
         776             :   real(r8), intent(out) :: nscngtot(mgncol,nlev)        ! change n conversion to graupel due to collection droplets by snow (lcldm)
         777             :   real(r8), intent(out) :: ngracstot(mgncol,nlev)       ! change n conversion to graupel due to collection rain by snow (precipf)
         778             :   real(r8), intent(out) :: nmultgtot(mgncol,nlev)       ! ice mult due to acc droplets by graupel  (lcldm)
         779             :   real(r8), intent(out) :: nmultrgtot(mgncol,nlev)      ! ice mult due to acc rain by graupel  (precipf)
         780             :   real(r8), intent(out) :: npsacwgtot(mgncol,nlev)      ! change n collection droplets by graupel (lcldm?)
         781             :   real(r8), intent(out) :: nrout(mgncol,nlev)        ! rain number concentration (1/m3)
         782             :   real(r8), intent(out) :: nsout(mgncol,nlev)        ! snow number concentration (1/m3)
         783             :   real(r8), intent(out) :: refl(mgncol,nlev)         ! analytic radar reflectivity
         784             :   real(r8), intent(out) :: arefl(mgncol,nlev)        ! average reflectivity will zero points outside valid range
         785             :   real(r8), intent(out) :: areflz(mgncol,nlev)       ! average reflectivity in z.
         786             :   real(r8), intent(out) :: frefl(mgncol,nlev)        ! fractional occurrence of radar reflectivity
         787             :   real(r8), intent(out) :: csrfl(mgncol,nlev)        ! cloudsat reflectivity
         788             :   real(r8), intent(out) :: acsrfl(mgncol,nlev)       ! cloudsat average
         789             :   real(r8), intent(out) :: fcsrfl(mgncol,nlev)       ! cloudsat fractional occurrence of radar reflectivity
         790             :   real(r8), intent(out) :: rercld(mgncol,nlev)       ! effective radius calculation for rain + cloud
         791             :   real(r8), intent(out) :: ncai(mgncol,nlev)         ! output number conc of ice nuclei available (1/m3)
         792             :   real(r8), intent(out) :: ncal(mgncol,nlev)         ! output number conc of CCN (1/m3)
         793             :   real(r8), intent(out) :: qrout2(mgncol,nlev)       ! copy of qrout as used to compute drout2
         794             :   real(r8), intent(out) :: qsout2(mgncol,nlev)       ! copy of qsout as used to compute dsout2
         795             :   real(r8), intent(out) :: nrout2(mgncol,nlev)       ! copy of nrout as used to compute drout2
         796             :   real(r8), intent(out) :: nsout2(mgncol,nlev)       ! copy of nsout as used to compute dsout2
         797             :   real(r8), intent(out) :: drout2(mgncol,nlev)       ! mean rain particle diameter (m)
         798             :   real(r8), intent(out) :: dsout2(mgncol,nlev)       ! mean snow particle diameter (m)
         799             :   real(r8), intent(out) :: freqs(mgncol,nlev)        ! fractional occurrence of snow
         800             :   real(r8), intent(out) :: freqr(mgncol,nlev)        ! fractional occurrence of rain
         801             :   real(r8), intent(out) :: nfice(mgncol,nlev)        ! fractional occurrence of ice
         802             :   real(r8), intent(out) :: qcrat(mgncol,nlev)        ! limiter for qc process rates (1=no limit --> 0. no qc)
         803             :   real(r8), intent(out) :: qgout(mgncol,nlev)        ! graupel/hail mixing ratio (kg/kg)
         804             :   real(r8), intent(out) :: dgout(mgncol,nlev)        ! graupel/hail diameter (m)
         805             :   real(r8), intent(out) :: ngout(mgncol,nlev)        ! graupel/hail number concentration (1/m3)
         806             :   real(r8), intent(out) :: qgout2(mgncol,nlev)       ! copy of qgout as used to compute dgout2
         807             :   real(r8), intent(out) :: ngout2(mgncol,nlev)       ! copy of ngout as used to compute dgout2
         808             :   real(r8), intent(out) :: dgout2(mgncol,nlev)       ! mean graupel/hail particle diameter (m)
         809             :   real(r8), intent(out) :: freqg(mgncol,nlev)        ! fractional occurrence of graupel
         810             :
         811             :   real(r8), intent(out) :: prer_evap(mgncol,nlev)
         812             :
         813             :   character(128),   intent(out) :: errstring  ! output status (non-blank for error return)
         814             :
         815             :   ! Tendencies calculated by external schemes that can replace MG's native
         816             :   ! process tendencies.
         817             :
         818             :   ! Used with CARMA cirrus microphysics
         819             :   ! (or similar external microphysics model)
         820             :   real(r8), intent(in) :: tnd_qsnow(:,:) ! snow mass tendency (kg/kg/s)
         821             :   real(r8), intent(in) :: tnd_nsnow(:,:) ! snow number tendency (#/kg/s)
         822             :   real(r8), intent(in) :: re_ice(:,:)    ! ice effective radius (m)
         823             :
         824             :   ! From external ice nucleation.
         825             :   real(r8), intent(in) :: frzimm(:,:) ! Number tendency due to immersion freezing (1/cm3)
         826             :   real(r8), intent(in) :: frzcnt(:,:) ! Number tendency due to contact freezing (1/cm3)
         827             :   real(r8), intent(in) :: frzdep(:,:) ! Number tendency due to deposition nucleation (1/cm3)
         828             :
         829             :   ! local workspace
         830             :   ! all units mks unless otherwise stated
         831             :
         832             :   ! local copies of input variables
         833           0 :   real(r8) :: qc(mgncol,nlev)      ! cloud liquid mixing ratio (kg/kg)
         834           0 :   real(r8) :: qi(mgncol,nlev)      ! cloud ice mixing ratio (kg/kg)
         835           0 :   real(r8) :: nc(mgncol,nlev)      ! cloud liquid number concentration (1/kg)
         836           0 :   real(r8) :: ni(mgncol,nlev)      ! cloud liquid number concentration (1/kg)
         837           0 :   real(r8) :: qr(mgncol,nlev)      ! rain mixing ratio (kg/kg)
         838           0 :   real(r8) :: qs(mgncol,nlev)      ! snow mixing ratio (kg/kg)
         839           0 :   real(r8) :: nr(mgncol,nlev)      ! rain number concentration (1/kg)
         840           0 :   real(r8) :: ns(mgncol,nlev)      ! snow number concentration (1/kg)
         841           0 :   real(r8) :: qg(mgncol,nlev)      ! graupel mixing ratio (kg/kg)
         842           0 :   real(r8) :: ng(mgncol,nlev)      ! graupel number concentration (1/kg)
         843             :   real(r8) :: rhogtmp              ! hail or graupel density (kg m-3)
         844             :
         845             :   ! general purpose variables
         846             :   real(r8) :: deltat            ! sub-time step (s)
         847             :   real(r8) :: mtime             ! the assumed ice nucleation timescale
         848             :   real(r8) :: rdeltat           ! reciprocal of sub-time step (1/s)
         849             :
         850             :   ! physical properties of the air at a given point
         851           0 :   real(r8) :: rho(mgncol,nlev)    ! density (kg m-3)
         852           0 :   real(r8) :: dv(mgncol,nlev)     ! diffusivity of water vapor
         853           0 :   real(r8) :: mu(mgncol,nlev)     ! viscosity
         854           0 :   real(r8) :: sc(mgncol,nlev)     ! schmidt number
         855           0 :   real(r8) :: rhof(mgncol,nlev)   ! density correction factor for fallspeed
         856             :
         857             :   ! cloud fractions
         858           0 :   real(r8) :: precip_frac(mgncol,nlev) ! precip fraction assuming maximum overlap
         859           0 :   real(r8) :: cldm(mgncol,nlev)   ! cloud fraction
         860           0 :   real(r8) :: icldm(mgncol,nlev)  ! ice cloud fraction
         861           0 :   real(r8) :: lcldm(mgncol,nlev)  ! liq cloud fraction
         862           0 :   real(r8) :: qsfm(mgncol,nlev)   ! subgrid cloud water saturation scaling factor
         863             :
         864             :   ! mass mixing ratios
         865           0 :   real(r8) :: qcic(mgncol,nlev)   ! in-cloud cloud liquid
         866           0 :   real(r8) :: qiic(mgncol,nlev)   ! in-cloud cloud ice
         867           0 :   real(r8) :: qsic(mgncol,nlev)   ! in-precip snow
         868           0 :   real(r8) :: qric(mgncol,nlev)   ! in-precip rain
         869           0 :   real(r8) :: qgic(mgncol,nlev)   ! in-precip graupel/hail
         870             :
         871             :   ! number concentrations
         872           0 :   real(r8) :: ncic(mgncol,nlev)   ! in-cloud droplet
         873           0 :   real(r8) :: niic(mgncol,nlev)   ! in-cloud cloud ice
         874           0 :   real(r8) :: nsic(mgncol,nlev)   ! in-precip snow
         875           0 :   real(r8) :: nric(mgncol,nlev)   ! in-precip rain
         876           0 :   real(r8) :: ngic(mgncol,nlev)   ! in-precip graupel/hail
         877             :
         878             :   ! maximum allowed ni value
         879           0 :   real(r8) :: nimax(mgncol,nlev)
         880             :
         881             :   ! Size distribution parameters for:
         882             :   ! cloud ice
         883           0 :   real(r8) :: lami(mgncol,nlev)   ! slope
         884           0 :   real(r8) :: n0i(mgncol,nlev)    ! intercept
         885             :   ! cloud liquid
         886           0 :   real(r8) :: lamc(mgncol,nlev)   ! slope
         887           0 :   real(r8) :: pgam(mgncol,nlev)   ! spectral width parameter
         888             :   ! snow
         889           0 :   real(r8) :: lams(mgncol,nlev)   ! slope
         890           0 :   real(r8) :: n0s(mgncol,nlev)    ! intercept
         891             :   ! rain
         892           0 :   real(r8) :: lamr(mgncol,nlev)   ! slope
         893           0 :   real(r8) :: n0r(mgncol,nlev)    ! intercept
         894             :   ! graupel/hail
         895           0 :   real(r8) :: lamg(mgncol,nlev)   ! slope
         896           0 :   real(r8) :: n0g(mgncol,nlev)    ! intercept
         897             :   real(r8) :: bgtmp               ! tmp fall speed parameter
         898             :
         899             :   ! Rates/tendencies due to:
         900             :
         901             :   ! Instantaneous snow melting
         902           0 :   real(r8) :: minstsm(mgncol,nlev)    ! mass mixing ratio
         903           0 :   real(r8) :: ninstsm(mgncol,nlev)    ! number concentration
         904             :   ! Instantaneous graupel melting
         905           0 :   real(r8) :: minstgm(mgncol,nlev)    ! mass mixing ratio
         906           0 :   real(r8) :: ninstgm(mgncol,nlev)    ! number concentration
         907             :
         908             :   ! Instantaneous rain freezing
         909           0 :   real(r8) :: minstrf(mgncol,nlev)    ! mass mixing ratio
         910           0 :   real(r8) :: ninstrf(mgncol,nlev)    ! number concentration
         911             :
         912             :   ! deposition of cloud ice
         913           0 :   real(r8) :: vap_dep(mgncol,nlev)    ! deposition from vapor to ice PMC 12/3/12
         914             :   ! sublimation of cloud ice
         915           0 :   real(r8) :: ice_sublim(mgncol,nlev) ! sublimation from ice to vapor PMC 12/3/12
         916             :   ! ice nucleation
         917           0 :   real(r8) :: nnuccd(mgncol,nlev) ! number rate from deposition/cond.-freezing
         918           0 :   real(r8) :: mnuccd(mgncol,nlev) ! mass mixing ratio
         919             :   ! freezing of cloud water
         920           0 :   real(r8) :: mnuccc(mgncol,nlev) ! mass mixing ratio
         921           0 :   real(r8) :: nnuccc(mgncol,nlev) ! number concentration
         922             :   ! contact freezing of cloud water
         923           0 :   real(r8) :: mnucct(mgncol,nlev) ! mass mixing ratio
         924           0 :   real(r8) :: nnucct(mgncol,nlev) ! number concentration
         925             :   ! deposition nucleation in mixed-phase clouds (from external scheme)
         926           0 :   real(r8) :: mnudep(mgncol,nlev) ! mass mixing ratio
         927           0 :   real(r8) :: nnudep(mgncol,nlev) ! number concentration
         928             :   ! ice multiplication
         929           0 :   real(r8) :: msacwi(mgncol,nlev) ! mass mixing ratio
         930           0 :   real(r8) :: nsacwi(mgncol,nlev) ! number concentration
         931             :   ! autoconversion of cloud droplets
         932           0 :   real(r8) :: prc(mgncol,nlev)    ! mass mixing ratio
         933           0 :   real(r8) :: nprc(mgncol,nlev)   ! number concentration (rain)
         934           0 :   real(r8) :: nprc1(mgncol,nlev)  ! number concentration (cloud droplets)
         935             :   ! self-aggregation of snow
         936           0 :   real(r8) :: nsagg(mgncol,nlev)  ! number concentration
         937             :   ! self-collection of rain
         938           0 :   real(r8) :: nragg(mgncol,nlev)  ! number concentration
         939             :   ! collection of droplets by snow
         940           0 :   real(r8) :: psacws(mgncol,nlev)     ! mass mixing ratio
         941           0 :   real(r8) :: npsacws(mgncol,nlev)    ! number concentration
         942             :   ! collection of rain by snow
         943           0 :   real(r8) :: pracs(mgncol,nlev)  ! mass mixing ratio
         944           0 :   real(r8) :: npracs(mgncol,nlev) ! number concentration
         945             :   ! freezing of rain
         946           0 :   real(r8) :: mnuccr(mgncol,nlev) ! mass mixing ratio
         947           0 :   real(r8) :: nnuccr(mgncol,nlev) ! number concentration
         948             :   ! freezing of rain to form ice (mg add 4/26/13)
         949           0 :   real(r8) :: mnuccri(mgncol,nlev)    ! mass mixing ratio
         950           0 :   real(r8) :: nnuccri(mgncol,nlev)    ! number concentration
         951             :   ! accretion of droplets by rain
         952           0 :   real(r8) :: pra(mgncol,nlev)    ! mass mixing ratio
         953           0 :   real(r8) :: npra(mgncol,nlev)   ! number concentration
         954             :   ! autoconversion of cloud ice to snow
         955           0 :   real(r8) :: prci(mgncol,nlev)   ! mass mixing ratio
         956           0 :   real(r8) :: nprci(mgncol,nlev)  ! number concentration
         957             :   ! accretion of cloud ice by snow
         958           0 :   real(r8) :: prai(mgncol,nlev)   ! mass mixing ratio
         959           0 :   real(r8) :: nprai(mgncol,nlev)  ! number concentration
         960             :   ! evaporation of rain
         961           0 :   real(r8) :: pre(mgncol,nlev)    ! mass mixing ratio
         962             :   ! sublimation of snow
         963           0 :   real(r8) :: prds(mgncol,nlev)   ! mass mixing ratio
         964             :   ! number evaporation
         965           0 :   real(r8) :: nsubi(mgncol,nlev)  ! cloud ice
         966           0 :   real(r8) :: nsubc(mgncol,nlev)  ! droplet
         967           0 :   real(r8) :: nsubs(mgncol,nlev)  ! snow
         968           0 :   real(r8) :: nsubr(mgncol,nlev)  ! rain
         969             :   ! bergeron process
         970           0 :   real(r8) :: berg(mgncol,nlev)   ! mass mixing ratio (cloud ice)
         971           0 :   real(r8) :: bergs(mgncol,nlev)  ! mass mixing ratio (snow)
         972             :
         973             :   !graupel/hail processes
         974           0 :   real(r8) :: npracg(mgncol,nlev)  ! change n collection rain by graupel  (precipf)
         975           0 :   real(r8) :: nscng(mgncol,nlev)   ! change n conversion to graupel due to collection droplets by snow (lcldm)
         976           0 :   real(r8) :: ngracs(mgncol,nlev)  ! change n conversion to graupel due to collection rain by snow (precipf)
         977           0 :   real(r8) :: nmultg(mgncol,nlev)  ! ice mult due to acc droplets by graupel  (lcldm)
         978           0 :   real(r8) :: nmultrg(mgncol,nlev) ! ice mult due to acc rain by graupel  (precipf)
         979           0 :   real(r8) :: npsacwg(mgncol,nlev) ! change n collection droplets by graupel (lcldm)
         980             :
         981           0 :   real(r8) :: psacr(mgncol,nlev)   ! conversion due to coll of snow by rain (precipf)
         982           0 :   real(r8) :: pracg(mgncol,nlev)   ! change in q collection rain by graupel  (precipf)
         983           0 :   real(r8) :: psacwg(mgncol,nlev)  ! change in q collection droplets by graupel (lcldm)
         984           0 :   real(r8) :: pgsacw(mgncol,nlev)  ! conversion q to graupel due to collection droplets by snow  (lcldm)
         985           0 :   real(r8) :: pgracs(mgncol,nlev)  ! conversion q to graupel due to collection rain by snow (precipf)
         986           0 :   real(r8) :: prdg(mgncol,nlev)    ! dep of graupel (precipf)
         987           0 :   real(r8) :: qmultg(mgncol,nlev)  ! change q due to ice mult droplets/graupel  (lcldm)
         988           0 :   real(r8) :: qmultrg(mgncol,nlev) ! change q due to ice mult rain/graupel (precipf)
         989             :
         990             :
         991             :   ! fallspeeds
         992             :   ! number-weighted
         993           0 :   real(r8) :: uns(mgncol,nlev)    ! snow
         994           0 :   real(r8) :: unr(mgncol,nlev)    ! rain
         995           0 :   real(r8) :: ung(mgncol,nlev)    ! graupel/hail
         996             :
         997             :   ! air density corrected fallspeed parameters
         998           0 :   real(r8) :: arn(mgncol,nlev)    ! rain
         999           0 :   real(r8) :: asn(mgncol,nlev)    ! snow
        1000           0 :   real(r8) :: agn(mgncol,nlev)    ! graupel
        1001           0 :   real(r8) :: acn(mgncol,nlev)    ! cloud droplet
        1002           0 :   real(r8) :: ain(mgncol,nlev)    ! cloud ice
        1003           0 :   real(r8) :: ajn(mgncol,nlev)    ! cloud small ice
        1004             :
        1005             :   ! Mass of liquid droplets used with external heterogeneous freezing.
        1006           0 :   real(r8) :: mi0l(mgncol,nlev)
        1007             :
        1008             :   ! saturation vapor pressures
        1009           0 :   real(r8) :: esl(mgncol,nlev)    ! liquid
        1010           0 :   real(r8) :: esi(mgncol,nlev)    ! ice
        1011           0 :   real(r8) :: esnA(mgncol,nlev)   ! checking for RH after rain evap
        1012             :
        1013             :   ! saturation vapor mixing ratios
        1014           0 :   real(r8) :: qvl(mgncol,nlev)    ! liquid
        1015           0 :   real(r8) :: qvi(mgncol,nlev)    ! ice
        1016           0 :   real(r8) :: qvnA(mgncol,nlev), qvnAI(mgncol,nlev) ! checking for RH after rain evap
        1017             :
        1018             :   ! relative humidity
        1019           0 :   real(r8) :: relhum(mgncol,nlev)
        1020             :
        1021             :   ! parameters for cloud water and cloud ice sedimentation calculations
        1022           0 :   real(r8) :: fc(mgncol,nlev)
        1023           0 :   real(r8) :: fnc(mgncol,nlev)
        1024           0 :   real(r8) :: fi(mgncol,nlev)
        1025           0 :   real(r8) :: fni(mgncol,nlev)
        1026           0 :   real(r8) :: fg(mgncol,nlev)
        1027           0 :   real(r8) :: fng(mgncol,nlev)
        1028           0 :   real(r8) :: fr(mgncol,nlev)
        1029           0 :   real(r8) :: fnr(mgncol,nlev)
        1030           0 :   real(r8) :: fs(mgncol,nlev)
        1031           0 :   real(r8) :: fns(mgncol,nlev)
        1032             :
        1033             :   real(r8) :: faloutc(nlev)
        1034             :   real(r8) :: faloutnc(nlev)
        1035             :   real(r8) :: falouti(nlev)
        1036             :   real(r8) :: faloutni(nlev)
        1037             :
        1038             :   real(r8) :: faloutr(nlev)
        1039             :   real(r8) :: faloutnr(nlev)
        1040             :   real(r8) :: falouts(nlev)
        1041             :   real(r8) :: faloutns(nlev)
        1042             :
        1043           0 :   real(r8) :: rainrt(mgncol,nlev)     ! rain rate for reflectivity calculation
        1044             :
        1045             :   ! dummy variables
        1046             :   real(r8) :: dum, dum1, dum2, dum3, dum4, qtmp
        1047           0 :   real(r8) :: dum1A(mgncol,nlev), dum2A(mgncol,nlev), dum3A(mgncol,nlev)
        1048           0 :   real(r8) :: dumni0, dumni0A2D(mgncol,nlev)
        1049           0 :   real(r8) :: dumns0, dumns0A2D(mgncol,nlev)
        1050             :   ! dummies for checking RH
        1051           0 :   real(r8) :: ttmpA(mgncol,nlev), qtmpAI(mgncol,nlev)
        1052             :   ! dummies for conservation check
        1053             :   real(r8) :: ratio
        1054             :   real(r8) :: tmpfrz
        1055             :   ! dummies for in-cloud variables
        1056           0 :   real(r8) :: dumc(mgncol,nlev)   ! qc
        1057           0 :   real(r8) :: dumnc(mgncol,nlev)  ! nc
        1058           0 :   real(r8) :: dumi(mgncol,nlev)   ! qi
        1059           0 :   real(r8) :: dumni(mgncol,nlev)  ! ni
        1060           0 :   real(r8) :: dumr(mgncol,nlev)   ! rain mixing ratio
        1061           0 :   real(r8) :: dumnr(mgncol,nlev)  ! rain number concentration
        1062           0 :   real(r8) :: dums(mgncol,nlev)   ! snow mixing ratio
        1063           0 :   real(r8) :: dumns(mgncol,nlev)  ! snow number concentration
        1064           0 :   real(r8) :: dumg(mgncol,nlev)   ! graupel mixing ratio
        1065           0 :   real(r8) :: dumng(mgncol,nlev)  ! graupel number concentration
        1066             :   ! Array dummy variable
        1067           0 :   real(r8) :: dum_2D(mgncol,nlev)
        1068           0 :   real(r8) :: pdel_inv(mgncol,nlev)
        1069             :
        1070             :   ! loop array variables
        1071             :   ! "i" and "k" are column/level iterators for internal (MG) variables
        1072             :   ! "n" is used for other looping (currently just sedimentation)
        1073             :   integer i, k, n
        1074             :
        1075             :   ! number of sub-steps for loops over "n" (for sedimentation)
        1076             :   integer nstep
        1077             :   integer mdust
        1078             :   integer :: precip_frac_method
        1079             :
        1080             :   ! Varaibles to scale fall velocity between small and regular ice regimes.
        1081             :   real(r8) :: irad
        1082             :   real(r8) :: ifrac
        1083             :
        1084             :   real(r8) :: nimey  !meyers ice nucleation
        1085           0 :   real(r8) :: niact(mgncol,nlev) ! dummy for modified activation
        1086             :
        1087             :   !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        1088             :
        1089             :   ! Return error message
        1090           0 :   errstring = ' '
        1091             :
        1092             :   ! Process inputs
        1093             :
        1094             :   ! assign variable deltat to deltatin
        1095           0 :   deltat  = deltatin
        1096           0 :   rdeltat = 1._r8 / deltat
        1097             :
        1098           0 :   if (trim(micro_mg_precip_frac_method) == 'in_cloud') then
        1099             :      precip_frac_method = MG_PRECIP_FRAC_INCLOUD
        1100           0 :   else if(trim(micro_mg_precip_frac_method) == 'max_overlap') then
        1101             :      precip_frac_method = MG_PRECIP_FRAC_OVERLAP
        1102             :   endif
        1103             :
        1104             :   !===============================================
        1105             :   ! Set ice nucleation timescale to deltat before microphysics loop
        1106           0 :   mtime=deltat
        1107             :
        1108             :   !......................................................................
        1109             :   !       graupel/hail density set (Hail = 400, Graupel = 500 from M2005)
        1110           0 :   bgtmp=0._r8
        1111           0 :   rhogtmp=0._r8
        1112           0 :   if (do_hail) then
        1113           0 :      bgtmp = bh
        1114           0 :      rhogtmp = rhoh
        1115             :   end if
        1116           0 :   if (do_graupel) then
        1117           0 :      bgtmp = bg
        1118           0 :      rhogtmp = rhog
        1119             :   end if
        1120             :
        1121             :   ! set mdust as the number of dust bins for use later in contact freezing subroutine
        1122           0 :   mdust = size(rndst,3)
        1123             :
        1124             :   !$acc data copyin  (t,q,qcn,qin,ncn,nin,qrn,qsn,nrn,nsn,qgr,ngr,relvar, &
        1125             :   !$acc               accre_enhan,p,pdel,cldn,liqcldf,icecldf,qsatfac,    &
        1126             :   !$acc               naai,npccn,rndst,nacon,tnd_qsnow,tnd_nsnow,re_ice,  &
        1127             :   !$acc               frzimm,frzcnt,frzdep,mg_liq_props,mg_ice_props,     &
        1128             :   !$acc               mg_rain_props,mg_graupel_props,mg_hail_props,       &
        1129             :   !$acc               mg_snow_props)                                      &
        1130             :   !$acc      copyout (qcsinksum_rate1ord,tlat,qvlat,qctend,qitend,nctend, &
        1131             :   !$acc               nitend,qrtend,qstend,nrtend,nstend,qgtend,ngtend,   &
        1132             :   !$acc               effc,effc_fn,effi,sadice,sadsnow,prect,preci,       &
        1133             :   !$acc               nevapr,evapsnow,am_evp_st,prain,prodsnow,cmeout,    &
        1134             :   !$acc               deffi,pgamrad,lamcrad,qsout,dsout,lflx,iflx,rflx,   &
        1135             :   !$acc               sflx,gflx,qrout,reff_rain,reff_snow,reff_grau,      &
        1136             :   !$acc               qcsevap,qisevap,qvres,cmeitot,vtrmc,vtrmi,umr,ums,  &
        1137             :   !$acc               umg,qgsedten,qcsedten,qisedten,qrsedten,qssedten,   &
        1138             :   !$acc               pratot,prctot,mnuccctot,mnuccttot,msacwitot,        &
        1139             :   !$acc               psacwstot,bergstot,bergtot,melttot,meltstot,        &
        1140             :   !$acc               meltgtot,mnudeptot,homotot,qcrestot,prcitot,        &
        1141             :   !$acc               praitot,qirestot,mnuccrtot,mnuccritot,pracstot,     &
        1142             :   !$acc               meltsdttot,frzrdttot,mnuccdtot,pracgtot,psacwgtot,  &
        1143             :   !$acc               pgsacwtot,pgracstot,prdgtot,qmultgtot,qmultrgtot,   &
        1144             :   !$acc               psacrtot,npracgtot,nscngtot,ngracstot,nmultgtot,    &
        1145             :   !$acc               nmultrgtot,npsacwgtot,nrout,nsout,refl,arefl,       &
        1146             :   !$acc               areflz,frefl,csrfl,acsrfl,fcsrfl,rercld,ncai,ncal,  &
        1147             :   !$acc               qrout2,qsout2,nrout2,nsout2,drout2,dsout2,freqs,    &
        1148             :   !$acc               freqr,nfice,qcrat,qgout,dgout,ngout,qgout2,ngout2,  &
        1149             :   !$acc               dgout2,freqg,prer_evap)                             &
        1150             :   !$acc      create  (qc,qi,nc,ni,qr,qs,nr,ns,qg,ng,rho,dv,mu,sc,rhof,    &
        1151             :   !$acc               precip_frac,cldm,icldm,lcldm,qsfm,qcic,qiic,qsic,   &
        1152             :   !$acc               qric,qgic,ncic,niic,nsic,nric,ngic,nimax,lami,n0i,  &
        1153             :   !$acc               lamc,pgam,lams,n0s,lamr,n0r,lamg,n0g,minstsm,       &
        1154             :   !$acc               ninstsm,minstgm,ninstgm,minstrf,ninstrf,vap_dep,    &
        1155             :   !$acc               ice_sublim,nnuccd,mnuccd,mnuccc,nnuccc,mnucct,      &
        1156             :   !$acc               nnucct,mnudep,nnudep,msacwi,nsacwi,prc,nprc,nprc1,  &
        1157             :   !$acc               nsagg,nragg,psacws,npsacws,pracs,npracs,mnuccr,     &
        1158             :   !$acc               nnuccr,mnuccri,nnuccri,pra,npra,prci,nprci,prai,    &
        1159             :   !$acc               nprai,pre,prds,nsubi,nsubc,nsubs,nsubr,berg,bergs,  &
        1160             :   !$acc               npracg,nscng,ngracs,nmultg,nmultrg,npsacwg,psacr,   &
        1161             :   !$acc               pracg,psacwg,pgsacw,pgracs,prdg,qmultg,qmultrg,uns, &
        1162             :   !$acc               unr,ung,arn,asn,agn,acn,ain,ajn,mi0l,esl,esi,esnA,  &
        1163             :   !$acc               qvl,qvi,qvnA,qvnAI,relhum,fc,fnc,fi,fni,fg,fng,fr,  &
        1164             :   !$acc               fnr,fs,fns,faloutc,faloutnc,falouti,faloutni,       &
        1165             :   !$acc               faloutr,faloutnr,falouts,faloutns,rainrt,dum1A,     &
        1166             :   !$acc               dum2A,dum3A,dumni0A2D,dumns0A2D,ttmpA,qtmpAI,dumc,  &
        1167             :   !$acc               dumnc,dumi,dumni,dumr,dumnr,dums,dumns,dumg,dumng,  &
        1168             :   !$acc               dum_2D,pdel_inv,niact)
        1169             :
        1170             :   ! Copies of input concentrations that may be changed internally.
        1171             :
        1172             :   !$acc parallel vector_length(VLENS) default(present)
        1173             :   !$acc loop gang vector collapse(2)
        1174           0 :   do k = 1,nlev
        1175           0 :      do i = 1,mgncol
        1176           0 :         qc(i,k) = qcn(i,k)
        1177           0 :         nc(i,k) = ncn(i,k)
        1178           0 :         qi(i,k) = qin(i,k)
        1179           0 :         ni(i,k) = nin(i,k)
        1180           0 :         qr(i,k) = qrn(i,k)
        1181           0 :         nr(i,k) = nrn(i,k)
        1182           0 :         qs(i,k) = qsn(i,k)
        1183           0 :         ns(i,k) = nsn(i,k)
        1184           0 :         qg(i,k) = qgr(i,k)
        1185           0 :         ng(i,k) = ngr(i,k)
        1186             :      end do
        1187             :   end do
        1188             :
        1189             :   ! cldn: used to set cldm, unused for subcolumns
        1190             :   ! liqcldf: used to set lcldm, unused for subcolumns
        1191             :   ! icecldf: used to set icldm, unused for subcolumns
        1192             :
        1193           0 :   if (microp_uniform) then
        1194             :      ! subcolumns, set cloud fraction variables to one
        1195             :      ! if cloud water or ice is present, if not present
        1196             :      ! set to mincld (mincld used instead of zero, to prevent
        1197             :      ! possible division by zero errors).
        1198             :
        1199             :      !$acc loop gang vector collapse(2)
        1200           0 :      do k=1,nlev
        1201           0 :        do i=1,mgncol
        1202           0 :           if (qc(i,k) >= qsmall) then
        1203           0 :              lcldm(i,k) = 1._r8
        1204             :           else
        1205           0 :              lcldm(i,k) = mincld
        1206             :           end if
        1207             :
        1208           0 :           if (qi(i,k) >= qsmall) then
        1209           0 :              icldm(i,k) = 1._r8
        1210             :           else
        1211           0 :              icldm(i,k) = mincld
        1212             :           end if
        1213             :
        1214           0 :           cldm(i,k) = max(icldm(i,k), lcldm(i,k))
        1215           0 :           qsfm(i,k) = 1._r8
        1216             :         end do
        1217             :      end do
        1218             :   else
        1219             :      ! get cloud fraction, check for minimum
        1220             :
        1221             :      !$acc loop gang vector collapse(2)
        1222           0 :      do k=1,nlev
        1223           0 :         do i=1,mgncol
        1224           0 :           cldm(i,k) = max(cldn(i,k),mincld)
        1225           0 :           lcldm(i,k) = max(liqcldf(i,k),mincld)
        1226           0 :           icldm(i,k) = max(icecldf(i,k),mincld)
        1227           0 :           qsfm(i,k) = qsatfac(i,k)
        1228             :         end do
        1229             :      end do
        1230             :   end if
        1231             :
        1232             :   ! Initialize local variables
        1233             :
        1234             :   ! local physical properties
        1235             :
        1236             :   !$acc loop gang vector collapse(2)
        1237           0 :   do k=1,nlev
        1238           0 :      do i=1,mgncol
        1239           0 :         rho(i,k) = p(i,k)/(r*t(i,k))
        1240           0 :         dv(i,k) = 8.794E-5_r8 * t(i,k)**1.81_r8 / p(i,k)
        1241           0 :         mu(i,k) = 1.496E-6_r8 * t(i,k)**1.5_r8 / (t(i,k) + 120._r8)
        1242           0 :         sc(i,k) = mu(i,k)/(rho(i,k)*dv(i,k))
        1243             :
        1244             :         ! air density adjustment for fallspeed parameters
        1245             :         ! includes air density correction factor to the
        1246             :         ! power of 0.54 following Heymsfield and Bansemer 2007
        1247             :
        1248           0 :         rhof(i,k)=(rhosu/rho(i,k))**0.54_r8
        1249             :
        1250           0 :         arn(i,k)=ar*rhof(i,k)
        1251           0 :         asn(i,k)=as*rhof(i,k)
        1252             :         ! Hail use ah*rhof graupel use ag*rhof
        1253             :         ! Note that do_hail and do_graupel can't both be true
        1254           0 :         if (do_hail) then
        1255           0 :            agn(i,k) = ah*rhof(i,k)
        1256             :         end if
        1257           0 :         if (do_graupel) then
        1258           0 :            agn(i,k) = ag*rhof(i,k)
        1259             :         end if
        1260           0 :         acn(i,k)=g*rhow/(18._r8*mu(i,k))
        1261           0 :         ain(i,k)=ai*(rhosu/rho(i,k))**0.35_r8
        1262           0 :         ajn(i,k)=aj*(rhosu/rho(i,k))**0.35_r8
        1263             :      end do
        1264             :   end do
        1265             :   !$acc end parallel
        1266             :
        1267             :   !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        1268             :   ! Get humidity and saturation vapor pressures
        1269             :
        1270           0 :   call qsat_water(t, p, esl, qvl, mgncol*nlev)
        1271           0 :   call qsat_ice(t, p, esi, qvi, mgncol*nlev)
        1272             :
        1273             :   !$acc parallel vector_length(VLENS) default(present)
        1274             :   !$acc loop gang vector collapse(2)
        1275           0 :   do k=1,nlev
        1276           0 :      do i=1,mgncol
        1277             :         ! make sure when above freezing that esi=esl, not active yet
        1278           0 :         if (t(i,k) >= tmelt) then
        1279           0 :            esi(i,k)=esl(i,k)
        1280           0 :            qvi(i,k)=qvl(i,k)
        1281             :         else
        1282             :            ! Scale the water saturation values to reflect subgrid scale
        1283             :            ! ice cloud fraction, where ice clouds begin forming at a
        1284             :            ! gridbox average relative humidity of rhmini (not 1).
        1285             :            !
        1286             :            ! NOTE: For subcolumns and other non-subgrid clouds, qsfm willi
        1287             :            ! be 1.
        1288           0 :            qvi(i,k) = qsfm(i,k) * qvi(i,k)
        1289           0 :            esi(i,k) = qsfm(i,k) * esi(i,k)
        1290           0 :            qvl(i,k) = qsfm(i,k) * qvl(i,k)
        1291           0 :            esl(i,k) = qsfm(i,k) * esl(i,k)
        1292             :         end if
        1293           0 :         relhum(i,k) = q(i,k) / max(qvl(i,k), qsmall)
        1294             :
        1295             : !  Adjust NAAI with meyers ice nucleation for 0 > T > -37
        1296           0 :         if ((t(i,k).gt.tmelt-37._r8.and.t(i,k).lt.tmelt).and.icenuc_use_meyers) then
        1297           0 :            nimey=1.e3_r8*exp(12.96_r8*(esl(i,k)-esi(i,k))/esi(i,k) - 0.639_r8)
        1298           0 :            niact(i,k)=max(naai(i,k),nimey/rho(i,k))
        1299             :         else
        1300           0 :            niact(i,k)=naai(i,k)
        1301             :         end if
        1302             :
        1303             :      end do
        1304             :   end do
        1305             :
        1306             :   ! initialize microphysics output
        1307             :
        1308             :   !$acc loop gang vector collapse(2)
        1309           0 :   do k=1,nlev
        1310           0 :      do i=1,mgncol
        1311           0 :         qcsevap(i,k)            = 0._r8
        1312           0 :         qisevap(i,k)            = 0._r8
        1313           0 :         qvres(i,k)              = 0._r8
        1314           0 :         cmeitot(i,k)            = 0._r8
        1315           0 :         vtrmc(i,k)              = 0._r8
        1316           0 :         vtrmi(i,k)              = 0._r8
        1317           0 :         qcsedten(i,k)           = 0._r8
        1318           0 :         qisedten(i,k)           = 0._r8
        1319           0 :         qrsedten(i,k)           = 0._r8
        1320           0 :         qssedten(i,k)           = 0._r8
        1321           0 :         qgsedten(i,k)           = 0._r8
        1322             :
        1323           0 :         pratot(i,k)             = 0._r8
        1324           0 :         prctot(i,k)             = 0._r8
        1325           0 :         mnuccctot(i,k)          = 0._r8
        1326           0 :         mnuccttot(i,k)          = 0._r8
        1327           0 :         msacwitot(i,k)          = 0._r8
        1328           0 :         psacwstot(i,k)          = 0._r8
        1329           0 :         bergstot(i,k)           = 0._r8
        1330           0 :         bergtot(i,k)            = 0._r8
        1331           0 :         melttot(i,k)            = 0._r8
        1332             :
        1333           0 :         mnudeptot(i,k)          = 0._r8
        1334           0 :         meltstot(i,k)           = 0._r8
        1335           0 :         meltgtot(i,k)           = 0._r8
        1336           0 :         homotot(i,k)            = 0._r8
        1337           0 :         qcrestot(i,k)           = 0._r8
        1338           0 :         prcitot(i,k)            = 0._r8
        1339           0 :         praitot(i,k)            = 0._r8
        1340           0 :         qirestot(i,k)           = 0._r8
        1341           0 :         mnuccrtot(i,k)          = 0._r8
        1342           0 :         mnuccritot(i,k)         = 0._r8
        1343           0 :         pracstot(i,k)           = 0._r8
        1344           0 :         meltsdttot(i,k)         = 0._r8
        1345           0 :         frzrdttot(i,k)          = 0._r8
        1346           0 :         mnuccdtot(i,k)          = 0._r8
        1347           0 :         psacrtot(i,k)           = 0._r8
        1348           0 :         pracgtot(i,k)           = 0._r8
        1349           0 :         psacwgtot(i,k)          = 0._r8
        1350           0 :         pgsacwtot(i,k)          = 0._r8
        1351           0 :         pgracstot(i,k)          = 0._r8
        1352           0 :         prdgtot(i,k)            = 0._r8
        1353           0 :         qmultgtot(i,k)          = 0._r8
        1354           0 :         qmultrgtot(i,k)         = 0._r8
        1355           0 :         npracgtot(i,k)          = 0._r8
        1356           0 :         nscngtot(i,k)           = 0._r8
        1357           0 :         ngracstot(i,k)          = 0._r8
        1358           0 :         nmultgtot(i,k)          = 0._r8
        1359           0 :         nmultrgtot(i,k)         = 0._r8
        1360           0 :         npsacwgtot(i,k)         = 0._r8
        1361             :
        1362             : !need to zero these out to be totally switchable (for conservation)
        1363           0 :         psacr(i,k)              = 0._r8
        1364           0 :         pracg(i,k)              = 0._r8
        1365           0 :         psacwg(i,k)             = 0._r8
        1366           0 :         pgsacw(i,k)             = 0._r8
        1367           0 :         pgracs(i,k)             = 0._r8
        1368           0 :         prdg(i,k)               = 0._r8
        1369           0 :         qmultg(i,k)             = 0._r8
        1370           0 :         qmultrg(i,k)            = 0._r8
        1371           0 :         npracg(i,k)             = 0._r8
        1372           0 :         nscng(i,k)              = 0._r8
        1373           0 :         ngracs(i,k)             = 0._r8
        1374           0 :         nmultg(i,k)             = 0._r8
        1375           0 :         nmultrg(i,k)            = 0._r8
        1376           0 :         npsacwg(i,k)            = 0._r8
        1377             :      end do
        1378             :   end do
        1379             :
        1380             :   !$acc loop gang vector collapse(2)
        1381           0 :   do k=1,nlev+1
        1382           0 :      do i=1,mgncol
        1383           0 :         rflx(i,k)               = 0._r8
        1384           0 :         sflx(i,k)               = 0._r8
        1385           0 :         lflx(i,k)               = 0._r8
        1386           0 :         iflx(i,k)               = 0._r8
        1387           0 :         gflx(i,k)               = 0._r8
        1388             :      end do
        1389             :   end do
        1390             :
        1391             :   !$acc loop gang vector collapse(2)
        1392           0 :   do k=1,nlev
        1393           0 :      do i=1,mgncol
        1394             :         ! initialize precip output
        1395           0 :         qrout(i,k)              = 0._r8
        1396           0 :         qsout(i,k)              = 0._r8
        1397           0 :         nrout(i,k)              = 0._r8
        1398           0 :         nsout(i,k)              = 0._r8
        1399           0 :         qgout(i,k)              = 0._r8
        1400           0 :         ngout(i,k)              = 0._r8
        1401             :
        1402             :         ! for refl calc
        1403           0 :         rainrt(i,k)             = 0._r8
        1404             :
        1405             :         ! initialize rain size
        1406           0 :         rercld(i,k)             = 0._r8
        1407             :
        1408           0 :         qcsinksum_rate1ord(i,k) = 0._r8
        1409             :
        1410             :         ! initialize variables for trop_mozart
        1411           0 :         nevapr(i,k)             = 0._r8
        1412           0 :         prer_evap(i,k)          = 0._r8
        1413           0 :         evapsnow(i,k)           = 0._r8
        1414           0 :         am_evp_st(i,k)          = 0._r8
        1415           0 :         prain(i,k)              = 0._r8
        1416           0 :         prodsnow(i,k)           = 0._r8
        1417           0 :         cmeout(i,k)             = 0._r8
        1418             :
        1419           0 :         precip_frac(i,k)        = mincld
        1420           0 :         lamc(i,k)               = 0._r8
        1421           0 :         lamg(i,k)               = 0._r8
        1422             :
        1423             :         ! initialize microphysical tendencies
        1424           0 :         tlat(i,k)               = 0._r8
        1425           0 :         qvlat(i,k)              = 0._r8
        1426           0 :         qctend(i,k)             = 0._r8
        1427           0 :         qitend(i,k)             = 0._r8
        1428           0 :         qstend(i,k)             = 0._r8
        1429           0 :         qrtend(i,k)             = 0._r8
        1430           0 :         nctend(i,k)             = 0._r8
        1431           0 :         nitend(i,k)             = 0._r8
        1432           0 :         nrtend(i,k)             = 0._r8
        1433           0 :         nstend(i,k)             = 0._r8
        1434           0 :         qgtend(i,k)             = 0._r8
        1435           0 :         ngtend(i,k)             = 0._r8
        1436             :
        1437             :         ! initialize in-cloud and in-precip quantities to zero
        1438           0 :         qcic(i,k)               = 0._r8
        1439           0 :         qiic(i,k)               = 0._r8
        1440           0 :         qsic(i,k)               = 0._r8
        1441           0 :         qric(i,k)               = 0._r8
        1442           0 :         qgic(i,k)               = 0._r8
        1443             :
        1444           0 :         ncic(i,k)               = 0._r8
        1445           0 :         niic(i,k)               = 0._r8
        1446           0 :         nsic(i,k)               = 0._r8
        1447           0 :         nric(i,k)               = 0._r8
        1448           0 :         ngic(i,k)               = 0._r8
        1449             :      end do
        1450             :   end do
        1451             :
        1452             :   ! initialize precip at surface
        1453             :
        1454             :   !$acc loop gang vector
        1455           0 :   do i=1,mgncol
        1456           0 :      prect(i)                   = 0._r8
        1457           0 :      preci(i)                   = 0._r8
        1458             :   end do
        1459             :
        1460             :   !$acc loop gang vector collapse(2)
        1461           0 :   do k=1,nlev
        1462           0 :      do i=1,mgncol
        1463             :         ! initialize precip fallspeeds to zero
        1464           0 :         ums(i,k)                = 0._r8
        1465           0 :         uns(i,k)                = 0._r8
        1466           0 :         umr(i,k)                = 0._r8
        1467           0 :         unr(i,k)                = 0._r8
        1468           0 :         umg(i,k)                = 0._r8
        1469           0 :         ung(i,k)                = 0._r8
        1470             :
        1471             :         ! initialize limiter for output
        1472           0 :         qcrat(i,k)              = 1._r8
        1473             :
        1474             :         ! Many outputs have to be initialized here at the top to work around
        1475             :         ! ifort problems, even if they are always overwritten later.
        1476           0 :         effc(i,k)               = 10._r8
        1477           0 :         lamcrad(i,k)            = 0._r8
        1478           0 :         pgamrad(i,k)            = 0._r8
        1479           0 :         effc_fn(i,k)            = 10._r8
        1480           0 :         effi(i,k)               = 25._r8
        1481           0 :         effi(i,k) = effi(i,k)*micro_mg_effi_factor
        1482           0 :         sadice(i,k)             = 0._r8
        1483           0 :         sadsnow(i,k)            = 0._r8
        1484           0 :         deffi(i,k)              = 50._r8
        1485             :
        1486           0 :         qrout2(i,k)             = 0._r8
        1487           0 :         nrout2(i,k)             = 0._r8
        1488           0 :         drout2(i,k)             = 0._r8
        1489           0 :         qsout2(i,k)             = 0._r8
        1490           0 :         nsout2(i,k)             = 0._r8
        1491           0 :         dsout(i,k)              = 0._r8
        1492           0 :         dsout2(i,k)             = 0._r8
        1493           0 :         qgout2(i,k)             = 0._r8
        1494           0 :         ngout2(i,k)             = 0._r8
        1495           0 :         freqg(i,k)              = 0._r8
        1496           0 :         freqr(i,k)              = 0._r8
        1497           0 :         freqs(i,k)              = 0._r8
        1498             :
        1499           0 :         reff_rain(i,k)          = 0._r8
        1500           0 :         reff_snow(i,k)          = 0._r8
        1501           0 :         reff_grau(i,k)          = 0._r8
        1502             :
        1503           0 :         refl(i,k)               = -9999._r8
        1504           0 :         arefl(i,k)              = 0._r8
        1505           0 :         areflz(i,k)             = 0._r8
        1506           0 :         frefl(i,k)              = 0._r8
        1507           0 :         csrfl(i,k)              = 0._r8
        1508           0 :         acsrfl(i,k)             = 0._r8
        1509           0 :         fcsrfl(i,k)             = 0._r8
        1510             :
        1511           0 :         ncal(i,k)               = 0._r8
        1512           0 :         ncai(i,k)               = 0._r8
        1513           0 :         nfice(i,k)              = 0._r8
        1514             :      end do
        1515             :   end do
        1516             :
        1517             :   !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        1518             :   ! droplet activation
        1519             :   ! get provisional droplet number after activation. This is used for
        1520             :   ! all microphysical process calculations, for consistency with update of
        1521             :   ! droplet mass before microphysics
        1522             :
        1523             :   ! calculate potential for droplet activation if cloud water is present
        1524             :   ! tendency from activation (npccn) is read in from companion routine
        1525             :
        1526             :   ! output activated liquid and ice (convert from #/kg -> #/m3)
        1527             :   !--------------------------------------------------
        1528             :
        1529             :   !$acc loop gang vector collapse(2)
        1530           0 :   do k=1,nlev
        1531           0 :      do i=1,mgncol
        1532           0 :         if (qc(i,k) >= qsmall) then
        1533           0 :            nc(i,k) = max(nc(i,k) + npccn(i,k)*deltat, 0._r8)
        1534           0 :            ncal(i,k) = nc(i,k)*rho(i,k)/lcldm(i,k) ! sghan minimum in #/cm3
        1535             :         else
        1536           0 :            ncal(i,k) = 0._r8
        1537             :         end if
        1538             :
        1539           0 :         if (t(i,k) < icenuct) then
        1540           0 :            ncai(i,k) = naai(i,k)*rho(i,k)
        1541             :         else
        1542           0 :            ncai(i,k) = 0._r8
        1543             :         end if
        1544             :     end do
        1545             :   end do
        1546             :
        1547             :   !===============================================
        1548             :
        1549             :   ! ice nucleation if activated nuclei exist at t<-5C AND rhmini + 5%
        1550             :   !
        1551             :   ! NOTE: If using gridbox average values, condensation will not occur until rh=1,
        1552             :   ! so the threshold seems like it should be 1.05 and not rhmini + 0.05. For subgrid
        1553             :   ! clouds (using rhmini and qsfacm), the relhum has already been adjusted, and thus
        1554             :   ! the nucleation threshold should also be 1.05 and not rhmini + 0.05.
        1555             :   !-------------------------------------------------------
        1556             :
        1557           0 :   if (do_cldice) then
        1558           0 :      if (icenuc_rh_off) then
        1559             :         !$acc loop gang vector collapse(2)
        1560           0 :         do k=1,nlev
        1561           0 :            do i=1,mgncol
        1562           0 :               if (niact(i,k) > 0._r8 .and. t(i,k) < icenuct) then
        1563             :                  !if NAAI > 0. then set numice = naai (as before)
        1564             :                  !note: this is gridbox averaged
        1565           0 :                  nnuccd(i,k) = (niact(i,k)-ni(i,k)/icldm(i,k))/mtime*icldm(i,k)
        1566           0 :                  nnuccd(i,k) = max(nnuccd(i,k),0._r8)
        1567           0 :                  nimax(i,k) = naai(i,k)*icldm(i,k)
        1568             :
        1569             :                  !Calc mass of new particles using new crystal mass...
        1570             :                  !also this will be multiplied by mtime as nnuccd is...
        1571           0 :                  mnuccd(i,k) = nnuccd(i,k) * mi0
        1572             :               else
        1573           0 :                  nnuccd(i,k) = 0._r8
        1574           0 :                  nimax(i,k)  = 0._r8
        1575           0 :                  mnuccd(i,k) = 0._r8
        1576             :               end if
        1577             :            end do
        1578             :         end do
        1579             :      else
        1580             :         !$acc loop gang vector collapse(2)
        1581           0 :         do k=1,nlev
        1582           0 :            do i=1,mgncol
        1583           0 :               if (naai(i,k) > 0._r8 .and. t(i,k) < icenuct .and. &
        1584           0 :                  relhum(i,k)*esl(i,k)/esi(i,k) > 1.05_r8) then
        1585             :                  !if NAAI > 0. then set numice = naai (as before)
        1586             :                  !note: this is gridbox averaged
        1587           0 :                  nnuccd(i,k) = (naai(i,k)-ni(i,k)/icldm(i,k))/mtime*icldm(i,k)
        1588           0 :                  nnuccd(i,k) = max(nnuccd(i,k),0._r8)
        1589           0 :                  nimax(i,k) = naai(i,k)*icldm(i,k)
        1590             :
        1591             :                  !Calc mass of new particles using new crystal mass...
        1592             :                  !also this will be multiplied by mtime as nnuccd is...
        1593           0 :                  mnuccd(i,k) = nnuccd(i,k) * mi0
        1594             :               else
        1595           0 :                  nnuccd(i,k) = 0._r8
        1596           0 :                  nimax(i,k)  = 0._r8
        1597           0 :                  mnuccd(i,k) = 0._r8
        1598             :               end if
        1599             :            end do
        1600             :         end do
        1601             :     end if
        1602             : !--ag
        1603             :   end if
        1604             :
        1605             :   !=============================================================================
        1606             :   !$acc loop gang vector collapse(2) private(dum,dum1)
        1607           0 :   do k=1,nlev
        1608           0 :      do i=1,mgncol
        1609             :         ! calculate instantaneous precip processes (melting and homogeneous freezing)
        1610             :         ! melting of snow at +2 C
        1611           0 :         if (t(i,k) > snowmelt) then
        1612           0 :            if (qs(i,k) > 0._r8) then
        1613             :               ! make sure melting snow doesn't reduce temperature below threshold
        1614           0 :               dum = -xlf/cpp*qs(i,k)
        1615           0 :               if (t(i,k)+dum < snowmelt) then
        1616           0 :                  dum = (t(i,k)-snowmelt)*cpp/xlf
        1617           0 :                  dum = dum/qs(i,k)
        1618           0 :                  dum = max(0._r8,dum)
        1619           0 :                  dum = min(1._r8,dum)
        1620             :               else
        1621             :                  dum = 1._r8
        1622             :               end if
        1623             :
        1624           0 :               minstsm(i,k) = dum*qs(i,k)
        1625           0 :               ninstsm(i,k) = dum*ns(i,k)
        1626             :
        1627           0 :               dum1=-xlf*minstsm(i,k)*rdeltat
        1628           0 :               tlat(i,k)=tlat(i,k)+dum1
        1629           0 :               meltsdttot(i,k)=meltsdttot(i,k) + dum1
        1630           0 :               meltstot(i,k)=minstsm(i,k)*rdeltat
        1631             :
        1632           0 :               qs(i,k) = max(qs(i,k) - minstsm(i,k), 0._r8)
        1633           0 :               ns(i,k) = max(ns(i,k) - ninstsm(i,k), 0._r8)
        1634           0 :               qr(i,k) = max(qr(i,k) + minstsm(i,k), 0._r8)
        1635           0 :               nr(i,k) = max(nr(i,k) + ninstsm(i,k), 0._r8)
        1636             :            end if
        1637             :         end if
        1638             :
        1639             :      end do
        1640             :   end do
        1641             :
        1642             :   ! melting of graupel at +2 C
        1643             :
        1644             :   !$acc loop gang vector collapse(2) private(dum,dum1)
        1645           0 :   do k=1,nlev
        1646           0 :      do i=1,mgncol
        1647             :
        1648           0 :         if (t(i,k) > snowmelt) then
        1649           0 :            if (qg(i,k) > 0._r8) then
        1650             :
        1651             :               ! make sure melting graupel doesn't reduce temperature below threshold
        1652           0 :               dum = -xlf/cpp*qg(i,k)
        1653           0 :               if (t(i,k)+dum < snowmelt) then
        1654           0 :                  dum = (t(i,k)-snowmelt)*cpp/xlf
        1655           0 :                  dum = dum/qg(i,k)
        1656           0 :                  dum = max(0._r8,dum)
        1657           0 :                  dum = min(1._r8,dum)
        1658             :               else
        1659             :                  dum = 1._r8
        1660             :               end if
        1661             :
        1662           0 :               minstgm(i,k) = dum*qg(i,k)
        1663           0 :               ninstgm(i,k) = dum*ng(i,k)
        1664             :
        1665           0 :               dum1=-xlf*minstgm(i,k)*rdeltat
        1666           0 :               tlat(i,k)=tlat(i,k)+dum1
        1667           0 :               meltsdttot(i,k)=meltsdttot(i,k) + dum1
        1668           0 :               meltgtot(i,k)=minstgm(i,k)*rdeltat
        1669             :
        1670           0 :               qg(i,k) = max(qg(i,k) - minstgm(i,k), 0._r8)
        1671           0 :               ng(i,k) = max(ng(i,k) - ninstgm(i,k), 0._r8)
        1672           0 :               qr(i,k) = max(qr(i,k) + minstgm(i,k), 0._r8)
        1673           0 :               nr(i,k) = max(nr(i,k) + ninstgm(i,k), 0._r8)
        1674             :            end if
        1675             :         end if
        1676             :
        1677             :      end do
        1678             :   end do
        1679             :
        1680             :   !$acc loop gang vector collapse(2) private(dum,dum1)
        1681           0 :   do k=1,nlev
        1682           0 :     do i=1,mgncol
        1683             :         ! freezing of rain at -5 C
        1684             :
        1685           0 :         if (t(i,k) < rainfrze) then
        1686             :
        1687           0 :            if (qr(i,k) > 0._r8) then
        1688             :
        1689             :               ! make sure freezing rain doesn't increase temperature above threshold
        1690           0 :               dum = xlf/cpp*qr(i,k)
        1691           0 :               if (t(i,k)+dum > rainfrze) then
        1692           0 :                  dum = -(t(i,k)-rainfrze)*cpp/xlf
        1693           0 :                  dum = dum/qr(i,k)
        1694           0 :                  dum = max(0._r8,dum)
        1695           0 :                  dum = min(1._r8,dum)
        1696             :               else
        1697             :                  dum = 1._r8
        1698             :               end if
        1699             :
        1700           0 :               minstrf(i,k) = dum*qr(i,k)
        1701           0 :               ninstrf(i,k) = dum*nr(i,k)
        1702             :
        1703             :               ! heating tendency
        1704           0 :               dum1 = xlf*minstrf(i,k)*rdeltat
        1705           0 :               tlat(i,k)=tlat(i,k)+dum1
        1706           0 :               frzrdttot(i,k)=frzrdttot(i,k) + dum1
        1707             :
        1708           0 :               qr(i,k) = max(qr(i,k) - minstrf(i,k), 0._r8)
        1709           0 :               nr(i,k) = max(nr(i,k) - ninstrf(i,k), 0._r8)
        1710             :
        1711             :               ! freeze rain to graupel not snow.
        1712           0 :               if(do_hail.or.do_graupel) then
        1713           0 :                  qg(i,k) = max(qg(i,k) + minstrf(i,k), 0._r8)
        1714           0 :                  ng(i,k) = max(ng(i,k) + ninstrf(i,k), 0._r8)
        1715             :               else
        1716           0 :                  qs(i,k) = max(qs(i,k) + minstrf(i,k), 0._r8)
        1717           0 :                  ns(i,k) = max(ns(i,k) + ninstrf(i,k), 0._r8)
        1718             :               end if
        1719             :            end if
        1720             :         end if
        1721             :      end do
        1722             :   end do
        1723             :
        1724             :   !$acc loop gang vector collapse(2)
        1725           0 :   do k=1,nlev
        1726           0 :     do i=1,mgncol
        1727             :         ! obtain in-cloud values of cloud water/ice mixing ratios and number concentrations
        1728             :         !-------------------------------------------------------
        1729             :         ! for microphysical process calculations
        1730             :         ! units are kg/kg for mixing ratio, 1/kg for number conc
        1731             :
        1732           0 :         if (qc(i,k).ge.qsmall) then
        1733             :            ! limit in-cloud values to 0.005 kg/kg
        1734           0 :            qcic(i,k)=min(qc(i,k)/lcldm(i,k),5.e-3_r8)
        1735           0 :            ncic(i,k)=max(nc(i,k)/lcldm(i,k),0._r8)
        1736             :
        1737             :            ! specify droplet concentration
        1738           0 :            if (nccons) then
        1739           0 :               ncic(i,k)=ncnst/rho(i,k)
        1740             :            end if
        1741             :         else
        1742           0 :            qcic(i,k)=0._r8
        1743           0 :            ncic(i,k)=0._r8
        1744             :         end if
        1745             :
        1746           0 :         if (qi(i,k).ge.qsmall) then
        1747             :            ! limit in-cloud values to 0.005 kg/kg
        1748           0 :            qiic(i,k)=min(qi(i,k)/icldm(i,k),5.e-3_r8)
        1749           0 :            niic(i,k)=max(ni(i,k)/icldm(i,k),0._r8)
        1750             :
        1751             :            ! switch for specification of cloud ice number
        1752           0 :            if (nicons) then
        1753           0 :               niic(i,k)=ninst/rho(i,k)
        1754             :            end if
        1755             :         else
        1756           0 :            qiic(i,k)=0._r8
        1757           0 :            niic(i,k)=0._r8
        1758             :         end if
        1759             :
        1760             :      end do
        1761             :   end do
        1762             :
        1763             :   !========================================================================
        1764             :
        1765             :   ! for sub-columns cldm has already been set to 1 if cloud
        1766             :   ! water or ice is present, so precip_frac will be correctly set below
        1767             :   ! and nothing extra needs to be done here
        1768             :
        1769             :   !$acc loop gang vector collapse(2)
        1770           0 :   do k=1,nlev
        1771           0 :      do i=1,mgncol
        1772           0 :         precip_frac(i,k) = cldm(i,k)
        1773             :      end do
        1774             :   end do
        1775             :   !$acc end parallel
        1776             :
        1777             :   !$acc parallel vector_length(VLENS) default(present)
        1778           0 :   if (precip_frac_method == MG_PRECIP_FRAC_INCLOUD) then
        1779             :      !$acc loop seq
        1780           0 :      do k=2,nlev
        1781             :         !$acc loop gang vector
        1782           0 :         do i=1,mgncol
        1783           0 :            if (qc(i,k) < qsmall .and. qi(i,k) < qsmall) then
        1784           0 :               precip_frac(i,k) = precip_frac(i,k-1)
        1785             :            end if
        1786             :        end do
        1787             :      end do
        1788           0 :   else if (precip_frac_method == MG_PRECIP_FRAC_OVERLAP) then
        1789             :      ! calculate precip fraction based on maximum overlap assumption
        1790             :
        1791             :      ! if rain or snow mix ratios are smaller than threshold,
        1792             :      ! then leave precip_frac as cloud fraction at current level
        1793             :
        1794             :      !$acc loop seq
        1795           0 :      do k=2,nlev
        1796             :         !$acc loop gang vector
        1797           0 :         do i=1,mgncol
        1798           0 :            if (qr(i,k-1) >= qsmall .or. qs(i,k-1) >= qsmall .or. qg(i,k-1) >= qsmall) then
        1799           0 :               precip_frac(i,k)=max(precip_frac(i,k-1),precip_frac(i,k))
        1800             :            end if
        1801             :         end do
        1802             :      end do
        1803             :
        1804             :   end if
        1805             :   !$acc end parallel
        1806             :
        1807             :   !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        1808             :   ! get size distribution parameters based on in-cloud cloud water
        1809             :   ! these calculations also ensure consistency between number and mixing ratio
        1810             :   !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        1811             :
        1812             :   ! cloud liquid
        1813             :   !-------------------------------------------
        1814           0 :   call size_dist_param_liq(mg_liq_props, qcic, ncic, rho, pgam, lamc, mgncol, nlev)
        1815             :
        1816             :   !========================================================================
        1817             :   ! autoconversion of cloud liquid water to rain
        1818             :   ! formula from Khrouditnov and Kogan (2000), modified for sub-grid distribution of qc
        1819             :   ! minimum qc of 1 x 10^-8 prevents floating point error
        1820             :
        1821           0 :   if (.not. do_sb_physics) then
        1822           0 :     call kk2000_liq_autoconversion(microp_uniform, qcic, ncic, rho, relvar, prc, nprc, nprc1, micro_mg_autocon_fact, micro_mg_autocon_nd_exp, micro_mg_autocon_lwp_exp, mgncol*nlev)
        1823             :   end if
        1824             :
        1825             :   !$acc parallel vector_length(VLENS) default(present)
        1826             :   !$acc loop gang vector collapse(2)
        1827           0 :   do k=1,nlev
        1828           0 :      do i=1,mgncol
        1829             :         ! assign qric based on prognostic qr, using assumed precip fraction
        1830             :         ! note: this could be moved above for consistency with qcic and qiic calculations
        1831           0 :         qric(i,k) = qr(i,k)/precip_frac(i,k)
        1832           0 :         nric(i,k) = nr(i,k)/precip_frac(i,k)
        1833             :
        1834             :         ! limit in-precip mixing ratios to 10 g/kg
        1835           0 :         qric(i,k)=min(qric(i,k),0.01_r8)
        1836             :
        1837             :         ! add autoconversion to precip from above to get provisional rain mixing ratio
        1838             :         ! and number concentration (qric and nric)
        1839             :
        1840           0 :         if (qric(i,k).lt.qsmall) then
        1841           0 :            qric(i,k)=0._r8
        1842           0 :            nric(i,k)=0._r8
        1843             :         end if
        1844             :
        1845             :         ! make sure number concentration is a positive number to avoid
        1846             :         ! taking root of negative later
        1847             :
        1848           0 :         nric(i,k)=max(nric(i,k),0._r8)
        1849             :      end do
        1850             :   end do
        1851             :   !$acc end parallel
        1852             :
        1853             :   ! Get size distribution parameters for cloud ice
        1854           0 :   call size_dist_param_basic(mg_ice_props, qiic, niic, lami, mgncol, nlev, n0=n0i)
        1855             :
        1856             :   ! Alternative autoconversion
        1857           0 :   if (do_sb_physics) then
        1858           0 :      call sb2001v2_liq_autoconversion(pgam, qcic, ncic, qric, rho, relvar, prc, nprc, nprc1, mgncol*nlev)
        1859             :   end if
        1860             :
        1861             :   !.......................................................................
        1862             :   ! Autoconversion of cloud ice to snow
        1863             :   ! similar to Ferrier (1994)
        1864           0 :   if (do_cldice) then
        1865           0 :      call ice_autoconversion(t, qiic, lami, n0i, dcs, prci, nprci, mgncol*nlev)
        1866             :   else
        1867             :      ! Add in the particles that we have already converted to snow, and
        1868             :      ! don't do any further autoconversion of ice.
        1869             :
        1870             :      !$acc parallel vector_length(VLENS) default(present)
        1871             :      !$acc loop gang vector collapse(2)
        1872           0 :      do k=1,nlev
        1873           0 :         do i=1,mgncol
        1874           0 :            prci(i,k)  = tnd_qsnow(i,k) / cldm(i,k)
        1875           0 :            nprci(i,k) = tnd_nsnow(i,k) / cldm(i,k)
        1876             :         end do
        1877             :      end do
        1878             :      !$acc end parallel
        1879             :   end if
        1880             :
        1881             :   ! note, currently we don't have this
        1882             :   ! inside the do_cldice block, should be changed later
        1883             :   ! assign qsic based on prognostic qs, using assumed precip fraction
        1884             :
        1885             :   !$acc parallel vector_length(VLENS) default(present)
        1886             :   !$acc loop gang vector collapse(2)
        1887           0 :   do k=1,nlev
        1888           0 :      do i=1,mgncol
        1889           0 :         qsic(i,k) = qs(i,k)/precip_frac(i,k)
        1890           0 :         nsic(i,k) = ns(i,k)/precip_frac(i,k)
        1891             :
        1892             :         ! limit in-precip mixing ratios to 10 g/kg
        1893           0 :         qsic(i,k)=min(qsic(i,k),0.01_r8)
        1894             :
        1895             :         ! if precip mix ratio is zero so should number concentration
        1896           0 :         if (qsic(i,k) < qsmall) then
        1897           0 :            qsic(i,k)=0._r8
        1898           0 :            nsic(i,k)=0._r8
        1899             :         end if
        1900             :
        1901             :         ! make sure number concentration is a positive number to avoid
        1902             :         ! taking root of negative later
        1903           0 :         nsic(i,k)=max(nsic(i,k),0._r8)
        1904             :
        1905             :         ! also do this for graupel, which is assumed to be 'precip_frac'
        1906           0 :         qgic(i,k) = qg(i,k)/precip_frac(i,k)
        1907           0 :         ngic(i,k) = ng(i,k)/precip_frac(i,k)
        1908             :
        1909             :         ! limit in-precip mixing ratios to 10 g/kg
        1910           0 :         qgic(i,k)=min(qgic(i,k),0.01_r8)
        1911             :
        1912             :         ! if precip mix ratio is zero so should number concentration
        1913           0 :         if (qgic(i,k) < qsmall) then
        1914           0 :            qgic(i,k)=0._r8
        1915           0 :            ngic(i,k)=0._r8
        1916             :         end if
        1917             :
        1918             :         ! make sure number concentration is a positive number to avoid
        1919             :         ! taking root of negative later
        1920           0 :         ngic(i,k)=max(ngic(i,k),0._r8)
        1921             :      end do
        1922             :   end do
        1923             :   !$acc end parallel
        1924             :
        1925             :   !.......................................................................
        1926             :   ! get size distribution parameters for precip
        1927             :   !......................................................................
        1928             :   ! rain
        1929           0 :   call size_dist_param_basic(mg_rain_props, qric, nric, lamr, mgncol, nlev, n0=n0r)
        1930             :
        1931             :   !$acc parallel vector_length(VLENS) default(present)
        1932             :   !$acc loop gang vector collapse(2)
        1933           0 :   do k=1,nlev
        1934           0 :      do i=1,mgncol
        1935           0 :         if (lamr(i,k) >= qsmall) then
        1936           0 :            dum_2D(i,k)= lamr(i,k)**br
        1937             :            ! provisional rain number and mass weighted mean fallspeed (m/s)
        1938           0 :            unr(i,k) = min(arn(i,k)*gamma_br_plus1/dum_2D(i,k),9.1_r8*rhof(i,k))
        1939           0 :            umr(i,k) = min(arn(i,k)*gamma_br_plus4/(6._r8*dum_2D(i,k)),9.1_r8*rhof(i,k))
        1940             :         else
        1941           0 :            umr(i,k) = 0._r8
        1942           0 :            unr(i,k) = 0._r8
        1943             :         end if
        1944             :      end do
        1945             :   end do
        1946             :   !$acc end parallel
        1947             :
        1948             :   !......................................................................
        1949             :   ! snow
        1950           0 :   call size_dist_param_basic(mg_snow_props, qsic, nsic, lams, mgncol, nlev, n0=n0s)
        1951             :
        1952             :   !$acc parallel vector_length(VLENS) default(present)
        1953             :   !$acc loop gang vector collapse(2)
        1954           0 :   do k=1,nlev
        1955           0 :      do i=1,mgncol
        1956           0 :         if (ifs_sed) then
        1957           0 :            if (lams(i,k) > 0._r8) then
        1958           0 :               ums(i,k) = 1._r8
        1959           0 :               uns(i,k) = 1._r8
        1960             :            else
        1961           0 :               ums(i,k) = 0._r8
        1962           0 :               uns(i,k) = 0._r8
        1963             :            end if
        1964             :         else
        1965           0 :            if (lams(i,k) > 0._r8) then
        1966           0 :               dum_2D(i,k) = lams(i,k)**bs
        1967             :               ! provisional snow number and mass weighted mean fallspeed (m/s)
        1968           0 :               ums(i,k) = min(asn(i,k)*gamma_bs_plus4/(6._r8*dum_2D(i,k)),1.2_r8*rhof(i,k))
        1969           0 :               ums(i,k) = ums(i,k)*micro_mg_vtrmi_factor
        1970           0 :               uns(i,k) = min(asn(i,k)*gamma_bs_plus1/dum_2D(i,k),1.2_r8*rhof(i,k))
        1971             :            else
        1972           0 :               ums(i,k) = 0._r8
        1973           0 :               uns(i,k) = 0._r8
        1974             :            end if
        1975             :         end if
        1976             :      end do
        1977             :   end do
        1978             :   !$acc end parallel
        1979             :
        1980             :   !  graupel/hail size distributions and properties
        1981             :
        1982           0 :   if (do_hail) then
        1983           0 :      call size_dist_param_basic(mg_hail_props, qgic, ngic, lamg, mgncol, nlev, n0=n0g)
        1984             :   end if
        1985           0 :   if (do_graupel) then
        1986           0 :      call size_dist_param_basic(mg_graupel_props, qgic, ngic, lamg, mgncol, nlev, n0=n0g)
        1987             :   end if
        1988             :
        1989             :   !$acc parallel vector_length(VLENS) default(present)
        1990             :   !$acc loop gang vector collapse(2)
        1991           0 :   do k=1,nlev
        1992           0 :      do i=1,mgncol
        1993           0 :         if (lamg(i,k) > 0._r8) then
        1994           0 :            dum_2D(i,k) = lamg(i,k)**bgtmp
        1995             :            ! provisional graupel/hail number and mass weighted mean fallspeed (m/s)
        1996           0 :            umg(i,k) = min(agn(i,k)*gamma_bg_plus4/(6._r8*dum_2D(i,k)),20._r8*rhof(i,k))
        1997           0 :            ung(i,k) = min(agn(i,k)*gamma_bg_plus1/dum_2D(i,k),20._r8*rhof(i,k))
        1998             :         else
        1999           0 :            umg(i,k) = 0._r8
        2000           0 :            ung(i,k) = 0._r8
        2001             :         end if
        2002             :      end do
        2003             :   end do
        2004             :   !$acc end parallel
        2005             :
        2006           0 :   if (do_cldice) then
        2007           0 :      if (.not. use_hetfrz_classnuc) then
        2008             :         ! heterogeneous freezing of cloud water
        2009             :         !----------------------------------------------
        2010           0 :         call immersion_freezing(microp_uniform, t, pgam, lamc, qcic, ncic, relvar, mnuccc, nnuccc, mgncol*nlev)
        2011             :
        2012             :         ! make sure number of droplets frozen does not exceed available ice nuclei concentration
        2013             :         ! this prevents 'runaway' droplet freezing
        2014             :
        2015             :         !$acc parallel vector_length(VLENS) default(present)
        2016             :         !$acc loop gang vector collapse(2)
        2017           0 :         do k=1,nlev
        2018           0 :            do i=1,mgncol
        2019           0 :               if (qcic(i,k).ge.qsmall .and. t(i,k).lt.269.15_r8 .and. &
        2020             :                   nnuccc(i,k)*lcldm(i,k).gt.nnuccd(i,k)) then
        2021             :                   ! scale mixing ratio of droplet freezing with limit
        2022           0 :                   mnuccc(i,k)=mnuccc(i,k)*(nnuccd(i,k)/(nnuccc(i,k)*lcldm(i,k)))
        2023           0 :                   nnuccc(i,k)=nnuccd(i,k)/lcldm(i,k)
        2024             :               end if
        2025           0 :               mnudep(i,k)=0._r8
        2026           0 :               nnudep(i,k)=0._r8
        2027             :            end do
        2028             :         end do
        2029             :         !$acc end parallel
        2030             :
        2031             :         call contact_freezing(microp_uniform, t, p, rndst, nacon, pgam, lamc, qcic, ncic, &
        2032           0 :                               relvar, mnucct, nnucct, mgncol*nlev, mdust)
        2033             :      else
        2034             :         ! Mass of droplets frozen is the average droplet mass, except
        2035             :         ! with two limiters: concentration must be at least 1/cm^3, and
        2036             :         ! mass must be at least the minimum defined above.
        2037             :
        2038             :         !$acc parallel vector_length(VLENS) default(present)
        2039             :         !$acc loop gang vector collapse(2)
        2040           0 :         do k=1,nlev
        2041           0 :            do i=1,mgncol
        2042           0 :               mi0l(i,k) = qcic(i,k)/max(ncic(i,k), 1.0e6_r8/rho(i,k))
        2043           0 :               mi0l(i,k) = max(mi0l_min, mi0l(i,k))
        2044           0 :               if (qcic(i,k) >= qsmall) then
        2045           0 :                  nnuccc(i,k) = frzimm(i,k)*1.0e6_r8/rho(i,k)
        2046           0 :                  mnuccc(i,k) = nnuccc(i,k)*mi0l(i,k)
        2047           0 :                  nnucct(i,k) = frzcnt(i,k)*1.0e6_r8/rho(i,k)
        2048           0 :                  mnucct(i,k) = nnucct(i,k)*mi0l(i,k)
        2049           0 :                  nnudep(i,k) = frzdep(i,k)*1.0e6_r8/rho(i,k)
        2050           0 :                  mnudep(i,k) = nnudep(i,k)*mi0
        2051             :               else
        2052           0 :                  nnuccc(i,k) = 0._r8
        2053           0 :                  mnuccc(i,k) = 0._r8
        2054           0 :                  nnucct(i,k) = 0._r8
        2055           0 :                  mnucct(i,k) = 0._r8
        2056           0 :                  nnudep(i,k) = 0._r8
        2057           0 :                  mnudep(i,k) = 0._r8
        2058             :               end if
        2059             :            end do
        2060             :         end do
        2061             :         !$acc end parallel
        2062             :      end if
        2063             :   else
        2064             :      !$acc parallel vector_length(VLENS) default(present)
        2065             :      !$acc loop gang vector collapse(2)
        2066           0 :      do k=1,nlev
        2067           0 :         do i=1,mgncol
        2068           0 :            mnuccc(i,k)=0._r8
        2069           0 :            nnuccc(i,k)=0._r8
        2070           0 :            mnucct(i,k)=0._r8
        2071           0 :            nnucct(i,k)=0._r8
        2072           0 :            mnudep(i,k)=0._r8
        2073           0 :            nnudep(i,k)=0._r8
        2074             :         end do
        2075             :      end do
        2076             :      !$acc end parallel
        2077             :   end if
        2078             :
        2079           0 :   call snow_self_aggregation(t, rho, asn, rhosn, qsic, nsic, nsagg, mgncol*nlev)
        2080             :
        2081             :   call accrete_cloud_water_snow(t, rho, asn, uns, mu, qcic, ncic, qsic, pgam, &
        2082           0 :                                 lamc, lams, n0s, psacws, npsacws, mgncol*nlev)
        2083             :
        2084           0 :   psacws = psacws*micro_mg_iaccr_factor
        2085           0 :   npsacws = npsacws*micro_mg_iaccr_factor
        2086           0 :   if (do_cldice) then
        2087           0 :      call secondary_ice_production(t, psacws, msacwi, nsacwi, mgncol*nlev)
        2088             :   else
        2089             :      !$acc parallel vector_length(VLENS) default(present)
        2090             :      !$acc loop gang vector collapse(2)
        2091           0 :      do k=1,nlev
        2092           0 :         do i=1,mgncol
        2093           0 :            nsacwi(i,k) = 0.0_r8
        2094           0 :            msacwi(i,k) = 0.0_r8
        2095             :         end do
        2096             :      end do
        2097             :      !$acc end parallel
        2098             :   end if
        2099             :
        2100             :   call accrete_rain_snow(t, rho, umr, ums, unr, uns, qric, qsic, lamr, &
        2101           0 :                          n0r, lams, n0s, pracs, npracs, mgncol*nlev)
        2102             :
        2103           0 :   call heterogeneous_rain_freezing(t, qric, nric, lamr, mnuccr, nnuccr, mgncol*nlev)
        2104             :
        2105           0 :   if (do_sb_physics) then
        2106           0 :      call sb2001v2_accre_cld_water_rain(qcic, ncic, qric, rho, relvar, pra, npra, mgncol*nlev)
        2107             :   else
        2108           0 :      call accrete_cloud_water_rain(microp_uniform, qric, qcic, ncic, relvar, accre_enhan, pra, npra, mgncol*nlev)
        2109             :   endif
        2110             :
        2111           0 :   pra = pra*micro_mg_accre_enhan_fact
        2112           0 :   npra = npra*micro_mg_accre_enhan_fact
        2113             :
        2114           0 :   call self_collection_rain(rho, qric, nric, nragg, mgncol*nlev)
        2115             :
        2116           0 :   if (do_cldice) then
        2117           0 :      call accrete_cloud_ice_snow(t, rho, asn, qiic, niic, qsic, lams, n0s, prai, nprai, mgncol*nlev)
        2118             :   else
        2119             :      !$acc parallel vector_length(VLENS) default(present)
        2120             :      !$acc loop gang vector collapse(2)
        2121           0 :      do k=1,nlev
        2122           0 :         do i=1,mgncol
        2123           0 :            prai(i,k) = 0._r8
        2124           0 :            nprai(i,k) = 0._r8
        2125             :         end do
        2126             :      end do
        2127             :      !$acc end parallel
        2128             :   end if
        2129             :
        2130           0 :   call bergeron_process_snow(t, rho, dv, mu, sc, qvl, qvi, asn, qcic, qsic, lams, n0s, bergs, mgncol*nlev)
        2131             :   !$acc parallel vector_length(VLENS) default(present)
        2132             :   !$acc loop gang vector collapse(2)
        2133           0 :   do k=1,nlev
        2134           0 :      do i=1,mgncol
        2135           0 :         bergs(i,k)=bergs(i,k)*micro_mg_berg_eff_factor
        2136             :      end do
        2137             :   end do
        2138             :   !$acc end parallel
        2139             :
        2140           0 :   if (do_cldice) then
        2141             :      call ice_deposition_sublimation(t, q, qi, ni, icldm, rho, dv, qvl, qvi, &
        2142           0 :                                      berg, vap_dep, ice_sublim, mgncol*nlev)
        2143             :      !$acc parallel vector_length(VLENS) default(present)
        2144             :      !$acc loop gang vector collapse(2)
        2145           0 :      do k=1,nlev
        2146           0 :         do i=1,mgncol
        2147           0 :            berg(i,k)=berg(i,k)*micro_mg_berg_eff_factor
        2148           0 :            if (ice_sublim(i,k) < 0._r8 .and. qi(i,k) > qsmall .and. icldm(i,k) > mincld) then
        2149           0 :               nsubi(i,k) = sublim_factor*ice_sublim(i,k) / qi(i,k) * ni(i,k) / icldm(i,k)
        2150             :            else
        2151           0 :               nsubi(i,k) = 0._r8
        2152             :            end if
        2153             :
        2154             :            ! bergeron process should not reduce nc unless
        2155             :            ! all ql is removed (which is handled elsewhere)
        2156             :            !in fact, nothing in this entire file makes nsubc nonzero.
        2157           0 :            nsubc(i,k) = 0._r8
        2158             :
        2159             :         end do
        2160             :      end do
        2161             :      !$acc end parallel
        2162             :   end if !do_cldice
        2163             :
        2164             : ! Process rate calls for graupel
        2165             : !===================================================================
        2166             :
        2167           0 :   if (do_hail.or.do_graupel) then
        2168           0 :      call graupel_collecting_snow(qsic, qric, umr, ums, rho, lamr, n0r, lams, n0s, psacr, mgncol*nlev)
        2169             :
        2170           0 :      call graupel_collecting_cld_water(qgic, qcic, ncic, rho, n0g, lamg, bgtmp, agn, psacwg, npsacwg, mgncol*nlev)
        2171             :
        2172           0 :      psacwg = psacwg*micro_mg_iaccr_factor
        2173           0 :      npsacwg = npsacwg*micro_mg_iaccr_factor
        2174             :
        2175             :      call graupel_riming_liquid_snow(psacws, qsic, qcic, nsic, rho, rhosn, rhogtmp, asn, &
        2176           0 :                                      lams, n0s, deltat, pgsacw, nscng, mgncol*nlev)
        2177             :
        2178             :      call graupel_collecting_rain(qric, qgic, umg, umr, ung, unr, rho, n0r, &
        2179           0 :                                   lamr, n0g, lamg, pracg, npracg, mgncol*nlev)
        2180             :
        2181           0 :      pracg = pracg*micro_mg_iaccr_factor
        2182           0 :      npracg = npracg*micro_mg_iaccr_factor
        2183             :
        2184             : !AG note: Graupel rain riming snow changes
        2185             : !    pracs, npracs, (accretion of rain by snow)  psacr (collection of snow by rain)
        2186             :
        2187             :      call graupel_rain_riming_snow(pracs, npracs, psacr, qsic, qric, nric, nsic, &
        2188           0 :                                    n0s, lams, n0r, lamr, deltat, pgracs, ngracs, mgncol*nlev)
        2189             :
        2190           0 :      call graupel_rime_splintering(t, qcic, qric, qgic, psacwg, pracg, qmultg, nmultg, qmultrg, nmultrg,mgncol*nlev)
        2191             :
        2192             :
        2193             :      call evaporate_sublimate_precip_graupel(t, rho, dv, mu, sc, q, qvl, qvi, lcldm, precip_frac, arn, asn, agn, &
        2194             :                                              bgtmp, qcic, qiic, qric, qsic, qgic, lamr, n0r, lams, n0s, lamg, n0g, &
        2195           0 :                                              pre, prds, prdg, am_evp_st, mgncol*nlev, evap_rhthrsh_ifs)
        2196             :   else
        2197             :      ! Routine without Graupel (original)
        2198             :      call evaporate_sublimate_precip(t, rho, dv, mu, sc, q, qvl, qvi, lcldm, precip_frac, arn, asn, qcic, qiic, &
        2199           0 :                                      qric, qsic, lamr, n0r, lams, n0s, pre, prds, am_evp_st, mgncol*nlev, evap_rhthrsh_ifs)
        2200             :   end if ! end do_graupel/hail loop
        2201             :
        2202             : ! scale precip evaporation to match IFS 'new' version (option 2)
        2203           0 :   if (evap_scl_ifs) then
        2204             :      !$acc parallel vector_length(VLENS) default(present)
        2205             :      !$acc loop gang vector collapse(2)
        2206           0 :      do k=1,nlev
        2207           0 :         do i=1,mgncol
        2208           0 :            pre(i,k)= 0.15_r8 * pre(i,k)
        2209             :         end do
        2210             :      end do
        2211             :      !$acc end parallel
        2212             :   end if
        2213             :
        2214             :   !$acc parallel vector_length(VLENS) default(present)
        2215             :   !$acc loop gang vector collapse(2) private(dum,ratio)
        2216           0 :   do k=1,nlev
        2217           0 :      do i=1,mgncol
        2218             :         ! conservation to ensure no negative values of cloud water/precipitation
        2219             :         ! in case microphysical process rates are large
        2220             :         !===================================================================
        2221             :
        2222             :         ! note: for check on conservation, processes are multiplied by omsm
        2223             :         ! to prevent problems due to round off error
        2224             :
        2225             :         ! conservation of qc
        2226             :         !-------------------------------------------------------------------
        2227           0 :         dum = ((prc(i,k)+pra(i,k)+mnuccc(i,k)+mnucct(i,k)+msacwi(i,k)+ &
        2228             :              psacws(i,k)+bergs(i,k)+qmultg(i,k)+psacwg(i,k)+pgsacw(i,k))*lcldm(i,k)+ &
        2229           0 :              berg(i,k))*deltat
        2230           0 :         if (dum.gt.qc(i,k)) then
        2231             :            ratio = qc(i,k)*rdeltat/((prc(i,k)+pra(i,k)+mnuccc(i,k)+mnucct(i,k)+ &
        2232             :                 msacwi(i,k)+psacws(i,k)+bergs(i,k)+qmultg(i,k)+psacwg(i,k)+pgsacw(i,k))*lcldm(i,k)+&
        2233           0 :                 berg(i,k))*omsm
        2234           0 :            qmultg(i,k)=qmultg(i,k)*ratio
        2235           0 :            psacwg(i,k)=psacwg(i,k)*ratio
        2236           0 :            pgsacw(i,k)=pgsacw(i,k)*ratio
        2237           0 :            prc(i,k) = prc(i,k)*ratio
        2238           0 :            pra(i,k) = pra(i,k)*ratio
        2239           0 :            mnuccc(i,k) = mnuccc(i,k)*ratio
        2240           0 :            mnucct(i,k) = mnucct(i,k)*ratio
        2241           0 :            msacwi(i,k) = msacwi(i,k)*ratio
        2242           0 :            psacws(i,k) = psacws(i,k)*ratio
        2243           0 :            bergs(i,k) = bergs(i,k)*ratio
        2244           0 :            berg(i,k) = berg(i,k)*ratio
        2245           0 :            qcrat(i,k) = ratio
        2246             :         else
        2247           0 :            qcrat(i,k) = 1._r8
        2248             :         end if
        2249             :         !PMC 12/3/12: ratio is also frac of step w/ liquid.
        2250             :         !thus we apply berg for "ratio" of timestep and vapor
        2251             :         !deposition for the remaining frac of the timestep.
        2252           0 :         if (qc(i,k) >= qsmall) then
        2253           0 :            vap_dep(i,k) = vap_dep(i,k)*(1._r8-qcrat(i,k))
        2254             :         end if
        2255             :      end do
        2256             :   end do
        2257             :
        2258             :   !$acc loop gang vector collapse(2) private(dum,dum1)
        2259           0 :   do k=1,nlev
        2260           0 :      do i=1,mgncol
        2261             :         !=================================================================
        2262             :         ! apply limiter to ensure that ice/snow sublimation and rain evap
        2263             :         ! don't push conditions into supersaturation, and ice deposition/nucleation don't
        2264             :         ! push conditions into sub-saturation
        2265             :         ! note this is done after qc conservation since we don't know how large
        2266             :         ! vap_dep is before then
        2267             :         ! estimates are only approximate since other process terms haven't been limited
        2268             :         ! for conservation yet
        2269             :
        2270             :         ! first limit ice deposition/nucleation vap_dep + mnuccd
        2271           0 :         dum1 = vap_dep(i,k) + mnuccd(i,k)
        2272           0 :         if (dum1 > 1.e-20_r8) then
        2273           0 :            dum = (q(i,k)-qvi(i,k))/(1._r8 + xxls_squared*qvi(i,k)/(cpp*rv*t(i,k)**2))*rdeltat
        2274           0 :            dum = max(dum,0._r8)
        2275           0 :            if (dum1 > dum) then
        2276             :               ! Allocate the limited "dum" tendency to mnuccd and vap_dep
        2277             :               ! processes. Don't divide by cloud fraction; these are grid-
        2278             :               ! mean rates.
        2279           0 :               dum1 = mnuccd(i,k) / (vap_dep(i,k)+mnuccd(i,k))
        2280           0 :               mnuccd(i,k) = dum*dum1
        2281           0 :               vap_dep(i,k) = dum - mnuccd(i,k)
        2282             :            end if
        2283             :         end if
        2284             :      end do
        2285             :   end do
        2286             :
        2287             :   !$acc loop gang vector collapse(2) private(dum,ratio)
        2288           0 :   do k=1,nlev
        2289           0 :      do i=1,mgncol
        2290             :         !===================================================================
        2291             :         ! conservation of nc
        2292             :         !-------------------------------------------------------------------
        2293           0 :         dum = (nprc1(i,k)+npra(i,k)+nnuccc(i,k)+nnucct(i,k)+ &
        2294           0 :                npsacws(i,k)-nsubc(i,k)+npsacwg(i,k))*lcldm(i,k)*deltat
        2295           0 :         if (dum.gt.nc(i,k)) then
        2296             :            ratio = nc(i,k)*rdeltat/((nprc1(i,k)+npra(i,k)+nnuccc(i,k)+nnucct(i,k)+&
        2297           0 :                    npsacws(i,k)-nsubc(i,k)+npsacwg(i,k))*lcldm(i,k))*omsm
        2298           0 :            npsacwg(i,k) = npsacwg(i,k)*ratio
        2299           0 :            nprc1(i,k)   = nprc1(i,k)*ratio
        2300           0 :            npra(i,k)    = npra(i,k)*ratio
        2301           0 :            nnuccc(i,k)  = nnuccc(i,k)*ratio
        2302           0 :            nnucct(i,k)  = nnucct(i,k)*ratio
        2303           0 :            npsacws(i,k) = npsacws(i,k)*ratio
        2304           0 :            nsubc(i,k)   = nsubc(i,k)*ratio
        2305             :         end if
        2306           0 :         mnuccri(i,k)=0._r8
        2307           0 :         nnuccri(i,k)=0._r8
        2308             :
        2309           0 :         if (do_cldice) then
        2310             :            ! freezing of rain to produce ice if mean rain size is smaller than Dcs
        2311           0 :            if (lamr(i,k) > qsmall) then
        2312           0 :               if (1._r8/lamr(i,k) < Dcs) then
        2313           0 :                  mnuccri(i,k)=mnuccr(i,k)
        2314           0 :                  nnuccri(i,k)=nnuccr(i,k)
        2315           0 :                  mnuccr(i,k)=0._r8
        2316           0 :                  nnuccr(i,k)=0._r8
        2317             :               end if
        2318             :            end if
        2319             :         end if
        2320             :      end do
        2321             :   end do
        2322             :
        2323             :   !$acc loop gang vector collapse(2) private(dum,ratio)
        2324           0 :   do k=1,nlev
        2325           0 :      do i=1,mgncol
        2326             :         ! conservation of rain mixing ratio
        2327             :         !-------------------------------------------------------------------
        2328           0 :         dum = ((-pre(i,k)+pracs(i,k)+mnuccr(i,k)+mnuccri(i,k) &
        2329             :              +qmultrg(i,k)+pracg(i,k)+pgracs(i,k))*precip_frac(i,k)- &
        2330           0 :              (pra(i,k)+prc(i,k))*lcldm(i,k))*deltat
        2331             :         ! note that qrtend is included below because of instantaneous freezing/melt
        2332           0 :         if (dum.gt.qr(i,k).and. &
        2333           0 :              (-pre(i,k)+pracs(i,k)+mnuccr(i,k)+mnuccri(i,k)+qmultrg(i,k)+pracg(i,k)+pgracs(i,k)).ge.qsmall) then
        2334             :            ratio = (qr(i,k)*rdeltat+(pra(i,k)+prc(i,k))*lcldm(i,k))/   &
        2335             :                 precip_frac(i,k)/(-pre(i,k)+pracs(i,k)+mnuccr(i,k)+mnuccri(i,k) &
        2336           0 :                 +qmultrg(i,k)+pracg(i,k)+pgracs(i,k))*omsm
        2337           0 :            qmultrg(i,k)= qmultrg(i,k)*ratio
        2338           0 :            pracg(i,k)=pracg(i,k)*ratio
        2339           0 :            pgracs(i,k)=pgracs(i,k)*ratio
        2340           0 :            pre(i,k)=pre(i,k)*ratio
        2341           0 :            pracs(i,k)=pracs(i,k)*ratio
        2342           0 :            mnuccr(i,k)=mnuccr(i,k)*ratio
        2343           0 :            mnuccri(i,k)=mnuccri(i,k)*ratio
        2344             :         end if
        2345             :      end do
        2346             :   end do
        2347             :
        2348             :   !$acc loop gang vector collapse(2)
        2349           0 :   do k=1,nlev
        2350           0 :      do i=1,mgncol
        2351             :         ! conservation of rain number
        2352             :         !-------------------------------------------------------------------
        2353             :         ! Add evaporation of rain number.
        2354           0 :         if (pre(i,k) < 0._r8) then
        2355           0 :            nsubr(i,k) = pre(i,k)*nr(i,k)/qr(i,k)
        2356             :         else
        2357           0 :            nsubr(i,k) = 0._r8
        2358             :         end if
        2359             :      end do
        2360             :   end do
        2361             :
        2362             :   !$acc loop gang vector collapse(2) private(dum,ratio)
        2363           0 :   do k=1,nlev
        2364           0 :      do i=1,mgncol
        2365           0 :         dum = ((-nsubr(i,k)+npracs(i,k)+nnuccr(i,k)+nnuccri(i,k)-nragg(i,k)+npracg(i,k)+ngracs(i,k)) &
        2366           0 :              *precip_frac(i,k)- nprc(i,k)*lcldm(i,k))*deltat
        2367           0 :         if (dum.gt.nr(i,k)) then
        2368             :            ratio = (nr(i,k)*rdeltat+nprc(i,k)*lcldm(i,k))/precip_frac(i,k)/ &
        2369           0 :                 (-nsubr(i,k)+npracs(i,k)+nnuccr(i,k)+nnuccri(i,k)-nragg(i,k)+npracg(i,k)+ngracs(i,k))*omsm
        2370           0 :            npracg(i,k)=npracg(i,k)*ratio
        2371           0 :            ngracs(i,k)=ngracs(i,k)*ratio
        2372           0 :            nragg(i,k)=nragg(i,k)*ratio
        2373           0 :            npracs(i,k)=npracs(i,k)*ratio
        2374           0 :            nnuccr(i,k)=nnuccr(i,k)*ratio
        2375           0 :            nsubr(i,k)=nsubr(i,k)*ratio
        2376           0 :            nnuccri(i,k)=nnuccri(i,k)*ratio
        2377             :         end if
        2378             :      end do
        2379             :   end do
        2380             :
        2381           0 :   if (do_cldice) then
        2382             :      !$acc loop gang vector collapse(2) private(dum,ratio)
        2383           0 :      do k=1,nlev
        2384           0 :         do i=1,mgncol
        2385             :            ! conservation of qi
        2386             :            !-------------------------------------------------------------------
        2387           0 :            dum = ((-mnuccc(i,k)-mnucct(i,k)-mnudep(i,k)-msacwi(i,k)-qmultg(i,k))*lcldm(i,k)+(prci(i,k)+ &
        2388             :                 prai(i,k))*icldm(i,k)+(-qmultrg(i,k)-mnuccri(i,k))*precip_frac(i,k) &
        2389           0 :                 -ice_sublim(i,k)-vap_dep(i,k)-berg(i,k)-mnuccd(i,k))*deltat
        2390           0 :            if (dum.gt.qi(i,k)) then
        2391             :               ratio = (qi(i,k)*rdeltat+vap_dep(i,k)+berg(i,k)+mnuccd(i,k)+ &
        2392             :                    (mnuccc(i,k)+mnucct(i,k)+mnudep(i,k)+msacwi(i,k)+qmultg(i,k))*lcldm(i,k)+ &
        2393             :                    (qmultrg(i,k)+mnuccri(i,k))*precip_frac(i,k))/ &
        2394           0 :                    ((prci(i,k)+prai(i,k))*icldm(i,k)-ice_sublim(i,k))*omsm
        2395           0 :               prci(i,k) = prci(i,k)*ratio
        2396           0 :               prai(i,k) = prai(i,k)*ratio
        2397           0 :               ice_sublim(i,k) = ice_sublim(i,k)*ratio
        2398             :            end if
        2399             :         end do
        2400             :      end do
        2401             :   end if
        2402             :
        2403           0 :   if (do_cldice) then
        2404             :      !$acc loop gang vector collapse(2) private(dum,ratio,tmpfrz)
        2405           0 :      do k=1,nlev
        2406           0 :         do i=1,mgncol
        2407             :            ! conservation of ni
        2408             :            !-------------------------------------------------------------------
        2409           0 :            if (use_hetfrz_classnuc) then
        2410           0 :               tmpfrz = nnuccc(i,k)
        2411             :            else
        2412             :               tmpfrz = 0._r8
        2413             :            end if
        2414           0 :            dum = ((-nnucct(i,k)-tmpfrz-nnudep(i,k)-nsacwi(i,k)-nmultg(i,k))*lcldm(i,k)+(nprci(i,k)+ &
        2415             :                 nprai(i,k)-nsubi(i,k))*icldm(i,k)+(-nmultrg(i,k)-nnuccri(i,k))*precip_frac(i,k)- &
        2416           0 :                 nnuccd(i,k))*deltat
        2417           0 :            if (dum.gt.ni(i,k)) then
        2418             :               ratio = (ni(i,k)*rdeltat+nnuccd(i,k)+ &
        2419             :                  (nnucct(i,k)+tmpfrz+nnudep(i,k)+nsacwi(i,k)+nmultg(i,k))*lcldm(i,k)+ &
        2420             :                  (nnuccri(i,k)+nmultrg(i,k))*precip_frac(i,k))/ &
        2421           0 :                  ((nprci(i,k)+nprai(i,k)-nsubi(i,k))*icldm(i,k))*omsm
        2422           0 :               nprci(i,k) = nprci(i,k)*ratio
        2423           0 :               nprai(i,k) = nprai(i,k)*ratio
        2424           0 :               nsubi(i,k) = nsubi(i,k)*ratio
        2425             :            end if
        2426             :         end do
        2427             :      end do
        2428             :   end if
        2429             :
        2430             :   !$acc loop gang vector collapse(2) private(dum,ratio)
        2431           0 :   do k=1,nlev
        2432           0 :      do i=1,mgncol
        2433             :         ! conservation of snow mixing ratio
        2434             :         !-------------------------------------------------------------------
        2435           0 :         if (do_hail .or. do_graupel) then
        2436             :         ! NOTE: mnuccr is moved to graupel when active
        2437             :         ! psacr is a positive value, but a loss for snow
        2438             :         !HM: psacr is positive in dum (two negatives)
        2439           0 :            dum = (-(prds(i,k)+pracs(i,k)-psacr(i,k))*precip_frac(i,k)-(prai(i,k)+prci(i,k))*icldm(i,k) &
        2440           0 :              -(bergs(i,k)+psacws(i,k))*lcldm(i,k))*deltat
        2441             :         else
        2442           0 :            dum = (-(prds(i,k)+pracs(i,k)+mnuccr(i,k))*precip_frac(i,k)-(prai(i,k)+prci(i,k))*icldm(i,k) &
        2443           0 :              -(bergs(i,k)+psacws(i,k))*lcldm(i,k))*deltat
        2444             :         end if
        2445           0 :         if (dum.gt.qs(i,k).and.(psacr(i,k)-prds(i,k)).ge.qsmall) then
        2446           0 :            if (do_hail .or. do_graupel) then
        2447             :               ratio = (qs(i,k)*rdeltat+(prai(i,k)+prci(i,k))*icldm(i,k)+ &
        2448             :                    (bergs(i,k)+psacws(i,k))*lcldm(i,k)+pracs(i,k)*precip_frac(i,k))/ &
        2449           0 :                    precip_frac(i,k)/(psacr(i,k)-prds(i,k))*omsm
        2450           0 :               psacr(i,k)=psacr(i,k)*ratio
        2451             :            else
        2452             :               ratio = (qs(i,k)*rdeltat+(prai(i,k)+prci(i,k))*icldm(i,k)+ &
        2453             :                    (bergs(i,k)+psacws(i,k))*lcldm(i,k)+(pracs(i,k)+mnuccr(i,k))*precip_frac(i,k))/ &
        2454           0 :                    precip_frac(i,k)/(-prds(i,k))*omsm
        2455             :            end if
        2456           0 :            prds(i,k)=prds(i,k)*ratio
        2457             :         end if
        2458             :      end do
        2459             :   end do
        2460             :
        2461             :   !$acc loop gang vector collapse(2) private(dum,ratio)
        2462           0 :   do k=1,nlev
        2463           0 :      do i=1,mgncol
        2464             :         ! conservation of snow number
        2465             :         !-------------------------------------------------------------------
        2466             :         ! calculate loss of number due to sublimation
        2467             :         ! for now neglect sublimation of ns
        2468           0 :         nsubs(i,k)=0._r8
        2469           0 :         if (do_hail .or. do_graupel) then
        2470           0 :            dum = ((-nsagg(i,k)-nsubs(i,k)+ngracs(i,k))*precip_frac(i,k)-nprci(i,k)*icldm(i,k)+nscng(i,k)*lcldm(i,k))*deltat
        2471             :         else
        2472           0 :            dum = ((-nsagg(i,k)-nsubs(i,k)-nnuccr(i,k))*precip_frac(i,k)-nprci(i,k)*icldm(i,k))*deltat
        2473             :         end if
        2474           0 :         if (dum.gt.ns(i,k)) then
        2475           0 :            if (do_hail .or. do_graupel) then
        2476             :               ratio = (ns(i,k)*rdeltat+nprci(i,k)*icldm(i,k))/precip_frac(i,k)/ &
        2477           0 :                    (-nsubs(i,k)-nsagg(i,k)+ngracs(i,k)+lcldm(i,k)/precip_frac(i,k)*nscng(i,k))*omsm
        2478           0 :               nscng(i,k)=nscng(i,k)*ratio
        2479           0 :               ngracs(i,k)=ngracs(i,k)*ratio
        2480             :            else
        2481             :               ratio = (ns(i,k)*rdeltat+nnuccr(i,k)* &
        2482             :                    precip_frac(i,k)+nprci(i,k)*icldm(i,k))/precip_frac(i,k)/ &
        2483           0 :                    (-nsubs(i,k)-nsagg(i,k))*omsm
        2484             :            endif
        2485           0 :            nsubs(i,k)=nsubs(i,k)*ratio
        2486           0 :            nsagg(i,k)=nsagg(i,k)*ratio
        2487             :         end if
        2488             :      end do
        2489             :   end do
        2490             :
        2491             : ! Graupel Conservation Checks
        2492             : !-------------------------------------------------------------------
        2493             :
        2494           0 :   if (do_hail.or.do_graupel) then
        2495             :      ! conservation of graupel mass
        2496             :      !-------------------------------------------------------------------
        2497             :
        2498             :      !$acc loop gang vector collapse(2) private(dum,ratio)
        2499           0 :      do k=1,nlev
        2500           0 :         do i=1,mgncol
        2501           0 :            dum= ((-pracg(i,k)-pgracs(i,k)-prdg(i,k)-psacr(i,k)-mnuccr(i,k))*precip_frac(i,k) &
        2502           0 :                 + (-psacwg(i,k)-pgsacw(i,k))*lcldm(i,k))*deltat
        2503           0 :            if (dum.gt.qg(i,k)) then
        2504             :               ! note: prdg is always negative (like prds), so it needs to be subtracted in ratio
        2505             :               ratio = (qg(i,k)*rdeltat + (pracg(i,k)+pgracs(i,k)+psacr(i,k)+mnuccr(i,k))*precip_frac(i,k) &
        2506           0 :                        + (psacwg(i,k)+pgsacw(i,k))*lcldm(i,k)) / ((-prdg(i,k))*precip_frac(i,k)) * omsm
        2507           0 :               prdg(i,k)= prdg(i,k)*ratio
        2508             :            end if
        2509             :         end do
        2510             :      end do
        2511             :      ! conservation of graupel number: not needed, no sinks
        2512             :      !-------------------------------------------------------------------
        2513             :   end if
        2514             :
        2515             :   !$acc loop gang vector collapse(2)
        2516           0 :   do k=1,nlev
        2517           0 :      do i=1,mgncol
        2518             :         ! next limit ice and snow sublimation and rain evaporation
        2519             :         ! get estimate of q and t at end of time step
        2520             :         ! don't include other microphysical processes since they haven't
        2521             :         ! been limited via conservation checks yet
        2522           0 :         qtmpAI(i,k)=q(i,k)-(ice_sublim(i,k)+vap_dep(i,k)+mnuccd(i,k)+ &
        2523           0 :                 (pre(i,k)+prds(i,k)+prdg(i,k))*precip_frac(i,k))*deltat
        2524             :         ttmpA(i,k)=t(i,k)+((pre(i,k)*precip_frac(i,k))*xxlv+ &
        2525           0 :              ((prds(i,k)+prdg(i,k))*precip_frac(i,k)+vap_dep(i,k)+ice_sublim(i,k)+mnuccd(i,k))*xxls)*deltat/cpp
        2526             :      end do
        2527             :   end do
        2528             :   !$acc end parallel
        2529             :
        2530             :   ! use rhw to allow ice supersaturation
        2531           0 :   call qsat_water(ttmpA, p, esnA, qvnAI, mgncol*nlev)
        2532             :
        2533             :   !$acc parallel vector_length(VLENS) default(present)
        2534             :   !$acc loop gang vector collapse(2)
        2535           0 :   do k=1,nlev
        2536           0 :      do i=1,mgncol
        2537           0 :         if ((pre(i,k)+prds(i,k)+prdg(i,k))*precip_frac(i,k)+ice_sublim(i,k) < -1.e-20_r8) then
        2538             :            ! modify ice/precip evaporation rate if q > qsat
        2539           0 :            if (qtmpAI(i,k) > qvnAI(i,k)) then
        2540           0 :               dum1A(i,k)=pre(i,k)*precip_frac(i,k)/((pre(i,k)+prds(i,k)+prdg(i,k))*precip_frac(i,k)+ice_sublim(i,k))
        2541           0 :               dum2A(i,k)=prds(i,k)*precip_frac(i,k)/((pre(i,k)+prds(i,k)+prdg(i,k))*precip_frac(i,k)+ice_sublim(i,k))
        2542           0 :               dum3A(i,k)=prdg(i,k)*precip_frac(i,k)/((pre(i,k)+prds(i,k)+prdg(i,k))*precip_frac(i,k)+ice_sublim(i,k))
        2543             :               ! recalculate q and t after vap_dep and mnuccd but without evap or sublim
        2544           0 :               ttmpA(i,k)=t(i,k)+((vap_dep(i,k)+mnuccd(i,k))*xxls)*deltat/cpp
        2545           0 :               dum_2D(i,k)=q(i,k)-(vap_dep(i,k)+mnuccd(i,k))*deltat
        2546             :            end if
        2547             :         end if
        2548             :      end do
        2549             :   end do
        2550             :   !$acc end parallel
        2551             :
        2552             :   ! use rhw to allow ice supersaturation
        2553           0 :   call qsat_water(ttmpA, p, esnA, qvnA, mgncol*nlev)
        2554             :
        2555             :   !$acc parallel vector_length(VLENS) default(present)
        2556             :   !$acc loop gang vector collapse(2) private(dum)
        2557           0 :   do k=1,nlev
        2558           0 :      do i=1,mgncol
        2559           0 :         if ((pre(i,k)+prds(i,k)+prdg(i,k))*precip_frac(i,k)+ice_sublim(i,k) < -1.e-20_r8) then
        2560             :            ! modify ice/precip evaporation rate if q > qsat
        2561           0 :            if (qtmpAI(i,k) > qvnAI(i,k)) then
        2562           0 :               dum=(dum_2D(i,k)-qvnA(i,k))/(1._r8 + xxlv_squared*qvnA(i,k)/(cpp*rv*ttmpA(i,k)**2))
        2563           0 :               dum=min(dum,0._r8)
        2564             :               ! modify rates if needed, divide by precip_frac to get local (in-precip) value
        2565           0 :               pre(i,k)=dum*dum1A(i,k)*rdeltat/precip_frac(i,k)
        2566             :            end if
        2567             :         end if
        2568             :      end do
        2569             :   end do
        2570             :   !$acc end parallel
        2571             :
        2572             :   ! do separately using RHI for prds and ice_sublim
        2573           0 :   call qsat_ice(ttmpA, p, esnA, qvnA, mgncol*nlev)
        2574             :
        2575             :   !$acc parallel vector_length(VLENS) default(present)
        2576             :   !$acc loop gang vector collapse(2) private(dum)
        2577           0 :   do k=1,nlev
        2578           0 :      do i=1,mgncol
        2579           0 :         if ((pre(i,k)+prds(i,k)+prdg(i,k))*precip_frac(i,k)+ice_sublim(i,k) < -1.e-20_r8) then
        2580             :            ! modify ice/precip evaporation rate if q > qsat
        2581           0 :            if (qtmpAI(i,k) > qvnAI(i,k)) then
        2582           0 :               dum=(dum_2D(i,k)-qvnA(i,k))/(1._r8 + xxls_squared*qvnA(i,k)/(cpp*rv*ttmpA(i,k)**2))
        2583           0 :               dum=min(dum,0._r8)
        2584             :               ! modify rates if needed, divide by precip_frac to get local (in-precip) value
        2585           0 :               prds(i,k) = dum*dum2A(i,k)*rdeltat/precip_frac(i,k)
        2586           0 :               prdg(i,k) = dum*dum3A(i,k)*rdeltat/precip_frac(i,k)
        2587             :               ! don't divide ice_sublim by cloud fraction since it is grid-averaged
        2588           0 :               dum1A(i,k) = (1._r8-dum1A(i,k)-dum2A(i,k)-dum3A(i,k))
        2589           0 :               ice_sublim(i,k) = dum*dum1A(i,k)*rdeltat
        2590             :            end if
        2591             :         end if
        2592             :      end do
        2593             :   end do
        2594             :
        2595             :   ! Big "administration" loop enforces conservation, updates variables
        2596             :   ! that accumulate over substeps, and sets output variables.
        2597             :
        2598             :   !$acc loop gang vector collapse(2)
        2599           0 :   do k=1,nlev
        2600           0 :      do i=1,mgncol
        2601             :         ! get tendencies due to microphysical conversion processes
        2602             :         !==========================================================
        2603             :         ! note: tendencies are multiplied by appropriate cloud/precip
        2604             :         ! fraction to get grid-scale values
        2605             :         ! note: vap_dep is already grid-average values
        2606             :
        2607             :         ! The net tendencies need to be added to rather than overwritten,
        2608             :         ! because they may have a value already set for instantaneous
        2609             :         ! melting/freezing.
        2610           0 :         qvlat(i,k) = qvlat(i,k)-(pre(i,k)+prds(i,k))*precip_frac(i,k)-&
        2611             :              vap_dep(i,k)-ice_sublim(i,k)-mnuccd(i,k)-mnudep(i,k)*lcldm(i,k) &
        2612           0 :              -prdg(i,k)*precip_frac(i,k)
        2613             :         tlat(i,k) = tlat(i,k)+((pre(i,k)*precip_frac(i,k))*xxlv+ &
        2614             :              ((prds(i,k)+prdg(i,k))*precip_frac(i,k)+vap_dep(i,k)+ice_sublim(i,k)+ &
        2615             :                  mnuccd(i,k)+mnudep(i,k)*lcldm(i,k))*xxls+ &
        2616             :              ((bergs(i,k)+psacws(i,k)+mnuccc(i,k)+mnucct(i,k)+msacwi(i,k)+psacwg(i,k)+ &
        2617             :                   qmultg(i,k)+pgsacw(i,k))*lcldm(i,k)+ &
        2618             :              (mnuccr(i,k)+pracs(i,k)+mnuccri(i,k)+pracg(i,k)+pgracs(i,k)+qmultrg(i,k))*precip_frac(i,k)+ &
        2619           0 :                   berg(i,k))*xlf)
        2620             :         qctend(i,k) = qctend(i,k)+ &
        2621             :              (-pra(i,k)-prc(i,k)-mnuccc(i,k)-mnucct(i,k)-msacwi(i,k)- &
        2622           0 :              psacws(i,k)-bergs(i,k)-qmultg(i,k)-psacwg(i,k)-pgsacw(i,k))*lcldm(i,k)-berg(i,k)
        2623             :
        2624           0 :         if (do_cldice) then
        2625             :            qitend(i,k) = qitend(i,k)+ &
        2626             :               (mnuccc(i,k)+mnucct(i,k)+mnudep(i,k)+msacwi(i,k)+qmultg(i,k))*lcldm(i,k)+(-prci(i,k)- &
        2627             :               prai(i,k))*icldm(i,k)+vap_dep(i,k)+berg(i,k)+ice_sublim(i,k)+ &
        2628           0 :               mnuccd(i,k)+(mnuccri(i,k)+qmultrg(i,k))*precip_frac(i,k)
        2629             :         end if
        2630             :
        2631             :         qrtend(i,k) = qrtend(i,k)+ &
        2632             :              (pra(i,k)+prc(i,k))*lcldm(i,k)+(pre(i,k)-pracs(i,k)- &
        2633           0 :              mnuccr(i,k)-mnuccri(i,k)-qmultrg(i,k)-pracg(i,k)-pgracs(i,k))*precip_frac(i,k)
        2634             :
        2635           0 :         if (do_hail.or.do_graupel) then
        2636             :            qgtend(i,k) = qgtend(i,k) + (pracg(i,k)+pgracs(i,k)+prdg(i,k)+psacr(i,k)+mnuccr(i,k))*precip_frac(i,k) &
        2637           0 :                 + (psacwg(i,k)+pgsacw(i,k))*lcldm(i,k)
        2638             :            qstend(i,k) = qstend(i,k)+ &
        2639             :                 (prai(i,k)+prci(i,k))*icldm(i,k)+(psacws(i,k)+bergs(i,k))*lcldm(i,k)+(prds(i,k)+ &
        2640           0 :                 pracs(i,k)-psacr(i,k))*precip_frac(i,k)
        2641             :         else
        2642             :            !necessary since mnuccr moved to graupel
        2643             :            qstend(i,k) = qstend(i,k)+ &
        2644             :                 (prai(i,k)+prci(i,k))*icldm(i,k)+(psacws(i,k)+bergs(i,k))*lcldm(i,k)+(prds(i,k)+ &
        2645           0 :                 pracs(i,k)+mnuccr(i,k))*precip_frac(i,k)
        2646             :         end if
        2647             :      end do
        2648             :   end do
        2649             :
        2650             :   !$acc loop gang vector collapse(2)
        2651           0 :   do k=1,nlev
        2652           0 :      do i=1,mgncol
        2653           0 :         cmeout(i,k) = vap_dep(i,k) + ice_sublim(i,k) + mnuccd(i,k)
        2654             :         ! add output for cmei (accumulate)
        2655           0 :         cmeitot(i,k) = vap_dep(i,k) + ice_sublim(i,k) + mnuccd(i,k)
        2656             :         !-------------------------------------------------------------------
        2657             :         ! evaporation/sublimation is stored here as positive term
        2658             :         ! Add to evapsnow via prdg
        2659           0 :         evapsnow(i,k) = (-prds(i,k)-prdg(i,k))*precip_frac(i,k)
        2660           0 :         nevapr(i,k) = -pre(i,k)*precip_frac(i,k)
        2661           0 :         prer_evap(i,k) = -pre(i,k)*precip_frac(i,k)
        2662             :         ! change to make sure prain is positive: do not remove snow from
        2663             :         ! prain used for wet deposition
        2664             :         prain(i,k) = (pra(i,k)+prc(i,k))*lcldm(i,k)+(-pracs(i,k)- &
        2665           0 :              mnuccr(i,k)-mnuccri(i,k))*precip_frac(i,k)
        2666           0 :         if (do_hail .or. do_graupel) then
        2667             :            prodsnow(i,k) = (prai(i,k)+prci(i,k))*icldm(i,k)+(psacws(i,k)+bergs(i,k))*lcldm(i,k)+(&
        2668           0 :                 pracs(i,k))*precip_frac(i,k)
        2669             :         else
        2670             :            prodsnow(i,k) = (prai(i,k)+prci(i,k))*icldm(i,k)+(psacws(i,k)+bergs(i,k))*lcldm(i,k)+(&
        2671           0 :                 pracs(i,k)+mnuccr(i,k))*precip_frac(i,k)
        2672             :         end if
        2673             :         ! following are used to calculate 1st order conversion rate of cloud water
        2674             :         !    to rain and snow (1/s), for later use in aerosol wet removal routine
        2675             :         ! previously, wetdepa used (prain/qc) for this, and the qc in wetdepa may be smaller than the qc
        2676             :         !    used to calculate pra, prc, ... in this routine
        2677             :         ! qcsinksum_rate1ord = { rate of direct transfer of cloud water to rain & snow }
        2678             :         !                      (no cloud ice or bergeron terms)
        2679           0 :         qcsinksum_rate1ord(i,k) = (pra(i,k)+prc(i,k)+psacws(i,k)+psacwg(i,k)+pgsacw(i,k))*lcldm(i,k)
        2680             :         ! Avoid zero/near-zero division.
        2681             :         qcsinksum_rate1ord(i,k) = qcsinksum_rate1ord(i,k) / &
        2682           0 :              max(qc(i,k),1.0e-30_r8)
        2683             :      end do
        2684             :   end do
        2685             :
        2686             :   !$acc loop gang vector collapse(2)
        2687           0 :   do k=1,nlev
        2688           0 :      do i=1,mgncol
        2689             :         ! microphysics output, note this is grid-averaged
        2690           0 :         pratot(i,k)     = pra(i,k)*lcldm(i,k)
        2691           0 :         prctot(i,k)     = prc(i,k)*lcldm(i,k)
        2692           0 :         mnuccctot(i,k)  = mnuccc(i,k)*lcldm(i,k)
        2693           0 :         mnudeptot(i,k)  = mnudep(i,k)*lcldm(i,k)
        2694           0 :         mnuccttot(i,k)  = mnucct(i,k)*lcldm(i,k)
        2695           0 :         msacwitot(i,k)  = msacwi(i,k)*lcldm(i,k)
        2696           0 :         psacwstot(i,k)  = psacws(i,k)*lcldm(i,k)
        2697           0 :         bergstot(i,k)   = bergs(i,k)*lcldm(i,k)
        2698           0 :         bergtot(i,k)    = berg(i,k)
        2699           0 :         prcitot(i,k)    = prci(i,k)*icldm(i,k)
        2700           0 :         praitot(i,k)    = prai(i,k)*icldm(i,k)
        2701           0 :         mnuccdtot(i,k)  = mnuccd(i,k)*icldm(i,k)
        2702           0 :         pracstot(i,k)   = pracs(i,k)*precip_frac(i,k)
        2703           0 :         mnuccrtot(i,k)  = mnuccr(i,k)*precip_frac(i,k)
        2704           0 :         mnuccritot(i,k) = mnuccri(i,k)*precip_frac(i,k)
        2705             :      end do
        2706             :   end do
        2707             :
        2708             :   !$acc loop gang vector collapse(2)
        2709           0 :   do k=1,nlev
        2710           0 :      do i=1,mgncol
        2711           0 :         psacrtot(i,k)   = psacr(i,k)*precip_frac(i,k)
        2712           0 :         pracgtot(i,k)   = pracg(i,k)*precip_frac(i,k)
        2713           0 :         psacwgtot(i,k)  = psacwg(i,k)*lcldm(i,k)
        2714           0 :         pgsacwtot(i,k)  = pgsacw(i,k)*lcldm(i,k)
        2715           0 :         pgracstot(i,k)  = pgracs(i,k)*precip_frac(i,k)
        2716           0 :         prdgtot(i,k)    = prdg(i,k)*precip_frac(i,k)
        2717           0 :         qmultgtot(i,k)  = qmultg(i,k)*lcldm(i,k)
        2718           0 :         qmultrgtot(i,k) = qmultrg(i,k)*precip_frac(i,k)
        2719           0 :         npracgtot(i,k)  = npracg(i,k)*precip_frac(i,k)
        2720           0 :         nscngtot(i,k)   = nscng(i,k)*lcldm(i,k)
        2721           0 :         ngracstot(i,k)  = ngracs(i,k)*precip_frac(i,k)
        2722           0 :         nmultgtot(i,k)  = nmultg(i,k)*lcldm(i,k)
        2723           0 :         nmultrgtot(i,k) = nmultrg(i,k)*precip_frac(i,k)
        2724           0 :         npsacwgtot(i,k) = npsacwg(i,k)*lcldm(i,k)
        2725             :      end do
        2726             :   end do
        2727             :
        2728             :   !$acc loop gang vector collapse(2) private(tmpfrz)
        2729           0 :   do k=1,nlev
        2730           0 :      do i=1,mgncol
        2731           0 :         nctend(i,k) = nctend(i,k)+&
        2732             :              (-nnuccc(i,k)-nnucct(i,k)-npsacws(i,k)+nsubc(i,k) &
        2733           0 :              -npra(i,k)-nprc1(i,k)-npsacwg(i,k))*lcldm(i,k)
        2734             :
        2735           0 :         if (do_cldice) then
        2736           0 :            if (use_hetfrz_classnuc) then
        2737           0 :               tmpfrz = nnuccc(i,k)
        2738             :            else
        2739             :               tmpfrz = 0._r8
        2740             :            end if
        2741             :            nitend(i,k) = nitend(i,k)+ nnuccd(i,k)+ &
        2742             :                 (nnucct(i,k)+tmpfrz+nnudep(i,k)+nsacwi(i,k)+nmultg(i,k))*lcldm(i,k)+(nsubi(i,k)-nprci(i,k)- &
        2743           0 :                 nprai(i,k))*icldm(i,k)+(nnuccri(i,k)+nmultrg(i,k))*precip_frac(i,k)
        2744             :         end if
        2745             :
        2746           0 :         if(do_graupel.or.do_hail) then
        2747             :            nstend(i,k) = nstend(i,k)+(nsubs(i,k)+ &
        2748           0 :                 nsagg(i,k)-ngracs(i,k))*precip_frac(i,k)+nprci(i,k)*icldm(i,k)-nscng(i,k)*lcldm(i,k)
        2749           0 :            ngtend(i,k) = ngtend(i,k)+nscng(i,k)*lcldm(i,k)+(ngracs(i,k)+nnuccr(i,k))*precip_frac(i,k)
        2750             :         else
        2751             :            !necessary since mnuccr moved to graupel
        2752             :            nstend(i,k) = nstend(i,k)+(nsubs(i,k)+ &
        2753           0 :                 nsagg(i,k)+nnuccr(i,k))*precip_frac(i,k)+nprci(i,k)*icldm(i,k)
        2754             :         end if
        2755             :
        2756             :         nrtend(i,k) = nrtend(i,k)+ &
        2757             :              nprc(i,k)*lcldm(i,k)+(nsubr(i,k)-npracs(i,k)-nnuccr(i,k) &
        2758           0 :              -nnuccri(i,k)+nragg(i,k)-npracg(i,k)-ngracs(i,k))*precip_frac(i,k)
        2759             :
        2760             :         ! make sure that ni at advanced time step does not exceed
        2761             :         ! maximum (existing N + source terms*dt), which is possible if mtime < deltat
        2762             :         ! note that currently mtime = deltat
        2763             :         !================================================================
        2764           0 :         if (do_cldice .and. nitend(i,k).gt.0._r8.and.ni(i,k)+nitend(i,k)*deltat.gt.nimax(i,k)) then
        2765           0 :            nitend(i,k)=max(0._r8,(nimax(i,k)-ni(i,k))*rdeltat)
        2766             :         end if
        2767             :      end do
        2768             :   end do ! end k loop
        2769             :   ! End of "administration" loop
        2770             :
        2771             :   !-----------------------------------------------------
        2772             :   ! convert rain/snow q and N for output to history, note,
        2773             :   ! output is for gridbox average
        2774             :
        2775             :   !$acc loop gang vector collapse(2)
        2776           0 :   do k=1,nlev
        2777           0 :      do i=1,mgncol
        2778           0 :         qrout(i,k) = qr(i,k)
        2779           0 :         nrout(i,k) = nr(i,k) * rho(i,k)
        2780           0 :         qsout(i,k) = qs(i,k)
        2781           0 :         nsout(i,k) = ns(i,k) * rho(i,k)
        2782           0 :         qgout(i,k) = qg(i,k)
        2783           0 :         ngout(i,k) = ng(i,k) * rho(i,k)
        2784             :      end do
        2785             :   end do
        2786             :   !$acc end parallel
        2787             :
        2788             :   ! calculate n0r and lamr from rain mass and number
        2789             :   ! divide by precip fraction to get in-precip (local) values of
        2790             :   ! rain mass and number, divide by rhow to get rain number in kg^-1
        2791           0 :   call size_dist_param_basic(mg_rain_props, qric, nric, lamr, mgncol, nlev, n0=n0r)
        2792             :
        2793             :   ! Calculate rercld
        2794             :   ! calculate mean size of combined rain and cloud water
        2795           0 :   call calc_rercld(lamr, n0r, lamc, pgam, qric, qcic, ncic, rercld, mgncol*nlev)
        2796             :
        2797             :   ! Assign variables back to start-of-timestep values
        2798             :   ! Some state variables are changed before the main microphysics loop
        2799             :   ! to make "instantaneous" adjustments. Afterward, we must move those changes
        2800             :   ! back into the tendencies.
        2801             :   ! These processes:
        2802             :   !  - Droplet activation (npccn, impacts nc)
        2803             :   !  - Instantaneous snow melting  (minstsm/ninstsm, impacts qr/qs/nr/ns)
        2804             :   !  - Instantaneous rain freezing (minstfr/ninstrf, impacts qr/qs/nr/ns)
        2805             :   !================================================================================
        2806             :   ! Re-apply droplet activation tendency
        2807             :
        2808             :   !$acc parallel vector_length(VLENS) default(present)
        2809             :   !$acc loop gang vector collapse(2)
        2810           0 :   do k=1,nlev
        2811           0 :      do i=1,mgncol
        2812           0 :         nc(i,k) = ncn(i,k)
        2813           0 :         nctend(i,k) = nctend(i,k) + npccn(i,k)
        2814             :         ! Re-apply rain freezing and snow melting.
        2815           0 :         dum_2D(i,k) = qs(i,k)
        2816           0 :         qs(i,k)     = qsn(i,k)
        2817           0 :         qstend(i,k) = qstend(i,k) + (dum_2D(i,k)-qs(i,k))*rdeltat
        2818             :
        2819           0 :         dum_2D(i,k) = ns(i,k)
        2820           0 :         ns(i,k)     = nsn(i,k)
        2821           0 :         nstend(i,k) = nstend(i,k) + (dum_2D(i,k)-ns(i,k))*rdeltat
        2822             :
        2823           0 :         dum_2D(i,k) = qr(i,k)
        2824           0 :         qr(i,k)     = qrn(i,k)
        2825           0 :         qrtend(i,k) = qrtend(i,k) + (dum_2D(i,k)-qr(i,k))*rdeltat
        2826             :
        2827           0 :         dum_2D(i,k) = nr(i,k)
        2828           0 :         nr(i,k)     = nrn(i,k)
        2829           0 :         nrtend(i,k) = nrtend(i,k) + (dum_2D(i,k)-nr(i,k))*rdeltat
        2830             :
        2831             :         ! Re-apply graupel freezing/melting
        2832           0 :         dum_2D(i,k) = qg(i,k)
        2833           0 :         qg(i,k)     = qgr(i,k)
        2834           0 :         qgtend(i,k) = qgtend(i,k) + (dum_2D(i,k)-qg(i,k))*rdeltat
        2835             :
        2836           0 :         dum_2D(i,k) = ng(i,k)
        2837           0 :         ng(i,k)     = ngr(i,k)
        2838           0 :         ngtend(i,k) = ngtend(i,k) + (dum_2D(i,k)-ng(i,k))*rdeltat
        2839             :         !.............................................................................
        2840             :         !================================================================================
        2841             :         ! modify to include snow. in prain & evap (diagnostic here: for wet dep)
        2842           0 :         nevapr(i,k) = nevapr(i,k) + evapsnow(i,k)
        2843           0 :         prain(i,k) = prain(i,k) + prodsnow(i,k)
        2844             :      end do
        2845             :   end do
        2846             :
        2847             :   !$acc loop gang vector collapse(2)
        2848           0 :   do k=1,nlev
        2849           0 :      do i=1,mgncol
        2850             :         ! calculate sedimentation for cloud water and ice
        2851             :         ! and Graupel (mg3)
        2852             :         !================================================================================
        2853             :         ! update in-cloud cloud mixing ratio and number concentration
        2854             :         ! with microphysical tendencies to calculate sedimentation, assign to dummy vars
        2855             :         ! note: these are in-cloud values***, hence we divide by cloud fraction
        2856           0 :         dumc(i,k)  = (qc(i,k)+qctend(i,k)*deltat)/lcldm(i,k)
        2857           0 :         dumi(i,k)  = (qi(i,k)+qitend(i,k)*deltat)/icldm(i,k)
        2858           0 :         dumnc(i,k) = max((nc(i,k)+nctend(i,k)*deltat)/lcldm(i,k),0._r8)
        2859           0 :         dumni(i,k) = max((ni(i,k)+nitend(i,k)*deltat)/icldm(i,k),0._r8)
        2860             :
        2861           0 :         dumr(i,k)  = (qr(i,k)+qrtend(i,k)*deltat)/precip_frac(i,k)
        2862           0 :         dumnr(i,k) = max((nr(i,k)+nrtend(i,k)*deltat)/precip_frac(i,k),0._r8)
        2863           0 :         dums(i,k)  = (qs(i,k)+qstend(i,k)*deltat)/precip_frac(i,k)
        2864           0 :         dumns(i,k) = max((ns(i,k)+nstend(i,k)*deltat)/precip_frac(i,k),0._r8)
        2865             :
        2866           0 :         dumg(i,k)  = (qg(i,k)+qgtend(i,k)*deltat)/precip_frac(i,k)
        2867           0 :         dumng(i,k) = max((ng(i,k)+ngtend(i,k)*deltat)/precip_frac(i,k),0._r8)
        2868             :
        2869             :         ! switch for specification of droplet and crystal number
        2870           0 :         if (ngcons) then
        2871           0 :            dumng(i,k)=ngnst/rho(i,k)
        2872             :         end if
        2873             :
        2874             :         ! switch for specification of droplet and crystal number
        2875           0 :         if (nccons) then
        2876           0 :            dumnc(i,k)=ncnst/rho(i,k)
        2877             :         end if
        2878             :
        2879             :         ! switch for specification of cloud ice number
        2880           0 :         if (nicons) then
        2881           0 :            dumni(i,k)=ninst/rho(i,k)
        2882             :         end if
        2883             :
        2884             :         ! switch for specification of constant number
        2885           0 :         if (nscons) then
        2886           0 :             dumns(i,k)=nsnst/rho(i,k)
        2887             :         end if
        2888             :
        2889             :         ! switch for specification of constant number
        2890           0 :         if (nrcons) then
        2891           0 :             dumnr(i,k)=nrnst/rho(i,k)
        2892             :         end if
        2893             :      end do
        2894             :   end do
        2895             :   !$acc end parallel
        2896             :
        2897             :   ! obtain new slope parameter to avoid possible singularity
        2898           0 :   call size_dist_param_basic(mg_ice_props, dumi, dumni, lami, mgncol, nlev)
        2899           0 :   call size_dist_param_liq(mg_liq_props, dumc, dumnc, rho, pgam, lamc, mgncol, nlev)
        2900             :
        2901             :   !$acc parallel vector_length(VLENS) default(present)
        2902             :   !$acc loop gang vector collapse(2) private(dum1,dum2,dum3,dum4)
        2903           0 :   do k=1,nlev
        2904           0 :      do i=1,mgncol
        2905             :         ! calculate number and mass weighted fall velocity for droplets and cloud ice
        2906             :         !-------------------------------------------------------------------
        2907           0 :         if (dumc(i,k).ge.qsmall) then
        2908           0 :            dum1 = 4._r8+bc+pgam(i,k)
        2909           0 :            dum2 = pgam(i,k)+4._r8
        2910           0 :            vtrmc(i,k)=acn(i,k)*gamma(dum1)/(lamc(i,k)**bc*gamma(dum2))
        2911             :            ! Following ifs, no condensate sedimentation
        2912           0 :            if (ifs_sed) then
        2913           0 :               fc(i,k)  = 0._r8
        2914           0 :               fnc(i,k) = 0._r8
        2915             :            else
        2916           0 :               dum3     = 1._r8+bc+pgam(i,k)
        2917           0 :               dum4     = pgam(i,k)+1._r8
        2918           0 :               fc(i,k)  = g*rho(i,k)*vtrmc(i,k)
        2919             :               fnc(i,k) = g*rho(i,k)* &
        2920             :                    acn(i,k)*gamma(dum3)/ &
        2921           0 :                    (lamc(i,k)**bc*gamma(dum4))
        2922             :            end if
        2923             :         else
        2924           0 :            fc(i,k) = 0._r8
        2925           0 :            fnc(i,k)= 0._r8
        2926             :         end if
        2927             :      end do
        2928             :   end do
        2929             :
        2930             :   !$acc loop gang vector collapse(2) private(irad,ifrac)
        2931           0 :   do k=1,nlev
        2932           0 :      do i=1,mgncol
        2933             :         ! calculate number and mass weighted fall velocity for cloud ice
        2934           0 :         if (dumi(i,k).ge.qsmall) then
        2935             :            vtrmi(i,k)=min(ain(i,k)*gamma_bi_plus4/(6._r8*lami(i,k)**bi), &
        2936           0 :                 1.2_r8*rhof(i,k))
        2937           0 :            vtrmi(i,k)=vtrmi(i,k)*micro_mg_vtrmi_factor
        2938             :
        2939           0 :            fi(i,k) = g*rho(i,k)*vtrmi(i,k)
        2940             :            fni(i,k) = g*rho(i,k)* &
        2941           0 :                 min(ain(i,k)*gamma_bi_plus1/lami(i,k)**bi,1.2_r8*rhof(i,k))
        2942             :
        2943             :            ! adjust the ice fall velocity for smaller (r < 20 um) ice
        2944             :            ! particles (blend over 8-20 um)
        2945           0 :            irad = 1.5_r8 / lami(i,k) * 1e6_r8
        2946           0 :            ifrac = min(1._r8, max(0._r8, (irad - 18._r8) / 2._r8))
        2947             :
        2948           0 :            if (ifrac .lt. 1._r8) then
        2949             :               vtrmi(i,k) = ifrac * vtrmi(i,k) + &
        2950             :                  (1._r8 - ifrac) * &
        2951             :                  min(ajn(i,k)*gamma_bj_plus4/(6._r8*lami(i,k)**bj), &
        2952           0 :                  1.2_r8*rhof(i,k))
        2953           0 :               vtrmi(i,k)=vtrmi(i,k)*micro_mg_vtrmi_factor
        2954             :
        2955           0 :               fi(i,k)  = g*rho(i,k)*vtrmi(i,k)
        2956             :               fni(i,k) = ifrac * fni(i,k) + &
        2957             :                  (1._r8 - ifrac) * &
        2958             :                  g*rho(i,k)* &
        2959           0 :                  min(ajn(i,k)*gamma_bj_plus1/lami(i,k)**bj,1.2_r8*rhof(i,k))
        2960             :            end if
        2961             :
        2962             :            ! Fix ice fall speed following IFS microphysics
        2963           0 :            if (ifs_sed) then
        2964           0 :               fi(i,k)=g*rho(i,k)*0.1_r8
        2965           0 :               fni(i,k)=g*rho(i,k)*0.1_r8
        2966             :            end if
        2967             :         else
        2968           0 :            fi(i,k) = 0._r8
        2969           0 :            fni(i,k)= 0._r8
        2970             :         end if
        2971             :      end do
        2972             :   end do
        2973             :   !$acc end parallel
        2974             :
        2975             :   ! fallspeed for rain
        2976           0 :   call size_dist_param_basic(mg_rain_props, dumr, dumnr, lamr, mgncol, nlev)
        2977             :   ! fallspeed for snow
        2978           0 :   call size_dist_param_basic(mg_snow_props, dums, dumns, lams, mgncol, nlev)
        2979             :   ! fallspeed for graupel
        2980           0 :   if (do_hail) then
        2981           0 :      call size_dist_param_basic(mg_hail_props, dumg, dumng, lamg, mgncol, nlev)
        2982             :   end if
        2983           0 :   if (do_graupel) then
        2984           0 :      call size_dist_param_basic(mg_graupel_props, dumg, dumng, lamg, mgncol, nlev)
        2985             :   end if
        2986             :
        2987             :   !$acc parallel vector_length(VLENS) default(present)
        2988             :   !$acc loop gang vector collapse(2) private(qtmp)
        2989           0 :   do k=1,nlev
        2990           0 :      do i=1,mgncol
        2991           0 :         if (lamr(i,k).ge.qsmall) then
        2992           0 :            qtmp = lamr(i,k)**br
        2993             :            ! 'final' values of number and mass weighted mean fallspeed for rain (m/s)
        2994           0 :            unr(i,k) = min(arn(i,k)*gamma_br_plus1/qtmp,9.1_r8*rhof(i,k))
        2995           0 :            fnr(i,k) = g*rho(i,k)*unr(i,k)
        2996           0 :            umr(i,k) = min(arn(i,k)*gamma_br_plus4/(6._r8*qtmp),9.1_r8*rhof(i,k))
        2997           0 :            fr(i,k) = g*rho(i,k)*umr(i,k)
        2998             :         else
        2999           0 :            fr(i,k)=0._r8
        3000           0 :            fnr(i,k)=0._r8
        3001             :         end if
        3002             :
        3003             :         ! Fallspeed correction to ensure non-zero if rain in the column
        3004             :         ! from updated Morrison (WRFv3.3) and P3 schemes
        3005             :         ! If fallspeed exists at a higher level, apply it below to eliminate
        3006           0 :         if (precip_fall_corr) then
        3007           0 :            if (k.gt.2) then
        3008           0 :               if (fr(i,k).lt.1.e-10_r8) then
        3009           0 :                  fr(i,k)=fr(i,k-1)
        3010           0 :                  fnr(i,k)=fnr(i,k-1)
        3011             :               end if
        3012             :            end if
        3013             :         end if
        3014             :
        3015           0 :         if (lams(i,k).ge.qsmall) then
        3016           0 :            qtmp = lams(i,k)**bs
        3017             :            ! 'final' values of number and mass weighted mean fallspeed for snow (m/s)
        3018           0 :            ums(i,k) = min(asn(i,k)*gamma_bs_plus4/(6._r8*qtmp),1.2_r8*rhof(i,k))
        3019           0 :            ums(:,k)=ums(:,k)*micro_mg_vtrmi_factor
        3020             :
        3021           0 :            fs(i,k)  = g*rho(i,k)*ums(i,k)
        3022           0 :            uns(i,k) = min(asn(i,k)*gamma_bs_plus1/qtmp,1.2_r8*rhof(i,k))
        3023           0 :            fns(i,k) = g*rho(i,k)*uns(i,k)
        3024             :            ! Fix fallspeed for snow
        3025           0 :            if (ifs_sed) then
        3026           0 :               ums(i,k) = 1._r8
        3027           0 :               uns(i,k) = 1._r8
        3028             :             end if
        3029             :         else
        3030           0 :            fs(i,k)=0._r8
        3031           0 :            fns(i,k)=0._r8
        3032             :         end if
        3033             :
        3034           0 :         if (precip_fall_corr) then
        3035           0 :            if (k.gt.2) then
        3036           0 :               if (fs(i,k).lt.1.e-10_r8) then
        3037           0 :                  fs(i,k)=fs(i,k-1)
        3038           0 :                  fns(i,k)=fns(i,k-1)
        3039             :               end if
        3040             :            end if
        3041             :         end if
        3042             :
        3043           0 :         if (lamg(i,k).ge.qsmall) then
        3044           0 :            qtmp = lamg(i,k)**bgtmp
        3045             :            ! 'final' values of number and mass weighted mean fallspeed for graupel (m/s)
        3046           0 :            umg(i,k) = min(agn(i,k)*gamma_bg_plus4/(6._r8*qtmp),20._r8*rhof(i,k))
        3047           0 :            fg(i,k) = g*rho(i,k)*umg(i,k)
        3048           0 :            ung(i,k) = min(agn(i,k)*gamma_bg_plus1/qtmp,20._r8*rhof(i,k))
        3049           0 :            fng(i,k) = g*rho(i,k)*ung(i,k)
        3050             :         else
        3051           0 :            fg(i,k)=0._r8
        3052           0 :            fng(i,k)=0._r8
        3053             :         end if
        3054             :
        3055           0 :         if (precip_fall_corr) then
        3056           0 :            if (k.gt.2) then
        3057           0 :               if (fg(i,k).lt.1.e-10_r8) then
        3058           0 :                  fg(i,k)=fg(i,k-1)
        3059           0 :                  fng(i,k)=fng(i,k-1)
        3060             :               end if
        3061             :            end if
        3062             :         end if
        3063           0 :         pdel_inv(i,k) = 1._r8/pdel(i,k)
        3064             :      end do
        3065             :   end do
        3066             :
        3067             :   !$acc loop gang vector collapse(2)
        3068           0 :   do k=1,nlev
        3069           0 :      do i=1,mgncol
        3070             :         ! redefine dummy variables - sedimentation is calculated over grid-scale
        3071             :         ! quantities to ensure conservation
        3072           0 :         dumc(i,k) = (qc(i,k)+qctend(i,k)*deltat)
        3073           0 :         dumnc(i,k) = max((nc(i,k)+nctend(i,k)*deltat),0._r8)
        3074           0 :         dumi(i,k) = (qi(i,k)+qitend(i,k)*deltat)
        3075           0 :         dumni(i,k) = max((ni(i,k)+nitend(i,k)*deltat),0._r8)
        3076           0 :         dumr(i,k) = (qr(i,k)+qrtend(i,k)*deltat)
        3077           0 :         dumnr(i,k) = max((nr(i,k)+nrtend(i,k)*deltat),0._r8)
        3078           0 :         dums(i,k) = (qs(i,k)+qstend(i,k)*deltat)
        3079           0 :         dumns(i,k) = max((ns(i,k)+nstend(i,k)*deltat),0._r8)
        3080           0 :         dumg(i,k) = (qg(i,k)+qgtend(i,k)*deltat)
        3081           0 :         dumng(i,k) = max((ng(i,k)+ngtend(i,k)*deltat),0._r8)
        3082           0 :         if (dumc(i,k).lt.qsmall) dumnc(i,k)=0._r8
        3083           0 :         if (dumi(i,k).lt.qsmall) dumni(i,k)=0._r8
        3084           0 :         if (dumr(i,k).lt.qsmall) dumnr(i,k)=0._r8
        3085           0 :         if (dums(i,k).lt.qsmall) dumns(i,k)=0._r8
        3086           0 :         if (dumg(i,k).lt.qsmall) dumng(i,k)=0._r8
        3087             :      end do
        3088             :   end do
        3089             :   !$acc end parallel
        3090             :
        3091             :   ! begin sedimentation
        3092             :   ! ice
        3093             :   call Sedimentation(mgncol,nlev,do_cldice,deltat,fi,fni,pdel_inv, &
        3094             :                        qitend,nitend,qisedten,dumi,dumni,prect,iflx, &
        3095             :                        xxlx=xxls,qxsevap=qisevap,tlat=tlat,qvlat=qvlat, &
        3096           0 :                        xcldm=icldm,preci=preci)
        3097             :   ! liq
        3098             :   call Sedimentation(mgncol,nlev,.TRUE.,deltat,fc,fnc,pdel_inv, &
        3099             :                        qctend,nctend,qcsedten,dumc,dumnc,prect,lflx, &
        3100           0 :                        xxlx=xxlv,qxsevap=qcsevap,tlat=tlat,qvlat=qvlat,xcldm=lcldm)
        3101             :   ! rain
        3102             :   call Sedimentation(mgncol,nlev,.TRUE.,deltat,fr,fnr,pdel_inv, &
        3103           0 :                        qrtend,nrtend,qrsedten,dumr,dumnr,prect,rflx)
        3104             :   ! snow
        3105             :   call Sedimentation(mgncol,nlev,.TRUE.,deltat,fs,fns,pdel_inv, &
        3106           0 :                        qstend,nstend,qssedten,dums,dumns,prect,sflx,preci=preci)
        3107             :   ! graupel
        3108             :   call Sedimentation(mgncol,nlev,.TRUE.,deltat,fg,fng,pdel_inv, &
        3109           0 :                        qgtend,ngtend,qgsedten,dumg,dumng,prect,gflx,preci=preci)
        3110             :   ! end sedimentation
        3111             :
        3112             :   !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        3113             :
        3114             :   ! get new update for variables that includes sedimentation tendency
        3115             :   ! note : here dum variables are grid-average, NOT in-cloud
        3116             :
        3117             :   !$acc parallel vector_length(VLENS) default(present)
        3118             :   !$acc loop gang vector collapse(2)
        3119           0 :   do k=1,nlev
        3120           0 :      do i=1,mgncol
        3121           0 :         dumc(i,k)  = max(qc(i,k)+qctend(i,k)*deltat,0._r8)
        3122           0 :         dumi(i,k)  = max(qi(i,k)+qitend(i,k)*deltat,0._r8)
        3123           0 :         dumnc(i,k) = max(nc(i,k)+nctend(i,k)*deltat,0._r8)
        3124           0 :         dumni(i,k) = max(ni(i,k)+nitend(i,k)*deltat,0._r8)
        3125             :
        3126           0 :         dumr(i,k)  = max(qr(i,k)+qrtend(i,k)*deltat,0._r8)
        3127           0 :         dumnr(i,k) = max(nr(i,k)+nrtend(i,k)*deltat,0._r8)
        3128           0 :         dums(i,k)  = max(qs(i,k)+qstend(i,k)*deltat,0._r8)
        3129           0 :         dumns(i,k) = max(ns(i,k)+nstend(i,k)*deltat,0._r8)
        3130           0 :         dumg(i,k)  = max(qg(i,k)+qgtend(i,k)*deltat,0._r8)
        3131           0 :         dumng(i,k) = max(ng(i,k)+ngtend(i,k)*deltat,0._r8)
        3132             :
        3133             :         ! switch for specification of droplet and crystal number
        3134           0 :         if (nccons) then
        3135           0 :            dumnc(i,k)=ncnst/rho(i,k)*lcldm(i,k)
        3136             :         end if
        3137             :
        3138             :         ! switch for specification of cloud ice number
        3139           0 :         if (nicons) then
        3140           0 :            dumni(i,k)=ninst/rho(i,k)*icldm(i,k)
        3141             :         end if
        3142             :
        3143             :         ! switch for specification of graupel number
        3144           0 :         if (ngcons) then
        3145           0 :            dumng(i,k)=ngnst/rho(i,k)*precip_frac(i,k)
        3146             :         end if
        3147             :
        3148             :         ! switch for specification of constant snow number
        3149           0 :         if (nscons) then
        3150           0 :             dumns(i,k)=nsnst/rho(i,k)
        3151             :         end if
        3152             :
        3153             :         ! switch for specification of constant rain number
        3154           0 :         if (nrcons) then
        3155           0 :             dumnr(i,k)=nrnst/rho(i,k)
        3156             :         end if
        3157             :
        3158           0 :         if (dumc(i,k).lt.qsmall) dumnc(i,k)=0._r8
        3159           0 :         if (dumi(i,k).lt.qsmall) dumni(i,k)=0._r8
        3160           0 :         if (dumr(i,k).lt.qsmall) dumnr(i,k)=0._r8
        3161           0 :         if (dums(i,k).lt.qsmall) dumns(i,k)=0._r8
        3162           0 :         if (dumg(i,k).lt.qsmall) dumng(i,k)=0._r8
        3163             :      end do
        3164             :   end do
        3165             :
        3166             :   ! calculate instantaneous processes (melting, homogeneous freezing)
        3167             :   !====================================================================
        3168             :   ! melting of snow at +2 C
        3169             :
        3170             :   !$acc loop gang vector collapse(2) private(dum,dum1)
        3171           0 :   do k=1,nlev
        3172           0 :      do i=1,mgncol
        3173           0 :         if (t(i,k)+tlat(i,k)/cpp*deltat > snowmelt) then
        3174           0 :            if (dums(i,k) > 0._r8) then
        3175             :               ! make sure melting snow doesn't reduce temperature below threshold
        3176           0 :               dum = -xlf/cpp*dums(i,k)
        3177           0 :               if (t(i,k)+tlat(i,k)/cpp*deltat+dum.lt. snowmelt) then
        3178           0 :                  dum = (t(i,k)+tlat(i,k)/cpp*deltat-snowmelt)*cpp/xlf
        3179           0 :                  dum = dum/dums(i,k)
        3180           0 :                  dum = max(0._r8,dum)
        3181           0 :                  dum = min(1._r8,dum)
        3182             :               else
        3183             :                  dum = 1._r8
        3184             :               end if
        3185             :
        3186           0 :               qstend(i,k)=qstend(i,k)-dum*dums(i,k)*rdeltat
        3187           0 :               nstend(i,k)=nstend(i,k)-dum*dumns(i,k)*rdeltat
        3188           0 :               qrtend(i,k)=qrtend(i,k)+dum*dums(i,k)*rdeltat
        3189           0 :               nrtend(i,k)=nrtend(i,k)+dum*dumns(i,k)*rdeltat
        3190             :
        3191           0 :               dum1=-xlf*dum*dums(i,k)*rdeltat
        3192           0 :               tlat(i,k)=tlat(i,k)+dum1
        3193           0 :               meltsdttot(i,k)=meltsdttot(i,k) + dum1
        3194             :
        3195             : !STOPPED FIX FOR SNOW NUMBER
        3196             : !ensure that snow... number does not go negative with constant number set
        3197             : !necessary because dumng is updated above.
        3198           0 :               if (nscons .and. ((ns(i,k)+nstend(i,k)*deltat) .lt. 0._r8)) then
        3199           0 :                  nstend(i,k)=-ns(i,k)*rdeltat
        3200             :               end if
        3201             :            end if
        3202             :         end if
        3203             :      end do
        3204             :   end do
        3205             :
        3206             :   ! melting of graupel at +2 C
        3207             :
        3208             :   !$acc loop gang vector collapse(2) private(dum,dum1)
        3209           0 :   do k=1,nlev
        3210           0 :      do i=1,mgncol
        3211           0 :         if (t(i,k)+tlat(i,k)/cpp*deltat > snowmelt) then
        3212           0 :            if (dumg(i,k) > 0._r8) then
        3213             :               ! make sure melting graupel doesn't reduce temperature below threshold
        3214           0 :               dum = -xlf/cpp*dumg(i,k)
        3215           0 :               if (t(i,k)+tlat(i,k)/cpp*deltat+dum .lt. snowmelt) then
        3216           0 :                  dum = (t(i,k)+tlat(i,k)/cpp*deltat-snowmelt)*cpp/xlf
        3217           0 :                  dum = dum/dumg(i,k)
        3218           0 :                  dum = max(0._r8,dum)
        3219           0 :                  dum = min(1._r8,dum)
        3220             :               else
        3221             :                  dum = 1._r8
        3222             :               end if
        3223             :
        3224           0 :               qgtend(i,k)=qgtend(i,k)-dum*dumg(i,k)*rdeltat
        3225           0 :               ngtend(i,k)=ngtend(i,k)-dum*dumng(i,k)*rdeltat
        3226           0 :               qrtend(i,k)=qrtend(i,k)+dum*dumg(i,k)*rdeltat
        3227           0 :               nrtend(i,k)=nrtend(i,k)+dum*dumng(i,k)*rdeltat
        3228             :
        3229           0 :               dum1=-xlf*dum*dumg(i,k)*rdeltat
        3230           0 :               tlat(i,k)=tlat(i,k)+dum1
        3231           0 :               meltsdttot(i,k)=meltsdttot(i,k) + dum1
        3232             :
        3233             : !ensure that graupel number does not go negative with constant number set
        3234             : !necessary because dumng is updated above.
        3235           0 :               if (ngcons .and. ((ng(i,k)+ngtend(i,k)*deltat) .lt. 0._r8)) then
        3236           0 :                  ngtend(i,k)=-ng(i,k)*rdeltat
        3237             :               end if
        3238             :            end if
        3239             :         end if
        3240             :      end do
        3241             :   end do
        3242             :   !$acc end parallel
        3243             :
        3244             :   ! get mean size of rain = 1/lamr, add frozen rain to either snow or cloud ice
        3245             :   ! depending on mean rain size
        3246             :   ! add to graupel if using that option....
        3247           0 :   call size_dist_param_basic(mg_rain_props, dumr, dumnr, lamr, mgncol, nlev)
        3248             :
        3249             :   !$acc parallel vector_length(VLENS) default(present)
        3250             :   !$acc loop gang vector collapse(2) private(dum,dum1)
        3251           0 :   do k=1,nlev
        3252           0 :      do i=1,mgncol
        3253             :         ! freezing of rain at -5 C
        3254           0 :         if (t(i,k)+tlat(i,k)/cpp*deltat < rainfrze) then
        3255           0 :            if (dumr(i,k) > 0._r8) then
        3256             :               ! make sure freezing rain doesn't increase temperature above threshold
        3257           0 :               dum = xlf/cpp*dumr(i,k)
        3258           0 :               if (t(i,k)+tlat(i,k)/cpp*deltat+dum.gt.rainfrze) then
        3259           0 :                  dum = -(t(i,k)+tlat(i,k)/cpp*deltat-rainfrze)*cpp/xlf
        3260           0 :                  dum = dum/dumr(i,k)
        3261           0 :                  dum = max(0._r8,dum)
        3262           0 :                  dum = min(1._r8,dum)
        3263             :               else
        3264             :                  dum = 1._r8
        3265             :               end if
        3266             :
        3267           0 :               qrtend(i,k)=qrtend(i,k)-dum*dumr(i,k)*rdeltat
        3268           0 :               nrtend(i,k)=nrtend(i,k)-dum*dumnr(i,k)*rdeltat
        3269             :
        3270           0 :               if (lamr(i,k) < 1._r8/Dcs) then
        3271           0 :                  if (do_hail.or.do_graupel) then
        3272           0 :                     qgtend(i,k)=qgtend(i,k)+dum*dumr(i,k)*rdeltat
        3273           0 :                     ngtend(i,k)=ngtend(i,k)+dum*dumnr(i,k)*rdeltat
        3274             :                  else
        3275           0 :                     qstend(i,k)=qstend(i,k)+dum*dumr(i,k)*rdeltat
        3276           0 :                     nstend(i,k)=nstend(i,k)+dum*dumnr(i,k)*rdeltat
        3277             :                  end if
        3278             :               else
        3279           0 :                  qitend(i,k)=qitend(i,k)+dum*dumr(i,k)*rdeltat
        3280           0 :                  nitend(i,k)=nitend(i,k)+dum*dumnr(i,k)*rdeltat
        3281             :               end if
        3282             :
        3283             :               ! heating tendency
        3284           0 :               dum1 = xlf*dum*dumr(i,k)*rdeltat
        3285           0 :               frzrdttot(i,k)=frzrdttot(i,k) + dum1
        3286           0 :               tlat(i,k)=tlat(i,k)+dum1
        3287             :            end if
        3288             :         end if
        3289             :       end do
        3290             :    end do
        3291             :    !$acc end parallel
        3292             :
        3293           0 :    if (do_cldice) then
        3294             :       !$acc parallel vector_length(VLENS) default(present)
        3295             :       !$acc loop gang vector collapse(2) private(dum)
        3296           0 :       do k=1,nlev
        3297           0 :         do i=1,mgncol
        3298           0 :            if (t(i,k)+tlat(i,k)/cpp*deltat > tmelt) then
        3299           0 :               if (dumi(i,k) > 0._r8) then
        3300             :                  ! limit so that melting does not push temperature below freezing
        3301             :                  !-----------------------------------------------------------------
        3302           0 :                  dum = -dumi(i,k)*xlf/cpp
        3303           0 :                  if (t(i,k)+tlat(i,k)/cpp*deltat+dum.lt.tmelt) then
        3304           0 :                     dum = (t(i,k)+tlat(i,k)/cpp*deltat-tmelt)*cpp/xlf
        3305           0 :                     dum = dum/dumi(i,k)
        3306           0 :                     dum = max(0._r8,dum)
        3307           0 :                     dum = min(1._r8,dum)
        3308             :                  else
        3309             :                     dum = 1._r8
        3310             :                  end if
        3311             :
        3312           0 :                  qctend(i,k)=qctend(i,k)+dum*dumi(i,k)*rdeltat
        3313             :
        3314             :                  ! for output
        3315           0 :                  melttot(i,k)=dum*dumi(i,k)*rdeltat
        3316             :
        3317             :                  ! assume melting ice produces droplet
        3318             :                  ! mean volume radius of 8 micron
        3319             :                  nctend(i,k)=nctend(i,k)+3._r8*dum*dumi(i,k)*rdeltat/ &
        3320           0 :                       (4._r8*pi*5.12e-16_r8*rhow)
        3321             :
        3322           0 :                  qitend(i,k)=((1._r8-dum)*dumi(i,k)-qi(i,k))*rdeltat
        3323           0 :                  nitend(i,k)=((1._r8-dum)*dumni(i,k)-ni(i,k))*rdeltat
        3324           0 :                  tlat(i,k)=tlat(i,k)-xlf*dum*dumi(i,k)*rdeltat
        3325             :               end if
        3326             :            end if
        3327             :         end do
        3328             :      end do
        3329             :
        3330             :      ! homogeneously freeze droplets at -40 C
        3331             :      !-----------------------------------------------------------------
        3332             :
        3333             :      !$acc loop gang vector collapse(2) private(dum)
        3334           0 :      do k=1,nlev
        3335           0 :         do i=1,mgncol
        3336           0 :            if (t(i,k)+tlat(i,k)/cpp*deltat < 233.15_r8) then
        3337           0 :               if (dumc(i,k) > 0._r8) then
        3338             :                  ! limit so that freezing does not push temperature above threshold
        3339           0 :                  dum = dumc(i,k)*xlf/cpp
        3340           0 :                  if (t(i,k)+tlat(i,k)/cpp*deltat+dum.gt.233.15_r8) then
        3341           0 :                     dum = -(t(i,k)+tlat(i,k)/cpp*deltat-233.15_r8)*cpp/xlf
        3342           0 :                     dum = dum/dumc(i,k)
        3343           0 :                     dum = max(0._r8,dum)
        3344           0 :                     dum = min(1._r8,dum)
        3345             :                  else
        3346             :                     dum = 1._r8
        3347             :                  end if
        3348             :
        3349           0 :                  qitend(i,k)=qitend(i,k)+dum*dumc(i,k)*rdeltat
        3350             :                  ! for output
        3351           0 :                  homotot(i,k)=dum*dumc(i,k)*rdeltat
        3352             :
        3353             :                  ! assume 25 micron mean volume radius of homogeneously frozen droplets
        3354             :                  ! consistent with size of detrained ice in stratiform.F90
        3355           0 :                  nitend(i,k)=nitend(i,k)+dum*3._r8*dumc(i,k)/(4._r8*3.14_r8*micro_mg_homog_size**3._r8*500._r8)*rdeltat
        3356             :
        3357           0 :                  qctend(i,k)=((1._r8-dum)*dumc(i,k)-qc(i,k))*rdeltat
        3358           0 :                  nctend(i,k)=((1._r8-dum)*dumnc(i,k)-nc(i,k))*rdeltat
        3359           0 :                  tlat(i,k)=tlat(i,k)+xlf*dum*dumc(i,k)*rdeltat
        3360             :               end if
        3361             :            end if
        3362             :         end do
        3363             :      end do
        3364             :
        3365             :      ! ice number limiter
        3366           0 :      do k=1,nlev
        3367           0 :         do i=1,mgncol
        3368           0 :             if (do_cldice .and. nitend(i,k).gt.0._r8.and.ni(i,k)+nitend(i,k)*deltat.gt.micro_mg_max_nicons/rho(i,k)) then
        3369           0 :                nitend(i,k)=max(0._r8,(micro_mg_max_nicons/rho(i,k)-ni(i,k))/deltat)
        3370             :             end if
        3371             :         end do
        3372             :      end do
        3373             :
        3374             :      ! remove any excess over-saturation, which is possible due to non-linearity when adding
        3375             :      ! together all microphysical processes
        3376             :      !-----------------------------------------------------------------
        3377             :      ! follow code similar to old CAM scheme
        3378             :
        3379             :      !$acc loop gang vector collapse(2)
        3380           0 :      do k=1,nlev
        3381           0 :         do i=1,mgncol
        3382           0 :            dum_2D(i,k)=q(i,k)+qvlat(i,k)*deltat
        3383           0 :            ttmpA(i,k)=t(i,k)+tlat(i,k)/cpp*deltat
        3384             :         end do
        3385             :      end do
        3386             :      !$acc end parallel
        3387             :
        3388             :      ! use rhw to allow ice supersaturation
        3389           0 :      call qsat_water(ttmpA, p, esnA, qvnA, mgncol*nlev)
        3390             :
        3391             :      !$acc parallel vector_length(VLENS) default(present)
        3392             :      !$acc loop gang vector collapse(2) private(dum,dum1)
        3393           0 :      do k=1,nlev
        3394           0 :         do i=1,mgncol
        3395           0 :            if (dum_2D(i,k) > qvnA(i,k) .and. qvnA(i,k) > 0 .and. remove_supersat) then
        3396             :               ! expression below is approximate since there may be ice deposition
        3397           0 :               dum = (dum_2D(i,k)-qvnA(i,k))/(1._r8+xxlv_squared*qvnA(i,k)/(cpp*rv*ttmpA(i,k)**2))*rdeltat
        3398             :               ! add to output cme
        3399           0 :               cmeout(i,k) = cmeout(i,k)+dum
        3400             :               ! now add to tendencies, partition between liquid and ice based on temperature
        3401           0 :               if (ttmpA(i,k) > 268.15_r8) then
        3402           0 :                  dum1=0.0_r8
        3403             :                  ! now add to tendencies, partition between liquid and ice based on te
        3404             :                  !-------------------------------------------------------
        3405           0 :               else if (ttmpA(i,k) < 238.15_r8) then
        3406           0 :                  dum1=1.0_r8
        3407             :               else
        3408           0 :                  dum1=(268.15_r8-ttmpA(i,k))/30._r8
        3409             :               end if
        3410             :               dum = (dum_2D(i,k)-qvnA(i,k))/(1._r8+(xxls*dum1+xxlv*(1._r8-dum1))**2 &
        3411           0 :                     *qvnA(i,k)/(cpp*rv*ttmpA(i,k)**2))*rdeltat
        3412           0 :               qctend(i,k)=qctend(i,k)+dum*(1._r8-dum1)
        3413             :               ! for output
        3414           0 :               qcrestot(i,k)=dum*(1._r8-dum1)
        3415           0 :               qitend(i,k)=qitend(i,k)+dum*dum1
        3416           0 :               qirestot(i,k)=dum*dum1
        3417           0 :               qvlat(i,k)=qvlat(i,k)-dum
        3418             :               ! for output
        3419           0 :               qvres(i,k)=-dum
        3420           0 :               tlat(i,k)=tlat(i,k)+dum*(1._r8-dum1)*xxlv+dum*dum1*xxls
        3421             :            end if
        3422             :         end do
        3423             :      end do
        3424             :      !$acc end parallel
        3425             :   end if
        3426             :
        3427             :   ! calculate effective radius for pass to radiation code
        3428             :   !=========================================================
        3429             :   ! if no cloud water, default value is 10 micron for droplets,
        3430             :   ! 25 micron for cloud ice
        3431             :
        3432             :   ! update cloud variables after instantaneous processes to get effective radius
        3433             :   ! variables are in-cloud to calculate size dist parameters
        3434             :
        3435             :   !$acc parallel vector_length(VLENS) default(present)
        3436             :   !$acc loop gang vector collapse(2)
        3437           0 :   do k=1,nlev
        3438           0 :      do i=1,mgncol
        3439           0 :         dumc(i,k) = max(qc(i,k)+qctend(i,k)*deltat,0._r8)/lcldm(i,k)
        3440           0 :         dumi(i,k) = max(qi(i,k)+qitend(i,k)*deltat,0._r8)/icldm(i,k)
        3441           0 :         dumnc(i,k) = max(nc(i,k)+nctend(i,k)*deltat,0._r8)/lcldm(i,k)
        3442           0 :         dumni(i,k) = max(ni(i,k)+nitend(i,k)*deltat,0._r8)/icldm(i,k)
        3443             :
        3444           0 :         dumr(i,k) = max(qr(i,k)+qrtend(i,k)*deltat,0._r8)/precip_frac(i,k)
        3445           0 :         dumnr(i,k) = max(nr(i,k)+nrtend(i,k)*deltat,0._r8)/precip_frac(i,k)
        3446           0 :         dums(i,k) = max(qs(i,k)+qstend(i,k)*deltat,0._r8)/precip_frac(i,k)
        3447           0 :         dumns(i,k) = max(ns(i,k)+nstend(i,k)*deltat,0._r8)/precip_frac(i,k)
        3448           0 :         dumg(i,k) = max(qg(i,k)+qgtend(i,k)*deltat,0._r8)
        3449           0 :         dumng(i,k) = max(ng(i,k)+ngtend(i,k)*deltat,0._r8)
        3450             :
        3451             :         ! switch for specification of droplet and crystal number
        3452           0 :         if (nccons) then
        3453           0 :            dumnc(i,k)=ncnst/rho(i,k)
        3454             :         end if
        3455             :
        3456             :         ! switch for specification of cloud ice number
        3457           0 :         if (nicons) then
        3458           0 :            dumni(i,k)=ninst/rho(i,k)
        3459             :         end if
        3460             :
        3461             :         ! switch for specification of graupel number
        3462           0 :         if (ngcons) then
        3463           0 :            dumng(i,k)=ngnst/rho(i,k)*precip_frac(i,k)
        3464             :         end if
        3465             :
        3466             :         ! switch for specification of constant snow number
        3467           0 :         if (nscons) then
        3468           0 :             dumns(i,k)=nsnst/rho(i,k)
        3469             :         end if
        3470             :
        3471             :         ! switch for specification of constant rain number
        3472           0 :         if (nrcons) then
        3473           0 :             dumnr(i,k)=nrnst/rho(i,k)
        3474             :         end if
        3475             :
        3476             :         ! limit in-cloud mixing ratio to reasonable value of 5 g kg-1
        3477           0 :         dumc(i,k)=min(dumc(i,k),5.e-3_r8)
        3478           0 :         dumi(i,k)=min(dumi(i,k),5.e-3_r8)
        3479             :         ! limit in-precip mixing ratios
        3480           0 :         dumr(i,k)=min(dumr(i,k),10.e-3_r8)
        3481           0 :         dums(i,k)=min(dums(i,k),10.e-3_r8)
        3482           0 :         dumg(i,k)=min(dumg(i,k),10.e-3_r8)
        3483             :      end do
        3484             :   end do
        3485             :   !$acc end parallel
        3486             :
        3487             :   ! cloud ice effective radius
        3488             :   !-----------------------------------------------------------------
        3489           0 :   if (do_cldice) then
        3490             :      !$acc parallel vector_length(VLENS) default(present)
        3491             :      !$acc loop gang vector collapse(2)
        3492           0 :      do k=1,nlev
        3493           0 :         do i=1,mgncol
        3494           0 :            dum_2D(i,k) = dumni(i,k)
        3495             :         end do
        3496             :      end do
        3497             :      !$acc end parallel
        3498             :
        3499           0 :      call size_dist_param_basic(mg_ice_props, dumi, dumni, lami, mgncol, nlev, n0=dumni0A2D)
        3500             :
        3501             :      !$acc parallel vector_length(VLENS) default(present)
        3502             :      !$acc loop gang vector collapse(2)
        3503           0 :      do k=1,nlev
        3504           0 :         do i=1,mgncol
        3505           0 :            if (dumi(i,k).ge.qsmall) then
        3506           0 :               if (dumni(i,k) /=dum_2D(i,k)) then
        3507             :                  ! adjust number conc if needed to keep mean size in reasonable range
        3508           0 :                  nitend(i,k)=(dumni(i,k)*icldm(i,k)-ni(i,k))*rdeltat
        3509             :               end if
        3510           0 :               effi(i,k)   = 1.5_r8/lami(i,k)*1.e6_r8
        3511           0 :               effi(i,k) = effi(i,k)*micro_mg_effi_factor
        3512             :
        3513           0 :               sadice(i,k) = 2._r8*pi*(lami(i,k)**(-3))*dumni0A2D(i,k)*rho(i,k)*1.e-2_r8  ! m2/m3 -> cm2/cm3
        3514             :            else
        3515           0 :               effi(i,k)   = 25._r8
        3516           0 :               effi(i,k) = effi(i,k)*micro_mg_effi_factor
        3517             :
        3518           0 :               sadice(i,k) = 0._r8
        3519             :            end if
        3520             :            ! ice effective diameter for david mitchell's optics
        3521           0 :            deffi(i,k)=effi(i,k)*rhoi/rhows*2._r8
        3522             :         end do
        3523             :      end do
        3524             :      !$acc end parallel
        3525             :   else
        3526             :      !$acc parallel vector_length(VLENS) default(present)
        3527             :      !acc loop gang vector collapse(2)
        3528           0 :      do k=1,nlev
        3529           0 :         do i=1,mgncol
        3530             :            ! NOTE: If CARMA is doing the ice microphysics, then the ice effective
        3531             :            ! radius has already been determined from the size distribution.
        3532           0 :            effi(i,k)   = re_ice(i,k) * 1.e6_r8      ! m -> um
        3533           0 :            effi(i,k) = effi(i,k)*micro_mg_effi_factor
        3534             :
        3535           0 :            deffi(i,k)  = effi(i,k) * 2._r8
        3536           0 :            sadice(i,k) = 4._r8*pi*(effi(i,k)**2)*ni(i,k)*rho(i,k)*1e-2_r8
        3537             :         end do
        3538             :      end do
        3539             :      !$acc end parallel
        3540             :   end if
        3541             :
        3542             :   ! cloud droplet effective radius
        3543             :   !-----------------------------------------------------------------
        3544             :
        3545             :   !$acc parallel vector_length(VLENS) default(present)
        3546             :   !$acc loop gang vector collapse(2)
        3547           0 :   do k=1,nlev
        3548           0 :      do i=1,mgncol
        3549           0 :         dum_2D(i,k) = dumnc(i,k)
        3550             :      end do
        3551             :   end do
        3552             :   !$acc end parallel
        3553             :
        3554           0 :   call size_dist_param_liq(mg_liq_props, dumc, dumnc, rho, pgam, lamc, mgncol, nlev)
        3555             :
        3556             :   !$acc parallel vector_length(VLENS) default(present)
        3557             :   !$acc loop gang vector collapse(2)
        3558           0 :   do k=1,nlev
        3559           0 :      do i=1,mgncol
        3560           0 :         if (dumc(i,k).ge.qsmall) then
        3561             :            ! switch for specification of droplet and crystal number
        3562           0 :            if (nccons) then
        3563             :               ! make sure nc is consistence with the constant N by adjusting tendency, need
        3564             :               ! to multiply by cloud fraction
        3565             :               ! note that nctend may be further adjusted below if mean droplet size is
        3566             :               ! out of bounds
        3567           0 :               nctend(i,k)=(ncnst/rho(i,k)*lcldm(i,k)-nc(i,k))*rdeltat
        3568             :            end if
        3569           0 :            if (dum_2D(i,k) /= dumnc(i,k)) then
        3570             :               ! adjust number conc if needed to keep mean size in reasonable range
        3571           0 :               nctend(i,k)=(dumnc(i,k)*lcldm(i,k)-nc(i,k))*rdeltat
        3572             :            end if
        3573             :
        3574           0 :            effc(i,k) = (pgam(i,k)+3._r8)/lamc(i,k)/2._r8*1.e6_r8
        3575             :            !assign output fields for shape here
        3576           0 :            lamcrad(i,k)=lamc(i,k)
        3577           0 :            pgamrad(i,k)=pgam(i,k)
        3578             :
        3579             :            ! recalculate effective radius for constant number, in order to separate
        3580             :            ! first and second indirect effects
        3581             :            !======================================
        3582             :            ! assume constant number of 10^8 kg-1
        3583           0 :            dumnc(i,k)=1.e8_r8
        3584             :         end if
        3585             :      end do
        3586             :   end do
        3587             :   !$acc end parallel
        3588             :
        3589             :   ! Pass in "false" adjust flag to prevent number from being changed within
        3590             :   ! size distribution subroutine.
        3591           0 :   call size_dist_param_liq(mg_liq_props, dumc, dumnc, rho, pgam, lamc, mgncol, nlev)
        3592             :
        3593             :   !$acc parallel vector_length(VLENS) default(present)
        3594             :   !$acc loop gang vector collapse(2)
        3595           0 :   do k =1,nlev
        3596           0 :      do i=1,mgncol
        3597           0 :         if (dumc(i,k).ge.qsmall) then
        3598           0 :            effc_fn(i,k) = (pgam(i,k)+3._r8)/lamc(i,k)/2._r8*1.e6_r8
        3599             :         else
        3600           0 :            effc(i,k) = 10._r8
        3601           0 :            lamcrad(i,k)=0._r8
        3602           0 :            pgamrad(i,k)=0._r8
        3603           0 :            effc_fn(i,k) = 10._r8
        3604             :         end if
        3605             :      end do
        3606             :   end do
        3607             :
        3608             :   ! recalculate 'final' rain size distribution parameters
        3609             :   ! to ensure that rain size is in bounds, adjust rain number if needed
        3610             :
        3611             :   !$acc loop gang vector collapse(2)
        3612           0 :   do k=1,nlev
        3613           0 :      do i=1,mgncol
        3614           0 :         dum_2D(i,k) = dumnr(i,k)
        3615             :      end do
        3616             :   end do
        3617             :   !$acc end parallel
        3618             :
        3619           0 :   call size_dist_param_basic(mg_rain_props, dumr, dumnr, lamr, mgncol, nlev)
        3620             :
        3621             :   !$acc parallel vector_length(VLENS) default(present)
        3622             :   !$acc loop gang vector collapse(2)
        3623           0 :   do k=1,nlev
        3624           0 :      do i=1,mgncol
        3625           0 :         if (dumr(i,k).ge.qsmall) then
        3626           0 :            if (dum_2D(i,k) /= dumnr(i,k)) then
        3627             :               ! adjust number conc if needed to keep mean size in reasonable range
        3628           0 :               nrtend(i,k)=(dumnr(i,k)*precip_frac(i,k)-nr(i,k))*rdeltat
        3629             :            end if
        3630             :
        3631             :         end if
        3632             :      end do
        3633             :   end do
        3634             :
        3635             :   ! recalculate 'final' snow size distribution parameters
        3636             :   ! to ensure that snow size is in bounds, adjust snow number if needed
        3637             :
        3638             :   !$acc loop gang vector collapse(2)
        3639           0 :   do k=1,nlev
        3640           0 :      do i=1,mgncol
        3641           0 :         dum_2D(i,k) = dumns(i,k)
        3642             :      end do
        3643             :   end do
        3644             :   !$acc end parallel
        3645             :
        3646           0 :   call size_dist_param_basic(mg_snow_props, dums, dumns, lams, mgncol, nlev, n0=dumns0A2D)
        3647             :
        3648             :   !$acc parallel vector_length(VLENS) default(present)
        3649             :   !$acc loop gang vector collapse(2)
        3650           0 :   do k=1,nlev
        3651           0 :      do i=1,mgncol
        3652           0 :         if (dums(i,k).ge.qsmall) then
        3653             :
        3654           0 :            if (dum_2D(i,k) /= dumns(i,k)) then
        3655             :               ! adjust number conc if needed to keep mean size in reasonable range
        3656           0 :               nstend(i,k)=(dumns(i,k)*precip_frac(i,k)-ns(i,k))*rdeltat
        3657             :            end if
        3658             :
        3659           0 :            sadsnow(i,k) = 2._r8*pi*(lams(i,k)**(-3))*dumns0A2D(i,k)*rho(i,k)*1.e-2_r8  ! m2/m3 -> cm2/cm3
        3660             :
        3661             :         end if
        3662             :      end do
        3663             :   end do
        3664             :
        3665             :   ! recalculate 'final' graupel size distribution parameters
        3666             :   ! to ensure that  size is in bounds, addjust number if needed
        3667             :
        3668             :   !$acc loop gang vector collapse(2)
        3669           0 :   do k=1,nlev
        3670           0 :      do i=1,mgncol
        3671           0 :         dum_2D(i,k) = dumng(i,k)
        3672             :      end do
        3673             :   end do
        3674             :   !$acc end parallel
        3675             :
        3676           0 :   if (do_hail) then
        3677           0 :      call size_dist_param_basic(mg_hail_props, dumg, dumng, lamg, mgncol, nlev)
        3678             :   end if
        3679           0 :   if (do_graupel) then
        3680           0 :      call size_dist_param_basic(mg_graupel_props, dumg, dumng, lamg, mgncol, nlev)
        3681             :   end if
        3682             :
        3683             :   !$acc parallel vector_length(VLENS) default(present)
        3684             :   !$acc loop gang vector collapse(2)
        3685           0 :   do k=1,nlev
        3686           0 :      do i=1,mgncol
        3687           0 :         if (dumg(i,k).ge.qsmall) then
        3688           0 :            if (dum_2D(i,k) /= dumng(i,k)) then
        3689             :               ! adjust number conc if needed to keep mean size in reasonable range
        3690           0 :               ngtend(i,k)=(dumng(i,k)*precip_frac(i,k)-ng(i,k))*rdeltat
        3691             :            end if
        3692             :         end if
        3693             :      end do
        3694             :   end do
        3695             :
        3696             :   !$acc loop gang vector collapse(2)
        3697           0 :   do k=1,nlev
        3698           0 :      do i=1,mgncol
        3699             :         ! if updated q (after microphysics) is zero, then ensure updated n is also zero
        3700             :         !=================================================================================
        3701           0 :         if (qc(i,k)+qctend(i,k)*deltat.lt.qsmall) nctend(i,k)=-nc(i,k)*rdeltat
        3702           0 :         if (do_cldice .and. qi(i,k)+qitend(i,k)*deltat.lt.qsmall) nitend(i,k)=-ni(i,k)*rdeltat
        3703           0 :         if (qr(i,k)+qrtend(i,k)*deltat.lt.qsmall) nrtend(i,k)=-nr(i,k)*rdeltat
        3704           0 :         if (qs(i,k)+qstend(i,k)*deltat.lt.qsmall) nstend(i,k)=-ns(i,k)*rdeltat
        3705           0 :         if (qg(i,k)+qgtend(i,k)*deltat.lt.qsmall) ngtend(i,k)=-ng(i,k)*rdeltat
        3706             :      end do
        3707             :   end do
        3708             :
        3709             :   ! DO STUFF FOR OUTPUT:
        3710             :   !==================================================
        3711             :   ! qc and qi are only used for output calculations past here,
        3712             :   ! so add qctend and qitend back in one more time
        3713             :
        3714             :   !$acc loop gang vector collapse(2)
        3715           0 :   do k=1,nlev
        3716           0 :      do i=1,mgncol
        3717           0 :         qc(i,k) = qc(i,k) + qctend(i,k)*deltat
        3718           0 :         qi(i,k) = qi(i,k) + qitend(i,k)*deltat
        3719             :      end do
        3720             :   end do
        3721             :
        3722             :   ! averaging for snow and rain number and diameter
        3723             :   !--------------------------------------------------
        3724             :
        3725             :   ! drout2/dsout2:
        3726             :   ! diameter of rain and snow
        3727             :   ! dsout:
        3728             :   ! scaled diameter of snow (passed to radiation in CAM)
        3729             :   ! reff_rain/reff_snow:
        3730             :   ! calculate effective radius of rain and snow in microns for COSP using Eq. 9 of COSP v1.3 manual
        3731             :
        3732             :   ! avoid divide by zero in avg_diameter_vec
        3733             :
        3734             :   !$acc loop gang vector collapse(2)
        3735           0 :   do k=1,nlev
        3736           0 :      do i=1,mgncol
        3737           0 :         if (nrout(i,k) .eq. 0._r8) nrout(i,k)=1.e-34_r8
        3738             :      end do
        3739             :   end do
        3740             :   !$acc end parallel
        3741             :
        3742             :   ! The avg_diameter_vec call does the actual calculation; other diameter
        3743             :   ! outputs are just drout2 times constants.
        3744           0 :   call avg_diameter_vec(qrout,nrout,rho,rhow,drout2,mgncol*nlev)
        3745             :
        3746             :   !$acc parallel vector_length(VLENS) default(present)
        3747             :   !$acc loop gang vector collapse(2)
        3748           0 :   do k=1,nlev
        3749           0 :      do i=1,mgncol
        3750           0 :         if (qrout(i,k) .gt. 1.e-7_r8 .and. nrout(i,k) .gt. 0._r8) then
        3751           0 :            qrout2(i,k) = qrout(i,k) * precip_frac(i,k)
        3752           0 :            nrout2(i,k) = nrout(i,k) * precip_frac(i,k)
        3753           0 :            freqr(i,k) = precip_frac(i,k)
        3754           0 :            reff_rain(i,k)=1.5_r8*drout2(i,k)*1.e6_r8
        3755             :         else
        3756           0 :            qrout2(i,k) = 0._r8
        3757           0 :            nrout2(i,k) = 0._r8
        3758           0 :            drout2(i,k) = 0._r8
        3759           0 :            freqr(i,k) = 0._r8
        3760           0 :            reff_rain(i,k) = 0._r8
        3761             :         end if
        3762             :      end do
        3763             :   end do
        3764             :
        3765             :   ! avoid divide by zero in avg_diameter_vec
        3766             :
        3767             :   !$acc loop gang vector collapse(2)
        3768           0 :   do k=1,nlev
        3769           0 :      do i=1,mgncol
        3770           0 :         if (nsout(i,k) .eq. 0._r8) nsout(i,k) = 1.e-34_r8
        3771             :      end do
        3772             :   end do
        3773             :   !$acc end parallel
        3774             :
        3775             :   ! The avg_diameter_vec call does the actual calculation; other diameter
        3776             :   ! outputs are just dsout2 times constants.
        3777           0 :   call avg_diameter_vec(qsout, nsout, rho, rhosn,dsout2,mgncol*nlev)
        3778             :
        3779             :   !$acc parallel vector_length(VLENS) default(present)
        3780             :   !$acc loop gang vector collapse(2)
        3781           0 :   do k=1,nlev
        3782           0 :      do i=1,mgncol
        3783           0 :         if (qsout(i,k) .gt. 1.e-7_r8 .and. nsout(i,k) .gt. 0._r8) then
        3784           0 :            qsout2(i,k) = qsout(i,k) * precip_frac(i,k)
        3785           0 :            nsout2(i,k) = nsout(i,k) * precip_frac(i,k)
        3786           0 :            freqs(i,k) = precip_frac(i,k)
        3787           0 :            dsout(i,k)=3._r8*rhosn/rhows*dsout2(i,k)
        3788           0 :            reff_snow(i,k)=1.5_r8*dsout2(i,k)*1.e6_r8
        3789             :         else
        3790           0 :            dsout(i,k)  = 0._r8
        3791           0 :            qsout2(i,k) = 0._r8
        3792           0 :            nsout2(i,k) = 0._r8
        3793           0 :            dsout2(i,k) = 0._r8
        3794           0 :            freqs(i,k)  = 0._r8
        3795           0 :            reff_snow(i,k)=0._r8
        3796             :         end if
        3797             :      end do
        3798             :   end do
        3799             :
        3800             :   ! avoid divide by zero in avg_diameter_vec
        3801             :
        3802             :   !$acc loop gang vector collapse(2)
        3803           0 :   do k=1,nlev
        3804           0 :      do i=1,mgncol
        3805           0 :         if (ngout(i,k) .eq. 0._r8) ngout(i,k) = 1.e-34_r8
        3806             :      end do
        3807             :   end do
        3808             :   !$acc end parallel
        3809             :
        3810             :   ! The avg_diameter_vec call does the actual calculation; other diameter
        3811             :   ! outputs are just dgout2 times constants.
        3812           0 :   if (do_hail .or. do_graupel) then
        3813           0 :      call avg_diameter_vec(qgout, ngout, rho, rhogtmp, dgout2, mgncol*nlev)
        3814             :   else
        3815             :      ! need this if statement for MG2, where rhogtmp = 0
        3816             :
        3817             :      !$acc parallel vector_length(VLENS) default(present)
        3818             :      !$acc loop gang vector collapse(2)
        3819           0 :      do k=1,nlev
        3820           0 :         do i=1,mgncol
        3821           0 :            dgout2(i,k) = 0._r8
        3822             :         end do
        3823             :      end do
        3824             :      !$acc end parallel
        3825             :   end if
        3826             :
        3827             :   !$acc parallel vector_length(VLENS) default(present)
        3828             :   !$acc loop gang vector collapse(2)
        3829           0 :   do k=1,nlev
        3830           0 :      do i=1,mgncol
        3831           0 :         if (qgout(i,k) .gt. 1.e-7_r8 .and. ngout(i,k) .gt. 0._r8) then
        3832           0 :            qgout2(i,k) = qgout(i,k) * precip_frac(i,k)
        3833           0 :            ngout2(i,k) = ngout(i,k) * precip_frac(i,k)
        3834           0 :            freqg(i,k) = precip_frac(i,k)
        3835           0 :            dgout(i,k)=3._r8*rhogtmp/rhows*dgout2(i,k)
        3836           0 :            reff_grau(i,k)=1.5_r8*dgout2(i,k)*1.e6_r8
        3837             :         else
        3838           0 :            dgout(i,k)  = 0._r8
        3839           0 :            qgout2(i,k) = 0._r8
        3840           0 :            ngout2(i,k) = 0._r8
        3841           0 :            dgout2(i,k) = 0._r8
        3842           0 :            freqg(i,k)  = 0._r8
        3843           0 :            reff_grau(i,k)=0._r8
        3844             :         end if
        3845             :      end do
        3846             :   end do
        3847             :
        3848             :   ! analytic radar reflectivity
        3849             :   !--------------------------------------------------
        3850             :   ! formulas from Matthew Shupe, NOAA/CERES
        3851             :   ! *****note: radar reflectivity is local (in-precip average)
        3852             :   ! units of mm^6/m^3
        3853             :
        3854             :   !$acc loop gang vector collapse(2) private(dum,dum1)
        3855           0 :   do k=1,nlev
        3856           0 :      do i=1,mgncol
        3857           0 :         if (qc(i,k).ge.qsmall .and. (nc(i,k)+nctend(i,k)*deltat).gt.10._r8) then
        3858             :            dum=(qc(i,k)/lcldm(i,k)*rho(i,k)*1000._r8)**2 &
        3859           0 :                 /(0.109_r8*(nc(i,k)+nctend(i,k)*deltat)/lcldm(i,k)*rho(i,k)/1.e6_r8)*lcldm(i,k)/precip_frac(i,k)
        3860             :         else
        3861             :            dum=0._r8
        3862             :         end if
        3863           0 :         if (qi(i,k).ge.qsmall) then
        3864           0 :            dum1=(qi(i,k)*rho(i,k)/icldm(i,k)*1000._r8/0.1_r8)**(1._r8/0.63_r8)*icldm(i,k)/precip_frac(i,k)
        3865             :         else
        3866           0 :            dum1=0._r8
        3867             :         end if
        3868           0 :         if (qsout(i,k).ge.qsmall) then
        3869           0 :            dum1=dum1+(qsout(i,k)*rho(i,k)*1000._r8/0.1_r8)**(1._r8/0.63_r8)
        3870             :         end if
        3871           0 :         refl(i,k)=dum+dum1
        3872             :
        3873             :         ! add rain rate, but for 37 GHz formulation instead of 94 GHz
        3874             :         ! formula approximated from data of Matrasov (2007)
        3875             :         ! rainrt is the rain rate in mm/hr
        3876             :         ! reflectivity (dum) is in DBz
        3877           0 :         if (rainrt(i,k).ge.0.001_r8) then
        3878           0 :            dum=log10(rainrt(i,k)**6._r8)+16._r8
        3879             :            ! convert from DBz to mm^6/m^3
        3880           0 :            dum = 10._r8**(dum/10._r8)
        3881             :         else
        3882             :            ! don't include rain rate in R calculation for values less than 0.001 mm/hr
        3883             :            dum=0._r8
        3884             :         end if
        3885             :
        3886             :         ! add to refl
        3887           0 :         refl(i,k)=refl(i,k)+dum
        3888             :
        3889             :         !output reflectivity in Z.
        3890           0 :         areflz(i,k)=refl(i,k) * precip_frac(i,k)
        3891             :
        3892             :         ! convert back to DBz
        3893           0 :         if (refl(i,k).gt.minrefl) then
        3894           0 :            refl(i,k)=10._r8*log10(refl(i,k))
        3895             :         else
        3896           0 :            refl(i,k)=-9999._r8
        3897             :         end if
        3898             :
        3899             :         !set averaging flag
        3900           0 :         if (refl(i,k).gt.mindbz) then
        3901           0 :            arefl(i,k)=refl(i,k) * precip_frac(i,k)
        3902           0 :            frefl(i,k)=precip_frac(i,k)
        3903             :         else
        3904           0 :            arefl(i,k)=0._r8
        3905           0 :            areflz(i,k)=0._r8
        3906           0 :            frefl(i,k)=0._r8
        3907             :         end if
        3908             :
        3909             :         ! bound cloudsat reflectivity
        3910           0 :         csrfl(i,k)=min(csmax,refl(i,k))
        3911             :
        3912             :         !set averaging flag
        3913           0 :         if (csrfl(i,k).gt.csmin) then
        3914           0 :            acsrfl(i,k)=refl(i,k) * precip_frac(i,k)
        3915           0 :            fcsrfl(i,k)=precip_frac(i,k)
        3916             :         else
        3917           0 :            acsrfl(i,k)=0._r8
        3918           0 :            fcsrfl(i,k)=0._r8
        3919             :         end if
        3920             :      end do
        3921             :   end do
        3922             :
        3923             :   !redefine fice here....
        3924             :
        3925             :   !$acc loop gang vector collapse(2)
        3926           0 :   do k=1,nlev
        3927           0 :      do i=1,mgncol
        3928           0 :         dum_2D(i,k) = qsout(i,k) + qrout(i,k) + qc(i,k) + qi(i,k)
        3929           0 :         dumi(i,k)   = qsout(i,k) + qi(i,k)
        3930           0 :         if (dumi(i,k) .gt. qsmall .and. dum_2D(i,k) .gt. qsmall) then
        3931           0 :            nfice(i,k) = min(dumi(i,k)/dum_2D(i,k),1._r8)
        3932             :         else
        3933           0 :            nfice(i,k) = 0._r8
        3934             :         end if
        3935             :      end do
        3936             :   end do
        3937             :   !$acc end parallel
        3938             :
        3939             :   !$acc end data
        3940             :
        3941           0 : end subroutine micro_pumas_tend
        3942             :
        3943             : !========================================================================
        3944             : !OUTPUT CALCULATIONS
        3945             : !========================================================================
        3946             :
        3947           0 : subroutine calc_rercld(lamr, n0r, lamc, pgam, qric, qcic, ncic, rercld, vlen)
        3948             :   integer,                   intent(in) :: vlen
        3949             :   real(r8), dimension(vlen), intent(in) :: lamr          ! rain size parameter (slope)
        3950             :   real(r8), dimension(vlen), intent(in) :: n0r           ! rain size parameter (intercept)
        3951             :   real(r8), dimension(vlen), intent(in) :: lamc          ! size distribution parameter (slope)
        3952             :   real(r8), dimension(vlen), intent(in) :: pgam          ! droplet size parameter
        3953             :   real(r8), dimension(vlen), intent(in) :: qric          ! in-cloud rain mass mixing ratio
        3954             :   real(r8), dimension(vlen), intent(in) :: qcic          ! in-cloud cloud liquid
        3955             :   real(r8), dimension(vlen), intent(in) :: ncic          ! in-cloud droplet number concentration
        3956             :
        3957             :   real(r8), dimension(vlen), intent(inout) :: rercld     ! effective radius calculation for rain + cloud
        3958             :
        3959             :   ! combined size of precip & cloud drops
        3960           0 :   real(r8) :: Atmp,tmp(vlen), pgamp1(vlen)
        3961             :
        3962             :   integer :: i
        3963             :
        3964             :   !$acc data present (rercld,lamr,n0r,lamc,pgam,qric,qcic,ncic) &
        3965             :   !$acc      create  (Atmp,tmp,pgamp1)
        3966             :
        3967             :   !$acc parallel vector_length(VLENS) default(present)
        3968             :   !$acc loop gang vector
        3969           0 :   do i=1,vlen
        3970           0 :      pgamp1(i) = pgam(i)+1._r8
        3971             :   end do
        3972             :   !$acc end parallel
        3973             :
        3974           0 :   call rising_factorial(pgamp1, 2, tmp, vlen)
        3975             :
        3976             :   !$acc parallel vector_length(VLENS) default(present)
        3977             :   !$acc loop gang vector private(Atmp)
        3978           0 :   do i=1,vlen
        3979             :      ! Rain drops
        3980           0 :      if (lamr(i) > 0._r8) then
        3981           0 :         Atmp = n0r(i) * pi / (2._r8 * lamr(i)**3._r8)
        3982             :      else
        3983             :         Atmp = 0._r8
        3984             :      end if
        3985             :      ! Add cloud drops
        3986           0 :      if (lamc(i) > 0._r8) then
        3987             :         Atmp = Atmp + &
        3988           0 :              ncic(i) * pi * tmp(i) / (4._r8 * lamc(i)**2._r8)
        3989             :      end if
        3990           0 :      if (Atmp > 0._r8) then
        3991           0 :         rercld(i) = rercld(i) + 3._r8 *(qric(i) + qcic(i)) / (4._r8 * rhow * Atmp)
        3992             :      end if
        3993             :   end do
        3994             :   !$acc end parallel
        3995             :
        3996             :   !$acc end data
        3997           0 : end subroutine calc_rercld
        3998             :
        3999             : !========================================================================
        4000             : !2020-09-15: Follow John Dennis's version to generate a new interface
        4001             : !            to update tendency in the sedimentation loop
        4002             : !========================================================================
        4003           0 : subroutine Sedimentation(mgncol,nlev,do_cldice,deltat,fx,fnx,pdel_inv,qxtend,nxtend, &
        4004           0 :                             qxsedten,dumx,dumnx,prect,xflx,xxlx,qxsevap,xcldm,tlat,qvlat,preci)
        4005             :
        4006             :    integer, intent(in)               :: mgncol,nlev
        4007             :    logical, intent(in)               :: do_cldice
        4008             :    real(r8),intent(in)               :: deltat
        4009             :    real(r8), intent(in)              :: fx(mgncol,nlev)
        4010             :    real(r8), intent(in)              :: fnx(mgncol,nlev)
        4011             :    real(r8), intent(in)              :: pdel_inv(mgncol,nlev)
        4012             :    real(r8), intent(inout)           :: qxtend(mgncol,nlev)
        4013             :    real(r8), intent(inout)           :: nxtend(mgncol,nlev)
        4014             :    real(r8), intent(inout)           :: qxsedten(mgncol,nlev)
        4015             :    real(r8), intent(inout)           :: dumx(mgncol,nlev)
        4016             :    real(r8), intent(inout)           :: dumnx(mgncol,nlev)
        4017             :    real(r8), intent(inout)           :: prect(mgncol)
        4018             :    real(r8), intent(inout)           :: xflx(mgncol,nlev+1)
        4019             :    real(r8), intent(in)   , optional :: xxlx
        4020             :    real(r8), intent(inout), optional :: qxsevap(mgncol,nlev)
        4021             :    real(r8), intent(in)   , optional :: xcldm(mgncol,nlev)
        4022             :    real(r8), intent(inout), optional :: tlat(mgncol,nlev)
        4023             :    real(r8), intent(inout), optional :: qvlat(mgncol,nlev)
        4024             :    real(r8), intent(inout), optional :: preci(mgncol)
        4025             :    integer  :: i,k,n,nstep
        4026             :    real(r8) :: faltndx,faltndnx,rnstep,faltndqxe
        4027           0 :    real(r8) :: dum1(mgncol,nlev),faloutx(mgncol,0:nlev),faloutnx(mgncol,0:nlev)
        4028             :    logical  :: present_tlat,present_qvlat, present_xcldm,present_qxsevap, present_preci
        4029             :
        4030           0 :    present_tlat    = present(tlat)
        4031           0 :    present_qvlat   = present(qvlat)
        4032           0 :    present_xcldm   = present(xcldm)
        4033           0 :    present_qxsevap = present(qxsevap)
        4034           0 :    present_preci   = present(preci)
        4035             :
        4036             :    ! loop over sedimentation sub-time step to ensure stability
        4037             :    !==============================================================
        4038             :
        4039             :    !$acc data present (fx,fnx,pdel_inv,qxtend,nxtend,qxsedten,dumx,dumnx) &
        4040             :    !$acc      present (prect,xflx,xxlx,qxsevap,xcldm,tlat,qvlat,preci) &
        4041             :    !$acc      create  (faloutx,faloutnx,dum1)
        4042             :
        4043             :    !$acc parallel vector_length(VLENS) default(present)
        4044             :    !$acc loop gang vector private(faltndx,faltndnx,faltndqxe,nstep,rnstep)
        4045           0 :    do i = 1,mgncol
        4046           0 :       nstep   = 1 + int( max( maxval( fx(i,:)*pdel_inv(i,:) ), &
        4047             :                               maxval( fnx(i,:)*pdel_inv(i,:) ) ) &
        4048           0 :                               * deltat )
        4049           0 :       rnstep  = 1._r8/real(nstep)
        4050             :
        4051           0 :       dum1(i,1) = 0._r8
        4052           0 :       if (present_xcldm) then
        4053             :          !$acc loop vector
        4054           0 :          do k = 2,nlev
        4055           0 :             dum1(i,k) = xcldm(i,k)/xcldm(i,k-1)
        4056           0 :             dum1(i,k) = min(dum1(i,k),1._r8)
        4057             :          end do
        4058             :       else
        4059             :          !$acc loop vector
        4060           0 :          do k=2,nlev
        4061           0 :             dum1(i,k) = 1._r8
        4062             :          end do
        4063             :       end if
        4064             :
        4065             :       !$acc loop seq
        4066           0 :       do n = 1,nstep
        4067           0 :          faloutx(i,0)  = 0._r8
        4068           0 :          faloutnx(i,0) = 0._r8
        4069           0 :          if (do_cldice) then
        4070             :             !$acc loop vector
        4071           0 :             do k=1,nlev
        4072           0 :                faloutx(i,k)  = fx(i,k)  * dumx(i,k)
        4073           0 :                faloutnx(i,k) = fnx(i,k) * dumnx(i,k)
        4074             :             end do
        4075             :          else
        4076             :             !$acc loop vector
        4077           0 :             do k=1,nlev
        4078           0 :                faloutx(i,k)  = 0._r8
        4079           0 :                faloutnx(i,k) = 0._r8
        4080             :             end do
        4081             :          end if
        4082             :
        4083             :          !$acc loop vector
        4084           0 :          do k = 1,nlev
        4085             :             ! for cloud liquid and ice, if cloud fraction increases with height
        4086             :             ! then add flux from above to both vapor and cloud water of current level
        4087             :             ! this means that flux entering clear portion of cell from above evaporates
        4088             :             ! instantly
        4089             :             ! note: this is not an issue with precip, since we assume max overlap
        4090           0 :             faltndx       = (faloutx(i,k)  - dum1(i,k)*faloutx(i,k-1))*pdel_inv(i,k)
        4091           0 :             faltndnx      = (faloutnx(i,k) - dum1(i,k)*faloutnx(i,k-1))*pdel_inv(i,k)
        4092             :             ! add fallout terms to eulerian tendencies
        4093           0 :             qxtend(i,k)   = qxtend(i,k)-faltndx*rnstep
        4094           0 :             nxtend(i,k)   = nxtend(i,k)-faltndnx*rnstep
        4095             :             ! sedimentation tendency for output
        4096           0 :             qxsedten(i,k) = qxsedten(i,k)-faltndx*rnstep
        4097             :             ! add terms to to evap/sub of cloud water
        4098           0 :             dumx(i,k)     = dumx(i,k)  - faltndx*deltat*rnstep
        4099           0 :             dumnx(i,k)    = dumnx(i,k) - faltndnx*deltat*rnstep
        4100             :
        4101           0 :             if (k>1) then
        4102           0 :                faltndqxe     = (faloutx(i,k)-faloutx(i,k-1))*pdel_inv(i,k)
        4103             :                ! for output
        4104           0 :                if(present_qxsevap) qxsevap(i,k)= qxsevap(i,k) - (faltndqxe-faltndx)*rnstep
        4105           0 :                if(present_qvlat)   qvlat(i,k)  = qvlat(i,k) - (faltndqxe-faltndx)*rnstep
        4106           0 :                if(present_tlat)    tlat(i,k)   = tlat(i,k) + (faltndqxe-faltndx)*xxlx*rnstep
        4107             :             end if
        4108             :
        4109           0 :             xflx(i,k+1) = xflx(i,k+1) + faloutx(i,k) / g * rnstep
        4110             :          end do
        4111             :
        4112             :          ! units below are m/s
        4113             :          ! sedimentation flux at surface is added to precip flux at surface
        4114             :          ! to get total precip (cloud + precip water) rate
        4115           0 :          prect(i) = prect(i)+faloutx(i,nlev)/g*rnstep/1000._r8
        4116           0 :          if(present_preci) preci(i) = preci(i)+faloutx(i,nlev)/g*rnstep/1000._r8
        4117             :       end do  ! n loop of 1, nstep
        4118             :    end do  ! i loop of 1, mgncol
        4119             :    !$acc end parallel
        4120             :
        4121             :    !$acc end data
        4122           0 : end subroutine Sedimentation
        4123             :
        4124             : !========================================================================
        4125             : !UTILITIES
        4126             : !========================================================================
        4127             :
        4128             :
        4129           0 : pure subroutine micro_pumas_get_cols(ncol, nlev, top_lev, mgncol, mgcols, &
        4130           0 :      qcn, qin, qrn, qsn, qgr)
        4131             :
        4132             :   ! Determines which columns microphysics should operate over by
        4133             :   ! checking for non-zero cloud water/ice.
        4134             :
        4135             :   integer, intent(in) :: ncol      ! Number of columns with meaningful data
        4136             :   integer, intent(in) :: nlev      ! Number of levels to use
        4137             :   integer, intent(in) :: top_lev   ! Top level for microphysics
        4138             :   integer, intent(out) :: mgncol   ! Number of columns MG will use
        4139             :   integer, allocatable, intent(out) :: mgcols(:) ! column indices
        4140             :
        4141             :   real(r8), intent(in) :: qcn(:,:) ! cloud water mixing ratio (kg/kg)
        4142             :   real(r8), intent(in) :: qin(:,:) ! cloud ice mixing ratio (kg/kg)
        4143             :   real(r8), intent(in) :: qrn(:,:) ! rain mixing ratio (kg/kg)
        4144             :   real(r8), intent(in) :: qsn(:,:) ! snow mixing ratio (kg/kg)
        4145             :   real(r8), optional, intent(in) :: qgr(:,:) ! graupel mixing ratio (kg/kg)
        4146             :
        4147             :   integer :: lev_offset  ! top_lev - 1 (defined here for consistency)
        4148           0 :   logical :: ltrue(ncol) ! store tests for each column
        4149             :
        4150             :   integer :: i, ii ! column indices
        4151             :
        4152           0 :   if (allocated(mgcols)) deallocate(mgcols)
        4153             :
        4154           0 :   lev_offset = top_lev - 1
        4155             :
        4156             :   ! Using "any" along dimension 2 collapses across levels, but
        4157             :   ! not columns, so we know if water is present at any level
        4158             :   ! in each column.
        4159             :
        4160           0 :   ltrue = any(qcn(:ncol,top_lev:(nlev+lev_offset)) >= qsmall, 2)
        4161           0 :   ltrue = ltrue .or. any(qin(:ncol,top_lev:(nlev+lev_offset)) >= qsmall, 2)
        4162           0 :   ltrue = ltrue .or. any(qrn(:ncol,top_lev:(nlev+lev_offset)) >= qsmall, 2)
        4163           0 :   ltrue = ltrue .or. any(qsn(:ncol,top_lev:(nlev+lev_offset)) >= qsmall, 2)
        4164             :
        4165           0 :   if(present(qgr)) ltrue = ltrue .or. any(qgr(:ncol,top_lev:(nlev+lev_offset)) >= qsmall, 2)
        4166             :
        4167             :   ! Scan for true values to get a usable list of indices.
        4168             :
        4169           0 :   mgncol = count(ltrue)
        4170           0 :   allocate(mgcols(mgncol))
        4171           0 :   i = 0
        4172           0 :   do ii = 1,ncol
        4173           0 :      if (ltrue(ii)) then
        4174           0 :         i = i + 1
        4175           0 :         mgcols(i) = ii
        4176             :      end if
        4177             :   end do
        4178             :
        4179           0 : end subroutine micro_pumas_get_cols
        4180             :
        4181             : end module micro_pumas_v1

Generated by: LCOV version 1.14