V3FIT
woflam.f
1  SUBROUTINE woflam(trigsu, trigsv, ifaxu, ifaxv, irho)
2 
3 c evaluate the lambda-dependent part of the functions w1(lambda)
4 c for all lambda.
5 c Calculates the SUM over n and m of (mR+nS)/(m-nq) * dalphamn/dlambda
6 c * <V||/V (mn)>/<V||/V>*lambda*(-2).
7 c Here, m is the poloidal mode number and n is the toroidal mode
8 c number/periods.
9 C-----------------------------------------------
10 C M o d u l e s
11 C-----------------------------------------------
12  USE parambs
13  USE vmec0
14  IMPLICIT NONE
15 C-----------------------------------------------
16 C D u m m y A r g u m e n t s
17 C-----------------------------------------------
18  INTEGER :: irho
19  INTEGER, DIMENSION(13) :: ifaxu, ifaxv
20  REAL(rprec), DIMENSION(3*nthetah/2 + 1) :: trigsu
21  REAL(rprec), DIMENSION(2*nzetah) :: trigsv
22 C-----------------------------------------------
23 C L o c a l P a r a m e t e r s
24 C-----------------------------------------------
25  REAL(rprec), PARAMETER :: zero = 0, one = 1,
26  1 d18 = 1.0e-18_dp, xlam = 0.96_dp
27  INTEGER, PARAMETER :: n_lam_coarse = 97, n_lam = 137
28 
29 c as discussed below, the mesh is split with 96 + 40 intervals
30 
31 C-----------------------------------------------
32 C L o c a l V a r i a b l e s
33 C-----------------------------------------------
34  INTEGER :: l, n, m
35  REAL(rprec), DIMENSION(:,:), ALLOCATABLE :: a11
36  complex(rprec), DIMENSION(:,:), ALLOCATABLE ::
37  1 alphamn, vmn
38  REAL(rprec) :: qn, numer, avg_vpov, denom
39  REAL(rprec), DIMENSION(n_lam) :: xlam_f
40  REAL(rprec), DIMENSION(n_lam-1) :: xlam_h
41 C-----------------------------------------------
42 C E x t e r n a l F u n c t i o n s
43 C-----------------------------------------------
44  REAL(rprec) , EXTERNAL :: sumit
45 C-----------------------------------------------
46 c
47 c
48  ALLOCATE (a11(nthetah+2,nzetah),
49  5 alphamn(-mbuse:mbuse,0:nbuse),
50  8 vmn(-mbuse:mbuse,0:nbuse), stat=l)
51 
52  IF (l .ne. 0) stop 'allocation error in woflam'
53 
54 c-----------------------------------------------------------------------
55 c///////////////////////////////////////////////////////////////////////
56 c-----------------------------------------------------------------------
57 
58 c meshing for lamda integral. the lambda integral has a near singular
59 c value at lambda = 1. to evaluate, construct dalpha/dlamda on full
60 c mesh. evaluate the rest on the full mesh. becasue of the (almost)
61 c singularity, split the mesh into coarse (0 to 0.96) and fine (.96 to
62 c 1.0) components. this was tested on both qas and qos configuratons
63 c and gives results that are a percent or better on the integral and
64 c ten to twenty times better for the total current. to save storage,
65 c there is only one loop on lambda. all grids are local.
66 
67 c-----------------------------------------------------------------------
68 c///////////////////////////////////////////////////////////////////////
69 c-----------------------------------------------------------------------
70 
71 
72  IF (isymm0 .ne. 0) THEN
73  aiterm1(irho) = zero
74  RETURN
75  ENDIF
76 
77 c First contruct the meshes. This scheme is hardwired.
78 
79 
80  xlam_f(1) = zero
81  xlam_h(1) = 0.005_dp
82 
83  DO l = 2, n_lam_coarse-1
84  xlam_f(l) = xlam_f(l-1) + 0.01_dp
85  xlam_h(l) = xlam_h(l-1) + 0.01_dp
86  ENDDO
87 
88  xlam_f(n_lam_coarse) = 0.96_dp !!End coarse mesh/Begin fine mesh
89  xlam_h(n_lam_coarse) = 0.9605_dp !!Begin of fine half-mesh
90 
91  DO l = n_lam_coarse+1, n_lam-1
92  xlam_f(l) = xlam_f(l-1) + 0.001_dp
93  xlam_h(l) = xlam_h(l-1) + 0.001_dp
94  ENDDO
95 
96  xlam_f(n_lam) = one
97 
98 c Loop over LAMBDA = 0 to 1, calculating alpha on the full mesh, and
99 c the rest of the integral on the half mesh. The jacobian for alpha
100 c puts in a b2avg but one of the b values in the denominator cancels the
101 c b/bmax in alpha. This form of the expression paralles that in the
102 c hamada paper but with flux surface averages for the alpha. This results
103 c in a form that does not have the beta found in the boozer paper
104 c see refereces at beginning.
105 
106  DO l = 1, n_lam
107  a11(:nthetah,:nzetah) =
108  2 sqrt(abs(one - xlam_f(l)*bfield(:nthetah,:nzetah))+d18)
109  3 *b2avg(irho)/bmax1(irho)**2/bfield(:nthetah,:nzetah)
110  a11(nthetah+1,:nzetah) = zero
111  a11(nthetah+2,:nzetah) = zero
112  CALL do_fft (a11, fmn, trigsu, trigsv, ifaxu, ifaxv, nthetah,
113  1 nzetah, mbuse, nbuse)
114  IF(l .eq. 1) THEN
115  alphamn = fmn
116  cycle !need to calculate to alphas and difference
117  END IF !before a term for the integral can be evaluated.
118 
119 
120 c SAVE alpha at lambda = 1
121 
122  IF(l .eq. n_lam) alpha1mn = fmn
123 
124 c form d_alphalmn at this point we have the value of alpha for l in fmn
125 c and the previous value of alpha (l-1) in alphamn. Store the difference in
126 c fmn using vmn as a temp variable to hold alpha, THEN update alpha for
127 c the next CYCLE.
128 
129  vmn = fmn
130  fmn = fmn - alphamn
131  alphamn = vmn
132 
133 c- Now DO the fft to get <EXP(i(m*theta-n*zeta)V||/V> on the half mesh
134 
135  a11(:nthetah,:nzetah) =
136  1 sqrt(abs(one - xlam_h(l-1)*bfield(:nthetah,:nzetah)) + d18)
137  2 *(b2avg(irho)/bfield(:nthetah,:nzetah)*bmax1(irho))**2
138 
139  a11(nthetah+1,:nzetah) = zero
140  a11(nthetah+2,:nzetah) = zero
141 
142  CALL do_fft (a11, vmn, trigsu, trigsv, ifaxu, ifaxv, nthetah,
143  1 nzetah, mbuse, nbuse)
144 
145 c This needs to be divided by <V||/V>. But this is exactly the REAL part of the
146 c 0,0 term of this transform.
147  avg_vpov = real(vmn(0,0))
148 
149 
150 c- for only those harmonics that are going to be used in the sum.
151 
152  qn = periods*qsafety(irho)*zetasign
153  rfmn(0,0) = zero
154  DO m = -mbuse, mbuse
155  DO n = 0, nbuse
156 
157  denom = m + n*qn
158  IF (n.ne.0 .or. m.ne.0) THEN
159 
160  numer = denom/(denom**2 + (damp_bs*m)**2)
161 
162  rfmn(m,n) = (m*capr(irho) + n*periods*caps(irho))*
163  1 real(fmn(m,n)*vmn(m,n))*(-2*numer)
164  ENDIF
165  END DO
166  END DO
167 c
168 c- the sum in w is the sum over all n and m (excluding (0,0)) of rfmn.
169 c but only the nonegative m have been calculated and stored. therefore
170 c the sum becomes
171 c 2 * sum over m = -mbuse to mbuse and n = 1 to nbuse of rfmn
172 c plus sum over n = -mbuse to mbuse of fn0
173 c minus f(0,0)
174 
175  w1(l-1) = xlam_h(l-1)*sumit(rfmn,mbuse,nbuse)/avg_vpov
176 
177 c- End loop over lambda.
178 
179  END DO
180 c evaluate integral
181 
182  aiterm1(irho) = sum(w1)
183  aiterm1(irho) = -0.75_dp*aiterm1(irho)*qsafety(irho)
184  1 /ftrapped(irho) * (one + aiogar(irho)/qsafety(irho))
185 
186  DEALLOCATE (a11, alphamn, vmn, stat=l)
187 
188  END SUBROUTINE woflam