review feedback

experiments/AMIP/components/atmosphere/climaatmos.jl

@@ -1,10 +1,10 @@
 # atmos_init: for ClimaAtmos interface
-import StaticArrays as SA
-import Statistics as Stats
-import LinearAlgebra as LinAlg
+import StaticArrays
+import Statistics
+import LinearAlgebra
 
 import ClimaAtmos as CA
-import ClimaCore
+import ClimaCore as CC
 import SurfaceFluxes as SF
 
 import ClimaCoupler.Interfacer: AtmosModelSimulation
@@ -59,9 +59,9 @@ function atmos_init(::Type{FT}, atmos_config_dict::Dict) where {FT}
     ᶜ3d_snow = integrator.p.precipitation.ᶜ3d_snow
 
     # set initial fluxes to zero
-    @. ρ_flux_h_tot = ClimaCore.Geometry.Covariant3Vector(FT(0.0))
-    @. ρ_flux_q_tot = ClimaCore.Geometry.Covariant3Vector(FT(0.0))
-    @. ᶠradiation_flux = ClimaCore.Geometry.WVector(FT(0))
+    @. ρ_flux_h_tot = CC.Geometry.Covariant3Vector(FT(0.0))
+    @. ρ_flux_q_tot = CC.Geometry.Covariant3Vector(FT(0.0))
+    @. ᶠradiation_flux = CC.Geometry.WVector(FT(0))
     ρ_flux_uₕ.components .= Ref(SMatrix{1, 2}([FT(0), FT(0)]))
     col_integrated_rain .= FT(0)
     col_integrated_snow .= FT(0)
@@ -96,27 +96,19 @@ function get_field(atmos_sim::ClimaAtmosSimulation, ::Val{:radiative_energy_flux
 
     if atmos_sim.integrator.p.radiation.radiation_model != nothing
         face_space = axes(atmos_sim.integrator.u.f)
-        nz_faces = length(ClimaCore.Spaces.vertical_topology(face_space).mesh.faces)
+        nz_faces = length(CC.Spaces.vertical_topology(face_space).mesh.faces)
 
         (; face_lw_flux_dn, face_lw_flux_up, face_sw_flux_dn, face_sw_flux_up) =
             atmos_sim.integrator.p.radiation.radiation_model
 
-        LWd_TOA = ClimaCore.Fields.level(
-            CA.RRTMGPI.array2field(FT.(face_lw_flux_dn), face_space),
-            nz_faces - ClimaCore.Utilities.half,
-        )
-        LWu_TOA = ClimaCore.Fields.level(
-            CA.RRTMGPI.array2field(FT.(face_lw_flux_up), face_space),
-            nz_faces - ClimaCore.Utilities.half,
-        )
-        SWd_TOA = ClimaCore.Fields.level(
-            CA.RRTMGPI.array2field(FT.(face_sw_flux_dn), face_space),
-            nz_faces - ClimaCore.Utilities.half,
-        )
-        SWu_TOA = ClimaCore.Fields.level(
-            CA.RRTMGPI.array2field(FT.(face_sw_flux_up), face_space),
-            nz_faces - ClimaCore.Utilities.half,
-        )
+        LWd_TOA =
+            CC.Fields.level(CA.RRTMGPI.array2field(FT.(face_lw_flux_dn), face_space), nz_faces - CC.Utilities.half)
+        LWu_TOA =
+            CC.Fields.level(CA.RRTMGPI.array2field(FT.(face_lw_flux_up), face_space), nz_faces - CC.Utilities.half)
+        SWd_TOA =
+            CC.Fields.level(CA.RRTMGPI.array2field(FT.(face_sw_flux_dn), face_space), nz_faces - CC.Utilities.half)
+        SWu_TOA =
+            CC.Fields.level(CA.RRTMGPI.array2field(FT.(face_sw_flux_up), face_space), nz_faces - CC.Utilities.half)
 
         return @. -(LWd_TOA + SWd_TOA - LWu_TOA - SWu_TOA)
     else
@@ -151,32 +143,31 @@ get_field(sim::ClimaAtmosSimulation, ::Val{:air_temperature}) =
     TD.air_temperature.(thermo_params, sim.integrator.p.precomputed.ᶜts)
 get_field(sim::ClimaAtmosSimulation, ::Val{:liquid_precipitation}) = sim.integrator.p.precipitation.col_integrated_rain
 get_field(sim::ClimaAtmosSimulation, ::Val{:radiative_energy_flux_sfc}) =
-    ClimaCore.Fields.level(sim.integrator.p.radiation.ᶠradiation_flux, ClimaCore.Utilities.half)
+    CC.Fields.level(sim.integrator.p.radiation.ᶠradiation_flux, CC.Utilities.half)
 get_field(sim::ClimaAtmosSimulation, ::Val{:snow_precipitation}) = sim.integrator.p.precipitation.col_integrated_snow
 get_field(sim::ClimaAtmosSimulation, ::Val{:turbulent_energy_flux}) =
-    ClimaCore.Geometry.WVector.(sim.integrator.p.precomputed.sfc_conditions.ρ_flux_h_tot)
+    CC.Geometry.WVector.(sim.integrator.p.precomputed.sfc_conditions.ρ_flux_h_tot)
 get_field(sim::ClimaAtmosSimulation, ::Val{:turbulent_moisture_flux}) =
-    ClimaCore.Geometry.WVector.(sim.integrator.p.precomputed.sfc_conditions.ρ_flux_q_tot)
-get_field(sim::ClimaAtmosSimulation, ::Val{:thermo_state_int}) =
-    ClimaCore.Spaces.level(sim.integrator.p.precomputed.ᶜts, 1)
+    CC.Geometry.WVector.(sim.integrator.p.precomputed.sfc_conditions.ρ_flux_q_tot)
+get_field(sim::ClimaAtmosSimulation, ::Val{:thermo_state_int}) = CC.Spaces.level(sim.integrator.p.precomputed.ᶜts, 1)
 get_field(atmos_sim::ClimaAtmosSimulation, ::Val{:water}) = atmos_sim.integrator.u.c.ρq_tot
 
 # extensions required by FluxCalculator (partitioned fluxes)
 get_field(sim::ClimaAtmosSimulation, ::Val{:height_int}) =
-    ClimaCore.Spaces.level(ClimaCore.Fields.coordinate_field(sim.integrator.u.c).z, 1)
+    CC.Spaces.level(CC.Fields.coordinate_field(sim.integrator.u.c).z, 1)
 get_field(sim::ClimaAtmosSimulation, ::Val{:height_sfc}) =
-    ClimaCore.Spaces.level(ClimaCore.Fields.coordinate_field(sim.integrator.u.f).z, ClimaCore.Utilities.half)
+    CC.Spaces.level(CC.Fields.coordinate_field(sim.integrator.u.f).z, CC.Utilities.half)
 function get_field(sim::ClimaAtmosSimulation, ::Val{:uv_int})
-    uₕ_int = ClimaCore.Geometry.UVVector.(ClimaCore.Spaces.level(sim.integrator.u.c.uₕ, 1))
-    return @. SA.SVector(uₕ_int.components.data.:1, uₕ_int.components.data.:2)
+    uₕ_int = CC.Geometry.UVVector.(CC.Spaces.level(sim.integrator.u.c.uₕ, 1))
+    return @. StaticArrays.SVector(uₕ_int.components.data.:1, uₕ_int.components.data.:2)
 end
 
 function update_field!(atmos_sim::ClimaAtmosSimulation, ::Val{:co2}, field)
     if atmos_sim.integrator.p.atmos.radiation_mode isa CA.RRTMGPI.GrayRadiation
         @warn("Gray radiation model initialized, skipping CO2 update", maxlog = 1)
         return
     else
-        atmos_sim.integrator.p.radiation.radiation_model.volume_mixing_ratio_co2 .= Stats.mean(parent(field))
+        atmos_sim.integrator.p.radiation.radiation_model.volume_mixing_ratio_co2 .= Statistics.mean(parent(field))
     end
 end
 # extensions required by the Interfacer
@@ -198,7 +189,7 @@ function update_field!(sim::ClimaAtmosSimulation, ::Val{:turbulent_fluxes}, fiel
     (; F_turb_energy, F_turb_moisture, F_turb_ρτxz, F_turb_ρτyz) = fields
 
     Y = sim.integrator.u
-    surface_local_geometry = ClimaCore.Fields.level(ClimaCore.Fields.local_geometry_field(Y.f), ClimaCore.Fields.half)
+    surface_local_geometry = CC.Fields.level(CC.Fields.local_geometry_field(Y.f), CC.Fields.half)
     surface_normal = @. CA.C3(CA.unit_basis_vector_data(CA.C3, surface_local_geometry))
 
     # get template objects for the contravariant components of the momentum fluxes (required by Atmos boundary conditions)
@@ -236,11 +227,11 @@ function update_sim!(atmos_sim::ClimaAtmosSimulation, csf, turbulent_fluxes)
 end
 
 """
-    calculate_surface_air_density(atmos_sim::ClimaAtmosSimulation, T_S::ClimaCore.Fields.Field)
+    calculate_surface_air_density(atmos_sim::ClimaAtmosSimulation, T_S::CC.Fields.Field)
 
 Extension for this function to calculate surface density.
 """
-function calculate_surface_air_density(atmos_sim::ClimaAtmosSimulation, T_S::ClimaCore.Fields.Field)
+function calculate_surface_air_density(atmos_sim::ClimaAtmosSimulation, T_S::CC.Fields.Field)
     thermo_params = get_thermo_params(atmos_sim)
     ts_int = get_field(atmos_sim, Val(:thermo_state_int))
     extrapolate_ρ_to_sfc.(Ref(thermo_params), ts_int, swap_space!(ones(axes(ts_int.ρ)), T_S))
@@ -309,7 +300,7 @@ struct CoupledMoninObukhov end
 """
     coupler_surface_setup(::CoupledMoninObukhov, p, csf_sfc = (; T = nothing, z0m = nothing, z0b = nothing, beta = nothing, q_vap = nothing))
 
-Sets up `surface_setup` as a `ClimaCore.Fields.Field` of `SurfaceState`s.
+Sets up `surface_setup` as a `CC.Fields.Field` of `SurfaceState`s.
 """
 function coupler_surface_setup(
     ::CoupledMoninObukhov,
@@ -397,20 +388,20 @@ function dss_state!(sim::ClimaAtmosSimulation)
     Y = sim.integrator.u
     for key in propertynames(Y)
         field = getproperty(Y, key)
-        buffer = ClimaCore.Spaces.create_dss_buffer(field)
-        ClimaCore.Spaces.weighted_dss!(field, buffer)
+        buffer = CC.Spaces.create_dss_buffer(field)
+        CC.Spaces.weighted_dss!(field, buffer)
     end
 end
 
 """
-    water_albedo_from_atmosphere!(atmos_sim::ClimaAtmosSimulation, direct_albedo::ClimaCore.Fields.Field, diffuse_albedo::ClimaCore.Fields.Field)
+    water_albedo_from_atmosphere!(atmos_sim::ClimaAtmosSimulation, direct_albedo::CC.Fields.Field, diffuse_albedo::CC.Fields.Field)
 
 Extension to calculate the water surface albedo from wind speed and insolation. It can be used for prescribed ocean and lakes.
 """
 function water_albedo_from_atmosphere!(
     atmos_sim::ClimaAtmosSimulation,
-    direct_albedo::ClimaCore.Fields.Field,
-    diffuse_albedo::ClimaCore.Fields.Field,
+    direct_albedo::CC.Fields.Field,
+    diffuse_albedo::CC.Fields.Field,
 )
 
     Y = atmos_sim.integrator.u
@@ -424,14 +415,14 @@ function water_albedo_from_atmosphere!(
     # update for current zenith angle
     CA.set_insolation_variables!(Y, p, t)
 
-    bottom_coords = ClimaCore.Fields.coordinate_field(ClimaCore.Spaces.level(Y.c, 1))
+    bottom_coords = CC.Fields.coordinate_field(CC.Spaces.level(Y.c, 1))
     μ = CA.RRTMGPI.array2field(radiation_model.cos_zenith, axes(bottom_coords))
     FT = eltype(atmos_sim.integrator.u)
     α_model = CA.RegressionFunctionAlbedo{FT}()
 
 
     # set the direct and diffuse surface albedos
-    direct_albedo .= CA.surface_albedo_direct(α_model).(λ, μ, LinAlg.norm.(ClimaCore.Fields.level(Y.c.uₕ, 1)))
-    diffuse_albedo .= CA.surface_albedo_diffuse(α_model).(λ, μ, LinAlg.norm.(ClimaCore.Fields.level(Y.c.uₕ, 1)))
+    direct_albedo .= CA.surface_albedo_direct(α_model).(λ, μ, LinearAlgebra.norm.(CC.Fields.level(Y.c.uₕ, 1)))
+    diffuse_albedo .= CA.surface_albedo_diffuse(α_model).(λ, μ, LinearAlgebra.norm.(CC.Fields.level(Y.c.uₕ, 1)))
 
 end

experiments/AMIP/components/atmosphere/climaatmos_extra_diags.jl

@@ -23,7 +23,7 @@ add_diagnostic_variable!(
         (; ᶜts) = cache.precomputed
         c_space = axes(state.c)
         thermo_params = CAP.thermodynamics_params(params)
-        e_pot = CAP.grav(params) .* ClimaCore.Fields.coordinate_field(c_space).z
+        e_pot = CAP.grav(params) .* CC.Fields.coordinate_field(c_space).z
         if isnothing(out)
             return TD.moist_static_energy.(thermo_params, ᶜts, e_pot)
         else
@@ -44,9 +44,9 @@ add_diagnostic_variable!(
         thermo_params = CAP.thermodynamics_params(params)
         ᶜT = @. TD.air_temperature(thermo_params, ᶜts)
         if isnothing(out)
-            return ClimaCore.Geometry.WVector.(CAD.ᶜgradᵥ.(CAD.ᶠinterp.(ᶜT))).components.data.:1
+            return CC.Geometry.WVector.(CAD.ᶜgradᵥ.(CAD.ᶠinterp.(ᶜT))).components.data.:1
         else
-            out .= ClimaCore.Geometry.WVector.(CAD.ᶜgradᵥ.(CAD.ᶠinterp.(ᶜT))).components.data.:1
+            out .= CC.Geometry.WVector.(CAD.ᶜgradᵥ.(CAD.ᶠinterp.(ᶜT))).components.data.:1
         end
 
     end,
@@ -61,11 +61,11 @@ add_diagnostic_variable!(
     compute! = (out, state, cache, time) -> begin
         (; ᶜp) = cache.precomputed
         (; C_E) = cache.atmos.vert_diff
-        interior_uₕ = ClimaCore.Fields.level(state.c.uₕ, 1)
+        interior_uₕ = CC.Fields.level(state.c.uₕ, 1)
         ᶠp = ᶠK_E = cache.scratch.ᶠtemp_scalar
-        ClimaCore.Fields.bycolumn(axes(ᶜp)) do colidx
+        CC.Fields.bycolumn(axes(ᶜp)) do colidx
             @. ᶠp[colidx] = CAD.ᶠinterp(ᶜp[colidx])
-            ᶜΔz_surface = ClimaCore.Fields.Δz_field(interior_uₕ)
+            ᶜΔz_surface = CC.Fields.Δz_field(interior_uₕ)
             @. ᶠK_E[colidx] = CA.eddy_diffusivity_coefficient(
                 C_E,
                 CA.norm(interior_uₕ[colidx]),

experiments/AMIP/components/land/climaland_bucket.jl

@@ -1,12 +1,12 @@
 import Dates
 import Thermodynamics as TD
 
-import ClimaCore
+import ClimaCore as CC
 import ClimaTimeSteppers as CTS
 import ClimaParams
 
 import ClimaLand as CL
-import CL.Bucket:
+import ClimaLand.Bucket:
     BucketModel,
     BucketModelParameters,
     AbstractAtmosphericDrivers,
@@ -158,7 +158,7 @@ Extension of Interfacer.get_field that provides the total energy contained in th
 function get_field(bucket_sim::BucketSimulation, ::Val{:energy})
     # required by ConservationChecker
     e_per_area = zeros(axes(bucket_sim.integrator.u.bucket.W))
-    ClimaCore.Fields.bycolumn(axes(bucket_sim.integrator.u.bucket.T)) do colidx
+    CC.Fields.bycolumn(axes(bucket_sim.integrator.u.bucket.T)) do colidx
         e_per_area[colidx] .=
             bucket_sim.model.parameters.ρc_soil .* mean(bucket_sim.integrator.u.bucket.T[colidx]) .*
             bucket_sim.domain.soil_depth
@@ -209,7 +209,7 @@ step!(sim::BucketSimulation, t) = step!(sim.integrator, t - sim.integrator.t, tr
 reinit!(sim::BucketSimulation) = reinit!(sim.integrator)
 
 # extensions required by FluxCalculator (partitioned fluxes)
-function update_turbulent_fluxes_point!(sim::BucketSimulation, fields::NamedTuple, colidx::ClimaCore.Fields.ColumnIndex)
+function update_turbulent_fluxes_point!(sim::BucketSimulation, fields::NamedTuple, colidx::CC.Fields.ColumnIndex)
     (; F_turb_energy, F_turb_moisture) = fields
     turbulent_fluxes = sim.integrator.p.bucket.turbulent_fluxes
     turbulent_fluxes.shf[colidx] .= F_turb_energy
@@ -223,7 +223,7 @@ function surface_thermo_state(
     sim::BucketSimulation,
     thermo_params::TD.Parameters.ThermodynamicsParameters,
     thermo_state_int,
-    colidx::ClimaCore.Fields.ColumnIndex,
+    colidx::CC.Fields.ColumnIndex,
 )
 
     T_sfc = get_field(sim, Val(:surface_temperature), colidx)
@@ -290,38 +290,37 @@ temp_anomaly_amip(coord) = 40 * cosd(coord.lat)^4
 
 """
     make_land_domain(
-        atmos_boundary_space::ClimaCore.Spaces.SpectralElementSpace2D,
+        atmos_boundary_space::CC.Spaces.SpectralElementSpace2D,
         zlim::Tuple{FT, FT},
         nelements_vert::Int,) where {FT}
 
 Creates the land model domain from the horizontal space of the atmosphere, and information
 about the number of elements and extent of the vertical domain.
 """
 function make_land_domain(
-    atmos_boundary_space::ClimaCore.Spaces.SpectralElementSpace2D,
+    atmos_boundary_space::CC.Spaces.SpectralElementSpace2D,
     zlim::Tuple{FT, FT},
     nelements_vert::Int,
 ) where {FT}
     @assert zlim[1] < zlim[2]
     depth = zlim[2] - zlim[1]
 
-    mesh = ClimaCore.Spaces.topology(atmos_boundary_space).mesh
+    mesh = CC.Spaces.topology(atmos_boundary_space).mesh
 
     radius = mesh.domain.radius
     nelements_horz = mesh.ne
-    npolynomial =
-        ClimaCore.Spaces.Quadratures.polynomial_degree(ClimaCore.Spaces.quadrature_style(atmos_boundary_space))
+    npolynomial = CC.Spaces.Quadratures.polynomial_degree(CC.Spaces.quadrature_style(atmos_boundary_space))
     nelements = (nelements_horz, nelements_vert)
-    vertdomain = ClimaCore.Domains.IntervalDomain(
-        ClimaCore.Geometry.ZPoint(FT(zlim[1])),
-        ClimaCore.Geometry.ZPoint(FT(zlim[2]));
+    vertdomain = CC.Domains.IntervalDomain(
+        CC.Geometry.ZPoint(FT(zlim[1])),
+        CC.Geometry.ZPoint(FT(zlim[2]));
         boundary_names = (:bottom, :top),
     )
 
-    vertmesh = ClimaCore.Meshes.IntervalMesh(vertdomain, ClimaCore.Meshes.Uniform(), nelems = nelements[2])
-    verttopology = ClimaCore.Topologies.IntervalTopology(vertmesh)
-    vert_center_space = ClimaCore.Spaces.CenterFiniteDifferenceSpace(verttopology)
-    subsurface_space = ClimaCore.Spaces.ExtrudedFiniteDifferenceSpace(atmos_boundary_space, vert_center_space)
+    vertmesh = CC.Meshes.IntervalMesh(vertdomain, CC.Meshes.Uniform(), nelems = nelements[2])
+    verttopology = CC.Topologies.IntervalTopology(vertmesh)
+    vert_center_space = CC.Spaces.CenterFiniteDifferenceSpace(verttopology)
+    subsurface_space = CC.Spaces.ExtrudedFiniteDifferenceSpace(atmos_boundary_space, vert_center_space)
     space = (; surface = atmos_boundary_space, subsurface = subsurface_space)
 
     return CL.Domains.SphericalShell{FT}(radius, depth, nothing, nelements, npolynomial, space)

experiments/AMIP/components/ocean/eisenman_seaice.jl

@@ -1,6 +1,6 @@
 import SciMLBase
 
-import ClimaCore
+import ClimaCore as CC
 import ClimaTimeSteppers as CTS
 
 import ClimaCoupler: FluxCalculator
@@ -79,7 +79,7 @@ function eisenman_seaice_init(
         area_fraction = area_fraction,
         ice_area_fraction = zeros(space),
         thermo_params = thermo_params,
-        dss_buffer = ClimaCore.Spaces.create_dss_buffer(ClimaCore.Fields.zeros(space)),
+        dss_buffer = CC.Spaces.create_dss_buffer(CC.Fields.zeros(space)),
     )
     problem = SciMLBase.ODEProblem(ode_function, Y, Float64.(tspan), cache)
     integrator = SciMLBase.init(problem, ode_algo, dt = Float64(dt), saveat = Float64(saveat), adaptive = false)
@@ -142,7 +142,7 @@ end
 function update_field!(sim::EisenmanIceSimulation, ::Val{:air_density}, field)
     parent(sim.integrator.p.Ya.ρ_sfc) .= parent(field)
 end
-function update_field!(sim::EisenmanIceSimulation, ::Val{:area_fraction}, field::ClimaCore.Fields.Field)
+function update_field!(sim::EisenmanIceSimulation, ::Val{:area_fraction}, field::CC.Fields.Field)
     sim.integrator.p.area_fraction .= field
 end
 function update_field!(sim::EisenmanIceSimulation, ::Val{:∂F_turb_energy∂T_sfc}, field, colidx)
@@ -160,11 +160,7 @@ step!(sim::EisenmanIceSimulation, t) = step!(sim.integrator, t - sim.integrator.
 reinit!(sim::EisenmanIceSimulation) = reinit!(sim.integrator)
 
 # extensions required by FluxCalculator (partitioned fluxes)
-function update_turbulent_fluxes_point!(
-    sim::EisenmanIceSimulation,
-    fields::NamedTuple,
-    colidx::ClimaCore.Fields.ColumnIndex,
-)
+function update_turbulent_fluxes_point!(sim::EisenmanIceSimulation, fields::NamedTuple, colidx::CC.Fields.ColumnIndex)
     (; F_turb_energy) = fields
     @. sim.integrator.p.Ya.F_turb[colidx] = F_turb_energy
 end
@@ -219,20 +215,20 @@ end
 
 Initialize the state vectors for the Eisenman-Zhang sea ice model.
 """
-function state_init(p::EisenmanIceParameters, space::ClimaCore.Spaces.AbstractSpace)
-    Y = ClimaCore.Fields.FieldVector(
+function state_init(p::EisenmanIceParameters, space::CC.Spaces.AbstractSpace)
+    Y = CC.Fields.FieldVector(
         T_sfc = ones(space) .* p.T_freeze,
         h_ice = zeros(space),
         T_ml = ones(space) .* 277,
-        q_sfc = ClimaCore.Fields.zeros(space),
+        q_sfc = CC.Fields.zeros(space),
     )
-    Ya = ClimaCore.Fields.FieldVector(
-        F_turb = ClimaCore.Fields.zeros(space),
-        ∂F_turb_energy∂T_sfc = ClimaCore.Fields.zeros(space),
-        F_rad = ClimaCore.Fields.zeros(space),
-        e_base = ClimaCore.Fields.zeros(space),
-        ocean_qflux = ClimaCore.Fields.zeros(space),
-        ρ_sfc = ClimaCore.Fields.zeros(space),
+    Ya = CC.Fields.FieldVector(
+        F_turb = CC.Fields.zeros(space),
+        ∂F_turb_energy∂T_sfc = CC.Fields.zeros(space),
+        F_rad = CC.Fields.zeros(space),
+        e_base = CC.Fields.zeros(space),
+        ocean_qflux = CC.Fields.zeros(space),
+        ρ_sfc = CC.Fields.zeros(space),
     )
     return Y, Ya
 end
@@ -352,7 +348,7 @@ function ∑tendencies(dY, Y, cache, _)
     @. dY.T_sfc = -Y.T_sfc / Δt
     @. dY.q_sfc = -Y.q_sfc / Δt
 
-    ClimaCore.Fields.bycolumn(axes(Y.T_sfc)) do colidx
+    CC.Fields.bycolumn(axes(Y.T_sfc)) do colidx
         solve_eisenman_model!(Y[colidx], Ya[colidx], p, thermo_params, Δt)
     end