5#include <PB3D_macros.h>
27 logical :: BR_normalization_provided(2)
50 module procedure calc_eq_1
52 module procedure calc_eq_2
92 module procedure print_output_eq_1
94 module procedure print_output_eq_2
105 module procedure redistribute_output_eq_1
107 module procedure redistribute_output_eq_2
128 module procedure calc_rzl_ind
130 module procedure calc_rzl_arr
141 module procedure calc_g_c_ind
143 module procedure calc_g_c_arr
158 module procedure calc_g_v_ind
160 module procedure calc_g_v_arr
171 module procedure calc_g_h_ind
173 module procedure calc_g_h_arr
185 module procedure calc_g_f_ind
187 module procedure calc_g_f_arr
198 module procedure calc_jac_v_ind
200 module procedure calc_jac_v_arr
215 module procedure calc_jac_h_ind
217 module procedure calc_jac_h_arr
234 module procedure calc_jac_f_ind
236 module procedure calc_jac_f_arr
247 module procedure calc_t_vc_ind
249 module procedure calc_t_vc_arr
260 module procedure calc_t_vf_ind
262 module procedure calc_t_vf_arr
273 module procedure calc_t_hf_ind
275 module procedure calc_t_hf_arr
280 integer function calc_eq_1(grid_eq,eq)
result(ierr)
287 character(*),
parameter :: rout_name =
'calc_eq_1'
290 type(
grid_type),
intent(inout) :: grid_eq
294 character(len=max_str_ln) :: err_msg
295 character(len=max_str_ln) :: file_name
296 character(len=max_str_ln) :: plot_title
302 call writo(
'Start setting up flux equilibrium quantities')
307 call eq%init(grid_eq)
314 call calc_flux_q_vmec()
316 call calc_flux_q_hel()
320 err_msg =
'Check out the FAQ at &
321 &https://pb3d.github.io/Doxygen/html/page_faq.html'
324 if (grid_eq%r_F(grid_eq%n(3)) .le. grid_eq%r_F(1))
then
326 plot_title =
'Flux normal coordinate'
327 call print_ex_2d(plot_title,file_name,grid_eq%r_F,draw=.false.)
328 call draw_ex([plot_title],file_name,1,1,.false.)
330 call writo(
'The code has not been tested satisfactorily for &
331 &decreasing normal coordinate. See plot '// trim(file_name),&
337 if (maxval(eq%q_saf_E(:,1))*minval(eq%q_saf_E(:,1)) .lt. 0._dp)
then
339 plot_title =
'derivative of safety factor'
340 call print_ex_2d(plot_title,file_name,eq%q_saf_E(:,1),draw=.false.)
341 call draw_ex([plot_title],file_name,1,1,.false.)
343 call writo(
'The code has to be adapted for reversed-shear regions &
344 &still. See plot ' // trim(file_name),alert=.true.)
349 grid_eq%loc_r_E = grid_eq%r_E(grid_eq%i_min:grid_eq%i_max)
350 grid_eq%loc_r_F = grid_eq%r_F(grid_eq%i_min:grid_eq%i_max)
363 err_msg =
'No density style associated with '//&
373 call writo(
'Done setting up flux equilibrium quantities')
380 subroutine calc_flux_q_vmec()
385 eq%flux_p_E =
flux_p_v(grid_eq%i_min:grid_eq%i_max,:)
386 eq%flux_t_E =
flux_t_v(grid_eq%i_min:grid_eq%i_max,:)
387 eq%pres_E =
pres_v(grid_eq%i_min:grid_eq%i_max,:)
388 eq%q_saf_E =
q_saf_v(grid_eq%i_min:grid_eq%i_max,:)
389 eq%rot_t_E =
rot_t_v(grid_eq%i_min:grid_eq%i_max,:)
405 end subroutine calc_flux_q_vmec
410 subroutine calc_flux_q_hel()
415 eq%flux_p_E =
flux_p_h(grid_eq%i_min:grid_eq%i_max,:)
416 eq%flux_t_E =
flux_t_h(grid_eq%i_min:grid_eq%i_max,:)
417 eq%pres_E =
pres_h(grid_eq%i_min:grid_eq%i_max,:)
418 eq%q_saf_E =
q_saf_h(grid_eq%i_min:grid_eq%i_max,:)
419 eq%rot_t_E =
rot_t_h(grid_eq%i_min:grid_eq%i_max,:)
431 end subroutine calc_flux_q_hel
432 end function calc_eq_1
434 integer function calc_eq_2(grid_eq,eq_1,eq_2,dealloc_vars)
result(ierr)
446 character(*),
parameter :: rout_name =
'calc_eq_2'
449 type(
grid_type),
intent(inout) :: grid_eq
450 type(eq_1_type),
intent(in) :: eq_1
451 type(eq_2_type),
intent(inout) :: eq_2
452 logical,
intent(in),
optional :: dealloc_vars
455 integer :: id, jd, kd
457 logical :: dealloc_vars_loc
458 character(len=max_str_ln) :: err_msg
464 call writo(
'Start setting up metric equilibrium quantities')
469 dealloc_vars_loc = .false.
470 if (
present(dealloc_vars)) dealloc_vars_loc = dealloc_vars
473 call eq_2%init(grid_eq)
482 call writo(
'Calculate R, Z, λ, ...')
483 if (.not.
allocated(grid_eq%trigon_factors))
then
485 &grid_eq%trigon_factors)
502 call writo(
'Calculate g_C')
509 call writo(
'Calculate T_VC')
516 call writo(
'Calculate g_V')
523 call writo(
'Calculate jac_V')
531 call writo(
'Test calculation of g_V?')
533 ierr = test_g_v(grid_eq,eq_2)
537 call writo(
'Test calculation of jac_V?')
539 ierr = test_jac_v(grid_eq,eq_2)
547 call writo(
'Calculate T_VF')
557 if (dealloc_vars_loc)
then
568 do kd = 1,grid_eq%loc_n_r
569 do jd = 1,grid_eq%n(2)
570 do id = 1,grid_eq%n(1)
571 if (abs(grid_eq%theta_E(id,jd,kd)-
chi_h(id)).gt.&
574 err_msg =
'theta_E is not identical to chi_H'
577 if (abs(grid_eq%zeta_E(id,jd,kd)-0._dp).gt.&
580 err_msg =
'zeta_E is not identical to 0'
588 call writo(
'Test consistency of metric factors?')
589 if(
get_log(.false.,ind=.true.))
then
595 call writo(
'Test harmonic content of HELENA output?')
596 if(
get_log(.false.,ind=.true.))
then
605 call writo(
'Calculate jac_H')
612 call writo(
'Calculate g_H')
619 call writo(
'Calculate h_H')
627 call writo(
'Test calculation of D1 D2 h_H?')
629 ierr = test_d12h_h(grid_eq,eq_2)
638 call writo(
'Exporting HELENA equilibrium for VMEC porting')
640 ierr = create_vmec_input(grid_eq,eq_1)
643 call writo(
'Done exporting')
647 call writo(
'Calculate T_HF')
657 if (dealloc_vars_loc)
then
664 call writo(
'Test calculation of T_EF?')
666 ierr = test_t_ef(grid_eq,eq_1,eq_2)
674 call writo(
'Calculate T_FE')
683 call writo(
'Calculate g_F')
690 call writo(
'Calculate h_F')
697 call writo(
'Calculate jac_F')
704 eq_2%jac_F(:,:,:,0,0,0) = sign(&
705 &max(
tol_zero,abs(eq_2%jac_F(:,:,:,0,0,0))),eq_2%jac_F(:,:,:,0,0,0))
706 eq_2%jac_E(:,:,:,0,0,0) = sign(&
707 &max(
tol_zero,abs(eq_2%jac_E(:,:,:,0,0,0))),eq_2%jac_E(:,:,:,0,0,0))
710 if (dealloc_vars_loc)
then
719 if (dealloc_vars_loc)
then
720 deallocate(eq_2%g_F,eq_2%h_F,eq_2%jac_F)
721 deallocate(eq_2%det_T_EF,eq_2%det_T_FE)
730 call writo(
'Test Jacobian in Flux coordinates?')
732 ierr = test_jac_f(grid_eq,eq_1,eq_2)
736 call writo(
'Test calculation of B_F?')
738 ierr = test_b_f(grid_eq,eq_1,eq_2)
742 call writo(
'Test consistency with given pressure?')
744 ierr = test_p(grid_eq,eq_1,eq_2)
753 call writo(
'Done setting up metric equilibrium quantities')
754 end function calc_eq_2
783 integer function create_vmec_input(grid_eq,eq_1)
result(ierr)
804 character(*),
parameter :: rout_name =
'create_VMEC_input'
811 type(ezspline2_r8) :: f_spl
812 integer :: id, jd, kd, ld
815 integer :: tot_nr_pert
818 integer :: n_loc, m_loc
819 integer :: pert_map_n_loc
821 integer :: n_pert_map(2)
822 integer :: plot_dim(2)
824 integer :: pert_style
826 integer :: max_n_b_output
827 integer :: n_prop_b_tor
831 integer :: m_range(2)
832 integer,
allocatable :: n_pert(:)
833 integer,
allocatable :: n_pert_copy(:)
834 integer,
allocatable :: piv(:)
835 character(len=1) :: pm(3)
836 character(len=8) :: flux_name(2)
837 character(len=6) :: path_prefix =
'../../'
838 character(len=max_str_ln) :: hel_pert_file_name
839 character(len=max_str_ln) :: err_msg
840 character(len=max_str_ln) :: file_name
841 character(len=max_str_ln) :: plot_name(2)
842 character(len=max_str_ln) :: plot_title(2)
843 character(len=max_str_ln) :: prop_b_tor_file_name
844 real(
dp) :: eq_vert_shift
845 real(
dp) :: pert_map_vert_shift
846 real(
dp) :: max_pert_on_axis
847 real(
dp) :: m_tol = 1.e-7_dp
848 real(
dp) :: delta_loc(2)
849 real(
dp) :: plot_lims(2,2)
852 real(
dp) :: prop_b_tor_smooth
853 real(
dp) :: rz_b_0(2)
855 real(
dp) :: pres_vj(99,0:1)
856 real(
dp) :: rot_t_v(99)
858 real(
dp) :: ffp_vj(99)
859 real(
dp) :: flux_j(99)
860 real(
dp) :: q_saf_vj(99)
861 real(
dp) :: i_tor_v(99)
862 real(
dp) :: i_tor_j(99)
863 real(
dp) :: norm_j(3)
864 real(
dp),
allocatable :: psi_t(:)
865 real(
dp),
allocatable :: ffp(:)
866 real(
dp),
allocatable :: i_tor(:)
867 real(
dp),
allocatable :: norm_transf(:,:)
868 real(
dp),
allocatable :: i_tor_int(:)
869 real(
dp),
allocatable :: rrint(:)
870 real(
dp),
allocatable :: rrint_loc(:)
871 real(
dp),
allocatable :: r_plot(:,:,:)
872 real(
dp),
allocatable :: z_plot(:,:,:)
873 real(
dp),
allocatable :: pert_map_r(:)
874 real(
dp),
allocatable :: pert_map_z(:)
875 real(
dp),
allocatable :: pert_map(:,:)
876 real(
dp),
allocatable :: pert_map_interp(:)
877 real(
dp),
allocatable :: pert_map_interp_f(:,:)
878 real(
dp),
allocatable :: r_h_loc(:,:)
879 real(
dp),
allocatable :: z_h_loc(:,:)
880 real(
dp),
allocatable :: delta(:,:,:)
881 real(
dp),
allocatable :: delta_copy(:,:,:)
882 real(
dp),
allocatable :: bh_0(:,:)
883 real(
dp),
allocatable :: b_f(:,:,:,:)
884 real(
dp),
allocatable :: b_f_copy(:,:,:,:)
885 real(
dp),
allocatable :: b_f_dum(:,:)
886 real(
dp),
allocatable :: b_f_dum2(:,:)
887 real(
dp),
allocatable :: b_f_dum3(:,:)
888 real(
dp),
allocatable :: theta_b(:)
889 real(
dp),
allocatable :: theta_geo(:,:,:)
890 real(
dp),
allocatable :: prop_b_tor(:,:)
891 real(
dp),
allocatable :: prop_b_tor_ord(:,:)
892 real(
dp),
allocatable :: prop_b_tor_interp(:)
893 real(
dp),
allocatable :: prop_b_tor_interp_f(:,:)
894 real(
dp),
allocatable :: xyz_plot(:,:,:,:)
895 logical :: zero_n_pert
898 logical :: change_max_n_b_output
901 real(
dp),
allocatable :: bh_0_alt(:,:)
910 err_msg =
'This routine should be run with the full normal &
916 if (.not.all(br_normalization_provided,1)) &
917 &
call writo(
'No normalization factors were provided. &
918 &Are you sure this is correct?',warning=.true.)
932 call writo(
'Was the equilibrium shifted vertically?')
933 eq_vert_shift = 0._dp
935 call writo(
'How much higher [m] should the &
936 &equilibrium be in the real world?')
938 call writo(
'The equilibrium will be shifted by '//&
939 &trim(
r2strt(eq_vert_shift))//
'm to compensate for this')
940 z_h_loc = z_h_loc + eq_vert_shift
946 allocate(rrint_loc(
nchi))
953 rrint(kd) = rrint_loc(
nchi)
955 if (
ias.eq.0) rrint = 2*rrint
973 allocate(norm_transf(
n_r_eq,2))
976 i_tor = - (eq_1%pres_E(:,1)*rrint*eq_1%q_saf_E(:,0)/
rbphi_h(:,0) + &
980 allocate(i_tor_int(
size(i_tor)))
981 ierr =
calc_int(i_tor*norm_transf(:,1),&
987 call writo(
'Initialize boundary')
996 allocate(bh_0(n_b,2))
997 allocate(theta_b(n_b))
1006 bh_0(:,1) = r_h_loc(1:n_b,
n_r_eq)
1007 bh_0(:,2) = z_h_loc(1:n_b,
n_r_eq)
1009 rz_b_0(1) = sum(r_h_loc(:,1))/
size(r_h_loc,1)
1010 rz_b_0(2) = sum(z_h_loc(:,1))/
size(z_h_loc,1)
1011 theta_b = atan2(bh_0(:,2)-rz_b_0(2),bh_0(:,1)-rz_b_0(1))
1012 where (theta_b.lt.0) theta_b = theta_b + 2*
pi
1013 call writo(
'Magnetic axis used for geometrical coordinates:')
1015 call writo(
'('//trim(
r2str(rz_b_0(1)))//
'm, '//&
1016 &trim(
r2str(rz_b_0(2)))//
'm)')
1021 call writo(
'Test with circular tokamak?')
1024 call writo(
'Testing with circular tokamak:')
1026 call writo(
'R = 3/2 + 1/2 cos(θ)')
1027 call writo(
'Z = 1/2 sin(θ)')
1028 call writo(
'with θ equidistant')
1030 theta_b = [((id-1._dp)/n_b*2*
pi,id=1,n_b)]
1031 bh_0(:,1) = 1.5_dp + 0.5*cos(theta_b)
1032 bh_0(:,2) = 0.5*sin(theta_b)
1039 allocate(theta_geo(grid_eq%n(1),1,grid_eq%loc_n_r))
1040 do kd = 1,grid_eq%loc_n_r
1041 theta_geo(:,1,kd) = atan2(z_h_loc(:,kd)-rz_b_0(2),&
1042 &r_h_loc(:,kd)-rz_b_0(1))
1044 where (theta_geo.lt.0._dp) theta_geo = theta_geo + 2*
pi
1045 allocate(xyz_plot(grid_eq%n(1),1,grid_eq%loc_n_r,3))
1047 &xyz_plot(:,:,:,2),xyz_plot(:,:,:,3))
1049 call plot_hdf5(
'theta_geo',
'theta_geo',theta_geo,&
1050 &tot_dim=[grid_eq%n(1),1,grid_eq%n(3),3],&
1051 &loc_offset=[0,0,grid_eq%i_min-1,0],&
1052 &x=xyz_plot(:,:,:,1),y=xyz_plot(:,:,:,2),z=xyz_plot(:,:,:,3),&
1053 &descr=
'geometric poloidal angle used to create the VMEC &
1055 deallocate(xyz_plot)
1059 plot_title = [
'R_H',
'Z_H']
1060 call print_ex_2d(plot_title(1:2),plot_name(1),bh_0,&
1061 &x=reshape(theta_b,[n_b,1])/
pi,&
1063 call draw_ex(plot_title(1:2),plot_name(1),2,1,.false.)
1066 flux_name = [
'poloidal',
'toroidal']
1072 plot_lims(:,1) = [0.0_dp,2.0_dp]
1073 plot_lims(:,2) = [0.5_dp,2.0_dp]
1074 call writo(
'Change plot properties from defaults?')
1077 call writo(trim(
i2str(plot_dim(id)))//
' geometrical '//&
1078 &flux_name(id)//
' points on range '//&
1079 &trim(
r2strt(plot_lims(1,id)))//
' pi .. '//&
1080 &trim(
r2strt(plot_lims(2,id)))//
' pi')
1086 call writo(
'number of '//flux_name(id)//
' points?')
1087 plot_dim(id) =
get_int(lim_lo=2)
1088 call writo(
'lower limit of '//flux_name(id)//&
1091 call writo(
'upper limit of '//flux_name(id)//&
1099 call writo(
'VMEC current style?')
1101 call writo(
'0: prescribe iota')
1102 call writo(
'1: prescribe toroidal current')
1103 ncurr =
get_int(lim_lo=0,lim_hi=1)
1109 call writo(
'Do you want to perturb the equilibrium?')
1112 zero_n_pert = .false.
1115 call writo(
'Which perturbation do you want to describe?')
1117 call writo(
'1: plasma boundary position')
1118 call writo(
'2: B_tor magnetic ripple with fixed N')
1120 pert_style =
get_int(lim_lo=1,lim_hi=2)
1123 if (pert_style.eq.2)
then
1126 do while (.not.found)
1127 call writo(
'Proportionality file name?')
1129 read(*,*,iostat=ierr) prop_b_tor_file_name
1130 err_msg =
'failed to read prop_B_tor_file_name'
1135 call writo(
'Trying "'//path_prefix(1:kd*3)//&
1136 &trim(prop_b_tor_file_name)//
'"')
1138 &trim(prop_b_tor_file_name),iostat=ierr,&
1142 call writo(
"Success")
1150 call writo(
'Could not open file "'//&
1151 &trim(prop_b_tor_file_name)//
'"')
1158 call writo(
'Analyzing perturbation proportionality file "'&
1159 &//trim(prop_b_tor_file_name)//
'"')
1162 allocate(prop_b_tor(n_prop_b_tor,2))
1164 &file_name=prop_b_tor_file_name)
1166 do id = 1,n_prop_b_tor
1172 &
' poloidal points '//&
1173 &trim(
r2strt(minval(prop_b_tor(:,1))))//
'..'//&
1174 &trim(
r2strt(maxval(prop_b_tor(:,1)))))
1175 call writo(
'The proportionality factor should be tabulated &
1176 &in a geometrical angle that has the SAME origin as &
1177 &the one used here',alert=.true.)
1183 &tol=0.5_dp*
pi/n_prop_b_tor)
1185 deallocate(prop_b_tor)
1186 allocate(prop_b_tor_interp(n_b))
1187 ierr =
spline(prop_b_tor_ord(:,1),prop_b_tor_ord(:,2),&
1191 deallocate(prop_b_tor_ord)
1192 if (grid_eq%n(2).ne.1)
call writo(
'There should be only 1 &
1193 &geodesical position, but there are '//&
1194 &trim(
i2str(grid_eq%n(2))),warning=.true.)
1197 prop_b_tor_smooth = 1._dp
1198 call writo(
'Do you want to smooth the perturbation?')
1199 if (
get_log(.false.)) prop_b_tor_smooth = &
1200 &
get_real(lim_lo=0._dp,lim_hi=1._dp)
1203 ierr =
nufft(theta_b,prop_b_tor_interp,prop_b_tor_interp_f)
1205 deallocate(prop_b_tor_interp)
1211 call writo(
'How do you want to prescribe the perturbation?')
1213 call writo(
'1: through a file with Fourier modes in &
1214 &geometrical coordinates')
1216 call writo(
'It should provide N M delta_cos delta_sin in four &
1218 call writo(
'Rows starting with "#" are ignored')
1220 call writo(
'2: same as 1 but interactively')
1221 call writo(
'3: through a 2-D map in R, Z for constant &
1222 &geometrical coordinates')
1224 call writo(
'It should provide R, Z and delta')
1225 call writo(
'Rows starting with "#" are ignored')
1227 pert_type =
get_int(lim_lo=1,lim_hi=3)
1230 select case (pert_type)
1233 do while (.not.found)
1235 call writo(
'File name?')
1237 read(*,*,iostat=ierr) hel_pert_file_name
1238 err_msg =
'failed to read HEL_pert_file_name'
1243 call writo(
'Trying "'//path_prefix(1:kd*3)//&
1244 &trim(hel_pert_file_name)//
'"')
1246 &trim(hel_pert_file_name),iostat=ierr,&
1250 call writo(
"Success")
1258 call writo(
'Could not open file "'//&
1259 &trim(hel_pert_file_name)//
'"')
1265 call writo(
'Parsing "'//trim(hel_pert_file_name)//
'"')
1268 if (pert_type.eq.1)
then
1271 else if (pert_type.eq.3)
then
1273 call writo(
'Was the equilibrium shifted &
1275 if (abs(eq_vert_shift).gt.0._dp)
then
1277 call writo(
'Note: If you already shifted just &
1278 &now, don''t do it again!',alert=.true.)
1281 pert_map_vert_shift = 0._dp
1283 call writo(
'How much higher [m] should the &
1284 &equilibrium be in the real world?')
1286 call writo(
'The perturbation map will be &
1288 &trim(
r2strt(-pert_map_vert_shift))//&
1289 &
'm to compensate for this')
1294 &file_name=hel_pert_file_name)
1297 allocate(pert_map_r(n_pert_map(1)))
1298 do kd = 1,n_pert_map(1)
1302 &file_name=hel_pert_file_name)
1305 allocate(pert_map_z(n_pert_map(2)))
1306 do kd = 1,n_pert_map(2)
1309 pert_map_z = pert_map_z - pert_map_vert_shift
1311 &file_name=hel_pert_file_name)
1313 allocate(pert_map(n_pert_map(1),n_pert_map(2)))
1314 do kd = 1,n_pert_map(1)
1319 allocate(r_plot(n_pert_map(1),n_pert_map(2),1))
1320 allocate(z_plot(n_pert_map(1),n_pert_map(2),1))
1321 do kd = 1,n_pert_map(2)
1322 r_plot(:,kd,1) = pert_map_r
1324 do kd = 1,n_pert_map(1)
1325 z_plot(kd,:,1) = pert_map_z
1328 &reshape(pert_map,[n_pert_map,1]),&
1329 &x=r_plot,y=z_plot,descr=&
1330 &
'perturbation map for '//&
1331 &trim(hel_pert_file_name))
1332 deallocate(r_plot,z_plot)
1335 if (minval(r_h_loc).lt.minval(pert_map_r) .or. &
1336 &maxval(r_h_loc).gt.maxval(pert_map_r) .or. &
1337 &minval(z_h_loc).lt.minval(pert_map_z) .or. &
1338 &maxval(z_h_loc).gt.maxval(pert_map_z))
then
1340 call writo(
'Are you sure you have specified &
1341 &a normalization factor R_0?',alert=.true.)
1342 err_msg =
'R and Z are not contained in &
1348 allocate(pert_map_interp(n_b))
1351 call ezspline_init(f_spl,n_pert_map(1),n_pert_map(2),&
1352 &bcs(:,1),bcs(:,2),ierr)
1353 call ezspline_error(ierr)
1357 f_spl%x1 = pert_map_r
1358 f_spl%x2 = pert_map_z
1361 call ezspline_setup(f_spl,pert_map,ierr,&
1363 call ezspline_error(ierr)
1368 call ezspline_interp(f_spl,bh_0(id,1),bh_0(id,2),&
1369 &pert_map_interp(id),ierr)
1370 call ezspline_error(ierr)
1374 &pert_map_interp,x=theta_b,draw=.false.)
1375 call draw_ex([
'ripple'],
'pert_map_interp',&
1379 ierr =
nufft(theta_b,pert_map_interp,&
1382 tot_nr_pert =
size(pert_map_interp_f,1)*2
1385 call writo(
'toiroidal period for 2-D map?')
1386 pert_map_n_loc =
get_int(lim_lo=0)
1392 call writo(
'Prescribe interactively')
1393 call writo(
'How many different combinations of &
1394 &toroidal and poloidal mode numbers (N,M)?')
1395 tot_nr_pert =
get_int(lim_lo=1)
1399 allocate(n_pert(tot_nr_pert+1))
1400 allocate(delta(tot_nr_pert+1,0:0,2))
1403 do jd = 1,tot_nr_pert
1404 select case (pert_type)
1407 &file_name=hel_pert_file_name)
1409 read(
hel_pert_i,*,iostat=ierr) n_loc, m_loc, &
1411 err_msg =
'Could not read file '//&
1412 &trim(hel_pert_file_name)
1415 call writo(
'For mode '//trim(
i2str(jd))//
'/'//&
1416 &trim(
i2str(tot_nr_pert))//
':')
1420 call writo(
'Toroidal mode number N?')
1422 call writo(
'Poloidal mode number M?')
1424 call writo(
'Perturbation strength ~ cos('//&
1425 &trim(
i2str(m_loc))//
' θ - '//&
1426 &trim(
i2str(n_loc))//
' ζ)?')
1428 call writo(
'Perturbation strength ~ sin('//&
1429 &trim(
i2str(m_loc))//
' θ - '//&
1430 &trim(
i2str(n_loc))//
' ζ)?')
1436 delta_loc = pert_map_interp_f(m_loc+1,:)*0.5_dp
1437 n_loc = pert_map_n_loc
1438 if (mod(jd,2).eq.0) n_loc = -n_loc
1442 if (n_loc.eq.0) zero_n_pert = .true.
1447 if (n_loc.eq.n_pert(id)) n_id = id
1449 if (n_id.gt.nr_n)
then
1451 n_pert(n_id) = n_loc
1453 m_range = [lbound(delta,2),ubound(delta,2)]
1454 if (m_loc.gt.m_range(1) .or. m_loc.lt.m_range(2))
then
1455 allocate(delta_copy(tot_nr_pert+1,&
1456 &m_range(1):m_range(2),2))
1459 allocate(delta(tot_nr_pert+1,&
1460 &min(m_range(1),m_loc):max(m_range(2),m_loc),2))
1462 delta(:,m_range(1):m_range(2),:) = delta_copy
1463 deallocate(delta_copy)
1465 delta(n_id,m_loc,:) = delta(n_id,m_loc,:) + delta_loc
1469 call writo(
'perturbation '//trim(
i2str(jd))//&
1470 &
' at position '//trim(
i2str(n_id))//
', adding ('//&
1471 &trim(
r2strt(delta(n_id,m_loc,1)))//
', '//&
1472 &trim(
r2strt(delta(n_id,m_loc,2)))//
') at (n,m) = ('//&
1473 &trim(
i2str(n_loc))//
', '//trim(
i2str(m_loc))//
')')
1475 do kd = lbound(delta,2),ubound(delta,2)
1476 call writo(
'delta ('//trim(
i2str(kd))//
') = ('//&
1477 &trim(
r2strt(delta(n_id,kd,1)))//
', '//&
1478 &trim(
r2strt(delta(n_id,kd,2)))//
')')
1485 if (pert_type.eq.1 .or. pert_type.eq.3)
close(
hel_pert_i)
1490 max_pert_on_axis = sum(sqrt(sum(delta**2,3)))
1491 select case (pert_style)
1493 call writo(
'Maximum absolute perturbation on axis: '//&
1494 &trim(
r2strt(100*max_pert_on_axis))//
'cm')
1495 call writo(
'Rescale to some value [cm]?')
1497 call writo(
'Maximum relative perturbation on axis: '//&
1498 &trim(
r2strt(100*max_pert_on_axis))//
'%')
1499 call writo(
'Rescale to some value [%]?')
1501 mult_fac = max_pert_on_axis
1504 mult_fac = mult_fac * 0.01_dp
1506 delta = delta * mult_fac / max_pert_on_axis
1510 if (.not.zero_n_pert)
then
1519 call writo(
'Set up axisymmetric data')
1523 nr_m_max = (n_b-1)/2
1525 nr_m_max = max(nr_m_max, max(-lbound(delta,2),ubound(delta,2)))
1526 if (pert_style.eq.2) &
1527 &nr_m_max = max(nr_m_max,
size(prop_b_tor_interp_f,1)-1)
1529 plot_name(1) =
'R_F'
1530 plot_name(2) =
'Z_F'
1531 allocate(b_f(nr_n,-nr_m_max:nr_m_max,2,2))
1536 call writo(
'analyzing '//trim(plot_name(kd)))
1540 ierr =
nufft(theta_b,bh_0(:,kd),b_f_dum,plot_name(kd))
1542 b_f(1,0:
size(b_f_dum,1)-1,:,kd) = b_f_dum
1546 allocate(bh_0_alt(
size(bh_0,1),
size(b_f_dum,1)))
1548 call writo(
'Comparing R or Z with reconstruction &
1549 &through Fourier coefficients')
1550 bh_0_alt(:,1) = b_f_dum(1,1)
1551 do id = 1,
size(b_f_dum,1)-1
1552 bh_0_alt(:,id+1) = bh_0_alt(:,id) + &
1553 &b_f_dum(id+1,1)*cos(id*theta_b) + &
1554 &b_f_dum(id+1,2)*sin(id*theta_b)
1557 &reshape([bh_0(:,kd),bh_0_alt(:,
size(b_f_dum,1))],&
1558 &[
size(bh_0,1),2]),x=&
1559 &reshape([theta_b],[
size(bh_0,1),1]))
1561 call writo(
'Plotting Fourier approximation')
1562 call print_ex_2d([
'alt BH'],
'TEST_'//trim(plot_name(kd))//&
1563 &
'_F_series',bh_0_alt,&
1564 &x=reshape([theta_b],[
size(bh_0,1),1]))
1565 call draw_ex([
'alt BH'],
'TEST_'//trim(plot_name(kd))//&
1566 &
'_F_series',
size(b_f_dum,1),1,.false.)
1568 call writo(
'Making animation')
1569 call print_ex_2d([
'alt BH'],
'TEST_'//trim(plot_name(kd))//&
1570 &
'_F_series_anim',bh_0_alt,&
1571 &x=reshape([theta_b],[
size(bh_0,1),1]))
1572 call draw_ex([
'alt BH'],
'TEST_'//trim(plot_name(kd))//&
1573 &
'_F_series_anim',
size(b_f_dum,1),1,.false.,&
1574 &is_animated=.true.)
1576 deallocate(bh_0_alt)
1586 call plot_boundary(b_f(1:1,:,:,:),[-nr_m_max,nr_m_max],[0],&
1587 &
'last_fs',plot_dim,plot_lims)
1594 call writo(
'Set up perturbed data')
1598 allocate(n_pert_copy(nr_n))
1599 allocate(delta_copy(nr_n,-nr_m_max:nr_m_max,2))
1601 n_pert_copy = n_pert(1:nr_n)
1602 delta_copy(:,lbound(delta,2):ubound(delta,2),:) = delta(1:nr_n,:,:)
1603 deallocate(n_pert);
allocate(n_pert(nr_n)); n_pert = n_pert_copy
1604 deallocate(delta);
allocate(delta(nr_n,-nr_m_max:nr_m_max,2))
1605 deallocate(n_pert_copy)
1609 delta(jd,:,:) = delta_copy(piv(jd),:,:)
1614 id_n_0 = minloc(abs(n_pert),1)
1615 if (id_n_0.gt.1)
then
1616 b_f(id_n_0,:,:,:) = b_f(1,:,:,:)
1617 b_f(1,:,:,:) = 0._dp
1621 call writo(
'Summary of perturbation form:')
1624 do id = -nr_m_max,nr_m_max
1625 if (maxval(abs(delta(jd,id,:))).lt.
tol_zero) cycle
1626 if (n_pert(jd).ge.0)
then
1631 if (delta(jd,id,1).ge.0)
then
1636 if (jd.eq.1 .and. id.eq.1) pm(2) =
''
1637 if (delta(jd,id,2).ge.0)
then
1643 &pm(2)//
' '//trim(
r2strt(abs(delta(jd,id,1))))//&
1644 &
' cos('//trim(
i2str(id))//
' θ '//pm(1)//
' '//&
1645 &trim(
i2str(abs(n_pert(jd))))//
' ζ) '//&
1646 &pm(3)//
' '//trim(
r2strt(abs(delta(jd,id,2))))//&
1647 &
' sin('//trim(
i2str(id))//
' θ '//pm(1)//
' '//&
1648 &trim(
i2str(abs(n_pert(jd))))//
' ζ)')
1685 m_range = [-nr_m_max,nr_m_max]
1686 pert_n:
do jd = 1,nr_n
1688 allocate(b_f_dum(-nr_m_max:nr_m_max,2))
1689 allocate(b_f_dum2(-nr_m_max:nr_m_max,2))
1690 allocate(b_f_dum3(-nr_m_max:nr_m_max,2))
1695 b_f_dum(1,kd) = 1._dp
1696 if (pert_style.eq.2)
then
1697 b_f_dum2(0:ubound(prop_b_tor_interp_f,1),:) = &
1698 &prop_b_tor_interp_f
1699 call shift_f(m_range,m_range,m_range,&
1700 &b_f_dum,b_f_dum2,b_f_dum3)
1705 b_f_dum2(:,1) = delta(jd,:,1)
1706 b_f_dum2(:,2) = delta(jd,:,2)
1707 call shift_f(m_range,m_range,m_range,&
1708 &b_f_dum,b_f_dum2,b_f_dum3)
1709 b_f_dum3 = b_f_dum3*0.5_dp
1711 b_f(jd,:,1,kd) = b_f(jd,:,1,kd) + b_f_dum3(:,1)
1713 b_f(jd,:,2,kd) = b_f(jd,:,2,kd) + b_f_dum3(:,2)
1715 b_f(jd,:,1,kd) = b_f(jd,:,1,kd) + &
1716 &b_f_dum3(nr_m_max:-nr_m_max:-1,1)
1718 b_f(jd,:,2,kd) = b_f(jd,:,2,kd) - &
1719 &b_f_dum3(nr_m_max:-nr_m_max:-1,2)
1722 b_f_dum2(:,1) = -delta(jd,:,2)
1723 b_f_dum2(:,2) = delta(jd,:,1)
1724 call shift_f(m_range,m_range,m_range,&
1725 &b_f_dum,b_f_dum2,b_f_dum3)
1726 b_f_dum3 = b_f_dum3*0.5_dp
1728 b_f(jd,:,1,kd) = b_f(jd,:,1,kd) + b_f_dum3(:,2)
1730 b_f(jd,:,2,kd) = b_f(jd,:,2,kd) - b_f_dum3(:,1)
1732 b_f(jd,:,1,kd) = b_f(jd,:,1,kd) - &
1733 &b_f_dum3(nr_m_max:-nr_m_max:-1,2)
1735 b_f(jd,:,2,kd) = b_f(jd,:,2,kd) - &
1736 &b_f_dum3(nr_m_max:-nr_m_max:-1,1)
1740 deallocate(b_f_dum2)
1741 deallocate(b_f_dum3)
1744 allocate(b_f_copy(2*nr_n,0:nr_m_max,2,2))
1745 allocate(n_pert_copy(2*nr_n))
1749 allocate(b_f_dum(0:nr_m_max,2))
1750 pert_n_convert:
do jd = 1,
size(b_f,1)
1759 b_f_dum(0:nr_m_max,1) = b_f(jd,0:-nr_m_max:-1,1,ld)
1760 b_f_dum(0:nr_m_max,2) = -b_f(jd,0:-nr_m_max:-1,2,ld)
1764 b_f_dum(0:nr_m_max,:) = b_f(jd,0:nr_m_max,:,ld)
1766 b_f_dum(0,:) = b_f_dum(0,:)*0.5_dp
1768 if (maxval(abs(b_f_dum)).gt.
tol_zero)
then
1772 if (n_loc.eq.n_pert_copy(id)) n_id = id
1774 if (n_id.gt.nr_n)
then
1776 n_pert_copy(n_id) = n_loc
1780 call writo(
'putting n = '//trim(
i2str(n_loc))//&
1781 &
' in index '//trim(
i2str(n_id))//
':')
1783 call writo(
'(negative range)')
1785 call writo(
'(positive range)')
1788 call writo(
'(for R)')
1790 call writo(
'(for Z)')
1794 if (maxval(abs(b_f_dum(id,:)))&
1798 &trim(
r2strt(b_f_dum(id,1)))//&
1800 &trim(
r2strt(b_f_dum(id,2)))//
')')
1808 b_f_copy(n_id,:,:,ld) = b_f_copy(n_id,:,:,ld) + &
1813 end do pert_n_convert
1818 allocate(b_f(nr_n,0:nr_m_max,2,2))
1819 allocate(n_pert(nr_n))
1820 b_f = b_f_copy(1:nr_n,:,:,:)
1821 n_pert = n_pert_copy(1:nr_n)
1822 deallocate(b_f_copy)
1823 deallocate(n_pert_copy)
1827 call plot_boundary(b_f,[0,nr_m_max],n_pert(1:nr_n),
'last_fs_pert',&
1828 &plot_dim,plot_lims)
1837 nfp =
gcd(nfp,n_pert(jd))
1847 allocate(b_f_copy(1,0:nr_m_max,2,2))
1849 b_f_copy(1,id,1,:) = b_f(1,id,1,:) + b_f(1,-id,1,:)
1850 b_f_copy(1,id,2,:) = b_f(1,id,2,:) - b_f(1,-id,2,:)
1852 b_f_copy(1,0,:,:) = b_f_copy(1,0,:,:)*0.5_dp
1854 allocate(b_f(1,0:nr_m_max,2,2))
1856 deallocate(b_f_copy)
1868 file_name =
"input."//trim(
eq_name)
1870 select case (pert_type)
1872 file_name = trim(file_name)//
'_'//&
1873 &trim(hel_pert_file_name)
1876 file_name = trim(file_name)//
'_N'//&
1877 &trim(
i2str(n_pert(jd)))
1879 file_name = trim(file_name)//
'M'//&
1884 file_name = trim(file_name)//
'_2DMAP'
1887 call writo(
'Generate VMEC input file "'//trim(file_name)//
'"')
1888 call writo(
'This can be used for VMEC porting')
1893 if (maxval(abs(b_f(:,:,2,1)))/maxval(abs(b_f(:,:,1,1))) &
1895 maxval(abs(b_f(:,:,1,2)))/maxval(abs(b_f(:,:,2,2))) &
1899 call writo(
"The equilibrium configuration has stellarator &
1902 call writo(
"The equilibrium configuration does not have &
1903 &stellarator symmetry")
1908 norm_b_h = maxval(abs(b_f))
1910 if (maxval(abs(b_f(:,id,:,1)/norm_b_h)).gt.m_tol .or. &
1911 &maxval(abs(b_f(:,id,:,2)/norm_b_h)).gt.m_tol) &
1914 call writo(
"Detected recommended number of poloidal modes: "//&
1915 &trim(
i2str(rec_min_m)))
1918 max_n_b_output = rec_min_m
1919 call writo(
'Do you want to change the maximum number of modes &
1920 &from current '//trim(
i2str(max_n_b_output))//
'?')
1921 change_max_n_b_output =
get_log(.false.)
1922 if (change_max_n_b_output)
then
1923 call writo(
'Maximum number of modes to output?')
1924 max_n_b_output =
get_int(lim_lo=1)
1928 s_vj = [((kd-1._dp)/(
size(s_vj)-1),kd=1,
size(s_vj))]
1929 allocate(psi_t(grid_eq%n(3)))
1930 psi_t = eq_1%flux_t_E(:,0)/eq_1%flux_t_E(grid_eq%n(3),0)
1932 ierr =
spline(psi_t,eq_1%pres_E(:,id),s_vj,pres_vj(:,id),&
1938 pres_vj(:,1) = pres_vj(:,1) /
psi_0
1951 ierr =
spline(psi_t,eq_1%q_saf_E(:,0),s_vj,q_saf_vj,&
1954 ierr =
spline(psi_t,-eq_1%rot_t_E(:,0),s_vj,rot_t_v,&
1957 ierr =
spline(psi_t,i_tor*norm_transf(:,1),s_vj,i_tor_v,&
1960 ierr =
spline(psi_t,i_tor*norm_transf(:,2),s_vj,i_tor_j,&
1967 open(
hel_export_i,status=
'replace',file=trim(file_name),&
1969 chckerr(
'Failed to open file')
1971 write(
hel_export_i,
"(A)")
"!----- General Parameters -----"
1977 write(
hel_export_i,
"(A)")
"NS_ARRAY = 19, 39, 79, 159, 319"
1984 write(
hel_export_i,
"(A)")
"FTOL_ARRAY = 1.E-6, 1.E-6, 1.E-8, 1.E-10, &
1991 &trim(
r2str(-eq_1%flux_t_E(grid_eq%n(3),0)*
psi_0))
1994 &trim(
r2str(-eq_1%flux_t_E(grid_eq%n(3),0)))
1997 &trim(
r2str(i_tor_int(
size(i_tor_int))))
2002 write(
hel_export_i,
"(A)")
"!----- Pressure Parameters -----"
2003 write(
hel_export_i,
"(A)")
"PMASS_TYPE = 'cubic_spline'"
2004 write(
hel_export_i,
"(A)")
"GAMMA = 0.00000000000000E+00"
2007 do kd = 1,
size(s_vj)
2012 do kd = 1,
size(pres_vj,1)
2020 &
"!----- Current/Iota Parameters -----"
2025 write(
hel_export_i,
"(A)")
"PIOTA_TYPE = 'Cubic_spline'"
2028 do kd = 1,
size(s_vj)
2033 do kd = 1,
size(rot_t_v)
2039 write(
hel_export_i,
"(A)")
"PCURR_TYPE = 'cubic_spline_Ip'"
2042 do kd = 1,
size(s_vj)
2047 do kd = 1,
size(rot_t_v)
2055 &
"!----- Boundary Shape Parameters -----"
2056 ierr = print_mode_numbers(
hel_export_i,b_f(id_n_0,:,:,:),0,&
2061 if (jd.eq.id_n_0) cycle
2063 &b_f(jd,:,:,:),n_pert(jd),max_n_b_output)
2067 write(
hel_export_i,
"(A,ES23.16)")
"RAXIS = ", rz_b_0(1)
2068 write(
hel_export_i,
"(A,ES23.16)")
"ZAXIS = ", rz_b_0(2)
2074 file_name =
"flux_and_boundary_quantities_"//trim(
eq_name)//
'.jorek'
2075 open(
hel_export_i,status=
'replace',file=trim(file_name),&
2077 chckerr(
'Failed to open file')
2093 norm_j = [1._dp,1._dp,1._dp]
2095 write(
hel_export_i,
"('! ',A)")
'normalization factors:'
2096 write(
hel_export_i,
"('R_geo = ',ES23.16,' ! [m]')") norm_j(1)
2097 write(
hel_export_i,
"('Z_geo = ',ES23.16,' ! [m]')") norm_j(2)
2098 write(
hel_export_i,
"('F0 = ',ES23.16,' ! [Tm]')") norm_j(3)
2100 write(
hel_export_i,
"('! ',7(A23,' '))")
'pol. flux [Tm^2]', &
2101 &
'normalized pol. flux []',
'mu_0 pressure [T^2]',
'F [Tm]', &
2102 &
'-FFp [T^2 m^2/(Wb/rad)]',
'safety factor []',
'I_tor [A]'
2103 do kd = 1,
size(flux_j)
2104 write(
hel_export_i,
"(' ',7(ES23.16,' '))") flux_j(kd), &
2105 &flux_j(kd)/flux_j(
size(flux_j)),
mu_0_original*pres_vj(kd,0), &
2106 &f_vj(kd), -ffp_vj(kd), q_saf_vj(kd), i_tor_j(kd)
2109 write(
hel_export_i,*)
'! R [m] and Z [m] of boundary at psi [ ]'
2112 &', Z_boundary(',I4,') = ',ES23.16,&
2113 &', psi_boundary(',I4,') = ',ES23.16)") &
2114 &kd, bh_0(kd,1), kd, bh_0(kd,2), kd, 1._dp
2120 file_name =
"jorek_ffprime"
2121 open(
hel_export_i,status=
'replace',file=trim(file_name),&
2123 chckerr(
'Failed to open file')
2125 do kd = 1,
size(flux_j)
2127 &flux_j(kd)/flux_j(
size(flux_j)), -ffp_vj(kd)
2130 file_name =
"jorek_temperature"
2131 open(
hel_export_i,status=
'replace',file=trim(file_name),&
2133 chckerr(
'Failed to open file')
2136 do kd = 1,
size(flux_j)
2138 &flux_j(kd)/flux_j(
size(flux_j)),
mu_0_original*pres_vj(kd,0)
2152 subroutine plot_boundary(B,m_lims,n,plot_name,plot_dim,plot_lims)
2154 real(
dp),
intent(in) :: b(1:,0:,1:,1:)
2155 integer,
intent(in) :: m_lims(2)
2156 integer,
intent(in) :: n(:)
2157 character(len=*),
intent(in) :: plot_name
2158 integer,
intent(in) :: plot_dim(2)
2159 real(
dp),
intent(in) :: plot_lims(2,2)
2162 real(
dp),
allocatable :: ang_plot(:,:,:)
2163 real(
dp),
allocatable :: xyz_plot(:,:,:)
2168 allocate(ang_plot(plot_dim(1),plot_dim(2),2))
2169 allocate(xyz_plot(plot_dim(1),plot_dim(2),3))
2173 do id = 1,plot_dim(1)
2174 ang_plot(id,:,1) =
pi * (plot_lims(1,1) + (id-1._dp)/&
2175 &(plot_dim(1)-1)*(plot_lims(2,1)-plot_lims(1,1)))
2177 do kd = 1,plot_dim(2)
2178 ang_plot(:,kd,2) =
pi * (plot_lims(1,2) + (kd-1._dp)/&
2179 &(plot_dim(2)-1)*(plot_lims(2,2)-plot_lims(1,2)))
2184 do kd = 0,m_lims(2)-m_lims(1)
2188 xyz_plot(:,:,1) = xyz_plot(:,:,1) + &
2189 &b(id,kd,1,1)*cos(
m*ang_plot(:,:,1)-&
2190 &n(id)*ang_plot(:,:,2)) + &
2191 &b(id,kd,2,1)*sin(
m*ang_plot(:,:,1)-&
2192 &n(id)*ang_plot(:,:,2))
2194 xyz_plot(:,:,3) = xyz_plot(:,:,3) + &
2195 &b(id,kd,1,2)*cos(
m*ang_plot(:,:,1)-&
2196 &n(id)*ang_plot(:,:,2)) + &
2197 &b(id,kd,2,2)*sin(
m*ang_plot(:,:,1)-&
2198 &n(id)*ang_plot(:,:,2))
2203 call print_ex_2d([
'cross_section_'//trim(plot_name)],&
2204 &
'cross_section_'//trim(plot_name),xyz_plot(:,:,3),&
2205 &x=xyz_plot(:,:,1),draw=.false.)
2206 call draw_ex([
'cross_section_'//trim(plot_name)],&
2207 &
'cross_section_'//trim(plot_name),plot_dim(2),1,.false.)
2208 call print_ex_2d([
'cross_section_'//trim(plot_name)//
'_inv'],&
2209 &
'cross_section_'//trim(plot_name)//
'_inv',&
2210 &transpose(xyz_plot(:,:,3)),x=transpose(xyz_plot(:,:,1)),&
2212 call draw_ex([
'cross_section_'//trim(plot_name)//
'_inv'],&
2213 &
'cross_section_'//trim(plot_name)//
'_inv',plot_dim(1),1,&
2217 xyz_plot(:,:,2) = xyz_plot(:,:,1)*sin(ang_plot(:,:,2))
2218 xyz_plot(:,:,1) = xyz_plot(:,:,1)*cos(ang_plot(:,:,2))
2221 call print_ex_3d(
'plasma boundary',trim(plot_name),xyz_plot(:,:,3),&
2222 &x=xyz_plot(:,:,1),y=xyz_plot(:,:,2),draw=.false.)
2223 call draw_ex([
'plasma boundary'],trim(plot_name),1,2,.false.)
2224 end subroutine plot_boundary
2228 integer function print_mode_numbers(file_i,B,n_pert,max_n_B_output) &
2230 character(*),
parameter :: rout_name =
'print_mode_numbers'
2233 integer,
intent(in) :: file_i
2234 real(
dp),
intent(in) :: b(:,:,:)
2235 integer,
intent(in) :: n_pert
2236 integer,
intent(in) :: max_n_b_output
2241 character :: var_name
2242 character(len=2*max_str_ln) :: temp_output_str
2258 do m = 0,min(
size(b,1)-1,max_n_b_output)
2259 write(temp_output_str,
"(' ',A,'BC(',I4,', ',I4,') = ',&
2260 &ES23.16,' ',A,'BS(',I4,', ',I4,') = ',ES23.16)") &
2261 &var_name,n_pert,
m,b(
m+1,1,kd), &
2262 &var_name,n_pert,
m,b(
m+1,2,kd)
2263 write(unit=file_i,fmt=
'(A)',iostat=ierr) &
2264 &trim(temp_output_str)
2265 chckerr(
'Failed to write')
2267 write(unit=file_i,fmt=
"(A)",iostat=ierr)
""
2268 chckerr(
'Failed to write')
2270 end function print_mode_numbers
2271 end function create_vmec_input
2274 integer function print_output_eq_1(grid_eq,eq,data_name)
result(ierr)
2281 character(*),
parameter :: rout_name =
'print_output_eq_1'
2284 type(grid_type),
intent(in) :: grid_eq
2285 type(eq_1_type),
intent(in) :: eq
2286 character(len=*),
intent(in) :: data_name
2289 integer :: norm_id(2)
2290 type(
var_1d_type),
allocatable,
target :: eq_1d(:)
2292 type(grid_type) :: grid_trim
2300 call writo(
'Write flux equilibrium variables to output file')
2304 ierr =
trim_grid(grid_eq,grid_trim,norm_id)
2314 eq_1d_loc => eq_1d(id); id = id+1
2315 eq_1d_loc%var_name =
'pres_FD'
2316 allocate(eq_1d_loc%tot_i_min(2),eq_1d_loc%tot_i_max(2))
2317 allocate(eq_1d_loc%loc_i_min(2),eq_1d_loc%loc_i_max(2))
2318 eq_1d_loc%tot_i_min = [1,0]
2319 eq_1d_loc%tot_i_max = [grid_trim%n(3),
size(eq%pres_FD,2)-1]
2320 eq_1d_loc%loc_i_min = [grid_trim%i_min,0]
2321 eq_1d_loc%loc_i_max = [grid_trim%i_max,
size(eq%pres_FD,2)-1]
2322 loc_size =
size(eq%pres_FD(norm_id(1):norm_id(2),:))
2323 allocate(eq_1d_loc%p(loc_size))
2324 eq_1d_loc%p = reshape(eq%pres_FD(norm_id(1):norm_id(2),:),[loc_size])
2327 eq_1d_loc => eq_1d(id); id = id+1
2328 eq_1d_loc%var_name =
'q_saf_FD'
2329 allocate(eq_1d_loc%tot_i_min(2),eq_1d_loc%tot_i_max(2))
2330 allocate(eq_1d_loc%loc_i_min(2),eq_1d_loc%loc_i_max(2))
2331 eq_1d_loc%tot_i_min = [1,0]
2332 eq_1d_loc%tot_i_max = [grid_trim%n(3),
size(eq%q_saf_FD,2)-1]
2333 eq_1d_loc%loc_i_min = [grid_trim%i_min,0]
2334 eq_1d_loc%loc_i_max = [grid_trim%i_max,
size(eq%q_saf_FD,2)-1]
2335 loc_size =
size(eq%q_saf_FD(norm_id(1):norm_id(2),:))
2336 allocate(eq_1d_loc%p(loc_size))
2337 eq_1d_loc%p = reshape(eq%q_saf_FD(norm_id(1):norm_id(2),:),[loc_size])
2340 eq_1d_loc => eq_1d(id); id = id+1
2341 eq_1d_loc%var_name =
'rot_t_FD'
2342 allocate(eq_1d_loc%tot_i_min(2),eq_1d_loc%tot_i_max(2))
2343 allocate(eq_1d_loc%loc_i_min(2),eq_1d_loc%loc_i_max(2))
2344 eq_1d_loc%tot_i_min = [1,0]
2345 eq_1d_loc%tot_i_max = [grid_trim%n(3),
size(eq%rot_t_FD,2)-1]
2346 eq_1d_loc%loc_i_min = [grid_trim%i_min,0]
2347 eq_1d_loc%loc_i_max = [grid_trim%i_max,
size(eq%rot_t_FD,2)-1]
2348 loc_size =
size(eq%rot_t_FD(norm_id(1):norm_id(2),:))
2349 allocate(eq_1d_loc%p(loc_size))
2350 eq_1d_loc%p = reshape(eq%rot_t_FD(norm_id(1):norm_id(2),:),[loc_size])
2353 eq_1d_loc => eq_1d(id); id = id+1
2354 eq_1d_loc%var_name =
'flux_p_FD'
2355 allocate(eq_1d_loc%tot_i_min(2),eq_1d_loc%tot_i_max(2))
2356 allocate(eq_1d_loc%loc_i_min(2),eq_1d_loc%loc_i_max(2))
2357 eq_1d_loc%tot_i_min = [1,0]
2358 eq_1d_loc%tot_i_max = [grid_trim%n(3),
size(eq%flux_p_FD,2)-1]
2359 eq_1d_loc%loc_i_min = [grid_trim%i_min,0]
2360 eq_1d_loc%loc_i_max = [grid_trim%i_max,
size(eq%flux_p_FD,2)-1]
2361 loc_size =
size(eq%flux_p_FD(norm_id(1):norm_id(2),:))
2362 allocate(eq_1d_loc%p(loc_size))
2363 eq_1d_loc%p = reshape(eq%flux_p_FD(norm_id(1):norm_id(2),:),[loc_size])
2366 eq_1d_loc => eq_1d(id); id = id+1
2367 eq_1d_loc%var_name =
'flux_t_FD'
2368 allocate(eq_1d_loc%tot_i_min(2),eq_1d_loc%tot_i_max(2))
2369 allocate(eq_1d_loc%loc_i_min(2),eq_1d_loc%loc_i_max(2))
2370 eq_1d_loc%tot_i_min = [1,0]
2371 eq_1d_loc%tot_i_max = [grid_trim%n(3),
size(eq%flux_t_FD,2)-1]
2372 eq_1d_loc%loc_i_min = [grid_trim%i_min,0]
2373 eq_1d_loc%loc_i_max = [grid_trim%i_max,
size(eq%flux_t_FD,2)-1]
2374 loc_size =
size(eq%flux_t_FD(norm_id(1):norm_id(2),:))
2375 allocate(eq_1d_loc%p(loc_size))
2376 eq_1d_loc%p = reshape(eq%flux_t_FD(norm_id(1):norm_id(2),:),[loc_size])
2379 eq_1d_loc => eq_1d(id); id = id+1
2380 eq_1d_loc%var_name =
'rho'
2381 allocate(eq_1d_loc%tot_i_min(1),eq_1d_loc%tot_i_max(1))
2382 allocate(eq_1d_loc%loc_i_min(1),eq_1d_loc%loc_i_max(1))
2383 eq_1d_loc%tot_i_min = 1
2384 eq_1d_loc%tot_i_max = grid_trim%n(3)
2385 eq_1d_loc%loc_i_min = grid_trim%i_min
2386 eq_1d_loc%loc_i_max = grid_trim%i_max
2387 loc_size =
size(eq%rho(norm_id(1):norm_id(2)))
2388 allocate(eq_1d_loc%p(loc_size))
2389 eq_1d_loc%p = eq%rho(norm_id(1):norm_id(2))
2394 &ind_print=.not.grid_trim%divided)
2398 call grid_trim%dealloc()
2404 end function print_output_eq_1
2406 integer function print_output_eq_2(grid_eq,eq,data_name,rich_lvl,par_div,&
2407 &dealloc_vars)
result(ierr)
2414 character(*),
parameter :: rout_name =
'print_output_eq_2'
2417 type(grid_type),
intent(in) :: grid_eq
2418 type(eq_2_type),
intent(inout) :: eq
2419 character(len=*),
intent(in) :: data_name
2420 integer,
intent(in),
optional :: rich_lvl
2421 logical,
intent(in),
optional :: par_div
2422 logical,
intent(in),
optional :: dealloc_vars
2426 integer :: par_id(2)
2427 integer :: norm_id(2)
2428 type(
var_1d_type),
allocatable,
target :: eq_1d(:)
2430 type(grid_type) :: grid_trim
2432 logical :: par_div_loc
2434 logical :: dealloc_vars_loc
2440 call writo(
'Write metric equilibrium variables to output file')
2444 ierr =
trim_grid(grid_eq,grid_trim,norm_id)
2448 par_div_loc = .false.
2449 if (
present(par_div)) par_div_loc = par_div
2453 par_id = [1,n_tot(1)]
2454 if (grid_trim%n(1).gt.0 .and. par_div_loc)
then
2460 dealloc_vars_loc = .false.
2461 if (
present(dealloc_vars)) dealloc_vars_loc = dealloc_vars
2470 eq_1d_loc => eq_1d(id); id = id+1
2471 eq_1d_loc%var_name =
'g_FD'
2472 allocate(eq_1d_loc%tot_i_min(7),eq_1d_loc%tot_i_max(7))
2473 allocate(eq_1d_loc%loc_i_min(7),eq_1d_loc%loc_i_max(7))
2474 eq_1d_loc%tot_i_min = [1,1,1,1,0,0,0]
2475 eq_1d_loc%tot_i_max = [n_tot,6,
size(eq%g_FD,5)-1,
size(eq%g_FD,6)-1,&
2477 eq_1d_loc%loc_i_min = [par_id(1),1,grid_trim%i_min,1,0,0,0]
2478 eq_1d_loc%loc_i_max = [par_id(2),n_tot(2),grid_trim%i_max,6,&
2479 &
size(eq%g_FD,5)-1,
size(eq%g_FD,6)-1,
size(eq%g_FD,7)-1]
2480 loc_size =
size(eq%g_FD(:,:,norm_id(1):norm_id(2),:,:,:,:))
2481 allocate(eq_1d_loc%p(loc_size))
2482 eq_1d_loc%p = reshape(eq%g_FD(:,:,norm_id(1):norm_id(2),:,:,:,:),&
2484 if (dealloc_vars_loc)
deallocate(eq%g_FD)
2487 eq_1d_loc => eq_1d(id); id = id+1
2488 eq_1d_loc%var_name =
'h_FD'
2489 allocate(eq_1d_loc%tot_i_min(7),eq_1d_loc%tot_i_max(7))
2490 allocate(eq_1d_loc%loc_i_min(7),eq_1d_loc%loc_i_max(7))
2491 eq_1d_loc%tot_i_min = [1,1,1,1,0,0,0]
2492 eq_1d_loc%tot_i_max = [n_tot,6,
size(eq%h_FD,5)-1,
size(eq%h_FD,6)-1,&
2494 eq_1d_loc%loc_i_min = [par_id(1),1,grid_trim%i_min,1,0,0,0]
2495 eq_1d_loc%loc_i_max = [par_id(2),n_tot(2),grid_trim%i_max,6,&
2496 &
size(eq%h_FD,5)-1,
size(eq%h_FD,6)-1,
size(eq%h_FD,7)-1]
2497 loc_size =
size(eq%h_FD(:,:,norm_id(1):norm_id(2),:,:,:,:))
2498 allocate(eq_1d_loc%p(loc_size))
2499 eq_1d_loc%p = reshape(eq%h_FD(:,:,norm_id(1):norm_id(2),:,:,:,:),&
2501 if (dealloc_vars_loc)
deallocate(eq%h_FD)
2504 eq_1d_loc => eq_1d(id); id = id+1
2505 eq_1d_loc%var_name =
'jac_FD'
2506 allocate(eq_1d_loc%tot_i_min(6),eq_1d_loc%tot_i_max(6))
2507 allocate(eq_1d_loc%loc_i_min(6),eq_1d_loc%loc_i_max(6))
2508 eq_1d_loc%tot_i_min = [1,1,1,0,0,0]
2509 eq_1d_loc%tot_i_max = [n_tot,
size(eq%jac_FD,4)-1,
size(eq%jac_FD,5)-1,&
2510 &
size(eq%jac_FD,6)-1]
2511 eq_1d_loc%loc_i_min = [par_id(1),1,grid_trim%i_min,0,0,0]
2512 eq_1d_loc%loc_i_max = [par_id(2),n_tot(2),grid_trim%i_max,&
2513 &
size(eq%jac_FD,4)-1,
size(eq%jac_FD,5)-1,
size(eq%jac_FD,6)-1]
2514 loc_size =
size(eq%jac_FD(:,:,norm_id(1):norm_id(2),:,:,:))
2515 allocate(eq_1d_loc%p(loc_size))
2516 eq_1d_loc%p = reshape(eq%jac_FD(:,:,norm_id(1):norm_id(2),:,:,:),&
2518 if (dealloc_vars_loc)
deallocate(eq%jac_FD)
2521 eq_1d_loc => eq_1d(id); id = id+1
2522 eq_1d_loc%var_name =
'S'
2523 allocate(eq_1d_loc%tot_i_min(3),eq_1d_loc%tot_i_max(3))
2524 allocate(eq_1d_loc%loc_i_min(3),eq_1d_loc%loc_i_max(3))
2525 eq_1d_loc%tot_i_min = [1,1,1]
2526 eq_1d_loc%tot_i_max = n_tot
2527 eq_1d_loc%loc_i_min = [par_id(1),1,grid_trim%i_min]
2528 eq_1d_loc%loc_i_max = [par_id(2),n_tot(2),grid_trim%i_max]
2529 loc_size =
size(eq%S(:,:,norm_id(1):norm_id(2)))
2530 allocate(eq_1d_loc%p(loc_size))
2531 eq_1d_loc%p = reshape(eq%S(:,:,norm_id(1):norm_id(2)),&
2533 if (dealloc_vars_loc)
deallocate(eq%S)
2536 eq_1d_loc => eq_1d(id); id = id+1
2537 eq_1d_loc%var_name =
'kappa_n'
2538 allocate(eq_1d_loc%tot_i_min(3),eq_1d_loc%tot_i_max(3))
2539 allocate(eq_1d_loc%loc_i_min(3),eq_1d_loc%loc_i_max(3))
2540 eq_1d_loc%tot_i_min = [1,1,1]
2541 eq_1d_loc%tot_i_max = n_tot
2542 eq_1d_loc%loc_i_min = [par_id(1),1,grid_trim%i_min]
2543 eq_1d_loc%loc_i_max = [par_id(2),n_tot(2),grid_trim%i_max]
2544 loc_size =
size(eq%kappa_n(:,:,norm_id(1):norm_id(2)))
2545 allocate(eq_1d_loc%p(loc_size))
2546 eq_1d_loc%p = reshape(eq%kappa_n(:,:,norm_id(1):norm_id(2)),&
2548 if (dealloc_vars_loc)
deallocate(eq%kappa_n)
2551 eq_1d_loc => eq_1d(id); id = id+1
2552 eq_1d_loc%var_name =
'kappa_g'
2553 allocate(eq_1d_loc%tot_i_min(3),eq_1d_loc%tot_i_max(3))
2554 allocate(eq_1d_loc%loc_i_min(3),eq_1d_loc%loc_i_max(3))
2555 eq_1d_loc%tot_i_min = [1,1,1]
2556 eq_1d_loc%tot_i_max = n_tot
2557 eq_1d_loc%loc_i_min = [par_id(1),1,grid_trim%i_min]
2558 eq_1d_loc%loc_i_max = [par_id(2),n_tot(2),grid_trim%i_max]
2559 loc_size =
size(eq%kappa_g(:,:,norm_id(1):norm_id(2)))
2560 allocate(eq_1d_loc%p(loc_size))
2561 eq_1d_loc%p = reshape(eq%kappa_g(:,:,norm_id(1):norm_id(2)),&
2563 if (dealloc_vars_loc)
deallocate(eq%kappa_g)
2566 eq_1d_loc => eq_1d(id); id = id+1
2567 eq_1d_loc%var_name =
'sigma'
2568 allocate(eq_1d_loc%tot_i_min(3),eq_1d_loc%tot_i_max(3))
2569 allocate(eq_1d_loc%loc_i_min(3),eq_1d_loc%loc_i_max(3))
2570 eq_1d_loc%tot_i_min = [1,1,1]
2571 eq_1d_loc%tot_i_max = n_tot
2572 eq_1d_loc%loc_i_min = [par_id(1),1,grid_trim%i_min]
2573 eq_1d_loc%loc_i_max = [par_id(2),n_tot(2),grid_trim%i_max]
2574 loc_size =
size(eq%sigma(:,:,norm_id(1):norm_id(2)))
2575 allocate(eq_1d_loc%p(loc_size))
2576 eq_1d_loc%p = reshape(eq%sigma(:,:,norm_id(1):norm_id(2)),&
2578 if (dealloc_vars_loc)
deallocate(eq%sigma)
2582 &rich_lvl=rich_lvl,ind_print=.not.grid_trim%divided)
2586 call grid_trim%dealloc()
2592 end function print_output_eq_2
2595 integer function redistribute_output_eq_1(grid,grid_out,eq,eq_out) &
2599 character(*),
parameter :: rout_name =
'redistribute_output_eq_1'
2602 type(grid_type),
intent(in) :: grid
2603 type(grid_type),
intent(in) :: grid_out
2604 type(eq_1_type),
intent(in) :: eq
2605 type(eq_1_type),
intent(inout) :: eq_out
2614 call writo(
'Redistribute flux equilibrium variables')
2618 if (grid%n(1).ne.grid_out%n(1) .or. grid%n(2).ne.grid_out%n(2))
then
2620 chckerr(
'grid and grid_out are not compatible')
2624 call eq_out%init(grid_out,setup_e=.false.,setup_f=.true.)
2627 do id = 0,max_deriv+1
2630 &[grid%i_min,grid%i_max],[grid_out%i_min,grid_out%i_max])
2635 &[grid%i_min,grid%i_max],[grid_out%i_min,grid_out%i_max])
2640 &[grid%i_min,grid%i_max],[grid_out%i_min,grid_out%i_max])
2645 &[grid%i_min,grid%i_max],[grid_out%i_min,grid_out%i_max])
2650 &[grid%i_min,grid%i_max],[grid_out%i_min,grid_out%i_max])
2656 &[grid%i_min,grid%i_max],[grid_out%i_min,grid_out%i_max])
2661 end function redistribute_output_eq_1
2663 integer function redistribute_output_eq_2(grid,grid_out,eq,eq_out) &
2667 character(*),
parameter :: rout_name =
'redistribute_output_eq_2'
2670 type(grid_type),
intent(in) :: grid
2671 type(grid_type),
intent(in) :: grid_out
2672 type(eq_2_type),
intent(in) :: eq
2673 type(eq_2_type),
intent(inout) :: eq_out
2676 integer :: id, jd, kd, ld
2677 integer :: lims(2), lims_dis(2)
2678 integer :: siz(3), siz_dis(3)
2679 real(dp),
allocatable :: temp_var(:)
2685 call writo(
'Redistribute metric equilibrium variables')
2689 if (grid%n(1).ne.grid_out%n(1) .or. grid%n(2).ne.grid_out%n(2))
then
2691 chckerr(
'grid and grid_out are not compatible')
2695 call eq_out%init(grid_out,setup_e=.false.,setup_f=.true.)
2698 lims(1) = product(grid%n(1:2))*(grid%i_min-1)+1
2699 lims(2) = product(grid%n(1:2))*grid%i_max
2700 lims_dis(1) = product(grid%n(1:2))*(grid_out%i_min-1)+1
2701 lims_dis(2) = product(grid%n(1:2))*grid_out%i_max
2702 siz = [grid%n(1:2),grid%loc_n_r]
2703 siz_dis = [grid_out%n(1:2),grid_out%loc_n_r]
2704 allocate(temp_var(product(siz_dis)))
2710 do id = 1,
size(eq%g_FD,4)
2713 &eq%g_FD(:,:,:,id,jd,kd,ld),[product(siz)]),&
2714 &temp_var,lims,lims_dis)
2716 eq_out%g_FD(:,:,:,id,jd,kd,ld) = &
2717 &reshape(temp_var,siz_dis)
2721 &eq%h_FD(:,:,:,id,jd,kd,ld),[product(siz)]),&
2722 &temp_var,lims,lims_dis)
2724 eq_out%h_FD(:,:,:,id,jd,kd,ld) = &
2725 &reshape(temp_var,siz_dis)
2729 &eq%jac_FD(:,:,:,jd,kd,ld),[product(siz)]),&
2730 &temp_var,lims,lims_dis)
2732 eq_out%jac_FD(:,:,:,jd,kd,ld) = &
2733 &reshape(temp_var,siz_dis)
2742 eq_out%S = reshape(temp_var,siz_dis)
2748 eq_out%kappa_n = reshape(temp_var,siz_dis)
2754 eq_out%kappa_g = reshape(temp_var,siz_dis)
2760 eq_out%sigma = reshape(temp_var,siz_dis)
2764 end function redistribute_output_eq_2
2767 integer function calc_rzl_ind(grid_eq,eq,deriv)
result(ierr)
2771 character(*),
parameter :: rout_name =
'calc_RZL_ind'
2774 type(grid_type),
intent(in) :: grid_eq
2775 type(eq_2_type),
intent(inout) :: eq
2776 integer,
intent(in) :: deriv(3)
2783 ierr = check_deriv(deriv,max_deriv+1,
'calc_RZL')
2789 &
r_v_s(:,grid_eq%i_min:grid_eq%i_max,deriv(1)),&
2790 &grid_eq%trigon_factors,eq%R_E(:,:,:,deriv(1),deriv(2),deriv(3)),&
2791 &sym=[.true.,is_asym_v],deriv=[deriv(2),deriv(3)])
2794 &
z_v_s(:,grid_eq%i_min:grid_eq%i_max,deriv(1)),&
2795 &grid_eq%trigon_factors,eq%Z_E(:,:,:,deriv(1),deriv(2),deriv(3)),&
2796 &sym=[is_asym_v,.true.],deriv=[deriv(2),deriv(3)])
2799 &
l_v_s(:,grid_eq%i_min:grid_eq%i_max,deriv(1)),&
2800 &grid_eq%trigon_factors,eq%L_E(:,:,:,deriv(1),deriv(2),deriv(3)),&
2801 &sym=[is_asym_v,.true.],deriv=[deriv(2),deriv(3)])
2803 end function calc_rzl_ind
2805 integer function calc_rzl_arr(grid_eq,eq,deriv)
result(ierr)
2806 character(*),
parameter :: rout_name =
'calc_RZL_arr'
2809 type(grid_type),
intent(in) :: grid_eq
2810 type(eq_2_type),
intent(inout) :: eq
2811 integer,
intent(in) :: deriv(:,:)
2819 do id = 1,
size(deriv,2)
2820 ierr = calc_rzl_ind(grid_eq,eq,deriv(:,id))
2823 end function calc_rzl_arr
2826 integer function calc_g_c_ind(eq,deriv)
result(ierr)
2829 character(*),
parameter :: rout_name =
'calc_g_C_ind'
2832 type(eq_2_type),
intent(inout) :: eq
2833 integer,
intent(in) :: deriv(:)
2845 eq%g_C(:,:,:,:,deriv(1),deriv(2),deriv(3)) = 0._dp
2848 if (sum(deriv).eq.0)
then
2849 eq%g_C(:,:,:,
c([1,1],.true.),deriv(1),deriv(2),deriv(3)) = 1.0_dp
2850 eq%g_C(:,:,:,
c([3,3],.true.),deriv(1),deriv(2),deriv(3)) = 1.0_dp
2853 &eq%g_C(:,:,:,
c([2,2],.true.),deriv(1),deriv(2),deriv(3)),deriv)
2855 end function calc_g_c_ind
2857 integer function calc_g_c_arr(eq,deriv)
result(ierr)
2858 character(*),
parameter :: rout_name =
'calc_g_C_arr'
2861 type(eq_2_type),
intent(inout) :: eq
2862 integer,
intent(in) :: deriv(:,:)
2870 do id = 1,
size(deriv,2)
2871 ierr = calc_g_c_ind(eq,deriv(:,id))
2874 end function calc_g_c_arr
2877 integer function calc_g_v_ind(eq,deriv)
result(ierr)
2880 character(*),
parameter :: rout_name =
'calc_g_V_ind'
2883 type(eq_2_type),
intent(inout) :: eq
2884 integer,
intent(in) :: deriv(:)
2891 ierr = check_deriv(deriv,max_deriv,
'calc_g_V')
2895 ierr =
calc_g(eq%g_C,eq%T_VC,eq%g_E,deriv,max_deriv)
2897 end function calc_g_v_ind
2899 integer function calc_g_v_arr(eq,deriv)
result(ierr)
2900 character(*),
parameter :: rout_name =
'calc_g_V_arr'
2903 type(eq_2_type),
intent(inout) :: eq
2904 integer,
intent(in) :: deriv(:,:)
2912 do id = 1,
size(deriv,2)
2913 ierr = calc_g_v_ind(eq,deriv(:,id))
2916 end function calc_g_v_arr
2919 integer function calc_g_h_ind(grid_eq,eq,deriv)
result(ierr)
2926 character(*),
parameter :: rout_name =
'calc_g_H_ind'
2929 type(grid_type),
intent(in) :: grid_eq
2930 type(eq_2_type),
intent(inout) :: eq
2931 integer,
intent(in) :: deriv(:)
2934 type(ezspline2_r8) :: f_spl(2)
2935 character(len=max_str_ln) :: err_msg
2936 integer :: id, jd, kd, ld
2939 real(
dp) :: bcs_val(2,3)
2940 real(
dp),
allocatable :: rchi(:,:), rpsi(:,:)
2941 real(
dp),
allocatable :: zchi(:,:), zpsi(:,:)
2964 bc_ld = [1,2,0,1,0,1]
2967 eq%g_E(:,:,:,:,deriv(1),deriv(2),deriv(3)) = 0._dp
2969 if (sum(deriv).eq.0)
then
2971 allocate(rchi(
nchi,grid_eq%loc_n_r),rpsi(
nchi,grid_eq%loc_n_r))
2972 allocate(zchi(
nchi,grid_eq%loc_n_r),zpsi(
nchi,grid_eq%loc_n_r))
2975 do id = 1,grid_eq%n(1)
2976 ierr =
spline(grid_eq%loc_r_E,&
2977 &
r_h(id,grid_eq%i_min:grid_eq%i_max),grid_eq%loc_r_E,&
2979 &bcs_val=bcs_val(:,3))
2981 ierr =
spline(grid_eq%loc_r_E,&
2982 &
z_h(id,grid_eq%i_min:grid_eq%i_max),grid_eq%loc_r_E,&
2984 &bcs_val=bcs_val(:,3))
2989 do kd = 1,grid_eq%loc_n_r
2992 &bcs_val=bcs_val(:,1))
2996 &bcs_val=bcs_val(:,2))
3001 do jd = 1,grid_eq%n(2)
3002 eq%g_E(:,jd,:,
c([1,1],.true.),0,0,0) = rpsi*rpsi+zpsi*zpsi
3003 eq%g_E(:,jd,:,
c([1,2],.true.),0,0,0) = rpsi*rchi+zpsi*zchi
3004 eq%g_E(:,jd,:,
c([2,2],.true.),0,0,0) = rpsi*rchi+zpsi*zchi
3005 eq%g_E(:,jd,:,
c([2,2],.true.),0,0,0) = rchi*rchi+zchi*zchi
3006 eq%g_E(:,jd,:,
c([3,3],.true.),0,0,0) = &
3007 &1._dp/
h_h_33(:,grid_eq%i_min:grid_eq%i_max)
3011 deallocate(zchi,rchi,zpsi,rpsi)
3012 else if (deriv(3).ne.0)
then
3014 else if (deriv(1).le.2 .and. deriv(2).le.2)
then
3017 call ezspline_init(f_spl(ld),grid_eq%n(1),grid_eq%loc_n_r,&
3018 &bcs(:,ld),bcs(:,3),ierr)
3019 call ezspline_error(ierr)
3023 f_spl(ld)%x1 =
chi_h
3024 f_spl(ld)%x2 = grid_eq%loc_r_E
3027 f_spl(ld)%bcval1min = bcs_val(1,ld)
3028 f_spl(ld)%bcval1max = bcs_val(2,ld)
3029 f_spl(ld)%bcval2min = bcs_val(1,3)
3030 f_spl(ld)%bcval2max = bcs_val(2,3)
3034 do jd = 1,grid_eq%n(2)
3035 if (bc_ld(ld).eq.0) cycle
3038 call ezspline_setup(f_spl(bc_ld(ld)),&
3039 &eq%g_E(:,jd,:,ld,0,0,0),ierr,exact_dim=.true.)
3040 call ezspline_error(ierr)
3046 call ezspline_derivative(f_spl(bc_ld(ld)),deriv(2),&
3047 &deriv(1),grid_eq%n(1),grid_eq%loc_n_r,
chi_h,&
3048 &grid_eq%loc_r_E,eq%g_E(:,jd,:,ld,deriv(1),deriv(2),0),&
3050 call ezspline_error(ierr)
3057 call ezspline_free(f_spl(ld),ierr)
3058 call ezspline_error(ierr)
3062 err_msg =
'Derivative of order ('//trim(
i2str(deriv(1)))//
','//&
3063 &trim(
i2str(deriv(2)))//
','//trim(
i2str(deriv(3)))//
'&
3067 end function calc_g_h_ind
3069 integer function calc_g_h_arr(grid_eq,eq,deriv)
result(ierr)
3070 character(*),
parameter :: rout_name =
'calc_g_H_arr'
3073 type(grid_type),
intent(in) :: grid_eq
3074 type(eq_2_type),
intent(inout) :: eq
3075 integer,
intent(in) :: deriv(:,:)
3083 do id = 1,
size(deriv,2)
3084 ierr = calc_g_h_ind(grid_eq,eq,deriv(:,id))
3087 end function calc_g_h_arr
3090 integer function calc_g_f_ind(eq,deriv)
result(ierr)
3093 character(*),
parameter :: rout_name =
'calc_g_F_ind'
3096 type(eq_2_type),
intent(inout) :: eq
3097 integer,
intent(in) :: deriv(:)
3104 ierr = check_deriv(deriv,max_deriv,
'calc_g_F')
3109 ierr =
calc_g(eq%g_E,eq%T_FE,eq%g_F,deriv,max_deriv)
3111 end function calc_g_f_ind
3113 integer function calc_g_f_arr(eq,deriv)
result(ierr)
3114 character(*),
parameter :: rout_name =
'calc_g_F_arr'
3117 type(eq_2_type),
intent(inout) :: eq
3118 integer,
intent(in) :: deriv(:,:)
3126 do id = 1,
size(deriv,2)
3127 ierr = calc_g_f_ind(eq,deriv(:,id))
3130 end function calc_g_f_arr
3133 integer function calc_jac_v_ind(grid,eq,deriv)
result(ierr)
3137 character(*),
parameter :: rout_name =
'calc_jac_V_ind'
3140 type(grid_type),
intent(in) :: grid
3141 type(eq_2_type),
intent(inout) :: eq
3142 integer,
intent(in) :: deriv(:)
3149 ierr = check_deriv(deriv,max_deriv,
'calc_J_V')
3154 eq%jac_E(:,:,:,deriv(1),deriv(2),deriv(3)) = 0.0_dp
3158 &
jac_v_s(:,grid%i_min:grid%i_max,deriv(1)),grid%trigon_factors,&
3159 &eq%jac_E(:,:,:,deriv(1),deriv(2),deriv(3)),&
3160 &sym=[.true.,is_asym_v],deriv=[deriv(2),deriv(3)])
3162 end function calc_jac_v_ind
3164 integer function calc_jac_v_arr(grid,eq,deriv)
result(ierr)
3165 character(*),
parameter :: rout_name =
'calc_jac_V_arr'
3168 type(grid_type),
intent(in) :: grid
3169 type(eq_2_type),
intent(inout) :: eq
3170 integer,
intent(in) :: deriv(:,:)
3178 do id = 1,
size(deriv,2)
3179 ierr = calc_jac_v_ind(grid,eq,deriv(:,id))
3182 end function calc_jac_v_arr
3185 integer function calc_jac_h_ind(grid_eq,eq_1,eq_2,deriv)
result(ierr)
3190 character(*),
parameter :: rout_name =
'calc_jac_H_ind'
3193 type(grid_type),
intent(in) :: grid_eq
3194 type(eq_1_type),
intent(in) :: eq_1
3195 type(eq_2_type),
intent(inout) :: eq_2
3196 integer,
intent(in) :: deriv(:)
3199 character(len=max_str_ln) :: err_msg
3200 type(ezspline2_r8) :: f_spl
3203 real(dp) :: bcs_val(2,2)
3210 ierr = check_deriv(deriv,max_deriv,
'calc_jac_H')
3224 if (sum(deriv).eq.0)
then
3225 do kd = 1,grid_eq%loc_n_r
3226 do jd = 1,grid_eq%n(2)
3227 eq_2%jac_E(:,jd,kd,0,0,0) = eq_1%q_saf_E(kd,0)/&
3228 &(
h_h_33(:,grid_eq%i_min-1+kd)*&
3229 &
rbphi_h(grid_eq%i_min-1+kd,0))
3232 else if (deriv(3).ne.0)
then
3233 eq_2%jac_E(:,:,:,deriv(1),deriv(2),deriv(3)) = 0.0_dp
3234 else if (deriv(1).le.2 .and. deriv(2).le.2)
then
3236 call ezspline_init(f_spl,grid_eq%n(1),grid_eq%loc_n_r,&
3237 &bcs(:,1),bcs(:,2),ierr)
3238 call ezspline_error(ierr)
3244 f_spl%x2 = grid_eq%loc_r_E
3247 f_spl%bcval1min = bcs_val(1,1)
3248 f_spl%bcval1max = bcs_val(2,1)
3249 f_spl%bcval2min = bcs_val(1,2)
3250 f_spl%bcval2max = bcs_val(2,2)
3252 do jd = 1,grid_eq%n(2)
3254 call ezspline_setup(f_spl,eq_2%jac_E(:,jd,:,0,0,0),ierr,&
3256 call ezspline_error(ierr)
3262 call ezspline_derivative(f_spl,deriv(2),deriv(1),grid_eq%n(1),&
3263 &grid_eq%loc_n_r,
chi_h,grid_eq%loc_r_E,&
3264 &eq_2%jac_E(:,jd,:,deriv(1),deriv(2),0),ierr)
3265 call ezspline_error(ierr)
3270 call ezspline_free(f_spl,ierr)
3271 call ezspline_error(ierr)
3274 err_msg =
'Derivative of order ('//trim(
i2str(deriv(1)))//
','//&
3275 &trim(
i2str(deriv(2)))//
','//trim(
i2str(deriv(3)))//
'&
3279 end function calc_jac_h_ind
3281 integer function calc_jac_h_arr(grid_eq,eq_1,eq_2,deriv)
result(ierr)
3282 character(*),
parameter :: rout_name =
'calc_jac_H_arr'
3285 type(grid_type),
intent(in) :: grid_eq
3286 type(eq_1_type),
intent(in) :: eq_1
3287 type(eq_2_type),
intent(inout) :: eq_2
3288 integer,
intent(in) :: deriv(:,:)
3296 do id = 1,
size(deriv,2)
3297 ierr = calc_jac_h_ind(grid_eq,eq_1,eq_2,deriv(:,id))
3300 end function calc_jac_h_arr
3303 integer function calc_jac_f_ind(eq,deriv)
result(ierr)
3306 character(*),
parameter :: rout_name =
'calc_jac_F_ind'
3309 type(eq_2_type),
intent(inout) :: eq
3310 integer,
intent(in) :: deriv(:)
3322 eq%jac_F(:,:,:,deriv(1),deriv(2),deriv(3)) = 0.0_dp
3326 &eq%jac_F(:,:,:,deriv(1),deriv(2),deriv(3)),deriv)
3328 end function calc_jac_f_ind
3330 integer function calc_jac_f_arr(eq,deriv)
result(ierr)
3331 character(*),
parameter :: rout_name =
'calc_jac_F_arr'
3334 type(eq_2_type),
intent(inout) :: eq
3335 integer,
intent(in) :: deriv(:,:)
3343 do id = 1,
size(deriv,2)
3344 ierr = calc_jac_f_ind(eq,deriv(:,id))
3347 end function calc_jac_f_arr
3350 integer function calc_t_vc_ind(eq,deriv)
result(ierr)
3353 character(*),
parameter :: rout_name =
'calc_T_VC_ind'
3356 type(eq_2_type),
intent(inout) :: eq
3357 integer,
intent(in) :: deriv(:)
3369 eq%T_VC(:,:,:,:,deriv(1),deriv(2),deriv(3)) = 0._dp
3372 eq%T_VC(:,:,:,
c([1,1],.false.),deriv(1),deriv(2),deriv(3)) = &
3373 &eq%R_E(:,:,:,deriv(1)+1,deriv(2),deriv(3))
3375 eq%T_VC(:,:,:,
c([1,3],.false.),deriv(1),deriv(2),deriv(3)) = &
3376 &eq%Z_E(:,:,:,deriv(1)+1,deriv(2),deriv(3))
3377 eq%T_VC(:,:,:,
c([2,1],.false.),deriv(1),deriv(2),deriv(3)) = &
3378 &eq%R_E(:,:,:,deriv(1),deriv(2)+1,deriv(3))
3380 eq%T_VC(:,:,:,
c([2,3],.false.),deriv(1),deriv(2),deriv(3)) = &
3381 &eq%Z_E(:,:,:,deriv(1),deriv(2)+1,deriv(3))
3382 eq%T_VC(:,:,:,
c([3,1],.false.),deriv(1),deriv(2),deriv(3)) = &
3383 &eq%R_E(:,:,:,deriv(1),deriv(2),deriv(3)+1)
3384 if (sum(deriv).eq.0)
then
3385 eq%T_VC(:,:,:,
c([3,2],.false.),deriv(1),deriv(2),deriv(3)) = 1.0_dp
3389 eq%T_VC(:,:,:,
c([3,3],.false.),deriv(1),deriv(2),deriv(3)) = &
3390 &eq%Z_E(:,:,:,deriv(1),deriv(2),deriv(3)+1)
3391 end function calc_t_vc_ind
3393 integer function calc_t_vc_arr(eq,deriv)
result(ierr)
3394 character(*),
parameter :: rout_name =
'calc_T_VC_arr'
3397 type(eq_2_type),
intent(inout) :: eq
3398 integer,
intent(in) :: deriv(:,:)
3406 do id = 1,
size(deriv,2)
3407 ierr = calc_t_vc_ind(eq,deriv(:,id))
3410 end function calc_t_vc_arr
3413 integer function calc_t_vf_ind(grid_eq,eq_1,eq_2,deriv)
result(ierr)
3417 character(*),
parameter :: rout_name =
'calc_T_VF_ind'
3420 type(grid_type),
intent(in) :: grid_eq
3421 type(eq_1_type),
intent(in) :: eq_1
3422 type(eq_2_type),
intent(inout) :: eq_2
3423 integer,
intent(in) :: deriv(:)
3427 real(
dp),
allocatable :: theta_s(:,:,:,:,:,:)
3428 real(
dp),
allocatable :: zeta_s(:,:,:,:,:,:)
3442 dims = [grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r]
3445 eq_2%T_EF(:,:,:,:,deriv(1),deriv(2),deriv(3)) = 0._dp
3448 allocate(theta_s(dims(1),dims(2),dims(3),&
3449 &0:deriv(1)+1,0:deriv(2)+1,0:deriv(3)+1))
3452 theta_s(:,:,:,0,0,0) = grid_eq%theta_E
3453 theta_s(:,:,:,0,1,0) = 1.0_dp
3455 theta_s = theta_s + &
3456 &eq_2%L_E(:,:,:,0:deriv(1)+1,0:deriv(2)+1,0:deriv(3)+1)
3462 &eq_2%T_EF(:,:,:,
c([1,1],.false.),deriv(1),deriv(2),deriv(3)),&
3465 ierr =
add_arr_mult(eq_2%L_E(:,:,:,1:,0:,0:),eq_1%q_saf_E(:,0:),&
3466 &eq_2%T_EF(:,:,:,
c([1,1],.false.),deriv(1),deriv(2),deriv(3)),&
3470 if (deriv(2).eq.0 .and. deriv(3).eq.0)
then
3472 eq_2%T_EF(:,:,kd,
c([1,2],.false.),deriv(1),0,0) = &
3473 &eq_1%flux_p_E(kd,deriv(1)+1)/(2*
pi)
3480 eq_2%T_EF(:,:,:,
c([1,3],.false.),deriv(1),deriv(2),deriv(3)) = &
3481 &eq_2%L_E(:,:,:,deriv(1)+1,deriv(2),deriv(3))
3483 ierr =
add_arr_mult(theta_s(:,:,:,0:,1:,0:),eq_1%q_saf_E,&
3484 &eq_2%T_EF(:,:,:,
c([2,1],.false.),deriv(1),deriv(2),deriv(3)),&
3491 eq_2%T_EF(:,:,:,
c([2,3],.false.),deriv(1),deriv(2),deriv(3)) = &
3492 &theta_s(:,:,:,deriv(1),deriv(2)+1,deriv(3))
3494 if (sum(deriv).eq.0)
then
3495 eq_2%T_EF(:,:,:,
c([3,1],.false.),0,0,0) = -1.0_dp
3497 ierr =
add_arr_mult(eq_2%L_E(:,:,:,0:,0:,1:),eq_1%q_saf_E,&
3498 &eq_2%T_EF(:,:,:,
c([3,1],.false.),deriv(1),deriv(2),deriv(3)),&
3505 eq_2%T_EF(:,:,:,
c([3,3],.false.),deriv(1),deriv(2),deriv(3)) = &
3506 &eq_2%L_E(:,:,:,deriv(1),deriv(2),deriv(3)+1)
3509 eq_2%det_T_EF(:,:,:,deriv(1),deriv(2),deriv(3)) = 0.0_dp
3511 &eq_1%flux_p_E(:,1:)/(2*
pi),&
3512 &eq_2%det_T_EF(:,:,:,deriv(1),deriv(2),deriv(3)),deriv)
3516 allocate(zeta_s(dims(1),dims(2),dims(3),0:deriv(1)+1,&
3517 &0:deriv(2)+1,0:deriv(3)+1))
3520 zeta_s(:,:,:,0,0,0) = grid_eq%zeta_E
3521 zeta_s(:,:,:,0,0,1) = 1.0_dp
3525 eq_2%T_EF(:,:,:,
c([1,1],.false.),deriv(1),deriv(2),deriv(3)) = &
3526 &-theta_s(:,:,:,deriv(1)+1,deriv(2),deriv(3))
3528 &eq_2%T_EF(:,:,:,
c([1,1],.false.),deriv(1),deriv(2),deriv(3)),&
3532 if (deriv(2).eq.0 .and. deriv(3).eq.0)
then
3534 eq_2%T_EF(:,:,kd,
c([1,2],.false.),deriv(1),0,0) = &
3535 &-eq_1%flux_t_E(kd,deriv(1)+1)/(2*
pi)
3545 eq_2%T_EF(:,:,:,
c([2,1],.false.),deriv(1),deriv(2),deriv(3)) = &
3546 &-theta_s(:,:,:,deriv(1),deriv(2)+1,deriv(3))
3554 if (deriv(2).eq.0 .and. deriv(3).eq.0)
then
3556 eq_2%T_EF(:,:,kd,
c([3,1],.false.),deriv(1),0,0) = &
3557 &eq_1%rot_t_E(kd,deriv(1))
3560 c1 =
c([3,1],.false.)
3561 eq_2%T_EF(:,:,:,c1,deriv(1),deriv(2),deriv(3)) = &
3562 &eq_2%T_EF(:,:,:,c1,deriv(1),deriv(2),deriv(3)) &
3563 &-theta_s(:,:,:,deriv(1),deriv(2),deriv(3)+1)
3568 if (sum(deriv).eq.0)
then
3569 eq_2%T_EF(:,:,:,
c([3,3],.false.),0,0,0) = -1.0_dp
3576 eq_2%det_T_EF(:,:,:,deriv(1),deriv(2),deriv(3)) = 0.0_dp
3578 &eq_1%flux_t_E(:,1:)/(2*
pi),&
3579 &eq_2%det_T_EF(:,:,:,deriv(1),deriv(2),deriv(3)),deriv)
3582 end function calc_t_vf_ind
3584 integer function calc_t_vf_arr(grid_eq,eq_1,eq_2,deriv)
result(ierr)
3585 character(*),
parameter :: rout_name =
'calc_T_VF_arr'
3588 type(grid_type),
intent(in) :: grid_eq
3589 type(eq_1_type),
intent(in) :: eq_1
3590 type(eq_2_type),
intent(inout) :: eq_2
3591 integer,
intent(in) :: deriv(:,:)
3599 do id = 1,
size(deriv,2)
3600 ierr = calc_t_vf_ind(grid_eq,eq_1,eq_2,deriv(:,id))
3603 end function calc_t_vf_arr
3606 integer function calc_t_hf_ind(grid_eq,eq_1,eq_2,deriv)
result(ierr)
3610 character(*),
parameter :: rout_name =
'calc_T_HF_ind'
3613 type(grid_type),
intent(in) :: grid_eq
3614 type(eq_1_type),
intent(in) :: eq_1
3615 type(eq_2_type),
intent(inout) :: eq_2
3616 integer,
intent(in) :: deriv(:)
3620 real(
dp),
allocatable :: theta_s(:,:,:,:,:,:)
3621 real(
dp),
allocatable :: zeta_s(:,:,:,:,:,:)
3634 dims = [grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r]
3637 eq_2%T_EF(:,:,:,:,deriv(1),deriv(2),deriv(3)) = 0._dp
3640 allocate(theta_s(dims(1),dims(2),dims(3),&
3641 &0:deriv(1)+1,0:deriv(2)+1,0:deriv(3)+1))
3643 theta_s(:,:,:,0,0,0) = grid_eq%theta_E
3644 theta_s(:,:,:,0,1,0) = 1.0_dp
3650 &eq_2%T_EF(:,:,:,
c([1,1],.false.),deriv(1),deriv(2),deriv(3)),&
3654 if (sum(deriv).eq.0)
then
3655 eq_2%T_EF(:,:,:,
c([1,2],.false.),0,0,0) = 1._dp
3664 if (deriv(2).eq.0 .and. deriv(3).eq.0)
then
3666 eq_2%T_EF(:,:,kd,
c([2,1],.false.),deriv(1),0,0) = &
3667 &-eq_1%q_saf_E(kd,deriv(1))
3677 if (sum(deriv).eq.0)
then
3678 eq_2%T_EF(:,:,:,
c([2,3],.false.),0,0,0) = 1.0_dp
3684 if (sum(deriv).eq.0)
then
3685 eq_2%T_EF(:,:,:,
c([3,1],.false.),0,0,0) = 1.0_dp
3698 eq_2%det_T_EF(:,:,:,deriv(1),deriv(2),deriv(3)) = 0.0_dp
3699 if (sum(deriv).eq.0)
then
3700 eq_2%det_T_EF(:,:,:,deriv(1),0,0) = 1._dp
3706 allocate(zeta_s(dims(1),dims(2),dims(3),&
3707 &0:deriv(1)+1,0:deriv(2)+1,0:deriv(3)+1))
3710 zeta_s(:,:,:,0,0,0) = grid_eq%zeta_E
3711 zeta_s(:,:,:,0,0,1) = 1.0_dp
3716 &eq_2%T_EF(:,:,:,
c([1,1],.false.),deriv(1),deriv(2),deriv(3)),&
3720 if (deriv(2).eq.0 .and. deriv(3).eq.0)
then
3722 eq_2%T_EF(:,:,kd,
c([1,2],.false.),deriv(1),0,0) = &
3723 &eq_1%q_saf_E(kd,deriv(1))
3733 if (sum(deriv).eq.0)
then
3734 eq_2%T_EF(:,:,:,
c([2,1],.false.),0,0,0) = -1.0_dp
3746 if (deriv(2).eq.0 .and. deriv(3).eq.0)
then
3748 eq_2%T_EF(:,:,kd,
c([3,1],.false.),deriv(1),0,0) = &
3749 &eq_1%rot_t_E(kd,deriv(1))
3759 if (sum(deriv).eq.0)
then
3760 eq_2%T_EF(:,:,:,
c([3,3],.false.),deriv(1),0,0) = 1.0_dp
3767 eq_2%det_T_EF(:,:,:,deriv(1),deriv(2),deriv(3)) = 0.0_dp
3768 if (deriv(2).eq.0 .and. deriv(3).eq.0)
then
3770 eq_2%det_T_EF(:,:,kd,deriv(1),0,0) = &
3771 &eq_1%q_saf_E(kd,deriv(1))
3777 end function calc_t_hf_ind
3779 integer function calc_t_hf_arr(grid_eq,eq_1,eq_2,deriv)
result(ierr)
3780 character(*),
parameter :: rout_name =
'calc_T_HF_arr'
3783 type(grid_type),
intent(in) :: grid_eq
3784 type(eq_1_type),
intent(in) :: eq_1
3785 type(eq_2_type),
intent(inout) :: eq_2
3786 integer,
intent(in) :: deriv(:,:)
3794 do id = 1,
size(deriv,2)
3795 ierr = calc_t_hf_ind(grid_eq,eq_1,eq_2,deriv(:,id))
3798 end function calc_t_hf_arr
3815 character(*),
parameter :: rout_name =
'flux_q_plot'
3818 type(grid_type),
intent(in) :: grid_eq
3822 integer :: norm_id(2)
3824 integer :: n_vars = 5
3825 character(len=max_str_ln),
allocatable :: plot_titles(:)
3826 type(grid_type) :: grid_trim
3827 real(
dp),
allocatable :: x_plot_2d(:,:)
3828 real(
dp),
allocatable :: y_plot_2d(:,:)
3837 call writo(
'Plotting flux quantities')
3842 ierr =
trim_grid(grid_eq,grid_trim,norm_id)
3846 allocate(plot_titles(n_vars))
3847 plot_titles(1) =
'safety factor []'
3848 plot_titles(2) =
'rotational transform []'
3849 plot_titles(3) =
'pressure [pa]'
3850 plot_titles(4) =
'poloidal flux [Tm^2]'
3851 plot_titles(5) =
'toroidal flux [Tm^2]'
3854 ierr = flux_q_plot_hdf5()
3858 ierr = flux_q_plot_ex()
3862 call grid_trim%dealloc()
3868 integer function flux_q_plot_hdf5()
result(ierr)
3874 character(*),
parameter :: rout_name =
'flux_q_plot_HDF5'
3878 real(
dp),
allocatable :: x_plot(:,:,:,:)
3879 real(
dp),
allocatable :: y_plot(:,:,:,:)
3880 real(
dp),
allocatable :: z_plot(:,:,:,:)
3881 real(
dp),
allocatable :: f_plot(:,:,:,:)
3882 integer :: plot_dim(4)
3883 integer :: plot_offset(4)
3884 type(grid_type) :: grid_plot
3885 character(len=max_str_ln) :: file_name
3891 file_name =
'flux_quantities'
3894 allocate(y_plot_2d(grid_trim%loc_n_r,n_vars))
3896 y_plot_2d(:,1) = eq%q_saf_FD(norm_id(1):norm_id(2),0)
3897 y_plot_2d(:,2) = eq%rot_t_FD(norm_id(1):norm_id(2),0)
3898 y_plot_2d(:,3) = eq%pres_FD(norm_id(1):norm_id(2),0)
3899 y_plot_2d(:,4) = eq%flux_p_FD(norm_id(1):norm_id(2),0)
3900 y_plot_2d(:,5) = eq%flux_t_FD(norm_id(1):norm_id(2),0)
3909 &grid_plot%trigon_factors)
3914 plot_dim = [grid_plot%n(1),grid_plot%n(2),grid_plot%n(3),n_vars]
3915 plot_offset = [0,0,grid_plot%i_min-1,n_vars]
3918 allocate(x_plot(grid_plot%n(1),grid_plot%n(2),grid_plot%loc_n_r,&
3920 allocate(y_plot(grid_plot%n(1),grid_plot%n(2),grid_plot%loc_n_r,&
3922 allocate(z_plot(grid_plot%n(1),grid_plot%n(2),grid_plot%loc_n_r,&
3924 allocate(f_plot(grid_plot%n(1),grid_plot%n(2),grid_plot%loc_n_r,&
3929 &y_plot(:,:,:,1),z_plot(:,:,:,1))
3932 x_plot(:,:,:,id) = x_plot(:,:,:,1)
3933 y_plot(:,:,:,id) = y_plot(:,:,:,1)
3934 z_plot(:,:,:,id) = z_plot(:,:,:,1)
3937 do kd = 1,grid_plot%loc_n_r
3938 f_plot(:,:,kd,1) = y_plot_2d(kd,1)
3939 f_plot(:,:,kd,2) = y_plot_2d(kd,2)
3940 f_plot(:,:,kd,3) = y_plot_2d(kd,3)
3941 f_plot(:,:,kd,4) = y_plot_2d(kd,4)
3942 f_plot(:,:,kd,5) = y_plot_2d(kd,5)
3947 f_plot(:,:,:,3) = f_plot(:,:,:,3)*
pres_0
3948 f_plot(:,:,:,4) = f_plot(:,:,:,4)*
psi_0
3949 f_plot(:,:,:,5) = f_plot(:,:,:,5)*
psi_0
3953 call plot_hdf5(plot_titles,file_name,f_plot,plot_dim,plot_offset,&
3954 &x_plot,y_plot,z_plot,col=1,descr=
'Flux quantities')
3957 deallocate(y_plot_2d)
3958 deallocate(x_plot,y_plot,z_plot,f_plot)
3959 call grid_plot%dealloc()
3960 end function flux_q_plot_hdf5
3964 integer function flux_q_plot_ex()
result(ierr)
3967 character(*),
parameter :: rout_name =
'flux_q_plot_ex'
3970 character(len=max_str_ln),
allocatable :: file_name(:)
3971 real(dp),
allocatable :: ser_var_loc(:)
3978 allocate(x_plot_2d(grid_trim%n(3),n_vars))
3979 allocate(y_plot_2d(grid_trim%n(3),n_vars))
3989 ierr =
get_ser_var(eq%pres_FD(norm_id(1):norm_id(2),0),ser_var_loc)
3991 if (rank.eq.0) y_plot_2d(:,3) = ser_var_loc
3992 ierr =
get_ser_var(eq%flux_p_FD(norm_id(1):norm_id(2),0),&
3995 if (rank.eq.0) y_plot_2d(:,4) = ser_var_loc
3996 ierr =
get_ser_var(eq%flux_t_FD(norm_id(1):norm_id(2),0),&
3999 if (rank.eq.0) y_plot_2d(:,5) = ser_var_loc
4002 if (use_normalization .and. rank.eq.0)
then
4003 y_plot_2d(:,3) = y_plot_2d(:,3)*
pres_0
4004 y_plot_2d(:,4) = y_plot_2d(:,4)*
psi_0
4005 y_plot_2d(:,5) = y_plot_2d(:,5)*
psi_0
4011 deallocate(ser_var_loc)
4014 allocate(file_name(2))
4015 file_name(1) =
'pres'
4016 file_name(2) =
'flux'
4019 x_plot_2d(:,1) = grid_trim%r_F*2*pi/
max_flux_f
4021 x_plot_2d(:,id) = x_plot_2d(:,1)
4026 call print_ex_2d(plot_titles(3),file_name(1),&
4027 &y_plot_2d(:,3),x_plot_2d(:,3),draw=.false.)
4029 call print_ex_2d(plot_titles(4:5),file_name(2),&
4030 &y_plot_2d(:,4:5),x_plot_2d(:,4:5),draw=.false.)
4033 call draw_ex(plot_titles(3:3),file_name(1),1,1,.false.)
4034 call draw_ex(plot_titles(4:5),file_name(2),2,1,.false.)
4037 call writo(
'Safety factor and rotational transform are plotted &
4038 &using "plot_resonance"')
4041 deallocate(x_plot_2d,y_plot_2d)
4043 end function flux_q_plot_ex
4292 character(*),
parameter :: rout_name =
'calc_derived_q'
4295 type(grid_type),
intent(in) :: grid_eq
4297 type(
eq_2_type),
intent(inout),
target :: eq_2
4300 integer :: id, jd, kd, ld
4303 real(
dp) :: bcs_val(2,2)
4304 real(
dp),
allocatable :: d1_epar(:,:,:,:)
4305 real(
dp),
allocatable :: d3_epar(:,:,:,:)
4306 real(
dp),
allocatable :: b_n(:,:,:,:)
4307 real(
dp),
allocatable :: b_g(:,:,:,:)
4308 real(
dp),
allocatable :: de(:,:,:,:,:,:)
4309 real(
dp),
allocatable :: d_de(:,:,:,:,:)
4310 real(
dp),
allocatable :: rchi(:,:,:)
4311 real(
dp),
allocatable :: zchi(:,:,:)
4317 call writo(
'Set up derived quantities in flux coordinates')
4322 allocate(de(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r,2,2,3))
4323 allocate(d_de(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r,2,3))
4324 allocate(b_n(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r,3))
4325 allocate(b_g(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r,3))
4326 allocate(d1_epar(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r,3))
4327 allocate(d3_epar(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r,3))
4332 call calc_derived_s_direct(eq_2,eq_2%S)
4336 call calc_derived_de_epar_vmec(grid_eq,eq_2,de,d_de,b_n,b_g)
4339 call calc_derived_sigma_hel(grid_eq,eq_2,eq_2%sigma)
4359 allocate(rchi(grid_eq%n(1),grid_eq%loc_n_r,0:2))
4360 allocate(zchi(grid_eq%n(1),grid_eq%loc_n_r,0:2))
4361 rchi(:,:,0) =
r_h(:,grid_eq%i_min:grid_eq%i_max)
4362 zchi(:,:,0) =
z_h(:,grid_eq%i_min:grid_eq%i_max)
4363 do kd = 1,grid_eq%loc_n_r
4364 kd_h = grid_eq%i_min-1+kd
4367 &ord=3,deriv=id,bcs=bcs(:,1),bcs_val=bcs_val(:,1))
4370 &ord=3,deriv=id,bcs=bcs(:,2),bcs_val=bcs_val(:,2))
4377 call calc_derived_de_epar_hel(grid_eq,eq_1,rchi,zchi,&
4383 call calc_derived_dc_epar(de,d_de,eq_2%T_FE(:,:,:,:,0,0,0),&
4390 call calc_derived_kappa_from_e(grid_eq,eq_1,eq_2,b_n,b_g,&
4391 eq_2%kappa_n,eq_2%kappa_g)
4396 call calc_derived_sigma_from_e(eq_2,b_n,eq_2%sigma)
4399 ierr = calc_derived_s_from_deriv_hel(grid_eq,eq_2,bcs(:,1),&
4400 &bcs_val(:,1),eq_2%S)
4410 ierr = plot_derived_q(grid_eq,eq_2)
4413 call writo(
'Testing consistency of derived quantities')
4417 ierr = test_sigma_with_kappa_g(grid_eq,eq_1,eq_2)
4421 ierr = test_kappa(grid_eq,eq_1,eq_2)
4425 ierr = test_sigma_vmec(grid_eq,eq_1,eq_2)
4428 call test_s_hel(grid_eq,eq_2,rchi,zchi)
4432 call writo(
'Testing done')
4439 subroutine calc_derived_s_direct(eq_2,S)
4444 type(
eq_2_type),
intent(in),
target :: eq_2
4445 real(
dp),
intent(out) :: s(:,:,:)
4448 real(
dp),
pointer :: j(:,:,:) => null()
4449 real(
dp),
pointer :: h12(:,:,:) => null()
4450 real(
dp),
pointer :: d3h12(:,:,:) => null()
4451 real(
dp),
pointer :: h22(:,:,:) => null()
4452 real(
dp),
pointer :: d3h22(:,:,:) => null()
4453 real(
dp),
allocatable :: h22_corrected(:,:,:)
4456 j => eq_2%jac_FD(:,:,:,0,0,0)
4457 h12 => eq_2%h_FD(:,:,:,
c([1,2],.true.),0,0,0)
4458 d3h12 => eq_2%h_FD(:,:,:,
c([1,2],.true.),0,0,1)
4459 h22 => eq_2%h_FD(:,:,:,
c([2,2],.true.),0,0,0)
4460 d3h22 => eq_2%h_FD(:,:,:,
c([2,2],.true.),0,0,1)
4461 allocate(h22_corrected(
size(s,1),
size(s,2),
size(s,3)))
4462 h22_corrected = sign(max(abs(h22),
tol_zero),h22)
4465 s = -(d3h12/h22_corrected - d3h22*h12/h22_corrected**2)/j
4468 nullify(j,h12,d3h12,h22,d3h22)
4469 end subroutine calc_derived_s_direct
4472 integer function calc_derived_s_from_deriv_hel(grid_eq,eq_2,bcs,&
4473 &bcs_val,S)
result(ierr)
4478 character(*),
parameter :: rout_name = &
4479 &
'calc_derived_S_from_deriv_HEL'
4482 type(grid_type),
intent(in) :: grid_eq
4483 type(
eq_2_type),
intent(in),
target :: eq_2
4485 real(dp) :: bcs_val(2)
4486 real(dp),
intent(out) :: s(:,:,:)
4490 real(dp),
pointer :: j(:,:,:) => null()
4491 real(dp),
pointer :: h12(:,:,:) => null()
4492 real(dp),
pointer :: h22(:,:,:) => null()
4498 j => eq_2%jac_FD(:,:,:,0,0,0)
4499 h12 => eq_2%h_FD(:,:,:,
c([1,2],.true.),0,0,0)
4500 h22 => eq_2%h_FD(:,:,:,
c([2,2],.true.),0,0,0)
4503 do kd = 1,grid_eq%loc_n_r
4504 do jd = 1,grid_eq%n(2)
4506 &s(:,jd,kd),ord=3,deriv=1,bcs=bcs,bcs_val=bcs_val)
4514 end function calc_derived_s_from_deriv_hel
4519 subroutine calc_derived_s_from_sigma_hel(grid_eq,eq_2,Rchi,Zchi,S)
4524 type(grid_type),
intent(in) :: grid_eq
4525 type(
eq_2_type),
intent(in),
target :: eq_2
4526 real(dp),
intent(in) :: rchi(:,:,0:)
4527 real(dp),
intent(in) :: zchi(:,:,0:)
4528 real(dp),
intent(out) :: s(:,:,:)
4533 real(dp),
allocatable :: k(:,:,:)
4534 real(dp),
pointer :: j(:,:,:) => null()
4535 real(dp),
pointer :: g33(:,:,:) => null()
4536 real(dp),
pointer :: h22(:,:,:) => null()
4539 j => eq_2%jac_FD(:,:,:,0,0,0)
4540 g33 => eq_2%g_FD(:,:,:,
c([3,3],.true.),0,0,0)
4541 h22 => eq_2%h_FD(:,:,:,
c([2,2],.true.),0,0,0)
4544 allocate(k(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r))
4545 do jd = 1,grid_eq%n(2)
4546 k(:,jd,:) = zchi(:,:,1)/rchi(:,:,0) + &
4547 &(zchi(:,:,1)*rchi(:,:,2)-rchi(:,:,1)*zchi(:,:,2))/&
4548 &(rchi(:,:,1)**2+zchi(:,:,1)**2)
4549 do kd = 1,grid_eq%loc_n_r
4550 kd_h = grid_eq%i_min-1+kd
4551 k(:,jd,kd) = -2._dp*
rbphi_h(kd_h,0)/rchi(:,kd,0) * &
4564 end subroutine calc_derived_s_from_sigma_hel
4568 subroutine calc_derived_dc_epar(de,D_de,T_FE,D1_epar,D3_epar)
4572 real(dp),
intent(in) :: de(:,:,:,:,:,:)
4573 real(dp),
intent(in) :: d_de(:,:,:,:,:)
4574 real(dp),
intent(in) :: t_fe(:,:,:,:)
4575 real(dp),
intent(out) :: d1_epar(:,:,:,:)
4576 real(dp),
intent(out) :: d3_epar(:,:,:,:)
4587 d1_epar(:,:,:,ld) = d1_epar(:,:,:,ld) + &
4588 &de(:,:,:,id,jd,ld) * &
4589 &t_fe(:,:,:,
c([1,1+id],.false.)) * &
4590 &t_fe(:,:,:,
c([3,1+jd],.false.))
4592 d3_epar(:,:,:,ld) = d3_epar(:,:,:,ld) + &
4593 &de(:,:,:,id,jd,ld) * &
4594 &t_fe(:,:,:,
c([3,1+id],.false.)) * &
4595 &t_fe(:,:,:,
c([3,1+jd],.false.))
4598 d1_epar(:,:,:,ld) = d1_epar(:,:,:,ld) + &
4599 &d_de(:,:,:,id,ld) * &
4600 &t_fe(:,:,:,
c([1,1+id],.false.))
4602 d3_epar(:,:,:,ld) = d3_epar(:,:,:,ld) + &
4603 &d_de(:,:,:,id,ld) * &
4604 &t_fe(:,:,:,
c([3,1+id],.false.))
4607 end subroutine calc_derived_dc_epar
4613 subroutine calc_derived_de_epar_vmec(grid_eq,eq_2,de,D_de,b_n,b_g)
4618 type(grid_type),
intent(in) :: grid_eq
4620 real(
dp),
intent(out) :: de(:,:,:,:,:,:)
4621 real(
dp),
intent(out) :: d_de(:,:,:,:,:)
4622 real(
dp),
intent(out) :: b_n(:,:,:,:)
4623 real(
dp),
intent(out) :: b_g(:,:,:,:)
4635 de(:,:,:,1,1,1) = eq_2%R_E(:,:,:,0,2,0)
4637 de(:,:,:,1,1,3) = eq_2%Z_E(:,:,:,0,2,0)
4640 de(:,:,:,2,1,1) = eq_2%R_E(:,:,:,0,1,1)
4641 de(:,:,:,2,1,2) = eq_2%R_E(:,:,:,0,1,0)/eq_2%R_E(:,:,:,0,0,0)
4642 de(:,:,:,2,1,3) = eq_2%Z_E(:,:,:,0,1,1)
4645 de(:,:,:,1,2,1) = eq_2%R_E(:,:,:,0,1,1)
4646 de(:,:,:,1,2,2) = eq_2%R_E(:,:,:,0,1,0)/eq_2%R_E(:,:,:,0,0,0)
4647 de(:,:,:,1,2,3) = eq_2%Z_E(:,:,:,0,1,1)
4650 de(:,:,:,2,2,1) = eq_2%R_E(:,:,:,0,0,2)-eq_2%R_E(:,:,:,0,0,0)
4651 de(:,:,:,2,2,2) = 2*eq_2%R_E(:,:,:,0,0,1)/eq_2%R_E(:,:,:,0,0,0)
4652 de(:,:,:,2,2,3) = eq_2%Z_E(:,:,:,0,0,2)
4658 d_de(:,:,:,1,1) = eq_2%R_E(:,:,:,0,1,0)*&
4659 &eq_2%T_FE(:,:,:,
c([3,2],.false.),0,1,0) + &
4660 &eq_2%R_E(:,:,:,0,0,1)*&
4661 &eq_2%T_FE(:,:,:,
c([3,3],.false.),0,1,0)
4662 d_de(:,:,:,1,2) = eq_2%T_FE(:,:,:,
c([3,3],.false.),0,1,0)
4663 d_de(:,:,:,1,3) = eq_2%Z_E(:,:,:,0,1,0)*&
4664 &eq_2%T_FE(:,:,:,
c([3,2],.false.),0,1,0) + &
4665 &eq_2%Z_E(:,:,:,0,0,1)*&
4666 &eq_2%T_FE(:,:,:,
c([3,3],.false.),0,1,0)
4669 d_de(:,:,:,2,1) = eq_2%R_E(:,:,:,0,1,0)*&
4670 &eq_2%T_FE(:,:,:,
c([3,2],.false.),0,0,1) + &
4671 &eq_2%R_E(:,:,:,0,0,1)*&
4672 &eq_2%T_FE(:,:,:,
c([3,3],.false.),0,0,1)
4673 d_de(:,:,:,2,2) = eq_2%T_FE(:,:,:,
c([3,3],.false.),0,0,1)
4674 d_de(:,:,:,2,3) = eq_2%Z_E(:,:,:,0,1,0)*&
4675 &eq_2%T_FE(:,:,:,
c([3,2],.false.),0,0,1) + &
4676 &eq_2%Z_E(:,:,:,0,0,1)*&
4677 &eq_2%T_FE(:,:,:,
c([3,3],.false.),0,0,1)
4683 b_n(:,:,:,2) = (eq_2%R_E(:,:,:,0,0,1)*eq_2%Z_E(:,:,:,0,1,0) - &
4684 &eq_2%R_E(:,:,:,0,1,0)*eq_2%Z_E(:,:,:,0,0,1))
4686 &-eq_2%jac_FD(:,:,:,0,0,0) * (1+eq_2%L_E(:,:,:,0,1,0)) / &
4687 &max(
tol_zero, ( (b_n(:,:,:,2)/eq_2%R_E(:,:,:,0,0,0))**2 + &
4688 &eq_2%R_E(:,:,:,0,1,0)**2 + eq_2%Z_E(:,:,:,0,1,0)**2 )) / &
4689 &eq_2%R_E(:,:,:,0,0,0)
4690 b_n(:,:,:,2) = b_n(:,:,:,1) * b_n(:,:,:,2)
4691 b_n(:,:,:,3) = b_n(:,:,:,1) * eq_2%R_E(:,:,:,0,1,0)
4692 b_n(:,:,:,1) = -b_n(:,:,:,1) * eq_2%Z_E(:,:,:,0,1,0)
4695 do kd = 1,grid_eq%loc_n_r
4696 b_g(:,:,kd,2) = eq_2%g_FD(:,:,kd,
c([3,3],.true.),0,0,0) - &
4697 &eq_2%g_FD(:,:,kd,
c([1,3],.true.),0,0,0) * &
4698 &eq_1%q_saf_FD(kd,0)
4699 b_g(:,:,kd,1) = (b_g(:,:,kd,2)*eq_2%L_E(:,:,kd,0,0,1)-&
4700 &eq_2%g_FD(:,:,kd,
c([1,3],.true.),0,0,0)) / &
4701 &(1+eq_2%L_E(:,:,kd,0,1,0))
4703 b_g(:,:,:,3) = b_g(:,:,:,1)*eq_2%Z_E(:,:,:,0,1,0) - &
4704 &b_g(:,:,:,2)*eq_2%Z_E(:,:,:,0,0,1)
4705 b_g(:,:,:,1) = b_g(:,:,:,1)*eq_2%R_E(:,:,:,0,1,0) - &
4706 &b_g(:,:,:,2)*eq_2%R_E(:,:,:,0,0,1)
4707 b_g(:,:,:,2) = -b_g(:,:,:,2)*eq_2%R_E(:,:,:,0,0,0)**2
4709 b_g(:,:,:,ld) = b_g(:,:,:,ld) / &
4710 &eq_2%g_FD(:,:,:,
c([3,3],.true.),0,0,0)
4712 end subroutine calc_derived_de_epar_vmec
4713 subroutine calc_derived_de_epar_hel(grid_eq,eq_1,Rchi,Zchi,&
4719 type(grid_type),
intent(in) :: grid_eq
4721 real(dp),
intent(in) :: rchi(:,:,0:)
4722 real(dp),
intent(in) :: zchi(:,:,0:)
4723 real(dp),
intent(out) :: de(:,:,:,:,:,:)
4724 real(dp),
intent(out) :: d_de(:,:,:,:,:)
4725 real(dp),
intent(out) :: b_n(:,:,:,:)
4726 real(dp),
intent(out) :: b_g(:,:,:,:)
4730 real(dp),
allocatable :: dum1(:,:)
4731 real(dp),
allocatable :: dum2(:,:)
4743 de(:,1,:,1,1,1) = rchi(:,:,2)
4745 de(:,1,:,1,1,3) = zchi(:,:,2)
4749 de(:,1,:,2,1,2) = -rchi(:,:,1)/rchi(:,:,0)
4754 de(:,1,:,1,2,2) = -rchi(:,:,1)/rchi(:,:,0)
4758 de(:,1,:,2,2,1) = -rchi(:,:,0)
4764 do jd = 1,grid_eq%n(2)
4766 allocate(dum1(grid_eq%n(1),grid_eq%loc_n_r))
4767 allocate(dum2(grid_eq%n(1),grid_eq%loc_n_r))
4768 dum1 = rchi(:,:,1)**2 + zchi(:,:,1)**2
4769 do kd = 1,grid_eq%loc_n_r
4770 dum2(:,kd) = dum1(:,kd) + &
4771 &(eq_1%q_saf_FD(kd,0)*rchi(:,kd,0))**2
4775 b_n(:,jd,:,1) = rchi(:,:,0)*zchi(:,:,1)/dum1
4776 b_n(:,jd,:,3) = -rchi(:,:,0)*rchi(:,:,1)/dum1
4779 b_g(:,jd,:,1) = -rchi(:,:,0)**2 * rchi(:,:,1)/dum2
4780 b_g(:,jd,:,2) = -rchi(:,:,0)**2 * dum1/dum2
4781 b_g(:,jd,:,3) = -rchi(:,:,0)**2 * zchi(:,:,1)/dum2
4784 deallocate(dum1,dum2)
4786 do kd = 1,grid_eq%loc_n_r
4787 kd_h = grid_eq%i_min-1+kd
4788 b_n(:,:,kd,:) = b_n(:,:,kd,:) * &
4789 &eq_1%q_saf_FD(kd,0)/
rbphi_h(kd_h,0)
4790 b_g(:,:,kd,1:3:2) = b_g(:,:,kd,1:3:2) * eq_1%q_saf_FD(kd,0)
4792 end subroutine calc_derived_de_epar_hel
4796 subroutine calc_derived_kappa_from_e(grid_eq,eq_1,eq_2,b_n,b_g,&
4802 type(grid_type),
intent(in) :: grid_eq
4805 real(dp),
intent(in) :: b_n(:,:,:,:)
4806 real(dp),
intent(in) :: b_g(:,:,:,:)
4807 real(dp),
intent(out) :: kappa_n(:,:,:)
4808 real(dp),
intent(out) :: kappa_g(:,:,:)
4814 kappa_n = kappa_n + d3_epar(:,:,:,ld) * b_n(:,:,:,ld)
4815 kappa_g = kappa_g + d3_epar(:,:,:,ld) * b_g(:,:,:,ld)
4819 kappa_n = kappa_n / eq_2%g_FD(:,:,:,
c([3,3],.true.),0,0,0)
4820 kappa_g = kappa_g / eq_2%g_FD(:,:,:,
c([3,3],.true.),0,0,0)
4823 if (.not.use_pol_flux_f)
then
4824 do kd = 1,grid_eq%loc_n_r
4825 kappa_n(:,:,kd) = kappa_n(:,:,kd) * eq_1%q_saf_E(kd,0)
4826 kappa_g(:,:,kd) = kappa_g(:,:,kd) / eq_1%q_saf_E(kd,0)
4829 end subroutine calc_derived_kappa_from_e
4832 subroutine calc_derived_sigma_from_e(eq_2,b_n,sigma)
4837 real(dp),
intent(in) :: b_n(:,:,:,:)
4838 real(dp),
intent(out) :: sigma(:,:,:)
4843 sigma = sigma + 2._dp * b_n(:,:,:,ld) * &
4844 &(d3_epar(:,:,:,ld) * &
4845 &eq_2%g_FD(:,:,:,
c([1,3],.true.),0,0,0) - &
4846 &d1_epar(:,:,:,ld) * &
4847 &eq_2%g_FD(:,:,:,
c([3,3],.true.),0,0,0)) / &
4848 &eq_2%jac_FD(:,:,:,0,0,0)**2
4852 sigma = -(sigma + eq_2%jac_FD(:,:,:,0,0,0) * &
4853 &eq_2%h_FD(:,:,:,
c([2,2],.true.),0,0,0)*eq_2%S) * &
4854 &eq_2%jac_FD(:,:,:,0,0,0)/&
4855 &(
vac_perm*eq_2%g_FD(:,:,:,
c([3,3],.true.),0,0,0))
4856 end subroutine calc_derived_sigma_from_e
4859 subroutine calc_derived_sigma_hel(grid_eq,eq_2,sigma)
4864 type(grid_type),
intent(in) :: grid_eq
4865 type(
eq_2_type),
intent(in),
target :: eq_2
4866 real(dp),
intent(out) :: sigma(:,:,:)
4871 real(dp),
pointer :: j(:,:,:) => null()
4872 real(dp),
pointer :: g13(:,:,:) => null()
4873 real(dp),
pointer :: g33(:,:,:) => null()
4876 j => eq_2%jac_FD(:,:,:,0,0,0)
4877 g13 => eq_2%g_FD(:,:,:,
c([1,3],.true.),0,0,0)
4878 g33 => eq_2%g_FD(:,:,:,
c([3,3],.true.),0,0,0)
4881 do kd = 1,grid_eq%loc_n_r
4882 kd_h = grid_eq%i_min-1+kd
4884 &eq_1%pres_FD(kd,1)*j(:,:,kd)*g13(:,:,kd)/g33(:,:,kd)
4889 end subroutine calc_derived_sigma_hel
4893 integer function plot_derived_q(grid_eq,eq_2)
result(ierr)
4900 character(*),
parameter :: rout_name =
'plot_derived_q'
4908 real(
dp) :: norm_factors(4)
4909 real(
dp),
allocatable :: x_plot(:,:,:)
4910 real(
dp),
allocatable :: y_plot(:,:,:)
4911 real(
dp),
allocatable :: z_plot(:,:,:)
4912 real(
dp),
pointer :: ang_par_f(:,:,:) => null()
4913 integer :: norm_id(2)
4914 integer :: plot_dim(3)
4915 integer :: plot_offset(3)
4916 logical :: cont_plot
4921 call writo(
'Plotting derived equilibrium quantities')
4925 ierr =
trim_grid(grid_eq,grid_trim,norm_id)
4929 allocate(x_plot(grid_trim%n(1),grid_trim%n(2),grid_trim%loc_n_r))
4930 allocate(y_plot(grid_trim%n(1),grid_trim%n(2),grid_trim%loc_n_r))
4931 allocate(z_plot(grid_trim%n(1),grid_trim%n(2),grid_trim%loc_n_r))
4935 ang_par_f => grid_eq%theta_F
4937 ang_par_f => grid_eq%zeta_F
4945 x_plot = ang_par_f(:,:,norm_id(1):norm_id(2))
4946 do jd = 1,grid_trim%n(2)
4947 y_plot(:,jd,:) =
alpha(jd)
4949 do kd = 1,grid_trim%loc_n_r
4959 ierr =
calc_xyz_grid(grid_eq,grid_trim,x_plot,y_plot,z_plot)
4972 plot_dim = grid_trim%n
4973 plot_offset = [0,0,grid_trim%i_min-1]
4978 norm_factors = 1._dp
4981 norm_factors(2) = 1._dp/(
r_0**3)
4982 norm_factors(3) = 1._dp/(
r_0*
b_0)
4983 norm_factors(4) = 1._dp
4987 call plot_hdf5(
'sigma',
'TEST_sigma',&
4988 &eq_2%sigma(:,:,norm_id(1):norm_id(2))*norm_factors(1),&
4989 &tot_dim=plot_dim,loc_offset=plot_offset,cont_plot=cont_plot,&
4990 &x=x_plot,y=y_plot,z=z_plot)
4993 call plot_hdf5(
'shear',
'TEST_S',&
4994 &eq_2%S(:,:,norm_id(1):norm_id(2))*norm_factors(2),&
4995 &tot_dim=plot_dim,loc_offset=plot_offset,cont_plot=cont_plot,&
4996 &x=x_plot,y=y_plot,z=z_plot)
4999 call plot_hdf5(
'kappa_n',
'TEST_kappa_n',&
5000 &eq_2%kappa_n(:,:,norm_id(1):norm_id(2))*norm_factors(3),&
5001 &tot_dim=plot_dim,loc_offset=plot_offset,cont_plot=cont_plot,&
5002 &x=x_plot,y=y_plot,z=z_plot)
5005 call plot_hdf5(
'kappa_g',
'TEST_kappa_g',&
5006 &eq_2%kappa_g(:,:,norm_id(1):norm_id(2))*norm_factors(4),&
5007 &tot_dim=plot_dim,loc_offset=plot_offset,cont_plot=cont_plot,&
5008 &x=x_plot,y=y_plot,z=z_plot)
5012 call grid_trim%dealloc()
5015 end function plot_derived_q
5018 integer function test_sigma_with_kappa_g(grid_eq,eq_1,eq_2)
result(ierr)
5022 character(*),
parameter :: rout_name =
'test_sigma_with_kappa_g'
5025 type(grid_type),
intent(in) :: grid_eq
5026 type(eq_1_type),
intent(in) :: eq_1
5027 type(eq_2_type),
intent(in) :: eq_2
5030 real(
dp),
allocatable :: sigma_alt(:,:,:)
5031 real(
dp),
allocatable :: d3sigma(:,:,:)
5032 real(
dp),
allocatable :: d3sigma_alt(:,:,:)
5033 real(
dp),
pointer :: ang_par_f(:,:,:) => null()
5038 call writo(
'Testing whether -2 p'' J kappa_g = D3sigma')
5043 ang_par_f => grid_eq%theta_F
5045 ang_par_f => grid_eq%zeta_F
5049 allocate(d3sigma(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r))
5050 allocate(d3sigma_alt(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r))
5051 allocate(sigma_alt(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r))
5054 do kd = 1,grid_eq%loc_n_r
5055 do jd = 1,grid_eq%n(2)
5056 ierr =
spline(ang_par_f(:,jd,kd),eq_2%sigma(:,jd,kd),&
5057 &ang_par_f(:,jd,kd),d3sigma(:,jd,kd),&
5064 do kd = 1,grid_eq%loc_n_r
5065 d3sigma_alt(:,:,kd) = -2*eq_1%pres_FD(kd,1)*&
5066 &eq_2%kappa_g(:,:,kd)*eq_2%jac_FD(:,:,kd,0,0,0)
5072 do kd = 1,grid_eq%loc_n_r
5073 do jd = 1,grid_eq%n(2)
5074 ierr =
calc_int(d3sigma_alt(:,jd,kd),ang_par_f(:,jd,kd),&
5075 &sigma_alt(:,jd,kd))
5078 sigma_alt(:,:,kd) = eq_2%sigma(1,1,kd) + sigma_alt(:,:,kd)
5082 call plot_diff_hdf5(d3sigma,d3sigma_alt,&
5083 &
'TEST_diff_D3sigma_through_kappa_g',&
5084 &grid_eq%n,[0,0,grid_eq%i_min-1],&
5085 &descr=
'To test whether -2 p'' J kappa_g = D3sigma',&
5086 &output_message=.true.)
5087 call plot_diff_hdf5(eq_2%sigma,sigma_alt,&
5088 &
'TEST_diff_sigma_through_kappa_g',&
5089 &grid_eq%n,[0,0,grid_eq%i_min-1],descr=
'To test whether &
5090 &int(-2 p'' J kappa_g) = sigma',output_message=.true.)
5096 end function test_sigma_with_kappa_g
5099 integer function test_kappa(grid_eq,eq_1,eq_2)
result(ierr)
5102 character(*),
parameter :: rout_name =
'test_kappa'
5105 type(grid_type),
intent(in) :: grid_eq
5106 type(eq_1_type),
intent(in) :: eq_1
5107 type(eq_2_type),
intent(in),
target :: eq_2
5110 real(dp),
allocatable :: kappa_alt(:,:,:,:)
5111 real(dp),
pointer :: j(:,:,:) => null()
5112 real(dp),
pointer :: d1j(:,:,:) => null()
5113 real(dp),
pointer :: d2j(:,:,:) => null()
5114 real(dp),
pointer :: d3j(:,:,:) => null()
5115 real(dp),
pointer :: g13(:,:,:) => null()
5116 real(dp),
pointer :: g33(:,:,:) => null()
5117 real(dp),
pointer :: d1g33(:,:,:) => null()
5118 real(dp),
pointer :: d2g33(:,:,:) => null()
5119 real(dp),
pointer :: d3g33(:,:,:) => null()
5120 real(dp),
pointer :: h12(:,:,:) => null()
5121 real(dp),
pointer :: h22(:,:,:) => null()
5122 real(dp),
pointer :: h23(:,:,:) => null()
5127 call writo(
'Testing whether kappa agrees with naive implementation')
5131 j => eq_2%jac_FD(:,:,:,0,0,0)
5132 d1j => eq_2%jac_FD(:,:,:,1,0,0)
5133 d2j => eq_2%jac_FD(:,:,:,0,1,0)
5134 d3j => eq_2%jac_FD(:,:,:,0,0,1)
5135 g13 => eq_2%g_FD(:,:,:,
c([1,3],.true.),0,0,0)
5136 g33 => eq_2%g_FD(:,:,:,
c([3,3],.true.),0,0,0)
5137 d1g33 => eq_2%g_FD(:,:,:,
c([3,3],.true.),1,0,0)
5138 d2g33 => eq_2%g_FD(:,:,:,
c([3,3],.true.),0,1,0)
5139 d3g33 => eq_2%g_FD(:,:,:,
c([3,3],.true.),0,0,1)
5140 h12 => eq_2%h_FD(:,:,:,
c([1,2],.true.),0,0,0)
5141 h22 => eq_2%h_FD(:,:,:,
c([2,2],.true.),0,0,0)
5142 h23 => eq_2%h_FD(:,:,:,
c([2,3],.true.),0,0,0)
5145 allocate(kappa_alt(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r,2))
5149 do kd = 1,grid_eq%loc_n_r
5150 kappa_alt(:,:,kd,1) = &
5151 &vac_perm*j(:,:,kd)**2*eq_1%pres_FD(kd,1)/g33(:,:,kd) + &
5152 &1._dp/(2*h22(:,:,kd)) * ( &
5153 &h12(:,:,kd) * ( d1g33(:,:,kd)/g33(:,:,kd) - &
5154 &2*d1j(:,:,kd)/j(:,:,kd) ) + &
5155 &h22(:,:,kd) * ( d2g33(:,:,kd)/g33(:,:,kd) - &
5156 &2*d2j(:,:,kd)/j(:,:,kd) ) + &
5157 &h23(:,:,kd) * ( d3g33(:,:,kd)/g33(:,:,kd) - &
5158 &2*d3j(:,:,kd)/j(:,:,kd) ) )
5162 kappa_alt(:,:,:,2) = (0.5*d1g33/g33 - d1j/j) - &
5163 &g13/g33*(0.5*d3g33/g33 - d3j/j)
5165 call plot_diff_hdf5(eq_2%kappa_n,kappa_alt(:,:,:,1),&
5166 &
'TEST_diff_kappa_n',grid_eq%n,[0,0,grid_eq%i_min-1],&
5167 &descr=
'To test whether kappa_n agrees with naive calculation',&
5168 &output_message=.true.)
5169 call plot_diff_hdf5(eq_2%kappa_g,kappa_alt(:,:,:,2),&
5170 &
'TEST_diff_kappa_g',grid_eq%n,[0,0,grid_eq%i_min-1],&
5171 &descr=
'To test whether kappa_g agrees with naive calculation',&
5172 &output_message=.true.)
5183 nullify(j,d1j,d2j,d3j,g13,g33,d1g33,d2g33,d3g33,h12,h22,h23)
5186 end function test_kappa
5190 integer function test_sigma_vmec(grid_eq,eq_1,eq_2)
result(ierr)
5196 character(*),
parameter :: rout_name =
'test_sigma_VMEC'
5199 type(grid_type),
intent(in) :: grid_eq
5200 type(eq_1_type),
intent(in) :: eq_1
5201 type(eq_2_type),
intent(in),
target :: eq_2
5204 real(
dp),
allocatable :: sigma_alt(:,:,:)
5205 real(
dp),
allocatable :: b_v(:,:,:,:)
5206 real(
dp),
allocatable :: b_alpha(:,:,:)
5207 real(
dp),
pointer :: j(:,:,:) => null()
5208 real(
dp),
pointer :: d1j(:,:,:) => null()
5209 real(
dp),
pointer :: g13(:,:,:) => null()
5210 real(
dp),
pointer :: g23(:,:,:) => null()
5211 real(
dp),
pointer :: d1g23(:,:,:) => null()
5212 real(
dp),
pointer :: g33(:,:,:) => null()
5217 call writo(
'Testing whether sigma agrees with naive implementation')
5221 j => eq_2%jac_FD(:,:,:,0,0,0)
5222 d1j => eq_2%jac_FD(:,:,:,1,0,0)
5223 g13 => eq_2%g_FD(:,:,:,
c([1,3],.true.),0,0,0)
5224 g23 => eq_2%g_FD(:,:,:,
c([2,3],.true.),0,0,0)
5225 d1g23 => eq_2%g_FD(:,:,:,
c([2,3],.true.),1,0,0)
5226 g33 => eq_2%g_FD(:,:,:,
c([3,3],.true.),0,0,0)
5229 allocate(sigma_alt(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r))
5233 allocate(b_v(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r,2:3))
5236 &
b_v_sub_c(:,grid_eq%i_min:grid_eq%i_max,id),&
5237 &
b_v_sub_s(:,grid_eq%i_min:grid_eq%i_max,id),&
5238 &grid_eq%trigon_factors,b_v(:,:,:,id),&
5239 &[.true.,is_asym_v])
5244 allocate(b_alpha(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r))
5247 b_alpha = b_alpha + b_v(:,:,:,kd) * &
5248 &eq_2%T_FE(:,:,:,
c([1,kd],.false.),0,0,0)
5255 do jd = 1,grid_eq%n(2)
5256 do id = 1,grid_eq%n(1)
5257 ierr =
spline(grid_eq%loc_r_F,b_alpha(id,jd,:),&
5258 &grid_eq%loc_r_F,sigma_alt(id,jd,:),&
5265 sigma_alt = (d1g23 - g23*d1j/j)/j - sigma_alt
5269 do kd = 1,grid_eq%loc_n_r
5270 sigma_alt(:,:,kd) = sigma_alt(:,:,kd) / vac_perm - &
5271 &eq_1%pres_FD(kd,1)*j(:,:,kd)*g13(:,:,kd)/g33(:,:,kd)
5274 call plot_diff_hdf5(eq_2%sigma,sigma_alt,
'TEST_diff_sigma',&
5275 &grid_eq%n,[0,0,grid_eq%i_min-1],descr=
'To test whether &
5276 &sigma agrees with naive calculation',output_message=.true.)
5284 nullify(j,d1j,g13,g23,d1g23,g33)
5287 end function test_sigma_vmec
5291 subroutine test_s_hel(grid_eq,eq_2,Rchi,Zchi)
5293 type(grid_type),
intent(in) :: grid_eq
5294 type(eq_2_type),
intent(in) :: eq_2
5295 real(dp),
intent(in) :: rchi(:,:,0:)
5296 real(dp),
intent(in) :: zchi(:,:,0:)
5299 real(dp),
allocatable :: s_alt(:,:,:)
5301 call writo(
'Testing whether shear agrees with alternative &
5302 &implementation using sigma')
5305 allocate(s_alt(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r))
5308 call calc_derived_s_from_sigma_hel(grid_eq,eq_2,rchi,zchi,s_alt)
5309 call plot_diff_hdf5(eq_2%S,s_alt,
'TEST_diff_S_from_sigma',&
5310 &grid_eq%n,[0,0,grid_eq%i_min-1],descr=
'To test whether &
5311 &sigma agrees with naive calculation',output_message=.true.)
5314 call calc_derived_s_direct(eq_2,s_alt)
5315 call plot_diff_hdf5(eq_2%S,s_alt,
'TEST_diff_S_naive',&
5316 &grid_eq%n,[0,0,grid_eq%i_min-1],descr=
'To test whether &
5317 &sigma agrees with naive calculation',output_message=.true.)
5320 end subroutine test_s_hel
5323 integer function plot_diff_for_paper(r, A, B, title)
result(ierr)
5326 character(*),
parameter :: rout_name =
'plot_diff_for_paper'
5329 real(
dp),
intent(in) :: r(:)
5330 real(
dp),
intent(in) :: a(:,:,:)
5331 real(
dp),
intent(in) :: b(:,:,:)
5332 character(len=*),
intent(in) :: title
5338 character(len=max_str_ln) :: plot_title(3)
5345 chckerr(
'Need 1 process')
5352 plot_title(1) = trim(title)//
'_GOOD'
5353 plot_title(2) = trim(title)//
'_BAD'
5354 plot_title(3) = trim(title)//
'_DIFF'
5357 plot_title(kd) = trim(plot_title(kd))//
'_'//&
5361 call print_ex_2d([
''],plot_title(1),transpose(a(:,jd,:)),&
5362 &x=reshape(r*2*
pi/max_flux_f,[n_r,1]),draw=.false.)
5363 call draw_ex([
''],plot_title(1),n_par,1,.false.)
5364 call print_ex_2d([
''],plot_title(2),transpose(b(:,jd,:)),&
5365 &x=reshape(r*2*
pi/max_flux_f,[n_r,1]),draw=.false.)
5366 call draw_ex([
''],plot_title(2),n_par,1,.false.)
5367 call print_ex_2d([
''],plot_title(3),&
5368 &transpose(abs(b(:,jd,:)-a(:,jd,:))/&
5369 &(abs(b(:,jd,:))+abs(a(:,jd,:)))),&
5370 &x=reshape(r*2*
pi/max_flux_f,[n_r,1]),draw=.false.)
5371 call draw_ex([
''],plot_title(3),n_par,1,.false.)
5373 end function plot_diff_for_paper
5409 character(*),
parameter :: rout_name =
'calc_normalization_const'
5412 integer :: nr_overridden_const
5413 character(len=max_str_ln) :: err_msg
5419 br_normalization_provided = [.false.,.false.]
5423 call writo(
'Calculating the normalization constants')
5427 nr_overridden_const = 0
5436 call calc_normalization_const_vmec
5438 call calc_normalization_const_hel
5443 if (nr_overridden_const.gt.0) &
5444 &
call writo(trim(
i2str(nr_overridden_const))//&
5445 &
' constants were overridden by user. Consistency is NOT &
5446 &checked!',warning=.true.)
5453 if (
t_0.lt.0._dp)
then
5455 err_msg =
'Alfven time is negative. Are you sure the &
5456 &equilibrium is consistent?'
5462 call writo(
'Normalization constants calculated')
5465 call writo(
'Normalization not used')
5470 subroutine calc_normalization_const_vmec
5496 call writo(
'Using MISHKA normalization')
5500 if (
r_0.ge.huge(1._dp))
then
5503 nr_overridden_const = nr_overridden_const + 1
5507 if (
b_0.ge.huge(1._dp))
then
5510 nr_overridden_const = nr_overridden_const + 1
5514 if (
pres_0.ge.huge(1._dp))
then
5517 nr_overridden_const = nr_overridden_const + 1
5521 if (
psi_0.ge.huge(1._dp))
then
5524 nr_overridden_const = nr_overridden_const + 1
5528 call writo(
'Using COBRA normalization')
5531 if (
r_0.ge.huge(1._dp))
then
5532 r_0 = 0.5_dp*(rmin_surf+rmax_surf)
5534 nr_overridden_const = nr_overridden_const + 1
5538 if (
pres_0.ge.huge(1._dp))
then
5541 nr_overridden_const = nr_overridden_const + 1
5545 if (
b_0.ge.huge(1._dp))
then
5549 nr_overridden_const = nr_overridden_const + 1
5553 if (
psi_0.ge.huge(1._dp))
then
5557 nr_overridden_const = nr_overridden_const + 1
5561 call writo(
'Exact COBRA normalization is substituted by &
5562 &"pure", modified version',warning=.true.)
5563 call writo(
'This leaves the equations unmodified',&
5565 call writo(
'To translate to exact COBRA, manually multiply &
5566 &PB3D Eigenvalues by',alert=.true.)
5568 call writo(trim(
r2str(beta_v/2)),alert=.true.)
5575 if (
t_0.ge.huge(1._dp))
then
5578 nr_overridden_const = nr_overridden_const + 1
5580 end subroutine calc_normalization_const_vmec
5584 subroutine calc_normalization_const_hel
5586 call writo(
'Using MISHKA normalization')
5588 if (
r_0.ge.huge(1._dp))
then
5591 br_normalization_provided(1) = .true.
5593 if (
b_0.ge.huge(1._dp))
then
5596 br_normalization_provided(2) = .true.
5598 if (
pres_0.ge.huge(1._dp))
then
5601 nr_overridden_const = nr_overridden_const + 1
5603 if (
psi_0.ge.huge(1._dp))
then
5606 nr_overridden_const = nr_overridden_const + 1
5609 end subroutine calc_normalization_const_hel
5613 subroutine print_normalization_const(R_0,rho_0,B_0,pres_0,psi_0,&
5616 real(dp),
intent(in),
optional ::
r_0
5617 real(dp),
intent(in),
optional ::
rho_0
5618 real(dp),
intent(in),
optional ::
b_0
5619 real(dp),
intent(in),
optional ::
pres_0
5620 real(dp),
intent(in),
optional ::
psi_0
5621 real(dp),
intent(in),
optional :: mu_0
5622 real(dp),
intent(in),
optional ::
t_0
5633 if (
present(mu_0))
call writo(
'mu_0 = '//&
5634 &trim(
r2str(mu_0))//
' Tm/A')
5636 end subroutine print_normalization_const
5650 call writo(
'Start normalizing the input variables')
5664 call writo(
'HELENA input is already normalized with MISHKA &
5670 call writo(
'Normalization done')
5698 integer function b_plot(grid_eq,eq_1,eq_2,rich_lvl,plot_fluxes,XYZ) &
5707 character(*),
parameter :: rout_name =
'B_plot'
5710 type(grid_type),
intent(inout) :: grid_eq
5713 integer,
intent(in),
optional :: rich_lvl
5714 logical,
intent(in),
optional :: plot_fluxes
5715 real(
dp),
intent(in),
optional :: xyz(:,:,:,:)
5719 real(
dp),
allocatable :: b_com(:,:,:,:,:)
5720 real(
dp),
allocatable :: b_mag(:,:,:)
5721 real(
dp),
allocatable,
save :: b_flux_tor(:,:), b_flux_pol(:,:)
5722 character(len=10) :: base_name
5723 logical :: plot_fluxes_loc
5729 if (
eq_style.eq.1 .and. .not.
allocated(grid_eq%trigon_factors))
then
5731 chckerr(
'trigonometric factors not allocated')
5735 plot_fluxes_loc = .false.
5736 if (
present(plot_fluxes)) plot_fluxes_loc = plot_fluxes
5739 allocate(b_com(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r,3,2))
5740 allocate(b_mag(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r))
5744 b_com(:,:,:,id,1) = eq_2%g_FD(:,:,:,
c([3,id],.true.),0,0,0)/&
5745 &eq_2%jac_FD(:,:,:,0,0,0)
5747 b_com(:,:,:,3,2) = 1._dp/eq_2%jac_FD(:,:,:,0,0,0)
5754 if (
present(rich_lvl))
then
5755 if (rich_lvl.gt.0) base_name = trim(base_name)//
'_R_'//&
5756 &trim(
i2str(rich_lvl))
5760 if (plot_fluxes_loc)
then
5763 &v_flux_tor=b_flux_tor,v_flux_pol=b_flux_pol,xyz=xyz)
5801 integer function j_plot(grid_eq,eq_1,eq_2,rich_lvl,plot_fluxes,XYZ) &
5815 character(*),
parameter :: rout_name =
'J_plot'
5818 type(grid_type),
intent(inout) :: grid_eq
5821 integer,
intent(in),
optional :: rich_lvl
5822 logical,
intent(in),
optional :: plot_fluxes
5823 real(
dp),
intent(in),
optional :: xyz(:,:,:,:)
5827 real(
dp),
allocatable :: j_com(:,:,:,:,:)
5828 real(
dp),
allocatable :: j_mag(:,:,:)
5829 character(len=10) :: base_name
5830 real(
dp),
allocatable,
save :: j_flux_tor(:,:), j_flux_pol(:,:)
5831 logical :: plot_fluxes_loc
5832 character(len=max_str_ln) :: plot_name
5833 character(len=max_str_ln) :: plot_titles(2)
5835 real(
dp),
allocatable :: j_v_sup_int2(:,:)
5842 if (
eq_style.eq.1 .and. .not.
allocated(grid_eq%trigon_factors))
then
5844 chckerr(
'trigonometric factors not allocated')
5848 plot_fluxes_loc = .false.
5849 if (
present(plot_fluxes)) plot_fluxes_loc = plot_fluxes
5852 allocate(j_com(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r,3,2))
5853 allocate(j_mag(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r))
5858 j_com(:,:,:,1,2) = eq_2%g_FD(:,:,:,
c([3,3],.true.),0,1,0) - &
5859 &eq_2%g_FD(:,:,:,
c([3,3],.true.),0,0,0)*&
5860 &eq_2%jac_FD(:,:,:,0,1,0)/&
5861 &eq_2%jac_FD(:,:,:,0,0,0) - &
5862 &eq_2%g_FD(:,:,:,
c([2,3],.true.),0,0,1) + &
5863 &eq_2%g_FD(:,:,:,
c([2,3],.true.),0,0,0)*&
5864 &eq_2%jac_FD(:,:,:,0,0,1)/&
5865 &eq_2%jac_FD(:,:,:,0,0,0)
5866 j_com(:,:,:,1,2) = j_com(:,:,:,1,2)/eq_2%jac_FD(:,:,:,0,0,0)**2
5867 j_com(:,:,:,3,2) = eq_2%g_FD(:,:,:,
c([2,3],.true.),1,0,0) - &
5868 &eq_2%g_FD(:,:,:,
c([2,3],.true.),0,0,0)*&
5869 &eq_2%jac_FD(:,:,:,1,0,0)/&
5870 &eq_2%jac_FD(:,:,:,0,0,0) - &
5871 &eq_2%g_FD(:,:,:,
c([1,3],.true.),0,1,0) + &
5872 &eq_2%g_FD(:,:,:,
c([1,3],.true.),0,0,0)*&
5873 &eq_2%jac_FD(:,:,:,0,1,0)/&
5874 &eq_2%jac_FD(:,:,:,0,0,0)
5875 j_com(:,:,:,3,2) = j_com(:,:,:,3,2)/eq_2%jac_FD(:,:,:,0,0,0)**2
5878 do kd = 1,grid_eq%loc_n_r
5879 j_com(:,:,kd,1,2) = -eq_1%pres_FD(kd,1)
5880 j_com(:,:,kd,3,2) = eq_2%sigma(:,:,kd)/&
5881 &eq_2%jac_FD(:,:,kd,0,0,0) + eq_1%pres_FD(kd,1)*&
5882 &eq_2%g_FD(:,:,kd,
c([1,3],.true.),0,0,0) / &
5883 &eq_2%g_FD(:,:,kd,
c([3,3],.true.),0,0,0)
5889 j_com(:,:,:,id,1) = &
5890 &j_com(:,:,:,1,2)*eq_2%g_FD(:,:,:,
c([1,id],.true.),0,0,0) + &
5891 &j_com(:,:,:,3,2)*eq_2%g_FD(:,:,:,
c([3,id],.true.),0,0,0)
5899 if (
present(rich_lvl))
then
5900 if (rich_lvl.gt.0) base_name = trim(base_name)//
'_R_'//&
5901 &trim(
i2str(rich_lvl))
5905 if (plot_fluxes_loc)
then
5908 &v_flux_tor=j_flux_tor,v_flux_pol=j_flux_pol,xyz=xyz)
5922 plot_name =
'J_V_sup_int'
5923 plot_titles = [
'J^theta_V',
'J^phi_V ']
5924 call print_ex_2d(plot_titles,plot_name,
j_v_sup_int,&
5926 &[
size(grid_eq%r_F),1]),draw=.false.)
5927 call draw_ex(plot_titles,plot_name,2,1,.false.)
5933 &j_v_sup_int2(:,id))
5936 plot_name =
'J_V_sup_int2'
5937 plot_titles = [
'J^theta_V',
'J^phi_V ']
5938 call print_ex_2d(plot_titles,plot_name,j_v_sup_int2,&
5940 &[
size(grid_eq%r_F),1]),draw=.false.)
5941 call draw_ex(plot_titles,plot_name,2,1,.false.)
5985 character(*),
parameter :: rout_name =
'kappa_plot'
5988 type(grid_type),
intent(inout) :: grid_eq
5991 integer,
intent(in),
optional :: rich_lvl
5992 real(
dp),
intent(in),
optional :: xyz(:,:,:,:)
5995 real(
dp),
allocatable :: k_com(:,:,:,:,:)
5996 real(
dp),
allocatable :: k_mag(:,:,:)
5997 character(len=15) :: base_name
5999 type(grid_type) :: grid_trim
6001 integer :: norm_id(2)
6002 integer :: plot_dim(4)
6003 integer :: plot_offset(4)
6004 real(
dp),
allocatable :: xyz_loc(:,:,:,:,:)
6005 real(
dp),
allocatable :: k_com_inv(:,:,:,:)
6006 logical,
save :: asked_for_testing = .false.
6007 logical,
save :: testing = .false.
6014 if (
eq_style.eq.1 .and. .not.
allocated(grid_eq%trigon_factors))
then
6016 chckerr(
'trigonometric factors not allocated')
6020 allocate(k_com(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r,3,2))
6021 allocate(k_mag(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r))
6030 k_com(:,:,:,1,1) = eq_2%kappa_g
6031 k_com(:,:,:,2,1) = eq_2%kappa_n - &
6033 &eq_2%h_FD(:,:,:,
c([1,2],.true.),0,0,0)/&
6034 &eq_2%h_FD(:,:,:,
c([2,2],.true.),0,0,0)
6035 k_com(:,:,:,3,1) = 0._dp
6042 k_com(:,:,:,1,2) = eq_2%kappa_n * &
6043 &eq_2%h_FD(:,:,:,
c([1,2],.true.),0,0,0) + &
6045 &eq_2%g_FD(:,:,:,
c([3,3],.true.),0,0,0) / &
6046 &(eq_2%jac_FD(:,:,:,0,0,0)**2*&
6047 &eq_2%h_FD(:,:,:,
c([2,2],.true.),0,0,0))
6048 k_com(:,:,:,2,2) = eq_2%kappa_n * &
6049 &eq_2%h_FD(:,:,:,
c([2,2],.true.),0,0,0)
6050 k_com(:,:,:,3,2) = eq_2%kappa_n * &
6051 &eq_2%h_FD(:,:,:,
c([3,2],.true.),0,0,0) - &
6053 &eq_2%g_FD(:,:,:,
c([3,1],.true.),0,0,0) / &
6054 &(eq_2%jac_FD(:,:,:,0,0,0)**2*&
6055 &eq_2%h_FD(:,:,:,
c([2,2],.true.),0,0,0))
6062 if (
present(rich_lvl))
then
6063 if (rich_lvl.gt.0) base_name = trim(base_name)//
'_R_'//&
6064 &trim(
i2str(rich_lvl))
6074 if (.not.asked_for_testing)
then
6075 call writo(
'Do you want to plot the inversion in kappa?')
6078 asked_for_testing = .true.
6082 call writo(
'Every point in the grid is transformed by &
6083 &displacing it in the direction of the curvature,')
6084 call writo(
'by a distance equal to the inverse of the &
6085 &curvature. I.e. the point is displaced to the center &
6090 ierr =
trim_grid(grid_eq,grid_trim,norm_id)
6094 plot_dim = [grid_trim%n,3]
6095 plot_offset = [0,0,grid_trim%i_min-1,0]
6105 allocate(k_com_inv(grid_trim%n(1),grid_trim%n(2),&
6106 &grid_trim%loc_n_r,3))
6108 k_com_inv(:,:,:,id) = &
6109 &k_com(:,:,norm_id(1):norm_id(2),id,1)/&
6110 &(k_mag(:,:,norm_id(1):norm_id(2))**2)
6116 &grid_trim%zeta_E,grid_trim%trigon_factors)
6121 allocate(xyz_loc(grid_trim%n(1),grid_trim%n(2),&
6122 &grid_trim%loc_n_r,3,3))
6124 &xyz_loc(:,:,:,1,2),xyz_loc(:,:,:,1,3))
6128 call plot_hdf5([
'cen_of_curv'],
'TEST_cen_of_curv_vec',&
6129 &k_com_inv,tot_dim=plot_dim,loc_offset=plot_offset,&
6130 &x=xyz_loc(:,:,:,:,1),y=xyz_loc(:,:,:,:,2),&
6131 &z=xyz_loc(:,:,:,:,3),col=4,cont_plot=
eq_job_nr.gt.1,&
6132 &descr=
'center of curvature')
6136 xyz_loc(:,:,:,1,id) = xyz_loc(:,:,:,1,id) + &
6137 &k_com_inv(:,:,:,id)
6139 xyz_loc(:,:,:,2,:) = xyz_loc(:,:,:,1,:)
6140 xyz_loc(:,:,:,3,:) = xyz_loc(:,:,:,3,:)
6143 call plot_hdf5([
'cen_of_curv_inv'],
'TEST_cen_of_curv_inv_vec',&
6144 &-k_com_inv,tot_dim=plot_dim,loc_offset=plot_offset,&
6145 &x=xyz_loc(:,:,:,:,1),y=xyz_loc(:,:,:,:,2),&
6146 &z=xyz_loc(:,:,:,:,3),col=4,cont_plot=
eq_job_nr.gt.1,&
6147 &descr=
'center of curvature')
6150 call grid_trim%dealloc()
6184 character(*),
parameter :: rout_name =
'delta_r_plot'
6187 type(grid_type),
intent(inout) :: grid_eq
6190 real(
dp),
intent(in) :: xyz(:,:,:,:)
6191 integer,
intent(in),
optional :: rich_lvl
6195 integer :: norm_id(2)
6196 integer :: r_lim_2d(2)
6198 real(
dp),
allocatable :: xyz_loc(:,:,:,:)
6199 real(
dp),
allocatable :: b_com(:,:,:,:,:)
6200 real(
dp),
allocatable :: delta_b_tor(:,:,:)
6201 real(
dp),
allocatable :: prop_b_tor_tot(:,:)
6202 real(
dp),
allocatable :: prop_b_tor(:,:,:)
6203 real(
dp),
allocatable :: var_tot_loc(:)
6204 real(
dp),
allocatable :: x_2d(:,:)
6205 real(
dp),
allocatable :: delta_r(:,:,:)
6206 real(
dp),
allocatable :: theta_geo(:,:,:)
6207 real(
dp),
allocatable :: r_geo(:,:,:)
6208 real(
dp),
allocatable :: prop_b_tor_plot(:,:)
6209 real(
dp),
allocatable :: delta_r_f(:,:)
6210 logical :: new_file_found
6211 character(len=25) :: base_name
6212 character(len=25) :: plot_name
6213 character(len=max_str_ln) :: prop_b_tor_file_name
6214 type(grid_type) :: grid_trim
6220 if (
eq_style.eq.1 .and. .not.
allocated(grid_eq%trigon_factors))
then
6222 chckerr(
'trigonometric factors not allocated')
6225 call writo(
'Calculate normal displacement')
6229 ierr =
trim_grid(grid_eq,grid_trim,norm_id)
6234 allocate(xyz_loc(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r,4))
6235 ierr =
calc_xyz_grid(grid_eq,grid_eq,xyz_loc(:,:,:,1),xyz_loc(:,:,:,2),&
6236 &xyz_loc(:,:,:,3),r=xyz_loc(:,:,:,4))
6240 allocate(theta_geo(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r))
6241 allocate(r_geo(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r))
6242 theta_geo = atan2(xyz_loc(:,:,:,3)-
rz_0(2),&
6243 &xyz_loc(:,:,:,4)-
rz_0(1))
6244 where (theta_geo.lt.0._dp) theta_geo = theta_geo + 2*
pi
6245 r_geo = sqrt((xyz_loc(:,:,:,3)-
rz_0(2))**2 + &
6246 &(xyz_loc(:,:,:,4)-
rz_0(1))**2)
6249 call plot_hdf5(
'theta_geo',
'theta_geo',&
6250 &theta_geo(:,:,norm_id(1):norm_id(2)),&
6251 &tot_dim=[grid_trim%n(1),grid_trim%n(2),grid_trim%n(3)],&
6252 &loc_offset=[0,0,grid_trim%i_min-1],&
6253 &x=xyz(:,:,norm_id(1):norm_id(2),1),&
6254 &y=xyz(:,:,norm_id(1):norm_id(2),2),&
6255 &z=xyz(:,:,norm_id(1):norm_id(2),3),cont_plot=
eq_job_nr.gt.1,&
6256 &descr=
'geometrical poloidal angle')
6257 call plot_hdf5(
'r_geo',
'r_geo',&
6258 &r_geo(:,:,norm_id(1):norm_id(2)),&
6259 &tot_dim=[grid_trim%n(1),grid_trim%n(2),grid_trim%n(3)],&
6260 &loc_offset=[0,0,grid_trim%i_min-1],&
6261 &x=xyz(:,:,norm_id(1):norm_id(2),1),&
6262 &y=xyz(:,:,norm_id(1):norm_id(2),2),&
6263 &z=xyz(:,:,norm_id(1):norm_id(2),3),cont_plot=
eq_job_nr.gt.1,&
6264 &descr=
'geometrical radius')
6268 &absolute=.true.,r=grid_eq%loc_r_F)
6270 allocate(delta_r(grid_eq%n(1),1,grid_eq%loc_n_r))
6271 delta_r(:,1,:) = r_geo(:,2,:)*0.5_dp
6275 base_name =
'delta_r'
6276 plot_name =
'delta_r'
6277 if (
present(rich_lvl))
then
6278 if (rich_lvl.gt.0) base_name = trim(base_name)//
'_R_'//&
6279 &trim(
i2str(rich_lvl))
6283 call plot_hdf5(trim(plot_name),trim(base_name),&
6284 &delta_r(:,:,norm_id(1):norm_id(2)),&
6285 &tot_dim=[grid_trim%n(1),1,grid_trim%n(3)],&
6286 &loc_offset=[0,0,grid_trim%i_min-1],&
6287 &x=xyz(:,2:2,norm_id(1):norm_id(2),1),&
6288 &y=xyz(:,2:2,norm_id(1):norm_id(2),2),&
6289 &z=xyz(:,2:2,norm_id(1):norm_id(2),3),&
6291 &descr=
'plasma position displacement')
6296 call print_ex_2d(trim(plot_name),trim(base_name),&
6297 &delta_r(1:grid_eq%n(1)-1,1,grid_eq%loc_n_r),&
6298 &x=theta_geo(1:grid_eq%n(1)-1,2,grid_eq%loc_n_r),&
6300 call draw_ex([trim(plot_name)],trim(base_name),1,1,.false.)
6303 ierr =
nufft(theta_geo(1:grid_eq%n(1)-1,2,grid_eq%loc_n_r),&
6304 &delta_r(1:grid_eq%n(1)-1,1,grid_eq%loc_n_r),delta_r_f)
6306 base_name = trim(base_name)//
'_F'
6307 call print_ex_2d([
'delta_r cos',
'delta_r sin'],trim(base_name),&
6308 &delta_r_f,draw=.false.)
6309 call draw_ex([
'delta_r cos',
'delta_r sin'],trim(base_name),2,1,&
6313 new_file_found = .false.
6314 prop_b_tor_file_name = base_name
6316 do while (.not.new_file_found)
6317 open(
prop_b_tor_i,file=trim(prop_b_tor_file_name)//
'_'//&
6318 &trim(
i2str(kd))//
'.dat',iostat=ierr,status=
'old')
6322 prop_b_tor_file_name = trim(prop_b_tor_file_name)//&
6323 &
'_'//trim(
i2str(kd))//
'.dat'
6324 new_file_found = .true.
6326 &iostat=ierr,status=
'new')
6327 chckerr(
'Failed to open file')
6334 &
'N',
'M',
'delta_c',
'delta_s'
6335 do id = 1,min(
size(delta_r_f,1),max_m+1)
6337 &18, id-1, delta_r_f(id,:)
6346 call writo(
'Calculate magnetic field ripple')
6350 allocate(b_com(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r,3,2))
6353 b_com(:,:,:,id,1) = eq_2%g_FD(:,:,:,
c([3,id],.true.),0,0,0)/&
6354 &eq_2%jac_FD(:,:,:,0,0,0)
6356 b_com(:,:,:,3,2) = 1._dp/eq_2%jac_FD(:,:,:,0,0,0)
6366 allocate(delta_b_tor(grid_eq%n(1),1,grid_eq%loc_n_r))
6367 delta_b_tor(:,1,:) = 0.5_dp*sum(b_com(:,2,:,2,:),3)
6371 base_name =
'delta_B_tor'
6372 plot_name =
'delta_B_tor'
6373 if (
present(rich_lvl))
then
6374 if (rich_lvl.gt.0) base_name = trim(base_name)//
'_R_'//&
6375 &trim(
i2str(rich_lvl))
6379 call plot_hdf5(trim(plot_name),trim(base_name),&
6380 &delta_b_tor(:,:,norm_id(1):norm_id(2)),&
6381 &tot_dim=[grid_trim%n(1),1,grid_trim%n(3)],&
6382 &loc_offset=[0,0,grid_trim%i_min-1],&
6383 &x=xyz(:,2:2,norm_id(1):norm_id(2),1),&
6384 &y=xyz(:,2:2,norm_id(1):norm_id(2),2),&
6385 &z=xyz(:,2:2,norm_id(1):norm_id(2),3),cont_plot=
eq_job_nr.gt.1,&
6386 &descr=
'plasma position displacement')
6390 call writo(
'Calculate proportionality factor')
6393 allocate(prop_b_tor(grid_eq%n(1),1,grid_eq%loc_n_r))
6395 where (abs(delta_b_tor).gt.
tol_zero) prop_b_tor = delta_r / delta_b_tor
6398 base_name =
'prop_B_tor'
6399 plot_name = trim(base_name)
6400 if (
present(rich_lvl))
then
6401 if (rich_lvl.gt.0) base_name = trim(base_name)//
'_R_'//&
6402 &trim(
i2str(rich_lvl))
6406 call plot_hdf5(plot_name,trim(base_name),&
6407 &prop_b_tor(:,:,norm_id(1):norm_id(2)),&
6408 &tot_dim=[grid_trim%n(1),1,grid_trim%n(3)],&
6409 &loc_offset=[0,0,grid_trim%i_min-1],&
6410 &x=xyz(:,2:2,norm_id(1):norm_id(2),1),&
6411 &y=xyz(:,2:2,norm_id(1):norm_id(2),2),&
6412 &z=xyz(:,2:2,norm_id(1):norm_id(2),3),&
6414 &descr=
'delta_r divided by delta_B_tor')
6417 r_lim_2d = [1,grid_trim%n(3)]
6419 &max(r_lim_2d(1),r_lim_2d(2)-127+1)
6421 allocate(x_2d(grid_trim%n(1),grid_trim%n(3)))
6423 allocate(prop_b_tor_tot(grid_trim%n(1),grid_trim%n(3)))
6424 do id = 1,grid_trim%n(1)
6426 ierr =
get_ser_var(prop_b_tor(id,1,norm_id(1):norm_id(2)),&
6427 &var_tot_loc,scatter=.true.)
6429 prop_b_tor_tot(id,:) = max(min(var_tot_loc,2._dp),-2._dp)
6430 deallocate(var_tot_loc)
6433 ierr =
get_ser_var(theta_geo(id,2,norm_id(1):norm_id(2)),&
6436 if (
rank.eq.0) x_2d(id,:) = var_tot_loc/
pi
6437 deallocate(var_tot_loc)
6440 call print_ex_2d([trim(plot_name)],trim(base_name),&
6441 &prop_b_tor_tot(:,r_lim_2d(1):r_lim_2d(2)),&
6442 &x=x_2d(:,r_lim_2d(1):r_lim_2d(2)),draw=.false.)
6443 call draw_ex([trim(plot_name)],trim(base_name),&
6444 &r_lim_2d(2)-r_lim_2d(1)+1,1,.false.)
6449 call writo(
'Output to file as function of poloidal flux angle')
6455 new_file_found = .false.
6456 prop_b_tor_file_name =
'prop_B_tor'
6458 do while (.not.new_file_found)
6459 open(
prop_b_tor_i,file=trim(prop_b_tor_file_name)//
'_'//&
6460 &trim(
i2str(kd))//
'.dat',iostat=ierr,status=
'old')
6464 prop_b_tor_file_name = trim(prop_b_tor_file_name)//&
6465 &
'_'//trim(
i2str(kd))//
'.dat'
6466 new_file_found = .true.
6468 &iostat=ierr,status=
'new')
6469 chckerr(
'Failed to open file')
6474 call writo(
'Save toroidal field proportionality factor in &
6475 &file "'//trim(prop_b_tor_file_name)//
'"',persistent=.true.)
6480 &theta_geo(1:grid_trim%n(1)-1,1,grid_trim%loc_n_r),&
6481 &prop_b_tor_tot(1:grid_trim%n(1)-1,grid_trim%n(3))],&
6482 &[grid_trim%n(1)-1,2]),prop_b_tor_plot,0)
6487 &
'pol. angle [ ]',
'prop. factor [ ]'
6488 do id = 1,grid_trim%n(1)-1
6490 &prop_b_tor_plot(id,:)
6500 call grid_trim%dealloc()
6565 &n_par_X_base,range_name)
result(ierr)
6574 character(*),
parameter :: rout_name =
'divide_eq_jobs'
6577 integer,
intent(in) ::
n_par_x
6578 integer,
intent(in) :: arr_size(2)
6579 integer,
intent(inout) :: n_div
6580 integer,
intent(in),
optional :: n_div_max
6581 integer,
intent(in),
optional :: n_par_x_base
6582 character(len=*),
intent(in),
optional :: range_name
6585 real(
dp) :: mem_size(2)
6586 integer :: n_par_x_base_loc
6587 integer :: max_mem_req
6588 integer :: n_div_max_loc
6589 integer :: n_par_range
6590 character(len=max_str_ln) :: range_message
6591 character(len=max_str_ln) :: err_msg
6592 character(len=max_str_ln) :: range_name_loc
6598 call writo(
'Dividing the equilibrium jobs')
6610 n_par_x_base_loc = 0
6611 if (
present(n_par_x_base)) n_par_x_base_loc = n_par_x_base
6614 range_name_loc =
'parallel points'
6615 if (
present(range_name)) range_name_loc = trim(range_name)
6619 mem_size = huge(1._dp)*0.49_dp
6621 n_div_max_loc = (
n_par_x-1)/fund_n_par
6623 n_div_max_loc =
n_par_x/fund_n_par
6625 if (
present(n_div_max)) n_div_max_loc = n_div_max
6626 n_div_max_loc = max(1,n_div_max_loc)
6631 n_par_range = ceiling(
n_par_x*1._dp/n_div + n_par_x_base_loc)
6637 if (n_div.gt.n_div_max_loc)
then
6639 err_msg =
'The memory limit is too low, need more than '//&
6640 &trim(
i2str(max_mem_req))//
'MB'
6643 max_mem_req = ceiling(sum(mem_size))
6645 if (n_div.gt.1)
then
6647 &trim(range_name_loc)//
' are split into '//&
6648 &trim(
i2str(n_div))//
' and '//trim(
i2str(n_div))//&
6649 &
' collective jobs are done serially'
6652 &trim(range_name_loc)//
' can be done without splitting them'
6654 call writo(range_message)
6655 call writo(
'The total memory for all processes together is estimated &
6656 &to be about '//trim(
i2str(ceiling(sum(mem_size))))//
'MB')
6659 &MB, user specified)')
6666 call writo(
'In the perturbation phase, the equilibrium variables &
6667 &are not being operated on:')
6669 call writo(
'This translates to a scale factor 1/'//&
6671 call writo(
'Therefore, the memory left for the perturbation phase &
6678 call writo(
'Equilibrium jobs divided')
6704 character(*),
parameter :: rout_name =
'calc_eq_jobs_lims'
6707 integer,
intent(in) ::
n_par_x
6708 integer,
intent(in) :: n_div
6713 integer,
allocatable :: n_par(:)
6714 character(len=max_str_ln) :: err_msg
6722 allocate(n_par(n_div))
6755 if (n_div.gt.1)
then
6757 &fund_n_par * max(1,nint(&
6766 err_msg =
'Limits don''t match, try with more memory or lower &
6778 integer function test_t_ef(grid_eq,eq_1,eq_2)
result(ierr)
6784 character(*),
parameter :: rout_name =
'test_T_EF'
6787 type(grid_type),
intent(in) :: grid_eq
6788 type(eq_1_type),
intent(in) :: eq_1
6789 type(eq_2_type),
intent(in) :: eq_2
6793 integer :: norm_id(2)
6794 real(
dp),
allocatable :: res(:,:,:,:)
6795 character(len=max_str_ln) :: file_name
6796 character(len=max_str_ln) :: description
6797 integer :: tot_dim(3), loc_offset(3)
6798 type(grid_type) :: grid_trim
6804 call writo(
'Going to test whether T_EF complies with the theory')
6808 ierr =
trim_grid(grid_eq,grid_trim,norm_id)
6812 tot_dim = [grid_trim%n(1),grid_trim%n(2),grid_trim%n(3)]
6813 loc_offset = [0,0,grid_trim%i_min-1]
6816 allocate(res(grid_trim%n(1),grid_trim%n(2),grid_trim%loc_n_r,9))
6827 do kd = norm_id(1),norm_id(2)
6828 res(:,:,kd-norm_id(1)+1,1) = &
6829 &grid_eq%theta_F(:,:,kd)*eq_1%q_saf_E(kd,1) + &
6830 &eq_2%L_E(:,:,kd,1,0,0)*eq_1%q_saf_E(kd,0)
6833 do kd = norm_id(1),norm_id(2)
6834 res(:,:,kd-norm_id(1)+1,2) = &
6835 &(1._dp + eq_2%L_E(:,:,kd,0,1,0))*eq_1%q_saf_E(kd,0)
6838 do kd = norm_id(1),norm_id(2)
6839 res(:,:,kd-norm_id(1)+1,3) = -1._dp + &
6840 &eq_2%L_E(:,:,kd,0,0,1)*eq_1%q_saf_E(kd,0)
6843 do kd = norm_id(1),norm_id(2)
6844 res(:,:,kd-norm_id(1)+1,4) = eq_1%flux_p_E(kd,1)/(2*
pi)
6847 res(:,:,:,7) = eq_2%L_E(:,:,norm_id(1):norm_id(2),1,0,0)
6849 res(:,:,:,8) = 1._dp + &
6850 &eq_2%L_E(:,:,norm_id(1):norm_id(2),0,1,0)
6852 res(:,:,:,9) = eq_2%L_E(:,:,norm_id(1):norm_id(2),0,0,1)
6855 do kd = norm_id(1),norm_id(2)
6856 res(:,:,kd-norm_id(1)+1,1) = - eq_2%L_E(:,:,kd,1,0,0) &
6857 &+ grid_eq%zeta_E(:,:,kd)*eq_1%rot_t_E(kd,1)
6861 &- (1._dp + eq_2%L_E(:,:,norm_id(1):norm_id(2),0,1,0))
6863 do kd = norm_id(1),norm_id(2)
6864 res(:,:,kd-norm_id(1)+1,3) = &
6865 &eq_1%rot_t_E(kd,0) - eq_2%L_E(:,:,kd,0,0,1)
6868 do kd = norm_id(1),norm_id(2)
6869 res(:,:,kd-norm_id(1)+1,4) = &
6870 &- eq_1%flux_t_E(kd,1)/(2*
pi)
6873 res(:,:,:,9) = -1._dp
6878 do kd = norm_id(1),norm_id(2)
6879 res(:,:,kd-norm_id(1)+1,1) = &
6880 &- eq_1%q_saf_E(kd,1)*grid_eq%theta_E(:,:,kd)
6883 do kd = norm_id(1),norm_id(2)
6884 res(:,:,kd-norm_id(1)+1,2) = - eq_1%q_saf_E(kd,0)
6887 res(:,:,:,3) = 1._dp
6889 do kd = norm_id(1),norm_id(2)
6890 res(:,:,kd-norm_id(1)+1,4) = eq_1%flux_p_E(kd,1)/(2*
pi)
6893 res(:,:,:,8) = 1._dp
6896 do kd = norm_id(1),norm_id(2)
6897 res(:,:,kd-norm_id(1)+1,1) = &
6898 &eq_1%rot_t_E(kd,1)*grid_eq%zeta_E(:,:,kd)
6901 res(:,:,:,2) = -1._dp
6903 do kd = norm_id(1),norm_id(2)
6904 res(:,:,kd-norm_id(1)+1,3) = eq_1%rot_t_E(kd,0)
6907 do kd = norm_id(1),norm_id(2)
6908 res(:,:,kd-norm_id(1)+1,4) = eq_1%flux_t_E(kd,1)/(2*
pi)
6911 res(:,:,:,9) = 1._dp
6919 call writo(
'Testing T_EF('//trim(
i2str(kd))//
','//&
6920 &trim(
i2str(id))//
')')
6924 file_name =
'TEST_T_EF_'//trim(
i2str(kd))//
'_'//trim(
i2str(id))
6925 description =
'Testing calculated with given value for T_EF('//&
6930 &eq_2%T_EF(:,:,norm_id(1):norm_id(2),
c([kd,id],.false.),&
6931 &0,0,0),file_name,tot_dim,loc_offset,description,&
6932 &output_message=.true.)
6939 call grid_trim%dealloc()
6943 call writo(
'Test complete')
6944 end function test_t_ef
6952 integer function test_d12h_h(grid_eq,eq)
result(ierr)
6958 character(*),
parameter :: rout_name =
'test_D12h_H'
6961 type(grid_type),
intent(in) :: grid_eq
6962 type(eq_2_type),
intent(in) :: eq
6965 integer :: norm_id(2)
6966 integer :: id, jd, kd, ld
6969 real(
dp) :: bcs_val(2,3)
6970 real(
dp),
allocatable :: res(:,:,:,:)
6971 character(len=max_str_ln) :: file_name
6972 character(len=max_str_ln) :: description
6973 integer :: tot_dim(3), loc_offset(3)
6974 type(grid_type) :: grid_trim
6980 call writo(
'Going to test whether D1 D2 h_H is calculated correctly')
6984 ierr =
trim_grid(grid_eq,grid_trim,norm_id)
6988 allocate(res(grid_trim%n(1),grid_trim%n(2),grid_trim%loc_n_r,6))
6991 tot_dim = [grid_trim%n(1),grid_trim%n(2),grid_trim%n(3)]
6992 loc_offset = [0,0,grid_trim%i_min-1]
7005 bc_ld = [1,2,0,1,0,1]
7009 do kd = norm_id(1),norm_id(2)
7010 do jd = 1,grid_trim%n(2)
7011 if (bc_ld(ld).eq.0) cycle
7013 ierr =
spline(grid_eq%theta_E(:,jd,kd),&
7014 &eq%h_E(:,jd,kd,ld,0,1,0),grid_eq%theta_E(:,jd,kd),&
7016 &deriv=1,bcs=bcs(:,bc_ld(ld)),&
7017 &bcs_val=bcs_val(:,bc_ld(ld)))
7027 call writo(
'Testing h_H('//trim(
i2str(kd))//
','//&
7028 &trim(
i2str(id))//
')')
7032 file_name =
'TEST_D12h_H_'//trim(
i2str(kd))//
'_'//&
7034 description =
'Testing calculated with given value for D12h_H('&
7035 &//trim(
i2str(kd))//
','//trim(
i2str(id))//
')'
7038 call plot_diff_hdf5(res(:,:,:,
c([kd,id],.true.)),&
7039 &eq%h_E(:,:,norm_id(1):norm_id(2),&
7040 &
c([kd,id],.true.),1,1,0),file_name,tot_dim,loc_offset,&
7041 &description,output_message=.true.)
7048 call grid_trim%dealloc()
7052 call writo(
'Test complete')
7053 end function test_d12h_h
7063 integer function test_jac_f(grid_eq,eq_1,eq_2)
result(ierr)
7069 character(*),
parameter :: rout_name =
'test_jac_F'
7072 type(grid_type),
intent(in) :: grid_eq
7073 type(eq_1_type),
intent(in),
target :: eq_1
7074 type(eq_2_type),
intent(in) :: eq_2
7077 integer :: norm_id(2)
7078 real(
dp),
allocatable :: res(:,:,:)
7080 character(len=max_str_ln) :: file_name
7081 character(len=max_str_ln) :: description
7082 integer :: tot_dim(3), loc_offset(3)
7083 type(grid_type) :: grid_trim
7084 real(
dp),
pointer :: dflux(:) => null()
7091 call writo(
'Going to test the calculation of the Jacobian in Flux &
7096 ierr =
trim_grid(grid_eq,grid_trim,norm_id)
7100 allocate(res(grid_trim%n(1),grid_trim%n(2),grid_trim%loc_n_r))
7103 tot_dim = [grid_trim%n(1),grid_trim%n(2),grid_trim%n(3)]
7104 loc_offset = [0,0,grid_trim%i_min-1]
7109 ierr =
calc_det(res,eq_2%g_F(:,:,norm_id(1):norm_id(2),:,0,0,0),3)
7113 file_name =
'TEST_jac_F_1'
7114 description =
'Testing whether the Jacobian in Flux coordinates is &
7115 &consistent with determinant of metric matrix'
7118 call plot_diff_hdf5(res(:,:,:),&
7119 &eq_2%jac_F(:,:,norm_id(1):norm_id(2),0,0,0)**2,&
7120 &file_name,tot_dim,loc_offset,description,output_message=.true.)
7126 dflux => eq_1%flux_p_E(:,1)
7129 dflux => eq_1%flux_t_E(:,1)
7139 do kd = norm_id(1),norm_id(2)
7140 res(:,:,kd-norm_id(1)+1) = pmone * 2*
pi * &
7141 &eq_2%R_E(:,:,kd,0,0,0)*&
7142 &(eq_2%R_E(:,:,kd,1,0,0)*eq_2%Z_E(:,:,kd,0,1,0) - &
7143 &eq_2%Z_E(:,:,kd,1,0,0)*eq_2%R_E(:,:,kd,0,1,0)) / &
7144 &(dflux(kd)*(1+eq_2%L_E(:,:,kd,0,1,0)))
7147 do kd = norm_id(1),norm_id(2)
7148 do jd = 1,grid_trim%n(2)
7149 res(:,jd,kd-norm_id(1)+1) = eq_1%q_saf_E(kd,0)/&
7150 &(
h_h_33(:,kd+grid_eq%i_min-1)*&
7151 &
rbphi_h(kd+grid_eq%i_min-1,0))
7157 file_name =
'TEST_jac_F_2'
7158 description =
'Testing whether the Jacobian in Flux coordinates is &
7159 &consistent with explicit formula'
7162 call plot_diff_hdf5(res(:,:,:),&
7163 &eq_2%jac_F(:,:,norm_id(1):norm_id(2),0,0,0),file_name,tot_dim,&
7164 &loc_offset,description,output_message=.true.)
7168 call grid_trim%dealloc()
7172 call writo(
'Test complete')
7173 end function test_jac_f
7180 integer function test_g_v(grid_eq,eq)
result(ierr)
7184 character(*),
parameter :: rout_name =
'test_g_V'
7187 type(grid_type),
intent(in) :: grid_eq
7188 type(eq_2_type),
intent(in) :: eq
7191 integer :: norm_id(2)
7193 real(dp),
allocatable :: res(:,:,:,:)
7194 character(len=max_str_ln) :: file_name
7195 character(len=max_str_ln) :: description
7196 integer :: tot_dim(3), loc_offset(3)
7197 type(grid_type) :: grid_trim
7203 call writo(
'Going to test whether g_V is calculated correctly')
7207 ierr =
trim_grid(grid_eq,grid_trim,norm_id)
7211 allocate(res(grid_trim%n(1),grid_trim%n(2),grid_trim%loc_n_r,6))
7214 tot_dim = [grid_trim%n(1),grid_trim%n(2),grid_trim%n(3)]
7215 loc_offset = [0,0,grid_trim%i_min-1]
7218 res(:,:,:,1) = eq%R_E(:,:,norm_id(1):norm_id(2),1,0,0)**2 + &
7219 &eq%Z_E(:,:,norm_id(1):norm_id(2),1,0,0)**2
7221 res(:,:,:,2) = eq%R_E(:,:,norm_id(1):norm_id(2),1,0,0)*&
7222 &eq%R_E(:,:,norm_id(1):norm_id(2),0,1,0) + &
7223 &eq%Z_E(:,:,norm_id(1):norm_id(2),1,0,0)*&
7224 &eq%Z_E(:,:,norm_id(1):norm_id(2),0,1,0)
7226 res(:,:,:,3) = eq%R_E(:,:,norm_id(1):norm_id(2),1,0,0)*&
7227 &eq%R_E(:,:,norm_id(1):norm_id(2),0,0,1) + &
7228 &eq%Z_E(:,:,norm_id(1):norm_id(2),1,0,0)*&
7229 &eq%Z_E(:,:,norm_id(1):norm_id(2),0,0,1)
7231 res(:,:,:,4) = eq%R_E(:,:,norm_id(1):norm_id(2),0,1,0)**2 + &
7232 &eq%Z_E(:,:,norm_id(1):norm_id(2),0,1,0)**2
7234 res(:,:,:,5) = eq%R_E(:,:,norm_id(1):norm_id(2),0,1,0)*&
7235 &eq%R_E(:,:,norm_id(1):norm_id(2),0,0,1) + &
7236 &eq%Z_E(:,:,norm_id(1):norm_id(2),0,1,0)*&
7237 &eq%Z_E(:,:,norm_id(1):norm_id(2),0,0,1)
7239 res(:,:,:,6) = eq%R_E(:,:,norm_id(1):norm_id(2),0,0,1)**2 + &
7240 &eq%Z_E(:,:,norm_id(1):norm_id(2),0,0,1)**2 + &
7241 &eq%R_E(:,:,norm_id(1):norm_id(2),0,0,0)**2
7247 call writo(
'Testing g_V('//trim(
i2str(kd))//
','//&
7248 &trim(
i2str(id))//
')')
7252 file_name =
'TEST_g_V_'//trim(
i2str(kd))//
'_'//trim(
i2str(id))
7253 description =
'Testing calculated with given value for g_V('//&
7257 call plot_diff_hdf5(res(:,:,:,
c([kd,id],.true.)),&
7258 &eq%g_E(:,:,norm_id(1):norm_id(2),
c([kd,id],.true.),&
7259 &0,0,0),file_name,tot_dim,loc_offset,description,&
7260 &output_message=.true.)
7267 call grid_trim%dealloc()
7271 call writo(
'Test complete')
7272 end function test_g_v
7279 integer function test_jac_v(grid_eq,eq)
result(ierr)
7283 character(*),
parameter :: rout_name =
'test_jac_V'
7286 type(grid_type),
intent(in) :: grid_eq
7287 type(eq_2_type),
intent(in) :: eq
7290 integer :: norm_id(2)
7291 real(dp),
allocatable :: res(:,:,:)
7292 character(len=max_str_ln) :: file_name
7293 character(len=max_str_ln) :: description
7294 integer :: tot_dim(3), loc_offset(3)
7295 type(grid_type) :: grid_trim
7301 call writo(
'Going to test the calculation of the Jacobian in the VMEC &
7306 ierr =
trim_grid(grid_eq,grid_trim,norm_id)
7310 allocate(res(grid_trim%n(1),grid_trim%n(2),grid_trim%loc_n_r))
7313 tot_dim = [grid_trim%n(1),grid_trim%n(2),grid_trim%n(3)]
7314 loc_offset = [0,0,grid_trim%i_min-1]
7317 ierr =
calc_det(res,eq%g_E(:,:,norm_id(1):norm_id(2),:,0,0,0),3)
7321 file_name =
'TEST_jac_V'
7322 description =
'Testing whether the Jacobian in VMEC coordinates is &
7323 &consistent with determinant of metric matrix'
7326 call plot_diff_hdf5(-sqrt(res),&
7327 &eq%jac_E(:,:,norm_id(1):norm_id(2),0,0,0),&
7328 &file_name,tot_dim,loc_offset,description,output_message=.true.)
7331 call grid_trim%dealloc()
7335 call writo(
'Test complete')
7336 end function test_jac_v
7343 integer function test_b_f(grid_eq,eq_1,eq_2)
result(ierr)
7350 character(*),
parameter :: rout_name =
'test_B_F'
7353 type(grid_type),
intent(in) :: grid_eq
7354 type(eq_1_type),
intent(in) :: eq_1
7355 type(eq_2_type),
intent(in) :: eq_2
7358 integer :: norm_id(2)
7359 integer :: norm_id_f(2)
7361 real(
dp),
allocatable :: res(:,:,:,:)
7362 real(
dp),
allocatable :: res2(:,:,:,:)
7363 character(len=max_str_ln) :: file_name
7364 character(len=max_str_ln) :: description
7365 integer :: tot_dim(3), loc_offset(3)
7366 type(grid_type) :: grid_trim
7372 call writo(
'Going to test whether the magnetic field is calculated &
7377 ierr =
trim_grid(grid_eq,grid_trim,norm_id)
7381 norm_id_f = grid_eq%i_min+norm_id-1
7384 allocate(res(grid_trim%n(1),grid_trim%n(2),grid_trim%loc_n_r,4))
7385 allocate(res2(grid_trim%n(1),grid_trim%n(2),grid_trim%loc_n_r,4))
7388 tot_dim = [grid_trim%n(1),grid_trim%n(2),grid_trim%n(3)]
7389 loc_offset = [0,0,grid_trim%i_min-1]
7399 &
b_v_sub_c(:,norm_id_f(1):norm_id_f(2),id),&
7400 &
b_v_sub_s(:,norm_id_f(1):norm_id_f(2),id),&
7401 &grid_eq%trigon_factors(:,:,:,norm_id(1):&
7402 &norm_id(2),:),res(:,:,:,id))
7406 &
b_v_s(:,norm_id_f(1):norm_id_f(2)),&
7407 &grid_eq%trigon_factors(:,:,:,norm_id(1):norm_id(2),:),&
7408 &res(:,:,:,4),[.true.,is_asym_v])
7412 &eq_2%g_E(:,:,norm_id(1):norm_id(2),
c([1,2],.true.),0,0,0)
7414 &eq_2%g_E(:,:,norm_id(1):norm_id(2),
c([2,2],.true.),0,0,0)
7415 do kd = norm_id(1),norm_id(2)
7416 res(:,:,kd-norm_id(1)+1,3) = eq_1%q_saf_E(kd,0)*&
7417 &eq_2%g_E(:,:,kd,
c([3,3],.true.),0,0,0)
7418 res(:,:,kd-norm_id(1)+1,4) = &
7419 &sqrt(eq_2%g_E(:,:,kd,
c([2,2],.true.),0,0,0) + &
7420 &eq_1%q_saf_E(kd,0)**2*&
7421 &eq_2%g_E(:,:,kd,
c([3,3],.true.),0,0,0))
7424 do kd = norm_id(1),norm_id(2)
7425 res(:,:,kd-norm_id(1)+1,id) = &
7426 &res(:,:,kd-norm_id(1)+1,id)/&
7427 &eq_2%jac_E(:,:,kd,0,0,0)
7436 res2(:,:,:,id) = res2(:,:,:,id) + &
7437 &res(:,:,:,kd) * eq_2%T_FE(:,:,norm_id(1):&
7438 &norm_id(2),
c([id,kd],.false.),0,0,0)
7441 res2(:,:,:,4) = res(:,:,:,4)
7449 call writo(
'Testing B_F')
7455 file_name =
'TEST_B_F_'//trim(
i2str(id))
7456 description =
'Testing calculated with given value for B_F_'//&
7459 file_name =
'TEST_B_F'
7460 description =
'Testing calculated with given value for B_F'
7465 call plot_diff_hdf5(res2(:,:,:,id),&
7466 &eq_2%g_F(:,:,norm_id(1):norm_id(2),
c([3,id],.true.),0,0,0)&
7467 &/eq_2%jac_F(:,:,norm_id(1):norm_id(2),0,0,0),file_name,&
7468 &tot_dim,loc_offset,description,output_message=.true.)
7470 call plot_diff_hdf5(res2(:,:,:,id),sqrt(&
7471 &eq_2%g_F(:,:,norm_id(1):norm_id(2),
c([3,3],.true.),0,0,0))&
7472 &/eq_2%jac_F(:,:,norm_id(1):norm_id(2),0,0,0),file_name,&
7473 &tot_dim,loc_offset,description,output_message=.true.)
7480 call grid_trim%dealloc()
7484 call writo(
'Test complete')
7485 end function test_b_f
7501 integer function test_p(grid_eq,eq_1,eq_2)
result(ierr)
7508 character(*),
parameter :: rout_name =
'test_p'
7511 type(grid_type),
intent(in) :: grid_eq
7516 integer :: norm_id(2)
7517 integer :: norm_id_f(2)
7518 real(
dp),
allocatable :: res(:,:,:,:)
7519 integer :: id, jd, kd
7520 character(len=max_str_ln) :: file_name
7521 character(len=max_str_ln) :: description
7522 integer :: tot_dim(3), loc_offset(3)
7523 type(grid_type) :: grid_trim
7529 call writo(
'Going to test the consistency of the equilibrium variables &
7530 &with the given pressure')
7534 ierr =
trim_grid(grid_eq,grid_trim,norm_id)
7538 allocate(res(grid_trim%n(1),grid_trim%n(2),grid_trim%loc_n_r,2))
7541 call writo(
'Checking if mu_0 D2 p = 1/J (D3 B_2 - D2_B3)')
7545 tot_dim = [grid_trim%n(1),grid_trim%n(2),grid_trim%n(3)]
7546 loc_offset = [0,0,grid_trim%i_min-1]
7550 &(eq_2%g_FD(:,:,norm_id(1):norm_id(2),
c([2,3],.true.),0,0,1) - &
7551 &eq_2%g_FD(:,:,norm_id(1):norm_id(2),
c([3,3],.true.),0,1,0))/&
7552 &eq_2%jac_FD(:,:,norm_id(1):norm_id(2),0,0,0) - &
7553 &(eq_2%g_FD(:,:,norm_id(1):norm_id(2),
c([2,3],.true.),0,0,0)*&
7554 &eq_2%jac_FD(:,:,norm_id(1):norm_id(2),0,0,1) - &
7555 &eq_2%g_FD(:,:,norm_id(1):norm_id(2),
c([3,3],.true.),0,0,0)*&
7556 &eq_2%jac_FD(:,:,norm_id(1):norm_id(2),0,1,0))/ &
7557 &(eq_2%jac_FD(:,:,norm_id(1):norm_id(2),0,0,0)**2)
7558 res(:,:,:,1) = res(:,:,:,1)/eq_2%jac_FD(:,:,norm_id(1):norm_id(2),0,0,0)
7567 do kd = norm_id(1),norm_id(2)
7568 res(:,:,kd-norm_id(1)+1,2) =
vac_perm*eq_1%pres_FD(kd,1)
7572 file_name =
'TEST_D2p'
7573 description =
'Testing whether mu_0 D2 p = 1/J (D3 B_2 - D2_B3)'
7576 call plot_diff_hdf5(res(:,:,:,1),res(:,:,:,2),file_name,tot_dim,&
7577 &loc_offset,description,output_message=.true.)
7582 call writo(
'Checking if mu_0 J D1p = 0 => D3 B_1 = D2 B_3')
7587 &eq_2%g_FD(:,:,norm_id(1):norm_id(2),
c([1,3],.true.),0,0,1)/&
7588 &eq_2%jac_FD(:,:,norm_id(1):norm_id(2),0,0,0) - &
7589 &eq_2%g_FD(:,:,norm_id(1):norm_id(2),
c([1,3],.true.),0,0,0)*&
7590 &eq_2%jac_FD(:,:,norm_id(1):norm_id(2),0,0,1) / &
7591 &eq_2%jac_FD(:,:,norm_id(1):norm_id(2),0,0,0)**2
7595 &eq_2%g_FD(:,:,norm_id(1):norm_id(2),
c([3,3],.true.),1,0,0)/&
7596 &eq_2%jac_FD(:,:,norm_id(1):norm_id(2),0,0,0) - &
7597 &eq_2%g_FD(:,:,norm_id(1):norm_id(2),
c([3,3],.true.),0,0,0)*&
7598 &eq_2%jac_FD(:,:,norm_id(1):norm_id(2),1,0,0) / &
7599 &eq_2%jac_FD(:,:,norm_id(1):norm_id(2),0,0,0)**2
7602 file_name =
'TEST_D1p'
7603 description =
'Testing whether D3 B_1 = D1 B_3'
7606 call plot_diff_hdf5(res(:,:,:,1),res(:,:,:,2),file_name,tot_dim,&
7607 &loc_offset,description,output_message=.true.)
7614 norm_id_f = grid_eq%i_min+norm_id-1
7617 call writo(
'Checking if D2 B_3 |F = D1 (qF+qh_11/F) |H')
7622 &(eq_2%g_FD(:,:,norm_id(1):norm_id(2),
c([3,3],.true.),0,1,0)/&
7623 &eq_2%jac_FD(:,:,norm_id(1):norm_id(2),0,0,0) - &
7624 &eq_2%g_FD(:,:,norm_id(1):norm_id(2),
c([3,3],.true.),0,0,0)*&
7625 &eq_2%jac_FD(:,:,norm_id(1):norm_id(2),0,1,0)/ &
7626 &(eq_2%jac_FD(:,:,norm_id(1):norm_id(2),0,0,0)**2))
7627 do jd = 1,grid_trim%n(2)
7628 do id = 1,grid_trim%n(1)
7629 ierr =
spline(grid_trim%loc_r_E,&
7630 &eq_1%q_saf_E(norm_id(1):norm_id(2),0)*&
7631 &
rbphi_h(norm_id_f(1):norm_id_f(2),0)+&
7632 &eq_1%q_saf_E(norm_id(1):norm_id(2),0)*&
7633 &
h_h_11(id,norm_id_f(1):norm_id_f(2))/&
7634 &
rbphi_h(norm_id_f(1):norm_id_f(2),0),&
7635 &grid_trim%loc_r_E,res(id,jd,:,2),&
7642 file_name =
'TEST_D2B_3'
7643 description =
'Testing whether D2 B_3 |F = D1 (qF+qh_11/F) |H'
7646 call plot_diff_hdf5(res(:,:,:,1),res(:,:,:,2),file_name,tot_dim,&
7647 &loc_offset,description,output_message=.true.)
7652 call writo(
'Checking if D3 B_2 |F = -q/F D2 h_11 + F q'' |H')
7657 &eq_2%g_FD(:,:,norm_id(1):norm_id(2),
c([2,3],.true.),0,0,1)/&
7658 &eq_2%jac_FD(:,:,norm_id(1):norm_id(2),0,0,0) - &
7659 &eq_2%g_FD(:,:,norm_id(1):norm_id(2),
c([2,3],.true.),0,0,0)*&
7660 &eq_2%jac_FD(:,:,norm_id(1):norm_id(2),0,0,1)/ &
7661 &(eq_2%jac_FD(:,:,norm_id(1):norm_id(2),0,0,0)**2)
7662 do kd = norm_id(1),norm_id(2)
7663 do jd = 1,grid_trim%n(2)
7664 ierr =
spline(grid_eq%theta_E(:,jd,kd),&
7665 &
h_h_12(:,kd+grid_eq%i_min-1),grid_eq%theta_E(:,jd,kd),&
7666 &res(:,jd,kd-norm_id(1)+1,2),&
7670 res(:,:,kd-norm_id(1)+1,2) = &
7671 &-eq_1%q_saf_E(kd,0)/
rbphi_h(kd+grid_eq%i_min-1,0) * &
7672 &res(:,:,kd-norm_id(1)+1,2) + &
7673 &
rbphi_h(kd+grid_eq%i_min-1,0) * &
7678 file_name =
'TEST_D3B_2'
7679 description =
'Testing whether D3 B_2 |F = -D2 (qh_12/F) + f q'' |H'
7682 call plot_diff_hdf5(res(:,:,:,1),res(:,:,:,2),file_name,tot_dim,&
7683 &loc_offset,description,output_message=.true.)
7689 call grid_trim%dealloc()
7693 call writo(
'Test complete')
Calculate the equilibrium quantities on a grid determined by straight field lines.
Calculate the lower metric elements in the C(ylindrical) coordinate system.
Calculate the metric coefficients in the F(lux) coordinate system.
Calculate the lower metric coefficients in the equilibrium H(ELENA) coordinate system.
Calculate the lower metric coefficients in the equilibrium V(MEC) coordinate system.
Calculate , the jacobian in Flux coordinates.
Calculate , the jacobian in HELENA coordinates.
Calculate , the jacobian in V(MEC) coordinates.
Calculate , & and derivatives in VMEC coordinates.
Calculate , the transformation matrix between H(ELENA) and F(lux) coordinate systems.
Calculate , the transformation matrix between C(ylindrical) and V(mec) coordinate system.
Calculate , the transformation matrix between V(MEC) and F(lux) coordinate systems.
Print equilibrium quantities to an output file:
Redistribute the equilibrium variables, but only the Flux variables are saved.
Transforms derivatives of the equilibrium quantities in E coordinates to derivatives in the F coordin...
Calculate from and where and , according to weyens3d.
Calculates the toroidal difference for a magnitude calculated on three toroidal points: two extremiti...
Gather parallel variable in serial version on group master.
Add to an array (3) the product of arrays (1) and (2).
Sorting with the bubble sort routine.
Calculate determinant of a matrix.
Integrates a function using the trapezoidal rule.
Order a periodic function to include and an overlap.
Wrapper to the pspline library, making it easier to use for 1-D applications where speed is not the m...
Prints variables vars with names var_names in an HDF5 file with name c file_name and accompanying XDM...
Print 2-D output on a file.
Print 3-D output on a file.
Inverse Fourier transformation, from VMEC.
Operations on the equilibrium variables.
logical, public debug_create_vmec_input
plot debug information for create_vmec_input()
integer function, public calc_derived_q(grid_eq, eq_1, eq_2)
Calculates derived equilibrium quantities system.
subroutine, public normalize_input()
Normalize input quantities.
logical, public debug_calc_derived_q
plot debug information for calc_derived_q()
integer function, public delta_r_plot(grid_eq, eq_1, eq_2, xyz, rich_lvl)
Plots HALF of the change in the position vectors for 2 different toroidal positions,...
logical, public debug_j_plot
plot debug information for j_plot()
integer function, public kappa_plot(grid_eq, eq_1, eq_2, rich_lvl, xyz)
Plots the curvature.
integer function, public divide_eq_jobs(n_par_x, arr_size, n_div, n_div_max, n_par_x_base, range_name)
Divides the equilibrium jobs.
integer function, public calc_eq_jobs_lims(n_par_x, n_div)
Calculate eq_jobs_lims.
integer function, public calc_normalization_const()
Sets up normalization constants.
integer function, public j_plot(grid_eq, eq_1, eq_2, rich_lvl, plot_fluxes, xyz)
Plots the current.
integer function, public b_plot(grid_eq, eq_1, eq_2, rich_lvl, plot_fluxes, xyz)
Plots the magnetic fields.
integer function, public flux_q_plot(grid_eq, eq)
Plots the flux quantities in the solution grid.
Numerical utilities related to equilibrium variables.
integer function, public calc_memory_eq(arr_size, n_par, mem_size)
Calculate memory in MB necessary for variables in equilibrium job.
integer function, public calc_g(g_a, t_ba, g_b, deriv, max_deriv)
Calculate the metric coefficients in a coordinate system B ! using the.
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
Numerical utilities related to files.
integer function, public count_lines(file_i)
Count non-comment lines in a file.
integer function, public skip_comment(file_i, file_name)
Skips comment when reading a file.
Numerical utilities related to the grids and different coordinate systems.
integer function, public extend_grid_f(grid_in, grid_ext, grid_eq, n_theta_plot, n_zeta_plot, lim_theta_plot, lim_zeta_plot)
Extend a grid angularly.
integer function, public trim_grid(grid_in, grid_out, norm_id)
Trim a grid, removing any overlap between the different regions.
integer function, public calc_xyz_grid(grid_eq, grid_xyz, x, y, z, l, r)
Calculates , and on a grid grid_XYZ, determined through its E(quilibrium) coordinates.
integer function, public calc_vec_comp(grid, grid_eq, eq_1, eq_2, v_com, norm_disc_prec, v_mag, base_name, max_transf, v_flux_tor, v_flux_pol, xyz, compare_tor_pos)
Calculates contra- and covariant components of a vector in multiple coordinate systems.
integer function, public nufft(x, f, f_f, plot_name)
calculates the cosine and sine mode numbers of a function defined on a non-regular grid.
Variables pertaining to the different grids used.
real(dp), dimension(:), allocatable, public alpha
field line label alpha
integer, public n_r_eq
nr. of normal points in equilibrium 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 test_metrics_h()
Checks whether the metric elements provided by HELENA are consistent with a direct calculation using ...
integer function, public test_harm_cont_h()
Investaige harmonic content of the HELENA variables.
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), public bmtog_h
B_geo/B_mag.
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), public rmtog_h
R_geo/R_mag.
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 utilities related to MPI.
integer function, public redistribute_var(var, dis_var, lims, lims_dis)
Redistribute variables according to new limits.
integer function, public check_deriv(deriv, max_deriv, sr_name)
checks whether the derivatives requested for a certain subroutine are valid
recursive integer function, public gcd(u, v)
Returns least denominator using the GCD.
integer function, public c(ij, sym, n, lim_n)
Convert 2-D coordinates (i,j) to the storage convention used in matrices.
subroutine, public shift_f(al, bl, cl, a, b, c)
Calculate multiplication through shifting of fourier modes A and B into C.
integer function, dimension(:,:), allocatable, public derivs(order, dims)
Returns derivatives of certain order.
integer, dimension(:,:), allocatable, public m
1-D array indices of metric indices
Numerical variables used by most other modules.
integer, parameter, public dp
double precision
logical, public ltest
whether or not to call the testing routines
real(dp), parameter, public pi
real(dp), parameter, public mem_scale_fac
empirical scale factor of memory to calculate eq compared to just storing it
integer, public n_procs
nr. of MPI processes
character(len=max_str_ln), public eq_name
name of equilibrium file from VMEC or HELENA
integer, parameter, public max_str_ln
maximum length of strings
integer, public prog_style
program style (1: PB3D, 2: PB3D_POST)
real(dp), parameter, public mu_0_original
permeability of free space
logical, public use_normalization
whether to use normalization or not
integer, dimension(:,:), allocatable, public eq_jobs_lims
data about eq jobs: [ , ] for all jobs
integer, public rho_style
style for equilibrium density profile
integer, parameter, public max_deriv
highest derivatives for metric factors in Flux coords.
integer, public eq_style
either 1 (VMEC) or 2 (HELENA)
integer, parameter, public hel_pert_i
file number of HELENA equilibrium perturbation file
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
logical, public export_hel
export HELENA
integer, public ex_plot_style
external plot style (1: GNUPlot, 2: Bokeh for 2D, Mayavi for 3D)
real(dp), dimension(2), public rz_0
origin of geometrical poloidal coordinate
integer, public rank
MPI rank.
integer, public rich_restart_lvl
starting Richardson level (0: none [default])
integer, parameter, public hel_export_i
file number of output of HELENA equilibrium export file
integer, public norm_style
style for normalization
integer, parameter, public prop_b_tor_i
file number of proportionality factor file
logical, public use_pol_flux_e
whether poloidal flux is used in E coords.
integer, public eq_job_nr
nr. of eq job
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]
integer, public magn_int_style
style for magnetic integrals (1: trapezoidal, 2: Simpson 3/8)
real(dp), public max_x_mem
maximum memory for perturbation calculations for all processes [MB]
logical, public no_plots
no plots made
Operations concerning giving output, on the screen as well as in output files.
subroutine, public plot_diff_hdf5(a, b, file_name, tot_dim, loc_offset, descr, output_message)
Takes two input vectors and plots these as well as the relative and absolute difference in a HDF5 fil...
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.
Variables concerning Richardson extrapolation.
integer, public rich_lvl
current level of Richardson extrapolation
integer, public n_par_x
nr. of parallel points in field-aligned grid
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.
subroutine, public normalize_vmec
Normalizes VMEC input.
Numerical utilities related to the output of VMEC.
integer function, public calc_trigon_factors(theta, zeta, trigon_factors)
Calculate the trigonometric cosine and sine factors.
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).
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), public b_0_v
the magnitude of B at the magnetic axis,
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
Numerical utilities related to perturbation operations.
integer function, public calc_memory_x(ord, arr_size, n_mod, mem_size)
Calculate memory in MB necessary for X variables.
Variables pertaining to the perturbation quantities.
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.
1D equivalent of multidimensional variables, used for internal HDF5 storage.