Merge branch 'ln/bring-back-semtner' of https://github.com/CliMA/Clim…

…aCoupler.jl into ln/bring-back-semtner

experiments/AMIP/modular/components/ocean/README.md

@@ -0,0 +1,72 @@
+# **Simple Sea Ice Model**
+
+# Semtner's Zero-layer Model
+
+Sea ice is approximated as a model that does not occupy a particular gridpoint between the mixed layer (ocean) and the atmosphere. It is essentially implemented in each domain column as an ODE, with no horizontal transport between the columns. In the absence of ice, the sea-ice model reduces to the slab ocean formulation. The ice is assumed to have a negligible heat capacity (so there is no energy storage due to internal temperature changes of the ice). The only storage changes arise from the ice thickness changes, which result from temperature differences at the ice surface or at the ice base.
+
+We followed the FMS implementation as in [Zhang et al 22](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) and is itself a modification of the 0-layer model of [Semtner 1976](https://journals.ametsoc.org/view/journals/phoc/6/3/1520-0485_1976_006_0379_amfttg_2_0_co_2.xml) (appendix). Chronologically, the algorithm follows these steps, with all fluxes defined positive upward:
+
+## 1. Ice thickness, $h_i$
+$$
+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} \approx \lambda (-{T_{sfc}} - T_{atm})
+$$
+- where the latter approximation was used for testing in earlier prototypes. $F_{atm}$ will be obtained from the atmospheric model (via the coupler). The flux at the ice base from the mixed layer 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}$.
+
+## 2. Ocean mixed layer temperature, $T_{ml}$
+- $T_{ml}$ is the standard slab ocean formulation in ice-free conditions:
+$$
+\rho_w c_w h_{ml}\frac{dT_{ml}}{dt} = - F_{atm}
+$$
+- while ice-covered conditions require that:
+$$
+\rho_w c_w h_{ml}\frac{dT_{ml}}{dt} = - F_{base}
+$$
+
+## 3. Transitions between ice free and ice covered 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.
+
+## 4. Surface temperature ($T_s$)
+- $T_s$ is determined 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:
+$$
+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).
+- Where $T_s^{t+1} > T_{melt}$, we set $T_s^{t+1} = T_{melt}$. Where there is no ice $T_s^{t+1} = T_{ml}^{t+1}$.
+
+## 5. Update ice mask
+- `mask = 1` if ice and `mask = 0` if no ice
+
+# Q flux (optional)
+- We can add an additional flux to the RHS of the $T_{ml}$ equations in (2), which corresponds to a more realistic ocean heating, coarsely mimicking otherwise neglected ocean dynamics, such as lateral advection, convection and diffusion. This is especially needed to improve low latitude oceanic forcing.
+
+- An analytic formulation can be written as:
+$$
+Q = Q_0(1-2\phi^2/w_\phi^2) \frac{exp(- (\phi^2/w_\phi^2))}{cos(\phi)}
+$$
+- where $\phi$ is latitude in radians, $Q_0$ is the amplitude of the equatorial heating and $w_\phi$ the width of the heating in radians.
+
+# Alternatives
+
+## [Semtner 1976](https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2021MS002671) 3-layer model
+- Semtner (1976) developed a simple model for the evolution ice (2 layers) and snow (1 layer) temperatures, which is represented by a 1D diffusve process also accounting for radiation, melting and energy release from brine poskets and accumulating snow.
+- Main formation of sea ice - energy fluxes from vertical boundaries and heat storage in brine pockets. The vertical processes retain the central role.
+- Semtner 1976 found that the 0-layer model is broadly comparable in terms of accuracy as the 3-layer model
+
+## Ocean-based dynamical sea-ice model
+- to be coordinated after AMIP
+
+

experiments/AMIP/modular/components/ocean/semtner_seaice.jl

@@ -0,0 +1,290 @@
+# slab_ice
+
+"""
+    SemtnerIceSimulation{P, Y, D, I}
+
+Thermodynamic 0-layer sea ice model as specified by Semtner 1979  Ice concentration is prescribed, and we solve the following energy equation:
+
+    (h * ρ * c) d T_sfc dt = -(F_turb_energy + F_radiativead) + F_conductive
+
+    with
+    F_conductive = k_ice (T_base - T_sfc) / (h)
+
+    The surface temperature (`T_sfc`) is the prognostic variable which is being
+    modified by turbulent aerodynamic (`F_turb_energy`) and radiative (`F_turb_energy`) fluxes,
+    as well as a conductive flux that depends on the temperature difference
+    across the ice layer (with `T_base` being prescribed).
+
+In the current version, the sea ice has a prescribed thickness, and we assume that it is not
+sublimating. That contribution has been zeroed out in the atmos fluxes.
+
+"""
+struct SemtnerIceSimulation{P, Y, D, I} <: SeaIceModelSimulation
+    params::P
+    Y_init::Y
+    domain::D
+    integrator::I
+end
+name(::SemtnerIceSimulation) = "SemtnerIceSimulation"
+
+struct SemtnerParameters{FT <: AbstractFloat}
+    h::FT # sea ice height
+    ρ::FT
+    c::FT
+    T_init::FT
+    z0m::FT
+    z0b::FT
+    ρc_ml::FT    # density times heat transfer coefficient for mixed layer [J / m2 / K ]
+    F0_base::FT  # ice base transfer coefficient [W / m2 / K]
+    T_base::FT   # ice base temperature [K]
+    L_ice::FT    # latent heat coefficient for ice [J / m3]
+    h_ml::FT     # mixed layer depth [m]
+    T_freeze::FT # temperature at freezing point [K]
+    k_ice::FT   # thermal conductivity of ice [W / m / K]
+    α::FT # albedo
+    σ::FT # Stefan Boltzmann constant
+end
+
+# init simulation
+function slab_ice_space_init(::Type{FT}, space, p) where {FT}
+    Y = Fields.FieldVector(T_sfc = ones(space) .* p.T_freeze, h_ice = zeros(space), T_ml = ones(space) .* p.T_freeze)
+    return Y
+end
+
+"""
+    solve_ice!(dT_sfc, T_sfc, (parameters, F_accumulated), t)
+
+slab RHS with an implicit solve ice and explicit (forward Euler) solve for ocean
+
+"""
+function solve_ice!(integ, Δt)
+
+    Y = integ.u
+    Ya = integ.p.Ya
+    p = integ.p.Ya.params
+
+    # prognostic
+    T_ml = Y.T_ml
+    T_sfc = Y.T_sfc
+    h_ice = Y.h_ice
+
+    # auxiliary
+    F_atm = @. Ya.F_aero + Ya.F_rad
+    ∂F_atm∂T_sfc = get_∂F_atm∂T_sfc(p, T_sfc, Ya) # this will be passed from atmos/SF.jl
+    F_base = similar(T_sfc)
+
+    sm = semtner_zero_layer_model.(T_ml, T_sfc, h_ice, F_atm, ∂F_atm∂T_sfc, Ref(p), Ref(Δt))
+
+    T_ml .= sm.T_ml
+    h_ice .= sm.h_ice
+    T_sfc .= sm.T_sfc
+    F_base .= sm.F_base
+
+    # update state
+    @. Y.T_ml = T_ml
+    @. Y.h_ice = h_ice
+    @. Y.T_sfc = T_sfc
+
+    Ya.ice_mask .= get_ice_mask.(h_ice)
+    Ya.F_base .= F_base
+
+    return nothing
+end
+
+
+function semtner_zero_layer_model(T_ml, T_sfc, h_ice, F_atm, ∂F_atm∂T_sfc, p, Δt)
+
+    # local
+    ocean_qflux = FT(0)
+    F_base = FT(0)
+
+    # ice thickness and mixed layer temperature changes due to atmosphereic and ocean fluxes
+    if h_ice > 0 # ice-covered
+        F_base = p.F0_base * (T_ml - p.T_base)
+        ΔT_ml = -(F_base + ocean_qflux) * Δt / (p.h_ml * p.ρc_ml)
+        Δh_ice = (F_atm - F_base) * Δt / p.L_ice
+    else # ice-free
+        ΔT_ml = -(F_atm + ocean_qflux) * Δt / (p.h_ml * p.ρc_ml)
+        Δh_ice = 0
+    end
+
+    # adjust if transition to ice-covered
+    if (T_ml[1] + ΔT_ml[1] < p.T_freeze)
+        Δh_ice = Δh_ice - (T_ml + ΔT_ml - p.T_freeze) * (p.h_ml * p.ρc_ml) / p.L_ice
+        ΔT_ml = p.T_freeze - T_ml
+    end
+
+    # adjust if transition to ice-free
+    if ((h_ice[1] > 0) & (h_ice[1] + Δh_ice[1] <= 0))
+        ΔT_ml = ΔT_ml - (h_ice + Δh_ice) * p.L_ice / (p.h_ml * p.ρc_ml)
+        Δh_ice = -h_ice
+    end
+
+    # solve for T_sfc
+    if (h_ice[1] + Δh_ice[1] > 0) #  surface is 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.k_ice / (h_ice + Δh_ice) * (p.T_base - T_sfc)
+        ΔT_sfc = (-F_atm + F_conductive) / (p.k_ice / (h_ice + Δh_ice) + ∂F_atm∂T_sfc)
+        if (T_sfc[1] + ΔT_sfc[1] > p.T_freeze)
+            ΔT_sfc = p.T_freeze - T_sfc
+        end
+        # surface is ice-covered, so update T_sfc as ice surface temperature
+
+        T_sfc += FT(ΔT_sfc)
+    else # ice-free, so update T_sfc as mixed layer temperature
+        T_sfc = FT(T_ml + ΔT_ml)
+    end
+
+    T_ml += FT(ΔT_ml)
+    h_ice += FT(Δh_ice)
+
+    return (; T_ml = FT(T_ml), h_ice = FT(h_ice), F_base = FT(F_base), T_sfc = FT(T_sfc))
+end
+
+get_∂F_atm∂T_sfc(p, T_sfc, Ya) = @. FT(4) * (FT(1) - p.α) * p.σ * T_sfc^3 + Ya.∂F_aero∂T_sfc
+
+function ∑tendencies_ice_stub(du, u, p, _)
+    dY = du
+    Y = u
+    FT = eltype(dY)
+
+    if p.Ya.prescribed_sic_data !== nothing # Prognostic eqn for T_sfc
+
+        @show "pres SIC"
+        params = p.Ya.params
+        F_aero = p.Ya.F_aero
+        F_rad = p.Ya.F_rad
+        ice_mask = p.Ya.ice_mask
+
+        F_conductive = @. params.k_ice / (params.h) * (params.T_base - Y.T_sfc) .* FT(0)
+        rhs = @. (-F_aero - F_rad + F_conductive) / (params.h * params.ρ * params.c)
+
+        # ∂F_atm∂T_sfc = get_∂F_atm∂T_sfc(p, T_sfc, Ya)
+        # rhs = @. (-F_aero - F_rad + F_conductive) / (p.k_ice / params.h + ∂F_atm∂T_sfc)
+        parent(dY.T_sfc) .= apply_mask.(parent(ice_mask), >, parent(rhs), FT(0), FT(0))
+    else
+        solve_ice!((; u = u, p = p), p.Ya.Δt) # timestepping outside of DeffEq (but DeffEq still used here for saving vars in `integ.sol`)
+
+        @. dY.T_ml = FT(0)
+        @. dY.h_ice = FT(0)
+        @. dY.T_sfc = FT(0)
+    end
+end
+
+function semtner_ice_init(
+    ::Type{FT},
+    tspan,
+    ocean_params;
+    stepper = Euler(),
+    nelements = 6,
+    npolynomial = 4,
+    dt = 0.02,
+    saveat = 1.0e10,
+    space = nothing,
+    mask = nothing,
+    prescribed_sic_data = nothing,
+) where {FT}
+
+    params = SemtnerParameters(
+        FT(2),
+        FT(1500.0),
+        FT(800.0),
+        FT(280.0),
+        FT(1e-3),
+        FT(1e-5),
+        FT(1500.0 * 800.0), #rho c
+        FT(120),
+        FT(273.16),
+        FT(3e8),
+        FT(ocean_params.h), # h_ml
+        FT(273.16),
+        FT(2),
+        FT(0.38),
+        FT(5.67e-8),
+    ) # TODO: better interface, use CLIMAParameters
+
+    ice_mask =
+        prescribed_sic_data !== nothing ? get_ice_mask.(prescribed_sic_data .- FT(25)) : ClimaCore.Fields.zeros(space) # here 25% and lower is considered ice free # TODO: generalize to a smooth function of ice fraction
+
+    #ice_mask = ClimaCore.Fields.zeros(space)
+    Y = slab_ice_space_init(FT, space, params)
+    Ya = (;
+        params = params,
+        F_aero = ClimaCore.Fields.zeros(space),
+        ∂F_aero∂T_sfc = ClimaCore.Fields.zeros(space),
+        F_rad = ClimaCore.Fields.zeros(space),
+        mask = mask,
+        ice_mask = ice_mask,
+        Δt = dt,
+        prescribed_sic_data = prescribed_sic_data,
+        F_base = ClimaCore.Fields.zeros(space),
+    ) #auxiliary
+
+    problem = OrdinaryDiffEq.ODEProblem(∑tendencies_ice_stub, Y, tspan, (; Ya = Ya))
+    integrator = OrdinaryDiffEq.init(problem, stepper, dt = dt, saveat = saveat)
+
+
+    SemtnerIceSimulation(params, Y, space, integrator)
+end
+
+get_ice_mask(h_ice) = h_ice > FT(0) ? FT(1) : FT(0)
+
+get_ml_energy(slab_sim, T_sfc) = T_sfc .* slab_sim.params.h_ml .* slab_sim.params.ρc_ml #slab_sim.params.ρ .* slab_sim.params.c .* T_sfc .* slab_sim.params.h
+
+get_ice_energy(slab_sim, T_sfc) = T_sfc ./ slab_sim.params.h .* slab_sim.params.k_ice
+
+get_dyn_ice_energy(seaice_sim, h_ice) = .-h_ice .* seaice_sim.integrator.p.Ya.params.L_ice
+
+# file-specific
+"""
+    clean_sic(SIC, _info)
+Ensures that the space of the SIC struct matches that of the mask, and converts the units from area % to area fraction.
+"""
+clean_sic(SIC, _info) = swap_space!(zeros(axes(_info.land_fraction)), SIC) ./ float_type_bcf(_info)(100.0)
+
+# setting that SIC < 0.5 is counted as ocean if binary remapping.
+get_ice_fraction(h_ice::FT, mono::Bool, threshold = 0.5) where {FT} =
+    mono ? h_ice : Regridder.binary_mask(h_ice, threshold = FT(threshold))
+
+# extensions required by Interfacer
+get_field(sim::SemtnerIceSimulation, ::Val{:surface_temperature}) = sim.integrator.u.T_sfc
+get_field(sim::SemtnerIceSimulation, ::Val{:surface_humidity}) = sim.integrator.p.q_sfc
+get_field(sim::SemtnerIceSimulation, ::Val{:roughness_momentum}) = sim.integrator.p.params.z0m
+get_field(sim::SemtnerIceSimulation, ::Val{:roughness_buoyancy}) = sim.integrator.p.params.z0b
+get_field(sim::SemtnerIceSimulation, ::Val{:beta}) = convert(eltype(sim.integrator.u), 1.0)
+get_field(sim::SemtnerIceSimulation, ::Val{:albedo}) = sim.integrator.p.params.α
+get_field(sim::SemtnerIceSimulation, ::Val{:area_fraction}) = sim.integrator.p.area_fraction
+get_field(sim::SemtnerIceSimulation, ::Val{:air_density}) = sim.integrator.p.ρ_sfc
+
+function update_field!(sim::SemtnerIceSimulation, ::Val{:area_fraction}, field::Fields.Field)
+    sim.integrator.p.area_fraction .= field
+end
+
+function update_field!(sim::SemtnerIceSimulation, ::Val{:turbulent_energy_flux}, field)
+    parent(sim.integrator.p.F_turb_energy) .= parent(field)
+end
+function update_field!(sim::SemtnerIceSimulation, ::Val{:radiative_energy_flux}, field)
+    parent(sim.integrator.p.F_radiative) .= parent(field)
+end
+function update_field!(sim::SemtnerIceSimulation, ::Val{:air_density}, field)
+    parent(sim.integrator.p.ρ_sfc) .= parent(field)
+end
+
+# extensions required by FieldExchanger
+step!(sim::SemtnerIceSimulation, t) = step!(sim.integrator, t - sim.integrator.t, true)
+reinit!(sim::SemtnerIceSimulation) = reinit!(sim.integrator)
+
+# extensions required by FluxCalculator (partitioned fluxes)
+function update_turbulent_fluxes_point!(sim::SemtnerIceSimulation, fields::NamedTuple, colidx::Fields.ColumnIndex)
+    (; F_turb_energy) = fields
+    @. sim.integrator.p.F_turb_energy[colidx] = F_turb_energy
+end
+
+"""
+    get_model_state_vector(sim::SemtnerIceSimulation)
+
+Extension of Checkpointer.get_model_state_vector to get the model state.
+"""
+function get_model_state_vector(sim::SemtnerIceSimulation)
+    return sim.integrator.u
+end