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