diff --git a/esmvaltool/diag_scripts/aerosols/aod_aeronet_assess.py b/esmvaltool/diag_scripts/aerosols/aod_aeronet_assess.py
index 27ab6b2714..3866e3c51a 100644
--- a/esmvaltool/diag_scripts/aerosols/aod_aeronet_assess.py
+++ b/esmvaltool/diag_scripts/aerosols/aod_aeronet_assess.py
@@ -11,10 +11,10 @@
 import matplotlib.pyplot as plt
 import numpy as np
 import scipy
-from aero_utils import add_bounds, extract_pt
 from matplotlib import colors, gridspec
 from numpy import ma
 
+from esmvaltool.diag_scripts.aerosols.aero_utils import add_bounds, extract_pt
 from esmvaltool.diag_scripts.shared import group_metadata, run_diagnostic
 from esmvaltool.diag_scripts.shared._base import get_plot_filename
 
diff --git a/esmvaltool/diag_scripts/climate_patterns/climate_patterns.py b/esmvaltool/diag_scripts/climate_patterns/climate_patterns.py
index 7fdb98a293..bc265cda00 100644
--- a/esmvaltool/diag_scripts/climate_patterns/climate_patterns.py
+++ b/esmvaltool/diag_scripts/climate_patterns/climate_patterns.py
@@ -28,24 +28,24 @@
 """
 
 import logging
-from pathlib import Path
 import os
+from pathlib import Path
 
 import iris
 import iris.coord_categorisation
 import iris.cube
 import numpy as np
 import sklearn.linear_model
-
-import sub_functions as sf
-from plotting import (
-    plot_timeseries,
-    plot_patterns
-)
 from esmvalcore.preprocessor import (
     area_statistics,
+    climate_statistics,
     extract_time,
-    climate_statistics
+)
+
+import esmvaltool.diag_scripts.climate_patterns.sub_functions as sf
+from esmvaltool.diag_scripts.climate_patterns.plotting import (
+    plot_patterns,
+    plot_timeseries,
 )
 from esmvaltool.diag_scripts.shared import run_diagnostic
 
diff --git a/esmvaltool/diag_scripts/ensclus/ens_anom.py b/esmvaltool/diag_scripts/ensclus/ens_anom.py
index a3b664d11a..ef9c303786 100644
--- a/esmvaltool/diag_scripts/ensclus/ens_anom.py
+++ b/esmvaltool/diag_scripts/ensclus/ens_anom.py
@@ -1,12 +1,19 @@
 """Computation of ensemble anomalies based on a desired value."""
 
 import os
+
 import numpy as np
 from scipy import stats
 
 # User-defined packages
-from read_netcdf import read_iris, save_n_2d_fields
-from sel_season_area import sel_area, sel_season
+from esmvaltool.diag_scripts.ensclus.read_netcdf import (
+    read_iris,
+    save_n_2d_fields,
+)
+from esmvaltool.diag_scripts.ensclus.sel_season_area import (
+    sel_area,
+    sel_season,
+)
 
 
 def ens_anom(filenames, dir_output, name_outputs, varname, numens, season,
diff --git a/esmvaltool/diag_scripts/ensclus/ens_eof_kmeans.py b/esmvaltool/diag_scripts/ensclus/ens_eof_kmeans.py
index 2be17195b0..940af71278 100644
--- a/esmvaltool/diag_scripts/ensclus/ens_eof_kmeans.py
+++ b/esmvaltool/diag_scripts/ensclus/ens_eof_kmeans.py
@@ -10,8 +10,8 @@
 from sklearn.cluster import KMeans
 
 # User-defined libraries
-from eof_tool import eof_computation
-from read_netcdf import read_n_2d_fields
+from esmvaltool.diag_scripts.ensclus.eof_tool import eof_computation
+from esmvaltool.diag_scripts.ensclus.read_netcdf import read_n_2d_fields
 
 
 def ens_eof_kmeans(dir_output, name_outputs, numens, numpcs, perc, numclus):
diff --git a/esmvaltool/diag_scripts/ensclus/ens_plots.py b/esmvaltool/diag_scripts/ensclus/ens_plots.py
index 55ad420b08..1b06acd671 100644
--- a/esmvaltool/diag_scripts/ensclus/ens_plots.py
+++ b/esmvaltool/diag_scripts/ensclus/ens_plots.py
@@ -3,12 +3,12 @@
 import math
 import os
 
+import cartopy.crs as ccrs
 import matplotlib.pyplot as plt
 import numpy as np
-import cartopy.crs as ccrs
 
 # User-defined libraries
-from read_netcdf import read_n_2d_fields
+from esmvaltool.diag_scripts.ensclus.read_netcdf import read_n_2d_fields
 
 
 def ens_plots(dir_output, dir_plot, name_outputs, numclus,
diff --git a/esmvaltool/diag_scripts/ensclus/ensclus.py b/esmvaltool/diag_scripts/ensclus/ensclus.py
index a669a9a02d..df06ea411c 100644
--- a/esmvaltool/diag_scripts/ensclus/ensclus.py
+++ b/esmvaltool/diag_scripts/ensclus/ensclus.py
@@ -13,16 +13,21 @@
     20170710-mavilia_irene: routines written.
 """
 
-import os
 import logging
+import os
+
 import numpy as np
-from esmvaltool.diag_scripts.shared import group_metadata, run_diagnostic
-from esmvaltool.diag_scripts.shared import ProvenanceLogger, sorted_metadata
 
 # Import user diagnostic routines
-from ens_anom import ens_anom
-from ens_eof_kmeans import ens_eof_kmeans
-from ens_plots import ens_plots
+from esmvaltool.diag_scripts.ensclus.ens_anom import ens_anom
+from esmvaltool.diag_scripts.ensclus.ens_eof_kmeans import ens_eof_kmeans
+from esmvaltool.diag_scripts.ensclus.ens_plots import ens_plots
+from esmvaltool.diag_scripts.shared import (
+    ProvenanceLogger,
+    group_metadata,
+    run_diagnostic,
+    sorted_metadata,
+)
 
 logger = logging.getLogger(os.path.basename(__file__))
 
diff --git a/esmvaltool/diag_scripts/examples/correlate.py b/esmvaltool/diag_scripts/examples/correlate.py
index 171a24a51f..052f3d2bdc 100644
--- a/esmvaltool/diag_scripts/examples/correlate.py
+++ b/esmvaltool/diag_scripts/examples/correlate.py
@@ -6,7 +6,7 @@
 from iris.analysis import MEAN
 from iris.analysis.stats import pearsonr
 
-from diagnostic import plot_diagnostic
+from esmvaltool.diag_scripts.examples.diagnostic import plot_diagnostic
 from esmvaltool.diag_scripts.shared import group_metadata, run_diagnostic
 
 logger = logging.getLogger(os.path.basename(__file__))
diff --git a/esmvaltool/diag_scripts/iht_toa/single_model_diagnostics.py b/esmvaltool/diag_scripts/iht_toa/single_model_diagnostics.py
index fc917a8aa5..e56240c67a 100644
--- a/esmvaltool/diag_scripts/iht_toa/single_model_diagnostics.py
+++ b/esmvaltool/diag_scripts/iht_toa/single_model_diagnostics.py
@@ -16,8 +16,8 @@
 import numpy as np
 from iris import NameConstraint
 from matplotlib import gridspec, rcParams
-from poisson_solver import SphericalPoisson
 
+from esmvaltool.diag_scripts.iht_toa.poisson_solver import SphericalPoisson
 from esmvaltool.diag_scripts.shared import (
     group_metadata,
     run_diagnostic,
diff --git a/esmvaltool/diag_scripts/mpqb/mpqb_lineplot.py b/esmvaltool/diag_scripts/mpqb/mpqb_lineplot.py
index fe4cf2deb0..b08a2cd012 100644
--- a/esmvaltool/diag_scripts/mpqb/mpqb_lineplot.py
+++ b/esmvaltool/diag_scripts/mpqb/mpqb_lineplot.py
@@ -7,9 +7,9 @@
 import iris
 import matplotlib.dates as mdates
 import matplotlib.pyplot as plt
-from mpqb_utils import get_mpqb_cfg
 
 import esmvaltool.diag_scripts.shared.iris_helpers as ih
+from esmvaltool.diag_scripts.mpqb.mpqb_utils import get_mpqb_cfg
 from esmvaltool.diag_scripts.shared import group_metadata, run_diagnostic
 from esmvaltool.diag_scripts.shared._base import (
     ProvenanceLogger,
diff --git a/esmvaltool/diag_scripts/mpqb/mpqb_lineplot_anncyc.py b/esmvaltool/diag_scripts/mpqb/mpqb_lineplot_anncyc.py
index e2a5a662b2..d087499d60 100644
--- a/esmvaltool/diag_scripts/mpqb/mpqb_lineplot_anncyc.py
+++ b/esmvaltool/diag_scripts/mpqb/mpqb_lineplot_anncyc.py
@@ -6,8 +6,8 @@
 
 import iris
 import matplotlib.pyplot as plt
-from mpqb_utils import get_mpqb_cfg
 
+from esmvaltool.diag_scripts.mpqb.mpqb_utils import get_mpqb_cfg
 from esmvaltool.diag_scripts.shared import group_metadata, run_diagnostic
 from esmvaltool.diag_scripts.shared._base import (
     ProvenanceLogger,
diff --git a/esmvaltool/diag_scripts/mpqb/mpqb_lineplot_growthrate.py b/esmvaltool/diag_scripts/mpqb/mpqb_lineplot_growthrate.py
index a6ec136966..90662be72b 100644
--- a/esmvaltool/diag_scripts/mpqb/mpqb_lineplot_growthrate.py
+++ b/esmvaltool/diag_scripts/mpqb/mpqb_lineplot_growthrate.py
@@ -8,9 +8,9 @@
 import matplotlib.dates as mdates
 import matplotlib.pyplot as plt
 import numpy as np
-from mpqb_utils import get_mpqb_cfg
 
 import esmvaltool.diag_scripts.shared.iris_helpers as ih
+from esmvaltool.diag_scripts.mpqb.mpqb_utils import get_mpqb_cfg
 from esmvaltool.diag_scripts.shared import group_metadata, run_diagnostic
 from esmvaltool.diag_scripts.shared._base import (
     ProvenanceLogger,
diff --git a/esmvaltool/diag_scripts/weighting/climwip/calibrate_sigmas.py b/esmvaltool/diag_scripts/weighting/climwip/calibrate_sigmas.py
index 0998346b22..1d85e1f95a 100644
--- a/esmvaltool/diag_scripts/weighting/climwip/calibrate_sigmas.py
+++ b/esmvaltool/diag_scripts/weighting/climwip/calibrate_sigmas.py
@@ -6,7 +6,13 @@
 import matplotlib.pyplot as plt
 import numpy as np
 import xarray as xr
-from core_functions import (
+from scipy.optimize import brute
+
+from esmvaltool.diag_scripts.shared import (
+    get_diagnostic_filename,
+    get_plot_filename,
+)
+from esmvaltool.diag_scripts.weighting.climwip.core_functions import (
     area_weighted_mean,
     calculate_model_distances,
     calculate_weights,
@@ -14,17 +20,11 @@
     compute_overall_mean,
     weighted_quantile,
 )
-from io_functions import (
+from esmvaltool.diag_scripts.weighting.climwip.io_functions import (
     read_metadata,
     read_model_data,
     read_model_data_ancestor,
 )
-from scipy.optimize import brute
-
-from esmvaltool.diag_scripts.shared import (
-    get_diagnostic_filename,
-    get_plot_filename,
-)
 
 logger = logging.getLogger(os.path.basename(__file__))
 
diff --git a/esmvaltool/diag_scripts/weighting/climwip/main.py b/esmvaltool/diag_scripts/weighting/climwip/main.py
index 6b22399547..f46db62090 100644
--- a/esmvaltool/diag_scripts/weighting/climwip/main.py
+++ b/esmvaltool/diag_scripts/weighting/climwip/main.py
@@ -10,15 +10,22 @@
 import numpy as np
 import seaborn as sns
 import xarray as xr
-from calibrate_sigmas import calibrate_performance_sigma
-from core_functions import (
+
+from esmvaltool.diag_scripts.shared import (
+    get_diagnostic_filename,
+    get_plot_filename,
+    run_diagnostic,
+)
+from esmvaltool.diag_scripts.weighting.climwip.calibrate_sigmas import (
+    calibrate_performance_sigma, )
+from esmvaltool.diag_scripts.weighting.climwip.core_functions import (
     area_weighted_mean,
     calculate_model_distances,
     calculate_weights,
     combine_ensemble_members,
     compute_overall_mean,
 )
-from io_functions import (
+from esmvaltool.diag_scripts.weighting.climwip.io_functions import (
     log_provenance,
     read_metadata,
     read_model_data,
@@ -27,11 +34,6 @@
     read_observation_data_ancestor,
 )
 
-from esmvaltool.diag_scripts.shared import (
-    get_diagnostic_filename,
-    get_plot_filename,
-    run_diagnostic,
-)
 
 logger = logging.getLogger(os.path.basename(__file__))
 
diff --git a/esmvaltool/diag_scripts/weighting/plot_utilities.py b/esmvaltool/diag_scripts/weighting/plot_utilities.py
index 706dfd64c9..a071a3c7f0 100644
--- a/esmvaltool/diag_scripts/weighting/plot_utilities.py
+++ b/esmvaltool/diag_scripts/weighting/plot_utilities.py
@@ -3,7 +3,8 @@
 
 import xarray as xr
 
-from climwip.core_functions import weighted_quantile
+from esmvaltool.diag_scripts.weighting.climwip.core_functions import (
+    weighted_quantile, )
 
 
 def read_weights(filename: str) -> dict:
diff --git a/esmvaltool/diag_scripts/weighting/weighted_temperature_graph.py b/esmvaltool/diag_scripts/weighting/weighted_temperature_graph.py
index 5b871283d7..f5f1086e1e 100644
--- a/esmvaltool/diag_scripts/weighting/weighted_temperature_graph.py
+++ b/esmvaltool/diag_scripts/weighting/weighted_temperature_graph.py
@@ -10,13 +10,16 @@
 import matplotlib.pyplot as plt
 import numpy as np
 import xarray as xr
-from climwip.io_functions import log_provenance, read_model_data
 
 from esmvaltool.diag_scripts.shared import (
     get_diagnostic_filename,
     get_plot_filename,
     run_diagnostic,
 )
+from esmvaltool.diag_scripts.weighting.climwip.io_functions import (
+    log_provenance,
+    read_model_data,
+)
 from esmvaltool.diag_scripts.weighting.plot_utilities import (
     calculate_percentiles,
     read_metadata,
diff --git a/esmvaltool/diag_scripts/weighting/weighted_temperature_map.py b/esmvaltool/diag_scripts/weighting/weighted_temperature_map.py
index 7a161d2535..667a382d94 100644
--- a/esmvaltool/diag_scripts/weighting/weighted_temperature_map.py
+++ b/esmvaltool/diag_scripts/weighting/weighted_temperature_map.py
@@ -11,13 +11,16 @@
 import matplotlib.pyplot as plt
 import numpy as np
 from cartopy.mpl.ticker import LatitudeFormatter, LongitudeFormatter
-from climwip.io_functions import log_provenance, read_model_data
 
 from esmvaltool.diag_scripts.shared import (
     get_diagnostic_filename,
     get_plot_filename,
     run_diagnostic,
 )
+from esmvaltool.diag_scripts.weighting.climwip.io_functions import (
+    log_provenance,
+    read_model_data,
+)
 from esmvaltool.diag_scripts.weighting.plot_utilities import (
     calculate_percentiles,
     read_metadata,
diff --git a/esmvaltool/diag_scripts/zmnam/zmnam.py b/esmvaltool/diag_scripts/zmnam/zmnam.py
index d86ee48458..c0450c117a 100644
--- a/esmvaltool/diag_scripts/zmnam/zmnam.py
+++ b/esmvaltool/diag_scripts/zmnam/zmnam.py
@@ -20,9 +20,12 @@
 
 from esmvaltool.diag_scripts.shared import ProvenanceLogger, run_diagnostic
 # Import zmnam diagnostic routines
-from zmnam_calc import zmnam_calc
-from zmnam_plot import zmnam_plot
-from zmnam_preproc import (zmnam_preproc, zmnam_preproc_clean)
+from esmvaltool.diag_scripts.zmnam.zmnam_calc import zmnam_calc
+from esmvaltool.diag_scripts.zmnam.zmnam_plot import zmnam_plot
+from esmvaltool.diag_scripts.zmnam.zmnam_preproc import (
+    zmnam_preproc,
+    zmnam_preproc_clean,
+)
 
 logger = logging.getLogger(__name__)
 
diff --git a/setup.cfg b/setup.cfg
index c738c5d716..e28f8079a0 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -1,5 +1,6 @@
 [tool:pytest]
 addopts =
+    --import-mode=importlib
     --doctest-modules
     --ignore=doc/sphinx/source/conf.py
     --cov=esmvaltool