c     Copyright (C) 2008-2014 Vincent Eymet
c
c     KDISTRIBUTION is free software; you can redistribute it and/or modify
c     it under the terms of the GNU General Public License as published by
c     the Free Software Foundation; either version 3, or (at your option)
c     any later version.

c     KDISTRIBUTION is distributed in the hope that it will be useful,
c     but WITHOUT ANY WARRANTY; without even the implied warranty of
c     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
c     GNU General Public License for more details.

c     You should have received a copy of the GNU General Public License
c     along with KDISTRIBUTION; if not, see <http://www.gnu.org/licenses/>
c
      subroutine read_options(kind,kpts,Nq,fixed_g,alpha,beta,
     &     endpts,Np,epsilon,sit,quad_file)
      implicit none
      include 'max.inc'
      
      logical ex
      integer quad,i,kind,kpts,Nq,Np,Ni,ios
      double precision endpts(2)
      double precision fixed_g(2)
      double precision alpha,beta
      double precision epsilon
      integer sit
      character*(Nchar_mx) opfile,quad_file
      character*1 ch1,zch
      character*2 ch2
c     label
      integer strlen
      character*(Nchar_mx) label
      label='subroutine read_options'
 11   format(I1)
 12   format(I2)

      opfile='./options.in'
      open(10,file=opfile(1:strlen(opfile)),status='old',iostat=ios)
      if (ios.ne.0) then        ! file not found
         write(*,*) 'Error from routine ',
     &        label(1:strlen(label)),' :'
         write(*,*) 'File not found:'
         write(*,*) opfile(1:strlen(opfile))
         stop
      endif
      do i=1,3
         read(10,*)
      enddo
      read(10,*) quad
      read(10,*)
      read(10,*) Nq
      read(10,*)
      read(10,*) fixed_g(1)
      read(10,*)
      read(10,*) fixed_g(2)
      read(10,*)
      read(10,*) alpha
      read(10,*)
      read(10,*) beta
      read(10,*)
      read(10,*)
      read(10,*) Ni
      read(10,*)
      read(10,*) epsilon
      read(10,*)
      read(10,*) sit
      close(10)

      if ((quad.lt.0).or.(quad.gt.8)) then
         write(*,*) 'Error from routine ',
     &        label(1:strlen(label)),' :'
         write(*,*) 'Quadrature type is incorrect'
         stop
      endif

      if (Nq.gt.Nqmx) then
         call error(label)
         write(*,*) 'Nq=',Nq
         write(*,*) 'while Nqmx=',Nqmx
         stop
      endif

      kpts=0
      kind=quad
      if (quad.eq.0) then
         if (Nq.lt.10) then
            write(zch,11) 0
            write(ch1,11) Nq
            ch2=zch(1:strlen(zch))//ch1(1:strlen(ch1))
         else if (Nq.lt.100) then
            write(ch2,12) Nq
         else
            write(*,*) 'Error from routine ',
     &           label(1:strlen(label)),' :'
            write(*,*) 'quadrature order is too high'
            stop
         endif
         quad_file='./data/quad'//ch2(1:strlen(ch2))//'.in'
         inquire(file=quad_file(1:strlen(quad_file)),exist=ex)
         if (.not.ex) then
            write(*,*) 'Error from routine ',
     &           label(1:strlen(label)),' :'
            write(*,*) 'File could not be detected:'
            write(*,*) quad_file(1:strlen(quad_file))
            stop
         endif
      else if (quad.eq.1) then  ! Legendre
         kpts=0
      else if (quad.eq.2) then ! radau <=> g=0 is imposed
         kpts=1
      else if (quad.eq.3) then ! lobatto <=> g=0 and g=1 are imposed
         kpts=2
      endif

      if (Nq.le.0) then
         write(*,*) 'Error from routine ',
     &        label(1:strlen(label)),' :'
         write(*,*) 'Quadrature order is incorrect'
         stop
      endif

      if (((quad.eq.2).or.(quad.eq.3)).and.(fixed_g(1).lt.0.0D+0)) then
         write(*,*) 'Error from routine ',
     &        label(1:strlen(label)),' :'
         write(*,*) 'First g abscissa imposed is incorrect'
         stop
      endif
      if ((quad.eq.3).and.(fixed_g(2).gt.1.0D+0)) then
         write(*,*) 'Error from routine ',
     &        label(1:strlen(label)),' :'
         write(*,*) 'Seoncd g abscissa imposed is incorrect'
         stop
      endif

c     Modification of endpts
      do i=1,2
         endpts(i)=2.0D+0*fixed_g(i)-1.0D+0
      enddo

      Np=Ni+1
      if ((Np.lt.2).or.(Np.gt.Npmax)) then
         write(*,*) 'Error from routine ',
     &        label(1:strlen(label)),' :'
         write(*,*) 'Number of intervals is incorrect'
         write(*,*) 'Np=',Np,' Npmax=',Npmax
         stop
      endif

      if (epsilon.lt.0.0D+0) then
         write(*,*) 'Error from routine ',
     &        label(1:strlen(label)),' :'
         write(*,*) 'Value of epsilon is incorrect'
         stop
      endif

      if ((sit.ne.1).and.(sit.ne.2)) then
         write(*,*) 'Error from routine ',
     &        label(1:strlen(label)),' :'
         write(*,*) 'Incorrect value for:'
         write(*,*) 'Spline interpolation type'
         stop
      endif

      return
      end