From 1725f1eaf1b6d4e2f51536bb3967217376b4d688 Mon Sep 17 00:00:00 2001
From: ignis <ignis@postech.ac.kr>
Date: Sat, 31 Aug 2019 07:33:04 +0900
Subject: [PATCH] minor mods ysolve

---
 code/ysolve.f90 | 60 ++++++++++++++++++++++++++++---------------------
 1 file changed, 35 insertions(+), 25 deletions(-)

diff --git a/code/ysolve.f90 b/code/ysolve.f90
index 82be8ef..b123150 100644
--- a/code/ysolve.f90
+++ b/code/ysolve.f90
@@ -65,14 +65,15 @@
       INTEGER :: i,j,k,savenum
       REAL :: pflame,pflold,delf=0.
 
-      CALL ludcmp(nx,5,5,1,0,0)
+      CALL ludcmp(nx,100,100,1,0,0)
 
       CALL SET_IC
 
-      t_now=0.; t_uf=0.
+      t_now=0.
+      t_uf=0.
 
-      DO 
-       CALL SET_BC  
+      DO
+       CALL SET_BC
        CALL RK4(fns)
        IF(t_now.ge.tf) EXIT
        IF(t_uf.ge.dt_uf) THEN
@@ -249,7 +250,6 @@
      END SUBROUTINE write_pre
 
      SUBROUTINE SET_BC
-      
       y1(1,1)=1.-ctmp
      END SUBROUTINE SET_BC
 
@@ -257,10 +257,16 @@
       INTEGER :: i, ifl, si
       REAL :: xi
 
+      ! 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
@@ -274,8 +280,10 @@
        y1(i,1)=(1.-xi)   ! reactant mass fraction
       ENDDO
 
+      ! initialize temperature field
       y1(:,2) = 1.0 - y1(:,1)
 
+      ! calculate flame position
       pflame=0.
       DO i=1,nx
        pflame=pflame+y1(i,1)*hx
@@ -320,28 +328,29 @@
        END SUBROUTINE rhs
       END INTERFACE
 
-       REAL, DIMENSION(5), PARAMETER :: a=(/   970286171893./4311952581923.,   &
-                                              6584761158862./12103376702013.,  &
-                                              2251764453980./15575788980749.,  &
-                                             26877169314380./34165994151039.,  &
-                                                          0.                  /), &
-                                        b=(/  1153189308089./22510343858157.,  &
-                                              1772645290293./4653164025191.,   &
-                                             -1672844663538./4480602732383.,   &
-                                              2114624349019./3568978502595.,   &
-                                              5198255086312./14908931495163.  /)
-
-       INTEGER :: istage
-
-       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)
+       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(ri,r1,r2,f)
+     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
 
@@ -349,6 +358,7 @@
       REAL :: at,bt
 
        CALL rhs(ri,f)
+
        IF(istage<5) THEN
         at=a(istage)*dt
         bt=(b(istage)-a(istage))*dt