V3FIT
read_wout_mod.f90
1  MODULE read_wout_mod
2 !
3 ! USE READ_WOUT_MOD to include variables dynamically allocated
4 ! in this module
5 ! Call DEALLOCATE_READ_WOUT to free this memory when it is no longer needed
6 !
7 ! Reads in output from VMEC equilibrium code(s), contained in wout file
8 !
9 ! Contained subroutines:
10 !
11 ! read_wout_file wrapper alias called to read/open wout file
12 ! read_wout_text called by read_wout_file to read text file wout
13 ! read_wout_nc called by read_wout_file to read netcdf file wout
14 !
15 ! Post-processing routines
16 !
17 ! mse_pitch user-callable function to compute mse pitch angle
18 ! for the computed equilibrium
19 !
20 
21  USE vmec_input, ONLY: lrfp, lmove_axis, nbfld
22  USE mgrid_mod
23 
24  IMPLICIT NONE
25 #if defined(NETCDF)
26 !------------------------------------------------
27 ! L O C A L P A R A M E T E R S
28 !------------------------------------------------
29 ! Variable names (vn_...) : put eventually into library, used by read_wout too...
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', &
40 
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', &
61 
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', &
73  vn_fsqz = 'fsqz', &
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'
95 ! Long names (ln_...)
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', &
142 
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', &
147 
148  ln_chipf = 'd(chi)/ds: Poroidal flux deriv on full mesh', &
149 
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', &
169 
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', &
178 
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', &
183 
184  ln_bsupumnc = 'BSUPUmnc half', ln_bsupvmnc = 'BSUPVmnc half', &
185 
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', &
194 
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', &
199 
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', &
204 
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'
210 #endif
211 !-----------------------------------------------
212 ! L o c a l V a r i a b l e s
213 !-----------------------------------------------
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, &
217  mnmaxpot
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
248 
249  INTEGER, PARAMETER :: norm_term_flag=0, &
250  bad_jacobian_flag=1, more_iter_flag=2, jac75_flag=4
251 
252 ! OVERLOAD SUBROUTINE READ_WOUT_FILE TO ACCEPT BOTH UNIT NO. (OPENED EXTERNALLY)
253 ! OR FILENAME (HANDLE OPEN/CLOSE HERE)
254  INTERFACE read_wout_file
255  MODULE PROCEDURE readw_and_open, readw_only
256  END INTERFACE
257 
258 #if defined(NETCDF)
259  PRIVATE :: read_wout_text, read_wout_nc
260 #else
261  PRIVATE :: read_wout_text
262 #endif
263  PRIVATE :: norm_term_flag, bad_jacobian_flag, &
264  more_iter_flag, jac75_flag
265 
266  CONTAINS
267 
268  SUBROUTINE readw_and_open(file_or_extension, ierr, iopen)
269  USE safe_open_mod
270  IMPLICIT NONE
271 !------------------------------------------------
272 ! D u m m y A r g u m e n t s
273 !------------------------------------------------
274  INTEGER, INTENT(out) :: ierr
275  INTEGER, OPTIONAL :: iopen
276  CHARACTER(LEN=*), INTENT(in) :: file_or_extension
277 !------------------------------------------------
278 ! L o c a l V a r i a b l e s
279 !------------------------------------------------
280  INTEGER, PARAMETER :: iunit_init = 10
281  INTEGER :: iunit
282  LOGICAL :: isnc
283  CHARACTER(len=LEN_TRIM(file_or_extension)+10) :: filename
284 !------------------------------------------------
285 !
286 ! THIS SUBROUTINE READS THE WOUT FILE CREATED BY THE VMEC CODE
287 ! AND STORES THE DATA IN THE READ_WOUT MODULE
288 !
289 ! FIRST, CHECK IF THIS IS A FULLY-QUALIFIED PATH NAME
290 ! MAKE SURE wout IS NOT EMBEDDED IN THE NAME (PERVERSE USER...)
291 !
292  filename = 'wout'
293  CALL parse_extension(filename, file_or_extension, isnc)
294  CALL flush(6)
295 !SPH IF (.not.isnc) STOP 'ISNC ERR IN READ_WOUT_MOD'
296  IF (isnc) THEN
297 #if defined(NETCDF)
298  CALL read_wout_nc(filename, ierr)
299 #else
300  print *, "NETCDF wout file can not be opened on this platform"
301  ierr = -100
302 #endif
303  ELSE
304  iunit = iunit_init
305  CALL safe_open (iunit, ierr, filename, 'old', 'formatted')
306  IF (ierr .eq. 0) CALL read_wout_text(iunit, ierr)
307  CLOSE(unit=iunit)
308  END IF
309 
310  IF (PRESENT(iopen)) iopen = ierr
311  lwout_opened = (ierr .eq. 0)
312  ! WHEN READING A NETCDF FILE, A BAD RUN MAY PREVENT XN FROM BEING
313  ! READ, SUBSEQUENTLY WE MUST CHECK TO SEE IF XN HAS BEEN ALLOCATED
314  ! BEFORE DOING ANYTHING WITH IT OTHERWISE WE DEFAULT LTHREED TO
315  ! FALSE. - SAL 09/07/11
316  IF (ALLOCATED(xn)) THEN
317  lthreed = any(nint(xn) .ne. 0)
318  ELSE
319  lthreed = .false.
320  END IF
321 
322  END SUBROUTINE readw_and_open
323 
324 
325  SUBROUTINE readw_only(iunit, ierr, iopen)
326  IMPLICIT NONE
327 !------------------------------------------------
328 ! D u m m y A r g u m e n t s
329 !------------------------------------------------
330  INTEGER, INTENT(in) :: iunit
331  INTEGER, INTENT(out):: ierr
332  INTEGER, OPTIONAL :: iopen
333 !------------------------------------------------
334 ! L o c a l V a r i a b l e s
335 !------------------------------------------------
336  INTEGER :: istat
337  CHARACTER(LEN=256) :: vmec_version
338  LOGICAL :: exfile
339 !------------------------------------------------
340 !
341 ! User opened the file externally and has a unit number, iunit
342 !
343  ierr = 0
344 
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
350  ierr = -1
351  RETURN
352  ELSE
353  IF (PRESENT(iopen)) iopen = 0
354  END IF
355 
356  CALL read_wout_text(iunit, ierr)
357  lwout_opened = (ierr .eq. 0)
358  lthreed = any(nint(xn) .ne. 0)
359 
360  END SUBROUTINE readw_only
361 
362 
363  SUBROUTINE read_wout_text(iunit, ierr)
364  USE stel_constants, ONLY: mu0
365  IMPLICIT NONE
366 !------------------------------------------------
367 ! D u m m y A r g u m e n t s
368 !------------------------------------------------
369  INTEGER :: iunit, ierr
370 !------------------------------------------------
371 ! L o c a l P a r a m e t e r s
372 !------------------------------------------------
373  REAL(rprec), PARAMETER :: eps_w = 1.e-4_dp
374 !------------------------------------------------
375 ! L o c a l V a r i a b l e s
376 !------------------------------------------------
377  INTEGER :: istat(15), i, j, k, js, m, n, n1, mn, nparts_in
378  CHARACTER(LEN=256) :: vmec_version
379  LOGICAL :: lcurr
380 !------------------------------------------------
381 !
382 ! THIS SUBROUTINE READS THE TEXT FILE WOUT CREATED BY THE VMEC CODE
383 ! AND STORES THE INFORMATION IN THE read_WOUT MODULE
384 !
385 ! CALL read_wout_file - GENERIC INTERFACE - CAN BE CALLED WITH EITHER UNIT NO. OR FILENAME
386 !
387 ! RMNC, ZMNS: FULL-GRID
388 ! LMNS : HALF-GRID
389 !
390  istat = 0
391  ierr = 0
392  nextcur = 0
393 
394  READ (iunit, '(a)', iostat=istat(2), err=1000) vmec_version
395 
396  i = index(vmec_version,'=')
397  IF (i .ge. 0) THEN
398  READ(vmec_version(i+1:len_trim(vmec_version)),*) version_
399  ELSE
400  version_ = -1.0
401  END IF
402 
403  ierr_vmec = norm_term_flag
404 
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, &
408  ireconstruct
409  ELSE
410  IF (version_ .lt. 6.54) THEN
411  READ (iunit, *, iostat=istat(2), err=1000) wb, wp, gamma, &
412  pfac, rmax_surf, rmin_surf
413  ELSE
414  READ (iunit, *, iostat=istat(2), err=1000) wb, wp, gamma, &
415  pfac, rmax_surf, rmin_surf, zmax_surf
416  END IF
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
420  mnmax_nyq = mnmax
421  ELSE
422  READ (iunit, *, iostat=istat(2), err=1000) nfp, ns, mpol, &
423  ntor, mnmax, mnmax_nyq, itfsq, niter, iasym, ireconstruct, &
424  ierr_vmec
425  END IF
426  END IF
427 
428  lasym = (iasym .gt. 0)
429 
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
433  ELSE
434  READ (iunit, *, iostat=istat(1), err=1000)imse, itse, nbsets, &
435  nobd, nextcur
436  nstore_seq = 100
437  END IF
438 
439  IF (ierr_vmec.ne.norm_term_flag .and. &
440  ierr_vmec.ne.more_iter_flag) GOTO 1000
441 
442  IF (nextcur .gt. nigroup) istat(15) = -1
443 
444  IF (ALLOCATED(xm)) CALL read_wout_deallocate
445 
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), &
455  chipf(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))
461 
462  IF (nextcur .GT. 0) ALLOCATE(extcur(nextcur), curlabel(nextcur), &
463  stat = istat(6))
464 
465  IF (lasym) &
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), &
472  stat=istat(6))
473 
474  fsqt = 0; wdot = 0; raxis = 0; zaxis = 0
475 
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
479 
480  DO js = 1, ns
481  DO mn = 1, mnmax
482  IF(js .eq. 1) THEN
483  READ (iunit, *, iostat=istat(7), err=1000) m, n
484  xm(mn) = real(m,rprec)
485  xn(mn) = real(n,rprec)
486  END IF
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), &
492  currvmnc(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), &
498  currvmnc(mn,js)
499  ELSE
500  READ (iunit, *, iostat=istat(8), err=1000) &
501  rmnc(mn,js), zmns(mn,js), lmns(mn,js)
502  END IF
503  IF (lasym) THEN
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)
510  ELSE
511  READ (iunit, *, iostat=istat(8), err=1000) &
512  rmns(mn,js), zmnc(mn,js), lmnc(mn,js)
513  END IF
514  END IF
515  IF (js.eq.1 .and. m.eq.0) THEN
516  n1 = abs(n/nfp)
517  IF (n1 .le. ntor) THEN
518  raxis(n1,1) = rmnc(mn,1)
519  zaxis(n1,1) = zmns(mn,1)
520  IF (lasym) THEN
521  raxis(n1,2) = rmns(mn,1)
522  zaxis(n1,2) = zmnc(mn,1)
523  END IF
524  END IF
525  END IF
526  END DO
527 
528  IF (version_ .le. (8.0+eps_w)) cycle
529  DO mn = 1, mnmax_nyq
530  IF(js .eq. 1) THEN
531  READ (iunit, *, iostat=istat(7), err=1000) m, n
532  xm_nyq(mn) = real(m,rprec)
533  xn_nyq(mn) = real(n,rprec)
534  END IF
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)
539  IF (lasym) THEN
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)
544  END IF
545  END DO
546 
547  END DO
548 
549 ! Compute current coefficients on full mesh
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)
557  END IF
558 
559  mnyq = int(maxval(xm_nyq)); nnyq = int(maxval(abs(xn_nyq)))/nfp
560 
561 !
562 ! Read FULL AND HALF-MESH QUANTITIES
563 !
564 ! NOTE: In version_ <= 6.00, mass, press were written out in INTERNAL (VMEC) units
565 ! and are therefore multiplied here by 1/mu0 to transform to pascals. Same is true
566 ! for ALL the currents (jcuru, jcurv, jdotb). Also, in version_ = 6.10 and
567 ! above, PHI is the true (physical) toroidal flux (has the sign of jacobian correctly
568 ! built into it)
569 !
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
572  beta_vol(1) = 0
573 
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
595  ELSE
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
605  END IF
606 
607 
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
613  END IF
614 
615 !-----------------------------------------------
616 ! MERCIER CRITERION
617 !-----------------------------------------------
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)
626  END IF
627 
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)
632  ELSE
633  READ (iunit, *, iostat=istat(12), err=1000) &
634  (extcur(i),i=1,nextcur)
635  END IF
636  READ (iunit, *, iostat=istat(13)) lcurr
637  IF (lcurr) READ (iunit, *, iostat=istat(13), err=1000) &
638  (curlabel(i),i=1,nextcur)
639  END IF
640 
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)
644  ELSE
645  READ (iunit, *, iostat=istat(14)) &
646  (fsqt(i), wdot(i), i=1,nstore_seq)
647  END IF
648 
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)
656  ELSE
657  istat(14) = 0
658  END IF
659 
660  chipf = iotaf*phipf
661 !
662 ! CONVERT FROM INTERNAL UNITS TO PHYSICAL UNITS IF NEEDED
663 !
664  IF (version_ .le. (6.05+eps_w)) THEN
665  mass = mass/mu0
666  pres = pres/mu0
667  jcuru = jcuru/mu0
668  jcurv = jcurv/mu0
669  jdotb = jdotb/mu0
670  phi = -phi
671  END IF
672 
673 !-----------------------------------------------
674 ! DATA AND MSE FITS
675 !-----------------------------------------------
676  IF (ireconstruct .gt. 0) THEN
677 
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), &
689  i=1,isnodes)
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)
696  END IF
697 
698  IF (nobd .gt. 0) THEN
699  READ (iunit, *) (dsiext(i),plflux(i),dsiobt(i),i=1,nobd)
700  READ (iunit, *) flmwgt
701  END IF
702 
703  nbfldn = sum(nbfld(:nbsets))
704  IF (nbfldn .gt. 0) THEN
705  DO n = 1, nbsets
706  READ (iunit, *) (bcoil(i,n),plbfld(i,n),bbc(i,n), &
707  i=1,nbfld(n))
708  END DO
709  READ (iunit, *) bcwgt
710  END IF
711 
712  READ (iunit, *) phidiam, delphid
713 !
714 ! READ Limiter & Prout plotting specs
715 !
716  READ (iunit, *) nsets, nparts_in, nlim
717 
718  ALLOCATE (nsetsn(nsets))
719  READ (iunit, *) (nsetsn(i),i=1,nsets)
720 
721  n1 = maxval(nsetsn(:nsets))
722  ALLOCATE (pfcspec(nparts_in,n1,nsets), limitr(nlim))
723 
724  READ (iunit, *) (((pfcspec(i,j,k),i=1,nparts_in), &
725  j=1,nsetsn(k)),k=1,nsets)
726 
727  READ (iunit, *) (limitr(i), i=1,nlim)
728 
729  m = maxval(limitr(:nlim))
730  ALLOCATE (rlim(m,nlim), zlim(m,nlim))
731 
732  READ (iunit, *) ((rlim(i,j),zlim(i,j),i=1,limitr(j)), &
733  j=1,nlim)
734  READ (iunit, *) nrgrid, nzgrid
735  READ (iunit, *) tokid
736  READ (iunit, *) rx1, rx2, zy1, zy2, condif
737  READ (iunit, *) imatch_phiedge
738 
739  END IF
740 
741  1000 CONTINUE
742 
743  READ (iunit, iostat=ierr) mgrid_mode
744  IF (ierr .ne. 0) THEN
745  ierr = 0; mgrid_mode = 'N'
746  END IF
747 
748  IF (istat(2) .ne. 0) ierr_vmec = 1
749 
750  DO m = 1,15
751  IF (istat(m) .gt. 0) THEN
752  print *,' Error No. ',m,' in READ_WOUT, iostat = ',istat(m)
753  ierr = m
754  EXIT
755  END IF
756  END DO
757 
758 
759  720 FORMAT(8i10)
760  730 FORMAT(5e20.13)
761  740 FORMAT(a)
762  790 FORMAT(i5,/,(1p,3e12.4))
763 
764  END SUBROUTINE read_wout_text
765 
766 
767 #if defined(NETCDF)
768  SUBROUTINE read_wout_nc(filename, ierr)
769  USE ezcdf
770  USE stel_constants, ONLY: mu0
771  IMPLICIT NONE
772 !------------------------------------------------
773 ! D u m m y A r g u m e n t s
774 !------------------------------------------------
775  INTEGER, INTENT(out) :: ierr
776  CHARACTER(LEN=*), INTENT(in) :: filename
777 !------------------------------------------------
778 ! L o c a l V a r i a b l e s
779 !------------------------------------------------
780  INTEGER :: nwout, ierror
781  INTEGER, DIMENSION(3) :: dimlens
782 ! REAL(rprec) :: ohs
783  REAL(rprec), DIMENSION(:), ALLOCATABLE :: raxis_cc, raxis_cs, &
784  zaxis_cs, zaxis_cc
785 !------------------------------------------------
786 ! Open cdf File
787  CALL cdf_open(nwout,filename,'r', ierr)
788  IF (ierr .ne. 0) THEN
789  print *,' Error opening wout .nc file'
790  RETURN
791  END IF
792 
793 ! Be sure all arrays are deallocated
794  CALL read_wout_deallocate
795 
796 ! Read in scalar variables
797  CALL cdf_read(nwout, vn_error, ierr_vmec)
798 
799  IF (ierr_vmec.ne.norm_term_flag .and. &
800  ierr_vmec.ne.more_iter_flag) GOTO 1000
801 
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)
816  mnmaxpot=0
817  CALL cdf_read(nwout, vn_maxpot, mnmaxpot)
818  mnmax_nyq = -1
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)
823  IF (lasym) iasym = 1
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)
828 
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)
852  imse = -1
853  IF (lrecon) THEN
854  CALL cdf_read(nwout, vn_mse, imse)
855  CALL cdf_read(nwout, vn_thom, itse)
856  END IF
857  CALL cdf_read(nwout, vn_nextcur, nextcur)
858 
859  mgrid_mode = 'N'
860  CALL cdf_inquire(nwout, vn_mgmode, dimlens, ier=ierror)
861  IF (ierror.eq.0) CALL cdf_read(nwout, vn_mgmode, mgrid_mode)
862  IF (lfreeb) THEN
863  CALL cdf_read(nwout, vn_flp, nobser)
864  CALL cdf_read(nwout, vn_nobd, nobd)
865  CALL cdf_read(nwout, vn_nbset, nbsets)
866  END IF
867 
868 ! Inquire existence, dimensions of arrays for allocation
869 ! 1D Arrays
870  IF (lfreeb .and. nbsets.gt.0) THEN
871  CALL cdf_read(nwout, vn_nbfld, nbfld)
872  END IF
873 
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)
883  END IF
884 
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)
890  END IF
891 
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)
896  IF (lasym) THEN
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)
901  END IF
902 
903 ! Profile coefficients, dimensioned from 0
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)
910 
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)
923 
924  CALL cdf_inquire(nwout, vn_iotaf, dimlens)
925  ALLOCATE (iotaf(dimlens(1)), stat = ierror)
926  CALL cdf_inquire(nwout, vn_qfact, dimlens)
927 
928  ALLOCATE (qfact(dimlens(1)), stat = ierror)
929 
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)
935 
936  ALLOCATE (chi(dimlens(1)), stat = ierror)
937 
938  CALL cdf_inquire(nwout, vn_phipf, dimlens)
939  ALLOCATE (phipf(dimlens(1)), stat = ierror)
940 !OLD VERSION MAY NOT HAVE THIS! CALL cdf_inquire(nwout, vn_chipf, dimlens)
941  ALLOCATE (chipf(dimlens(1)), stat = ierror)
942  chipf = 1.e30_dp
943 
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)
968 
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)
975 
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)
988 
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)
993 
994  IF (nextcur .gt. 0) THEN
995  CALL cdf_inquire(nwout, vn_extcur, dimlens)
996  ALLOCATE (extcur(dimlens(1)), stat = ierror)
997 !NOTE: curlabel is an array of CHARACTER(30) strings - defined in mgrid_mod
998 ! so dimlens(1) == 30 (check this) and dimlens(2) is the number of strings in the array
999  CALL cdf_inquire(nwout, vn_curlab, dimlens)
1000  ALLOCATE (curlabel(dimlens(2)), stat = ierror)
1001  ENDIF
1002 
1003 ! 2D Arrays
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)
1020 
1021 ! ELIMINATE THESE EVENTUALLY: DON'T NEED THEM
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)
1026 
1027 ! The curr*mn* arrays have the same dimensions as the bsu**mn* arrays. No need
1028 ! to inquire about the dimension sizes.
1029  ALLOCATE (currumnc(dimlens(1),dimlens(2)), stat = ierror)
1030  ALLOCATE (currvmnc(dimlens(1),dimlens(2)), stat = ierror)
1031 
1032  IF (lfreeb) THEN
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)
1041  END IF
1042 
1043  IF (.NOT. lasym) GO TO 800
1044 
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)
1061 
1062 ! ELIMINATE THESE EVENTUALLY: DO NOT NEED THEM
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)
1067 
1068 ! The curr*mn* arrays have the same dimensions as the bsu**mn* arrays. No need
1069 ! to inquire about the dimension sizes.
1070  ALLOCATE (currumns(dimlens(1),dimlens(2)), stat = ierror)
1071  ALLOCATE (currvmns(dimlens(1),dimlens(2)), stat = ierror)
1072 
1073  IF (lfreeb) THEN
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)
1082  END IF
1083 
1084  800 CONTINUE
1085 
1086 ! Read Arrays
1087  CALL cdf_read(nwout, vn_pmod, xm)
1088  CALL cdf_read(nwout, vn_tmod, xn)
1089  IF (mnmax_nyq .le. 0) THEN
1090  mnmax_nyq = mnmax
1091  ALLOCATE (xm_nyq(mnmax_nyq), xn_nyq(mnmax_nyq), stat=ierror)
1092  xm_nyq = xm; xn_nyq = xn
1093  ELSE
1094  CALL cdf_read(nwout, vn_pmod_nyq, xm_nyq)
1095  CALL cdf_read(nwout, vn_tmod_nyq, xn_nyq)
1096  END IF
1097 
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) !includes nfp factor
1102  IF (lasym) CALL cdf_read(nwout, vn_potcos, potcos)
1103  END IF
1104 
1105  mnyq = int(maxval(xm_nyq)); nnyq = int(maxval(abs(xn_nyq)))/nfp
1106 
1107  CALL cdf_read(nwout, vn_racc, raxis_cc)
1108  CALL cdf_read(nwout, vn_zacs, zaxis_cs)
1109 
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)
1116 
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) !Half mesh
1121  CALL cdf_read(nwout, vn_bmnc, bmnc) !Half mesh
1122  CALL cdf_read(nwout, vn_bsubumnc, bsubumnc) !Half mesh
1123  CALL cdf_read(nwout, vn_bsubvmnc, bsubvmnc) !Half mesh
1124  CALL cdf_read(nwout, vn_bsubsmns, bsubsmns) !Full mesh
1125 ! ELIMINATE THESE EVENTUALLY: DON'T NEED THEM (can express in terms of lambdas)
1126  CALL cdf_read(nwout, vn_bsupumnc, bsupumnc)
1127  CALL cdf_read(nwout, vn_bsupvmnc, bsupvmnc)
1128 
1129  IF (version_ .ge. 9.0) THEN
1130  CALL cdf_read(nwout, vn_currumnc, currumnc)
1131  CALL cdf_read(nwout, vn_currvmnc, currvmnc)
1132  END IF
1133 
1134  IF (lfreeb) THEN
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)
1139  END IF
1140 
1141  IF (lasym) THEN
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)
1154 ! ELIMINATE THESE EVENTUALLY: DON'T NEED THEM
1155  CALL cdf_read(nwout, vn_bsupumns, bsupumns)
1156  CALL cdf_read(nwout, vn_bsupvmns, bsupvmns)
1157 
1158  IF (version_ .ge. 9.0) THEN
1159  CALL cdf_read(nwout, vn_currumns, currumns)
1160  CALL cdf_read(nwout, vn_currvmns, currvmns)
1161  END IF
1162 
1163  IF (lfreeb) THEN
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)
1168  END IF
1169  END IF
1170 
1171  CALL cdf_read(nwout, vn_am, am)
1172  CALL cdf_read(nwout, vn_ac, ac)
1173  CALL cdf_read(nwout, vn_ai, ai)
1174 
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)
1181 
1182  CALL cdf_read(nwout, vn_iotaf, iotaf)
1183  CALL cdf_read(nwout, vn_qfact, qfact)
1184 
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)
1189 
1190  CALL cdf_read(nwout, vn_chipf, chipf)
1191  IF (all(chipf .EQ. 1.e30_dp)) THEN
1192  chipf = iotaf*phipf
1193  END IF
1194 
1195  CALL cdf_read(nwout, vn_jcuru, jcuru)
1196  CALL cdf_read(nwout, vn_jcurv, jcurv)
1197 
1198 ! HALF-MESH quantities
1199 ! NOTE: jdotb is in units_of_A (1/mu0 incorporated in jxbforce...)
1200 ! prior to version 6.00, this was output in internal VMEC units...
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)
1213 
1214 ! MERCIER_CRITERION
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)
1221 
1222  CALL cdf_read(nwout, vn_fsq, fsqt)
1223  CALL cdf_read(nwout, vn_wdot, wdot)
1224 
1225  IF (nextcur .gt. 0) THEN
1226  CALL cdf_read(nwout, vn_extcur, extcur)
1227  CALL cdf_read(nwout, vn_curlab, curlabel)
1228  ENDIF
1229 
1230  1000 CONTINUE
1231 
1232  CALL cdf_close(nwout, ierr)
1233 
1234  IF (.not.ALLOCATED(bsubumnc)) RETURN !Moved this here because ns may not be set. SAL -09/07/11
1235 !
1236 ! COMPUTE CONTRAVARIANT CURRENT COMPONENTS IN AMPS
1237 ! ON THE FULL RADIAL MESH, WHERE JACOBIAN = SQRT(G)
1238 !
1239 ! CURRU = SQRT(G) * J dot grad(u)
1240 ! CURRV = SQRT(G) * J dot grad(v)
1241 !
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)
1250  END IF
1251  END IF
1252 
1253  IF (ierr .ne. 0) print *,'in read_wout_nc ierr=',ierr
1254  IF (ierror .ne. 0) print *,'in read_wout_nc ierror=',ierror
1255 
1256  END SUBROUTINE read_wout_nc
1257 #endif
1258 
1259  SUBROUTINE write_wout_text(filename, ierr)
1260  USE v3_utilities
1261  USE vsvd0, ONLY: nparts
1262  USE safe_open_mod
1263  USE stel_constants, ONLY: mu0
1264 
1265  IMPLICIT NONE
1266 !------------------------------------------------
1267 ! D u m m y A r g u m e n t s
1268 !------------------------------------------------
1269  CHARACTER (len=*) :: filename
1270  INTEGER, INTENT(out) :: ierr
1271 !------------------------------------------------
1272 ! L o c a l P a r a m e t e r s
1273 !------------------------------------------------
1274  REAL(rprec), PARAMETER :: eps_w = 1.e-4_dp
1275 !------------------------------------------------
1276 ! L o c a l V a r i a b l e s
1277 !------------------------------------------------
1278  INTEGER :: iounit, js, mn, i, j, k, m, n, iasymm
1279  LOGICAL :: lcurr
1280 !------------------------------------------------
1281 !
1282 ! THIS SUBROUTINE WRITES A TEXT FILE WOUT CREATED BY STORED THE INFORMATION
1283 ! IN THE read_WOUT MODULE. This routine can only be called if the wout has
1284 ! already been read in.
1285 
1286  iounit = 0
1287  ierr = 0
1288  CALL safe_open(iounit, ierr, &
1289  & 'wout_' // trim(filename) // '.txt', &
1290  & 'replace', 'formatted')
1291 
1292  CALL assert_eq(0, ierr, 'Error opening text wout file in ' // &
1293  & 'write_wout_text of read_wout_mod.')
1294 
1295 
1296 ! Write version info
1297  WRITE (iounit, '(a15,f4.2)') 'VMEC VERSION = ', version_
1298 
1299 ! Check version numbers since values change.
1300  IF (lasym) THEN
1301  iasymm = 1
1302  ELSE
1303  iasym = 0
1304  END IF
1305 
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
1309  ELSE
1310  IF (version_ .lt. 6.54) THEN
1311  WRITE (iounit, *) wb, wp, gamma, pfac, rmax_surf, rmin_surf
1312  ELSE
1313  WRITE (iounit, *) wb, wp, gamma, pfac, rmax_surf, rmin_surf, &
1314  & zmax_surf
1315  END IF
1316  IF (version_ .le. (8.0 + eps_w)) THEN
1317  WRITE (iounit, *) nfp, ns, mpol, ntor, mnmax, itfsq, niter, &
1318  & iasym, ireconstruct, ierr_vmec
1319  ELSE
1320  WRITE (iounit, *) nfp, ns, mpol, ntor, mnmax, mnmax_nyq, &
1321  & itfsq, niter, iasym, ireconstruct, &
1322  & ierr_vmec
1323  END IF
1324  END IF
1325 
1326  IF (version_ .gt. (6.20 + eps_w)) THEN
1327  WRITE (iounit, *) imse, itse, nbsets, nobd, nextcur, nstore_seq
1328  ELSE
1329  WRITE (iounit, *) imse, itse, nbsets, nobd, nextcur
1330  END IF
1331 
1332  IF (ierr_vmec .ne. norm_term_flag .and. &
1333  & ierr_vmec .ne. more_iter_flag) THEN
1334  GOTO 1000
1335  END IF
1336 
1337  IF (nbsets .gt. 0) THEN
1338  WRITE (iounit, *) nbfld(1:nbsets)
1339  END IF
1340  WRITE (iounit, *) trim(mgrid_file)
1341 
1342  DO js = 1, ns
1343  DO mn = 1, mnmax
1344  IF (js .eq. 1) THEN
1345  WRITE (iounit, *) nint(xm(mn)), nint(xn(mn)/nfp)
1346  END IF
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), &
1353  & currvmnc(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)
1360  ELSE
1361  WRITE (iounit, *) rmnc(mn,js), zmns(mn,js), lmns(mn,js)
1362  END IF
1363 
1364 ! Write asymmetric components.
1365  IF (lasym) THEN
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)
1372  ELSE
1373  WRITE (iounit, *) rmns(mn,js), zmnc(mn,js), &
1374  & lmnc(mn,js)
1375  END IF
1376  END IF
1377  END DO
1378 
1379  IF (version_ .le. (8.0 + eps_w)) THEN
1380  cycle
1381  END IF
1382 
1383  DO mn = 1, mnmax_nyq
1384  IF (js .eq. 1) THEN
1385  WRITE (iounit, *) nint(xm_nyq(mn)), &
1386  & nint(xn_nyq(mn)/nfp)
1387  END IF
1388  WRITE (iounit, *) bmnc(mn,js), gmnc(mn,js), &
1389  & bsubumnc(mn,js), bsubvmnc(mn,js), &
1390  & bsubsmns(mn,js), bsupumnc(mn,js), &
1391  & bsupvmnc(mn,js)
1392  IF (lasym) THEN
1393  WRITE (iounit, *) bmns(mn,js), gmns(mn,js), &
1394  & bsubumns(mn,js), bsubvmns(mn,js), &
1395  & bsubsmnc(mn,js), bsupumns(mn,js), &
1396  & bsupvmns(mn,js)
1397  END IF
1398  END DO
1399  END DO
1400 
1401 !
1402 ! Write FULL AND HALF-MESH QUANTITIES
1403 !
1404 ! NOTE: In version_ <= 6.00, mass, press were written out in INTERNAL (VMEC) units
1405 ! and are therefore multiplied here by 1/mu0 to transform to pascals. Same is true
1406 ! for ALL the currents (jcuru, jcurv, jdotb). Also, in version_ = 6.10 and
1407 ! above, PHI is the true (physical) toroidal flux (has the sign of jacobian correctly
1408 ! built into it)
1409 !
1410 
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), &
1422  & js=2, ns)
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), &
1429  & js=2, ns)
1430  WRITE (iounit, *) aspect, betatot, betapol, betaxis, b0
1431  ELSE
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), &
1437  & js = 2, ns)
1438  WRITE (iounit, *) aspect, betatot, betapol, betaxis, b0
1439  END IF
1440 
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
1446  END IF
1447 
1448 !-----------------------------------------------
1449 ! MERCIER CRITERION
1450 !-----------------------------------------------
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), &
1455  & js=2, ns - 1)
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), &
1459  & js=2, ns - 1)
1460  END IF
1461 
1462  IF (nextcur .gt. 0) THEN
1463  IF (version_ .le. (6.20 + eps_w)) THEN
1464  WRITE (iounit, 730) (extcur(js), js=1, nextcur)
1465  ELSE
1466  WRITE (iounit, *) (extcur(js), js=1, nextcur)
1467  END IF
1468 
1469  lcurr = len_trim(curlabel(1)) .gt. 0
1470  WRITE (iounit, *) lcurr
1471  IF (lcurr) THEN
1472  WRITE (iounit, *) (trim(curlabel(js)), js=1, nextcur)
1473  END IF
1474  END IF
1475 
1476  IF (version_ .le. (6.20 + eps_w)) THEN
1477  WRITE (iounit, 730) (fsqt(js), wdot(js), js = 1, nstore_seq)
1478  ELSE
1479  WRITE (iounit, *) (fsqt(js), wdot(js), js = 1, nstore_seq)
1480  END IF
1481 
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)
1487  END IF
1488 
1489 !-----------------------------------------------
1490 ! DATA AND MSE FITS
1491 !-----------------------------------------------
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), &
1503  & js=1, imse)
1504  WRITE (iounit, *) (rthom(js), datathom(i), js=1, itse)
1505  END IF
1506 
1507  IF (nobd .gt. 0) THEN
1508  WRITE (iounit, *) (dsiext(js), plflux(js), dsiobt(js), &
1509  & js=1, nobd)
1510  WRITE (iounit, *) flmwgt
1511  END IF
1512 
1513  IF (nbfldn .gt. 0) THEN
1514  DO n = 1, nbsets
1515  READ (iounit, *) (bcoil(i,n), plbfld(i,n), bbc(i,n), &
1516  & i=1,nbfld(n))
1517  END DO
1518  WRITE (iounit, *) bcwgt
1519  END IF
1520 
1521  WRITE (iounit, *) phidiam, delphid
1522 !
1523 ! READ Limiter & Prout plotting specs
1524 !
1525  WRITE (iounit, *) nsets, nparts, nlim
1526 
1527  WRITE (iounit, *) (nsetsn(js), js=1, nsets)
1528 
1529  WRITE (iounit, *) (((pfcspec(i,j,k), i=1, nparts), &
1530  & j=1, nsetsn(k)), k=1, nsets)
1531 
1532  WRITE (iounit, *) (limitr(i), i=1, nlim)
1533 
1534  WRITE (iounit, *) ((rlim(i,j), zlim(i,j), i=1, limitr(j)), &
1535  & j=1, nlim)
1536  WRITE (iounit, *) nrgrid, nzgrid
1537  WRITE (iounit, *) tokid
1538  WRITE (iounit, *) rx1, rx2, zy1, zy2, condif
1539  WRITE (iounit, *) imatch_phiedge
1540  END IF
1541 
1542 1000 CONTINUE
1543 
1544  WRITE (iounit, *) mgrid_mode
1545 
1546  730 FORMAT(5e20.13)
1547 
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.')
1551 
1552  END SUBROUTINE
1553 
1554  SUBROUTINE compute_currents(bsubsmnc_, bsubsmns_, &
1555  & bsubumnc_, bsubumns_, &
1556  & bsubvmnc_, bsubvmns_, &
1557  & xm_nyq_, xn_nyq_, mnmax_nyq_, &
1558  & lasym_, ns_, &
1559  & currumnc_, currvmnc_, &
1560  & currumns_, currvmns_)
1561  USE stel_constants, ONLY: mu0
1562  IMPLICIT NONE
1563 
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_
1570 
1571  REAL(rprec), DIMENSION(:), INTENT(in) :: xm_nyq_
1572  REAL(rprec), DIMENSION(:), INTENT(in) :: xn_nyq_
1573 
1574  INTEGER, INTENT(in) :: mnmax_nyq_
1575  LOGICAL, INTENT(in) :: lasym_
1576  INTEGER, INTENT(in) :: ns_
1577 
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_
1582 
1583 !-----------------------------------------------
1584 ! L o c a l V a r i a b l e s
1585 !-----------------------------------------------
1586  INTEGER :: js
1587  REAL(rprec) :: ohs, hs, shalf(ns_), sfull(ns_)
1588  REAL(rprec), DIMENSION(mnmax_nyq_) :: bu1, bu0, bv1, bv0, t1, t2, t3
1589 !-----------------------------------------------
1590 !
1591 ! Computes current harmonics for currXmn == sqrt(g)*JsupX, X = u,v
1592 ! [Corrected above "JsubX" to "JsupX", JDH 2010-08-16]
1593 
1594 ! NOTE: bsub(s,u,v)mn are on HALF radial grid
1595 ! (in earlier versions, bsubsmn was on FULL radial grid)
1596 
1597 ! NOTE: near the axis, b_s is dominated by the m=1 component of gsv ~ cos(u)/sqrt(s)
1598 ! we average it with a weight-factor of sqrt(s)
1599 !
1600 
1601  ohs = (ns_-1)
1602  hs = 1._dp/ohs
1603 
1604  DO js = 2, ns_
1605  shalf(js) = sqrt(hs*(js-1.5_dp))
1606  sfull(js) = sqrt(hs*(js-1))
1607  END DO
1608 
1609  DO js = 2, ns_-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)
1619  ELSEWHERE
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))
1623  ENDWHERE
1624  currumnc_(:,js) = -xn_nyq_(:)*t1 - t3
1625  currvmnc_(:,js) = -xm_nyq_(:)*t1 + t2
1626  END DO
1627 
1628  WHERE (xm_nyq_ .LE. 1)
1629  currvmnc_(:,1) = 2*currvmnc_(:,2) - currvmnc_(:,3)
1630  currumnc_(:,1) = 2*currumnc_(:,2) - currumnc_(:,3)
1631  ELSEWHERE
1632  currvmnc_(:,1) = 0
1633  currumnc_(:,1) = 0
1634  ENDWHERE
1635 
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
1639 
1640  IF (.NOT.lasym_) RETURN
1641 
1642  DO js = 2, ns_-1
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)
1652  ELSEWHERE
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))
1656  END WHERE
1657  currumns_(:,js) = xn_nyq_(:)*t1 - t3
1658  currvmns_(:,js) = xm_nyq_(:)*t1 + t2
1659  END DO
1660 
1661  WHERE (xm_nyq_ .LE. 1)
1662  currvmns_(:,1) = 2*currvmns_(:,2) - currvmns_(:,3)
1663  currumns_(:,1) = 2*currumns_(:,2) - currumns_(:,3)
1664  ELSEWHERE
1665  currvmns_(:,1) = 0
1666  currumns_(:,1) = 0
1667  END WHERE
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
1671 
1672  END SUBROUTINE compute_currents
1673 
1674  SUBROUTINE read_wout_deallocate
1675  IMPLICIT NONE
1676 !-----------------------------------------------
1677 ! L o c a l V a r i a b l e s
1678 !-----------------------------------------------
1679  INTEGER :: istat(10)
1680 !-----------------------------------------------
1681  istat=0
1682  lwout_opened=.false.
1683 
1684  IF (ALLOCATED(extcur)) DEALLOCATE (extcur, curlabel, &
1685  stat = istat(1))
1686  IF (ALLOCATED(overr)) DEALLOCATE (overr, stat = istat(2))
1687 
1688  IF (ALLOCATED(xm)) DEALLOCATE (xm, xn, xm_nyq, xn_nyq, &
1689  rmnc, zmns, lmns, bmnc, gmnc, bsubumnc, iotaf, presf, phipf, &
1690  chipf, &
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))
1695 
1696  IF (ALLOCATED(potsin)) DEALLOCATE (potsin)
1697  IF (ALLOCATED(potcos)) DEALLOCATE (potcos)
1698 
1699  IF (ALLOCATED(chipf)) DEALLOCATE (chipf, chi)
1700 
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))
1703 
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))
1709 
1710  IF (ALLOCATED(rmns)) DEALLOCATE (rmns, zmnc, lmnc, &
1711  bmns, gmns, bsubumns, bsubvmns, bsubsmnc, &
1712  bsupumns, bsupvmns, stat=istat(5))
1713 
1714  IF (ALLOCATED(bsubumnc_sur)) THEN
1715  DEALLOCATE(bsubumnc_sur, bsubvmnc_sur, &
1716  bsupumnc_sur, bsupvmnc_sur)
1717  END IF
1718 
1719  IF (ALLOCATED(bsubumns_sur)) THEN
1720  DEALLOCATE(bsubumns_sur, bsubvmns_sur, &
1721  bsupumns_sur, bsupvmns_sur)
1722  END IF
1723 
1724 ! Note currvmnc deallocated above.
1725  IF (ALLOCATED(currumnc)) DEALLOCATE (currumnc)
1726  IF (ALLOCATED(currumns)) DEALLOCATE (currumns, currvmns)
1727  IF (ALLOCATED(rzl_local)) DEALLOCATE (rzl_local)
1728 
1729  IF (any(istat .ne. 0)) &
1730  stop 'Deallocation error in read_wout_deallocate'
1731 
1732  END SUBROUTINE read_wout_deallocate
1733 
1734  SUBROUTINE tosuvspace (s_in, u_in, v_in, gsqrt, &
1735  bsupu, bsupv, jsupu, jsupv)
1736  USE stel_constants, ONLY: zero, one
1737  IMPLICIT NONE
1738 !------------------------------------------------
1739 ! D u m m y A r g u m e n t s
1740 !------------------------------------------------
1741  REAL(rprec), INTENT(in) :: s_in, u_in, v_in
1742  REAL(rprec), INTENT(out), OPTIONAL :: gsqrt, bsupu, bsupv, &
1743  jsupu, jsupv
1744 !------------------------------------------------
1745 ! L o c a l V a r i a b l e s
1746 !------------------------------------------------
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
1757 !------------------------------------------------
1758 !
1759 ! COMPUTE VARIOUS HALF/FULL-RADIAL GRID QUANTITIES AT THE INPUT POINT
1760 ! (S, U, V) , WHERE
1761 ! S = normalized toroidal flux (0 - 1),
1762 ! U = poloidal angle
1763 ! V = N*phi = toroidal angle * no. field periods
1764 !
1765 ! HALF-RADIAL GRID QUANTITIES
1766 ! gsqrt, bsupu, bsupv
1767 !
1768 ! FULL-RADIAL GRID QUANTITIES
1769 ! dbsubuds, dbsubvds, dbsubsdu, dbsubsdv
1770 !
1771 !------------------------------------------------
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'
1774  RETURN
1775  END IF
1776 
1777  IF (.not.lwout_opened) THEN
1778  WRITE(6, *)' tosuvspace can only be called AFTER opening wout file!'
1779  RETURN
1780  END IF
1781 
1782 !
1783 ! SETUP TRIG ARRAYS
1784 !
1785  cosu = cos(u_in); sinu = sin(u_in)
1786  cosv = cos(v_in); sinv = sin(v_in)
1787 
1788  cosmu(0) = 1; sinmu(0) = 0
1789  cosnv(0) = 1; sinnv(0) = 0
1790  DO m = 1, mnyq
1791  cosmu(m) = cosmu(m-1)*cosu - sinmu(m-1)*sinu
1792  sinmu(m) = sinmu(m-1)*cosu + cosmu(m-1)*sinu
1793  END DO
1794 
1795  DO n = 1, nnyq
1796  cosnv(n) = cosnv(n-1)*cosv - sinnv(n-1)*sinv
1797  sinnv(n) = sinnv(n-1)*cosv + cosnv(n-1)*sinv
1798  END DO
1799 
1800 
1801 !
1802 ! FIND INTERPOLATED s VALUE AND COMPUTE INTERPOLATION WEIGHTS wlo, whi
1803 ! RECALL THAT THESE QUANTITIES ARE ON THE HALF-RADIAL GRID...
1804 ! s-half(j) = (j-1.5)*hs, for j = 2,...ns
1805 !
1806  hs1 = one/(ns-1)
1807  jslo = int(c1p5 + s_in/hs1)
1808  jshi = jslo+1
1809  wlo = (hs1*(jshi-c1p5) - s_in)/hs1
1810  whi = 1 - wlo
1811  IF (jslo .eq. ns) THEN
1812 ! USE Xhalf(ns+1) = 2*Xhalf(ns) - Xhalf(ns-1) FOR "GHOST" POINT VALUE 1/2hs OUTSIDE EDGE
1813 ! THEN, X = wlo*Xhalf(ns) + whi*Xhalf(ns+1) == Xhalf(ns) + whi*(Xhalf(ns) - Xhalf(ns-1))
1814  jshi = jslo-1
1815  wlo = 1+whi; whi = -whi
1816  ELSE IF (jslo .eq. 1) THEN
1817  jslo = 2
1818  END IF
1819 
1820 !
1821 ! FOR ODD-m MODES X ~ SQRT(s), SO INTERPOLATE Xmn/SQRT(s)
1822 !
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)))
1826  ELSE
1827  wlo_odd = 0
1828  whi_odd = sqrt(s_in/(hs1*(jshi-c1p5)))
1829  END IF
1830 
1831  WHERE (mod(nint(xm_nyq(:)),2) .eq. 0)
1832  wmins = wlo
1833  wplus = whi
1834  ELSEWHERE
1835  wmins = wlo_odd
1836  wplus = whi_odd
1837  END WHERE
1838 
1839  ipresent = 0
1840  lgsqrt = PRESENT(gsqrt)
1841  IF (lgsqrt) THEN
1842  gsqrt = 0 ; ipresent = ipresent+1
1843  gmnc1 = wmins*gmnc(:,jslo) + wplus*gmnc(:,jshi)
1844  IF (lasym) gmns1 = wmins*gmns(:,jslo) + wplus*gmns(:,jshi)
1845  END IF
1846  lbsupu = PRESENT(bsupu)
1847  IF (lbsupu) THEN
1848  bsupu = 0 ; ipresent = ipresent+1
1849  bsupumnc1 = wmins*bsupumnc(:,jslo) + wplus*bsupumnc(:,jshi)
1850  IF (lasym) bsupumns1 = wmins*bsupumns(:,jslo) + wplus*bsupumns(:,jshi)
1851  END IF
1852  lbsupv = PRESENT(bsupv)
1853  IF (lbsupv) THEN
1854  bsupv = 0 ; ipresent = ipresent+1
1855  bsupvmnc1 = wmins*bsupvmnc(:,jslo) + wplus*bsupvmnc(:,jshi)
1856  IF (lasym) bsupvmns1 = wmins*bsupvmns(:,jslo) + wplus*bsupvmns(:,jshi)
1857  END IF
1858 
1859  IF (ipresent .eq. 0) GOTO 1000
1860 
1861 !
1862 ! COMPUTE GSQRT, ... IN REAL SPACE
1863 ! tcosmn = cos(mu - nv); tsinmn = sin(mu - nv)
1864 !
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
1872  END DO
1873 
1874  IF (.not.lasym) GOTO 1000
1875 
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
1883  END DO
1884 
1885  1000 CONTINUE
1886 
1887 ! FULL-MESH QUANTITIES
1888 !
1889 ! FIND INTERPOLATED s VALUE AND COMPUTE INTERPOLATION WEIGHTS wlo, whi
1890 ! RECALL THAT THESE QUANTITIES ARE ON THE FULL-RADIAL GRID...
1891 ! s-full(j) = (j-1)*hs, for j = 1,...ns
1892 !
1893  hs1 = one/(ns-1)
1894  jslo = 1+int(s_in/hs1)
1895  jshi = jslo+1
1896  IF (jslo .eq. ns) jshi = ns
1897  wlo = (hs1*(jshi-1) - s_in)/hs1
1898  whi = 1 - wlo
1899 !
1900 ! FOR ODD-m MODES X ~ SQRT(s), SO INTERPOLATE Xmn/SQRT(s)
1901 !
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)))
1905  ELSE
1906  wlo_odd = 0
1907  whi_odd = sqrt(s_in/(hs1*(jshi-1)))
1908  END IF
1909 
1910  WHERE (mod(nint(xm_nyq(:)),2) .eq. 0)
1911  wmins = wlo
1912  wplus = whi
1913  ELSEWHERE
1914  wmins = wlo_odd
1915  wplus = whi_odd
1916  END WHERE
1917 
1918  ipresent = 0
1919  ljsupu = PRESENT(jsupu)
1920  IF (ljsupu) THEN
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(:,jshi)
1925  END IF
1926 
1927  ljsupv = PRESENT(jsupv)
1928  IF (ljsupv) THEN
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(:,jshi)
1933  END IF
1934 
1935  IF (ipresent .eq. 0) RETURN
1936 
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
1943  END DO
1944 
1945  IF (.not.lasym) GOTO 2000
1946 
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
1953  END DO
1954 
1955  2000 CONTINUE
1956 
1957  IF (ljsupu) jsupu = jsupu/gsqrt
1958  IF (ljsupv) jsupv = jsupv/gsqrt
1959 
1960  END SUBROUTINE tosuvspace
1961 
1962  SUBROUTINE loadrzl
1963  IMPLICIT NONE
1964 !------------------------------------------------
1965 ! L o c a l V a r i a b l e s
1966 !------------------------------------------------
1967  INTEGER :: rcc, rss, zsc, zcs, rsc, rcs, zcc, zss
1968  INTEGER :: mpol1, mn, m, n, n1
1969  REAL(rprec) :: sgn
1970 !------------------------------------------------
1971 !
1972 ! Arrays must be stacked (and ns,ntor,mpol ordering imposed)
1973 ! as coefficients of cos(mu)*cos(nv), etc
1974 ! Only need R, Z components(not lambda, for now anyhow)
1975 !
1976  IF (ALLOCATED(rzl_local)) RETURN
1977 
1978  mpol1 = mpol-1
1979  rcc = 1; zsc = 1
1980  IF (.not.lasym) THEN
1981  IF (lthreed) THEN
1982  ntmax = 2
1983  rss = 2; zcs = 2
1984  ELSE
1985  ntmax = 1
1986  END IF
1987  ELSE
1988  IF (lthreed) THEN
1989  ntmax = 4
1990  rss = 2; rsc = 3; rcs = 4
1991  zcs = 2; zcc = 3; zss = 4
1992  ELSE
1993  ntmax = 2
1994  rsc = 2; zcc = 2
1995  END IF
1996  END IF
1997 
1998 ! really only need to ALLOCATE 2*ntmax (don't need lambdas)
1999 ! for consistency, we'll allocate 3*ntmax and set lambdas = 0
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'
2003  rzl_local = 0
2004 
2005  DO mn = 1, mnmax
2006  m = nint(xm(mn)); n = nint(xn(mn))/nfp; n1 = abs(n)
2007  sgn = sign(1, 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,:)
2010  IF (lthreed) THEN
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,:)
2013  END IF
2014  IF (lasym) THEN
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,:)
2017  IF (lthreed) THEN
2018  rzl_local(:,n1,m,rcs) = rzl_local(:,n1,m,rcs) &
2019  - sgn*rmns(mn,:)
2020  rzl_local(:,n1,m,zss) = rzl_local(:,n1,m,zss) &
2021  + sgn*zmnc(mn,:)
2022  END IF
2023  END IF
2024  END DO
2025 
2026  END SUBROUTINE loadrzl
2027 
2028  END MODULE read_wout_mod
2029 
v3_utilities::assert_eq
Definition: v3_utilities.f:62
Generated on Thu Mar 5 2020 15:49:23 for V3FIT by   doxygen 1.8.17