1 SUBROUTINE add_tension(amat, wten, hx, tens, tensv, fpoly,
8 INTEGER n, nb, ioff, nmat
9 REAL(rprec) tens, tensv, fpoly
10 REAL(rprec),
DIMENSION(nmat,*) :: amat
11 REAL(rprec),
DIMENSION(nmat) :: wten
12 REAL(rprec),
DIMENSION(n) :: hx
17 REAL(rprec),
DIMENSION(n) :: wten0, tenshx, work
18 REAL(rprec) :: delta_x, tension
37 IF (n + ioff > nmat) stop
'(n+ioff>nmat)'
38 IF (fpoly < 0.) stop
'(fpoly<0)'
39 IF (n < 1) stop
'(n < 1)'
45 delta_x = sum(hx(:n-1))
46 tension = 0.5*(delta_x/(n))**3
47 IF (fpoly.eq.0. .or. tens.eq.tensv)
THEN
51 tenshx(i) = tens + (tensv - tens)*(real(i - 1,rprec)/
57 tenshx(i) = tension * tenshx(i) / hx(i)
58 work(i) = hx(i)*(wten(i+ioff) + wten(i+ioff+1))
61 wten0(i) = 0.5 * ( work(i) + work(i-1) )/(hx(i) + hx(i-1))
64 wten0(n) = wten(n+ioff)
75 work(2) = tenshx(1)*wten0(2)
77 ELSE IF (k .eq. n)
THEN
78 work(n-1) = tenshx(n-1)*wten0(n-1)
80 work(k-1) = tenshx(k-1)*wten0(k-1)
81 work(k) = -(tenshx(k)+tenshx(k-1))*wten0(k)
82 work(k+1) = tenshx(k)*wten0(k+1)
84 IF (nb .eq. natur) work(1) = 0
89 CALL jacprod (work, hx, n, nb)
90 amat(koff,1+ioff:n+ioff) = amat(koff,1+ioff:n+ioff) + work(:n)
93 END SUBROUTINE add_tension