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