V3FIT
add_tension.f
1  SUBROUTINE add_tension(amat, wten, hx, tens, tensv, fpoly,
2  1 n, nb, ioff, nmat)
3  USE vspline
4  IMPLICIT NONE
5 C-----------------------------------------------
6 C D u m m y A r g u m e n t s
7 C-----------------------------------------------
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
13 C-----------------------------------------------
14 C L o c a l V a r i a b l e s
15 C-----------------------------------------------
16  INTEGER :: k, i, koff
17  REAL(rprec), DIMENSION(n) :: wten0, tenshx, work
18  REAL(rprec) :: delta_x, tension
19 C-----------------------------------------------
20 
21 !
22 ! THIS SUBROUTINE ADDS THE TENSION TERM WTEN*(DEL(K) - DEL(K-1))
23 ! TO THE REST OF THE CHI-SQ AMAT, WHERE DEL(K) = (G(K+1)-G(K))/H(K)
24 ! AND G(K) IS THE SECOND DERIVATIVE AT THE K-TH KNOT
25 !
26 ! IOFF ALLOWS FOR ADDING TENSION SEPARATELY TO IOTA (IOFF=0)
27 ! AND PRESSURE (IOFF=NIOTA) SPLINES
28 !
29 ! tens: spline tension (optionally, at the 1st pt only;see note)
30 ! tensv: vbl spline tension for n-1th point (see note)
31 ! fpoly: vbl spline tension form factor (note: if tens1<>tens2
32 ! then tension(i-th point) = tens+(tensv-tens)*(i/n-1))**fpoly)
33 
34 !
35 ! BOUNDS CHECKING
36 !
37  IF (n + ioff > nmat) stop '(n+ioff>nmat)'
38  IF (fpoly < 0.) stop '(fpoly<0)'
39  IF (n < 1) stop '(n < 1)'
40 
41 !
42 ! COMPUTE TENSION COEFFICIENTS
43 !
44 
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
48  tenshx(:n-1) = tens
49  ELSE
50  DO i = 1, n - 1
51  tenshx(i) = tens + (tensv - tens)*(real(i - 1,rprec)/
52  1 (n - 1))**fpoly
53  END DO
54  ENDIF
55 
56  DO i = 1,n-1
57  tenshx(i) = tension * tenshx(i) / hx(i)
58  work(i) = hx(i)*(wten(i+ioff) + wten(i+ioff+1))
59  ENDDO
60  DO i = 2,n-1
61  wten0(i) = 0.5 * ( work(i) + work(i-1) )/(hx(i) + hx(i-1))
62  ENDDO
63  wten0(1) = wten0(2)
64  wten0(n) = wten(n+ioff)
65 !
66 ! COMPUTE, FOR K = 1,N, B(K,L)*JACOBIAN(L,I) = W(K,I),
67 ! WHERE JACOBIAN = D[G]/D[F] and B is TRIDIAGONAL
68 ! SEE EQN(27) IN PHYS.PLASMAS 1, p 2277.
69 !
70  DO k = 1, n
71  koff = k + ioff
72  work(:n-1) = 0
73 ! SET UP COEFFICIENTS IN [G(K+1)-G(K)]/h(k) - [G(K)-G(K-1)]/h(k-1)
74  IF (k .eq. 1) THEN
75  work(2) = tenshx(1)*wten0(2)
76  work(1) = -work(2)
77  ELSE IF (k .eq. n) THEN
78  work(n-1) = tenshx(n-1)*wten0(n-1)
79  ELSE
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)
83  ENDIF
84  IF (nb .eq. natur) work(1) = 0
85  work(n) = 0
86 !
87 ! COMPUTE work(j) = work(i)*Jacobian(i,j) and add to amat(k,j)
88 !
89  CALL jacprod (work, hx, n, nb)
90  amat(koff,1+ioff:n+ioff) = amat(koff,1+ioff:n+ioff) + work(:n)
91  END DO
92 
93  END SUBROUTINE add_tension