Partially fixed SEB. Byun scheme not yet working.

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)
+

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

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