subroutine splintl(xint,yout) * * interpolates for liquid * implicit double precision(a-h,o-z) common g(650),x(650),rho(650),yliq(4,650),ray(4,650) common/savingl/k,kq common/setup/isetup common/mfac/m common/yfacs/y1(650),y2(650),y3(650),y4(650),y5(650) common/rs/r(650) c common/savel/m,c1(4,650),y1(650),c2(4,650),y2(650), c 1 c3(4,650),y3(650),c4(4,650),y4(650),c5(4,650), c 2 y5(650),c6(4,650) real*8 yout(6),y6(650) real*8 c1(4,650),c2(4,650),c3(4,650),c4(4,650),c5(4,650) real*8 c6(4,650) * * compute spline coefficients for a model if we haven't done so * if (isetup.eq.1)then do 180 i=1,m y6(i)=r(i) 180 continue call consl(x,y1,m,c1) call consl(x,y2,m,c2) call consl(x,y3,m,c3) call consl(x,y4,m,c4) call consl(x,y5,m,c5) call consl(x,y6,m,c6) isetup=0 endif mm=m-1 if (xint.ge.x(1).and.xint.lt.x(m)) go to 40 if (xint .lt. (1.+1.e-8)*x(1)) go to 110 if (xint .ge. x(m)) go to 20 k=1 go to 70 20 if (xint .gt. (1.+1.e-6)*x(m)) go to 110 k=mm go to 70 40 il=1 ir=m 50 k=il+((ir-il)/2) if (xint.ge.x(k)) go to 60 ir=k go to 50 60 if (xint.lt.x(k+1)) go to 70 il=k go to 50 70 continue x1=x(k+1)-xint xx=xint-x(k) x12=x1*x1 xx2=xx*xx yout(1)=x1*(c1(1,k)*x12 +c1(3,k))+xx*(c1(2,k)*xx2 +c1(4,k)) yout(2)=x1*(c2(1,k)*x12 +c2(3,k))+xx*(c2(2,k)*xx2 +c2(4,k)) yout(3)=x1*(c3(1,k)*x12 +c3(3,k))+xx*(c3(2,k)*xx2 +c3(4,k)) yout(4)=x1*(c4(1,k)*x12 +c4(3,k))+xx*(c4(2,k)*xx2 +c4(4,k)) yout(5)=x1*(c5(1,k)*x12 +c5(3,k))+xx*(c5(2,k)*xx2 +c5(4,k)) yout(6)=x1*(c6(1,k)*x12 +c6(3,k))+xx*(c6(2,k)*xx2 +c6(4,k)) return 110 continue return 700 continue stop end ************************************************************************