V3FIT
tok_fraction.f
1  SUBROUTINE tok_fraction(ft_tok)
2 c--
3 c This calculates the fraction passing and the fraction trapped
4 c for a tokamak using the Lin-Liu Miller approximation with Houlberg
5 c weighting factors. (Phys. Plas. 2, 1966 (1995).
6 c--
7 C-----------------------------------------------
8 C M o d u l e s
9 C-----------------------------------------------
10  USE parambs
11  USE vmec0
12  IMPLICIT NONE
13 C-----------------------------------------------
14 C D u m m y A r g u m e n t s
15 C-----------------------------------------------
16  REAL(rprec) :: ft_tok
17 C-----------------------------------------------
18 C L o c a l P a r a m e t e r s
19 C-----------------------------------------------
20  REAL(rprec), PARAMETER :: one = 1
21 C-----------------------------------------------
22 C L o c a l V a r i a b l e s
23 C-----------------------------------------------
24  INTEGER i
25  REAL(rprec) fub, flb, bavg, b2av, sum_gsqrt_tok
26  REAL(rprec) bmax
27 C-----------------------------------------------
28 
29  REAL(rprec), DIMENSION(:), ALLOCATABLE ::
30  1 bfield_tok, gsqrt_tok, b2obm_tok, one_b2
31  ALLOCATE(bfield_tok(nthetah), gsqrt_tok(nthetah),
32  1 b2obm_tok(nthetah), one_b2(nthetah), stat = i)
33  IF(i .ne. 0) stop 'allocation error of tokamak fields'
34 
35 
36 c make symetric copies of field quantites--average over zeta to extract
37 c symmetric component.
38 
39  DO i =1,nthetah
40  bfield_tok(i) = sum(bfield(i,:nzetah))/nzetah
41  ENDDO
42 
43 c renormalize to MAX tok field equivalent
44 
45  bmax = maxval(bfield_tok)
46  bfield_tok = bfield_tok/bmax
47  WHERE(bfield_tok .gt. one) bfield_tok = one
48 
49 c calculate 1/b**2, b**2
50 
51  one_b2 = one/bfield_tok**2
52  b2obm_tok = bfield_tok**2
53 
54 c jacobian only includes 1/b**2 component and is normalized to bmax tok
55 c integrate over jacobian to obtain normalization for averages
56  sum_gsqrt_tok= sum(one_b2)
57 
58 
59 c find average b, b**2
60 
61  bavg = sum(bfield_tok*one_b2)/sum_gsqrt_tok
62  b2av = sum(b2obm_tok*one_b2)/sum_gsqrt_tok
63 
64 c find upper bound
65 
66  fub = one-(one-sqrt(one-bavg)*(one+0.5_dp*bavg))*b2av/bavg**2
67 
68 c find lower bound
69 
70 
71 c minus <1/b**2>
72 
73  flb = - sum(one/b2obm_tok**2)/sum_gsqrt_tok
74 
75 c plus <SQRT(1-b)/b**2>
76 
77  flb = flb + sum(sqrt(one-bfield_tok)
78  1 /b2obm_tok**2)/sum_gsqrt_tok
79 
80 c plus <SQRT(1-b)/b)>/2
81 
82  flb = flb + 0.5_dp*sum(sqrt(one-bfield_tok)
83  1 /bfield_tok/b2obm_tok)/sum_gsqrt_tok
84 
85 c 1+(previous stuff)*<b**2>
86 
87 
88  flb = one + flb*b2av
89 
90 c finally combine upper and lower limits with Houlberg factors
91  ft_tok = 0.25_dp*flb + 0.75_dp*fub
92 
93  DEALLOCATE(bfield_tok, gsqrt_tok, b2obm_tok, one_b2)
94 
95  RETURN
96  END SUBROUTINE tok_fraction