1 SUBROUTINE spline_akima(x,y,xx,yy,npts,iflag)
8 integer,
intent(in) :: npts
9 integer,
intent(inout) :: iflag
10 real(rprec),
intent(in) :: x
11 real(rprec),
intent(out) :: y
12 real(rprec),
dimension(npts),
intent(in) :: xx, yy
16 real(rprec),
dimension(-1:size(xx)+2) :: a, b, c, d
17 real(rprec),
dimension(-1:size(xx)+2) :: xloc, yloc
18 real(rprec),
dimension(-1:size(xx)+2) :: m, t, dm, p, q
19 integer :: i, ix, iv, ivm1
20 real(rprec) :: cl,bl,cr,br,dx
34 > stop
'spline_akima: more active points requested than available'
35 if(ix /=
size(yy)) stop
'size mismatch of xx and yy!'
39 a = 0._dp ; b = 0._dp ; c = 0._dp ; d = 0._dp
40 xloc = 0._dp ; yloc = 0._dp
41 m = 0._dp ; t = 0._dp ; dm = 0._dp
45 xloc(-1)= 2*xloc(1)-xloc(3)
46 xloc( 0)= xloc(1)+xloc(2)-xloc(3)
47 xloc(iv+2)= 2*xloc(iv)-xloc(iv-2)
48 xloc(iv+1)= xloc(iv)+xloc(iv-1)-xloc(iv-2)
51 m(1:iv-1) = (yloc(2:iv)-yloc(1:iv-1))/
52 > (xloc(2:iv)-xloc(1:iv-1))
54 cl = (m(2)-m(1))/(xloc(3)-xloc(1))
55 bl = m(1) - cl*(xloc(2)-xloc(1))
56 cr = (m(iv-2)-m(iv-1))/(xloc(iv)-xloc(iv-2))
57 br = m(iv-2) - cr*(xloc(iv-1)-xloc(iv-2))
58 yloc( 0)=yloc(1)+bl*(xloc( 0)-xloc(1))+
59 > cl*(xloc( 0)-xloc(1))**2
60 yloc(-1)=yloc(1)+bl*(xloc(-1)-xloc(1))+
61 > cl*(xloc(-1)-xloc(1))**2
62 yloc(iv+1)=yloc(iv)+br*(xloc(iv+1)-xloc(iv))+
63 > cr*(xloc(iv+1)-xloc(iv))**2
64 yloc(iv+2)=yloc(iv)+br*(xloc(iv+2)-xloc(iv))+
65 > cr*(xloc(iv+2)-xloc(iv))**2
67 m(-1) = (yloc(0)-yloc(-1))/(xloc(0)-xloc(-1))
68 m( 0) = (yloc(1)-yloc( 0))/(xloc(1)-xloc( 0))
69 m(iv ) = (yloc(iv+1)-yloc(iv ))/(xloc(iv+1)-xloc(iv ))
70 m(iv+1) = (yloc(iv+2)-yloc(iv+1))/(xloc(iv+2)-xloc(iv+1))
72 dm(-1:iv)= abs(m(0:iv+1)-m(-1:iv))
74 p(1:iv) = dm(1:iv)/(dm(1:iv)+dm(-1:iv-2))
77 q(1:iv) = dm(-1:iv-2)/(dm(1:iv)+dm(-1:iv-2))
79 t(1:iv) = p(1:iv)*m(0:iv-1)+q(1:iv)*m(1:iv)
80 where ( p(1:iv)+q(1:iv) < tiny(1._dp))
81 t(1:iv) = 0.5_dp*m(0:iv-1)+0.5_dp*m(1:iv)
86 c(1:iv-1) = (3*m(1:iv-1)-t(2:iv)-2*t(1:iv-1))/
87 > (xloc(2:iv)-xloc(1:iv-1))
88 d(1:iv-1) = (t(2:iv)+t(1:iv-1)-2*m(1:iv-1))/
89 > (xloc(2:iv)-xloc(1:iv-1))**2
92 if(x<xloc(1) .or. x>xloc(iv))
then
97 if(x == xloc(iv))
then
103 if(x >= xloc(i) .and. x < xloc(i+1))
then
105 y=a(i)+dx*(b(i)+dx*(c(i)+d(i)*dx))
111 END SUBROUTINE spline_akima