temporarily fix viscosity ratio to 5

This reverts commit 0d771077f6.

Conflicts:
	code/ysolve.f90

.gitlab-ci.yml

@@ -16,4 +16,6 @@ test:
   script:
     - cd ${CI_PROJECT_DIR}
     - cd code
-    - make test
+    - make clean && make
+    - cd sample/4pi-IC1
+    - ../../ex

code/m_chemistry.f90

@@ -0,0 +1,154 @@
+module m_chemistry
+
+use m_parameters
+
+implicit none
+
+real, private :: coef(10)
+real, private :: lambda_onestep
+real, private :: lambda1_twostep
+real, private :: lambda2_twostep
+real, private :: beta1_twostep
+real, private :: hrp_twostep
+
+contains
+
+subroutine init_chemistry
+
+    character(len=40) :: nrxn_string
+
+    if (nrxn == 1) then
+       reaction_type = "onestep"
+    else if (nrxn == 2) then
+       reaction_type = "twostep"
+    else
+       write(nrxn_string, *) nrxn
+       reaction_type = trim(nrxn_string) // "-step"
+    end if
+
+    if ( reaction_type == "onestep" ) then
+        lambda_onestep = pre * exp ( - beta / hrp )
+    else if ( reaction_type == "twostep" ) then
+        lambda1_twostep = lambda1
+        lambda2_twostep = lambda2
+        beta1_twostep = beta1
+        hrp_twostep = hrp
+    else
+        WRITE(*,*) 'ERROR, UNDEFINED REACTION TYPE ', reaction_type
+        stop
+    end if
+
+end subroutine init_chemistry
+
+
+subroutine update_chemistry (t)
+
+    real :: t
+
+    real :: factor
+
+    real :: relax_duration = 60.
+
+    if ( reaction_type == "onestep" ) then
+        lambda_onestep = pre * exp ( - beta / hrp )
+    else if ( reaction_type == "twostep" ) then
+        if (t < relax_duration) then
+            factor = (relax_duration + t) / relax_duration / 2.
+        else
+            factor = 1.
+        end if
+
+        lambda1_twostep = factor * lambda1
+        lambda2_twostep = factor * lambda2
+
+    else
+        stop
+    end if
+
+end subroutine update_chemistry
+
+
+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
+
+
+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 = lambda1_twostep * y1 * y2 * & 
+                  exp (-(beta1_twostep*(1. - t_reduce))/(1. - hrp_twostep*(1. - t_reduce)))
+
+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 = lambda2_twostep * yx * yx
+
+end function rate2_2step
+
+end module m_chemistry

code/m_compact.f90

@@ -1,14 +1,13 @@
-    MODULE Compact
+    MODULE m_compact
      IMPLICIT NONE
-     PRIVATE
+
-     REAL, DIMENSION(:), ALLOCATABLE :: lxf,lxs,wxf,wxs, & 
+     REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: lxf,lxs,wxf,wxs, &
-                                        lyf,lys,wyf,wys, &
+                                                lyf,lys,wyf,wys, &
-                                        lzf,lzs,wzf,wzs
+                                                lzf,lzs,wzf,wzs
-!                                           lyzf,lyzs,wyzf,wyzs
+
      INTEGER :: nxc,nyc,nzc
-     REAL, PARAMETER :: ezero = 1.0e-14
+
-      
+     REAL(KIND=8), PARAMETER :: ezero = 1.0e-14
-     PUBLIC :: ludcmp,dfnonp,d2fnonp,dfp,d2fp
 
     CONTAINS
 
@@ -20,6 +19,22 @@
        nxc=nx
        nyc=ny
        nzc=nz
+
+       CALL ludcmp_allocate(nx,ny,nz,xp,yp,zp)
+
+       CALL ludcmp_calculate(nx,ny,nz,xp,yp,zp)
+
+     END SUBROUTINE ludcmp
+
+     SUBROUTINE ludcmp_allocate(nx,ny,nz,xp,yp,zp)
+      INTEGER, INTENT(IN) :: nx,ny,nz
+      INTEGER, INTENT(IN) :: xp,yp,zp
+      INTEGER :: ierr
+
+       nxc=nx
+       nyc=ny
+       nzc=nz
+
 !       IF(nyc /= nzc) PRINT*,'ny should be equal nz'
 
 !      xp, yp, zp = 0 : periodic
@@ -32,9 +47,6 @@
          IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
          ALLOCATE(wxs(nxc),STAT=ierr)
          IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
-         CALL p_lud(1,nxc)
-       ELSE
-         CALL nonp_lud(1,nxc)
        ENDIF
 
        ALLOCATE(lyf(nyc),STAT=ierr)
@@ -46,9 +58,6 @@
          IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
          ALLOCATE(wys(nyc),STAT=ierr)
          IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
-         CALL p_lud(2,nyc)
-       ELSE
-         call nonp_lud(2,nyc)
        ENDIF
 
        ALLOCATE(lzf(nzc),STAT=ierr)
@@ -60,19 +69,140 @@
          IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
          ALLOCATE(wzs(nzc),STAT=ierr)
          IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
+       ENDIF
+
+     END SUBROUTINE ludcmp_allocate
+
+     SUBROUTINE ludcmp_deallocate(xp,yp,zp)
+      INTEGER, INTENT(IN) :: xp,yp,zp
+
+!       IF(nyc /= nzc) PRINT*,'ny should be equal nz'
+
+!      xp, yp, zp = 0 : periodic
+       DEALLOCATE(lxf)
+       DEALLOCATE(lxs)
+       IF(xp.eq.0) THEN
+         DEALLOCATE(wxf)
+         DEALLOCATE(wxs)
+       ENDIF
+
+       DEALLOCATE(lyf)
+       DEALLOCATE(lys)
+       IF(yp.eq.0) THEN
+         DEALLOCATE(wyf)
+         DEALLOCATE(wys)
+       ENDIF
+
+       DEALLOCATE(lzf)
+       DEALLOCATE(lzs)
+       IF(zp.eq.0) THEN
+         DEALLOCATE(wzf)
+         DEALLOCATE(wzs)
+       ENDIF
+
+     END SUBROUTINE ludcmp_deallocate
+
+     SUBROUTINE ludcmp_testalloc
+
+       IF (.not. ALLOCATED(lxf)) print *, "lxf not allocated"
+       IF (.not. ALLOCATED(lxs)) print *, "lxs not allocated"
+       IF (.not. ALLOCATED(wxf)) print *, "wxf not allocated"
+       IF (.not. ALLOCATED(wxs)) print *, "wxs not allocated"
+
+       IF (.not. ALLOCATED(lyf)) print *, "lyf not allocated"
+       IF (.not. ALLOCATED(lys)) print *, "lys not allocated"
+       IF (.not. ALLOCATED(wyf)) print *, "wyf not allocated"
+       IF (.not. ALLOCATED(wys)) print *, "wys not allocated"
+
+       IF (.not. ALLOCATED(lzf)) print *, "lzf not allocated"
+       IF (.not. ALLOCATED(lzs)) print *, "lzs not allocated"
+       IF (.not. ALLOCATED(wzf)) print *, "wzf not allocated"
+       IF (.not. ALLOCATED(wzs)) print *, "wzs not allocated"
+
+     END SUBROUTINE ludcmp_testalloc
+
+     SUBROUTINE ludcmp_calculate(nx,ny,nz,xp,yp,zp)
+      INTEGER, INTENT(IN) :: nx,ny,nz
+      INTEGER, INTENT(IN) :: xp,yp,zp
+      INTEGER :: ierr
+
+       nxc=nx
+       nyc=ny
+       nzc=nz
+
+       ! CALL ludcmp_testalloc
+
+!       IF(nyc /= nzc) PRINT*,'ny should be equal nz'
+
+!      xp, yp, zp = 0 : periodic
+       IF(xp.eq.0) THEN
+         CALL p_lud(1,nxc)
+       ELSE
+         CALL nonp_lud(1,nxc)
+       ENDIF
+
+       IF(yp.eq.0) THEN
+         CALL p_lud(2,nyc)
+       ELSE
+         call nonp_lud(2,nyc)
+       ENDIF
+
+       IF(zp.eq.0) THEN
          CALL p_lud(3,nzc)
        ELSE 
          call nonp_lud(3,nzc)
        ENDIF
 
-!       CALL x_lud
+     END SUBROUTINE ludcmp_calculate
-!       CALL yz_lud
+
+
+     SUBROUTINE test_nonp_lud1(xx, coef)
+      INTEGER :: xx
+      REAL(KIND=8), DIMENSION(xx) :: aa
+      REAL(KIND=8), DIMENSION(xx), INTENT(OUT) :: coef
+       aa=3.
+       aa(1)=0.5 ; aa(2)=4.
+       aa(xx-1)=4. ; aa(xx)=0.5
+
+       CALL stdlu(aa,xx,coef)
+     END SUBROUTINE test_nonp_lud1
+
+
+     SUBROUTINE test_nonp_lud2(xx, coef)
+      INTEGER :: xx
+      REAL(KIND=8), DIMENSION(xx) :: aa
+      REAL(KIND=8), DIMENSION(xx), INTENT(OUT) :: coef
+       aa=5.5
+       aa(1)=2./11. ; aa(2)=10.
+       aa(xx-1)=10. ; aa(xx)=2./11.
+
+       CALL stdlu(aa,xx,coef)
+     END SUBROUTINE test_nonp_lud2
+
+
+     SUBROUTINE test_p_lud1(xx, coef1, coef2)
+      INTEGER :: xx
+      REAL(KIND=8) :: a
+      REAL(KIND=8), DIMENSION(xx), INTENT(OUT) :: coef1, coef2
+      a=3. ! first derivative
+
+      CALL ptdlu(a,xx,coef1,coef2) ! x-direction
+     END SUBROUTINE test_p_lud1
+
+
+     SUBROUTINE test_p_lud2(xx, coef1, coef2)
+      INTEGER :: xx
+      REAL(KIND=8) :: a
+      REAL(KIND=8), DIMENSION(xx), INTENT(OUT) :: coef1, coef2
+      a=11./2. ! second derivative
+
+      CALL ptdlu(a,xx,coef1,coef2) ! x-direction
+     END SUBROUTINE test_p_lud2
 
-     END SUBROUTINE ludcmp
 
      SUBROUTINE nonp_lud(xyz,xx)
       INTEGER :: i,xyz,xx
-      REAL, DIMENSION(xx) :: aa
+      REAL(KIND=8), DIMENSION(xx) :: aa
        aa=3.
        aa(1)=0.5 ; aa(2)=4.
        aa(xx-1)=4. ; aa(xx)=0.5
@@ -92,7 +222,7 @@
 
      SUBROUTINE p_lud(xyz,xx)
       INTEGER :: i,xyz,xx
-      REAL :: a
+      REAL(KIND=8) :: a
       a=3. ! first derivative
       IF (xyz.eq.1) CALL ptdlu(a,xx,lxf,wxf) ! x-direction
       IF (xyz.eq.2) CALL ptdlu(a,xx,lyf,wyf) ! y-direction
@@ -105,21 +235,23 @@
 
      SUBROUTINE stdlu(a,n,l)
       INTEGER :: n
-      REAL :: a(n),l(n)
+      REAL(KIND=8), INTENT(IN)  :: a(n)
-      REAL :: d
+      REAL(KIND=8), INTENT(OUT) :: l(n)
+      REAL(KIND=8) :: d
       INTEGER :: i
-      l(1)=1.0/a(1)
+      l(1)=1.0d0/a(1)
       DO i=2,n
         d=a(i)-l(i-1)
-        l(i)=1.0/d
+        l(i)=1.0d0/d
       ENDDO
      END SUBROUTINE stdlu
 
      SUBROUTINE ptdlu(a,n,l,w)
       INTEGER :: n
-      REAL :: a,l(n),w(n)
+      REAL(KIND=8), INTENT(OUT) :: a
+      REAL(KIND=8), INTENT(OUT) :: l(n),w(n)
       INTEGER :: i
-      REAL :: aa(n),d
+      REAL(KIND=8) :: aa(n),d
 
       DO i=1,n-1
         aa(i)=a
@@ -141,18 +273,19 @@
       l(n)=1./d
      END SUBROUTINE ptdlu
 
-     SUBROUTINE dfnonp(n,h,x,dx,nd,dir)
+
-      INTEGER,INTENT(IN) :: n,nd,dir
+     SUBROUTINE rhs1np(n,h,x,dx,nd)
-      REAL,INTENT(IN) :: h
+      INTEGER,INTENT(IN) :: n,nd
-      REAL,INTENT(IN),DIMENSION(nd,n) :: x
+      REAL(KIND=8),INTENT(IN) :: h
-      REAL,INTENT(OUT),DIMENSION(nd,n) :: dx
+      REAL(KIND=8),INTENT(IN),DIMENSION(nd,n) :: x
+      REAL(KIND=8),INTENT(OUT),DIMENSION(nd,n) :: dx
       INTEGER :: i,j
-      REAL :: r1,r2,r3,a,b,c,h1,t1,t2,t3,t4
+      REAL(KIND=8) :: r1,r2,r3,a,b,c,h1,t1,t2,t3,t4
 
-      h1=1./h
+      h1=1.d0/h
 
-      r1=7./3.
+      r1=7.d0/3.d0
-      r2=1./12.
+      r2=1.d0/12.d0
       r3=3.
       a=-1.25
       b=1.
@@ -178,6 +311,20 @@
         dx(j,i)=h1*(r1*t1+r2*t2)
        ENDDO
       ENDDO
+
+     END SUBROUTINE rhs1np
+
+
+     SUBROUTINE dfnonp(n,h,x,dx,nd,dir)
+      INTEGER,INTENT(IN) :: n,nd,dir
+      REAL(KIND=8),INTENT(IN) :: h
+      REAL(KIND=8),INTENT(IN),DIMENSION(nd,n) :: x
+      REAL(KIND=8),INTENT(OUT),DIMENSION(nd,n) :: dx
+      INTEGER :: i,j
+      REAL(KIND=8) :: r1,r2,r3,a,b,c,h1,t1,t2,t3,t4
+
+      CALL rhs1np (n,h,x,dx,nd)
+
       IF (dir.eq.1) CALL tdslv(dx,n,lxf,nd) ! x-direction
       IF (dir.eq.2) CALL tdslv(dx,n,lyf,nd) ! y-direction
       IF (dir.eq.3) CALL tdslv(dx,n,lzf,nd) ! z-direction
@@ -185,11 +332,14 @@
 
      SUBROUTINE dfp(n,h,x,dx,nd,dir)
       INTEGER,INTENT(IN) :: n,nd,dir
-      REAL,INTENT(IN) :: h
+      REAL(KIND=8),INTENT(IN) :: h
-      REAL,INTENT(IN),DIMENSION(nd,n) :: x
+      REAL(KIND=8),INTENT(IN),DIMENSION(nd,n) :: x
-      REAL,INTENT(OUT),DIMENSION(nd,n) :: dx
+      REAL(KIND=8),INTENT(OUT),DIMENSION(nd,n) :: dx
       INTEGER :: i,j
-      REAL :: r1,r2,h1
+      REAL(KIND=8) :: r1,r2,h1
+
+
+      ! print *, "dfnonp received (nd,n)", nd, n
 
       h1=1./h
       r1=7./3.
@@ -221,10 +371,10 @@
 
      SUBROUTINE ptdslv(r,n,l,w,nd)
       INTEGER,INTENT(IN) :: n,nd
-      REAL,INTENT(INOUT),DIMENSION(nd,n) :: r
+      REAL(KIND=8),INTENT(INOUT),DIMENSION(nd,n) :: r
-      REAL,INTENT(IN),DIMENSION(:) :: l,w
+      REAL(KIND=8),INTENT(IN),DIMENSION(n) :: l,w
       INTEGER i,j
-      REAL, DIMENSION(nd) :: sum
+      REAL(KIND=8), DIMENSION(nd) :: sum
       DO j=1,nd
         sum(j)=w(1)*r(j,1)
         r(j,1)=r(j,1)*l(1)
@@ -249,11 +399,13 @@
 
      SUBROUTINE d2fp(n,h,x,dx,nd,dir)
       INTEGER,INTENT(IN) :: n,nd,dir
-      REAL,INTENT(IN) :: h
+      REAL(KIND=8),INTENT(IN) :: h
-      REAL,INTENT(IN),DIMENSION(nd,n) :: x
+      REAL(KIND=8),INTENT(IN),DIMENSION(nd,n) :: x
-      REAL,INTENT(OUT),DIMENSION(nd,n) :: dx
+      REAL(KIND=8),INTENT(OUT),DIMENSION(nd,n) :: dx
       INTEGER :: i,j
-      REAL :: h2,r1,r2,t1,t2
+      REAL(KIND=8) :: h2,r1,r2,t1,t2
+
+
       h2=1./(h*h)
       r1=6.
       r2=3./8.
@@ -312,10 +464,10 @@
 
      SUBROUTINE tdslv(r,n,l,nd)
       INTEGER,INTENT(IN) :: n,nd
-      REAL,INTENT(INOUT),DIMENSION(nd,n) :: r
+      REAL(KIND=8),INTENT(INOUT),DIMENSION(nd,n) :: r
-      REAL,INTENT(IN),DIMENSION(:) :: l
+      REAL(KIND=8),INTENT(IN),DIMENSION(n) :: l
       INTEGER i,j
-      REAL t1
+      REAL(KIND=8) t1
       DO j=1,nd
         r(j,1)=r(j,1)*l(1)
       ENDDO
@@ -334,11 +486,12 @@
 
      SUBROUTINE d2fnonp(n,h,x,dx,nd,dir)
       INTEGER,INTENT(IN) :: n,nd,dir
-      REAL,INTENT(IN) :: h
+      REAL(KIND=8),INTENT(IN) :: h
-      REAL,INTENT(IN),DIMENSION(nd,n) :: x
+      REAL(KIND=8),INTENT(IN),DIMENSION(nd,n) :: x
-      REAL,INTENT(OUT),DIMENSION(nd,n) :: dx
+      REAL(KIND=8),INTENT(OUT),DIMENSION(nd,n) :: dx
       INTEGER :: i,j
-      REAL :: h2,r1,r2,r3,a,b,c,e,t1,t2
+      REAL(KIND=8) :: h2,r1,r2,r3,a,b,c,e,t1,t2
+
 
       h2=1./(h*h)
       r1=6.
@@ -383,4 +536,4 @@
 
      END SUBROUTINE d2fnonp
 
-    END MODULE Compact
+    END MODULE m_compact

code/m_parameters.f90

@@ -0,0 +1,157 @@
+module m_parameters
+
+  implicit none
+
+     ! Domain Parameter
+     integer :: nxp,nyp,nzp
+     integer :: nx
+     real :: hxp,hyp,hzp
+     real :: l_0
+     real :: hx
+
+     ! Transport Properties
+     real :: vis,sc,diff
+     real :: d_turb
+     integer :: d_mode
+     real :: scp,prp,lep,vis0p,rod
+     real :: lewis, le_a, le_x
+
+     ! Chemistry Properties
+     real :: prof_wr,min_wr,min_fsd,min_c,refwr
+     real :: pre,ac,bc,c_cut,c_ref
+     real :: lambda1, lambda2, beta1, hrp, beta
+     character(100) :: reaction_type
+
+
+     ! Constants
+     real, parameter :: pi=3.14159265358979323846d0
+     ! real, parameter :: pi=acos(-1.d0)
+     real, parameter :: me=1.00e-20
+
+     ! Flame Control
+     real :: minf,tar_lo,u0,ctmp,lo_flm=0.
+     real :: pflame,pflold,oldu
+
+     ! Solver Control
+     integer :: nsp, nrxn
+     integer :: ncyc=0,int_pr
+     integer :: fctrl_species=1
+     real :: absolute_tolerence=1e-8
+     real :: dt,tf,t_now,t_uf,dt_uf
+
+     ! Input File
+     LOGICAL :: read_itape, read_stdin
+     CHARACTER(100) :: itape_name
+
+     integer :: table_size
+     real :: rate_relaxation (2,100)
+
+     CHARACTER(LEN=8) :: cdum
+     INTEGER :: itape=300, otape=301
+
+contains
+
+     SUBROUTINE SET_CHEMISTRY
+
+       reaction_type = "twostep"
+       lambda1 = 100000.
+       lambda2 = lambda1 * 63.6
+       beta1 = 9.1
+       hrp = 0.8136765
+
+     END SUBROUTINE SET_CHEMISTRY
+
+     SUBROUTINE READ_INTRO
+
+     IMPLICIT NONE
+
+     INTERFACE READ_PARAMETER
+
+     SUBROUTINE READ_INT (x)
+        integer :: x
+     END SUBROUTINE READ_INT
+
+     SUBROUTINE READ_REAL (x)
+        real :: x
+     END SUBROUTINE READ_REAL
+
+     END INTERFACE
+
+       INTEGER :: i
+
+       CALL SET_CHEMISTRY
+
+       IF (read_stdin) THEN
+         itape=5
+       ELSE
+         OPEN(itape,FILE=itape_name)
+       END IF
+
+       OPEN(otape,FILE='otape')
+        CALL READ_PARAMETER (nx)       ! n grid points
+        CALL READ_PARAMETER (l_0)      ! domain length
+        CALL READ_PARAMETER (int_pr)   ! print interval
+        CALL READ_PARAMETER (tar_lo)   ! target flame location
+        CALL READ_PARAMETER (dt)       ! time step size
+        CALL READ_PARAMETER (sc)       ! Schmidt Number
+        CALL READ_PARAMETER (vis)      ! viscosity
+
+        CALL READ_PARAMETER (pre)      ! 1-step chemistry
+        CALL READ_PARAMETER (ac)       ! 1-step chemistry
+        CALL READ_PARAMETER (bc)       ! 1-step chemistry
+
+        CALL READ_PARAMETER (u0)       ! inlet velocity
+        CALL READ_PARAMETER (tf)       ! end time
+        CALL READ_PARAMETER (dt_uf)    ! flame control time step
+        CALL READ_PARAMETER (ctmp)     ! unburned c value
+
+        CALL READ_PARAMETER (c_cut)    ! reaction rate profile parameter
+        CALL READ_PARAMETER (c_ref)    ! reaction rate profile parameter
+        CALL READ_PARAMETER (min_wr)   ! reaction rate profile parameter
+        CALL READ_PARAMETER (prof_wr)  ! reaction rate profile parameter
+
+        CALL READ_PARAMETER (lewis)    ! Lewis number
+
+        CALL READ_PARAMETER (lambda1)  ! 2-step chemistry
+        CALL READ_PARAMETER (lambda2)  ! 2-step chemistry
+        CALL READ_PARAMETER (beta1)    ! 2-step chemistry
+        CALL READ_PARAMETER (hrp)      ! 2-step chemistry
+        CALL READ_PARAMETER (le_a)     ! 2-step chemistry
+        CALL READ_PARAMETER (le_x)     ! 2-step chemistry
+
+        CALL READ_PARAMETER (nsp)      ! number of species
+        CALL READ_PARAMETER (nrxn)     ! number of reactions
+
+        CALL READ_PARAMETER (d_turb)   ! turbulent diffusivity
+        CALL READ_PARAMETER (d_mode)   ! diffusivity formula selector
+
+        ! READ(itape,*) cdum,table_size
+        ! WRITE(otape,*) cdum,table_size
+
+        ! do i = 1:table_size
+        ! READ(itape,*) cdum,rate_relaxation(1,i),rate_relaxation(2,i)
+        ! WRITE(otape,*) cdum,rate_relaxation(1,i),rate_relaxation(2,i)
+        ! end do
+
+       IF (.not.read_stdin) THEN
+         CLOSE(itape)
+       END IF
+
+       l_0=l_0*pi
+
+       hx=l_0/REAL(nx)
+       cdum='hx'
+       WRITE(otape,*) cdum,hx
+       cdum='Ta/Tu'
+       WRITE(otape,*) cdum,ac
+       cdum='Tb/Tu'
+       WRITE(otape,*) cdum,bc+1
+       diff=vis/sc
+       cdum='diff'
+       WRITE(otape,*) cdum,diff
+
+       refwr=pre*1.*exp(-ac/(1.+bc*c_ref))
+
+     END SUBROUTINE READ_INTRO
+
+end module m_parameters

code/makefile

@@ -1,18 +1,37 @@
-flags = -Wall -O3 -fdefault-integer-8 -fdefault-double-8 -fdefault-real-8  -march=native 
+BLOCKSIZE?=32
+
+flags = -Wall -O3 -cpp -DBLOCKSIZE=$(BLOCKSIZE) -fdefault-integer-8 -fdefault-double-8 -fdefault-real-8  -march=native 
+
+ifdef CHECK
+flags += -fcheck=all
+endif
+
+ifdef DEBUG
+flags += -g
+endif
 
 compiler = gfortran
 
-ex  : test.o ysolve.mod compact.mod
+ex  : test.o read_parameter.o ysolve.o m_chemistry.o m_parameters.o m_compact.o
-	${compiler} -o ex test.o ysolve.o Compact.o 
+	${compiler} -o ex test.o read_parameter.o ysolve.o m_chemistry.o m_parameters.o m_compact.o
 
-test.o  : test.f90 ysolve.mod
+test.o  : test.f90 ysolve.mod m_compact.mod 
 	${compiler} -c ${flags} test.f90
 
-ysolve.mod: ysolve.f90 compact.mod
+read_parameter.o  : read_parameter.f90 m_parameters.o
+	${compiler} -c ${flags} read_parameter.f90
+
+ysolve.o ysolve.mod : ysolve.f90 m_compact.mod m_chemistry.mod m_parameters.mod
 	${compiler} -c ${flags} ysolve.f90
 
-compact.mod: Compact.f90
+m_compact.o m_compact.mod : m_compact.f90
-	${compiler} -c ${flags} Compact.f90
+	${compiler} -c ${flags} m_compact.f90
+
+m_parameters.o m_parameters.mod : m_parameters.f90
+	${compiler} -c ${flags} m_parameters.f90
+
+m_chemistry.o m_chemistry.mod : m_chemistry.f90 m_parameters.mod
+	${compiler} -c ${flags} m_chemistry.f90
 
 test: ex
 	sh test.sh

code/read_parameter.f90

@@ -0,0 +1,32 @@
+subroutine read_int (x)
+  use m_parameters
+
+  implicit none
+
+  integer :: x
+  READ(itape,*) cdum,x
+  WRITE(otape,*) cdum,x
+
+end subroutine read_int
+
+subroutine read_real (x)
+  use m_parameters
+
+  implicit none
+
+  real :: x
+  READ(itape,*) cdum,x
+  WRITE(otape,*) cdum,x
+
+end subroutine read_real
+
+subroutine read_string (x)
+  use m_parameters
+
+  implicit none
+
+  character(len=*) :: x
+  READ(itape,*) cdum,x
+  WRITE(otape,*) cdum,x
+
+end subroutine read_string

code/run.sh

@@ -0,0 +1 @@
+ (TSTAMP=`date -Iseconds`; python sample/itape.py --default sample/4pi-IC1/itape | tee ${TSTAMP}-input.txt | ./ex - | tee ${TSTAMP}-output.txt)

code/sample/4pi-IC1/itape

@@ -2,18 +2,24 @@
        l_0     4.   !3[pi]   
     int_pr     400
     tar_lo     0.60
-        dt     0.001 
+        dt     0.0005 
         sc     0.75
        vis     0.020           !2pi(162)ref  !2pi(256)mid  !2pi(256)Max !2p(162)Low
        pre     2.10E+4       !185.75          !185.75       !316.88       !440.5       !102.69
         ac     26.7            !22.250          !22.25        !22.56        !22.3        !22.39
         bc     3.              !7.000           !7.0          !7.0          !7.0         !7.0
         u0     0.21
-        tf     80.
+        tf     200.
      dt_uf     0.4
       ctmp     0.        !1.0e-14
      c_cut     0.001 !0.012     ! c < c_cut -> wrate = 0.
      c_ref     0.01  !0.003     ! 
     min_wr     0.        ! 5.0e-14
    prof_wr     1.0       
-
+     lewis     1.0
+       lambda1    6.25E+4       7.7E+4       5E+4      
+       lambda2    318E+4   445.2E+4   489.72E+4   
+       beta1   5
+       hrp   0.75
+       le_a  1.
+       le_x  0.20

code/sample/itape.py

@@ -0,0 +1,89 @@
+#!/bin/python
+
+import os, sys
+
+import argparse
+
+from pprint import pprint
+
+parser = argparse.ArgumentParser()
+file_parser = argparse.ArgumentParser()
+
+# Default values for all parameters
+raw_defaults = (
+'''\
+nx        512
+l_0       4.
+int_pr    400
+tar_lo    0.60
+dt        0.0005
+sc        0.75
+vis       0.020
+pre       2.10E+4
+ac        26.7
+bc        3.
+u0        0.21
+tf        100.
+dt_uf     0.4
+ctmp      0.
+c_cut     0.001
+c_ref     0.01
+min_wr    0.
+prof_wr   1.0
+lewis     1.0
+lambda1   6.25E+4
+lambda2   318E+4
+beta1     5
+hrp       0.75
+le_a      1.
+le_x      0.20
+nsp       2
+nrxn      1
+d_turb    0
+d_mode    0
+'''
+).split()
+
+ordered_args = raw_defaults[::2]
+
+ordered_values = []
+for v in raw_defaults[1::2]:
+    try:
+        ordered_values.append(int(v))
+    except ValueError:
+        ordered_values.append(float(v))
+
+param_dict = dict(zip(ordered_args, ordered_values))
+
+
+file_parser.add_argument("--default", help="file containing override default values")
+
+file_args, others = file_parser.parse_known_args()
+
+if file_args.default is not None:
+    with open(file_args.default) as df:
+        try:
+            for line in df:
+                k, v = line.split()[:2]
+                if k in param_dict:
+                    param_dict[k] = v
+        except ValueError:
+          pass
+
+
+for k, v in param_dict.iteritems():
+    parser.add_argument("--"+k, default=v)
+
+args = parser.parse_args(others)
+
+params = vars(args)
+
+line_new = '{:>12}     {}'
+
+for a in ordered_args:
+    print line_new.format(a, params[a])
+
+# for i in range(15,60,5):
+    # with open("itape", 'w') as itape:
+        # itape.write(format.format(float(i)/10.0))
+    # os.system("../../ex &>"+str(i)+".txt")

code/test.f90

@@ -6,5 +6,12 @@
     PROGRAM test
      USE ysolve
      IMPLICIT NONE
+     CALL parse
+
+     CALL init_solver
+
      CALL solve
+
+     CALL finalize_solver
+
     END PROGRAM test

code/test.sh

@@ -1,4 +1,4 @@
 cd sample/4pi-IC1
 ls
-../../ex
+../../ex itape
-diff sfield.dat sfield.diff
+diff -y sfield.dat sfield.diff

code/ysolve.f90

@@ -1,76 +1,202 @@
     MODULE ysolve
-     USE Compact
+     USE m_compact
+     USE m_parameters
+     USE m_chemistry
+
      IMPLICIT NONE
-     PRIVATE
+
-     REAL, PARAMETER :: pi=3.14159265358979323846
+     REAL, DIMENSION(:), ALLOCATABLE :: u
-     REAL :: hx,dt,vis,sc,diff,pre,ac,bc,tf,t_now,t_uf,dt_uf
+     REAL, DIMENSION(:), ALLOCATABLE :: inletbc
-     REAL :: c_cut,c_ref,refwr,minf,tar_lo,u0,ctmp,l_0,lo_flm=0.
+     REAL, DIMENSION(:,:), ALLOCATABLE :: y1,y2,yf,yold
-     REAL :: er_lof=0., erdot=0.,min_wr,prof_wr
-     REAL :: pflame,pflold,oldu
-     INTEGER :: ncyc=0,int_pr,nx
-     REAL, DIMENSION(:,:,:), ALLOCATABLE :: u,y1,y2,yf
      REAL, DIMENSION(:), ALLOCATABLE :: uxt,duxt
+     REAL, DIMENSION(:), ALLOCATABLE :: dm
 
-     INTEGER :: istage
-     REAL, DIMENSION(5) :: a=(/ 970286171893./4311952581923.,    &
-                                6584761158862./12103376702013.,  &
-                                2251764453980./15575788980749.,  &
-                                26877169314380./34165994151039., &
-                                0. /),                           &
-                           b=(/ 1153189308089./22510343858157.,  &
-                                1772645290293./4653164025191.,   &
-                               -1672844663538./4480602732383.,   &
-                                2114624349019./3568978502595.,   &
-                                5198255086312./14908931495163. /)
-
-     PUBLIC :: solve
     CONTAINS
 !------------------------------------------------------------------------
+     SUBROUTINE parse
+        CHARACTER(100) :: num1char
+
+        ! Good. One command line argument provided.
+        IF(COMMAND_ARGUMENT_COUNT().EQ.1)THEN
+          CALL GET_COMMAND_ARGUMENT(1,num1char)
+
+        ! Bad. No command line argument provided. Print help message.
+        ELSE IF(COMMAND_ARGUMENT_COUNT().EQ.0)THEN
+          WRITE(*,*)'ERROR, NO COMMAND-LINE ARGUMENT'
+          num1char = "-h"
+
+        ! Bad. Multiple command line arguments provided. Print help message.
+        ELSE
+          WRITE(*,*)'ERROR, TOO MANY COMMAND-LINE ARGUMENTS( > 2 )'
+          num1char = "-h"
+        ENDIF
+
+
+        ! Option or STDIN
+        IF(num1char(1:1) == "-") THEN
+
+          if(num1char=="-h") then
+            write(*,'(a)') "usage: ex [-h] -|input_file"
+            write(*,'(a)') ""
+            write(*,'(a)') "positional arguments:"
+            write(*,'(a)') " -                 read from std_in"
+            write(*,'(a)') " input_file        input file name"
+            write(*,'(a)') ""
+            write(*,'(a)') "optional arguments:"
+            write(*,'(a)') " -h                show this help message and exit"
+            stop
+
+          else if (num1char=="-") then
+            read_stdin = .true.
+
+          else
+            WRITE(*,*)'ERROR, UNSUPPORTED OPTION ', trim(num1char), '. STOPPING'
+            STOP
+          end if
+
+        END IF
+
+
+        IF(read_stdin)THEN
+          WRITE(*,*) "Read from STDIN"
+        ELSE 
+          itape_name = num1char
+          WRITE(*,*) "Read from file " // num1char
+        END IF
+
+     END SUBROUTINE parse
+
+!------------------------------------------------------------------------
+
+     REAL FUNCTION residual (x0, x1)
+
+       REAL, DIMENSION(:,:) :: x0, x1
+
+       residual = sum(abs(x0-x1))
+
+     END FUNCTION residual
+
+     LOGICAL FUNCTION converged (x0, x1)
+
+       REAL, DIMENSION(:,:) :: x0, x1
+       REAL :: r = 0.
+
+       r = residual(x0 ,x1)
+       converged = r < absolute_tolerence
+
+     END FUNCTION converged
+
+!------------------------------------------------------------------------
+
      SUBROUTINE solve
-      INTEGER :: i,j,k,savenum
+       IF (d_mode == 0) THEN
+         CALL solve_with_diffusivity (update_dm)
+       ELSE IF (d_mode == 1) THEN
+         CALL solve_with_diffusivity (update_dm_cd)
+       ELSE IF (d_mode == 2) THEN
+         CALL solve_with_diffusivity (update_dm_sutherland)
+       ELSE
+         WRITE(*,*)'ERROR, UNSUPPORTED DIFFUSIVITY ', d_mode, '. STOPPING'
+         STOP
+       END IF
+     END SUBROUTINE solve
+
+     SUBROUTINE solve_with_diffusivity (calc_diff)
+      INTERFACE
+       SUBROUTINE calc_diff(r1)
+        REAL, INTENT(IN),  DIMENSION(:,:) :: r1
+       END SUBROUTINE calc_diff
+      END INTERFACE
+
+       IF (reaction_type == "onestep") THEN
+         CALL solve_ (fonestep, SET_IC_ONESTEP, calc_diff)
+       ELSE IF (reaction_type == "twostep") THEN
+         CALL solve_ (fns, SET_IC_TWOSTEP, calc_diff)
+       ELSE
+         WRITE(*,*)'ERROR, UNSUPPORTED REACTION ', trim(reaction_type), '. STOPPING'
+         STOP
+       END IF
+     END SUBROUTINE solve_with_diffusivity
+
+     SUBROUTINE solve_ (rhs, SET_IC, calc_diff)
+
+      INTERFACE
+       SUBROUTINE rhs(r1,f)
+        REAL, INTENT(IN),  DIMENSION(:,:) :: r1
+        REAL, INTENT(OUT), DIMENSION(:,:) :: f
+       END SUBROUTINE rhs
+
+       SUBROUTINE SET_IC
+       END SUBROUTINE SET_IC
+
+       SUBROUTINE calc_diff(r1)
+        REAL, INTENT(IN),  DIMENSION(:,:) :: r1
+       END SUBROUTINE calc_diff
+      END INTERFACE
+
+      INTEGER :: i
       REAL :: pflame,pflold,delf=0.
-
+      REAL :: residue = 0.
-      CALL READ_INTRO
-
-      CALL ludcmp(nx,5,5,1,0,0)
 
       CALL SET_IC
 
-      t_now=0.; t_uf=0.
+      pflame=0.
+      DO i=1,nx
+       pflame=pflame+y1(i,fctrl_species)*hx
+      ENDDO
 
-      DO 
-       CALL SET_BC  
-       CALL RK4
-       IF(t_now.ge.tf) EXIT
-       IF(t_uf.ge.dt_uf) THEN
-        pflold=pflame
-        pflame=0.
-        DO i=1,nx
-         pflame=pflame+y1(i,1,1)*hx
-         uxt(i)=y1(i,1,1)
-        ENDDO
-        CALL dfnonp(nx,hx,uxt,duxt,1,1)
-        delf=1./MAXVAL(ABS(duxt))
-        oldu=u(1,1,1)
-        u=u+0.5*(hx*REAL(nx)*tar_lo-pflame)+0.5*(pflold-pflame)
-        t_uf=0.
-        WRITE(*,'(a3,f8.3,a10,f6.3,a10,f6.3,a10,f7.4,a10,f7.4,a8,f7.4)') &
-         ' T:',t_now,' // Tar_L:',l_0*tar_lo,' // cur_L:',pflame/hx/REAL(nx)*l_0, &
-         ' // Old_U:',oldu,' // New_U:',u(1,1,1),' // L_f:',delf
 
-!        WRITE(*,'(a7,f7.4,a7,f7.4,a10,f7.4)') ' cur_U:',oldu,' // dU:',u(1,1,1)-oldu,' // new_U:',u(1,1,1)
+      t_now=0.
-!        WRITE(*,*)
+      t_uf=0.
-       ENDIF
+
+      DO
+
+       IF (t_now.ge.tf) EXIT
+
+       residue = residual(yold, y1)
+
+       IF (converged(yold, y1)) EXIT
+
        ncyc=ncyc+1
        t_uf=t_uf+dt
        t_now=t_now+dt
+
+       yold = y1
+
+       CALL update_chemistry(t_now)
+
+       CALL SET_BC
+
+       CALL RK4(rhs, calc_diff)
+
+       IF (t_uf.ge.dt_uf) THEN
+         pflold=pflame
+         pflame=0.
+         DO i=1,nx
+           pflame=pflame+y1(i,fctrl_species)*hx
+           uxt(i)= 1. - y1(i,fctrl_species)
+         ENDDO
+
+         CALL dfnonp(nx,hx,uxt,duxt,1,1)
+         delf=1./MAXVAL(ABS(duxt))
+
+         oldu=u(1)
+         u=u+0.5*(hx*REAL(nx)*tar_lo-pflame)+0.5*(pflold-pflame)
+         t_uf=0.
+
+         WRITE(*,'(a3,f8.3,a10,f6.3,a10,f6.3,a10,f7.4,a10,f7.4,a8,f7.4,a6,e10.4)') &
+          ' T:',t_now,' // Tar_L:',l_0*tar_lo,' // cur_L:',pflame/hx/REAL(nx)*l_0, &
+          ' // Old_U:',oldu,' // New_U:',u(1),' // L_f:',delf,' // R:',residue
+
+       ENDIF
+
        IF (MOD(ncyc,int_pr).eq.0) THEN
 !        WRITE(*,'(a2,f8.3,a9,f10.7,a11,i6,a7,f9.5)') &
 !           'T:',t_now,'  //  dT:',dt,'  //  NCYC:',ncyc,'  // U:',u(1,1,1)
-        ENDIF
+       ENDIF
       ENDDO
 
-      CALL write_sd
+      ! CALL write_sd
       CALL write_pre
       CALL save_final_field
 
@@ -78,40 +204,41 @@
       WRITE(*,*) 'Fin.'
       WRITE(*,*)
 
-     END SUBROUTINE solve
+     END SUBROUTINE solve_
 
      SUBROUTINE write_sd
-      REAL :: c,yr,wrate,dely,sdr,sdd,sd,uu,onelw,dd
+      REAL :: c,yr,theta,wrate,dely,sdr,sdd,sd,uu,onelw,dd
       INTEGER :: i,j,k,nd
-      REAL, DIMENSION(1,nx) :: ux,dux,d2ux
+      REAL, DIMENSION(2,nx) :: ux,dux,d2ux
       REAL, DIMENSION(10,nx) :: sav
-        
+
       sav=0.
 
 !      refwr=pre*1.*exp(-ac/(1.+bc*c_cut))
 !      minf=exp((c_ref-c_cut)*prof_wr)
       refwr=pre*1.*exp(-ac/(1.+bc*c_ref))   
-     
+
-               
+
       WRITE(500,*) 'VARIABLES = "X","Yr","C","U","Wrate","|DEL(Y)|","Sdr"'
       WRITE(500,*) ' "Sdd","Sd","(1/C)/(dC/dx)","DIV(rho*Dmu*Gra(C))"'
 
       DO i=1,nx
-       ux(1,i)=y1(i,1,1)                               ! Yr
+       ux(1,i)=y1(i,1)                               ! Yr
+       ux(2,i)=y1(i,2)                               ! T
        IF (ux(1,i).gt.1.) ux(1,i)=1.                                  
       ENDDO
-          
+
-      nd=1
+      nd=2
-      CALL dfnonp(nx,hx,ux(1,:),dux(1,:),nd,1)  
+      CALL dfnonp(nx,hx,ux(:,:),dux(:,:),nd,1)
-      nd=1 
+      CALL d2fnonp(nx,hx,ux(:,:),d2ux(:,:),nd,1)
-      CALL d2fnonp(nx,hx,ux(1,:),d2ux(1,:),nd,1) 
 
       DO i=1,nx
        yr=ux(1,i)
+       theta=ux(2,i)
        c=1.-yr
        IF (c.lt.0.) c=0.
 
-       wrate=pre*yr*exp(-ac/(1.+bc*(1.-yr)))   !wrate
+       wrate=pre*yr*exp(-ac/(1.+bc*(theta)))   !wrate
 
 !       IF (c.le.c_cut) THEN
 !        wrate=min_wr
@@ -119,22 +246,22 @@
 !         ((exp((c-c_cut)*prof_wr)-minf)/(1.-minf)*(refwr-min_wr))+min_wr
 !       ENDIF
 
-       IF (c.le.c_ref) THEN
+       IF (theta.le.c_ref) THEN
          wrate=min_wr
-        IF (c.gt.c_cut) wrate=((refwr-min_wr)*exp(prof_wr*(c-c_ref))+ &
+        IF (theta.gt.c_cut) wrate=((refwr-min_wr)*exp(prof_wr*(theta-c_ref))+ &
         min_wr-refwr*exp(prof_wr*(c_cut-c_ref)))/(1.-exp(prof_wr*(c_cut-c_ref)))
        ENDIF
 
        dely=ABS(dux(1,i))
        sdr=wrate/dely
-       sdd=-diff*d2ux(1,i)/dely
+       sdd=-dm(i)*d2ux(1,i)/dely
        IF (dely.eq.0.) THEN
         sdr=0.; sdd=0.
        ENDIF
        sd=sdr+sdd
-       uu=u(1,1,1)
+       uu=u(1)
        onelw=(-dux(1,i))/c
-       dd=-diff*d2ux(1,i)
+       dd=-dm(i)*d2ux(1,i)
        if (c.eq.0.) onelw=0.
 
        sav(1,i)=sav(1,i)+yr
@@ -155,48 +282,63 @@
      END SUBROUTINE write_sd
 
      SUBROUTINE save_final_field
-      INTEGER :: i,j,k
+      INTEGER :: i
 
       OPEN (305,FILE='sfield.bin',form='unformatted',status='unknown')
       DO i=1,nx
-       WRITE (305) y1(i,1,1)
+       WRITE (305) y1(i,1)
       ENDDO
       CLOSE (305)
 
       OPEN (305,FILE='sfield.dat')
       DO i=1,nx
-       WRITE (305,'(e30.20)') y1(i,1,1)
+       WRITE (305,*) (/ i*hx, y1(i,:) /)
       ENDDO
       CLOSE (305)
 
      END SUBROUTINE save_final_field
 
      SUBROUTINE write_pre
-      REAL :: yr,c,dy,maxdy=0.,del_f
+      REAL :: theta,yr,c,dy,maxdy=0.,del_f
-      REAL :: S_L=0.,wrate
+      REAL :: S_L=0.,wrate,wrate1,wrate2
+      REAL :: ya,yx
       INTEGER :: i
       REAL, DIMENSION(1,nx) :: ux, dux
-      
-      DO i=1,nx
-       yr=y1(i,1,1)
-       c=1.-yr
-       ux(1,i)=yr
 
-       wrate=pre*yr*exp(-ac/(1.+bc*(1.-yr)))   !wrate
+      if ( reaction_type == "onestep" ) then
-!       IF (c.le.c_cut) THEN
+
-!        wrate=min_wr
+      DO i=1,nx
-!        IF (c.gt.c_ref) wrate= &
+       yr=y1(i,1)
-!         ((exp((c-c_cut)*prof_wr)-minf)/(1.-minf)*(refwr-min_wr))+min_wr
+       theta=y1(i,2)
-!       ENDIF
+
-       IF (c.le.c_ref) THEN
+       ux(1,i)=theta
-         wrate=min_wr
+
-        IF (c.gt.c_cut) wrate=((refwr-min_wr)*exp(prof_wr*(c-c_ref))+ &
+        wrate=rate_1step(yr, theta)
-        min_wr-refwr*exp(prof_wr*(c_cut-c_ref)))/(1.-exp(prof_wr*(c_cut-c_ref)))
-       ENDIF
 
        S_L=S_L+wrate*hx
       ENDDO
-      
+
+      else if ( reaction_type == "twostep" ) then
+
+      DO i=1,nx
+       ya=y1(i,1)
+       yx=y1(i,2)
+       theta=y1(i,3)
+
+       ux(1,i)=theta
+
+        wrate1=rate1_2step(ya, yx, theta)
+
+        wrate2=rate2_2step(yx, theta)
+
+       S_L=S_L+wrate1*hx
+      ENDDO
+
+      else
+          WRITE(*,*) 'ERROR, UNDEFINED REACTION TYPE ', reaction_type
+          stop
+      end if
+
       WRITE(*,'(a31,e14.8)') ' INTEGRAL( Wrate x dx ) => Sc :',S_L
 
       CALL dfnonp(nx,hx,ux(1,:),dux(1,:),1,1)
@@ -206,253 +348,338 @@
        maxdy=MAX(ABS(dy),maxdy)
       ENDDO
       del_f=1./maxdy
-       
+
-      WRITE(*,'(a13,e14.8,a25,e14.8)') ' Grid size : ',hx,' / Laminar flame speed : ',u(1,1,1)
+      WRITE(*,'(a13,e14.8,a25,e14.8)') ' Grid size : ',hx,' / Laminar flame speed : ',u(1)
       WRITE(*,'(a19,e14.8,a3,f9.5,a20)') &
        ' Flame thickness : ',del_f,' / ',del_f/hx,' grids in the flame.'
       WRITE(*,*)
      END SUBROUTINE write_pre
 
      SUBROUTINE SET_BC
-      
+      INTEGER :: i
-      y1(1,1,1)=1.-ctmp
+
+      DO i = 1, nsp
+       y1(1,i) = inletbc(i)
+      END DO
      END SUBROUTINE SET_BC
 
-     SUBROUTINE SET_IC
+     SUBROUTINE BASE_FLAME_PROFILE(x)
       INTEGER :: i, ifl, si
-      REAL :: xi
+      REAL :: x(nx)
 
+      ! initial flame thickness (0.2 pi n grids)
       ifl=INT(2.0*pi*0.1/hx)
 
-      u=u0; si=INT(nx*(1.-tar_lo))
+      ! initialize velocity field
+      u=u0;
 
+      ! initial flame center (grid index)
+      si=INT(nx*(1.-tar_lo))
+
+      ! initialize Yr field
       DO i=1,nx
        IF(i< nx-(si+ifl/2)) THEN
-        xi=0.+ctmp
+        x(i)=0.+ctmp
        ELSE IF(i> nx-(si-ifl/2)) THEN
-        xi=1.
+        x(i)=1.
        ELSE
-        xi=0.5+REAL(i-nx+si)/REAL(ifl)
+        x(i)=0.5+REAL(i-nx+si)/REAL(ifl)
        ENDIF
+      END DO
 
-!       y1(i,1,1)=(1.-xi)*1.
+      ! initialize species field
-       y1(i,1,1)=(1.-xi)   ! reactant mass fraction
+      ! max_ysum = maxval(y1(:,1) + y1(:,2))
-      ENDDO
+      ! y1(:,1) = y1(:,1) / maxval(y1(:,1) + y1(:,2))
+      ! y1(:,2) = y1(:,2) / maxval(y1(:,1) + y1(;,2))
+      ! y1(:,1) = y1(:,1) / (y1(1,1) + y1(1,2))
+      ! y1(:,2) = y1(:,2) / (y1(1,1) + y1(1,2))
 
-      pflame=0.
+     END SUBROUTINE BASE_FLAME_PROFILE
-      DO i=1,nx
-       pflame=pflame+y1(i,1,1)*hx
-      ENDDO
 
-     END SUBROUTINE SET_IC
+     SUBROUTINE SET_IC_ONESTEP
+      INTEGER :: i
+      REAL :: xi
+      REAL :: x(nx)
 
-     SUBROUTINE READ_INTRO
+      CALL BASE_FLAME_PROFILE(x)
-       CHARACTER(LEN=8) :: cdum
-       INTEGER :: itape=300,otape=301,ierr
 
+        DO i=1,nx
+         xi = x(i)
 
-       OPEN(itape,FILE='itape')
+         y1(i,1)=(1.-xi)   ! reactant mass fraction
-       OPEN(otape,FILE='otape')
+         y1(i,2)=(xi)   ! reactant mass fraction
-        READ(itape,*) cdum,nx
+        ENDDO
-        WRITE(otape,*) cdum,nx
+     END SUBROUTINE SET_IC_ONESTEP
-        READ(itape,*) cdum,l_0
-        WRITE(otape,*) cdum,l_0
-        READ(itape,*) cdum,int_pr
-        WRITE(otape,*) cdum,int_pr
-        READ(itape,*) cdum,tar_lo
-        WRITE(otape,*) cdum,tar_lo
-        READ(itape,*) cdum,dt
-        WRITE(otape,*) cdum,dt
-        READ(itape,*) cdum,sc
-        WRITE(otape,*) cdum,sc
-        READ(itape,*) cdum,vis
-        WRITE(otape,*) cdum,vis
-        READ(itape,*) cdum,pre
-        WRITE(otape,*) cdum,pre
-        READ(itape,*) cdum,ac
-        WRITE(otape,*) cdum,ac
-        READ(itape,*) cdum,bc
-        WRITE(otape,*) cdum,bc
-        READ(itape,*) cdum,u0
-        WRITE(otape,*) cdum,u0
-        READ(itape,*) cdum,tf
-        WRITE(otape,*) cdum,tf
-        READ(itape,*) cdum,dt_uf
-        WRITE(otape,*) cdum,dt_uf
-        READ(itape,*) cdum,ctmp
-        WRITE(otape,*) cdum,ctmp
-        READ(itape,*) cdum,c_cut
-        WRITE(otape,*) cdum,c_cut
-!        READ(itape,*) cdum,cs
-!        WRITE(otape,*) cdum,cs
-        READ(itape,*) cdum,c_ref
-        WRITE(otape,*) cdum,c_ref
-        READ(itape,*) cdum,min_wr
-        WRITE(otape,*) cdum,min_wr
-        READ(itape,*) cdum,prof_wr
-        WRITE(otape,*) cdum,prof_wr
 
-       CLOSE(itape)
+     SUBROUTINE SET_IC_TWOSTEP
+      INTEGER :: i
+      REAL :: xi
+      REAL :: x(nx)
 
-!       hx=l_0*pi/REAL(nx)
+      CALL BASE_FLAME_PROFILE(x)
-       hx=l_0*pi/REAL(nx-1)
-       cdum='hx'
-       WRITE(otape,*) cdum,hx
-       cdum='Ta/Tu'
-       WRITE(otape,*) cdum,ac
-       cdum='Tb/Tu'
-       WRITE(otape,*) cdum,bc+1
-       diff=vis/sc
-       cdum='diff'
-       WRITE(otape,*) cdum,diff
 
-!       refwr=pre*1.*exp(-ac/(1.+bc*c_cut))          
+        DO i=1,nx
-!       minf=exp((c_ref-c_cut)*prof_wr)
+         xi = x(i)
 
-       refwr=pre*1.*exp(-ac/(1.+bc*c_ref))          
+         y1(i,1)=(1.-xi)   ! reactant mass fraction
-       l_0=l_0*pi
 
-       ALLOCATE(u(nx,1,1),STAT=ierr)   ; u=0.
+         y1(i,2) = (1./2.) * (lambda1/lambda2) * y1(i,1) * &
-       ALLOCATE(y1(nx,1,1),STAT=ierr)  ; y1=0.
+           exp (-(beta1*(1. - xi))/(1. - hrp*(1. - xi)))
-       ALLOCATE(y2(nx,1,1),STAT=ierr)  ; y2=0.
+
-       ALLOCATE(yf(nx,1,1),STAT=ierr)  ; yf=0.
+         y1(i,3) = xi
+        ENDDO
+     END SUBROUTINE SET_IC_TWOSTEP
+
+     SUBROUTINE init_solver
+
+       INTEGER :: ierr
+
+       CALL READ_INTRO
+
+       CALL init_chemistry
+
+       ALLOCATE(   u(nx),STAT=ierr)    ; u=0.
+       ALLOCATE(  inletbc(nsp),STAT=ierr)  ; inletbc=0.
+       ALLOCATE(  y1(nx,nsp),STAT=ierr)  ; y1=0.
+       ALLOCATE(  y2(nx,nsp),STAT=ierr)  ; y2=0.
+       ALLOCATE(  yf(nx,nsp),STAT=ierr)  ; yf=0.
+       ALLOCATE(yold(nx,nsp),STAT=ierr)  ; yold=0.
        ALLOCATE(uxt(nx),STAT=ierr)     ; uxt=0.
        ALLOCATE(duxt(nx),STAT=ierr)    ; duxt=0.
+       ALLOCATE(dm(nx),STAT=ierr)      ; dm=diff
 
-     END SUBROUTINE READ_INTRO
+       if ( reaction_type == "onestep" ) then
+         inletbc(1) = 1.
+       else if ( reaction_type == "twostep" ) then
+         inletbc(1) = 1.
+         inletbc(2) = 0.
+         inletbc(3) = 0.
+       else
+           WRITE(*,*) 'ERROR, UNDEFINED REACTION TYPE ', reaction_type
+           stop
+       end if
 
-     SUBROUTINE RK4
-       istage=1; CALL substep(y1,y1,y2,yf)
-       istage=2; CALL substep(y1,y2,y1,yf)
-       istage=3; CALL substep(y2,y1,y2,yf)
-       istage=4; CALL substep(y1,y2,y1,yf)
-       istage=5; CALL substep(y2,y1,y2,yf)
-     END SUBROUTINE RK4
-!------------------------------------------------------------------------
-     SUBROUTINE substep(ri,r1,r2,f)
-      REAL, INTENT(INOUT),DIMENSION(:,:,:) :: ri,r1,r2
-      REAL, INTENT(OUT),DIMENSION(:,:,:) :: f
 
-      INTEGER :: i,j,k
+
+
+       CALL ludcmp(nx,100,100,1,0,0)
+
+     END SUBROUTINE init_solver
+
+     SUBROUTINE finalize_solver
+
+       DEALLOCATE( u)
+       DEALLOCATE(inletbc)
+       DEALLOCATE(y1)
+       DEALLOCATE(y2)
+       DEALLOCATE(yf)
+       DEALLOCATE(yold)
+       DEALLOCATE(uxt)
+       DEALLOCATE(duxt)
+       DEALLOCATE(dm)
+
+     END SUBROUTINE finalize_solver
+
+     SUBROUTINE RK4(rhs, calc_diff)
+
+      INTERFACE
+       SUBROUTINE rhs(r1,f)
+        REAL, INTENT(IN),  DIMENSION(:,:) :: r1
+        REAL, INTENT(OUT), DIMENSION(:,:) :: f
+       END SUBROUTINE rhs
+
+       SUBROUTINE calc_diff(r1)
+        REAL, INTENT(IN),  DIMENSION(:,:) :: r1
+       END SUBROUTINE calc_diff
+      END INTERFACE
+
+       CALL substep(1,y1,y1,y2,yf)
+       CALL substep(2,y1,y2,y1,yf)
+       CALL substep(3,y2,y1,y2,yf)
+       CALL substep(4,y1,y2,y1,yf)
+       CALL substep(5,y2,y1,y2,yf)
+
+       CONTAINS
+
+     SUBROUTINE substep(istage, ri,r1,r2,f)
+
+       REAL, DIMENSION(5), PARAMETER :: a=(/   970286171893.d0/4311952581923.,      &
+                                              6584761158862.d0/12103376702013.,     &
+                                              2251764453980.d0/15575788980749.,     &
+                                             26877169314380.d0/34165994151039.,     &
+                                                          0.d0                  /), &
+                                        b=(/  1153189308089.d0/22510343858157.,  &
+                                              1772645290293.d0/4653164025191.,   &
+                                             -1672844663538.d0/4480602732383.,   &
+                                              2114624349019.d0/3568978502595.,   &
+                                              5198255086312.d0/14908931495163.  /)
+
+       INTEGER :: istage
+
+      REAL, INTENT(INOUT),DIMENSION(:,:) :: ri,r1,r2
+      REAL, INTENT(OUT),DIMENSION(:,:) :: f
+
       REAL :: at,bt
 
-       CALL fns(ri,f)
+       CALL calc_diff(ri)
+
+       CALL rhs(ri,f)
+
        IF(istage<5) THEN
         at=a(istage)*dt
         bt=(b(istage)-a(istage))*dt
-        DO k=1,1  ! nz
+        r1=r1+at*f
-        DO j=1,1  ! ny
+        r2=r1+bt*f
-        DO i=1,nx
-        r1(i,j,k)=r1(i,j,k)+at*f(i,j,k)
-        r2(i,j,k)=r1(i,j,k)+bt*f(i,j,k)
-        ENDDO
-        ENDDO
-        ENDDO
        ELSE
         bt=b(istage)*dt
-        DO k=1,1  ! nz
+        r1=r1+bt*f
-        DO j=1,1  ! ny
-        DO i=1,nx
-        r1(i,j,k)=r1(i,j,k)+bt*f(i,j,k)
-        ENDDO
-        ENDDO
-        ENDDO
        ENDIF
      END SUBROUTINE substep
-!------------------------------------------------------------------------
-     SUBROUTINE fns(r1,f)
-      REAL, INTENT(IN),DIMENSION(:,:,:) :: r1
-      REAL, INTENT(OUT),DIMENSION(:,:,:) :: f
-      REAL, DIMENSION(3,nx) :: ux,dux,d2ux
-      INTEGER :: i,j,k
-      REAL :: wrate,Ly,Dy
 
-      DO k=1,1   !nz
+     END SUBROUTINE RK4
-! x-direction
-       DO j=1,1   !ny
-        DO i=1,nx
-         ux(1,i)=u(i,j,k)*r1(i,j,k)                     ! u*Y
-         ux(2,i)=u(i,j,k)                               ! u
-         ux(3,i)=r1(i,j,k)                              ! Y
-        ENDDO
 
-        CALL dfnonp(nx,hx,ux(1:3,:),dux(1:3,:),3,1)
+     SUBROUTINE update_dm(r1)
-        CALL d2fnonp(nx,hx,ux(3:3,:),d2ux(1,:),1,1)
+      REAL, INTENT(IN),DIMENSION(:,:) :: r1
+      INTEGER :: i
 
-        DO i=1,nx
+      DO i=1,nx
-         wrate=pre*ux(3,i)*exp(-ac/(1.+bc*(1.-ux(3,i))))   !wrate
+        dm(i) = diff + d_turb
-!         IF ((1.-ux(3,i)).le.c_cut) THEN
-!          wrate=min_wr
-!          IF ((1.-ux(3,i)).gt.c_ref) wrate= &
-!           ((exp(((1.-ux(3,i))-c_cut)*prof_wr)-minf)/(1.-minf)*(refwr-min_wr))+min_wr
-!         ENDIF
-         IF ((1.-ux(3,i)).le.c_ref) THEN
-           wrate=min_wr
-          IF ((1.-ux(3,i)).gt.c_cut) wrate=((refwr-min_wr)*exp(prof_wr*(1.-ux(3,i)-c_ref))+ &
-          min_wr-refwr*exp(prof_wr*(c_cut-c_ref)))/(1.-exp(prof_wr*(c_cut-c_ref)))
-         ENDIF
-
-! -0.5*( d(u*Y)/dx + u*dY/dx + Y*du/dx ) + D*d2Y/d2x
-         f(i,j,k)=-0.5*( dux(1,i) + ux(2,i)*dux(3,i) + &
-          ux(3,i)*dux(2,i) ) &
-                  + diff*d2ux(1,i) - wrate
-         IF (i.eq.nx) THEN
-          Ly=ux(2,nx)*dux(3,nx)                           ! Ly = u*dYr/dx
-          f(nx,1,1)=-wrate
-         ENDIF
-        ENDDO
-       ENDDO
-!! y-direction
-!       DO i=1,nx
-!        DO j=1,ny
-!         uy(1,j)=v(i,j,k)*r1(i,j,k)                     ! v*Y
-!         uy(2,j)=v(i,j,k)                               ! v
-!         uy(3,j)=r1(i,j,k)                              ! Y
-!        ENDDO
-!        CALL dfyz3(ny,hy,uy(1:3,:),duy)
-!        CALL d2fyz1(ny,hy,uy(3:3,:),d2uy)
-!        DO j=1,ny
-!        f(i,j,k)=f(i,j,k) + &
-!        -0.5*( duy(1,j) + uy(2,j)*duy(3,j) + uy(3,j)*duy(2,j) ) &
-!                   + diff*d2uy(j)
-!        ENDDO
-!       ENDDO
       ENDDO
-!      f(1,1,1)=1.0
+     END SUBROUTINE update_dm
-!      f(nx,1,1)=-0.5*( dux(1,nx) + ux(2,nx)*dux(3,nx) + ux(3,nx)*dux(2,nx) )
-!      f(nx,1,1)=0.0
-!! z-direction
-!      DO j=1,ny
-!       DO i=1,nx
-!        DO k=1,nz
-!         uy(1,k)=w(i,j,k)*r1(i,j,k)                     ! w*Y
-!         uy(2,k)=w(i,j,k)                               ! w
-!         uy(3,k)=r1(i,j,k)                              ! Y
-!        ENDDO
-!        CALL dfyz3(ny,hy,uy(1:3,:),duy)
-!        CALL d2fyz1(ny,hy,uy(3:3,:),d2uy)
-!        DO k=1,nz
-!        f(i,j,k)=f(i,j,k) + &
-!        -0.5*( duy(1,k) + uy(2,k)*duy(3,k) + uy(3,k)*duy(2,k) ) &
-!                        + diff*d2uy(k)
-!        ENDDO
-!       ENDDO
-!      ENDDO
 
-!      DO j=1,ny
+     SUBROUTINE update_dm_cd(r1)
-!       DO k=1,nz
+      REAL, INTENT(IN),DIMENSION(:,:) :: r1
-!         f(1,j,k)=0.0
+      REAL :: conv_rxn_boundary
-!         f(nx,j,k)=-0.5*(dux(1,nx)+ux(2,nx)*dux(3,nx)+ux(3,nx)*&
+      REAL :: beta, tutb
-!           dux(2,nx))
+      INTEGER :: i
-!       ENDDO
-!      ENDDO
 
-!    Boundary conditionS
+      tutb = 1. / (bc+1.)
-      f(1,1,1)=0.
 
-      Dy=Ly
+      ! \beta = Ta/Tb * (Tb-Tu)/Tb
-      f(nx,1,1)=f(nx,1,1)-Dy
+      beta = ac * tutb * (1. - tutb)
+
+      conv_rxn_boundary = 1./beta
+
+      DO i=1,nx
+        IF (r1(i,fctrl_species) > conv_rxn_boundary) THEN
+          dm(i) = diff + d_turb
+        ELSE
+          dm(i) = diff
+        END IF
+      ENDDO
+     END SUBROUTINE update_dm_cd
+
+     SUBROUTINE update_dm_sutherland(r1)
+      REAL, INTENT(IN),DIMENSION(:,:) :: r1
+      INTEGER :: i
+
+      DO i=1,nx
+        dm(i) = diff * diffusivity_sutherland(r1(i,nsp)) + d_turb
+      ENDDO
+     END SUBROUTINE update_dm_sutherland
+
+     SUBROUTINE fonestep(r1,f)
+      REAL, INTENT(IN),DIMENSION(:,:) :: r1
+      REAL, INTENT(OUT),DIMENSION(:,:) :: f
+      REAL, DIMENSION(nsp,nx) :: ux, dux, d2ux
+      REAL, DIMENSION(nx) :: dxdm
+      INTEGER :: i
+      REAL :: wrate,Ly,Dy,Lt,Dt
+
+! x-direction
+      DO i=1,nx
+         ux(1,i)=r1(i,1)                              ! Y
+         ux(2,i)=r1(i,2)                              ! T
+      ENDDO
+
+      CALL dfnonp(nx,hx,ux(:,:),dux(:,:),nsp,1)
+
+      CALL dfnonp(nx,hx,dm,dxdm,1,1)
+
+      CALL d2fnonp(nx,hx,ux(:,:),d2ux(:,:),nsp,1)
+
+      DO i=1,nx
+        wrate=rate_1step(ux(1,i), ux(2,i))
+
+        ! - u*dY/dx + D*d2Y/d2x
+        f(i,1) = - ( u(i)*dux(1,i) ) + dm(i) * d2ux(1,i) + dxdm(i) * dux(1,i) - wrate
+
+        ! - u*dY/dx + D*d2Y/d2x
+        f(i,2) = - ( u(i)*dux(2,i) ) + dm(i) * d2ux(2,i) + dxdm(i) * dux(2,i) + wrate
+
+        ! Boundary conditions
+        IF (i.eq.nx) THEN
+            f(nx,1) = -wrate - u(nx)*dux(1,nx)
+            f(nx,2) = wrate - u(nx)*dux(2,nx)
+        ENDIF
+      ENDDO
+
+      ! Boundary conditions
+      f(1,1)=0.
+      f(1,2)=0.
+
+     END SUBROUTINE fonestep
+
+     SUBROUTINE fns(r1,f)
+      REAL, INTENT(IN),DIMENSION(:,:) :: r1
+      REAL, INTENT(OUT),DIMENSION(:,:) :: f
+      REAL, DIMENSION(3,nx) :: ux, dux, d2ux
+      INTEGER :: i
+      REAL :: Ly,Dy,Lt,Dt
+      REAL :: wrate1, wrate2
+
+! x-direction
+      DO i=1,nx
+         ux(1,i)=r1(i,1)                              ! Y
+         ux(2,i)=r1(i,2)                              ! Y
+         ux(3,i)=r1(i,3)                              ! T
+      ENDDO
+
+      CALL dfnonp(nx,hx,ux(:,:),dux(:,:),3,1)
+
+      CALL d2fnonp(nx,hx,ux(:,:),d2ux(:,:),3,1)
+
+      DO i=1,nx
+        wrate1=rate1_2step(ux(1,i), ux(2,i), ux(3,i))
+        wrate2=rate2_2step(ux(2,i), ux(3,i))
+
+        ! - u*dY/dx + D*d2Y/d2x
+        f(i,1) = - ( u(i)*dux(1,i) ) + (diff/le_a) * d2ux(1,i) - wrate1
+
+        ! - u*dY/dx + D*d2Y/d2x
+        f(i,2) = - ( u(i)*dux(2,i) ) + (diff/le_x) * d2ux(2,i) + wrate1 - 2.*wrate2
+
+        ! -  u*dT/dx + D*d2T/d2x
+        f(i,3) = - ( u(i)*dux(3,i) ) + (diff) * d2ux(3,i) + 2.*wrate2
+
+        ! Boundary conditions
+        IF (i.eq.nx) THEN
+            f(nx,1) = -wrate1 - u(nx)*dux(1,nx)
+            f(nx,2) = wrate1 - 2.*wrate2 - u(nx)*dux(2,nx)
+            f(nx,3) = 2.*wrate2 - u(nx)*dux(3,nx)
+        ENDIF
+      ENDDO
+
+      ! Boundary conditions
+      f(1,1)=0.
+      f(1,2)=0.
+      f(1,3)=0.
 
      END SUBROUTINE fns
-!------------------------------------------------------------------------
+
+     REAL FUNCTION diffusivity_sutherland(c)
+         REAL, INTENT(IN) :: c
+         REAL :: theta, As, Ts, T0, T1
+         REAL, PARAMETER :: rvis=5.0d0
+         T0 = 1.0
+         T1 = (1.0 + bc)
+
+         Ts = (rvis*T1 - (T1**(3./2.))) / (T1**(3./2.) - rvis)
+         As = (T0 + Ts)
+
+         theta = (1.0 + bc * c)
+
+         diffusivity_sutherland = As * sqrt(theta) / (1. + Ts/theta)
+     END FUNCTION diffusivity_sutherland
+
     END MODULE ysolve