29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49#include "implicit_f.inc"
50
51
52
53 INTEGER PXI, IDXI
55 my_real,
DIMENSION(*),
INTENT(IN) :: kxi
56 my_real,
DIMENSION(*),
INTENT(OUT) :: ders1, ders2
57
58
59
60 INTEGER J, K, L, KR, KP, J1, J2, NDERS, LS1, LS2
62 my_real,
DIMENSION(PXI+1) :: aleft, right
63 my_real,
DIMENSION(2,PXI+1) :: ders, a
64 my_real,
DIMENSION(PXI+1,PXI+1) :: andu
65
66 nders=1
67 andu(1,1)=one
68
69 DO j = 1,pxi
70 aleft(j+1) = xi - kxi(idxi+1-j)
71 right(j+1) = kxi(idxi+j) - xi
72 saved = zero
73 DO l = 0,j-1
74 andu(j+1,l+1) = right(l+2) + aleft(j-l+1)
75 temp = andu(l+1,j)/andu(j+1,l+1)
76 andu(l+1,j+1) = saved + right(l+2)*temp
77 saved = aleft(j-l+1)*temp
78 ENDDO
79 andu(j+1,j+1) = saved
80 ENDDO
81
82
83 DO j = 0,pxi
84 ders(1,j+1) = andu(j+1,pxi+1)
85 ENDDO
86
87
88 DO l = 0,pxi
89 ls1 = 0
90 ls2 = 1
91 a(1,1) = one
92
93 DO k = 1,nders
94 d = zero
95 kr = l-k
96 kp = pxi-k
97 IF (l >= k) THEN
98 a(ls2+1,1) = a(ls1+1,1)/andu(kp+2,kr+1)
99 d = a(ls2+1,1)*andu(kr+1,kp+1)
100 ENDIF
101 IF (kr >= -1) THEN
102 j1 = 1
103 ELSE
104 j1 = -kr
105 ENDIF
106 IF ((l-1) <= kp) THEN
107 j2 = k-1
108 ELSE
109 j2 = pxi-l
110 ENDIF
111 DO j = j1,j2
112 a(ls2+1,j+1) = (a(ls1+1,j+1) - a(ls1+1,j))/andu(kp+2,kr+j+1)
113 d = d + a(ls2+1,j+1)*andu(kr+j+1,kp+1)
114 ENDDO
115 IF (l <= kp) THEN
116 a(ls2+1,k+1) = -a(ls1+1,k)/andu(kp+2,l+1)
117 d = d + a(ls2+1,k+1)*andu(l+1,kp+1)
118 ENDIF
119 ders(k+1,l+1) = d
120 j = ls1
121 ls1 = ls2
122 ls2 = j
123 ENDDO
124 ENDDO
125
126
127
128 l = pxi
129 DO k = 1,nders
130 DO j = 0,pxi
131 ders(k+1,j+1) = ders(k+1,j+1)*l
132 ENDDO
133 l = l*(pxi-k)
134 ENDDO
135
136 DO j = 1,pxi+1
137 ders1(j) = ders(1,j)
138 ders2(j) = ders(2,j)
139 ENDDO
140
141 RETURN