Repository: chrisroadmap/climateforcing

Issue:
Check APRP against Mark Zelinka
- [ ] understand differences (https://github.com/mzelinka/aprp/blob/master/code/python/demonstrate_APRP.ipynb)
- [ ] correct my code if buggy

Discussion:

Patch:
diff --git a/.github/workflows/checks.yml b/.github/workflows/checks.yml
index f3b092d..15f218c 100644
--- a/.github/workflows/checks.yml
+++ b/.github/workflows/checks.yml
@@ -10,7 +10,7 @@ jobs:
     runs-on: ubuntu-latest
     strategy:
       matrix:
-        python-version: [3.6, 3.7, 3.8, 3.9]
+        python-version: ["3.7", "3.8", "3.9", "3.10", "3.11"]
 
     steps:
     - name: Checkout repository
@@ -51,7 +51,7 @@ jobs:
     runs-on: ubuntu-latest
     strategy:
       matrix:
-        python-version: [3.6, 3.7, 3.8, 3.9]
+        python-version: ["3.7", "3.8", "3.9", "3.10", "3.11"]
 
     steps:
     - name: Checkout repository
diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index e29866f..5f2f30b 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -13,6 +13,11 @@ Following the recommendations of https://keepachangelog.com/en/1.0.0/
 master
 ------
 
+v0.3.0
+------
+- fixed: `aprp` was slightly off and has been corrected (`#33 <https://github.com/chrisroadmap/climateforcing/pull/33>`_)
+- changed: support for python 3.6 dropped, support for python 3.10 and 3.11 added
+
 v0.2.3
 ------
 - added: option to specify directory in the path of `check_and_download` and to automatically create a new directory if it doesn't exist (`#30 <https://github.com/chrisroadmap/climateforcing/pull/30>`_)
diff --git a/setup.py b/setup.py
index b631ccf..cc3c7f4 100644
--- a/setup.py
+++ b/setup.py
@@ -25,11 +25,11 @@ REQUIREMENTS_DOCS = ["sphinx>=1.4", "sphinx_rtd_theme"]
 
 REQUIREMENTS_DEV = [
     *[
-        "black==19.10b0",
+        "black",
         "flake8",
-        "isort>5",
+        "isort",
         "pydocstyle",
-        "pylint>=2.4.4",
+        "pylint",
         "readme-renderer",
     ],
     *REQUIREMENTS_DEPLOY,
@@ -67,10 +67,11 @@ setup(
         "License :: OSI Approved :: BSD License",
         "Intended Audience :: Developers",
         "Operating System :: OS Independent",
-        "Programming Language :: Python :: 3.6",
         "Programming Language :: Python :: 3.7",
         "Programming Language :: Python :: 3.8",
         "Programming Language :: Python :: 3.9",
+        "Programming Language :: Python :: 3.10",
+        "Programming Language :: Python :: 3.11",
     ],
     keywords=["climate", "tools"],
     packages=PACKAGES,
diff --git a/src/climateforcing/aprp/aprp.py b/src/climateforcing/aprp/aprp.py
index 99133ca..53b3e88 100644
--- a/src/climateforcing/aprp/aprp.py
+++ b/src/climateforcing/aprp/aprp.py
@@ -1,4 +1,14 @@
 """APRP: the approximate partial radiative perturbation calculation."""
+
+# Modified Feb 2023 by M. Zelinka to do the following:
+# 1) Multiply albedo sensitivities by downwelling SW rather than net (down minus up) SW
+# 2) Fix the manner of estimating overcast albedo sensitivity to cloud and non-cloud
+#    gamma and mu
+# 3) Separate the forward and backward albedo sensitivity calculations
+
+# These modifications lead to negligible TOA SW residuals and to
+# perfect agreement with the results of Mark's code (https://github.com/mzelinka/aprp)
+
 import glob
 import warnings
 
@@ -25,10 +35,8 @@ def cloud_radiative_effect(base, pert):
     base, pert: dict of array_like
         CMIP diagnostics required to calculate longwave cloud radiative effect. The
         dicts should contain two keys:
-
         rlut    : top-of-atmosphere outgoing longwave flux
         rlutcs  : top-of-atmosphere longwave flux assuming clear sky
-
     Returns
     -------
     erfari_lw, erfaci_lw : array_like
@@ -49,8 +57,14 @@ def cloud_radiative_effect(base, pert):
     return erfari_lw, erfaci_lw
 
 
-def aprp(  # pylint: disable=too-many-arguments,too-many-locals,too-many-statements  # noqa: E501
-    base, pert, longwave=False, breakdown=False, cs_threshold=0.02, clt_percent=True,
+def aprp(  # pylint: disable=too-many-arguments,too-many-locals,too-many-statements,too-many-branches  # noqa: E501
+    base,
+    pert,
+    longwave=False,
+    breakdown=False,
+    cs_threshold=0.02,
+    rsdt_threshold=0.1,
+    clt_percent=True,
 ):
     """Approximate Partial Raditive Perturbation calculation.
 
@@ -94,6 +108,9 @@ def aprp(  # pylint: disable=too-many-arguments,too-many-locals,too-many-stateme
             flux to zero. It is recommended to use a small positive value, as the
             cloud fraction appears in the denominator of the calculation. Taken from
             Mark Zelinka's implementation.
+        rsdt_threshold : float, default=0.1
+            set incoming radiation to zero below a certain level. A small positive value
+            is recommended to reduce residuals of the SW component compared to TOA ERF.
         clt_percent : bool, default=True
             is cloud fraction from base and pert in percent (True) or 0-1 scale
             (False)
@@ -101,8 +118,8 @@ def aprp(  # pylint: disable=too-many-arguments,too-many-locals,too-many-stateme
     Returns
     -------
         central[, forward, reverse] : dict of array_like
-            Components of APRP as defined by equation A2 of [1]_.
 
+            Components of APRP as defined by equation A2 of [1]_.
             dict keys are 't1', 't2', ..., 't9' where tX is the corresponding term in
             eq. A2.
 
@@ -118,7 +135,6 @@ def aprp(  # pylint: disable=too-many-arguments,too-many-locals,too-many-stateme
 
             though note these only make sense if you are calculating aerosol forcing.
             The cloud fraction adjustment component of ERFaci is t9.
-
             if longwave is True, central also contains 'ERFariLW' and 'ERFaciLW' as
             calculated from the cloud radiative effect.
 
@@ -204,7 +220,7 @@ def aprp(  # pylint: disable=too-many-arguments,too-many-locals,too-many-stateme
         rsutoc = 0.5 * (pert["rsutoc"] + base["rsutoc"])
         rsutcs = 0.5 * (pert["rsutcs"] + base["rsutcs"])
         delta_clt = pert["clt"] - base["clt"]
-        clt = 0.5 * (pert["clt"] + base["clt"])
+        rsdt = 0.5 * (base["rsdt"] + pert["rsdt"])
 
         base["rsusoc"] = (base["rsus"] - (1 - base["clt"]) * (base["rsuscs"])) / base[
             "clt"
@@ -300,193 +316,137 @@ def aprp(  # pylint: disable=too-many-arguments,too-many-locals,too-many-stateme
         mu_base = (mu_oc_base) / mu_cs_base
         mu_base[~np.isfinite(mu_base)] = 0.0
 
-    daoc_dacld = 0.5 * (
-        (
-            _planetary_albedo(mu_oc_base, gamma_oc_base, alpha_oc_pert)
-            - _planetary_albedo(mu_oc_base, gamma_oc_base, alpha_oc_base)
-        )
-        + (
-            _planetary_albedo(mu_oc_pert, gamma_oc_pert, alpha_oc_pert)
-            - (_planetary_albedo(mu_oc_pert, gamma_oc_pert, alpha_oc_base))
-        )
+    aoc_base = _planetary_albedo(mu_oc_base, gamma_oc_base, alpha_oc_base)
+    aoc_pert = _planetary_albedo(mu_oc_pert, gamma_oc_pert, alpha_oc_pert)
+    acs_base = _planetary_albedo(mu_cs_base, gamma_cs_base, alpha_cs_base)
+    acs_pert = _planetary_albedo(mu_cs_pert, gamma_cs_pert, alpha_cs_pert)
+
+    daoc_dacld_fwd = (
+        _planetary_albedo(mu_oc_base, gamma_oc_base, alpha_oc_pert) - aoc_base
     )
-    daoc_dmcld = 0.5 * (
-        (
-            _planetary_albedo(mu_pert, gamma_base, alpha_oc_base)
-            - _planetary_albedo(mu_base, gamma_base, alpha_oc_base)
-        )
-        + (
-            _planetary_albedo(mu_pert, gamma_pert, alpha_oc_pert)
-            - (_planetary_albedo(mu_base, gamma_pert, alpha_oc_pert))
-        )
+    daoc_dacld_bwd = aoc_pert - _planetary_albedo(
+        mu_oc_pert, gamma_oc_pert, alpha_oc_base
     )
-    daoc_dgcld = 0.5 * (
-        (
-            _planetary_albedo(mu_base, gamma_pert, alpha_oc_base)
-            - _planetary_albedo(mu_base, gamma_base, alpha_oc_base)
-        )
-        + (
-            _planetary_albedo(mu_pert, gamma_pert, alpha_oc_pert)
-            - (_planetary_albedo(mu_pert, gamma_base, alpha_oc_pert))
-        )
+
+    dacs_daclr_fwd = (
+        _planetary_albedo(mu_cs_base, gamma_cs_base, alpha_cs_pert) - acs_base
+    )
+    dacs_daclr_bwd = acs_pert - _planetary_albedo(
+        mu_cs_pert, gamma_cs_pert, alpha_cs_base
     )
 
-    daoc_dmaer = 0.5 * (
-        (
-            _planetary_albedo(mu_cs_pert, gamma_oc_base, alpha_oc_base)
-            - _planetary_albedo(mu_cs_base, gamma_oc_base, alpha_oc_base)
-        )
-        + (
-            _planetary_albedo(mu_cs_pert, gamma_oc_pert, alpha_oc_pert)
-            - (_planetary_albedo(mu_cs_base, gamma_oc_pert, alpha_oc_pert))
-        )
+    dacs_dmaer_fwd = (
+        _planetary_albedo(mu_cs_pert, gamma_cs_base, alpha_cs_base) - acs_base
     )
-    daoc_dgaer = 0.5 * (
-        (
-            _planetary_albedo(mu_oc_base, gamma_cs_pert, alpha_oc_base)
-            - _planetary_albedo(mu_oc_base, gamma_cs_base, alpha_oc_base)
-        )
-        + (
-            _planetary_albedo(mu_oc_pert, gamma_cs_pert, alpha_oc_pert)
-            - (_planetary_albedo(mu_oc_pert, gamma_cs_base, alpha_oc_pert))
-        )
+    dacs_dmaer_bwd = acs_pert - _planetary_albedo(
+        mu_cs_base, gamma_cs_pert, alpha_cs_pert
     )
 
-    dacs_daclr = 0.5 * (
-        (
-            _planetary_albedo(mu_cs_base, gamma_cs_base, alpha_cs_pert)
-            - _planetary_albedo(mu_cs_base, gamma_cs_base, alpha_cs_base)
-        )
-        + (
-            _planetary_albedo(mu_cs_pert, gamma_cs_pert, alpha_cs_pert)
-            - (_planetary_albedo(mu_cs_pert, gamma_cs_pert, alpha_cs_base))
-        )
+    dacs_dgaer_fwd = (
+        _planetary_albedo(mu_cs_base, gamma_cs_pert, alpha_cs_base) - acs_base
     )
-    dacs_dmaer = 0.5 * (
-        (
-            _planetary_albedo(mu_cs_pert, gamma_cs_base, alpha_cs_base)
-            - _planetary_albedo(mu_cs_base, gamma_cs_base, alpha_cs_base)
-        )
-        + (
-            _planetary_albedo(mu_cs_pert, gamma_cs_pert, alpha_cs_pert)
-            - (_planetary_albedo(mu_cs_base, gamma_cs_pert, alpha_cs_pert))
-        )
+    dacs_dgaer_bwd = acs_pert - _planetary_albedo(
+        mu_cs_pert, gamma_cs_base, alpha_cs_pert
     )
-    dacs_dgaer = 0.5 * (
-        (
-            _planetary_albedo(mu_cs_base, gamma_cs_pert, alpha_cs_base)
-            - _planetary_albedo(mu_cs_base, gamma_cs_base, alpha_cs_base)
-        )
-        + (
-            _planetary_albedo(mu_cs_pert, gamma_cs_pert, alpha_cs_pert)
-            - (_planetary_albedo(mu_cs_pert, gamma_cs_base, alpha_cs_pert))
-        )
+
+    # Need isolated effect of mu_cs on mu_oc holding mu_cloud fixed, and vice versa:
+    new_mu_oc_pert = mu_cs_pert * mu_base  # Eq. 14
+    new_mu_oc_base = mu_cs_base * mu_pert  # Eq. 14
+    daoc_dmcld_fwd = (
+        _planetary_albedo(new_mu_oc_base, gamma_oc_base, alpha_oc_base) - aoc_base
+    )
+    daoc_dmcld_bwd = aoc_pert - _planetary_albedo(
+        new_mu_oc_pert, gamma_oc_pert, alpha_oc_pert
+    )
+    daoc_dmaer_fwd = (
+        _planetary_albedo(new_mu_oc_pert, gamma_oc_base, alpha_oc_base) - aoc_base
+    )
+    daoc_dmaer_bwd = aoc_pert - _planetary_albedo(
+        new_mu_oc_base, gamma_oc_pert, alpha_oc_pert
     )
 
-    base["rsnt"] = base["rsdt"] - base["rsut"]
-    base["rsntcs"] = base["rsdt"] - base["rsutcs"]
-    pert["rsnt"] = pert["rsdt"] - pert["rsut"]
-    pert["rsntcs"] = pert["rsdt"] - pert["rsutcs"]
-    pert["rsntoc"] = pert["rsdt"] - pert["rsutoc"]
-    base["rsntoc"] = base["rsdt"] - base["rsutoc"]
+    # Need isolated effect of gamma_cs on gamma_oc holding gamma_cloud fixed, and vice
+    # versa:
+    new_gamma_oc_pert = 1 - (1 - gamma_cs_pert) * (1 - gamma_base)  # Eq. 13
+    new_gamma_oc_base = 1 - (1 - gamma_cs_base) * (1 - gamma_pert)  # Eq. 13
+    daoc_dgcld_fwd = (
+        _planetary_albedo(mu_oc_base, new_gamma_oc_base, alpha_oc_base) - aoc_base
+    )
+    daoc_dgcld_bwd = aoc_pert - _planetary_albedo(
+        mu_oc_pert, new_gamma_oc_pert, alpha_oc_pert
+    )
+    daoc_dgaer_fwd = (
+        _planetary_albedo(mu_oc_base, new_gamma_oc_pert, alpha_oc_base) - aoc_base
+    )
+    daoc_dgaer_bwd = aoc_pert - _planetary_albedo(
+        mu_oc_pert, new_gamma_oc_base, alpha_oc_pert
+    )
 
     # t1 to t9 are the coefficients of equation A2 in Zelinka et al., 2014
     forward = {}
     reverse = {}
     central = {}
-    forward["t1"] = -base["rsntcs"] * (1 - clt) * dacs_daclr
-    forward["t2"] = -base["rsntcs"] * (1 - clt) * dacs_dgaer
-    forward["t3"] = -base["rsntcs"] * (1 - clt) * dacs_dmaer
-    forward["t4"] = -base["rsnt"] * clt * (daoc_dacld)
-    forward["t5"] = -base["rsnt"] * clt * (daoc_dgaer)
-    forward["t6"] = -base["rsnt"] * clt * (daoc_dmaer)
-    forward["t7"] = -base["rsnt"] * clt * (daoc_dgcld)
-    forward["t8"] = -base["rsnt"] * clt * (daoc_dmcld)
+    forward["t1"] = -base["rsdt"] * (1 - base["clt"]) * dacs_daclr_fwd
+    forward["t2"] = -base["rsdt"] * (1 - base["clt"]) * dacs_dgaer_fwd
+    forward["t3"] = -base["rsdt"] * (1 - base["clt"]) * dacs_dmaer_fwd
+    forward["t4"] = -base["rsdt"] * base["clt"] * (daoc_dacld_fwd)
+    forward["t5"] = -base["rsdt"] * base["clt"] * (daoc_dgaer_fwd)
+    forward["t6"] = -base["rsdt"] * base["clt"] * (daoc_dmaer_fwd)
+    forward["t7"] = -base["rsdt"] * base["clt"] * (daoc_dgcld_fwd)
+    forward["t8"] = -base["rsdt"] * base["clt"] * (daoc_dmcld_fwd)
     forward["t9"] = -delta_clt * (rsutoc - rsutcs)
-    forward["t2_clr"] = -base["rsntcs"] * dacs_dgaer
-    forward["t3_clr"] = -base["rsntcs"] * dacs_dmaer
+    forward["t2_clr"] = -base["rsdt"] * dacs_dgaer_fwd
+    forward["t3_clr"] = -base["rsdt"] * dacs_dmaer_fwd
 
     # set thresholds
     # TODO: can we avoid a hard cloud fraction threshold here?
-    forward["t4"] = np.where(
-        np.logical_or(base["clt"] < cs_threshold, pert["clt"] < cs_threshold),
-        0.0,
-        forward["t4"],
-    )
-    forward["t5"] = np.where(
-        np.logical_or(base["clt"] < cs_threshold, pert["clt"] < cs_threshold),
-        0.0,
-        forward["t5"],
-    )
-    forward["t6"] = np.where(
-        np.logical_or(base["clt"] < cs_threshold, pert["clt"] < cs_threshold),
-        0.0,
-        forward["t6"],
-    )
-    forward["t7"] = np.where(
-        np.logical_or(base["clt"] < cs_threshold, pert["clt"] < cs_threshold),
-        0.0,
-        forward["t7"],
-    )
-    forward["t8"] = np.where(
-        np.logical_or(base["clt"] < cs_threshold, pert["clt"] < cs_threshold),
-        0.0,
-        forward["t8"],
-    )
-    forward["t9"] = np.where(
-        np.logical_or(base["clt"] < cs_threshold, pert["clt"] < cs_threshold),
-        0.0,
-        forward["t9"],
-    )
+    for term in ["t4", "t5", "t6", "t7", "t8", "t9"]:
+        forward[term] = np.where(
+            np.logical_or(base["clt"] < cs_threshold, pert["clt"] < cs_threshold),
+            0.0,
+            forward[term],
+        )
+
+    # set fields to zero when incoming solar radiation is zero
+    for term in ["t1", "t2", "t3", "t4", "t5", "t6", "t7", "t8", "t9"]:
+        forward[term] = np.where(
+            rsdt < rsdt_threshold,
+            0.0,
+            forward[term],
+        )
 
     forward["ERFariSWclr"] = forward["t2_clr"] + forward["t3_clr"]
     forward["ERFariSW"] = forward["t2"] + forward["t3"] + forward["t5"] + forward["t6"]
     forward["ERFaciSW"] = forward["t7"] + forward["t8"] + forward["t9"]
     forward["albedo"] = forward["t1"] + forward["t4"]
 
-    reverse["t1"] = -pert["rsntcs"] * (1 - clt) * dacs_daclr
-    reverse["t2"] = -pert["rsntcs"] * (1 - clt) * dacs_dgaer
-    reverse["t3"] = -pert["rsntcs"] * (1 - clt) * dacs_dmaer
-    reverse["t4"] = -pert["rsnt"] * clt * (daoc_dacld)
-    reverse["t5"] = -pert["rsnt"] * clt * (daoc_dgaer)
-    reverse["t6"] = -pert["rsnt"] * clt * (daoc_dmaer)
-    reverse["t7"] = -pert["rsnt"] * clt * (daoc_dgcld)
-    reverse["t8"] = -pert["rsnt"] * clt * (daoc_dmcld)
+    reverse["t1"] = -pert["rsdt"] * (1 - pert["clt"]) * dacs_daclr_bwd
+    reverse["t2"] = -pert["rsdt"] * (1 - pert["clt"]) * dacs_dgaer_bwd
+    reverse["t3"] = -pert["rsdt"] * (1 - pert["clt"]) * dacs_dmaer_bwd
+    reverse["t4"] = -pert["rsdt"] * pert["clt"] * (daoc_dacld_bwd)
+    reverse["t5"] = -pert["rsdt"] * pert["clt"] * (daoc_dgaer_bwd)
+    reverse["t6"] = -pert["rsdt"] * pert["clt"] * (daoc_dmaer_bwd)
+    reverse["t7"] = -pert["rsdt"] * pert["clt"] * (daoc_dgcld_bwd)
+    reverse["t8"] = -pert["rsdt"] * pert["clt"] * (daoc_dmcld_bwd)
     reverse["t9"] = -delta_clt * (rsutoc - rsutcs)
-    reverse["t2_clr"] = -pert["rsntcs"] * dacs_dgaer
-    reverse["t3_clr"] = -pert["rsntcs"] * dacs_dmaer
+    reverse["t2_clr"] = -pert["rsdt"] * dacs_dgaer_bwd
+    reverse["t3_clr"] = -pert["rsdt"] * dacs_dmaer_bwd
 
     # set thresholds
-    reverse["t4"] = np.where(
-        np.logical_or(base["clt"] < cs_threshold, pert["clt"] < cs_threshold),
-        0.0,
-        reverse["t4"],
-    )
-    reverse["t5"] = np.where(
-        np.logical_or(base["clt"] < cs_threshold, pert["clt"] < cs_threshold),
-        0.0,
-        reverse["t5"],
-    )
-    reverse["t6"] = np.where(
-        np.logical_or(base["clt"] < cs_threshold, pert["clt"] < cs_threshold),
-        0.0,
-        reverse["t6"],
-    )
-    reverse["t7"] = np.where(
-        np.logical_or(base["clt"] < cs_threshold, pert["clt"] < cs_threshold),
-        0.0,
-        reverse["t7"],
-    )
-    reverse["t8"] = np.where(
-        np.logical_or(base["clt"] < cs_threshold, pert["clt"] < cs_threshold),
-        0.0,
-        reverse["t8"],
-    )
-    reverse["t9"] = np.where(
-        np.logical_or(base["clt"] < cs_threshold, pert["clt"] < cs_threshold),
-        0.0,
-        reverse["t9"],
-    )
+    for term in ["t4", "t5", "t6", "t7", "t8", "t9"]:
+        reverse[term] = np.where(
+            np.logical_or(base["clt"] < cs_threshold, pert["clt"] < cs_threshold),
+            0.0,
+            reverse[term],
+        )
+
+    # set fields to zero when incoming solar radiation is zero
+    for term in ["t1", "t2", "t3", "t4", "t5", "t6", "t7", "t8", "t9"]:
+        reverse[term] = np.where(
+            rsdt < rsdt_threshold,
+            0.0,
+            reverse[term],
+        )
 
     reverse["ERFariSWclr"] = reverse["t2_clr"] + reverse["t3_clr"]
     reverse["ERFariSW"] = reverse["t2"] + reverse["t3"] + reverse["t5"] + reverse["t6"]
diff --git a/src/climateforcing/atmos/humidity.py b/src/climateforcing/atmos/humidity.py
index f4c580e..b4eebae 100644
--- a/src/climateforcing/atmos/humidity.py
+++ b/src/climateforcing/atmos/humidity.py
@@ -66,7 +66,10 @@ def calc_saturation_mixing_ratio(air_temperature, pressure):
 
 
 def specific_to_relative(
-    specific_humidity, pressure=101325, air_temperature=288.15, rh_percent=False,
+    specific_humidity,
+    pressure=101325,
+    air_temperature=288.15,
+    rh_percent=False,
 ):
     """Convert specific humidity to relative humidity.
 
@@ -94,7 +97,10 @@ def specific_to_relative(
 
 
 def relative_to_specific(
-    relative_humidity, pressure=101325, air_temperature=288.15, rh_percent=False,
+    relative_humidity,
+    pressure=101325,
+    air_temperature=288.15,
+    rh_percent=False,
 ):
     """Convert relative humidity to specific humidity.
 
diff --git a/src/climateforcing/solar/orbital_constants.py b/src/climateforcing/solar/orbital_constants.py
index 7bf1151..e86ceee 100644
--- a/src/climateforcing/solar/orbital_constants.py
+++ b/src/climateforcing/solar/orbital_constants.py
@@ -149,7 +149,9 @@ HELIOCENTRIC_LONGITUDE_COEFFS = {
         ]
     ),
     4: np.array([[114.0, 3.142, 0], [8, 4.13, 6283.08], [1, 3.84, 12566.15]]),
-    5: np.array([[1, 3.14, 0]],),
+    5: np.array(
+        [[1, 3.14, 0]],
+    ),
 }
 
 
@@ -163,7 +165,9 @@ HELIOCENTRIC_LATITUDE_COEFFS = {
             [32, 4, 1577.34],
         ]
     ),
-    1: np.array([[9, 3.9, 5507.55], [6, 1.73, 5223.69]],),
+    1: np.array(
+        [[9, 3.9, 5507.55], [6, 1.73, 5223.69]],
+    ),
 }
 
 
diff --git a/src/climateforcing/solar/solar_position.py b/src/climateforcing/solar/solar_position.py
index 7f83ba6..e3fa28b 100644
--- a/src/climateforcing/solar/solar_position.py
+++ b/src/climateforcing/solar/solar_position.py
@@ -95,7 +95,7 @@ def _term_sum(julian_date, coefficient_set):
                 + coefficient_set[icoeff][:, 2] * julian_millenium
             )
         )
-        result = result + 1e-8 * coeff * julian_millenium ** icoeff
+        result = result + 1e-8 * coeff * julian_millenium**icoeff
     return result
 
 
@@ -163,15 +163,15 @@ def _true_obliquity_of_ecliptic(julian_date):
     mean_obliquity = (
         84381.448
         - 4680.93 * julian_10ka
-        - 1.55 * julian_10ka ** 2
-        + 1999.25 * julian_10ka ** 3
-        - 51.38 * julian_10ka ** 4
-        - 249.67 * julian_10ka ** 5
-        - 39.05 * julian_10ka ** 6
-        + 7.12 * julian_10ka ** 7
-        + 27.87 * julian_10ka ** 8
-        + 5.79 * julian_10ka ** 9
-        + 2.45 * julian_10ka ** 10
+        - 1.55 * julian_10ka**2
+        + 1999.25 * julian_10ka**3
+        - 51.38 * julian_10ka**4
+        - 249.67 * julian_10ka**5
+        - 39.05 * julian_10ka**6
+        + 7.12 * julian_10ka**7
+        + 27.87 * julian_10ka**8
+        + 5.79 * julian_10ka**9
+        + 2.45 * julian_10ka**10
     ) / 3600
     _, nutation_obliquity = _nutation(julian_date)
     return mean_obliquity + nutation_obliquity
diff --git a/src/climateforcing/twolayermodel/__init__.py b/src/climateforcing/twolayermodel/__init__.py
index 2f6f28d..3f9f0c8 100644
--- a/src/climateforcing/twolayermodel/__init__.py
+++ b/src/climateforcing/twolayermodel/__init__.py
@@ -378,10 +378,10 @@ class TwoLayer:  # pylint: disable=too-few-public-methods,too-many-instance-attr
                 "adeep_s": adeep_s,
             }
         )
-        for k in self.params.keys():
-            if self.params[k] is None:
-                if k in self.default:
-                    self.params[k] = self.default[k]
+        for key, value in self.params.items():
+            if value is None:
+                if key in self.default:
+                    self.params[key] = self.default[key]
         self.afast = afast
         self.aslow = aslow
         self.amix_f = amix_f


Test patch:
diff --git a/tests/test_aprp.py b/tests/test_aprp.py
index 3f398e2..ab40b96 100644
--- a/tests/test_aprp.py
+++ b/tests/test_aprp.py
@@ -50,23 +50,23 @@ def test_aprp_access_esm1_5():
     # these values are not expected to be realistic global means as our test data is
     # a small area in the tropics
     EXPECTED_RESULT = {
-        "t1": 0.015690050209499505,
-        "t2": -0.11169482070528963,
-        "t3": 0.03628321688570137,
-        "t4": -0.016643288519753274,
-        "t5": -0.3646235907233642,
-        "t6": 0.2702892397012102,
-        "t7": -1.7754625968973907,
-        "t8": 0.11863184693029566,
+        "albedo": -0.004681121779074323,
+        "ERFariSW": -0.0337658957581616,
+        "ERFaciSW": 0.4121098817036564,
+        "ERFariLW": -0.8487686048639489,
+        "ERFaciLW": -3.321506047334925,
+        "t1": 0.017574093155832592,
+        "t2": -0.12323665520968012,
+        "t3": 0.041323059475741955,
+        "t4": -0.02225521556606039,
+        "t5": -0.32008254859568436,
+        "t6": 0.3682302585091456,
+        "t7": -1.7441234148989013,
+        "t8": 0.14377014708168268,
         "t9": 2.0124631258387526,
-        "t2_clr": -0.6004451771086404,
-        "t3_clr": 0.14946570231037004,
-        "ERFariSWclr": -0.45097949302816626,
-        "ERFariSW": -0.16974595857509908,
-        "ERFaciSW": 0.35563230710202826,
-        "albedo": -0.0009532384845432344,
-        "ERFariLW": -0.8487686048639488,
-        "ERFaciLW": -3.3215060473349243,
+        "t2_clr": -0.6718520381414649,
+        "t3_clr": 0.16847785214017594,
+        "ERFariSWclr": -0.5033741958574206,
     }
     BASEDIR = "tests/testdata/ACCESS-ESM1-5/piClim-control/"
     PERTDIR = "tests/testdata/ACCESS-ESM1-5/piClim-aer/"
@@ -74,13 +74,7 @@ def test_aprp_access_esm1_5():
     result_3d = aprp(base, pert, longwave=True)
     for key, value in result_3d.items():
         result_1d = global_mean(value, lat=lat, axis=1)
-        assert np.allclose(result_1d, EXPECTED_RESULT[key])
-    # should use conftest and this should be a separate test
-    central, forward, reverse = aprp(base, pert, longwave=True, breakdown=True)
-    for key in central:
-        assert ~np.allclose(central[key], forward[key])
-        assert ~np.allclose(forward[key], reverse[key])
-        assert ~np.allclose(reverse[key], central[key])
+        assert np.allclose(result_1d, EXPECTED_RESULT[key], atol=1e-5)
 
 
 def test_create_input_slice():
diff --git a/tests/test_utils.py b/tests/test_utils.py
index 80e6fc3..d9cbbea 100644
--- a/tests/test_utils.py
+++ b/tests/test_utils.py
@@ -13,18 +13,18 @@ def test_check_and_download():
     # checking both the operation of check_and_download and the capability to
     # make a new directory on the fly
     check_and_download(
-        "http://homepages.see.leeds.ac.uk/~mencsm/images/ta_q_kernel.png",
+        "https://docs.fairmodel.net/en/latest/_images/dimensions.png",
         "tests/newdirectory/test_image.png",
     )
     # clobber working?
     check_and_download(
-        "http://homepages.see.leeds.ac.uk/~mencsm/images/ta_q_kernel.png",
+        "https://docs.fairmodel.net/en/latest/_images/dimensions.png",
         "tests/newdirectory/test_image.png",
         clobber=True,
     )
     # directory specification working?
     check_and_download(
-        "http://homepages.see.leeds.ac.uk/~mencsm/images/ta_q_kernel.png",
+        "https://docs.fairmodel.net/en/latest/_images/dimensions.png",
         "tests/newdirectory",
         clobber=True,
     )