diff --git a/src/Physics/SEB/SEB.jl b/src/Physics/SEB/SEB.jl
index be02727a8ad2e93ec15f5b518b465b28f183c0c7..5327fb28aefbbbd17473e7a2bcb39d0c4e4c7ea5 100644
--- a/src/Physics/SEB/SEB.jl
+++ b/src/Physics/SEB/SEB.jl
@@ -12,6 +12,8 @@ using Parameters
import CryoGrid: BoundaryProcess, BoundaryStyle, Neumann, Top
import CryoGrid: initialcondition!, variables, boundaryvalue
+export Businger, HøgstrømSHEBA, Iterative, Analytical, Numerical
+
abstract type SolutionScheme end
"""
Byun 1990
@@ -42,8 +44,6 @@ SHEBA, Uttal et al., 2002, Grachev et al. 2007 (stable conditions)
"""
struct HøgstrømSHEBA <: StabilityFunctions end
-#abstract type SEBParams{TSolution,TStabFun} <: Params end
-
@with_kw struct SEBParams{TSolution,TStabFun} <: Params
# surface properties --> should be associated with the Stratigraphy and maybe made state variables
α::Float"1" = 0.2xu"1" # surface albedo [-]
@@ -71,30 +71,30 @@ struct HøgstrømSHEBA <: StabilityFunctions end
# type-dependet parameters
solscheme::TSolution = Iterative()
stabfun::TStabFun = HøgstrømSHEBA()
-
-# function SEBParams( ss::SolutionScheme, sf::StabilityFunctions )
-# solscheme = ss;
-# stabfun = sf;
-# new{typeof(solscheme),typeof(stabfun)}(solscheme, stabfun)
-# end
end
-
+"""
+Surface energy balance upper boundary condition.
+"""
struct SurfaceEnergyBalance{TSolution,TStabFun,F} <: BoundaryProcess
- forcing::F
+ forcings::F
sebparams::SEBParams{TSolution,TStabFun}
- function SurfaceEnergyBalance(Tair::TimeSeriesForcing, p::TimeSeriesForcing, q::TimeSeriesForcing,
- wind::TimeSeriesForcing, Lin::TimeSeriesForcing,
- Sin::TimeSeriesForcing, z::Float"m";
- kwargs...)
- forcing = (Tair = Tair, p = p, q = q, wind = wind, Lin = Lin, Sin = Sin, z = z);
+ function SurfaceEnergyBalance(
+ Tair::TimeSeriesForcing{Float"°C"},
+ p::TimeSeriesForcing,
+ q::TimeSeriesForcing,
+ wind::TimeSeriesForcing,
+ Lin::TimeSeriesForcing,
+ Sin::TimeSeriesForcing,
+ z::Float"m";
+ kwargs...
+ )
+ forcings = (Tair = Tair, p = p, q = q, wind = wind, Lin = Lin, Sin = Sin, z = z);
sebparams = SEBParams(;kwargs...);
- new{typeof(sebparams.solscheme),typeof(sebparams.stabfun),typeof(forcing)}(forcing, sebparams)
+ new{typeof(sebparams.solscheme),typeof(sebparams.stabfun),typeof(forcings)}(forcings, sebparams)
end
end
include("seb_heat.jl")
-export SEBParams, Businger, HøgstrømSHEBA, Iterative, Analytical, Numerical
-
end
diff --git a/src/Physics/SEB/seb_heat.jl b/src/Physics/SEB/seb_heat.jl
index 9264a78cbb4e51049d28f36635ffd29bc1aeaff0..c365443ba38ab9df1d8006c5c4652c48d0973e6d 100644
--- a/src/Physics/SEB/seb_heat.jl
+++ b/src/Physics/SEB/seb_heat.jl
@@ -25,23 +25,23 @@ BoundaryStyle(::Type{<:SurfaceEnergyBalance}) = Neumann()
"""
Top interaction, ground heat flux from surface energy balance. (no snow, no water body, no infiltration)
"""
-function boundaryvalue(seb::SurfaceEnergyBalance, ::Soil, ::Heat, stop, ssoil)
+function boundaryvalue(seb::SurfaceEnergyBalance, ::Top, ::Heat, ::Soil, stop, ssoil)
# TODO (optimize): pre-compute all forcings at time t here, then pass to functions
# 1. calculate radiation budget
# outgoing shortwave radiation as reflected
- @setscalar stop.Sout = let α = seb.sebparams.α, Sin = seb.forcing.Sin(stop.t);
+ @setscalar stop.Sout = let α = seb.sebparams.α, Sin = seb.forcings.Sin(stop.t);
-α * Sin # Eq. (2) in Westermann et al. (2016)
end
# outgoing longwave radiation composed of emitted and reflected radiation
- @setscalar stop.Lout = let ϵ = seb.sebparams.ϵ, σ = seb.sebparams.σ, T₀ = ssoil.T[1], Lin = seb.forcing.Lin(stop.t);
+ @setscalar stop.Lout = let ϵ = seb.sebparams.ϵ, σ = seb.sebparams.σ, T₀ = ssoil.T[1], Lin = seb.forcings.Lin(stop.t);
-ϵ * σ * T₀^4 - (1 - ϵ) * Lin # Eq. (3) in Westermann et al. (2016)
end
# net radiation budget
- @setscalar stop.Qnet = let Sin = seb.forcing.Sin(stop.t), Lin = seb.forcing.Lin(stop.t), Sout = getscalar(stop.Sout), Lout = getscalar(stop.Lout);
+ @setscalar stop.Qnet = let Sin = seb.forcings.Sin(stop.t), Lin = seb.forcings.Lin(stop.t), Sout = getscalar(stop.Sout), Lout = getscalar(stop.Lout);
Sin + Sout + Lin + Lout
end
@@ -93,8 +93,8 @@ Friction velocity according to Monin-Obukhov theory
"""
function ustar(seb::SurfaceEnergyBalance, stop)
let κ = seb.sebparams.κ,
- uz = seb.forcing.wind(stop.t), # wind speed at height z
- z = seb.forcing.z, # height z of wind forcing
+ uz = seb.forcings.wind(stop.t), # wind speed at height z
+ z = seb.forcings.z, # height z of wind forcing
zâ‚€ = seb.sebparams.zâ‚€, # aerodynamic roughness length [m]
Lstar = stop.Lstar |> getscalar;
κ * uz ./ (log(z / z₀) - Ψ_M(seb, z / Lstar, z₀ / Lstar)) # Eq. (7) in Westermann et al. (2016)
@@ -110,8 +110,8 @@ function Lstar(seb::SurfaceEnergyBalance{Iterative}, stop, ssoil)
g = seb.sebparams.g,
Râ‚ = seb.sebparams.Râ‚,
câ‚š = seb.sebparams.câ‚ / seb.sebparams.Ïâ‚, # specific heat capacity of air at constant pressure
- Tâ‚• = seb.forcing.Tair(stop.t), # air temperature at height z over surface
- p = seb.forcing.p(stop.t), # atmospheric pressure at surface
+ Tâ‚• = seb.forcings.Tair(stop.t), # air temperature at height z over surface
+ p = seb.forcings.p(stop.t), # atmospheric pressure at surface
ustar = stop.ustar |> getscalar,
Qâ‚‘ = stop.Qe |> getscalar,
Qâ‚• = stop.Qh |> getscalar,
@@ -133,12 +133,12 @@ function Lstar(seb::SurfaceEnergyBalance{Analytical}, stop, ssoil)
res = let g = seb.sebparams.g,
Râ‚ = seb.sebparams.Râ‚,
câ‚š = seb.sebparams.câ‚ / seb.sebparams.Ïâ‚, # specific heat capacity of air at constant pressure
- Tâ‚• = seb.forcing.Tair(stop.t), # air temperature at height z over surface
+ Tâ‚• = seb.forcings.Tair(stop.t), # air temperature at height z over surface
Tâ‚€ = ssoil.T[1], # surface temperature
- p = seb.forcing.p(stop.t), # atmospheric pressure at surface (height z)
- pâ‚€ = seb.forcing.p(stop.t), # normal pressure (for now assumed to be equal to p)
- uz = seb.forcing.wind(stop.t), # wind speed at height z
- z = seb.forcing.z, # height z of wind forcing
+ p = seb.forcings.p(stop.t), # atmospheric pressure at surface (height z)
+ pâ‚€ = seb.forcings.p(stop.t), # normal pressure (for now assumed to be equal to p)
+ uz = seb.forcings.wind(stop.t), # wind speed at height z
+ z = seb.forcings.z, # height z of wind forcing
zâ‚€ = seb.sebparams.zâ‚€, # aerodynamic roughness length [m]
Prâ‚€ = seb.sebparams.Prâ‚€, # turbulent Prandtl number
γₕ = seb.sebparams.γₕ,
@@ -186,13 +186,13 @@ Sensible heat flux, defined as positive if it is a flux towards the surface
function Q_H(seb::SurfaceEnergyBalance, stop, ssoil)
let κ = seb.sebparams.κ,
Râ‚ = seb.sebparams.Râ‚,
- Tâ‚• = seb.forcing.Tair(stop.t), # air temperature
+ Tâ‚• = seb.forcings.Tair(stop.t), # air temperature
Tâ‚€ = ssoil.T[1], # surface temperature
câ‚š = seb.sebparams.câ‚ / seb.sebparams.Ïâ‚, # specific heat capacity of air at constant pressure
- z = seb.forcing.z, # height at which forcing data are provided
+ z = seb.forcings.z, # height at which forcing data are provided
Lstar = stop.Lstar |> getscalar,
ustar = stop.ustar |> getscalar,
- p = seb.forcing.p(stop.t),
+ p = seb.forcings.p(stop.t),
zâ‚€ = seb.sebparams.zâ‚€,
Ïâ‚ = density_air(seb, Tâ‚•, p); # density of air at surface air temperature and surface pressure [kg/m^3]
@@ -212,18 +212,18 @@ function Q_E(seb::SurfaceEnergyBalance, stop, ssoil)
let κ = seb.sebparams.κ,
γ = seb.sebparams.γ,
Râ‚ = seb.sebparams.Râ‚,
- Tâ‚• = seb.forcing.Tair(stop.t), # air temperature at height z over surface
+ Tâ‚• = seb.forcings.Tair(stop.t), # air temperature at height z over surface
Tâ‚€ = ssoil.T[1], # surface temperature
- p = seb.forcing.p(stop.t), # atmospheric pressure at surface
- qâ‚• = seb.forcing.q(stop.t), # specific humidity at height h over surface
- z = seb.forcing.z, # height at which forcing data are provided
+ p = seb.forcings.p(stop.t), # atmospheric pressure at surface
+ qâ‚• = seb.forcings.q(stop.t), # specific humidity at height h over surface
+ z = seb.forcings.z, # height at which forcing data are provided
râ‚› = seb.sebparams.râ‚›, # surface resistance against evapotranspiration / sublimation [1/m]
Lstar = stop.Lstar |> getscalar,
ustar = stop.ustar |> getscalar,
Llg = L_lg(ssoil.T[1]),
Lsg = L_sg(ssoil.T[1]),
zâ‚€ = seb.sebparams.zâ‚€, # aerodynamic roughness length [m]
- Ïâ‚ = density_air(seb, seb.forcing.Tair(stop.t), seb.forcing.p(stop.t)); # density of air at surface air temperature and surface pressure [kg/m^3]
+ Ïâ‚ = density_air(seb, seb.forcings.Tair(stop.t), seb.forcings.p(stop.t)); # density of air at surface air temperature and surface pressure [kg/m^3]
q₀ = γ * estar(T₀) / p # saturation pressure of water/ice at the surface; Eq. (B1) in Westermann et al (2016)
râ‚ᵂ = (κ * ustar)^-1 * (log(z / zâ‚€) - Ψ_HW(seb, z / Lstar, zâ‚€ / Lstar)) # aerodynamic resistance Eq. (6) in Westermann et al. (2016)