V3FIT
•All Classes Namespaces Files Functions Variables Enumerations Macros Pages
makegrid.f
1 !******************************************************************************
2 ! File makegrid.f
3 !*******************************************************************************
4 !**********************************************************************
5 !** **
6 !** MAIN PROGRAM: Make Grid **
7 !** **
8 !** **
9 !** PROGRAM DESCRIPTION: **
10 !** MAKEGRID reads in a coils-dot file, and generates **
11 !** an MGRID file. **
12 !** **
13 !** Author: Steve Hirshman **
14 !** James Hanson **
15 !** Jonathan Hebert **
16 !** **
17 !** REFERENCES: **
18 !** (1) **
19 !** **
20 !**********************************************************************
21 !-------------------------------------------------------------------------------
22 ! DEPENDENCIES
23 !-------------------------------------------------------------------------------
24 !
25 ! This file uses the following modules:
26 ! stel_kinds
27 ! located in LIBSTELL/Sources/Modules
28 !
29 ! stel_constants
30 ! located in LIBSTELL/Sources/Modules
31 !
32 !
33 !-------------------------------------------------------------------------------
34 ! CHANGE HISTORY
35 !-------------------------------------------------------------------------------
36 !
37 ! See Section V, at end of file
38 !
39 !-------------------------------------------------------------------------------
40 ! USAGE
41 !-------------------------------------------------------------------------------
42 !
43 ! Executing 'xgrid -h' will printout a help message
44 !
45 ! INPUT FILES
46 !
47 !
48 ! OUTPUT FILES
49 !
50 !-------------------------------------------------------------------------------
51 ! COMMENTS
52 !-------------------------------------------------------------------------------
53 !
54 !
55 !*******************************************************************************
56 ! CODE MAKEGRID
57 !
58 ! SECTION I. Main Program
59 !
60 ! SECTION II. Initialization Subroutines
61 ! interactive_input
62 ! namelist_input
63 !
64 ! SECTION III. TASK SUBROUTINES
65 ! task_mgrid
66 ! task_mgrid_rs
67 ! task_circ_tor_grid
68 !
69 ! SECTION IV. AUXILIARY SUBROUTINES
70 ! coil_group_report
71 ! fft_local
72 !
73 ! SECTION V. SUBROUTINES FOR TESTNG
74 !
75 ! SECTION VI. COMMENTS FOR DIFFERENT REVISIONS
76 !*******************************************************************************
77 
78 !*******************************************************************************
79 ! SECTION I. MAIN PROGRAM
80 !*******************************************************************************
81 
82  PROGRAM makegrid
83 !
84 ! THIS CODE (MAKEGRID) GENERATES B-FIELD COMPONENTS ON R,Z, PHI GRID
85 ! AND COMPUTES THE POLOIDAL FLUX AT NOBS OBSERVATION POINTS
86 !
87 ! NOTE TO USER: EXPERIENCE SHOWS THAT A GRID SPACING OF
88 ! DEL R = DEL Z <= 1/50 WORKS WELL. HERE, DEL R = rmax-rmin.
89 ! GRID SPACING MUCH LARGER THAN THIS MAY ADVERSELY EFFECT CONVERGENCE
90 ! OF VACUUM CODE.
91 !
92 ! BOX DIMENSIONS: rmin <= R <= rmax, zmin <= Z <= zmax
93 ! KP = NO. TOROIDAL PLANES/FIELD PERIOD (MUST EQUAL VMEC VALUE)
94 !-------------------------------------------------------------------------------
95 !
96 ! grid dimensions:
97 ! ir = no. radial (r) points in box
98 ! jz = no. z points in box
99 !
100 ! suggest choosing hr == (rmax-rmin)/(ir-1) equal to
101 ! hz == (zmax-zmin)/(jz-1)
102 !
103 !
104 ! THIS CAN BE RUN EITHER FROM THE COMMAND LINE (INTERACTIVELY) OR FROM A COMMAND FILE:
105 !
106 ! xgrid < filename.cmd
107 !
108 ! WHERE the command file (filename.cmd) has entries:
109 ! coils file extension
110 ! stell_sym (T/F)
111 ! rmin value
112 ! rmax value
113 ! zmin value
114 ! zmax value
115 !
116 
117 !
118 ! SET UP GRID DIMENSIONS. USE EITHER INTERACTIVE (OR DRIVER FILE)
119 ! OR COMMAND LINE ARGUMENT LIST
120 !
121 !-------------------------------------------------------------------------------
122 ! CHANGE LOG
123 !-------------------------------------------------------------------------------
124 ! 1-18-2010 JDHe
125 ! Added task structure to makegrid (does not change functionality of old
126 ! code)
127 ! Added namelist input support (usable with all tasks)
128 ! Added help message
129 ! Added task to calculate vector potential on inside of vacuum vessel
130 ! for use with NIMROD (task CIRC_TOR_GRID)
131 ! Previous functionality available as called before (no command line
132 ! inputs prompts for mgrid data, 10 command line inputs runs mgrid with
133 ! those inputs) as well as through the namelist (task MGRID)
134 ! Added 2 arguments to command line input to reflect interactive input
135 !---------------------------------------------------------------------------
136 ! Use statements are followed by variables and subroutines used in this file
137 ! (subroutines are followed by parentheses)
138 !---------------------------------------------------------------------------
139 
140  USE write_mgrid, only: mgrid_ext, mgrid_mode, lstell_sym, &
141  & rmin, rmax, zmin, zmax, kp, ir, jz
142 
143  USE makegrid_global, only: task
144 
145  IMPLICIT NONE
146 !-----------------------------------------------
147 ! L o c a l V a r i a b l e s
148 !-----------------------------------------------
149  INTEGER :: numargs, i
150  CHARACTER(LEN=100) :: arg1
151  CHARACTER(LEN=20) :: task_use
152 
153  CHARACTER(LEN=57), DIMENSION(23), PARAMETER :: help_message=
154  & (/' Program makegrid command line help ',
155  & ' -h Display this message. ',
156  & ' Number of Command Line Arguments: ',
157  & ' 1 (not -h): ',
158  & ' Run makegrid using Namelist input file supplied in ',
159  & ' argument ',
160  & ' 10: ',
161  & ' Use command line arguments as input to makegrid in ',
162  & ' order (only for task MGRID): ',
163  & ' mgrid_ext ',
164  & ' mgrid_mode ',
165  & ' lstell_sym ',
166  & ' rmin ',
167  & ' rmax ',
168  & ' zmin ',
169  & ' zmax ',
170  & ' kp ',
171  & ' ir (optional,default=121) ',
172  & ' jz (optional,default=121) ',
173  & ' 0: ',
174  & ' Get command line input interactively (only for task ',
175  & ' MGRID) ',
176  & ' END HELP MESSAGE '/)
177 
178 !-----------------------------------------------
179 ! Start of Executable Code
180 !-----------------------------------------------
181 
182  CALL getcarg(1, arg1, numargs)
183 
184  SELECT CASE(numargs)
185 
186  CASE(1)
187 ! INTERACTIVE (USE REDIRECTED DRIVER FILE ALSO, XGRID < DRIVER)
188  IF (arg1 .EQ. '-h' .OR. arg1 .EQ. '--h' .OR. arg1 .EQ. '-H'
189  1 .OR. arg1 .EQ. '--H') THEN
190  DO i=1, SIZE(help_message)
191  WRITE(6,*) help_message(i)
192  END DO
193  stop
194  ELSE
195  CALL namelist_input(arg1)
196  END IF
197 
198  CASE(0)
199 ! Only used for task 'MGRID'
200  CALL interactive_input() !sets task='MGRID'
201 
202  CASE(8:10)
203 ! Only used for task 'MGRID'
204  arg1 = adjustl(arg1)
205  mgrid_ext = trim(arg1)
206  CALL getcarg(2, arg1, numargs)
207  IF (arg1(1:1) == 'R' .or. arg1(1:1) == 'r')
208  1 mgrid_mode = 'R'
209  CALL getcarg(3, arg1, numargs)
210  IF (arg1(1:1) == 'Y' .or. arg1(1:1) == 'y' .or.
211  1 arg1(1:1) == 'T' .or. arg1(1:1) == 't')
212  1 lstell_sym = .true.
213  CALL getcarg(4, arg1, numargs)
214  READ (arg1, *) rmin
215  CALL getcarg(5, arg1, numargs)
216  READ (arg1, *) rmax
217  CALL getcarg(6, arg1, numargs)
218  READ (arg1, *) zmin
219  CALL getcarg(7, arg1, numargs)
220  READ (arg1, *) zmax
221  CALL getcarg(8, arg1, numargs)
222  READ (arg1, *) kp
223  IF (numargs .GE. 9) THEN
224  CALL getcarg(9, arg1, numargs)
225  READ (arg1, *) ir
226  END IF
227  IF (numargs .EQ. 10) THEN
228  CALL getcarg(10, arg1, numargs)
229  READ (arg1, *) jz
230  END IF
231  task='MGRID'
232 
233  CASE DEFAULT
234  stop 'Unknown number of arguments, type "xgrid -h" for help.'
235 
236  END SELECT
237 
238 !-----------------------------------------------------------------------
239 ! TASKS:
240 ! 'makegrid' can now run two separate tasks:
241 ! 'B_GRID': find the magnetic field on a grid of points inside the
242 ! plasma volume (for VMEC, V3FIT, etc.)
243 ! 'CIRC_TOR_GRID': find the vector potential on the outer-most surface
244 ! (for NIMROD)
245 !-----------------------------------------------------------------------
246 
247 ! Convert to lower case, avoid simple misinterpretations.
248  task_use = task
249  CALL tolower(task_use)
250 
251  SELECT CASE(trim(adjustl(task_use)))
252 
253  CASE('mgrid')
254 
255  CALL task_mgrid()
256 
257  CASE('mgrid_rs')
258 
259  CALL task_mgrid_rs()
260 
261  CASE('circ_tor_mgrid')
262 
263  CALL task_circ_tor_grid()
264 
265  CASE DEFAULT
266 
267  WRITE(*,*) 'Unknown task ', trim(adjustl(task_use))
268 
269  END SELECT
270 
271  END PROGRAM makegrid
272 
273 
274 !*******************************************************************************
275 ! SECTION II. Initialization Subroutines
276 !*******************************************************************************
277 !
278 !----------------------------------------------------------------------
279 !*******************************************************************************
280 !----------------------------------------------------------------------
281 !
282  SUBROUTINE namelist_input(arg1)
283 ! namelist_input
284 ! Read a namelist input file for either task.
285 
286  USE stel_constants
287 
288  USE write_mgrid, only: mgrid_ext, mgrid_mode, lstell_sym, &
289  & rmin, rmax, zmin, zmax, kp, ir, jz
290 
291  USE makegrid_global, only: task, rmajor, aminor, nphi, ntheta, &
292  & extcur_mgrid, &
293  & cg_shift_1, cg_shift_2, cg_rot_xcent, cg_rot_theta, &
294  & cg_rot_phi, cg_rot_angle, l_rot_coil_center
295 
296  USE sym_check, ONLY: sym_ir, sym_jz, sym_kp, sym_perform_tests
297 
298  IMPLICIT NONE
299 !-----------------------------------------------
300 ! L o c a l V a r i a b l e s
301 !-----------------------------------------------
302 
303  CHARACTER(LEN=100), INTENT(IN) :: arg1
304  INTEGER :: ferror
305  INTEGER, PARAMETER :: iou_nli=12
306 
307 !-----------------------------------------------
308 ! ** Namelist Variables **
309 !
310 ! Used in all tasks:
311 ! task String to select which mode 'makegrid' to run. See comment
312 ! in program 'makegrid' above. (makegrid_global)
313 ! mgrid_ext String to concatenate to all file names. Must be the extension
314 ! of the 'coils-dot' file to be used. (write_mgrid)
315 !
316 ! Used in task 'mgrid' and 'mgrid_rs'
317 ! mgrid_mode Character to choose wheter to run in raw or scaled mode (write_mgrid)
318 ! lstell_sym Logical of whether or not to assume stellarator symmetry (write_mgrid)
319 ! rmin Minimum radial position on grid (write_mgrid)
320 ! rmax Maximum radial position on grid (write_mgrid)
321 ! zmin Minimum vertical position on grid (write_mgrid)
322 ! zmax Maximum vertical position on grid (write_mgrid)
323 ! kp Number of toroidal planes per period (write_mgrid)
324 ! ir Number of radial mesh points (write_mgrid)
325 ! jz Number of vertical mesh points (write_mgrid)
326 !
327 ! Used in task 'mgrid_rs'.
328 ! (All variables are declared in module makegrid_global)
329 ! ( : indicates dimension nextcur_dim - parameter in makegrid_global)
330 ! cg_shift_1(:,3) Vector to shift all the coils. (Before rotation)
331 ! cg_rot_theta(:) Spherical polar angle to specify axis of rotation
332 ! cg_rot_phi(:) Spherical azimuthal angle to specify axis of rotation
333 ! cg_rot_angle(:) Angle to rotate about axis of rotation.
334 ! NB - LEFT HAND convention. Put left thumb along
335 ! axis of rotation, fingers indicate direction of
336 ! positive rotation.
337 ! cg_rot_xcent(:,3) Position of center of rotation
338 ! l_rot_coil_center(:) Logical. True - use current-averaged center of
339 ! coil-group for center of rotation
340 ! False - use position specified in cg_rot_xcent
341 ! for center of rotation
342 ! cg_shift_2(:,3) Vector to shift all the coils. (After rotation)
343 !
344 ! Used in task 'circ_tor_grid':
345 ! rmajor Major radius of circular torus (makegrid_global)
346 ! aminor Minor radius of circular torus (makegrid_global)
347 ! nphi Number of toroidal mesh points (makegrid_global)
348 ! ntheta Number of poloidal mesh points (makegrid_global)
349 ! extcur_mgrid Current in each coil of 'coils-dot' file (makegrid_global)
350 !-----------------------------------------------
351 
352  namelist /mgrid_nli/ task, &
353  & mgrid_ext, mgrid_mode, lstell_sym, &
354  & rmin, rmax, zmin, zmax, kp, ir, jz, &
355  & rmajor, aminor, nphi, ntheta, extcur_mgrid, &
356  & cg_shift_1, cg_shift_2, cg_rot_xcent, cg_rot_theta, &
357  & cg_rot_phi, cg_rot_angle, l_rot_coil_center, &
358  & sym_ir, sym_jz, sym_kp, sym_perform_tests
359 
360 !-----------------------------------------------
361 ! Start of Executable Code
362 !-----------------------------------------------
363 ! Initialize namelist
364 
365  task='MGRID'
366  mgrid_ext='dummy.cth.m.d12c.f5ss'
367  mgrid_mode='R'
368  lstell_sym=.true.
369  rmin=.45
370  rmax=1.05
371  zmin=-.3
372  zmax=.3
373  kp=5
374 ! Commented out below, so that ir and jz have default values as specified in
375 ! module write_mgrid
376 ! ir=20
377 ! jz=20
378  rmajor=1
379  aminor=0.5
380  nphi=3
381  ntheta=3
382  extcur_mgrid=1.
383  sym_ir=5
384  sym_jz=5
385  sym_kp=4
386  sym_perform_tests=.false.
387 
388  cg_shift_1 = zero
389  cg_shift_2 = zero
390  cg_rot_xcent = zero
391  cg_rot_theta = zero
392  cg_rot_phi = zero
393  cg_rot_angle = zero
394  l_rot_coil_center = .true.
395 
396  WRITE(*,*) ' Running makegrid using NLI file ',
397  1 arg1
398  OPEN(iou_nli, file=trim(adjustl(arg1)),status='OLD',iostat=ferror)
399  IF(ferror .NE. 0) THEN
400  WRITE(*,*) 'Could not open NLI file ', arg1
401  WRITE(*,*) 'Open failed with error status ', ferror
402  END IF
403 
404  READ(iou_nli,mgrid_nli)
405 
406  CLOSE(iou_nli)
407 
408  END SUBROUTINE namelist_input
409 !
410 !----------------------------------------------------------------------
411 !*******************************************************************************
412 !----------------------------------------------------------------------
413 !
414  SUBROUTINE interactive_input()
415 ! interactive_input
416 ! For task MGRID, input parameters when prompted. This is for
417 ! backwards compatibility with the earlier version of makegrid.
418 
419  USE write_mgrid, only: mgrid_ext, mgrid_mode, lstell_sym, &
420  & rmin, rmax, zmin, zmax, kp, ir, jz, &
421  & mgrid_file
422 
423  USE makegrid_global, only: task
424 
425  IMPLICIT NONE
426 !-----------------------------------------------
427 ! L o c a l V a r i a b l e s
428 !-----------------------------------------------
429  INTEGER :: numargs, i
430  CHARACTER(LEN=100) :: arg1
431 
432 !-----------------------------------------------
433 ! Start of Executable Code
434 !-----------------------------------------------
435 
436  WRITE (6, 220, advance='no') &
437  & ' Enter extension of "coils" file : '
438  READ (*, *) mgrid_ext
439 
440  WRITE (6, '(a,/,a)', advance='no') &
441  & ' Scale (S) bfield to unit current/turn OR', &
442  & ' use raw (R) currents from coils file: '
443  READ (*, *) mgrid_file
444  IF (mgrid_file(1:1) == 'R' .or. mgrid_file(1:1) == 'r')
445  & mgrid_mode = 'R'
446 
447  WRITE (6, 220, advance='no') &
448  & ' Assume stellarator symmetry (Y/N)? : '
449  READ (*, *) mgrid_file
450  IF (mgrid_file(1:1) == 'Y' .or. mgrid_file(1:1) == 'y')
451  & lstell_sym = .true.
452 
453  WRITE (6, 220, advance='no') &
454  & ' Enter rmin (min radial grid dimension) : '
455  READ (*, *) rmin
456 
457  WRITE (6, 220, advance='no') &
458  & ' Enter rmax (max radial grid dimension) : '
459  READ (*, *) rmax
460 
461  WRITE (6, 220, advance='no') &
462  & ' Enter zmin (min vertical grid dimension): '
463  READ (*, *) zmin
464 
465  WRITE (6, 220, advance='no') &
466  & ' Enter zmax (max vertical grid dimension): '
467  READ (*, *) zmax
468 
469  WRITE (6, 220, advance='no') &
470  & ' Enter number of toroidal planes/period : '
471  READ (*, *) kp
472 
473  WRITE (6, 220, advance='no') &
474  & ' Enter number of r (radial) mesh points : '
475  READ (*, *) ir
476 
477  WRITE (6, 220, advance='no') &
478  & ' Enter number of z mesh points : '
479  READ (*, *) jz
480  task='MGRID'
481 
482  220 FORMAT(a)
483 
484  END SUBROUTINE interactive_input
485 !
486 !----------------------------------------------------------------------
487 !*******************************************************************************
488 !----------------------------------------------------------------------
489 !
490 
491 !*******************************************************************************
492 ! SECTION III. TASK SUBROUTINES
493 !*******************************************************************************
494 !
495 !----------------------------------------------------------------------
496 !*******************************************************************************
497 !----------------------------------------------------------------------
498 !
499  SUBROUTINE task_mgrid()
500 ! task_mgrid
501 ! By this point, the parameters (mgrid_ext, mgrid_mode, lstell_sym, rmin,
502 ! rmax, zmin, zmax, kp, ir, jz) have been input. This subroutine evaluates
503 ! the magnetic field in the plasma volume.
504 
505  USE stel_kinds
506 
507  USE write_mgrid
508  ! , only: mgrid_ext, mgrid_mode, lstell_sym,
509  ! rmin, rmax, zmin, zmax, kp, ir, jz,
510  ! br, bz, bp, kp2, jz2, kp_odd, jz_odd, coil_file, mgrid_file,
511  ! iextc, nfp, nextcur, extcur, fperiod, delr, delp, delz
512  ! subroutines: write_mgrid_nc, cleanup_biotsavart
513  ! access to bsc_b in module bsc.
514 
515  USE biotsavart ! variables nfp_bs, coil_group
516  ! subroutines parse_coils_file
517 
518  USE makegrid_global, only: task
519 
520  USE safe_open_mod !safe_open()
521 
522  USE sym_check, ONLY: init_symmetry, cleanup_symmetry
523 
524  IMPLICIT NONE
525 !-----------------------------------------------
526 ! L o c a l V a r i a b l e s
527 !-----------------------------------------------
528  INTEGER :: istat
529  CHARACTER(LEN=100) :: extcur_file
530  REAL(rprec) :: time_start, time_end1, time_end2
531 
532  REAL(rprec) :: r_ave, z_ave
533  REAL(rprec), DIMENSION(:), ALLOCATABLE :: phi_array
534  INTEGER :: kmax, k, icoll
535 
536 !-----------------------------------------------
537 ! Start of Executable Code
538 !-----------------------------------------------
539 
540 
541  WRITE(*,*) ' Running makegrid with the following parameters:'
542  WRITE(*,*) ' task = ', task
543  WRITE(*,*) ' mgrid_ext = ', mgrid_ext
544  WRITE(*,*) ' mgrid_mode = ', mgrid_mode
545  WRITE(*,*) ' lstell_sym = ', lstell_sym
546  WRITE(*,*) ' rmin = ', rmin
547  WRITE(*,*) ' rmax = ', rmax
548  WRITE(*,*) ' zmin = ', zmin
549  WRITE(*,*) ' zmax = ', zmax
550  WRITE(*,*) ' kp = ', kp
551  WRITE(*,*) ' ir = ', ir
552  WRITE(*,*) ' jz = ', jz
553  WRITE(*,*)
554 
555  IF (rmin .lt. 0.) stop ' rmin must be > 0 in xgrid'
556  IF (rmax .le. rmin) stop ' rmax must be > rmin in xgrid'
557  IF (zmax .le. zmin) stop ' zmax must be > zmin in xgrid'
558  IF (kp .le. 0) stop 'kp must be > 0 in xgrid'
559 
560  ALLOCATE (br(ir,jz,kp), bz(ir,jz,kp), bp(ir,jz,kp), stat=istat)
561  IF (istat .ne. 0) stop ' allocation error in xgrid'
562  ALLOCATE (ar(ir,jz,kp), az(ir,jz,kp), ap(ir,jz,kp), stat=istat)
563  IF (istat .ne. 0) stop ' allocation error in xgrid'
564 
565  IF (lstell_sym) THEN
566  kp2 = kp/2; jz2 = jz/2
567  kp_odd = mod(kp,2)
568  jz_odd = mod(jz,2)
569 !
570 ! Must be sure zmax = -zmin
571 !
572  IF (abs(zmax) > abs(zmin)) THEN
573  zmax = abs(zmax)
574  zmin = -zmax
575  ELSE
576  zmin = -abs(zmin)
577  zmax = -zmin
578  END IF
579  ELSE
580  kp2 = kp; jz2 = jz
581  kp_odd = 0; jz_odd = 0
582  END IF
583 
584  coil_file = 'coils.' // trim(mgrid_ext)
585  mgrid_file = 'mgrid_' // trim(mgrid_ext)
586  extcur_file = 'extcur.' // trim(mgrid_ext)
587  IF (lstell_sym) THEN
588  WRITE (6,*) 'Stellarator symmetry IS assumed'
589  ELSE
590  WRITE (6,*) 'Stellarator symmetry IS NOT assumed'
591  END IF
592  WRITE (6, *) 'rmin = ', rmin,' rmax = ', rmax
593  WRITE (6, *) 'zmin = ', zmin,' zmax = ', zmax
594  WRITE (6, *) 'kp = ', kp, ' ir = ',ir,' jz = ',jz
595  print *
596  WRITE (6, *) 'Input file: ',trim(coil_file)
597  WRITE (6, *) 'Mgrid file: ',trim(mgrid_file)
598  WRITE (6, *) 'Extcur file: ',trim(extcur_file)
599 
600 
601  iextc = 100
602  CALL safe_open(iextc, istat, trim(extcur_file),
603  1 'replace', 'formatted')
604  IF (istat .ne. 0) THEN
605  WRITE (6,*) 'XGRID could not create ', trim(extcur_file)
606  WRITE (6,*) 'IOSTAT = ', istat,' IUNIT = ', iextc
607  stop 25
608  END IF
609 
610 !-----------------------------------------------
611 !
612 ! PARSE FILAMENT FILE FOR NUMBER OF FIELD PERIODS
613 ! SPLIT INTO COIL GROUPS. DETERMINE NEXTCUR
614 ! COMING OUT, IGROUP+100=UNIT NO IS OPENED AND READY TO READ
615 ! AFTER REWINDING
616 !
617  CALL second0(time_start)
618 
619 ! parse_coils_file is in module biotsavart, made available through the
620 ! USE write_mgrid
621  CALL parse_coils_file (coil_file)
622  nfp = nfp_bs
623  nextcur = SIZE(coil_group)
624 
625  ALLOCATE (extcur(nextcur))
626 
627  CALL second0(time_end1)
628 
629  fperiod = (8*atan(one))/nfp
630  delr = (rmax-rmin)/(ir-1)
631  delz = (zmax-zmin)/(jz-1)
632  delp = fperiod/kp
633 
634 #if defined(NETCDF)
635  CALL write_mgrid_nc
636 #else
637  CALL write_mgrid_bin
638 #endif
639 
640  CALL second0(time_end2)
641 
642  WRITE (*, '(2(/,a,f8.3,a))')
643  1 ' TIME IN PARSER = ', time_end1-time_start,' SECONDS',
644  2 ' TIME IN BFIELD = ', time_end2-time_end1,' SECONDS'
645 
646 !-----------------------------------------------
647 ! Print out information about each coil group
648 ! NB coil_group is an allocatable array of bsc_coilcoll declared in module
649 ! biotsavart
650  WRITE(*,'(/a/)') "Extra Information about coil groups:"
651  kmax = max(1,kp2 + kp_odd) ! logic from write_mgrid
652  IF ((kp_odd == 0) .and. lstell_sym) THEN
653  kmax = max(kmax,kp2 + 1)
654  ENDIF
655  ALLOCATE (phi_array(kmax))
656  r_ave = (rmax + rmin) / 2.
657  z_ave = (zmax + zmin) / 2.
658  DO k = 1, kmax
659  phi_array(k) = (k-1) * delp
660  END DO
661  CALL init_symmetry
662  DO icoll = 1, SIZE(coil_group)
663  CALL coil_group_report(coil_group(icoll),icoll,r_ave,z_ave &
664  & ,kmax,phi_array(1:kmax))
665  END DO
666  DEALLOCATE (phi_array)
667  WRITE(*,'(/80("*"))')
668 
669 !-----------------------------------------------
670 ! clean up and finish
671  CALL cleanup_symmetry
672  CALL cleanup_biotsavart
673 
674  DEALLOCATE (br, bp, bz)
675  DEALLOCATE (ar, ap, az)
676 
677  IF (mgrid_mode .eq. 'R') THEN
678  WRITE (6, 205)
679  WRITE (iextc, 205)
680  ELSE
681  WRITE (6, 200) extcur_file
682  END IF
683 
684  CLOSE (iextc)
685 
686  DEALLOCATE (extcur)
687 
688  200 FORMAT(/,
689  1 ' THE BFIELDS HAVE BEEN STORED IN THE MGRID FILE IN SCALED',
690  2 ' MODE. THE EXTERNAL',/,' CURRENTS CORRESPONDING TO THOSE',
691  3 ' IN THE COILS-DOT FILE',/,' ARE GIVEN IN THE EXTCUR ARRAY IN',
692  4 ' THE FILE',1x,a,'. THEY SHOULD BE ENTERED INTO THE VMEC',
693  5 ' INPUT (INDATA) FILE.'/)
694  205 FORMAT(/,
695  1 ' THE BFIELDS HAVE BEEN STORED IN THE MGRID FILE IN RAW MODE.',
696  2 /' THE USER MUST PROVIDE THE EXTCUR ARRAY VALUES FOR THE',
697  3 ' VMEC INPUT (INDATA) FILE.',
698  4 /' CHOOSING EXTCUR=1 CORRESPONDS TO THE CURRENTS PRESENT',
699  5 ' IN THE COILS FILE.')
700 
701 
702  END SUBROUTINE task_mgrid
703 !
704 !----------------------------------------------------------------------
705 !*******************************************************************************
706 !----------------------------------------------------------------------
707 !
708  SUBROUTINE task_mgrid_rs()
709 ! Task: mgrid_rs, generate an MGRID file with Rotated, Shifted coil groups
710 
711  USE stel_kinds
712 
713  USE bsc_t
714  ! Derived type bsc_rs
715  ! subroutines bsc_spill_coil, bsc_mean_r, bsc_construct_rs,
716  ! bsc_rot_shift
717 
718  USE write_mgrid
719  ! , only: mgrid_ext, mgrid_mode, lstell_sym,
720  ! rmin, rmax, zmin, zmax, kp, ir, jz,
721  ! br, bz, bp, kp2, jz2, kp_odd, jz_odd, coil_file, mgrid_file,
722  ! iextc, nfp, nextcur, extcur, fperiod, delr, delp, delz
723  ! subroutines: write_mgrid_nc, cleanup_biotsavart
724 
725  USE biotsavart ! variables nfp_bs, coil_group
726  ! subroutines parse_coils_file
727 
728  USE makegrid_global, only: task, nextcur_dim, &
729  & cg_shift_1, cg_shift_2, cg_rot_xcent, cg_rot_theta, &
730  & cg_rot_phi, cg_rot_angle, l_rot_coil_center
731 
732  USE safe_open_mod !safe_open()
733 
734  USE sym_check, ONLY: init_symmetry, cleanup_symmetry
735 
736  IMPLICIT NONE
737 !-----------------------------------------------
738 ! L o c a l V a r i a b l e s
739 !-----------------------------------------------
740  INTEGER :: istat
741  CHARACTER(LEN=100) :: extcur_file
742  REAL(rprec) :: time_start, time_end1, time_end2
743 
744  INTEGER :: iextcur, icoil
745  REAL(rprec), DIMENSION(3) :: cgxcent, cg_rot_xcent_use, mean_r
746  REAL(rprec), DIMENSION(3) :: zero_a3 = (/zero,zero,zero/)
747  REAL(rprec) :: cur_coil, cur_total
748 
749  REAL(rprec) :: r_ave, z_ave
750  REAL(rprec), DIMENSION(:), ALLOCATABLE :: phi_array
751  INTEGER :: kmax, k, icoll
752 
753  TYPE (bsc_rs) :: this_bsc_rs ! Derived type from bsc, for rotation and shift
754 
755 !-----------------------------------------------
756 ! Start of Executable Code
757 !-----------------------------------------------
758 
759  WRITE(*,*) ' Running makegrid with the following parameters:'
760  WRITE(*,*) ' task = ', task
761  WRITE(*,*) ' mgrid_ext = ', mgrid_ext
762  WRITE(*,*) ' mgrid_mode = ', mgrid_mode
763  WRITE(*,*) ' lstell_sym = ', lstell_sym
764  WRITE(*,*) ' rmin = ', rmin
765  WRITE(*,*) ' rmax = ', rmax
766  WRITE(*,*) ' zmin = ', zmin
767  WRITE(*,*) ' zmax = ', zmax
768  WRITE(*,*) ' kp = ', kp
769  WRITE(*,*) ' ir = ', ir
770  WRITE(*,*) ' jz = ', jz
771  WRITE(*,*)
772 
773  IF (rmin .lt. 0.) stop ' rmin must be > 0 in xgrid'
774  IF (rmax .le. rmin) stop ' rmax must be > rmin in xgrid'
775  IF (zmax .le. zmin) stop ' zmax must be > zmin in xgrid'
776  IF (kp .le. 0) stop 'kp must be > 0 in xgrid'
777 
778  ALLOCATE (br(ir,jz,kp), bz(ir,jz,kp), bp(ir,jz,kp), stat=istat)
779  IF (istat .ne. 0) stop ' allocation error in xgrid'
780  ALLOCATE (ar(ir,jz,kp), az(ir,jz,kp), ap(ir,jz,kp), stat=istat)
781  IF (istat .ne. 0) stop ' allocation error in xgrid'
782 
783  IF (lstell_sym) THEN
784  kp2 = kp/2; jz2 = jz/2
785  kp_odd = mod(kp,2)
786  jz_odd = mod(jz,2)
787 !
788 ! Must be sure zmax = -zmin
789 !
790  IF (abs(zmax) > abs(zmin)) THEN
791  zmax = abs(zmax)
792  zmin = -zmax
793  ELSE
794  zmin = -abs(zmin)
795  zmax = -zmin
796  END IF
797  ELSE
798  kp2 = kp; jz2 = jz
799  kp_odd = 0; jz_odd = 0
800  END IF
801 
802  coil_file = 'coils.' // trim(mgrid_ext)
803  mgrid_file = 'mgrid_' // trim(mgrid_ext)
804  extcur_file = 'extcur.' // trim(mgrid_ext)
805  IF (lstell_sym) THEN
806  WRITE (6,*) 'Stellarator symmetry IS assumed'
807  ELSE
808  WRITE (6,*) 'Stellarator symmetry IS NOT assumed'
809  END IF
810  WRITE (6, *) 'rmin = ', rmin,' rmax = ', rmax
811  WRITE (6, *) 'zmin = ', zmin,' zmax = ', zmax
812  WRITE (6, *) 'kp = ', kp, ' ir = ',ir,' jz = ',jz
813  print *
814  WRITE (6, *) 'Input file: ',trim(coil_file)
815  WRITE (6, *) 'Mgrid file: ',trim(mgrid_file)
816  WRITE (6, *) 'Extcur file: ',trim(extcur_file)
817 
818  iextc = 100
819  CALL safe_open(iextc, istat, trim(extcur_file),
820  1 'replace', 'formatted')
821  IF (istat .ne. 0) THEN
822  WRITE (6,*) 'XGRID could not create ', trim(extcur_file)
823  WRITE (6,*) 'IOSTAT = ', istat,' IUNIT = ', iextc
824  stop 25
825  END IF
826 !
827 ! PARSE FILAMENT FILE FOR NUMBER OF FIELD PERIODS
828 ! SPLIT INTO COIL GROUPS. DETERMINE NEXTCUR
829 ! COMING OUT, IGROUP+100=UNIT NO IS OPENED AND READY TO READ
830 ! AFTER REWINDING
831 !
832  CALL second0(time_start)
833 
834 ! parse_coils_file is in module biotsavart, made available through the
835 ! USE write_mgrid
836  CALL parse_coils_file (coil_file)
837  nfp = nfp_bs
838  nextcur = SIZE(coil_group)
839 
840 ! Test to make sure that nextcur <= nextcur_dim
841  IF(nextcur .gt. nextcur_dim) THEN
842  WRITE(*,*) 'Number of coils greater than default number of ', &
843  & 'currents.'
844  stop
845  END IF
846 
847  ALLOCATE (extcur(nextcur))
848 
849 !-----------------------------------------------
850 ! Rotate and shift the coil groups
851 !-----------------------------------------------
852 ! Coils are stored in array of bsc_coilcoll named coil_group,
853 ! declared in module biotsavart
854 ! Loop over coil groups
855  WRITE(*,*)
856  WRITE(*,*) ' Rotate and Shift of the Coil Groups'
857  DO iextcur = 1,nextcur
858  WRITE(*,*)
859  WRITE(*,*) ' Coil Group ', iextcur,' with s_name ',
860  & coil_group(iextcur) % s_name
861 
862 ! Debug/test - spill the first coil in the coil group. Comment out if unneeded.
863 ! CALL bsc_spill_coil(coil_group(iextcur) % coils(1), &
864 ! & 'coil(1) before rs:')
865 ! WRITE(*,*)
866 
867 ! Compute current-averaged center of coil group (cgxcent)
868  cgxcent(1:3) = zero
869  cur_total = zero
870  DO icoil = 1, coil_group(iextcur) % ncoil
871  cur_coil = coil_group(iextcur) % coils(icoil) % current
872  cur_total = cur_total + cur_coil
873  CALL bsc_mean_r(coil_group(iextcur) % coils(icoil), &
874  & mean_r)
875  cgxcent(1:3) = cgxcent(1:3) + mean_r(1:3) * cur_coil
876  END DO
877  IF (cur_total .ne. 0) cgxcent(1:3) = cgxcent(1:3) / cur_total
878  IF (l_rot_coil_center(iextcur)) THEN
879  cg_rot_xcent_use = cgxcent
880  ELSE
881  cg_rot_xcent_use = cg_rot_xcent(iextcur,1:3)
882  ENDIF
883 
884 ! Generate bsc_rs for first shift, and apply it
885  CALL bsc_construct_rs(this_bsc_rs,zero,zero,zero,zero_a3, &
886  & cg_shift_1(iextcur,1:3))
887  CALL bsc_rot_shift(coil_group(iextcur),this_bsc_rs)
888 
889 ! Generate bsc_rs for rotation and second shift, and apply it
890  CALL bsc_construct_rs(this_bsc_rs,cg_rot_theta(iextcur), &
891  & cg_rot_phi(iextcur),cg_rot_angle(iextcur), &
892  & cg_rot_xcent_use(1:3),cg_shift_2(iextcur,1:3))
893  CALL bsc_rot_shift(coil_group(iextcur),this_bsc_rs)
894 
895 
896  WRITE(*,1000) ' Current-Averaged center of cg = ', &
897  & cgxcent(1:3)
898  WRITE(*,1000) ' First Shift = ', cg_shift_1(iextcur,1:3)
899  WRITE(*,1000) ' Center of Rotation Used = ', &
900  & cg_rot_xcent_use
901  WRITE(*,1000) ' Rotation theta, phi, angle = ', &
902  & cg_rot_theta(iextcur), cg_rot_phi(iextcur), &
903  & cg_rot_angle(iextcur)
904  WRITE(*,1000) ' Second Shift = ', cg_shift_2(iextcur,1:3)
905 1000 FORMAT(a34,3(2x,es12.5))
906 
907 ! Debug/test - spill the first coil in the coil group. Comment out if unneeded.
908 ! CALL bsc_spill_coil(coil_group(iextcur) % coils(1), &
909 ! & 'coil(1) after rs:')
910 
911  END DO
912 ! End loop over coil groups
913 
914  CALL second0(time_end1)
915 
916  fperiod = (8*atan(one))/nfp
917  delr = (rmax-rmin)/(ir-1)
918  delz = (zmax-zmin)/(jz-1)
919  delp = fperiod/kp
920 
921 #if defined(NETCDF)
922  CALL write_mgrid_nc
923 #else
924  CALL write_mgrid_bin
925 #endif
926 
927  CALL second0(time_end2)
928 
929  WRITE (*, '(2(/,a,f8.3,a))')
930  1 ' TIME IN PARSER = ', time_end1-time_start,' SECONDS',
931  2 ' TIME IN BFIELD = ', time_end2-time_end1,' SECONDS'
932 
933 !-----------------------------------------------
934 ! Print out information about each coil group
935 ! NB coil_group is an allocatable array of bsc_coilcoll declared in module
936 ! biotsavart
937  WRITE(*,'(/a/)') "Extra Information about coil groups:"
938  kmax = max(1,kp2 + kp_odd) ! logic from write_mgrid
939  IF ((kp_odd == 0) .and. lstell_sym) THEN
940  kmax = max(kmax,kp2 + 1)
941  ENDIF
942  ALLOCATE (phi_array(kmax))
943  r_ave = (rmax + rmin) / 2.
944  z_ave = (zmax + zmin) / 2.
945  DO k = 1, kmax
946  phi_array(k) = (k-1) * delp
947  END DO
948  CALL init_symmetry
949  DO icoll = 1, SIZE(coil_group)
950  CALL coil_group_report(coil_group(icoll),icoll,r_ave,z_ave &
951  & ,kmax,phi_array(1:kmax))
952  END DO
953  DEALLOCATE (phi_array)
954  WRITE(*,'(/80("*"))')
955 
956 !-----------------------------------------------
957 ! Clean up and finish
958  CALL cleanup_symmetry
959  CALL cleanup_biotsavart
960 
961  DEALLOCATE (br, bp, bz)
962  DEALLOCATE (ar, ap, az)
963 
964  IF (mgrid_mode .eq. 'R') THEN
965  WRITE (6, 205)
966  WRITE (iextc, 205)
967  ELSE
968  WRITE (6, 200) extcur_file
969  END IF
970 
971  CLOSE (iextc)
972 
973  DEALLOCATE (extcur)
974 
975  200 FORMAT(/,
976  1 ' THE BFIELDS HAVE BEEN STORED IN THE MGRID FILE IN SCALED',
977  2 ' MODE. THE EXTERNAL',/,' CURRENTS CORRESPONDING TO THOSE',
978  3 ' IN THE COILS-DOT FILE',/,' ARE GIVEN IN THE EXTCUR ARRAY IN',
979  4 ' THE FILE',1x,a,'. THEY SHOULD BE ENTERED INTO THE VMEC',
980  5 ' INPUT (INDATA) FILE.'/)
981  205 FORMAT(/,
982  1 ' THE BFIELDS HAVE BEEN STORED IN THE MGRID FILE IN RAW MODE.',
983  2 /' THE USER MUST PROVIDE THE EXTCUR ARRAY VALUES FOR THE',
984  3 ' VMEC INPUT (INDATA) FILE.',
985  4 /' CHOOSING EXTCUR=1 CORRESPONDS TO THE CURRENTS PRESENT',
986  5 ' IN THE COILS FILE.')
987 
988 
989  END SUBROUTINE task_mgrid_rs
990 !
991 !----------------------------------------------------------------------
992 !*******************************************************************************
993 !----------------------------------------------------------------------
994 !
995  SUBROUTINE task_circ_tor_grid()
996 ! task_circ_tor_grid
997 ! This subroutine evaluates the vector potential and
998 ! magnetic field on the surface of the (circular) torus described by rmajor
999 ! and aminor at nphi toroidal slices, ntheta times poloidally per slice.
1000 
1001 ! By this point, the parameters (mgrid_ext, rmajor, aminor, nphi, ntheta)
1002 ! have been input.
1003 
1004  USE stel_kinds
1005 
1006  USE stel_constants
1007 
1008  USE write_mgrid, only: mgrid_ext, mgrid_mode, lstell_sym, &
1009  & rmin, rmax, zmin, zmax, kp, ir, jz, &
1010  & coil_file, nextcur
1011 
1012  USE makegrid_global, only: task, rmajor, aminor, nphi, ntheta, &
1013  & extcur_mgrid, nextcur_dim
1014 
1015  USE biotsavart
1016  ! , only: coil_group
1017  ! subroutines: parse_coils_file
1018  ! access to bsc_b in module bsc.
1019 
1020  IMPLICIT NONE
1021 
1022 !-----------------------------------------------
1023 ! L o c a l V a r i a b l e s
1024 !-----------------------------------------------
1025 
1026  REAL(rprec), DIMENSION(3) :: x
1027  REAL(rprec), DIMENSION(3) :: btot, bdum
1028  REAL(rprec) :: tfrac, pfrac, brr, bpp, brout, btheta
1029  REAL(rprec) :: cost, sint, cosp, sinp, norm
1030  REAL(rprec) :: imagr, imagp, imagt, realr, realp, realt
1031  REAL(rprec), DIMENSION(nphi,ntheta,2) :: brf, bpf, bthetaf
1032  CHARACTER(LEN=80) :: ctg_file
1033  INTEGER :: i, j, k, nmode, mmode
1034  INTEGER, PARAMETER :: ctg_iou=17, fp=5
1035  !fp=5, do one field period and multiply phi mode numbers by 5
1036  !fp=1, do all field periods
1037 
1038 !-----------------------------------------------
1039 ! Subroutine Interface
1040 !-----------------------------------------------
1041 
1042  INTERFACE
1043  SUBROUTINE fft_local(dat)
1044  USE stel_kinds
1045  REAL(rprec), DIMENSION(:,:,:) :: dat
1046  END SUBROUTINE fft_local
1047  END INTERFACE
1048 
1049 !-----------------------------------------------
1050 ! Start of Executable Code
1051 !-----------------------------------------------
1052 
1053  WRITE(*,*) ' Running makegrid with the following parameters:'
1054  WRITE(*,*) ' task = ', task
1055  WRITE(*,*) ' mgrid_ext = ', mgrid_ext
1056  WRITE(*,*) ' rmajor = ', rmajor
1057  WRITE(*,*) ' aminor = ', aminor
1058  WRITE(*,*) ' nphi = ', nphi
1059  WRITE(*,*) ' ntheta = ', ntheta
1060  WRITE(*,*)
1061 
1062  coil_file = 'coils.' // trim(mgrid_ext)
1063  WRITE (6, *) 'Input file: ',trim(coil_file)
1064 
1065  ctg_file = 'ctg.'// trim(mgrid_ext)
1066  WRITE(*,*) 'CTG file: ', trim(adjustl(ctg_file))
1067  OPEN(unit=ctg_iou, file=trim(adjustl(ctg_file)))
1068 
1069  CALL parse_coils_file(coil_file)
1070  nextcur = SIZE(coil_group)
1071 
1072  x=(/0.0,0.0,0.0/)
1073 
1074  WRITE(ctg_iou,*) 'Fourier coefficients of magnetic field.'
1075  WRITE(ctg_iou,*) 'This magnetic geometry is 5 fold periodic in ', &
1076  & 'phi. Only modes which preserve this ', &
1077  & 'periodicity are kept.'
1078  WRITE(ctg_iou,*) 'If the field is stellarator symmetric, the ', &
1079  & 'Fourier transform can be expressed in terms ', &
1080  & 'of only sines or only cosines.'
1081  WRITE(ctg_iou,*)
1082  WRITE(ctg_iou,*) 'Major radius: R0 = ', rmajor
1083  WRITE(ctg_iou,*) 'Minor radius: a = ', aminor
1084  WRITE(ctg_iou,*) 'nphi (number of points in phi) = ', nphi
1085  WRITE(ctg_iou,*) 'ntheta (number of points in theta) = ', ntheta
1086  WRITE(ctg_iou,*)
1087  WRITE(ctg_iou,*) 'Coordinate system:'
1088  WRITE(ctg_iou,*) 'r = SQRT((SQRT(x**2+y**2)-R0**2)**2+z**2)'
1089  WRITE(ctg_iou,*) 'THETA = ARCSIN(z/(SQRT(x**2+y**2)-R0)'
1090  WRITE(ctg_iou,*) 'PHI = ARCTAN(y/x)'
1091  WRITE(ctg_iou,*) 'Br = B.r_hat'
1092  WRITE(ctg_iou,*) 'Btheta = B.theta_hat'
1093  WRITE(ctg_iou,*) 'Bphi = B.phi_hat'
1094  WRITE(ctg_iou,*)
1095  WRITE(ctg_iou,*) 'The magnetic field is given by'
1096  WRITE(ctg_iou,*) 'Br(theta,phi) = 1/SQRT(ntheta*nphi) Sum_over_', &
1097  & 'm_and_n (BR_smn * (SIN(m*theta+n*phi))'
1098  WRITE(ctg_iou,*) 'Bphi(theta,phi) = 1/SQRT(ntheta*nphi) Sum_over_' &
1099  & ,'m_and_n (BPHI_cmn * (COS(m*theta+n*phi))'
1100  WRITE(ctg_iou,*) 'Btheta(theta,phi) = 1/SQRT(ntheta*nphi) Sum_', &
1101  & 'over_m_and_n(BTHETA_cmn * (COS(m*theta+n*phi))'
1102  WRITE(ctg_iou,*) 'Where the subscripts "cmn" and"smn" is ', &
1103  & 'shorthand for the trig function involved ("c" ', &
1104  & 'or "s") and the mode number ("m", "n)'
1105 
1106  WRITE(ctg_iou,*)
1107  WRITE(ctg_iou,245) "N","M","BR_smn","BPHI_cmn","BTHETA_cmn"
1108 
1109 ! Test to make sure that nextcur <= nextcur_dim
1110 
1111  IF(nextcur .gt. nextcur_dim) THEN
1112  WRITE(*,*) 'Number of coils greater than default number of ', &
1113  & 'currents.'
1114  stop
1115  END IF
1116 ! nphi, toroidal
1117 ! ntheta, poloidal
1118  DO i=1, nphi
1119  pfrac=(i-1)*1.0/nphi
1120  DO j=1, ntheta
1121  tfrac=(j-1)*1.0/ntheta
1122 
1123 ! Define spatial position to calculate field
1124  cost=cos(twopi*tfrac)
1125  sint=sin(twopi*tfrac)
1126  cosp=cos(twopi*pfrac/fp)
1127  sinp=sin(-twopi*pfrac/fp)
1128  x(1)=(rmajor+aminor*cost)*cosp
1129  x(2)=(rmajor+aminor*cost)*sinp
1130  x(3)=aminor*sint
1131 
1132 ! Calculate field
1133  btot=0
1134  DO k=1, nextcur
1135  CALL bsc_b(coil_group(k),x,bdum)
1136  btot=btot+bdum*extcur_mgrid(k)
1137  END DO
1138 
1139 ! Transform field to local normal coordinates
1140  brr=btot(1)*cosp+btot(2)*sinp
1141  brout=cost*brr+sint*btot(3)
1142  btheta=-sint*brr+cost*btot(3)
1143  bpp=-btot(1)*sinp+btot(2)*cosp
1144 
1145 ! Arrange fields for Fourier transform
1146  brf(i,j,1)=brout
1147  brf(i,j,2)=0
1148  bpf(i,j,1)=bpp
1149  bpf(i,j,2)=0
1150  bthetaf(i,j,1)=btheta
1151  bthetaf(i,j,2)=0
1152 
1153  END DO
1154  END DO
1155 
1156 !Fourier transform the data in theta and phi
1157  CALL fft_local(brf)
1158  CALL fft_local(bpf)
1159  CALL fft_local(bthetaf)
1160 
1161 !Arrange the data so that it runs from negative modes to positive modes and
1162 !write to output file
1163  DO i=nphi/2+1, nphi
1164  DO j=ntheta/2+1, ntheta
1165 
1166  norm=sqrt(nphi*one)*sqrt(ntheta*one)
1167  realr=brf(i,j,1)/norm
1168  imagr=brf(i,j,2)/norm
1169  realp=bpf(i,j,1)/norm
1170  imagp=bpf(i,j,2)/norm
1171  realt=bthetaf(i,j,1)/norm
1172  imagt=bthetaf(i,j,2)/norm
1173 
1174  mmode=j-ntheta-1
1175  nmode=fp*(i-nphi-1)
1176 
1177  WRITE(ctg_iou,240) nmode, mmode, -imagr, realp, realt
1178  END DO
1179  DO j=1, ntheta/2
1180 
1181  norm=sqrt(nphi*one)*sqrt(ntheta*one)
1182  realr=brf(i,j,1)/norm
1183  imagr=brf(i,j,2)/norm
1184  realp=bpf(i,j,1)/norm
1185  imagp=bpf(i,j,2)/norm
1186  realt=bthetaf(i,j,1)/norm
1187  imagt=bthetaf(i,j,2)/norm
1188 
1189  mmode=j-1
1190  nmode=fp*(i-nphi-1)
1191 
1192  WRITE(ctg_iou,240) nmode, mmode, -imagr, realp, realt
1193  END DO
1194  END DO
1195 
1196 
1197  DO i=1, nphi/2
1198  DO j=ntheta/2+1, ntheta
1199 
1200  norm=sqrt(nphi*one)*sqrt(ntheta*one)
1201  realr=brf(i,j,1)/norm
1202  imagr=brf(i,j,2)/norm
1203  realp=bpf(i,j,1)/norm
1204  imagp=bpf(i,j,2)/norm
1205  realt=bthetaf(i,j,1)/norm
1206  imagt=bthetaf(i,j,2)/norm
1207 
1208  mmode=j-ntheta-1
1209  nmode=fp*(i-1)
1210 
1211  WRITE(ctg_iou,240) nmode, mmode, -imagr, realp, realt
1212  END DO
1213  DO j=1, ntheta/2
1214 
1215  norm=sqrt(nphi*one)*sqrt(ntheta*one)
1216  realr=brf(i,j,1)/norm
1217  imagr=brf(i,j,2)/norm
1218  realp=bpf(i,j,1)/norm
1219  imagp=bpf(i,j,2)/norm
1220  realt=bthetaf(i,j,1)/norm
1221  imagt=bthetaf(i,j,2)/norm
1222 
1223  mmode=j-1
1224  nmode=fp*(i-1)
1225 
1226  WRITE(ctg_iou,240) nmode, mmode, -imagr, realp, realt
1227  END DO
1228  END DO
1229 
1230  240 FORMAT(2i6,3es20.12)
1231  245 FORMAT(2a6,3a20)
1232  END SUBROUTINE task_circ_tor_grid
1233 
1234 !*******************************************************************************
1235 ! SECTION IV. AUXILIARY SUBROUTINES
1236 !*******************************************************************************
1237 !----------------------------------------------------------------------
1238 !*******************************************************************************
1239 !----------------------------------------------------------------------
1240 !
1241  SUBROUTINE coil_group_report(cg,igroup,rR,zZ,nphi,phi_array)
1242 ! Subroutine to print out information about a coil group
1243 ! (Stored in a bsc_coilcoll)
1244 
1245 ! JDH 2011-08-15 - First Version
1246 
1247  USE stel_kinds
1248  USE stel_constants
1249 
1250  USE bsc_t
1251  ! Access to derived types, magnetic field computation
1252 
1253  USE sym_check, ONLY: check_symmetry
1254 
1255  IMPLICIT NONE
1256 !-----------------------------------------------
1257 ! A R G U M E N T V a r i a b l e s
1258 !-----------------------------------------------
1259 
1260  TYPE (bsc_coilcoll), INTENT(inout) :: cg
1261  INTEGER, INTENT(in) :: igroup, nphi
1262  REAL(rprec), INTENT(in) :: rR, zZ
1263  REAL(rprec), DIMENSION(1:nphi), INTENT(in) :: phi_array
1264 
1265 ! cg A bsc_coilcoll (hodling a coil group) to report on
1266 ! igroup Coil group number (merely reported here)
1267 ! rR Cylindrical R at which to report B
1268 ! zZ Cylindrical Z at which to report B
1269 ! nphi Number of phi value at which to report B (length of phi_array)
1270 ! phi_array Array of cylindrical phi values at which to report B
1271 !-----------------------------------------------
1272 ! L o c a l V a r i a b l e s
1273 !-----------------------------------------------
1274  REAL(rprec), PARAMETER :: big = 1.e40_rprec
1275  INTEGER, PARAMETER :: int_big = 987654321
1276 
1277  INTEGER :: n_coil_fl, n_coil_fc, n_coil_other, nfil_max_fl, &
1278  & nfil_min_fl, nfil_sum_fl, nfil_sumsq_fl
1279 
1280  INTEGER :: ic, ncoil, nfil, iphi, i
1281 
1282  REAL(rprec) :: cur_max_fl, cur_min_fl, cur_sum_fl, &
1283  & cur_sumsq_fl, cur_max_fc, cur_min_fc, cur_sum_fc, &
1284  & cur_sumsq_fc, ave_nfil_fl, sd_nfil_fl, cur_ave_fl, &
1285  & cur_sd_fl, cur_ave_fc, cur_sd_fc
1286 
1287  TYPE(bsc_coil), POINTER :: coil
1288 
1289  REAL(rprec), DIMENSION(3) :: x, b
1290 
1291  REAL(rprec) :: c, s, br, bphi, bz, cur, phi
1292 
1293 !-----------------------------------------------
1294 ! Start of Executable Code
1295 !-----------------------------------------------
1296 
1297 ! Find out how many coils
1298  ncoil = cg % ncoil
1299  IF (ncoil .le. 0) THEN
1300  WRITE(*,4000) igroup, cg % s_name
1301 4000 FORMAT(/80("*")/"No coils in coil group # ",i4,", group ID ",a30)
1302  RETURN
1303  ELSE
1304  WRITE(*,4100) igroup, cg % s_name
1305 4100 FORMAT(/80("*")/"Coil group # ",i4,", with group ID ",a30)
1306  ENDIF
1307 
1308 ! Initialize variables before loop over coils
1309  n_coil_fl = 0
1310  n_coil_fc = 0
1311  n_coil_other = 0
1312  nfil_max_fl = 0
1313  nfil_min_fl = int_big
1314  nfil_sum_fl = 0
1315  nfil_sumsq_fl = 0
1316  cur_max_fl = - big
1317  cur_min_fl = big
1318  cur_sum_fl = zero
1319  cur_sumsq_fl = zero
1320  cur_max_fc = - big
1321  cur_min_fc = big
1322  cur_sum_fc = zero
1323  cur_sumsq_fc = zero
1324 
1325 ! Loop over coils
1326  DO ic = 1,ncoil
1327  coil => cg % coils(ic)
1328  cur = coil % current
1329 
1330 ! Different coding, depending on c_type
1331  SELECT CASE (coil % c_type)
1332  CASE ('fil_loop','floop')
1333  n_coil_fl = n_coil_fl + 1
1334 
1335  nfil = SIZE(coil % xnod,2) - 1
1336  nfil_max_fl = max(nfil_max_fl,nfil)
1337  nfil_min_fl = min(nfil_min_fl,nfil)
1338  nfil_sum_fl = nfil_sum_fl + nfil
1339  nfil_sumsq_fl = nfil_sumsq_fl + nfil * nfil
1340 
1341  cur_max_fl = max(cur_max_fl,cur)
1342  cur_min_fl = min(cur_min_fl,cur)
1343  cur_sum_fl = cur_sum_fl + cur
1344  cur_sumsq_fl = cur_sumsq_fl + cur * cur
1345 
1346  CASE ('fil_circ','fcirc')
1347  n_coil_fc = n_coil_fc + 1
1348  cur_max_fc = max(cur_max_fc,cur)
1349  cur_min_fc = min(cur_min_fc,cur)
1350  cur_sum_fc = cur_sum_fc + cur
1351  cur_sumsq_fc = cur_sumsq_fc + cur * cur
1352 
1353  CASE DEFAULT
1354  n_coil_other = n_coil_other + 1
1355 
1356  END SELECT
1357  END DO
1358 
1359  ave_nfil_fl = nfil_sum_fl * one / max(n_coil_fl,1)
1360  sd_nfil_fl = sqrt(nfil_sumsq_fl * one / max(n_coil_fl,1) - &
1361  & ave_nfil_fl ** 2)
1362  cur_ave_fl = cur_sum_fl / max(n_coil_fl,1)
1363  cur_sd_fl = sqrt(cur_sumsq_fl / max(n_coil_fl,1) - &
1364  & cur_ave_fl ** 2)
1365  cur_ave_fc = cur_sum_fc / max(n_coil_fc,1)
1366  cur_sd_fc = sqrt(cur_sumsq_fc / max(n_coil_fc,1) - &
1367  & cur_ave_fc ** 2)
1368 
1369 ! Write out results for all coils
1370  WRITE(*,5000)
1371  WRITE(*,5100) n_coil_fl, n_coil_fc, n_coil_other, ncoil
1372 
1373  IF (n_coil_fl .gt. 0) THEN
1374  WRITE(*,5200) nfil_min_fl, nfil_max_fl, ave_nfil_fl, &
1375  & sd_nfil_fl
1376  WRITE(*,5300) cur_min_fl, cur_max_fl, cur_ave_fl, cur_sd_fl
1377  ENDIF
1378 
1379  IF (n_coil_fc .gt. 0) THEN
1380  WRITE(*,5400) cur_min_fc, cur_max_fc, cur_ave_fc, cur_sd_fc
1381  ENDIF
1382 
1383 5000 FORMAT(/" Loops Circles Other Total")
1384 5100 FORMAT(" # coils ",4(3x,i7))
1385 5200 FORMAT(/" Filamentary loops, number of filaments:"/ &
1386  & " Min Max Average SD "/ &
1387  & i8,3x,i8,2(2x,f11.1))
1388 5300 FORMAT(/" Filamentary loops, current:"/ &
1389  & " Min Max Average SD "/ &
1390  & 3(2x,es12.5),2x,es9.2)
1391 5400 FORMAT(/" Filamentary circles, current:"/ &
1392  & " Min Max Average SD "/ &
1393  & 3(2x,es12.5),2x,es9.2)
1394 
1395 ! Loop over phi values
1396  WRITE(*,6000) rr, zz
1397  DO iphi = 1,nphi
1398  phi = phi_array(iphi)
1399  c = cos(phi)
1400  s = sin(phi)
1401  x(1) = rr * c
1402  x(2) = rr * s
1403  x(3) = zz
1404  CALL bsc_b(cg,x,b)
1405  br = b(1) * c + b(2) * s
1406  bphi = - b(1) * s + b(2) * c
1407  bz = b(3)
1408  WRITE(*,6100) phi, br, bphi, bz
1409  END DO
1410 
1411  CALL check_symmetry(igroup)
1412 
1413 6000 FORMAT(/"Magnetic field with unit current multiplier, at R,Z=", &
1414  & 2(2x,es12.5)/t7,"phi",t24,"B . rhat",t38,"B . phihat", &
1415  & t52,"B . zhat")
1416 6100 FORMAT(4(3x,es12.5))
1417 
1418  END SUBROUTINE coil_group_report
1419 !
1420 !----------------------------------------------------------------------
1421 !*******************************************************************************
1422 !----------------------------------------------------------------------
1423 !
1424  SUBROUTINE fft_local(dat)
1425 ! fft_local
1426 ! Subroutine to interface with 'cftfax_g' and 'cfft99', Fourier transform
1427 ! subroutines found in LIBSTELL/FFTpack. 'cftfax_g' conditions the vectors 'trigs' and
1428 ! 'ifax', while 'cfft99' does the actual transform. See those subroutines for more
1429 ! information.
1430 
1431 ! Performs a 2D (two 1D) Fourier transform on the input data. The data
1432 ! should have the form dat(N1, N2, 2) where N1 and N2 are the number of
1433 ! data points in the first and second direction, respectively, and
1434 ! the final index separates the real (1) and imaginary (2) parts of the
1435 ! data.
1436 
1437  USE stel_kinds, only: rprec
1438 
1439  IMPLICIT NONE
1440 !-----------------------------------------------
1441 ! L o c a l V a r i a b l e s
1442 !-----------------------------------------------
1443 
1444  REAL(rprec), DIMENSION(:,:,:) :: dat
1445  REAL(rprec), DIMENSION(:), ALLOCATABLE :: a, trigs, work
1446  INTEGER :: N, i, inc, jump, lot, isign, j
1447  INTEGER, DIMENSION(13) :: ifax
1448 
1449 !-----------------------------------------------
1450 ! Start of Executable Code
1451 !-----------------------------------------------
1452 
1453  n=SIZE(dat,1)
1454  isign=-1 !Sign of exponent in transform
1455  lot=SIZE(dat,2)
1456  jump=n !Number to skip between data vectors. Since we only do
1457  !one transform at a time, it shouldn't matter, but it
1458  !seems to have problems for values larger than N
1459  inc=1 !Number of array spaces between value pairs
1460 
1461  ALLOCATE(a(2*n))
1462  ALLOCATE(trigs(2*n))
1463  ALLOCATE(work(lot*n))
1464 
1465  DO j=1, lot
1466  DO i=1, n
1467 ! The values are stored in a 1-d array with alternating real and imaginary
1468 ! parts.
1469  a(2*i-1)=dat(i,j,1)
1470  a(2*i)=dat(i,j,2)
1471  END DO
1472 
1473 ! Do the transform
1474  CALL cftfax_g(n, ifax, trigs)
1475  CALL cfft99(a,work,trigs, ifax, inc, jump, n, 1, isign)
1476 
1477 ! Put the new values back into the input/output array
1478  DO i=1, n
1479  dat(i,j,1)=a(2*i-1)
1480  dat(i,j,2)=a(2*i)
1481  END DO
1482  END DO
1483 
1484 ! Reallocate to do the transform in the other index
1485  IF(ALLOCATED(a)) DEALLOCATE(a)
1486  ALLOCATE(a(2*lot))
1487  IF(ALLOCATED(trigs)) DEALLOCATE(trigs)
1488  ALLOCATE(trigs(2*lot))
1489  IF(ALLOCATED(work)) DEALLOCATE(work)
1490  ALLOCATE(work(lot*n))
1491 
1492  DO i=1, n
1493  DO j=1, lot
1494  a(2*j-1)=dat(i,j,1)
1495  a(2*j)=dat(i,j,2)
1496  END DO
1497 
1498  CALL cftfax_g(lot, ifax, trigs)
1499  CALL cfft99(a,work,trigs, ifax, inc, jump, lot, 1, isign)
1500 
1501  DO j=1, lot
1502  dat(i,j,1)=a(2*j-1)
1503  dat(i,j,2)=a(2*j)
1504  END DO
1505  END DO
1506 
1507  END SUBROUTINE fft_local
1508 
1509 !*******************************************************************************
1510 ! SECTION VI. COMMENTS FOR DIFFERENT REVISIONS
1511 !*******************************************************************************
1512 !
1513 ! JDH 2011-08-15
1514 ! Added coil_group_report subroutine, and coding to call it.
bsc_t::bsc_b
Definition: bsc_T.f:213
bsc_t::bsc_rot_shift
Definition: bsc_T.f:197
Generated on Thu Mar 5 2020 15:49:25 for V3FIT by   doxygen 1.8.17