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 768 : 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 768 : if (masterproc) then
97 2 : open( newunit=unitn, file=trim(nlfile), status='old' )
98 2 : call find_group_name(unitn, 'aero_wetdep_nl', status=ierr)
99 2 : if (ierr == 0) then
100 2 : read(unitn, aero_wetdep_nl, iostat=ierr)
101 2 : if (ierr /= 0) then
102 0 : call endrun(subname // ':: ERROR reading namelist')
103 : end if
104 : end if
105 2 : close(unitn)
106 :
107 : ! ============================
108 : ! Log namelist options
109 : ! ============================
110 2 : write(iulog,*) subname,' namelist settings: '
111 2 : write(iulog,*) ' sol_facti_cloud_borne : ',sol_facti_cloud_borne
112 2 : write(iulog,*) ' sol_factb_interstitial : ',sol_factb_interstitial
113 2 : write(iulog,*) ' sol_factic_interstitial: ',sol_factic_interstitial
114 : end if
115 :
116 : ! ============================
117 : ! Broadcast namelist variables
118 : ! ============================
119 768 : call mpi_bcast(sol_facti_cloud_borne, 1, mpi_real8, masterprocid, mpicom, ierr)
120 768 : if (ierr/=mpi_success) then
121 0 : call endrun(subname//': MPI_BCAST ERROR: sol_facti_cloud_borne')
122 : end if
123 768 : call mpi_bcast(sol_factb_interstitial, 1, mpi_real8, masterprocid, mpicom, ierr)
124 768 : if (ierr/=mpi_success) then
125 0 : call endrun(subname//': MPI_BCAST ERROR: sol_factb_interstitial')
126 : end if
127 768 : call mpi_bcast(sol_factic_interstitial, 1, mpi_real8, masterprocid, mpicom, ierr)
128 768 : if (ierr/=mpi_success) then
129 0 : call endrun(subname//': MPI_BCAST ERROR: sol_factic_interstitial')
130 : end if
131 :
132 768 : call mpi_bcast(nwetdep, 1, mpi_integer, masterprocid, mpicom, ierr)
133 768 : if (ierr/=mpi_success) then
134 0 : call endrun(subname//': MPI_BCAST ERROR: nwetdep')
135 : end if
136 :
137 768 : wetdep_active = .true. !nwetdep>0
138 :
139 768 : if (masterproc) then
140 2 : write(iulog,*) subname,' wetdep_active = ',wetdep_active,' nwetdep = ',nwetdep
141 : endif
142 :
143 768 : call aero_convproc_readnl(nlfile)
144 :
145 768 : end subroutine aero_wetdep_readnl
146 :
147 : !------------------------------------------------------------------------------
148 : !------------------------------------------------------------------------------
149 768 : 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 768 : fracis_idx = pbuf_get_index('FRACIS')
164 768 : rprddp_idx = pbuf_get_index('RPRDDP')
165 768 : rprdsh_idx = pbuf_get_index('RPRDSH')
166 768 : nevapr_shcu_idx = pbuf_get_index('NEVAPR_SHCU')
167 768 : nevapr_dpcu_idx = pbuf_get_index('NEVAPR_DPCU')
168 :
169 768 : 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 768 : call rad_cnst_get_info(0, nmodes=nmodes, nbins=nbins)
176 :
177 768 : if (nmodes>0) then
178 768 : aero_props => modal_aerosol_properties()
179 768 : 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 768 : nele_tot = aero_props%ncnst_tot()
192 :
193 7680 : allocate(aero_cnst_lq(aero_props%nbins(),0:maxval(aero_props%nspecies())), stat=astat)
194 768 : if (astat/=0) then
195 0 : call endrun(subrname//' : not able to allocate aero_cnst_lq array')
196 : end if
197 33024 : aero_cnst_lq(:,:) = .false.
198 :
199 7680 : allocate(aero_cnst_id(aero_props%nbins(),0:maxval(aero_props%nspecies())), stat=astat)
200 768 : if (astat/=0) then
201 0 : call endrun(subrname//' : not able to allocate aero_cnst_id array')
202 : end if
203 33024 : aero_cnst_id(:,:) = -1
204 :
205 768 : wetdep_lq = .false.
206 :
207 4608 : do m = 1, aero_props%nbins()
208 3840 : write(binstr,'(i2.2)') m
209 7680 : call addfld('SOLFACTB'//binstr, (/ 'lev' /), 'A', '1', 'below cld sol fact')
210 :
211 20736 : do l = 0, aero_props%nspecies(m)
212 :
213 16128 : if (l == 0) then ! number
214 3840 : call aero_props%num_names( m, tmpname, tmpname_cw)
215 : else
216 12288 : call aero_props%mmr_names( m,l, tmpname, tmpname_cw)
217 : end if
218 :
219 16128 : call cnst_get_ind(tmpname, id, abort=.false.)
220 16128 : aero_cnst_id(m,l) = id
221 16128 : aero_cnst_lq(m,l) = id > 0
222 16128 : if (id > 0) then
223 16128 : wetdep_lq(id) = .true.
224 : end if
225 :
226 : ! units --
227 16128 : if (l==0) then
228 3840 : unit_basename = ' 1' ! for num
229 : else
230 12288 : unit_basename = 'kg'
231 : endif
232 :
233 16128 : call add_hist_fields(tmpname, unit_basename)
234 16128 : call add_hist_fields(tmpname_cw, unit_basename)
235 :
236 : call addfld( trim(tmpname_cw)//'RSPTD', (/ 'lev' /), 'A', unit_basename//'/kg/s', &
237 36096 : trim(tmpname_cw)//' resuspension tendency')
238 :
239 : end do
240 : end do
241 :
242 2304 : allocate(scavimptblnum(nimptblgrow_mind:nimptblgrow_maxd, aero_props%nbins()), stat=astat)
243 768 : if (astat/=0) then
244 0 : call endrun(subrname//' : not able to allocate scavimptblnum array')
245 : end if
246 2304 : allocate(scavimptblvol(nimptblgrow_mind:nimptblgrow_maxd, aero_props%nbins()), stat=astat)
247 768 : if (astat/=0) then
248 0 : call endrun(subrname//' : not able to allocate scavimptblvol array')
249 : end if
250 768 : scavimptblnum = nan
251 768 : scavimptblvol = nan
252 :
253 768 : call wetdep_init()
254 :
255 4608 : nspec_max = maxval(aero_props%nspecies()) + 2
256 :
257 768 : call init_bcscavcoef()
258 :
259 3072 : if (convproc_do_aer) then
260 768 : call aero_convproc_init(aero_props)
261 : end if
262 :
263 : contains
264 :
265 32256 : 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 32256 : horiz_only, 'A',baseunits//'/m2/s ','Wet deposition flux at surface')
271 : call addfld (trim(name)//'SFSIC', &
272 32256 : horiz_only, 'A',baseunits//'/m2/s ','Wet deposition flux (incloud, convective) at surface')
273 : call addfld (trim(name)//'SFSIS', &
274 32256 : horiz_only, 'A',baseunits//'/m2/s ','Wet deposition flux (incloud, stratiform) at surface')
275 : call addfld (trim(name)//'SFSBC', &
276 32256 : horiz_only, 'A',baseunits//'/m2/s ','Wet deposition flux (belowcloud, convective) at surface')
277 : call addfld (trim(name)//'SFSBS', &
278 32256 : horiz_only, 'A',baseunits//'/m2/s ','Wet deposition flux (belowcloud, stratiform) at surface')
279 :
280 32256 : if (convproc_do_aer) then
281 : call addfld (trim(name)//'SFSEC', &
282 32256 : horiz_only, 'A',unit_basename//'/m2/s','Wet deposition flux (precip evap, convective) at surface')
283 : call addfld (trim(name)//'SFSES', &
284 32256 : horiz_only, 'A',unit_basename//'/m2/s','Wet deposition flux (precip evap, stratiform) at surface')
285 : call addfld (trim(name)//'SFSBD', &
286 32256 : horiz_only, 'A',unit_basename//'/m2/s','Wet deposition flux (belowcloud, deep convective) at surface')
287 : call addfld (trim(name)//'WETC', &
288 64512 : (/ 'lev' /), 'A',unit_basename//'/kg/s ','wet deposition tendency')
289 : call addfld (trim(name)//'CONU', &
290 64512 : (/ 'lev' /), 'A',unit_basename//'/kg ','updraft mixing ratio')
291 : end if
292 :
293 64512 : call addfld (trim(name)//'WET',(/ 'lev' /), 'A',baseunits//'/kg/s ','wet deposition tendency')
294 64512 : 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 64512 : trim(name)//' ic wet deposition')
298 : call addfld (trim(name)//'SIS',(/ 'lev' /), 'A',baseunits//'/kg/s ', &
299 64512 : trim(name)//' is wet deposition')
300 : call addfld (trim(name)//'SBC',(/ 'lev' /), 'A',baseunits//'/kg/s ', &
301 64512 : trim(name)//' bc wet deposition')
302 : call addfld (trim(name)//'SBS',(/ 'lev' /), 'A',baseunits//'/kg/s ', &
303 64512 : trim(name)//' bs wet deposition')
304 :
305 32256 : if ( history_aerosol .or. history_chemistry ) then
306 32256 : call add_default (trim(name)//'SFWET', 1, ' ')
307 : endif
308 32256 : 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 32256 : end subroutine add_hist_fields
321 :
322 : end subroutine aero_wetdep_init
323 :
324 : !------------------------------------------------------------------------------
325 : !------------------------------------------------------------------------------
326 2491776 : 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 24192 : 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 48384 : real(r8) :: diam_wet(state%ncol, pver)
364 : logical :: isprx(pcols,pver) ! true if precipation
365 : real(r8) :: prec(pcols) ! precipitation rate
366 :
367 48384 : real(r8) :: rtscavt(pcols, pver, 0:nspec_max)
368 :
369 : integer :: ncol, lchnk, m, ndx,mm, l
370 : integer :: i,k
371 :
372 24192 : 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 48384 : type(ptr2d_t) :: raer(nele_tot)
377 48384 : type(ptr2d_t) :: qqcw(nele_tot)
378 :
379 : real(r8) :: sflx(pcols)
380 : character(len=aero_name_len) :: aname, cname, name
381 :
382 48384 : 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 48384 : real(r8) :: dcondt_resusp3d(nele_tot,pcols,pver)
389 :
390 : integer, parameter :: nsrflx_mzaer2cnvpr = 2
391 48384 : real(r8) :: qsrflx_mzaer2cnvpr(pcols,nele_tot,nsrflx_mzaer2cnvpr)
392 :
393 24192 : real(r8), pointer :: rprddp(:,:) ! rain production, deep convection
394 24192 : real(r8), pointer :: rprdsh(:,:) ! rain production, shallow convection
395 24192 : real(r8), pointer :: evapcdp(:,:) ! Evaporation rate of deep convective precipitation >=0.
396 24192 : 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 24192 : nullify(aero_state)
412 :
413 24192 : if (.not.wetdep_active) return
414 :
415 1110195072 : dcondt_resusp3d(:,:,:) = 0._r8
416 :
417 24192 : if (nmodes>0) then
418 24192 : aero_state => modal_aerosol_state(state,pbuf)
419 24192 : 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 24192 : lchnk = state%lchnk
432 24192 : ncol = state%ncol
433 :
434 24192 : call physics_ptend_init(ptend, state%psetcols, subrname, lq=wetdep_lq)
435 :
436 24192 : call wetdep_inputs_set( state, pbuf, dep_inputs )
437 :
438 24192 : call pbuf_get_field(pbuf, fracis_idx, fracis)
439 :
440 24192 : call aero_state%get_states( aero_props, raer, qqcw )
441 :
442 17345664 : qsrflx_mzaer2cnvpr(:,:,:) = 0.0_r8
443 24192 : aerdepwetis(:,:) = 0.0_r8
444 24192 : aerdepwetcw(:,:) = 0.0_r8
445 :
446 24192 : 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 314496 : prec(:ncol)=0._r8
460 3169152 : do k=1,pver
461 40884480 : where (prec(:ncol) >= 1.e-7_r8)
462 : isprx(:ncol,k) = .true.
463 : elsewhere
464 3144960 : 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 40908672 : *state%pdel(:ncol,k)/gravit
468 : end do
469 :
470 24192 : f_act_conv = 0._r8
471 24192 : scavcoefnv = nan
472 24192 : qqcw_sav = nan
473 :
474 24192 : if (convproc_do_aer) then
475 :
476 24192 : call t_startf('aero_convproc')
477 : call aero_convproc_intr( aero_props, aero_state, state, ptend, pbuf, dt, &
478 24192 : nsrflx_mzaer2cnvpr, qsrflx_mzaer2cnvpr, aerdepwetis, dcondt_resusp3d )
479 :
480 24192 : if (convproc_do_evaprain_atonce) then
481 :
482 145152 : do m = 1,aero_props%nbins()
483 653184 : do l = 0,aero_props%nspecies(m)
484 508032 : mm = aero_props%indexer(m,l)
485 :
486 508032 : if (l == 0) then ! number
487 120960 : call aero_props%num_names(m, aname, cname)
488 : else
489 387072 : call aero_props%mmr_names(m,l, aname, cname)
490 : end if
491 :
492 859082112 : call outfld( trim(cname)//'RSPTD', dcondt_resusp3d(mm,:ncol,:), ncol, lchnk )
493 :
494 66673152 : do k = 1,pver
495 859082112 : do i = 1,ncol
496 858574080 : 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 24192 : call t_stopf('aero_convproc')
504 :
505 : end if
506 :
507 145152 : bins_loop: do m = 1,aero_props%nbins()
508 :
509 387072 : phase_loop: do lphase = strt_loop,end_loop, stride_loop ! loop over interstitial (1) and cloud-borne (2) forms
510 :
511 241920 : cldbrn = lphase==2
512 :
513 241920 : sol_factb = nan
514 241920 : sol_facti = nan
515 241920 : sol_factic = nan
516 :
517 241920 : if (lphase == 1) then ! interstial aerosol
518 :
519 120960 : sol_facti = 0.0_r8 ! strat in-cloud scav totally OFF for institial
520 :
521 267442560 : sol_factic = sol_factic_interstitial
522 :
523 : else ! cloud-borne aerosol (borne by stratiform cloud drops)
524 :
525 120960 : sol_factb = 0.0_r8 ! all below-cloud scav OFF (anything cloud-borne is located "in-cloud")
526 267442560 : sol_facti = sol_facti_cloud_borne ! strat in-cloud scav cloud-borne tuning factor
527 120960 : sol_factic = 0.0_r8 ! conv in-cloud scav OFF (having this on would mean
528 : ! that conv precip collects strat droplets)
529 120960 : f_act_conv = 0.0_r8 ! conv in-cloud scav OFF (having this on would mean
530 :
531 : end if
532 241920 : 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 120960 : sol_factic = 0.0_r8
540 : endif
541 :
542 409086720 : diam_wet = aero_state%wet_diameter(m,ncol,pver)
543 :
544 241920 : scavcoefnv = 0.0_r8
545 :
546 241920 : if (lphase == 1) then ! interstial aerosol
547 120960 : call get_bcscavcoefs( m, ncol, isprx, diam_wet, scavcoefnv(:,:,1), scavcoefnv(:,:,2) )
548 :
549 120960 : if ( sol_factb_interstitial /= NOTSET ) then
550 204543360 : 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 120960 : write(binstr,'(i2.2)') m
556 120960 : call outfld('SOLFACTB'//binstr,sol_factb, pcols, lchnk)
557 :
558 : end if
559 :
560 1378944 : elem_loop: do l = 0,aero_props%nspecies(m)
561 :
562 1016064 : ndx = aero_cnst_id(m,l)
563 1016064 : if (ndx<1) cycle elem_loop
564 :
565 1016064 : if (.not. cldbrn .and. ndx>0) then
566 508032 : insolfr_ptr => fracis(:,:,ndx)
567 : else
568 508032 : insolfr_ptr => fracis_nadv
569 : endif
570 :
571 1016064 : mm = aero_props%indexer(m,l)
572 :
573 1016064 : if (l == 0) then ! number
574 241920 : call aero_props%num_names( m, aname, cname)
575 : else
576 774144 : call aero_props%mmr_names( m,l, aname, cname)
577 : end if
578 :
579 1016064 : if (cldbrn) then
580 859082112 : q_tmp(1:ncol,:) = qqcw(mm)%fld(1:ncol,:)
581 508032 : jnv = 0
582 508032 : if (convproc_do_aer) then
583 859082112 : qqcw_sav(:ncol,:,l) = q_tmp(1:ncol,:)
584 : endif
585 508032 : name = cname
586 508032 : qqcw_in = nan
587 508032 : f_act_conv = nan
588 : else ! interstial aerosol
589 859082112 : q_tmp(1:ncol,:) = raer(mm)%fld(1:ncol,:) + ptend%q(1:ncol,:,ndx)*dt
590 508032 : if (l==0) then
591 : jnv = 1
592 : else
593 387072 : jnv = 2
594 : end if
595 508032 : if(convproc_do_aer) then
596 : !Feed in the saved cloudborne mixing ratios from phase 2
597 859082112 : qqcw_in(:ncol,:) = qqcw_sav(:ncol,:,l)
598 : else
599 0 : qqcw_in(:ncol,:) = qqcw(mm)%fld(:ncol,:)
600 : end if
601 :
602 859082112 : f_act_conv(:ncol,:) = aero_state%convcld_actfrac( m, l, ncol, pver)
603 508032 : name = aname
604 : end if
605 :
606 1718164224 : 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 1016064 : bergso_in=dep_inputs%bergso )
622 :
623 1016064 : if(convproc_do_aer) then
624 1016064 : if(cldbrn) then
625 : ! save resuspension of cloudborne species
626 859082112 : 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 859082112 : 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 859082112 : dqdt_tmp(1:ncol,:) = dqdt_tmp(1:ncol,:) + rtscavt(1:ncol,:,l)
633 : end if
634 : endif
635 :
636 1016064 : if (cldbrn) then
637 66552192 : do k = 1,pver
638 859082112 : do i = 1,ncol
639 858574080 : 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 859082112 : qqcw(mm)%fld(1:ncol,:) = qqcw(mm)%fld(1:ncol,:) + dqdt_tmp(1:ncol,:) * dt
646 :
647 : else
648 859082112 : ptend%q(1:ncol,:,ndx) = ptend%q(1:ncol,:,ndx) + dqdt_tmp(1:ncol,:)
649 : end if
650 :
651 1016064 : call outfld( trim(name)//'WET', dqdt_tmp(:,:), pcols, lchnk)
652 1016064 : call outfld( trim(name)//'SIC', icscavt, pcols, lchnk)
653 1016064 : call outfld( trim(name)//'SIS', isscavt, pcols, lchnk)
654 1016064 : call outfld( trim(name)//'SBC', bcscavt, pcols, lchnk)
655 1016064 : call outfld( trim(name)//'SBS', bsscavt, pcols, lchnk)
656 :
657 1016064 : call outfld( trim(name)//'INS', insolfr_ptr, pcols, lchnk)
658 :
659 1016064 : sflx(:)=0._r8
660 133104384 : do k=1,pver
661 1718164224 : do i=1,ncol
662 1717148160 : sflx(i)=sflx(i)+dqdt_tmp(i,k)*state%pdel(i,k)/gravit
663 : enddo
664 : enddo
665 1016064 : if (cldbrn) then
666 508032 : call outfld( trim(name)//'SFWET', sflx, pcols, lchnk)
667 6604416 : if (ndx>0) aerdepwetcw(:ncol,ndx) = sflx(:ncol)
668 : else
669 508032 : if (.not.convproc_do_aer) call outfld( trim(name)//'SFWET', sflx, pcols, lchnk)
670 6604416 : if (ndx>0) aerdepwetis(:ncol,ndx) = aerdepwetis(:ncol,ndx) + sflx(:ncol)
671 : end if
672 :
673 1016064 : sflx(:)=0._r8
674 133104384 : do k=1,pver
675 1718164224 : do i=1,ncol
676 1717148160 : sflx(i)=sflx(i)+icscavt(i,k)*state%pdel(i,k)/gravit
677 : enddo
678 : enddo
679 1016064 : if (cldbrn) then
680 508032 : call outfld( trim(name)//'SFSIC', sflx, pcols, lchnk)
681 : else
682 508032 : if (.not.convproc_do_aer) call outfld( trim(name)//'SFSIC', sflx, pcols, lchnk)
683 508032 : if (convproc_do_aer) sflxic = sflx
684 : end if
685 :
686 1016064 : sflx(:)=0._r8
687 133104384 : do k=1,pver
688 1718164224 : do i=1,ncol
689 1717148160 : sflx(i)=sflx(i)+isscavt(i,k)*state%pdel(i,k)/gravit
690 : enddo
691 : enddo
692 1016064 : call outfld( trim(name)//'SFSIS', sflx, pcols, lchnk)
693 :
694 1016064 : sflx(:)=0._r8
695 133104384 : do k=1,pver
696 1718164224 : do i=1,ncol
697 1717148160 : sflx(i)=sflx(i)+bcscavt(i,k)*state%pdel(i,k)/gravit
698 : enddo
699 : enddo
700 1016064 : call outfld( trim(name)//'SFSBC', sflx, pcols, lchnk)
701 1016064 : if (convproc_do_aer) sflxbc = sflx
702 :
703 1016064 : sflx(:)=0._r8
704 133104384 : do k=1,pver
705 1718164224 : do i=1,ncol
706 1717148160 : sflx(i)=sflx(i)+bsscavt(i,k)*state%pdel(i,k)/gravit
707 : enddo
708 : enddo
709 1016064 : call outfld( trim(name)//'SFSBS', sflx, pcols, lchnk)
710 :
711 1257984 : if(convproc_do_aer) then
712 :
713 1016064 : sflx(:)=0.0_r8
714 133104384 : do k=1,pver
715 1718164224 : do i=1,ncol
716 1717148160 : sflx(i)=sflx(i)+rsscavt(i,k)*state%pdel(i,k)/gravit
717 : enddo
718 : enddo
719 1016064 : call outfld( trim(name)//'SFSES', sflx, pcols, lchnk)
720 :
721 1016064 : sflx(:)=0._r8
722 133104384 : do k=1,pver
723 1718164224 : do i=1,ncol
724 1717148160 : sflx(i)=sflx(i)+rcscavt(i,k)*state%pdel(i,k)/gravit
725 : enddo
726 : enddo
727 1016064 : if (.not.convproc_do_aer) call outfld( trim(name)//'SFSEC', sflx, pcols, lchnk)
728 1016064 : sflxec = sflx
729 :
730 1016064 : 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 508032 : if ( deepconv_wetdep_history) then
739 :
740 508032 : call pbuf_get_field(pbuf, rprddp_idx, rprddp )
741 508032 : call pbuf_get_field(pbuf, rprdsh_idx, rprdsh )
742 508032 : call pbuf_get_field(pbuf, nevapr_dpcu_idx, evapcdp )
743 508032 : call pbuf_get_field(pbuf, nevapr_shcu_idx, evapcsh )
744 :
745 508032 : rprddpsum(:) = 0.0_r8
746 508032 : rprdshsum(:) = 0.0_r8
747 508032 : evapcdpsum(:) = 0.0_r8
748 508032 : evapcshsum(:) = 0.0_r8
749 :
750 66552192 : do k = 1, pver
751 858574080 : rprddpsum(:ncol) = rprddpsum(:ncol) + rprddp(:ncol,k)*state%pdel(:ncol,k)/gravit
752 858574080 : rprdshsum(:ncol) = rprdshsum(:ncol) + rprdsh(:ncol,k)*state%pdel(:ncol,k)/gravit
753 858574080 : evapcdpsum(:ncol) = evapcdpsum(:ncol) + evapcdp(:ncol,k)*state%pdel(:ncol,k)/gravit
754 859082112 : evapcshsum(:ncol) = evapcshsum(:ncol) + evapcsh(:ncol,k)*state%pdel(:ncol,k)/gravit
755 : end do
756 :
757 6604416 : do i = 1, ncol
758 6096384 : rprddpsum(i) = max( rprddpsum(i), 1.0e-35_r8 )
759 6096384 : rprdshsum(i) = max( rprdshsum(i), 1.0e-35_r8 )
760 6096384 : evapcdpsum(i) = max( evapcdpsum(i), 0.1e-35_r8 )
761 6096384 : 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 6096384 : tmpa = rprddpsum(i) / (rprddpsum(i) + rprdshsum(i))
765 6096384 : tmpa = max( 0.0_r8, min( 1.0_r8, tmpa ) )
766 6096384 : sflxicdp(i) = sflxic(i)*tmpa
767 6096384 : sflxbcdp(i) = sflxbc(i)*tmpa
768 :
769 : ! assume that resuspension is proportional to (wet removal)*[(precip evap)/(precip production)]
770 6096384 : tmp_resudp = tmpa * min( (evapcdpsum(i)/rprddpsum(i)), 1.0_r8 )
771 6096384 : tmp_resush = (1.0_r8 - tmpa) * min( (evapcshsum(i)/rprdshsum(i)), 1.0_r8 )
772 6096384 : tmpb = max( tmp_resudp, 1.0e-35_r8 ) / max( (tmp_resudp+tmp_resush), 1.0e-35_r8 )
773 6096384 : tmpb = max( 0.0_r8, min( 1.0_r8, tmpb ) )
774 6604416 : sflxecdp(i) = sflxec(i)*tmpb
775 : end do
776 508032 : 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 6604416 : qsrflx_mzaer2cnvpr(1:ncol,mm,1) = sflxec( 1:ncol)
787 6604416 : 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 24192 : if (associated(aero_state)) then
798 48384 : 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 72576 : if (.not. aerodep_flx_prescribed()) then
805 24192 : call aero_deposition_cam_setwet(aerdepwetis, aerdepwetcw, cam_out)
806 : endif
807 :
808 : contains
809 :
810 : ! below cloud impaction scavenging coefs
811 120960 : 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 15845760 : do k = 1, pver
822 204543360 : do i = 1, ncol
823 :
824 : ! do only if no precip
825 204422400 : 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 11252425 : dumdgratio = diam_wet(i,k)/aero_props%scav_diam(m)
830 11252425 : if ((dumdgratio >= 0.99_r8) .and. (dumdgratio <= 1.01_r8)) then
831 108858 : scavimpvol = scavimptblvol(0,m)
832 108858 : scavimpnum = scavimptblnum(0,m)
833 : else
834 11143567 : xgrow = log( dumdgratio ) / dlndg_nimptblgrow
835 11143567 : jgrow = int( xgrow )
836 11143567 : if (xgrow < 0._r8) jgrow = jgrow - 1
837 11143567 : if (jgrow < nimptblgrow_mind) then
838 : jgrow = nimptblgrow_mind
839 : xgrow = jgrow
840 : else
841 11143567 : jgrow = min( jgrow, nimptblgrow_maxd-1 )
842 : end if
843 :
844 11143567 : dumfhi = xgrow - jgrow
845 11143567 : dumflo = 1._r8 - dumfhi
846 :
847 0 : scavimpvol = dumflo*scavimptblvol(jgrow,m) + &
848 11143567 : dumfhi*scavimptblvol(jgrow+1,m)
849 0 : scavimpnum = dumflo*scavimptblnum(jgrow,m) + &
850 11143567 : dumfhi*scavimptblnum(jgrow+1,m)
851 :
852 : end if
853 :
854 : ! impaction scavenging removal amount for volume
855 11252425 : scavcoefvol(i,k) = exp( scavimpvol )
856 : ! impaction scavenging removal amount to number
857 11252425 : scavcoefnum(i,k) = exp( scavimpnum )
858 :
859 : else
860 177445175 : scavcoefvol(i,k) = 0._r8
861 177445175 : scavcoefnum(i,k) = 0._r8
862 : end if
863 :
864 : end do
865 : end do
866 :
867 24192 : end subroutine get_bcscavcoefs
868 :
869 : end subroutine aero_wetdep_tend
870 :
871 : !------------------------------------------------------------------------------
872 : !------------------------------------------------------------------------------
873 768 : 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 768 : lunerr = iulog
905 768 : dlndg_nimptblgrow = log( 1.25_r8 )
906 :
907 : ! bin model: main loop over aerosol bins
908 :
909 4608 : 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 3840 : call aero_props%get(m,1,density=rhodryaero)
914 :
915 3840 : dg0_base = aero_props%scav_diam(m)
916 :
917 3840 : logsig = aero_props%alogsig(m)
918 :
919 81408 : growloop: do jgrow = nimptblgrow_mind, nimptblgrow_maxd
920 :
921 76800 : wetdiaratio = exp( jgrow*dlndg_nimptblgrow )
922 76800 : dg0 = dg0_base*wetdiaratio
923 :
924 : wetvolratio = exp( jgrow*dlndg_nimptblgrow*3._r8 )
925 76800 : rhowetaero = 1.0_r8 + (rhodryaero-1.0_r8)/wetvolratio
926 76800 : rhowetaero = min( rhowetaero, rhodryaero )
927 :
928 : !
929 : ! compute impaction scavenging rates at 1 temp-press pair and save
930 : !
931 76800 : nnfit = 0
932 :
933 76800 : temp = 273.16_r8
934 76800 : press = 0.75e6_r8 ! dynes/cm2
935 76800 : rhowetaero = rhodryaero
936 :
937 76800 : dg0_cgs = dg0*1.0e2_r8 ! m to cm
938 :
939 76800 : 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 76800 : scavratenum, scavratevol, lunerr )
944 :
945 76800 : 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 76800 : xxfitnum(1,nnfit) = 1._r8
953 76800 : yyfitnum(nnfit) = log( scavratenum )
954 :
955 76800 : xxfitvol(1,nnfit) = 1._r8
956 76800 : yyfitvol(nnfit) = log( scavratevol )
957 :
958 : !depends on both bins and different species
959 76800 : scavimptblnum(jgrow,m) = yyfitnum(1)
960 157440 : scavimptblvol(jgrow,m) = yyfitvol(1)
961 :
962 : enddo growloop
963 : enddo modeloop
964 :
965 : contains
966 :
967 : !===============================================================================
968 76800 : 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 76800 : rlo = .005_r8
1018 76800 : rhi = .250_r8
1019 76800 : dr = 0.005_r8
1020 76800 : 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 76800 : precipmmhr = 1.0_r8
1029 76800 : precip = precipmmhr/36000._r8
1030 :
1031 76800 : ag0 = dg0/2._r8
1032 76800 : sx = logsig
1033 76800 : xg0 = log( ag0 )
1034 76800 : xg3 = xg0 + 3._r8*sx*sx
1035 :
1036 76800 : xlo = xg3 - 4._r8*sx
1037 76800 : xhi = xg3 + 4._r8*sx
1038 76800 : dx = 0.2_r8*sx
1039 :
1040 76800 : dx = max( 0.2_r8*sx, 0.01_r8 )
1041 76800 : xlo = xg3 - max( 4._r8*sx, 2._r8*dx )
1042 76800 : xhi = xg3 + max( 4._r8*sx, 2._r8*dx )
1043 :
1044 76800 : na = 1 + nint( (xhi-xlo)/dx )
1045 76800 : 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 76800 : cair = press/(rgas*temp)
1054 : ! air mass density
1055 76800 : rhoair = 28.966_r8*cair
1056 : ! molecular freepath
1057 76800 : freepath = 2.8052e-10_r8/cair
1058 : ! air dynamic viscosity
1059 : airdynvisc = 1.8325e-4_r8 * (416.16_r8/(temp+120._r8)) * &
1060 76800 : ((temp/296.16_r8)**1.5_r8)
1061 : ! air kinemaic viscosity
1062 76800 : airkinvisc = airdynvisc/rhoair
1063 : ! ratio of water viscosity to air viscosity (from Slinn)
1064 76800 : 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 76800 : precipsum = 0._r8
1074 3916800 : do i = 1, nr
1075 3840000 : r = rlo + (i-1)*dr
1076 3840000 : rrainsv(i) = r
1077 3840000 : xnumrainsv(i) = exp( -r/2.7e-2_r8 )
1078 :
1079 3840000 : d = 2._r8*r
1080 3840000 : if (d <= 0.007_r8) then
1081 0 : vfallstp = 2.88e5_r8 * d**2._r8
1082 3840000 : else if (d <= 0.025_r8) then
1083 153600 : vfallstp = 2.8008e4_r8 * d**1.528_r8
1084 3686400 : else if (d <= 0.1_r8) then
1085 614400 : vfallstp = 4104.9_r8 * d**1.008_r8
1086 3072000 : else if (d <= 0.25_r8) then
1087 1152000 : vfallstp = 1812.1_r8 * d**0.638_r8
1088 : else
1089 1920000 : vfallstp = 1069.8_r8 * d**0.235_r8
1090 : end if
1091 :
1092 3840000 : vfall = vfallstp * sqrt(1.204e-3_r8/rhoair)
1093 3840000 : vfallrainsv(i) = vfall
1094 3916800 : precipsum = precipsum + vfall*(r**3)*xnumrainsv(i)
1095 : end do
1096 76800 : precipsum = precipsum*pi*1.333333_r8
1097 :
1098 76800 : rnumsum = 0._r8
1099 3916800 : do i = 1, nr
1100 3840000 : xnumrainsv(i) = xnumrainsv(i)*(precip/precipsum)
1101 3916800 : 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 3225600 : do i = 1, na
1113 3148800 : x = xlo + (i-1)*dx
1114 3148800 : a = exp( x )
1115 3148800 : aaerosv(i) = a
1116 3148800 : dum = (x - xg0)/sx
1117 3148800 : ynumaerosv(i) = exp( -0.5_r8*dum*dum )
1118 3148800 : yvolaerosv(i) = ynumaerosv(i)*1.3333_r8*pi*a*a*a
1119 3148800 : anumsum = anumsum + ynumaerosv(i)
1120 3225600 : avolsum = avolsum + yvolaerosv(i)
1121 : end do
1122 :
1123 3225600 : do i = 1, na
1124 3148800 : ynumaerosv(i) = ynumaerosv(i)/anumsum
1125 3225600 : 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 3916800 : jr_loop: do jr = 1, nr
1137 :
1138 3840000 : r = rrainsv(jr)
1139 3840000 : vfall = vfallrainsv(jr)
1140 :
1141 3840000 : reynolds = r * vfall / airkinvisc
1142 3840000 : sqrtreynolds = sqrt( reynolds )
1143 :
1144 : !
1145 : ! inner loop for aerosol particle radius
1146 : !
1147 3840000 : scavsumnumbb = 0._r8
1148 3840000 : scavsumvolbb = 0._r8
1149 :
1150 161280000 : ja_loop: do ja = 1, na
1151 :
1152 157440000 : a = aaerosv(ja)
1153 :
1154 157440000 : chi = a/r
1155 :
1156 157440000 : dum = freepath/a
1157 157440000 : dumfuchs = 1._r8 + 1.246_r8*dum + 0.42_r8*dum*exp(-0.87_r8/dum)
1158 157440000 : taurelax = 2._r8*rhoaero*a*a*dumfuchs/(9._r8*rhoair*airkinvisc)
1159 :
1160 157440000 : aeromass = 4._r8*pi*a*a*a*rhoaero/3._r8
1161 157440000 : aerodiffus = boltz_cgs*temp*taurelax/aeromass
1162 :
1163 157440000 : schmidt = airkinvisc/aerodiffus
1164 157440000 : stokes = vfall*taurelax/r
1165 :
1166 : ebrown = 4._r8*(1._r8 + 0.4_r8*sqrtreynolds*(schmidt**0.3333333_r8)) / &
1167 157440000 : (reynolds*schmidt)
1168 :
1169 : dum = (1._r8 + 2._r8*xmuwaterair*chi) / &
1170 157440000 : (1._r8 + xmuwaterair/sqrtreynolds)
1171 157440000 : eintercept = 4._r8*chi*(chi + dum)
1172 :
1173 157440000 : dum = log( 1._r8 + reynolds )
1174 157440000 : sstar = (1.2_r8 + dum/12._r8) / (1._r8 + dum)
1175 157440000 : eimpact = 0._r8
1176 157440000 : if (stokes > sstar) then
1177 44023296 : dum = stokes - sstar
1178 44023296 : eimpact = (dum/(dum+0.6666667_r8)) ** 1.5_r8
1179 : end if
1180 :
1181 157440000 : etotal = ebrown + eintercept + eimpact
1182 157440000 : etotal = min( etotal, 1.0_r8 )
1183 :
1184 157440000 : rainsweepout = xnumrainsv(jr)*4._r8*pi*r*r*vfall
1185 :
1186 157440000 : scavsumnumbb = scavsumnumbb + rainsweepout*etotal*ynumaerosv(ja)
1187 161280000 : scavsumvolbb = scavsumvolbb + rainsweepout*etotal*yvolaerosv(ja)
1188 :
1189 : enddo ja_loop
1190 :
1191 3840000 : scavsumnum = scavsumnum + scavsumnumbb
1192 3916800 : scavsumvol = scavsumvol + scavsumvolbb
1193 :
1194 : enddo jr_loop
1195 :
1196 76800 : scavratenum = scavsumnum*3600._r8
1197 76800 : scavratevol = scavsumvol*3600._r8
1198 :
1199 76800 : end subroutine calc_1_impact_rate
1200 :
1201 : end subroutine init_bcscavcoef
1202 :
1203 : end module aero_wetdep_cam
|
