Support Robin boundary conditions for finite difference derivative operators

**Is your feature request related to a problem? Please describe.**
The finite difference derivative operators support `SetValue` and `SetGradient` boundary operators.
These allow to solve PDE problems with Neumann and Dirichlet boundary conditions (BC), see for example [this tutorial](https://clima.github.io/ClimaCore.jl/dev/tutorials/introduction/#3.1-Heat-equation-using-finite-differences).
A third, standard boundary condition for PDEs is the Robin boundary condition; a linear combination of Neumann and Dirichlet BCs given by:
$\left.(a \nabla u \cdot n + b u)\right|_{\partial\Omega} = g$,
with $a,b$ some given constants, $u$ our unknown on $\Omega$, $n$ the normal vector, $\partial\Omega$ the boundary of the domain, and $g$ a function on said boundary.

Right now, as far as I understand it is not possible (or at least not documented and not obvious to me) how one could enforce such a boundary condition.
Robin-type boundary conditions appear in climate modeling, especially in simplified models suitable for testing and analysis.[^1]

**Describe the solution you'd like**
It should be possible to provide such a boundary condition for a diffusion-like PDE.
Since it consists of a combination of Dirichlet and Neumann BC, and since Dirichlet BC have to be implemented via `SetValue` for gradient operators, I think it makes sense to implement this for the gradient operators.
I personally need this for the combination `GradientC2F` x `DivergenceF2C`.

With the stencils used, for the lower boundary this will take the form:
```
a*(u[1] - u[0]) + b u[0] = g
a*u[1] + (b-a) * u[0] = g
=> u[0] = (g - a*u[1]) / (b-a)
```
and thus can be implemented as
```
G(1/2) = u[1] - (g - a*u[1]) / (b-a) (1)
```

Note that the Dirichlet/Neumann cases are recovered for `a=0,b=1`, `a=1,b=0`, respectively.

**Describe alternatives you've considered**
These are not really alternatives but here are some thoughts I have on this and what make it a bit complicated to get running:
- Implementing this as `SetRobin` would deviate in spirit from the existing finite difference boundary conditions which use the terms "value" and "gradient" instead of "Dirichlet" and "Neumann". `SetLinearCombination` might be an alternative but it is not obvious to understand
- `SetGradient` expects a vector, `SetValue` a scalar. This is part of why I don't think it is obvious how I could implement a Robin-type BC with the current feature set. Note that Neumann BC and Robin BC typically explicitly set not the gradient but $\nabla u \cdot n$, which, again, is a scalar. In the finite difference framework in ClimaCoupler, as far as I understand from #202, one implicitly assumes to live in the vertical dimension, so $n=(0,0,\pm 1)^T$.
- Note that $b-a=0$ for $b=a$. It might be reasonable to rescale $\tilde{g}=g/a, \tilde{b}=b/a, \tilde{a}=1$ to prevent division by zero in (1).

[^1]: Examples: Bulk conditions at the atmosphere-ocean interface (e.g., Sockwell et al., [2020](https://doi.org/10.1016/j.rinam.2020.100110)) or the surface boundary condition for vertical sea ice thermodynamics (e.g., Maykut & Untersteiner, [1971](https://doi.org/10.1029/JC076i006p01550))

Start by reading the existing finite difference derivative operators and their SetValue and SetGradient boundary-condition APIs, then review the linked heat-equation tutorial. Work through the proposed Robin formulation for GradientC2F with DivergenceF2C, including the b=a case, and consider how the API should represent scalar boundary data. Done means Robin conditions can be applied to a diffusion-like PDE while retaining the existing Dirichlet and Neumann cases.

Written by the indexing model from the issue text.