UPP (develop)
Loading...
Searching...
No Matches
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