1 SUBROUTINE vacuum_par (rmnc, rmns, zmns, zmnc, xm, xn,
2 & plascur, rbtor, wint, ns, ivac_skip,
3 & ivac, mnmax, ier_flag, lscreen)
5 USE vparams,
ONLY: nthreed, zero, one, mu0
6 USE vmec_params,
ONLY: norm_term_flag, phiedge_error_flag
7 USE vmec_input,
ONLY: lrfp
8 USE vmec_main,
ONLY: nznt, irst
9 USE parallel_include_module
15 INTEGER :: ns, ivac_skip, ivac, mnmax, ier_flag
16 REAL(dp) :: plascur, rbtor
17 REAL(dp),
DIMENSION(mnmax),
INTENT(in) ::
18 & rmnc, rmns, zmns, zmnc, xm, xn
19 REAL(dp),
DIMENSION(nuv3),
INTENT(in) :: wint
24 INTEGER :: mn, n, n1, m, i, info=0, j
25 REAL(dp),
DIMENSION(:),
POINTER :: potcos, potsin
26 REAL(dp),
ALLOCATABLE :: potu(:), potv(:)
27 REAL(dp),
ALLOCATABLE :: amatrix(:)
28 REAL(dp):: dn2, dm2, cosmn, sinmn, huv, hvv,
29 & det, bsubuvac, fac, ton, toff
30 REAL(dp) :: tmp1(2), tmp2(2)
43 IF (.NOT.
ALLOCATED(potvac)) stop
'POTVAC not ALLOCATED in VACCUM'
45 ALLOCATE (amatrix(mnpd2*mnpd2), potu(nuv3), potv(nuv3), stat=i)
46 IF (i .NE. 0) stop
'Allocation error in vacuum'
48 potsin => potvac(1:mnpd)
49 potcos => potvac(1+mnpd:)
51 ALLOCATE (bexu(nuv3), bexv(nuv3), bexn(nuv3),
52 & bexni(nuv3), r1b(nuv), rub(nuv3), rvb(nuv3),
53 & z1b(nuv), zub(nuv3), zvb(nuv3), auu(nuv3), auv(nuv3),
54 & avv(nuv3), snr(nuv3), snv(nuv3), snz(nuv3), drv(nuv3),
55 & guu_b(nuv3), guv_b(nuv3), gvv_b(nuv3), rzb2(nuv),
56 & rcosuv(nuv), rsinuv(nuv), stat=i)
57 IF (i .NE. 0) stop
'Allocation error in vacuum'
84 IF (.NOT.
ALLOCATED(tanu))
CALL precal (wint)
85 CALL surface (rmnc, rmns, zmns, zmnc, xm, xn, mnmax)
87 timer_vac(tsurf) = timer_vac(tsurf) + (toff-ton)
90 CALL bextern (plascur, wint, lscreen)
92 timer_vac(tbext) = timer_vac(tbext)+(toff-ton)
97 CALL scalpot (potvac, amatrix, ivac_skip)
99 timer_vac(tscal) = timer_vac(tscal) + (toff-ton)
102 CALL solver (amatrix, potvac, mnpd2, 1, info)
104 timer_vac(tsolver) = timer_vac(tsolver) + (toff-ton)
105 solver_time = solver_time + (toff - ton)
107 IF (info .NE. 0) stop
'Error in solver in VACUUM'
112 potu(nuv3min:nuv3max) = 0; potv(nuv3min:nuv3max) = 0
122 DO i = nuv3min, nuv3max
124 cosmn = potsin(mn)*cosmni(j,mn)/(pi2*pi2*wint(i))
125 potu(i) = potu(i) + dm2*cosmn
126 potv(i) = potv(i) + dn2*cosmn
128 IF (.NOT.lasym) cycle
130 DO i = nuv3min, nuv3max
132 sinmn = potcos(mn)*sinmni(j,mn)/(pi2*pi2*wint(i))
133 potu(i) = potu(i) - dm2*sinmn
134 potv(i) = potv(i) - dn2*sinmn
139 DO i = nuv3min, nuv3max
140 bsubu_sur(i) = potu(i) + bexu(i)
141 bsubv_sur(i) = potv(i) + bexv(i)
142 huv = p5*guv_b(i)*(nfper)
143 hvv = gvv_b(i)*(nfper*nfper)
144 det = one/(guu_b(i)*hvv-huv*huv)
145 bsupu_sur(i) = (hvv*bsubu_sur(i)-huv*bsubv_sur(i))*det
146 bsupv_sur(i) = ((-huv*bsubu_sur(i))+guu_b(i)*bsubv_sur(i))*det
147 bsqvac(i) = p5*(bsubu_sur(i)*bsupu_sur(i)
148 & + bsubv_sur(i)*bsupv_sur(i))
149 brv(i) = rub(i)*bsupu_sur(i) + rvb(i)*bsupv_sur(i)
150 bphiv(i) = r1b(i)*bsupv_sur(i)
151 bzv(i) = zub(i)*bsupu_sur(i) + zvb(i)*bsupv_sur(i)
156 IF (ivac .EQ. 0)
THEN
158 IF (vrank .EQ. 0)
THEN
159 IF (lscreen)
WRITE (*, 200) nfper, mf, nf, nu, nv
160 WRITE (nthreed, 200) nfper, mf, nf, nu, nv
162 200
FORMAT(/,2x,
'In VACUUM, np =',i3,2x,
'mf =',i3,2x,
'nf =',i3,
163 &
' nu =',i3,2x,
'nv = ',i4)
167 DO i = nuv3min, nuv3max
168 bsubuvac = bsubuvac + bsubu_sur(i)*wint(i)
169 bsubvvac = bsubvvac + bsubv_sur(i)*wint(i)
176 CALL mpi_allreduce(tmp1, tmp2, 2, mpi_real8, mpi_sum,
181 allreduce_time = allreduce_time + (toff - ton)
182 bsubuvac = tmp2(1); bsubvvac = tmp2(2)
183 bsubuvac = bsubuvac*signgs*pi2
186 IF (vrank .EQ. 0)
THEN
188 WRITE (*,1000) bsubuvac*fac, plascur*fac, bsubvvac, rbtor
190 WRITE (nthreed, 1000) bsubuvac*fac, plascur*fac, bsubvvac,
193 1000
FORMAT(2x,
'2*pi * a * -BPOL(vac) = ',1p,e10.2,
194 &
' TOROIDAL CURRENT = ',e10.2,/,2x,
'R * BTOR(vac) = ',
195 & e10.2,
' R * BTOR(plasma) = ',e10.2)
200 IF (rbtor*bsubvvac .LT. zero)
THEN
202 IF (vrank .EQ. 0)
THEN
203 IF (lscreen)
WRITE(*,1100)
207 ier_flag = phiedge_error_flag
210 IF (abs((plascur - bsubuvac)/rbtor) .GT. 5.e-2_dp)
THEN
212 IF (vrank .EQ. 0)
THEN
213 IF (lscreen)
WRITE(*,1200)
222 1100
FORMAT(
'lrfp is TRUE. Ignore phiedge sign problem')
223 1200
FORMAT(
'lrfp is TRUE. Proceed with convergence')
226 &
DEALLOCATE (bexu, bexv, bexn, bexni, r1b, rub, rvb, z1b, zub,
227 & zvb, auu, auv, avv, snr, snv, snz, drv, guu_b, guv_b, gvv_b,
228 & rzb2, rcosuv, rsinuv, stat=i)
229 IF (i .NE. 0) stop
'Deallocation error in vacuum'
231 DEALLOCATE (amatrix, potu, potv, stat=i)
232 IF (i .NE. 0) stop
'Deallocation error in vacuum'
237 CALL mpi_allgatherv(mpi_in_place, numjs_vac, mpi_real8, bsqvac,
238 & counts_vac, disps_vac, mpi_real8, vac_comm,
243 timer_vac(tallgv) = timer_vac(tallgv) + (toff-ton)
246 vacuum_time = vacuum_time + (tvacoff - tvacon)
248 END SUBROUTINE vacuum_par
251 SUBROUTINE vacuum(rmnc, rmns, zmns, zmnc, xm, xn,
252 & plascur, rbtor, wint, ns, ivac_skip, ivac,
253 & mnmax, ier_flag, lscreen)
255 USE vparams,
ONLY: nthreed, zero, one, mu0
256 USE vmec_params,
ONLY: norm_term_flag, phiedge_error_flag
257 USE vmec_input,
ONLY: lrfp
263 INTEGER :: ns, ivac_skip, ivac, mnmax, ier_flag
264 REAL(dp) :: plascur, rbtor
265 REAL(dp),
DIMENSION(mnmax),
INTENT(in) ::
266 & rmnc, rmns, zmns, zmnc, xm, xn
267 REAL(dp),
DIMENSION(ns, nuv3),
INTENT(in) :: wint
273 REAL(dp),
ALLOCATABLE :: tmpwint(:)
275 ALLOCATE (tmpwint(nuv3), stat=i)
276 IF (i .NE. 0) stop
'Allocation error in vacuum'
278 tmpwint(:) = wint(ns, :)
280 CALL vacuum_par (rmnc, rmns, zmns, zmnc, xm, xn,
281 & plascur, rbtor, tmpwint, ns, ivac_skip,
282 & ivac, mnmax, ier_flag, lscreen)
284 DEALLOCATE (tmpwint, stat = i)
286 END SUBROUTINE vacuum