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,:)
403 grid_eq%r_F = -
flux_t_v(:,0)/(2*pi)
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
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)
478 select case (eq_style)
482 call writo(
'Calculate R, Z, λ, ...')
483 if (.not.
allocated(grid_eq%trigon_factors))
then
485 &grid_eq%trigon_factors)
488 do id = 0,max_deriv+1
499 select case (eq_style)
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?')
537 call writo(
'Test calculation of jac_V?')
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?')
638 call writo(
'Exporting HELENA equilibrium for VMEC porting')
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?')
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?')
736 call writo(
'Test calculation of B_F?')
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
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.)
925 if (use_normalization)
then
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))
948 &
flux_p_h(:,0)/(2*pi),ffp,ord=norm_disc_prec_eq,deriv=1)
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),&
983 if (use_normalization) i_tor_int = i_tor_int *
r_0*
b_0/mu_0_original
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)//
'"')
1137 open(prop_b_tor_i,file=path_prefix(1:kd*3)//&
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
1167 read(prop_b_tor_i,*,iostat=ierr) prop_b_tor(id,:)
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),&
1188 &theta_b,prop_b_tor_interp,ord=norm_disc_prec_eq,&
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)//
'"')
1245 open(hel_pert_i,file=path_prefix(1:kd*3)//&
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)
1296 read(hel_pert_i,*) n_pert_map(1)
1297 allocate(pert_map_r(n_pert_map(1)))
1298 do kd = 1,n_pert_map(1)
1299 read(hel_pert_i,*) pert_map_r(kd)
1302 &file_name=hel_pert_file_name)
1304 read(hel_pert_i,*) n_pert_map(2)
1305 allocate(pert_map_z(n_pert_map(2)))
1306 do kd = 1,n_pert_map(2)
1307 read(hel_pert_i,*) pert_map_z(kd)
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)
1315 read(hel_pert_i,*) pert_map(kd,:)
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),&
1933 &ord=norm_disc_prec_eq)
1936 if (use_normalization)
then
1938 pres_vj(:,1) = pres_vj(:,1) /
psi_0
1940 ierr =
spline(psi_t,
rbphi_h(:,0),s_vj,f_vj,ord=norm_disc_prec_eq)
1942 if (use_normalization) f_vj = f_vj *
r_0*
b_0
1943 ierr =
spline(psi_t,ffp,s_vj,ffp_vj,ord=norm_disc_prec_eq)
1945 if (use_normalization)
then
1948 ierr =
spline(psi_t,
flux_p_h(:,0),s_vj,flux_j,ord=norm_disc_prec_eq)
1950 if (use_normalization) flux_j = flux_j*
psi_0
1951 ierr =
spline(psi_t,eq_1%q_saf_E(:,0),s_vj,q_saf_vj,&
1952 &ord=norm_disc_prec_eq)
1954 ierr =
spline(psi_t,-eq_1%rot_t_E(:,0),s_vj,rot_t_v,&
1955 &ord=norm_disc_prec_eq)
1957 ierr =
spline(psi_t,i_tor*norm_transf(:,1),s_vj,i_tor_v,&
1958 &ord=norm_disc_prec_eq)
1960 ierr =
spline(psi_t,i_tor*norm_transf(:,2),s_vj,i_tor_j,&
1961 &ord=norm_disc_prec_eq)
1963 if (use_normalization) i_tor_v = i_tor_v *
r_0*
b_0/mu_0_original
1964 if (use_normalization) i_tor_j = i_tor_j *
r_0*
b_0/mu_0_original
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 -----"
1972 write(hel_export_i,
"(A)")
"&INDATA"
1973 write(hel_export_i,
"(A)")
"MGRID_FILE = 'NONE',"
1974 write(hel_export_i,
"(A)")
"PRECON_TYPE = 'GMRES'"
1975 write(hel_export_i,
"(A)")
"PREC2D_THRESHOLD = 1.E-9"
1976 write(hel_export_i,
"(A)")
"DELT = 1.00E+00,"
1977 write(hel_export_i,
"(A)")
"NS_ARRAY = 19, 39, 79, 159, 319"
1978 write(hel_export_i,
"(A)")
"LRFP = F"
1979 write(hel_export_i,
"(A,L1)")
"LASYM = ", .not.stel_sym
1980 write(hel_export_i,
"(A)")
"LFREEB = F"
1981 write(hel_export_i,
"(A)")
"NTOR = "//trim(
i2str(nr_n-1))
1982 write(hel_export_i,
"(A)")
"MPOL = "//trim(
i2str(max_n_b_output))
1983 write(hel_export_i,
"(A)")
"TCON0 = 1"
1984 write(hel_export_i,
"(A)")
"FTOL_ARRAY = 1.E-6, 1.E-6, 1.E-8, 1.E-10, &
1986 write(hel_export_i,
"(A)")
"NITER = 100000,"
1987 write(hel_export_i,
"(A)")
"NSTEP = 200,"
1988 write(hel_export_i,
"(A)")
"NFP = "//trim(
i2str(nfp))
1989 if (use_normalization)
then
1990 write(hel_export_i,
"(A)")
"PHIEDGE = "//&
1991 &trim(
r2str(-eq_1%flux_t_E(grid_eq%n(3),0)*
psi_0))
1993 write(hel_export_i,
"(A)")
"PHIEDGE = "//&
1994 &trim(
r2str(-eq_1%flux_t_E(grid_eq%n(3),0)))
1996 write(hel_export_i,
"(A)")
"CURTOR = "//&
1997 &trim(
r2str(i_tor_int(
size(i_tor_int))))
1999 write(hel_export_i,
"(A)")
""
2000 write(hel_export_i,
"(A)")
""
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"
2005 write(hel_export_i,
"(A)")
""
2006 write(hel_export_i,
"(A)",advance=
"no")
"AM_AUX_S ="
2007 do kd = 1,
size(s_vj)
2008 write(hel_export_i,
"(A1,ES23.16)")
" ", s_vj(kd)
2010 write(hel_export_i,
"(A)")
""
2011 write(hel_export_i,
"(A)",advance=
"no")
"AM_AUX_F ="
2012 do kd = 1,
size(pres_vj,1)
2013 write(hel_export_i,
"(A1,ES23.16)")
" ", pres_vj(kd,0)
2016 write(hel_export_i,
"(A)")
""
2017 write(hel_export_i,
"(A)")
""
2019 write(hel_export_i,
"(A)") &
2020 &
"!----- Current/Iota Parameters -----"
2021 write(hel_export_i,
"(A)")
"NCURR = "//trim(
i2str(ncurr))
2022 write(hel_export_i,
"(A)")
""
2025 write(hel_export_i,
"(A)")
"PIOTA_TYPE = 'Cubic_spline'"
2026 write(hel_export_i,
"(A)")
""
2027 write(hel_export_i,
"(A)",advance=
"no")
"AI_AUX_S ="
2028 do kd = 1,
size(s_vj)
2029 write(hel_export_i,
"(A1,ES23.16)")
" ", s_vj(kd)
2031 write(hel_export_i,
"(A)")
""
2032 write(hel_export_i,
"(A)",advance=
"no")
"AI_AUX_F ="
2033 do kd = 1,
size(rot_t_v)
2034 write(hel_export_i,
"(A1,ES23.16)")
" ", rot_t_v(kd)
2036 write(hel_export_i,
"(A)")
""
2039 write(hel_export_i,
"(A)")
"PCURR_TYPE = 'cubic_spline_Ip'"
2040 write(hel_export_i,
"(A)")
""
2041 write(hel_export_i,
"(A)",advance=
"no")
"AC_AUX_S ="
2042 do kd = 1,
size(s_vj)
2043 write(hel_export_i,
"(A1,ES23.16)")
" ", s_vj(kd)
2045 write(hel_export_i,
"(A)")
""
2046 write(hel_export_i,
"(A)",advance=
"no")
"AC_AUX_F ="
2047 do kd = 1,
size(rot_t_v)
2048 write(hel_export_i,
"(A1,ES23.16)")
" ", i_tor_v(kd)
2051 write(hel_export_i,
"(A)")
""
2052 write(hel_export_i,
"(A)")
""
2054 write(hel_export_i,
"(A)") &
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
2062 ierr = print_mode_numbers(hel_export_i,&
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)
2069 write(hel_export_i,
"(A)")
"&END"
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')
2090 if (use_normalization)
then
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)
2099 write(hel_export_i,*)
''
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)
2108 write(hel_export_i,*)
''
2109 write(hel_export_i,*)
'! R [m] and Z [m] of boundary at psi [ ]'
2111 write(hel_export_i,
"('R_boundary(',I4,') = ',ES23.16,&
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)
2126 write(hel_export_i,
"(' ',2(ES23.16,' '))") &
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)
2137 write(hel_export_i,
"(' ',2(ES23.16,' '))") &
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
2274 integer function print_output_eq_1(grid_eq,eq,data_name)
result(ierr)
2281 character(*),
parameter :: rout_name =
'print_output_eq_1'
2286 character(len=*),
intent(in) :: data_name
2289 integer :: norm_id(2)
2290 type(
var_1d_type),
allocatable,
target :: eq_1d(:)
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'
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(:)
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'
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.)
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'
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'
2776 integer,
intent(in) :: deriv(3)
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'
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'
2833 integer,
intent(in) :: deriv(:)
2840 ierr = check_deriv(deriv,
max_deriv,
'calc_g_C')
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'
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'
2884 integer,
intent(in) :: 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'
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'
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(:,:)
2948 ierr = check_deriv(deriv,max_deriv,
'calc_g_H')
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'
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'
3097 integer,
intent(in) :: 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'
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'
3142 integer,
intent(in) :: deriv(:)
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'
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'
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)
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'
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'
3310 integer,
intent(in) :: deriv(:)
3317 ierr = check_deriv(deriv,
max_deriv,
'calc_J_F')
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'
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'
3357 integer,
intent(in) :: deriv(:)
3364 ierr = check_deriv(deriv,
max_deriv,
'calc_T_VC')
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'
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'
3423 integer,
intent(in) :: deriv(:)
3427 real(dp),
allocatable :: theta_s(:,:,:,:,:,:)
3428 real(dp),
allocatable :: zeta_s(:,:,:,:,:,:)
3437 ierr = check_deriv(deriv,max_deriv,
'calc_T_VF')
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'
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'
3616 integer,
intent(in) :: deriv(:)
3620 real(dp),
allocatable :: theta_s(:,:,:,:,:,:)
3621 real(dp),
allocatable :: zeta_s(:,:,:,:,:,:)
3629 ierr = check_deriv(deriv,max_deriv,
'calc_T_HF')
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'
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
3809 integer function flux_q_plot(grid_eq,eq)
result(ierr)
3815 character(*),
parameter :: rout_name =
'flux_q_plot'
3822 integer :: norm_id(2)
3824 integer :: n_vars = 5
3825 character(len=max_str_ln),
allocatable :: plot_titles(:)
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)
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)
3946 if (use_normalization)
then
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
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)
4027 &y_plot_2d(:,3),x_plot_2d(:,3),draw=.false.)
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'
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))
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'
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)
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)
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
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)
4796 subroutine calc_derived_kappa_from_e(grid_eq,eq_1,eq_2,b_n,b_g,&
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)
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)
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))
4934 if (use_pol_flux_f)
then
4935 ang_par_f => grid_eq%theta_F
4937 ang_par_f => grid_eq%zeta_F
4943 select case (eq_style)
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)
4964 select case (eq_style)
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
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)
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'
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')
5042 if (use_pol_flux_f)
then
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)
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.)
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'
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)
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.)
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'
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)
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)
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)
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)
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.)
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 ']
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.)
5969 integer function kappa_plot(grid_eq,eq_1,eq_2,rich_lvl,XYZ) &
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
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))
6058 if (use_normalization) k_com = k_com /
r_0
6062 if (
present(rich_lvl))
then
6063 if (rich_lvl.gt.0) base_name = trim(base_name)//
'_R_'//&
6064 &trim(
i2str(rich_lvl))
6068 ierr = calc_vec_comp(grid_eq,grid_eq,eq_1,eq_2,k_com,&
6069 &norm_disc_prec_eq,v_mag=k_mag,base_name=base_name,xyz=xyz)
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)
6114 if (eq_style.eq.1)
then
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()
6171 integer function delta_r_plot(grid_eq,eq_1,eq_2,XYZ,rich_lvl) &
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
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')
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()
6564 integer function divide_eq_jobs(n_par_X,arr_size,n_div,n_div_max,&
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
6646 range_message =
'The '//trim(
i2str(n_par_x))//
' '//&
6647 &trim(range_name_loc)//
' are split into '//&
6648 &trim(
i2str(n_div))//
' and '//trim(
i2str(n_div))//&
6649 &
' collective jobs are done serially'
6651 range_message =
'The '//trim(
i2str(n_par_x))//
' '//&
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))
6723 n_par = n_par_x/n_div
6724 n_par(1:mod(n_par_x,n_div)) = n_par(1:mod(n_par_x,n_div)) + 1
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'
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)
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')
6952 integer function test_d12h_h(grid_eq,eq)
result(ierr)
6958 character(*),
parameter :: rout_name =
'test_D12h_H'
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)
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))//
')'
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')
7063 integer function test_jac_f(grid_eq,eq_1,eq_2)
result(ierr)
7069 character(*),
parameter :: rout_name =
'test_jac_F'
7073 type(
eq_1_type),
intent(in),
target :: eq_1
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)
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'
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'
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')
7180 integer function test_g_v(grid_eq,eq)
result(ierr)
7184 character(*),
parameter :: rout_name =
'test_g_V'
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)
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('//&
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')
7279 integer function test_jac_v(grid_eq,eq)
result(ierr)
7283 character(*),
parameter :: rout_name =
'test_jac_V'
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)
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'
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')
7343 integer function test_b_f(grid_eq,eq_1,eq_2)
result(ierr)
7350 character(*),
parameter :: rout_name =
'test_B_F'
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)
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'
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.)
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')
7501 integer function test_p(grid_eq,eq_1,eq_2)
result(ierr)
7508 character(*),
parameter :: rout_name =
'test_p'
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)
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')