Update Simple FCT: lin combination of UpwindBiasedProductC2F + Upwind3rdOrderBiasedProductC2F convergence test to use same order approximation

I spoke with @akshaysridhar today about this [FCT test](https://github.com/CliMA/ClimaCore.jl/blob/a932c51ac0bf780586539629de2c0cf7a411bb5b/test/Operators/finitedifference/convergence_column.jl#L691-L760), and plotted the computed vs exact solutions:

Here's the full modified script in case it's helpful:

```julia
#=
julia --check-bounds=yes --project=.buildkite
using Revise; include("../cc_fd_biased_prod.jl")
=#
using Test
using ClimaCorePlots
using LazyBroadcast: @lazy
using Plots
using StaticArrays, IntervalSets, LinearAlgebra

using ClimaComms
ClimaComms.@import_required_backends
import ClimaCore: slab, Domains, Meshes, Topologies, Spaces, Fields, Operators
import ClimaCore.Domains: Geometry
import ClimaCore.DataLayouts: vindex

device = ClimaComms.device()

function stencil_location(space, bc, idx)
(li, lw, rw, ri) = Operators.window_bounds(space, bc)
if !Topologies.isperiodic(Spaces.vertical_topology(space))
if idx in li:(lw - 1)
return Operators.LeftBoundaryWindow{Spaces.left_boundary_name(space)}()
end
end
if idx in lw:rw
return Interior()
end
if !Topologies.isperiodic(Spaces.vertical_topology(space))
if idx in (rw + 1):ri
return Operators.RightBoundaryWindow{Spaces.right_boundary_name(space)}()
end
end
error("Uncaught case")
end

"""
convergence_rate(err, Δh)

Estimate convergence rate given vectors `err` and `Δh`

err = C Δh^p+ H.O.T
err_k ≈ C Δh_k^p
err_k/err_m ≈ Δh_k^p/Δh_m^p
log(err_k/err_m) ≈ log((Δh_k/Δh_m)^p)
log(err_k/err_m) ≈ p*log(Δh_k/Δh_m)
log(err_k/err_m)/log(Δh_k/Δh_m) ≈ p

"""
convergence_rate(err, Δh) =
[log(err[i] / err[i - 1]) / log(Δh[i] / Δh[i - 1]) for i in 2:length(Δh)]

FT = Float64
n_elems_seq = 2 .^ (4, 6, 8, 10)
stretch_fns = (Meshes.Uniform(), Meshes.ExponentialStretching(1.0))
device = ClimaComms.device()

i = 1
stretch_fn = stretch_fns[1]
err_adv_wc = zeros(FT, length(n_elems_seq))
Δh = zeros(FT, length(n_elems_seq))
Plots.plot()
for (k, n) in enumerate(n_elems_seq)
domain = Domains.IntervalDomain(
Geometry.ZPoint{FT}(-pi),
Geometry.ZPoint{FT}(pi);
boundary_names = (:bottom, :top),
)
mesh = Meshes.IntervalMesh(domain, stretch_fn; nelems = n)

cs = Spaces.CenterFiniteDifferenceSpace(device, mesh)
fs = Spaces.FaceFiniteDifferenceSpace(cs)

zc = Fields.coordinate_field(cs).z
C = FT(1.0) # flux-correction coefficient (falling back to third-order upwinding)

# UpwindBiasedProductC2F & Upwind3rdOrderBiasedProductC2F Center -> Face operator
# Unitary, constant advective velocity
w = Geometry.WVector.(ones(fs))
# parent(w) .= parent(w) * -1
Δz = FT(2pi / n)
c = sin.(zc)

zcontra3 = Geometry.Contravariant3Vector(FT(0.0))

first_order_fluxᶠ = Operators.UpwindBiasedProductC2F(
bottom = Operators.Extrapolate(),
top = Operators.Extrapolate(),
)
third_order_fluxᶠ = Operators.Upwind3rdOrderBiasedProductC2F(
bottom = Operators.FirstOrderOneSided(),
top = Operators.FirstOrderOneSided(),
)

divf2c = Operators.DivergenceF2C(
bottom = Operators.SetValue(
zcontra3,
),
top = Operators.SetValue(
zcontra3,
),
)

adv_wc = @lazy @. divf2c.(third_order_fluxᶠ(w, c))
zero_bc = Operators.SetBoundaryOperator(;
bottom = Operators.SetValue(zcontra3), top = Operators.SetValue(zcontra3)
)
zbc = zero_bc.(third_order_fluxᶠ.(w, c)).components.data.:1
zbc .= zbc ./ maximum(zbc) .* maximum(c)
Δh[k] = Spaces.local_geometry_data(fs).J[vindex(1)]

exact = cos.(zc)
err = parent(Fields.field_values(abs.(Base.materialize(adv_wc) .- exact)))

if k==1
Plots.plot!(cos.(zc); label = "exact")
end
Plots.plot!(Base.materialize(adv_wc); label = "computed, dh=$(round(Δh[k]; digits=4))")
Plots.plot!(c; label = "c, dh=$(round(Δh[k]; digits=4))")
Plots.plot!(zbc; label = "zbc, dh=$(round(Δh[k]; digits=4))")

# Error
err_adv_wc[k] = norm(Base.materialize(adv_wc) .- cos.(zc))
end
Plots.plot!(; legend = :topleft)
Plots.png("convergence_ubp.png")
```

Start with test/Operators/finitedifference/convergence_column.jl at lines 691-760 and compare it with the modified script in the issue. Run the FCT convergence test, inspect the computed versus exact solutions and convergence behavior, then update the test so the linear combination uses the same order approximation. Done means the test reflects the intended approximation and its convergence results are consistent.

Written by the indexing model from the issue text.