Line data Source code
1 : module aero_wetdep_cam
2 : use shr_kind_mod, only: r8 => shr_kind_r8
3 : use physics_types, only: physics_state, physics_ptend, physics_ptend_init
4 : use camsrfexch, only: cam_out_t
5 : use physics_buffer,only: physics_buffer_desc, pbuf_get_index, pbuf_set_field, pbuf_get_field
6 : use constituents, only: pcnst, cnst_name, cnst_get_ind
7 : use phys_control, only: phys_getopts
8 : use ppgrid, only: pcols, pver
9 : use physconst, only: gravit
10 :
11 : use cam_abortutils,only: endrun
12 : use cam_logfile, only: iulog
13 : use spmd_utils, only: masterproc
14 : use infnan, only: nan, assignment(=)
15 :
16 : use cam_history, only: addfld, add_default, horiz_only, outfld
17 : use wetdep, only: wetdep_init
18 :
19 : use rad_constituents, only: rad_cnst_get_info
20 :
21 : use aerosol_properties_mod, only: aero_name_len
22 : use aerosol_properties_mod, only: aerosol_properties
23 : use modal_aerosol_properties_mod, only: modal_aerosol_properties
24 :
25 : use aerosol_state_mod, only: aerosol_state, ptr2d_t
26 : use modal_aerosol_state_mod, only: modal_aerosol_state
27 :
28 : use aero_convproc, only: aero_convproc_readnl, aero_convproc_init, aero_convproc_intr
29 : use aero_convproc, only: convproc_do_evaprain_atonce
30 : use aero_convproc, only: deepconv_wetdep_history
31 :
32 : use infnan, only: nan, assignment(=)
33 : use perf_mod, only: t_startf, t_stopf
34 :
35 : implicit none
36 : private
37 :
38 : public :: aero_wetdep_readnl
39 : public :: aero_wetdep_init
40 : public :: aero_wetdep_tend
41 :
42 : real(r8), parameter :: NOTSET = -huge(1._r8)
43 : real(r8) :: sol_facti_cloud_borne = NOTSET
44 : real(r8) :: sol_factb_interstitial = NOTSET
45 : real(r8) :: sol_factic_interstitial = NOTSET
46 :
47 : integer :: fracis_idx = -1
48 : integer :: rprddp_idx = -1
49 : integer :: rprdsh_idx = -1
50 : integer :: nevapr_shcu_idx = -1
51 : integer :: nevapr_dpcu_idx = -1
52 :
53 : logical :: wetdep_active = .false.
54 : integer :: nwetdep = 0
55 : logical :: convproc_do_aer = .false.
56 : logical,allocatable :: aero_cnst_lq(:,:)
57 : integer,allocatable :: aero_cnst_id(:,:)
58 : logical, public, protected :: wetdep_lq(pcnst) ! set flags true for constituents with non-zero tendencies
59 :
60 : ! variables for table lookup of aerosol impaction/interception scavenging rates
61 : integer, parameter :: nimptblgrow_mind=-7, nimptblgrow_maxd=12
62 : real(r8) :: dlndg_nimptblgrow
63 : real(r8),allocatable :: scavimptblnum(:,:)
64 : real(r8),allocatable :: scavimptblvol(:,:)
65 :
66 : integer :: nmodes=0
67 : integer :: nspec_max=0
68 : integer :: nele_tot ! total number of aerosol elements
69 : class(aerosol_properties), pointer :: aero_props=>null()
70 :
71 : contains
72 :
73 : !------------------------------------------------------------------------------
74 : !------------------------------------------------------------------------------
75 0 : subroutine aero_wetdep_readnl(nlfile)
76 : use namelist_utils, only: find_group_name
77 : use spmd_utils, only: mpicom, masterprocid, mpi_character, mpi_real8, mpi_integer, mpi_success
78 : use spmd_utils, only: mpi_logical
79 :
80 : character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input
81 :
82 : integer :: unitn, ierr
83 : character(len=*), parameter :: subname = 'aero_wetdep_readnl'
84 :
85 : ! ===================
86 : ! Namelist definition
87 : ! ===================
88 : namelist /aero_wetdep_nl/ sol_facti_cloud_borne, sol_factb_interstitial, sol_factic_interstitial
89 :
90 : ! =============
91 : ! Read namelist
92 : ! =============
93 0 : if (masterproc) then
94 0 : open( newunit=unitn, file=trim(nlfile), status='old' )
95 0 : call find_group_name(unitn, 'aero_wetdep_nl', status=ierr)
96 0 : if (ierr == 0) then
97 0 : read(unitn, aero_wetdep_nl, iostat=ierr)
98 0 : if (ierr /= 0) then
99 0 : call endrun(subname // ':: ERROR reading namelist')
100 : end if
101 : end if
102 0 : close(unitn)
103 :
104 : ! ============================
105 : ! Log namelist options
106 : ! ============================
107 0 : write(iulog,*) subname,' namelist settings: '
108 0 : write(iulog,*) ' sol_facti_cloud_borne : ',sol_facti_cloud_borne
109 0 : write(iulog,*) ' sol_factb_interstitial : ',sol_factb_interstitial
110 0 : write(iulog,*) ' sol_factic_interstitial: ',sol_factic_interstitial
111 : end if
112 :
113 : ! ============================
114 : ! Broadcast namelist variables
115 : ! ============================
116 0 : call mpi_bcast(sol_facti_cloud_borne, 1, mpi_real8, masterprocid, mpicom, ierr)
117 0 : if (ierr/=mpi_success) then
118 0 : call endrun(subname//': MPI_BCAST ERROR: sol_facti_cloud_borne')
119 : end if
120 0 : call mpi_bcast(sol_factb_interstitial, 1, mpi_real8, masterprocid, mpicom, ierr)
121 0 : if (ierr/=mpi_success) then
122 0 : call endrun(subname//': MPI_BCAST ERROR: sol_factb_interstitial')
123 : end if
124 0 : call mpi_bcast(sol_factic_interstitial, 1, mpi_real8, masterprocid, mpicom, ierr)
125 0 : if (ierr/=mpi_success) then
126 0 : call endrun(subname//': MPI_BCAST ERROR: sol_factic_interstitial')
127 : end if
128 :
129 0 : call mpi_bcast(nwetdep, 1, mpi_integer, masterprocid, mpicom, ierr)
130 0 : if (ierr/=mpi_success) then
131 0 : call endrun(subname//': MPI_BCAST ERROR: nwetdep')
132 : end if
133 :
134 0 : wetdep_active = .true. !nwetdep>0
135 :
136 0 : if (masterproc) then
137 0 : write(iulog,*) subname,' wetdep_active = ',wetdep_active,' nwetdep = ',nwetdep
138 : endif
139 :
140 0 : call aero_convproc_readnl(nlfile)
141 :
142 0 : end subroutine aero_wetdep_readnl
143 :
144 : !------------------------------------------------------------------------------
145 : !------------------------------------------------------------------------------
146 0 : subroutine aero_wetdep_init( )
147 :
148 : character(len=*), parameter :: subrname = 'aero_wetdep_init'
149 :
150 : character(len=2) :: unit_basename ! Units 'kg' or '1'
151 : character(len=aero_name_len) :: tmpname
152 : character(len=aero_name_len) :: tmpname_cw
153 :
154 : logical :: history_aerosol ! Output MAM or SECT aerosol tendencies
155 : logical :: history_chemistry
156 :
157 : integer :: l,m, id, astat
158 : character(len=2) :: binstr
159 :
160 0 : fracis_idx = pbuf_get_index('FRACIS')
161 0 : rprddp_idx = pbuf_get_index('RPRDDP')
162 0 : rprdsh_idx = pbuf_get_index('RPRDSH')
163 0 : nevapr_shcu_idx = pbuf_get_index('NEVAPR_SHCU')
164 0 : nevapr_dpcu_idx = pbuf_get_index('NEVAPR_DPCU')
165 :
166 0 : if (.not.wetdep_active) return
167 :
168 : call phys_getopts(history_aerosol_out = history_aerosol, &
169 : history_chemistry_out=history_chemistry, &
170 : convproc_do_aer_out = convproc_do_aer)
171 :
172 0 : call rad_cnst_get_info(0, nmodes=nmodes)
173 :
174 0 : if (nmodes>0) then
175 0 : aero_props => modal_aerosol_properties()
176 0 : if (.not.associated(aero_props)) then
177 0 : call endrun(subrname//' : construction of aero_props modal_aerosol_properties object failed')
178 : end if
179 : else
180 0 : call endrun(subrname//' : cannot determine aerosol model')
181 : endif
182 :
183 0 : nele_tot = aero_props%ncnst_tot()
184 :
185 0 : allocate(aero_cnst_lq(aero_props%nbins(),0:maxval(aero_props%nmasses())), stat=astat)
186 0 : if (astat/=0) then
187 0 : call endrun(subrname//' : not able to allocate aero_cnst_lq array')
188 : end if
189 0 : aero_cnst_lq(:,:) = .false.
190 :
191 0 : allocate(aero_cnst_id(aero_props%nbins(),0:maxval(aero_props%nmasses())), stat=astat)
192 0 : if (astat/=0) then
193 0 : call endrun(subrname//' : not able to allocate aero_cnst_id array')
194 : end if
195 0 : aero_cnst_id(:,:) = -1
196 :
197 0 : wetdep_lq = .false.
198 :
199 0 : do m = 1, aero_props%nbins()
200 0 : write(binstr,'(i2.2)') m
201 0 : call addfld('SOLFACTB'//binstr, (/ 'lev' /), 'A', '1', 'below cld sol fact')
202 :
203 0 : do l = 0, aero_props%nmasses(m)
204 :
205 0 : if (l == 0) then ! number
206 0 : call aero_props%num_names( m, tmpname, tmpname_cw)
207 : else
208 0 : call aero_props%mmr_names( m,l, tmpname, tmpname_cw)
209 : end if
210 :
211 0 : call cnst_get_ind(tmpname, id, abort=.false.)
212 0 : aero_cnst_id(m,l) = id
213 0 : aero_cnst_lq(m,l) = id > 0
214 0 : if (id > 0) then
215 0 : wetdep_lq(id) = .true.
216 : end if
217 :
218 : ! units --
219 0 : if (l==0) then
220 0 : unit_basename = ' 1' ! for num
221 : else
222 0 : unit_basename = 'kg'
223 : endif
224 :
225 0 : call add_hist_fields(tmpname, unit_basename)
226 0 : call add_hist_fields(tmpname_cw, unit_basename)
227 :
228 : call addfld( trim(tmpname_cw)//'RSPTD', (/ 'lev' /), 'A', unit_basename//'/kg/s', &
229 0 : trim(tmpname_cw)//' resuspension tendency')
230 :
231 : end do
232 : end do
233 :
234 0 : allocate(scavimptblnum(nimptblgrow_mind:nimptblgrow_maxd, aero_props%nbins()), stat=astat)
235 0 : if (astat/=0) then
236 0 : call endrun(subrname//' : not able to allocate scavimptblnum array')
237 : end if
238 0 : allocate(scavimptblvol(nimptblgrow_mind:nimptblgrow_maxd, aero_props%nbins()), stat=astat)
239 0 : if (astat/=0) then
240 0 : call endrun(subrname//' : not able to allocate scavimptblvol array')
241 : end if
242 0 : scavimptblnum = nan
243 0 : scavimptblvol = nan
244 :
245 0 : call wetdep_init()
246 :
247 0 : nspec_max = maxval(aero_props%nspecies()) + 2
248 :
249 0 : call init_bcscavcoef()
250 :
251 0 : if (convproc_do_aer) then
252 0 : call aero_convproc_init(aero_props)
253 : end if
254 :
255 : contains
256 :
257 0 : subroutine add_hist_fields(name,baseunits)
258 : character(len=*), intent(in) :: name
259 : character(len=*), intent(in) :: baseunits
260 :
261 : call addfld (trim(name)//'SFWET', &
262 0 : horiz_only, 'A',baseunits//'/m2/s ','Wet deposition flux at surface')
263 : call addfld (trim(name)//'SFSIC', &
264 0 : horiz_only, 'A',baseunits//'/m2/s ','Wet deposition flux (incloud, convective) at surface')
265 : call addfld (trim(name)//'SFSIS', &
266 0 : horiz_only, 'A',baseunits//'/m2/s ','Wet deposition flux (incloud, stratiform) at surface')
267 : call addfld (trim(name)//'SFSBC', &
268 0 : horiz_only, 'A',baseunits//'/m2/s ','Wet deposition flux (belowcloud, convective) at surface')
269 : call addfld (trim(name)//'SFSBS', &
270 0 : horiz_only, 'A',baseunits//'/m2/s ','Wet deposition flux (belowcloud, stratiform) at surface')
271 :
272 0 : if (convproc_do_aer) then
273 : call addfld (trim(name)//'SFSEC', &
274 0 : horiz_only, 'A',unit_basename//'/m2/s','Wet deposition flux (precip evap, convective) at surface')
275 : call addfld (trim(name)//'SFSES', &
276 0 : horiz_only, 'A',unit_basename//'/m2/s','Wet deposition flux (precip evap, stratiform) at surface')
277 : call addfld (trim(name)//'SFSBD', &
278 0 : horiz_only, 'A',unit_basename//'/m2/s','Wet deposition flux (belowcloud, deep convective) at surface')
279 : call addfld (trim(name)//'WETC', &
280 0 : (/ 'lev' /), 'A',unit_basename//'/kg/s ','wet deposition tendency')
281 : call addfld (trim(name)//'CONU', &
282 0 : (/ 'lev' /), 'A',unit_basename//'/kg ','updraft mixing ratio')
283 : end if
284 :
285 0 : call addfld (trim(name)//'WET',(/ 'lev' /), 'A',baseunits//'/kg/s ','wet deposition tendency')
286 0 : call addfld (trim(name)//'INS',(/ 'lev' /), 'A',baseunits//'/kg/s ','insol frac')
287 :
288 : call addfld (trim(name)//'SIC',(/ 'lev' /), 'A',baseunits//'/kg/s ', &
289 0 : trim(name)//' ic wet deposition')
290 : call addfld (trim(name)//'SIS',(/ 'lev' /), 'A',baseunits//'/kg/s ', &
291 0 : trim(name)//' is wet deposition')
292 : call addfld (trim(name)//'SBC',(/ 'lev' /), 'A',baseunits//'/kg/s ', &
293 0 : trim(name)//' bc wet deposition')
294 : call addfld (trim(name)//'SBS',(/ 'lev' /), 'A',baseunits//'/kg/s ', &
295 0 : trim(name)//' bs wet deposition')
296 :
297 0 : if ( history_aerosol .or. history_chemistry ) then
298 0 : call add_default (trim(name)//'SFWET', 1, ' ')
299 : endif
300 0 : if ( history_aerosol ) then
301 0 : call add_default (trim(name)//'SFSEC', 1, ' ')
302 0 : call add_default (trim(name)//'SFSIC', 1, ' ')
303 0 : call add_default (trim(name)//'SFSIS', 1, ' ')
304 0 : call add_default (trim(name)//'SFSBC', 1, ' ')
305 0 : call add_default (trim(name)//'SFSBS', 1, ' ')
306 0 : if (convproc_do_aer) then
307 0 : call add_default (trim(name)//'SFSES', 1, ' ')
308 0 : call add_default (trim(name)//'SFSBD', 1, ' ')
309 : end if
310 : endif
311 :
312 0 : end subroutine add_hist_fields
313 :
314 : end subroutine aero_wetdep_init
315 :
316 : !------------------------------------------------------------------------------
317 : !------------------------------------------------------------------------------
318 0 : subroutine aero_wetdep_tend( state, dt, dlf, cam_out, ptend, pbuf)
319 : use wetdep, only: wetdepa_v2, wetdep_inputs_set, wetdep_inputs_t
320 : use aerodep_flx, only: aerodep_flx_prescribed
321 : use aero_deposition_cam, only: aero_deposition_cam_setwet
322 :
323 : type(physics_state), target, intent(in) :: state ! Physics state variables
324 : real(r8), intent(in) :: dt ! time step
325 : real(r8), intent(in) :: dlf(:,:) ! shallow+deep convective detrainment [kg/kg/s]
326 : type(cam_out_t), intent(inout) :: cam_out ! export state
327 : type(physics_ptend), intent(out) :: ptend ! indivdual parameterization tendencies
328 : type(physics_buffer_desc), pointer :: pbuf(:)
329 :
330 : character(len=*), parameter :: subrname = 'aero_wetdep_tend'
331 : type(wetdep_inputs_t) :: dep_inputs
332 0 : real(r8), pointer :: fracis(:,:,:) ! fraction of transported species that are insoluble (pcols, pver, pcnst)
333 : real(r8), target :: fracis_nadv(pcols,pver) ! fraction of not-transported aerosols
334 :
335 : real(r8) :: scavcoefnv(pcols,pver,0:2) ! Dana and Hales coefficient (/mm) for
336 : ! cloud-borne num & vol (0),
337 : ! interstitial num (1), interstitial vol (2)
338 : integer :: jnv ! index for scavcoefnv 3rd dimension
339 : integer :: lphase ! index for interstitial / cloudborne aerosol
340 : integer :: strt_loop, end_loop, stride_loop !loop indices for the lphase loop
341 :
342 : real(r8) :: sol_factb(pcols, pver)
343 : real(r8) :: sol_facti(pcols, pver)
344 : real(r8) :: sol_factic(pcols,pver)
345 :
346 : real(r8) :: dqdt_tmp(pcols,pver) ! temporary array to hold tendency for 1 species
347 : real(r8) :: rcscavt(pcols, pver)
348 : real(r8) :: rsscavt(pcols, pver)
349 : real(r8) :: iscavt(pcols, pver)
350 : real(r8) :: icscavt(pcols, pver)
351 : real(r8) :: isscavt(pcols, pver)
352 : real(r8) :: bcscavt(pcols, pver)
353 : real(r8) :: bsscavt(pcols, pver)
354 :
355 0 : real(r8) :: diam_wet(state%ncol, pver)
356 : logical :: isprx(pcols,pver) ! true if precipation
357 : real(r8) :: prec(pcols) ! precipitation rate
358 :
359 0 : real(r8) :: rtscavt(pcols, pver, 0:nspec_max)
360 :
361 : integer :: ncol, lchnk, m, ndx,mm, l
362 : integer :: i,k
363 :
364 0 : real(r8), pointer :: insolfr_ptr(:,:)
365 : real(r8) :: q_tmp(pcols,pver) ! temporary array to hold "most current" mixing ratio for 1 species
366 : logical :: cldbrn
367 :
368 0 : type(ptr2d_t) :: raer(nele_tot)
369 0 : type(ptr2d_t) :: qqcw(nele_tot)
370 :
371 : real(r8) :: sflx(pcols)
372 : character(len=aero_name_len) :: aname, cname, name
373 :
374 0 : real(r8) :: qqcw_in(pcols,pver), qqcw_sav(pcols,pver,0:nspec_max)
375 : real(r8) :: f_act_conv(pcols,pver) ! prescribed aerosol activation fraction for convective cloud ! rce 2010/05/01
376 :
377 : character(len=2) :: binstr
378 : real(r8) :: aerdepwetcw(pcols,pcnst) ! aerosol wet deposition (cloud water)
379 : real(r8) :: aerdepwetis(pcols,pcnst) ! aerosol wet deposition (interstitial)
380 0 : real(r8) :: dcondt_resusp3d(nele_tot,pcols,pver)
381 :
382 : integer, parameter :: nsrflx_mzaer2cnvpr = 2
383 0 : real(r8) :: qsrflx_mzaer2cnvpr(pcols,nele_tot,nsrflx_mzaer2cnvpr)
384 :
385 0 : real(r8), pointer :: rprddp(:,:) ! rain production, deep convection
386 0 : real(r8), pointer :: rprdsh(:,:) ! rain production, shallow convection
387 0 : real(r8), pointer :: evapcdp(:,:) ! Evaporation rate of deep convective precipitation >=0.
388 0 : real(r8), pointer :: evapcsh(:,:) ! Evaporation rate of shallow convective precipitation >=0.
389 :
390 : real(r8) :: rprddpsum(pcols)
391 : real(r8) :: rprdshsum(pcols)
392 : real(r8) :: evapcdpsum(pcols)
393 : real(r8) :: evapcshsum(pcols)
394 :
395 : real(r8) :: tmp_resudp, tmp_resush
396 : real(r8) :: tmpa, tmpb
397 : real(r8) :: sflxec(pcols), sflxecdp(pcols) ! deposition flux
398 : real(r8) :: sflxic(pcols), sflxicdp(pcols) ! deposition flux
399 : real(r8) :: sflxbc(pcols), sflxbcdp(pcols) ! deposition flux
400 :
401 : class(aerosol_state), pointer :: aero_state
402 :
403 0 : nullify(aero_state)
404 :
405 0 : if (.not.wetdep_active) return
406 :
407 0 : dcondt_resusp3d(:,:,:) = 0._r8
408 :
409 0 : if (nmodes>0) then
410 0 : aero_state => modal_aerosol_state(state,pbuf)
411 0 : if (.not.associated(aero_state)) then
412 0 : call endrun(subrname//' : construction of aero_state modal_aerosol_state object failed')
413 : end if
414 : else
415 0 : call endrun(subrname//' : cannot determine aerosol model')
416 : endif
417 :
418 0 : lchnk = state%lchnk
419 0 : ncol = state%ncol
420 :
421 0 : call physics_ptend_init(ptend, state%psetcols, subrname, lq=wetdep_lq)
422 :
423 0 : call wetdep_inputs_set( state, pbuf, dep_inputs )
424 :
425 0 : call pbuf_get_field(pbuf, fracis_idx, fracis)
426 :
427 0 : call aero_state%get_states( aero_props, raer, qqcw )
428 :
429 0 : qsrflx_mzaer2cnvpr(:,:,:) = 0.0_r8
430 0 : aerdepwetis(:,:) = 0.0_r8
431 0 : aerdepwetcw(:,:) = 0.0_r8
432 :
433 0 : if (convproc_do_aer) then
434 : !Do cloudborne first for unified convection scheme so that the resuspension of cloudborne
435 : !can be saved then applied to interstitial
436 : strt_loop = 2
437 : end_loop = 1
438 : stride_loop = -1
439 : else
440 : ! Counters for "without" unified convective treatment (i.e. default case)
441 0 : strt_loop = 1
442 0 : end_loop = 2
443 0 : stride_loop = 1
444 : endif
445 :
446 0 : prec(:ncol)=0._r8
447 0 : do k=1,pver
448 0 : where (prec(:ncol) >= 1.e-7_r8)
449 : isprx(:ncol,k) = .true.
450 : elsewhere
451 0 : isprx(:ncol,k) = .false.
452 : endwhere
453 0 : prec(:ncol) = prec(:ncol) + (dep_inputs%prain(:ncol,k) + dep_inputs%cmfdqr(:ncol,k) - dep_inputs%evapr(:ncol,k)) &
454 0 : *state%pdel(:ncol,k)/gravit
455 : end do
456 :
457 0 : f_act_conv = 0._r8
458 0 : scavcoefnv = nan
459 0 : qqcw_sav = nan
460 :
461 0 : if (convproc_do_aer) then
462 :
463 0 : call t_startf('aero_convproc')
464 : call aero_convproc_intr( aero_props, aero_state, state, ptend, pbuf, dt, &
465 0 : nsrflx_mzaer2cnvpr, qsrflx_mzaer2cnvpr, aerdepwetis, dcondt_resusp3d )
466 :
467 0 : if (convproc_do_evaprain_atonce) then
468 :
469 0 : do m = 1,aero_props%nbins()
470 0 : do l = 0,aero_props%nmasses(m)
471 0 : mm = aero_props%indexer(m,l)
472 :
473 0 : if (l == 0) then ! number
474 0 : call aero_props%num_names(m, aname, cname)
475 : else
476 0 : call aero_props%mmr_names(m,l, aname, cname)
477 : end if
478 :
479 0 : call outfld( trim(cname)//'RSPTD', dcondt_resusp3d(mm,:ncol,:), ncol, lchnk )
480 :
481 0 : do k = 1,pver
482 0 : do i = 1,ncol
483 0 : qqcw(mm)%fld(i,k) = max(0._r8, qqcw(mm)%fld(i,k) + dcondt_resusp3d(mm,i,k)*dt)
484 : end do
485 : end do
486 :
487 : end do
488 : end do
489 : end if
490 0 : call t_stopf('aero_convproc')
491 :
492 : end if
493 :
494 0 : bins_loop: do m = 1,aero_props%nbins()
495 :
496 0 : phase_loop: do lphase = strt_loop,end_loop, stride_loop ! loop over interstitial (1) and cloud-borne (2) forms
497 :
498 0 : cldbrn = lphase==2
499 :
500 0 : sol_factb = nan
501 0 : sol_facti = nan
502 0 : sol_factic = nan
503 :
504 0 : if (lphase == 1) then ! interstial aerosol
505 :
506 0 : sol_facti = 0.0_r8 ! strat in-cloud scav totally OFF for institial
507 :
508 0 : sol_factic = sol_factic_interstitial
509 :
510 : else ! cloud-borne aerosol (borne by stratiform cloud drops)
511 :
512 0 : sol_factb = 0.0_r8 ! all below-cloud scav OFF (anything cloud-borne is located "in-cloud")
513 0 : sol_facti = sol_facti_cloud_borne ! strat in-cloud scav cloud-borne tuning factor
514 0 : sol_factic = 0.0_r8 ! conv in-cloud scav OFF (having this on would mean
515 : ! that conv precip collects strat droplets)
516 0 : f_act_conv = 0.0_r8 ! conv in-cloud scav OFF (having this on would mean
517 :
518 : end if
519 0 : if (convproc_do_aer .and. lphase == 1) then
520 : ! if modal aero convproc is turned on for aerosols, then
521 : ! turn off the convective in-cloud removal for interstitial aerosols
522 : ! (but leave the below-cloud on, as convproc only does in-cloud)
523 : ! and turn off the outfld SFWET, SFSIC, SFSID, SFSEC, and SFSED calls
524 : ! for (stratiform)-cloudborne aerosols, convective wet removal
525 : ! (all forms) is zero, so no action is needed
526 0 : sol_factic = 0.0_r8
527 : endif
528 :
529 0 : diam_wet = aero_state%wet_diameter(m,ncol,pver)
530 :
531 0 : scavcoefnv = 0.0_r8
532 :
533 0 : if (lphase == 1) then ! interstial aerosol
534 0 : call get_bcscavcoefs( m, ncol, isprx, diam_wet, scavcoefnv(:,:,1), scavcoefnv(:,:,2) )
535 :
536 0 : if ( sol_factb_interstitial /= NOTSET ) then
537 0 : sol_factb(:ncol,:) = sol_factb_interstitial ! all below-cloud scav
538 : else
539 0 : sol_factb(:ncol,:) = aero_state%sol_factb_interstitial( m, ncol, pver, aero_props )
540 : end if
541 :
542 0 : write(binstr,'(i2.2)') m
543 0 : call outfld('SOLFACTB'//binstr,sol_factb, pcols, lchnk)
544 :
545 : end if
546 :
547 0 : masses_loop: do l = 0,aero_props%nmasses(m)
548 :
549 0 : ndx = aero_cnst_id(m,l)
550 :
551 0 : if (.not. cldbrn .and. ndx>0) then
552 0 : insolfr_ptr => fracis(:,:,ndx)
553 : else
554 0 : insolfr_ptr => fracis_nadv
555 : endif
556 :
557 0 : mm = aero_props%indexer(m,l)
558 :
559 0 : if (l == 0) then ! number
560 0 : call aero_props%num_names( m, aname, cname)
561 : else
562 0 : call aero_props%mmr_names( m,l, aname, cname)
563 : end if
564 :
565 0 : if (cldbrn) then
566 0 : q_tmp(1:ncol,:) = qqcw(mm)%fld(1:ncol,:)
567 0 : jnv = 0
568 0 : if (convproc_do_aer) then
569 0 : qqcw_sav(:ncol,:,l) = q_tmp(1:ncol,:)
570 : endif
571 0 : name = cname
572 0 : qqcw_in = nan
573 0 : f_act_conv = nan
574 : else ! interstial aerosol
575 0 : q_tmp(1:ncol,:) = raer(mm)%fld(1:ncol,:) + ptend%q(1:ncol,:,ndx)*dt
576 0 : if (l==0) then
577 : jnv = 1
578 : else
579 0 : jnv = 2
580 : end if
581 0 : if(convproc_do_aer) then
582 : !Feed in the saved cloudborne mixing ratios from phase 2
583 0 : qqcw_in(:ncol,:) = qqcw_sav(:ncol,:,l)
584 : else
585 0 : qqcw_in(:ncol,:) = qqcw(mm)%fld(:ncol,:)
586 : end if
587 :
588 0 : f_act_conv(:ncol,:) = aero_state%convcld_actfrac( m, l, ncol, pver)
589 0 : name = aname
590 : end if
591 :
592 0 : dqdt_tmp(1:ncol,:) = 0.0_r8
593 :
594 : call wetdepa_v2(state%pmid, state%q(:,:,1), state%pdel, &
595 : dep_inputs%cldt, dep_inputs%cldcu, dep_inputs%cmfdqr, &
596 : dep_inputs%evapc, dep_inputs%conicw, dep_inputs%prain, dep_inputs%qme, &
597 : dep_inputs%evapr, dep_inputs%totcond, q_tmp, dt, &
598 : dqdt_tmp, iscavt, dep_inputs%cldvcu, dep_inputs%cldvst, &
599 : dlf, insolfr_ptr, sol_factb, ncol, &
600 : scavcoefnv(:,:,jnv), &
601 : is_strat_cloudborne=cldbrn, &
602 : qqcw=qqcw_in(:,:), f_act_conv=f_act_conv, &
603 : icscavt=icscavt, isscavt=isscavt, bcscavt=bcscavt, bsscavt=bsscavt, &
604 : convproc_do_aer=convproc_do_aer, rcscavt=rcscavt, rsscavt=rsscavt, &
605 : sol_facti_in=sol_facti, sol_factic_in=sol_factic, &
606 : convproc_do_evaprain_atonce_in=convproc_do_evaprain_atonce, &
607 0 : bergso_in=dep_inputs%bergso )
608 :
609 0 : if(convproc_do_aer) then
610 0 : if(cldbrn) then
611 : ! save resuspension of cloudborne species
612 0 : rtscavt(1:ncol,:,l) = rcscavt(1:ncol,:) + rsscavt(1:ncol,:)
613 : ! wetdepa_v2 adds the resuspension of cloudborne to the dqdt of cloudborne (as a source)
614 : ! undo this, so the resuspension of cloudborne can be added to the dqdt of interstitial (above)
615 0 : dqdt_tmp(1:ncol,:) = dqdt_tmp(1:ncol,:) - rtscavt(1:ncol,:,l)
616 : else
617 : ! add resuspension of cloudborne species to dqdt of interstitial species
618 0 : dqdt_tmp(1:ncol,:) = dqdt_tmp(1:ncol,:) + rtscavt(1:ncol,:,l)
619 : end if
620 : endif
621 :
622 0 : if (cldbrn .or. ndx<0) then
623 0 : do k = 1,pver
624 0 : do i = 1,ncol
625 0 : if ( (qqcw(mm)%fld(i,k) + dqdt_tmp(i,k) * dt) .lt. 0.0_r8 ) then
626 0 : dqdt_tmp(i,k) = - qqcw(mm)%fld(i,k) / dt
627 : end if
628 : end do
629 : end do
630 :
631 0 : qqcw(mm)%fld(1:ncol,:) = qqcw(mm)%fld(1:ncol,:) + dqdt_tmp(1:ncol,:) * dt
632 :
633 : else
634 0 : ptend%q(1:ncol,:,ndx) = ptend%q(1:ncol,:,ndx) + dqdt_tmp(1:ncol,:)
635 : end if
636 :
637 0 : call outfld( trim(name)//'WET', dqdt_tmp(:,:), pcols, lchnk)
638 0 : call outfld( trim(name)//'SIC', icscavt, pcols, lchnk)
639 0 : call outfld( trim(name)//'SIS', isscavt, pcols, lchnk)
640 0 : call outfld( trim(name)//'SBC', bcscavt, pcols, lchnk)
641 0 : call outfld( trim(name)//'SBS', bsscavt, pcols, lchnk)
642 :
643 0 : call outfld( trim(name)//'INS', insolfr_ptr, pcols, lchnk)
644 :
645 0 : sflx(:)=0._r8
646 0 : do k=1,pver
647 0 : do i=1,ncol
648 0 : sflx(i)=sflx(i)+dqdt_tmp(i,k)*state%pdel(i,k)/gravit
649 : enddo
650 : enddo
651 0 : if (cldbrn) then
652 0 : call outfld( trim(name)//'SFWET', sflx, pcols, lchnk)
653 0 : if (ndx>0) aerdepwetcw(:ncol,ndx) = sflx(:ncol)
654 : else
655 0 : if (.not.convproc_do_aer) call outfld( trim(name)//'SFWET', sflx, pcols, lchnk)
656 0 : if (ndx>0) aerdepwetis(:ncol,ndx) = aerdepwetis(:ncol,ndx) + sflx(:ncol)
657 : end if
658 :
659 0 : sflx(:)=0._r8
660 0 : do k=1,pver
661 0 : do i=1,ncol
662 0 : sflx(i)=sflx(i)+icscavt(i,k)*state%pdel(i,k)/gravit
663 : enddo
664 : enddo
665 0 : if (cldbrn) then
666 0 : call outfld( trim(name)//'SFSIC', sflx, pcols, lchnk)
667 : else
668 0 : if (.not.convproc_do_aer) call outfld( trim(name)//'SFSIC', sflx, pcols, lchnk)
669 0 : if (convproc_do_aer) sflxic = sflx
670 : end if
671 :
672 0 : sflx(:)=0._r8
673 0 : do k=1,pver
674 0 : do i=1,ncol
675 0 : sflx(i)=sflx(i)+isscavt(i,k)*state%pdel(i,k)/gravit
676 : enddo
677 : enddo
678 0 : call outfld( trim(name)//'SFSIS', sflx, pcols, lchnk)
679 :
680 0 : sflx(:)=0._r8
681 0 : do k=1,pver
682 0 : do i=1,ncol
683 0 : sflx(i)=sflx(i)+bcscavt(i,k)*state%pdel(i,k)/gravit
684 : enddo
685 : enddo
686 0 : call outfld( trim(name)//'SFSBC', sflx, pcols, lchnk)
687 0 : if (convproc_do_aer) sflxbc = sflx
688 :
689 0 : sflx(:)=0._r8
690 0 : do k=1,pver
691 0 : do i=1,ncol
692 0 : sflx(i)=sflx(i)+bsscavt(i,k)*state%pdel(i,k)/gravit
693 : enddo
694 : enddo
695 0 : call outfld( trim(name)//'SFSBS', sflx, pcols, lchnk)
696 :
697 0 : if(convproc_do_aer) then
698 :
699 0 : sflx(:)=0.0_r8
700 0 : do k=1,pver
701 0 : do i=1,ncol
702 0 : sflx(i)=sflx(i)+rsscavt(i,k)*state%pdel(i,k)/gravit
703 : enddo
704 : enddo
705 0 : call outfld( trim(name)//'SFSES', sflx, pcols, lchnk)
706 :
707 0 : sflx(:)=0._r8
708 0 : do k=1,pver
709 0 : do i=1,ncol
710 0 : sflx(i)=sflx(i)+rcscavt(i,k)*state%pdel(i,k)/gravit
711 : enddo
712 : enddo
713 0 : if (.not.convproc_do_aer) call outfld( trim(name)//'SFSEC', sflx, pcols, lchnk)
714 0 : sflxec = sflx
715 :
716 0 : if(.not.cldbrn) then
717 :
718 : ! apportion convective surface fluxes to deep and shallow conv
719 : ! this could be done more accurately in subr wetdepa
720 : ! since deep and shallow rarely occur simultaneously, and these
721 : ! fields are just diagnostics, this approximate method is adequate
722 : ! only do this for interstitial aerosol, because conv clouds to not
723 : ! affect the stratiform-cloudborne aerosol
724 0 : if ( deepconv_wetdep_history) then
725 :
726 0 : call pbuf_get_field(pbuf, rprddp_idx, rprddp )
727 0 : call pbuf_get_field(pbuf, rprdsh_idx, rprdsh )
728 0 : call pbuf_get_field(pbuf, nevapr_dpcu_idx, evapcdp )
729 0 : call pbuf_get_field(pbuf, nevapr_shcu_idx, evapcsh )
730 :
731 0 : rprddpsum(:) = 0.0_r8
732 0 : rprdshsum(:) = 0.0_r8
733 0 : evapcdpsum(:) = 0.0_r8
734 0 : evapcshsum(:) = 0.0_r8
735 :
736 0 : do k = 1, pver
737 0 : rprddpsum(:ncol) = rprddpsum(:ncol) + rprddp(:ncol,k)*state%pdel(:ncol,k)/gravit
738 0 : rprdshsum(:ncol) = rprdshsum(:ncol) + rprdsh(:ncol,k)*state%pdel(:ncol,k)/gravit
739 0 : evapcdpsum(:ncol) = evapcdpsum(:ncol) + evapcdp(:ncol,k)*state%pdel(:ncol,k)/gravit
740 0 : evapcshsum(:ncol) = evapcshsum(:ncol) + evapcsh(:ncol,k)*state%pdel(:ncol,k)/gravit
741 : end do
742 :
743 0 : do i = 1, ncol
744 0 : rprddpsum(i) = max( rprddpsum(i), 1.0e-35_r8 )
745 0 : rprdshsum(i) = max( rprdshsum(i), 1.0e-35_r8 )
746 0 : evapcdpsum(i) = max( evapcdpsum(i), 0.1e-35_r8 )
747 0 : evapcshsum(i) = max( evapcshsum(i), 0.1e-35_r8 )
748 :
749 : ! assume that in- and below-cloud removal are proportional to column precip production
750 0 : tmpa = rprddpsum(i) / (rprddpsum(i) + rprdshsum(i))
751 0 : tmpa = max( 0.0_r8, min( 1.0_r8, tmpa ) )
752 0 : sflxicdp(i) = sflxic(i)*tmpa
753 0 : sflxbcdp(i) = sflxbc(i)*tmpa
754 :
755 : ! assume that resuspension is proportional to (wet removal)*[(precip evap)/(precip production)]
756 0 : tmp_resudp = tmpa * min( (evapcdpsum(i)/rprddpsum(i)), 1.0_r8 )
757 0 : tmp_resush = (1.0_r8 - tmpa) * min( (evapcshsum(i)/rprdshsum(i)), 1.0_r8 )
758 0 : tmpb = max( tmp_resudp, 1.0e-35_r8 ) / max( (tmp_resudp+tmp_resush), 1.0e-35_r8 )
759 0 : tmpb = max( 0.0_r8, min( 1.0_r8, tmpb ) )
760 0 : sflxecdp(i) = sflxec(i)*tmpb
761 : end do
762 0 : call outfld( trim(name)//'SFSBD', sflxbcdp, pcols, lchnk)
763 : else
764 0 : sflxec(1:ncol) = 0.0_r8
765 0 : sflxecdp(1:ncol) = 0.0_r8
766 : end if
767 :
768 : ! when ma_convproc_intr is used, convective in-cloud wet removal is done there
769 : ! the convective (total and deep) precip-evap-resuspension includes in- and below-cloud
770 : ! contributions
771 : ! so pass the below-cloud contribution to ma_convproc_intr
772 0 : qsrflx_mzaer2cnvpr(1:ncol,mm,1) = sflxec( 1:ncol)
773 0 : qsrflx_mzaer2cnvpr(1:ncol,mm,2) = sflxecdp(1:ncol)
774 :
775 : end if
776 : end if
777 :
778 : end do masses_loop
779 : end do phase_loop
780 :
781 : end do bins_loop
782 :
783 0 : if (associated(aero_state)) then
784 0 : deallocate(aero_state)
785 : nullify(aero_state)
786 : end if
787 :
788 : ! if the user has specified prescribed aerosol dep fluxes then
789 : ! do not set cam_out dep fluxes according to the prognostic aerosols
790 0 : if (.not. aerodep_flx_prescribed()) then
791 0 : call aero_deposition_cam_setwet(aerdepwetis, aerdepwetcw, cam_out)
792 : endif
793 :
794 : contains
795 :
796 : ! below cloud impaction scavenging coefs
797 0 : subroutine get_bcscavcoefs( m, ncol, isprx, diam_wet, scavcoefnum, scavcoefvol )
798 :
799 : integer,intent(in) :: m, ncol
800 : logical,intent(in):: isprx(:,:)
801 : real(r8), intent(in) :: diam_wet(:,:)
802 : real(r8), intent(out) :: scavcoefnum(:,:), scavcoefvol(:,:)
803 :
804 : integer i, k, jgrow
805 : real(r8) dumdgratio, xgrow, dumfhi, dumflo, scavimpvol, scavimpnum
806 :
807 0 : do k = 1, pver
808 0 : do i = 1, ncol
809 :
810 : ! do only if no precip
811 0 : if ( isprx(i,k) .and. diam_wet(i,k)>0.0_r8) then
812 : !
813 : ! interpolate table values using log of (actual-wet-size)/(base-dry-size)
814 :
815 0 : dumdgratio = diam_wet(i,k)/aero_props%scav_diam(m)
816 0 : if ((dumdgratio >= 0.99_r8) .and. (dumdgratio <= 1.01_r8)) then
817 0 : scavimpvol = scavimptblvol(0,m)
818 0 : scavimpnum = scavimptblnum(0,m)
819 : else
820 0 : xgrow = log( dumdgratio ) / dlndg_nimptblgrow
821 0 : jgrow = int( xgrow )
822 0 : if (xgrow < 0._r8) jgrow = jgrow - 1
823 0 : if (jgrow < nimptblgrow_mind) then
824 : jgrow = nimptblgrow_mind
825 : xgrow = jgrow
826 : else
827 0 : jgrow = min( jgrow, nimptblgrow_maxd-1 )
828 : end if
829 :
830 0 : dumfhi = xgrow - jgrow
831 0 : dumflo = 1._r8 - dumfhi
832 :
833 0 : scavimpvol = dumflo*scavimptblvol(jgrow,m) + &
834 0 : dumfhi*scavimptblvol(jgrow+1,m)
835 0 : scavimpnum = dumflo*scavimptblnum(jgrow,m) + &
836 0 : dumfhi*scavimptblnum(jgrow+1,m)
837 :
838 : end if
839 :
840 : ! impaction scavenging removal amount for volume
841 0 : scavcoefvol(i,k) = exp( scavimpvol )
842 : ! impaction scavenging removal amount to number
843 0 : scavcoefnum(i,k) = exp( scavimpnum )
844 :
845 : else
846 0 : scavcoefvol(i,k) = 0._r8
847 0 : scavcoefnum(i,k) = 0._r8
848 : end if
849 :
850 : end do
851 : end do
852 :
853 0 : end subroutine get_bcscavcoefs
854 :
855 : end subroutine aero_wetdep_tend
856 :
857 : !------------------------------------------------------------------------------
858 : !------------------------------------------------------------------------------
859 0 : subroutine init_bcscavcoef( )
860 : !-----------------------------------------------------------------------
861 : !
862 : ! Purpose:
863 : ! Computes lookup table for aerosol impaction/interception scavenging rates
864 : !
865 : ! Authors: R. Easter
866 : ! Simone Tilmes Nov 2021
867 : ! added modifications for bin model, assuming sigma = 1.
868 : !
869 : !-----------------------------------------------------------------------
870 :
871 : use mo_constants, only: pi
872 :
873 : ! local variables
874 : integer nnfit_maxd
875 : parameter (nnfit_maxd=27)
876 :
877 : integer m, jgrow, nnfit
878 : integer lunerr
879 :
880 : real(r8) dg0, dg0_cgs, press, dg0_base, &
881 : rhodryaero, rhowetaero, rhowetaero_cgs, &
882 : scavratenum, scavratevol, logsig, &
883 : temp, wetdiaratio, wetvolratio
884 :
885 : real(r8) :: xxfitnum(1,nnfit_maxd), yyfitnum(nnfit_maxd)
886 : real(r8) :: xxfitvol(1,nnfit_maxd), yyfitvol(nnfit_maxd)
887 :
888 : character(len=*), parameter :: subname = 'aero_wetdep_cam::init_bcscavcoef'
889 :
890 0 : lunerr = iulog
891 0 : dlndg_nimptblgrow = log( 1.25_r8 )
892 :
893 : ! bin model: main loop over aerosol bins
894 :
895 0 : modeloop: do m = 1, aero_props%nbins()
896 :
897 : ! for setting up the lookup table, use the dry density of the first species
898 : ! -- assume the first species of the mode/bin is the dominate species
899 0 : call aero_props%get(m,1,density=rhodryaero)
900 :
901 0 : dg0_base = aero_props%scav_diam(m)
902 :
903 0 : logsig = aero_props%alogsig(m)
904 :
905 0 : growloop: do jgrow = nimptblgrow_mind, nimptblgrow_maxd
906 :
907 0 : wetdiaratio = exp( jgrow*dlndg_nimptblgrow )
908 0 : dg0 = dg0_base*wetdiaratio
909 :
910 : wetvolratio = exp( jgrow*dlndg_nimptblgrow*3._r8 )
911 0 : rhowetaero = 1.0_r8 + (rhodryaero-1.0_r8)/wetvolratio
912 0 : rhowetaero = min( rhowetaero, rhodryaero )
913 :
914 : !
915 : ! compute impaction scavenging rates at 1 temp-press pair and save
916 : !
917 0 : nnfit = 0
918 :
919 0 : temp = 273.16_r8
920 0 : press = 0.75e6_r8 ! dynes/cm2
921 0 : rhowetaero = rhodryaero
922 :
923 0 : dg0_cgs = dg0*1.0e2_r8 ! m to cm
924 :
925 0 : rhowetaero_cgs = rhowetaero*1.0e-3_r8 ! kg/m3 to g/cm3
926 :
927 : call calc_1_impact_rate( &
928 : dg0_cgs, logsig, rhowetaero_cgs, temp, press, &
929 0 : scavratenum, scavratevol, lunerr )
930 :
931 0 : nnfit = nnfit + 1
932 : if (nnfit > nnfit_maxd) then
933 : write(lunerr,9110)
934 : call endrun(subname//' : nnfit > nnfit_maxd')
935 : end if
936 : 9110 format( '*** subr. init_bcscavcoef -- nnfit too big' )
937 :
938 0 : xxfitnum(1,nnfit) = 1._r8
939 0 : yyfitnum(nnfit) = log( scavratenum )
940 :
941 0 : xxfitvol(1,nnfit) = 1._r8
942 0 : yyfitvol(nnfit) = log( scavratevol )
943 :
944 : !depends on both bins and different species
945 0 : scavimptblnum(jgrow,m) = yyfitnum(1)
946 0 : scavimptblvol(jgrow,m) = yyfitvol(1)
947 :
948 : enddo growloop
949 : enddo modeloop
950 :
951 : contains
952 :
953 : !===============================================================================
954 0 : subroutine calc_1_impact_rate( &
955 : dg0, logsig, rhoaero, temp, press, &
956 : scavratenum, scavratevol, lunerr )
957 : !
958 : ! routine computes a single impaction scavenging rate
959 : ! for precipitation rate of 1 mm/h
960 : !
961 : ! dg0 = geometric mean diameter of aerosol number size distrib. (cm)
962 : ! sigmag = geometric standard deviation of size distrib.
963 : ! rhoaero = density of aerosol particles (g/cm^3)
964 : ! temp = temperature (K)
965 : ! press = pressure (dyne/cm^2)
966 : ! scavratenum = number scavenging rate (1/h)
967 : ! scavratevol = volume or mass scavenging rate (1/h)
968 : ! lunerr = logical unit for error message
969 : !
970 : use mo_constants, only: boltz_cgs, pi, rhowater => rhoh2o_cgs, rgas => rgas_cgs
971 :
972 : implicit none
973 :
974 : ! subr. parameters
975 : integer, intent(in) :: lunerr
976 : real(r8), intent(in) :: dg0, logsig, rhoaero, temp, press
977 : real(r8), intent(out) :: scavratenum, scavratevol
978 :
979 : ! local variables
980 : integer nrainsvmax
981 : parameter (nrainsvmax=50)
982 : real(r8) rrainsv(nrainsvmax), xnumrainsv(nrainsvmax),&
983 : vfallrainsv(nrainsvmax)
984 :
985 : integer naerosvmax
986 : parameter (naerosvmax=51)
987 : real(r8) aaerosv(naerosvmax), &
988 : ynumaerosv(naerosvmax), yvolaerosv(naerosvmax)
989 :
990 : integer i, ja, jr, na, nr
991 : real(r8) a, aerodiffus, aeromass, ag0, airdynvisc, airkinvisc
992 : real(r8) anumsum, avolsum, cair, chi
993 : real(r8) d, dr, dum, dumfuchs, dx
994 : real(r8) ebrown, eimpact, eintercept, etotal, freepath
995 : real(r8) precip, precipmmhr, precipsum
996 : real(r8) r, rainsweepout, reynolds, rhi, rhoair, rlo, rnumsum
997 : real(r8) scavsumnum, scavsumnumbb
998 : real(r8) scavsumvol, scavsumvolbb
999 : real(r8) schmidt, sqrtreynolds, sstar, stokes, sx
1000 : real(r8) taurelax, vfall, vfallstp
1001 : real(r8) x, xg0, xg3, xhi, xlo, xmuwaterair
1002 :
1003 0 : rlo = .005_r8
1004 0 : rhi = .250_r8
1005 0 : dr = 0.005_r8
1006 0 : nr = 1 + nint( (rhi-rlo)/dr )
1007 : if (nr > nrainsvmax) then
1008 : write(lunerr,9110)
1009 : call endrun(subname//' : nr > nrainsvmax')
1010 : end if
1011 :
1012 : 9110 format( '*** subr. calc_1_impact_rate -- nr > nrainsvmax' )
1013 :
1014 0 : precipmmhr = 1.0_r8
1015 0 : precip = precipmmhr/36000._r8
1016 :
1017 0 : ag0 = dg0/2._r8
1018 0 : sx = logsig
1019 0 : xg0 = log( ag0 )
1020 0 : xg3 = xg0 + 3._r8*sx*sx
1021 :
1022 0 : xlo = xg3 - 4._r8*sx
1023 0 : xhi = xg3 + 4._r8*sx
1024 0 : dx = 0.2_r8*sx
1025 :
1026 0 : dx = max( 0.2_r8*sx, 0.01_r8 )
1027 0 : xlo = xg3 - max( 4._r8*sx, 2._r8*dx )
1028 0 : xhi = xg3 + max( 4._r8*sx, 2._r8*dx )
1029 :
1030 0 : na = 1 + nint( (xhi-xlo)/dx )
1031 0 : if (na > naerosvmax) then
1032 0 : write(lunerr,9120)
1033 0 : call endrun(subname//' : na > naerosvmax')
1034 : end if
1035 :
1036 : 9120 format( '*** subr. calc_1_impact_rate -- na > naerosvmax' )
1037 :
1038 : ! air molar density
1039 0 : cair = press/(rgas*temp)
1040 : ! air mass density
1041 0 : rhoair = 28.966_r8*cair
1042 : ! molecular freepath
1043 0 : freepath = 2.8052e-10_r8/cair
1044 : ! air dynamic viscosity
1045 : airdynvisc = 1.8325e-4_r8 * (416.16_r8/(temp+120._r8)) * &
1046 0 : ((temp/296.16_r8)**1.5_r8)
1047 : ! air kinemaic viscosity
1048 0 : airkinvisc = airdynvisc/rhoair
1049 : ! ratio of water viscosity to air viscosity (from Slinn)
1050 0 : xmuwaterair = 60.0_r8
1051 :
1052 : !
1053 : ! compute rain drop number concentrations
1054 : ! rrainsv = raindrop radius (cm)
1055 : ! xnumrainsv = raindrop number concentration (#/cm^3)
1056 : ! (number in the bin, not number density)
1057 : ! vfallrainsv = fall velocity (cm/s)
1058 : !
1059 0 : precipsum = 0._r8
1060 0 : do i = 1, nr
1061 0 : r = rlo + (i-1)*dr
1062 0 : rrainsv(i) = r
1063 0 : xnumrainsv(i) = exp( -r/2.7e-2_r8 )
1064 :
1065 0 : d = 2._r8*r
1066 0 : if (d <= 0.007_r8) then
1067 0 : vfallstp = 2.88e5_r8 * d**2._r8
1068 0 : else if (d <= 0.025_r8) then
1069 0 : vfallstp = 2.8008e4_r8 * d**1.528_r8
1070 0 : else if (d <= 0.1_r8) then
1071 0 : vfallstp = 4104.9_r8 * d**1.008_r8
1072 0 : else if (d <= 0.25_r8) then
1073 0 : vfallstp = 1812.1_r8 * d**0.638_r8
1074 : else
1075 0 : vfallstp = 1069.8_r8 * d**0.235_r8
1076 : end if
1077 :
1078 0 : vfall = vfallstp * sqrt(1.204e-3_r8/rhoair)
1079 0 : vfallrainsv(i) = vfall
1080 0 : precipsum = precipsum + vfall*(r**3)*xnumrainsv(i)
1081 : end do
1082 0 : precipsum = precipsum*pi*1.333333_r8
1083 :
1084 0 : rnumsum = 0._r8
1085 0 : do i = 1, nr
1086 0 : xnumrainsv(i) = xnumrainsv(i)*(precip/precipsum)
1087 0 : rnumsum = rnumsum + xnumrainsv(i)
1088 : end do
1089 :
1090 : !
1091 : ! compute aerosol concentrations
1092 : ! aaerosv = particle radius (cm)
1093 : ! fnumaerosv = fraction of total number in the bin (--)
1094 : ! fvolaerosv = fraction of total volume in the bin (--)
1095 : !
1096 : anumsum = 0._r8
1097 : avolsum = 0._r8
1098 0 : do i = 1, na
1099 0 : x = xlo + (i-1)*dx
1100 0 : a = exp( x )
1101 0 : aaerosv(i) = a
1102 0 : dum = (x - xg0)/sx
1103 0 : ynumaerosv(i) = exp( -0.5_r8*dum*dum )
1104 0 : yvolaerosv(i) = ynumaerosv(i)*1.3333_r8*pi*a*a*a
1105 0 : anumsum = anumsum + ynumaerosv(i)
1106 0 : avolsum = avolsum + yvolaerosv(i)
1107 : end do
1108 :
1109 0 : do i = 1, na
1110 0 : ynumaerosv(i) = ynumaerosv(i)/anumsum
1111 0 : yvolaerosv(i) = yvolaerosv(i)/avolsum
1112 : end do
1113 :
1114 : !
1115 : ! compute scavenging
1116 : !
1117 : scavsumnum = 0._r8
1118 : scavsumvol = 0._r8
1119 : !
1120 : ! outer loop for rain drop radius
1121 : !
1122 0 : jr_loop: do jr = 1, nr
1123 :
1124 0 : r = rrainsv(jr)
1125 0 : vfall = vfallrainsv(jr)
1126 :
1127 0 : reynolds = r * vfall / airkinvisc
1128 0 : sqrtreynolds = sqrt( reynolds )
1129 :
1130 : !
1131 : ! inner loop for aerosol particle radius
1132 : !
1133 0 : scavsumnumbb = 0._r8
1134 0 : scavsumvolbb = 0._r8
1135 :
1136 0 : ja_loop: do ja = 1, na
1137 :
1138 0 : a = aaerosv(ja)
1139 :
1140 0 : chi = a/r
1141 :
1142 0 : dum = freepath/a
1143 0 : dumfuchs = 1._r8 + 1.246_r8*dum + 0.42_r8*dum*exp(-0.87_r8/dum)
1144 0 : taurelax = 2._r8*rhoaero*a*a*dumfuchs/(9._r8*rhoair*airkinvisc)
1145 :
1146 0 : aeromass = 4._r8*pi*a*a*a*rhoaero/3._r8
1147 0 : aerodiffus = boltz_cgs*temp*taurelax/aeromass
1148 :
1149 0 : schmidt = airkinvisc/aerodiffus
1150 0 : stokes = vfall*taurelax/r
1151 :
1152 : ebrown = 4._r8*(1._r8 + 0.4_r8*sqrtreynolds*(schmidt**0.3333333_r8)) / &
1153 0 : (reynolds*schmidt)
1154 :
1155 : dum = (1._r8 + 2._r8*xmuwaterair*chi) / &
1156 0 : (1._r8 + xmuwaterair/sqrtreynolds)
1157 0 : eintercept = 4._r8*chi*(chi + dum)
1158 :
1159 0 : dum = log( 1._r8 + reynolds )
1160 0 : sstar = (1.2_r8 + dum/12._r8) / (1._r8 + dum)
1161 0 : eimpact = 0._r8
1162 0 : if (stokes > sstar) then
1163 0 : dum = stokes - sstar
1164 0 : eimpact = (dum/(dum+0.6666667_r8)) ** 1.5_r8
1165 : end if
1166 :
1167 0 : etotal = ebrown + eintercept + eimpact
1168 0 : etotal = min( etotal, 1.0_r8 )
1169 :
1170 0 : rainsweepout = xnumrainsv(jr)*4._r8*pi*r*r*vfall
1171 :
1172 0 : scavsumnumbb = scavsumnumbb + rainsweepout*etotal*ynumaerosv(ja)
1173 0 : scavsumvolbb = scavsumvolbb + rainsweepout*etotal*yvolaerosv(ja)
1174 :
1175 : enddo ja_loop
1176 :
1177 0 : scavsumnum = scavsumnum + scavsumnumbb
1178 0 : scavsumvol = scavsumvol + scavsumvolbb
1179 :
1180 : enddo jr_loop
1181 :
1182 0 : scavratenum = scavsumnum*3600._r8
1183 0 : scavratevol = scavsumvol*3600._r8
1184 :
1185 0 : end subroutine calc_1_impact_rate
1186 :
1187 : end subroutine init_bcscavcoef
1188 :
1189 : end module aero_wetdep_cam
|
