subroutine fit_dat(filenames)

	character*(*) filenames

	integer iret, fit_dat_opt

	iret=fit_dat_opt(filenames, 0, 0.0) ! option 0: standard
	end


	subroutine fit_link(filenames)

	character*(*) filenames

	integer iret, fit_dat_opt

	iret=fit_dat_opt(filenames, 1, 0.0) ! option 1: link to existing data
	end


	subroutine fit_dat_merge(filenames, step)

	character*(*) filenames
	real step

	integer iret, fit_dat_opt

	iret=fit_dat_opt(filenames, 2, step) ! option 2: merge read data
	end


	subroutine fit_open(filenames)

	character*(*) filenames

	integer iret, fit_dat_opt

	iret=fit_dat_opt(filenames, 4, 0.0) ! option 4: <return> reads last.fit3
	end


	subroutine fit_dat_next(filenames, flag)

	character*(*) filenames
	integer flag  ! input: 0: start read, 1: read next file, 2: read next range
	              ! output: 0: end of list reached, 1: o.k., no range, 2: o.k., next range


	call fit_dat_next_opt(filenames, flag, 0, 0.0)
	end


	subroutine fit_dat_next_opt(filenames, flag, opt, step)

	character*(*) filenames
	integer flag  ! input: 0: start read, 1: read next file, 2: read next range
	              ! output: 0: end of list reached, 1: o.k., no range, 2: o.k., next range
	integer opt  ! 0: DAT, 1: LINK, 2: DAT&MERGE
	real step

	integer lopt, fit_dat_opt
	integer iret
	integer idx/0/
	
	logical silent, silent0
	common/fit_dat_com/silent

	if (flag .eq. 0) then
	  call dat_set_index(1)
	  lopt=8  ! read 1st
	else if (flag .eq. 1) then
	  call dat_set_index(1)
	  lopt=16 ! read next file
	else
	  idx=idx+1
	  lopt=24 ! read again
	endif
	silent0=silent
	iret=fit_dat_opt(filenames, lopt+opt, step)
	if (iret .eq. 0 .and. flag .eq. 2) then
	   call dat_set_index(1)
	   silent=silent0
	   iret=fit_dat_opt(filenames, 16+opt, step) ! read next file
	endif
	if (iret .eq. 0) then
	  flag=0
	else
	  call dat_next_index(flag)
	endif
	end


	integer function fit_dat_opt(filenames, opt, step)
!	--------------------------------------------------

	implicit none

	include 'fit.inc'

	character filenames*(*)
	integer opt
	real step

	integer n_mul
	parameter (n_mul=256)
	character filelist*(maxflen), sample*80
	integer i,j,l,ll,nread,n1
	integer merge, list_mode
	integer fidx, sets/0/
	logical init, init_dat/.true./, purged
	integer nmax, ncol, nset0
	external fit_putval, fit_fit3_read
	real stp, t1, dt1, w1, tmean, wmean, ds, upar

	logical silent
	common/fit_dat_com/silent
	data silent/.false./

! setup different flags
	init=mod(opt,2) .eq. 0     ! bit0=0: purge old data
	merge=mod(opt,4)/2         ! bit1=1: merge data ?
	if (merge .ne. 0) stp=step ! step for merge 
	load_state=mod(opt,8)/4	   ! bit2: do a load for FitSave files
	list_mode=mod(opt,32)/8    ! bit3+4: 1: read first file, 2: read next file (ignore filenames argument), 3: read again

	fit_dat_opt=1
	if (list_mode .gt. 0) then
	  if (list_mode .ne. 3) then
	    call dat_init_list
	    if (list_mode .eq. 1) then
	      call dat_filelist(filenames)
	    endif
	  endif
	  goto 15
	endif
	fidx=1
11	if (filenames(fidx:fidx) .eq. '+') then
	  fidx=fidx+1
	  if (load_state .eq. 0) init=.false.
	  goto 11
	elseif (filenames(fidx:fidx) .eq. '&') then
	  fidx=fidx+1
	  merge=1
	  goto 11
	endif

	if (init_dat) then
	  call dat_init
	  call dat_fit3_replace(fit_fit3_read)
	  init_dat=.false.
	endif
	if (filenames(fidx:) .eq. ' ') then
	  if (load_state .eq. 1) then
	    call dat_ask_filelist(filelist
	1   , 'Hit <RETURN> to continue with last data')
	  else
	    call dat_ask_filelist(filelist, ' ')
	  endif
	  if (filelist .eq. ' ') then
	    if (load_state .eq. 1) then
	      filelist='last.fit3'
	    else
	      fit_dat_opt=0
	      goto 91
	    endif
	  endif
	  fidx=1
12	  if (filelist(fidx:fidx) .eq. '+') then
	    fidx=fidx+1
	    if (load_state .eq. 0) init=.false.
	    goto 12
	  elseif (filelist(fidx:fidx) .eq. '&') then
	    fidx=fidx+1
	    merge=1
	    goto 12
	  endif
	  call str_trim(filelist,filelist(fidx:),ll)
	else
	  call str_trim(filelist,filenames(fidx:),ll)
	endif	  


	if (init) then
	  sets=0
	  call dat_init_list
	else
	  call dat_put_filelist(fillis)
	endif

	call dat_filelist(filelist(1:ll))
	filesave=filelist(1:ll)

15	if (init) then
	  npkt=0
	  call sym_purge(1)
	  temp=0
	  dtemp=0
	  wtemp=0
	endif

	purged=.false.
	n1=npkt+1
30	if (silent) call dat_silent
	call fit_dat_init_table
	call fit_ncurves(0)
	if (list_mode .ne. 3) then
	  call dat_read_next(fit_putval, maxdat-npkt, nread
	1   , xval(npkt+1), YVAL(npkt+1), sig(npkt+1), rmon(npkt+1))
	else
	  call dat_read_again(fit_putval, maxdat-npkt, nread
	1   , xval(npkt+1), YVAL(npkt+1), sig(npkt+1), rmon(npkt+1))
	endif
	if (nread .gt. 0) then
	  if (npkt .ne. 0) then
	    if (ymon0 .ne. ymon) then
	      if (ymon .eq. 0) then
	        ymon=ymon0
	      elseif (ymon0 .eq. 0) then
	        if (ymon .ne. 0) then
	          do i=npkt+1,npkt+nread
	            rmon(i)=ymon
	          enddo
	        endif
	      else
	        do i=npkt+1,npkt+nread
	          YVAL(i)=YVAL(i)*ymon/ymon0
	          sig(i)=sig(i)*ymon/ymon0
	        enddo
	      endif
	    endif
	  else
	    nset=0
	    ymon=ymon0
	  endif

	  w1=0
	  j=1
	  nset0=nset
	  nmax=npkt+nread
	  do while (npkt .lt. nmax)
	    ncol=cols(j)
	    do i=0,min(nmax-npkt,rows(j)*ncol-1)
	      if (npkt .lt. nmax) then
	        npkt=npkt+1
	        iset(npkt)=nset+mod(i,ncol)+1
	        w1=w1+rmon(npkt) ! weight for temperature stddev
	      endif
	    enddo
	    nset=nset+ncol
	    j=min(nmult,j+1)
	  enddo
	  t1=0
	  dt1=0
	  call sym_get_real('Temp', t1)
	  if (t1 .eq. 0) then
	    call sym_get_real('sTemp', t1) ! this is for old DMC/HRPT files
	  endif
	  call sym_get_real('dTemp', dt1)

	  if (t1 .ne. 0 .and. w1 .ne. 0) then ! calculate merged temperature stddev
	    wmean=wtemp+w1
	    tmean=(temp*wtemp+t1*w1)/wmean
	    dtemp=sqrt((wtemp*(dtemp**2+(temp-tmean)**2)
	1                + w1*(  dt1**2+  (t1-tmean)**2))/wmean )
	    temp=tmean
	    wtemp=wmean
	  endif

	  if (nset .gt. maxset) goto 31

	  do i=1,npd
	    upar=none
	    call sym_get_real(userpar(usernp+i), upar)
	    if (upar .eq. none) then
	      if (pdpar(i, nset) .eq. none .and.
	1         userpar(usernp+i) .ne. 'Numor' .and.
	1         userpar(usernp+i) .ne. 'Temp') then 
	        print *,'parameter ',userpar(usernp+i),' not found'
	        pdpar(i, nset)=0.0
	      endif
	    else
	      if (userpar(usernp+i) .eq. 'Numor') then
	        if (nset0+1 .eq. nset) then
	          ds=0.0
	          call sym_get_real('Pal', ds)
	        else
	          ds=0.1
	        endif
	        do j=nset0+1,nset
	          pdpar(i, j)=upar+ds
	          ds=ds+0.1
	        enddo
	      else
	        do j=nset0+1,nset
	          pdpar(i, j)=upar
	        enddo
	      endif
	    endif
	  enddo

31	  continue  
	  ISW(1) = 0
	  ISW(2) = 0
	  ISW(3) = 0
	  NFCN = 1
	  SIGMA = 0.

          if (ififu .eq. 0) ififu=8       ! set to plot only mode
c	  print '(x,a,i7)','Npkt=',nread
	  call sym_show(1)
	  call meta_put('Npkt', float(nread))
	  call meta_set_format('Npkt', 5003)
	  call meta_set_format('Numor', 5003)
	  call meta_set_format('Counts', 5003)
	  if (.not. silent) call sym_list(isyswr,2,1,' ')
	  sets=sets+1
	  if (merge .ne. 0 .and. sets .gt. 1) then
	    if (npkt .gt. n1 .and.
	1       (npkt .gt. maxdat/2 .or. nset .ge. maxset)) then ! probably not enough space -> merge
	      nxmin=n1
	      nxmax=npkt
	      call fit_merge(stp)
	      merge=1
	    else
	      merge=2
	    endif
	  endif

	  if (keepregion) then
	    if (npkt .gt. n1 .and. npkt .gt. maxdat/2) then ! probably not enough space -> purge
	      call fit_restore_region(.true.)
	      purged=.true.
	    endif
	  endif

	  if (list_mode .eq. 0) then
	    if (npkt .lt. maxdat) goto 30
	    print *,'no more datapoints accepted'
	  endif
	else if (list_mode .ne. 0) then
	  fit_dat_opt=0
	endif

	if (merge .eq. 2 .and. npkt .ge. n1) then ! merge (if not yet done)
	  nxmin=n1
	  nxmax=npkt
	  call fit_merge(stp)
	endif
	if (keepregion) then
	  call fit_restore_region(purged)
	  if (purged) then
	    print *,'points outside window/region were purged'
	  endif
	else
	  call fit_win(1.0,1.0) ! disable window
	endif

	if (load_state .gt. 1) then
	  if (load_state .eq. 3) then ! for FitSave 3.2 and older
	    call dat_init_list
	    call dat_filelist(filnam)
	    load_state=0
	    goto 30
	  endif
	else
	  call dat_get_filelist(fillis)
	  call str_trim(fillis, fillis, l)
	  if (.not. silent) then
	    print *
	    if (list_mode .eq. 0) then
	      if (fillis(1:l) .eq. ' ') then
	        print *,'no file read'
	      else
	        print *,'Files read: ',fillis(1:l)
	      endif
	    endif
	  endif

	  wavlen=0
	  call sym_get_real('lambda', wavlen)
	  call sym_get_str('Title', ll, itit)
	  sample=' '
	  call sym_get_str('sample', l, sample)
	  if (sample(1:l) .ne. ' '
	1     .and. index(itit(1:ll),sample(1:l)) .eq. 0) then
	    if (ll+l .lt. len(itit)) then
	      if (ll+l+7 .lt. len(itit)) then
	        itit(ll+1:)=' Sample='
	        ll=ll+8
	      else
	        ll=ll+1
	      endif
	      itit(ll+1:)=sample
	    endif
	  endif

	  call sym_show(1)
	  call sym_put_real('Temp', temp)
	  call sym_put_real('dTemp', dtemp)
	  call sym_put_real('Monitor', ymon)

	  if (npkt .gt. 0) call fit_set(0,0.0,0.0,0.0,0.0)
	endif
91	silent=.false.
	end


	subroutine fit_dat_silent

	logical silent
	common/fit_dat_com/silent

	silent=.true.
	end


	subroutine fit_datraw
!	---------------------
	print *,'Subroutine FIT_DATRAW no longer available'
	end


      	subroutine fit_spec
!       -------------------
	print *,'Subroutine FIT_SPEC is no longer available'
	END



	subroutine fit_dat_put(mode, xx, nx, yy, ny, ss, ns, ww, nw)  !!
!!	------------------------------------------------------------

! preconditions: (ny>0) and (nx=ny or nx=1 or nx=2) and (ns=ny or ns=1) and (nw=ny or nw=1)
! nx=1: x-values are xx(1),xx(1)+1,xx(1)+2,...,xx(1)+(nx-1)
! nx=2: x-values are xx(1),xx(2),...,xx(1)+(xx(2)-xx(1))*(nx-1)
! ns=1: ss(1) <> 0.0: sigma values are ss(1)
!       ss(1) =  0.0: sigma values are sqrt(max(1.0,yy(i)))
! nw=1: weights are ww(1)

	include 'fit.inc'

	integer mode !! mode=0: purge before, mode=1: link new dataset, mode=2: link to existing dataset
	integer nx, ny, ns, nw !!
	real xx(nx), yy(ny), ss(ns), ww(nw) !!

	integer i,j,nread
	real x0,xs

	if (mode .eq. 0) then
	  npkt=0
	  call sym_purge(1)
	endif

	if (ny .lt. 0) stop 'FIT_DAT_INT: illegal ny'

	nread=ny
	if (npkt+ny .gt. maxdat) then
	  print *,'Too many datapoints: truncated'
	  nread=maxdat-npkt
	endif

	j=npkt
	do i=1,nread
	  j=j+1
	  YVAL(j)=yy(i)
	enddo

	j=npkt
	if (nx .eq. ny) then
	  do i=1,nread
	    j=j+1
	    xval(j)=xx(i)
	  enddo
	else
	  x0=xx(1)
	  if (nx .eq. 1) then
	    xs=1
	  elseif (nx .eq. 2) then
	    xs=(xval(nx)-x0)/(nx-1)
	  else
	    stop 'FIT_DAT_INT: illegal nx'
	  endif
	  x0=x0-xs
	  do i=nx+1,nread
	    j=j+1
	    xval(j)=x0+xs*i
	  enddo
	endif

	j=npkt
	if (ns .eq. ny) then
	  do i=1,nread
	    j=j+1
	    sig(j)=ss(i)
	  enddo
	elseif (ns .eq. 1) then
	  xs=ss(1)
	  if (xs .eq. 0.0) then
	    do i=1,nread
	      j=j+1
	      sig(j)=sqrt(max(1.0,YVAL(j)))
	    enddo
	  else
	    do i=1,nread
	      j=j+1
	      sig(j)=xs
	    enddo
	  endif
	else
	  stop 'FIT_DAT_INT: illegal ns'
	endif

	j=npkt
	if (nw .eq. ny) then
	  ymon0=ww(1)
	  do i=1,nread
	    j=j+1
	    rmon(j)=ww(i)
	  enddo
	elseif (nw .eq. 1) then
	  ymon0=ww(1)
	  do i=1,nread
	    j=j+1
	    rmon(j)=ymon0
	  enddo
	else
	  stop 'FIT_DAT_INT: illegal nw'
	endif

	if (npkt .ne. 0) then
	  if (mode .eq. 1) nset=nset+1
	  if (ymon0 .ne. ymon) then
	    if (ymon .eq. 0) then
	      ymon=ymon0
	    elseif (ymon0 .eq. 0) then
	      if (ymon .ne. 0) then
	        do i=npkt+1,npkt+nread
	          rmon(i)=ymon
	        enddo
	      endif
	    else
	      do i=npkt+1,npkt+nread
	        YVAL(i)=YVAL(i)*ymon/ymon0
	        sig(i)=sig(i)*ymon/ymon0
	      enddo
	    endif
	  endif
	else
	  nset=1
	  ymon=ymon0
	endif
	do i=npkt+1,npkt+nread
	  iset(i)=nset
	enddo
	npkt=npkt+nread

	ISW(1) = 0
	ISW(2) = 0
	ISW(3) = 0
	NFCN = 1
	SIGMA = 0.

        nxmin=1
        nxmax=npkt
	end


	subroutine fit_putval(str, val)

	character*(*) str
	real val

	include 'fit.inc'

	if (str .eq. 'ShowLevel') then
	  call sym_show(nint(val))
	  return
	elseif (str .eq. 'Monitor' .or. str .eq. 'monitor') then
	  ymon0=val
	endif
	call meta_put(str, val)
	end


	subroutine fit_fit3_read(lun,forced,nread,putval,nmax,xx,yy,ss,ww)
!	------------------------------------------------------------------
	implicit none

	include 'fit.inc'

	integer lun, forced, nread, nmax
	external putval
	real xx(nmax), yy(nmax), ss(nmax), ww(nmax)

	character str*8

	if (load_state .eq. 1) then
	  read(lun, '(a)') str
	  rewind lun
	  if (str .eq. 'FitSave') then
	    load_state=2
	    call fit_load_nopn(lun)
	    nread=0
	    return
	  endif
	endif
	call dat_fit3_read(lun, forced, nread, putval, nmax,xx,yy,ss,ww)
	end


	subroutine fit_dat_options(options)

	include 'fit.inc'

	character*(*) options

	character typ*16, desc*256
	integer done
	integer l

	if (options .eq. '?' .or. options .eq. ' ') then
	  call dat_gettyp(typ)
	  call str_trim(typ, typ, l)
	  print *
          done = 1
	  call dat_desc_opt(done, typ(1:l))
	  if (done .eq. 0) then
	    print '(x,2a)','No options available for filetype ',typ(1:l)
	  else
	    call dat_get_def_options(desc)
	    if (desc .ne. ' ') then
	      call str_trim(desc, desc, l)
	      print '(x,2a)','Actual options: ',desc(1:l)
	    endif
	    if (options .eq. ' ') then
	      print *
	      print '(x,a)','Syntax: option1=value1,option2=value2 ...'
	      print '(x/x,a,$)','Enter default options: '
	      read(*,'(a)',end=9,err=9) desc
	      call dat_def_options(desc)
	    endif
	  endif
	else
	  call dat_def_options(options)
	endif
9	end


	subroutine fit_dat_table(itbl, icols, irows)
! to be called within dat_..._read routines
! table number itbl has icols columns and irows rows
	include 'fit.inc'

	integer itbl, irows, icols
	integer i

	if (itbl .gt. nmult) then
	  do i=nmult+1,itbl
	    rows(i)=0
	    cols(i)=1
	  enddo
	  nmult=itbl
	endif
	rows(itbl)=irows
	cols(itbl)=icols
	end


	subroutine fit_dat_init_table

	include 'fit.inc'

	nmult=1
	rows(1)=maxdat
	cols(1)=1
	end


	subroutine fit_dat_pdp_idx(name, ipdp)
	! to be called within dat_..._read routines
	include 'fit.inc'

	character name*(*)
	integer ipdp

	integer i
	character unam*32, upar*32

	call str_upcase(unam, name)
	do i=1,npd
	  call str_upcase(upar, userpar(usernp+i))
	  if (upar .eq. unam) then
	    ipdp=i
	    return
	  endif
	enddo
	if (usernp .lt. maxext) then
	  call fit_userpdp(name)
	  ipdp=npd
	else
	  ipdp=0
	endif
	end


	subroutine fit_dat_pdp_set(ipdp, idx, value)
	! to be called within dat_..._read routines

	include 'fit.inc'

	integer ipdp, idx
	real value

	if (ipdp .gt. 0 .and. ipdp .le. npd .and.
	1   idx .gt. 0 .and. nset+idx .le. maxset) then
	  pdpar(ipdp, nset+idx)=value
	endif
	end