From 078e7d82e75816ea8b0661e27773278d61022b89 Mon Sep 17 00:00:00 2001
From: ignis <ignis@postech.ac.kr>
Date: Wed, 26 Feb 2020 11:49:02 +0900
Subject: [PATCH] stubs for dm modification

---
 code/ysolve.f90 | 33 +++++++++++++++++++++++++++------
 1 file changed, 27 insertions(+), 6 deletions(-)

diff --git a/code/ysolve.f90 b/code/ysolve.f90
index f2e2e70..ede43c3 100644
--- a/code/ysolve.f90
+++ b/code/ysolve.f90
@@ -89,6 +89,16 @@
 !------------------------------------------------------------------------
 
      SUBROUTINE solve
+       IF (reaction_type == "onestep") THEN
+         CALL solve_
+       ELSE IF (reaction_type == "twostep") THEN
+         CALL solve_
+       ELSE
+         CALL solve_
+       END IF
+     END SUBROUTINE solve
+
+     SUBROUTINE solve_
 
       INTEGER :: i,j,k,savenum
       REAL :: pflame,pflold,delf=0.
@@ -161,7 +171,7 @@
       WRITE(*,*) 'Fin.'
       WRITE(*,*)
 
-     END SUBROUTINE solve
+     END SUBROUTINE solve_
 
      SUBROUTINE write_sd
       REAL :: c,yr,theta,wrate,dely,sdr,sdd,sd,uu,onelw,dd
@@ -475,6 +485,8 @@
       INTEGER :: i,j,k
       REAL :: at,bt
 
+       CALL update_dm
+
        CALL rhs(ri,f)
 
        IF(istage<5) THEN
@@ -489,11 +501,20 @@
      END SUBROUTINE substep
 
      END SUBROUTINE RK4
-!------------------------------------------------------------------------
+
+     SUBROUTINE update_dm
+      INTEGER :: i
+
+      DO i=1,nx
+        dm(i) = diff + d_turb
+      ENDDO
+     END SUBROUTINE update_dm
+
      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,j,k
       REAL :: wrate,Ly,Dy,Lt,Dt
       REAL :: wrate1, wrate2
@@ -506,16 +527,18 @@
 
       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) ) + (diff) * d2ux(1,i) - wrate
+        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) ) + (diff) * d2ux(2,i) + wrate
+        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
@@ -578,6 +601,4 @@
 
      END SUBROUTINE fns
 
-!------------------------------------------------------------------------
-
     END MODULE ysolve