diff --git a/pyaerocom/helpers.py b/pyaerocom/helpers.py
index 326b86990..e04fd0b8c 100644
--- a/pyaerocom/helpers.py
+++ b/pyaerocom/helpers.py
@@ -6,7 +6,7 @@
 from cf_units import Unit
 from collections import Counter
 from datetime import MINYEAR, datetime, date
-import iris
+import iris, iris.analysis, iris.cube, iris.coords
 import math as ma
 import numpy as np
 import pandas as pd
diff --git a/pyaerocom/io/iris_io.py b/pyaerocom/io/iris_io.py
index b32df7870..b6a2ddf0d 100644
--- a/pyaerocom/io/iris_io.py
+++ b/pyaerocom/io/iris_io.py
@@ -9,7 +9,7 @@
 
 import cf_units
 from datetime import datetime
-import iris
+import iris, iris.coords, iris.util
 try:
     # as of iris version 3
     from iris.util  import equalise_attributes
@@ -546,11 +546,7 @@ def _check_correct_dtypes_timedim_cube_list(cubes):
 def concatenate_iris_cubes(cubes, error_on_mismatch=True):
     """Concatenate list of :class:`iris.Cube` instances cubes into single Cube
 
-    Helper method for concatenating list of cubes and that helps
-    with handling the fact that the corresponding iris method is not well
-    defined in the sense of what it returns (i.e. instance of
-    :class:`Cube` or :class:`CubeList`, depending on whether all cubes
-    could be concatenated or not...)
+    Helper method for concatenating list of cubes
 
     This method is not supposed to be called directly but rather
     :func:`concatenate_cubes` (which ALWAYS returns instance of
@@ -559,7 +555,7 @@ def concatenate_iris_cubes(cubes, error_on_mismatch=True):
 
     Parameters
     ----------
-    cubes : CubeList
+    cubes : CubeList or list(Cubes)
         list of individual cubes
     error_on_mismatch
         boolean specifying whether an Exception is supposed to be raised
@@ -567,7 +563,7 @@ def concatenate_iris_cubes(cubes, error_on_mismatch=True):
 
     Returns
     -------
-    :obj:`Cube` or :obj:`CubeList`
+    :obj:`Cube`
         result of concatenation
 
     Raises
@@ -577,6 +573,7 @@ def concatenate_iris_cubes(cubes, error_on_mismatch=True):
         into a single instance of :class:`iris.Cube` class.
 
     """
+    cubes = iris.cube.CubeList(cubes)
     var_name = cubes[0].var_name
     if const.GRID_IO.EQUALISE_METADATA:
         meta_init = cubes[0].metadata
@@ -595,14 +592,15 @@ def concatenate_iris_cubes(cubes, error_on_mismatch=True):
     #unify time units
     iris.util.unify_time_units(cubes)
     #now concatenate the cube list to one cube
+
     try:
-        cubes_concat = iris._concatenate.concatenate(cubes, error_on_mismatch)
-    except Exception as e: #
+        cubes_concat = iris.cube.CubeList.concatenate_cube(cubes, error_on_mismatch)
+    except Exception as e:
         if _check_correct_dtypes_timedim_cube_list(cubes):
-            cubes_concat = iris._concatenate.concatenate(cubes,
-                                                         error_on_mismatch)
+            cubes_concat = iris.cube.CubeList.concatenate_cube(cubes,
+                                                               error_on_mismatch)
         else:
             raise
 
-    return cubes_concat[0]
+    return cubes_concat
 
diff --git a/pyaerocom_env.yml b/pyaerocom_env.yml
index 103619fca..6bafb632d 100644
--- a/pyaerocom_env.yml
+++ b/pyaerocom_env.yml
@@ -3,7 +3,7 @@
 channels:
   - conda-forge
 dependencies:
-  - iris >=3.0.1,<3.1.0
+  - iris >=3.1.0
   - xarray >=0.16.0
   - cartopy >=0.16.0
   - matplotlib >=3.0.1
diff --git a/setup.py b/setup.py
index 142c4c02f..7720ca135 100644
--- a/setup.py
+++ b/setup.py
@@ -48,8 +48,7 @@
 
     python_requires = '>=3.7,<4',
     install_requires = [
-        'scitools-iris>=3.0.1,<3.1.0',
-        'scitools-pyke>=1.1.1',
+        'scitools-iris>=3.1.0',
         'xarray>=0.16.0',
         'cartopy>=0.16.0',
         'matplotlib>=3.0.1',