module global_props
!This module is used to declare the fluid property storage variables for the user-defined
!materials
implicit none
save
  double precision,allocatable:: fprop(:,:),db(:,:)
  integer::fl_bounds(25), fl_flag, fl_ct, LuFlEr
  double precision Gjsav,ttime
end module

!################################################################################

subroutine readFluidPropFile(LU_FL)
use global_props
!This subroutine opens and reads the fluid property file specified by 
!..the user, and allocates space for the property array in the 
!..global_props module
implicit none

integer,intent(in):: LU_FL
integer::ios,ntot,nufl,i,j,m,actb
character::tc,checkname*300

!If user-defined fluid properties have not already been loaded by another unit of this type,
!..do so now.  
if(.not.allocated(fprop)) then
    inquire(unit=LU_FL,name=checkname) !Find the fluid property file name
    open(unit=LU_fL,file=trim(checkname),iostat=ios,status="OLD",position="REWIND")  !Open the file

    if(ios.ne.0) then !there was a problem opening the file
        allocate(fprop(1,1))
        close(LU_FL)
        goto 100  !If there was a problem opening the file, skip out and continue
    endif 
    ntot=0 ; nufl=0 ; fl_bounds=0 ; ios=0;
    !First check to see how many fluids are included, and how many total lines there are            
    do
        read(LU_FL,fmt='(A1)',iostat=ios) tc
        if(ios.lt.0) exit  !negative ios means end of file
        if(tc.ne.'&') then
            ntot=ntot+1 !increment the number of entries
        else 
            nufl=nufl+1 !increment the number of user fluids
            fl_bounds(nufl)=ntot+1 !Specify the lower bounds
        endif
    enddo
    ! Always set the upper bound (+1)
    fl_bounds(nufl+1)=ntot+1
    
    if(nufl.le.0) then !There are currently no user fluids, so skip the rest
        allocate(fprop(1,1))
        close(LU_FL)
        goto 100  
    endif
    
    !Allocate and initiate fprop array
    allocate(fprop(7,ntot)) ; fprop=0. ;
    
    !Go back and read in the values, keep track of where each fluid's entries end
    rewind(LU_FL) ; ios=0;
    m=1; actb=fl_bounds(m); i=1;  !initialize values
    do
        if(i.ne.actb) then  !If the index i is NOT equal to one of the stored indices known to be a new fluid marker
            read(LU_FL,fmt=*,iostat=ios) (fprop(j,i),j=1,7)  !Read in the fluid properties
            i=i+1
            if(i.gt.ntot) exit !We've read through the file
        else  !If the current line is known to be a new fluid marker
            read(LU_FL,fmt='(A1)') tc
            m=m+1  !Increment to the next marker number
            actb=fl_bounds(m) !Set the active boundary to the next marker
        endif
        
        if(ios.lt.0) exit !end of file
    enddo


    close(LU_FL)
endif            
100 continue      


end subroutine

!**************************************************************************************************
!**  The following functions are used to interpolate user-defined material property tables.      **
!**  Code is adapted from Numerical Recipes, 3rd Edition.                                        **
!**************************************************************************************************

subroutine interp(x,n,xx,yy,jsav,y)

    !********************************************************************************************
    !* Given a value x, return an interpolated value y, using data xx and yy of size n. Save   **
    !* previous table location as jsav, and return it to the subroutine when calling.          **
    !********************************************************************************************
    integer,intent(in)::n
    integer,intent(inout)::jsav
    integer  cor, j, dj
    double precision x,xx(n),yy(n)
    double precision,intent(out)::y

    !Estimated limit for using hunt vs locate methods
    dj=min(1,nint(n**.25))  

    !Decide which interpolation method to use
    if (cor.eq.1) then  
        call hunt(x,n,xx,yy,dj,cor,j,jsav)  !We anticipate a location close to the previous
    else  
        call locate(x,n,xx,yy,dj,cor,j,jsav)  !The anticipated location is not near the previous
    endif

    !Calculate the interpolated value
    y = yy(j)+((x-xx(j))/(xx(j+1)-xx(j)))*(yy(j+1)-yy(j))

end subroutine


!--------------------------------------------------------------------------------------------


subroutine locate(x,n,xx,yy,dj,cor,j,jsav)

    !********************************************************************************************
    !* Given a value x, return a value j such that x is (insofar as possible) centered in the   *
    !* subrange xx[j..j+mm-1], where xx is the stored pointer.  The values in xx must be mono-  *
    !* tonic, either increasing or decreasing.  The returned value is not less than 0, nor      *
    !* greater than n-1.                                                                        *
    !********************************************************************************************
    implicit none

    integer,intent(in):: n, dj
    integer,intent(out):: cor, j
    integer,intent(inout):: jsav
    double precision,intent(in):: x,xx(*),yy(*) 
    integer ju, jm, jl, mm
    logical ascnd

    !mm is the number of abscissas used to bracket the interpolation value
    !note for linear interpolation, mm is always equal to 2
    mm=2

    if ((n.lt.2).or.(mm.lt.2.).or.(mm.gt.n)) &
        print *,"LOCATE size error"

    ascnd=.FALSE.
    if (xx(n).ge.xx(1)) ascnd=.TRUE. !True if ascending order table, false otherwise

    jl=1; ju=n;           !Initialize lower, upper bound
                        
    do 
        if((ju-jl).le.1) exit
        jm=(ju+jl)/2            !Compute a midpoint
        if((x.ge.xx(jm)).eq.ascnd) then  !Replace either the lower limit or upper limit
            jl=jm
        else
            ju=jm
        endif
    enddo                     !Repeat until the test condition is satisfied

    !Set the flag that determines which interpolation method to use in the next time step
    cor=1;  if(abs(jl-jsav).gt.dj) cor=0; 

    jsav=jl  !Save the lower bound

    j=max(1,min(n-mm+1,jl-((mm-2)/2))) !find the lower bound, must be between 1 and n-1

end subroutine

!--------------------------------------------------------------------------------------------

subroutine hunt(x,n,xx,yy,dj,cor,j,jsav)

    !********************************************************************************************
    !* Given a value x, return a value j such that x is (insofar as possible) centered in the  **
    !* subrange xx(j..j+mm-1), where xx is the stored data array.  The values in xx must be    **
    !* monatonic, either increasing or decreasing.  The returned value is not less than 0, nor **
    !* greater than n-1.                                                                       **
    !********************************************************************************************
    implicit none

    integer,intent(in):: n, dj
    integer,intent(out):: cor, j
    integer,intent(inout):: jsav
    double precision,intent(in):: x, xx(n), yy(n)
    integer ju, jm, jl, mm, inc
    logical ascnd

    inc=1  !Set the initial increment

    !mm is the number of abscissas used to bracket the interpolation value
    !note for linear interpolation, mm is always equal to 2
    mm=2

    if ((n.lt.2).or.(mm.lt.2.).or.(mm.gt.n)) &
     print *,"LOCATE size error"

    ascnd=.FALSE.
    if (xx(n).ge.xx(1)) ascnd=.TRUE. !True if ascending order table, false otherwise

    if((jl.lt.1).or.(jl.gt.n)) then  !Check to see if the guess values are useful
        jl=1  !if not, define the lower and upper bounds as the min/max abscissas of the array
        ju=n
        else
        if(x.ge.xx(jl).eq.ascnd) then
            do
                ju=jl+inc
                if(ju.ge.n) then
                    ju=n
                    exit    !off end of table
                elseif((x.lt.xx(ju)).eq.ascnd) then
                    exit    !Found bracket
                else
                    jl=ju
                    inc=inc+inc  !Double increment and try again
                endif
            enddo
        else
            ju=jl
            do
                jl=jl-inc
                if(jl.le.1) then
                    jl=1
                    exit  !Off end of table
                elseif((x.ge.xx(jl)).eq.ascnd) then
                    exit  !Found bracket
                else
                    ju=jl
                    inc=inc+inc  !Double increment and try again
                endif
            enddo
        endif
    endif
    do while((ju-jl).gt.1)
        jm=(ju+jl)/2
        if((x.ge.xx(jm)).eq.ascnd) then
            jl=jm
        else
            ju=jm
        endif
    enddo

    !Set the flag that determines which interpolation method to use in the next time step
    cor=1;  if(abs(jl-jsav).gt.dj) cor=0; 

    jsav=jl  !Save the lower bound

    j=max(1,min(n-mm+1,jl-((mm-2)/2))) !find the lower bound, must be between 1 and n-1

end subroutine


integer function t_warn(X,XLO,XHI,substance)
use global_props
implicit none

double precision,intent(in)::X,XHI,XLO
character(len=*),intent(in)::substance

    fl_ct=fl_ct+1.
    if((fl_ct.gt.50.).and.(fl_flag.eq.0.)) then
    
        write (LuFlEr,fmt=100) trim(substance),X,XLO,XHI,ttime
        fl_flag = 1.  !Note that the warning has now been written. Don't write again.
        
    endif
    100   format("WARNING: The property range for ",A,&
                 " has been exceeded more than 50 times. The value provided was ",F6.1,&
                 ", and the lower and upper limits are ",&
                 F6.1," and ",F6.1,", respectively. (time: ",F6.1,")")
    t_warn=0
end function
