Add new surface energy balance parameterizations

New implementation of the SEB module which allows to choose between an analytical solution of the Monin-Obukhov theory (Byun et al. 1990) and an iterative one which was used in CryoGrid3 (Westermann et al. 2013). In addition, one can choose between different stability functions (Businger 1971, Högström 1988, SHEBA) to compute the turbulent heat fluxes.

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