diff --git a/m_fdm_calc.f90 b/m_fdm_calc.f90
index 6281fd8..5909692 100644
--- a/m_fdm_calc.f90
+++ b/m_fdm_calc.f90
@@ -30,6 +30,8 @@ module m_fdm_calc
 
       integer, parameter :: neq = 2
 
+      logical :: fors
+
 !===========================================================================
 !===========================================================================
 
@@ -217,6 +219,8 @@ module m_fdm_calc
          return
       endif
 
+      fors = .true.
+
       allocate(u_(nx,ny,nz))
       allocate(v_(nx,ny,nz))
       allocate(w_(nx,ny,nz))
@@ -655,6 +659,80 @@ module m_fdm_calc
       return
       END SUBROUTINE fns
 
+
+      subroutine EE1(xx,yy,zz,uu_,vv_,ww_,yy1,rhs1)
+
+      implicit none
+
+      integer :: xx,yy,zz
+      real*8 :: uu_(xx,yy,zz),vv_(xx,yy,zz),ww_(xx,yy,zz)
+      real*8 :: yy1(xx,yy,zz,neq),rhs1(xx,yy,zz,neq)
+
+      integer :: i, j, k
+
+      CALL fns(yy1,rhs1,xx,yy,zz,uu_,vv_,ww_)
+
+      yy1 = yy1 + fdmdt * rhs1
+
+!==========rho=1 treatment
+      yy1(:,:,:,2) = yy1(:,:,:,2)/yy1(:,:,:,1)
+      yy1(:,:,:,1) = 1.
+
+      DO k = 1,zz
+       DO j = 1,yy
+        DO i = 1,xx
+
+!==========Max Yr=1 treatment
+         yy1(i,j,k,2)=MIN(in_yr,yy1(i,j,k,2))
+
+!==========Min Yr=0 treatment
+!          yy1(i,j,k,2)=MAX(out_yr,yy1(i,j,k,2))
+
+        ENDDO
+       ENDDO
+      ENDDO
+
+      return
+      END SUBROUTINE EE1
+
+
+      subroutine AB2(xx,yy,zz,uu_,vv_,ww_,yy1,rhs0,rhs1)
+
+      implicit none
+
+      integer :: xx,yy,zz
+      real*8 :: uu_(xx,yy,zz),vv_(xx,yy,zz),ww_(xx,yy,zz)
+      real*8 :: yy1(xx,yy,zz,neq),rhs0(xx,yy,zz,neq),rhs1(xx,yy,zz,neq)
+
+      integer :: i, j, k
+
+      CALL fns(yy1,rhs1,xx,yy,zz,uu_,vv_,ww_)
+
+      yy1 = yy1 + fdmdt * ( 1.5d0 * rhs1 - 0.5d0 * rhs0 )
+      rhs0 = rhs1
+
+!==========rho=1 treatment
+      yy1(:,:,:,2) = yy1(:,:,:,2)/yy1(:,:,:,1)
+      yy1(:,:,:,1) = 1.
+
+      DO k = 1,zz
+       DO j = 1,yy
+        DO i = 1,xx
+
+!==========Max Yr=1 treatment
+         yy1(i,j,k,2)=MIN(in_yr,yy1(i,j,k,2))
+
+!==========Min Yr=0 treatment
+!          yy1(i,j,k,2)=MAX(out_yr,yy1(i,j,k,2))
+
+        ENDDO
+       ENDDO
+      ENDDO
+
+      return
+      END SUBROUTINE AB2
+
+
       subroutine RK4(xx,yy,zz,uu_,vv_,ww_,yy1,yy2,yyf)
 
       implicit none
@@ -681,8 +759,14 @@ module m_fdm_calc
        real*8 :: uu_(xx,yy,zz),vv_(xx,yy,zz),ww_(xx,yy,zz)
        real*8 :: yy1(xx,yy,zz,neq),yy2(xx,yy,zz,neq),yyf(xx,yy,zz,neq)
 
-
-       CALL RK4(xx,yy,zz,uu_,vv_,ww_,yy1,yy2,yyf)
+       ! advance reacting scalars - either Euler or Adams-Bashforth
+       if (fors) then
+          call EE1(xx,yy,zz,uu_,vv_,ww_,yy1,yy2)
+          fors = .false.
+       else
+          call AB2(xx,yy,zz,uu_,vv_,ww_,yy1,yy2,yyf)
+          yy2 = yyf
+       end if
 
        return
       END SUBROUTINE solve