FUNCTION tau,lnr,ee,xmu
	@common_cdata
	@common_cmesh
	; First load table: eos_table,'table name'

	s=size(lnr)
	my=s(2)
	print,'absorption ...'
	lnab=lnr
	for j=0,my-1 do lnab(*,j,*)=lookup(lnr(*,j,*),ee(*,j,*),iv=1)
	ab=exp(lnab)
	dlnab=fltarr(nx,my,nz)
	print,'dlnab ...'
	for j=1,my-2 do dlnab(*,j,*)=0.5*(lnab(*,j+1,*)-lnab(*,j-1,*))
	dlnab[*,0,*]=(lnab[*,1,*]-lnab[*,0,*])
	dlnab[*,my-1,*]=(lnab[*,my-1,*]-lnab[*,my-2,*])
	tau1=fltarr(nx,my,nz)
	tau1(*,0,*)=0.
	print,'tau ...'
	for j=1,my-1 do begin
	  dy=(ym(j)-ym(j-1))
	  dy12=dy/12.
	  ds=dy/xmu
	  tau1(*,j,*)=tau1(*,j-1,*)+ds* $
	    (ab(*,j-1,*)*((0.5+dlnab(*,j-1,*)*dy12))>0+ $
	     ab(*,j  ,*)*((0.5-dlnab(*,j  ,*)*dy12))>0)
	endfor
	return,tau1
end