UPP  V11.0.0
 All Data Structures Files Functions Pages
SPLINE.f
1 !
2 !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
3  SUBROUTINE spline(JTB,NOLD,XOLD,YOLD,Y2,NNEW,XNEW,YNEW,P,Q)
4 ! ******************************************************************
5 ! * *
6 ! * THIS IS A ONE-DIMENSIONAL CUBIC SPLINE FITTING ROUTINE *
7 ! * PROGRAMED FOR A SMALL SCALAR MACHINE. *
8 ! * *
9 ! * PROGRAMER Z. JANJIC, YUGOSLAV FED. HYDROMET. INST., BEOGRAD *
10 ! * *
11 ! * *
12 ! * *
13 ! * NOLD - NUMBER OF GIVEN VALUES OF THE FUNCTION. MUST BE GE 3. *
14 ! * XOLD - LOCATIONS OF THE POINTS AT WHICH THE VALUES OF THE *
15 ! * FUNCTION ARE GIVEN. MUST BE IN ASCENDING ORDER. *
16 ! * YOLD - THE GIVEN VALUES OF THE FUNCTION AT THE POINTS XOLD. *
17 ! * Y2 - THE SECOND DERIVATIVES AT THE POINTS XOLD. IF NATURAL *
18 ! * SPLINE IS FITTED Y2(1)=0. AND Y2(NOLD)=0. MUST BE *
19 ! * SPECIFIED. *
20 ! * NNEW - NUMBER OF VALUES OF THE FUNCTION TO BE CALCULATED. *
21 ! * XNEW - LOCATIONS OF THE POINTS AT WHICH THE VALUES OF THE *
22 ! * FUNCTION ARE CALCULATED. XNEW(K) MUST BE GE XOLD(1) *
23 ! * AND LE XOLD(NOLD). *
24 ! * YNEW - THE VALUES OF THE FUNCTION TO BE CALCULATED. *
25 ! * P, Q - AUXILIARY VECTORS OF THE LENGTH NOLD-2. *
26 ! * *
27 ! ******************************************************************
28 !
29 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
30  implicit none
31 !
32  integer,intent(in) :: jtb,nold,nnew
33  real,dimension(JTB),intent(in) :: xold,yold,xnew
34  real,dimension(JTB),intent(inout) :: p,q,y2
35  real,dimension(JTB),intent(out) :: ynew
36 !
37  integer noldm1,k,k1,k2,kold
38  real dxl,dxr,dydxl,dydxr,rtdxc,dxc,den,xk,y2k,y2kp1,dx,rdx, &
39  ak,bk,ck,x,xsq
40 !-----------------------------------------------------------------------
41  noldm1=nold-1
42 !
43  dxl=xold(2)-xold(1)
44  dxr=xold(3)-xold(2)
45  dydxl=(yold(2)-yold(1))/dxl
46  dydxr=(yold(3)-yold(2))/dxr
47  rtdxc=.5/(dxl+dxr)
48 !
49  p(1)= rtdxc*(6.*(dydxr-dydxl)-dxl*y2(1))
50  q(1)=-rtdxc*dxr
51 !
52  IF(nold==3) go to 700
53 !-----------------------------------------------------------------------
54  k=3
55 !
56  100 dxl=dxr
57  dydxl=dydxr
58  dxr=xold(k+1)-xold(k)
59  dydxr=(yold(k+1)-yold(k))/dxr
60  dxc=dxl+dxr
61  den=1./(dxl*q(k-2)+dxc+dxc)
62 !
63  p(k-1)= den*(6.*(dydxr-dydxl)-dxl*p(k-2))
64  q(k-1)=-den*dxr
65 !
66  k=k+1
67  IF(k<nold) go to 100
68 !-----------------------------------------------------------------------
69  700 k=noldm1
70 !
71  200 y2(k)=p(k-1)+q(k-1)*y2(k+1)
72 !
73  k=k-1
74  IF(k>1) go to 200
75 !-----------------------------------------------------------------------
76  k1=1
77 !
78  300 xk=xnew(k1)
79 !
80  DO 400 k2=2,nold
81  IF(xold(k2)<=xk) go to 400
82  kold=k2-1
83  go to 450
84  400 CONTINUE
85  ynew(k1)=yold(nold)
86  go to 600
87 !
88  450 IF(k1==1) go to 500
89  IF(k==kold) go to 550
90 !
91  500 k=kold
92 !
93  y2k=y2(k)
94  y2kp1=y2(k+1)
95  dx=xold(k+1)-xold(k)
96  rdx=1./dx
97 !
98 !VVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVV
99 ! WRITE(6,5000) K,Y2K,Y2KP1,DX,RDX,YOLD(K),YOLD(K+1)
100 !5000 FORMAT(' K=',I4,' Y2K=',E12.4,' Y2KP1=',E12.4,' DX=',E12.4,' RDX='
101 ! 2,E12.4,' YOK=',E12.4,' YOP1=',E12.4)
102 !AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
103  ak=.1666667*rdx*(y2kp1-y2k)
104  bk=.5*y2k
105  ck=rdx*(yold(k+1)-yold(k))-.1666667*dx*(y2kp1+y2k+y2k)
106 !
107  550 x=xk-xold(k)
108  xsq=x*x
109 !
110  ynew(k1)=ak*xsq*x+bk*xsq+ck*x+yold(k)
111 !
112  600 k1=k1+1
113  IF(k1<=nnew) go to 300
114 !-----------------------------------------------------------------------
115  RETURN
116  END