Line data Source code
1 : module modal_aerosol_properties_mod
2 : use shr_kind_mod, only: r8 => shr_kind_r8
3 : use physconst, only: pi
4 : use aerosol_properties_mod, only: aerosol_properties, aero_name_len
5 : use rad_constituents, only: rad_cnst_get_info, rad_cnst_get_mode_props, rad_cnst_get_aer_props
6 :
7 : implicit none
8 :
9 : private
10 :
11 : public :: modal_aerosol_properties
12 :
13 : type, extends(aerosol_properties) :: modal_aerosol_properties
14 : private
15 : real(r8), allocatable :: exp45logsig_(:)
16 : real(r8), allocatable :: voltonumblo_(:)
17 : real(r8), allocatable :: voltonumbhi_(:)
18 : integer, allocatable :: sulfate_mode_ndxs_(:)
19 : integer, allocatable :: dust_mode_ndxs_(:)
20 : integer, allocatable :: ssalt_mode_ndxs_(:)
21 : integer, allocatable :: ammon_mode_ndxs_(:)
22 : integer, allocatable :: nitrate_mode_ndxs_(:)
23 : integer, allocatable :: msa_mode_ndxs_(:)
24 : integer, allocatable :: bcarbon_mode_ndxs_(:,:)
25 : integer, allocatable :: porganic_mode_ndxs_(:,:)
26 : integer, allocatable :: sorganic_mode_ndxs_(:,:)
27 : integer :: num_soa_ = 0
28 : integer :: num_poa_ = 0
29 : integer :: num_bc_ = 0
30 : contains
31 : procedure :: number_transported
32 : procedure :: get
33 : procedure :: amcube
34 : procedure :: actfracs
35 : procedure :: num_names
36 : procedure :: mmr_names
37 : procedure :: amb_num_name
38 : procedure :: amb_mmr_name
39 : procedure :: species_type
40 : procedure :: icenuc_updates_num
41 : procedure :: icenuc_updates_mmr
42 : procedure :: apply_number_limits
43 : procedure :: hetfrz_species
44 : procedure :: optics_params
45 : procedure :: nbins_rlist
46 : procedure :: nspecies_per_bin_rlist
47 : procedure :: alogsig_rlist
48 : procedure :: soluble
49 : procedure :: min_mass_mean_rad
50 : procedure :: bin_name
51 : procedure :: scav_diam
52 : procedure :: resuspension_resize
53 : procedure :: rebin_bulk_fluxes
54 : procedure :: hydrophilic
55 :
56 : final :: destructor
57 : end type modal_aerosol_properties
58 :
59 : interface modal_aerosol_properties
60 : procedure :: constructor
61 : end interface modal_aerosol_properties
62 :
63 : logical, parameter :: debug = .false.
64 :
65 : contains
66 :
67 : !------------------------------------------------------------------------------
68 : !------------------------------------------------------------------------------
69 0 : function constructor() result(newobj)
70 :
71 : type(modal_aerosol_properties), pointer :: newobj
72 :
73 : integer :: l, m, nmodes, ncnst_tot, mm
74 : real(r8) :: dgnumlo
75 : real(r8) :: dgnumhi
76 0 : integer,allocatable :: nspecies(:)
77 0 : real(r8),allocatable :: sigmag(:)
78 0 : real(r8),allocatable :: alogsig(:)
79 0 : real(r8),allocatable :: f1(:)
80 0 : real(r8),allocatable :: f2(:)
81 : integer :: ierr
82 :
83 : character(len=aero_name_len) :: spectype
84 :
85 : integer :: npoa, nsoa, nbc
86 :
87 0 : allocate(newobj,stat=ierr)
88 0 : if( ierr /= 0 ) then
89 0 : nullify(newobj)
90 0 : return
91 : end if
92 :
93 0 : call rad_cnst_get_info(0, nmodes=nmodes)
94 :
95 0 : allocate(nspecies(nmodes),stat=ierr)
96 0 : if( ierr /= 0 ) then
97 0 : nullify(newobj)
98 : return
99 : end if
100 0 : allocate(alogsig(nmodes),stat=ierr)
101 0 : if( ierr /= 0 ) then
102 0 : nullify(newobj)
103 : return
104 : end if
105 0 : allocate( f1(nmodes),stat=ierr )
106 0 : if( ierr /= 0 ) then
107 0 : nullify(newobj)
108 : return
109 : end if
110 0 : allocate( f2(nmodes),stat=ierr )
111 0 : if( ierr /= 0 ) then
112 0 : nullify(newobj)
113 : return
114 : end if
115 :
116 0 : allocate(sigmag(nmodes),stat=ierr)
117 0 : if( ierr /= 0 ) then
118 0 : nullify(newobj)
119 : return
120 : end if
121 0 : allocate(newobj%exp45logsig_(nmodes),stat=ierr)
122 0 : if( ierr /= 0 ) then
123 0 : nullify(newobj)
124 : return
125 : end if
126 0 : allocate(newobj%voltonumblo_(nmodes),stat=ierr)
127 0 : if( ierr /= 0 ) then
128 0 : nullify(newobj)
129 : return
130 : end if
131 0 : allocate(newobj%voltonumbhi_(nmodes),stat=ierr)
132 0 : if( ierr /= 0 ) then
133 0 : nullify(newobj)
134 : return
135 : end if
136 :
137 0 : ncnst_tot = 0
138 :
139 0 : do m = 1, nmodes
140 0 : call rad_cnst_get_info(0, m, nspec=nspecies(m))
141 :
142 0 : ncnst_tot = ncnst_tot + nspecies(m) + 1
143 :
144 : call rad_cnst_get_mode_props(0, m, sigmag=sigmag(m), &
145 0 : dgnumhi=dgnumhi, dgnumlo=dgnumlo )
146 :
147 0 : alogsig(m) = log(sigmag(m))
148 :
149 0 : newobj%exp45logsig_(m) = exp(4.5_r8*alogsig(m)*alogsig(m))
150 :
151 0 : f1(m) = 0.5_r8*exp(2.5_r8*alogsig(m)*alogsig(m))
152 0 : f2(m) = 1._r8 + 0.25_r8*alogsig(m)
153 :
154 0 : newobj%voltonumblo_(m) = 1._r8 / ( (pi/6._r8)* &
155 0 : (dgnumlo**3._r8)*exp(4.5_r8*alogsig(m)**2._r8) )
156 0 : newobj%voltonumbhi_(m) = 1._r8 / ( (pi/6._r8)* &
157 0 : (dgnumhi**3._r8)*exp(4.5_r8*alogsig(m)**2._r8) )
158 :
159 : end do
160 :
161 0 : call newobj%initialize(nmodes,ncnst_tot,nspecies,nspecies,alogsig,f1,f2,ierr)
162 :
163 0 : npoa = 0
164 0 : nsoa = 0
165 0 : nbc = 0
166 :
167 0 : m = 1
168 0 : do l = 1,newobj%nspecies(m)
169 0 : mm = newobj%indexer(m,l)
170 0 : call newobj%species_type(m, l, spectype)
171 0 : select case ( trim(spectype) )
172 : case('p-organic')
173 0 : npoa = npoa + 1
174 : case('s-organic')
175 0 : nsoa = nsoa + 1
176 : case('black-c')
177 0 : nbc = nbc + 1
178 : end select
179 : end do
180 :
181 0 : newobj%num_soa_ = nsoa
182 0 : newobj%num_poa_ = npoa
183 0 : newobj%num_bc_ = nbc
184 :
185 0 : allocate(newobj%sulfate_mode_ndxs_(newobj%nbins()),stat=ierr)
186 0 : if( ierr /= 0 ) then
187 0 : nullify(newobj)
188 : return
189 : end if
190 0 : allocate(newobj%dust_mode_ndxs_(newobj%nbins()),stat=ierr)
191 0 : if( ierr /= 0 ) then
192 0 : nullify(newobj)
193 : return
194 : end if
195 0 : allocate(newobj%ssalt_mode_ndxs_(newobj%nbins()),stat=ierr)
196 0 : if( ierr /= 0 ) then
197 0 : nullify(newobj)
198 : return
199 : end if
200 0 : allocate(newobj%ammon_mode_ndxs_(newobj%nbins()),stat=ierr)
201 0 : if( ierr /= 0 ) then
202 0 : nullify(newobj)
203 : return
204 : end if
205 0 : allocate(newobj%nitrate_mode_ndxs_(newobj%nbins()),stat=ierr)
206 0 : if( ierr /= 0 ) then
207 0 : nullify(newobj)
208 : return
209 : end if
210 0 : allocate(newobj%msa_mode_ndxs_(newobj%nbins()),stat=ierr)
211 0 : if( ierr /= 0 ) then
212 0 : nullify(newobj)
213 : return
214 : end if
215 :
216 0 : newobj%sulfate_mode_ndxs_ = 0
217 0 : newobj%dust_mode_ndxs_ = 0
218 0 : newobj%ssalt_mode_ndxs_ = 0
219 0 : newobj%ammon_mode_ndxs_ = 0
220 0 : newobj%nitrate_mode_ndxs_ = 0
221 0 : newobj%msa_mode_ndxs_ = 0
222 :
223 0 : allocate(newobj%porganic_mode_ndxs_(newobj%nbins(),npoa),stat=ierr)
224 0 : if( ierr /= 0 ) then
225 0 : nullify(newobj)
226 : return
227 : end if
228 0 : allocate(newobj%sorganic_mode_ndxs_(newobj%nbins(),nsoa),stat=ierr)
229 0 : if( ierr /= 0 ) then
230 0 : nullify(newobj)
231 : return
232 : end if
233 0 : allocate(newobj%bcarbon_mode_ndxs_(newobj%nbins(),nbc),stat=ierr)
234 0 : if( ierr /= 0 ) then
235 0 : nullify(newobj)
236 : return
237 : end if
238 :
239 0 : newobj%porganic_mode_ndxs_ = 0._r8
240 0 : newobj%sorganic_mode_ndxs_ = 0._r8
241 0 : newobj%bcarbon_mode_ndxs_ = 0._r8
242 :
243 0 : do m = 1,newobj%nbins()
244 0 : npoa = 0
245 0 : nsoa = 0
246 0 : nbc = 0
247 :
248 0 : do l = 1,newobj%nspecies(m)
249 0 : mm = newobj%indexer(m,l)
250 0 : call newobj%species_type(m, l, spectype)
251 :
252 0 : select case ( trim(spectype) )
253 : case('sulfate')
254 0 : newobj%sulfate_mode_ndxs_(m) = mm
255 : case('dust')
256 0 : newobj%dust_mode_ndxs_(m) = mm
257 : case('nitrate')
258 0 : newobj%nitrate_mode_ndxs_(m) = mm
259 : case('ammonium')
260 0 : newobj%ammon_mode_ndxs_(m) = mm
261 : case('seasalt')
262 0 : newobj%ssalt_mode_ndxs_(m) = mm
263 : case('msa')
264 0 : newobj%msa_mode_ndxs_(m) = mm
265 : case('p-organic')
266 0 : npoa = npoa + 1
267 0 : newobj%porganic_mode_ndxs_(m,npoa) = mm
268 : case('s-organic')
269 0 : nsoa = nsoa + 1
270 0 : newobj%sorganic_mode_ndxs_(m,nsoa) = mm
271 : case('black-c')
272 0 : nbc = nbc + 1
273 0 : newobj%bcarbon_mode_ndxs_(m,nbc) = mm
274 : end select
275 :
276 : end do
277 : end do
278 :
279 0 : if( ierr /= 0 ) then
280 0 : nullify(newobj)
281 : return
282 : end if
283 0 : deallocate(nspecies)
284 0 : deallocate(alogsig)
285 0 : deallocate(sigmag)
286 0 : deallocate(f1)
287 0 : deallocate(f2)
288 :
289 0 : end function constructor
290 :
291 : !------------------------------------------------------------------------------
292 : !------------------------------------------------------------------------------
293 0 : subroutine destructor(self)
294 : type(modal_aerosol_properties), intent(inout) :: self
295 :
296 0 : if (allocated(self%exp45logsig_)) then
297 0 : deallocate(self%exp45logsig_)
298 : end if
299 0 : if (allocated(self%voltonumblo_)) then
300 0 : deallocate(self%voltonumblo_)
301 : end if
302 0 : if (allocated(self%voltonumbhi_)) then
303 0 : deallocate(self%voltonumbhi_)
304 : end if
305 :
306 0 : if (allocated(self%sulfate_mode_ndxs_)) then
307 0 : deallocate(self%sulfate_mode_ndxs_)
308 : end if
309 0 : if (allocated(self%dust_mode_ndxs_)) then
310 0 : deallocate(self%dust_mode_ndxs_)
311 : end if
312 0 : if (allocated(self%ssalt_mode_ndxs_)) then
313 0 : deallocate(self%ssalt_mode_ndxs_)
314 : end if
315 0 : if (allocated(self%ammon_mode_ndxs_)) then
316 0 : deallocate(self%ammon_mode_ndxs_)
317 : end if
318 0 : if (allocated(self%nitrate_mode_ndxs_)) then
319 0 : deallocate(self%nitrate_mode_ndxs_)
320 : end if
321 0 : if (allocated(self%msa_mode_ndxs_)) then
322 0 : deallocate(self%msa_mode_ndxs_)
323 : end if
324 0 : if (allocated(self%porganic_mode_ndxs_)) then
325 0 : deallocate(self%porganic_mode_ndxs_)
326 : end if
327 0 : if (allocated(self%sorganic_mode_ndxs_)) then
328 0 : deallocate(self%sorganic_mode_ndxs_)
329 : end if
330 0 : if (allocated(self%bcarbon_mode_ndxs_)) then
331 0 : deallocate(self%bcarbon_mode_ndxs_)
332 : end if
333 :
334 0 : call self%final()
335 :
336 0 : end subroutine destructor
337 :
338 : !------------------------------------------------------------------------------
339 : ! returns number of transported aerosol constituents
340 : !------------------------------------------------------------------------------
341 0 : integer function number_transported(self)
342 : class(modal_aerosol_properties), intent(in) :: self
343 : ! to be implemented later
344 0 : number_transported = -1
345 0 : end function number_transported
346 :
347 : !------------------------------------------------------------------------
348 : ! returns aerosol properties:
349 : ! density
350 : ! hygroscopicity
351 : ! species type
352 : ! species name
353 : ! short wave species refractive indices
354 : ! long wave species refractive indices
355 : ! species morphology
356 : !------------------------------------------------------------------------
357 0 : subroutine get(self, bin_ndx, species_ndx, list_ndx, density, hygro, &
358 : spectype, specname, specmorph, refindex_sw, refindex_lw)
359 :
360 : class(modal_aerosol_properties), intent(in) :: self
361 : integer, intent(in) :: bin_ndx ! bin index
362 : integer, intent(in) :: species_ndx ! species index
363 : integer, optional, intent(in) :: list_ndx ! climate or a diagnostic list number
364 : real(r8), optional, intent(out) :: density ! density (kg/m3)
365 : real(r8), optional, intent(out) :: hygro ! hygroscopicity
366 : character(len=*), optional, intent(out) :: spectype ! species type
367 : character(len=*), optional, intent(out) :: specname ! species name
368 : character(len=*), optional, intent(out) :: specmorph ! species morphology
369 : complex(r8), pointer, optional, intent(out) :: refindex_sw(:) ! short wave species refractive indices
370 : complex(r8), pointer, optional, intent(out) :: refindex_lw(:) ! long wave species refractive indices
371 :
372 : integer :: ilist
373 :
374 0 : if (present(list_ndx)) then
375 0 : ilist = list_ndx
376 : else
377 0 : ilist = 0
378 : end if
379 :
380 : call rad_cnst_get_aer_props(ilist, bin_ndx, species_ndx, &
381 : density_aer=density, hygro_aer=hygro, spectype=spectype, &
382 0 : refindex_aer_sw=refindex_sw, refindex_aer_lw=refindex_lw)
383 :
384 0 : if (present(specname)) then
385 0 : call rad_cnst_get_info(ilist, bin_ndx, species_ndx, spec_name=specname)
386 : end if
387 :
388 0 : if (present(specmorph)) then
389 0 : specmorph = 'UNKNOWN'
390 : end if
391 :
392 0 : end subroutine get
393 :
394 : !------------------------------------------------------------------------
395 : ! returns optics type and table parameters
396 : !------------------------------------------------------------------------
397 0 : subroutine optics_params(self, list_ndx, bin_ndx, opticstype, extpsw, abspsw, asmpsw, absplw, &
398 : refrtabsw, refitabsw, refrtablw, refitablw, ncoef, prefr, prefi, sw_hygro_ext_wtp, &
399 : sw_hygro_ssa_wtp, sw_hygro_asm_wtp, lw_hygro_ext_wtp, wgtpct, nwtp, &
400 : sw_hygro_coreshell_ext, sw_hygro_coreshell_ssa, sw_hygro_coreshell_asm, lw_hygro_coreshell_ext, &
401 : corefrac, bcdust, kap, relh, nfrac, nbcdust, nkap, nrelh )
402 :
403 : class(modal_aerosol_properties), intent(in) :: self
404 : integer, intent(in) :: bin_ndx ! bin index
405 : integer, intent(in) :: list_ndx ! rad climate/diags list
406 :
407 : character(len=*), optional, intent(out) :: opticstype
408 :
409 : ! refactive index table parameters
410 : real(r8), optional, pointer :: extpsw(:,:,:,:) ! short wave specific extinction
411 : real(r8), optional, pointer :: abspsw(:,:,:,:) ! short wave specific absorption
412 : real(r8), optional, pointer :: asmpsw(:,:,:,:) ! short wave asymmetry factor
413 : real(r8), optional, pointer :: absplw(:,:,:,:) ! long wave specific absorption
414 : real(r8), optional, pointer :: refrtabsw(:,:) ! table of short wave real refractive indices for aerosols
415 : real(r8), optional, pointer :: refitabsw(:,:) ! table of short wave imaginary refractive indices for aerosols
416 : real(r8), optional, pointer :: refrtablw(:,:) ! table of long wave real refractive indices for aerosols
417 : real(r8), optional, pointer :: refitablw(:,:) ! table of long wave imaginary refractive indices for aerosols
418 : integer, optional, intent(out) :: ncoef ! number of chebychev polynomials
419 : integer, optional, intent(out) :: prefr ! number of real refractive indices in table
420 : integer, optional, intent(out) :: prefi ! number of imaginary refractive indices in table
421 :
422 : ! hygrowghtpct table parameters
423 : real(r8), optional, pointer :: sw_hygro_ext_wtp(:,:) ! short wave extinction table
424 : real(r8), optional, pointer :: sw_hygro_ssa_wtp(:,:) ! short wave single-scatter albedo table
425 : real(r8), optional, pointer :: sw_hygro_asm_wtp(:,:) ! short wave asymmetry table
426 : real(r8), optional, pointer :: lw_hygro_ext_wtp(:,:) ! long wave absorption table
427 : real(r8), optional, pointer :: wgtpct(:) ! weight precent of H2SO4/H2O solution
428 : integer, optional, intent(out) :: nwtp ! number of weight precent values
429 :
430 : ! hygrocoreshell table parameters
431 : real(r8), optional, pointer :: sw_hygro_coreshell_ext(:,:,:,:,:) ! short wave extinction table
432 : real(r8), optional, pointer :: sw_hygro_coreshell_ssa(:,:,:,:,:) ! short wave single-scatter albedo table
433 : real(r8), optional, pointer :: sw_hygro_coreshell_asm(:,:,:,:,:) ! short wave asymmetry table
434 : real(r8), optional, pointer :: lw_hygro_coreshell_ext(:,:,:,:,:) ! long wave absorption table
435 : real(r8), optional, pointer :: corefrac(:) ! core fraction dimension values
436 : real(r8), optional, pointer :: bcdust(:) ! bc/(bc + dust) fraction dimension values
437 : real(r8), optional, pointer :: kap(:) ! hygroscopicity dimension values
438 : real(r8), optional, pointer :: relh(:) ! relative humidity dimension values
439 : integer, optional, intent(out) :: nfrac ! core fraction dimension size
440 : integer, optional, intent(out) :: nbcdust ! bc/(bc + dust) fraction dimension size
441 : integer, optional, intent(out) :: nkap ! hygroscopicity dimension size
442 : integer, optional, intent(out) :: nrelh ! relative humidity dimension size
443 :
444 : ! refactive index table parameters
445 : call rad_cnst_get_mode_props(list_ndx, bin_ndx, &
446 : opticstype=opticstype, &
447 : extpsw=extpsw, &
448 : abspsw=abspsw, &
449 : asmpsw=asmpsw, &
450 : absplw=absplw, &
451 : refrtabsw=refrtabsw, &
452 : refitabsw=refitabsw, &
453 : refrtablw=refrtablw, &
454 : refitablw=refitablw, &
455 : ncoef=ncoef, &
456 : prefr=prefr, &
457 0 : prefi=prefi)
458 :
459 : ! hygrowghtpct table parameters
460 0 : if (present(sw_hygro_ext_wtp)) then
461 0 : nullify(sw_hygro_ext_wtp)
462 : end if
463 0 : if (present(sw_hygro_ssa_wtp)) then
464 0 : nullify(sw_hygro_ssa_wtp)
465 : end if
466 0 : if (present(sw_hygro_asm_wtp)) then
467 0 : nullify(sw_hygro_asm_wtp)
468 : end if
469 0 : if (present(lw_hygro_ext_wtp)) then
470 0 : nullify(lw_hygro_ext_wtp)
471 : end if
472 0 : if (present(wgtpct)) then
473 0 : nullify(wgtpct)
474 : end if
475 0 : if (present(nwtp)) then
476 0 : nwtp = -1
477 : end if
478 :
479 : ! hygrocoreshell table parameters
480 0 : if (present(sw_hygro_coreshell_ext)) then
481 0 : nullify(sw_hygro_coreshell_ext)
482 : end if
483 0 : if (present(sw_hygro_coreshell_ssa)) then
484 0 : nullify(sw_hygro_coreshell_ssa)
485 : end if
486 0 : if (present(sw_hygro_coreshell_asm)) then
487 0 : nullify(sw_hygro_coreshell_asm)
488 : end if
489 0 : if (present(lw_hygro_coreshell_ext)) then
490 0 : nullify(lw_hygro_coreshell_ext)
491 : end if
492 0 : if (present(corefrac)) then
493 0 : nullify(corefrac)
494 : end if
495 0 : if (present(bcdust)) then
496 0 : nullify(bcdust)
497 : end if
498 0 : if (present(kap)) then
499 0 : nullify(kap)
500 : end if
501 0 : if (present(relh)) then
502 0 : nullify(relh)
503 : end if
504 0 : if (present(nfrac)) then
505 0 : nfrac = -1
506 : end if
507 0 : if (present(nbcdust)) then
508 0 : nbcdust = -1
509 : end if
510 0 : if (present(nkap)) then
511 0 : nkap = -1
512 : end if
513 0 : if (present(nrelh)) then
514 0 : nrelh = -1
515 : end if
516 :
517 0 : end subroutine optics_params
518 :
519 : !------------------------------------------------------------------------------
520 : ! returns radius^3 (m3) of a given bin number
521 : !------------------------------------------------------------------------------
522 0 : pure elemental real(r8) function amcube(self, bin_ndx, volconc, numconc)
523 :
524 : class(modal_aerosol_properties), intent(in) :: self
525 : integer, intent(in) :: bin_ndx ! bin number
526 : real(r8), intent(in) :: volconc ! volume conc (m3/m3)
527 : real(r8), intent(in) :: numconc ! number conc (1/m3)
528 :
529 0 : amcube = (3._r8*volconc/(4._r8*pi*self%exp45logsig_(bin_ndx)*numconc))
530 :
531 0 : end function amcube
532 :
533 : !------------------------------------------------------------------------------
534 : ! returns mass and number activation fractions
535 : !------------------------------------------------------------------------------
536 0 : subroutine actfracs(self, bin_ndx, smc, smax, fn, fm )
537 : use shr_spfn_mod, only: erf => shr_spfn_erf
538 : class(modal_aerosol_properties), intent(in) :: self
539 : integer, intent(in) :: bin_ndx ! bin index
540 : real(r8),intent(in) :: smc ! critical supersaturation for particles of bin radius
541 : real(r8),intent(in) :: smax ! maximum supersaturation for multiple competing aerosols
542 : real(r8),intent(out) :: fn ! activation fraction for aerosol number
543 : real(r8),intent(out) :: fm ! activation fraction for aerosol mass
544 :
545 : real(r8) :: x,y
546 : real(r8), parameter :: twothird = 2._r8/3._r8
547 : real(r8), parameter :: sq2 = sqrt(2._r8)
548 :
549 0 : x=twothird*(log(smc)-log(smax))/(sq2*self%alogsig(bin_ndx))
550 0 : y=x-1.5_r8*sq2*self%alogsig(bin_ndx)
551 :
552 0 : fn = 0.5_r8*(1._r8-erf(x))
553 0 : fm = 0.5_r8*(1._r8-erf(y))
554 :
555 0 : end subroutine actfracs
556 :
557 : !------------------------------------------------------------------------
558 : ! returns constituents names of aerosol number mixing ratios
559 : !------------------------------------------------------------------------
560 0 : subroutine num_names(self, bin_ndx, name_a, name_c)
561 : class(modal_aerosol_properties), intent(in) :: self
562 : integer, intent(in) :: bin_ndx ! bin number
563 : character(len=*), intent(out) :: name_a ! constituent name of ambient aerosol number dens
564 : character(len=*), intent(out) :: name_c ! constituent name of cloud-borne aerosol number dens
565 :
566 0 : call rad_cnst_get_info(0,bin_ndx, num_name=name_a, num_name_cw=name_c)
567 0 : end subroutine num_names
568 :
569 : !------------------------------------------------------------------------
570 : ! returns constituents names of aerosol mass mixing ratios
571 : !------------------------------------------------------------------------
572 0 : subroutine mmr_names(self, bin_ndx, species_ndx, name_a, name_c)
573 : class(modal_aerosol_properties), intent(in) :: self
574 : integer, intent(in) :: bin_ndx ! bin number
575 : integer, intent(in) :: species_ndx ! species number
576 : character(len=*), intent(out) :: name_a ! constituent name of ambient aerosol MMR
577 : character(len=*), intent(out) :: name_c ! constituent name of cloud-borne aerosol MMR
578 :
579 0 : call rad_cnst_get_info(0, bin_ndx, species_ndx, spec_name=name_a, spec_name_cw=name_c)
580 0 : end subroutine mmr_names
581 :
582 : !------------------------------------------------------------------------
583 : ! returns constituent name of ambient aerosol number mixing ratios
584 : !------------------------------------------------------------------------
585 0 : subroutine amb_num_name(self, bin_ndx, name)
586 : class(modal_aerosol_properties), intent(in) :: self
587 : integer, intent(in) :: bin_ndx ! bin number
588 : character(len=*), intent(out) :: name ! constituent name of ambient aerosol number dens
589 :
590 0 : call rad_cnst_get_info(0,bin_ndx, num_name=name)
591 :
592 0 : end subroutine amb_num_name
593 :
594 : !------------------------------------------------------------------------
595 : ! returns constituent name of ambient aerosol mass mixing ratios
596 : !------------------------------------------------------------------------
597 0 : subroutine amb_mmr_name(self, bin_ndx, species_ndx, name)
598 : class(modal_aerosol_properties), intent(in) :: self
599 : integer, intent(in) :: bin_ndx ! bin number
600 : integer, intent(in) :: species_ndx ! species number
601 : character(len=*), intent(out) :: name ! constituent name of ambient aerosol MMR
602 :
603 0 : call rad_cnst_get_info(0, bin_ndx, species_ndx, spec_name=name)
604 :
605 0 : end subroutine amb_mmr_name
606 :
607 : !------------------------------------------------------------------------
608 : ! returns species type
609 : !------------------------------------------------------------------------
610 0 : subroutine species_type(self, bin_ndx, species_ndx, spectype)
611 : class(modal_aerosol_properties), intent(in) :: self
612 : integer, intent(in) :: bin_ndx ! bin number
613 : integer, intent(in) :: species_ndx ! species number
614 : character(len=*), intent(out) :: spectype ! species type
615 :
616 0 : call rad_cnst_get_info(0, bin_ndx, species_ndx, spec_type=spectype)
617 :
618 0 : end subroutine species_type
619 :
620 : !------------------------------------------------------------------------------
621 : ! returns TRUE if Ice Nucleation tendencies are applied to given aerosol bin number
622 : !------------------------------------------------------------------------------
623 0 : function icenuc_updates_num(self, bin_ndx) result(res)
624 : class(modal_aerosol_properties), intent(in) :: self
625 : integer, intent(in) :: bin_ndx ! bin number
626 :
627 : logical :: res
628 :
629 : character(len=aero_name_len) :: spectype
630 : character(len=aero_name_len) :: modetype
631 : integer :: spc_ndx
632 :
633 0 : res = .false.
634 :
635 0 : call rad_cnst_get_info(0, bin_ndx, mode_type=modetype)
636 0 : if (.not.(modetype=='coarse' .or. modetype=='coarse_dust')) then
637 : return
638 : end if
639 :
640 0 : do spc_ndx = 1, self%nspecies(bin_ndx)
641 0 : call self%species_type( bin_ndx, spc_ndx, spectype)
642 0 : if (spectype=='dust') res = .true.
643 : end do
644 :
645 0 : end function icenuc_updates_num
646 :
647 : !------------------------------------------------------------------------------
648 : ! returns TRUE if Ice Nucleation tendencies are applied to a given species within a bin
649 : !------------------------------------------------------------------------------
650 0 : function icenuc_updates_mmr(self, bin_ndx, species_ndx) result(res)
651 : class(modal_aerosol_properties), intent(in) :: self
652 : integer, intent(in) :: bin_ndx ! bin number
653 : integer, intent(in) :: species_ndx ! species number
654 :
655 : logical :: res
656 :
657 : character(len=32) :: spectype
658 : character(len=32) :: modetype
659 :
660 0 : res = .false.
661 :
662 0 : if (species_ndx>0) then
663 :
664 0 : call rad_cnst_get_info(0, bin_ndx, mode_type=modetype)
665 0 : if (.not.(modetype=='coarse' .or. modetype=='coarse_dust')) then
666 : return
667 : end if
668 :
669 0 : call self%species_type( bin_ndx, species_ndx, spectype)
670 0 : if (spectype=='dust') res = .true.
671 : end if
672 :
673 : end function icenuc_updates_mmr
674 :
675 : !------------------------------------------------------------------------------
676 : ! apply max / min to number concentration
677 : !------------------------------------------------------------------------------
678 0 : subroutine apply_number_limits( self, naerosol, vaerosol, istart, istop, m )
679 : class(modal_aerosol_properties), intent(in) :: self
680 : real(r8), intent(inout) :: naerosol(:) ! number conc (1/m3)
681 : real(r8), intent(in) :: vaerosol(:) ! volume conc (m3/m3)
682 : integer, intent(in) :: istart ! start column index (1 <= istart <= istop <= pcols)
683 : integer, intent(in) :: istop ! stop column index
684 : integer, intent(in) :: m ! mode or bin index
685 :
686 : integer :: i
687 :
688 : ! adjust number so that dgnumlo < dgnum < dgnumhi
689 : ! -- the diameter falls within the lower and upper limits which are
690 : ! represented by voltonumhi and voltonumblo values, respectively
691 0 : do i = istart, istop
692 0 : naerosol(i) = max(naerosol(i), vaerosol(i)*self%voltonumbhi_(m))
693 0 : naerosol(i) = min(naerosol(i), vaerosol(i)*self%voltonumblo_(m))
694 : end do
695 :
696 0 : end subroutine apply_number_limits
697 :
698 : !------------------------------------------------------------------------------
699 : ! returns TRUE if species `spc_ndx` in aerosol subset `bin_ndx` contributes to
700 : ! the particles' ability to act as heterogeneous freezing nuclei
701 : !------------------------------------------------------------------------------
702 0 : function hetfrz_species(self, bin_ndx, spc_ndx) result(res)
703 : class(modal_aerosol_properties), intent(in) :: self
704 : integer, intent(in) :: bin_ndx ! bin number
705 : integer, intent(in) :: spc_ndx ! species number
706 :
707 : logical :: res
708 :
709 : character(len=aero_name_len) :: mode_name, species_type
710 :
711 0 : res = .false.
712 :
713 0 : call rad_cnst_get_info(0, bin_ndx, mode_type=mode_name)
714 :
715 0 : if ((trim(mode_name)/='aitken')) then
716 :
717 0 : call self%species_type(bin_ndx, spc_ndx, species_type)
718 :
719 0 : if ((trim(species_type)=='black-c').or.(trim(species_type)=='dust')) then
720 :
721 0 : res = .true.
722 :
723 : end if
724 :
725 : end if
726 :
727 0 : end function hetfrz_species
728 :
729 : !------------------------------------------------------------------------------
730 : ! returns TRUE if soluble
731 : !------------------------------------------------------------------------------
732 0 : logical function soluble(self,bin_ndx)
733 : class(modal_aerosol_properties), intent(in) :: self
734 : integer, intent(in) :: bin_ndx ! bin number
735 :
736 : character(len=aero_name_len) :: mode_name
737 :
738 0 : call rad_cnst_get_info(0, bin_ndx, mode_type=mode_name)
739 :
740 0 : soluble = trim(mode_name)/='primary_carbon'
741 :
742 0 : end function soluble
743 :
744 : !------------------------------------------------------------------------------
745 : ! returns minimum mass mean radius (meters)
746 : !------------------------------------------------------------------------------
747 0 : function min_mass_mean_rad(self,bin_ndx,species_ndx) result(minrad)
748 : class(modal_aerosol_properties), intent(in) :: self
749 : integer, intent(in) :: bin_ndx ! bin number
750 : integer, intent(in) :: species_ndx ! species number
751 :
752 : real(r8) :: minrad ! meters
753 :
754 : integer :: nmodes
755 : character(len=aero_name_len) :: species_type, mode_type
756 :
757 0 : call self%species_type(bin_ndx, species_ndx, spectype=species_type)
758 0 : select case ( trim(species_type) )
759 : case('dust')
760 0 : call rad_cnst_get_info(0, bin_ndx, mode_type=mode_type)
761 0 : select case ( trim(mode_type) )
762 : case ('accum','fine_dust')
763 0 : minrad = 0.258e-6_r8
764 : case ('coarse','coarse_dust')
765 0 : minrad = 1.576e-6_r8
766 : case default
767 0 : minrad = -huge(1._r8)
768 : end select
769 : case('black-c')
770 0 : call rad_cnst_get_info(0, nmodes=nmodes)
771 0 : if (nmodes==3) then
772 : minrad = 0.04e-6_r8
773 : else
774 0 : minrad = 0.067e-6_r8 ! from emission size
775 : endif
776 : case default
777 0 : minrad = -huge(1._r8)
778 : end select
779 :
780 0 : end function min_mass_mean_rad
781 :
782 : !------------------------------------------------------------------------------
783 : ! returns the total number of bins for a given radiation list index
784 : !------------------------------------------------------------------------------
785 0 : function nbins_rlist(self, list_ndx) result(res)
786 : class(modal_aerosol_properties), intent(in) :: self
787 : integer, intent(in) :: list_ndx ! radiation list number
788 :
789 : integer :: res
790 :
791 0 : call rad_cnst_get_info(list_ndx, nmodes=res)
792 :
793 0 : end function nbins_rlist
794 :
795 : !------------------------------------------------------------------------------
796 : ! returns number of species in a bin for a given radiation list index
797 : !------------------------------------------------------------------------------
798 0 : function nspecies_per_bin_rlist(self, list_ndx, bin_ndx) result(res)
799 : class(modal_aerosol_properties), intent(in) :: self
800 : integer, intent(in) :: list_ndx ! radiation list number
801 : integer, intent(in) :: bin_ndx ! bin number
802 :
803 : integer :: res
804 :
805 0 : call rad_cnst_get_info(list_ndx, bin_ndx, nspec=res)
806 :
807 0 : end function nspecies_per_bin_rlist
808 :
809 : !------------------------------------------------------------------------------
810 : ! returns the natural log of geometric standard deviation of the number
811 : ! distribution for radiation list number and aerosol bin
812 : !------------------------------------------------------------------------------
813 0 : function alogsig_rlist(self, list_ndx, bin_ndx) result(res)
814 : class(modal_aerosol_properties), intent(in) :: self
815 : integer, intent(in) :: list_ndx ! radiation list number
816 : integer, intent(in) :: bin_ndx ! bin number
817 :
818 : real(r8) :: res
819 :
820 : real(r8) :: sig
821 :
822 0 : call rad_cnst_get_mode_props(list_ndx, bin_ndx, sigmag=sig)
823 0 : res = log(sig)
824 :
825 0 : end function alogsig_rlist
826 :
827 : !------------------------------------------------------------------------------
828 : ! returns name for a given radiation list number and aerosol bin
829 : !------------------------------------------------------------------------------
830 0 : function bin_name(self, list_ndx, bin_ndx) result(name)
831 : class(modal_aerosol_properties), intent(in) :: self
832 : integer, intent(in) :: list_ndx ! radiation list number
833 : integer, intent(in) :: bin_ndx ! bin number
834 :
835 : character(len=32) name
836 :
837 0 : call rad_cnst_get_info(list_ndx, bin_ndx, mode_type=name)
838 :
839 0 : end function bin_name
840 :
841 : !------------------------------------------------------------------------------
842 : ! returns scavenging diameter (cm) for a given aerosol bin number
843 : !------------------------------------------------------------------------------
844 0 : function scav_diam(self, bin_ndx) result(diam)
845 : use modal_aero_data, only: dgnum_amode
846 :
847 : class(modal_aerosol_properties), intent(in) :: self
848 : integer, intent(in) :: bin_ndx ! bin number
849 :
850 : real(r8) :: diam
851 :
852 0 : diam = dgnum_amode(bin_ndx)
853 :
854 0 : end function scav_diam
855 :
856 : !------------------------------------------------------------------------------
857 : ! adjust aerosol concentration tendencies to create larger sizes of aerosols
858 : ! during resuspension
859 : !------------------------------------------------------------------------------
860 0 : subroutine resuspension_resize(self, dcondt)
861 :
862 0 : use modal_aero_data, only: mode_size_order
863 :
864 : class(modal_aerosol_properties), intent(in) :: self
865 : real(r8), intent(inout) :: dcondt(:)
866 :
867 : integer :: i
868 : character(len=4) :: spcstr
869 :
870 0 : call accumulate_to_larger_mode( 'SO4', self%sulfate_mode_ndxs_, dcondt )
871 0 : call accumulate_to_larger_mode( 'DUST',self%dust_mode_ndxs_,dcondt )
872 0 : call accumulate_to_larger_mode( 'NACL',self%ssalt_mode_ndxs_,dcondt )
873 0 : call accumulate_to_larger_mode( 'MSA', self%msa_mode_ndxs_, dcondt )
874 0 : call accumulate_to_larger_mode( 'NH4', self%ammon_mode_ndxs_, dcondt )
875 0 : call accumulate_to_larger_mode( 'NO3', self%nitrate_mode_ndxs_, dcondt )
876 :
877 0 : spcstr = ' '
878 0 : do i = 1,self%num_soa_
879 0 : write(spcstr,'(i4)') i
880 0 : call accumulate_to_larger_mode( 'SOA'//adjustl(spcstr), self%sorganic_mode_ndxs_(:,i), dcondt )
881 : enddo
882 0 : spcstr = ' '
883 0 : do i = 1,self%num_poa_
884 0 : write(spcstr,'(i4)') i
885 0 : call accumulate_to_larger_mode( 'POM'//adjustl(spcstr), self%porganic_mode_ndxs_(:,i), dcondt )
886 : enddo
887 0 : spcstr = ' '
888 0 : do i = 1,self%num_bc_
889 0 : write(spcstr,'(i4)') i
890 0 : call accumulate_to_larger_mode( 'BC'//adjustl(spcstr), self%bcarbon_mode_ndxs_(:,i), dcondt )
891 : enddo
892 :
893 : contains
894 :
895 : !------------------------------------------------------------------------------
896 0 : subroutine accumulate_to_larger_mode( spc_name, lptr, prevap )
897 :
898 0 : use cam_logfile, only: iulog
899 : use spmd_utils, only: masterproc
900 :
901 : character(len=*), intent(in) :: spc_name
902 : integer, intent(in) :: lptr(:)
903 : real(r8), intent(inout) :: prevap(:)
904 :
905 : integer :: m,n, nl,ns
906 :
907 : logical, parameter :: debug = .false.
908 :
909 : ! find constituent index of the largest mode for the species
910 0 : loop1: do m = 1,self%nbins()-1
911 0 : nl = lptr(mode_size_order(m))
912 0 : if (nl>0) exit loop1
913 : end do loop1
914 :
915 0 : if (.not. nl>0) return
916 :
917 : ! accumulate the smaller modes into the largest mode
918 0 : do n = m+1,self%nbins()
919 0 : ns = lptr(mode_size_order(n))
920 0 : if (ns>0) then
921 0 : prevap(nl) = prevap(nl) + prevap(ns)
922 0 : prevap(ns) = 0._r8
923 : if (masterproc .and. debug) then
924 : write(iulog,'(a,i3,a,i3)') trim(spc_name)//' mode number accumulate ',ns,'->',nl
925 : endif
926 : endif
927 : end do
928 :
929 : end subroutine accumulate_to_larger_mode
930 : !------------------------------------------------------------------------------
931 :
932 : end subroutine resuspension_resize
933 :
934 : !------------------------------------------------------------------------------
935 : ! returns bulk deposition fluxes of the specified species type
936 : ! rebinned to specified diameter limits
937 : !------------------------------------------------------------------------------
938 0 : subroutine rebin_bulk_fluxes(self, bulk_type, dep_fluxes, diam_edges, bulk_fluxes, &
939 : error_code, error_string)
940 : use infnan, only: nan, assignment(=)
941 :
942 : class(modal_aerosol_properties), intent(in) :: self
943 : character(len=*),intent(in) :: bulk_type ! aerosol type to rebin
944 : real(r8), intent(in) :: dep_fluxes(:) ! kg/m2
945 : real(r8), intent(in) :: diam_edges(:) ! meters
946 : real(r8), intent(out) :: bulk_fluxes(:) ! kg/m2
947 : integer, intent(out) :: error_code ! error code (0 if no error)
948 : character(len=*), intent(out) :: error_string ! error string
949 :
950 : real(r8) :: dns_dst ! kg/m3
951 0 : real(r8) :: sigma_g, vmd, tmp, massfrac_bin(size(bulk_fluxes))
952 : real(r8) :: Ntype, Mtype, Mtotal, Ntot
953 : integer :: k,l,m,mm, nbulk
954 : logical :: has_type, type_not_found
955 :
956 : character(len=aero_name_len) :: spectype
957 : character(len=aero_name_len) :: modetype
958 :
959 : real(r8), parameter :: sqrtwo = sqrt(2._r8)
960 : real(r8), parameter :: onethrd = 1._r8/3._r8
961 :
962 0 : error_code = 0
963 0 : error_string = ' '
964 :
965 0 : type_not_found = .true.
966 :
967 0 : nbulk = size(bulk_fluxes)
968 :
969 0 : bulk_fluxes(:) = 0.0_r8
970 :
971 0 : do m = 1,self%nbins()
972 0 : Mtype = 0._r8
973 0 : Mtotal = 0._r8
974 0 : mm = self%indexer(m,0)
975 0 : Ntot = dep_fluxes(mm) ! #/m2
976 :
977 0 : has_type = .false.
978 :
979 0 : do l = 1,self%nspecies(m)
980 0 : mm = self%indexer(m,l)
981 0 : call self%get(m,l, spectype=spectype, density=dns_dst) ! kg/m3
982 0 : if (spectype==bulk_type) then
983 0 : Mtype = dep_fluxes(mm) ! kg/m2
984 0 : has_type = .true.
985 0 : type_not_found = .false.
986 : end if
987 0 : Mtotal = Mtotal + dep_fluxes(mm) ! kg/m2
988 : end do
989 0 : mode_has_type: if (has_type) then
990 0 : call rad_cnst_get_info(0, m, mode_type=modetype)
991 0 : if (Ntot>1.e-40_r8 .and. Mtype>1.e-40_r8 .and. Mtotal>1.e-40_r8) then
992 :
993 0 : call rad_cnst_get_mode_props(0, m, sigmag=sigma_g)
994 0 : tmp = sqrtwo*log(sigma_g)
995 :
996 : ! type number concentration
997 0 : Ntype = Ntot * Mtype/Mtotal ! #/m2
998 :
999 : ! volume median diameter (meters)
1000 0 : vmd = (6._r8*Mtype/(pi*Ntype*dns_dst))**onethrd * exp(1.5_r8*(log(sigma_g))**2)
1001 :
1002 0 : massfrac_bin = 0._r8
1003 :
1004 0 : do k = 1,nbulk
1005 0 : massfrac_bin(k) = 0.5_r8*( erf((log(diam_edges(k+1)/vmd))/tmp) &
1006 0 : - erf((log(diam_edges(k )/vmd))/tmp) )
1007 0 : bulk_fluxes(k) = bulk_fluxes(k) + massfrac_bin(k) * Mtype
1008 : end do
1009 :
1010 : if (debug) then
1011 : if (abs(1._r8-sum(massfrac_bin)) > 1.e-6_r8) then
1012 : write(*,*) 'rebin_bulk_fluxes WARNING mode-num, massfrac_bin, sum(massfrac_bin) = ', &
1013 : m, massfrac_bin, sum(massfrac_bin)
1014 : end if
1015 : end if
1016 :
1017 : end if
1018 : end if mode_has_type
1019 : end do
1020 :
1021 0 : if (type_not_found) then
1022 0 : bulk_fluxes(:) = nan
1023 0 : error_code = 1
1024 0 : write(error_string,*) 'aerosol_properties::rebin_bulk_fluxes ERROR : ',trim(bulk_type),' not found'
1025 : end if
1026 :
1027 0 : end subroutine rebin_bulk_fluxes
1028 :
1029 : !------------------------------------------------------------------------------
1030 : ! Returns TRUE if bin is hydrophilic, otherwise FALSE
1031 : !------------------------------------------------------------------------------
1032 0 : logical function hydrophilic(self, bin_ndx)
1033 : class(modal_aerosol_properties), intent(in) :: self
1034 : integer, intent(in) :: bin_ndx ! bin number
1035 :
1036 : character(len=aero_name_len) :: modetype
1037 :
1038 0 : call rad_cnst_get_info(0, bin_ndx, mode_type=modetype)
1039 :
1040 0 : hydrophilic = (trim(modetype) == 'accum')
1041 :
1042 0 : end function hydrophilic
1043 :
1044 0 : end module modal_aerosol_properties_mod
|
