Fix more miscellaneous issues

examples/heat_freeW_lake_lite_implicit.jl

@@ -6,7 +6,12 @@ Plots.plotly()
 
 CryoGrid.debug(true)
 
-forcings = loadforcings(CryoGrid.Forcings.Samoylov_ERA_MkL3_CCSM4_long_term);
+raw_forcings = loadforcings(CryoGrid.Forcings.Samoylov_ERA_MkL3_CCSM4_long_term);
+Tair = raw_forcings.data.Tair
+Ptot = uconvert.(u"m/s", raw_forcings.data.Ptot)
+rainfall = Ptot.*(Tair .> 0u"°C")
+snowfall = Ptot.*(Tair .<= 0u"°C")
+forcings = rebuild(raw_forcings; Tair, rainfall, snowfall);
 soilprofile = SoilProfile(
     0.0u"m" => SimpleSoil(por=0.80,sat=0.9,org=0.75),
     0.1u"m" => SimpleSoil(por=0.80,sat=1.0,org=0.25),

examples/heat_freeW_lite_implicit.jl

@@ -7,11 +7,17 @@
 using CryoGrid
 using CryoGrid.LiteImplicit
 
-# Load forcings and build stratigraphy like before.
-forcings = loadforcings(CryoGrid.Forcings.Samoylov_ERA_MkL3_CCSM4_long_term);
-forcings = Base.rename(forcings, :Ptot => :precip)
-precip_values = forcings.precip.interpolant.coefs
-precip_values .*= 2
+# First we prepare the forcing data. This data includes only air temperature and
+# *total* precipitation, so we need to separate this into rainfall and snowfall
+# using air temperature.
+raw_forcings = loadforcings(CryoGrid.Forcings.Samoylov_ERA_MkL3_CCSM4_long_term);
+Tair = raw_forcings.data.Tair
+Ptot = uconvert.(u"m/s", raw_forcings.data.Ptot)
+rainfall = Ptot.*(Tair .> 0u"°C")
+snowfall = Ptot.*(Tair .<= 0u"°C")
+forcings = rebuild(raw_forcings; Tair, rainfall, snowfall);
+
+# Next we define the boundary conditions and stratigraphy.
 z_top = -2.0u"m"
 z_bot = 1000.0u"m"
 upperbc = WaterHeatBC(
@@ -29,7 +35,7 @@ soilprofile = SoilProfile(
 )
 heat = HeatBalance(heatop)
 water = WaterBalance()
-soil_layers = map(para -> Ground(para; heat, water), soilprofile)
+soil_layers = map(para -> Ground(para; heat, water), soilprofile);
 strat = @Stratigraphy(
     z_top => Top(upperbc),
     z_top => Snowpack(Snow.LiteGridded(); heat),
@@ -40,11 +46,11 @@ modelgrid = Grid(vcat(z_top:0.02u"m":-0.02u"m", CryoGrid.DefaultGrid_2cm))
 tile = Tile(strat, modelgrid, forcings, ssinit);
 
 # Since the solver can take daily timesteps, we can easily specify longer simulation time spans at minimal cost.
-# Here we specify a time span of 30 years.
-tspan = (DateTime(1990,10,1), DateTime(2020,10,1))
+# Here we specify a time span of 20 years.
+tspan = (DateTime(2000,10,1), DateTime(2020,10,1))
 u0, du0 = initialcondition!(tile, tspan);
-prob = CryoGridProblem(tile, u0, tspan, saveat=24*3600.0, savevars=(:T,:θw,:θwi,:kc,:ρsn))
-sol = @time solve(prob, LiteImplicitEuler())
+prob = CryoGridProblem(tile, u0, tspan, saveat=24*3600, savevars=(:T,:θw))
+sol = @time solve(prob, LiteImplicitEuler(), dt=24*3600)
 out = CryoGridOutput(sol)
 
 # Plot the results!
@@ -54,20 +60,6 @@ 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.

examples/heat_sfcc_richardseq_samoylov.jl

@@ -5,6 +5,8 @@
 # pressure gradients.
 
 using CryoGrid
+
+
 forcings = loadforcings(CryoGrid.Forcings.Samoylov_ERA_obs_fitted_1979_2014_spinup_extended_2044);
 tspan = (DateTime(2011,1,1),DateTime(2011,12,31))
 T0 = forcings.Tair(tspan[1])

src/IO/InputOutput.jl

@@ -67,7 +67,7 @@ export Input, InputProvider, InputFunctionProvider
 export inputs
 include("input.jl")
 
-export ForcingFormat, ForcingFormatJSON, ForcingFormatNCD
+export ForcingFormat, ForcingFormatJSON, ForcingFormatNCD, Interpolated1D
 export loadforcings
 include("forcings/forcings.jl")
 

src/IO/forcings/providers.jl

@@ -41,6 +41,11 @@ function Base.show(io::IO, mimetype::MIME"text/plain", data::Interpolated1D)
     show(io, mimetype, getfield(data, :data))
 end
 
+function DimensionalData.rebuild(data::Interpolated1D, default_interp_mode=Linear(); stripunits=true, interp_modes=(;), datasets...)
+    newdata = DimStack((; datasets...))
+    return Interpolated1D(newdata, default_interp_mode; stripunits, interp_modes...)
+end
+
 function inputmap(data::Interpolated1D)
     return map(_wrap_interpolant, getfield(data, :interpolants))
 end

src/Physics/Heat/heat_water.jl

@@ -18,7 +18,6 @@ advectiveflux(jw, T₁, T₂, cw, L) = jw*(cw*T₁*(jw > zero(jw)) + cw*T₂*(jw
 Adds advective energy fluxes for all internal grid cell faces.
 """
 function water_energy_advection!(jH, jw, T, cw::Real, L::Real)
-    water, heat = ps
     @inbounds for i in 2:length(jw)-1
         let jw = jw[i],
             T₁ = T[i-1],
@@ -86,18 +85,23 @@ function CryoGrid.interact!(
     end
 end
 
-# Flux calculation
 function CryoGrid.computeprognostic!(sub::SubSurface, ps::Coupled(WaterBalance, HeatBalance), state)
     water, heat = ps
-    cw = heatcapacitywater(sub, state)
-    L = heat.prop.L
     CryoGrid.computeprognostic!(sub, water, state)
     if heat.advection
+        cw = heatcapacitywater(sub, state)
+        L = heat.prop.L
         water_energy_advection!(state.jH, state.jw, state.T, cw, L)
     end
     CryoGrid.computeprognostic!(sub, heat, state)
 end
 
+function CryoGrid.computeprognostic!(sub::SubSurface, ps::Coupled(WaterBalance{NoFlow}, HeatBalance), state)
+    water, heat = ps
+    CryoGrid.computeprognostic!(sub, water, state)
+    CryoGrid.computeprognostic!(sub, heat, state)
+end
+
 function heatcapacitywater(sub::SubSurface, state)
     if hasproperty(state, :ch_w)
         # TODO: this is a temporary solution that works so long as the

src/Tiles/stratigraphy.jl

@@ -93,7 +93,9 @@ Base.NamedTuple(strat::Stratigraphy) = layers(strat)
 function Base.show(io::IO, ::MIME"text/plain", strat::Stratigraphy)
     print(io, "Stratigraphy:\n")
     for (k,v,b) in zip(keys(strat), values(strat), boundaries(strat))
-        print(repeat(" ", Base.indent_width), "$b: $k :: $(nameof(typeof(v)))\n")
+        procs = isa(processes(v), Tuple) ? processes(v) : (processes(v),)
+        procnames = map(p -> string(nameof(typeof(p))), procs)
+        print(repeat(" ", Base.indent_width), "$(b*u"m"): $k :: $(nameof(typeof(v))) with processes $(join(procnames, ", "))\n")
     end
 end
 

src/Tiles/tile.jl

@@ -34,7 +34,13 @@ Flatten.flattenable(::Type{<:Tile}, ::Type{Val{:inits}}) = true
 Flatten.flattenable(::Type{<:Tile}, ::Type{Val{:events}}) = true
 # all other fields are non-flattenable by default
 Flatten.flattenable(::Type{<:Tile}, ::Type{Val{name}}) where name = false
-Base.show(io::IO, ::MIME"text/plain", tile::Tile{TStrat,TGrid,TStates,TInits,TEvents,iip}) where {TStrat,TGrid,TStates,TInits,TEvents,iip} = print(io, "Tile (iip=$iip) with layers $(keys(tile.strat)), $TGrid, $TStrat")
+
+function Base.show(io::IO, mime::MIME"text/plain", tile::Tile{TStrat,TGrid,TStates,TInits,TEvents,TInputs,iip}) where {TStrat,TGrid,TStates,TInits,TEvents,TInputs,iip}
+    println(io, "Tile (iip=$iip) with layers $(keys(tile.strat))")
+    print("  "); show(io, mime, tile.grid); println();
+    print("  "); show(io, mime, tile.strat);
+    println("  metadata: $(tile.metadata)")
+end
 
 """
     Tile(