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. iget(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