Line data Source code
1 : !===============================================================================
2 : ! Modal Aerosol Model
3 : !===============================================================================
4 : module aero_model
5 : use shr_kind_mod, only: r8 => shr_kind_r8
6 : use constituents, only: pcnst, cnst_name, cnst_get_ind
7 : use ppgrid, only: pcols, pver, pverp
8 : use phys_control, only: phys_getopts, cam_physpkg_is
9 : use cam_abortutils, only: endrun
10 : use cam_logfile, only: iulog
11 : use perf_mod, only: t_startf, t_stopf
12 : use camsrfexch, only: cam_in_t, cam_out_t
13 : use aerodep_flx, only: aerodep_flx_prescribed
14 : use physics_types, only: physics_state, physics_ptend, physics_ptend_init
15 : use physics_buffer, only: physics_buffer_desc
16 : use physics_buffer, only: pbuf_get_field, pbuf_get_index, pbuf_set_field
17 : use physconst, only: gravit, rair, rhoh2o
18 : use spmd_utils, only: masterproc
19 : use infnan, only: nan, assignment(=)
20 :
21 : use cam_history, only: outfld, fieldname_len
22 : use chem_mods, only: gas_pcnst, adv_mass
23 : use mo_tracname, only: solsym
24 :
25 : use modal_aero_data,only: cnst_name_cw, lptr_so4_cw_amode
26 : use modal_aero_data,only: ntot_amode, modename_amode, nspec_max
27 :
28 : use ref_pres, only: top_lev => clim_modal_aero_top_lev
29 :
30 : use modal_aero_wateruptake, only: modal_strat_sulfate
31 : use mo_setsox, only: setsox, has_sox
32 :
33 : implicit none
34 : private
35 :
36 : public :: aero_model_readnl
37 : public :: aero_model_register
38 : public :: aero_model_init
39 : public :: aero_model_gasaerexch ! create, grow, change, and shrink aerosols.
40 : public :: aero_model_drydep ! aerosol dry deposition and sediment
41 : public :: aero_model_wetdep ! aerosol wet removal
42 : public :: aero_model_emissions ! aerosol emissions
43 : public :: aero_model_surfarea ! tropopspheric aerosol wet surface area for chemistry
44 : public :: aero_model_strat_surfarea ! stratospheric aerosol wet surface area for chemistry
45 :
46 : public :: calc_1_impact_rate
47 : public :: nimptblgrow_mind, nimptblgrow_maxd
48 :
49 : ! Accessor functions
50 : public :: get_scavimptblvol, get_scavimptblnum, get_dlndg_nimptblgrow
51 :
52 : ! Misc private data
53 :
54 : ! number of modes
55 : integer :: nmodes
56 : integer :: pblh_idx = 0
57 : integer :: dgnum_idx = 0
58 : integer :: dgnumwet_idx = 0
59 : integer :: rate1_cw2pr_st_idx = 0
60 :
61 : integer :: wetdens_ap_idx = 0
62 : integer :: qaerwat_idx = 0
63 :
64 : integer :: fracis_idx = 0
65 : integer :: prain_idx = 0
66 : integer :: rprddp_idx = 0
67 : integer :: rprdsh_idx = 0
68 : integer :: nevapr_shcu_idx = 0
69 : integer :: nevapr_dpcu_idx = 0
70 :
71 : integer :: sulfeq_idx = -1
72 :
73 : integer :: nh3_ndx = 0
74 : integer :: nh4_ndx = 0
75 :
76 : ! variables for table lookup of aerosol impaction/interception scavenging rates
77 : integer, parameter :: nimptblgrow_mind=-7, nimptblgrow_maxd=12
78 : real(r8) :: dlndg_nimptblgrow
79 : real(r8),allocatable :: scavimptblnum(:,:)
80 : real(r8),allocatable :: scavimptblvol(:,:)
81 :
82 : ! for surf_area_dens
83 : integer,allocatable :: num_idx(:)
84 : integer,allocatable :: index_tot_mass(:,:)
85 : integer,allocatable :: index_chm_mass(:,:)
86 :
87 : integer :: ndx_h2so4
88 : character(len=fieldname_len), allocatable :: dgnum_name(:), dgnumwet_name(:)
89 :
90 : ! Namelist variables
91 : character(len=16) :: wetdep_list(pcnst) = ' '
92 : character(len=16) :: drydep_list(pcnst) = ' '
93 : real(r8) :: sol_facti_cloud_borne = 1._r8
94 : real(r8) :: sol_factb_interstitial = 0.1_r8
95 : real(r8) :: sol_factic_interstitial = 0.4_r8
96 : real(r8) :: seasalt_emis_scale
97 :
98 : integer :: ndrydep = 0
99 : integer,allocatable :: drydep_indices(:)
100 : integer :: nwetdep = 0
101 : integer,allocatable :: wetdep_indices(:)
102 : logical :: drydep_lq(pcnst)
103 : logical :: wetdep_lq(pcnst)
104 :
105 : logical :: modal_accum_coarse_exch = .false.
106 :
107 : logical :: convproc_do_aer
108 :
109 : contains
110 :
111 : !=============================================================================
112 : ! reads aerosol namelist options
113 : !=============================================================================
114 1536 : subroutine aero_model_readnl(nlfile)
115 :
116 : use namelist_utils, only: find_group_name
117 : use units, only: getunit, freeunit
118 : use mpishorthand
119 : use modal_aero_convproc, only: ma_convproc_readnl
120 :
121 : character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input
122 :
123 : ! Local variables
124 : integer :: unitn, ierr
125 : character(len=*), parameter :: subname = 'aero_model_readnl'
126 :
127 : ! Namelist variables
128 : character(len=16) :: aer_wetdep_list(pcnst) = ' '
129 : character(len=16) :: aer_drydep_list(pcnst) = ' '
130 :
131 : namelist /aerosol_nl/ aer_wetdep_list, aer_drydep_list, sol_facti_cloud_borne, &
132 : sol_factb_interstitial, sol_factic_interstitial, modal_strat_sulfate, modal_accum_coarse_exch, seasalt_emis_scale
133 :
134 : !-----------------------------------------------------------------------------
135 :
136 : ! Read namelist
137 1536 : if (masterproc) then
138 2 : unitn = getunit()
139 2 : open( unitn, file=trim(nlfile), status='old' )
140 2 : call find_group_name(unitn, 'aerosol_nl', status=ierr)
141 2 : if (ierr == 0) then
142 2 : read(unitn, aerosol_nl, iostat=ierr)
143 2 : if (ierr /= 0) then
144 0 : call endrun(subname // ':: ERROR reading namelist')
145 : end if
146 : end if
147 2 : close(unitn)
148 2 : call freeunit(unitn)
149 : end if
150 :
151 : #ifdef SPMD
152 : ! Broadcast namelist variables
153 1536 : call mpibcast(aer_wetdep_list, len(aer_wetdep_list(1))*pcnst, mpichar, 0, mpicom)
154 1536 : call mpibcast(aer_drydep_list, len(aer_drydep_list(1))*pcnst, mpichar, 0, mpicom)
155 1536 : call mpibcast(sol_facti_cloud_borne, 1, mpir8, 0, mpicom)
156 1536 : call mpibcast(sol_factb_interstitial, 1, mpir8, 0, mpicom)
157 1536 : call mpibcast(sol_factic_interstitial, 1, mpir8, 0, mpicom)
158 1536 : call mpibcast(modal_strat_sulfate, 1, mpilog, 0, mpicom)
159 1536 : call mpibcast(seasalt_emis_scale, 1, mpir8, 0, mpicom)
160 1536 : call mpibcast(modal_accum_coarse_exch, 1, mpilog, 0, mpicom)
161 : #endif
162 :
163 64512 : wetdep_list = aer_wetdep_list
164 64512 : drydep_list = aer_drydep_list
165 :
166 1536 : call ma_convproc_readnl(nlfile)
167 :
168 1536 : end subroutine aero_model_readnl
169 :
170 : !=============================================================================
171 : !=============================================================================
172 1536 : subroutine aero_model_register()
173 1536 : use modal_aero_data, only: modal_aero_data_reg
174 :
175 1536 : call modal_aero_data_reg()
176 :
177 1536 : end subroutine aero_model_register
178 :
179 : !=============================================================================
180 : !=============================================================================
181 1536 : subroutine aero_model_init( pbuf2d )
182 :
183 : use mo_chem_utls, only: get_inv_ndx
184 : use cam_history, only: addfld, add_default, horiz_only
185 1536 : use mo_chem_utls, only: get_spc_ndx
186 : use modal_aero_data, only: cnst_name_cw
187 : use modal_aero_data, only: modal_aero_data_init
188 : use rad_constituents,only: rad_cnst_get_info
189 : use dust_model, only: dust_init, dust_names, dust_active, dust_nbin, dust_nnum
190 : use seasalt_model, only: seasalt_init, seasalt_names, seasalt_active,seasalt_nbin
191 : use aer_drydep_mod, only: inidrydep
192 : use wetdep, only: wetdep_init
193 :
194 : use modal_aero_calcsize, only: modal_aero_calcsize_init
195 : use modal_aero_coag, only: modal_aero_coag_init
196 : use modal_aero_deposition, only: modal_aero_deposition_init
197 : use modal_aero_gasaerexch, only: modal_aero_gasaerexch_init
198 : use modal_aero_newnuc, only: modal_aero_newnuc_init
199 : use modal_aero_rename, only: modal_aero_rename_init
200 : use modal_aero_convproc, only: ma_convproc_init
201 :
202 : ! args
203 : type(physics_buffer_desc), pointer :: pbuf2d(:,:)
204 :
205 : ! local vars
206 : character(len=*), parameter :: subrname = 'aero_model_init'
207 : integer :: m, n, id
208 : character(len=20) :: dummy
209 :
210 : logical :: history_aerosol ! Output MAM or SECT aerosol tendencies
211 : logical :: history_chemistry, history_cesm_forcing, history_dust
212 :
213 : integer :: l
214 : character(len=6) :: test_name
215 : character(len=64) :: errmes
216 :
217 : character(len=2) :: unit_basename ! Units 'kg' or '1'
218 : integer :: errcode
219 : character(len=fieldname_len) :: field_name
220 :
221 : character(len=32) :: spec_name
222 : character(len=32) :: spec_type
223 : character(len=32) :: mode_type
224 : integer :: nspec
225 :
226 1536 : dgnum_idx = pbuf_get_index('DGNUM')
227 1536 : dgnumwet_idx = pbuf_get_index('DGNUMWET')
228 1536 : fracis_idx = pbuf_get_index('FRACIS')
229 1536 : prain_idx = pbuf_get_index('PRAIN')
230 1536 : rprddp_idx = pbuf_get_index('RPRDDP')
231 1536 : rprdsh_idx = pbuf_get_index('RPRDSH')
232 1536 : nevapr_shcu_idx = pbuf_get_index('NEVAPR_SHCU')
233 1536 : nevapr_dpcu_idx = pbuf_get_index('NEVAPR_DPCU')
234 1536 : sulfeq_idx = pbuf_get_index('MAMH2SO4EQ',errcode)
235 :
236 : call phys_getopts(history_aerosol_out = history_aerosol, &
237 : history_chemistry_out=history_chemistry, &
238 : history_cesm_forcing_out=history_cesm_forcing, &
239 : history_dust_out=history_dust, &
240 1536 : convproc_do_aer_out = convproc_do_aer)
241 :
242 1536 : call rad_cnst_get_info(0, nmodes=nmodes)
243 :
244 1536 : call modal_aero_data_init(pbuf2d)
245 1536 : call modal_aero_bcscavcoef_init()
246 :
247 1536 : call modal_aero_rename_init( modal_accum_coarse_exch )
248 : ! calcsize call must follow rename call
249 1536 : call modal_aero_calcsize_init( pbuf2d )
250 1536 : call modal_aero_gasaerexch_init
251 : ! coag call must follow gasaerexch call
252 1536 : call modal_aero_coag_init
253 1536 : call modal_aero_newnuc_init
254 :
255 : ! call modal_aero_deposition_init only if the user has not specified
256 : ! prescribed aerosol deposition fluxes
257 1536 : if (.not.aerodep_flx_prescribed()) then
258 1536 : call modal_aero_deposition_init
259 : endif
260 :
261 1536 : if (convproc_do_aer) then
262 1536 : call ma_convproc_init()
263 : endif
264 :
265 1536 : call dust_init()
266 1536 : call seasalt_init(seasalt_emis_scale)
267 1536 : call wetdep_init()
268 :
269 1536 : nwetdep = 0
270 1536 : ndrydep = 0
271 :
272 64512 : count_species: do m = 1,pcnst
273 62976 : if ( len_trim(wetdep_list(m)) /= 0 ) then
274 29184 : nwetdep = nwetdep+1
275 : endif
276 64512 : if ( len_trim(drydep_list(m)) /= 0 ) then
277 29184 : ndrydep = ndrydep+1
278 : endif
279 : enddo count_species
280 :
281 1536 : if (nwetdep>0) &
282 4608 : allocate(wetdep_indices(nwetdep))
283 1536 : if (ndrydep>0) &
284 4608 : allocate(drydep_indices(ndrydep))
285 :
286 30720 : do m = 1,ndrydep
287 29184 : call cnst_get_ind ( drydep_list(m), id, abort=.false. )
288 29184 : if (id>0) then
289 29184 : drydep_indices(m) = id
290 : else
291 0 : call endrun(subrname//': invalid drydep species: '//trim(drydep_list(m)) )
292 : endif
293 :
294 30720 : if (masterproc) then
295 38 : write(iulog,*) subrname//': '//drydep_list(m)//' will have drydep applied'
296 : endif
297 : enddo
298 30720 : do m = 1,nwetdep
299 29184 : call cnst_get_ind ( wetdep_list(m), id, abort=.false. )
300 29184 : if (id>0) then
301 29184 : wetdep_indices(m) = id
302 : else
303 0 : call endrun(subrname//': invalid wetdep species: '//trim(wetdep_list(m)) )
304 : endif
305 :
306 30720 : if (masterproc) then
307 38 : write(iulog,*) subrname//': '//wetdep_list(m)//' will have wet removal'
308 : endif
309 : enddo
310 :
311 1536 : if (ndrydep>0) then
312 :
313 1536 : call inidrydep(rair, gravit)
314 :
315 1536 : dummy = 'RAM1'
316 1536 : call addfld (dummy,horiz_only, 'A','frac','RAM1')
317 1536 : if ( history_aerosol ) then
318 0 : call add_default (dummy, 1, ' ')
319 : endif
320 1536 : dummy = 'airFV'
321 1536 : call addfld (dummy,horiz_only, 'A','frac','FV')
322 1536 : if ( history_aerosol ) then
323 0 : call add_default (dummy, 1, ' ')
324 : endif
325 :
326 : endif
327 :
328 1536 : if (dust_active) then
329 : ! emissions diagnostics ....
330 :
331 10752 : do m = 1, dust_nbin+dust_nnum
332 9216 : dummy = trim(dust_names(m)) // 'SF'
333 9216 : call addfld (dummy,horiz_only, 'A','kg/m2/s',trim(dust_names(m))//' dust surface emission')
334 10752 : if (history_aerosol.or.history_chemistry) then
335 9216 : call add_default (dummy, 1, ' ')
336 : endif
337 : enddo
338 :
339 1536 : dummy = 'DSTSFMBL'
340 1536 : call addfld (dummy,horiz_only, 'A','kg/m2/s','Mobilization flux at surface')
341 1536 : if (history_aerosol .or. history_dust) then
342 0 : call add_default (dummy, 1, ' ')
343 : endif
344 :
345 1536 : dummy = 'LND_MBL'
346 1536 : call addfld (dummy,horiz_only, 'A','frac','Soil erodibility factor')
347 1536 : if (history_aerosol) then
348 0 : call add_default (dummy, 1, ' ')
349 : endif
350 :
351 : endif
352 :
353 1536 : if (seasalt_active) then
354 :
355 1536 : dummy = 'SSTSFMBL'
356 1536 : call addfld (dummy,horiz_only, 'A','kg/m2/s','Mobilization flux at surface')
357 1536 : if (history_aerosol) then
358 0 : call add_default (dummy, 1, ' ')
359 : endif
360 :
361 6144 : do m = 1, seasalt_nbin
362 4608 : dummy = trim(seasalt_names(m)) // 'SF'
363 4608 : call addfld (dummy,horiz_only, 'A','kg/m2/s',trim(seasalt_names(m))//' seasalt surface emission')
364 6144 : if (history_aerosol.or.history_chemistry) then
365 4608 : call add_default (dummy, 1, ' ')
366 : endif
367 : enddo
368 :
369 : endif
370 :
371 :
372 : ! set flags for drydep tendencies
373 1536 : drydep_lq(:) = .false.
374 30720 : do m=1,ndrydep
375 29184 : id = drydep_indices(m)
376 30720 : drydep_lq(id) = .true.
377 : enddo
378 :
379 : ! set flags for wetdep tendencies
380 1536 : wetdep_lq(:) = .false.
381 30720 : do m=1,nwetdep
382 29184 : id = wetdep_indices(m)
383 30720 : wetdep_lq(id) = .true.
384 : enddo
385 :
386 1536 : wetdens_ap_idx = pbuf_get_index('WETDENS_AP')
387 1536 : qaerwat_idx = pbuf_get_index('QAERWAT')
388 1536 : pblh_idx = pbuf_get_index('pblh')
389 :
390 1536 : rate1_cw2pr_st_idx = pbuf_get_index('RATE1_CW2PR_ST')
391 1536 : call pbuf_set_field(pbuf2d, rate1_cw2pr_st_idx, 0.0_r8)
392 :
393 30720 : do m = 1,ndrydep
394 :
395 : ! units
396 29184 : if (drydep_list(m)(1:3) == 'num') then
397 6144 : unit_basename = ' 1'
398 : else
399 23040 : unit_basename = 'kg'
400 : endif
401 :
402 : call addfld (trim(drydep_list(m))//'DDF', horiz_only, 'A',unit_basename//'/m2/s ', &
403 29184 : trim(drydep_list(m))//' dry deposition flux at bottom (grav + turb)')
404 : call addfld (trim(drydep_list(m))//'TBF', horiz_only, 'A',unit_basename//'/m2/s', &
405 29184 : trim(drydep_list(m))//' turbulent dry deposition flux')
406 : call addfld (trim(drydep_list(m))//'GVF', horiz_only, 'A',unit_basename//'/m2/s ', &
407 29184 : trim(drydep_list(m))//' gravitational dry deposition flux')
408 : call addfld (trim(drydep_list(m))//'DTQ', (/ 'lev' /), 'A',unit_basename//'/kg/s ', &
409 58368 : trim(drydep_list(m))//' dry deposition')
410 : call addfld (trim(drydep_list(m))//'DDV', (/ 'lev' /), 'A','m/s', &
411 58368 : trim(drydep_list(m))//' deposition velocity')
412 :
413 29184 : if ( history_aerosol.or.history_chemistry ) then
414 29184 : call add_default (trim(drydep_list(m))//'DDF', 1, ' ')
415 : endif
416 30720 : if ( history_aerosol ) then
417 0 : call add_default (trim(drydep_list(m))//'TBF', 1, ' ')
418 0 : call add_default (trim(drydep_list(m))//'GVF', 1, ' ')
419 : endif
420 :
421 : enddo
422 :
423 30720 : do m = 1,nwetdep
424 :
425 : ! units
426 29184 : if (wetdep_list(m)(1:3) == 'num') then
427 6144 : unit_basename = ' 1'
428 : else
429 23040 : unit_basename = 'kg'
430 : endif
431 :
432 : call addfld (trim(wetdep_list(m))//'SFWET', &
433 29184 : horiz_only, 'A',unit_basename//'/m2/s ','Wet deposition flux at surface')
434 : call addfld (trim(wetdep_list(m))//'SFSIC', &
435 29184 : horiz_only, 'A',unit_basename//'/m2/s ','Wet deposition flux (incloud, convective) at surface')
436 : call addfld (trim(wetdep_list(m))//'SFSIS', &
437 29184 : horiz_only, 'A',unit_basename//'/m2/s ','Wet deposition flux (incloud, stratiform) at surface')
438 : call addfld (trim(wetdep_list(m))//'SFSBC', &
439 29184 : horiz_only, 'A',unit_basename//'/m2/s ','Wet deposition flux (belowcloud, convective) at surface')
440 : call addfld (trim(wetdep_list(m))//'SFSBS', &
441 29184 : horiz_only, 'A',unit_basename//'/m2/s ','Wet deposition flux (belowcloud, stratiform) at surface')
442 :
443 29184 : if (convproc_do_aer) then
444 : call addfld (trim(wetdep_list(m))//'SFSES', &
445 29184 : horiz_only, 'A','kg/m2/s','Wet deposition flux (precip evap, stratiform) at surface')
446 : call addfld (trim(wetdep_list(m))//'SFSBD', &
447 29184 : horiz_only, 'A','kg/m2/s','Wet deposition flux (belowcloud, deep convective) at surface')
448 : end if
449 :
450 58368 : call addfld (trim(wetdep_list(m))//'WET',(/ 'lev' /), 'A',unit_basename//'/kg/s ','wet deposition tendency')
451 : call addfld (trim(wetdep_list(m))//'SIC',(/ 'lev' /), 'A',unit_basename//'/kg/s ', &
452 58368 : trim(wetdep_list(m))//' ic wet deposition')
453 : call addfld (trim(wetdep_list(m))//'SIS',(/ 'lev' /), 'A',unit_basename//'/kg/s ', &
454 58368 : trim(wetdep_list(m))//' is wet deposition')
455 : call addfld (trim(wetdep_list(m))//'SBC',(/ 'lev' /), 'A',unit_basename//'/kg/s ', &
456 58368 : trim(wetdep_list(m))//' bc wet deposition')
457 : call addfld (trim(wetdep_list(m))//'SBS',(/ 'lev' /), 'A',unit_basename//'/kg/s ', &
458 58368 : trim(wetdep_list(m))//' bs wet deposition')
459 :
460 29184 : if ( history_aerosol .or. history_chemistry ) then
461 29184 : call add_default (trim(wetdep_list(m))//'SFWET', 1, ' ')
462 : endif
463 30720 : if ( history_aerosol ) then
464 0 : call add_default (trim(wetdep_list(m))//'SFSIC', 1, ' ')
465 0 : call add_default (trim(wetdep_list(m))//'SFSIS', 1, ' ')
466 0 : call add_default (trim(wetdep_list(m))//'SFSBC', 1, ' ')
467 0 : call add_default (trim(wetdep_list(m))//'SFSBS', 1, ' ')
468 0 : if (convproc_do_aer) then
469 0 : call add_default (trim(wetdep_list(m))//'SFSES', 1, ' ')
470 0 : call add_default (trim(wetdep_list(m))//'SFSBD', 1, ' ')
471 : end if
472 : endif
473 :
474 : enddo
475 :
476 49152 : do m = 1,gas_pcnst
477 :
478 47616 : if ( solsym(m)(1:3) == 'num') then
479 6144 : unit_basename = ' 1' ! Units 'kg' or '1'
480 : else
481 41472 : unit_basename = 'kg' ! Units 'kg' or '1'
482 : end if
483 :
484 : call addfld( 'GS_'//trim(solsym(m)),horiz_only, 'A', unit_basename//'/m2/s ', &
485 47616 : trim(solsym(m))//' gas chemistry/wet removal (for gas species)')
486 : call addfld( 'AQ_'//trim(solsym(m)),horiz_only, 'A', unit_basename//'/m2/s ', &
487 47616 : trim(solsym(m))//' aqueous chemistry (for gas species)')
488 49152 : if ( history_aerosol ) then
489 0 : call add_default( 'AQ_'//trim(solsym(m)), 1, ' ')
490 : endif
491 :
492 : enddo
493 64512 : do n = 1,pcnst
494 64512 : if( .not. (cnst_name_cw(n) == ' ') ) then
495 :
496 29184 : if (cnst_name_cw(n)(1:3) == 'num') then
497 6144 : unit_basename = ' 1'
498 : else
499 23040 : unit_basename = 'kg'
500 : endif
501 :
502 : call addfld( cnst_name_cw(n), (/ 'lev' /), 'A', unit_basename//'/kg ', &
503 58368 : trim(cnst_name_cw(n))//' in cloud water')
504 29184 : call addfld (trim(cnst_name_cw(n))//'SFWET', horiz_only, 'A', unit_basename//'/m2/s ', &
505 58368 : trim(cnst_name_cw(n))//' wet deposition flux at surface')
506 29184 : call addfld (trim(cnst_name_cw(n))//'SFSIC', horiz_only, 'A', unit_basename//'/m2/s ', &
507 58368 : trim(cnst_name_cw(n))//' wet deposition flux (incloud, convective) at surface')
508 29184 : call addfld (trim(cnst_name_cw(n))//'SFSIS', horiz_only, 'A', unit_basename//'/m2/s ', &
509 58368 : trim(cnst_name_cw(n))//' wet deposition flux (incloud, stratiform) at surface')
510 29184 : call addfld (trim(cnst_name_cw(n))//'SFSBC', horiz_only, 'A', unit_basename//'/m2/s ', &
511 58368 : trim(cnst_name_cw(n))//' wet deposition flux (belowcloud, convective) at surface')
512 29184 : call addfld (trim(cnst_name_cw(n))//'SFSBS', horiz_only, 'A', unit_basename//'/m2/s ', &
513 58368 : trim(cnst_name_cw(n))//' wet deposition flux (belowcloud, stratiform) at surface')
514 29184 : call addfld (trim(cnst_name_cw(n))//'DDF', horiz_only, 'A', unit_basename//'/m2/s ', &
515 58368 : trim(cnst_name_cw(n))//' dry deposition flux at bottom (grav + turb)')
516 29184 : call addfld (trim(cnst_name_cw(n))//'TBF', horiz_only, 'A', unit_basename//'/m2/s ', &
517 58368 : trim(cnst_name_cw(n))//' turbulent dry deposition flux')
518 29184 : call addfld (trim(cnst_name_cw(n))//'GVF', horiz_only, 'A', unit_basename//'/m2/s ', &
519 58368 : trim(cnst_name_cw(n))//' gravitational dry deposition flux')
520 :
521 29184 : if (convproc_do_aer) then
522 29184 : call addfld (trim(cnst_name_cw(n))//'SFSEC', &
523 58368 : horiz_only, 'A','kg/m2/s','Wet deposition flux (precip evap, convective) at surface')
524 29184 : call addfld (trim(cnst_name_cw(n))//'SFSES', &
525 58368 : horiz_only, 'A','kg/m2/s','Wet deposition flux (precip evap, stratiform) at surface')
526 29184 : call addfld (trim(cnst_name_cw(n))//'SFSBD', &
527 58368 : horiz_only, 'A','kg/m2/s','Wet deposition flux (belowcloud, deep convective) at surface')
528 : end if
529 :
530 :
531 29184 : if ( history_aerosol.or. history_chemistry ) then
532 29184 : call add_default( cnst_name_cw(n), 1, ' ' )
533 29184 : call add_default (trim(cnst_name_cw(n))//'SFWET', 1, ' ')
534 : endif
535 29184 : if ( history_aerosol ) then
536 0 : call add_default (trim(cnst_name_cw(n))//'GVF', 1, ' ')
537 0 : call add_default (trim(cnst_name_cw(n))//'TBF', 1, ' ')
538 0 : call add_default (trim(cnst_name_cw(n))//'DDF', 1, ' ')
539 0 : call add_default (trim(cnst_name_cw(n))//'SFSBS', 1, ' ')
540 0 : call add_default (trim(cnst_name_cw(n))//'SFSIC', 1, ' ')
541 0 : call add_default (trim(cnst_name_cw(n))//'SFSBC', 1, ' ')
542 0 : call add_default (trim(cnst_name_cw(n))//'SFSIS', 1, ' ')
543 0 : if (convproc_do_aer) then
544 0 : call add_default (trim(cnst_name_cw(n))//'SFSEC', 1, ' ')
545 0 : call add_default (trim(cnst_name_cw(n))//'SFSES', 1, ' ')
546 0 : call add_default (trim(cnst_name_cw(n))//'SFSBD', 1, ' ')
547 : end if
548 : endif
549 : endif
550 : enddo
551 :
552 6144 : allocate(dgnum_name(ntot_amode), dgnumwet_name(ntot_amode))
553 7680 : do n=1,ntot_amode
554 6144 : dgnum_name(n) = ' '
555 6144 : dgnumwet_name(n) = ' '
556 6144 : write(dgnum_name(n),fmt='(a,i1)') 'dgnum',n
557 6144 : write(dgnumwet_name(n),fmt='(a,i1)') 'dgnumwet',n
558 12288 : call addfld( dgnum_name(n), (/ 'lev' /), 'I', 'm', 'Aerosol mode dry diameter' )
559 12288 : call addfld( dgnumwet_name(n), (/ 'lev' /), 'I', 'm', 'Aerosol mode wet diameter' )
560 6144 : if ( history_aerosol ) then
561 0 : call add_default( dgnum_name(n), 1, ' ' )
562 0 : call add_default( dgnumwet_name(n), 1, ' ' )
563 : endif
564 6144 : if ( history_cesm_forcing .and. n<4 ) then
565 0 : call add_default( dgnumwet_name(n), 8, ' ' )
566 : endif
567 :
568 7680 : if (modal_strat_sulfate) then
569 0 : field_name = ' '
570 0 : write(field_name,fmt='(a,i1)') 'wtpct_a',n
571 0 : call addfld( field_name, (/ 'lev' /), 'I', '%', 'Aerosol mode weight percent H2SO4' )
572 0 : if ( history_aerosol ) then
573 0 : call add_default (field_name, 0, 'I')
574 : endif
575 :
576 0 : field_name = ' '
577 0 : write(field_name,fmt='(a,i1)') 'sulfeq_a',n
578 0 : call addfld( field_name, (/ 'lev' /), 'I', 'kg/kg', 'H2SO4 equilibrium mixing ratio' )
579 0 : if ( history_aerosol ) then
580 0 : call add_default (field_name, 0, 'I')
581 : endif
582 :
583 0 : field_name = ' '
584 0 : write(field_name,fmt='(a,i1)') 'sulden_a',n
585 0 : call addfld( field_name, (/ 'lev' /), 'I', 'g/cm3', 'Sulfate aerosol particle mass density' )
586 0 : if ( history_aerosol ) then
587 0 : call add_default (field_name, 0, 'I')
588 : endif
589 :
590 : end if
591 : end do
592 :
593 1536 : ndx_h2so4 = get_spc_ndx('H2SO4')
594 1536 : nh3_ndx = get_spc_ndx('NH3')
595 1536 : nh4_ndx = get_spc_ndx('NH4')
596 :
597 4608 : allocate(num_idx(ntot_amode))
598 7680 : num_idx = -1
599 :
600 : ! for aero_model_surfarea called from mo_usrrxt
601 7680 : do l=1,ntot_amode
602 6144 : test_name = ' '
603 6144 : write(test_name,fmt='(a5,i1)') 'num_a',l
604 6144 : num_idx(l) = get_spc_ndx( trim(test_name) )
605 7680 : if (num_idx(l) < 0) then
606 0 : write(errmes,fmt='(a,i1)') 'usrrxt_inti: cannot find MAM num_idx ',l
607 0 : write(iulog,*) errmes
608 0 : call endrun(errmes)
609 : endif
610 : end do
611 :
612 6144 : allocate(index_tot_mass(nmodes,nspec_max))
613 4608 : allocate(index_chm_mass(nmodes,nspec_max))
614 47616 : index_tot_mass = -1
615 47616 : index_chm_mass = -1
616 :
617 : ! for surf_area_dens
618 : ! define indices associated with the various aerosol types
619 7680 : do n = 1,nmodes
620 6144 : call rad_cnst_get_info(0, n, mode_type=mode_type, nspec=nspec)
621 7680 : if ( trim(mode_type) /= 'primary_carbon') then ! ignore the primary_carbon mode
622 24576 : do l = 1, nspec
623 19968 : call rad_cnst_get_info(0, n, l, spec_type=spec_type, spec_name=spec_name)
624 19968 : index_tot_mass(n,l) = get_spc_ndx(spec_name)
625 : if ( trim(spec_type) == 'sulfate' .or. &
626 : trim(spec_type) == 's-organic' .or. &
627 : trim(spec_type) == 'p-organic' .or. &
628 19968 : trim(spec_type) == 'black-c' .or. &
629 4608 : trim(spec_type) == 'ammonium') then
630 10752 : index_chm_mass(n,l) = get_spc_ndx(spec_name)
631 : endif
632 : enddo
633 : endif
634 : enddo
635 :
636 1536 : if (has_sox) then
637 7680 : do m = 1, ntot_amode
638 :
639 6144 : l = lptr_so4_cw_amode(m)
640 7680 : if (l > 0) then
641 : call addfld (&
642 4608 : trim(cnst_name_cw(l))//'AQSO4',horiz_only, 'A','kg/m2/s', &
643 9216 : trim(cnst_name_cw(l))//' aqueous phase chemistry')
644 : call addfld (&
645 4608 : trim(cnst_name_cw(l))//'AQH2SO4',horiz_only, 'A','kg/m2/s', &
646 9216 : trim(cnst_name_cw(l))//' aqueous phase chemistry')
647 4608 : if ( history_aerosol ) then
648 0 : call add_default (trim(cnst_name_cw(l))//'AQSO4', 1, ' ')
649 0 : call add_default (trim(cnst_name_cw(l))//'AQH2SO4', 1, ' ')
650 : endif
651 : end if
652 :
653 : end do
654 :
655 3072 : call addfld( 'XPH_LWC', (/ 'lev' /), 'A','kg/kg', 'pH value multiplied by lwc')
656 1536 : call addfld ('AQSO4_H2O2', horiz_only, 'A','kg/m2/s', 'SO4 aqueous phase chemistry due to H2O2')
657 1536 : call addfld ('AQSO4_O3', horiz_only, 'A','kg/m2/s', 'SO4 aqueous phase chemistry due to O3')
658 :
659 1536 : if ( history_aerosol ) then
660 0 : call add_default ('XPH_LWC', 1, ' ')
661 0 : call add_default ('AQSO4_H2O2', 1, ' ')
662 0 : call add_default ('AQSO4_O3', 1, ' ')
663 : endif
664 : endif
665 :
666 3072 : end subroutine aero_model_init
667 :
668 : !=============================================================================
669 : !=============================================================================
670 2529432 : subroutine aero_model_drydep ( state, pbuf, obklen, ustar, cam_in, dt, cam_out, ptend )
671 :
672 1536 : use dust_sediment_mod, only: dust_sediment_tend
673 : use aer_drydep_mod, only: d3ddflux, calcram
674 : use modal_aero_data, only: qqcw_get_field
675 : use modal_aero_data, only: cnst_name_cw
676 : use modal_aero_data, only: alnsg_amode
677 : use modal_aero_data, only: sigmag_amode
678 : use modal_aero_data, only: nspec_amode
679 : use modal_aero_data, only: numptr_amode
680 : use modal_aero_data, only: numptrcw_amode
681 : use modal_aero_data, only: lmassptr_amode
682 : use modal_aero_data, only: lmassptrcw_amode
683 : use modal_aero_deposition, only: set_srf_drydep
684 :
685 : ! args
686 : type(physics_state), intent(in) :: state ! Physics state variables
687 : real(r8), intent(in) :: obklen(:)
688 : real(r8), intent(in) :: ustar(:) ! sfc fric vel
689 : type(cam_in_t), target, intent(in) :: cam_in ! import state
690 : real(r8), intent(in) :: dt ! time step
691 : type(cam_out_t), intent(inout) :: cam_out ! export state
692 : type(physics_ptend), intent(out) :: ptend ! indivdual parameterization tendencies
693 : type(physics_buffer_desc), pointer :: pbuf(:)
694 :
695 : ! local vars
696 58824 : real(r8), pointer :: landfrac(:) ! land fraction
697 58824 : real(r8), pointer :: icefrac(:) ! ice fraction
698 58824 : real(r8), pointer :: ocnfrac(:) ! ocean fraction
699 58824 : real(r8), pointer :: fvin(:) !
700 58824 : real(r8), pointer :: ram1in(:) ! for dry dep velocities from land model for progseasalts
701 :
702 : real(r8) :: fv(pcols) ! for dry dep velocities, from land modified over ocean & ice
703 : real(r8) :: ram1(pcols) ! for dry dep velocities, from land modified over ocean & ice
704 :
705 : integer :: lchnk ! chunk identifier
706 : integer :: ncol ! number of atmospheric columns
707 : integer :: jvlc ! index for last dimension of vlc_xxx arrays
708 : integer :: lphase ! index for interstitial / cloudborne aerosol
709 : integer :: lspec ! index for aerosol number / chem-mass / water-mass
710 : integer :: m ! aerosol mode index
711 : integer :: mm ! tracer index
712 : integer :: i
713 :
714 : real(r8) :: tvs(pcols,pver)
715 : real(r8) :: rho(pcols,pver) ! air density in kg/m3
716 : real(r8) :: sflx(pcols) ! deposition flux
717 : real(r8) :: dep_trb(pcols) !kg/m2/s
718 : real(r8) :: dep_grv(pcols) !kg/m2/s (total of grav and trb)
719 : real(r8) :: pvmzaer(pcols,pverp) ! sedimentation velocity in Pa
720 : real(r8) :: dqdt_tmp(pcols,pver) ! temporary array to hold tendency for 1 species
721 :
722 : real(r8) :: rad_drop(pcols,pver)
723 : real(r8) :: dens_drop(pcols,pver)
724 : real(r8) :: sg_drop(pcols,pver)
725 : real(r8) :: rad_aer(pcols,pver)
726 : real(r8) :: dens_aer(pcols,pver)
727 : real(r8) :: sg_aer(pcols,pver)
728 :
729 : real(r8) :: vlc_dry(pcols,pver,4) ! dep velocity
730 : real(r8) :: vlc_grv(pcols,pver,4) ! dep velocity
731 : real(r8):: vlc_trb(pcols,4) ! dep velocity
732 : real(r8) :: aerdepdryis(pcols,pcnst) ! aerosol dry deposition (interstitial)
733 : real(r8) :: aerdepdrycw(pcols,pcnst) ! aerosol dry deposition (cloud water)
734 58824 : real(r8), pointer :: fldcw(:,:)
735 58824 : real(r8), pointer :: dgncur_awet(:,:,:)
736 58824 : real(r8), pointer :: wetdens(:,:,:)
737 58824 : real(r8), pointer :: qaerwat(:,:,:)
738 :
739 58824 : landfrac => cam_in%landfrac(:)
740 58824 : icefrac => cam_in%icefrac(:)
741 58824 : ocnfrac => cam_in%ocnfrac(:)
742 58824 : fvin => cam_in%fv(:)
743 58824 : ram1in => cam_in%ram1(:)
744 :
745 58824 : lchnk = state%lchnk
746 58824 : ncol = state%ncol
747 :
748 : ! calc ram and fv over ocean and sea ice ...
749 : call calcram( ncol,landfrac,icefrac,ocnfrac,obklen,&
750 : ustar,ram1in,ram1,state%t(:,pver),state%pmid(:,pver),&
751 58824 : state%pdel(:,pver),fvin,fv)
752 :
753 58824 : call outfld( 'airFV', fv(:), pcols, lchnk )
754 58824 : call outfld( 'RAM1', ram1(:), pcols, lchnk )
755 :
756 : ! note that tendencies are not only in sfc layer (because of sedimentation)
757 : ! and that ptend is updated within each subroutine for different species
758 :
759 58824 : call physics_ptend_init(ptend, state%psetcols, 'aero_model_drydep', lq=drydep_lq)
760 :
761 235296 : call pbuf_get_field(pbuf, dgnumwet_idx, dgncur_awet, start=(/1,1,1/), kount=(/pcols,pver,nmodes/) )
762 235296 : call pbuf_get_field(pbuf, wetdens_ap_idx, wetdens, start=(/1,1,1/), kount=(/pcols,pver,nmodes/) )
763 235296 : call pbuf_get_field(pbuf, qaerwat_idx, qaerwat, start=(/1,1,1/), kount=(/pcols,pver,nmodes/) )
764 :
765 91405656 : tvs(:ncol,:) = state%t(:ncol,:)!*(1+state%q(:ncol,k)
766 91405656 : rho(:ncol,:)= state%pmid(:ncol,:)/(rair*state%t(:ncol,:))
767 :
768 : !
769 : ! calc settling/deposition velocities for cloud droplets (and cloud-borne aerosols)
770 : !
771 : ! *** mean drop radius should eventually be computed from ndrop and qcldwtr
772 93059568 : rad_drop(:,:) = 5.0e-6_r8
773 93059568 : dens_drop(:,:) = rhoh2o
774 93059568 : sg_drop(:,:) = 1.46_r8
775 58824 : jvlc = 3
776 : call modal_aero_depvel_part( ncol,state%t(:,:), state%pmid(:,:), ram1, fv, &
777 : vlc_dry(:,:,jvlc), vlc_trb(:,jvlc), vlc_grv(:,:,jvlc), &
778 58824 : rad_drop(:,:), dens_drop(:,:), sg_drop(:,:), 0, lchnk)
779 58824 : jvlc = 4
780 : call modal_aero_depvel_part( ncol,state%t(:,:), state%pmid(:,:), ram1, fv, &
781 : vlc_dry(:,:,jvlc), vlc_trb(:,jvlc), vlc_grv(:,:,jvlc), &
782 58824 : rad_drop(:,:), dens_drop(:,:), sg_drop(:,:), 3, lchnk)
783 :
784 :
785 :
786 294120 : do m = 1, ntot_amode ! main loop over aerosol modes
787 :
788 764712 : do lphase = 1, 2 ! loop over interstitial / cloud-borne forms
789 :
790 470592 : if (lphase == 1) then ! interstial aerosol - calc settling/dep velocities of mode
791 :
792 : ! rad_aer = volume mean wet radius (m)
793 : ! dgncur_awet = geometric mean wet diameter for number distribution (m)
794 470592 : rad_aer(1:ncol,:) = 0.5_r8*dgncur_awet(1:ncol,:,m) &
795 366093216 : *exp(1.5_r8*(alnsg_amode(m)**2))
796 : ! dens_aer(1:ncol,:) = wet density (kg/m3)
797 365622624 : dens_aer(1:ncol,:) = wetdens(1:ncol,:,m)
798 365622624 : sg_aer(1:ncol,:) = sigmag_amode(m)
799 :
800 235296 : jvlc = 1
801 : call modal_aero_depvel_part( ncol, state%t(:,:), state%pmid(:,:), ram1, fv, &
802 : vlc_dry(:,:,jvlc), vlc_trb(:,jvlc), vlc_grv(:,:,jvlc), &
803 235296 : rad_aer(:,:), dens_aer(:,:), sg_aer(:,:), 0, lchnk)
804 235296 : jvlc = 2
805 : call modal_aero_depvel_part( ncol, state%t(:,:), state%pmid(:,:), ram1, fv, &
806 : vlc_dry(:,:,jvlc), vlc_trb(:,jvlc), vlc_grv(:,:,jvlc), &
807 235296 : rad_aer(:,:), dens_aer(:,:), sg_aer(:,:), 3, lchnk)
808 : end if
809 :
810 3411792 : do lspec = 0, nspec_amode(m)+1 ! loop over number + constituents + water
811 :
812 2705904 : if (lspec == 0) then ! number
813 470592 : if (lphase == 1) then
814 235296 : mm = numptr_amode(m)
815 235296 : jvlc = 1
816 : else
817 235296 : mm = numptrcw_amode(m)
818 235296 : jvlc = 3
819 : endif
820 2235312 : else if (lspec <= nspec_amode(m)) then ! non-water mass
821 1764720 : if (lphase == 1) then
822 882360 : mm = lmassptr_amode(lspec,m)
823 882360 : jvlc = 2
824 : else
825 882360 : mm = lmassptrcw_amode(lspec,m)
826 882360 : jvlc = 4
827 : endif
828 : else ! water mass
829 : ! bypass dry deposition of aerosol water
830 : cycle
831 : if (lphase == 1) then
832 : mm = 0
833 : ! mm = lwaterptr_amode(m)
834 : jvlc = 2
835 : else
836 : mm = 0
837 : jvlc = 4
838 : endif
839 : endif
840 :
841 :
842 2235312 : if (mm <= 0) cycle
843 :
844 : ! if (lphase == 1) then
845 2705904 : if ((lphase == 1) .and. (lspec <= nspec_amode(m))) then
846 1117656 : ptend%lq(mm) = .TRUE.
847 :
848 : ! use pvprogseasalts instead (means making the top level 0)
849 18662256 : pvmzaer(:ncol,1)=0._r8
850 1736707464 : pvmzaer(:ncol,2:pverp) = vlc_dry(:ncol,:,jvlc)
851 :
852 1117656 : call outfld( trim(cnst_name(mm))//'DDV', pvmzaer(:,2:pverp), pcols, lchnk )
853 :
854 : if(.true.) then ! use phil's method
855 : ! convert from meters/sec to pascals/sec
856 : ! pvprogseasalts(:,1) is assumed zero, use density from layer above in conversion
857 1736707464 : pvmzaer(:ncol,2:pverp) = pvmzaer(:ncol,2:pverp) * rho(:ncol,:)*gravit
858 :
859 : ! calculate the tendencies and sfc fluxes from the above velocities
860 : call dust_sediment_tend( &
861 : ncol, dt, state%pint(:,:), state%pmid, state%pdel, state%t , &
862 1117656 : state%q(:,:,mm), pvmzaer, ptend%q(:,:,mm), sflx )
863 : else !use charlie's method
864 : call d3ddflux( ncol, vlc_dry(:,:,jvlc), state%q(:,:,mm), state%pmid, &
865 : state%pdel, tvs, sflx, ptend%q(:,:,mm), dt )
866 : endif
867 :
868 : ! apportion dry deposition into turb and gravitational settling for tapes
869 1117656 : dep_trb = 0._r8
870 1117656 : dep_grv = 0._r8
871 18662256 : do i=1,ncol
872 18662256 : if (vlc_dry(i,pver,jvlc) /= 0._r8) then
873 17544600 : dep_trb(i)=sflx(i)*vlc_trb(i,jvlc)/vlc_dry(i,pver,jvlc)
874 17544600 : dep_grv(i)=sflx(i)*vlc_grv(i,pver,jvlc)/vlc_dry(i,pver,jvlc)
875 : end if
876 : enddo
877 :
878 1117656 : call outfld( trim(cnst_name(mm))//'DDF', sflx, pcols, lchnk)
879 1117656 : call outfld( trim(cnst_name(mm))//'TBF', dep_trb, pcols, lchnk )
880 1117656 : call outfld( trim(cnst_name(mm))//'GVF', dep_grv, pcols, lchnk )
881 1117656 : call outfld( trim(cnst_name(mm))//'DTQ', ptend%q(:,:,mm), pcols, lchnk)
882 18662256 : aerdepdryis(:ncol,mm) = sflx(:ncol)
883 :
884 1117656 : else if ((lphase == 1) .and. (lspec == nspec_amode(m)+1)) then ! aerosol water
885 : ! use pvprogseasalts instead (means making the top level 0)
886 0 : pvmzaer(:ncol,1)=0._r8
887 0 : pvmzaer(:ncol,2:pverp) = vlc_dry(:ncol,:,jvlc)
888 :
889 : if(.true.) then ! use phil's method
890 : ! convert from meters/sec to pascals/sec
891 : ! pvprogseasalts(:,1) is assumed zero, use density from layer above in conversion
892 0 : pvmzaer(:ncol,2:pverp) = pvmzaer(:ncol,2:pverp) * rho(:ncol,:)*gravit
893 :
894 : ! calculate the tendencies and sfc fluxes from the above velocities
895 : call dust_sediment_tend( &
896 : ncol, dt, state%pint(:,:), state%pmid, state%pdel, state%t , &
897 0 : qaerwat(:,:,mm), pvmzaer, dqdt_tmp(:,:), sflx )
898 : else !use charlie's method
899 : call d3ddflux( ncol, vlc_dry(:,:,jvlc), qaerwat(:,:,mm), state%pmid, &
900 : state%pdel, tvs, sflx, dqdt_tmp(:,:), dt )
901 : endif
902 :
903 : ! apportion dry deposition into turb and gravitational settling for tapes
904 0 : dep_trb = 0._r8
905 0 : dep_grv = 0._r8
906 0 : do i=1,ncol
907 0 : if (vlc_dry(i,pver,jvlc) /= 0._r8) then
908 0 : dep_trb(i)=sflx(i)*vlc_trb(i,jvlc)/vlc_dry(i,pver,jvlc)
909 0 : dep_grv(i)=sflx(i)*vlc_grv(i,pver,jvlc)/vlc_dry(i,pver,jvlc)
910 : end if
911 : enddo
912 :
913 0 : qaerwat(1:ncol,:,mm) = qaerwat(1:ncol,:,mm) + dqdt_tmp(1:ncol,:) * dt
914 :
915 : else ! lphase == 2
916 : ! use pvprogseasalts instead (means making the top level 0)
917 18662256 : pvmzaer(:ncol,1)=0._r8
918 1736707464 : pvmzaer(:ncol,2:pverp) = vlc_dry(:ncol,:,jvlc)
919 1117656 : fldcw => qqcw_get_field(pbuf, mm,lchnk)
920 :
921 : if(.true.) then ! use phil's method
922 : ! convert from meters/sec to pascals/sec
923 : ! pvprogseasalts(:,1) is assumed zero, use density from layer above in conversion
924 1736707464 : pvmzaer(:ncol,2:pverp) = pvmzaer(:ncol,2:pverp) * rho(:ncol,:)*gravit
925 :
926 : ! calculate the tendencies and sfc fluxes from the above velocities
927 : call dust_sediment_tend( &
928 : ncol, dt, state%pint(:,:), state%pmid, state%pdel, state%t , &
929 1117656 : fldcw(:,:), pvmzaer, dqdt_tmp(:,:), sflx )
930 : else !use charlie's method
931 : call d3ddflux( ncol, vlc_dry(:,:,jvlc), fldcw(:,:), state%pmid, &
932 : state%pdel, tvs, sflx, dqdt_tmp(:,:), dt )
933 : endif
934 :
935 : ! apportion dry deposition into turb and gravitational settling for tapes
936 1117656 : dep_trb = 0._r8
937 1117656 : dep_grv = 0._r8
938 18662256 : do i=1,ncol
939 18662256 : if (vlc_dry(i,pver,jvlc) /= 0._r8) then
940 17544600 : dep_trb(i)=sflx(i)*vlc_trb(i,jvlc)/vlc_dry(i,pver,jvlc)
941 17544600 : dep_grv(i)=sflx(i)*vlc_grv(i,pver,jvlc)/vlc_dry(i,pver,jvlc)
942 : end if
943 : enddo
944 :
945 1736707464 : fldcw(1:ncol,:) = fldcw(1:ncol,:) + dqdt_tmp(1:ncol,:) * dt
946 :
947 1117656 : call outfld( trim(cnst_name_cw(mm))//'DDF', sflx, pcols, lchnk)
948 1117656 : call outfld( trim(cnst_name_cw(mm))//'TBF', dep_trb, pcols, lchnk )
949 1117656 : call outfld( trim(cnst_name_cw(mm))//'GVF', dep_grv, pcols, lchnk )
950 18662256 : aerdepdrycw(:ncol,mm) = sflx(:ncol)
951 :
952 : endif
953 :
954 : enddo ! lspec = 0, nspec_amode(m)+1
955 : enddo ! lphase = 1, 2
956 : enddo ! m = 1, ntot_amode
957 :
958 : ! if the user has specified prescribed aerosol dep fluxes then
959 : ! do not set cam_out dep fluxes according to the prognostic aerosols
960 58824 : if (.not.aerodep_flx_prescribed()) then
961 58824 : call set_srf_drydep(aerdepdryis, aerdepdrycw, cam_out)
962 : endif
963 :
964 117648 : endsubroutine aero_model_drydep
965 :
966 : !=============================================================================
967 : !=============================================================================
968 2529432 : subroutine aero_model_wetdep( state, dt, dlf, cam_out, ptend, pbuf)
969 :
970 58824 : use modal_aero_deposition, only: set_srf_wetdep
971 : use wetdep, only: wetdepa_v2, wetdep_inputs_set, wetdep_inputs_t
972 : use modal_aero_data
973 : use modal_aero_calcsize, only: modal_aero_calcsize_sub
974 : use modal_aero_wateruptake,only: modal_aero_wateruptake_dr
975 : use modal_aero_convproc, only: deepconv_wetdep_history, ma_convproc_intr, convproc_do_evaprain_atonce
976 :
977 : ! args
978 :
979 : type(physics_state), intent(in) :: state ! Physics state variables
980 : real(r8), intent(in) :: dt ! time step
981 : real(r8), intent(in) :: dlf(:,:) ! shallow+deep convective detrainment [kg/kg/s]
982 : type(cam_out_t), intent(inout) :: cam_out ! export state
983 : type(physics_ptend), intent(out) :: ptend ! indivdual parameterization tendencies
984 : type(physics_buffer_desc), pointer :: pbuf(:)
985 :
986 : ! local vars
987 :
988 : integer :: m ! tracer index
989 :
990 : integer :: lchnk ! chunk identifier
991 : integer :: ncol ! number of atmospheric columns
992 :
993 : real(r8) :: iscavt(pcols, pver)
994 :
995 : integer :: mm
996 : integer :: i,k
997 :
998 : real(r8) :: icscavt(pcols, pver)
999 : real(r8) :: isscavt(pcols, pver)
1000 : real(r8) :: bcscavt(pcols, pver)
1001 : real(r8) :: bsscavt(pcols, pver)
1002 : real(r8) :: sol_factb, sol_facti
1003 : real(r8) :: sol_factic(pcols,pver)
1004 :
1005 : real(r8) :: sflx(pcols) ! deposition flux
1006 :
1007 : integer :: jnv ! index for scavcoefnv 3rd dimension
1008 : integer :: lphase ! index for interstitial / cloudborne aerosol
1009 : integer :: strt_loop, end_loop, stride_loop !loop indices for the lphase loop
1010 : integer :: lspec ! index for aerosol number / chem-mass / water-mass
1011 : integer :: lcoardust, lcoarnacl ! indices for coarse mode dust and seasalt masses
1012 : real(r8) :: dqdt_tmp(pcols,pver) ! temporary array to hold tendency for 1 species
1013 : real(r8) :: f_act_conv(pcols,pver) ! prescribed aerosol activation fraction for convective cloud ! rce 2010/05/01
1014 : real(r8) :: f_act_conv_coarse(pcols,pver) ! similar but for coarse mode ! rce 2010/05/02
1015 : real(r8) :: f_act_conv_coarse_dust, f_act_conv_coarse_nacl ! rce 2010/05/02
1016 : real(r8) :: fracis_cw(pcols,pver)
1017 : real(r8) :: hygro_sum_old(pcols,pver) ! before removal [sum of (mass*hydro/dens)]
1018 : real(r8) :: hygro_sum_del(pcols,pver) ! removal change to [sum of (mass*hydro/dens)]
1019 : real(r8) :: hygro_sum_old_ik, hygro_sum_new_ik
1020 : real(r8) :: prec(pcols) ! precipitation rate
1021 : real(r8) :: q_tmp(pcols,pver) ! temporary array to hold "most current" mixing ratio for 1 species
1022 : real(r8) :: scavcoefnv(pcols,pver,0:2) ! Dana and Hales coefficient (/mm) for
1023 : ! cloud-borne num & vol (0),
1024 : ! interstitial num (1), interstitial vol (2)
1025 : real(r8) :: tmpa, tmpb
1026 : real(r8) :: tmpdust, tmpnacl
1027 : real(r8) :: water_old, water_new ! temporary old/new aerosol water mix-rat
1028 : logical :: isprx(pcols,pver) ! true if precipation
1029 : real(r8) :: aerdepwetis(pcols,pcnst) ! aerosol wet deposition (interstitial)
1030 : real(r8) :: aerdepwetcw(pcols,pcnst) ! aerosol wet deposition (cloud water)
1031 :
1032 : ! For unified convection scheme
1033 : logical, parameter :: do_aero_water_removal = .false. ! True if aerosol water reduction by wet removal is to be calculated
1034 : ! (this has not been fully tested, so best to leave it off)
1035 : logical :: do_hygro_sum_del, do_lphase1, do_lphase2
1036 :
1037 58824 : real(r8), pointer :: rprddp(:,:) ! rain production, deep convection
1038 58824 : real(r8), pointer :: rprdsh(:,:) ! rain production, shallow convection
1039 58824 : real(r8), pointer :: evapcdp(:,:) ! Evaporation rate of deep convective precipitation >=0.
1040 58824 : real(r8), pointer :: evapcsh(:,:) ! Evaporation rate of shallow convective precipitation >=0.
1041 :
1042 : real(r8) :: rprddpsum(pcols)
1043 : real(r8) :: rprdshsum(pcols)
1044 : real(r8) :: evapcdpsum(pcols)
1045 : real(r8) :: evapcshsum(pcols)
1046 :
1047 : real(r8) :: tmp_resudp, tmp_resush
1048 :
1049 : real(r8) :: sflxec(pcols), sflxecdp(pcols) ! deposition flux
1050 : real(r8) :: sflxic(pcols), sflxicdp(pcols) ! deposition flux
1051 : real(r8) :: sflxbc(pcols), sflxbcdp(pcols) ! deposition flux
1052 : real(r8) :: rcscavt(pcols, pver)
1053 : real(r8) :: rsscavt(pcols, pver)
1054 117648 : real(r8) :: qqcw_in(pcols,pver), qqcw_sav(pcols,pver,0:nspec_max) ! temporary array to hold qqcw for the current mode
1055 117648 : real(r8) :: rtscavt(pcols, pver, 0:nspec_max)
1056 :
1057 : integer, parameter :: nsrflx_mzaer2cnvpr = 2
1058 : real(r8) :: qsrflx_mzaer2cnvpr(pcols,pcnst,nsrflx_mzaer2cnvpr)
1059 : ! End unified convection scheme
1060 :
1061 58824 : real(r8), pointer :: fldcw(:,:)
1062 :
1063 58824 : real(r8), pointer :: dgnumwet(:,:,:)
1064 58824 : real(r8), pointer :: qaerwat(:,:,:) ! aerosol water
1065 :
1066 58824 : real(r8), pointer :: fracis(:,:,:) ! fraction of transported species that are insoluble
1067 :
1068 : type(wetdep_inputs_t) :: dep_inputs
1069 :
1070 : real(r8) :: dcondt_resusp3d(2*pcnst,pcols, pver)
1071 :
1072 58824 : lchnk = state%lchnk
1073 58824 : ncol = state%ncol
1074 :
1075 58824 : dcondt_resusp3d(:,:,:) = 0._r8
1076 :
1077 58824 : call physics_ptend_init(ptend, state%psetcols, 'aero_model_wetdep', lq=wetdep_lq)
1078 :
1079 : ! Do calculations of mode radius and water uptake if:
1080 : ! 1) modal aerosols are affecting the climate, or
1081 : ! 2) prognostic modal aerosols are enabled
1082 :
1083 58824 : call t_startf('calcsize')
1084 : ! for prognostic modal aerosols the transfer of mass between aitken and accumulation
1085 : ! modes is done in conjunction with the dry radius calculation
1086 58824 : call modal_aero_calcsize_sub(state, ptend, dt, pbuf)
1087 58824 : call t_stopf('calcsize')
1088 :
1089 58824 : call t_startf('wateruptake')
1090 58824 : call modal_aero_wateruptake_dr(state, pbuf)
1091 58824 : call t_stopf('wateruptake')
1092 :
1093 58824 : if (nwetdep<1) return
1094 :
1095 58824 : call wetdep_inputs_set( state, pbuf, dep_inputs )
1096 :
1097 235296 : call pbuf_get_field(pbuf, dgnumwet_idx, dgnumwet, start=(/1,1,1/), kount=(/pcols,pver,nmodes/) )
1098 235296 : call pbuf_get_field(pbuf, qaerwat_idx, qaerwat, start=(/1,1,1/), kount=(/pcols,pver,nmodes/) )
1099 58824 : call pbuf_get_field(pbuf, fracis_idx, fracis, start=(/1,1,1/), kount=(/pcols, pver, pcnst/) )
1100 :
1101 982224 : prec(:ncol)=0._r8
1102 5529456 : do k=1,pver
1103 91346832 : where (prec(:ncol) >= 1.e-7_r8)
1104 : isprx(:ncol,k) = .true.
1105 : elsewhere
1106 5470632 : isprx(:ncol,k) = .false.
1107 : endwhere
1108 0 : prec(:ncol) = prec(:ncol) + (dep_inputs%prain(:ncol,k) + dep_inputs%cmfdqr(:ncol,k) - dep_inputs%evapr(:ncol,k)) &
1109 91405656 : *state%pdel(:ncol,k)/gravit
1110 : end do
1111 :
1112 58824 : if(convproc_do_aer) then
1113 58824 : qsrflx_mzaer2cnvpr(:,:,:) = 0.0_r8
1114 58824 : aerdepwetis(:,:) = 0.0_r8
1115 58824 : aerdepwetcw(:,:) = 0.0_r8
1116 : else
1117 0 : qsrflx_mzaer2cnvpr(:,:,:) = nan
1118 0 : aerdepwetis(:,:) = nan
1119 0 : aerdepwetcw(:,:) = nan
1120 : endif
1121 :
1122 : ! calculate the mass-weighted sol_factic for coarse mode species
1123 : ! sol_factic_coarse(:,:) = 0.30_r8 ! tuned 1/4
1124 93059568 : f_act_conv_coarse(:,:) = 0.60_r8 ! rce 2010/05/02
1125 58824 : f_act_conv_coarse_dust = 0.40_r8 ! rce 2010/05/02
1126 58824 : f_act_conv_coarse_nacl = 0.80_r8 ! rce 2010/05/02
1127 58824 : if (modeptr_coarse > 0) then
1128 58824 : lcoardust = lptr_dust_a_amode(modeptr_coarse)
1129 58824 : lcoarnacl = lptr_nacl_a_amode(modeptr_coarse)
1130 58824 : if ((lcoardust > 0) .and. (lcoarnacl > 0)) then
1131 5529456 : do k = 1, pver
1132 91405656 : do i = 1, ncol
1133 85876200 : tmpdust = max( 0.0_r8, state%q(i,k,lcoardust) + ptend%q(i,k,lcoardust)*dt )
1134 85876200 : tmpnacl = max( 0.0_r8, state%q(i,k,lcoarnacl) + ptend%q(i,k,lcoarnacl)*dt )
1135 91346832 : if ((tmpdust+tmpnacl) > 1.0e-30_r8) then
1136 : ! sol_factic_coarse(i,k) = (0.2_r8*tmpdust + 0.4_r8*tmpnacl)/(tmpdust+tmpnacl) ! tuned 1/6
1137 79839942 : f_act_conv_coarse(i,k) = (f_act_conv_coarse_dust*tmpdust &
1138 79839942 : + f_act_conv_coarse_nacl*tmpnacl)/(tmpdust+tmpnacl) ! rce 2010/05/02
1139 : end if
1140 : end do
1141 : end do
1142 : end if
1143 : end if
1144 :
1145 93059568 : scavcoefnv(:,:,0) = 0.0_r8 ! below-cloud scavcoef = 0.0 for cloud-borne species
1146 :
1147 : ! Counters for "without" unified convective treatment (i.e. default case)
1148 58824 : strt_loop = 1
1149 58824 : end_loop = 2
1150 58824 : stride_loop = 1
1151 58824 : if (convproc_do_aer) then
1152 : !Do cloudborne first for unified convection scheme so that the resuspension of cloudborne
1153 : !can be saved then applied to interstitial
1154 58824 : strt_loop = 2
1155 58824 : end_loop = 1
1156 58824 : stride_loop = -1
1157 : endif
1158 :
1159 294120 : do m = 1, ntot_amode ! main loop over aerosol modes
1160 :
1161 764712 : do lphase = strt_loop,end_loop, stride_loop ! loop over interstitial (1) and cloud-borne (2) forms
1162 :
1163 : ! sol_factb and sol_facti values
1164 : ! sol_factb - currently this is basically a tuning factor
1165 : ! sol_facti & sol_factic - currently has a physical basis, and reflects activation fraction
1166 : !
1167 : ! 2008-mar-07 rce - sol_factb (interstitial) changed from 0.3 to 0.1
1168 : ! - sol_factic (interstitial, dust modes) changed from 1.0 to 0.5
1169 : ! - sol_factic (cloud-borne, pcarb modes) no need to set it to 0.0
1170 : ! because the cloud-borne pcarbon == 0 (no activation)
1171 : !
1172 : ! rce 2010/05/02
1173 : ! prior to this date, sol_factic was used for convective in-cloud wet removal,
1174 : ! and its value reflected a combination of an activation fraction (which varied between modes)
1175 : ! and a tuning factor
1176 : ! from this date forward, two parameters are used for convective in-cloud wet removal
1177 : ! f_act_conv is the activation fraction
1178 : ! note that "non-activation" of aerosol in air entrained into updrafts should
1179 : ! be included here
1180 : ! eventually we might use the activate routine (with w ~= 1 m/s) to calculate
1181 : ! this, but there is still the entrainment issue
1182 : ! sol_factic is strictly a tuning factor
1183 : !
1184 470592 : if (lphase == 1) then ! interstial aerosol
1185 235296 : hygro_sum_old(:,:) = 0.0_r8
1186 235296 : hygro_sum_del(:,:) = 0.0_r8
1187 : call modal_aero_bcscavcoef_get( m, ncol, isprx, dgnumwet, &
1188 235296 : scavcoefnv(:,:,1), scavcoefnv(:,:,2) )
1189 :
1190 235296 : sol_factb = sol_factb_interstitial ! all below-cloud scav ON (0.1 "tuning factor")
1191 :
1192 235296 : sol_facti = 0.0_r8 ! strat in-cloud scav totally OFF for institial
1193 :
1194 372238272 : sol_factic = sol_factic_interstitial
1195 :
1196 235296 : if (m == modeptr_pcarbon) then
1197 : ! sol_factic = 0.0_r8 ! conv in-cloud scav OFF (0.0 activation fraction)
1198 58824 : f_act_conv = 0.0_r8 ! rce 2010/05/02
1199 176472 : else if ((m == modeptr_finedust) .or. (m == modeptr_coardust)) then
1200 : ! sol_factic = 0.2_r8 ! conv in-cloud scav ON (0.5 activation fraction) ! tuned 1/4
1201 0 : f_act_conv = 0.4_r8 ! rce 2010/05/02
1202 : else
1203 : ! sol_factic = 0.4_r8 ! conv in-cloud scav ON (1.0 activation fraction) ! tuned 1/4
1204 279178704 : f_act_conv = 0.8_r8 ! rce 2010/05/02
1205 : end if
1206 :
1207 : else ! cloud-borne aerosol (borne by stratiform cloud drops)
1208 :
1209 235296 : sol_factb = 0.0_r8 ! all below-cloud scav OFF (anything cloud-borne is located "in-cloud")
1210 235296 : sol_facti = sol_facti_cloud_borne ! strat in-cloud scav cloud-borne tuning factor
1211 235296 : sol_factic = 0.0_r8 ! conv in-cloud scav OFF (having this on would mean
1212 : ! that conv precip collects strat droplets)
1213 235296 : f_act_conv = 0.0_r8 ! conv in-cloud scav OFF (having this on would mean
1214 :
1215 : end if
1216 470592 : if (convproc_do_aer .and. lphase == 1) then
1217 : ! if modal aero convproc is turned on for aerosols, then
1218 : ! turn off the convective in-cloud removal for interstitial aerosols
1219 : ! (but leave the below-cloud on, as convproc only does in-cloud)
1220 : ! and turn off the outfld SFWET, SFSIC, SFSID, SFSEC, and SFSED calls
1221 : ! for (stratiform)-cloudborne aerosols, convective wet removal
1222 : ! (all forms) is zero, so no action is needed
1223 235296 : sol_factic = 0.0_r8
1224 : endif
1225 :
1226 : !
1227 : ! rce 2010/05/03
1228 : ! wetdepa has "sol_fact" parameters:
1229 : ! sol_facti, sol_factic, sol_factb for liquid cloud
1230 :
1231 3411792 : do lspec = 0, nspec_amode(m)+1 ! loop over number + chem constituents + water
1232 :
1233 2705904 : if (lspec == 0) then ! number
1234 470592 : if (lphase == 1) then
1235 235296 : mm = numptr_amode(m)
1236 235296 : jnv = 1
1237 : else
1238 235296 : mm = numptrcw_amode(m)
1239 235296 : jnv = 0
1240 : endif
1241 2235312 : else if (lspec <= nspec_amode(m)) then ! non-water mass
1242 1764720 : if (lphase == 1) then
1243 882360 : mm = lmassptr_amode(lspec,m)
1244 882360 : jnv = 2
1245 : else
1246 882360 : mm = lmassptrcw_amode(lspec,m)
1247 882360 : jnv = 0
1248 : endif
1249 : else ! water mass
1250 : ! bypass wet removal of aerosol water
1251 : if(convproc_do_aer) then
1252 : if ( .not. do_aero_water_removal ) cycle
1253 : else
1254 : cycle
1255 : endif
1256 : if (lphase == 1) then
1257 : mm = 0
1258 : ! mm = lwaterptr_amode(m)
1259 : jnv = 2
1260 : else
1261 : mm = 0
1262 : jnv = 0
1263 : endif
1264 : endif
1265 :
1266 2235312 : if (mm <= 0) cycle
1267 :
1268 :
1269 : ! set f_act_conv for interstitial (lphase=1) coarse mode species
1270 : ! for the convective in-cloud, we conceptually treat the coarse dust and seasalt
1271 : ! as being externally mixed, and apply f_act_conv = f_act_conv_coarse_dust/nacl to dust/seasalt
1272 : ! number and sulfate are conceptually partitioned to the dust and seasalt
1273 : ! on a mass basis, so the f_act_conv for number and sulfate are
1274 : ! mass-weighted averages of the values used for dust/seasalt
1275 2235312 : if ((lphase == 1) .and. (m == modeptr_coarse)) then
1276 : ! sol_factic = sol_factic_coarse
1277 235296 : f_act_conv = f_act_conv_coarse ! rce 2010/05/02
1278 235296 : if (lspec > 0) then
1279 176472 : if (lmassptr_amode(lspec,m) == lptr_dust_a_amode(m)) then
1280 : ! sol_factic = 0.2_r8 ! tuned 1/4
1281 93059568 : f_act_conv = f_act_conv_coarse_dust ! rce 2010/05/02
1282 117648 : else if (lmassptr_amode(lspec,m) == lptr_nacl_a_amode(m)) then
1283 : ! sol_factic = 0.4_r8 ! tuned 1/6
1284 93059568 : f_act_conv = f_act_conv_coarse_nacl ! rce 2010/05/02
1285 : end if
1286 : end if
1287 : end if
1288 :
1289 2705904 : if ((lphase == 1) .and. (lspec <= nspec_amode(m))) then
1290 1117656 : ptend%lq(mm) = .TRUE.
1291 1117656 : dqdt_tmp(:,:) = 0.0_r8
1292 : ! q_tmp reflects changes from modal_aero_calcsize and is the "most current" q
1293 1736707464 : q_tmp(1:ncol,:) = state%q(1:ncol,:,mm) + ptend%q(1:ncol,:,mm)*dt
1294 1117656 : if(convproc_do_aer) then
1295 : !Feed in the saved cloudborne mixing ratios from phase 2
1296 1768131792 : qqcw_in(:,:) = qqcw_sav(:,:,lspec)
1297 : else
1298 0 : fldcw => qqcw_get_field(pbuf, mm,lchnk)
1299 0 : qqcw_in(:,:) = fldcw(:,:)
1300 : endif
1301 :
1302 : call wetdepa_v2( state%pmid, state%q(:,:,1), state%pdel, &
1303 : dep_inputs%cldt, dep_inputs%cldcu, dep_inputs%cmfdqr, &
1304 : dep_inputs%evapc, dep_inputs%conicw, dep_inputs%prain, dep_inputs%qme, &
1305 : dep_inputs%evapr, dep_inputs%totcond, q_tmp, dt, &
1306 : dqdt_tmp, iscavt, dep_inputs%cldvcu, dep_inputs%cldvst, &
1307 : dlf, fracis(:,:,mm), sol_factb, ncol, &
1308 : scavcoefnv(:,:,jnv), &
1309 : is_strat_cloudborne=.false., &
1310 : qqcw=qqcw_in(:,:), &
1311 : f_act_conv=f_act_conv, &
1312 : icscavt=icscavt, isscavt=isscavt, bcscavt=bcscavt, bsscavt=bsscavt, &
1313 : convproc_do_aer=convproc_do_aer, rcscavt=rcscavt, rsscavt=rsscavt, &
1314 : sol_facti_in=sol_facti, sol_factic_in=sol_factic, &
1315 1117656 : convproc_do_evaprain_atonce_in=convproc_do_evaprain_atonce )
1316 :
1317 1117656 : do_hygro_sum_del = .false.
1318 1117656 : if ( lspec > 0 ) do_hygro_sum_del = .true.
1319 :
1320 1117656 : if(convproc_do_aer) then
1321 1117656 : do_hygro_sum_del = .false.
1322 : ! add resuspension of cloudborne species to dqdt of interstitial species
1323 1736707464 : dqdt_tmp(1:ncol,:) = dqdt_tmp(1:ncol,:) + rtscavt(1:ncol,:,lspec)
1324 : if ( (lspec > 0) .and. do_aero_water_removal ) then
1325 : do_hygro_sum_del = .true.
1326 : endif
1327 : endif
1328 :
1329 1736707464 : ptend%q(1:ncol,:,mm) = ptend%q(1:ncol,:,mm) + dqdt_tmp(1:ncol,:)
1330 :
1331 1117656 : call outfld( trim(cnst_name(mm))//'WET', dqdt_tmp(:,:), pcols, lchnk)
1332 1117656 : call outfld( trim(cnst_name(mm))//'SIC', icscavt, pcols, lchnk)
1333 1117656 : call outfld( trim(cnst_name(mm))//'SIS', isscavt, pcols, lchnk)
1334 1117656 : call outfld( trim(cnst_name(mm))//'SBC', bcscavt, pcols, lchnk)
1335 1117656 : call outfld( trim(cnst_name(mm))//'SBS', bsscavt, pcols, lchnk)
1336 :
1337 1117656 : sflx(:)=0._r8
1338 105059664 : do k=1,pver
1339 1736707464 : do i=1,ncol
1340 1735589808 : sflx(i)=sflx(i)+dqdt_tmp(i,k)*state%pdel(i,k)/gravit
1341 : enddo
1342 : enddo
1343 1117656 : if (.not.convproc_do_aer) call outfld( trim(cnst_name(mm))//'SFWET', sflx, pcols, lchnk)
1344 18662256 : aerdepwetis(:ncol,mm) = sflx(:ncol)
1345 :
1346 1117656 : sflx(:)=0._r8
1347 105059664 : do k=1,pver
1348 1736707464 : do i=1,ncol
1349 1735589808 : sflx(i)=sflx(i)+icscavt(i,k)*state%pdel(i,k)/gravit
1350 : enddo
1351 : enddo
1352 1117656 : if (.not.convproc_do_aer) call outfld( trim(cnst_name(mm))//'SFSIC', sflx, pcols, lchnk)
1353 1117656 : if (convproc_do_aer) sflxic = sflx
1354 :
1355 1117656 : sflx(:)=0._r8
1356 105059664 : do k=1,pver
1357 1736707464 : do i=1,ncol
1358 1735589808 : sflx(i)=sflx(i)+isscavt(i,k)*state%pdel(i,k)/gravit
1359 : enddo
1360 : enddo
1361 1117656 : call outfld( trim(cnst_name(mm))//'SFSIS', sflx, pcols, lchnk)
1362 :
1363 1117656 : sflx(:)=0._r8
1364 105059664 : do k=1,pver
1365 1736707464 : do i=1,ncol
1366 1735589808 : sflx(i)=sflx(i)+bcscavt(i,k)*state%pdel(i,k)/gravit
1367 : enddo
1368 : enddo
1369 1117656 : call outfld( trim(cnst_name(mm))//'SFSBC', sflx, pcols, lchnk)
1370 1117656 : if (convproc_do_aer)sflxbc = sflx
1371 :
1372 1117656 : sflx(:)=0._r8
1373 105059664 : do k=1,pver
1374 1736707464 : do i=1,ncol
1375 1735589808 : sflx(i)=sflx(i)+bsscavt(i,k)*state%pdel(i,k)/gravit
1376 : enddo
1377 : enddo
1378 1117656 : call outfld( trim(cnst_name(mm))//'SFSBS', sflx, pcols, lchnk)
1379 :
1380 1117656 : if (convproc_do_aer) then
1381 :
1382 1117656 : sflx(:)=0._r8
1383 105059664 : do k=1,pver
1384 1736707464 : do i=1,ncol
1385 1735589808 : sflx(i)=sflx(i)+rcscavt(i,k)*state%pdel(i,k)/gravit
1386 : enddo
1387 : enddo
1388 1117656 : sflxec = sflx
1389 :
1390 1117656 : sflx(:)=0._r8
1391 105059664 : do k=1,pver
1392 1736707464 : do i=1,ncol
1393 1735589808 : sflx(i)=sflx(i)+rsscavt(i,k)*state%pdel(i,k)/gravit
1394 : enddo
1395 : enddo
1396 1117656 : call outfld( trim(cnst_name(mm))//'SFSES', sflx, pcols, lchnk)
1397 :
1398 : ! apportion convective surface fluxes to deep and shallow conv
1399 : ! this could be done more accurately in subr wetdepa
1400 : ! since deep and shallow rarely occur simultaneously, and these
1401 : ! fields are just diagnostics, this approximate method is adequate
1402 : ! only do this for interstitial aerosol, because conv clouds to not
1403 : ! affect the stratiform-cloudborne aerosol
1404 1117656 : if ( deepconv_wetdep_history) then
1405 :
1406 1117656 : call pbuf_get_field(pbuf, rprddp_idx, rprddp )
1407 1117656 : call pbuf_get_field(pbuf, rprdsh_idx, rprdsh )
1408 1117656 : call pbuf_get_field(pbuf, nevapr_dpcu_idx, evapcdp )
1409 1117656 : call pbuf_get_field(pbuf, nevapr_shcu_idx, evapcsh )
1410 :
1411 1117656 : rprddpsum(:) = 0.0_r8
1412 1117656 : rprdshsum(:) = 0.0_r8
1413 1117656 : evapcdpsum(:) = 0.0_r8
1414 1117656 : evapcshsum(:) = 0.0_r8
1415 :
1416 105059664 : do k = 1, pver
1417 1735589808 : rprddpsum(:ncol) = rprddpsum(:ncol) + rprddp(:ncol,k)*state%pdel(:ncol,k)/gravit
1418 1735589808 : rprdshsum(:ncol) = rprdshsum(:ncol) + rprdsh(:ncol,k)*state%pdel(:ncol,k)/gravit
1419 1735589808 : evapcdpsum(:ncol) = evapcdpsum(:ncol) + evapcdp(:ncol,k)*state%pdel(:ncol,k)/gravit
1420 1736707464 : evapcshsum(:ncol) = evapcshsum(:ncol) + evapcsh(:ncol,k)*state%pdel(:ncol,k)/gravit
1421 : end do
1422 :
1423 18662256 : do i = 1, ncol
1424 17544600 : rprddpsum(i) = max( rprddpsum(i), 1.0e-35_r8 )
1425 17544600 : rprdshsum(i) = max( rprdshsum(i), 1.0e-35_r8 )
1426 17544600 : evapcdpsum(i) = max( evapcdpsum(i), 0.1e-35_r8 )
1427 17544600 : evapcshsum(i) = max( evapcshsum(i), 0.1e-35_r8 )
1428 :
1429 : ! assume that in- and below-cloud removal are proportional to column precip production
1430 17544600 : tmpa = rprddpsum(i) / (rprddpsum(i) + rprdshsum(i))
1431 17544600 : tmpa = max( 0.0_r8, min( 1.0_r8, tmpa ) )
1432 17544600 : sflxicdp(i) = sflxic(i)*tmpa
1433 17544600 : sflxbcdp(i) = sflxbc(i)*tmpa
1434 :
1435 : ! assume that resuspension is proportional to (wet removal)*[(precip evap)/(precip production)]
1436 17544600 : tmp_resudp = tmpa * min( (evapcdpsum(i)/rprddpsum(i)), 1.0_r8 )
1437 17544600 : tmp_resush = (1.0_r8 - tmpa) * min( (evapcshsum(i)/rprdshsum(i)), 1.0_r8 )
1438 17544600 : tmpb = max( tmp_resudp, 1.0e-35_r8 ) / max( (tmp_resudp+tmp_resush), 1.0e-35_r8 )
1439 17544600 : tmpb = max( 0.0_r8, min( 1.0_r8, tmpb ) )
1440 18662256 : sflxecdp(i) = sflxec(i)*tmpb
1441 : end do
1442 1117656 : call outfld( trim(cnst_name(mm))//'SFSBD', sflxbcdp, pcols, lchnk)
1443 : else
1444 0 : sflxec(1:ncol) = 0.0_r8
1445 0 : sflxecdp(1:ncol) = 0.0_r8
1446 : end if
1447 :
1448 : ! when ma_convproc_intr is used, convective in-cloud wet removal is done there
1449 : ! the convective (total and deep) precip-evap-resuspension includes in- and below-cloud
1450 : ! contributions
1451 : ! so pass the below-cloud contribution to ma_convproc_intr
1452 18662256 : qsrflx_mzaer2cnvpr(1:ncol,mm,1) = sflxec( 1:ncol)
1453 18662256 : qsrflx_mzaer2cnvpr(1:ncol,mm,2) = sflxecdp(1:ncol)
1454 :
1455 : endif
1456 :
1457 1117656 : if (do_hygro_sum_del) then
1458 0 : tmpa = spechygro(lspec,m)/ &
1459 0 : specdens_amode(lspec,m)
1460 0 : tmpb = tmpa*dt
1461 0 : hygro_sum_old(1:ncol,:) = hygro_sum_old(1:ncol,:) &
1462 0 : + tmpa*q_tmp(1:ncol,:)
1463 : hygro_sum_del(1:ncol,:) = hygro_sum_del(1:ncol,:) &
1464 0 : + tmpb*dqdt_tmp(1:ncol,:)
1465 : end if
1466 :
1467 1117656 : else if ((lphase == 1) .and. (lspec == nspec_amode(m)+1)) then
1468 0 : do_lphase1 = .true.
1469 0 : if(convproc_do_aer) then
1470 : do_lphase1 = .false.
1471 : if(do_aero_water_removal)do_lphase1 = .true.
1472 : endif
1473 : if(do_lphase1) then
1474 : ! aerosol water -- because of how wetdepa treats evaporation of stratiform
1475 : ! precip, it is not appropriate to apply wetdepa to aerosol water
1476 : ! instead, "hygro_sum" = [sum of (mass*hygro/dens)] is calculated before and
1477 : ! after wet removal, and new water is calculated using
1478 : ! new_water = old_water*min(10,(hygro_sum_new/hygro_sum_old))
1479 : ! the "min(10,...)" is to avoid potential problems when hygro_sum_old ~= 0
1480 : ! also, individual wet removal terms (ic,is,bc,bs) are not output to history
1481 : ! ptend%lq(mm) = .TRUE.
1482 : ! dqdt_tmp(:,:) = 0.0_r8
1483 0 : do k = 1, pver
1484 0 : do i = 1, ncol
1485 : ! water_old = max( 0.0_r8, state%q(i,k,mm)+ptend%q(i,k,mm)*dt )
1486 0 : water_old = max( 0.0_r8, qaerwat(i,k,mm) )
1487 0 : hygro_sum_old_ik = max( 0.0_r8, hygro_sum_old(i,k) )
1488 0 : hygro_sum_new_ik = max( 0.0_r8, hygro_sum_old_ik+hygro_sum_del(i,k) )
1489 0 : if (hygro_sum_new_ik >= 10.0_r8*hygro_sum_old_ik) then
1490 0 : water_new = 10.0_r8*water_old
1491 : else
1492 0 : water_new = water_old*(hygro_sum_new_ik/hygro_sum_old_ik)
1493 : end if
1494 : ! dqdt_tmp(i,k) = (water_new - water_old)/dt
1495 0 : qaerwat(i,k,mm) = water_new
1496 : end do
1497 : end do
1498 :
1499 : ! ptend%q(1:ncol,:,mm) = ptend%q(1:ncol,:,mm) + dqdt_tmp(1:ncol,:)
1500 :
1501 : ! call outfld( trim(cnst_name(mm))
1502 :
1503 : ! sflx(:)=0._r8
1504 : ! do k=1,pver
1505 : ! do i=1,ncol
1506 : ! sflx(i)=sflx(i)+dqdt_tmp(i,k)*state%pdel(i,k)/gravit
1507 : ! enddo
1508 : ! enddo
1509 : ! call outfld( trim(cnst_name(mm))
1510 : endif
1511 :
1512 1117656 : elseif (lphase == 2) then
1513 :
1514 1117656 : do_lphase2 = .true.
1515 1117656 : if (convproc_do_aer) then
1516 1117656 : do_lphase2 = .false.
1517 1117656 : if (lspec <= nspec_amode(m)) do_lphase2 = .true.
1518 : endif
1519 :
1520 : if (do_lphase2) then
1521 :
1522 1117656 : dqdt_tmp(:,:) = 0.0_r8
1523 :
1524 1117656 : if (convproc_do_aer) then
1525 1117656 : fldcw => qqcw_get_field(pbuf,mm,lchnk)
1526 1736707464 : qqcw_sav(1:ncol,:,lspec) = fldcw(1:ncol,:)
1527 : else
1528 0 : fldcw => qqcw_get_field(pbuf, mm,lchnk)
1529 : endif
1530 :
1531 : call wetdepa_v2(state%pmid, state%q(:,:,1), state%pdel, &
1532 : dep_inputs%cldt, dep_inputs%cldcu, dep_inputs%cmfdqr, &
1533 : dep_inputs%evapc, dep_inputs%conicw, dep_inputs%prain, dep_inputs%qme, &
1534 : dep_inputs%evapr, dep_inputs%totcond, fldcw, dt, &
1535 : dqdt_tmp, iscavt, dep_inputs%cldvcu, dep_inputs%cldvst, &
1536 : dlf, fracis_cw, sol_factb, ncol, &
1537 : scavcoefnv(:,:,jnv), &
1538 : is_strat_cloudborne=.true., &
1539 : icscavt=icscavt, isscavt=isscavt, bcscavt=bcscavt, bsscavt=bsscavt, &
1540 : convproc_do_aer=convproc_do_aer, rcscavt=rcscavt, rsscavt=rsscavt, &
1541 : sol_facti_in=sol_facti, sol_factic_in=sol_factic, &
1542 : convproc_do_evaprain_atonce_in=convproc_do_evaprain_atonce, &
1543 1117656 : bergso_in=dep_inputs%bergso )
1544 :
1545 1117656 : if(convproc_do_aer) then
1546 : ! save resuspension of cloudborne species
1547 1736707464 : rtscavt(1:ncol,:,lspec) = rcscavt(1:ncol,:) + rsscavt(1:ncol,:)
1548 : ! wetdepa_v2 adds the resuspension of cloudborne to the dqdt of cloudborne (as a source)
1549 : ! undo this, so the resuspension of cloudborne can be added to the dqdt of interstitial (above)
1550 1736707464 : dqdt_tmp(1:ncol,:) = dqdt_tmp(1:ncol,:) - rtscavt(1:ncol,:,lspec)
1551 : endif
1552 :
1553 :
1554 1736707464 : fldcw(1:ncol,:) = fldcw(1:ncol,:) + dqdt_tmp(1:ncol,:) * dt
1555 :
1556 1117656 : sflx(:)=0._r8
1557 105059664 : do k=1,pver
1558 1736707464 : do i=1,ncol
1559 1735589808 : sflx(i)=sflx(i)+dqdt_tmp(i,k)*state%pdel(i,k)/gravit
1560 : enddo
1561 : enddo
1562 1117656 : call outfld( trim(cnst_name_cw(mm))//'SFWET', sflx, pcols, lchnk)
1563 18662256 : aerdepwetcw(:ncol,mm) = sflx(:ncol)
1564 :
1565 1117656 : sflx(:)=0._r8
1566 105059664 : do k=1,pver
1567 1736707464 : do i=1,ncol
1568 1735589808 : sflx(i)=sflx(i)+icscavt(i,k)*state%pdel(i,k)/gravit
1569 : enddo
1570 : enddo
1571 1117656 : call outfld( trim(cnst_name_cw(mm))//'SFSIC', sflx, pcols, lchnk)
1572 1117656 : sflx(:)=0._r8
1573 105059664 : do k=1,pver
1574 1736707464 : do i=1,ncol
1575 1735589808 : sflx(i)=sflx(i)+isscavt(i,k)*state%pdel(i,k)/gravit
1576 : enddo
1577 : enddo
1578 1117656 : call outfld( trim(cnst_name_cw(mm))//'SFSIS', sflx, pcols, lchnk)
1579 1117656 : sflx(:)=0._r8
1580 105059664 : do k=1,pver
1581 1736707464 : do i=1,ncol
1582 1735589808 : sflx(i)=sflx(i)+bcscavt(i,k)*state%pdel(i,k)/gravit
1583 : enddo
1584 : enddo
1585 1117656 : call outfld( trim(cnst_name_cw(mm))//'SFSBC', sflx, pcols, lchnk)
1586 1117656 : sflx(:)=0._r8
1587 105059664 : do k=1,pver
1588 1736707464 : do i=1,ncol
1589 1735589808 : sflx(i)=sflx(i)+bsscavt(i,k)*state%pdel(i,k)/gravit
1590 : enddo
1591 : enddo
1592 1117656 : call outfld( trim(cnst_name_cw(mm))//'SFSBS', sflx, pcols, lchnk)
1593 :
1594 1117656 : if(convproc_do_aer) then
1595 1117656 : sflx(:)=0.0_r8
1596 105059664 : do k=1,pver
1597 1736707464 : sflx(1:ncol)=sflx(1:ncol)+rcscavt(1:ncol,k)*state%pdel(1:ncol,k)/gravit
1598 : enddo
1599 1117656 : call outfld( trim(cnst_name_cw(mm))//'SFSEC', sflx, pcols, lchnk)
1600 :
1601 1117656 : sflx(:)=0.0_r8
1602 105059664 : do k=1,pver
1603 1736707464 : sflx(1:ncol)=sflx(1:ncol)+rsscavt(1:ncol,k)*state%pdel(1:ncol,k)/gravit
1604 : enddo
1605 1117656 : call outfld( trim(cnst_name_cw(mm))//'SFSES', sflx, pcols, lchnk)
1606 : endif
1607 : endif
1608 : endif
1609 :
1610 : enddo ! lspec = 0, nspec_amode(m)+1
1611 : enddo ! lphase = 1, 2
1612 : enddo ! m = 1, ntot_amode
1613 :
1614 58824 : if (convproc_do_aer) then
1615 58824 : call t_startf('ma_convproc')
1616 : call ma_convproc_intr( state, ptend, pbuf, dt, &
1617 : nsrflx_mzaer2cnvpr, qsrflx_mzaer2cnvpr, aerdepwetis, &
1618 58824 : dcondt_resusp3d)
1619 :
1620 58824 : if (convproc_do_evaprain_atonce) then
1621 294120 : do m = 1, ntot_amode ! main loop over aerosol modes
1622 764712 : do lphase = strt_loop,end_loop, stride_loop
1623 : ! loop over interstitial (1) and cloud-borne (2) forms
1624 3411792 : do lspec = 0, nspec_amode(m)+1 ! loop over number + chem constituents + water
1625 2705904 : if (lspec == 0) then ! number
1626 470592 : if (lphase == 1) then
1627 235296 : mm = numptr_amode(m)
1628 : else
1629 235296 : mm = numptrcw_amode(m)
1630 : endif
1631 2235312 : else if (lspec <= nspec_amode(m)) then ! non-water mass
1632 1764720 : if (lphase == 1) then
1633 882360 : mm = lmassptr_amode(lspec,m)
1634 : else
1635 882360 : mm = lmassptrcw_amode(lspec,m)
1636 : endif
1637 : endif
1638 3176496 : if (lphase == 2) then
1639 1352952 : fldcw => qqcw_get_field(pbuf, mm,lchnk)
1640 2102330088 : fldcw(:ncol,:) = fldcw(:ncol,:) + dcondt_resusp3d(mm,:ncol,:)*dt
1641 : end if
1642 : end do ! loop over number + chem constituents + water
1643 : end do ! lphase
1644 : end do ! m aerosol modes
1645 : end if
1646 :
1647 58824 : call t_stopf('ma_convproc')
1648 : endif
1649 :
1650 : ! if the user has specified prescribed aerosol dep fluxes then
1651 : ! do not set cam_out dep fluxes according to the prognostic aerosols
1652 58824 : if (.not. aerodep_flx_prescribed()) then
1653 58824 : call set_srf_wetdep(aerdepwetis, aerdepwetcw, cam_out)
1654 : endif
1655 :
1656 117648 : endsubroutine aero_model_wetdep
1657 :
1658 : !-------------------------------------------------------------------------
1659 : ! provides wet tropospheric aerosol surface area info for modal aerosols
1660 : ! called from mo_usrrxt
1661 : !-------------------------------------------------------------------------
1662 0 : subroutine aero_model_surfarea( &
1663 0 : mmr, radmean, relhum, pmid, temp, strato_sad, sulfate, rho, ltrop, &
1664 0 : dlat, het1_ndx, pbuf, ncol, sfc, dm_aer, sad_trop, reff_trop )
1665 :
1666 : ! dummy args
1667 : real(r8), intent(in) :: pmid(:,:)
1668 : real(r8), intent(in) :: temp(:,:)
1669 : real(r8), intent(in) :: mmr(:,:,:)
1670 : real(r8), intent(in) :: radmean ! mean radii in cm
1671 : real(r8), intent(in) :: strato_sad(:,:)
1672 : integer, intent(in) :: ncol
1673 : integer, intent(in) :: ltrop(:)
1674 : real(r8), intent(in) :: dlat(:) ! degrees latitude
1675 : integer, intent(in) :: het1_ndx
1676 : real(r8), intent(in) :: relhum(:,:)
1677 : real(r8), intent(in) :: rho(:,:) ! total atm density (/cm^3)
1678 : real(r8), intent(in) :: sulfate(:,:)
1679 : type(physics_buffer_desc), pointer :: pbuf(:)
1680 :
1681 : real(r8), intent(inout) :: sfc(:,:,:)
1682 : real(r8), intent(inout) :: dm_aer(:,:,:)
1683 : real(r8), intent(inout) :: sad_trop(:,:)
1684 : real(r8), intent(out) :: reff_trop(:,:)
1685 :
1686 : ! local vars
1687 0 : real(r8), pointer, dimension(:,:,:) :: dgnumwet
1688 0 : integer :: beglev(ncol)
1689 0 : integer :: endlev(ncol)
1690 : integer :: i,k
1691 :
1692 0 : call pbuf_get_field(pbuf, dgnumwet_idx, dgnumwet )
1693 :
1694 0 : beglev(:ncol)=ltrop(:ncol)+1
1695 0 : endlev(:ncol)=pver
1696 0 : call surf_area_dens( ncol, mmr, pmid, temp, dgnumwet, beglev, endlev, sad_trop, reff_trop, sfc=sfc )
1697 :
1698 0 : do i = 1,ncol
1699 0 : do k = ltrop(i)+1,pver
1700 0 : dm_aer(i,k,:) = dgnumwet(i,k,:) * 1.e2_r8 ! convert m to cm
1701 : enddo
1702 : enddo
1703 :
1704 58824 : end subroutine aero_model_surfarea
1705 :
1706 : !-------------------------------------------------------------------------
1707 : ! provides WET stratospheric aerosol surface area info for modal aerosols
1708 : ! if modal_strat_sulfate = TRUE -- called from mo_gas_phase_chemdr
1709 : !-------------------------------------------------------------------------
1710 0 : subroutine aero_model_strat_surfarea( ncol, mmr, pmid, temp, ltrop, pbuf, strato_sad, reff_strat )
1711 :
1712 : ! dummy args
1713 : integer, intent(in) :: ncol
1714 : real(r8), intent(in) :: mmr(:,:,:)
1715 : real(r8), intent(in) :: pmid(:,:)
1716 : real(r8), intent(in) :: temp(:,:)
1717 : integer, intent(in) :: ltrop(:) ! tropopause level indices
1718 : type(physics_buffer_desc), pointer :: pbuf(:)
1719 : real(r8), intent(out) :: strato_sad(:,:)
1720 : real(r8), intent(out) :: reff_strat(:,:)
1721 :
1722 : ! local vars
1723 0 : real(r8), pointer, dimension(:,:,:) :: dgnumwet
1724 0 : integer :: beglev(ncol)
1725 0 : integer :: endlev(ncol)
1726 :
1727 0 : reff_strat = 0._r8
1728 0 : strato_sad = 0._r8
1729 :
1730 0 : if (.not.modal_strat_sulfate) return
1731 :
1732 0 : call pbuf_get_field(pbuf, dgnumwet_idx, dgnumwet )
1733 :
1734 0 : beglev(:ncol)=top_lev
1735 0 : endlev(:ncol)=ltrop(:ncol)
1736 0 : call surf_area_dens( ncol, mmr, pmid, temp, dgnumwet, beglev, endlev, strato_sad, reff_strat )
1737 :
1738 0 : end subroutine aero_model_strat_surfarea
1739 :
1740 : !=============================================================================
1741 : !=============================================================================
1742 58824 : subroutine aero_model_gasaerexch( loffset, ncol, lchnk, troplev, delt, reaction_rates, &
1743 58824 : tfld, pmid, pdel, mbar, relhum, &
1744 117648 : zm, qh2o, cwat, cldfr, cldnum, &
1745 117648 : airdens, invariants, del_h2so4_gasprod, &
1746 58824 : vmr0, vmr, pbuf )
1747 :
1748 : use time_manager, only : get_nstep
1749 : use modal_aero_coag, only : modal_aero_coag_sub
1750 : use modal_aero_gasaerexch, only : modal_aero_gasaerexch_sub
1751 : use modal_aero_newnuc, only : modal_aero_newnuc_sub
1752 : use modal_aero_data, only : cnst_name_cw, qqcw_get_field
1753 :
1754 : !-----------------------------------------------------------------------
1755 : ! ... dummy arguments
1756 : !-----------------------------------------------------------------------
1757 : integer, intent(in) :: loffset ! offset applied to modal aero "pointers"
1758 : integer, intent(in) :: ncol ! number columns in chunk
1759 : integer, intent(in) :: lchnk ! chunk index
1760 : integer, intent(in) :: troplev(pcols)
1761 : real(r8), intent(in) :: delt ! time step size (sec)
1762 : real(r8), intent(in) :: reaction_rates(:,:,:) ! reaction rates
1763 : real(r8), intent(in) :: tfld(:,:) ! temperature (K)
1764 : real(r8), intent(in) :: pmid(:,:) ! pressure at model levels (Pa)
1765 : real(r8), intent(in) :: pdel(:,:) ! pressure thickness of levels (Pa)
1766 : real(r8), intent(in) :: mbar(:,:) ! mean wet atmospheric mass ( amu )
1767 : real(r8), intent(in) :: relhum(:,:) ! relative humidity
1768 : real(r8), intent(in) :: airdens(:,:) ! total atms density (molec/cm**3)
1769 : real(r8), intent(in) :: invariants(:,:,:)
1770 : real(r8), intent(in) :: del_h2so4_gasprod(:,:)
1771 : real(r8), intent(in) :: zm(:,:)
1772 : real(r8), intent(in) :: qh2o(:,:)
1773 : real(r8), intent(in) :: cwat(:,:) ! cloud liquid water content (kg/kg)
1774 : real(r8), intent(in) :: cldfr(:,:)
1775 : real(r8), intent(in) :: cldnum(:,:) ! droplet number concentration (#/kg)
1776 : real(r8), intent(in) :: vmr0(:,:,:) ! initial mixing ratios (before gas-phase chem changes)
1777 : real(r8), intent(inout) :: vmr(:,:,:) ! mixing ratios ( vmr )
1778 :
1779 : type(physics_buffer_desc), pointer :: pbuf(:)
1780 :
1781 : ! local vars
1782 :
1783 : integer :: n, m
1784 : integer :: i,k,l
1785 : integer :: nstep
1786 :
1787 117648 : real(r8) :: del_h2so4_aeruptk(ncol,pver)
1788 :
1789 58824 : real(r8), pointer :: dgnum(:,:,:), dgnumwet(:,:,:), wetdens(:,:,:)
1790 58824 : real(r8), pointer :: pblh(:) ! pbl height (m)
1791 :
1792 117648 : real(r8), dimension(ncol) :: wrk
1793 : character(len=32) :: name
1794 117648 : real(r8) :: dvmrcwdt(ncol,pver,gas_pcnst)
1795 117648 : real(r8) :: dvmrdt(ncol,pver,gas_pcnst)
1796 117648 : real(r8) :: vmrcw(ncol,pver,gas_pcnst) ! cloud-borne aerosol (vmr)
1797 :
1798 117648 : real(r8) :: aqso4(ncol,ntot_amode) ! aqueous phase chemistry
1799 117648 : real(r8) :: aqh2so4(ncol,ntot_amode) ! aqueous phase chemistry
1800 117648 : real(r8) :: aqso4_h2o2(ncol) ! SO4 aqueous phase chemistry due to H2O2
1801 117648 : real(r8) :: aqso4_o3(ncol) ! SO4 aqueous phase chemistry due to O3
1802 117648 : real(r8) :: xphlwc(ncol,pver) ! pH value multiplied by lwc
1803 117648 : real(r8) :: nh3_beg(ncol,pver)
1804 58824 : real(r8), pointer :: fldcw(:,:)
1805 58824 : real(r8), pointer :: sulfeq(:,:,:)
1806 :
1807 : logical :: is_spcam_m2005
1808 : !
1809 : ! ... initialize nh3
1810 : !
1811 58824 : if ( nh3_ndx > 0 ) then
1812 0 : nh3_beg = vmr(1:ncol,:,nh3_ndx)
1813 : end if
1814 : !
1815 58824 : is_spcam_m2005 = cam_physpkg_is('spcam_m2005')
1816 :
1817 58824 : call pbuf_get_field(pbuf, dgnum_idx, dgnum)
1818 58824 : call pbuf_get_field(pbuf, dgnumwet_idx, dgnumwet )
1819 58824 : call pbuf_get_field(pbuf, wetdens_ap_idx, wetdens )
1820 58824 : call pbuf_get_field(pbuf, pblh_idx, pblh)
1821 :
1822 294120 : do n=1,ntot_amode
1823 235296 : call outfld(dgnum_name(n), dgnum(1:ncol,1:pver,n), ncol, lchnk )
1824 294120 : call outfld(dgnumwet_name(n), dgnumwet(1:ncol,1:pver,n), ncol, lchnk )
1825 : end do
1826 :
1827 : ! do gas-aerosol exchange (h2so4, msa, nh3 condensation)
1828 :
1829 58824 : nstep = get_nstep()
1830 :
1831 : ! calculate tendency due to gas phase chemistry and processes
1832 2833634160 : dvmrdt(:ncol,:,:) = (vmr(:ncol,:,:) - vmr0(:ncol,:,:)) / delt
1833 1882368 : do m = 1, gas_pcnst
1834 30448944 : wrk(:) = 0._r8
1835 171413136 : do k = 1,pver
1836 2833575336 : wrk(:ncol) = wrk(:ncol) + dvmrdt(:ncol,k,m)*adv_mass(m)/mbar(:ncol,k)*pdel(:ncol,k)/gravit
1837 : end do
1838 1823544 : name = 'GS_'//trim(solsym(m))
1839 1882368 : call outfld( name, wrk(:ncol), ncol, lchnk )
1840 : enddo
1841 :
1842 : !
1843 : ! Aerosol processes ...
1844 : !
1845 58824 : call qqcw2vmr( lchnk, vmrcw, mbar, ncol, loffset, pbuf )
1846 :
1847 58824 : if (.not. is_spcam_m2005) then ! regular CAM
1848 2833634160 : dvmrdt(:ncol,:,:) = vmr(:ncol,:,:)
1849 2833634160 : dvmrcwdt(:ncol,:,:) = vmrcw(:ncol,:,:)
1850 :
1851 : ! aqueous chemistry ...
1852 :
1853 58824 : if( has_sox ) then
1854 : call setsox( &
1855 : ncol, &
1856 : lchnk, &
1857 : loffset, &
1858 : delt, &
1859 : pmid, &
1860 : pdel, &
1861 : tfld, &
1862 : mbar, &
1863 : cwat, &
1864 : cldfr, &
1865 : cldnum, &
1866 : airdens, &
1867 : invariants, &
1868 : vmrcw, &
1869 : vmr, &
1870 : xphlwc, &
1871 : aqso4, &
1872 : aqh2so4, &
1873 : aqso4_h2o2, &
1874 : aqso4_o3 &
1875 58824 : )
1876 :
1877 294120 : do n = 1, ntot_amode
1878 235296 : l = lptr_so4_cw_amode(n)
1879 294120 : if (l > 0) then
1880 176472 : call outfld( trim(cnst_name_cw(l))//'AQSO4', aqso4(:ncol,n), ncol, lchnk)
1881 176472 : call outfld( trim(cnst_name_cw(l))//'AQH2SO4', aqh2so4(:ncol,n), ncol, lchnk)
1882 : end if
1883 : end do
1884 :
1885 58824 : call outfld( 'AQSO4_H2O2', aqso4_h2o2(:ncol), ncol, lchnk)
1886 58824 : call outfld( 'AQSO4_O3', aqso4_o3(:ncol), ncol, lchnk)
1887 58824 : call outfld( 'XPH_LWC', xphlwc(:ncol,:), ncol, lchnk )
1888 :
1889 : endif
1890 :
1891 : ! Tendency due to aqueous chemistry
1892 2833692984 : dvmrdt = (vmr - dvmrdt) / delt
1893 2833634160 : dvmrcwdt = (vmrcw - dvmrcwdt) / delt
1894 1882368 : do m = 1, gas_pcnst
1895 30448944 : wrk(:) = 0._r8
1896 171413136 : do k = 1,pver
1897 2833575336 : wrk(:ncol) = wrk(:ncol) + dvmrdt(:ncol,k,m) * adv_mass(m)/mbar(:ncol,k)*pdel(:ncol,k)/gravit
1898 : end do
1899 1823544 : name = 'AQ_'//trim(solsym(m))
1900 1882368 : call outfld( name, wrk(:ncol), ncol, lchnk )
1901 : enddo
1902 :
1903 : else if (is_spcam_m2005) then ! SPCAM ECPP
1904 : ! when ECPP is used, aqueous chemistry is done in ECPP,
1905 : ! and not updated here.
1906 : ! Minghuai Wang, 2010-02 (Minghuai.Wang@pnl.gov)
1907 : !
1908 0 : dvmrdt = 0.0_r8
1909 0 : dvmrcwdt = 0.0_r8
1910 : endif
1911 :
1912 : ! do gas-aerosol exchange (h2so4, msa, nh3 condensation)
1913 :
1914 58824 : if (ndx_h2so4 > 0) then
1915 91405656 : del_h2so4_aeruptk(1:ncol,:) = vmr(1:ncol,:,ndx_h2so4)
1916 : else
1917 0 : del_h2so4_aeruptk(:,:) = 0.0_r8
1918 : endif
1919 :
1920 58824 : call t_startf('modal_gas-aer_exchng')
1921 :
1922 58824 : if ( sulfeq_idx>0 ) then
1923 0 : call pbuf_get_field( pbuf, sulfeq_idx, sulfeq )
1924 : else
1925 58824 : nullify( sulfeq )
1926 : endif
1927 :
1928 : call modal_aero_gasaerexch_sub( &
1929 : lchnk, ncol, nstep, &
1930 : loffset, delt, &
1931 : tfld, pmid, pdel, &
1932 : qh2o, troplev, &
1933 : vmr, vmrcw, &
1934 : dvmrdt, dvmrcwdt, &
1935 : dgnum, dgnumwet, &
1936 58824 : sulfeq )
1937 :
1938 58824 : if (ndx_h2so4 > 0) then
1939 91405656 : del_h2so4_aeruptk(1:ncol,:) = vmr(1:ncol,:,ndx_h2so4) - del_h2so4_aeruptk(1:ncol,:)
1940 : endif
1941 :
1942 58824 : call t_stopf('modal_gas-aer_exchng')
1943 :
1944 58824 : call t_startf('modal_nucl')
1945 :
1946 : ! do aerosol nucleation (new particle formation)
1947 : call modal_aero_newnuc_sub( &
1948 : lchnk, ncol, nstep, &
1949 : loffset, delt, &
1950 : tfld, pmid, pdel, &
1951 : zm, pblh, &
1952 : qh2o, cldfr, &
1953 : vmr, &
1954 21474864 : del_h2so4_gasprod, del_h2so4_aeruptk )
1955 :
1956 58824 : call t_stopf('modal_nucl')
1957 :
1958 58824 : call t_startf('modal_coag')
1959 :
1960 : ! do aerosol coagulation
1961 : call modal_aero_coag_sub( &
1962 : lchnk, ncol, nstep, &
1963 : loffset, delt, &
1964 : tfld, pmid, pdel, &
1965 : vmr, &
1966 : dgnum, dgnumwet, &
1967 58824 : wetdens )
1968 :
1969 58824 : call t_stopf('modal_coag')
1970 :
1971 58824 : call vmr2qqcw( lchnk, vmrcw, mbar, ncol, loffset, pbuf )
1972 :
1973 : ! diagnostics for cloud-borne aerosols...
1974 2470608 : do n = 1,pcnst
1975 2411784 : fldcw => qqcw_get_field(pbuf,n,lchnk,errorhandle=.true.)
1976 2470608 : if(associated(fldcw)) then
1977 1117656 : call outfld( cnst_name_cw(n), fldcw(:,:), pcols, lchnk )
1978 : endif
1979 : end do
1980 : !
1981 : ! ... put missing NH3 into NH4
1982 : !
1983 58824 : if ( nh3_ndx > 0 .and. nh4_ndx > 0 ) then
1984 0 : vmr(1:ncol,:,nh4_ndx) = vmr(1:ncol,:,nh4_ndx) + (nh3_beg-vmr(1:ncol,:,nh3_ndx))
1985 0 : vmr(1:ncol,:,nh4_ndx) = max(0._r8,vmr(1:ncol,:,nh4_ndx))
1986 : end if
1987 :
1988 117648 : end subroutine aero_model_gasaerexch
1989 :
1990 : !=============================================================================
1991 : !=============================================================================
1992 58824 : subroutine aero_model_emissions( state, cam_in )
1993 58824 : use seasalt_model, only: seasalt_emis, seasalt_names, seasalt_indices, seasalt_active,seasalt_nbin
1994 : use dust_model, only: dust_emis, dust_names, dust_indices, dust_active,dust_nbin, dust_nnum
1995 : use physics_types, only: physics_state
1996 :
1997 : ! Arguments:
1998 :
1999 : type(physics_state), intent(in) :: state ! Physics state variables
2000 : type(cam_in_t), intent(inout) :: cam_in ! import state
2001 :
2002 : ! local vars
2003 :
2004 : integer :: lchnk, ncol
2005 : integer :: m, mm
2006 : real(r8) :: soil_erod_tmp(pcols)
2007 : real(r8) :: sflx(pcols) ! accumulate over all bins for output
2008 : real(r8) :: u10cubed(pcols)
2009 : real (r8), parameter :: z0=0.0001_r8 ! m roughness length over oceans--from ocean model
2010 :
2011 58824 : lchnk = state%lchnk
2012 58824 : ncol = state%ncol
2013 :
2014 58824 : if (dust_active) then
2015 :
2016 58824 : call dust_emis( ncol, lchnk, cam_in%dstflx, cam_in%cflx, soil_erod_tmp )
2017 :
2018 : ! some dust emis diagnostics ...
2019 58824 : sflx(:)=0._r8
2020 411768 : do m=1,dust_nbin+dust_nnum
2021 352944 : mm = dust_indices(m)
2022 3123144 : if (m<=dust_nbin) sflx(:ncol)=sflx(:ncol)+cam_in%cflx(:ncol,mm)
2023 411768 : call outfld(trim(dust_names(m))//'SF',cam_in%cflx(:,mm),pcols, lchnk)
2024 : enddo
2025 58824 : call outfld('DSTSFMBL',sflx(:),pcols,lchnk)
2026 58824 : call outfld('LND_MBL',soil_erod_tmp(:),pcols, lchnk )
2027 : endif
2028 :
2029 58824 : if (seasalt_active) then
2030 982224 : u10cubed(:ncol)=sqrt(state%u(:ncol,pver)**2+state%v(:ncol,pver)**2)
2031 : ! move the winds to 10m high from the midpoint of the gridbox:
2032 : ! follows Tie and Seinfeld and Pandis, p.859 with math.
2033 :
2034 982224 : u10cubed(:ncol)=u10cubed(:ncol)*log(10._r8/z0)/log(state%zm(:ncol,pver)/z0)
2035 :
2036 : ! we need them to the 3.41 power, according to Gong et al., 1997:
2037 982224 : u10cubed(:ncol)=u10cubed(:ncol)**3.41_r8
2038 :
2039 58824 : sflx(:)=0._r8
2040 :
2041 58824 : call seasalt_emis( u10cubed, cam_in%sst, cam_in%ocnfrac, ncol, cam_in%cflx )
2042 :
2043 235296 : do m=1,seasalt_nbin
2044 176472 : mm = seasalt_indices(m)
2045 2946672 : sflx(:ncol)=sflx(:ncol)+cam_in%cflx(:ncol,mm)
2046 235296 : call outfld(trim(seasalt_names(m))//'SF',cam_in%cflx(:,mm),pcols,lchnk)
2047 : enddo
2048 58824 : call outfld('SSTSFMBL',sflx(:),pcols,lchnk)
2049 : endif
2050 :
2051 58824 : end subroutine aero_model_emissions
2052 :
2053 : !===============================================================================
2054 : ! private methods
2055 :
2056 :
2057 : !=============================================================================
2058 : !=============================================================================
2059 0 : subroutine surf_area_dens( ncol, mmr, pmid, temp, diam, beglev, endlev, sad, reff, sfc )
2060 58824 : use mo_constants, only : pi
2061 : use modal_aero_data, only : nspec_amode, alnsg_amode
2062 :
2063 : ! dummy args
2064 : integer, intent(in) :: ncol
2065 : real(r8), intent(in) :: mmr(:,:,:)
2066 : real(r8), intent(in) :: pmid(:,:)
2067 : real(r8), intent(in) :: temp(:,:)
2068 : real(r8), intent(in) :: diam(:,:,:)
2069 : integer, intent(in) :: beglev(:)
2070 : integer, intent(in) :: endlev(:)
2071 : real(r8), intent(out) :: sad(:,:)
2072 : real(r8), intent(out) :: reff(:,:)
2073 : real(r8),optional, intent(out) :: sfc(:,:,:)
2074 :
2075 : ! local vars
2076 0 : real(r8) :: sad_mode(pcols,pver,ntot_amode),radeff(pcols,pver)
2077 0 : real(r8) :: vol(pcols,pver),vol_mode(pcols,pver,ntot_amode)
2078 : real(r8) :: rho_air
2079 : integer :: i,k,l,m
2080 : real(r8) :: chm_mass, tot_mass
2081 :
2082 : !
2083 : ! Compute surface aero for each mode.
2084 : ! Total over all modes as the surface area for chemical reactions.
2085 : !
2086 :
2087 0 : sad = 0._r8
2088 0 : sad_mode = 0._r8
2089 0 : vol = 0._r8
2090 0 : vol_mode = 0._r8
2091 0 : reff = 0._r8
2092 :
2093 0 : do i = 1,ncol
2094 0 : do k = beglev(i),endlev(i)
2095 0 : rho_air = pmid(i,k)/(temp(i,k)*287.04_r8)
2096 0 : do l=1,ntot_amode
2097 : !
2098 : ! compute a mass weighting of the number
2099 : !
2100 0 : tot_mass = 0._r8
2101 0 : chm_mass = 0._r8
2102 0 : do m=1,nspec_amode(l)
2103 0 : if ( index_tot_mass(l,m) > 0 ) &
2104 0 : tot_mass = tot_mass + mmr(i,k,index_tot_mass(l,m))
2105 0 : if ( index_chm_mass(l,m) > 0 ) &
2106 0 : chm_mass = chm_mass + mmr(i,k,index_chm_mass(l,m))
2107 : end do
2108 0 : if ( tot_mass > 0._r8 ) then
2109 : ! surface area density
2110 0 : sad_mode(i,k,l) = chm_mass/tot_mass &
2111 0 : * mmr(i,k,num_idx(l))*rho_air*pi*diam(i,k,l)**2._r8 &
2112 0 : * exp(2._r8*alnsg_amode(l)**2._r8) ! m^2/m^3
2113 0 : sad_mode(i,k,l) = 1.e-2_r8 * sad_mode(i,k,l) ! cm^2/cm^3
2114 :
2115 : ! volume calculation, for use in effective radius calculation
2116 : vol_mode(i,k,l) = chm_mass/tot_mass &
2117 0 : * mmr(i,k,num_idx(l))*rho_air*pi/6._r8*diam(i,k,l)**3._r8 &
2118 0 : * exp(4.5_r8*alnsg_amode(l)**2._r8) ! m^3/m^3 = cm^3/cm^3
2119 : else
2120 0 : sad_mode(i,k,l) = 0._r8
2121 0 : vol_mode(i,k,l) = 0._r8
2122 : end if
2123 : end do
2124 0 : sad(i,k) = sum(sad_mode(i,k,:))
2125 0 : vol(i,k) = sum(vol_mode(i,k,:))
2126 0 : reff(i,k) = 3._r8*vol(i,k)/sad(i,k)
2127 :
2128 : enddo
2129 : enddo
2130 :
2131 0 : if (present(sfc)) then
2132 0 : sfc(:,:,:) = sad_mode(:,:,:)
2133 : endif
2134 :
2135 0 : end subroutine surf_area_dens
2136 :
2137 : !===============================================================================
2138 : !===============================================================================
2139 1536 : subroutine modal_aero_bcscavcoef_init
2140 : !-----------------------------------------------------------------------
2141 : !
2142 : ! Purpose:
2143 : ! Computes lookup table for aerosol impaction/interception scavenging rates
2144 : !
2145 : ! Authors: R. Easter
2146 : !
2147 : !-----------------------------------------------------------------------
2148 :
2149 0 : use shr_kind_mod, only: r8 => shr_kind_r8
2150 : use modal_aero_data
2151 : use cam_abortutils, only: endrun
2152 :
2153 : implicit none
2154 :
2155 :
2156 : ! local variables
2157 : integer nnfit_maxd
2158 : parameter (nnfit_maxd=27)
2159 :
2160 : integer i, jgrow, jdens, jpress, jtemp, mode, nnfit
2161 : integer lunerr
2162 :
2163 : real(r8) dg0, dg0_cgs, press, &
2164 : rhodryaero, rhowetaero, rhowetaero_cgs, rmserr, &
2165 : scavratenum, scavratevol, sigmag, &
2166 : temp, wetdiaratio, wetvolratio
2167 : real(r8) aafitnum(1), xxfitnum(1,nnfit_maxd), yyfitnum(nnfit_maxd)
2168 : real(r8) aafitvol(1), xxfitvol(1,nnfit_maxd), yyfitvol(nnfit_maxd)
2169 :
2170 4608 : allocate(scavimptblnum(nimptblgrow_mind:nimptblgrow_maxd, ntot_amode))
2171 3072 : allocate(scavimptblvol(nimptblgrow_mind:nimptblgrow_maxd, ntot_amode))
2172 :
2173 1536 : lunerr = 6
2174 1536 : dlndg_nimptblgrow = log( 1.25_r8 )
2175 :
2176 7680 : modeloop: do mode = 1, ntot_amode
2177 :
2178 6144 : sigmag = sigmag_amode(mode)
2179 :
2180 6144 : rhodryaero = specdens_amode(1,mode)
2181 :
2182 130560 : growloop: do jgrow = nimptblgrow_mind, nimptblgrow_maxd
2183 :
2184 122880 : wetdiaratio = exp( jgrow*dlndg_nimptblgrow )
2185 122880 : dg0 = dgnum_amode(mode)*wetdiaratio
2186 :
2187 : wetvolratio = exp( jgrow*dlndg_nimptblgrow*3._r8 )
2188 122880 : rhowetaero = 1.0_r8 + (rhodryaero-1.0_r8)/wetvolratio
2189 122880 : rhowetaero = min( rhowetaero, rhodryaero )
2190 :
2191 : !
2192 : ! compute impaction scavenging rates at 1 temp-press pair and save
2193 : !
2194 122880 : nnfit = 0
2195 :
2196 122880 : temp = 273.16_r8
2197 122880 : press = 0.75e6_r8 ! dynes/cm2
2198 122880 : rhowetaero = rhodryaero
2199 :
2200 122880 : dg0_cgs = dg0*1.0e2_r8 ! m to cm
2201 122880 : rhowetaero_cgs = rhowetaero*1.0e-3_r8 ! kg/m3 to g/cm3
2202 : call calc_1_impact_rate( &
2203 : dg0_cgs, sigmag, rhowetaero_cgs, temp, press, &
2204 122880 : scavratenum, scavratevol, lunerr )
2205 :
2206 122880 : nnfit = nnfit + 1
2207 : if (nnfit > nnfit_maxd) then
2208 : write(lunerr,9110)
2209 : call endrun()
2210 : end if
2211 : 9110 format( '*** subr. modal_aero_bcscavcoef_init -- nnfit too big' )
2212 :
2213 122880 : xxfitnum(1,nnfit) = 1._r8
2214 122880 : yyfitnum(nnfit) = log( scavratenum )
2215 :
2216 122880 : xxfitvol(1,nnfit) = 1._r8
2217 122880 : yyfitvol(nnfit) = log( scavratevol )
2218 :
2219 : 5900 continue
2220 :
2221 : !
2222 : ! skip mlinfit stuff because scav table no longer has dependencies on
2223 : ! air temp, air press, and particle wet density
2224 : ! just load the log( scavrate--- ) values
2225 : !
2226 : !!
2227 : !! do linear regression
2228 : !! log(scavrate) = a1 + a2*log(wetdens)
2229 : !!
2230 : ! call mlinft( xxfitnum, yyfitnum, aafitnum, nnfit, 1, 1, rmserr )
2231 : ! call mlinft( xxfitvol, yyfitvol, aafitvol, nnfit, 1, 1, rmserr )
2232 : !
2233 : ! scavimptblnum(jgrow,mode) = aafitnum(1)
2234 : ! scavimptblvol(jgrow,mode) = aafitvol(1)
2235 :
2236 122880 : scavimptblnum(jgrow,mode) = yyfitnum(1)
2237 129024 : scavimptblvol(jgrow,mode) = yyfitvol(1)
2238 :
2239 : enddo growloop
2240 : enddo modeloop
2241 1536 : return
2242 1536 : end subroutine modal_aero_bcscavcoef_init
2243 :
2244 : !===============================================================================
2245 : !===============================================================================
2246 588240 : subroutine modal_aero_depvel_part( ncol, t, pmid, ram1, fv, vlc_dry, vlc_trb, vlc_grv, &
2247 : radius_part, density_part, sig_part, moment, lchnk )
2248 :
2249 : ! calculates surface deposition velocity of particles
2250 : ! L. Zhang, S. Gong, J. Padro, and L. Barrie
2251 : ! A size-seggregated particle dry deposition scheme for an atmospheric aerosol module
2252 : ! Atmospheric Environment, 35, 549-560, 2001.
2253 : !
2254 : ! Authors: X. Liu
2255 :
2256 : !
2257 : ! !USES
2258 : !
2259 1536 : use physconst, only: pi,boltz, gravit, rair
2260 : use mo_drydep, only: n_land_type, fraction_landuse
2261 :
2262 : ! !ARGUMENTS:
2263 : !
2264 : implicit none
2265 : !
2266 : real(r8), intent(in) :: t(pcols,pver) !atm temperature (K)
2267 : real(r8), intent(in) :: pmid(pcols,pver) !atm pressure (Pa)
2268 : real(r8), intent(in) :: fv(pcols) !friction velocity (m/s)
2269 : real(r8), intent(in) :: ram1(pcols) !aerodynamical resistance (s/m)
2270 : real(r8), intent(in) :: radius_part(pcols,pver) ! mean (volume/number) particle radius (m)
2271 : real(r8), intent(in) :: density_part(pcols,pver) ! density of particle material (kg/m3)
2272 : real(r8), intent(in) :: sig_part(pcols,pver) ! geometric standard deviation of particles
2273 : integer, intent(in) :: moment ! moment of size distribution (0 for number, 2 for surface area, 3 for volume)
2274 : integer, intent(in) :: ncol
2275 : integer, intent(in) :: lchnk
2276 :
2277 : real(r8), intent(out) :: vlc_trb(pcols) !Turbulent deposn velocity (m/s)
2278 : real(r8), intent(out) :: vlc_grv(pcols,pver) !grav deposn velocity (m/s)
2279 : real(r8), intent(out) :: vlc_dry(pcols,pver) !dry deposn velocity (m/s)
2280 : !------------------------------------------------------------------------
2281 :
2282 : !------------------------------------------------------------------------
2283 : ! Local Variables
2284 : integer :: m,i,k,ix !indices
2285 : real(r8) :: rho !atm density (kg/m**3)
2286 : real(r8) :: vsc_dyn_atm(pcols,pver) ![kg m-1 s-1] Dynamic viscosity of air
2287 : real(r8) :: vsc_knm_atm(pcols,pver) ![m2 s-1] Kinematic viscosity of atmosphere
2288 : real(r8) :: shm_nbr ![frc] Schmidt number
2289 : real(r8) :: stk_nbr ![frc] Stokes number
2290 : real(r8) :: mfp_atm(pcols,pver) ![m] Mean free path of air
2291 : real(r8) :: dff_aer ![m2 s-1] Brownian diffusivity of particle
2292 : real(r8) :: slp_crc(pcols,pver) ![frc] Slip correction factor
2293 : real(r8) :: rss_trb ![s m-1] Resistance to turbulent deposition
2294 : real(r8) :: rss_lmn ![s m-1] Quasi-laminar layer resistance
2295 : real(r8) :: brownian ! collection efficiency for Browning diffusion
2296 : real(r8) :: impaction ! collection efficiency for impaction
2297 : real(r8) :: interception ! collection efficiency for interception
2298 : real(r8) :: stickfrac ! fraction of particles sticking to surface
2299 : real(r8) :: radius_moment(pcols,pver) ! median radius (m) for moment
2300 : real(r8) :: lnsig ! ln(sig_part)
2301 : real(r8) :: dispersion ! accounts for influence of size dist dispersion on bulk settling velocity
2302 : ! assuming radius_part is number mode radius * exp(1.5 ln(sigma))
2303 :
2304 : integer :: lt
2305 : real(r8) :: lnd_frc
2306 : real(r8) :: wrk1, wrk2, wrk3
2307 :
2308 : ! constants
2309 : real(r8) gamma(11) ! exponent of schmidt number
2310 : ! data gamma/0.54d+00, 0.56d+00, 0.57d+00, 0.54d+00, 0.54d+00, &
2311 : ! 0.56d+00, 0.54d+00, 0.54d+00, 0.54d+00, 0.56d+00, &
2312 : ! 0.50d+00/
2313 : data gamma/0.56e+00_r8, 0.54e+00_r8, 0.54e+00_r8, 0.56e+00_r8, 0.56e+00_r8, &
2314 : 0.56e+00_r8, 0.50e+00_r8, 0.54e+00_r8, 0.54e+00_r8, 0.54e+00_r8, &
2315 : 0.54e+00_r8/
2316 : save gamma
2317 :
2318 : real(r8) alpha(11) ! parameter for impaction
2319 : ! data alpha/50.00d+00, 0.95d+00, 0.80d+00, 1.20d+00, 1.30d+00, &
2320 : ! 0.80d+00, 50.00d+00, 50.00d+00, 2.00d+00, 1.50d+00, &
2321 : ! 100.00d+00/
2322 : data alpha/1.50e+00_r8, 1.20e+00_r8, 1.20e+00_r8, 0.80e+00_r8, 1.00e+00_r8, &
2323 : 0.80e+00_r8, 100.00e+00_r8, 50.00e+00_r8, 2.00e+00_r8, 1.20e+00_r8, &
2324 : 50.00e+00_r8/
2325 : save alpha
2326 :
2327 : real(r8) radius_collector(11) ! radius (m) of surface collectors
2328 : ! data radius_collector/-1.00d+00, 5.10d-03, 3.50d-03, 3.20d-03, 10.00d-03, &
2329 : ! 5.00d-03, -1.00d+00, -1.00d+00, 10.00d-03, 10.00d-03, &
2330 : ! -1.00d+00/
2331 : data radius_collector/10.00e-03_r8, 3.50e-03_r8, 3.50e-03_r8, 5.10e-03_r8, 2.00e-03_r8, &
2332 : 5.00e-03_r8, -1.00e+00_r8, -1.00e+00_r8, 10.00e-03_r8, 3.50e-03_r8, &
2333 : -1.00e+00_r8/
2334 : save radius_collector
2335 :
2336 : integer :: iwet(11) ! flag for wet surface = 1, otherwise = -1
2337 : ! data iwet/1, -1, -1, -1, -1, &
2338 : ! -1, -1, -1, 1, -1, &
2339 : ! 1/
2340 : data iwet/-1, -1, -1, -1, -1, &
2341 : -1, 1, -1, 1, -1, &
2342 : -1/
2343 : save iwet
2344 :
2345 :
2346 588240 : vlc_trb = 0._r8
2347 588240 : vlc_grv = 0._r8
2348 588240 : vlc_dry = 0._r8
2349 :
2350 : !------------------------------------------------------------------------
2351 55294560 : do k=top_lev,pver ! radius_part is not defined above top_lev
2352 914056560 : do i=1,ncol
2353 :
2354 858762000 : lnsig = log(sig_part(i,k))
2355 : ! use a maximum radius of 50 microns when calculating deposition velocity
2356 : radius_moment(i,k) = min(50.0e-6_r8,radius_part(i,k))* &
2357 858762000 : exp((float(moment)-1.5_r8)*lnsig*lnsig)
2358 858762000 : dispersion = exp(2._r8*lnsig*lnsig)
2359 :
2360 858762000 : rho=pmid(i,k)/rair/t(i,k)
2361 :
2362 : ! Quasi-laminar layer resistance: call rss_lmn_get
2363 : ! Size-independent thermokinetic properties
2364 : vsc_dyn_atm(i,k) = 1.72e-5_r8 * ((t(i,k)/273.0_r8)**1.5_r8) * 393.0_r8 / &
2365 858762000 : (t(i,k)+120.0_r8) ![kg m-1 s-1] RoY94 p. 102
2366 : mfp_atm(i,k) = 2.0_r8 * vsc_dyn_atm(i,k) / & ![m] SeP97 p. 455
2367 858762000 : (pmid(i,k)*sqrt(8.0_r8/(pi*rair*t(i,k))))
2368 858762000 : vsc_knm_atm(i,k) = vsc_dyn_atm(i,k) / rho ![m2 s-1] Kinematic viscosity of air
2369 :
2370 : slp_crc(i,k) = 1.0_r8 + mfp_atm(i,k) * &
2371 : (1.257_r8+0.4_r8*exp(-1.1_r8*radius_moment(i,k)/(mfp_atm(i,k)))) / &
2372 858762000 : radius_moment(i,k) ![frc] Slip correction factor SeP97 p. 464
2373 : vlc_grv(i,k) = (4.0_r8/18.0_r8) * radius_moment(i,k)*radius_moment(i,k)*density_part(i,k)* &
2374 858762000 : gravit*slp_crc(i,k) / vsc_dyn_atm(i,k) ![m s-1] Stokes' settling velocity SeP97 p. 466
2375 858762000 : vlc_grv(i,k) = vlc_grv(i,k) * dispersion
2376 :
2377 913468320 : vlc_dry(i,k)=vlc_grv(i,k)
2378 : enddo
2379 : enddo
2380 588240 : k=pver ! only look at bottom level for next part
2381 9822240 : do i=1,ncol
2382 9234000 : dff_aer = boltz * t(i,k) * slp_crc(i,k) / & ![m2 s-1]
2383 9234000 : (6.0_r8*pi*vsc_dyn_atm(i,k)*radius_moment(i,k)) !SeP97 p.474
2384 9234000 : shm_nbr = vsc_knm_atm(i,k) / dff_aer ![frc] SeP97 p.972
2385 :
2386 9234000 : wrk2 = 0._r8
2387 9234000 : wrk3 = 0._r8
2388 110808000 : do lt = 1,n_land_type
2389 101574000 : lnd_frc = fraction_landuse(i,lt,lchnk)
2390 110808000 : if ( lnd_frc /= 0._r8 ) then
2391 22542930 : brownian = shm_nbr**(-gamma(lt))
2392 22542930 : if (radius_collector(lt) > 0.0_r8) then
2393 : ! vegetated surface
2394 10266270 : stk_nbr = vlc_grv(i,k) * fv(i) / (gravit*radius_collector(lt))
2395 10266270 : interception = 2.0_r8*(radius_moment(i,k)/radius_collector(lt))**2.0_r8
2396 : else
2397 : ! non-vegetated surface
2398 12276660 : stk_nbr = vlc_grv(i,k) * fv(i) * fv(i) / (gravit*vsc_knm_atm(i,k)) ![frc] SeP97 p.965
2399 12276660 : interception = 0.0_r8
2400 : endif
2401 22542930 : impaction = (stk_nbr/(alpha(lt)+stk_nbr))**2.0_r8
2402 :
2403 22542930 : if (iwet(lt) > 0) then
2404 : stickfrac = 1.0_r8
2405 : else
2406 13968990 : stickfrac = exp(-sqrt(stk_nbr))
2407 13968990 : if (stickfrac < 1.0e-10_r8) stickfrac = 1.0e-10_r8
2408 : endif
2409 22542930 : rss_lmn = 1.0_r8 / (3.0_r8 * fv(i) * stickfrac * (brownian+interception+impaction))
2410 22542930 : rss_trb = ram1(i) + rss_lmn + ram1(i)*rss_lmn*vlc_grv(i,k)
2411 :
2412 22542930 : wrk1 = 1.0_r8 / rss_trb
2413 22542930 : wrk2 = wrk2 + lnd_frc*( wrk1 )
2414 22542930 : wrk3 = wrk3 + lnd_frc*( wrk1 + vlc_grv(i,k) )
2415 : endif
2416 : enddo ! n_land_type
2417 9234000 : vlc_trb(i) = wrk2
2418 9822240 : vlc_dry(i,k) = wrk3
2419 : enddo !ncol
2420 :
2421 588240 : return
2422 588240 : end subroutine modal_aero_depvel_part
2423 :
2424 : !===============================================================================
2425 470592 : subroutine modal_aero_bcscavcoef_get( m, ncol, isprx, dgn_awet, scavcoefnum, scavcoefvol )
2426 :
2427 588240 : use modal_aero_data
2428 : !-----------------------------------------------------------------------
2429 : implicit none
2430 :
2431 : integer,intent(in) :: m, ncol
2432 : logical,intent(in):: isprx(pcols,pver)
2433 : real(r8), intent(in) :: dgn_awet(pcols,pver,ntot_amode)
2434 : real(r8), intent(out) :: scavcoefnum(pcols,pver), scavcoefvol(pcols,pver)
2435 :
2436 : integer i, k, jgrow
2437 : real(r8) dumdgratio, xgrow, dumfhi, dumflo, scavimpvol, scavimpnum
2438 :
2439 :
2440 22117824 : do k = 1, pver
2441 365622624 : do i = 1, ncol
2442 :
2443 : ! do only if no precip
2444 365387328 : if ( isprx(i,k) ) then
2445 : !
2446 : ! interpolate table values using log of (actual-wet-size)/(base-dry-size)
2447 :
2448 56378516 : dumdgratio = dgn_awet(i,k,m)/dgnum_amode(m)
2449 :
2450 56378516 : if ((dumdgratio >= 0.99_r8) .and. (dumdgratio <= 1.01_r8)) then
2451 458255 : scavimpvol = scavimptblvol(0,m)
2452 458255 : scavimpnum = scavimptblnum(0,m)
2453 : else
2454 55920261 : xgrow = log( dumdgratio ) / dlndg_nimptblgrow
2455 55920261 : jgrow = int( xgrow )
2456 55920261 : if (xgrow < 0._r8) jgrow = jgrow - 1
2457 55920261 : if (jgrow < nimptblgrow_mind) then
2458 : jgrow = nimptblgrow_mind
2459 : xgrow = jgrow
2460 : else
2461 55920220 : jgrow = min( jgrow, nimptblgrow_maxd-1 )
2462 : end if
2463 :
2464 55920261 : dumfhi = xgrow - jgrow
2465 55920261 : dumflo = 1._r8 - dumfhi
2466 :
2467 0 : scavimpvol = dumflo*scavimptblvol(jgrow,m) + &
2468 55920261 : dumfhi*scavimptblvol(jgrow+1,m)
2469 0 : scavimpnum = dumflo*scavimptblnum(jgrow,m) + &
2470 55920261 : dumfhi*scavimptblnum(jgrow+1,m)
2471 :
2472 : end if
2473 :
2474 : ! impaction scavenging removal amount for volume
2475 56378516 : scavcoefvol(i,k) = exp( scavimpvol )
2476 : ! impaction scavenging removal amount to number
2477 56378516 : scavcoefnum(i,k) = exp( scavimpnum )
2478 :
2479 : ! scavcoef = impaction scav rate (1/h) for precip = 1 mm/h
2480 : ! scavcoef = impaction scav rate (1/s) for precip = pfx_inrain
2481 : ! (scavcoef/3600) = impaction scav rate (1/s) for precip = 1 mm/h
2482 : ! (pfx_inrain*3600) = in-rain-area precip rate (mm/h)
2483 : ! impactrate = (scavcoef/3600) * (pfx_inrain*3600)
2484 : else
2485 287126284 : scavcoefvol(i,k) = 0._r8
2486 287126284 : scavcoefnum(i,k) = 0._r8
2487 : end if
2488 :
2489 : end do
2490 : end do
2491 :
2492 235296 : return
2493 235296 : end subroutine modal_aero_bcscavcoef_get
2494 :
2495 : !===============================================================================
2496 122880 : subroutine calc_1_impact_rate( &
2497 : dg0, sigmag, rhoaero, temp, press, &
2498 : scavratenum, scavratevol, lunerr )
2499 : !
2500 : ! routine computes a single impaction scavenging rate
2501 : ! for precipitation rate of 1 mm/h
2502 : !
2503 : ! dg0 = geometric mean diameter of aerosol number size distrib. (cm)
2504 : ! sigmag = geometric standard deviation of size distrib.
2505 : ! rhoaero = density of aerosol particles (g/cm^3)
2506 : ! temp = temperature (K)
2507 : ! press = pressure (dyne/cm^2)
2508 : ! scavratenum = number scavenging rate (1/h)
2509 : ! scavratevol = volume or mass scavenging rate (1/h)
2510 : ! lunerr = logical unit for error message
2511 : !
2512 235296 : use shr_kind_mod, only: r8 => shr_kind_r8
2513 : use mo_constants, only: boltz_cgs, pi, rhowater => rhoh2o_cgs, &
2514 : gravity => gravity_cgs, rgas => rgas_cgs
2515 :
2516 : implicit none
2517 :
2518 : ! subr. parameters
2519 : integer lunerr
2520 : real(r8) dg0, sigmag, rhoaero, temp, press, scavratenum, scavratevol
2521 :
2522 : ! local variables
2523 : integer nrainsvmax
2524 : parameter (nrainsvmax=50)
2525 : real(r8) rrainsv(nrainsvmax), xnumrainsv(nrainsvmax),&
2526 : vfallrainsv(nrainsvmax)
2527 :
2528 : integer naerosvmax
2529 : parameter (naerosvmax=51)
2530 : real(r8) aaerosv(naerosvmax), &
2531 : ynumaerosv(naerosvmax), yvolaerosv(naerosvmax)
2532 :
2533 : integer i, ja, jr, na, nr
2534 : real(r8) a, aerodiffus, aeromass, ag0, airdynvisc, airkinvisc
2535 : real(r8) anumsum, avolsum, cair, chi
2536 : real(r8) d, dr, dum, dumfuchs, dx
2537 : real(r8) ebrown, eimpact, eintercept, etotal, freepath
2538 : real(r8) precip, precipmmhr, precipsum
2539 : real(r8) r, rainsweepout, reynolds, rhi, rhoair, rlo, rnumsum
2540 : real(r8) scavsumnum, scavsumnumbb
2541 : real(r8) scavsumvol, scavsumvolbb
2542 : real(r8) schmidt, sqrtreynolds, sstar, stokes, sx
2543 : real(r8) taurelax, vfall, vfallstp
2544 : real(r8) x, xg0, xg3, xhi, xlo, xmuwaterair
2545 :
2546 :
2547 122880 : rlo = .005_r8
2548 122880 : rhi = .250_r8
2549 122880 : dr = 0.005_r8
2550 122880 : nr = 1 + nint( (rhi-rlo)/dr )
2551 : if (nr > nrainsvmax) then
2552 : write(lunerr,9110)
2553 : call endrun()
2554 : end if
2555 :
2556 : 9110 format( '*** subr. calc_1_impact_rate -- nr > nrainsvmax' )
2557 :
2558 122880 : precipmmhr = 1.0_r8
2559 122880 : precip = precipmmhr/36000._r8
2560 :
2561 122880 : ag0 = dg0/2._r8
2562 122880 : sx = log( sigmag )
2563 122880 : xg0 = log( ag0 )
2564 122880 : xg3 = xg0 + 3._r8*sx*sx
2565 :
2566 122880 : xlo = xg3 - 4._r8*sx
2567 122880 : xhi = xg3 + 4._r8*sx
2568 122880 : dx = 0.2_r8*sx
2569 :
2570 122880 : dx = max( 0.2_r8*sx, 0.01_r8 )
2571 122880 : xlo = xg3 - max( 4._r8*sx, 2._r8*dx )
2572 122880 : xhi = xg3 + max( 4._r8*sx, 2._r8*dx )
2573 :
2574 122880 : na = 1 + nint( (xhi-xlo)/dx )
2575 122880 : if (na > naerosvmax) then
2576 0 : write(lunerr,9120)
2577 0 : call endrun()
2578 : end if
2579 :
2580 : 9120 format( '*** subr. calc_1_impact_rate -- na > naerosvmax' )
2581 :
2582 : ! air molar density
2583 122880 : cair = press/(rgas*temp)
2584 : ! air mass density
2585 122880 : rhoair = 28.966_r8*cair
2586 : ! molecular freepath
2587 122880 : freepath = 2.8052e-10_r8/cair
2588 : ! air dynamic viscosity
2589 : airdynvisc = 1.8325e-4_r8 * (416.16_r8/(temp+120._r8)) * &
2590 122880 : ((temp/296.16_r8)**1.5_r8)
2591 : ! air kinemaic viscosity
2592 122880 : airkinvisc = airdynvisc/rhoair
2593 : ! ratio of water viscosity to air viscosity (from Slinn)
2594 122880 : xmuwaterair = 60.0_r8
2595 :
2596 : !
2597 : ! compute rain drop number concentrations
2598 : ! rrainsv = raindrop radius (cm)
2599 : ! xnumrainsv = raindrop number concentration (#/cm^3)
2600 : ! (number in the bin, not number density)
2601 : ! vfallrainsv = fall velocity (cm/s)
2602 : !
2603 122880 : precipsum = 0._r8
2604 6266880 : do i = 1, nr
2605 6144000 : r = rlo + (i-1)*dr
2606 6144000 : rrainsv(i) = r
2607 6144000 : xnumrainsv(i) = exp( -r/2.7e-2_r8 )
2608 :
2609 6144000 : d = 2._r8*r
2610 6144000 : if (d <= 0.007_r8) then
2611 0 : vfallstp = 2.88e5_r8 * d**2._r8
2612 6144000 : else if (d <= 0.025_r8) then
2613 245760 : vfallstp = 2.8008e4_r8 * d**1.528_r8
2614 5898240 : else if (d <= 0.1_r8) then
2615 983040 : vfallstp = 4104.9_r8 * d**1.008_r8
2616 4915200 : else if (d <= 0.25_r8) then
2617 1843200 : vfallstp = 1812.1_r8 * d**0.638_r8
2618 : else
2619 3072000 : vfallstp = 1069.8_r8 * d**0.235_r8
2620 : end if
2621 :
2622 6144000 : vfall = vfallstp * sqrt(1.204e-3_r8/rhoair)
2623 6144000 : vfallrainsv(i) = vfall
2624 6266880 : precipsum = precipsum + vfall*(r**3)*xnumrainsv(i)
2625 : end do
2626 122880 : precipsum = precipsum*pi*1.333333_r8
2627 :
2628 122880 : rnumsum = 0._r8
2629 6266880 : do i = 1, nr
2630 6144000 : xnumrainsv(i) = xnumrainsv(i)*(precip/precipsum)
2631 6266880 : rnumsum = rnumsum + xnumrainsv(i)
2632 : end do
2633 :
2634 : !
2635 : ! compute aerosol concentrations
2636 : ! aaerosv = particle radius (cm)
2637 : ! fnumaerosv = fraction of total number in the bin (--)
2638 : ! fvolaerosv = fraction of total volume in the bin (--)
2639 : !
2640 : anumsum = 0._r8
2641 : avolsum = 0._r8
2642 5160960 : do i = 1, na
2643 5038080 : x = xlo + (i-1)*dx
2644 5038080 : a = exp( x )
2645 5038080 : aaerosv(i) = a
2646 5038080 : dum = (x - xg0)/sx
2647 5038080 : ynumaerosv(i) = exp( -0.5_r8*dum*dum )
2648 5038080 : yvolaerosv(i) = ynumaerosv(i)*1.3333_r8*pi*a*a*a
2649 5038080 : anumsum = anumsum + ynumaerosv(i)
2650 5160960 : avolsum = avolsum + yvolaerosv(i)
2651 : end do
2652 :
2653 5160960 : do i = 1, na
2654 5038080 : ynumaerosv(i) = ynumaerosv(i)/anumsum
2655 5160960 : yvolaerosv(i) = yvolaerosv(i)/avolsum
2656 : end do
2657 :
2658 :
2659 : !
2660 : ! compute scavenging
2661 : !
2662 : scavsumnum = 0._r8
2663 : scavsumvol = 0._r8
2664 : !
2665 : ! outer loop for rain drop radius
2666 : !
2667 6266880 : jr_loop: do jr = 1, nr
2668 :
2669 6144000 : r = rrainsv(jr)
2670 6144000 : vfall = vfallrainsv(jr)
2671 :
2672 6144000 : reynolds = r * vfall / airkinvisc
2673 6144000 : sqrtreynolds = sqrt( reynolds )
2674 :
2675 : !
2676 : ! inner loop for aerosol particle radius
2677 : !
2678 6144000 : scavsumnumbb = 0._r8
2679 6144000 : scavsumvolbb = 0._r8
2680 :
2681 258048000 : ja_loop: do ja = 1, na
2682 :
2683 251904000 : a = aaerosv(ja)
2684 :
2685 251904000 : chi = a/r
2686 :
2687 251904000 : dum = freepath/a
2688 251904000 : dumfuchs = 1._r8 + 1.246_r8*dum + 0.42_r8*dum*exp(-0.87_r8/dum)
2689 251904000 : taurelax = 2._r8*rhoaero*a*a*dumfuchs/(9._r8*rhoair*airkinvisc)
2690 :
2691 251904000 : aeromass = 4._r8*pi*a*a*a*rhoaero/3._r8
2692 251904000 : aerodiffus = boltz_cgs*temp*taurelax/aeromass
2693 :
2694 251904000 : schmidt = airkinvisc/aerodiffus
2695 251904000 : stokes = vfall*taurelax/r
2696 :
2697 : ebrown = 4._r8*(1._r8 + 0.4_r8*sqrtreynolds*(schmidt**0.3333333_r8)) / &
2698 251904000 : (reynolds*schmidt)
2699 :
2700 : dum = (1._r8 + 2._r8*xmuwaterair*chi) / &
2701 251904000 : (1._r8 + xmuwaterair/sqrtreynolds)
2702 251904000 : eintercept = 4._r8*chi*(chi + dum)
2703 :
2704 251904000 : dum = log( 1._r8 + reynolds )
2705 251904000 : sstar = (1.2_r8 + dum/12._r8) / (1._r8 + dum)
2706 251904000 : eimpact = 0._r8
2707 251904000 : if (stokes > sstar) then
2708 61905408 : dum = stokes - sstar
2709 61905408 : eimpact = (dum/(dum+0.6666667_r8)) ** 1.5_r8
2710 : end if
2711 :
2712 251904000 : etotal = ebrown + eintercept + eimpact
2713 251904000 : etotal = min( etotal, 1.0_r8 )
2714 :
2715 251904000 : rainsweepout = xnumrainsv(jr)*4._r8*pi*r*r*vfall
2716 :
2717 251904000 : scavsumnumbb = scavsumnumbb + rainsweepout*etotal*ynumaerosv(ja)
2718 258048000 : scavsumvolbb = scavsumvolbb + rainsweepout*etotal*yvolaerosv(ja)
2719 :
2720 : enddo ja_loop
2721 :
2722 6144000 : scavsumnum = scavsumnum + scavsumnumbb
2723 6266880 : scavsumvol = scavsumvol + scavsumvolbb
2724 :
2725 : enddo jr_loop
2726 :
2727 122880 : scavratenum = scavsumnum*3600._r8
2728 122880 : scavratevol = scavsumvol*3600._r8
2729 :
2730 122880 : return
2731 : end subroutine calc_1_impact_rate
2732 :
2733 : !=============================================================================
2734 : !=============================================================================
2735 117648 : subroutine qqcw2vmr(lchnk, vmr, mbar, ncol, im, pbuf)
2736 : use modal_aero_data, only : qqcw_get_field
2737 : use physics_buffer, only : physics_buffer_desc
2738 : !-----------------------------------------------------------------
2739 : ! ... Xfrom from mass to volume mixing ratio
2740 : !-----------------------------------------------------------------
2741 :
2742 : use chem_mods, only : adv_mass, gas_pcnst
2743 :
2744 : implicit none
2745 :
2746 : !-----------------------------------------------------------------
2747 : ! ... Dummy args
2748 : !-----------------------------------------------------------------
2749 : integer, intent(in) :: lchnk, ncol, im
2750 : real(r8), intent(in) :: mbar(ncol,pver)
2751 : real(r8), intent(inout) :: vmr(ncol,pver,gas_pcnst)
2752 : type(physics_buffer_desc), pointer :: pbuf(:)
2753 :
2754 : !-----------------------------------------------------------------
2755 : ! ... Local variables
2756 : !-----------------------------------------------------------------
2757 : integer :: k, m
2758 58824 : real(r8), pointer :: fldcw(:,:)
2759 :
2760 1882368 : do m=1,gas_pcnst
2761 1882368 : if( adv_mass(m) /= 0._r8 ) then
2762 1823544 : fldcw => qqcw_get_field(pbuf, m+im,lchnk,errorhandle=.true.)
2763 1823544 : if(associated(fldcw)) then
2764 105059664 : do k=1,pver
2765 1736707464 : vmr(:ncol,k,m) = mbar(:ncol,k) * fldcw(:ncol,k) / adv_mass(m)
2766 : end do
2767 : else
2768 1096867872 : vmr(:,:,m) = 0.0_r8
2769 : end if
2770 : end if
2771 : end do
2772 117648 : end subroutine qqcw2vmr
2773 :
2774 :
2775 : !=============================================================================
2776 : !=============================================================================
2777 117648 : subroutine vmr2qqcw( lchnk, vmr, mbar, ncol, im, pbuf )
2778 : !-----------------------------------------------------------------
2779 : ! ... Xfrom from volume to mass mixing ratio
2780 : !-----------------------------------------------------------------
2781 :
2782 58824 : use m_spc_id
2783 : use chem_mods, only : adv_mass, gas_pcnst
2784 : use modal_aero_data, only : qqcw_get_field
2785 : use physics_buffer, only : physics_buffer_desc
2786 :
2787 : implicit none
2788 :
2789 : !-----------------------------------------------------------------
2790 : ! ... Dummy args
2791 : !-----------------------------------------------------------------
2792 : integer, intent(in) :: lchnk, ncol, im
2793 : real(r8), intent(in) :: mbar(ncol,pver)
2794 : real(r8), intent(in) :: vmr(ncol,pver,gas_pcnst)
2795 : type(physics_buffer_desc), pointer :: pbuf(:)
2796 :
2797 : !-----------------------------------------------------------------
2798 : ! ... Local variables
2799 : !-----------------------------------------------------------------
2800 : integer :: k, m
2801 58824 : real(r8), pointer :: fldcw(:,:)
2802 : !-----------------------------------------------------------------
2803 : ! ... The non-group species
2804 : !-----------------------------------------------------------------
2805 1882368 : do m = 1,gas_pcnst
2806 1823544 : fldcw => qqcw_get_field(pbuf, m+im,lchnk,errorhandle=.true.)
2807 1882368 : if( adv_mass(m) /= 0._r8 .and. associated(fldcw)) then
2808 105059664 : do k = 1,pver
2809 1736707464 : fldcw(:ncol,k) = adv_mass(m) * vmr(:ncol,k,m) / mbar(:ncol,k)
2810 : end do
2811 : end if
2812 : end do
2813 :
2814 117648 : end subroutine vmr2qqcw
2815 :
2816 0 : function get_dlndg_nimptblgrow() result (dlndg_nimptblgrow_ret)
2817 : real(r8) :: dlndg_nimptblgrow_ret
2818 0 : dlndg_nimptblgrow_ret = dlndg_nimptblgrow
2819 58824 : end function get_dlndg_nimptblgrow
2820 :
2821 0 : function get_scavimptblvol() result (scavimptblvol_ret)
2822 : real(r8) :: scavimptblvol_ret(nimptblgrow_mind:nimptblgrow_maxd, ntot_amode)
2823 0 : scavimptblvol_ret = scavimptblvol
2824 0 : end function get_scavimptblvol
2825 :
2826 0 : function get_scavimptblnum() result (scavimptblnum_ret)
2827 : real(r8) :: scavimptblnum_ret(nimptblgrow_mind:nimptblgrow_maxd, ntot_amode)
2828 0 : scavimptblnum_ret = scavimptblnum
2829 0 : end function get_scavimptblnum
2830 :
2831 : end module aero_model
|
