Skip to content

Commit

Permalink
step works
Browse files Browse the repository at this point in the history 
revamp model + unit tests

runs in AMIP driver

clean

step works

ice runs from driver with dss

name switch

delete redundant directory

add ocean params

add to CI pipeline

sea ice animation

format

unit test fixes

test for differentiate_turbulent_fluxes!

add artifacts to buildkite

dss q_sfc, solving tropical instability

shorten t_end for ci

format

fix test

dep fix

revs
  • Loading branch information
@LenkaNovak
LenkaNovak committed Sep 18, 2023
1 parent c435636 commit 85e1afe
Show file tree
Hide file tree
Showing 33 changed files with 966 additions and 3,445 deletions.
15 changes: 12 additions & 3 deletions .buildkite/pipeline.yml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -114,8 +114,12 @@ steps:
slurm_nodes: 3
slurm_tasks_per_node: 1

- label: "Component model tests"
command: "julia --color=yes --project=test/ test/component_model_tests.jl --run_name component_test --job_id component_test"
- label: "Component model tests: Slab sea ice"
command: "julia --color=yes --project=test/ test/component_model_tests/slab_seaice_tests.jl --run_name component_test_slab_seaice --job_id component_test_slab_seaice"
timeout_in_minutes: 5

- label: "Component model tests: Eisenman sea ice"
command: "julia --color=yes --project=test/ test/component_model_tests/eisenman_seaice_tests.jl --run_name component_test_eisenman_seaice --job_id component_test_eisenman_seaice"
timeout_in_minutes: 5

- group: "Integration Tests"
Expand Down Expand Up @@ -146,7 +150,6 @@ steps:
agents:
slurm_mem: 20GB

# Test: non-monotonous remapping for land mask
- label: "Slabplanet: non-monotonous surface remap"
key: "slabplanet_non-monotonous"
command: "julia --color=yes --project=experiments/AMIP/modular/ experiments/AMIP/modular/coupler_driver_modular.jl --run_name slabplanet_nonmono --FLOAT_TYPE Float64 --coupled true --surface_setup PrescribedSurface --moist equil --vert_diff true --rad gray --energy_check true --mode_name slabplanet --t_end 10days --dt_save_to_sol 9days --dt_cpl 200 --dt 200secs --mono_surface false --h_elem 4 --precip_model 0M --anim true --job_id slabplanet_nonmono"
Expand All @@ -161,6 +164,12 @@ steps:
agents:
slurm_mem: 20GB

- label: "Slabplanet: eisenman sea ice"
key: "slabplanet_eisenman"
command: "julia --color=yes --project=experiments/AMIP/modular/ experiments/AMIP/modular/coupler_driver_modular.jl --run_name slabplanet_eisenman --FLOAT_TYPE Float64 --coupled true --surface_setup PrescribedSurface --moist equil --vert_diff true --rad gray --energy_check true --mode_name slabplanet_eisenman --t_end 5days --dt_save_to_sol 9days --dt_cpl 200 --dt 200secs --mono_surface true --h_elem 4 --precip_model 0M --anim true --job_id slabplanet_eisenman"
artifact_paths: "experiments/AMIP/modular/output/slabplanet_eisenman/slabplanet_eisenman_artifacts/*"
agents:
slurm_mem: 20GB

# AMIP

Expand Down
296 changes: 296 additions & 0 deletions experiments/AMIP/modular/components/ocean/eisenman_seaice.jl
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,296 @@
import ClimaCoupler: FluxCalculator
import ClimaTimeSteppers as CTS
import DiffEqBase: ODEProblem, init, step!
import ClimaCoupler.FluxCalculator: update_turbulent_fluxes_point!, differentiate_turbulent_fluxes!
import ClimaCoupler.Interfacer: get_field, update_field!

# TODO: move to ClimaCore
import Base.getindex
Base.getindex(f::Fields.FieldVector, colidx) =
Fields.FieldVector(; zip(propertynames(f), map(x -> getproperty(f, x)[colidx], propertynames(f)))...)

"""
EisenmanIceSimulation{P, Y, D, I}
Thermodynamic 0-layer, based on the Semtner 1979 model and later refined by
Eisenmen 2009 and Zhang et al 2021.
"""
struct EisenmanIceSimulation{P, Y, D, I} <: SeaIceModelSimulation
params_ice::P
Y_init::Y
domain::D
integrator::I
end
name(::EisenmanIceSimulation) = "EisenmanIceSimulation"

Base.@kwdef struct EisenmanIceParameters{FT <: AbstractFloat}
z0m::FT = 1e-3 # roughness length for momentum [m]
z0b::FT = 1e-5 # roughness length for buoyancy [m]
C0_base::FT = 120 # ice base transfer coefficient [W / m2 / K]
T_base::FT = 273.16 # ice base temperature [K]
L_ice::FT = 3e8 # latent heat coefficient for ice [J / m3]
T_freeze::FT = 273.16 # temperature at freezing point [K]
k_ice::FT = 2 # thermal conductivity of ice [W / m / K]
α::FT = 0.70 # albedo
σ::FT = 5.67e-8 # Stefan Boltzmann constant [W / m^2 / K^4]
end

Base.@kwdef struct EisenmanOceanParameters{FT <: AbstractFloat}
h::FT = 1 # mixed layer depth [m]
ρ::FT = 1020 # density of the mixed layer [kg / m^3]
c::FT = 4000 # mixed layer heat capacity [J / kg / K ]
z0m::FT = 1e-3 # roughness length for momentum [m]
z0b::FT = 1e-5 # roughness length for buoyancy [m]
α::FT = 0.1 # albedo
end

"""
state_init(p::EisenmanIceParameters, space)
Initialize the state vectors for the Eisenman-Zhang sea ice model.
"""
function state_init(p::EisenmanIceParameters, space::Spaces.AbstractSpace)
Y = Fields.FieldVector(
T_sfc = ones(space) .* p.T_freeze,
h_ice = zeros(space),
T_ml = ones(space) .* 277,
q_sfc = ClimaCore.Fields.zeros(space),
)
Ya = Fields.FieldVector(
F_turb = ClimaCore.Fields.zeros(space),
∂F_turb_energy∂T_sfc = ClimaCore.Fields.zeros(space),
F_rad = ClimaCore.Fields.zeros(space),
F_base = ClimaCore.Fields.zeros(space),
ocean_qflux = ClimaCore.Fields.zeros(space),
ρ_sfc = ClimaCore.Fields.zeros(space),
)
return Y, Ya
end

"""
get_∂F_rad_energy∂T_sfc(T_sfc, p)
Calculate the derivative of the radiative flux with respect to the surface temperature.
"""
function get_∂F_rad_energy∂T_sfc(T_sfc, p)
FT = eltype(T_sfc)
@. FT(4) * (FT(1) - p.α) * p.σ * T_sfc^3
end

"""
solve_eisenman_model!(Y, Ya, p, Δt::FT)
Solve the Eisenman-Zhang sea ice model for one timestep.
"""
function solve_eisenman_model!(Y, Ya, p, thermo_params, Δt)

# model parameter sets
(; p_i, p_o) = p

# prognostic
(; T_sfc, h_ice, T_ml, q_sfc) = Y

# auxiliary
(; F_turb, ∂F_turb_energy∂T_sfc, F_rad, ocean_qflux, F_base) = Ya

# local
F_atm = @. F_turb + F_rad
∂F_atmo∂T_sfc = get_∂F_rad_energy∂T_sfc.(T_sfc, Ref(p_i)) .+ ∂F_turb_energy∂T_sfc

# ice thickness and mixed layer temperature changes due to atmosphereic and ocean fluxes
ice_covered = parent(h_ice)[1] > 0
if ice_covered # ice-covered
@. F_base = p_i.C0_base * (T_ml - p_i.T_base)
ΔT_ml = @. -(F_base - ocean_qflux) * Δt / (p_o.h * p_o.ρ * p_o.c)
Δh_ice = @. (F_atm - F_base - ocean_qflux) * Δt / p_i.L_ice
else # ice-free
ΔT_ml = @. -(F_atm - ocean_qflux) * Δt / (p_o.h * p_o.ρ * p_o.c)
Δh_ice = 0
end

# T_ml is not allowed to be below freezing
frazil_ice_formation = parent(T_ml .+ ΔT_ml)[1] < p_i.T_freeze
if frazil_ice_formation
# Note that ice formation, which requires T_ml < T_freeze, implies that Δh_ice increases
Δh_ice = @. Δh_ice - (T_ml + ΔT_ml - p_i.T_freeze) * (p_o.h * p_o.ρ * p_o.c) / p_i.L_ice
ΔT_ml = @. p_i.T_freeze - T_ml
end

# adjust ocean temperature if transition to ice-free
transition_to_icefree = (parent(h_ice)[1] > 0) & (parent(h_ice .+ Δh_ice)[1] <= 0)
if transition_to_icefree
ΔT_ml = @. ΔT_ml - (h_ice + Δh_ice) * p_i.L_ice / (p_o.h * p_o.ρ * p_o.c)
Δh_ice = @. -h_ice
end

# solve for T_sfc
remains_ice_covered = (parent(h_ice .+ Δh_ice)[1] > 0)
if remains_ice_covered
# if ice covered, solve implicity (for now one Newton iteration: ΔT_s = - F(T_s) / dF(T_s)/dT_s )
F_conductive = @. p_i.k_ice / (h_ice + Δh_ice) * (p_i.T_base - T_sfc)
ΔT_sfc = @. (-F_atm + F_conductive) / (p_i.k_ice / (h_ice + Δh_ice) + ∂F_atmo∂T_sfc)
surface_melting = (parent(T_sfc .+ ΔT_sfc)[1] > p_i.T_freeze)
if surface_melting
ΔT_sfc = @. p_i.T_freeze - T_sfc # NB: T_sfc not storing energy
end
# surface is ice-covered, so update T_sfc as ice surface temperature
T_sfc .+= ΔT_sfc
# update surface humidity
@. q_sfc = TD.q_vap_saturation_generic.(thermo_params, T_sfc, Ya.ρ_sfc, TD.Ice())
else # ice-free, so update T_sfc as mixed layer temperature
T_sfc .= T_ml .+ ΔT_ml
# update surface humidity
@. q_sfc = TD.q_vap_saturation_generic.(thermo_params, T_sfc, Ya.ρ_sfc, TD.Liquid())
end

Y.T_ml .+= ΔT_ml
Y.h_ice .+= Δh_ice
Y.T_sfc .= T_sfc
Y.q_sfc .= q_sfc

return Y, Ya
end

"""
∑tendencies(dY, Y, cache, _)
Calculate the tendencies for the Eisenman-Zhang sea ice model.
"""
function ∑tendencies(dY, Y, cache, _)
FT = eltype(dY)
Δt = cache.Δt
Ya = cache.Ya
thermo_params = cache.thermo_params
p = cache.params

@. dY.T_ml = -Y.T_ml / Δt
@. dY.h_ice = -Y.h_ice / Δt
@. dY.T_sfc = -Y.T_sfc / Δt
@. dY.q_sfc = -Y.q_sfc / Δt

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

# Get dY/dt
@. dY.T_ml += Y.T_ml / Δt
@. dY.h_ice += Y.h_ice / Δt
@. dY.T_sfc += Y.T_sfc / Δt
@. dY.q_sfc = -Y.q_sfc / Δt

# update ice area fraction (binary mask for now)
cache.ice_area_fraction .= Regridder.binary_mask.(Y.h_ice, threshold = eps())

end

"""
eisenman_seaice_init(::Type{FT}, tspan; space = nothing, area_fraction = nothing, thermo_params = nothing, stepper = CTS.RK4(), dt = 0.02, saveat = 1.0e10)
Initialize the Eisenman-Zhang sea ice model and simulation.
"""
function eisenman_seaice_init(
::Type{FT},
tspan;
space = nothing,
area_fraction = nothing,
thermo_params = nothing,
stepper = CTS.RK4(),
dt = 0.02,
saveat = 1.0e10,
) where {FT}

params_ice = EisenmanIceParameters{FT}()
params_ocean = EisenmanOceanParameters{FT}()
params = (; p_i = params_ice, p_o = params_ocean)

# initiate prognostic variables
Y, Ya = state_init(params_ice, space)

ode_algo = CTS.ExplicitAlgorithm(stepper)
ode_function = CTS.ClimaODEFunction(T_exp! = ∑tendencies, dss! = weighted_dss_slab!)

cache = (;
Ya = Ya,
Δt = dt,
params = params,
area_fraction = area_fraction,
ice_area_fraction = zeros(space),
thermo_params = thermo_params,
dss_buffer = ClimaCore.Spaces.create_dss_buffer(ClimaCore.Fields.zeros(space)),
)
problem = ODEProblem(ode_function, Y, FT.(tspan), cache)
integrator = init(problem, ode_algo, dt = FT(dt), saveat = FT(saveat), adaptive = false)

sim = EisenmanIceSimulation(params, Y, space, integrator)
@warn name(sim) *
" assumes gray radiation, no snow coverage, and PartitionedStateFluxes for the surface flux calculation."
return sim
end

# extensions required by Interfacer
get_field(sim::EisenmanIceSimulation, ::Val{:surface_temperature}) = sim.integrator.u.T_sfc
get_field(sim::EisenmanIceSimulation, ::Val{:surface_humidity}) = sim.integrator.u.q_sfc
get_field(sim::EisenmanIceSimulation, ::Val{:roughness_momentum}) =
@. sim.integrator.p.params.p_i.z0m * (sim.integrator.p.ice_area_fraction) +
sim.integrator.p.params.p_o.z0m .* (1 - sim.integrator.p.ice_area_fraction)
get_field(sim::EisenmanIceSimulation, ::Val{:roughness_buoyancy}) =
@. sim.integrator.p.params.p_i.z0b * (sim.integrator.p.ice_area_fraction) +
sim.integrator.p.params.p_o.z0b .* (1 - sim.integrator.p.ice_area_fraction)
get_field(sim::EisenmanIceSimulation, ::Val{:beta}) = convert(eltype(sim.integrator.u), 1.0)
get_field(sim::EisenmanIceSimulation, ::Val{:albedo}) =
@. sim.integrator.p.params.p_i.α * (sim.integrator.p.ice_area_fraction) +
sim.integrator.p.params.p_o.α .* (1 - sim.integrator.p.ice_area_fraction)
get_field(sim::EisenmanIceSimulation, ::Val{:area_fraction}) = sim.integrator.p.area_fraction
get_field(sim::EisenmanIceSimulation, ::Val{:air_density}) = sim.integrator.p.Ya.ρ_sfc

function update_field!(sim::EisenmanIceSimulation, ::Val{:area_fraction}, field::Fields.Field)
sim.integrator.p.area_fraction .= field
end
function update_field!(sim::EisenmanIceSimulation, ::Val{:turbulent_energy_flux}, field)
parent(sim.integrator.p.Ya.F_turb) .= parent(field)
end
function update_field!(sim::EisenmanIceSimulation, ::Val{:radiative_energy_flux}, field)
parent(sim.integrator.p.Ya.F_rad) .= parent(field)
end
function update_field!(sim::EisenmanIceSimulation, ::Val{:air_density}, field)
parent(sim.integrator.p.Ya.ρ_sfc) .= parent(field)
end

# extensions required by FieldExchanger
step!(sim::EisenmanIceSimulation, t) = step!(sim.integrator, t - sim.integrator.t, true)
reinit!(sim::EisenmanIceSimulation) = reinit!(sim.integrator)

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

"""
get_model_state_vector(sim::EisenmanIceSimulation)
Extension of Checkpointer.get_model_state_vector to get the model state.
"""
function get_model_state_vector(sim::EisenmanIceSimulation)
return sim.integrator.u
end

"""
differentiate_turbulent_fluxes!(sim::EisenmanIceSimulation, args)
Extension of differentiate_turbulent_fluxes! from FluxCalculator to get the turbulent fluxes.
"""
differentiate_turbulent_fluxes!(sim::EisenmanIceSimulation, args) =
differentiate_turbulent_fluxes!(sim::EisenmanIceSimulation, args..., ΔT_sfc = 0.1)

function update_field!(sim::EisenmanIceSimulation, ::Val{:∂F_turb_energy∂T_sfc}, field, colidx)
sim.integrator.p.Ya.∂F_turb_energy∂T_sfc[colidx] .= field
end

function get_total_energy(sim::EisenmanIceSimulation, _)
Y = sim.integrator.u
(; p_i, p_o) = sim.integrator.p.params
heat_ml = @. p_o.h * p_o.ρ * p_o.c * Y.T_ml
phase = @. p_i.L_ice .* sim.integrator.u.h_ice
return heat_ml .+ phase
end
59 changes: 59 additions & 0 deletions experiments/AMIP/modular/components/ocean/eisenman_seaice.md
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,59 @@
# Eisenman-Zhang Model
Thermodynamic 0-layer, based on the Semtner 1979 model and later refined by sea ice model as specified by
Eisenman & Wettlaufer 2009 and Zhang et al 2021.

There are three prognostic variables, the height of ice (`h_i`), the ocean mixed layer depth (`T_ml`) and the
surface air temperature (`T_s`).

## Formulation
$$
L_i \frac{dh_i}{dt} = F_{atm} - F_{base}
$$
with the latent heat of fusion, $L_i=3 \times 10^8$ J m$^{-3}$, and where the (upward-pointing) flux into the atmosphere is
$$
F_{atm} = F_{rad} + F_{SH} + F_{LH}
$$
where $F\_{rad}$ is the net radiative flux, $F_{SH}$ is the sensible heat flux and $F_{LH}$ is the latent heat flux

The basal flux is
$$
F_{base} = F_0(T_{ml} - T_{melt})
$$
where $T_{melt} = 273.16$ K is the freezing temperature, and the basal heat coefficient $F_0 = 120$ W m$^{-2}$ K$^{-1}$.

The surface temperature is solved implicitly using a balance between $F_{atm}(T_s)$ and the conductive heat flux through the ice slab, $F_{ice}$:
$$
F_{atm} = F_{ice} = k_i \frac{T_{melt} - T_s}{h_i}
$$
where $k_i = 2$ W m$^{-2}$ k$^{-1}$ is the thermal conductivity of ice.
Currently the implicit solve is implemented as one Newton iteration (sufficient for the current spatial and temporal resolution - see Semtner 1976):
$$
T_s^{t+1} = T_s + \frac{F}{dF /d T_s} = T_s^{t} + \frac{- F_{atm}^t + F_{ice}^{t+1}}{k_i/h_i^{t+1} + d F_{atm}^t / d T_s^t}
$$
where $h_i^{t+1}$ is the updated $h^i$ from the previous section, and $d F_{atm}^t / d T_s^t$ needs to be supplied from the atmosphere model (or crudely calculated in the coupler, given $T_s$, turbulent diffusivities and transfer coefficients, and atmos state).

### Warm surface
Where $T_s^{t+1} > T_{melt}$, we set $T_s^{t+1} = T_{melt}$.
The oceal temperature $T_{ml}$ uses the standard slab model representation in ice-free conditions:
$$
\rho_w c_w h_{ml}\frac{dT_{ml}}{dt} = - F_{atm}
$$
In this case of no ice $T_s^{t+1} = T_{ml}^{t+1}$, while ice-covered conditions require that:
$$
\rho_w c_w h_{ml}\frac{dT_{ml}}{dt} = - F_{base}
$$
with $T_{ml}^{t+1} = T_{melt}$ in the absence of oceanic lateral flux (Q-flux).

### Transitions between ice-covered and ice-free conditions
- If the updated $T_{ml}^{t+1} < T_{melt}$, set $T_{ml}^{t+1} = T_{melt}$ and grow ice ($h_i^{t+1}$) due to the corresponding energy deficit.
- If the updated $h_i^{t+1} <= 0$ from a non-zero $h_i^t$, adjust $h_i^{t+1} = 0$ and use the surplus energy to warm the mixed layer,
where ice is present the $T_{ml}$ = $T_{melt}$ and $F_{base}$ = 0 during the transition

## Potential extensions
- area `ice_area_fraction` adjustment (e.g., assuming a minimal thickness)
- add a skin temperature equation

# References
- [Semtner 1976](https://journals.ametsoc.org/view/journals/phoc/6/3/1520-0485_1976_006_0379_amfttg_2_0_co_2.xml)
- [Eisenman & Wettlaufer 2009](https://www.pnas.org/doi/full/10.1073/pnas.0806887106)
- [Zhang et al 2021](https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2021MS002671), which can be found on [GitHub](https://github.com/sally-xiyue/fms-idealized/blob/sea_ice_v1.0/exp/sea_ice/srcmods/mixed_layer.f90)
Loading

0 comments on commit 85e1afe

Please sign in to comment.