MODULE ysolve
     USE Compact
     IMPLICIT NONE
     PRIVATE
     REAL, PARAMETER :: pi=3.14159265358979323846
     REAL :: hx,dt,vis,sc,diff,pre,ac,bc,tf,t_now,t_uf,dt_uf
     REAL :: rvis,cstar,sigw
     REAL :: c_cut,c_ref,refwr,minf,tar_lo,u0,ctmp,l_0,lo_flm=0.
     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. /)

     LOGICAL :: read_itape, read_stdin
     CHARACTER(100) :: itape_name
     PUBLIC :: solve, parse
    CONTAINS
!------------------------------------------------------------------------
     SUBROUTINE parse
        CHARACTER(100) :: num1char

        !First, make sure the right number of inputs have been provided
        IF(COMMAND_ARGUMENT_COUNT().EQ.0)THEN
          ! read from itape
          read_itape = .true.
        ELSE IF(COMMAND_ARGUMENT_COUNT().EQ.1)THEN
          CALL GET_COMMAND_ARGUMENT(1,num1char)   !first, read in the two values
        ELSE
          WRITE(*,*)'ERROR, TOO MANY COMMAND-LINE ARGUMENTS(MORE THAN ONE). STOPPING'
          STOP
        ENDIF

        read_stdin = (num1char=="-")

        IF(read_itape)THEN
          itape_name = "itape"
          WRITE(*,*) "Read from itape"
        ELSE IF(read_stdin)THEN
          WRITE(*,*) "Read from STDIN"
        ELSE 
          itape_name = num1char
          WRITE(*,*) "Read from " // num1char
        END IF

     END SUBROUTINE parse
!------------------------------------------------------------------------
     SUBROUTINE solve
      INTEGER :: i,j,k,savenum
      REAL :: pflame,pflold,delf=0.

      CALL READ_INTRO

      CALL ludcmp(nx,5,5,1,0,0)

      CALL SET_IC

      t_now=0.; t_uf=0.

      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)
!        WRITE(*,*)
       ENDIF
       ncyc=ncyc+1
       t_uf=t_uf+dt
       t_now=t_now+dt
       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,wrate,dely,sdr,sdd,sd,uu,onelw,dd
      INTEGER :: i,j,k,nd
      REAL, DIMENSION(1,nx) :: ux,dux,d2ux
      REAL, DIMENSION(1,nx) :: duy
      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
       IF (ux(1,i).gt.1.) ux(1,i)=1.                                  
      ENDDO
          
      dm = (/ (diff * diffusivity_sutherland(1.-y1(i,1,1)), i=1,nx) /)

      nd=1
      CALL dfnonp(nx,hx,ux(1,:),dux(1,:),nd,1)  
      CALL dfnonp(nx,hx,dm,duy(1,:),nd,1)
      nd=1 
      CALL d2fnonp(nx,hx,ux(1,:),d2ux(1,:),nd,1) 


      DO i=1,nx
       yr=ux(1,i)
       c=1.-yr
       IF (c.lt.0.) c=0.

       wrate=pre*yr*exp(-ac/(1.+bc*(1.-yr)))   !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 (c.le.c_ref) THEN
         wrate=min_wr
        IF (c.gt.c_cut) wrate=((refwr-min_wr)*exp(prof_wr*(c-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)+duy(1,i)*dux(1,i))/dely
       IF (dely.eq.0.) THEN
        sdr=0.; sdd=0.
       ENDIF
       sd=sdr+sdd
       uu=u(1,1,1)
       onelw=(-dux(1,i))/c
       dd=-(dm(i)*d2ux(1,i)+duy(1,i)*dux(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,1)
      ENDDO
      CLOSE (305)

      OPEN (305,FILE='sfield.dat')
      DO i=1,nx
       WRITE (305,'(e30.20)') y1(i,1,1)
      ENDDO
      CLOSE (305)

     END SUBROUTINE save_final_field

     SUBROUTINE write_pre
      REAL :: yr,c,dy,maxdy=0.,del_f
      REAL :: S_L=0.,wrate
      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 (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 (c.le.c_ref) THEN
         wrate=min_wr
        IF (c.gt.c_cut) wrate=((refwr-min_wr)*exp(prof_wr*(c-c_ref))+ &
        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
      
      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,1,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
      
      y1(1,1,1)=1.-ctmp
     END SUBROUTINE SET_BC

     SUBROUTINE SET_IC
      INTEGER :: i, ifl, si
      REAL :: xi

      ifl=INT(2.0*pi*0.1/hx)

      u=u0; si=INT(nx*(1.-tar_lo))

      DO i=1,nx
       IF(i< nx-(si+ifl/2)) THEN
        xi=0.+ctmp
       ELSE IF(i> nx-(si-ifl/2)) THEN
        xi=1.
       ELSE
        xi=0.5+REAL(i-nx+si)/REAL(ifl)
       ENDIF

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

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

     END SUBROUTINE SET_IC

     SUBROUTINE READ_INTRO
       CHARACTER(LEN=8) :: cdum
       INTEGER :: itape=300,otape=301,ierr
       INTEGER :: i

       IF (read_stdin) THEN
         itape=5
       ELSE
         OPEN(itape,FILE=itape_name)
       END IF

       OPEN(otape,FILE='otape')
        READ(itape,*) cdum,nx
        WRITE(otape,*) cdum,nx
        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
        READ(itape,*) cdum,rvis
        WRITE(otape,*) cdum,rvis

       IF (.not.read_stdin) THEN
         CLOSE(itape)
       END IF

!       hx=l_0*pi/REAL(nx)
       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))          
!       minf=exp((c_ref-c_cut)*prof_wr)

       refwr=pre*1.*exp(-ac/(1.+bc*c_ref))          
       l_0=l_0*pi

       ALLOCATE(u(nx,1,1),STAT=ierr)   ; u=0.
       ALLOCATE(y1(nx,1,1),STAT=ierr)  ; y1=0.
       ALLOCATE(y2(nx,1,1),STAT=ierr)  ; y2=0.
       ALLOCATE(yf(nx,1,1),STAT=ierr)  ; yf=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

     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
      REAL :: at,bt

       CALL fns(ri,f)
       IF(istage<5) THEN
        at=a(istage)*dt
        bt=(b(istage)-a(istage))*dt
        DO k=1,1  ! nz
        DO j=1,1  ! ny
        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
        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
      REAL, DIMENSION(1,nx) :: duy
      INTEGER :: i,j,k
      REAL :: wrate,Ly,Dy

      DO k=1,1   !nz
! 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

      dm = (/ (diff * diffusivity_sutherland(1.-r1(i,1,1)), i=1,nx) /)

        CALL dfnonp(nx,hx,ux(1:3,:),dux(1:3,:),3,1)
        CALL dfnonp(nx,hx,dm,duy(1,:),1,1)
        CALL d2fnonp(nx,hx,ux(3:3,:),d2ux(1,:),1,1)

        DO i=1,nx
         wrate=pre*ux(3,i)*exp(-ac/(1.+bc*(1.-ux(3,i))))   !wrate
!         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) ) &
                  + dm(i)*d2ux(1,i) + duy(1,i)*dux(3,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
      ENDDO

!    Boundary conditionS
      f(1,1,1)=0.

      Dy=Ly
      f(nx,1,1)=f(nx,1,1)-Dy

     END SUBROUTINE fns
!------------------------------------------------------------------------
     REAL FUNCTION sigmoid(x)
         REAL, INTENT(IN) :: x
         sigmoid = 1. / (1 + exp(-x))
     END FUNCTION sigmoid
!------------------------------------------------------------------------
     REAL FUNCTION diffusivity_sigmoid(c)
         REAL, INTENT(IN) :: c
         diffusivity_sigmoid = 1. + (rvis - 1.) * sigmoid(sigw*(c-cstar))
     END FUNCTION diffusivity_sigmoid
!------------------------------------------------------------------------
     REAL FUNCTION diffusivity_T076(c)
         REAL, INTENT(IN) :: c
         diffusivity_T076 = ((1.0 + bc * c) ** 0.76)
     END FUNCTION diffusivity_T076
!------------------------------------------------------------------------
     REAL FUNCTION diffusivity_sutherland(c)
         REAL, INTENT(IN) :: c
         REAL :: theta, As, Ts, T0, T1
         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