diff --git a/src/FieldArray.jl b/src/FieldArray.jl
index 1edfec0..ce0b426 100644
--- a/src/FieldArray.jl
+++ b/src/FieldArray.jl
@@ -1,5 +1,10 @@
 """
-    FieldArray
+    FieldArray(
+        name::String, values::Array;
+        dims::Union{Nothing, Vector{<:AbstractString}, Vector{<:Pair}}=nothing,
+        coords::Vector{<:FieldArray}=FieldArray[],
+        attributes=nothing,
+    )
 
 A generic [xarray](https://xarray.pydata.org/en/stable/index.html)-like or 
 [IRIS](https://scitools-iris.readthedocs.io/en/latest/)-like 
@@ -7,25 +12,111 @@ Array with named dimensions and optional coordinates.
 
 NB: this aims to be simple and generic, not efficient !!! Intended for representing model output,
 not for numerically-intensive calculations.
+
+# Examples
+
+2-dimensional FieldArray with no supplied dimension names or coordinates:
+
+    julia> fa = PALEOmodel.FieldArray("vec", [1 2 3; 4 5 6])
+    FieldArray(name="vec", eltype=Int64)
+      dims: (NamedDimension(name=dim_1, size=2, coords=String[]), NamedDimension(name=dim_2, size=3, coords=String[]))
+
+2-dimensional FieldArray with supplied dimension names, no coordinates:
+
+    julia> fa = PALEOmodel.FieldArray("vec", [1 2 3; 4 5 6]; dims=["x", "y"])
+    FieldArray(name="vec", eltype=Int64)
+      dims: (NamedDimension(name=x, size=2, coords=String[]), NamedDimension(name=y, size=3, coords=String[]))
+
+2-dimensional FieldArray with supplied dimension names and coordinates:
+
+    julia> fa = PALEOmodel.FieldArray(
+                "vec", [1 2 3; 4 5 6]; 
+                dims=[
+                    "y"=>("ymid",), 
+                    "x"=>("xmid",)
+                ], 
+                coords=[
+                    PALEOmodel.FieldArray("xmid", [0.5, 1.5, 2.5], dims=["x"]), 
+                    PALEOmodel.FieldArray("ymid", [0.5, 1.5], dims=["y"])
+                ],
+            )
+
+    FieldArray(name="vec", eltype=Int64)
+      dims: (NamedDimension(name=y, size=2, coords=["ymid"]), NamedDimension(name=x, size=3, coords=["xmid"]))
+      coords:
+        FieldArray(name="xmid", eltype=Float64, dims=(NamedDimension(name=x, size=3, coords=String[]),))
+        FieldArray(name="ymid", eltype=Float64, dims=(NamedDimension(name=y, size=2, coords=String[]),))
+
 """
 struct FieldArray{T, N}
     name::String
     values::T
     dims::NTuple{N, PB.NamedDimension}
+    coords::Vector{FieldArray}
     attributes::Union{Dict, Nothing}
 end
 
 
+function FieldArray(
+    name::String, values::Array;
+    dims::Union{Nothing, Vector{<:AbstractString}, Vector{<:Pair}}=nothing,
+    coords::Vector{<:FieldArray}=FieldArray[],
+    attributes=nothing,
+)
+    if isnothing(dims) 
+        dims = [PB.NamedDimension("dim_$i", s, []) for (i, s) in enumerate(size(values))]
+        
+    elseif dims isa Vector{<:AbstractString}
+        tmpdims = []
+        for (d, s) in PB.IteratorUtils.zipstrict(dims, size(values); 
+                errmsg="values has $(length(size(values))) dimensions but $(length(dims)) dims supplied")
+            push!(tmpdims, PB.NamedDimension(d, s, []))
+        end
+        dims=tmpdims
+    elseif dims isa Vector{<:Pair}
+        tmpdims = []
+        for ((dname, dcoords), s) in PB.IteratorUtils.zipstrict(dims, size(values); 
+            errmsg="values has $(length(size(values))) dimensions but $(length(dims)) dims supplied")
+            push!(tmpdims, PB.NamedDimension(dname, s, [dc for dc in dcoords]))
+        end
+        dims=tmpdims
+    else
+        error("invalid dims: ", dims)
+    end
+
+    # check supplied coords have same named dimensions of same size
+    tmpcoords = FieldArray[] # narrow type
+    for co in coords
+        for codim in co.dims
+            arraydimidx = findfirst(d->d.name == codim.name, dims)
+            !isnothing(arraydimidx) || error("coord '$(co.name)' dimension '$(codim.name)' not present in array dimensions $([d.name for d in dims])")
+            arraydim = dims[arraydimidx]
+            codim.size == arraydim.size || error("coord '$(co.name)' dimension '$(codim.name)': size $(codim.size) != array dimension size $(arraydim.size)")
+        end
+        push!(tmpcoords, co)
+    end
+
+    return FieldArray(name, values, Tuple(dims), tmpcoords, attributes)
+end
+
+
 function Base.show(io::IO, fa::FieldArray)
     print(io, 
-        "FieldArray(name=\"", fa.name, "\", eltype=", eltype(fa.values),", size=", size(fa.values), ", dims=", fa.dims, ")"
+        "FieldArray(name=\"", fa.name, "\", eltype=", eltype(fa.values), ", dims=", fa.dims, ")"
     )
 end
 
 function Base.show(io::IO, ::MIME"text/plain", fa::FieldArray)
     println(io, "FieldArray(name=\"", fa.name, "\", eltype=", eltype(fa.values),")")
     println(io, "  dims: ", fa.dims)
-    println(io, "  attributes: ", fa.attributes)
+   
+    if !isempty(fa.coords)
+        println(io, "  coords:")
+        for co in fa.coords
+            println(io, "    ", co)
+        end
+    end
+    # println(io, "  attributes: ", fa.attributes)
 end
 
 # basic arithmetic operations
@@ -52,6 +143,84 @@ function default_fieldarray_name(attributes::Dict)
     return name
 end
 
+"""
+    get_dimension(fa::FieldArray, dimname::AbstractString) -> nd::NamedDimension
+
+Get dimension `dimname` (error if not present)
+"""
+function get_dimension(fa::FieldArray, dimname::AbstractString)
+    arraydimidx = findfirst(d->d.name == dimname, fa.dims)
+    !isnothing(arraydimidx) || error("dimension '$dimname' not present")
+    return dims[arraydimidx]
+end
+
+"""
+    get_coord_arrays(fa::FieldArray, dimname::AbstractString) -> coord_arrays::Vector{FieldArray}
+
+Get coordinates for dimension `dimname` (empty Vector if no coordinates defined)
+"""
+function get_coord_arrays(fa::FieldArray, dimname::AbstractString)
+    
+    nd = get_dimension(fa, dimname)
+    coord_names = nd.coords
+    coord_arrays = FieldArray[]
+    for cname in coord_names
+        cnidx = findfirst(c->c.name == cname, fa.coords)
+        !isnothing(cnidx) || error("dimension '$dimname' coordinate '$cname' not present")
+        push!(coord_arrays, fa.coords[cnidx])
+    end
+
+    return coord_arrays
+end
+
+"""
+    update_coordinates(varray::FieldArray, vec_coord_arrays::AbstractVector) -> FieldArray
+
+Replace or add coordinates `vec_coord_arrays` to `varray`.
+
+`new_coord_arrays` is a Vector of Pairs of "dim_name"=>(var1::FieldArray, var2::FieldArray, ...)
+
+Example: to replace a 1D column default pressure coordinate with a z coordinate:
+ 
+    coords=["z"=>(zmid::FieldArray, zlower::FieldArray, atm.zupper::FieldArray)]
+"""
+function update_coordinates(varray::FieldArray, vec_coord_arrays::AbstractVector)
+
+    check_coords_argument(vec_coord_arrays) || 
+        error("argument vec_coords_arrays should be a Vector of Pairs of \"dim_name\"=>(var1::FieldArray, var2::FieldArray, ...), eg: [\"z\"=>(zmid::FieldArray, zlower::FieldArray, atm.zupper::FieldArray)]")
+
+    # generate Vector of NamedDimensions to use as new coordinates
+    named_dimensions = PB.NamedDimension[]
+    for (coord_name, coord_arrays) in vec_coord_arrays
+        fixed_coords = []
+        for coord_array in coord_arrays
+            push!(fixed_coords, PB.FixedCoord(get(coord_array.attributes, :var_name, ""), coord_array.values, coord_array.attributes))
+        end
+        push!(
+            named_dimensions, PB.NamedDimension(
+                coord_name,
+                length(first(fixed_coords).values),
+                fixed_coords,
+            )
+        )
+    end
+
+    # replace coordinates
+    varray_newcoords = FieldArray(varray.name, varray.values, Tuple(named_dimensions), varray.attributes)
+
+    return varray_newcoords
+end
+
+# check 'coords' of form [] or ["z"=>[ ... ], ] or ["z"=>(...),]
+check_coords_argument(coords) =
+    isa(coords, AbstractVector) && (
+        isempty(coords) || (
+            isa(coords, AbstractVector{<:Pair}) &&
+            isa(first(first(coords)), AbstractString) &&
+            isa(last(first(coords)), Union{AbstractVector, Tuple})
+        )
+    )
+
 #############################################################
 # Create from PALEO objects
 #############################################################
@@ -184,50 +353,154 @@ function get_array(
     return varray
 end
 
-# check 'coords' of form [] or ["z"=>[ ... ], ] or ["z"=>(...),]
-check_coords_argument(coords) =
-    isa(coords, AbstractVector) && (
-        isempty(coords) || (
-            isa(coords, AbstractVector{<:Pair}) &&
-            isa(first(first(coords)), AbstractString) &&
-            isa(last(first(coords)), Union{AbstractVector, Tuple})
-        )
-    )
+
+
+
+
+#########################################
+# get_region moved from PALEOboxes.jl
+#########################################
 
 """
-    update_coordinates(varray::FieldArray, vec_coord_arrays::AbstractVector) -> FieldArray
+    get_region(grid::Union{PB.AbstractMesh, Nothing}, values; selectargs...) -> values_subset, (dim_subset::NamedDimension, ...)
 
-Replace or add coordinates `vec_coord_arrays` to `varray`.
+Return the subset of `values` given by `selectargs` (Grid-specific keywords eg cell=, column=, ...)
+and corresponding dimensions (with attached coordinates).
+"""
+function get_region(grid::Union{PB.AbstractMesh, Nothing}, values) end
 
-`new_coord_arrays` is a Vector of Pairs of "coord_name"=>(var1::FieldArray, var2::FieldArray, ...)
+"""
+    get_region(grid::Nothing, values) -> values[]
 
-Example: to replace a 1D column default pressure coordinate with a z coordinate:
- 
-    coords=["z"=>(zmid::FieldArray, zlower::FieldArray, atm.zupper::FieldArray)]
+Fallback for Domain with no grid, assumed 1 cell
 """
-function update_coordinates(varray::FieldArray, vec_coord_arrays::AbstractVector)
+function get_region(grid::Nothing, values)
+    length(values) == 1 ||
+        throw(ArgumentError("grid==Nothing and length(values) != 1"))
+    return values[], ()
+end
+   
+"""
+    get_region(grid::UnstructuredVectorGrid, values; cell) -> 
+        values_subset, (dim_subset::NamedDimension, ...)
 
-    check_coords_argument(vec_coord_arrays) || 
-        error("argument vec_coords_arrays should be a Vector of Pairs of \"coord_name\"=>(var1::FieldArray, var2::FieldArray, ...), eg: [\"z\"=>(zmid::FieldArray, zlower::FieldArray, atm.zupper::FieldArray)]")
+# Keywords for region selection:
+- `cell::Union{Int, Symbol}`: an Int, or a Symbol to look up in `cellnames`
+"""
+function get_region(grid::PB.Grids.UnstructuredVectorGrid, values; cell::Union{Int, Symbol})
+    if cell isa Int
+        idx = cell
+    else
+        idx = get(grid.cellnames, cell, nothing)
+        !isnothing(idx) ||
+            throw(ArgumentError("cell ':$cell' not present in  grid.cellnames=$(grid.cellnames)"))
+    end
 
-    # generate Vector of NamedDimensions to use as new coordinates
-    named_dimensions = PB.NamedDimension[]
-    for (coord_name, coord_arrays) in vec_coord_arrays
-        fixed_coords = []
-        for coord_array in coord_arrays
-            push!(fixed_coords, PB.FixedCoord(get(coord_array.attributes, :var_name, ""), coord_array.values, coord_array.attributes))
-        end
-        push!(
-            named_dimensions, PB.NamedDimension(
-                coord_name,
-                length(first(fixed_coords).values),
-                fixed_coords,
-            )
+    return (
+        values[idx],
+        (),  # no dimensions (ie squeeze out a dimension length 1 for single cell)
+    )
+end
+
+"""
+    get_region(grid::UnstructuredColumnGrid, values; column, [cell=nothing]) -> 
+        values_subset, (dim_subset::NamedDimension, ...)
+
+# Keywords for region selection:
+- `column::Union{Int, Symbol}`: (may be an Int, or a Symbol to look up in `columnames`)
+- `cell::Int`: optional cell index within `column`, highest cell is cell 1
+"""
+function get_region(grid::PB.Grids.UnstructuredColumnGrid, values; column, cell::Union{Nothing, Int}=nothing)
+
+    if column isa Int
+        column in 1:length(grid.Icolumns) ||
+            throw(ArgumentError("column index $column out of range"))
+        colidx = column
+    else
+        colidx = findfirst(isequal(column), grid.columnnames)
+        !isnothing(colidx) || 
+            throw(ArgumentError("columnname '$column' not present in  grid.columnnames=$(grid.columnnames)"))
+    end
+
+    if isnothing(cell)
+        indices = grid.Icolumns[colidx]
+        return (
+            values[indices],
+            (PB.NamedDimension("z", length(indices), PB.get_region(grid.z_coords, indices)), ),
+        )
+    else
+        # squeeze out z dimension
+        idx = grid.Icolumns[colidx][cell]
+        return (
+            values[idx],
+            (),  # no dimensions (ie squeeze out a dimension length 1 for single cell)
         )
     end
+    
+end
 
-    # replace coordinates
-    varray_newcoords = FieldArray(varray.name, varray.values, Tuple(named_dimensions), varray.attributes)
+"""
+    get_region(grid::Union{CartesianLinearGrid{2}, CartesianArrayGrid{2}} , internalvalues; [i=i_idx], [j=j_idx]) ->
+        arrayvalues_subset, (dim_subset::NamedDimension, ...)
 
-    return varray_newcoords
-end
\ No newline at end of file
+# Keywords for region selection:
+- `i::Int`: optional, slice along first dimension
+- `j::Int`: optional, slice along second dimension
+
+`internalvalues` are transformed if needed from internal Field representation as a Vector length `ncells`, to
+an Array (2D if neither i, j arguments present, 1D if i or j present, 0D ie one cell if both present)
+"""
+function get_region(
+    grid::Union{PB.Grids.CartesianLinearGrid{2}, PB.Grids.CartesianArrayGrid{2}}, internalvalues; 
+    i::Union{Integer, Colon}=Colon(), j::Union{Integer, Colon}=Colon()
+)
+    return _get_region(grid, internalvalues, [i, j])
+end
+
+"""
+    get_region(grid::Union{CartesianLinearGrid{3}, CartesianArrayGrid{3}}, internalvalues; [i=i_idx], [j=j_idx]) ->
+        arrayvalues_subset, (dim_subset::NamedDimension, ...)
+
+# Keywords for region selection:
+- `i::Int`: optional, slice along first dimension
+- `j::Int`: optional, slice along second dimension
+- `k::Int`: optional, slice along third dimension
+
+`internalvalues` are transformed if needed from internal Field representation as a Vector length `ncells`, to
+an Array (3D if neither i, j, k arguments present, 2D if one of i, j or k present, 1D if two present,
+0D ie one cell if i, j, k all specified).
+"""
+function get_region(
+    grid::Union{PB.Grids.CartesianLinearGrid{3}, PB.Grids.CartesianArrayGrid{3}}, internalvalues;
+    i::Union{Integer, Colon}=Colon(), j::Union{Integer, Colon}=Colon(), k::Union{Integer, Colon}=Colon()
+)
+    return _get_region(grid, internalvalues, [i, j, k])
+end
+
+function _get_region(
+    grid::Union{PB.Grids.CartesianLinearGrid, PB.Grids.CartesianArrayGrid}, internalvalues, indices
+)
+    if !isempty(grid.coords) && !isempty(grid.coords_edges)
+        dims = [
+            PB.NamedDimension(grid.dimnames[idx], grid.coords[idx], grid.coords_edges[idx])
+            for (idx, ind) in enumerate(indices) if isa(ind, Colon)
+        ]
+    elseif !isempty(grid.coords)
+        dims = [
+            PB.NamedDimension(grid.dimnames[idx], grid.coords[idx])
+            for (idx, ind) in enumerate(indices) if isa(ind, Colon)
+        ]
+    else
+        dims = [
+            PB.NamedDimension(grid.dimnames[idx])
+            for (idx, ind) in enumerate(indices) if isa(ind, Colon)
+        ]
+    end
+
+    values = internal_to_cartesian(grid, internalvalues)
+    if !all(isequal(Colon()), indices)
+        values = values[indices...]
+    end
+
+    return values, Tuple(dims)    
+end
diff --git a/src/FieldRecord.jl b/src/FieldRecord.jl
index 16d5e77..bdb6d84 100644
--- a/src/FieldRecord.jl
+++ b/src/FieldRecord.jl
@@ -1,3 +1,20 @@
+mutable struct RecordCoord
+    name::String
+    values::Vector{Float64}
+    attributes::Dict{Symbol, Any}
+end
+
+
+
+function get_region(fc::RecordCoord, indices::AbstractVector)
+    return FixedCoord(fc.name, fc.values[indices], fc.attributes)
+end
+
+function get_region(fcv::Vector{RecordCoord}, indices::AbstractVector)
+    return [RecordCoord(fc.name, fc.values[indices], fc.attributes) for fc in fcv]
+end
+
+
 """
     FieldRecord{D <: AbstractData, S <: AbstractSpace, V, N, M, R}
     FieldRecord(
@@ -17,15 +34,17 @@ Fields with single `values` (`field_single_element` true) are stored as a Vector
 struct FieldRecord{D <: PB.AbstractData, S <: PB.AbstractSpace, V, N, M, R}
     records::Vector{R}
     data_dims::NTuple{N, PB.NamedDimension}
+    record_dim::PB.NamedDimension
     mesh::M
+    coords::Vector{FieldArray}
     attributes::Dict{Symbol, Any}
-    coords_record::Vector{PB.FixedCoord} # coordinates attached to record dimension
+    # coords_record::Vector{RecordCoord} # coordinates attached to record dimension
 end
 
 function Base.show(io::IO, fr::FieldRecord)
     print(io, 
-        "FieldRecord(eltype=", eltype(fr),", length=", length(fr), 
-        ", attributes=", fr.attributes, ", coords_record=", fr.coords_record, ")"
+        "FieldRecord(eltype=", eltype(fr),", length=", length(fr), ", record_dim=", fr.record_dim, 
+        ", attributes=", fr.attributes, ")"
     )
 end
 
@@ -223,7 +242,7 @@ function _get_array(
             # find ridx corresponding to a coordinate
             for cr in fr.coords_record
                 if String(k) == cr.name
-                    ridx, cvalue = PB.find_indices(cr.values, v)
+                    ridx, cvalue = find_indices(cr.values, v)
                     selectargs_records=NamedTuple((k=>cvalue,))
                 end
             end            
@@ -231,7 +250,8 @@ function _get_array(
             selectargs_region = merge(selectargs_region, NamedTuple((k=>v,)))
         end
     end
-    records_dim = PB.NamedDimension("records", length(ridx), PB.get_region(fr.coords_record, ridx))
+    records_coord = 
+    records_dim = PB.NamedDimension("records", length(ridx), get_region(fr.coords_record, ridx))
 
     # add attributes for selection used
     attributes = copy(fr.attributes)
@@ -312,3 +332,28 @@ function _get_array(
 
 end
 
+"find indices of coord from first before range[1] to first after range[2]"
+function find_indices(coord::AbstractVector, range)
+    length(range) == 2 ||
+        throw(ArgumentError("find_indices: length(range) != 2  $range"))
+
+    idxstart = findlast(t -> t<=range[1], coord)
+    isnothing(idxstart) && (idxstart = 1)
+
+    idxend = findfirst(t -> t>=range[2], coord)
+    isnothing(idxend) && (idxend = length(coord))
+
+    return idxstart:idxend, (coord[idxstart], coord[idxend])
+end
+
+"find indices of coord nearest val"
+function find_indices(coord::AbstractVector, val::Real)
+    idx = 1
+    for i in 1:length(coord)
+        if abs(coord[i] - val) < abs(coord[idx] - val)
+            idx = i
+        end
+    end
+
+    return [idx], coord[idx]
+end
\ No newline at end of file
diff --git a/src/OutputWriters.jl b/src/OutputWriters.jl
index 9302709..50fcb04 100644
--- a/src/OutputWriters.jl
+++ b/src/OutputWriters.jl
@@ -440,7 +440,7 @@ function PB.get_field(odom::OutputMemoryDomain, varname)
     fr = PALEOmodel.wrap_fieldrecord(
         vdata, field_data, Tuple(data_dims), missing, space, grid, attributes,
         coords_record=[
-            PB.FixedCoord(
+            PALEOmodel.RecordCoord(
                 String(odom.coords_record),
                 # df[!, odom.coords_record],
                 use_all_records ? df[!, odom.coords_record] : df[1:odom._nrecs, odom.coords_record],
diff --git a/src/PlotRecipes.jl b/src/PlotRecipes.jl
index e3cc46c..b2958e7 100644
--- a/src/PlotRecipes.jl
+++ b/src/PlotRecipes.jl
@@ -364,6 +364,12 @@ RecipesBase.@recipe function f(
         end
     end
 
+    # guess coordinate edges from midpoints, assuming uniform spacing
+    function guess_coords_edges(x_midpoints)
+        first_x = x_midpoints[1] - 0.5*(x_midpoints[2] - x_midpoints[1])
+        last_x = x_midpoints[end] + 0.5*(x_midpoints[end] - x_midpoints[end-1])
+        return [first_x; 0.5.*(x_midpoints[1:end-1] .+ x_midpoints[2:end]); last_x]
+    end
 
     values = fa.values
     name = fa.name
@@ -386,20 +392,20 @@ RecipesBase.@recipe function f(
         else
             label --> labelprefix*name
         end
-        f_label = PB.append_units(get_attribute(fa, :var_name, ""), fa.attributes)
+        f_label = append_units(get_attribute(fa, :var_name, ""), fa.attributes)
         coords_vec = fa.dims[1].coords
 
         if length(coords_vec) == 1 || length(coords_vec) > 3
             co = first(coords_vec)            
             return swap_coords_values(
                 co.values, values,
-                PB.append_units(co.name, co.attributes), f_label,
+                append_units(co.name, co.attributes), f_label,
             )
             
         elseif length(coords_vec) in (2, 3)
             co_lower = coords_vec[end-1]
             co_upper = coords_vec[end]
-            c_label = PB.append_units(co_lower.name*", "*co_upper.name, co_lower.attributes)           
+            c_label = append_units(co_lower.name*", "*co_upper.name, co_lower.attributes)           
             return swap_coords_values(
                 create_stepped(co_lower.values, co_upper.values, values)...,
                 c_label, f_label,
@@ -412,19 +418,19 @@ RecipesBase.@recipe function f(
         end
 
     elseif length(fa.dims) == 2
-        title --> PB.append_units(fa.name, fa.attributes)        
+        title --> append_units(fa.name, fa.attributes)        
        
-        i_values, i_label = PB.build_coords_edges(fa.dims[1])
-        j_values, j_label = PB.build_coords_edges(fa.dims[2])
+        i_values, i_label = build_coords_edges(fa, fa.dims[1].name)
+        j_values, j_label = build_coords_edges(fa, fa.dims[2].name)
 
         # workaround for Plots.jl heatmap: needs x, y to either both be midpoints, or both be edges
         if length(i_values) == size(values, 1) && length(j_values) == size(values, 2)+1
             @warn "$(fa.name) guessing '$i_label' coordinate edges from midpoints assuming uniform spacing"
-            i_values = PB.guess_coords_edges(i_values)
+            i_values = guess_coords_edges(i_values)
         end
         if length(i_values) == size(values, 1)+1 && length(j_values) == size(values, 2)
             @warn "$(fa.name) guessing '$j_label' coordinate edges from midpoints assuming uniform spacing"
-            j_values = PB.guess_coords_edges(j_values)
+            j_values = guess_coords_edges(j_values)
         end
 
         x_values, y_values, f_values = swap_coords_xy(
@@ -527,3 +533,79 @@ function broadcast_dict(dv::Vector{<:Dict})
         return dvout
     end
 end
+
+"""
+    build_coords_edges(nd::NamedDimension) -> Vector{Float64}
+
+Call [`build_coords_edges`](@ref)(nd.coords), or fallback to just returning indices
+if no coords present.
+"""
+function build_coords_edges(nd::PB.NamedDimension)
+    if !isempty(nd.coords)
+        return build_coords_edges(nd.coords)
+    else
+        @warn "no coords for NamedDimension $(nd.name), returning indices"
+        return collect(1:nd.size), nd.name*" (indices)"
+    end
+end
+
+"""
+    build_coords_edges(coords_vec::Vector{FixedCoord}) -> Vector{Float64}
+
+Build a vector of coordinate edges (length `n+1``) from `coords_vec`, assuming the PALEO
+convention that `coords_vec` contains three elements with 
+cell midpoints, lower edges, upper edges each of length `n`, in that order.
+
+Falls back to just returning the first entry in `coords_vec` for other cases.
+"""
+function build_coords_edges(fa::FieldArray, dimname::AbstractString)
+    
+    coords_vec = get_coord_arrays(fa, dimname)
+
+    if isempty(coords_vec)
+        nd = get_dimension(fa, dimname)
+        @warn "no coords for dimension $(nd), returning indices"
+        return collect(1:nd.size), nd.name*" (indices)"
+    elseif length(coords_vec) == 1 || length(coords_vec) > 3
+        # 1 coordinate or something we don't understand - take first
+        co = first(coords_vec)
+        co_values = co.values
+        co_label = append_units(co.name, co.attributes)
+    elseif length(coords_vec) in (2, 3)
+        # 2 coordinates assume lower, upper edges
+        # 3 coordinates assume mid, lower, upper
+        co_lower = coords_vec[end-1]
+        co_upper = coords_vec[end]
+        co_label = append_units(co_lower.name*", "*co_upper.name, co_lower.attributes)
+        first(co_lower.values) < first(co_upper.values) ||
+            @warn "build_coords_edges: $co_label co_lower is > co_upper - check model grid"            
+        if co_lower.values[end] > co_lower.values[1] # ascending order
+            co_lower.values[2:end] == co_upper.values[1:end-1] ||
+                @warn "build_coords_edges: $co_label lower and upper edges don't match"
+            co_values = [co_lower.values; co_upper.values[end]]
+        else # descending order
+            co_lower.values[1:end-1] == co_upper.values[2:end] ||
+                @warn "build_coords_edges: $co_label lower and upper edges don't match"
+            co_values = [co_upper.values[1]; co_lower.values]
+        end
+
+    end
+
+    return co_values, co_label
+end
+
+"""
+    append_units(name::AbstractString, attributes) -> "name (units)"
+
+Utility function to append variable units string to a variable name for display.
+"""
+function append_units(name::AbstractString, attributes::Dict{Symbol, Any})
+    units = get(attributes, :units, "")
+    if isempty(units)
+        return name
+    else
+        return name*" ($units)"
+    end
+end
+
+append_units(name::AbstractString, attributes::Nothing) = name