PB3D [2.47]
Ideal linear high-n MHD stability in 3-D
Loading...
Searching...
No Matches
input_ops.f90
Go to the documentation of this file.
1!------------------------------------------------------------------------------!
2!> Operations concerning giving input.
3!------------------------------------------------------------------------------!
5#include <PB3D_macros.h>
7 use output_ops
8 use messages
9 use num_vars, only: dp, max_str_ln
10
11 implicit none
12 private
14
15contains
16 !> Reads input options from user-provided input file.
17 !!
18 !! This is done by the master process and these variables will then be
19 !! passed on through MPI to other processes in broadcast_input_opts().
20 !!
21 !! \return ierr
22 integer function read_input_opts() result(ierr)
23 use num_vars, only: &
31 &bc_style, tol_norm, tol_slepc_loc => tol_slepc, max_it_slepc, &
41 use eq_vars, only: rho_0, r_0, pres_0, b_0, psi_0, t_0
48
49 character(*), parameter :: rout_name = 'read_input_opts'
50
51 ! local variables
52 character(len=max_str_ln) :: err_msg !< error message
53 integer :: pb3d_rich_lvl !< Richardson level to post-process (for POST)
54 integer, parameter :: max_size_tol_slepc = 100 !< maximum size of tol_SLEPC
55 real(dp) :: alpha !< field line label (for alpha_style 1)
56 real(dp) :: tol_slepc(max_size_tol_slepc) !< tol_SLEPC
57 real(dp), parameter :: min_tol = 1.e-12_dp !< minimum general tolerance
58 real(dp), parameter :: max_tol = 1.e-3_dp !< maximum general tolerance
59
60 ! input options
61 namelist /inputdata_pb3d/ min_par_x, max_par_x, alpha, min_r_sol, &
75 namelist /inputdata_post/ n_sol_plotted, n_theta_plot, n_zeta_plot, &
78 &pb3d_rich_lvl, max_it_zero, tol_zero, min_theta_plot, &
84
85 ! initialize ierr
86 ierr = 0
87
88 if (rank.eq.0) then ! only master
89 call writo('Setting up user-provided input "'//trim(input_name)//&
90 &'"')
91 call lvl_ud(1)
92
93 ! initialize input variables (optionally overwritten by user later)
94 call writo('Initialize all the inputs with default values')
95 call lvl_ud(1)
96
97 ! common variables for all program styles
98 sol_n_procs = n_procs ! use equal amount of processes for solution as for rest of program
99 max_tot_mem = 6000_dp ! count with 6GB total
100 plot_size = [10,5] ! size of plot in inch
101 min_r_plot = 0._dp ! minimal plot normal range
102 max_r_plot = 1._dp ! maximal plot normal range
103 min_theta_plot = 1 ! starting from pi gives nicer plots
105 select case(eq_style)
106 case (1) ! VMEC
107 n_theta_plot = 201 ! nr. poloidal points in plot
108 n_zeta_plot = 51 ! nr. toroidal points in plot
109 min_zeta_plot = 0 ! min. toroidal plot angle [pi]
110 max_zeta_plot = 2 ! max. toroidal plot angle [pi]
111 case (2) ! HELENA
112 n_theta_plot = 501 ! nr. poloidal points in plot
113 n_zeta_plot = 1 ! nr. toroidal points in plot
114 min_zeta_plot = 0 ! min. toroidal plot angle [pi]
115 max_zeta_plot = min_zeta_plot ! max. toroidal plot angle [pi]
116 end select
117 plot_grid_style = 0 ! normal plots on 3D geometry
118 max_nr_backtracks_hh = 20 ! standard nr. of backtracks
119 ex_plot_style = 1 ! GNUPlot
120 norm_disc_prec_eq = 3 ! cubic precision for normal discretization of equilibrium
121 norm_disc_prec_x = 3 ! cubic precision for normal discretization of perturbation
122 norm_disc_prec_sol = 3 ! cubic precision for normal discretization of solution
123
124 ! select depending on program style
125 select case (prog_style)
126 case(1) ! PB3D
127 call default_input_pb3d()
128 case(2) ! POST
129 call default_input_post()
130 end select
131
132 call lvl_ud(-1)
133
134 ! read user input
135 ! select depending on program style
136 select case (prog_style)
137 case(1) ! PB3D
138 read(unit=input_i,nml=inputdata_pb3d,iostat=ierr,&
139 &iomsg=err_msg) ! read input data
140
141 ! check input if successful read
142 if (ierr.eq.0) then ! input file succesfully read
143 call writo('Overwriting with user-provided file "'&
144 &//trim(input_name) // '"')
145
146 call writo('Checking user-provided file')
147
148 call lvl_ud(1)
149
150 ! adapt MPI variables if needed
151 call adapt_mpi()
152
153 ! adapt run-time variables if needed
154 ierr = adapt_run_pb3d()
155 chckerr('')
156
157 ! adapt plotting variables if needed
158 ierr = adapt_plot()
159 chckerr('')
160
161 ! adapt perturbation modes
162 ierr = adapt_x_modes()
163 chckerr('')
164
165 ! adapt min_n_par_X if needed
166 call adapt_min_n_par_x
167
168 ! adapt alpha if needed
169 ierr = adapt_alpha()
170 chckerr('')
171
172 ! adapt solution grid
173 ierr = adapt_sol_grid(min_r_sol,max_r_sol,'sol')
174 chckerr('')
175
176 ! adapt normalization variables if needed
177 ierr = adapt_normalization()
178 chckerr('')
179
180 ! adapt Richardson variables if needed
181 call adapt_rich
182
183 ! adapt input / output variables if needed
184 ierr = adapt_inoutput_pb3d()
185 chckerr('')
186
187 ! adapt Householder variables if needed
188 call adapt_zero
189
190 ! adapt tolerances if needed
191 call adapt_tols
192
193 call lvl_ud(-1)
194 else ! cannot read input data
195 ierr = 1
196 call writo('Cannot read user-provided file "'&
197 &//trim(input_name)//'", error message:',&
198 &error=.true.)
199 call lvl_ud(1)
200 call writo('"'//trim(err_msg)//'"')
201 call lvl_ud(-1)
202 chckerr("")
203 end if
204 case(2) ! POST
205 read(unit=input_i,nml=inputdata_post,iostat=ierr,&
206 &iomsg=err_msg) ! read input data
207
208 ! check input if successful read
209 if (ierr.eq.0) then ! input file succesfully read
210 call writo('Overwriting with user-provided file "'&
211 &//trim(input_name) // '"')
212
213 call writo('Checking user-provided file')
214
215 call lvl_ud(1)
216
217 ! adapt run-time variables if needed
218 ierr = adapt_run_post()
219 chckerr('')
220
221 ! adapt plotting variables if needed
222 ierr = adapt_plot()
223 chckerr('')
224
225 ! adapt plot grid and save to solution variable
226 ierr = adapt_sol_grid(min_r_plot,max_r_plot,'plot')
227 chckerr('')
228
229 ! adapt input / output variables if needed
230 ierr = adapt_inoutput_post()
231 chckerr('')
232
233 call lvl_ud(-1)
234 else ! cannot read input data
235 ierr = 1
236 call writo('Cannot read user-provided file "'&
237 &//trim(input_name)//'", error message:',&
238 &error=.true.)
239 call lvl_ud(1)
240 call writo('"'//trim(err_msg)//'"')
241 call lvl_ud(-1)
242 chckerr("")
243 end if
244 end select
245
246 ! user output
247 call writo('Close input file')
248 close(input_i)
249
250 call lvl_ud(-1)
251 call writo('Input values set')
252 end if
253 contains
254 ! PB3D version
255 !> \private
256 subroutine default_input_pb3d
257 !use num_vars, only: use_pol_flux_E
258
259 ! concerning solution
260 n_r_sol = 100 ! 100 points in solution grid
261 min_r_sol = 0.1_dp ! minimum normal range
262 max_r_sol = 1.0_dp ! maximum normal range
263 tol_norm = 0.05 ! tolerance for normal range
264 ev_style = 1 ! slepc solver for EV problem
265
266 ! concerning field line
267 alpha = 0._dp ! field line based in outboard
268 min_alpha = 0.0_dp ! minimum field-line label [pi]
269 max_alpha = 2.0_dp ! maximum field-line label [pi]
270 select case (eq_style)
271 case (1) ! VMEC
272 alpha_style = 2 ! multiple field-lines, single turns
273 n_alpha = 10
274 case (2) ! HELENA
275 alpha_style = 1 ! single field-line, multiple turns
276 end select
277
278 ! concerning perturbation
279 min_n_par_x = req_min_n_par_x ! nonsensible value to check for user overwriting
280 min_par_x = 0.0_dp ! minimum parallel angle [pi]
281 max_par_x = 2.0_dp ! maximum parallel angle [pi]
282 prim_x = 20 ! main mode number of perturbation
283 min_sec_x = huge(1) ! nonsensible value to check for user overwriting
284 max_sec_x = huge(1) ! nonsensible value to check for user overwriting
285 n_mod_x = huge(1) ! nonsensible value to check for user overwriting
286 !use_pol_flux_F = use_pol_flux_E ! use same normal flux coordinate as the equilibrium
287 use_pol_flux_f = .true. ! use poloidail flux as toroidal flux is untested
288
289 ! concerning normalization
290 rho_0 = 1e-7_dp ! for fusion, particle density of around 5E19 for detuerium
291 r_0 = huge(1._dp) ! nonsensible value to check for user overwriting
292 pres_0 = huge(1._dp) ! nonsensible value to check for user overwriting
293 psi_0 = huge(1._dp) ! nonsensible value to check for user overwriting
294 b_0 = huge(1._dp) ! nonsensible value to check for user overwriting
295 t_0 = huge(1._dp) ! nonsensible value to check for user overwriting
296
297 ! concerning input / output
298 n_sol_requested = 1 ! request solution closest to guess
299 retain_all_sol = .false. ! don't retain faulty ones
300 plot_resonance = .false. ! do not plot the q-profile with nq-m = 0
301 plot_magn_grid = .false. ! do not plot the magnetic grid
302 plot_kappa = .false. ! do not plot the curvature
303 plot_b = .false. ! do not plot the magnetic field
304 plot_j = .false. ! do not plot the current
305 plot_flux_q = .false. ! do not plot the flux quantities
306
307 ! concerning Richardson extrapolation
308 max_it_rich = 1 ! by default no Richardson extrapolation
309 tol_rich = 1.e-4_dp ! tolerance of 1E-4
310 rich_restart_lvl = 1 ! don't restart
311
312 ! concerning finding zeros
313 max_it_zero = 100 ! maximum 100 iterations
314 tol_zero = 1.0e-10_dp ! very low tolerance for calculation of field lines
315
316 ! runtime variables
317 use_normalization = .true. ! use normalization for the variables
318 norm_disc_style_sol = 2 ! left finite differences
319 magn_int_style = 1 ! trapezoidal rule
320 max_it_slepc = 1000 ! max. nr. of iterations for SLEPC
321 ev_bc = 1._dp ! use 1 as artificial EV for the Boundary Conditions
322 ev_guess = -0.3_dp ! sensible guess
323 tol_slepc = huge(1._dp) ! nonsensible value to check for user overwriting
324 rho_style = 1 ! constant pressure profile, equal to rho_0
325 u_style = 3 ! full expression for U, up to order 3
326 k_style = 1 ! perpendicular kinetic energy normalized
327 norm_style = 1 ! MISHKA normalization
328 bc_style = [1,2] ! left BC zeroed and right BC through minimization of energy
329 x_style = 2 ! fast style: mode numbers optimized in normal coordinate
330 solver_slepc_style = 1 ! Krylov-Schur
331 matrix_slepc_style = 1 ! sparse matrix storage
332 x_grid_style = huge(1) ! nonsensible value to check for user overwriting
333 max_njq_change = 0.49_dp ! maximum change of just under 0.5 for n q (pol. flux) or n iota (tor. flux)
334 end subroutine default_input_pb3d
335
336 ! POST version
337 !> \private
338 subroutine default_input_post()
339 ! concerning finding zeros
340 max_it_zero = 200 ! more iterations than PB3D
341 tol_zero = 1.0e-8_dp ! less relative error than PB3D
342
343 ! runtime variables
344 plot_resonance = .true. ! plot the q-profile with nq-m = 0
345 plot_flux_q = .true. ! plot the flux quantities
346 plot_kappa = .true. ! plot the curvature
347 plot_magn_grid = .true. ! plot the magnetic grid
348 plot_b = .true. ! plot the magnetic field
349 plot_j = .true. ! plot the current
350 plot_sol_xi = .true. ! plot plasma perturbation of solution
351 plot_sol_q = .true. ! plot magnetic perturbation of solution
352 plot_e_rec = .true. ! plot energy reconstruction
353 plot_vac_pot = .true. ! plot vacuum potential
354 post_style = 1 ! process on extended plot grid
355 write(*,*) '!!!!!! min and max need to take into account normalization'
356 min_rvac_plot = 0.1*r_0 ! minimum R of vacuum plot
357 max_rvac_plot = 2*r_0 ! maximum R of vacuum plot
358 min_zvac_plot = -r_0 ! minimum R of vacuum plot
359 max_zvac_plot = r_0 ! maximum R of vacuum plot
360 n_vac_plot = [100,100] ! number of points in R and Z of vacuum plot
361
362 ! variables concerning input / output
363 pert_mult_factor_post = 0._dp ! factor by which to XYZ is perturbed in POST
364 pb3d_rich_lvl = huge(1) ! don't restart
365
366 ! variables concerning output
367 n_sol_plotted = -1 ! plot all solutions
368 end subroutine default_input_post
369
370 ! Checks whether the variables concerning MPI are chosen correctly:
371 ! sol_n_procs cannot be greater than n_procs. If it is lower than 1,
372 ! the maximum amount of n_procs is used.
373 !> \private
374 subroutine adapt_mpi()
375 use num_vars, only: n_procs
376
377 if (sol_n_procs.lt.1) then
378 call writo('Using the maximum number of MPI processes '//&
379 &trim(i2str(n_procs))//' for SLEPC')
381 end if
382 if (sol_n_procs.gt.n_procs) then
383 call writo('Cannot use more than '//trim(i2str(n_procs))//&
384 &' MPI processes for SLEPC',warning=.true.)
386 end if
387 end subroutine adapt_mpi
388
389 ! Checks whether the variables concerning run-time are chosen correctly:
390 ! rho_style has to be 1 (constant rho = rho_0),
391 ! matrix_SLEPC_style has to be 0 or 1,
392 ! max_it_SLEPC has to be at least 1,
393 ! magnetic integral style has to be 1..2,
394 ! perturbation grid style has to be 1..2,
395 ! tensorial perturbation variables interpolation style has to be 1..2
396 ! (for style 2, a warning should be displayed)
397 ! for HELENA (eq_style 1), only poloidal flux can be used.
398 ! norm_disc_prec_eq has to be 3 because cubics are the maximum order
399 ! for splines and derivatives of order 3 are needed for equilibrium
400 ! quantities.
401 ! The other normal discretization precisions can be lower, but not
402 ! higher.
403 !> \private
404 integer function adapt_run_pb3d() result(ierr)
405 use num_vars, only: eq_style
406
407 character(*), parameter :: rout_name = 'adapt_run_PB3D'
408
409 ! initialize ierr
410 ierr = 0
411
412 if (rho_style.ne.1) then
413 ierr = 1
414 err_msg = 'rho_style has to be 1 (constant)'
415 chckerr(err_msg)
416 end if
417 if (eq_style.eq.2 .and. .not.use_normalization) then
418 use_normalization = .true.
419 call writo('normalization is always used for HELENA',&
420 &warning=.true.)
421 end if
422 if (matrix_slepc_style.lt.1 .or. matrix_slepc_style.gt.2) then
423 ierr = 1
424 err_msg = 'matrix_SLEPC_style has to be 1 (sparse) or 2 (shell)'
425 chckerr(err_msg)
426 end if
427 if (max_it_slepc.le.0) then
428 max_it_slepc = 1
429 call writo('max_it_SLEPC has to be positive and is set to 1',&
430 &warning=.true.)
431 end if
432 if (magn_int_style.lt.1 .or. magn_int_style.gt.2) then
433 ierr = 1
434 err_msg = 'magn_int_style has to be 1 (trapezoidal) or 2 &
435 &(Simpson 3/8 rule)'
436 chckerr(err_msg)
437 end if
438 if (eq_style.eq.2 .and. .not.use_pol_flux_f) then
439 use_pol_flux_f = .true.
440 call writo('can only use poloidal flux for HELENA',&
441 &warning=.true.)
442 end if
443 if (norm_disc_prec_eq.ne.3) then
444 ierr =1
445 err_msg = 'can only use norm_disc_prec_eq = 3 currently'
446 chckerr(err_msg)
447 end if
448 if (norm_disc_prec_x.lt.1 .or. norm_disc_prec_x.gt.3) then
449 ierr =1
450 err_msg = 'norm_disc_prec_X has to be between 1 and 3'
451 chckerr(err_msg)
452 end if
453 if (norm_disc_prec_sol.lt.1 .or. norm_disc_prec_sol.gt.3) then
454 ierr = 1
455 err_msg = 'norm_disc_prec_sol has to be between 1 and 3'
456 chckerr(err_msg)
457 end if
458 end function adapt_run_pb3d
459
460 ! Checks whether the variables concerning run-time are chosen correctly:
461 ! POST_style should be 1 (plotting on extended plot grid) or 2
462 ! (plotting on field-aligned grid also used in PB3D).
463 !> \private
464 integer function adapt_run_post() result(ierr)
465 character(*), parameter :: rout_name = 'adapt_run_POST'
466
467 ! initialize ierr
468 ierr = 0
469
470 if (post_style.lt.1 .or. post_style.gt.2) then
471 ierr = 1
472 err_msg = 'POST_style has to be 1 (extended grid) or 2 &
473 &(field-aligned grid)'
474 chckerr(err_msg)
475 end if
476 end function adapt_run_post
477
478 ! Checks whether the variables concerning output are chosen correctly:
479 ! PB3D_rich_lvl can be at most the maximally found level
480 ! POST_output_full is true if non-flux quantities are plot
481 ! POST_output_sol is true if solution quantities are plot
482 !> \private
483 integer function adapt_inoutput_post() result(ierr)
484 use num_vars, only: compare_tor_pos
485
486 character(*), parameter :: rout_name = 'adapt_inoutput_POST'
487
488 ! local variables
489 integer :: max_pb3d_rich_lvl ! maximum Richardson level
490
491 ! initialize ierr
492 ierr = 0
493
494 ! find maximum levl
495 ierr = find_max_rich_lvl(pb3d_rich_lvl,max_pb3d_rich_lvl)
496 chckerr('')
497 if (max_pb3d_rich_lvl.le.0) then
498 call writo('No suitable Richardson level found, so only &
499 &equilibrium output will be done.',alert=.true.)
500 plot_sol_xi = .false.
501 plot_sol_q = .false.
502 plot_e_rec = .false.
503 pb3d_rich_lvl = 1
504 else
505 if (pb3d_rich_lvl.eq.huge(1)) then
506 pb3d_rich_lvl = max_pb3d_rich_lvl
507 call writo('Treating maximum Richardson level found: '//&
508 &trim(i2str(max_pb3d_rich_lvl)))
509 else if (pb3d_rich_lvl.le.max_pb3d_rich_lvl) then
510 call writo('Treating Richardson level requested: '//&
511 &trim(i2str(pb3d_rich_lvl)))
512 if (pb3d_rich_lvl.lt.max_pb3d_rich_lvl) then
513 call lvl_ud(1)
514 call writo('Less than maximum Richardson level: '//&
515 &trim(i2str(max_pb3d_rich_lvl)),warning=.true.)
516 call lvl_ud(-1)
517 end if
518 else if (pb3d_rich_lvl.gt.max_pb3d_rich_lvl) then
519 ierr = 1
520 err_msg = 'PB3D_rich_lvl cannot be higher than '//&
521 &trim(i2str(max_pb3d_rich_lvl))
522 chckerr(err_msg)
523 else if (pb3d_rich_lvl.lt.1) then
524 call writo('Richardson level requested '//&
525 &trim(i2str(pb3d_rich_lvl))//', so only equilibrium &
526 &output will be done.',alert=.true.)
527 plot_sol_xi = .false.
528 plot_sol_q = .false.
529 plot_e_rec = .false.
530 pb3d_rich_lvl = 1
531 end if
532 end if
533 max_it_rich = pb3d_rich_lvl
534 rich_lvl = pb3d_rich_lvl
535
536 ! set POST_output_full and POST_output_sol
540 end function adapt_inoutput_post
541
542 ! Checks whether the variables concerning output are chosen correctly:
543 ! n_sol_requested has to be at least one,
544 ! rich_restart_lvl can be at most one more than the maximally found
545 ! level, nor can it be higher than max_it_rich.
546 !> \private
547 integer function adapt_inoutput_pb3d() result(ierr)
548 character(*), parameter :: rout_name = 'adapt_inoutput_PB3D'
549
550 ! initialize ierr
551 ierr = 0
552
553 if (n_sol_requested.lt.1) then
554 n_sol_requested = 1
555 call writo('n_sol_requested has been increased to '&
556 &//trim(i2str(n_sol_requested)),warning=.true.)
557 end if
558 if (rich_restart_lvl.lt.1) then
559 call writo('rich_restart_lvl was '//&
560 &trim(i2str(rich_restart_lvl))//&
561 &' but cannot be lower than 1, so it was reset to 1',&
562 &warning=.true.)
563 rich_restart_lvl = 1
564 end if
565 if (rich_restart_lvl.gt.max_it_rich) then
566 ierr = 1
567 err_msg = 'rich_restart_lvl not within 1..max_it_rich = 1..'//&
568 &trim(i2str(max_it_rich))
569 chckerr(err_msg)
570 end if
571 if (rich_restart_lvl.gt.1) then
572 ierr = find_max_rich_lvl(pb3d_rich_lvl,pb3d_rich_lvl)
573 chckerr('')
574 if (rich_restart_lvl.gt.pb3d_rich_lvl+1) then
575 ierr = 1
576 if (pb3d_rich_lvl.gt.0) then
577 err_msg = 'The highest Richardson level found was '//&
578 &trim(i2str(pb3d_rich_lvl))
579 else
580 err_msg = 'No Richardson level found'
581 end if
582 chckerr(err_msg)
583 end if
584 end if
585 end function adapt_inoutput_pb3d
586
587 ! Checks whether the variables concerning plotting are chosen correctly:
588 ! n_theta_plot and n_zeta_plot have to be positive,
589 ! if n_theta_plot or n_zeta_plot is equal to 1, the minimum value is
590 ! chosen.
591 ! if comparing different toroidal positions, only 2 values can be
592 ! used.
593 ! ex_plot_style should be 1 (GNUPlot) or 2 (Bokeh / Mayavi)
594 ! if toroidal positions are compared, theta_plot limits should contain
595 ! 0..2pi and max_r_plot should be 1.
596 !> \private
597 integer function adapt_plot() result(ierr)
598 use num_vars, only: compare_tor_pos
599
600 character(*), parameter :: rout_name = 'adapt_plot'
601
602 ! local variables
603 real(dp), parameter :: tol = 1.e-6_dp ! tolerance for checks
604 character(len=max_str_ln) :: err_msg ! error message
605
606 ! initialize ierr
607 ierr = 0
608
609 if (n_theta_plot.lt.1) then
610 n_theta_plot = 1
611 call writo('n_theta_plot has to be positive and is set to '//&
612 &trim(i2str(n_theta_plot)),warning=.true.)
613 end if
614 if (n_zeta_plot.lt.1) then
615 n_zeta_plot = 1
616 call writo('n_zeta_plot has to be positive and is set to '//&
617 &trim(i2str(n_zeta_plot)),warning=.true.)
618 end if
619 if (n_theta_plot.eq.1 .and. &
620 &abs(max_theta_plot-min_theta_plot).gt.tol) then
621 call writo('For n_theta_plot = 1, max_theta_plot is changed to &
622 &min_theta_plot = '//trim(r2strt(min_theta_plot)),&
623 &warning=.true.)
625 end if
626 if (n_zeta_plot.eq.1 .and. &
627 &abs(max_zeta_plot-min_zeta_plot).gt.tol) then
628 call writo('For n_zeta_plot = 1, min_zeta_plot and &
629 &max_zeta_plot are changed to their average = '//&
631 &warning=.true.)
634 end if
635 if (ex_plot_style.lt.1 .or. ex_plot_style.gt.2) then
636 call writo('external plot style should be')
637 call lvl_ud(1)
638 call writo('1: GNUPlot')
639 call writo('2: Bokeh (2-D) / Mayavi (3-D)')
640 call lvl_ud(-1)
641 call writo('Default (1: GNUPlot) chosen',warning=.true.)
642 ex_plot_style = 1
643 end if
644 if (prog_style.eq.2) then ! only for POST
645 if (compare_tor_pos) then
646 if (n_zeta_plot.ne.3) then
647 ierr = 1
648 err_msg = 'For compare_tor_pos, you can only use &
649 &n_zeta = 3 (one point in the middle)'
650 chckerr(err_msg)
651 end if
652 if (post_style.ne.1) then
653 ierr = 1
654 err_msg = 'For compare_tor_pos, you can only use &
655 &POST_style = 1 (extended grid)'
656 chckerr(err_msg)
657 end if
658 if (min_r_plot.lt.tol_zero) then
659 ierr = 1
660 err_msg = 'min_r_plot has to be greater than 0.'
661 chckerr(err_msg)
662 end if
663 if (max_r_plot.lt.1._dp) then
664 ierr = 1
665 err_msg = 'max_r_plot has to be 1.'
666 chckerr(err_msg)
667 end if
668 if (abs(min_theta_plot-0).gt.tol_zero .or. &
669 &abs(max_theta_plot-2).gt.tol_zero) then
670 ierr = 1
671 err_msg = 'theta limits of plot need to be 0..2pi'
672 chckerr(err_msg)
673 end if
674 end if
675 if (min_rvac_plot.lt.0) then
676 call writo('The minimum value for Rvac_plot has to be &
677 &greater than zero',warning=.true.)
678 call lvl_ud(1)
679 min_rvac_plot = 0.1*r_0
680 call writo('It has been set to '//&
681 &trim(r2str(min_rvac_plot)))
682 call lvl_ud(-1)
683 end if
684 if (max_rvac_plot.lt.min_rvac_plot) then
685 ierr = 1
686 err_msg = 'The maximum value for Rvac_plot has to be &
687 &greater than the minimum value'
688 chckerr(err_msg)
689 end if
690 if (eq_style.eq.1 .and. plot_vac_pot) then
691 plot_vac_pot = .false.
692 call writo('Not plotting vacuum for VMEC yet',&
693 &warning=.true.)
694 end if
695 end if
696 end function adapt_plot
697
698 ! Checks whether variables concerning perturbation modes are correct:
699 ! the absolute value of prim_X has to be at least min_nm_X,
700 ! deduces the X style,
701 ! min_sec_X and max_sec_X, pertaining to X style 1 and n_sec_X,
702 ! pertaining to X style 2 cannot be both given,
703 ! if X style 1 (prescribed), max_sec_X has to be greater than
704 ! min_sec_X in absolute value,
705 ! if X style 2 (fast), n_mod_X has to be a number greater than 0,
706 ! for every X style, some checks are made and some variables set.
707 !> \private
708 integer function adapt_x_modes() result(ierr)
709 use x_vars, only: min_nm_x
710
711 character(*), parameter :: rout_name = 'adapt_X_modes'
712
713 ! local variables
714 character(len=max_str_ln) :: err_msg ! error message
715
716 ! initialize ierr
717 ierr = 0
718
719 ! check prim_X: common for all X styles
720 if (abs(prim_x).lt.min_nm_x) then
721 ierr = 1
722 err_msg = 'prim_X should be at least '//trim(i2str(min_nm_x))//&
723 &', and probably a bit higher still'
724 chckerr(err_msg)
725 end if
726
727 ! set X_style
728 if (min_sec_x.ne.huge(1) .and. min_sec_x.ne.huge(1)) then ! min_sec_X and max_sec_X prescribed
729 if (n_mod_x.eq.huge(1)) then ! n_mod_X not prescribed
730 x_style = 1 ! prescribed style
731 else ! n_mod_X also prescribed
732 ierr = 1
733 err_msg = 'Cannot prescribe both min and max of secondary &
734 &X and nr. of modes.'
735 chckerr(err_msg)
736 end if
737 else ! min_sec_X and max_sec_X not prescribed
738 if (n_mod_x.eq.huge(1)) then ! n_mod_X not prescribed either: default
739 x_style = 2 ! X style 2: fast
740 n_mod_x = 10 ! 10 resonating modes at every place
741 else ! but n_mod_X prescribed
742 x_style = 2
743 end if
744 end if
745
746 ! set X_grid_style if not overwritten by user
747 if (x_grid_style.ge.huge(1)) then
748 select case (x_style)
749 case (1) ! prescribed
750 x_grid_style = 1 ! equilibrium
751 case (2) ! fast
752 x_grid_style = 3 ! enriched
753 end select
754 end if
755
756 ! check X_grid style
757 if (x_grid_style.lt.1 .or. x_grid_style.gt.3) then
758 ierr = 1
759 err_msg = 'X_grid_style has to be 1 (equilibrium), 2 &
760 &(solution) or 3 (enriched)'
761 chckerr(err_msg)
762 end if
763
764 ! checks for each X style
765 select case (x_style)
766 case (1) ! prescribed
767 ! check min_sec_X, max_sec_X
768 if (abs(max_sec_x).lt.abs(min_sec_x)) then
769 ierr = 1
770 err_msg = 'max_sec_X has to be larger or equal to &
771 &min_sec_X in absolute value'
772 chckerr(err_msg)
773 end if
774 if (max_sec_x.ge.huge(1)) then
775 ierr = 1
776 err_msg = 'max_sec_X not set'
777 chckerr(err_msg)
778 end if
779 if (min_sec_x.ge.huge(1)) then
780 ierr = 1
781 err_msg = 'min_sec_X not set'
782 chckerr(err_msg)
783 end if
784 if (x_grid_style.eq.3) then
785 ierr = 1
786 err_msg = 'X_grid_style 3 (enriched) does not make &
787 &sense for X_style 1 (prescribed)'
788 chckerr(err_msg)
789 end if
790
791 ! set n_mod_X
792 n_mod_x = abs(max_sec_x)-(min_sec_x)+1
793 case (2) ! fast
794 if (n_mod_x.lt.1) then
795 ierr = 1
796 err_msg = 'the number of resonating modes has to be at &
797 &least one'
798 chckerr(err_msg)
799 end if
800 end select
801 end function adapt_x_modes
802
803 ! Checks whether min_n_par_X has been chosen correctly:
804 ! It should be larger than req_min_n_par_X (arbitrary) if an integral
805 ! is to be calculated. However, there are additional requirements for
806 ! the different magnetic integral styles: for style 1 min_par_X has to
807 ! be of the form 2+k and for style 2 of the form 4+3k.
808 !> \private
809 subroutine adapt_min_n_par_x
810 ! local variables
811 integer :: fund_n_par ! fundamental interval width
812
813 ! check req_min_n_par_X
814 if (min_n_par_x.lt.req_min_n_par_x) then
816 call writo('min_n_par_X has been increased to '//&
817 &trim(i2str(min_n_par_x))//' to have enough parallel &
818 &resolution',warning=.true.)
819 end if
820
821 ! check for individual magnetic integral style
822 select case (magn_int_style)
823 case (1) ! Trapezoidal rule
824 fund_n_par = 1
825 case (2) ! Simpson's 3/8 rule
826 fund_n_par = 3
827 end select
828 if (mod(min_n_par_x-1,fund_n_par).ne.0) then ! there is a remainder
829 min_n_par_x = 1 + fund_n_par * &
830 &((min_n_par_x-1)/fund_n_par + 1)
831 call writo('min_n_par_X has been increased to '//&
832 &trim(i2str(min_n_par_x))//' to be a of the form '//&
833 &trim(i2str(fund_n_par+1))//' + '//&
834 &trim(i2str(fund_n_par))//&
835 &'k for magnetic integral style '//&
836 &trim(i2str(magn_int_style)),warning=.true.)
837 end if
838 end subroutine adapt_min_n_par_x
839
840 ! Checks whether variables related to alpha are chosen correctly :
841 ! alpha_style has to be either 0 or 1
842 ! for 1, n_alpha = 1 and min/max_alpha = alpha
843 ! for 2, max_par_X - min_par_X = 2
844 !> \private
845 integer function adapt_alpha() result(ierr)
846 character(*), parameter :: rout_name = 'adapt_alpha'
847
848 ! local variables
849 character(len=max_str_ln) :: err_msg ! error message
850
851 ! initialize ierr
852 ierr = 0
853
854 ! check alpha style
855 if (alpha_style.lt.1 .or. alpha_style.gt.2) then
856 ierr = 1
857 err_msg = 'Alpha style needs to be 1 or 2 [def]'
858 chckerr(err_msg)
859 end if
860
861 ! adapt variables pertaining to particular alpha style
862 select case (alpha_style)
863 case (1) ! single field line, multiple turns
864 if (n_alpha.ne.1) then
865 call writo('For alpha style 1, n_alpha can only be 1',&
866 &warning=.true.)
867 n_alpha = 1
868 end if
871 case (2) ! multiple field lines, single turns
872 if (n_alpha.lt.1) then
873 ierr = 1
874 err_msg = 'For alpha style 2, n_alpha needs to be &
875 &positive'
876 chckerr(err_msg)
877 end if
878 if (abs(max_par_x - min_par_x - 2) .gt. tol_zero) then
879 min_par_x = 0._dp
880 max_par_x = 2._dp
881 call writo('min/max_par_X has been set to ['//&
882 &trim(r2strt(min_par_x))//':'//&
883 &trim(r2strt(max_par_x))//']',warning=.true.)
884 end if
885 end select
886
887 ! for HELENA, warn about axisymmetry
888 if (eq_style.eq.2) then
889 if (abs(max_par_x - min_par_x - 2) .gt. tol_zero) then
890 call writo('HELENA equilibria are axisymmetric, so only &
891 &max - min_par_X = 2 needed',warning=.true.)
892 end if
893 if (n_alpha .gt. 1) then
894 call writo('HELENA equilibria are axisymmetric, so only &
895 &n_alpha = 1 needed',warning=.true.)
896 end if
897 end if
898 end function adapt_alpha
899
900 ! Checks whether variables concerning the solution grid are correct:
901 ! min_r should be greater than 0
902 ! max_r should be lesser than 1
903 ! min_r should be lesser than max_r
904 ! The number of solution points should be large enough for the numerical
905 ! discretization scheme, and there should be enough points per process.
906 !> \private
907 integer function adapt_sol_grid(min_r,max_r,var_name) result(ierr)
908 use num_vars, only: n_procs
909
910 character(*), parameter :: rout_name = 'adapt_sol_grid'
911
912 ! input / output
913 real(dp), intent(inout) :: min_r, max_r ! min and max of norma range
914 character(len=*), intent(in) :: var_name ! name of variable
915
916 ! local variables
917 character(len=max_str_ln) :: err_msg ! error message
918
919 ! initialize ierr
920 ierr = 0
921
922 ! check min_r
923 if (min_r.lt.0._dp) then
924 min_r = 0._dp
925 call writo('min_r_'//trim(var_name)//' has been increased to '&
926 &//trim(r2strt(min_r)),warning=.true.)
927 end if
928
929 ! check if max_r is not greater than 1
930 if (max_r.gt.1._dp) then
931 max_r = 1._dp
932 call writo('max_r_'//trim(var_name)//' has been decreased to '&
933 &//trim(r2strt(max_r)),warning=.true.)
934 end if
935
936 ! check if min_r < max_r
937 if (min_r.ge.max_r) then
938 ierr = 1
939 err_msg = 'min_r_'//trim(var_name)//' should be lower than &
940 & max_r_'//trim(var_name)
941 chckerr(err_msg)
942 end if
943
944 ! check n_r_sol
945 if (n_r_sol.lt.6*norm_disc_prec_sol+2) then
947 call writo('n_r_sol has been increased to '//&
948 &trim(i2str(n_r_sol))//' to have enough points for the &
949 &normal discretization precision',warning=.true.)
950 end if
951 if (n_r_sol.lt.n_procs*(2*norm_disc_prec_sol+1)) then
953 call writo('n_r_sol has been increased to '//&
954 &trim(i2str(n_r_sol))//' to have enough points per &
955 &process',warning=.true.)
956 end if
957
958 ! warning if too close to axis for VMEC
959 if (eq_style.eq.1) then
960 if (min_r.le.tol_zero) &
961 &call writo('For VMEC, you should not go all the way to &
962 &the magnetic axis.',warning=.true.)
963 end if
964 end function adapt_sol_grid
965
966 ! Checks whether the variables concerning Richardson extrapolation are
967 ! correct:
968 ! max_it_rich has to be at least 1.
969 !> \private
970 subroutine adapt_rich
971 if (max_it_rich.lt.1) then
972 max_it_rich = 1
973 call writo('max_it_rich has been increased to 1',warning=.true.)
974 end if
975 end subroutine adapt_rich
976
977 ! Checks whether the variables concerning finding zeros are correct:
978 ! max_it_zero has to be at least 2,
979 !> \private
980 subroutine adapt_zero
981 if (max_it_zero.lt.1) then
982 max_it_zero = 2
983 call writo('max_it_zero has been increased to 2',warning=.true.)
984 end if
985 end subroutine adapt_zero
986
987 ! Checks whether tolerances are correct:
988 ! tol_norm needs to be 0..1,
989 ! tol_zero needs to be min_tol..max_tol,
990 ! tol_rich needs to be min_tol..max_tol,
991 ! tol_SLEPC needs to be min_tol..max_tol.
992 ! Also sets local tol_SLEPC.
993 !> \private
994 subroutine adapt_tols
995 ! local variables
996 integer :: id ! counter
997 integer :: max_id ! maximum of user-set id
998 real(dp), parameter :: tol_slepc_def = 1.e-5_dp ! default tol_SLEPC
999
1000 ! check tol_norm
1001 if (tol_norm.lt.0._dp) then
1002 call writo('tol_norm has been increased to 0',warning=.true.)
1003 tol_norm = 0._dp
1004 end if
1005 if (tol_norm.gt.1._dp) then
1006 call writo('tol_norm has been decreased to 1',warning=.true.)
1007 tol_norm = 1._dp
1008 end if
1009
1010 ! check tol_rich
1011 if (tol_rich.lt.min_tol) then
1012 call writo('tol_rich has been increased to '//&
1013 &trim(r2str(min_tol)),warning=.true.)
1014 tol_rich = min_tol
1015 end if
1016 if (tol_rich.gt.max_tol) then
1017 call writo('tol_rich has been decreased to '//&
1018 &trim(r2str(max_tol)),warning=.true.)
1019 tol_rich = max_tol
1020 end if
1021
1022 ! check tol_zero
1023 if (tol_zero.lt.min_tol) then
1024 call writo('tol_zero has been increased to '//&
1025 &trim(r2str(min_tol)),warning=.true.)
1026 tol_zero = min_tol
1027 end if
1028 if (tol_zero.gt.max_tol) then
1029 call writo('tol_zero has been decreased to '//&
1030 &trim(r2str(max_tol)),warning=.true.)
1031 tol_zero = max_tol
1032 end if
1033
1034 ! check and setup tol_SLEPC
1035 allocate(tol_slepc_loc(max_it_rich))
1036 ! get maximum index set by user and correct
1037 max_id = 0
1038 do id = 1,max_size_tol_slepc
1039 if (tol_slepc(id).lt.huge(1._dp)) then ! was overwritten by user
1040 max_id = id
1041 if (tol_slepc(id).lt.min_tol) then
1042 call writo('tol_SLEPC('//trim(i2str(id))//&
1043 &') has been increased to '//trim(r2str(min_tol)),&
1044 &warning=.true.)
1045 tol_slepc(id) = min_tol
1046 end if
1047 if (tol_slepc(id).gt.max_tol) then
1048 call writo('tol_SLEPC('//trim(i2str(id))//&
1049 &') has been decreased to '//trim(r2str(max_tol)),&
1050 &warning=.true.)
1051 tol_slepc(id) = max_tol
1052 end if
1053 end if
1054 end do
1055 ! set local variable
1056 if (max_id.gt.0) then ! was overwritten by user
1057 if (max_id.gt.max_it_rich) then ! too many values given
1058 call writo(trim(i2str(max_id-max_it_rich))//&
1059 &' last values discarded for tol_SLEPC',warning=.true.)
1060 tol_slepc_loc = tol_slepc(1:max_it_rich)
1061 else if (max_id.eq.max_it_rich) then ! exactly the right amount of values given
1062 tol_slepc_loc = tol_slepc(1:max_it_rich)
1063 else ! not enough values given
1064 call writo('tol_SLEPC('//trim(i2str(max_id+1))//':'//&
1065 &trim(i2str(max_it_rich))//') has been set to '//&
1066 &trim(r2str(tol_slepc_def))//' [default]',&
1067 &warning=.true.)
1068 tol_slepc_loc(1:max_id) = tol_slepc(1:max_id)
1069 tol_slepc_loc(max_id+1:max_it_rich) = tol_slepc_def
1070 end if
1071 else ! was not overwritten by user
1072 tol_slepc_loc = tol_slepc_def ! all tolerances equal to default
1073 end if
1074 end subroutine adapt_tols
1075
1076 ! Checks whether normalization variables are chosen correctly:
1077 ! rho_0 has to be positive.
1078 !> \private
1079 integer function adapt_normalization() result(ierr)
1080 character(*), parameter :: rout_name = 'adapt_normalization'
1081
1082 ! local variables
1083 character(len=max_str_ln) :: err_msg ! error message
1084
1085 ! initialize ierr
1086 ierr = 0
1087
1088 if (rho_0.le.0) then
1089 ierr = 1
1090 err_msg = 'rho_0 has to be positive'
1091 chckerr(err_msg)
1092 end if
1093 end function adapt_normalization
1094 end function read_input_opts
1095
1096 !> Reads the equilibrium input file if no Richardson restart.
1097 !!
1098 !! These variables will be passed on through the HDF5 data system.
1099 !!
1100 !! \see reconstruct_pb3d_in().
1101 !!
1102 !! \return ierr
1103 integer function read_input_eq() result(ierr)
1105 use vmec_ops, only: read_vmec
1106 use helena_ops, only: read_hel
1107 use grid_vars, only: n_r_in
1108
1109 character(*), parameter :: rout_name = 'read_input_eq'
1110
1111
1112 ! initialize ierr
1113 ierr = 0
1114
1115 ! choose which equilibrium style is being used:
1116 ! 1: VMEC
1117 ! 2: HELENA
1118 select case (eq_style)
1119 case (1) ! VMEC
1121 chckerr('')
1122 case (2) ! HELENA
1124 chckerr('')
1125 end select
1126
1127 ! close equilibrium file
1128 close(eq_i)
1129 end function read_input_eq
1130
1131 !> Print input quantities to an output file.
1132 !!
1133 !! Variables printed:
1134 !! - \c misc_in: \c prog_version, \c eq_style, \c rho_style, \c R_0,
1135 !! \c pres_0, \c B_0, \c psi_0 \c rho_0, \c T_0, \c vac_perm,
1136 !! \c use_pol_flux_F, \c use_pol_flux_E, \c use_normalization, \c
1137 !! norm_disc_prec_eq, \c norm_disc_style_sol, \c n_r_in, \c n_r_eq, \c
1138 !! n_r_sol, \c debug_version
1139 !! - \c misc_in_V: \c is_asym_V, \c is_freeb_V, \c mnmax_V, \c mpol_V, \c
1140 !! ntor_V, \c gam_V
1141 !! - \c flux_t_V
1142 !! - \c flux_p_V
1143 !! - \c pres_V
1144 !! - \c rot_t_V
1145 !! - \c q_saf_V
1146 !! - \c mn_V
1147 !! - \c R_V_c
1148 !! - \c R_V_s
1149 !! - \c Z_V_c
1150 !! - \c Z_V_s
1151 !! - \c L_V_c
1152 !! - \c L_V_s
1153 !! - \c B_V_sub: \c B_V_sub_c, \c B_V_sub_s
1154 !! - \c B_V: \c B_V_c, \c B_V_s
1155 !! - \c J_V_sup_int
1156 !! - \c jac_V: \c jac_V_c, \c jac_V_s
1157 !! - \c misc_in_H: \c ias, \c nchi
1158 !! - \c RZ_H: \c R_H, \c Z_H
1159 !! - \c chi_H
1160 !! - \c q_saf_H
1161 !! - \c rot_t_H
1162 !! - \c RBphi_H
1163 !! - \c pres_H
1164 !! - \c flux_p_H
1165 !! - \c flux_t_H
1166 !! - \c misc_X: \c prim_X, \c n_mod_X, \c min_sec_X, \c max_sec_X,
1167 !! \c norm_disc_prec_X, \c norm_style, \c U_style, \c X_style, \c
1168 !! matrix_SLEPC_style, \c K_style, \c alpha_style
1169 !! - \c misc_sol: \c min_r_sol, \c max_r_sol, \c norm_disc_prec_sol, \c
1170 !! BC_style, \c EV_BC, \c EV_BC
1171 !!
1172 !! \return ierr
1173 integer function print_output_in(data_name,remove_previous_arrs) &
1174 &result(ierr)
1175
1182 use eq_vars, only: r_0, pres_0, b_0, psi_0, rho_0, t_0, vac_perm, &
1185 &max_alpha
1186 use grid_ops, only: calc_norm_range
1188 &n_mod_x
1189 use hdf5_ops, only: print_hdf5_arrs
1190 use hdf5_vars, only: dealloc_var_1d, var_1d_type, &
1192 use helena_vars, only: chi_h, flux_p_h, flux_t_h, r_h, z_h, nchi, ias, &
1194 use vmec_vars, only: is_freeb_v, mnmax_v, mpol_v, ntor_v, is_asym_v, &
1195 &gam_v, r_v_c, r_v_s, z_v_c, z_v_s, l_v_c, l_v_s, jac_v_c, &
1196 &jac_v_s, mnmax_v, mn_v, rot_t_v, q_saf_v, pres_v, flux_t_v, &
1197 &flux_p_v, nfp_v, b_v_sub_c, b_v_sub_s
1198 use helena_vars, only: h_h_11, h_h_12, h_h_33
1199#if ldebug
1200 use vmec_vars, only: b_v_c, b_v_s, j_v_sup_int
1201#endif
1202
1203 character(*), parameter :: rout_name = 'print_output_in'
1204
1205 ! input / output
1206 character(len=*), intent(in) :: data_name !< name under which to store
1207 logical, intent(in), optional :: remove_previous_arrs !< remove previous variables if present
1208
1209 ! local variables
1210 type(var_1d_type), allocatable, target :: in_1d(:) ! 1D equivalent of input variables
1211 type(var_1d_type), pointer :: in_1d_loc => null() ! local element in in_1D
1212 integer :: id ! counter
1213 integer :: in_limits(2) ! min. and max. index of input variable grid of this process
1214 integer :: loc_size ! local size
1215
1216 ! initialize ierr
1217 ierr = 0
1218
1219 ! user output
1220 call writo('Write input variables to output file')
1221 call lvl_ud(1)
1222
1223 ! plot modes if requested and VMEC
1224 if (plot_vmec_modes) then
1225 call writo('Plotting decay of VMEC modes')
1226 call print_ex_2d('R_V_c','R_V_c',&
1227 &log10(maxval(abs(r_v_c(:,:,0)),2)),draw=.false.)
1228 call draw_ex(['R_V_c'],'R_V_c',1,1,.false.,ex_plot_style=1)
1229 call print_ex_2d('R_V_s','R_V_s',&
1230 &log10(maxval(abs(r_v_s(:,:,0)),2)),draw=.false.)
1231 call draw_ex(['R_V_s'],'R_V_s',1,1,.false.,ex_plot_style=1)
1232 call print_ex_2d('Z_V_c','Z_V_c',&
1233 &log10(maxval(abs(z_v_c(:,:,0)),2)),draw=.false.)
1234 call draw_ex(['Z_V_c'],'Z_V_c',1,1,.false.,ex_plot_style=1)
1235 call print_ex_2d('Z_V_s','Z_V_s',&
1236 &log10(maxval(abs(z_v_s(:,:,0)),2)),draw=.false.)
1237 call draw_ex(['Z_V_s'],'Z_V_s',1,1,.false.,ex_plot_style=1)
1238 call print_ex_2d('L_V_c','L_V_c',&
1239 &log10(maxval(abs(l_v_c(:,:,0)),2)),draw=.false.)
1240 call draw_ex(['L_V_c'],'L_V_c',1,1,.false.,ex_plot_style=1)
1241 call print_ex_2d('L_V_s','L_V_s',&
1242 &log10(maxval(abs(l_v_s(:,:,0)),2)),draw=.false.)
1243 call draw_ex(['L_V_s'],'L_V_s',1,1,.false.,ex_plot_style=1)
1244 call print_ex_2d('jac_V_c','jac_V_c',&
1245 &log10(maxval(abs(jac_v_c(:,:,0)),2)),draw=.false.)
1246 call draw_ex(['jac_V_c'],'jac_V_c',1,1,.false.,ex_plot_style=1)
1247 call print_ex_2d('jac_V_s','jac_V_s',&
1248 &log10(maxval(abs(jac_v_s(:,:,0)),2)),draw=.false.)
1249 call draw_ex(['jac_V_s'],'jac_V_s',1,1,.false.,ex_plot_style=1)
1250#if ldebug
1251 call print_ex_2d('B_V_c','B_V_c',&
1252 &log10(maxval(abs(b_v_c),2)),draw=.false.)
1253 call draw_ex(['B_V_c'],'B_V_c',1,1,.false.,ex_plot_style=1)
1254 call print_ex_2d('B_V_s','B_V_s',&
1255 &log10(maxval(abs(b_v_s),2)),draw=.false.)
1256 call draw_ex(['B_V_s'],'B_V_s',1,1,.false.,ex_plot_style=1)
1257#endif
1258 end if
1259
1260 ! calculate limits of input range
1261 ierr = calc_norm_range('PB3D_in',in_limits=in_limits)
1262 chckerr('')
1263 n_r_eq = in_limits(2)-in_limits(1)+1
1264 if (in_limits(1).gt.1 .or. in_limits(2).lt.n_r_in) then
1265 call writo('Only the normal range '//trim(i2str(in_limits(1)))//&
1266 &'..'//trim(i2str(in_limits(2)))//' is written')
1267 else
1268 call writo('The entire normal range '//trim(i2str(in_limits(1)))//&
1269 &'..'//trim(i2str(in_limits(2)))//' is written')
1270 end if
1271 call writo('(relevant to solution range '//&
1272 &trim(r2strt(min_r_sol))//'..'//trim(r2strt(max_r_sol))//')')
1273
1274 ! set up 1D equivalents of input variables
1275 allocate(in_1d(max_dim_var_1d))
1276
1277 ! Set up common variables in_1D
1278 id = 1
1279
1280 ! misc_in
1281 in_1d_loc => in_1d(id); id = id+1
1282 in_1d_loc%var_name = 'misc_in'
1283 allocate(in_1d_loc%tot_i_min(1),in_1d_loc%tot_i_max(1))
1284 allocate(in_1d_loc%loc_i_min(1),in_1d_loc%loc_i_max(1))
1285 in_1d_loc%loc_i_min = [1]
1286 in_1d_loc%loc_i_max = [20]
1287 in_1d_loc%tot_i_min = in_1d_loc%loc_i_min
1288 in_1d_loc%tot_i_max = in_1d_loc%loc_i_max
1289 allocate(in_1d_loc%p(20))
1290 in_1d_loc%p = [prog_version,eq_style*1._dp,rho_style*1._dp,&
1291 &r_0,pres_0,b_0,psi_0,rho_0,t_0,vac_perm,-1._dp,-1._dp,&
1292 &-1._dp,norm_disc_prec_eq*1._dp,n_r_in*1._dp,n_r_eq*1._dp,&
1293 &n_r_sol*1._dp,max_flux_e,max_flux_f,-1._dp]
1294 if (use_pol_flux_e) in_1d_loc%p(11) = 1._dp
1295 if (use_pol_flux_f) in_1d_loc%p(12) = 1._dp
1296 if (use_normalization) in_1d_loc%p(13) = 1._dp
1297 if (debug_version) in_1d_loc%p(20) = 1._dp
1298
1299 ! misc_in_V or misc_in_H, depending on equilibrium style
1300 select case (eq_style)
1301 case (1) ! VMEC
1302 ! misc_in_V
1303 in_1d_loc => in_1d(id); id = id+1
1304 in_1d_loc%var_name = 'misc_in_V'
1305 allocate(in_1d_loc%tot_i_min(1),&
1306 &in_1d_loc%tot_i_max(1))
1307 allocate(in_1d_loc%loc_i_min(1),&
1308 &in_1d_loc%loc_i_max(1))
1309 in_1d_loc%loc_i_min = [1]
1310 in_1d_loc%loc_i_max = [7]
1311 in_1d_loc%tot_i_min = in_1d_loc%loc_i_min
1312 in_1d_loc%tot_i_max = in_1d_loc%loc_i_max
1313 allocate(in_1d_loc%p(7))
1314 in_1d_loc%p = [-1._dp,-1._dp,mnmax_v*1._dp,mpol_v*1._dp,&
1315 &ntor_v*1._dp,nfp_v*1._dp,gam_v]
1316 if (is_asym_v) in_1d_loc%p(1) = 1._dp
1317 if (is_freeb_v) in_1d_loc%p(2) = 1._dp
1318
1319 loc_size = n_r_eq*size(flux_t_v,2)
1320
1321 ! flux_t_V
1322 in_1d_loc => in_1d(id); id = id+1
1323 in_1d_loc%var_name = 'flux_t_V'
1324 allocate(in_1d_loc%tot_i_min(2),in_1d_loc%tot_i_max(2))
1325 allocate(in_1d_loc%loc_i_min(2),in_1d_loc%loc_i_max(2))
1326 in_1d_loc%loc_i_min = [1,0]
1327 in_1d_loc%loc_i_max = [n_r_eq,size(flux_t_v,2)-1]
1328 in_1d_loc%tot_i_min = in_1d_loc%loc_i_min
1329 in_1d_loc%tot_i_max = in_1d_loc%loc_i_max
1330 allocate(in_1d_loc%p(loc_size))
1331 in_1d_loc%p = reshape(flux_t_v(in_limits(1):in_limits(2),:),&
1332 &[loc_size])
1333
1334 ! flux_p_V
1335 in_1d_loc => in_1d(id); id = id+1
1336 in_1d_loc%var_name = 'flux_p_V'
1337 allocate(in_1d_loc%tot_i_min(2),in_1d_loc%tot_i_max(2))
1338 allocate(in_1d_loc%loc_i_min(2),in_1d_loc%loc_i_max(2))
1339 in_1d_loc%loc_i_min = [1,0]
1340 in_1d_loc%loc_i_max = [n_r_eq,size(flux_p_v,2)-1]
1341 in_1d_loc%tot_i_min = in_1d_loc%loc_i_min
1342 in_1d_loc%tot_i_max = in_1d_loc%loc_i_max
1343 allocate(in_1d_loc%p(loc_size))
1344 in_1d_loc%p = reshape(flux_p_v(in_limits(1):in_limits(2),:),&
1345 &[loc_size])
1346
1347 loc_size = n_r_eq*size(pres_v,2)
1348
1349 ! pres_V
1350 in_1d_loc => in_1d(id); id = id+1
1351 in_1d_loc%var_name = 'pres_V'
1352 allocate(in_1d_loc%tot_i_min(2),in_1d_loc%tot_i_max(2))
1353 allocate(in_1d_loc%loc_i_min(2),in_1d_loc%loc_i_max(2))
1354 in_1d_loc%loc_i_min = [1,0]
1355 in_1d_loc%loc_i_max = [n_r_eq,size(pres_v,2)-1]
1356 in_1d_loc%tot_i_min = in_1d_loc%loc_i_min
1357 in_1d_loc%tot_i_max = in_1d_loc%loc_i_max
1358 allocate(in_1d_loc%p(loc_size))
1359 in_1d_loc%p = reshape(pres_v(in_limits(1):in_limits(2),:),&
1360 &[loc_size])
1361
1362 ! rot_t_V
1363 in_1d_loc => in_1d(id); id = id+1
1364 in_1d_loc%var_name = 'rot_t_V'
1365 allocate(in_1d_loc%tot_i_min(2),in_1d_loc%tot_i_max(2))
1366 allocate(in_1d_loc%loc_i_min(2),in_1d_loc%loc_i_max(2))
1367 in_1d_loc%loc_i_min = [1,0]
1368 in_1d_loc%loc_i_max = [n_r_eq,size(rot_t_v,2)-1]
1369 in_1d_loc%tot_i_min = in_1d_loc%loc_i_min
1370 in_1d_loc%tot_i_max = in_1d_loc%loc_i_max
1371 allocate(in_1d_loc%p(loc_size))
1372 in_1d_loc%p = reshape(rot_t_v(in_limits(1):in_limits(2),:),&
1373 &[loc_size])
1374
1375 ! q_saf_V
1376 in_1d_loc => in_1d(id); id = id+1
1377 in_1d_loc%var_name = 'q_saf_V'
1378 allocate(in_1d_loc%tot_i_min(2),in_1d_loc%tot_i_max(2))
1379 allocate(in_1d_loc%loc_i_min(2),in_1d_loc%loc_i_max(2))
1380 in_1d_loc%loc_i_min = [1,0]
1381 in_1d_loc%loc_i_max = [n_r_eq,size(q_saf_v,2)-1]
1382 in_1d_loc%tot_i_min = in_1d_loc%loc_i_min
1383 in_1d_loc%tot_i_max = in_1d_loc%loc_i_max
1384 allocate(in_1d_loc%p(loc_size))
1385 in_1d_loc%p = reshape(q_saf_v(in_limits(1):in_limits(2),:),&
1386 &[loc_size])
1387
1388 ! mn_V
1389 in_1d_loc => in_1d(id); id = id+1
1390 in_1d_loc%var_name = 'mn_V'
1391 allocate(in_1d_loc%tot_i_min(2),in_1d_loc%tot_i_max(2))
1392 allocate(in_1d_loc%loc_i_min(2),in_1d_loc%loc_i_max(2))
1393 in_1d_loc%loc_i_min = [1,1]
1394 in_1d_loc%loc_i_max = [mnmax_v,2]
1395 in_1d_loc%tot_i_min = in_1d_loc%loc_i_min
1396 in_1d_loc%tot_i_max = in_1d_loc%loc_i_max
1397 allocate(in_1d_loc%p(size(mn_v)))
1398 in_1d_loc%p = reshape(mn_v,[size(mn_v)])
1399
1400 ! RZL_V
1401 in_1d_loc => in_1d(id); id = id+1
1402 in_1d_loc%var_name = 'RZL_V'
1403 allocate(in_1d_loc%tot_i_min(4),in_1d_loc%tot_i_max(4))
1404 allocate(in_1d_loc%loc_i_min(4),in_1d_loc%loc_i_max(4))
1405 in_1d_loc%loc_i_min = [1,1,0,1]
1406 in_1d_loc%loc_i_max = [mnmax_v,n_r_eq,size(r_v_c,3)-1,6]
1407 in_1d_loc%tot_i_min = in_1d_loc%loc_i_min
1408 in_1d_loc%tot_i_max = in_1d_loc%loc_i_max
1409 allocate(in_1d_loc%p(6*mnmax_v*n_r_eq*size(r_v_c,3)))
1410 in_1d_loc%p = reshape([r_v_c(:,in_limits(1):in_limits(2),:),&
1411 &r_v_s(:,in_limits(1):in_limits(2),:),&
1412 &z_v_c(:,in_limits(1):in_limits(2),:),&
1413 &z_v_s(:,in_limits(1):in_limits(2),:),&
1414 &l_v_c(:,in_limits(1):in_limits(2),:),&
1415 &l_v_s(:,in_limits(1):in_limits(2),:)],&
1416 &[6*mnmax_v*n_r_eq*size(r_v_c,3)])
1417
1418 ! jac_V
1419 in_1d_loc => in_1d(id); id = id+1
1420 in_1d_loc%var_name = 'jac_V'
1421 allocate(in_1d_loc%tot_i_min(4),in_1d_loc%tot_i_max(4))
1422 allocate(in_1d_loc%loc_i_min(4),in_1d_loc%loc_i_max(4))
1423 in_1d_loc%loc_i_min = [1,1,0,1]
1424 in_1d_loc%loc_i_max = [mnmax_v,n_r_eq,size(jac_v_c,3)-1,2]
1425 in_1d_loc%tot_i_min = in_1d_loc%loc_i_min
1426 in_1d_loc%tot_i_max = in_1d_loc%loc_i_max
1427 allocate(in_1d_loc%p(2*mnmax_v*n_r_eq*size(jac_v_c,3)))
1428 in_1d_loc%p = reshape([jac_v_c(:,in_limits(1):in_limits(2),:),&
1429 &jac_v_s(:,in_limits(1):in_limits(2),:)],&
1430 &[2*mnmax_v*n_r_eq*size(jac_v_c,3)])
1431
1432 ! B_V_sub
1433 in_1d_loc => in_1d(id); id = id+1
1434 in_1d_loc%var_name = 'B_V_sub'
1435 allocate(in_1d_loc%tot_i_min(4),in_1d_loc%tot_i_max(4))
1436 allocate(in_1d_loc%loc_i_min(4),in_1d_loc%loc_i_max(4))
1437 in_1d_loc%loc_i_min = [1,1,1,1]
1438 in_1d_loc%loc_i_max = [mnmax_v,n_r_eq,3,2]
1439 in_1d_loc%tot_i_min = in_1d_loc%loc_i_min
1440 in_1d_loc%tot_i_max = in_1d_loc%loc_i_max
1441 allocate(in_1d_loc%p(2*mnmax_v*n_r_eq*3))
1442 in_1d_loc%p = reshape([&
1443 &b_v_sub_c(:,in_limits(1):in_limits(2),:),&
1444 &b_v_sub_s(:,in_limits(1):in_limits(2),:)],&
1445 &[2*mnmax_v*n_r_eq*3])
1446
1447#if ldebug
1448 ! B_V
1449 in_1d_loc => in_1d(id); id = id+1
1450 in_1d_loc%var_name = 'B_V'
1451 allocate(in_1d_loc%tot_i_min(3),in_1d_loc%tot_i_max(3))
1452 allocate(in_1d_loc%loc_i_min(3),in_1d_loc%loc_i_max(3))
1453 in_1d_loc%loc_i_min = [1,1,1]
1454 in_1d_loc%loc_i_max = [mnmax_v,n_r_eq,2]
1455 in_1d_loc%tot_i_min = in_1d_loc%loc_i_min
1456 in_1d_loc%tot_i_max = in_1d_loc%loc_i_max
1457 allocate(in_1d_loc%p(2*mnmax_v*n_r_eq))
1458 in_1d_loc%p = reshape([b_v_c(:,in_limits(1):in_limits(2)),&
1459 &b_v_s(:,in_limits(1):in_limits(2))],[2*mnmax_v*n_r_eq])
1460
1461 ! J_V_sup_int
1462 in_1d_loc => in_1d(id); id = id+1
1463 in_1d_loc%var_name = 'J_V_sup_int'
1464 allocate(in_1d_loc%tot_i_min(2),in_1d_loc%tot_i_max(2))
1465 allocate(in_1d_loc%loc_i_min(2),in_1d_loc%loc_i_max(2))
1466 in_1d_loc%loc_i_min = [1,1]
1467 in_1d_loc%loc_i_max = [n_r_eq,2]
1468 in_1d_loc%tot_i_min = in_1d_loc%loc_i_min
1469 in_1d_loc%tot_i_max = in_1d_loc%loc_i_max
1470 allocate(in_1d_loc%p(2*n_r_eq))
1471 in_1d_loc%p = reshape(j_v_sup_int(in_limits(1):in_limits(2),:),&
1472 &[2*n_r_eq])
1473
1474#endif
1475 case (2) ! HELENA
1476 ! misc_in_H
1477 in_1d_loc => in_1d(id); id = id+1
1478 in_1d_loc%var_name = 'misc_in_H'
1479 allocate(in_1d_loc%tot_i_min(1),&
1480 &in_1d_loc%tot_i_max(1))
1481 allocate(in_1d_loc%loc_i_min(1),&
1482 &in_1d_loc%loc_i_max(1))
1483 in_1d_loc%loc_i_min = [1]
1484 in_1d_loc%loc_i_max = [2]
1485 in_1d_loc%tot_i_min = in_1d_loc%loc_i_min
1486 in_1d_loc%tot_i_max = in_1d_loc%loc_i_max
1487 allocate(in_1d_loc%p(2))
1488 in_1d_loc%p = [ias*1._dp,nchi*1._dp]
1489
1490 ! RZ_H
1491 in_1d_loc => in_1d(id); id = id+1
1492 in_1d_loc%var_name = 'RZ_H'
1493 allocate(in_1d_loc%tot_i_min(3),in_1d_loc%tot_i_max(3))
1494 allocate(in_1d_loc%loc_i_min(3),in_1d_loc%loc_i_max(3))
1495 in_1d_loc%loc_i_min = [1,1,1]
1496 in_1d_loc%loc_i_max = [nchi,n_r_eq,2]
1497 in_1d_loc%tot_i_min = in_1d_loc%loc_i_min
1498 in_1d_loc%tot_i_max = in_1d_loc%loc_i_max
1499 allocate(in_1d_loc%p(2*nchi*n_r_eq))
1500 in_1d_loc%p = reshape([r_h(:,in_limits(1):in_limits(2)),&
1501 &z_h(:,in_limits(1):in_limits(2))],[2*nchi*n_r_eq])
1502
1503 ! chi_H
1504 in_1d_loc => in_1d(id); id = id+1
1505 in_1d_loc%var_name = 'chi_H'
1506 allocate(in_1d_loc%tot_i_min(1),in_1d_loc%tot_i_max(1))
1507 allocate(in_1d_loc%loc_i_min(1),in_1d_loc%loc_i_max(1))
1508 in_1d_loc%loc_i_min = [1]
1509 in_1d_loc%loc_i_max = [nchi]
1510 in_1d_loc%tot_i_min = in_1d_loc%loc_i_min
1511 in_1d_loc%tot_i_max = in_1d_loc%loc_i_max
1512 allocate(in_1d_loc%p(nchi))
1513 in_1d_loc%p = chi_h
1514
1515 loc_size = n_r_eq*size(flux_p_h,2)
1516
1517 ! flux_p_H
1518 in_1d_loc => in_1d(id); id = id+1
1519 in_1d_loc%var_name = 'flux_p_H'
1520 allocate(in_1d_loc%tot_i_min(2),in_1d_loc%tot_i_max(2))
1521 allocate(in_1d_loc%loc_i_min(2),in_1d_loc%loc_i_max(2))
1522 in_1d_loc%loc_i_min = [1,0]
1523 in_1d_loc%loc_i_max = [n_r_eq,size(flux_p_h,2)-1]
1524 in_1d_loc%tot_i_min = in_1d_loc%loc_i_min
1525 in_1d_loc%tot_i_max = in_1d_loc%loc_i_max
1526 allocate(in_1d_loc%p(loc_size))
1527 in_1d_loc%p = reshape(flux_p_h(in_limits(1):in_limits(2),:),&
1528 &[loc_size])
1529
1530 ! flux_t_H
1531 in_1d_loc => in_1d(id); id = id+1
1532 in_1d_loc%var_name = 'flux_t_H'
1533 allocate(in_1d_loc%tot_i_min(2),in_1d_loc%tot_i_max(2))
1534 allocate(in_1d_loc%loc_i_min(2),in_1d_loc%loc_i_max(2))
1535 in_1d_loc%loc_i_min = [1,0]
1536 in_1d_loc%loc_i_max = [n_r_eq,size(flux_t_h,2)-1]
1537 in_1d_loc%tot_i_min = in_1d_loc%loc_i_min
1538 in_1d_loc%tot_i_max = in_1d_loc%loc_i_max
1539 allocate(in_1d_loc%p(loc_size))
1540 in_1d_loc%p = reshape(flux_t_h(in_limits(1):in_limits(2),:),&
1541 &[loc_size])
1542
1543 ! q_saf_H
1544 in_1d_loc => in_1d(id); id = id+1
1545 in_1d_loc%var_name = 'q_saf_H'
1546 allocate(in_1d_loc%tot_i_min(2),in_1d_loc%tot_i_max(2))
1547 allocate(in_1d_loc%loc_i_min(2),in_1d_loc%loc_i_max(2))
1548 in_1d_loc%loc_i_min = [1,0]
1549 in_1d_loc%loc_i_max = [n_r_eq,size(q_saf_h,2)-1]
1550 in_1d_loc%tot_i_min = in_1d_loc%loc_i_min
1551 in_1d_loc%tot_i_max = in_1d_loc%loc_i_max
1552 allocate(in_1d_loc%p(loc_size))
1553 in_1d_loc%p = reshape(q_saf_h(in_limits(1):in_limits(2),:),&
1554 &[loc_size])
1555
1556 ! rot_t_H
1557 in_1d_loc => in_1d(id); id = id+1
1558 in_1d_loc%var_name = 'rot_t_H'
1559 allocate(in_1d_loc%tot_i_min(2),in_1d_loc%tot_i_max(2))
1560 allocate(in_1d_loc%loc_i_min(2),in_1d_loc%loc_i_max(2))
1561 in_1d_loc%loc_i_min = [1,0]
1562 in_1d_loc%loc_i_max = [n_r_eq,size(rot_t_h,2)-1]
1563 in_1d_loc%tot_i_min = in_1d_loc%loc_i_min
1564 in_1d_loc%tot_i_max = in_1d_loc%loc_i_max
1565 allocate(in_1d_loc%p(loc_size))
1566 in_1d_loc%p = reshape(rot_t_h(in_limits(1):in_limits(2),:),&
1567 &[loc_size])
1568
1569 ! pres_H
1570 in_1d_loc => in_1d(id); id = id+1
1571 in_1d_loc%var_name = 'pres_H'
1572 allocate(in_1d_loc%tot_i_min(2),in_1d_loc%tot_i_max(2))
1573 allocate(in_1d_loc%loc_i_min(2),in_1d_loc%loc_i_max(2))
1574 in_1d_loc%loc_i_min = [1,0]
1575 in_1d_loc%loc_i_max = [n_r_eq,size(pres_h,2)-1]
1576 in_1d_loc%tot_i_min = in_1d_loc%loc_i_min
1577 in_1d_loc%tot_i_max = in_1d_loc%loc_i_max
1578 allocate(in_1d_loc%p(loc_size))
1579 in_1d_loc%p = reshape(pres_h(in_limits(1):in_limits(2),:),&
1580 &[loc_size])
1581
1582 ! RBphi_H
1583 in_1d_loc => in_1d(id); id = id+1
1584 in_1d_loc%var_name = 'RBphi_H'
1585 allocate(in_1d_loc%tot_i_min(2),in_1d_loc%tot_i_max(2))
1586 allocate(in_1d_loc%loc_i_min(2),in_1d_loc%loc_i_max(2))
1587 in_1d_loc%loc_i_min = [1,0]
1588 in_1d_loc%loc_i_max = [n_r_eq,size(rbphi_h,2)-1]
1589 in_1d_loc%tot_i_min = in_1d_loc%loc_i_min
1590 in_1d_loc%tot_i_max = in_1d_loc%loc_i_max
1591 allocate(in_1d_loc%p(loc_size))
1592 in_1d_loc%p = reshape(rbphi_h(in_limits(1):in_limits(2),:),&
1593 &[loc_size])
1594
1595 ! h_H
1596 in_1d_loc => in_1d(id); id = id+1
1597 in_1d_loc%var_name = 'h_H'
1598 allocate(in_1d_loc%tot_i_min(3),in_1d_loc%tot_i_max(3))
1599 allocate(in_1d_loc%loc_i_min(3),in_1d_loc%loc_i_max(3))
1600 in_1d_loc%loc_i_min = [1,1,1]
1601 in_1d_loc%loc_i_max = [nchi,n_r_eq,3]
1602 in_1d_loc%tot_i_min = in_1d_loc%loc_i_min
1603 in_1d_loc%tot_i_max = in_1d_loc%loc_i_max
1604 allocate(in_1d_loc%p(3*nchi*n_r_eq))
1605 in_1d_loc%p = reshape([h_h_11(:,in_limits(1):in_limits(2)),&
1606 &h_h_12(:,in_limits(1):in_limits(2)),&
1607 &h_h_33(:,in_limits(1):in_limits(2))],&
1608 &[3*nchi*n_r_eq])
1609 end select
1610
1611 ! misc_X
1612 in_1d_loc => in_1d(id); id = id+1
1613 in_1d_loc%var_name = 'misc_X'
1614 allocate(in_1d_loc%tot_i_min(1),in_1d_loc%tot_i_max(1))
1615 allocate(in_1d_loc%loc_i_min(1),in_1d_loc%loc_i_max(1))
1616 in_1d_loc%loc_i_min = [1]
1617 in_1d_loc%loc_i_max = [17]
1618 in_1d_loc%tot_i_min = in_1d_loc%loc_i_min
1619 in_1d_loc%tot_i_max = in_1d_loc%loc_i_max
1620 allocate(in_1d_loc%p(17))
1621 in_1d_loc%p = [prim_x*1._dp,n_mod_x*1._dp,min_sec_x*1._dp,&
1622 &max_sec_x*1._dp,norm_disc_prec_x*1._dp,norm_style*1._dp,&
1623 &u_style*1._dp,x_style*1._dp,matrix_slepc_style*1._dp,&
1624 &magn_int_style*1._dp,k_style*1._dp,alpha_style*1._dp,&
1627
1628 ! misc_sol
1629 in_1d_loc => in_1d(id); id = id+1
1630 in_1d_loc%var_name = 'misc_sol'
1631 allocate(in_1d_loc%tot_i_min(1),in_1d_loc%tot_i_max(1))
1632 allocate(in_1d_loc%loc_i_min(1),in_1d_loc%loc_i_max(1))
1633 in_1d_loc%loc_i_min = [1]
1634 in_1d_loc%loc_i_max = [8]
1635 in_1d_loc%tot_i_min = in_1d_loc%loc_i_min
1636 in_1d_loc%tot_i_max = in_1d_loc%loc_i_max
1637 allocate(in_1d_loc%p(8))
1638 in_1d_loc%p = [min_r_sol,max_r_sol,norm_disc_prec_sol*1._dp,&
1639 &norm_disc_style_sol*1._dp,bc_style(1)*1._dp,bc_style(2)*1._dp,&
1640 &ev_style*1._dp,ev_bc]
1641
1642 ! write
1643 ierr = print_hdf5_arrs(in_1d(1:id-1),pb3d_name,trim(data_name),&
1644 &remove_previous_arrs=remove_previous_arrs)
1645 chckerr('')
1646
1647 ! clean up
1648 call dealloc_var_1d(in_1d)
1649 nullify(in_1d_loc)
1650
1651 ! user output
1652 call lvl_ud(-1)
1653 end function print_output_in
1654end module input_ops
Deallocates 1D variable.
Definition HDF5_vars.f90:68
Print 2-D output on a file.
Variables that have to do with equilibrium quantities and the grid used in the calculations:
Definition eq_vars.f90:27
real(dp), public r_0
independent normalization constant for nondimensionalization
Definition eq_vars.f90:42
real(dp), public psi_0
derived normalization constant for nondimensionalization
Definition eq_vars.f90:46
real(dp), public max_flux_f
max. flux in Flux coordinates, set in calc_norm_range_PB3D_in
Definition eq_vars.f90:50
real(dp), public t_0
derived normalization constant for nondimensionalization
Definition eq_vars.f90:47
real(dp), public max_flux_e
max. flux in Equilibrium coordinates, set in calc_norm_range_PB3D_in
Definition eq_vars.f90:49
real(dp), public rho_0
independent normalization constant for nondimensionalization
Definition eq_vars.f90:44
real(dp), public pres_0
independent normalization constant for nondimensionalization
Definition eq_vars.f90:43
real(dp), public vac_perm
either usual mu_0 (default) or normalized
Definition eq_vars.f90:48
real(dp), public b_0
derived normalization constant for nondimensionalization
Definition eq_vars.f90:45
Operations that have to do with the grids and different coordinate systems.
Definition grid_ops.f90:4
integer function, public calc_norm_range(style, in_limits, eq_limits, x_limits, sol_limits, r_f_eq, r_f_x, r_f_sol, jq)
Calculates normal range for the input grid, the equilibrium grid and/or the solution grid.
Definition grid_ops.f90:46
Variables pertaining to the different grids used.
Definition grid_vars.f90:4
real(dp), dimension(:), allocatable, public alpha
field line label alpha
Definition grid_vars.f90:28
real(dp), public min_par_x
min. of parallel coordinate [ ] in field-aligned grid
Definition grid_vars.f90:24
integer, public n_alpha
nr. of field-lines
Definition grid_vars.f90:23
integer, public n_r_eq
nr. of normal points in equilibrium grid
Definition grid_vars.f90:20
real(dp), public max_par_x
max. of parallel coordinate [ ] in field-aligned grid
Definition grid_vars.f90:25
integer, public n_r_in
nr. of normal points in input grid
Definition grid_vars.f90:19
real(dp), public max_alpha
max. of field-line label [ ] in field-aligned grid
Definition grid_vars.f90:27
real(dp), public min_alpha
min. of field-line label [ ] in field-aligned grid
Definition grid_vars.f90:26
integer, public n_r_sol
nr. of normal points in solution grid
Definition grid_vars.f90:22
Operations on HDF5 and XDMF variables.
Definition HDF5_ops.f90:27
integer function, public print_hdf5_arrs(vars, pb3d_name, head_name, rich_lvl, disp_info, ind_print, remove_previous_arrs)
Prints a series of arrays, in the form of an array of pointers, to an HDF5 file.
Variables pertaining to HDF5 and XDMF.
Definition HDF5_vars.f90:4
integer, parameter, public max_dim_var_1d
maximum dimension of var_1D
Definition HDF5_vars.f90:21
Operations on HELENA variables.
Definition HELENA_ops.f90:4
integer function, public read_hel(n_r_in, use_pol_flux_h)
Reads the HELENA equilibrium data.
Variables that have to do with HELENA quantities.
integer, public ias
0 if top-bottom symmetric, 1 if not
integer, public nchi
nr. of poloidal points
real(dp), dimension(:,:), allocatable, public r_h
major radius (xout)
real(dp), dimension(:,:), allocatable, public h_h_33
upper metric factor (1 / gem12)
real(dp), dimension(:,:), allocatable, public rbphi_h
real(dp), dimension(:,:), allocatable, public pres_h
pressure profile
real(dp), dimension(:,:), allocatable, public z_h
height (yout)
real(dp), dimension(:,:), allocatable, public q_saf_h
safety factor
real(dp), dimension(:), allocatable, public chi_h
poloidal angle
real(dp), dimension(:,:), allocatable, public h_h_12
upper metric factor (gem12)
real(dp), dimension(:,:), allocatable, public flux_p_h
poloidal flux
real(dp), dimension(:,:), allocatable, public h_h_11
upper metric factor (gem11)
real(dp), dimension(:,:), allocatable, public rot_t_h
rotational transform
real(dp), dimension(:,:), allocatable, public flux_t_h
toroidal flux
Operations concerning giving input.
Definition input_ops.f90:4
integer function, public read_input_opts()
Reads input options from user-provided input file.
Definition input_ops.f90:23
integer function, public read_input_eq()
Reads the equilibrium input file if no Richardson restart.
integer function, public print_output_in(data_name, remove_previous_arrs)
Print input quantities to an output file.
Numerical utilities related to giving output.
Definition messages.f90:4
subroutine, public lvl_ud(inc)
Increases/decreases lvl of output.
Definition messages.f90:254
subroutine, public writo(input_str, persistent, error, warning, alert)
Write output to file identified by output_i.
Definition messages.f90:275
Numerical variables used by most other modules.
Definition num_vars.f90:4
integer, public sol_n_procs
nr. of MPI processes for solution with SLEPC
Definition num_vars.f90:70
integer, parameter, public dp
double precision
Definition num_vars.f90:46
integer, public post_style
style for POST (1: extended grid, 2: B-aligned grid)
Definition num_vars.f90:98
logical, public plot_b
whether to plot the magnetic field in real coordinates
Definition num_vars.f90:104
integer, public n_theta_plot
nr. of poloidal points in plot
Definition num_vars.f90:156
integer, public plot_grid_style
style for POST plot grid (0: 3-D plots, 1: slab plots, 2: slab plots with folding,...
Definition num_vars.f90:169
integer, public norm_disc_prec_sol
precision for normal discretization for solution
Definition num_vars.f90:122
logical, public plot_kappa
whether to plot curvature
Definition num_vars.f90:107
character(len=max_str_ln), public input_name
will hold the full name of the input file
Definition num_vars.f90:174
real(dp), public max_r_plot
max. of r_plot
Definition num_vars.f90:164
real(dp), public min_theta_plot
min. of theta_plot
Definition num_vars.f90:159
real(dp), public max_zeta_plot
max. of zeta_plot
Definition num_vars.f90:162
logical, public compare_tor_pos
compare quantities at toroidal positions (only for POST)
Definition num_vars.f90:151
real(dp), public max_njq_change
maximum change of prim. mode number times saf. fac. / rot. transf. when using X_style 2 (fast)
Definition num_vars.f90:119
real(dp), public tol_norm
tolerance for normal range (normalized to 0..1)
Definition num_vars.f90:134
integer, public n_procs
nr. of MPI processes
Definition num_vars.f90:69
logical, public plot_flux_q
whether to plot flux quantities in real coordinates
Definition num_vars.f90:106
integer, dimension(2), public plot_size
size of plot in inches
Definition num_vars.f90:172
integer, public n_sol_requested
how many solutions requested
Definition num_vars.f90:170
integer, parameter, public max_str_ln
maximum length of strings
Definition num_vars.f90:50
logical, public plot_sol_q
whether to plot magnetic perturbation of solution in POST
Definition num_vars.f90:109
real(dp), public ev_guess
first guess for eigenvalue
Definition num_vars.f90:117
real(dp), public ev_bc
value of artificial Eigenvalue for boundary condition
Definition num_vars.f90:116
logical, public post_output_sol
POST has outputs of solution.
Definition num_vars.f90:147
integer, public prog_style
program style (1: PB3D, 2: PB3D_POST)
Definition num_vars.f90:53
integer, dimension(2), public n_vac_plot
nr. of points in R and Z in vacuum
Definition num_vars.f90:158
logical, public use_normalization
whether to use normalization or not
Definition num_vars.f90:115
integer, public norm_disc_style_sol
style for normal discretization for solution (1: central fin. diff., 2: left fin. diff....
Definition num_vars.f90:123
integer, public u_style
style for calculation of U (1: ord.2, 2: ord.1, 1: ord.0)
Definition num_vars.f90:91
real(dp), public max_zvac_plot
max. of Z in which to plot vacuum
Definition num_vars.f90:168
logical, public debug_version
debug version used
Definition num_vars.f90:62
logical, public plot_vac_pot
whether to plot vacuum potential in POST
Definition num_vars.f90:110
integer, public n_zeta_plot
nr. of toroidal points in plot
Definition num_vars.f90:157
real(dp), public max_rvac_plot
max. of R in which to plot vacuum
Definition num_vars.f90:166
logical, public plot_resonance
whether to plot the q-profile or iota-profile with resonances
Definition num_vars.f90:102
integer, public x_grid_style
style for normal component of X grid (1: eq, 2: sol, 3: enriched)
Definition num_vars.f90:99
real(dp), public min_zeta_plot
min. of zeta_plot
Definition num_vars.f90:161
integer, public rho_style
style for equilibrium density profile
Definition num_vars.f90:90
real(dp), parameter, public prog_version
version number
Definition num_vars.f90:59
integer, parameter, public eq_i
file number of equilibrium file from VMEC or HELENA
Definition num_vars.f90:183
logical, public post_output_full
POST has output on full grids.
Definition num_vars.f90:146
integer, public eq_style
either 1 (VMEC) or 2 (HELENA)
Definition num_vars.f90:89
integer, public max_it_slepc
maximum nr. of iterations for SLEPC
Definition num_vars.f90:101
logical, public plot_vmec_modes
plot VMEC modes
Definition num_vars.f90:143
character(len=max_str_ln), public pb3d_name
name of PB3D output file
Definition num_vars.f90:139
integer, public norm_disc_prec_eq
precision for normal discretization for equilibrium
Definition num_vars.f90:120
integer, dimension(4), public n_sol_plotted
how many solutions to be plot (first unstable, last unstable, first stable, last stable)
Definition num_vars.f90:171
real(dp), public min_rvac_plot
min. of R in which to plot vacuum
Definition num_vars.f90:165
real(dp), public tol_rich
tolerance for Richardson extrapolation
Definition num_vars.f90:128
integer, public alpha_style
style for alpha (1: one field line, many turns, 2: many field lines, one turn)
Definition num_vars.f90:100
real(dp), public pert_mult_factor_post
factor with which to multiply perturbation strength for POST
Definition num_vars.f90:176
logical, public plot_j
whether to plot the current in real coordinates
Definition num_vars.f90:105
integer, public max_nr_backtracks_hh
maximum number of backtracks for Householder, relax. factors
Definition num_vars.f90:132
integer, dimension(2), public bc_style
style for BC left and right
Definition num_vars.f90:94
integer, public ex_plot_style
external plot style (1: GNUPlot, 2: Bokeh for 2D, Mayavi for 3D)
Definition num_vars.f90:175
integer, public ev_style
determines the method used for solving an EV problem
Definition num_vars.f90:88
real(dp), public min_r_plot
min. of r_plot
Definition num_vars.f90:163
integer, public rank
MPI rank.
Definition num_vars.f90:68
real(dp), public min_zvac_plot
min. of Z in which to plot vacuum
Definition num_vars.f90:167
integer, parameter, public input_i
file number of input file
Definition num_vars.f90:182
integer, public rich_restart_lvl
starting Richardson level (0: none [default])
Definition num_vars.f90:173
integer, public k_style
style for kinetic energy
Definition num_vars.f90:93
logical, public plot_sol_xi
whether to plot plasma perturbation of solution in POST
Definition num_vars.f90:108
integer, public norm_disc_prec_x
precision for normal discretization for perturbation
Definition num_vars.f90:121
real(dp), dimension(:), allocatable, public tol_slepc
tolerance for SLEPC for different Richardson levels
Definition num_vars.f90:118
integer, public norm_style
style for normalization
Definition num_vars.f90:92
integer, public max_it_zero
maximum number of iterations to find zeros
Definition num_vars.f90:131
logical, public plot_magn_grid
whether to plot the grid in real coordinates
Definition num_vars.f90:103
integer, public x_style
style for secondary mode numbers (1: prescribed, 2: fast)
Definition num_vars.f90:95
logical, public retain_all_sol
retain also faulty solutions
Definition num_vars.f90:152
integer, public matrix_slepc_style
style for matrix storage (1: sparse, 2: shell)
Definition num_vars.f90:96
real(dp), public max_theta_plot
max. of theta_plot
Definition num_vars.f90:160
logical, public use_pol_flux_e
whether poloidal flux is used in E coords.
Definition num_vars.f90:113
logical, public use_pol_flux_f
whether poloidal flux is used in F coords.
Definition num_vars.f90:114
real(dp), public tol_zero
tolerance for zeros
Definition num_vars.f90:133
real(dp), public max_tot_mem
maximum total memory for all processes [MB]
Definition num_vars.f90:74
logical, public plot_e_rec
whether to plot energy reconstruction in POST
Definition num_vars.f90:111
integer, public magn_int_style
style for magnetic integrals (1: trapezoidal, 2: Simpson 3/8)
Definition num_vars.f90:124
integer, public max_it_rich
number of levels for Richardson extrapolation
Definition num_vars.f90:127
integer, public solver_slepc_style
style for solver (1: Krylov-Schur, 2: GD)
Definition num_vars.f90:97
Operations concerning giving output, on the screen as well as in output files.
Definition output_ops.f90:5
subroutine, public draw_ex(var_names, draw_name, nplt, draw_dim, plot_on_screen, ex_plot_style, data_name, draw_ops, extra_ops, is_animated, ranges, delay, persistent)
Use external program to draw a plot.
Operations concerning Richardson extrapolation.
Definition rich_ops.f90:4
integer function, public find_max_rich_lvl(lvl_rich_req, max_lvl_rich_file)
Probe to find out which Richardson levels are available.
Definition rich_ops.f90:539
Variables concerning Richardson extrapolation.
Definition rich_vars.f90:4
integer, public min_n_par_x
min. of n_par_X (e.g. first value in Richardson loop)
Definition rich_vars.f90:22
integer, public rich_lvl
current level of Richardson extrapolation
Definition rich_vars.f90:19
integer, parameter, public req_min_n_par_x
required minimum n_par_X
Definition rich_vars.f90:21
Operations on strings.
elemental character(len=max_str_ln) function, public i2str(k)
Convert an integer to string.
elemental character(len=max_str_ln) function, public r2str(k)
Convert a real (double) to string.
elemental character(len=max_str_ln) function, public r2strt(k)
Convert a real (double) to string.
Operations that concern the output of VMEC.
Definition VMEC_ops.f90:4
integer function, public read_vmec(n_r_in, use_pol_flux_v)
Reads the VMEC equilibrium data.
Definition VMEC_ops.f90:34
Variables that concern the output of VMEC.
Definition VMEC_vars.f90:4
real(dp), dimension(:,:), allocatable, public q_saf_v
safety factor
Definition VMEC_vars.f90:38
real(dp), dimension(:,:,:), allocatable, public jac_v_c
Coeff. of in sine series (HM and FM) and norm. deriv.
Definition VMEC_vars.f90:45
real(dp), dimension(:,:,:), allocatable, public b_v_sub_c
Coeff. of B_i in sine series (r,theta,phi) (FM).
Definition VMEC_vars.f90:47
integer, dimension(:,:), allocatable, public mn_v
m and n of modes
Definition VMEC_vars.f90:32
real(dp), dimension(:,:,:), allocatable, public l_v_s
Coeff. of in cosine series (HM) and norm. deriv.
Definition VMEC_vars.f90:44
real(dp), dimension(:,:,:), allocatable, public z_v_c
Coeff. of in sine series (FM) and norm. deriv.
Definition VMEC_vars.f90:41
real(dp), dimension(:,:), allocatable, public rot_t_v
rotational transform
Definition VMEC_vars.f90:37
real(dp), dimension(:,:,:), allocatable, public jac_v_s
Coeff. of in cosine series (HM and FM) and norm. deriv.
Definition VMEC_vars.f90:46
real(dp), dimension(:,:,:), allocatable, public r_v_c
Coeff. of in sine series (FM) and norm. deriv.
Definition VMEC_vars.f90:39
real(dp), dimension(:,:), allocatable, public pres_v
pressure
Definition VMEC_vars.f90:36
real(dp), dimension(:,:), allocatable, public j_v_sup_int
Integrated poloidal and toroidal current (FM).
Definition VMEC_vars.f90:52
real(dp), dimension(:,:,:), allocatable, public b_v_sub_s
Coeff. of B_i in cosine series (r,theta,phi) (FM).
Definition VMEC_vars.f90:48
real(dp), dimension(:,:), allocatable, public b_v_s
Coeff. of magnitude of B in cosine series (HM and FM).
Definition VMEC_vars.f90:51
real(dp), dimension(:,:), allocatable, public flux_t_v
toroidal flux
Definition VMEC_vars.f90:34
real(dp), dimension(:,:,:), allocatable, public z_v_s
Coeff. of in cosine series (FM) and norm. deriv.
Definition VMEC_vars.f90:42
real(dp), dimension(:,:,:), allocatable, public r_v_s
Coeff. of in cosine series (FM) and norm. deriv.
Definition VMEC_vars.f90:40
real(dp), dimension(:,:,:), allocatable, public l_v_c
Coeff. of in sine series (HM) and norm. deriv.
Definition VMEC_vars.f90:43
real(dp), dimension(:,:), allocatable, public b_v_c
Coeff. of magnitude of B in sine series (HM and FM).
Definition VMEC_vars.f90:50
real(dp), dimension(:,:), allocatable, public flux_p_v
poloidal flux
Definition VMEC_vars.f90:35
Variables pertaining to the perturbation quantities.
Definition X_vars.f90:4
integer, public min_sec_x
m_X (pol. flux) or n_X (tor. flux) (only for X style 1)
Definition X_vars.f90:127
real(dp), public max_r_sol
max. normal range for pert.
Definition X_vars.f90:136
integer, public n_mod_x
size of m_X (pol. flux) or n_X (tor. flux)
Definition X_vars.f90:129
real(dp), public min_r_sol
min. normal range for pert.
Definition X_vars.f90:135
integer, public min_nm_x
minimum for the high-n theory (debable)
Definition X_vars.f90:130
integer, public max_sec_x
m_X (pol. flux) or n_X (tor. flux) (only for\ c X style 1)
Definition X_vars.f90:128
integer, public prim_x
n_X (pol. flux) or m_X (tor. flux)
Definition X_vars.f90:126
1D equivalent of multidimensional variables, used for internal HDF5 storage.
Definition HDF5_vars.f90:48