diff --git a/src/chemistry/aerosol/aerosol_properties_mod.F90 b/src/chemistry/aerosol/aerosol_properties_mod.F90 index e7cea68ad4..127a8758dc 100644 --- a/src/chemistry/aerosol/aerosol_properties_mod.F90 +++ b/src/chemistry/aerosol/aerosol_properties_mod.F90 @@ -74,6 +74,7 @@ module aerosol_properties_mod procedure(aero_resuspension_resize), deferred :: resuspension_resize procedure(aero_rebin_bulk_fluxes), deferred :: rebin_bulk_fluxes procedure(aero_hydrophilic), deferred :: hydrophilic + procedure :: is_bulk procedure :: final=>aero_props_final end type aerosol_properties @@ -101,7 +102,7 @@ end function aero_number_transported ! species morphology !------------------------------------------------------------------------ subroutine aero_props_get(self, bin_ndx, species_ndx, list_ndx, density, hygro, & - spectype, specname, specmorph, refindex_sw, refindex_lw) + spectype, specname, specmorph, refindex_sw, refindex_lw) import :: aerosol_properties, r8 class(aerosol_properties), intent(in) :: self integer, intent(in) :: bin_ndx ! bin index @@ -124,7 +125,10 @@ subroutine aero_optics_params(self, list_ndx, bin_ndx, opticstype, extpsw, absps refrtabsw, refitabsw, refrtablw, refitablw, ncoef, prefr, prefi, sw_hygro_ext_wtp, & sw_hygro_ssa_wtp, sw_hygro_asm_wtp, lw_hygro_ext_wtp, wgtpct, nwtp, & sw_hygro_coreshell_ext, sw_hygro_coreshell_ssa, sw_hygro_coreshell_asm, lw_hygro_coreshell_ext, & - corefrac, bcdust, kap, relh, nfrac, nbcdust, nkap, nrelh ) + corefrac, bcdust, kap, relh, nfrac, nbcdust, nkap, nrelh, & + sw_hygroscopic_ext, sw_hygroscopic_ssa, sw_hygroscopic_asm, lw_hygroscopic_ext, & + sw_insoluble_ext, sw_insoluble_ssa, sw_insoluble_asm, lw_insoluble_ext, & + r_sw_ext, r_sw_scat, r_sw_ascat, r_mu, r_lw_abs ) import :: aerosol_properties, r8 @@ -169,6 +173,25 @@ subroutine aero_optics_params(self, list_ndx, bin_ndx, opticstype, extpsw, absps integer, optional, intent(out) :: nkap ! hygroscopicity dimension size integer, optional, intent(out) :: nrelh ! relative humidity dimension size + ! hygroscopic + real(r8), optional, pointer :: sw_hygroscopic_ext(:,:) ! short wave extinction table + real(r8), optional, pointer :: sw_hygroscopic_ssa(:,:) ! short wave single-scatter albedo table + real(r8), optional, pointer :: sw_hygroscopic_asm(:,:) ! short wave asymmetry table + real(r8), optional, pointer :: lw_hygroscopic_ext(:,:) ! long wave absorption table + + ! non-hygroscopic (insoluble) + real(r8), optional, pointer :: sw_insoluble_ext(:) ! short wave extinction table + real(r8), optional, pointer :: sw_insoluble_ssa(:) ! short wave single-scatter albedo table + real(r8), optional, pointer :: sw_insoluble_asm(:) ! short wave asymmetry table + real(r8), optional, pointer :: lw_insoluble_ext(:) ! long wave absorption table + + ! volcanic radius + real(r8), optional, pointer :: r_sw_ext(:,:) + real(r8), optional, pointer :: r_sw_scat (:,:) + real(r8), optional, pointer :: r_sw_ascat(:,:) + real(r8), optional, pointer :: r_mu(:) + real(r8), optional, pointer :: r_lw_abs(:,:) + end subroutine aero_optics_params !------------------------------------------------------------------------ @@ -375,12 +398,12 @@ end function aero_alogsig_rlist ! returns name for a given radiation list number and aerosol bin !------------------------------------------------------------------------------ function aero_bin_name(self, list_ndx, bin_ndx) result(name) - import :: aerosol_properties, r8 + import :: aerosol_properties, r8, aero_name_len class(aerosol_properties), intent(in) :: self integer, intent(in) :: list_ndx ! radiation list number integer, intent(in) :: bin_ndx ! bin number - character(len=32) name + character(len=aero_name_len) :: name end function aero_bin_name @@ -710,4 +733,14 @@ pure real(r8) function pom_equivso4_factor(self) end function pom_equivso4_factor + !------------------------------------------------------------------------------ + ! returns TRUE if bulk aerosol representation + !------------------------------------------------------------------------------ + pure logical function is_bulk(self) + class(aerosol_properties), intent(in) :: self + + is_bulk = .false. + + end function is_bulk + end module aerosol_properties_mod diff --git a/src/chemistry/aerosol/bulk_aerosol_properties_mod.F90 b/src/chemistry/aerosol/bulk_aerosol_properties_mod.F90 new file mode 100644 index 0000000000..c5f6304b20 --- /dev/null +++ b/src/chemistry/aerosol/bulk_aerosol_properties_mod.F90 @@ -0,0 +1,722 @@ +!-------------------------------------------------------------------------------- +! For bulk aerosol representation. +! Here each aerosol is treated as a separate bin. +!-------------------------------------------------------------------------------- +module bulk_aerosol_properties_mod + use shr_kind_mod, only: r8 => shr_kind_r8 + use cam_abortutils, only: endrun + use string_utils, only : to_lower + + use aerosol_properties_mod, only: aerosol_properties, aero_name_len + + use rad_constituents, only: rad_cnst_get_info, rad_cnst_get_aer_props, rad_cnst_get_aer_mmr + use infnan, only: nan, assignment(=) + + implicit none + + private + + public :: bulk_aerosol_properties + + type, extends(aerosol_properties) :: bulk_aerosol_properties + + private + + contains + + procedure :: number_transported + procedure :: get + procedure :: amcube + procedure :: actfracs + procedure :: num_names + procedure :: mmr_names + procedure :: amb_num_name + procedure :: amb_mmr_name + procedure :: species_type + procedure :: icenuc_updates_num + procedure :: icenuc_updates_mmr + procedure :: apply_number_limits + procedure :: hetfrz_species + procedure :: optics_params + procedure :: nbins_rlist + procedure :: nspecies_per_bin_rlist + procedure :: alogsig_rlist + procedure :: soluble + procedure :: min_mass_mean_rad + procedure :: bin_name + procedure :: scav_diam + procedure :: resuspension_resize + procedure :: rebin_bulk_fluxes + procedure :: hydrophilic + procedure :: is_bulk + + final :: destructor + + end type bulk_aerosol_properties + + interface bulk_aerosol_properties + procedure :: constructor + end interface bulk_aerosol_properties + +contains + + !------------------------------------------------------------------------------ + !------------------------------------------------------------------------------ + function constructor() result(newobj) + + type(bulk_aerosol_properties), pointer :: newobj + + integer,allocatable :: nspecies(:) + real(r8),allocatable :: alogsig(:) + real(r8),allocatable :: f1(:) + integer :: ierr, naero + + allocate(newobj,stat=ierr) + if( ierr /= 0 ) then + nullify(newobj) + return + end if + + call rad_cnst_get_info(0, naero=naero) + + ! Here treat each aerosol as a separate bin + allocate( nspecies(naero),stat=ierr ) + if( ierr /= 0 ) then + nullify(newobj) + return + end if + allocate( alogsig(naero),stat=ierr ) + if( ierr /= 0 ) then + nullify(newobj) + return + end if + allocate( f1(naero),stat=ierr ) + if( ierr /= 0 ) then + nullify(newobj) + return + end if + + ! Bulk aerosols have 1 chemical species in each bin + nspecies(:) = 1 + + ! Taken from CARMA -- not sure if it will be used for our purposes + alogsig(:) = log(2._r8) + f1(:) = 1._r8 + + ! For bulk aerosols, the number of bins and total number of constituents are + ! the same (naero) -- one constituent (species and mass) per bin. + call newobj%initialize(nbin=naero, ncnst=naero, nspec=nspecies, nmasses=nspecies, & + alogsig=alogsig, f1=f1, f2=f1, ierr=ierr) + + deallocate(nspecies) + deallocate(alogsig) + deallocate(f1) + + if( ierr /= 0 ) then + nullify(newobj) + return + end if + + end function constructor + + !------------------------------------------------------------------------------ + !------------------------------------------------------------------------------ + subroutine destructor(self) + type(bulk_aerosol_properties), intent(inout) :: self + + end subroutine destructor + + !------------------------------------------------------------------------------ + ! returns number of transported aerosol constituents + !------------------------------------------------------------------------------ + integer function number_transported(self) + class(bulk_aerosol_properties), intent(in) :: self + ! to be implemented later + number_transported = -1 + end function number_transported + + + !------------------------------------------------------------------------ + ! returns aerosol properties: + ! density + ! hygroscopicity + ! species type + ! species name + ! short wave species refractive indices + ! long wave species refractive indices + ! species morphology + !------------------------------------------------------------------------ + subroutine get(self, bin_ndx, species_ndx, list_ndx, density, hygro, & + spectype, specname, specmorph, refindex_sw, refindex_lw) + + class(bulk_aerosol_properties), intent(in) :: self + integer, intent(in) :: bin_ndx ! bin index + integer, intent(in) :: species_ndx ! species index + integer, optional, intent(in) :: list_ndx ! climate or a diagnostic list number + real(r8), optional, intent(out) :: density ! density (kg/m3) + real(r8), optional, intent(out) :: hygro ! hygroscopicity + character(len=*), optional, intent(out) :: spectype ! species type + character(len=*), optional, intent(out) :: specname ! species name + character(len=*), optional, intent(out) :: specmorph ! species morphology + complex(r8), pointer, optional, intent(out) :: refindex_sw(:) ! short wave species refractive indices + complex(r8), pointer, optional, intent(out) :: refindex_lw(:) ! long wave species refractive indices + + integer :: ilist + character(len=20) :: aername + + if (present(list_ndx)) then + ilist = list_ndx + else + ilist = 0 + end if + + if (present(density)) then + call rad_cnst_get_aer_props(ilist, bin_ndx, density_aer=density) + end if + + if (present(hygro)) then + call rad_cnst_get_aer_props(ilist, bin_ndx, hygro_aer=hygro) + end if + if (present(spectype)) then + + call rad_cnst_get_aer_props(ilist, bin_ndx, aername=aername) + + select case ( to_lower( aername(:4) ) ) + case('dust') + spectype = 'dust' + case('sulf','volc') + spectype = 'sulfate' + case('bcar','bcph') + spectype = 'black-c' + case('ocar','ocph') + spectype = 'p-organic' + case('sslt','seas','ssam','sscm') + spectype = 'seasalt' + case default + spectype = 'UNKNOWN' + call endrun('ERROR: bulk_aerosol_properties_mod%get aername not recognized : '//aername) + end select + + end if + if (present(specmorph)) then + call endrun('ERROR: bulk_aerosol_properties_mod%get specmorph not yet implemented') + end if + if (present(specname)) then + call rad_cnst_get_aer_props(ilist, bin_ndx, aername=specname) + end if + if (present(refindex_sw)) then + call rad_cnst_get_aer_props(ilist, bin_ndx, refindex_aer_sw=refindex_sw) + end if + if (present(refindex_lw)) then + call rad_cnst_get_aer_props(ilist, bin_ndx, refindex_aer_lw=refindex_lw) + end if + + end subroutine get + + !------------------------------------------------------------------------ + ! returns optics type and table parameters + !------------------------------------------------------------------------ + subroutine optics_params(self, list_ndx, bin_ndx, opticstype, extpsw, abspsw, asmpsw, absplw, & + refrtabsw, refitabsw, refrtablw, refitablw, ncoef, prefr, prefi, sw_hygro_ext_wtp, & + sw_hygro_ssa_wtp, sw_hygro_asm_wtp, lw_hygro_ext_wtp, wgtpct, nwtp, & + sw_hygro_coreshell_ext, sw_hygro_coreshell_ssa, sw_hygro_coreshell_asm, lw_hygro_coreshell_ext, & + corefrac, bcdust, kap, relh, nfrac, nbcdust, nkap, nrelh, & + sw_hygroscopic_ext, sw_hygroscopic_ssa, sw_hygroscopic_asm, lw_hygroscopic_ext, & + sw_insoluble_ext, sw_insoluble_ssa, sw_insoluble_asm, lw_insoluble_ext, & + r_sw_ext, r_sw_scat, r_sw_ascat, r_mu, r_lw_abs ) + + + class(bulk_aerosol_properties), intent(in) :: self + integer, intent(in) :: bin_ndx ! bin index + integer, intent(in) :: list_ndx ! rad climate/diags list + + character(len=*), optional, intent(out) :: opticstype + + ! refactive index table parameters + real(r8), optional, pointer :: extpsw(:,:,:,:) ! short wave specific extinction + real(r8), optional, pointer :: abspsw(:,:,:,:) ! short wave specific absorption + real(r8), optional, pointer :: asmpsw(:,:,:,:) ! short wave asymmetry factor + real(r8), optional, pointer :: absplw(:,:,:,:) ! long wave specific absorption + real(r8), optional, pointer :: refrtabsw(:,:) ! table of short wave real refractive indices for aerosols + real(r8), optional, pointer :: refitabsw(:,:) ! table of short wave imaginary refractive indices for aerosols + real(r8), optional, pointer :: refrtablw(:,:) ! table of long wave real refractive indices for aerosols + real(r8), optional, pointer :: refitablw(:,:) ! table of long wave imaginary refractive indices for aerosols + integer, optional, intent(out) :: ncoef ! number of chebychev polynomials + integer, optional, intent(out) :: prefr ! number of real refractive indices in table + integer, optional, intent(out) :: prefi ! number of imaginary refractive indices in table + + ! hygrowghtpct table parameters + real(r8), optional, pointer :: sw_hygro_ext_wtp(:,:) ! short wave extinction table + real(r8), optional, pointer :: sw_hygro_ssa_wtp(:,:) ! short wave single-scatter albedo table + real(r8), optional, pointer :: sw_hygro_asm_wtp(:,:) ! short wave asymmetry table + real(r8), optional, pointer :: lw_hygro_ext_wtp(:,:) ! long wave absorption table + real(r8), optional, pointer :: wgtpct(:) ! weight precent of H2SO4/H2O solution + integer, optional, intent(out) :: nwtp ! number of weight precent values + + ! hygrocoreshell table parameters + real(r8), optional, pointer :: sw_hygro_coreshell_ext(:,:,:,:,:) ! short wave extinction table + real(r8), optional, pointer :: sw_hygro_coreshell_ssa(:,:,:,:,:) ! short wave single-scatter albedo table + real(r8), optional, pointer :: sw_hygro_coreshell_asm(:,:,:,:,:) ! short wave asymmetry table + real(r8), optional, pointer :: lw_hygro_coreshell_ext(:,:,:,:,:) ! long wave absorption table + real(r8), optional, pointer :: corefrac(:) ! core fraction dimension values + real(r8), optional, pointer :: bcdust(:) ! bc/(bc + dust) fraction dimension values + real(r8), optional, pointer :: kap(:) ! hygroscopicity dimension values + real(r8), optional, pointer :: relh(:) ! relative humidity dimension values + integer, optional, intent(out) :: nfrac ! core fraction dimension size + integer, optional, intent(out) :: nbcdust ! bc/(bc + dust) fraction dimension size + integer, optional, intent(out) :: nkap ! hygroscopicity dimension size + integer, optional, intent(out) :: nrelh ! relative humidity dimension size + + ! hygroscopic + real(r8), optional, pointer :: sw_hygroscopic_ext(:,:) ! short wave extinction table + real(r8), optional, pointer :: sw_hygroscopic_ssa(:,:) ! short wave single-scatter albedo table + real(r8), optional, pointer :: sw_hygroscopic_asm(:,:) ! short wave asymmetry table + real(r8), optional, pointer :: lw_hygroscopic_ext(:,:) ! long wave absorption table + + ! non-hygroscopic (insoluble) + real(r8), optional, pointer :: sw_insoluble_ext(:) ! short wave extinction table + real(r8), optional, pointer :: sw_insoluble_ssa(:) ! short wave single-scatter albedo table + real(r8), optional, pointer :: sw_insoluble_asm(:) ! short wave asymmetry table + real(r8), optional, pointer :: lw_insoluble_ext(:) ! long wave absorption table + + ! volcanic radius + real(r8), optional, pointer :: r_sw_ext(:,:) + real(r8), optional, pointer :: r_sw_scat (:,:) + real(r8), optional, pointer :: r_sw_ascat(:,:) + real(r8), optional, pointer :: r_mu(:) + real(r8), optional, pointer :: r_lw_abs(:,:) + + ! refactive index table parameters + call rad_cnst_get_aer_props(list_ndx, bin_ndx, & + opticstype=opticstype, & + sw_hygro_ext=sw_hygroscopic_ext, & + sw_hygro_ssa=sw_hygroscopic_ssa, & + sw_hygro_asm=sw_hygroscopic_asm, & + lw_hygro_ext=lw_hygroscopic_ext, & + sw_nonhygro_ext=sw_insoluble_ext, & + sw_nonhygro_ssa=sw_insoluble_ssa, & + sw_nonhygro_asm=sw_insoluble_asm, & + lw_ext=lw_insoluble_ext, & + r_sw_ext=r_sw_ext, r_sw_scat=r_sw_scat, r_sw_ascat=r_sw_ascat, & + r_lw_abs=r_lw_abs, mu=r_mu ) + + if (present(extpsw)) then + nullify(extpsw) + endif + + if (present(abspsw)) then + nullify(abspsw) + endif + if (present(asmpsw)) then + nullify(asmpsw) + endif + if (present(absplw)) then + nullify(absplw) + endif + if (present(refrtabsw)) then + nullify(refrtabsw) + endif + if (present(refitabsw)) then + nullify(refitabsw) + endif + if (present(refrtablw)) then + nullify(refrtablw) + endif + if (present(refitablw)) then + nullify(refitablw) + endif + if (present(ncoef)) then + ncoef = -huge(1) + endif + if (present(prefr)) then + prefr = -huge(1) + endif + if (present(prefi)) then + prefi = -huge(1) + endif + if (present(sw_hygro_ext_wtp)) then + nullify(sw_hygro_ext_wtp) + endif + if (present(sw_hygro_ssa_wtp)) then + nullify(sw_hygro_ssa_wtp) + endif + if (present(sw_hygro_asm_wtp)) then + nullify(sw_hygro_asm_wtp) + endif + if (present(lw_hygro_ext_wtp)) then + nullify(lw_hygro_ext_wtp) + endif + if (present(wgtpct)) then + nullify(wgtpct) + endif + if (present(nwtp)) then + nwtp = -huge(1) + endif + + if (present(sw_hygro_coreshell_ext)) then + nullify(sw_hygro_coreshell_ext) + endif + if (present(sw_hygro_coreshell_ssa)) then + nullify(sw_hygro_coreshell_ssa) + endif + if (present(sw_hygro_coreshell_asm)) then + nullify(sw_hygro_coreshell_asm) + endif + if (present(lw_hygro_coreshell_ext)) then + nullify(lw_hygro_coreshell_ext) + endif + if (present(corefrac)) then + nullify(corefrac) + endif + if (present(bcdust)) then + nullify(bcdust) + endif + if (present(kap)) then + nullify(kap) + endif + if (present(relh)) then + nullify(relh) + endif + + if (present(nfrac)) then + nfrac = -huge(1) + endif + if (present(nbcdust)) then + nbcdust = -huge(1) + endif + if (present(nkap)) then + nkap = -huge(1) + endif + if (present(nrelh)) then + nrelh = -huge(1) + endif + + end subroutine optics_params + + !------------------------------------------------------------------------------ + ! returns radius^3 (m3) of a given bin number + !------------------------------------------------------------------------------ + pure elemental real(r8) function amcube(self, bin_ndx, volconc, numconc) + + class(bulk_aerosol_properties), intent(in) :: self + integer, intent(in) :: bin_ndx ! bin number + real(r8), intent(in) :: volconc ! volume conc (m3/m3) + real(r8), intent(in) :: numconc ! number conc (1/m3) + + amcube = nan ! to be implemented later if needed + + end function amcube + + !------------------------------------------------------------------------------ + ! returns mass and number activation fractions + !------------------------------------------------------------------------------ + subroutine actfracs(self, bin_ndx, smc, smax, fn, fm ) + + class(bulk_aerosol_properties), intent(in) :: self + integer, intent(in) :: bin_ndx ! bin index + real(r8),intent(in) :: smc ! critical supersaturation for particles of bin radius + real(r8),intent(in) :: smax ! maximum supersaturation for multiple competing aerosols + real(r8),intent(out) :: fn ! activation fraction for aerosol number + real(r8),intent(out) :: fm ! activation fraction for aerosol mass + + ! to be implemented later if needed + call endrun('ERROR: bulk_aerosol_properties_mod%actfracs not yet implemented') + + end subroutine actfracs + + !------------------------------------------------------------------------ + ! returns constituents names of aerosol number mixing ratios + !------------------------------------------------------------------------ + subroutine num_names(self, bin_ndx, name_a, name_c) + class(bulk_aerosol_properties), intent(in) :: self + integer, intent(in) :: bin_ndx ! bin number + character(len=*), intent(out) :: name_a ! constituent name of ambient aerosol number dens + character(len=*), intent(out) :: name_c ! constituent name of cloud-borne aerosol number dens + + ! to be implemented later if needed + call endrun('ERROR: bulk_aerosol_properties_mod%num_names not yet implemented') + + end subroutine num_names + + !------------------------------------------------------------------------ + ! returns constituents names of aerosol mass mixing ratios + !------------------------------------------------------------------------ + subroutine mmr_names(self, bin_ndx, species_ndx, name_a, name_c) + class(bulk_aerosol_properties), intent(in) :: self + integer, intent(in) :: bin_ndx ! bin number + integer, intent(in) :: species_ndx ! species number + character(len=*), intent(out) :: name_a ! constituent name of ambient aerosol MMR + character(len=*), intent(out) :: name_c ! constituent name of cloud-borne aerosol MMR + + ! to be implemented later if needed + call endrun('ERROR: bulk_aerosol_properties_mod%mmr_names not yet implemented') + + end subroutine mmr_names + + !------------------------------------------------------------------------ + ! returns constituent name of ambient aerosol number mixing ratios + !------------------------------------------------------------------------ + subroutine amb_num_name(self, bin_ndx, name) + class(bulk_aerosol_properties), intent(in) :: self + integer, intent(in) :: bin_ndx ! bin number + character(len=*), intent(out) :: name ! constituent name of ambient aerosol number dens + + ! to be implemented later if needed + call endrun('ERROR: bulk_aerosol_properties_mod%amb_num_name not yet implemented') + + end subroutine amb_num_name + + !------------------------------------------------------------------------ + ! returns constituent name of ambient aerosol mass mixing ratios + !------------------------------------------------------------------------ + subroutine amb_mmr_name(self, bin_ndx, species_ndx, name) + class(bulk_aerosol_properties), intent(in) :: self + integer, intent(in) :: bin_ndx ! bin number + integer, intent(in) :: species_ndx ! species number + character(len=*), intent(out) :: name ! constituent name of ambient aerosol MMR + + ! to be implemented later if needed + call endrun('ERROR: bulk_aerosol_properties_mod%amb_mmr_name not yet implemented') + + end subroutine amb_mmr_name + + !------------------------------------------------------------------------ + ! returns species type + !------------------------------------------------------------------------ + subroutine species_type(self, bin_ndx, species_ndx, spectype) + class(bulk_aerosol_properties), intent(in) :: self + integer, intent(in) :: bin_ndx ! bin number + integer, intent(in) :: species_ndx ! species number + character(len=*), intent(out) :: spectype ! species type + + call self%get(bin_ndx, species_ndx, spectype=spectype) + + end subroutine species_type + + !------------------------------------------------------------------------------ + ! returns TRUE if Ice Nucleation tendencies are applied to given aerosol bin number + !------------------------------------------------------------------------------ + function icenuc_updates_num(self, bin_ndx) result(res) + class(bulk_aerosol_properties), intent(in) :: self + integer, intent(in) :: bin_ndx ! bin number + logical :: res + + ! to be implemented later if needed + res = .false. + + end function icenuc_updates_num + + !------------------------------------------------------------------------------ + ! returns TRUE if Ice Nucleation tendencies are applied to a given species within a bin + !------------------------------------------------------------------------------ + function icenuc_updates_mmr(self, bin_ndx, species_ndx) result(res) + class(bulk_aerosol_properties), intent(in) :: self + integer, intent(in) :: bin_ndx ! bin number + integer, intent(in) :: species_ndx ! species number + logical :: res + + ! to be implemented later if needed + res = .false. + end function icenuc_updates_mmr + + !------------------------------------------------------------------------------ + ! apply max / min to number concentration + !------------------------------------------------------------------------------ + subroutine apply_number_limits( self, naerosol, vaerosol, istart, istop, m ) + class(bulk_aerosol_properties), intent(in) :: self + real(r8), intent(inout) :: naerosol(:) ! number conc (1/m3) + real(r8), intent(in) :: vaerosol(:) ! volume conc (m3/m3) + integer, intent(in) :: istart ! start column index (1 <= istart <= istop <= pcols) + integer, intent(in) :: istop ! stop column index + integer, intent(in) :: m ! mode or bin index + + call endrun('ERROR: bulk_aerosol_properties_mod%apply_number_limits not yet implemented') + + end subroutine apply_number_limits + + !------------------------------------------------------------------------------ + ! returns TRUE if species `spc_ndx` in aerosol subset `bin_ndx` contributes to + ! the particles' ability to act as heterogeneous freezing nuclei + !------------------------------------------------------------------------------ + function hetfrz_species(self, bin_ndx, spc_ndx) result(res) + class(bulk_aerosol_properties), intent(in) :: self + integer, intent(in) :: bin_ndx ! bin number + integer, intent(in) :: spc_ndx ! species number + + logical :: res + + ! to be implemented later if needed + res = .false. + end function hetfrz_species + + !------------------------------------------------------------------------------ + ! returns TRUE if soluble + !------------------------------------------------------------------------------ + logical function soluble(self,bin_ndx) + class(bulk_aerosol_properties), intent(in) :: self + integer, intent(in) :: bin_ndx ! bin number + + character(len=20) :: aername + logical :: primary_carbon ! primary carbons (CB1 and OC1) are hydrophobic + + call rad_cnst_get_aer_props(0, bin_ndx, aername=aername) + + aername = to_lower(aername) + + primary_carbon = (aername=='bcpho') .or. (aername=='ocpho') + soluble = .not. primary_carbon + + end function soluble + + !------------------------------------------------------------------------------ + ! returns minimum mass mean radius (meters) + !------------------------------------------------------------------------------ + function min_mass_mean_rad(self,bin_ndx,species_ndx) result(minrad) + class(bulk_aerosol_properties), intent(in) :: self + integer, intent(in) :: bin_ndx ! bin number + integer, intent(in) :: species_ndx ! species number + + real(r8) :: minrad ! meters + + minrad = 0._r8 + + end function min_mass_mean_rad + + !------------------------------------------------------------------------------ + ! returns the total number of bins for a given radiation list index + !------------------------------------------------------------------------------ + function nbins_rlist(self, list_ndx) result(res) + class(bulk_aerosol_properties), intent(in) :: self + integer, intent(in) :: list_ndx ! radiation list number + + integer :: res + + call rad_cnst_get_info(list_ndx, naero=res) + + end function nbins_rlist + + !------------------------------------------------------------------------------ + ! returns number of species in a bin for a given radiation list index + !------------------------------------------------------------------------------ + function nspecies_per_bin_rlist(self, list_ndx, bin_ndx) result(res) + class(bulk_aerosol_properties), intent(in) :: self + integer, intent(in) :: list_ndx ! radiation list number + integer, intent(in) :: bin_ndx ! bin number + + integer :: res + + res = 1 + + end function nspecies_per_bin_rlist + + !------------------------------------------------------------------------------ + ! returns the natural log of geometric standard deviation of the number + ! distribution for radiation list number and aerosol bin + !------------------------------------------------------------------------------ + function alogsig_rlist(self, list_ndx, bin_ndx) result(res) + class(bulk_aerosol_properties), intent(in) :: self + integer, intent(in) :: list_ndx ! radiation list number + integer, intent(in) :: bin_ndx ! bin number + + real(r8) :: res + + ! to be implemented later if needed + res = nan + + end function alogsig_rlist + + !------------------------------------------------------------------------------ + ! returns name for a given radiation list number and aerosol bin + !------------------------------------------------------------------------------ + function bin_name(self, list_ndx, bin_ndx) result(name) + class(bulk_aerosol_properties), intent(in) :: self + integer, intent(in) :: list_ndx ! radiation list number + integer, intent(in) :: bin_ndx ! bin number + + character(len=aero_name_len) :: name + character(len=64), allocatable :: names(:) + integer :: naer, astat + + + call rad_cnst_get_info(list_ndx, naero=naer) + + allocate( names(naer), stat=astat) + if( astat/= 0 ) call endrun('bulk_aerosol_properties_mod%bin_name: names allocate error') + + call rad_cnst_get_info(list_ndx, aernames=names) + + name = names(bin_ndx) + + deallocate(names) + + end function bin_name + + !------------------------------------------------------------------------------ + ! returns scavenging diameter (cm) for a given aerosol bin number + !------------------------------------------------------------------------------ + function scav_diam(self, bin_ndx) result(diam) + + class(bulk_aerosol_properties), intent(in) :: self + integer, intent(in) :: bin_ndx ! bin number + + real(r8) :: diam + + diam = nan ! to be implemented later if needed + + end function scav_diam + + !------------------------------------------------------------------------------ + ! adjust aerosol concentration tendencies to create larger sizes of aerosols + ! during resuspension + !------------------------------------------------------------------------------ + subroutine resuspension_resize(self, dcondt) + + class(bulk_aerosol_properties), intent(in) :: self + real(r8), intent(inout) :: dcondt(:) + + dcondt = nan ! to be implemented later if needed + + end subroutine resuspension_resize + + !------------------------------------------------------------------------------ + ! returns bulk deposition fluxes of the specified species type + ! rebinned to specified diameter limits + !------------------------------------------------------------------------------ + subroutine rebin_bulk_fluxes(self, bulk_type, dep_fluxes, diam_edges, bulk_fluxes, & + error_code, error_string) + + class(bulk_aerosol_properties), intent(in) :: self + character(len=*),intent(in) :: bulk_type ! aerosol type to rebin + real(r8), intent(in) :: dep_fluxes(:) ! kg/m2 + real(r8), intent(in) :: diam_edges(:) ! meters + real(r8), intent(out) :: bulk_fluxes(:) ! kg/m2 + integer, intent(out) :: error_code ! error code (0 if no error) + character(len=*), intent(out) :: error_string ! error string + + ! to be implemented later if needed + call endrun('ERROR: bulk_aerosol_properties_mod%rebin_bulk_fluxes not yet implemented') + + end subroutine rebin_bulk_fluxes + + !------------------------------------------------------------------------------ + ! Returns TRUE if bin is hydrophilic, otherwise FALSE + !------------------------------------------------------------------------------ + logical function hydrophilic(self, bin_ndx) + class(bulk_aerosol_properties), intent(in) :: self + integer, intent(in) :: bin_ndx ! bin number + + hydrophilic = self%soluble(bin_ndx) + + end function hydrophilic + + !------------------------------------------------------------------------------ + ! returns TRUE if bulk aerosol representation + !------------------------------------------------------------------------------ + pure logical function is_bulk(self) + class(bulk_aerosol_properties), intent(in) :: self + + is_bulk = .true. + + end function is_bulk + +end module bulk_aerosol_properties_mod diff --git a/src/chemistry/aerosol/bulk_aerosol_state_mod.F90 b/src/chemistry/aerosol/bulk_aerosol_state_mod.F90 new file mode 100644 index 0000000000..4289adb086 --- /dev/null +++ b/src/chemistry/aerosol/bulk_aerosol_state_mod.F90 @@ -0,0 +1,426 @@ +module bulk_aerosol_state_mod + use shr_kind_mod, only: r8 => shr_kind_r8 + use rad_constituents, only: rad_cnst_get_aer_mmr + use cam_abortutils, only: endrun + + use physics_buffer, only: physics_buffer_desc + use physics_types, only: physics_state + + use aerosol_state_mod, only: aerosol_state, ptr2d_t + use aerosol_properties_mod, only: aerosol_properties + + implicit none + + type, extends(aerosol_state) :: bulk_aerosol_state + private + + type(physics_state), pointer :: state => null() + type(physics_buffer_desc), pointer :: pbuf(:) => null() + + contains + + procedure :: get_transported + procedure :: set_transported + procedure :: ambient_total_bin_mmr + procedure :: get_ambient_mmr_0list + procedure :: get_ambient_mmr_rlist + procedure :: get_cldbrne_mmr + procedure :: get_ambient_num + procedure :: get_cldbrne_num + procedure :: get_states + procedure :: icenuc_size_wght_arr + procedure :: icenuc_size_wght_val + procedure :: icenuc_type_wght + procedure :: update_bin + procedure :: hetfrz_size_wght + procedure :: hygroscopicity + procedure :: water_uptake + procedure :: dry_volume + procedure :: wet_volume + procedure :: water_volume + procedure :: wet_diameter + procedure :: convcld_actfrac + procedure :: wgtpct + + final :: destructor + + end type bulk_aerosol_state + + interface bulk_aerosol_state + procedure :: constructor + end interface bulk_aerosol_state + +contains + + !------------------------------------------------------------------------------ + !------------------------------------------------------------------------------ + function constructor(state,pbuf) result(newobj) + type(physics_state), target :: state + type(physics_buffer_desc), pointer :: pbuf(:) + type(bulk_aerosol_state), pointer :: newobj + + integer :: ierr + + allocate(newobj,stat=ierr) + if( ierr /= 0 ) then + nullify(newobj) + return + end if + + newobj%state => state + newobj%pbuf => pbuf + + end function constructor + + !------------------------------------------------------------------------------ + !------------------------------------------------------------------------------ + subroutine destructor(self) + type(bulk_aerosol_state), intent(inout) :: self + + nullify(self%state) + nullify(self%pbuf) + + end subroutine destructor + + !------------------------------------------------------------------------------ + ! sets transported components + ! This aerosol model with the state of the transported aerosol constituents + ! (mass mixing ratios or number mixing ratios) + !------------------------------------------------------------------------------ + subroutine set_transported( self, transported_array ) + class(bulk_aerosol_state), intent(inout) :: self + real(r8), intent(in) :: transported_array(:,:,:) + ! to be implemented later + end subroutine set_transported + + !------------------------------------------------------------------------------ + ! returns transported components + ! This returns to current state of the transported aerosol constituents + ! (mass mixing ratios or number mixing ratios) + !------------------------------------------------------------------------------ + subroutine get_transported( self, transported_array ) + class(bulk_aerosol_state), intent(in) :: self + real(r8), intent(out) :: transported_array(:,:,:) + ! to be implemented later + end subroutine get_transported + + !------------------------------------------------------------------------ + ! Total aerosol mass mixing ratio for a bin in a given grid box location (column and layer) + !------------------------------------------------------------------------ + function ambient_total_bin_mmr(self, aero_props, bin_ndx, col_ndx, lyr_ndx) result(mmr_tot) + class(bulk_aerosol_state), intent(in) :: self + class(aerosol_properties), intent(in) :: aero_props + integer, intent(in) :: bin_ndx ! bin index + integer, intent(in) :: col_ndx ! column index + integer, intent(in) :: lyr_ndx ! vertical layer index + + real(r8) :: mmr_tot ! mass mixing ratios totaled for all species + real(r8), pointer :: mmr(:,:) ! mass mixing ratios (ncol,nlev) + + call self%get_ambient_mmr(1, bin_ndx, mmr) + + mmr_tot = mmr(col_ndx, lyr_ndx) + + end function ambient_total_bin_mmr + + !------------------------------------------------------------------------------ + ! returns ambient aerosol mass mixing ratio for a given species index and bin index + !------------------------------------------------------------------------------ + subroutine get_ambient_mmr_0list(self, species_ndx, bin_ndx, mmr) + class(bulk_aerosol_state), intent(in) :: self + integer, intent(in) :: species_ndx ! species index + integer, intent(in) :: bin_ndx ! bin index + real(r8), pointer :: mmr(:,:) ! mass mixing ratios (ncol,nlev) + + call self%get_ambient_mmr(0, species_ndx, bin_ndx, mmr) + + end subroutine get_ambient_mmr_0list + + !------------------------------------------------------------------------------ + ! returns ambient aerosol mass mixing ratio for a given radiation diagnostics + ! list index, species index and bin index + !------------------------------------------------------------------------------ + subroutine get_ambient_mmr_rlist(self, list_ndx, species_ndx, bin_ndx, mmr) + class(bulk_aerosol_state), intent(in) :: self + integer, intent(in) :: list_ndx ! rad climate list index + integer, intent(in) :: species_ndx ! species index + integer, intent(in) :: bin_ndx ! bin index + real(r8), pointer :: mmr(:,:) ! mass mixing ratios (ncol,nlev) + + call rad_cnst_get_aer_mmr(list_ndx, bin_ndx, self%state, self%pbuf, mmr) + + end subroutine get_ambient_mmr_rlist + + !------------------------------------------------------------------------------ + ! returns cloud-borne aerosol number mixing ratio for a given species index and bin index + !------------------------------------------------------------------------------ + subroutine get_cldbrne_mmr(self, species_ndx, bin_ndx, mmr) + class(bulk_aerosol_state), intent(in) :: self + integer, intent(in) :: species_ndx ! species index + integer, intent(in) :: bin_ndx ! bin index + real(r8), pointer :: mmr(:,:) ! mass mixing ratios (ncol,nlev) + + call endrun('ERROR: bulk_aerosol_state_mod%get_cldbrne_mmr not yet implemented') + + end subroutine get_cldbrne_mmr + + !------------------------------------------------------------------------------ + ! returns ambient aerosol number mixing ratio for a given species index and bin index + !------------------------------------------------------------------------------ + subroutine get_ambient_num(self, bin_ndx, num) + class(bulk_aerosol_state), intent(in) :: self + integer, intent(in) :: bin_ndx ! bin index + real(r8), pointer :: num(:,:) ! number densities + + nullify(num) + + call endrun('ERROR: bulk_aerosol_state_mod%get_ambient_num not yet implemented') + + end subroutine get_ambient_num + + !------------------------------------------------------------------------------ + ! returns cloud-borne aerosol number mixing ratio for a given species index and bin index + !------------------------------------------------------------------------------ + subroutine get_cldbrne_num(self, bin_ndx, num) + class(bulk_aerosol_state), intent(in) :: self + integer, intent(in) :: bin_ndx ! bin index + real(r8), pointer :: num(:,:) + + nullify(num) + + call endrun('ERROR: bulk_aerosol_state_mod%get_cldbrne_num not yet implemented') + + end subroutine get_cldbrne_num + + !------------------------------------------------------------------------------ + ! returns interstitial and cloud-borne aerosol states + !------------------------------------------------------------------------------ + subroutine get_states( self, aero_props, raer, qqcw ) + class(bulk_aerosol_state), intent(in) :: self + class(aerosol_properties), intent(in) :: aero_props + type(ptr2d_t), intent(out) :: raer(:) + type(ptr2d_t), intent(out) :: qqcw(:) + + call endrun('ERROR: bulk_aerosol_state_mod%get_states not yet implemented') + + end subroutine get_states + + !------------------------------------------------------------------------------ + ! return aerosol bin size weights for a given bin + !------------------------------------------------------------------------------ + subroutine icenuc_size_wght_arr(self, bin_ndx, ncol, nlev, species_type, use_preexisting_ice, wght) + class(bulk_aerosol_state), intent(in) :: self + integer, intent(in) :: bin_ndx ! bin number + integer, intent(in) :: ncol ! number of columns + integer, intent(in) :: nlev ! number of vertical levels + character(len=*), intent(in) :: species_type ! species type + logical, intent(in) :: use_preexisting_ice ! pre-existing ice flag + real(r8), intent(out) :: wght(:,:) + + call endrun('ERROR: bulk_aerosol_state_mod%icenuc_size_wght_arr not yet implemented') + + end subroutine icenuc_size_wght_arr + + !------------------------------------------------------------------------------ + ! return aerosol bin size weights for a given bin, column and vertical layer + !------------------------------------------------------------------------------ + subroutine icenuc_size_wght_val(self, bin_ndx, col_ndx, lyr_ndx, species_type, use_preexisting_ice, wght) + class(bulk_aerosol_state), intent(in) :: self + integer, intent(in) :: bin_ndx ! bin number + integer, intent(in) :: col_ndx ! column index + integer, intent(in) :: lyr_ndx ! vertical layer index + character(len=*), intent(in) :: species_type ! species type + logical, intent(in) :: use_preexisting_ice ! pre-existing ice flag + real(r8), intent(out) :: wght + + call endrun('ERROR: bulk_aerosol_state_mod%icenuc_size_wght_val not yet implemented') + + end subroutine icenuc_size_wght_val + + !------------------------------------------------------------------------------ + ! returns aerosol type weights for a given aerosol type and bin + !------------------------------------------------------------------------------ + subroutine icenuc_type_wght(self, bin_ndx, ncol, nlev, species_type, aero_props, rho, wght, cloud_borne) + + class(bulk_aerosol_state), intent(in) :: self + integer, intent(in) :: bin_ndx ! bin number + integer, intent(in) :: ncol ! number of columns + integer, intent(in) :: nlev ! number of vertical levels + character(len=*), intent(in) :: species_type ! species type + class(aerosol_properties), intent(in) :: aero_props ! aerosol properties object + real(r8), intent(in) :: rho(:,:) ! air density (kg m-3) + real(r8), intent(out) :: wght(:,:) ! type weights + logical, optional, intent(in) :: cloud_borne ! if TRUE cloud-borne aerosols are used + ! otherwise ambient aerosols are used + + call endrun('ERROR: bulk_aerosol_state_mod%icenuc_type_wght not yet implemented') + + end subroutine icenuc_type_wght + + !------------------------------------------------------------------------------ + !------------------------------------------------------------------------------ + subroutine update_bin( self, bin_ndx, col_ndx, lyr_ndx, delmmr_sum, delnum_sum, tnd_ndx, dtime, tend ) + class(bulk_aerosol_state), intent(in) :: self + integer, intent(in) :: bin_ndx ! bin number + integer, intent(in) :: col_ndx ! column index + integer, intent(in) :: lyr_ndx ! vertical layer index + real(r8),intent(in) :: delmmr_sum ! mass mixing ratio change summed over all species in bin + real(r8),intent(in) :: delnum_sum ! number mixing ratio change summed over all species in bin + integer, intent(in) :: tnd_ndx ! tendency index + real(r8),intent(in) :: dtime ! time step size (sec) + real(r8),intent(inout) :: tend(:,:,:) ! tendency + + call endrun('ERROR: bulk_aerosol_state_mod%update_bin not yet implemented') + + end subroutine update_bin + + !------------------------------------------------------------------------------ + ! returns the volume-weighted fractions of aerosol subset `bin_ndx` that can act + ! as heterogeneous freezing nuclei + !------------------------------------------------------------------------------ + function hetfrz_size_wght(self, bin_ndx, ncol, nlev) result(wght) + class(bulk_aerosol_state), intent(in) :: self + integer, intent(in) :: bin_ndx ! bin number + integer, intent(in) :: ncol ! number of columns + integer, intent(in) :: nlev ! number of vertical levels + real(r8) :: wght(ncol,nlev) + + call endrun('ERROR: bulk_aerosol_state_mod%hetfrz_size_wght not yet implemented') + + end function hetfrz_size_wght + + !------------------------------------------------------------------------------ + ! returns hygroscopicity for a given radiation diagnostic list number and + ! bin number + !------------------------------------------------------------------------------ + subroutine hygroscopicity(self, list_ndx, bin_ndx, kappa) + class(bulk_aerosol_state), intent(in) :: self + integer, intent(in) :: list_ndx ! rad climate list number + integer, intent(in) :: bin_ndx ! bin number + real(r8), intent(out) :: kappa(:,:) ! hygroscopicity (ncol,nlev) + + call endrun('ERROR: bulk_aerosol_state_mod%hygroscopicity not yet implemented') + + end subroutine hygroscopicity + + !------------------------------------------------------------------------------ + ! returns aerosol wet diameter and aerosol water concentration for a given + ! radiation diagnostic list number and bin number + !------------------------------------------------------------------------------ + subroutine water_uptake(self, aero_props, list_idx, bin_idx, ncol, nlev, dgnumwet, qaerwat) + + class(bulk_aerosol_state), intent(in) :: self + class(aerosol_properties), intent(in) :: aero_props + integer, intent(in) :: list_idx ! rad climate/diags list number + integer, intent(in) :: bin_idx ! bin number + integer, intent(in) :: ncol ! number of columns + integer, intent(in) :: nlev ! number of levels + real(r8),intent(out) :: dgnumwet(ncol,nlev) ! aerosol wet diameter (m) + real(r8),intent(out) :: qaerwat(ncol,nlev) ! aerosol water concentration (g/g) + + call endrun('ERROR: bulk_aerosol_state_mod%water_uptake not yet implemented') + + end subroutine water_uptake + + !------------------------------------------------------------------------------ + ! aerosol dry volume (m3/kg) for given radiation diagnostic list number and bin number + !------------------------------------------------------------------------------ + function dry_volume(self, aero_props, list_idx, bin_idx, ncol, nlev) result(vol) + + class(bulk_aerosol_state), intent(in) :: self + class(aerosol_properties), intent(in) :: aero_props + + integer, intent(in) :: list_idx ! rad climate/diags list number + integer, intent(in) :: bin_idx ! bin number + integer, intent(in) :: ncol ! number of columns + integer, intent(in) :: nlev ! number of levels + + real(r8) :: vol(ncol,nlev) ! m3/kg + + vol = -huge(1._r8) + + end function dry_volume + + !------------------------------------------------------------------------------ + ! aerosol wet volume (m3/kg) for given radiation diagnostic list number and bin number + !------------------------------------------------------------------------------ + function wet_volume(self, aero_props, list_idx, bin_idx, ncol, nlev) result(vol) + + class(bulk_aerosol_state), intent(in) :: self + class(aerosol_properties), intent(in) :: aero_props + + integer, intent(in) :: list_idx ! rad climate/diags list number + integer, intent(in) :: bin_idx ! bin number + integer, intent(in) :: ncol ! number of columns + integer, intent(in) :: nlev ! number of levels + + real(r8) :: vol(ncol,nlev) ! m3/kg + + vol = -huge(1._r8) + + end function wet_volume + + !------------------------------------------------------------------------------ + ! aerosol water volume (m3/kg) for given radiation diagnostic list number and bin number + !------------------------------------------------------------------------------ + function water_volume(self, aero_props, list_idx, bin_idx, ncol, nlev) result(vol) + + class(bulk_aerosol_state), intent(in) :: self + class(aerosol_properties), intent(in) :: aero_props + + integer, intent(in) :: list_idx ! rad climate/diags list number + integer, intent(in) :: bin_idx ! bin number + integer, intent(in) :: ncol ! number of columns + integer, intent(in) :: nlev ! number of levels + + real(r8) :: vol(ncol,nlev) ! m3/kg + + vol = -huge(1._r8) + + end function water_volume + + !------------------------------------------------------------------------------ + ! aerosol wet diameter + !------------------------------------------------------------------------------ + function wet_diameter(self, bin_idx, ncol, nlev) result(diam) + class(bulk_aerosol_state), intent(in) :: self + integer, intent(in) :: bin_idx ! bin number + integer, intent(in) :: ncol ! number of columns + integer, intent(in) :: nlev ! number of levels + + real(r8) :: diam(ncol,nlev) + + call endrun('ERROR: bulk_aerosol_state_mod%wet_diameter not yet implemented') + + end function wet_diameter + + !------------------------------------------------------------------------------ + ! prescribed aerosol activation fraction for convective cloud + !------------------------------------------------------------------------------ + function convcld_actfrac(self, ibin, ispc, ncol, nlev) result(frac) + + class(bulk_aerosol_state), intent(in) :: self + integer, intent(in) :: ibin ! bin index + integer, intent(in) :: ispc ! species index + integer, intent(in) :: ncol ! number of columns + integer, intent(in) :: nlev ! number of vertical levels + + real(r8) :: frac(ncol,nlev) + + call endrun('ERROR: bulk_aerosol_state_mod%convcld_actfrac not yet implemented') + + end function convcld_actfrac + + !------------------------------------------------------------------------------ + ! aerosol weight precent of H2SO4/H2O solution + !------------------------------------------------------------------------------ + function wgtpct(self, ncol, nlev) result(wtp) + class(bulk_aerosol_state), intent(in) :: self + integer, intent(in) :: ncol, nlev + real(r8) :: wtp(ncol,nlev) ! weight precent of H2SO4/H2O solution for given icol, ilev + + wtp = -huge(1._r8) + + end function wgtpct + +end module bulk_aerosol_state_mod diff --git a/src/chemistry/aerosol/carma_aerosol_properties_mod.F90 b/src/chemistry/aerosol/carma_aerosol_properties_mod.F90 index 222ea38c34..1f0c4d9287 100644 --- a/src/chemistry/aerosol/carma_aerosol_properties_mod.F90 +++ b/src/chemistry/aerosol/carma_aerosol_properties_mod.F90 @@ -318,7 +318,10 @@ subroutine optics_params(self, list_ndx, bin_ndx, opticstype, extpsw, abspsw, as refrtabsw, refitabsw, refrtablw, refitablw, ncoef, prefr, prefi, sw_hygro_ext_wtp, & sw_hygro_ssa_wtp, sw_hygro_asm_wtp, lw_hygro_ext_wtp, wgtpct, nwtp, & sw_hygro_coreshell_ext, sw_hygro_coreshell_ssa, sw_hygro_coreshell_asm, lw_hygro_coreshell_ext, & - corefrac, bcdust, kap, relh, nfrac, nbcdust, nkap, nrelh ) + corefrac, bcdust, kap, relh, nfrac, nbcdust, nkap, nrelh, & + sw_hygroscopic_ext, sw_hygroscopic_ssa, sw_hygroscopic_asm, lw_hygroscopic_ext, & + sw_insoluble_ext, sw_insoluble_ssa, sw_insoluble_asm, lw_insoluble_ext, & + r_sw_ext, r_sw_scat, r_sw_ascat, r_mu, r_lw_abs ) class(carma_aerosol_properties), intent(in) :: self integer, intent(in) :: bin_ndx ! bin index @@ -361,6 +364,25 @@ subroutine optics_params(self, list_ndx, bin_ndx, opticstype, extpsw, abspsw, as integer, optional, intent(out) :: nkap ! hygroscopicity dimension size integer, optional, intent(out) :: nrelh ! relative humidity dimension size + ! hygroscopic + real(r8), optional, pointer :: sw_hygroscopic_ext(:,:) ! short wave extinction table + real(r8), optional, pointer :: sw_hygroscopic_ssa(:,:) ! short wave single-scatter albedo table + real(r8), optional, pointer :: sw_hygroscopic_asm(:,:) ! short wave asymmetry table + real(r8), optional, pointer :: lw_hygroscopic_ext(:,:) ! long wave absorption table + + ! non-hygroscopic (insoluble) + real(r8), optional, pointer :: sw_insoluble_ext(:) ! short wave extinction table + real(r8), optional, pointer :: sw_insoluble_ssa(:) ! short wave single-scatter albedo table + real(r8), optional, pointer :: sw_insoluble_asm(:) ! short wave asymmetry table + real(r8), optional, pointer :: lw_insoluble_ext(:) ! long wave absorption table + + ! volcanic radius + real(r8), optional, pointer :: r_sw_ext(:,:) + real(r8), optional, pointer :: r_sw_scat (:,:) + real(r8), optional, pointer :: r_sw_ascat(:,:) + real(r8), optional, pointer :: r_mu(:) + real(r8), optional, pointer :: r_lw_abs(:,:) + if (present(extpsw)) then nullify(extpsw) end if @@ -416,6 +438,52 @@ subroutine optics_params(self, list_ndx, bin_ndx, opticstype, extpsw, abspsw, as nrelh=nrelh, & nfrac=nfrac ) + + ! hygroscopic + if (present(sw_hygroscopic_ext)) then + nullify(sw_hygroscopic_ext) + end if + if (present(sw_hygroscopic_ssa)) then + nullify(sw_hygroscopic_ssa) + end if + if (present(sw_hygroscopic_asm)) then + nullify(sw_hygroscopic_asm) + end if + if (present(lw_hygroscopic_ext)) then + nullify(lw_hygroscopic_ext) + end if + + ! non-hygroscopic (insoluble) + if (present(sw_insoluble_ext)) then + nullify(sw_insoluble_ext) + end if + if (present(sw_insoluble_ssa)) then + nullify(sw_insoluble_ssa) + end if + if (present(sw_insoluble_asm)) then + nullify(sw_insoluble_asm) + end if + if (present(lw_insoluble_ext)) then + nullify(lw_insoluble_ext) + end if + + ! volcanic radius + if (present(r_sw_ext)) then + nullify(r_sw_ext) + end if + if (present(r_sw_scat)) then + nullify(r_sw_scat) + end if + if (present(r_sw_ascat)) then + nullify(r_sw_ascat) + end if + if (present(r_lw_abs)) then + nullify(r_lw_abs) + end if + if (present(r_mu)) then + nullify(r_mu) + end if + end subroutine optics_params !------------------------------------------------------------------------------ diff --git a/src/chemistry/aerosol/hygro_aerosol_optics_mod.F90 b/src/chemistry/aerosol/hygro_aerosol_optics_mod.F90 new file mode 100644 index 0000000000..d6bf3d4c4e --- /dev/null +++ b/src/chemistry/aerosol/hygro_aerosol_optics_mod.F90 @@ -0,0 +1,160 @@ +!------------------------------------------------------------------------------- +! Short-wave hygroscopic aerosol, Long-wave non-hygroscopic (insoluble) +! aerosol optical properties +!------------------------------------------------------------------------------- +module hygro_aerosol_optics_mod + use shr_kind_mod, only: r8 => shr_kind_r8 + + use aerosol_optics_mod, only: aerosol_optics + use aerosol_properties_mod, only: aerosol_properties + use aerosol_state_mod, only: aerosol_state + + implicit none + + private + + public :: hygro_aerosol_optics + + type, extends(aerosol_optics) :: hygro_aerosol_optics + + ! aerosol optics properties tables (from physprops files) + real(r8), pointer :: ext_sw(:,:) => null() + real(r8), pointer :: ssa_sw(:,:) => null() + real(r8), pointer :: asm_sw(:,:) => null() + real(r8), pointer :: abs_lw(:) => null() + + ! from state + real(r8), allocatable :: wrh(:,:) ! (-) weighting on left side values ! (pcols,pver) + integer , allocatable :: krh(:,:) ! index into rh mesh + + ! aerosol mass mixing ratio + real(r8), pointer :: mmr(:,:) + + contains + + procedure :: sw_props + procedure :: lw_props + + final :: destructor + + end type hygro_aerosol_optics + + interface hygro_aerosol_optics + procedure :: constructor + end interface hygro_aerosol_optics + +contains + + !------------------------------------------------------------------------------ + !------------------------------------------------------------------------------ + function constructor(aero_props, aero_state, ilist, ibin, ncols, nlevs, numrh, relhum) & + result(newobj) + class(aerosol_properties),intent(in) :: aero_props ! aerosol_properties object + class(aerosol_state), intent(in) :: aero_state ! aerosol_state object + integer, intent(in) :: ilist ! climate or a diagnostic list number + integer, intent(in) :: ibin ! bin number + integer, intent(in) :: ncols, nlevs, numrh + real(r8),intent(in) :: relhum(ncols,nlevs) + + type(hygro_aerosol_optics), pointer :: newobj + + real(r8) :: rhtrunc(ncols,nlevs) + integer :: ierr + + allocate(newobj, stat=ierr) + if (ierr/=0) then + nullify(newobj) + return + end if + + allocate(newobj%wrh(ncols,nlevs), stat=ierr) + if (ierr/=0) then + nullify(newobj) + return + end if + + allocate(newobj%krh(ncols,nlevs), stat=ierr) + if (ierr/=0) then + nullify(newobj) + return + end if + +! NOTE should try to use table_interp_mod utility !!! + rhtrunc(1:ncols,1:nlevs) = min(relhum(1:ncols,1:nlevs),1._r8) + newobj%krh(1:ncols,1:nlevs) = min(floor( rhtrunc(1:ncols,1:nlevs) * numrh ) + 1, numrh - 1) ! index into rh mesh + newobj%wrh(1:ncols,1:nlevs) = rhtrunc(1:ncols,1:nlevs) * numrh - newobj%krh(1:ncols,1:nlevs) ! (-) weighting on left side values + + ! optical properties tables + call aero_props%optics_params(ilist, ibin, & + sw_hygroscopic_ext=newobj%ext_sw, & + sw_hygroscopic_ssa=newobj%ssa_sw, & + sw_hygroscopic_asm=newobj%asm_sw, & + lw_insoluble_ext=newobj%abs_lw ) + + call aero_state%get_ambient_mmr(ilist, species_ndx=1, bin_ndx=ibin, mmr=newobj%mmr) + + end function constructor + + !------------------------------------------------------------------------------ + ! returns short wave aerosol optics properties + !------------------------------------------------------------------------------ + subroutine sw_props(self, ncol, ilev, iwav, pext, pabs, palb, pasm) + + class(hygro_aerosol_optics), intent(in) :: self + + integer, intent(in) :: ncol ! number of columns + integer, intent(in) :: ilev ! vertical level index + integer, intent(in) :: iwav ! wave length index + real(r8),intent(out) :: pext(ncol) ! parameterized specific extinction (m2/kg) + real(r8),intent(out) :: pabs(ncol) ! parameterized specific absorption (m2/kg) + real(r8),intent(out) :: palb(ncol) ! parameterized single scattering albedo + real(r8),intent(out) :: pasm(ncol) ! parameterized asymmetry factor + + integer :: icol + + ! interpolate the properties tables + do icol = 1, ncol + pext(icol) = (1._r8 + self%wrh(icol,ilev)) * self%ext_sw(self%krh(icol,ilev)+1,iwav) & + - self%wrh(icol,ilev) * self%ext_sw(self%krh(icol,ilev), iwav) + palb(icol) = (1._r8 + self%wrh(icol,ilev)) * self%ssa_sw(self%krh(icol,ilev)+1,iwav) & + - self%wrh(icol,ilev) * self%ssa_sw(self%krh(icol,ilev), iwav) + pasm(icol) = (1._r8 + self%wrh(icol,ilev)) * self%asm_sw(self%krh(icol,ilev)+1,iwav) & + - self%wrh(icol,ilev) * self%asm_sw(self%krh(icol,ilev), iwav) + + pext(icol) = pext(icol) * self%mmr(icol,ilev) + + pabs(icol) = pext(icol) * ( 1._r8 - palb(icol) ) + + end do + + end subroutine sw_props + + !------------------------------------------------------------------------------ + ! returns long wave aerosol optics properties + !------------------------------------------------------------------------------ + subroutine lw_props(self, ncol, ilev, iwav, pabs) + + class(hygro_aerosol_optics), intent(in) :: self + integer, intent(in) :: ncol ! number of columns + integer, intent(in) :: ilev ! vertical level index + integer, intent(in) :: iwav ! wave length index + real(r8),intent(out) :: pabs(ncol) ! parameterized specific absorption (m2/kg) + + integer :: icol + + pabs(:ncol) = self%abs_lw(iwav) * self%mmr(:ncol,ilev) + + end subroutine lw_props + + !------------------------------------------------------------------------------ + !------------------------------------------------------------------------------ + subroutine destructor(self) + + type(hygro_aerosol_optics), intent(inout) :: self + + deallocate(self%wrh) + deallocate(self%krh) + + end subroutine destructor + +end module hygro_aerosol_optics_mod diff --git a/src/chemistry/aerosol/hygroscopic_aerosol_optics_mod.F90 b/src/chemistry/aerosol/hygroscopic_aerosol_optics_mod.F90 new file mode 100644 index 0000000000..3fa1d01f0e --- /dev/null +++ b/src/chemistry/aerosol/hygroscopic_aerosol_optics_mod.F90 @@ -0,0 +1,165 @@ +!------------------------------------------------------------------------------- +! Short-wave and long-wave hygroscopic aerosol properties +!------------------------------------------------------------------------------- +module hygroscopic_aerosol_optics_mod + use shr_kind_mod, only: r8 => shr_kind_r8 + + use aerosol_optics_mod, only: aerosol_optics + use aerosol_properties_mod, only: aerosol_properties + use aerosol_state_mod, only: aerosol_state + + implicit none + + private + + public :: hygroscopic_aerosol_optics + + type, extends(aerosol_optics) :: hygroscopic_aerosol_optics + + ! aerosol optics properties tables (from physprops files) + real(r8), pointer :: ext_sw(:,:) => null() + real(r8), pointer :: ssa_sw(:,:) => null() + real(r8), pointer :: asm_sw(:,:) => null() + real(r8), pointer :: abs_lw(:,:) => null() + + ! from state + real(r8), allocatable :: wrh(:,:) ! (-) weighting on left side values ! (pcols,pver) + integer , allocatable :: krh(:,:) ! index into rh mesh + + ! aerosol mass mixing ratio + real(r8), pointer :: mmr(:,:) + + contains + + procedure :: sw_props + procedure :: lw_props + + final :: destructor + + end type hygroscopic_aerosol_optics + + interface hygroscopic_aerosol_optics + procedure :: constructor + end interface hygroscopic_aerosol_optics + +contains + + !------------------------------------------------------------------------------ + !------------------------------------------------------------------------------ + function constructor(aero_props, aero_state, ilist, ibin, ncols, nlevs, numrh, relhum) & + result(newobj) + class(aerosol_properties),intent(in) :: aero_props ! aerosol_properties object + class(aerosol_state), intent(in) :: aero_state ! aerosol_state object + integer, intent(in) :: ilist ! climate or a diagnostic list number + integer, intent(in) :: ibin ! bin number + integer, intent(in) :: ncols, nlevs, numrh + real(r8),intent(in) :: relhum(ncols,nlevs) + + type(hygroscopic_aerosol_optics), pointer :: newobj + + real(r8) :: rhtrunc(ncols,nlevs) + integer :: ierr + + allocate(newobj, stat=ierr) + if (ierr/=0) then + nullify(newobj) + return + end if + + allocate(newobj%wrh(ncols,nlevs), stat=ierr) + if (ierr/=0) then + nullify(newobj) + return + end if + + allocate(newobj%krh(ncols,nlevs), stat=ierr) + if (ierr/=0) then + nullify(newobj) + return + end if + +! NOTE should try to use table_interp_mod utility !!! + rhtrunc(1:ncols,1:nlevs) = min(relhum(1:ncols,1:nlevs),1._r8) + newobj%krh(1:ncols,1:nlevs) = min(floor( rhtrunc(1:ncols,1:nlevs) * numrh ) + 1, numrh - 1) ! index into rh mesh + newobj%wrh(1:ncols,1:nlevs) = rhtrunc(1:ncols,1:nlevs) * numrh - newobj%krh(1:ncols,1:nlevs) ! (-) weighting on left side values + + ! optical properties tables + call aero_props%optics_params(ilist, ibin, & + sw_hygroscopic_ext=newobj%ext_sw, & + sw_hygroscopic_ssa=newobj%ssa_sw, & + sw_hygroscopic_asm=newobj%asm_sw, & + lw_hygroscopic_ext=newobj%abs_lw ) + + call aero_state%get_ambient_mmr(ilist, species_ndx=1, bin_ndx=ibin, mmr=newobj%mmr) + + end function constructor + + !------------------------------------------------------------------------------ + ! returns short wave aerosol optics properties + !------------------------------------------------------------------------------ + subroutine sw_props(self, ncol, ilev, iwav, pext, pabs, palb, pasm) + + class(hygroscopic_aerosol_optics), intent(in) :: self + + integer, intent(in) :: ncol ! number of columns + integer, intent(in) :: ilev ! vertical level index + integer, intent(in) :: iwav ! wave length index + real(r8),intent(out) :: pext(ncol) ! parameterized specific extinction (m2/kg) + real(r8),intent(out) :: pabs(ncol) ! parameterized specific absorption (m2/kg) + real(r8),intent(out) :: palb(ncol) ! parameterized single scattering albedo + real(r8),intent(out) :: pasm(ncol) ! parameterized asymmetry factor + + integer :: icol + + ! interpolate the properties tables + do icol = 1, ncol + pext(icol) = (1._r8 + self%wrh(icol,ilev)) * self%ext_sw(self%krh(icol,ilev)+1,iwav) & + - self%wrh(icol,ilev) * self%ext_sw(self%krh(icol,ilev), iwav) + palb(icol) = (1._r8 + self%wrh(icol,ilev)) * self%ssa_sw(self%krh(icol,ilev)+1,iwav) & + - self%wrh(icol,ilev) * self%ssa_sw(self%krh(icol,ilev), iwav) + pasm(icol) = (1._r8 + self%wrh(icol,ilev)) * self%asm_sw(self%krh(icol,ilev)+1,iwav) & + - self%wrh(icol,ilev) * self%asm_sw(self%krh(icol,ilev), iwav) + + pext(icol) = pext(icol) * self%mmr(icol,ilev) + + pabs(icol) = pext(icol) * ( 1._r8 - palb(icol) ) + + end do + + end subroutine sw_props + + !------------------------------------------------------------------------------ + ! returns long wave aerosol optics properties + !------------------------------------------------------------------------------ + subroutine lw_props(self, ncol, ilev, iwav, pabs) + + class(hygroscopic_aerosol_optics), intent(in) :: self + integer, intent(in) :: ncol ! number of columns + integer, intent(in) :: ilev ! vertical level index + integer, intent(in) :: iwav ! wave length index + real(r8),intent(out) :: pabs(ncol) ! parameterized specific absorption (m2/kg) + + integer :: icol + + ! interpolate the properties tables + do icol = 1, ncol + pabs(icol) = (1._r8 + self%wrh(icol,ilev)) * self%abs_lw(self%krh(icol,ilev)+1,iwav) & + - self%wrh(icol,ilev) * self%abs_lw(self%krh(icol,ilev), iwav) + + pabs(icol) = pabs(icol) * self%mmr(icol,ilev) + end do + + end subroutine lw_props + + !------------------------------------------------------------------------------ + !------------------------------------------------------------------------------ + subroutine destructor(self) + + type(hygroscopic_aerosol_optics), intent(inout) :: self + + deallocate(self%wrh) + deallocate(self%krh) + + end subroutine destructor + +end module hygroscopic_aerosol_optics_mod diff --git a/src/chemistry/aerosol/insoluble_aerosol_optics_mod.F90 b/src/chemistry/aerosol/insoluble_aerosol_optics_mod.F90 new file mode 100644 index 0000000000..ac7475dc33 --- /dev/null +++ b/src/chemistry/aerosol/insoluble_aerosol_optics_mod.F90 @@ -0,0 +1,118 @@ +!------------------------------------------------------------------------------- +! Insoluble (non-hygroscopic) aerosol optical properties +!------------------------------------------------------------------------------- +module insoluble_aerosol_optics_mod + use shr_kind_mod, only: r8 => shr_kind_r8 + + use aerosol_optics_mod, only: aerosol_optics + use aerosol_properties_mod, only: aerosol_properties + use aerosol_state_mod, only: aerosol_state + + implicit none + + private + + public :: insoluble_aerosol_optics + + type, extends(aerosol_optics) :: insoluble_aerosol_optics + real(r8), pointer :: lw_abs(:) + real(r8), pointer :: sw_ext(:) + real(r8), pointer :: sw_ssa(:) + real(r8), pointer :: sw_asm(:) + + ! aerosol mass mixing ratio + real(r8), pointer :: mmr(:,:) + + contains + + procedure :: sw_props + procedure :: lw_props + + final :: destructor + + end type insoluble_aerosol_optics + + interface insoluble_aerosol_optics + procedure :: constructor + end interface insoluble_aerosol_optics + +contains + + !------------------------------------------------------------------------------ + !------------------------------------------------------------------------------ + function constructor(aero_props, aero_state, ilist, ibin) result(newobj) + + class(aerosol_properties),intent(in) :: aero_props ! aerosol_properties object + class(aerosol_state), intent(in) :: aero_state ! aerosol_state object + integer, intent(in) :: ilist ! climate or a diagnostic list number + integer, intent(in) :: ibin ! bin number + + type(insoluble_aerosol_optics), pointer :: newobj + + integer :: ierr + + allocate(newobj, stat=ierr) + if (ierr/=0) then + nullify(newobj) + return + end if + + ! get mode properties + call aero_props%optics_params(ilist, ibin, & + sw_insoluble_ext=newobj%sw_ext, & + sw_insoluble_ssa=newobj%sw_ssa, & + sw_insoluble_asm=newobj%sw_asm, & + lw_insoluble_ext=newobj%lw_abs ) + + call aero_state%get_ambient_mmr(ilist, species_ndx=1, bin_ndx=ibin, mmr=newobj%mmr) + + end function constructor + + + !------------------------------------------------------------------------------ + !------------------------------------------------------------------------------ + subroutine destructor(self) + + type(insoluble_aerosol_optics), intent(inout) :: self + + end subroutine destructor + + !------------------------------------------------------------------------------ + ! returns short wave aerosol optics properties + !------------------------------------------------------------------------------ + subroutine sw_props(self, ncol, ilev, iwav, pext, pabs, palb, pasm) + + class(insoluble_aerosol_optics), intent(in) :: self + + integer, intent(in) :: ncol ! number of columns + integer, intent(in) :: ilev ! vertical level index + integer, intent(in) :: iwav ! wave length index + real(r8),intent(out) :: pext(ncol) ! parameterized specific extinction (m2/kg) + real(r8),intent(out) :: pabs(ncol) ! parameterized specific absorption (m2/kg) + real(r8),intent(out) :: palb(ncol) ! parameterized single scattering albedo + real(r8),intent(out) :: pasm(ncol) ! parameterized asymmetry factor + + pext(:ncol) = self%sw_ext(iwav) * self%mmr(:ncol,ilev) + palb(:ncol) = self%sw_ssa(iwav) + pasm(:ncol) = self%sw_asm(iwav) + + pabs(:ncol) = pext(:ncol) * ( 1._r8 - palb(:ncol) ) + + end subroutine sw_props + + !------------------------------------------------------------------------------ + ! returns long wave aerosol optics properties + !------------------------------------------------------------------------------ + subroutine lw_props(self, ncol, ilev, iwav, pabs) + + class(insoluble_aerosol_optics), intent(in) :: self + integer, intent(in) :: ncol ! number of columns + integer, intent(in) :: ilev ! vertical level index + integer, intent(in) :: iwav ! wave length index + real(r8),intent(out) :: pabs(ncol) ! parameterized specific absorption (m2/kg) + + pabs(:ncol) = self%lw_abs(iwav) * self%mmr(:ncol,ilev) + + end subroutine lw_props + +end module insoluble_aerosol_optics_mod diff --git a/src/chemistry/aerosol/modal_aerosol_properties_mod.F90 b/src/chemistry/aerosol/modal_aerosol_properties_mod.F90 index 828b54ed99..5a92b8df5a 100644 --- a/src/chemistry/aerosol/modal_aerosol_properties_mod.F90 +++ b/src/chemistry/aerosol/modal_aerosol_properties_mod.F90 @@ -398,7 +398,10 @@ subroutine optics_params(self, list_ndx, bin_ndx, opticstype, extpsw, abspsw, as refrtabsw, refitabsw, refrtablw, refitablw, ncoef, prefr, prefi, sw_hygro_ext_wtp, & sw_hygro_ssa_wtp, sw_hygro_asm_wtp, lw_hygro_ext_wtp, wgtpct, nwtp, & sw_hygro_coreshell_ext, sw_hygro_coreshell_ssa, sw_hygro_coreshell_asm, lw_hygro_coreshell_ext, & - corefrac, bcdust, kap, relh, nfrac, nbcdust, nkap, nrelh ) + corefrac, bcdust, kap, relh, nfrac, nbcdust, nkap, nrelh, & + sw_hygroscopic_ext, sw_hygroscopic_ssa, sw_hygroscopic_asm, lw_hygroscopic_ext, & + sw_insoluble_ext, sw_insoluble_ssa, sw_insoluble_asm, lw_insoluble_ext, & + r_sw_ext, r_sw_scat, r_sw_ascat, r_mu, r_lw_abs ) class(modal_aerosol_properties), intent(in) :: self integer, intent(in) :: bin_ndx ! bin index @@ -441,6 +444,25 @@ subroutine optics_params(self, list_ndx, bin_ndx, opticstype, extpsw, abspsw, as integer, optional, intent(out) :: nkap ! hygroscopicity dimension size integer, optional, intent(out) :: nrelh ! relative humidity dimension size + ! hygroscopic + real(r8), optional, pointer :: sw_hygroscopic_ext(:,:) ! short wave extinction table + real(r8), optional, pointer :: sw_hygroscopic_ssa(:,:) ! short wave single-scatter albedo table + real(r8), optional, pointer :: sw_hygroscopic_asm(:,:) ! short wave asymmetry table + real(r8), optional, pointer :: lw_hygroscopic_ext(:,:) ! long wave absorption table + + ! non-hygroscopic (insoluble) + real(r8), optional, pointer :: sw_insoluble_ext(:) ! short wave extinction table + real(r8), optional, pointer :: sw_insoluble_ssa(:) ! short wave single-scatter albedo table + real(r8), optional, pointer :: sw_insoluble_asm(:) ! short wave asymmetry table + real(r8), optional, pointer :: lw_insoluble_ext(:) ! long wave absorption table + + ! volcanic radius + real(r8), optional, pointer :: r_sw_ext(:,:) + real(r8), optional, pointer :: r_sw_scat (:,:) + real(r8), optional, pointer :: r_sw_ascat(:,:) + real(r8), optional, pointer :: r_mu(:) + real(r8), optional, pointer :: r_lw_abs(:,:) + ! refactive index table parameters call rad_cnst_get_mode_props(list_ndx, bin_ndx, & opticstype=opticstype, & @@ -514,6 +536,51 @@ subroutine optics_params(self, list_ndx, bin_ndx, opticstype, extpsw, abspsw, as nrelh = -1 end if + ! hygroscopic + if (present(sw_hygroscopic_ext)) then + nullify(sw_hygroscopic_ext) + end if + if (present(sw_hygroscopic_ssa)) then + nullify(sw_hygroscopic_ssa) + end if + if (present(sw_hygroscopic_asm)) then + nullify(sw_hygroscopic_asm) + end if + if (present(lw_hygroscopic_ext)) then + nullify(lw_hygroscopic_ext) + end if + + ! non-hygroscopic (insoluble) + if (present(sw_insoluble_ext)) then + nullify(sw_insoluble_ext) + end if + if (present(sw_insoluble_ssa)) then + nullify(sw_insoluble_ssa) + end if + if (present(sw_insoluble_asm)) then + nullify(sw_insoluble_asm) + end if + if (present(lw_insoluble_ext)) then + nullify(lw_insoluble_ext) + end if + + ! volcanic radius + if (present(r_sw_ext)) then + nullify(r_sw_ext) + end if + if (present(r_sw_scat)) then + nullify(r_sw_scat) + end if + if (present(r_sw_ascat)) then + nullify(r_sw_ascat) + end if + if (present(r_lw_abs)) then + nullify(r_lw_abs) + end if + if (present(r_mu)) then + nullify(r_mu) + end if + end subroutine optics_params !------------------------------------------------------------------------------ diff --git a/src/chemistry/aerosol/volcrad_aerosol_optics_mod.F90 b/src/chemistry/aerosol/volcrad_aerosol_optics_mod.F90 new file mode 100644 index 0000000000..45227dd363 --- /dev/null +++ b/src/chemistry/aerosol/volcrad_aerosol_optics_mod.F90 @@ -0,0 +1,202 @@ +!------------------------------------------------------------------------------- +! Geometric mean radius parameterized optical properties for volcanic +! stratospheric aerosols +!------------------------------------------------------------------------------- +module volcrad_aerosol_optics_mod + use shr_kind_mod, only: r8 => shr_kind_r8 + + use aerosol_optics_mod, only: aerosol_optics + use aerosol_properties_mod, only: aerosol_properties + use aerosol_state_mod, only: aerosol_state + + implicit none + + private + + public :: volcrad_aerosol_optics + + type, extends(aerosol_optics) :: volcrad_aerosol_optics + + ! aerosol optics properties tables (from physprops files) + real(r8), pointer :: r_sw_ext(:,:) => null() + real(r8), pointer :: r_sw_scat(:,:) => null() + real(r8), pointer :: r_sw_ascat(:,:) => null() + real(r8), pointer :: r_lw_abs(:,:) => null() + real(r8), pointer :: r_mu(:) + + ! from state + real(r8), allocatable :: wmu(:,:) ! (-) weighting on left side values ! (pcols,pver) + integer , allocatable :: kmu(:,:) ! index into rh mesh + + ! aerosol mass mixing ratio + real(r8), pointer :: mmr(:,:) + + contains + + procedure :: sw_props + procedure :: lw_props + + final :: destructor + + end type volcrad_aerosol_optics + + interface volcrad_aerosol_optics + procedure :: constructor + end interface volcrad_aerosol_optics + +contains + + !------------------------------------------------------------------------------ + !------------------------------------------------------------------------------ + function constructor(aero_props, aero_state, ilist, ibin, ncols, nlevs, geometric_radius) & + result(newobj) + class(aerosol_properties),intent(in) :: aero_props ! aerosol_properties object + class(aerosol_state), intent(in) :: aero_state ! aerosol_state object + integer, intent(in) :: ilist ! climate or a diagnostic list number + integer, intent(in) :: ibin ! bin number + integer, intent(in) :: ncols, nlevs + real(r8),intent(in) :: geometric_radius(ncols,nlevs) + + type(volcrad_aerosol_optics), pointer :: newobj + + integer :: ierr, nmu, i, k + real(r8) :: r_mu_min, r_mu_max, mutrunc, mu + + allocate(newobj, stat=ierr) + if (ierr/=0) then + nullify(newobj) + return + end if + + allocate(newobj%wmu(ncols,nlevs), stat=ierr) + if (ierr/=0) then + nullify(newobj) + return + end if + + allocate(newobj%kmu(ncols,nlevs), stat=ierr) + if (ierr/=0) then + nullify(newobj) + return + end if + + ! optical properties tables + call aero_props%optics_params(ilist, ibin, & + r_sw_ext=newobj%r_sw_ext, & + r_sw_scat=newobj%r_sw_scat, & + r_sw_ascat=newobj%r_sw_ascat, & + r_lw_abs=newobj%r_lw_abs, & + r_mu=newobj%r_mu ) + + call aero_state%get_ambient_mmr(ilist, species_ndx=1, bin_ndx=ibin, mmr=newobj%mmr) + + +! NOTE should try to use table_interp_mod utility !!! + + nmu = size(newobj%r_mu) + r_mu_max = newobj%r_mu(nmu) + r_mu_min = newobj%r_mu(1) + + do i = 1, ncols + do k = 1, nlevs + if(geometric_radius(i,k) > 0._r8) then + mu = log(geometric_radius(i,k)) + else + mu = 0._r8 + endif + + ASSOCIATE ( kmu=>newobj%kmu(i,k), wmu=>newobj%wmu(i,k), r_mu=>newobj%r_mu ) + + mutrunc = max(min(mu,r_mu_max),r_mu_min) + kmu = max(min(1 + (mutrunc-r_mu_min)/(r_mu_max-r_mu_min)*(nmu-1),nmu-1._r8),1._r8) + wmu = max(min( (mutrunc -r_mu(kmu)) / (r_mu(kmu+1) - r_mu(kmu)) ,1._r8),0._r8) + + END ASSOCIATE + + end do + end do + + end function constructor + + !------------------------------------------------------------------------------ + ! returns short wave aerosol optics properties + !------------------------------------------------------------------------------ + subroutine sw_props(self, ncol, ilev, iwav, pext, pabs, palb, pasm) + + class(volcrad_aerosol_optics), intent(in) :: self + + integer, intent(in) :: ncol ! number of columns + integer, intent(in) :: ilev ! vertical level index + integer, intent(in) :: iwav ! wave length index + real(r8),intent(out) :: pext(ncol) ! parameterized specific extinction (m2/kg) + real(r8),intent(out) :: pabs(ncol) ! parameterized specific absorption (m2/kg) + real(r8),intent(out) :: palb(ncol) ! parameterized single scattering albedo + real(r8),intent(out) :: pasm(ncol) ! parameterized asymmetry factor + + real(r8) :: scat(ncol) + real(r8) :: ascat(ncol) + integer :: icol + + ! interpolate the properties tables + do icol = 1, ncol + + pext(icol) = ((1._r8 - self%wmu(icol,ilev)) * self%r_sw_ext(iwav, self%kmu(icol,ilev) ) + & + (self%wmu(icol,ilev)) * self%r_sw_ext(iwav, self%kmu(icol,ilev)+1)) + + scat(icol) = ((1._r8 - self%wmu(icol,ilev)) * self%r_sw_scat(iwav, self%kmu(icol,ilev) ) + & + (self%wmu(icol,ilev)) * self%r_sw_scat(iwav, self%kmu(icol,ilev)+1)) + + ascat(icol) = ((1._r8 - self%wmu(icol,ilev)) * self%r_sw_ascat(iwav, self%kmu(icol,ilev) ) + & + (self%wmu(icol,ilev)) * self%r_sw_ascat(iwav, self%kmu(icol,ilev)+1)) + + + palb(icol) = scat(icol) / pext(icol) + + if (scat(icol)>0._r8) then + pasm(icol) = ascat(icol) / scat(icol) + else + pasm(icol) = 0._r8 + end if + + pext(icol) = pext(icol) * self%mmr(icol,ilev) + + pabs(icol) = pext(icol) * ( 1._r8 - palb(icol) ) + + end do + + end subroutine sw_props + + !------------------------------------------------------------------------------ + ! returns long wave aerosol optics properties + !------------------------------------------------------------------------------ + subroutine lw_props(self, ncol, ilev, iwav, pabs) + + class(volcrad_aerosol_optics), intent(in) :: self + integer, intent(in) :: ncol ! number of columns + integer, intent(in) :: ilev ! vertical level index + integer, intent(in) :: iwav ! wave length index + real(r8),intent(out) :: pabs(ncol) ! parameterized specific absorption (m2/kg) + + integer :: icol + + ! interpolate the properties tables + do icol = 1, ncol + pabs(icol) = ((1._r8 - self%wmu(icol,ilev)) * self%r_lw_abs(iwav, self%kmu(icol,ilev) ) + & + (self%wmu(icol,ilev)) * self%r_lw_abs(iwav, self%kmu(icol,ilev)+1)) + pabs(icol) = pabs(icol) * self%mmr(icol,ilev) + end do + + end subroutine lw_props + + !------------------------------------------------------------------------------ + !------------------------------------------------------------------------------ + subroutine destructor(self) + + type(volcrad_aerosol_optics), intent(inout) :: self + + deallocate(self%wmu) + deallocate(self%kmu) + + end subroutine destructor + +end module volcrad_aerosol_optics_mod diff --git a/src/physics/cam/aer_rad_props.F90 b/src/physics/cam/aer_rad_props.F90 index 33ecc9056f..435c94f7b1 100644 --- a/src/physics/cam/aer_rad_props.F90 +++ b/src/physics/cam/aer_rad_props.F90 @@ -19,25 +19,27 @@ module aer_rad_props use aerosol_optics_cam,only: aerosol_optics_cam_init, aerosol_optics_cam_sw, aerosol_optics_cam_lw use cam_history, only: fieldname_len, addfld, outfld, add_default, horiz_only use cam_history_support, only : fillvalue -! Placed here due to PGI bug. -use ref_pres, only: clim_modal_aero_top_lev use cam_abortutils, only: endrun +use aer_vis_diag_mod, only: aer_vis_diag_init, aer_vis_diag_out implicit none private save -integer :: top_lev = 1 - public :: & aer_rad_props_init, & aer_rad_props_sw, & ! return SW optical props of aerosols aer_rad_props_lw ! return LW optical props of aerosols -! Private data -character(len=fieldname_len), pointer :: odv_names(:) ! outfld names for visible OD +! debug diags +character(len=32) :: sw_ta_name(nswbands) +character(len=32) :: sw_wa_name(nswbands) +character(len=32) :: sw_ga_name(nswbands) +character(len=32) :: sw_fa_name(nswbands) +character(len=32) :: lw_ta_name(nlwbands) +logical :: bam_debug = .false. !============================================================================== contains @@ -46,26 +48,20 @@ module aer_rad_props subroutine aer_rad_props_init() use phys_control, only: phys_getopts - - integer :: i integer :: numaerosols ! number of aerosols - character(len=64), pointer :: aernames(:) ! aerosol names logical :: history_amwg ! output the variables used by the AMWG diag package logical :: history_aero_optics ! Output aerosol optics diagnostics logical :: history_dust ! Output dust diagnostics - logical :: prog_modal_aero ! Prognostic modal aerosols present integer :: nmodes ! number of aerosol modes integer :: nbins ! number of aerosol bins + character(len=2) :: numch + integer :: i !---------------------------------------------------------------------------- call phys_getopts( history_aero_optics_out = history_aero_optics, & history_amwg_out = history_amwg, & - history_dust_out = history_dust, & - prog_modal_aero_out = prog_modal_aero ) - - ! Limit modal aerosols with top_lev here. - if (prog_modal_aero) top_lev = clim_modal_aero_top_lev + history_dust_out = history_dust ) call addfld ('AEROD_v ', horiz_only, 'A', '1', & 'Total Aerosol Optical Depth in visible band', flag_xyfill=.true.) @@ -76,39 +72,50 @@ subroutine aer_rad_props_init() ! Contributions to AEROD_v from individual aerosols (climate species). ! number of bulk aerosols in climate list - call rad_cnst_get_info(0, naero=numaerosols) + call rad_cnst_get_info(0, naero=numaerosols, nmodes=nmodes, nbins=nbins) - ! get names of bulk aerosols - allocate(aernames(numaerosols)) - call rad_cnst_get_info(0, aernames=aernames, nmodes=nmodes, nbins=nbins) - - ! diagnostic output for bulk aerosols - ! create outfld names for visible OD - allocate(odv_names(numaerosols)) - do i = 1, numaerosols - odv_names(i) = 'ODV_'//trim(aernames(i)) - call addfld (odv_names(i), horiz_only, 'A', '1', & - trim(aernames(i))//' optical depth in visible band', flag_xyfill=.true.) - end do + call aer_vis_diag_init() ! Determine default fields - if (history_amwg .or. history_dust ) then + if (history_amwg .or. history_dust .or. history_aero_optics) then call add_default ('AEROD_v', 1, ' ') endif - if ( history_aero_optics ) then - call add_default ('AEROD_v', 1, ' ') - do i = 1, numaerosols - odv_names(i) = 'ODV_'//trim(aernames(i)) - call add_default (odv_names(i), 1, ' ') - end do - endif - - if (nmodes>0 .or. nbins>0) then + if (numaerosols>0 .or. nmodes>0 .or. nbins>0) then call aerosol_optics_cam_init() - end if + endif + + bam_debug = numaerosols>0 .and. history_aero_optics + + ! for BAM debugging + if (bam_debug) then + call add_default ('AEROD_v', 2, ' ') + call add_default ('AODvstrt', 2, ' ') + + ! debug diags + do i = 1,nswbands + write(numch,'(I2.2)') i + sw_ta_name(i) = 'TASW'//numch + sw_wa_name(i) = 'WASW'//numch + sw_ga_name(i) = 'GASW'//numch + sw_fa_name(i) = 'FASW'//numch + call addfld(sw_ta_name(i), (/'lev'/), 'A',' ', 'SW TAU for wave band '//numch, flag_xyfill=.true.) + call addfld(sw_wa_name(i), (/'lev'/), 'A',' ', 'SW WA for wave band '//numch, flag_xyfill=.true.) + call addfld(sw_ga_name(i), (/'lev'/), 'A',' ', 'SW GA for wave band '//numch, flag_xyfill=.true.) + call addfld(sw_fa_name(i), (/'lev'/), 'A',' ', 'SW FA for wave band '//numch, flag_xyfill=.true.) + call add_default (sw_ta_name(i), 2, ' ') + call add_default (sw_wa_name(i), 2, ' ') + call add_default (sw_ga_name(i), 2, ' ') + call add_default (sw_fa_name(i), 2, ' ') + end do + do i = 1,nlwbands + write(numch,'(I2.2)') i + lw_ta_name(i) = 'TALW'//numch + call addfld(lw_ta_name(i), (/'lev'/), 'A',' ', 'LW TAU for wave band '//numch, flag_xyfill=.true.) + call add_default (lw_ta_name(i), 2, ' ') + end do - deallocate(aernames) + end if end subroutine aer_rad_props_init @@ -124,7 +131,6 @@ subroutine aer_rad_props_sw(list_idx, state, pbuf, nnite, idxnite, & ! Arguments integer, intent(in) :: list_idx ! index of the climate or a diagnostic list type(physics_state), intent(in), target :: state - type(physics_buffer_desc), pointer :: pbuf(:) integer, intent(in) :: nnite ! number of night columns integer, intent(in) :: idxnite(:) ! local column indices of night columns @@ -138,97 +144,31 @@ subroutine aer_rad_props_sw(list_idx, state, pbuf, nnite, idxnite, & integer :: ncol integer :: lchnk - integer :: k ! index integer :: troplev(pcols) - ! optical props for each aerosol - ! hygroscopic - real(r8), pointer :: h_ext(:,:) - real(r8), pointer :: h_ssa(:,:) - real(r8), pointer :: h_asm(:,:) - ! non-hygroscopic - real(r8), pointer :: n_ext(:) - real(r8), pointer :: n_ssa(:) - real(r8), pointer :: n_asm(:) - real(r8), pointer :: n_scat(:) - real(r8), pointer :: n_ascat(:) - ! radius-dependent - real(r8), pointer :: r_ext(:,:) ! radius-dependent mass-specific extinction - real(r8), pointer :: r_scat(:,:) - real(r8), pointer :: r_ascat(:,:) - real(r8), pointer :: r_mu(:) ! log(radius) domain variable for r_ext, r_scat, r_ascat - - ! radiative properties for each aerosol - real(r8) :: ta (pcols,pver,nswbands) - real(r8) :: tw (pcols,pver,nswbands) - real(r8) :: twf(pcols,pver,nswbands) - real(r8) :: twg(pcols,pver,nswbands) - - ! aerosol masses - real(r8), pointer :: aermmr(:,:) ! mass mixing ratio of aerosols - real(r8) :: mmr_to_mass(pcols,pver) ! conversion factor for mmr to mass - real(r8) :: aermass(pcols,pver) ! mass of aerosols - - ! for table lookup into rh grid - real(r8) :: es(pcols,pver) ! saturation vapor pressure - real(r8) :: qs(pcols,pver) ! saturation specific humidity - real(r8) :: rh(pcols,pver) - real(r8) :: rhtrunc(pcols,pver) - real(r8) :: wrh(pcols,pver) - integer :: krh(pcols,pver) - integer :: numaerosols ! number of bulk aerosols in climate/diagnostic list integer :: nmodes ! number of aerosol modes in climate/diagnostic list integer :: nbins ! number of aerosol bins in climate/diagnostic list - integer :: iaerosol ! index into bulk aerosol list - - character(len=ot_length) :: opticstype ! hygro or nonhygro - character(len=16) :: pbuf_fld + integer :: i !----------------------------------------------------------------------------- ncol = state%ncol lchnk = state%lchnk - ! compute mixing ratio to mass conversion - do k = 1, pver - mmr_to_mass(:ncol,k) = rga * state%pdeldry(:ncol,k) - enddo - - ! initialize to conditions that would cause failure - tau (:,:,:) = -100._r8 - tau_w (:,:,:) = -100._r8 - tau_w_g (:,:,:) = -100._r8 - tau_w_f (:,:,:) = -100._r8 - ! top layer (ilev = 0) has no aerosol (ie tau = 0) - ! also initialize rest of layers to accumulate od's + ! also initialize rest of layers to accumulate optical depths tau (1:ncol,:,:) = 0._r8 tau_w (1:ncol,:,:) = 0._r8 tau_w_g(1:ncol,:,:) = 0._r8 tau_w_f(1:ncol,:,:) = 0._r8 - ! calculate relative humidity for table lookup into rh grid - do k = 1, pver - call qsat(state%t(1:ncol,k), state%pmid(1:ncol,k), es(1:ncol,k), qs(1:ncol,k), ncol) - end do - rh(1:ncol,1:pver) = state%q(1:ncol,1:pver,1) / qs(1:ncol,1:pver) - - rhtrunc(1:ncol,1:pver) = min(rh(1:ncol,1:pver),1._r8) - krh(1:ncol,1:pver) = min(floor( rhtrunc(1:ncol,1:pver) * nrh ) + 1, nrh - 1) ! index into rh mesh - wrh(1:ncol,1:pver) = rhtrunc(1:ncol,1:pver) * nrh - krh(1:ncol,1:pver) ! (-) weighting on left side values - ! get number of bulk aerosols and number of modes in current list call rad_cnst_get_info(list_idx, naero=numaerosols, nmodes=nmodes, nbins=nbins) ! Contributions from modal and bin aerosols. - if (nmodes>0 .or. nbins>0) then + if (numaerosols>0 .or. nmodes>0 .or. nbins>0) then call aerosol_optics_cam_sw(list_idx, state, pbuf, nnite, idxnite, & tau, tau_w, tau_w_g, tau_w_f) - else - tau (1:ncol,:,:) = 0._r8 - tau_w (1:ncol,:,:) = 0._r8 - tau_w_g(1:ncol,:,:) = 0._r8 - tau_w_f(1:ncol,:,:) = 0._r8 end if !REMOVECAM - no longer need this when CAM is retired and pcols no longer exists @@ -236,92 +176,30 @@ subroutine aer_rad_props_sw(list_idx, state, pbuf, nnite, idxnite, & !REMOVECAM_END call tropopause_find_cam(state, troplev) - ! Contributions from bulk aerosols. - do iaerosol = 1, numaerosols - - ! get bulk aerosol mass mixing ratio - call rad_cnst_get_aer_mmr(list_idx, iaerosol, state, pbuf, aermmr) - aermass(1:ncol,1:top_lev-1) = 0._r8 - aermass(1:ncol,top_lev:pver) = aermmr(1:ncol,top_lev:pver) * mmr_to_mass(1:ncol,top_lev:pver) - - ! get optics type - call rad_cnst_get_aer_props(list_idx, iaerosol, opticstype=opticstype) - - select case (trim(opticstype)) - case('hygro','hygroscopic','hygroscopi') - ! get optical properties for hygroscopic aerosols - call rad_cnst_get_aer_props(list_idx, iaerosol, sw_hygro_ext=h_ext, sw_hygro_ssa=h_ssa, sw_hygro_asm=h_asm) - call get_hygro_rad_props(ncol, krh, wrh, aermass, h_ext, h_ssa, h_asm, ta, tw, twg, twf) - tau (1:ncol,1:pver,:) = tau (1:ncol,1:pver,:) + ta (1:ncol,:,:) - tau_w (1:ncol,1:pver,:) = tau_w (1:ncol,1:pver,:) + tw (1:ncol,:,:) - tau_w_g(1:ncol,1:pver,:) = tau_w_g(1:ncol,1:pver,:) + twg(1:ncol,:,:) - tau_w_f(1:ncol,1:pver,:) = tau_w_f(1:ncol,1:pver,:) + twf(1:ncol,:,:) - - case('nonhygro','insoluble ') - ! get optical properties for non-hygroscopic aerosols - call rad_cnst_get_aer_props(list_idx, iaerosol, sw_nonhygro_ext=n_ext, sw_nonhygro_ssa=n_ssa, & - sw_nonhygro_asm=n_asm) - - call get_nonhygro_rad_props(ncol, aermass, n_ext, n_ssa, n_asm, ta, tw, twg, twf) - tau (1:ncol,1:pver,:) = tau (1:ncol,1:pver,:) + ta (1:ncol,:,:) - tau_w (1:ncol,1:pver,:) = tau_w (1:ncol,1:pver,:) + tw (1:ncol,:,:) - tau_w_g(1:ncol,1:pver,:) = tau_w_g(1:ncol,1:pver,:) + twg(1:ncol,:,:) - tau_w_f(1:ncol,1:pver,:) = tau_w_f(1:ncol,1:pver,:) + twf(1:ncol,:,:) - - case('volcanic') - ! get optical properties for volcanic aerosols - call rad_cnst_get_aer_props(list_idx, iaerosol, sw_nonhygro_ext=n_ext, sw_nonhygro_scat=n_scat, & - sw_nonhygro_ascat=n_ascat) - - call get_volcanic_rad_props(ncol, aermass, n_ext, n_scat, n_ascat, ta, tw, twg, twf) - tau (1:ncol,1:pver,:) = tau (1:ncol,1:pver,:) + ta (1:ncol,:,:) - tau_w (1:ncol,1:pver,:) = tau_w (1:ncol,1:pver,:) + tw (1:ncol,:,:) - tau_w_g(1:ncol,1:pver,:) = tau_w_g(1:ncol,1:pver,:) + twg(1:ncol,:,:) - tau_w_f(1:ncol,1:pver,:) = tau_w_f(1:ncol,1:pver,:) + twf(1:ncol,:,:) - - case('volcanic_radius','volcanic_radius1','volcanic_radius2','volcanic_radius3') - pbuf_fld = 'VOLC_RAD_GEOM ' - if (len_trim(opticstype)>15) then - pbuf_fld = trim(pbuf_fld)//opticstype(16:16) - endif - ! get optical properties for volcanic aerosols - call rad_cnst_get_aer_props(list_idx, iaerosol, r_sw_ext=r_ext, r_sw_scat=r_scat, r_sw_ascat=r_ascat, mu=r_mu) - call get_volcanic_radius_rad_props(ncol, aermass, pbuf_fld, pbuf, r_ext, r_scat, r_ascat, r_mu, ta, tw, twg, twf) - tau (1:ncol,1:pver,:) = tau (1:ncol,1:pver,:) + ta (1:ncol,:,:) - tau_w (1:ncol,1:pver,:) = tau_w (1:ncol,1:pver,:) + tw (1:ncol,:,:) - tau_w_g(1:ncol,1:pver,:) = tau_w_g(1:ncol,1:pver,:) + twg(1:ncol,:,:) - tau_w_f(1:ncol,1:pver,:) = tau_w_f(1:ncol,1:pver,:) + twf(1:ncol,:,:) - - case('zero') - ! no effect of "zero" aerosols, so update nothing - case default - call endrun('aer_rad_props_sw: unsupported opticstype :'//trim(opticstype)//':') - end select - - ! diagnostic output of individual aerosol optical properties - ! currently implemented for climate list only - call aer_vis_diag_out(lchnk, ncol, nnite, idxnite, iaerosol, ta(:,:,idx_sw_diag), list_idx, troplev) - - enddo - ! diagnostic output of total aerosol optical properties ! currently implemented for climate list only call aer_vis_diag_out(lchnk, ncol, nnite, idxnite, 0, tau(:,:,idx_sw_diag), list_idx, troplev) + ! debug diags + if (bam_debug) then + do i = 1,nswbands + call outfld(sw_ta_name(i), tau(1:ncol,1:pver,i), ncol, lchnk ) + call outfld(sw_wa_name(i), tau_w(1:ncol,1:pver,i), ncol, lchnk ) + call outfld(sw_ga_name(i), tau_w_g(1:ncol,1:pver,i), ncol, lchnk ) + call outfld(sw_fa_name(i), tau_w_f(1:ncol,1:pver,i), ncol, lchnk ) + end do + end if + end subroutine aer_rad_props_sw !============================================================================== -subroutine aer_rad_props_lw(list_idx, state, pbuf, odap_aer) +subroutine aer_rad_props_lw(list_idx, state, pbuf, odap_aer) ! Purpose: Compute aerosol transmissions needed in absorptivity/ ! emissivity calculations - ! lw extinction is the same representation for all - ! species. If this changes, this routine will need to do something - ! similar to the sw with routines like get_hygro_lw_abs - - use physics_buffer, only : pbuf_get_field, pbuf_get_index, physics_buffer_desc + use physics_buffer, only : physics_buffer_desc ! Arguments integer, intent(in) :: list_idx ! index of the climate or a diagnostic list @@ -332,407 +210,32 @@ subroutine aer_rad_props_lw(list_idx, state, pbuf, odap_aer) ! Local variables - integer :: bnd_idx ! LW band index - integer :: i ! column index - integer :: k ! lev index - integer :: ncol ! number of columns integer :: numaerosols ! number of bulk aerosols in climate/diagnostic list integer :: nmodes ! number of aerosol modes in climate/diagnostic list integer :: nbins ! number of aerosol bins in climate/diagnostic list - integer :: iaerosol ! index into bulk aerosol list - character(len=ot_length) :: opticstype ! hygro or nonhygro - - ! optical props for each aerosol - real(r8), pointer :: lw_abs(:) - real(r8), pointer :: lw_hygro_abs(:,:) - real(r8), pointer :: geometric_radius(:,:) - - ! volcanic lookup table - real(r8), pointer :: r_lw_abs(:,:) ! radius dependent mass-specific absorption coefficient - real(r8), pointer :: r_mu(:) ! log(geometric_mean_radius) domain samples of r_lw_abs(:,:) - integer :: idx ! index to pbuf for geometric radius - real(r8) :: mu(pcols,pver) ! log(geometric_radius) - real(r8) :: r_mu_min, r_mu_max, wmu, mutrunc - integer :: nmu, kmu - - ! for table lookup into rh grid - real(r8) :: es(pcols,pver) ! saturation vapor pressure - real(r8) :: qs(pcols,pver) ! saturation specific humidity - real(r8) :: rh(pcols,pver) - real(r8) :: rhtrunc(pcols,pver) - real(r8) :: wrh(pcols,pver) - integer :: krh(pcols,pver) - - ! aerosol (vertical) mass path and extinction - ! aerosol masses - real(r8), pointer :: aermmr(:,:) ! mass mixing ratio of aerosols - real(r8) :: mmr_to_mass(pcols,pver) ! conversion factor for mmr to mass - real(r8) :: aermass(pcols,pver) ! mass of aerosols - - character(len=16) :: pbuf_fld + integer :: i, ncol !----------------------------------------------------------------------------- - ncol = state%ncol - ! get number of bulk aerosols and number of modes in current list call rad_cnst_get_info(list_idx, naero=numaerosols, nmodes=nmodes, nbins=nbins) + odap_aer = 0._r8 + ! Contributions from modal and sectional aerosols. - if (nmodes>0 .or. nbins>0) then + if (numaerosols>0 .or. nmodes>0 .or. nbins>0) then call aerosol_optics_cam_lw(list_idx, state, pbuf, odap_aer) - else - odap_aer = 0._r8 end if - ! Contributions from bulk aerosols. - if (numaerosols > 0) then - - ! compute mixing ratio to mass conversion - do k = 1, pver - mmr_to_mass(:ncol,k) = rga * state%pdeldry(:ncol,k) - end do - - ! calculate relative humidity for table lookup into rh grid - do k = 1, pver - call qsat(state%t(1:ncol,k), state%pmid(1:ncol,k), es(1:ncol,k), qs(1:ncol,k), ncol) + ! debug diags + if (bam_debug) then + ncol = state%ncol + do i = 1,nlwbands + call outfld(lw_ta_name(i), odap_aer(1:ncol,1:pver,i), ncol, state%lchnk) end do - rh(1:ncol,1:pver) = state%q(1:ncol,1:pver,1) / qs(1:ncol,1:pver) - - rhtrunc(1:ncol,1:pver) = min(rh(1:ncol,1:pver),1._r8) - krh(1:ncol,1:pver) = min(floor( rhtrunc(1:ncol,1:pver) * nrh ) + 1, nrh - 1) ! index into rh mesh - wrh(1:ncol,1:pver) = rhtrunc(1:ncol,1:pver) * nrh - krh(1:ncol,1:pver) ! (-) weighting on left side values - end if - ! Loop over bulk aerosols in list. - do iaerosol = 1, numaerosols - - ! get aerosol mass mixing ratio - call rad_cnst_get_aer_mmr(list_idx, iaerosol, state, pbuf, aermmr) - aermass(1:ncol,1:top_lev-1) = 0._r8 - aermass(1:ncol,top_lev:pver) = aermmr(1:ncol,top_lev:pver) * mmr_to_mass(1:ncol,top_lev:pver) - - ! get optics type - call rad_cnst_get_aer_props(list_idx, iaerosol, opticstype=opticstype) - select case (trim(opticstype)) - case('hygroscopic') - ! get optical properties for hygroscopic aerosols - call rad_cnst_get_aer_props(list_idx, iaerosol, lw_hygro_ext=lw_hygro_abs) - do bnd_idx = 1, nlwbands - do k = 1, pver - do i = 1, ncol - odap_aer(i, k, bnd_idx) = odap_aer(i, k, bnd_idx) + & - aermass(i, k) * & - ((1 + wrh(i,k)) * lw_hygro_abs(krh(i,k)+1,bnd_idx) & - - wrh(i,k) * lw_hygro_abs(krh(i,k), bnd_idx)) - end do - end do - end do - case('insoluble','nonhygro','hygro','volcanic') - ! get optical properties for hygroscopic aerosols - call rad_cnst_get_aer_props(list_idx, iaerosol, lw_ext=lw_abs) - do bnd_idx = 1, nlwbands - do k = 1, pver - do i = 1, ncol - odap_aer(i,k,bnd_idx) = odap_aer(i,k,bnd_idx) + lw_abs(bnd_idx)*aermass(i,k) - end do - end do - end do - - case('volcanic_radius','volcanic_radius1','volcanic_radius2','volcanic_radius3') - pbuf_fld = 'VOLC_RAD_GEOM ' - if (len_trim(opticstype)>15) then - pbuf_fld = trim(pbuf_fld)//opticstype(16:16) - endif - - ! get optical properties for hygroscopic aerosols - call rad_cnst_get_aer_props(list_idx, iaerosol, r_lw_abs=r_lw_abs, mu=r_mu) - ! get microphysical properties for volcanic aerosols - idx = pbuf_get_index(pbuf_fld) - call pbuf_get_field(pbuf, idx, geometric_radius ) - - ! interpolate in radius - ! caution: clip the table with no warning when outside bounds - nmu = size(r_mu) - r_mu_max = r_mu(nmu) - r_mu_min = r_mu(1) - do i = 1, ncol - do k = 1, pver - if(geometric_radius(i,k) > 0._r8) then - mu(i,k) = log(geometric_radius(i,k)) - else - mu(i,k) = 0._r8 - endif - mutrunc = max(min(mu(i,k),r_mu_max),r_mu_min) - kmu = max(min(1 + (mutrunc-r_mu_min)/(r_mu_max-r_mu_min)*(nmu-1),nmu-1._r8),1._r8) - wmu = max(min( (mutrunc -r_mu(kmu)) / (r_mu(kmu+1) - r_mu(kmu)) ,1._r8),0._r8) - do bnd_idx = 1, nlwbands - odap_aer(i,k,bnd_idx) = odap_aer(i,k,bnd_idx) + & - aermass(i,k) * & - ((1._r8 - wmu) * r_lw_abs(bnd_idx, kmu ) + & - (wmu) * r_lw_abs(bnd_idx, kmu+1)) - end do - end do - end do - - case('zero') - ! zero aerosols types have no optical effect, so do nothing. - case default - call endrun('aer_rad_props_lw: unsupported opticstype: '//trim(opticstype)) - end select - end do - end subroutine aer_rad_props_lw -!============================================================================== -! Private methods -!============================================================================== - -subroutine get_hygro_rad_props(ncol, krh, wrh, mass, ext, ssa, asm, & - tau, tau_w, tau_w_g, tau_w_f) - - ! Arguments - integer, intent(in) :: ncol - integer, intent(in) :: krh(pcols,pver) ! index for linear interpolation of optics on rh - real(r8), intent(in) :: wrh(pcols,pver) ! weight for linear interpolation of optics on rh - real(r8), intent(in) :: mass(pcols,pver) - real(r8), intent(in) :: ext(:,:) - real(r8), intent(in) :: ssa(:,:) - real(r8), intent(in) :: asm(:,:) - - real(r8), intent(out) :: tau (pcols,pver,nswbands) - real(r8), intent(out) :: tau_w (pcols,pver,nswbands) - real(r8), intent(out) :: tau_w_g(pcols,pver,nswbands) - real(r8), intent(out) :: tau_w_f(pcols,pver,nswbands) - - ! Local variables - real(r8) :: ext1, ssa1, asm1 - integer :: icol, ilev, iswband - !----------------------------------------------------------------------------- - - do iswband = 1, nswbands - do icol = 1, ncol - do ilev = 1, pver - ext1 = (1 + wrh(icol,ilev)) * ext(krh(icol,ilev)+1,iswband) & - - wrh(icol,ilev) * ext(krh(icol,ilev), iswband) - ssa1 = (1 + wrh(icol,ilev)) * ssa(krh(icol,ilev)+1,iswband) & - - wrh(icol,ilev) * ssa(krh(icol,ilev), iswband) - asm1 = (1 + wrh(icol,ilev)) * asm(krh(icol,ilev)+1,iswband) & - - wrh(icol,ilev) * asm(krh(icol,ilev), iswband) - - tau (icol, ilev, iswband) = mass(icol, ilev) * ext1 - tau_w (icol, ilev, iswband) = mass(icol, ilev) * ext1 * ssa1 - tau_w_g(icol, ilev, iswband) = mass(icol, ilev) * ext1 * ssa1 * asm1 - tau_w_f(icol, ilev, iswband) = mass(icol, ilev) * ext1 * ssa1 * asm1 * asm1 - enddo - enddo - enddo - -end subroutine get_hygro_rad_props - -!============================================================================== - -subroutine get_nonhygro_rad_props(ncol, mass, ext, ssa, asm, & - tau, tau_w, tau_w_g, tau_w_f) - - ! Arguments - integer, intent(in) :: ncol - real(r8), intent(in) :: mass(pcols, pver) - real(r8), intent(in) :: ext(:) - real(r8), intent(in) :: ssa(:) - real(r8), intent(in) :: asm(:) - - real(r8), intent(out) :: tau (pcols, pver, nswbands) - real(r8), intent(out) :: tau_w (pcols, pver, nswbands) - real(r8), intent(out) :: tau_w_g(pcols, pver, nswbands) - real(r8), intent(out) :: tau_w_f(pcols, pver, nswbands) - - ! Local variables - integer :: iswband - real(r8) :: ext1, ssa1, asm1 - !----------------------------------------------------------------------------- - - do iswband = 1, nswbands - ext1 = ext(iswband) - ssa1 = ssa(iswband) - asm1 = asm(iswband) - tau (1:ncol,1:pver,iswband) = mass(1:ncol,1:pver) * ext1 - tau_w (1:ncol,1:pver,iswband) = mass(1:ncol,1:pver) * ext1 * ssa1 - tau_w_g(1:ncol,1:pver,iswband) = mass(1:ncol,1:pver) * ext1 * ssa1 * asm1 - tau_w_f(1:ncol,1:pver,iswband) = mass(1:ncol,1:pver) * ext1 * ssa1 * asm1 * asm1 - enddo - -end subroutine get_nonhygro_rad_props - -!============================================================================== - -subroutine get_volcanic_radius_rad_props(ncol, mass, pbuf_radius_name, pbuf, r_ext, r_scat, r_ascat, r_mu, & - tau, tau_w, tau_w_g, tau_w_f) - - - use physics_buffer, only : pbuf_get_field, pbuf_get_index - - ! Arguments - integer, intent(in) :: ncol - real(r8), intent(in) :: mass(pcols, pver) - character(len=*) :: pbuf_radius_name - type(physics_buffer_desc), pointer :: pbuf(:) - real(r8), intent(in) :: r_ext(:,:) - real(r8), intent(in) :: r_scat(:,:) - real(r8), intent(in) :: r_ascat(:,:) - real(r8), intent(in) :: r_mu(:) ! log(radius) domain of mass-specific optics - - real(r8), intent(out) :: tau (pcols, pver, nswbands) - real(r8), intent(out) :: tau_w (pcols, pver, nswbands) - real(r8), intent(out) :: tau_w_g(pcols, pver, nswbands) - real(r8), intent(out) :: tau_w_f(pcols, pver, nswbands) - - ! Local variables - integer :: iswband - real(r8) :: g - - integer :: idx ! index to radius in physics buffer - real(r8), pointer :: geometric_radius(:,:) ! geometric mean radius of volcanic aerosol - real(r8) :: mu(pcols,pver) ! log(geometric mean radius of volcanic aerosol) - integer :: kmu, nmu - real(r8) :: wmu, mutrunc, r_mu_max, r_mu_min - - ! interpolated values from table - real(r8) :: ext(nswbands) - real(r8) :: scat(nswbands) - real(r8) :: ascat(nswbands) - - integer :: i, k ! column level iterator - !----------------------------------------------------------------------------- - - tau =0._r8 - tau_w =0._r8 - tau_w_g=0._r8 - tau_w_f=0._r8 - - ! get microphysical properties for volcanic aerosols - idx = pbuf_get_index(pbuf_radius_name) - call pbuf_get_field(pbuf, idx, geometric_radius ) - - ! interpolate in radius - ! caution: clip the table with no warning when outside bounds - nmu = size(r_mu) - r_mu_max = r_mu(nmu) - r_mu_min = r_mu(1) - do i = 1, ncol - do k = 1, pver - if(geometric_radius(i,k) > 0._r8) then - mu(i,k) = log(geometric_radius(i,k)) - else - mu(i,k) = 0._r8 - endif - mutrunc = max(min(mu(i,k),r_mu_max),r_mu_min) - kmu = max(min(1 + (mutrunc-r_mu_min)/(r_mu_max-r_mu_min)*(nmu-1),nmu-1._r8),1._r8) - wmu = max(min( (mutrunc -r_mu(kmu)) / (r_mu(kmu+1) - r_mu(kmu)) ,1._r8),0._r8) - do iswband = 1, nswbands - ext(iswband) = & - ((1._r8 - wmu) * r_ext(iswband, kmu ) + & - (wmu) * r_ext(iswband, kmu+1)) - scat(iswband) = & - ((1._r8 - wmu) * r_scat(iswband, kmu ) + & - (wmu) * r_scat(iswband, kmu+1)) - ascat(iswband) = & - ((1._r8 - wmu) * r_ascat(iswband, kmu ) + & - (wmu) * r_ascat(iswband, kmu+1)) - if (scat(iswband).gt.0._r8) then - g = ascat(iswband)/scat(iswband) - else - g=0._r8 - endif - tau (i,k,iswband) = mass(i,k) * ext(iswband) - tau_w (i,k,iswband) = mass(i,k) * scat(iswband) - tau_w_g(i,k,iswband) = mass(i,k) * ascat(iswband) - tau_w_f(i,k,iswband) = mass(i,k) * g * ascat(iswband) - end do - enddo - enddo - -end subroutine get_volcanic_radius_rad_props - -!============================================================================== - -subroutine get_volcanic_rad_props(ncol, mass, ext, scat, ascat, & - tau, tau_w, tau_w_g, tau_w_f) - - ! Arguments - integer, intent(in) :: ncol - real(r8), intent(in) :: mass(pcols, pver) - real(r8), intent(in) :: ext(:) - real(r8), intent(in) :: scat(:) - real(r8), intent(in) :: ascat(:) - - real(r8), intent(out) :: tau (pcols, pver, nswbands) - real(r8), intent(out) :: tau_w (pcols, pver, nswbands) - real(r8), intent(out) :: tau_w_g(pcols, pver, nswbands) - real(r8), intent(out) :: tau_w_f(pcols, pver, nswbands) - - ! Local variables - integer :: iswband - real(r8) :: g - !----------------------------------------------------------------------------- - - do iswband = 1, nswbands - if (scat(iswband).gt.0._r8) then - g = ascat(iswband)/scat(iswband) - else - g=0._r8 - endif - tau (1:ncol,1:pver,iswband) = mass(1:ncol,1:pver) * ext(iswband) - tau_w (1:ncol,1:pver,iswband) = mass(1:ncol,1:pver) * scat(iswband) - tau_w_g(1:ncol,1:pver,iswband) = mass(1:ncol,1:pver) * ascat(iswband) - tau_w_f(1:ncol,1:pver,iswband) = mass(1:ncol,1:pver) * g * ascat(iswband) - enddo - -end subroutine get_volcanic_rad_props - -!============================================================================== - -subroutine aer_vis_diag_out(lchnk, ncol, nnite, idxnite, iaer, tau, diag_idx, troplev) - - ! output aerosol optical depth for the visible band - - integer, intent(in) :: lchnk - integer, intent(in) :: ncol ! number of columns - integer, intent(in) :: nnite ! number of night columns - integer, intent(in) :: idxnite(:) ! local column indices of night columns - integer, intent(in) :: iaer ! aerosol index -- if 0 then tau is a total for all aerosols - real(r8), intent(in) :: tau(:,:) ! aerosol optical depth for the visible band - integer, intent(in) :: diag_idx ! identifies whether the aerosol optics - ! is for the climate calc or a diagnostic calc - integer, intent(in) :: troplev(:) ! tropopause level - - ! Local variables - integer :: i - real(r8) :: tmp(pcols), tmp2(pcols) - !----------------------------------------------------------------------------- - - ! currently only implemented for climate calc - if (diag_idx > 0) return - - ! compute total column aerosol optical depth - tmp(1:ncol) = sum(tau(1:ncol,:), 2) - ! use fillvalue to indicate night columns - do i = 1, nnite - tmp(idxnite(i)) = fillvalue - end do - - if (iaer > 0) then - call outfld(odv_names(iaer), tmp, pcols, lchnk) - else - call outfld('AEROD_v', tmp, pcols, lchnk) - do i = 1, ncol - tmp2(i) = sum(tau(i,:troplev(i))) - end do - call outfld('AODvstrt', tmp2, pcols, lchnk) - end if - -end subroutine aer_vis_diag_out - !============================================================================== end module aer_rad_props diff --git a/src/physics/cam/aer_vis_diag_mod.F90 b/src/physics/cam/aer_vis_diag_mod.F90 new file mode 100644 index 0000000000..da96377864 --- /dev/null +++ b/src/physics/cam/aer_vis_diag_mod.F90 @@ -0,0 +1,115 @@ +!------------------------------------------------------------------------------- +! output aerosol optical depths for the visible band +!------------------------------------------------------------------------------- +module aer_vis_diag_mod + use shr_kind_mod, only: r8 => shr_kind_r8 + use cam_history, only: fieldname_len, addfld, outfld, add_default, horiz_only, hist_fld_active + use cam_history_support, only : fillvalue + use rad_constituents, only: rad_cnst_get_info + use ppgrid, only: pcols, pver + use phys_control, only: phys_getopts + use cam_abortutils, only: endrun + + implicit none + + private + public :: aer_vis_diag_init + public :: aer_vis_diag_out + + integer :: numaerosols=0 + character(len=fieldname_len), pointer :: odv_names(:) ! outfld names for visible OD + +contains + + !============================================================================== + subroutine aer_vis_diag_init() + + integer :: i, astat + character(len=64), allocatable :: aernames(:) + logical :: history_aero_optics ! Output aerosol optics diagnostics + + ! number of bulk aerosols in climate list + call rad_cnst_get_info(0, naero=numaerosols) + + if (numaerosols<1) return + + ! get names of bulk aerosols + allocate(aernames(numaerosols),stat=astat) + if( astat/= 0 ) call endrun('aer_vis_diag_init: aernames allocate error') + + call rad_cnst_get_info(0, aernames=aernames) + + call phys_getopts( history_aero_optics_out = history_aero_optics ) + + ! diagnostic output for bulk aerosols + ! create outfld names for visible OD + allocate(odv_names(numaerosols),stat=astat) + if( astat/= 0 ) call endrun('aer_vis_diag_init: odv_names allocate error') + + do i = 1, numaerosols + odv_names(i) = 'ODV_'//trim(aernames(i)) + call addfld (odv_names(i), horiz_only, 'A', '1', & + trim(aernames(i))//' optical depth in visible band', flag_xyfill=.true.) + call add_default(odv_names(i), 1, ' ') + end do + + end subroutine aer_vis_diag_init + + !============================================================================== + subroutine aer_vis_diag_out(lchnk, ncol, nnite, idxnite, iaer, tau, diag_idx, troplev) + + ! output aerosol optical depth for the visible band + + integer, intent(in) :: lchnk + integer, intent(in) :: ncol ! number of columns + integer, intent(in) :: nnite ! number of night columns + integer, intent(in) :: idxnite(:) ! local column indices of night columns + integer, intent(in) :: iaer ! aerosol index -- if 0 then tau is a total for all aerosols + real(r8), intent(in) :: tau(:,:) ! aerosol optical depth for the visible band + integer, intent(in) :: diag_idx ! identifies whether the aerosol optics + ! is for the climate calc or a diagnostic calc + integer, intent(in) :: troplev(:) ! tropopause level + + ! Local variables + integer :: i + real(r8) :: tmp(pcols), tmp2(pcols) + logical :: do_calc + !----------------------------------------------------------------------------- + + ! currently only implemented for climate calc + if (diag_idx > 0) return + + do_calc = .false. + if (iaer > 0) then + do_calc = hist_fld_active(odv_names(iaer)) + else + do_calc = hist_fld_active('AEROD_v') + end if + + if (do_calc) then + ! compute total column aerosol optical depth + tmp(1:ncol) = sum(tau(1:ncol,:), 2) + ! use fillvalue to indicate night columns + do i = 1, nnite + tmp(idxnite(i)) = fillvalue + end do + end if + + if (iaer > 0) then + call outfld(odv_names(iaer), tmp, pcols, lchnk) + else + call outfld('AEROD_v', tmp, pcols, lchnk) + + if (hist_fld_active('AODvstrt')) then + do i = 1, ncol + tmp2(i) = sum(tau(i,:troplev(i))) + end do + call outfld('AODvstrt', tmp2, pcols, lchnk) + end if + end if + + end subroutine aer_vis_diag_out + + !============================================================================== + +end module aer_vis_diag_mod diff --git a/src/physics/cam/aerosol_optics_cam.F90 b/src/physics/cam/aerosol_optics_cam.F90 index 9ea02e5d9f..1c5b7b201f 100644 --- a/src/physics/cam/aerosol_optics_cam.F90 +++ b/src/physics/cam/aerosol_optics_cam.F90 @@ -4,9 +4,9 @@ module aerosol_optics_cam use cam_logfile, only: iulog use radconstants, only: nswbands, nlwbands, idx_sw_diag, idx_uv_diag, idx_nir_diag use radconstants, only: get_lw_spectral_boundaries - use phys_prop, only: ot_length + use phys_prop, only: ot_length, numrh=>nrh use physics_types,only: physics_state - use physics_buffer,only: physics_buffer_desc + use physics_buffer,only: physics_buffer_desc, pbuf_get_field, pbuf_get_index use ppgrid, only: pcols, pver use physconst, only: rga, rair use cam_abortutils, only: endrun @@ -15,21 +15,30 @@ module aerosol_optics_cam use cam_history, only: addfld, add_default, outfld, horiz_only, fieldname_len use cam_history_support, only: fillvalue + use ref_pres, only: clim_modal_aero_top_lev use tropopause, only : tropopause_findChemTrop use wv_saturation, only: qsat use aerosol_properties_mod, only: aerosol_properties use modal_aerosol_properties_mod, only: modal_aerosol_properties use carma_aerosol_properties_mod, only: carma_aerosol_properties + use bulk_aerosol_properties_mod, only: bulk_aerosol_properties use aerosol_state_mod, only: aerosol_state use modal_aerosol_state_mod,only: modal_aerosol_state use carma_aerosol_state_mod,only: carma_aerosol_state + use bulk_aerosol_state_mod, only: bulk_aerosol_state use aerosol_optics_mod, only: aerosol_optics use refractive_aerosol_optics_mod, only: refractive_aerosol_optics use hygrocoreshell_aerosol_optics_mod, only: hygrocoreshell_aerosol_optics use hygrowghtpct_aerosol_optics_mod, only: hygrowghtpct_aerosol_optics + use hygroscopic_aerosol_optics_mod, only: hygroscopic_aerosol_optics + use hygro_aerosol_optics_mod, only: hygro_aerosol_optics + use insoluble_aerosol_optics_mod, only: insoluble_aerosol_optics + use volcrad_aerosol_optics_mod, only: volcrad_aerosol_optics + + use aer_vis_diag_mod, only: aer_vis_diag_out implicit none @@ -59,6 +68,8 @@ module aerosol_optics_cam logical :: carma_active = .false. logical :: modal_active = .false. + logical :: bulk_active = .false. + integer :: num_aero_models = 0 integer :: lw10um_indx = -1 ! wavelength index corresponding to 10 microns real(r8), parameter :: lw10um = 10._r8 ! microns @@ -76,6 +87,8 @@ module aerosol_optics_cam type(out_name), allocatable :: aodbindn_fields(:) type(out_name), allocatable :: aoddustdn_fields(:) + integer :: top_lev = 1 + contains !=============================================================================== @@ -132,28 +145,39 @@ subroutine aerosol_optics_cam_init character(len=*), parameter :: prefix = 'aerosol_optics_cam_init: ' integer :: nmodes=0, nbins=0, iaermod, istat, ilist, i + integer :: nbulk_aerosols=0 logical :: call_list(0:n_diag) real(r8) :: lwavlen_lo(nlwbands), lwavlen_hi(nlwbands) - integer :: m, n + integer :: m, n, cnt character(len=fieldname_len) :: fldname character(len=128) :: lngname logical :: history_aero_optics ! output aerosol optics diagnostics logical :: history_amwg ! output the variables used by the AMWG diag package logical :: history_dust ! output dust diagnostics + logical :: prog_modal_aero ! prognostic modal aerosols present character(len=cl) :: locfile call phys_getopts(history_amwg_out = history_amwg, & history_aero_optics_out = history_aero_optics, & - history_dust_out = history_dust ) + history_dust_out = history_dust, & + prog_modal_aero_out = prog_modal_aero ) + + ! Limit modal aerosols with top_lev here. + if (prog_modal_aero) then + top_lev = clim_modal_aero_top_lev + else + top_lev = 1 + endif num_aero_models = 0 - call rad_cnst_get_info(0, nmodes=nmodes, nbins=nbins) + call rad_cnst_get_info(0, nmodes=nmodes, nbins=nbins, naero=nbulk_aerosols) modal_active = nmodes>0 carma_active = nbins>0 + bulk_active = nbulk_aerosols>0 ! count aerosol models if (modal_active) then @@ -162,6 +186,9 @@ subroutine aerosol_optics_cam_init if (carma_active) then num_aero_models = num_aero_models+1 end if + if (bulk_active) then + num_aero_models = num_aero_models+1 + end if if (num_aero_models>0) then allocate(aero_props(num_aero_models), stat=istat) @@ -180,9 +207,15 @@ subroutine aerosol_optics_cam_init iaermod = iaermod+1 aero_props(iaermod)%obj => carma_aerosol_properties() end if + if (bulk_active) then + iaermod = iaermod+1 + aero_props(iaermod)%obj => bulk_aerosol_properties() + end if if (water_refindex_file=='NONE') then - call endrun(prefix//'water_refindex_file must be specified') + if (modal_active .or. carma_active) then + call endrun(prefix//'water_refindex_file must be specified') + end if else call getfil(water_refindex_file, locfile) call read_water_refindex(locfile) @@ -302,6 +335,8 @@ subroutine aerosol_optics_cam_init call endrun(prefix//'array allocation error: aoddustdn_fields') end if + cnt=0 + do n = 1,num_aero_models allocate(burden_fields(n)%name(aero_props(n)%obj%nbins()), stat=istat) @@ -332,7 +367,9 @@ subroutine aerosol_optics_cam_init do m = 1, aero_props(n)%obj%nbins() - write(fldname,'(a,i2.2)') 'BURDEN', m + cnt = cnt+1 + + write(fldname,'(a,i2.2)') 'BURDEN', cnt burden_fields(n)%name(m) = fldname write(lngname,'(a,i2.2)') 'Aerosol burden bin ', m call addfld (fldname, horiz_only, 'A', 'kg/m2', lngname, flag_xyfill=.true.) @@ -348,7 +385,7 @@ subroutine aerosol_optics_cam_init call add_default (fldname, 1, ' ') end if - write(fldname,'(a,i2.2)') 'AODDUST', m + write(fldname,'(a,i2.2)') 'AODDUST', cnt aoddust_fields(n)%name(m) = fldname write(lngname,'(a,i2,a)') 'Aerosol optical depth, day only, 550 nm mode ',m,' from dust' call addfld (aoddust_fields(n)%name(m), horiz_only, 'A', ' ', lngname, flag_xyfill=.true.) @@ -356,7 +393,7 @@ subroutine aerosol_optics_cam_init call add_default (fldname, 1, ' ') end if - write(fldname,'(a,i2.2)') 'BURDENdn', m + write(fldname,'(a,i2.2)') 'BURDENdn', cnt burdendn_fields(n)%name(m) = fldname write(lngname,'(a,i2)') 'Aerosol burden, day night, bin ', m call addfld (burdendn_fields(n)%name(m), horiz_only, 'A', 'kg/m2', lngname, flag_xyfill=.true.) @@ -372,7 +409,7 @@ subroutine aerosol_optics_cam_init call add_default (fldname, 1, ' ') end if - write(fldname,'(a,i2.2)') 'AODdnDUST', m + write(fldname,'(a,i2.2)') 'AODdnDUST', cnt aoddustdn_fields(n)%name(m) = fldname write(lngname,'(a,i2,a)') 'Aerosol optical depth 550 nm, day night, bin ',m,' from dust' call addfld (aoddustdn_fields(n)%name(m), horiz_only, 'A', ' ', lngname, flag_xyfill=.true.) @@ -554,6 +591,7 @@ subroutine aerosol_optics_cam_sw(list_idx, state, pbuf, nnite, idxnite, tauxar, real(r8), intent(inout) :: ga(pcols,0:pver,nswbands) ! asymmetry factor real(r8), intent(inout) :: fa(pcols,0:pver,nswbands) ! forward scattered fraction + real(r8) :: taubam(pcols,0:pver) character(len=*), parameter :: prefix = 'aerosol_optics_cam_sw: ' integer :: ibin, nbins @@ -659,7 +697,11 @@ subroutine aerosol_optics_cam_sw(list_idx, state, pbuf, nnite, idxnite, tauxar, real(r8) :: ssavis(pcols) integer :: troplev(pcols) - integer :: i, k + integer :: bam_cnt + + real(r8), pointer :: geometric_radius(:,:) + integer :: idx ! index to pbuf for geometric radius + character(len=16) :: pbuf_fld nullify(aero_optics) @@ -681,6 +723,9 @@ subroutine aerosol_optics_cam_sw(list_idx, state, pbuf, nnite, idxnite, tauxar, absorb = 0._r8 aodtot = 0._r8 tauxar = 0._r8 + wa = 0._r8 + ga = 0._r8 + fa = 0._r8 extinct = 0._r8 extinctnir = 0._r8 extinctuv = 0._r8 @@ -722,6 +767,10 @@ subroutine aerosol_optics_cam_sw(list_idx, state, pbuf, nnite, idxnite, tauxar, iaermod = iaermod+1 aero_state(iaermod)%obj => carma_aerosol_state( state, pbuf ) end if + if (bulk_active) then + iaermod = iaermod+1 + aero_state(iaermod)%obj => bulk_aerosol_state( state, pbuf ) + end if allocate(pext(ncol), stat=istat) if (istat/=0) then @@ -740,6 +789,13 @@ subroutine aerosol_optics_cam_sw(list_idx, state, pbuf, nnite, idxnite, tauxar, call endrun(prefix//'array allocation error: pasm') end if + ! calculate relative humidity for table lookup into rh grid + call qsat(state%t(:ncol,:), state%pmid(:ncol,:), sate(:ncol,:), satq(:ncol,:), ncol, pver) + relh(:ncol,:) = state%q(1:ncol,:,1) / satq(:ncol,:) + relh(:ncol,:) = max(1.e-20_r8,relh(:ncol,:)) + + bam_cnt = 0 + aeromodel: do iaermod = 1,num_aero_models aeroprops => aero_props(iaermod)%obj @@ -755,6 +811,7 @@ subroutine aerosol_optics_cam_sw(list_idx, state, pbuf, nnite, idxnite, tauxar, dustaodbin(:) = 0._r8 burden(:) = 0._r8 aodbin(:) = 0.0_r8 + taubam(:,:) = 0._r8 call aeroprops%optics_params(list_idx, ibin, opticstype=opticstype) @@ -763,17 +820,38 @@ subroutine aerosol_optics_cam_sw(list_idx, state, pbuf, nnite, idxnite, tauxar, aero_optics=>refractive_aerosol_optics(aeroprops, aerostate, list_idx, ibin, & ncol, pver, nswbands, nlwbands, crefwsw, crefwlw) case('hygroscopic_coreshell') - ! calculate relative humidity for table lookup into rh grid - call qsat(state%t(:ncol,:), state%pmid(:ncol,:), sate(:ncol,:), satq(:ncol,:), ncol, pver) - relh(:ncol,:) = state%q(1:ncol,:,1) / satq(:ncol,:) - relh(:ncol,:) = max(1.e-20_r8,relh(:ncol,:)) aero_optics=>hygrocoreshell_aerosol_optics(aeroprops, aerostate, list_idx, & ibin, ncol, pver, relh(:ncol,:)) case('hygroscopic_wtp') aero_optics=>hygrowghtpct_aerosol_optics(aeroprops, aerostate, list_idx, & ibin, ncol, pver, sulfwtpct(:ncol,:)) + case('hygro') + ! Short-wave hygroscopic aerosol, Long-wave non-hygroscopic + ! aerosol optical properties + aero_optics=>hygro_aerosol_optics(aeroprops, aerostate, list_idx, & + ibin, ncol, pver, numrh, relh(:ncol,:)) + case('hygroscopic') + ! Short-wave and long-wave hygroscopic aerosol properties + aero_optics=>hygroscopic_aerosol_optics(aeroprops, aerostate, list_idx, & + ibin, ncol, pver, numrh, relh(:ncol,:)) + + case('nonhygro', 'insoluble') + aero_optics=>insoluble_aerosol_optics(aeroprops, aerostate, list_idx, ibin) + + case('volcanic_radius','volcanic_radius1','volcanic_radius2','volcanic_radius3') + pbuf_fld = 'VOLC_RAD_GEOM ' + if (len_trim(opticstype)>15) then + pbuf_fld = trim(pbuf_fld)//opticstype(16:16) + endif + ! get microphysical properties for volcanic aerosols + idx = pbuf_get_index(pbuf_fld) + call pbuf_get_field(pbuf, idx, geometric_radius ) + + aero_optics=>volcrad_aerosol_optics(aeroprops, aerostate, list_idx, & + ibin, ncol, pver, geometric_radius(:ncol,:)) + case default - call endrun(prefix//'optics method not recognized') + call endrun(prefix//'optics method not recognized: '//trim(opticstype)) end select if (associated(aero_optics)) then @@ -783,24 +861,31 @@ subroutine aerosol_optics_cam_sw(list_idx, state, pbuf, nnite, idxnite, tauxar, wavelength: do iwav = 1, nswbands - vertical: do ilev = 1, pver + vertical: do ilev = top_lev, pver call aero_optics%sw_props(ncol, ilev, iwav, pext, pabs, palb, pasm ) call init_diags column: do icol = 1,ncol - dopaer(icol) = pext(icol)*mass(icol,ilev) + dopaer(icol) = pext(icol) * mass(icol,ilev) tauxar(icol,ilev,iwav) = tauxar(icol,ilev,iwav) + dopaer(icol) wa(icol,ilev,iwav) = wa(icol,ilev,iwav) + dopaer(icol)*palb(icol) ga(icol,ilev,iwav) = ga(icol,ilev,iwav) + dopaer(icol)*palb(icol)*pasm(icol) fa(icol,ilev,iwav) = fa(icol,ilev,iwav) + dopaer(icol)*palb(icol)*pasm(icol)*pasm(icol) - call update_diags + if (.not.aeroprops%is_bulk()) then + call update_diags() + end if end do column + if (aeroprops%is_bulk().and.iwav==idx_sw_diag) then + taubam(:ncol,ilev) = dopaer(:ncol) + end if + end do vertical + end do wavelength else @@ -810,7 +895,13 @@ subroutine aerosol_optics_cam_sw(list_idx, state, pbuf, nnite, idxnite, tauxar, deallocate(aero_optics) nullify(aero_optics) - call output_bin_diags + if (aeroprops%is_bulk()) then + bam_cnt = bam_cnt+1 + call aer_vis_diag_out(lchnk, ncol, nnite, idxnite, bam_cnt, taubam, & + list_idx, troplev) + else + call output_bin_diags() + end if end do binloop end do aeromodel @@ -1190,6 +1281,10 @@ subroutine aerosol_optics_cam_lw(list_idx, state, pbuf, tauxar) character(len=32) :: opticstype integer :: iaermod + real(r8), pointer :: geometric_radius(:,:) + integer :: idx ! index to pbuf for geometric radius + character(len=16) :: pbuf_fld + real(r8) :: lwabs(pcols,pver) lwabs = 0._r8 tauxar = 0._r8 @@ -1210,6 +1305,10 @@ subroutine aerosol_optics_cam_lw(list_idx, state, pbuf, tauxar) iaermod = iaermod+1 aero_state(iaermod)%obj => carma_aerosol_state( state, pbuf ) end if + if (bulk_active) then + iaermod = iaermod+1 + aero_state(iaermod)%obj => bulk_aerosol_state( state, pbuf ) + end if ncol = state%ncol @@ -1220,6 +1319,11 @@ subroutine aerosol_optics_cam_lw(list_idx, state, pbuf, tauxar) call endrun(prefix//'array allocation error: pabs') end if + ! calculate relative humidity for table lookup into rh grid + call qsat(state%t(:ncol,:), state%pmid(:ncol,:), sate(:ncol,:), satq(:ncol,:), ncol, pver) + relh(:ncol,:) = state%q(1:ncol,:,1) / satq(:ncol,:) + relh(:ncol,:) = max(1.e-20_r8,relh(:ncol,:)) + aeromodel: do iaermod = 1,num_aero_models aeroprops => aero_props(iaermod)%obj @@ -1238,17 +1342,37 @@ subroutine aerosol_optics_cam_lw(list_idx, state, pbuf, tauxar) aero_optics=>refractive_aerosol_optics(aeroprops, aerostate, list_idx, ibin, & ncol, pver, nswbands, nlwbands, crefwsw, crefwlw) case('hygroscopic_coreshell') - ! calculate relative humidity for table lookup into rh grid - call qsat(state%t(:ncol,:), state%pmid(:ncol,:), sate(:ncol,:), satq(:ncol,:), ncol, pver) - relh(:ncol,:) = state%q(1:ncol,:,1) / satq(:ncol,:) - relh(:ncol,:) = max(1.e-20_r8,relh(:ncol,:)) aero_optics=>hygrocoreshell_aerosol_optics(aeroprops, aerostate, list_idx, & ibin, ncol, pver, relh(:ncol,:)) case('hygroscopic_wtp') aero_optics=>hygrowghtpct_aerosol_optics(aeroprops, aerostate, list_idx, & ibin, ncol, pver, sulfwtpct(:ncol,:)) + + case('hygroscopic') + aero_optics=>hygroscopic_aerosol_optics(aeroprops, aerostate, list_idx, ibin, & + ncol, pver, numrh, relh(:ncol,:)) + + case('hygro') + aero_optics=>hygro_aerosol_optics(aeroprops, aerostate, list_idx, ibin, & + ncol, pver, numrh, relh(:ncol,:)) + + case('nonhygro', 'insoluble') + aero_optics=>insoluble_aerosol_optics(aeroprops, aerostate, list_idx, ibin) + + case('volcanic_radius','volcanic_radius1','volcanic_radius2','volcanic_radius3') + pbuf_fld = 'VOLC_RAD_GEOM ' + if (len_trim(opticstype)>15) then + pbuf_fld = trim(pbuf_fld)//opticstype(16:16) + endif + ! get microphysical properties for volcanic aerosols + idx = pbuf_get_index(pbuf_fld) + call pbuf_get_field(pbuf, idx, geometric_radius ) + + aero_optics=>volcrad_aerosol_optics(aeroprops, aerostate, list_idx, & + ibin, ncol, pver, geometric_radius(:ncol,:)) + case default - call endrun(prefix//'optics method not recognized') + call endrun(prefix//'optics method not recognized: '//trim(opticstype)) end select if (associated(aero_optics)) then @@ -1259,7 +1383,7 @@ subroutine aerosol_optics_cam_lw(list_idx, state, pbuf, tauxar) call aero_optics%lw_props(ncol, ilev, iwav, pabs ) column: do icol = 1, ncol - dopaer(icol) = pabs(icol)*mass(icol,ilev) + dopaer(icol) = pabs(icol) * mass(icol,ilev) tauxar(icol,ilev,iwav) = tauxar(icol,ilev,iwav) + dopaer(icol) lwabs(icol,ilev) = lwabs(icol,ilev) + pabs(icol) end do column