	function vol(tl,th,cspt,ntpt,tpt)
c     returns integral between tl and th using trapezoidal rule on
c     ntpt points at tpt with values cspt
c
	dimension cspt(1),tpt(1)
c
	ilt=0
	iht=0
	do 1 iv=1,ntpt
	if (tl.gt.tpt(iv)) ilt=iv
	if (th.gt.tpt(iv)) iht=iv
1		continue
	if (ilt.eq.0) goto 4
	if (ilt.eq.ntpt) goto 7
	ilt=ilt+1
c     t bins ilt to iht are inside t limits
	if (iht.lt.ilt) goto 8
	vol=.5*(tpt(ilt)-tl)*(cspt(ilt)+((tpt(ilt)-tl)*cspt(ilt-1)+
& (tl-tpt(ilt-1))*cspt(ilt))/(tpt(ilt)-tpt(ilt-1)))
6		if (iht.eq.ntpt) goto 10
	vol=vol+.5*(th-tpt(iht))*(cspt(iht)+((tpt(iht+1)-th)*cspt(iht)+
& (th-tpt(iht))*cspt(iht+1))/(tpt(iht+1)-tpt(iht)))
3		iht1=iht-1
	if (iht1.lt.ilt) return
	do 2 i=ilt,iht1
2		vol=vol+.5*(tpt(i+1)-tpt(i))*(cspt(i)+cspt(i+1))
	return
4		if (iht.gt.0) goto 5
c     extrapolation using first 2 points
	vol=.5*(th-tl)*((tpt(2)-tl)*cspt(1)+(tl-tpt(1))*cspt(2)+
& (tpt(2)-th)*cspt(1)+(th-tpt(1))*cspt(2))/(tpt(2)-tpt(1))
	return
c
c     extrapolation up to tpt(1)
5		vol=.5*(tpt(1)-tl)*(cspt(1)+((tpt(2)-tl)*cspt(1)+(tl-tpt(1))
& *cspt(2))/(tpt(2)-tpt(1)))
	ilt=ilt+1
	goto 6
c
c     extrapolation using last 2 points
7		vol=.5*(th-tl)*((tpt(ntpt)-tl)*cspt(ntpt-1)+(tl-tpt(ntpt-1))
& *cspt(ntpt)+(tpt(ntpt)-th)*cspt(ntpt-1)+(th-tpt(ntpt-1))
& *cspt(ntpt))/(tpt(ntpt)-tpt(ntpt-1))
	return
c
c     tl and th are between iht(lower) and ilt(higher) by 1
8		vol=.5*(th-tl)*((tl-tpt(iht))*cspt(ilt)+(tpt(ilt)-tl)*cspt(iht)+
& (th-tpt(iht))*cspt(ilt)+(tpt(ilt)-th)*cspt(iht))/
& (tpt(ilt)-tpt(iht))
	return
10		vol=vol+.5*(th-tpt(iht))*(cspt(iht)+((tpt(iht)-th)*cspt(iht-1)
& +(th-tpt(iht-1))*cspt(iht))/(tpt(iht)-tpt(iht-1)))
	goto 3
	end
