Bug in NetCDF output for curvilinear immersed grids (#3666)

Co-authored-by: Ali Ramadhan <[email protected]>

src/Grids/latitude_longitude_grid.jl

@@ -182,7 +182,7 @@ function LatitudeLongitudeGrid(architecture::AbstractArchitecture = CPU(),
                                precompute_metrics = true,
                                halo = nothing)
 
-    if architecture == GPU() && !has_cuda() 
+    if architecture == GPU() && !has_cuda()
         throw(ArgumentError("Cannot create a GPU grid. No CUDA-enabled GPU was detected!"))
     end
 
@@ -193,7 +193,7 @@ function LatitudeLongitudeGrid(architecture::AbstractArchitecture = CPU(),
     Hλ, Hφ, Hz = halo
 
     # Calculate all direction (which might be stretched)
-    # A direction is regular if the domain passed is a Tuple{<:Real, <:Real}, 
+    # A direction is regular if the domain passed is a Tuple{<:Real, <:Real},
     # it is stretched if being passed is a function or vector (as for the VerticallyStretchedRectilinearGrid)
     TX, TY, TZ = topology
 
@@ -281,7 +281,7 @@ function validate_lat_lon_grid_args(topology, size, halo, FT, latitude, longitud
     φ₂ <= 90  || throw(ArgumentError("The northern latitude cannot be less than -90 degrees."))
     φ₁ <= φ₂  || throw(ArgumentError("Latitudes must increase south to north."))
 
-    if TX == Flat || TY == Flat 
+    if TX == Flat || TY == Flat
         precompute_metrics = false
     end
 
@@ -313,7 +313,7 @@ function Base.show(io::IO, grid::LatitudeLongitudeGrid, withsummary=true)
     Ωφ = domain(TY(), size(grid, 2), grid.φᵃᶠᵃ)
     Ωz = domain(TZ(), size(grid, 3), grid.zᵃᵃᶠ)
 
-    x_summary = domain_summary(TX(), "λ", Ωλ) 
+    x_summary = domain_summary(TX(), "λ", Ωλ)
     y_summary = domain_summary(TY(), "φ", Ωφ)
     z_summary = domain_summary(TZ(), "z", Ωz)
 
@@ -457,20 +457,20 @@ end
 @inline Azᶜᶜᵃ(i, j, k, grid::XRegularLLG) = @inbounds grid.radius^2 * deg2rad(grid.Δλᶜᵃᵃ) * (hack_sind(grid.φᵃᶠᵃ[j+1]) - hack_sind(grid.φᵃᶠᵃ[j]))
 
 #####
-##### Utilities to precompute metrics 
+##### Utilities to precompute metrics
 #####
 
-@inline metrics_precomputed(::LatitudeLongitudeGrid{<:Any, <:Any, <:Any, <:Any, Nothing}) = false 
+@inline metrics_precomputed(::LatitudeLongitudeGrid{<:Any, <:Any, <:Any, <:Any, Nothing}) = false
 @inline metrics_precomputed(::LatitudeLongitudeGrid) = true
 
 #####
 ##### Kernels that precompute the z- and x-metric
 #####
 
-@inline metric_worksize(grid::LatitudeLongitudeGrid)  = (length(grid.Δλᶜᵃᵃ), length(grid.φᵃᶠᵃ) - 2) 
-@inline metric_workgroup(grid::LatitudeLongitudeGrid) = (16, 16) 
+@inline metric_worksize(grid::LatitudeLongitudeGrid)  = (length(grid.Δλᶜᵃᵃ), length(grid.φᵃᶠᵃ) - 2)
+@inline metric_workgroup(grid::LatitudeLongitudeGrid) = (16, 16)
 
-@inline metric_worksize(grid::XRegularLLG)  = length(grid.φᵃᶠᵃ) - 2 
+@inline metric_worksize(grid::XRegularLLG)  = length(grid.φᵃᶠᵃ) - 2
 @inline metric_workgroup(grid::XRegularLLG) = 16
 
 @kernel function compute_Δx_Az!(grid::LatitudeLongitudeGrid, Δxᶠᶜ, Δxᶜᶠ, Δxᶠᶠ, Δxᶜᶜ, Azᶠᶜ, Azᶜᶠ, Azᶠᶠ, Azᶜᶜ)
@@ -534,7 +534,7 @@ function allocate_metrics(grid::LatitudeLongitudeGrid)
     FT = eltype(grid)
 
     arch = grid.architecture
-    
+
     if grid isa XRegularLLG
         offsets     = grid.φᵃᶜᵃ.offsets[1]
         metric_size = length(grid.φᵃᶜᵃ)
@@ -561,7 +561,7 @@ function allocate_metrics(grid::LatitudeLongitudeGrid)
         Δyᶠᶜ    = OffsetArray(on_architecture(arch, parentC), grid.Δφᵃᶜᵃ.offsets[1])
         Δyᶜᶠ    = OffsetArray(on_architecture(arch, parentF), grid.Δφᵃᶜᵃ.offsets[1])
     end
-    
+
     return Δxᶠᶜ, Δxᶜᶠ, Δxᶠᶠ, Δxᶜᶜ, Δyᶠᶜ, Δyᶜᶠ, Azᶠᶜ, Azᶜᶠ, Azᶠᶠ, Azᶜᶜ
 end
 
@@ -616,7 +616,7 @@ function nodes(grid::LLG, ℓx, ℓy, ℓz; reshape=false, with_halos=false)
         # consideration...
         #
         # See also `nodes` for `RectilinearGrid`.
-        
+
         Nλ = isnothing(λ) ? 1 : length(λ)
         Nφ = isnothing(φ) ? 1 : length(φ)
         Nz = isnothing(z) ? 1 : length(z)
@@ -655,6 +655,15 @@ end
 @inline xnodes(grid::LLG, ℓx, ℓy, ℓz; with_halos=false) = xnodes(grid, ℓx, ℓy; with_halos)
 @inline ynodes(grid::LLG, ℓx, ℓy, ℓz; with_halos=false) = ynodes(grid, ℓy; with_halos)
 
+# Generalized coordinates
+@inline ξnodes(grid::LLG, ℓx; kwargs...) = λnodes(grid, ℓx; kwargs...)
+@inline ηnodes(grid::LLG, ℓy; kwargs...) = φnodes(grid, ℓy; kwargs...)
+@inline rnodes(grid::LLG, ℓz; kwargs...) = znodes(grid, ℓz; kwargs...)
+
+@inline ξnodes(grid::LLG, ℓx, ℓy, ℓz; kwargs...) = λnodes(grid, ℓx; kwargs...)
+@inline ηnodes(grid::LLG, ℓx, ℓy, ℓz; kwargs...) = φnodes(grid, ℓy; kwargs...)
+@inline rnodes(grid::LLG, ℓx, ℓy, ℓz; kwargs...) = znodes(grid, ℓz; kwargs...)
+
 #####
 ##### Grid spacings in x, y, z (in meters)
 #####
@@ -670,11 +679,11 @@ end
 @inline xspacings(grid::LLG, ℓx::F, ℓy::C; with_halos=false) = _property(grid.Δxᶠᶜᵃ, ℓx, ℓy,
                                                                          topology(grid, 1), topology(grid, 2),
                                                                          size(grid, 1), size(grid, 2), with_halos)
- 
+
 @inline xspacings(grid::LLG, ℓx::F, ℓy::F; with_halos=false) = _property(grid.Δxᶠᶠᵃ, ℓx, ℓy,
                                                                          topology(grid, 1), topology(grid, 2),
                                                                          size(grid, 1), size(grid, 2), with_halos)
- 
+
 @inline xspacings(grid::HRegularLLG, ℓx::C, ℓy::C; with_halos=false) = _property(grid.Δxᶜᶜᵃ, ℓy, topology(grid, 2),
                                                                                  size(grid, 2), with_halos)
 
@@ -737,4 +746,3 @@ end
 
 @inline λspacing(i, j, k, grid::LLG, ℓx, ℓy, ℓz) = λspacing(i, grid, ℓx)
 @inline φspacing(i, j, k, grid::LLG, ℓx, ℓy, ℓz) = φspacing(j, grid, ℓy)
-

src/Grids/rectilinear_grid.jl

@@ -38,7 +38,7 @@ struct RectilinearGrid{FT, TX, TY, TZ, FX, FY, FZ, VX, VY, VZ, Arch} <: Abstract
                                                                  FX, VX, FY,
                                                                  VY, FZ, VZ} =
         new{FT, TX, TY, TZ, FX, FY, FZ, VX, VY, VZ, Arch}(arch, Nx, Ny, Nz,
-                                                          Hx, Hy, Hz, Lx, Ly, Lz, 
+                                                          Hx, Hy, Hz, Lx, Ly, Lz,
                                                           Δxᶠᵃᵃ, Δxᶜᵃᵃ, xᶠᵃᵃ, xᶜᵃᵃ,
                                                           Δyᵃᶠᵃ, Δyᵃᶜᵃ, yᵃᶠᵃ, yᵃᶜᵃ,
                                                           Δzᵃᵃᶠ, Δzᵃᵃᶜ, zᵃᵃᶠ, zᵃᵃᶜ)
@@ -261,7 +261,7 @@ function RectilinearGrid(architecture::AbstractArchitecture = CPU(),
                          extent = nothing,
                          topology = (Periodic, Periodic, Bounded))
 
-    if architecture == GPU() && !has_cuda() 
+    if architecture == GPU() && !has_cuda()
         throw(ArgumentError("Cannot create a GPU grid. No CUDA-enabled GPU was detected!"))
     end
 
@@ -388,7 +388,7 @@ function constructor_arguments(grid::RectilinearGrid)
                   :x => cpu_face_constructor_x(grid),
                   :y => cpu_face_constructor_y(grid),
                   :z => cpu_face_constructor_z(grid),
-                  :topology => topo) 
+                  :topology => topo)
 
     return args, kwargs
 end
@@ -479,7 +479,7 @@ function nodes(grid::RectilinearGrid, ℓx, ℓy, ℓz; reshape=false, with_halo
         # might be to omit the `nothing` nodes in the `reshape`. In other words,
         # if `TX === Flat`, then we should return `(x, z)`. This is for future
         # consideration...
-        # 
+        #
         # See also `nodes` for `LatitudeLongitudeGrid`.
 
         Nx = isnothing(x) ? 1 : length(x)
@@ -509,6 +509,15 @@ const C = Center
 @inline ynodes(grid::RG, ℓx, ℓy, ℓz; with_halos=false) = ynodes(grid, ℓy; with_halos)
 @inline znodes(grid::RG, ℓx, ℓy, ℓz; with_halos=false) = znodes(grid, ℓz; with_halos)
 
+# Generalized coordinates
+@inline ξnodes(grid::RG, ℓx; kwargs...) = xnodes(grid, ℓx; kwargs...)
+@inline ηnodes(grid::RG, ℓy; kwargs...) = ynodes(grid, ℓy; kwargs...)
+@inline rnodes(grid::RG, ℓz; kwargs...) = znodes(grid, ℓz; kwargs...)
+
+@inline ξnodes(grid::RG, ℓx, ℓy, ℓz; kwargs...) = xnodes(grid, ℓx; kwargs...)
+@inline ηnodes(grid::RG, ℓx, ℓy, ℓz; kwargs...) = ynodes(grid, ℓy; kwargs...)
+@inline rnodes(grid::RG, ℓx, ℓy, ℓz; kwargs...) = znodes(grid, ℓz; kwargs...)
+
 #####
 ##### Grid spacings
 #####
@@ -525,4 +534,3 @@ const C = Center
 @inline zspacings(grid::RG, ℓx, ℓy, ℓz; kwargs...) = zspacings(grid, ℓz; kwargs...)
 
 @inline isrectilinear(::RG) = true
-

src/ImmersedBoundaries/ImmersedBoundaries.jl

@@ -1,7 +1,7 @@
 module ImmersedBoundaries
 
 export ImmersedBoundaryGrid, GridFittedBoundary, GridFittedBottom, PartialCellBottom, ImmersedBoundaryCondition
-       
+
 using Adapt
 
 using Oceananigans.Grids
@@ -52,6 +52,7 @@ import Oceananigans.Grids: architecture, on_architecture, with_halo, inflate_hal
                            ξnode, ηnode, rnode,
                            ξname, ηname, rname, node_names,
                            xnodes, ynodes, znodes, λnodes, φnodes, nodes,
+                           ξnodes, ηnodes, rnodes,
                            inactive_cell
 
 import Oceananigans.Coriolis: φᶠᶠᵃ
@@ -145,15 +146,15 @@ const IBG = ImmersedBoundaryGrid
 @inline z_domain(ibg::IBG) = z_domain(ibg.underlying_grid)
 
 Adapt.adapt_structure(to, ibg::IBG{FT, TX, TY, TZ}) where {FT, TX, TY, TZ} =
-    ImmersedBoundaryGrid{TX, TY, TZ}(adapt(to, ibg.underlying_grid), 
-                                     adapt(to, ibg.immersed_boundary), 
-                                     nothing, 
+    ImmersedBoundaryGrid{TX, TY, TZ}(adapt(to, ibg.underlying_grid),
+                                     adapt(to, ibg.immersed_boundary),
+                                     nothing,
                                      nothing)
 
 with_halo(halo, ibg::ImmersedBoundaryGrid) =
     ImmersedBoundaryGrid(with_halo(halo, ibg.underlying_grid), ibg.immersed_boundary)
 
-# ImmersedBoundaryGrids require an extra halo point to check the "inactivity" of a `Face` node at N + H 
+# ImmersedBoundaryGrids require an extra halo point to check the "inactivity" of a `Face` node at N + H
 # (which requires checking `Center` nodes at N + H and N + H + 1)
 inflate_halo_size_one_dimension(req_H, old_H, _, ::IBG)            = max(req_H + 1, old_H)
 inflate_halo_size_one_dimension(req_H, old_H, ::Type{Flat}, ::IBG) = 0
@@ -198,7 +199,7 @@ is not part of the prognostic state.
 
 Return `true` if the tracer cell at `i, j, k` either (i) lies outside the `Bounded` domain
 or (ii) lies within the immersed region of `ImmersedBoundaryGrid`.
-        
+
 Example
 =======
 
@@ -215,7 +216,7 @@ Consider the configuration
  × === ∘ === ×     ∘     ×
  | ========= |           |
 
-i-1          i          
+i-1          i
  f           f           f
 ```
 
@@ -248,11 +249,12 @@ const f = Face()
 @inline Base.zero(ibg::IBG) = zero(ibg.underlying_grid)
 @inline φᶠᶠᵃ(i, j, k, ibg::IBG) = φᶠᶠᵃ(i, j, k, ibg.underlying_grid)
 
-@inline λnode(i, j, k, ibg::IBG, ℓ...) = λnode(i, j, k, ibg.underlying_grid, ℓ...)
-@inline φnode(i, j, k, ibg::IBG, ℓ...) = φnode(i, j, k, ibg.underlying_grid, ℓ...)
-@inline xnode(i, j, k, ibg::IBG, ℓ...) = xnode(i, j, k, ibg.underlying_grid, ℓ...)
-@inline ynode(i, j, k, ibg::IBG, ℓ...) = ynode(i, j, k, ibg.underlying_grid, ℓ...)
-@inline znode(i, j, k, ibg::IBG, ℓ...) = znode(i, j, k, ibg.underlying_grid, ℓ...)
+@inline xnode(i, j, k, ibg::IBG, ℓx, ℓy, ℓz) = xnode(i, j, k, ibg.underlying_grid, ℓx, ℓy, ℓz)
+@inline ynode(i, j, k, ibg::IBG, ℓx, ℓy, ℓz) = ynode(i, j, k, ibg.underlying_grid, ℓx, ℓy, ℓz)
+@inline znode(i, j, k, ibg::IBG, ℓx, ℓy, ℓz) = znode(i, j, k, ibg.underlying_grid, ℓx, ℓy, ℓz)
+
+@inline λnode(i, j, k, ibg::IBG, ℓx, ℓy, ℓz) = λnode(i, j, k, ibg.underlying_grid, ℓx, ℓy, ℓz)
+@inline φnode(i, j, k, ibg::IBG, ℓx, ℓy, ℓz) = φnode(i, j, k, ibg.underlying_grid, ℓx, ℓy, ℓz)
 
 @inline ξnode(i, j, k, ibg::IBG, ℓx, ℓy, ℓz) = ξnode(i, j, k, ibg.underlying_grid, ℓx, ℓy, ℓz)
 @inline ηnode(i, j, k, ibg::IBG, ℓx, ℓy, ℓz) = ηnode(i, j, k, ibg.underlying_grid, ℓx, ℓy, ℓz)
@@ -263,18 +265,28 @@ const f = Face()
 nodes(ibg::IBG, ℓx, ℓy, ℓz; kwargs...) = nodes(ibg.underlying_grid, ℓx, ℓy, ℓz; kwargs...)
 nodes(ibg::IBG, (ℓx, ℓy, ℓz); kwargs...) = nodes(ibg, ℓx, ℓy, ℓz; kwargs...)
 
-λnodes(ibg::IBG, loc; kwargs...) = λnodes(ibg.underlying_grid, loc; kwargs...)
-φnodes(ibg::IBG, loc; kwargs...) = φnodes(ibg.underlying_grid, loc; kwargs...)
 xnodes(ibg::IBG, loc; kwargs...) = xnodes(ibg.underlying_grid, loc; kwargs...)
 ynodes(ibg::IBG, loc; kwargs...) = ynodes(ibg.underlying_grid, loc; kwargs...)
 znodes(ibg::IBG, loc; kwargs...) = znodes(ibg.underlying_grid, loc; kwargs...)
 
-λnodes(ibg::IBG, ℓx, ℓy, ℓz; kwargs...) = λnodes(ibg.underlying_grid, ℓx, ℓy, ℓz; kwargs...)
-φnodes(ibg::IBG, ℓx, ℓy, ℓz; kwargs...) = φnodes(ibg.underlying_grid, ℓx, ℓy, ℓz; kwargs...)
+λnodes(ibg::IBG, loc; kwargs...) = λnodes(ibg.underlying_grid, loc; kwargs...)
+φnodes(ibg::IBG, loc; kwargs...) = φnodes(ibg.underlying_grid, loc; kwargs...)
+
+ξnodes(ibg::IBG, loc; kwargs...) = ξnodes(ibg.underlying_grid, loc; kwargs...)
+ηnodes(ibg::IBG, loc; kwargs...) = ηnodes(ibg.underlying_grid, loc; kwargs...)
+rnodes(ibg::IBG, loc; kwargs...) = rnodes(ibg.underlying_grid, loc; kwargs...)
+
 xnodes(ibg::IBG, ℓx, ℓy, ℓz; kwargs...) = xnodes(ibg.underlying_grid, ℓx, ℓy, ℓz; kwargs...)
 ynodes(ibg::IBG, ℓx, ℓy, ℓz; kwargs...) = ynodes(ibg.underlying_grid, ℓx, ℓy, ℓz; kwargs...)
 znodes(ibg::IBG, ℓx, ℓy, ℓz; kwargs...) = znodes(ibg.underlying_grid, ℓx, ℓy, ℓz; kwargs...)
 
+λnodes(ibg::IBG, ℓx, ℓy, ℓz; kwargs...) = λnodes(ibg.underlying_grid, ℓx, ℓy, ℓz; kwargs...)
+φnodes(ibg::IBG, ℓx, ℓy, ℓz; kwargs...) = φnodes(ibg.underlying_grid, ℓx, ℓy, ℓz; kwargs...)
+
+ξnodes(ibg::IBG, ℓx, ℓy, ℓz; kwargs...) = ξnodes(ibg.underlying_grid, ℓx, ℓy, ℓz; kwargs...)
+ηnodes(ibg::IBG, ℓx, ℓy, ℓz; kwargs...) = ηnodes(ibg.underlying_grid, ℓx, ℓy, ℓz; kwargs...)
+rnodes(ibg::IBG, ℓx, ℓy, ℓz; kwargs...) = rnodes(ibg.underlying_grid, ℓx, ℓy, ℓz; kwargs...)
+
 @inline cpu_face_constructor_x(ibg::IBG) = cpu_face_constructor_x(ibg.underlying_grid)
 @inline cpu_face_constructor_y(ibg::IBG) = cpu_face_constructor_y(ibg.underlying_grid)
 @inline cpu_face_constructor_z(ibg::IBG) = cpu_face_constructor_z(ibg.underlying_grid)
@@ -292,9 +304,9 @@ end
 
 isrectilinear(ibg::IBG) = isrectilinear(ibg.underlying_grid)
 
-@inline fractional_x_index(x, locs, grid::ImmersedBoundaryGrid) = fractional_x_index(x, locs, grid.underlying_grid) 
-@inline fractional_y_index(x, locs, grid::ImmersedBoundaryGrid) = fractional_y_index(x, locs, grid.underlying_grid) 
-@inline fractional_z_index(x, locs, grid::ImmersedBoundaryGrid) = fractional_z_index(x, locs, grid.underlying_grid) 
+@inline fractional_x_index(x, locs, grid::ImmersedBoundaryGrid) = fractional_x_index(x, locs, grid.underlying_grid)
+@inline fractional_y_index(x, locs, grid::ImmersedBoundaryGrid) = fractional_y_index(x, locs, grid.underlying_grid)
+@inline fractional_z_index(x, locs, grid::ImmersedBoundaryGrid) = fractional_z_index(x, locs, grid.underlying_grid)
 
 #####
 ##### Diffusivities (for VerticallyImplicit)
@@ -331,4 +343,3 @@ include("immersed_reductions.jl")
 include("distributed_immersed_boundaries.jl")
 
 end # module
-

src/OutputWriters/netcdf_output_writer.jl

@@ -4,7 +4,7 @@ using Dates: AbstractTime, now
 
 using Oceananigans.Fields
 
-using Oceananigans.Grids: AbstractCurvilinearGrid, RectilinearGrid, topology, halo_size, parent_index_range
+using Oceananigans.Grids: AbstractCurvilinearGrid, RectilinearGrid, topology, halo_size, parent_index_range, ξnodes, ηnodes, rnodes
 using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid
 using Oceananigans.Utils: versioninfo_with_gpu, oceananigans_versioninfo, prettykeys
 using Oceananigans.TimeSteppers: float_or_date_time
@@ -52,46 +52,19 @@ netcdf_spatial_dimensions(::AbstractField{LX, LY, LZ}) where {LX, LY, LZ} =
 function native_dimensions_for_netcdf_output(grid, indices, TX, TY, TZ, Hx, Hy, Hz)
     with_halos = true
 
-    xC = xnodes(grid, c; with_halos)
-    xF = xnodes(grid, f; with_halos)
-    yC = ynodes(grid, c; with_halos)
-    yF = ynodes(grid, f; with_halos)
-    zC = znodes(grid, c; with_halos)
-    zF = znodes(grid, f; with_halos)
+    xC = ξnodes(grid, c; with_halos)
+    xF = ξnodes(grid, f; with_halos)
+    yC = ηnodes(grid, c; with_halos)
+    yF = ηnodes(grid, f; with_halos)
+    zC = rnodes(grid, c; with_halos)
+    zF = rnodes(grid, f; with_halos)
 
-    xC = isnothing(xC) ? [0.0] : parent(xC)  
-    xF = isnothing(xF) ? [0.0] : parent(xF) 
-    yC = isnothing(yC) ? [0.0] : parent(yC) 
-    yF = isnothing(yF) ? [0.0] : parent(yF) 
-    zC = isnothing(zC) ? [0.0] : parent(zC) 
-    zF = isnothing(zF) ? [0.0] : parent(zF) 
-
-    dims = Dict("xC" => xC[parent_index_range(indices["xC"][1], c, TX(), Hx)],
-                "xF" => xF[parent_index_range(indices["xF"][1], f, TX(), Hx)],
-                "yC" => yC[parent_index_range(indices["yC"][2], c, TY(), Hy)],
-                "yF" => yF[parent_index_range(indices["yF"][2], f, TY(), Hy)],
-                "zC" => zC[parent_index_range(indices["zC"][3], c, TZ(), Hz)],
-                "zF" => zF[parent_index_range(indices["zF"][3], f, TZ(), Hz)])
-
-    return dims
-end
-
-function native_dimensions_for_netcdf_output(grid::AbstractCurvilinearGrid, indices, TX, TY, TZ, Hx, Hy, Hz)
-    with_halos = true
-
-    xC = λnodes(grid, c; with_halos)
-    xF = λnodes(grid, f; with_halos)
-    yC = φnodes(grid, c; with_halos)
-    yF = φnodes(grid, f; with_halos)
-    zC = znodes(grid, c; with_halos)
-    zF = znodes(grid, f; with_halos)
-
-    xC = isnothing(xC) ? [0.0] : parent(xC)  
-    xF = isnothing(xF) ? [0.0] : parent(xF) 
-    yC = isnothing(yC) ? [0.0] : parent(yC) 
-    yF = isnothing(yF) ? [0.0] : parent(yF) 
-    zC = isnothing(zC) ? [0.0] : parent(zC) 
-    zF = isnothing(zF) ? [0.0] : parent(zF) 
+    xC = isnothing(xC) ? [0.0] : parent(xC)
+    xF = isnothing(xF) ? [0.0] : parent(xF)
+    yC = isnothing(yC) ? [0.0] : parent(yC)
+    yF = isnothing(yF) ? [0.0] : parent(yF)
+    zC = isnothing(zC) ? [0.0] : parent(zC)
+    zF = isnothing(zF) ? [0.0] : parent(zF)
 
     dims = Dict("xC" => xC[parent_index_range(indices["xC"][1], c, TX(), Hx)],
                 "xF" => xF[parent_index_range(indices["xF"][1], f, TX(), Hx)],
@@ -269,7 +242,7 @@ Keyword arguments
 
 - `file_splitting`: Schedule for splitting the output file. The new files will be suffixed with
           `_part1`, `_part2`, etc. For example `file_splitting = FileSizeLimit(sz)` will
-          split the output file when its size exceeds `sz`. Another example is 
+          split the output file when its size exceeds `sz`. Another example is
           `file_splitting = TimeInterval(30days)`, which will split files every 30 days of
           simulation time. The default incurs no splitting (`NoFileSplitting()`).
 
@@ -721,7 +694,7 @@ function start_next_file(model, ow::NetCDFOutputWriter)
     verbose && @info "Now writing to: $(ow.filepath)"
 
     initialize_nc_file!(ow, model)
-    
+
     return nothing
 end
 
@@ -785,7 +758,7 @@ function initialize_nc_file!(filepath,
 
         for (name, output) in outputs
             attributes = try output_attributes[name]; catch; Dict(); end
-            materialized = materialize_output(output, model) 
+            materialized = materialize_output(output, model)
             define_output_variable!(dataset,
                                     materialized,
                                     name,