diff --git a/docs/index.rst b/docs/index.rst
index c6be79b5..f5d2ed56 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -11,6 +11,7 @@ Contents:
    qcforward
    qcproperties
    qcreset
+   properties
    create_rft_ertobs
    rename_rms_scripts
    rms_import_localmodule
diff --git a/docs/pdf/sw_calc.pdf b/docs/pdf/sw_calc.pdf
new file mode 100755
index 00000000..561d4f0f
Binary files /dev/null and b/docs/pdf/sw_calc.pdf differ
diff --git a/docs/properties.rst b/docs/properties.rst
new file mode 100644
index 00000000..bf2ca3df
--- /dev/null
+++ b/docs/properties.rst
@@ -0,0 +1,79 @@
+The properties package
+======================
+
+The fmu.tools ``properties`` package provides functionality (methods) to special calculations
+in e.g. RMS. Currently a library for computing water saturation is present.
+
+Generic water saturation calculations
+-------------------------------------
+
+Water saturation in RMS user interface is limited to Leverett J functions and lookup functions.
+However there is a need to extend this to Leverett J function which has an additional constant term,
+and to BVW and Brooks-Corey methods. The function offers several options:
+
+* Using direct methods (cell mid-point) or integrate across cell thickness
+* Look at the center line (vertically) through cells, og use cells corners.
+* Both normal and inverse formulation of Leverett J is supported.
+* Normalization is supported.
+
+The theory is given in the PDF file provided here: :download:`pdf <pdf/sw_calc.pdf>`.
+
+Note that this library will not give exact same values as RMS built-in Sw GUI will give.
+This is partly due to computing height across a cell is a bit different implemented.
+
+The library input is documented here: :meth:`.SwFunction`
+
+
+Example: Using simple Sw Leverett J
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+The Leverett J is given as $S_w = a J^b$, where $J$ is $h\sqrt(k/\phi)$. In this example, the
+factors are constant numbers:
+
+.. code-block:: python
+
+    import numpy as np
+    import xtgeo
+    from fmu.tools.properties import SwFunction
+
+    GRIDNAME = "Geogrid"
+    A = 1
+    B = -2
+    FWL = 1700
+
+    PORONAME = "PHIT"
+    PERMNAME = "KLOGH"
+
+    CALC = "integrated"
+    METHOD = "cell_corners_above_ffl"
+
+    SW_RESULT = "SwJ"
+
+
+    def compute_saturation():
+        grid = xtgeo.grid_from_roxar(project, GRIDNAME)
+        poro = xtgeo.gridproperty_from_roxar(project, GRIDNAME, PORONAME)
+        perm = xtgeo.gridproperty_from_roxar(project, GRIDNAME, PERMNAME)
+
+        # avoid zero porosities in the division below
+        poro.values[poro.values == 0.0] = 0.01
+        x = poro.copy()
+        x.values = np.sqrt(np.divide(perm.values, poro.values))
+
+        sw_func = SwFunction(
+            grid=grid,
+            x=x,
+            a=A,
+            b=B,
+            ffl=FWL,
+            invert=True,
+            method=METHOD,
+        )
+        sw = sw_func.compute(CALC)
+        sw.values[poro.values < 0.05] = 1.0  # 100% water when low porosity
+        sw.to_roxar(project, GRIDNAME, SW_RESULT)
+
+
+    if __name__ == "__main__":
+        compute_saturation()
+        print("Done")
diff --git a/src/fmu/tools/_common.py b/src/fmu/tools/_common.py
index e3414ee8..6c56a281 100644
--- a/src/fmu/tools/_common.py
+++ b/src/fmu/tools/_common.py
@@ -2,9 +2,9 @@
 This private module is used for some common functions and testing
 """
 
+import os
 import sys
 from functools import wraps
-import os
 
 
 class _QCCommon(object):
diff --git a/src/fmu/tools/properties/__init__.py b/src/fmu/tools/properties/__init__.py
new file mode 100644
index 00000000..90bd6d6b
--- /dev/null
+++ b/src/fmu/tools/properties/__init__.py
@@ -0,0 +1,5 @@
+from .swfunction import SwFunction
+
+__all__ = [
+    "SwFunction",
+]
diff --git a/src/fmu/tools/properties/swfunction.py b/src/fmu/tools/properties/swfunction.py
new file mode 100644
index 00000000..294554f0
--- /dev/null
+++ b/src/fmu/tools/properties/swfunction.py
@@ -0,0 +1,324 @@
+"""Sw J or BVW or Brooks-Corey general functions, by direct or integration.
+
+Solve a system on form:
+
+sw = a(m + x*h)^b
+
+The "x" can e.g. be the J term in the Leverett equation, but without the height term.
+
+In many cases m is zero.
+
+Note:
+The public library `xtgeo` must be available.
+"""
+
+from __future__ import annotations
+
+import logging
+from dataclasses import dataclass, field
+from typing import Any
+
+import numpy as np
+import xtgeo
+
+# custom logger for this module
+logger = logging.getLogger(__name__)
+
+ALLOWED_METHODS = [
+    "cell_center",
+    "cell_center_above_ffl",
+    "cell_corners_above_ffl",
+]
+
+
+@dataclass
+class SwFunction:
+    """Generic Sw calc for formulation Sw = a * (m + x * h)^b.
+
+    Some theory for this is shown here: :download:`pdf <pdf/sw_calc.pdf>`.
+
+    Args:
+        grid: The xtgeo grid object
+        x: The ``x`` term in the generic saturation. For a simplified J function, this
+            will be sqrt(perm/poro).
+        a: The ``a`` term in the generic saturation.
+        b: The ``b`` term in the generic saturation.
+        ffl: The free fluid level, either as a number or a 3D xtgeo gridproperty
+            object.
+        m: The ``m`` term in the equation above, either as a number or a 3D xtgeo
+            gridproperty object. Defaults to zero.
+        swira: The absolute irreducable satuation applied as asymptote in the equation.
+            Defaults to zero.
+        swmax: The maximum applied as asymptote in the equation.
+            Defaults to 1.
+        method: How to look at cell geometry. There are three options: 'cell_center',
+            'cell_center-above_ffl' and 'cell_corners_above_ffl'.
+        gridname: Optional - only needed when debug keyword is True.
+        project: Optional. Must be included when working in RMS.
+        invert: Optional, will in case invert ``a`` and ``b`` in case the J function
+            (or similar) is formulated as Sw = (J/A)^(1/B). This is the case in e.g.
+            RMS
+        zdepth: Optional. May speed up computation if provided, in case the function
+            is called several time for the same grid.
+        hcenter: Optional. May speed up computation if provided, in case the function
+            is called several time for the same grid, or hcenter is pre-computed in a
+            special way. The hcenter will work with direct method only, not integration.
+        htop: Optional. May speed up computation if provided, in case the function
+            is called several time for the same grid.
+        hbot: Optional. May speed up computation if provided, in case the function
+            is called several time for the same grid.
+        debug: If True, several "check parameters" will we given, either as grid
+            properties (if working in RMS) or as ROFF files (if working outside RMS).
+        tag: Optional string to identify debug parameters.
+
+
+    """
+
+    grid: xtgeo.Grid
+    # all input are constants or xtgeo grid or gridproperty
+    x: float | xtgeo.GridProperty = 0.0
+    a: float | xtgeo.GridProperty = 1.0
+    b: float | xtgeo.GridProperty = -1.0
+    ffl: float | xtgeo.GridProperty = 999.0
+
+    m: float | xtgeo.GridProperty = 0.0
+
+    swira: float | xtgeo.GridProperty = 0.0
+    swmax: float | xtgeo.GridProperty = 1.0
+
+    method: str = "cell_center_above_ffl"
+    gridname: str = ""  # only used when debugging in RMS
+    project: Any | None = None
+    invert: bool = False  # if the SwJ function is on "reverse" form
+
+    # if None they will be computed here; otherwise they can be given explicitly:
+    hcenter: xtgeo.GridProperty = None
+    htop: xtgeo.GridProperty = None
+    hbot: xtgeo.GridProperty = None
+
+    # debug flag for additional parameters. Given in RMS if project is given; otherwise
+    # as files on your working folder
+    debug: bool = False
+    tag: str = ""  # identification tag to add to debug params
+
+    # derived and internal
+    _swval: np.ma.MaskedArray = field(init=False)
+    _sw: Any = field(init=False)
+
+    def __post_init__(self) -> None:
+        """Post init."""
+
+        if self.tag and not self.tag.startswith("_"):
+            self.tag = "_" + self.tag
+
+        self._process_input()
+
+        if all([self.hcenter, self.htop, self.hbot]):
+            logger.info("Essential geometries are pre-computed")
+        else:
+            self._compute_htop_hbot()
+
+        if self.debug:
+            if self.project:
+                self.grid.to_roxar(self.project, "DEBUG_" + self.gridname)
+            else:
+                self.grid.to_file("debug_grid.roff")
+
+    def _process_input(self) -> None:
+        """Work with a, b, x, ie. inversing and convert from float."""
+        logger.info("Process a, b, etc")
+
+        if self.method not in ALLOWED_METHODS:
+            raise ValueError(
+                f"The method <{self.method}> is invalid. Allowed: {ALLOWED_METHODS}"
+            )
+
+        if isinstance(self.a, (int, float)):
+            self.a = xtgeo.GridProperty(self.grid, values=self.a)
+        if isinstance(self.b, (int, float)):
+            self.b = xtgeo.GridProperty(self.grid, values=self.b)
+        if isinstance(self.x, (int, float)):
+            self.x = xtgeo.GridProperty(self.grid, values=self.x)
+        if isinstance(self.m, (int, float)):
+            self.m = xtgeo.GridProperty(self.grid, values=self.m)
+        if isinstance(self.ffl, (int, float)):
+            self.ffl = xtgeo.GridProperty(self.grid, values=self.ffl)
+        if isinstance(self.swira, (int, float)):
+            self.swira = xtgeo.GridProperty(self.grid, values=self.swira)
+        if isinstance(self.swmax, (int, float)):
+            self.swmax = xtgeo.GridProperty(self.grid, values=self.swmax)
+
+        if self.invert:
+            orig_a = self.a.values.mean()
+            orig_b = self.b.values.mean()
+            self.b.values = 1.0 / self.b.values
+            self.a.values = np.ma.power(1.0 / self.a.values, self.b.values)
+            new_a = self.a.values.mean()
+            new_b = self.b.values.mean()
+            logger.info(
+                "Inverted a and b; original average was %s %s, new is %s %s",
+                orig_a,
+                orig_b,
+                new_a,
+                new_b,
+            )
+
+    def _compute_htop_hbot(self) -> None:
+        """Setting geometries for 'bottom' and 'top' cell.
+
+        All these values are relative to the given contact, as "height above". For
+        cells under the contact, zero or a negative number will be applied.
+        """
+        assert isinstance(self.ffl, xtgeo.GridProperty)  # mypy
+
+        htop, hbot, hmid = self.grid.get_heights_above_ffl(self.ffl, option=self.method)
+
+        self.htop = htop
+        self.hbot = hbot
+        self.hcenter = hmid
+        logger.info("Use method %s", self.method)
+
+        if self.debug:
+            htop.name = "TOP_SW" + self.tag
+            hbot.name = "BOT_SW" + self.tag
+            hmid.name = "CENTER_SW" + self.tag
+            if self.project:
+                logger.debug("TAG is ", self.tag)
+                htop.to_roxar(self.project, self.gridname, "HTOP" + self.tag)
+                hbot.to_roxar(self.project, self.gridname, "HBOT" + self.tag)
+                hmid.to_roxar(self.project, self.gridname, "HCENTER" + self.tag)
+            else:
+                htop.to_file(f"debug_htop{self.tag}.roff")
+                hbot.to_file(f"debug_hbot{self.tag}.roff")
+                hmid.to_file(f"debug_hcenter{self.tag}.roff")
+
+    def _sw_function_direct(self) -> None:
+        """Use function on form Sw = A*(M + X*h)^B; generic function!"""
+        assert isinstance(self.a, xtgeo.GridProperty)  # mypy
+        assert isinstance(self.b, xtgeo.GridProperty)  # mypy
+        assert isinstance(self.x, xtgeo.GridProperty)  # mypy
+        assert isinstance(self.m, xtgeo.GridProperty)  # mypy
+        assert isinstance(self.swira, xtgeo.GridProperty)  # mypy
+        assert isinstance(self.swmax, xtgeo.GridProperty)  # mypy
+
+        # the direct function is mostly used to compare with integrated approach, as QC
+        height = self.hcenter.values
+        xh_val = self.x.values * height
+
+        self._swval = self.a.values * np.ma.power(self.m.values + xh_val, self.b.values)
+
+        # normalize before any limitations
+        self._swval = (
+            self.swira.values + (self.swmax.values - self.swira.values) * self._swval
+        )
+
+        self._swval = np.ma.where(height <= 0.0, 1, self._swval)
+        self._swval = np.ma.where(self._swval > 1.0, 1.0, self._swval)
+
+        # store result
+        self._sw = xtgeo.GridProperty(self.grid, values=self._swval)
+
+    def _sw_function_integrate_w_mterm(self) -> None:
+        """Integrate a function on form Sw = A*(M + JX*H)^B; generic function!
+
+        It assumes here:
+        * A: self.a
+        * B: self.b
+        * X: self.x
+        * M: self.m
+        * self.ffl is the free fluid level for the phase under investigation
+        """
+        assert isinstance(self.a, xtgeo.GridProperty)  # mypy
+        assert isinstance(self.b, xtgeo.GridProperty)  # mypy
+        assert isinstance(self.x, xtgeo.GridProperty)  # mypy
+        assert isinstance(self.m, xtgeo.GridProperty)  # mypy
+        assert isinstance(self.swira, xtgeo.GridProperty)  # mypy
+        assert isinstance(self.swmax, xtgeo.GridProperty)  # mypy
+        ht = (
+            ((1.0 / self.a.values) ** (1.0 / self.b.values)) - self.m.values
+        ) / self.x.values  # threshold height
+
+        h2 = self.htop.values.copy()  # h_top or H2 in integration
+        h1 = self.hbot.values.copy()  # h_bot or H1 in integration
+
+        if self.debug:
+            tmp = xtgeo.GridProperty(self.grid, values=ht, name="HT" + self.tag)
+            if self.project:
+                tmp.to_roxar(self.project, self.gridname, "THRESHOLD HEIGHT" + self.tag)
+            else:
+                tmp.to_file(f"debug_ht{self.tag}.roff")
+
+        water = h2 * 0.0
+        water = np.ma.where(h2 < ht, 1.0, water)
+        water = np.ma.where(h2 <= 0.0, 1.0, water)  # may occur for negative ht
+
+        # get all possible corner cases when a cell is close to free fluid level
+        # some of these are probably hitting twice or more; the point is to integrate
+        # over correct interval
+        h1x = np.ma.where(h1 <= 0.0, 0.0, h1)
+        h2x = np.ma.where(h2 <= h1x, h1x, h2)
+        htx = ht - h1x
+        htx = np.ma.where(htx <= 0.0, 0.0, htx)
+        htx = np.ma.where(h2x <= 0.0, 0.0, htx)
+        hxx = np.ma.where(ht > h1x, ht, h1x)
+
+        dh = h2x - h1x
+        dh = np.ma.where(dh < 0.001, 0.001, dh)  # avoid zero division
+
+        self._swval = (
+            htx
+            + self.a.values
+            / (self.x.values * (self.b.values + 1))
+            * (
+                np.ma.power(self.m.values + self.x.values * h2x, self.b.values + 1)
+                - np.ma.power(self.m.values + self.x.values * hxx, self.b.values + 1)
+            )
+        ) / dh
+
+        # mark collapsed cells, where direct calulation is needed
+        mark = (h2x - h1x) * 0.0
+        mark = np.ma.where(dh <= 0.001, 1.0, mark)
+        jval = self.x.values * h2x
+        coll_swval = self.a.values * np.ma.power(self.m.values + jval, self.b.values)
+        self._swval = np.ma.where(mark > 0.5, coll_swval, self._swval)
+
+        # normalize before any limitations
+        self._swval = (
+            self.swira.values + (self.swmax.values - self.swira.values) * self._swval
+        )
+
+        # limit
+        self._swval = np.ma.where(water > 0.5, 1.0, self._swval)
+        self._swval = np.ma.where(self._swval > 1.0, 1.0, self._swval)
+        self._swval = np.ma.where(self._swval < 0.00099, 0.00099, self._swval)
+
+        self._sw = xtgeo.GridProperty(self.grid, values=self._swval)
+
+    def _compute_integrated(self) -> None:
+        """Compute Sw, integrated version."""
+        logger.info("Integration, do integration over height of cell...")
+        self._sw_function_integrate_w_mterm()
+
+        if self._sw.values.max() > 1.0:
+            raise RuntimeError(f"SW max out of range: {self._sw.values.max()}")
+        if self._sw.values.min() < 0.0:
+            raise RuntimeError(f"SW min out of range: {self._sw.values.min()}")
+
+    def _compute_direct(self) -> None:
+        """Compute Sw, direct version (less precise wiht thick cells)."""
+        logger.info("Direct, no integration over height...")
+        self._sw_function_direct()
+
+        if self._sw.values.max() > 1.0:
+            raise RuntimeError(f"SW max out of range: {self._sw.values.max()}")
+        if self._sw.values.min() < 0.0:
+            raise RuntimeError(f"SW min out of range: {self._sw.values.min()}")
+
+    def compute(self, compute_method: str = "integrated") -> xtgeo.GridProperty:
+        """Common compute function for saturation, and returns the Sw property"""
+        if compute_method == "integrated":
+            self._compute_integrated()
+        else:
+            self._compute_direct()
+
+        return self._sw
diff --git a/tests/properties/test_swfunction.py b/tests/properties/test_swfunction.py
new file mode 100644
index 00000000..2eb8293d
--- /dev/null
+++ b/tests/properties/test_swfunction.py
@@ -0,0 +1,303 @@
+"""Testing properties Sw"""
+
+# for tests vs RMS cf /private/jriv/work/testing/swfunc/test_swfunc.rms13.1.2
+import math
+from pathlib import Path
+
+import numpy as np
+import pytest
+import xtgeo
+from fmu.tools.properties import SwFunction
+
+TPATH = Path("../xtgeo-testdata")
+
+GRIDDATA1 = TPATH / "3dgrids/reek/reek_sim_grid.roff"
+PORODATA1 = TPATH / "3dgrids/reek/reek_sim_poro.roff"
+PERMDATA1 = TPATH / "3dgrids/reek/reek_sim_permx.roff"
+
+
+@pytest.mark.parametrize(
+    "avalue, bvalue, ffl, direct, cellmethod, expected_mean",
+    [
+        (1, -2, 13, True, "cell_center_above_ffl", 0.408742),
+        (1, -2, 13, False, "cell_center_above_ffl", 0.412536),
+        (0.3, -5, 12.5, True, "cell_center_above_ffl", 0.547251),
+        (0.3, -5, 12.5, False, "cell_center_above_ffl", 0.549613),
+        (0.3, -5, 12.5, False, "cell_corners_above_ffl", 0.549613),
+    ],
+)
+def test_swj_simple(avalue, bvalue, ffl, direct, cellmethod, expected_mean):
+    """Test a simple SwJ setup, expected mean are checked with RMS Sw job"""
+
+    grid1 = xtgeo.create_box_grid((10, 10, 20), increment=(1, 1, 1), origin=(0, 0, 0))
+    poro = 0.3
+    perm = 300
+
+    xvalue = math.sqrt(perm / poro)
+
+    sw_obj = SwFunction(
+        grid=grid1,
+        a=avalue,
+        b=bvalue,
+        x=xvalue,
+        ffl=ffl,
+        invert=True,
+        method=cellmethod,
+    )
+    sw = sw_obj.compute("direct" if direct else "integrated")
+    assert sw.values.mean() == pytest.approx(expected_mean, rel=0.01)
+
+
+def test_swj_simple_threshold_2grids():
+    """Test a simple SwJ setup, expected mean are checked with RMS Sw job.
+
+    In this case, the threshold height will be approx 5.42 meters. Values for assertion
+    are checked in RMS.
+    """
+
+    grid1 = xtgeo.create_box_grid((10, 10, 20), increment=(1, 1, 1), origin=(0, 0, 0))
+    grid2 = xtgeo.create_box_grid(
+        (10, 10, 200), increment=(1, 1, 0.1), origin=(12, 0, 0)
+    )
+    poro = 0.3
+    perm = 1.0
+    avalue = 10.0
+    bvalue = -5.0
+    ffl = 12.5
+    cellmethod = "cell_corners_above_ffl"
+
+    xvalue = math.sqrt(perm / poro)
+
+    sw_obj = SwFunction(
+        grid=grid1,
+        a=avalue,
+        b=bvalue,
+        x=xvalue,
+        ffl=ffl,
+        invert=True,
+        method=cellmethod,
+    )
+    sw = sw_obj.compute("integrated")
+    assert sw.values.mean() == pytest.approx(0.9689, rel=0.01)
+
+    sw_obj = SwFunction(
+        grid=grid2,
+        a=avalue,
+        b=bvalue,
+        x=xvalue,
+        ffl=ffl,
+        invert=True,
+        method=cellmethod,
+    )
+    sw = sw_obj.compute("integrated")
+
+    assert sw.values.mean() == pytest.approx(0.9689, rel=0.01)
+    assert float(sw.values[0, 0, 70]) == pytest.approx(1, rel=0.0001)
+
+
+@pytest.mark.parametrize(
+    "direct, cellmethod, expected_mean, exp_cell1",
+    [
+        (True, "cell_center_above_ffl", 0.7057, 0.046719),  # n/a vs RMS
+        (False, "cell_center_above_ffl", 0.70736, 0.046724),  # n/a vs RMS
+        (False, "cell_corners_above_ffl", 0.674485, 0.046791),  # n/a vs RMS
+    ],
+)
+def test_swj_simple_reek(direct, cellmethod, expected_mean, exp_cell1):
+    """Test a simple SwJ setup, expected mean are checked with RMS Sw job"""
+
+    avalue = 1
+    bvalue = -0.5
+    ffl = 1700
+
+    grid = xtgeo.grid_from_file(GRIDDATA1)
+    poro = xtgeo.gridproperty_from_file(PORODATA1)
+    perm = xtgeo.gridproperty_from_file(PERMDATA1)
+
+    xval = perm.copy()
+
+    xval.values = np.ma.sqrt(perm.values / poro.values)
+
+    sw_obj = SwFunction(
+        grid=grid,
+        a=avalue,
+        b=bvalue,
+        x=xval,
+        ffl=ffl,
+        invert=False,
+        method=cellmethod,
+    )
+
+    sw = sw_obj.compute("direct" if direct else "integrated")
+
+    assert sw.values.mean() == pytest.approx(expected_mean, rel=0.01)
+
+    assert float(sw.values[32, 35, 6]) == pytest.approx(exp_cell1, abs=0.0001)
+
+
+def test_sw_bvw():
+    """Test a simple BVW setup, expected values are checked from spreadsheet.
+
+    In BVW, the Sw = A * P^B * poro^C
+
+    Here A, B, C can be derived from regressions
+
+    A = a1*poro + a2
+    B = b1*poro + b2
+    C = c1*poro + c2  ...or...  c1*poro + c2 - 1
+
+    The P kan be normalized capillary pressure, and for gas above an oil zone with
+    thickness H_oil, this will be
+
+    P = Pcn = H_oil * (gradw - grado)/ado  +  h * (gradw - gradg)/adg
+
+    """
+    poro = 0.3
+    a1 = -0.5561
+    a2 = 0.3385
+    b1 = 0.01985
+    b2 = -0.3953
+    c1 = 5.648
+    c2 = 0.4987
+
+    rw = 0.1002
+    rg = 0.0113
+    ro = 0.0785
+
+    adg = 42.7
+    ado = 26.0
+    h_oil = 3
+
+    ffl = 30.5  # FOL in this example
+
+    grid = xtgeo.create_box_grid((1, 1, 200))
+
+    a_term = a1 * poro + a2
+    b_term = b1 * poro + b2
+    c_term = c1 * poro + c2 - 1.0
+
+    assert a_term == pytest.approx(0.171670)
+    assert b_term == pytest.approx(-0.389345)
+    assert c_term == pytest.approx(2.193100 - 1.0)
+    # was done in input
+
+    avalue = a_term * poro**c_term
+    bvalue = b_term
+    mvalue = h_oil * (rw - ro) / ado
+    xvalue = (rw - rg) / adg
+
+    # PCN at 10m above FFL:
+
+    # in practice, this will always be the first term(?)
+    m10 = min([h_oil * (rw - ro) / ado, max([0, 10 + h_oil * (rw - ro) / ado])])
+
+    # in practice, this will always be the second term(?)
+    x10 = max([0, 10 * (rw - rg) / adg])
+    assert m10 + x10 == pytest.approx(0.0233235)
+
+    manual10 = a_term * ((m10 + x10) ** b_term) * (poro**c_term)
+    assert manual10 == pytest.approx(0.176335192)
+
+    sw_obj = SwFunction(
+        grid=grid,
+        a=avalue,
+        b=bvalue,
+        m=mvalue,
+        x=xvalue,
+        ffl=ffl,
+        method="cell_center_above_ffl",
+        debug=False,
+    )
+    sw = sw_obj.compute("direct")
+
+    sw10 = float(sw.values[:, :, 20])  # 10 meter above contact
+    assert sw10 == pytest.approx(manual10)
+
+    sw20 = float(sw.values[:, :, 10])  # 20 meter above contact
+    assert sw20 == pytest.approx(0.13755086)
+
+    sw = sw_obj.compute("integrated")
+    sw10_i = float(sw.values[:, :, 20])
+    assert sw10_i == pytest.approx(sw10, abs=0.0001)
+
+
+def test_sw_brooks_corey():
+    """Test a simple Brooks-Corey setup, expected values are checked from spreadsheet.
+
+    In BVW, the Sw = (PCNe / Pcn)^(1/N)
+
+    Sw,final = Swi + (1-Swi) * Sw
+
+    Perhaps some regressions are made:
+    Swi = a1*poro + a2
+    Pcne = b1*poro^b2
+    N = c1*poro + c2
+
+
+    The Pcn is be normalized capillary pressure, and for gas above an oil zone with
+    thickness H_oil, this will be
+
+    Pcn = H_oil * (gradw - grado)/ado  +  h * (gradw - gradg)/adg
+
+    a = Pcne^(1/N)
+    b = (1/N)
+    x = Pcn
+
+
+    """
+    poro = 0.3
+    a1 = 0.0
+    a2 = 0.3
+    b1 = 4.4e-06
+    b2 = -5.7
+    c1 = -28.45
+    c2 = 10.42
+
+    rw = 0.1002
+    rg = 0.0113
+    ro = 0.0785
+
+    adg = 42.7
+    ado = 26.0
+    h_oil = 3
+
+    ffl = 30.5  # FOL in this example
+
+    grid = xtgeo.create_box_grid((1, 1, 200))
+
+    swi = a1 * poro + a2
+    pcne = b1 * poro**b2
+    n = c1 * poro + c2
+
+    assert swi == pytest.approx(0.3)
+    assert n == pytest.approx(1.885)
+    assert pcne == pytest.approx(0.004205925)
+
+    avalue = pcne ** (1.0 / n)
+    bvalue = -1.0 / n
+    mvalue = h_oil * (rw - ro) / ado
+    xvalue = (rw - rg) / adg
+
+    sw_obj = SwFunction(
+        grid=grid,
+        a=avalue,
+        b=bvalue,
+        m=mvalue,
+        x=xvalue,
+        ffl=ffl,
+        method="cell_center_above_ffl",
+        debug=False,
+    )
+    sw = sw_obj.compute("direct")
+
+    sw10 = float(sw.values[:, :, 20])  # 10 meter above contact
+    assert sw10 == pytest.approx(0.40303321)
+    sw10final = swi + (1 - swi) * sw10
+    assert sw10final == pytest.approx(0.582123)
+
+    sw20 = float(sw.values[:, :, 10])  # 20 meter above contact
+    assert sw20 == pytest.approx(0.28731234)
+
+    sw = sw_obj.compute("integrated")
+    sw10_i = float(sw.values[:, :, 20])
+    assert sw10_i == pytest.approx(sw10, abs=0.001)
diff --git a/tests/rms/test_create_rft_ertobs.py b/tests/rms/test_create_rft_ertobs.py
index fd5f6e5c..7890329b 100644
--- a/tests/rms/test_create_rft_ertobs.py
+++ b/tests/rms/test_create_rft_ertobs.py
@@ -9,9 +9,9 @@
 import numpy as np
 import pandas as pd
 import pytest
+from fmu.tools._common import preserve_cwd
 from fmu.tools.rms import create_rft_ertobs
 from fmu.tools.rms.create_rft_ertobs import check_and_parse_config
-from fmu.tools._common import preserve_cwd
 
 logging.basicConfig(level=logging.INFO)
 
diff --git a/tests/rms/test_import_localmodule.py b/tests/rms/test_import_localmodule.py
index a1c8a32c..84db3005 100644
--- a/tests/rms/test_import_localmodule.py
+++ b/tests/rms/test_import_localmodule.py
@@ -3,7 +3,6 @@
 import os
 import sys
 
-import fmu.tools
 import fmu.tools.rms as toolsrms
 import pytest
 from fmu.tools._common import preserve_cwd
diff --git a/tests/rms/test_rename_rms_scripts.py b/tests/rms/test_rename_rms_scripts.py
index aff8006d..ccc03b3a 100644
--- a/tests/rms/test_rename_rms_scripts.py
+++ b/tests/rms/test_rename_rms_scripts.py
@@ -6,8 +6,8 @@
 from textwrap import dedent
 
 import pytest
-from fmu.tools.rms.rename_rms_scripts import PythonCompMaster, main
 from fmu.tools._common import preserve_cwd
+from fmu.tools.rms.rename_rms_scripts import PythonCompMaster, main
 
 logging.basicConfig(level=logging.INFO)
 
diff --git a/tests/sensitivities/test_designmatrix.py b/tests/sensitivities/test_designmatrix.py
index 786c205a..b20a5e0a 100644
--- a/tests/sensitivities/test_designmatrix.py
+++ b/tests/sensitivities/test_designmatrix.py
@@ -6,9 +6,9 @@
 
 import pandas as pd
 import pytest
+from fmu.tools._common import preserve_cwd
 from fmu.tools.sensitivities import DesignMatrix
 from packaging import version
-from fmu.tools._common import preserve_cwd
 
 TESTDATA = Path(__file__).parent / "data"
 
diff --git a/tests/sensitivities/test_excel2dict.py b/tests/sensitivities/test_excel2dict.py
index 5caa4d2f..a0680059 100644
--- a/tests/sensitivities/test_excel2dict.py
+++ b/tests/sensitivities/test_excel2dict.py
@@ -5,10 +5,10 @@
 import numpy as np
 import pandas as pd
 import pytest
+from fmu.tools._common import preserve_cwd
 from fmu.tools.sensitivities import excel2dict_design, inputdict_to_yaml
 from fmu.tools.sensitivities._excel2dict import _has_value
 from packaging import version
-from fmu.tools._common import preserve_cwd
 
 MOCK_GENERAL_INPUT = pd.DataFrame(
     data=[
diff --git a/tests/sensitivities/test_use_cases.py b/tests/sensitivities/test_use_cases.py
index bbf2efac..8c24c28e 100644
--- a/tests/sensitivities/test_use_cases.py
+++ b/tests/sensitivities/test_use_cases.py
@@ -2,9 +2,9 @@
 
 import pandas as pd
 import pytest
+from fmu.tools._common import preserve_cwd
 from fmu.tools.sensitivities import DesignMatrix, excel2dict_design
 from packaging import version
-from fmu.tools._common import preserve_cwd
 
 
 @preserve_cwd
diff --git a/tests/test_ensembles/test_ensemble_well_props.py b/tests/test_ensembles/test_ensemble_well_props.py
index 66187ec7..e577938e 100644
--- a/tests/test_ensembles/test_ensemble_well_props.py
+++ b/tests/test_ensembles/test_ensemble_well_props.py
@@ -7,8 +7,8 @@
 import pytest
 import xtgeo
 import yaml
-from fmu.tools.ensembles import ensemble_well_props
 from fmu.tools._common import preserve_cwd
+from fmu.tools.ensembles import ensemble_well_props
 
 SOURCE = pathlib.Path(__file__).absolute().parent.parent.parent