Fix NaNs in nonlinear solver

[simone-silvestri]

When h = 0, the internal flux is zero and does not depend on the input temperature.
For this reason, a root solver cannot find a solution for T, which ends up being NaN.

[glwagner]

Are we sure this is valid for all kinds of thermodynamics models we might implement?

[glwagner]

I think maybe rather than putting this check in, we should avoid solving for surface temperature if the thickness is 0. How about that?

[glwagner]

But we already have

ClimaSeaIce.jl/src/SeaIceThermodynamics/slab_thermodynamics_tendencies.jl

Line 29 in 2dee0f4

consolidated_ice = hᵢ >= hᶜ

and

ClimaSeaIce.jl/src/SeaIceThermodynamics/slab_thermodynamics_tendencies.jl

Lines 35 to 40 in 2dee0f4

	if consolidated_ice # slab is consolidated and has an independent surface temperature
	Tu⁻ = @inbounds Tu[i, j, k]
	Tuⁿ = top_surface_temperature(i, j, grid, top_heat_bc, Tu⁻, Qi, Qu, clock, model_fields)
	else # slab is unconsolidated and does not have an independent surface temperature
	Tuⁿ = bottom_temperature(i, j, grid, bottom_heat_bc, liquidus)
	end

Can you start by adding a test that shows what the problem is?