From ed5c6e997a8c4bad6bf7a8c3b89f0ebe9a0da871 Mon Sep 17 00:00:00 2001
From: Jan Nitzbon <jan.nitzbon@posteo.de>
Date: Tue, 2 Aug 2022 16:09:03 +0200
Subject: [PATCH] Partially fixed SEB. Byun scheme not yet working.
---
examples/heat_vgfc_seb_samoylov_custom.jl | 23 ++++++++---------------
src/Physics/SEB/SEB.jl | 1 +
src/Physics/SEB/seb_heat.jl | 18 ++++++++----------
3 files changed, 17 insertions(+), 25 deletions(-)
diff --git a/examples/heat_vgfc_seb_samoylov_custom.jl b/examples/heat_vgfc_seb_samoylov_custom.jl
index f4e924f3..4a3605f6 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,7 +28,7 @@ tspan = (DateTime(2010,1,1), DateTime(2011,1,1))
soilprofile, tempprofile = CryoGrid.Presets.SamoylovDefault
initT = initializer(:T, tempprofile)
strat = Stratigraphy(
- -z*u"m" => top(SurfaceEnergyBalance(Tair,pr,q,wind,Lin,Sin,z)),
+ -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")),
);
@@ -47,11 +39,12 @@ tspan = (DateTime(2010,10,30),DateTime(2011,10,30))
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 f6ea2e06..00bedece 100644
--- a/src/Physics/SEB/SEB.jl
+++ b/src/Physics/SEB/SEB.jl
@@ -13,6 +13,7 @@ import CryoGrid: initialcondition!, variables, boundaryvalue
using Unitful
+export SurfaceEnergyBalance, SEBParams
export Businger, HøgstrømSHEBA, Iterative, Analytical, Numerical
abstract type SolutionScheme end
diff --git a/src/Physics/SEB/seb_heat.jl b/src/Physics/SEB/seb_heat.jl
index 2c6e3390..71bfa08d 100644
--- a/src/Physics/SEB/seb_heat.jl
+++ b/src/Physics/SEB/seb_heat.jl
@@ -37,7 +37,7 @@ function boundaryvalue(seb::SurfaceEnergyBalance, ::Top, ::Heat, ::Soil, stop, s
# outgoing longwave radiation composed of emitted and reflected radiation
@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)
+ -ϵ * σ * normalize_temperature(T₀)^4 - (1 - ϵ) * Lin # Eq. (3) in Westermann et al. (2016)
end
# net radiation budget
@@ -68,11 +68,10 @@ 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)) : 611.2 * exp(22.46 * T / (272.62 + T)); # Eq. (B3) in Westermann et al. (2016)
@@ -110,13 +109,14 @@ 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.forcings.Tair(stop.t), # air temperature at height z over 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]
+ Ïâ‚ = 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
@@ -133,8 +133,8 @@ 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.forcings.Tair(stop.t), # air temperature at height z over surface
- Tâ‚€ = ssoil.T[1], # surface temperature
+ 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
@@ -223,7 +223,7 @@ function Q_E(seb::SurfaceEnergyBalance, stop, ssoil)
Llg = L_lg(ssoil.T[1]),
Lsg = L_sg(ssoil.T[1]),
zâ‚€ = seb.sebparams.zâ‚€, # aerodynamic roughness length [m]
- Ïâ‚ = 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]
+ Ïâ‚ = 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(seb, z / Lstar, zâ‚€ / Lstar)) # aerodynamic resistance Eq. (6) in Westermann et al. (2016)
@@ -323,5 +323,3 @@ function Ψ_M(seb::SurfaceEnergyBalance{T,Businger},ζâ‚::Float64, ζ₂::Float
end
end
end
-
-export SurfaceEnergyBalance
--
GitLab