FUNCTION ftau,lnr,tt,AA5000=AA5000,LINEAR=LINEAR,ab=ab
	common cdat,x,y,z,nx,ny,nz,nw,ntmax,date0,time0
	common opacs,kaparr,nrho,nT

	if (N_ELEMENTS(kaparr) EQ 0) then begin
		nrho = 20
		nT   = 70
		kaparr = dblarr(2, nT, nrho)
;		datdir = 'opac/c38/'
;		hosts_dir, scrdir=datdir
;		openr, 5, datdir+'kap_solar'+'.dat'
		openr, 5, 'kap.dat'
		readu, 5, kaparr
		close, 5
	endif
	s=size(lnr)
	itype = s(0)+1
	if KEYWORD_SET(AA5000) then ivk=1   else ivk=0
	if KEYWORD_SET(LINEAR) then i12=0.0 else i12 = 1./12.
  if (s(0) EQ 3) then begin
	mz=s(3)
	nx=s(1) & ny=s(2) & nz=s(3) & nw=nx*ny
	print,'absorption ...'
	lnab=lnr
;	for n=0,nz-1 do begin
	irho = (lnr/alog(1e1)-7.+14.)/11.*(nrho-1.0)
	itt  = (alog10(tt)-3.2)/1.4*(nT-1.0)
	for n=0,nz-1 do begin
;		lnab(*,*,n)=(bispl(reform(kaparr(ivk,*,*)), $
;			reform(itt(*,*,n)), reform(irho(*,*,n)) )+1.)*alog(1e1) + lnr(*,*,n)
		lnab(*,*,n)=(interpolate(reform(kaparr(ivk,*,*)), $
			reform(itt(*,*,n)), reform(irho(*,*,n)) )+1.)*alog(1e1) + lnr(*,*,n)
	endfor
	ab=exp(lnab)
	dlnab=make_array(nx,ny,mz,type=s(itype))
	print,'dlnab ...'
; dlnab is twice the log der
	for n=1,mz-2 do dlnab(*,*,n)=lnab(*,*,n+1)-lnab(*,*,n-1)
	dlnab(*,*,0)=2.*(lnab(*,*,1)-lnab(*,*,0))
	dlnab(*,*,mz-1)=2.*(lnab(*,*,mz-1)-lnab(*,*,mz-2))
;	dlnab = zderz(lnab,z)
	tau1=make_array(nx,ny,mz,type=s(itype))
	tau1(*,*,0)=0.
	print,'tau ...'
	for n=1,mz-1 do begin
	  dz=z(n)-z(n-1)
	  tau1(*,*,n)=tau1(*,*,n-1)+dz* $
	    (ab(*,*,n-1)*0.5*(1.0+i12*dlnab(*,*,n-1))+ $
	     ab(*,*,n  )*0.5*(1.0-i12*dlnab(*,*,n  )))
	endfor
	tau1(*,*,0)=1.e1^(alog10(tau1(*,*,1)) - $
			(alog10(tau1(*,*,2))-alog10(tau1(*,*,1)))/(z(2)-z(1))*(z(1)-z(0)))
  endif
  if (s(0) EQ 1) then begin
	mz=s(1)
	print,'absorption ...'
	lnab=lnr
	irho = (lnr/alog(1e1)-7.+14.)/11.*(nrho-1.0)
	itt  = (alog10(tt)-3.2)/1.4*(nT-1.0)
	for n=0,nz-1 do begin
;		lnab=(bispl(reform(kaparr(ivk,*,*)),itt,irho)+1.)*alog(1e1) + lnr
		lnab=(interpolate(reform(kaparr(ivk,*,*)),itt,irho)+1.)*alog(1e1) + lnr
	endfor
	ab=exp(lnab)
	dlnab=make_array(mz,type=s(itype))
	print,'dlnab ...'
	; dlnab is twice the log der
	for n=1,mz-2 do dlnab(n)=lnab(n+1)-lnab(n-1)
	dlnab(0)=2.*(lnab(1)-lnab(0))
	dlnab(mz-1)=2.*(lnab(mz-1)-lnab(mz-2))
	tau1=make_array(mz,type=s(itype))
	tau1(0)=0.
	print,'tau ...'
	for n=1,mz-1 do begin
	  dz=z(n)-z(n-1)
	  tau1(n)=tau1(n-1)+dz* $
	    (ab(n-1)*0.5*(1.0+i12*dlnab(n-1))+ $
	     ab(n  )*0.5*(1.0-i12*dlnab(n  )))
	endfor
	tau1(0)=1.e1^(alog10(tau1(1)) - $
			(alog10(tau1(2))-alog10(tau1(1)))/(z(2)-z(1))*(z(1)-z(0)))
  endif
	return,tau1
end