Skip to content
Snippets Groups Projects
Commit 87605a50 authored by Brian Groenke's avatar Brian Groenke
Browse files

Minor clean up

parent ecf81d1b
No related branches found
No related tags found
1 merge request!34Add new surface energy balance parameterizations
This commit is part of merge request !34. Comments created here will be created in the context of that merge request.
Hide whitespace changes
Inline Side-by-side
src/Physics/SEB/SEB.jl
+ 18
− 18
View file @ 87605a50
......@@ -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
src/Physics/SEB/seb_heat.jl
+ 21
− 21
View file @ 87605a50
......@@ -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)
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment