4 USE vmec_input,
ONLY: nbfld, nflxs, lfreeb, lrecon
5 USE vsvd0,
ONLY: nigroup, nparts, npfcoil, nbcoilsp, nfloops,
9 INTEGER,
PARAMETER :: nlimset = 2
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',
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'
78 INTEGER :: nr0b, np0b, nfper0, nz0b
79 INTEGER :: nobd, nobser, nextcur, nbfldn, nbsets, nbcoilsn
80 INTEGER :: nbvac, nbcoil_max, nlim, nlim_max, nsets,
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
106 PRIVATE :: read_mgrid_bin, read_mgrid_nc
108 PRIVATE :: read_mgrid_bin
113 SUBROUTINE read_mgrid (mgrid_file, extcur, nv, nfp, lscreen,
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
138 CHARACTER(LEN=*),
PARAMETER :: mgrid_defarea=
'vmec$:[makegrid]'
140 CHARACTER(LEN=*),
PARAMETER :: mgrid_defarea=
'$HOME/vmec/MAKEGRID'
149 CHARACTER(LEN=200) :: home_dir
150 LOGICAL :: lgrid_exist, lfind
155 mpi_comm = mpi_comm_world
156 IF (
PRESENT(comm)) mpi_comm = comm
159 mgrid_path = trim(mgrid_file)
161 IF ((mgrid_path .eq. trim(mgrid_path_old)) .and.
162 1
ALLOCATED(curlabel))
THEN
163 print *,
' mgrid file previously parsed!'
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:)
179 home_dir = trim(home_dir) // mgrid_defarea(istat+5:)
182 home_dir = mgrid_defarea
184 mgrid_path = trim(home_dir) // mgrid_file(ii+1:)
185 INQUIRE (file=mgrid_path,exist=lgrid_exist,iostat=istat)
188 mgrid_path_old = mgrid_path
192 IF (lgrid_exist)
THEN
193 IF (lscreen) print
'(2x,2a)',
194 1
'Opening vacuum field file: ', trim(mgrid_file)
198 ii = len_trim(mgrid_path) - 2
199 lfind = (mgrid_path(ii:ii+2) ==
'.nc')
202 CALL read_mgrid_nc (mgrid_path, extcur, nv, nfp,
203 1 ier_flag, lscreen, comm)
205 lgrid_exist = .false.
208 CALL read_mgrid_bin (mgrid_path, extcur, nv, nfp,
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'
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
225 IF (ier_flag .ne. 0)
RETURN
227 IF (.not.lgrid_exist .or. ier_flag.ne.0)
THEN
231 print *,
' Error opening/reading mgrid file in 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'
240 END SUBROUTINE read_mgrid
243 SUBROUTINE read_mgrid_bin (filename, extcur, nv, nfp, ier_flag,
253 INTEGER,
INTENT(in) :: nv, nfp
254 CHARACTER(LEN=*),
INTENT(in) :: filename
255 REAL(rprec),
INTENT(in) :: extcur(:)
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
265 CALL safe_open(iunit, istat, filename,
'old',
'unformatted')
266 IF (istat .ne. 0)
THEN
271 READ (iunit,iostat=istat) nr0b, nz0b, np0b, nfper0, nextcur
272 IF (istat .ne. 0) ier_flag = 9
274 IF (nfper0.NE.nfp .OR. mod(np0b, nv).NE.0)
RETURN
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
281 IF (nextcur .eq. 0)
THEN
282 print *,
' NEXTCUR = 0 IN READING MGRID FILE'
284 ELSE IF (nextcur .gt. nigroup)
THEN
285 print *,
' NEXTCUR > NIGROUP IN MGRID FILE'
289 IF (ier_flag .ne. 0)
RETURN
291 ALLOCATE (curlabel(5*(nextcur/5+1)), stat=istat)
293 READ(iunit,iostat=istat) (curlabel(n),n=1,nextcur)
294 IF (istat .ne. 0)
THEN
295 print *,
' reading mgrid file failed (curlabel)'
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))
310 ALLOCATE (brtemp(nr0b,nz0b,np0b), bptemp(nr0b,nz0b,np0b),
311 1 bztemp(nr0b,nz0b,np0b), stat=istat)
313 IF (istat .ne. 0)
THEN
314 print *,
' allocation for b-vector storage failed'
322 IF (lstyle_2000)
THEN
323 READ(iunit, iostat=istat) brtemp, bptemp, bztemp
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)
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)
338 DEALLOCATE (brtemp, bztemp, bptemp)
341 CALL assign_bptrs(bvac)
343 IF (istat .ne. 0)
THEN
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'
363 READ(iunit,iostat=istat) nobser, nobd, nbsets
368 IF (lscreen) print *,
' No observation data in mgrid data'
372 nbfldn = sum(nbfld(:nbsets))
373 ALLOCATE (nbcoils(nbsets), stat=istat)
374 READ(iunit) (nbcoils(n),n=1,nbsets)
376 nbcoil_max = maxval(nbcoils(:nbsets))
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
389 1 print *,
' allocation error for xobser: istat = ',istat
394 IF (nobser .gt. nfloops)
THEN
395 print *,
'NOBSER>NFLOOPS'
398 IF (nobd .gt. nfloops)
THEN
399 print *,
'NOBD>NFLOOPS'
402 IF (nflxs .gt. nfloops)
THEN
403 print *,
'NFLXS>NFLOOPS'
406 IF (nbfldn .gt. nbctotp)
THEN
407 print *,
'NBFLDN>NBCTOTP'
410 IF (nbcoil_max .gt. nbcoilsp)
THEN
411 print *,
'NBCOIL_max>NBCOILSP'
415 IF (ier_flag .ne. 0)
RETURN
417 IF (nobser+nobd .gt. 0) iconnect(:4,:nobser+nobd) = 0
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)
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
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),
436 dbcoil(:nbcoil_max,:nbsets,:nextcur) = 0
440 READ(iunit) (unpsiext(n,ig),n=1,nobser)
442 READ(iunit) (dbcoil(m,n,ig),m=1,nbcoils(n))
449 ALLOCATE (limitr(nlimset), nsetsn(nigroup))
451 READ (iunit,iostat=istat) nlim,(limitr(i),i=1,nlim)
452 IF (istat .ne. 0)
then
454 IF (lscreen) print *,
' No limiter data in mgrid file'
458 nlim_max = maxval(limitr)
460 IF (nlim .gt. nlimset)
THEN
461 print *,
'nlim>nlimset'
466 ALLOCATE( rlim(nlim_max,nlim), zlim(nlim_max,nlim),
467 1 reslim(nlim_max,nlim) ,seplim(nlim_max,nlim),
469 IF (istat .ne. 0)
THEN
470 print *,
'rlim istat!=0'
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)
479 IF (nsets .gt. nigroup)
THEN
480 print *,
'nsets>nigroup'
483 ELSE IF (istat .ne. 0)
THEN
488 nsets_max = maxval(nsetsn)
490 IF (nsets_max .gt. npfcoil)
THEN
491 print *,
'nsetsn>npfcoil'
496 ALLOCATE (pfcspec(nparts,nsets_max,nsets), stat=istat)
500 READ(iunit, iostat=istat) (((pfcspec(i,j,k),i=1,nparts),
501 1 j=1,nsetsn(k)), k=1,nsets)
503 DEALLOCATE (limitr, nsetsn)
505 READ(iunit, iostat=istat) rx1,rx2,zy1,zy2,condif,
506 1 nrgrid,nzgrid,tokid
508 IF (istat .ne. 0)
THEN
513 IF (nobser .gt. 0) xobsqr(:nobser) = sqrt(xobser(:nobser))
517 nbcoilsn = sum(nbcoils(:nbsets))
520 rbcoilsqr(:nbcoils(n),n) = sqrt(rbcoil(:nbcoils(n),n))
527 delrb = (rmaxb-rminb)/(nr0b-1)
528 delzb = (zmaxb-zminb)/(nz0b-1)
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
540 1 psiext(:nobser) = psiext(:nobser) +
541 2 extcur(ig)*unpsiext(:nobser,ig)
544 bcoil(:n1,n) = bcoil(:n1,n) +
545 1 extcur(ig)*dbcoil(:n1,n,ig)
550 END SUBROUTINE read_mgrid_bin
557 SUBROUTINE read_mgrid_nc (filename, extcur, nv, nfp,
558 1 ier_flag, lscreen, comm)
565 CHARACTER(LEN=*),
INTENT(in) :: filename
566 INTEGER,
INTENT(in) :: nv, nfp
567 REAL(rprec),
INTENT(in) :: extcur(:)
568 INTEGER,
INTENT(in) :: comm
572 REAL(rprec),
DIMENSION(:,:,:),
ALLOCATABLE ::
573 1 brtemp, bztemp, bptemp
574 INTEGER :: ier_flag, ngrid
577 INTEGER,
DIMENSION(3) :: dimlens
578 CHARACTER(LEN=100) :: temp
580 INTEGER :: mpi_rank, mpi_size, lMPIInit, MPI_ERR
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)
587 mpi_rank = 0; mpi_size = 1
591 call cdf_open(ngrid, filename,
'r', istat)
592 IF (istat .ne. 0)
THEN
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)
605 IF (nfper0.NE.nfp .OR. mod(np0b, nv).NE.0)
RETURN
607 CALL cdf_read(ngrid, vn_nextcur, nextcur)
609 IF (nextcur .eq. 0)
THEN
610 print *,
' NEXTCUR = 0 IN READING MGRID FILE'
613 ELSE IF (nextcur .gt. nigroup)
THEN
614 print *,
' NEXTCUR > NIGROUP IN MGRID FILE'
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)
624 delrb = (rmaxb-rminb)/(nr0b-1)
625 delzb = (zmaxb-zminb)/(nz0b-1)
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)
635 IF (nextcur .eq. 1)
THEN
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))
642 IF (istat .ne. 0) stop
'Error allocating CURLABEL in mgrid_mod'
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)
652 IF (istat .ne. 0) stop
'Error allocating bvac in mgrid_mod'
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'
665 DO ig = mpi_rank + 1, nextcur, mpi_size
669 WRITE (temp, 1000) vn_br0, ig
670 CALL cdf_read(ngrid, temp, brtemp)
672 WRITE (temp, 1000) vn_bp0, ig
673 CALL cdf_read(ngrid, temp, bptemp)
675 WRITE (temp, 1000) vn_bz0, ig
676 CALL cdf_read(ngrid, temp, bztemp)
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)
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')
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)
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)
720 CALL cdf_close(ngrid)
722 IF (
ALLOCATED(brtemp))
723 1
DEALLOCATE (brtemp, bptemp, bztemp)
725 CALL assign_bptrs(bvac)
727 1000
FORMAT(a,
'_',i3.3)
729 END SUBROUTINE read_mgrid_nc
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)
740 bfield = bfield + cur*bf_add(:,1:np0b:nskip)
742 END SUBROUTINE sum_bfield
744 SUBROUTINE assign_bptrs(bptr)
746 REAL(rprec),
TARGET,
INTENT(in) :: bptr(nr0b,nz0b,np0b,3)
748 brvac => bptr(:,:,:,1)
749 bpvac => bptr(:,:,:,2)
750 bzvac => bptr(:,:,:,3)
752 END SUBROUTINE assign_bptrs
754 SUBROUTINE free_mgrid (istat)
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)
768 1
DEALLOCATE (rlim,zlim, reslim,seplim,stat=istat)
773 END SUBROUTINE free_mgrid