diff --git a/tools/tokamak_grids/gridue_to_bout b/tools/tokamak_grids/gridue_to_bout
new file mode 100755
index 0000000000..17f07b7821
--- /dev/null
+++ b/tools/tokamak_grids/gridue_to_bout
@@ -0,0 +1,791 @@
+#!/usr/bin/env python3
+#
+# Convert UEDGE grid files (gridue) to BOUT++ grids
+#
+# Parts adapted from INGRID https://github.com/LLNL/INGRID
+#  Copyright (c) 2020, Lawrence Livermore National Security, LLC
+# Parts adapted from Hypnotoad https://github.com/boutproject/hypnotoad/
+#  Copyright 2019 J.T. Omotani
+#
+
+import numpy as np
+import matplotlib
+import matplotlib.pyplot as plt
+
+from scipy import linalg
+
+
+def _importBody(gridue_settings, f):
+    next(f)
+    BodyItems = ["rm", "zm", "psi", "br", "bz", "bpol", "bphi", "b"]
+    Str = {i: [] for i in BodyItems}
+    k = iter(Str.keys())
+    Key = next(k)
+
+    for line in f:
+        if line == "iogridue\n":
+            continue
+        if line == "\n":
+            try:
+                Key = next(k)
+            except:
+                continue
+        else:
+            Str[Key].append(line)
+    f.close()
+    nx = gridue_settings["nxm"] + 2
+    ny = gridue_settings["nym"] + 2
+
+    for k, v in Str.items():
+        L = ("".join(v).replace("\n", "").replace("D", "e")).split()
+        _l = iter(L)
+        vv = next(_l)
+
+        data_ = np.zeros((nx, ny, 5))
+        for n in range(5):
+            for j in range(ny):
+                for i in range(nx):
+
+                    data_[i][j][n] = float(vv)
+
+                    try:
+                        vv = next(_l)
+                    except:
+                        continue
+                    gridue_settings[k] = data_
+    return gridue_settings
+
+
+def _importSN(values, f):
+    """
+    Import a single null file
+    """
+    HeaderItems = ["nxm", "nym", "ixpt1", "ixpt2", "iyseparatrix1"]
+    gridue_settings = dict(zip(HeaderItems, values))
+
+    # Add indices that are present in DN but not SN
+    gridue_settings.update(
+        {
+            "iyseparatrix2": gridue_settings["nym"] + 10,  # Outside grid
+            "ix_cut1": gridue_settings["ixpt1"],
+            "ix_cut2": gridue_settings["nxm"] // 2,
+            "ix_inner": gridue_settings["nxm"] // 2,
+            "ix_cut3": gridue_settings["nxm"] // 2,
+            "ix_cut4": gridue_settings["ixpt2"],
+        }
+    )
+
+    return _importBody(gridue_settings, f)
+
+
+def _importDN(values, f):
+
+    gridue_settings = dict(zip(["nxm", "nym"], values))
+
+    header_rows = [
+        ["iyseparatrix1", "iyseparatrix2"],
+        ["ix_plate1", "ix_cut1", "_FILLER_", "ix_cut2", "ix_plate2"],
+        ["iyseparatrix3", "iyseparatrix4"],
+        ["ix_plate3", "ix_cut3", "_FILLER_", "ix_cut4", "ix_plate4"],
+    ]
+
+    for row in header_rows:
+        values = [int(x) for x in next(f).split()]
+        if len(values) != len(row):
+            raise ValueError(
+                "Expected row with {} integers, found {}".format(len(row), len(values))
+            )
+        gridue_settings.update(zip(row, values))
+
+    gridue_settings.update({"ix_inner": gridue_settings["ix_plate2"]})
+
+    return _importBody(gridue_settings, f)
+
+
+def importGridue(fname: str = "gridue") -> dict:
+    """
+    Import UEDGE grid file as dictionary.
+
+    Parameters
+    ----------
+    fname : str, optional
+        Path/file name to gridue formatted file.
+
+    Returns
+    -------
+    A dict containing header and body information from the gridue file.
+
+    """
+    f = open(fname, mode="r")
+    values = [int(x) for x in next(f).split()]
+
+    if len(values) == 5:
+        return _importSN(values, f)
+    if len(values) == 2:
+        return _importDN(values, f)
+
+    raise ValueError("Unrecognised gridue format")
+
+
+def plot(GridueParams: dict, edgecolor="black", ax: object = None, show=True):
+    """
+    Plot UEDGE grid from 'dict' obtained from method 'ImportGridue'
+
+    Parameters
+    ----------
+    GridueParams : dict
+        Gridue header and body information as a dictionary.
+
+    edgecolor : str, optional
+        Color of grid.
+
+    ax : object, optional
+        Matplotlib axes to plot on.
+
+    """
+    r = GridueParams["rm"]
+    z = GridueParams["zm"]
+    Nx = len(r)
+    Ny = len(r[0])
+    patches = []
+    plt.figure(figsize=(6, 10))
+    if ax is None:
+        ax = plt.gca()
+    idx = [np.array([1, 2, 4, 3, 1])]
+    for i in range(Nx):
+        for j in range(Ny):
+            p = matplotlib.patches.Polygon(
+                np.concatenate((r[i][j][idx], z[i][j][idx])).reshape(2, 5).T,
+                fill=False,
+                closed=True,
+                edgecolor=edgecolor,
+            )
+            ax.add_patch(p)  # draw the contours
+    ax.set_aspect("equal", adjustable="box")
+    ax.set_xlabel("R")
+    ax.set_ylabel("Z")
+    ax.set_ylim(z.min(), z.max())
+    ax.set_xlim(r.min(), r.max())
+    if show:
+        plt.show()
+    return ax
+
+
+def calcHy(g: dict, dy: float):
+    """
+    Calculate poloidal arc length metric from gridue dictionary
+    """
+    rm = g["rm"]
+    zm = g["zm"]
+
+    hy = np.zeros((nx, ny))
+    for i in range(nx):
+        for j in range(ny):
+            r = rm[j, i, :]
+            z = zm[j, i, :]
+            # Find intersection with (1) -- (3)
+            R1 = 0.5 * (r[1] + r[3])
+            Z1 = 0.5 * (z[1] + z[3])
+            # Find intersection with (2) -- (4)
+            R2 = 0.5 * (r[2] + r[4])
+            Z2 = 0.5 * (z[2] + z[4])
+
+            # Distance from one side to the other
+            dl = np.sqrt((R2 - R1) ** 2 + (Z2 - Z1) ** 2)
+            hy[i, j] = dl / dy
+    return hy
+
+
+def calcGridAngle(g: dict):
+    Brxy = g["br"][:, :, 0].T
+    Bzxy = g["bz"][:, :, 0].T
+    Bpxy = g["bpol"][:, :, 0].T
+
+    psi = g["psi"]
+    rm = g["rm"]
+    zm = g["zm"]
+
+    delta_y = [Brxy / Bpxy, Bzxy / Bpxy]  # Unit vector along e_y
+
+    R_xlow = 0.5 * (rm[:, :, 1] + rm[:, :, 2]).T
+    Z_xlow = 0.5 * (zm[:, :, 1] + zm[:, :, 2]).T
+    R_xhigh = 0.5 * (rm[:, :, 3] + rm[:, :, 4]).T
+    Z_xhigh = 0.5 * (zm[:, :, 3] + zm[:, :, 4]).T
+    dR = R_xhigh - R_xlow
+    dZ = Z_xhigh - Z_xlow
+    dl = np.sqrt(dR**2 + dZ**2)
+    delta_x = [dR / dl, dZ / dl]  # unit vector along e_x
+
+    # Calculate angle between x and y unit vectors.
+    # sin(beta) = cos(pi/2 - beta) = e_x_hat.e_y_hat = delta_x.delta_y
+    sinBeta = delta_x[0] * delta_y[0] + delta_x[1] * delta_y[1]
+
+    # Rotate delta_y by 90 degrees anticlockwise to get unit vector in psi direction
+    delta_psi = [-delta_y[1], delta_y[0]]
+
+    # cosBeta = delta_x.delta_psi
+    cosBeta = delta_x[0] * delta_psi[0] + delta_x[1] * delta_psi[1]
+
+    return sinBeta, cosBeta
+
+
+def calcRZderivs(R, Z, var):
+    """
+    Calculate derivatives of var w.r.t R and Z, using 5 points
+
+    Inputs R, Z and var should be 3D [poloidal, radial, 5]
+    as stored in gridue files.
+
+    Returns (dvar/dR, dvar/dZ) as 2D arrays [poloidal, radial]
+    at the location of point index 0
+    """
+    npol, nr, _ = R.shape
+    dR = np.zeros((npol, nr))
+    dZ = np.zeros((npol, nr))
+    for i in range(npol):
+        for j in range(nr):
+            A = np.zeros((5, 5))
+            for p in range(5):
+                # Constraint on derivatives for this point
+                A[p, 0] = 1.0  # Constant
+                A[p, 1] = R[i, j, p] - R[i, j, 0]  # dR
+                A[p, 2] = Z[i, j, p] - Z[i, j, 0]  # dZ
+                A[p, 3] = (R[i, j, p] - R[i, j, 0]) ** 2  # dR^2
+                A[p, 4] = (Z[i, j, p] - Z[i, j, 0]) ** 2  # dR^2
+            # Values of the function being fitted
+            b = var[i, j, :].squeeze()
+            # Invert using LU for stability
+            lu, piv = linalg.lu_factor(A)
+            x = linalg.lu_solve((lu, piv), b)
+            dR[i, j] = x[1]
+            dZ[i, j] = x[2]
+    return dR, dZ  # dvar/dR, dvar/dZ at cell center
+
+
+def calcMetric(grd: dict, bpsign, verbose=False, ignore_checks=False):
+    """
+    Calculate metric tensor given BOUT++ cell centers
+    """
+    Rxy = grd["Rxy"]
+    Zxy = grd["Zxy"]
+    Brxy = grd["Brxy"]
+    Bzxy = grd["Bzxy"]
+    Btxy = grd["Btxy"]
+    Bpxy = grd["Bpxy"]
+    Bxy = grd["Bxy"]
+    hy = grd["hy"]
+    cosBeta = grd["cosBeta"]
+    tanBeta = grd["tanBeta"]
+    curl_bOverB_Rhat = grd["curl_bOverB_Rhat"]
+    curl_bOverB_Zhat = grd["curl_bOverB_Zhat"]
+    curl_bOverB_zetahat = grd["curl_bOverB_zetahat"]
+
+    if verbose:
+        for var, name in [
+            (cosBeta, "cos(beta)"),
+            (tanBeta, "tan(beta)"),
+        ]:
+            print(
+                "{} min {}, mean {}, max {}".format(
+                    name, np.amin(var), np.mean(var), np.amax(var)
+                )
+            )
+
+    dphidy = hy * Btxy / (Bpxy * Rxy)
+
+    I = np.zeros(Rxy.shape)
+
+    g11 = (Rxy * Bpxy) ** 2
+    g22 = 1.0 / (hy * cosBeta) ** 2
+    g33 = (
+        1.0 / Rxy**2
+        + (Rxy * Bpxy * I) ** 2
+        + (dphidy / (hy * cosBeta)) ** 2
+        + 2.0 * Rxy * Bpxy * I * dphidy * tanBeta / hy
+    )
+    g12 = Rxy * np.abs(Bpxy) * tanBeta / hy
+    g13 = -Rxy * Bpxy * dphidy * tanBeta / hy - I * (Rxy * Bpxy) ** 2
+    g23 = -bpsign * dphidy / (hy * cosBeta) ** 2 - Rxy * np.abs(Bpxy) * I * tanBeta / hy
+
+    J = hy / Bpxy
+
+    g_11 = 1.0 / (Rxy * Bpxy * cosBeta) ** 2 + (I * Rxy) ** 2
+    g_22 = hy**2 + (dphidy * Rxy) ** 2
+    g_33 = Rxy**2
+    g_12 = bpsign * I * dphidy * Rxy**2 - hy * tanBeta / (Rxy * np.abs(Bpxy))
+    g_13 = I * Rxy**2
+    g_23 = bpsign * dphidy * Rxy**2
+
+    Jcheck = (
+        bpsign
+        * 1.0
+        / np.sqrt(
+            g11 * g22 * g33
+            + 2.0 * g12 * g13 * g23
+            - g11 * g23**2
+            - g22 * g13**2
+            - g33 * g12**2
+        )
+    )
+
+    rel_error = (J - Jcheck) / J
+
+    if verbose:
+        for var, name in [
+            (hy, "hy"),
+            (J, "J"),
+            (Jcheck, "Jcheck"),
+            (J - Jcheck, "J - Jcheck"),
+            (rel_error, "(J - Jcheck)/J"),
+        ]:
+            print(
+                "{} min {}, mean {}, max {}".format(
+                    name, np.amin(var), np.mean(var), np.amax(var)
+                )
+            )
+
+    if np.max(np.abs(rel_error)) > 1e-6:
+        if ignore_checks:
+            print("WARNING: Relative error in Jacobian too large.")
+        else:
+            raise ValueError("Relative error in Jacobian too large.")
+
+    # We want to output contravariant components of Curl(b/B) in the
+    # locally field-aligned coordinate system.
+    # The contravariant components of an arbitrary vector A are
+    # A^x = A.Grad(x)
+    # A^y = A.Grad(y)
+    # A^z = A.Grad(z)
+
+    # Grad in cylindrical coordinates is
+    # Grad(f) = df/dR Rhat + 1/R df/dzeta zetahat + df/dZ Zhat
+    # https://en.wikipedia.org/wiki/Del_in_cylindrical_and_spherical_coordinates,
+
+    # x = psi - psi_min
+    # dpsi/dR = -R*BZ
+    # dpsi/dZ = R*BR
+    # => Grad(x) = (dpsi/dR, 0, dpsi/dZ).(Rhat, zetahat, Zhat)
+    # => Grad(x) = (-R BZ, 0, R BR).(Rhat, zetahat, Zhat)
+    curl_bOverB_x = -Rxy * Bzxy * curl_bOverB_Rhat + Rxy * Brxy * curl_bOverB_Zhat
+
+    # Grad(y) = (d_Z, 0, -d_R)/(hy*cosBeta)
+    #         = (BR*cosBeta-BZ*sinBeta, 0, BZ*cosBeta+BR*sinBeta)
+    #           /(Bp*hy*cosBeta)
+    #         = (BR-BZ*tanBeta, 0, BZ+BR*tanBeta)/(Bp*hy)
+    curl_bOverB_y = (
+        curl_bOverB_Rhat * (Brxy - Bzxy * tanBeta)
+        + curl_bOverB_Zhat * (Bzxy + Brxy * tanBeta)
+    ) / (Bpxy * hy)
+
+    # Grad(z) = Grad(zeta) - Bt*hy/(Bp*R)*Grad(y) - I*Grad(x)
+    # Grad(z) = (0, 1/R, 0) - Bt*hy/(Bp*R)*Grad(y) - I*Grad(x)
+    curl_bOverB_z = (
+        curl_bOverB_zetahat / Rxy
+        - Btxy * hy / (Bpxy * Rxy) * curl_bOverB_y
+        - I * curl_bOverB_x
+    )
+
+    bxcvx = Bxy / 2.0 * curl_bOverB_x
+    bxcvy = Bxy / 2.0 * curl_bOverB_y
+    bxcvz = Bxy / 2.0 * curl_bOverB_z
+
+    if verbose:
+        for var, name in [
+            (bxcvx, "bxcvx"),
+            (bxcvy, "bxcvy"),
+            (bxcvz, "bxcvz"),
+        ]:
+            print(
+                "{} min {}, mean {}, max {}".format(
+                    name, np.amin(var), np.mean(var), np.amax(var)
+                )
+            )
+
+    return {
+        "dphidy": dphidy,
+        # Metric tensor
+        "g11": g11,
+        "g22": g22,
+        "g33": g33,
+        "g12": g12,
+        "g13": g13,
+        "g23": g23,
+        # Inverse metric tensor
+        "g_11": g_11,
+        "g_22": g_22,
+        "g_33": g_33,
+        "g_12": g_12,
+        "g_13": g_13,
+        "g_23": g_23,
+        # Jacobian
+        "J": J,
+        # Integrated shear
+        "sinty": I,
+        # Curvature
+        "curl_bOverB_x": curl_bOverB_x,
+        "curl_bOverB_y": curl_bOverB_y,
+        "curl_bOverB_z": curl_bOverB_z,
+        "bxcvx": bxcvx,
+        "bxcvy": bxcvy,
+        "bxcvz": bxcvz,
+    }
+
+
+def calcRZCurvature(g: dict):
+    """
+    Calculate curvature in (R, Z, zeta)
+    Returns 2D arrays [radial, poloidal]
+    """
+    # Variables in UEDGE format [poloidal, radial, 5]
+    BR = g["br"]
+    BZ = g["bz"]
+    Bzeta = g["bphi"]
+    B2 = g["b"] ** 2
+    R = g["rm"]
+    Z = g["zm"]
+
+    dBzetadR, dBzetadZ = calcRZderivs(R, Z, Bzeta)
+    dBRdR, dBRdZ = calcRZderivs(R, Z, BR)
+    dBZdR, dBZdZ = calcRZderivs(R, Z, BZ)
+    dB2dR, dB2dZ = calcRZderivs(R, Z, B2)
+
+    # Select point at centre of cell
+    BR = BR[:, :, 0]
+    BZ = BZ[:, :, 0]
+    Bzeta = Bzeta[:, :, 0]
+    B2 = B2[:, :, 0]
+    R = R[:, :, 0]
+    Z = Z[:, :, 0]
+
+    # In cylindrical coords
+    # curl(A) = (1/R*d(AZ)/dzeta - d(Azeta)/dZ)  * Rhat
+    #           + 1/R*(d(R Azeta)/dR - d(AR)/dzeta) * Zhat
+    #           + (d(AR)/dZ - d(AZ)/dR) * zetahat
+    # Where AR, AZ and Azeta are the components on a basis of unit vectors,
+    # i.e. AR = A.Rhat; AZ = A.Zhat; Azeta = A.zetahat
+    # https://en.wikipedia.org/wiki/Del_in_cylindrical_and_spherical_coordinates,
+    #
+    # curl(b/B) = curl((BR/B2), (BZ/B2), (Bzeta/B2))
+    # curl(b/B)_Rhat = 1/R d(BZ/B2)/dzeta - d(Bzeta/B2)/dZ
+    #                = 1/(R*B2)*d(BZ)/dzeta - BZ/(R*B4)*d(B2)/dzeta
+    #                  - 1/B2*d(Bzeta)/dZ + Bzeta/B4*d(B2)/dZ
+    #                = -1/B2*d(Bzeta)/dZ + Bzeta/B4*d(B2)/dZ
+    # curl(b/B)_Zhat = 1/R * (d(R Bzeta/B2)/dR - d(BR/B2)/dzeta)
+    #                = Bzeta/(R*B2) + 1/B2*d(Bzeta)/dR - Bzeta/B4*d(B2)/dR
+    #                  - 1/(R*B2)*d(BR)/dzeta + BR/(R*B4)*d(B2)/dzeta
+    #                = Bzeta/(R*B2) + 1/B2*d(Bzeta)/dR - Bzeta/B4*d(B2)/dR
+    # curl(b/B)_zetahat = d(BR/B2)/dZ - d(BZ/B2)/dR
+    #                   = 1/B2*d(BR)/dZ - BR/B4*d(B2)/dZ
+    #                     - 1/B2*d(BZ)/dR + BZ/B4*d(B2)/dR
+    # remembering d/dzeta=0 for axisymmetric equilibrium
+
+    curl_bOverB_Rhat = -dBzetadZ / B2 + Bzeta / B2**2 * dB2dZ
+    curl_bOverB_Zhat = Bzeta / (R * B2) + dBzetadR / B2 - Bzeta / B2**2 * dB2dR
+    curl_bOverB_zetahat = (
+        dBRdZ / B2 - BR / B2**2 * dB2dZ - dBZdR / B2 + BZ / B2**2 * dB2dR
+    )
+
+    # Return as [radial, poloidal]
+    return curl_bOverB_Rhat.T, curl_bOverB_Zhat.T, curl_bOverB_zetahat.T
+
+
+if __name__ == "__main__":
+
+    import argparse
+
+    parser = argparse.ArgumentParser(
+        description="""Converts UEDGE grid files (gridue) into BOUT++ grids.
+Note that in most cases these grids are non-orthogonal."""
+    )
+    parser.add_argument("gridue_file", type=str)
+    parser.add_argument("-o", "--output", default="bout.grd.nc")
+    parser.add_argument("-p", "--plot", action="store_true", default=False)
+    parser.add_argument("-v", "--verbose", action="store_true", default=False)
+    parser.add_argument("-i", "--ignore-checks", action="store_true", default=False)
+    try:
+        import argcomplete
+
+        argcomplete.autocomplete(parser)
+    except ImportError:
+        pass
+
+    args = parser.parse_args()
+    gridue_file = args.gridue_file
+    output_filename = args.output
+    plotting = args.plot
+    verbose = args.verbose
+    ignore_checks = args.ignore_checks
+
+    g = importGridue(gridue_file)
+
+    if plotting:
+        plot(g)
+
+    psi = g["psi"]
+    rm = g["rm"]
+    zm = g["zm"]
+
+    Rxy = rm[:, :, 0].T
+    Zxy = zm[:, :, 0].T
+    Brxy = g["br"][:, :, 0].T
+    Bzxy = g["bz"][:, :, 0].T
+    Bpxy = g["bpol"][:, :, 0].T
+    Btxy = g["bphi"][:, :, 0].T
+    Bxy = g["b"][:, :, 0].T
+    psixy = psi[:, :, 0].T
+    nx, ny = Rxy.shape
+
+    # Ordering
+    # (1) -- (3)
+    #  |      |
+    #  |  (0) |  -> Radial, BOUT++ "x"
+    #  |      |
+    # (2) -- (4)
+
+    # calculate change in psi across cell -> dx
+
+    dx = np.zeros((nx, ny))
+    for i in range(nx):
+        for j in range(ny):
+            if i > 1 and i < nx - 2:
+                dx[i, j] = 0.5 * (psi[j, i + 1, 0] - psi[j, i - 1, 0])
+            else:
+                dx[i, j] = 0.5 * (
+                    psi[j, i, 3] + psi[j, i, 4] - psi[j, i, 1] - psi[j, i, 2]
+                )
+    # Note: psi and dx are not necessarily constant on flux surfaces
+
+    # Note: UEDGE grids have narrow cells on the radial
+    # boundaries. BOUT++ applies boundary conditions half-way between
+    # cells, and usually expects two boundary cells.
+
+    # Calculate direction of magnetic field Bp dot Grad(y)
+    Bp_dot_grady = Brxy[1, 1] * (Rxy[1, 1] - Rxy[1, 0]) + Bzxy[1, 1] * (
+        Zxy[1, 1] - Zxy[1, 0]
+    )
+
+    if Bp_dot_grady * Bpxy[1, 1] < 0.0:
+        raise ValueError("Bp_dot_grady and Bpxy have opposite signs")
+    bpsign = np.sign(Bpxy[1, 1])  # Sign of the poloidal magnetic field
+
+    # Choose angle dy across poloidal cell. This is somewhat arbitrary,
+    # but it is helpful if the y angle changes by 2pi for each poloidal transit of the core
+
+    dy = 2 * np.pi / ny
+
+    # Calculate hy, the arc length along the flux surface passing through
+    # the center of each cell.
+    hy = calcHy(g, dy)
+
+    # Calculate angle between x and y coordinates. sinBeta = 0, cosBeta = 1 for an orthogonal mesh
+    sinBeta, cosBeta = calcGridAngle(g)
+    tanBeta = sinBeta / cosBeta
+
+    # Calculate curvature
+    curl_bOverB_Rhat, curl_bOverB_Zhat, curl_bOverB_zetahat = calcRZCurvature(g)
+
+    # Collect 2D variables for output
+    grd = {
+        "Rxy": Rxy,
+        "Zxy": Zxy,
+        "psixy": psixy,
+        "dx": dx,
+        "dy": np.full((nx, ny), dy),
+        # Magnetic field components
+        "Brxy": Brxy,
+        "Bzxy": Bzxy,
+        "Bpxy": Bpxy,
+        "Btxy": Btxy,
+        "Bxy": Bxy,
+        # Poloidal arc length
+        "hy": hy,
+        "hthe": hy,
+        # Curvature
+        "curl_bOverB_Rhat": curl_bOverB_Rhat,
+        "curl_bOverB_Zhat": curl_bOverB_Zhat,
+        "curl_bOverB_zetahat": curl_bOverB_zetahat,
+        # Grid angles
+        "cosBeta": cosBeta,
+        "tanBeta": tanBeta,
+    }
+
+    # Remove Y (poloidal) boundary cells
+    for name in grd:
+        grd[name] = grd[name][:, 1:-1]
+
+    if g["ix_cut2"] != g["ix_cut3"]:
+        # Double null -> Remove upper Y guard cells
+        ny_inner = g["ix_inner"]
+        for name in grd:
+            var = grd[name]
+            nx, ny = var.shape
+            newvar = np.zeros((nx, ny - 2))
+            newvar[:, :ny_inner] = var[:, :ny_inner]
+            newvar[:, ny_inner:] = var[:, (ny_inner + 2) :]
+            grd[name] = newvar
+        g["ix_cut2"] = g["ix_cut2"] - 1
+        g["ix_cut3"] = g["ix_cut3"] - 3
+        g["ix_cut4"] = g["ix_cut4"] - 2
+
+    # Extrapolate X (radial) boundary cells
+    # Removing one cell, adding two on each X boundary
+    for name in grd:
+        var = grd[name]
+        nx, ny = var.shape
+        newvar = np.zeros((nx + 2, ny))
+        newvar[2:-2, :] = var[1:-1, :]
+        if name in ["Rxy", "Zxy", "psixy"]:
+            # Linear extrapolation
+            newvar[1, :] = 2.0 * newvar[2, :] - newvar[3, :]
+            newvar[0, :] = 2.0 * newvar[1, :] - newvar[2, :]
+            newvar[-2, :] = 2.0 * newvar[-3, :] - newvar[-4, :]
+            newvar[-1, :] = 2.0 * newvar[-2, :] - newvar[-3, :]
+        else:
+            # Constant extrapolation
+            newvar[1, :] = newvar[2, :]
+            newvar[0, :] = newvar[2, :]
+            newvar[-2, :] = newvar[-3, :]
+            newvar[-1, :] = newvar[-3, :]
+        grd[name] = newvar
+
+    # Calculate metric tensor
+    grd.update(calcMetric(grd, bpsign, verbose, ignore_checks))
+
+    dphidy = grd["dphidy"]
+    Rxy = grd["Rxy"]
+    Zxy = grd["Zxy"]
+    nx, ny = Rxy.shape
+
+    # Grid indices
+    ixseps1 = g["iyseparatrix1"] + 2  # Lower X-point separatrix
+    ixseps2 = min(g["iyseparatrix2"] + 2, nx)  # Upper X-point separatrix
+    jyseps1_1 = g["ix_cut1"] - 1
+    jyseps2_1 = g["ix_cut2"]
+    ny_inner = g["ix_inner"]
+    jyseps1_2 = g["ix_cut3"]
+    jyseps2_2 = g["ix_cut4"] - 1
+
+    # Calculate zShift and ShiftAngle
+    zShift = np.zeros((nx, ny))
+
+    # Inner core region
+    zShift[:, jyseps1_1 + 1] = 0.5 * dphidy[:, jyseps1_1 + 1] * dy
+    # Note: This goes to jyseps2_1 + 1, the first cell in the upper inner leg (including upper PF)
+    for jy in range(jyseps1_1 + 2, jyseps2_1 + 2):
+        zShift[:, jy] = (
+            zShift[:, jy - 1] + 0.5 * (dphidy[:, jy] + dphidy[:, jy - 1]) * dy
+        )
+
+    # Outer core
+    # Note: ixseps2 points are connected from inner to outer core
+    # If single null, ixseps2 = nx, jyseps1_2 = jyseps2_1
+    zShift[:, jyseps1_2 + 1] = (
+        zShift[:, jyseps2_1]
+        + 0.5 * (dphidy[:, jyseps2_1] + dphidy[:, jyseps1_2 + 1]) * dy
+    )
+    for jy in range(jyseps1_2 + 2, jyseps2_2 + 2):
+        zShift[:, jy] = (
+            zShift[:, jy - 1] + 0.5 * (dphidy[:, jy] + dphidy[:, jy - 1]) * dy
+        )
+
+    if jyseps2_1 != jyseps1_2:
+        # Double null, with upper X-point
+
+        # Upper inner leg. Note that upper PF region set from core at y index jyseps2_1 + 1
+        for jy in range(jyseps2_1 + 2, ny_inner):
+            zShift[:, jy] = (
+                zShift[:, jy - 1] + 0.5 * (dphidy[:, jy] + dphidy[:, jy - 1]) * dy
+            )
+
+        # Upper outer leg. Joins onto upper inner leg (jyseps2_1 + 1) for x < ixseps2
+        zShift[:ixseps2, jyseps1_2] = (
+            zShift[:ixseps2, jyseps2_1 + 1]
+            - 0.5 * (dphidy[:ixseps2, jyseps2_1 + 1] + dphidy[:ixseps2, jyseps1_2]) * dy
+        )
+        # joins outer core/SOL region for x >= ixseps2
+        zShift[ixseps2:, jyseps1_2] = (
+            zShift[ixseps2:, jyseps1_2 + 1]
+            - 0.5 * (dphidy[ixseps2:, jyseps1_2 + 1] + dphidy[ixseps2:, jyseps1_2]) * dy
+        )
+        # Iterate backwards along upper outer leg from X-point towards target
+        for jy in range(jyseps1_2 - 1, ny_inner - 1, -1):
+            zShift[:, jy] = (
+                zShift[:, jy + 1] - 0.5 * (dphidy[:, jy] + dphidy[:, jy - 1]) * dy
+            )
+
+    # Lower outer leg
+    zShift[:ixseps1, jyseps2_2 + 1] = 0.5 * dphidy[:ixseps1, jyseps2_2 + 1] * dy
+    for jy in range(jyseps2_2 + 2, ny):
+        zShift[:, jy] = (
+            zShift[:, jy - 1] + 0.5 * (dphidy[:, jy] + dphidy[:, jy - 1]) * dy
+        )
+    # Lower inner leg. Going backwards in Y toward the plate
+    zShift[:ixseps1, jyseps1_1] = -0.5 * dphidy[:ixseps1, jyseps1_1] * dy
+    for jy in range(jyseps1_1 - 1, -1, -1):
+        zShift[:, jy] = (
+            zShift[:, jy + 1] - 0.5 * (dphidy[:, jy] + dphidy[:, jy + 1]) * dy
+        )
+
+    ShiftAngle = np.zeros(nx)
+    ShiftAngle[:ixseps1] = np.sum(
+        dphidy[:ixseps1, (jyseps1_1 + 1) : (jyseps2_1 + 1)] * dy, axis=1  # Inner core
+    ) + np.sum(
+        dphidy[:ixseps1, (jyseps1_2 + 1) : (jyseps2_2 + 1)] * dy, axis=1  # Outer core
+    )
+
+    if verbose:
+        print(
+            "Safety factor: min {}, mean {}, max {}".format(
+                np.amin(ShiftAngle[:ixseps1]) / (2 * np.pi),
+                np.mean(ShiftAngle[:ixseps1]) / (2 * np.pi),
+                np.amax(ShiftAngle[:ixseps1]) / (2 * np.pi),
+            )
+        )
+
+    if plotting:
+        plt.plot(Rxy, Zxy, "x")
+        plt.plot(Rxy[ixseps1, :], Zxy[ixseps1, :], color="b", label="ixseps1")
+        if ixseps2 < nx:
+            plt.plot(Rxy[ixseps2, :], Zxy[ixseps2, :], color="r", label="ixseps2")
+
+        plt.plot(Rxy[:, jyseps1_1], Zxy[:, jyseps1_1], color="k", label="jyseps1_1")
+        plt.plot(Rxy[:, jyseps1_2], Zxy[:, jyseps1_2], color="b", label="jyseps1_2")
+        plt.plot(Rxy[:, ny_inner], Zxy[:, ny_inner], color="c", label="ny_inner")
+        plt.plot(Rxy[:, jyseps2_1], Zxy[:, jyseps2_1], color="g", label="jyseps2_1")
+        plt.plot(Rxy[:, jyseps2_2], Zxy[:, jyseps2_2], color="r", label="jyseps2_2")
+
+        ax = plt.gca()
+        ax.set_aspect("equal", adjustable="box")
+        ax.set_xlabel("R")
+        ax.set_ylabel("Z")
+        plt.legend()
+        plt.show()
+
+    from boututils.datafile import DataFile
+
+    if verbose:
+        print("Saving to " + output_filename)
+
+    with DataFile(output_filename, create=True, format="NETCDF4") as f:
+        # Save unique ID for grid file
+        import uuid
+
+        f.write_file_attribute("grid_id", str(uuid.uuid1()))
+        f.write_file_attribute("gridue", str(gridue_file))
+
+        f.write("nx", nx)
+        f.write("ny", ny)
+        f.write("ixseps1", ixseps1)
+        f.write("ixseps2", ixseps2)
+        f.write("jyseps1_1", jyseps1_1)
+        f.write("jyseps2_1", jyseps2_1)
+        f.write("ny_inner", ny_inner)
+        f.write("jyseps1_2", jyseps1_2)
+        f.write("jyseps2_2", jyseps2_2)
+
+        # 2D fields
+        for name in grd:
+            f.write(name, grd[name])
+
+        f.write("zShift", zShift)
+        f.write("ShiftAngle", ShiftAngle)