Line data Source code
1 : !===============================================================================
2 : ! Bulk 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 cam_abortutils, only: endrun
9 : use cam_logfile, only: iulog
10 : use perf_mod, only: t_startf, t_stopf
11 : use camsrfexch, only: cam_in_t, cam_out_t
12 : use aerodep_flx, only: aerodep_flx_prescribed
13 : use physics_types, only: physics_state, physics_ptend, physics_ptend_init
14 : use physics_buffer, only: physics_buffer_desc
15 : use physconst, only: gravit, rair
16 : use dust_model, only: dust_active, dust_names, dust_nbin
17 : use seasalt_model, only: sslt_active=>seasalt_active, seasalt_names, seasalt_nbin
18 : use spmd_utils, only: masterproc
19 : use physics_buffer, only: pbuf_get_field, pbuf_get_index
20 : use cam_history, only: outfld
21 : use infnan, only: nan, assignment(=)
22 :
23 : implicit none
24 : private
25 :
26 : public :: aero_model_readnl
27 : public :: aero_model_register
28 : public :: aero_model_init
29 : public :: aero_model_gasaerexch ! create, grow, change, and shrink aerosols.
30 : public :: aero_model_drydep ! aerosol dry deposition and sediment
31 : public :: aero_model_wetdep ! aerosol wet removal
32 : public :: aero_model_emissions ! aerosol emissions
33 : public :: aero_model_surfarea ! tropospheric aerosol wet surface area for chemistry
34 : public :: aero_model_strat_surfarea ! stub
35 :
36 : public :: wetdep_lq
37 :
38 : ! Misc private data
39 :
40 : integer :: so4_ndx, cb2_ndx, oc2_ndx, nit_ndx
41 : integer :: soa_ndx, soai_ndx, soam_ndx, soab_ndx, soat_ndx, soax_ndx
42 :
43 : ! Namelist variables
44 : character(len=16) :: wetdep_list(pcnst) = ' '
45 : character(len=16) :: drydep_list(pcnst) = ' '
46 :
47 : integer :: ndrydep = 0
48 : integer,allocatable :: drydep_indices(:)
49 : integer :: nwetdep = 0
50 : integer,allocatable :: wetdep_indices(:)
51 : logical :: drydep_lq(pcnst)
52 : logical, protected :: wetdep_lq(pcnst)
53 :
54 : integer :: fracis_idx = 0
55 :
56 : real(r8) :: aer_sol_facti(pcnst) ! in-cloud solubility factor
57 : real(r8) :: aer_sol_factb(pcnst) ! below-cloud solubility factor
58 : real(r8) :: aer_scav_coef(pcnst)
59 :
60 : contains
61 :
62 : !=============================================================================
63 : ! reads aerosol namelist options
64 : !=============================================================================
65 0 : subroutine aero_model_readnl(nlfile)
66 :
67 : use namelist_utils, only: find_group_name
68 : use units, only: getunit, freeunit
69 : use mpishorthand
70 :
71 : character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input
72 :
73 : ! Local variables
74 : integer :: unitn, ierr
75 : character(len=*), parameter :: subname = 'aero_model_readnl'
76 :
77 : ! Namelist variables
78 : character(len=16) :: aer_wetdep_list(pcnst) = ' '
79 : character(len=16) :: aer_drydep_list(pcnst) = ' '
80 :
81 : namelist /aerosol_nl/ aer_wetdep_list, aer_drydep_list
82 : namelist /aerosol_nl/ aer_sol_facti, aer_sol_factb, aer_scav_coef
83 : !-----------------------------------------------------------------------------
84 0 : aer_sol_facti = nan
85 0 : aer_sol_factb = nan
86 0 : aer_scav_coef = nan
87 :
88 : ! Read namelist
89 0 : if (masterproc) then
90 0 : unitn = getunit()
91 0 : open( unitn, file=trim(nlfile), status='old' )
92 0 : call find_group_name(unitn, 'aerosol_nl', status=ierr)
93 0 : if (ierr == 0) then
94 0 : read(unitn, aerosol_nl, iostat=ierr)
95 0 : if (ierr /= 0) then
96 0 : call endrun(subname // ':: ERROR reading namelist')
97 : end if
98 : end if
99 0 : close(unitn)
100 0 : call freeunit(unitn)
101 : end if
102 :
103 : #ifdef SPMD
104 : ! Broadcast namelist variables
105 0 : call mpibcast(aer_wetdep_list, len(aer_wetdep_list(1))*pcnst, mpichar, 0, mpicom)
106 0 : call mpibcast(aer_drydep_list, len(aer_drydep_list(1))*pcnst, mpichar, 0, mpicom)
107 0 : call mpibcast(aer_sol_facti, pcnst, mpir8, 0, mpicom)
108 0 : call mpibcast(aer_sol_factb, pcnst, mpir8, 0, mpicom)
109 0 : call mpibcast(aer_scav_coef, pcnst, mpir8, 0, mpicom)
110 : #endif
111 :
112 0 : wetdep_list = aer_wetdep_list
113 0 : drydep_list = aer_drydep_list
114 :
115 0 : end subroutine aero_model_readnl
116 :
117 : !=============================================================================
118 : !=============================================================================
119 1536 : subroutine aero_model_register()
120 : use mo_setsoa, only : soa_register
121 :
122 1536 : call soa_register()
123 1536 : end subroutine aero_model_register
124 :
125 : !=============================================================================
126 : !=============================================================================
127 1536 : subroutine aero_model_init( pbuf2d )
128 :
129 : use mo_chem_utls, only: get_inv_ndx, get_spc_ndx
130 : use cam_history, only: addfld, add_default, horiz_only
131 : use phys_control, only: phys_getopts
132 : use mo_aerosols, only: aerosols_inti
133 : use mo_setsoa, only: soa_inti
134 : use dust_model, only: dust_init
135 : use seasalt_model, only: seasalt_init
136 : use aer_drydep_mod, only: inidrydep
137 : use wetdep, only: wetdep_init
138 : use mo_setsox, only: has_sox
139 :
140 : ! args
141 : type(physics_buffer_desc), pointer :: pbuf2d(:,:)
142 :
143 : ! local vars
144 : character(len=12), parameter :: subrname = 'aero_model_init'
145 : integer :: m, id
146 : character(len=20) :: dummy
147 : logical :: history_aerosol ! Output MAM or SECT aerosol tendencies
148 : logical :: history_dust ! Output dust
149 :
150 : call phys_getopts( history_aerosol_out = history_aerosol,&
151 1536 : history_dust_out = history_dust )
152 :
153 1536 : call aerosols_inti()
154 1536 : call soa_inti(pbuf2d)
155 1536 : call dust_init()
156 1536 : call seasalt_init()
157 1536 : call wetdep_init()
158 :
159 1536 : fracis_idx = pbuf_get_index('FRACIS')
160 :
161 1536 : nwetdep = 0
162 1536 : ndrydep = 0
163 :
164 18432 : count_species: do m = 1,pcnst
165 16896 : if ( len_trim(wetdep_list(m)) /= 0 ) then
166 0 : nwetdep = nwetdep+1
167 : endif
168 18432 : if ( len_trim(drydep_list(m)) /= 0 ) then
169 0 : ndrydep = ndrydep+1
170 : endif
171 : enddo count_species
172 :
173 1536 : if (nwetdep>0) &
174 0 : allocate(wetdep_indices(nwetdep))
175 1536 : if (ndrydep>0) &
176 0 : allocate(drydep_indices(ndrydep))
177 :
178 1536 : do m = 1,ndrydep
179 0 : call cnst_get_ind ( drydep_list(m), id, abort=.false. )
180 0 : if (id>0) then
181 0 : drydep_indices(m) = id
182 : else
183 0 : call endrun(subrname//': invalid drydep species: '//trim(drydep_list(m)) )
184 : endif
185 :
186 1536 : if (masterproc) then
187 0 : write(iulog,*) subrname//': '//drydep_list(m)//' will have drydep applied'
188 : endif
189 : enddo
190 1536 : do m = 1,nwetdep
191 0 : call cnst_get_ind ( wetdep_list(m), id, abort=.false. )
192 0 : if (id>0) then
193 0 : wetdep_indices(m) = id
194 : else
195 0 : call endrun(subrname//': invalid wetdep species: '//trim(wetdep_list(m)) )
196 : endif
197 :
198 1536 : if (masterproc) then
199 0 : write(iulog,*) subrname//': '//wetdep_list(m)//' will have wet removal'
200 : endif
201 : enddo
202 :
203 : ! set flags for drydep tendencies
204 1536 : drydep_lq(:) = .false.
205 1536 : do m=1,ndrydep
206 0 : id = drydep_indices(m)
207 1536 : drydep_lq(id) = .true.
208 : enddo
209 :
210 : ! set flags for wetdep tendencies
211 1536 : wetdep_lq(:) = .false.
212 1536 : do m=1,nwetdep
213 0 : id = wetdep_indices(m)
214 1536 : wetdep_lq(id) = .true.
215 : enddo
216 :
217 1536 : do m = 1,ndrydep
218 :
219 0 : dummy = trim(drydep_list(m)) // 'TB'
220 0 : call addfld (dummy,horiz_only, 'A','kg/m2/s',trim(drydep_list(m))//' turbulent dry deposition flux')
221 0 : if ( history_aerosol ) then
222 0 : call add_default (dummy, 1, ' ')
223 : endif
224 0 : dummy = trim(drydep_list(m)) // 'GV'
225 0 : call addfld (dummy,horiz_only, 'A','kg/m2/s',trim(drydep_list(m)) //' gravitational dry deposition flux')
226 0 : if ( history_aerosol ) then
227 0 : call add_default (dummy, 1, ' ')
228 : endif
229 0 : dummy = trim(drydep_list(m)) // 'DD'
230 0 : call addfld (dummy,horiz_only, 'A','kg/m2/s',trim(drydep_list(m)) //' dry deposition flux at bottom (grav + turb)')
231 0 : if ( history_aerosol ) then
232 0 : call add_default (dummy, 1, ' ')
233 : endif
234 0 : dummy = trim(drydep_list(m)) // 'DT'
235 0 : call addfld (dummy,(/ 'lev' /), 'A','kg/kg/s',trim(drydep_list(m))//' dry deposition')
236 0 : if ( history_aerosol ) then
237 0 : call add_default (dummy, 1, ' ')
238 : endif
239 0 : dummy = trim(drydep_list(m)) // 'DV'
240 0 : call addfld (dummy,(/ 'lev' /), 'A','m/s',trim(drydep_list(m))//' deposition velocity')
241 1536 : if ( history_aerosol ) then
242 0 : call add_default (dummy, 1, ' ')
243 : endif
244 :
245 : enddo
246 :
247 1536 : if (ndrydep>0) then
248 :
249 0 : call inidrydep(rair, gravit)
250 :
251 0 : dummy = 'RAM1'
252 0 : call addfld (dummy,horiz_only, 'A','frac','RAM1')
253 0 : if ( history_aerosol ) then
254 0 : call add_default (dummy, 1, ' ')
255 : endif
256 0 : dummy = 'airFV'
257 0 : call addfld (dummy,horiz_only, 'A','frac','FV')
258 0 : if ( history_aerosol ) then
259 0 : call add_default (dummy, 1, ' ')
260 : endif
261 :
262 0 : if (sslt_active) then
263 0 : dummy = 'SSTSFDRY'
264 0 : call addfld (dummy,horiz_only, 'A','kg/m2/s','Sea salt deposition flux at surface')
265 0 : if ( history_aerosol ) then
266 0 : call add_default (dummy, 1, ' ')
267 : endif
268 : endif
269 0 : if (dust_active) then
270 0 : dummy = 'DSTSFDRY'
271 0 : call addfld (dummy,horiz_only, 'A','kg/m2/s','Dust deposition flux at surface')
272 0 : if ( history_aerosol ) then
273 0 : call add_default (dummy, 1, ' ')
274 : endif
275 : endif
276 :
277 : endif
278 :
279 1536 : do m = 1,nwetdep
280 :
281 0 : call addfld (trim(wetdep_list(m))//'SFWET', horiz_only, 'A','kg/m2/s', &
282 0 : 'Wet deposition flux at surface')
283 : call addfld (trim(wetdep_list(m))//'SFSIC', horiz_only, 'A','kg/m2/s', &
284 0 : 'Wet deposition flux (incloud, convective) at surface')
285 : call addfld (trim(wetdep_list(m))//'SFSIS', horiz_only, 'A','kg/m2/s', &
286 0 : 'Wet deposition flux (incloud, stratiform) at surface')
287 : call addfld (trim(wetdep_list(m))//'SFSBC', horiz_only, 'A','kg/m2/s', &
288 0 : 'Wet deposition flux (belowcloud, convective) at surface')
289 : call addfld (trim(wetdep_list(m))//'SFSBS', horiz_only, 'A','kg/m2/s', &
290 0 : 'Wet deposition flux (belowcloud, stratiform) at surface')
291 : call addfld (trim(wetdep_list(m))//'WET', (/ 'lev' /), 'A','kg/kg/s', &
292 0 : 'wet deposition tendency')
293 : call addfld (trim(wetdep_list(m))//'SIC', (/ 'lev' /), 'A','kg/kg/s', &
294 0 : trim(wetdep_list(m))//' ic wet deposition')
295 : call addfld (trim(wetdep_list(m))//'SIS', (/ 'lev' /), 'A','kg/kg/s', &
296 0 : trim(wetdep_list(m))//' is wet deposition')
297 : call addfld (trim(wetdep_list(m))//'SBC', (/ 'lev' /), 'A','kg/kg/s', &
298 0 : trim(wetdep_list(m))//' bc wet deposition')
299 : call addfld (trim(wetdep_list(m))//'SBS', (/ 'lev' /), 'A','kg/kg/s', &
300 1536 : trim(wetdep_list(m))//' bs wet deposition')
301 : enddo
302 :
303 1536 : if (nwetdep>0) then
304 0 : if (sslt_active) then
305 0 : dummy = 'SSTSFWET'
306 0 : call addfld (dummy,horiz_only, 'A','kg/m2/s','Sea salt wet deposition flux at surface')
307 0 : if ( history_aerosol ) then
308 0 : call add_default (dummy, 1, ' ')
309 : endif
310 : endif
311 0 : if (dust_active) then
312 0 : dummy = 'DSTSFWET'
313 0 : call addfld (dummy,horiz_only, 'A','kg/m2/s','Dust wet deposition flux at surface')
314 0 : if ( history_aerosol ) then
315 0 : call add_default (dummy, 1, ' ')
316 : endif
317 : endif
318 : endif
319 :
320 1536 : if (dust_active) then
321 : ! emissions diagnostics ....
322 :
323 0 : do m = 1, dust_nbin
324 0 : dummy = trim(dust_names(m)) // 'SF'
325 0 : call addfld (dummy,horiz_only, 'A','kg/m2/s',trim(dust_names(m))//' dust surface emission')
326 0 : if (history_aerosol) then
327 0 : call add_default (dummy, 1, ' ')
328 : endif
329 : enddo
330 :
331 0 : dummy = 'DSTSFMBL'
332 0 : call addfld (dummy,horiz_only, 'A','kg/m2/s','Mobilization flux at surface')
333 0 : if (history_aerosol .or. history_dust) then
334 0 : call add_default (dummy, 1, ' ')
335 : endif
336 :
337 0 : dummy = 'LND_MBL'
338 0 : call addfld (dummy,horiz_only, 'A','frac','Soil erodibility factor')
339 0 : if (history_aerosol) then
340 0 : call add_default (dummy, 1, ' ')
341 : endif
342 :
343 : endif
344 :
345 1536 : if (sslt_active) then
346 :
347 0 : dummy = 'SSTSFMBL'
348 0 : call addfld (dummy,horiz_only, 'A','kg/m2/s','Mobilization flux at surface')
349 0 : if (history_aerosol) then
350 0 : call add_default (dummy, 1, ' ')
351 : endif
352 :
353 0 : do m = 1, seasalt_nbin
354 0 : dummy = trim(seasalt_names(m)) // 'SF'
355 0 : call addfld (dummy,horiz_only, 'A','kg/m2/s',trim(seasalt_names(m))//' seasalt surface emission')
356 0 : if (history_aerosol) then
357 0 : call add_default (dummy, 1, ' ')
358 : endif
359 : enddo
360 :
361 : endif
362 :
363 1536 : if( has_sox ) then
364 3072 : call addfld( 'XPH_LWC',(/ 'lev' /), 'A','kg/kg', 'pH value multiplied by lwc')
365 :
366 1536 : if ( history_aerosol ) then
367 0 : call add_default ('XPH_LWC', 1, ' ')
368 : endif
369 : endif
370 :
371 1536 : so4_ndx = get_spc_ndx( 'SO4' )
372 1536 : soa_ndx = get_spc_ndx( 'SOA' )
373 1536 : soai_ndx = get_spc_ndx( 'SOAI' )
374 1536 : soam_ndx = get_spc_ndx( 'SOAM' )
375 1536 : soab_ndx = get_spc_ndx( 'SOAB' )
376 1536 : soat_ndx = get_spc_ndx( 'SOAT' )
377 1536 : soax_ndx = get_spc_ndx( 'SOAX' )
378 1536 : cb2_ndx = get_spc_ndx( 'CB2' )
379 1536 : oc2_ndx = get_spc_ndx( 'OC2' )
380 1536 : nit_ndx = get_spc_ndx( 'NH4NO3' )
381 :
382 1536 : end subroutine aero_model_init
383 :
384 : !=============================================================================
385 : !=============================================================================
386 19359288 : subroutine aero_model_drydep ( state, pbuf, obklen, ustar, cam_in, dt, cam_out, ptend )
387 :
388 1536 : use dust_sediment_mod, only: dust_sediment_tend
389 : use aer_drydep_mod, only: d3ddflux, calcram
390 : use dust_model, only: dust_depvel, dust_nbin, dust_names
391 : use seasalt_model, only: sslt_depvel=>seasalt_depvel, sslt_nbin=>seasalt_nbin, sslt_names=>seasalt_names
392 :
393 : ! args
394 : type(physics_state), intent(in) :: state ! Physics state variables
395 : real(r8), intent(in) :: obklen(:)
396 : real(r8), intent(in) :: ustar(:) ! sfc fric vel
397 : type(cam_in_t), target, intent(in) :: cam_in ! import state
398 : real(r8), intent(in) :: dt ! time step
399 : type(cam_out_t), intent(inout) :: cam_out ! export state
400 : type(physics_ptend), intent(out) :: ptend ! indivdual parameterization tendencies
401 : type(physics_buffer_desc), pointer :: pbuf(:)
402 :
403 : ! local vars
404 1489176 : real(r8), pointer :: landfrac(:) ! land fraction
405 1489176 : real(r8), pointer :: icefrac(:) ! ice fraction
406 1489176 : real(r8), pointer :: ocnfrac(:) ! ocean fraction
407 1489176 : real(r8), pointer :: fvin(:) !
408 1489176 : real(r8), pointer :: ram1in(:) ! for dry dep velocities from land model for progseasalts
409 :
410 : real(r8) :: fv(pcols) ! for dry dep velocities, from land modified over ocean & ice
411 : real(r8) :: ram1(pcols) ! for dry dep velocities, from land modified over ocean & ice
412 :
413 : ! local decarations
414 :
415 : integer, parameter :: naero = sslt_nbin+dust_nbin
416 : integer, parameter :: begslt = 1
417 : integer, parameter :: endslt = sslt_nbin
418 : integer, parameter :: begdst = sslt_nbin+1
419 : integer, parameter :: enddst = sslt_nbin+dust_nbin
420 :
421 : integer :: ncol, lchnk
422 :
423 : character(len=6) :: aeronames(naero) ! = (/ sslt_names, dust_names /)
424 :
425 : real(r8) :: vlc_trb(pcols,naero) !Turbulent deposn velocity (m/s)
426 : real(r8) :: vlc_grv(pcols,pver,naero) !grav deposn velocity (m/s)
427 : real(r8) :: vlc_dry(pcols,pver,naero) !dry deposn velocity (m/s)
428 :
429 : real(r8) :: dep_trb(pcols) !kg/m2/s
430 : real(r8) :: dep_grv(pcols) !kg/m2/s (total of grav and trb)
431 :
432 : real(r8) :: tsflx_dst(pcols)
433 : real(r8) :: tsflx_slt(pcols)
434 : real(r8) :: pvaeros(pcols,pverp) ! sedimentation velocity in Pa
435 : real(r8) :: sflx(pcols)
436 :
437 : real(r8) :: tvs(pcols,pver)
438 : real(r8) :: rho(pcols,pver) ! air density in kg/m3
439 :
440 : integer :: m,mm, i, im
441 :
442 1489176 : if (ndrydep<1) return
443 :
444 0 : landfrac => cam_in%landfrac(:)
445 0 : icefrac => cam_in%icefrac(:)
446 0 : ocnfrac => cam_in%ocnfrac(:)
447 0 : fvin => cam_in%fv(:)
448 0 : ram1in => cam_in%ram1(:)
449 :
450 0 : lchnk = state%lchnk
451 0 : ncol = state%ncol
452 :
453 : ! calc ram and fv over ocean and sea ice ...
454 : call calcram( ncol,landfrac,icefrac,ocnfrac,obklen,&
455 : ustar,ram1in,ram1,state%t(:,pver),state%pmid(:,pver),&
456 0 : state%pdel(:,pver),fvin,fv)
457 :
458 0 : call outfld( 'airFV', fv(:), pcols, lchnk )
459 0 : call outfld( 'RAM1', ram1(:), pcols, lchnk )
460 :
461 : ! note that tendencies are not only in sfc layer (because of sedimentation)
462 : ! and that ptend is updated within each subroutine for different species
463 :
464 0 : call physics_ptend_init(ptend, state%psetcols, 'aero_model_drydep', lq=drydep_lq)
465 :
466 0 : aeronames(:sslt_nbin) = sslt_names(:)
467 0 : aeronames(sslt_nbin+1:) = dust_names(:)
468 :
469 0 : lchnk = state%lchnk
470 0 : ncol = state%ncol
471 :
472 0 : tvs(:ncol,:) = state%t(:ncol,:)
473 0 : rho(:ncol,:) = state%pmid(:ncol,:)/(rair*state%t(:ncol,:))
474 :
475 : ! compute dep velocities for sea salt and dust...
476 0 : if (sslt_active) then
477 : call sslt_depvel( state%t(:,:), state%pmid(:,:), state%q(:,:,1), ram1, fv, ncol, lchnk, &
478 0 : vlc_dry(:,:,begslt:endslt), vlc_trb(:,begslt:endslt), vlc_grv(:,:,begslt:endslt))
479 : endif
480 0 : if (dust_active) then
481 : call dust_depvel( state%t(:,:), state%pmid(:,:), ram1, fv, ncol, &
482 0 : vlc_dry(:,:,begdst:enddst), vlc_trb(:,begdst:enddst), vlc_grv(:,:,begdst:enddst) )
483 : endif
484 :
485 0 : tsflx_dst(:)=0._r8
486 0 : tsflx_slt(:)=0._r8
487 :
488 : ! do drydep for each of the bins of dust and seasalt
489 0 : do m=1,ndrydep
490 :
491 0 : mm = drydep_indices(m)
492 0 : findindex: do im = 1,naero
493 0 : if (trim(cnst_name(mm))==trim(aeronames(im))) exit findindex
494 : enddo findindex
495 :
496 0 : pvaeros(:ncol,1)=0._r8
497 0 : pvaeros(:ncol,2:pverp) = vlc_dry(:ncol,:,im)
498 :
499 0 : call outfld( trim(cnst_name(mm))//'DV', pvaeros(:,2:pverp), pcols, lchnk )
500 :
501 : if(.true.) then ! use phil's method
502 : ! convert from meters/sec to pascals/sec
503 : ! pvaeros(:,1) is assumed zero, use density from layer above in conversion
504 0 : pvaeros(:ncol,2:pverp) = pvaeros(:ncol,2:pverp) * rho(:ncol,:)*gravit
505 :
506 : ! calculate the tendencies and sfc fluxes from the above velocities
507 : call dust_sediment_tend( &
508 : ncol, dt, state%pint(:,:), state%pmid, state%pdel, state%t , &
509 0 : state%q(:,:,mm) , pvaeros , ptend%q(:,:,mm), sflx )
510 : else !use charlie's method
511 : call d3ddflux(ncol, vlc_dry(:,:,im), state%q(:,:,mm),state%pmid,state%pdel, tvs,sflx,ptend%q(:,:,mm),dt)
512 : endif
513 : ! apportion dry deposition into turb and gravitational settling for tapes
514 0 : do i=1,ncol
515 0 : dep_trb(i)=sflx(i)*vlc_trb(i,im)/vlc_dry(i,pver,im)
516 0 : dep_grv(i)=sflx(i)*vlc_grv(i,pver,im)/vlc_dry(i,pver,im)
517 : enddo
518 :
519 0 : if ( any( sslt_names(:)==trim(cnst_name(mm)) ) ) &
520 0 : tsflx_slt(:ncol)=tsflx_slt(:ncol)+sflx(:ncol)
521 0 : if ( any( dust_names(:)==trim(cnst_name(mm)) ) ) &
522 0 : tsflx_dst(:ncol)=tsflx_dst(:ncol)+sflx(:ncol)
523 :
524 : ! if the user has specified prescribed aerosol dep fluxes then
525 : ! do not set cam_out dep fluxes according to the prognostic aerosols
526 0 : if (.not. aerodep_flx_prescribed()) then
527 : ! set deposition in export state
528 0 : if (im==begdst) then
529 0 : cam_out%dstdry1(:ncol) = max(sflx(:ncol), 0._r8)
530 0 : elseif(im==begdst+1) then
531 0 : cam_out%dstdry2(:ncol) = max(sflx(:ncol), 0._r8)
532 0 : elseif(im==begdst+2) then
533 0 : cam_out%dstdry3(:ncol) = max(sflx(:ncol), 0._r8)
534 0 : elseif(im==begdst+3) then
535 0 : cam_out%dstdry4(:ncol) = max(sflx(:ncol), 0._r8)
536 : endif
537 : endif
538 :
539 0 : call outfld( trim(cnst_name(mm))//'DD', sflx, pcols, lchnk)
540 0 : call outfld( trim(cnst_name(mm))//'TB', dep_trb, pcols, lchnk )
541 0 : call outfld( trim(cnst_name(mm))//'GV', dep_grv, pcols, lchnk )
542 0 : call outfld( trim(cnst_name(mm))//'DT', ptend%q(:,:,mm), pcols, lchnk)
543 :
544 : end do
545 :
546 : ! output the total dry deposition
547 0 : if (sslt_active) then
548 0 : call outfld( 'SSTSFDRY', tsflx_slt, pcols, lchnk)
549 : endif
550 0 : if (dust_active) then
551 0 : call outfld( 'DSTSFDRY', tsflx_dst, pcols, lchnk)
552 : endif
553 :
554 1489176 : endsubroutine aero_model_drydep
555 :
556 : !=============================================================================
557 : !=============================================================================
558 17870112 : subroutine aero_model_wetdep( state, dt, dlf, cam_out, ptend, pbuf)
559 :
560 : use wetdep, only : wetdepa_v1, wetdep_inputs_set, wetdep_inputs_t
561 : use dust_model, only : dust_names
562 : use seasalt_model, only : sslt_names=>seasalt_names
563 :
564 : ! args
565 :
566 : type(physics_state), intent(in) :: state ! Physics state variables
567 : real(r8), intent(in) :: dt ! time step
568 : real(r8), intent(in) :: dlf(:,:) ! shallow+deep convective detrainment [kg/kg/s]
569 : type(cam_out_t), intent(inout) :: cam_out ! export state
570 : type(physics_ptend), intent(out) :: ptend ! indivdual parameterization tendencies
571 : type(physics_buffer_desc), pointer :: pbuf(:)
572 :
573 : ! local vars
574 :
575 : integer :: ncol ! number of atmospheric columns
576 : integer :: lchnk ! chunk identifier
577 : integer :: m,mm, i,k
578 :
579 : real(r8) :: sflx_tot_dst(pcols)
580 : real(r8) :: sflx_tot_slt(pcols)
581 :
582 : real(r8) :: iscavt(pcols, pver)
583 : real(r8) :: scavt(pcols, pver)
584 : real(r8) :: scavcoef(pcols,pver) ! Dana and Hales coefficient (/mm) (0.1)
585 : real(r8) :: sflx(pcols) ! deposition flux
586 :
587 : real(r8) :: icscavt(pcols, pver)
588 : real(r8) :: isscavt(pcols, pver)
589 : real(r8) :: bcscavt(pcols, pver)
590 : real(r8) :: bsscavt(pcols, pver)
591 :
592 : real(r8) :: sol_factb, sol_facti
593 :
594 : real(r8) :: rainmr(pcols,pver) ! mixing ratio of rain within cloud volume
595 : real(r8) :: cldv(pcols,pver) ! cloudy volume undergoing scavenging
596 : real(r8) :: cldvcu(pcols,pver) ! Convective precipitation area at the top interface of current layer
597 : real(r8) :: cldvst(pcols,pver) ! Stratiform precipitation area at the top interface of current layer
598 :
599 1489176 : real(r8), pointer :: fracis(:,:,:) ! fraction of transported species that are insoluble
600 :
601 : type(wetdep_inputs_t) :: dep_inputs ! obj that contains inputs to wetdepa routine
602 :
603 1489176 : if (nwetdep<1) return
604 :
605 0 : call pbuf_get_field(pbuf, fracis_idx, fracis, start=(/1,1,1/), kount=(/pcols, pver, pcnst/) )
606 :
607 0 : call physics_ptend_init(ptend, state%psetcols, 'aero_model_wetdep', lq=wetdep_lq)
608 :
609 0 : call wetdep_inputs_set( state, pbuf, dep_inputs )
610 :
611 0 : lchnk = state%lchnk
612 0 : ncol = state%ncol
613 :
614 0 : sflx_tot_dst(:) = 0._r8
615 0 : sflx_tot_slt(:) = 0._r8
616 :
617 0 : do m = 1, nwetdep
618 :
619 0 : mm = wetdep_indices(m)
620 :
621 0 : sol_factb = aer_sol_factb(m)
622 0 : sol_facti = aer_sol_facti(m)
623 :
624 0 : scavcoef(:ncol,:) = aer_scav_coef(m)
625 :
626 : call wetdepa_v1( state%t, state%pmid, state%q(:,:,1), state%pdel, &
627 : dep_inputs%cldt, dep_inputs%cldcu, dep_inputs%cmfdqr, &
628 : dep_inputs%conicw, dep_inputs%prain, dep_inputs%qme, &
629 : dep_inputs%evapr, dep_inputs%totcond, state%q(:,:,mm), dt, &
630 : scavt, iscavt, dep_inputs%cldv, &
631 : fracis(:,:,mm), sol_factb, ncol, &
632 : scavcoef, &
633 : sol_facti_in=sol_facti, &
634 0 : icscavt=icscavt, isscavt=isscavt, bcscavt=bcscavt, bsscavt=bsscavt )
635 :
636 0 : ptend%q(:ncol,:,mm)=scavt(:ncol,:)
637 :
638 0 : call outfld( trim(cnst_name(mm))//'WET', ptend%q(:,:,mm), pcols, lchnk)
639 0 : call outfld( trim(cnst_name(mm))//'SIC', icscavt , pcols, lchnk)
640 0 : call outfld( trim(cnst_name(mm))//'SIS', isscavt, pcols, lchnk)
641 0 : call outfld( trim(cnst_name(mm))//'SBC', bcscavt, pcols, lchnk)
642 0 : call outfld( trim(cnst_name(mm))//'SBS', bsscavt, pcols, lchnk)
643 :
644 0 : sflx(:)=0._r8
645 :
646 0 : do k=1,pver
647 0 : do i=1,ncol
648 0 : sflx(i)=sflx(i)+ptend%q(i,k,mm)*state%pdel(i,k)/gravit
649 : enddo
650 : enddo
651 0 : call outfld( trim(cnst_name(mm))//'SFWET', sflx, pcols, lchnk)
652 :
653 0 : if ( any( sslt_names(:)==trim(cnst_name(mm)) ) ) &
654 0 : sflx_tot_slt(:ncol) = sflx_tot_slt(:ncol) + sflx(:ncol)
655 0 : if ( any( dust_names(:)==trim(cnst_name(mm)) ) ) &
656 0 : sflx_tot_dst(:ncol) = sflx_tot_dst(:ncol) + sflx(:ncol)
657 :
658 : ! if the user has specified prescribed aerosol dep fluxes then
659 : ! do not set cam_out dep fluxes according to the prognostic aerosols
660 0 : if (.not.aerodep_flx_prescribed()) then
661 : ! export deposition fluxes to coupler ??? why "-" sign ???
662 0 : if (trim(cnst_name(mm))=='CB2') then
663 0 : cam_out%bcphiwet(:) = max(-sflx(:), 0._r8)
664 0 : elseif (trim(cnst_name(mm))=='OC2') then
665 0 : cam_out%ocphiwet(:) = max(-sflx(:), 0._r8)
666 0 : elseif (trim(cnst_name(mm))==trim(dust_names(1))) then
667 0 : cam_out%dstwet1(:) = max(-sflx(:), 0._r8)
668 0 : elseif (trim(cnst_name(mm))==trim(dust_names(2))) then
669 0 : cam_out%dstwet2(:) = max(-sflx(:), 0._r8)
670 0 : elseif (trim(cnst_name(mm))==trim(dust_names(3))) then
671 0 : cam_out%dstwet3(:) = max(-sflx(:), 0._r8)
672 0 : elseif (trim(cnst_name(mm))==trim(dust_names(4))) then
673 0 : cam_out%dstwet4(:) = max(-sflx(:), 0._r8)
674 : endif
675 : endif
676 :
677 : enddo
678 :
679 0 : if (sslt_active) then
680 0 : call outfld( 'SSTSFWET', sflx_tot_slt, pcols, lchnk)
681 : endif
682 0 : if (dust_active) then
683 0 : call outfld( 'DSTSFWET', sflx_tot_dst, pcols, lchnk)
684 : endif
685 :
686 1489176 : endsubroutine aero_model_wetdep
687 :
688 : !-------------------------------------------------------------------------
689 : ! provides aerosol surface area info for sectional aerosols
690 : ! called from mo_usrrxt
691 : !-------------------------------------------------------------------------
692 0 : subroutine aero_model_surfarea( &
693 0 : mmr, radmean, relhum, pmid, temp, strato_sad, sulfate, m, ltrop, &
694 0 : dlat, het1_ndx, pbuf, ncol, sfc, dm_aer, sad_total, reff_trop )
695 :
696 : use mo_constants, only : pi, avo => avogadro
697 :
698 : ! dummy args
699 : real(r8), intent(in) :: pmid(:,:)
700 : real(r8), intent(in) :: temp(:,:)
701 : real(r8), intent(in) :: mmr(:,:,:)
702 : real(r8), intent(in) :: radmean ! mean radii in cm
703 : real(r8), intent(in) :: strato_sad(:,:)
704 : integer, intent(in) :: ncol
705 : integer, intent(in) :: ltrop(:)
706 : real(r8), intent(in) :: dlat(:) ! degrees latitude
707 : integer, intent(in) :: het1_ndx
708 : real(r8), intent(in) :: relhum(:,:)
709 : real(r8), intent(in) :: m(:,:) ! total atm density (/cm^3)
710 : real(r8), intent(in) :: sulfate(:,:)
711 : type(physics_buffer_desc), pointer :: pbuf(:)
712 :
713 : real(r8), intent(inout) :: sfc(:,:,:)
714 : real(r8), intent(inout) :: dm_aer(:,:,:)
715 : real(r8), intent(inout) :: sad_total(:,:)
716 : real(r8), intent(out) :: reff_trop(:,:)
717 :
718 : ! local vars
719 :
720 : integer :: i,k
721 : real(r8) :: rho_air
722 : real(r8) :: v, n, n_exp, r_rd, r_sd
723 : real(r8) :: dm_sulf, dm_sulf_wet, log_sd_sulf, sfc_sulf, sfc_nit
724 : real(r8) :: dm_orgc, dm_orgc_wet, log_sd_orgc, sfc_oc, sfc_soa
725 : real(r8) :: sfc_soai, sfc_soam, sfc_soab, sfc_soat, sfc_soax
726 : real(r8) :: dm_bc, dm_bc_wet, log_sd_bc, sfc_bc
727 : real(r8) :: rxt_sulf, rxt_nit, rxt_oc, rxt_soa
728 : real(r8) :: c_n2o5, c_ho2, c_no2, c_no3
729 : real(r8) :: s_exp
730 :
731 : !-----------------------------------------------------------------
732 : ! ... parameters for log-normal distribution by number
733 : ! references:
734 : ! Chin et al., JAS, 59, 461, 2003
735 : ! Liao et al., JGR, 108(D1), 4001, 2003
736 : ! Martin et al., JGR, 108(D3), 4097, 2003
737 : !-----------------------------------------------------------------
738 : real(r8), parameter :: rm_sulf = 6.95e-6_r8 ! mean radius of sulfate particles (cm) (Chin)
739 : real(r8), parameter :: sd_sulf = 2.03_r8 ! standard deviation of radius for sulfate (Chin)
740 : real(r8), parameter :: rho_sulf = 1.7e3_r8 ! density of sulfate aerosols (kg/m3) (Chin)
741 :
742 : real(r8), parameter :: rm_orgc = 2.12e-6_r8 ! mean radius of organic carbon particles (cm) (Chin)
743 : real(r8), parameter :: sd_orgc = 2.20_r8 ! standard deviation of radius for OC (Chin)
744 : real(r8), parameter :: rho_orgc = 1.8e3_r8 ! density of OC aerosols (kg/m3) (Chin)
745 :
746 : real(r8), parameter :: rm_bc = 1.18e-6_r8 ! mean radius of soot/BC particles (cm) (Chin)
747 : real(r8), parameter :: sd_bc = 2.00_r8 ! standard deviation of radius for BC (Chin)
748 : real(r8), parameter :: rho_bc = 1.0e3_r8 ! density of BC aerosols (kg/m3) (Chin)
749 :
750 : real(r8), parameter :: mw_so4 = 98.e-3_r8 ! so4 molecular wt (kg/mole)
751 :
752 : integer :: irh, rh_l, rh_u
753 : real(r8) :: factor, rfac_sulf, rfac_oc, rfac_bc, rfac_ss
754 : logical :: zero_aerosols
755 :
756 : !-----------------------------------------------------------------
757 : ! ... table for hygroscopic growth effect on radius (Chin et al)
758 : ! (no growth effect for mineral dust)
759 : !-----------------------------------------------------------------
760 : real(r8), dimension(7) :: table_rh, table_rfac_sulf, table_rfac_bc, table_rfac_oc, table_rfac_ss
761 :
762 : data table_rh(1:7) / 0.0_r8, 0.5_r8, 0.7_r8, 0.8_r8, 0.9_r8, 0.95_r8, 0.99_r8/
763 : data table_rfac_sulf(1:7) / 1.0_r8, 1.4_r8, 1.5_r8, 1.6_r8, 1.8_r8, 1.9_r8, 2.2_r8/
764 : data table_rfac_oc(1:7) / 1.0_r8, 1.2_r8, 1.4_r8, 1.5_r8, 1.6_r8, 1.8_r8, 2.2_r8/
765 : data table_rfac_bc(1:7) / 1.0_r8, 1.0_r8, 1.0_r8, 1.2_r8, 1.4_r8, 1.5_r8, 1.9_r8/
766 : data table_rfac_ss(1:7) / 1.0_r8, 1.6_r8, 1.8_r8, 2.0_r8, 2.4_r8, 2.9_r8, 4.8_r8/
767 :
768 : !-----------------------------------------------------------------
769 : ! ... exponent for calculating number density
770 : !-----------------------------------------------------------------
771 0 : n_exp = exp( -4.5_r8*log(sd_sulf)*log(sd_sulf) )
772 :
773 0 : dm_sulf = 2._r8 * rm_sulf
774 0 : dm_orgc = 2._r8 * rm_orgc
775 0 : dm_bc = 2._r8 * rm_bc
776 :
777 0 : log_sd_sulf = log(sd_sulf)
778 0 : log_sd_orgc = log(sd_orgc)
779 0 : log_sd_bc = log(sd_bc)
780 :
781 0 : reff_trop(:,:) = 0._r8
782 :
783 0 : ver_loop: do k = 1,pver
784 0 : col_loop: do i = 1,ncol
785 : !-------------------------------------------------------------------------
786 : ! ... air density (kg/m3)
787 : !-------------------------------------------------------------------------
788 0 : rho_air = pmid(i,k)/(temp(i,k)*287.04_r8)
789 : !-------------------------------------------------------------------------
790 : ! ... aerosol growth interpolated from M.Chin's table
791 : !-------------------------------------------------------------------------
792 0 : if (relhum(i,k) >= table_rh(7)) then
793 0 : rfac_sulf = table_rfac_sulf(7)
794 0 : rfac_oc = table_rfac_oc(7)
795 0 : rfac_bc = table_rfac_bc(7)
796 : else
797 0 : do irh = 2,7
798 0 : if (relhum(i,k) <= table_rh(irh)) then
799 : exit
800 : end if
801 : end do
802 0 : rh_l = irh-1
803 0 : rh_u = irh
804 :
805 0 : factor = (relhum(i,k) - table_rh(rh_l))/(table_rh(rh_u) - table_rh(rh_l))
806 :
807 0 : rfac_sulf = table_rfac_sulf(rh_l) + factor*(table_rfac_sulf(rh_u) - table_rfac_sulf(rh_l))
808 0 : rfac_oc = table_rfac_oc(rh_u) + factor*(table_rfac_oc(rh_u) - table_rfac_oc(rh_l))
809 0 : rfac_bc = table_rfac_bc(rh_u) + factor*(table_rfac_bc(rh_u) - table_rfac_bc(rh_l))
810 : end if
811 :
812 0 : dm_sulf_wet = dm_sulf * rfac_sulf
813 0 : dm_orgc_wet = dm_orgc * rfac_oc
814 0 : dm_bc_wet = dm_bc * rfac_bc
815 :
816 0 : dm_bc_wet = min(dm_bc_wet ,50.e-6_r8) ! maximum size is 0.5 micron (Chin)
817 0 : dm_orgc_wet = min(dm_orgc_wet,50.e-6_r8) ! maximum size is 0.5 micron (Chin)
818 :
819 :
820 : !-------------------------------------------------------------------------
821 : ! ... sulfate aerosols
822 : !-------------------------------------------------------------------------
823 0 : zero_aerosols = k < ltrop(i)
824 0 : if ( abs( dlat(i) ) > 50._r8 ) then
825 0 : zero_aerosols = pmid(i,k) < 30000._r8
826 : endif
827 : !-------------------------------------------------------------------------
828 : ! ... use ubvals climatology for stratospheric sulfate surface area density
829 : !-------------------------------------------------------------------------
830 0 : if( zero_aerosols ) then
831 0 : sfc_sulf = strato_sad(i,k)
832 0 : if ( het1_ndx > 0 ) then
833 0 : sfc_sulf = 0._r8 ! reaction already taken into account in mo_strato_rates.F90
834 : end if
835 : sfc_nit = 0._r8
836 : sfc_soa = 0._r8
837 : sfc_oc = 0._r8
838 : sfc_bc = 0._r8
839 : else
840 :
841 0 : if( so4_ndx > 0 ) then
842 : !-------------------------------------------------------------------------
843 : ! convert mass mixing ratio of aerosol to cm3/cm3 (cm^3_aerosol/cm^3_air)
844 : ! v=volume density (m^3/m^3)
845 : ! rho_aer=density of aerosol (kg/m^3)
846 : ! v=m*rho_air/rho_aer [kg/kg * (kg/m3)_air/(kg/m3)_aer]
847 : !-------------------------------------------------------------------------
848 0 : v = mmr(i,k,so4_ndx) * rho_air/rho_sulf
849 : !-------------------------------------------------------------------------
850 : ! calculate the number density of aerosol (aerosols/cm3)
851 : ! assuming a lognormal distribution
852 : ! n = (aerosols/cm3)
853 : ! dm = geometric mean diameter
854 : !
855 : ! because only the dry mass of the aerosols is known, we
856 : ! use the mean dry radius
857 : !-------------------------------------------------------------------------
858 0 : n = v * (6._r8/pi)*(1._r8/(dm_sulf**3._r8))*n_exp
859 : !-------------------------------------------------------------------------
860 : ! find surface area of aerosols using dm_wet, log_sd
861 : ! (increase of sd due to RH is negligible)
862 : ! and number density calculated above as distribution
863 : ! parameters
864 : ! sfc = surface area of wet aerosols (cm^2/cm^3)
865 : !-------------------------------------------------------------------------
866 0 : s_exp = exp(2._r8*log_sd_sulf*log_sd_sulf)
867 0 : sfc_sulf = n * pi * (dm_sulf_wet**2._r8) * s_exp
868 :
869 : else
870 : !-------------------------------------------------------------------------
871 : ! if so4 not simulated, use off-line sulfate and calculate as above
872 : ! convert sulfate vmr to volume density of aerosol (cm^3_aerosol/cm^3_air)
873 : !-------------------------------------------------------------------------
874 0 : v = sulfate(i,k) * m(i,k) * mw_so4 / (avo * rho_sulf) *1.e6_r8
875 0 : n = v * (6._r8/pi)*(1._r8/(dm_sulf**3._r8))*n_exp
876 0 : s_exp = exp(2._r8*log_sd_sulf*log_sd_sulf)
877 0 : sfc_sulf = n * pi * (dm_sulf_wet**2._r8) * s_exp
878 :
879 : end if
880 :
881 : !-------------------------------------------------------------------------
882 : ! ammonium nitrate (follow same procedure as sulfate, using size and density of sulfate)
883 : !-------------------------------------------------------------------------
884 0 : if( nit_ndx > 0 ) then
885 0 : v = mmr(i,k,nit_ndx) * rho_air/rho_sulf
886 0 : n = v * (6._r8/pi)*(1._r8/(dm_sulf**3._r8))*n_exp
887 0 : s_exp = exp(2._r8*log_sd_sulf*log_sd_sulf)
888 0 : sfc_nit = n * pi * (dm_sulf_wet**2._r8) * s_exp
889 : else
890 : sfc_nit = 0._r8
891 : end if
892 :
893 : !-------------------------------------------------------------------------
894 : ! hydrophylic organic carbon (follow same procedure as sulfate)
895 : !-------------------------------------------------------------------------
896 0 : if( oc2_ndx > 0 ) then
897 0 : v = mmr(i,k,oc2_ndx) * rho_air/rho_orgc
898 0 : n = v * (6._r8/pi)*(1._r8/(dm_orgc**3))*n_exp
899 0 : s_exp = exp(2._r8*log_sd_orgc*log_sd_orgc)
900 0 : sfc_oc = n * pi * (dm_orgc_wet**2._r8) * s_exp
901 : else
902 : sfc_oc = 0._r8
903 : end if
904 :
905 : !-------------------------------------------------------------------------
906 : ! secondary organic carbon (follow same procedure as sulfate)
907 : !-------------------------------------------------------------------------
908 0 : if( soa_ndx > 0 ) then
909 0 : v = mmr(i,k,soa_ndx) * rho_air/rho_orgc
910 0 : n = v * (6._r8/pi)*(1._r8/(dm_orgc**3._r8))*n_exp
911 0 : s_exp = exp(2._r8*log_sd_orgc*log_sd_orgc)
912 0 : sfc_soa = n * pi * (dm_orgc_wet**2._r8) * s_exp
913 : else
914 : sfc_soa = 0._r8
915 : end if
916 :
917 : !-------------------------------------------------------------------------
918 : ! black carbon (follow same procedure as sulfate)
919 : !-------------------------------------------------------------------------
920 0 : if( cb2_ndx > 0 ) then
921 0 : v = mmr(i,k,cb2_ndx) * rho_air/rho_bc
922 0 : n = v * (6._r8/pi)*(1._r8/(dm_bc**3._r8))*n_exp
923 0 : s_exp = exp(2._r8*log_sd_bc*log_sd_bc)
924 0 : sfc_bc = n * pi * (dm_bc_wet**2._r8) * s_exp
925 : else
926 : sfc_bc = 0._r8
927 : end if
928 0 : if( soai_ndx > 0 ) then
929 0 : v = mmr(i,k,soai_ndx) * rho_air/rho_orgc
930 0 : n = v * (6._r8/pi)*(1._r8/(dm_orgc**3._r8))*n_exp
931 0 : s_exp = exp(2._r8*log_sd_orgc*log_sd_orgc)
932 0 : sfc_soai = n * pi * (dm_orgc_wet**2._r8) * s_exp
933 : else
934 : sfc_soai = 0._r8
935 : end if
936 0 : if( soam_ndx > 0 ) then
937 0 : v = mmr(i,k,soam_ndx) * rho_air/rho_orgc
938 0 : n = v * (6._r8/pi)*(1._r8/(dm_orgc**3._r8))*n_exp
939 0 : s_exp = exp(2._r8*log_sd_orgc*log_sd_orgc)
940 0 : sfc_soam = n * pi * (dm_orgc_wet**2._r8) * s_exp
941 : else
942 : sfc_soam = 0._r8
943 : end if
944 0 : if( soab_ndx > 0 ) then
945 0 : v = mmr(i,k,soab_ndx) * rho_air/rho_orgc
946 0 : n = v * (6._r8/pi)*(1._r8/(dm_orgc**3._r8))*n_exp
947 0 : s_exp = exp(2._r8*log_sd_orgc*log_sd_orgc)
948 0 : sfc_soab = n * pi * (dm_orgc_wet**2._r8) * s_exp
949 : else
950 : sfc_soab = 0._r8
951 : end if
952 0 : if( soat_ndx > 0 ) then
953 0 : v = mmr(i,k,soat_ndx) * rho_air/rho_orgc
954 0 : n = v * (6._r8/pi)*(1._r8/(dm_orgc**3._r8))*n_exp
955 0 : s_exp = exp(2._r8*log_sd_orgc*log_sd_orgc)
956 0 : sfc_soat = n * pi * (dm_orgc_wet**2._r8) * s_exp
957 : else
958 : sfc_soat = 0._r8
959 : end if
960 0 : if( soax_ndx > 0 ) then
961 0 : v = mmr(i,k,soax_ndx) * rho_air/rho_orgc
962 0 : n = v * (6._r8/pi)*(1._r8/(dm_orgc**3._r8))*n_exp
963 0 : s_exp = exp(2._r8*log_sd_orgc*log_sd_orgc)
964 0 : sfc_soax = n * pi * (dm_orgc_wet**2._r8) * s_exp
965 : else
966 : sfc_soax = 0._r8
967 : end if
968 0 : sfc_soa = sfc_soa + sfc_soai + sfc_soam + sfc_soab + sfc_soat + sfc_soax
969 :
970 : end if
971 :
972 0 : sfc(i,k,:) = (/ sfc_sulf, sfc_nit, sfc_oc, sfc_soa, sfc_bc /)
973 0 : dm_aer(i,k,:) = (/ dm_sulf_wet,dm_sulf_wet,dm_orgc_wet,dm_orgc_wet,dm_bc_wet /)
974 :
975 : !-------------------------------------------------------------------------
976 : ! ... add up total surface area density for output
977 : !-------------------------------------------------------------------------
978 0 : sad_total(i,k) = sfc_sulf + sfc_nit + sfc_oc + sfc_soa + sfc_bc
979 :
980 : enddo col_loop
981 : enddo ver_loop
982 :
983 0 : end subroutine aero_model_surfarea
984 :
985 : !-------------------------------------------------------------------------
986 : ! stub
987 : !-------------------------------------------------------------------------
988 0 : subroutine aero_model_strat_surfarea( ncol, mmr, pmid, temp, ltrop, pbuf, strato_sad, reff_strat )
989 :
990 : ! dummy args
991 : integer, intent(in) :: ncol
992 : real(r8), intent(in) :: mmr(:,:,:)
993 : real(r8), intent(in) :: pmid(:,:)
994 : real(r8), intent(in) :: temp(:,:)
995 : integer, intent(in) :: ltrop(:) ! tropopause level indices
996 : type(physics_buffer_desc), pointer :: pbuf(:)
997 : real(r8), intent(out) :: strato_sad(:,:)
998 : real(r8), intent(out) :: reff_strat(:,:)
999 :
1000 0 : strato_sad(:,:) = 0._r8
1001 0 : reff_strat(:,:) = 0._r8
1002 :
1003 0 : end subroutine aero_model_strat_surfarea
1004 :
1005 : !=============================================================================
1006 : !=============================================================================
1007 0 : subroutine aero_model_gasaerexch( loffset, ncol, lchnk, troplev, delt, reaction_rates, &
1008 0 : tfld, pmid, pdel, mbar, relhum, &
1009 0 : zm, qh2o, cwat, cldfr, cldnum, &
1010 0 : airdens, invariants, del_h2so4_gasprod, &
1011 0 : vmr0, vmr, pbuf )
1012 :
1013 : use chem_mods, only : gas_pcnst
1014 : use mo_aerosols, only : aerosols_formation, has_aerosols
1015 : use mo_setsox, only : setsox, has_sox
1016 : use mo_setsoa, only : setsoa, has_soa
1017 :
1018 : !-----------------------------------------------------------------------
1019 : ! ... dummy arguments
1020 : !-----------------------------------------------------------------------
1021 : integer, intent(in) :: loffset ! offset applied to modal aero "pointers"
1022 : integer, intent(in) :: ncol ! number columns in chunk
1023 : integer, intent(in) :: lchnk ! chunk index
1024 : integer, intent(in) :: troplev(:)
1025 : real(r8), intent(in) :: delt ! time step size (sec)
1026 : real(r8), intent(in) :: reaction_rates(:,:,:) ! reaction rates
1027 : real(r8), intent(in) :: tfld(:,:) ! temperature (K)
1028 : real(r8), intent(in) :: pmid(:,:) ! pressure at model levels (Pa)
1029 : real(r8), intent(in) :: pdel(:,:) ! pressure thickness of levels (Pa)
1030 : real(r8), intent(in) :: mbar(:,:) ! mean wet atmospheric mass ( amu )
1031 : real(r8), intent(in) :: relhum(:,:) ! relative humidity
1032 : real(r8), intent(in) :: airdens(:,:) ! total atms density (molec/cm**3)
1033 : real(r8), intent(in) :: invariants(:,:,:)
1034 : real(r8), intent(in) :: del_h2so4_gasprod(:,:)
1035 : real(r8), intent(in) :: zm(:,:)
1036 : real(r8), intent(in) :: qh2o(:,:)
1037 : real(r8), intent(in) :: cwat(:,:) ! cloud liquid water content (kg/kg)
1038 : real(r8), intent(in) :: cldfr(:,:)
1039 : real(r8), intent(in) :: cldnum(:,:) ! droplet number concentration (#/kg)
1040 : real(r8), intent(in) :: vmr0(:,:,:) ! initial mixing ratios (before gas-phase chem changes)
1041 : real(r8), intent(inout) :: vmr(:,:,:) ! mixing ratios ( vmr )
1042 :
1043 : type(physics_buffer_desc), pointer :: pbuf(:)
1044 :
1045 : ! local vars
1046 :
1047 0 : real(r8) :: vmrcw(ncol,pver,gas_pcnst) ! cloud-borne aerosol (vmr)
1048 :
1049 0 : real(r8) :: aqso4(ncol,1) ! aqueous phase chemistry
1050 0 : real(r8) :: aqh2so4(ncol,1) ! aqueous phase chemistry
1051 0 : real(r8) :: aqso4_h2o2(ncol) ! SO4 aqueous phase chemistry due to H2O2
1052 0 : real(r8) :: aqso4_o3(ncol) ! SO4 aqueous phase chemistry due to O3
1053 0 : real(r8) :: xphlwc(ncol,pver) ! pH value multiplied by lwc
1054 :
1055 :
1056 : ! aqueous chemistry ...
1057 :
1058 0 : if( has_sox ) then
1059 : call setsox( &
1060 : ncol, &
1061 : lchnk, &
1062 : loffset, &
1063 : delt, &
1064 : pmid, &
1065 : pdel, &
1066 : tfld, &
1067 : mbar, &
1068 : cwat, &
1069 : cldfr, &
1070 : cldnum, &
1071 : airdens, &
1072 : invariants, &
1073 : vmrcw, &
1074 : vmr, &
1075 : xphlwc, &
1076 : aqso4, &
1077 : aqh2so4, &
1078 : aqso4_h2o2,&
1079 : aqso4_o3 &
1080 0 : )
1081 0 : call outfld( 'XPH_LWC',xphlwc(:ncol,:), ncol , lchnk )
1082 : endif
1083 :
1084 0 : if( has_soa ) then
1085 0 : call setsoa( ncol, lchnk, delt, reaction_rates, tfld, airdens, vmr, pbuf)
1086 : endif
1087 :
1088 0 : if( has_aerosols ) then
1089 0 : call aerosols_formation( ncol, lchnk, tfld, relhum, vmr )
1090 : endif
1091 :
1092 :
1093 0 : end subroutine aero_model_gasaerexch
1094 :
1095 : !=============================================================================
1096 : !=============================================================================
1097 0 : subroutine aero_model_emissions( state, cam_in )
1098 : use seasalt_model, only: seasalt_emis, seasalt_indices
1099 : use dust_model, only: dust_emis, dust_indices
1100 : use physics_types, only: physics_state
1101 :
1102 : ! Arguments:
1103 :
1104 : type(physics_state), intent(in) :: state ! Physics state variables
1105 : type(cam_in_t), intent(inout) :: cam_in ! import state
1106 :
1107 : ! local vars
1108 :
1109 : integer :: lchnk, ncol
1110 : integer :: m, mm
1111 : real(r8) :: soil_erod_tmp(pcols)
1112 : real(r8) :: sflx(pcols) ! accumulate over all bins for output
1113 : real(r8) :: u10cubed(pcols)
1114 : real (r8), parameter :: z0=0.0001_r8 ! m roughness length over oceans--from ocean model
1115 :
1116 0 : lchnk = state%lchnk
1117 0 : ncol = state%ncol
1118 :
1119 0 : if (dust_active) then
1120 :
1121 0 : call dust_emis( ncol, lchnk, cam_in%dstflx, cam_in%cflx, soil_erod_tmp )
1122 :
1123 : ! some dust emis diagnostics ...
1124 0 : sflx(:)=0._r8
1125 0 : do m=1,dust_nbin
1126 0 : mm = dust_indices(m)
1127 0 : sflx(:ncol)=sflx(:ncol)+cam_in%cflx(:ncol,mm)
1128 0 : call outfld(trim(dust_names(m))//'SF',cam_in%cflx(:,mm),pcols, lchnk)
1129 : enddo
1130 0 : call outfld('DSTSFMBL',sflx(:),pcols,lchnk)
1131 0 : call outfld('LND_MBL',soil_erod_tmp(:),pcols, lchnk )
1132 : endif
1133 :
1134 0 : if (sslt_active) then
1135 0 : u10cubed(:ncol)=sqrt(state%u(:ncol,pver)**2+state%v(:ncol,pver)**2)
1136 : ! move the winds to 10m high from the midpoint of the gridbox:
1137 : ! follows Tie and Seinfeld and Pandis, p.859 with math.
1138 :
1139 0 : u10cubed(:ncol)=u10cubed(:ncol)*log(10._r8/z0)/log(state%zm(:ncol,pver)/z0)
1140 :
1141 : ! we need them to the 3.41 power, according to Gong et al., 1997:
1142 0 : u10cubed(:ncol)=u10cubed(:ncol)**3.41_r8
1143 :
1144 0 : sflx(:)=0._r8
1145 :
1146 0 : call seasalt_emis( u10cubed, cam_in%sst, cam_in%ocnfrac, ncol, cam_in%cflx )
1147 :
1148 0 : do m=1,seasalt_nbin
1149 0 : mm = seasalt_indices(m)
1150 0 : sflx(:ncol)=sflx(:ncol)+cam_in%cflx(:ncol,mm)
1151 0 : call outfld(trim(seasalt_names(m))//'SF',cam_in%cflx(:,mm),pcols,lchnk)
1152 : enddo
1153 0 : call outfld('SSTSFMBL',sflx(:),pcols,lchnk)
1154 : endif
1155 :
1156 0 : end subroutine aero_model_emissions
1157 :
1158 : end module aero_model
|
