Add 5-bin sea salt aerosol tracer#4713
Conversation
|
📖 Docs preview for this PR: https://clima.github.io/ClimaAtmos.jl/previews/PR4713/ |
| center_dust1_radius = 0.73, | ||
| center_dust2_radius = 1.4, | ||
| center_dust3_radius = 2.4, | ||
| center_dust4_radius = 4.5, | ||
| center_dust5_radius = 8, | ||
| center_ss1_radius = 0.55, | ||
| center_ss2_radius = 1.4, | ||
| center_ss3_radius = 2.4, | ||
| center_ss4_radius = 4.5, | ||
| center_ss5_radius = 8, | ||
| center_ss1_radius = aero_params.SSLT01_radius * FT(1e6), | ||
| center_ss2_radius = aero_params.SSLT02_radius * FT(1e6), | ||
| center_ss3_radius = aero_params.SSLT03_radius * FT(1e6), | ||
| center_ss4_radius = aero_params.SSLT04_radius * FT(1e6), | ||
| center_ss5_radius = aero_params.SSLT05_radius * FT(1e6), |
There was a problem hiding this comment.
Let's not change the radius here. It depends on how RRTMGP creates the lookup table for aerosol optics, and I need to double check it. But I think what we have before is correct.
There was a problem hiding this comment.
Changed dust radius based on this MERRA document (https://gmao.gsfc.nasa.gov/media/gmao_products/IqPeIBhIWLGrweILfVPj9Vmdoc/Dust_z72405udf3g1la8jwxmyspohi.pdf), given the rest of the values matched, but don't know how the RRTMGP lookup tables work, so can revert. Should these also be placed in ClimaParams?
There was a problem hiding this comment.
I think we can keep these values hard-coded and move them to parameters after finalizing aerosol optics. Where is e.g. SSLT01_radius defined now?
| precomputed = precomputed_quantities(Y, atmos) | ||
| tracers = tracer_cache( | ||
| Y, | ||
| prescribed_aerosol_names, |
There was a problem hiding this comment.
Do we want to support prescribed and prognostic aerosol at the same time? To allow having prognostic sea salt but prescribed everything else? This is appealing for testing but may result in complicated code. The alternative would be to allow for either fully prescribed or fully interactive.
| prescribed_aerosol_names, | ||
| time_varying_trace_gas_names, | ||
| start_date, | ||
| _interactive_aerosol_names(atmos.interactive_aerosols), |
There was a problem hiding this comment.
I would add it right after prescribed aerosol names, drop the underscore and name in prognostic aerosol names. "Interactive" is a new way to say prognostic and will be confusing if we start to use it interchangeably in the model.
There was a problem hiding this comment.
Can refactor, originally had as "prognostic" but chickened out after reading the Wikipedia page for "prognostic variable" since we're not integrating anything to get these
| prescribed_aerosols: | ||
| help: "Which aerosols to add. List of keys from the data file (e.g., CB1, CB2)." | ||
| value: [] | ||
| interactive_aerosols: |
| ORGW, | ||
| RAD, | ||
| TRAC, | ||
| OFRAC, |
There was a problem hiding this comment.
Not for this PR, but we should re-think this pattern of adding new fields to this struct whenever we add new functionality to the model
| start_date, | ||
| _interactive_aerosol_names(atmos.interactive_aerosols), | ||
| ) | ||
| ocean_fraction = ones(axes(Fields.level(Y.f, Fields.half))) # Populated via Interfacer.update_field! before first timestep |
There was a problem hiding this comment.
What to do in purely Atmos runs?
There was a problem hiding this comment.
Also, do we really have to cache it?
There was a problem hiding this comment.
Currently in purely Atmos runs I just treat everything as ocean. I presume somewhere there's an Artifact or something with an ocean mask strictly in Atmos but don't know where, will discuss.
Don't know if we need to cache since it only gets consumed by set_sea_salt_emission_flux!. That currently takes state and cache, is there another good way to pass it?
|
|
||
| set_precomputed_quantities!(Y, precomputing_arguments, FT(0)) | ||
|
|
||
| radiation_aerosol_names = ( |
There was a problem hiding this comment.
Is this used anywhere else? Maybe we can add both categories directly to radiation_args without specifying another entity
There was a problem hiding this comment.
It's not used outside creating a radiation cache, we could pass both--I think the right approach will depend on how we handle the prescribed vs prognostic aerosol model
| non_orographic_gravity_wave = non_orographic_gravity_wave_cache(Y, atmos) | ||
| orographic_gravity_wave = orographic_gravity_wave_cache(Y, atmos) | ||
| radiation = radiation_model_cache(Y, atmos, radiation_args...) | ||
| tracers = tracer_cache(Y, aerosol_names, time_varying_trace_gas_names, start_date) |
There was a problem hiding this comment.
We don't need the trace_gas_names here anyomre?
There was a problem hiding this comment.
Was not deleted, moved upwards so it can be placed in precomputing_arguments which I think was necessary for a previous EDMF wiring. I assume it can be moved back down with the other caches since passive tracer EDMF now runs through automatic inclusion, will experiment.
| SurfaceConditions.update_surface_conditions!(Y, p, FT(t)) | ||
| end | ||
|
|
||
| set_sea_salt_emission_flux!(Y, p) |
There was a problem hiding this comment.
Do we want to always apply it, or only under some user defined option?
| aod_sw_scattering = NaN, | ||
| # assuming fixed aerosol radius | ||
| center_dust1_radius = 0.55, | ||
| center_dust1_radius = 0.73, |
There was a problem hiding this comment.
Why this needs to be changed?
| center_ss2_radius = aero_params.SSLT02_radius * FT(1e6), | ||
| center_ss3_radius = aero_params.SSLT03_radius * FT(1e6), | ||
| center_ss4_radius = aero_params.SSLT04_radius * FT(1e6), | ||
| center_ss5_radius = aero_params.SSLT05_radius * FT(1e6), |
There was a problem hiding this comment.
We should not change those numbers by accident. They are consistent with RRTMGP assumptions, I think
| turbconv_model::PrognosticEDMFX, | ||
| microphysics_model, | ||
| chemistry_model, | ||
| interactive_aerosols, |
| # tracer_diagnostics.jl | ||
| "loadss", "mmrbcpi", "mmrbcpo", "mmrdust", "mmrocpi", "mmrocpo", | ||
| "mmrso4", "mmrss", "o3", | ||
| "mmrdst01", "mmrdst02", "mmrdst03", "mmrdst04", "mmrdst05", |
There was a problem hiding this comment.
Do we document anywhere what thise acronyms stand for?
| # at the cost of water vapor. | ||
| tracer_nonnegativity_vapor_tendency!(Yₜ, Y, p, t, microphysics_model) | ||
|
|
||
| sea_salt_deposition_tendency!(Yₜ, Y, p, t) |
There was a problem hiding this comment.
Should we dispatch on it based on the aerosol_model type?
| end | ||
| end | ||
|
|
||
| sea_salt_emission_tendency!(Yₜ, Y, p, t) |
There was a problem hiding this comment.
Consider dispatching based on aerosol model type?
| # assuming one updraft | ||
| ᶜχʲₜ = MatrixFields.get_field(Yₜ.c, get_χʲ_name_from_ρχ_name(ρχ_name)) | ||
| @. ᶜχʲₜ -= specific(btt, p.precomputed.ᶜρʲs.:(1)) | ||
| ᶜχʲ_name = get_χʲ_name_from_ρχ_name(ρχ_name) |
There was a problem hiding this comment.
Could you explain why this needs the additional if now?
There was a problem hiding this comment.
Yeah, I don't think it does need the if given we're sourcing
| function get_tracers(parsed_args) | ||
| aerosol_names = Tuple(parsed_args["prescribed_aerosols"]) | ||
| prescribed_aerosol_names = Tuple(parsed_args["prescribed_aerosols"]) | ||
| interactive_aerosol_names = Tuple(parsed_args["interactive_aerosols"]) |
| species is interactive; otherwise the prescribed (climatology) concentration | ||
| is used. | ||
| """ | ||
| function update_radiation_aerosol_concentrations!((; u, p, t)::I) where {I} |
There was a problem hiding this comment.
I think the naming is a bit confusing here. We have the update_prescribed_aerosol_concentrations! followed by update_radiation_aerosol_concentrations! To me this reads as prescribed aerosol do not affect radiation. There is a separate group of aerosol that do affect radiation that gets updated later.
There was a problem hiding this comment.
When it was just prescribed aerosols, this script both updated the prescribed aerosol fields forward in time, and then fed them to RRTMGP. To add prognostic, I split it into two functions, so we have separate functions that step froward in time, and then one unified function to compute the aerosol column density from both prescribed and prognostic sources and feed it to RRTMGP. I agree I should change the naming
|
|
||
| # BIN DEFNITIONS | ||
| const SEA_SALT_BIN_BOUNDS = ( | ||
| (0.03, 0.1), # μm, rep size 0.079 μm, SSLT01, |
There was a problem hiding this comment.
Are we ok with hard coding those, or could this be a parameter?
There was a problem hiding this comment.
Can convert to ClimaParam! Just wasn't sure what belongs in parameters.toml, and what data types to use (i.e. for a bin with a low_bound, high_bound, and representative size, should those be three separate entries or a tuple)
| UF.dimensionless_profile(uf_params, z_target, ζ_target, z₀, UF.MomentumTransport()) | ||
| F_anchor = | ||
| UF.dimensionless_profile(uf_params, z_anchor, ζ_anchor, z₀, UF.MomentumTransport()) | ||
| return max(u_anchor * F_target / F_anchor, FT(0)) |
There was a problem hiding this comment.
Is it possible for F_anchor to be zero?
| Precomputable radius-dependent part of the Gong (2003) integrand | ||
| """ | ||
| function _gong2003_r_integrand(r, theta) | ||
| A = 4.7 * (1 + theta * r)^(-0.017 * r^(-1.44)) |
There was a problem hiding this comment.
Those should not be hardcoded. We should put them in ClimaParams
There was a problem hiding this comment.
Yes, and it should be written in a way that is invariant under changes in units. That is, r -> r/r_gong_ref, with a reference scale r_gong_ref, with the exponents correspondingly adjusted. We do not want to take fractional powers of dimensional quantities.
There was a problem hiding this comment.
Also, here and throughout, the expressions promote to f64. Make type stable.
| return 1.373 * r^(-A) * (1 + 0.057 * r^3.45) * 10^(1.607 * exp(-B^2)) | ||
| end | ||
|
|
||
| function _gong_bin_integral(r_lo, r_hi, theta = 30.0, N = 512) |
There was a problem hiding this comment.
Theta a the shape parameter set in Gong, which controls the shape of the particle size distribution in the <0.2
\mu m size scale. The paper doesn't make an effort to fit it, given a lack of measurements at that size, and all implementations I've seen default to theta=30.
N is the number of quadrature points used in computing the integral.
| function sea_salt_emission_flux(u_10, T_sfc, bin_index; SST_adj = false) | ||
| FT = typeof(u_10) | ||
| bin_integral = FT(SEA_SALT_BIN_R_INTEGRALS[bin_index]) | ||
| number_flux = bin_integral * abs(u_10)^FT(3.41) |
There was a problem hiding this comment.
Those parameters should be part of Clima Params
There was a problem hiding this comment.
And, as before, never take a fractional power of a dimensional quantity. Introduce a reference scale and adjust constants to match results.
| ) | ||
|
|
||
| """ | ||
| sea_salt_emission_flux(u_10, T_sfc, bin_index) |
There was a problem hiding this comment.
Could you add a short entry to ClimaAtmos documentation?
|
|
||
| const SEA_SALT_BIN_R_INTEGRALS = ntuple( | ||
| i -> _gong_bin_integral(SEA_SALT_BIN_BOUNDS[i]...), | ||
| Val(5), |
There was a problem hiding this comment.
Do we ever want to support more than 5 bins?
|
|
||
| Called during `set_explicit_precomputed_quantities!'. `sea_salt_emission_tendency!` reads the cached values rather than recomputing. | ||
| """ | ||
| function set_sea_salt_emission_flux!(Y, p) |
There was a problem hiding this comment.
You could dispatch based on the aerosol_model and avoid the isempty. I think its a cleaner pattern. Thats what we do in most other tendencies
There was a problem hiding this comment.
As in having different AerosolModel structs, one for prescribed and one for prognostic, and dispatch based on the input struct?
| surface_fluxes_params = CAP.surface_fluxes_params(p.params) | ||
| uf_params = SFP.uf_params(surface_fluxes_params) | ||
| κ = SFP.von_karman_const(surface_fluxes_params) | ||
| ustar_p = parent(sfc_conditions.ustar) |
There was a problem hiding this comment.
Why do we need the parent here?
| roughness_spec = SF.COARE3RoughnessParams{FT}() | ||
|
|
||
| z0_eff = p.scratch.temp_field_level | ||
| z₀_p = parent(z0_eff) |
There was a problem hiding this comment.
I don't think the parent is needed here?
|
|
||
| z0_eff = p.scratch.temp_field_level | ||
| z₀_p = parent(z0_eff) | ||
| z₀_p .= SF.momentum_roughness.(roughness_spec, ustar_p, surface_fluxes_params, nothing) |
There was a problem hiding this comment.
Maybe
| z₀_p .= SF.momentum_roughness.(roughness_spec, ustar_p, surface_fluxes_params, nothing) | |
| @. z₀_p = SF.momentum_roughness(roughness_spec, ustar_p, surface_fluxes_params, nothing) |
| z₀_p = parent(z0_eff) | ||
| z₀_p .= SF.momentum_roughness.(roughness_spec, ustar_p, surface_fluxes_params, nothing) | ||
|
|
||
| z_c1_p = parent(Fields.level(Fields.coordinate_field(axes(Y.c)).z, 1)) |
| ) | ||
| mass_per_particle = ntuple( | ||
| i -> FT(4 / 3 * π * bin_radii[i]^3 * aero_params.seasalt_density), | ||
| Val(5), |
There was a problem hiding this comment.
I think its safer to make the number of bins a parameter. If we ever change it it'll be hard to trace all the places it is being hardcoded
| aero_params.SSLT05_radius, | ||
| ) | ||
| mass_per_particle = ntuple( | ||
| i -> FT(4 / 3 * π * bin_radii[i]^3 * aero_params.seasalt_density), |
There was a problem hiding this comment.
I suspect we don't need FT here
tapios
left a comment
tapios
left a comment
There was a problem hiding this comment.
Avoid magic numbers; make sure we only have nondimensional parameters (i.e., nondimensionalize the dimensional quantities appearing by introducing reference scales and and making them parameters).
| `u_anchor` at height `z_anchor` (m) via the MOST profile ratio. | ||
| TODO: This MOST computation assumes center-face layer 1 wind is is point value, not layer average | ||
| """ | ||
| function monin_obukhov_wind_extrapolated(z_target, z_anchor, u_anchor, L, uf_params, κ, z₀) |
There was a problem hiding this comment.
Why do you need this? Can't you just use 'compute_profile_value` to get the wind wherever you need it?
| Precomputable radius-dependent part of the Gong (2003) integrand | ||
| """ | ||
| function _gong2003_r_integrand(r, theta) | ||
| A = 4.7 * (1 + theta * r)^(-0.017 * r^(-1.44)) |
There was a problem hiding this comment.
Yes, and it should be written in a way that is invariant under changes in units. That is, r -> r/r_gong_ref, with a reference scale r_gong_ref, with the exponents correspondingly adjusted. We do not want to take fractional powers of dimensional quantities.
| """ | ||
| function _gong2003_r_integrand(r, theta) | ||
| A = 4.7 * (1 + theta * r)^(-0.017 * r^(-1.44)) | ||
| B = (0.433 - log10(r)) / 0.433 |
There was a problem hiding this comment.
Same here: use r/r_gong_ref and adjust constants to match result. Make type-stable.
| function _gong2003_r_integrand(r, theta) | ||
| A = 4.7 * (1 + theta * r)^(-0.017 * r^(-1.44)) | ||
| B = (0.433 - log10(r)) / 0.433 | ||
| return 1.373 * r^(-A) * (1 + 0.057 * r^3.45) * 10^(1.607 * exp(-B^2)) |
| Precomputable radius-dependent part of the Gong (2003) integrand | ||
| """ | ||
| function _gong2003_r_integrand(r, theta) | ||
| A = 4.7 * (1 + theta * r)^(-0.017 * r^(-1.44)) |
There was a problem hiding this comment.
Also, here and throughout, the expressions promote to f64. Make type stable.
| function sea_salt_emission_flux(u_10, T_sfc, bin_index; SST_adj = false) | ||
| FT = typeof(u_10) | ||
| bin_integral = FT(SEA_SALT_BIN_R_INTEGRALS[bin_index]) | ||
| number_flux = bin_integral * abs(u_10)^FT(3.41) |
There was a problem hiding this comment.
And, as before, never take a fractional power of a dimensional quantity. Introduce a reference scale and adjust constants to match results.
| number_flux = bin_integral * abs(u_10)^FT(3.41) | ||
| if SST_adj | ||
| SST_factor = | ||
| FT(0.3) + FT(0.1) * T_sfc - FT(0.0076) * T_sfc^2 + FT(0.00021) * T_sfc^3 |
There was a problem hiding this comment.
Same here: Introduce reference scales and make sure end results match by adjusting coefficients. Also double check: here temperatures are in K. Are they in the original reference?
| FT = eltype(Y) | ||
| (; sfc_conditions) = p.precomputed | ||
| surface_fluxes_params = CAP.surface_fluxes_params(p.params) | ||
| uf_params = SFP.uf_params(surface_fluxes_params) | ||
| κ = SFP.von_karman_const(surface_fluxes_params) | ||
| ustar_p = parent(sfc_conditions.ustar) | ||
| L_p = parent(sfc_conditions.obukhov_length) | ||
| roughness_spec = SF.COARE3RoughnessParams{FT}() | ||
|
|
||
| z0_eff = p.scratch.temp_field_level | ||
| z₀_p = parent(z0_eff) | ||
| z₀_p .= SF.momentum_roughness.(roughness_spec, ustar_p, surface_fluxes_params, nothing) | ||
|
|
||
| z_c1_p = parent(Fields.level(Fields.coordinate_field(axes(Y.c)).z, 1)) | ||
| u_z1_p = parent(norm.(Fields.level(Y.c.uₕ, 1))) | ||
|
|
||
| u_10 = p.scratch.ᶠtemp_field_level | ||
| parent(u_10) .= | ||
| monin_obukhov_wind_extrapolated.(FT(10), z_c1_p, u_z1_p, L_p, uf_params, κ, z₀_p) |
There was a problem hiding this comment.
You can just use 'compute_profile_value` from SurfaceFluxes.
| ) | ||
| mass_per_particle = ntuple( | ||
| i -> FT(4 / 3 * π * bin_radii[i]^3 * aero_params.seasalt_density), | ||
| Val(5), |
|
|
||
| FT = eltype(Y) | ||
|
|
||
| half_life = FT(0.55 * 86400) # AeroCOM SS residence is ~0.55 days |
4 participants
PR Dependency
New feature
ocean_fractionof cellsset_sea_salt_emission_flux!(aerosols/sea_salt.jl) duringset_explicit_precomputed_quantities!, using a 10 m wind extrapolated from the lowest model level via MOST, and applied as a bottom boundary flux bysea_salt_emission_tendency!insurface_flux_tendency!horizontal_tracer_advection_tendency!and vertical transport viavertical_transportinexplicit_vertical_advection_tendency!(both inprognostic_equations/advection.jl)sea_salt_deposition_tendency!,called from
additional_tendency!), a placeholder until weionaltracer_diagnostics.jlfor individual aerosol size bins (mmrsslt01–mmrsslt05, and analogousmmrdst0Xfor dust); new diagnostics for SSA emissions: total fluxemissand per-binemisslt01–emisslt05Config changes
interactive_aerosolsconfig option, takes list of (currently sea salt aerosol) i.e. ["SSLT01", ...], defaults empty. Stored onAtmosModelas aValof names so tracer sets are inferrable at compile time.Wiring changes
ClimaCoupler.jlviaIntcluded inAtmosSimulation` structtracer_cacheintracer_cache.jlinteractive_aerosol_sgs_variablesinprognostic_variables.jl), zero-initialized following wiring of ClimaAtmos #4629radiation/update_inputs.jl:update_aerosol_concentrations!split intoupdate_prescribed_aerosol_concentrations!(evaluates MERRA-2 time-varying inputs) andupdate_radiation_aerosol_concentrations!(packages per-layer column mass densities for RRTMGP). For each species, the interactive (prognostic) mass is used when active, falling back to the prescribed climatology otherwise, so interactive sea salt feeds aerosol radiation directly. The RRTMGP↔MERRA-2 name mapping is hoisted into a sharedRRTMGP_AEROSOL_NAMES_PAIRconstant.radiation/radiation.jl: sea-salt bin radii passed to RRTMGP now come from ClimaParams (MERRA2_seasalt_aerosol_bin0X_radius). Aerosol check inradiation_model_cacheis derived fromRRTMGP_AEROSOL_NAMES_PAIRrather than a duplicated name list, and now accepts interactive as well as prescribed aerosol names.Renames
AtmosSimulationkwargaerosol_names->prescribed_aerosol_names.Future Work