diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 7783738..eb3a265 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -4,13 +4,19 @@ on:
   # Trigger the workflow on push to master or develop, except tag creation
   push:
     branches:
-    - 'main'
-    - 'develop'
+    - "main"
+    - "develop"
     tags-ignore:
-    - '**'
+    - "**"
+    paths-ignore:
+    - "docs/**"
+    - "README.md"
 
   # Trigger the workflow on pull request
   pull_request:
+    paths-ignore:
+    - "docs/**"
+    - "README.md"
 
   # Trigger the workflow manually
   workflow_dispatch:
@@ -18,6 +24,9 @@ on:
   # Trigger after public PR approved for CI
   pull_request_target:
     types: [labeled]
+    paths-ignore:
+    - "docs/**"
+    - "README.md"
 
 jobs:
   # Run CI including downstream packages on self-hosted runners
@@ -28,10 +37,9 @@ jobs:
     with:
       earthkit-meteo: ecmwf/earthkit-meteo@${{ github.event.pull_request.head.sha || github.sha }}
       codecov_upload: true
-      python_qa: true
+      python_qa: false
     secrets: inherit
 
-
   # Build downstream packages on HPC
   downstream-ci-hpc:
     name: downstream-ci-hpc
diff --git a/.github/workflows/legacy-ci.yml b/.github/workflows/legacy-ci.yml
index b9619d8..77b015a 100644
--- a/.github/workflows/legacy-ci.yml
+++ b/.github/workflows/legacy-ci.yml
@@ -24,18 +24,6 @@ defaults:
     shell: bash -l {0}
 
 jobs:
-  pre-commit:
-    if: ${{ !github.event.pull_request.head.repo.fork && github.event.action != 'labeled' || github.event.label.name == 'approved-for-ci' }}
-    runs-on: ubuntu-latest
-    steps:
-    - uses: actions/checkout@v3
-      with:
-        ref: ${{ github.event.pull_request.head.sha || github.ref }}
-    - uses: actions/setup-python@v4
-      with:
-        python-version: 3.x
-    - uses: pre-commit/action@v3.0.0
-
   documentation:
     if: ${{ !github.event.pull_request.head.repo.fork && github.event.action != 'labeled' || github.event.label.name == 'approved-for-ci' }}
     runs-on: ubuntu-latest
@@ -48,7 +36,7 @@ jobs:
       with:
         ref: ${{ github.event.pull_request.head.sha || github.ref }}
     - name: Install Conda environment with Micromamba
-      uses: mamba-org/provision-with-micromamba@v12
+      uses: mamba-org/setup-micromamba@v1
       with:
         environment-file: tests/environment-unit-tests.yml
         environment-name: DEVELOP
diff --git a/.github/workflows/python-pill-request.yml b/.github/workflows/python-pill-request.yml
new file mode 100644
index 0000000..2c38063
--- /dev/null
+++ b/.github/workflows/python-pill-request.yml
@@ -0,0 +1,12 @@
+name: Code Quality checks for PRs
+
+on:
+  push:
+  pull_request:
+    types: [opened, synchronize, reopened]
+
+jobs:
+  quality:
+    uses: ecmwf-actions/reusable-workflows/.github/workflows/qa-precommit-run.yml@v2
+    with:
+      skip-hooks: "no-commit-to-branch"
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 5ffa512..260d26e 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -1,32 +1,48 @@
 repos:
 - repo: https://github.com/pre-commit/pre-commit-hooks
-  rev: v4.4.0
+  rev: v5.0.0
   hooks:
-  - id: trailing-whitespace
-  - id: end-of-file-fixer
+  - id: trailing-whitespace     # Trailing whitespace checker
+  - id: end-of-file-fixer     # Ensure files end in a newline
   - id: check-json
-  - id: check-yaml
+  - id: check-yaml     # Check YAML files for syntax errors only
+    args: [--unsafe, --allow-multiple-documents]
   - id: check-toml
-  # - id: check-added-large-files
-  - id: debug-statements
+      # - id: check-added-large-files
+  - id: debug-statements     # Check for debugger imports and py37+ breakpoint()
   - id: mixed-line-ending
+  - id: no-commit-to-branch     # Prevent committing to main / master
+  - id: check-merge-conflict     # Check for files that contain merge conflict
+    exclude: /README\.rst$|^docs/.*\.rst$
 - repo: https://github.com/PyCQA/isort
-  rev: 5.12.0
+  rev: 5.13.2
   hooks:
   - id: isort
+    args:
+    - -l 110
+    - --force-single-line-imports
+    - --profile black
 - repo: https://github.com/psf/black
-  rev: 23.9.1
+  rev: 24.8.0
   hooks:
   - id: black
+    args: [--line-length=110]
 - repo: https://github.com/keewis/blackdoc
   rev: v0.3.8
   hooks:
   - id: blackdoc
     additional_dependencies: [black==23.3.0]
-- repo: https://github.com/PyCQA/flake8
-  rev: 6.1.0
-  hooks:
-  - id: flake8
+    exclude: xr_engine_profile_rst\.py
+- repo: https://github.com/astral-sh/ruff-pre-commit
+  rev: v0.6.9
+  hooks:
+  - id: ruff
+    exclude: '(dev/.*|.*_)\.py$'
+    args:
+    - --line-length=110
+    - --fix
+    - --exit-non-zero-on-fix
+    - --preview
 - repo: https://github.com/executablebooks/mdformat
   rev: 0.7.14
   hooks:
@@ -37,11 +53,11 @@ repos:
   hooks:
   - id: pretty-format-yaml
     args: [--autofix, --preserve-quotes]
-  - id: pretty-format-toml
-    args: [--autofix]
-- repo: https://github.com/PyCQA/pydocstyle.git
-  rev: 6.1.1
-  hooks:
-  - id: pydocstyle
-    additional_dependencies: [toml]
-    exclude: tests|docs
+- repo: https://github.com/sphinx-contrib/sphinx-lint
+  rev: v1.0.0
+  hooks:
+  - id: sphinx-lint
+- repo: https://github.com/tox-dev/pyproject-fmt
+  rev: "2.2.4"
+  hooks:
+  - id: pyproject-fmt
diff --git a/docs/conf.py b/docs/conf.py
index 4854c59..fa0b391 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -36,6 +36,7 @@
     "sphinx.ext.autodoc",
     "sphinx.ext.napoleon",
     "autoapi.extension",
+    "sphinx_issues",
 ]
 
 # autodoc configuration
@@ -62,6 +63,9 @@
 napoleon_numpy_docstring = True
 napoleon_preprocess_types = True
 
+# Path to GitHub repo {group}/{project}  (note that `group` is the GitHub user or organization)
+issues_github_path = "ecmwf/earthkit-meteo"
+
 # Add any paths that contain templates here, relative to this directory.
 templates_path = ["_templates"]
 
diff --git a/docs/release_notes/index.rst b/docs/release_notes/index.rst
index 37ca6f6..3e54796 100644
--- a/docs/release_notes/index.rst
+++ b/docs/release_notes/index.rst
@@ -4,4 +4,5 @@ Release notes
 .. toctree::
     :maxdepth: 1
 
+    version_0.2_updates
     version_0.1_updates
diff --git a/docs/release_notes/version_0.2_updates.rst b/docs/release_notes/version_0.2_updates.rst
new file mode 100644
index 0000000..70dc830
--- /dev/null
+++ b/docs/release_notes/version_0.2_updates.rst
@@ -0,0 +1,15 @@
+Version 0.2 Updates
+/////////////////////////
+
+
+Version 0.2.0
+===============
+
+New features
++++++++++++++++
+
+- Added the ``eps`` and ``out`` keywords arguments to control zero humidity in the following methods (:pr:`23`):
+
+  - :py:meth:`dewpoint_from_specific_humidity <meteo.thermo.array.dewpoint_from_specific_humidity>`
+  - :py:meth:`dewpoint_from_relative_humidity <meteo.thermo.array.dewpoint_from_relative_humidity>`
+  - :py:meth:`temperature_from_saturation_vapour_pressure <meteo.thermo.array.temperature_from_saturation_vapour_pressure>`
diff --git a/docs/requirements.txt b/docs/requirements.txt
index bcfdc7d..ae01e54 100644
--- a/docs/requirements.txt
+++ b/docs/requirements.txt
@@ -7,3 +7,4 @@ sphinx-rtd-theme
 nbsphinx
 setuptools
 sphinx-autoapi
+sphinx-issues
diff --git a/environment.yml b/environment.yml
index ec9248d..6511388 100644
--- a/environment.yml
+++ b/environment.yml
@@ -17,6 +17,7 @@ dependencies:
 - sphinx-autoapi
 - sphinx_rtd_theme
 - sphinxcontrib-apidoc
+- sphinx-issues
 - nbformat
 - nbconvert
 - nbsphinx
diff --git a/pyproject.toml b/pyproject.toml
index 7279832..ad4cc62 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,58 +1,56 @@
 [build-system]
-requires = ["setuptools>=61", "setuptools-scm>=8.0"]
+requires = [ "setuptools>=61", "setuptools-scm>=8" ]
 
 [project]
+name = "earthkit-meteo"
+description = "Meteorological computations"
+readme = "README.md"
+license = { text = "Apache License Version 2.0" }
 authors = [
-  {name = "European Centre for Medium-Range Weather Forecasts (ECMWF)", email = "software.support@ecmwf.int"}
+  { name = "European Centre for Medium-Range Weather Forecasts (ECMWF)", email = "software.support@ecmwf.int" },
 ]
+requires-python = ">=3.8"
+
 classifiers = [
   "Development Status :: 2 - Pre-Alpha",
   "Intended Audience :: Developers",
   "License :: OSI Approved :: Apache Software License",
-  "Programming Language :: Python :: 3",
+  "Operating System :: OS Independent",
+  "Programming Language :: Python :: 3 :: Only",
   "Programming Language :: Python :: 3.8",
   "Programming Language :: Python :: 3.9",
   "Programming Language :: Python :: 3.10",
   "Programming Language :: Python :: 3.11",
+  "Programming Language :: Python :: 3.12",
   "Programming Language :: Python :: Implementation :: CPython",
   "Programming Language :: Python :: Implementation :: PyPy",
-  "Operating System :: OS Independent"
 ]
+dynamic = [ "version" ]
 dependencies = [
-  "numpy"
+  "numpy",
 ]
-description = "Meteorological computations"
-dynamic = ["version"]
-license = {text = "Apache License Version 2.0"}
-name = "earthkit-meteo"
-readme = "README.md"
-requires-python = ">= 3.8"
-
-[project.optional-dependencies]
-test = [
+optional-dependencies.test = [
   "pytest",
-  "pytest-cov"
+  "pytest-cov",
 ]
+urls.Documentation = "https://earthkit-meteo.readthedocs.io/"
+urls.Homepage = "https://github.com/ecmwf/earthkit-meteo/"
+urls.Issues = "https://github.com/ecmwf/earthkit-meteo.issues"
+urls.Repository = "https://github.com/ecmwf/earthkit-meteo/"
 
-[project.urls]
-Documentation = "https://earthkit-meteo.readthedocs.io/"
-Homepage = "https://github.com/ecmwf/earthkit-meteo/"
-Issues = "https://github.com/ecmwf/earthkit-meteo.issues"
-Repository = "https://github.com/ecmwf/earthkit-meteo/"
+[tool.setuptools.packages.find]
+include = [ "earthkit.meteo" ]
+where = [ "src/" ]
 
-[tool.coverage.run]
-branch = "true"
+[tool.setuptools_scm]
+version_file = "src/earthkit/meteo/_version.py"
 
 [tool.isort]
 profile = "black"
 
+[tool.coverage.run]
+branch = "true"
+
 [tool.pydocstyle]
-add_ignore = ["D1", "D200", "D205", "D400", "D401"]
+add_ignore = [ "D1", "D200", "D205", "D400", "D401" ]
 convention = "numpy"
-
-[tool.setuptools.packages.find]
-include = ["earthkit.meteo"]
-where = ["src/"]
-
-[tool.setuptools_scm]
-version_file = "src/earthkit/meteo/_version.py"
diff --git a/src/earthkit/meteo/extreme/array/cpf.py b/src/earthkit/meteo/extreme/array/cpf.py
index 16596de..1a737a9 100644
--- a/src/earthkit/meteo/extreme/array/cpf.py
+++ b/src/earthkit/meteo/extreme/array/cpf.py
@@ -64,10 +64,9 @@ def cpf(clim, ens, sort_clim=True, sort_ens=True):
                     idx = (qv_f < qv_c) & (qv_c_2 < qv_c)
 
                     # intersection between two lines
-                    tau_i = (
-                        tau_c * (qv_c_2[idx] - qv_f[idx])
-                        + tau_c_2 * (qv_f[idx] - qv_c[idx])
-                    ) / (qv_c_2[idx] - qv_c[idx])
+                    tau_i = (tau_c * (qv_c_2[idx] - qv_f[idx]) + tau_c_2 * (qv_f[idx] - qv_c[idx])) / (
+                        qv_c_2[idx] - qv_c[idx]
+                    )
 
                     # populate matrix, no values below 0
                     cpf[idx] = np.maximum(tau_i, 0)
@@ -85,10 +84,9 @@ def cpf(clim, ens, sort_clim=True, sort_ens=True):
 
                     idx = (qv_f > qv_c) & (qv_c_2 > qv_c) & (~mask)
 
-                    tau_i = (
-                        tau_c * (qv_c_2[idx] - qv_f[idx])
-                        + tau_c_2 * (qv_f[idx] - qv_c[idx])
-                    ) / (qv_c_2[idx] - qv_c[idx])
+                    tau_i = (tau_c * (qv_c_2[idx] - qv_f[idx]) + tau_c_2 * (qv_f[idx] - qv_c[idx])) / (
+                        qv_c_2[idx] - qv_c[idx]
+                    )
 
                     # populate matrix, no values above 1
                     cpf[idx] = np.minimum(tau_i, 1)
diff --git a/src/earthkit/meteo/extreme/array/efi.py b/src/earthkit/meteo/extreme/array/efi.py
index 455222b..6c08911 100644
--- a/src/earthkit/meteo/extreme/array/efi.py
+++ b/src/earthkit/meteo/extreme/array/efi.py
@@ -31,9 +31,7 @@ def efi(clim, ens, eps=-0.1):
         EFI values
     """
     # locate missing values
-    missing_mask = np.logical_or(
-        np.sum(np.isnan(clim), axis=0), np.sum(np.isnan(ens), axis=0)
-    )
+    missing_mask = np.logical_or(np.sum(np.isnan(clim), axis=0), np.sum(np.isnan(ens), axis=0))
 
     # Compute fraction of the forecast below climatology
     nclim, npoints = clim.shape
@@ -65,9 +63,7 @@ def efi(clim, ens, eps=-0.1):
             mask = clim[icl + 1, :] > eps
             dEFI = np.where(
                 mask,
-                (2.0 * frac[icl, :] - 1.0) * acosdiff[icl]
-                + acoef[icl] * dFdp[icl, :]
-                - proddiff[icl],
+                (2.0 * frac[icl, :] - 1.0) * acosdiff[icl] + acoef[icl] * dFdp[icl, :] - proddiff[icl],
                 0.0,
             )
             defimax = np.where(mask, -acosdiff[icl] - proddiff[icl], 0.0)
@@ -77,11 +73,7 @@ def efi(clim, ens, eps=-0.1):
         efi /= efimax
     else:
         for icl in range(nclim - 1):
-            dEFI = (
-                (2.0 * frac[icl, :] - 1.0) * acosdiff[icl]
-                + acoef[icl] * dFdp[icl, :]
-                - proddiff[icl]
-            )
+            dEFI = (2.0 * frac[icl, :] - 1.0) * acosdiff[icl] + acoef[icl] * dFdp[icl, :] - proddiff[icl]
             efi += dEFI
         efi *= 2.0 / np.pi
     ##################################
diff --git a/src/earthkit/meteo/extreme/array/sot.py b/src/earthkit/meteo/extreme/array/sot.py
index f556afa..8a7d933 100644
--- a/src/earthkit/meteo/extreme/array/sot.py
+++ b/src/earthkit/meteo/extreme/array/sot.py
@@ -39,9 +39,7 @@ def sot_func(qc_tail, qc, qf, eps=-1e-4, lower_bound=-10, upper_bound=10):
     err = np.seterr(divide="ignore", invalid="ignore")
 
     min_den = np.fmax(eps, 0)
-    sot = np.where(
-        np.fabs(qc_tail - qc) > min_den, (qf - qc_tail) / (qc_tail - qc), np.nan
-    )
+    sot = np.where(np.fabs(qc_tail - qc) > min_den, (qf - qc_tail) / (qc_tail - qc), np.nan)
 
     # revert to original error state
     np.seterr(**err)
@@ -78,20 +76,12 @@ def sot(clim, ens, perc, eps=-1e4):
     numpy array (npoints)
         SOT values
     """
-    if not (isinstance(perc, int) or isinstance(perc, np.int64)) or (
-        perc < 2 or perc > 98
-    ):
-        raise Exception(
-            "Percentile value should be and Integer between 2 and 98, is {}".format(
-                perc
-            )
-        )
+    if not (isinstance(perc, int) or isinstance(perc, np.int64)) or (perc < 2 or perc > 98):
+        raise Exception("Percentile value should be and Integer between 2 and 98, is {}".format(perc))
 
     if clim.shape[0] != 101:
         raise Exception(
-            "Climatology array should contain 101 percentiles, it has {} values".format(
-                clim.shape
-            )
+            "Climatology array should contain 101 percentiles, it has {} values".format(clim.shape)
         )
 
     qc = clim[perc]
@@ -136,20 +126,12 @@ def sot_unsorted(clim, ens, perc, eps=-1e4):
     numpy array (npoints)
         SOT values
     """
-    if not (isinstance(perc, int) or isinstance(perc, np.int64)) or (
-        perc < 2 or perc > 98
-    ):
-        raise Exception(
-            "Percentile value should be and Integer between 2 and 98, is {}".format(
-                perc
-            )
-        )
+    if not (isinstance(perc, int) or isinstance(perc, np.int64)) or (perc < 2 or perc > 98):
+        raise Exception("Percentile value should be and Integer between 2 and 98, is {}".format(perc))
 
     if clim.shape[0] != 101:
         raise Exception(
-            "Climatology array should contain 101 percentiles, it has {} values".format(
-                clim.shape
-            )
+            "Climatology array should contain 101 percentiles, it has {} values".format(clim.shape)
         )
 
     if eps > 0:
diff --git a/src/earthkit/meteo/score/array/crps.py b/src/earthkit/meteo/score/array/crps.py
index 9595094..dcce923 100644
--- a/src/earthkit/meteo/score/array/crps.py
+++ b/src/earthkit/meteo/score/array/crps.py
@@ -48,9 +48,7 @@ def crps(x, y):
     alpha[-1] = np.fmax(-diffxy[-1], 0)  # y-x(n)
     beta[-1] = 0
     # else
-    alpha[1:-1] = np.fmin(
-        diffxx, np.fmax(-diffxy[:-1], 0)
-    )  # x(i+1)-x(i) or y-x(i) or 0
+    alpha[1:-1] = np.fmin(diffxx, np.fmax(-diffxy[:-1], 0))  # x(i+1)-x(i) or y-x(i) or 0
     beta[1:-1] = np.fmin(diffxx, np.fmax(diffxy[1:], 0))  # 0 or x(i+1)-y or x(i+1)-x(i)
 
     # compute crps
diff --git a/src/earthkit/meteo/stats/array/quantiles.py b/src/earthkit/meteo/stats/array/quantiles.py
index 0f9109a..b590a1a 100644
--- a/src/earthkit/meteo/stats/array/quantiles.py
+++ b/src/earthkit/meteo/stats/array/quantiles.py
@@ -7,7 +7,9 @@
 # nor does it submit to any jurisdiction.
 #
 
-from typing import Iterable, List, Union
+from typing import Iterable
+from typing import List
+from typing import Union
 
 import numpy as np
 
@@ -41,9 +43,7 @@ def iter_quantiles(
         Quantiles, in increasing order if `which` is an `int`, otherwise in the order specified
     """
     if method not in ("sort", "numpy_bulk", "numpy"):
-        raise ValueError(
-            f"Invalid method {method!r}, expected 'sort', 'numpy_bulk', or 'numpy'"
-        )
+        raise ValueError(f"Invalid method {method!r}, expected 'sort', 'numpy_bulk', or 'numpy'")
 
     if isinstance(which, int):
         n = which
diff --git a/src/earthkit/meteo/thermo/array/thermo.py b/src/earthkit/meteo/thermo/array/thermo.py
index f0c2b7a..4155d1a 100644
--- a/src/earthkit/meteo/thermo/array/thermo.py
+++ b/src/earthkit/meteo/thermo/array/thermo.py
@@ -13,6 +13,10 @@
 from earthkit.meteo import constants
 
 
+def _valid_number(x):
+    return x is not None and not np.isnan(x)
+
+
 def celsius_to_kelvin(t):
     """Converts temperature values from Celsius to Kelvin.
 
@@ -182,9 +186,7 @@ def specific_humidity_from_vapour_pressure(e, p, eps=1e-4):
 
     """
     if eps <= 0:
-        raise ValueError(
-            f"specific_humidity_from_vapour_pressure(): eps={eps} must be > 0"
-        )
+        raise ValueError(f"specific_humidity_from_vapour_pressure(): eps={eps} must be > 0")
 
     v = np.asarray(p + (constants.epsilon - 1) * e)
     v[np.asarray(p - e) < eps] = np.nan
@@ -261,9 +263,21 @@ def compute_slope(self, t, phase):
         elif phase == "ice":
             return self._es_ice_slope(t)
 
-    def t_from_es(self, es):
-        v = np.log(es / self.c1)
-        return (v * self.c4w - self.c3w * self.t0) / (v - self.c3w)
+    def t_from_es(self, es, eps=1e-8, out=None):
+        def _comp(x):
+            x = x / self.c1
+            x = np.log(x)
+            return (x * self.c4w - self.c3w * self.t0) / (x - self.c3w)
+
+        if out is not None:
+            v = np.asarray(es)
+            z_mask = v < eps
+            v_mask = ~z_mask
+            v[v_mask] = _comp(v[v_mask])
+            v[z_mask] = out
+            return v
+        else:
+            return _comp(es)
 
     def _es_water(self, t):
         return self.c1 * np.exp(self.c3w * (t - self.t0) / (t - self.c4w))
@@ -294,24 +308,14 @@ def _es_mixed(self, t):
         # mixed range
         m_mask = ~(i_mask | w_mask)
         alpha = np.square(t[m_mask] - self.ti) / np.square(self.t0 - self.ti)
-        svp[m_mask] = alpha * self._es_water(t[m_mask]) + (1.0 - alpha) * self._es_ice(
-            t[m_mask]
-        )
+        svp[m_mask] = alpha * self._es_water(t[m_mask]) + (1.0 - alpha) * self._es_ice(t[m_mask])
         return svp
 
     def _es_water_slope(self, t):
-        return (
-            self._es_water(t)
-            * (self.c3w * (self.t0 - self.c4w))
-            / np.square(t - self.c4w)
-        )
+        return self._es_water(t) * (self.c3w * (self.t0 - self.c4w)) / np.square(t - self.c4w)
 
     def _es_ice_slope(self, t):
-        return (
-            self._es_ice(t)
-            * (self.c3i * (self.t0 - self.c4i))
-            / np.square(t - self.c4i)
-        )
+        return self._es_ice(t) * (self.c3i * (self.t0 - self.c4i)) / np.square(t - self.c4i)
 
     def _es_mixed_slope(self, t):
         t = np.asarray(t)
@@ -467,9 +471,7 @@ def saturation_vapour_pressure_slope(t, phase="mixed"):
     return _EsComp().compute_slope(t, phase)
 
 
-def saturation_mixing_ratio_slope(
-    t, p, es=None, es_slope=None, phase="mixed", eps=1e-4
-):
+def saturation_mixing_ratio_slope(t, p, es=None, es_slope=None, phase="mixed", eps=1e-4):
     r"""Computes the slope of saturation mixing ratio with respect to temperature.
 
     Parameters
@@ -519,9 +521,7 @@ def saturation_mixing_ratio_slope(
     return constants.epsilon * es_slope * p / np.square(v)
 
 
-def saturation_specific_humidity_slope(
-    t, p, es=None, es_slope=None, phase="mixed", eps=1e-4
-):
+def saturation_specific_humidity_slope(t, p, es=None, es_slope=None, phase="mixed", eps=1e-4):
     r"""Computes the slope of saturation specific humidity with respect to temperature.
 
     Parameters
@@ -572,13 +572,22 @@ def saturation_specific_humidity_slope(
     return constants.epsilon * es_slope * p / v
 
 
-def temperature_from_saturation_vapour_pressure(es):
-    r"""Computes the temperature from saturation vapour pressure.
+def temperature_from_saturation_vapour_pressure(es, eps=1e-8, out=None):
+    r"""Compute the temperature from saturation vapour pressure.
 
     Parameters
     ----------
     es: ndarray
         :func:`saturation_vapour_pressure` (Pa)
+    eps: number
+        If ``out`` is not None, return ``out`` when ``es`` < ``eps``. If out
+        is None, ``eps`` is ignored and return np.nan for ``es`` values very
+        close to zero. *New in version 0.2.0*
+    out: number or None
+        If not None, return ``out`` when ``es`` < ``eps``. If None, ``eps`` is
+        ignored and return np.nan for ``es`` values very close to zero.
+        *New in version 0.2.0*
+
 
     Returns
     -------
@@ -591,7 +600,7 @@ def temperature_from_saturation_vapour_pressure(es):
     phase ``es`` was computed to.
 
     """
-    return _EsComp().t_from_es(es)
+    return _EsComp().t_from_es(es, eps=eps, out=out)
 
 
 def relative_humidity_from_dewpoint(t, td):
@@ -766,8 +775,8 @@ def specific_humidity_from_relative_humidity(t, r, p):
     return specific_humidity_from_vapour_pressure(e, p)
 
 
-def dewpoint_from_relative_humidity(t, r):
-    r"""Computes the dewpoint temperature from relative humidity.
+def dewpoint_from_relative_humidity(t, r, eps=1e-8, out=None):
+    r"""Compute the dewpoint temperature from relative humidity.
 
     Parameters
     ----------
@@ -775,11 +784,19 @@ def dewpoint_from_relative_humidity(t, r):
         Temperature (K)
     r: ndarray
         Relative humidity (%)
+    eps: number
+        If ``out`` is not None, return ``out`` when ``r`` < ``eps``.
+        If out is None, ``eps`` is ignored and return np.nan for ``r``
+        values very close to zero. *New in version 0.2.0*
+    out: number or None
+        If not None, return ``out`` when ``r`` < ``eps``. If None, ``eps`` is
+        ignored and return np.nan for ``r`` values very close to zero.
+        *New in version 0.2.0*
 
     Returns
     -------
     ndarray
-        Dewpoint (K)
+        Dewpoint temperature (K)
 
 
     The computation starts with determining the the saturation vapour pressure over
@@ -798,12 +815,25 @@ def dewpoint_from_relative_humidity(t, r):
     equations used in :func:`saturation_vapour_pressure`.
 
     """
+    # by masking upfront we avoid RuntimeWarning when calling log() in
+    # the computation of td when r is very small
+    if out is not None:
+        r = np.asarray(r)
+        mask = r < eps
+        r[mask] = np.nan
+
     es = saturation_vapour_pressure(t, phase="water") * r / 100.0
-    return temperature_from_saturation_vapour_pressure(es)
+    td = temperature_from_saturation_vapour_pressure(es)
+
+    if out is not None and not np.isnan(out):
+        td = np.asarray(td)
+        td[mask] = out
 
+    return td
 
-def dewpoint_from_specific_humidity(q, p):
-    r"""Computes the dewpoint temperature from specific humidity.
+
+def dewpoint_from_specific_humidity(q, p, eps=1e-8, out=None):
+    r"""Compute the dewpoint temperature from specific humidity.
 
     Parameters
     ----------
@@ -811,11 +841,19 @@ def dewpoint_from_specific_humidity(q, p):
         Specific humidity (kg/kg)
     p: ndarray
         Pressure (Pa)
+    eps: number
+        If ``out`` is not None, return ``out`` when ``q`` < ``eps``.
+        If out is None, ``eps`` is ignored and return np.nan for ``q``
+        values very close to zero. *New in version 0.2.0*
+    out: number or None
+        If not None, return ``out`` when ``q`` < ``eps``. If None, ``eps`` is
+        ignored and return np.nan for ``q`` values very close to zero.
+        *New in version 0.2.0*
 
     Returns
     -------
     ndarray
-        Dewpoint (K)
+        Dewpoint temperature (K)
 
 
     The computation starts with determining the the saturation vapour pressure over
@@ -835,9 +873,21 @@ def dewpoint_from_specific_humidity(q, p):
     used in :func:`saturation_vapour_pressure`.
 
     """
-    return temperature_from_saturation_vapour_pressure(
-        vapour_pressure_from_specific_humidity(q, p)
-    )
+    # by masking upfront we avoid RuntimeWarning when calling log() in
+    # the computation of td when q is very small
+    if out is not None:
+        q = np.asarray(q)
+        mask = q < eps
+        q[mask] = np.nan
+
+    es = vapour_pressure_from_specific_humidity(q, p)
+    td = temperature_from_saturation_vapour_pressure(es)
+
+    if out is not None and not np.isnan(out):
+        td = np.asarray(td)
+        td[mask] = out
+
+    return td
 
 
 def virtual_temperature(t, q):
@@ -846,7 +896,7 @@ def virtual_temperature(t, q):
     Parameters
     ----------
     t: number or ndarray
-        Temperature (K)
+        Temperature (K)s
     q: number or ndarray
         Specific humidity (kg/kg)
 
@@ -1057,9 +1107,7 @@ def lcl_temperature(t, td, method="davies"):
     """
     # Davies-Jones formula
     if method == "davies":
-        t_lcl = td - (
-            0.212 + 1.571e-3 * (td - constants.T0) - 4.36e-4 * (t - constants.T0)
-        ) * (t - td)
+        t_lcl = td - (0.212 + 1.571e-3 * (td - constants.T0) - 4.36e-4 * (t - constants.T0)) * (t - td)
         return t_lcl
     # Bolton formula
     elif method == "bolton":
@@ -1144,10 +1192,7 @@ def compute_wbpt(self, ept):
         x = ept / t0
         a = [7.101574, -20.68208, 16.11182, 2.574631, -5.205688]
         b = [1.0, -3.552497, 3.781782, -0.6899655, -0.5929340]
-        return ept - np.exp(
-            np.polynomial.polynomial.polyval(x, a)
-            / np.polynomial.polynomial.polyval(x, b)
-        )
+        return ept - np.exp(np.polynomial.polynomial.polyval(x, a) / np.polynomial.polynomial.polyval(x, b))
 
     def compute_t_on_ma_stipanuk(self, ept, p):
         if isinstance(p, np.ndarray):
@@ -1161,10 +1206,7 @@ def compute_t_on_ma_stipanuk(self, ept, p):
         for _ in range(max_iter):
             ths = _ThermoState(t=t, p=p)
             dt /= 2.0
-            t += (
-                np.sign(ept * np.exp(self._G_sat(ths, scale=-1.0)) - self._th_sat(ths))
-                * dt
-            )
+            t += np.sign(ept * np.exp(self._G_sat(ths, scale=-1.0)) - self._th_sat(ths)) * dt
         # ths.t = t
         # return ept - self._th_sat(ths) * np.exp(self._G_sat(ths))
         return t
@@ -1211,11 +1253,7 @@ def _D(p):
         tw[mask] = (_k1(pp[mask]) - 1.21) - (_k2(pp[mask]) - 1.21) * c_te[mask]
 
         mask = c_te < 0.4
-        tw[mask] = (
-            (_k1(pp[mask]) - 2.66)
-            - (_k2(pp[mask]) - 1.21) * c_te[mask]
-            + 0.58 / c_te[mask]
-        )
+        tw[mask] = (_k1(pp[mask]) - 2.66) - (_k2(pp[mask]) - 1.21) * c_te[mask] + 0.58 / c_te[mask]
         # tw has to be converted to Kelvin
         tw = celsius_to_kelvin(tw)
 
@@ -1306,9 +1344,7 @@ def _ept(self, ths):
     def _th_sat(self, ths):
         if ths.ws is None:
             ths.ws = saturation_mixing_ratio(ths.t, ths.p)
-        return ths.t * np.power(
-            constants.p0 / ths.p, constants.kappa * (1 - self.K3 * ths.ws)
-        )
+        return ths.t * np.power(constants.p0 / ths.p, constants.kappa * (1 - self.K3 * ths.ws))
 
     def _G_sat(self, ths, scale=1.0):
         if ths.ws is None:
@@ -1334,9 +1370,7 @@ def _f(self, ths):
     def _d_lnf(self, ths):
         return -self.c_lambda * (
             1 / ths.t
-            + self.K3
-            * np.log(ths.p / constants.p0)
-            * saturation_vapour_pressure_slope(ths.t)
+            + self.K3 * np.log(ths.p / constants.p0) * saturation_vapour_pressure_slope(ths.t)
             + self._d_G_sat(ths)
         )
 
@@ -1364,9 +1398,7 @@ def _ept(self, ths):
             w = mixing_ratio_from_specific_humidity(ths.q)
 
         e = vapour_pressure_from_mixing_ratio(w, ths.p)
-        th = potential_temperature(ths.t, ths.p - e) * np.power(
-            ths.t / t_lcl, self.K4 * w
-        )
+        th = potential_temperature(ths.t, ths.p - e) * np.power(ths.t / t_lcl, self.K4 * w)
         return th * np.exp((self.K0 / t_lcl - self.K1) * w * (1.0 + self.K2 * w))
 
     def _th_sat(self, ths):
@@ -1382,36 +1414,24 @@ def _G_sat(self, ths, scale=1.0):
         # print(f" es={ths.es}")
         ws = mixing_ratio_from_vapour_pressure(ths.es, ths.p)
         # print(f" ws={ws}")
-        return (
-            ((scale * self.K0) / ths.t - (scale * self.K1)) * ws * (1.0 + self.K2 * ws)
-        )
+        return ((scale * self.K0) / ths.t - (scale * self.K1)) * ws * (1.0 + self.K2 * ws)
 
     def _d_G_sat(self, ths):
         # print(f" d_ws={saturation_mixing_ratio_slope(ths.t, ths.p)}")
-        return -self.K0 * (ths.ws + self.K2 * np.square(ths.ws)) / (
-            np.square(ths.t)
-        ) + (self.K0 / ths.t - self.K1) * (
-            1 + (2 * self.K2) * ths.ws
-        ) * saturation_mixing_ratio_slope(
-            ths.t, ths.p
-        )
+        return -self.K0 * (ths.ws + self.K2 * np.square(ths.ws)) / (np.square(ths.t)) + (
+            self.K0 / ths.t - self.K1
+        ) * (1 + (2 * self.K2) * ths.ws) * saturation_mixing_ratio_slope(ths.t, ths.p)
 
     def _f(self, ths):
         # print(f" c_tw={ths.c_tw}")
         # print(f" es_frac={ths.es / ths.p}")
         # print(f" exp={self._G_sat(ths, scale=-self.c_lambda)}")
-        return (
-            ths.c_tw
-            * (1 - ths.es / ths.p)
-            * np.exp(self._G_sat(ths, scale=-self.c_lambda))
-        )
+        return ths.c_tw * (1 - ths.es / ths.p) * np.exp(self._G_sat(ths, scale=-self.c_lambda))
 
     def _d_lnf(self, ths):
         return -self.c_lambda * (
             1 / ths.t
-            + constants.kappa
-            * saturation_vapour_pressure_slope(ths.t)
-            / (ths.p - ths.es)
+            + constants.kappa * saturation_vapour_pressure_slope(ths.t) / (ths.p - ths.es)
             + self._d_G_sat(ths)
         )
 
@@ -1606,9 +1626,7 @@ def temperature_on_moist_adiabat(ept, p, ept_method="ifs", t_method="bisect"):
     elif t_method == "newton":
         return cm.compute_t_on_ma_davies(ept, p)
     else:
-        raise ValueError(
-            f"temperature_on_moist_adiabat: invalid t_method={t_method} specified!"
-        )
+        raise ValueError(f"temperature_on_moist_adiabat: invalid t_method={t_method} specified!")
 
 
 def wet_bulb_temperature_from_dewpoint(t, td, p, ept_method="ifs", t_method="bisect"):
@@ -1648,14 +1666,10 @@ def wet_bulb_temperature_from_dewpoint(t, td, p, ept_method="ifs", t_method="bis
 
     """
     ept = ept_from_dewpoint(t, td, p, method=ept_method)
-    return temperature_on_moist_adiabat(
-        ept, p, ept_method=ept_method, t_method=t_method
-    )
+    return temperature_on_moist_adiabat(ept, p, ept_method=ept_method, t_method=t_method)
 
 
-def wet_bulb_temperature_from_specific_humidity(
-    t, q, p, ept_method="ifs", t_method="bisect"
-):
+def wet_bulb_temperature_from_specific_humidity(t, q, p, ept_method="ifs", t_method="bisect"):
     r"""Computes the pseudo adiabatic wet bulb temperature from specific humidity.
 
     Parameters
@@ -1693,14 +1707,10 @@ def wet_bulb_temperature_from_specific_humidity(
 
     """
     ept = ept_from_specific_humidity(t, q, p, method=ept_method)
-    return temperature_on_moist_adiabat(
-        ept, p, ept_method=ept_method, t_method=t_method
-    )
+    return temperature_on_moist_adiabat(ept, p, ept_method=ept_method, t_method=t_method)
 
 
-def wet_bulb_potential_temperature_from_dewpoint(
-    t, td, p, ept_method="ifs", t_method="direct"
-):
+def wet_bulb_potential_temperature_from_dewpoint(t, td, p, ept_method="ifs", t_method="direct"):
     r"""Computes the pseudo adiabatic wet bulb potential temperature from dewpoint.
 
     Parameters
@@ -1741,14 +1751,10 @@ def wet_bulb_potential_temperature_from_dewpoint(
     if t_method == "direct":
         return _EptComp.make(ept_method).compute_wbpt(ept)
     else:
-        return temperature_on_moist_adiabat(
-            ept, constants.p0, ept_method=ept_method, t_method=t_method
-        )
+        return temperature_on_moist_adiabat(ept, constants.p0, ept_method=ept_method, t_method=t_method)
 
 
-def wet_bulb_potential_temperature_from_specific_humidity(
-    t, q, p, ept_method="ifs", t_method="direct"
-):
+def wet_bulb_potential_temperature_from_specific_humidity(t, q, p, ept_method="ifs", t_method="direct"):
     r"""Computes the pseudo adiabatic wet bulb potential temperature from specific humidity.
 
     Parameters
@@ -1786,6 +1792,4 @@ def wet_bulb_potential_temperature_from_specific_humidity(
     if t_method == "direct":
         return _EptComp.make(ept_method).compute_wbpt(ept)
     else:
-        return temperature_on_moist_adiabat(
-            ept, constants.p0, ept_method=ept_method, t_method=t_method
-        )
+        return temperature_on_moist_adiabat(ept, constants.p0, ept_method=ept_method, t_method=t_method)
diff --git a/src/earthkit/meteo/vertical/__init__.py b/src/earthkit/meteo/vertical/__init__.py
new file mode 100644
index 0000000..5552f22
--- /dev/null
+++ b/src/earthkit/meteo/vertical/__init__.py
@@ -0,0 +1,19 @@
+# (C) Copyright 2021 ECMWF.
+#
+# This software is licensed under the terms of the Apache Licence Version 2.0
+# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+# In applying this licence, ECMWF does not waive the privileges and immunities
+# granted to it by virtue of its status as an intergovernmental organisation
+# nor does it submit to any jurisdiction.
+#
+
+"""
+Vertical computation functions.
+
+The API is split into two levels. The low level functions are in the ``array`` submodule and they
+can be used to operate on numpy arrays. The high level functions are still to be developed and
+planned to work with objects like *earthkit.data FieldLists* or *xarray DataSets*.
+"""
+
+
+from .vertical import *  # noqa
diff --git a/src/earthkit/meteo/vertical/array/__init__.py b/src/earthkit/meteo/vertical/array/__init__.py
new file mode 100644
index 0000000..9884ca5
--- /dev/null
+++ b/src/earthkit/meteo/vertical/array/__init__.py
@@ -0,0 +1,14 @@
+# (C) Copyright 2021 ECMWF.
+#
+# This software is licensed under the terms of the Apache Licence Version 2.0
+# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+# In applying this licence, ECMWF does not waive the privileges and immunities
+# granted to it by virtue of its status as an intergovernmental organisation
+# nor does it submit to any jurisdiction.
+#
+
+"""
+Vertical computation functions operating on numpy arrays.
+"""
+
+from .vertical import *  # noqa
diff --git a/src/earthkit/meteo/vertical/array/vertical.py b/src/earthkit/meteo/vertical/array/vertical.py
new file mode 100644
index 0000000..975426c
--- /dev/null
+++ b/src/earthkit/meteo/vertical/array/vertical.py
@@ -0,0 +1,8 @@
+# (C) Copyright 2021 ECMWF.
+#
+# This software is licensed under the terms of the Apache Licence Version 2.0
+# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+# In applying this licence, ECMWF does not waive the privileges and immunities
+# granted to it by virtue of its status as an intergovernmental organisation
+# nor does it submit to any jurisdiction.
+#
diff --git a/src/earthkit/meteo/vertical/vertical.py b/src/earthkit/meteo/vertical/vertical.py
new file mode 100644
index 0000000..5a1d5be
--- /dev/null
+++ b/src/earthkit/meteo/vertical/vertical.py
@@ -0,0 +1,10 @@
+# (C) Copyright 2021 ECMWF.
+#
+# This software is licensed under the terms of the Apache Licence Version 2.0
+# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+# In applying this licence, ECMWF does not waive the privileges and immunities
+# granted to it by virtue of its status as an intergovernmental organisation
+# nor does it submit to any jurisdiction.
+#
+
+# from . import array
diff --git a/src/earthkit/meteo/wind/array/wind.py b/src/earthkit/meteo/wind/array/wind.py
index f611ca9..a20a529 100644
--- a/src/earthkit/meteo/wind/array/wind.py
+++ b/src/earthkit/meteo/wind/array/wind.py
@@ -279,13 +279,9 @@ def windrose(speed, direction, sectors=16, speed_bins=[], percent=True):
     speed = np.atleast_1d(speed)
     direction = np.atleast_1d(direction)
     dir_step = 360.0 / sectors
-    dir_bins = np.linspace(
-        int(-dir_step / 2), int(360 + dir_step / 2), int(360 / dir_step) + 2
-    )
+    dir_bins = np.linspace(int(-dir_step / 2), int(360 + dir_step / 2), int(360 / dir_step) + 2)
 
-    res = np.histogram2d(speed, direction, bins=[speed_bins, dir_bins], density=False)[
-        0
-    ]
+    res = np.histogram2d(speed, direction, bins=[speed_bins, dir_bins], density=False)[0]
 
     # unify the north bins
     res[:, 0] = res[:, 0] + res[:, -1]
diff --git a/tests/environment-unit-tests.yml b/tests/environment-unit-tests.yml
index 17eefa1..78e3148 100644
--- a/tests/environment-unit-tests.yml
+++ b/tests/environment-unit-tests.yml
@@ -16,6 +16,7 @@ dependencies:
 - sphinx-autoapi
 - sphinx_rtd_theme
 - sphinxcontrib-apidoc
+- sphinx-issues
 - nbformat
 - nbconvert
 - nbsphinx
diff --git a/tests/score/test_score.py b/tests/score/test_score.py
index aa80c04..ffa8cc1 100644
--- a/tests/score/test_score.py
+++ b/tests/score/test_score.py
@@ -43,13 +43,9 @@ def crps_quaver2(x, y):
         lcond = esarr[0, :] > anarr
         aa[0, lcond] = 1.0
         bb[0, :] = np.where(lcond, esarr[0, :] - anarr, 0.0)
-        aa[1:-1, :] = np.where(
-            esarr[1:, :] <= anarr, esarr[1:, :] - esarr[:-1, :], anarr - esarr[:-1, :]
-        )
+        aa[1:-1, :] = np.where(esarr[1:, :] <= anarr, esarr[1:, :] - esarr[:-1, :], anarr - esarr[:-1, :])
         aa[1:-1, :][esarr[: n_ens - 1, :] > anarr] = 0.0  # this would be hard in xarray
-        bb[1:-1, :] = np.where(
-            esarr[:-1, :] > anarr, esarr[1:, :] - esarr[:-1, :], esarr[1:, :] - anarr
-        )
+        bb[1:-1, :] = np.where(esarr[:-1, :] > anarr, esarr[1:, :] - esarr[:-1, :], esarr[1:, :] - anarr)
         bb[1:-1, :][esarr[1:, :] <= anarr] = 0.0
         lcond = anarr > esarr[-1, :]
         aa[-1, :] = np.where(lcond, anarr - esarr[-1, :], 0.0)
diff --git a/tests/solar/test_solar.py b/tests/solar/test_solar.py
index bb362b8..64d00f4 100644
--- a/tests/solar/test_solar.py
+++ b/tests/solar/test_solar.py
@@ -58,9 +58,7 @@ def test_cos_solar_zenith_angle_integrated():
     latitudes = np.array([40.0])
     longitudes = np.array([18.0])
 
-    v = solar.cos_solar_zenith_angle_integrated(
-        begin_date, end_date, latitudes, longitudes
-    )
+    v = solar.cos_solar_zenith_angle_integrated(begin_date, end_date, latitudes, longitudes)
     assert np.isclose(v[0], 0.3110738757)
 
 
diff --git a/tests/stats/test_stats.py b/tests/stats/test_stats.py
index ba7aea4..03aa327 100644
--- a/tests/stats/test_stats.py
+++ b/tests/stats/test_stats.py
@@ -26,9 +26,7 @@ def test_nanaverage():
     # replace nan values
     target_result[:, 1] = np.nansum(data, axis=-1)[:, 1]
     print(target_result)
-    np.testing.assert_equal(
-        stats.nanaverage(data, axis=-1, weights=weights), target_result
-    )
+    np.testing.assert_equal(stats.nanaverage(data, axis=-1, weights=weights), target_result)
 
 
 @pytest.mark.parametrize("method", ["sort", "numpy_bulk", "numpy"])
@@ -71,10 +69,6 @@ def test_quantiles_nans():
     arr = np.random.rand(100, 100, 100)
     arr.ravel()[np.random.choice(arr.size, 100000, replace=False)] = np.nan
     qs = [0.0, 0.25, 0.5, 0.75, 1.0]
-    sort = [
-        quantile for quantile in stats.iter_quantiles(arr.copy(), qs, method="sort")
-    ]
-    numpy = [
-        quantile for quantile in stats.iter_quantiles(arr.copy(), qs, method="numpy")
-    ]
+    sort = [quantile for quantile in stats.iter_quantiles(arr.copy(), qs, method="sort")]
+    numpy = [quantile for quantile in stats.iter_quantiles(arr.copy(), qs, method="numpy")]
     assert np.all(np.isclose(sort, numpy, equal_nan=True))
diff --git a/tests/thermo/test_thermo_array.py b/tests/thermo/test_thermo_array.py
index aa6a792..a87fcf8 100644
--- a/tests/thermo/test_thermo_array.py
+++ b/tests/thermo/test_thermo_array.py
@@ -10,6 +10,7 @@
 import os
 
 import numpy as np
+import pytest
 
 from earthkit.meteo import thermo
 
@@ -302,9 +303,7 @@ def test_saturation_specific_humidity_slope():
     )
 
     for phase in phases:
-        svp = thermo.array.saturation_specific_humidity_slope(
-            d["t"], d["p"], phase=phase
-        )
+        svp = thermo.array.saturation_specific_humidity_slope(d["t"], d["p"], phase=phase)
         np.testing.assert_allclose(svp, d[phase])
 
     v_ref = np.array([np.nan])
@@ -322,7 +321,7 @@ def test_saturation_specific_humidity_slope():
         np.testing.assert_allclose(svp, v_ref[i])
 
 
-def test_temperature_from_saturation_vapour_pressure():
+def test_temperature_from_saturation_vapour_pressure_1():
     ref_file = "sat_vp.csv"
     d = np.genfromtxt(
         data_file(ref_file),
@@ -331,7 +330,32 @@ def test_temperature_from_saturation_vapour_pressure():
     )
 
     t = thermo.array.temperature_from_saturation_vapour_pressure(d["water"])
-    np.testing.assert_allclose(t, d["t"])
+    assert np.allclose(t, d["t"], equal_nan=True)
+
+
+@pytest.mark.parametrize(
+    "es,kwargs,expected_values",
+    [
+        (4.2, {}, 219.7796336743947),
+        ([4.2, 0, 0.001, np.nan], {"eps": 1e-2, "out": 100}, [219.7796336743947, 100, 100, np.nan]),
+        ([4.2, 0, 0.001, np.nan], {"eps": 1e-2, "out": np.nan}, [219.7796336743947, np.nan, np.nan, np.nan]),
+        (0, {}, np.nan),
+        (0.001, {"eps": 1e-2, "out": 100}, 100.0),
+        (0.001, {"eps": 1e-2, "out": np.nan}, np.nan),
+    ],
+)
+def test_temperature_from_saturation_vapour_pressure_2(es, kwargs, expected_values):
+
+    multi = isinstance(es, list)
+    if multi:
+        es = np.array(es)
+        expected_values = np.array(expected_values)
+
+    t = thermo.array.temperature_from_saturation_vapour_pressure(es, **kwargs)
+    if multi:
+        np.testing.assert_allclose(t, expected_values, equal_nan=True)
+    else:
+        assert np.isclose(t, expected_values, equal_nan=True)
 
 
 def test_relative_humidity_from_dewpoint():
@@ -355,9 +379,7 @@ def test_relative_humidity_from_dewpoint():
 
 
 def test_relative_humidity_from_specific_humidity():
-    t = thermo.array.celsius_to_kelvin(
-        np.array([-29.2884, -14.4118, -5.9235, 9.72339, 18.4514])
-    )
+    t = thermo.array.celsius_to_kelvin(np.array([-29.2884, -14.4118, -5.9235, 9.72339, 18.4514]))
     p = np.array([300, 400, 500, 700, 850]) * 100.0
     q = np.array(
         [
@@ -401,9 +423,7 @@ def test_specific_humidity_from_dewpoint():
 
 
 def test_specific_humidity_from_relative_humidity():
-    t = thermo.array.celsius_to_kelvin(
-        np.array([-29.2884, -14.4118, -5.9235, 9.72339, 18.4514])
-    )
+    t = thermo.array.celsius_to_kelvin(np.array([-29.2884, -14.4118, -5.9235, 9.72339, 18.4514]))
     p = np.array([300, 400, 500, 700, 850]) * 100.0
     r = np.array(
         [
@@ -427,43 +447,122 @@ def test_specific_humidity_from_relative_humidity():
     np.testing.assert_allclose(q, v_ref)
 
 
-def test_dewpoint_from_relative_humidity():
-    t = thermo.array.celsius_to_kelvin(np.array([20.0, 20, 0, 35, 5, -15, 25]))
-    r = np.array(
-        [
-            100.0000000000,
-            52.5224541378,
-            46.8714823296,
-            84.5391163313,
-            21.9244774232,
-            46.1081101229,
-            15.4779832381,
-        ]
-    )
-    v_ref = thermo.array.celsius_to_kelvin(np.array([20.0, 10, -10, 32, -15, -24, -3]))
+@pytest.mark.parametrize(
+    "t,r,kwargs,expected_values",
+    [
+        (
+            [20.0, 20, 0, 35, 5, -15, 25, 25],
+            [
+                100.0000000000,
+                52.5224541378,
+                46.8714823296,
+                84.5391163313,
+                21.9244774232,
+                46.1081101229,
+                15.4779832381,
+                0,
+            ],
+            {},
+            [20.0, 10, -10, 32, -15, -24, -3, np.nan],
+        ),
+        (
+            [20.0, 20.0, 20.0],
+            [
+                52.5224541378,
+                0.0,
+                0.000001,
+            ],
+            {"eps": 1e-3, "out": thermo.array.celsius_to_kelvin(100)},
+            [10, 100, 100],
+        ),
+        (
+            [20.0, 20.0, 20.0],
+            [
+                52.5224541378,
+                0.0,
+                0.000001,
+            ],
+            {"eps": 1e-3, "out": np.nan},
+            [10, np.nan, np.nan],
+        ),
+        (20.0, 52.5224541378, {}, 10.0),
+        (20.0, 0.0, {"eps": 1e-3, "out": thermo.array.celsius_to_kelvin(100)}, 100),
+        (20.0, 0.000001, {"eps": 1e-3, "out": thermo.array.celsius_to_kelvin(100)}, 100),
+        (20.0, 0.0, {"eps": 1e-3, "out": np.nan}, np.nan),
+        (20, 0.000001, {"eps": 1e-3, "out": np.nan}, np.nan),
+    ],
+)
+def test_dewpoint_from_relative_humidity(t, r, kwargs, expected_values):
     # reference was tested with an online relhum calculator at:
     # https://bmcnoldy.rsmas.miami.edu/Humidity.html
-    td = thermo.array.dewpoint_from_relative_humidity(t, r)
-    np.testing.assert_allclose(td, v_ref)
 
-
-def test_dewpoint_from_specific_humidity():
-    p = np.array([967.5085, 936.3775, 872.248, 756.1647, 649.157, 422.4207]) * 100
-    q = np.array(
-        [
-            0.0169461501,
-            0.0155840075,
-            0.0134912382,
-            0.0083409720,
-            0.0057268584,
-            0.0025150791,
-        ]
-    )
-    v_ref = thermo.array.celsius_to_kelvin(
-        np.array([21.78907, 19.90885, 16.50236, 7.104064, -0.3548709, -16.37916])
-    )
-    td = thermo.array.dewpoint_from_specific_humidity(q, p)
-    np.testing.assert_allclose(td, v_ref)
+    multi = isinstance(t, list)
+    if multi:
+        t = np.array(t)
+        r = np.array(r)
+        expected_values = np.array(expected_values)
+
+    t = thermo.array.celsius_to_kelvin(t)
+    v_ref = thermo.array.celsius_to_kelvin(expected_values)
+
+    td = thermo.array.dewpoint_from_relative_humidity(t, r, **kwargs)
+    if multi:
+        assert np.allclose(td, v_ref, equal_nan=True)
+    else:
+        assert np.isclose(td, v_ref, equal_nan=True)
+
+
+@pytest.mark.parametrize(
+    "q,p,kwargs,expected_values",
+    [
+        (
+            [0.0169461501, 0.0155840075, 0.0134912382, 0.0083409720, 0.0057268584, 0.0025150791, 0],
+            [967.5085, 936.3775, 872.248, 756.1647, 649.157, 422.4207, 422.4207],
+            {},
+            [21.78907, 19.90885, 16.50236, 7.104064, -0.3548709, -16.37916, np.nan],
+        ),
+        (
+            [
+                0.0169461501,
+                0.0,
+                0.000001,
+            ],
+            [967.5085, 967.5085, 967.5085],
+            {"eps": 1e-3, "out": thermo.array.celsius_to_kelvin(100)},
+            [21.78907, 100, 100],
+        ),
+        (
+            [
+                0.0169461501,
+                0.0,
+                0.000001,
+            ],
+            [967.5085, 967.5085, 967.5085],
+            {"eps": 1e-3, "out": np.nan},
+            [21.78907, np.nan, np.nan],
+        ),
+        (0.0169461501, 967.508, {}, 21.78907),
+        (0.0, 967.5085, {"eps": 1e-3, "out": thermo.array.celsius_to_kelvin(100)}, 100),
+        (0.000001, 967.5085, {"eps": 1e-3, "out": thermo.array.celsius_to_kelvin(100)}, 100),
+        (0.0, 967.5085, {"eps": 1e-3, "out": np.nan}, np.nan),
+        (0.000001, 967.5085, {"eps": 1e-3, "out": np.nan}, np.nan),
+    ],
+)
+def test_dewpoint_from_specific_humidity(q, p, kwargs, expected_values):
+    multi = isinstance(q, list)
+    if multi:
+        q = np.array(q)
+        p = np.array(p)
+        expected_values = np.array(expected_values)
+
+    p = p * 100.0
+    v_ref = thermo.array.celsius_to_kelvin(expected_values)
+
+    td = thermo.array.dewpoint_from_specific_humidity(q, p, **kwargs)
+    if multi:
+        assert np.allclose(td, v_ref, equal_nan=True)
+    else:
+        assert np.isclose(td, v_ref, equal_nan=True)
 
 
 def test_virtual_temperature():
@@ -638,9 +737,7 @@ def test_temperature_on_moist_adiabat():
             pt = thermo.array.temperature_on_moist_adiabat(
                 ref["ept"], ref["p"], ept_method=m_ept, t_method=m_t
             )
-            np.testing.assert_allclose(
-                pt, ref[f"{m_ept}_{m_t}"], err_msg=f"method={m_ept}_{m_t}"
-            )
+            np.testing.assert_allclose(pt, ref[f"{m_ept}_{m_t}"], err_msg=f"method={m_ept}_{m_t}")
 
 
 def test_wet_bulb_temperature():
diff --git a/tests/wind/test_wind.py b/tests/wind/test_wind.py
index 09f1b1d..3648065 100644
--- a/tests/wind/test_wind.py
+++ b/tests/wind/test_wind.py
@@ -42,9 +42,7 @@ def test_direction():
     # meteo
     u = np.array([0, 1, 1, 1, 0, -1, -1, -1, 0, np.nan, 1, np.nan])
     v = np.array([1, 1, 0, -1, -1, -1, 0, 1, 0, 1, np.nan, np.nan])
-    v_ref = np.array(
-        [180.0, 225, 270, 315, 0, 45, 90, 135, 270, np.nan, np.nan, np.nan]
-    )
+    v_ref = np.array([180.0, 225, 270, 315, 0, 45, 90, 135, 270, np.nan, np.nan, np.nan])
     d = wind.direction(u, v, convention="meteo")
     np.testing.assert_allclose(d, v_ref)
 
@@ -122,9 +120,7 @@ def test_xy_to_polar():
             np.nan,
         ]
     )
-    d_ref = np.array(
-        [180.0, 225, 270, 315, 0, 45, 90, 135, 270, np.nan, np.nan, np.nan]
-    )
+    d_ref = np.array([180.0, 225, 270, 315, 0, 45, 90, 135, 270, np.nan, np.nan, np.nan])
     sp, d = wind.xy_to_polar(u, v)
     np.testing.assert_allclose(sp, sp_ref)
     np.testing.assert_allclose(d, d_ref)
@@ -172,12 +168,8 @@ def test_coriolis():
 
 
 def test_windrose():
-    sp = np.array(
-        [3.5, 1, 1.1, 2.1, 0.1, 0.0, 2.4, 1.9, 1.7, 3.9, 3.1, 2.1, np.nan, np.nan]
-    )
-    d = np.array(
-        [1.0, 29, 31, 93.0, 121, 171, 189, 245, 240.11, 311, 359.1, np.nan, 11, np.nan]
-    )
+    sp = np.array([3.5, 1, 1.1, 2.1, 0.1, 0.0, 2.4, 1.9, 1.7, 3.9, 3.1, 2.1, np.nan, np.nan])
+    d = np.array([1.0, 29, 31, 93.0, 121, 171, 189, 245, 240.11, 311, 359.1, np.nan, 11, np.nan])
     sp_bins = [0, 1, 2, 3, 4]
 
     # count