21 USE vmec_input,
ONLY: lrfp, lmove_axis, nbfld
30 CHARACTER(LEN=*),
PARAMETER :: vn_version =
'version_', &
31 vn_extension =
'input_extension', vn_mgrid =
'mgrid_file',
32 vn_magen =
'wb', vn_therm =
'wp', vn_gam =
'gamma',
33 vn_maxr =
'rmax_surf', vn_minr =
'rmin_surf',
34 vn_maxz =
'zmax_surf', vn_fp =
'nfp',
35 vn_radnod =
'ns', vn_polmod =
'mpol', vn_tormod =
'ntor',
36 vn_maxmod =
'mnmax', vn_maxit =
'niter', vn_actit =
'itfsq',
37 vn_maxpot =
'mnmaxpot', vn_potsin =
'potsin', vn_potcos =
'potcos'
38 vn_xmpot =
'xmpot', vn_xnpot=
'xnpot',
39 vn_asym =
'lasym', vn_recon =
'lrecon',
41 vn_free =
'lfreeb', vn_moveaxis =
'lmove_axis',
42 vn_error =
'ier_flag', vn_aspect =
'aspect', vn_rfp =
'lrfp',
43 vn_maxmod_nyq =
'mnmax_nyq',
44 vn_beta =
'betatotal', vn_pbeta =
'betapol',
45 vn_tbeta =
'betator', vn_abeta =
'betaxis',
46 vn_b0 =
'b0', vn_rbt0 =
'rbtor0', vn_rbt1 =
'rbtor',
47 vn_sgs =
'signgs', vn_lar =
'IonLarmor', vn_modb =
'volavgB',
48 vn_ctor =
'ctor', vn_amin =
'Aminor_p', vn_rmaj =
'Rmajor_p',
49 vn_vol =
'volume_p', vn_am =
'am', vn_ai =
'ai', vn_ac =
'ac',
50 vn_ah =
'hot particle fraction', vn_atuname =
'T-perp/T-par',
51 vn_pmass_type =
'pmass_type', vn_piota_type =
'piota_type',
52 vn_pcurr_type =
'pcurr_type',
53 vn_am_aux_s =
'am_aux_s', vn_am_aux_f =
'am_aux_f',
54 vn_ai_aux_s =
'ai_aux_s', vn_ai_aux_f =
'ai_aux_f',
55 vn_ac_aux_s =
'ac_aux_s', vn_ac_aux_f =
'ac_aux_f',
56 vn_mse =
'imse', vn_thom =
'itse',
57 vn_pmod =
'xm', vn_tmod =
'xn', vn_pmod_nyq =
'xm_nyq',
58 vn_tmod_nyq =
'xn_nyq',
59 vn_racc =
'raxis_cc', vn_zacs =
'zaxis_cs',
60 vn_racs =
'raxis_cs', vn_zacc =
'zaxis_cc', vn_iotaf =
'iotaf',
62 vn_qfact=
'q-factor', vn_chi=
'chi', vn_chipf=
'chipf',
63 vn_presf =
'presf', vn_phi =
'phi', vn_phipf =
'phipf',
64 vn_jcuru =
'jcuru', vn_jcurv =
'jcurv', vn_iotah =
'iotas',
65 vn_mass =
'mass', vn_presh =
'pres', vn_betah =
'beta_vol',
66 vn_buco =
'buco', vn_bvco =
'bvco', vn_vp =
'vp',
67 vn_specw =
'specw', vn_phip =
'phips', vn_jdotb =
'jdotb',
68 vn_bdotb =
'bdotb', vn_overr =
'over_r',
69 vn_bgrv =
'bdotgradv', vn_merc =
'DMerc', vn_mshear =
'DShear',
70 vn_mwell =
'DWell', vn_mcurr =
'DCurr', vn_mgeo =
'DGeod',
71 vn_equif =
'equif', vn_fsq =
'fsqt', vn_wdot =
'wdot',
72 vn_ftolv =
'ftolv', vn_fsql=
'fsql', vn_fsqr =
'fsqr',
74 vn_extcur =
'extcur', vn_curlab =
'curlabel', vn_rmnc =
'rmnc',
75 vn_zmns =
'zmns', vn_lmns =
'lmns', vn_gmnc =
'gmnc',
76 vn_bmnc =
'bmnc', vn_bsubumnc =
'bsubumnc',
77 vn_bsubvmnc =
'bsubvmnc', vn_bsubsmns =
'bsubsmns',
78 vn_bsupumnc =
'bsupumnc', vn_bsupvmnc =
'bsupvmnc',
79 vn_rmns =
'rmns', vn_zmnc =
'zmnc',
80 vn_lmnc =
'lmnc', vn_gmns =
'gmns', vn_bmns =
'bmns',
81 vn_bsubumns =
'bsubumns', vn_bsubvmns =
'bsubvmns',
82 vn_bsubsmnc =
'bsubsmnc', vn_bsupumns =
'bsupumns',
83 vn_currumnc =
'currumnc', vn_currumns =
'currumns',
84 vn_currvmnc =
'currvmnc', vn_currvmns =
'currvmns',
85 vn_bsupvmns =
'bsupvmns',
86 vn_bsubumnc_sur =
'bsubumnc_sur',
87 vn_bsubvmnc_sur =
'bsubvmnc_sur',
88 vn_bsupumnc_sur =
'bsupumnc_sur',
89 vn_bsupvmnc_sur =
'bsupvmnc_sur',
90 vn_bsubumns_sur =
'bsubumns_sur',
91 vn_bsubvmns_sur =
'bsubvmns_sur',
92 vn_bsupumns_sur =
'bsupumns_sur',
93 vn_bsupvmns_sur =
'bsupvmns_sur',
94 vn_rbc =
'rbc', vn_zbs =
'zbs', vn_rbs =
'rbs', vn_zbc =
'zbc'
96 CHARACTER(LEN=*),
PARAMETER :: ln_version =
'VMEC Version', &
97 ln_extension =
'Input file extension', ln_mgrid =
'MGRID file',
98 ln_magen =
'Magnetic Energy', ln_therm =
'Thermal Energy',
99 ln_gam =
'Gamma', ln_maxr =
'Maximum R', ln_minr =
'Minimum R',
100 ln_maxz =
'Maximum Z', ln_fp =
'Field Periods',
101 ln_radnod =
'Radial nodes', ln_polmod =
'Poloidal modes',
102 ln_tormod =
'Toroidal modes', ln_maxmod =
'Fourier modes',
103 ln_maxmod_nyq =
'Max # Fourier modes (Nyquist)',
104 ln_maxpot =
'Max # Fourier modes (vacuum potential)',
105 ln_potsin =
'Vacuum potential sin modes',
106 ln_potcos =
'Vacuum potential cos modes',
107 ln_xmpot =
'Vacuum potential poloidal modes',
108 ln_xnpot =
'Vacuum potential toroidal modes',
109 ln_maxit =
'Max iterations', ln_actit =
'Actual iterations',
110 ln_asym =
'Asymmetry', ln_recon =
'Reconstruction',
111 ln_free =
'Free boundary',
112 ln_error =
'Error flag', ln_aspect =
'Aspect ratio',
113 ln_beta =
'Total beta', ln_pbeta =
'Poloidal beta',
114 ln_tbeta =
'Toroidal beta', ln_abeta =
'Beta axis',
115 ln_b0 =
'RB-t over R axis', ln_rbt0 =
'RB-t axis',
116 ln_rbt1 =
'RB-t edge', ln_sgs =
'Sign jacobian',
117 ln_lar =
'Ion Larmor radius', ln_modb =
'avg mod B',
118 ln_ctor =
'Toroidal current', ln_amin =
'minor radius',
119 ln_rmaj =
'major radius', ln_vol =
'Plasma volume',
120 ln_mse =
'Number of MSE points',
121 ln_thom =
'Number of Thompson scattering points',
122 ln_am =
'Specification parameters for mass(s)',
123 ln_ac =
'Specification parameters for <J>(s)',
124 ln_ai =
'Specification parameters for iota(s)',
125 ln_pmass_type =
'Profile type specifier for mass(s)',
126 ln_pcurr_type =
'Profile type specifier for <J>(s)',
127 ln_piota_type =
'Profile type specifier for iota(s)',
128 ln_am_aux_s =
'Auxiliary-s parameters for mass(s)',
129 ln_am_aux_f =
'Auxiliary-f parameters for mass(s)',
130 ln_ac_aux_s =
'Auxiliary-s parameters for <J>(s)',
131 ln_ac_aux_f =
'Auxiliary-f parameters for <J>(s)',
132 ln_ai_aux_s =
'Auxiliary-s parameters for iota(s)',
133 ln_ai_aux_f =
'Auxiliary-f parameters for iota(s)',
134 ln_pmod =
'Poloidal mode numbers',
135 ln_tmod =
'Toroidal mode numbers',
136 ln_pmod_nyq =
'Poloidal mode numbers (Nyquist)',
137 ln_tmod_nyq =
'Toroidal mode numbers (Nyquist)',
138 ln_racc =
'raxis (cosnv)', ln_racs =
'raxis (sinnv)',
139 ln_zacs =
'zaxis (sinnv)', ln_zacc =
'zaxis (cosnv)',
140 ln_iotaf =
'iota on full mesh',
141 ln_qfact =
'q-factor on full mesh',
143 ln_presf =
'pressure on full mesh',
144 ln_phi =
'Toroidal flux on full mesh',
145 ln_phipf =
'd(phi)/ds: Toroidal flux deriv on full mesh',
146 ln_chi =
'Poloidal flux on full mesh',
148 ln_chipf =
'd(chi)/ds: Poroidal flux deriv on full mesh',
150 ln_jcuru =
'j dot gradu full',
151 ln_jcurv =
'j dot gradv full', ln_iotah =
'iota half',
152 ln_mass =
'mass half', ln_presh =
'pressure half',
153 ln_betah =
'beta half', ln_buco =
'bsubu half',
154 ln_bvco =
'bsubv half', ln_vp =
'volume deriv half',
155 ln_specw =
'Spectral width half',
156 ln_phip =
'tor flux deriv over 2pi half',
157 ln_jdotb =
'J dot B', ln_bgrv =
'B dot grad v',
158 ln_bdotb =
'B dot B',
159 ln_merc =
'Mercier criterion', ln_mshear =
'Shear Mercier',
160 ln_mwell =
'Well Mercier', ln_mcurr =
'Current Mercier',
161 ln_mgeo =
'Geodesic Mercier', ln_equif=
'Average force balance',
162 ln_fsq =
'Residual decay',
163 ln_wdot =
'Wdot decay', ln_extcur =
'External coil currents',
164 ln_fsqr =
'Residual decay - radial',
165 ln_fsqz =
'Residual decay - vertical',
166 ln_fsql =
'Residual decay - hoop',
167 ln_ftolv =
'Residual decay - requested',
168 ln_curlab =
'External current names',
170 ln_rmnc =
'cosmn component of cylindrical R, full mesh',
171 ln_zmns =
'sinmn component of cylindrical Z, full mesh',
172 ln_lmns =
'sinmn component of lambda, half mesh',
173 ln_gmnc =
'cosmn component of jacobian, half mesh',
174 ln_bmnc =
'cosmn component of mod-B, half mesh',
175 ln_bsubumnc =
'cosmn covariant u-component of B, half mesh',
176 ln_bsubvmnc =
'cosmn covariant v-component of B, half mesh',
177 ln_bsubsmns =
'sinmn covariant s-component of B, half mesh',
179 ln_bsubumnc_sur =
'cosmn bsubu of B, surface',
180 ln_bsubvmnc_sur =
'cosmn bsubv of B, surface',
181 ln_bsupumnc_sur =
'cosmn bsupu of B, surface',
182 ln_bsupvmnc_sur =
'cosmn bsupv of B, surface',
184 ln_bsupumnc =
'BSUPUmnc half', ln_bsupvmnc =
'BSUPVmnc half',
186 ln_rmns =
'sinmn component of cylindrical R, full mesh',
187 ln_zmnc =
'cosmn component of cylindrical Z, full mesh',
188 ln_lmnc =
'cosmn component of lambda, half mesh',
189 ln_gmns =
'sinmn component of jacobian, half mesh',
190 ln_bmns =
'sinmn component of mod-B, half mesh',
191 ln_bsubumns =
'sinmn covariant u-component of B, half mesh',
192 ln_bsubvmns =
'sinmn covariant v-component of B, half mesh',
193 ln_bsubsmnc =
'cosmn covariant s-component of B, half mesh',
195 ln_currumnc =
'cosmn covariant u-component of J, full mesh',
196 ln_currumns =
'sinmn covariant u-component of J, full mesh',
197 ln_currvmnc =
'cosmn covariant v-component of J, full mesh',
198 ln_currvmns =
'sinmn covariant v-component of J, full mesh',
200 ln_bsubumns_sur =
'sinmn bsubu of B, surface',
201 ln_bsubvmns_sur =
'sinmn bsubv of B, surface',
202 ln_bsupumns_sur =
'sinmn bsupu of B, surface',
203 ln_bsupvmns_sur =
'sinmn bsupv of B, surface',
205 ln_bsupumns =
'BSUPUmns half', ln_bsupvmns =
'BSUPVmns half',
206 ln_rbc =
'Initial boundary R cos(mu-nv) coefficients',
207 ln_zbs =
'Initial boundary Z sin(mu-nv) coefficients',
208 ln_rbs =
'Initial boundary R sin(mu-nv) coefficients',
209 ln_zbc =
'Initial boundary Z cos(mu-nv) coefficients'
214 INTEGER :: nfp, ns, mpol, ntor, mnmax, mnmax_nyq, itfsq, niter, &
215 iasym, ireconstruct, ierr_vmec, imse, itse, nstore_seq,
216 isnodes, ipnodes, imatch_phiedge, isigng, mnyq, nnyq, ntmax,
218 REAL(rprec) :: wb, wp, gamma, pfac, rmax_surf, rmin_surf, &
219 zmax_surf, aspect, betatot, betapol, betator, betaxis, b0,
220 tswgt, msewgt, flmwgt, bcwgt, phidiam, version_,
221 delphid, ionlarmor, volavgb,
222 fsql, fsqr, fsqz, ftolv,
223 aminor, rmajor, volume, rbtor, rbtor0, itor
224 REAL(rprec),
ALLOCATABLE :: rzl_local(:,:,:,:)
225 REAL(rprec),
DIMENSION(:,:),
ALLOCATABLE :: &
226 rmnc, zmns, lmns, rmns, zmnc, lmnc, bmnc, gmnc, bsubumnc,
227 bsubvmnc, bsubsmns, bsupumnc, bsupvmnc, currvmnc,
228 currumnc, bbc, raxis, zaxis
229 REAL(rprec),
DIMENSION(:,:),
ALLOCATABLE :: &
230 bmns, gmns, bsubumns, bsubvmns, bsubsmnc,
231 bsupumns, bsupvmns, currumns, currvmns
232 REAL(rprec),
DIMENSION(:),
ALLOCATABLE :: &
233 bsubumnc_sur, bsubumns_sur, bsubvmnc_sur, bsubvmns_sur,
234 bsupumnc_sur, bsupumns_sur, bsupvmnc_sur, bsupvmns_sur
235 REAL(rprec),
DIMENSION(:),
ALLOCATABLE :: &
236 iotas, iotaf, presf, phipf, mass, pres, beta_vol, xm, xn,
237 potsin, potcos, xmpot, xnpot, qfact, chipf, phi, chi,
238 xm_nyq, xn_nyq, phip, buco, bvco, vp, overr, jcuru, jcurv,
239 specw, jdotb, bdotb, bdotgradv, fsqt, wdot, am, ac, ai,
240 am_aux_s, am_aux_f, ac_aux_s, ac_aux_f, ai_aux_s, ai_aux_f,
241 dmerc, dshear, dwell, dcurr, dgeod, equif, extcur,
242 sknots, ystark, y2stark, pknots, ythom, y2thom,
243 anglemse, rmid, qmid, shear, presmid, alfa, curmid, rstark,
244 qmeas, datastark, rthom, datathom, dsiobt
245 LOGICAL :: lasym, lthreed, lwout_opened=.false.
246 CHARACTER :: mgrid_file*200, input_extension*100
247 CHARACTER :: pmass_type*20, pcurr_type*20, piota_type*20
249 INTEGER,
PARAMETER :: norm_term_flag=0, &
250 bad_jacobian_flag=1, more_iter_flag=2, jac75_flag=4
254 INTERFACE read_wout_file
255 MODULE PROCEDURE readw_and_open, readw_only
259 PRIVATE :: read_wout_text, read_wout_nc
261 PRIVATE :: read_wout_text
263 PRIVATE :: norm_term_flag, bad_jacobian_flag,
264 more_iter_flag, jac75_flag
268 SUBROUTINE readw_and_open(file_or_extension, ierr, iopen)
274 INTEGER,
INTENT(out) :: ierr
275 INTEGER,
OPTIONAL :: iopen
276 CHARACTER(LEN=*),
INTENT(in) :: file_or_extension
280 INTEGER,
PARAMETER :: iunit_init = 10
283 CHARACTER(len=LEN_TRIM(file_or_extension)+10) :: filename
293 CALL parse_extension(filename, file_or_extension, isnc)
298 CALL read_wout_nc(filename, ierr)
300 print *,
"NETCDF wout file can not be opened on this platform"
305 CALL safe_open (iunit, ierr, filename,
'old',
'formatted')
306 IF (ierr .eq. 0)
CALL read_wout_text(iunit, ierr)
310 IF (
PRESENT(iopen)) iopen = ierr
311 lwout_opened = (ierr .eq. 0)
316 IF (
ALLOCATED(xn))
THEN
317 lthreed = any(nint(xn) .ne. 0)
322 END SUBROUTINE readw_and_open
325 SUBROUTINE readw_only(iunit, ierr, iopen)
330 INTEGER,
INTENT(in) :: iunit
331 INTEGER,
INTENT(out):: ierr
332 INTEGER,
OPTIONAL :: iopen
337 CHARACTER(LEN=256) :: vmec_version
345 INQUIRE(unit=iunit, exist=exfile, name=vmec_version,iostat=istat)
346 IF (istat.ne.0 .or. .not.exfile)
THEN
347 print *,
' In READ_WOUT_FILE, Unit = ',iunit,
348 ' File = ',trim(vmec_version),
' DOES NOT EXIST'
349 IF (
PRESENT(iopen)) iopen = -1
353 IF (
PRESENT(iopen)) iopen = 0
356 CALL read_wout_text(iunit, ierr)
357 lwout_opened = (ierr .eq. 0)
358 lthreed = any(nint(xn) .ne. 0)
360 END SUBROUTINE readw_only
363 SUBROUTINE read_wout_text(iunit, ierr)
364 USE stel_constants,
ONLY: mu0
369 INTEGER :: iunit, ierr
373 REAL(rprec),
PARAMETER :: eps_w = 1.e-4_dp
377 INTEGER :: istat(15), i, j, k, js, m, n, n1, mn, nparts_in
378 CHARACTER(LEN=256) :: vmec_version
394 READ (iunit,
'(a)', iostat=istat(2), err=1000) vmec_version
396 i = index(vmec_version,
'=')
398 READ(vmec_version(i+1:len_trim(vmec_version)),*) version_
403 ierr_vmec = norm_term_flag
405 IF (version_ .le. (5.10 + eps_w))
THEN
406 READ (iunit, *, iostat=istat(2), err=1000) wb, wp, gamma,
407 pfac, nfp, ns, mpol, ntor, mnmax, itfsq, niter, iasym,
410 IF (version_ .lt. 6.54)
THEN
411 READ (iunit, *, iostat=istat(2), err=1000) wb, wp, gamma,
412 pfac, rmax_surf, rmin_surf
414 READ (iunit, *, iostat=istat(2), err=1000) wb, wp, gamma,
415 pfac, rmax_surf, rmin_surf, zmax_surf
417 IF (version_ .le. (8.0+eps_w))
THEN
418 READ (iunit, *, iostat=istat(2), err=1000) nfp, ns, mpol,
419 ntor, mnmax, itfsq, niter, iasym, ireconstruct, ierr_vmec
422 READ (iunit, *, iostat=istat(2), err=1000) nfp, ns, mpol,
423 ntor, mnmax, mnmax_nyq, itfsq, niter, iasym, ireconstruct,
428 lasym = (iasym .gt. 0)
430 IF (version_ .gt. (6.20+eps_w))
THEN
431 READ (iunit, *, iostat=istat(1), err=1000)imse, itse, nbsets,
432 nobd, nextcur, nstore_seq
434 READ (iunit, *, iostat=istat(1), err=1000)imse, itse, nbsets,
439 IF (ierr_vmec.ne.norm_term_flag .and.
440 ierr_vmec.ne.more_iter_flag)
GOTO 1000
442 IF (nextcur .gt. nigroup) istat(15) = -1
444 IF (
ALLOCATED(xm))
CALL read_wout_deallocate
446 ALLOCATE (xm(mnmax), xn(mnmax), xm_nyq(mnmax_nyq),
447 xn_nyq(mnmax_nyq),rmnc(mnmax,ns), zmns(mnmax,ns),
448 lmns(mnmax,ns), bmnc(mnmax_nyq,ns), gmnc(mnmax_nyq,ns),
449 bsubumnc(mnmax_nyq,ns), bsubvmnc(mnmax_nyq,ns),
450 bsubsmns(mnmax_nyq,ns), bsupumnc(mnmax_nyq,ns),
451 bsupvmnc(mnmax_nyq,ns), currvmnc(mnmax_nyq,ns),
452 currumnc(mnmax_nyq,ns),
453 iotas(ns), mass(ns), pres(ns), beta_vol(ns), phip(ns),
454 buco(ns), bvco(ns), phi(ns), iotaf(ns), presf(ns), phipf(ns),
456 vp(ns), overr(ns), jcuru(ns), jcurv(ns), specw(ns), dmerc(ns),
457 dshear(ns), dwell(ns), dcurr(ns), dgeod(ns), equif(ns),
458 raxis(0:ntor,2), zaxis(0:ntor,2), jdotb(ns), bdotb(ns),
459 bdotgradv(ns), am(0:20), ac(0:20), ai(0:20),
460 fsqt(nstore_seq), wdot(nstore_seq), stat = istat(6))
462 IF (nextcur .GT. 0)
ALLOCATE(extcur(nextcur), curlabel(nextcur),
466 ALLOCATE (rmns(mnmax,ns), zmnc(mnmax,ns), lmnc(mnmax,ns),
467 bmns(mnmax_nyq,ns), gmns(mnmax_nyq,ns),
468 bsubumns(mnmax_nyq,ns),
469 bsubvmns(mnmax_nyq,ns), bsubsmnc(mnmax_nyq,ns),
470 bsupumns(mnmax_nyq,ns), bsupvmns(mnmax_nyq,ns),
471 currumns(mnmax_nyq,ns), currvmns(mnmax_nyq,ns),
474 fsqt = 0; wdot = 0; raxis = 0; zaxis = 0
476 IF (nbsets .gt. 0)
READ (iunit, *, iostat=istat(4), err=1000)
477 (nbfld(i),i=1,nbsets)
478 READ (iunit,
'(a)', iostat=istat(5), err=1000) mgrid_file
483 READ (iunit, *, iostat=istat(7), err=1000) m, n
484 xm(mn) = real(m,rprec)
485 xn(mn) = real(n,rprec)
487 IF (version_ .le. (6.20+eps_w))
THEN
488 READ (iunit, 730, iostat=istat(8), err=1000)
489 rmnc(mn,js), zmns(mn,js), lmns(mn,js), bmnc(mn,js),
490 gmnc(mn,js), bsubumnc(mn,js), bsubvmnc(mn,js),
491 bsubsmns(mn,js), bsupumnc(mn,js), bsupvmnc(mn,js),
493 ELSE IF (version_ .le. (8.0+eps_w))
THEN
494 READ (iunit, *, iostat=istat(8), err=1000)
495 rmnc(mn,js), zmns(mn,js), lmns(mn,js), bmnc(mn,js),
496 gmnc(mn,js), bsubumnc(mn,js), bsubvmnc(mn,js),
497 bsubsmns(mn,js), bsupumnc(mn,js), bsupvmnc(mn,js),
500 READ (iunit, *, iostat=istat(8), err=1000)
501 rmnc(mn,js), zmns(mn,js), lmns(mn,js)
504 IF (version_ .le. (8.0+eps_w))
THEN
505 READ (iunit, *, iostat=istat(8), err=1000)
506 rmns(mn,js), zmnc(mn,js), lmnc(mn,js),
507 bmns(mn,js), gmns(mn,js), bsubumns(mn,js),
508 bsubvmns(mn,js), bsubsmnc(mn,js),
509 bsupumns(mn,js), bsupvmns(mn,js)
511 READ (iunit, *, iostat=istat(8), err=1000)
512 rmns(mn,js), zmnc(mn,js), lmnc(mn,js)
515 IF (js.eq.1 .and. m.eq.0)
THEN
517 IF (n1 .le. ntor)
THEN
518 raxis(n1,1) = rmnc(mn,1)
519 zaxis(n1,1) = zmns(mn,1)
521 raxis(n1,2) = rmns(mn,1)
522 zaxis(n1,2) = zmnc(mn,1)
528 IF (version_ .le. (8.0+eps_w)) cycle
531 READ (iunit, *, iostat=istat(7), err=1000) m, n
532 xm_nyq(mn) = real(m,rprec)
533 xn_nyq(mn) = real(n,rprec)
535 READ (iunit, *, iostat=istat(8), err=1000)
536 bmnc(mn,js), gmnc(mn,js), bsubumnc(mn,js),
537 bsubvmnc(mn,js), bsubsmns(mn,js),
538 bsupumnc(mn,js), bsupvmnc(mn,js)
540 READ (iunit, *, iostat=istat(8), err=1000)
541 bmns(mn,js), gmns(mn,js), bsubumns(mn,js),
542 bsubvmns(mn,js), bsubsmnc(mn,js),
543 bsupumns(mn,js), bsupvmns(mn,js)
550 IF (version_ .gt. (8.0+eps_w))
THEN
551 CALL compute_currents(bsubsmnc, bsubsmns,
552 & bsubumnc, bsubumns,
553 & bsubvmnc, bsubvmns,
554 & xm_nyq, xn_nyq, mnmax_nyq, lasym, ns,
555 & currumnc, currvmnc,
556 & currumns, currvmns)
559 mnyq = int(maxval(xm_nyq)); nnyq = int(maxval(abs(xn_nyq)))/nfp
570 iotas(1) = 0; mass(1) = 0; pres(1) = 0; phip(1) = 0;
571 buco(1) = 0; bvco(1) = 0; vp(1) = 0; overr(1) = 0; specw(1) = 1
574 IF (version_ .le. (6.05+eps_w))
THEN
575 READ (iunit, 730, iostat=istat(9), err=1000)
576 (iotas(js), mass(js), pres(js),
577 phip(js), buco(js), bvco(js), phi(js), vp(js), overr(js),
578 jcuru(js), jcurv(js), specw(js),js=2,ns)
579 READ (iunit, 730, iostat=istat(10), err=1000)
580 aspect, betatot, betapol, betator, betaxis, b0
581 ELSE IF (version_ .le. (6.20+eps_w))
THEN
582 READ (iunit, 730, iostat=istat(9), err=1000)
583 (iotas(js), mass(js), pres(js), beta_vol(js),
584 phip(js), buco(js), bvco(js), phi(js), vp(js), overr(js),
585 jcuru(js), jcurv(js), specw(js),js=2,ns)
586 READ (iunit, 730, iostat=istat(10), err=1000)
587 aspect, betatot, betapol, betator, betaxis, b0
588 ELSE IF (version_ .le. (6.95+eps_w))
THEN
589 READ (iunit, *, iostat=istat(9), err=1000)
590 (iotas(js), mass(js), pres(js), beta_vol(js),
591 phip(js), buco(js), bvco(js), phi(js), vp(js), overr(js),
592 jcuru(js), jcurv(js), specw(js),js=2,ns)
593 READ (iunit, *, iostat=istat(10), err=1000)
594 aspect, betatot, betapol, betator, betaxis, b0
596 READ (iunit, *, iostat=istat(9), err=1000)
597 (iotaf(js), presf(js), phipf(js), phi(js),
598 jcuru(js), jcurv(js), js=1,ns)
599 READ (iunit, *, iostat=istat(9), err=1000)
600 (iotas(js), mass(js), pres(js),
601 beta_vol(js), phip(js), buco(js), bvco(js), vp(js),
602 overr(js), specw(js),js=2,ns)
603 READ (iunit, *, iostat=istat(10), err=1000)
604 aspect, betatot, betapol, betator, betaxis, b0
608 IF (version_ .gt. (6.10+eps_w))
THEN
609 READ (iunit, *, iostat=istat(10), err=1000) isigng
610 READ (iunit, *, iostat=istat(10), err=1000) input_extension
611 READ (iunit, *, iostat=istat(10), err=1000) ionlarmor,
612 volavgb, rbtor0, rbtor, itor, aminor, rmajor, volume
618 IF (version_.gt.(5.10+eps_w) .and. version_.lt.(6.20-eps_w))
THEN
619 READ (iunit, 730, iostat=istat(11), err=1000)
620 (dmerc(js), dshear(js), dwell(js), dcurr(js),
621 dgeod(js), equif(js), js=2,ns-1)
622 ELSE IF (version_ .ge. (6.20-eps_w))
THEN
623 READ (iunit, *, iostat=istat(11), err=1000)
624 (dmerc(js), dshear(js), dwell(js), dcurr(js),
625 dgeod(js), equif(js), js=2,ns-1)
628 IF (nextcur .gt. 0)
THEN
629 IF (version_ .le. (6.20+eps_w))
THEN
630 READ (iunit, 730, iostat=istat(12), err=1000)
631 (extcur(i),i=1,nextcur)
633 READ (iunit, *, iostat=istat(12), err=1000)
634 (extcur(i),i=1,nextcur)
636 READ (iunit, *, iostat=istat(13)) lcurr
637 IF (lcurr)
READ (iunit, *, iostat=istat(13), err=1000)
638 (curlabel(i),i=1,nextcur)
641 IF (version_ .le. (6.20+eps_w))
THEN
642 READ (iunit, 730, iostat=istat(14))
643 (fsqt(i), wdot(i), i=1,nstore_seq)
645 READ (iunit, *, iostat=istat(14))
646 (fsqt(i), wdot(i), i=1,nstore_seq)
649 IF ((version_.ge.6.20-eps_w) .and. (version_ .lt. (6.50-eps_w))
650 .and. (istat(14).eq.0))
THEN
651 READ (iunit, 730, iostat=istat(14), err=1000)
652 (jdotb(js), bdotgradv(js), bdotb(js), js=1,ns)
653 ELSE IF (version_ .ge. (6.50-eps_w))
THEN
654 READ (iunit, *, iostat=istat(14), err=1000)
655 (jdotb(js), bdotgradv(js), bdotb(js), js=1,ns)
664 IF (version_ .le. (6.05+eps_w))
THEN
676 IF (ireconstruct .gt. 0)
THEN
678 n1 = maxval(nbfld(:nbsets))
679 ALLOCATE (sknots(isnodes), ystark(isnodes), y2stark(isnodes),
680 pknots(ipnodes), ythom(ipnodes), y2thom(ipnodes),
681 anglemse(2*ns), rmid(2*ns), qmid(2*ns), shear(2*ns),
682 presmid(2*ns), alfa(2*ns), curmid(2*ns), rstark(imse),
683 datastark(imse), rthom(itse), datathom(itse),
684 dsiext(nobd), plflux(nobd), dsiobt(nobd), bcoil(n1,nbsets),
685 plbfld(n1,nbsets), bbc(n1,nbsets))
686 IF (imse.ge.2 .or. itse.gt.0)
THEN
687 READ (iunit, *) tswgt, msewgt
688 READ (iunit, *) isnodes, (sknots(i),ystark(i),y2stark(i),
690 READ (iunit, *) ipnodes, (pknots(i), ythom(i),
691 y2thom(i),i=1,ipnodes)
692 READ(iunit, *)(anglemse(i),rmid(i),qmid(i),shear(i),
693 presmid(i),alfa(i),curmid(i),i=1,2*ns-1)
694 READ(iunit, *)(rstark(i),datastark(i),qmeas(i),i=1,imse)
695 READ(iunit, *)(rthom(i),datathom(i),i=1,itse)
698 IF (nobd .gt. 0)
THEN
699 READ (iunit, *) (dsiext(i),plflux(i),dsiobt(i),i=1,nobd)
700 READ (iunit, *) flmwgt
703 nbfldn = sum(nbfld(:nbsets))
704 IF (nbfldn .gt. 0)
THEN
706 READ (iunit, *) (bcoil(i,n),plbfld(i,n),bbc(i,n),
709 READ (iunit, *) bcwgt
712 READ (iunit, *) phidiam, delphid
716 READ (iunit, *) nsets, nparts_in, nlim
718 ALLOCATE (nsetsn(nsets))
719 READ (iunit, *) (nsetsn(i),i=1,nsets)
721 n1 = maxval(nsetsn(:nsets))
722 ALLOCATE (pfcspec(nparts_in,n1,nsets), limitr(nlim))
724 READ (iunit, *) (((pfcspec(i,j,k),i=1,nparts_in),
725 j=1,nsetsn(k)),k=1,nsets)
727 READ (iunit, *) (limitr(i), i=1,nlim)
729 m = maxval(limitr(:nlim))
730 ALLOCATE (rlim(m,nlim), zlim(m,nlim))
732 READ (iunit, *) ((rlim(i,j),zlim(i,j),i=1,limitr(j)),
734 READ (iunit, *) nrgrid, nzgrid
735 READ (iunit, *) tokid
736 READ (iunit, *) rx1, rx2, zy1, zy2, condif
737 READ (iunit, *) imatch_phiedge
743 READ (iunit, iostat=ierr) mgrid_mode
744 IF (ierr .ne. 0)
THEN
745 ierr = 0; mgrid_mode =
'N'
748 IF (istat(2) .ne. 0) ierr_vmec = 1
751 IF (istat(m) .gt. 0)
THEN
752 print *,
' Error No. ',m,
' in READ_WOUT, iostat = ',istat(m)
762 790
FORMAT(i5,/,(1p,3e12.4))
764 END SUBROUTINE read_wout_text
768 SUBROUTINE read_wout_nc(filename, ierr)
770 USE stel_constants,
ONLY: mu0
775 INTEGER,
INTENT(out) :: ierr
776 CHARACTER(LEN=*),
INTENT(in) :: filename
780 INTEGER :: nwout, ierror
781 INTEGER,
DIMENSION(3) :: dimlens
783 REAL(rprec),
DIMENSION(:),
ALLOCATABLE :: raxis_cc, raxis_cs, &
787 CALL cdf_open(nwout,filename,
'r', ierr)
788 IF (ierr .ne. 0)
THEN
789 print *,
' Error opening wout .nc file'
794 CALL read_wout_deallocate
797 CALL cdf_read(nwout, vn_error, ierr_vmec)
799 IF (ierr_vmec.ne.norm_term_flag .and.
800 ierr_vmec.ne.more_iter_flag)
GOTO 1000
802 CALL cdf_read(nwout, vn_version, version_)
803 CALL cdf_read(nwout, vn_extension, input_extension)
804 CALL cdf_read(nwout, vn_mgrid, mgrid_file)
805 CALL cdf_read(nwout, vn_magen, wb)
806 CALL cdf_read(nwout, vn_therm, wp)
807 CALL cdf_read(nwout, vn_gam, gamma)
808 CALL cdf_read(nwout, vn_maxr, rmax_surf)
809 CALL cdf_read(nwout, vn_minr, rmin_surf)
810 CALL cdf_read(nwout, vn_maxz, zmax_surf)
811 CALL cdf_read(nwout, vn_fp, nfp)
812 CALL cdf_read(nwout, vn_radnod, ns)
813 CALL cdf_read(nwout, vn_polmod, mpol)
814 CALL cdf_read(nwout, vn_tormod, ntor)
815 CALL cdf_read(nwout, vn_maxmod, mnmax)
817 CALL cdf_read(nwout, vn_maxpot, mnmaxpot)
819 CALL cdf_read(nwout, vn_maxmod_nyq, mnmax_nyq)
820 CALL cdf_read(nwout, vn_maxit, niter)
821 CALL cdf_read(nwout, vn_actit, itfsq)
822 CALL cdf_read(nwout, vn_asym, lasym)
824 CALL cdf_read(nwout, vn_recon, lrecon)
825 IF (lrecon) ireconstruct = 1
826 CALL cdf_read(nwout, vn_free, lfreeb)
827 CALL cdf_read(nwout, vn_moveaxis, lmove_axis)
829 CALL cdf_read(nwout, vn_rfp, lrfp)
830 CALL cdf_read(nwout, vn_aspect, aspect)
831 CALL cdf_read(nwout, vn_beta, betatot)
832 CALL cdf_read(nwout, vn_pbeta, betapol)
833 CALL cdf_read(nwout, vn_tbeta, betator)
834 CALL cdf_read(nwout, vn_abeta, betaxis)
835 CALL cdf_read(nwout, vn_b0, b0)
836 CALL cdf_read(nwout, vn_rbt0, rbtor0)
837 CALL cdf_read(nwout, vn_rbt1, rbtor)
838 CALL cdf_read(nwout, vn_sgs, isigng)
839 CALL cdf_read(nwout, vn_lar, ionlarmor)
840 CALL cdf_read(nwout, vn_modb, volavgb)
841 CALL cdf_read(nwout, vn_ctor, itor)
842 CALL cdf_read(nwout, vn_amin, aminor)
843 CALL cdf_read(nwout, vn_rmaj, rmajor)
844 CALL cdf_read(nwout, vn_vol, volume)
845 CALL cdf_read(nwout, vn_ftolv, ftolv)
846 CALL cdf_read(nwout, vn_fsqr, fsqr)
847 CALL cdf_read(nwout, vn_fsqz, fsqz)
848 CALL cdf_read(nwout, vn_fsql, fsql)
849 CALL cdf_read(nwout, vn_pcurr_type, pcurr_type)
850 CALL cdf_read(nwout, vn_piota_type, piota_type)
851 CALL cdf_read(nwout, vn_pmass_type, pmass_type)
854 CALL cdf_read(nwout, vn_mse, imse)
855 CALL cdf_read(nwout, vn_thom, itse)
857 CALL cdf_read(nwout, vn_nextcur, nextcur)
860 CALL cdf_inquire(nwout, vn_mgmode, dimlens, ier=ierror)
861 IF (ierror.eq.0)
CALL cdf_read(nwout, vn_mgmode, mgrid_mode)
863 CALL cdf_read(nwout, vn_flp, nobser)
864 CALL cdf_read(nwout, vn_nobd, nobd)
865 CALL cdf_read(nwout, vn_nbset, nbsets)
870 IF (lfreeb .and. nbsets.gt.0)
THEN
871 CALL cdf_read(nwout, vn_nbfld, nbfld)
874 CALL cdf_inquire(nwout, vn_pmod, dimlens)
875 ALLOCATE (xm(dimlens(1)), stat = ierror)
876 CALL cdf_inquire(nwout, vn_tmod, dimlens)
877 ALLOCATE (xn(dimlens(1)), stat = ierror)
878 IF (mnmax_nyq .gt. 0)
THEN
879 CALL cdf_inquire(nwout, vn_pmod_nyq, dimlens)
880 ALLOCATE (xm_nyq(dimlens(1)), stat = ierror)
881 CALL cdf_inquire(nwout, vn_tmod_nyq, dimlens)
882 ALLOCATE (xn_nyq(dimlens(1)), stat = ierror)
885 IF (mnmaxpot .gt. 0)
THEN
886 CALL cdf_inquire(nwout, vn_potsin, dimlens)
887 ALLOCATE (potsin(dimlens(1)), xmpot(dimlens(1)),
888 xnpot(dimlens(1)), stat = ierror)
889 IF (lasym)
ALLOCATE (potcos(dimlens(1)), stat = ierror)
892 CALL cdf_inquire(nwout, vn_racc, dimlens)
893 ALLOCATE (raxis_cc(0:dimlens(1)-1), stat = ierror)
894 CALL cdf_inquire(nwout, vn_zacs, dimlens)
895 ALLOCATE (zaxis_cs(0:dimlens(1)-1), stat = ierror)
897 CALL cdf_inquire(nwout, vn_racs, dimlens)
898 ALLOCATE (raxis_cs(0:dimlens(1)-1), stat = ierror)
899 CALL cdf_inquire(nwout, vn_zacc, dimlens)
900 ALLOCATE (zaxis_cc(0:dimlens(1)-1), stat = ierror)
904 CALL cdf_inquire(nwout, vn_am, dimlens)
905 ALLOCATE (am(0:dimlens(1)-1), stat = ierror)
906 CALL cdf_inquire(nwout, vn_ac, dimlens)
907 ALLOCATE (ac(0:dimlens(1)-1), stat = ierror)
908 CALL cdf_inquire(nwout, vn_ai, dimlens)
909 ALLOCATE (ai(0:dimlens(1)-1), stat = ierror)
911 CALL cdf_inquire(nwout, vn_ac_aux_s, dimlens)
912 ALLOCATE (ac_aux_s(dimlens(1)), stat = ierror)
913 CALL cdf_inquire(nwout, vn_ac_aux_f, dimlens)
914 ALLOCATE (ac_aux_f(dimlens(1)), stat = ierror)
915 CALL cdf_inquire(nwout, vn_ai_aux_s, dimlens)
916 ALLOCATE (ai_aux_s(dimlens(1)), stat = ierror)
917 CALL cdf_inquire(nwout, vn_ai_aux_f, dimlens)
918 ALLOCATE (ai_aux_f(dimlens(1)), stat = ierror)
919 CALL cdf_inquire(nwout, vn_am_aux_s, dimlens)
920 ALLOCATE (am_aux_s(dimlens(1)), stat = ierror)
921 CALL cdf_inquire(nwout, vn_am_aux_f, dimlens)
922 ALLOCATE (am_aux_f(dimlens(1)), stat = ierror)
924 CALL cdf_inquire(nwout, vn_iotaf, dimlens)
925 ALLOCATE (iotaf(dimlens(1)), stat = ierror)
926 CALL cdf_inquire(nwout, vn_qfact, dimlens)
928 ALLOCATE (qfact(dimlens(1)), stat = ierror)
930 CALL cdf_inquire(nwout, vn_presf, dimlens)
931 ALLOCATE (presf(dimlens(1)), stat = ierror)
932 CALL cdf_inquire(nwout, vn_phi, dimlens)
933 ALLOCATE (phi(dimlens(1)), stat = ierror)
934 CALL cdf_inquire(nwout, vn_chi, dimlens)
936 ALLOCATE (chi(dimlens(1)), stat = ierror)
938 CALL cdf_inquire(nwout, vn_phipf, dimlens)
939 ALLOCATE (phipf(dimlens(1)), stat = ierror)
941 ALLOCATE (chipf(dimlens(1)), stat = ierror)
944 CALL cdf_inquire(nwout, vn_jcuru, dimlens)
945 ALLOCATE (jcuru(dimlens(1)), stat = ierror)
946 CALL cdf_inquire(nwout, vn_jcurv, dimlens)
947 ALLOCATE (jcurv(dimlens(1)), stat = ierror)
948 CALL cdf_inquire(nwout, vn_iotah, dimlens)
949 ALLOCATE (iotas(dimlens(1)), stat = ierror)
950 CALL cdf_inquire(nwout, vn_mass, dimlens)
951 ALLOCATE (mass(dimlens(1)), stat = ierror)
952 CALL cdf_inquire(nwout, vn_presh, dimlens)
953 ALLOCATE (pres(dimlens(1)), stat = ierror)
954 CALL cdf_inquire(nwout, vn_betah, dimlens)
955 ALLOCATE (beta_vol(dimlens(1)), stat = ierror)
956 CALL cdf_inquire(nwout, vn_buco, dimlens)
957 ALLOCATE (buco(dimlens(1)), stat = ierror)
958 CALL cdf_inquire(nwout, vn_bvco, dimlens)
959 ALLOCATE (bvco(dimlens(1)), stat = ierror)
960 CALL cdf_inquire(nwout, vn_vp, dimlens)
961 ALLOCATE (vp(dimlens(1)), stat = ierror)
962 CALL cdf_inquire(nwout, vn_specw, dimlens)
963 ALLOCATE (specw(dimlens(1)), stat = ierror)
964 CALL cdf_inquire(nwout, vn_phip, dimlens)
965 ALLOCATE (phip(dimlens(1)), stat = ierror)
966 CALL cdf_inquire(nwout, vn_overr, dimlens)
967 ALLOCATE (overr(dimlens(1)), stat = ierror)
969 CALL cdf_inquire(nwout, vn_jdotb, dimlens)
970 ALLOCATE (jdotb(dimlens(1)), stat = ierror)
971 CALL cdf_inquire(nwout, vn_bdotb, dimlens)
972 ALLOCATE (bdotb(dimlens(1)), stat = ierror)
973 CALL cdf_inquire(nwout, vn_bgrv, dimlens)
974 ALLOCATE (bdotgradv(dimlens(1)), stat = ierror)
976 CALL cdf_inquire(nwout, vn_merc, dimlens)
977 ALLOCATE (dmerc(dimlens(1)), stat = ierror)
978 CALL cdf_inquire(nwout, vn_mshear, dimlens)
979 ALLOCATE (dshear(dimlens(1)), stat = ierror)
980 CALL cdf_inquire(nwout, vn_mwell, dimlens)
981 ALLOCATE (dwell(dimlens(1)), stat = ierror)
982 CALL cdf_inquire(nwout, vn_mcurr, dimlens)
983 ALLOCATE (dcurr(dimlens(1)), stat = ierror)
984 CALL cdf_inquire(nwout, vn_mgeo, dimlens)
985 ALLOCATE (dgeod(dimlens(1)), stat = ierror)
986 CALL cdf_inquire(nwout, vn_equif, dimlens)
987 ALLOCATE (equif(dimlens(1)), stat = ierror)
989 CALL cdf_inquire(nwout, vn_fsq, dimlens)
990 ALLOCATE (fsqt(dimlens(1)), stat = ierror)
991 CALL cdf_inquire(nwout, vn_wdot, dimlens)
992 ALLOCATE (wdot(dimlens(1)), stat = ierror)
994 IF (nextcur .gt. 0)
THEN
995 CALL cdf_inquire(nwout, vn_extcur, dimlens)
996 ALLOCATE (extcur(dimlens(1)), stat = ierror)
999 CALL cdf_inquire(nwout, vn_curlab, dimlens)
1000 ALLOCATE (curlabel(dimlens(2)), stat = ierror)
1004 CALL cdf_inquire(nwout, vn_rmnc, dimlens)
1005 ALLOCATE (rmnc(dimlens(1),dimlens(2)), stat = ierror)
1006 CALL cdf_inquire(nwout, vn_zmns, dimlens)
1007 ALLOCATE (zmns(dimlens(1),dimlens(2)), stat = ierror)
1008 CALL cdf_inquire(nwout, vn_lmns, dimlens)
1009 ALLOCATE (lmns(dimlens(1),dimlens(2)), stat = ierror)
1010 CALL cdf_inquire(nwout, vn_gmnc, dimlens)
1011 ALLOCATE (gmnc(dimlens(1),dimlens(2)), stat = ierror)
1012 CALL cdf_inquire(nwout, vn_bmnc, dimlens)
1013 ALLOCATE (bmnc(dimlens(1),dimlens(2)), stat = ierror)
1014 CALL cdf_inquire(nwout, vn_bsubumnc, dimlens)
1015 ALLOCATE (bsubumnc(dimlens(1),dimlens(2)), stat = ierror)
1016 CALL cdf_inquire(nwout, vn_bsubvmnc, dimlens)
1017 ALLOCATE (bsubvmnc(dimlens(1),dimlens(2)), stat = ierror)
1018 CALL cdf_inquire(nwout, vn_bsubsmns, dimlens)
1019 ALLOCATE (bsubsmns(dimlens(1),dimlens(2)), stat = ierror)
1022 CALL cdf_inquire(nwout, vn_bsupumnc, dimlens)
1023 ALLOCATE (bsupumnc(dimlens(1),dimlens(2)), stat = ierror)
1024 CALL cdf_inquire(nwout, vn_bsupvmnc, dimlens)
1025 ALLOCATE (bsupvmnc(dimlens(1),dimlens(2)), stat = ierror)
1029 ALLOCATE (currumnc(dimlens(1),dimlens(2)), stat = ierror)
1030 ALLOCATE (currvmnc(dimlens(1),dimlens(2)), stat = ierror)
1033 CALL cdf_inquire(nwout, vn_bsubumnc_sur, dimlens)
1034 ALLOCATE (bsubumnc_sur(dimlens(1)), stat = ierror)
1035 CALL cdf_inquire(nwout, vn_bsubvmnc_sur, dimlens)
1036 ALLOCATE (bsubvmnc_sur(dimlens(1)), stat = ierror)
1037 CALL cdf_inquire(nwout, vn_bsupumnc_sur, dimlens)
1038 ALLOCATE (bsupumnc_sur(dimlens(1)), stat = ierror)
1039 CALL cdf_inquire(nwout, vn_bsupvmnc_sur, dimlens)
1040 ALLOCATE (bsupvmnc_sur(dimlens(1)), stat = ierror)
1043 IF (.NOT. lasym)
GO TO 800
1045 CALL cdf_inquire(nwout, vn_rmns, dimlens)
1046 ALLOCATE (rmns(dimlens(1),dimlens(2)), stat = ierror)
1047 CALL cdf_inquire(nwout, vn_zmnc, dimlens)
1048 ALLOCATE (zmnc(dimlens(1),dimlens(2)), stat = ierror)
1049 CALL cdf_inquire(nwout, vn_lmnc, dimlens)
1050 ALLOCATE (lmnc(dimlens(1),dimlens(2)), stat = ierror)
1051 CALL cdf_inquire(nwout, vn_gmns, dimlens)
1052 ALLOCATE (gmns(dimlens(1),dimlens(2)), stat = ierror)
1053 CALL cdf_inquire(nwout, vn_bmns, dimlens)
1054 ALLOCATE (bmns(dimlens(1),dimlens(2)), stat = ierror)
1055 CALL cdf_inquire(nwout, vn_bsubumns, dimlens)
1056 ALLOCATE (bsubumns(dimlens(1),dimlens(2)), stat = ierror)
1057 CALL cdf_inquire(nwout, vn_bsubvmns, dimlens)
1058 ALLOCATE (bsubvmns(dimlens(1),dimlens(2)), stat = ierror)
1059 CALL cdf_inquire(nwout, vn_bsubsmnc, dimlens)
1060 ALLOCATE (bsubsmnc(dimlens(1),dimlens(2)), stat = ierror)
1063 CALL cdf_inquire(nwout, vn_bsupumns, dimlens)
1064 ALLOCATE (bsupumns(dimlens(1),dimlens(2)), stat = ierror)
1065 CALL cdf_inquire(nwout, vn_bsupvmns, dimlens)
1066 ALLOCATE (bsupvmns(dimlens(1),dimlens(2)), stat = ierror)
1070 ALLOCATE (currumns(dimlens(1),dimlens(2)), stat = ierror)
1071 ALLOCATE (currvmns(dimlens(1),dimlens(2)), stat = ierror)
1074 CALL cdf_inquire(nwout, vn_bsubumns_sur, dimlens)
1075 ALLOCATE (bsubumns_sur(dimlens(1)), stat = ierror)
1076 CALL cdf_inquire(nwout, vn_bsubvmns_sur, dimlens)
1077 ALLOCATE (bsubvmns_sur(dimlens(1)), stat = ierror)
1078 CALL cdf_inquire(nwout, vn_bsupumns_sur, dimlens)
1079 ALLOCATE (bsupumns_sur(dimlens(1)), stat = ierror)
1080 CALL cdf_inquire(nwout, vn_bsupvmns_sur, dimlens)
1081 ALLOCATE (bsupvmns_sur(dimlens(1)), stat = ierror)
1087 CALL cdf_read(nwout, vn_pmod, xm)
1088 CALL cdf_read(nwout, vn_tmod, xn)
1089 IF (mnmax_nyq .le. 0)
THEN
1091 ALLOCATE (xm_nyq(mnmax_nyq), xn_nyq(mnmax_nyq), stat=ierror)
1092 xm_nyq = xm; xn_nyq = xn
1094 CALL cdf_read(nwout, vn_pmod_nyq, xm_nyq)
1095 CALL cdf_read(nwout, vn_tmod_nyq, xn_nyq)
1098 IF (mnmaxpot .GT. 0)
THEN
1099 CALL cdf_read(nwout, vn_potsin, potsin)
1100 CALL cdf_read(nwout, vn_xmpot, xmpot)
1101 CALL cdf_read(nwout, vn_xnpot, xnpot)
1102 IF (lasym)
CALL cdf_read(nwout, vn_potcos, potcos)
1105 mnyq = int(maxval(xm_nyq)); nnyq = int(maxval(abs(xn_nyq)))/nfp
1107 CALL cdf_read(nwout, vn_racc, raxis_cc)
1108 CALL cdf_read(nwout, vn_zacs, zaxis_cs)
1110 IF (
SIZE(raxis_cc) .ne. ntor+1)
1111 stop
'WRONG SIZE(raxis_cc) in READ_WOUT_NC'
1112 ALLOCATE (raxis(0:ntor,2), zaxis(0:ntor,2), stat=ierror)
1113 raxis(:,1) = raxis_cc(0:ntor); zaxis(:,1) = zaxis_cs(0:ntor)
1114 raxis(:,2) = 0; zaxis(:,2) = 0
1115 DEALLOCATE (raxis_cc, zaxis_cs, stat=ierror)
1117 CALL cdf_read(nwout, vn_rmnc, rmnc)
1118 CALL cdf_read(nwout, vn_zmns, zmns)
1119 CALL cdf_read(nwout, vn_lmns, lmns)
1120 CALL cdf_read(nwout, vn_gmnc, gmnc)
1121 CALL cdf_read(nwout, vn_bmnc, bmnc)
1122 CALL cdf_read(nwout, vn_bsubumnc, bsubumnc)
1123 CALL cdf_read(nwout, vn_bsubvmnc, bsubvmnc)
1124 CALL cdf_read(nwout, vn_bsubsmns, bsubsmns)
1126 CALL cdf_read(nwout, vn_bsupumnc, bsupumnc)
1127 CALL cdf_read(nwout, vn_bsupvmnc, bsupvmnc)
1129 IF (version_ .ge. 9.0)
THEN
1130 CALL cdf_read(nwout, vn_currumnc, currumnc)
1131 CALL cdf_read(nwout, vn_currvmnc, currvmnc)
1135 CALL cdf_read(nwout, vn_bsubumnc_sur, bsubumnc_sur)
1136 CALL cdf_read(nwout, vn_bsubvmnc_sur, bsubvmnc_sur)
1137 CALL cdf_read(nwout, vn_bsupumnc_sur, bsupumnc_sur)
1138 CALL cdf_read(nwout, vn_bsupvmnc_sur, bsupvmnc_sur)
1142 CALL cdf_read(nwout, vn_racs, raxis_cs)
1143 CALL cdf_read(nwout, vn_zacc, zaxis_cc)
1144 raxis(:,2) = raxis_cs; zaxis(:,2) = zaxis_cc
1145 DEALLOCATE (raxis_cs, zaxis_cc, stat=ierror)
1146 CALL cdf_read(nwout, vn_rmns, rmns)
1147 CALL cdf_read(nwout, vn_zmnc, zmnc)
1148 CALL cdf_read(nwout, vn_lmnc, lmnc)
1149 CALL cdf_read(nwout, vn_gmns, gmns)
1150 CALL cdf_read(nwout, vn_bmns, bmns)
1151 CALL cdf_read(nwout, vn_bsubumns, bsubumns)
1152 CALL cdf_read(nwout, vn_bsubvmns, bsubvmns)
1153 CALL cdf_read(nwout, vn_bsubsmnc, bsubsmnc)
1155 CALL cdf_read(nwout, vn_bsupumns, bsupumns)
1156 CALL cdf_read(nwout, vn_bsupvmns, bsupvmns)
1158 IF (version_ .ge. 9.0)
THEN
1159 CALL cdf_read(nwout, vn_currumns, currumns)
1160 CALL cdf_read(nwout, vn_currvmns, currvmns)
1164 CALL cdf_read(nwout, vn_bsubumns_sur, bsubumns_sur)
1165 CALL cdf_read(nwout, vn_bsubvmns_sur, bsubvmns_sur)
1166 CALL cdf_read(nwout, vn_bsupumns_sur, bsupumns_sur)
1167 CALL cdf_read(nwout, vn_bsupvmns_sur, bsupvmns_sur)
1171 CALL cdf_read(nwout, vn_am, am)
1172 CALL cdf_read(nwout, vn_ac, ac)
1173 CALL cdf_read(nwout, vn_ai, ai)
1175 CALL cdf_read(nwout, vn_am_aux_s, am_aux_s)
1176 CALL cdf_read(nwout, vn_am_aux_f, am_aux_f)
1177 CALL cdf_read(nwout, vn_ac_aux_s, ac_aux_s)
1178 CALL cdf_read(nwout, vn_ac_aux_f, ac_aux_f)
1179 CALL cdf_read(nwout, vn_ai_aux_s, ai_aux_s)
1180 CALL cdf_read(nwout, vn_ai_aux_f, ai_aux_f)
1182 CALL cdf_read(nwout, vn_iotaf, iotaf)
1183 CALL cdf_read(nwout, vn_qfact, qfact)
1185 CALL cdf_read(nwout, vn_presf, presf)
1186 CALL cdf_read(nwout, vn_phi, phi)
1187 CALL cdf_read(nwout, vn_phipf, phipf)
1188 CALL cdf_read(nwout, vn_chi, chi)
1190 CALL cdf_read(nwout, vn_chipf, chipf)
1191 IF (all(chipf .EQ. 1.e30_dp))
THEN
1195 CALL cdf_read(nwout, vn_jcuru, jcuru)
1196 CALL cdf_read(nwout, vn_jcurv, jcurv)
1201 CALL cdf_read(nwout, vn_iotah, iotas)
1202 CALL cdf_read(nwout, vn_mass, mass)
1203 CALL cdf_read(nwout, vn_presh, pres)
1204 CALL cdf_read(nwout, vn_betah, beta_vol)
1205 CALL cdf_read(nwout, vn_buco, buco)
1206 CALL cdf_read(nwout, vn_bvco, bvco)
1207 CALL cdf_read(nwout, vn_vp, vp)
1208 CALL cdf_read(nwout, vn_specw, specw)
1209 CALL cdf_read(nwout, vn_phip, phip)
1210 CALL cdf_read(nwout, vn_jdotb, jdotb)
1211 CALL cdf_read(nwout, vn_bdotb, bdotb)
1212 CALL cdf_read(nwout, vn_bgrv, bdotgradv)
1215 CALL cdf_read(nwout, vn_merc, dmerc)
1216 CALL cdf_read(nwout, vn_mshear, dshear)
1217 CALL cdf_read(nwout, vn_mwell, dwell)
1218 CALL cdf_read(nwout, vn_mcurr, dcurr)
1219 CALL cdf_read(nwout, vn_mgeo, dgeod)
1220 CALL cdf_read(nwout, vn_equif, equif)
1222 CALL cdf_read(nwout, vn_fsq, fsqt)
1223 CALL cdf_read(nwout, vn_wdot, wdot)
1225 IF (nextcur .gt. 0)
THEN
1226 CALL cdf_read(nwout, vn_extcur, extcur)
1227 CALL cdf_read(nwout, vn_curlab, curlabel)
1232 CALL cdf_close(nwout, ierr)
1234 IF (.not.
ALLOCATED(bsubumnc))
RETURN
1242 IF (version_ .lt. 9.0)
THEN
1243 IF (ierror .eq. 0)
THEN
1244 CALL compute_currents(bsubsmnc, bsubsmns,
1245 & bsubumnc, bsubumns,
1246 & bsubvmnc, bsubvmns,
1247 & xm_nyq, xn_nyq, mnmax_nyq, lasym, ns,
1248 & currumnc, currvmnc,
1249 & currumns, currvmns)
1253 IF (ierr .ne. 0) print *,
'in read_wout_nc ierr=',ierr
1254 IF (ierror .ne. 0) print *,
'in read_wout_nc ierror=',ierror
1256 END SUBROUTINE read_wout_nc
1259 SUBROUTINE write_wout_text(filename, ierr)
1261 USE vsvd0,
ONLY: nparts
1263 USE stel_constants,
ONLY: mu0
1269 CHARACTER (len=*) :: filename
1270 INTEGER,
INTENT(out) :: ierr
1274 REAL(rprec),
PARAMETER :: eps_w = 1.e-4_dp
1278 INTEGER :: iounit, js, mn, i, j, k, m, n, iasymm
1288 CALL safe_open(iounit, ierr,
1289 &
'wout_' // trim(filename) //
'.txt',
1290 &
'replace',
'formatted')
1292 CALL assert_eq(0, ierr,
'Error opening text wout file in ' //
1293 &
'write_wout_text of read_wout_mod.')
1297 WRITE (iounit,
'(a15,f4.2)')
'VMEC VERSION = ', version_
1306 IF (version_ .le. (5.10 + eps_w))
THEN
1307 WRITE (iounit, *) wb, wp, gamma, pfac, nfp, ns, mpol, ntor,
1308 & mnmax, itfsq, niter, iasymm, ireconstruct
1310 IF (version_ .lt. 6.54)
THEN
1311 WRITE (iounit, *) wb, wp, gamma, pfac, rmax_surf, rmin_surf
1313 WRITE (iounit, *) wb, wp, gamma, pfac, rmax_surf, rmin_surf,
1316 IF (version_ .le. (8.0 + eps_w))
THEN
1317 WRITE (iounit, *) nfp, ns, mpol, ntor, mnmax, itfsq, niter,
1318 & iasym, ireconstruct, ierr_vmec
1320 WRITE (iounit, *) nfp, ns, mpol, ntor, mnmax, mnmax_nyq,
1321 & itfsq, niter, iasym, ireconstruct,
1326 IF (version_ .gt. (6.20 + eps_w))
THEN
1327 WRITE (iounit, *) imse, itse, nbsets, nobd, nextcur, nstore_seq
1329 WRITE (iounit, *) imse, itse, nbsets, nobd, nextcur
1332 IF (ierr_vmec .ne. norm_term_flag .and.
1333 & ierr_vmec .ne. more_iter_flag)
THEN
1337 IF (nbsets .gt. 0)
THEN
1338 WRITE (iounit, *) nbfld(1:nbsets)
1340 WRITE (iounit, *) trim(mgrid_file)
1345 WRITE (iounit, *) nint(xm(mn)), nint(xn(mn)/nfp)
1347 IF (version_ .le. (6.20 + eps_w))
THEN
1348 WRITE (iounit, 730) rmnc(mn,js), zmns(mn,js),
1349 & lmns(mn,js), bmnc(mn,js),
1350 & gmnc(mn,js), bsubumnc(mn,js),
1351 & bsubvmnc(mn,js), bsubsmns(mn,js),
1352 & bsupumnc(mn,js), bsupvmnc(mn,js),
1354 ELSE IF (version_ .le. (8.0 + eps_w))
THEN
1355 WRITE (iounit, *) rmnc(mn,js), zmns(mn,js), lmns(mn,js),
1356 & bmnc(mn,js), gmnc(mn,js),
1357 & bsubumnc(mn,js), bsubvmnc(mn,js),
1358 & bsubsmns(mn,js), bsupumnc(mn,js),
1359 & bsupvmnc(mn,js), currvmnc(mn,js)
1361 WRITE (iounit, *) rmnc(mn,js), zmns(mn,js), lmns(mn,js)
1366 IF (version_ .le. (8.0 + eps_w))
THEN
1367 WRITE (iounit, *) rmns(mn,js), zmnc(mn,js),
1368 & lmnc(mn,js), bmns(mn,js),
1369 & gmns(mn,js), bsubumns(mn,js),
1370 & bsubvmns(mn,js), bsubsmnc(mn,js),
1371 & bsupumns(mn,js), bsubvmns(mn,js)
1373 WRITE (iounit, *) rmns(mn,js), zmnc(mn,js),
1379 IF (version_ .le. (8.0 + eps_w))
THEN
1383 DO mn = 1, mnmax_nyq
1385 WRITE (iounit, *) nint(xm_nyq(mn)),
1386 & nint(xn_nyq(mn)/nfp)
1388 WRITE (iounit, *) bmnc(mn,js), gmnc(mn,js),
1389 & bsubumnc(mn,js), bsubvmnc(mn,js),
1390 & bsubsmns(mn,js), bsupumnc(mn,js),
1393 WRITE (iounit, *) bmns(mn,js), gmns(mn,js),
1394 & bsubumns(mn,js), bsubvmns(mn,js),
1395 & bsubsmnc(mn,js), bsupumns(mn,js),
1411 IF (version_ .le. (6.05 + eps_w))
THEN
1412 WRITE (iounit, 730) (iotas(js), mass(js)*mu0, pres(js)*mu0,
1413 & phip(js), buco(js), bvco(js), -phi(js),
1414 & vp(js), overr(js), jcuru(js)*mu0,
1415 & jcurv(js)*mu0, specw(js), js=2, ns)
1416 WRITE (iounit, 730) aspect, betatot, betapol, betaxis, b0
1417 ELSE IF (version_ .le. (6.20 + eps_w))
THEN
1418 WRITE (iounit, 730) (iotas(js), mass(js), pres(js),
1419 & beta_vol(js), phip(js), buco(js),
1420 & bvco(js), phi(js), vp(js), overr(js),
1421 & jcuru(js), jcurv(js), specw(js),
1423 WRITE (iounit, 730) aspect, betatot, betapol, betaxis, b0
1424 ELSE IF (version_ .le. (6.95 + eps_w))
THEN
1425 WRITE (iounit, *) (iotas(js), mass(js), pres(js),
1426 & beta_vol(js), phip(js), buco(js),
1427 & bvco(js), phi(js), vp(js), overr(js),
1428 & jcuru(js), jcurv(js), specw(js),
1430 WRITE (iounit, *) aspect, betatot, betapol, betaxis, b0
1432 WRITE (iounit, *) (iotaf(js), presf(js), phipf(js), phi(js),
1433 & jcuru(js), jcurv(js), js=1, ns)
1434 WRITE (iounit, *) (iotas(js), mass(js), pres(js),
1435 & beta_vol(js), phip(js), buco(js),
1436 & bvco(js), vp(js), overr(js), specw(js),
1438 WRITE (iounit, *) aspect, betatot, betapol, betaxis, b0
1441 IF (version_ .gt. (6.10 + eps_w))
THEN
1442 WRITE (iounit, *) isigng
1443 WRITE (iounit, *) trim(input_extension)
1444 WRITE (iounit, *) ionlarmor, volavgb, rbtor0, rbtor, itor,
1445 & aminor, rmajor, volume
1451 IF (version_ .gt. (5.10 + eps_w) .and.
1452 & version_ .lt. (6.20 - eps_w))
THEN
1453 WRITE (iounit, 730) (dmerc(js), dshear(js), dwell(js),
1454 & dcurr(js), dgeod(js), equif(js),
1456 ELSE IF (version_ .ge. (6.20 - eps_w))
THEN
1457 WRITE (iounit, *) (dmerc(js), dshear(js), dwell(js),
1458 & dcurr(js), dgeod(js), equif(js),
1462 IF (nextcur .gt. 0)
THEN
1463 IF (version_ .le. (6.20 + eps_w))
THEN
1464 WRITE (iounit, 730) (extcur(js), js=1, nextcur)
1466 WRITE (iounit, *) (extcur(js), js=1, nextcur)
1469 lcurr = len_trim(curlabel(1)) .gt. 0
1470 WRITE (iounit, *) lcurr
1472 WRITE (iounit, *) (trim(curlabel(js)), js=1, nextcur)
1476 IF (version_ .le. (6.20 + eps_w))
THEN
1477 WRITE (iounit, 730) (fsqt(js), wdot(js), js = 1, nstore_seq)
1479 WRITE (iounit, *) (fsqt(js), wdot(js), js = 1, nstore_seq)
1482 IF (version_ .ge. (6.20 - eps_w) .and.
1483 & version_ .lt. (6.50 - eps_w))
THEN
1484 WRITE (iounit, 730) (jdotb(js), bdotgradv(js), js=1, ns)
1485 ELSE IF (version_ .ge. (6.50 - eps_w))
THEN
1486 WRITE (iounit, *) (jdotb(js), bdotgradv(js), js=1, ns)
1492 IF (ireconstruct .gt. 0)
THEN
1493 IF (imse .ge. 2 .or. itse .gt. 0)
THEN
1494 WRITE (iounit, *) tswgt, msewgt
1495 WRITE (iounit, *) isnodes, (sknots(js), ystark(js),
1496 & y2stark(js), js=1, isnodes)
1497 WRITE (iounit, *) ipnodes, (pknots(js), ythom(js),
1498 & y2thom(js), js=1, ipnodes)
1499 WRITE (iounit, *) (anglemse(js), rmid(js), qmid(js),
1500 & shear(js), presmid(js), alfa(js),
1501 & curmid(js), js=1, 2*ns - 1)
1502 WRITE (iounit, *) (rstark(js), datastark(js), qmeas(js),
1504 WRITE (iounit, *) (rthom(js), datathom(i), js=1, itse)
1507 IF (nobd .gt. 0)
THEN
1508 WRITE (iounit, *) (dsiext(js), plflux(js), dsiobt(js),
1510 WRITE (iounit, *) flmwgt
1513 IF (nbfldn .gt. 0)
THEN
1515 READ (iounit, *) (bcoil(i,n), plbfld(i,n), bbc(i,n),
1518 WRITE (iounit, *) bcwgt
1521 WRITE (iounit, *) phidiam, delphid
1525 WRITE (iounit, *) nsets, nparts, nlim
1527 WRITE (iounit, *) (nsetsn(js), js=1, nsets)
1529 WRITE (iounit, *) (((pfcspec(i,j,k), i=1, nparts),
1530 & j=1, nsetsn(k)), k=1, nsets)
1532 WRITE (iounit, *) (limitr(i), i=1, nlim)
1534 WRITE (iounit, *) ((rlim(i,j), zlim(i,j), i=1, limitr(j)),
1536 WRITE (iounit, *) nrgrid, nzgrid
1537 WRITE (iounit, *) tokid
1538 WRITE (iounit, *) rx1, rx2, zy1, zy2, condif
1539 WRITE (iounit, *) imatch_phiedge
1544 WRITE (iounit, *) mgrid_mode
1548 CLOSE (iounit, iostat = ierr)
1549 CALL assert_eq(0, ierr,
'Error closing text wout file in ' //
1550 &
'write_wout_text of read_wout_mod.')
1554 SUBROUTINE compute_currents(bsubsmnc_, bsubsmns_, &
1555 & bsubumnc_, bsubumns_, &
1556 & bsubvmnc_, bsubvmns_, &
1557 & xm_nyq_, xn_nyq_, mnmax_nyq_, &
1559 & currumnc_, currvmnc_, &
1560 & currumns_, currvmns_)
1561 USE stel_constants,
ONLY: mu0
1564 REAL(rprec),
DIMENSION(:,:),
INTENT(in) :: bsubsmnc_
1565 REAL(rprec),
DIMENSION(:,:),
INTENT(in) :: bsubsmns_
1566 REAL(rprec),
DIMENSION(:,:),
INTENT(in) :: bsubumnc_
1567 REAL(rprec),
DIMENSION(:,:),
INTENT(in) :: bsubumns_
1568 REAL(rprec),
DIMENSION(:,:),
INTENT(in) :: bsubvmnc_
1569 REAL(rprec),
DIMENSION(:,:),
INTENT(in) :: bsubvmns_
1571 REAL(rprec),
DIMENSION(:),
INTENT(in) :: xm_nyq_
1572 REAL(rprec),
DIMENSION(:),
INTENT(in) :: xn_nyq_
1574 INTEGER,
INTENT(in) :: mnmax_nyq_
1575 LOGICAL,
INTENT(in) :: lasym_
1576 INTEGER,
INTENT(in) :: ns_
1578 REAL(rprec),
DIMENSION(:,:),
INTENT(out) :: currumnc_
1579 REAL(rprec),
DIMENSION(:,:),
INTENT(out) :: currvmnc_
1580 REAL(rprec),
DIMENSION(:,:),
INTENT(out) :: currumns_
1581 REAL(rprec),
DIMENSION(:,:),
INTENT(out) :: currvmns_
1587 REAL(rprec) :: ohs, hs, shalf(ns_), sfull(ns_)
1588 REAL(rprec),
DIMENSION(mnmax_nyq_) :: bu1, bu0, bv1, bv0, t1, t2, t3
1605 shalf(js) = sqrt(hs*(js-1.5_dp))
1606 sfull(js) = sqrt(hs*(js-1))
1610 WHERE (mod(int(xm_nyq_),2) .EQ. 1)
1611 t1 = 0.5_dp*(shalf(js+1)*bsubsmns_(:,js+1)
1612 + shalf(js) *bsubsmns_(:,js)) /sfull(js)
1613 bu0 = bsubumnc_(:,js )/shalf(js)
1614 bu1 = bsubumnc_(:,js+1)/shalf(js+1)
1615 t2 = ohs*(bu1-bu0)*sfull(js)+0.25_dp*(bu0+bu1)/sfull(js)
1616 bv0 = bsubvmnc_(:,js )/shalf(js)
1617 bv1 = bsubvmnc_(:,js+1)/shalf(js+1)
1618 t3 = ohs*(bv1-bv0)*sfull(js)+0.25_dp*(bv0+bv1)/sfull(js)
1620 t1 = 0.5_dp*(bsubsmns_(:,js+1)+bsubsmns_(:,js))
1621 t2 = ohs*(bsubumnc_(:,js+1)-bsubumnc_(:,js))
1622 t3 = ohs*(bsubvmnc_(:,js+1)-bsubvmnc_(:,js))
1624 currumnc_(:,js) = -xn_nyq_(:)*t1 - t3
1625 currvmnc_(:,js) = -xm_nyq_(:)*t1 + t2
1628 WHERE (xm_nyq_ .LE. 1)
1629 currvmnc_(:,1) = 2*currvmnc_(:,2) - currvmnc_(:,3)
1630 currumnc_(:,1) = 2*currumnc_(:,2) - currumnc_(:,3)
1636 currumnc_(:,ns_) = 2*currumnc_(:,ns_-1) - currumnc_(:,ns_-2)
1637 currvmnc_(:,ns_) = 2*currvmnc_(:,ns_-1) - currvmnc_(:,ns_-2)
1638 currumnc_ = currumnc_ /mu0; currvmnc_ = currvmnc_/mu0
1640 IF (.NOT.lasym_)
RETURN
1643 WHERE (mod(int(xm_nyq_),2) .EQ. 1)
1644 t1 = 0.5_dp*(shalf(js+1)*bsubsmnc_(:,js+1)
1645 + shalf(js) *bsubsmnc_(:,js)) / sfull(js)
1646 bu0 = bsubumns_(:,js )/shalf(js+1)
1647 bu1 = bsubumns_(:,js+1)/shalf(js+1)
1648 t2 = ohs*(bu1-bu0)*sfull(js) + 0.25_dp*(bu0+bu1)/sfull(js)
1649 bv0 = bsubvmns_(:,js )/shalf(js)
1650 bv1 = bsubvmns_(:,js+1)/shalf(js+1)
1651 t3 = ohs*(bv1-bv0)*sfull(js)+0.25_dp*(bv0+bv1)/sfull(js)
1653 t1 = 0.5_dp*(bsubsmnc_(:,js+1) + bsubsmnc_(:,js))
1654 t2 = ohs*(bsubumns_(:,js+1)-bsubumns_(:,js))
1655 t3 = ohs*(bsubvmns_(:,js+1)-bsubvmns_(:,js))
1657 currumns_(:,js) = xn_nyq_(:)*t1 - t3
1658 currvmns_(:,js) = xm_nyq_(:)*t1 + t2
1661 WHERE (xm_nyq_ .LE. 1)
1662 currvmns_(:,1) = 2*currvmns_(:,2) - currvmns_(:,3)
1663 currumns_(:,1) = 2*currumns_(:,2) - currumns_(:,3)
1668 currumns_(:,ns_) = 2*currumns_(:,ns_-1) - currumns_(:,ns_-2)
1669 currvmns_(:,ns_) = 2*currvmns_(:,ns_-1) - currvmns_(:,ns_-2)
1670 currumns_ = currumns_/mu0; currvmns_ = currvmns_/mu0
1672 END SUBROUTINE compute_currents
1674 SUBROUTINE read_wout_deallocate
1679 INTEGER :: istat(10)
1682 lwout_opened=.false.
1684 IF (
ALLOCATED(extcur))
DEALLOCATE (extcur, curlabel,
1686 IF (
ALLOCATED(overr))
DEALLOCATE (overr, stat = istat(2))
1688 IF (
ALLOCATED(xm))
DEALLOCATE (xm, xn, xm_nyq, xn_nyq,
1689 rmnc, zmns, lmns, bmnc, gmnc, bsubumnc, iotaf, presf, phipf,
1691 bsubvmnc, bsubsmns, bsupumnc, bsupvmnc, currvmnc, iotas, mass,
1692 pres, beta_vol, phip, buco, bvco, phi, vp, jcuru, am, ac, ai,
1693 jcurv, specw, dmerc, dshear, dwell, dcurr, dgeod, equif, jdotb,
1694 bdotb, bdotgradv, raxis, zaxis, fsqt, wdot, stat=istat(3))
1696 IF (
ALLOCATED(potsin))
DEALLOCATE (potsin)
1697 IF (
ALLOCATED(potcos))
DEALLOCATE (potcos)
1699 IF (
ALLOCATED(chipf))
DEALLOCATE (chipf, chi)
1701 IF (
ALLOCATED(am_aux_s))
DEALLOCATE (am_aux_s, am_aux_f,
1702 ac_aux_s, ac_aux_f, ai_aux_s, ai_aux_f, stat=istat(6))
1704 IF (ireconstruct.gt.0 .and.
ALLOCATED(sknots))
DEALLOCATE (
1705 ystark, y2stark, pknots, anglemse, rmid, qmid, shear,
1706 presmid, alfa, curmid, rstark, datastark, rthom, datathom,
1707 ythom, y2thom, plflux, dsiobt, bcoil, plbfld, bbc, sknots,
1708 pfcspec, limitr, rlim, zlim, nsetsn, stat = istat(4))
1710 IF (
ALLOCATED(rmns))
DEALLOCATE (rmns, zmnc, lmnc,
1711 bmns, gmns, bsubumns, bsubvmns, bsubsmnc,
1712 bsupumns, bsupvmns, stat=istat(5))
1714 IF (
ALLOCATED(bsubumnc_sur))
THEN
1715 DEALLOCATE(bsubumnc_sur, bsubvmnc_sur,
1716 bsupumnc_sur, bsupvmnc_sur)
1719 IF (
ALLOCATED(bsubumns_sur))
THEN
1720 DEALLOCATE(bsubumns_sur, bsubvmns_sur,
1721 bsupumns_sur, bsupvmns_sur)
1725 IF (
ALLOCATED(currumnc))
DEALLOCATE (currumnc)
1726 IF (
ALLOCATED(currumns))
DEALLOCATE (currumns, currvmns)
1727 IF (
ALLOCATED(rzl_local))
DEALLOCATE (rzl_local)
1729 IF (any(istat .ne. 0))
1730 stop
'Deallocation error in read_wout_deallocate'
1732 END SUBROUTINE read_wout_deallocate
1734 SUBROUTINE tosuvspace (s_in, u_in, v_in, gsqrt, &
1735 bsupu, bsupv, jsupu, jsupv)
1736 USE stel_constants,
ONLY: zero, one
1741 REAL(rprec),
INTENT(in) :: s_in, u_in, v_in
1742 REAL(rprec),
INTENT(out),
OPTIONAL :: gsqrt, bsupu, bsupv, &
1747 REAL(rprec),
PARAMETER :: c1p5 = 1.5_dp
1748 INTEGER :: m, n, n1, mn, ipresent, jslo, jshi
1749 REAL(rprec) :: hs1, wlo, whi, wlo_odd, whi_odd
1750 REAL(rprec),
DIMENSION(mnmax_nyq) :: gmnc1, gmns1, bsupumnc1, &
1751 bsupumns1, bsupvmnc1, bsupvmns1, jsupumnc1, jsupumns1, &
1752 jsupvmnc1, jsupvmns1, wmins, wplus
1753 REAL(rprec) :: cosu, sinu, cosv, sinv, tcosmn, tsinmn, sgn
1754 REAL(rprec) :: cosmu(0:mnyq), sinmu(0:mnyq), &
1755 cosnv(0:nnyq), sinnv(0:nnyq)
1756 LOGICAL :: lgsqrt, lbsupu, lbsupv, ljsupu, ljsupv
1772 IF (s_in.lt.zero .or. s_in.gt.one)
THEN
1773 WRITE(6, *)
' In tosuvspace, s(flux) must be between 0 and 1'
1777 IF (.not.lwout_opened)
THEN
1778 WRITE(6, *)
' tosuvspace can only be called AFTER opening wout file!'
1785 cosu = cos(u_in); sinu = sin(u_in)
1786 cosv = cos(v_in); sinv = sin(v_in)
1788 cosmu(0) = 1; sinmu(0) = 0
1789 cosnv(0) = 1; sinnv(0) = 0
1791 cosmu(m) = cosmu(m-1)*cosu - sinmu(m-1)*sinu
1792 sinmu(m) = sinmu(m-1)*cosu + cosmu(m-1)*sinu
1796 cosnv(n) = cosnv(n-1)*cosv - sinnv(n-1)*sinv
1797 sinnv(n) = sinnv(n-1)*cosv + cosnv(n-1)*sinv
1807 jslo = int(c1p5 + s_in/hs1)
1809 wlo = (hs1*(jshi-c1p5) - s_in)/hs1
1811 IF (jslo .eq. ns)
THEN
1815 wlo = 1+whi; whi = -whi
1816 ELSE IF (jslo .eq. 1)
THEN
1823 whi_odd = whi*sqrt(s_in/(hs1*(jshi-c1p5)))
1824 IF (jslo .ne. 1)
THEN
1825 wlo_odd = wlo*sqrt(s_in/(hs1*(jslo-c1p5)))
1828 whi_odd = sqrt(s_in/(hs1*(jshi-c1p5)))
1831 WHERE (mod(nint(xm_nyq(:)),2) .eq. 0)
1840 lgsqrt =
PRESENT(gsqrt)
1842 gsqrt = 0 ; ipresent = ipresent+1
1843 gmnc1 = wmins*gmnc(:,jslo) + wplus*gmnc(:,jshi)
1844 IF (lasym) gmns1 = wmins*gmns(:,jslo) + wplus*gmns(:,jshi)
1846 lbsupu =
PRESENT(bsupu)
1848 bsupu = 0 ; ipresent = ipresent+1
1849 bsupumnc1 = wmins*bsupumnc(:,jslo) + wplus*bsupumnc(:,jshi)
1850 IF (lasym) bsupumns1 = wmins*bsupumns(:,jslo) + wplus*bsupumns(
1852 lbsupv =
PRESENT(bsupv)
1854 bsupv = 0 ; ipresent = ipresent+1
1855 bsupvmnc1 = wmins*bsupvmnc(:,jslo) + wplus*bsupvmnc(:,jshi)
1856 IF (lasym) bsupvmns1 = wmins*bsupvmns(:,jslo) + wplus*bsupvmns(
1859 IF (ipresent .eq. 0)
GOTO 1000
1865 DO mn = 1, mnmax_nyq
1866 m = nint(xm_nyq(mn)); n = nint(xn_nyq(mn))/nfp
1867 n1 = abs(n); sgn = sign(1,n)
1868 tcosmn = cosmu(m)*cosnv(n1) + sgn*sinmu(m)*sinnv(n1)
1869 IF (lgsqrt) gsqrt = gsqrt + gmnc1(mn)*tcosmn
1870 IF (lbsupu) bsupu = bsupu + bsupumnc1(mn)*tcosmn
1871 IF (lbsupv) bsupv = bsupv + bsupvmnc1(mn)*tcosmn
1874 IF (.not.lasym)
GOTO 1000
1876 DO mn = 1, mnmax_nyq
1877 m = nint(xm_nyq(mn)); n = nint(xn_nyq(mn))/nfp
1878 n1 = abs(n); sgn = sign(1,n)
1879 tsinmn = sinmu(m)*cosnv(n1) - sgn*cosmu(m)*sinnv(n1)
1880 IF (lgsqrt) gsqrt = gsqrt + gmns1(mn)*tsinmn
1881 IF (lbsupu) bsupu = bsupu + bsupumns1(mn)*tsinmn
1882 IF (lbsupv) bsupv = bsupv + bsupvmns1(mn)*tsinmn
1894 jslo = 1+int(s_in/hs1)
1896 IF (jslo .eq. ns) jshi = ns
1897 wlo = (hs1*(jshi-1) - s_in)/hs1
1902 whi_odd = whi*sqrt(s_in/(hs1*(jshi-1)))
1903 IF (jslo .ne. 1)
THEN
1904 wlo_odd = wlo*sqrt(s_in/(hs1*(jslo-1)))
1907 whi_odd = sqrt(s_in/(hs1*(jshi-1)))
1910 WHERE (mod(nint(xm_nyq(:)),2) .eq. 0)
1919 ljsupu =
PRESENT(jsupu)
1921 IF (.not.lgsqrt) stop
'MUST compute gsqrt for jsupu'
1922 jsupu = 0 ; ipresent = ipresent+1
1923 jsupumnc1 = wmins*currumnc(:,jslo) + wplus*currumnc(:,jshi)
1924 IF (lasym) jsupumns1 = wmins*currumns(:,jslo) + wplus*currumns(
1927 ljsupv =
PRESENT(jsupv)
1929 IF (.not.lgsqrt) stop
'MUST compute gsqrt for jsupv'
1930 jsupv = 0 ; ipresent = ipresent+1
1931 jsupvmnc1 = wmins*currvmnc(:,jslo) + wplus*currvmnc(:,jshi)
1932 IF (lasym) jsupvmns1 = wmins*currvmns(:,jslo) + wplus*currvmns(
1935 IF (ipresent .eq. 0)
RETURN
1937 DO mn = 1, mnmax_nyq
1938 m = nint(xm_nyq(mn)); n = nint(xn_nyq(mn))/nfp
1939 n1 = abs(n); sgn = sign(1,n)
1940 tcosmn = cosmu(m)*cosnv(n1) + sgn*sinmu(m)*sinnv(n1)
1941 IF (ljsupu) jsupu = jsupu + jsupumnc1(mn)*tcosmn
1942 IF (ljsupv) jsupv = jsupv + jsupvmnc1(mn)*tcosmn
1945 IF (.not.lasym)
GOTO 2000
1947 DO mn = 1, mnmax_nyq
1948 m = nint(xm_nyq(mn)); n = nint(xn_nyq(mn))/nfp
1949 n1 = abs(n); sgn = sign(1,n)
1950 tsinmn = sinmu(m)*cosnv(n1) - sgn*cosmu(m)*sinnv(n1)
1951 IF (ljsupu) jsupu = jsupu + jsupumns1(mn)*tsinmn
1952 IF (ljsupv) jsupv = jsupv + jsupvmns1(mn)*tsinmn
1957 IF (ljsupu) jsupu = jsupu/gsqrt
1958 IF (ljsupv) jsupv = jsupv/gsqrt
1960 END SUBROUTINE tosuvspace
1967 INTEGER :: rcc, rss, zsc, zcs, rsc, rcs, zcc, zss
1968 INTEGER :: mpol1, mn, m, n, n1
1976 IF (
ALLOCATED(rzl_local))
RETURN
1980 IF (.not.lasym)
THEN
1990 rss = 2; rsc = 3; rcs = 4
1991 zcs = 2; zcc = 3; zss = 4
2000 zsc = 1+ntmax; zcs = zcs+ntmax; zcc = zcc+ntmax; zss = zss+ntmax
2001 ALLOCATE(rzl_local(ns,0:ntor,0:mpol1,3*ntmax), stat=n)
2002 IF (n .ne. 0) stop
'Allocation error in LoadRZL'
2006 m = nint(xm(mn)); n = nint(xn(mn))/nfp; n1 = abs(n)
2008 rzl_local(:,n1,m,rcc) = rzl_local(:,n1,m,rcc) + rmnc(mn,:)
2009 rzl_local(:,n1,m,zsc) = rzl_local(:,n1,m,zsc) + zmns(mn,:)
2011 rzl_local(:,n1,m,rss) = rzl_local(:,n1,m,rss) + sgn*rmnc(mn,
2012 rzl_local(:,n1,m,zcs) = rzl_local(:,n1,m,zcs) - sgn*zmns(mn,
2015 rzl_local(:,n1,m,rsc) = rzl_local(:,n1,m,rsc) + rmns(mn,:)
2016 rzl_local(:,n1,m,zcc) = rzl_local(:,n1,m,zcc) + zmnc(mn,:)
2018 rzl_local(:,n1,m,rcs) = rzl_local(:,n1,m,rcs)
2020 rzl_local(:,n1,m,zss) = rzl_local(:,n1,m,zss)
2026 END SUBROUTINE loadrzl
2028 END MODULE read_wout_mod