Fix cglite snow depth not going to zero

Project.toml

 name = "CryoGrid"
 uuid = "a535b82e-5f3d-4d97-8b0b-d6483f5bebd5"
 authors = ["Brian Groenke <brian.groenke@awi.de>", "Jan Nitzbon <jan.nitzbon@awi.de>", "Moritz Langer <moritz.langer@awi.de>"]
-version = "0.22.1"
+version = "0.22.2"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"

examples/heat_freeW_lite_implicit.jl

@@ -10,6 +10,8 @@ using CryoGrid.LiteImplicit
 # Load forcings and build stratigraphy like before.
 forcings = loadforcings(CryoGrid.Presets.Forcings.Samoylov_ERA_MkL3_CCSM4_long_term);
 forcings = Base.rename(forcings, :Ptot => :precip)
+precip_values = forcings.precip.interpolant.coefs
+precip_values .*= 2
 z_top = -2.0u"m"
 z_bot = 1000.0u"m"
 upperbc = WaterHeatBC(
@@ -52,6 +54,20 @@ cg = Plots.cgrad(:copper,rev=true);
 Plots.plot(out.T[Z(Near(zs))], color=cg[LinRange(0.0,1.0,length(zs))]', ylabel="Temperature", title="", leg=false, dpi=150)
 Plots.plot(out.snowpack.dsn)
 
+integrator = init(prob, LiteImplicitEuler())
+for i in integrator
+    state = getstate(integrator)
+    last_above_ground_idx = findlast(<(0), cells(state.grid))
+    if state.snowpack.dsn[1] > 0.0 && state.snowpack.swe[last_above_ground_idx] <= 0 && state.snowpack.snowfall[1] <= 0
+        println("stopping")
+        break
+    end
+end
+
+CryoGrid.debug(true)
+
+step!(integrator)
+
 # CryoGridLite can also be embedded into integrators from OrdinaryDiffEq.jl via the `NLCGLite` nonlinear solver interface.
 # Note that these sovers generally will not be faster (in execution time) but may be more stable in some cases. Adaptive timestepping can be employed by
 # removing the `adaptive=false` argument.

src/Physics/Snow/snow_lite.jl

@@ -28,8 +28,8 @@ accumulation!(::Top, ::SnowBC, ::LiteSnowpack, ::SnowMassBalance, stop, ssnow) =
 
 # Declare snow mass balance variables for Lite scheme
 CryoGrid.variables(::LiteSnowpack, ::SnowMassBalance) = (
-    Prognostic(:swe, OnGrid(Cells), u"m"),
     Prognostic(:dsn, Scalar, u"m"),
+    Diagnostic(:swe, OnGrid(Cells), u"m"),
     Diagnostic(:ρsn, OnGrid(Cells), u"kg/m^3"),
     Diagnostic(:snowfall, Scalar, u"m/s"),
     Diagnostic(:T_ub, Scalar, u"°C"),
@@ -126,11 +126,16 @@ function CryoGrid.diagnosticstep!(snowpack::LiteSnowpack, state)
     new_dsn, new_ubc_idx = lite_snow!(state.θw, state.θwi, grid, dsn, para.dsn_max, sf, state.ρsn, para.ρsn_max, para.ρsn_k1, para.ρsn_k2, ρw, date)
     state.dsn .= new_dsn
     state.ubc_idx .= new_ubc_idx
+    total_swe = 0.0
     # calculate swe from updated water/ice contents
     @inbounds for i in eachindex(state.θwi)
-        state.swe[i] = (state.θwi[i] - state.θw[i])*min(dz[i], abs(dsn - grid[i+1]))
+        state.swe[i] = state.θwi[i]*dz[i]
+        total_swe += state.swe[i]
     end
     @. state.ρsn = ρw*state.swe / dz
+    if state.T_ub[1] > 1.0 && total_swe <= 0.0 && state.dsn[1] > 0.0
+        state.dsn .= 0.0
+    end
     return false
 end
 
@@ -138,7 +143,7 @@ function CryoGrid.initialcondition!(snowpack::LiteSnowpack, ::SnowMassBalance, s
     @. state.ρsn = snowpack.para.ρsn_0
     @. state.swe = 0.0
     @. state.dsn = 0.0
-    @. state.ubc_idx = length(state.T)
+    state.ubc_idx .= length(state.T)
 end
 
 function Hydrology.watercontent!(snow::LiteSnowpack, ::WaterBalance, state)