revert to old branch status

src/standalone/Snow/snow_parameterizations.jl

@@ -42,20 +42,19 @@ end
         p,
     ) where {FT}
 
-Returns the surface height of the `Snow` model; surface (ground) elevation
-and ignores snow depth (CLM).
+Returns the surface height of the `Snow` model; this returns the depth
+of the snow and hence implicitly treats the surface (ground) elevation as
+at zero.
 
 Once topography or land ice is incorporated, this will need to change to
-z_sfc + land_ice_depth. Note that land ice can
+z_sfc + land_ice_depth + snow_depth. Note that land ice can 
 be ~1-3 km thick on Greenland/
 
-In order to compute surface fluxes, this cannot be large than the
-height of the atmosphere measurement location (z_atmos > z_land always).
-
-This assumes that the surface elevation is zero.
+In order to compute surface fluxes, this cannot be large than the 
+height of the atmosphere measurement location (z_atmos > z_land always). 
 """
-function ClimaLand.surface_height(model::SnowModel{FT}, Y, p) where {FT}
-    return FT(0)
+function ClimaLand.surface_height(model::SnowModel, Y, p)
+    return snow_depth(model.parameters.density, Y, p, model.parameters)
 end
 
 
@@ -80,7 +79,7 @@ end
 """
     ClimaLand.surface_temperature(model::SnowModel, Y, p)
 
-a helper function which returns the surface temperature for the snow
+a helper function which returns the surface temperature for the snow 
 model, which is stored in the aux state.
 """
 function ClimaLand.surface_temperature(model::SnowModel, Y, p, t)
@@ -164,23 +163,15 @@ Computes the thermal conductivity, given the density
 of the snow, according to Equation
 5.33 from Bonan's textbook, which in turn is taken from
 Jordan (1991).
-
-We have adjusted the original equation to make the coefficients
-non-dimensional by multiplying by the first by x = ρ_ice/ρ_ice
-and the second by x², with ρ_ice in kg/m³.
-
-When ρ_snow = ρ_ice, we recover κ_snow = κ_ice.
 """
 function snow_thermal_conductivity(
     ρ_snow::FT,
     parameters::SnowParameters{FT},
 ) where {FT}
     _κ_air = FT(LP.K_therm(parameters.earth_param_set))
-    _ρ_ice = FT(LP.ρ_cloud_ice(parameters.earth_param_set))
     κ_ice = parameters.κ_ice
     return _κ_air +
-           (FT(0.07) * (ρ_snow / _ρ_ice) + FT(0.93) * (ρ_snow / _ρ_ice)^2) *
-           (κ_ice - _κ_air)
+           (FT(7.75e-5) * ρ_snow + FT(1.105e-6) * ρ_snow^2) * (κ_ice - _κ_air)
 end
 
 """
@@ -341,7 +332,7 @@ end
 """
     energy_from_q_l_and_swe(S::FT, q_l::FT, parameters) where {FT}
 
-A helper function for compute the snow energy per unit area, given snow
+A helper function for compute the snow energy per unit area, given snow 
 water equivalent S, liquid fraction q_l,  and snow model parameters.
 
 Note that liquid water can only exist at the freezing point in this model,
@@ -363,7 +354,7 @@ end
 """
     energy_from_T_and_swe(S::FT, T::FT, parameters) where {FT}
 
-A helper function for compute the snow energy per unit area, given snow
+A helper function for compute the snow energy per unit area, given snow 
 water equivalent S, bulk temperature T, and snow model parameters.
 
 If T = T_freeze, we return the energy as if q_l = 0.
@@ -516,10 +507,6 @@ function dzdt(density::Anderson1976{FT}, model::SnowModel{FT}, Y, p) where {FT}
             Ref(density),
         )
     Δz = snowfall ./ ρ_newsnow
-    if sum(isnan.(parent(Δz))) > 0
-        print("ΔZ NAN: ", snowfall, ρ_newsnow, "\n")
-        error()
-    end
 
     #Estimate resulting change from compaction effects of old snow (and the incoming snow on top):
     W_ice = (1 .- p.snow.q_l) .* Y.snow.S
@@ -530,26 +517,6 @@ function dzdt(density::Anderson1976{FT}, model::SnowModel{FT}, Y, p) where {FT}
     parent(ρ_ice)[parent(W_ice) .<= eps(FT)] .= FT(0.1) #any nonzero value for no-snowpack to avoid NaN and return 0.0
     ρ_est = compact_density.(W_use, sliq, ρ_ice, Δt_t, Tsnow, Ref(density))
     dz_compaction = (W_ice ./ ρ_est) .- Y.snow.Z
-    if sum(isnan.(parent(dz_compaction))) > 0
-        print(
-            "dz_compaction NAN: \n",
-            ρ_est,
-            "\n",
-            W_ice,
-            "\n",
-            W_use,
-            "\n",
-            sliq,
-            "\n",
-            ρ_ice,
-            "\n",
-            Y.snow.Z,
-            "\n",
-            Y.snow.S,
-            "\n***\n",
-        )
-        error()
-    end
 
     return (dz_compaction .+ Δz) ./ model.parameters.Δt
 end
@@ -583,10 +550,4 @@ function update_density_prog!(density::Anderson1976, model::SnowModel, dY, Y, p)
             dzdt(density, model, Y, p),
             model.parameters.Δt,
         )
-    newz = parent(dY.snow.Z)[1] * model.parameters.Δt + parent(Y.snow.Z)[1]
-    newswe = parent(dY.snow.S)[1] * model.parameters.Δt + parent(Y.snow.S)[1]
-    #if newz < newswe
-    #    print("BAD Z < SWE\n", dY.snow.Z, "\n", dY.snow.S, "\n", Y.snow.Z, "\n", Y.snow.S, "\n", dzdt(density, model, Y, p), "\n", newz, " | ", newswe)
-    #    error()
-    #end
 end