Skip to content

Commit

Permalink
rm bycolumn
Browse files Browse the repository at this point in the history
  • Loading branch information
juliasloan25 committed May 21, 2024
1 parent ab243e2 commit 6f571c3
Show file tree
Hide file tree
Showing 13 changed files with 123 additions and 238 deletions.
2 changes: 1 addition & 1 deletion docs/src/interfacer.md
Original file line number Diff line number Diff line change
Expand Up @@ -151,7 +151,7 @@ following properties:
| `turbulent_moisture_flux` | aerodynamic turbulent surface fluxes of energy (evaporation) | kg m^-2 s^-1 |

### SurfaceModelSimulation - optional functions
- `update_turbulent_fluxes_point!(::ComponentModelSimulation, fields::NamedTuple, colidx)`:
- `update_turbulent_fluxes_point!(::ComponentModelSimulation, fields::NamedTuple)`:
This function updates the turbulent fluxes of the component model simulation
at this point in horizontal space. The values are updated using the energy
and moisture turbulent fluxes stored in fields which are calculated by the
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -65,16 +65,9 @@ CAD.add_diagnostic_variable!(
(; C_E) = cache.atmos.vert_diff
interior_uₕ = CC.Fields.level(state.c.uₕ, 1)
ᶠp = ᶠK_E = cache.scratch.ᶠtemp_scalar
CC.Fields.bycolumn(axes(ᶜp)) do colidx
@. ᶠp[colidx] = CAD.ᶠinterp(ᶜp[colidx])
ᶜΔz_surface = CC.Fields.Δz_field(interior_uₕ)
@. ᶠK_E[colidx] = CA.eddy_diffusivity_coefficient(
C_E,
CA.norm(interior_uₕ[colidx]),
ᶜΔz_surface[colidx] / 2,
ᶠp[colidx],
)
end
@. ᶠp = CAD.ᶠinterp(ᶜp)
ᶜΔz_surface = CC.Fields.Δz_field(interior_uₕ)
@. ᶠK_E = CA.eddy_diffusivity_coefficient(C_E, CA.norm(interior_uₕ), ᶜΔz_surface / 2, ᶠp)
if isnothing(out)
return CAD.ᶜinterp.(ᶠK_E)
else
Expand Down
21 changes: 8 additions & 13 deletions experiments/ClimaEarth/components/land/climaland_bucket.jl
Original file line number Diff line number Diff line change
Expand Up @@ -29,8 +29,8 @@ Interfacer.name(::BucketSimulation) = "BucketSimulation"
"""
get_new_cache(p, Y, energy_check)
Returns a new `p` with an updated field to store e_per_area if energy conservation
checks are turned on.
Returns a new `p` with an updated field to store e_per_area if energy conservation
checks are turned on.
"""
function get_new_cache(p, Y, energy_check)
if energy_check
Expand Down Expand Up @@ -216,16 +216,12 @@ Interfacer.step!(sim::BucketSimulation, t) = Interfacer.step!(sim.integrator, t
Interfacer.reinit!(sim::BucketSimulation) = Interfacer.reinit!(sim.integrator)

# extensions required by FluxCalculator (partitioned fluxes)
function FluxCalculator.update_turbulent_fluxes_point!(
sim::BucketSimulation,
fields::NamedTuple,
colidx::CC.Fields.ColumnIndex,
)
function FluxCalculator.update_turbulent_fluxes_point!(sim::BucketSimulation, fields::NamedTuple)
(; F_turb_energy, F_turb_moisture) = fields
turbulent_fluxes = sim.integrator.p.bucket.turbulent_fluxes
turbulent_fluxes.shf[colidx] .= F_turb_energy
turbulent_fluxes.shf .= F_turb_energy
earth_params = sim.model.parameters.earth_param_set
turbulent_fluxes.vapor_flux[colidx] .= F_turb_moisture ./ LP.ρ_cloud_liq(earth_params)
turbulent_fluxes.vapor_flux .= F_turb_moisture ./ LP.ρ_cloud_liq(earth_params)
return nothing
end

Expand All @@ -234,15 +230,14 @@ function FluxCalculator.surface_thermo_state(
sim::BucketSimulation,
thermo_params::TD.Parameters.ThermodynamicsParameters,
thermo_state_int,
colidx::CC.Fields.ColumnIndex,
)

T_sfc = Interfacer.get_field(sim, Val(:surface_temperature), colidx)
T_sfc = Interfacer.get_field(sim, Val(:surface_temperature))
# Note that the surface air density, ρ_sfc, is computed using the atmospheric state at the first level and making ideal gas
# and hydrostatic balance assumptions. The land model does not compute the surface air density so this is
# a reasonable stand-in.
ρ_sfc = Interfacer.get_field(sim, Val(:air_density), colidx)
q_sfc = Interfacer.get_field(sim, Val(:surface_humidity), colidx) # already calculated in rhs! (cache)
ρ_sfc = Interfacer.get_field(sim, Val(:air_density))
q_sfc = Interfacer.get_field(sim, Val(:surface_humidity)) # already calculated in rhs! (cache)
@. TD.PhaseEquil_ρTq.(thermo_params, ρ_sfc, T_sfc, q_sfc)
end

Expand Down
23 changes: 8 additions & 15 deletions experiments/ClimaEarth/components/ocean/eisenman_seaice.jl
Original file line number Diff line number Diff line change
Expand Up @@ -140,8 +140,8 @@ end
function Interfacer.update_field!(sim::EisenmanIceSimulation, ::Val{:area_fraction}, field::CC.Fields.Field)
sim.integrator.p.area_fraction .= field
end
function Interfacer.update_field!(sim::EisenmanIceSimulation, ::Val{:∂F_turb_energy∂T_sfc}, field, colidx)
sim.integrator.p.Ya.∂F_turb_energy∂T_sfc[colidx] .= field
function Interfacer.update_field!(sim::EisenmanIceSimulation, ::Val{:∂F_turb_energy∂T_sfc}, field)
sim.integrator.p.Ya.∂F_turb_energy∂T_sfc .= field
end
function Interfacer.update_field!(sim::EisenmanIceSimulation, ::Val{:radiative_energy_flux_sfc}, field)
parent(sim.integrator.p.Ya.F_rad) .= parent(field)
Expand All @@ -155,13 +155,9 @@ Interfacer.step!(sim::EisenmanIceSimulation, t) = Interfacer.step!(sim.integrato
Interfacer.reinit!(sim::EisenmanIceSimulation) = Interfacer.reinit!(sim.integrator)

# extensions required by FluxCalculator (partitioned fluxes)
function FluxCalculator.update_turbulent_fluxes_point!(
sim::EisenmanIceSimulation,
fields::NamedTuple,
colidx::CC.Fields.ColumnIndex,
)
function FluxCalculator.update_turbulent_fluxes_point!(sim::EisenmanIceSimulation, fields::NamedTuple)
(; F_turb_energy) = fields
@. sim.integrator.p.Ya.F_turb[colidx] = F_turb_energy
@. sim.integrator.p.Ya.F_turb = F_turb_energy
end

"""
Expand Down Expand Up @@ -194,9 +190,8 @@ function FluxCalculator.differentiate_turbulent_fluxes!(
fluxes;
δT_sfc = 0.1,
)
(; thermo_state_int, surface_params, surface_scheme, colidx) = input_args
thermo_state_sfc_dT =
FluxCalculator.surface_thermo_state(sim, thermo_params, thermo_state_int, colidx, δT_sfc = δT_sfc)
(; thermo_state_int, surface_params, surface_scheme) = input_args
thermo_state_sfc_dT = FluxCalculator.surface_thermo_state(sim, thermo_params, thermo_state_int, δT_sfc = δT_sfc)
input_args = merge(input_args, (; thermo_state_sfc = thermo_state_sfc_dT))

# set inputs based on whether the surface_scheme is `MoninObukhovScheme` or `BulkScheme`
Expand All @@ -210,7 +205,7 @@ function FluxCalculator.differentiate_turbulent_fluxes!(
# calculate the derivative
∂F_turb_energy∂T_sfc = @. (F_shf_δT_sfc + F_lhf_δT_sfc - F_shf - F_lhf) / δT_sfc

Interfacer.update_field!(sim, Val(:∂F_turb_energy∂T_sfc), ∂F_turb_energy∂T_sfc, colidx)
Interfacer.update_field!(sim, Val(:∂F_turb_energy∂T_sfc), ∂F_turb_energy∂T_sfc)
end

###
Expand Down Expand Up @@ -354,9 +349,7 @@ function ∑tendencies(dY, Y, cache, _)
@. dY.T_sfc = -Y.T_sfc / Δt
@. dY.q_sfc = -Y.q_sfc / Δt

CC.Fields.bycolumn(axes(Y.T_sfc)) do colidx
solve_eisenman_model!(Y[colidx], Ya[colidx], p, thermo_params, Δt)
end
solve_eisenman_model!(Y, Ya, p, thermo_params, Δt)

# Get dY/dt
@. dY.T_ml += Y.T_ml / Δt
Expand Down
8 changes: 2 additions & 6 deletions experiments/ClimaEarth/components/ocean/prescr_seaice.jl
Original file line number Diff line number Diff line change
Expand Up @@ -131,13 +131,9 @@ Interfacer.step!(sim::PrescribedIceSimulation, t) = Interfacer.step!(sim.integra
Interfacer.reinit!(sim::PrescribedIceSimulation) = Interfacer.reinit!(sim.integrator)

# extensions required by FluxCalculator (partitioned fluxes)
function FluxCalculator.update_turbulent_fluxes_point!(
sim::PrescribedIceSimulation,
fields::NamedTuple,
colidx::CC.Fields.ColumnIndex,
)
function FluxCalculator.update_turbulent_fluxes_point!(sim::PrescribedIceSimulation, fields::NamedTuple)
(; F_turb_energy) = fields
@. sim.integrator.p.F_turb_energy[colidx] = F_turb_energy
@. sim.integrator.p.F_turb_energy = F_turb_energy
end

"""
Expand Down
8 changes: 2 additions & 6 deletions experiments/ClimaEarth/components/ocean/slab_ocean.jl
Original file line number Diff line number Diff line change
Expand Up @@ -143,13 +143,9 @@ Interfacer.step!(sim::SlabOceanSimulation, t) = Interfacer.step!(sim.integrator,
Interfacer.reinit!(sim::SlabOceanSimulation) = Interfacer.reinit!(sim.integrator)

# extensions required by FluxCalculator (partitioned fluxes)
function FluxCalculator.update_turbulent_fluxes_point!(
sim::SlabOceanSimulation,
fields::NamedTuple,
colidx::CC.Fields.ColumnIndex,
)
function FluxCalculator.update_turbulent_fluxes_point!(sim::SlabOceanSimulation, fields::NamedTuple)
(; F_turb_energy) = fields
@. sim.integrator.p.F_turb_energy[colidx] = F_turb_energy
@. sim.integrator.p.F_turb_energy = F_turb_energy
end

"""
Expand Down
156 changes: 71 additions & 85 deletions src/FluxCalculator.jl
Original file line number Diff line number Diff line change
Expand Up @@ -142,70 +142,65 @@ function partitioned_turbulent_fluxes!(
csf.F_turb_energy .*= FT(0)
csf.F_turb_moisture .*= FT(0)

# iterate over all columns (when regridding, this will need to change)
CC.Fields.bycolumn(boundary_space) do colidx
# atmos state of center level 1
z_int = Interfacer.get_field(atmos_sim, Val(:height_int), colidx)
uₕ_int = Interfacer.get_field(atmos_sim, Val(:uv_int), colidx)
thermo_state_int = Interfacer.get_field(atmos_sim, Val(:thermo_state_int), colidx)
z_sfc = Interfacer.get_field(atmos_sim, Val(:height_sfc), colidx)

for sim in model_sims
# iterate over all surface models with non-zero area fractions
if sim isa Interfacer.SurfaceModelSimulation
# get area fraction (min = 0, max = 1)
area_fraction = Interfacer.get_field(sim, Val(:area_fraction), colidx)
# get area mask [0, 1], where area_mask = 1 if area_fraction > 0
area_mask = Regridder.binary_mask.(area_fraction)

if !iszero(parent(area_mask))

thermo_state_sfc = surface_thermo_state(sim, thermo_params, thermo_state_int, colidx)

# set inputs based on whether the surface_scheme is `MoninObukhovScheme` or `BulkScheme`
surface_params = get_surface_params(atmos_sim)
scheme_properties = get_scheme_properties(surface_scheme, sim, colidx)
input_args = (;
thermo_state_sfc,
thermo_state_int,
uₕ_int,
z_int,
z_sfc,
scheme_properties,
surface_params,
surface_scheme,
colidx,
)
inputs = surface_inputs(surface_scheme, input_args)

# calculate the surface fluxes
fluxes = get_surface_fluxes_point!(inputs, surface_params)
(; F_turb_ρτxz, F_turb_ρτyz, F_shf, F_lhf, F_turb_moisture) = fluxes

# perform additional diagnostics if required
differentiate_turbulent_fluxes!(sim, (thermo_params, input_args, fluxes))

# update fluxes in the coupler
fields = (;
F_turb_ρτxz = F_turb_ρτxz,
F_turb_ρτyz = F_turb_ρτyz,
F_turb_energy = F_shf .+ F_lhf,
F_turb_moisture = F_turb_moisture,
)

# update the fluxes of this surface model
update_turbulent_fluxes_point!(sim, fields, colidx)

# add the flux contributing from this surface to the coupler field
@. csf.F_turb_ρτxz[colidx] += F_turb_ρτxz * area_fraction * area_mask
@. csf.F_turb_ρτyz[colidx] += F_turb_ρτyz * area_fraction * area_mask
@. csf.F_turb_energy[colidx] += (F_shf .+ F_lhf) * area_fraction * area_mask
@. csf.F_turb_moisture[colidx] += F_turb_moisture * area_fraction * area_mask

end
# atmos state of center level 1
z_int = Interfacer.get_field(atmos_sim, Val(:height_int))
uₕ_int = Interfacer.get_field(atmos_sim, Val(:uv_int))
thermo_state_int = Interfacer.get_field(atmos_sim, Val(:thermo_state_int))
z_sfc = Interfacer.get_field(atmos_sim, Val(:height_sfc))

for sim in model_sims
# iterate over all surface models with non-zero area fractions
if sim isa Interfacer.SurfaceModelSimulation
# get area fraction (min = 0, max = 1)
area_fraction = Interfacer.get_field(sim, Val(:area_fraction))
# get area mask [0, 1], where area_mask = 1 if area_fraction > 0
area_mask = Regridder.binary_mask.(area_fraction)

if !iszero(parent(area_mask))

thermo_state_sfc = surface_thermo_state(sim, thermo_params, thermo_state_int)

# set inputs based on whether the surface_scheme is `MoninObukhovScheme` or `BulkScheme`
surface_params = get_surface_params(atmos_sim)
scheme_properties = get_scheme_properties(surface_scheme, sim)
input_args = (;
thermo_state_sfc,
thermo_state_int,
uₕ_int,
z_int,
z_sfc,
scheme_properties,
surface_params,
surface_scheme,
)
inputs = surface_inputs(surface_scheme, input_args)

# calculate the surface fluxes
fluxes = get_surface_fluxes_point!(inputs, surface_params)
(; F_turb_ρτxz, F_turb_ρτyz, F_shf, F_lhf, F_turb_moisture) = fluxes

# perform additional diagnostics if required
differentiate_turbulent_fluxes!(sim, (thermo_params, input_args, fluxes))

# update fluxes in the coupler
fields = (;
F_turb_ρτxz = F_turb_ρτxz,
F_turb_ρτyz = F_turb_ρτyz,
F_turb_energy = F_shf .+ F_lhf,
F_turb_moisture = F_turb_moisture,
)

# update the fluxes of this surface model
update_turbulent_fluxes_point!(sim, fields)

# add the flux contributing from this surface to the coupler field
@. csf.F_turb_ρτxz += F_turb_ρτxz * area_fraction * area_mask
@. csf.F_turb_ρτyz += F_turb_ρτyz * area_fraction * area_mask
@. csf.F_turb_energy += (F_shf .+ F_lhf) * area_fraction * area_mask
@. csf.F_turb_moisture += F_turb_moisture * area_fraction * area_mask

end
end

end

# TODO: add allowable bounds here, check explicitly that all fluxes are equal
Expand All @@ -217,25 +212,21 @@ struct BulkScheme <: AbstractSurfaceFluxScheme end
struct MoninObukhovScheme <: AbstractSurfaceFluxScheme end

"""
get_scheme_properties(scheme::AbstractSurfaceFluxScheme, sim::Interfacer.SurfaceModelSimulation, colidx::CC.Fields.ColumnIndex)
get_scheme_properties(scheme::AbstractSurfaceFluxScheme, sim::Interfacer.SurfaceModelSimulation)
Returns the scheme-specific properties for the surface model simulation `sim`.
"""
function get_scheme_properties(::BulkScheme, sim::Interfacer.SurfaceModelSimulation, colidx::CC.Fields.ColumnIndex)
Ch = Interfacer.get_field(sim, Val(:heat_transfer_coefficient), colidx)
Cd = Interfacer.get_field(sim, Val(:beta), colidx)
beta = Interfacer.get_field(sim, Val(:drag_coefficient), colidx)
function get_scheme_properties(::BulkScheme, sim::Interfacer.SurfaceModelSimulation)
Ch = Interfacer.get_field(sim, Val(:heat_transfer_coefficient))
Cd = Interfacer.get_field(sim, Val(:beta))
beta = Interfacer.get_field(sim, Val(:drag_coefficient))
FT = eltype(Ch)
return (; z0b = FT(0), z0m = FT(0), Ch = Ch, Cd = Cd, beta = beta, gustiness = FT(1))
end
function get_scheme_properties(
::MoninObukhovScheme,
sim::Interfacer.SurfaceModelSimulation,
colidx::CC.Fields.ColumnIndex,
)
z0m = Interfacer.get_field(sim, Val(:roughness_momentum), colidx)
z0b = Interfacer.get_field(sim, Val(:roughness_buoyancy), colidx)
beta = Interfacer.get_field(sim, Val(:beta), colidx)
function get_scheme_properties(::MoninObukhovScheme, sim::Interfacer.SurfaceModelSimulation)
z0m = Interfacer.get_field(sim, Val(:roughness_momentum))
z0b = Interfacer.get_field(sim, Val(:roughness_buoyancy))
beta = Interfacer.get_field(sim, Val(:beta))
FT = eltype(z0m)
return (; z0b = z0b, z0m = z0m, Ch = FT(0), Cd = FT(0), beta = beta, gustiness = FT(1))
end
Expand Down Expand Up @@ -287,21 +278,20 @@ function surface_inputs(::MoninObukhovScheme, input_args::NamedTuple)
end

"""
surface_thermo_state(sim::Interfacer.SurfaceModelSimulation, thermo_params::TD.Parameters.ThermodynamicsParameters, thermo_state_int, colidx::CC.Fields.ColumnIndex)
surface_thermo_state(sim::Interfacer.SurfaceModelSimulation, thermo_params::TD.Parameters.ThermodynamicsParameters, thermo_state_int)
Returns the surface parameters for the surface model simulation `sim`. The default is assuming saturated surfaces, unless an extension is defined for the given `SurfaceModelSimulation`.
"""
function surface_thermo_state(
sim::Interfacer.SurfaceModelSimulation,
thermo_params::TD.Parameters.ThermodynamicsParameters,
thermo_state_int,
colidx::CC.Fields.ColumnIndex;
thermo_state_int;
δT_sfc = 0,
)
FT = eltype(parent(thermo_state_int))
@warn("Simulation " * Interfacer.name(sim) * " uses the default thermo (saturated) surface state", maxlog = 10)
# get surface temperature (or perturbed surface temperature for differentiation)
T_sfc = Interfacer.get_field(sim, Val(:surface_temperature), colidx) .+ FT(δT_sfc)
T_sfc = Interfacer.get_field(sim, Val(:surface_temperature)) .+ FT(δT_sfc)
ρ_sfc = extrapolate_ρ_to_sfc.(thermo_params, thermo_state_int, T_sfc) # ideally the # calculate elsewhere, here just getter...
q_sfc = TD.q_vap_saturation_generic.(thermo_params, T_sfc, ρ_sfc, TD.Liquid()) # default: saturated liquid surface
@. TD.PhaseEquil_ρTq.(thermo_params, ρ_sfc, T_sfc, q_sfc)
Expand Down Expand Up @@ -367,15 +357,11 @@ function get_surface_params(atmos_sim::Interfacer.AtmosModelSimulation)
end

"""
update_turbulent_fluxes_point!(sim::Interfacer.SurfaceModelSimulation, fields::NamedTuple, colidx::CC.Fields.ColumnIndex)
update_turbulent_fluxes_point!(sim::Interfacer.SurfaceModelSimulation, fields::NamedTuple)
Updates the fluxes in the surface model simulation `sim` with the fluxes in `fields`.
"""
function update_turbulent_fluxes_point!(
sim::Interfacer.SurfaceModelSimulation,
fields::NamedTuple,
colidx::CC.Fields.ColumnIndex,
)
function update_turbulent_fluxes_point!(sim::Interfacer.SurfaceModelSimulation, fields::NamedTuple)
return error(
"update_turbulent_fluxes_point! is required to be dispatched on" *
Interfacer.name(sim) *
Expand Down
13 changes: 0 additions & 13 deletions src/Interfacer.jl
Original file line number Diff line number Diff line change
Expand Up @@ -168,19 +168,6 @@ get_field(sim::ComponentModelSimulation, val::Val) = get_field_error(sim, val)

get_field_error(sim, val::Val{X}) where {X} = error("undefined field `$X` for " * name(sim))

"""
get_field(::ComponentModelSimulation, ::Val, colidx::CC.Fields.ColumnIndex)
Extension of `get_field(::ComponentModelSimulation, ::Val)`, indexing into the specified colum index.
"""
function get_field(sim::ComponentModelSimulation, val::Val, colidx::CC.Fields.ColumnIndex)
if get_field(sim, val) isa AbstractFloat
get_field(sim, val)
else
get_field(sim, val)[colidx]
end
end

"""
update_field!(::AtmosModelSimulation, ::Val, _...)
Expand Down
Loading

0 comments on commit 6f571c3

Please sign in to comment.