integrate neural network into Snow.jl architecture

[a-charbon]

Purpose

Integrate the neural snow depth model into the Snow model architecture. This involves small changes to Snow.jl and SnowParameterizations.jl as well as expanding the NeuralSnow extension to allow the neural model as a type of AbstractDensityModel. "CSV" was also added to the buildkite dependencies to permit testing of the NeuralSnow extension, and BSON has been appropriately added to the ClimaLand weak dependencies to enable the expanded functionality of the extension.

To-do

Questions are commented in the changed files about different possible approaches. With this PR merged, I can also add a slightly faster and more abstract/streamlined form of the neural model itself (changes to ModelTools, etc).

Content

Expanded NeuralSnow extension and changed some of the syntax of the Snow model.

---

• I have read and checked the items on the review checklist.

[kmdeck]

you dont need to set this, it is set in set_initial_cache!

[a-charbon]

Note to myself to check a model simulation without this line, I think I ran into problems without it

[kmdeck]

our style is more to have the download link download the file, managed by ClimaArtifacts, and then to load from file.

Perhaps you could switch to that style (read from local file), and open an issue in ClimaArtifacts to add this? we can resolve that later, though

[kmdeck]

we have neural snow model extension tests currently:
https://github.com/CliMA/ClimaLand.jl/blob/main/test/standalone/Snow/tool_tests.jl

it would be good to add some tests for these functions

[kmdeck]

yes, anything with a dot (broadcast) over a field needs to be in an expression that involves x .= ...

so that there are no allocations.

You can do this like:

    @. dY.snow.P_avg = β * (abs(p.drivers.P_snow) .- Y.snow.P_avg)
    @. dY.snow.T_avg = β * (p.drivers.T - model.parameters.earth_param_set.T_freeze - Y.snow.T_avg)
    @. dY.snow.R_avg = β * (p.drivers.SW_d - Y.snow.R_avg) # SW_d is positive by definition
    @. dY.snow.Qrel_avg = β * (Thermodynamics.relative_humidity(
            Ref(model.boundary_conditions.atmos.thermo_params),
            p.drivers.thermal_state,
        ) - Y.snow.Qrel_avg)
    @. dY.snow.u_avg = β * (wind_speed - Y.snow.u_avg)

the way it is written, it is not possible to do this for dY.snow.Z. But you could add in an aux var

 p.snow.predicted_dzdt .=  dzdt(density, Y)

which is updated here (if you like) and then compute

@. dY.snow.Z = clip_dZdt(swe, z, dswedt, predicted_dzdt, dt)