V3FIT
eqfor.f
1  SUBROUTINE eqfor(br, bz, bsubu, bsubv, tau, rzl_array, ier_flag)
2  USE vmec_main
3  USE vmec_params
4  USE realspace
5  USE vforces, r12 => armn_o, bsupu => crmn_e, bsupv => czmn_e,
6  1 gsqrt => azmn_o, bsq => bzmn_o, izeta => azmn_e,
7  2 brho => bzmn_e, bphi => czmn_o, curtheta => brmn_e
8 #ifdef _ANIMEC
9  3 ,tau_an => brmn_o
10 #endif
11  USE vacmod
12  USE vspline
13  USE csplinx
14  USE vmec_io
15  USE mgrid_mod
16  USE fbal
17  USE v3f_vmec_comm
18  USE stel_constants, ONLY: pi
19 
20  USE xstuff, ONLY: xc, pxc
21  USE safe_open_mod
22  USE timer_sub
23  IMPLICIT NONE
24 C-----------------------------------------------
25 C D u m m y A r g u m e n t s
26 C-----------------------------------------------
27  INTEGER :: ier_flag
28  REAL(dp), DIMENSION(ns,nznt,0:1), INTENT(in) :: bsubu, bsubv
29  REAL(dp), DIMENSION(nrzt), INTENT(out) :: br, bz
30  REAL(dp), DIMENSION(nrzt), INTENT(out) :: tau
31  REAL(dp), DIMENSION(ns,0:ntor,0:mpol1,3*ntmax), TARGET,
32  & INTENT(in) :: rzl_array
33 C-----------------------------------------------
34 C L o c a l P a r a m e t e r s
35 C-----------------------------------------------
36  REAL(dp), PARAMETER :: c1p5=1.5_dp
37 C-----------------------------------------------
38 C L o c a l V a r i a b l e s
39 C-----------------------------------------------
40  INTEGER :: i, icount, itheta, js, js1, l, loff,
41  & lpi, lt, n, n1, nchicur, nchiiota0, noff,
42  & nout, nsort, iv, iu, lk, nplanes
43  REAL(dp), DIMENSION(:), POINTER ::
44  & rmags, zmags, rmaga, zmaga
45  REAL(dp), DIMENSION(:,:,:), POINTER :: rmncc,zmnsc
46  REAL(dp), DIMENSION(ns) :: phi1, chi1, jPS2
47  REAL(dp) :: modb(nznt)
48  REAL(dp), DIMENSION(:), ALLOCATABLE ::
49  & btor_vac, btor1, dbtor, phat, t12u, guu_1u, surf_area,
50  & r3v, redge, rbps1u, bpol2vac, phipf_loc
51  REAL(dp) :: aminr1, !aminr2,
52  & aminr2in, anorm,
53  & aspectratio, betai, betstr, scaling_ratio,
54  & bminz2, bminz2in, btor, iotamax, musubi,
55  & bzcalc, bzin, chisq, chiwgt, cur0,
56  & delphid_exact, delta1, delta2, delta3, denwgt, lambda,
57  & dlogidsi, !dmusubi_meas,
58  & er, es, fac, facnorm, factor, fgeo,
59  & fmax, fmin, flao, fpsi0, pavg, pitchc, pitchm,
60  & pprime, qedge, qmin1, qmin2, qmin3, qzero,
61  & raxis0, rcalc, rcen,
62  & rcenin, rgeo, rjs,
63  & rjs1, rlao, rqmin1, rqmin2, rshaf, rshaf1, rshaf2, s11, s12,
64  & s13, s2, s3, sigr0, sigr1, sigz1, smaleli,
65  & splintx, splints, sqmin, sumbpol, sumbtot, sumbtor, sump,
66  & sump2, sump20, t1, tz, jpar_perp=0, jparps_perp=0,
67  & tol, toroidal_flux, vnorm, vprime, wght0, xmax,
68  & xmida, xmidb, xmin, rzmax, rzmin, zxmax, zxmin, zaxis0,
69  & zmax, zmin, yr1u, yz1u, waist(2), height(2)
70  REAL(dp) :: d_of_kappa, tmpxc, rmssum
71  INTEGER :: istat1, OFU, j, k
72 C-----------------------------------------------
73 C E x t e r n a l F u n c t i o n s
74 C-----------------------------------------------
75  EXTERNAL splintx, splints
76 C-----------------------------------------------
77 !
78 ! POINTER ASSOCIATIONS
79 !
80  rmags => rzl_array(1,:,0,rcc)
81  zmags => rzl_array(1,:,0,zcs+ntmax)
82  rmncc => rzl_array(:,:,:,rcc)
83  zmnsc => rzl_array(:,:,:,zsc+ntmax)
84  IF (lasym) THEN
85  rmaga => rzl_array(1,:,0,rcs)
86  zmaga => rzl_array(1,:,0,zcc+ntmax)
87  END IF
88 
89 ! crmn_o => bss on half grid
90  CALL bss(r12, bzmn, brmn, azmn, armn, crmn_o, bsupu, bsupv,
91  & br, bphi, bz)
92 
93 !
94 ! STORE EDGE VALUES OF B-FIELD
95 !
96  IF ((lfreeb .and. ivac .gt. 1) .or. ledge_dump) THEN
97  IF (ALLOCATED(bredge)) THEN
98  DEALLOCATE(bredge, bpedge, bzedge)
99  END IF
100  ALLOCATE(bredge(2*nznt), bpedge(2*nznt), bzedge(2*nznt),
101  & stat=i)
102  IF (i .ne. 0) stop 'Error in EQFOR allocating bredge'
103  DO iv = 1,nzeta
104  DO iu = 1, ntheta3
105  lk = iv + nzeta*(iu - 1)
106  n1 = ns*lk
107  bredge(lk) = 1.5_dp*br(n1) - p5*br(n1 - 1)
108  bpedge(lk) = 1.5_dp*bphi(n1) - p5*bphi(n1 - 1)
109  bzedge(lk) = 1.5_dp*bz(n1) - p5*bz(n1 - 1)
110  END DO
111  END DO
112  END IF
113 
114 #ifdef _ANIMEC
115 !EVALUATE MIRROR STABILITY CRITERION; BSQ ==> MAGNETIC PRESSURE
116  CALL mirror_crit(tau_an, bsq)
117 #endif
118 !
119 ! NOTE: JXBFORCE ROUTINE MUST BE CALLED TO COMPUTE IZETA, JDOTB
120 ! ON OUTPUT, J, IZETA, JDOTB ARE IN MKS UNITS (1/MU0 FACTOR)
121 !
122 ! CAUTION: THIS CALL WILL WRITE OVER br, bz
123 !
124 
125  CALL jxbforce(bsupu, bsupv, bsubu, bsubv, crmn_o, rcon, zcon,
126  & gsqrt, bsq, curtheta, izeta, brho, sigma_an,
127  & ier_flag
128 #ifdef _ANIMEC
129  & , pp1, pp2, ppar, onembc
130 #endif
131  & )
132 
133 !
134 ! HALF-MESH VOLUME-AVERAGED BETA
135 !
136 
137  tau(1) = 0
138  tau(2:nrzt) = signgs*wint(2:nrzt)*gsqrt(2:nrzt)
139  DO i = 2, ns
140  s2 = sum(bsq(i:nrzt:ns)*tau(i:nrzt:ns))/vp(i) - pres(i)
141  overr(i) = sum(tau(i:nrzt:ns)/r12(i:nrzt:ns)) / vp(i)
142  beta_vol(i) = pres(i)/s2
143  END DO
144 
145  betaxis = c1p5*beta_vol(2) - p5*beta_vol(3)
146 
147  IF (grank.EQ. 0) WRITE (nthreed, 5)
148  5 FORMAT(/,' NOTE: S=normalized toroidal flux (0 - 1)',/,
149  1 ' U=poloidal angle (0 - 2*pi)',/,
150  1 ' V=geometric toroidal angle (0 - 2*pi)',/,
151  1 ' <RADIAL FORCE> = d(Ipol)/dPHI',
152  1 ' - IOTA*d(Itor)/dPHI - dp/dPHI * d(VOL)/dPHI',/,
153  1 ' = d(VOL)/dPHI*[<JSUPU>',
154  1 ' - IOTA*<JSUPV> - SIGN(JAC)*dp/dPHI]',/,
155  1 ' (NORMED TO SUM OF INDIVIDUAL TERMS)',//,
156  1 ' S <RADIAL TOROIDAL IOTA ',
157  1 ' <JSUPU> <JSUPV> d(VOL)/',
158  2 ' d(PRES)/ <M> PRESF <BSUBU> <BSUBV>',
159  3 ' <J.B> <B.B>',/,
160  4 ' FORCE> FLUX ',
161  5 ' d(PHI) ',
162  6 ' d(PHI) ',/,148('-'),/)
163  ALLOCATE (phipf_loc(ns))
164 
165  phipf_loc(1) = twopi*signgs*(c1p5*phip(2) - p5*phip(3))
166  presf(1) = c1p5*pres(2) - p5*pres(3)
167  DO i = 2,ns1
168  presf(i) = p5*(pres(i) + pres(i + 1))
169  phipf_loc(i) = p5*twopi*signgs*(phip(i) + phip(i + 1))
170  END DO
171  presf(ns) = c1p5*pres(ns) - p5*pres(ns - 1)
172  phipf_loc(ns) = twopi*signgs*(c1p5*phip(ns) - p5*phip(ns1))
173 
174  phi1(1) = zero
175  chi1(1) = zero
176  DO i = 2, ns
177  phi1(i) = phi1(i - 1) + hs*phip(i)
178  chi1(i) = chi1(i - 1) + hs*(phip(i)*iotas(i))
179  END DO
180 
181  chi = twopi*chi1
182 
183 ! WRITE (36, 201) (i, phi1(i), chi1(i), iotas(i), i=1, ns)
184 !201 FORMAT (i4, 1p,3e14.6)
185 
186  CALL calc_fbal(bsubu, bsubv)
187 
188  bucof(1) = 0
189  bvcof(1) = c1p5*bvco(2) - p5*bvco(3)
190 !
191 ! NOTE: jcuru, jcurv on FULL radial mesh coming out of calc_fbal
192 ! They are local (surface-averaged) current densities (NOT integrated in s)
193 ! jcurX = (dV/ds)/twopi**2 <JsupX> for X=u,v
194 !
195  DO i = 2, ns1
196  equif(i) = equif(i)*vpphi(i)/(abs(jcurv(i)*chipf(i)) +
197  & abs(jcuru(i)*phipf(i)) +
198  & abs(presgrad(i)*vpphi(i)))
199  bucof(i) = p5*(buco(i) + buco(i + 1))
200  bvcof(i) = p5*(bvco(i) + bvco(i + 1))
201  END DO
202 
203  bucof(ns) = c1p5*buco(ns) - p5*buco(ns1)
204  bvcof(ns) = c1p5*bvco(ns) - p5*bvco(ns1)
205 
206  equif(1) = two*equif(2) - equif(3)
207  jcuru(1) = two*jcuru(2) - jcuru(3)
208  jcurv(1) = two*jcurv(2) - jcurv(3)
209  presgrad(1) = two*presgrad(2) - presgrad(3)
210  presgrad(ns) = two*presgrad(ns1) - presgrad(ns1-1)
211  vpphi(1) = two*vpphi(2) - vpphi(3)
212  vpphi(ns) = two*vpphi(ns1) - vpphi(ns1-1)
213  equif(ns) = two*equif(ns1) - equif(ns1-1)
214  jcuru(ns) = two*jcuru(ns1) - jcuru(ns1-1)
215  jcurv(ns) = two*jcurv(ns1) - jcurv(ns1-1)
216 ! NOTE: phipf = phipf_loc/(twopi), phipf_loc ACTUAL (twopi factor) Toroidal flux derivative
217 ! SPH/JDH (060211): remove twopi factors from <JSUPU,V> (agree with output in JXBOUT file)
218  fac = twopi*signgs
219  DO js = 1, ns
220  es = (js - 1)*hs
221  cur0 = fac*vpphi(js)*twopi !==dV/ds = dV/dPHI * d(PHI/ds) (V=actual volume)
222  IF (rank .EQ. 0) THEN
223  WRITE (nthreed, 30) es, equif(js), fac*phi1(js), iotaf(js),
224  & jcuru(js)/vpphi(js)/mu0,
225  & jcurv(js)/vpphi(js)/mu0,
226  & cur0/phipf_loc(js),
227  & presgrad(js)/phipf_loc(js)/mu0,
228  & specw(js), presf(js)/mu0, bucof(js),
229  & bvcof(js), jdotb(js), bdotb(js)
230  END IF
231  END DO
232  30 FORMAT(1p,2e10.2,2e12.4,4e11.3,0p,f7.3,1p,5e11.3)
233 
234  DEALLOCATE(phipf_loc)
235 
236 !
237 ! MAKE SURE WOUT FILE DOES NOT REQUIRE ANY STUFF COMPUTED BELOW....
238 !
239  IF (ier_flag .NE. successful_term_flag) RETURN
240 
241 !
242 ! Calculate mean (toroidally averaged) poloidal cross section area & toroidal flux
243 !
244  anorm = twopi*hs
245  vnorm = twopi*anorm
246  toroidal_flux = anorm*sum(bsupv(2:nrzt)*tau(2:nrzt))
247 
248 !
249 ! Calculate poloidal circumference and normal surface area and aspect ratio
250 ! Normal is | dr/du X dr/dv | = SQRT [R**2 guu + (RuZv - RvZu)**2]
251 !
252  ALLOCATE(guu_1u(nznt), surf_area(nznt))
253  guu_1u(:nznt) = ru0(ns:nrzt:ns)*ru0(ns:nrzt:ns) +
254  1 zu0(ns:nrzt:ns)*zu0(ns:nrzt:ns)
255  surf_area(:nznt) = wint(ns:nrzt:ns)*sqrt(guu_1u(:nznt))
256  circum_p = twopi*sum(surf_area(:nznt))
257  surf_area(:nznt) = wint(ns:nrzt:ns)*sqrt(
258  1 + (r1(ns:nrzt:ns,0) + r1(ns:nrzt:ns,1))**2*guu_1u(:nznt)
259  2 +((rv(ns:nrzt:ns,0) + rv(ns:nrzt:ns,1))*zu0(ns:nrzt:ns)
260  3 - (zv(ns:nrzt:ns,0) + zv(ns:nrzt:ns,1))*ru0(ns:nrzt:ns))**2)
261  surf_area_p = twopi**2*sum(surf_area(:nznt))
262  DEALLOCATE (guu_1u)
263 
264  aspect = aspectratio()
265 
266 ! Also, estimate mean elongation of plasma from the following relations
267 ! for an axisymmetric torus with elliptical cross section and semi-axes
268 ! a and a * kappa (kappa >= 1)
269 !
270 ! surf_area _p = 2*pi*R * 2*pi*a ctwiddle(kappa_p)
271 ! volume_p = 2*pi*R * pi*a ** 2 * kappa_p
272 ! cross_area _p = pi*a ** 2 * kappa_p
273 !
274 ! The cirumference of an ellipse of semi-axes a and a * kappa_p is
275 ! 2 * pi * a ctwiddle(kappa_p)
276 ! The exact form for ctwiddle is 4 E(1 - kappa_p^2) / (2 pi), where
277 ! E is the complete elliptic integral of the second kind
278 ! (with parameter argument m, not modulus argument k)
279 !
280 ! The coding below implements an approximate inverse of the function
281 ! d(kappa) = ctwiddle(kappa) / sqrt(kappa)
282 ! The approximate inverse is
283 ! kappa = 1 + (pi^2/8) * (d^2+sqrt(d^4-1)-1)
284 ! Note that the variable aminor_p, for an elliptic cross section,
285 ! would be a * sqrt(kappa)
286 !
287  d_of_kappa = surf_area_p * aminor_p / ( 2 * volume_p)
288  kappa_p = 1 + (pi * pi / 8) * (d_of_kappa ** 2 + sqrt &
289  & (abs(d_of_kappa ** 4 - 1)) -1)
290  vvc_kappa_p = kappa_p ! Save result for v3fit.
291 
292  aminr1 = 2*volume_p/surf_area_p
293 
294 !
295 !
296 ! OUTPUT BETAS, INDUCTANCES, SAFETY FACTORS, ETC.
297 ! (EXTRACTED FROM FQ-CODE, 9-10-92)
298 !
299 ! b poloidals (cylindrical estimates)
300 !
301 ! rcen = p5*(router + rinner) !geometric center
302  n = 0
303  n1 = n + 1
304  rcenin = dot_product(rmncc(ns,n1,:mpol1+1:2),
305  1 mscale(:mpol1:2)*nscale(n))
306 
307  l = (mpol1+1)/2
308  ALLOCATE (t12u(l))
309  t12u(:l) = mscale(1:mpol1:2)*nscale(n)
310  aminr2in = dot_product(rmncc(ns,n1,2:mpol1+1:2),t12u(:l))
311  bminz2in = dot_product(zmnsc(ns,n1,2:mpol1+1:2),t12u(:l))
312  bminz2 = dot_product(zmnsc(ns,n1,2:mpol1+1:2),t12u(:l))
313  DEALLOCATE (t12u)
314  aminr1 = sqrt(two*volume_p/(twopi*twopi*r00)) !vol av minor radius
315 ! aminr2 = p5*(router - rinner) !geometric plasma radius
316 !
317 ! cylindrical estimates for beta poloidal
318 !
319  sump = vnorm*sum(vp(2:ns)*pres(2:ns))
320  pavg = sump/volume_p
321 ! ppeak = presf(1)
322  factor = 2*pavg
323 !
324 ! delphid_exact = Integral[ (Bvac - B) * dSphi ]
325 ! rshaf [= RT in Eq.(12), Phys Fluids B 5 (1993) 3119]
326 !
327 ! Note: tau = |gsqrt|*wint
328 !
329  ALLOCATE (btor_vac(nznt), btor1(nznt), dbtor(nznt),
330  1 phat(nznt), redge(nznt))
331  delphid_exact = zero !Eq. 20 in Shafranov
332  musubi = zero
333  rshaf1 = zero
334  rshaf2 = zero
335  DO js = 2, ns
336  btor_vac(:nznt) = rbtor/r12(js:nrzt:ns)
337  btor1(:nznt) = r12(js:nrzt:ns)*bsupv(js:nrzt:ns)
338  delphid_exact = delphid_exact + sum(
339  1 (btor_vac(:nznt)/r12(js:nrzt:ns) -
340  2 bsupv(js:nrzt:ns))*tau(js:nrzt:ns))
341  dbtor(:nznt) = btor1(:nznt)**2 - btor_vac(:nznt)**2
342  musubi = musubi - sum(dbtor(:nznt)*tau(js:nrzt:ns))
343  phat(:nznt) = bsq(js:nrzt:ns) - p5*btor_vac(:nznt)**2
344  phat(:nznt) = (phat(:nznt) - dbtor(:nznt))*tau(js:nrzt:ns)
345  rshaf1 = rshaf1 + sum(phat(:nznt))
346  rshaf2 = rshaf2 + sum(phat(:nznt)/r12(js:nrzt:ns))
347  END DO
348 
349  redge(:nznt) = r1(ns:nrzt:ns,0) + r1(ns:nrzt:ns,1)
350  IF (lfreeb .and. ivac.gt.1) THEN
351  phat = bsqvac(:) - p5*(bsubvvac/redge(:))**2
352  ELSE
353  phat = c1p5*bsq(ns:nrzt:ns) - p5*bsq(ns-1:nrzt:ns)
354  1 - p5*(rbtor/redge(:))**2
355  END IF
356 
357  DEALLOCATE (btor_vac, btor1, dbtor)
358 
359  delphid_exact = anorm*delphid_exact
360  rshaf = rshaf1/rshaf2
361  fpsi0 = c1p5*bvco(2) - p5*bvco(3)
362  b0 = fpsi0/r00
363 
364  rmax_surf = maxval(r1(ns:nrzt:ns,0)+r1(ns:nrzt:ns,1))
365  rmin_surf = minval(r1(ns:nrzt:ns,0)+r1(ns:nrzt:ns,1))
366  zmax_surf = maxval(abs(z1(ns:nrzt:ns,0)+z1(ns:nrzt:ns,1)))
367 
368  DO js = 2, ns
369  modb(:nznt) = sqrt(two*(bsq(js:nrzt:ns)-pres(js)))
370  CALL bextrema (modb, bmin(1,js), bmax(1,js), nzeta, ntheta2)
371  END DO
372 
373 !
374 ! output geometrical, |B| quantities
375 !
376  CALL elongation (r1, z1, waist, height)
377 
378  IF(rank.EQ.0) WRITE (nthreed, 75) bmin(1,ns), bmax(1,ns),
379  1 bmin(ntheta2,ns), bmax(ntheta2,ns)
380  75 FORMAT(/
381  1 ' Magnetic field modulation (averaged over toroidal angle)',/,
382  2 1x,71('-')/,' BMIN(u=0) = ',f14.6/
383  3 ' BMAX(u=0) = ',f14.6/' BMIN(u=pi) = ',
384  4 f14.6/' BMAX(u=pi) = ',f14.6/)
385 
386  sumbtot = 2*(vnorm*sum(bsq(2:nrzt)*tau(2:nrzt)) - sump)
387  sumbtor = vnorm*sum(tau(2:nrzt)*(r12(2:nrzt)*bsupv(2:nrzt))**2)
388  sumbpol = sumbtot - sumbtor
389  betapol = 2*sump/sumbpol
390  sump20 = 2*sump
391  sump2 = sum(pres(2:ns)*pres(2:ns)*vp(2:ns)*vnorm)
392  betatot = sump20/sumbtot
393  betator = sump20/sumbtor
394  volavgb = sqrt(abs(sumbtot/volume_p))
395  ionlarmor = 0.0032_dp/volavgb
396  jps2(2:ns1) = (jpar2(2:ns1) - jdotb(2:ns1)**2/bdotb(2:ns1))
397  jpar_perp = sum(jpar2(2:ns1)*(vp(2:ns1) + vp(3:ns)))
398  jparps_perp = sum(jps2(2:ns1)*(vp(2:ns1) + vp(3:ns)))
399  s2 = sum(jperp2(2:ns1)*(vp(2:ns1) + vp(3:ns)))
400  IF (s2 .ne. zero) THEN
401  jpar_perp = jpar_perp/s2
402  jparps_perp = jparps_perp/s2
403  END IF
404  IF (ntor .gt. 1) THEN
405  IF(rank.EQ.0) WRITE (nthreed, 80) aspect, kappa_p, volume_p,
406  1 cross_area_p, surf_area_p, circum_p, rmajor_p, aminor_p,
407  2 rmin_surf, rmax_surf, zmax_surf, waist(1), height(1),
408  3 waist(2), height(2)
409  ELSE
410  IF(rank.EQ.0) WRITE (nthreed, 80) aspect, kappa_p, volume_p,
411  1 cross_area_p, surf_area_p, circum_p, rmajor_p, aminor_p,
412  2 rmin_surf, rmax_surf, zmax_surf, waist(1), height(1)
413  END IF
414  80 FORMAT(/,' Geometric and Magnetic Quantities',/,1x,71('-')/,
415  1 ' Aspect Ratio = ',f14.6, /
416  1 ' Mean Elongation = ',f14.6, /
417  1 ' Plasma Volume = ',f14.6,' [M**3]',/
418  2 ' Cross Sectional Area = ',f14.6,' [M**2]',/
419  3 ' Normal Surface Area = ',f14.6,' [M**2]',/
420  4 ' Poloidal Circumference= ',f14.6,' [M]',/
421  5 ' Major Radius = ',f14.6,' [M]',
422  6 ' (from Volume and Cross Section)',/
423  7 ' Minor Radius = ',f14.6,' [M]',
424  8 ' (from Cross Section)',/
425  9 ' Minimum (inboard) R = ',f14.6,' [M]',/
426  a ' Maximum (outboard) R = ',f14.6,' [M]',/
427  a ' Maximum height Z = ',f14.6,' [M]',/
428  b ' Waist (v = 0) in R = ',f14.6,' [M]',/
429  b ' Full Height(v = 0) = ',f14.6,' [M]',:,/
430  b ' Waist (v = pi) in R = ',f14.6,' [M]',:,/
431  b ' Full Height(v = pi) = ',f14.6,' [M]')
432  IF (rank.EQ.0) WRITE (nthreed, 85) toroidal_flux,
433  1 1.e-6_dp*ctor/mu0, rbtor, rbtor0, volavgb, ionlarmor,
434  2 jpar_perp, jparps_perp
435  85 FORMAT(
436  1 ' Toroidal Flux = ',f14.6,' [Wb]',/
437  1 ' Toroidal Current = ',f14.6,' [MA]',/
438  1 ' RBtor(s=1) = ',f14.6,' [T-m]',/
439  1 ' RBtor(s=0) = ',f14.6,' [T-m]',/
440  1 ' Volume Average B = ',f14.6,' [T]',/
441  2 ' Ion Larmor Radius = ',f14.6,' [M] X Ti(keV)**0.5',/
442  3 ' <J||**2>/<J-perp**2> = ',f14.6,' (Vol. Averaged)',/
443  4 ' <JPS**2>/<J-perp**2> = ',f14.6,' (Vol. Averaged)',/)
444 
445  IF(rank.EQ.0) WRITE (nthreed, 90)
446  90 FORMAT(/,71('-'),/,' MORE GEOMETRIC AND PHYSICS QUANTITIES',/,
447  1 71('-'),/,' Toroidal Plane: Phi = 0',/,
448  1 5x,'j',3x,'psi-psiaxis',9x,'a [M]',3x,'ellipticity',3x,
449  2 'indentation',7x,'d-shape',4x,'rel. shift',6x,'<J||**2>/',4x,
450  3 '<JPS**2>/',/,95x,
451  4 '<J-perp**2>',3x,'<J-perp**2>'/,' -----',8(2x,12('-')))
452 
453  fac = twopi*hs*signgs
454  psi(1) = zero
455  ALLOCATE (r3v(ns-1))
456  r3v(:ns-1) = fac*phip(2:ns)*iotas(2:ns)
457  DO i = 1, ns - 1
458  psi(1+i) = psi(i) + r3v(i)
459  END DO
460  DEALLOCATE (r3v)
461 
462 ! nphi-plane, noff = 1,....,nzeta
463  DO nplanes = 1, 2
464  IF (nplanes .eq. 1) THEN
465  noff = 1
466  ELSE
467  IF (nzeta .eq. 1) EXIT
468  IF(rank.EQ.0) WRITE (nthreed, 95)
469  noff = 1 + nzeta/2
470  END IF
471 
472  ygeo(1) = zero
473  DO js = 2, ns
474  zmin = huge(zmin)
475  zmax = -huge(zmax)
476  xmin = huge(xmin)
477  xmax = -huge(xmax)
478  rzmax = zero
479 
480 ! Theta = 0 to pi in upper half of X-Z plane
481  DO icount = 1,2
482  n1 = noff !nphi-plane, n1 = noff,...,nzeta
483  IF (icount .eq. 2)
484  1 n1 = mod(nzeta + 1 - noff,nzeta) + 1 !(twopi-v), reflected plane
485  loff = js + ns*(n1-1)
486  t1 = one
487  IF (icount .eq. 2) t1 = -one
488  DO itheta = 1,ntheta2
489  yr1u = r1(loff,0) + sqrts(js)*r1(loff,1)
490  yz1u = z1(loff,0) + sqrts(js)*z1(loff,1)
491  yz1u = t1*yz1u
492  IF (yz1u .ge. zmax) THEN
493  zmax = abs(yz1u)
494  rzmax = yr1u
495  ELSE IF (yz1u .le. zmin) THEN
496  zmin = yz1u
497  rzmin = yr1u
498  END IF
499  IF (yr1u .ge. xmax) THEN
500  xmax = yr1u
501  zxmax = yz1u
502  ELSE IF (yr1u .le. xmin) THEN
503  xmin = yr1u
504  zxmin = yz1u
505  END IF
506  loff = loff + ns*nzeta
507  END DO
508  END DO
509 
510 
511  lpi = ns*((noff-1) + nzeta*(ntheta2 - 1))
512  lt = ns*((noff-1))
513  xmida = r1(js+lpi,0) + sqrts(js)*r1(js+lpi,1)
514  xmidb = r1(js+lt,0) + sqrts(js)*r1(js+lt,1)
515 
516  rgeo = p5*(xmidb + xmida) !Geometric major radius
517  ygeo(js) = p5*(xmidb - xmida) !Geometric minor radius
518 
519  yinden(js) = (xmida - xmin)/(xmax - xmin) !Geometric indentation
520  yellip(js) = (zmax - zmin)/(xmax - xmin) !Geometric ellipticity
521 
522  ytrian(js) = (rgeo - rzmax)/(xmax - xmin) !Geometric triangularity
523  yshift(js) = (r1(1+lt,0)-rgeo)/(xmax - xmin) !Geometric shift
524 
525  IF (jperp2(js) .eq. zero) jperp2(js) = epsilon(jperp2(js))
526  jpar_perp = jpar2(js)/jperp2(js)
527  IF (js .lt. ns) THEN
528  jparps_perp = jps2(js)/jperp2(js)
529  ELSE
530  jparps_perp = zero
531  END IF
532 
533  IF (nplanes .eq. 1) THEN
534  IF(rank.EQ.0) WRITE (nthreed, 120) js, psi(js), ygeo(js),
535  1 yellip(js), yinden(js), ytrian(js), yshift(js), jpar_perp,
536  2 jparps_perp
537  ELSE
538  IF(rank.EQ.0) WRITE (nthreed, 120) js, psi(js), ygeo(js),
539  1 yellip(js), yinden(js), ytrian(js), yshift(js)
540  END IF
541 
542  END DO
543  END DO
544 
545  95 FORMAT(/,71('-'),/,' Toroidal Plane: Phi = 180/Nfp',/,71('-'),/)
546  120 FORMAT(1x,i5,6f14.5,1p,3e14.2)
547 
548  IF(rank.EQ.0) WRITE (nthreed, 130)
549  130 FORMAT(//,' Magnetic Fields and Pressure',/,1x,71('-'))
550  fac = p5/mu0
551  IF(rank.EQ.0) WRITE (nthreed, 140) sump/mu0, pavg/mu0,
552  1 fac*sumbpol, fac*sumbpol/volume_p, fac*sumbtor, fac*sumbtor/
553  2 volume_p, fac*sumbtot, fac*sumbtot/volume_p, c1p5*sump/mu0,
554  3 c1p5*pavg/mu0
555  140 FORMAT(' Volume Integrals (Joules) and Volume ',
556  1 'Averages (Pascals)',/,24x,'Integral',6x,'Average',/,
557  2 ' pressure = ',1p,2e14.6,/,' bpol**2 /(2 mu0) = ',
558  3 2e14.6,/,' btor**2/(2 mu0) = ',2e14.6,/,
559  4 ' b**2/(2 mu0) = ',2e14.6,/,' EKIN (3/2p) = ',
560  5 2e14.6,/)
561 
562  IF(rank.EQ.0) WRITE (nthreed, 800)
563  800 FORMAT(/,' MAGNETIC AXIS COEFFICIENTS'/,
564  1 ' n rmag zmag rmag zmag',/,
565  2 ' (cos nv) (sin nv) (sin nv) (cos nv)',/)
566  loff = lbound(rmags,1)
567  IF (rank. eq. 0) THEN
568  DO n = 0, ntor
569  n1 = n + loff
570  t1 = mscale(0)*nscale(n)
571  tz = t1
572  IF (.NOT.lthreed) tz = 0
573  IF (lasym) THEN
574  WRITE (nthreed, 820) n, t1*rmags(n1),
575  1 tz*zmags(n1), tz*rmaga(n1), t1*zmaga(n1)
576  ELSE
577  WRITE (nthreed, 820) n, t1*rmags(n1),
578  1 tz*zmags(n1)
579  END IF
580  END DO
581  END IF
582  820 FORMAT(i5,1p,4e12.4)
583 
584  betstr = two*sqrt(sump2/volume_p)/(sumbtot/volume_p)
585 
586  IF(rank.EQ.0) WRITE (nthreed, 150) betatot, betapol, betator
587  150 FORMAT(/,' From volume averages over plasma, betas are',/,
588  1 ' beta total = ',f14.6,/,' beta poloidal = ',f14.6,/,
589  2 ' beta toroidal = ',f14.6,/)
590 
591  IF(rank.EQ.0) WRITE (nthreed, 160) rbtor, betaxis, betstr
592  160 FORMAT(' R * Btor-vac = ',f14.6,' [Wb/M]',/,
593  1 ' Peak Beta = ',f14.6,/,
594  2 ' Beta-star = ',f14.6,/)
595 
596 !
597 !
598 ! Shafranov surface integrals s1,s2
599 ! Plasma Physics vol 13, pp 757-762 (1971)
600 ! Also, s3 = .5*S3, defined in Lao, Nucl. Fusion 25, p.1421 (1985)
601 ! Note: if ctor = 0, use Int(Bsupu*Bsubu dV) for ctor*ctor/R
602 ! Phys. Fluids B, Vol 5 (1993) p 3121, Eq. 9a-9d
603 !
604 
605  ALLOCATE (rbps1u(nznt), bpol2vac(nznt))
606  IF (lfreeb .and. ivac.gt.1) THEN
607  bpol2vac = 2*bsqvac - bphiv*bphiv
608  ELSE
609  bpol2vac = 2*(c1p5*bsq(ns:nrzt:ns) - p5*bsq(ns-1:nrzt:ns))
610  1 - ((c1p5*bsupv(ns:nrzt:ns) - p5*bsupv(ns-1:nrzt:ns))
611  2 * redge)**2
612  END IF
613 
614 ! Compute current-like norm (factor) in Eq.(8), <a> * int(Bpol**2 * dA)
615 ! where <a> == 2*pi*Rs in Eq. 8 is the effective minor radius = Vol/Asurf
616 ! (corrects wrong description of Rs in paper, which is NOT the major radius)
617 ! This aminr1 = 1/2 the "correct" aminr1
618  aminr1 = volume_p/surf_area_p
619  factor = twopi**2*aminr1*sum(bpol2vac*surf_area)
620  factor = one/factor
621  facnorm = factor*twopi**2
622 
623 ! Lao's definition of normalization factor
624  scaling_ratio = (mu0*curtor/circum_p)**2*volume_p
625  scaling_ratio = scaling_ratio*factor
626 
627  rbps1u(:nznt) = facnorm*redge(:nznt)*phat(:nznt)
628  1 *wint(ns:nznt*ns:ns)
629  sigr0 = sum(rbps1u(:nznt)*zu0(ns:nrzt:ns))
630  sigr1 = sum(rbps1u(:nznt)*zu0(ns:nrzt:ns)*redge(:nznt))
631  sigz1 =-sum(rbps1u(:nznt)*ru0(ns:nrzt:ns)*
632  1 (z1(ns:nrzt:ns,0) + z1(ns:nrzt:ns,1)))
633  DEALLOCATE (redge, phat, rbps1u, bpol2vac, surf_area)
634 
635  er = sigr1 + sigz1
636  rlao = volume_p/(twopi*cross_area_p) !LAO, NUCL.FUS.25(1985)1421
637  flao = rshaf/rlao
638 ! fgeo = rshaf/rcen
639 
640  smaleli = factor*sumbpol
641  vvc_smaleli = smaleli ! Save result for v3fit.
642  betai = 2*factor*sump
643  musubi = vnorm*factor*musubi
644 ! dmusubi_meas = 2*twopi*factor*delphid*rbtor
645  lambda = p5*smaleli + betai
646  s11 = er - rshaf*sigr0 !Shafranov def. based on RT, Eq.(12)
647  s12 = er - rcen*sigr0 !R = Rgeometric
648  s13 = er - rlao*sigr0 !R = RLao
649  s2 = sigr0*rshaf
650  s3 = sigz1 !1/2 S3 in Eq.(14c)
651  delta1 = zero
652  delta2 = one - fgeo
653  delta3 = one - flao
654  IF(rank.EQ.0) WRITE (nthreed, 168)
655  IF(rank.EQ.0) WRITE (nthreed, 170) rshaf, rcen, rlao,
656  1 scaling_ratio, s3, smaleli, musubi, betai, lambda
657 ! IF (lrecon.AND.rank.EQ.0) WRITE (nthreed, 172) dmusubi_meas
658  IF(rank.EQ.0) WRITE (nthreed, 174) delta1, delta2, delta3,
659  1 s11, s12, s13, s2, s2/fgeo, s2/flao,
660  5 musubi + s11,musubi + s12,
661  6 musubi + s13,
662  8 p5*s11 + s2, p5*s12 + s2/fgeo, p5*s13 + s2/flao,
663  a p5*(3*betai+smaleli-musubi)/(s11+s2) - one,
664  b p5*(3*betai+smaleli-musubi)/(s12+s2/fgeo) - one,
665  c p5*(3*betai+smaleli-musubi)/(s13+s2/flao) - one,
666  d p5*(betai+smaleli+musubi)/s2 - one,
667  e p5*fgeo*(betai+smaleli+musubi)/s2 - one,
668  f p5*flao*(betai+smaleli+musubi)/s2 - one
669 
670  168 FORMAT(' Shafranov Surface Integrals',/
671  1 ' Ref: S. P. Hirshman, Phys. Fluids B, 5, (1993) 3119',/,
672  2 ' Note: s1 = S1/2, s2 = S2/2, where ',
673  3 ' s1,s2 are the Shafranov definitions,',/,
674  4 ' and s3 = S3/2, where S3 is Lao''s definition.',/,
675  5 ' The quantity lsubi gives the ratio of volume poloidal',
676  6 /,' field energy to the field energy estimated from the',
677  7 /,' surface integral in Eq.8.',/,1x,22('-'),/)
678 
679  170 FORMAT(
680  5 ' RT (Pressure-weighted) = ',f14.6,' [M]',/,
681  6 ' RG (Geometric) = ',f14.6,' [M]',/,
682  7 ' RL (Vol/2*pi*Area-Lao) = ',f14.6,' [M]',/,
683  8 ' Poloidal Field Energy',/,
684  8 ' Normalization Ratio = ',f14.6,' (Lao/Hirshman)',//,
685  8 ' s3 = ',f14.6,/,
686  9 ' lsubi = ',f14.6,/,
687  6 ' musubi = ',f14.6,/,
688  7 ' betai = ',f14.6,/,
689  9 ' lambda = ',f14.6,/)
690  172 FORMAT(' musubi (diamagnetism) = ',f14.6)
691  174 FORMAT(/,32x,'R = RT',12x,'R = RG',12x,'R = RL',/,
692  1 20x,3(10x,8('-')),/,
693  1 ' delta = 1 - RT/R = ',3(f14.6,4x),/,
694  2 ' s1 = ',3(f14.6,4x),/,
695  3 ' s2 = ',3(f14.6,4x),/,
696  8 ' betai (Mui + s1) = ',3(f14.6,4x),/,
697  a ' lambda (s1/2 + s2) = ',3(f14.6,4x),/,
698  b ' 1st Shafr''v relation = ',3(f14.6,4x),/,
699  c ' (3*Betai + Li - Mui)/[2*(s1+s2)] - 1',/,
700  d ' Radial force balance = ',3(f14.6,4x),/,
701  e ' (Betai + Li + Mui)/(2*s2) - 1',/)
702 
703  END SUBROUTINE eqfor
Generated on Thu Mar 5 2020 15:49:27 for V3FIT by   doxygen 1.8.17