diff --git a/.JuliaFormatter.toml b/.JuliaFormatter.toml
new file mode 100644
index 0000000..aa7f0eb
--- /dev/null
+++ b/.JuliaFormatter.toml
@@ -0,0 +1,27 @@
+align_assignment = true
+align_conditional = true
+align_matrix = true
+align_pair_arrow = true
+align_struct_field = true
+always_for_in = true
+always_use_return = true
+annotate_untyped_fields_with_any = false
+conditional_to_if = true
+for_in_replacement = "∈"
+format_docstrings = true
+import_to_using = true
+indent = 4
+indent_submodule = true
+join_lines_based_on_source = true
+long_to_short_function_def = true
+margin = 88
+normalize_line_endings = "unix"
+pipe_to_function_call = true
+remove_extra_newlines = true
+separate_kwargs_with_semicolon = true
+surround_whereop_typeparameters = true
+trailing_comma = true
+whitespace_in_kwargs = false
+whitespace_ops_in_indices = false
+whitespace_typedefs = true
+yas_style_nesting = true
diff --git a/.github/workflows/format_check.yml b/.github/workflows/format_check.yml
new file mode 100644
index 0000000..4e46a44
--- /dev/null
+++ b/.github/workflows/format_check.yml
@@ -0,0 +1,36 @@
+name: Format Check
+
+on:
+  push:
+    branches: ["master", "dev", "format"]
+  pull_request:
+    branches: ["master", "dev"]
+jobs:
+  check:
+    runs-on: ${{ matrix.os }}
+    strategy:
+      matrix:
+        julia-version: [1.9.3]
+        julia-arch: [x86]
+        os: [ubuntu-latest]
+    steps:
+      - uses: julia-actions/setup-julia@latest
+        with:
+          version: ${{ matrix.julia-version }}
+
+      - uses: actions/checkout@v1
+      - name: Install JuliaFormatter and format
+        run: |
+          julia  -e 'using Pkg; Pkg.add(PackageSpec(name="JuliaFormatter"))'
+          julia  -e 'using JuliaFormatter; format(".", verbose=true)'
+      - name: Format check
+        run: |
+          julia -e '
+          out = Cmd(`git diff --name-only`) |> read |> String
+          if out == ""
+              exit(0)
+          else
+              @error "Some files have not been formatted !!!"
+              write(stdout, out)
+              exit(1)
+          end'
diff --git a/src/SOLPS2IMAS.jl b/src/SOLPS2IMAS.jl
index 2e0f3a0..130001b 100644
--- a/src/SOLPS2IMAS.jl
+++ b/src/SOLPS2IMAS.jl
@@ -1,7 +1,7 @@
 module SOLPS2IMAS
 
 using Revise
-import OMAS
+using OMAS: OMAS
 using NCDatasets: Dataset, dimnames
 using YAML: load_file as YAML_load_file
 export try_omas
@@ -10,7 +10,6 @@ export read_b2_output
 export search_points
 export solps2imas
 
-
 function try_omas()
     ids = OMAS.dd()
     resize!(ids.equilibrium.time_slice, 1)
@@ -18,7 +17,6 @@ function try_omas()
     return nothing
 end
 
-
 function generate_test_data(nx=94, ny=38, sep=15, lcut=25, rcut=69)
     """This doesn't work well"""
     center_r = 1.7
@@ -27,9 +25,9 @@ function generate_test_data(nx=94, ny=38, sep=15, lcut=25, rcut=69)
     inner_minor = 0.3
     outer_minor = 0.8
     elongation = 1.8
-    θₓ = 8 * π/6
-    θᵢ = θₓ - π/4
-    θₒ = θₓ + π/4
+    θₓ = 8 * π / 6
+    θᵢ = θₓ - π / 4
+    θₒ = θₓ + π / 4
     inner_leg_length = 0.5
     outer_leg_length = 0.6
     dx_il = inner_leg_length / lcut
@@ -42,8 +40,8 @@ function generate_test_data(nx=94, ny=38, sep=15, lcut=25, rcut=69)
     cry = Array{Float64}(undef, 4, ny, nx)
     xpoint_r = center_r + aminor * cos(θₓ)
     xpoint_z = center_z + aminor * sin(θₓ) * elongation
-    for ix::Int = 1:nx
-        for iy::Int = 1:ny
+    for ix::Int ∈ 1:nx
+        for iy::Int ∈ 1:ny
             if iy < sep
                 dy = dr_in
             else
@@ -57,11 +55,15 @@ function generate_test_data(nx=94, ny=38, sep=15, lcut=25, rcut=69)
                 cell_centers_r[iy, ix] = cos(θ) * a + center_r
                 cell_centers_z[iy, ix] = sin(θ) * a * elongation + center_z
             elseif ix < lcut
-                cell_centers_r[iy, ix] = xpoint_r + cos(θᵢ) * dx_il * (lcut - ix) + sin(θᵢ) * dy * (iy - sep)
-                cell_centers_z[iy, ix] = xpoint_z + sin(θᵢ) * dx_il * (lcut - ix) + cos(θᵢ) * dy * (iy - sep)
+                cell_centers_r[iy, ix] =
+                    xpoint_r + cos(θᵢ) * dx_il * (lcut - ix) + sin(θᵢ) * dy * (iy - sep)
+                cell_centers_z[iy, ix] =
+                    xpoint_z + sin(θᵢ) * dx_il * (lcut - ix) + cos(θᵢ) * dy * (iy - sep)
             else
-                cell_centers_r[iy, ix] = xpoint_r + cos(θₒ) * dx_ol * (ix - rcut) + sin(θₒ) * dy * (iy - sep)
-                cell_centers_z[iy, ix] = xpoint_z + sin(θₒ) * dx_ol * (ix - rcut) + cos(θₒ) * dy * (iy - sep)
+                cell_centers_r[iy, ix] =
+                    xpoint_r + cos(θₒ) * dx_ol * (ix - rcut) + sin(θₒ) * dy * (iy - sep)
+                cell_centers_z[iy, ix] =
+                    xpoint_z + sin(θₒ) * dx_ol * (ix - rcut) + cos(θₒ) * dy * (iy - sep)
             end
         end
     end
@@ -70,7 +72,7 @@ function generate_test_data(nx=94, ny=38, sep=15, lcut=25, rcut=69)
     rmax = 3.0
     zmin = -3.0
     zmax = 3.0
-    for i::Int = 1:nx*ny
+    for i::Int ∈ 1:nx*ny
         # if mod(i, nx) == 0
         #     print(cell_centers_r[i], " ", cell_centers_r[i] > rmax, " ", rmax)
         # end
@@ -85,7 +87,6 @@ function generate_test_data(nx=94, ny=38, sep=15, lcut=25, rcut=69)
     return (cell_centers_r, cell_centers_z)
 end
 
-
 function read_b2time_output(filename)
     dim_order = (
         "time",
@@ -98,13 +99,17 @@ function read_b2time_output(filename)
     )
     ret_dict = Dict("dim" => Dict(), "data" => Dict())
     ds = Dataset(filename)
-    for key in keys(ds.dim)
+    for key ∈ keys(ds.dim)
         ret_dict["dim"][key] = ds.dim[key]
     end
-    for key in keys(ds)
+    for key ∈ keys(ds)
         if key != "ntstep"
             d = dimnames(ds[key])
-            permute = [y for y in [findfirst(x->x==dimord, d) for dimord in dim_order] if y !== nothing]
+            permute = [
+                y for
+                y ∈ [findfirst(x -> x == dimord, d) for dimord ∈ dim_order] if
+                y !== nothing
+            ]
             try
                 ret_dict["data"][key] = permutedims(Array(ds[key]), permute)
             catch e
@@ -117,13 +122,12 @@ function read_b2time_output(filename)
     return ret_dict
 end
 
-
 function read_b2mn_output(filename)
     lines = open(filename) do f
-        readlines(f)
+        return readlines(f)
     end
     contents = Dict()
-    for line in lines
+    for line ∈ lines
         if startswith(line, "'")
             splits = split(line, "'"; keepempty=false)
             contents[splits[1]] = parse(Float64, splits[end])
@@ -132,10 +136,7 @@ function read_b2mn_output(filename)
     return contents
 end
 
-
 function read_b2_output(filename)
-    # x = readdlm(filename, ' ', Float64, '\n', header=true)  # Doesn't handle the text lines at the start of each array
-
     if cmp(splitext(filename)[2], ".nc") == 0
         return read_b2time_output(filename)
     end
@@ -144,7 +145,7 @@ function read_b2_output(filename)
     array_sizes = Dict()
     ret_dict = Dict()
     lines = open(filename) do f
-        readlines(f)
+        return readlines(f)
     end
     nx = 0
     ny = 0
@@ -153,7 +154,7 @@ function read_b2_output(filename)
     arraysize = 0
     arraytype = nothing
     j = 1
-    for l in lines
+    for l ∈ lines
         if startswith(l, "*cf:")
             j = 1  # Reset intra-array element counter
             _, arraytype, arraysize, tag = split(l)
@@ -168,10 +169,10 @@ function read_b2_output(filename)
             array_sizes[tag] = arraysize
         elseif tag != ""
             if arraytype == "int"
-                array_line = [parse(Int, ss) for ss in split(l)]
+                array_line = [parse(Int, ss) for ss ∈ split(l)]
                 array_inc = size(array_line)[1]
             elseif arraytype == "real"
-                array_line = [parse(Float64, ss) for ss in split(l)]
+                array_line = [parse(Float64, ss) for ss ∈ split(l)]
                 array_inc = size(array_line)[1]
             else
                 array_line = l
@@ -190,12 +191,14 @@ function read_b2_output(filename)
     elseif "nx,ny,ns" ∈ keys(contents)
         return extract_state_quantities(contents)
     else
-        throw(DomainError(keys(contents),
-              "nx,ny (b2fgmtry) or nx,ny,ns (b2fstate) must be present in b2 output file."))
+        throw(
+            DomainError(keys(contents),
+                "nx,ny (b2fgmtry) or nx,ny,ns (b2fstate) must be present in b2 output file.",
+            ),
+        )
     end
 end
 
-
 function extract_geometry(gmtry)
     ret_dict = Dict("dim" => Dict(), "data" => Dict())
     ret_dict["dim"]["nx_no_guard"], ret_dict["dim"]["ny_no_guard"] = gmtry["nx,ny"]
@@ -208,12 +211,13 @@ function extract_geometry(gmtry)
     # Adding placeholder timestamp
     ret_dict["dim"]["time"] = 1
     ret_dict["data"]["timesa"] = [0.0]
-    for k in keys(gmtry)
+    for k ∈ keys(gmtry)
         # The 4 fields of bb are poloidal, radial, toroidal, and total magnetic field
         # according to page 212 of D. Coster, "SOLPS-ITER [manual]" (2019)
         # The 4 fields in crx and cry are the corners of each grid cell.
         if k ∈ ["crx", "cry", "bb"]
-            ret_dict["data"][k] = permutedims(reshape(gmtry[k], (nx, ny, 4, 1)), (4, 3, 2, 1))
+            ret_dict["data"][k] =
+                permutedims(reshape(gmtry[k], (nx, ny, 4, 1)), (4, 3, 2, 1))
         elseif k ∈ ["leftcut", "bottomcut", "rightcut", "topcut"]
             ret_dict["data"][k] = Array([gmtry[k][1]])
         elseif k ∈ ["leftcut2", "bottomcut2", "rightcut2", "topcut2"]
@@ -227,10 +231,11 @@ function extract_geometry(gmtry)
     return ret_dict
 end
 
-
 function extract_state_quantities(state)
     ret_dict = Dict("dim" => Dict(), "data" => Dict())
-    ret_dict["dim"]["nx_no_guard"], ret_dict["dim"]["ny_no_guard"], ret_dict["dim"]["ns"] = state["nx,ny,ns"]
+    ret_dict["dim"]["nx_no_guard"],
+    ret_dict["dim"]["ny_no_guard"],
+    ret_dict["dim"]["ns"] = state["nx,ny,ns"]
     # includes guard cells
     nx = ret_dict["dim"]["nx"] = ret_dict["dim"]["nx_no_guard"] + 2
     ny = ret_dict["dim"]["ny"] = ret_dict["dim"]["ny_no_guard"] + 2
@@ -239,16 +244,19 @@ function extract_state_quantities(state)
     # Adding placeholder timestamp
     ret_dict["dim"]["time"] = 1
     ret_dict["data"]["timesa"] = [0.0]
-    for k in keys(state)
+    for k ∈ keys(state)
         l = length(state[k])
         if l == nx * ny
             ret_dict["data"][k] = permutedims(reshape(state[k], (nx, ny, 1)), (3, 2, 1))
         elseif l == nx * ny * ns
-            ret_dict["data"][k] = permutedims(reshape(state[k], (nx, ny, ns, 1)), (4, 3, 2, 1))
+            ret_dict["data"][k] =
+                permutedims(reshape(state[k], (nx, ny, ns, 1)), (4, 3, 2, 1))
         elseif l == nx * ny * ndir
-            ret_dict["data"][k] = permutedims(reshape(state[k], (nx, ny, ndir, 1)), (4, 3, 2, 1))
+            ret_dict["data"][k] =
+                permutedims(reshape(state[k], (nx, ny, ndir, 1)), (4, 3, 2, 1))
         elseif l == nx * ny * ndir * ns
-            ret_dict["data"][k] = permutedims(reshape(state[k], (nx, ny, ndir, ns, 1)), (5, 4, 3, 2, 1))
+            ret_dict["data"][k] =
+                permutedims(reshape(state[k], (nx, ny, ndir, ns, 1)), (5, 4, 3, 2, 1))
         elseif l == ns
             ret_dict["data"][k] = permutedims(reshape(state[k], (ns, 1)), (2, 1))
         elseif k ∉ keys(ret_dict["dim"])
@@ -258,7 +266,6 @@ function extract_state_quantities(state)
     return ret_dict
 end
 
-
 """
     xytoc(ix, iy; nx)
 
@@ -270,7 +277,6 @@ function xytoc(ix, iy; nx)
     return ic
 end
 
-
 """
     ctoxy(ic; nx)
 
@@ -282,9 +288,8 @@ function ctoxy(ic; nx)
     return ix, iy
 end
 
-
 """
-   in_core(; ix, iy, topcut, bottomcut, leftcut, rightcut)
+in_core(; ix, iy, topcut, bottomcut, leftcut, rightcut)
 
 Returns true if cell indexed ix, iy lie inside the core
 """
@@ -292,16 +297,12 @@ function in_core(; ix, iy, topcut, bottomcut, leftcut, rightcut)
     return bottomcut + 1 < iy < topcut + 2 && leftcut + 1 < ix < rightcut + 2
 end
 
-
 """
     in_sol(; iy, topcut, kwargs...)
 
 Returns true if cell indexed ix, iy lie inside the SOL
 """
-function in_sol(; iy, topcut, kwargs...)
-    return topcut + 1 < iy
-end
-
+in_sol(; iy, topcut, kwargs...) = topcut + 1 < iy
 
 """
     in_idr(; ix, iy, topcut, bottomcut, leftcut, kwargs...)
@@ -312,7 +313,6 @@ function in_idr(; ix, iy, topcut, bottomcut, leftcut, kwargs...)
     return bottomcut + 1 < iy < topcut + 2 && ix < leftcut + 2
 end
 
-
 """
     in_odr(; ix, iy, topcut, bottomcut, rightcut, kwargs...)
 
@@ -328,52 +328,66 @@ end
 # Thus, x-edges will have odd boundary index and y_edges will have even
 # List of tuples (boundary_ind, (corner pair forming edge))
 chosen_edge_order = [(1, (1, 2)),
-                     (2, (2, 4)),
-                     (3, (4, 3)),
-                     (4, (3, 1))]
-
+    (2, (2, 4)),
+    (3, (4, 3)),
+    (4, (3, 1))]
 
 """
     is_x_aligned(;boundary_ind)
 
 y_aligned edges will have odd boundary_ind based on chosen order of numbering them
 """
-function is_x_aligned(;boundary_ind)
-    return mod(boundary_ind, 2) == 1
-end
-
+is_x_aligned(; boundary_ind) = mod(boundary_ind, 2) == 1
 
 """
     is_y_aligned(; boundary_ind)
 
 y_aligned edges will have even boundary_ind based on chosen order of numbering them
 """
-function is_y_aligned(; boundary_ind)
-    return mod(boundary_ind, 2) == 0
-end
+is_y_aligned(; boundary_ind) = mod(boundary_ind, 2) == 0
 
 """
-   is_core_cut(; ix, iy, nx, cells, boundary_ind, topcut, bottomcut, leftcut, rightcut)
+is_core_cut(; ix, iy, nx, cells, boundary_ind, topcut, bottomcut, leftcut, rightcut)
 
 Returns true if boundary_ind of a cell at ix, iy is on core_cut (Y-aliged edge)
 """
-function is_core_cut(; ix, iy, cells, nx, boundary_ind, topcut, bottomcut, leftcut, rightcut)
+function is_core_cut(;
+    ix,
+    iy,
+    cells,
+    nx,
+    boundary_ind,
+    topcut,
+    bottomcut,
+    leftcut,
+    rightcut,
+)
     if bottomcut + 1 < iy < topcut + 2 && ix == leftcut + 2 && mod(boundary_ind, 2) == 0
         ixr = rightcut + 1
         this_cell = cells[xytoc(ix, iy; nx=nx)]
-        other_side_cell_ind = xytoc(ixr, iy; nx=nx)  # Cell on the other side of core cut
+        # Cell on the other side of core cut
+        other_side_cell_ind = xytoc(ixr, iy; nx=nx)
         return other_side_cell_ind ∈ this_cell.boundary[boundary_ind].neighbours
     end
     return false
 end
 
-
 """
-   is_pfr_cut(; ix, iy, nx, cells, boundary_ind, topcut, bottomcut, leftcut, rightcut)
+is_pfr_cut(; ix, iy, nx, cells, boundary_ind, topcut, bottomcut, leftcut, rightcut)
 
 Returns true if boundary_ind of a cell at ix, iy is on core_cut (Y-aliged edge)
 """
-function is_pfr_cut(; ix, iy,  cells, nx, boundary_ind, topcut, bottomcut, leftcut, rightcut)
+function is_pfr_cut(;
+    ix,
+    iy,
+    cells,
+    nx,
+    boundary_ind,
+    topcut,
+    bottomcut,
+    leftcut,
+    rightcut,
+)
     if bottomcut + 1 < iy < topcut + 2 && ix == leftcut + 1 && mod(boundary_ind, 2) == 0
         ixr = rightcut + 2
         this_cell = cells[xytoc(ix, iy; nx=nx)]
@@ -383,7 +397,6 @@ function is_pfr_cut(; ix, iy,  cells, nx, boundary_ind, topcut, bottomcut, leftc
     return false
 end
 
-
 """
     is_outer_throat(; ix, iy, boundary_ind, topcut, rightcut, kwargs...)
 
@@ -393,7 +406,6 @@ function is_outer_throat(; ix, iy, boundary_ind, topcut, rightcut, kwargs...)
     return topcut + 1 < iy && ix == rightcut + 1 && boundary_ind == 2
 end
 
-
 """
     is_inner_throat(; ix, iy, boundary_ind, topcut, leftcut, kwargs...)
 
@@ -403,7 +415,6 @@ function is_inner_throat(; ix, iy, boundary_ind, topcut, leftcut, kwargs...)
     return topcut + 1 < iy && ix == leftcut + 2 && boundary_ind == 4
 end
 
-
 """
     is_outer_midplane(; ix, jxa, boundary_ind)
 
@@ -414,7 +425,6 @@ function is_outer_midplane(; ix, iy, jxa, boundary_ind, topcut, kwargs...)
     return ix == jxa && boundary_ind == 2
 end
 
-
 """
     is_inner_midplane(; ix, jxa, boundary_ind)
 
@@ -425,37 +435,47 @@ function is_inner_midplane(; ix, iy, jxi, boundary_ind, topcut, kwargs...)
     return ix == jxi && boundary_ind == 4
 end
 
-
 """
-    is_outer_target(; ix, nx, boundary_ind, kwargs...)
+    is_outer_target(; ix, nx, boundary_ind)
 
 Returns true if boundary_ind of a cell at ix, iy is on outer target
 """
-function is_outer_target(; ix, nx, boundary_ind)
-    return ix == nx && boundary_ind == 2
-end
-
+is_outer_target(; ix, nx, boundary_ind) = ix == nx && boundary_ind == 2
 
 """
     is_inner_target(; ix, boundary_ind, kwargs...)
 
 Returns true if boundary_ind of a cell at ix, iy is on inner target
 """
-function is_inner_target(; ix, boundary_ind)
-    return ix == 1 && boundary_ind == 4
-end
-
+is_inner_target(; ix, boundary_ind) = ix == 1 && boundary_ind == 4
 
 """
-    is_core_boundary(; ix, iy, boundary_ind, bottomcut, leftcut, rightcut, kwargs...)
-
-Returns true if boundary_ind of a cell at ix, iy is on core boundary (central blank spot boundary)
+    is_core_boundary(;
+    ix,
+    iy,
+    boundary_ind,
+    bottomcut,
+    leftcut,
+    rightcut,
+    kwargs...,
+
+)
+
+Returns true if boundary_ind of a cell at ix, iy is on core boundary (central blank
+spot boundary)
 """
-function is_core_boundary(; ix, iy, boundary_ind, bottomcut, leftcut, rightcut, kwargs...)
+function is_core_boundary(;
+    ix,
+    iy,
+    boundary_ind,
+    bottomcut,
+    leftcut,
+    rightcut,
+    kwargs...,
+)
     return bottomcut + 2 == iy && leftcut + 1 < ix < rightcut + 2 && boundary_ind == 1
 end
 
-
 """
     is_separatix(; iy, boundary_ind, topcut, kwargs...)
 
@@ -465,7 +485,6 @@ function is_separatix(; iy, boundary_ind, topcut, kwargs...)
     return topcut + 2 == iy && boundary_ind == 1
 end
 
-
 """
     add_subset_element!(subset, sn, dim, index::Int, in_subset=(x...)->true; kwargs...)
 
@@ -474,7 +493,14 @@ with space number sn, dimension dim, index. To determine,
 if the element should be added or not, a function in_subset can be provided that gets the arguments
 (kwargs...). These functions will be in_core, in_sol etc as difined above.
 """
-function add_subset_element!(subset, sn, dim, index::Int, in_subset=(x...)->true; kwargs...)
+function add_subset_element!(
+    subset,
+    sn,
+    dim,
+    index::Int,
+    in_subset=(x...) -> true;
+    kwargs...,
+)
     if in_subset(; kwargs...)
         dd_ind = length(subset.element) + 1
         resize!(subset.element, dd_ind)
@@ -485,18 +511,32 @@ function add_subset_element!(subset, sn, dim, index::Int, in_subset=(x...)->true
     end
 end
 
-
 """
-    add_subset_element!(subset, sn, dim, index::Vector{Int}, in_subset=(x...)->true; kwargs...)
+    add_subset_element!(
+    subset,
+    sn,
+    dim,
+    index::Vector{Int},
+    in_subset=(x...) -> true;
+    kwargs...,
+
+)
 
 Overloaded to work differently (faster) with list of indices to be added.
 """
-function add_subset_element!(subset, sn, dim, index::Vector{Int}, in_subset=(x...)->true; kwargs...)
+function add_subset_element!(
+    subset,
+    sn,
+    dim,
+    index::Vector{Int},
+    in_subset=(x...) -> true;
+    kwargs...,
+)
     if in_subset(; kwargs...)
         dd_start_ind = length(subset.element) + 1
         resize!(subset.element, length(subset.element) + length(index))
         dd_stop_ind = length(subset.element)
-        for (ii, dd_ind) in enumerate(dd_start_ind:dd_stop_ind)
+        for (ii, dd_ind) ∈ enumerate(dd_start_ind:dd_stop_ind)
             resize!(subset.element[dd_ind].object, 1)
             subset.element[dd_ind].object[1].space = sn
             subset.element[dd_ind].object[1].dimension = dim
@@ -505,7 +545,6 @@ function add_subset_element!(subset, sn, dim, index::Vector{Int}, in_subset=(x..
     end
 end
 
-
 """
     get_xpoint_nodes(gmtry)
 
@@ -517,23 +556,27 @@ function get_xpoint_nodes(gmtry)
     nt = gmtry["dim"]["time"]
     # Find cells around x-point
     xpcells = [(gmtry["data"]["topcut"][1] + 1, gmtry["data"]["leftcut"][1] + 1),
-               (gmtry["data"]["topcut"][1] + 1, gmtry["data"]["leftcut"][1] + 2),
-               (gmtry["data"]["topcut"][1] + 2, gmtry["data"]["leftcut"][1] + 1),
-               (gmtry["data"]["topcut"][1] + 2, gmtry["data"]["leftcut"][1] + 2),
-               (gmtry["data"]["topcut"][1] + 1, gmtry["data"]["rightcut"][1] + 1),
-               (gmtry["data"]["topcut"][1] + 1, gmtry["data"]["rightcut"][1] + 2),
-               (gmtry["data"]["topcut"][1] + 2, gmtry["data"]["rightcut"][1] + 1),
-               (gmtry["data"]["topcut"][1] + 2, gmtry["data"]["rightcut"][1] + 2)]
+        (gmtry["data"]["topcut"][1] + 1, gmtry["data"]["leftcut"][1] + 2),
+        (gmtry["data"]["topcut"][1] + 2, gmtry["data"]["leftcut"][1] + 1),
+        (gmtry["data"]["topcut"][1] + 2, gmtry["data"]["leftcut"][1] + 2),
+        (gmtry["data"]["topcut"][1] + 1, gmtry["data"]["rightcut"][1] + 1),
+        (gmtry["data"]["topcut"][1] + 1, gmtry["data"]["rightcut"][1] + 2),
+        (gmtry["data"]["topcut"][1] + 2, gmtry["data"]["rightcut"][1] + 1),
+        (gmtry["data"]["topcut"][1] + 2, gmtry["data"]["rightcut"][1] + 2)]
     # Get list of all nodes in these cells
     candidate_nodes = []
     resize!(candidate_nodes, nt)
-    for it in 1:nt
-        candidate_nodes[it] = [[[crx[it, icorner, iy, ix], cry[it, icorner, iy, ix]] for icorner in 1:4] for (iy, ix) in xpcells]
+    for it ∈ 1:nt
+        candidate_nodes[it] = [
+            [
+                [crx[it, icorner, iy, ix], cry[it, icorner, iy, ix]] for icorner ∈ 1:4
+            ] for (iy, ix) ∈ xpcells
+        ]
     end
     xpoint_nodes = []
     resize!(xpoint_nodes, nt)
     # Find the node that is common among all the cells
-    for it in 1:nt
+    for it ∈ 1:nt
         xpoint_nodes[it] = intersect(candidate_nodes[it]...)
     end
     return xpoint_nodes
@@ -546,9 +589,10 @@ function search_points(ids, r, z)
     space_number = 1
     subset_idx_node = 1
     # If an index remains at 0, it means the point in question was not found
-    nodes = ids.edge_profiles.grid_ggd[grid_number].space[space_number].objects_per_dimension[subset_idx_node].object
-    for j = 1:n
-        for i in eachindex(nodes)
+    nodes =
+        ids.edge_profiles.grid_ggd[grid_number].space[space_number].objects_per_dimension[subset_idx_node].object
+    for j ∈ 1:n
+        for i ∈ eachindex(nodes)
             rn = nodes[i].geometry[1]
             zn = nodes[i].geometry[2]
             if (rn == r[j]) && (zn == z[j])
@@ -560,36 +604,34 @@ function search_points(ids, r, z)
     return indices
 end
 
-
 """
     search_edges(edges, edge_nodes)
 
 search if an edge with nodes as edge_nodes already exists
 """
 function search_edges(edges, edge_nodes)
-    for ii in eachindex(edges)
+    for ii ∈ eachindex(edges)
         if edge_nodes[2] == edges[ii].nodes[1] && edge_nodes[1] == edges[ii].nodes[2]
             return ii
-        elseif edge_nodes[2] == edges[ii].nodes[1] && edge_nodes[1] == edges[ii].nodes[2]
+        elseif edge_nodes[2] == edges[ii].nodes[1] &&
+               edge_nodes[1] == edges[ii].nodes[2]
             return ii
         end
     end
     return 0
 end
 
-
 """
     distance_between_nodes(nodes, node_inds)
 
 Return distance between two node indices
 """
 function distance_between_nodes(nodes, node_inds)
-    return √(sum((nodes[node_inds[1]].geometry - nodes[node_inds[2]].geometry).^2))
+    return √(sum((nodes[node_inds[1]].geometry - nodes[node_inds[2]].geometry) .^ 2))
 end
 
-
 function neighbour_inds(ic; nx, ny, leftcut, rightcut, topcut, bottomcut)
-    ix, iy = ctoxy(ic, nx=nx)
+    ix, iy = ctoxy(ic; nx=nx)
     neighbour_x_inds = []
     neighbour_y_inds = []
     if ix > 1
@@ -624,25 +666,24 @@ function neighbour_inds(ic; nx, ny, leftcut, rightcut, topcut, bottomcut)
     if iy < ny
         append!(neighbour_y_inds, iy + 1)
     end
-        
+
     neighbour_inds = []
-    for x_ind in neighbour_x_inds
+    for x_ind ∈ neighbour_x_inds
         append!(neighbour_inds, xytoc(x_ind, iy; nx=nx))
     end
-    for y_ind in neighbour_y_inds
+    for y_ind ∈ neighbour_y_inds
         append!(neighbour_inds, xytoc(ix, y_ind; nx=nx))
     end
     return neighbour_inds
 end
 
-
 function get_neighbour_inds(ic, gmtry, it)
     nx = gmtry["dim"]["nx"]
     ny = gmtry["dim"]["ny"]
-    ix, iy = ctoxy(ic, nx=nx)
+    ix, iy = ctoxy(ic; nx=nx)
     neighbour_inds = []
     # println(ix, ", ", iy)
-    for neighbour in ["left", "right", "top", "bottom"]
+    for neighbour ∈ ["left", "right", "top", "bottom"]
         nix = gmtry["data"][neighbour*"ix"][it, iy, ix] + 2
         niy = gmtry["data"][neighbour*"iy"][it, iy, ix] + 2
         # println(neighbour, ": ", nix, ", ", niy)
@@ -653,30 +694,34 @@ function get_neighbour_inds(ic, gmtry, it)
     return neighbour_inds
 end
 
-
 function attach_neightbours(cells, edges, gmtry, it)
-    for (ic, cell) in enumerate(cells)
-        for neighbour_ind in get_neighbour_inds(ic, gmtry, it)
-            for boundary in cell.boundary
-                for neighbour_boundary in cells[neighbour_ind].boundary
-                    if boundary.index == neighbour_boundary.index && neighbour_ind ∉ boundary.neighbours
+    for (ic, cell) ∈ enumerate(cells)
+        for neighbour_ind ∈ get_neighbour_inds(ic, gmtry, it)
+            for boundary ∈ cell.boundary
+                for neighbour_boundary ∈ cells[neighbour_ind].boundary
+                    if boundary.index == neighbour_boundary.index &&
+                       neighbour_ind ∉ boundary.neighbours
                         append!(boundary.neighbours, neighbour_ind)
                     end
                 end
             end
         end
     end
-    for (ic, cell) in enumerate(cells)
-        for edge_ind in [bnd.index for bnd in cell.boundary]
-            neighbour_edge_inds = [bnd.index for bnd in cell.boundary]
-            for neighbour_ind in get_neighbour_inds(ic, gmtry, it)
-                union!(neighbour_edge_inds, [bnd.index for bnd in cells[neighbour_ind].boundary])
+    for (ic, cell) ∈ enumerate(cells)
+        for edge_ind ∈ [bnd.index for bnd ∈ cell.boundary]
+            neighbour_edge_inds = [bnd.index for bnd ∈ cell.boundary]
+            for neighbour_ind ∈ get_neighbour_inds(ic, gmtry, it)
+                union!(
+                    neighbour_edge_inds,
+                    [bnd.index for bnd ∈ cells[neighbour_ind].boundary],
+                )
             end
             setdiff!(neighbour_edge_inds, edge_ind)
-            for neighbour_edge_ind in neighbour_edge_inds
-                for edge_bnd in edges[edge_ind].boundary
-                    for neighbour_edge_bnd in edges[neighbour_edge_ind].boundary
-                        if edge_bnd.index == neighbour_edge_bnd.index && neighbour_edge_ind ∉ edge_bnd.neighbours
+            for neighbour_edge_ind ∈ neighbour_edge_inds
+                for edge_bnd ∈ edges[edge_ind].boundary
+                    for neighbour_edge_bnd ∈ edges[neighbour_edge_ind].boundary
+                        if edge_bnd.index == neighbour_edge_bnd.index &&
+                           neighbour_edge_ind ∉ edge_bnd.neighbours
                             append!(edge_bnd.neighbours, neighbour_edge_ind)
                         end
                     end
@@ -686,14 +731,13 @@ function attach_neightbours(cells, edges, gmtry, it)
     end
 end
 
-
 """
     dict2prop!(obj, dict)
 
 Copies grid_ggd and space description in dict format to the data structure recursively.
 """
-function dict2prop!(obj, dict)
-    for (key, prop) in dict
+dict2prop!(obj, dict) =
+    for (key, prop) ∈ dict
         if isa(key, Int)
             if length(obj) < key
                 resize!(obj, key)
@@ -707,8 +751,6 @@ function dict2prop!(obj, dict)
             end
         end
     end
-end
-
 
 """
     path_to_obj(obj, path)
@@ -724,7 +766,7 @@ Note: If integer is -1, the array of field is resized to
 increase by 1 in length and last element is returned.
 """
 function path_to_obj(obj, path)
-    for ele in path
+    for ele ∈ path
         if typeof(ele) == String
             obj = getfield(obj, Symbol(ele))
         elseif typeof(ele) == Int
@@ -740,32 +782,28 @@ function path_to_obj(obj, path)
     return obj
 end
 
-
 isint(x) = typeof(x) == Int
 
-
 """
     val_obj(var, ggd, grid_ggd_index)
 
 Given SOLPS variable name (var), return the field to write values on from
-ggd object. 
+ggd object.
 """
 solps_var_to_imas = YAML_load_file("$(@__DIR__)/solps_var_to_imas.yml")
-function val_obj(ggd, var, grid_ggd_index)
+val_obj(ggd, var, grid_ggd_index) =
     if var ∉ keys(solps_var_to_imas)
         return nothing
     else
         path = solps_var_to_imas[var]
         gsi_ind = findlast(isint, path)
         path_to_prop = path[1:gsi_ind]
-        prop_to_obj = path[gsi_ind + 1:end]
+        prop_to_obj = path[gsi_ind+1:end]
         prop = path_to_obj(ggd, path_to_prop)
         prop.grid_index = grid_ggd_index
         prop.grid_subset_index = gsi_ind
         return path_to_obj(prop, prop_to_obj)
     end
-end
-
 
 """
     find_subset_index()
@@ -778,7 +816,7 @@ We'll call them dd_index (the place in the DD) and ggd_index (the type of subset
 """
 function find_subset_index(gsdesc, ggd_index)
     subsets = length(gsdesc["grid_subset"])
-    for dd_index = 1:subsets
+    for dd_index ∈ 1:subsets
         if gsdesc["grid_subset"][dd_index]["identifier"]["index"] == ggd_index
             return dd_index
         end
@@ -786,14 +824,13 @@ function find_subset_index(gsdesc, ggd_index)
     return 0  # Indicates failure
 end
 
-
 """
     get_grid_subset_with_index(grid_subset_index)
 
 Returns the grid_subset in a grid_ggd with the matching grid_subset_index
 """
 function get_grid_subset_with_index(grid_ggd, grid_subset_index)
-    for subset in grid_ggd.grid_subset
+    for subset ∈ grid_ggd.grid_subset
         if subset.identifier.index == grid_subset_index
             return subset
         end
@@ -801,24 +838,27 @@ function get_grid_subset_with_index(grid_ggd, grid_subset_index)
     return 0  # Indicates failure
 end
 
-
-function get_subset_boundary_inds(space::OMAS.edge_profiles__grid_ggd___space, subset::OMAS.edge_profiles__grid_ggd___grid_subset)
+function get_subset_boundary_inds(
+    space::OMAS.edge_profiles__grid_ggd___space,
+    subset::OMAS.edge_profiles__grid_ggd___grid_subset,
+)
     nD = subset.element[1].object[1].dimension
     if nD > 0
-        nD_objects = space.objects_per_dimension[nD + 1].object
-        elements = [nD_objects[ele.object[1].index] for ele in subset.element]
+        nD_objects = space.objects_per_dimension[nD+1].object
+        elements = [nD_objects[ele.object[1].index] for ele ∈ subset.element]
         boundary_inds = Int[]
-        for ele in elements
-            symdiff!(boundary_inds, [bnd.index for bnd in ele.boundary])
-
+        for ele ∈ elements
+            symdiff!(boundary_inds, [bnd.index for bnd ∈ ele.boundary])
         end
         return boundary_inds
     end
     return []
 end
 
-
-function get_subset_boundary(space::OMAS.edge_profiles__grid_ggd___space, subset::OMAS.edge_profiles__grid_ggd___grid_subset)
+function get_subset_boundary(
+    space::OMAS.edge_profiles__grid_ggd___space,
+    subset::OMAS.edge_profiles__grid_ggd___grid_subset,
+)
     ret_subset = OMAS.edge_profiles__grid_ggd___grid_subset()
     boundary_inds = get_subset_boundary_inds(space, subset)
     bnd_dim = subset.element[1].object[1].dimension - 1
@@ -827,33 +867,41 @@ function get_subset_boundary(space::OMAS.edge_profiles__grid_ggd___space, subset
     return ret_subset.element
 end
 
-
 function get_subset_space(space::OMAS.edge_profiles__grid_ggd___space,
-                          elements::Vector{OMAS.edge_profiles__grid_ggd___grid_subset___element})
+    elements::Vector{OMAS.edge_profiles__grid_ggd___grid_subset___element})
     nD = elements[1].object[1].dimension
     nD_objects = space.objects_per_dimension[nD+1].object
-    return [nD_objects[ele.object[1].index] for ele in elements]
+    return [nD_objects[ele.object[1].index] for ele ∈ elements]
 end
 
 """
-    subset_do(set_operator, itrs::Vararg{Vector{OMAS.edge_profiles__grid_ggd___grid_subset___element}};
-              space::OMAS.edge_profiles__grid_ggd___space=OMAS.edge_profiles__grid_ggd___space(),
-              use_nodes=false)
-
-Function to perform any set operation (intersect, union, setdiff etc.) on subset.element to
-generate a list of elements to go to subset object. If use_nodes is true, the set operation will be
-applied on the set of nodes from subset.element, space argument is required for this.
-Note: that the arguments are subset.element (not the subset itself). Similarly, the return object is a
-list of OMAS.edge_profiles__grid_ggd___grid_subset___element.
+    subset_do(set_operator,
+    itrs::Vararg{Vector{OMAS.edge_profiles__grid_ggd___grid_subset___element}};
+    space::OMAS.edge_profiles__grid_ggd___space=OMAS.edge_profiles__grid_ggd___space(),
+    use_nodes=false)
+
+Function to perform any set operation (intersect, union, setdiff etc.) on
+subset.element to generate a list of elements to go to subset object. If use_nodes is
+true, the set operation will be applied on the set of nodes from subset.element, space
+argument is required for this.
+Note: that the arguments are subset.element (not the subset itself). Similarly, the
+return object is a list of OMAS.edge_profiles__grid_ggd___grid_subset___element.
 """
-function subset_do(set_operator, itrs::Vararg{Vector{OMAS.edge_profiles__grid_ggd___grid_subset___element}};
-                   space::OMAS.edge_profiles__grid_ggd___space=OMAS.edge_profiles__grid_ggd___space(),
-                   use_nodes=false)
+function subset_do(set_operator,
+    itrs::Vararg{Vector{OMAS.edge_profiles__grid_ggd___grid_subset___element}};
+    space::OMAS.edge_profiles__grid_ggd___space=OMAS.edge_profiles__grid_ggd___space(),
+    use_nodes=false)
     if use_nodes
-        ele_inds = set_operator([union([obj.nodes for obj in get_subset_space(space, set_elements)]...) for set_elements in itrs]...)
+        ele_inds = set_operator(
+            [
+                union([obj.nodes for obj ∈ get_subset_space(space, set_elements)]...) for set_elements ∈ itrs
+            ]...,
+        )
         dim = 0
     else
-        ele_inds = set_operator([[ele.object[1].index for ele in set_elements] for set_elements in itrs]...)
+        ele_inds = set_operator(
+            [[ele.object[1].index for ele ∈ set_elements] for set_elements ∈ itrs]...,
+        )
         dim = itrs[1][1].object[1].dimension
     end
     ret_subset = OMAS.edge_profiles__grid_ggd___grid_subset()
@@ -862,7 +910,6 @@ function subset_do(set_operator, itrs::Vararg{Vector{OMAS.edge_profiles__grid_gg
     return ret_subset.element
 end
 
-
 """
     solps2imas(b2gmtry, b2output, gsdesc; load_bb=false)
 
@@ -897,62 +944,61 @@ function solps2imas(b2gmtry, b2output, gsdesc, b2mn=nothing; load_bb=false)
     cut_keys = ["leftcut", "rightcut", "bottomcut", "topcut"]
     cuts_found = cut_keys ⊆ keys(gmtry["data"])
     if cuts_found
-        cuts = Dict([(Symbol(key), gmtry["data"][key][1]) for key in cut_keys])
+        cuts = Dict([(Symbol(key), gmtry["data"][key][1]) for key ∈ cut_keys])
         xpoints_nodes = get_xpoint_nodes(gmtry)
     end
 
-
     if typeof(gsdesc) == String
         gsdesc = YAML_load_file(gsdesc)
     end
 
     # Add grid_ggd array equal to number of time steps
     resize!(ids.edge_profiles.grid_ggd, gmtry["dim"]["time"])
-    for it in 1:gmtry["dim"]["time"]
-        
+    for it ∈ 1:gmtry["dim"]["time"]
         grid_ggd = ids.edge_profiles.grid_ggd[it]
         grid_ggd.time = Float64.(gmtry["data"]["timesa"][it])
         dict2prop!(grid_ggd, gsdesc)
-        for sn in keys(gsdesc["space"])
+        for sn ∈ keys(gsdesc["space"])
             space = grid_ggd.space[sn]
-            # Assuming following to be standard for now. We can add this info through YAML as well
+            # Assuming following to be standard for now.
+            # We can add this info through YAML as well
             resize!(space.objects_per_dimension, 4)
             o0 = space.objects_per_dimension[1]  # 0D objects
             o1 = space.objects_per_dimension[2]  # 1D objects
             o2 = space.objects_per_dimension[3]  # 2D objects
             o3 = space.objects_per_dimension[4]  # 3D objects
 
-            subset_nodes = get_grid_subset_with_index(grid_ggd, 1)  # nodes have index 1
-            subset_faces = get_grid_subset_with_index(grid_ggd, 2)  # faces (edges with elongation in third dimension) have index 2
-            subset_xedges = get_grid_subset_with_index(grid_ggd, 3)  # faces (edges with elongation in third dimension) have index 2
-            subset_yedges = get_grid_subset_with_index(grid_ggd, 4)  # faces (edges with elongation in third dimension) have index 2
-            subset_cells = get_grid_subset_with_index(grid_ggd, 5)  # cells (cell in 2D grid, volume with elongation) have index 5
+            subset_nodes = get_grid_subset_with_index(grid_ggd, 1)
+            subset_faces = get_grid_subset_with_index(grid_ggd, 2)
+            subset_xedges = get_grid_subset_with_index(grid_ggd, 3)
+            subset_yedges = get_grid_subset_with_index(grid_ggd, 4)
+            subset_cells = get_grid_subset_with_index(grid_ggd, 5)
             if cuts_found
-                subset_core = get_grid_subset_with_index(grid_ggd, 22)  # core cells have index 22
-                subset_sol = get_grid_subset_with_index(grid_ggd, 23)   # sol cells have index 23
-                subset_odr = get_grid_subset_with_index(grid_ggd, 24)   # odr cells have index 24
-                subset_idr = get_grid_subset_with_index(grid_ggd, 25)   # idr cells have index 25
-                subset_xp = get_grid_subset_with_index(grid_ggd, 6)     # x_points have index 6
-                subset_corecut = get_grid_subset_with_index(grid_ggd, 7)     # core_cut have index 7
-                subset_pfrcut = get_grid_subset_with_index(grid_ggd, 8)      # pfr_cut have index 8
-                subset_othroat = get_grid_subset_with_index(grid_ggd, 9)     # outer_throat index 9
-                subset_ithroat = get_grid_subset_with_index(grid_ggd, 10)    # inner_throat index 10
+                subset_core = get_grid_subset_with_index(grid_ggd, 22)
+                subset_sol = get_grid_subset_with_index(grid_ggd, 23)
+                subset_odr = get_grid_subset_with_index(grid_ggd, 24)
+                subset_idr = get_grid_subset_with_index(grid_ggd, 25)
+                subset_xp = get_grid_subset_with_index(grid_ggd, 6)
+                subset_corecut = get_grid_subset_with_index(grid_ggd, 7)
+                subset_pfrcut = get_grid_subset_with_index(grid_ggd, 8)
+                subset_othroat = get_grid_subset_with_index(grid_ggd, 9)
+                subset_ithroat = get_grid_subset_with_index(grid_ggd, 10)
                 if !isnothing(jxa)
-                    subset_omp = get_grid_subset_with_index(grid_ggd, 11)     # outer midplane 11
-                    subset_ompsep = get_grid_subset_with_index(grid_ggd, 101) # outer midplane node at separatix 101
+                    subset_omp = get_grid_subset_with_index(grid_ggd, 11)
+                    subset_ompsep = get_grid_subset_with_index(grid_ggd, 101)
                 end
                 if !isnothing(jxi)
-                    subset_imp = get_grid_subset_with_index(grid_ggd, 12)     # inner midplane 12
-                    subset_impsep = get_grid_subset_with_index(grid_ggd, 102) # inner midplane node at separatix 102
+                    subset_imp = get_grid_subset_with_index(grid_ggd, 12)
+                    subset_impsep = get_grid_subset_with_index(grid_ggd, 102)
                 end
-                subset_corebnd = get_grid_subset_with_index(grid_ggd, 15)    # core inner boundary index 15
-                subset_separatix = get_grid_subset_with_index(grid_ggd, 16)  # separatix index 16
-                subset_otsep = get_grid_subset_with_index(grid_ggd, 103)     # outer target separatix meeting point 103
-                subset_itsep = get_grid_subset_with_index(grid_ggd, 104)     # inner target separatix meeting point 104
+                subset_corebnd = get_grid_subset_with_index(grid_ggd, 15)
+                subset_separatix = get_grid_subset_with_index(grid_ggd, 16)
+                subset_otsep = get_grid_subset_with_index(grid_ggd, 103)
+                subset_itsep = get_grid_subset_with_index(grid_ggd, 104)
             end
-            
-            subset_otarget = get_grid_subset_with_index(grid_ggd, 13)     # outer_target index 13
-            subset_itarget = get_grid_subset_with_index(grid_ggd, 14)     # inner_target index 14
+
+            subset_otarget = get_grid_subset_with_index(grid_ggd, 13)
+            subset_itarget = get_grid_subset_with_index(grid_ggd, 14)
 
             # Resizing objects to hold cell geometry data
             # Should be fewer than this many points, but this way we won't under-fill
@@ -962,19 +1008,19 @@ function solps2imas(b2gmtry, b2output, gsdesc, b2mn=nothing; load_bb=false)
             resize!(o3.object, ncell)  # Volumes
 
             # Initialize geometry for 0D objects(nodes), nodes for 1D objects(edges)
-            for i = 1:(ncell * 4)
+            for i ∈ 1:(ncell*4)
                 nodes[i].geometry = [0.0, 0.0]
                 edges[i].nodes = [0, 0]
                 resize!(edges[i].boundary, 2)
-                for bnd in edges[i].boundary
+                for bnd ∈ edges[i].boundary
                     bnd.neighbours = Int64[]
                 end
             end
             # Initialize nodes and boundaries for cells
-            for i = 1:(ncell)
+            for i ∈ 1:(ncell)
                 cells[i].nodes = [0, 0, 0, 0]
                 resize!(cells[i].boundary, 4)
-                for bnd in cells[i].boundary
+                for bnd ∈ cells[i].boundary
                     bnd.neighbours = Int64[]
                 end
             end
@@ -982,19 +1028,25 @@ function solps2imas(b2gmtry, b2output, gsdesc, b2mn=nothing; load_bb=false)
             j = 1
             edge_ind = 1
             # Setting up space with nodes, edges and cells
-            for iy = 1:ny
-                for ix = 1:nx
+            for iy ∈ 1:ny
+                for ix ∈ 1:nx
                     ic::Int = (iy - 1) * nx + ix
-                    # Adding node positions data to grid_ggd[grid_number].space[space_number].objects_per_dimension[0].object[:].geometry
-                    # Adding cell corners data to grid_ggd[grid_number].space[space_number].objects_per_dimension[2].object[:].nodes[1:4]
-                    for icorner = 1:4
-                        # Have to search to see if the node is already added and then record its index
-                        # If not already listed, then list it under new index and record that
-                        # Note that time index has been fixed to 1 here. Only handling fixed grid geometry
-                        # through the run cases.
-                        i_existing = search_points(ids, crx[1, icorner, iy, ix], cry[1, icorner, iy, ix])[1]
+                    # Adding node positions and cell corners data
+                    for icorner ∈ 1:4
+                        # Have to search to see if the node is already added and then
+                        # record its index
+                        # If not already listed, then list it under new index and
+                        # record that
+                        # Note that time index has been fixed to 1 here. Only handling
+                        # fixed grid geometry through the run cases.
+                        i_existing = search_points(
+                            ids,
+                            crx[1, icorner, iy, ix],
+                            cry[1, icorner, iy, ix],
+                        )[1]
                         if i_existing == 0
-                            nodes[j].geometry = [crx[1, icorner, iy, ix], cry[1, icorner, iy, ix]]
+                            nodes[j].geometry =
+                                [crx[1, icorner, iy, ix], cry[1, icorner, iy, ix]]
                             cells[ic].nodes[icorner] = j
                             add_subset_element!(subset_nodes, sn, 0, j)
                             if cuts_found && xpoints_nodes[it][1] == nodes[j].geometry
@@ -1005,21 +1057,36 @@ function solps2imas(b2gmtry, b2output, gsdesc, b2mn=nothing; load_bb=false)
                             cells[ic].nodes[icorner] = i_existing[1]
                         end
                     end
-                    # Adding edges (faced with toroidal elongation) to grid_ggd[grid_number].space[space_number].objects_per_dimension[1].object[:].nodes[1:2]
-                    # Adding same edges as boundary to  grid_ggd[grid_number].space[space_number].objects_per_dimension[2].object[:].boundary[1:4]
-                    for (boundary_ind, edge_pair) in chosen_edge_order
-                        edge_nodes = [cells[ic].nodes[icorner] for icorner in edge_pair]
+                    # Adding edges (faces with toroidal elongation)
+                    # Adding same edges as boundary to cells
+                    for (boundary_ind, edge_pair) ∈ chosen_edge_order
+                        edge_nodes = [cells[ic].nodes[icorner] for icorner ∈ edge_pair]
                         existing_edge_ind = search_edges(edges, edge_nodes)
                         if existing_edge_ind == 0
                             edges[edge_ind].nodes = edge_nodes
-                            for (ii, edge_bnd) in enumerate(edges[edge_ind].boundary)
+                            for (ii, edge_bnd) ∈ enumerate(edges[edge_ind].boundary)
                                 edge_bnd.index = edge_nodes[ii]
                             end
-                            edges[edge_ind].measure = distance_between_nodes(nodes, edge_nodes)
+                            edges[edge_ind].measure =
+                                distance_between_nodes(nodes, edge_nodes)
                             cells[ic].boundary[boundary_ind].index = edge_ind
                             add_subset_element!(subset_faces, sn, 1, edge_ind)
-                            add_subset_element!(subset_xedges, sn, 1, edge_ind, is_x_aligned; boundary_ind)
-                            add_subset_element!(subset_yedges, sn, 1, edge_ind, is_y_aligned; boundary_ind)
+                            add_subset_element!(
+                                subset_xedges,
+                                sn,
+                                1,
+                                edge_ind,
+                                is_x_aligned;
+                                boundary_ind,
+                            )
+                            add_subset_element!(
+                                subset_yedges,
+                                sn,
+                                1,
+                                edge_ind,
+                                is_y_aligned;
+                                boundary_ind,
+                            )
                             edge_ind += 1
                         else
                             cells[ic].boundary[boundary_ind].index = existing_edge_ind
@@ -1027,11 +1094,37 @@ function solps2imas(b2gmtry, b2output, gsdesc, b2mn=nothing; load_bb=false)
                     end
                     add_subset_element!(subset_cells, sn, 2, ic)
                     if cuts_found
-                        # add_subset_element!(subset, sn, dim, index, ix, iy, in_subset=(x...) -> true; kwargs...)
-                        add_subset_element!(subset_core, sn, 2, ic, in_core; ix, iy, cuts...)
+                        add_subset_element!(
+                            subset_core,
+                            sn,
+                            2,
+                            ic,
+                            in_core;
+                            ix,
+                            iy,
+                            cuts...,
+                        )
                         add_subset_element!(subset_sol, sn, 2, ic, in_sol; iy, cuts...)
-                        add_subset_element!(subset_idr, sn, 2, ic, in_idr; ix, iy, cuts...)
-                        add_subset_element!(subset_odr, sn, 2, ic, in_odr; ix, iy, cuts...)
+                        add_subset_element!(
+                            subset_idr,
+                            sn,
+                            2,
+                            ic,
+                            in_idr;
+                            ix,
+                            iy,
+                            cuts...,
+                        )
+                        add_subset_element!(
+                            subset_odr,
+                            sn,
+                            2,
+                            ic,
+                            in_odr;
+                            ix,
+                            iy,
+                            cuts...,
+                        )
                     end
                 end
             end
@@ -1039,24 +1132,93 @@ function solps2imas(b2gmtry, b2output, gsdesc, b2mn=nothing; load_bb=false)
                 # Add boundaries
                 attach_neightbours(cells, edges, gmtry, it)
                 # Adding edges to subsets
-                for iy = 1:ny
-                    for ix = 1:nx
-                        for boundary_ind = 1:4
-                            edge_ind = cells[xytoc(ix, iy; nx)].boundary[boundary_ind].index
-                            add_subset_element!(subset_corecut, sn, 1, edge_ind, is_core_cut; ix, iy, cells, nx, boundary_ind, cuts...)
-                            # add_subset_element!(subset_pfrcut, sn, 1, edge_ind, is_pfr_cut; ix, iy, cells, nx, boundary_ind, cuts...)
+                for iy ∈ 1:ny
+                    for ix ∈ 1:nx
+                        for boundary_ind ∈ 1:4
+                            edge_ind =
+                                cells[xytoc(ix, iy; nx)].boundary[boundary_ind].index
+                            add_subset_element!(
+                                subset_corecut,
+                                sn,
+                                1,
+                                edge_ind,
+                                is_core_cut;
+                                ix,
+                                iy,
+                                cells,
+                                nx,
+                                boundary_ind,
+                                cuts...,
+                            )
                             if !isnothing(jxa)
-                                add_subset_element!(subset_omp, sn, 1, edge_ind, is_outer_midplane; ix, iy, jxa, boundary_ind, cuts...)
+                                add_subset_element!(
+                                    subset_omp,
+                                    sn,
+                                    1,
+                                    edge_ind,
+                                    is_outer_midplane;
+                                    ix,
+                                    iy,
+                                    jxa,
+                                    boundary_ind,
+                                    cuts...,
+                                )
                             end
                             if !isnothing(jxi)
-                                add_subset_element!(subset_imp, sn, 1, edge_ind, is_inner_midplane; ix, iy, jxi, boundary_ind, cuts...)
+                                add_subset_element!(
+                                    subset_imp,
+                                    sn,
+                                    1,
+                                    edge_ind,
+                                    is_inner_midplane;
+                                    ix,
+                                    iy,
+                                    jxi,
+                                    boundary_ind,
+                                    cuts...,
+                                )
                             end
-                            add_subset_element!(subset_othroat, sn, 1, edge_ind, is_outer_throat; ix, iy, boundary_ind, cuts...)
-                            add_subset_element!(subset_ithroat, sn, 1, edge_ind, is_inner_throat; ix, iy, boundary_ind, cuts...)
-                            add_subset_element!(subset_otarget, sn, 1, edge_ind, is_outer_target; ix, nx, boundary_ind)
-                            add_subset_element!(subset_itarget, sn, 1, edge_ind, is_inner_target; ix, boundary_ind)
-                            # add_subset_element!(subset_corebnd, sn, 1, edge_ind, is_core_boundary; ix, iy, boundary_ind, cuts...)
-                            # add_subset_element!(subset_separatix, sn, 1, edge_ind, is_separatix; iy, boundary_ind, cuts...)
+                            add_subset_element!(
+                                subset_othroat,
+                                sn,
+                                1,
+                                edge_ind,
+                                is_outer_throat;
+                                ix,
+                                iy,
+                                boundary_ind,
+                                cuts...,
+                            )
+                            add_subset_element!(
+                                subset_ithroat,
+                                sn,
+                                1,
+                                edge_ind,
+                                is_inner_throat;
+                                ix,
+                                iy,
+                                boundary_ind,
+                                cuts...,
+                            )
+                            add_subset_element!(
+                                subset_otarget,
+                                sn,
+                                1,
+                                edge_ind,
+                                is_outer_target;
+                                ix,
+                                nx,
+                                boundary_ind,
+                            )
+                            add_subset_element!(
+                                subset_itarget,
+                                sn,
+                                1,
+                                edge_ind,
+                                is_inner_target;
+                                ix,
+                                boundary_ind,
+                            )
                         end
                     end
                 end
@@ -1064,20 +1226,46 @@ function solps2imas(b2gmtry, b2output, gsdesc, b2mn=nothing; load_bb=false)
                 sol_boundary_elements = get_subset_boundary(space, subset_sol)
                 idr_boundary_elements = get_subset_boundary(space, subset_idr)
                 odr_boundary_elements = get_subset_boundary(space, subset_odr)
-                subset_pfrcut.element = subset_do(intersect, idr_boundary_elements, odr_boundary_elements)
-                subset_corebnd.element = subset_do(setdiff, core_boundary_elements, sol_boundary_elements)
+                subset_pfrcut.element =
+                    subset_do(intersect, idr_boundary_elements, odr_boundary_elements)
+                subset_corebnd.element =
+                    subset_do(setdiff, core_boundary_elements, sol_boundary_elements)
                 subset_separatix.element = subset_do(intersect, sol_boundary_elements,
-                                                                   subset_do(union, core_boundary_elements,
-                                                                                       odr_boundary_elements,
-                                                                                       idr_boundary_elements))
+                    subset_do(union, core_boundary_elements,
+                        odr_boundary_elements,
+                        idr_boundary_elements))
                 if !isnothing(jxa)
-                    subset_ompsep.element = subset_do(intersect, subset_separatix.element, subset_omp.element; space, use_nodes=true)
+                    subset_ompsep.element = subset_do(
+                        intersect,
+                        subset_separatix.element,
+                        subset_omp.element;
+                        space,
+                        use_nodes=true,
+                    )
                 end
                 if !isnothing(jxi)
-                    subset_impsep.element = subset_do(intersect, subset_separatix.element, subset_imp.element; space, use_nodes=true)
+                    subset_impsep.element = subset_do(
+                        intersect,
+                        subset_separatix.element,
+                        subset_imp.element;
+                        space,
+                        use_nodes=true,
+                    )
                 end
-                subset_otsep.element = subset_do(intersect, subset_separatix.element, subset_otarget.element; space, use_nodes=true)
-                subset_itsep.element = subset_do(intersect, subset_separatix.element, subset_itarget.element; space, use_nodes=true)
+                subset_otsep.element = subset_do(
+                    intersect,
+                    subset_separatix.element,
+                    subset_otarget.element;
+                    space,
+                    use_nodes=true,
+                )
+                subset_itsep.element = subset_do(
+                    intersect,
+                    subset_separatix.element,
+                    subset_itarget.element;
+                    space,
+                    use_nodes=true,
+                )
             end
         end  # End of setting up space
     end
@@ -1086,15 +1274,15 @@ function solps2imas(b2gmtry, b2output, gsdesc, b2mn=nothing; load_bb=false)
     b2 = read_b2_output(b2output)
     # Add grid_ggd array equal to number of time steps
     resize!(ids.edge_profiles.ggd, b2["dim"]["time"])
-    for it in 1:b2["dim"]["time"]
+    for it ∈ 1:b2["dim"]["time"]
         ggd = ids.edge_profiles.ggd[it]
         ggd.time = Float64.(b2["data"]["timesa"][it])
-        for (key, data) in b2["data"]
+        for (key, data) ∈ b2["data"]
             obj = val_obj(ggd, key, gsdesc["identifier"]["index"])
             if !isnothing(obj)
                 resize!(obj, ncell)
-                for iy = 1:ny
-                    for ix = 1:nx
+                for iy ∈ 1:ny
+                    for ix ∈ 1:nx
                         ic::Int = (iy - 1) * nx + ix
                         obj[ic] = data[it, iy, ix]
                     end
@@ -1111,7 +1299,7 @@ function solps2imas(b2gmtry, b2output, gsdesc, b2mn=nothing; load_bb=false)
             resize!(ids.equilibrium.time_slice, gmtry["dim"]["time"])
             ids.equilibrium.time = gmtry["data"]["timesa"]
         end
-        for it in 1:gmtry["dim"]["time"]
+        for it ∈ 1:gmtry["dim"]["time"]
             # Note
             # Ideally, equilibrium keeps separate grids_ggd object for each time step
             # But since we have already created them in edge_profiles.grid_ggd, we
@@ -1128,13 +1316,14 @@ function solps2imas(b2gmtry, b2output, gsdesc, b2mn=nothing; load_bb=false)
             b_r = time_slice.ggd[1].b_field_r[1]
             b_z = time_slice.ggd[1].b_field_z[1]
 
-            b_z.grid_index = b_r.grid_index = b_t.grid_index = gsdesc["identifier"]["index"]
+            b_z.grid_index =
+                b_r.grid_index = b_t.grid_index = gsdesc["identifier"]["index"]
             b_z.grid_subset_index = b_r.grid_subset_index = b_t.grid_subset_index = 5
             resize!(b_z.values, ncell)
             resize!(b_r.values, ncell)
             resize!(b_t.values, ncell)
-            for iy = 1:ny
-                for ix = 1:nx
+            for iy ∈ 1:ny
+                for ix ∈ 1:nx
                     ic::Int = (iy - 1) * nx + ix
                     b_z.values[ic] = bb[it, 1, iy, ix]
                     b_r.values[ic] = bb[it, 2, iy, ix]
diff --git a/test/runtests.jl b/test/runtests.jl
index 46278c6..c4e0f12 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,4 +1,4 @@
-import SOLPS2IMAS
+using SOLPS2IMAS: SOLPS2IMAS
 using Test
 using YAML: load_file as YAML_load_file
 
@@ -12,18 +12,27 @@ function test_ind_conversion()
     nx = 92
     ny = 38
     success = true
-    for iy = 1:ny
-        for ix = 1:nx
+    for iy ∈ 1:ny
+        for ix ∈ 1:nx
             cxy = SOLPS2IMAS.ctoxy(SOLPS2IMAS.xytoc(ix, iy; nx=nx); nx=nx)
             if cxy != (ix, iy)
                 ic = SOLPS2IMAS.xytoc(ix, iy; nx=nx)
-                println("ic: ", ic, ", (ix, iy): (", ix, ", ", iy, "), converted_ix,iy: ", cxy)
+                println(
+                    "ic: ",
+                    ic,
+                    ", (ix, iy): (",
+                    ix,
+                    ", ",
+                    iy,
+                    "), converted_ix,iy: ",
+                    cxy,
+                )
                 success = false
             end
         end
     end
 
-    for ic = 1:nx*ny
+    for ic ∈ 1:nx*ny
         cic = SOLPS2IMAS.xytoc(SOLPS2IMAS.ctoxy(ic; nx=nx)...; nx=nx)
         if cic != ic
             ix, iy = SOLPS2IMAS.ctoxy(ic; nx=nx)
@@ -69,7 +78,7 @@ function test_read_b2_output()
     @assert(size(contents["data"]["fl3dr"]) == (nt, nybr))
     @assert(size(contents["data"]["fc3dr"]) == (nt, nybr))
     @assert(size(contents["data"]["fna3da"]) == (nt, ns, nybr))
-    @assert(size(contents["data"]["tmhacore"]) == (nt, ))
+    @assert(size(contents["data"]["tmhacore"]) == (nt,))
     @assert(size(contents["data"]["tmhasol"]) == (nt,))
     @assert(size(contents["data"]["tmhadiv"]) == (nt,))
     return true
@@ -88,29 +97,37 @@ function test_solps2imas()
     it = 3
     iy = 4
     ix = 5
-    @assert(b2t["data"]["ne2d"][3, iy, ix] == dd.edge_profiles.ggd[3].electrons.density[5].values[(iy - 1) * nx + ix])
+    @assert(
+        b2t["data"]["ne2d"][3, iy, ix] ==
+        dd.edge_profiles.ggd[3].electrons.density[5].values[(iy-1)*nx+ix]
+    )
     # Checking if subsets obtained from set operations are same
     # as using a brute force definition which is too dependent
     # on the correct ordering of nodes in SOLPS data files.
     cut_keys = ["leftcut", "rightcut", "bottomcut", "topcut"]
     gmtry = SOLPS2IMAS.read_b2_output(b2gmtry)
     ny = gmtry["dim"]["ny"]
-    cuts = Dict([(Symbol(key), gmtry["data"][key][1]) for key in cut_keys])
-    subset_pfrcut = SOLPS2IMAS.get_grid_subset_with_index(dd.edge_profiles.grid_ggd[1], 8)
-    subset_corebnd = SOLPS2IMAS.get_grid_subset_with_index(dd.edge_profiles.grid_ggd[1], 15)
-    subset_separatix= SOLPS2IMAS.get_grid_subset_with_index(dd.edge_profiles.grid_ggd[1], 16)
+    cuts = Dict([(Symbol(key), gmtry["data"][key][1]) for key ∈ cut_keys])
+    subset_pfrcut =
+        SOLPS2IMAS.get_grid_subset_with_index(dd.edge_profiles.grid_ggd[1], 8)
+    subset_corebnd =
+        SOLPS2IMAS.get_grid_subset_with_index(dd.edge_profiles.grid_ggd[1], 15)
+    subset_separatix =
+        SOLPS2IMAS.get_grid_subset_with_index(dd.edge_profiles.grid_ggd[1], 16)
     cells = dd.edge_profiles.grid_ggd[1].space[1].objects_per_dimension[3].object
-    subset_pfrcut_element_list = [ele.object[1].index for ele in subset_pfrcut.element]
-    subset_corebnd_element_list = [ele.object[1].index for ele in subset_corebnd.element]
-    subset_separatix_element_list = [ele.object[1].index for ele in subset_separatix.element]
+    subset_pfrcut_element_list = [ele.object[1].index for ele ∈ subset_pfrcut.element]
+    subset_corebnd_element_list = [ele.object[1].index for ele ∈ subset_corebnd.element]
+    subset_separatix_element_list =
+        [ele.object[1].index for ele ∈ subset_separatix.element]
     brute_force_pfrcut_list = []
     brute_force_corebnd_list = []
     brute_force_separatix_list = []
-    for iy = 1:ny
-        for ix = 1:nx
-            for boundary_ind = 1:4
-                edge_ind = cells[SOLPS2IMAS.xytoc(ix, iy; nx)].boundary[boundary_ind].index
-                if SOLPS2IMAS.is_pfr_cut(;ix, iy, cells, nx, boundary_ind, cuts...)
+    for iy ∈ 1:ny
+        for ix ∈ 1:nx
+            for boundary_ind ∈ 1:4
+                edge_ind =
+                    cells[SOLPS2IMAS.xytoc(ix, iy; nx)].boundary[boundary_ind].index
+                if SOLPS2IMAS.is_pfr_cut(; ix, iy, cells, nx, boundary_ind, cuts...)
                     append!(brute_force_pfrcut_list, edge_ind)
                 elseif SOLPS2IMAS.is_core_boundary(; ix, iy, boundary_ind, cuts...)
                     append!(brute_force_corebnd_list, edge_ind)
@@ -126,7 +143,6 @@ function test_solps2imas()
     return true
 end
 
-
 @testset "omasstuff" begin
     @test SOLPS2IMAS.try_omas() === nothing
     @test test_ind_conversion()