2 USE stel_kinds,
ONLY: dp
3 REAL(dp),
DIMENSION(:),
ALLOCATABLE :: rzu_fac, rru_fac,
8 SUBROUTINE calc_fbal_par(bsubu, bsubv)
9 USE vmec_main,
ONLY: buco, bvco, equif, iequi,
10 1 jcurv, jcuru, chipf, vp, pres,
11 2 phipf, vpphi, presgrad, ohs
12 USE vmec_params,
ONLY: signgs
13 USE vmec_dim,
ONLY: ns, nrzt, nznt, ns1
14 USE realspace,
ONLY: pwint, phip
15 USE vmec_input,
ONLY: nzeta
16 USE vmec_dim,
ONLY: ntheta3
17 USE parallel_include_module
20 REAL(dp),
INTENT(in) :: bsubu(nznt,ns),
26 INTEGER :: nsmin, nsmax
29 nsmin=t1lglob; nsmax=t1rglob
31 buco(js) = sum(bsubu(:,js)*pwint(:,js))
32 bvco(js) = sum(bsubv(:,js)*pwint(:,js))
35 CALL gather1xarray(bvco)
36 CALL gather1xarray(buco)
41 nsmin=max(2,tlglob); nsmax=min(trglob,ns-1)
43 jcurv(js) = (signgs*ohs)*(buco(js+1) - buco(js))
44 jcuru(js) =-(signgs*ohs)*(bvco(js+1) - bvco(js))
45 vpphi(js) = (vp(js+1) + vp(js))/2
46 presgrad(js) = (pres(js+1) - pres(js))*ohs
47 equif(js) = (-phipf(js)*jcuru(js) + chipf(js)*jcurv(js))
48 1 /vpphi(js) + presgrad(js)
55 END SUBROUTINE calc_fbal_par
57 SUBROUTINE calc_fbal(bsubu, bsubv)
58 USE vmec_main,
ONLY: buco, bvco, equif,
59 1 jcurv, jcuru, chipf, vp, pres,
60 2 phipf, vpphi, presgrad, ohs
62 3 ,pmap, pd, phot, tpotb, zero
64 USE vmec_params,
ONLY: signgs
65 USE vmec_dim,
ONLY: ns, nrzt, nznt, ns1
66 USE realspace,
ONLY: wint, phip
68 1 ,pperp, ppar, onembc, sigma_an,
70 USE vforces, gsqrt => azmn_o
75 REAL(dp),
INTENT(in) :: bsubu(1:nrzt), bsubv(1:nrzt)
80 INTEGER :: nsmin, nsmax
86 buco(js) = sum(bsubu(js:nrzt:ns)*wint(js:nrzt:ns))
87 bvco(js) = sum(bsubv(js:nrzt:ns)*wint(js:nrzt:ns))
93 jcurv(js) = (signgs*ohs)*(buco(js+1) - buco(js))
94 jcuru(js) =-(signgs*ohs)*(bvco(js+1) - bvco(js))
96 vpphi(js) = (vp(js+1) + vp(js))/2
97 presgrad(js) = (pres(js+1) - pres(js))*ohs
98 equif(js) = (-phipf(js)*jcuru(js) + chipf(js)*jcurv(js))
99 1 /vpphi(js) + presgrad(js)
104 END SUBROUTINE calc_fbal
107 SUBROUTINE bimax_ppargrad(pp1, pp2, pp3, ppar, onembc, pres, phot,
109 USE vmec_main,
ONLY: zero
110 USE vmec_dim,
ONLY: ns, nrzt, nznt, ns1
114 REAL(dp),
INTENT(in) :: pp3(nrzt), ppar(nrzt), phot(ns),
115 1 onembc(nrzt) ,pres(ns), tpotb(ns)
116 REAL(dp),
INTENT(out):: pp1(nrzt), pp2(nrzt)
121 REAL(dp),
ALLOCATABLE :: tmp0(:), tmp2(:)
131 ALLOCATE (tmp0(nrzt), tmp2(nrzt))
133 pp1(1:nrzt:ns) = 0; pp2(1:nrzt:ns) = 0
135 tmp0(js:nrzt:ns) = pp3(js:nrzt:ns)*(ppar(js:nrzt:ns)-pres(js))
136 tmp2(js:nrzt:ns) = tmp0(js:nrzt:ns)/(pres(js)*phot(js)+eps)
140 pp2(l) = 0.5_dp*(tmp2(l)+tmp2(l+1))
143 tmp2(l) = pp3(l)*(1-onembc(l))
149 IF (onembc(l) <= zero)
THEN
151 & (1 - tpotb(js) * onembc(l)+eps)
153 pp1(l)=(2*tpotb(js )*tmp0(l)*onembc(l)
154 & + pres(js ) * phot(js ) * tmp2(l)*
155 & ( 1 - 5*(tpotb(js )*onembc(l))**1.5_dp))
156 & / (1 - (tpotb(js )*onembc(l))**2 + eps)
163 pp1(l) = 0.5_dp*(pp1(l)+pp1(l+1))
166 DEALLOCATE (tmp0, tmp2)
168 END SUBROUTINE bimax_ppargrad
170 SUBROUTINE mirror_crit(taumir, bsq)
171 USE stel_kinds,
ONLY: rprec, dp
172 USE realspace,
ONLY: sigma_an, pperp, ppar, onembc
174 USE vmec_main,
ONLY: tpotb, pppr, pmap, pres, papr, vp, wb,
175 & wp, wpar, wper, ns, nznt, nrzt,
184 REAL(dp),
DIMENSION(nrzt),
INTENT(out) :: taumir
185 REAL(dp),
DIMENSION(nrzt),
INTENT(inout) :: bsq
189 INTEGER :: js , lk , l
190 REAL(dp) :: whpar, whper, bhotmx, taumin, ppprht, parpht,
191 & sigmin, bhotc, betath, betaht, betato, betapa,
192 & betape, betadi, betaew, es, betprc, betprx
220 bhotmx = -huge(bhotmx)
221 taumin = huge(taumin)
222 sigmin = huge(sigmin)
227 taumir(l) = one + (pperp(l) - pres(js)) / bsq(l)*
228 & (one - tpotb(js)) / (one - tpotb(js) * onembc(l))
230 if(onembc(l) .gt. zero)
231 & taumir(l) = taumir(l)
232 & + (pperp(l) - pres(js)) / bsq(l) * 0.25_dp*tpotb(js)
233 & * (one-onembc(l))*sqrt(tpotb(js)*onembc(l))
234 & / (one + tpotb(js) * onembc(l))
235 & * (15 - 5*tpotb(js) * onembc(l)
236 & - 7*(tpotb(js)*onembc(l))**2
237 & - 3*(tpotb(js)*onembc(l))**3 )
238 & / ((one + tpotb(js)*onembc(l))**2
239 & - (tpotb(js)*onembc(l))**1.5_dp * (5
240 & - (tpotb(js)*onembc(l))**2))
246 taumin = min(taumir(l),taumin)
247 sigmin = min(sigma_an(l),sigmin)
248 bhotc = (2*pperp(l)+ppar(l)-3*pres(js))/bsq(l)
249 bhotmx = max(bhotc,bhotmx)
252 betprx = huge(betprx)
254 betprc = 0.76786580_dp - 0.29708540_dp * tpotb(js)
255 & + 0.054249860_dp * tpotb(js)**2
256 & - 0.0054254849_dp * tpotb(js)**3
257 & + 0.00030947525_dp * tpotb(js)**4
258 & - 9.7144781e-6_dp * tpotb(js)**5
259 & + 1.3844199e-7_dp * tpotb(js)**6
260 & - 1.4328673e-11_dp * (one + tpotb(js))*tpotb(js)**7
261 betprx=min(betprx,betprc)
264 write (nthreed,101) taumin, sigmin, bhotmx
265 101
format(
" taumin=",1pe15.6,
" sigmin=",1pe15.6,
266 &
" peak hot beta=",1pe15.6)
268 betaht = (2*whper+whpar)/(3*wb)
269 betato = (2*wper+wpar)/(3*wb)
272 betadi = betath + betape
273 betaew = betath + 0.5*(betapa+betape)
274 write (nthreed,102) betath,betaht,betato
275 102
format(
" thermal beta =",1pe13.4,
" beta-hot =",1pe13.4,
276 &
" total beta =",1pe13.4)
277 write (nthreed,103) betapa,betape,betaew-betath
278 103
format(
' hot parallel beta =',1pe13.4,
279 &
' hot perpendicular beta =',1pe13.4,
' beta-hot(EW) =',1pe13.4)
280 write (nthreed,104) betadi,betaew,betprx
281 104
format(
' diamagnetic beta =',1pe13.4,
' beta(equal weighting) ='
282 & , 1pe13.4,
' maximum hot perp. beta_c =',1pe10.3)
284 IF (abs(betaht) .GT. 1.e-12_dp)
THEN
286 105
format(/,
' js',8x,
's',9x,
'th. press',3x,
'par. pres',2x,
289 es = real(js - 1.5,dp) / (ns-1)
290 ppprht = pppr(js) - pres(js)
291 parpht = papr(js) - pres(js)
292 WRITE (nthreed,106) js, es, pres(js), parpht, ppprht
295 106
format (i3,1p,1e15.6,1p,3e12.3)
297 END SUBROUTINE mirror_crit