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