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 : use modal_aerosol_properties_mod, only: modal_aerosol_properties
33 :
34 : implicit none
35 : private
36 :
37 : public :: aero_model_readnl
38 : public :: aero_model_register
39 : public :: aero_model_init
40 : public :: aero_model_gasaerexch ! create, grow, change, and shrink aerosols.
41 : public :: aero_model_drydep ! aerosol dry deposition and sediment
42 : public :: aero_model_wetdep ! aerosol wet removal
43 : public :: aero_model_emissions ! aerosol emissions
44 : public :: aero_model_surfarea ! tropopspheric aerosol wet surface area for chemistry
45 : public :: aero_model_strat_surfarea ! stratospheric aerosol wet surface area for chemistry
46 :
47 : public :: calc_1_impact_rate
48 : public :: nimptblgrow_mind, nimptblgrow_maxd
49 :
50 : ! Accessor functions
51 : public :: get_scavimptblvol, get_scavimptblnum, get_dlndg_nimptblgrow
52 :
53 : ! Misc private data
54 :
55 : ! number of modes
56 : integer :: nmodes
57 : integer :: pblh_idx = 0
58 : integer :: dgnum_idx = 0
59 : integer :: dgnumwet_idx = 0
60 : integer :: rate1_cw2pr_st_idx = 0
61 :
62 : integer :: wetdens_ap_idx = 0
63 : integer :: qaerwat_idx = 0
64 :
65 : integer :: fracis_idx = 0
66 : integer :: prain_idx = 0
67 : integer :: rprddp_idx = 0
68 : integer :: rprdsh_idx = 0
69 : integer :: nevapr_shcu_idx = 0
70 : integer :: nevapr_dpcu_idx = 0
71 :
72 : integer :: sulfeq_idx = -1
73 :
74 : integer :: nh3_ndx = 0
75 : integer :: nh4_ndx = 0
76 :
77 : ! variables for table lookup of aerosol impaction/interception scavenging rates
78 : integer, parameter :: nimptblgrow_mind=-7, nimptblgrow_maxd=12
79 : real(r8) :: dlndg_nimptblgrow
80 : real(r8),allocatable :: scavimptblnum(:,:)
81 : real(r8),allocatable :: scavimptblvol(:,:)
82 :
83 : ! for surf_area_dens
84 : integer,allocatable :: num_idx(:)
85 : integer,allocatable :: index_tot_mass(:,:)
86 : integer,allocatable :: index_chm_mass(:,:)
87 :
88 : integer :: ndx_h2so4
89 : character(len=fieldname_len), allocatable :: dgnum_name(:), dgnumwet_name(:)
90 :
91 : ! Namelist variables
92 : character(len=16) :: drydep_list(pcnst) = ' '
93 : real(r8) :: seasalt_emis_scale
94 :
95 : integer :: ndrydep = 0
96 : integer,allocatable :: drydep_indices(:)
97 : logical :: drydep_lq(pcnst)
98 :
99 : logical :: modal_accum_coarse_exch = .false.
100 :
101 : type(modal_aerosol_properties), pointer :: aero_props=>null()
102 :
103 : contains
104 :
105 : !=============================================================================
106 : ! reads aerosol namelist options
107 : !=============================================================================
108 2304 : subroutine aero_model_readnl(nlfile)
109 :
110 : use namelist_utils, only: find_group_name
111 : use units, only: getunit, freeunit
112 : use mpishorthand
113 : use aero_wetdep_cam, only: aero_wetdep_readnl
114 : use dust_model, only: dust_readnl
115 :
116 : character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input
117 :
118 : ! Local variables
119 : integer :: unitn, ierr
120 : character(len=*), parameter :: subname = 'aero_model_readnl'
121 :
122 : ! Namelist variables
123 : character(len=16) :: aer_drydep_list(pcnst) = ' '
124 :
125 : namelist /aerosol_nl/ aer_drydep_list, modal_strat_sulfate, modal_accum_coarse_exch, seasalt_emis_scale
126 :
127 : !-----------------------------------------------------------------------------
128 :
129 : ! Read namelist
130 2304 : if (masterproc) then
131 3 : unitn = getunit()
132 3 : open( unitn, file=trim(nlfile), status='old' )
133 3 : call find_group_name(unitn, 'aerosol_nl', status=ierr)
134 3 : if (ierr == 0) then
135 3 : read(unitn, aerosol_nl, iostat=ierr)
136 3 : if (ierr /= 0) then
137 0 : call endrun(subname // ':: ERROR reading namelist')
138 : end if
139 : end if
140 3 : close(unitn)
141 3 : call freeunit(unitn)
142 : end if
143 :
144 : #ifdef SPMD
145 : ! Broadcast namelist variables
146 2304 : call mpibcast(aer_drydep_list, len(aer_drydep_list(1))*pcnst, mpichar, 0, mpicom)
147 2304 : call mpibcast(modal_strat_sulfate, 1, mpilog, 0, mpicom)
148 2304 : call mpibcast(seasalt_emis_scale, 1, mpir8, 0, mpicom)
149 2304 : call mpibcast(modal_accum_coarse_exch, 1, mpilog, 0, mpicom)
150 : #endif
151 :
152 96768 : drydep_list = aer_drydep_list
153 :
154 2304 : call aero_wetdep_readnl(nlfile)
155 2304 : call dust_readnl(nlfile)
156 :
157 2304 : end subroutine aero_model_readnl
158 :
159 : !=============================================================================
160 : !=============================================================================
161 2304 : subroutine aero_model_register()
162 2304 : use modal_aero_data, only: modal_aero_data_reg
163 :
164 2304 : call modal_aero_data_reg()
165 :
166 2304 : end subroutine aero_model_register
167 :
168 : !=============================================================================
169 : !=============================================================================
170 2304 : subroutine aero_model_init( pbuf2d )
171 :
172 : use mo_chem_utls, only: get_inv_ndx
173 : use cam_history, only: addfld, add_default, horiz_only
174 2304 : use mo_chem_utls, only: get_spc_ndx
175 : use modal_aero_data, only: cnst_name_cw
176 : use modal_aero_data, only: modal_aero_data_init
177 : use rad_constituents,only: rad_cnst_get_info
178 : use dust_model, only: dust_init, dust_names, dust_active, dust_nbin, dust_nnum
179 : use seasalt_model, only: seasalt_init, seasalt_names, seasalt_active,seasalt_nbin
180 : use aer_drydep_mod, only: inidrydep
181 : use aero_wetdep_cam, only: aero_wetdep_init
182 : use mo_setsox, only: sox_inti
183 :
184 : use modal_aero_calcsize, only: modal_aero_calcsize_init
185 : use modal_aero_coag, only: modal_aero_coag_init
186 : use aero_deposition_cam, only: aero_deposition_cam_init
187 : use modal_aero_gasaerexch, only: modal_aero_gasaerexch_init
188 : use modal_aero_newnuc, only: modal_aero_newnuc_init
189 : use modal_aero_rename, only: modal_aero_rename_init
190 :
191 : ! args
192 : type(physics_buffer_desc), pointer :: pbuf2d(:,:)
193 :
194 : ! local vars
195 : character(len=*), parameter :: subrname = 'aero_model_init'
196 : integer :: m, n, id
197 : character(len=20) :: dummy
198 :
199 : logical :: history_aerosol ! Output MAM or SECT aerosol tendencies
200 : logical :: history_chemistry, history_cesm_forcing, history_dust
201 :
202 : integer :: l
203 : character(len=6) :: test_name
204 : character(len=64) :: errmes
205 :
206 : character(len=2) :: unit_basename ! Units 'kg' or '1'
207 : integer :: errcode
208 : character(len=fieldname_len) :: field_name
209 :
210 : character(len=32) :: spec_name
211 : character(len=32) :: spec_type
212 : character(len=32) :: mode_type
213 : integer :: nspec
214 :
215 : ! aqueous chem initialization
216 2304 : call sox_inti()
217 :
218 2304 : dgnum_idx = pbuf_get_index('DGNUM')
219 2304 : dgnumwet_idx = pbuf_get_index('DGNUMWET')
220 2304 : fracis_idx = pbuf_get_index('FRACIS')
221 2304 : prain_idx = pbuf_get_index('PRAIN')
222 2304 : rprddp_idx = pbuf_get_index('RPRDDP')
223 2304 : rprdsh_idx = pbuf_get_index('RPRDSH')
224 2304 : nevapr_shcu_idx = pbuf_get_index('NEVAPR_SHCU')
225 2304 : nevapr_dpcu_idx = pbuf_get_index('NEVAPR_DPCU')
226 2304 : sulfeq_idx = pbuf_get_index('MAMH2SO4EQ',errcode)
227 :
228 : call phys_getopts(history_aerosol_out = history_aerosol, &
229 : history_chemistry_out=history_chemistry, &
230 : history_cesm_forcing_out=history_cesm_forcing, &
231 2304 : history_dust_out=history_dust)
232 :
233 2304 : call rad_cnst_get_info(0, nmodes=nmodes)
234 :
235 2304 : call modal_aero_data_init(pbuf2d)
236 2304 : call modal_aero_bcscavcoef_init()
237 :
238 2304 : call modal_aero_rename_init( modal_accum_coarse_exch )
239 : ! calcsize call must follow rename call
240 2304 : call modal_aero_calcsize_init( pbuf2d )
241 2304 : call modal_aero_gasaerexch_init
242 : ! coag call must follow gasaerexch call
243 2304 : call modal_aero_coag_init
244 2304 : call modal_aero_newnuc_init
245 :
246 : ! call aero_deposition_cam_init only if the user has not specified
247 : ! prescribed aerosol deposition fluxes
248 2304 : if (.not.aerodep_flx_prescribed()) then
249 2304 : aero_props => modal_aerosol_properties()
250 2304 : call aero_deposition_cam_init(aero_props)
251 : endif
252 :
253 2304 : call dust_init()
254 2304 : call seasalt_init(seasalt_emis_scale)
255 :
256 2304 : ndrydep = 0
257 :
258 96768 : count_species: do m = 1,pcnst
259 96768 : if ( len_trim(drydep_list(m)) /= 0 ) then
260 43776 : ndrydep = ndrydep+1
261 : endif
262 : enddo count_species
263 :
264 2304 : if (ndrydep>0) &
265 6912 : allocate(drydep_indices(ndrydep))
266 :
267 46080 : do m = 1,ndrydep
268 43776 : call cnst_get_ind ( drydep_list(m), id, abort=.false. )
269 43776 : if (id>0) then
270 43776 : drydep_indices(m) = id
271 : else
272 0 : call endrun(subrname//': invalid drydep species: '//trim(drydep_list(m)) )
273 : endif
274 :
275 46080 : if (masterproc) then
276 57 : write(iulog,*) subrname//': '//drydep_list(m)//' will have drydep applied'
277 : endif
278 : enddo
279 :
280 2304 : if (ndrydep>0) then
281 :
282 2304 : call inidrydep(rair, gravit)
283 :
284 2304 : dummy = 'RAM1'
285 2304 : call addfld (dummy,horiz_only, 'A','frac','RAM1')
286 2304 : if ( history_aerosol ) then
287 0 : call add_default (dummy, 1, ' ')
288 : endif
289 2304 : dummy = 'airFV'
290 2304 : call addfld (dummy,horiz_only, 'A','frac','FV')
291 2304 : if ( history_aerosol ) then
292 0 : call add_default (dummy, 1, ' ')
293 : endif
294 :
295 : endif
296 :
297 2304 : if (dust_active) then
298 : ! emissions diagnostics ....
299 :
300 16128 : do m = 1, dust_nbin+dust_nnum
301 13824 : dummy = trim(dust_names(m)) // 'SF'
302 13824 : call addfld (dummy,horiz_only, 'A','kg/m2/s',trim(dust_names(m))//' dust surface emission')
303 16128 : if (history_aerosol.or.history_chemistry) then
304 13824 : call add_default (dummy, 1, ' ')
305 : endif
306 : enddo
307 :
308 2304 : dummy = 'DSTSFMBL'
309 2304 : call addfld (dummy,horiz_only, 'A','kg/m2/s','Mobilization flux at surface')
310 2304 : if (history_aerosol .or. history_dust) then
311 0 : call add_default (dummy, 1, ' ')
312 : endif
313 :
314 2304 : dummy = 'LND_MBL'
315 2304 : call addfld (dummy,horiz_only, 'A','frac','Soil erodibility factor')
316 2304 : if (history_aerosol) then
317 0 : call add_default (dummy, 1, ' ')
318 : endif
319 :
320 : endif
321 :
322 2304 : if (seasalt_active) then
323 :
324 2304 : dummy = 'SSTSFMBL'
325 2304 : call addfld (dummy,horiz_only, 'A','kg/m2/s','Mobilization flux at surface')
326 2304 : if (history_aerosol) then
327 0 : call add_default (dummy, 1, ' ')
328 : endif
329 :
330 9216 : do m = 1, seasalt_nbin
331 6912 : dummy = trim(seasalt_names(m)) // 'SF'
332 6912 : call addfld (dummy,horiz_only, 'A','kg/m2/s',trim(seasalt_names(m))//' seasalt surface emission')
333 9216 : if (history_aerosol.or.history_chemistry) then
334 6912 : call add_default (dummy, 1, ' ')
335 : endif
336 : enddo
337 :
338 : endif
339 :
340 :
341 : ! set flags for drydep tendencies
342 2304 : drydep_lq(:) = .false.
343 46080 : do m=1,ndrydep
344 43776 : id = drydep_indices(m)
345 46080 : drydep_lq(id) = .true.
346 : enddo
347 :
348 2304 : wetdens_ap_idx = pbuf_get_index('WETDENS_AP')
349 2304 : qaerwat_idx = pbuf_get_index('QAERWAT')
350 2304 : pblh_idx = pbuf_get_index('pblh')
351 :
352 2304 : rate1_cw2pr_st_idx = pbuf_get_index('RATE1_CW2PR_ST')
353 2304 : call pbuf_set_field(pbuf2d, rate1_cw2pr_st_idx, 0.0_r8)
354 :
355 46080 : do m = 1,ndrydep
356 :
357 : ! units
358 43776 : if (drydep_list(m)(1:3) == 'num') then
359 9216 : unit_basename = ' 1'
360 : else
361 34560 : unit_basename = 'kg'
362 : endif
363 :
364 : call addfld (trim(drydep_list(m))//'DDF', horiz_only, 'A',unit_basename//'/m2/s ', &
365 43776 : trim(drydep_list(m))//' dry deposition flux at bottom (grav + turb)')
366 : call addfld (trim(drydep_list(m))//'TBF', horiz_only, 'A',unit_basename//'/m2/s', &
367 43776 : trim(drydep_list(m))//' turbulent dry deposition flux')
368 : call addfld (trim(drydep_list(m))//'GVF', horiz_only, 'A',unit_basename//'/m2/s ', &
369 43776 : trim(drydep_list(m))//' gravitational dry deposition flux')
370 : call addfld (trim(drydep_list(m))//'DTQ', (/ 'lev' /), 'A',unit_basename//'/kg/s ', &
371 87552 : trim(drydep_list(m))//' dry deposition')
372 : call addfld (trim(drydep_list(m))//'DDV', (/ 'lev' /), 'A','m/s', &
373 87552 : trim(drydep_list(m))//' deposition velocity')
374 :
375 43776 : if ( history_aerosol.or.history_chemistry ) then
376 43776 : call add_default (trim(drydep_list(m))//'DDF', 1, ' ')
377 : endif
378 46080 : if ( history_aerosol ) then
379 0 : call add_default (trim(drydep_list(m))//'TBF', 1, ' ')
380 0 : call add_default (trim(drydep_list(m))//'GVF', 1, ' ')
381 : endif
382 :
383 : enddo
384 :
385 73728 : do m = 1,gas_pcnst
386 :
387 71424 : if ( solsym(m)(1:3) == 'num') then
388 9216 : unit_basename = ' 1' ! Units 'kg' or '1'
389 : else
390 62208 : unit_basename = 'kg' ! Units 'kg' or '1'
391 : end if
392 :
393 : call addfld( 'GS_'//trim(solsym(m)),horiz_only, 'A', unit_basename//'/m2/s ', &
394 71424 : trim(solsym(m))//' gas chemistry/wet removal (for gas species)')
395 : call addfld( 'AQ_'//trim(solsym(m)),horiz_only, 'A', unit_basename//'/m2/s ', &
396 71424 : trim(solsym(m))//' aqueous chemistry (for gas species)')
397 73728 : if ( history_aerosol ) then
398 0 : call add_default( 'AQ_'//trim(solsym(m)), 1, ' ')
399 : endif
400 :
401 : enddo
402 96768 : do n = 1,pcnst
403 96768 : if( .not. (cnst_name_cw(n) == ' ') ) then
404 :
405 43776 : if (cnst_name_cw(n)(1:3) == 'num') then
406 9216 : unit_basename = ' 1'
407 : else
408 34560 : unit_basename = 'kg'
409 : endif
410 :
411 : call addfld( cnst_name_cw(n), (/ 'lev' /), 'A', unit_basename//'/kg ', &
412 87552 : trim(cnst_name_cw(n))//' in cloud water')
413 43776 : call addfld (trim(cnst_name_cw(n))//'DDF', horiz_only, 'A', unit_basename//'/m2/s ', &
414 87552 : trim(cnst_name_cw(n))//' dry deposition flux at bottom (grav + turb)')
415 43776 : call addfld (trim(cnst_name_cw(n))//'TBF', horiz_only, 'A', unit_basename//'/m2/s ', &
416 87552 : trim(cnst_name_cw(n))//' turbulent dry deposition flux')
417 43776 : call addfld (trim(cnst_name_cw(n))//'GVF', horiz_only, 'A', unit_basename//'/m2/s ', &
418 87552 : trim(cnst_name_cw(n))//' gravitational dry deposition flux')
419 :
420 43776 : if ( history_aerosol.or. history_chemistry ) then
421 43776 : call add_default( cnst_name_cw(n), 1, ' ' )
422 : endif
423 43776 : if ( history_aerosol ) then
424 0 : call add_default (trim(cnst_name_cw(n))//'GVF', 1, ' ')
425 0 : call add_default (trim(cnst_name_cw(n))//'TBF', 1, ' ')
426 0 : call add_default (trim(cnst_name_cw(n))//'DDF', 1, ' ')
427 : endif
428 : endif
429 : enddo
430 :
431 9216 : allocate(dgnum_name(ntot_amode), dgnumwet_name(ntot_amode))
432 11520 : do n=1,ntot_amode
433 9216 : dgnum_name(n) = ' '
434 9216 : dgnumwet_name(n) = ' '
435 9216 : write(dgnum_name(n),fmt='(a,i1)') 'dgnum',n
436 9216 : write(dgnumwet_name(n),fmt='(a,i1)') 'dgnumwet',n
437 18432 : call addfld( dgnum_name(n), (/ 'lev' /), 'I', 'm', 'Aerosol mode dry diameter' )
438 18432 : call addfld( dgnumwet_name(n), (/ 'lev' /), 'I', 'm', 'Aerosol mode wet diameter' )
439 9216 : if ( history_aerosol ) then
440 0 : call add_default( dgnum_name(n), 1, ' ' )
441 0 : call add_default( dgnumwet_name(n), 1, ' ' )
442 : endif
443 9216 : if ( history_cesm_forcing .and. n<4 ) then
444 0 : call add_default( dgnumwet_name(n), 8, ' ' )
445 : endif
446 :
447 11520 : if (modal_strat_sulfate) then
448 0 : field_name = ' '
449 0 : write(field_name,fmt='(a,i1)') 'wtpct_a',n
450 0 : call addfld( field_name, (/ 'lev' /), 'I', '%', 'Aerosol mode weight percent H2SO4' )
451 0 : if ( history_aerosol ) then
452 0 : call add_default (field_name, 0, 'I')
453 : endif
454 :
455 0 : field_name = ' '
456 0 : write(field_name,fmt='(a,i1)') 'sulfeq_a',n
457 0 : call addfld( field_name, (/ 'lev' /), 'I', 'kg/kg', 'H2SO4 equilibrium mixing ratio' )
458 0 : if ( history_aerosol ) then
459 0 : call add_default (field_name, 0, 'I')
460 : endif
461 :
462 0 : field_name = ' '
463 0 : write(field_name,fmt='(a,i1)') 'sulden_a',n
464 0 : call addfld( field_name, (/ 'lev' /), 'I', 'g/cm3', 'Sulfate aerosol particle mass density' )
465 0 : if ( history_aerosol ) then
466 0 : call add_default (field_name, 0, 'I')
467 : endif
468 :
469 : end if
470 : end do
471 :
472 2304 : ndx_h2so4 = get_spc_ndx('H2SO4')
473 2304 : nh3_ndx = get_spc_ndx('NH3')
474 2304 : nh4_ndx = get_spc_ndx('NH4')
475 :
476 6912 : allocate(num_idx(ntot_amode))
477 11520 : num_idx = -1
478 :
479 : ! for aero_model_surfarea called from mo_usrrxt
480 11520 : do l=1,ntot_amode
481 9216 : test_name = ' '
482 9216 : write(test_name,fmt='(a5,i1)') 'num_a',l
483 9216 : num_idx(l) = get_spc_ndx( trim(test_name) )
484 11520 : if (num_idx(l) < 0) then
485 0 : write(errmes,fmt='(a,i1)') 'usrrxt_inti: cannot find MAM num_idx ',l
486 0 : write(iulog,*) errmes
487 0 : call endrun(errmes)
488 : endif
489 : end do
490 :
491 9216 : allocate(index_tot_mass(nmodes,nspec_max))
492 6912 : allocate(index_chm_mass(nmodes,nspec_max))
493 71424 : index_tot_mass = -1
494 71424 : index_chm_mass = -1
495 :
496 : ! for surf_area_dens
497 : ! define indices associated with the various aerosol types
498 11520 : do n = 1,nmodes
499 9216 : call rad_cnst_get_info(0, n, mode_type=mode_type, nspec=nspec)
500 11520 : if ( trim(mode_type) /= 'primary_carbon') then ! ignore the primary_carbon mode
501 36864 : do l = 1, nspec
502 29952 : call rad_cnst_get_info(0, n, l, spec_type=spec_type, spec_name=spec_name)
503 29952 : index_tot_mass(n,l) = get_spc_ndx(spec_name)
504 : if ( trim(spec_type) == 'sulfate' .or. &
505 : trim(spec_type) == 's-organic' .or. &
506 : trim(spec_type) == 'p-organic' .or. &
507 29952 : trim(spec_type) == 'black-c' .or. &
508 6912 : trim(spec_type) == 'ammonium') then
509 16128 : index_chm_mass(n,l) = get_spc_ndx(spec_name)
510 : endif
511 : enddo
512 : endif
513 : enddo
514 :
515 2304 : if (has_sox) then
516 11520 : do m = 1, ntot_amode
517 :
518 9216 : l = lptr_so4_cw_amode(m)
519 11520 : if (l > 0) then
520 : call addfld (&
521 6912 : trim(cnst_name_cw(l))//'AQSO4',horiz_only, 'A','kg/m2/s', &
522 13824 : trim(cnst_name_cw(l))//' aqueous phase chemistry')
523 : call addfld (&
524 6912 : trim(cnst_name_cw(l))//'AQH2SO4',horiz_only, 'A','kg/m2/s', &
525 13824 : trim(cnst_name_cw(l))//' aqueous phase chemistry')
526 6912 : if ( history_aerosol ) then
527 0 : call add_default (trim(cnst_name_cw(l))//'AQSO4', 1, ' ')
528 0 : call add_default (trim(cnst_name_cw(l))//'AQH2SO4', 1, ' ')
529 : endif
530 : end if
531 :
532 : end do
533 :
534 4608 : call addfld( 'XPH_LWC', (/ 'lev' /), 'A','kg/kg', 'pH value multiplied by lwc')
535 2304 : call addfld ('AQSO4_H2O2', horiz_only, 'A','kg/m2/s', 'SO4 aqueous phase chemistry due to H2O2')
536 2304 : call addfld ('AQSO4_O3', horiz_only, 'A','kg/m2/s', 'SO4 aqueous phase chemistry due to O3')
537 :
538 2304 : if ( history_aerosol ) then
539 0 : call add_default ('XPH_LWC', 1, ' ')
540 0 : call add_default ('AQSO4_H2O2', 1, ' ')
541 0 : call add_default ('AQSO4_O3', 1, ' ')
542 : endif
543 : endif
544 :
545 2304 : call aero_wetdep_init()
546 :
547 4608 : end subroutine aero_model_init
548 :
549 : !=============================================================================
550 : !=============================================================================
551 3860712 : subroutine aero_model_drydep ( state, pbuf, obklen, ustar, cam_in, dt, cam_out, ptend )
552 :
553 2304 : use dust_sediment_mod, only: dust_sediment_tend
554 : use aer_drydep_mod, only: d3ddflux, calcram
555 : use modal_aero_data, only: qqcw_get_field
556 : use modal_aero_data, only: cnst_name_cw
557 : use modal_aero_data, only: alnsg_amode
558 : use modal_aero_data, only: sigmag_amode
559 : use modal_aero_data, only: nspec_amode
560 : use modal_aero_data, only: numptr_amode
561 : use modal_aero_data, only: numptrcw_amode
562 : use modal_aero_data, only: lmassptr_amode
563 : use modal_aero_data, only: lmassptrcw_amode
564 : use aero_deposition_cam,only: aero_deposition_cam_setdry
565 :
566 : ! args
567 : type(physics_state), intent(in) :: state ! Physics state variables
568 : real(r8), intent(in) :: obklen(:)
569 : real(r8), intent(in) :: ustar(:) ! sfc fric vel
570 : type(cam_in_t), target, intent(in) :: cam_in ! import state
571 : real(r8), intent(in) :: dt ! time step
572 : type(cam_out_t), intent(inout) :: cam_out ! export state
573 : type(physics_ptend), intent(out) :: ptend ! indivdual parameterization tendencies
574 : type(physics_buffer_desc), pointer :: pbuf(:)
575 :
576 : ! local vars
577 89784 : real(r8), pointer :: landfrac(:) ! land fraction
578 89784 : real(r8), pointer :: icefrac(:) ! ice fraction
579 89784 : real(r8), pointer :: ocnfrac(:) ! ocean fraction
580 89784 : real(r8), pointer :: fvin(:) !
581 89784 : real(r8), pointer :: ram1in(:) ! for dry dep velocities from land model for progseasalts
582 :
583 : real(r8) :: fv(pcols) ! for dry dep velocities, from land modified over ocean & ice
584 : real(r8) :: ram1(pcols) ! for dry dep velocities, from land modified over ocean & ice
585 :
586 : integer :: lchnk ! chunk identifier
587 : integer :: ncol ! number of atmospheric columns
588 : integer :: jvlc ! index for last dimension of vlc_xxx arrays
589 : integer :: lphase ! index for interstitial / cloudborne aerosol
590 : integer :: lspec ! index for aerosol number / chem-mass / water-mass
591 : integer :: m ! aerosol mode index
592 : integer :: mm ! tracer index
593 : integer :: i
594 :
595 : real(r8) :: tvs(pcols,pver)
596 : real(r8) :: rho(pcols,pver) ! air density in kg/m3
597 : real(r8) :: sflx(pcols) ! deposition flux
598 : real(r8) :: dep_trb(pcols) !kg/m2/s
599 : real(r8) :: dep_grv(pcols) !kg/m2/s (total of grav and trb)
600 : real(r8) :: pvmzaer(pcols,pverp) ! sedimentation velocity in Pa
601 : real(r8) :: dqdt_tmp(pcols,pver) ! temporary array to hold tendency for 1 species
602 :
603 : real(r8) :: rad_drop(pcols,pver)
604 : real(r8) :: dens_drop(pcols,pver)
605 : real(r8) :: sg_drop(pcols,pver)
606 : real(r8) :: rad_aer(pcols,pver)
607 : real(r8) :: dens_aer(pcols,pver)
608 : real(r8) :: sg_aer(pcols,pver)
609 :
610 : real(r8) :: vlc_dry(pcols,pver,4) ! dep velocity
611 : real(r8) :: vlc_grv(pcols,pver,4) ! dep velocity
612 : real(r8):: vlc_trb(pcols,4) ! dep velocity
613 : real(r8) :: aerdepdryis(pcols,pcnst) ! aerosol dry deposition (interstitial)
614 : real(r8) :: aerdepdrycw(pcols,pcnst) ! aerosol dry deposition (cloud water)
615 89784 : real(r8), pointer :: fldcw(:,:)
616 89784 : real(r8), pointer :: dgncur_awet(:,:,:)
617 89784 : real(r8), pointer :: wetdens(:,:,:)
618 89784 : real(r8), pointer :: qaerwat(:,:,:)
619 :
620 89784 : landfrac => cam_in%landfrac(:)
621 89784 : icefrac => cam_in%icefrac(:)
622 89784 : ocnfrac => cam_in%ocnfrac(:)
623 89784 : fvin => cam_in%fv(:)
624 89784 : ram1in => cam_in%ram1(:)
625 :
626 89784 : lchnk = state%lchnk
627 89784 : ncol = state%ncol
628 :
629 : ! calc ram and fv over ocean and sea ice ...
630 : call calcram( ncol,landfrac,icefrac,ocnfrac,obklen,&
631 : ustar,ram1in,ram1,state%t(:,pver),state%pmid(:,pver),&
632 89784 : state%pdel(:,pver),fvin,fv)
633 :
634 89784 : call outfld( 'airFV', fv(:), pcols, lchnk )
635 89784 : call outfld( 'RAM1', ram1(:), pcols, lchnk )
636 :
637 : ! note that tendencies are not only in sfc layer (because of sedimentation)
638 : ! and that ptend is updated within each subroutine for different species
639 :
640 89784 : call physics_ptend_init(ptend, state%psetcols, 'aero_model_drydep', lq=drydep_lq)
641 :
642 359136 : call pbuf_get_field(pbuf, dgnumwet_idx, dgncur_awet, start=(/1,1,1/), kount=(/pcols,pver,nmodes/) )
643 359136 : call pbuf_get_field(pbuf, wetdens_ap_idx, wetdens, start=(/1,1,1/), kount=(/pcols,pver,nmodes/) )
644 359136 : call pbuf_get_field(pbuf, qaerwat_idx, qaerwat, start=(/1,1,1/), kount=(/pcols,pver,nmodes/) )
645 :
646 139513896 : tvs(:ncol,:) = state%t(:ncol,:)!*(1+state%q(:ncol,k)
647 139513896 : rho(:ncol,:)= state%pmid(:ncol,:)/(rair*state%t(:ncol,:))
648 :
649 : !
650 : ! calc settling/deposition velocities for cloud droplets (and cloud-borne aerosols)
651 : !
652 : ! *** mean drop radius should eventually be computed from ndrop and qcldwtr
653 142038288 : rad_drop(:,:) = 5.0e-6_r8
654 142038288 : dens_drop(:,:) = rhoh2o
655 142038288 : sg_drop(:,:) = 1.46_r8
656 89784 : jvlc = 3
657 : call modal_aero_depvel_part( ncol,state%t(:,:), state%pmid(:,:), ram1, fv, &
658 : vlc_dry(:,:,jvlc), vlc_trb(:,jvlc), vlc_grv(:,:,jvlc), &
659 89784 : rad_drop(:,:), dens_drop(:,:), sg_drop(:,:), 0, lchnk)
660 89784 : jvlc = 4
661 : call modal_aero_depvel_part( ncol,state%t(:,:), state%pmid(:,:), ram1, fv, &
662 : vlc_dry(:,:,jvlc), vlc_trb(:,jvlc), vlc_grv(:,:,jvlc), &
663 89784 : rad_drop(:,:), dens_drop(:,:), sg_drop(:,:), 3, lchnk)
664 :
665 :
666 :
667 448920 : do m = 1, ntot_amode ! main loop over aerosol modes
668 :
669 1167192 : do lphase = 1, 2 ! loop over interstitial / cloud-borne forms
670 :
671 718272 : if (lphase == 1) then ! interstial aerosol - calc settling/dep velocities of mode
672 :
673 : ! rad_aer = volume mean wet radius (m)
674 : ! dgncur_awet = geometric mean wet diameter for number distribution (m)
675 718272 : rad_aer(1:ncol,:) = 0.5_r8*dgncur_awet(1:ncol,:,m) &
676 558773856 : *exp(1.5_r8*(alnsg_amode(m)**2))
677 : ! dens_aer(1:ncol,:) = wet density (kg/m3)
678 558055584 : dens_aer(1:ncol,:) = wetdens(1:ncol,:,m)
679 558055584 : sg_aer(1:ncol,:) = sigmag_amode(m)
680 :
681 359136 : jvlc = 1
682 : call modal_aero_depvel_part( ncol, state%t(:,:), state%pmid(:,:), ram1, fv, &
683 : vlc_dry(:,:,jvlc), vlc_trb(:,jvlc), vlc_grv(:,:,jvlc), &
684 359136 : rad_aer(:,:), dens_aer(:,:), sg_aer(:,:), 0, lchnk)
685 359136 : jvlc = 2
686 : call modal_aero_depvel_part( ncol, state%t(:,:), state%pmid(:,:), ram1, fv, &
687 : vlc_dry(:,:,jvlc), vlc_trb(:,jvlc), vlc_grv(:,:,jvlc), &
688 359136 : rad_aer(:,:), dens_aer(:,:), sg_aer(:,:), 3, lchnk)
689 : end if
690 :
691 5207472 : do lspec = 0, nspec_amode(m)+1 ! loop over number + constituents + water
692 :
693 4130064 : if (lspec == 0) then ! number
694 718272 : if (lphase == 1) then
695 359136 : mm = numptr_amode(m)
696 359136 : jvlc = 1
697 : else
698 359136 : mm = numptrcw_amode(m)
699 359136 : jvlc = 3
700 : endif
701 3411792 : else if (lspec <= nspec_amode(m)) then ! non-water mass
702 2693520 : if (lphase == 1) then
703 1346760 : mm = lmassptr_amode(lspec,m)
704 1346760 : jvlc = 2
705 : else
706 1346760 : mm = lmassptrcw_amode(lspec,m)
707 1346760 : jvlc = 4
708 : endif
709 : else ! water mass
710 : ! bypass dry deposition of aerosol water
711 : cycle
712 : if (lphase == 1) then
713 : mm = 0
714 : ! mm = lwaterptr_amode(m)
715 : jvlc = 2
716 : else
717 : mm = 0
718 : jvlc = 4
719 : endif
720 : endif
721 :
722 :
723 3411792 : if (mm <= 0) cycle
724 :
725 : ! if (lphase == 1) then
726 4130064 : if ((lphase == 1) .and. (lspec <= nspec_amode(m))) then
727 1705896 : ptend%lq(mm) = .TRUE.
728 :
729 : ! use pvprogseasalts instead (means making the top level 0)
730 28484496 : pvmzaer(:ncol,1)=0._r8
731 2650764024 : pvmzaer(:ncol,2:pverp) = vlc_dry(:ncol,:,jvlc)
732 :
733 1705896 : call outfld( trim(cnst_name(mm))//'DDV', pvmzaer(:,2:pverp), pcols, lchnk )
734 :
735 : if(.true.) then ! use phil's method
736 : ! convert from meters/sec to pascals/sec
737 : ! pvprogseasalts(:,1) is assumed zero, use density from layer above in conversion
738 2650764024 : pvmzaer(:ncol,2:pverp) = pvmzaer(:ncol,2:pverp) * rho(:ncol,:)*gravit
739 :
740 : ! calculate the tendencies and sfc fluxes from the above velocities
741 : call dust_sediment_tend( &
742 : ncol, dt, state%pint(:,:), state%pmid, state%pdel, state%t , &
743 1705896 : state%q(:,:,mm), pvmzaer, ptend%q(:,:,mm), sflx )
744 : else !use charlie's method
745 : call d3ddflux( ncol, vlc_dry(:,:,jvlc), state%q(:,:,mm), state%pmid, &
746 : state%pdel, tvs, sflx, ptend%q(:,:,mm), dt )
747 : endif
748 :
749 : ! apportion dry deposition into turb and gravitational settling for tapes
750 1705896 : dep_trb = 0._r8
751 1705896 : dep_grv = 0._r8
752 28484496 : do i=1,ncol
753 28484496 : if (vlc_dry(i,pver,jvlc) /= 0._r8) then
754 26778600 : dep_trb(i)=sflx(i)*vlc_trb(i,jvlc)/vlc_dry(i,pver,jvlc)
755 26778600 : dep_grv(i)=sflx(i)*vlc_grv(i,pver,jvlc)/vlc_dry(i,pver,jvlc)
756 : end if
757 : enddo
758 :
759 1705896 : call outfld( trim(cnst_name(mm))//'DDF', sflx, pcols, lchnk)
760 1705896 : call outfld( trim(cnst_name(mm))//'TBF', dep_trb, pcols, lchnk )
761 1705896 : call outfld( trim(cnst_name(mm))//'GVF', dep_grv, pcols, lchnk )
762 1705896 : call outfld( trim(cnst_name(mm))//'DTQ', ptend%q(:,:,mm), pcols, lchnk)
763 28484496 : aerdepdryis(:ncol,mm) = sflx(:ncol)
764 :
765 1705896 : else if ((lphase == 1) .and. (lspec == nspec_amode(m)+1)) then ! aerosol water
766 : ! use pvprogseasalts instead (means making the top level 0)
767 0 : pvmzaer(:ncol,1)=0._r8
768 0 : pvmzaer(:ncol,2:pverp) = vlc_dry(:ncol,:,jvlc)
769 :
770 : if(.true.) then ! use phil's method
771 : ! convert from meters/sec to pascals/sec
772 : ! pvprogseasalts(:,1) is assumed zero, use density from layer above in conversion
773 0 : pvmzaer(:ncol,2:pverp) = pvmzaer(:ncol,2:pverp) * rho(:ncol,:)*gravit
774 :
775 : ! calculate the tendencies and sfc fluxes from the above velocities
776 : call dust_sediment_tend( &
777 : ncol, dt, state%pint(:,:), state%pmid, state%pdel, state%t , &
778 0 : qaerwat(:,:,mm), pvmzaer, dqdt_tmp(:,:), sflx )
779 : else !use charlie's method
780 : call d3ddflux( ncol, vlc_dry(:,:,jvlc), qaerwat(:,:,mm), state%pmid, &
781 : state%pdel, tvs, sflx, dqdt_tmp(:,:), dt )
782 : endif
783 :
784 : ! apportion dry deposition into turb and gravitational settling for tapes
785 0 : dep_trb = 0._r8
786 0 : dep_grv = 0._r8
787 0 : do i=1,ncol
788 0 : if (vlc_dry(i,pver,jvlc) /= 0._r8) then
789 0 : dep_trb(i)=sflx(i)*vlc_trb(i,jvlc)/vlc_dry(i,pver,jvlc)
790 0 : dep_grv(i)=sflx(i)*vlc_grv(i,pver,jvlc)/vlc_dry(i,pver,jvlc)
791 : end if
792 : enddo
793 :
794 0 : qaerwat(1:ncol,:,mm) = qaerwat(1:ncol,:,mm) + dqdt_tmp(1:ncol,:) * dt
795 :
796 : else ! lphase == 2
797 : ! use pvprogseasalts instead (means making the top level 0)
798 28484496 : pvmzaer(:ncol,1)=0._r8
799 2650764024 : pvmzaer(:ncol,2:pverp) = vlc_dry(:ncol,:,jvlc)
800 1705896 : fldcw => qqcw_get_field(pbuf, mm,lchnk)
801 :
802 : if(.true.) then ! use phil's method
803 : ! convert from meters/sec to pascals/sec
804 : ! pvprogseasalts(:,1) is assumed zero, use density from layer above in conversion
805 2650764024 : pvmzaer(:ncol,2:pverp) = pvmzaer(:ncol,2:pverp) * rho(:ncol,:)*gravit
806 :
807 : ! calculate the tendencies and sfc fluxes from the above velocities
808 : call dust_sediment_tend( &
809 : ncol, dt, state%pint(:,:), state%pmid, state%pdel, state%t , &
810 1705896 : fldcw(:,:), pvmzaer, dqdt_tmp(:,:), sflx )
811 : else !use charlie's method
812 : call d3ddflux( ncol, vlc_dry(:,:,jvlc), fldcw(:,:), state%pmid, &
813 : state%pdel, tvs, sflx, dqdt_tmp(:,:), dt )
814 : endif
815 :
816 : ! apportion dry deposition into turb and gravitational settling for tapes
817 1705896 : dep_trb = 0._r8
818 1705896 : dep_grv = 0._r8
819 28484496 : do i=1,ncol
820 28484496 : if (vlc_dry(i,pver,jvlc) /= 0._r8) then
821 26778600 : dep_trb(i)=sflx(i)*vlc_trb(i,jvlc)/vlc_dry(i,pver,jvlc)
822 26778600 : dep_grv(i)=sflx(i)*vlc_grv(i,pver,jvlc)/vlc_dry(i,pver,jvlc)
823 : end if
824 : enddo
825 :
826 2650764024 : fldcw(1:ncol,:) = fldcw(1:ncol,:) + dqdt_tmp(1:ncol,:) * dt
827 :
828 1705896 : call outfld( trim(cnst_name_cw(mm))//'DDF', sflx, pcols, lchnk)
829 1705896 : call outfld( trim(cnst_name_cw(mm))//'TBF', dep_trb, pcols, lchnk )
830 1705896 : call outfld( trim(cnst_name_cw(mm))//'GVF', dep_grv, pcols, lchnk )
831 28484496 : aerdepdrycw(:ncol,mm) = sflx(:ncol)
832 :
833 : endif
834 :
835 : enddo ! lspec = 0, nspec_amode(m)+1
836 : enddo ! lphase = 1, 2
837 : enddo ! m = 1, ntot_amode
838 :
839 : ! if the user has specified prescribed aerosol dep fluxes then
840 : ! do not set cam_out dep fluxes according to the prognostic aerosols
841 89784 : if (.not.aerodep_flx_prescribed()) then
842 89784 : call aero_deposition_cam_setdry(aerdepdryis, aerdepdrycw, cam_out)
843 : endif
844 :
845 179568 : endsubroutine aero_model_drydep
846 :
847 : !=============================================================================
848 : !=============================================================================
849 3770928 : subroutine aero_model_wetdep( state, dt, dlf, cam_out, ptend, pbuf)
850 :
851 89784 : use aero_wetdep_cam, only: aero_wetdep_tend
852 :
853 : ! args
854 :
855 : type(physics_state), intent(in) :: state ! Physics state variables
856 : real(r8), intent(in) :: dt ! time step
857 : real(r8), intent(in) :: dlf(:,:) ! shallow+deep convective detrainment [kg/kg/s]
858 : type(cam_out_t), intent(inout) :: cam_out ! export state
859 : type(physics_ptend), intent(out) :: ptend ! indivdual parameterization tendencies
860 : type(physics_buffer_desc), pointer :: pbuf(:)
861 :
862 89784 : call aero_wetdep_tend(state, dt, dlf, cam_out, ptend, pbuf)
863 :
864 89784 : end subroutine aero_model_wetdep
865 :
866 : !-------------------------------------------------------------------------
867 : ! provides wet tropospheric aerosol surface area info for modal aerosols
868 : ! called from mo_usrrxt
869 : !-------------------------------------------------------------------------
870 0 : subroutine aero_model_surfarea( &
871 0 : state, mmr, radmean, relhum, pmid, temp, strato_sad, sulfate, rho, ltrop, &
872 0 : dlat, het1_ndx, pbuf, ncol, sfc, dm_aer, sad_trop, reff_trop )
873 :
874 : ! dummy args
875 : type(physics_state), intent(in) :: state ! Physics state variables
876 : real(r8), intent(in) :: pmid(:,:)
877 : real(r8), intent(in) :: temp(:,:)
878 : real(r8), intent(in) :: mmr(:,:,:)
879 : real(r8), intent(in) :: radmean ! mean radii in cm
880 : real(r8), intent(in) :: strato_sad(:,:)
881 : integer, intent(in) :: ncol
882 : integer, intent(in) :: ltrop(:)
883 : real(r8), intent(in) :: dlat(:) ! degrees latitude
884 : integer, intent(in) :: het1_ndx
885 : real(r8), intent(in) :: relhum(:,:)
886 : real(r8), intent(in) :: rho(:,:) ! total atm density (/cm^3)
887 : real(r8), intent(in) :: sulfate(:,:)
888 : type(physics_buffer_desc), pointer :: pbuf(:)
889 :
890 : real(r8), intent(inout) :: sfc(:,:,:)
891 : real(r8), intent(inout) :: dm_aer(:,:,:)
892 : real(r8), intent(inout) :: sad_trop(:,:)
893 : real(r8), intent(out) :: reff_trop(:,:)
894 :
895 : ! local vars
896 0 : real(r8), pointer, dimension(:,:,:) :: dgnumwet
897 0 : integer :: beglev(ncol)
898 0 : integer :: endlev(ncol)
899 : integer :: i,k
900 :
901 0 : call pbuf_get_field(pbuf, dgnumwet_idx, dgnumwet )
902 :
903 0 : beglev(:ncol)=ltrop(:ncol)+1
904 0 : endlev(:ncol)=pver
905 0 : call surf_area_dens( ncol, mmr, pmid, temp, dgnumwet, beglev, endlev, sad_trop, reff_trop, sfc=sfc )
906 :
907 0 : do i = 1,ncol
908 0 : do k = ltrop(i)+1,pver
909 0 : dm_aer(i,k,:) = dgnumwet(i,k,:) * 1.e2_r8 ! convert m to cm
910 : enddo
911 : enddo
912 :
913 89784 : end subroutine aero_model_surfarea
914 :
915 : !-------------------------------------------------------------------------
916 : ! provides WET stratospheric aerosol surface area info for modal aerosols
917 : ! if modal_strat_sulfate = TRUE -- called from mo_gas_phase_chemdr
918 : !-------------------------------------------------------------------------
919 0 : subroutine aero_model_strat_surfarea( state, ncol, mmr, pmid, temp, ltrop, pbuf, strato_sad, reff_strat )
920 :
921 : ! dummy args
922 : type(physics_state), intent(in) :: state ! Physics state variables
923 : integer, intent(in) :: ncol
924 : real(r8), intent(in) :: mmr(:,:,:)
925 : real(r8), intent(in) :: pmid(:,:)
926 : real(r8), intent(in) :: temp(:,:)
927 : integer, intent(in) :: ltrop(:) ! tropopause level indices
928 : type(physics_buffer_desc), pointer :: pbuf(:)
929 : real(r8), intent(out) :: strato_sad(:,:)
930 : real(r8), intent(out) :: reff_strat(:,:)
931 :
932 : ! local vars
933 0 : real(r8), pointer, dimension(:,:,:) :: dgnumwet
934 0 : integer :: beglev(ncol)
935 0 : integer :: endlev(ncol)
936 :
937 0 : reff_strat = 0._r8
938 0 : strato_sad = 0._r8
939 :
940 0 : if (.not.modal_strat_sulfate) return
941 :
942 0 : call pbuf_get_field(pbuf, dgnumwet_idx, dgnumwet )
943 :
944 0 : beglev(:ncol)=top_lev
945 0 : endlev(:ncol)=ltrop(:ncol)
946 0 : call surf_area_dens( ncol, mmr, pmid, temp, dgnumwet, beglev, endlev, strato_sad, reff_strat )
947 :
948 0 : end subroutine aero_model_strat_surfarea
949 :
950 : !=============================================================================
951 : !=============================================================================
952 89784 : subroutine aero_model_gasaerexch( state, loffset, ncol, lchnk, troplev, delt, reaction_rates, &
953 89784 : tfld, pmid, pdel, mbar, relhum, &
954 179568 : zm, qh2o, cwat, cldfr, cldnum, &
955 179568 : airdens, invariants, del_h2so4_gasprod, &
956 89784 : vmr0, vmr, pbuf )
957 :
958 : use time_manager, only : get_nstep
959 : use modal_aero_coag, only : modal_aero_coag_sub
960 : use modal_aero_gasaerexch, only : modal_aero_gasaerexch_sub
961 : use modal_aero_newnuc, only : modal_aero_newnuc_sub
962 : use modal_aero_data, only : cnst_name_cw, qqcw_get_field
963 :
964 : !-----------------------------------------------------------------------
965 : ! ... dummy arguments
966 : !-----------------------------------------------------------------------
967 : type(physics_state), intent(in) :: state ! Physics state variables
968 : integer, intent(in) :: loffset ! offset applied to modal aero "pointers"
969 : integer, intent(in) :: ncol ! number columns in chunk
970 : integer, intent(in) :: lchnk ! chunk index
971 : integer, intent(in) :: troplev(pcols)
972 : real(r8), intent(in) :: delt ! time step size (sec)
973 : real(r8), intent(in) :: reaction_rates(:,:,:) ! reaction rates
974 : real(r8), intent(in) :: tfld(:,:) ! temperature (K)
975 : real(r8), intent(in) :: pmid(:,:) ! pressure at model levels (Pa)
976 : real(r8), intent(in) :: pdel(:,:) ! pressure thickness of levels (Pa)
977 : real(r8), intent(in) :: mbar(:,:) ! mean wet atmospheric mass ( amu )
978 : real(r8), intent(in) :: relhum(:,:) ! relative humidity
979 : real(r8), intent(in) :: airdens(:,:) ! total atms density (molec/cm**3)
980 : real(r8), intent(in) :: invariants(:,:,:)
981 : real(r8), intent(in) :: del_h2so4_gasprod(:,:)
982 : real(r8), intent(in) :: zm(:,:)
983 : real(r8), intent(in) :: qh2o(:,:)
984 : real(r8), intent(in) :: cwat(:,:) ! cloud liquid water content (kg/kg)
985 : real(r8), intent(in) :: cldfr(:,:)
986 : real(r8), intent(in) :: cldnum(:,:) ! droplet number concentration (#/kg)
987 : real(r8), intent(in) :: vmr0(:,:,:) ! initial mixing ratios (before gas-phase chem changes)
988 : real(r8), intent(inout) :: vmr(:,:,:) ! mixing ratios ( vmr )
989 :
990 : type(physics_buffer_desc), pointer :: pbuf(:)
991 :
992 : ! local vars
993 :
994 : integer :: n, m
995 : integer :: i,k,l
996 : integer :: nstep
997 :
998 179568 : real(r8) :: del_h2so4_aeruptk(ncol,pver)
999 :
1000 89784 : real(r8), pointer :: dgnum(:,:,:), dgnumwet(:,:,:), wetdens(:,:,:)
1001 89784 : real(r8), pointer :: pblh(:) ! pbl height (m)
1002 :
1003 179568 : real(r8), dimension(ncol) :: wrk
1004 : character(len=32) :: name
1005 179568 : real(r8) :: dvmrcwdt(ncol,pver,gas_pcnst)
1006 179568 : real(r8) :: dvmrdt(ncol,pver,gas_pcnst)
1007 179568 : real(r8) :: vmrcw(ncol,pver,gas_pcnst) ! cloud-borne aerosol (vmr)
1008 :
1009 179568 : real(r8) :: aqso4(ncol,ntot_amode) ! aqueous phase chemistry
1010 179568 : real(r8) :: aqh2so4(ncol,ntot_amode) ! aqueous phase chemistry
1011 179568 : real(r8) :: aqso4_h2o2(ncol) ! SO4 aqueous phase chemistry due to H2O2
1012 179568 : real(r8) :: aqso4_o3(ncol) ! SO4 aqueous phase chemistry due to O3
1013 179568 : real(r8) :: xphlwc(ncol,pver) ! pH value multiplied by lwc
1014 179568 : real(r8) :: nh3_beg(ncol,pver)
1015 89784 : real(r8), pointer :: fldcw(:,:)
1016 89784 : real(r8), pointer :: sulfeq(:,:,:)
1017 :
1018 : !
1019 : ! ... initialize nh3
1020 : !
1021 89784 : if ( nh3_ndx > 0 ) then
1022 0 : nh3_beg = vmr(1:ncol,:,nh3_ndx)
1023 : end if
1024 : !
1025 :
1026 89784 : call pbuf_get_field(pbuf, dgnum_idx, dgnum)
1027 89784 : call pbuf_get_field(pbuf, dgnumwet_idx, dgnumwet )
1028 89784 : call pbuf_get_field(pbuf, wetdens_ap_idx, wetdens )
1029 89784 : call pbuf_get_field(pbuf, pblh_idx, pblh)
1030 :
1031 448920 : do n=1,ntot_amode
1032 359136 : call outfld(dgnum_name(n), dgnum(1:ncol,1:pver,n), ncol, lchnk )
1033 448920 : call outfld(dgnumwet_name(n), dgnumwet(1:ncol,1:pver,n), ncol, lchnk )
1034 : end do
1035 :
1036 : ! do gas-aerosol exchange (h2so4, msa, nh3 condensation)
1037 :
1038 89784 : nstep = get_nstep()
1039 :
1040 : ! calculate tendency due to gas phase chemistry and processes
1041 4325020560 : dvmrdt(:ncol,:,:) = (vmr(:ncol,:,:) - vmr0(:ncol,:,:)) / delt
1042 2873088 : do m = 1, gas_pcnst
1043 46474704 : wrk(:) = 0._r8
1044 261630576 : do k = 1,pver
1045 4324930776 : wrk(:ncol) = wrk(:ncol) + dvmrdt(:ncol,k,m)*adv_mass(m)/mbar(:ncol,k)*pdel(:ncol,k)/gravit
1046 : end do
1047 2783304 : name = 'GS_'//trim(solsym(m))
1048 2873088 : call outfld( name, wrk(:ncol), ncol, lchnk )
1049 : enddo
1050 :
1051 : !
1052 : ! Aerosol processes ...
1053 : !
1054 89784 : call qqcw2vmr( lchnk, vmrcw, mbar, ncol, loffset, pbuf )
1055 :
1056 4325020560 : dvmrdt(:ncol,:,:) = vmr(:ncol,:,:)
1057 4325020560 : dvmrcwdt(:ncol,:,:) = vmrcw(:ncol,:,:)
1058 :
1059 : ! aqueous chemistry ...
1060 :
1061 89784 : if( has_sox ) then
1062 : call setsox( state, &
1063 : ncol, &
1064 : lchnk, &
1065 : loffset, &
1066 : delt, &
1067 : pmid, &
1068 : pdel, &
1069 : tfld, &
1070 : mbar, &
1071 : cwat, &
1072 : cldfr, &
1073 : cldnum, &
1074 : airdens, &
1075 : invariants, &
1076 : vmrcw, &
1077 : vmr, &
1078 : xphlwc, &
1079 : aqso4, &
1080 : aqh2so4, &
1081 : aqso4_h2o2, &
1082 : aqso4_o3 &
1083 89784 : )
1084 :
1085 448920 : do n = 1, ntot_amode
1086 359136 : l = lptr_so4_cw_amode(n)
1087 448920 : if (l > 0) then
1088 269352 : call outfld( trim(cnst_name_cw(l))//'AQSO4', aqso4(:ncol,n), ncol, lchnk)
1089 269352 : call outfld( trim(cnst_name_cw(l))//'AQH2SO4', aqh2so4(:ncol,n), ncol, lchnk)
1090 : end if
1091 : end do
1092 :
1093 89784 : call outfld( 'AQSO4_H2O2', aqso4_h2o2(:ncol), ncol, lchnk)
1094 89784 : call outfld( 'AQSO4_O3', aqso4_o3(:ncol), ncol, lchnk)
1095 89784 : call outfld( 'XPH_LWC', xphlwc(:ncol,:), ncol, lchnk )
1096 :
1097 : endif
1098 :
1099 : ! Tendency due to aqueous chemistry
1100 4325110344 : dvmrdt = (vmr - dvmrdt) / delt
1101 4325020560 : dvmrcwdt = (vmrcw - dvmrcwdt) / delt
1102 2873088 : do m = 1, gas_pcnst
1103 46474704 : wrk(:) = 0._r8
1104 261630576 : do k = 1,pver
1105 4324930776 : wrk(:ncol) = wrk(:ncol) + dvmrdt(:ncol,k,m) * adv_mass(m)/mbar(:ncol,k)*pdel(:ncol,k)/gravit
1106 : end do
1107 2783304 : name = 'AQ_'//trim(solsym(m))
1108 2873088 : call outfld( name, wrk(:ncol), ncol, lchnk )
1109 : enddo
1110 :
1111 : ! do gas-aerosol exchange (h2so4, msa, nh3 condensation)
1112 :
1113 89784 : if (ndx_h2so4 > 0) then
1114 139513896 : del_h2so4_aeruptk(1:ncol,:) = vmr(1:ncol,:,ndx_h2so4)
1115 : else
1116 0 : del_h2so4_aeruptk(:,:) = 0.0_r8
1117 : endif
1118 :
1119 89784 : call t_startf('modal_gas-aer_exchng')
1120 :
1121 89784 : if ( sulfeq_idx>0 ) then
1122 0 : call pbuf_get_field( pbuf, sulfeq_idx, sulfeq )
1123 : else
1124 89784 : nullify( sulfeq )
1125 : endif
1126 :
1127 : call modal_aero_gasaerexch_sub( &
1128 : lchnk, ncol, nstep, &
1129 : loffset, delt, &
1130 : tfld, pmid, pdel, &
1131 : qh2o, troplev, &
1132 : vmr, vmrcw, &
1133 : dvmrdt, dvmrcwdt, &
1134 : dgnum, dgnumwet, &
1135 89784 : sulfeq )
1136 :
1137 89784 : if (ndx_h2so4 > 0) then
1138 139513896 : del_h2so4_aeruptk(1:ncol,:) = vmr(1:ncol,:,ndx_h2so4) - del_h2so4_aeruptk(1:ncol,:)
1139 : endif
1140 :
1141 89784 : call t_stopf('modal_gas-aer_exchng')
1142 :
1143 89784 : call t_startf('modal_nucl')
1144 :
1145 : ! do aerosol nucleation (new particle formation)
1146 : call modal_aero_newnuc_sub( &
1147 : lchnk, ncol, nstep, &
1148 : loffset, delt, &
1149 : tfld, pmid, pdel, &
1150 : zm, pblh, &
1151 : qh2o, cldfr, &
1152 : vmr, &
1153 32777424 : del_h2so4_gasprod, del_h2so4_aeruptk )
1154 :
1155 89784 : call t_stopf('modal_nucl')
1156 :
1157 89784 : call t_startf('modal_coag')
1158 :
1159 : ! do aerosol coagulation
1160 : call modal_aero_coag_sub( &
1161 : lchnk, ncol, nstep, &
1162 : loffset, delt, &
1163 : tfld, pmid, pdel, &
1164 : vmr, &
1165 : dgnum, dgnumwet, &
1166 89784 : wetdens )
1167 :
1168 89784 : call t_stopf('modal_coag')
1169 :
1170 89784 : call vmr2qqcw( lchnk, vmrcw, mbar, ncol, loffset, pbuf )
1171 :
1172 : ! diagnostics for cloud-borne aerosols...
1173 3770928 : do n = 1,pcnst
1174 3681144 : fldcw => qqcw_get_field(pbuf,n,lchnk,errorhandle=.true.)
1175 3770928 : if(associated(fldcw)) then
1176 1705896 : call outfld( cnst_name_cw(n), fldcw(:,:), pcols, lchnk )
1177 : endif
1178 : end do
1179 : !
1180 : ! ... put missing NH3 into NH4
1181 : !
1182 89784 : if ( nh3_ndx > 0 .and. nh4_ndx > 0 ) then
1183 0 : vmr(1:ncol,:,nh4_ndx) = vmr(1:ncol,:,nh4_ndx) + (nh3_beg-vmr(1:ncol,:,nh3_ndx))
1184 0 : vmr(1:ncol,:,nh4_ndx) = max(0._r8,vmr(1:ncol,:,nh4_ndx))
1185 : end if
1186 :
1187 179568 : end subroutine aero_model_gasaerexch
1188 :
1189 : !=============================================================================
1190 : !=============================================================================
1191 89784 : subroutine aero_model_emissions( state, cam_in )
1192 89784 : use seasalt_model, only: seasalt_emis, seasalt_names, seasalt_indices, seasalt_active,seasalt_nbin
1193 : use dust_model, only: dust_emis, dust_names, dust_indices, dust_active,dust_nbin, dust_nnum
1194 : use physics_types, only: physics_state
1195 :
1196 : ! Arguments:
1197 :
1198 : type(physics_state), intent(in) :: state ! Physics state variables
1199 : type(cam_in_t), intent(inout) :: cam_in ! import state
1200 :
1201 : ! local vars
1202 :
1203 : integer :: lchnk, ncol
1204 : integer :: m, mm
1205 : real(r8) :: soil_erod_tmp(pcols)
1206 : real(r8) :: sflx(pcols) ! accumulate over all bins for output
1207 : real(r8) :: u10cubed(pcols)
1208 : real (r8), parameter :: z0=0.0001_r8 ! m roughness length over oceans--from ocean model
1209 :
1210 89784 : lchnk = state%lchnk
1211 89784 : ncol = state%ncol
1212 :
1213 89784 : if (dust_active) then
1214 :
1215 89784 : call dust_emis( ncol, lchnk, cam_in%dstflx, cam_in%cflx, soil_erod_tmp )
1216 :
1217 : ! some dust emis diagnostics ...
1218 89784 : sflx(:)=0._r8
1219 628488 : do m=1,dust_nbin+dust_nnum
1220 538704 : mm = dust_indices(m)
1221 4766904 : if (m<=dust_nbin) sflx(:ncol)=sflx(:ncol)+cam_in%cflx(:ncol,mm)
1222 628488 : call outfld(trim(dust_names(m))//'SF',cam_in%cflx(:,mm),pcols, lchnk)
1223 : enddo
1224 89784 : call outfld('DSTSFMBL',sflx(:),pcols,lchnk)
1225 89784 : call outfld('LND_MBL',soil_erod_tmp(:),pcols, lchnk )
1226 : endif
1227 :
1228 89784 : if (seasalt_active) then
1229 1499184 : u10cubed(:ncol)=sqrt(state%u(:ncol,pver)**2+state%v(:ncol,pver)**2)
1230 : ! move the winds to 10m high from the midpoint of the gridbox:
1231 : ! follows Tie and Seinfeld and Pandis, p.859 with math.
1232 :
1233 1499184 : u10cubed(:ncol)=u10cubed(:ncol)*log(10._r8/z0)/log(state%zm(:ncol,pver)/z0)
1234 :
1235 : ! we need them to the 3.41 power, according to Gong et al., 1997:
1236 1499184 : u10cubed(:ncol)=u10cubed(:ncol)**3.41_r8
1237 :
1238 89784 : sflx(:)=0._r8
1239 :
1240 89784 : call seasalt_emis( u10cubed, cam_in%sst, cam_in%ocnfrac, ncol, cam_in%cflx )
1241 :
1242 359136 : do m=1,seasalt_nbin
1243 269352 : mm = seasalt_indices(m)
1244 4497552 : sflx(:ncol)=sflx(:ncol)+cam_in%cflx(:ncol,mm)
1245 359136 : call outfld(trim(seasalt_names(m))//'SF',cam_in%cflx(:,mm),pcols,lchnk)
1246 : enddo
1247 89784 : call outfld('SSTSFMBL',sflx(:),pcols,lchnk)
1248 : endif
1249 :
1250 89784 : end subroutine aero_model_emissions
1251 :
1252 : !===============================================================================
1253 : ! private methods
1254 :
1255 :
1256 : !=============================================================================
1257 : !=============================================================================
1258 0 : subroutine surf_area_dens( ncol, mmr, pmid, temp, diam, beglev, endlev, sad, reff, sfc )
1259 89784 : use mo_constants, only : pi
1260 : use modal_aero_data, only : nspec_amode, alnsg_amode
1261 :
1262 : ! dummy args
1263 : integer, intent(in) :: ncol
1264 : real(r8), intent(in) :: mmr(:,:,:)
1265 : real(r8), intent(in) :: pmid(:,:)
1266 : real(r8), intent(in) :: temp(:,:)
1267 : real(r8), intent(in) :: diam(:,:,:)
1268 : integer, intent(in) :: beglev(:)
1269 : integer, intent(in) :: endlev(:)
1270 : real(r8), intent(out) :: sad(:,:)
1271 : real(r8), intent(out) :: reff(:,:)
1272 : real(r8),optional, intent(out) :: sfc(:,:,:)
1273 :
1274 : ! local vars
1275 0 : real(r8) :: sad_mode(pcols,pver,ntot_amode),radeff(pcols,pver)
1276 0 : real(r8) :: vol(pcols,pver),vol_mode(pcols,pver,ntot_amode)
1277 : real(r8) :: rho_air
1278 : integer :: i,k,l,m
1279 : real(r8) :: chm_mass, tot_mass
1280 :
1281 : !
1282 : ! Compute surface aero for each mode.
1283 : ! Total over all modes as the surface area for chemical reactions.
1284 : !
1285 :
1286 0 : sad = 0._r8
1287 0 : sad_mode = 0._r8
1288 0 : vol = 0._r8
1289 0 : vol_mode = 0._r8
1290 0 : reff = 0._r8
1291 :
1292 0 : do i = 1,ncol
1293 0 : do k = beglev(i),endlev(i)
1294 0 : rho_air = pmid(i,k)/(temp(i,k)*287.04_r8)
1295 0 : do l=1,ntot_amode
1296 : !
1297 : ! compute a mass weighting of the number
1298 : !
1299 0 : tot_mass = 0._r8
1300 0 : chm_mass = 0._r8
1301 0 : do m=1,nspec_amode(l)
1302 0 : if ( index_tot_mass(l,m) > 0 ) &
1303 0 : tot_mass = tot_mass + mmr(i,k,index_tot_mass(l,m))
1304 0 : if ( index_chm_mass(l,m) > 0 ) &
1305 0 : chm_mass = chm_mass + mmr(i,k,index_chm_mass(l,m))
1306 : end do
1307 0 : if ( tot_mass > 0._r8 ) then
1308 : ! surface area density
1309 0 : sad_mode(i,k,l) = chm_mass/tot_mass &
1310 0 : * mmr(i,k,num_idx(l))*rho_air*pi*diam(i,k,l)**2._r8 &
1311 0 : * exp(2._r8*alnsg_amode(l)**2._r8) ! m^2/m^3
1312 0 : sad_mode(i,k,l) = 1.e-2_r8 * sad_mode(i,k,l) ! cm^2/cm^3
1313 :
1314 : ! volume calculation, for use in effective radius calculation
1315 : vol_mode(i,k,l) = chm_mass/tot_mass &
1316 0 : * mmr(i,k,num_idx(l))*rho_air*pi/6._r8*diam(i,k,l)**3._r8 &
1317 0 : * exp(4.5_r8*alnsg_amode(l)**2._r8) ! m^3/m^3 = cm^3/cm^3
1318 : else
1319 0 : sad_mode(i,k,l) = 0._r8
1320 0 : vol_mode(i,k,l) = 0._r8
1321 : end if
1322 : end do
1323 0 : sad(i,k) = sum(sad_mode(i,k,:))
1324 0 : vol(i,k) = sum(vol_mode(i,k,:))
1325 0 : reff(i,k) = 3._r8*vol(i,k)/sad(i,k)
1326 :
1327 : enddo
1328 : enddo
1329 :
1330 0 : if (present(sfc)) then
1331 0 : sfc(:,:,:) = sad_mode(:,:,:)
1332 : endif
1333 :
1334 0 : end subroutine surf_area_dens
1335 :
1336 : !===============================================================================
1337 : !===============================================================================
1338 2304 : subroutine modal_aero_bcscavcoef_init
1339 : !-----------------------------------------------------------------------
1340 : !
1341 : ! Purpose:
1342 : ! Computes lookup table for aerosol impaction/interception scavenging rates
1343 : !
1344 : ! Authors: R. Easter
1345 : !
1346 : !-----------------------------------------------------------------------
1347 :
1348 0 : use shr_kind_mod, only: r8 => shr_kind_r8
1349 : use modal_aero_data
1350 : use cam_abortutils, only: endrun
1351 :
1352 : implicit none
1353 :
1354 :
1355 : ! local variables
1356 : integer nnfit_maxd
1357 : parameter (nnfit_maxd=27)
1358 :
1359 : integer i, jgrow, jdens, jpress, jtemp, mode, nnfit
1360 : integer lunerr
1361 :
1362 : real(r8) dg0, dg0_cgs, press, &
1363 : rhodryaero, rhowetaero, rhowetaero_cgs, rmserr, &
1364 : scavratenum, scavratevol, sigmag, &
1365 : temp, wetdiaratio, wetvolratio
1366 : real(r8) aafitnum(1), xxfitnum(1,nnfit_maxd), yyfitnum(nnfit_maxd)
1367 : real(r8) aafitvol(1), xxfitvol(1,nnfit_maxd), yyfitvol(nnfit_maxd)
1368 :
1369 6912 : allocate(scavimptblnum(nimptblgrow_mind:nimptblgrow_maxd, ntot_amode))
1370 4608 : allocate(scavimptblvol(nimptblgrow_mind:nimptblgrow_maxd, ntot_amode))
1371 :
1372 2304 : lunerr = 6
1373 2304 : dlndg_nimptblgrow = log( 1.25_r8 )
1374 :
1375 11520 : modeloop: do mode = 1, ntot_amode
1376 :
1377 9216 : sigmag = sigmag_amode(mode)
1378 :
1379 9216 : rhodryaero = specdens_amode(1,mode)
1380 :
1381 195840 : growloop: do jgrow = nimptblgrow_mind, nimptblgrow_maxd
1382 :
1383 184320 : wetdiaratio = exp( jgrow*dlndg_nimptblgrow )
1384 184320 : dg0 = dgnum_amode(mode)*wetdiaratio
1385 :
1386 : wetvolratio = exp( jgrow*dlndg_nimptblgrow*3._r8 )
1387 184320 : rhowetaero = 1.0_r8 + (rhodryaero-1.0_r8)/wetvolratio
1388 184320 : rhowetaero = min( rhowetaero, rhodryaero )
1389 :
1390 : !
1391 : ! compute impaction scavenging rates at 1 temp-press pair and save
1392 : !
1393 184320 : nnfit = 0
1394 :
1395 184320 : temp = 273.16_r8
1396 184320 : press = 0.75e6_r8 ! dynes/cm2
1397 184320 : rhowetaero = rhodryaero
1398 :
1399 184320 : dg0_cgs = dg0*1.0e2_r8 ! m to cm
1400 184320 : rhowetaero_cgs = rhowetaero*1.0e-3_r8 ! kg/m3 to g/cm3
1401 : call calc_1_impact_rate( &
1402 : dg0_cgs, sigmag, rhowetaero_cgs, temp, press, &
1403 184320 : scavratenum, scavratevol, lunerr )
1404 :
1405 184320 : nnfit = nnfit + 1
1406 : if (nnfit > nnfit_maxd) then
1407 : write(lunerr,9110)
1408 : call endrun()
1409 : end if
1410 : 9110 format( '*** subr. modal_aero_bcscavcoef_init -- nnfit too big' )
1411 :
1412 184320 : xxfitnum(1,nnfit) = 1._r8
1413 184320 : yyfitnum(nnfit) = log( scavratenum )
1414 :
1415 184320 : xxfitvol(1,nnfit) = 1._r8
1416 184320 : yyfitvol(nnfit) = log( scavratevol )
1417 :
1418 : 5900 continue
1419 :
1420 : !
1421 : ! skip mlinfit stuff because scav table no longer has dependencies on
1422 : ! air temp, air press, and particle wet density
1423 : ! just load the log( scavrate--- ) values
1424 : !
1425 : !!
1426 : !! do linear regression
1427 : !! log(scavrate) = a1 + a2*log(wetdens)
1428 : !!
1429 : ! call mlinft( xxfitnum, yyfitnum, aafitnum, nnfit, 1, 1, rmserr )
1430 : ! call mlinft( xxfitvol, yyfitvol, aafitvol, nnfit, 1, 1, rmserr )
1431 : !
1432 : ! scavimptblnum(jgrow,mode) = aafitnum(1)
1433 : ! scavimptblvol(jgrow,mode) = aafitvol(1)
1434 :
1435 184320 : scavimptblnum(jgrow,mode) = yyfitnum(1)
1436 193536 : scavimptblvol(jgrow,mode) = yyfitvol(1)
1437 :
1438 : enddo growloop
1439 : enddo modeloop
1440 2304 : return
1441 2304 : end subroutine modal_aero_bcscavcoef_init
1442 :
1443 : !===============================================================================
1444 : !===============================================================================
1445 897840 : subroutine modal_aero_depvel_part( ncol, t, pmid, ram1, fv, vlc_dry, vlc_trb, vlc_grv, &
1446 : radius_part, density_part, sig_part, moment, lchnk )
1447 :
1448 : ! calculates surface deposition velocity of particles
1449 : ! L. Zhang, S. Gong, J. Padro, and L. Barrie
1450 : ! A size-seggregated particle dry deposition scheme for an atmospheric aerosol module
1451 : ! Atmospheric Environment, 35, 549-560, 2001.
1452 : !
1453 : ! Authors: X. Liu
1454 :
1455 : !
1456 : ! !USES
1457 : !
1458 2304 : use physconst, only: pi,boltz, gravit, rair
1459 : use mo_drydep, only: n_land_type, fraction_landuse
1460 :
1461 : ! !ARGUMENTS:
1462 : !
1463 : implicit none
1464 : !
1465 : real(r8), intent(in) :: t(pcols,pver) !atm temperature (K)
1466 : real(r8), intent(in) :: pmid(pcols,pver) !atm pressure (Pa)
1467 : real(r8), intent(in) :: fv(pcols) !friction velocity (m/s)
1468 : real(r8), intent(in) :: ram1(pcols) !aerodynamical resistance (s/m)
1469 : real(r8), intent(in) :: radius_part(pcols,pver) ! mean (volume/number) particle radius (m)
1470 : real(r8), intent(in) :: density_part(pcols,pver) ! density of particle material (kg/m3)
1471 : real(r8), intent(in) :: sig_part(pcols,pver) ! geometric standard deviation of particles
1472 : integer, intent(in) :: moment ! moment of size distribution (0 for number, 2 for surface area, 3 for volume)
1473 : integer, intent(in) :: ncol
1474 : integer, intent(in) :: lchnk
1475 :
1476 : real(r8), intent(out) :: vlc_trb(pcols) !Turbulent deposn velocity (m/s)
1477 : real(r8), intent(out) :: vlc_grv(pcols,pver) !grav deposn velocity (m/s)
1478 : real(r8), intent(out) :: vlc_dry(pcols,pver) !dry deposn velocity (m/s)
1479 : !------------------------------------------------------------------------
1480 :
1481 : !------------------------------------------------------------------------
1482 : ! Local Variables
1483 : integer :: m,i,k,ix !indices
1484 : real(r8) :: rho !atm density (kg/m**3)
1485 : real(r8) :: vsc_dyn_atm(pcols,pver) ![kg m-1 s-1] Dynamic viscosity of air
1486 : real(r8) :: vsc_knm_atm(pcols,pver) ![m2 s-1] Kinematic viscosity of atmosphere
1487 : real(r8) :: shm_nbr ![frc] Schmidt number
1488 : real(r8) :: stk_nbr ![frc] Stokes number
1489 : real(r8) :: mfp_atm(pcols,pver) ![m] Mean free path of air
1490 : real(r8) :: dff_aer ![m2 s-1] Brownian diffusivity of particle
1491 : real(r8) :: slp_crc(pcols,pver) ![frc] Slip correction factor
1492 : real(r8) :: rss_trb ![s m-1] Resistance to turbulent deposition
1493 : real(r8) :: rss_lmn ![s m-1] Quasi-laminar layer resistance
1494 : real(r8) :: brownian ! collection efficiency for Browning diffusion
1495 : real(r8) :: impaction ! collection efficiency for impaction
1496 : real(r8) :: interception ! collection efficiency for interception
1497 : real(r8) :: stickfrac ! fraction of particles sticking to surface
1498 : real(r8) :: radius_moment(pcols,pver) ! median radius (m) for moment
1499 : real(r8) :: lnsig ! ln(sig_part)
1500 : real(r8) :: dispersion ! accounts for influence of size dist dispersion on bulk settling velocity
1501 : ! assuming radius_part is number mode radius * exp(1.5 ln(sigma))
1502 :
1503 : integer :: lt
1504 : real(r8) :: lnd_frc
1505 : real(r8) :: wrk1, wrk2, wrk3
1506 :
1507 : ! constants
1508 : real(r8) gamma(11) ! exponent of schmidt number
1509 : ! data gamma/0.54d+00, 0.56d+00, 0.57d+00, 0.54d+00, 0.54d+00, &
1510 : ! 0.56d+00, 0.54d+00, 0.54d+00, 0.54d+00, 0.56d+00, &
1511 : ! 0.50d+00/
1512 : data gamma/0.56e+00_r8, 0.54e+00_r8, 0.54e+00_r8, 0.56e+00_r8, 0.56e+00_r8, &
1513 : 0.56e+00_r8, 0.50e+00_r8, 0.54e+00_r8, 0.54e+00_r8, 0.54e+00_r8, &
1514 : 0.54e+00_r8/
1515 : save gamma
1516 :
1517 : real(r8) alpha(11) ! parameter for impaction
1518 : ! data alpha/50.00d+00, 0.95d+00, 0.80d+00, 1.20d+00, 1.30d+00, &
1519 : ! 0.80d+00, 50.00d+00, 50.00d+00, 2.00d+00, 1.50d+00, &
1520 : ! 100.00d+00/
1521 : data alpha/1.50e+00_r8, 1.20e+00_r8, 1.20e+00_r8, 0.80e+00_r8, 1.00e+00_r8, &
1522 : 0.80e+00_r8, 100.00e+00_r8, 50.00e+00_r8, 2.00e+00_r8, 1.20e+00_r8, &
1523 : 50.00e+00_r8/
1524 : save alpha
1525 :
1526 : real(r8) radius_collector(11) ! radius (m) of surface collectors
1527 : ! data radius_collector/-1.00d+00, 5.10d-03, 3.50d-03, 3.20d-03, 10.00d-03, &
1528 : ! 5.00d-03, -1.00d+00, -1.00d+00, 10.00d-03, 10.00d-03, &
1529 : ! -1.00d+00/
1530 : data radius_collector/10.00e-03_r8, 3.50e-03_r8, 3.50e-03_r8, 5.10e-03_r8, 2.00e-03_r8, &
1531 : 5.00e-03_r8, -1.00e+00_r8, -1.00e+00_r8, 10.00e-03_r8, 3.50e-03_r8, &
1532 : -1.00e+00_r8/
1533 : save radius_collector
1534 :
1535 : integer :: iwet(11) ! flag for wet surface = 1, otherwise = -1
1536 : ! data iwet/1, -1, -1, -1, -1, &
1537 : ! -1, -1, -1, 1, -1, &
1538 : ! 1/
1539 : data iwet/-1, -1, -1, -1, -1, &
1540 : -1, 1, -1, 1, -1, &
1541 : -1/
1542 : save iwet
1543 :
1544 :
1545 897840 : vlc_trb = 0._r8
1546 897840 : vlc_grv = 0._r8
1547 897840 : vlc_dry = 0._r8
1548 :
1549 : !------------------------------------------------------------------------
1550 84396960 : do k=top_lev,pver ! radius_part is not defined above top_lev
1551 1395138960 : do i=1,ncol
1552 :
1553 1310742000 : lnsig = log(sig_part(i,k))
1554 : ! use a maximum radius of 50 microns when calculating deposition velocity
1555 : radius_moment(i,k) = min(50.0e-6_r8,radius_part(i,k))* &
1556 1310742000 : exp((float(moment)-1.5_r8)*lnsig*lnsig)
1557 1310742000 : dispersion = exp(2._r8*lnsig*lnsig)
1558 :
1559 1310742000 : rho=pmid(i,k)/rair/t(i,k)
1560 :
1561 : ! Quasi-laminar layer resistance: call rss_lmn_get
1562 : ! Size-independent thermokinetic properties
1563 : vsc_dyn_atm(i,k) = 1.72e-5_r8 * ((t(i,k)/273.0_r8)**1.5_r8) * 393.0_r8 / &
1564 1310742000 : (t(i,k)+120.0_r8) ![kg m-1 s-1] RoY94 p. 102
1565 : mfp_atm(i,k) = 2.0_r8 * vsc_dyn_atm(i,k) / & ![m] SeP97 p. 455
1566 1310742000 : (pmid(i,k)*sqrt(8.0_r8/(pi*rair*t(i,k))))
1567 1310742000 : vsc_knm_atm(i,k) = vsc_dyn_atm(i,k) / rho ![m2 s-1] Kinematic viscosity of air
1568 :
1569 : slp_crc(i,k) = 1.0_r8 + mfp_atm(i,k) * &
1570 : (1.257_r8+0.4_r8*exp(-1.1_r8*radius_moment(i,k)/(mfp_atm(i,k)))) / &
1571 1310742000 : radius_moment(i,k) ![frc] Slip correction factor SeP97 p. 464
1572 : vlc_grv(i,k) = (4.0_r8/18.0_r8) * radius_moment(i,k)*radius_moment(i,k)*density_part(i,k)* &
1573 1310742000 : gravit*slp_crc(i,k) / vsc_dyn_atm(i,k) ![m s-1] Stokes' settling velocity SeP97 p. 466
1574 1310742000 : vlc_grv(i,k) = vlc_grv(i,k) * dispersion
1575 :
1576 1394241120 : vlc_dry(i,k)=vlc_grv(i,k)
1577 : enddo
1578 : enddo
1579 897840 : k=pver ! only look at bottom level for next part
1580 14991840 : do i=1,ncol
1581 14094000 : dff_aer = boltz * t(i,k) * slp_crc(i,k) / & ![m2 s-1]
1582 14094000 : (6.0_r8*pi*vsc_dyn_atm(i,k)*radius_moment(i,k)) !SeP97 p.474
1583 14094000 : shm_nbr = vsc_knm_atm(i,k) / dff_aer ![frc] SeP97 p.972
1584 :
1585 14094000 : wrk2 = 0._r8
1586 14094000 : wrk3 = 0._r8
1587 169128000 : do lt = 1,n_land_type
1588 155034000 : lnd_frc = fraction_landuse(i,lt,lchnk)
1589 169128000 : if ( lnd_frc /= 0._r8 ) then
1590 34407630 : brownian = shm_nbr**(-gamma(lt))
1591 34407630 : if (radius_collector(lt) > 0.0_r8) then
1592 : ! vegetated surface
1593 15669570 : stk_nbr = vlc_grv(i,k) * fv(i) / (gravit*radius_collector(lt))
1594 15669570 : interception = 2.0_r8*(radius_moment(i,k)/radius_collector(lt))**2.0_r8
1595 : else
1596 : ! non-vegetated surface
1597 18738060 : stk_nbr = vlc_grv(i,k) * fv(i) * fv(i) / (gravit*vsc_knm_atm(i,k)) ![frc] SeP97 p.965
1598 18738060 : interception = 0.0_r8
1599 : endif
1600 34407630 : impaction = (stk_nbr/(alpha(lt)+stk_nbr))**2.0_r8
1601 :
1602 34407630 : if (iwet(lt) > 0) then
1603 : stickfrac = 1.0_r8
1604 : else
1605 21321090 : stickfrac = exp(-sqrt(stk_nbr))
1606 21321090 : if (stickfrac < 1.0e-10_r8) stickfrac = 1.0e-10_r8
1607 : endif
1608 34407630 : rss_lmn = 1.0_r8 / (3.0_r8 * fv(i) * stickfrac * (brownian+interception+impaction))
1609 34407630 : rss_trb = ram1(i) + rss_lmn + ram1(i)*rss_lmn*vlc_grv(i,k)
1610 :
1611 34407630 : wrk1 = 1.0_r8 / rss_trb
1612 34407630 : wrk2 = wrk2 + lnd_frc*( wrk1 )
1613 34407630 : wrk3 = wrk3 + lnd_frc*( wrk1 + vlc_grv(i,k) )
1614 : endif
1615 : enddo ! n_land_type
1616 14094000 : vlc_trb(i) = wrk2
1617 14991840 : vlc_dry(i,k) = wrk3
1618 : enddo !ncol
1619 :
1620 897840 : return
1621 897840 : end subroutine modal_aero_depvel_part
1622 :
1623 : !===============================================================================
1624 : subroutine modal_aero_bcscavcoef_get( m, ncol, isprx, dgn_awet, scavcoefnum, scavcoefvol )
1625 :
1626 897840 : use modal_aero_data
1627 : !-----------------------------------------------------------------------
1628 : implicit none
1629 :
1630 : integer,intent(in) :: m, ncol
1631 : logical,intent(in):: isprx(pcols,pver)
1632 : real(r8), intent(in) :: dgn_awet(pcols,pver,ntot_amode)
1633 : real(r8), intent(out) :: scavcoefnum(pcols,pver), scavcoefvol(pcols,pver)
1634 :
1635 : integer i, k, jgrow
1636 : real(r8) dumdgratio, xgrow, dumfhi, dumflo, scavimpvol, scavimpnum
1637 :
1638 :
1639 : do k = 1, pver
1640 : do i = 1, ncol
1641 :
1642 : ! do only if no precip
1643 : if ( isprx(i,k) ) then
1644 : !
1645 : ! interpolate table values using log of (actual-wet-size)/(base-dry-size)
1646 :
1647 : dumdgratio = dgn_awet(i,k,m)/dgnum_amode(m)
1648 :
1649 : if ((dumdgratio >= 0.99_r8) .and. (dumdgratio <= 1.01_r8)) then
1650 : scavimpvol = scavimptblvol(0,m)
1651 : scavimpnum = scavimptblnum(0,m)
1652 : else
1653 : xgrow = log( dumdgratio ) / dlndg_nimptblgrow
1654 : jgrow = int( xgrow )
1655 : if (xgrow < 0._r8) jgrow = jgrow - 1
1656 : if (jgrow < nimptblgrow_mind) then
1657 : jgrow = nimptblgrow_mind
1658 : xgrow = jgrow
1659 : else
1660 : jgrow = min( jgrow, nimptblgrow_maxd-1 )
1661 : end if
1662 :
1663 : dumfhi = xgrow - jgrow
1664 : dumflo = 1._r8 - dumfhi
1665 :
1666 : scavimpvol = dumflo*scavimptblvol(jgrow,m) + &
1667 : dumfhi*scavimptblvol(jgrow+1,m)
1668 : scavimpnum = dumflo*scavimptblnum(jgrow,m) + &
1669 : dumfhi*scavimptblnum(jgrow+1,m)
1670 :
1671 : end if
1672 :
1673 : ! impaction scavenging removal amount for volume
1674 : scavcoefvol(i,k) = exp( scavimpvol )
1675 : ! impaction scavenging removal amount to number
1676 : scavcoefnum(i,k) = exp( scavimpnum )
1677 :
1678 : ! scavcoef = impaction scav rate (1/h) for precip = 1 mm/h
1679 : ! scavcoef = impaction scav rate (1/s) for precip = pfx_inrain
1680 : ! (scavcoef/3600) = impaction scav rate (1/s) for precip = 1 mm/h
1681 : ! (pfx_inrain*3600) = in-rain-area precip rate (mm/h)
1682 : ! impactrate = (scavcoef/3600) * (pfx_inrain*3600)
1683 : else
1684 : scavcoefvol(i,k) = 0._r8
1685 : scavcoefnum(i,k) = 0._r8
1686 : end if
1687 :
1688 : end do
1689 : end do
1690 :
1691 : return
1692 : end subroutine modal_aero_bcscavcoef_get
1693 :
1694 : !===============================================================================
1695 184320 : subroutine calc_1_impact_rate( &
1696 : dg0, sigmag, rhoaero, temp, press, &
1697 : scavratenum, scavratevol, lunerr )
1698 : !
1699 : ! routine computes a single impaction scavenging rate
1700 : ! for precipitation rate of 1 mm/h
1701 : !
1702 : ! dg0 = geometric mean diameter of aerosol number size distrib. (cm)
1703 : ! sigmag = geometric standard deviation of size distrib.
1704 : ! rhoaero = density of aerosol particles (g/cm^3)
1705 : ! temp = temperature (K)
1706 : ! press = pressure (dyne/cm^2)
1707 : ! scavratenum = number scavenging rate (1/h)
1708 : ! scavratevol = volume or mass scavenging rate (1/h)
1709 : ! lunerr = logical unit for error message
1710 : !
1711 : use shr_kind_mod, only: r8 => shr_kind_r8
1712 : use mo_constants, only: boltz_cgs, pi, rhowater => rhoh2o_cgs, &
1713 : gravity => gravity_cgs, rgas => rgas_cgs
1714 :
1715 : implicit none
1716 :
1717 : ! subr. parameters
1718 : integer lunerr
1719 : real(r8) dg0, sigmag, rhoaero, temp, press, scavratenum, scavratevol
1720 :
1721 : ! local variables
1722 : integer nrainsvmax
1723 : parameter (nrainsvmax=50)
1724 : real(r8) rrainsv(nrainsvmax), xnumrainsv(nrainsvmax),&
1725 : vfallrainsv(nrainsvmax)
1726 :
1727 : integer naerosvmax
1728 : parameter (naerosvmax=51)
1729 : real(r8) aaerosv(naerosvmax), &
1730 : ynumaerosv(naerosvmax), yvolaerosv(naerosvmax)
1731 :
1732 : integer i, ja, jr, na, nr
1733 : real(r8) a, aerodiffus, aeromass, ag0, airdynvisc, airkinvisc
1734 : real(r8) anumsum, avolsum, cair, chi
1735 : real(r8) d, dr, dum, dumfuchs, dx
1736 : real(r8) ebrown, eimpact, eintercept, etotal, freepath
1737 : real(r8) precip, precipmmhr, precipsum
1738 : real(r8) r, rainsweepout, reynolds, rhi, rhoair, rlo, rnumsum
1739 : real(r8) scavsumnum, scavsumnumbb
1740 : real(r8) scavsumvol, scavsumvolbb
1741 : real(r8) schmidt, sqrtreynolds, sstar, stokes, sx
1742 : real(r8) taurelax, vfall, vfallstp
1743 : real(r8) x, xg0, xg3, xhi, xlo, xmuwaterair
1744 :
1745 :
1746 184320 : rlo = .005_r8
1747 184320 : rhi = .250_r8
1748 184320 : dr = 0.005_r8
1749 184320 : nr = 1 + nint( (rhi-rlo)/dr )
1750 : if (nr > nrainsvmax) then
1751 : write(lunerr,9110)
1752 : call endrun()
1753 : end if
1754 :
1755 : 9110 format( '*** subr. calc_1_impact_rate -- nr > nrainsvmax' )
1756 :
1757 184320 : precipmmhr = 1.0_r8
1758 184320 : precip = precipmmhr/36000._r8
1759 :
1760 184320 : ag0 = dg0/2._r8
1761 184320 : sx = log( sigmag )
1762 184320 : xg0 = log( ag0 )
1763 184320 : xg3 = xg0 + 3._r8*sx*sx
1764 :
1765 184320 : xlo = xg3 - 4._r8*sx
1766 184320 : xhi = xg3 + 4._r8*sx
1767 184320 : dx = 0.2_r8*sx
1768 :
1769 184320 : dx = max( 0.2_r8*sx, 0.01_r8 )
1770 184320 : xlo = xg3 - max( 4._r8*sx, 2._r8*dx )
1771 184320 : xhi = xg3 + max( 4._r8*sx, 2._r8*dx )
1772 :
1773 184320 : na = 1 + nint( (xhi-xlo)/dx )
1774 184320 : if (na > naerosvmax) then
1775 0 : write(lunerr,9120)
1776 0 : call endrun()
1777 : end if
1778 :
1779 : 9120 format( '*** subr. calc_1_impact_rate -- na > naerosvmax' )
1780 :
1781 : ! air molar density
1782 184320 : cair = press/(rgas*temp)
1783 : ! air mass density
1784 184320 : rhoair = 28.966_r8*cair
1785 : ! molecular freepath
1786 184320 : freepath = 2.8052e-10_r8/cair
1787 : ! air dynamic viscosity
1788 : airdynvisc = 1.8325e-4_r8 * (416.16_r8/(temp+120._r8)) * &
1789 184320 : ((temp/296.16_r8)**1.5_r8)
1790 : ! air kinemaic viscosity
1791 184320 : airkinvisc = airdynvisc/rhoair
1792 : ! ratio of water viscosity to air viscosity (from Slinn)
1793 184320 : xmuwaterair = 60.0_r8
1794 :
1795 : !
1796 : ! compute rain drop number concentrations
1797 : ! rrainsv = raindrop radius (cm)
1798 : ! xnumrainsv = raindrop number concentration (#/cm^3)
1799 : ! (number in the bin, not number density)
1800 : ! vfallrainsv = fall velocity (cm/s)
1801 : !
1802 184320 : precipsum = 0._r8
1803 9400320 : do i = 1, nr
1804 9216000 : r = rlo + (i-1)*dr
1805 9216000 : rrainsv(i) = r
1806 9216000 : xnumrainsv(i) = exp( -r/2.7e-2_r8 )
1807 :
1808 9216000 : d = 2._r8*r
1809 9216000 : if (d <= 0.007_r8) then
1810 0 : vfallstp = 2.88e5_r8 * d**2._r8
1811 9216000 : else if (d <= 0.025_r8) then
1812 368640 : vfallstp = 2.8008e4_r8 * d**1.528_r8
1813 8847360 : else if (d <= 0.1_r8) then
1814 1474560 : vfallstp = 4104.9_r8 * d**1.008_r8
1815 7372800 : else if (d <= 0.25_r8) then
1816 2764800 : vfallstp = 1812.1_r8 * d**0.638_r8
1817 : else
1818 4608000 : vfallstp = 1069.8_r8 * d**0.235_r8
1819 : end if
1820 :
1821 9216000 : vfall = vfallstp * sqrt(1.204e-3_r8/rhoair)
1822 9216000 : vfallrainsv(i) = vfall
1823 9400320 : precipsum = precipsum + vfall*(r**3)*xnumrainsv(i)
1824 : end do
1825 184320 : precipsum = precipsum*pi*1.333333_r8
1826 :
1827 184320 : rnumsum = 0._r8
1828 9400320 : do i = 1, nr
1829 9216000 : xnumrainsv(i) = xnumrainsv(i)*(precip/precipsum)
1830 9400320 : rnumsum = rnumsum + xnumrainsv(i)
1831 : end do
1832 :
1833 : !
1834 : ! compute aerosol concentrations
1835 : ! aaerosv = particle radius (cm)
1836 : ! fnumaerosv = fraction of total number in the bin (--)
1837 : ! fvolaerosv = fraction of total volume in the bin (--)
1838 : !
1839 : anumsum = 0._r8
1840 : avolsum = 0._r8
1841 7741440 : do i = 1, na
1842 7557120 : x = xlo + (i-1)*dx
1843 7557120 : a = exp( x )
1844 7557120 : aaerosv(i) = a
1845 7557120 : dum = (x - xg0)/sx
1846 7557120 : ynumaerosv(i) = exp( -0.5_r8*dum*dum )
1847 7557120 : yvolaerosv(i) = ynumaerosv(i)*1.3333_r8*pi*a*a*a
1848 7557120 : anumsum = anumsum + ynumaerosv(i)
1849 7741440 : avolsum = avolsum + yvolaerosv(i)
1850 : end do
1851 :
1852 7741440 : do i = 1, na
1853 7557120 : ynumaerosv(i) = ynumaerosv(i)/anumsum
1854 7741440 : yvolaerosv(i) = yvolaerosv(i)/avolsum
1855 : end do
1856 :
1857 :
1858 : !
1859 : ! compute scavenging
1860 : !
1861 : scavsumnum = 0._r8
1862 : scavsumvol = 0._r8
1863 : !
1864 : ! outer loop for rain drop radius
1865 : !
1866 9400320 : jr_loop: do jr = 1, nr
1867 :
1868 9216000 : r = rrainsv(jr)
1869 9216000 : vfall = vfallrainsv(jr)
1870 :
1871 9216000 : reynolds = r * vfall / airkinvisc
1872 9216000 : sqrtreynolds = sqrt( reynolds )
1873 :
1874 : !
1875 : ! inner loop for aerosol particle radius
1876 : !
1877 9216000 : scavsumnumbb = 0._r8
1878 9216000 : scavsumvolbb = 0._r8
1879 :
1880 387072000 : ja_loop: do ja = 1, na
1881 :
1882 377856000 : a = aaerosv(ja)
1883 :
1884 377856000 : chi = a/r
1885 :
1886 377856000 : dum = freepath/a
1887 377856000 : dumfuchs = 1._r8 + 1.246_r8*dum + 0.42_r8*dum*exp(-0.87_r8/dum)
1888 377856000 : taurelax = 2._r8*rhoaero*a*a*dumfuchs/(9._r8*rhoair*airkinvisc)
1889 :
1890 377856000 : aeromass = 4._r8*pi*a*a*a*rhoaero/3._r8
1891 377856000 : aerodiffus = boltz_cgs*temp*taurelax/aeromass
1892 :
1893 377856000 : schmidt = airkinvisc/aerodiffus
1894 377856000 : stokes = vfall*taurelax/r
1895 :
1896 : ebrown = 4._r8*(1._r8 + 0.4_r8*sqrtreynolds*(schmidt**0.3333333_r8)) / &
1897 377856000 : (reynolds*schmidt)
1898 :
1899 : dum = (1._r8 + 2._r8*xmuwaterair*chi) / &
1900 377856000 : (1._r8 + xmuwaterair/sqrtreynolds)
1901 377856000 : eintercept = 4._r8*chi*(chi + dum)
1902 :
1903 377856000 : dum = log( 1._r8 + reynolds )
1904 377856000 : sstar = (1.2_r8 + dum/12._r8) / (1._r8 + dum)
1905 377856000 : eimpact = 0._r8
1906 377856000 : if (stokes > sstar) then
1907 92858112 : dum = stokes - sstar
1908 92858112 : eimpact = (dum/(dum+0.6666667_r8)) ** 1.5_r8
1909 : end if
1910 :
1911 377856000 : etotal = ebrown + eintercept + eimpact
1912 377856000 : etotal = min( etotal, 1.0_r8 )
1913 :
1914 377856000 : rainsweepout = xnumrainsv(jr)*4._r8*pi*r*r*vfall
1915 :
1916 377856000 : scavsumnumbb = scavsumnumbb + rainsweepout*etotal*ynumaerosv(ja)
1917 387072000 : scavsumvolbb = scavsumvolbb + rainsweepout*etotal*yvolaerosv(ja)
1918 :
1919 : enddo ja_loop
1920 :
1921 9216000 : scavsumnum = scavsumnum + scavsumnumbb
1922 9400320 : scavsumvol = scavsumvol + scavsumvolbb
1923 :
1924 : enddo jr_loop
1925 :
1926 184320 : scavratenum = scavsumnum*3600._r8
1927 184320 : scavratevol = scavsumvol*3600._r8
1928 :
1929 184320 : return
1930 : end subroutine calc_1_impact_rate
1931 :
1932 : !=============================================================================
1933 : !=============================================================================
1934 179568 : subroutine qqcw2vmr(lchnk, vmr, mbar, ncol, im, pbuf)
1935 : use modal_aero_data, only : qqcw_get_field
1936 : use physics_buffer, only : physics_buffer_desc
1937 : !-----------------------------------------------------------------
1938 : ! ... Xfrom from mass to volume mixing ratio
1939 : !-----------------------------------------------------------------
1940 :
1941 : use chem_mods, only : adv_mass, gas_pcnst
1942 :
1943 : implicit none
1944 :
1945 : !-----------------------------------------------------------------
1946 : ! ... Dummy args
1947 : !-----------------------------------------------------------------
1948 : integer, intent(in) :: lchnk, ncol, im
1949 : real(r8), intent(in) :: mbar(ncol,pver)
1950 : real(r8), intent(inout) :: vmr(ncol,pver,gas_pcnst)
1951 : type(physics_buffer_desc), pointer :: pbuf(:)
1952 :
1953 : !-----------------------------------------------------------------
1954 : ! ... Local variables
1955 : !-----------------------------------------------------------------
1956 : integer :: k, m
1957 89784 : real(r8), pointer :: fldcw(:,:)
1958 :
1959 2873088 : do m=1,gas_pcnst
1960 2873088 : if( adv_mass(m) /= 0._r8 ) then
1961 2783304 : fldcw => qqcw_get_field(pbuf, m+im,lchnk,errorhandle=.true.)
1962 2783304 : if(associated(fldcw)) then
1963 160354224 : do k=1,pver
1964 2650764024 : vmr(:ncol,k,m) = mbar(:ncol,k) * fldcw(:ncol,k) / adv_mass(m)
1965 : end do
1966 : else
1967 1674166752 : vmr(:,:,m) = 0.0_r8
1968 : end if
1969 : end if
1970 : end do
1971 179568 : end subroutine qqcw2vmr
1972 :
1973 :
1974 : !=============================================================================
1975 : !=============================================================================
1976 179568 : subroutine vmr2qqcw( lchnk, vmr, mbar, ncol, im, pbuf )
1977 : !-----------------------------------------------------------------
1978 : ! ... Xfrom from volume to mass mixing ratio
1979 : !-----------------------------------------------------------------
1980 :
1981 89784 : use m_spc_id
1982 : use chem_mods, only : adv_mass, gas_pcnst
1983 : use modal_aero_data, only : qqcw_get_field
1984 : use physics_buffer, only : physics_buffer_desc
1985 :
1986 : implicit none
1987 :
1988 : !-----------------------------------------------------------------
1989 : ! ... Dummy args
1990 : !-----------------------------------------------------------------
1991 : integer, intent(in) :: lchnk, ncol, im
1992 : real(r8), intent(in) :: mbar(ncol,pver)
1993 : real(r8), intent(in) :: vmr(ncol,pver,gas_pcnst)
1994 : type(physics_buffer_desc), pointer :: pbuf(:)
1995 :
1996 : !-----------------------------------------------------------------
1997 : ! ... Local variables
1998 : !-----------------------------------------------------------------
1999 : integer :: k, m
2000 89784 : real(r8), pointer :: fldcw(:,:)
2001 : !-----------------------------------------------------------------
2002 : ! ... The non-group species
2003 : !-----------------------------------------------------------------
2004 2873088 : do m = 1,gas_pcnst
2005 2783304 : fldcw => qqcw_get_field(pbuf, m+im,lchnk,errorhandle=.true.)
2006 2873088 : if( adv_mass(m) /= 0._r8 .and. associated(fldcw)) then
2007 160354224 : do k = 1,pver
2008 2650764024 : fldcw(:ncol,k) = adv_mass(m) * vmr(:ncol,k,m) / mbar(:ncol,k)
2009 : end do
2010 : end if
2011 : end do
2012 :
2013 179568 : end subroutine vmr2qqcw
2014 :
2015 0 : function get_dlndg_nimptblgrow() result (dlndg_nimptblgrow_ret)
2016 : real(r8) :: dlndg_nimptblgrow_ret
2017 0 : dlndg_nimptblgrow_ret = dlndg_nimptblgrow
2018 89784 : end function get_dlndg_nimptblgrow
2019 :
2020 0 : function get_scavimptblvol() result (scavimptblvol_ret)
2021 : real(r8) :: scavimptblvol_ret(nimptblgrow_mind:nimptblgrow_maxd, ntot_amode)
2022 0 : scavimptblvol_ret = scavimptblvol
2023 0 : end function get_scavimptblvol
2024 :
2025 0 : function get_scavimptblnum() result (scavimptblnum_ret)
2026 : real(r8) :: scavimptblnum_ret(nimptblgrow_mind:nimptblgrow_maxd, ntot_amode)
2027 0 : scavimptblnum_ret = scavimptblnum
2028 0 : end function get_scavimptblnum
2029 :
2030 : end module aero_model
|
