diff --git a/code/m_chemistry.f90 b/code/m_chemistry.f90
index 28d11a2..4dd3dce 100644
--- a/code/m_chemistry.f90
+++ b/code/m_chemistry.f90
@@ -10,34 +10,82 @@ real function rate_1step (yr, theta)
 
     real, intent(in)  :: yr
     real, intent(in)  :: theta
-    
+
     real :: y
     real :: t_reduce
-    
+
     y=yr
-    if(yr.lt.0.) y=0.
-    if(yr.gt.1.) y=1.
-    
+    ! if(yr.lt.0.) y=0.
+    ! if(yr.gt.1.) y=1.
+
     t_reduce=theta
-    if(theta.lt.0.) t_reduce=0.
-    if(theta.gt.1.) t_reduce=1.
-    
+    ! if(theta.lt.0.) t_reduce=0.
+    ! if(theta.gt.1.) t_reduce=1.
+
     if (t_reduce.gt.c_ref) then
-    
+
       rate_1step = pre*y*exp(-ac/(1.+bc*t_reduce))
-    
+
     else if (t_reduce.le.c_cut) then
-    
+
       rate_1step = min_wr
-    
+
     else
-    
+
       rate_1step = &
         ((refwr-min_wr)*exp(prof_wr*(t_reduce-c_ref)) + min_wr - refwr*exp(prof_wr*(c_cut-c_ref))) &
         / (1.-exp(prof_wr*(c_cut-c_ref)))
-    
+
    endif
 
 end function rate_1step
 
+
+real function rate1_2step (ya, yx, theta)
+
+    real, intent(in)  :: ya
+    real, intent(in)  :: yx
+    real, intent(in)  :: theta
+
+    real :: y1
+    real :: y2
+    real :: t_reduce
+
+    y1=ya
+    if(ya.lt.0.) y1=0.
+    if(ya.gt.1.) y1=1.
+
+    y2=yx
+    if(yx.lt.0.) y2=0.
+    if(yx.gt.1.) y2=1.
+
+    t_reduce=theta
+    if(theta.lt.0.) t_reduce=0.
+    if(theta.gt.1.) t_reduce=1.
+
+    rate1_2step = 0.
+
+end function rate1_2step
+
+
+real function rate2_2step (yx, theta)
+
+    real, intent(in)  :: yx
+    real, intent(in)  :: theta
+
+    real :: y
+    real :: t_reduce
+
+    y=yx
+    if(yx.lt.0.) y=0.
+    if(yx.gt.1.) y=1.
+
+    t_reduce=theta
+    if(theta.lt.0.) t_reduce=0.
+    if(theta.gt.1.) t_reduce=1.
+
+    rate2_2step = 0.
+
+end function rate2_2step
+
 end module m_chemistry