From fea70aa4cb375da81aafd5c61b46f613ffc5c367 Mon Sep 17 00:00:00 2001
From: ignis <ignis@postech.ac.kr>
Date: Sun, 1 Sep 2019 14:36:17 +0900
Subject: [PATCH] refactoring, put rk substeps into rk4

---
 m_fdm_calc.f90 | 107 ++++++++++++++++++++++---------------------------
 1 file changed, 47 insertions(+), 60 deletions(-)

diff --git a/m_fdm_calc.f90 b/m_fdm_calc.f90
index 53aa4ad..e428f4e 100644
--- a/m_fdm_calc.f90
+++ b/m_fdm_calc.f90
@@ -693,6 +693,27 @@ module m_fdm_calc
       END SUBROUTINE fns
 
 
+      SUBROUTINE solve(xx,yy,zz,uu_,vv_,ww_,yy1,yy2,yyf)
+
+       IMPLICIT NONE
+
+       integer :: i,j,k,xx,yy,zz
+       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)
+
+       ! 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
+
+
       subroutine EE1(xx,yy,zz,uu_,vv_,ww_,yy1,rhs1)
 
       implicit none
@@ -774,48 +795,7 @@ 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 rotarysubstep(1, xx, yy, zz, uu_, vv_, ww_, yy1, yy1, yy2)
-      CALL rotarysubstep(2, xx, yy, zz, uu_, vv_, ww_, yy1, yy2, yyf)
-      CALL rotarysubstep(3, xx, yy, zz, uu_, vv_, ww_, yy2, yyf, yy1)
-      CALL rotarysubstep(4, xx, yy, zz, uu_, vv_, ww_, yyf, yy1, yy2)
-      CALL lastsubstep     (xx, yy, zz, uu_, vv_, ww_, yy1, yy2, yyf)
-
-      return
-      END SUBROUTINE RK4
-
-      SUBROUTINE solve(xx,yy,zz,uu_,vv_,ww_,yy1,yy2,yyf)
-
-       IMPLICIT NONE
-
-       integer :: i,j,k,xx,yy,zz
-       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)
-
-       ! 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
-
-      subroutine rotarysubstep(istage,xx,yy,zz,uu_,vv_,ww_,r0,r1,r2)
-
-      implicit none
-
-      integer :: xx,yy,zz,istage
-      real*8 :: uu_(xx,yy,zz),vv_(xx,yy,zz),ww_(xx,yy,zz)
-      real*8 :: r0(xx,yy,zz,neq),r1(xx,yy,zz,neq),r2(xx,yy,zz,neq)
-
       real*8 :: a(5),b(5)
-      real*8 :: at, bt
-
-      integer :: i,j,k,l
 
       a(1)= 970286171893./4311952581923.
       a(2)= 6584761158862./12103376702013.
@@ -829,6 +809,29 @@ module m_fdm_calc
       b(4)= 2114624349019./3568978502595.
       b(5)= 5198255086312./14908931495163.
 
+
+      CALL rotarysubstep(1, xx, yy, zz, uu_, vv_, ww_, yy1, yy1, yy2)
+      CALL rotarysubstep(2, xx, yy, zz, uu_, vv_, ww_, yy1, yy2, yyf)
+      CALL rotarysubstep(3, xx, yy, zz, uu_, vv_, ww_, yy2, yyf, yy1)
+      CALL rotarysubstep(4, xx, yy, zz, uu_, vv_, ww_, yyf, yy1, yy2)
+      CALL lastsubstep     (xx, yy, zz, uu_, vv_, ww_, yy1, yy2, yyf)
+
+      return
+
+      contains
+
+      subroutine rotarysubstep(istage,xx,yy,zz,uu_,vv_,ww_,r0,r1,r2)
+
+      implicit none
+
+      integer :: xx,yy,zz,istage
+      real*8 :: uu_(xx,yy,zz),vv_(xx,yy,zz),ww_(xx,yy,zz)
+      real*8 :: r0(xx,yy,zz,neq),r1(xx,yy,zz,neq),r2(xx,yy,zz,neq)
+
+      real*8 :: at, bt
+
+      integer :: i,j,k,l
+
        CALL fns(r0,r2,xx,yy,zz,uu_,vv_,ww_)
 
         at=a(istage)*fdmdt
@@ -860,15 +863,11 @@ module m_fdm_calc
       real*8 :: r0(xx,yy,zz,neq), r1(xx,yy,zz,neq), r2(xx,yy,zz,neq)
       real*8 :: uu_(xx,yy,zz), vv_(xx,yy,zz), ww_(xx,yy,zz)
 
-      real*8 :: b
-
       integer :: i, j, k
 
-      b = 5198255086312./14908931495163.
-
        CALL fns(r0,r2,xx,yy,zz,uu_,vv_,ww_)
 
-        bt = b*fdmdt
+        bt = b(5)*fdmdt
 
         r0 = r1 + bt*r2
 
@@ -901,23 +900,9 @@ module m_fdm_calc
       integer :: i,j,k,xx,yy,zz,istage
       real*8 :: at,bt , wrate , yr
       real*8 :: ri(xx,yy,zz,neq),r1(xx,yy,zz,neq),r2(xx,yy,zz,neq),f(xx,yy,zz,neq)
-      real*8 :: a(5),b(5)
       real*8 :: uu_(xx,yy,zz),vv_(xx,yy,zz),ww_(xx,yy,zz)
       integer :: nfinal, iscr, mspec, mpict, msave, nmindt, nv
 
-
-      a(1)= 970286171893./4311952581923.
-      a(2)= 6584761158862./12103376702013.
-      a(3)= 2251764453980./15575788980749.
-      a(4)= 26877169314380./34165994151039.
-      a(5)=0.
-
-      b(1)= 1153189308089./22510343858157.
-      b(2)= 1772645290293./4653164025191.
-      b(3)= -1672844663538./4480602732383.
-      b(4)= 2114624349019./3568978502595.
-      b(5)= 5198255086312./14908931495163.
-
        CALL fns(ri,f,xx,yy,zz,uu_,vv_,ww_)
 
       IF(istage<5) THEN
@@ -957,6 +942,8 @@ module m_fdm_calc
       return
       END SUBROUTINE substep
 
+      END SUBROUTINE RK4
+
   subroutine tp2 (a, b, n1, n2)
     ! a = transpose(b)