source:
LMDZ4/trunk/libf/dyn3dpar/spline.F
@
936
Last change on this file since 936 was 774, checked in by , 17 years ago | |
---|---|
|
|
File size: 1.6 KB |
Rev | Line | |
---|---|---|
[630] | 1 | ! |
2 | ! $Header$ | |
3 | ! | |
4 | subroutine spline(x,y,n,yp1,ypn,y2) | |
5 | ||
6 | c | |
7 | ||
8 | c Routine to set up the interpolating function for a cubic spline | |
9 | ||
10 | c interpolation (see "Numerical Recipes" for details). | |
11 | ||
12 | c | |
13 | implicit real (a-h,o-z) | |
14 | implicit integer (i-n) | |
15 | ||
16 | parameter(nllm=4096) | |
17 | ||
18 | dimension x(n),y(n),y2(n),u(nllm) | |
19 | ||
20 | c | |
21 | c write(6,*)(x(i),i=1,n) | |
22 | c write(6,*)(y(i),i=1,n) | |
23 | ||
24 | if(yp1.gt.0.99E30) then | |
25 | c the lower boundary condition is set | |
26 | y2(1)=0. | |
27 | c either to be "natural" | |
28 | u(1)=0. | |
29 | ||
30 | else | |
31 | c or else to have a specified first | |
32 | y2(1)=-0.5 | |
33 | c derivative | |
34 | u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1) | |
35 | ||
36 | end if | |
37 | ||
38 | do 11 i=2,n-1 | |
39 | c decomposition loop of the tridiagonal | |
40 | sig=(x(i)-x(i-1))/(x(i+1)-x(i-1)) | |
41 | c algorithm. Y2 and U are used | |
42 | p=sig*y2(i-1)+2. | |
43 | c for temporary storage of the decompo- | |
44 | y2(i)=(sig-1.)/p | |
45 | c sed factors | |
46 | u(i)=(6.*((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1)) | |
47 | ||
48 | . /(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*u(i-1))/p | |
49 | ||
50 | 11 continue | |
51 | ||
52 | if(ypn.gt.0.99E30) then | |
53 | c the upper boundary condition is set | |
54 | qn=0. | |
55 | c either to be "natural" | |
56 | un=0. | |
57 | ||
58 | else | |
59 | c or else to have a specified first | |
60 | qn=0.5 | |
61 | c derivative | |
62 | un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1))) | |
63 | ||
64 | end if | |
65 | ||
66 | y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.) | |
67 | ||
68 | do 12 k=n-1,1,-1 | |
69 | c this is the backsubstitution loop of | |
70 | y2(k)=y2(k)*y2(k+1)+u(k) | |
71 | c the tridiagonal algorithm | |
72 | 12 continue | |
73 | ||
74 | c | |
75 | ||
76 | return | |
77 | ||
78 | end | |
79 |
Note: See TracBrowser
for help on using the repository browser.