        subroutine inittab()
C	Read the data tables of fragmentation functions and interpolate.
        implicit real*8 (a-h, o-z)

	parameter (iqbin=77, ixbin=70)
        parameter (meqn=28, mxgrd=70)!actually need only 6*4, not 7*4

C	For the interpolations:
	real*8 q2bin(77), xbin(70)
	integer NA(2)

	real*8 Fur(77, 70), Fdr(77, 70), Fsr(77, 70), Fcr(77, 70),
     1			Fbr(77, 70), Fgr(77, 70) ! rho
	real*8 Fuk(77, 70), Fdk(77, 70), Fsk(77, 70), Fck(77, 70),
     1			Fbk(77, 70), Fgk(77, 70) ! kstar
	real*8 Fuo(77, 70), Fdo(77, 70), Fso(77, 70), Fco(77, 70),
     1			Fbo(77, 70), Fgo(77, 70) ! omega
	real*8 Fup(77, 70), Fdp(77, 70), Fsp(77, 70), Fcp(77, 70),
     1			Fbp(77, 70), Fgp(77, 70) ! phi

 	common / bin /q2bin, xbin
	common / rhoff/  Fur, Fdr, Fsr, Fcr, Fbr, Fgr
	common / kstarff/  Fuk, Fdk, Fsk, Fck, Fbk, Fgk
	common / omegaff/  Fuo, Fdo, Fso, Fco, Fbo, Fgo
	common / phiff/  Fup, Fdp, Fsp, Fcp, Fbp, Fgp
	common / NAA/  NA

C	To use the tables to determine the fits to the data:

C	Suggest setting min of (x,Q2) = (10^-3, 1.6129) (charm)
C	Here mxgrd is the number of x values in the tables and meqn goes
C	over V,u,d,s,c,b,g=7 where V, valence, has not been computed
C	because observed data are for (h+hbar) hadrons.

C	Read the tables.
	open(12, file = 'rho.tab')
	open(13, file = 'kstar.tab')
	open(14, file = 'omega.tab')
	open(16, file = 'phi.tab')

C	(Q2, x) values at which fragmentation functions are defined:
	NA(1) = iqbin ! q2
	NA(2) = ixbin ! x
C	print*, 'Q2, x-bins', NA(1), NA(2)
	do i = 1, NA(1)
	  do j = 1, NA(2)
 	   read(12, *) q2bin(i), xbin(j)
C 	   write(*, *) 'q2, xbin', i, j, q2bin(i), xbin(j)
	  end do
 	  read(12, *)  ! gnuplot blanks
 	  read(12, *)  ! "
	end do
	rewind(12) ! so that can read data from it

C	Read (Q2, x) data for all hadrons:
	do i = 1, NA(1)
 	  do j = 1, NA(2)
	    read(12, *) qq, xx, Fur(i,j), Fdr(i,j), Fsr(i,j), Fcr(i,j),	
     1				 Fbr(i,j), Fgr(i,j)
	    read(13, *) qq, xx, Fuk(i,j), Fdk(i,j), Fsk(i,j), Fck(i,j),	
     1				 Fbk(i,j), Fgk(i,j)
	    read(14, *) qq, xx, Fuo(i,j), Fdo(i,j), Fso(i,j), Fco(i,j),	
     1				 Fbo(i,j), Fgo(i,j)
	    read(16, *) qq, xx, Fup(i,j), Fdp(i,j), Fsp(i,j), Fcp(i,j),	
     1				 Fbp(i,j), Fgp(i,j)

	  end do
 	  read(12, *)  ! gnuplot blanks
 	  read(12, *)  ! "
 	  read(13, *)  ! gnuplot blanks
 	  read(13, *)  ! "
 	  read(14, *)  ! gnuplot blanks
 	  read(14, *)  ! "
 	  read(16, *)  ! gnuplot blanks
 	  read(16, *)  ! "
	end do

C	Now have 2D tables of all ffs for all hadrons. Use these to
C	generate the D(x) at interpolated Q2 for each data point.
	
	close(12)
	close(13)
	close(14)
	close(16)

	return
	end

