V3FIT
jxbforce.f
1  SUBROUTINE jxbforce(bsupu, bsupv, bsubu, bsubv, bsubs, bsubsu,
2  1 bsubsv, gsqrt, bsq, itheta, izeta, brho, sigma_an, ier_flag
3 #ifdef _ANIMEC
4  2 ,pp1, pp2, ppar, onembc
5 #endif
6  3 )
7  USE safe_open_mod
8  USE vmec_main
9  USE vmec_params, ONLY: mscale, nscale, signgs, mnyq, nnyq,
10  1 successful_term_flag
11  USE realspace, ONLY: shalf, wint, guu, guv, gvv, r1, ru, rv,
12  1 zu , zv , phip
13 #ifdef _ANIMEC
14  2 ,pp3
15  USE fbal, ONLY: bimax_ppargrad
16 #endif
17 #ifdef NETCDF
18  USE ezcdf
19 #endif
20  USE xstuff, ONLY: xc
21  USE parallel_include_module
22  IMPLICIT NONE
23 !-----------------------------------------------
24 ! D u m m y A r g u m e n t s
25 !-----------------------------------------------
26  REAL(rprec), DIMENSION(ns,nznt), INTENT(in) ::
27  1 bsupu, bsupv, bsq, gsqrt, sigma_an
28 #ifdef _ANIMEC
29  2 ,ppar, onembc
30 #endif
31  REAL(rprec), DIMENSION(ns,nznt,0:1), TARGET, INTENT(inout) ::
32  1 bsubu, bsubv
33  REAL(rprec), DIMENSION(ns,nznt), INTENT(inout), TARGET :: bsubs
34  REAL(rprec), DIMENSION(ns,nznt), INTENT(out) ::
35  1 itheta, brho, izeta
36 #ifdef _ANIMEC
37  1 ,pp1, pp2
38 #endif
39  REAL(rprec), DIMENSION(ns,nznt,0:1) :: bsubsu, bsubsv
40  INTEGER, INTENT(in) :: ier_flag
41 !-----------------------------------------------
42 ! L o c a l P a r a m e t e r s
43 !-----------------------------------------------
44 !RESET lbsubs DEFAULT FLAG TO FALSE TO CAPTURE CURRENT SHEETS!
45 ! LOGICAL, PARAMETER :: lbsubs = .false. !!False to use (correct) bsubs calculation (from metrics)
46  !!True to use (modified) bsubs calculation (from mag. diff. eq.)
47 ! J Hanson 2014-01-12. Commented out above line. Variable is now declared
48 ! in module vmec_input, available here through module vmec_main.
49 ! lbsubs is now a namelist input variable, so user can change.
50 ! LOGICAL, PARAMETER :: lbsubs = .true. !!True to use NEW bsubs calculation (from mag. diff. eq.)
51 ! !!False to use OLD bsubs calculation (from metrics)
52  LOGICAL, PARAMETER :: lprint = .false. !!Prints out bsubs spectrum to fort.33
53  REAL(rprec), PARAMETER :: two=2, p5=0.5_dp, c1p5=1.5_dp
54 !-----------------------------------------------
55 ! L o c a l V a r i a b l e s
56 !-----------------------------------------------
57  INTEGER lk, lz, lt, k, m, js, j, n, injxbout, mparity, nznt1
58  INTEGER :: njxbout = jxbout0, nmin, info
59  INTEGER, PARAMETER :: ns_skip = 1, nu_skip = 1, nv_skip = 1
60  REAL(rprec), DIMENSION(:,:), ALLOCATABLE ::
61  1 bdotk, bsubuv, bsubvu
62  REAL(rprec), DIMENSION(:,:,:,:), ALLOCATABLE :: bsubsmn
63  REAL(rprec), DIMENSION(:,:,:), ALLOCATABLE :: brhomn,
64  1 bsubs3, bsubv3, bsubu3, jxb_gradp, jcrossb, sqrtg3,
65  2 bsupv3, bsupu3, jsups3, jsupv3, jsupu3, jdotb_sqrtg
66  REAL(rprec), POINTER :: bs1(:), bu1(:,:), bv1(:,:)
67  REAL(rprec), DIMENSION(:), ALLOCATABLE :: kperpu, kperpv,
68  2 sqgb2, sqrtg, kp2, jxb, jxb2, bsupu1, bsupv1, bsubu1, bsubv1,
69  3 avforce, aminfor, amaxfor, toroidal_angle, phin, pprim, pprime
70  REAL(rprec), DIMENSION(:,:), ALLOCATABLE :: bsubua, bsubva
71  REAL(rprec) ::
72  1 bsubsmn1, bsubsmn2, bsubvmn1, bsubvmn2, bsubumn1, bsubumn2,
73  1 bsubsmn3, bsubsmn4, bsubvmn3, bsubvmn4, bsubumn3, bsubumn4,
74  2 dnorm1, tcos1, tcos2, tsini1, tsini2, tcosi1, tcosi2,
75  3 tcosm1, tcosm2, tcosn1, tcosn2, tsinm1, tsinm2, tsin1, tsin2,
76  4 tsinn1, tsinn2, tjnorm, ovp, pnorm, brho00(ns)
77  REAL(rprec), DIMENSION(:,:), ALLOCATABLE ::
78  1 bsubu_s, bsubu_a, bsubv_s, bsubv_a
79  REAL(rprec), DIMENSION(:), ALLOCATABLE ::
80  1 bsubs_s, bsubs_a
81  CHARACTER(LEN=100) :: jxbout_file
82  CHARACTER(LEN=100) :: legend(13)
83  LOGICAL :: lprint_flag
84 !-----------------------------------------------
85 #ifdef NETCDF
86  CHARACTER(LEN=*), PARAMETER ::
87  1 vn_legend = 'legend',
88  1 vn_radial_surfaces = 'radial_surfaces',
89  1 vn_poloidal_grid_points = 'poloidal_grid_points',
90  1 vn_toroidal_grid_points = 'toroidal_grid_points',
91  1 vn_mpol = 'mpol',
92  1 vn_ntor = 'ntor',
93  1 vn_phin = 'phin',
94  1 vn_toroidal_angle = 'toroidal_angle',
95  1 vn_avforce = 'avforce',
96  1 vn_jdotb = 'surf_av_jdotb',
97  1 vn_sqg_bdotk = 'sqrt(g)*bdotk',
98  1 vn_sqrtg = 'sqrt(g)',
99  1 vn_bdotgradv = 'bdotgradv',
100  1 vn_amaxfor = 'amaxfor',
101  1 vn_aminfor = 'aminfor',
102  1 vn_pprime = 'pprime',
103  1 vn_jsupu = 'jsupu',
104  1 vn_jsupv = 'jsupv',
105  1 vn_jsups = 'jsups',
106  1 vn_bsupu = 'bsupu',
107  1 vn_bsupv = 'bsupv',
108  1 vn_jcrossb = 'jcrossb',
109  1 vn_jxb_gradp = 'jxb_gradp',
110  1 vn_bsubu = 'bsubu',
111  1 vn_bsubv = 'bsubv',
112  1 vn_bsubs = 'bsubs'
113 !-----------------------------------------------
114 #endif
115  lprint_flag = (ier_flag.eq.successful_term_flag)
116  IF (lprint_flag) THEN
117 #ifdef NETCDF
118  jxbout_file = 'jxbout_'//trim(input_extension)//'.nc'
119 
120  CALL cdf_open(njxbout,jxbout_file,'w',injxbout)
121 #else
122  jxbout_file = 'jxbout.'//trim(input_extension)//'.txt'
123  CALL safe_open(njxbout, injxbout, jxbout_file, 'replace',
124  1 'formatted')
125 #endif
126 #if !defined(_ANIMEC)
127  IF (any(sigma_an .NE. one)) stop 'SIGMA_AN != 1'
128 #endif
129  IF (injxbout .ne. 0) THEN
130  print *,' Error opening JXBOUT file in jxbforce'
131  RETURN
132  END IF
133 
134 ! PROGRAM FOR COMPUTING LOCAL KXB = grad-p FORCE BALANCE
135 !
136 ! sigma_an = one + (pperp-ppar)/(2*bsq(1:nrzt)) (=1 for isotropic pressure)
137 ! K = CURL(sigma_an B) is the "effective" current (=J for isotropic pressure)
138 ! Compute u (=theta), v (=zeta) derivatives of B sub s
139 !
140  legend(1) = " S = normalized toroidal flux (0 - 1)"
141  IF (lasym) THEN
142  legend(2) = " U = VMEC poloidal angle (0 - 2*pi, FULL period)"
143  ELSE
144  legend(2) = " U = VMEC poloidal angle (0 - pi, HALF a period)"
145  END IF
146  legend(3) = " V = VMEC (geometric) toroidal angle (0 - 2*pi)"
147  legend(4) = " SQRT(g') = |SQRT(g-VMEC)| / VOL':" //
148  1 " Cylindrical-to-s,u,v Jacobian normed to volume derivative"
149  legend(5) = " VOL = Int_s'=0,s Int_u Int_v |SQRT(g_VMEC)| :" //
150  1 " plasma volume enclosed by surface s'=s"
151  legend(6) = " VOL' = d(VOL)/ds: differential volume element"
152  legend(7) = " Es = SQRT(g') [grad(U) X grad(V)] : covariant" //
153  1 " radial unit vector (based on volume radial coordinate)"
154  legend(8) = " BSUP{U,V} = sigma_an B DOT GRAD{U,V}:" //
155  1 " contravariant components of B"
156  legend(9) = " JSUP{U,V} = SQRT(g') J DOT GRAD{U,V}"
157  legend(10)=
158  1 " K X B = Es DOT [K X B]: covariant component of K X B force"
159  legend(11)= " K * B = K DOT B * SQRT(g')"
160  legend(12)= " p' = dp/d(VOL): pressure gradient (based on" //
161  1 " volume radial coordinate)"
162  legend(13)= " <KSUP{U,V}> = Int_u Int_v [KSUP{U,V}]/dV/ds"
163 
164 #if !defined(NETCDF)
165  WRITE (njxbout,5) (ns1-1)/ns_skip, ntheta3/nu_skip, nzeta/nv_skip,
166  1 mpol, ntor, phiedge
167  5 FORMAT(/,' Radial surfaces = ',i3, ' Poloidal grid points = ',i3,
168  1 ' Toroidal grid points = ',i3,/,
169  2 ' Poloidal modes = ',i3,' Toroidal modes = ', i3,
170  3 ' Toroidal Flux = ',1pe12.3)
171  WRITE (njxbout, 6) (legend(j), j=1,13)
172  6 FORMAT(/,100('-'),/,' LEGEND:',/,100('-'),/,
173  1 2(3(a,/),/),5(a,/),/,2(a,/),100('-'),//)
174 #endif
175  ENDIF
176 
177  nznt1 = nzeta*ntheta2
178  ALLOCATE (avforce(ns),aminfor(ns),amaxfor(ns))
179  ALLOCATE (bdotk(ns,nznt), bsubuv(ns,nznt),
180  1 bsubvu(ns,nznt), kperpu(nznt), kperpv(nznt),
181  2 sqgb2(nznt), brhomn(0:mnyq,-nnyq:nnyq,0:1),kp2(nznt),
182  3 jxb(nznt), jxb2(nznt), bsupu1(nznt),
183  3 bsubua(nznt1,0:1), bsubva(nznt1,0:1),
184  4 bsupv1(nznt), bsubu1(nznt), bsubv1(nznt),
185  5 bsubsmn(ns,0:mnyq,-nnyq:nnyq,0:1),
186  6 bsubs_s(nznt), bsubs_a(nznt), sqrtg(nznt),
187  7 bsubu_s(nznt1,0:1), bsubu_a(nznt1,0:1),
188  8 bsubv_s(nznt1,0:1), bsubv_a(nznt1,0:1), stat=j)
189  IF (j .ne. 0) stop 'Allocation error in jxbforce'
190 
191 !
192 ! NOTE: bsubuv, bsubvu are used to compute the radial current (should be zero)
193 !
194  bsubsu = 0; bsubsv = 0; bsubuv = 0; bsubvu = 0; bdotk = 0
195  bsubs(1,:) = 0; bsubsmn = 0
196 
197  radial: DO js = 1, ns
198 !
199 ! Put bsubs on full mesh
200 !
201  IF (js.gt.1 .and. js.lt.ns) THEN
202  bsubs(js,:) = p5*(bsubs(js,:) + bsubs(js+1,:))
203  END IF
204 
205  bsubu(js,:,1) = bsubu(js,:,1)/shalf(js)
206  bsubv(js,:,1) = bsubv(js,:,1)/shalf(js)
207  bsubua = 0; bsubva = 0
208 
209 ! _s: symmetric in u,v _a: antisymmetric in u,v on half (ntheta2) interval
210  IF (lasym) THEN
211  bs1=>bsubs(js,:); bu1=>bsubu(js,:,:); bv1=>bsubv(js,:,:)
212  CALL fsym_fft (bs1, bu1, bv1, bsubs_s, bsubu_s, bsubv_s,
213  1 bsubs_a, bsubu_a, bsubv_a)
214  ELSE
215  bsubs_s(:) = bsubs(js,:)
216  bsubu_s = bsubu(js,:,:); bsubv_s = bsubv(js,:,:)
217  END IF
218 
219 !
220 ! FOURIER LOW-PASS FILTER bsubX
221 
222  DO m = 0, mpol1
223  mparity = mod(m, 2)
224  DO n = 0, ntor
225 !
226 ! FOURIER TRANSFORM
227 !
228  dnorm1 = one/r0scale**2
229  IF (m .eq. mnyq) dnorm1 = p5*dnorm1
230  IF (n.eq.nnyq .and. n.ne.0) dnorm1 = p5*dnorm1
231  bsubsmn1 = 0; bsubsmn2 = 0
232  IF (lasym) THEN
233  dnorm1 = 2*dnorm1 !SPH012314 in FixAray
234  bsubsmn3 = 0; bsubsmn4 = 0
235  END IF
236  bsubumn1 = 0; bsubumn2 = 0; bsubvmn1 = 0; bsubvmn2 = 0
237  IF (lasym) THEN
238  bsubumn3 = 0; bsubumn4 = 0; bsubvmn3 = 0; bsubvmn4 = 0
239  END IF
240 
241  DO k = 1, nzeta
242  lk = k
243  DO j = 1, ntheta2
244  tsini1 = sinmui(j,m)*cosnv(k,n)*dnorm1
245  tsini2 = cosmui(j,m)*sinnv(k,n)*dnorm1
246  tcosi1 = cosmui(j,m)*cosnv(k,n)*dnorm1
247  tcosi2 = sinmui(j,m)*sinnv(k,n)*dnorm1
248  bsubsmn1 = bsubsmn1 + tsini1*bsubs_s(lk)
249  bsubsmn2 = bsubsmn2 + tsini2*bsubs_s(lk)
250  bsubvmn1 = bsubvmn1 + tcosi1*bsubv_s(lk, mparity)
251  bsubvmn2 = bsubvmn2 + tcosi2*bsubv_s(lk, mparity)
252  bsubumn1 = bsubumn1 + tcosi1*bsubu_s(lk, mparity)
253  bsubumn2 = bsubumn2 + tcosi2*bsubu_s(lk, mparity)
254 
255  IF (lasym) THEN
256  bsubsmn3 = bsubsmn3 + tcosi1*bsubs_a(lk)
257  bsubsmn4 = bsubsmn4 + tcosi2*bsubs_a(lk)
258  bsubvmn3 = bsubvmn3 + tsini1*bsubv_a(lk, mparity)
259  bsubvmn4 = bsubvmn4 + tsini2*bsubv_a(lk, mparity)
260  bsubumn3 = bsubumn3 + tsini1*bsubu_a(lk, mparity)
261  bsubumn4 = bsubumn4 + tsini2*bsubu_a(lk, mparity)
262  END IF
263 
264  lk = lk + nzeta
265 
266  END DO
267  END DO
268 
269 !
270 ! FOURIER INVERSE TRANSFORM
271 ! Compute on u-v grid (must add symmetric, antisymmetric parts for lasym=T)
272 !
273  DO k = 1, nzeta
274  lk = k
275  DO j = 1, ntheta2
276  tcos1 = cosmu(j,m)*cosnv(k,n)
277  tcos2 = sinmu(j,m)*sinnv(k,n)
278  bsubua(lk,0) = bsubua(lk,0) + tcos1*bsubumn1 +
279  1 tcos2*bsubumn2
280  bsubva(lk,0) = bsubva(lk,0) + tcos1*bsubvmn1 +
281  1 tcos2*bsubvmn2
282 
283  tcosm1 = cosmum(j,m)*cosnv(k,n)
284  tcosm2 = sinmum(j,m)*sinnv(k,n)
285  bsubsu(js,lk,0) = bsubsu(js,lk,0) +
286  1 tcosm1*bsubsmn1 + tcosm2*bsubsmn2
287  tcosn1 = sinmu(j,m)*sinnvn(k,n)
288  tcosn2 = cosmu(j,m)*cosnvn(k,n)
289  bsubsv(js,lk,0) = bsubsv(js,lk,0) +
290  1 tcosn1*bsubsmn1 + tcosn2*bsubsmn2
291  bsubvu(js,lk) = bsubvu(js,lk) +
292  1 sinmum(j,m)*cosnv(k,n)*bsubvmn1 +
293  2 cosmum(j,m)*sinnv(k,n)*bsubvmn2
294  bsubuv(js,lk) = bsubuv(js,lk) +
295  1 cosmu(j,m)*sinnvn(k,n)*bsubumn1 +
296  2 sinmu(j,m)*cosnvn(k,n)*bsubumn2
297 
298  IF (lasym) THEN
299  tsin1 = sinmu(j,m)*cosnv(k,n)
300  tsin2 = cosmu(j,m)*sinnv(k,n)
301  bsubua(lk,1) = bsubua(lk,1) + tsin1*bsubumn3 +
302  1 tsin2*bsubumn4
303  bsubva(lk,1) = bsubva(lk,1) + tsin1*bsubvmn3 +
304  1 tsin2*bsubvmn4
305 
306  tsinm1 = sinmum(j,m)*cosnv(k,n)
307  tsinm2 = cosmum(j,m)*sinnv(k,n)
308  bsubsu(js,lk,1) = bsubsu(js,lk,1) +
309  1 tsinm1*bsubsmn3 + tsinm2*bsubsmn4
310  tsinn1 = cosmu(j,m)*sinnvn(k,n)
311  tsinn2 = sinmu(j,m)*cosnvn(k,n)
312  bsubsv(js,lk,1) = bsubsv(js,lk,1) +
313  1 tsinn1*bsubsmn3 + tsinn2*bsubsmn4
314  bsubvu(js,lk) = bsubvu(js,lk) +
315  1 cosmum(j,m)*cosnv(k,n)*bsubvmn3 +
316  2 sinmum(j,m)*sinnv(k,n)*bsubvmn4
317  bsubuv(js,lk) = bsubuv(js,lk) +
318  1 sinmu(j,m)*sinnvn(k,n)*bsubumn3 +
319  2 cosmu(j,m)*cosnvn(k,n)*bsubumn4
320  END IF
321 
322  lk = lk + nzeta
323 
324  END DO
325  END DO
326 
327 !
328 ! bsubsmn: coefficients of sin(mu)cos(nv), n>=0, cos(mu)sin(nv), n<0 (type=0)
329 ! cos(mu)cos(nv), n>=0, sin(mu)sin(nv), n<0 (type=1, nonzero only for lasym=T)
330 !
331  IF (.not.lprint) cycle !Don't need these except for comparison
332 
333  bsubsmn(js,m,n,0) = bsubsmn1
334  IF (n .gt. 0) bsubsmn(js,m,-n,0) = bsubsmn2
335 
336  IF (.not.lasym) cycle
337 
338  bsubsmn(js,m,n,1) = bsubsmn3
339  IF (n .gt. 0) bsubsmn(js,m,-n,0) = bsubsmn4
340 
341  END DO
342  END DO
343 
344  IF (lasym) THEN
345 ! EXTEND FILTERED bsubu, bsubv TO NTHETA3 MESH
346 ! NOTE: INDEX 0 - COS(mu-nv) SYMMETRY; 1 - SIN(mu-nv) SYMMETRY
347  CALL fext_fft (bsubu(js,:,0), bsubua(:,0), bsubua(:,1))
348  CALL fext_fft (bsubv(js,:,0), bsubva(:,0), bsubva(:,1))
349  ELSE
350  bsubu(js,:,0) = bsubua(:,0)
351  bsubv(js,:,0) = bsubva(:,0)
352  END IF
353 
354  END DO radial
355 
356  DEALLOCATE (bsubua, bsubva)
357 
358 ! EXTEND bsubsu, bsubsv TO NTHETA3 MESH
359  IF (lasym) CALL fsym_invfft (bsubsu, bsubsv)
360 
361 #ifdef _ANIMEC
362  CALL bimax_ppargrad(pp1, pp2, gsqrt, ppar, onembc, pres,
363  1 phot,tpotb)
364 #endif
365 
366 ! SKIPS Bsubs Correction - uses Bsubs from metric elements
367  IF (.not.lbsubs) GOTO 1500
368 
369 !
370 ! Compute corrected Bsubs coefficients (brhomn) (impacts currents)
371 ! by solving es dot (KXB - gradp_parallel) = 0 equation for brhomn in REAL SPACE
372 ! Can be written Bsupu d(bs)/du + Bsupv d(bs)/dv = RHS (jxb below), bs==bsubs
373 ! brho==sigma B_s, pp1 and pp2 are the Jacobian times the hot particle parallel
374 ! pressure radial gradient Amplitudes on the full integer mesh
375 !
376  correct_bsubs: DO js = 2, ns-1
377  jxb(:) = p5*(gsqrt(js,:) + gsqrt(js+1,:))
378  bsupu1(:) = p5*(bsupu(js,:)*gsqrt(js,:)
379  1 + bsupu(js+1,:)*gsqrt(js+1,:))
380  bsupv1(:) = p5*(bsupv(js,:)*gsqrt(js,:)
381  1 + bsupv(js+1,:)*gsqrt(js+1,:))
382  brho(js,:) = ohs*
383  1 ( bsupu1(:)*(bsubu(js+1,:,0) - bsubu(js,:,0))
384  2 + bsupv1(:)*(bsubv(js+1,:,0) - bsubv(js,:,0)))
385  3 + (pres(js+1) - pres(js))*ohs*jxb(:)
386 #ifdef _ANIMEC
387 !WAC Last two lines of brho contain hot particle parallel pressure gradients
388  4 + ohs*((pres(js+1)*phot(js+1) - pres(js)*phot(js)) * pp2(js,:)
389  5 + (tpotb(js+1) - tpotb(js) ) * pp1(js,:))
390 #endif
391 !
392 ! SUBTRACT FLUX-SURFACE AVERAGE FORCE BALANCE FROM brho, OTHERWISE
393 ! LOCAL FORCE BALANCE EQUATION B dot grad(Bs) = brho CAN'T BE SOLVED
394 !
395  brho00(js) = sum(brho(js,:)*wint(js:nrzt:ns))
396  brho(js,:) = brho(js,:) - signgs*jxb(:)*brho00(js)/
397  1 (p5*(vp(js) + vp(js+1)))
398 
399  jxb(:) = brho(js,:)
400  CALL getbsubs (brhomn, jxb, bsupu1, bsupv1, mnyq, nnyq, info)
401  IF (info .ne. 0) THEN
402  print *, 'Error in GETBRHO: info= ',info, ' js= ',js
403  ELSE IF (lprint) THEN
404  WRITE (33, *) ' JS = ', js
405  IF (lasym) THEN
406  WRITE (33, '(a)')
407  1 ' M N BSUBS(old) BSUBS(new)' //
408  2 ' BSUBS(old) BSUBS(new)'
409  ELSE
410  WRITE (33, *) ' M N BSUBS(old) BSUBS(new)'
411  END IF
412  DO m = 0, mpol1
413  DO n = -ntor, ntor
414  IF (lasym) THEN
415  WRITE(33,1223) m, n, bsubsmn(js,m,n,0), brhomn(m,n,0),
416  1 bsubsmn(js,m,n,1), brhomn(m,n,1)
417  ELSE
418  WRITE(33,1224) m, n, bsubsmn(js,m,n,0), brhomn(m,n,0)
419  END IF
420  END DO
421  END DO
422  END IF
423  1223 FORMAT (i4,1x,i4,4(6x,1p,e12.3))
424  1224 FORMAT (i4,1x,i4,2(6x,1p,e12.3))
425 
426 !
427 ! Recompute bsubsu,v now using corrected bsubs
428 ! Store old values (itheta,izeta) for checking force balance later
429 !
430  itheta(js,:) = bsubsu(js,:,0); izeta(js,:) = bsubsv(js,:,0)
431 
432  IF (info .ne. 0) cycle
433  bsubsu(js,:,:) = 0; bsubsv(js,:,:) = 0; bsubs_s = 0
434  IF (lasym) bsubs_a = 0;
435 
436  DO m = 0, mnyq
437  DO n = 0, nnyq
438  IF (n .eq. 0) THEN
439  bsubsmn1 = brhomn(m,0,0)
440  bsubsmn2 = 0
441  ELSE
442  bsubsmn1 = brhomn(m,n,0)
443  bsubsmn2 = brhomn(m,-n,0)
444  END IF
445 
446  IF (lasym) THEN
447  IF (n .eq. 0) THEN
448  bsubsmn3 = brhomn(m,0,1)
449  bsubsmn4 = 0
450  ELSE
451  bsubsmn3 = brhomn(m,n,1)
452  bsubsmn4 = brhomn(m,-n,1)
453  END IF
454  END IF
455 
456  DO k = 1, nzeta
457  lk = k
458  DO j = 1, ntheta2
459  tsin1 = sinmu(j,m)*cosnv(k,n)
460  tsin2 = cosmu(j,m)*sinnv(k,n)
461  bsubs_s(lk) = bsubs_s(lk) + tsin1*bsubsmn1
462  1 + tsin2*bsubsmn2
463  tcosm1 = cosmum(j,m)*cosnv(k,n)
464  tcosm2 = sinmum(j,m)*sinnv(k,n)
465  bsubsu(js,lk,0) = bsubsu(js,lk,0) +
466  1 tcosm1*bsubsmn1 + tcosm2*bsubsmn2
467  tcosn1 = sinmu(j,m)*sinnvn(k,n)
468  tcosn2 = cosmu(j,m)*cosnvn(k,n)
469  bsubsv(js,lk,0) = bsubsv(js,lk,0) +
470  1 tcosn1*bsubsmn1 + tcosn2*bsubsmn2
471 
472  IF (lasym) THEN
473  tcos1 = cosmu(j,m)*cosnv(k,n)
474  tcos2 = sinmu(j,m)*sinnv(k,n)
475  bsubs_a(lk) = bsubs_a(lk) + tcos1*bsubsmn3
476  1 + tcos2*bsubsmn4
477  tsinm1 = sinmum(j,m)*cosnv(k,n)
478  tsinm2 = cosmum(j,m)*sinnv(k,n)
479  bsubsu(js,lk,1) = bsubsu(js,lk,1) +
480  1 tsinm1*bsubsmn3 + tsinm2*bsubsmn4
481  tsinn1 = cosmu(j,m)*sinnvn(k,n)
482  tsinn2 = sinmu(j,m)*cosnvn(k,n)
483  bsubsv(js,lk,1) = bsubsv(js,lk,1) +
484  1 tsinn1*bsubsmn3 + tsinn2*bsubsmn4
485 
486  END IF
487 
488  lk = lk + nzeta
489 
490  END DO
491  END DO
492  END DO
493  END DO
494 
495  IF (lasym) THEN
496 ! EXTEND TO FULL (ntheta3) u-GRID
497  bs1 => bsubs(js,:)
498  CALL fext_fft (bs1, bsubs_a, bsubs_s)
499  ELSE
500  bsubs(js,:) = bsubs_s(:)
501  END IF
502 
503  END DO correct_bsubs
504 
505 ! EXTEND bsubsu, bsubsv TO NTHETA3 MESH
506  IF (lasym) CALL fsym_invfft (bsubsu, bsubsv)
507 
508 !
509 ! CHECK FORCE BALANCE: SQRT(g)*(bsupu*bsubsu + bsupv*bsubsv) = brho
510 !
511  IF (.not.lprint) GOTO 1500
512 
513  WRITE (33, '(/,2a,/)') 'ANGLE INDEX B*grad(Bs) Frhs',
514  1 ' Fold'
515  check_fb: DO js = 2, ns-1
516  bsupu1(:) = p5*(bsupu(js,:)*gsqrt(js,:)
517  1 + bsupu(js+1,:)*gsqrt(js+1,:))
518  bsupv1(:) = p5*(bsupv(js,:)*gsqrt(js,:)
519  1 + bsupv(js+1,:)*gsqrt(js+1,:))
520  kp2(:) = bsupu1(:)*bsubsu(js,:,0) + bsupv1(:)*bsubsv(js,:,0)
521  jxb(:) = bsupu1(:)*itheta(js,:) + bsupv1(:)*izeta(js,:)
522 
523  WRITE (33, '(/,a,i4)') 'JS = ',js
524  DO lk = 1, nznt
525  WRITE(33,1230) lk, brho(js,lk), kp2(lk), jxb(lk)
526  END DO
527 
528  END DO check_fb
529 
530  1230 FORMAT (i9,5x, 1p,3e14.4)
531 
532  1500 CONTINUE
533 
534  DEALLOCATE (bsubs_s, bsubs_a, bsubu_s,
535  1 bsubu_a, bsubv_s, bsubv_a, stat=lk)
536 
537 !
538 ! Compute end point values for bsubs
539 !
540  bsubs(1,:) = 2*bsubs(2,:) - bsubs(3,:)
541  bsubs(ns,:) = 2*bsubs(ns,:) - bsubs(ns-1,:)
542 !
543 ! Now compute currents on the FULL radial mesh
544 ! Here:
545 !
546 ! Itheta = sqrt(g) * Ksupu
547 ! Izeta = sqrt(g) * Ksupv
548 ! Ksupx = K dot grad(x) x=(u,v)
549 ! jxb = (K X B) dot (grad-u X grad-v) sqrt(g)
550 ! bdotk = sigma*sqrt(g)*K dot B
551 ! kperpx = (B X gradp) dot grad(x) / |B|**2 x=(u,v)
552 ! sqgb2 = sigma*sqrt(g)*|B|**2
553 ! sqrtg = sqrt(g)
554 ! pprime = d(p||)/dV
555 !
556 ! kp2 == |k-perp|**2 = kperpu**2 * guu + 2*kperpu*kperpv*guv + kperpv**2 * gvv
557 ! This was compared to the alternative expression (agreed very well):
558 ! |j-perp|**2 = |grad-s|**2 * (dp/ds)**2 / |B|**2
559 !
560 ! Note: Multiply currents, pressure by 1/mu0 to get in mks units!
561 ! TWOPI*TWOPI factor incorporated in vp (thru ovp factor below), so V' = (2pi)**2*vp
562 !
563 #ifdef NETCDF
564  ALLOCATE(
565  1 bsubs3(ns,nzeta,ntheta3), bsubv3(ns,nzeta,ntheta3),
566  2 bsubu3(ns,nzeta,ntheta3), jxb_gradp(ns,nzeta,ntheta3),
567  3 jcrossb(ns,nzeta,ntheta3), bsupv3(ns,nzeta,ntheta3),
568  4 bsupu3(ns,nzeta,ntheta3), jsups3(ns,nzeta,ntheta3),
569  5 jsupv3(ns,nzeta,ntheta3), jsupu3(ns,nzeta,ntheta3),
570  6 jdotb_sqrtg(ns,nzeta,ntheta3), sqrtg3(ns,nzeta,ntheta3),
571  7 phin(ns), toroidal_angle(nzeta), stat=j)
572 
573  bsubs3=0; bsubv3=0; bsubu3=0; jxb_gradp=0
574  jcrossb=0 ; bsupv3=0; bsupu3=0; jsups3=0
575  jsupv3=0; jsupu3=0; phin=0; phin(ns)=1
576  jdotb_sqrtg=0; sqrtg3=0
577 #endif
578 
579  ALLOCATE (pprime(nznt), pprim(ns),stat=j)
580  pprim=0
581 
582  avforce=0; aminfor=0; amaxfor=0
583  dnorm1 = twopi*twopi
584 
585  DO js = 2, ns1
586  ovp = two/(vp(js+1) + vp(js))/dnorm1
587  tjnorm = ovp*signgs
588  sqgb2(:nznt) = sigma_an(js+1,:)*gsqrt(js+1,:)*
589  1 (bsq(js+1,:)- pres(js+1))
590  2 + sigma_an(js,:)*gsqrt(js,:) *
591  3 (bsq(js,:) - pres(js))
592 #ifdef _ANIMEC
593  sqgb2(:nznt) = sigma_an(js+1,:nznt)*gsqrt(js+1,:nznt)
594  1 * bsq(js+1,:nznt)
595  2 + sigma_an(js ,:nznt)*gsqrt(js ,:nznt)
596  3 * bsq(js ,:nznt)
597 #elif defined _FLOW
598 ! sqgb2(:nznt) = sigma_an(js+1,:nznt)*gsqrt(js+1,:nznt)
599 ! 1 * (bsq(js+1,:nznt)-prot(js+1,:nznt))
600 ! 2 + sigma_an(js ,:nznt)*gsqrt(js ,:nznt)
601 ! 3 * (bsq(js ,:nznt)-prot(js ,:nznt))
602 #else
603  sqgb2(:nznt) = sigma_an(js+1,:nznt)*gsqrt(js+1,:nznt)
604  1 * (bsq(js+1,:nznt)-pres(js+1))
605  2 + sigma_an(js ,:nznt)*gsqrt(js ,:nznt)
606  3 * (bsq(js ,:nznt)-pres(js ))
607 #endif
608 ! TAKE THIS OUT: MAY BE POORLY CONVERGED AT THIS POINT....
609 ! IF (ANY(sqgb2(:nznt)*signgs .le. zero))
610 ! 1 STOP ' SQGB2 <= 0 in JXBFORCE'
611  pprime(:) = ohs*(pres(js+1)-pres(js))/mu0 !dp/ds here
612 #ifdef _ANIMEC
613 !WAC Last two lines of 'pprime' contain the hot particle parallel pressure
614  1 + ohs*((pres(js+1)*phot(js+1) - pres(js)*phot(js))*pp2(js,:nznt)
615  2 + (tpotb(js+1) - tpotb(js) ) *pp1(js,:nznt))
616  3 / mu0
617 #endif
618  kperpu(:nznt) = p5*(bsubv(js+1,:nznt,0) + bsubv(js,:nznt,0))*
619  1 pprime(:)/sqgb2
620  kperpv(:nznt) =-p5*(bsubu(js+1,:nznt,0) + bsubu(js,:nznt,0))*
621  1 pprime(:)/sqgb2
622  kp2(:nznt)=p5*(kperpu**2*(guu(js+1:nrzt:ns) + guu(js:nrzt:ns))
623  1 + 2*kperpu*kperpv*(guv(js+1:nrzt:ns) + guv(js:nrzt:ns))
624  2 + kperpv**2*(gvv(js+1:nrzt:ns) + gvv(js:nrzt:ns)))
625  itheta(js,:nznt) = bsubsv(js,:nznt,0) - ohs*
626  1 (bsubv(js+1,:nznt,0) - bsubv(js,:nznt,0))
627  izeta(js,:nznt) = -bsubsu(js,:nznt,0) + ohs*
628  1 (bsubu(js+1,:nznt,0) - bsubu(js,:nznt,0))
629  itheta(js,:nznt) = itheta(js,:nznt)/mu0
630  izeta(js,:nznt) = izeta(js,:nznt)/mu0
631  sqrtg(:) = p5*(gsqrt(js,:) + gsqrt(js+1,:))
632  bsupu1(:nznt) = p5*(bsupu(js+1,:nznt)*gsqrt(js+1,:)
633  1 + bsupu(js,:nznt) *gsqrt(js,:)) / sqrtg(:)
634  bsupv1(:nznt) = p5*(bsupv(js+1,:nznt)*gsqrt(js+1,:)
635  1 + bsupv(js,:nznt) *gsqrt(js,:)) / sqrtg(:)
636  bsubu1(:nznt) = p5*(bsubu(js+1,:nznt,0) + bsubu(js,:nznt,0))
637  bsubv1(:nznt) = p5*(bsubv(js+1,:nznt,0) + bsubv(js,:nznt,0))
638  jxb(:nznt) = ovp*(itheta(js,:nznt) * bsupv1(:nznt)
639  1 - izeta(js,:nznt) * bsupu1(:nznt))
640  bdotk(js,:nznt) = itheta(js,:nznt) * bsubu1(:nznt) +
641  1 izeta(js,:nznt) * bsubv1(:nznt)
642  pprime(:nznt) = ovp*pprime(:nznt)
643  pnorm = one/(abs(pprime(1)) + epsilon(pprime(1)))
644  amaxfor(js) = maxval(jxb(:nznt)-pprime(:))*pnorm
645  aminfor(js) = minval(jxb(:nznt)-pprime(:))*pnorm
646  avforce(js) = sum(wint(2:nrzt:ns)*(jxb(:nznt) - pprime(:)))
647  amaxfor(js) = 100*min(amaxfor(js),9.999_dp)
648  aminfor(js) = 100*max(aminfor(js),-9.999_dp)
649  pprim(js) = sum(wint(js:nrzt:ns)*pprime(:))
650 ! Compute <K dot B>, <B sup v> = signgs*phip
651 ! jpar2 = <j||**2>, jperp2 = <j-perp**2>, with <...> = flux surface average
652 
653  jdotb(js) = dnorm1*tjnorm*sum(bdotk(js,:nznt)*wint(2:nrzt:ns)
654  1 / sigma_an(js,:nznt))
655  bdotb(js) = dnorm1*tjnorm*sum(sqgb2(:nznt)*wint(2:nrzt:ns)
656  1 / sigma_an(js,:nznt))
657 
658  bdotgradv(js) = p5*dnorm1*tjnorm*(phip(js) + phip(js+1))
659  jpar2(js) = dnorm1*tjnorm*
660  1 sum(bdotk(js,:nznt)**2*wint(2:nrzt:ns)
661  2 /(sigma_an(js,:nznt)*sqgb2(:nznt)))
662  jperp2(js)= dnorm1*tjnorm*
663  1 sum(kp2(:nznt)*wint(2:nrzt:ns)*sqrtg(:nznt))
664 
665  IF (mod(js,ns_skip) .eq. 0 .and. lprint_flag) THEN
666 #ifdef NETCDF
667  phin(js) = phi(js)/phi(ns)
668  DO lz = 1, nzeta
669  toroidal_angle(lz)=real(360*(lz-1),rprec)/nzeta
670  DO lt = 1, ntheta3
671  lk = lz + nzeta*(lt-1)
672 C lu (js,lz,lt ) = lt
673  jsupu3(js,lz,lt) = ovp*itheta(js,lk)
674  jsupv3(js,lz,lt) = ovp*izeta(js,lk)
675  jsups3(js,lz,lt) = ovp*(bsubuv(js,lk)
676  1 - bsubvu(js,lk))/mu0
677  bsupu3(js,lz,lt) = bsupu1(lk)
678  bsupv3(js,lz,lt) = bsupv1(lk)
679  jcrossb(js,lz,lt) = jxb(lk)
680  jxb_gradp(js,lz,lt) = (jxb(lk) - pprime(lk))
681  jdotb_sqrtg(js,lz,lt) = ovp*bdotk(js,lk)
682  sqrtg3(js,lz,lt) = sqrtg(lk)*ovp
683  bsubu3(js,lz,lt) = bsubu(js,lk,0)
684  bsubv3(js,lz,lt) = bsubv(js,lk,0)
685  bsubs3(js,lz,lt) = bsubs(js,lk)
686  END DO
687  END DO
688 #else
689  WRITE (njxbout, 200) phi(js), avforce(js), jdotb(js),
690  1 bdotgradv(js), pprime(1), one/ovp,
691  2 (twopi**2)*tjnorm*sum(itheta(js,:)*wint(js:nrzt:ns)),
692  3 (twopi**2)*tjnorm*sum(izeta(js,:)*wint(js:nrzt:ns)),
693  4 amaxfor(js), aminfor(js)
694  WRITE (njxbout, 90)
695  DO lz = 1, nzeta, nv_skip
696  WRITE (njxbout, 100) real(360*(lz-1),rprec)/nzeta, lz
697  DO lt = 1, ntheta3, nu_skip
698  lk = lz + nzeta*(lt - 1)
699  WRITE (njxbout, 110) lt, tjnorm*itheta(js,lk),
700  1 tjnorm*izeta(js,lk), ovp*(bsubuv(js,lk) -
701  2 bsubvu(js,lk))/mu0, bsupu1(lk), bsupv1(lk),
702  3 sqrtg(lk)*ovp, jxb(lk), jxb(lk) - pprime(lk),
703  4 ovp*bdotk(js,lk), bsubu(js,lk,0),
704  5 bsubv(js,lk,0), bsubs(js,lk)
705  END DO
706  END DO
707 #endif
708  ENDIF
709  END DO
710 
711  izeta(1,:nznt) = two*izeta(2,:nznt) - izeta(3,:nznt) !!Need in wrout
712  izeta(ns,:nznt)= two*izeta(ns-1,:nznt) - izeta(ns-2,:nznt) !!Need in wrout
713  jdotb(1) = two*jdotb(2) - jdotb(3)
714  jdotb(ns) = two*jdotb(ns-1) - jdotb(ns-2)
715  bdotb(1) = two*bdotb(3) - bdotb(2)
716  bdotb(ns) = two*bdotb(ns-1) - bdotb(ns-2)
717  bdotgradv(1) = two*bdotgradv(2) - bdotgradv(3)
718  bdotgradv(ns) = two*bdotgradv(ns-1) - bdotgradv(ns-2)
719  jpar2(1) = 0; jpar2(ns) = 0; jperp2(1) = 0; jperp2(ns) = 0
720  pprim(1) = 2*pprim(ns-1) - pprim(ns-2)
721  pprim(ns) = 2*pprim(ns-1) - pprim(ns-2)
722 
723  IF (lprint_flag) THEN
724 #ifdef NETCDF
725  CALL cdf_define(njxbout,vn_legend,legend)
726  CALL cdf_define(njxbout,vn_mpol,mpol)
727  CALL cdf_define(njxbout,vn_ntor,ntor)
728  CALL cdf_define(njxbout,vn_phin,phin)
729  CALL cdf_define(njxbout,vn_radial_surfaces,ns)
730  CALL cdf_define(njxbout,vn_poloidal_grid_points,ntheta3)
731  CALL cdf_define(njxbout,vn_toroidal_grid_points,nzeta)
732  CALL cdf_define(njxbout,vn_avforce,avforce)
733  CALL cdf_define(njxbout,vn_jdotb,jdotb)
734 
735  CALL cdf_define(njxbout,vn_sqg_bdotk,jdotb_sqrtg)
736  CALL cdf_define(njxbout,vn_sqrtg,sqrtg3)
737 
738  CALL cdf_define(njxbout,vn_bdotgradv,bdotgradv)
739  CALL cdf_define(njxbout,vn_pprime,pprim)
740  CALL cdf_define(njxbout,vn_aminfor,aminfor)
741  CALL cdf_define(njxbout,vn_amaxfor,amaxfor)
742  CALL cdf_define(njxbout,vn_jsupu,jsupu3)
743  CALL cdf_define(njxbout,vn_jsupv,jsupv3)
744  CALL cdf_define(njxbout,vn_jsups,jsups3)
745  CALL cdf_define(njxbout,vn_bsupu,bsupu3)
746  CALL cdf_define(njxbout,vn_bsupv,bsupv3)
747  CALL cdf_define(njxbout,vn_jcrossb,jcrossb)
748  CALL cdf_define(njxbout,vn_jxb_gradp,jxb_gradp)
749  CALL cdf_define(njxbout,vn_bsubu,bsubu3)
750  CALL cdf_define(njxbout,vn_bsubv,bsubv3)
751  CALL cdf_define(njxbout,vn_bsubs,bsubs3)
752 
753  CALL cdf_write(njxbout,vn_legend,legend)
754  CALL cdf_write(njxbout,vn_mpol,mpol)
755  CALL cdf_write(njxbout,vn_ntor,ntor)
756  CALL cdf_write(njxbout,vn_phin,phin)
757  CALL cdf_write(njxbout,vn_radial_surfaces,ns)
758  CALL cdf_write(njxbout,vn_poloidal_grid_points,ntheta3)
759  CALL cdf_write(njxbout,vn_toroidal_grid_points,nzeta)
760  CALL cdf_write(njxbout,vn_avforce,avforce)
761  CALL cdf_write(njxbout,vn_jdotb,jdotb)
762 
763  CALL cdf_write(njxbout,vn_sqg_bdotk,jdotb_sqrtg)
764  CALL cdf_write(njxbout,vn_sqrtg,sqrtg3)
765 
766  CALL cdf_write(njxbout,vn_bdotgradv,bdotgradv)
767  CALL cdf_write(njxbout,vn_pprime,pprim)
768  CALL cdf_write(njxbout,vn_aminfor,aminfor)
769  CALL cdf_write(njxbout,vn_amaxfor,amaxfor)
770  CALL cdf_write(njxbout,vn_jsupu,jsupu3)
771  CALL cdf_write(njxbout,vn_jsupv,jsupv3)
772  CALL cdf_write(njxbout,vn_jsups,jsups3)
773  CALL cdf_write(njxbout,vn_bsupu,bsupu3)
774  CALL cdf_write(njxbout,vn_bsupv,bsupv3)
775  CALL cdf_write(njxbout,vn_jcrossb,jcrossb)
776  CALL cdf_write(njxbout,vn_jxb_gradp,jxb_gradp)
777  CALL cdf_write(njxbout,vn_bsubu,bsubu3)
778  CALL cdf_write(njxbout,vn_bsubv,bsubv3)
779  CALL cdf_write(njxbout,vn_bsubs,bsubs3)
780 
781  CALL cdf_close(njxbout)
782 
783  DEALLOCATE(
784  1 bsubs3, bsubv3, bsubu3, jxb_gradp, jcrossb, bsupv3,
785  2 bsupu3, jsups3, jsupv3, jsupu3, jdotb_sqrtg, phin,
786  3 toroidal_angle, sqrtg3, stat=j)
787 
788 #else
789  CLOSE (njxbout)
790 
791  90 FORMAT(/" LU JSUPU JSUPV JSUPS BSUPU",
792  1 " BSUPV SQRT(g') J X B J X B - p' J * B",
793  2 " BSUBU BSUBV BSUBS "/)
794  100 FORMAT( " TOROIDAL ANGLE (PER PERIOD) = ", f8.3," DEGREES",
795  1 " (PLANE #", i3,")")
796  110 FORMAT(i5,1p,12e11.3)
797  200 FORMAT(/" TOROIDAL FLUX = ",1p,e12.3,3x,"<J X B - p'> = ",
798  1 e12.3,3x,"<J DOT B> = ",e12.3,3x,
799  2 "<B DOT GRAD(V)> = ",e12.3,/,
800  2 " dp/d(VOL) [p'] = ",e12.3,3x,'d(VOL)/ds = ',e12.3,3x,
801  2 "<JSUPU> = ",e12.3,3x,"<JSUPV> = ",e12.3,/,
802  3 " MAXIMUM FORCE DEVIATIONS (RELATIVE TO p'): ",sp,0p,f7.2,"%",
803  4 3x,f7.2,"%")
804 #endif
805 
806  END IF
807 
808  DEALLOCATE (kperpu, kperpv, sqgb2, sqrtg, kp2, brhomn, bsubsmn,
809  1 jxb, jxb2, bsupu1, bsupv1, bsubu1, bsubv1, avforce, aminfor,
810  2 amaxfor, pprim, stat=j)
811 !
812 ! COMPUTE MERCIER CRITERION
813 !
814  bdotk = mu0*bdotk
815  CALL mercier(gsqrt,bsq,bdotk,iotas,wint,r1,ru,rv,zu,zv,bsubu,
816  1 vp,phips,pres,ns,nznt)
817 
818  DEALLOCATE (bdotk, bsubuv, bsubvu, pprime, stat=j)
819 
820  END SUBROUTINE jxbforce
Generated on Thu Mar 5 2020 15:49:27 for V3FIT by   doxygen 1.8.17