source: LMDZ.3.3/trunk/libf/dyn3d/spline.F @ 5306

Last change on this file since 5306 was 2, checked in by lmdz, 25 years ago

Initial revision

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