3 USE mgrid_mod,
ONLY: vn_nextcur, vn_mgmode, vn_ir,
4 1 vn_jz, vn_kp, vn_nfp, vn_rmin, vn_rmax,
5 2 vn_zmin, vn_zmax, vn_coilgrp, vn_coilcur,
6 3 vn_br0, vn_bz0, vn_bp0, vn_ar0, vn_az0, vn_ap0
11 REAL(rprec),
PARAMETER :: one = 1
12 CHARACTER(LEN=*),
PARAMETER ::
13 1 coildim(2) = (/
'stringsize ',
14 2
'external_coil_groups'/),
15 3 groupdim(1)= (/
'external_coils'/),
16 4 cylcoord(3)= (/
'rad',
'zee',
'phi'/)
20 INTEGER :: ir = 121, jz = 121
21 INTEGER :: iextc, nfp, nextcur, istat
22 INTEGER :: kp, kp2, jz2, kp_odd, jz_odd
23 REAL(rprec),
ALLOCATABLE,
DIMENSION(:, :, :) :: br, bz, bp
24 REAL(rprec),
ALLOCATABLE,
DIMENSION(:, :, :) :: ar, az, ap
25 REAL(rprec),
ALLOCATABLE,
DIMENSION(:) :: extcur
26 REAL(rprec) :: rmin, rmax, zmin, zmax
27 REAL(rprec) :: fperiod, delr, delz, delp
28 LOGICAL :: lstell_sym=.false.
29 CHARACTER(LEN=1) :: mgrid_mode=
'S'
30 CHARACTER(LEN=70) :: mgrid_file, coil_file
31 CHARACTER(LEN=60) :: mgrid_ext
38 SUBROUTINE write_mgrid_nc
45 CHARACTER(LEN=100) :: temp
46 CHARACTER(LEN=100),
ALLOCATABLE,
DIMENSION(:) ::
48 CHARACTER(LEN=100),
ALLOCATABLE,
DIMENSION(:) ::
51 mgrid_file = trim(mgrid_file) //
'.nc'
53 CALL cdf_open(ngrid,mgrid_file,
'w',istat)
54 if (istat .ne. 0) stop
'Error opening WOUT.nc file VMEC WROUT'
59 CALL cdf_define(ngrid, vn_ir, ir)
60 CALL cdf_define(ngrid, vn_jz, jz)
61 CALL cdf_define(ngrid, vn_kp, kp)
62 CALL cdf_define(ngrid, vn_nfp, nfp)
63 CALL cdf_define(ngrid, vn_nextcur, nextcur)
64 CALL cdf_define(ngrid, vn_rmin, rmin)
65 CALL cdf_define(ngrid, vn_zmin, zmin)
66 CALL cdf_define(ngrid, vn_rmax, rmax)
67 CALL cdf_define(ngrid, vn_zmax, zmax)
68 IF (nextcur .eq. 1)
THEN
69 CALL cdf_define(ngrid, vn_coilgrp,coil_group(1)%s_name)
71 CALL cdf_define(ngrid, vn_coilgrp,coil_group(1:nextcur)%s_name,
74 CALL cdf_define(ngrid, vn_mgmode, mgrid_mode)
75 CALL cdf_define(ngrid, vn_coilcur, extcur(1:nextcur),
80 ALLOCATE (vn_br(nextcur), vn_bz(nextcur), vn_bp(nextcur))
81 ALLOCATE (vn_ar(nextcur), vn_az(nextcur), vn_ap(nextcur))
84 write (temp,
'(a,i3.3)')
"_",ig
85 vn_br(ig) = vn_br0 // temp
86 vn_bp(ig) = vn_bp0 // temp
87 vn_bz(ig) = vn_bz0 // temp
88 CALL cdf_define(ngrid, vn_br(ig), br, dimname=cylcoord)
89 CALL cdf_define(ngrid, vn_bp(ig), bp, dimname=cylcoord)
90 CALL cdf_define(ngrid, vn_bz(ig), bz, dimname=cylcoord)
92 vn_ar(ig) = vn_ar0 // temp
93 vn_ap(ig) = vn_ap0 // temp
94 vn_az(ig) = vn_az0 // temp
95 CALL cdf_define(ngrid, vn_ar(ig), ar, dimname=cylcoord)
96 CALL cdf_define(ngrid, vn_ap(ig), ap, dimname=cylcoord)
97 CALL cdf_define(ngrid, vn_az(ig), az, dimname=cylcoord)
104 CALL cdf_write(ngrid, vn_ir, ir)
105 CALL cdf_write(ngrid, vn_jz, jz)
106 CALL cdf_write(ngrid, vn_kp, kp)
107 CALL cdf_write(ngrid, vn_nfp, nfp)
108 CALL cdf_write(ngrid, vn_nextcur, nextcur)
109 CALL cdf_write(ngrid, vn_rmin, rmin)
110 CALL cdf_write(ngrid, vn_zmin, zmin)
111 CALL cdf_write(ngrid, vn_rmax, rmax)
112 CALL cdf_write(ngrid, vn_zmax, zmax)
113 IF (nextcur .eq. 1)
THEN
114 CALL cdf_write(ngrid, vn_coilgrp, coil_group(1)%s_name)
116 CALL cdf_write(ngrid, vn_coilgrp, coil_group(1:nextcur)%s_name)
124 groups:
DO ig = 1,nextcur
126 CALL compute_bfield(ig)
128 CALL cdf_write(ngrid, vn_br(ig), br)
129 CALL cdf_write(ngrid, vn_bp(ig), bp)
130 CALL cdf_write(ngrid, vn_bz(ig), bz)
132 CALL cdf_write(ngrid, vn_ar(ig), ar)
133 CALL cdf_write(ngrid, vn_ap(ig), ap)
134 CALL cdf_write(ngrid, vn_az(ig), az)
138 CALL cdf_write(ngrid, vn_mgmode, mgrid_mode)
139 CALL cdf_write(ngrid, vn_coilcur, extcur(1:nextcur))
141 CALL cdf_close(ngrid)
143 DEALLOCATE (vn_br, vn_bz, vn_bp)
144 DEALLOCATE (vn_ar, vn_az, vn_ap)
146 END SUBROUTINE write_mgrid_nc
149 SUBROUTINE write_mgrid_bin
155 INTEGER,
PARAMETER :: igrid0 = 10
159 INTEGER :: igrid, ig, i
163 mgrid_file = trim(mgrid_file) //
'.bin'
164 CALL safe_open (igrid, istat, trim(mgrid_file),
165 1
'replace',
'unformatted')
167 IF (istat .ne. 0)
THEN
168 WRITE (6,*)
'XGRID could not create ', trim(mgrid_file)
169 WRITE (6,*)
'IOSTAT = ', istat,
' IUNIT = ', igrid
177 WRITE(igrid) ir, jz, kp, nfp, -nextcur
178 WRITE(igrid) rmin, zmin, rmax, zmax
179 WRITE(igrid) (coil_group(i) % s_name, i=1,nextcur)
186 groups:
DO ig = 1,nextcur
188 CALL compute_bfield (ig)
191 WRITE (igrid, iostat=istat) br, bp, bz
192 WRITE (igrid, iostat=istat) ar, ap, az
195 IF (istat .ne. 0) stop
' Error writing bfield components'
199 WRITE (igrid) mgrid_mode
200 WRITE (igrid) extcur(1:nextcur)
208 END SUBROUTINE write_mgrid_bin
211 SUBROUTINE compute_bfield (ig)
213 INTEGER,
INTENT(IN) :: ig
217 INTEGER :: icoil, numcoils, k, j, i, numfils, curindex(1)
218 REAL(rprec) :: rad, zee, phi, extcur_ig, extcur_ig_inv
220 numcoils = coil_group(ig) % ncoil
221 curindex = maxloc(abs(coil_group(ig)%coils(1:numcoils)%current))
222 extcur_ig = coil_group(ig)%coils(curindex(1))%current
224 extcur(ig) = extcur_ig
225 IF (extcur_ig .eq. zero)
WRITE (6, 100) ig
226 100
FORMAT (
'In COILS file, current(s) vanished for coil group ',i4)
229 DO icoil = 1, numcoils
231 & max(
SIZE(coil_group(ig) % coils(icoil) % xnod,2),2) - 1
246 WRITE (6,
'(/,2a)')
' COIL GROUP : ',
247 1 trim(coil_group(ig) % s_name)
248 WRITE (6,
'(a,i6,a,i6)')
' TOTAL COILS IN GROUP: ', numcoils,
249 1
' TOTAL FILAMENTS: ', numfils
256 zee = zmin + (j-1)*delz
258 rad = rmin + (i-1)*delr
259 CALL bfield (rad, phi, zee, br(i,j,k),
260 1 bp(i,j,k), bz(i,j,k), ig)
263 br(i,jz+1-j,k) = -br(i,j,k)
264 bz(i,jz+1-j,k) = bz(i,j,k)
265 bp(i,jz+1-j,k) = bp(i,j,k)
268 CALL afield (rad, phi, zee, ar(i,j,k),
269 1 ap(i,j,k), az(i,j,k), ig)
272 ar(i,jz+1-j,k) = -ar(i,j,k)
273 az(i,jz+1-j,k) = az(i,j,k)
274 ap(i,jz+1-j,k) = ap(i,j,k)
281 1
WRITE (6,
'(a,i4,a,i4,a)')
' K = ',k,
' (OUT OF ',kp,
')'
287 zee = zmin + (j-1)*delz
289 rad = rmin + (i-1)*delr
290 CALL bfield (rad, phi, zee, br(i,j,k),
291 1 bp(i,j,k), bz(i,j,k), ig)
293 br(i,jz+1-j,kp+2-k) = -br(i,j,k)
294 bz(i,jz+1-j,kp+2-k) = bz(i,j,k)
295 bp(i,jz+1-j,kp+2-k) = bp(i,j,k)
298 CALL afield (rad, phi, zee, ar(i,j,k),
299 1 ap(i,j,k), az(i,j,k), ig)
301 ar(i,jz+1-j,kp+2-k) = -ar(i,j,k)
302 az(i,jz+1-j,kp+2-k) = az(i,j,k)
303 ap(i,jz+1-j,kp+2-k) = ap(i,j,k)
307 WRITE (6,
'(a,i4)')
' K = ',k
311 IF ((kp_odd == 0) .and. lstell_sym)
THEN
316 zee = zmin + (j-1)*delz
318 rad = rmin + (i-1)*delr
319 CALL bfield (rad, phi, zee, br(i,j,k),
320 1 bp(i,j,k), bz(i,j,k), ig)
321 br(i,jz+1-j,k) = -br(i,j,k)
322 bz(i,jz+1-j,k) = bz(i,j,k)
323 bp(i,jz+1-j,k) = bp(i,j,k)
325 CALL afield (rad, phi, zee, ar(i,j,k),
326 1 ap(i,j,k), az(i,j,k), ig)
327 ar(i,jz+1-j,k) = -ar(i,j,k)
328 az(i,jz+1-j,k) = az(i,j,k)
329 ap(i,jz+1-j,k) = ap(i,j,k)
333 WRITE (6,
'(a,i4)')
' K = ',k
336 IF (mgrid_mode .eq.
'R')
THEN
339 IF (ig .lt. 10)
WRITE (iextc, 220) ig, extcur_ig, extcur_ig
340 IF (ig .ge. 10)
WRITE (iextc, 225) ig, extcur_ig, extcur_ig
343 IF (abs(extcur_ig) .gt. 1.e-100_rprec)
THEN
344 extcur_ig_inv = 1. / extcur_ig
346 extcur_ig_inv = 1.e100_rprec
349 br = br*extcur_ig_inv; bp = bp*extcur_ig_inv;
350 bz = bz*extcur_ig_inv
352 ar = ar*extcur_ig_inv; ap = ap*extcur_ig_inv;
353 az = az*extcur_ig_inv
355 IF (ig .lt. 10)
WRITE (iextc, 210) ig, extcur_ig
356 IF (ig .ge. 10)
WRITE (iextc, 215) ig, extcur_ig
360 210
FORMAT(
'EXTCUR(', i1,
') = ', 1p,e22.14)
361 215
FORMAT(
'EXTCUR(', i2,
') = ', 1p,e22.14)
362 220
FORMAT(
'EXTCUR(', i1,
') = ', 1p,e22.14,
'/',e22.14)
363 225
FORMAT(
'EXTCUR(', i2,
') = ', 1p,e22.14,
'/',e22.14)
365 END SUBROUTINE compute_bfield
367 END MODULE write_mgrid