5#include <PB3D_macros.h>
49 character(*),
parameter :: rout_name =
'read_input_opts'
52 character(len=max_str_ln) :: err_msg
53 integer :: pb3d_rich_lvl
54 integer,
parameter :: max_size_tol_slepc = 100
57 real(
dp),
parameter :: min_tol = 1.e-12_dp
58 real(
dp),
parameter :: max_tol = 1.e-3_dp
94 call writo(
'Initialize all the inputs with default values')
127 call default_input_pb3d()
129 call default_input_post()
138 read(unit=
input_i,nml=inputdata_pb3d,iostat=ierr,&
143 call writo(
'Overwriting with user-provided file "'&
146 call writo(
'Checking user-provided file')
154 ierr = adapt_run_pb3d()
162 ierr = adapt_x_modes()
166 call adapt_min_n_par_x
177 ierr = adapt_normalization()
184 ierr = adapt_inoutput_pb3d()
196 call writo(
'Cannot read user-provided file "'&
200 call writo(
'"'//trim(err_msg)//
'"')
205 read(unit=
input_i,nml=inputdata_post,iostat=ierr,&
210 call writo(
'Overwriting with user-provided file "'&
213 call writo(
'Checking user-provided file')
218 ierr = adapt_run_post()
230 ierr = adapt_inoutput_post()
236 call writo(
'Cannot read user-provided file "'&
240 call writo(
'"'//trim(err_msg)//
'"')
247 call writo(
'Close input file')
251 call writo(
'Input values set')
256 subroutine default_input_pb3d
334 end subroutine default_input_pb3d
338 subroutine default_input_post()
355 write(*,*)
'!!!!!! min and max need to take into account normalization'
364 pb3d_rich_lvl = huge(1)
368 end subroutine default_input_post
374 subroutine adapt_mpi()
378 call writo(
'Using the maximum number of MPI processes '//&
384 &
' MPI processes for SLEPC',warning=.true.)
387 end subroutine adapt_mpi
404 integer function adapt_run_pb3d()
result(ierr)
407 character(*),
parameter :: rout_name =
'adapt_run_PB3D'
414 err_msg =
'rho_style has to be 1 (constant)'
419 call writo(
'normalization is always used for HELENA',&
424 err_msg =
'matrix_SLEPC_style has to be 1 (sparse) or 2 (shell)'
429 call writo(
'max_it_SLEPC has to be positive and is set to 1',&
434 err_msg =
'magn_int_style has to be 1 (trapezoidal) or 2 &
440 call writo(
'can only use poloidal flux for HELENA',&
445 err_msg =
'can only use norm_disc_prec_eq = 3 currently'
450 err_msg =
'norm_disc_prec_X has to be between 1 and 3'
455 err_msg =
'norm_disc_prec_sol has to be between 1 and 3'
458 end function adapt_run_pb3d
464 integer function adapt_run_post()
result(ierr)
465 character(*),
parameter :: rout_name =
'adapt_run_POST'
470 if (post_style.lt.1 .or. post_style.gt.2)
then
472 err_msg =
'POST_style has to be 1 (extended grid) or 2 &
473 &(field-aligned grid)'
476 end function adapt_run_post
483 integer function adapt_inoutput_post()
result(ierr)
486 character(*),
parameter :: rout_name =
'adapt_inoutput_POST'
489 integer :: max_pb3d_rich_lvl
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.)
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
514 call writo(
'Less than maximum Richardson level: '//&
515 &trim(
i2str(max_pb3d_rich_lvl)),warning=.true.)
518 else if (pb3d_rich_lvl.gt.max_pb3d_rich_lvl)
then
520 err_msg =
'PB3D_rich_lvl cannot be higher than '//&
521 &trim(
i2str(max_pb3d_rich_lvl))
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.)
540 end function adapt_inoutput_post
547 integer function adapt_inoutput_pb3d()
result(ierr)
548 character(*),
parameter :: rout_name =
'adapt_inoutput_PB3D'
553 if (n_sol_requested.lt.1)
then
555 call writo(
'n_sol_requested has been increased to '&
556 &//trim(
i2str(n_sol_requested)),warning=.true.)
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',&
565 if (rich_restart_lvl.gt.max_it_rich)
then
567 err_msg =
'rich_restart_lvl not within 1..max_it_rich = 1..'//&
568 &trim(
i2str(max_it_rich))
571 if (rich_restart_lvl.gt.1)
then
574 if (rich_restart_lvl.gt.pb3d_rich_lvl+1)
then
576 if (pb3d_rich_lvl.gt.0)
then
577 err_msg =
'The highest Richardson level found was '//&
578 &trim(
i2str(pb3d_rich_lvl))
580 err_msg =
'No Richardson level found'
585 end function adapt_inoutput_pb3d
597 integer function adapt_plot()
result(ierr)
600 character(*),
parameter :: rout_name =
'adapt_plot'
603 real(
dp),
parameter :: tol = 1.e-6_dp
604 character(len=max_str_ln) :: err_msg
611 call writo(
'n_theta_plot has to be positive and is set to '//&
616 call writo(
'n_zeta_plot has to be positive and is set to '//&
621 call writo(
'For n_theta_plot = 1, max_theta_plot is changed to &
628 call writo(
'For n_zeta_plot = 1, min_zeta_plot and &
629 &max_zeta_plot are changed to their average = '//&
636 call writo(
'external plot style should be')
638 call writo(
'1: GNUPlot')
639 call writo(
'2: Bokeh (2-D) / Mayavi (3-D)')
641 call writo(
'Default (1: GNUPlot) chosen',warning=.true.)
648 err_msg =
'For compare_tor_pos, you can only use &
649 &n_zeta = 3 (one point in the middle)'
654 err_msg =
'For compare_tor_pos, you can only use &
655 &POST_style = 1 (extended grid)'
660 err_msg =
'min_r_plot has to be greater than 0.'
665 err_msg =
'max_r_plot has to be 1.'
671 err_msg =
'theta limits of plot need to be 0..2pi'
676 call writo(
'The minimum value for Rvac_plot has to be &
677 &greater than zero',warning=.true.)
680 call writo(
'It has been set to '//&
686 err_msg =
'The maximum value for Rvac_plot has to be &
687 &greater than the minimum value'
692 call writo(
'Not plotting vacuum for VMEC yet',&
696 end function adapt_plot
708 integer function adapt_x_modes()
result(ierr)
711 character(*),
parameter :: rout_name =
'adapt_X_modes'
714 character(len=max_str_ln) :: err_msg
722 err_msg =
'prim_X should be at least '//trim(
i2str(
min_nm_x))//&
723 &
', and probably a bit higher still'
733 err_msg =
'Cannot prescribe both min and max of secondary &
734 &X and nr. of modes.'
747 if (x_grid_style.ge.huge(1))
then
748 select case (x_style)
757 if (x_grid_style.lt.1 .or. x_grid_style.gt.3)
then
759 err_msg =
'X_grid_style has to be 1 (equilibrium), 2 &
760 &(solution) or 3 (enriched)'
765 select case (x_style)
770 err_msg =
'max_sec_X has to be larger or equal to &
771 &min_sec_X in absolute value'
776 err_msg =
'max_sec_X not set'
781 err_msg =
'min_sec_X not set'
784 if (x_grid_style.eq.3)
then
786 err_msg =
'X_grid_style 3 (enriched) does not make &
787 &sense for X_style 1 (prescribed)'
796 err_msg =
'the number of resonating modes has to be at &
801 end function adapt_x_modes
809 subroutine adapt_min_n_par_x
811 integer :: fund_n_par
816 call writo(
'min_n_par_X has been increased to '//&
818 &resolution',warning=.true.)
822 select case (magn_int_style)
831 call writo(
'min_n_par_X has been increased to '//&
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.)
838 end subroutine adapt_min_n_par_x
845 integer function adapt_alpha()
result(ierr)
846 character(*),
parameter :: rout_name =
'adapt_alpha'
849 character(len=max_str_ln) :: err_msg
855 if (alpha_style.lt.1 .or. alpha_style.gt.2)
then
857 err_msg =
'Alpha style needs to be 1 or 2 [def]'
862 select case (alpha_style)
865 call writo(
'For alpha style 1, n_alpha can only be 1',&
874 err_msg =
'For alpha style 2, n_alpha needs to be &
881 call writo(
'min/max_par_X has been set to ['//&
888 if (eq_style.eq.2)
then
890 call writo(
'HELENA equilibria are axisymmetric, so only &
891 &max - min_par_X = 2 needed',warning=.true.)
894 call writo(
'HELENA equilibria are axisymmetric, so only &
895 &n_alpha = 1 needed',warning=.true.)
898 end function adapt_alpha
907 integer function adapt_sol_grid(min_r,max_r,var_name)
result(ierr)
910 character(*),
parameter :: rout_name =
'adapt_sol_grid'
913 real(
dp),
intent(inout) :: min_r, max_r
914 character(len=*),
intent(in) :: var_name
917 character(len=max_str_ln) :: err_msg
923 if (min_r.lt.0._dp)
then
925 call writo(
'min_r_'//trim(var_name)//
' has been increased to '&
926 &//trim(
r2strt(min_r)),warning=.true.)
930 if (max_r.gt.1._dp)
then
932 call writo(
'max_r_'//trim(var_name)//
' has been decreased to '&
933 &//trim(
r2strt(max_r)),warning=.true.)
937 if (min_r.ge.max_r)
then
939 err_msg =
'min_r_'//trim(var_name)//
' should be lower than &
940 & max_r_'//trim(var_name)
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.)
953 call writo(
'n_r_sol has been increased to '//&
955 &process',warning=.true.)
961 &
call writo(
'For VMEC, you should not go all the way to &
962 &the magnetic axis.',warning=.true.)
964 end function adapt_sol_grid
970 subroutine adapt_rich
971 if (max_it_rich.lt.1)
then
973 call writo(
'max_it_rich has been increased to 1',warning=.true.)
975 end subroutine adapt_rich
980 subroutine adapt_zero
981 if (max_it_zero.lt.1)
then
983 call writo(
'max_it_zero has been increased to 2',warning=.true.)
985 end subroutine adapt_zero
994 subroutine adapt_tols
998 real(dp),
parameter :: tol_slepc_def = 1.e-5_dp
1001 if (tol_norm.lt.0._dp)
then
1002 call writo(
'tol_norm has been increased to 0',warning=.true.)
1005 if (tol_norm.gt.1._dp)
then
1006 call writo(
'tol_norm has been decreased to 1',warning=.true.)
1011 if (tol_rich.lt.min_tol)
then
1012 call writo(
'tol_rich has been increased to '//&
1013 &trim(
r2str(min_tol)),warning=.true.)
1016 if (tol_rich.gt.max_tol)
then
1017 call writo(
'tol_rich has been decreased to '//&
1018 &trim(
r2str(max_tol)),warning=.true.)
1023 if (tol_zero.lt.min_tol)
then
1024 call writo(
'tol_zero has been increased to '//&
1025 &trim(
r2str(min_tol)),warning=.true.)
1028 if (tol_zero.gt.max_tol)
then
1029 call writo(
'tol_zero has been decreased to '//&
1030 &trim(
r2str(max_tol)),warning=.true.)
1035 allocate(tol_slepc_loc(max_it_rich))
1038 do id = 1,max_size_tol_slepc
1039 if (tol_slepc(id).lt.huge(1._dp))
then
1041 if (tol_slepc(id).lt.min_tol)
then
1043 &
') has been increased to '//trim(
r2str(min_tol)),&
1045 tol_slepc(id) = min_tol
1047 if (tol_slepc(id).gt.max_tol)
then
1049 &
') has been decreased to '//trim(
r2str(max_tol)),&
1051 tol_slepc(id) = max_tol
1056 if (max_id.gt.0)
then
1057 if (max_id.gt.max_it_rich)
then
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
1062 tol_slepc_loc = tol_slepc(1:max_it_rich)
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]',&
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
1072 tol_slepc_loc = tol_slepc_def
1074 end subroutine adapt_tols
1079 integer function adapt_normalization()
result(ierr)
1080 character(*),
parameter :: rout_name =
'adapt_normalization'
1083 character(len=max_str_ln) :: err_msg
1088 if (
rho_0.le.0)
then
1090 err_msg =
'rho_0 has to be positive'
1093 end function adapt_normalization
1109 character(*),
parameter :: rout_name =
'read_input_eq'
1194 use vmec_vars,
only: is_freeb_v, mnmax_v, mpol_v, ntor_v, is_asym_v, &
1203 character(*),
parameter :: rout_name =
'print_output_in'
1206 character(len=*),
intent(in) :: data_name
1207 logical,
intent(in),
optional :: remove_previous_arrs
1210 type(
var_1d_type),
allocatable,
target :: in_1d(:)
1213 integer :: in_limits(2)
1220 call writo(
'Write input variables to output file')
1225 call writo(
'Plotting decay of VMEC modes')
1227 &log10(maxval(abs(
r_v_c(:,:,0)),2)),draw=.false.)
1230 &log10(maxval(abs(
r_v_s(:,:,0)),2)),draw=.false.)
1233 &log10(maxval(abs(
z_v_c(:,:,0)),2)),draw=.false.)
1236 &log10(maxval(abs(
z_v_s(:,:,0)),2)),draw=.false.)
1239 &log10(maxval(abs(
l_v_c(:,:,0)),2)),draw=.false.)
1242 &log10(maxval(abs(
l_v_s(:,:,0)),2)),draw=.false.)
1245 &log10(maxval(abs(
jac_v_c(:,:,0)),2)),draw=.false.)
1248 &log10(maxval(abs(
jac_v_s(:,:,0)),2)),draw=.false.)
1252 &log10(maxval(abs(
b_v_c),2)),draw=.false.)
1255 &log10(maxval(abs(
b_v_s),2)),draw=.false.)
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')
1268 call writo(
'The entire normal range '//trim(
i2str(in_limits(1)))//&
1269 &
'..'//trim(
i2str(in_limits(2)))//
' is written')
1271 call writo(
'(relevant to solution range '//&
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))
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
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]
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),:),&
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]
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),:),&
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]
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),:),&
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]
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),:),&
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]
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),:),&
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)])
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),:)],&
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]
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
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),:)],&
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),:)],&
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])
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),:),&
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]
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]
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
1500 in_1d_loc%p = reshape([
r_h(:,in_limits(1):in_limits(2)),&
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))
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]
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),:),&
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]
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),:),&
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]
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),:),&
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]
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),:),&
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]
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),:),&
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]
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),:),&
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]
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
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))],&
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))
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))
1644 &remove_previous_arrs=remove_previous_arrs)
Print 2-D output on a file.
Variables that have to do with equilibrium quantities and the grid used in the calculations:
real(dp), public r_0
independent normalization constant for nondimensionalization
real(dp), public psi_0
derived normalization constant for nondimensionalization
real(dp), public max_flux_f
max. flux in Flux coordinates, set in calc_norm_range_PB3D_in
real(dp), public t_0
derived normalization constant for nondimensionalization
real(dp), public max_flux_e
max. flux in Equilibrium coordinates, set in calc_norm_range_PB3D_in
real(dp), public rho_0
independent normalization constant for nondimensionalization
real(dp), public pres_0
independent normalization constant for nondimensionalization
real(dp), public vac_perm
either usual mu_0 (default) or normalized
real(dp), public b_0
derived normalization constant for nondimensionalization
Operations that have to do with the grids and different coordinate systems.
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.
Variables pertaining to the different grids used.
real(dp), dimension(:), allocatable, public alpha
field line label alpha
real(dp), public min_par_x
min. of parallel coordinate [ ] in field-aligned grid
integer, public n_alpha
nr. of field-lines
integer, public n_r_eq
nr. of normal points in equilibrium grid
real(dp), public max_par_x
max. of parallel coordinate [ ] in field-aligned grid
integer, public n_r_in
nr. of normal points in input grid
real(dp), public max_alpha
max. of field-line label [ ] in field-aligned grid
real(dp), public min_alpha
min. of field-line label [ ] in field-aligned grid
integer, public n_r_sol
nr. of normal points in solution grid
Operations on HDF5 and XDMF variables.
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.
integer, parameter, public max_dim_var_1d
maximum dimension of var_1D
Operations on HELENA variables.
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
Numerical utilities related to giving output.
subroutine, public lvl_ud(inc)
Increases/decreases lvl of output.
subroutine, public writo(input_str, persistent, error, warning, alert)
Write output to file identified by output_i.
Numerical variables used by most other modules.
integer, public sol_n_procs
nr. of MPI processes for solution with SLEPC
integer, parameter, public dp
double precision
integer, public post_style
style for POST (1: extended grid, 2: B-aligned grid)
logical, public plot_b
whether to plot the magnetic field in real coordinates
integer, public n_theta_plot
nr. of poloidal points in plot
integer, public plot_grid_style
style for POST plot grid (0: 3-D plots, 1: slab plots, 2: slab plots with folding,...
integer, public norm_disc_prec_sol
precision for normal discretization for solution
logical, public plot_kappa
whether to plot curvature
character(len=max_str_ln), public input_name
will hold the full name of the input file
real(dp), public max_r_plot
max. of r_plot
real(dp), public min_theta_plot
min. of theta_plot
real(dp), public max_zeta_plot
max. of zeta_plot
logical, public compare_tor_pos
compare quantities at toroidal positions (only for POST)
real(dp), public max_njq_change
maximum change of prim. mode number times saf. fac. / rot. transf. when using X_style 2 (fast)
real(dp), public tol_norm
tolerance for normal range (normalized to 0..1)
integer, public n_procs
nr. of MPI processes
logical, public plot_flux_q
whether to plot flux quantities in real coordinates
integer, dimension(2), public plot_size
size of plot in inches
integer, public n_sol_requested
how many solutions requested
integer, parameter, public max_str_ln
maximum length of strings
logical, public plot_sol_q
whether to plot magnetic perturbation of solution in POST
real(dp), public ev_guess
first guess for eigenvalue
real(dp), public ev_bc
value of artificial Eigenvalue for boundary condition
logical, public post_output_sol
POST has outputs of solution.
integer, public prog_style
program style (1: PB3D, 2: PB3D_POST)
integer, dimension(2), public n_vac_plot
nr. of points in R and Z in vacuum
logical, public use_normalization
whether to use normalization or not
integer, public norm_disc_style_sol
style for normal discretization for solution (1: central fin. diff., 2: left fin. diff....
integer, public u_style
style for calculation of U (1: ord.2, 2: ord.1, 1: ord.0)
real(dp), public max_zvac_plot
max. of Z in which to plot vacuum
logical, public debug_version
debug version used
logical, public plot_vac_pot
whether to plot vacuum potential in POST
integer, public n_zeta_plot
nr. of toroidal points in plot
real(dp), public max_rvac_plot
max. of R in which to plot vacuum
logical, public plot_resonance
whether to plot the q-profile or iota-profile with resonances
integer, public x_grid_style
style for normal component of X grid (1: eq, 2: sol, 3: enriched)
real(dp), public min_zeta_plot
min. of zeta_plot
integer, public rho_style
style for equilibrium density profile
real(dp), parameter, public prog_version
version number
integer, parameter, public eq_i
file number of equilibrium file from VMEC or HELENA
logical, public post_output_full
POST has output on full grids.
integer, public eq_style
either 1 (VMEC) or 2 (HELENA)
integer, public max_it_slepc
maximum nr. of iterations for SLEPC
logical, public plot_vmec_modes
plot VMEC modes
character(len=max_str_ln), public pb3d_name
name of PB3D output file
integer, public norm_disc_prec_eq
precision for normal discretization for equilibrium
integer, dimension(4), public n_sol_plotted
how many solutions to be plot (first unstable, last unstable, first stable, last stable)
real(dp), public min_rvac_plot
min. of R in which to plot vacuum
real(dp), public tol_rich
tolerance for Richardson extrapolation
integer, public alpha_style
style for alpha (1: one field line, many turns, 2: many field lines, one turn)
real(dp), public pert_mult_factor_post
factor with which to multiply perturbation strength for POST
logical, public plot_j
whether to plot the current in real coordinates
integer, public max_nr_backtracks_hh
maximum number of backtracks for Householder, relax. factors
integer, dimension(2), public bc_style
style for BC left and right
integer, public ex_plot_style
external plot style (1: GNUPlot, 2: Bokeh for 2D, Mayavi for 3D)
integer, public ev_style
determines the method used for solving an EV problem
real(dp), public min_r_plot
min. of r_plot
integer, public rank
MPI rank.
real(dp), public min_zvac_plot
min. of Z in which to plot vacuum
integer, parameter, public input_i
file number of input file
integer, public rich_restart_lvl
starting Richardson level (0: none [default])
integer, public k_style
style for kinetic energy
logical, public plot_sol_xi
whether to plot plasma perturbation of solution in POST
integer, public norm_disc_prec_x
precision for normal discretization for perturbation
real(dp), dimension(:), allocatable, public tol_slepc
tolerance for SLEPC for different Richardson levels
integer, public norm_style
style for normalization
integer, public max_it_zero
maximum number of iterations to find zeros
logical, public plot_magn_grid
whether to plot the grid in real coordinates
integer, public x_style
style for secondary mode numbers (1: prescribed, 2: fast)
logical, public retain_all_sol
retain also faulty solutions
integer, public matrix_slepc_style
style for matrix storage (1: sparse, 2: shell)
real(dp), public max_theta_plot
max. of theta_plot
logical, public use_pol_flux_e
whether poloidal flux is used in E coords.
logical, public use_pol_flux_f
whether poloidal flux is used in F coords.
real(dp), public tol_zero
tolerance for zeros
real(dp), public max_tot_mem
maximum total memory for all processes [MB]
logical, public plot_e_rec
whether to plot energy reconstruction in POST
integer, public magn_int_style
style for magnetic integrals (1: trapezoidal, 2: Simpson 3/8)
integer, public max_it_rich
number of levels for Richardson extrapolation
integer, public solver_slepc_style
style for solver (1: Krylov-Schur, 2: GD)
Operations concerning giving output, on the screen as well as in output files.
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.
integer function, public find_max_rich_lvl(lvl_rich_req, max_lvl_rich_file)
Probe to find out which Richardson levels are available.
Variables concerning Richardson extrapolation.
integer, public min_n_par_x
min. of n_par_X (e.g. first value in Richardson loop)
integer, public rich_lvl
current level of Richardson extrapolation
integer, parameter, public req_min_n_par_x
required minimum n_par_X
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.
integer function, public read_vmec(n_r_in, use_pol_flux_v)
Reads the VMEC equilibrium data.
Variables that concern the output of VMEC.
real(dp), dimension(:,:), allocatable, public q_saf_v
safety factor
real(dp), dimension(:,:,:), allocatable, public jac_v_c
Coeff. of in sine series (HM and FM) and norm. deriv.
real(dp), dimension(:,:,:), allocatable, public b_v_sub_c
Coeff. of B_i in sine series (r,theta,phi) (FM).
integer, dimension(:,:), allocatable, public mn_v
m and n of modes
real(dp), dimension(:,:,:), allocatable, public l_v_s
Coeff. of in cosine series (HM) and norm. deriv.
real(dp), dimension(:,:,:), allocatable, public z_v_c
Coeff. of in sine series (FM) and norm. deriv.
real(dp), dimension(:,:), allocatable, public rot_t_v
rotational transform
real(dp), dimension(:,:,:), allocatable, public jac_v_s
Coeff. of in cosine series (HM and FM) and norm. deriv.
real(dp), dimension(:,:,:), allocatable, public r_v_c
Coeff. of in sine series (FM) and norm. deriv.
real(dp), dimension(:,:), allocatable, public pres_v
pressure
real(dp), dimension(:,:), allocatable, public j_v_sup_int
Integrated poloidal and toroidal current (FM).
real(dp), dimension(:,:,:), allocatable, public b_v_sub_s
Coeff. of B_i in cosine series (r,theta,phi) (FM).
real(dp), dimension(:,:), allocatable, public b_v_s
Coeff. of magnitude of B in cosine series (HM and FM).
real(dp), dimension(:,:), allocatable, public flux_t_v
toroidal flux
real(dp), dimension(:,:,:), allocatable, public z_v_s
Coeff. of in cosine series (FM) and norm. deriv.
real(dp), dimension(:,:,:), allocatable, public r_v_s
Coeff. of in cosine series (FM) and norm. deriv.
real(dp), dimension(:,:,:), allocatable, public l_v_c
Coeff. of in sine series (HM) and norm. deriv.
real(dp), dimension(:,:), allocatable, public b_v_c
Coeff. of magnitude of B in sine series (HM and FM).
real(dp), dimension(:,:), allocatable, public flux_p_v
poloidal flux
Variables pertaining to the perturbation quantities.
integer, public min_sec_x
m_X (pol. flux) or n_X (tor. flux) (only for X style 1)
real(dp), public max_r_sol
max. normal range for pert.
integer, public n_mod_x
size of m_X (pol. flux) or n_X (tor. flux)
real(dp), public min_r_sol
min. normal range for pert.
integer, public min_nm_x
minimum for the high-n theory (debable)
integer, public max_sec_x
m_X (pol. flux) or n_X (tor. flux) (only for\ c X style 1)
integer, public prim_x
n_X (pol. flux) or m_X (tor. flux)
1D equivalent of multidimensional variables, used for internal HDF5 storage.