MODULE ysolve
     USE m_compact
     USE m_parameters
     USE m_chemistry

     IMPLICIT NONE

     REAL, DIMENSION(:), ALLOCATABLE :: u
     REAL, DIMENSION(:), ALLOCATABLE :: inletbc
     REAL, DIMENSION(:,:), ALLOCATABLE :: y1,y2,yf,yold
     REAL, DIMENSION(:), ALLOCATABLE :: uxt,duxt
     REAL, DIMENSION(:), ALLOCATABLE :: dm

    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
       IF (d_mode == 0) THEN
         CALL solve_with_diffusivity (update_dm)
       ELSE IF (d_mode == 1) THEN
         CALL solve_with_diffusivity (update_dm)
       ELSE
         WRITE(*,*)'ERROR, UNSUPPORTED DIFFUSIVITY ', d_mode, '. STOPPING'
         STOP
       END IF
     END SUBROUTINE solve

     SUBROUTINE solve_with_diffusivity (calc_diff)
      INTERFACE
       SUBROUTINE calc_diff
       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
       END SUBROUTINE calc_diff
      END INTERFACE

      INTEGER :: i,j,k,savenum
      REAL :: pflame,pflold,delf=0.
      REAL :: residue = 0.

      CALL SET_IC

      pflame=0.
      DO i=1,nx
       pflame=pflame+y1(i,fctrl_species)*hx
      ENDDO


      t_now=0.
      t_uf=0.

      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)

       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
      ENDDO

      ! CALL write_sd
      CALL write_pre
      CALL save_final_field

      WRITE(*,*)
      WRITE(*,*) 'Fin.'
      WRITE(*,*)

     END SUBROUTINE solve_

     SUBROUTINE write_sd
      REAL :: c,yr,theta,wrate,dely,sdr,sdd,sd,uu,onelw,dd
      INTEGER :: i,j,k,nd
      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)                               ! Yr
       ux(2,i)=y1(i,2)                               ! T
       IF (ux(1,i).gt.1.) ux(1,i)=1.                                  
      ENDDO

      nd=2
      CALL dfnonp(nx,hx,ux(:,:),dux(:,:),nd,1)
      CALL d2fnonp(nx,hx,ux(:,:),d2ux(:,:),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*(theta)))   !wrate

!       IF (c.le.c_cut) THEN
!        wrate=min_wr
!        IF (c.gt.c_ref) wrate= &
!         ((exp((c-c_cut)*prof_wr)-minf)/(1.-minf)*(refwr-min_wr))+min_wr
!       ENDIF

       IF (theta.le.c_ref) THEN
         wrate=min_wr
        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=-dm(i)*d2ux(1,i)/dely
       IF (dely.eq.0.) THEN
        sdr=0.; sdd=0.
       ENDIF
       sd=sdr+sdd
       uu=u(1)
       onelw=(-dux(1,i))/c
       dd=-dm(i)*d2ux(1,i)
       if (c.eq.0.) onelw=0.

       sav(1,i)=sav(1,i)+yr
       sav(2,i)=sav(2,i)+c
       sav(3,i)=sav(3,i)+uu
       sav(4,i)=sav(4,i)+wrate
       sav(5,i)=sav(5,i)+dely
       sav(6,i)=sav(6,i)+sdr
       sav(7,i)=sav(7,i)+sdd
       sav(8,i)=sav(8,i)+sd
       sav(9,i)=sav(9,i)+onelw
       sav(10,i)=sav(10,i)+dd
      ENDDO

      DO i=1,nx
       WRITE(500,'(37e30.20)') (i-1)*hx,sav(1:10,i)
      ENDDO
     END SUBROUTINE write_sd

     SUBROUTINE save_final_field
      INTEGER :: i,j,k

      OPEN (305,FILE='sfield.bin',form='unformatted',status='unknown')
      DO i=1,nx
       WRITE (305) y1(i,1)
      ENDDO
      CLOSE (305)

      OPEN (305,FILE='sfield.dat')
      DO i=1,nx
       WRITE (305,*) (/ i*hx, y1(i,:) /)
      ENDDO
      CLOSE (305)

     END SUBROUTINE save_final_field

     SUBROUTINE write_pre
      REAL :: theta,yr,c,dy,maxdy=0.,del_f
      REAL :: S_L=0.,wrate,wrate1,wrate2
      REAL :: ya,yx
      INTEGER :: i
      REAL, DIMENSION(1,nx) :: ux, dux

      if ( reaction_type == "onestep" ) then

      DO i=1,nx
       yr=y1(i,1)
       theta=y1(i,2)

       ux(1,i)=theta

        wrate=rate_1step(yr, theta)

       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)

      DO i=(nx/10),nx-(nx/10)
       dy=dux(1,i)
       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)
      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

      DO i = 1, nsp
       y1(1,i) = inletbc(i)
      END DO
     END SUBROUTINE SET_BC

     SUBROUTINE BASE_FLAME_PROFILE(x)
      INTEGER :: i, ifl, si
      REAL :: xi
      REAL :: x(nx)
      REAL :: max_ysum

      ! initial flame thickness (0.2 pi n grids)
      ifl=INT(2.0*pi*0.1/hx)

      ! 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
        x(i)=0.+ctmp
       ELSE IF(i> nx-(si-ifl/2)) THEN
        x(i)=1.
       ELSE
        x(i)=0.5+REAL(i-nx+si)/REAL(ifl)
       ENDIF
      END DO

      ! initialize species field
      ! max_ysum = maxval(y1(:,1) + y1(:,2))
      ! 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))

     END SUBROUTINE BASE_FLAME_PROFILE

     SUBROUTINE SET_IC_ONESTEP
      INTEGER :: i
      REAL :: xi
      REAL :: x(nx)

      CALL BASE_FLAME_PROFILE(x)

        DO i=1,nx
         xi = x(i)

         y1(i,1)=(1.-xi)   ! reactant mass fraction
         y1(i,2)=(xi)   ! reactant mass fraction
        ENDDO
     END SUBROUTINE SET_IC_ONESTEP

     SUBROUTINE SET_IC_TWOSTEP
      INTEGER :: i
      REAL :: xi
      REAL :: x(nx)

      CALL BASE_FLAME_PROFILE(x)

        DO i=1,nx
         xi = x(i)

         y1(i,1)=(1.-xi)   ! reactant mass fraction

         y1(i,2) = (1./2.) * (lambda1/lambda2) * y1(i,1) * &
           exp (-(beta1*(1. - xi))/(1. - hrp*(1. - xi)))

         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

       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




       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)

      INTERFACE
       SUBROUTINE rhs(r1,f)
        REAL, INTENT(IN),  DIMENSION(:,:) :: r1
        REAL, INTENT(OUT), DIMENSION(:,:) :: f
       END SUBROUTINE rhs
      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

      INTEGER :: i,j,k
      REAL :: at,bt

       CALL update_dm

       CALL rhs(ri,f)

       IF(istage<5) THEN
        at=a(istage)*dt
        bt=(b(istage)-a(istage))*dt
        r1=r1+at*f
        r2=r1+bt*f
       ELSE
        bt=b(istage)*dt
        r1=r1+bt*f
       ENDIF
     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

! 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) = -wrate1 - u(nx)*dux(1,nx)
            f(nx,2) = wrate1 - 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,j,k
      REAL :: wrate,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
        wrate=rate_1step(ux(1,i), ux(3,i))
        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

    END MODULE ysolve