diff --git a/examples/heat_vgfc_seb_samoylov_custom.jl b/examples/heat_vgfc_seb_samoylov_custom.jl
index 4364f781921b471dd84ecbb3b5a0cad4d0c09fa0..4a3605f6dd93c4f295a3977e305f60834036039b 100644
--- a/examples/heat_vgfc_seb_samoylov_custom.jl
+++ b/examples/heat_vgfc_seb_samoylov_custom.jl
@@ -12,17 +12,9 @@ soilprofile = SoilProfile(
0.4u"m" => soilparameters(xic=0.30,por=0.55,sat=1.0,org=0.25), #(θwi=0.80,θm=0.15,θo=0.05,ϕ=0.55),
3.0u"m" => soilparameters(xic=0.0,por=0.50,sat=1.0,org=0.0), #(θwi=0.50,θm=0.50,θo=0.0,ϕ=0.50),
10.0u"m" => soilparameters(xic=0.0,por=0.30,sat=1.0,org=0.0), #(θwi=0.30,θm=0.70,θo=0.0,ϕ=0.30),
-),
+);
# mid-winter temperature profile
-tempprofile = TemperatureProfile(
- 0.01u"m" => -15.9u"°C",
- 0.05u"m" => -15.49u"°C",
- 0.10u"m" => -15.32u"°C",
- 0.20u"m" => -14.44u"°C",
- 0.30u"m" => -14.18u"°C",
- 0.40u"m" => -13.50u"°C",
- 1000.0u"m" => 10.2u"°C",
-)
+tempprofile = CryoGrid.Presets.SamoylovDefault.tempprofile
forcings = loadforcings(CryoGrid.Presets.Forcings.Samoylov_ERA_obs_fitted_1979_2014_spinup_extended_2044, :Tair => u"°C");
Tair = TimeSeriesForcing(ustrip.(forcings.data.Tair), forcings.timestamps, :Tair);
# assume other forcings don't (yet) have units
@@ -36,22 +28,23 @@ tspan = (DateTime(2010,1,1), DateTime(2011,1,1))
soilprofile, tempprofile = CryoGrid.Presets.SamoylovDefault
initT = initializer(:T, tempprofile)
strat = Stratigraphy(
- -2.0u"m" => top(SurfaceEnergyBalance(Tair,pr,q,wind,Lin,Sin,z)),
- Tuple(z => subsurface(:soil1, Soil(para=para), Heat(:H, freezecurve=SFCC(DallAmico()))) for (z,para) in soilprofile),
+ -z*u"m" => top(SurfaceEnergyBalance(Tair,pr,q,wind,Lin,Sin,z,solscheme=Analytical(),stabfun=Businger())),
+ Tuple(knot.depth => subsurface(Symbol(:soil,i), Soil(para=knot.value), Heat(:H, freezecurve=SFCC(DallAmico()))) for (i,knot) in enumerate(soilprofile)),
1000.0u"m" => bottom(GeothermalHeatFlux(0.053u"J/s/m^2")),
);
grid = Grid(gridvals);
-model = Tile(strat, grid, initT);
+tile = Tile(strat, grid, initT);
# define time span
tspan = (DateTime(2010,10,30),DateTime(2011,10,30))
-p = parameters(model)
-u0, du0 = initialcondition!(model, tspan, p)
+p = parameters(tile)
+u0, du0 = initialcondition!(tile, tspan, p)
# CryoGrid front-end for ODEProblem
-prob = CryoGridProblem(model,u0,tspan,p,savevars=(:T,))
+prob = CryoGridProblem(tile,u0,tspan,p,step_limiter=nothing,savevars=(:T,))
# solve with forward Euler (w/ CFL) and construct CryoGridOutput from solution
out = @time solve(prob, Euler(), dt=2*60.0, saveat=24*3600.0, progress=true) |> CryoGridOutput;
# Plot it!
-zs = [1:10...,20:10:100...]
+zs = [1,5,10,15,20:10:100...]
cg = Plots.cgrad(:copper,rev=true);
-plot(out.H[Z(zs)], color=cg[LinRange(0.0,1.0,length(zs))]', ylabel="Enthalpy", leg=false, dpi=150)
-plot(out.T[Z(zs)], color=cg[LinRange(0.0,1.0,length(zs))]', ylabel="Temperature", leg=false, size=(800,500), dpi=150)
+plot(ustrip.(out.H[Z(zs)]), color=cg[LinRange(0.0,1.0,length(zs))]', ylabel="Enthalpy", leg=false, dpi=150)
+plot(ustrip.(out.T[Z(zs)]), color=cg[LinRange(0.0,1.0,length(zs))]', ylabel="Temperature", leg=false, size=(800,500), dpi=150)
+
diff --git a/src/Physics/SEB/SEB.jl b/src/Physics/SEB/SEB.jl
index 03fefd8ad596196f1a27e7810ce9ebf3216ff696..00bedeceb81186647de439010c6ec8a8488110ba 100644
--- a/src/Physics/SEB/SEB.jl
+++ b/src/Physics/SEB/SEB.jl
@@ -1,19 +1,52 @@
module SEB
-using ..HeatConduction: Heat
-using ..Physics
-using ..Boundaries
using CryoGrid.InputOutput: Forcing
+using CryoGrid.Physics
+using CryoGrid.Physics.Boundaries
+using CryoGrid.Physics.HeatConduction
using CryoGrid.Physics.Soils
using CryoGrid.Numerics
using CryoGrid.Utils
import CryoGrid: BoundaryProcess, BoundaryStyle, Neumann, Top
-import CryoGrid: initialcondition!, variables
+import CryoGrid: initialcondition!, variables, boundaryvalue
using Unitful
-Base.@kwdef struct SEBParams <: IterableStruct
+export SurfaceEnergyBalance, SEBParams
+export Businger, HøgstrømSHEBA, Iterative, Analytical, Numerical
+
+abstract type SolutionScheme end
+"""
+Byun 1990
+"""
+struct Analytical <: SolutionScheme end
+
+"""
+equation by Westermann2016, but use Newton solver
+- not implemented yet
+"""
+struct Numerical <: SolutionScheme end
+
+"""
+Westermann2016, use info from last time step
+"""
+struct Iterative <: SolutionScheme end
+
+abstract type StabilityFunctions end
+
+"""
+Businger 1971
+"""
+struct Businger <: StabilityFunctions end
+
+"""
+Høgstrøm 1988 (unstable conditions)
+SHEBA, Uttal et al., 2002, Grachev et al. 2007 (stable conditions)
+"""
+struct HøgstrømSHEBA <: StabilityFunctions end
+
+Base.@kwdef struct SEBParams{TSolution,TStabFun}
# surface properties --> should be associated with the Stratigraphy and maybe made state variables
α::Float"1" = 0.2xu"1" # surface albedo [-]
ϵ::Float"1" = 0.97xu"1" # surface emissivity [-]
@@ -30,28 +63,43 @@ Base.@kwdef struct SEBParams <: IterableStruct
# material properties (assumed to be constant)
Ïâ‚::Float"kg/m^3" = 1.293xu"kg/m^3" # density of air at standard pressure and 0°C [kg/m^3]
câ‚::Float"J/(m^3*K)" = 1005.7xu"J/(kg*K)" * Ïâ‚ # volumetric heat capacity of dry air at standard pressure and 0°C [J/(m^3*K)]
+
+ Prâ‚€::Float64 = 0.74 # turbulent Prandtl number
+ βₘ::Float64 = 4.7
+ βₕ::Float64 = βₘ/Pr₀
+ γₘ::Float64 = 15.0
+ γₕ::Float64 = 9.0
+
+ # type-dependent parameters
+ solscheme::TSolution = Iterative()
+ stabfun::TStabFun = HøgstrømSHEBA()
end
-struct SurfaceEnergyBalance{F} <: BoundaryProcess{Heat}
- forcing::F
- sebparams::SEBParams
+"""
+ SurfaceEnergyBalance{TSolution,TStabFun,F} <: BoundaryProcess{Heat}
+
+Surface energy balance upper boundary condition.
+"""
+struct SurfaceEnergyBalance{TSolution,TStabFun,F} <: BoundaryProcess{Heat}
+ forcings::F
+ sebparams::SEBParams{TSolution,TStabFun}
+ SurfaceEnergyBalance(forcings::NamedTuple, sebparams::SEBParams) = new{typeof(sebparams.solscheme),typeof(sebparams.stabfun),typeof(forcings)}(forcings, sebparams)
function SurfaceEnergyBalance(
- Tair::Forcing,
+ Tair::Forcing{Float"°C"},
p::Forcing,
q::Forcing,
wind::Forcing,
Lin::Forcing,
Sin::Forcing,
- z::Float"m",
- params::SEBParams=SEBParams()
+ z::Float"m";
+ kwargs...
)
- forcing = (Tair = Tair, p = p, q = q, wind = wind, Lin = Lin, Sin = Sin, z = z);
- new{typeof(forcing)}(forcing, params)
+ 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(forcings)}(forcings, sebparams)
end
end
-BoundaryStyle(::Type{<:SurfaceEnergyBalance}) = Neumann()
-
include("seb_heat.jl")
-end
\ No newline at end of file
+end
diff --git a/src/Physics/SEB/seb_heat.jl b/src/Physics/SEB/seb_heat.jl
index 96f113db1a840b660d32656fe5ce99bb09f32710..71bfa08d41e5bf61a628299722c3256da22adc78 100644
--- a/src/Physics/SEB/seb_heat.jl
+++ b/src/Physics/SEB/seb_heat.jl
@@ -20,33 +20,35 @@ function initialcondition!(top::Top, seb::SurfaceEnergyBalance, state)
@setscalar state.ustar = 10.;
end
+BoundaryStyle(::Type{<:SurfaceEnergyBalance}) = Neumann()
+
"""
Top interaction, ground heat flux from surface energy balance. (no snow, no water body, no infiltration)
"""
-function (seb::SurfaceEnergyBalance)(top::Top, soil::Soil, heat::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);
- -α * Sin # Eq. (2) in Westermann et al. (2016)
+ @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);
- -ϵ * σ * T₀^4 - (1 - ϵ) * Lin # Eq. (3) in Westermann et al. (2016)
+ @setscalar stop.Lout = let ϵ = seb.sebparams.ϵ, σ = seb.sebparams.σ, T₀ = ssoil.T[1], Lin = seb.forcings.Lin(stop.t);
+ -ϵ * σ * normalize_temperature(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
# 2. calcuate turbulent heat flux budget
# determine atmospheric stability conditions
- @setscalar stop.ustar = ustar(seb, stop);
@setscalar stop.Lstar = Lstar(seb, stop, ssoil);
+ @setscalar stop.ustar = ustar(seb, stop);
# sensible heat flux
@setscalar stop.Qh = Q_H(seb, stop, ssoil);
@@ -66,15 +68,12 @@ end
"""
Density of air at given tempeature and pressure
"""
-density_air(seb::SurfaceEnergyBalance,T::Real"°C",p::Real"Pa") = p / (T * seb.sebparams.Râ‚);
+density_air(seb::SurfaceEnergyBalance,T::Real"°C",p::Real"Pa") = p / (normalize_temperature(T) * seb.sebparams.Râ‚);
"""
Saturation pressure of water/ice according to the empirical August-Roche-Magnus formula
-Note: T is passed [K] and converted to [°C]
"""
-estar(T::Real"°C") = ( (T > 0) ? 611.2 * exp(17.62 * T / (243.12 + T))# Eq. (B3) in Westermann et al. (2016)
-: 611.2 * exp(22.46 * T / (272.62 + T)) ;
-)
+estar(T::Real"°C") = (T > 0) ? 611.2 * exp(17.62 * T / (243.12 + T)) : 611.2 * exp(22.46 * T / (272.62 + T)); # Eq. (B3) in Westermann et al. (2016)
"""
Latent heat of evaporation/condensation of water in [J/kg]
@@ -93,30 +92,88 @@ 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(z / Lstar, z₀ / Lstar)) # Eq. (7) in Westermann et al. (2016)
+ κ * uz ./ (log(z / z₀) - Ψ_M(seb, z / Lstar, z₀ / Lstar)) # Eq. (7) in Westermann et al. (2016)
end
end
"""
-Obukhov length according to Monin-Obukhov theory
+Obukhov length according to Monin-Obukhov theory, iterative determination as in CryoGrid3 / Westermann et al. 2016
+- uses the turubulent fluxes Qe and Qh as well as the friction velocity of the previous time step
"""
-function Lstar(seb::SurfaceEnergyBalance, stop, ssoil)
+function Lstar(seb::SurfaceEnergyBalance{Iterative}, stop, ssoil)
res = let κ = seb.sebparams.κ,
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
+ Tair = seb.forcings.Tair(stop.t),
+ Tâ‚• = normalize_temperature(Tair), # 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,
Llg = L_lg(ssoil.T[1]),
- Ïâ‚ = density_air(seb, Tâ‚•, p); # density of air at surface air temperature and surface pressure [kg/m^3]
- -Ïâ‚ * câ‚š * Tâ‚• / (κ * g) * ustar^3 / (Qâ‚• + 0.61 * câ‚š / Llg * Tâ‚• * Qâ‚‘) # Eq. (8) in Westermann et al. (2016)
+ Ïâ‚ = density_air(seb, Tair, p); # density of air at surface air temperature and surface pressure [kg/m^3]
+ -Ïâ‚ * câ‚š * Tâ‚• / (κ * g) * ustar^3 / (Qâ‚• + 0.61 * câ‚š / Llg * Tâ‚• * Qâ‚‘) # Eq. (8) in Westermann et al. (2016)
+ end
+ # upper and lower limits for Lstar
+ res = (abs(res) < 1e-7) ? sign(res) * 1e-7 : res
+ res = (abs(res) > 1e+7) ? sign(res) * 1e+7 : res
+end
+
+"""
+Obukhov length, analytical solution of Monin-Obukhov theory according to Byun 1990
+- implicitly assumes the Businger 1971 stability functions
+- only uses the current surface temperature as input
+"""
+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â‚• = normalize_temperature(seb.forcings.Tair(stop.t)), # air temperature at height z over surface
+ Tâ‚€ = normalize_temperature(ssoil.T[1]), # surface temperature
+ 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.γₕ,
+ γₘ = seb.sebparams.γₘ,
+ βₕ = seb.sebparams.βₕ,
+ βₘ = seb.sebparams.βₘ;
+
+ Θₕ = Tâ‚• * (pâ‚€/p)^(Râ‚/câ‚š); # potential temperature (for now identical to actual temperature)
+ Θ₀ = Tâ‚€ * (pâ‚€/p)^(Râ‚/câ‚š);
+
+ # calcuate bulk Richardson number
+ Ri_b = g / Θ₀ * (Θₕ - Θ₀) * (z - z₀) / uz^2; # bulk Richardson number, eq. (9) in Byun 1990
+
+ # calulate ζ
+ a = (z / (z-zâ‚€)) * log(z/zâ‚€);
+ if Ri_b>0 # stable conditions
+ b = 2 * βₕ * (βₘ * Ri_b - 1);
+ c = -(2 * βₕ * Ri_b - 1) - (1 + (4 * (βₕ - βₘ) * Ri_b) / Pr₀ )^0.5
+ ζ = a / b * c; # eq. (19) in Byun 1990
+ else # unstable conditions
+ s_b = Ri_b / Prâ‚€; # eq. (30) in Byun 1990
+ Q_b = 1/9 * ( 1/γₘ^2 + 3 * γₕ/γₘ * s_b^2 ); # eq. (31) in Byun 1990
+ P_b = 1/54 * ( -2/γₘ^3 + 9/γₘ * (3 - γₕ/γₘ) * s_b^2 ); # eq. (32) in Byun 1990
+ θ_b = acos( P_b / Q_b^(3/2) ); # eq. (33) in Byun 1990
+ T_b = ( (P_b^2 - Q_b^3)^0.5 + abs(P_b) )^(1/3); # eq. (34) in Byun 1990
+ if Q_b^3-P_b^2 >= 0
+ ζ = a * (-2 * Q_b^0.5 * cos(θ_b/3) + 1/(3*γₘ) ) # eq. (28) in Byun 1990
+ else
+ ζ = a * (-(T_b + Q_b/T_b ) + 1/(3*γₘ) ); # eq. (29) in Byun 1990
+ end
+ end
+
+ # calculate Lstar
+ z / ζ;
+
end
# upper and lower limits for Lstar
res = (abs(res) < 1e-7) ? sign(res) * 1e-7 : res
@@ -129,17 +186,17 @@ 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]
- râ‚á´´ = (κ * ustar)^-1 * (log(z / zâ‚€) - Ψ_HW(z / Lstar, zâ‚€ / Lstar)) # Eq. (6) in Westermann et al. (2016)
+ râ‚á´´ = (κ * ustar)^-1 * (log(z / zâ‚€) - Ψ_HW(seb, z / Lstar, zâ‚€ / Lstar)) # Eq. (6) in Westermann et al. (2016)
# calculate Q_H
-Ïâ‚ * câ‚š * (Tâ‚• - Tâ‚€) / râ‚á´´ # Eq. (4) in Westermann et al. (2016)
@@ -155,28 +212,28 @@ 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, 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(z / Lstar, zâ‚€ / Lstar)) # aerodynamic resistance Eq. (6) in Westermann et al. (2016)
- L = (Tâ‚€ <= 0.0) ? Lsg : Llg # latent heat of sublimation/resublimation or evaporation/condensation [J/kg]
+ râ‚ᵂ = (κ * ustar)^-1 * (log(z / zâ‚€) - Ψ_HW(seb, z / Lstar, zâ‚€ / Lstar)) # aerodynamic resistance Eq. (6) in Westermann et al. (2016)
+ L = (Tâ‚€ <= 0.0) ? Lsg : Llg # latent heat of sublimation/resublimation or evaporation/condensation [J/kg]
# calculate Q_E
res = (qâ‚• > qâ‚€) ? -Ïâ‚ * L * (qâ‚• - qâ‚€) / (râ‚ᵂ) : # Eq. (5) in Westermann et al. (2016) # condensation / deposition (no aerodynamics resistance)
-Ïâ‚ * L * (qâ‚• - qâ‚€) / (râ‚ᵂ + râ‚›) ; # evaporation / sublimation (account for surface resistance against evapotranspiration/sublimation)
- res = (Tâ‚€ <= 0.0) ? 0 : res ; # for now: set sublimation and deposition to zero.
+ res = (Tâ‚€ <= 0.0) ? 0 : res ; # for now: set sublimation and deposition to zero.
end
end
@@ -185,7 +242,7 @@ Integrated stability function for heat/water transport
Högström, 1988
SHEBA, Uttal et al., 2002, Grachev et al. 2007
"""
-function Ψ_HW(ζâ‚::Real, ζ₂::Real)
+function Ψ_HW(seb::SurfaceEnergyBalance{T,HøgstrømSHEBA}, ζâ‚::Float64, ζ₂::Float64) where T
if ζ₠<= 0 # neutral and unstable conditions (according to Høgstrøm, 1988)
# computed using WolframAlpha command "Integrate[ (1-(0.95*(1-11.6x)^(-1/2)))/x ] assuming x<0"
res = real( log(Complex(ζâ‚)) + 1.9 * atanh(Complex((1 - 11.6 * ζâ‚)^0.5)) -
@@ -205,8 +262,10 @@ end
"""
Integrated stability function for momentum transport
+ Högström, 1988
+ SHEBA, Uttal et al., 2002, Grachev et al. 2007
"""
-function Ψ_M(ζâ‚::Real, ζ₂::Real)
+function Ψ_M(seb::SurfaceEnergyBalance{T,HøgstrømSHEBA}, ζâ‚::Float64, ζ₂::Float64) where T
if ζ₠<= 0 # neutral and unstable conditions (according to Høgstrøm, 1988)
# computed using WolframAlpha command "Integrate[ (1-(1-19.3x)^(-1/4))/x ] assuming x<0"
# log(ζâ‚) - 2*atan((1-19.3*ζâ‚)^(1/4)) + 2*atanh((1-19.3*ζâ‚)^(1/4)) -
@@ -230,4 +289,37 @@ function Ψ_M(ζâ‚::Real, ζ₂::Real)
end
end
-export SurfaceEnergyBalance
+
+"""
+Integrated stability function for heat/water transport
+ Businger 1971
+"""
+function Ψ_HW(seb::SurfaceEnergyBalance{T,Businger},ζâ‚::Float64, ζ₂::Float64) where T
+ let γₕ = seb.sebparams.γₕ,
+ βₕ = seb.sebparams.βₕ;
+ if ζ₠<= 0 # neutral and unstable conditions (according to Businger, 1971)
+ return 2*log( ( (1-γₕ*ζâ‚)^0.5 + 1 ) / ( (1-γₕ*ζ₂)^0.5 + 1 ) ); # eq. (15) in Byun 1990
+ else # stable stratification (according to Busigner, 1971)
+ return -βₕ*(ζâ‚-ζ₂); # eq. (13) in Byun 1990
+ end
+ end
+end
+
+"""
+Integrated stability function for momentum transport
+ Busigner 1971
+"""
+function Ψ_M(seb::SurfaceEnergyBalance{T,Businger},ζâ‚::Float64, ζ₂::Float64) where T
+ let γₘ = seb.sebparams.γₘ,
+ βₘ = seb.sebparams.βₘ;
+ if ζ₠<= 0 # neutral and unstable conditions (according to Businger, 1971)
+ x=(1-γₘ*ζâ‚)^0.25;
+ x₀=(1-γₘ*ζ₂)^0.25;
+ return 2*log( (1+x)/(1+xâ‚€) ) + log( (1+x^2) / (1+xâ‚€^2) ) -
+ 2*atan(x) + 2*atan(x) ; # eq. (14) in Byun 1990
+
+ else # stable stratification (according to Busigner, 1971)
+ return -βₘ*(ζâ‚-ζ₂); # eq. (12) in Byun 1990
+ end
+ end
+end