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 . INTENT(IN) :: xi
57 . DIMENSION(*), INTENT(IN) :: kxi
59
60
61
62 INTEGER J, JJ, K, NDERS
63 my_real saved, temp, aleft, right
65 my_real,
DIMENSION(PXI+1,PXI+1) :: andu
67
68 nders=1
69 andu(:,:)=zero
70
71
72 DO j=0,pxi
73 IF ((xi>=kxi(idxi+j)).AND.(xi<kxi(idxi+j+1))) THEN
74 andu(j+1,1) = one
75 ELSE
76 andu(j+1,1) = zero
77 ENDIF
78 ENDDO
79
80 DO k=1,pxi
81 IF (andu(1,k) == 0) THEN
82 saved = zero
83 ELSE
84 saved = ((xi-kxi(idxi))*andu(1,k))/(kxi(idxi+k)-kxi(idxi))
85 ENDIF
86 DO j=0,pxi-k
87 aleft = kxi(idxi+j+1)
88 right = kxi(idxi+j+k+1)
89 IF (andu(j+2,k) == 0) THEN
90 andu(j+1,k+1) = saved
91 saved = zero
92 ELSE
93 temp = andu(j+2,k)/(right-aleft)
94 andu(j+1,k+1) = saved+(right-xi)*temp
95 saved = (xi-aleft)*temp
96 ENDIF
97 ENDDO
98 ENDDO
99
100 ders(1) = andu(1,pxi+1)
101
102 DO k=1,nders
103 DO j=1,k+1
104 nd(j) = andu(j,pxi-k+1)
105 ENDDO
106 DO jj=1,k
107 IF (nd(1) == 0) THEN
108 saved = zero
109 ELSE
110 saved = nd(1)/(kxi(idxi+pxi-k+jj)-kxi(idxi))
111 ENDIF
112 DO j=1,k-jj+1
113 aleft = kxi(idxi+j)
114 right = kxi(idxi+j+pxi+jj-1)
115
116 IF (nd(j+1) == 0) THEN
117 nd(j) = (pxi-k+jj)*saved
118 saved = zero
119 ELSE
120 temp = nd(j+1)/(right-aleft)
121 nd(j) = (pxi-k+jj)*(saved-temp)
122 saved = temp
123 ENDIF
124 ENDDO
125 ENDDO
126 ders(2) = nd(1)
127 ENDDO
128
129 ders1 = ders(1)
130 ders2 = ders(2)
131
132 RETURN