source: LMDZ4/branches/unlabeled-1.1.1/libf/dyn3d/spline.F @ 4667

Last change on this file since 4667 was 524, checked in by lmdzadmin, 20 years ago

Initial revision

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 1.6 KB
Line 
1!
2! $Header$
3!
4      subroutine spline(x,y,n,yp1,ypn,y2)
5     
6c
7     
8c     Routine to set up the interpolating function for a cubic spline
9     
10c     interpolation (see "Numerical Recipes" for details).
11     
12c
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     
20c
21c       write(6,*)(x(i),i=1,n)
22c       write(6,*)(y(i),i=1,n)
23     
24      if(yp1.gt.0.99E30) then
25c the lower boundary condition is set
26       y2(1)=0.
27c either to be "natural"
28       u(1)=0.
29     
30      else
31c or else to have a specified first
32       y2(1)=-0.5
33c 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
39c decomposition loop of the tridiagonal
40       sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
41c algorithm. Y2 and U are used
42       p=sig*y2(i-1)+2.
43c for temporary storage of the decompo-
44       y2(i)=(sig-1.)/p
45c 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
53c the upper boundary condition is set
54       qn=0.
55c either to be "natural"
56       un=0.
57     
58      else
59c or else to have a specified first
60       qn=0.5
61c 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
69c this is the backsubstitution loop of
70       y2(k)=y2(k)*y2(k+1)+u(k)
71c the tridiagonal algorithm
72 12   continue
73     
74c
75     
76      return
77     
78      end
79     
Note: See TracBrowser for help on using the repository browser.