diff --git a/.gitignore b/.gitignore
index e041081..8d05456 100644
--- a/.gitignore
+++ b/.gitignore
@@ -169,6 +169,8 @@ examples/permz_*,1,*.png
 examples/poro_*,1,*.png
 examples/porv_*,1,*.png
 examples/satnum_*,1,*.png
+examples/sgas_*,1,*.png
+examples/fgip.png
 tests/generic_deck
 
 # Python environment
diff --git a/MANIFEST.in b/MANIFEST.in
index b539b81..98fd363 100755
--- a/MANIFEST.in
+++ b/MANIFEST.in
@@ -1,3 +1,4 @@
 include README.md
 include LICENSE
-include src/plopm/core/*
\ No newline at end of file
+include src/plopm/core/*
+include src/plopm/utils/*
\ No newline at end of file
diff --git a/src/plopm/core/plopm.py b/src/plopm/core/plopm.py
index 3e0f701..b2ea3c5 100755
--- a/src/plopm/core/plopm.py
+++ b/src/plopm/core/plopm.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 # SPDX-FileCopyrightText: 2024 NORCE
 # SPDX-License-Identifier: GPL-3.0
 # pylint: disable=R1702, W0123, W1401
@@ -8,26 +7,20 @@
 """
 
 import argparse
-import csv
-import os
-from io import StringIO
 import numpy as np
-import matplotlib
-import matplotlib.pyplot as plt
-from mpl_toolkits.axes_grid1 import make_axes_locatable
-
-try:
-    from opm.io.ecl import EclFile as OpmFile
-    from opm.io.ecl import EGrid as OpmGrid
-    from opm.io.ecl import ESmry as OpmSummary
-except ImportError:
-    print("The Python package opm was not found, using resdata")
-try:
-    from resdata.resfile import ResdataFile
-    from resdata.summary import Summary
-    from resdata.grid import Grid
-except ImportError:
-    print("The resdata Python package was not found, using opm")
+from plopm.utils.initialization import ini_dic, ini_properties, ini_readers
+from plopm.utils.readers import (
+    get_kws,
+    get_wells,
+    get_xycoords_resdata,
+    get_xycoords_opm,
+    get_xzcoords_resdata,
+    get_xzcoords_opm,
+    get_yzcoords_resdata,
+    get_yzcoords_opm,
+)
+from plopm.utils.mapping import map_xycoords, map_xzcoords, map_yzcoords
+from plopm.utils.write import make_summary, make_2dmaps
 
 
 def plopm():
@@ -36,129 +29,26 @@ def plopm():
     dic = ini_dic(cmdargs)
     ini_properties(dic)
     ini_readers(dic)
-    if not os.path.exists(dic["output"]):
-        os.system(f"mkdir {dic['output']}")
-    if dic["slide"][0] > -1:
-        dic["mx"], dic["my"] = 2 * dic["ny"] - 1, 2 * dic["nz"] - 1
-        dic["xmeaning"], dic["ymeaning"] = "y", "z"
-    elif dic["slide"][1] > -1:
-        dic["mx"], dic["my"] = 2 * dic["nx"] - 1, 2 * dic["nz"] - 1
-        dic["xmeaning"], dic["ymeaning"] = "x", "z"
-    else:
-        dic["mx"], dic["my"] = 2 * dic["nx"] - 1, 2 * dic["ny"] - 1
-        dic["xmeaning"], dic["ymeaning"] = "x", "y"
-    dic["wx"], dic["wy"] = 2 * dic["nx"] - 1, 2 * dic["ny"] - 1
     get_kws(dic)
     if len(dic["vsum"]) > 0:
         make_summary(dic)
         return
-    get_mesh(dic)
-    if dic["wells"] == 1:
+    if dic["well"] == 1:
         get_wells(dic)
-        get_xy_wells(dic)
-        make_wells(dic)
-        return
+    if dic["grid"] == 1:
+        dic["grida"] = np.ones((dic["mx"]) * (dic["my"])) * np.nan
     if dic["slide"][0] > -1:
-        dic["tslide"] = f", slide i={dic['slide'][0]+1}"
-        dic["nslide"] = f"{dic['slide'][0]+1},*,*"
-        get_yzcoords(dic)
-        map_yzcoords(dic)
+        handle_slide_x(dic)
     elif dic["slide"][1] > -1:
-        dic["tslide"] = f", slide j={dic['slide'][1]+1}"
-        dic["nslide"] = f"*,{dic['slide'][1]+1},*"
-        get_xzcoords(dic)
-        map_xzcoords(dic)
+        handle_slide_y(dic)
     else:
-        dic["tslide"] = f", slide k={dic['slide'][2]+1}"
-        dic["nslide"] = f"*,*,{dic['slide'][2]+1}"
-        get_xycoords(dic)
-        map_xycoords(dic)
+        handle_slide_z(dic)
     make_2dmaps(dic)
 
 
-def make_summary(dic):
+def handle_slide_x(dic):
     """
-    Plot the summary variable
-
-    Args:
-        cmdargs (dict): Command arguments
-
-    Returns:
-        None
-
-    """
-    plt.rcParams.update({"axes.grid": True})
-    fig, axis = plt.subplots()
-    if len(dic["cmaps"]) == 1:
-        axis.step(dic["time"], dic["vsum"], color=dic["cmaps"][0])
-    else:
-        axis.step(dic["time"], dic["vsum"], color="b")
-    axis.set_ylabel(dic["props"][0])
-    axis.set_xlabel("Time [s]")
-    axis.set_xlim([min(dic["time"]), max(dic["time"])])
-    axis.set_ylim([min(dic["vsum"]), max(dic["vsum"])])
-    axis.set_xticks(
-        np.linspace(
-            min(dic["time"]),
-            max(dic["time"]),
-            4,
-        )
-    )
-    axis.set_yticks(
-        np.linspace(
-            min(dic["vsum"]),
-            max(dic["vsum"]),
-            4,
-        )
-    )
-    fig.savefig(f"{dic['output']}/{dic['variable']}.png", bbox_inches="tight", dpi=300)
-    plt.close()
-
-
-def ini_dic(cmdargs):
-    """
-    Initialize the global dictionary
-
-    Args:
-        cmdargs (dict): Command arguments
-
-    Returns:
-        dic (dict): Modified global dictionary
-
-    """
-    dic = {"name": cmdargs["input"].strip()}
-    dic["coords"] = ["x", "y", "z"]
-    dic["scale"] = cmdargs["scale"].strip()
-    dic["use"] = cmdargs["use"].strip()
-    dic["variable"] = cmdargs["variable"].strip()
-    dic["size"] = float(cmdargs["size"])
-    dic["legend"] = cmdargs["legend"].strip()
-    dic["numbers"] = cmdargs["numbers"].strip()
-    dic["colormap"] = cmdargs["colormap"].strip()
-    dic["output"] = cmdargs["output"].strip()
-    dic["xlim"], dic["ylim"], dic["wellsij"], dic["vsum"] = [], [], [], []
-    dic["dtitle"] = ""
-    dic["restart"] = int(cmdargs["restart"])
-    dic["wells"] = int(cmdargs["wells"])
-    if dic["restart"] == -1:
-        dic["restart"] = 0
-    dic["slide"] = (
-        np.genfromtxt(StringIO(cmdargs["slide"]), delimiter=",", dtype=int) - 1
-    )
-    if cmdargs["xlim"]:
-        dic["xlim"] = np.genfromtxt(
-            StringIO(cmdargs["xlim"]), delimiter=",", dtype=float
-        )
-    if cmdargs["ylim"]:
-        dic["ylim"] = np.genfromtxt(
-            StringIO(cmdargs["ylim"]), delimiter=",", dtype=float
-        )
-    return dic
-
-
-def ini_readers(dic):
-    """
-    Set the classes for reading the properties
+    Processing the selected yz slide to obtain the grid properties
 
     Args:
         dic (dict): Global dictionary
@@ -167,208 +57,21 @@ def ini_readers(dic):
         dic (dict): Modified global dictionary
 
     """
+    dic["tslide"] = f", slide i={dic['slide'][0]+1}"
+    dic["nslide"] = f"{dic['slide'][0]+1},*,*"
+    for well in reversed(dic["wells"]):
+        if well[0] != dic["slide"][0]:
+            dic["wells"].remove(well)
     if dic["use"] == "resdata":
-        dic["egrid"] = Grid(f"{dic['name']}.EGRID")
-        dic["init"] = ResdataFile(f"{dic['name']}.INIT")
-        dic["nx"] = dic["egrid"].nx
-        dic["ny"] = dic["egrid"].ny
-        dic["nz"] = dic["egrid"].nz
-        if os.path.isfile(f"{dic['name']}.UNRST"):
-            dic["unrst"] = ResdataFile(f"{dic['name']}.UNRST")
-        if os.path.isfile(f"{dic['name']}.SMSPEC"):
-            dic["summary"] = Summary(f"{dic['name']}.SMSPEC")
-    else:
-        dic["egrid"] = OpmGrid(f"{dic['name']}.EGRID")
-        dic["init"] = OpmFile(f"{dic['name']}.INIT")
-        dic["nx"] = dic["egrid"].dimension[0]
-        dic["ny"] = dic["egrid"].dimension[1]
-        dic["nz"] = dic["egrid"].dimension[2]
-        if os.path.isfile(f"{dic['name']}.UNRST"):
-            dic["unrst"] = OpmFile(f"{dic['name']}.UNRST")
-        if os.path.isfile(f"{dic['name']}.SMSPEC"):
-            dic["summary"] = OpmSummary(f"{dic['name']}.SMSPEC")
-
-
-def ini_properties(dic):
-    """
-    Define the properties to plot
-
-    Args:
-        dic (dict): Global dictionary
-
-    Returns:
-        dic (dict): Modified global dictionary
-
-    """
-    if dic["variable"]:
-        dic["props"] = [dic["variable"]]
-        if dic["variable"].lower() in ["disperc", "depth", "dx", "dy", "dz"]:
-            dic["units"] = [" [m]"]
-            dic["cmaps"] = ["jet"]
-            dic["format"] = [lambda x, _: f"{x:.2e}"]
-        elif dic["variable"].lower() in ["porv"]:
-            dic["units"] = [r" [m$^3$]"]
-            dic["cmaps"] = ["terrain"]
-            dic["format"] = [lambda x, _: f"{x:.2e}"]
-        elif dic["variable"].lower() in ["permx", "permy", "permz"]:
-            dic["units"] = [" [mD]"]
-            dic["cmaps"] = ["turbo"]
-            dic["format"] = [lambda x, _: f"{x:.0f}"]
-        elif "num" in dic["variable"].lower():
-            dic["units"] = [" [-]"]
-            dic["cmaps"] = ["tab20b"]
-            dic["format"] = [lambda x, _: f"{x:.2f}"]
-        else:
-            dic["units"] = [" [-]"]
-            dic["cmaps"] = ["gnuplot"]
-            dic["format"] = [lambda x, _: f"{x:.0f}"]
-        if dic["legend"]:
-            dic["units"] = [f" {dic['legend']}"]
-        if dic["numbers"]:
-            dic["format"] = [eval(dic["numbers"])]
-        if dic["colormap"]:
-            dic["cmaps"] = [dic["colormap"]]
+        get_yzcoords_resdata(dic)
     else:
-        dic["props"] = [
-            "porv",
-            "permx",
-            "permz",
-            "poro",
-            "fipnum",
-            "satnum",
-        ]
-        dic["units"] = [
-            r" [m$^3$]",
-            " [mD]",
-            " [mD]",
-            " [-]",
-            " [-]",
-            " [-]",
-        ]
-        dic["format"] = [
-            lambda x, _: f"{x:.2e}",
-            lambda x, _: f"{x:.0f}",
-            lambda x, _: f"{x:.0f}",
-            lambda x, _: f"{x:.1f}",
-            lambda x, _: f"{x:.0f}",
-            lambda x, _: f"{x:.0f}",
-        ]
-        dic["cmaps"] = [
-            "terrain",
-            "turbo",
-            "turbo",
-            "gnuplot",
-            "tab20b",
-            "tab20b",
-        ]
-    font = {"family": "normal", "weight": "normal", "size": dic["size"]}
-    matplotlib.rc("font", **font)
-    plt.rcParams.update(
-        {
-            "text.usetex": True,
-            "font.family": "monospace",
-            "legend.columnspacing": 0.9,
-            "legend.handlelength": 3.5,
-            "legend.fontsize": dic["size"],
-            "lines.linewidth": 4,
-            "axes.titlesize": dic["size"],
-            "axes.grid": False,
-            "figure.figsize": (8, 16),
-        }
-    )
+        get_yzcoords_opm(dic)
+    map_yzcoords(dic)
 
 
-def get_yzcoords(dic):
+def handle_slide_y(dic):
     """
-    Handle the coordinates from the OPM Grid to the 2D yz-mesh
-
-    Args:
-        dic (dict): Global dictionary
-
-    Returns:
-        dic (dict): Modified global dictionary
-
-    """
-    for j in range(dic["nz"]):
-        dic["xc"].append([])
-        dic["yc"].append([])
-        for k, n in zip(["x", "x", "y", "y"], [1, 7, 2, 8]):
-            dic[f"{k}c"][-1].append(
-                dic["mesh"][(dic["nz"] - j - 1) * dic["nx"] * dic["ny"]][n]
-            )
-        for i in range(dic["ny"] - 1):
-            for k, n in zip(["x", "x", "y", "y"], [1, 7, 2, 8]):
-                dic[f"{k}c"][-1].append(
-                    dic["mesh"][
-                        (i + 1) * dic["nx"]
-                        + (dic["nz"] - j - 1) * dic["nx"] * dic["ny"]
-                    ][n]
-                )
-        dic["xc"].append([])
-        dic["yc"].append([])
-        for k, n in zip(["x", "x", "y", "y"], [13, 19, 14, 20]):
-            dic[f"{k}c"][-1].append(
-                dic["mesh"][(dic["nz"] - j - 1) * dic["nx"] * dic["ny"]][n]
-            )
-        for i in range(dic["ny"] - 1):
-            for k, n in zip(["x", "x", "y", "y"], [13, 19, 14, 20]):
-                dic[f"{k}c"][-1].append(
-                    dic["mesh"][
-                        (i + 1) * dic["nx"]
-                        + (dic["nz"] - j - 1) * dic["nx"] * dic["ny"]
-                    ][n]
-                )
-
-
-def map_yzcoords(dic):
-    """
-    Map the properties from the simulations to the 2D slide
-
-    Args:
-        dic (dict): Global dictionary
-
-    Returns:
-        dic (dict): Modified global dictionary
-
-    """
-    if dic["use"] == "resdata":
-        for cell in dic["egrid"].cells():
-            if cell.active and cell.i == dic["slide"][0]:
-                for name in dic["props"]:
-                    dic[name + "a"][
-                        2 * cell.j + 2 * (dic["nz"] - cell.k - 1) * dic["mx"]
-                    ] = dic[name][cell.active_index]
-                if "porv" in dic["props"]:
-                    dic["porva"][
-                        2 * cell.j + 2 * (dic["nz"] - cell.k - 1) * dic["mx"]
-                    ] = dic["porv"][cell.global_index]
-    else:
-        for k in range(dic["nz"]):
-            for j in range(dic["ny"]):
-                for i in range(dic["nx"]):
-                    if (
-                        dic["porv"][i + j * dic["nx"] + k * dic["nx"] * dic["ny"]] > 0
-                        and i == dic["slide"][0]
-                    ):
-                        for name in dic["props"]:
-                            dic[name + "a"][
-                                2 * j + 2 * (dic["nz"] - k - 1) * dic["mx"]
-                            ] = dic[name][
-                                dic["actind"][
-                                    i + j * dic["nx"] + k * dic["nx"] * dic["ny"]
-                                ]
-                            ]
-                        if "porv" in dic["props"]:
-                            dic["porva"][
-                                2 * j + 2 * (dic["nz"] - k - 1) * dic["mx"]
-                            ] = dic["porv"][
-                                i + j * dic["nx"] + k * dic["nx"] * dic["ny"]
-                            ]
-
-
-def get_mesh(dic):
-    """
-    Read the coordinates using either opm or resdata
+    Processing the selected xz slide to obtain the grid properties
 
     Args:
         dic (dict): Global dictionary
@@ -377,62 +80,21 @@ def get_mesh(dic):
         dic (dict): Modified global dictionary
 
     """
+    dic["tslide"] = f", slide j={dic['slide'][1]+1}"
+    dic["nslide"] = f"*,{dic['slide'][1]+1},*"
+    for well in reversed(dic["wells"]):
+        if well[1] != dic["slide"][1]:
+            dic["wells"].remove(well)
     if dic["use"] == "resdata":
-        dic["mesh"] = dic["egrid"].export_corners(dic["egrid"].export_index())
+        get_xzcoords_resdata(dic)
     else:
-        dic["mesh"] = []
-        for k in range(dic["nz"]):
-            for j in range(dic["ny"]):
-                for i in range(dic["nx"]):
-                    dic["mesh"].append([])
-                    for n in range(8):
-                        dic["mesh"][-1] += [
-                            dic["egrid"].xyz_from_ijk(i, j, k)[0][n],
-                            dic["egrid"].xyz_from_ijk(i, j, k)[1][n],
-                            dic["egrid"].xyz_from_ijk(i, j, k)[2][n],
-                        ]
-    dic["xc"], dic["yc"] = [], []
-
-
-def get_xzcoords(dic):
-    """
-    Handle the coordinates from the OPM Grid to the 2D xz-mesh
-
-    Args:
-        dic (dict): Global dictionary
-
-    Returns:
-        dic (dict): Modified global dictionary
-
-    """
-    for j in range(dic["nz"]):
-        dic["xc"].append([])
-        dic["yc"].append([])
-        for k, n in zip(["x", "x", "y", "y"], [0, 3, 2, 5]):
-            dic[f"{k}c"][-1].append(
-                dic["mesh"][(dic["nz"] - j - 1) * dic["nx"] * dic["ny"]][n]
-            )
-        for i in range(dic["nx"] - 1):
-            for k, n in zip(["x", "x", "y", "y"], [0, 3, 2, 5]):
-                dic[f"{k}c"][-1].append(
-                    dic["mesh"][1 + i + (dic["nz"] - j - 1) * dic["nx"] * dic["ny"]][n]
-                )
-        dic["xc"].append([])
-        dic["yc"].append([])
-        for k, n in zip(["x", "x", "y", "y"], [12, 15, 14, 17]):
-            dic[f"{k}c"][-1].append(
-                dic["mesh"][(dic["nz"] - j - 1) * dic["nx"] * dic["ny"]][n]
-            )
-        for i in range(dic["nx"] - 1):
-            for k, n in zip(["x", "x", "y", "y"], [12, 15, 14, 17]):
-                dic[f"{k}c"][-1].append(
-                    dic["mesh"][1 + i + (dic["nz"] - j - 1) * dic["nx"] * dic["ny"]][n]
-                )
+        get_xzcoords_opm(dic)
+    map_xzcoords(dic)
 
 
-def map_xzcoords(dic):
+def handle_slide_z(dic):
     """
-    Map the properties from the simulations to the 2D slide
+    Processing the selected xy slide to obtain the grid properties
 
     Args:
         dic (dict): Global dictionary
@@ -441,399 +103,16 @@ def map_xzcoords(dic):
         dic (dict): Modified global dictionary
 
     """
+    dic["tslide"] = f", slide k={dic['slide'][2]+1}"
+    dic["nslide"] = f"*,*,{dic['slide'][2]+1}"
+    for well in reversed(dic["wells"]):
+        if dic["slide"][2] not in range(well[2], well[3] + 1):
+            dic["wells"].remove(well)
     if dic["use"] == "resdata":
-        for cell in dic["egrid"].cells():
-            if cell.active and cell.j == dic["slide"][1]:
-                for name in dic["props"]:
-                    dic[name + "a"][
-                        2 * cell.i + 2 * (dic["nz"] - cell.k - 1) * dic["mx"]
-                    ] = dic[name][cell.active_index]
-                if "porv" in dic["props"]:
-                    dic["porva"][
-                        2 * cell.i + 2 * (dic["nz"] - cell.k - 1) * dic["mx"]
-                    ] = dic["porv"][cell.global_index]
+        get_xycoords_resdata(dic)
     else:
-        for k in range(dic["nz"]):
-            for j in range(dic["ny"]):
-                for i in range(dic["nx"]):
-                    if (
-                        dic["porv"][i + j * dic["nx"] + k * dic["nx"] * dic["ny"]] > 0
-                        and j == dic["slide"][1]
-                    ):
-                        for name in dic["props"]:
-                            dic[name + "a"][
-                                2 * i + 2 * (dic["nz"] - k - 1) * dic["mx"]
-                            ] = dic[name][
-                                dic["actind"][
-                                    i + j * dic["nx"] + k * dic["nx"] * dic["ny"]
-                                ]
-                            ]
-                        if "porv" in dic["props"]:
-                            dic["porva"][
-                                2 * i + 2 * (dic["nz"] - k - 1) * dic["mx"]
-                            ] = dic["porv"][
-                                i + j * dic["nx"] + k * dic["nx"] * dic["ny"]
-                            ]
-
-
-def get_xycoords(dic):
-    """
-    Handle the coordinates from the OPM Grid to the 2D xy-mesh
-
-    Args:
-        dic (dict): Global dictionary
-
-    Returns:
-        dic (dict): Modified global dictionary
-
-    """
-    for j in range(dic["ny"]):
-        dic["xc"].append([])
-        dic["yc"].append([])
-        for k, n in zip(["x", "x", "y", "y"], [0, 3, 1, 4]):
-            dic[f"{k}c"][-1].append(dic["mesh"][j * dic["nx"]][n])
-        for i in range(dic["nx"] - 1):
-            for k, n in zip(["x", "x", "y", "y"], [0, 3, 1, 4]):
-                dic[f"{k}c"][-1].append(dic["mesh"][1 + i + j * dic["nx"]][n])
-        dic["xc"].append([])
-        dic["yc"].append([])
-        for k, n in zip(["x", "x", "y", "y"], [6, 9, 7, 10]):
-            dic[f"{k}c"][-1].append(dic["mesh"][j * dic["nx"]][n])
-        for i in range(dic["nx"] - 1):
-            for k, n in zip(["x", "x", "y", "y"], [6, 9, 7, 10]):
-                dic[f"{k}c"][-1].append(dic["mesh"][1 + i + j * dic["nx"]][n])
-
-
-def map_xycoords(dic):
-    """
-    Map the properties from the simulations to the 2D slide
-
-    Args:
-        dic (dict): Global dictionary
-
-    Returns:
-        dic (dict): Modified global dictionary
-
-    """
-    if dic["use"] == "resdata":
-        for cell in dic["egrid"].cells():
-            if cell.active and cell.k == dic["slide"][2]:
-                for name in dic["props"]:
-                    dic[name + "a"][2 * cell.i + 2 * cell.j * dic["mx"]] = dic[name][
-                        cell.active_index
-                    ]
-                if "porv" in dic["props"]:
-                    dic["porva"][2 * cell.i + 2 * cell.j * dic["mx"]] = dic["porv"][
-                        cell.global_index
-                    ]
-    else:
-        for k in range(dic["nz"]):
-            for j in range(dic["ny"]):
-                for i in range(dic["nx"]):
-                    if (
-                        dic["porv"][i + j * dic["nx"] + k * dic["nx"] * dic["ny"]] > 0
-                        and k == dic["slide"][2]
-                    ):
-                        for name in dic["props"]:
-                            dic[name + "a"][2 * i + 2 * j * dic["mx"]] = dic[name][
-                                dic["actind"][
-                                    i + j * dic["nx"] + k * dic["nx"] * dic["ny"]
-                                ]
-                            ]
-                        if "porv" in dic["props"]:
-                            dic["porva"][2 * i + 2 * j * dic["mx"]] = dic["porv"][
-                                i + j * dic["nx"] + k * dic["nx"] * dic["ny"]
-                            ]
-
-
-def get_kws(dic):
-    """
-    Get the static properties from the INIT file using ResdataFile
-
-    Args:
-        dic (dict): Global dictionary
-
-    Returns:
-        dic (dict): Modified global dictionary
-
-    """
-    for name in dic["props"]:
-        if dic["use"] == "resdata":
-            if dic["init"].has_kw(name.upper()):
-                dic[name] = np.array(dic["init"].iget_kw(name.upper())[0])
-            elif dic["unrst"].has_kw(name.upper()):
-                dic[name] = np.array(
-                    dic["unrst"].iget_kw(name.upper())[dic["restart"] - 1]
-                )
-                ntot = len(dic["unrst"].iget_kw(name.upper()))
-                dic["dtitle"] = (
-                    f", rst {ntot if dic['restart']==0 else dic['restart']} out of {ntot}"
-                )
-            elif dic["summary"].has_key(name.upper()):
-                dic["vsum"] = dic["summary"][name.upper()].values
-                dic["time"] = dic["summary"]["TIME"].values
-                return
-        else:
-            if dic["init"].count(name.upper()):
-                dic[name] = np.array(dic["init"][name.upper()])
-            elif dic["unrst"].count(name.upper()):
-                nrst = (
-                    dic["restart"] - 1
-                    if dic["restart"] > 0
-                    else dic["unrst"].count(name.upper()) - 1
-                )
-                dic[name] = np.array(dic["unrst"][name.upper(), nrst])
-                ntot = dic["unrst"].count(name.upper())
-                dic["dtitle"] = (
-                    f", rst {ntot if dic['restart']==0 else dic['restart']} out of {ntot}"
-                )
-            elif name.upper() in dic["summary"].keys():
-                dic["vsum"] = dic["summary"][name.upper()]
-                dic["time"] = dic["summary"]["TIME"]
-                return
-        dic[name + "a"] = np.ones(dic["mx"] * dic["my"]) * np.nan
-    if dic["use"] == "opm":
-        dic["porv"] = np.array(dic["init"]["PORV"])
-        dic["actind"] = np.cumsum([1 if porv > 0 else 0 for porv in dic["porv"]]) - 1
-
-
-def get_xy_wells(dic):
-    """
-    Get the top x,y coordinates from the geological model for the well figure
-
-    Args:
-        dic (dict): Global dictionary
-
-    Returns:
-        dic (dict): Modified global dictionary
-
-    """
-    dic["mesh"] = dic["egrid"].export_corners(dic["egrid"].export_index())
-    (
-        dic["xw"],
-        dic["yw"],
-    ) = (
-        [],
-        [],
-    )
-    for j in range(dic["ny"]):
-        dic["xw"].append([])
-        dic["yw"].append([])
-        for k, n in zip(["x", "x", "y", "y"], [0, 3, 1, 4]):
-            dic[f"{k}w"][-1].append(dic["mesh"][j * dic["nx"]][n])
-        for i in range(dic["nx"] - 1):
-            for k, n in zip(["x", "x", "y", "y"], [0, 3, 1, 4]):
-                dic[f"{k}w"][-1].append(dic["mesh"][1 + i + j * dic["nx"]][n])
-        dic["xw"].append([])
-        dic["yw"].append([])
-        for k, n in zip(["x", "x", "y", "y"], [6, 9, 7, 10]):
-            dic[f"{k}w"][-1].append(dic["mesh"][j * dic["nx"]][n])
-        for i in range(dic["nx"] - 1):
-            for k, n in zip(["x", "x", "y", "y"], [6, 9, 7, 10]):
-                dic[f"{k}w"][-1].append(dic["mesh"][1 + i + j * dic["nx"]][n])
-    for cell in dic["egrid"].cells():
-        if cell.active:
-            dic["wellsa"][2 * cell.i + 2 * cell.j * dic["wx"]] = len(dic["wellsij"])
-    for i, well in enumerate(dic["wellsij"]):
-        dic["wellsa"][2 * well[1] + 2 * well[2] * dic["wx"]] = i
-
-
-def get_wells(dic):
-    """
-    Using the input deck (.DATA) to read the i,j well locations
-
-    Args:
-        dic (dict): Global dictionary
-
-    Returns:
-        dic (dict): Modified global dictionary
-
-    """
-    dic["wellsa"] = np.ones((dic["wx"]) * (dic["wy"])) * np.nan
-    wells, sources = False, False
-    with open(f"{dic['name']}.DATA", "r", encoding="utf8") as file:
-        for row in csv.reader(file):
-            nrwo = str(row)[2:-2]
-            if wells:
-                if len(nrwo.split()) < 2:
-                    break
-                dic["wellsij"].append(
-                    [
-                        nrwo.split()[0],
-                        int(nrwo.split()[2]) - 1,
-                        int(nrwo.split()[3]) - 1,
-                    ]
-                )
-            if sources:
-                if len(nrwo.split()) < 2:
-                    break
-                dic["wellsij"].append(
-                    [
-                        f"SOURCE{len(dic['wellsij'])+1}",
-                        int(nrwo.split()[0]) - 1,
-                        int(nrwo.split()[1]) - 1,
-                    ]
-                )
-            if nrwo == "WELSPECS":
-                wells = True
-            if nrwo == "SOURCE":
-                sources = True
-
-
-def make_2dmaps(dic):
-    """
-    Method to create the 2d maps using pcolormesh
-
-    Args:
-        dic (dict): Global dictionary
-
-    Returns:
-        dic (dict): Modified global dictionary
-
-    """
-    maps = np.ones([dic["my"], dic["mx"]]) * np.nan
-    for n, name in enumerate(dic["props"]):
-        fig, axis = plt.subplots()
-        for i in np.arange(0, dic["my"]):
-            maps[i, :] = dic[name + "a"][i * (dic["mx"]) : (i + 1) * (dic["mx"])]
-        imag = axis.pcolormesh(
-            dic["xc"],
-            dic["yc"],
-            maps,
-            shading="flat",
-            cmap=dic["cmaps"][n],
-        )
-        if dic["scale"] == "yes":
-            axis.axis("scaled")
-        axis.set_xlabel(f"{dic['xmeaning']} [m]")
-        axis.set_ylabel(f"{dic['ymeaning']} [m]")
-        extra = ""
-        if name == "porv":
-            extra = f" (sum={sum(dic[name]):.3e})"
-        axis.set_title(name + dic["tslide"] + dic["dtitle"] + extra)
-        minc = dic[name][~np.isnan(dic[name])].min()
-        maxc = dic[name][~np.isnan(dic[name])].max()
-        ntick = 5
-        if maxc == minc:
-            ntick = 1
-        elif name in ["fipnum", "satnum"]:
-            ntick = maxc
-        divider = make_axes_locatable(axis)
-        vect = np.linspace(
-            minc,
-            maxc,
-            ntick,
-            endpoint=True,
-        )
-        fig.colorbar(
-            imag,
-            cax=divider.append_axes("right", size="5%", pad=0.05),
-            orientation="vertical",
-            ticks=vect,
-            label=dic["units"][n],
-            format=dic["format"][n],
-        )
-        imag.set_clim(
-            minc,
-            maxc,
-        )
-        if dic["slide"][2] == -2:
-            axis.invert_yaxis()
-        axis.set_xticks(
-            np.linspace(
-                min(min(dic["xc"])),
-                max(max(dic["xc"])),
-                4,
-            )
-        )
-        axis.set_yticks(
-            np.linspace(
-                min(min(dic["yc"])),
-                max(max(dic["yc"])),
-                4,
-            )
-        )
-        if len(dic["xlim"]) > 1:
-            axis.set_xlim(dic["xlim"])
-        if len(dic["ylim"]) > 1:
-            axis.set_ylim(dic["ylim"])
-        fig.savefig(
-            f"{dic['output']}/{name}_{dic['nslide']}.png", bbox_inches="tight", dpi=300
-        )
-        plt.close()
-
-
-def make_wells(dic):
-    """
-    If there are any well, we only plot them on the top xy view
-
-     Args:
-         dic (dict): Global dictionary
-
-     Returns:
-         dic (dict): Modified global dictionary
-
-    """
-    maps = np.ones([dic["wy"], dic["wx"]]) * np.nan
-    fig, axis = plt.subplots()
-    for i in np.arange(0, dic["wy"]):
-        maps[i, :] = dic["wellsa"][i * (dic["wx"]) : (i + 1) * (dic["wx"])]
-    imag = axis.pcolormesh(
-        dic["xw"],
-        dic["yw"],
-        maps,
-        shading="flat",
-        cmap="nipy_spectral",
-    )
-    minc = 0
-    maxc = len(dic["wellsij"])
-    divider = make_axes_locatable(axis)
-    vect = np.linspace(
-        minc,
-        maxc,
-        5,
-        endpoint=True,
-    )
-    fig.colorbar(
-        imag,
-        cax=divider.append_axes("right", size="0%", pad=0.05),
-        orientation="vertical",
-        ticks=vect,
-        format=lambda x, _: "",
-    )
-    imag.set_clim(
-        minc,
-        maxc,
-    )
-    cmap = matplotlib.colormaps["nipy_spectral"]
-    colors = cmap(np.linspace(0, 1, len(dic["wellsij"]) + 1))
-    for i, well in enumerate(dic["wellsij"]):
-        plt.text(0, i, f"{i}-({well[1]+1},{well[2]+1})", c=colors[i], fontweight="bold")
-    if dic["scale"] == "yes":
-        axis.axis("scaled")
-    axis.set_title("Well's location, top view xy (k=1)")
-    axis.set_xlabel(f"{dic['xmeaning']} [m]")
-    axis.set_ylabel(f"{dic['ymeaning']} [m]")
-    axis.set_xticks(
-        np.linspace(
-            min(min(dic["xw"])),
-            max(max(dic["xw"])),
-            4,
-        )
-    )
-    axis.set_yticks(
-        np.linspace(
-            min(min(dic["yw"])),
-            max(max(dic["yw"])),
-            4,
-        )
-    )
-    if dic["xlim"]:
-        axis.set_xlim(dic["xlim"])
-    if dic["ylim"]:
-        axis.set_ylim(dic["ylim"])
-    fig.savefig(f"{dic['output']}/wells.png", bbox_inches="tight", dpi=300)
-    plt.close()
+        get_xycoords_opm(dic)
+    map_xycoords(dic)
 
 
 def load_parser():
@@ -930,7 +209,14 @@ def load_parser():
         "-w",
         "--wells",
         default=0,
-        help="Plot the xz-position of the wells or sources ('0' by default).",
+        help="Plot the positions of the wells or sources ('0' by default).",
+    )
+    parser.add_argument(
+        "-g",
+        "--grid",
+        default=0,
+        help="Plot information about the number of cells in the x, y, and z "
+        " directions and number of active grid cells ('0' by default).",
     )
     return vars(parser.parse_known_args()[0])
 
diff --git a/src/plopm/utils/initialization.py b/src/plopm/utils/initialization.py
new file mode 100755
index 0000000..69438e6
--- /dev/null
+++ b/src/plopm/utils/initialization.py
@@ -0,0 +1,209 @@
+# SPDX-FileCopyrightText: 2024 NORCE
+# SPDX-License-Identifier: GPL-3.0
+# pylint: disable=W0123
+
+"""
+Utiliy functions to set the requiried input values by plopm.
+"""
+
+import os
+from io import StringIO
+import numpy as np
+import matplotlib
+import matplotlib.pyplot as plt
+
+try:
+    from opm.io.ecl import EclFile as OpmFile
+    from opm.io.ecl import EGrid as OpmGrid
+    from opm.io.ecl import ESmry as OpmSummary
+except ImportError:
+    print("The Python package opm was not found, using resdata")
+try:
+    from resdata.resfile import ResdataFile
+    from resdata.summary import Summary
+    from resdata.grid import Grid
+except ImportError:
+    print("The resdata Python package was not found, using opm")
+
+
+def ini_dic(cmdargs):
+    """
+    Initialize the global dictionary
+
+    Args:
+        cmdargs (dict): Command arguments
+
+    Returns:
+        dic (dict): Modified global dictionary
+
+    """
+    dic = {"name": cmdargs["input"].strip()}
+    dic["coords"] = ["x", "y", "z"]
+    dic["scale"] = cmdargs["scale"].strip()
+    dic["use"] = cmdargs["use"].strip()
+    dic["variable"] = cmdargs["variable"].strip()
+    dic["size"] = float(cmdargs["size"])
+    dic["legend"] = cmdargs["legend"].strip()
+    dic["numbers"] = cmdargs["numbers"].strip()
+    dic["colormap"] = cmdargs["colormap"].strip()
+    dic["output"] = cmdargs["output"].strip()
+    dic["xlim"], dic["ylim"], dic["wells"], dic["vsum"] = [], [], [], []
+    dic["grid"] = []
+    dic["dtitle"] = ""
+    dic["restart"] = int(cmdargs["restart"])
+    dic["well"] = int(cmdargs["wells"])
+    dic["grid"] = int(cmdargs["grid"])
+    if dic["restart"] == -1:
+        dic["restart"] = 0
+    dic["slide"] = (
+        np.genfromtxt(StringIO(cmdargs["slide"]), delimiter=",", dtype=int) - 1
+    )
+    if cmdargs["xlim"]:
+        dic["xlim"] = np.genfromtxt(
+            StringIO(cmdargs["xlim"]), delimiter=",", dtype=float
+        )
+    if cmdargs["ylim"]:
+        dic["ylim"] = np.genfromtxt(
+            StringIO(cmdargs["ylim"]), delimiter=",", dtype=float
+        )
+    return dic
+
+
+def ini_readers(dic):
+    """
+    Set the classes for reading the properties
+
+    Args:
+        dic (dict): Global dictionary
+
+    Returns:
+        dic (dict): Modified global dictionary
+
+    """
+    if dic["use"] == "resdata":
+        dic["egrid"] = Grid(f"{dic['name']}.EGRID")
+        dic["init"] = ResdataFile(f"{dic['name']}.INIT")
+        dic["nx"] = dic["egrid"].nx
+        dic["ny"] = dic["egrid"].ny
+        dic["nz"] = dic["egrid"].nz
+        if os.path.isfile(f"{dic['name']}.UNRST"):
+            dic["unrst"] = ResdataFile(f"{dic['name']}.UNRST")
+        if os.path.isfile(f"{dic['name']}.SMSPEC"):
+            dic["summary"] = Summary(f"{dic['name']}.SMSPEC")
+    else:
+        dic["egrid"] = OpmGrid(f"{dic['name']}.EGRID")
+        dic["init"] = OpmFile(f"{dic['name']}.INIT")
+        dic["nx"] = dic["egrid"].dimension[0]
+        dic["ny"] = dic["egrid"].dimension[1]
+        dic["nz"] = dic["egrid"].dimension[2]
+        if os.path.isfile(f"{dic['name']}.UNRST"):
+            dic["unrst"] = OpmFile(f"{dic['name']}.UNRST")
+        if os.path.isfile(f"{dic['name']}.SMSPEC"):
+            dic["summary"] = OpmSummary(f"{dic['name']}.SMSPEC")
+    if not os.path.exists(dic["output"]):
+        os.system(f"mkdir {dic['output']}")
+    if dic["slide"][0] > -1:
+        dic["mx"], dic["my"] = 2 * dic["ny"] - 1, 2 * dic["nz"] - 1
+        dic["xmeaning"], dic["ymeaning"] = "y", "z"
+    elif dic["slide"][1] > -1:
+        dic["mx"], dic["my"] = 2 * dic["nx"] - 1, 2 * dic["nz"] - 1
+        dic["xmeaning"], dic["ymeaning"] = "x", "z"
+    else:
+        dic["mx"], dic["my"] = 2 * dic["nx"] - 1, 2 * dic["ny"] - 1
+        dic["xmeaning"], dic["ymeaning"] = "x", "y"
+
+
+def ini_properties(dic):
+    """
+    Define the properties to plot
+
+    Args:
+        dic (dict): Global dictionary
+
+    Returns:
+        dic (dict): Modified global dictionary
+
+    """
+    if dic["variable"]:
+        dic["props"] = [dic["variable"]]
+        if dic["variable"].lower() in ["disperc", "depth", "dx", "dy", "dz"]:
+            dic["units"] = [" [m]"]
+            dic["cmaps"] = ["jet"]
+            dic["format"] = [lambda x, _: f"{x:.2e}"]
+        elif dic["variable"].lower() in ["porv"]:
+            dic["units"] = [r" [m$^3$]"]
+            dic["cmaps"] = ["terrain"]
+            dic["format"] = [lambda x, _: f"{x:.2e}"]
+        elif dic["variable"].lower() in ["permx", "permy", "permz"]:
+            dic["units"] = [" [mD]"]
+            dic["cmaps"] = ["turbo"]
+            dic["format"] = [lambda x, _: f"{x:.0f}"]
+        elif "num" in dic["variable"].lower():
+            dic["units"] = [" [-]"]
+            dic["cmaps"] = ["tab20b"]
+            dic["format"] = [lambda x, _: f"{x:.2f}"]
+        else:
+            dic["units"] = [" [-]"]
+            dic["cmaps"] = ["gnuplot"]
+            dic["format"] = [lambda x, _: f"{x:.0f}"]
+        if dic["legend"]:
+            dic["units"] = [f" {dic['legend']}"]
+        if dic["numbers"]:
+            dic["format"] = [eval(dic["numbers"])]
+        if dic["colormap"]:
+            dic["cmaps"] = [dic["colormap"]]
+    else:
+        dic["props"] = [
+            "porv",
+            "permx",
+            "permz",
+            "poro",
+            "fipnum",
+            "satnum",
+        ]
+        dic["units"] = [
+            r" [m$^3$]",
+            " [mD]",
+            " [mD]",
+            " [-]",
+            " [-]",
+            " [-]",
+        ]
+        dic["format"] = [
+            lambda x, _: f"{x:.2e}",
+            lambda x, _: f"{x:.0f}",
+            lambda x, _: f"{x:.0f}",
+            lambda x, _: f"{x:.1f}",
+            lambda x, _: f"{x:.0f}",
+            lambda x, _: f"{x:.0f}",
+        ]
+        dic["cmaps"] = [
+            "terrain",
+            "turbo",
+            "turbo",
+            "gnuplot",
+            "tab20b",
+            "tab20b",
+        ]
+    if dic["well"] == 1 or dic["grid"] == 1:
+        dic["props"] = []
+        dic["grid"] = 1
+        dic["units"] = [" [-]"]
+        dic["cmaps"] = ["nipy_spectral"]
+        dic["format"] = [lambda x, _: f"{x:.0f}"]
+    font = {"family": "normal", "weight": "normal", "size": dic["size"]}
+    matplotlib.rc("font", **font)
+    plt.rcParams.update(
+        {
+            "text.usetex": True,
+            "font.family": "monospace",
+            "legend.columnspacing": 0.9,
+            "legend.handlelength": 3.5,
+            "legend.fontsize": dic["size"],
+            "lines.linewidth": 4,
+            "axes.titlesize": dic["size"],
+            "axes.grid": False,
+            "figure.figsize": (8, 16),
+        }
+    )
+    dic["xc"], dic["yc"] = [], []
diff --git a/src/plopm/utils/mapping.py b/src/plopm/utils/mapping.py
new file mode 100644
index 0000000..1ce7b2f
--- /dev/null
+++ b/src/plopm/utils/mapping.py
@@ -0,0 +1,103 @@
+# SPDX-FileCopyrightText: 2024 NORCE
+# SPDX-License-Identifier: GPL-3.0
+
+"""
+Utiliy function for the grid and locations in the geological models.
+"""
+
+
+def map_yzcoords(dic):
+    """
+    Map the properties from the simulations to the 2D slide
+
+    Args:
+        dic (dict): Global dictionary
+
+    Returns:
+        dic (dict): Modified global dictionary
+
+    """
+    for k in range(dic["nz"]):
+        for j in range(dic["ny"]):
+            ind = dic["slide"][0] + j * dic["nx"] + k * dic["nx"] * dic["ny"]
+            if dic["porv"][ind] > 0:
+                for name in dic["props"]:
+                    dic[name + "a"][2 * j + 2 * (dic["nz"] - k - 1) * dic["mx"]] = dic[
+                        name
+                    ][dic["actind"][ind]]
+                if "porv" in dic["props"]:
+                    dic["porva"][2 * j + 2 * (dic["nz"] - k - 1) * dic["mx"]] = dic[
+                        "porv"
+                    ][ind]
+                if len(dic["wells"]) > 0:
+                    dic["wellsa"][2 * j + 2 * (dic["nz"] - k - 1) * dic["mx"]] = len(
+                        dic["wells"]
+                    )
+                if dic["grid"] == 1:
+                    dic["grida"][2 * j + 2 * (dic["nz"] - k - 1) * dic["mx"]] = 1
+    for i, well in enumerate(dic["wells"]):
+        for k in range(well[2], well[3] + 1):
+            dic["wellsa"][2 * well[1] + 2 * (dic["nz"] - k - 1) * dic["mx"]] = i
+
+
+def map_xzcoords(dic):
+    """
+    Map the properties from the simulations to the 2D slide
+
+    Args:
+        dic (dict): Global dictionary
+
+    Returns:
+        dic (dict): Modified global dictionary
+
+    """
+    for k in range(dic["nz"]):
+        for i in range(dic["nx"]):
+            ind = i + dic["slide"][1] * dic["nx"] + k * dic["nx"] * dic["ny"]
+            if dic["porv"][ind] > 0:
+                for name in dic["props"]:
+                    dic[name + "a"][2 * i + 2 * (dic["nz"] - k - 1) * dic["mx"]] = dic[
+                        name
+                    ][dic["actind"][ind]]
+                if "porv" in dic["props"]:
+                    dic["porva"][2 * i + 2 * (dic["nz"] - k - 1) * dic["mx"]] = dic[
+                        "porv"
+                    ][ind]
+                if len(dic["wells"]) > 0:
+                    dic["wellsa"][2 * i + 2 * (dic["nz"] - k - 1) * dic["mx"]] = len(
+                        dic["wells"]
+                    )
+                if dic["grid"] == 1:
+                    dic["grida"][2 * i + 2 * (dic["nz"] - k - 1) * dic["mx"]] = 1
+    for i, well in enumerate(dic["wells"]):
+        for k in range(well[2], well[3] + 1):
+            dic["wellsa"][2 * well[0] + 2 * (dic["nz"] - k - 1) * dic["mx"]] = i
+
+
+def map_xycoords(dic):
+    """
+    Map the properties from the simulations to the 2D slide
+
+    Args:
+        dic (dict): Global dictionary
+
+    Returns:
+        dic (dict): Modified global dictionary
+
+    """
+    for j in range(dic["ny"]):
+        for i in range(dic["nx"]):
+            ind = i + j * dic["nx"] + dic["slide"][2] * dic["nx"] * dic["ny"]
+            if dic["porv"][ind] > 0:
+                for name in dic["props"]:
+                    dic[name + "a"][2 * i + 2 * j * dic["mx"]] = dic[name][
+                        dic["actind"][ind]
+                    ]
+                if "porv" in dic["props"]:
+                    dic["porva"][2 * i + 2 * j * dic["mx"]] = dic["porv"][ind]
+                if len(dic["wells"]) > 0:
+                    dic["wellsa"][2 * i + 2 * j * dic["mx"]] = len(dic["wells"])
+                if dic["grid"] == 1:
+                    dic["grida"][2 * i + 2 * j * dic["mx"]] = 1
+    for i, well in enumerate(dic["wells"]):
+        dic["wellsa"][2 * well[0] + 2 * well[1] * dic["mx"]] = i
diff --git a/src/plopm/utils/readers.py b/src/plopm/utils/readers.py
new file mode 100644
index 0000000..316cef7
--- /dev/null
+++ b/src/plopm/utils/readers.py
@@ -0,0 +1,341 @@
+# SPDX-FileCopyrightText: 2024 NORCE
+# SPDX-License-Identifier: GPL-3.0
+
+"""
+Utiliy functions to read the OPM Flow simulator type output files.
+"""
+
+import csv
+import numpy as np
+
+
+def get_yzcoords_resdata(dic):
+    """
+    Handle the coordinates from the OPM Grid to the 2D yz-mesh using resdata
+
+    Args:
+        dic (dict): Global dictionary
+
+    Returns:
+        dic (dict): Modified global dictionary
+
+    """
+    dic["mesh"] = dic["egrid"].export_corners(dic["egrid"].export_index())
+    s_l = dic["slide"][0]
+    for j in range(dic["nz"]):
+        dic["xc"].append([])
+        dic["yc"].append([])
+        for k, n in zip(["x", "x", "y", "y"], [1, 7, 2, 8]):
+            dic[f"{k}c"][-1].append(
+                dic["mesh"][(dic["nz"] - j - 1) * dic["nx"] * dic["ny"] + s_l][n]
+            )
+        for i in range(dic["ny"] - 1):
+            for k, n in zip(["x", "x", "y", "y"], [1, 7, 2, 8]):
+                dic[f"{k}c"][-1].append(
+                    dic["mesh"][
+                        (i + 1) * dic["nx"]
+                        + (dic["nz"] - j - 1) * dic["nx"] * dic["ny"]
+                        + s_l
+                    ][n]
+                )
+        dic["xc"].append([])
+        dic["yc"].append([])
+        for k, n in zip(["x", "x", "y", "y"], [13, 19, 14, 20]):
+            dic[f"{k}c"][-1].append(
+                dic["mesh"][(dic["nz"] - j - 1) * dic["nx"] * dic["ny"] + s_l][n]
+            )
+        for i in range(dic["ny"] - 1):
+            for k, n in zip(["x", "x", "y", "y"], [13, 19, 14, 20]):
+                dic[f"{k}c"][-1].append(
+                    dic["mesh"][
+                        (i + 1) * dic["nx"]
+                        + (dic["nz"] - j - 1) * dic["nx"] * dic["ny"]
+                        + s_l
+                    ][n]
+                )
+
+
+def get_yzcoords_opm(dic):
+    """
+    Handle the coordinates from the OPM Grid to the 2D yz-mesh using opm
+
+    Args:
+        dic (dict): Global dictionary
+
+    Returns:
+        dic (dict): Modified global dictionary
+
+    """
+    for j in range(dic["nz"]):
+        dic["xc"].append([])
+        dic["yc"].append([])
+        for k, c, p in zip(["x", "x", "y", "y"], [1, 1, 2, 2], [0, 2, 0, 2]):
+            dic[f"{k}c"][-1].append(
+                dic["egrid"].xyz_from_ijk(dic["slide"][0], 0, dic["nz"] - j - 1)[c][p]
+            )
+        for i in range(dic["ny"] - 1):
+            for k, c, p in zip(["x", "x", "y", "y"], [1, 1, 2, 2], [0, 2, 0, 2]):
+                dic[f"{k}c"][-1].append(
+                    dic["egrid"].xyz_from_ijk(
+                        dic["slide"][0], i + 1, dic["nz"] - j - 1
+                    )[c][p]
+                )
+        dic["xc"].append([])
+        dic["yc"].append([])
+        for k, c, p in zip(["x", "x", "y", "y"], [1, 1, 2, 2], [4, 6, 4, 6]):
+            dic[f"{k}c"][-1].append(
+                dic["egrid"].xyz_from_ijk(dic["slide"][0], 0, dic["nz"] - j - 1)[c][p]
+            )
+        for i in range(dic["ny"] - 1):
+            for k, c, p in zip(["x", "x", "y", "y"], [1, 1, 2, 2], [4, 6, 4, 6]):
+                dic[f"{k}c"][-1].append(
+                    dic["egrid"].xyz_from_ijk(
+                        dic["slide"][0], i + 1, dic["nz"] - j - 1
+                    )[c][p]
+                )
+
+
+def get_xzcoords_resdata(dic):
+    """
+    Handle the coordinates from the OPM Grid to the 2D xz-mesh using resdata
+
+    Args:
+        dic (dict): Global dictionary
+
+    Returns:
+        dic (dict): Modified global dictionary
+
+    """
+    dic["mesh"] = dic["egrid"].export_corners(dic["egrid"].export_index())
+    s_l = dic["slide"][1] * dic["nx"]
+    for j in range(dic["nz"]):
+        dic["xc"].append([])
+        dic["yc"].append([])
+        for k, n in zip(["x", "x", "y", "y"], [0, 3, 2, 5]):
+            dic[f"{k}c"][-1].append(
+                dic["mesh"][(dic["nz"] - j - 1) * dic["nx"] * dic["ny"] + s_l][n]
+            )
+        for i in range(dic["nx"] - 1):
+            for k, n in zip(["x", "x", "y", "y"], [0, 3, 2, 5]):
+                dic[f"{k}c"][-1].append(
+                    dic["mesh"][
+                        1 + i + (dic["nz"] - j - 1) * dic["nx"] * dic["ny"] + s_l
+                    ][n]
+                )
+        dic["xc"].append([])
+        dic["yc"].append([])
+        for k, n in zip(["x", "x", "y", "y"], [12, 15, 14, 17]):
+            dic[f"{k}c"][-1].append(
+                dic["mesh"][(dic["nz"] - j - 1) * dic["nx"] * dic["ny"] + s_l][n]
+            )
+        for i in range(dic["nx"] - 1):
+            for k, n in zip(["x", "x", "y", "y"], [12, 15, 14, 17]):
+                dic[f"{k}c"][-1].append(
+                    dic["mesh"][
+                        1 + i + (dic["nz"] - j - 1) * dic["nx"] * dic["ny"] + s_l
+                    ][n]
+                )
+
+
+def get_xzcoords_opm(dic):
+    """
+    Handle the coordinates from the OPM Grid to the 2D xz-mesh using opm
+
+    Args:
+        dic (dict): Global dictionary
+
+    Returns:
+        dic (dict): Modified global dictionary
+
+    """
+    for j in range(dic["nz"]):
+        dic["xc"].append([])
+        dic["yc"].append([])
+        for k, c, p in zip(["x", "x", "y", "y"], [0, 0, 2, 2], [0, 1, 0, 1]):
+            dic[f"{k}c"][-1].append(
+                dic["egrid"].xyz_from_ijk(0, dic["slide"][1], dic["nz"] - j - 1)[c][p]
+            )
+        for i in range(dic["nx"] - 1):
+            for k, c, p in zip(["x", "x", "y", "y"], [0, 0, 2, 2], [0, 1, 0, 1]):
+                dic[f"{k}c"][-1].append(
+                    dic["egrid"].xyz_from_ijk(
+                        i + 1, dic["slide"][1], dic["nz"] - j - 1
+                    )[c][p]
+                )
+        dic["xc"].append([])
+        dic["yc"].append([])
+        for k, c, p in zip(["x", "x", "y", "y"], [0, 0, 2, 2], [4, 5, 4, 5]):
+            dic[f"{k}c"][-1].append(
+                dic["egrid"].xyz_from_ijk(0, dic["slide"][1], dic["nz"] - j - 1)[c][p]
+            )
+        for i in range(dic["nx"] - 1):
+            for k, c, p in zip(["x", "x", "y", "y"], [0, 0, 2, 2], [4, 5, 4, 5]):
+                dic[f"{k}c"][-1].append(
+                    dic["egrid"].xyz_from_ijk(
+                        1 + i, dic["slide"][1], dic["nz"] - j - 1
+                    )[c][p]
+                )
+
+
+def get_xycoords_resdata(dic):
+    """
+    Handle the coordinates from the OPM Grid to the 2D xy-mesh using resdata
+
+    Args:
+        dic (dict): Global dictionary
+
+    Returns:
+        dic (dict): Modified global dictionary
+
+    """
+    dic["mesh"] = dic["egrid"].export_corners(dic["egrid"].export_index())
+    s_l = dic["slide"][2] * dic["nx"] * dic["ny"]
+    for j in range(dic["ny"]):
+        dic["xc"].append([])
+        dic["yc"].append([])
+        for k, n in zip(["x", "x", "y", "y"], [0, 3, 1, 4]):
+            dic[f"{k}c"][-1].append(dic["mesh"][j * dic["nx"] + s_l][n])
+        for i in range(dic["nx"] - 1):
+            for k, n in zip(["x", "x", "y", "y"], [0, 3, 1, 4]):
+                dic[f"{k}c"][-1].append(dic["mesh"][1 + i + j * dic["nx"] + s_l][n])
+        dic["xc"].append([])
+        dic["yc"].append([])
+        for k, n in zip(["x", "x", "y", "y"], [6, 9, 7, 10]):
+            dic[f"{k}c"][-1].append(dic["mesh"][j * dic["nx"] + s_l][n])
+        for i in range(dic["nx"] - 1):
+            for k, n in zip(["x", "x", "y", "y"], [6, 9, 7, 10]):
+                dic[f"{k}c"][-1].append(dic["mesh"][1 + i + j * dic["nx"] + s_l][n])
+
+
+def get_xycoords_opm(dic):
+    """
+    Handle the coordinates from the OPM Grid to the 2D xy-mesh using opm
+
+    Args:
+        dic (dict): Global dictionary
+
+    Returns:
+        dic (dict): Modified global dictionary
+
+    """
+    for j in range(dic["ny"]):
+        dic["xc"].append([])
+        dic["yc"].append([])
+        for k, c, p in zip(["x", "x", "y", "y"], [0, 0, 1, 1], [0, 1, 0, 1]):
+            dic[f"{k}c"][-1].append(
+                dic["egrid"].xyz_from_ijk(0, j, dic["slide"][2])[c][p]
+            )
+        for i in range(dic["nx"] - 1):
+            for k, c, p in zip(["x", "x", "y", "y"], [0, 0, 1, 1], [0, 1, 0, 1]):
+                dic[f"{k}c"][-1].append(
+                    dic["egrid"].xyz_from_ijk(i + 1, j, dic["slide"][2])[c][p]
+                )
+        dic["xc"].append([])
+        dic["yc"].append([])
+        for k, c, p in zip(["x", "x", "y", "y"], [0, 0, 1, 1], [2, 3, 2, 3]):
+            dic[f"{k}c"][-1].append(
+                dic["egrid"].xyz_from_ijk(0, j, dic["slide"][2])[c][p]
+            )
+        for i in range(dic["nx"] - 1):
+            for k, c, p in zip(["x", "x", "y", "y"], [0, 0, 1, 1], [2, 3, 2, 3]):
+                dic[f"{k}c"][-1].append(
+                    dic["egrid"].xyz_from_ijk(i + 1, j, dic["slide"][2])[c][p]
+                )
+
+
+def get_kws(dic):
+    """
+    Get the static properties from the INIT file using ResdataFile
+
+    Args:
+        dic (dict): Global dictionary
+
+    Returns:
+        dic (dict): Modified global dictionary
+
+    """
+    for name in dic["props"]:
+        if dic["use"] == "resdata":
+            if dic["init"].has_kw(name.upper()):
+                dic[name] = np.array(dic["init"].iget_kw(name.upper())[0])
+            elif dic["unrst"].has_kw(name.upper()):
+                dic[name] = np.array(
+                    dic["unrst"].iget_kw(name.upper())[dic["restart"] - 1]
+                )
+                ntot = len(dic["unrst"].iget_kw(name.upper()))
+                dic["dtitle"] = (
+                    f", rst {ntot if dic['restart']==0 else dic['restart']} out of {ntot}"
+                )
+            elif dic["summary"].has_key(name.upper()):
+                dic["vsum"] = dic["summary"][name.upper()].values
+                dic["time"] = dic["summary"]["TIME"].values
+                return
+        else:
+            if dic["init"].count(name.upper()):
+                dic[name] = np.array(dic["init"][name.upper()])
+            elif dic["unrst"].count(name.upper()):
+                nrst = (
+                    dic["restart"] - 1
+                    if dic["restart"] > 0
+                    else dic["unrst"].count(name.upper()) - 1
+                )
+                dic[name] = np.array(dic["unrst"][name.upper(), nrst])
+                ntot = dic["unrst"].count(name.upper())
+                dic["dtitle"] = (
+                    f", rst {ntot if dic['restart']==0 else dic['restart']} out of {ntot}"
+                )
+            elif name.upper() in dic["summary"].keys():
+                dic["vsum"] = dic["summary"][name.upper()]
+                dic["time"] = dic["summary"]["TIME"]
+                return
+        dic[name + "a"] = np.ones(dic["mx"] * dic["my"]) * np.nan
+    if dic["use"] == "opm":
+        dic["porv"] = np.array(dic["init"]["PORV"])
+    else:
+        dic["porv"] = np.array(dic["init"].iget_kw("PORV")[0])
+    dic["porva"] = np.ones(dic["mx"] * dic["my"]) * np.nan
+    dic["actind"] = np.cumsum([1 if porv > 0 else 0 for porv in dic["porv"]]) - 1
+
+
+def get_wells(dic):
+    """
+    Using the input deck (.DATA) to read the i,j well locations
+
+    Args:
+        dic (dict): Global dictionary
+
+    Returns:
+        dic (dict): Modified global dictionary
+
+    """
+    dic["wellsa"] = np.ones((dic["mx"]) * (dic["my"])) * np.nan
+    wells, sources = False, False
+    with open(f"{dic['name']}.DATA", "r", encoding="utf8") as file:
+        for row in csv.reader(file):
+            nrwo = str(row)[2:-2]
+            if wells:
+                if len(nrwo.split()) < 2:
+                    break
+                dic["wells"].append(
+                    [
+                        int(nrwo.split()[1]) - 1,
+                        int(nrwo.split()[2]) - 1,
+                        int(nrwo.split()[3]) - 1,
+                        int(nrwo.split()[4]) - 1,
+                    ]
+                )
+            if sources:
+                if len(nrwo.split()) < 2:
+                    break
+                dic["wells"].append(
+                    [
+                        int(nrwo.split()[0]) - 1,
+                        int(nrwo.split()[1]) - 1,
+                        int(nrwo.split()[2]) - 1,
+                        int(nrwo.split()[2]) - 1,
+                    ]
+                )
+            if nrwo == "COMPDAT":
+                wells = True
+            if nrwo == "SOURCE":
+                sources = True
diff --git a/src/plopm/utils/write.py b/src/plopm/utils/write.py
new file mode 100644
index 0000000..45d05bb
--- /dev/null
+++ b/src/plopm/utils/write.py
@@ -0,0 +1,230 @@
+# SPDX-FileCopyrightText: 2024 NORCE
+# SPDX-License-Identifier: GPL-3.0
+
+"""
+Utiliy functions to write the PNGs figures.
+"""
+
+import numpy as np
+import matplotlib
+import matplotlib.pyplot as plt
+from mpl_toolkits.axes_grid1 import make_axes_locatable
+
+
+def make_summary(dic):
+    """
+    Plot the summary variable
+
+    Args:
+        cmdargs (dict): Command arguments
+
+    Returns:
+        None
+
+    """
+    plt.rcParams.update({"axes.grid": True})
+    fig, axis = plt.subplots()
+    if len(dic["cmaps"]) == 1:
+        axis.step(dic["time"], dic["vsum"], color=dic["cmaps"][0])
+    else:
+        axis.step(dic["time"], dic["vsum"], color="b")
+    axis.set_ylabel(dic["props"][0])
+    axis.set_xlabel("Time [s]")
+    axis.set_xlim([min(dic["time"]), max(dic["time"])])
+    axis.set_ylim([min(dic["vsum"]), max(dic["vsum"])])
+    axis.set_xticks(
+        np.linspace(
+            min(dic["time"]),
+            max(dic["time"]),
+            4,
+        )
+    )
+    axis.set_yticks(
+        np.linspace(
+            min(dic["vsum"]),
+            max(dic["vsum"]),
+            4,
+        )
+    )
+    fig.savefig(f"{dic['output']}/{dic['variable']}.png", bbox_inches="tight", dpi=300)
+    plt.close()
+
+
+def make_2dmaps(dic):
+    """
+    Method to create the 2d maps using pcolormesh
+
+    Args:
+        dic (dict): Global dictionary
+
+    Returns:
+        dic (dict): Modified global dictionary
+
+    """
+    maps = np.ones([dic["my"], dic["mx"]]) * np.nan
+    if dic["well"] == 1:
+        dic["props"] += ["wells"]
+    elif dic["grid"] == 1:
+        dic["props"] += ["grid"]
+    for n, name in enumerate(dic["props"]):
+        fig, axis = plt.subplots()
+        for i in np.arange(0, dic["my"]):
+            maps[i, :] = dic[name + "a"][i * (dic["mx"]) : (i + 1) * (dic["mx"])]
+        imag = axis.pcolormesh(
+            dic["xc"],
+            dic["yc"],
+            maps,
+            shading="flat",
+            cmap=dic["cmaps"][n],
+        )
+        ntick = 5
+        if dic["well"] != 1 and dic["grid"] != 1:
+            minc = dic[name][~np.isnan(dic[name])].min()
+            maxc = dic[name][~np.isnan(dic[name])].max()
+            if maxc == minc:
+                ntick = 1
+            elif name in ["fipnum", "satnum"]:
+                ntick = maxc
+        else:
+            minc = 0
+            maxc = len(dic["wells"])
+        divider = make_axes_locatable(axis)
+        vect = np.linspace(
+            minc,
+            maxc,
+            ntick,
+            endpoint=True,
+        )
+        if dic["well"] != 1 and dic["grid"] != 1:
+            fig.colorbar(
+                imag,
+                cax=divider.append_axes("right", size="5%", pad=0.05),
+                orientation="vertical",
+                ticks=vect,
+                label=dic["units"][n],
+                format=dic["format"][n],
+            )
+        else:
+            handle_well_or_grid(dic, fig, imag, divider, vect)
+        imag.set_clim(
+            minc,
+            maxc,
+        )
+        handle_axis(dic, axis, name)
+        fig.savefig(
+            f"{dic['output']}/{name}_{dic['nslide']}.png", bbox_inches="tight", dpi=300
+        )
+        plt.close()
+
+
+def handle_axis(dic, axis, name):
+    """
+    Method to handle the figure axis
+
+    Args:
+        dic (dict): Global dictionary\n
+        axis (class): Axis object\n
+        name (str): Property to plot
+
+    Returns:
+        axis (class): Modified axis object
+
+    """
+    if dic["scale"] == "yes":
+        axis.axis("scaled")
+    axis.set_xlabel(f"{dic['xmeaning']} [m]")
+    axis.set_ylabel(f"{dic['ymeaning']} [m]")
+    extra = ""
+    if name == "porv":
+        extra = f" (sum={sum(dic[name]):.3e})"
+    axis.set_title(name + dic["tslide"] + dic["dtitle"] + extra)
+    if name == "grid":
+        axis.set_title(
+            f"Grid = [{dic['nx']},{dic['ny']},{dic['nz']}], "
+            + f"Total no. active cells = {max(dic['actind'])+1}"
+        )
+    if dic["slide"][2] == -2:
+        axis.invert_yaxis()
+    axis.set_xticks(
+        np.linspace(
+            min(min(dic["xc"])),
+            max(max(dic["xc"])),
+            4,
+        )
+    )
+    axis.set_yticks(
+        np.linspace(
+            min(min(dic["yc"])),
+            max(max(dic["yc"])),
+            4,
+        )
+    )
+    if len(dic["xlim"]) > 1:
+        axis.set_xlim(dic["xlim"])
+    if len(dic["ylim"]) > 1:
+        axis.set_ylim(dic["ylim"])
+
+
+def handle_well_or_grid(dic, fig, imag, divider, vect):
+    """
+    Method to create the 2d maps using pcolormesh
+
+    Args:
+        dic (dict): Global dictionary\n
+        fig (class): Figure object\n
+        imag (class): Actual plot object\n
+        divider (class): Object for the color bar axis\n
+        vect (array): Floats for the labels in the color bar
+
+    Returns:
+        fig (class): Modified figure object\n
+        plt (class): Modified plotting object\n
+
+    """
+    fig.colorbar(
+        imag,
+        cax=divider.append_axes("right", size="0%", pad=0.05),
+        orientation="vertical",
+        ticks=vect,
+        format=lambda x, _: "",
+    )
+    cmap = matplotlib.colormaps["nipy_spectral"]
+    colors = cmap(np.linspace(0, 1, len(dic["wells"]) + 1))
+    if len(dic["wells"]) < 20:
+        for i, well in enumerate(dic["wells"]):
+            if well[2] != well[3]:
+                plt.text(
+                    0,
+                    i,
+                    f"{i}-({well[0]+1},{well[1]+1},{well[2]+1}-{well[3]+1})",
+                    c=colors[i],
+                    fontweight="bold",
+                )
+            else:
+                plt.text(
+                    0,
+                    i,
+                    f"{i}-({well[0]+1},{well[1]+1},{well[2]+1})",
+                    c=colors[i],
+                    fontweight="bold",
+                )
+    else:
+        for i, well in zip(
+            [0, len(dic["wells"]) - 1], [dic["wells"][0], dic["wells"][-1]]
+        ):
+            if well[2] != well[3]:
+                plt.text(
+                    0,
+                    i,
+                    f"{i}-({well[0]+1},{well[1]+1},{well[2]+1}-{well[3]+1})",
+                    c=colors[i],
+                    fontweight="bold",
+                )
+            else:
+                plt.text(
+                    0,
+                    i,
+                    f"{i}-({well[0]+1},{well[1]+1},{well[2]+1})",
+                    c=colors[i],
+                    fontweight="bold",
+                )
diff --git a/tests/test_wells.py b/tests/test_wells.py
index 1f52deb..022cff8 100644
--- a/tests/test_wells.py
+++ b/tests/test_wells.py
@@ -12,5 +12,5 @@ def test_wells():
         ["plopm", "-i", "SPE10_MODEL2", "-o", ".", "-w", "1"],
         check=True,
     )
-    assert os.path.exists(f"{cwd}/tests/generic_deck/wells.png")
+    assert os.path.exists(f"{cwd}/tests/generic_deck/wells_*,1,*.png")
     os.chdir(cwd)