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 0 : 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 0 : res = gamma(x+n)/gamma(x)
432 :
433 0 : 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 0 : 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 0 : res = 1._r8
465 0 : factor = x
466 :
467 : !$acc loop seq
468 0 : do i = 1, n
469 0 : res = res * factor
470 0 : factor = factor + 1._r8
471 : end do
472 :
473 0 : end subroutine rising_factorial_integer
474 :
475 772882344 : 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 42009204744 : do i=1,vlen
487 41236322400 : res(i) = 1._r8
488 41236322400 : factor = x(i)
489 :
490 42009204744 : 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 5890903200 : elseif (n == 2) then
497 5890903200 : res(i) = res(i) * factor
498 5890903200 : factor = factor + 1._r8
499 5890903200 : 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 772882344 : end subroutine rising_factorial_integer_vec
511 :
512 : ! Calculate correction due to latent heat for evaporation/sublimation
513 2649323946 : 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 2649323946 : dqsdt = xxl*qv / (rv * t**2)
523 2649323946 : ab = 1._r8 + dqsdt*xxl/cpp
524 :
525 2649323946 : 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 2622210112 : pgam(i,k) = 1.0_r8 - 0.7_r8 * exp(-0.008_r8*1.e-6_r8*ncic(i,k)*rho(i,k))
624 2622210112 : pgam(i,k) = 1._r8/(pgam(i,k)**2) - 1._r8
625 2622210112 : pgam(i,k) = max(pgam(i,k), 2._r8)
626 2622210112 : pgamp1(i,k) = pgam(i,k)+1._r8
627 : else
628 20941402688 : 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 2622210112 : shapeC(i,k) = pi / 6._r8 * props%rho * tmp(i,k)
648 : ! Limit to between 2 and 50 microns mean size.
649 2622210112 : lbnd(i,k) = pgamp1(i,k)*1._r8/50.e-6_r8
650 2622210112 : ubnd(i,k) = pgamp1(i,k)*1._r8/2.e-6_r8
651 : else
652 20941402688 : shapeC(i,k) = 0._r8
653 20941402688 : lbnd(i,k) = 0._r8
654 20941402688 : 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 20941402688 : pgam(i,k) = -100._r8
671 20941402688 : 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 1495329478 : pgam(i) = 1.0_r8 - 0.7_r8 * exp(-0.008_r8*1.e-6_r8*ncic(i)*rho(i))
705 1495329478 : pgam(i) = 1._r8/(pgam(i)**2) - 1._r8
706 1495329478 : pgam(i) = max(pgam(i), 2._r8)
707 1495329478 : pgamp1(i) = pgam(i)+1._r8
708 : else
709 10286476922 : 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 1495329478 : shapeC(i) = pi / 6._r8 * props%rho * tmp(i)
727 : ! Limit to between 2 and 50 microns mean size.
728 1495329478 : lbnd(i) = pgamp1(i)*1._r8/50.e-6_r8
729 1495329478 : ubnd(i) = pgamp1(i)*1._r8/2.e-6_r8
730 : else
731 10286476922 : shapeC(i) = 0._r8
732 10286476922 : lbnd(i) = 0._r8
733 10286476922 : 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 10286476922 : pgam(i) = -100._r8
748 10286476922 : 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 508303159 : 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 508303159 : present_n0 = present(n0)
767 :
768 508303159 : if (qic > qsmall) then
769 : ! add upper limit to in-cloud number concentration to prevent
770 : ! numerical error
771 508303159 : 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 508303159 : lam = (props%shape_coef * nic/qic)**(1._r8/props%eff_dim)
777 :
778 : ! check for slope
779 : ! adjust vars
780 508303159 : if (lam < props%lambda_bounds(1)) then
781 778871 : lam = props%lambda_bounds(1)
782 778871 : nic = lam**(props%eff_dim) * qic/props%shape_coef
783 507524288 : else if (lam > props%lambda_bounds(2)) then
784 2299 : lam = props%lambda_bounds(2)
785 2299 : nic = lam**(props%eff_dim) * qic/props%shape_coef
786 : end if
787 : else
788 0 : lam = 0._r8
789 : end if
790 :
791 508303159 : if (present_n0) n0 = nic * lam
792 :
793 508303159 : 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 3157427111 : if (limiterActive) then
822 3157427111 : nic(i) = min(nic(i), qic(i) / minMass)
823 : end if
824 :
825 : ! lambda = (c n/q)^(1/d)
826 3157427111 : lam(i) = (shapeCoef * nic(i)/qic(i))**(1._r8/effDim)
827 :
828 : ! check for slope
829 : ! adjust vars
830 3157427111 : if (lam(i) < lbnd) then
831 0 : lam(i) = lbnd
832 0 : nic(i) = lam(i)**(effDim) * qic(i)/shapeCoef
833 3157427111 : else if (lam(i) > ubnd) then
834 13174539 : lam(i) = ubnd
835 13174539 : nic(i) = lam(i)**(effDim) * qic(i)/shapeCoef
836 : end if
837 : else
838 8624379289 : 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 19101136092 : if (limiterActive) then
875 5091714925 : nic(i,k) = min(nic(i,k), qic(i,k) / minMass)
876 : end if
877 :
878 : ! lambda = (c n/q)^(1/d)
879 19101136092 : lam(i,k) = (shapeCoef * nic(i,k)/qic(i,k))**(1._r8/effDim)
880 :
881 : ! check for slope
882 : ! adjust vars
883 19101136092 : if (lam(i,k) < lbnd) then
884 1363261480 : lam(i,k) = lbnd
885 1363261480 : nic(i,k) = lam(i,k)**(effDim) * qic(i,k)/shapeCoef
886 17737874612 : else if (lam(i,k) > ubnd) then
887 1863527578 : lam(i,k) = ubnd
888 1863527578 : nic(i,k) = lam(i,k)**(effDim) * qic(i,k)/shapeCoef
889 : end if
890 : else
891 63371508708 : 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 4117539590 : if (limiterActive) then
929 4117539590 : nic(i) = min(nic(i), qic(i) / minMass)
930 : end if
931 : ! lambda = (c n/q)^(1/d)
932 :
933 4117539590 : lam(i) = (shapeC(i) * nic(i)/qic(i))**(1._r8/effDim)
934 : ! check for slope
935 : ! adjust vars
936 :
937 4117539590 : if (lam(i) < lbnd(i)) then
938 571520655 : lam(i) = lbnd(i)
939 571520655 : nic(i) = lam(i)**(effDim) * qic(i)/shapeC(i)
940 3546018935 : else if (lam(i) > ubnd(i)) then
941 895048242 : lam(i) = ubnd(i)
942 895048242 : nic(i) = lam(i)**(effDim) * qic(i)/shapeC(i)
943 : end if
944 : else
945 31227879610 : 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 1369952840 : 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 1369952840 : avg_diameter = (q*rho_air/(pi * rho_sub * n))**(1._r8/3._r8)
964 :
965 1369952840 : 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 0 : 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 0 : call rising_factorial(relvar, a, res)
997 0 : res = res / relvar**a
998 :
999 0 : 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 0 : 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 0 : real(r8) :: tmp(vlen)
1048 :
1049 : !$acc data create (tmp)
1050 :
1051 0 : call rising_factorial(relvar, a,tmp,vlen)
1052 :
1053 : !$acc parallel vector_length(VLENS) default(present)
1054 : !$acc loop gang vector
1055 0 : do i=1,vlen
1056 0 : res(i) = tmp(i) / relvar(i)**a
1057 : end do
1058 : !$acc end parallel
1059 :
1060 : !$acc end data
1061 0 : 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 4528047249 : qiic(i) = 0._r8
1119 4528047249 : niic(i) = 0._r8
1120 4528047249 : 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 1362855951 : epsi = 2._r8*pi*n0i(i)*rho(i)*Dv(i)/(lami(i)*lami(i))
1134 :
1135 : !Compute deposition/sublimation
1136 1362855951 : vap_dep(i) = epsi/ab(i)*(qv(i) - qvi(i))
1137 :
1138 : !Make this a grid-averaged quantity
1139 1362855951 : vap_dep(i) = vap_dep(i)*icldm(i)
1140 :
1141 : !Split into deposition or sublimation.
1142 1362855951 : if (t(i) < tmelt .and. vap_dep(i) > 0._r8) then
1143 557245533 : ice_sublim(i) = 0._r8
1144 : else
1145 : ! make ice_sublim negative for consistency with other evap/sub processes
1146 805610418 : ice_sublim(i) = min(vap_dep(i),0._r8)
1147 805610418 : vap_dep(i) = 0._r8
1148 : end if
1149 :
1150 : !sublimation occurs @ any T. Not so for berg.
1151 1362855951 : if (t(i) < tmelt) then
1152 :
1153 : !Compute bergeron rate assuming cloud for whole step.
1154 1235238825 : berg(i) = max(epsi/ab(i)*(qvl(i) - qvi(i)), 0._r8)
1155 : else !T>frz
1156 127617126 : berg(i)=0._r8
1157 : end if !T<frz
1158 : else !where qi<qsmall
1159 4528047249 : berg(i) = 0._r8
1160 4528047249 : vap_dep(i) = 0._r8
1161 4528047249 : 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 508303159 : qsic = qs(i)/precip_frac(i)
1310 508303159 : nsic = ns(i)/precip_frac(i)
1311 :
1312 : !Compute linearized condensational heating correction
1313 508303159 : call calc_ab(t(i), qvi(i), xxls, ab)
1314 : !Get slope and intercept of gamma distn for ice.
1315 508303159 : 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 508303159 : (lams**(5._r8/2._r8+bs/2._r8)))
1322 :
1323 : !Compute deposition/sublimation
1324 508303159 : vap_deps(i) = eps/ab*(qv(i) - qvi(i))
1325 :
1326 : !Make this a grid-averaged quantity
1327 508303159 : vap_deps(i)=vap_deps(i)*precip_frac(i)
1328 :
1329 : !Only want deposition
1330 508303159 : if (t(i) < tmelt .and. vap_deps(i)>0._r8) then
1331 449445060 : vap_deps(i)=max(vap_deps(i),0._r8)
1332 : else
1333 58858099 : 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 536556896 : (ncic(i)*1.e-6_r8*rho(i))**(micro_mg_autocon_nd_exp)
1398 :
1399 536556896 : nprc(i) = prc(i) * (1._r8/droplet_mass_25um)
1400 536556896 : nprc1(i) = prc(i)*ncic(i)/qcic(i)
1401 :
1402 : else
1403 5354346304 : prc(i) = 0._r8
1404 5354346304 : nprc(i) = 0._r8
1405 5354346304 : 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 1235238825 : d_rat = lami(i)*dcs
1560 :
1561 : ! Rate of ice particle conversion (number).
1562 1235238825 : nprci(i) = n0i(i)/(lami(i)*ac_time)*exp(-d_rat)
1563 1235238825 : 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 1235238825 : (((d_rat + 3._r8)*d_rat + 6._r8)*d_rat + 6._r8)
1570 : else
1571 4655664375 : prci(i) = 0._r8
1572 4655664375 : 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 0 : subroutine immersion_freezing(microp_uniform, t, pgam, lamc, &
1583 0 : 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 0 : real(r8), dimension(vlen) :: dum
1608 : integer :: i
1609 : real(r8) :: tmp
1610 :
1611 : !$acc data create (dum)
1612 :
1613 0 : if (.not. microp_uniform) then
1614 0 : 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 0 : do i=1,vlen
1627 0 : if (qcic(i) >= qsmall .and. t(i) < 269.15_r8) then
1628 0 : call rising_factorial(pgam(i)+1._r8, 3, tmp)
1629 : nnuccc(i) = &
1630 : pi/6._r8*ncic(i)*tmp* &
1631 0 : bimm*(exp(aimm*(tmelt - t(i)))-1._r8)/lamc(i)**3
1632 :
1633 0 : call rising_factorial(pgam(i)+4._r8, 3, tmp)
1634 : mnuccc(i) = dum(i) * nnuccc(i) * &
1635 : pi/6._r8*rhow* &
1636 0 : tmp/lamc(i)**3
1637 : else
1638 0 : mnuccc(i) = 0._r8
1639 0 : 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 0 : 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 0 : subroutine contact_freezing (microp_uniform, t, p, rndst, nacon, &
1652 0 : 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 0 : real(r8) :: nslip(size(rndst,2)) ! slip correction factors
1685 0 : 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 0 : do i = 1,vlen
1700 0 : if (qcic(i) >= qsmall .and. t(i) < 269.15_r8) then
1701 0 : if (microp_uniform) then
1702 0 : dum = 1._r8
1703 0 : dum1 = 1._r8
1704 : else
1705 0 : call var_coef(relvar(i), 4._r8/3._r8, dum)
1706 0 : call var_coef(relvar(i), 1._r8/3._r8, dum1)
1707 : end if
1708 :
1709 0 : tcnt=(270.16_r8-t(i))**1.3_r8
1710 0 : 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 0 : (p(i)*sqrt( 8.0_r8*28.96e-3_r8/(pi*8.314409_r8*t(i)) ))
1713 :
1714 0 : do j = 1, mdust
1715 : ! Note that these two are vectors.
1716 0 : nslip(j) = 1.0_r8+(mfp/rndst(i,j))*(1.257_r8+(0.4_r8*exp(-(1.1_r8*rndst(i,j)/mfp))))! Slip correction factor
1717 :
1718 0 : 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 0 : contact_factor = dot_product(ndfaer,nacon(i,:)*tcnt) * pi * &
1722 0 : ncic(i) * (pgam(i) + 1._r8) / lamc(i)
1723 0 : call rising_factorial(pgam(i)+2._r8, 3, tmp)
1724 0 : mnucct(i) = dum * contact_factor * pi/3._r8*rhow*tmp/lamc(i)**3
1725 0 : nnucct(i) = dum1 * 2._r8 * contact_factor
1726 : else
1727 0 : mnucct(i) = 0._r8
1728 0 : 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 0 : 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 1544718626 : ((qsic(i)/nsic(i))*(1._r8/(rhosn*pi)))**((bs-1._r8)/3._r8)
1765 : else
1766 4346184574 : 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 190178600 : dc0 = (pgam(i)+1._r8)/lamc(i)
1865 190178600 : dum = dc0*dc0*uns(i)*rhow*lams(i)/(9._r8*mu(i))
1866 190178600 : eci = dum*dum/((dum+0.4_r8)*(dum+0.4_r8))
1867 190178600 : eci = max(eci,0._r8)
1868 190178600 : eci = min(eci,1._r8)
1869 :
1870 : ! no impact of sub-grid distribution of qc since psacws
1871 : ! is linear in qc
1872 190178600 : accrete_rate = pi/4._r8*asn(i)*rho(i)*n0s(i)*eci*gamma_bs_plus3 / lams(i)**(bs+3._r8)
1873 190178600 : psacws(i) = accrete_rate*qcic(i)
1874 190178600 : npsacws(i) = accrete_rate*ncic(i)
1875 : else
1876 5700724600 : psacws(i) = 0._r8
1877 5700724600 : 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 74355375 : nsacwi(i) = 3.5e8_r8*(270.16_r8-t(i))/2.0_r8*psacws(i)
1944 5816547825 : else if((t(i) < 268.16_r8) .and. (t(i) >= 265.16_r8)) then
1945 111011270 : nsacwi(i) = 3.5e8_r8*(t(i)-265.16_r8)/3.0_r8*psacws(i)
1946 : else
1947 5705536555 : 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 240753561 : common_factor = pi*ecr*rho(i)*n0r(i)*n0s(i)/(lamr(i)**3 * lams(i))
2000 240753561 : 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 240753561 : ((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 240753561 : ((d_rat + 1._r8)*d_rat + 1._r8)
2007 : else
2008 5650149639 : pracs(i) = 0._r8
2009 5650149639 : 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 414625232 : (exp(aimm*(tmelt - t(i)))-1._r8)/lamr(i)**3
2078 414625232 : mnuccr(i) = nnuccr(i) * 20._r8*pi*rhow/lamr(i)**3
2079 : else
2080 5476277968 : mnuccr(i) = 0._r8
2081 5476277968 : 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 411116956 : pra(i) = pra_coef(i) * 67._r8*(qcic(i)*qric(i))**1.15_r8
2138 411116956 : npra(i) = pra(i)*ncic(i)/qcic(i)
2139 : else
2140 5479786244 : pra(i) = 0._r8
2141 5479786244 : 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 1513923157 : nragg(i) = -8._r8*nric(i)*qric(i)*rho(i)
2169 : else
2170 4376980043 : 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 1180121467 : lams(i)**(bs+3._r8)
2209 1180121467 : prai(i) = accrete_rate * qiic(i)
2210 1180121467 : nprai(i) = accrete_rate * niic(i)
2211 : else
2212 4710781733 : prai(i) = 0._r8
2213 4710781733 : 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 4810229004 : dum(i) = 0._r8
2512 : else
2513 1080674196 : 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 4810229004 : if (qric(i) >= qsmall .or. qsic(i) >= qsmall .or. qgic(i) >= qsmall) then
2524 1461752642 : am_evp_st(i) = precip_frac(i) - dum(i)
2525 : ! calculate q for out-of-cloud region
2526 1461752642 : qclr=(q(i)-dum(i)*qvl(i))/(1._r8-dum(i))
2527 : end if
2528 :
2529 : ! evaporation of rain
2530 4810229004 : if (qric(i) >= qsmall) then
2531 777846309 : if (.not.evap_rhthrsh.or.(evap_rhthrsh.and.qclr/qvl(i).le.0.9_r8)) then
2532 777846309 : 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 777846309 : (lamr(i)**(5._r8/2._r8+br/2._r8)))
2538 777846309 : pre(i) = eps*(qclr-qvl(i))/abr
2539 : ! only evaporate in out-of-cloud region
2540 : ! and distribute across precip_frac
2541 777846309 : pre(i) = min(pre(i)*am_evp_st(i),0._r8)
2542 777846309 : pre(i) = pre(i)/precip_frac(i)
2543 : else
2544 0 : pre(i) = 0._r8
2545 : end if
2546 : else
2547 4032382695 : pre(i) = 0._r8
2548 : end if
2549 :
2550 : ! sublimation of snow
2551 4810229004 : if (qsic(i) >= qsmall) then
2552 863465824 : 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 863465824 : (lams(i)**(5._r8/2._r8+bs/2._r8)))
2558 863465824 : prds(i) = eps*(qclr-qvi(i))/ab
2559 : ! only sublimate in out-of-cloud region and distribute over precip_frac
2560 863465824 : prds(i) = min(prds(i)*am_evp_st(i),0._r8)
2561 863465824 : prds(i) = prds(i)/precip_frac(i)
2562 : else
2563 3946763180 : prds(i) = 0._r8
2564 : end if
2565 :
2566 : ! ADD GRAUPEL, do Same with prdg.
2567 4810229004 : if (qgic(i).ge.qsmall) then
2568 309530054 : 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 309530054 : (lamg(i)**(5._r8/2._r8+bs/2._r8)))
2574 309530054 : prdg(i) = eps*(qclr-qvi(i))/abg
2575 : ! only sublimate in out-of-cloud region and distribute over precip_frac
2576 309530054 : prdg(i) = min(prdg(i)*am_evp_st(i),0._r8)
2577 309530054 : prdg(i) = prdg(i)/precip_frac(i)
2578 : else
2579 4500698950 : prdg(i) = 0._r8
2580 : end if
2581 :
2582 : else
2583 1080674196 : prds(i) = 0._r8
2584 1080674196 : pre(i) = 0._r8
2585 1080674196 : 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 190178600 : 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 190178600 : (lams(i)**(5._r8/2._r8+bs/2._r8)))
2751 190178600 : bergs(i) = eps*(qvl(i)-qvi(i))/ab
2752 : else
2753 5700724600 : 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 35001461 : 0.5_r8/(lams(i)*lamr(i)**3)))
2796 : else
2797 5855901739 : 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 52339648 : lamg(i)**(bg+3._r8)
2840 : npsacwg(i) = cons*agn(i)*ncic(i)*rho(i)* &
2841 : n0g(i)/ &
2842 52339648 : lamg(i)**(bg+3._r8)
2843 : else
2844 5838563552 : psacwg(i) = 0._r8
2845 5838563552 : 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 4783776 : (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 4783776 : dum = max(rhosn/(rhog-rhosn)*pgsacw(i),0._r8)
2905 :
2906 : ! Number concentraiton of embryo graupel from riming of snow
2907 4783776 : nscng(i) = dum/mg0*rho(i)
2908 : ! Limit max number converted to snow number (dtime here correct? We think yes)
2909 4783776 : nscng(i) = min(nscng(i),nsic(i)/dtime)
2910 :
2911 : ! Portion of riming left for snow
2912 4783776 : psacws(i) = psacws(i) - pgsacw(i)
2913 : else
2914 5886119424 : pgsacw(i) = 0._r8
2915 5886119424 : 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 116432012 : 0.5_r8/(lamr(i)*lamg(i)**3._r8)))
2969 :
2970 : ! assume 1 mm drops are shed, get number shed per sec
2971 116432012 : 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 116432012 : 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 116432012 : npracg(i) = npracg(i) - dum
2982 : else
2983 5774471188 : npracg(i) = 0._r8
2984 5774471188 : 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 20439644 : cons19*(4._r8/lamr(i))**3*(4._r8/lamr(i))**3)
3037 20439644 : dum = min(dum,1._r8)
3038 20439644 : dum = max(dum,0._r8)
3039 20439644 : pgracs(i) = (1._r8-dum)*pracs(i)
3040 20439644 : ngracs(i) = (1._r8-dum)*npracs(i)
3041 : ! limit max number converted to min of either rain or snow number concentration
3042 20439644 : ngracs(i) = min(ngracs(i),nric(i)/dtime)
3043 20439644 : ngracs(i) = min(ngracs(i),nsic(i)/dtime)
3044 :
3045 : ! amount left for snow production
3046 20439644 : pracs(i) = pracs(i) - pgracs(i)
3047 20439644 : npracs(i) = npracs(i) - ngracs(i)
3048 :
3049 : ! conversion to graupel due to collection of snow by rain
3050 20439644 : psacr(i)=psacr(i)*(1._r8-dum)
3051 : else
3052 5870463556 : pgracs(i) = 0._r8
3053 5870463556 : 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 36871539 : if (qcic(i).ge.0.5e-3_r8.or.qric(i).ge.0.1e-3_r8) then
3103 24635142 : if (psacwg(i).gt.0._r8.or.pracg(i).gt.0._r8) then
3104 24635142 : if (t(i).lt.tm_3 .and. t(i).gt.tm_8) then
3105 4038875 : if (t(i).gt.tm_3) then
3106 : fmult = 0._r8
3107 4038875 : else if (t(i).le.tm_3.and.t(i).gt.tm_5) then
3108 2143470 : fmult = (tm_3-t(i))/2._r8
3109 1895405 : else if (t(i).ge.tm_8.and.t(i).le.tm_5) then
3110 1895405 : 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 4038875 : if (psacwg(i).gt.0._r8) then
3117 928216 : nmultg(i) = 35.e4_r8*psacwg(i)*fmult*1000._r8
3118 928216 : 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 928216 : qmultg(i) = min(qmultg(i),psacwg(i))
3122 928216 : 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 4038875 : if (pracg(i).gt.0._r8) then
3130 3996849 : nmultrg(i) = 35.e4_r8*pracg(i)*fmult*1000._r8
3131 3996849 : 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 3996849 : qmultrg(i) = min(qmultrg(i),pracg(i))
3135 3996849 : 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 1719003247 : pure function limiter_is_on(lim)
3437 : !$acc routine seq
3438 :
3439 : real(r8), intent(in) :: lim
3440 : logical :: limiter_is_on
3441 :
3442 1719003247 : limiter_is_on = transfer(lim, limiter_off) /= limiter_off
3443 :
3444 1719003247 : end function limiter_is_on
3445 :
3446 0 : end module micro_pumas_utils
|
