program f90_SPL implicit none integer,parameter::NDmax=100 integer::nd,md,i,io real(8)::xdat(1:NDmax),ydat(1:NDmax) real(8),allocatable::datax(:),datay(:) character(len=50)::fnameR character(len=20)::dummy call getarg(1,fnameR) !Input data file call getarg(2,dummy) read(dummy,*) md open(11,file=fnameR,status='old') do i=1,NDmax read(11,*,iostat=io) xdat(i),ydat(i) if(io<0)exit end do close(11) nd=i-1 allocate(datax(1:nd),datay(1:nd)) do i=1,nd datax(i)=xdat(i) datay(i)=ydat(i) end do call SPL3(nd,datax,datay,md) stop contains subroutine SPL3(nd,datax,datay,md) integer,intent(in)::nd,md real(8),intent(in)::datax(1:nd),datay(1:nd) integer::i,j,nn,mm real(8)::a(1:nd),b(1:nd),c(1:nd),d(1:nd) real(8)::h0,h1,hj real(8)::dx,xx,sj real(8),allocatable::ary(:,:),vec(:) nn=nd-2 mm=2 allocate(ary(1:nn,1:mm),vec(1:nn)) do i=1,nn do j=1,mm ary(i,j)=0.0D0 end do vec(i)=0.0D0 end do do i=1,nd a(i)=datay(i) b(i)=0.0D0 c(i)=0.0D0 d(i)=0.0D0 end do do i=2,nd-1 h0=datax(i)-datax(i-1) h1=datax(i+1)-datax(i) ary(i-1,1)=2.0D0*(h0+h1) ary(i-1,2)=h1 vec(i-1)=3.0D0*(a(i+1)-a(i))/h1-3.0D0*(a(i)-a(i-1))/h0 end do call CHBAND(nn,mm,ary,vec) c(1)=0.0D0 do i=2,nd-1 c(i)=vec(i-1) end do c(nd)=0.0D0 do i=1,nd-1 hj=datax(i+1)-datax(i) d(i)=(c(i+1)-c(i))/(3.0D0*hj) b(i)=(a(i+1)-a(i))/hj-(2.0D0*c(i)+c(i+1))*hj/3.0D0 end do do i=1,nd-1 dx=(datax(i+1)-datax(i))/dble(md) do j=1,md xx=dble(j-1)*dx sj=a(i)+b(i)*xx+c(i)*xx**2.0D0+d(i)*xx**3.0D0 write(6,*) datax(i)+xx,sj end do end do i=nd-1 j=md+1 xx=dble(j-1)*dx sj=a(i)+b(i)*xx+c(i)*xx**2.0D0+d(i)*xx**3.0D0 write(6,*) datax(i)+xx,sj end subroutine SPL3 subroutine CHBAND(n1,m1,array,vector) integer,intent(in)::n1,m1 real(8),intent(inout)::array(1:n1,1:m1) real(8),intent(inout)::vector(1:n1) integer :: i,j,k,mm,mn real(8) :: z array(1,1)=sqrt(array(1,1)) do j=2,m1 array(1,j)=array(1,j)/array(1,1) end do do i=2,n1 if(n1-i+1>m1)then mm=m1 else mm=n1-i+1 end if do j=1,mm z=array(i,j) if(j/=m1)then if(m1-j+11)then array(i,j)=z/array(i,1) else array(i,1)=sqrt(z) end if else array(i,j)=z/array(i,1) end if end do end do vector(1)=vector(1)/array(1,1) do i=2,n1 z=vector(i) if(i>m1)then mm=m1 else mm=i end if do j=2,mm z=z-array(i-j+1,j)*vector(i-j+1) end do vector(i)=z/array(i,1) end do vector(n1)=vector(n1)/array(n1,1) do i=n1-1,1,-1 z=vector(i) do j=2,m1 if(i+j-1>n1)exit z=z-array(i,j)*vector(i+j-1) end do vector(i)=z/array(i,1) end do end subroutine CHBAND end program f90_SPL