From 6eb10824ba57f6580a7f000313799613c3e8c876 Mon Sep 17 00:00:00 2001
From: Brian Groenke <brian.groenke@awi.de>
Date: Fri, 13 May 2022 17:26:39 +0200
Subject: [PATCH] Fix post-merge issues

---
 examples/heat_vgfc_seb_samoylov_custom.jl | 20 +++++++--------
 src/Physics/SEB/SEB.jl                    | 31 +++++++++++++----------
 2 files changed, 27 insertions(+), 24 deletions(-)

diff --git a/examples/heat_vgfc_seb_samoylov_custom.jl b/examples/heat_vgfc_seb_samoylov_custom.jl
index a1bf4f7b..fba30648 100644
--- a/examples/heat_vgfc_seb_samoylov_custom.jl
+++ b/examples/heat_vgfc_seb_samoylov_custom.jl
@@ -12,7 +12,7 @@ soilprofile = SoilProfile(
     0.4u"m" => CharacteristicFractions(xic=0.30,por=0.55,sat=1.0,org=0.25), #(θw=0.80,θm=0.15,θo=0.05,ϕ=0.55),
     3.0u"m" => CharacteristicFractions(xic=0.0,por=0.50,sat=1.0,org=0.0), #(θw=0.50,θm=0.50,θo=0.0,ϕ=0.50),
     10.0u"m" => CharacteristicFractions(xic=0.0,por=0.30,sat=1.0,org=0.0), #(θw=0.30,θm=0.70,θo=0.0,ϕ=0.30),
-),
+)
 # mid-winter temperature profile
 tempprofile = TemperatureProfile(
     0.01u"m" => -15.9u"°C",
@@ -36,22 +36,20 @@ tspan = (DateTime(2010,1,1), DateTime(2011,1,1))
 soilprofile, tempprofile = CryoGrid.Presets.SamoylovDefault
 initT = initializer(:T, tempprofile)
 strat = Stratigraphy(
-    -z => top(SurfaceEnergyBalance(Tair,pr,q,wind,Lin,Sin,z,solscheme=Analytical(),stabfun=Businger())),
-    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)),
+    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,))
-# solve with forward Euler (w/ CFL) and construct CryoGridOutput from solution
-out = @time solve(prob, Euler(), dt=2*60.0, callback=CFLStepLimiter(model), saveat=24*3600.0, progress=true) |> CryoGridOutput;
+prob = CryoGridProblem(tile, u0, tspan, p, savevars=(:T,))
+# solve with forward Euler and construct CryoGridOutput from solution
+out = @time solve(prob, Euler(), dt=120.0, saveat=24*3600.0, progress=true) |> CryoGridOutput;
 # Plot it!
 zs = [1:10...,20:10:100...]
 cg = Plots.cgrad(:copper,rev=true);
diff --git a/src/Physics/SEB/SEB.jl b/src/Physics/SEB/SEB.jl
index 14d3ee1f..f6ea2e06 100644
--- a/src/Physics/SEB/SEB.jl
+++ b/src/Physics/SEB/SEB.jl
@@ -1,9 +1,9 @@
 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
@@ -11,6 +11,8 @@ using CryoGrid.Utils
 import CryoGrid: BoundaryProcess, BoundaryStyle, Neumann, Top
 import CryoGrid: initialcondition!, variables, boundaryvalue
 
+using Unitful
+
 export Businger, HøgstrømSHEBA, Iterative, Analytical, Numerical
 
 abstract type SolutionScheme end
@@ -43,7 +45,7 @@ SHEBA, Uttal et al., 2002, Grachev et al. 2007 (stable conditions)
 """
 struct HøgstrømSHEBA <: StabilityFunctions end
 
-Base.@kwdef struct SEBParams
+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 [-]
@@ -59,32 +61,35 @@ Base.@kwdef struct SEBParams
 
     # 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)]
+    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
+    Prâ‚€::Float64 = 0.74                             # turbulent Prandtl number
     βₘ::Float64 = 4.7
     βₕ::Float64 = βₘ/Pr₀
     γₘ::Float64 = 15.0
     γₕ::Float64 = 9.0
 
-    # type-dependet parameters
+    # type-dependent parameters
     solscheme::TSolution = Iterative()
     stabfun::TStabFun = HøgstrømSHEBA()
 end
 
 """
+    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::TimeSeriesForcing{Float"°C"},
-        p::TimeSeriesForcing,
-        q::TimeSeriesForcing,
-        wind::TimeSeriesForcing,
-        Lin::TimeSeriesForcing,
-        Sin::TimeSeriesForcing,
+        Tair::Forcing{Float"°C"},
+        p::Forcing,
+        q::Forcing,
+        wind::Forcing,
+        Lin::Forcing,
+        Sin::Forcing,
         z::Float"m";
         kwargs...
     )
-- 
GitLab