diff --git a/code/ysolve.f90 b/code/ysolve.f90
index 6e88986..c2dba97 100644
--- a/code/ysolve.f90
+++ b/code/ysolve.f90
@@ -380,13 +380,7 @@
         CALL d2fnonp(nx,hx,ux(:,:),d2ux(:,:),2,1)
 
         DO i=1,nx
-         wrate=pre*ux(1,i)*exp(-ac/(1.+bc*(ux(2,i))))   !wrate
-
-         IF ((ux(2,i)).le.c_ref) THEN
-           wrate=min_wr
-          IF ((ux(2,i)).gt.c_cut) wrate=((refwr-min_wr)*exp(prof_wr*(ux(2,i)-c_ref))+ &
-          min_wr-refwr*exp(prof_wr*(c_cut-c_ref)))/(1.-exp(prof_wr*(c_cut-c_ref)))
-         ENDIF
+         wrate=rate_1step(ux(1,i), ux(2,i))
 
 ! -0.5*( d(u*Y)/dx + u*dY/dx + Y*du/dx ) + D*d2Y/d2x
          f(i,1) = - ( u(i)*dux(1,i) ) + diff * d2ux(1,i) - wrate
@@ -416,5 +410,41 @@
       f(nx,2)=f(nx,2)-Dt
 
      END SUBROUTINE fns
+
 !------------------------------------------------------------------------
+
+     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.
+
+        t_reduce=theta
+        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
+
     END MODULE ysolve