diff --git a/2dtestcases.py b/2dtestcases.py
index 334aa9a48465e0d4281f565ccdc0b0a93a6cc22d..70bf2a5e209e2aaa4c54efcefb57099b6305a5b5 100755
--- a/2dtestcases.py
+++ b/2dtestcases.py
@@ -486,7 +486,13 @@ chemIntegrator = None
 if chemtype == "ozoneco":
     l = l + mechbib.leightonco_reacs()
     mech = mechanism(l)
-    chemIntegrator = "EulerBackward"
+    mech.define_integrator('EulerBackward',
+                           errortarget_relative=1.0e-6,
+                           nitersmax=10000
+                           )
+else:
+    mech = mechanism(l)
+    mech.define_integrator('null')
 #############################################################
 
 
@@ -554,13 +560,10 @@ myctm = ctm(
     vts=vts,
     dt=dt,
     MassFluxFunc=MassFluxFunc,
-    chemIntegrator=chemIntegrator,
     split=split,
     label=label,
     prescset=defaults.prescset,
     ndtchem=ndtchem,
-    ebitolerance=1.0e-6,
-    ebinitersmax=10000,
 )
 #############################################################
 
diff --git a/CItest-MRconservation.py b/CItest-MRconservation.py
index e6a76d7ad98b5d12e469a5d0dc5ac8ac5ab15c53..5651f9058b33566075cdd3ec7c77626933c7dc75 100755
--- a/CItest-MRconservation.py
+++ b/CItest-MRconservation.py
@@ -109,6 +109,7 @@ mrdiag = True
 
 l = TRC_reacs()
 mech = mechanism(l)
+mech.define_integrator('null')
 
 # Grid creation
 maxlat = 45.0
@@ -163,8 +164,7 @@ for tsname in tslist:
             vts=defaults.vts,
             prescset=defaults.prescset,
             dt=dt,
-            MassFluxFunc=MassFluxFunc,
-            chemIntegrator=None,
+            MassFluxFunc=MassFluxFunc
         )
         masstest, mrtest = test(myctm)
         if masstest and mrtest:
diff --git a/CItest-PrescEmi.py b/CItest-PrescEmi.py
index aea9efd26353afb1509601d390f783b686dd03b0..21cb9756b8e4ddf5c24cb36e70ce62b209f3cd6e 100755
--- a/CItest-PrescEmi.py
+++ b/CItest-PrescEmi.py
@@ -144,6 +144,7 @@ airdensity = 2.69e25
 
 l = trop_ozone_reacs() + CO_hotwo_reacs() + TRC_reacs()
 mech = mechanism(l)
+mech.define_integrator('EulerBackward')
 
 # Initial concentrations
 
@@ -204,8 +205,7 @@ myctm = ctm(
     ts=ts,
     vts=defaults.vts,
     dt=1200.0,
-    MassFluxFunc=MassFluxFunc,
-    chemIntegrator="EulerBackward",
+    MassFluxFunc=MassFluxFunc
 )
 
 # Some screen outputs to see if things go fine
diff --git a/CItest-chem.py b/CItest-chem.py
index 340f138edee2befc2891d92d49167b6deb215bad..8aba0f937bbaed01823df173b41810bf601d9ddb 100755
--- a/CItest-chem.py
+++ b/CItest-chem.py
@@ -75,6 +75,7 @@ ts_b_dict = {}
 ts_c_dict = {}
 
 for exp in experimentlist:
+    mech.define_integrator(exp[0])
     myctm = ctm(
         mech=mech,
         finaltime=finaltime,
@@ -85,7 +86,6 @@ for exp in experimentlist:
         tdgrid=defaults.tdgrid,
         MassFluxFunc=defaults.MassFluxFunc,
         dt=exp[1],
-        chemIntegrator=exp[0],
     )
     ts_a = ((myctm.t, myctm.activeConc_ppb("A")[0, 0, 0]),)
     ts_b = ((myctm.t, myctm.activeConc_ppb("B")[0, 0, 0]),)
diff --git a/CItest-compressible.py b/CItest-compressible.py
index d5fdcb9aab79998afe092862342eb40d654d8ae2..51597d75faf6ca19a3c5bf0f2569d044d69b6000 100755
--- a/CItest-compressible.py
+++ b/CItest-compressible.py
@@ -279,6 +279,8 @@ tdgrid = Threedimgrid("3d grid", myGrid, dzvector, vts.nneigh, vertperiod=False)
 # Define mechanism
 l = TRC_reacs()
 mech = mechanism(l)
+mech.define_integrator('null')
+
 #############################################################
 
 #############################################################
@@ -311,7 +313,6 @@ myctm = ctm(
     vts=vts,
     dt=dt,
     MassFluxFunc=MassFluxFunc,
-    chemIntegrator=None,
     prescset=defaults.prescset,
 )
 #############################################################
diff --git a/CItest-isopleth.py b/CItest-isopleth.py
index ddb119f37333d05c086d68240840ee01ab1f99f1..0d3dae48b34a27c34f3e88f991c8ab804ce2dda9 100755
--- a/CItest-isopleth.py
+++ b/CItest-isopleth.py
@@ -118,6 +118,7 @@ def add_standard_emi(dict, **kwargs):
 
 l = trop_ozone_reacs() + hotwo_reacs() + NMVOC_reacs()
 mech = mechanism(l)
+mech.define_integrator('EulerBackward')
 
 tdgrid = defaults.tdgrid
 MassFluxFunc = defaults.MassFluxFunc
@@ -164,7 +165,6 @@ for nox in noxvec:
             ts=defaults.ts,
             vts=defaults.vts,
             dt=100.0,
-            chemIntegrator="EulerBackward",
             MassFluxFunc=MassFluxFunc,
         )
 
diff --git a/CItest-leighton-radicals.py b/CItest-leighton-radicals.py
index efe19bde792d819ab9d0350febc7047ac6535a8c..fee9b89cbdf9e99012b2314b1252eeaa9b83b728 100755
--- a/CItest-leighton-radicals.py
+++ b/CItest-leighton-radicals.py
@@ -48,6 +48,9 @@ Mh2o = 0.018
 # Chemical scheme creation
 
 mech = mechanism(mechbib.leightonco_reacs())
+mech.define_integrator('EulerBackward',
+                       nitersmax = 10000,
+                       errortarget_relative=1.0e-6)
 
 # Initial concentrations
 inidict = {}
@@ -95,8 +98,6 @@ myctm = ctm(
     ts=defaults.ts,
     vts=defaults.vts,
     prescset=defaults.prescset,
-    ebitolerance=1.0e-6,
-    ebinitersmax=10000,
 )
 o3 = myctm.centerConc("O3")
 no = myctm.centerConc("NO")
diff --git a/CItest-xz.py b/CItest-xz.py
index 687ccb97165eb6103ff62144395ec26ca33b69e5..10ad8c202c063270b351fe19fbdec6e6f0055b47 100755
--- a/CItest-xz.py
+++ b/CItest-xz.py
@@ -473,6 +473,7 @@ print("dt", dt)
 
 l = tracer_reacs()
 mech = mechanism(l)  # These two lines create fake chemistry with one inert tracer
+mech.define_integrator('null')
 
 # No emissions
 emidict = {}
@@ -642,8 +643,7 @@ myctm = ctm(
     prescset=defaults.prescset,
     inidict=inidict,
     lagdict=lagdict,
-    label=label,
-    chemIntegrator=None,
+    label=label
 )
 
 if myctm.vts.name == "ppmSL":
diff --git a/toyctm/toyctm.py b/toyctm/toyctm.py
index 451c26c6277ae40deb31573c73c09dfe1c713f45..40f0a509ab8351b6fc781409172f67783ea37c20 100644
--- a/toyctm/toyctm.py
+++ b/toyctm/toyctm.py
@@ -62,8 +62,6 @@ class ctm:
     - MassFluxFunc   : Instance of the MassFluxGenerator class (examples in threedmassfluxes.py)
     - lagdict        : defaults to {}
     - label          : string tha gives a name to the simulation
-    - chemIntegrator : time scheme for chemistry. So far available,
-                       EulerBackward (strongly recommended) and EulerForward
     """
 
         # read Mandatory arguments
@@ -94,6 +92,14 @@ class ctm:
                 "toyctm.py ERROR: Please provide a chemical mechanism with mech= when calling ctm()"
             )
             sys.exit(1)
+        if self.mech.integrator_name is None:
+            print(
+                """toyctm.py ERROR: Please initialize the chemical integrator 
+                   with mech.define_integrator(). If no chemistry integration
+                   is needed please use mech.define_integrator('null')"""
+            )
+            sys.exit(1)
+            
         if self.tdgrid is None:
             print(
                 "toyctm.py ERROR: Please provide a grid with tdgrid= when calling ctm()"
@@ -133,10 +139,18 @@ class ctm:
             sys.exit(1)
 
         # read optional arguments
-        self.chemIntegrator = kwargs.get(
-            "chemIntegrator", "EulerBackward"
-        )  # Chemical integrator. 'EulerBackward', 'EulerForward' or None are supported.
-           # 'EulerBackward' is strongly recommended if you have chemistry.
+        
+        # Deprecated arguments produce an error
+        if 'chemIntegrator' in kwargs :
+            print('ERROR: chemIntegrator is deprecated, please use mech.define_integrator')
+            sys.exit(1)
+        if 'ebitolerance' in kwargs :
+            print('ERROR: ebitolerance is deprecated, please use mech.define_integrator')
+            sys.exit(1)
+        if 'ebinitersmax' in kwargs :
+            print('ERROR: ebinitersmax is deprecated, please use mech.define_integrator')
+            sys.exit(1)
+        
         inidict = kwargs.get(
             "inidict", {}
         )  # Dictionnary associating an instance of the initialization class to each model species.
@@ -166,8 +180,6 @@ class ctm:
             "split", "lie"
         )  # splitting strategy, so far 'strang' or 'lie'
         self.ndtchem = kwargs.get("ndtchem", None)
-        self.ebitolerance = kwargs.get("ebitolerance", 1.0e-2)
-        self.ebinitersmax = kwargs.get("ebinitersmax", 1000)
 
         # End of arguments read
 
@@ -216,11 +228,6 @@ class ctm:
         self.zdconc = np.zeros(self.conc.shape)
         self.mdconc = np.zeros(self.conc.shape)
 
-        # Initialize the chemistry solver
-        if self.chemIntegrator is not None : 
-            self.mech.define_integrator(self.chemIntegrator,
-                                        errortarget_relative = self.ebitolerance,
-                                        nitersmax            = self.ebinitersmax)
         # Misc initializations
         self.niters = -999
         self.vdconc = 0.
@@ -693,7 +700,7 @@ class ctm:
             self.emi()
 
         # dtchem for chemistry
-        if self.chemIntegrator is not None:
+        if self.mech.integrator is not None:
             # CHEMISTRY
             self.chemstep()
             for sp in self.species:
@@ -706,7 +713,7 @@ class ctm:
     def chemupdate_reverse(self):
         """ Everything except transport on time dtchem in reverse order
         during time dtchem"""
-        if self.chemIntegrator is not None:
+        if self.mech.integrator is not None:
             # CHEMISTRY
             self.chemstep()
             for sp in self.species:
diff --git a/toyctm_addons/chemistry.py b/toyctm_addons/chemistry.py
index 51459ae1b24ecfc7446ce051df8dcf443184f387..04f108f9e8c15311564159e50bc865a569b307b7 100644
--- a/toyctm_addons/chemistry.py
+++ b/toyctm_addons/chemistry.py
@@ -309,6 +309,8 @@ class mechanism:
             self.errortarget_absolute = self.minimalconcentration
         elif integrator_name == 'EulerForward' : 
             self.integrator = self.eulerForwardTimestep
+        elif integrator_name == 'null' :
+            self.integrator = None
         else : 
             print('ERROR : integrator', integrator_name, 'not available')
             sys.exit(1)