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