V3FIT
wrout.f
1  SUBROUTINE wrout(bsq, gsqrt, bsubu, bsubv, bsubs, bsupv, bsupu,
2  1 rzl_array, gc_array, ier_flag, lwrite
3 #ifdef _ANIMEC
4  2 ,tau_an, sigma_an, ppar, pperp, onembc, pbprim,
5  3 ppprim, densit
6 #endif
7  4 )
8 ! ... from SPH 2009-10-05; changes for modB sine-harmonics included
9  USE vmec_input, ONLY: ns_array, ftol_array, lwouttxt
10  USE vmec_params
11  USE vmec_main
12  USE vmercier
13  USE vmec_persistent
14  USE vparams, p5 => cp5, two => c2p0
15  USE vac_persistent
16  USE vspline
17  USE xstuff
18  USE vmec_io
19  USE realspace, ONLY: phip, chip, gsqrta=>z1, z1=>z1
20  USE totzsp_mod
21  USE vforces, ONLY: bsupua=>brmn_e, bsupva=>czmn_o, bsqa=>bzmn_e,
22  1 bsubsa=>armn_e, bsubua=>azmn_e, bsubva=>armn_o
23 #ifdef _VACUUM2
24  USE vac2_vacmod, ONLY: potvac, mnpd, xmpot, xnpot
25 #else
26  USE vacmod, ONLY: potvac, mnpd, xmpot, xnpot !added for diagno, J.Geiger
27 #endif
28 #ifdef _HBANGLE
29  USE angle_constraints, ONLY: getrz
30 #endif
31 !undef NETCDF IF TXT DESIRED
32 #ifdef NETCDF
33  USE ezcdf
34  USE read_wout_mod, ONLY: compute_currents,
35  1 vn_version, vn_extension, vn_mgrid,
36  1 vn_magen, vn_therm, vn_gam, vn_maxr, vn_minr, vn_maxz, vn_fp,
37  2 vn_radnod, vn_polmod, vn_tormod, vn_maxmod, vn_maxit, vn_actit,
38  3 vn_asym, vn_recon, vn_free, vn_error, vn_aspect, vn_beta,
39  4 vn_pbeta, vn_tbeta, vn_abeta, vn_b0, vn_rbt0, vn_maxmod_nyq,
40  5 vn_rbt1, vn_sgs, vn_lar, vn_modb, vn_ctor, vn_amin, vn_rmaj,
41  5 vn_potsin, vn_potcos, vn_maxpot, vn_xmpot, vn_xnpot, !diagno/extender output (SPH071414)
42  6 vn_vol, vn_mse, vn_thom, vn_ac, vn_ai, vn_am, vn_rfp,
43  6 vn_pmass_type, vn_pcurr_type, vn_piota_type,
44  6 vn_am_aux_s, vn_am_aux_f, vn_ac_aux_s, vn_ac_aux_f,
45  6 vn_ai_aux_s, vn_ai_aux_f,
46  6 vn_ftolv, vn_fsqr, vn_fsqz, vn_fsql,
47  7 vn_pmod, vn_tmod, vn_pmod_nyq, vn_tmod_nyq,
48  7 vn_racc, vn_zacs, vn_racs, vn_zacc, vn_iotaf, vn_qfact,
49  8 vn_presf, vn_phi, vn_phipf, vn_jcuru, vn_jcurv, vn_iotah,
50  8 vn_chi, vn_chipf,
51  9 vn_mass, vn_presh, vn_betah, vn_buco, vn_bvco, vn_vp, vn_specw,
52  a vn_phip, vn_jdotb, vn_bdotb, vn_overr, vn_bgrv, vn_merc,
53  b vn_mshear, vn_mwell, vn_mcurr, vn_mgeo, vn_equif, vn_fsq,
54  c vn_wdot, vn_extcur, vn_curlab, vn_rmnc, vn_zmns, vn_lmns,
55  d vn_gmnc, vn_bmnc, vn_bsubumnc, vn_bsubvmnc, vn_bsubsmns,
56  e vn_bsupumnc, vn_bsupvmnc, vn_rmns, vn_zmnc, vn_lmnc, vn_gmns,
57  f vn_bmns, vn_bsubumns, vn_bsubvmns, vn_bsubsmnc, vn_bsupumns,
58  g vn_bsupvmns, vn_rbc, vn_zbs, vn_rbs, vn_zbc,
59  h ln_version, ln_extension, ln_mgrid,
60 
61  & vn_bsubumnc_sur, vn_bsubvmnc_sur, !MRC 10-15-15
62  & vn_bsupumnc_sur, vn_bsupvmnc_sur,
63  & vn_bsubumns_sur, vn_bsubvmns_sur,
64  & vn_bsupumns_sur, vn_bsupvmns_sur,
65  & vn_currumnc, vn_currumns, vn_currvmnc, vn_currvmns, !MRC 8-12-16
66 
67  1 ln_magen, ln_therm, ln_gam, ln_maxr, ln_minr, ln_maxz, ln_fp,
68  2 ln_radnod, ln_polmod, ln_tormod, ln_maxmod, ln_maxit, ln_actit,
69  2 ln_maxpot, ln_potsin, ln_potcos,
70  3 ln_asym, ln_recon, ln_free, ln_error, ln_aspect, ln_beta,
71  4 ln_pbeta, ln_tbeta, ln_abeta, ln_b0, ln_rbt0, ln_maxmod_nyq,
72  5 ln_rbt1, ln_sgs, ln_lar, ln_modb, ln_ctor, ln_amin, ln_rmaj,
73  6 ln_mse, ln_thom, ln_flp, ln_nobd, ln_nbset, ln_next, ln_nbfld,
74  7 ln_pmod, ln_tmod, ln_pmod_nyq, ln_tmod_nyq, ln_racc, ln_zacs,
75  7 ln_racs, ln_zacc, ln_iotaf, ln_qfact, ln_am, ln_ac, ln_ai,
76  7 ln_pmass_type, ln_pcurr_type, ln_piota_type,
77  7 ln_am_aux_s, ln_am_aux_f, ln_ac_aux_s, ln_ac_aux_f,
78  7 ln_ai_aux_s, ln_ai_aux_f, ln_chi, ln_chipf,
79  8 ln_presf, ln_phi, ln_phipf, ln_jcuru, ln_jcurv, ln_iotah,
80  9 ln_mass, ln_presh, ln_betah, ln_buco, ln_bvco, ln_vp, ln_specw,
81  a ln_vol, ln_phip, ln_jdotb, ln_bdotb, ln_bgrv, ln_merc,
82  b ln_mshear, ln_mwell, ln_mcurr, ln_mgeo, ln_equif, ln_fsq,
83  c ln_wdot, ln_extcur, ln_curlab, ln_rmnc, ln_zmns, ln_lmns,
84  d ln_gmnc, ln_bmnc, ln_bsubumnc, ln_bsubvmnc, ln_bsubsmns,
85  e ln_bsupumnc, ln_bsupvmnc, ln_rmns, ln_zmnc, ln_lmnc, ln_gmns,
86  f ln_bmns, ln_bsubumns, ln_bsubvmns, ln_bsubsmnc, ln_bsupumns,
87  g ln_bsupvmns, ln_rbc, ln_zbs, ln_rbs, ln_zbc,
88 
89  & ln_bsubumnc_sur, ln_bsubvmnc_sur, !MRC 10-15-15
90  & ln_bsupumnc_sur, ln_bsupvmnc_sur,
91  & ln_bsubumns_sur, ln_bsubvmns_sur,
92  & ln_bsupumns_sur, ln_bsupvmns_sur,
93  & ln_currumnc, ln_currumns, ln_currvmnc, ln_currvmns
94 
95 !------------------DEC$ ELSE !to use safe_open_mod in any case (J.Geiger)
96 #endif
97  USE safe_open_mod
98  USE mgrid_mod
99 
100  IMPLICIT NONE
101 !-----------------------------------------------
102 ! D u m m y A r g u m e n t s
103 !-----------------------------------------------
104  INTEGER, INTENT(in) :: ier_flag
105  REAL(dp), DIMENSION(mnmax,ns,3*MAX(ntmax/2,1)), !reverse ns, mnmax for backwards compatibility
106  1 INTENT(inout), TARGET :: rzl_array, gc_array
107  REAL(dp), DIMENSION(ns,nznt), INTENT(inout) ::
108  1 bsq, gsqrt, bsubu, bsubv, bsubs, bsupv, bsupu
109 #ifdef _ANIMEC
110  2 ,tau_an, ppar, pperp, onembc, sigma_an
111  REAL(dp), DIMENSION(ns,nznt), INTENT(out) ::
112  1 densit, pbprim, ppprim
113 #endif
114  REAL(dp) :: qfact(ns)
115  LOGICAL :: lwrite
116 !-----------------------------------------------
117 ! L o c a l P a r a m e t e r s
118 !-----------------------------------------------
119  REAL(dp), PARAMETER :: c1p5 = 1.5_dp
120  LOGICAL :: lnyquist = .true. !=false, suppress nyquist stuff
121 #ifdef NETCDF
122  CHARACTER(LEN=*), PARAMETER, DIMENSION(1) ::
123  1 r1dim = (/'radius'/), mn1dim = (/'mn_mode'/),
124  2 mn2dim = (/'mn_mode_nyq'/),
125  2 mnpotdim = (/'mn_mode_pot'/),
126  3 currg = (/'ext_current'/),
127  4 currl = (/'current_label'/)
128  CHARACTER(LEN=*), DIMENSION(2), PARAMETER ::
129  1 r2dim = (/'mn_mode','radius '/),
130  1 r3dim = (/'mn_mode_nyq','radius '/)
131 #endif
132 !-----------------------------------------------
133 ! L o c a l V a r i a b l e s
134 !-----------------------------------------------
135  INTEGER :: j, js, jlk, mn, lk, iasym,
136  1 m, n, k, iwout0, n1, nwout, istat, i, indx1(1),
137  2 mnmax_nyq0, mnyq0, nnyq0, nwout2 ! nwout2 by J.Geiger
138  3 ,isgn, js2, nfort !for diagno 1.5
139  REAL(dp) :: dmult, tcosi, tsini, vversion, sgn, tmult,
140  1 presfactor, ftolx1, d_bsupumn, d_bsupvmn ! diagno 1.5
141 #ifdef _ANIMEC
142  2 ,hotdam, omtbc, optbc, pdh, pmh, pde, pme, eps
143 #endif
144  REAL(dp), POINTER, DIMENSION(:,:) :: rmnc, rmns, zmns,
145  1 zmnc, lmns, lmnc
146  REAL(dp), ALLOCATABLE, DIMENSION(:,:) ::
147  1 gmnc, bmnc, gmns, bmns,
148  2 bsubumnc, bsubvmnc, bsubsmns, bsubumns, bsubvmns, bsubsmnc,
149  3 currumnc, currvmnc, currumns, currvmns
150 #ifdef _ANIMEC
151  3 ,sigmnc , taumnc , pparmnc , ppermnc , pbprmnc , ppprmnc ,
152  4 hotdmnc , hotdmns ,
153  5 sigmns , taumns , pparmns , ppermns , pbprmns , ppprmns
154  REAL(dp), DIMENSION(:,:), ALLOCATABLE :: sigma_ana, tau_ana,
155  1 ppara, pperpa, pbprima, ppprima, densita
156 #endif
157  REAL(dp), DIMENSION(mnmax) :: rmnc1, zmns1, lmns1,
158  1 rmns1, zmnc1, lmnc1, bmodmn, bmodmn1
159  REAL(dp), DIMENSION(:), ALLOCATABLE :: gmn, bmn,
160  1 bsubumn, bsubvmn, bsubsmn, bsupumn, bsupvmn
161 #ifdef _ANIMEC
162  2 ,sigmn , taumn , pparmn , ppermn , pbprmn , ppprmn , hotdmn
163 #endif
164  REAL(dp), DIMENSION(:), ALLOCATABLE :: bsubumnc_sur !MRC 10-15-15
165  REAL(dp), DIMENSION(:), ALLOCATABLE :: bsubvmnc_sur
166  REAL(dp), DIMENSION(:), ALLOCATABLE :: bsupumnc_sur
167  REAL(dp), DIMENSION(:), ALLOCATABLE :: bsupvmnc_sur
168  REAL(dp), DIMENSION(:), ALLOCATABLE :: bsubumns_sur
169  REAL(dp), DIMENSION(:), ALLOCATABLE :: bsubvmns_sur
170  REAL(dp), DIMENSION(:), ALLOCATABLE :: bsupumns_sur
171  REAL(dp), DIMENSION(:), ALLOCATABLE :: bsupvmns_sur
172  REAL(dp), DIMENSION(:), ALLOCATABLE :: bsubua_sur, bsubva_sur
173  REAL(dp), DIMENSION(:), ALLOCATABLE :: bsupua_sur, bsupva_sur
174 
175  CHARACTER(LEN=120) :: wout_file, wout2_file ! wout2_file by J.Geiger
176  CHARACTER(LEN=120) :: fort_file ! fort_file for diagno 1.5
177  REAL(dp), DIMENSION(:), ALLOCATABLE :: xfinal
178  REAL(dp), DIMENSION(:), POINTER :: xm_nyq0, xn_nyq0
179 ! ELIMINATE THESE EVENTUALLY
180  REAL(dp), ALLOCATABLE, DIMENSION(:,:) ::
181  1 bsupumnc, bsupumns, bsupvmnc, bsupvmns
182 
183  LOGICAL :: lcurr
184  INTEGER :: nmin0 ! J Geiger: Added for diagno-file
185 
186 !-----------------------------------------------
187  CALL second0 (twouton)
188 !
189 ! Pointer assignments for storage arrays
190 !
191  n1 = max(1,ntmax/2)
192  rmnc => rzl_array(:,:,1) !!store COS(mu-nv) components
193  zmns => rzl_array(:,:,1+n1) !!store SIN(mu-nv)
194  lmns => rzl_array(:,:,1+2*n1) !!store SIN(mu-nv)
195 
196  IF (lasym) THEN
197  rmns => gc_array(:,:,1) !!store SIN(mu-nv)
198  zmnc => gc_array(:,:,1+n1) !!store COS(mu-nv)
199  lmnc => gc_array(:,:,1+2*n1) !!store COS(mu-nv)
200  END IF
201 
202 !
203 ! THIS SUBROUTINE CREATES THE FILE WOUT.IT CONTAINS THE CYLINDRICAL COORDINATE SPECTRAL
204 ! COEFFICIENTS RMN,ZMN (full), LMN (half_mesh - CONVERTED FROM
205 ! INTERNAL full REPRESENTATION), AS WELL AS COEFFICIENTS (ON NYQ MESH) FOR COMPUTED
206 ! QUANTITIES:
207 !
208 ! BSQ, BSUPU,V, BSUBU,V, GSQRT (HALF); BSUBS (FULL-CONVERTED IN JXBFORCE)
209 !
210  IF (lnyquist) THEN
211  mnmax_nyq0 = mnmax_nyq
212  mnyq0 = mnyq
213  nnyq0 = nnyq
214  xm_nyq0 => xm_nyq; xn_nyq0 => xn_nyq
215  ELSE
216  mnmax_nyq0 = mnmax
217  mnyq0 = mpol1
218  nnyq0 = ntor
219  xm_nyq0 => xm; xn_nyq0 => xn
220  END IF
221 
222  ALLOCATE (gmn(mnmax_nyq0), bmn(mnmax_nyq0),
223  1 bsubumn(mnmax_nyq0), bsubvmn(mnmax_nyq0), bsubsmn(mnmax_nyq0),
224  2 bsupumn(mnmax_nyq0), bsupvmn(mnmax_nyq0),
225 #ifdef _ANIMEC
226  3 sigmn(mnmax_nyq0) ,
227  4 taumn(mnmax_nyq0) , pparmn(mnmax_nyq0) , ppermn(mnmax_nyq0) ,
228  5 pbprmn(mnmax_nyq0) , ppprmn(mnmax_nyq0) , hotdmn(mnmax_nyq0) ,
229 #endif
230  6 stat=istat)
231 
232  IF (lfreeb) THEN !MRC 10-15-15
233  ALLOCATE (bsubua_sur(nzeta*ntheta2), bsubva_sur(nzeta*ntheta2))
234  ALLOCATE (bsupua_sur(nzeta*ntheta2), bsupva_sur(nzeta*ntheta2))
235 
236  ALLOCATE (bsubumnc_sur(mnmax_nyq0), bsubvmnc_sur(mnmax_nyq0))
237  ALLOCATE (bsupumnc_sur(mnmax_nyq0), bsupvmnc_sur(mnmax_nyq0))
238  IF (lasym) THEN
239  ALLOCATE (bsubumns_sur(mnmax_nyq0), &
240  & bsubvmns_sur(mnmax_nyq0))
241  ALLOCATE (bsupumns_sur(mnmax_nyq0), &
242  & bsupvmns_sur(mnmax_nyq0))
243  END IF
244  END IF
245 
246  ALLOCATE (gmnc(mnmax_nyq0,ns), bmnc(mnmax_nyq0,ns),
247  1 bsubumnc(mnmax_nyq0,ns), bsubvmnc(mnmax_nyq0,ns),
248  2 bsubsmns(mnmax_nyq0,ns), bsupumnc(mnmax_nyq0,ns),
249  3 bsupvmnc(mnmax_nyq0,ns),
250  4 currumnc(mnmax_nyq0,ns), currvmnc(mnmax_nyq0,ns),
251 #ifdef _ANIMEC
252  5 sigmnc(mnmax_nyq0,ns) ,
253  6 taumnc(mnmax_nyq0,ns) , pparmnc(mnmax_nyq0,ns) ,
254  7 ppermnc(mnmax_nyq0,ns) , pbprmnc(mnmax_nyq0,ns) ,
255  8 ppprmnc(mnmax_nyq0,ns) , hotdmnc(mnmax_nyq0,ns) ,
256 #endif
257  9 stat=istat)
258  IF (lasym) THEN
259  ALLOCATE (gmns(mnmax_nyq0,ns), bmns(mnmax_nyq0,ns),
260  1 bsubumns(mnmax_nyq0,ns), bsubvmns(mnmax_nyq0,ns),
261  2 bsubsmnc(mnmax_nyq0,ns), bsupumns(mnmax_nyq0,ns),
262  3 bsupvmns(mnmax_nyq0,ns),
263  4 currumns(mnmax_nyq0,ns), currvmns(mnmax_nyq0,ns),
264 #ifdef _ANIMEC
265  4 sigmns(mnmax_nyq0,ns) ,
266  5 taumns(mnmax_nyq0,ns) , pparmns(mnmax_nyq0,ns) ,
267  6 ppermns(mnmax_nyq0,ns) , pbprmns(mnmax_nyq0,ns) ,
268  7 ppprmns(mnmax_nyq0,ns) , hotdmns(mnmax_nyq0,ns) ,
269 #endif
270  8 stat=istat)
271 #ifdef _ANIMEC
272  ALLOCATE (sigma_ana(ns,nznt) ,tau_ana(ns,nznt) ,densita(ns,nznt),
273  1 ppara(ns,nznt) ,pperpa(ns,nznt) ,pbprima(ns,nznt),
274  2 ppprima(ns,nznt), stat=istat)
275 #endif
276  END IF
277  IF (istat .ne. 0) stop 'Error allocating arrays in VMEC WROUT'
278 
279 ! IF (nextcur .eq. 0) THEN
280 ! DO j = SIZE(extcur), 1, -1
281 ! IF (extcur(j) .ne. zero) THEN
282 ! nextcur = j
283 ! EXIT
284 ! END IF
285 ! END DO
286 ! END IF
287 
288 ! ftol info evaluated here!
289  indx1=maxloc(ns_array)
290  ftolx1=ftol_array(indx1(1))
291 
292 ! NYQUIST FREQUENCY REQUIRES FACTOR OF 1/2
293  IF (lnyquist) THEN
294  IF (mnyq .ne. 0) cosmui(:,mnyq) = p5*cosmui(:,mnyq)
295  IF (nnyq .ne. 0) cosnv(:,nnyq) = p5*cosnv(:,nnyq)
296  END IF
297  wout_file = version_
298  READ (wout_file, *) vversion
299 
300 #ifdef NETCDF
301  wout_file = 'wout_' // trim(input_extension) // '.nc'
302  CALL cdf_open(nwout,wout_file,'w',iwout0)
303  IF (iwout0 .ne. 0) stop 'Error opening wout.nc file VMEC WROUT'
304 
305 !================================
306 ! Define Variables
307 !================================
308 ! Scalars
309  CALL cdf_define(nwout, vn_version, vversion)
310  CALL cdf_define(nwout, vn_extension, input_extension)
311  CALL cdf_define(nwout, vn_mgrid, mgrid_file)
312  CALL cdf_define(nwout, vn_pcurr_type, pcurr_type)
313  CALL cdf_define(nwout, vn_pmass_type, pmass_type)
314  CALL cdf_define(nwout, vn_piota_type, piota_type)
315  CALL cdf_define(nwout, vn_magen, wb)
316  CALL cdf_define(nwout, vn_therm, wp)
317  CALL cdf_define(nwout, vn_gam, gamma)
318  CALL cdf_define(nwout, vn_maxr, rmax_surf)
319  CALL cdf_define(nwout, vn_minr, rmin_surf)
320  CALL cdf_define(nwout, vn_maxz, zmax_surf)
321  CALL cdf_define(nwout, vn_fp, nfp)
322  CALL cdf_define(nwout, vn_radnod, ns)
323  CALL cdf_define(nwout, vn_polmod, mpol)
324  CALL cdf_define(nwout, vn_tormod, ntor)
325  CALL cdf_define(nwout, vn_maxmod, mnmax)
326  CALL cdf_define(nwout, vn_maxmod_nyq, mnmax_nyq0)
327  CALL cdf_define(nwout, vn_maxit, iter2)
328  CALL cdf_define(nwout, vn_actit, itfsq)
329  CALL cdf_define(nwout, vn_asym, lasym)
330  CALL cdf_define(nwout, vn_recon, lrecon)
331  CALL cdf_define(nwout, vn_free, lfreeb)
332  CALL cdf_define(nwout, vn_rfp, lrfp)
333  CALL cdf_define(nwout, vn_error, ier_flag)
334  CALL cdf_define(nwout, vn_aspect, aspect)
335  CALL cdf_define(nwout, vn_beta, betatot)
336  CALL cdf_define(nwout, vn_pbeta, betapol)
337  CALL cdf_define(nwout, vn_tbeta, betator)
338  CALL cdf_define(nwout, vn_abeta, betaxis)
339  CALL cdf_define(nwout, vn_b0, b0)
340  CALL cdf_define(nwout, vn_rbt0, rbtor0)
341  CALL cdf_define(nwout, vn_rbt1, rbtor)
342  CALL cdf_define(nwout, vn_sgs, nint(signgs))
343  CALL cdf_define(nwout, vn_lar, ionlarmor)
344  CALL cdf_define(nwout, vn_modb, volavgb)
345  CALL cdf_define(nwout, vn_ctor, ctor)
346  CALL cdf_define(nwout, vn_amin, aminor_p)
347  CALL cdf_define(nwout, vn_rmaj, rmajor_p)
348  CALL cdf_define(nwout, vn_vol, volume_p)
349  CALL cdf_define(nwout, vn_ftolv, ftolx1)
350  CALL cdf_define(nwout, vn_fsql, fsql)
351  CALL cdf_define(nwout, vn_fsqr, fsqr)
352  CALL cdf_define(nwout, vn_fsqz, fsqz)
353 
354  CALL cdf_define(nwout, vn_nextcur, nextcur)
355  CALL cdf_define(nwout, vn_extcur, extcur(1:nextcur),
356  1 dimname=currg)
357  CALL cdf_define(nwout, vn_mgmode, mgrid_mode)
358  IF (lfreeb) THEN
359  CALL cdf_define(nwout, vn_maxpot, mnpd)
360  CALL cdf_define(nwout, vn_flp, nobser)
361  CALL cdf_define(nwout, vn_nobd, nobd)
362  CALL cdf_define(nwout, vn_nbset, nbsets)
363  IF (nbsets .gt. 0)
364  1 CALL cdf_define(nwout,vn_nbfld,nbfld(1:nbsets))
365  END IF
366 
367  IF (.not.lwrite) GO TO 800
368 
369 ! 1D Arrays
370 
371  CALL cdf_define(nwout, vn_pmod, xm, dimname=mn1dim)
372  CALL cdf_setatt(nwout, vn_pmod, ln_pmod)
373  CALL cdf_define(nwout, vn_tmod, xn, dimname=mn1dim)
374  CALL cdf_setatt(nwout, vn_tmod, ln_tmod)
375  CALL cdf_define(nwout, vn_pmod_nyq, xm_nyq0, dimname=mn2dim)
376  CALL cdf_setatt(nwout, vn_pmod_nyq, ln_pmod_nyq)
377  CALL cdf_define(nwout, vn_tmod_nyq, xn_nyq0, dimname=mn2dim)
378  CALL cdf_setatt(nwout, vn_tmod_nyq, ln_tmod_nyq)
379 
380  CALL cdf_define(nwout, vn_racc, raxis_cc(0:ntor),
381  1 dimname=(/'n_tor'/))
382  CALL cdf_setatt(nwout, vn_racc, ln_racc)
383  CALL cdf_define(nwout, vn_zacs, zaxis_cs(0:ntor),
384  1 dimname=(/'n_tor'/))
385  CALL cdf_setatt(nwout, vn_zacs, ln_zacs)
386  IF (lasym) THEN
387  CALL cdf_define(nwout, vn_racs, raxis_cs(0:ntor),
388  1 dimname=(/'n_tor'/))
389  CALL cdf_setatt(nwout, vn_racs, ln_racs)
390  CALL cdf_define(nwout, vn_zacc, zaxis_cc(0:ntor),
391  1 dimname=(/'n_tor'/))
392  CALL cdf_setatt(nwout, vn_zacc, ln_zacc)
393  END IF
394 
395  j = SIZE(am)-1
396  CALL cdf_define(nwout, vn_am, am(0:j),
397  1 dimname=(/'preset'/))
398  j = SIZE(ac)-1
399  CALL cdf_define(nwout, vn_ac, ac(0:j),
400  1 dimname=(/'preset'/))
401  j = SIZE(ai)-1
402  CALL cdf_define(nwout, vn_ai, ai(0:j),
403  1 dimname=(/'preset'/))
404 
405  j = SIZE(am_aux_s)
406  CALL cdf_define(nwout, vn_am_aux_s, am_aux_s(1:j),
407  1 dimname=(/'ndfmax'/))
408  j = SIZE(am_aux_f)
409  CALL cdf_define(nwout, vn_am_aux_f, am_aux_f(1:j),
410  1 dimname=(/'ndfmax'/))
411  j = SIZE(ai_aux_s)
412  CALL cdf_define(nwout, vn_ai_aux_s, ai_aux_s(1:j),
413  1 dimname=(/'ndfmax'/))
414  j = SIZE(ai_aux_f)
415  CALL cdf_define(nwout, vn_ai_aux_f, ai_aux_f(1:j),
416  1 dimname=(/'ndfmax'/))
417  j = SIZE(ac_aux_s)
418  CALL cdf_define(nwout, vn_ac_aux_s, ac_aux_s(1:j),
419  1 dimname=(/'ndfmax'/))
420  j = SIZE(ac_aux_f)
421  CALL cdf_define(nwout, vn_ac_aux_f, ac_aux_f(1:j),
422  1 dimname=(/'ndfmax'/))
423 
424 
425  CALL cdf_define(nwout, vn_iotaf, iotaf(1:ns),
426  1 dimname=r1dim)
427  CALL cdf_setatt(nwout, vn_iotaf, ln_iotaf)
428 
429  qfact=huge(qfact)
430  WHERE (iotaf(1:ns) .NE. zero) qfact=one/iotaf(1:ns)
431 
432  CALL cdf_define(nwout, vn_qfact, qfact(1:ns),
433  1 dimname=r1dim)
434  CALL cdf_setatt(nwout, vn_iotaf, ln_qfact)
435  CALL cdf_define(nwout, vn_presf, presf,
436  1 dimname=r1dim)
437  CALL cdf_setatt(nwout, vn_presf, ln_presf, units='Pa')
438  CALL cdf_define(nwout, vn_phi, phi,
439  1 dimname=r1dim)
440  CALL cdf_setatt(nwout, vn_phi, ln_phi, units='wb')
441  CALL cdf_define(nwout, vn_phipf,
442  1 phipf, dimname=r1dim)
443  CALL cdf_setatt(nwout, vn_phipf, ln_phipf)
444  CALL cdf_define(nwout, vn_chi, chi,
445  1 dimname=r1dim)
446  CALL cdf_setatt(nwout, vn_chi, ln_chi, units='wb')
447  CALL cdf_define(nwout, vn_chipf,
448  1 phipf, dimname=r1dim)
449  CALL cdf_setatt(nwout, vn_chipf, ln_chipf)
450  CALL cdf_define(nwout, vn_jcuru,
451  1 jcuru, dimname=r1dim)
452  CALL cdf_define(nwout, vn_jcurv,
453  1 jcurv, dimname=r1dim)
454 
455  CALL cdf_define(nwout, vn_iotah, iotas(1:ns),
456  1 dimname=r1dim)
457  CALL cdf_setatt(nwout, vn_iotah, ln_iotah)
458  CALL cdf_define(nwout, vn_mass, mass,
459  1 dimname=r1dim)
460  CALL cdf_setatt(nwout, vn_mass, ln_mass)
461  CALL cdf_define(nwout, vn_presh, pres(1:ns),
462  1 dimname=r1dim)
463  CALL cdf_setatt(nwout, vn_presh, ln_presh, units='Pa')
464  CALL cdf_define(nwout, vn_betah, beta_vol,
465  1 dimname=r1dim)
466  CALL cdf_define(nwout, vn_buco, buco,
467  1 dimname=r1dim)
468  CALL cdf_define(nwout, vn_bvco, bvco,
469  1 dimname=r1dim)
470  CALL cdf_define(nwout, vn_vp, vp(1:ns),
471  1 dimname=r1dim)
472  CALL cdf_define(nwout, vn_specw, specw,
473  1 dimname=r1dim)
474  CALL cdf_define(nwout, vn_phip,
475  1 phips(1:ns), dimname=r1dim)
476  CALL cdf_define(nwout, vn_overr,
477  2 overr(1:ns), dimname=r1dim)
478 
479  CALL cdf_define(nwout, vn_jdotb, jdotb,
480  1 dimname=r1dim)
481  CALL cdf_define(nwout, vn_bdotb, bdotb,
482  1 dimname=r1dim)
483  CALL cdf_define(nwout, vn_bgrv, bdotgradv,
484  1 dimname=r1dim)
485 
486  CALL cdf_define(nwout, vn_merc, dmerc,
487  1 dimname=r1dim)
488  CALL cdf_define(nwout, vn_mshear, dshear,
489  1 dimname=r1dim)
490  CALL cdf_define(nwout, vn_mwell, dwell,
491  1 dimname=r1dim)
492  CALL cdf_define(nwout, vn_mcurr, dcurr,
493  1 dimname=r1dim)
494  CALL cdf_define(nwout, vn_mgeo,
495  1 dgeod, dimname=r1dim)
496  CALL cdf_define(nwout, vn_equif,
497  1 equif, dimname=r1dim)
498 
499  CALL cdf_define(nwout, vn_fsq, fsqt(1:nstore_seq),
500  1 dimname=(/'time'/))
501  CALL cdf_define(nwout, vn_wdot, wdot(1:nstore_seq),
502  1 dimname=(/'time'/))
503 
504  IF (lfreeb) THEN
505  CALL cdf_define(nwout, vn_potsin, potvac(1:mnpd),
506  1 dimname=mnpotdim)
507  CALL cdf_setatt(nwout, vn_potsin, ln_potsin)
508  CALL cdf_define(nwout, vn_xmpot, xmpot(1:mnpd),
509  1 dimname=mnpotdim)
510  CALL cdf_define(nwout, vn_xnpot, xnpot(1:mnpd),
511  1 dimname=mnpotdim)
512  IF (lasym) THEN
513  CALL cdf_define(nwout, vn_potcos,
514  1 potvac(1+mnpd:2*mnpd), dimname=mnpotdim)
515  CALL cdf_setatt(nwout, vn_potcos, ln_potcos)
516  END IF
517 
518  IF (nextcur.gt.0 .and. ALLOCATED(curlabel)) THEN
519  CALL cdf_define(nwout, vn_curlab,
520  1 curlabel(1:nextcur), dimname=currl)
521  END IF
522  ENDIF
523 
524 ! 2D Arrays
525  CALL cdf_define(nwout, vn_rmnc, rmnc, dimname=r2dim)
526  CALL cdf_setatt(nwout, vn_rmnc, ln_rmnc, units='m')
527  CALL cdf_define(nwout, vn_zmns, zmns, dimname=r2dim)
528  CALL cdf_setatt(nwout, vn_zmns, ln_zmns, units='m')
529  CALL cdf_define(nwout, vn_lmns, lmns, dimname=r2dim)
530  CALL cdf_setatt(nwout, vn_lmns, ln_lmns)
531  CALL cdf_define(nwout, vn_gmnc, gmnc, dimname=r3dim)
532  CALL cdf_setatt(nwout, vn_gmnc, ln_gmnc)
533  CALL cdf_define(nwout, vn_bmnc, bmnc, dimname=r3dim)
534  CALL cdf_setatt(nwout, vn_bmnc, ln_bmnc)
535  CALL cdf_define(nwout, vn_bsubumnc, bsubumnc, dimname=r3dim)
536  CALL cdf_setatt(nwout, vn_bsubumnc, ln_bsubumnc)
537  CALL cdf_define(nwout, vn_bsubvmnc, bsubvmnc, dimname=r3dim)
538  CALL cdf_setatt(nwout, vn_bsubvmnc, ln_bsubvmnc)
539  CALL cdf_define(nwout, vn_bsubsmns, bsubsmns, dimname=r3dim)
540  CALL cdf_setatt(nwout, vn_bsubsmns, ln_bsubsmns)
541 
542  CALL cdf_define(nwout, vn_currumnc, currumnc, dimname=r3dim) !MRC 8-12-16
543  CALL cdf_setatt(nwout, vn_currumnc, ln_currumnc)
544  CALL cdf_define(nwout, vn_currvmnc, currvmnc, dimname=r3dim)
545  CALL cdf_setatt(nwout, vn_currvmnc, ln_currvmnc)
546 
547  IF (lfreeb) THEN
548  CALL cdf_define(nwout, vn_bsubumnc_sur, bsubumnc_sur, &
549  & dimname=mn2dim)
550  CALL cdf_setatt(nwout, vn_bsubumnc_sur, ln_bsubumnc_sur)
551  CALL cdf_define(nwout, vn_bsubvmnc_sur, bsubvmnc_sur, &
552  & dimname=mn2dim)
553  CALL cdf_setatt(nwout, vn_bsubvmnc_sur, ln_bsubvmnc_sur)
554  CALL cdf_define(nwout, vn_bsupumnc_sur, bsupumnc_sur, &
555  & dimname=mn2dim)
556  CALL cdf_setatt(nwout, vn_bsupumnc_sur, ln_bsupumnc_sur)
557  CALL cdf_define(nwout, vn_bsupvmnc_sur, bsupvmnc_sur, &
558  & dimname=mn2dim)
559  CALL cdf_setatt(nwout, vn_bsupvmnc_sur, ln_bsupvmnc_sur)
560  END IF
561 
562 ! ELIMINATE THESE EVENTUALLY: DON'T NEED THEM - CAN COMPUTE FROM GSQRT
563  CALL cdf_define(nwout, vn_bsupumnc, bsupumnc, dimname=r3dim)
564  CALL cdf_define(nwout, vn_bsupvmnc, bsupvmnc, dimname=r3dim)
565 ! IF (lfreeb) THEN
566 ! CALL cdf_define(nwout, vn_rbc, rbc,
567 ! 1 dimname=(/'n_mode','m_mode'/))
568 ! CALL cdf_setatt(nwout, vn_rbc, ln_rbc, units='m')
569 ! CALL cdf_define(nwout, vn_zbs, zbs,
570 ! 1 dimname=(/'n_mode','m_mode'/))
571 ! CALL cdf_setatt(nwout, vn_zbs, ln_zbs, units='m')
572 ! IF (lasym) THEN
573 ! CALL cdf_define(nwout, vn_rbs, rbs,
574 ! 1 dimname=(/'n_mode','m_mode'/))
575 ! CALL cdf_define(nwout, vn_zbc, zbc,
576 ! 1 dimname=(/'n_mode','m_mode'/))
577 ! END IF
578 ! END IF
579 
580  IF (.NOT. lasym) GO TO 800
581 
582  CALL cdf_define(nwout, vn_rmns, rmns, dimname=r2dim)
583  CALL cdf_setatt(nwout, vn_rmns, ln_rmns, units='m')
584  CALL cdf_define(nwout, vn_zmnc, zmnc, dimname=r2dim)
585  CALL cdf_setatt(nwout, vn_zmnc, ln_zmnc, units='m')
586  CALL cdf_define(nwout, vn_lmnc, lmnc, dimname=r2dim)
587  CALL cdf_setatt(nwout, vn_lmnc, ln_lmnc)
588  CALL cdf_define(nwout, vn_gmns, gmns, dimname=r3dim)
589  CALL cdf_setatt(nwout, vn_gmns, ln_gmns)
590  CALL cdf_define(nwout, vn_bmns, bmns, dimname=r3dim)
591  CALL cdf_setatt(nwout, vn_bmns, ln_bmns)
592  CALL cdf_define(nwout, vn_bsubumns, bsubumns, dimname=r3dim)
593  CALL cdf_setatt(nwout, vn_bsubumns, ln_bsubumns)
594  CALL cdf_define(nwout, vn_bsubvmns, bsubvmns, dimname=r3dim)
595  CALL cdf_setatt(nwout, vn_bsubvmns, ln_bsubvmns)
596  CALL cdf_define(nwout, vn_bsubsmnc, bsubsmnc, dimname=r3dim)
597  CALL cdf_setatt(nwout, vn_bsubsmnc, ln_bsubsmnc)
598 
599  CALL cdf_define(nwout, vn_currumns, currumns, dimname=r3dim)
600  CALL cdf_setatt(nwout, vn_currumns, ln_currumns)
601  CALL cdf_define(nwout, vn_currvmns, currvmns, dimname=r3dim)
602  CALL cdf_setatt(nwout, vn_currvmns, ln_currvmns)
603 
604  IF (lfreeb) THEN
605  CALL cdf_define(nwout, vn_bsubumns_sur, bsubumns_sur, &
606  & dimname=mn2dim)
607  CALL cdf_setatt(nwout, vn_bsubumns_sur, ln_bsubumns_sur)
608  CALL cdf_define(nwout, vn_bsubvmns_sur, bsubvmns_sur, &
609  & dimname=mn2dim)
610  CALL cdf_setatt(nwout, vn_bsubvmns_sur, ln_bsubvmns_sur)
611  CALL cdf_define(nwout, vn_bsupumns_sur, bsupumns_sur, &
612  & dimname=mn2dim)
613  CALL cdf_setatt(nwout, vn_bsupumns_sur, ln_bsupumns_sur)
614  CALL cdf_define(nwout, vn_bsupvmns_sur, bsupvmns_sur, &
615  & dimname=mn2dim)
616  CALL cdf_setatt(nwout, vn_bsupvmns_sur, ln_bsupvmns_sur)
617  END IF
618 
619 ! ELIMINATE THESE EVENTUALLY: DON'T NEED THEM
620  CALL cdf_define(nwout, vn_bsupumns, bsupumns, dimname=r3dim)
621  CALL cdf_define(nwout, vn_bsupvmns, bsupvmns, dimname=r3dim)
622 
623  800 CONTINUE
624 
625 !================================
626 ! Write Variables
627 !================================
628 
629 ! Scalars
630  CALL cdf_write(nwout, vn_version, vversion)
631  CALL cdf_write(nwout, vn_extension, input_extension)
632  CALL cdf_write(nwout, vn_mgrid, mgrid_file)
633  CALL cdf_write(nwout, vn_pcurr_type, pcurr_type)
634  CALL cdf_write(nwout, vn_piota_type, piota_type)
635  CALL cdf_write(nwout, vn_pmass_type, pmass_type)
636  CALL cdf_write(nwout, vn_magen, wb)
637  CALL cdf_write(nwout, vn_therm, wp)
638  CALL cdf_write(nwout, vn_gam, gamma)
639  CALL cdf_write(nwout, vn_maxr, rmax_surf)
640  CALL cdf_write(nwout, vn_minr, rmin_surf)
641  CALL cdf_write(nwout, vn_maxz, zmax_surf)
642  CALL cdf_write(nwout, vn_fp, nfp)
643  CALL cdf_write(nwout, vn_radnod, ns)
644  CALL cdf_write(nwout, vn_polmod, mpol)
645  CALL cdf_write(nwout, vn_tormod, ntor)
646  CALL cdf_write(nwout, vn_maxmod, mnmax)
647  CALL cdf_write(nwout, vn_maxmod_nyq, mnmax_nyq0)
648  CALL cdf_write(nwout, vn_maxit, iter2)
649  CALL cdf_write(nwout, vn_actit, itfsq)
650  CALL cdf_write(nwout, vn_asym, lasym)
651  CALL cdf_write(nwout, vn_recon, lrecon)
652  CALL cdf_write(nwout, vn_free, lfreeb)
653  CALL cdf_write(nwout, vn_rfp, lrfp)
654  CALL cdf_write(nwout, vn_error, ier_flag)
655 !
656  CALL cdf_write(nwout, vn_aspect, aspect)
657  CALL cdf_write(nwout, vn_beta, betatot)
658  CALL cdf_write(nwout, vn_pbeta, betapol)
659  CALL cdf_write(nwout, vn_tbeta, betator)
660  CALL cdf_write(nwout, vn_abeta, betaxis)
661  CALL cdf_write(nwout, vn_b0, b0)
662  CALL cdf_write(nwout, vn_rbt0, rbtor0)
663  CALL cdf_write(nwout, vn_rbt1, rbtor)
664  CALL cdf_write(nwout, vn_sgs, nint(signgs))
665  CALL cdf_write(nwout, vn_lar, ionlarmor)
666  CALL cdf_write(nwout, vn_modb, volavgb)
667  CALL cdf_write(nwout, vn_ctor, ctor/mu0)
668  CALL cdf_write(nwout, vn_amin, aminor_p)
669  CALL cdf_write(nwout, vn_rmaj, rmajor_p)
670  CALL cdf_write(nwout, vn_vol, volume_p)
671  CALL cdf_write(nwout, vn_ftolv, ftolx1)
672  CALL cdf_write(nwout, vn_fsql, fsql)
673  CALL cdf_write(nwout, vn_fsqr, fsqr)
674  CALL cdf_write(nwout, vn_fsqz, fsqz)
675 
676  CALL cdf_write(nwout, vn_nextcur, nextcur)
677  IF (nextcur .gt. 0) THEN
678  CALL cdf_write(nwout, vn_extcur, extcur(1:nextcur))
679  CALL cdf_write(nwout, vn_mgmode, mgrid_mode)
680  ENDIF
681  IF (lfreeb) THEN
682  CALL cdf_write(nwout, vn_flp, nobser)
683  CALL cdf_write(nwout, vn_maxpot, mnpd)
684  CALL cdf_write(nwout, vn_nobd, nobd)
685  CALL cdf_write(nwout, vn_nbset, nbsets)
686  IF (nextcur.gt.0 .and. ALLOCATED(curlabel))
687  1 CALL cdf_write(nwout, vn_curlab, curlabel(1:nextcur))
688  END IF
689 
690 ! 1D Arrays
691  IF (nbsets .gt. 0) CALL cdf_write(nwout,vn_nbfld,nbfld(1:nbsets))
692 
693  IF (.not.lwrite) GO TO 940 !change 970 to 940, J.Geiger
694  !skip the next cdf_write-calls
695  CALL cdf_write(nwout, vn_pmod, xm)
696  CALL cdf_write(nwout, vn_tmod, xn)
697  CALL cdf_write(nwout, vn_pmod_nyq, xm_nyq0)
698  CALL cdf_write(nwout, vn_tmod_nyq, xn_nyq0)
699 
700  IF (lfreeb) THEN
701  CALL cdf_write(nwout, vn_potsin, potvac(1:mnpd))
702  IF (lasym) &
703  & CALL cdf_write(nwout, vn_potcos, potvac(1+mnpd:2*mnpd))
704  CALL cdf_write(nwout, vn_xmpot, xmpot)
705  CALL cdf_write(nwout, vn_xnpot, xnpot)
706  END IF
707 
708 940 CONTINUE ! before closing, write the initial part of the wouttxt-file
709 #endif
710  IF (lwouttxt) THEN
711  wout2_file = 'wout_'//trim(input_extension) // '.txt'
712  nwout2 = nwout0
713  CALL safe_open(nwout2, iwout0, wout2_file,
714  1 'replace', 'formatted')
715  IF (iwout0 .ne. 0) stop 'Error opening WOUT.txt file in WROUT'
716 
717  IF (lasym) THEN
718  iasym = 1 ! asymmetric mode
719  ELSE
720  iasym = 0
721  END IF
722 
723 !
724 ! Insert version information into wout file. This will be parsed in
725 ! read_wout_file to return the real value version_ to check the version number.
726 !
727 #if defined(_ANIMEC)
728  WRITE (nwout2, '(a15,a,a)') 'VMEC VERSION = ', version_,
729  1 '_ANIMEC'
730 #elif defined(_FLOW)
731  WRITE (nwout2, '(a15,a,a)') 'VMEC VERSION = ', version_,'_FLOW'
732 #else
733  WRITE (nwout2, '(a15,a)') 'VMEC VERSION = ', version_
734 #endif
735 
736 #ifdef _ANIMEC
737  WRITE (nwout2, *) wb, wpar, gamma, pfac,
738 #else
739  WRITE (nwout2, *) wb, wp, gamma, 1,
740 #endif
741  1 rmax_surf, rmin_surf, zmax_surf
742 
743  WRITE (nwout2, *) nfp, ns, mpol, ntor, mnmax, mnmax_nyq0,
744  1 itfsq, iter2, iasym, 0, ier_flag
745 
746  WRITE (nwout2, *) 0, 0, nbsets, nobd, nextcur, nstore_seq
747  IF (nbsets .gt. 0) WRITE (nwout2, *) (nbfld(i),i=1,nbsets)
748  WRITE (nwout2, '(a)') mgrid_file
749 
750  IF (.not. lwrite) GOTO 950 ! J Geiger: At this point only closing of
751  ! txt- and nc-file is needed
752  pres(1) = pres(2)
753  END IF
754 !---------------------DEC$ ENDIF
755 
756  IF (.not. lwrite) GOTO 970 ! J Geiger: in case lwouttxt is not true
757  ! jump to close nc-file
758  ALLOCATE (xfinal(neqs), stat=js)
759  IF (js .NE. 0) stop 'Allocation error for xfinal in WROUT!'
760  xfinal = xc
761 #ifdef _HBANGLE
762  CALL getrz(xfinal)
763 #else
764 !
765 ! MUST CONVERT m=1 MODES... FROM INTERNAL TO PHYSICAL FORM
766 ! Extrapolation of m=0 Lambda (cs) modes, which are not evolved at j=1, done in CONVERT
767 !
768  lk = ns*ntor1
769  IF (lthreed) CALL convert_sym (xfinal(1+mns*(rss-1)+lk),
770  1 xfinal(1+irzloff+mns*(zcs-1)+lk))
771  IF (lasym) CALL convert_asym (xfinal(1+mns*(rsc-1)+lk),
772  1 xfinal(1+irzloff+mns*(zcc-1)+lk))
773 #endif
774 !
775 ! CONVERT TO rmnc, zmns, lmns, etc EXTERNAL representation (without internal mscale, nscale)
776 ! IF B^v ~ phip + lamu, MUST DIVIDE BY phipf(js) below to maintain old-style format
777 ! THIS COULD BE A PROBLEM FOR RFP WHERE PHIPF->0 INSIDE THE PLASMA!
778 !
779  radius1: DO js = 1, ns
780 
781  CALL convert (rmnc1, zmns1, lmns1, rmns1, zmnc1, lmnc1,
782  1 xfinal, js)
783 
784  rmnc(:,js) = rmnc1(:)
785  zmns(:,js) = zmns1(:)
786  lmns(:,js) = (lmns1(:)/phipf(js)) * lamscale
787  IF (lasym) THEN
788  rmns(:,js) = rmns1(:)
789  zmnc(:,js) = zmnc1(:)
790  lmnc(:,js) = (lmnc1(:)/phipf(js)) * lamscale
791  END IF
792 
793  END DO radius1
794 
795  DEALLOCATE (xfinal)
796 
797 !
798 ! INTERPOLATE LAMBDA ONTO HALF-MESH FOR BACKWARDS CONSISTENCY WITH EARLIER VERSIONS OF VMEC
799 ! AND SMOOTHS POSSIBLE UNPHYSICAL "WIGGLE" ON RADIAL MESH
800 !
801 
802  WHERE (nint(xm) .le. 1) lmns(:,1) = lmns(:,2)
803  DO js = ns,2,-1
804  WHERE (mod(nint(xm),2) .eq. 0)
805  lmns(:,js) = p5*(lmns(:,js) + lmns(:,js-1))
806  ELSEWHERE
807  lmns(:,js) = p5*(sm(js)*lmns(:,js) + sp(js-1)*lmns(:,js-1))
808  END WHERE
809  END DO
810 
811  lmns(:,1) = 0
812  raxis_cc(0:ntor) = rmnc(1:ntor+1,1)
813  zaxis_cs(0:ntor) = zmns(1:ntor+1,1)
814 
815  IF (.NOT.lasym) GOTO 900
816 
817  WHERE (nint(xm) .le. 1) lmnc(:,1) = lmnc(:,2)
818  DO js = ns,2,-1
819  WHERE (mod(nint(xm),2) .eq. 0)
820  lmnc(:,js) = p5*(lmnc(:,js) + lmnc(:,js-1))
821  ELSEWHERE
822  lmnc(:,js) = p5*(sm(js)*lmnc(:,js) + sp(js-1)*lmnc(:,js-1))
823  END WHERE
824  END DO
825 
826  lmnc(:,1) = 0;
827  raxis_cs(0:ntor) = rmns(1:ntor+1,1)
828  zaxis_cc(0:ntor) = zmnc(1:ntor+1,1)
829 
830  900 CONTINUE
831 
832 #ifdef _ANIMEC
833 !... CALCULATE RADIAL DERIVATIVES OF HOT PARTICLE PRESSURE TERMS
834 !... STORE IN ARRAYS pm AND pd PREVIOUSLY USED IN PRESSURE AND EQFOR
835  eps = epsilon(eps)
836  DO js=2,ns-1
837  pd(js) = ohs * (pres(js+1) * phot(js+1) - pres(js) * phot(js))
838  pmap(js) = ohs * (tpotb(js+1) - tpotb(js))
839  END DO
840 !... INTERPOLATE (EXTRAPOLATE) TO HALF INTEGER MESH
841  pdh = c1p5 * pd(2) - p5 * pd(3)
842  pmh = c1p5 * pmap(2) - p5 * pmap(3)
843  pde = c1p5 * pd(ns-1) - p5 * pd(ns-2)
844  pme = c1p5 * pmap(ns-1) - p5 * pmap(ns-2)
845  DO js=ns-2,2,-1
846  pd(js+1) = p5*(pd(js+1) + pd(js)) / (pres(js+1)*phot(js+1)+eps)
847  pmap(js+1) = p5 * (pmap(js+1) + pmap(js))
848  END DO
849  pd(2) = pdh / (pres(2)*phot(2)+eps)
850  pd(ns) = pde / (pres(ns)*phot(ns)+eps)
851  pmap(2) = pmh
852  pmap(ns) = pme
853 !ALTERNATE EXTRAPOLATION
854  pd(2) = 2*pd(3) - pd(4)
855  pd(ns) = 2*pd(ns-1) - pd(ns-2)
856 
857 !CALCULATE HOT PARTICLE PARALLEL AND PERPENDICULAR PRESSURE GRADIENT; DENSITY
858  DO 20 js = 2, ns
859  hotdam = pres(js) * phot(js) / sqrt(tpotb(js)+eps)
860  DO 10 lk = 1, nznt
861 !
862  omtbc = one - tpotb(js) * onembc(js,lk)
863  optbc = one + tpotb(js) * onembc(js,lk)
864  IF (onembc(js,lk) <= zero) THEN
865  densit(js,lk)= (ppar(js,lk) - pres(js))*hotdam /
866  & (pres(js)*phot(js)+eps)
867  pbprim(js,lk) = (ppar(js,lk) -pres(js)) *
868  & (pd(js) + onembc(js,lk) * pmap(js) / (omtbc+eps))
869  ppprim(js,lk) = (pperp(js,lk)-pres(js)) *
870  & (pd(js) + optbc * pmap(js) / (omtbc * tpotb(js)+eps))
871  ELSE
872  densit(js,lk) = hotdam * (one - onembc(js,lk)) *
873  & (optbc - 2*(tpotb(js)*onembc(js,lk))**c1p5) / (omtbc*optbc+eps)
874  pbprim(js,lk) = (ppar(js,lk) -pres(js)) * pd(js) +
875  & ( 2 * tpotb(js) * onembc(js,lk)**2 * (ppar(js,lk)-pres(js))
876  & + pres(js)*phot(js)*(one-onembc(js,lk))*onembc(js,lk)*(one -5
877  & *(tpotb(js)*onembc(js,lk))**c1p5))* pmap(js) / (omtbc*optbc+eps)
878  ppprim(js,lk) = (pperp(js,lk)-pres(js)) * pd(js) +
879  & ((pperp(js,lk)-pres(js))*(one+3*(tpotb(js)*onembc(js,lk))**2) /
880  & (tpotb(js)+eps)+ pres(js)*phot(js)*tpotb(js)
881  & *(one-onembc(js,lk))**2
882  & * onembc(js,lk)*(two*optbc-sqrt(tpotb(js)*onembc(js,lk))*(7.5
883  & - 3.5_dp*(tpotb(js)*onembc(js,lk))**2))/(omtbc*optbc+eps))
884  & * pmap(js)/ (omtbc * optbc + eps)
885  END IF
886  10 END DO
887  20 END DO
888 #endif
889 !SPH100209: COMPUTE |B| = SQRT(|B|**2) and store in bsq, bsqa
890  DO js = 2, ns
891  bsq(js,:nznt) = sqrt(2*abs(bsq(js,:nznt)-pres(js)))
892  END DO
893 
894  tmult = p5/r0scale**2
895 !SPH: FIXED THIS 03-05-07 TO CALL symmetrization routine
896  IF (lasym) THEN
897 !Changed integration norm in fixaray, SPH012314
898  tmult = 2*tmult
899  bsubs(1,:) = 0
900  CALL symoutput (bsq, gsqrt, bsubu, bsubv, bsupu,
901  1 bsupv, bsubs,
902 #ifdef _ANIMEC
903  2 ppar , pperp , densit ,
904  3 sigma_an , tau_an , pbprim , ppprim ,
905 #endif
906  4 bsqa, gsqrta, bsubua, bsubva, bsupua,
907  5 bsupva, bsubsa
908 #ifdef _ANIMEC
909  6 ,ppara , pperpa , densita,
910  7 sigma_ana, tau_ana, pbprima, ppprima
911 #endif
912  8 )
913 
914  IF (lfreeb) THEN !MRC 10-15-15
915  CALL symoutput_sur(bsubu_sur, bsubv_sur, &
916  & bsupu_sur, bsupv_sur, &
917  & bsubua_sur, bsubva_sur, &
918  & bsupua_sur, bsupva_sur)
919  END IF
920  END IF
921 
922 ! DO js = 2, ns
923 ! WRITE (200, *) 'JS: ', js, 'BSUBU, BSUBV'
924 ! WRITE (200, '(1p,6e12.4)') bsubu(js,:), bsubv(js,:)
925 ! END DO
926 
927  radius2: DO js = 2, ns
928  gmn = 0
929  bmn = 0
930  bsubumn = 0
931  bsubvmn = 0
932  bsubsmn = 0
933  bsupumn = 0
934  bsupvmn = 0
935 
936  mn2: DO mn = 1, mnmax_nyq0
937  n = nint(xn_nyq0(mn))/nfp
938  m = nint(xm_nyq0(mn))
939  n1 = abs(n)
940  dmult = mscale(m)*nscale(n1)*tmult
941  IF (m.eq.0 .or. n.eq.0) dmult = 2*dmult
942  sgn = sign(1, n)
943  lk = 0
944  DO j = 1, ntheta2
945  DO k = 1, nzeta
946  lk = lk + 1
947  tcosi = dmult*(cosmui(j,m)*cosnv(k,n1) +
948  1 sgn*sinmui(j,m)*sinnv(k,n1)) !cos(mu - nv)
949  tsini = dmult*(sinmui(j,m)*cosnv(k,n1) -
950  1 sgn*cosmui(j,m)*sinnv(k,n1)) !sin(mu - nv)
951  bmn(mn) = bmn(mn) + tcosi*bsq(js,lk)
952  gmn(mn) = gmn(mn) + tcosi*gsqrt(js,lk)
953  bsubumn(mn) = bsubumn(mn) + tcosi*bsubu(js,lk)
954  bsubvmn(mn) = bsubvmn(mn) + tcosi*bsubv(js,lk)
955  bsubsmn(mn) = bsubsmn(mn) + tsini*bsubs(js,lk)
956  bsupumn(mn) = bsupumn(mn) + tcosi*bsupu(js,lk)
957  bsupvmn(mn) = bsupvmn(mn) + tcosi*bsupv(js,lk)
958 #ifdef _ANIMEC
959  pparmn(mn) = pparmn(mn) + tcosi*
960  1 (ppar(js,lk)-pres(js))
961  ppermn(mn) = ppermn(mn) + tcosi*
962  1 (pperp(js,lk)-pres(js))
963  sigmn(mn) = sigmn(mn) + tcosi*sigma_an(js,lk)
964  taumn(mn) = taumn(mn) + tcosi*tau_an(js,lk)
965  pbprmn(mn) = pbprmn(mn) + tcosi*pbprim(js,lk)
966  ppprmn(mn) = ppprmn(mn) + tcosi*ppprim(js,lk)
967  hotdmn(mn) = hotdmn(mn) + tcosi*densit(js,lk)
968 #endif
969  END DO
970  END DO
971  END DO mn2
972 
973  IF (js .eq. ns/2) bmodmn = bmn(1:mnmax)
974  IF (js .eq. ns) bmodmn1 = bmn(1:mnmax)
975  gmnc(:,js) = gmn(:)
976  bmnc(:,js) = bmn(:)
977  bsubumnc(:,js) = bsubumn(:)
978  bsubvmnc(:,js) = bsubvmn(:)
979  bsubsmns(:,js) = bsubsmn(:)
980  bsupumnc(:,js) = bsupumn(:)
981  bsupvmnc(:,js) = bsupvmn(:)
982 #ifdef _ANIMEC
983  pparmnc(:,js) = pparmn(:)
984  ppermnc(:,js) = ppermn(:)
985  sigmnc(:,js) = sigmn(:)
986  taumnc(:,js) = taumn(:)
987  pbprmnc(:,js) = pbprmn(:)
988  ppprmnc(:,js) = ppprmn(:)
989  hotdmnc(:,js) = hotdmn(:)
990 #endif
991  END DO radius2
992 
993  IF (lfreeb) THEN !MRC 10-15-15
994  bsubumnc_sur = 0
995  bsubvmnc_sur = 0
996  bsupumnc_sur = 0
997  bsupvmnc_sur = 0
998  DO mn = 1, mnmax_nyq0
999  n = nint(xn_nyq0(mn))/nfp
1000  m = nint(xm_nyq0(mn))
1001  n1 = abs(n)
1002  dmult = mscale(m)*nscale(n1)*tmult
1003  IF (m.eq.0 .or. n.eq.0) dmult = 2*dmult
1004  sgn = sign(1, n)
1005  lk = 0
1006  DO j = 1, ntheta2
1007  DO k = 1, nzeta
1008  lk = lk + 1
1009  tcosi = dmult*(cosmui(j,m)*cosnv(k,n1) +
1010  1 sgn*sinmui(j,m)*sinnv(k,n1))
1011  bsubumnc_sur(mn) = bsubumnc_sur(mn) &
1012  & + tcosi*bsubu_sur(lk)
1013  bsubvmnc_sur(mn) = bsubvmnc_sur(mn) &
1014  & + tcosi*bsubv_sur(lk)
1015  bsupumnc_sur(mn) = bsupumnc_sur(mn) &
1016  & + tcosi*bsupu_sur(lk)
1017  bsupvmnc_sur(mn) = bsupvmnc_sur(mn) &
1018  & + tcosi*bsupv_sur(lk)
1019  END DO
1020  END DO
1021  END DO
1022  END IF
1023 
1024  gmnc(:,1) = 0; bmnc(:,1) = 0;
1025  bsubumnc(:,1) = 0
1026  bsubvmnc(:,1) = 0
1027  bsubsmns(:,1) = 2*bsubsmns(:,2) - bsubsmns(:,3)
1028  bsupumnc(:,1) = 0; bsupvmnc(:,1) = 0
1029 
1030 #ifdef _ANIMEC
1031  hotdmnc(:,1) = 0; pparmnc(:,1) = 0; ppermnc(:,1) = 0
1032  pbprmnc(:,1) = 0; ppprmnc(:,1) = 0
1033  sigmnc(:,1) = 0; taumnc(:,1) = 0
1034 #endif
1035 
1036  IF (.not.lasym) GO TO 200
1037 
1038  radius3: DO js = 2, ns
1039  gmn = 0
1040  bmn = 0
1041  bsubumn = 0
1042  bsubvmn = 0
1043  bsubsmn = 0
1044  bsupumn = 0
1045  bsupvmn = 0
1046 #ifdef _ANIMEC
1047  pparmn = 0
1048  ppermn = 0
1049  sigmn = 0
1050  taumn = 0
1051  pbprmn = 0
1052  ppprmn = 0
1053  hotdmn = 0
1054 #endif
1055  mn3: DO mn = 1, mnmax_nyq0
1056  n = nint(xn_nyq0(mn))/nfp
1057  m = nint(xm_nyq0(mn))
1058  n1 = abs(n)
1059  dmult = mscale(m)*nscale(n1)*tmult
1060  IF (m.eq.0 .or. n.eq.0) dmult = 2*dmult
1061  sgn = sign(1, n)
1062  lk = 0
1063  jlk = js
1064  DO j = 1, ntheta2
1065  DO k = 1, nzeta
1066  lk = lk + 1
1067  tcosi = dmult*(cosmui(j,m)*cosnv(k,n1) +
1068  1 sgn*sinmui(j,m)*sinnv(k,n1))
1069  tsini = dmult*(sinmui(j,m)*cosnv(k,n1) -
1070  1 sgn*cosmui(j,m)*sinnv(k,n1))
1071  bmn(mn) = bmn(mn) + tsini*bsqa(jlk)
1072  gmn(mn) = gmn(mn) + tsini*gsqrta(jlk,0)
1073  bsubumn(mn) = bsubumn(mn) + tsini*bsubua(jlk)
1074  bsubvmn(mn) = bsubvmn(mn) + tsini*bsubva(jlk)
1075  bsubsmn(mn) = bsubsmn(mn) + tcosi*bsubsa(jlk)
1076  bsupumn(mn) = bsupumn(mn) + tsini*bsupua(jlk)
1077  bsupvmn(mn) = bsupvmn(mn) + tsini*bsupva(jlk)
1078 
1079 #ifdef _ANIMEC
1080  pparmn(mn) = pparmn(mn) + tsini*ppara(js,lk)
1081  ppermn(mn) = ppermn(mn) + tsini*pperpa(js,lk)
1082  sigmn(mn) = sigmn(mn) + tsini*sigma_ana(js,lk)
1083  taumn(mn) = taumn(mn) + tsini*tau_ana(js,lk)
1084  pbprmn(mn) = pbprmn(mn) + tsini*pbprima(js,lk)
1085  ppprmn(mn) = ppprmn(mn) + tsini*ppprima(js,lk)
1086  hotdmn(mn) = hotdmn(mn) + tsini*densita(js,lk)
1087 #endif
1088  jlk = jlk+ns
1089  END DO
1090  END DO
1091  END DO mn3
1092 
1093  gmns(:,js) = gmn(:)
1094  bmns(:,js) = bmn(:)
1095  bsubumns(:,js) = bsubumn(:)
1096  bsubvmns(:,js) = bsubvmn(:)
1097  bsubsmnc(:,js) = bsubsmn(:)
1098  bsupumns(:,js) = bsupumn(:)
1099  bsupvmns(:,js) = bsupvmn(:)
1100 #ifdef _ANIMEC
1101  pparmns(:,js) = pparmn(:)
1102  ppermns(:,js) = ppermn(:)
1103  sigmns(:,js) = sigmn(:)
1104  taumns(:,js) = taumn(:)
1105  pbprmns(:,js) = pbprmn(:)
1106  ppprmns(:,js) = ppprmn(:)
1107  hotdmns(:,js) = hotdmn(:)
1108 #endif
1109  END DO radius3
1110 
1111  gmns(:,1) = 0; bmns(:,1) = 0
1112  bsubumns(:,1) = 0
1113  bsubvmns(:,1) = 0
1114  bsubsmnc(:,1) = 2*bsubsmnc(:,2) - bsubsmnc(:,3)
1115  bsupumns(:,1) = 0; bsupvmns(:,1) = 0
1116 #ifdef _ANIMEC
1117  hotdmns(:,1) = 0; pparmns(:,1) = 0; ppermns(:,1) = 0
1118  pbprmns(:,1) = 0; ppprmns(:,1) = 0
1119  sigmns(:,1) = 0; taumns(:,1) = 0
1120 #endif
1121 
1122  IF (lfreeb) THEN !MRC 10-15-15
1123  bsubumns_sur = 0
1124  bsubvmns_sur = 0
1125  bsupumns_sur = 0
1126  bsupvmns_sur = 0
1127 
1128  DO mn = 1, mnmax_nyq0
1129  n = nint(xn_nyq0(mn))/nfp
1130  m = nint(xm_nyq0(mn))
1131  n1 = abs(n)
1132  dmult = mscale(m)*nscale(n1)*tmult
1133  IF (m.eq.0 .or. n.eq.0) dmult = 2*dmult
1134  sgn = sign(1, n)
1135  lk = 0
1136  DO j = 1, ntheta2
1137  DO k = 1, nzeta
1138  lk = lk + 1
1139  tsini = dmult*(sinmui(j,m)*cosnv(k,n1) -
1140  1 sgn*cosmui(j,m)*sinnv(k,n1))
1141  bsubumns_sur(mn) = bsubumns_sur(mn) &
1142  & + tsini*bsubua_sur(lk)
1143  bsubvmns_sur(mn) = bsubvmns_sur(mn) &
1144  & + tsini*bsubva_sur(lk)
1145  bsupumns_sur(mn) = bsupumns_sur(mn) &
1146  & + tsini*bsupua_sur(lk)
1147  bsupvmns_sur(mn) = bsupvmns_sur(mn) &
1148  & + tsini*bsupva_sur(lk)
1149  END DO
1150  END DO
1151  END DO
1152  END IF
1153 
1154  200 CONTINUE
1155 
1156  CALL compute_currents(bsubsmnc, bsubsmns, bsubumnc, bsubumns, &
1157  & bsubvmnc, bsubvmns, &
1158  & xm_nyq0, xn_nyq0, mnmax_nyq0, lasym, ns, &
1159  & currumnc, currvmnc, currumns, currvmns)
1160 
1161 #ifdef _DEBUG
1162  WRITE (333, *) ' JS M*B_S GRAD(B_U) J^V'
1163  DO mn = 1, mnmax_nyq0
1164  WRITE (333,'(2(a,i4))') 'm=', int(xm_nyq0(mn)),
1165  1 ' n=', int(xn_nyq0(mn))/nfp
1166  DO js = 2,ns-1
1167  tmult=-xm_nyq0(mn)*bsubsmns(mn,js) +
1168  1 ohs*(bsubumnc(mn,js+1)-bsubumnc(mn,js))
1169  WRITE (333,'(i6,1p,3e12.4)') js,
1170  1 bsubsmns(mn,js)*xm_nyq0(mn),
1171  2 ohs*(bsubumnc(mn,js+1)-bsubumnc(mn,js)),
1172  3 tmult
1173  END DO
1174  END DO
1175 
1176  WRITE(333,*) version_
1177  IF (lasym) THEN
1178  WRITE(333,2002) 'mn', 'rmnc', 'rmns', 'zmnc', 'zmns', &
1179  & 'lmnc', 'lmns', 'gmnc', 'gmns', &
1180  & 'bmnc', 'bmns', &
1181  & 'bsubumnc', 'bsubumns', &
1182  & 'bsubvmnc', 'bsubvmns', &
1183  & 'bsubsmnc', 'bsubsmns', &
1184  & 'bsupumnc', 'bsupumns', &
1185  & 'bsupvmnc', 'bsupvmns'
1186  ELSE
1187  WRITE(333,2000) 'mn', 'rmnc', 'lmns', 'gmnc', 'bmnc', &
1188  & 'bsubumnc', 'bsubvmnc', &
1189  & 'bsubsmns', &
1190  & 'bsupumnc', 'bsupvmnc'
1191  END IF
1192  DO mn = 1, mnmax
1193  IF (lasym) THEN
1194  WRITE(333,2003) mn, rmnc(mn,ns/2), rmns(mn,ns/2), &
1195  & zmnc(mn,ns/2), zmns(mn,ns/2), &
1196  & lmnc(mn,ns/2), lmns(mn,ns/2), &
1197  & gmnc(mn,ns/2), gmns(mn,ns/2), &
1198  & bmnc(mn,ns/2), bmns(mn,ns/2), &
1199  & bsubumnc(mn,ns/2), bsubumns(mn,ns/2), &
1200  & bsubvmnc(mn,ns/2), bsubvmns(mn,ns/2), &
1201  & bsubsmnc(mn,ns/2), bsubsmns(mn,ns/2), &
1202  & bsupumnc(mn,ns/2), bsupumns(mn,ns/2), &
1203  & bsupvmnc(mn,ns/2), bsupvmns(mn,ns/2)
1204  ELSE
1205  WRITE(333,2001) mn, rmnc(mn,ns/2), lmns(mn,ns/2), &
1206  & gmnc(mn,ns/2), bmnc(mn,ns/2), &
1207  & bsubumnc(mn,ns/2), bsubvmnc(mn,ns/2), &
1208  & bsubsmns(mn,ns/2), &
1209  & bsupumnc(mn,ns/2), bsupvmnc(mn,ns/2)
1210  END IF
1211  END DO
1212 2000 FORMAT(a2,10(2x,a12))
1213 2001 FORMAT(i2,10(2x,e12.5))
1214 2002 FORMAT(a2,20(2x,a12))
1215 2003 FORMAT(i2,20(2x,es12.5))
1216 #endif
1217 !
1218 ! WRITE OUT ARRAYS
1219 !
1220 #ifdef NETCDF
1221  CALL cdf_write(nwout, vn_racc, raxis_cc(0:ntor))
1222  CALL cdf_write(nwout, vn_zacs, zaxis_cs(0:ntor))
1223  CALL cdf_write(nwout, vn_rmnc, rmnc)
1224  CALL cdf_write(nwout, vn_zmns, zmns)
1225  CALL cdf_write(nwout, vn_lmns, lmns)
1226  CALL cdf_write(nwout, vn_gmnc, gmnc) !Half mesh
1227  CALL cdf_write(nwout, vn_bmnc, bmnc) !Half mesh
1228  CALL cdf_write(nwout, vn_bsubumnc, bsubumnc) !Half mesh
1229  CALL cdf_write(nwout, vn_bsubvmnc, bsubvmnc) !Half mesh
1230  CALL cdf_write(nwout, vn_bsubsmns, bsubsmns) !Full mesh
1231 
1232  CALL cdf_write(nwout, vn_currumnc, currumnc) !MRK 8-12-16
1233  CALL cdf_write(nwout, vn_currvmnc, currvmnc)
1234 
1235 ! GET RID OF THESE EVENTUALLY: DON'T NEED THEM (can express in terms of lambdas)
1236  CALL cdf_write(nwout, vn_bsupumnc, bsupumnc)
1237  CALL cdf_write(nwout, vn_bsupvmnc, bsupvmnc)
1238 
1239  IF (lfreeb) THEN !MRC 10-15-15
1240  CALL cdf_write(nwout, vn_bsubumnc_sur, bsubumnc_sur)
1241  CALL cdf_write(nwout, vn_bsubvmnc_sur, bsubvmnc_sur)
1242  CALL cdf_write(nwout, vn_bsupumnc_sur, bsupumnc_sur)
1243  CALL cdf_write(nwout, vn_bsupvmnc_sur, bsupvmnc_sur)
1244  END IF
1245 
1246 ! FULL-MESH quantities
1247 ! NOTE: jdotb is in units_of_A (1/mu0 incorporated in jxbforce...)
1248 ! prior to version 6.00, this was output in internal VMEC units...
1249 
1250  j = SIZE(am)-1
1251  CALL cdf_write(nwout, vn_am, am(0:j))
1252  j = SIZE(ac)-1
1253  CALL cdf_write(nwout, vn_ac, ac(0:j))
1254  j = SIZE(ai)-1
1255  CALL cdf_write(nwout, vn_ai, ai(0:j))
1256 
1257  j = SIZE(am_aux_s)
1258  CALL cdf_write(nwout, vn_am_aux_s, am_aux_s(1:j))
1259  j = SIZE(am_aux_f)
1260  CALL cdf_write(nwout, vn_am_aux_f, am_aux_f(1:j))
1261  j = SIZE(ac_aux_s)
1262  CALL cdf_write(nwout, vn_ac_aux_s, ac_aux_s(1:j))
1263  j = SIZE(ac_aux_f)
1264  CALL cdf_write(nwout, vn_ac_aux_f, ac_aux_f(1:j))
1265  j = SIZE(ai_aux_s)
1266  CALL cdf_write(nwout, vn_ai_aux_s, ai_aux_s(1:j))
1267  j = SIZE(ai_aux_f)
1268  CALL cdf_write(nwout, vn_ai_aux_f, ai_aux_f(1:j))
1269 
1270  CALL cdf_write(nwout, vn_iotaf, iotaf(1:ns))
1271  CALL cdf_write(nwout, vn_qfact, qfact(1:ns))
1272  CALL cdf_write(nwout, vn_presf, presf/mu0)
1273  CALL cdf_write(nwout, vn_phi, phi)
1274  CALL cdf_write(nwout, vn_phipf, twopi*signgs*phipf)
1275  CALL cdf_write(nwout, vn_chi, chi)
1276  CALL cdf_write(nwout, vn_chipf, twopi*signgs*chipf)
1277  CALL cdf_write(nwout, vn_jcuru, jcuru/mu0)
1278  CALL cdf_write(nwout, vn_jcurv, jcurv/mu0)
1279  CALL cdf_write(nwout, vn_jdotb, jdotb)
1280  CALL cdf_write(nwout, vn_bdotb, bdotb)
1281  CALL cdf_write(nwout, vn_bgrv, bdotgradv)
1282 
1283 ! HALF-MESH quantities
1284  iotas(1) = 0; mass(1) = 0; pres(1) = 0; phip(1) = 0;
1285  buco(1) = 0; bvco(1) = 0; vp(1) = 0; overr(1) = 0; specw(1) = 1
1286  beta_vol(1) = 0
1287  CALL cdf_write(nwout, vn_iotah, iotas(1:ns))
1288  CALL cdf_write(nwout, vn_mass, mass/mu0)
1289  CALL cdf_write(nwout, vn_presh, pres(1:ns)/mu0)
1290  CALL cdf_write(nwout, vn_betah, beta_vol)
1291  CALL cdf_write(nwout, vn_buco, buco)
1292  CALL cdf_write(nwout, vn_bvco, bvco)
1293  CALL cdf_write(nwout, vn_vp, vp(1:ns))
1294  CALL cdf_write(nwout, vn_specw, specw)
1295  CALL cdf_write(nwout, vn_phip, phips(1:ns))
1296  CALL cdf_write(nwout, vn_overr, overr(1:ns))
1297 
1298 ! MERCIER_CRITERION
1299  CALL cdf_write(nwout, vn_merc, dmerc)
1300  CALL cdf_write(nwout, vn_mshear, dshear)
1301  CALL cdf_write(nwout, vn_mwell, dwell)
1302  CALL cdf_write(nwout, vn_mcurr, dcurr)
1303  CALL cdf_write(nwout, vn_mgeo, dgeod)
1304  CALL cdf_write(nwout, vn_equif, equif)
1305 
1306  CALL cdf_write(nwout, vn_fsq, fsqt(1:nstore_seq))
1307  CALL cdf_write(nwout, vn_wdot, wdot(1:nstore_seq))
1308 
1309 !-----------------------------------------------
1310 ! DATA AND MSE FITS : HAVE NOT CONVERTED TO NETCDF
1311 ! SINCE THIS WILL BE REPLACED SOON
1312 !-----------------------------------------------
1313 ! IF (.not.lrecon) GOTO 925
1314 ! 925 CONTINUE
1315 
1316 #endif
1317  IF (lwouttxt) THEN
1318  DO js = 1, ns
1319  WRITE(nwout2, *) "JS: ", js
1320  mn1_out: DO mn = 1, mnmax
1321  IF (js .eq. 1) THEN
1322  WRITE (nwout2, *) nint(xm(mn)), nint(xn(mn))
1323  END IF
1324 
1325  WRITE (nwout2, *) rmnc(mn,js), zmns(mn,js), lmns(mn,js)
1326  IF (lasym) THEN
1327  WRITE (nwout2, *)rmns(mn,js),zmnc(mn,js),lmnc(mn,js)
1328  ENDIF
1329  END DO mn1_out
1330 
1331  mn2_out: DO mn = 1, mnmax_nyq0
1332  IF (js .eq. 1) THEN
1333  WRITE (nwout2, *) nint(xm_nyq0(mn)), nint(xn_nyq0(mn))
1334  END IF
1335  WRITE (nwout2, *) bmnc(mn,js), gmnc(mn,js),
1336  1 bsubumnc(mn,js), bsubvmnc(mn,js), bsubsmns(mn,js),
1337  2 bsupumnc(mn,js), bsupvmnc(mn,js)
1338 #ifdef _ANIMEC
1339  3 ,pparmnc(mn,js), ppermnc(mn,js), hotdmnc(mn,js),
1340  4 pbprmnc(mn,js), ppprmnc(mn,js), sigmnc(mn,js),
1341  5 taumnc(mn,js)
1342 #endif
1343  IF (lasym) THEN
1344  WRITE (nwout2, *) bmns(mn,js), gmns(mn,js),
1345  1 bsubumns(mn,js), bsubvmns(mn,js), bsubsmnc(mn,js),
1346  2 bsupumns(mn,js), bsupvmns(mn,js)
1347 #ifdef _ANIMEC
1348  3 ,pparmns(mn,js), ppermns(mn,js), hotdmns(mn,js),
1349  4 pbprmns(mn,js), ppprmns(mn,js), sigmns(mn,js),
1350  5 taumns(mn,js)
1351 #endif
1352  ENDIF
1353  END DO mn2_out
1354  END DO
1355 
1356 !
1357 ! HALF-MESH QUANTITIES (except phi, jcuru, jcurv which are FULL MESH - computed in eqfor)
1358 ! NOTE: jcuru, jcurv are local current densities, NOT integrated in s and normed to twopi
1359 ! NOTE: In version <= 6.00, mass, press are written out in INTERNAL units
1360 ! and should be multiplied by 1/mu0 to transform to pascals. In version > 6.00,
1361 ! the pressure, mass are in correct (physical) units
1362 !
1363 
1364 ! NOTE: phipf has a twopi * signgs factor compared with phips...
1365 
1366 
1367  WRITE (nwout2, *) (iotaf(js), presf(js)/mu0,
1368  1 twopi*signgs*phipf(js),
1369  2 phi(js), jcuru(js)/mu0, jcurv(js)/mu0, js=1,ns)
1370  WRITE (nwout2, *) (iotas(js), mass(js)/mu0, pres(js)/mu0,
1371  1 beta_vol(js), phip(js), buco(js), bvco(js), vp(js),
1372  2 overr(js), specw(js),js=2,ns)
1373 !-----------------------------------------------
1374 
1375  WRITE (nwout2, *) aspect, betatot, betapol, betator, betaxis,
1376  1 b0
1377 
1378 !-----------------------------------------------
1379 ! New output added to version 6.20
1380 !-----------------------------------------------
1381  WRITE (nwout2, *) nint(signgs)
1382  WRITE (nwout2, '(a)') input_extension
1383  WRITE (nwout2, *) ionlarmor, volavgb, rbtor0, rbtor, ctor/mu0,
1384  1 aminor_p, rmajor_p, volume_p
1385 !-----------------------------------------------
1386 ! MERCIER CRITERION
1387 !-----------------------------------------------
1388  WRITE (nwout2, *) (dmerc(js), dshear(js), dwell(js), dcurr(js),
1389  1 dgeod(js), equif(js), js=2,ns-1)
1390 
1391  IF (nextcur.gt.0) THEN
1392  WRITE (nwout2, *) (extcur(i),i=1,nextcur)
1393  lcurr = ALLOCATED(curlabel) .and. lfreeb
1394  WRITE (nwout2, *) lcurr
1395  IF (lcurr) WRITE (nwout2, *) (curlabel(i),i=1,nextcur)
1396  ENDIF
1397 
1398 !-----------------------------------------------
1399 ! NOTE: jdotb is in units of A (1/mu0 incorporated in jxbforce...)
1400 ! prior to version 6.00, this was output in internal VMEC units...
1401 !-----------------------------------------------
1402  WRITE (nwout2, *) (fsqt(i),wdot(i),i=1,nstore_seq)
1403  WRITE (nwout2, *) (jdotb(js),bdotgradv(js),bdotb(js),js=1,ns)
1404 
1405 !-----------------------------------------------
1406 ! Modification to obtain old fort.8 file (depracated)
1407 ! Write out only the stellarator symmetric parts
1408 ! Only kept for old codes. (J. Geiger)
1409 !-----------------------------------------------
1410  IF (loldout) THEN
1411  WRITE (nfort8, '(f10.3,7i6)')
1412  1 gamma, nfp, ns, mpol, ntor, mnmax, itfsq, iter2/100+1
1413  DO js = 1, ns
1414  mn = 0
1415  DO m = 0, mpol1
1416  nmin0 = -ntor
1417  IF (m .eq. 0) nmin0 = 0
1418  DO n = nmin0, ntor
1419  mn = mn + 1
1420  IF (js .eq. 1)
1421  1 WRITE (nfort8,'(2i10)') nint(xm(mn)),nint(xn(mn))
1422  WRITE (nfort8,'(5e20.13)')
1423  1 rmnc(mn,js),zmns(mn,js),lmns(mn,js),
1424  2 bmnc(mn,js),gmnc(mn,js),
1425  3 bsubumnc(mn,js),bsubvmnc(mn,js),bsubsmns(mn,js),
1426  4 bsupumnc(mn,js),bsupvmnc(mn,js)
1427  END DO
1428  END DO
1429  END DO
1430  WRITE (nfort8,'(5e20.13)') (iotas(js),mass(js),pres(js),
1431  1 phips(js),buco(js),bvco(js),phi(js),vp(js),
1432  2 jcuru(js)/mu0,jcurv(js)/mu0,specw(js),js=2,ns)
1433  WRITE (nfort8,'(5e20.13)') (fsqt(i),wdot(i),i=1,100)
1434  CLOSE(nfort8) !last write to nfort8
1435  END IF
1436 !-----------------------------------------------
1437 ! Write diagno file (J.Geiger)
1438 !-----------------------------------------------
1439  IF ((.not.lasym).and. ldiagno .and.lfreeb) THEN
1440  open(unit=21,file='diagno_input.data',status='unknown',
1441  1 action='write')
1442  write(21,'(a8)') "vmec2000"
1443  write(21,*) "nfp mpol ntor"
1444  write(21,*) nfp, mpol, ntor
1445  write(21,*) "rmnc"
1446  write(21,'(1p,e24.16)') (rmnc(mn,ns),mn=1,mnmax)
1447  write(21,*) "zmns"
1448  write(21,'(1p,e24.16)') (zmns(mn,ns),mn=1,mnmax)
1449  write(21,*) "potsin"
1450  DO mn = 1, mnpd
1451  write(21,'(1p,e24.16)') potvac(mn)
1452  END DO
1453  write(21,*) "phiedge"
1454  write(21,*) phiedge
1455  write(21,*) "nextcur"
1456  write(21,*) nextcur
1457  write(21,*) "external currents"
1458  write(21,*) extcur(1:nextcur)
1459  write(21,*) "plasma current"
1460  write(21,*) ctor
1461  write(21,*) "plasma current filament fc R"
1462  write(21,*) rmnc(1:ntor+1,1)
1463  write(21,*) "plasma current filament fc z"
1464  write(21,*) zmns(1:ntor+1,1)
1465  close(unit=21)
1466  END IF
1467 !-----------------------------------------------
1468 ! for diagno version 1.5 written by Sam Lazerson (SAL) start
1469 !-----------------------------------------------
1470  IF(ldiagno)THEN
1471  IF(lfreeb .and. (.not.lasym))THEN
1472  nfort = 21
1473  fort_file = 'diagno1.5_in.'//input_extension
1474  call safe_open(nfort,iwout0,fort_file,'replace',
1475  1 'formatted')
1476  if (iwout0 .ne. 0)
1477  1 stop 'Error writing diagno_in. file in VMEC WROUT'
1478 
1479  write(21,'(a)') "vmec2000_B"
1480  write(21,*) "nfp mpol ntor"
1481  write(21,*) nfp, mpol, ntor
1482  write(21,*) "rmnc"
1483  write(21,'(1p,e24.16)') (rmnc(mn,ns),mn=1,mnmax)
1484  write(21,*) "zmns"
1485  write(21,'(1p,e24.16)') (zmns(mn,ns),mn=1,mnmax)
1486 
1487  write(21,*) "potsin"
1488  DO i = 1, mnpd
1489  write(21,'(1p,e24.16)') potvac(i)
1490  END DO
1491 
1492 !----- Added by SAL 11/2010
1493  write(21,*) "bsupu"
1494  js=ns
1495  js2=ns-1
1496  do m = 0, mpol1
1497  nmin0 = -ntor
1498  if (m .eq. 0) nmin0 = 0
1499  do n = nmin0, ntor
1500  dmult = two/(mscale(m)*nscale(abs(n)))
1501  if (m .eq. 0 .and. n .eq. 0) dmult = p5*dmult
1502  n1 = abs(n)
1503  isgn = sign(1, n)
1504  d_bsupumn = 0
1505  do j = 1, ntheta2
1506  do k = 1, nzeta
1507  lk = k + nzeta*(j - 1)
1508  tcosi = dmult*(cosmui(j,m)*cosnv(k,n1) +
1509  1 isgn*sinmui(j,m)*sinnv(k,n1))
1510  d_bsupumn = d_bsupumn + 1.5*tcosi*bsupu(js,lk)
1511  1 - 0.5*tcosi*bsupu(js2,lk)
1512  end do
1513  end do
1514  write (21,'(1p,e24.16)') d_bsupumn
1515  end do
1516  end do
1517 !----- Added by SAL 11/2010
1518  write(21,*) "bsupv"
1519  js=ns
1520  js2=ns-1
1521  do m = 0, mpol1
1522  nmin0 = -ntor
1523  if (m .eq. 0) nmin0 = 0
1524  do n = nmin0, ntor
1525  dmult = two/(mscale(m)*nscale(abs(n)))
1526  if (m .eq. 0 .and. n .eq. 0) dmult = p5*dmult
1527  n1 = abs(n)
1528  isgn = sign(1, n)
1529  d_bsupvmn = 0
1530  do j = 1, ntheta2
1531  do k = 1, nzeta
1532  lk = k + nzeta*(j - 1)
1533  tcosi = dmult*(cosmui(j,m)*cosnv(k,n1) +
1534  1 isgn*sinmui(j,m)*sinnv(k,n1))
1535  d_bsupvmn = d_bsupvmn + 1.5*tcosi*bsupv(js,lk)
1536  1 - 0.5*tcosi*bsupv(js2,lk)
1537  end do
1538  end do
1539  write (21,'(1p,e24.16)') d_bsupvmn
1540  end do
1541  end do
1542 
1543  write(21,*) "phiedge"
1544  write(21,*) phiedge
1545  write(21,*) "nextcur"
1546  write(21,*) nextcur
1547  write(21,*) "external currents"
1548  write(21,*) extcur(1:nextcur)
1549 
1550  write(21,*) "plasma current"
1551  write(21,*) ctor
1552  write(21,*) "plasma current filament fc R"
1553  write(21,*) rmnc(1:ntor+1,1)
1554  write(21,*) "plasma current filament fc z"
1555  write(21,*) zmns(1:ntor+1,1)
1556 
1557  close(unit=21)
1558  ELSE
1559  write(6,*)"Diagno-file request not completed!"
1560  write(6,*)"VMEC2000 not running in free-boundary mode!"
1561  write(6,*)"-or- LASYM = .true. !"
1562  write(6,*)"LASYM = ",lasym
1563  write(6,*)"LFREEB = ",lfreeb
1564  write(6,*)"Check mgrid-file and namelist!"
1565  ENDIF
1566  ENDIF !added for diagno version 1.5 end
1567 
1568  ENDIF
1569 
1570  950 CONTINUE
1571 
1572  IF (lwouttxt) CLOSE (unit=nwout2) !J Geiger: Close only if open, i.e. lwouttxt==true
1573 !--------------------DEC$ ENDIF
1574  IF (.not. lwrite) GOTO 970 ! J Geiger: in case lwouttxt is not true and netcdf-write is finished
1575 #ifdef NETCDF
1576  IF (lasym) THEN
1577  CALL cdf_write(nwout, vn_racs, raxis_cs(0:ntor))
1578  CALL cdf_write(nwout, vn_zacc, zaxis_cc(0:ntor))
1579  CALL cdf_write(nwout, vn_rmns, rmns)
1580  CALL cdf_write(nwout, vn_zmnc, zmnc)
1581  CALL cdf_write(nwout, vn_lmnc, lmnc)
1582  CALL cdf_write(nwout, vn_gmns, gmns)
1583  CALL cdf_write(nwout, vn_bmns, bmns)
1584  CALL cdf_write(nwout, vn_bsubumns, bsubumns)
1585  CALL cdf_write(nwout, vn_bsubvmns, bsubvmns)
1586  CALL cdf_write(nwout, vn_bsubsmnc, bsubsmnc)
1587 
1588  CALL cdf_write(nwout, vn_currumns, currumns) !MRC 8-12-16
1589  CALL cdf_write(nwout, vn_currvmns, currvmns)
1590 
1591 ! GET RID OF THESE EVENTUALLY: DON'T NEED THEM
1592  CALL cdf_write(nwout, vn_bsupumns, bsupumns)
1593  CALL cdf_write(nwout, vn_bsupvmns, bsupvmns)
1594 
1595  IF (lfreeb) THEN !MRC 10-15-15
1596  CALL cdf_write(nwout, vn_bsubumns_sur, bsubumns_sur)
1597  CALL cdf_write(nwout, vn_bsubvmns_sur, bsubvmns_sur)
1598  CALL cdf_write(nwout, vn_bsupumns_sur, bsupumns_sur)
1599  CALL cdf_write(nwout, vn_bsupvmns_sur, bsupvmns_sur)
1600  END IF
1601  END IF
1602 #endif
1603  970 CONTINUE ! J Geiger: need to keep label 970 out of NETCDF defines.
1604 
1605 #ifdef NETCDF
1606  CALL cdf_close(nwout)
1607 #endif
1608 !
1609 ! RESTORE nyq ENDPOINT VALUES
1610 !
1611  IF (lnyquist) THEN
1612  IF (mnyq .ne. 0) cosmui(:,mnyq) = 2*cosmui(:,mnyq)
1613  IF (nnyq .ne. 0) cosnv(:,nnyq) = 2*cosnv(:,nnyq)
1614  END IF
1615 
1616 !
1617 ! DEALLOCATIONS ! J Geiger: these have been moved downwards.
1618 !
1619  IF (ALLOCATED(gmnc)) DEALLOCATE(gmnc, bmnc, bsubumnc, bsubvmnc,
1620  1 bsubsmns, bsupumnc, bsupvmnc
1621 #ifdef _ANIMEC
1622  2 ,sigmnc,taumnc,pparmnc,ppermnc,pbprmnc,ppprmnc,hotdmnc
1623 #endif
1624  3 )
1625  IF (ALLOCATED(gmns)) DEALLOCATE(gmns, bmns, bsubumns, bsubvmns,
1626  1 bsubsmnc, bsupumns, bsupvmns
1627 #ifdef _ANIMEC
1628  2 ,sigmns,taumns,pparmns,ppermns,pbprmns,ppprmns,hotdmns
1629 #endif
1630  3 )
1631 #ifdef _ANIMEC
1632  IF (ALLOCATED(tau_ana)) DEALLOCATE(sigma_ana, tau_ana, ppara,
1633  1 pperpa, pbprima, ppprima, densita)
1634 #endif
1635 ! J Geiger: check also for allocation.
1636  IF (ALLOCATED(gmn)) DEALLOCATE (gmn, bmn, bsubumn, bsubvmn,
1637  1 bsubsmn, bsupumn, bsupvmn,
1638 #ifdef _ANIMEC
1639  2 sigmn, taumn, pparmn, ppermn, pbprmn, ppprmn,
1640  3 hotdmn,
1641 #endif
1642  4 stat=istat)
1643 
1644  IF (ALLOCATED(bsubumnc_sur)) THEN
1645  DEALLOCATE(bsubumnc_sur, bsubvmnc_sur)
1646  DEALLOCATE(bsupumnc_sur, bsupvmnc_sur)
1647  END IF
1648  IF (ALLOCATED(bsubumns_sur)) THEN
1649  DEALLOCATE(bsubumns_sur, bsubvmns_sur)
1650  DEALLOCATE(bsupumns_sur, bsupvmns_sur)
1651  END IF
1652  IF (ALLOCATED(bsubua_sur)) THEN
1653  DEALLOCATE(bsubua_sur, bsubva_sur)
1654  DEALLOCATE(bsupua_sur, bsupva_sur)
1655  END IF
1656 
1657 !-----------------------------------------------
1658 ! FREE BOUNDARY DATA
1659 !-----------------------------------------------
1660  IF (lwrite )
1661  1 CALL freeb_data(rmnc1, zmns1, rmns1, zmnc1, bmodmn, bmodmn1)
1662  1000 CONTINUE
1663 
1664  rzl_array = 0
1665 
1666  CALL second0 (twoutoff)
1667  timer(twout) = timer(twout) + twoutoff - twouton
1668  fo_wrout_time = timer(twout)
1669 
1670  END SUBROUTINE wrout
Generated on Thu Mar 5 2020 15:49:27 for V3FIT by   doxygen 1.8.17