def exceed_ns_icpm(cln_min, cln_max, cls_min, cls_max, dep_n, dep_s):
    """Calculates exceedances based on the methodology outlined by Max
        Posch in the ICP Mapping manual (section VII.4):

        https://www.umweltbundesamt.de/sites/default/files/medien/4292/dokumente/ch7-mapman-2016-04-26.pdf

        NB: All units should be in eq/l.

    Args:
        cln_min: Float. Parameter to define "critical load function" (see PDF)
        cln_max: Float. Parameter to define "critical load function" (see PDF)
        cls_min: Float. Parameter to define "critical load function" (see PDF)
        cls_max: Float. Parameter to define "critical load function" (see PDF)
        dep_n:   Float. Total N deposition
        dep_s:   Float. Total (non-marine) S deposition

    Returns:
        Tuple (ex_n, ex_s, reg_id)

        ex_n and ex_s are the exceedances for N and S depositions dep_n and dep_s
        and the CLF defined by (cln_min, cls_max) and (cln_max, cls_min). The
        overall exceedance is (ex_n + ex_s).

        reg_id is an integer region ID, as defined in Figure VII.3 of the PDF.
    """
    # Check inputs
    assert (dep_n >= 0) and (dep_s >= 0), "Deposition cannot be negative."

    # Make sure floats
    cln_min = float(cln_min)
    cln_max = float(cln_max)
    cls_min = float(cls_min)
    cls_max = float(cls_max)
    dep_n = float(dep_n)
    dep_s = float(dep_s)

    # Handle edge cases
    # CLF pars < 0
    if (cln_min < 0) or (cln_max < 0) or (cls_min < 0) or (cls_max < 0):
        # Pars not valid
        # Updated 07.11.2020. Values < 0  do not make sense, so were originally set to -1.
        # This change is equivalent to setting values less than zero back to zero
        # return (-1, -1, -1)
        return (dep_n, dep_s, -1)

    # CL = 0
    if (cls_max == 0) and (cln_max == 0):
        # All dep is above CL
        return (dep_n, dep_s, 9)

    # Otherwise, we're somewhere on Fig. VII.3
    dn = cln_min - cln_max
    ds = cls_max - cls_min

    if (
        (dep_s <= cls_max)
        and (dep_n <= cln_max)
        and ((dep_n - cln_max) * ds <= (dep_s - cls_min) * dn)
    ):
        # Non-exceedance
        return (0, 0, 0)

    elif dep_s <= cls_min:
        # Region 1
        ex_s = 0
        ex_n = dep_n - cln_max

        return (ex_n, ex_s, 1)

    elif dep_n <= cln_min:
        # Region 5
        ex_s = dep_s - cls_max
        ex_n = 0

        return (ex_n, ex_s, 5)

    elif -(dep_n - cln_max) * dn >= (dep_s - cls_min) * ds:
        # Region 2
        ex_n = dep_n - cln_max
        ex_s = dep_s - cls_min

        return (ex_n, ex_s, 2)

    elif -(dep_n - cln_min) * dn <= (dep_s - cls_max) * ds:
        # Region 4
        ex_n = dep_n - cln_min
        ex_s = dep_s - cls_max

        return (ex_n, ex_s, 4)

    else:
        # Region 3
        dd = dn ** 2 + ds ** 2
        s = dep_n * dn + dep_s * ds
        v = cln_max * ds - cls_min * dn
        xf = (dn * s + ds * v) / dd
        yf = (ds * s - dn * v) / dd
        ex_n = dep_n - xf
        ex_s = dep_s - yf

        return (ex_n, ex_s, 3)