DivergenceC2F Boundary Condition Bug

## Describe the bug
When constructing a `DivergenceC2F` operator, if an entire field is given for the boundary condition instead of a single column, the operator will just use the first value in the field instead of failing outright.

This was found in ClimaAtmos.jl, where the following code did not return an error, but just used the first column regardless of the set `colidx`
```
ᶜdivᵥ_ρe_tot = Operators.DivergenceF2C(
top = Operators.SetValue(CT3(FT(0))),
bottom = Operators.SetValue(sfc_conditions.ρ_flux_h_tot),
)
@. Yₜ.c.ρe_tot[colidx] -= ᶜdivᵥ_ρe_tot(
-(
ᶠρK_E[colidx] *
ᶠgradᵥ((Y.c.ρe_tot[colidx] + ᶜp[colidx]) / Y.c.ρ[colidx])
),
)
```
The actual operator should have been the following:
```
ᶜdivᵥ_ρe_tot = Operators.DivergenceF2C(
top = Operators.SetValue(CT3(FT(0))),
bottom = Operators.SetValue(sfc_conditions.ρ_flux_h_tot[colidx]),
)
```

@trontrytel I am happy to talk this over more, I wasn't able to reproduce this bug without running Atmos

Start by tracing construction and boundary-value handling in the DivergenceC2F operator, using the ClimaAtmos.jl example in the issue as the reproduction context. Confirm that passing a full field where one column is required fails instead of silently using the first value; done means the invalid input is rejected while the indexed form remains valid.

Written by the indexing model from the issue text.