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