diff --git a/opendrift/models/basemodel/environment.py b/opendrift/models/basemodel/environment.py
index fcd186bdd..d5844a587 100644
--- a/opendrift/models/basemodel/environment.py
+++ b/opendrift/models/basemodel/environment.py
@@ -639,10 +639,10 @@ def get_environment(self, variables, time, lon, lat, z, profiles=None, profiles_
                         profiles_from_reader = None
                     env_tmp, env_profiles_tmp = \
                         reader.get_variables_interpolated(
-                            variable_group, profiles_from_reader,
-                            profiles_depth, time,
-                            lon[missing_indices], lat[missing_indices],
-                            z[missing_indices], self.proj_latlon)
+                            variable_group, profiles = profiles_from_reader,
+                            profiles_depth = profiles_depth, time = time,
+                            lon=lon[missing_indices], lat=lat[missing_indices],
+                            z=z[missing_indices], rotate_to_proj=self.proj_latlon)
 
                 except NotCoveredError as e:
                     logger.info(e)
diff --git a/opendrift/readers/operators/ops.py b/opendrift/readers/operators/ops.py
index 4e673a1c3..a3e83f24d 100644
--- a/opendrift/readers/operators/ops.py
+++ b/opendrift/readers/operators/ops.py
@@ -5,7 +5,9 @@
 import xarray as xr
 import numpy as np
 class Combine:
-
+    """Combining two readers into a third one. You can use usual operators,
+    but also more complex ones such as gaussian combining.
+    """
 
     def __add__(self, other):
         from .readerops import Combined as ReaderCombined
@@ -38,6 +40,10 @@ def __sub__(self, other):
         return self + (-1 * other)
 
     def combine_gaussian(self, measurement_reader, std):
+        """Mix two readers with a gaussian, whose std is the one given as an argument.
+        The measurment reader have to be of type timeseries, with a lon and lat
+        attributes that are taken as the center of the measure.
+        """
         from .readerops import Combined as ReaderCombined
 
         def gaussian_melting(a, b, lon, lat, lon_center, lat_center, std):
@@ -46,21 +52,9 @@ def gaussian_melting(a, b, lon, lat, lon_center, lat_center, std):
                 lon_center = lon_center * np.ones(len(lon))
             if ( isinstance(lat_center, float) or len(lat_center) == 1 ) and len(lat) != 1 :
                 lat_center = lat_center * np.ones(len(lat))
-
-            print("lon length: {}".format(len(lon)))
-            print("lat length: {}".format(len(lat)))
-            print("lon_center length: {}".format(len(lon_center)))
-            print("lat_center length: {}".format(len(lat_center)))
-
             _, _, distance = geod.inv(lon, lat, lon_center, lat_center)
             exponential_factor = np.exp( -np.power(distance/std, 2.) / 2)
             res = a * (1 - exponential_factor) + b * exponential_factor
-            print("==================================")
-            print("Gaussian calculation")
-            print("Lon, lat, lon_center, lat_center = {}".format((lon, lat, lon_center, lat_center)))
-            print("Distance used for calculation: {}".format(distance))
-            print("a = {}, b = {}, c = {}".format(a, b, res))
-            print("=================================")
             return res
 
         return ReaderCombined(a = self, b = measurement_reader, op = gaussian_melting, op_type= "combine_gaussian", external_params = std)
diff --git a/opendrift/readers/operators/readerops.py b/opendrift/readers/operators/readerops.py
index 7618db506..1c276ab66 100644
--- a/opendrift/readers/operators/readerops.py
+++ b/opendrift/readers/operators/readerops.py
@@ -20,7 +20,7 @@ class Combined(BaseReader):
     def __init__(self, a, b, op, op_type = "easy", external_params = None):
         '''Combine two readers a and b followinf the operator op. If needed,
         you can ad an op_type that will enable you to use the external parameters
-        you want in your op. '''
+        you want in your op between lines 63 and 73.'''
         self.a = a
         self.b = b
         self.op = op
@@ -49,16 +49,11 @@ def covers_positions(self, lon, lat):
     def covers_time(self, time):
         return self.a.covers_time(time) and self.b.covers_time(time)
 
-    def get_variables_interpolated(self, variables, *args, **kwargs):
+    def get_variables_interpolated(self, variables, **kwargs):
         assert set(variables).issubset(self.variables), f"{variables} is not subset of {self.variables}"
 
-        debug=True
-        if debug:
-            print("args passed are: {}".format(args))
-            print("kwargs passed are: {}".format(kwargs))
-
-        env_a, env_profiles_a = self.a.get_variables_interpolated(variables, *args, **kwargs)
-        env_b, env_profiles_b = self.b.get_variables_interpolated(variables, *args, **kwargs)
+        env_a, env_profiles_a = self.a.get_variables_interpolated(variables, **kwargs)
+        env_b, env_profiles_b = self.b.get_variables_interpolated(variables, **kwargs)
 
         variables = [
             var for var in env_a.keys() if var not in ['x', 'y', 'z', 'time']