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

              Line data    Source code
           1             : module micro_pumas_utils
           2             :
           3             : !--------------------------------------------------------------------------
           4             : !
           5             : ! This module contains process rates and utility functions used by the MG
           6             : ! microphysics.
           7             : !
           8             : ! Original MG authors: Andrew Gettelman, Hugh Morrison
           9             : ! Contributions from: Peter Caldwell, Xiaohong Liu and Steve Ghan
          10             : !
          11             : ! Separated from MG 1.5 by B. Eaton.
          12             : ! Separated module switched to MG 2.0 and further changes by S. Santos.
          13             : !
          14             : ! for questions contact Hugh Morrison, Andrew Gettelman
          15             : ! e-mail: morrison@ucar.edu, andrew@ucar.edu
          16             : !
          17             : !--------------------------------------------------------------------------
          18             : !
          19             : ! List of required external functions that must be supplied:
          20             : !   gamma --> standard mathematical gamma function (if gamma is an
          21             : !       intrinsic, define HAVE_GAMMA_INTRINSICS)
          22             : !
          23             : !--------------------------------------------------------------------------
          24             : !
          25             : ! Constants that must be specified in the "init" method (module variables):
          26             : !
          27             : ! kind            kind of reals (to verify correct linkage only) -
          28             : ! gravit          acceleration due to gravity                    m s-2
          29             : ! rair            dry air gas constant for air                   J kg-1 K-1
          30             : ! rh2o            gas constant for water vapor                   J kg-1 K-1
          31             : ! cpair           specific heat at constant pressure for dry air J kg-1 K-1
          32             : ! tmelt           temperature of melting point for water         K
          33             : ! latvap          latent heat of vaporization                    J kg-1
          34             : ! latice          latent heat of fusion                          J kg-1
          35             : !
          36             : !--------------------------------------------------------------------------
          37             :
          38             : #ifndef HAVE_GAMMA_INTRINSICS
          39             : use shr_spfn_mod, only: gamma => shr_spfn_gamma
          40             : #endif
          41             :
          42             : implicit none
          43             : private
          44             : save
          45             :
          46             : public :: &
          47             :      micro_pumas_utils_init, &
          48             :      size_dist_param_liq, &
          49             :      size_dist_param_basic, &
          50             :      avg_diameter, &
          51             :      rising_factorial, &
          52             :      ice_deposition_sublimation, &
          53             :      ice_deposition_sublimation_mg4, &  !! ktc
          54             :      sb2001v2_liq_autoconversion,&
          55             :      sb2001v2_accre_cld_water_rain,&
          56             :      kk2000_liq_autoconversion, &
          57             :      ice_autoconversion, &
          58             :      immersion_freezing, &
          59             :      contact_freezing, &
          60             :      snow_self_aggregation, &
          61             :      ice_self_aggregation, &  !! mg4 added for using instead of snow_self_agg
          62             :      accrete_cloud_water_snow, &
          63             :      accrete_cloud_water_ice, &  !! mg4 added for using instead of cloud_water_snow
          64             :      secondary_ice_production, &
          65             :      accrete_rain_snow, &
          66             :      accrete_rain_ice, &  !! mg4 added for using instead of rain_snow
          67             :      heterogeneous_rain_freezing, &
          68             :      accrete_cloud_water_rain, &
          69             :      self_collection_rain, &
          70             :      accrete_cloud_ice_snow, &
          71             :      evaporate_sublimate_precip, &
          72             :      bergeron_process_snow, &
          73             :      graupel_collecting_snow, &
          74             :      graupel_collecting_rain, &
          75             :      graupel_collecting_cld_water, &
          76             :      graupel_riming_liquid_snow, &
          77             :      graupel_rain_riming_snow, &
          78             :      graupel_rime_splintering, &
          79             :      evaporate_sublimate_precip_graupel, &
          80             :      evaporate_sublimate_precip_mg4, &
          81             :      evaporate_sublimate_precip_graupel_mg4, &
          82             :      access_lookup_table, & !! mg4
          83             :      access_lookup_table_coll, & !! mg4
          84             :      init_lookup_table, &      !! mg4
          85             :      avg_diameter_vec
          86             :
          87             : ! 8 byte real and integer
          88             : integer, parameter, public :: r8 = selected_real_kind(12)
          89             : integer, parameter, public :: i8 = selected_int_kind(18)
          90             : integer, parameter         :: VLENS = 128  ! vector length of a GPU compute kernel
          91             :
          92             : public :: MGHydrometeorProps
          93             :
          94             : type :: MGHydrometeorProps
          95             :    ! Density (kg/m^3)
          96             :    real(r8) :: rho
          97             :    ! Information for size calculations.
          98             :    ! Basic calculation of mean size is:
          99             :    !     lambda = (shape_coef*nic/qic)^(1/eff_dim)
         100             :    ! Then lambda is constrained by bounds.
         101             :    real(r8) :: eff_dim
         102             :    real(r8) :: shape_coef
         103             :    real(r8) :: lambda_bounds(2)
         104             :    ! Minimum average particle mass (kg).
         105             :    ! Limit is applied at the beginning of the size distribution calculations.
         106             :    real(r8) :: min_mean_mass
         107             : end type MGHydrometeorProps
         108             :
         109             : interface MGHydrometeorProps
         110             :    module procedure NewMGHydrometeorProps
         111             : end interface
         112             :
         113             : type(MGHydrometeorProps), public :: mg_liq_props
         114             : type(MGHydrometeorProps), public :: mg_ice_props
         115             : type(MGHydrometeorProps), public :: mg_rain_props
         116             : type(MGHydrometeorProps), public :: mg_snow_props
         117             : type(MGHydrometeorProps), public :: mg_graupel_props
         118             : type(MGHydrometeorProps), public :: mg_hail_props
         119             :
         120             : interface size_dist_param_liq
         121             :   module procedure size_dist_param_liq_2D
         122             :   module procedure size_dist_param_liq_vect
         123             :   module procedure size_dist_param_liq_line
         124             : end interface
         125             : interface size_dist_param_basic
         126             :   module procedure size_dist_param_basic_2D
         127             :   module procedure size_dist_param_basic_vect
         128             :   module procedure size_dist_param_basic_line
         129             : end interface
         130             : interface calc_ab
         131             :   module procedure calc_ab_line
         132             :   module procedure calc_ab_vect
         133             : end interface
         134             :
         135             : !=================================================
         136             : ! Public module parameters (mostly for MG itself)
         137             : !=================================================
         138             :
         139             : ! Pi to 20 digits; more than enough to reach the limit of double precision.
         140             : real(r8), parameter, public :: pi = 3.14159265358979323846_r8
         141             :
         142             : ! "One minus small number": number near unity for round-off issues.
         143             : real(r8), parameter, public :: omsm   = 1._r8 - 1.e-5_r8
         144             :
         145             : ! Smallest mixing ratio considered in microphysics.
         146             : real(r8), parameter, public :: qsmall = 1.e-18_r8
         147             :
         148             : ! minimum allowed cloud fraction
         149             : real(r8), parameter, public :: mincld = 0.0001_r8
         150             :
         151             : real(r8), parameter, public :: rhosn = 250._r8  ! bulk density snow
         152             : real(r8), parameter, public :: rhoi = 500._r8   ! bulk density ice
         153             : real(r8), parameter, public :: rhow = 1000._r8  ! bulk density liquid
         154             : real(r8), parameter, public :: rhows = 917._r8  ! bulk density water solid
         155             :
         156             : !Hail and Graupel (set in MG3)
         157             : real(r8), parameter, public :: rhog = 500._r8
         158             : real(r8), parameter, public :: rhoh = 900._r8
         159             :
         160             : ! fall speed parameters, V = aD^b (V is in m/s)
         161             : ! droplets
         162             : real(r8), parameter, public :: ac = 3.e7_r8
         163             : real(r8), parameter, public :: bc = 2._r8
         164             : ! snow
         165             : real(r8), parameter, public :: as = 11.72_r8
         166             : real(r8), parameter, public :: bs = 0.41_r8
         167             : ! cloud ice
         168             : real(r8), parameter, public :: ai = 700._r8
         169             : real(r8), parameter, public :: bi = 1._r8
         170             : ! small cloud ice (r< 10 um) - sphere, bulk density
         171             : real(r8), parameter, public :: aj = ac*((rhoi/rhows)**(bc/3._r8))*rhows/rhow
         172             : real(r8), parameter, public :: bj = bc
         173             : ! rain
         174             : real(r8), parameter, public :: ar = 841.99667_r8
         175             : real(r8), parameter, public :: br = 0.8_r8
         176             : ! graupel
         177             : real(r8), parameter, public :: ag = 19.3_r8
         178             : real(r8), parameter, public :: bg = 0.37_r8
         179             : ! hail
         180             : real(r8), parameter, public :: ah = 114.5_r8
         181             : real(r8), parameter, public :: bh = 0.5_r8
         182             :
         183             : ! mass of new crystal due to aerosol freezing and growth (kg)
         184             : ! Make this consistent with the lower bound, to support UTLS and
         185             : ! stratospheric ice, and the smaller ice size limit.
         186             : real(r8), parameter, public :: mi0 = 4._r8/3._r8*pi*rhoi*(1.e-6_r8)**3
         187             :
         188             : ! mass of new graupel particle  (assume same as mi0 for now, may want to make bigger?)
         189             : !real(r8), parameter, public :: mg0 = 4._r8/3._r8*pi*rhoi*(1.e-6_r8)**3
         190             : !or set based on M2005:
         191             : real(r8), parameter, public :: mg0 = 1.6e-10_r8
         192             : ! radius of contact nuclei
         193             : real(r8), parameter, public :: mmult = 4._r8/3._r8*pi*rhoi*(5.e-6_r8)**3
         194             :
         195             : !=================================================
         196             : ! Private module parameters
         197             : !=================================================
         198             :
         199             : ! Signaling NaN bit pattern that represents a limiter that's turned off.
         200             : integer(i8), parameter :: limiter_off = int(Z'7FF1111111111111', i8)
         201             :
         202             : ! alternate threshold used for some in-cloud mmr
         203             : real(r8), parameter, public :: icsmall = 1.e-8_r8
         204             :
         205             : ! particle mass-diameter relationship
         206             : ! currently we assume spherical particles for cloud ice/snow
         207             : ! m = cD^d
         208             : ! exponent
         209             : real(r8), parameter :: dsph = 3._r8
         210             :
         211             : ! Bounds for mean diameter for different constituents.
         212             : real(r8), parameter :: lam_bnd_rain(2) = 1._r8/[500.e-6_r8, 20.e-6_r8]
         213             : real(r8), parameter :: lam_bnd_snow(2) = 1._r8/[2000.e-6_r8, 10.e-6_r8]
         214             :
         215             : ! Minimum average mass of particles.
         216             : real(r8), parameter :: min_mean_mass_liq = 1.e-20_r8
         217             : real(r8), parameter :: min_mean_mass_ice = 1.e-20_r8
         218             :
         219             : ! ventilation parameters
         220             : ! for snow
         221             : real(r8), parameter :: f1s = 0.86_r8
         222             : real(r8), parameter :: f2s = 0.28_r8
         223             : ! for rain
         224             : real(r8), parameter :: f1r = 0.78_r8
         225             : real(r8), parameter :: f2r = 0.308_r8
         226             :
         227             : ! collection efficiencies
         228             : ! aggregation of cloud ice and snow
         229             : real(r8), parameter :: eii = 0.5_r8
         230             : ! collection efficiency, ice-droplet collisions
         231             : real(r8), parameter, public :: ecid = 0.7_r8
         232             : ! collection efficiency between droplets/rain and snow/rain
         233             : real(r8), parameter, public :: ecr = 1.0_r8
         234             :
         235             : ! immersion freezing parameters, bigg 1953
         236             : real(r8), parameter :: bimm = 100._r8
         237             : real(r8), parameter :: aimm = 0.66_r8
         238             :
         239             : ! Mass of each raindrop created from autoconversion.
         240             : real(r8), parameter :: droplet_mass_25um = 4._r8/3._r8*pi*rhow*(25.e-6_r8)**3
         241             : real(r8), parameter :: droplet_mass_40um = 4._r8/3._r8*pi*rhow*(40.e-6_r8)**3
         242             :
         243             : !!!!!!
         244             : ! MG4 look up table parameters
         245             : !!!!!!
         246             : ! ice microphysics lookup table array dimensions
         247             : integer, parameter, public :: tsize     = 3
         248             : integer, parameter, public :: isize     = 20
         249             : integer, parameter, public :: jsize     = 20
         250             : integer, parameter, public :: rcollsize = 30
         251             :
         252             : ! number of ice microphysical quantities used from lookup table
         253             : integer, parameter, public :: tabsize   = 14
         254             :
         255             : ! number of ice-rain collection microphysical quantities used from lookup table
         256             : integer, parameter, public :: colltabsize = 2
         257             :
         258             :
         259             : !ice lookup table values
         260             : real(r8) :: itab(tsize,isize,jsize,tabsize)
         261             :
         262             : !ice lookup table values for ice-rain collision/collection
         263             : real(r8) :: itabcoll(tsize,isize,jsize,rcollsize,colltabsize)
         264             :
         265             : private :: itab,itabcoll
         266             :
         267             : !=========================================================
         268             : ! Constants set in initialization
         269             : !=========================================================
         270             :
         271             : ! Set using arguments to micro_pumas_init
         272             : real(r8) :: rv          ! water vapor gas constant
         273             : real(r8) :: cpp         ! specific heat of dry air
         274             : real(r8) :: tmelt       ! freezing point of water (K)
         275             :
         276             : real(r8) :: ra        ! dry air gas constant
         277             :
         278             : ! latent heats of:
         279             : real(r8) :: xxlv        ! vaporization
         280             : real(r8) :: xlf         ! freezing
         281             : real(r8) :: xxls        ! sublimation
         282             :
         283             : ! additional constants to help speed up code
         284             : real(r8) :: gamma_bs_plus3
         285             : real(r8) :: gamma_half_br_plus5
         286             : real(r8) :: gamma_half_bs_plus5
         287             : real(r8) :: gamma_2bs_plus2
         288             :
         289             : !$acc declare create (rv,cpp,tmelt,xxlv,xxls,gamma_bs_plus3,   &
         290             : !$acc                 gamma_half_br_plus5,gamma_half_bs_plus5, &
         291             : !$acc                 gamma_2bs_plus2)
         292             :
         293             : !=========================================================
         294             : ! Utilities that are cheaper if the compiler knows that
         295             : ! some argument is an integer.
         296             : !=========================================================
         297             :
         298             : interface rising_factorial
         299             :    module procedure rising_factorial_r8
         300             :    module procedure rising_factorial_r8_vec
         301             :    module procedure rising_factorial_integer
         302             :    module procedure rising_factorial_integer_vec
         303             : end interface rising_factorial
         304             :
         305             : interface var_coef
         306             :    module procedure var_coef_r8
         307             :    module procedure var_coef_r8_vect
         308             :    module procedure var_coef_integer
         309             :    module procedure var_coef_integer_vect
         310             : end interface var_coef
         311             :
         312             : !==========================================================================
         313             : contains
         314             : !==========================================================================
         315             :
         316             : ! Initialize module variables.
         317             : !
         318             : ! "kind" serves no purpose here except to check for unlikely linking
         319             : ! issues; always pass in the kind for a double precision real.
         320             : !
         321             : ! "errstring" is the only output; it is blank if there is no error, or set
         322             : ! to a message if there is an error.
         323             : !
         324             : ! Check the list at the top of this module for descriptions of all other
         325             : ! arguments.
         326           0 : subroutine micro_pumas_utils_init( kind, rair, rh2o, cpair, tmelt_in, latvap, &
         327           0 :      latice, dcs, errstring)
         328             :
         329             :   integer,  intent(in)  :: kind
         330             :   real(r8), intent(in)  :: rair
         331             :   real(r8), intent(in)  :: rh2o
         332             :   real(r8), intent(in)  :: cpair
         333             :   real(r8), intent(in)  :: tmelt_in
         334             :   real(r8), intent(in)  :: latvap
         335             :   real(r8), intent(in)  :: latice
         336             :   real(r8), intent(in)  :: dcs
         337             :
         338             :   character(128), intent(out) :: errstring
         339             :
         340             :   ! Name this array to workaround an XLF bug (otherwise could just use the
         341             :   ! expression that sets it).
         342             :   real(r8) :: ice_lambda_bounds(2)
         343             :
         344             :   !-----------------------------------------------------------------------
         345             :
         346           0 :   errstring = ' '
         347             :
         348           0 :   if( kind .ne. r8 ) then
         349           0 :      errstring = 'micro_pumas_init: KIND of reals does not match'
         350           0 :      return
         351             :   endif
         352             :
         353             :   ! declarations for MG code (transforms variable names)
         354             :
         355           0 :   rv    = rh2o          ! water vapor gas constant
         356           0 :   cpp   = cpair         ! specific heat of dry air
         357           0 :   tmelt = tmelt_in
         358           0 :   ra    = rair          ! dry air gas constant
         359             :
         360             :   ! latent heats
         361             :
         362           0 :   xxlv = latvap         ! latent heat vaporization
         363           0 :   xlf  = latice         ! latent heat freezing
         364           0 :   xxls = xxlv + xlf     ! latent heat of sublimation
         365             :
         366             :   ! Define constants to help speed up code (this limits calls to gamma function)
         367           0 :   gamma_bs_plus3      = gamma(3._r8+bs)
         368           0 :   gamma_half_br_plus5 = gamma(5._r8/2._r8+br/2._r8)
         369           0 :   gamma_half_bs_plus5 = gamma(5._r8/2._r8+bs/2._r8)
         370           0 :   gamma_2bs_plus2     = gamma(2._r8*bs+2._r8)
         371             :
         372             :   ! Don't specify lambda bounds for cloud liquid, as they are determined by
         373             :   ! pgam dynamically.
         374             :   mg_liq_props = MGHydrometeorProps(rhow, dsph, &
         375           0 :        min_mean_mass=min_mean_mass_liq)
         376             :
         377             :   ! Mean ice diameter can not grow bigger than twice the autoconversion
         378             :   ! threshold for snow.
         379           0 :   ice_lambda_bounds = 1._r8/[2._r8*dcs, 1.e-6_r8]
         380             :
         381             :   mg_ice_props = MGHydrometeorProps(rhoi, dsph, &
         382           0 :        ice_lambda_bounds, min_mean_mass_ice)
         383             :
         384           0 :   mg_rain_props = MGHydrometeorProps(rhow, dsph, lam_bnd_rain)
         385           0 :   mg_snow_props = MGHydrometeorProps(rhosn, dsph, lam_bnd_snow)
         386           0 :   mg_graupel_props = MGHydrometeorProps(rhog, dsph, lam_bnd_snow)
         387           0 :   mg_hail_props = MGHydrometeorProps(rhoh, dsph, lam_bnd_snow)
         388             :
         389             :   !$acc update device (rv,cpp,tmelt,xxlv,xxls,gamma_bs_plus3,   &
         390             :   !$acc                gamma_half_br_plus5,gamma_half_bs_plus5, &
         391             :   !$acc                gamma_2bs_plus2)
         392             :
         393             : end subroutine micro_pumas_utils_init
         394             :
         395             : ! Constructor for a constituent property object.
         396           0 : function NewMGHydrometeorProps(rho, eff_dim, lambda_bounds, min_mean_mass) &
         397             :      result(res)
         398             :   real(r8), intent(in) :: rho, eff_dim
         399             :   real(r8), intent(in), optional :: lambda_bounds(2), min_mean_mass
         400             :   type(MGHydrometeorProps) :: res
         401             :
         402           0 :   res%rho = rho
         403           0 :   res%eff_dim = eff_dim
         404           0 :   if (present(lambda_bounds)) then
         405           0 :      res%lambda_bounds = lambda_bounds
         406             :   else
         407           0 :      res%lambda_bounds = no_limiter()
         408             :   end if
         409           0 :   if (present(min_mean_mass)) then
         410           0 :      res%min_mean_mass = min_mean_mass
         411             :   else
         412           0 :      res%min_mean_mass = no_limiter()
         413             :   end if
         414             :
         415           0 :   res%shape_coef = rho*pi*gamma(eff_dim+1._r8)/6._r8
         416             :
         417           0 : end function NewMGHydrometeorProps
         418             :
         419             : !========================================================================
         420             : !FORMULAS
         421             : !========================================================================
         422             :
         423             : ! Use gamma function to implement rising factorial extended to the reals.
         424           0 : subroutine rising_factorial_r8(x, n, res)
         425             :   !$acc routine seq
         426             :   real(r8), intent(in) :: x, n
         427             :   real(r8), intent(out) :: res
         428             :
         429             :   !$acc data present (x,res)
         430             :
         431           0 :   res = gamma(x+n)/gamma(x)
         432             :
         433             :   !$acc end data
         434           0 : end subroutine rising_factorial_r8
         435             :
         436           0 : subroutine rising_factorial_r8_vec(x, n, res,vlen)
         437             :   integer, intent(in)   :: vlen
         438             :   real(r8), intent(in)  :: x(vlen), n
         439             :   real(r8), intent(out) :: res(vlen)
         440             :   integer :: i
         441             :   real(r8) :: tmp
         442             :
         443             :   !$acc data present (x,res)
         444             :
         445             :   !$acc parallel vector_length(VLENS) default(present)
         446             :   !$acc loop gang vector private(tmp)
         447           0 :   do i=1,vlen
         448           0 :      tmp = x(i)+n
         449           0 :      res(i) = gamma(tmp)
         450           0 :      tmp = gamma(x(i))
         451           0 :      res(i) = res(i)/tmp
         452             :   end do
         453             :   !$acc end parallel
         454             :
         455             :   !$acc end data
         456             :
         457           0 : end subroutine rising_factorial_r8_vec
         458             :
         459             : ! Rising factorial can be performed much cheaper if n is a small integer.
         460           0 : subroutine rising_factorial_integer(x, n, res)
         461             :   !$acc routine seq
         462             :   real(r8), intent(in)  :: x
         463             :   integer,  intent(in)  :: n
         464             :   real(r8), intent(out) :: res
         465             :
         466             :   integer :: i
         467             :   real(r8) :: factor
         468             :
         469             :   !$acc data present (x,res)
         470             :
         471           0 :   res = 1._r8
         472           0 :   factor = x
         473             :
         474             :   !$acc loop seq
         475           0 :   do i = 1, n
         476           0 :      res = res * factor
         477           0 :      factor = factor + 1._r8
         478             :   end do
         479             :
         480             :   !$acc end data
         481           0 : end subroutine rising_factorial_integer
         482             :
         483           0 : subroutine rising_factorial_integer_vec(x, n, res,vlen)
         484             :   integer,  intent(in)  :: vlen
         485             :   real(r8), intent(in)  :: x(vlen)
         486             :   integer,  intent(in)  :: n
         487             :   real(r8), intent(out) :: res(vlen)
         488             :
         489             :   integer  :: i,j
         490           0 :   real(r8) :: factor(vlen)
         491             :
         492             :   !$acc data present (x,res) &
         493             :   !$acc      create  (factor)
         494             :
         495             :   !$acc parallel vector_length(VLENS) default(present)
         496             :   !$acc loop gang vector
         497           0 :   do i=1,vlen
         498           0 :      res(i)    = 1._r8
         499           0 :      factor(i) = x(i)
         500             :   end do
         501             :
         502           0 :   if (n == 3) then
         503             :     !$acc loop gang vector
         504           0 :     do i=1,vlen
         505           0 :        res(i)    = res(i) * factor(i)
         506           0 :        factor(i) = factor(i) + 1._r8
         507           0 :        res(i)    = res(i) * factor(i)
         508           0 :        factor(i) = factor(i) + 1._r8
         509           0 :        res(i)    = res(i) * factor(i)
         510             :     end do
         511           0 :   elseif (n == 2) then
         512             :     !$acc loop gang vector
         513           0 :     do i=1,vlen
         514           0 :        res(i)    = res(i) * factor(i)
         515           0 :        factor(i) = factor(i) + 1._r8
         516           0 :        res(i)    = res(i) * factor(i)
         517             :     end do
         518             :   else
         519             :     !$acc loop seq
         520           0 :     do j = 1, n
         521             :        !$acc loop gang vector
         522           0 :        do i = 1, vlen
         523           0 :           res(i)    = res(i) * factor(i)
         524           0 :           factor(i) = factor(i) + 1._r8
         525             :        end do
         526             :     end do
         527             :   end if
         528             :   !$acc end parallel
         529             :
         530             :   !$acc end data
         531           0 : end subroutine rising_factorial_integer_vec
         532             :
         533             : ! Calculate correction due to latent heat for evaporation/sublimation
         534           0 : subroutine calc_ab_line(t, qv, xxl, ab)
         535             :   !$acc routine seq
         536             :   real(r8), intent(in)  :: t     ! Temperature
         537             :   real(r8), intent(in)  :: qv    ! Saturation vapor pressure
         538             :   real(r8), intent(in)  :: xxl   ! Latent heat
         539             :   real(r8), intent(out) :: ab
         540             :
         541             :   real(r8) :: dqsdt
         542             :
         543             :   !$acc data present (t,qv,xxl,ab)
         544             :
         545           0 :   dqsdt = xxl*qv / (rv * t**2)
         546           0 :   ab = 1._r8 + dqsdt*xxl/cpp
         547             :
         548             :   !$acc end data
         549           0 : end subroutine calc_ab_line
         550             :
         551             : ! Calculate correction due to latent heat for evaporation/sublimation
         552           0 : subroutine calc_ab_vect(t, qv, xxl, ab, vlen)
         553             :   integer,  intent(in) :: vlen
         554             :   real(r8), intent(in) :: t(vlen)     ! Temperature
         555             :   real(r8), intent(in) :: qv(vlen)    ! Saturation vapor pressure
         556             :   real(r8), intent(in) :: xxl         ! Latent heat
         557             :
         558             :   real(r8), intent(out) :: ab(vlen)
         559             :   real(r8) :: dqsdt
         560             :   integer :: i
         561             :
         562             :   !$acc data present (t,qv,xxl,ab)
         563             :
         564             :   !$acc parallel vector_length(VLENS) default(present)
         565             :   !$acc loop gang vector private(dqsdt)
         566           0 :   do i=1,vlen
         567           0 :      dqsdt = xxl*qv(i) / (rv * t(i)**2)
         568           0 :      ab(i) = 1._r8 + dqsdt*xxl/cpp
         569             :   end do
         570             :   !$acc end parallel
         571             :
         572             :   !$acc end data
         573           0 : end subroutine calc_ab_vect
         574             :
         575             : ! get cloud droplet size distribution parameters
         576           0 : subroutine size_dist_param_liq_line(props, qcic, ncic, rho, pgam, lamc)
         577             :   !$acc routine seq
         578             :   type(MGHydrometeorProps), intent(in)    :: props
         579             :   real(r8),                 intent(in)    :: qcic
         580             :   real(r8),                 intent(inout) :: ncic
         581             :   real(r8),                 intent(in)    :: rho
         582             :   real(r8),                 intent(out)   :: pgam
         583             :   real(r8),                 intent(out)   :: lamc
         584             :
         585             :   ! local variables
         586             :   type(MGHydrometeorProps)                :: props_loc
         587             :   real(r8)                                :: tmp
         588             :
         589             :   !$acc data present (props,qcic,ncic,rho,pgam,lamc)
         590             :
         591           0 :   if (qcic > qsmall) then
         592             :      ! Local copy of properties that can be modified.
         593             :      ! (Elemental routines that operate on arrays can't modify scalar
         594             :      ! arguments.)
         595           0 :      props_loc = props
         596             :
         597             :      ! Get pgam from fit Rotstayn and Liu 2003 (changed from Martin 1994 for CAM6)
         598           0 :      pgam = 1.0_r8 - 0.7_r8 * exp(-0.008_r8*1.e-6_r8*ncic*rho)
         599           0 :      pgam = 1._r8/(pgam**2) - 1._r8
         600           0 :      pgam = max(pgam, 2._r8)
         601             :
         602             :      ! Set coefficient for use in size_dist_param_basic.
         603             :      ! The 3D case is so common and optimizable that we specialize it:
         604           0 :      if (props_loc%eff_dim == 3._r8) then
         605           0 :         call rising_factorial(pgam+1._r8, 3, tmp)
         606           0 :         props_loc%shape_coef = pi / 6._r8 * props_loc%rho * tmp
         607             :      else
         608           0 :         call rising_factorial(pgam+1._r8, props_loc%eff_dim, tmp)
         609           0 :         props_loc%shape_coef = pi / 6._r8 * props_loc%rho * tmp
         610             :      end if
         611             :
         612             :      ! Limit to between 2 and 50 microns mean size.
         613           0 :      props_loc%lambda_bounds = (pgam+1._r8)*1._r8/[50.e-6_r8, 2.e-6_r8]
         614             :
         615           0 :      call size_dist_param_basic(props_loc, qcic, ncic, lamc)
         616             :   else
         617             :      ! pgam not calculated in this case, so set it to a value likely to
         618             :      ! cause an error if it is accidentally used
         619             :      ! (gamma function undefined for negative integers)
         620           0 :      pgam = -100._r8
         621           0 :      lamc = 0._r8
         622             :   end if
         623             :
         624             :   !$acc end data
         625           0 : end subroutine size_dist_param_liq_line
         626             :
         627             : ! get cloud droplet size distribution parameters
         628             :
         629           0 : subroutine size_dist_param_liq_2D(props, qcic, ncic, rho, pgam, lamc, dim1, dim2)
         630             :
         631             :   type(mghydrometeorprops),       intent(in)    :: props
         632             :   integer,                        intent(in)    :: dim1, dim2
         633             :   real(r8), dimension(dim1,dim2), intent(in)    :: qcic
         634             :   real(r8), dimension(dim1,dim2), intent(inout) :: ncic
         635             :   real(r8), dimension(dim1,dim2), intent(in)    :: rho
         636             :   real(r8), dimension(dim1,dim2), intent(out)   :: pgam
         637             :   real(r8), dimension(dim1,dim2), intent(out)   :: lamc
         638             :
         639             :   ! local variables
         640             :   integer  :: i, k
         641           0 :   real(r8) :: tmp(dim1,dim2),pgamp1(dim1,dim2)
         642           0 :   real(r8) :: shapeC(dim1,dim2),lbnd(dim1,dim2),ubnd(dim1,dim2)
         643             :
         644             :   !$acc data present (props,qcic,ncic,rho,pgam,lamc) &
         645             :   !$acc      create  (tmp,pgamp1,shapeC,lbnd,ubnd)
         646             :
         647             :   !$acc parallel vector_length(VLENS) default(present)
         648             :   !$acc loop gang vector collapse(2)
         649           0 :   do k = 1, dim2
         650           0 :      do i = 1, dim1
         651           0 :         if (qcic(i,k) > qsmall) then
         652             :            ! Get pgam from fit Rotstayn and Liu 2003 (changed from Martin 1994 for CAM6)
         653           0 :            pgam(i,k) = 1.0_r8 - 0.7_r8 * exp(-0.008_r8*1.e-6_r8*ncic(i,k)*rho(i,k))
         654           0 :            pgam(i,k) = 1._r8/(pgam(i,k)**2) - 1._r8
         655           0 :            pgam(i,k) = max(pgam(i,k), 2._r8)
         656           0 :            pgamp1(i,k) = pgam(i,k)+1._r8
         657             :         else
         658           0 :            pgamp1(i,k) = 0._r8
         659             :         end if
         660             :      end do
         661             :   end do
         662             :   !$acc end parallel
         663             :
         664             :   ! Set coefficient for use in size_dist_param_basic.
         665             :   ! The 3D case is so common and optimizable that we specialize it:
         666           0 :   if (props%eff_dim == 3._r8) then
         667           0 :      call rising_factorial_integer_vec(pgamp1,3,tmp,dim1*dim2)
         668             :   else
         669           0 :      call rising_factorial_r8_vec(pgamp1, props%eff_dim,tmp,dim1*dim2)
         670             :   end if
         671             :
         672             :   !$acc parallel vector_length(VLENS) default(present)
         673             :   !$acc loop gang vector collapse(2)
         674           0 :   do k = 1, dim2
         675           0 :      do i = 1, dim1
         676           0 :         if (qcic(i,k) > qsmall) then
         677           0 :            shapeC(i,k) = pi / 6._r8 * props%rho * tmp(i,k)
         678             :            ! Limit to between 2 and 50 microns mean size.
         679           0 :            lbnd(i,k)   = pgamp1(i,k)*1._r8/50.e-6_r8
         680           0 :            ubnd(i,k)   = pgamp1(i,k)*1._r8/2.e-6_r8
         681             :         else
         682           0 :            shapeC(i,k) = 0._r8
         683           0 :            lbnd(i,k) = 0._r8
         684           0 :            ubnd(i,k) = 0._r8
         685             :         end if
         686             :      end do
         687             :   end do
         688             :   !$acc end parallel
         689             :
         690           0 :   call size_dist_param_basic_vect2(props, qcic, ncic, shapeC, lbnd, ubnd, lamc, dim1*dim2)
         691             :
         692             :   !$acc parallel vector_length(VLENS) default(present)
         693             :   !$acc loop gang vector collapse(2)
         694           0 :   do k = 1, dim2
         695           0 :      do i = 1, dim1
         696           0 :         if (qcic(i,k) <= qsmall) then
         697             :            ! pgam not calculated in this case, so set it to a value likely to
         698             :            ! cause an error if it is accidentally used
         699             :            ! (gamma function undefined for negative integers)
         700           0 :            pgam(i,k) = -100._r8
         701           0 :            lamc(i,k) = 0._r8
         702             :         end if
         703             :      end do
         704             :   end do
         705             :   !$acc end parallel
         706             :
         707             :   !$acc end data
         708           0 : end subroutine size_dist_param_liq_2D
         709             :
         710             : ! get cloud droplet size distribution parameters
         711             :
         712           0 : subroutine size_dist_param_liq_vect(props, qcic, ncic, rho, pgam, lamc, vlen)
         713             :
         714             :   type(mghydrometeorprops),  intent(in)    :: props
         715             :   integer,                   intent(in)    :: vlen
         716             :   real(r8), dimension(vlen), intent(in)    :: qcic
         717             :   real(r8), dimension(vlen), intent(inout) :: ncic
         718             :   real(r8), dimension(vlen), intent(in)    :: rho
         719             :   real(r8), dimension(vlen), intent(out)   :: pgam
         720             :   real(r8), dimension(vlen), intent(out)   :: lamc
         721             :
         722             :   ! local variables
         723             :   integer  :: i
         724           0 :   real(r8) :: tmp(vlen),pgamp1(vlen)
         725           0 :   real(r8) :: shapeC(vlen),lbnd(vlen),ubnd(vlen)
         726             :
         727             :   !$acc data present (props,qcic,ncic,rho,pgam,lamc) &
         728             :   !$acc      create  (tmp,pgamp1,shapeC,lbnd,ubnd)
         729             :
         730             :   !$acc parallel vector_length(VLENS) default(present)
         731             :   !$acc loop gang vector
         732           0 :   do i = 1, vlen
         733           0 :      if (qcic(i) > qsmall) then
         734             :         ! Get pgam from fit Rotstayn and Liu 2003 (changed from Martin 1994 for CAM6)
         735           0 :         pgam(i) = 1.0_r8 - 0.7_r8 * exp(-0.008_r8*1.e-6_r8*ncic(i)*rho(i))
         736           0 :         pgam(i) = 1._r8/(pgam(i)**2) - 1._r8
         737           0 :         pgam(i) = max(pgam(i), 2._r8)
         738           0 :         pgamp1(i) = pgam(i)+1._r8
         739             :      else
         740           0 :         pgamp1(i) = 0._r8
         741             :      end if
         742             :   end do
         743             :   !$acc end parallel
         744             :
         745             :   ! Set coefficient for use in size_dist_param_basic.
         746             :   ! The 3D case is so common and optimizable that we specialize it:
         747           0 :   if (props%eff_dim == 3._r8) then
         748           0 :      call rising_factorial_integer_vec(pgamp1,3,tmp,vlen)
         749             :   else
         750           0 :      call rising_factorial_r8_vec(pgamp1, props%eff_dim,tmp,vlen)
         751             :   end if
         752             :
         753             :   !$acc parallel vector_length(VLENS) default(present)
         754             :   !$acc loop gang vector
         755           0 :   do i = 1, vlen
         756           0 :      if (qcic(i) > qsmall) then
         757           0 :         shapeC(i) = pi / 6._r8 * props%rho * tmp(i)
         758             :         ! Limit to between 2 and 50 microns mean size.
         759           0 :         lbnd(i)   = pgamp1(i)*1._r8/50.e-6_r8
         760           0 :         ubnd(i)   = pgamp1(i)*1._r8/2.e-6_r8
         761             :      else
         762           0 :         shapeC(i) = 0._r8
         763           0 :         lbnd(i) = 0._r8
         764           0 :         ubnd(i) = 0._r8
         765             :      end if
         766             :   end do
         767             :   !$acc end parallel
         768             :
         769           0 :   call size_dist_param_basic_vect2(props, qcic, ncic, shapeC, lbnd, ubnd, lamc, vlen)
         770             :
         771             :   !$acc parallel vector_length(VLENS) default(present)
         772             :   !$acc loop gang vector
         773           0 :   do i = 1, vlen
         774           0 :      if (qcic(i) <= qsmall) then
         775             :         ! pgam not calculated in this case, so set it to a value likely to
         776             :         ! cause an error if it is accidentally used
         777             :         ! (gamma function undefined for negative integers)
         778           0 :         pgam(i) = -100._r8
         779           0 :         lamc(i) = 0._r8
         780             :      end if
         781             :   end do
         782             :   !$acc end parallel
         783             :
         784             :   !$acc end data
         785           0 : end subroutine size_dist_param_liq_vect
         786             :
         787             : ! Basic routine for getting size distribution parameters.
         788           0 : subroutine size_dist_param_basic_line(props, qic, nic, lam, n0)
         789             :   !$acc routine seq
         790             :   type(MGHydrometeorProps), intent(in)    :: props
         791             :   real(r8),                 intent(in)    :: qic
         792             :   real(r8),                 intent(inout) :: nic
         793             :   real(r8),                 intent(out)           :: lam
         794             :   real(r8),                 intent(out), optional :: n0
         795             :
         796             :   logical :: present_n0
         797           0 :   present_n0 = present(n0)
         798             :
         799             :   !$acc data present (props,qic,nic,lam,n0)
         800             :
         801           0 :   if (qic > qsmall) then
         802             :      ! add upper limit to in-cloud number concentration to prevent
         803             :      ! numerical error
         804           0 :      if (limiter_is_on(props%min_mean_mass)) then
         805           0 :         nic = min(nic, qic / props%min_mean_mass)
         806             :      end if
         807             :
         808             :      ! lambda = (c n/q)^(1/d)
         809           0 :      lam = (props%shape_coef * nic/qic)**(1._r8/props%eff_dim)
         810             :
         811             :      ! check for slope
         812             :      ! adjust vars
         813           0 :      if (lam < props%lambda_bounds(1)) then
         814           0 :         lam = props%lambda_bounds(1)
         815           0 :         nic = lam**(props%eff_dim) * qic/props%shape_coef
         816           0 :      else if (lam > props%lambda_bounds(2)) then
         817           0 :         lam = props%lambda_bounds(2)
         818           0 :         nic = lam**(props%eff_dim) * qic/props%shape_coef
         819             :      end if
         820             :   else
         821           0 :      lam = 0._r8
         822             :   end if
         823             :
         824           0 :   if (present_n0) n0 = nic * lam
         825             :
         826             :   !$acc end data
         827           0 : end subroutine size_dist_param_basic_line
         828             :
         829           0 : subroutine size_dist_param_basic_vect(props, qic, nic, lam, vlen, n0)
         830             :
         831             :   type (mghydrometeorprops), intent(in)    :: props
         832             :   integer,                   intent(in)    :: vlen
         833             :   real(r8), dimension(vlen), intent(in)    :: qic
         834             :   real(r8), dimension(vlen), intent(inout) :: nic
         835             :   real(r8), dimension(vlen), intent(out)   :: lam
         836             :   real(r8), dimension(vlen), intent(out), optional :: n0
         837             :   integer  :: i
         838             :   logical  :: limiterActive, present_n0
         839             :   real(r8) :: effDim,shapeCoef,ubnd,lbnd, minMass
         840             :
         841             :   !$acc data present (props,qic,nic,lam,n0)
         842             :
         843           0 :   limiterActive = limiter_is_on(props%min_mean_mass)
         844           0 :   effDim    = props%eff_dim
         845           0 :   shapeCoef = props%shape_coef
         846           0 :   lbnd      = props%lambda_bounds(1)
         847           0 :   ubnd      = props%lambda_bounds(2)
         848           0 :   minMass   = props%min_mean_mass
         849           0 :   present_n0 = present(n0)
         850             :
         851             :   !$acc parallel vector_length(VLENS) default(present)
         852             :   !$acc loop gang vector
         853           0 :   do i = 1, vlen
         854           0 :      if (qic(i) > qsmall) then
         855             :         ! add upper limit to in-cloud number concentration to prevent
         856             :         ! numerical error
         857           0 :         if (limiterActive) then
         858           0 :            nic(i) = min(nic(i), qic(i) / minMass)
         859             :         end if
         860             :
         861             :         ! lambda = (c n/q)^(1/d)
         862           0 :         lam(i) = (shapeCoef * nic(i)/qic(i))**(1._r8/effDim)
         863             :
         864             :         ! check for slope
         865             :         ! adjust vars
         866           0 :         if (lam(i) < lbnd) then
         867           0 :            lam(i) = lbnd
         868           0 :            nic(i) = lam(i)**(effDim) * qic(i)/shapeCoef
         869           0 :         else if (lam(i) > ubnd) then
         870           0 :            lam(i) = ubnd
         871           0 :            nic(i) = lam(i)**(effDim) * qic(i)/shapeCoef
         872             :         end if
         873             :
         874             :      else
         875           0 :         lam(i) = 0._r8
         876             :      end if
         877             :   end do
         878             :
         879           0 :   if (present_n0) then
         880             :      !$acc loop gang vector
         881           0 :      do i = 1, vlen
         882           0 :         n0(i) = nic(i) * lam(i)
         883             :      end do
         884             :   end if
         885             :   !$acc end parallel
         886             :
         887             :   !$acc end data
         888           0 : end subroutine size_dist_param_basic_vect
         889             :
         890           0 : subroutine size_dist_param_basic_2D(props, qic, nic, lam, dim1, dim2, n0)
         891             :
         892             :   type (mghydrometeorprops),      intent(in)    :: props
         893             :   integer,                        intent(in)    :: dim1, dim2
         894             :   real(r8), dimension(dim1,dim2), intent(in)    :: qic
         895             :   real(r8), dimension(dim1,dim2), intent(inout) :: nic
         896             :   real(r8), dimension(dim1,dim2), intent(out)   :: lam
         897             :   real(r8), dimension(dim1,dim2), intent(out), optional :: n0
         898             :   integer  :: i, k
         899             :   logical  :: limiterActive, present_n0
         900             :   real(r8) :: effDim,shapeCoef,ubnd,lbnd, minMass
         901             :
         902           0 :   limiterActive = limiter_is_on(props%min_mean_mass)
         903           0 :   effDim    = props%eff_dim
         904           0 :   shapeCoef = props%shape_coef
         905           0 :   lbnd      = props%lambda_bounds(1)
         906           0 :   ubnd      = props%lambda_bounds(2)
         907           0 :   minMass   = props%min_mean_mass
         908           0 :   present_n0 = present(n0)
         909             :
         910             :   !$acc data present (props,qic,nic,lam,n0)
         911             :
         912             :   !$acc parallel vector_length(VLENS) default(present)
         913             :   !$acc loop gang vector collapse(2)
         914           0 :   do k = 1, dim2
         915           0 :      do i = 1, dim1
         916           0 :         if (qic(i,k) > qsmall) then
         917             :            ! add upper limit to in-cloud number concentration to prevent
         918             :            ! numerical error
         919           0 :            if (limiterActive) then
         920           0 :               nic(i,k) = min(nic(i,k), qic(i,k) / minMass)
         921             :            end if
         922             :
         923             :            ! lambda = (c n/q)^(1/d)
         924           0 :            lam(i,k) = (shapeCoef * nic(i,k)/qic(i,k))**(1._r8/effDim)
         925             :
         926             :            ! check for slope
         927             :            ! adjust vars
         928           0 :            if (lam(i,k) < lbnd) then
         929           0 :               lam(i,k) = lbnd
         930           0 :               nic(i,k) = lam(i,k)**(effDim) * qic(i,k)/shapeCoef
         931           0 :            else if (lam(i,k) > ubnd) then
         932           0 :               lam(i,k) = ubnd
         933           0 :               nic(i,k) = lam(i,k)**(effDim) * qic(i,k)/shapeCoef
         934             :            end if
         935             :
         936             :         else
         937           0 :            lam(i,k) = 0._r8
         938             :         end if
         939             :      end do
         940             :   end do
         941             :
         942           0 :   if (present_n0) then
         943             :      !$acc loop gang vector collapse(2)
         944           0 :      do k = 1, dim2
         945           0 :         do i = 1, dim1
         946           0 :            n0(i,k) = nic(i,k) * lam(i,k)
         947             :         end do
         948             :      end do
         949             :   end if
         950             :   !$acc end parallel
         951             :
         952             :   !$acc end data
         953           0 : end subroutine size_dist_param_basic_2D
         954             :
         955           0 : subroutine size_dist_param_basic_vect2(props, qic, nic, shapeC, lbnd, ubnd, lam, vlen, n0)
         956             :
         957             :   type (mghydrometeorprops), intent(in)    :: props
         958             :   integer,                   intent(in)    :: vlen
         959             :   real(r8), dimension(vlen), intent(in)    :: qic
         960             :   real(r8), dimension(vlen), intent(inout) :: nic
         961             :   real(r8), dimension(vlen), intent(in)    :: shapeC,lbnd,ubnd
         962             :   real(r8), dimension(vlen), intent(out)   :: lam
         963             :   real(r8), dimension(vlen), intent(out), optional :: n0
         964             :   integer  :: i
         965             :   integer  :: cnt
         966             :   logical  :: limiterActive, present_n0
         967             :   real(r8) :: effDim,shapeCoef, minMass
         968             :
         969           0 :   limiterActive = limiter_is_on(props%min_mean_mass)
         970           0 :   effDim        = props%eff_dim
         971           0 :   minMass       = props%min_mean_mass
         972           0 :   present_n0    = present(n0)
         973             :
         974             :   !$acc data present (props,qic,nic,shapeC,lbnd,ubnd,lam,n0)
         975             :
         976             :   !$acc parallel vector_length(VLENS) default(present)
         977             :   !$acc loop gang vector
         978           0 :   do i=1,vlen
         979             :
         980           0 :      if (qic(i) > qsmall) then
         981             :         ! add upper limit to in-cloud number concentration to prevent
         982             :         ! numerical error
         983             :
         984           0 :         if (limiterActive) then
         985           0 :            nic(i) = min(nic(i), qic(i) / minMass)
         986             :         end if
         987             :         ! lambda = (c n/q)^(1/d)
         988             :
         989           0 :         lam(i) = (shapeC(i) * nic(i)/qic(i))**(1._r8/effDim)
         990             :         ! check for slope
         991             :         ! adjust vars
         992             :
         993           0 :         if (lam(i) < lbnd(i)) then
         994           0 :            lam(i) = lbnd(i)
         995           0 :            nic(i) = lam(i)**(effDim) * qic(i)/shapeC(i)
         996           0 :         else if (lam(i) > ubnd(i)) then
         997           0 :            lam(i) = ubnd(i)
         998           0 :            nic(i) = lam(i)**(effDim) * qic(i)/shapeC(i)
         999             :         end if
        1000             :      else
        1001           0 :         lam(i) = 0._r8
        1002             :      end if
        1003             :
        1004             :   end do
        1005             :
        1006           0 :   if (present_n0) then
        1007             :      !$acc loop gang vector
        1008           0 :      do i = 1,vlen
        1009           0 :         n0(i) = nic(i) * lam(i)
        1010             :      end do
        1011             :   end if
        1012             :   !$acc end parallel
        1013             :
        1014             :   !$acc end data
        1015           0 : end subroutine size_dist_param_basic_vect2
        1016             :
        1017           0 : real(r8) elemental function avg_diameter(q, n, rho_air, rho_sub)
        1018             :   ! Finds the average diameter of particles given their density, and
        1019             :   ! mass/number concentrations in the air.
        1020             :   ! Assumes that diameter follows an exponential distribution.
        1021             :   real(r8), intent(in) :: q         ! mass mixing ratio
        1022             :   real(r8), intent(in) :: n         ! number concentration (per volume)
        1023             :   real(r8), intent(in) :: rho_air   ! local density of the air
        1024             :   real(r8), intent(in) :: rho_sub   ! density of the particle substance
        1025             :
        1026           0 :   avg_diameter = (q*rho_air/(pi * rho_sub * n))**(1._r8/3._r8)
        1027             :
        1028           0 : end function avg_diameter
        1029             :
        1030           0 : subroutine avg_diameter_vec (q, n, rho_air, rho_sub, avg_diameter, vlen)
        1031             :    ! Finds the average diameter of particles given their density, and
        1032             :    ! mass/number concentrations in the air.
        1033             :    ! Assumes that diameter follows an exponential distribution.
        1034             :    integer,  intent(in)  :: vlen
        1035             :    real(r8), intent(in)  :: q(vlen)         ! mass mixing ratio
        1036             :    real(r8), intent(in)  :: n(vlen)         ! number concentration (per volume)
        1037             :    real(r8), intent(in)  :: rho_air(vlen)   ! local density of the air
        1038             :    real(r8), intent(in)  :: rho_sub   ! density of the particle substance
        1039             :    real(r8), intent(out) :: avg_diameter(vlen)
        1040             :    integer :: i
        1041             :
        1042             :    !$acc data present (q,n,rho_air,avg_diameter)
        1043             :
        1044             :    !$acc parallel vector_length(VLENS) default(present)
        1045             :    !$acc loop gang vector
        1046           0 :    do i=1,vlen
        1047           0 :       avg_diameter(i) = (q(i)*rho_air(i)/(pi * rho_sub * n(i)))**(1._r8/3._r8)
        1048             :    end do
        1049             :    !$acc end parallel
        1050             :
        1051             :    !$acc end data
        1052           0 : end subroutine avg_diameter_vec
        1053             :
        1054           0 : subroutine var_coef_r8(relvar, a, res)
        1055             :   !$acc routine seq
        1056             :
        1057             :   ! Finds a coefficient for process rates based on the relative variance
        1058             :   ! of cloud water.
        1059             :   real(r8), intent(in)  :: relvar
        1060             :   real(r8), intent(in)  :: a
        1061             :   real(r8), intent(out) :: res
        1062             :
        1063             :   !$acc data present (relvar,res)
        1064             :
        1065           0 :   call rising_factorial(relvar, a, res)
        1066           0 :   res = res / relvar**a
        1067             :
        1068             :   !$acc end data
        1069           0 : end subroutine var_coef_r8
        1070             :
        1071           0 : subroutine var_coef_r8_vect(relvar, a, res, vlen)
        1072             :   ! Finds a coefficient for process rates based on the relative variance
        1073             :   ! of cloud water.
        1074             :   integer,  intent(in)  :: vlen
        1075             :   real(r8), intent(in)  :: relvar(vlen)
        1076             :   real(r8), intent(in)  :: a
        1077             :   real(r8), intent(out) :: res(vlen)
        1078             :   integer  :: i
        1079           0 :   real(r8) :: tmpA(vlen)
        1080             :
        1081             :   !$acc data present (relvar,res) &
        1082             :   !$acc      create  (tmpA)
        1083             :
        1084           0 :   call rising_factorial(relvar,a,tmpA,vlen)
        1085             :   !$acc parallel vector_length(VLENS) default(present)
        1086             :   !$acc loop gang vector
        1087           0 :   do i=1,vlen
        1088           0 :      res(i) = tmpA(i)/relvar(i)**a
        1089             :   end do
        1090             :   !$acc end parallel
        1091             :
        1092             :   !$acc end data
        1093           0 : end subroutine var_coef_r8_vect
        1094             :
        1095             : subroutine var_coef_integer(relvar, a, res)
        1096             :   !$acc routine seq
        1097             :
        1098             :   ! Finds a coefficient for process rates based on the relative variance
        1099             :   ! of cloud water.
        1100             :   real(r8), intent(in) :: relvar
        1101             :   integer, intent(in) :: a
        1102             :   real(r8), intent(out) :: res
        1103             :
        1104             :   !$acc data present (relvar,res)
        1105             :
        1106             :   call rising_factorial(relvar, a, res)
        1107             :   res = res / relvar**a
        1108             :
        1109             :   !$acc end data
        1110             : end subroutine var_coef_integer
        1111             :
        1112           0 : subroutine var_coef_integer_vect(relvar, a, res, vlen)
        1113             :   ! Finds a coefficient for process rates based on the relative variance
        1114             :   ! of cloud water.
        1115             :   integer, intent(in)   :: vlen
        1116             :   real(r8), intent(in)  :: relvar(vlen)
        1117             :   integer, intent(in)   :: a
        1118             :   real(r8), intent(out) :: res(vlen)
        1119             :   integer  :: i
        1120           0 :   real(r8) :: tmp(vlen)
        1121             :
        1122             :   !$acc data present (relvar,res) &
        1123             :   !$acc      create  (tmp)
        1124             :
        1125           0 :   call rising_factorial(relvar, a,tmp,vlen)
        1126             :   !$acc parallel vector_length(VLENS) default(present)
        1127             :   !$acc loop gang vector
        1128           0 :   do i=1,vlen
        1129           0 :      res(i) = tmp(i) / relvar(i)**a
        1130             :   end do
        1131             :   !$acc end parallel
        1132             :
        1133             :   !$acc end data
        1134           0 : end subroutine var_coef_integer_vect
        1135             :
        1136             : !========================================================================
        1137             : !MICROPHYSICAL PROCESS CALCULATIONS
        1138             : !========================================================================
        1139             : !========================================================================
        1140             : ! Initial ice deposition and sublimation loop.
        1141             : ! Run before the main loop
        1142             : ! This subroutine written by Peter Caldwell
        1143             :
        1144           0 : subroutine ice_deposition_sublimation(t, qv, qi, ni, &
        1145           0 :                                       icldm, rho, dv,qvl, qvi, &
        1146           0 :                                       berg, vap_dep, ice_sublim, vlen)
        1147             :
        1148             :   !INPUT VARS:
        1149             :   !===============================================
        1150             :   integer,  intent(in)                  :: vlen
        1151             :   real(r8), dimension(vlen), intent(in) :: t
        1152             :   real(r8), dimension(vlen), intent(in) :: qv
        1153             :   real(r8), dimension(vlen), intent(in) :: qi
        1154             :   real(r8), dimension(vlen), intent(in) :: ni
        1155             :   real(r8), dimension(vlen), intent(in) :: icldm
        1156             :   real(r8), dimension(vlen), intent(in) :: rho
        1157             :   real(r8), dimension(vlen), intent(in) :: dv
        1158             :   real(r8), dimension(vlen), intent(in) :: qvl
        1159             :   real(r8), dimension(vlen), intent(in) :: qvi
        1160             :
        1161             :   !OUTPUT VARS:
        1162             :   !===============================================
        1163             :   real(r8), dimension(vlen), intent(out) :: vap_dep !ice deposition (cell-ave value)
        1164             :   real(r8), dimension(vlen), intent(out) :: ice_sublim !ice sublimation (cell-ave value)
        1165             :   real(r8), dimension(vlen), intent(out) :: berg !bergeron enhancement (cell-ave value)
        1166             :
        1167             :   !INTERNAL VARS:
        1168             :   !===============================================
        1169           0 :   real(r8) :: ab(vlen)
        1170             :   real(r8) :: epsi
        1171           0 :   real(r8) :: qiic(vlen)
        1172           0 :   real(r8) :: niic(vlen)
        1173           0 :   real(r8) :: lami(vlen)
        1174           0 :   real(r8) :: n0i(vlen)
        1175             :   integer :: i
        1176             :
        1177             :   !$acc data present (t,qv,qi,ni,icldm,rho,dv,qvl) &
        1178             :   !$acc      present (qvi,vap_dep,ice_sublim,berg) &
        1179             :   !$acc      create  (ab,qiic,niic,lami,n0i)
        1180             :
        1181             :   !GET IN-CLOUD qi, ni
        1182             :   !===============================================
        1183             :
        1184             :   !Compute linearized condensational heating correction
        1185           0 :   call calc_ab(t, qvi, xxls, ab, vlen)
        1186             :
        1187             :   !$acc parallel vector_length(VLENS) default(present)
        1188             :   !$acc loop gang vector
        1189           0 :   do i = 1,vlen
        1190           0 :      qiic(i) = qi(i)/icldm(i)
        1191           0 :      niic(i) = ni(i)/icldm(i)
        1192           0 :      if (qi(i)<qsmall) then
        1193           0 :         qiic(i) = 0._r8
        1194           0 :         niic(i) = 0._r8
        1195           0 :         ab(i)   = 0._r8
        1196             :      end if
        1197             :   end do
        1198             :   !$acc end parallel
        1199             :
        1200             :   !Get slope and intercept of gamma distn for ice.
        1201           0 :   call size_dist_param_basic_vect(mg_ice_props, qiic, niic, lami, vlen, n0i)
        1202             :
        1203             :   !$acc parallel vector_length(VLENS) default(present)
        1204             :   !$acc loop gang vector private(epsi)
        1205           0 :   do i=1,vlen
        1206           0 :      if (qi(i)>=qsmall) then
        1207             :         !Get depletion timescale=1/eps
        1208           0 :         epsi = 2._r8*pi*n0i(i)*rho(i)*Dv(i)/(lami(i)*lami(i))
        1209             :
        1210             :         !Compute deposition/sublimation
        1211           0 :         vap_dep(i) = epsi/ab(i)*(qv(i) - qvi(i))
        1212             :
        1213             :         !Make this a grid-averaged quantity
        1214           0 :         vap_dep(i) = vap_dep(i)*icldm(i)
        1215             :
        1216             :         !Split into deposition or sublimation.
        1217           0 :         if (t(i) < tmelt .and. vap_dep(i) > 0._r8) then
        1218           0 :            ice_sublim(i) = 0._r8
        1219             :         else
        1220             :         ! make ice_sublim negative for consistency with other evap/sub processes
        1221           0 :            ice_sublim(i) = min(vap_dep(i),0._r8)
        1222           0 :            vap_dep(i)    = 0._r8
        1223             :         end if
        1224             :
        1225             :         !sublimation occurs @ any T. Not so for berg.
        1226           0 :         if (t(i) < tmelt) then
        1227             :
        1228             :            !Compute bergeron rate assuming cloud for whole step.
        1229           0 :            berg(i) = max(epsi/ab(i)*(qvl(i) - qvi(i)), 0._r8)
        1230             :         else   !T>frz
        1231           0 :            berg(i)=0._r8
        1232             :         end if !T<frz
        1233             :      else      !where qi<qsmall
        1234           0 :         berg(i)       = 0._r8
        1235           0 :         vap_dep(i)    = 0._r8
        1236           0 :         ice_sublim(i) = 0._r8
        1237             :      end if    !qi>qsmall
        1238             :   end do
        1239             :   !$acc end parallel
        1240             :
        1241             :   !$acc end data
        1242           0 : end subroutine ice_deposition_sublimation
        1243             :
        1244           0 : subroutine ice_deposition_sublimation_mg4(t, qv, qi, niic, &
        1245           0 :                                       icldm, rho, dv,qvl, qvi, &
        1246           0 :                                       berg, vap_dep, ice_sublim, &
        1247           0 :                                       af1pr5, af1pr14, rhof, mu, sc, &
        1248             :                                       vlen)
        1249             :
        1250             :   !INPUT VARS:
        1251             :   !===============================================
        1252             :   integer,  intent(in)                  :: vlen
        1253             :   real(r8), dimension(vlen), intent(in) :: t
        1254             :   real(r8), dimension(vlen), intent(in) :: qv
        1255             :   real(r8), dimension(vlen), intent(in) :: qi
        1256             :   real(r8), dimension(vlen), intent(in) :: niic ! trude: input nicc as other routines
        1257             :   real(r8), dimension(vlen), intent(in) :: icldm
        1258             :   real(r8), dimension(vlen), intent(in) :: rho
        1259             :   real(r8), dimension(vlen), intent(in) :: dv
        1260             :   real(r8), dimension(vlen), intent(in) :: qvl
        1261             :   real(r8), dimension(vlen), intent(in) :: qvi
        1262             :   real(r8), dimension(vlen), intent(in) :: af1pr5, af1pr14
        1263             :   real(r8), dimension(vlen), intent(in) :: rhof
        1264             :   real(r8), dimension(vlen), intent(in) :: mu
        1265             :   real(r8), dimension(vlen), intent(in) :: sc
        1266             :
        1267             :   !OUTPUT VARS:
        1268             :   !===============================================
        1269             :   real(r8), dimension(vlen), intent(out) :: vap_dep !ice deposition (cell-ave value)
        1270             :   real(r8), dimension(vlen), intent(out) :: ice_sublim !ice sublimation (cell-ave value)
        1271             :   real(r8), dimension(vlen), intent(out) :: berg !bergeron enhancement (cell-ave value)
        1272             :
        1273             :   !INTERNAL VARS:
        1274             :   !===============================================
        1275           0 :   real(r8) :: ab(vlen)
        1276             :   real(r8) :: epsi
        1277           0 :   real(r8) :: qiic(vlen)
        1278             :   real(r8) :: lami(vlen)
        1279             :   real(r8) :: n0i(vlen)
        1280             :   real(r8), parameter :: thrd = 1._r8/3._r8
        1281             :   integer :: i
        1282             :
        1283             :   !$acc data present (t,qv,qi,icldm,rho,dv,qvl) &
        1284             :   !$acc      present (qvi,vap_dep,ice_sublim,berg) &
        1285             :   !$acc      present (niic,af1pr5,af1pr14,rhof,mu,sc) &
        1286             :   !$acc      create  (ab,qiic,lami,n0i)
        1287             :
        1288             :   !GET IN-CLOUD qi, ni
        1289             :   !===============================================
        1290             :   !$acc parallel vector_length(VLENS) default(present)
        1291             :   !$acc loop gang vector
        1292           0 :   do i = 1,vlen
        1293           0 :      if (qi(i)>=qsmall) then
        1294           0 :         qiic(i) = qi(i)/icldm(i)
        1295             :         !Compute linearized condensational heating correction
        1296             :         call calc_ab(t(i), qvi(i), xxls, ab(i))
        1297             :      end if
        1298             :   end do
        1299             :   !$acc end parallel
        1300             :
        1301             :   !$acc parallel vector_length(VLENS) default(present)
        1302             :   !$acc loop gang vector private(epsi)
        1303           0 :   do i=1,vlen
        1304           0 :      if (qi(i)>=qsmall) then
        1305           0 :         if ( t(i) .lt. tmelt ) then
        1306             :            epsi = (af1pr5(i)+af1pr14(i)*sc(i)**thrd*(rhof(i)*rho(i)/mu(i))**0.5_r8)*2._r8*pi* &
        1307           0 :              rho(i)*dv(i)
        1308             :         else
        1309             :            epsi = 0._r8
        1310             :         end if
        1311             :
        1312             :         !Compute deposition/sublimation
        1313           0 :         vap_dep(i) = epsi/ab(i)*(qv(i) - qvi(i))
        1314             :
        1315             :         !Make this a grid-averaged quantity
        1316           0 :         vap_dep(i) = vap_dep(i)*icldm(i)
        1317             :
        1318             :         !Split into deposition or sublimation.
        1319           0 :         if (t(i) < tmelt .and. vap_dep(i) > 0._r8) then
        1320           0 :            ice_sublim(i)=0._r8
        1321             :         else
        1322             :            ! make ice_sublim negative for consistency with other evap/sub processes
        1323           0 :            ice_sublim(i)=min(vap_dep(i),0._r8)
        1324           0 :            vap_dep(i)=0._r8
        1325             :         end if
        1326             :
        1327             :         !sublimation occurs @ any T. Not so for berg.
        1328           0 :         if (t(i) < tmelt) then
        1329             :            !Compute bergeron rate assuming cloud for whole step.
        1330           0 :            berg(i) = max(epsi/ab(i)*(qvl(i) - qvi(i)), 0._r8)
        1331             :         else   !T>frz
        1332           0 :            berg(i) = 0._r8
        1333             :         end if !T<frz
        1334             :      else      !where qi<qsmall
        1335           0 :         berg(i)       = 0._r8
        1336           0 :         vap_dep(i)    = 0._r8
        1337           0 :         ice_sublim(i) = 0._r8
        1338             :      end if    !qi>qsmall
        1339             :   end do
        1340             :   !$acc end parallel
        1341             :
        1342             :   !$acc end data
        1343           0 : end subroutine ice_deposition_sublimation_mg4
        1344             :
        1345             : !========================================================================
        1346             : ! autoconversion of cloud liquid water to rain
        1347             : ! formula from Khrouditnov and Kogan (2000), modified for sub-grid distribution of qc
        1348             : ! minimum qc of 1 x 10^-8 prevents floating point error
        1349             :
        1350           0 : subroutine kk2000_liq_autoconversion(microp_uniform, qcic, &
        1351           0 :      ncic, rho, relvar, prc, nprc, nprc1, micro_mg_autocon_fact, &
        1352             :      micro_mg_autocon_nd_exp, micro_mg_autocon_lwp_exp, vlen)
        1353             :
        1354             :   integer, intent(in) :: vlen
        1355             :   logical, intent(in) :: microp_uniform
        1356             :
        1357             :   real(r8), dimension(vlen), intent(in) :: qcic
        1358             :   real(r8), dimension(vlen), intent(in) :: ncic
        1359             :   real(r8), dimension(vlen), intent(in) :: rho
        1360             :
        1361             :   real(r8), dimension(vlen), intent(in) :: relvar
        1362             :
        1363             :   real(r8), dimension(vlen), intent(out) :: prc
        1364             :   real(r8), dimension(vlen), intent(out) :: nprc
        1365             :   real(r8), dimension(vlen), intent(out) :: nprc1
        1366             :
        1367           0 :   real(r8), dimension(vlen) :: prc_coef
        1368             :   integer :: i
        1369             :
        1370             :   real(r8), intent(in) :: micro_mg_autocon_fact
        1371             :   real(r8), intent(in) :: micro_mg_autocon_nd_exp
        1372             :   real(r8), intent(in) :: micro_mg_autocon_lwp_exp
        1373             :
        1374             :   !$acc data present (qcic,ncic,rho,relvar,prc,nprc,nprc1) &
        1375             :   !$acc      create  (prc_coef)
        1376             :
        1377             :   ! Take variance into account, or use uniform value.
        1378           0 :   if (.not. microp_uniform) then
        1379           0 :      call var_coef(relvar, micro_mg_autocon_lwp_exp, prc_coef,vlen)
        1380             :   else
        1381             :      !$acc parallel vector_length(VLENS) default(present)
        1382             :      !$acc loop gang vector
        1383           0 :      do i = 1,vlen
        1384           0 :         prc_coef(i) = 1._r8
        1385             :      end do
        1386             :      !$acc end parallel
        1387             :   end if
        1388             :
        1389             :   !$acc parallel vector_length(VLENS) default(present)
        1390             :   !$acc loop gang vector
        1391           0 :   do i=1,vlen
        1392           0 :      if (qcic(i) >= icsmall) then
        1393             :         ! nprc is increase in rain number conc due to autoconversion
        1394             :         ! nprc1 is decrease in cloud droplet conc due to autoconversion
        1395             :
        1396             :         ! assume exponential sub-grid distribution of qc, resulting in additional
        1397             :         ! factor related to qcvar below
        1398             :         ! switch for sub-columns, don't include sub-grid qc
        1399             :         prc(i)   = prc_coef(i) * &
        1400             :              micro_mg_autocon_fact * 1350._r8 * qcic(i)**micro_mg_autocon_lwp_exp * &
        1401           0 :              (ncic(i)*1.e-6_r8*rho(i))**(micro_mg_autocon_nd_exp)
        1402             :
        1403           0 :         nprc(i)  = prc(i) * (1._r8/droplet_mass_25um)
        1404           0 :         nprc1(i) = prc(i)*ncic(i)/qcic(i)
        1405             :
        1406             :      else
        1407           0 :         prc(i)   = 0._r8
        1408           0 :         nprc(i)  = 0._r8
        1409           0 :         nprc1(i) = 0._r8
        1410             :      end if
        1411             :   end do
        1412             :   !$acc end parallel
        1413             :
        1414             :   !$acc end data
        1415           0 : end subroutine kk2000_liq_autoconversion
        1416             :
        1417             : !========================================================================
        1418             :
        1419           0 : subroutine sb2001v2_liq_autoconversion(pgam,qc,nc,qr,rho,relvar,au,nprc,nprc1,vlen)
        1420             :   !
        1421             :   ! ---------------------------------------------------------------------
        1422             :   ! AUTO_SB:  calculates the evolution of mass- and number mxg-ratio for
        1423             :   ! drizzle drops due to autoconversion. The autoconversion rate assumes
        1424             :   ! f(x)=A*x**(nu_c)*exp(-Bx) in drop MASS x.
        1425             :
        1426             :   ! Code from Hugh Morrison, Sept 2014
        1427             :
        1428             :   ! autoconversion
        1429             :   ! use simple lookup table of dnu values to get mass spectral shape parameter
        1430             :   ! equivalent to the size spectral shape parameter pgam
        1431             :
        1432             :   integer,                   intent(in)     :: vlen
        1433             :   real(r8), dimension(vlen), intent (in)    :: pgam
        1434             :   real(r8), dimension(vlen), intent (in)    :: qc  ! = qc (cld water mixing ratio)
        1435             :   real(r8), dimension(vlen), intent (in)    :: nc  ! = nc (cld water number conc /kg)
        1436             :   real(r8), dimension(vlen), intent (in)    :: qr  ! = qr (rain water mixing ratio)
        1437             :   real(r8), dimension(vlen), intent (in)    :: rho ! = rho : density profile
        1438             :   real(r8), dimension(vlen), intent (in)    :: relvar
        1439             :
        1440             :   real(r8), dimension(vlen), intent (out)   :: au ! = prc autoconversion rate
        1441             :   real(r8), dimension(vlen), intent (out)   :: nprc1 ! = number tendency
        1442             :   real(r8), dimension(vlen), intent (out)   :: nprc ! = number tendency fixed size for rain
        1443             :
        1444             :   ! parameters for droplet mass spectral shape,
        1445             :   !used by Seifert and Beheng (2001)
        1446             :   ! warm rain scheme only (iparam = 1)
        1447             :   real(r8), parameter :: dnu(16) = [0._r8,-0.557_r8,-0.430_r8,-0.307_r8, &
        1448             :      -0.186_r8,-0.067_r8,0.050_r8,0.167_r8,0.282_r8,0.397_r8,0.512_r8, &
        1449             :      0.626_r8,0.739_r8,0.853_r8,0.966_r8,0.966_r8]
        1450             :
        1451             :   ! parameters for Seifert and Beheng (2001) autoconversion/accretion
        1452             :   real(r8), parameter :: kc = 9.44e9_r8
        1453             :   real(r8), parameter :: kr = 5.78e3_r8
        1454             :   real(r8) :: dum, dum1, nu
        1455             :   integer  :: dumi, i
        1456             :
        1457             :   !$acc data present (pgam,qc,nc,qr,rho,relvar,au,nprc1,nprc)
        1458             :
        1459             :   !$acc parallel vector_length(VLENS) default(present)
        1460             :   !$acc loop gang vector private(dumi,nu,dum,dum1)
        1461           0 :   do i=1,vlen
        1462           0 :      if (qc(i) > qsmall) then
        1463           0 :        dumi = int(pgam(i))
        1464           0 :        nu   = dnu(dumi)+(dnu(dumi+1)-dnu(dumi))* &
        1465           0 :                (pgam(i)-dumi)
        1466             :
        1467           0 :        dum  = 1._r8-qc(i)/(qc(i)+qr(i))
        1468           0 :        dum1 = 600._r8*dum**0.68_r8*(1._r8-dum**0.68_r8)**3
        1469             :
        1470             :        au(i) = kc/(20._r8*2.6e-7_r8)* &
        1471             :          (nu+2._r8)*(nu+4._r8)/(nu+1._r8)**2._r8* &
        1472             :          (rho(i)*qc(i)/1000._r8)**4._r8/(rho(i)*nc(i)/1.e6_r8)**2._r8* &
        1473           0 :          (1._r8+dum1/(1._r8-dum)**2)*1000._r8 / rho(i)
        1474             :
        1475           0 :        nprc1(i) = au(i)*2._r8/2.6e-7_r8*1000._r8
        1476           0 :        nprc(i)  = au(i)/droplet_mass_40um
        1477             :      else
        1478           0 :        au(i)    = 0._r8
        1479           0 :        nprc1(i) = 0._r8
        1480           0 :        nprc(i)  = 0._r8
        1481             :      end if
        1482             :   end do
        1483             :   !$acc end parallel
        1484             :
        1485             :   !$acc end data
        1486           0 : end subroutine sb2001v2_liq_autoconversion
        1487             :
        1488             : !========================================================================
        1489             : !SB2001 Accretion V2
        1490             :
        1491           0 : subroutine sb2001v2_accre_cld_water_rain(qc,nc,qr,rho,relvar,pra,npra,vlen)
        1492             :   !
        1493             :   ! ---------------------------------------------------------------------
        1494             :   ! ACCR_SB calculates the evolution of mass mxng-ratio due to accretion
        1495             :   ! and self collection following Seifert & Beheng (2001).
        1496             :   !
        1497             :
        1498             :   integer, intent(in) :: vlen
        1499             :
        1500             :   real(r8), dimension(vlen), intent (in)    :: qc  ! = qc (cld water mixing ratio)
        1501             :   real(r8), dimension(vlen), intent (in)    :: nc  ! = nc (cld water number conc /kg)
        1502             :   real(r8), dimension(vlen), intent (in)    :: qr  ! = qr (rain water mixing ratio)
        1503             :   real(r8), dimension(vlen), intent (in)    :: rho ! = rho : density profile
        1504             :   real(r8), dimension(vlen), intent (in)    :: relvar
        1505             :
        1506             :   ! Output tendencies
        1507             :   real(r8), dimension(vlen), intent(out) :: pra  ! MMR
        1508             :   real(r8), dimension(vlen), intent(out) :: npra ! Number
        1509             :
        1510             :   ! parameters for Seifert and Beheng (2001) autoconversion/accretion
        1511             :   real(r8), parameter :: kc = 9.44e9_r8
        1512             :   real(r8), parameter :: kr = 5.78e3_r8
        1513             :
        1514             :   real(r8) :: dum, dum1
        1515             :   integer :: i
        1516             :
        1517             :   ! accretion
        1518             :
        1519             :   !$acc data present (qc,nc,qr,rho,relvar,pra,npra)
        1520             :
        1521             :   !$acc parallel vector_length(VLENS) default(present)
        1522             :   !$acc loop gang vector private(dum,dum1)
        1523           0 :   do i = 1,vlen
        1524           0 :     if (qc(i) > qsmall) then
        1525           0 :       dum     = 1._r8-qc(i)/(qc(i)+qr(i))
        1526           0 :       dum1    = (dum/(dum+5.e-4_r8))**4._r8
        1527           0 :       pra(i)  = kr*rho(i)*0.001_r8*qc(i)*qr(i)*dum1
        1528             :       npra(i) = pra(i)*rho(i)*0.001_r8*(nc(i)*rho(i)*1.e-6_r8)/ &
        1529           0 :            (qc(i)*rho(i)*0.001_r8)*1.e6_r8 / rho(i)
        1530             :     else
        1531           0 :       pra(i)  = 0._r8
        1532           0 :       npra(i) = 0._r8
        1533             :     end if
        1534             :   end do
        1535             :   !$acc end parallel
        1536             :
        1537             :   !$acc end data
        1538           0 : end subroutine sb2001v2_accre_cld_water_rain
        1539             :
        1540             : !========================================================================
        1541             : ! Autoconversion of cloud ice to snow
        1542             : ! similar to Ferrier (1994)
        1543             :
        1544           0 : subroutine ice_autoconversion(t, qiic, lami, n0i, dcs, prci, nprci, vlen)
        1545             :
        1546             :   integer, intent(in) :: vlen
        1547             :   real(r8), dimension(vlen), intent(in) :: t
        1548             :   real(r8), dimension(vlen), intent(in) :: qiic
        1549             :   real(r8), dimension(vlen), intent(in) :: lami
        1550             :   real(r8), dimension(vlen), intent(in) :: n0i
        1551             :   real(r8),                    intent(in) :: dcs
        1552             :
        1553             :   real(r8), dimension(vlen), intent(out) :: prci
        1554             :   real(r8), dimension(vlen), intent(out) :: nprci
        1555             :
        1556             :   ! Assume autoconversion timescale of 180 seconds.
        1557             :   real(r8), parameter :: ac_time = 180._r8
        1558             :
        1559             :   ! Average mass of an ice particle.
        1560             :   real(r8) :: m_ip
        1561             :   ! Ratio of autoconversion diameter to average diameter.
        1562             :   real(r8) :: d_rat
        1563             :   integer :: i
        1564             :
        1565             :   !$acc data present (t,qiic,lami,n0i,prci,nprci)
        1566             :
        1567             :   !$acc parallel vector_length(VLENS) default(present)
        1568             :   !$acc loop gang vector private(d_rat,m_ip)
        1569           0 :   do i=1,vlen
        1570           0 :      if (t(i) <= tmelt .and. qiic(i) >= qsmall) then
        1571           0 :         d_rat = lami(i)*dcs
        1572             :
        1573             :         ! Rate of ice particle conversion (number).
        1574           0 :         nprci(i) = n0i(i)/(lami(i)*ac_time)*exp(-d_rat)
        1575           0 :         m_ip = (rhoi*pi/6._r8) / lami(i)**3
        1576             :
        1577             :         ! Rate of mass conversion.
        1578             :         ! Note that this is:
        1579             :         ! m n (d^3 + 3 d^2 + 6 d + 6)
        1580             :         prci(i) = m_ip * nprci(i) * &
        1581           0 :              (((d_rat + 3._r8)*d_rat + 6._r8)*d_rat + 6._r8)
        1582             :      else
        1583           0 :         prci(i)  = 0._r8
        1584           0 :         nprci(i) = 0._r8
        1585             :      end if
        1586             :   end do
        1587             :   !$acc end parallel
        1588             :
        1589             :   !$acc end data
        1590           0 : end subroutine ice_autoconversion
        1591             :
        1592             : ! immersion freezing (Bigg, 1953)
        1593             : !===================================
        1594             :
        1595           0 : subroutine immersion_freezing(microp_uniform, t, pgam, lamc, &
        1596           0 :      qcic, ncic, relvar, mnuccc, nnuccc, vlen)
        1597             :
        1598             :   integer, intent(in) :: vlen
        1599             :   logical, intent(in) :: microp_uniform
        1600             :
        1601             :   ! Temperature
        1602             :   real(r8), dimension(vlen), intent(in) :: t
        1603             :
        1604             :   ! Cloud droplet size distribution parameters
        1605             :   real(r8), dimension(vlen), intent(in) :: pgam
        1606             :   real(r8), dimension(vlen), intent(in) :: lamc
        1607             :
        1608             :   ! MMR and number concentration of in-cloud liquid water
        1609             :   real(r8), dimension(vlen), intent(in) :: qcic
        1610             :   real(r8), dimension(vlen), intent(in) :: ncic
        1611             :
        1612             :   ! Relative variance of cloud water
        1613             :   real(r8), dimension(vlen), intent(in) :: relvar
        1614             :
        1615             :   ! Output tendencies
        1616             :   real(r8), dimension(vlen), intent(out) :: mnuccc ! MMR
        1617             :   real(r8), dimension(vlen), intent(out) :: nnuccc ! Number
        1618             :
        1619             :   ! Coefficients that will be omitted for sub-columns
        1620           0 :   real(r8), dimension(vlen) :: dum
        1621             :   integer  :: i
        1622             :   real(r8) :: tmp
        1623             :
        1624             :   !$acc data present (t,pgam,lamc,qcic,ncic) &
        1625             :   !$acc      present (relvar,mnuccc,nnuccc)  &
        1626             :   !$acc      create  (dum)
        1627             :
        1628           0 :   if (.not. microp_uniform) then
        1629           0 :      call var_coef(relvar, 2, dum, vlen)
        1630             :   else
        1631             :      !$acc parallel vector_length(VLENS) default(present)
        1632             :      !$acc loop gang vector
        1633           0 :      do i =1,vlen
        1634           0 :         dum(i) = 1._r8
        1635             :      end do
        1636             :      !$acc end parallel
        1637             :   end if
        1638             :
        1639             :   !$acc parallel vector_length(VLENS) default(present)
        1640             :   !$acc loop gang vector private(tmp)
        1641           0 :   do i=1,vlen
        1642           0 :      if (qcic(i) >= qsmall .and. t(i) < 269.15_r8) then
        1643           0 :         call rising_factorial(pgam(i)+1._r8, 3, tmp)
        1644             :         nnuccc(i) = &
        1645             :              pi/6._r8*ncic(i)*tmp* &
        1646           0 :              bimm*(exp(aimm*(tmelt - t(i)))-1._r8)/lamc(i)**3
        1647             :
        1648           0 :         call rising_factorial(pgam(i)+4._r8, 3, tmp)
        1649             :         mnuccc(i) = dum(i) * nnuccc(i) * &
        1650             :              pi/6._r8*rhow* &
        1651           0 :              tmp/lamc(i)**3
        1652             :      else
        1653           0 :         mnuccc(i) = 0._r8
        1654           0 :         nnuccc(i) = 0._r8
        1655             :      end if ! qcic > qsmall and t < 4 deg C
        1656             :   end do
        1657             :   !$acc end parallel
        1658             :
        1659             :   !$acc end data
        1660           0 : end subroutine immersion_freezing
        1661             :
        1662             : ! contact freezing (-40<T<-3 C) (Young, 1974) with hooks into simulated dust
        1663             : !===================================================================
        1664             : ! dust size and number in multiple bins are read in from companion routine
        1665             :
        1666           0 : subroutine contact_freezing (microp_uniform, t, p, rndst, nacon, &
        1667           0 :      pgam, lamc, qcic, ncic, relvar, mnucct, nnucct, vlen, mdust)
        1668             :
        1669             :   logical, intent(in) :: microp_uniform
        1670             :
        1671             :   integer, intent(in) :: vlen
        1672             :   integer, intent(in) :: mdust
        1673             :
        1674             :   real(r8), dimension(vlen), intent(in) :: t            ! Temperature
        1675             :   real(r8), dimension(vlen), intent(in) :: p            ! Pressure
        1676             :   real(r8), dimension(vlen, mdust), intent(in) :: rndst ! Radius (for multiple dust bins)
        1677             :   real(r8), dimension(vlen, mdust), intent(in) :: nacon ! Number (for multiple dust bins)
        1678             :
        1679             :   ! Size distribution parameters for cloud droplets
        1680             :   real(r8), dimension(vlen), intent(in) :: pgam
        1681             :   real(r8), dimension(vlen), intent(in) :: lamc
        1682             :
        1683             :   ! MMR and number concentration of in-cloud liquid water
        1684             :   real(r8), dimension(vlen), intent(in) :: qcic
        1685             :   real(r8), dimension(vlen), intent(in) :: ncic
        1686             :
        1687             :   ! Relative cloud water variance
        1688             :   real(r8), dimension(vlen), intent(in) :: relvar
        1689             :
        1690             :   ! Output tendencies
        1691             :   real(r8), dimension(vlen), intent(out) :: mnucct ! MMR
        1692             :   real(r8), dimension(vlen), intent(out) :: nnucct ! Number
        1693             :
        1694             :   real(r8) :: tcnt                  ! scaled relative temperature
        1695             :   real(r8) :: viscosity             ! temperature-specific viscosity (kg/m/s)
        1696             :   real(r8) :: mfp                   ! temperature-specific mean free path (m)
        1697             :
        1698             :   ! Dimension these according to number of dust bins, inferred from rndst size
        1699           0 :   real(r8) :: nslip(size(rndst,2))  ! slip correction factors
        1700           0 :   real(r8) :: ndfaer(size(rndst,2)) ! aerosol diffusivities (m^2/sec)
        1701             :
        1702             :   ! Coefficients not used for subcolumns
        1703             :   real(r8) :: dum, dum1, tmp
        1704             :
        1705             :   ! Common factor between mass and number.
        1706             :   real(r8) :: contact_factor
        1707             :
        1708             :   integer  :: i, j
        1709             :
        1710             :   !$acc data present (t,p,rndst,nacon,pgam,lamc) &
        1711             :   !$acc      present (qcic,ncic,relvar,mnucct,nnucct) &
        1712             :   !$acc      create  (nslip,ndfaer)
        1713             :
        1714             :   !$acc parallel vector_length(VLENS) default(present)
        1715             :   !$acc loop gang vector private(dum,dum1,tcnt,viscosity,mfp,contact_factor,tmp)
        1716           0 :   do i = 1,vlen
        1717           0 :      if (qcic(i) >= qsmall .and. t(i) < 269.15_r8) then
        1718           0 :         if (microp_uniform) then
        1719           0 :            dum  = 1._r8
        1720           0 :            dum1 = 1._r8
        1721             :         else
        1722           0 :            call var_coef(relvar(i), 4._r8/3._r8, dum)
        1723           0 :            call var_coef(relvar(i), 1._r8/3._r8, dum1)
        1724             :         end if
        1725             :
        1726           0 :         tcnt=(270.16_r8-t(i))**1.3_r8
        1727           0 :         viscosity = 1.8e-5_r8*(t(i)/298.0_r8)**0.85_r8    ! Viscosity (kg/m/s)
        1728             :         mfp = 2.0_r8*viscosity/ &                         ! Mean free path (m)
        1729           0 :                      (p(i)*sqrt( 8.0_r8*28.96e-3_r8/(pi*8.314409_r8*t(i)) ))
        1730             :
        1731           0 :         do j = 1, mdust
        1732             :            ! Note that these two are vectors.
        1733           0 :            nslip(j) = 1.0_r8+(mfp/rndst(i,j))*(1.257_r8+(0.4_r8*exp(-(1.1_r8*rndst(i,j)/mfp))))! Slip correction factor
        1734             :
        1735           0 :            ndfaer(j) = 1.381e-23_r8*t(i)*nslip(j)/(6._r8*pi*viscosity*rndst(i,j))  ! aerosol diffusivity (m2/s)
        1736             :         end do
        1737             :
        1738           0 :         contact_factor = dot_product(ndfaer,nacon(i,:)*tcnt) * pi * &
        1739           0 :                                      ncic(i) * (pgam(i) + 1._r8) / lamc(i)
        1740           0 :         call rising_factorial(pgam(i)+2._r8, 3, tmp)
        1741           0 :         mnucct(i) = dum * contact_factor * pi/3._r8*rhow*tmp/lamc(i)**3
        1742           0 :         nnucct(i) = dum1 * 2._r8 * contact_factor
        1743             :      else
        1744           0 :         mnucct(i) = 0._r8
        1745           0 :         nnucct(i) = 0._r8
        1746             :      end if ! qcic > qsmall and t < 4 deg C
        1747             :   end do
        1748             :   !$acc end parallel
        1749             :
        1750             :   !$acc end data
        1751           0 : end subroutine contact_freezing
        1752             :
        1753             : ! snow self-aggregation from passarelli, 1978, used by reisner, 1998
        1754             : !===================================================================
        1755             : ! this is hard-wired for bs = 0.4 for now
        1756             : ! ignore self-collection of cloud ice
        1757             :
        1758           0 : subroutine snow_self_aggregation(t, rho, asn, rhosn, qsic, nsic, nsagg, vlen)
        1759             :
        1760             :   integer,                   intent(in) :: vlen
        1761             :
        1762             :   real(r8), dimension(vlen), intent(in) :: t     ! Temperature
        1763             :   real(r8), dimension(vlen), intent(in) :: rho   ! Density
        1764             :   real(r8), dimension(vlen), intent(in) :: asn   ! fall speed parameter for snow
        1765             :   real(r8),                  intent(in) :: rhosn ! density of snow
        1766             :
        1767             :   ! In-cloud snow
        1768             :   real(r8), dimension(vlen), intent(in) :: qsic ! MMR
        1769             :   real(r8), dimension(vlen), intent(in) :: nsic ! Number
        1770             :
        1771             :   ! Output number tendency
        1772             :   real(r8), dimension(vlen), intent(out) :: nsagg
        1773             :
        1774             :   integer :: i
        1775             :
        1776             :   !$acc data present (t,rho,asn,qsic,nsic,nsagg)
        1777             :
        1778             :   !$acc parallel vector_length(VLENS) default(present)
        1779             :   !$acc loop gang vector
        1780           0 :   do i=1,vlen
        1781           0 :      if (qsic(i) >= qsmall .and. t(i) <= tmelt) then
        1782             :         nsagg(i) = -1108._r8*eii/(4._r8*720._r8*rhosn)*asn(i)*qsic(i)*nsic(i)*rho(i)*&
        1783           0 :              ((qsic(i)/nsic(i))*(1._r8/(rhosn*pi)))**((bs-1._r8)/3._r8)
        1784             :      else
        1785           0 :         nsagg(i) = 0._r8
        1786             :      end if
        1787             :   end do
        1788             :   !$acc end parallel
        1789             :
        1790             :   !$acc end data
        1791           0 : end subroutine snow_self_aggregation
        1792             :
        1793             : ! ice self-aggregation used in mg4
        1794             : !===================================================================
        1795             :
        1796           0 : subroutine ice_self_aggregation(t, rho, rhof, af1pr3, qiic, nsagg, vlen)
        1797             :
        1798             :   integer,                   intent(in) :: vlen
        1799             :
        1800             :   real(r8), dimension(vlen), intent(in) :: t        ! Temperature
        1801             :   real(r8), dimension(vlen), intent(in) :: rho      ! Density
        1802             :   real(r8), dimension(vlen), intent(in) :: rhof     ! Density correction
        1803             :
        1804             :   real(r8), dimension(vlen), intent(in) :: af1pr3   ! process rate
        1805             :
        1806             :   ! In-cloud ice
        1807             :   real(r8), dimension(vlen), intent(in) :: qiic ! MMR
        1808             :
        1809             :   ! Collection efficiency
        1810             :   real(r8), parameter :: eii = 0.1_r8
        1811             :
        1812             :   ! Output number tendency
        1813             :   real(r8), dimension(vlen), intent(out) :: nsagg
        1814             :
        1815             :   integer :: i
        1816             :
        1817             :   !$acc data present (t,rho,rhof,af1pr3,qiic,nsagg)
        1818             :
        1819             :   !$acc parallel vector_length(VLENS) default(present)
        1820             :   !$acc loop gang vector
        1821           0 :   do i=1,vlen
        1822           0 :      if (qiic(i) >= qsmall .and. t(i) <= tmelt) then
        1823           0 :         nsagg(i) = af1pr3(i)*rho(i)*eii*rhof(i)
        1824             :      else
        1825           0 :         nsagg(i)=0._r8
        1826             :      end if
        1827             :   end do
        1828             :   !$acc end parallel
        1829             :
        1830             :   !$acc end data
        1831           0 : end subroutine ice_self_aggregation
        1832             :
        1833             : ! accretion of cloud droplets onto snow/graupel
        1834             : !===================================================================
        1835             : ! here use continuous collection equation with
        1836             : ! simple gravitational collection kernel
        1837             : ! ignore collisions between droplets/cloud ice
        1838             : ! since minimum size ice particle for accretion is 50 - 150 micron
        1839             :
        1840           0 : subroutine accrete_cloud_water_snow(t, rho, asn, uns, mu, qcic, ncic, qsic, &
        1841           0 :      pgam, lamc, lams, n0s, psacws, npsacws, vlen)
        1842             :
        1843             :   integer,                   intent(in) :: vlen
        1844             :   real(r8), dimension(vlen), intent(in) :: t   ! Temperature
        1845             :   real(r8), dimension(vlen), intent(in) :: rho ! Density
        1846             :   real(r8), dimension(vlen), intent(in) :: asn ! Fallspeed parameter (snow)
        1847             :   real(r8), dimension(vlen), intent(in) :: uns ! Current fallspeed   (snow)
        1848             :   real(r8), dimension(vlen), intent(in) :: mu  ! Viscosity
        1849             :
        1850             :   ! In-cloud liquid water
        1851             :   real(r8), dimension(vlen), intent(in) :: qcic ! MMR
        1852             :   real(r8), dimension(vlen), intent(in) :: ncic ! Number
        1853             :
        1854             :   ! In-cloud snow
        1855             :   real(r8), dimension(vlen), intent(in) :: qsic ! MMR
        1856             :
        1857             :   ! Cloud droplet size parameters
        1858             :   real(r8), dimension(vlen), intent(in) :: pgam
        1859             :   real(r8), dimension(vlen), intent(in) :: lamc
        1860             :
        1861             :   ! Snow size parameters
        1862             :   real(r8), dimension(vlen), intent(in) :: lams
        1863             :   real(r8), dimension(vlen), intent(in) :: n0s
        1864             :
        1865             :   ! Output tendencies
        1866             :   real(r8), dimension(vlen), intent(out) :: psacws  ! Mass mixing ratio
        1867             :   real(r8), dimension(vlen), intent(out) :: npsacws ! Number concentration
        1868             :
        1869             :   real(r8) :: dc0 ! Provisional mean droplet size
        1870             :   real(r8) :: dum
        1871             :   real(r8) :: eci ! collection efficiency for riming of snow by droplets
        1872             :
        1873             :   ! Fraction of cloud droplets accreted per second
        1874             :   real(r8) :: accrete_rate
        1875             :   integer :: i
        1876             :
        1877             :   ! ignore collision of snow with droplets above freezing
        1878             :
        1879             :   !$acc data present (t,rho,asn,uns,mu,qcic,ncic,qsic) &
        1880             :   !$acc      present (pgam,lamc,lams,n0s,psacws,npsacws)
        1881             :
        1882             :   !$acc parallel vector_length(VLENS) default(present)
        1883             :   !$acc loop gang vector private(dc0,dum,eci,accrete_rate)
        1884           0 :   do i=1,vlen
        1885           0 :      if (qsic(i) >= qsmall .and. t(i) <= tmelt .and. qcic(i) >= qsmall) then
        1886             :         ! put in size dependent collection efficiency
        1887             :         ! mean diameter of snow is area-weighted, since
        1888             :         ! accretion is function of crystal geometric area
        1889             :         ! collection efficiency is approximation based on stoke's law (Thompson et al. 2004)
        1890           0 :         dc0 = (pgam(i)+1._r8)/lamc(i)
        1891           0 :         dum = dc0*dc0*uns(i)*rhow*lams(i)/(9._r8*mu(i))
        1892           0 :         eci = dum*dum/((dum+0.4_r8)*(dum+0.4_r8))
        1893           0 :         eci = max(eci,0._r8)
        1894           0 :         eci = min(eci,1._r8)
        1895             :
        1896             :         ! no impact of sub-grid distribution of qc since psacws
        1897             :         ! is linear in qc
        1898           0 :         accrete_rate = pi/4._r8*asn(i)*rho(i)*n0s(i)*eci*gamma_bs_plus3 / lams(i)**(bs+3._r8)
        1899           0 :         psacws(i)  = accrete_rate*qcic(i)
        1900           0 :         npsacws(i) = accrete_rate*ncic(i)
        1901             :      else
        1902           0 :         psacws(i)  = 0._r8
        1903           0 :         npsacws(i) = 0._r8
        1904             :      end if
        1905             :   end do
        1906             :   !$acc end parallel
        1907             :
        1908             :   !$acc end data
        1909           0 : end subroutine accrete_cloud_water_snow
        1910             :
        1911             : ! Version used in MG4
        1912           0 : subroutine accrete_cloud_water_ice(t, rho, rhof, af1pr4, qcic, ncic, qiic, psacws, npsacws, vlen)
        1913             :
        1914             :   integer,                   intent(in) :: vlen
        1915             :   real(r8), dimension(vlen), intent(in) :: t        ! Temperature
        1916             :   real(r8), dimension(vlen), intent(in) :: rho      ! Density
        1917             :   real(r8), dimension(vlen), intent(in) :: rhof     ! Density correction factor
        1918             :   real(r8), dimension(vlen), intent(in) :: af1pr4   ! Process rate from look-up table
        1919             :   ! In-cloud liquid water
        1920             :   real(r8), dimension(vlen), intent(in) :: qcic     ! MMR
        1921             :   real(r8), dimension(vlen), intent(in) :: ncic     ! Number
        1922             :   ! In-cloud ice
        1923             :   real(r8), dimension(vlen), intent(in) :: qiic     ! MMR
        1924             :   ! Output tendencies
        1925             :   real(r8), dimension(vlen), intent(out) :: psacws  ! Mass mixing ratio
        1926             :   real(r8), dimension(vlen), intent(out) :: npsacws ! Number concentration
        1927             :
        1928             :   real(r8) :: eci     ! collection efficiency for riming of snow by droplets
        1929             :   integer  :: i
        1930             :
        1931             :   ! ignore collision of snow with droplets above freezing
        1932           0 :   eci = 0.5_r8 ! from p3 program
        1933             :
        1934             :   !$acc data present (t,rho,rhof,af1pr4,qcic) &
        1935             :   !$acc      present (ncic,qiic,psacws,npsacws)
        1936             :
        1937             :   !$acc parallel vector_length(VLENS) default(present)
        1938             :   !$acc loop gang vector
        1939           0 :   do i=1,vlen
        1940           0 :      if (qiic(i) >= qsmall .and. t(i) <= tmelt .and. qcic(i) >= qsmall) then
        1941           0 :          psacws(i)  = rhof(i)*af1pr4(i)*qcic(i)*eci*rho(i)
        1942           0 :          npsacws(i) = rhof(i)*af1pr4(i)*ncic(i)*eci*rho(i)
        1943             :      else
        1944           0 :          psacws(i)  = 0._r8
        1945           0 :          npsacws(i) = 0._r8
        1946             :      end if
        1947             :   end do
        1948             :   !$acc end parallel
        1949             :
        1950             :   !$acc end data
        1951           0 : end subroutine accrete_cloud_water_ice
        1952             :
        1953             : ! add secondary ice production due to accretion of droplets by snow
        1954             : !===================================================================
        1955             : ! (Hallet-Mossop process) (from Cotton et al., 1986)
        1956             :
        1957           0 : subroutine secondary_ice_production(t, psacws, msacwi, nsacwi, vlen)
        1958             :
        1959             :   integer,                   intent(in) :: vlen
        1960             :   real(r8), dimension(vlen), intent(in) :: t ! Temperature
        1961             :
        1962             :   ! Accretion of cloud water to snow tendencies
        1963             :   real(r8), dimension(vlen), intent(inout) :: psacws ! MMR
        1964             :
        1965             :   ! Output (ice) tendencies
        1966             :   real(r8), dimension(vlen), intent(out) :: msacwi ! MMR
        1967             :   real(r8), dimension(vlen), intent(out) :: nsacwi ! Number
        1968             :   integer :: i
        1969             :
        1970             :   !$acc data present (t,psacws,msacwi,nsacwi)
        1971             :
        1972             :   !$acc parallel vector_length(VLENS) default(present)
        1973             :   !$acc loop gang vector
        1974           0 :   do i=1,vlen
        1975           0 :      if((t(i) < 270.16_r8) .and. (t(i) >= 268.16_r8)) then
        1976           0 :         nsacwi(i) = 3.5e8_r8*(270.16_r8-t(i))/2.0_r8*psacws(i)
        1977           0 :      else if((t(i) < 268.16_r8) .and. (t(i) >= 265.16_r8)) then
        1978           0 :         nsacwi(i) = 3.5e8_r8*(t(i)-265.16_r8)/3.0_r8*psacws(i)
        1979             :      else
        1980           0 :         nsacwi(i) = 0.0_r8
        1981             :      endif
        1982             :   end do
        1983             :
        1984             :   !$acc loop gang vector
        1985           0 :   do i=1,vlen
        1986           0 :      msacwi(i) = min(nsacwi(i)*mi0, psacws(i))
        1987           0 :      psacws(i) = psacws(i) - msacwi(i)
        1988             :   end do
        1989             :   !$acc end parallel
        1990             :
        1991             :   !$acc end data
        1992           0 : end subroutine secondary_ice_production
        1993             :
        1994             : ! accretion of rain water by snow
        1995             : !===================================================================
        1996             : ! formula from ikawa and saito, 1991, used by reisner et al., 1998
        1997             :
        1998           0 : subroutine accrete_rain_snow(t, rho, umr, ums, unr, uns, qric, qsic, &
        1999           0 :      lamr, n0r, lams, n0s, pracs, npracs, vlen)
        2000             :
        2001             :   integer,                   intent(in) :: vlen
        2002             :   real(r8), dimension(vlen), intent(in) :: t   ! Temperature
        2003             :   real(r8), dimension(vlen), intent(in) :: rho ! Density
        2004             :   ! Fallspeeds
        2005             :   ! mass-weighted
        2006             :   real(r8), dimension(vlen), intent(in) :: umr ! rain
        2007             :   real(r8), dimension(vlen), intent(in) :: ums ! snow
        2008             :   ! number-weighted
        2009             :   real(r8), dimension(vlen), intent(in) :: unr ! rain
        2010             :   real(r8), dimension(vlen), intent(in) :: uns ! snow
        2011             :   ! In cloud MMRs
        2012             :   real(r8), dimension(vlen), intent(in) :: qric ! rain
        2013             :   real(r8), dimension(vlen), intent(in) :: qsic ! snow
        2014             :   ! Size distribution parameters
        2015             :   ! rain
        2016             :   real(r8), dimension(vlen), intent(in) :: lamr
        2017             :   real(r8), dimension(vlen), intent(in) :: n0r
        2018             :   ! snow
        2019             :   real(r8), dimension(vlen), intent(in) :: lams
        2020             :   real(r8), dimension(vlen), intent(in) :: n0s
        2021             :   ! Output tendencies
        2022             :   real(r8), dimension(vlen), intent(out) :: pracs  ! MMR
        2023             :   real(r8), dimension(vlen), intent(out) :: npracs ! Number
        2024             :   ! Collection efficiency for accretion of rain by snow
        2025             :   real(r8), parameter :: ecr = 1.0_r8
        2026             :   ! Ratio of average snow diameter to average rain diameter.
        2027             :   real(r8) :: d_rat
        2028             :   ! Common factor between mass and number expressions
        2029             :   real(r8) :: common_factor
        2030             :   integer  :: i
        2031             :
        2032             :   !$acc data present (t,rho,umr,ums,unr,uns,qric,qsic) &
        2033             :   !$acc      present (lamr,n0r,lams,n0s,pracs,npracs)
        2034             :
        2035             :   !$acc parallel vector_length(VLENS) default(present)
        2036             :   !$acc loop gang vector private(common_factor,d_rat)
        2037           0 :   do i=1,vlen
        2038           0 :      if (qric(i) >= icsmall .and. qsic(i) >= icsmall .and. t(i) <= tmelt) then
        2039           0 :         common_factor = pi*ecr*rho(i)*n0r(i)*n0s(i)/(lamr(i)**3 * lams(i))
        2040           0 :         d_rat = lamr(i)/lams(i)
        2041             :         pracs(i) = common_factor*pi*rhow* &
        2042             :              sqrt((1.2_r8*umr(i)-0.95_r8*ums(i))**2 + 0.08_r8*ums(i)*umr(i)) / lamr(i)**3 * &
        2043           0 :              ((0.5_r8*d_rat + 2._r8)*d_rat + 5._r8)
        2044             :         npracs(i) = common_factor*0.5_r8* &
        2045             :              sqrt(1.7_r8*(unr(i)-uns(i))**2 + 0.3_r8*unr(i)*uns(i)) * &
        2046           0 :              ((d_rat + 1._r8)*d_rat + 1._r8)
        2047             :      else
        2048           0 :         pracs(i) = 0._r8
        2049           0 :         npracs(i) = 0._r8
        2050             :      end if
        2051             :   end do
        2052             :   !$acc end parallel
        2053             :
        2054             :   !$acc end data
        2055           0 : end subroutine accrete_rain_snow
        2056             :
        2057             : ! Version used in MG4
        2058           0 : subroutine accrete_rain_ice(t, rho, rhof, af1pr8, af1pr7, qric, qiic, &
        2059           0 :       n0r, pracs, npracs, vlen )
        2060             :
        2061             :   integer,                   intent(in) :: vlen
        2062             :   real(r8), dimension(vlen), intent(in) :: rhof ! Density correction factor
        2063             :   real(r8), dimension(vlen), intent(in) :: rho  ! Density
        2064             :   real(r8), dimension(vlen), intent(in) :: t    ! Temperature
        2065             :   ! In cloud MMRs
        2066             :   real(r8), dimension(vlen), intent(in) :: qric ! rain
        2067             :   real(r8), dimension(vlen), intent(in) :: qiic ! combined snow and ice
        2068             :   ! Size distribution parameters
        2069             :   ! rain
        2070             :   real(r8), dimension(vlen), intent(in) :: n0r
        2071             :   ! Processrates
        2072             :   real(r8), dimension(vlen), intent(in) :: af1pr8
        2073             :   real(r8), dimension(vlen), intent(in) :: af1pr7
        2074             :   ! Output tendencies
        2075             :   real(r8), dimension(vlen), intent(out) :: pracs  ! MMR
        2076             :   real(r8), dimension(vlen), intent(out) :: npracs ! Number
        2077             :   ! Collection efficiency for accretion of rain by snow
        2078             :   real(r8), parameter :: ecr = 1.0_r8
        2079             :   integer :: i
        2080             :
        2081             :   !$acc data present (t,rho,rhof,qric,qiic,n0r) &
        2082             :   !$acc      present (af1pr8,af1pr7,pracs,npracs)
        2083             :
        2084             :   !$acc parallel vector_length(VLENS) default(present)
        2085             :   !$acc loop gang vector
        2086           0 :   do i=1,vlen
        2087           0 :      if (qric(i) >= icsmall .and. qiic(i) >= icsmall .and. t(i) <= tmelt) then
        2088           0 :         pracs(i)  = 10._r8**(af1pr8(i)+dlog10(n0r(i)))*rho(i)*rhof(i)*ecr
        2089           0 :         npracs(i) = 10._r8**(af1pr7(i)+dlog10(n0r(i)))*rho(i)*rhof(i)*ecr
        2090             :      else
        2091           0 :         pracs(i)  = 0._r8
        2092           0 :         npracs(i) = 0._r8
        2093             :      end if
        2094             :   end do
        2095             :   !$acc end parallel
        2096             :
        2097             :   !$acc end data
        2098           0 : end subroutine accrete_rain_ice
        2099             :
        2100             : ! heterogeneous freezing of rain drops
        2101             : !===================================================================
        2102             : ! follows from Bigg (1953)
        2103             :
        2104           0 : subroutine heterogeneous_rain_freezing(t, qric, nric, lamr, mnuccr, nnuccr, vlen)
        2105             :
        2106             :   integer,                   intent(in) :: vlen
        2107             :   real(r8), dimension(vlen), intent(in) :: t    ! Temperature
        2108             :   ! In-cloud rain
        2109             :   real(r8), dimension(vlen), intent(in) :: qric ! MMR
        2110             :   real(r8), dimension(vlen), intent(in) :: nric ! Number
        2111             :   real(r8), dimension(vlen), intent(in) :: lamr ! size parameter
        2112             :   ! Output tendencies
        2113             :   real(r8), dimension(vlen), intent(out) :: mnuccr ! MMR
        2114             :   real(r8), dimension(vlen), intent(out) :: nnuccr ! Number
        2115             :   integer :: i
        2116             :
        2117             :   !$acc data present (t,qric,nric,lamr,mnuccr,nnuccr)
        2118             :
        2119             :   !$acc parallel vector_length(VLENS) default(present)
        2120             :   !$acc loop gang vector
        2121           0 :   do i=1,vlen
        2122           0 :      if (t(i) < 269.15_r8 .and. qric(i) >= qsmall) then
        2123             :         nnuccr(i) = pi*nric(i)*bimm* &
        2124           0 :              (exp(aimm*(tmelt - t(i)))-1._r8)/lamr(i)**3
        2125           0 :         mnuccr(i) = nnuccr(i) * 20._r8*pi*rhow/lamr(i)**3
        2126             :      else
        2127           0 :         mnuccr(i) = 0._r8
        2128           0 :         nnuccr(i) = 0._r8
        2129             :      end if
        2130             :   end do
        2131             :   !$acc end parallel
        2132             :
        2133             :   !$acc end data
        2134           0 : end subroutine heterogeneous_rain_freezing
        2135             :
        2136             : ! accretion of cloud liquid water by rain
        2137             : !===================================================================
        2138             : ! formula from Khrouditnov and Kogan (2000)
        2139             : ! gravitational collection kernel, droplet fall speed neglected
        2140             :
        2141           0 : subroutine accrete_cloud_water_rain(microp_uniform, qric, qcic, &
        2142           0 :      ncic, relvar, accre_enhan, pra, npra, vlen)
        2143             :
        2144             :   logical, intent(in) :: microp_uniform
        2145             :   integer, intent(in) :: vlen
        2146             :   ! In-cloud rain
        2147             :   real(r8), dimension(vlen), intent(in) :: qric ! MMR
        2148             :   ! Cloud droplets
        2149             :   real(r8), dimension(vlen), intent(in) :: qcic ! MMR
        2150             :   real(r8), dimension(vlen), intent(in) :: ncic ! Number
        2151             :   ! SGS variability
        2152             :   real(r8), dimension(vlen), intent(in) :: relvar
        2153             :   real(r8), dimension(vlen), intent(in) :: accre_enhan
        2154             :   ! Output tendencies
        2155             :   real(r8), dimension(vlen), intent(out) :: pra  ! MMR
        2156             :   real(r8), dimension(vlen), intent(out) :: npra ! Number
        2157             :   ! Coefficient that varies for subcolumns
        2158           0 :   real(r8), dimension(vlen) :: pra_coef
        2159             :   integer :: i
        2160             :
        2161             :   !$acc data present (qric,qcic,ncic,relvar,accre_enhan,pra,npra) &
        2162             :   !$acc      create  (pra_coef)
        2163             :
        2164           0 :   if (.not. microp_uniform) then
        2165           0 :      call var_coef(relvar, 1.15_r8, pra_coef, vlen)
        2166             :      !$acc parallel vector_length(VLENS) default(present)
        2167             :      !$acc loop gang vector
        2168           0 :      do i = 1, vlen
        2169           0 :         pra_coef(i) = accre_enhan(i) * pra_coef(i)
        2170             :      end do
        2171             :      !$acc end parallel
        2172             :   else
        2173             :      !$acc parallel vector_length(VLENS) default(present)
        2174             :      !$acc loop gang vector
        2175           0 :      do i = 1, vlen
        2176           0 :         pra_coef(i) = 1._r8
        2177             :      end do
        2178             :      !$acc end parallel
        2179             :   end if
        2180             :
        2181             :   !$acc parallel vector_length(VLENS) default(present)
        2182             :   !$acc loop gang vector
        2183           0 :   do i=1,vlen
        2184           0 :      if (qric(i) >= qsmall .and. qcic(i) >= qsmall) then
        2185             :         ! include sub-grid distribution of cloud water
        2186           0 :         pra(i)  = pra_coef(i) * 67._r8*(qcic(i)*qric(i))**1.15_r8
        2187           0 :         npra(i) = pra(i)*ncic(i)/qcic(i)
        2188             :      else
        2189           0 :         pra(i)  = 0._r8
        2190           0 :         npra(i) = 0._r8
        2191             :      end if
        2192             :   end do
        2193             :   !$acc end parallel
        2194             :
        2195             :   !$acc end data
        2196           0 : end subroutine accrete_cloud_water_rain
        2197             :
        2198             : ! Self-collection of rain drops
        2199             : !===================================================================
        2200             : ! from Beheng(1994)
        2201             :
        2202           0 : subroutine self_collection_rain(rho, qric, nric, nragg, vlen)
        2203             :
        2204             :   integer,                   intent(in) :: vlen
        2205             :   real(r8), dimension(vlen), intent(in) :: rho  ! Air density
        2206             :   ! Rain
        2207             :   real(r8), dimension(vlen), intent(in) :: qric ! MMR
        2208             :   real(r8), dimension(vlen), intent(in) :: nric ! Number
        2209             :   ! Output number tendency
        2210             :   real(r8), dimension(vlen), intent(out) :: nragg
        2211             :   integer :: i
        2212             :
        2213             :   !$acc data present (rho,qric,nric,nragg)
        2214             :
        2215             :   !$acc parallel vector_length(VLENS) default(present)
        2216             :   !$acc loop gang vector
        2217           0 :   do i=1,vlen
        2218           0 :      if (qric(i) >= qsmall) then
        2219           0 :         nragg(i) = -8._r8*nric(i)*qric(i)*rho(i)
        2220             :      else
        2221           0 :         nragg(i) = 0._r8
        2222             :      end if
        2223             :   end do
        2224             :   !$acc end parallel
        2225             :
        2226             :   !$acc end data
        2227           0 : end subroutine self_collection_rain
        2228             :
        2229             : ! Accretion of cloud ice by snow
        2230             : !===================================================================
        2231             : ! For this calculation, it is assumed that the Vs >> Vi
        2232             : ! and Ds >> Di for continuous collection
        2233             :
        2234           0 : subroutine accrete_cloud_ice_snow(t, rho, asn, qiic, niic, qsic, &
        2235           0 :      lams, n0s, prai, nprai, vlen)
        2236             :
        2237             :   integer,                   intent(in) :: vlen
        2238             :   real(r8), dimension(vlen), intent(in) :: t    ! Temperature
        2239             :   real(r8), dimension(vlen), intent(in) :: rho   ! Density
        2240             :   real(r8), dimension(vlen), intent(in) :: asn  ! Snow fallspeed parameter
        2241             :   ! Cloud ice
        2242             :   real(r8), dimension(vlen), intent(in) :: qiic ! MMR
        2243             :   real(r8), dimension(vlen), intent(in) :: niic ! Number
        2244             :   real(r8), dimension(vlen), intent(in) :: qsic ! Snow MMR
        2245             :   ! Snow size parameters
        2246             :   real(r8), dimension(vlen), intent(in) :: lams
        2247             :   real(r8), dimension(vlen), intent(in) :: n0s
        2248             :   ! Output tendencies
        2249             :   real(r8), dimension(vlen), intent(out) :: prai ! MMR
        2250             :   real(r8), dimension(vlen), intent(out) :: nprai ! Number
        2251             :   ! Fraction of cloud ice particles accreted per second
        2252             :   real(r8) :: accrete_rate
        2253             :   integer  :: i
        2254             :
        2255             :   !$acc data present (t,rho,asn,qiic,niic,qsic,lams,n0s,prai,nprai)
        2256             :
        2257             :   !$acc parallel vector_length(VLENS) default(present)
        2258             :   !$acc loop gang vector private(accrete_rate)
        2259           0 :   do i=1,vlen
        2260           0 :      if (qsic(i) >= qsmall .and. qiic(i) >= qsmall .and. t(i) <= tmelt) then
        2261             :         accrete_rate = pi/4._r8 * eii * asn(i) * rho(i) * n0s(i) * gamma_bs_plus3/ &
        2262           0 :              lams(i)**(bs+3._r8)
        2263           0 :         prai(i)  = accrete_rate * qiic(i)
        2264           0 :         nprai(i) = accrete_rate * niic(i)
        2265             :      else
        2266           0 :         prai(i)  = 0._r8
        2267           0 :         nprai(i) = 0._r8
        2268             :      end if
        2269             :   end do
        2270             :   !$acc end parallel
        2271             :
        2272             :   !$acc end data
        2273           0 : end subroutine accrete_cloud_ice_snow
        2274             :
        2275             : ! calculate evaporation/sublimation of rain and snow
        2276             : !===================================================================
        2277             : ! note: evaporation/sublimation occurs only in cloud-free portion of grid cell
        2278             : ! in-cloud condensation/deposition of rain and snow is neglected
        2279             : ! except for transfer of cloud water to snow through bergeron process
        2280             :
        2281           0 : subroutine evaporate_sublimate_precip(t, rho, dv, mu, sc, q, qvl, qvi, &
        2282           0 :      lcldm, precip_frac, arn, asn, qcic, qiic, qric, qsic, lamr, n0r, lams, n0s, &
        2283           0 :      pre, prds, am_evp_st, vlen, evap_rhthrsh)
        2284             :
        2285             :   integer,                   intent(in) :: vlen
        2286             :   real(r8), dimension(vlen), intent(in) :: t           ! temperature
        2287             :   real(r8), dimension(vlen), intent(in) :: rho         ! air density
        2288             :   real(r8), dimension(vlen), intent(in) :: dv          ! water vapor diffusivity
        2289             :   real(r8), dimension(vlen), intent(in) :: mu          ! viscosity
        2290             :   real(r8), dimension(vlen), intent(in) :: sc          ! schmidt number
        2291             :   real(r8), dimension(vlen), intent(in) :: q           ! humidity
        2292             :   real(r8), dimension(vlen), intent(in) :: qvl         ! saturation humidity (water)
        2293             :   real(r8), dimension(vlen), intent(in) :: qvi         ! saturation humidity (ice)
        2294             :   real(r8), dimension(vlen), intent(in) :: lcldm       ! liquid cloud fraction
        2295             :   real(r8), dimension(vlen), intent(in) :: precip_frac ! precipitation fraction (maximum overlap)
        2296             :   ! fallspeed parameters
        2297             :   real(r8), dimension(vlen), intent(in) :: arn         ! rain
        2298             :   real(r8), dimension(vlen), intent(in) :: asn         ! snow
        2299             :   ! In-cloud MMRs
        2300             :   real(r8), dimension(vlen), intent(in) :: qcic        ! cloud liquid
        2301             :   real(r8), dimension(vlen), intent(in) :: qiic        ! cloud ice
        2302             :   real(r8), dimension(vlen), intent(in) :: qric        ! rain
        2303             :   real(r8), dimension(vlen), intent(in) :: qsic        ! snow
        2304             :   ! Size parameters
        2305             :   ! rain
        2306             :   real(r8), dimension(vlen), intent(in) :: lamr
        2307             :   real(r8), dimension(vlen), intent(in) :: n0r
        2308             :   ! snow
        2309             :   real(r8), dimension(vlen), intent(in) :: lams
        2310             :   real(r8), dimension(vlen), intent(in) :: n0s
        2311             :   ! Output tendencies
        2312             :   real(r8), dimension(vlen), intent(out) :: pre
        2313             :   real(r8), dimension(vlen), intent(out) :: prds
        2314             :   real(r8), dimension(vlen), intent(out) :: am_evp_st  ! Fractional area where rain evaporates.
        2315             :   ! Switch
        2316             :   logical, intent(in)                    :: evap_rhthrsh
        2317             :
        2318             :   real(r8) :: qclr   ! water vapor mixing ratio in clear air
        2319             :   real(r8) :: ab,abr ! correction to account for latent heat
        2320             :   real(r8) :: eps    ! 1/ sat relaxation timescale
        2321           0 :   real(r8), dimension(vlen) :: dum
        2322             :   integer :: i
        2323             :
        2324             :   !$acc data present (t,rho,dv,mu,sc,q,qvl,qvi,lcldm,precip_frac,arn,asn) &
        2325             :   !$acc      present (qcic,qiic,qric,qsic,lamr,n0r,lams,n0s,pre,prds,am_evp_st) &
        2326             :   !$acc      create  (dum)
        2327             :
        2328             :   ! set temporary cloud fraction to zero if cloud water + ice is very small
        2329             :   ! this will ensure that evaporation/sublimation of precip occurs over
        2330             :   ! entire grid cell, since min cloud fraction is specified otherwise
        2331             :
        2332             :   !$acc parallel vector_length(VLENS) default(present)
        2333             :   !$acc loop gang vector
        2334           0 :   do i=1,vlen
        2335           0 :      am_evp_st(i) = 0._r8
        2336           0 :      if (qcic(i)+qiic(i) < 1.e-6_r8) then
        2337           0 :         dum(i) = 0._r8
        2338             :      else
        2339           0 :         dum(i) = lcldm(i)
        2340             :      end if
        2341             :   end do
        2342             :
        2343             :   !$acc loop gang vector private(qclr,eps,abr,ab)
        2344           0 :   do i=1,vlen
        2345             :      ! only calculate if there is some precip fraction > cloud fraction
        2346           0 :      if (precip_frac(i) > dum(i)) then
        2347           0 :         if (qric(i) >= qsmall .or. qsic(i) >= qsmall) then
        2348           0 :            am_evp_st(i) = precip_frac(i) - dum(i)
        2349             :            ! calculate q for out-of-cloud region
        2350           0 :            qclr=(q(i)-dum(i)*qvl(i))/(1._r8-dum(i))
        2351             :         end if
        2352             :
        2353             :         ! evaporation of rain
        2354           0 :         if (qric(i) >= qsmall) then
        2355           0 :            if (.not.evap_rhthrsh.or.(evap_rhthrsh.and.qclr/qvl(i).le.0.9_r8)) then
        2356           0 :               call calc_ab(t(i), qvl(i), xxlv, abr)
        2357             :               eps = 2._r8*pi*n0r(i)*rho(i)*Dv(i)* &
        2358             :                     (f1r/(lamr(i)*lamr(i))+ &
        2359             :                     f2r*(arn(i)*rho(i)/mu(i))**0.5_r8* &
        2360             :                     sc(i)**(1._r8/3._r8)*gamma_half_br_plus5/ &
        2361           0 :                     (lamr(i)**(5._r8/2._r8+br/2._r8)))
        2362           0 :               pre(i) = eps*(qclr-qvl(i))/abr
        2363             :               ! only evaporate in out-of-cloud region
        2364             :               ! and distribute across precip_frac
        2365           0 :               pre(i)=min(pre(i)*am_evp_st(i),0._r8)
        2366           0 :               pre(i)=pre(i)/precip_frac(i)
        2367             :            else
        2368           0 :               pre(i) = 0._r8
        2369             :            end if
        2370             :         else
        2371           0 :            pre(i) = 0._r8
        2372             :         end if
        2373             :
        2374             :         ! sublimation of snow
        2375           0 :         if (qsic(i) >= qsmall) then
        2376           0 :            call calc_ab(t(i), qvi(i), xxls, ab)
        2377             :            eps = 2._r8*pi*n0s(i)*rho(i)*Dv(i)* &
        2378             :                 (f1s/(lams(i)*lams(i))+ &
        2379             :                 f2s*(asn(i)*rho(i)/mu(i))**0.5_r8* &
        2380             :                 sc(i)**(1._r8/3._r8)*gamma_half_bs_plus5/ &
        2381           0 :                 (lams(i)**(5._r8/2._r8+bs/2._r8)))
        2382           0 :            prds(i) = eps*(qclr-qvi(i))/ab
        2383             :            ! only sublimate in out-of-cloud region and distribute over precip_frac
        2384           0 :            prds(i) = min(prds(i)*am_evp_st(i),0._r8)
        2385           0 :            prds(i) = prds(i)/precip_frac(i)
        2386             :         else
        2387           0 :            prds(i) = 0._r8
        2388             :         end if
        2389             :      else
        2390           0 :         prds(i) = 0._r8
        2391           0 :         pre(i)  = 0._r8
        2392             :      end if
        2393             :   end do
        2394             :   !$acc end parallel
        2395             :
        2396             :   !$acc end data
        2397           0 : end subroutine evaporate_sublimate_precip
        2398             :
        2399             : ! calculate evaporation/sublimation of rain (no snow in mg4)
        2400             : !===================================================================
        2401             : ! note: evaporation/sublimation occurs only in cloud-free portion of grid cell
        2402             : ! in-cloud condensation/deposition of rain and snow is neglected
        2403             : ! except for transfer of cloud water to snow through bergeron process
        2404             :
        2405           0 : subroutine evaporate_sublimate_precip_mg4(t, rho, dv, mu, sc, q, qvl, qvi, &
        2406           0 :      lcldm, precip_frac, arn, qcic, qiic, qric, lamr, n0r, &
        2407           0 :      pre, am_evp_st, vlen, evap_rhthrsh)
        2408             :
        2409             :   integer,  intent(in)                   :: vlen
        2410             :   real(r8), dimension(vlen), intent(in)  :: t           ! temperature
        2411             :   real(r8), dimension(vlen), intent(in)  :: rho         ! air density
        2412             :   real(r8), dimension(vlen), intent(in)  :: dv          ! water vapor diffusivity
        2413             :   real(r8), dimension(vlen), intent(in)  :: mu          ! viscosity
        2414             :   real(r8), dimension(vlen), intent(in)  :: sc          ! schmidt number
        2415             :   real(r8), dimension(vlen), intent(in)  :: q           ! humidity
        2416             :   real(r8), dimension(vlen), intent(in)  :: qvl         ! saturation humidity (water)
        2417             :   real(r8), dimension(vlen), intent(in)  :: qvi         ! saturation humidity (ice)
        2418             :   real(r8), dimension(vlen), intent(in)  :: lcldm       ! liquid cloud fraction
        2419             :   real(r8), dimension(vlen), intent(in)  :: precip_frac ! precipitation fraction (maximum overlap)
        2420             :   ! fallspeed parameters
        2421             :   real(r8), dimension(vlen), intent(in)  :: arn         ! rain
        2422             :   ! In-cloud MMRs
        2423             :   real(r8), dimension(vlen), intent(in)  :: qcic        ! cloud liquid
        2424             :   real(r8), dimension(vlen), intent(in)  :: qiic        ! cloud ice
        2425             :   real(r8), dimension(vlen), intent(in)  :: qric        ! rain
        2426             :   ! Size parameters
        2427             :   ! rain
        2428             :   real(r8), dimension(vlen), intent(in)  :: lamr
        2429             :   real(r8), dimension(vlen), intent(in)  :: n0r
        2430             :   ! Output tendencies
        2431             :   real(r8), dimension(vlen), intent(out) :: pre
        2432             :   real(r8), dimension(vlen), intent(out) :: am_evp_st   ! Fractional area where rain evaporates.
        2433             :   ! Switch
        2434             :   logical, intent(in)                    :: evap_rhthrsh
        2435             :   real(r8) :: qclr   ! water vapor mixing ratio in clear air
        2436             :   real(r8) :: abr    ! correction to account for latent heat
        2437             :   real(r8) :: eps    ! 1/ sat relaxation timescale
        2438           0 :   real(r8), dimension(vlen) :: dum
        2439             :   integer :: i
        2440             :
        2441             :   !$acc data present (t,rho,dv,mu,sc,q,qvl,qvi,lcldm,precip_frac,arn) &
        2442             :   !$acc      present (qcic,qiic,qric,lamr,n0r,pre,am_evp_st) &
        2443             :   !$acc      create  (dum)
        2444             :
        2445             :   ! set temporary cloud fraction to zero if cloud water + ice is very small
        2446             :   ! this will ensure that evaporation/sublimation of precip occurs over
        2447             :   ! entire grid cell, since min cloud fraction is specified otherwise
        2448             :
        2449             :   !$acc parallel vector_length(VLENS) default(present)
        2450             :   !$acc loop gang vector
        2451           0 :   do i=1,vlen
        2452           0 :      am_evp_st(i) = 0._r8
        2453           0 :      if (qcic(i)+qiic(i) < 1.e-6_r8) then
        2454           0 :         dum(i) = 0._r8
        2455             :      else
        2456           0 :         dum(i) = lcldm(i)
        2457             :      end if
        2458             :   end do
        2459             :
        2460             :   !$acc loop gang vector private(qclr,eps,abr)
        2461           0 :   do i=1,vlen
        2462             :      ! only calculate if there is some precip fraction > cloud fraction
        2463           0 :      if (precip_frac(i) > dum(i)) then
        2464           0 :         if (qric(i) >= qsmall) then
        2465           0 :            am_evp_st(i) = precip_frac(i) - dum(i)
        2466             :            ! calculate q for out-of-cloud region
        2467           0 :            qclr=(q(i)-dum(i)*qvl(i))/(1._r8-dum(i))
        2468             :         end if
        2469             :
        2470             :         ! evaporation of rain
        2471           0 :         if (qric(i) >= qsmall) then
        2472           0 :            if (.not.evap_rhthrsh.or.(evap_rhthrsh.and.qclr/qvl(i).le.0.9_r8)) then
        2473           0 :               call calc_ab(t(i), qvl(i), xxlv, abr)
        2474             :               eps = 2._r8*pi*n0r(i)*rho(i)*Dv(i)* &
        2475             :                  (f1r/(lamr(i)*lamr(i))+ &
        2476             :                  f2r*(arn(i)*rho(i)/mu(i))**0.5_r8* &
        2477             :                  sc(i)**(1._r8/3._r8)*gamma_half_br_plus5/ &
        2478           0 :                  (lamr(i)**(5._r8/2._r8+br/2._r8)))
        2479           0 :               pre(i) = eps*(qclr-qvl(i))/abr
        2480             :               ! only evaporate in out-of-cloud region
        2481             :               ! and distribute across precip_frac
        2482           0 :               pre(i) = min(pre(i)*am_evp_st(i),0._r8)
        2483           0 :               pre(i) = pre(i)/precip_frac(i)
        2484             :            else
        2485           0 :               pre(i) = 0._r8
        2486             :            end if
        2487             :         else
        2488           0 :            pre(i) = 0._r8
        2489             :         end if
        2490             :      else
        2491           0 :         pre(i) = 0._r8
        2492             :      end if
        2493             :   end do
        2494             :   !$acc end parallel
        2495             :
        2496             :   !$acc end data
        2497           0 : end subroutine evaporate_sublimate_precip_mg4
        2498             :
        2499             : ! evaporation/sublimation of rain, snow and graupel
        2500             : !===================================================================
        2501             : ! note: evaporation/sublimation occurs only in cloud-free portion of grid cell
        2502             : ! in-cloud condensation/deposition of rain and snow is neglected
        2503             : ! except for transfer of cloud water to snow through bergeron process
        2504             :
        2505           0 : subroutine evaporate_sublimate_precip_graupel(t, rho, dv, mu, sc, q, qvl, qvi, &
        2506           0 :      lcldm, precip_frac, arn, asn, agn, bg, qcic, qiic, qric, qsic, qgic, lamr, n0r, lams, n0s, lamg, n0g, &
        2507           0 :      pre, prds, prdg, am_evp_st, vlen, evap_rhthrsh)
        2508             :
        2509             :   integer,                   intent(in)  :: vlen
        2510             :   real(r8), dimension(vlen), intent(in)  :: t           ! temperature
        2511             :   real(r8), dimension(vlen), intent(in)  :: rho         ! air density
        2512             :   real(r8), dimension(vlen), intent(in)  :: dv          ! water vapor diffusivity
        2513             :   real(r8), dimension(vlen), intent(in)  :: mu          ! viscosity
        2514             :   real(r8), dimension(vlen), intent(in)  :: sc          ! schmidt number
        2515             :   real(r8), dimension(vlen), intent(in)  :: q           ! humidity
        2516             :   real(r8), dimension(vlen), intent(in)  :: qvl         ! saturation humidity (water)
        2517             :   real(r8), dimension(vlen), intent(in)  :: qvi         ! saturation humidity (ice)
        2518             :   real(r8), dimension(vlen), intent(in)  :: lcldm       ! liquid cloud fraction
        2519             :   real(r8), dimension(vlen), intent(in)  :: precip_frac ! precipitation fraction (maximum overlap)
        2520             :   ! fallspeed parameters
        2521             :   real(r8), dimension(vlen), intent(in)  :: arn         ! rain
        2522             :   real(r8), dimension(vlen), intent(in)  :: asn         ! snow
        2523             :   real(r8), dimension(vlen), intent(in)  :: agn         ! graupel
        2524             :   real(r8),                  intent(in)  :: bg
        2525             :   ! In-cloud MMRs
        2526             :   real(r8), dimension(vlen), intent(in)  :: qcic        ! cloud liquid
        2527             :   real(r8), dimension(vlen), intent(in)  :: qiic        ! cloud ice
        2528             :   real(r8), dimension(vlen), intent(in)  :: qric        ! rain
        2529             :   real(r8), dimension(vlen), intent(in)  :: qsic        ! snow
        2530             :   real(r8), dimension(vlen), intent(in)  :: qgic        ! graupel
        2531             :   ! Size parameters
        2532             :   ! rain
        2533             :   real(r8), dimension(vlen), intent(in)  :: lamr
        2534             :   real(r8), dimension(vlen), intent(in)  :: n0r
        2535             :   ! snow
        2536             :   real(r8), dimension(vlen), intent(in)  :: lams
        2537             :   real(r8), dimension(vlen), intent(in)  :: n0s
        2538             :   ! graupel
        2539             :   real(r8), dimension(vlen), intent(in)  :: lamg
        2540             :   real(r8), dimension(vlen), intent(in)  :: n0g
        2541             :   ! Output tendencies
        2542             :   real(r8), dimension(vlen), intent(out) :: pre
        2543             :   real(r8), dimension(vlen), intent(out) :: prds
        2544             :   real(r8), dimension(vlen), intent(out) :: prdg
        2545             :   real(r8), dimension(vlen), intent(out) :: am_evp_st   ! Fractional area where rain evaporates.
        2546             :   ! Switch
        2547             :   logical, intent(in) :: evap_rhthrsh
        2548             :
        2549             :   real(r8) :: qclr   ! water vapor mixing ratio in clear air
        2550             :   real(r8) :: ab, abr, abg    ! correction to account for latent heat
        2551             :   real(r8) :: eps    ! 1/ sat relaxation timescale
        2552           0 :   real(r8), dimension(vlen) :: dum
        2553             :   integer :: i
        2554             :
        2555             :   !$acc data present (t,rho,dv,mu,sc,q,qvl,qvi,lcldm,precip_frac) &
        2556             :   !$acc      present (arn,asn,agn,qcic,qiic,qric,qsic,qgic,lamr)  &
        2557             :   !$acc      present (n0r,lams,n0s,lamg,n0g,pre,prds,prdg,am_evp_st) &
        2558             :   !$acc      create  (dum)
        2559             :
        2560             :   ! set temporary cloud fraction to zero if cloud water + ice is very small
        2561             :   ! this will ensure that evaporation/sublimation of precip occurs over
        2562             :   ! entire grid cell, since min cloud fraction is specified otherwise
        2563             :
        2564             :   !$acc parallel vector_length(VLENS) default(present)
        2565             :   !$acc loop gang vector
        2566           0 :   do i=1,vlen
        2567           0 :      am_evp_st(i) = 0._r8
        2568           0 :      if (qcic(i)+qiic(i) < 1.e-6_r8) then
        2569           0 :         dum(i) = 0._r8
        2570             :      else
        2571           0 :         dum(i) = lcldm(i)
        2572             :      end if
        2573             :   end do
        2574             :
        2575             :   !$acc loop gang vector private(qclr,eps,abr,ab,abg)
        2576           0 :   do i=1,vlen
        2577             :      ! only calculate if there is some precip fraction > cloud fraction
        2578           0 :      if (precip_frac(i) > dum(i)) then
        2579           0 :         if (qric(i) >= qsmall .or. qsic(i) >= qsmall .or. qgic(i) >= qsmall) then
        2580           0 :            am_evp_st(i) = precip_frac(i) - dum(i)
        2581             :            ! calculate q for out-of-cloud region
        2582           0 :            qclr=(q(i)-dum(i)*qvl(i))/(1._r8-dum(i))
        2583             :         end if
        2584             :
        2585             :         ! evaporation of rain
        2586           0 :         if (qric(i) >= qsmall) then
        2587           0 :            if (.not.evap_rhthrsh.or.(evap_rhthrsh.and.qclr/qvl(i).le.0.9_r8)) then
        2588           0 :               call calc_ab(t(i), qvl(i), xxlv, abr)
        2589             :               eps = 2._r8*pi*n0r(i)*rho(i)*Dv(i)* &
        2590             :                  (f1r/(lamr(i)*lamr(i))+ &
        2591             :                  f2r*(arn(i)*rho(i)/mu(i))**0.5_r8* &
        2592             :                  sc(i)**(1._r8/3._r8)*gamma_half_br_plus5/ &
        2593           0 :                  (lamr(i)**(5._r8/2._r8+br/2._r8)))
        2594           0 :               pre(i) = eps*(qclr-qvl(i))/abr
        2595             :               ! only evaporate in out-of-cloud region
        2596             :               ! and distribute across precip_frac
        2597           0 :               pre(i) = min(pre(i)*am_evp_st(i),0._r8)
        2598           0 :               pre(i) = pre(i)/precip_frac(i)
        2599             :            else
        2600           0 :               pre(i) = 0._r8
        2601             :            end if
        2602             :         else
        2603           0 :            pre(i) = 0._r8
        2604             :         end if
        2605             :
        2606             :         ! sublimation of snow
        2607           0 :         if (qsic(i) >= qsmall) then
        2608           0 :            call calc_ab(t(i), qvi(i), xxls, ab)
        2609             :            eps = 2._r8*pi*n0s(i)*rho(i)*Dv(i)* &
        2610             :                 (f1s/(lams(i)*lams(i))+ &
        2611             :                 f2s*(asn(i)*rho(i)/mu(i))**0.5_r8* &
        2612             :                 sc(i)**(1._r8/3._r8)*gamma_half_bs_plus5/ &
        2613           0 :                 (lams(i)**(5._r8/2._r8+bs/2._r8)))
        2614           0 :            prds(i) = eps*(qclr-qvi(i))/ab
        2615             :            ! only sublimate in out-of-cloud region and distribute over precip_frac
        2616           0 :            prds(i) = min(prds(i)*am_evp_st(i),0._r8)
        2617           0 :            prds(i) = prds(i)/precip_frac(i)
        2618             :         else
        2619           0 :            prds(i) = 0._r8
        2620             :         end if
        2621             :
        2622             :         ! ADD GRAUPEL, do Same with prdg.
        2623           0 :         if (qgic(i).ge.qsmall) then
        2624           0 :            call calc_ab(t(i), qvi(i), xxls, abg)
        2625             :            eps = 2._r8*pi*n0g(i)*rho(i)*Dv(i)*                    &
        2626             :                 (f1s/(lamg(i)*lamg(i))+                           &
        2627             :                 f2s*(agn(i)*rho(i)/mu(i))**0.5_r8*                &
        2628             :                 sc(i)**(1._r8/3._r8)*gamma(5._r8/2._r8+bg/2._r8)/ &
        2629           0 :                 (lamg(i)**(5._r8/2._r8+bs/2._r8)))
        2630           0 :            prdg(i) = eps*(qclr-qvi(i))/abg
        2631             :            ! only sublimate in out-of-cloud region and distribute over precip_frac
        2632           0 :            prdg(i) = min(prdg(i)*am_evp_st(i),0._r8)
        2633           0 :            prdg(i) = prdg(i)/precip_frac(i)
        2634             :         else
        2635           0 :            prdg(i) = 0._r8
        2636             :         end if
        2637             :
        2638             :      else
        2639           0 :         prds(i) = 0._r8
        2640           0 :         pre(i)  = 0._r8
        2641           0 :         prdg(i) = 0._r8
        2642             :      end if
        2643             :   end do
        2644             :   !$acc end parallel
        2645             :
        2646             :   !$acc end data
        2647           0 : end subroutine evaporate_sublimate_precip_graupel
        2648             :
        2649             : ! evaporation/sublimation of rain and graupel (no snow in mg4)
        2650             : !===================================================================
        2651             : ! note: evaporation/sublimation occurs only in cloud-free portion of grid cell
        2652             : ! in-cloud condensation/deposition of rain and graupel is neglected
        2653             : ! except for transfer of cloud water to snow through bergeron process
        2654             :
        2655           0 : subroutine evaporate_sublimate_precip_graupel_mg4(t, rho, dv, mu, sc, q, qvl, qvi, &
        2656           0 :      lcldm, precip_frac, arn, agn, bg, qcic, qiic, qric, qgic, lamr, n0r, lamg, n0g, &
        2657           0 :      pre, prdg, am_evp_st, vlen)
        2658             :
        2659             :   integer,                   intent(in) :: vlen
        2660             :   real(r8), dimension(vlen), intent(in) :: t            ! temperature
        2661             :   real(r8), dimension(vlen), intent(in) :: rho          ! air density
        2662             :   real(r8), dimension(vlen), intent(in) :: dv           ! water vapor diffusivity
        2663             :   real(r8), dimension(vlen), intent(in) :: mu           ! viscosity
        2664             :   real(r8), dimension(vlen), intent(in) :: sc           ! schmidt number
        2665             :   real(r8), dimension(vlen), intent(in) :: q            ! humidity
        2666             :   real(r8), dimension(vlen), intent(in) :: qvl          ! saturation humidity (water)
        2667             :   real(r8), dimension(vlen), intent(in) :: qvi          ! saturation humidity (ice)
        2668             :   real(r8), dimension(vlen), intent(in) :: lcldm        ! liquid cloud fraction
        2669             :   real(r8), dimension(vlen), intent(in) :: precip_frac  ! precipitation fraction (maximum overlap)
        2670             :   ! fallspeed parameters
        2671             :   real(r8), dimension(vlen), intent(in) :: arn          ! rain
        2672             :   real(r8), dimension(vlen), intent(in) :: agn          ! graupel
        2673             :   real(r8),                  intent(in) :: bg
        2674             :   ! In-cloud MMRs
        2675             :   real(r8), dimension(vlen), intent(in) :: qcic         ! cloud liquid
        2676             :   real(r8), dimension(vlen), intent(in) :: qiic         ! cloud ice
        2677             :   real(r8), dimension(vlen), intent(in) :: qric         ! rain
        2678             :   real(r8), dimension(vlen), intent(in) :: qgic         ! graupel
        2679             :   ! Size parameters
        2680             :   ! rain
        2681             :   real(r8), dimension(vlen), intent(in) :: lamr
        2682             :   real(r8), dimension(vlen), intent(in) :: n0r
        2683             :   ! graupel
        2684             :   real(r8), dimension(vlen), intent(in) :: lamg
        2685             :   real(r8), dimension(vlen), intent(in) :: n0g
        2686             :   ! Output tendencies
        2687             :   real(r8), dimension(vlen), intent(out) :: pre
        2688             :   real(r8), dimension(vlen), intent(out) :: prdg
        2689             :   real(r8), dimension(vlen), intent(out) :: am_evp_st   ! Fractional area where rain evaporates.
        2690             :
        2691             :   real(r8) :: qclr                                      ! water vapor mixing ratio in clear air
        2692             :   real(r8) :: abr, abg                                  ! correction to account for latent heat
        2693             :   real(r8) :: eps                                       ! 1/ sat relaxation timescale
        2694           0 :   real(r8), dimension(vlen) :: dum
        2695             :   integer :: i
        2696             :
        2697             :   !$acc data present (t,rho,dv,mu,sc,q,qvl,qvi,lcldm,precip_frac) &
        2698             :   !$acc      present (arn,agn,qcic,qiic,qric,qgic,lamr,n0r,lamg)  &
        2699             :   !$acc      present (n0g,pre,prdg,am_evp_st) &
        2700             :   !$acc      create  (dum)
        2701             :
        2702             :   ! set temporary cloud fraction to zero if cloud water + ice is very small
        2703             :   ! this will ensure that evaporation/sublimation of precip occurs over
        2704             :   ! entire grid cell, since min cloud fraction is specified otherwise
        2705             :   !$acc parallel vector_length(VLENS) default(present)
        2706             :   !$acc loop gang vector
        2707           0 :   do i=1,vlen
        2708           0 :      am_evp_st(i) = 0._r8
        2709           0 :      if (qcic(i)+qiic(i) < 1.e-6_r8) then
        2710           0 :         dum(i) = 0._r8
        2711             :      else
        2712           0 :         dum(i) = lcldm(i)
        2713             :      end if
        2714             :   end do
        2715             :
        2716             :   !$acc loop gang vector private(qclr,eps,abr,abg)
        2717           0 :   do i=1,vlen
        2718             :      ! only calculate if there is some precip fraction > cloud fraction
        2719           0 :      if (precip_frac(i) > dum(i)) then
        2720           0 :         if (qric(i) >= qsmall .or. qgic(i) >= qsmall) then
        2721           0 :            am_evp_st(i) = precip_frac(i) - dum(i)
        2722             :            ! calculate q for out-of-cloud region
        2723           0 :            qclr = (q(i)-dum(i)*qvl(i))/(1._r8-dum(i))
        2724             :         end if
        2725             :
        2726             :         ! evaporation of rain
        2727           0 :         if (qric(i) >= qsmall) then
        2728           0 :            call calc_ab(t(i), qvl(i), xxlv, abr)
        2729             :            eps = 2._r8*pi*n0r(i)*rho(i)*Dv(i)* &
        2730             :                 (f1r/(lamr(i)*lamr(i))+ &
        2731             :                 f2r*(arn(i)*rho(i)/mu(i))**0.5_r8* &
        2732             :                 sc(i)**(1._r8/3._r8)*gamma_half_br_plus5/ &
        2733           0 :                 (lamr(i)**(5._r8/2._r8+br/2._r8)))
        2734           0 :            pre(i) = eps*(qclr-qvl(i))/abr
        2735             :            ! only evaporate in out-of-cloud region
        2736             :            ! and distribute across precip_frac
        2737           0 :            pre(i) = min(pre(i)*am_evp_st(i),0._r8)
        2738           0 :            pre(i) = pre(i)/precip_frac(i)
        2739             :         else
        2740           0 :            pre(i) = 0._r8
        2741             :         end if
        2742             :
        2743             :         ! ADD GRAUPEL, do Same with prdg.
        2744           0 :         if (qgic(i).ge.qsmall) then
        2745           0 :            call calc_ab(t(i), qvi(i), xxls, abg)
        2746             :            eps = 2._r8*pi*n0g(i)*rho(i)*Dv(i)*                    &
        2747             :                 (f1s/(lamg(i)*lamg(i))+                           &
        2748             :                 f2s*(agn(i)*rho(i)/mu(i))**0.5_r8*                &
        2749             :                 sc(i)**(1._r8/3._r8)*gamma(5._r8/2._r8+bg/2._r8)/ &
        2750           0 :                 (lamg(i)**(5._r8/2._r8+bs/2._r8)))
        2751           0 :            prdg(i) = eps*(qclr-qvi(i))/abg
        2752             :            ! only sublimate in out-of-cloud region and distribute over precip_frac
        2753           0 :            prdg(i) = min(prdg(i)*am_evp_st(i),0._r8)
        2754           0 :            prdg(i) = prdg(i)/precip_frac(i)
        2755             :         else
        2756           0 :            prdg(i) = 0._r8
        2757             :         end if
        2758             :      else
        2759           0 :         pre(i) = 0._r8
        2760           0 :         prdg(i) = 0._r8
        2761             :      end if
        2762             :   end do
        2763             :   !$acc end parallel
        2764             :
        2765             :   !$acc end data
        2766           0 : end subroutine evaporate_sublimate_precip_graupel_mg4
        2767             :
        2768             : ! bergeron process - evaporation of droplets and deposition onto snow
        2769             : !===================================================================
        2770             :
        2771           0 : subroutine bergeron_process_snow(t, rho, dv, mu, sc, qvl, qvi, asn, &
        2772           0 :      qcic, qsic, lams, n0s, bergs, vlen)
        2773             :
        2774             :   integer,                   intent(in) :: vlen
        2775             :   real(r8), dimension(vlen), intent(in) :: t    ! temperature
        2776             :   real(r8), dimension(vlen), intent(in) :: rho  ! air density
        2777             :   real(r8), dimension(vlen), intent(in) :: dv   ! water vapor diffusivity
        2778             :   real(r8), dimension(vlen), intent(in) :: mu   ! viscosity
        2779             :   real(r8), dimension(vlen), intent(in) :: sc   ! schmidt number
        2780             :   real(r8), dimension(vlen), intent(in) :: qvl  ! saturation humidity (water)
        2781             :   real(r8), dimension(vlen), intent(in) :: qvi  ! saturation humidity (ice)
        2782             :   ! fallspeed parameter for snow
        2783             :   real(r8), dimension(vlen), intent(in) :: asn
        2784             :   ! In-cloud MMRs
        2785             :   real(r8), dimension(vlen), intent(in) :: qcic ! cloud liquid mixing ratio
        2786             :   real(r8), dimension(vlen), intent(in) :: qsic ! snow mixing ratio
        2787             :   ! Size parameters for snow
        2788             :   real(r8), dimension(vlen), intent(in) :: lams
        2789             :   real(r8), dimension(vlen), intent(in) :: n0s
        2790             :   ! Output tendencies
        2791             :   real(r8), dimension(vlen), intent(out) :: bergs
        2792             :
        2793             :   real(r8) :: ab     ! correction to account for latent heat
        2794             :   real(r8) :: eps    ! 1/ sat relaxation timescale
        2795             :   integer :: i
        2796             :
        2797             :   !$acc data present (t,rho,dv,mu,sc,qvl,qvi) &
        2798             :   !$acc      present (asn,qcic,qsic,lams,n0s,bergs)
        2799             :
        2800             :   !$acc parallel vector_length(VLENS) default(present)
        2801             :   !$acc loop gang vector private(eps,ab)
        2802           0 :   do i=1,vlen
        2803           0 :      if (qsic(i) >= qsmall.and. qcic(i) >= qsmall .and. t(i) < tmelt) then
        2804           0 :         call calc_ab(t(i), qvi(i), xxls, ab)
        2805             :         eps = 2._r8*pi*n0s(i)*rho(i)*Dv(i)* &
        2806             :              (f1s/(lams(i)*lams(i))+ &
        2807             :              f2s*(asn(i)*rho(i)/mu(i))**0.5_r8* &
        2808             :              sc(i)**(1._r8/3._r8)*gamma_half_bs_plus5/ &
        2809           0 :              (lams(i)**(5._r8/2._r8+bs/2._r8)))
        2810           0 :         bergs(i) = eps*(qvl(i)-qvi(i))/ab
        2811             :      else
        2812           0 :         bergs(i) = 0._r8
        2813             :      end if
        2814             :   end do
        2815             :   !$acc end parallel
        2816             :
        2817             :   !$acc end data
        2818           0 : end subroutine bergeron_process_snow
        2819             :
        2820             : !========================================================================
        2821             : ! Collection of snow by rain to form graupel
        2822             : !========================================================================
        2823             :
        2824           0 : subroutine graupel_collecting_snow(qsic,qric,umr,ums,rho,lamr,n0r,lams,n0s, &
        2825           0 :      psacr, vlen)
        2826             :
        2827             :   integer,                   intent(in)  :: vlen
        2828             :   ! In-cloud MMRs
        2829             :   real(r8), dimension(vlen), intent(in)  :: qsic  ! snow
        2830             :   real(r8), dimension(vlen), intent(in)  :: qric  ! rain
        2831             :   ! mass-weighted fall speeds
        2832             :   real(r8), dimension(vlen), intent(in)  :: umr   ! rain
        2833             :   real(r8), dimension(vlen), intent(in)  :: ums   ! snow
        2834             :   real(r8), dimension(vlen), intent(in)  :: rho   ! air density
        2835             :   ! Size parameters for rain
        2836             :   real(r8), dimension(vlen), intent(in)  :: lamr
        2837             :   real(r8), dimension(vlen), intent(in)  :: n0r
        2838             :   ! Size parameters for snow
        2839             :   real(r8), dimension(vlen), intent(in)  :: lams
        2840             :   real(r8), dimension(vlen), intent(in)  :: n0s
        2841             :   real(r8), dimension(vlen), intent(out) :: psacr ! conversion due to coll of snow by rain
        2842             :
        2843             :   real(r8), parameter :: cons31 = pi*pi*ecr*rhosn
        2844             :   integer :: i
        2845             :
        2846             :   !$acc data present (qsic,qric,umr,ums,rho,lamr) &
        2847             :   !$acc      present (n0r,lams,n0s,psacr)
        2848             :
        2849             :   !$acc parallel vector_length(VLENS) default(present)
        2850             :   !$acc loop gang vector
        2851           0 :   do i=1,vlen
        2852           0 :      if (qsic(i).ge.0.1e-3_r8 .and. qric(i).ge.0.1e-3_r8) then
        2853             :         psacr(i) = cons31*(((1.2_r8*umr(i)-0.95_r8*ums(i))**2+ &
        2854             :              0.08_r8*ums(i)*umr(i))**0.5_r8*rho(i)*            &
        2855             :              n0r(i)*n0s(i)/lams(i)**3*                         &
        2856             :              (5._r8/(lams(i)**3*lamr(i))+                      &
        2857             :              2._r8/(lams(i)**2*lamr(i)**2)+                    &
        2858           0 :              0.5_r8/(lams(i)*lamr(i)**3)))
        2859             :      else
        2860           0 :         psacr(i) = 0._r8
        2861             :      end if
        2862             :   end do
        2863             :   !$acc end parallel
        2864             :
        2865             :   !$acc end data
        2866           0 : end subroutine graupel_collecting_snow
        2867             :
        2868             : !========================================================================
        2869             : ! Collection of cloud water by graupel
        2870             : !========================================================================
        2871             :
        2872           0 : subroutine graupel_collecting_cld_water(qgic,qcic,ncic,rho,n0g,lamg,bg,agn, &
        2873           0 :      psacwg, npsacwg, vlen)
        2874             :
        2875             :   integer,                   intent(in) :: vlen
        2876             :   ! In-cloud MMRs
        2877             :   real(r8), dimension(vlen), intent(in) :: qgic ! graupel
        2878             :   real(r8), dimension(vlen), intent(in) :: qcic ! cloud water
        2879             :   real(r8), dimension(vlen), intent(in) :: ncic ! cloud water number conc
        2880             :   real(r8), dimension(vlen), intent(in) :: rho  ! air density
        2881             :   ! Size parameters for graupel
        2882             :   real(r8), dimension(vlen), intent(in) :: lamg
        2883             :   real(r8), dimension(vlen), intent(in) :: n0g
        2884             :   ! fallspeed parameters for graupel
        2885             :   ! Set AGN and BG  as input (in micro_pumas_v1.F90)  (select hail or graupel as appropriate)
        2886             :   real(r8),                  intent(in) :: bg
        2887             :   real(r8), dimension(vlen), intent(in) :: agn
        2888             :   ! Output tendencies
        2889             :   real(r8), dimension(vlen), intent(out) :: psacwg
        2890             :   real(r8), dimension(vlen), intent(out) :: npsacwg
        2891             :
        2892             :   real(r8) :: cons
        2893             :   integer  :: i
        2894             :
        2895           0 :   cons = gamma(bg + 3._r8)*pi/4._r8 * ecid
        2896             :
        2897             :   !$acc data present (qgic,qcic,ncic,rho,lamg,n0g,agn,psacwg,npsacwg)
        2898             :
        2899             :   !$acc parallel vector_length(VLENS) default(present)
        2900             :   !$acc loop gang vector
        2901           0 :   do i=1,vlen
        2902           0 :      if (qgic(i).ge.1.e-8_r8 .and. qcic(i).ge.qsmall) then
        2903             :         psacwg(i) = cons*agn(i)*qcic(i)*rho(i)*  &
        2904             :                n0g(i)/                           &
        2905           0 :                lamg(i)**(bg+3._r8)
        2906             :         npsacwg(i) = cons*agn(i)*ncic(i)*rho(i)* &
        2907             :                n0g(i)/                           &
        2908           0 :                lamg(i)**(bg+3._r8)
        2909             :      else
        2910           0 :         psacwg(i)  = 0._r8
        2911           0 :         npsacwg(i) = 0._r8
        2912             :      end if
        2913             :   end do
        2914             :   !$acc end parallel
        2915             :
        2916             :   !$acc end data
        2917           0 : end subroutine graupel_collecting_cld_water
        2918             :
        2919             : !========================================================================
        2920             : ! Conversion of rimed cloud water onto snow to graupel/hail
        2921             : !========================================================================
        2922             :
        2923           0 : subroutine graupel_riming_liquid_snow(psacws,qsic,qcic,nsic,rho,rhosn,rhog,asn,lams,n0s,dtime, &
        2924           0 :      pgsacw,nscng,vlen)
        2925             :
        2926             :   integer,                   intent(in)    :: vlen
        2927             :   ! Accretion of cloud water to snow tendency (modified)
        2928             :   real(r8), dimension(vlen), intent(inout) :: psacws
        2929             :   real(r8), dimension(vlen), intent(in)    :: qsic    ! snow mixing ratio
        2930             :   real(r8), dimension(vlen), intent(in)    :: qcic    ! cloud liquid mixing ratio
        2931             :   real(r8), dimension(vlen), intent(in)    :: nsic    ! snow number concentration
        2932             :   real(r8), dimension(vlen), intent(in)    :: rho     ! air density
        2933             :   real(r8),                  intent(in)    :: rhosn   ! snow density
        2934             :   real(r8),                  intent(in)    :: rhog    ! graupel density
        2935             :   real(r8), dimension(vlen), intent(in)    :: asn     ! fall speed parameter for snow
        2936             :   ! Size parameters for snow
        2937             :   real(r8), dimension(vlen), intent(in)    :: lams
        2938             :   real(r8), dimension(vlen), intent(in)    :: n0s
        2939             :   real(r8),                  intent(in)    :: dtime
        2940             :   !Output process rates
        2941             :   real(r8), dimension(vlen), intent(out)   :: pgsacw  ! dQ graupel due to collection droplets by snow
        2942             :   real(r8), dimension(vlen), intent(out)   :: nscng   ! dN graupel due to collection droplets by snow
        2943             :
        2944             :   real(r8) :: cons
        2945             :   real(r8) :: rhosu
        2946             :   real(r8) :: dum
        2947             :   integer  :: i
        2948             :
        2949             : !........................................................................
        2950             : !Input: PSACWS,qs,qc,n0s,rho,lams,rhos,rhog
        2951             : !Output:PSACWS,PGSACW,NSCNG
        2952             :
        2953           0 :   rhosu = 85000._r8/(ra * tmelt)    ! typical air density at 850 mb
        2954             :
        2955           0 :   cons  = 4._r8 *2._r8 *3._r8*rhosu*pi*ecid*ecid*gamma_2bs_plus2/(8._r8*(rhog-rhosn))
        2956             :
        2957             :   !$acc data present (psacws,qsic,qcic,nsic,rho,asn,lams,n0s,pgsacw,nscng)
        2958             :
        2959             :   !$acc parallel vector_length(VLENS) default(present)
        2960             :   !$acc loop gang vector private(dum)
        2961           0 :   do i=1,vlen
        2962           0 :      if (psacws(i).gt.0._r8 .and. qsic(i).GE.0.1e-3_r8 .AND. qcic(i).GE.0.5E-3_r8) then
        2963             :         ! Only allow conversion if qni > 0.1 and qc > 0.5 g/kg following Rutledge and Hobbs (1984)
        2964             :         ! if (qsic(i).GE.0.1e-3_r8 .AND. qcic(i).GE.0.5E-3_r8) then
        2965             :
        2966             :         ! portion of riming converted to graupel (Reisner et al. 1998, originally IS1991)
        2967             :         ! dtime here is correct.
        2968             :         pgsacw(i) = min(psacws(i),cons*dtime*n0s(i)*qcic(i)*qcic(i)* &
        2969             :              asn(i)*asn(i)/ &
        2970           0 :              (rho(i)*lams(i)**(2._r8*bs+2._r8)))
        2971             :
        2972             :         ! Mix rat converted into graupel as embryo (Reisner et al. 1998, orig M1990)
        2973           0 :         dum       = max(rhosn/(rhog-rhosn)*pgsacw(i),0._r8)
        2974             :
        2975             :         ! Number concentraiton of embryo graupel from riming of snow
        2976           0 :         nscng(i)  = dum/mg0*rho(i)
        2977             :         ! Limit max number converted to snow number  (dtime here correct? We think yes)
        2978           0 :         nscng(i)  = min(nscng(i),nsic(i)/dtime)
        2979             :
        2980             :         ! Portion of riming left for snow
        2981           0 :         psacws(i) = psacws(i) - pgsacw(i)
        2982             :      else
        2983           0 :         pgsacw(i) = 0._r8
        2984           0 :         nscng(i)  = 0._r8
        2985             :      end if
        2986             :   end do
        2987             :   !$acc end parallel
        2988             :
        2989             :   !$acc end data
        2990           0 : end subroutine graupel_riming_liquid_snow
        2991             :
        2992             : !========================================================================
        2993             : !CHANGE IN Q,N COLLECTION RAIN BY GRAUPEL
        2994             : !========================================================================
        2995             :
        2996           0 : subroutine graupel_collecting_rain(qric,qgic,umg,umr,ung,unr,rho,n0r,lamr,n0g,lamg,&
        2997           0 :      pracg,npracg,vlen)
        2998             :
        2999             :   integer,                   intent(in) :: vlen
        3000             :   !MMR
        3001             :   real(r8), dimension(vlen), intent(in) :: qric  !rain mixing ratio
        3002             :   real(r8), dimension(vlen), intent(in) :: qgic  !graupel mixing ratio
        3003             :   !Mass weighted Fall speeds
        3004             :   real(r8), dimension(vlen), intent(in) :: umg ! graupel fall speed
        3005             :   real(r8), dimension(vlen), intent(in) :: umr ! rain fall speed
        3006             :   !Number weighted fall speeds
        3007             :   real(r8), dimension(vlen), intent(in) :: ung ! graupel fall speed
        3008             :   real(r8), dimension(vlen), intent(in) :: unr ! rain fall speed
        3009             :   real(r8), dimension(vlen), intent(in) :: rho   ! air density
        3010             :   ! Size parameters for rain
        3011             :   real(r8), dimension(vlen), intent(in) :: n0r
        3012             :   real(r8), dimension(vlen), intent(in) :: lamr
        3013             :   ! Size parameters for graupel
        3014             :   real(r8), dimension(vlen), intent(in) :: n0g
        3015             :   real(r8), dimension(vlen), intent(in) :: lamg
        3016             :   !Output process rates
        3017             :   real(r8), dimension(vlen), intent(out) :: pracg   ! Q collection rain by graupel
        3018             :   real(r8), dimension(vlen), intent(out) :: npracg  ! N collection rain by graupel
        3019             :
        3020             : ! Add collection of graupel by rain above freezing
        3021             : ! assume all rain collection by graupel above freezing is shed
        3022             : ! assume shed drops are 1 mm in size
        3023             :   integer  :: i
        3024             :   real(r8), parameter :: cons41 = pi*pi*ecr*rhow
        3025             :   real(r8), parameter :: cons32 = pi/2._r8*ecr
        3026             :   real(r8) :: dum
        3027             :
        3028             :   !$acc data present (qric,qgic,umg,umr,ung,unr,rho) &
        3029             :   !$acc      present (n0r,lamr,n0g,lamg,pracg,npracg)
        3030             :
        3031             :   !$acc parallel vector_length(VLENS) default(present)
        3032             :   !$acc loop gang vector private(dum)
        3033           0 :   do i=1,vlen
        3034           0 :      if (qric(i).ge.1.e-8_r8.and.qgic(i).ge.1.e-8_r8) then
        3035             :         ! pracg is mixing ratio of rain per sec collected by graupel/hail
        3036             :         pracg(i) = cons41*(((1.2_r8*umr(i)-0.95_r8*umg(i))**2._r8+  &
        3037             :              0.08_r8*umg(i)*umr(i))**0.5_r8*rho(i)*                 &
        3038             :              n0r(i)*n0g(i)/lamr(i)**3._r8*                          &
        3039             :              (5._r8/(lamr(i)**3._r8*lamg(i))+                       &
        3040             :              2._r8/(lamr(i)**2._r8*lamg(i)**2._r8)+                 &
        3041           0 :              0.5_r8/(lamr(i)*lamg(i)**3._r8)))
        3042             :
        3043             : ! assume 1 mm drops are shed, get number shed per sec
        3044           0 :         dum = pracg(i)/5.2e-7_r8
        3045             :         npracg(i) = cons32*rho(i)*(1.7_r8*(unr(i)-ung(i))**2._r8+   &
        3046             :              0.3_r8*unr(i)*ung(i))**0.5_r8*n0r(i)*n0g(i)*           &
        3047             :              (1._r8/(lamr(i)**3._r8*lamg(i))+                       &
        3048             :              1._r8/(lamr(i)**2._r8*lamg(i)**2._r8)+                 &
        3049           0 :              1._r8/(lamr(i)*lamg(i)**3._r8))
        3050             :
        3051             : ! hm 7/15/13, remove limit so that the number of collected drops can smaller than
        3052             : ! number of shed drops
        3053             : !            NPRACG(K)=MAX(NPRACG(K)-DUM,0.)
        3054           0 :         npracg(i) = npracg(i) - dum
        3055             :      else
        3056           0 :         npracg(i) = 0._r8
        3057           0 :         pracg(i)  = 0._r8
        3058             :      end if
        3059             :   end do
        3060             :   !$acc end parallel
        3061             :
        3062             :   !$acc end data
        3063           0 : end subroutine graupel_collecting_rain
        3064             :
        3065             : !========================================================================
        3066             : ! Rain riming snow to graupel
        3067             : !========================================================================
        3068             : ! Conversion of rimed rainwater onto snow converted to graupel
        3069             :
        3070           0 : subroutine graupel_rain_riming_snow(pracs,npracs,psacr,qsic,qric,nric,nsic,n0s,lams,n0r,lamr,dtime,&
        3071           0 :      pgracs,ngracs,vlen)
        3072             :
        3073             :   integer,                   intent(in)    :: vlen
        3074             :   ! Accretion of rain by snow
        3075             :   real(r8), dimension(vlen), intent(inout) :: pracs
        3076             :   real(r8), dimension(vlen), intent(inout) :: npracs
        3077             :   real(r8), dimension(vlen), intent(inout) :: psacr  ! conversion due to coll of snow by rain
        3078             :   !MMR
        3079             :   real(r8), dimension(vlen), intent(in)    :: qsic   ! snow mixing ratio
        3080             :   real(r8), dimension(vlen), intent(in)    :: qric   ! rain mixing ratio
        3081             :   real(r8), dimension(vlen), intent(in)    :: nric   ! rain number concentration
        3082             :   real(r8), dimension(vlen), intent(in)    :: nsic   ! snow number concentration
        3083             :   ! Size parameters for snow
        3084             :   real(r8), dimension(vlen), intent(in)    :: n0s
        3085             :   real(r8), dimension(vlen), intent(in)    :: lams
        3086             :   ! Size parameters for rain
        3087             :   real(r8), dimension(vlen), intent(in)    :: n0r
        3088             :   real(r8), dimension(vlen), intent(in)    :: lamr
        3089             :   real(r8),                  intent(in)    :: dtime
        3090             :   !Output process rates
        3091             :   real(r8), dimension(vlen), intent(out)   :: pgracs  ! Q graupel due to collection rain by snow
        3092             :   real(r8), dimension(vlen), intent(out)   :: ngracs  ! N graupel due to collection rain by snow
        3093             :
        3094             : !Input: PRACS,NPRACS,PSACR,qni,qr,lams,lamr,nr,ns  Note: No PSACR in MG2
        3095             : !Output:PGRACS,NGRACS,PRACS,PSACR
        3096             :   integer  :: i
        3097             :   real(r8), parameter :: cons18 = rhosn*rhosn
        3098             :   real(r8), parameter :: cons19 = rhow*rhow
        3099             :   real(r8) :: dum
        3100             :
        3101             :   !$acc data present (pracs,npracs,psacr,qsic,qric,nric)  &
        3102             :   !$acc      present (nsic,n0s,lams,n0r,lamr,pgracs,ngracs)
        3103             :
        3104             :   !$acc parallel vector_length(VLENS) default(present)
        3105             :   !$acc loop gang vector private(dum)
        3106           0 :   do i=1,vlen
        3107           0 :      if (pracs(i).gt.0._r8.and.qsic(i).ge.0.1e-3_r8.and.qric(i).ge.0.1e-3_r8) then
        3108             :         ! only allow conversion if qs > 0.1 and qr > 0.1 g/kg following rutledge and hobbs (1984)
        3109             :         !if (qsic(i).ge.0.1e-3_r8.and.qric(i).ge.0.1e-3_r8) then
        3110             :            ! portion of collected rainwater converted to graupel (reisner et al. 1998)
        3111             :         dum = cons18*(4._r8/lams(i))**3*(4._r8/lams(i))**3 &
        3112             :              /(cons18*(4._r8/lams(i))**3*(4._r8/lams(i))**3+ &
        3113           0 :              cons19*(4._r8/lamr(i))**3*(4._r8/lamr(i))**3)
        3114           0 :         dum = min(dum,1._r8)
        3115           0 :         dum = max(dum,0._r8)
        3116           0 :         pgracs(i) = (1._r8-dum)*pracs(i)
        3117           0 :         ngracs(i) = (1._r8-dum)*npracs(i)
        3118             :         ! limit max number converted to min of either rain or snow number concentration
        3119           0 :         ngracs(i) = min(ngracs(i),nric(i)/dtime)
        3120           0 :         ngracs(i) = min(ngracs(i),nsic(i)/dtime)
        3121             :
        3122             :         ! amount left for snow production
        3123           0 :         pracs(i)  = pracs(i) - pgracs(i)
        3124           0 :         npracs(i) = npracs(i) - ngracs(i)
        3125             :
        3126             :         ! conversion to graupel due to collection of snow by rain
        3127           0 :         psacr(i)=psacr(i)*(1._r8-dum)
        3128             :      else
        3129           0 :         pgracs(i) = 0._r8
        3130           0 :         ngracs(i) = 0._r8
        3131             :      end if
        3132             :   end do
        3133             :   !$acc end parallel
        3134             :
        3135             :   !$acc end data
        3136           0 : end subroutine graupel_rain_riming_snow
        3137             :
        3138             : !========================================================================
        3139             : ! Rime Splintering
        3140             : !========================================================================
        3141           0 : subroutine graupel_rime_splintering(t,qcic,qric,qgic,psacwg,pracg,&
        3142           0 :      qmultg,nmultg,qmultrg,nmultrg,vlen)
        3143             :
        3144             :   integer,                   intent(in)    :: vlen
        3145             :   real(r8), dimension(vlen), intent(in)    :: t        ! temperature
        3146             :   !MMR
        3147             :   real(r8), dimension(vlen), intent(in)    :: qcic     ! liquid mixing ratio
        3148             :   real(r8), dimension(vlen), intent(in)    :: qric     ! rain mixing ratio
        3149             :   real(r8), dimension(vlen), intent(in)    :: qgic     ! graupel mixing ratio
        3150             :   ! Already calculated terms for collection
        3151             :   real(r8), dimension(vlen), intent(inout) :: psacwg   ! collection droplets by graupel
        3152             :   real(r8), dimension(vlen), intent(inout) :: pracg    ! collection rain by graupel
        3153             :   !Output process rates for splintering
        3154             :   real(r8), dimension(vlen), intent(out)   :: qmultg   ! Q ice mult of droplets/graupel
        3155             :   real(r8), dimension(vlen), intent(out)   :: nmultg   ! N ice mult of droplets/graupel
        3156             :   real(r8), dimension(vlen), intent(out)   :: qmultrg  ! Q ice mult of rain/graupel
        3157             :   real(r8), dimension(vlen), intent(out)   :: nmultrg  ! N ice mult of rain/graupel
        3158             :
        3159             : !Input: qg,qc,qr, PSACWG,PRACG,T
        3160             : !Output: NMULTG,QMULTG,NMULTRG,QMULTRG,PSACWG,PRACG
        3161             :   integer  :: i
        3162             :   real(r8) :: fmult
        3163             :   real(r8) :: tm_3, tm_5, tm_8
        3164             :
        3165           0 :   tm_3 = tmelt - 3._r8
        3166           0 :   tm_5 = tmelt - 5._r8
        3167           0 :   tm_8 = tmelt - 8._r8
        3168             :
        3169             :   !$acc data present (t,qcic,qric,qgic,psacwg,pracg) &
        3170             :   !$acc      present (qmultg,nmultg,qmultrg,nmultrg)
        3171             :
        3172             : !nmultg,qmultg                                                                             .
        3173             : !========================================================================
        3174             :
        3175             :   !$acc parallel vector_length(VLENS) default(present)
        3176             :   !$acc loop gang vector private(fmult)
        3177           0 :   do i=1,vlen
        3178           0 :      nmultrg(i) = 0._r8
        3179           0 :      qmultrg(i) = 0._r8
        3180           0 :      nmultg(i)  = 0._r8
        3181           0 :      qmultg(i)  = 0._r8
        3182           0 :      if (qgic(i).ge.0.1e-3_r8) then
        3183           0 :         if (qcic(i).ge.0.5e-3_r8.or.qric(i).ge.0.1e-3_r8) then
        3184           0 :            if (psacwg(i).gt.0._r8.or.pracg(i).gt.0._r8) then
        3185           0 :               if (t(i).lt.tm_3 .and. t(i).gt.tm_8) then
        3186           0 :                  if (t(i).gt.tm_3) then
        3187             :                     fmult = 0._r8
        3188           0 :                  else if (t(i).le.tm_3.and.t(i).gt.tm_5)  then
        3189           0 :                     fmult = (tm_3-t(i))/2._r8
        3190           0 :                  else if (t(i).ge.tm_8.and.t(i).le.tm_5)   then
        3191           0 :                     fmult = (t(i)-tm_8)/3._r8
        3192           0 :                  else if (t(i).lt.tm_8) then
        3193           0 :                     fmult = 0._r8
        3194             :                  end if
        3195             : ! 1000 is to convert from kg to g  (Check if needed for MG)
        3196             : ! splintering from droplets accreted onto graupel
        3197           0 :                  if (psacwg(i).gt.0._r8) then
        3198           0 :                     nmultg(i) = 35.e4_r8*psacwg(i)*fmult*1000._r8
        3199           0 :                     qmultg(i) = nmultg(i)*mmult
        3200             : ! constrain so that transfer of mass from graupel to ice cannot be more mass
        3201             : ! than was rimed onto graupel
        3202           0 :                     qmultg(i) = min(qmultg(i),psacwg(i))
        3203           0 :                     psacwg(i) = psacwg(i)-qmultg(i)
        3204             :                  end if
        3205             :
        3206             : !nmultrg,qmultrg                                                                             .
        3207             : !========================================================================
        3208             : ! riming and splintering from accreted raindrops
        3209             : ! Factor of 1000. again (Check)
        3210           0 :                  if (pracg(i).gt.0._r8) then
        3211           0 :                     nmultrg(i) = 35.e4_r8*pracg(i)*fmult*1000._r8
        3212           0 :                     qmultrg(i) = nmultrg(i)*mmult
        3213             : ! constrain so that transfer of mass from graupel to ice cannot be more mass
        3214             : ! than was rimed onto graupel
        3215           0 :                     qmultrg(i) = min(qmultrg(i),pracg(i))
        3216           0 :                     pracg(i)   = pracg(i)-qmultrg(i)
        3217             :                  end if
        3218             :               end if
        3219             :            end if
        3220             :         end if
        3221             :      end if
        3222             :   end do
        3223             :   !$acc end parallel
        3224             :
        3225             :   !$acc end data
        3226           0 : end subroutine graupel_rime_splintering
        3227             :
        3228             : !========================================================================
        3229             : ! Evaporation and sublimation of graupel
        3230             : !========================================================================
        3231             :
        3232             : !MERGE WITH RAIN AND SNOW EVAP
        3233             : !
        3234             : !subroutine graupel_sublimate_evap(t,q,qgic,rho,n0g,lamg,qvi,dv,mu,sc,bg,agn,&
        3235             : !     prdg,eprdg,vlen)
        3236             : !
        3237             : !  integer, intent(in) :: vlen
        3238             : !
        3239             : !  real(r8), dimension(vlen), intent(in) :: t  !temperature
        3240             : !  real(r8), dimension(vlen), intent(in) :: q  !specific humidity (mixing ratio)
        3241             : !
        3242             : !  !MMR
        3243             : !  real(r8), dimension(vlen), intent(in) :: qgic  !graupel mixing ratio
        3244             : !
        3245             : !  real(r8), dimension(vlen), intent(in) :: rho   ! air density
        3246             : !
        3247             : ! ! Size parameters for graupel
        3248             : !  real(r8), dimension(vlen), intent(in) :: n0g
        3249             : !  real(r8), dimension(vlen), intent(in) :: lamg
        3250             : !
        3251             : !  real(r8), dimension(vlen), intent(in) :: qvi  !saturation humidity (ice)
        3252             : !
        3253             : !  real(r8), dimension(vlen), intent(in) :: dv   ! water vapor diffusivity
        3254             : !  real(r8), dimension(vlen), intent(in) :: mu   ! viscosity
        3255             : !  real(r8), dimension(vlen), intent(in) :: sc   ! schmidt number
        3256             : !
        3257             : !  ! fallspeed parameters for graupel
        3258             : !  ! Set AGN and BG  as input (in micro_pumas_v1.F90)  (select hail or graupel as appropriate)
        3259             : !  real(r8),                    intent(in) :: bg
        3260             : !  real(r8), dimension(vlen), intent(in) :: agn
        3261             : !
        3262             : !  ! Output tendencies (sublimation or evaporation of graupel)
        3263             : !  real(r8), dimension(vlen), intent(out) :: prdg
        3264             : !  real(r8), dimension(vlen), intent(out) :: eprdg
        3265             : !
        3266             : !  real(r8) :: cons11,cons36
        3267             : !  real(r8) :: abi
        3268             : !  real(r8) :: epsg
        3269             : !  integer :: i
        3270             : !
        3271             : !  cons11=gamma(2.5_r8+bg/2._r8)  !bg will be different for graupel (bg) or hail (bh)
        3272             : !  cons36=(2.5_r8+bg/2._r8)
        3273             : !
        3274             : !
        3275             : !  do i=1,vlen
        3276             : !
        3277             : !     abi = calc_ab(t(i), qvi(i), xxls)
        3278             : !
        3279             : !      if (qgic(i).ge.qsmall) then
        3280             : !         epsg = 2._r8*pi*n0g(i)*rho(i)*dv(i)*                                &
        3281             : !              (f1s/(lamg(i)*lamg(i))+                               &
        3282             : !              f2s*(agn(i)*rho(i)/mu(i))**0.5_r8*                      &
        3283             : !              sc(i)**(1._r8/3._r8)*cons11/                   &
        3284             : !              (lamg(i)**cons36))
        3285             : !      else
        3286             : !         epsg = 0.
        3287             : !      end if
        3288             : !
        3289             : !! vapor dpeosition on graupel
        3290             : !      prdg(i) = epsg*(q(i)-qvi(i))/abi
        3291             : !
        3292             : !! make sure not pushed into ice supersat/subsat
        3293             : !! put this in main mg3 code…..check for it…
        3294             : !! formula from reisner 2 scheme
        3295             : !!
        3296             : !!          dum = (qv3d(k)-qvi(k))/dt
        3297             : !!
        3298             : !!           fudgef = 0.9999
        3299             : !!           sum_dep = prd(k)+prds(k)+mnuccd(k)+prdg(k)
        3300             : !!
        3301             : !!           if( (dum.gt.0. .and. sum_dep.gt.dum*fudgef) .or.                      &
        3302             : !!               (dum.lt.0. .and. sum_dep.lt.dum*fudgef) ) then
        3303             : !!             prdg(k) = fudgef*prdg(k)*dum/sum_dep
        3304             : !!           endif
        3305             : !
        3306             : !! if cloud ice/snow/graupel vap deposition is neg, then assign to sublimation processes
        3307             : !
        3308             : !      eprdg(i)=0._r8
        3309             : !
        3310             : !      if (prdg(i).lt.0._r8) then
        3311             : !         eprdg(i)=prdg(i)
        3312             : !         prdg(i)=0.
        3313             : !      end if
        3314             : !
        3315             : !   enddo
        3316             : !
        3317             : !end subroutine graupel_sublimate_evap
        3318             :
        3319             : !========================================================================
        3320             : !MG4 Lookup Table
        3321             : !========================================================================
        3322             :
        3323           0 : SUBROUTINE init_lookup_table(lkuptable_filename)
        3324             :
        3325             :   implicit none
        3326             :     !character(100) :: lkuptable_filename = "/home/katec/mg4work/lookup_table.dat"
        3327             :     !character(100) :: lkuptable_filename = "/glade/p/work/katec/mg3work/lookup_table.dat"
        3328             :     !character(100) :: lkuptable_filename = "./lookup_table.dat"
        3329             :     character(len=256), intent(in) :: lkuptable_filename
        3330             :     integer :: i,j,k,ii,jj,kk, tt
        3331             :     real(r8)    :: dum
        3332             :
        3333             :     ! ++ trude
        3334           0 :     open(unit=10,file=lkuptable_filename,status='old')
        3335             :
        3336             :     !------------------------------------------------------------------------------------------!
        3337             :     ! read in ice microphysics table
        3338           0 :      do tt = 1, tsize
        3339           0 :        do i = 1,isize
        3340           0 :           do k = 1,jsize
        3341           0 :              read(10,*) dum,dum,itab(tt,i,k,1),itab(tt,i,k,2),             &
        3342           0 :                   itab(tt,i,k,3),itab(tt,i,k,4),itab(tt,i,k,5),            &
        3343           0 :                   itab(tt,i,k,6),itab(tt,i,k,7),itab(tt,i,k,8),            &
        3344           0 :                   itab(tt,i,k,13),itab(tt,i,k,9),itab(tt,i,k,10),          &
        3345           0 :                   itab(tt,i,k,11),itab(tt,i,k,12),itab(tt,i,k,14)
        3346             :             enddo
        3347             :        enddo
        3348             :        ! read in table for ice-rain collection
        3349           0 :        do i = 1,isize
        3350           0 :           do k = 1,jsize
        3351           0 :              do j = 1,rcollsize
        3352           0 :                 read(10,*) dum,dum,dum,itabcoll(tt,i,k,j,1),                  &
        3353           0 :                      itabcoll(tt,i,k,j,2),dum
        3354           0 :                 itabcoll(tt,i,k,j,1) = dlog10(itabcoll(tt,i,k,j,1))
        3355           0 :                 itabcoll(tt,i,k,j,2) = dlog10(itabcoll(tt,i,k,j,2))
        3356             :              end do
        3357             :           end do
        3358             :        end do
        3359             :     end do
        3360             :
        3361           0 :     close(unit=10)
        3362             :
        3363           0 : END SUBROUTINE init_lookup_table
        3364             :
        3365           0 : pure SUBROUTINE access_lookup_table(dumii,dumi,dumk,index,dum1,dum2,dum4,proc)
        3366             :
        3367             :   implicit none
        3368             :
        3369             :   real(r8), intent(in) :: dum1,dum2,dum4
        3370             :   real(r8), intent(out) :: proc
        3371             :   real(r8) :: dproc1,dproc2,iproc1,gproc1,tmp1,tmp2
        3372             :   integer, intent(in) :: dumi,dumk,dumii,index
        3373             :
        3374             :
        3375             :   !(trude,  dumii could be the same as dumt)
        3376             :   ! get value at current density index
        3377             :   !***** without temperature dependence *******
        3378             :   ! first interpolate for current rimed fraction index
        3379             :   !   dproc1 = itab(dumi,dumk,index)+(dum1-real(dumi))*(itab(       &
        3380             :   !            dumi+1,dumk,index)-itab(dumi,dumk,index))
        3381             :
        3382             :   !   dproc2 = itab(dumi,dumk+1,index)+(dum1-real(dumi))*(itab(     &
        3383             :   !          dumi+1,dumk+1,index)-itab(dumi,dumk+1,index))
        3384             :
        3385             :   !   iproc1 = dproc1+(dum2-real(dumk))*(dproc2-dproc1)
        3386             :
        3387             :   !   proc    = iproc1
        3388             :   !! ***** end without temperature dependence ***********
        3389             :   ! ******************* from orig P3 ********************
        3390             :
        3391             :   ! trude comment, dumjj is for density index, dumii is for rime, which is subsituted by temperature. All dumjj should be removed.
        3392             :   ! first interpolate for current temperature
        3393           0 :   dproc1 = itab(dumii,dumi,dumk,index)+(dum1-real(dumi))*(itab(dumii,       &
        3394           0 :        dumi+1,dumk,index)-itab(dumii,dumi,dumk,index))
        3395             :
        3396           0 :   dproc2 = itab(dumii,dumi,dumk+1,index)+(dum1-real(dumi))*(itab(dumii,     &
        3397           0 :        dumi+1,dumk+1,index)-itab(dumii,dumi,dumk+1,index))
        3398             :
        3399           0 :   iproc1 = dproc1+(dum2-real(dumk))*(dproc2-dproc1)
        3400             :
        3401             :   ! linearly interpolate to get process rates for temperature index + 1
        3402           0 :   dproc1 = itab(dumii+1,dumi,dumk,index)+(dum1-real(dumi))*(itab(dumii+1,   &
        3403           0 :        dumi+1,dumk,index)-itab(dumii+1,dumi,dumk,index))
        3404             :
        3405             :   dproc2 = itab(dumii+1,dumi,dumk+1,index)+(dum1-real(dumi))*(itab(dumii+1, &
        3406           0 :        dumi+1,dumk+1,index)-itab(dumii+1,dumi,dumk+1,index))
        3407             :
        3408           0 :   gproc1 = dproc1+(dum2-real(dumk))*(dproc2-dproc1)
        3409           0 :   tmp1   = iproc1+(dum4-real(dumii))*(gproc1-iproc1)
        3410           0 :   proc   = tmp1
        3411             :
        3412             :   ! ********** end from orig P3 **************
        3413             :
        3414           0 : END SUBROUTINE access_lookup_table
        3415             :
        3416             : !------------------------------------------------------------------------------------------!
        3417           0 : pure SUBROUTINE access_lookup_table_coll(dumii,dumj,dumi,dumk,index,dum1,dum2,dum3,     &
        3418             :      dum4,proc)
        3419             :
        3420             :   implicit none
        3421             :
        3422             :   real(r8),intent(in)    :: dum1,dum2,dum3,dum4
        3423             :   real(r8), intent(out) :: proc
        3424             :   real(r8) :: dproc1,dproc2,iproc1,gproc1,tmp1,tmp2,dproc11, &
        3425             :        dproc12,dproc21,dproc22
        3426             :   integer, intent(in) :: dumii,dumj,dumi,dumk,index
        3427             :
        3428             :
        3429             :   ! This subroutine interpolates lookup table values for rain/ice collection processes
        3430             :   !!******for without temperature dependeny *****
        3431             :   !   dproc11 = itabcoll(dumi,dumk,dumj,index)+(dum1-real(dumi))*               &
        3432             :   !             (itabcoll(dumi+1,dumk,dumj,index)-itabcoll(dumi,    &
        3433             :   !             dumk,dumj,index))
        3434             :
        3435             :   !   dproc21 = itabcoll(dumi,dumk+1,dumj,index)+(dum1-real(dumi))*             &
        3436             :   !             (itabcoll(dumi+1,dumk+1,dumj,index)-itabcoll(dumi,  &
        3437             :   !             dumk+1,dumj,index))
        3438             :
        3439             :   !   dproc1  = dproc11+(dum2-real(dumk))*(dproc21-dproc11)
        3440             :
        3441             :   !   dproc12 = itabcoll(dumi,dumk,dumj+1,index)+(dum1-real(dumi))*             &
        3442             :   !             (itabcoll(dumi+1,dumk,dumj+1,index)-itabcoll(dumi,  &
        3443             :   !             dumk,dumj+1,index))
        3444             :
        3445             :   !   dproc22 = itabcoll(dumi,dumk+1,dumj+1,index)+(dum1-real(dumi))*           &
        3446             :   !             (itabcoll(dumi+1,dumk+1,dumj+1,index)-itabcoll(     &
        3447             :   !             dumi,dumk+1,dumj+1,index))
        3448             :
        3449             :   !   dproc2  = dproc12+(dum2-real(dumk))*(dproc22-dproc12)
        3450             :   !   iproc1  = dproc1+(dum3-real(dumj))*(dproc2-dproc1)
        3451             :
        3452             :   !   proc    = iproc1
        3453             :   !!******end for without temperature dependeny *****
        3454             :
        3455             :   ! ************ from original P3 ************
        3456             :   ! current density index
        3457             :
        3458             :   ! current temp index
        3459           0 :   dproc11 = itabcoll(dumii,dumi,dumk,dumj,index)+(dum1-real(dumi))*               &
        3460           0 :        (itabcoll(dumii,dumi+1,dumk,dumj,index)-itabcoll(dumii,dumi,    &
        3461           0 :        dumk,dumj,index))
        3462             :
        3463           0 :   dproc21 = itabcoll(dumii,dumi,dumk+1,dumj,index)+(dum1-real(dumi))*             &
        3464             :        (itabcoll(dumii,dumi+1,dumk+1,dumj,index)-itabcoll(dumii,dumi,  &
        3465           0 :        dumk+1,dumj,index))
        3466             :
        3467           0 :   dproc1  = dproc11+(dum2-real(dumk))*(dproc21-dproc11)
        3468           0 :   dproc12 = itabcoll(dumii,dumi,dumk,dumj+1,index)+(dum1-real(dumi))*             &
        3469             :        (itabcoll(dumii,dumi+1,dumk,dumj+1,index)-itabcoll(dumii,dumi,  &
        3470           0 :        dumk,dumj+1,index))
        3471             :
        3472             :   dproc22 = itabcoll(dumii,dumi,dumk+1,dumj+1,index)+(dum1-real(dumi))*           &
        3473             :        (itabcoll(dumii,dumi+1,dumk+1,dumj+1,index)-itabcoll(dumii,     &
        3474           0 :        dumi,dumk+1,dumj+1,index))
        3475             :
        3476           0 :   dproc2  = dproc12+(dum2-real(dumk))*(dproc22-dproc12)
        3477           0 :   iproc1  = dproc1+(dum3-real(dumj))*(dproc2-dproc1)
        3478             :
        3479             :   ! temp index + 1
        3480             :
        3481           0 :   dproc11 = itabcoll(dumii+1,dumi,dumk,dumj,index)+(dum1-real(dumi))*             &
        3482             :        (itabcoll(dumii+1,dumi+1,dumk,dumj,index)-itabcoll(dumii+1,     &
        3483           0 :        dumi,dumk,dumj,index))
        3484             :
        3485             :   dproc21 = itabcoll(dumii+1,dumi,dumk+1,dumj,index)+(dum1-real(dumi))*           &
        3486             :        (itabcoll(dumii+1,dumi+1,dumk+1,dumj,index)-itabcoll(dumii+1,   &
        3487           0 :        dumi,dumk+1,dumj,index))
        3488             :
        3489           0 :   dproc1  = dproc11+(dum2-real(dumk))*(dproc21-dproc11)
        3490             :
        3491             :   dproc12 = itabcoll(dumii+1,dumi,dumk,dumj+1,index)+(dum1-real(dumi))*           &
        3492             :        (itabcoll(dumii+1,dumi+1,dumk,dumj+1,index)-itabcoll(dumii+1,   &
        3493           0 :        dumi,dumk,dumj+1,index))
        3494             :
        3495             :   dproc22 = itabcoll(dumii+1,dumi,dumk+1,dumj+1,index)+(dum1-real(dumi))*         &
        3496             :        (itabcoll(dumii+1,dumi+1,dumk+1,dumj+1,index)-itabcoll(dumii+1, &
        3497           0 :        dumi,dumk+1,dumj+1,index))
        3498             :
        3499           0 :   dproc2  = dproc12+(dum2-real(dumk))*(dproc22-dproc12)
        3500             :
        3501           0 :   gproc1  = dproc1+(dum3-real(dumj))*(dproc2-dproc1)
        3502           0 :   tmp1    = iproc1+(dum4-real(dumii))*(gproc1-iproc1)
        3503           0 :   proc    = tmp1
        3504             :
        3505             :   ! ********** end orig P3 *******
        3506           0 : END SUBROUTINE access_lookup_table_coll
        3507             :
        3508             : !========================================================================
        3509             : !UTILITIES
        3510             : !========================================================================
        3511             :
        3512           0 : pure function no_limiter()
        3513             :   real(r8) :: no_limiter
        3514             :
        3515           0 :   no_limiter = transfer(limiter_off, no_limiter)
        3516             :
        3517           0 : end function no_limiter
        3518             :
        3519           0 : pure function limiter_is_on(lim)
        3520             :   !$acc routine seq
        3521             :
        3522             :   real(r8), intent(in) :: lim
        3523             :   logical :: limiter_is_on
        3524             :
        3525           0 :   limiter_is_on = transfer(lim, limiter_off) /= limiter_off
        3526             :
        3527           0 : end function limiter_is_on
        3528             :
        3529           0 : end module micro_pumas_utils

Generated by: LCOV version 1.14