spline interpolation function for advanced fortran

This subroutine will help you to spline the whole dimension of a vector. This can be easily modified to multi-dimensional matrix calculation.


subroutine spline(n,x,y,yp1,ypn,no,xo,yo)
cccccccccccccccccccccccccccccccccccccccccccccccccccc
c This subroutine will calculate the cubic-spline intepolated value of
c given any function for single variable of any dimension
c
c inputs/output:
c n : (input) [scalar] dimension of input array
c x : (input) [vector of length 'n'], independent variable of dimension 'n'
c y : (input)[vector of length 'n']' the function of x (an array of dimension 'n'
c yp1 : (input) [scalar] =1.e+30 :: the routine is signaled to set
c the corresponding boundary condition
c for a natural spline, with zero
c second derivative on that boundary
c ypn : (input) same to yp1 = 1.e+30
c no : (input) [scalar] dimension output array
c xo : (input) [vector of length 'no'], 'x' values at which you want to spline
c yo : (output) [vector of length 'no'], splined value
c
c Note: This works properly in any inter fortran compiler and f90-compiler
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c Tanmoy Das.
c
c Feb 13, 2008.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
integer::n,no,io
real,dimension(1:n)::x,y,y2,u
real,dimension(1:no)::xo,yo
c
if (yp1.gt.0.99e30) then
y2(1)=0.
u(1)=0.
else
y2(1)=-0.5
u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
endif
do 11 i=2,n-1
sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
p=sig*y2(i-1)+2.
y2(i)=(sig-1.)/p
u(i)=(6.*((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1))
& /(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*u(i-1))/p
11 continue
if (ypn.gt.0.99e30) then
qn=0.
un=0.
else
qn=0.5
un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
endif
y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.)
do 12 k=n-1,1,-1
y2(k)=y2(k)*y2(k+1)+u(k)
12 continue
cccccccccccccccccccccccccccccccccccccccccccccc
c now for any arbitrary x
cccccccccccccccccccccccccccccccccccccccccccccc
do io = 1,no
klo=1
khi=n
1 if (khi-klo.gt.1) then
k=(khi+klo)/2.
if (x(k).gt.xo(io)) then
khi=k
else
klo=k
endif
goto 1
endif
h=x(khi)-x(klo)
if (h.eq.0.) pause 'bad XA input'
a=(x(khi)-xo(io))/h
b=(xo(io)-x(klo))/h
yo(io)=a*y(klo)+b*y(khi)+
& (a*(a*a-1.)*y2(klo)+b*(b*b-1.)*y2(khi))*h*h/6.

enddo
c
return
end