V3FIT
mgrid_mod.f
1  MODULE mgrid_mod
2  USE v3_utilities
3  USE stel_kinds
4  USE vmec_input, ONLY: nbfld, nflxs, lfreeb, lrecon
5  USE vsvd0, ONLY: nigroup, nparts, npfcoil, nbcoilsp, nfloops,
6  1 nbctotp
7  IMPLICIT NONE
8 
9  INTEGER, PARAMETER :: nlimset = 2 !number of different limiters
10  CHARACTER(LEN=*), PARAMETER ::
11  1 vn_br0 = 'br', vn_bp0 = 'bp', vn_bz0 = 'bz',
12  2 vn_ar0 = 'ar', vn_ap0 = 'ap', vn_az0 = 'az',
13  3 vn_ir = 'ir', vn_jz = 'jz',
14  4 vn_kp = 'kp', vn_nfp = 'nfp',
15  5 vn_rmin='rmin', vn_rmax='rmax', vn_zmin='zmin',
16  6 vn_zmax='zmax', vn_coilgrp='coil_group'
17  CHARACTER(LEN=*), PARAMETER ::
18  1 vn_nextcur = 'nextcur', vn_mgmode='mgrid_mode',
19  2 vn_coilcur = 'raw_coil_cur',
20  3 vn_flp = 'nobser', vn_nobd = 'nobd', vn_nbset = 'nbsets',
21  4 vn_nbfld = 'nbfld',
22  2 ln_flp = 'flux loops', ln_nobd = 'Connected flux loops',
23  3 ln_nbset = 'B-coil loops', ln_next = 'External currents',
24  4 ln_nbfld = 'B-coil measurements'
25 
26 C-----------------------------------------------
27 C L o c a l V a r i a b l e s
28 C-----------------------------------------------
29 ! nr0b, np0b, nz0b:
30 ! : grid dimensions of magnetic field in mgrid file
31 ! nbvac : total number of grid points (nr0b*np0b*nz0b) in mgrid file
32 ! bvac(:,1): br (radial component of external magnetic field)
33 ! bvac(:,2): bp (toroidal component)
34 ! bvac(:,3) = bz (z-component)
35 ! rminb, rmaxb : min (max) radial dimension of grid in mgrid
36 ! zminb, zmaxb : min (max) vertical dimension of grid in mgrid
37 !
38 ! nextcur: no. of EXTERNAL current groups (eg., TF, PF, helical)
39 ! raw_coil_current array of raw currents for each coil group
40 ! mgrid_mode = 'S', scaled mode; = 'R', raw mode
41 ! curlabel: array of labels describing each current group
42 ! included in green''s FUNCTION BFIELD response
43 !
44 ! - - - - - - - - - - - - - - - - - -
45 ! FOR DIAGNOSTICS AND DATA ANALYSIS
46 ! (HERE,COILS ARE FOR MEASURING FIELDS, FLUXES)
47 ! iconnect: two-dimensional array describing electrical
48 ! needflx: =NEEDIT, loop required for flux match
49 ! >=ISYMCOIL, loop required but flux computed by
50 ! invoking symmetry in Z
51 ! =IDONTNEED, loop not required for flux match
52 ! needbfld: =NEEDIT, loop required for B-field match
53 ! =ISAMECOIL, loop at same position as previous loop
54 ! >=ISYMCOIL, loop required but B-field computed
55 ! by invoking symmetry in Z
56 ! =IDONTNEED, loop not required in B-field match
57 ! dsiext: connected flux loop signals due to EXTERNAL coils
58 ! plflux: array of measured (inferred) plasma contrib. to flux loops
59 ! plbfld: array of measured (inferred) plasma contrib. to B-loops
60 ! connection of up to four flux loops. Specifies
61 ! the sign and flux loop number of (up to) four
62 ! connected individual loops (indexing based on
63 ! xobser,zobser arrays).
64 ! nobd: number of connected flux loop measurements
65 ! nobser: number of individual flux loop positions
66 ! nbsets: number of B-coil sets defined in mgrid file
67 ! nbcoils(n): number of bfield coils in each set defined in mgrid file
68 ! nbcoilsn: total number of bfield coils defined in mgrid file
69 ! nbfldn: total number of EXTERNAL bfield measurements used in matching
70 ! bloopnames: array of labels describing b-field sets
71 ! dsilabel: array of labels describing connected flux loops
72 ! xobser: array of flux loop R-positions
73 ! zobser: array of flux loop Z-positions
74 ! rbcoil(m,n): R position of the m-th coil in the n-th set from mgrid file
75 ! zbcoil(m,n): Z position of the m-th coil in the n-th set from mgrid file
76 ! abcoil(m,n): orientation (surface normal wrt R axis; in radians)
77 !
78  INTEGER :: nr0b, np0b, nfper0, nz0b
79  INTEGER :: nobd, nobser, nextcur, nbfldn, nbsets, nbcoilsn
80  INTEGER :: nbvac, nbcoil_max, nlim, nlim_max, nsets,
81  1 nrgrid, nzgrid
82  INTEGER, DIMENSION(:), ALLOCATABLE :: needflx, nbcoils
83  INTEGER, DIMENSION(:), ALLOCATABLE :: limitr, nsetsn
84  INTEGER, DIMENSION(:,:), ALLOCATABLE :: iconnect, needbfld
85  REAL(rprec) :: rminb, zminb, rmaxb, zmaxb, delrb, delzb
86  REAL(rprec) ::rx1, rx2, zy1, zy2, condif
87  REAL(rprec), DIMENSION(:,:), ALLOCATABLE, TARGET :: bvac
88  REAL(rprec), DIMENSION(:,:,:), POINTER :: brvac, bzvac, bpvac
89  REAL(rprec), DIMENSION(:,:), ALLOCATABLE :: unpsiext,
90  1 plbfld, rbcoil, zbcoil, abcoil, bcoil, rbcoilsqr
91  REAL(rprec), DIMENSION(:), ALLOCATABLE :: raw_coil_current
92  REAL(rprec), DIMENSION(:), ALLOCATABLE :: xobser, zobser,
93  1 xobsqr, dsiext, psiext, plflux, b_chi
94  CHARACTER(LEN=300) :: mgrid_path
95  CHARACTER(LEN=300) :: mgrid_path_old = " "
96  CHARACTER(LEN=30), DIMENSION(:), ALLOCATABLE :: curlabel
97  CHARACTER(LEN=15), DIMENSION(:), ALLOCATABLE ::
98  1 dsilabel, bloopnames
99  CHARACTER(LEN=30) :: tokid
100  REAL(rprec), DIMENSION(:,:,:), ALLOCATABLE :: dbcoil, pfcspec
101  REAL(rprec), DIMENSION(:,:), ALLOCATABLE ::
102  1 rlim, zlim, reslim, seplim
103  CHARACTER(LEN=1) :: mgrid_mode
104 
105 #if defined(NETCDF)
106  PRIVATE :: read_mgrid_bin, read_mgrid_nc
107 #else
108  PRIVATE :: read_mgrid_bin
109 #endif
110 
111  CONTAINS
112 
113  SUBROUTINE read_mgrid (mgrid_file, extcur, nv, nfp, lscreen,
114  1 ier_flag, comm)
115  USE system_mod
116  USE mpi_inc
117  IMPLICIT NONE
118 C-----------------------------------------------
119 C D u m m y V a r i a b l e s
120 C-----------------------------------------------
121 !
122 ! mgrid_file: full path to mgrid file
123 ! lscreen : logical controlling output to screen
124 ! ier_flag ; error flag returned to caller
125 ! extcur(n) : external current multiplier for bfield(n) components
126 ! comm : Optional mpi communicator.
127 !
128  INTEGER, INTENT(out) :: ier_flag
129  INTEGER, INTENT(in) :: nv, nfp
130  LOGICAL, INTENT(in) :: lscreen
131  REAL(rprec), INTENT(in) :: extcur(:)
132  CHARACTER(len=*), INTENT(in) :: mgrid_file
133  INTEGER, INTENT(in), OPTIONAL :: comm
134 C-----------------------------------------------
135 C L o c a l P a r a m e t e r s
136 C-----------------------------------------------
137 #if defined(VMS)
138  CHARACTER(LEN=*), PARAMETER :: mgrid_defarea='vmec$:[makegrid]'
139 #else
140  CHARACTER(LEN=*), PARAMETER :: mgrid_defarea='$HOME/vmec/MAKEGRID'
141 #endif
142 C-----------------------------------------------
143 C L o c a l V a r i a b l e s
144 C-----------------------------------------------
145 !
146 ! lgrid_exist : logical set if mgrid file is found in given path
147 !
148  INTEGER :: istat, ii
149  CHARACTER(LEN=200) :: home_dir
150  LOGICAL :: lgrid_exist, lfind
151  INTEGER :: mpi_comm
152 C-----------------------------------------------
153 
154 #if defined(MPI_OPT)
155  mpi_comm = mpi_comm_world
156  IF (PRESENT(comm)) mpi_comm = comm
157 #endif
158 
159  mgrid_path = trim(mgrid_file)
160 
161  IF ((mgrid_path .eq. trim(mgrid_path_old)) .and.
162  1 ALLOCATED(curlabel)) THEN
163  print *,' mgrid file previously parsed!'
164  RETURN
165  END IF
166 
167  INQUIRE (file=mgrid_path,exist=lgrid_exist,iostat=istat)
168  IF (istat.ne.0 .or. .not.lgrid_exist) THEN
169  IF (lscreen) print *,' MGRID FILE NOT FOUND IN SPECIFIED ',
170  1 'PATH: SEARCHING DEFAULT AREA'
171  ii = index(mgrid_file,'/',back=.true.)
172  istat = index(mgrid_defarea, '$HOME')
173  IF (istat .ne. 0) THEN
174  CALL getenv('HOME', home_dir)
175  IF (istat .gt. 1) THEN
176  home_dir = mgrid_defarea(1:istat-1) // trim(home_dir)
177  1 // mgrid_defarea(istat+5:)
178  ELSE
179  home_dir = trim(home_dir) // mgrid_defarea(istat+5:)
180  END IF
181  ELSE
182  home_dir = mgrid_defarea
183  END IF
184  mgrid_path = trim(home_dir) // mgrid_file(ii+1:)
185  INQUIRE (file=mgrid_path,exist=lgrid_exist,iostat=istat)
186  END IF
187 
188  mgrid_path_old = mgrid_path
189 
190  ier_flag = 0
191 
192  IF (lgrid_exist) THEN
193  IF (lscreen) print '(2x,2a)',
194  1 'Opening vacuum field file: ', trim(mgrid_file)
195 !
196 ! Parse mgrid file name, look for .nc extension (netcdf format)
197 !
198  ii = len_trim(mgrid_path) - 2
199  lfind = (mgrid_path(ii:ii+2) == '.nc')
200  IF (lfind) THEN
201 #if defined(NETCDF)
202  CALL read_mgrid_nc (mgrid_path, extcur, nv, nfp,
203  1 ier_flag, lscreen, comm)
204 #else
205  lgrid_exist = .false.
206 #endif
207  ELSE
208  CALL read_mgrid_bin (mgrid_path, extcur, nv, nfp,
209  1 ier_flag, lscreen)
210  END IF
211 
212 !SPH060517 IF (np0b .ne. nv) THEN
213  IF (nv.EQ.0 .OR. mod(np0b, nv).NE.0) THEN
214  print *,' NZETA=',nv,
215  1 ' DOES NOT DIVIDE EVENLY INTO NP0B=',np0b,' IN MGRID FILE'
216  ier_flag = 9
217  ELSE IF (nfper0.ne.nfp) THEN
218  print *,' NFP(READ in) = ',nfp,' DOES NOT AGREE WITH ',
219  1 'NFPER (in vacuum field file) = ',nfper0
220  ier_flag = 9
221  END IF
222 
223  END IF
224 
225  IF (ier_flag .ne. 0) RETURN
226 
227  IF (.not.lgrid_exist .or. ier_flag.ne.0) THEN
228  lfreeb = .false.
229  lrecon = .false.
230  IF (lscreen) THEN
231  print *, ' Error opening/reading mgrid file in dir: ',
232  1 trim(home_dir)
233  print *, ' User must supply vacuum bfield in mgrid to ',
234  1 'run vmec in free-boundary mode!'
235  print *, ' Proceeding to run vmec in',
236  1 ' fixed boundary mode'
237  END IF
238  END IF
239 
240  END SUBROUTINE read_mgrid
241 
242 
243  SUBROUTINE read_mgrid_bin (filename, extcur, nv, nfp, ier_flag,
244  1 lscreen)
245  USE safe_open_mod
246  IMPLICIT NONE
247 C-----------------------------------------------
248 C D u m m y A r g u m e n t s
249 C-----------------------------------------------
250 !
251 ! lstyle_2000 : logical controlling ordering of magnetic field components read in
252 !
253  INTEGER, INTENT(in) :: nv, nfp
254  CHARACTER(LEN=*), INTENT(in) :: filename
255  REAL(rprec), INTENT(in) :: extcur(:)
256 C-----------------------------------------------
257 C L o c a l V a r i a b l e s
258 C-----------------------------------------------
259  REAL(rprec), DIMENSION(:,:,:), ALLOCATABLE ::
260  1 brtemp, bztemp, bptemp
261  INTEGER :: ier_flag, iunit = 50
262  integer :: istat, ig, i, j, n, n1, m, nsets_max, k
263  LOGICAL :: lscreen, lstyle_2000
264 C-----------------------------------------------
265  CALL safe_open(iunit, istat, filename, 'old', 'unformatted')
266  IF (istat .ne. 0) THEN
267  ier_flag = 9
268  RETURN
269  END IF
270 
271  READ (iunit,iostat=istat) nr0b, nz0b, np0b, nfper0, nextcur
272  IF (istat .ne. 0) ier_flag = 9
273 
274  IF (nfper0.NE.nfp .OR. mod(np0b, nv).NE.0) RETURN
275 
276  lstyle_2000 = (nextcur < 0)
277  nextcur = abs(nextcur)
278  READ(iunit,iostat=istat) rminb, zminb, rmaxb, zmaxb
279  IF (istat .ne. 0) ier_flag = 9
280 
281  IF (nextcur .eq. 0) THEN
282  print *,' NEXTCUR = 0 IN READING MGRID FILE'
283  ier_flag = 9
284  ELSE IF (nextcur .gt. nigroup) THEN
285  print *,' NEXTCUR > NIGROUP IN MGRID FILE'
286  ier_flag = 9
287  END IF
288 
289  IF (ier_flag .ne. 0) RETURN
290 
291  ALLOCATE (curlabel(5*(nextcur/5+1)), stat=istat) !MIN of 5 for printing
292  curlabel = " "
293  READ(iunit,iostat=istat) (curlabel(n),n=1,nextcur)
294  IF (istat .ne. 0) THEN
295  print *,' reading mgrid file failed (curlabel)'
296  ier_flag = 9
297  RETURN
298  END IF
299 
300 !
301 ! NOTE: ADD UP CONTRIBUTION TO BVAC DIRECTLY FOR ALL EXTERNAL CURRENT GROUPS
302 
303  nbvac = nr0b*nz0b*nv
304  IF (.NOT. ALLOCATED(bvac)) THEN
305  ALLOCATE (bvac(nbvac,3))
306  ELSE IF (SIZE(bvac,1) .ne. nbvac) THEN
307  DEALLOCATE (bvac); ALLOCATE(bvac(nbvac,3))
308  END IF
309 
310  ALLOCATE (brtemp(nr0b,nz0b,np0b), bptemp(nr0b,nz0b,np0b),
311  1 bztemp(nr0b,nz0b,np0b), stat=istat)
312 
313  IF (istat .ne. 0) THEN
314  print *,' allocation for b-vector storage failed'
315  ier_flag = 9
316  RETURN
317  END IF
318 
319  bvac = 0
320 
321  DO ig = 1,nextcur
322  IF (lstyle_2000) THEN
323  READ(iunit, iostat=istat) brtemp, bptemp, bztemp
324  ELSE
325  READ(iunit, iostat=istat) (((brtemp(i,j,k), bztemp(i,j,k),
326  1 bptemp(i,j,k), i= 1,nr0b),
327  2 j=1,nz0b), k=1,np0b)
328  END IF
329 !
330 ! STORE SUMMED BFIELD (OVER COIL GROUPS) IN BVAC
331 !
332  CALL sum_bfield(bvac(1,1), brtemp, extcur(ig), nv)
333  CALL sum_bfield(bvac(1,2), bptemp, extcur(ig), nv)
334  CALL sum_bfield(bvac(1,3), bztemp, extcur(ig), nv)
335 
336  END DO
337 
338  DEALLOCATE (brtemp, bztemp, bptemp)
339  np0b = nv
340 
341  CALL assign_bptrs(bvac)
342 
343  IF (istat .ne. 0) THEN
344  ier_flag = 9
345  RETURN
346  END IF
347 
348  IF (lstyle_2000) THEN
349  READ (iunit, iostat=istat) mgrid_mode
350  IF (istat .eq. 0) THEN
351  ALLOCATE (raw_coil_current(nextcur))
352  READ (iunit, iostat=istat) raw_coil_current(1:nextcur)
353  IF (istat .ne. 0) mgrid_mode = 'N'
354  END IF
355  ELSE
356  mgrid_mode = 'N' !Old-style, no mode info
357  END IF
358 
359 !
360 ! READ IN EXTERNAL POLOIDAL FLUX, FIELD MEASURMENT
361 ! LOOP COORDINATES AND LABELS
362 !
363  READ(iunit,iostat=istat) nobser, nobd, nbsets
364  IF (istat.ne.0) THEN
365  nobser = 0
366  nobd = 0
367  nbsets = 0
368  IF (lscreen) print *,' No observation data in mgrid data'
369  GOTO 900
370  END IF
371 
372  nbfldn = sum(nbfld(:nbsets))
373  ALLOCATE (nbcoils(nbsets), stat=istat)
374  READ(iunit) (nbcoils(n),n=1,nbsets)
375 
376  nbcoil_max = maxval(nbcoils(:nbsets))
377 
378  ALLOCATE (xobser(nobser), zobser(nobser), dsilabel(nobd),
379  1 iconnect(4,nobser+nobd), unpsiext(nobser,nextcur),
380  2 xobsqr(nobser), needflx(nobser), plflux(nobser+nobd),
381  3 dsiext(nobd), psiext(nobser), bloopnames(nbsets),
382  4 needbfld(nbcoil_max,nbsets), plbfld(nbcoil_max,nbsets),
383  5 rbcoil(nbcoil_max,nbsets), zbcoil(nbcoil_max,nbsets),
384  6 abcoil(nbcoil_max,nbsets), bcoil(nbcoil_max,nbsets),
385  7 rbcoilsqr(nbcoil_max,nbsets), b_chi(nbsets),
386  8 dbcoil(nbcoil_max,nbsets,nextcur), stat = istat)
387  IF (istat .ne. 0) THEN
388  IF (lscreen)
389  1 print *,' allocation error for xobser: istat = ',istat
390  ier_flag = 9
391  RETURN
392  END IF
393 
394  IF (nobser .gt. nfloops) THEN
395  print *, 'NOBSER>NFLOOPS'
396  ier_flag = 9
397  END IF
398  IF (nobd .gt. nfloops) THEN
399  print *, 'NOBD>NFLOOPS'
400  ier_flag = 9
401  END IF
402  IF (nflxs .gt. nfloops) THEN
403  print *, 'NFLXS>NFLOOPS'
404  ier_flag = 9
405  END IF
406  IF (nbfldn .gt. nbctotp) THEN
407  print *, 'NBFLDN>NBCTOTP'
408  ier_flag = 9
409  END IF
410  IF (nbcoil_max .gt. nbcoilsp) THEN
411  print *, 'NBCOIL_max>NBCOILSP'
412  ier_flag = 9
413  END IF
414 
415  IF (ier_flag .ne. 0) RETURN
416 
417  IF (nobser+nobd .gt. 0) iconnect(:4,:nobser+nobd) = 0
418 
419  READ(iunit) (xobser(n), zobser(n),n=1,nobser)
420  READ(iunit) (dsilabel(n),n=1,nobd)
421  READ(iunit) ((iconnect(j,n),j=1,4),n=1,nobd)
422 
423  IF (nbcoil_max.gt.0 .and. nbsets.gt.0) THEN
424  rbcoil(:nbcoil_max,:nbsets) = 0
425  zbcoil(:nbcoil_max,:nbsets) = 0
426  abcoil(:nbcoil_max,:nbsets) = 0
427 
428  DO n=1,nbsets
429  IF (nbcoils(n).gt.0) THEN
430  READ(iunit) n1,bloopnames(n1)
431  READ(iunit)(rbcoil(m,n),zbcoil(m,n),abcoil(m,n),
432  1 m=1,nbcoils(n))
433  ENDIF
434  ENDDO
435 
436  dbcoil(:nbcoil_max,:nbsets,:nextcur) = 0
437  END IF
438  DO ig = 1,nextcur
439  !un-connected coil fluxes
440  READ(iunit) (unpsiext(n,ig),n=1,nobser)
441  DO n = 1,nbsets
442  READ(iunit) (dbcoil(m,n,ig),m=1,nbcoils(n))
443  ENDDO
444  ENDDO
445 
446 !
447 ! READ LIMITER & PROUT PLOTTING SPECS
448 !
449  ALLOCATE (limitr(nlimset), nsetsn(nigroup))
450 
451  READ (iunit,iostat=istat) nlim,(limitr(i),i=1,nlim)
452  IF (istat .ne. 0)then
453  nlim = 0
454  IF (lscreen) print *,' No limiter data in mgrid file'
455  GOTO 900
456  END IF
457 
458  nlim_max = maxval(limitr)
459 
460  IF (nlim .gt. nlimset) THEN
461  print *, 'nlim>nlimset'
462  ier_flag = 9
463  RETURN
464  END IF
465 
466  ALLOCATE( rlim(nlim_max,nlim), zlim(nlim_max,nlim),
467  1 reslim(nlim_max,nlim) ,seplim(nlim_max,nlim),
468  2 stat=istat)
469  IF (istat .ne. 0) THEN
470  print *, 'rlim istat!=0'
471  ier_flag = 9
472  RETURN
473  END IF
474 
475  READ(iunit, iostat=istat)
476  1 ((rlim(i,j),zlim(i,j),i=1,limitr(j)),j=1,nlim)
477  READ(iunit, iostat=istat) nsets,(nsetsn(i), i=1,nsets)
478 
479  IF (nsets .gt. nigroup) THEN
480  print *, 'nsets>nigroup'
481  ier_flag = 9
482  RETURN
483  ELSE IF (istat .ne. 0) THEN
484  ier_flag = 9
485  RETURN
486  END IF
487 
488  nsets_max = maxval(nsetsn)
489 
490  IF (nsets_max .gt. npfcoil) THEN
491  print *, 'nsetsn>npfcoil'
492  ier_flag = 9
493  RETURN
494  END IF
495 
496  ALLOCATE (pfcspec(nparts,nsets_max,nsets), stat=istat)
497 
498 ! NOTE TO RMW: SHOULD READ IN NPARTS HERE (PUT INTO MGRID FILE)
499 
500  READ(iunit, iostat=istat) (((pfcspec(i,j,k),i=1,nparts),
501  1 j=1,nsetsn(k)), k=1,nsets)
502 
503  DEALLOCATE (limitr, nsetsn)
504 
505  READ(iunit, iostat=istat) rx1,rx2,zy1,zy2,condif,
506  1 nrgrid,nzgrid,tokid
507 
508  IF (istat .ne. 0) THEN
509  ier_flag = 9
510  RETURN
511  END IF
512 
513  IF (nobser .gt. 0) xobsqr(:nobser) = sqrt(xobser(:nobser))
514 !
515 ! PARTITION MGRID B-LOOPS INTO SETS
516 !
517  nbcoilsn = sum(nbcoils(:nbsets))
518 
519  DO n = 1,nbsets
520  rbcoilsqr(:nbcoils(n),n) = sqrt(rbcoil(:nbcoils(n),n))
521  ENDDO
522 
523  900 CONTINUE
524 
525  CLOSE (iunit)
526 
527  delrb = (rmaxb-rminb)/(nr0b-1)
528  delzb = (zmaxb-zminb)/(nz0b-1)
529 
530 !
531 ! SUM UP CONTRIBUTIONS FROM INDIVIDUAL COIL GROUPS
532 !
533  IF (lfreeb) THEN
534  IF (nobser .gt. 0) psiext(:nobser) = 0
535  IF (nbcoil_max.gt.0 .and. nbsets.gt.0)
536  1 bcoil(:nbcoil_max, :nbsets) = 0
537 
538  DO ig = 1,nextcur
539  IF (nobser .gt. 0)
540  1 psiext(:nobser) = psiext(:nobser) +
541  2 extcur(ig)*unpsiext(:nobser,ig)
542  DO n=1,nbsets
543  n1 = nbcoils(n)
544  bcoil(:n1,n) = bcoil(:n1,n) +
545  1 extcur(ig)*dbcoil(:n1,n,ig)
546  ENDDO
547  ENDDO
548  ENDIF !!IF LFREEB
549 
550  END SUBROUTINE read_mgrid_bin
551 
552 
553 #if defined(NETCDF)
554 !
555 ! PARALLEL MPI MODIFICATIONS ADDED BY Mark R. Cianciosa <cianciosamr@ornl.gov>, 092315
556 !
557  SUBROUTINE read_mgrid_nc (filename, extcur, nv, nfp,
558  1 ier_flag, lscreen, comm)
559  USE ezcdf
560  USE mpi_inc
561  IMPLICIT NONE
562 C-----------------------------------------------
563 C D u m m y A r g u m e n t s
564 C-----------------------------------------------
565  CHARACTER(LEN=*), INTENT(in) :: filename
566  INTEGER, INTENT(in) :: nv, nfp
567  REAL(rprec), INTENT(in) :: extcur(:)
568  INTEGER, INTENT(in) :: comm
569 C-----------------------------------------------
570 C L o c a l V a r i a b l e s
571 C-----------------------------------------------
572  REAL(rprec), DIMENSION(:,:,:), ALLOCATABLE ::
573  1 brtemp, bztemp, bptemp
574  INTEGER :: ier_flag, ngrid
575  INTEGER :: istat, ig
576  LOGICAL :: lscreen
577  INTEGER, DIMENSION(3) :: dimlens
578  CHARACTER(LEN=100) :: temp
579 #if defined(MPI_OPT)
580  INTEGER :: mpi_rank, mpi_size, lMPIInit, MPI_ERR
581 
582  CALL mpi_initialized(lmpiinit, mpi_err)
583  IF (lmpiinit.NE.0) THEN
584  CALL mpi_comm_rank(comm, mpi_rank, istat)
585  CALL mpi_comm_size(comm, mpi_size, istat)
586  ELSE
587  mpi_rank = 0; mpi_size = 1
588  END IF
589 #endif
590 C-----------------------------------------------
591  call cdf_open(ngrid, filename,'r', istat)
592  IF (istat .ne. 0) THEN
593  ier_flag = 9
594  RETURN
595  END IF
596 
597 !
598 ! READ IN DATA
599 !
600  CALL cdf_read(ngrid, vn_ir, nr0b)
601  CALL cdf_read(ngrid, vn_jz, nz0b)
602  CALL cdf_read(ngrid, vn_kp, np0b)
603  CALL cdf_read(ngrid, vn_nfp, nfper0)
604 
605  IF (nfper0.NE.nfp .OR. mod(np0b, nv).NE.0) RETURN
606 
607  CALL cdf_read(ngrid, vn_nextcur, nextcur)
608 
609  IF (nextcur .eq. 0) THEN
610  print *,' NEXTCUR = 0 IN READING MGRID FILE'
611  ier_flag = 9
612  RETURN
613  ELSE IF (nextcur .gt. nigroup) THEN
614  print *,' NEXTCUR > NIGROUP IN MGRID FILE'
615  ier_flag = 9
616  RETURN
617  END IF
618 
619  CALL cdf_read(ngrid, vn_rmin, rminb)
620  CALL cdf_read(ngrid, vn_zmin, zminb)
621  CALL cdf_read(ngrid, vn_rmax, rmaxb)
622  CALL cdf_read(ngrid, vn_zmax, zmaxb)
623 
624  delrb = (rmaxb-rminb)/(nr0b-1)
625  delzb = (zmaxb-zminb)/(nz0b-1)
626 
627  CALL cdf_inquire(ngrid, vn_coilgrp, dimlens)
628  IF (.NOT. ALLOCATED(curlabel)) THEN
629  ALLOCATE (curlabel(nextcur), stat=istat)
630  ELSE IF (SIZE(curlabel) .ne. nextcur) THEN
631  DEALLOCATE (curlabel)
632  ALLOCATE (curlabel(nextcur), stat=istat)
633  END IF
634 !THIS IS A GLITCH WITH cdf_read: must distinguish 1D char array from multi-D
635  IF (nextcur .eq. 1) THEN
636  IF (istat .eq. 0)
637  1 CALL cdf_read(ngrid, vn_coilgrp, curlabel(1))
638  ELSE IF (istat .eq. 0) THEN
639  CALL cdf_read(ngrid, vn_coilgrp, curlabel(1:nextcur))
640  END IF
641 
642  IF (istat .ne. 0) stop 'Error allocating CURLABEL in mgrid_mod'
643 !
644 ! READ 3D Br, Bp, Bz ARRAYS FOR EACH COIL GROUP
645 !
646  nbvac = nr0b*nz0b*nv
647  IF (.NOT. ALLOCATED(bvac)) THEN
648  ALLOCATE (bvac(nbvac,3), stat=istat)
649  ELSE IF (SIZE(bvac,1) .ne. nbvac) THEN
650  DEALLOCATE (bvac); ALLOCATE(bvac(nbvac,3), stat=istat)
651  END IF
652  IF (istat .ne. 0) stop 'Error allocating bvac in mgrid_mod'
653 
654  bvac = 0
655 
656  ALLOCATE (brtemp(nr0b, nz0b, np0b), &
657  & bptemp(nr0b, nz0b, np0b), &
658  & bztemp(nr0b, nz0b, np0b), stat=istat)
659  IF (istat .ne. 0)stop 'Error allocating bXtemp in mgrid_mod'
660  brtemp = 0
661  bptemp = 0
662  bztemp = 0
663 
664 #if defined(MPI_OPT)
665  DO ig = mpi_rank + 1, nextcur, mpi_size
666 #else
667  DO ig = 1, nextcur
668 #endif
669  WRITE (temp, 1000) vn_br0, ig
670  CALL cdf_read(ngrid, temp, brtemp)
671 
672  WRITE (temp, 1000) vn_bp0, ig
673  CALL cdf_read(ngrid, temp, bptemp)
674 
675  WRITE (temp, 1000) vn_bz0, ig
676  CALL cdf_read(ngrid, temp, bztemp)
677 
678 !
679 ! STORE SUMMED BFIELD (OVER COIL GROUPS) IN BVAC
680 !
681  CALL sum_bfield(bvac(1,1), brtemp, extcur(ig), nv)
682  CALL sum_bfield(bvac(1,2), bptemp, extcur(ig), nv)
683  CALL sum_bfield(bvac(1,3), bztemp, extcur(ig), nv)
684  END DO
685 
686  np0b = nv
687 
688 #if defined(MPI_OPT)
689  IF (lmpiinit.NE.0) THEN
690  CALL mpi_allreduce(mpi_in_place, bvac, SIZE(bvac), mpi_real8, &
691  & mpi_sum, comm, istat)
692  CALL assert_eq(istat,0,'MPI_ALLREDUCE failed in read_mgrid_nc')
693  END IF
694 #endif
695 
696 !
697 ! MUST ADD EXTERNAL LOOP STUFF LATER
698 ! MAY DECIDE TO WRITE THAT INFO INTO A SEPARATE FILE
699 ! FOR NOW, JUST SET DEFAULTS
700 !
701  nobser = 0
702  nobd = 0
703  nbsets = 0
704 
705  CALL cdf_inquire(ngrid, vn_mgmode, dimlens, ier=istat)
706  IF (istat .eq. 0) THEN
707  CALL cdf_read(ngrid, vn_mgmode, mgrid_mode)
708  ELSE
709  mgrid_mode = 'N'
710  END IF
711 
712  CALL cdf_inquire(ngrid, vn_coilcur, dimlens, ier=istat)
713  IF (istat .eq. 0) THEN
714  IF (ALLOCATED(raw_coil_current)) DEALLOCATE(raw_coil_current)
715  ALLOCATE (raw_coil_current(nextcur), stat=istat)
716  IF (istat .ne. 0) stop 'Error allocating RAW_COIL in mgrid_mod'
717  CALL cdf_read(ngrid, vn_coilcur, raw_coil_current)
718  END IF
719 
720  CALL cdf_close(ngrid)
721 
722  IF (ALLOCATED(brtemp))
723  1 DEALLOCATE (brtemp, bptemp, bztemp)
724 
725  CALL assign_bptrs(bvac)
726 
727 1000 FORMAT(a,'_',i3.3)
728 
729  END SUBROUTINE read_mgrid_nc
730 #endif
731 
732  SUBROUTINE sum_bfield(bfield, bf_add, cur, nv)
733  INTEGER, INTENT(IN) :: nv
734  REAL(rprec), INTENT(INOUT) :: bfield(nr0b*nz0b,nv)
735  REAL(rprec), INTENT(IN) :: bf_add(nr0b*nz0b,np0b)
736  INTEGER :: nskip
737  REAL(rprec) :: cur
738 
739  nskip = np0b/nv
740  bfield = bfield + cur*bf_add(:,1:np0b:nskip)
741 
742  END SUBROUTINE sum_bfield
743 
744  SUBROUTINE assign_bptrs(bptr)
745  IMPLICIT NONE
746  REAL(rprec), TARGET, INTENT(in) :: bptr(nr0b,nz0b,np0b,3)
747 
748  brvac => bptr(:,:,:,1)
749  bpvac => bptr(:,:,:,2)
750  bzvac => bptr(:,:,:,3)
751 
752  END SUBROUTINE assign_bptrs
753 
754  SUBROUTINE free_mgrid (istat)
755  INTEGER :: istat
756 
757  istat = 0
758 
759  IF (ALLOCATED(bvac)) DEALLOCATE (bvac,stat=istat)
760  IF (ALLOCATED(xobser))
761  1 DEALLOCATE (xobser, xobsqr, zobser, unpsiext, dsiext,
762  2 psiext,plflux, iconnect, needflx, needbfld, plbfld,
763  3 nbcoils, rbcoil, zbcoil, abcoil, bcoil, rbcoilsqr, dbcoil,
764  4 pfcspec,dsilabel, bloopnames, curlabel, b_chi, stat=istat)
765  IF (ALLOCATED(raw_coil_current)) DEALLOCATE(raw_coil_current)
766 
767  IF (ALLOCATED(rlim))
768  1 DEALLOCATE (rlim,zlim, reslim,seplim,stat=istat)
769 
770 ! Reset mgrid_path_old, so that can reread an mgrid file. SL, JDH 2012-07-16
771  mgrid_path_old = " "
772 
773  END SUBROUTINE free_mgrid
774 
775  END MODULE mgrid_mod
system_mod::getenv
Definition: system_mod.f:16
v3_utilities::assert_eq
Definition: v3_utilities.f:62
mpi_inc
Umbrella module avoid multiple inlcudes of the mpif.h header.
Definition: mpi_inc.f:11
Generated on Thu Mar 5 2020 15:49:23 for V3FIT by   doxygen 1.8.17