1 2! Copyright (C) 2002-2005 J. K. Dewhurst, S. Sharma and C. Ambrosch-Draxl. 3! This file is distributed under the terms of the GNU Lesser General Public 4! License. See the file COPYING for license details. 5 6!BOP 7! !ROUTINE: rfinterp 8! !INTERFACE: 9subroutine rfinterp(ni,xi,wci,fi,no,xo,fo) 10! !INPUT/OUTPUT PARAMETERS: 11! ni : number of input points (in,integer) 12! xi : input abscissa array (in,real(ni)) 13! wci : input spline coefficient weights (in,real(12,ni)) 14! fi : input data array (in,real(ni)) 15! no : number of output points (in,integer) 16! xo : output abscissa array (in,real(no)) 17! fo : output interpolated function (out,real(no)) 18! !DESCRIPTION: 19! Given a function defined on a set of input points, this routine uses a 20! clamped cubic spline to interpolate the function on a different set of 21! points. See routine {\tt spline}. 22! 23! !REVISION HISTORY: 24! Created January 2005 (JKD) 25! Arguments changed, April 2016 (JKD) 26!EOP 27!BOC 28implicit none 29! arguments 30integer, intent(in) :: ni 31real(8), intent(in) :: xi(ni),wci(12,ni),fi(ni) 32integer, intent(in) :: no 33real(8), intent(in) :: xo(no) 34real(8), intent(out) :: fo(no) 35! local variables 36integer i,j,k,l 37real(8) x,dx 38! automatic arrays 39real(8) cf(3,ni) 40if (ni.le.0) then 41 write(*,*) 42 write(*,'("Error(rfinterp): invalid number of input points : ",I8)') ni 43 write(*,*) 44 stop 45end if 46if (no.le.0) then 47 write(*,*) 48 write(*,'("Error(rfinterp): invalid number of output points : ",I8)') no 49 write(*,*) 50 stop 51end if 52if (ni.eq.1) then 53 fo(:)=fi(1) 54 return 55end if 56! compute the spline coefficients 57call splinew(ni,wci,fi,cf) 58! evaluate spline at output points 59i=1 60do l=1,no 61 x=xo(l) 62 if (i.ge.ni) i=1 63 if (x.lt.xi(i)) goto 10 64 if (x.le.xi(i+1)) goto 30 65! binary search 6610 continue 67 i=1 68 j=ni 6920 continue 70 k=(i+j)/2 71 if (x.lt.xi(k)) then 72 j=k 73 else 74 i=k 75 end if 76 if (j.gt.i+1) goto 20 7730 continue 78 dx=x-xi(i) 79 fo(l)=fi(i)+dx*(cf(1,i)+dx*(cf(2,i)+dx*cf(3,i))) 80end do 81end subroutine 82!EOC 83 84