LCOV - code coverage report
Current view: | top level - physics/pumas - micro_pumas_utils.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 562 | 822 | 68.4 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 43 | 57 | 75.4 %

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

Generated by: LCOV version 1.14