V3FIT
fbal.f
1  MODULE fbal
2  USE stel_kinds, ONLY: dp
3  REAL(dp), DIMENSION(:), ALLOCATABLE :: rzu_fac, rru_fac,
4  1 frcc_fac, fzsc_fac
5 
6  CONTAINS
7 
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
18  IMPLICIT NONE
19 !-----------------------------------------------
20  REAL(dp), INTENT(in) :: bsubu(nznt,ns),
21  1 bsubv(nznt,ns)
22 !-----------------------------------------------
23 ! L o c a l V a r i a b l e s
24 !-----------------------------------------------
25  INTEGER :: js, lk
26  INTEGER :: nsmin, nsmax
27 !-----------------------------------------------
28 
29  nsmin=t1lglob; nsmax=t1rglob
30  DO js = nsmin, nsmax
31  buco(js) = sum(bsubu(:,js)*pwint(:,js))
32  bvco(js) = sum(bsubv(:,js)*pwint(:,js))
33  END DO
34 
35  CALL gather1xarray(bvco)
36  CALL gather1xarray(buco)
37 
38 ! FROM AMPERE'S LAW, JcurX are angle averages of jac*JsupX, so
39 ! JcurX = (dV/ds)/twopi**2 <JsupX> where <...> is flux surface average
40  !nsmin=MAX(2,t1lglob); nsmax=MIN(t1rglob,ns-1)
41  nsmin=max(2,tlglob); nsmax=min(trglob,ns-1)
42  DO js = nsmin, nsmax
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)
49  END DO
50  equif(1) = 0
51  equif(ns) = 0
52 
53  !SKS-RANGE: All LHS's computed correctly in [t1lglob, trglob]
54 
55  END SUBROUTINE calc_fbal_par
56 
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
61 #ifdef _ANIMEC
62  3 ,pmap, pd, phot, tpotb, zero
63 #endif
64  USE vmec_params, ONLY: signgs
65  USE vmec_dim, ONLY: ns, nrzt, nznt, ns1
66  USE realspace, ONLY: wint, phip
67 #ifdef _ANIMEC
68  1 ,pperp, ppar, onembc, sigma_an,
69  2 pp1, pp2, pp3
70  USE vforces, gsqrt => azmn_o
71 #endif
72 
73  IMPLICIT NONE
74 !-----------------------------------------------
75  REAL(dp), INTENT(in) :: bsubu(1:nrzt), bsubv(1:nrzt)
76 !-----------------------------------------------
77 ! L o c a l V a r i a b l e s
78 !-----------------------------------------------
79  INTEGER :: js, lk
80  INTEGER :: nsmin, nsmax
81 #ifdef _ANIMEC
82  REAL(dp) :: t4, t5
83 #endif
84 !-----------------------------------------------
85  DO js = 2, ns
86  buco(js) = sum(bsubu(js:nrzt:ns)*wint(js:nrzt:ns))
87  bvco(js) = sum(bsubv(js:nrzt:ns)*wint(js:nrzt:ns))
88  END DO
89 
90 ! FROM AMPERE'S LAW, JcurX are angle averages of jac*JsupX, so
91 ! JcurX = (dV/ds)/twopi**2 <JsupX> where <...> is flux surface average
92  DO js = 2, ns1
93  jcurv(js) = (signgs*ohs)*(buco(js+1) - buco(js))
94  jcuru(js) =-(signgs*ohs)*(bvco(js+1) - bvco(js))
95 !FOR RFP vpphi(js) = (vp(js+1)/phip(js+1) + vp(js)/phip(js))/2
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)
100  END DO
101  equif(1) = 0
102  equif(ns) = 0
103 
104  END SUBROUTINE calc_fbal
105 
106 #ifdef _ANIMEC
107  SUBROUTINE bimax_ppargrad(pp1, pp2, pp3, ppar, onembc, pres, phot,
108  1 tpotb)
109  USE vmec_main, ONLY: zero
110  USE vmec_dim, ONLY: ns, nrzt, nznt, ns1
111 C-----------------------------------------------
112 C D u m m y A r g u m e n t s
113 C-----------------------------------------------
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)
117 C-----------------------------------------------
118 C L o c a l V a r i a b l e s
119 C-----------------------------------------------
120  INTEGER :: js, lk, l
121  REAL(dp), ALLOCATABLE :: tmp0(:), tmp2(:)
122  REAL(dp) :: eps
123 C-----------------------------------------------
124 !********0*********0*********0*********0*********0*********0*********0**
125 ! Model with Anisotropic Pressure from Bi-Maxwellian Distribution. *
126 ! EVALUATION OF AMPLITUDES OF partial p_parallel/partial s AT FIXED B.
127 ! Compute hot particle parallel pressure gradient at fixed B times PP3
128 ! PP3 can be sqrt(g), sqrt(g)/V' or unity
129 !*********0*********0*********0*********0*********0*********0*********0*
130  eps = epsilon(eps)
131  ALLOCATE (tmp0(nrzt), tmp2(nrzt))
132 
133  pp1(1:nrzt:ns) = 0; pp2(1:nrzt:ns) = 0
134  DO js = 2, ns
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)
137  END DO
138 
139  DO l = 1,nrzt-1
140  pp2(l) = 0.5_dp*(tmp2(l)+tmp2(l+1))
141  END DO
142  DO l = 1, nrzt
143  tmp2(l) = pp3(l)*(1-onembc(l))
144  END DO
145 
146  DO js = 2,ns
147  DO lk = 1,nznt
148  l = js +(lk-1)*ns
149  IF (onembc(l) <= zero) THEN
150  pp1(l)= tmp0(l)/
151  & (1 - tpotb(js) * onembc(l)+eps)
152  ELSE
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)
157  END IF
158  END DO
159  END DO
160 
161  pp1 = pp1 * onembc
162  DO l=1, nrzt-1
163  pp1(l) = 0.5_dp*(pp1(l)+pp1(l+1))
164  END DO
165 
166  DEALLOCATE (tmp0, tmp2)
167 
168  END SUBROUTINE bimax_ppargrad
169 
170  SUBROUTINE mirror_crit(taumir, bsq)
171  USE stel_kinds, ONLY: rprec, dp
172  USE realspace, ONLY: sigma_an, pperp, ppar, onembc
173 ! USE vforces, ONLY: bsq => bzmn_o
174  USE vmec_main, ONLY: tpotb, pppr, pmap, pres, papr, vp, wb,
175  & wp, wpar, wper, ns, nznt, nrzt,
176  & zero, one, nthreed
177 !
178 ! WAC (11/30/07): See description of anisotropic pressure below
179 !
180  IMPLICIT NONE
181 C-----------------------------------------------
182 C D u m m y A r g u m e n t s
183 C-----------------------------------------------
184  REAL(dp), DIMENSION(nrzt), INTENT(out) :: taumir
185  REAL(dp), DIMENSION(nrzt), INTENT(inout) :: bsq
186 C-----------------------------------------------
187 C L o c a l P a r a m e t e r s
188 C-----------------------------------------------
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
193 C-----------------------------------------------
194 ! CALCULATION OF MIRROR STABILITY CRITERION
195 !********0*********0*********0*********0*********0*********0*********0**
196 ! *
197 ! Anisotropic Pressure Model *
198 ! specific to case where: *
199 ! a Bi-Maxwellian distribution is considered (by J. Graves) *
200 ! p_parallel(s,B) = pth(1 + phot(s)*H(s,B)) *
201 ! H(s,B)=(B/B_crit)/[1-(T_perp/T_par)(1-B/B_crit)] for B>B_crit *
202 ! For B<B_crit, *
203 ! H(s,B)=H(s,B>B_crit){1-2[(T_perp/T_par)(1-B/B_crit)]^(5/2) / *
204 ! [1+(T_perp/T_par)(1-B/B_crit)]} *
205 ! *
206 !********0*********0*********0*********0*********0*********0*********0**
207 !
208 !********0*********0*********0*********0*********0*********0*********0**
209 ! 1. Change BSQ from total pressure to magnetic pressure *
210 !********0*********0*********0*********0*********0*********0*********0**
211 c
212  bsq = bsq - pperp
213 !
214 !********0*********0*********0*********0*********0*********0*********0**
215 ! 2. Compute Tau-minimum, Sigma-minimum, Peak Hot Particle Beta. *
216 !********0*********0*********0*********0*********0*********0*********0**
217 c
218  whpar = wpar - wp
219  whper = wper - wp
220  bhotmx = -huge(bhotmx)
221  taumin = huge(taumin)
222  sigmin = huge(sigmin)
223 
224  DO 57 js = 2,ns
225  DO 55 lk = 1,nznt
226  l = (lk-1) * ns + js
227  taumir(l) = one + (pperp(l) - pres(js)) / bsq(l)*
228  & (one - tpotb(js)) / (one - tpotb(js) * onembc(l))
229 ! taut = taumir(l) + pperp(l)/bsq(l) *0.25*tpotb(js)
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))
241 55 end do
242 57 end do
243  do 62 js = 2,ns
244  do 60 lk = 1,nznt
245  l = (lk-1) * ns + js
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)
250  60 end do
251  62 end do
252  betprx = huge(betprx)
253  do 64 js=2,ns
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)
262  64 end do
263  bhotmx = bhotmx/3
264  write (nthreed,101) taumin, sigmin, bhotmx
265  101 format(" taumin=",1pe15.6," sigmin=",1pe15.6,
266  & " peak hot beta=",1pe15.6)
267  betath = wp/wb
268  betaht = (2*whper+whpar)/(3*wb)
269  betato = (2*wper+wpar)/(3*wb)
270  betapa = whpar/wb
271  betape = whper/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)
283 
284  IF (abs(betaht) .GT. 1.e-12_dp) THEN
285  write (nthreed,105)
286  105 format(/,' js',8x,'s',9x,'th. press',3x,'par. pres',2x,
287  & 'perp. pres')
288  DO 79 js = 2,ns
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
293  79 END DO
294  END IF
295  106 format (i3,1p,1e15.6,1p,3e12.3)
296 
297  END SUBROUTINE mirror_crit
298 #endif
299 
300  END MODULE fbal
Generated on Thu Mar 5 2020 15:49:26 for V3FIT by   doxygen 1.8.17