V3FIT
readin.f
1  SUBROUTINE readin(input_file, iseq_count, ier_flag, lscreen)
2  USE vmec_main
3  USE vmec_params
4  USE vacmod
5  USE vspline
6  USE timer_sub
7  USE mgrid_mod, ONLY: nextcur, curlabel, nfper0, read_mgrid
8  USE init_geometry
9  USE parallel_include_module, ONLY: grank, mgrid_file_read_time,
10  & lprecond
11  USE parallel_vmec_module, ONLY: runvmec_comm_world
12  IMPLICIT NONE
13 C-----------------------------------------------
14 C D u m m y A r g u m e n t s
15 C-----------------------------------------------
16  INTEGER :: iseq_count, ier_flag
17  LOGICAL :: lscreen
18  CHARACTER(LEN=*) :: input_file
19 C-----------------------------------------------
20 C L o c a l V a r i a b l e s
21 C-----------------------------------------------
22  INTEGER :: iexit, ipoint, n, iunit, ier_flag_init,
23  & i, ni, m, nsmin, igrid, mj, isgn, ioff, joff,
24  & NonZeroLen
25  REAL(dp), DIMENSION(:,:), POINTER ::
26  & rbcc, rbss, rbcs, rbsc, zbcs, zbsc, zbcc, zbss
27  REAL(dp) :: rtest, ztest, tzc, trc, delta
28  REAL(dp), ALLOCATABLE :: temp(:)
29  CHARACTER(LEN=100) :: line, line2
30  CHARACTER(LEN=1) :: ch1, ch2
31  LOGICAL :: lwrite
32 C-----------------------------------------------
33 !
34 ! LOCAL VARIABLES
35 !
36 ! rbcc,rbss,rbcs,rbsc
37 ! boundary Fourier coefficient arrays for R (of cosu*cosv, etc)
38 ! zbcc,zbss,zbcs,zbsc
39 ! boundary Fourier coefficient arrays for Z
40 !
41 ! XCC*COS(MU)COS(NV), XCS*COS(MU)SIN(NV), ETC
42 !
43 ! STACKING ORDER DEPENDS ON LASYM AND LTHREED. EACH COMPONENT XCC, XSS, XSC, XCS
44 ! HAS SIZE = mns. (PHIFAC, MSE TAKE UP 1 INDEX EACH AT END OF ARRAY)
45 !
46 ! LTHREED=F, LTHREED=F, LTHREED=T, LTHREED=T
47 ! LASYM=F LASYM=T LASYM=F LASYM=T
48 !
49 ! rmncc rmncc rmncc rmncc
50 ! zmnsc rmnsc rmnss rmnss
51 ! lmnsc zmnsc zmnsc rmnsc
52 ! zmncc zmncs rmncs
53 ! lmnsc lmnsc zmnsc
54 ! lmncc lmncs zmncs
55 ! zmncc
56 ! zmnss
57 ! lmnsc
58 ! lmncs
59 ! lmncc
60 ! lmnss
61 !
62 !
63 ! STANDARD INPUT DATA AND RECOMMENDED VALUES
64 !
65 ! Plasma parameters (MKS units)
66 ! ai: expansion coefficients for iota (power series in s) used when ncurr=0
67 ! Interpretation changes with piota_type
68 ! am: mass or pressure (gamma=0) expansion coefficients (series in s)
69 ! in MKS units [NWT/M**2]
70 ! Interpretation changes with pmass_type
71 ! ac: expansion coefficients for the normalized (pcurr(s=1) = 1)
72 ! radial derivative of the flux-averaged toroidal current density
73 ! (power series in s) used when ncurr=1
74 ! Interpretation changes with pcurr_type
75 ! ai_aux_s: Auxiliary array for iota profile. Used for splines, s values
76 ! ai_aux_f: Auxiliary array for iota profile. Used for splines, function values
77 ! am_aux_s: Auxiliary array for mass profile. Used for splines, s values
78 ! am_aux_f: Auxiliary array for mass profile. Used for splines, function values
79 ! ac_aux_s: Auxiliary array for current profile. Used for splines, s values
80 ! ac_aux_f: Auxiliary array for current profile. Used for splines, function values
81 ! curtor: value of toroidal current [A]. Used if ncurr = 1 to specify
82 ! current profile, or IF in data reconstruction mode.
83 ! phiedge: toroidal flux enclosed by plasma at edge (in Wb)
84 ! extcur: array of currents in each external current group. Used to
85 ! multiply Green''s function for fields and loops read in from
86 ! MGRID file. Should use real current units (A).
87 ! gamma: value of compressibility index (gamma=0 => pressure prescribed)
88 ! nfp: number of toroidal field periods ( =1 for Tokamak)
89 ! rbc: boundary coefficients of COS(m*theta-n*zeta) for R [m]
90 ! zbs: boundary coefficients of SIN(m*theta-n*zeta) for Z [m]
91 ! rbs: boundary coefficients of SIN(m*theta-n*zeta) for R [m]
92 ! zbc: boundary coefficients of COS(m*theta-n*zeta) for Z [m]
93 !
94 !
95 ! Numerical and logical control parameters
96 ! ncurr: flux conserving (=0) or prescribed toroidal current (=1)
97 ! ns_array: array of radial mesh sizes to be used in multigrid sequence
98 ! nvacskip: number of iteration steps between accurate calculation of vacuum
99 ! response; use fast interpolation scheme in between
100 ! pres_scale: factor used to scale pressure profile (default value = 1)
101 ! useful so user can fix profile and change beta without having to change
102 ! all AM coefficients separately
103 ! tcon0: weight factor for constraint force (=1 by DEFAULT)
104 ! lasym: =T, run in asymmetric mode; =F, run in stellarator symmetry mode
105 ! lfreeb: =T, run in free boundary mode if mgrid_file exists
106 ! lforbal: =T, use non-variational forces to ensure <EQUIF> = 0;
107 ! =F, use variational form of forces, <EQUIF> ~ 0
108 !
109 ! Convergence control parameters
110 ! ftol_array: array of value of residual(s) at which each multigrid
111 ! iteration ends
112 ! niter_array: array of number of iterations (used to terminate run) at
113 ! each multigrid iteration
114 ! nstep: number of timesteps between printouts on screen
115 ! nvacskip: iterations skipped between full update of vacuum solution
116 !
117 ! Preconditioner control parameters (added 8/30/04)
118 ! precon_type: specifies type of 2D preconditioner to use ('default', diagonal in m,n,
119 ! tri-diagonal in s; 'conjugate-gradient', block tri-di, evolve using
120 ! cg method; 'gmres', block tri-di, generalized minimal residual method;
121 ! 'tfqmr', block tri-di, transpose-free quasi minimum residual
122 ! prec2d_threshold:
123 ! value of preconditioned force residuals at which block (2d) tri-di
124 ! solver is turned on, if requested via type_prec2d
125 !
126 ! Character parameters
127 ! mgrid_file: full path for vacuum Green''s function data
128 ! pcurr_type: Specifies parameterization type of pcurr function
129 ! 'power_series' - I'(s)=Sum[ ac(j) s ** j] - Default
130 ! 'gauss_trunc' - I'(s)=ac(0) (exp(-(s/ac(1)) ** 2) -
131 ! exp(-(1/ac(1)) ** 2))
132 ! others - see function pcurr
133 ! piota_type: Specifies parameterization type of piota function
134 ! 'power_series' - p(s)=Sum[ am(j) s ** j] - Default
135 ! others - see function piota
136 ! pmass_type: Specifies parameterization type of pmass function
137 ! 'power_series' - p(s)=Sum[ am(j) s ** j] - Default
138 ! 'gauss_trunc' - p(s)=am(0) (exp(-(s/am(1)) ** 2) -
139 ! exp(-(1/am(1)) ** 2))
140 ! others - see function pmass
141 
142 ! Equilibrium reconstruction parameters
143 ! phifac: factor scaling toroidal flux to match apres or limiter
144 ! datastark: pitch angle data from stark measurement
145 ! datathom: pressure data from Thompson, CHEERS (Pa)
146 ! imatch_ = 1 (default),match value of PHIEDGE in input file
147 ! phiedge: = 0, USE pressure profile width to determine PHIEDGE
148 ! = 2, USE LIMPOS data (in mgrid file) to find PHIEDGE
149 ! = 3, USE Ip to find PHIEDGE (fixed-boundary only)
150 ! imse: number of Motional Stark effect data points
151 ! >0, USE mse data to find iota; <=0, fixed iota profile ai
152 ! itse: number of pressure profile data points
153 ! = 0, no thompson scattering data to READ
154 ! isnodes: number of iota spline points (computed internally unless specified explicitly)
155 ! ipnodes: number of pressure spline points (computed internally unless specified explicitly)
156 ! lpofr: LOGICAL variable. =.true. IF pressure data are
157 ! prescribed in REAL space. =.false. IF data in flux space.
158 ! pknots: array of pressure knot values in SQRT(s) space
159 ! sknots: array of iota knot values in SQRT(s) space
160 ! tensp: spline tension for pressure profile
161 !
162 ! tensi: spline tension for iota
163 ! tensi2: vbl spline tension for iota
164 ! fpolyi: vbl spline tension form factor (note: IF tensi!=tensi2
165 ! THEN tension(i-th point) = tensi+(tensi2-tensi)*(i/n-1))**fpolyi
166 ! - - - - - - - - - - - - - - - - - -
167 ! mseangle_ uniform EXPerimental offset of MSE data
168 ! offset: (calibration offset) ... PLUS ...
169 ! mseangle_ multiplier on mseprof offset array
170 ! offsetM: (calibration offset)
171 ! mseprof: offset array from NAMELIST MSEPROFIL
172 ! so that the total offset on the i-th MSE data point is
173 ! taken to be
174 ! = mseangle_offset+mseangle_offsetM*mseprof(i)
175 ! - - - - - - - - - - - - - - - - - -
176 ! pres_offset: uniform arbitrary radial offset of pressure data
177 ! presfac: number by which Thomson scattering data is scaled
178 ! to get actual pressure
179 ! phidiam: diamagnetic toroidal flux (Wb)
180 ! dsiobt: measured flux loop signals corresponding to the
181 ! combination of signals in iconnect array
182 ! indxflx: array giving INDEX of flux measurement in iconnect array
183 ! indxbfld: array giving INDEX of bfield measurement used in matching
184 ! nobd: number of connected flux loop measurements
185 ! nobser: number of individual flux loop positions
186 ! nbsets: number of B-coil sets defined in mgrid file
187 ! nbcoils(n): number of bfield coils in each set defined in mgrid file
188 ! nbcoilsn: total number of bfield coils defined in mgrid file
189 ! bbc(m,n): measured magnetic field at rbcoil(m,n),zbcoil(m,n) at
190 ! the orientation br*COS(abcoil) + bz*SIN(abcoil)
191 ! rbcoil(m,n): R position of the m-th coil in the n-th set from mgrid file
192 ! zbcoil(m,n): Z position of the m-th coil in the n-th set from mgrid file
193 ! abcoil(m,n): orientation (surface normal wrt R axis; in radians)
194 ! of the m-th coil in the n-th set from mgrid file.
195 ! nflxs: number of flux loop measurements used in matching
196 ! nbfld(n): number of selected EXTERNAL bfield measurements in set n from nml file
197 ! nbfldn: total number of EXTERNAL bfield measurements used in matching
198 ! - - - - - - - - - - - - - - - - - -
199 ! NOTE: FOR STANDARD DEVIATIONS (sigma''s) < 0, INTERPRET
200 ! AS PERCENT OF RESPECTIVE MEASUREMENT
201 ! sigma_thom: standard deviation (Pa) for pressure profile data
202 ! sigma_stark: standard deviation (degrees) in MSE data
203 ! sigma_flux: standard deviaton (Wb) for EXTERNAL poloidal flux data
204 ! sigma_b: standard deviation (T) for EXTERNAL magnetic field data
205 !sigma_current: standard deviation (A) in toroidal current
206 !sigma_delphid: standard deviation (Wb) for diamagnetic match
207 !
208 !
209 ! THE (ABSOLUTE) CHI-SQ ERROR IS DEFINED AS FOLLOWS:
210 !
211 ! 2
212 ! CHI = SUM [ EQ(K,IOTA,PRESSURE) - DATA(K) ] ** 2
213 ! (K) -----------------------------------
214 ! SIGMA(K)**2
215 !
216 ! HERE, SIGMA IS THE STANDARD DEVIATION OF THE MEASURED DATA, AND
217 ! EQ(IOTA,PRESSURE) IS THE EQUILIBRIUM EXPRESSION FOR THE DATA TO BE
218 ! MATCHED:
219 !
220 ! EQ(I) = SUM [ W(I,J)*X(J) ]
221 ! (J)
222 !
223 ! WHERE W(I,J) ARE THE (LINEAR) MATRIX ELEMENTS AND X(J) REPRESENT
224 ! THE KNOT VALUES OF IOTA (AND/OR PRESSURE). THE RESULTING LEAST-SQUARES
225 ! MATRIX ELEMENTS AND DATA ARRAY CAN BE EXPRESSED AS FOLLOWS:
226 !
227 ! ALSQ(I,J) = SUM [ W(K,I) * W(K,J) / SIGMA(K) ** 2]
228 ! (K)
229 !
230 ! BLSQ(I) = SUM [ W(K,I) * DATA(K)/ SIGMA(K) ** 2]
231 ! (K)
232 !
233 ! THEREFORE, INTERNALLY IT IS CONVENIENT TO WORK WITH THE 'SCALED'
234 ! W'(K,I) = W(K,I)/SIGMA(K) AND DATA'(K) = DATA(K)/SIGMA(K)
235 !
236 ! ****! I - M - P - O - R - T - A - N - T N - O - T - E *****
237 !
238 ! THE INPUT DATA FILE WILL ACCEPT BOTH POSITIVE AND NEGATIVE
239 ! SIGMAS, WHICH IT INTERPRETS DIFFERENTLY. FOR SIGMA > 0, IT
240 ! TAKES SIGMA TO BE THE STANDARD DEVIATION FOR THAT MEASUREMENT
241 ! AS DESCRIBED ABOVE. FOR SIGMA < 0, SIGMA IS INTERPRETED AS
242 ! THE FRACTION OF THE MEASURED DATA NEEDED TO COMPUTE THE ABSOLUTE
243 ! SIGMA, I.E., (-SIGMA * DATA) = ACTUAL SIGMA USED IN CODE.
244 !
245 
246  CALL second0(treadon)
247 
248  lwrite = (grank .EQ. 0)
249  ier_flag_init = ier_flag
250  ier_flag = norm_term_flag
251  IF (ier_flag_init .EQ. more_iter_flag) GOTO 1000
252 
253 !
254 ! READ IN DATA FROM INDATA FILE
255 !
256  CALL read_indata(input_file, iunit, ier_flag)
257  IF (ier_flag .NE. norm_term_flag) RETURN
258 
259  IF (tensi2 .EQ. zero ) tensi2 = tensi
260 
261 !
262 ! Open output files here, print out heading to threed1 file
263 !
264 ! PRINT *,'IN READIN, LWRITE: ', lwrite
265  IF (lwrite) THEN
266  CALL heading(input_extension, time_slice,
267  & iseq_count, lmac, lscreen, lwrite)
268  END IF
269 
270 !
271 ! READ IN COMMENTS DEMARKED BY "!"
272 !
273  rewind(iunit, iostat=iexit)
274  IF (lwrite) THEN
275  DO WHILE(iexit .EQ. 0)
276  READ (iunit, '(a)', iostat=iexit) line
277  IF (iexit .NE. 0) EXIT
278  iexit = index(line,'INDATA')
279  iexit = iexit + index(line,'indata')
280  ipoint = index(line,'!')
281  IF (ipoint .EQ. 1) WRITE (nthreed, *) trim(line)
282  ENDDO
283  END IF
284  CLOSE (iunit)
285 
286 !
287 ! READ IN AND STORE (FOR SEQUENTIAL RUNNING) MAGNETIC FIELD DATA
288 ! FROM MGRID_FILE
289 !
290  IF (lfreeb) THEN
291  CALL second0(trc)
292  CALL read_mgrid (mgrid_file, extcur, nzeta, nfp,
293  & lscreen, ier_flag, comm = runvmec_comm_world)
294  CALL second0(tzc)
295  mgrid_file_read_time = mgrid_file_read_time + (tzc - trc)
296 
297  IF (lfreeb .AND. lscreen .AND. lwrite) THEN
298  WRITE (6,'(2x,a,1p,e10.2,a)') 'Time to read MGRID file: ',
299  & tzc - trc, ' s'
300  IF (ier_flag .ne. norm_term_flag) RETURN
301  IF (lwrite) WRITE (nthreed,20) nr0b, nz0b, np0b, rminb,
302  & rmaxb, zminb, zmaxb, trim(mgrid_file)
303  20 FORMAT(//,' VACUUM FIELD PARAMETERS:',/,1x,24('-'),/,
304  & ' nr-grid nz-grid np-grid rmin rmax zmin',
305  & ' zmax input-file',/,3i9,4f10.3,5x,a)
306  END IF
307  END IF
308 
309 !
310 ! PARSE NS_ARRAY
311 !
312  nsin = max(3, nsin)
313  multi_ns_grid = 1
314  IF (ns_array(1) .eq. 0) THEN !Old input style
315  ns_array(1) = min(nsin,nsd)
316  multi_ns_grid = 2
317  ns_array(multi_ns_grid) = ns_default !Run on 31-point mesh
318  ELSE
319  nsmin = 1
320  DO WHILE (ns_array(multi_ns_grid) .ge. nsmin .and.
321  & multi_ns_grid .lt. 100) ! .ge. previously .gt.
322  nsmin = max(nsmin, ns_array(multi_ns_grid))
323  IF (nsmin .le. nsd) THEN
324  multi_ns_grid = multi_ns_grid + 1
325  ELSE !Optimizer, Boozer code overflows otherwise
326  ns_array(multi_ns_grid) = nsd
327  nsmin = nsd
328  IF (lwrite) THEN
329  print *,' NS_ARRAY ELEMENTS CANNOT EXCEED ',nsd
330  print *,' CHANGING NS_ARRAY(',multi_ns_grid,') to ',
331  & nsd
332  END IF
333  END IF
334  END DO
335  multi_ns_grid = multi_ns_grid - 1
336  ENDIF
337  IF (ftol_array(1) .eq. zero) THEN
338  ftol_array(1) = 1.e-8_dp
339  IF (multi_ns_grid .eq. 1) ftol_array(1) = ftol
340  DO igrid = 2, multi_ns_grid
341  ftol_array(igrid) = 1.e-8_dp * (1.e8_dp * ftol)**
342  & ( real(igrid-1,dp)/(multi_ns_grid-1) )
343  END DO
344  ENDIF
345 
346  ns_maxval = nsmin
347 !
348 ! WRITE OUT DATA TO THREED1 FILE
349 !
350 
351 !SPH121912 - SCALING TO RENDER LAMSCALE=1
352 ! delta = twopi/phiedge !phiedge=>twopi
353 ! phiedge = phiedge*delta
354 ! bcrit = bcrit*delta
355 ! curtor = curtor*delta
356 ! extcur = extcur*delta
357 ! am = am*delta**2
358 
359  IF (nvacskip .LE. 0) nvacskip = nfp
360 
361  proc0: IF (lwrite) THEN
362  WRITE (nthreed,100)
363  & ns_array(multi_ns_grid),ntheta1,nzeta,mpol,ntor,nfp,
364 #ifdef _ANIMEC
365  & gamma,spres_ped,phiedge,curtor,bcrit,lrfp
366 #else
367  & gamma,spres_ped,phiedge,curtor,lrfp
368 #endif
369  100 FORMAT(/,' COMPUTATION PARAMETERS: (u = theta, v = zeta)'/,
370  & 1x,45('-'),/,
371  & ' ns nu nv mu mv',/,
372  & 5i7,//,' CONFIGURATION PARAMETERS:',/,1x,39('-'),/,
373  & ' nfp gamma spres_ped phiedge(wb)'
374 #ifdef _ANIMEC
375  & ' curtor(A) BCrit(T) lRFP',
376  & /,i7,1p,e11.3,2e15.3,2e14.3,l12/)
377 #else
378  & ' curtor(A) lRFP',
379  & /,i7,1p,e11.3,2e15.3,e14.3,l12/)
380 #endif
381  WRITE (nthreed,110) ncurr,niter_array(multi_ns_grid),
382  & ns_array(1),nstep,nvacskip,
383  & ftol_array(multi_ns_grid),tcon0,lasym,lforbal,lmove_axis,
384  & lconm1,mfilter_fbdy,nfilter_fbdy,lfull3d1out,
385  & max_main_iterations,lgiveup,fgiveup ! M Drevlak 20130114
386  110 FORMAT(' RUN CONTROL PARAMETERS:',/,1x,23('-'),/,
387  & ' ncurr niter nsin nstep nvacskip ftol tcon0',
388  & ' lasym lforbal lmove_axis lconm1',/,
389  & 4i7,i10,1p,2e10.2,4l9,/,
390  & ' mfilter_fbdy nfilter_fbdy lfull3d1out max_main_iterations', ! J Geiger 20120203
391  & ' lgiveup fgiveup',/, ! M Drevlak 20130114
392  & 2(6x,i7),l12,10x,i10,l8,e9.1,/) ! M Drevlak 20130114
393 
394  WRITE (nthreed,120) precon_type, prec2d_threshold
395  120 FORMAT(' PRECONDITIONER CONTROL PARAMETERS:',/,1x,34('-'),/,
396  & ' precon_type prec2d_threshold',/,2x,a10,1p,e20.2,/)
397 
398  IF (nextcur .gt. 0) THEN
399  WRITE(nthreed, "(' EXTERNAL CURRENTS',/,1x,17('-'))")
400  ni = 0
401  IF (ALLOCATED(curlabel)) THEN
402  ni = maxval(len_trim(curlabel(1:nextcur)))
403  END IF
404  ni = max(ni+4, 14)
405  WRITE (line, '(a,i2.2,a)') "(5a",ni,")"
406  WRITE (line2, '(a,i2.2,a)') "(5(",ni-12,"x,1p,e12.4))"
407  DO i = 1,nextcur,5
408  ni = min(i+4, nextcur)
409  IF (ALLOCATED(curlabel)) THEN
410  WRITE (nthreed, line, iostat=mj)
411  & (trim(curlabel(n)),n=i,ni)
412  WRITE (nthreed, line2,iostat=mj) (extcur(n), n=i,ni)
413  END IF
414  END DO
415  WRITE (nthreed, *)
416  END IF
417 
418  IF (bloat .ne. one) THEN
419  WRITE (nthreed,'(" Profile Bloat Factor: ",1pe11.4)') bloat
420  phiedge = phiedge*bloat
421  END IF
422 
423  IF (pres_scale .ne. one) THEN
424  WRITE (nthreed,121) pres_scale
425  END IF
426  121 FORMAT(' Pressure profile factor: ',1pe11.4,
427  & ' (multiplier for pressure)')
428 ! Print out am array
429  WRITE(nthreed,130)
430  WRITE(nthreed,131) trim(pmass_type)
431  WRITE(nthreed,132)
432  130 FORMAT(' MASS PROFILE COEFFICIENTS - newton/m**2',
433  & ' (EXPANSION IN NORMALIZED RADIUS):')
434  131 FORMAT(' PMASS parameterization type is ''', a,'''')
435  132 FORMAT(1x,35('-'))
436 ! WRITE(nthreed,135)(am(i-1),i=1, SIZE(am))
437  135 FORMAT(1p,6e12.3)
438 
439  SELECT CASE(trim(pmass_type))
440  CASE ('Akima_spline','cubic_spline')
441  WRITE(nthreed,"(' am_aux_s is' )")
442  n = nonzerolen(am_aux_s,SIZE(am_aux_s))
443  WRITE(nthreed,135)(am_aux_s(i),i=1, n)
444  n = nonzerolen(am_aux_f,SIZE(am_aux_f))
445  WRITE(nthreed,"(' am_aux_f is' )")
446  WRITE(nthreed,135)(am_aux_f(i),i=1, n)
447  CASE DEFAULT
448  n = nonzerolen(am,SIZE(am))
449  WRITE(nthreed,135)(am(i-1),i=1,n)
450  END SELECT
451 
452  IF (ncurr.eq.0) THEN
453  IF (lrfp) THEN
454  WRITE (nthreed,142)
455  ELSE
456  WRITE (nthreed,140)
457  END IF
458 ! Print out ai array
459 ! WRITE(nthreed,135)(ai(i-1),i=1, SIZE(ai))
460  WRITE(nthreed,143) trim(piota_type)
461  SELECT CASE(trim(piota_type))
462  CASE ('Akima_spline','cubic_spline')
463  n = nonzerolen(ai_aux_s,SIZE(ai_aux_s))
464  WRITE(nthreed,"(' ai_aux_s is' )")
465  WRITE(nthreed,135)(ai_aux_s(i),i=1, n)
466  n = nonzerolen(ai_aux_f,SIZE(ai_aux_f))
467  WRITE(nthreed,"(' ai_aux_f is' )")
468  WRITE(nthreed,135)(ai_aux_f(i),i=1, n)
469  CASE DEFAULT
470  n = nonzerolen(ai,SIZE(ai))
471  WRITE(nthreed,135)(ai(i-1),i=1, n)
472  END SELECT
473  ELSE
474 ! Print out ac array
475  WRITE(nthreed,145)
476  WRITE(nthreed,146) trim(pcurr_type)
477  WRITE(nthreed,147)
478 ! WRITE(nthreed,135)(ac(i-1),i=1, SIZE(ac))
479  SELECT CASE(trim(pcurr_type))
480  CASE ('Akima_spline_Ip','Akima_spline_I', &
481  & 'cubic_spline_Ip','cubic_spline_I')
482  n = nonzerolen(ac_aux_s,SIZE(ac_aux_s))
483  WRITE(nthreed,"(' ac_aux_s is' )")
484  WRITE(nthreed,135)(ac_aux_s(i),i=1, n)
485  n = nonzerolen(ac_aux_f,SIZE(ac_aux_f))
486  WRITE(nthreed,"(' ac_aux_f is' )")
487  WRITE(nthreed,135)(ac_aux_f(i),i=1, n)
488  CASE DEFAULT
489  n = nonzerolen(ac,SIZE(ac))
490  WRITE(nthreed,135)(ac(i-1),i=1, n)
491  END SELECT
492  END IF
493 
494  140 FORMAT(/' IOTA PROFILE COEFFICIENTS',
495  & ' (EXPANSION IN NORMALIZED RADIUS):',/,1x,35('-'))
496  142 FORMAT(/' SAFETY-FACTOR (q) PROFILE COEFFICIENTS ai',
497  & ' (EXPANSION IN NORMALIZED RADIUS):',/,1x,35('-'))
498  143 FORMAT(' PIOTA parameterization type is ''', a,'''')
499  145 FORMAT(/' TOROIDAL CURRENT DENSITY (*V'') COEFFICIENTS',
500  & ' ac (EXPANSION IN NORMALIZED RADIUS):')
501  146 FORMAT(' PCURR parameterization type is ''', a,'''')
502  147 FORMAT(1x,38('-'))
503 
504  WRITE(nthreed,150)
505  n = nonzerolen(aphi,SIZE(aphi))
506  WRITE(nthreed,135)(aphi(i),i=1,n)
507  150 FORMAT(/' NORMALIZED TOROIDAL FLUX COEFFICIENTS aphi',
508  & ' (EXPANSION IN S):',/,1x,35('-'))
509 #ifdef _ANIMEC
510  IF (any(ah .ne. zero)) THEN
511  WRITE(nthreed,160)
512  n = nonzerolen(ah,SIZE(ah))
513  WRITE(nthreed,135)(ah(i-1),i=1, n)
514  WRITE(nthreed,165)
515  n = nonzerolen(at,SIZE(at))
516  WRITE(nthreed,135)(at(i-1),i=1, n)
517  END IF
518 
519  160 FORMAT(' HOT PARTICLE PRESSURE COEFFICIENTS ah',
520  & ' (EXPANSION IN TOROIDAL FLUX):',/,1x,35('-'))
521  165 FORMAT(' HOT PARTICLE TPERP/T|| COEFFICIENTS at',
522  & ' (EXPANSION IN TOROIDAL FLUX):',/,1x,35('-'))
523 #endif
524 
525 ! Fourier Boundary Coefficients
526  WRITE(nthreed,180)
527  180 FORMAT(/,' R-Z FOURIER BOUNDARY COEFFICIENTS AND',
528  & ' MAGNETIC AXIS INITIAL GUESS',/,
529  & ' R = RBC*COS(m*u - n*v) + RBS*SIN(m*u - n*v),',
530  & ' Z = ZBC*COS(m*u - n*v) + ZBS*SIN(m*u-n*v)'/1x,86('-'),
531  & /,' nb mb rbc rbs zbc ',
532  & 'zbs ',
533  & ' raxis(c) raxis(s) zaxis(c) zaxis(s)')
534 
535  END IF proc0
536 
537 1000 CONTINUE
538 
539 !
540 ! CONVERT TO REPRESENTATION WITH RBS(m=1) = ZBC(m=1)
541 !
542  IF (lasym) THEN
543  delta = atan((rbs(0,1) - zbc(0,1))/
544  & (abs(rbc(0,1)) + abs(zbs(0,1))))
545  IF (delta .ne. zero) THEN
546  DO m = 0,mpol1
547  DO n = -ntor,ntor
548  trc = rbc(n,m)*cos(m*delta) + rbs(n,m)*sin(m*delta)
549  rbs(n,m) = rbs(n,m)*cos(m*delta)
550  & - rbc(n,m)*sin(m*delta)
551  rbc(n,m) = trc
552  tzc = zbc(n,m)*cos(m*delta) + zbs(n,m)*sin(m*delta)
553  zbs(n,m) = zbs(n,m)*cos(m*delta)
554  & - zbc(n,m)*sin(m*delta)
555  zbc(n,m) = tzc
556  END DO
557  END DO
558  END IF
559  END IF
560 
561 !
562 ! ALLOCATE MEMORY FOR NU, NV, MPOL, NTOR SIZED ARRAYS
563 !
564  CALL allocate_nunv
565 
566 !
567 ! CONVERT TO INTERNAL REPRESENTATION OF MODES
568 !
569 ! R = RBCC*COS(M*U)*COS(N*V) + RBSS*SIN(M*U)*SIN(N*V)
570 ! + RBCS*COS(M*U)*SIN(N*V) + RBSC*SIN(M*U)*COS(N*V)
571 ! Z = ZBCS*COS(M*U)*SIN(N*V) + ZBSC*SIN(M*U)*COS(N*V)
572 ! + ZBCC*COS(M*U)*COS(N*V) + ZBSS*SIN(M*U)*SIN(N*V)
573 !
574 !
575 ! POINTER ASSIGNMENTS (NOTE: INDICES START AT 1, NOT 0, FOR POINTERS, EVEN THOUGH
576 ! THEY START AT ZERO FOR RMN_BDY)
577 ! ARRAY STACKING ORDER DETERMINED HERE
578 !
579  rbcc => rmn_bdy(:,:,rcc)
580  zbsc => zmn_bdy(:,:,zsc)
581  IF (lthreed) THEN
582  rbss => rmn_bdy(:,:,rss)
583  zbcs => zmn_bdy(:,:,zcs)
584  END IF
585 
586  IF (lasym) THEN
587  rbsc => rmn_bdy(:,:,rsc)
588  zbcc => zmn_bdy(:,:,zcc)
589  IF (lthreed) THEN
590  rbcs => rmn_bdy(:,:,rcs)
591  zbss => zmn_bdy(:,:,zss)
592  END IF
593  ENDIF
594 
595  rmn_bdy = 0; zmn_bdy = 0
596 
597  ioff = lbound(rbcc,1)
598  joff = lbound(rbcc,2)
599 
600  DO m=0, mpol1
601  mj = m + joff
602  IF (lfreeb .and.
603  & (mfilter_fbdy.gt.1 .and. m.gt.mfilter_fbdy)) cycle
604  DO n = -ntor, ntor
605  IF (lfreeb .and.
606  & (nfilter_fbdy.gt.0 .and. abs(n).gt.nfilter_fbdy)) cycle
607  ni = abs(n) + ioff
608  IF (n .eq. 0) THEN
609  isgn = 0
610  ELSE IF (n .gt. 0) THEN
611  isgn = 1
612  ELSE
613  isgn = -1
614  END IF
615  rbcc(ni,mj) = rbcc(ni,mj) + rbc(n,m)
616  IF (m .gt. 0) zbsc(ni,mj) = zbsc(ni,mj) + zbs(n,m)
617 
618  IF (lthreed) THEN
619  IF (m .gt. 0) rbss(ni,mj) = rbss(ni,mj) + isgn*rbc(n,m)
620  zbcs(ni,mj) = zbcs(ni,mj) - isgn*zbs(n,m)
621  END IF
622 
623  IF (lasym) THEN
624  IF (m .gt. 0) rbsc(ni,mj) = rbsc(ni,mj) + rbs(n,m)
625  zbcc(ni,mj) = zbcc(ni,mj) + zbc(n,m)
626  IF (lthreed) THEN
627  rbcs(ni,mj) = rbcs(ni,mj) - isgn*rbs(n,m)
628  IF (m .gt. 0) zbss(ni,mj) = zbss(ni,mj) + isgn*zbc(n,m)
629  END IF
630  END IF
631 
632  IF (ier_flag_init .ne. norm_term_flag) cycle
633  trc = abs(rbc(n,m)) + abs(rbs(n,m))
634  & + abs(zbc(n,m)) + abs(zbs(n,m))
635  IF (m .eq. 0) THEN
636  IF (n .lt. 0) cycle
637  IF (trc.eq.zero .and. abs(raxis_cc(n)).eq.zero .and.
638  & abs(zaxis_cs(n)).eq.zero) cycle
639  IF (lwrite) WRITE (nthreed,195) n, m, rbc(n,m), rbs(n,m),
640  & zbc(n,m), zbs(n,m), raxis_cc(n), raxis_cs(n),
641  & zaxis_cc(n), zaxis_cs(n)
642  ELSE
643  IF (trc .eq. zero) cycle
644  IF (lwrite) WRITE (nthreed,195) n, m, rbc(n,m), rbs(n,m),
645  & zbc(n,m), zbs(n,m)
646  END IF
647  END DO
648  END DO
649  195 FORMAT(i5,i4,1p,8e12.4)
650 
651 !
652 ! CHECK SIGN OF JACOBIAN (SHOULD BE SAME AS SIGNGS)
653 !
654  m = 1
655  mj = m+joff
656  rtest = sum(rbcc(1:ntor1,mj))
657  ztest = sum(zbsc(1:ntor1,mj))
658  lflip=(rtest*ztest .lt. zero)
659  signgs = -1
660  IF (lflip) CALL flip_theta(rmn_bdy, zmn_bdy)
661 
662 !
663 ! CONVERT TO INTERNAL FORM FOR (CONSTRAINED) m=1 MODES
664 ! INTERNALLY, FOR m=1: XC(rss) = .5(RSS+ZCS), XC(zcs) = .5(RSS-ZCS)
665 ! WITH XC(zcs) -> 0 FOR POLAR CONSTRAINT
666 ! (see convert_sym, convert_asym in totzsp_mod file)
667 !
668 
669  IF (lconm1 .and. (lthreed .or. lasym)) THEN
670  ALLOCATE (temp(SIZE(rbcc,1)))
671  IF (lthreed) THEN
672  mj = 1+joff
673  temp = rbss(:,mj)
674  rbss(:,mj) = p5*(temp(:) + zbcs(:,mj))
675  zbcs(:,mj) = p5*(temp(:) - zbcs(:,mj))
676  END IF
677  IF (lasym) THEN
678  mj = 1+joff
679  temp = rbsc(:,mj)
680  rbsc(:,mj) = p5*(temp(:) + zbcc(:,mj))
681  zbcc(:,mj) = p5*(temp(:) - zbcc(:,mj))
682  END IF
683  IF (ALLOCATED(temp)) DEALLOCATE (temp)
684  END IF
685 
686 !
687 ! PARSE TYPE OF PRECONDITIONER
688 !
689  precon_type = trim(adjustl(precon_type))
690  itype_precon = 0 !default scalar tri-di preconditioner
691  lprecond = .false.
692  ch1 = precon_type(1:1); ch2 = precon_type(2:2)
693 
694 ! ALL THE FOLLOWING USE THE FULL 2D BLOCK-TRI PRECONDITIONER
695 ! BUT DIFFER IN THE WAY TIME-EVOLUTION IS HANDLED
696  SELECT CASE (ch1)
697  CASE ('c', 'C')
698 !conjugate gradient
699  IF (ch2 == 'g' .or. ch2 == 'G') itype_precon = 1
700  lprecond = .true.
701  CASE ('g', 'G')
702 !gmres or gmresr
703  IF (ch2 == 'm' .or. ch2 == 'M') itype_precon = 2
704  IF (len_trim(precon_type) == 6) itype_precon = 3
705  lprecond = .true.
706  CASE ('t', 'T')
707 !transpose free qmr
708  IF (ch2 == 'f' .or. ch2 == 'F') itype_precon = 4
709  lprecond = .true.
710  END SELECT
711 
712 
713  iresidue = -1
714 
715  currv = mu0*curtor !Convert to Internal units
716 
717  CALL second0(treadoff)
718  timer(tread) = timer(tread) + (treadoff-treadon)
719  CALL mpi_bcast(lprecond,1,mpi_logical,0,runvmec_comm_world, &
720  & mpi_err)
721  readin_time = timer(tread)
722 
723  END SUBROUTINE readin
724 
725  INTEGER FUNCTION nonzerolen(array, n)
726  USE vmec_main, ONLY: dp, zero
727  IMPLICIT NONE
728  INTEGER, INTENT(IN) :: n
729  REAL(dp), INTENT(IN) :: array(n)
730  INTEGER :: k
731 
732  DO k = n, 1, -1
733  IF (array(k) .NE. zero) EXIT
734  END DO
735 
736  nonzerolen = k
737 
738  END FUNCTION nonzerolen