module datafile_type
       type datafile
          private
          real, dimension(:,:), pointer :: dat
          integer :: m, n, nFile
          character*25, dimension(:), pointer :: headers
          logical :: contains_headers
          character*255 :: fname
       end type datafile

       interface get
          module procedure get_by_indices
          module procedure get_by_header
       end interface

       interface get_column
          module procedure get_column_by_indices
          module procedure get_column_by_header
       end interface

       interface get_header
          module procedure get_header_datafile
       end interface

       interface open
          module procedure open_file
       end interface

       interface isgood
          module procedure is_good_datafile
       end interface

       interface numrows
          module procedure numrows_datafile
       end interface

       interface numcols
          module procedure numcols_datafile
       end interface

       interface freeunit
          module procedure freeunit_
       end interface

       interface getfname
          module procedure getfname_
       end interface

       interface close
          module procedure close_datafile
       end interface

       interface dump
          module procedure dump_datafile
       end interface

       interface interpolate
          module procedure interpolate
       end interface
       
contains

       function getfname_(d)
              type(datafile), intent(in) :: d
              character*(len_trim(d%fname)) :: getfname_
              if (isgood(d)) then
                 getfname_ = d%fname
              else
                 stop "datafile_type::getfname: Uninitialized data file."
              end if
       end function getfname_

       function freeunit_()
              integer :: freeunit_
              logical :: openstatus

              freeunit_ = 10
              inquire(freeunit_, opened = openstatus)

              do while (openstatus)
                 freeunit_ = freeunit_ + 1
                 inquire(freeunit_, opened = openstatus)
              end do
       end function freeunit_

       function is_good_datafile(d)
              type(datafile), intent(in) :: d
              logical :: is_good_datafile

              is_good_datafile = .true.

              if (associated(d%dat) .and. d%m>0 .and. d%n>0) return

              is_good_datafile = .false.
       end function is_good_datafile

       function numrows_datafile(d)
              type(datafile), intent(in) :: d
              integer :: numrows_datafile

              numrows_datafile = d%m
       end function numrows_datafile

       function numcols_datafile(d)
              type(datafile), intent(in) :: d
              integer :: numcols_datafile

              numcols_datafile = d%n
       end function numcols_datafile

       subroutine get_column_by_header (d, v, head)
              integer :: j
              type(datafile), intent(in) :: d
              real, pointer, dimension(:) :: v
              character*(len(head)), intent(in) :: head

              if (.not. isgood(d)) then
                 print *,'datafile::get_column_by_header: Data array not allocated. Aborting.'
                 call dump(d)
                 stop
              end if

              allocate(v(numrows(d)))
              do j=1,numcols(d)
                 if (head .eq. d%headers(j)) then
                    v = d%dat(1:numrows(d),j)
                    return
                 end if
              end do

              print *,'datafile::get_column_by_header: Invalid header name "'//head//'". Aborting.'
              stop
       end subroutine get_column_by_header
       
       subroutine get_column_by_indices (d, v, j)
              integer, intent(in) :: j
              type(datafile), intent(in) :: d
              real, pointer, dimension(:) :: v

              if (j<1 .or. j>numcols(d)) then
                 print *,'datafile::get_column_by_indices: Invalid coordinate. Aborting.',j
                 call dump(d)
                 stop
              end if

              if (.not. isgood(d)) then
                 print *,'datafile::get_column_by_indices: Data array not allocated. Aborting.'
                 call dump(d)
                 stop
              end if

              allocate(v(numrows(d)))
              v = d%dat(1:numrows(d),j)
       end subroutine get_column_by_indices
       
       function get_by_indices(d, i, j)
              integer, intent(in) :: i, j
              type(datafile), intent(in) :: d
              real :: get_by_indices

              if (i<1 .or. i>numrows(d) .or. j<1 .or. j>numcols(d)) then
                 print *,'datafile::get_by_indices: Invalid coordinate. Aborting.',i,j
                 call dump(d)
                 stop
              end if

              if (.not. isgood(d)) then
                 print *,'datafile::get_by_indices: Data array not allocated. Aborting.'
                 call dump(d)
                 stop
              end if

              get_by_indices = d%dat(i,j)
       end function get_by_indices

       function get_by_header(d, head, i)
              character*15, intent(in) :: head
              integer, intent(in) :: i
              type(datafile), intent(in) :: d
              real :: get_by_header

              integer :: j = 1

              if (i<1 .or. i>numrows(d) .or. j<1 .or. j>numcols(d)) then
                 write(*,*)'datafile::get_by_header: Invalid coordinate. Aborting.',i,j
                 call dump(d)
                 stop
              end if

              if (.not. isgood(d)) then
                 print *,'datafile::get_by_header: Data array not allocated. Aborting.'
                 call dump(d)
                 stop
              end if

              do j=1,numcols(d)
                 if (trim(head) .eq. trim(d%headers(j))) then
                    get_by_header = d%dat(i,j)
                    return
                 end if
              end do

              print *,'datafile::get_by_header: Invalid header name. Aborting.'
              stop
       end function get_by_header

       function get_header_datafile(d, n)
              type(datafile), intent(in) :: d
              integer, intent(in) :: n
              character*25 :: get_header_datafile

              if (n<1 .or. n>numcols(d)) then
                 print *,'datafile::get_header_datafile: Invalid coordinate. Aborting.'
                 stop
              end if

              if (.not. isgood(d)) then
                 print *,'datafile::get_header_datafile: Data array not allocated. Aborting.'
                 stop
              end if

              get_header_datafile = d%headers(n)
       end function get_header_datafile

       subroutine open_file(d, filename, no_halt)
              type(datafile), intent(out) :: d
              character*(len(filename)), intent(in) :: filename
              logical, optional :: no_halt
              integer :: i, j, nStatus, nTemp, nStat2
              character :: chTemp*1
              real, allocatable :: rTemp(:)

              write(*,'(A)',advance='no')"Loading data file "//filename//"..."

              d%nFile = freeunit()
              open( unit = d%nFile, file = filename, status = "old", action = "read", &
                   position = "rewind", iostat = nStatus )

              if(nStatus > 0) then
                 if (present(no_halt) .and. no_halt) return
                 write (*,*) "datafile_type:open_file: failed to open input file ",filename
                 stop
              end if

              read(unit=d%nFile,fmt='(a1)',iostat=nStatus,advance='no') chTemp

              d%fname = filename

              if (chTemp==' ') then
                 do while ((chTemp==' ').and.(nStatus==0))
                    read(unit=d%nFile,fmt='(a1)',iostat=nStatus,advance='no') chTemp
                 end do
                 if (nStatus.ne.0) then
                    if (present(no_halt) .and. no_halt) return
                    write(*,*)"datafile::open_file: error reading first line of input file ",filename
                    stop
                 end if
              end if

              if (((chTemp>='0').and.(chTemp<='9')).or.(chTemp=='-').or.(chtemp=='+')) then
                 d%contains_headers = .false.
              else
                 d%contains_headers = .true.
              end if
              rewind(d%nFile)

              nTemp = 0
              nStatus=0

              read(unit=d%nFile,fmt='(a1)',iostat=nStatus,advance='no') chTemp
              do while ((chTemp.eq.' ').and.(nStatus==0))
                 read(unit=d%nFile,fmt='(a1)',iostat=nStatus,advance='no') chTemp
              end do
              do while(nStatus==0)
                 nTemp = nTemp + 1
                 do while ((chTemp.ne.' ').and.(nStatus==0))
                    read(unit=d%nFile,fmt='(a1)',iostat=nStatus&
                         ,advance='no') chTemp
                 end do
                 if (nStatus==0)  then
                    read(unit=d%nFile,fmt='(a1)',iostat=nStatus,advance='no') chTemp
                    do while (chTemp==' '.and.nStatus==0)
                       read(unit=d%nFile,fmt='(a1)',iostat=nStatus,advance='no') chTemp
                    end do
                 end if
              end do

              d%n=nTemp
              if (d%n==0) stop "datafile::open_file: no columns found"
              i = 0

              if(d%contains_headers) then
                 allocate(d%headers(numcols(d)))
                 do j=1,numcols(d)
                    d%headers(j)=''
                    nTemp=0
                    read(unit=d%nFile,fmt='(a1)',advance='no',iostat=nStat2) chTemp
                    do while (nStat2==0.and.chTemp==' ')
                       read(unit=d%nFile,fmt='(a1)',advance='no',iostat=nStat2) chTemp
                    end do
                    do while (nStat2==0.and.chTemp.ne.' ')
                       if (nTemp.gt.0) then
                          d%headers(j) = d%headers(j)(1:nTemp) // chTemp
                       else
                          d%headers(j) = chTemp
                       end if
                       nTemp = nTemp + 1
                       read(unit=d%nFile,fmt='(a1)',advance='no',iostat=nStat2) chTemp
                    end do
                 end do
              end if

              i=0
              rewind(d%nFile)
              if(d%contains_headers) read(d%nFile,*)
              allocate(rTemp(numcols(d)))
              read (unit=d%nFile,iostat=nStatus,fmt=*) (rTemp(j),j=1,numcols(d))
              do while (nStatus==0)
                 i=i+1
                 read (unit=d%nFile,iostat=nStatus,fmt=*) (rTemp(j),j=1,numcols(d))
              end do

              d%m=i
              allocate (d%dat(numrows(d),numcols(d)))

              rewind(unit=d%nFile)

              if(d%contains_headers) read(d%nFile,*)

              do i=1,numrows(d)
                 read (unit=d%nFile,fmt=*) (d%dat(i,j),j=1,numcols(d))
              end do

              if (.not. isgood(d)) then
                 write(*,*)"datafile_type:open_file: bad data file."
                 call dump_datafile(d)
                 stop
              end if
              write(*,*) "ok"
       end subroutine open_file

       subroutine close_datafile(d)
              type(datafile), intent(inout) :: d

              close(d%nFile)
       end subroutine close_datafile

       subroutine dump_datafile(d)
              type(datafile), intent(in) :: d
              integer :: i,j

              if (.not. isgood(d))  then
                 write(*,'(A)') "datafile_type::dump: Data array not allocated."
                 return
              else
                 write(*,'(/A80,/i5," x ",i3)') d%fname,d%m,d%n
              end if

              if (d%contains_headers) then
                 do j=1,numcols(d)
                    write(*,'(A7,1x)', advance='no') d%headers(j)
                 end do
                 write(*,'(/)')
              end if

              do i=1, numrows(d)
                 do j=1, numcols(d)
                    write(*,'(g7.2e2,1x)',advance='no') d%dat(i,j)
                 end do
                 write(*,'(/)',advance='no')
              end do
       end subroutine dump_datafile

       real function interpolate(d, E, col)
               integer, intent(in) :: col
               real, intent(in) :: E
               type(datafile) :: d

               integer :: i

               if (.not. isgood(d)) stop "datafile_type::interpolate: datafile not initialized"

               i = 1 
               do while (E>get(d,i,1) .and. i<numrows(d))
                  i= i + 1
               end do

               if (E>get(d,i,1) .or. i==1) stop "datafile_type::interpolate: invalid key"

               if (col>numcols(d)) stop "datafile_type::interpolate: invalid column"

               interpolate = (E-get(d,i-1,1))/(get(d,i,1)-get(d,i-1,1)) * (get(d,i,col)-get(d,i-1,col)) + get(d,i-1,col)
	end function interpolate
end module datafile_type