V3FIT
write_mgrid.f
1  MODULE write_mgrid
2  USE biotsavart
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
7  IMPLICIT NONE
8 C-----------------------------------------------
9 C L o c a l P a r a m e t e r s
10 C-----------------------------------------------
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'/)
17 C-----------------------------------------------
18 C L o c a l V a r i a b l e s
19 C-----------------------------------------------
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
32 
33  PRIVATE :: istat
34 C-----------------------------------------------
35 
36  CONTAINS
37 #if defined(NETCDF)
38  SUBROUTINE write_mgrid_nc
39  USE ezcdf
40  IMPLICIT NONE
41 C-----------------------------------------------
42 C L o c a l V a r i a b l e s
43 C-----------------------------------------------
44  INTEGER :: ngrid, ig
45  CHARACTER(LEN=100) :: temp
46  CHARACTER(LEN=100), ALLOCATABLE, DIMENSION(:) ::
47  1 vn_br, vn_bp, vn_bz
48  CHARACTER(LEN=100), ALLOCATABLE, DIMENSION(:) ::
49  1 vn_ar, vn_ap, vn_az
50 C-----------------------------------------------
51  mgrid_file = trim(mgrid_file) // '.nc'
52 
53  CALL cdf_open(ngrid,mgrid_file,'w',istat)
54  if (istat .ne. 0) stop 'Error opening WOUT.nc file VMEC WROUT'
55 
56 !
57 ! DEFINE DATA VARIABLES, DIMENSION NAMES
58 !
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)
70  ELSE
71  CALL cdf_define(ngrid, vn_coilgrp,coil_group(1:nextcur)%s_name,
72  1 dimname=coildim)
73  END IF
74  CALL cdf_define(ngrid, vn_mgmode, mgrid_mode)
75  CALL cdf_define(ngrid, vn_coilcur, extcur(1:nextcur),
76  1 dimname=groupdim)
77 !
78 ! STORED AS 3D ARRAYS (ACTUALLY 4D, BUT CUT THROUGH IG)
79 !
80  ALLOCATE (vn_br(nextcur), vn_bz(nextcur), vn_bp(nextcur))
81  ALLOCATE (vn_ar(nextcur), vn_az(nextcur), vn_ap(nextcur))
82 
83  DO ig = 1, 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)
91 
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)
98  END DO
99 
100 
101 !
102 ! WRITE OUT DATA
103 !
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)
115  ELSE
116  CALL cdf_write(ngrid, vn_coilgrp, coil_group(1:nextcur)%s_name)
117  END IF
118 
119 !
120 ! SET UP CYLINDRICAL COMPONENTS OF MAGNETIC FIELD ON GRID
121 ! SUM OVER CURRENT GROUPS IG = 1,NEXTCUR
122 ! NOTE: USER MUST SUPPLY SUBROUTINE " " TO COMPUTE THESE VALUES
123 !
124  groups: DO ig = 1,nextcur
125 
126  CALL compute_bfield(ig)
127 
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)
131 
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)
135 
136  END DO groups
137 
138  CALL cdf_write(ngrid, vn_mgmode, mgrid_mode)
139  CALL cdf_write(ngrid, vn_coilcur, extcur(1:nextcur))
140 
141  CALL cdf_close(ngrid)
142 
143  DEALLOCATE (vn_br, vn_bz, vn_bp)
144  DEALLOCATE (vn_ar, vn_az, vn_ap)
145 
146  END SUBROUTINE write_mgrid_nc
147 #endif
148 
149  SUBROUTINE write_mgrid_bin
150  USE safe_open_mod
151  IMPLICIT NONE
152 C-----------------------------------------------
153 C L o c a l P a r a m e t e r s
154 C-----------------------------------------------
155  INTEGER, PARAMETER :: igrid0 = 10
156 C-----------------------------------------------
157 C L o c a l V a r i a b l e s
158 C-----------------------------------------------
159  INTEGER :: igrid, ig, i
160 C-----------------------------------------------
161 
162  igrid = igrid0
163  mgrid_file = trim(mgrid_file) // '.bin'
164  CALL safe_open (igrid, istat, trim(mgrid_file),
165  1 'replace', 'unformatted')
166 
167  IF (istat .ne. 0) THEN
168  WRITE (6,*) 'XGRID could not create ', trim(mgrid_file)
169  WRITE (6,*) 'IOSTAT = ', istat,' IUNIT = ', igrid
170  stop 20
171  END IF
172 
173 !
174 ! MARK NEW-STYLE CODE WITH NEXTCUR < 0 SO READ_MGRID CAN DISTINGUISH
175 ! FROM OLD-STYLE MGRID FILE
176 !
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)
180 
181 !
182 ! SET UP CYLINDRICAL COMPONENTS OF MAGNETIC FIELD ON GRID
183 ! SUM OVER CURRENT GROUPS IG = 1,NEXTCUR
184 ! NOTE: USER MUST SUPPLY SUBROUTINE "BFIELD" TO COMPUTE THESE VALUES
185 !
186  groups: DO ig = 1,nextcur
187 
188  CALL compute_bfield (ig)
189 
190 
191  WRITE (igrid, iostat=istat) br, bp, bz
192  WRITE (igrid, iostat=istat) ar, ap, az
193 ! WRITE(igrid,'(1p,9e12.4)') (((br(i,j,k), bz(i,j,k), !OLD STYLE
194 ! 1 bp(i,j,k), i=1,ir), j=1,jz),k=1,kp)
195  IF (istat .ne. 0) stop ' Error writing bfield components'
196 
197  END DO groups
198 
199  WRITE (igrid) mgrid_mode
200  WRITE (igrid) extcur(1:nextcur)
201 
202 !
203 ! ADD RECONSTRUCTION/OBSERVATION STUFF HERE
204 !
205 
206  CLOSE (igrid)
207 
208  END SUBROUTINE write_mgrid_bin
209 
210 
211  SUBROUTINE compute_bfield (ig)
212  IMPLICIT NONE
213  INTEGER, INTENT(IN) :: ig
214 C-----------------------------------------------
215 C L o c a l V a r i a b l e s
216 C-----------------------------------------------
217  INTEGER :: icoil, numcoils, k, j, i, numfils, curindex(1)
218  REAL(rprec) :: rad, zee, phi, extcur_ig, extcur_ig_inv
219 C-----------------------------------------------
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
223 
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)
227 
228  numfils = 0
229  DO icoil = 1, numcoils
230  numfils = numfils + &
231  & max(SIZE(coil_group(ig) % coils(icoil) % xnod,2),2) - 1
232 !
233 ! Compute field for unit current
234 !
235 ! JDH 2011-08-16. Comment out below - don't want to change bsc_coil data
236 ! Adjust B values late on with Raw test
237 ! Will behave differently if extcur_ig == 0.
238 ! IF (extcur_ig .ne. zero) THEN
239 ! coil_group(ig) % coils(icoil) % current =
240 ! 1 coil_group(ig) % coils(icoil) % current /extcur_ig
241 ! ELSE
242 ! coil_group(ig) % coils(icoil) % current = 1
243 ! END IF
244  END DO
245 
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
250 
251  k = 1 ! this is always a symmetry plane
252  phi = (k-1)*delp
253 
254 !$omp parallel do private(i, j, zee, rad)
255  DO j=1,jz2 + jz_odd
256  zee = zmin + (j-1)*delz
257  DO i=1,ir
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)
261 
262  IF (lstell_sym) THEN
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)
266  END IF
267 
268  CALL afield (rad, phi, zee, ar(i,j,k),
269  1 ap(i,j,k), az(i,j,k), ig)
270 
271  IF (lstell_sym) THEN
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)
275  END IF
276 
277  END DO
278  END DO
279 !$omp end parallel do
280  IF (kp .ne. 1)
281  1 WRITE (6, '(a,i4,a,i4,a)') ' K = ',k,' (OUT OF ',kp,')'
282 
283 !$omp parallel do private(i, j, k, zee, rad, phi)
284  DO k=2,kp2+kp_odd
285  phi = (k-1)*delp
286  DO j=1,jz
287  zee = zmin + (j-1)*delz
288  DO i=1,ir
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)
292  IF (lstell_sym) THEN
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)
296  END IF
297 
298  CALL afield (rad, phi, zee, ar(i,j,k),
299  1 ap(i,j,k), az(i,j,k), ig)
300  IF (lstell_sym) THEN
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)
304  END IF
305  END DO
306  END DO
307  WRITE (6,'(a,i4)') ' K = ',k
308  END DO
309 !$omp end parallel do
310 
311  IF ((kp_odd == 0) .and. lstell_sym) THEN ! another symmetry plane
312  k = kp2 + 1
313  phi = (k-1)*delp
314 !$omp parallel do private(i, j, zee, rad)
315  DO j=1,jz2 + jz_odd
316  zee = zmin + (j-1)*delz
317  DO i=1,ir
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)
324 
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)
330  END DO
331  END DO
332 !$omp end parallel do
333  WRITE (6,'(a,i4)') ' K = ',k
334  END IF
335 
336  IF (mgrid_mode .eq. 'R') THEN
337 ! JDH 2011-080-16. Comment out below (1 line)
338 ! br = br*extcur_ig; bp = bp*extcur_ig; bz = bz*extcur_ig ! "Raw" fields
339  IF (ig .lt. 10) WRITE (iextc, 220) ig, extcur_ig, extcur_ig ! currents for "raw" mode
340  IF (ig .ge. 10) WRITE (iextc, 225) ig, extcur_ig, extcur_ig
341  ELSE
342 ! JDH 2011-08-16. Add below, scale b. 2011-08-24 - clean up extcur_ig_inv
343  IF (abs(extcur_ig) .gt. 1.e-100_rprec) THEN
344  extcur_ig_inv = 1. / extcur_ig
345  ELSE
346  extcur_ig_inv = 1.e100_rprec
347  ENDIF
348 ! extcur_ig_inv = 1. / MAX(extcur_ig, 1.E-100)
349  br = br*extcur_ig_inv; bp = bp*extcur_ig_inv;
350  bz = bz*extcur_ig_inv ! "Scaled" fields
351 
352  ar = ar*extcur_ig_inv; ap = ap*extcur_ig_inv;
353  az = az*extcur_ig_inv ! "Scaled" fields
354 ! JDH 2011-08-16. end addition
355  IF (ig .lt. 10) WRITE (iextc, 210) ig, extcur_ig ! currents for "scaled" mode
356  IF (ig .ge. 10) WRITE (iextc, 215) ig, extcur_ig
357  END IF
358 
359 
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)
364 
365  END SUBROUTINE compute_bfield
366 
367  END MODULE write_mgrid
Generated on Thu Mar 5 2020 15:49:25 for V3FIT by   doxygen 1.8.17