PB3D  [2.45]
Ideal linear high-n MHD stability in 3-D
grid_utilities.f90
Go to the documentation of this file.
1 !------------------------------------------------------------------------------!
3 !------------------------------------------------------------------------------!
4 module grid_utilities
5 #include <PB3D_macros.h>
6 #include <wrappers.h>
7  use str_utilities
8  use output_ops
9  use messages
10  use num_vars, only: dp, pi, max_str_ln, iu
11  use grid_vars, only: grid_type
12 
13  implicit none
14  private
15  public coord_f2e, coord_e2f, calc_xyz_grid, calc_eqd_grid, extend_grid_f, &
16  &calc_int_vol, copy_grid, trim_grid, untrim_grid, calc_n_par_x_rich, &
17  &calc_vec_comp, nufft, find_compr_range, calc_arc_angle, calc_tor_diff
18 #if ldebug
19  public debug_calc_int_vol, debug_calc_vec_comp
20 #endif
21 
22  ! global variables
23 #if ldebug
24 
25  logical :: debug_calc_int_vol = .false.
27  logical :: debug_calc_vec_comp = .false.
28 #endif
29 
30  ! interfaces
31 
47  interface coord_f2e
49  module procedure coord_f2e_r
51  module procedure coord_f2e_rtz
52  end interface
53 
66  interface coord_e2f
68  module procedure coord_e2f_r
70  module procedure coord_e2f_rtz
71  end interface
72 
75  interface calc_eqd_grid
77  module procedure calc_eqd_grid_1d
79  module procedure calc_eqd_grid_3d
80  end interface
81 
109  interface calc_tor_diff
111  module procedure calc_tor_diff_0d
113  module procedure calc_tor_diff_2d
114  end interface
115 
116 contains
118  integer function coord_f2e_rtz(grid_eq,r_F,theta_F,zeta_F,r_E,&
119  &theta_E,zeta_E,r_F_array,r_E_array,ord) result(ierr)
120  use num_vars, only: tol_zero, eq_style
121  use vmec_utilities, only: fourier2real
122  use vmec_vars, only: l_v_c, l_v_s, is_asym_v
123 
124  character(*), parameter :: rout_name = 'coord_F2E_rtz'
125 
126  ! input / output
127  type(grid_type), intent(in) :: grid_eq
128  real(dp), intent(in) :: r_f(:)
129  real(dp), intent(in) :: theta_f(:,:,:)
130  real(dp), intent(in) :: zeta_f(:,:,:)
131  real(dp), intent(inout) :: r_e(:)
132  real(dp), intent(inout) :: theta_e(:,:,:)
133  real(dp), intent(inout) :: zeta_e(:,:,:)
134  real(dp), intent(in), optional, target :: r_f_array(:)
135  real(dp), intent(in), optional, target :: r_e_array(:)
136  integer, intent(in), optional :: ord
137 
138  ! local variables (also used in child functions)
139  character(len=max_str_ln) :: err_msg ! error message
140  integer :: dims(3) ! dimensions of the grid
141  real(dp), allocatable :: l_v_c_loc(:,:) ! local version of L_V_c
142  real(dp), allocatable :: l_v_s_loc(:,:) ! local version of L_V_s
143  integer :: ord_loc ! local ord
144 
145  ! initialize ierr
146  ierr = 0
147 
148  ! set up array sizes
149  dims = shape(theta_e)
150 
151  ! set up local ord
152  ord_loc = 1
153  if (present(ord)) ord_loc = ord
154 
155  ! tests
156  if (dims(3).ne.size(r_f) .or. dims(3).ne.size(r_e)) then
157  ierr = 1
158  err_msg = 'r_F and r_E need to have the correct dimensions'
159  chckerr(err_msg)
160  end if
161  if (present(r_f_array).neqv.present(r_e_array)) then
162  ierr = 1
163  err_msg = 'both r_F_array and r_E_array have to be provided'
164  chckerr(err_msg)
165  end if
166 
167  ! convert normal position
168  ierr = coord_f2e_r(grid_eq,r_f,r_e,r_f_array,r_e_array)
169  chckerr('')
170 
171  ! choose which equilibrium style is being used:
172  ! 1: VMEC
173  ! 2: HELENA
174  select case (eq_style)
175  case (1) ! VMEC
176  ierr = coord_f2e_vmec()
177  chckerr('')
178  case (2) ! HELENA
179  ! trivial HELENA uses flux coordinates
180  theta_e = theta_f
181  zeta_e = zeta_f
182  end select
183  contains
184  ! VMEC version
186  integer function coord_f2e_vmec() result(ierr)
187  use num_vars, only: norm_disc_prec_eq, rank
188  use num_ops, only: calc_zero_hh
189  use vmec_vars, only: mnmax_v
190  use num_utilities, only: spline
191 
192  character(*), parameter :: rout_name = 'coord_F2E_VMEC'
193 
194  ! local variables
195  integer :: id ! counter
196  real(dp), pointer :: loc_r_f(:) => null() ! loc_r in F coords.
197  real(dp), allocatable :: theta_e_guess_3d(:,:,:) ! guess for theta_E
198  real(dp), allocatable :: lam(:,:,:,:) ! lambda and derivative
199 
200  ! initialize ierr
201  ierr = 0
202 
203  ! set up loc_r_F
204  if (present(r_f_array)) then
205  loc_r_f => r_f_array
206  else
207  loc_r_f => grid_eq%loc_r_F
208  end if
209 
210  ! the toroidal angle is trivial
211  zeta_e = - zeta_f ! conversion VMEC LH -> RH coord. system
212 
213  ! allocate local copies of L_V_c and L_V_s
214  allocate(l_v_c_loc(mnmax_v,dims(3)))
215  allocate(l_v_s_loc(mnmax_v,dims(3)))
216 
217  ! interpolate L_V_c and L_V_s at requested normal points r_F
218  do id = 1,mnmax_v
219  ierr = spline(loc_r_f,l_v_c(id,grid_eq%i_min:grid_eq%i_max,0),&
220  &r_f,l_v_c_loc(id,:),ord=norm_disc_prec_eq)
221  chckerr('')
222  ierr = spline(loc_r_f,l_v_s(id,grid_eq%i_min:grid_eq%i_max,0),&
223  &r_f,l_v_s_loc(id,:),ord=norm_disc_prec_eq)
224  chckerr('')
225  end do
226 
227  ! set up guess
228  allocate(theta_e_guess_3d(dims(1),dims(2),dims(3)))
229  allocate(lam(dims(1),dims(2),dims(3),0:1))
230 
231  ! set up guess: try theta_F - lambda(theta_F)/(1+Dlambda(theta_F)),
232  ! which is what you get when you approximate lambda(theta_V) as
233  ! lambda(theta_F) + Dlambda(theta_F) (theta_V-theta_F)
234  ! (lambda(theta_V) would be the exact solution)
235  do id = 0,1
236  ierr = fourier2real(l_v_c_loc,l_v_s_loc,theta_f,zeta_e,&
237  &lam(:,:,:,id),sym=[is_asym_v,.true.],deriv=[id,0])
238  chckerr('')
239  end do
240  theta_e_guess_3d = theta_f - lam(:,:,:,0)/(1+lam(:,:,:,1))
241 
242  ! the poloidal angle has to be found as the zero of
243  ! f = theta_F - theta_E - lambda
244  err_msg = calc_zero_hh(dims,theta_e,fun_pol,ord_loc,&
245  &theta_e_guess_3d,output=.true.)
246  if (err_msg.ne.'') then
247  ierr = 1
248  chckerr(err_msg)
249  end if
250 
251  ! do a check whether the result is indeed zero
252  if (maxval(abs(fun_pol(dims,theta_e,0))).gt.tol_zero*100) then
253  call plot_hdf5('f','ERROR_coord_F2E_VMEC_'//trim(i2str(rank)),&
254  &fun_pol(dims,theta_e,0))
255  err_msg = 'In coord_F2E_VMEC, calculating whether f=0 as a &
256  &check, yields a deviation from 0 of '//trim(r2strt(&
257  &maxval(abs(fun_pol(dims,theta_e,0)))))//'. See output.'
258  ierr = 1
259  chckerr(err_msg)
260  end if
261 
262  ! clean up
263  nullify(loc_r_f)
264  end function coord_f2e_vmec
265 
266  ! function that returns f = theta_F - theta_V - lambda or its
267  ! derivatives in the poloidal direction. It uses zeta_E (= zeta_V) from
268  ! the parent function.
270  function fun_pol(dims,theta_E_in,dpol)
271  character(*), parameter :: rout_name = 'fun_pol'
272 
273  ! input / output
274  integer, intent(in) :: dims(3)
275  real(dp), intent(in) :: theta_e_in(dims(1),dims(2),dims(3))
276  integer, intent(in) :: dpol
277  real(dp) :: fun_pol(dims(1),dims(2),dims(3))
278 
279  ! local variables
280  real(dp) :: lam(dims(1),dims(2),dims(3))
281 
282  ! initialize fun_pol
283  fun_pol = 0.0_dp
284 
285  ! calculate lambda
286  ierr = fourier2real(l_v_c_loc,l_v_s_loc,theta_e_in,zeta_e,lam,&
287  &sym=[is_asym_v,.true.],deriv=[dpol,0])
288  chckerr('')
289 
290  ! calculate the output function
291  if (dpol.eq.0) then
292  fun_pol = theta_f - theta_e_in - lam
293  else if (dpol.eq.1) then
294  fun_pol = -1 - lam
295  else if (dpol.gt.1) then
296  fun_pol = - lam
297  end if
298  end function fun_pol
299  end function coord_f2e_rtz
301  integer function coord_f2e_r(grid_eq,r_F,r_E,r_F_array,r_E_array) &
302  &result(ierr) ! version with only r
303  use num_vars, only: norm_disc_prec_eq
304  use num_utilities, only: spline
305 
306  character(*), parameter :: rout_name = 'coord_F2E_r'
307 
308  ! input / output
309  type(grid_type), intent(in) :: grid_eq
310  real(dp), intent(in) :: r_f(:)
311  real(dp), intent(inout) :: r_e(:)
312  real(dp), intent(in), optional, target :: r_f_array(:)
313  real(dp), intent(in), optional, target :: r_e_array(:)
314 
315  ! local variables
316  character(len=max_str_ln) :: err_msg ! error message
317  integer :: n_r ! dimension of the grid
318  real(dp), pointer :: loc_r_e(:) => null() ! loc_r in E coords.
319  real(dp), pointer :: loc_r_f(:) => null() ! loc_r in F coords.
320 
321  ! initialize ierr
322  ierr = 0
323 
324  ! set up array sizes
325  n_r = size(r_f)
326 
327  ! tests
328  if (n_r.ne.size(r_e)) then
329  ierr = 1
330  err_msg = 'r_F and r_E need to have the correct dimensions'
331  chckerr(err_msg)
332  end if
333  if (present(r_f_array).neqv.present(r_e_array)) then
334  ierr = 1
335  err_msg = 'both r_F_array and r_E_array have to be provided'
336  chckerr(err_msg)
337  end if
338 
339  ! set up loc_r_E and loc_r_F
340  if (present(r_f_array)) then
341  loc_r_e => r_e_array
342  loc_r_f => r_f_array
343  else
344  loc_r_e => grid_eq%loc_r_E
345  loc_r_f => grid_eq%loc_r_F
346  end if
347 
348  ! convert normal position
349  ierr = spline(loc_r_f,loc_r_e,r_f,r_e,ord=norm_disc_prec_eq)
350  chckerr('')
351 
352  ! clean up
353  nullify(loc_r_e,loc_r_f)
354  end function coord_f2e_r
355 
357  integer function coord_e2f_rtz(grid_eq,r_E,theta_E,zeta_E,r_F,&
358  &theta_F,zeta_F,r_E_array,r_F_array) result(ierr) ! version with r, theta and zeta
359  use num_vars, only: eq_style
360 
361  character(*), parameter :: rout_name = 'coord_E2F_rtz'
362 
363  ! input / output
364  type(grid_type), intent(in) :: grid_eq
365  real(dp), intent(in) :: r_e(:)
366  real(dp), intent(in) :: theta_e(:,:,:)
367  real(dp), intent(in) :: zeta_e(:,:,:)
368  real(dp), intent(inout) :: r_f(:)
369  real(dp), intent(inout) :: theta_f(:,:,:)
370  real(dp), intent(inout) :: zeta_f(:,:,:)
371  real(dp), intent(in), optional, target :: r_e_array(:)
372  real(dp), intent(in), optional, target :: r_f_array(:)
373 
374  ! local variables (also used in child functions)
375  character(len=max_str_ln) :: err_msg ! error message
376  integer :: dims(3) ! dimensions of the grid
377 
378  ! initialize ierr
379  ierr = 0
380 
381  ! set up array sizes
382  dims = shape(theta_f)
383 
384  ! tests
385  if (dims(3).ne.size(r_f) .or. dims(3).ne.size(r_e)) then
386  ierr = 1
387  err_msg = 'r_F and r_E need to have the correct dimensions'
388  chckerr(err_msg)
389  end if
390  if (present(r_f_array).neqv.present(r_e_array)) then
391  ierr = 1
392  err_msg = 'both r_F_array and r_E_array have to be provided'
393  chckerr(err_msg)
394  end if
395 
396  ! convert normal position
397  ierr = coord_e2f_r(grid_eq,r_e,r_f,r_e_array,r_f_array)
398  chckerr('')
399 
400  ! choose which equilibrium style is being used:
401  ! 1: VMEC
402  ! 2: HELENA
403  select case (eq_style)
404  case (1) ! VMEC
405  ierr = coord_e2f_vmec()
406  chckerr('')
407  case (2) ! HELENA
408  ! trivial HELENA uses flux coordinates
409  theta_f = theta_e
410  zeta_f = zeta_e
411  end select
412  contains
413  ! VMEC version
415  integer function coord_e2f_vmec() result(ierr)
416  use num_vars, only: norm_disc_prec_eq
417  use vmec_utilities, only: fourier2real
418  use vmec_vars, only: mnmax_v, l_v_c, l_v_s, is_asym_v
419  use num_utilities, only: spline
420 
421  character(*), parameter :: rout_name = 'coord_E2F_VMEC'
422 
423  ! local variables
424  integer :: id ! counter
425  real(dp), allocatable :: l_v_c_loc(:,:), l_v_s_loc(:,:) ! local version of L_V_c and L_V_s
426  real(dp), allocatable :: lam(:,:,:) ! lambda
427  real(dp), pointer :: loc_r_e(:) => null() ! loc_r in E coords.
428 
429  ! initialize ierr
430  ierr = 0
431 
432  ! set up loc_r_E
433  if (present(r_e_array)) then
434  loc_r_e => r_e_array
435  else
436  loc_r_e => grid_eq%loc_r_E
437  end if
438 
439  ! the toroidal angle is trivial
440  zeta_f = - zeta_e ! conversion VMEC LH -> RH coord. system
441 
442  ! allocate local copies of L_V_c and L_V_s and lambda
443  allocate(l_v_c_loc(mnmax_v,dims(3)))
444  allocate(l_v_s_loc(mnmax_v,dims(3)))
445  allocate(lam(dims(1),dims(2),dims(3)))
446 
447  ! interpolate L_V_c and L_V_s at requested normal points r_E
448  do id = 1,mnmax_v
449  ierr = spline(loc_r_e,l_v_c(id,grid_eq%i_min:grid_eq%i_max,0),&
450  &r_e,l_v_c_loc(id,:),ord=norm_disc_prec_eq)
451  chckerr('')
452  ierr = spline(loc_r_e,l_v_s(id,grid_eq%i_min:grid_eq%i_max,0),&
453  &r_e,l_v_s_loc(id,:),ord=norm_disc_prec_eq)
454  chckerr('')
455  end do
456 
457  ! calculate lambda
458  ierr = fourier2real(l_v_c_loc,l_v_s_loc,theta_e,zeta_e,lam,&
459  &sym=[is_asym_v,.true.])
460  chckerr('')
461 
462  ! the poloidal angle has to be found as
463  ! theta_F = theta_E + lambda
464  theta_f = theta_e + lam
465 
466  ! clean up
467  nullify(loc_r_e)
468  end function coord_e2f_vmec
469  end function coord_e2f_rtz
471  integer function coord_e2f_r(grid_eq,r_E,r_F,r_E_array,r_F_array) &
472  &result(ierr) ! version with only r
473  use num_vars, only: norm_disc_prec_eq
474  use num_utilities, only: spline
475 
476  character(*), parameter :: rout_name = 'coord_E2F_r'
477 
478  ! input / output
479  type(grid_type), intent(in) :: grid_eq
480  real(dp), intent(in) :: r_e(:)
481  real(dp), intent(inout) :: r_f(:)
482  real(dp), intent(in), optional, target :: r_e_array(:)
483  real(dp), intent(in), optional, target :: r_f_array(:)
484 
485  ! local variables
486  character(len=max_str_ln) :: err_msg ! error message
487  integer :: n_r ! dimension of the grid
488  real(dp), pointer :: loc_r_e(:) => null() ! loc_r in E coords.
489  real(dp), pointer :: loc_r_f(:) => null() ! loc_r in F coords.
490 
491  ! initialize ierr
492  ierr = 0
493 
494  ! set up array sizes
495  n_r = size(r_e)
496 
497  ! tests
498  if (n_r.ne.size(r_f)) then
499  ierr = 1
500  err_msg = 'r_E and r_F need to have the correct dimensions'
501  chckerr(err_msg)
502  end if
503  if (present(r_e_array).neqv.present(r_f_array)) then
504  ierr = 1
505  err_msg = 'both r_E_array and r_F_array have to be provided'
506  chckerr(err_msg)
507  end if
508 
509  ! set up loc_r_E and loc_r_F
510  if (present(r_f_array)) then
511  loc_r_e => r_e_array
512  loc_r_f => r_f_array
513  else
514  loc_r_e => grid_eq%loc_r_E
515  loc_r_f => grid_eq%loc_r_F
516  end if
517 
518  ! convert normal position
519  ierr = spline(loc_r_e,loc_r_f,r_e,r_f,ord=norm_disc_prec_eq)
520  chckerr('')
521 
522  ! clean up
523  nullify(loc_r_e,loc_r_f)
524  end function coord_e2f_r
525 
527  integer function calc_eqd_grid_3d(var,min_grid,max_grid,grid_dim,&
528  &excl_last) result(ierr)
529  character(*), parameter :: rout_name = 'calc_eqd_grid_3D'
530 
531  ! input and output
532  real(dp), intent(inout) :: var(:,:,:)
533  real(dp), intent(in) :: min_grid
534  real(dp), intent(in) :: max_grid
535  integer, intent(in) :: grid_dim
536  logical, intent(in), optional :: excl_last
537 
538  ! local variables
539  integer :: id ! counter
540  character(len=max_str_ln) :: err_msg ! error message
541  logical :: excl_last_loc ! local copy of excl_last
542  integer :: grid_size ! nr. of points
543  integer :: grid_size_mod ! grid_size or grid_size + 1
544 
545  ! initialize ierr
546  ierr = 0
547 
548  ! tests
549  if (grid_dim.lt.1 .or. grid_dim.gt.3) then
550  ierr = 1
551  err_msg = 'grid_dim has to point to a dimension going from 1 to 3'
552  chckerr(err_msg)
553  end if
554 
555  ! set up grid_size
556  grid_size = size(var,grid_dim)
557 
558  ! test some values
559  if (grid_size.lt.1) then
560  err_msg = 'The angular array has to have a length of at &
561  &least 1'
562  ierr = 1
563  chckerr(err_msg)
564  end if
565 
566  ! set up local excl_last
567  excl_last_loc = .false.
568  if (present(excl_last)) excl_last_loc = excl_last
569 
570  ! add 1 to modified grid size if last point is to be excluded
571  if (excl_last_loc) then
572  grid_size_mod = grid_size + 1
573  else
574  grid_size_mod = grid_size
575  end if
576 
577  ! initialize output vector
578  var = 0.0_dp
579 
580  ! calculate grid points
581  if (grid_size.eq.1) then
582  var = min_grid ! = max_grid
583  else
584  if (grid_dim.eq.1) then
585  do id = 1,grid_size
586  var(id,:,:) = min_grid + &
587  &(max_grid-min_grid)*(id-1)/(grid_size_mod-1)
588  end do
589  else if (grid_dim.eq.2) then
590  do id = 1,grid_size
591  var(:,id,:) = min_grid + &
592  &(max_grid-min_grid)*(id-1)/(grid_size_mod-1)
593  end do
594  else
595  do id = 1,grid_size
596  var(:,:,id) = min_grid + &
597  &(max_grid-min_grid)*(id-1)/(grid_size_mod-1)
598  end do
599  end if
600  end if
601  end function calc_eqd_grid_3d
603  integer function calc_eqd_grid_1d(var,min_grid,max_grid,excl_last) &
604  &result(ierr)
605  character(*), parameter :: rout_name = 'calc_eqd_grid_1D'
606 
607  ! input and output
608  real(dp), intent(inout) :: var(:)
609  real(dp), intent(in) :: min_grid
610  real(dp), intent(in) :: max_grid
611  logical, intent(in), optional :: excl_last
612 
613  ! local variables
614  real(dp), allocatable :: var_3d(:,:,:) ! 3D version of var
615 
616  ! initialize ierr
617  ierr = 0
618 
619  ! set up var_3D
620  allocate(var_3d(size(var),1,1))
621 
622  ! call 3D version
623  ierr = calc_eqd_grid_3d(var_3d,min_grid,max_grid,1,excl_last)
624  chckerr('')
625 
626  ! update var
627  var = var_3d(:,1,1)
628  end function calc_eqd_grid_1d
629 
631  integer function calc_tor_diff_2d(v_com,theta,norm_disc_prec,absolute,r) &
632  &result(ierr)
633  use num_vars, only: min_theta_plot, max_theta_plot, tol_zero
634  use num_utilities, only: spline, order_per_fun
635 
636  character(*), parameter :: rout_name = 'calc_tor_diff_2D'
637 
638  ! input / output
639  real(dp), intent(inout) :: v_com(:,:,:,:,:)
640  real(dp), intent(in) :: theta(:,:,:)
641  integer, intent(in) :: norm_disc_prec
642  logical, intent(in), optional :: absolute
643  real(dp), intent(in), optional :: r(:)
644 
645  ! local variables
646  integer :: id, jd, kd, ld ! counters
647  integer :: n_pol ! number of poloidal points to be used (1 less than total)
648  integer :: istat ! status
649  real(dp), allocatable :: theta_eqd(:) ! equidistant grid
650  real(dp), allocatable :: v_com_interp(:,:) ! interpolated local v_com for outer points
651  real(dp), allocatable :: theta_ord(:) ! ordered theta
652  real(dp), allocatable :: v_com_ord(:) ! ordered v_com
653  logical :: absolute_loc ! local absolute
654 
655  ! initialize ierr
656  ierr = 0
657 
658  ! set up variables
659  n_pol = size(theta,1)-1
660  allocate(theta_eqd(n_pol))
661  allocate(v_com_interp(n_pol,3))
662  absolute_loc = .false.
663  if (present(absolute)) absolute_loc = absolute
664 
665  norm: do kd = 1,size(v_com,3)
666  ! set up equidistant grid
667  ierr = calc_eqd_grid(theta_eqd,min_theta_plot*pi,max_theta_plot*pi,&
668  &excl_last=.true.)
669  chckerr('')
670  do id = 1,size(v_com,4)
671  do ld = 1,size(v_com,5)
672  ! interpolate the geometric poloidal angle for outer points
673  do jd = 1,3,2
674  ! order
675  ierr = order_per_fun(theta(1:n_pol,jd,kd),&
676  &v_com(1:n_pol,jd,kd,id,ld),theta_ord,v_com_ord,&
677  &norm_disc_prec)
678  chckerr('')
679 
680  istat = spline(theta_ord,v_com_ord,theta_eqd,&
681  &v_com_interp(:,jd),ord=min(norm_disc_prec,3),&
682  &deriv=0)
683  deallocate(theta_ord,v_com_ord)
684 
685  if (istat.ne.0) then
686  call display_interp_warning(r)
687  v_com(:,2,kd,:,:) = 0._dp
688  cycle norm
689  end if
690  end do
691 
692  ! calculate difference and save in middle point
693  v_com_interp(:,2) = &
694  &(v_com_interp(:,3)-v_com_interp(:,1))
695  if (.not.absolute_loc) & ! make it relative
696  &v_com_interp(:,2) = v_com_interp(:,2)/&
697  &sign(max(&
698  &tol_zero,abs(v_com_interp(:,3)+v_com_interp(:,1))),&
699  &v_com_interp(:,3)+v_com_interp(:,1))
700 
701  ! order
702  ierr = order_per_fun(theta_eqd,v_com_interp(:,2),&
703  &theta_ord,v_com_ord,norm_disc_prec)
704  chckerr('')
705 
706  ! interpolate back
707  istat = spline(theta_ord,v_com_ord,theta(:,2,kd),&
708  &v_com(:,2,kd,id,ld),min(norm_disc_prec,3),&
709  &deriv=0)
710  deallocate(theta_ord,v_com_ord)
711 
712  if (istat.ne.0) then
713  call display_interp_warning(r)
714  v_com(:,2,kd,:,:) = 0._dp
715  cycle norm
716  end if
717  end do
718  end do
719  end do norm
720  contains
721  ! Displays warning if no interpolation possible.
722  ! Uses variable kd from parent procedure.
724  subroutine display_interp_warning(r)
725  ! input / output
726  real(dp), intent(in), optional :: r(:) ! normal positions for theta
727 
728  if (present(r)) then
729  call writo('Cannot interpolate geometrical &
730  &theta at normal position '//&
731  &trim(r2str(r(kd))),warning=.true.)
732  else
733  call writo('Cannot interpolate geometrical &
734  &theta',warning=.true.)
735  end if
736  call lvl_ud(1)
737  call writo('Are you sure the origin is chosen &
738  &correctly?')
739  call writo('Skipping this normal position')
740  call lvl_ud(-1)
741  end subroutine display_interp_warning
742  end function calc_tor_diff_2d
744  integer function calc_tor_diff_0d(v_mag,theta,norm_disc_prec,absolute,r) &
745  &result(ierr)
746  character(*), parameter :: rout_name = 'calc_tor_diff_0D'
747 
748  ! input / output
749  real(dp), intent(inout) :: v_mag(:,:,:)
750  real(dp), intent(in) :: theta(:,:,:)
751  integer, intent(in) :: norm_disc_prec
752  logical, intent(in), optional :: absolute
753  real(dp), intent(in), optional :: r(:)
754 
755  ! local variable
756  real(dp), allocatable :: v_com_loc(:,:,:,:,:) ! local copy of v_mag
757 
758  ! initialize ierr
759  ierr = 0
760 
761  ! set up local copy of v_mag
762  allocate(v_com_loc(size(v_mag,1),size(v_mag,2),size(v_mag,3),1,1))
763  v_com_loc(:,:,:,1,1) = v_mag
764 
765  ! call 2-D version
766  ierr = calc_tor_diff_2d(v_com_loc,theta,norm_disc_prec,&
767  &absolute=absolute,r=r)
768  chckerr('')
769 
770  ! copy
771  v_mag = v_com_loc(:,:,:,1,1)
772  end function calc_tor_diff_0d
773 
798  integer function calc_xyz_grid(grid_eq,grid_XYZ,X,Y,Z,L,R) result(ierr)
799  use num_vars, only: eq_style, use_normalization
800  use num_utilities, only: round_with_tol
801  use eq_vars, only: r_0
802 
803  character(*), parameter :: rout_name = 'calc_XYZ_grid'
804 
805  ! input / output
806  type(grid_type), intent(in) :: grid_eq
807  type(grid_type), intent(in) :: grid_xyz
808  real(dp), intent(inout) :: x(:,:,:)
809  real(dp), intent(inout) :: y(:,:,:)
810  real(dp), intent(inout) :: z(:,:,:)
811  real(dp), intent(inout), optional :: l(:,:,:)
812  real(dp), intent(inout), optional :: r(:,:,:)
813 
814  ! local variables
815  character(len=max_str_ln) :: err_msg ! error message
816 
817  ! initialize ierr
818  ierr = 0
819 
820  ! test
821  if (size(x,1).ne.grid_xyz%n(1) .or. size(x,2).ne.grid_xyz%n(2) .or. &
822  &size(x,3).ne.grid_xyz%loc_n_r) then
823  ierr = 1
824  err_msg = 'X needs to have the correct dimensions'
825  chckerr(err_msg)
826  end if
827  if (size(y,1).ne.grid_xyz%n(1) .or. size(y,2).ne.grid_xyz%n(2) .or. &
828  &size(y,3).ne.grid_xyz%loc_n_r) then
829  ierr = 1
830  err_msg = 'Y needs to have the correct dimensions'
831  chckerr(err_msg)
832  end if
833  if (size(z,1).ne.grid_xyz%n(1) .or. size(z,2).ne.grid_xyz%n(2) .or. &
834  &size(z,3).ne.grid_xyz%loc_n_r) then
835  ierr = 1
836  err_msg = 'Z needs to have the correct dimensions'
837  chckerr(err_msg)
838  end if
839  if (present(l)) then
840  if (size(l,1).ne.grid_xyz%n(1) .or. size(l,2).ne.grid_xyz%n(2) &
841  &.or. size(l,3).ne.grid_xyz%loc_n_r) then
842  ierr = 1
843  err_msg = 'L needs to have the correct dimensions'
844  chckerr(err_msg)
845  end if
846  end if
847 
848  ! choose which equilibrium style is being used:
849  ! 1: VMEC
850  ! 2: HELENA
851  select case (eq_style)
852  case (1) ! VMEC
853  ierr = calc_xyz_grid_vmec(grid_eq,grid_xyz,x,y,z,l=l,r=r)
854  chckerr('')
855  case (2) ! HELENA
856  ierr = calc_xyz_grid_hel(grid_eq,grid_xyz,x,y,z,r=r)
857  chckerr('')
858  end select
859 
860  ! if normalized, return to physical value
861  if (use_normalization) then
862  x = x*r_0
863  y = y*r_0
864  z = z*r_0
865  if (present(r)) r = r*r_0
866  end if
867  contains
868  ! VMEC version
870  integer function calc_xyz_grid_vmec(grid_eq,grid_XYZ,X,Y,Z,L,R) &
871  &result(ierr)
872  use num_vars, only: norm_disc_prec_eq
873  use vmec_utilities, only: fourier2real
874  use vmec_vars, only: r_v_c, r_v_s, z_v_c, z_v_s, l_v_c, l_v_s, &
875  &mnmax_v, is_asym_v
876  use num_utilities, only: spline
877 
878  character(*), parameter :: rout_name = 'calc_XYZ_grid_VMEC'
879 
880  ! input / output
881  type(grid_type), intent(in) :: grid_eq ! equilibrium grid
882  type(grid_type), intent(in) :: grid_xyz ! grid for which to calculate X, Y, Z and optionally L
883  real(dp), intent(inout) :: x(:,:,:), y(:,:,:), z(:,:,:) ! X, Y and Z of grid
884  real(dp), intent(inout), optional :: l(:,:,:) ! lambda of grid
885  real(dp), intent(inout), optional :: r(:,:,:) ! R of grid
886 
887  ! local variables
888  integer :: id ! counter
889  real(dp), allocatable :: r_v_c_int(:,:), r_v_s_int(:,:) ! interpolated version of R_V_c and R_V_s
890  real(dp), allocatable :: z_v_c_int(:,:), z_v_s_int(:,:) ! interpolated version of Z_V_c and Z_V_s
891  real(dp), allocatable :: l_v_c_int(:,:), l_v_s_int(:,:) ! interpolated version of L_V_c and L_V_s
892  real(dp), allocatable :: r_loc(:,:,:) ! R in Cylindrical coordinates
893  character(len=max_str_ln) :: err_msg ! error message
894 
895  ! initialize ierr
896  ierr = 0
897 
898  ! test
899  if (.not.allocated(grid_xyz%trigon_factors)) then
900  ierr = 1
901  err_msg = 'trigon factors of grid_XYZ need to be allocated'
902  chckerr(err_msg)
903  end if
904 
905  ! set up interpolated R_V_c_int, ..
906  allocate(r_v_c_int(mnmax_v,grid_xyz%loc_n_r))
907  allocate(r_v_s_int(mnmax_v,grid_xyz%loc_n_r))
908  allocate(z_v_c_int(mnmax_v,grid_xyz%loc_n_r))
909  allocate(z_v_s_int(mnmax_v,grid_xyz%loc_n_r))
910  if (present(l)) then
911  allocate(l_v_c_int(mnmax_v,grid_xyz%loc_n_r))
912  allocate(l_v_s_int(mnmax_v,grid_xyz%loc_n_r))
913  end if
914 
915  ! interpolate VMEC tables
916  do id = 1,mnmax_v
917  ierr = spline(grid_eq%r_E,r_v_c(id,:,0),grid_xyz%loc_r_E,&
918  &r_v_c_int(id,:),ord=norm_disc_prec_eq)
919  chckerr('')
920  ierr = spline(grid_eq%r_E,r_v_s(id,:,0),grid_xyz%loc_r_E,&
921  &r_v_s_int(id,:),ord=norm_disc_prec_eq)
922  chckerr('')
923  ierr = spline(grid_eq%r_E,z_v_c(id,:,0),grid_xyz%loc_r_E,&
924  &z_v_c_int(id,:),ord=norm_disc_prec_eq)
925  chckerr('')
926  ierr = spline(grid_eq%r_E,z_v_s(id,:,0),grid_xyz%loc_r_E,&
927  &z_v_s_int(id,:),ord=norm_disc_prec_eq)
928  chckerr('')
929  if (present(l)) then
930  ierr = spline(grid_eq%r_E,l_v_c(id,:,0),grid_xyz%loc_r_E,&
931  &l_v_c_int(id,:),ord=norm_disc_prec_eq)
932  chckerr('')
933  ierr = spline(grid_eq%r_E,l_v_s(id,:,0),grid_xyz%loc_r_E,&
934  &l_v_s_int(id,:),ord=norm_disc_prec_eq)
935  chckerr('')
936  end if
937  end do
938 
939  ! allocate local R
940  allocate(r_loc(grid_xyz%n(1),grid_xyz%n(2),grid_xyz%loc_n_r))
941 
942  ! inverse fourier transform with trigonometric factors
943  ierr = fourier2real(r_v_c_int,r_v_s_int,grid_xyz%trigon_factors,&
944  &r_loc,sym=[.true.,is_asym_v])
945  chckerr('')
946  ierr = fourier2real(z_v_c_int,z_v_s_int,grid_xyz%trigon_factors,z,&
947  &sym=[is_asym_v,.true.])
948  chckerr('')
949  if (present(l)) then
950  ierr = fourier2real(l_v_c_int,l_v_s_int,&
951  &grid_xyz%trigon_factors,l,sym=[is_asym_v,.true.])
952  chckerr('')
953  end if
954  if (present(r)) r = r_loc
955 
956  ! transform cylindrical to cartesian
957  ! (the geometrical zeta is equal to the VMEC zeta)
958  x = r_loc*cos(grid_xyz%zeta_E)
959  y = r_loc*sin(grid_xyz%zeta_E)
960 
961  ! deallocate
962  deallocate(r_v_c_int,r_v_s_int,z_v_c_int,z_v_s_int)
963  if (present(l)) deallocate(l_v_c_int,l_v_s_int)
964  deallocate(r_loc)
965  end function calc_xyz_grid_vmec
966 
967  ! HELENA version
969  integer function calc_xyz_grid_hel(grid_eq,grid_XYZ,X,Y,Z,R) &
970  &result(ierr)
971  use helena_vars, only: r_h, z_h, chi_h, ias, nchi
972  use num_vars, only: norm_disc_prec_eq
973  use num_utilities, only: spline
974 
975  character(*), parameter :: rout_name = 'calc_XYZ_grid_HEL'
976 
977  ! input / output
978  type(grid_type), intent(in) :: grid_eq ! equilibrium grid
979  type(grid_type), intent(in) :: grid_xyz ! grid for which to calculate X, Y, Z and optionally L
980  real(dp), intent(inout) :: x(:,:,:), y(:,:,:), z(:,:,:) ! X, Y and Z of grid
981  real(dp), intent(inout), optional :: r(:,:,:) ! R of grid
982 
983  ! local variables
984  integer :: id, jd, kd ! counters
985  integer :: pmone ! plus or minus one
986  integer :: bcs(2,3) ! boundary conditions (theta(even), theta(odd), r)
987  real(dp) :: bcs_val(2,3) ! values for boundary conditions
988  real(dp), allocatable :: r_h_int(:,:), z_h_int(:,:) ! R and Z at interpolated normal value
989  real(dp), allocatable :: r_loc(:,:,:) ! R in Cylindrical coordinates
990  real(dp) :: theta_loc ! local copy of theta_E
991 
992  ! initialize ierr
993  ierr = 0
994 
995  ! set up boundary conditions
996  if (ias.eq.0) then ! top-bottom symmetric
997  bcs(:,1) = [1,1] ! theta(even): zero first derivative
998  bcs(:,2) = [2,2] ! theta(odd): zero first derivative
999  else
1000  bcs(:,1) = [-1,-1] ! theta(even): periodic
1001  bcs(:,2) = [-1,-1] ! theta(odd): periodic
1002  end if
1003  bcs(:,3) = [0,0] ! r: not-a-knot
1004  bcs_val = 0._dp
1005 
1006  ! allocate local R
1007  allocate(r_loc(grid_xyz%n(1),grid_xyz%n(2),grid_xyz%loc_n_r))
1008 
1009  ! set up interpolated R and Z
1010  allocate(r_h_int(nchi,grid_xyz%loc_n_r))
1011  allocate(z_h_int(nchi,grid_xyz%loc_n_r))
1012 
1013  ! interpolate HELENA output R_H and Z_H for every requested normal
1014  ! point
1015  do id = 1,nchi
1016  ierr = spline(grid_eq%r_E,r_h(id,:),grid_xyz%loc_r_E,&
1017  &r_h_int(id,:),ord=norm_disc_prec_eq,bcs=bcs(:,3),&
1018  &bcs_val=bcs_val(:,3))
1019  chckerr('')
1020  ierr = spline(grid_eq%r_E,z_h(id,:),grid_xyz%loc_r_E,&
1021  &z_h_int(id,:),ord=norm_disc_prec_eq,bcs=bcs(:,3),&
1022  &bcs_val=bcs_val(:,3))
1023  chckerr('')
1024  end do
1025 
1026  ! Note: R_H and Z_H are not adapted to the parallel grid, but
1027  ! tabulated in the original HELENA poloidal grid.
1028  ! loop over normal points
1029  do kd = 1,grid_xyz%loc_n_r ! loop over all normal points
1030  ! loop over toroidal points
1031  do jd = 1,grid_xyz%n(2)
1032  ! interpolate at the requested poloidal points
1033  do id = 1,grid_xyz%n(1)
1034  ! initialize local theta
1035  theta_loc = grid_xyz%theta_E(id,jd,kd)
1036 
1037  ! add or subtract 2pi to the parallel angle until it is
1038  ! at least 0 to get principal range 0..2pi
1039  if (theta_loc.lt.0._dp) then
1040  do while (theta_loc.lt.0._dp)
1041  theta_loc = theta_loc + 2*pi
1042  end do
1043  else if (theta_loc.gt.2*pi) then
1044  do while (theta_loc.gt.2._dp*pi)
1045  theta_loc = theta_loc - 2*pi
1046  end do
1047  end if
1048 
1049  ! take into account possible symmetry
1050  if (ias.eq.0 .and. theta_loc.gt.pi) then
1051  theta_loc = 2*pi-theta_loc
1052  pmone = -1
1053  else
1054  pmone = 1
1055  end if
1056 
1057  ! interpolate the HELENA variables poloidally
1058  ierr = spline(chi_h,r_h_int(:,kd),[theta_loc],&
1059  &r_loc(id:id,jd,kd),ord=norm_disc_prec_eq,&
1060  &bcs=bcs(:,1),bcs_val=bcs_val(:,1)) ! even
1061  chckerr('')
1062  ierr = spline(chi_h,pmone*z_h_int(:,kd),[theta_loc],&
1063  &z(id:id,jd,kd),ord=norm_disc_prec_eq,&
1064  &bcs=bcs(:,2),bcs_val=bcs_val(:,2)) ! odd
1065  chckerr('')
1066  end do
1067  end do
1068  end do
1069  if (present(r)) r = r_loc
1070 
1071  ! calculate X and Y, transforming cylindrical to cartesian
1072  ! (the geometrical zeta is the inverse of HELENA zeta)
1073  x = r_loc*cos(-grid_xyz%zeta_E)
1074  y = r_loc*sin(-grid_xyz%zeta_E)
1075 
1076  ! deallocate
1077  deallocate(r_loc)
1078  end function calc_xyz_grid_hel
1079  end function calc_xyz_grid
1080 
1094  integer function extend_grid_f(grid_in,grid_ext,grid_eq,n_theta_plot,&
1095  &n_zeta_plot,lim_theta_plot,lim_zeta_plot) result(ierr)
1096  use num_vars, only: n_theta => n_theta_plot, n_zeta => n_zeta_plot, &
1097  &min_theta_plot, max_theta_plot, min_zeta_plot, max_zeta_plot
1098  use mpi_utilities, only: get_ser_var
1099 
1100  character(*), parameter :: rout_name = 'extend_grid_F'
1101 
1102  ! input / output
1103  type(grid_type), intent(in) :: grid_in
1104  type(grid_type), intent(inout) :: grid_ext
1105  type(grid_type), intent(in), optional :: grid_eq
1106  integer, intent(in), optional :: n_theta_plot
1107  integer, intent(in), optional :: n_zeta_plot
1108  real(dp), intent(in), optional :: lim_theta_plot(2)
1109  real(dp), intent(in), optional :: lim_zeta_plot(2)
1110 
1111  ! local variables
1112  type(grid_type) :: grid_ext_trim ! trimmed grid_ext
1113  integer :: n_theta_plot_loc ! local n_theta_plot
1114  integer :: n_zeta_plot_loc ! local n_zeta_plot
1115  real(dp) :: lim_theta_plot_loc(2) ! local limits of theta_plot
1116  real(dp) :: lim_zeta_plot_loc(2) ! local limits of zeta_plot
1117 
1118  ! initialize ierr
1119  ierr = 0
1120 
1121  ! set up local variables
1122  n_theta_plot_loc = n_theta
1123  if (present(n_theta_plot)) n_theta_plot_loc = n_theta_plot
1124  n_zeta_plot_loc = n_zeta
1125  if (present(n_zeta_plot)) n_zeta_plot_loc = n_zeta_plot
1126  lim_theta_plot_loc = [min_theta_plot,max_theta_plot]
1127  if (present(lim_theta_plot)) lim_theta_plot_loc = lim_theta_plot
1128  lim_zeta_plot_loc = [min_zeta_plot,max_zeta_plot]
1129  if (present(lim_zeta_plot)) lim_zeta_plot_loc = lim_zeta_plot
1130 
1131  ! user ouput
1132  call writo('Theta limits: '//trim(r2strt(lim_theta_plot_loc(1)))//&
1133  &'pi .. '//trim(r2strt(lim_theta_plot_loc(2)))//'pi')
1134  call writo('Zeta limits: '//trim(r2strt(lim_zeta_plot_loc(1)))//&
1135  &'pi .. '//trim(r2strt(lim_zeta_plot_loc(2)))//'pi')
1136 
1137  ! creating equilibrium grid for the output that covers the whole
1138  ! geometry angularly in F coordinates
1139  ierr = grid_ext%init([n_theta_plot_loc,n_zeta_plot_loc,grid_in%n(3)],&
1140  &[grid_in%i_min,grid_in%i_max])
1141  chckerr('')
1142  grid_ext%r_F = grid_in%r_F
1143  grid_ext%loc_r_F = grid_in%loc_r_F
1144  ierr = calc_eqd_grid(grid_ext%theta_F,lim_theta_plot_loc(1)*pi,&
1145  &lim_theta_plot_loc(2)*pi,1)
1146  chckerr('')
1147  ierr = calc_eqd_grid(grid_ext%zeta_F,lim_zeta_plot_loc(1)*pi,&
1148  &lim_zeta_plot_loc(2)*pi,2)
1149  chckerr('')
1150 
1151  ! convert all F coordinates to E coordinates if requested
1152  if (present(grid_eq)) then
1153  call writo('convert F to E coordinates')
1154  call lvl_ud(1)
1155 
1156  ! set total r_E
1157  grid_ext%r_E = grid_in%r_E
1158 
1159  ! get local r_E and angular theta_E and zeta_E
1160  ierr = coord_f2e(grid_eq,&
1161  &grid_ext%loc_r_F,grid_ext%theta_F,grid_ext%zeta_F,&
1162  &grid_ext%loc_r_E,grid_ext%theta_E,grid_ext%zeta_E)
1163  chckerr('')
1164 
1165  ! trim external grid
1166  ierr = trim_grid(grid_ext,grid_ext_trim)
1167  chckerr('')
1168 
1169  ! clean up
1170  call grid_ext_trim%dealloc()
1171 
1172  call lvl_ud(-1)
1173  end if
1174  end function extend_grid_f
1175 
1188  integer function copy_grid(grid_A,grid_B,lims_B,i_lim) result(ierr)
1189  character(*), parameter :: rout_name = 'copy_grid'
1190 
1191  ! input / output
1192  class(grid_type), intent(in) :: grid_a
1193  class(grid_type), intent(inout) :: grid_b
1194  integer, intent(in), optional :: lims_b(3,2)
1195  integer, intent(in), optional :: i_lim(2)
1196 
1197  ! local variables
1198  integer :: n_b(3) ! n of grid B
1199  integer :: lims_b_loc(3,2) ! local lims_B
1200  integer :: i_lim_loc(2) ! local i_lim
1201 
1202  ! initialize ierr
1203  ierr = 0
1204 
1205  ! set up dimensions and of B
1206  lims_b_loc(1,:) = [1,grid_a%n(1)]
1207  lims_b_loc(2,:) = [1,grid_a%n(2)]
1208  lims_b_loc(3,:) = [1,grid_a%n(3)]
1209  if (present(lims_b)) lims_b_loc = lims_b
1210  n_b = lims_b_loc(:,2)-lims_b_loc(:,1)+1
1211  where (lims_b_loc(:,1).eq.lims_b_loc(:,2) &
1212  &.and. lims_b_loc(:,1).eq.[0,0,0]) n_b = 0
1213  i_lim_loc = lims_b_loc(3,:)
1214  if (present(i_lim)) i_lim_loc = lims_b_loc(3,1) - 1 + i_lim
1215  ! tests
1216  if (n_b(1).gt.grid_a%n(1) .or. n_b(2).gt.grid_a%n(2) .or. &
1217  &n_b(3).gt.grid_a%n(3)) then
1218  write(*,*) 'n_B' ,n_b
1219  write(*,*) 'grid_A', grid_a%n
1220  ierr = 1
1221  chckerr('lims_B is too large')
1222  end if
1223  if (i_lim_loc(1).gt.grid_a%n(3) .or. i_lim_loc(2).gt.grid_a%n(3)) then
1224  ierr = 1
1225  chckerr('normal limits too wide')
1226  end if
1227 
1228  ! allocate the new grid
1229  ierr = grid_b%init(n_b,i_lim)
1230  chckerr('')
1231 
1232  ! copy arrays, possibly subset
1233  grid_b%r_F = grid_a%r_F(lims_b_loc(3,1):lims_b_loc(3,2))
1234  grid_b%r_E = grid_a%r_E(lims_b_loc(3,1):lims_b_loc(3,2))
1235  grid_b%loc_r_F = grid_a%r_F(i_lim_loc(1):i_lim_loc(2))
1236  grid_b%loc_r_E = grid_a%r_E(i_lim_loc(1):i_lim_loc(2))
1237  if (n_b(1).ne.0 .and. n_b(2).ne.0) then
1238  if (grid_a%divided) then
1239  ierr = 1
1240  chckerr('grid_A cannot be divided')
1241  end if
1242  grid_b%theta_F = grid_a%theta_F(lims_b_loc(1,1):lims_b_loc(1,2),&
1243  &lims_b_loc(2,1):lims_b_loc(2,2),i_lim_loc(1):i_lim_loc(2))
1244  grid_b%theta_E = grid_a%theta_E(lims_b_loc(1,1):lims_b_loc(1,2),&
1245  &lims_b_loc(2,1):lims_b_loc(2,2),i_lim_loc(1):i_lim_loc(2))
1246  grid_b%zeta_F = grid_a%zeta_F(lims_b_loc(1,1):lims_b_loc(1,2),&
1247  &lims_b_loc(2,1):lims_b_loc(2,2),i_lim_loc(1):i_lim_loc(2))
1248  grid_b%zeta_E = grid_a%zeta_E(lims_b_loc(1,1):lims_b_loc(1,2),&
1249  &lims_b_loc(2,1):lims_b_loc(2,2),i_lim_loc(1):i_lim_loc(2))
1250  end if
1251  end function copy_grid
1252 
1319  integer function calc_int_vol(ang_1,ang_2,norm,J,f,f_int) result(ierr)
1320 #if ldebug
1321  use num_vars, only: rank, n_procs
1322 #endif
1323 
1324  character(*), parameter :: rout_name = 'calc_int_vol'
1325 
1326  ! input / output
1327  real(dp), intent(in) :: ang_1(:,:,:)
1328  real(dp), intent(in) :: ang_2(:,:,:)
1329  real(dp), intent(in) :: norm(:)
1330  real(dp), intent(in) :: j(:,:,:)
1331  complex(dp), intent(in) :: f(:,:,:,:)
1332  complex(dp), intent(inout) :: f_int(:)
1333 
1334  ! local variables
1335  integer :: id, jd, kd, ld ! counters
1336  integer :: nn_mod ! number of indices for V and V_int
1337  integer :: k_min ! minimum k
1338  integer :: dims(3) ! real dimensions
1339  real(dp), allocatable :: transf_j(:,:,:,:) ! comps. of transf. between J and Jxyz: theta_x, zeta_y, theta_y, zeta_x and r_F_z
1340  real(dp), allocatable :: transf_j_tot(:,:,:) ! transf. between J and Jxyz: (theta_x*zeta_y-theta_y*zeta_x)*r_F_z
1341  complex(dp), allocatable :: jf(:,:,:,:) ! J*f
1342  character(len=max_str_ln) :: err_msg ! error message
1343  integer :: i_a(3), i_b(3) ! limits of building blocks of intermediary variables
1344  logical :: dim_1(2) ! whether the dimension is equal to one
1345 #if ldebug
1346  character(len=max_str_ln), allocatable :: var_names(:,:) ! names of variables
1347  complex(dp), allocatable :: f_int_alt(:) ! alternative calculation for output
1348  complex(dp), allocatable :: f_int_alt_1d(:,:) ! intermediary step in f_int_ALT
1349  complex(dp), allocatable :: f_int_alt_2d(:,:,:) ! intermediary step in f_int_ALT
1350  complex(dp), allocatable :: f_int_alt_alt(:) ! another alternative calculation for output
1351  integer :: loc_norm(2) ! local normal index
1352 #endif
1353 
1354  ! initialize ierr
1355  ierr = 0
1356 
1357  ! set nn_mod
1358  nn_mod = size(f,4)
1359 
1360  ! tests
1361  if (size(f_int).ne.nn_mod) then
1362  ierr = 1
1363  err_msg = 'f and f_int need to have the same storage convention'
1364  chckerr(err_msg)
1365  end if
1366  if (size(ang_1).ne.size(ang_2) .or. size(ang_1,3).ne.size(norm)) then
1367  ierr = 1
1368  err_msg = 'The angular variables have to be compatible'
1369  chckerr(err_msg)
1370  end if
1371  if (size(norm).eq.1) then
1372  ierr = 1
1373  err_msg = 'The normal grid size has to be greater than one'
1374  chckerr(err_msg)
1375  end if
1376 
1377  ! set up dims
1378  dims = shape(ang_1)
1379 
1380  ! set up dim_1
1381  dim_1 = dims(1:2).eq.1
1382 
1383  ! set up Jf and transf_J
1384  allocate(jf(max(dims(1)-1,1),max(dims(2)-1,1),dims(3)-1,nn_mod))
1385  allocate(transf_j(max(dims(1)-1,1),max(dims(2)-1,1),dims(3)-1,5))
1386  allocate(transf_j_tot(max(dims(1)-1,1),max(dims(2)-1,1),dims(3)-1))
1387 
1388  ! set up k_min
1389  k_min = 1
1390 
1391  ! get intermediate Jf and components of transf_J:
1392  ! 1: ang_1_x
1393  ! 2: ang_2_y
1394  ! 3: ang_1_y
1395  ! 4: ang_2_x
1396  ! 5: r_F_z
1397  ! Note: the are always 8 terms in the summation, hence the factor 0.125
1398  jf = 0._dp
1399  transf_j = 0._dp
1400  do kd = 1,2
1401  if (kd.eq.1) then
1402  i_a(3) = 1
1403  i_b(3) = dims(3)-1
1404  else
1405  i_a(3) = 2
1406  i_b(3) = dims(3)
1407  end if
1408  do jd = 1,2
1409  if (dim_1(2)) then
1410  i_a(2) = 1
1411  i_b(2) = dims(2)
1412  else
1413  if (jd.eq.1) then
1414  i_a(2) = 1
1415  i_b(2) = dims(2)-1
1416  else
1417  i_a(2) = 2
1418  i_b(2) = dims(2)
1419  end if
1420  end if
1421  do id = 1,2
1422  if (dim_1(1)) then
1423  i_a(1) = 1
1424  i_b(1) = dims(1)
1425  else
1426  if (id.eq.1) then
1427  i_a(1) = 1
1428  i_b(1) = dims(1)-1
1429  else
1430  i_a(1) = 2
1431  i_b(1) = dims(1)
1432  end if
1433  end if
1434  do ld = 1,nn_mod
1435  ! Jf
1436  jf(:,:,:,ld) = jf(:,:,:,ld) + 0.125_dp * &
1437  &j(i_a(1):i_b(1),i_a(2):i_b(2),i_a(3):i_b(3)) * &
1438  &f(i_a(1):i_b(1),i_a(2):i_b(2),i_a(3):i_b(3),ld)
1439  end do
1440  if (dim_1(1)) then
1441  ! transf_J(1): ang_1_x
1442  transf_j(:,:,:,1) = 2*pi
1443  ! transf_J(4): ang_2_x
1444  transf_j(:,:,:,4) = 0._dp
1445  else
1446  do ld = 1,dims(1)-1
1447  ! transf_J(1): ang_1_x
1448  transf_j(ld,:,:,1) = transf_j(ld,:,:,1) + &
1449  &0.125_dp * ( &
1450  &ang_1(ld+1,i_a(2):i_b(2),i_a(3):i_b(3)) - &
1451  &ang_1(ld,i_a(2):i_b(2),i_a(3):i_b(3)) )
1452  ! transf_J(4): ang_2_x
1453  transf_j(ld,:,:,4) = transf_j(ld,:,:,4) + &
1454  &0.125_dp * ( &
1455  &ang_2(ld+1,i_a(2):i_b(2),i_a(3):i_b(3)) - &
1456  &ang_2(ld,i_a(2):i_b(2),i_a(3):i_b(3)) )
1457  end do
1458  end if
1459  if (dim_1(2)) then
1460  ! transf_J(2): ang_2_y
1461  transf_j(:,:,:,2) = 2*pi
1462  ! transf_J(3): ang_1_y
1463  transf_j(:,:,:,3) = 0._dp
1464  else
1465  do ld = 1,dims(2)-1
1466  ! transf_J(2): ang_2_y
1467  transf_j(:,ld,:,2) = transf_j(:,ld,:,2) + &
1468  &0.125_dp * ( &
1469  &ang_2(i_a(1):i_b(1),ld+1,i_a(3):i_b(3)) - &
1470  &ang_2(i_a(1):i_b(1),ld,i_a(3):i_b(3)) )
1471  ! transf_J(3): ang_1_y
1472  transf_j(:,ld,:,3) = transf_j(:,ld,:,3) + &
1473  &0.125_dp * ( &
1474  &ang_1(i_a(1):i_b(1),ld+1,i_a(3):i_b(3)) - &
1475  &ang_1(i_a(1):i_b(1),ld,i_a(3):i_b(3)) )
1476  end do
1477  end if
1478  do ld = 1,dims(3)-1
1479  ! transf_J(5): r_F_z
1480  transf_j(:,:,ld,5) = transf_j(:,:,ld,5) + 0.125_dp * ( &
1481  &norm(ld+1) - norm(ld) )
1482  end do
1483  end do
1484  end do
1485  end do
1486 
1487  ! set up total transf_J
1488  transf_j_tot = (transf_j(:,:,:,1)*transf_j(:,:,:,2)-&
1489  &transf_j(:,:,:,3)*transf_j(:,:,:,4))*transf_j(:,:,:,5)
1490 
1491  ! integrate term over ang_par_F for all equilibrium grid points, using
1492  ! dx = dy = dz = 1
1493  do ld = 1,nn_mod
1494  f_int(ld) = sum(jf(:,:,:,ld)*transf_j_tot)
1495  end do
1496 
1497 #if ldebug
1498  if (debug_calc_int_vol) then
1499  call writo('Testing whether volume integral is calculated &
1500  &correctly')
1501  call lvl_ud(1)
1502 
1503  allocate(var_names(nn_mod,2))
1504  var_names(:,1) = 'real part of integrand J f'
1505  var_names(:,2) = 'imaginary part of integrand J f'
1506  do ld = 1,nn_mod
1507  var_names(ld,1) = trim(var_names(ld,1))//' '//trim(i2str(ld))
1508  var_names(ld,2) = trim(var_names(ld,2))//' '//trim(i2str(ld))
1509  end do
1510 
1511  call plot_hdf5(var_names(:,1),'TEST_RE_Jf_'//trim(i2str(rank)),&
1512  &rp(jf))
1513  call plot_hdf5(var_names(:,2),'TEST_IM_Jf_'//trim(i2str(rank)),&
1514  &ip(jf))
1515  call plot_hdf5('transformation of Jacobians',&
1516  &'TEST_transf_J_'//trim(i2str(rank)),transf_j_tot)
1517 
1518  ! do alternative calculations
1519  ! assuming that coordinate i varies only in dimensions i!!!
1520  allocate(f_int_alt(nn_mod))
1521  allocate(f_int_alt_1d(size(f,3),nn_mod))
1522  allocate(f_int_alt_2d(size(f,2),size(f,3),nn_mod))
1523  allocate(f_int_alt_alt(nn_mod))
1524  f_int_alt = 0._dp
1525  f_int_alt_1d = 0._dp
1526  f_int_alt_2d = 0._dp
1527  f_int_alt_alt = 0._dp
1528 
1529  ! integrate in first coordinate
1530  do kd = 1,size(f,3)
1531  do jd = 1,size(f,2)
1532  if (size(f,1).gt.1) then
1533  do id = 2,size(f,1)
1534  f_int_alt_2d(jd,kd,:) = f_int_alt_2d(jd,kd,:) + &
1535  &0.5*(j(id,jd,kd)*f(id,jd,kd,:)+&
1536  &j(id-1,jd,kd)*f(id-1,jd,kd,:))*&
1537  &(ang_1(id,jd,kd)-ang_1(id-1,jd,kd))
1538  end do
1539  else
1540  f_int_alt_2d(jd,kd,:) = j(1,jd,kd)*f(1,jd,kd,:)
1541  end if
1542  end do
1543  end do
1544 
1545  ! integrate in second coordinate
1546  do kd = 1,size(f,3)
1547  if (size(f,2).gt.1) then
1548  do jd = 2,size(f,2)
1549  f_int_alt_1d(kd,:) = f_int_alt_1d(kd,:) + &
1550  &0.5*(f_int_alt_2d(jd,kd,:)+&
1551  &f_int_alt_2d(jd-1,kd,:))*&
1552  &(ang_2(1,jd,kd)-ang_2(1,jd-1,kd)) ! assuming that ang_2 is not dependent on dimension 1
1553  end do
1554  else
1555  f_int_alt_1d(kd,:) = f_int_alt_2d(1,kd,:)
1556  end if
1557  end do
1558 
1559  ! integrate in third coordinate
1560  if (size(f,3).gt.1) then
1561  do kd = 2,size(f,3)
1562  f_int_alt(:) = f_int_alt(:) + &
1563  &0.5*(f_int_alt_1d(kd,:)+f_int_alt_1d(kd-1,:))*&
1564  &(norm(kd)-norm(kd-1))
1565  end do
1566  else
1567  f_int_alt = f_int_alt_1d(1,:)
1568  end if
1569 
1570  ! second alternative, first order integration
1571  loc_norm = [1,size(f,3)]
1572  if (rank.lt.n_procs-1) loc_norm(2) = loc_norm(2)-1 ! no ghost region needed for this method
1573  do ld = 1,size(f,4)
1574  f_int_alt_alt(ld) = sum(j(:,:,loc_norm(1):loc_norm(2))*&
1575  &f(:,:,loc_norm(1):loc_norm(2),ld))
1576  end do
1577  if (size(f,1).gt.1) f_int_alt_alt = f_int_alt_alt*&
1578  &(ang_1(2,1,1)-ang_1(1,1,1))
1579  if (size(f,2).gt.1) f_int_alt_alt = f_int_alt_alt*&
1580  &(ang_2(1,2,1)-ang_2(1,1,1))
1581  if (size(f,3).gt.1) f_int_alt_alt = f_int_alt_alt*&
1582  &(norm(2)-norm(1))
1583 
1584  call writo('Resulting integral: ',persistent=.true.)
1585  call lvl_ud(1)
1586  do ld = 1,nn_mod
1587  call writo('rank '//trim(i2str(rank))//' integral '//&
1588  &trim(i2str(ld))//': '//trim(c2str(f_int(ld))),&
1589  &persistent=.true.)
1590  call writo('rank '//trim(i2str(rank))//' alternative '//&
1591  &trim(i2str(ld))//': '//trim(c2str(f_int_alt(ld))),&
1592  &persistent=.true.)
1593  call writo('rank '//trim(i2str(rank))//' other alt. '//&
1594  &trim(i2str(ld))//': '//trim(c2str(f_int_alt_alt(ld))),&
1595  &persistent=.true.)
1596  end do
1597  call lvl_ud(-1)
1598 
1599  call writo('(Note that alternative calculations are only valid if &
1600  &independent coordinates!)')
1601 
1602  call lvl_ud(-1)
1603  end if
1604 #endif
1605 
1606  ! deallocate local variables
1607  deallocate(jf,transf_j,transf_j_tot)
1608  end function calc_int_vol
1609 
1635  integer function trim_grid(grid_in,grid_out,norm_id) result(ierr)
1636  use num_vars, only: n_procs, rank
1637  use mpi_utilities, only: get_ser_var
1638 
1639  character(*), parameter :: rout_name = 'trim_grid'
1640 
1641  ! input / output
1642  type(grid_type), intent(in) :: grid_in
1643  type(grid_type), intent(inout) :: grid_out
1644  integer, intent(inout), optional :: norm_id(2)
1645 
1646  ! local variables
1647  integer, allocatable :: tot_i_min(:) ! i_min of grid of all processes
1648  integer, allocatable :: tot_i_max(:) ! i_max of grid of all processes
1649  integer :: i_lim_out(2) ! i_lim of output grid
1650  integer :: n_out(3) ! n of output grid
1651 
1652  ! initialize ierr
1653  ierr = 0
1654 
1655  ! detect whether grid divided
1656  if (grid_in%divided) then
1657  ! get min_i's of the grid_in
1658  ierr = get_ser_var([grid_in%i_min],tot_i_min,scatter=.true.)
1659  chckerr('')
1660 
1661  ! get max_i's of the grid_in
1662  ierr = get_ser_var([grid_in%i_max],tot_i_max,scatter=.true.)
1663  chckerr('')
1664 
1665  ! set i_lim of trimmed output grid (not yet shifted by first proc
1666  ! min)
1667  if (rank.gt.0) then
1668  i_lim_out(1) = max(tot_i_min(1),tot_i_min(rank+1)+&
1669  &ceiling((tot_i_max(rank)-tot_i_min(rank+1)+1._dp)/2))
1670  else
1671  i_lim_out(1) = tot_i_min(1)
1672  end if
1673  if (rank.lt.n_procs-1) then
1674  i_lim_out(2) = min(tot_i_max(n_procs),tot_i_max(rank+1)-&
1675  &floor((tot_i_max(rank+1)-tot_i_min(rank+2)+1._dp)/2))
1676  else
1677  i_lim_out(2) = tot_i_max(n_procs)
1678  end if
1679 
1680  ! limit to input limits (necessary if a process has only one point)
1681  i_lim_out(1) = max(min(i_lim_out(1),grid_in%i_max),grid_in%i_min)
1682  i_lim_out(2) = max(min(i_lim_out(2),grid_in%i_max),grid_in%i_min)
1683 
1684  ! get min_i's and max_i's of the grid_out, not shifted by min of
1685  ! first process
1686  ierr = get_ser_var([i_lim_out(1)],tot_i_min,scatter=.true.)
1687  chckerr('')
1688  ierr = get_ser_var([i_lim_out(2)],tot_i_max,scatter=.true.)
1689  chckerr('')
1690 
1691  ! set n of output grid
1692  n_out(1) = grid_in%n(1)
1693  n_out(2) = grid_in%n(2)
1694  n_out(3) = tot_i_max(n_procs)-tot_i_min(1)+1
1695 
1696  ! create new grid
1697  ierr = grid_out%init(n_out,i_lim_out-tot_i_min(1)+1,divided=.true.) ! limits shifted by min of first process
1698  chckerr('')
1699 
1700  ! recycle i_lim_out for shifted array indices, set norm_id if
1701  ! requested
1702  i_lim_out = i_lim_out - grid_in%i_min + 1
1703  if (present(norm_id)) norm_id = i_lim_out
1704  else
1705  ! set n of output grid
1706  n_out = grid_in%n
1707 
1708  ! set unshifted i_lim_out and norm_id if requested
1709  i_lim_out = [grid_in%i_min,grid_in%i_max]
1710  if (present(norm_id)) norm_id = i_lim_out
1711 
1712  ! shift i_lim_out by min.
1713  i_lim_out = i_lim_out-i_lim_out(1)+1
1714 
1715  ! create new grid
1716  ierr = grid_out%init(n_out,i_lim_out) ! grid not divided
1717  chckerr('')
1718  end if
1719 
1720  ! copy local arrays
1721  if (grid_in%n(1).ne.0 .and. grid_in%n(2).ne.0) then ! only if 3D grid
1722  grid_out%theta_E = grid_in%theta_E(:,:,i_lim_out(1):i_lim_out(2))
1723  grid_out%zeta_E = grid_in%zeta_E(:,:,i_lim_out(1):i_lim_out(2))
1724  grid_out%theta_F = grid_in%theta_F(:,:,i_lim_out(1):i_lim_out(2))
1725  grid_out%zeta_F = grid_in%zeta_F(:,:,i_lim_out(1):i_lim_out(2))
1726  end if
1727  if (grid_in%divided) then ! but if input grid divided, loc_r gets priority
1728  grid_out%loc_r_E = grid_in%loc_r_E(i_lim_out(1):i_lim_out(2))
1729  grid_out%loc_r_F = grid_in%loc_r_F(i_lim_out(1):i_lim_out(2))
1730  end if
1731 
1732  ! if divided, set total arrays
1733  if (grid_in%divided) then
1734  grid_out%r_E = grid_in%r_E(tot_i_min(1):tot_i_max(n_procs))
1735  grid_out%r_F = grid_in%r_F(tot_i_min(1):tot_i_max(n_procs))
1736  else
1737  grid_out%r_E = grid_in%r_E
1738  grid_out%r_F = grid_in%r_F
1739  end if
1740 
1741  ! if allocated, copy trigonometric factors
1742  if (allocated(grid_in%trigon_factors)) then
1743  allocate(grid_out%trigon_factors(&
1744  &size(grid_in%trigon_factors,1),&
1745  &size(grid_in%trigon_factors,2),&
1746  &size(grid_in%trigon_factors,3),&
1747  &grid_out%loc_n_r,2))
1748  grid_out%trigon_factors = &
1749  &grid_in%trigon_factors(:,:,:,i_lim_out(1):i_lim_out(2),:)
1750  end if
1751  end function trim_grid
1752 
1763  integer function untrim_grid(grid_in,grid_out,size_ghost) result(ierr)
1764  use num_vars, only: n_procs, rank
1765  use mpi_utilities, only: get_ghost_arr, get_ser_var
1766 
1767  character(*), parameter :: rout_name = 'untrim_grid'
1768 
1769  ! input / output
1770  type(grid_type), intent(in) :: grid_in
1771  type(grid_type), intent(inout) :: grid_out
1772  integer, intent(in) :: size_ghost
1773 
1774  ! local variables
1775  integer :: i_lim_in(2) ! limits of input grid
1776  integer :: i_lim_out(2) ! limits of ghosted grid
1777  integer, allocatable :: tot_loc_n_r(:) ! loc_n_r of all processes
1778  integer :: size_ghost_loc ! local size_ghost
1779 
1780  ! initialize ierr
1781  ierr = 0
1782 
1783  ! get array sizes of all processes
1784  ierr = get_ser_var([grid_in%loc_n_r],tot_loc_n_r,scatter=.true.)
1785  chckerr('')
1786 
1787  ! set local size_ghost
1788  size_ghost_loc = min(size_ghost,minval(tot_loc_n_r)-1)
1789 
1790  ! set i_lim_in and i_lim_out
1791  i_lim_in = [grid_in%i_min,grid_in%i_max]
1792  i_lim_out = i_lim_in
1793  if (rank.lt.n_procs-1) i_lim_out(2) = i_lim_out(2)+size_ghost_loc
1794 
1795  ! create grid
1796  ierr = grid_out%init(grid_in%n,i_lim_out)
1797  chckerr('')
1798 
1799  ! set ghosted variables
1800  if (grid_in%n(1).ne.0 .and. grid_in%n(2).ne.0) then ! only if 3d grid
1801  grid_out%theta_e(:,:,1:grid_in%loc_n_r) = grid_in%theta_e
1802  grid_out%zeta_e(:,:,1:grid_in%loc_n_r) = grid_in%zeta_e
1803  grid_out%theta_f(:,:,1:grid_in%loc_n_r) = grid_in%theta_f
1804  grid_out%zeta_f(:,:,1:grid_in%loc_n_r) = grid_in%zeta_f
1805  ierr = get_ghost_arr(grid_out%theta_e,size_ghost_loc)
1806  chckerr('')
1807  ierr = get_ghost_arr(grid_out%zeta_e,size_ghost_loc)
1808  chckerr('')
1809  ierr = get_ghost_arr(grid_out%theta_f,size_ghost_loc)
1810  chckerr('')
1811  ierr = get_ghost_arr(grid_out%zeta_f,size_ghost_loc)
1812  chckerr('')
1813  end if
1814  if (grid_in%divided) then ! but if input grid divided, loc_r gets priority
1815  grid_out%loc_r_e = grid_in%r_e(i_lim_out(1):i_lim_out(2))
1816  grid_out%loc_r_f = grid_in%r_f(i_lim_out(1):i_lim_out(2))
1817  end if
1818  grid_out%r_e = grid_in%r_e
1819  grid_out%r_f = grid_in%r_f
1820  end function untrim_grid
1821 
1856  integer function calc_vec_comp(grid,grid_eq,eq_1,eq_2,v_com,norm_disc_prec,&
1857  &v_mag,base_name,max_transf,v_flux_tor,v_flux_pol,XYZ,compare_tor_pos) &
1858  &result(ierr)
1859 
1860  use eq_vars, only: eq_1_type, eq_2_type, max_flux_f
1861  use eq_utilities, only: calc_inv_met
1862  use num_vars, only: eq_jobs_lims, eq_job_nr, use_pol_flux_f, eq_style, &
1863  &use_normalization, rank, tol_zero, rz_0, &
1864  &compare_tor_pos_glob => compare_tor_pos
1865  use num_utilities, only: c, calc_int, order_per_fun, spline
1866  use eq_vars, only: r_0, b_0, psi_0
1867  use vmec_utilities, only: calc_trigon_factors
1868  use mpi_utilities, only: get_ser_var
1869 
1870  character(*), parameter :: rout_name = 'calc_vec_comp'
1871 
1872  ! input / output
1873  type(grid_type), intent(in) :: grid
1874  type(grid_type), intent(in) :: grid_eq
1875  type(eq_1_type), intent(in) :: eq_1
1876  type(eq_2_type), intent(in) :: eq_2
1877  real(dp), intent(inout) :: v_com(:,:,:,:,:)
1878  integer, intent(in) :: norm_disc_prec
1879  real(dp), intent(inout), optional :: v_mag(:,:,:)
1880  character(len=*), intent(in), optional :: base_name
1881  integer, intent(in), optional :: max_transf
1882  real(dp), intent(inout), allocatable, optional :: v_flux_pol(:,:)
1883  real(dp), intent(inout), allocatable, optional :: v_flux_tor(:,:)
1884  real(dp), intent(in), optional :: xyz(:,:,:,:)
1885  logical, intent(in), optional :: compare_tor_pos
1886 
1887  ! local variables
1888  type(grid_type) :: grid_trim ! trimmed plot grid
1889  integer :: max_transf_loc ! local max_transf
1890  integer :: id, jd, kd, ld ! counters
1891  integer :: plot_dim(4) ! dimensions of plot
1892  integer :: plot_offset(4) ! local offset of plot
1893  integer :: c_loc ! local c
1894  integer :: tor_id(2) ! toroidal indices
1895  integer :: norm_id(2) ! untrimmed normal indices for trimmed grids
1896  logical :: cont_plot ! continued plot
1897  logical :: do_plot ! perform plotting
1898  logical :: compare_tor_pos_loc ! local compare_tor_pos
1899  character(len=max_str_ln) :: description(3) ! description of plots
1900  character(len=max_str_ln) :: file_names(3) ! plot file names
1901  character(len=max_str_ln) :: var_names(3,2) ! variable names
1902  character(len=5) :: coord_names(3) ! name of coordinates
1903  character(len=max_str_ln) :: err_msg ! error message
1904  real(dp) :: norm_len ! normalization factor for lengths, to cancel the one introduced in "calc_XYZ_grid"
1905  real(dp), allocatable :: q_saf(:,:) ! interpolated q_saf_FD and derivative
1906  real(dp), allocatable :: jac(:,:,:) ! interpolated jac_FD
1907  real(dp), allocatable :: xyzr(:,:,:,:) ! X, Y, Z and R of surface in cylindrical coordinates, untrimmed grid
1908  real(dp), allocatable :: theta_geo(:,:,:) ! geometrical theta
1909  real(dp), allocatable :: x(:,:,:,:), y(:,:,:,:), z(:,:,:,:) ! copy of X, Y and Z, trimmed grid
1910  real(dp), allocatable :: v_temp(:,:,:,:,:) ! temporary variable for v
1911  real(dp), allocatable :: v_ser_temp(:) ! temporary serial variable
1912  real(dp), allocatable :: v_ser_temp_int(:) ! temporary integrated serial variable
1913  real(dp), allocatable :: t_ba(:,:,:,:,:,:,:), t_ab(:,:,:,:,:,:,:) ! transformation matrices from A to B
1914  real(dp), allocatable :: d1r(:,:,:), d2r(:,:,:) ! dR/dpsi, dR/dtheta in HELENA coords.
1915  real(dp), allocatable :: d1z(:,:,:), d2z(:,:,:) ! dZ/dpsi, dZ/dtheta in HELENA coords.
1916 
1917  ! initialize ierr
1918  ierr = 0
1919 
1920  ! user output
1921  call writo('Prepare calculation of vector components')
1922  call lvl_ud(1)
1923 
1924  ! set up maximum level up to which to transform and whether to plot
1925  max_transf_loc = 5
1926  if (present(max_transf)) then
1927  if (max_transf.ge.1 .and. max_transf.le.5) &
1928  &max_transf_loc = max_transf
1929  end if
1930  do_plot = .false.
1931  if (present(base_name)) then
1932  if (trim(base_name).ne.'') do_plot = .true.
1933  end if
1934 
1935  ! set up continued plot
1936  cont_plot = eq_job_nr.gt.1
1937 
1938  ! set up normalization factor
1939  norm_len = 1._dp
1940  if (use_normalization) norm_len = r_0
1941 
1942  ! trim plot grid
1943  ierr = trim_grid(grid,grid_trim,norm_id)
1944  chckerr('')
1945 
1946  ! set up X, Y, Z and R in grid
1947  allocate(xyzr(grid%n(1),grid%n(2),grid%loc_n_r,4))
1948  ierr = calc_xyz_grid(grid_eq,grid,xyzr(:,:,:,1),xyzr(:,:,:,2),&
1949  &xyzr(:,:,:,3),r=xyzr(:,:,:,4))
1950  chckerr('')
1951 
1952  ! set up local compare_tor_pos
1953  compare_tor_pos_loc = compare_tor_pos_glob
1954  if (present(compare_tor_pos)) compare_tor_pos_loc = compare_tor_pos
1955 
1956  ! get geometrical poloidal angle if comparing toroidal position
1957  if (compare_tor_pos_loc) then
1958  allocate(theta_geo(grid%n(1),grid%n(2),grid%loc_n_r))
1959  theta_geo = atan2(xyzr(:,:,:,3)-rz_0(2),xyzr(:,:,:,4)-rz_0(1))
1960  where (theta_geo.lt.0) theta_geo = theta_geo + 2*pi
1961  end if
1962 
1963  ! if plotting is required
1964  if (do_plot) then
1965  ! set up plot dimensions and local dimensions
1966  plot_dim = [grid_trim%n,3]
1967  plot_offset = [0,0,grid_trim%i_min-1,0]
1968  tor_id = [1,size(v_com,2)]
1969  if (compare_tor_pos_loc) then
1970  if (plot_dim(2).ne.3) then
1971  ierr = 1
1972  err_msg = 'When comparing toroidal positions, need to &
1973  &have 3 toroidal points (one in the middle)'
1974  chckerr(err_msg)
1975  end if
1976  plot_dim(2) = 1
1977  tor_id = [2,2]
1978  end if
1979 
1980  ! possibly modify if multiple equilibrium parallel jobs
1981  if (size(eq_jobs_lims,2).gt.1) then
1982  plot_dim(1) = eq_jobs_lims(2,size(eq_jobs_lims,2)) - &
1983  &eq_jobs_lims(1,1) + 1
1984  plot_offset(1) = eq_jobs_lims(1,eq_job_nr) - 1
1985  if (compare_tor_pos_loc) then
1986  ierr = 1
1987  err_msg = 'When comparing toroidal positions, cannot have &
1988  &multiple equilibrium jobs'
1989  chckerr(err_msg)
1990  end if
1991  end if
1992 
1993  ! tests
1994  if (eq_style.eq.1 .and. .not.allocated(grid%trigon_factors)) then
1995  ierr = 1
1996  err_msg = 'trigonometric factors not allocated'
1997  chckerr(err_msg)
1998  end if
1999 
2000  ! user output
2001  if (compare_tor_pos_loc) call writo('Comparing toroidal positions')
2002 
2003  ! copy X, Y and Z to trimmed grid copies
2004  allocate(x(grid_trim%n(1),grid_trim%n(2),grid_trim%loc_n_r,1))
2005  allocate(y(grid_trim%n(1),grid_trim%n(2),grid_trim%loc_n_r,1))
2006  allocate(z(grid_trim%n(1),grid_trim%n(2),grid_trim%loc_n_r,1))
2007  if (present(xyz)) then
2008  x(:,:,:,1) = xyz(:,:,norm_id(1):norm_id(2),1)
2009  y(:,:,:,1) = xyz(:,:,norm_id(1):norm_id(2),2)
2010  z(:,:,:,1) = xyz(:,:,norm_id(1):norm_id(2),3)
2011  else
2012  x(:,:,:,1) = xyzr(:,:,norm_id(1):norm_id(2),1)
2013  y(:,:,:,1) = xyzr(:,:,norm_id(1):norm_id(2),2)
2014  z(:,:,:,1) = xyzr(:,:,norm_id(1):norm_id(2),3)
2015  end if
2016  end if
2017 
2018  ! set up temporal interpolated copy of v, T_BA and T_AB
2019  allocate(v_temp(grid%n(1),grid%n(2),grid%loc_n_r,3,2))
2020  allocate(t_ba(grid%n(1),grid%n(2),grid%loc_n_r,9,0:0,0:0,0:0))
2021  allocate(t_ab(grid%n(1),grid%n(2),grid%loc_n_r,9,0:0,0:0,0:0))
2022  allocate(q_saf(grid%loc_n_r,0:1))
2023  if ((present(v_flux_tor) .or. present(v_flux_pol)) .and. &
2024  &.not.compare_tor_pos_loc) then
2025  allocate(jac(grid%n(1),grid%n(2),grid%loc_n_r))
2026  if (rank.eq.0 .and. .not.cont_plot) then
2027  if (present(v_flux_tor)) then
2028  allocate(v_flux_tor(grid_trim%n(3),plot_dim(2)))
2029  v_flux_tor = 0._dp
2030  end if
2031  if (present(v_flux_pol)) then
2032  allocate(v_flux_pol(grid_trim%n(3),plot_dim(1)))
2033  v_flux_pol = 0._dp
2034  end if
2035  end if
2036  end if
2037 
2038  call lvl_ud(-1)
2039 
2040  ! 1. Flux coordinates (alpha,psi,theta)
2041  call writo('Flux coordinates')
2042  call lvl_ud(1)
2043  if (present(v_mag)) then
2044  v_mag = 0._dp
2045  do id = 1,3
2046  v_mag = v_mag + v_com(:,:,:,id,1)*v_com(:,:,:,id,2)
2047  end do
2048  v_mag = sqrt(v_mag)
2049  end if
2050 
2051  ! save temporary copy, normalized
2052  v_temp = v_com
2053 
2054  ! set up plot variables
2055  if (do_plot .and. .not.compare_tor_pos_loc) then ! comparison only works for periodic quantities
2056  coord_names(1) = 'alpha'
2057  coord_names(2) = 'psi'
2058  if (use_pol_flux_f) then
2059  coord_names(3) = 'theta'
2060  else
2061  coord_names(3) = 'zeta'
2062  end if
2063  var_names = trim(base_name)
2064  do id = 1,3
2065  var_names(id,1) = trim(var_names(id,1))//'_sub_'//&
2066  &trim(coord_names(id))
2067  var_names(id,2) = trim(var_names(id,2))//'_sup_'//&
2068  &trim(coord_names(id))
2069  end do
2070  description(1) = 'covariant components of the magnetic field in &
2071  &Flux coordinates'
2072  description(2) = 'contravariant components of the magnetic field &
2073  &in Flux coordinates'
2074  description(3) = 'magnitude of the magnetic field in Flux &
2075  &coordinates'
2076  file_names(1) = trim(base_name)//'_F_sub'
2077  file_names(2) = trim(base_name)//'_F_sup'
2078  file_names(3) = trim(base_name)//'_F_mag'
2079  if (use_normalization) then
2080  v_com(:,:,:,1,1) = v_com(:,:,:,1,1) * r_0 ! norm factor for e_alpha
2081  v_com(:,:,:,2,1) = v_com(:,:,:,2,1) / (r_0*b_0) ! norm factor for e_psi
2082  v_com(:,:,:,3,1) = v_com(:,:,:,3,1) * r_0 ! norm factor for e_theta
2083  v_com(:,:,:,1,2) = v_com(:,:,:,1,2) / r_0 ! norm factor for e^alpha
2084  v_com(:,:,:,2,2) = v_com(:,:,:,2,2) * (r_0*b_0) ! norm factor for e^psi
2085  v_com(:,:,:,3,2) = v_com(:,:,:,3,2) / r_0 ! norm factor for e^theta
2086  end if
2087  do id = 1,2
2088  call plot_hdf5(var_names(:,id),trim(file_names(id)),&
2089  &v_com(:,tor_id(1):tor_id(2),norm_id(1):norm_id(2),:,id),&
2090  &tot_dim=plot_dim,loc_offset=plot_offset,&
2091  &x=x(:,tor_id(1):tor_id(2),:,:),&
2092  &y=y(:,tor_id(1):tor_id(2),:,:),&
2093  &z=z(:,tor_id(1):tor_id(2),:,:),&
2094  &cont_plot=cont_plot,descr=description(id))
2095  end do
2096  if (present(v_mag)) then
2097  call plot_hdf5(trim(base_name),trim(file_names(3)),&
2098  &v_mag(:,tor_id(1):tor_id(2),norm_id(1):norm_id(2)),&
2099  &tot_dim=plot_dim(1:3),loc_offset=plot_offset(1:3),&
2100  &x=x(:,tor_id(1):tor_id(2),:,1),&
2101  &y=y(:,tor_id(1):tor_id(2),:,1),&
2102  &z=z(:,tor_id(1):tor_id(2),:,1),&
2103  &cont_plot=cont_plot,descr=description(3))
2104  end if
2105  end if
2106 
2107  call lvl_ud(-1)
2108 
2109  ! 2. Magnetic coordinates (phi,theta,zeta) by converting using
2110  ! transformation matrices:
2111  ! (v_i)_M = T_M^F (v_i)_F
2112  ! (v^i)_M = (v^i)_F T_F^M
2113  ! where
2114  ! ( -q' theta 1 0 )
2115  ! T_MF = ( -q 0 1 )
2116  ! ( 1 0 0 )
2117  ! if the poloidal flux is used, or
2118  ! ( -q'/q zeta q 0 )
2119  ! T_MF = ( -1 0 0 )
2120  ! ( 1/q 0 1 )
2121  ! if it is the toroidal flux. Its inverse can be calculated as well.
2122  call writo('magnetic coordinates')
2123  call lvl_ud(1)
2124  t_ba = 0._dp
2125  t_ab = 0._dp
2126  do id = 0,1
2127  ierr = spline(grid_eq%loc_r_F,eq_1%q_saf_FD(:,id),grid%loc_r_F,&
2128  &q_saf(:,id),ord=norm_disc_prec)
2129  chckerr('')
2130  end do
2131  c_loc = c([1,1],.false.)
2132  if (use_pol_flux_f) then
2133  t_ba(:,:,:,c_loc,0,0,0) = -grid%theta_F
2134  do kd = 1,grid%loc_n_r
2135  t_ba(:,:,kd,c_loc,0,0,0) = t_ba(:,:,kd,c_loc,0,0,0)*q_saf(kd,1)
2136  t_ba(:,:,kd,c([2,1],.false.),0,0,0) = -q_saf(kd,0)
2137  end do
2138  t_ba(:,:,:,c([3,1],.false.),0,0,0) = 1._dp
2139  t_ba(:,:,:,c([1,2],.false.),0,0,0) = 1._dp
2140  t_ba(:,:,:,c([2,3],.false.),0,0,0) = 1._dp
2141  else
2142  t_ba(:,:,:,c_loc,0,0,0) = -grid%zeta_F
2143  do kd = 1,grid%loc_n_r
2144  t_ba(:,:,kd,c_loc,0,0,0) = &
2145  &t_ba(:,:,kd,c_loc,0,0,0)*q_saf(kd,1)/q_saf(kd,0)
2146  t_ba(:,:,kd,c([2,1],.false.),0,0,0) = 1._dp/q_saf(kd,0)
2147  t_ba(:,:,kd,c([1,2],.false.),0,0,0) = q_saf(kd,0)
2148  end do
2149  t_ba(:,:,:,c([2,1],.false.),0,0,0) = -1._dp
2150  t_ba(:,:,:,c([3,3],.false.),0,0,0) = 1._dp
2151  end if
2152  ierr = calc_inv_met(t_ab,t_ba,[0,0,0])
2153  chckerr('')
2154  v_com = 0._dp
2155  do jd = 1,3
2156  do id = 1,3
2157  v_com(:,:,:,jd,1) = v_com(:,:,:,jd,1) + t_ba(:,:,:,&
2158  &c([jd,id],.false.),0,0,0) * v_temp(:,:,:,id,1)
2159  v_com(:,:,:,jd,2) = v_com(:,:,:,jd,2) + t_ab(:,:,:,&
2160  &c([id,jd],.false.),0,0,0) * v_temp(:,:,:,id,2)
2161  end do
2162  end do
2163  if (present(v_mag)) then
2164  v_mag = 0._dp
2165  do id = 1,3
2166  v_mag = v_mag + v_com(:,:,:,id,1)*v_com(:,:,:,id,2)
2167  end do
2168  v_mag = sqrt(v_mag)
2169  end if
2170 
2171  ! set up plot variables
2172  if (do_plot) then
2173  coord_names(1) = 'psi'
2174  coord_names(2) = 'theta'
2175  coord_names(3) = 'zeta'
2176  var_names = trim(base_name)
2177  do id = 1,3
2178  var_names(id,1) = trim(var_names(id,1))//'_sub_'//&
2179  &trim(coord_names(id))
2180  var_names(id,2) = trim(var_names(id,2))//'_sup_'//&
2181  &trim(coord_names(id))
2182  end do
2183  description(1) = 'covariant components of the magnetic field in &
2184  &Magnetic coordinates'
2185  description(2) = 'contravariant components of the magnetic field &
2186  &in Magnetic coordinates'
2187  description(3) = 'magnitude of the magnetic field in Magnetic &
2188  &coordinates'
2189  file_names(1) = trim(base_name)//'_M_sub'
2190  file_names(2) = trim(base_name)//'_M_sup'
2191  file_names(3) = trim(base_name)//'_M_mag'
2192  if (compare_tor_pos_loc) then
2193  do id = 1,3
2194  file_names(id) = trim(file_names(id))//'_COMP'
2195  end do
2196  end if
2197  if (use_normalization) then
2198  v_com(:,:,:,1,1) = v_com(:,:,:,1,1) / (r_0*b_0) ! norm factor for e_psi
2199  v_com(:,:,:,2,1) = v_com(:,:,:,2,1) * r_0 ! norm factor for e_theta
2200  v_com(:,:,:,3,1) = v_com(:,:,:,3,1) * r_0 ! norm factor for e_zeta
2201  v_com(:,:,:,1,2) = v_com(:,:,:,1,2) * (r_0*b_0) ! norm factor for e^psi
2202  v_com(:,:,:,2,2) = v_com(:,:,:,2,2) / r_0 ! norm factor for e^theta
2203  v_com(:,:,:,3,2) = v_com(:,:,:,3,2) / r_0 ! norm factor for e^zeta
2204  end if
2205  if (compare_tor_pos_loc) then
2206  ierr = calc_tor_diff(v_com,theta_geo,norm_disc_prec)
2207  chckerr('')
2208  end if
2209  do id = 1,2
2210  call plot_hdf5(var_names(:,id),trim(file_names(id)),&
2211  &v_com(:,tor_id(1):tor_id(2),norm_id(1):norm_id(2),:,id),&
2212  &tot_dim=plot_dim,loc_offset=plot_offset,&
2213  &x=x(:,tor_id(1):tor_id(2),:,:),&
2214  &y=y(:,tor_id(1):tor_id(2),:,:),&
2215  &z=z(:,tor_id(1):tor_id(2),:,:),&
2216  &cont_plot=cont_plot,descr=description(id))
2217  end do
2218  if (present(v_mag)) then
2219  if (compare_tor_pos_loc) then
2220  ierr = calc_tor_diff(v_mag,theta_geo,norm_disc_prec)
2221  chckerr('')
2222  end if
2223  call plot_hdf5(trim(base_name),trim(file_names(3)),&
2224  &v_mag(:,tor_id(1):tor_id(2),norm_id(1):norm_id(2)),&
2225  &tot_dim=plot_dim(1:3),loc_offset=plot_offset(1:3),&
2226  &x=x(:,tor_id(1):tor_id(2),:,1),&
2227  &y=y(:,tor_id(1):tor_id(2),:,1),&
2228  &z=z(:,tor_id(1):tor_id(2),:,1),&
2229  &cont_plot=cont_plot,descr=description(3))
2230  end if
2231  end if
2232 
2233  if ((present(v_flux_tor) .or. present(v_flux_pol)) .and. &
2234  &.not.compare_tor_pos_loc) then
2235  ! tests
2236  if (maxval(grid%theta_F(grid%n(1),:,:)).lt.&
2237  &minval(grid%theta_F(1,:,:))) then
2238  call writo('theta of the grid monotomically decreases in Flux &
2239  &quantities.',alert=.true.)
2240  call lvl_ud(1)
2241  call writo('This inverts the sign of the toroidal flux.')
2242  call writo('Remember that the grid limits are prescribed &
2243  &in Flux quantities.')
2244  call lvl_ud(-1)
2245  end if
2246  if (maxval(grid%zeta_F(:,grid%n(2),:)).lt.&
2247  &minval(grid%zeta_F(:,1,:))) then
2248  call writo('zeta of the grid monotomically decreases in Flux &
2249  &quantities.',alert=.true.)
2250  call lvl_ud(1)
2251  call writo('This inverts the sign of the poloidal flux.')
2252  call writo('Remember that the grid limits are prescribed &
2253  &in Flux quantities.')
2254  if (eq_style.eq.1) call writo('For VMEC, these are inverse.')
2255  call lvl_ud(-1)
2256  end if
2257  if (grid_trim%r_F(grid_trim%n(3)).lt.grid_trim%r_F(1)) then
2258  call writo('r of the grid monotomically decreases in Flux &
2259  &quantities.',alert=.true.)
2260  call lvl_ud(1)
2261  call writo('This inverts the sign of the fluxes.')
2262  call writo('Remember that the grid limits are prescribed &
2263  &in Flux quantities.')
2264  call lvl_ud(-1)
2265  end if
2266  if (abs(minval(grid_trim%r_F)).gt.tol_zero) then
2267  call writo('r of the grid does not start at zero.',alert=.true.)
2268  call lvl_ud(1)
2269  call writo('This leaves out part of the fluxes.')
2270  call lvl_ud(-1)
2271  end if
2272 
2273  ! initialize temporary variable
2274  allocate(v_ser_temp_int(grid_trim%n(3)))
2275 
2276  ! set up plot variables
2277  var_names(1,2) = 'integrated poloidal flux of '//trim(base_name)
2278  var_names(2,2) = 'integrated toroidal flux of '//trim(base_name)
2279  file_names(1) = trim(base_name)//'_flux_pol_int'
2280  file_names(2) = trim(base_name)//'_flux_tor_int'
2281 
2282 #if ldebug
2283  if (debug_calc_vec_comp) then
2284  ! user output
2285  call writo('Instead of calculating fluxes, returning &
2286  &integrals in the angular variables.',alert=.true.)
2287  call lvl_ud(1)
2288  call writo('Useful to check whether Maxwell holds.')
2289  call writo('i.e. whether loop integral of J returns &
2290  &B flux')
2291  call writo('Note that the J-variables now refer to B and &
2292  &vice-versa')
2293  call writo('Don''t forget the contribution of the toroidal &
2294  &field B_zeta on axis, times 2piR.')
2295  call writo('Don''t forget the vacuum permeability!')
2296  call lvl_ud(-1)
2297  else
2298 #endif
2299  ! multiply angular contravariant coordinates by Jacobian
2300  do jd = 1,grid_eq%n(2)
2301  do id = 1,grid_eq%n(1)
2302  ierr = spline(grid_eq%loc_r_F,&
2303  &eq_2%jac_FD(id,jd,:,0,0,0),grid%loc_r_F,&
2304  &jac(id,jd,:),ord=norm_disc_prec)
2305  chckerr('')
2306  end do
2307  end do
2308  if (use_normalization) jac = jac*r_0/b_0
2309  v_com(:,:,:,2,2) = v_com(:,:,:,2,2)*jac
2310  v_com(:,:,:,3,2) = v_com(:,:,:,3,2)*jac
2311  deallocate(jac)
2312 #if ldebug
2313  end if
2314 #endif
2315 
2316  ! integrate toroidal flux
2317  if (present(v_flux_tor) .and. grid_trim%n(1).gt.1 .and. &
2318  &.not.compare_tor_pos_loc) then ! integrate poloidally and normally
2319  ! loop over all toroidal points
2320  do jd = 1,grid_trim%n(2)
2321 #if ldebug
2322  if (debug_calc_vec_comp) then
2323  ! for all normal points, integrate poloidally the
2324  ! covariant poloidal quantity
2325  do kd = norm_id(1),norm_id(2)
2326  ierr = calc_int(v_com(:,jd,kd,2,1),&
2327  &grid%theta_F(:,jd,kd),v_com(:,jd,kd,2,2)) ! save in contravariant variable
2328  chckerr('')
2329  end do
2330 
2331  ! gather on master
2332  ierr = get_ser_var(v_com(grid_trim%n(1),jd,&
2333  &norm_id(1):norm_id(2),2,2),v_ser_temp)
2334  chckerr('')
2335 
2336  ! update integrals if master
2337  if (rank.eq.0) then
2338  if (use_pol_flux_f) then
2339  v_flux_tor(:,jd) = v_flux_tor(:,jd) + v_ser_temp
2340  else
2341  v_flux_tor(:,jd+plot_offset(1)) = v_ser_temp
2342  end if
2343  end if
2344  else
2345 #endif
2346  ! for all normal points, integrate poloidally
2347  do kd = norm_id(1),norm_id(2)
2348  ierr = calc_int(v_com(:,jd,kd,3,2),&
2349  &grid%theta_F(:,jd,kd),v_com(:,jd,kd,3,1)) ! save in covariant variable
2350  chckerr('')
2351  end do
2352 
2353  ! gather on master
2354  ierr = get_ser_var(v_com(grid_trim%n(1),jd,&
2355  &norm_id(1):norm_id(2),3,1),v_ser_temp)
2356  chckerr('')
2357 
2358  ! integrate normally and update integrals if master
2359  if (rank.eq.0) then
2360  ierr = calc_int(v_ser_temp,&
2361  &grid_trim%r_F(norm_id(1):),v_ser_temp_int)
2362  chckerr('')
2363  if (use_normalization) v_ser_temp_int = &
2364  &v_ser_temp_int*psi_0
2365  if (use_pol_flux_f) then
2366  v_flux_tor(:,jd) = v_flux_tor(:,jd) + &
2367  &v_ser_temp_int
2368  else
2369  v_flux_tor(:,jd+plot_offset(1)) = v_ser_temp_int
2370  end if
2371  end if
2372 #if ldebug
2373  end if
2374 #endif
2375  end do
2376 
2377  ! plot output
2378  if (rank.eq.0 .and. do_plot .and. &
2379  &eq_job_nr.eq.size(eq_jobs_lims,2)) then
2380  call print_ex_2d([var_names(2,2)],&
2381  &trim(file_names(2)),v_flux_tor,&
2382  &x=reshape(grid_trim%r_F(norm_id(1):)*2*pi/max_flux_f,&
2383  &[grid_trim%n(3),1]),draw=.false.)
2384  call draw_ex([var_names(2,2)],trim(file_names(2)),&
2385  &plot_dim(2),1,.false.)
2386  end if
2387  end if
2388 
2389  ! integrate poloidal flux
2390  if (present(v_flux_pol) .and. grid_trim%n(2).gt.1 .and. &
2391  &.not.compare_tor_pos_loc) then ! integrate toroidally and normally
2392  ! loop over all poloidal points
2393  do id = 1,grid_trim%n(1)
2394 #if ldebug
2395  if (debug_calc_vec_comp) then
2396  ! for all normal points, integrate toroidally the
2397  ! covariant toroidal quantity
2398  do kd = norm_id(1),norm_id(2)
2399  ierr = calc_int(v_com(id,:,kd,3,1),&
2400  &grid%zeta_F(id,:,kd),v_com(id,:,kd,3,2)) ! save in contravariant variable
2401  chckerr('')
2402  end do
2403 
2404  ! gather on master
2405  ierr = get_ser_var(v_com(id,grid_trim%n(2),&
2406  &norm_id(1):norm_id(2),3,2),v_ser_temp)
2407  chckerr('')
2408 
2409  ! update integrals if master
2410  if (rank.eq.0) then
2411  if (use_pol_flux_f) then
2412  v_flux_pol(:,id+plot_offset(1)) = v_ser_temp
2413  else
2414  v_flux_pol(:,id) = v_flux_pol(:,id) + v_ser_temp
2415  end if
2416  end if
2417  else
2418 #endif
2419  ! for all normal points, integrate toroidally
2420  do kd = norm_id(1),norm_id(2)
2421  ierr = calc_int(v_com(id,:,kd,2,2),&
2422  &grid%zeta_F(id,:,kd),v_com(id,:,kd,2,1)) ! save in covariant variable
2423  chckerr('')
2424  end do
2425 
2426  ! gather on master
2427  ierr = get_ser_var(v_com(id,grid_trim%n(2),&
2428  &norm_id(1):norm_id(2),2,1),v_ser_temp)
2429  chckerr('')
2430 
2431  ! integrate normally and update integrals if master
2432  if (rank.eq.0) then
2433  ierr = calc_int(v_ser_temp,&
2434  &grid_trim%r_F(norm_id(1):),v_ser_temp_int)
2435  chckerr('')
2436  if (use_normalization) v_ser_temp_int = &
2437  &v_ser_temp_int*psi_0
2438  if (use_pol_flux_f) then
2439  v_flux_pol(:,id+plot_offset(1)) = v_ser_temp_int
2440  else
2441  v_flux_pol(:,id) = v_flux_pol(:,id) + &
2442  &v_ser_temp_int
2443  end if
2444  end if
2445 #if ldebug
2446  end if
2447 #endif
2448  end do
2449 
2450  ! plot output
2451  if (rank.eq.0 .and. do_plot .and. &
2452  &eq_job_nr.eq.size(eq_jobs_lims,2)) then
2453  call print_ex_2d([var_names(1,2)],&
2454  &trim(file_names(1)),v_flux_pol,&
2455  &x=reshape(grid_trim%r_F(norm_id(1):)*2*pi/max_flux_f,&
2456  &[grid_trim%n(3),1]),draw=.false.)
2457  call draw_ex([var_names(1,2)],trim(file_names(1)),&
2458  &plot_dim(1),1,.false.)
2459  end if
2460  end if
2461 
2462  ! clean up
2463  deallocate(v_ser_temp)
2464  if (rank.eq.0) deallocate(v_ser_temp_int)
2465  end if
2466 
2467  call lvl_ud(-1)
2468 
2469  if (max_transf_loc.eq.2) return
2470 
2471  ! 3. HELENA coordinates (psi,theta,zeta) or VMEC coordinates
2472  ! (r,theta,phi), by converting using transformation matrices
2473  ! (v_i)_H = T_H^F (v_i)_F
2474  ! (v^i)_H = (v^i)_F T_F^H
2475  ! Note: This is done directly from step 1; The results from step 2 are
2476  ! skipped, i.e. v_temp is not overwritten after step 2.
2477  call writo('Equilibrium coordinates')
2478  call lvl_ud(1)
2479  do ld = 1,9
2480  do jd = 1,grid_eq%n(2)
2481  do id = 1,grid_eq%n(1)
2482  ierr = spline(grid_eq%loc_r_F,&
2483  &eq_2%T_FE(id,jd,:,ld,0,0,0),grid%loc_r_F,&
2484  &t_ab(id,jd,:,ld,0,0,0),ord=norm_disc_prec)
2485  chckerr('')
2486  ierr = spline(grid_eq%loc_r_F,&
2487  &eq_2%T_EF(id,jd,:,ld,0,0,0),grid%loc_r_F,&
2488  &t_ba(id,jd,:,ld,0,0,0),ord=norm_disc_prec)
2489  chckerr('')
2490  end do
2491  end do
2492  end do
2493  v_com = 0._dp
2494  do jd = 1,3
2495  do id = 1,3
2496  v_com(:,:,:,jd,1) = v_com(:,:,:,jd,1) + &
2497  &t_ba(:,:,:,c([jd,id],.false.),0,0,0) * &
2498  &v_temp(:,:,:,id,1)
2499  v_com(:,:,:,jd,2) = v_com(:,:,:,jd,2) + &
2500  &t_ab(:,:,:,c([id,jd],.false.),0,0,0) * &
2501  &v_temp(:,:,:,id,2)
2502  end do
2503  end do
2504  if (present(v_mag)) then
2505  v_mag = 0._dp
2506  do id = 1,3
2507  v_mag = v_mag + v_com(:,:,:,id,1)*v_com(:,:,:,id,2)
2508  end do
2509  v_mag = sqrt(v_mag)
2510  end if
2511 
2512  ! save temporary copy, normalized
2513  v_temp = v_com
2514 
2515  ! set up plot variables
2516  if (do_plot) then
2517  select case (eq_style)
2518  case (1) ! VMEC
2519  coord_names(1) = 'r'
2520  coord_names(2) = 'theta'
2521  coord_names(3) = 'phi'
2522  description(1) = 'covariant components of the magnetic &
2523  &field in VMEC coordinates'
2524  description(2) = 'contravariant components of the magnetic &
2525  &field in VMEC coordinates'
2526  description(3) = 'magnitude of the magnetic field in VMEC &
2527  &coordinates'
2528  file_names(1) = trim(base_name)//'_V_sub'
2529  file_names(2) = trim(base_name)//'_V_sup'
2530  file_names(3) = trim(base_name)//'_V_mag'
2531  case (2) ! HELENA
2532  coord_names(1) = 'psi'
2533  coord_names(2) = 'theta'
2534  coord_names(3) = 'phi'
2535  description(1) = 'covariant components of the magnetic &
2536  &field in HELENA coordinates'
2537  description(2) = 'contravariant components of the magnetic &
2538  &field in HELENA coordinates'
2539  description(3) = 'magnitude of the magnetic field in &
2540  &HELENA coordinates'
2541  file_names(1) = trim(base_name)//'_H_sub'
2542  file_names(2) = trim(base_name)//'_H_sup'
2543  file_names(3) = trim(base_name)//'_H_mag'
2544  end select
2545  if (compare_tor_pos_loc) then
2546  do id = 1,3
2547  file_names(id) = trim(file_names(id))//'_COMP'
2548  end do
2549  end if
2550  var_names = trim(base_name)
2551  do id = 1,3
2552  var_names(id,1) = trim(var_names(id,1))//'_sub_'//&
2553  &trim(coord_names(id))
2554  var_names(id,2) = trim(var_names(id,2))//'_sup_'//&
2555  &trim(coord_names(id))
2556  end do
2557  if (use_normalization) then
2558  v_com(:,:,:,1,1) = v_com(:,:,:,1,1) * r_0 ! norm factor for e_r or e_psi
2559  v_com(:,:,:,2,1) = v_com(:,:,:,2,1) * r_0 ! norm factor for e_theta
2560  v_com(:,:,:,3,1) = v_com(:,:,:,3,1) * r_0 ! norm factor for e_zeta or e_phi
2561  v_com(:,:,:,1,2) = v_com(:,:,:,1,2) / r_0 ! norm factor for e^r or e^psi
2562  v_com(:,:,:,2,2) = v_com(:,:,:,2,2) / r_0 ! norm factor for e^theta
2563  v_com(:,:,:,3,2) = v_com(:,:,:,3,2) / r_0 ! norm factor for e^zeta or e^phi
2564  end if
2565  if (compare_tor_pos_loc) then
2566  ierr = calc_tor_diff(v_com,theta_geo,norm_disc_prec)
2567  chckerr('')
2568  end if
2569  do id = 1,2
2570  call plot_hdf5(var_names(:,id),trim(file_names(id)),&
2571  &v_com(:,tor_id(1):tor_id(2),norm_id(1):norm_id(2),:,id),&
2572  &tot_dim=plot_dim,loc_offset=plot_offset,&
2573  &x=x(:,tor_id(1):tor_id(2),:,:),&
2574  &y=y(:,tor_id(1):tor_id(2),:,:),&
2575  &z=z(:,tor_id(1):tor_id(2),:,:),&
2576  &cont_plot=cont_plot,descr=description(id))
2577  end do
2578  if (present(v_mag)) then
2579  if (compare_tor_pos_loc) then
2580  ierr = calc_tor_diff(v_mag,theta_geo,norm_disc_prec)
2581  chckerr('')
2582  end if
2583  call plot_hdf5(trim(base_name),trim(file_names(3)),&
2584  &v_mag(:,tor_id(1):tor_id(2),norm_id(1):norm_id(2)),&
2585  &tot_dim=plot_dim(1:3),loc_offset=plot_offset(1:3),&
2586  &x=x(:,tor_id(1):tor_id(2),:,1),&
2587  &y=y(:,tor_id(1):tor_id(2),:,1),&
2588  &z=z(:,tor_id(1):tor_id(2),:,1),&
2589  &cont_plot=cont_plot,descr=description(3))
2590  end if
2591  end if
2592 
2593  call lvl_ud(-1)
2594 
2595  if (max_transf_loc.eq.3) return
2596 
2597  ! 4. cylindrical coordinates (R,phi,Z) by converting using
2598  ! transformation matrices:
2599  ! (v_i)_C = T_C^E (v_i)_E
2600  ! (v^i)_C = (v^i)_E T_E^C
2601  ! where for VMEC (E=V), the transformation matrix T_VC is tabulated, and
2602  ! its inverse is calculate, and for HELENA (E=H), the transformation
2603  ! matrix T_HC is given by
2604  ! ( R' 0 Z' )
2605  ! T_HC = ( R_theta 0 Z_theta )
2606  ! ( 0 -1 0 )
2607  ! and its inverse can be calculated as well.
2608  call writo('Cylindrical coordinates')
2609  call lvl_ud(1)
2610  t_ba = 0._dp
2611  t_ab = 0._dp
2612  select case (eq_style)
2613  case (1) ! VMEC
2614  do ld = 1,9
2615  do jd = 1,grid_eq%n(2)
2616  do id = 1,grid_eq%n(1)
2617  ierr = spline(grid_eq%loc_r_F,&
2618  &eq_2%T_VC(id,jd,:,ld,0,0,0),grid%loc_r_F,&
2619  &t_ab(id,jd,:,ld,0,0,0),ord=norm_disc_prec)
2620  chckerr('')
2621  end do
2622  end do
2623  end do
2624  case (2) ! HELENA
2625  ! setup D1R and D2R, and same thing for Z
2626  ! (use untrimmed grid, with ghost region)
2627  allocate(d1r(grid%n(1),grid%n(2),grid%loc_n_r))
2628  allocate(d2r(grid%n(1),grid%n(2),grid%loc_n_r))
2629  allocate(d1z(grid%n(1),grid%n(2),grid%loc_n_r))
2630  allocate(d2z(grid%n(1),grid%n(2),grid%loc_n_r))
2631  do jd = 1,grid%n(2)
2632  do id = 1,grid%n(1)
2633  ierr = spline(grid%loc_r_E,xyzr(id,jd,:,4)/norm_len,&
2634  &grid%loc_r_E,d1r(id,jd,:),ord=norm_disc_prec,&
2635  &deriv=1)
2636  chckerr('')
2637  ierr = spline(grid%loc_r_E,xyzr(id,jd,:,3)/norm_len,&
2638  &grid%loc_r_E,d1z(id,jd,:),ord=norm_disc_prec,&
2639  &deriv=1)
2640  chckerr('')
2641  end do
2642  end do
2643  do kd = 1,grid%loc_n_r
2644  do jd = 1,grid%n(2)
2645  ierr = spline(grid%theta_E(:,jd,kd),&
2646  &xyzr(:,jd,kd,4)/norm_len,grid%theta_E(:,jd,kd),&
2647  &d2r(:,jd,kd),ord=norm_disc_prec,deriv=1)
2648  chckerr('')
2649  ierr = spline(grid%theta_E(:,jd,kd),&
2650  &xyzr(:,jd,kd,3)/norm_len,grid%theta_E(:,jd,kd),&
2651  &d2z(:,jd,kd),ord=norm_disc_prec,deriv=1)
2652  chckerr('')
2653  end do
2654  end do
2655  t_ab(:,:,:,c([1,1],.false.),0,0,0) = d1r
2656  t_ab(:,:,:,c([1,3],.false.),0,0,0) = d1z
2657  t_ab(:,:,:,c([2,1],.false.),0,0,0) = d2r
2658  t_ab(:,:,:,c([2,3],.false.),0,0,0) = d2z
2659  t_ab(:,:,:,c([3,2],.false.),0,0,0) = -1._dp
2660  deallocate(d1r,d2r,d1z,d2z)
2661  end select
2662  ierr = calc_inv_met(t_ba,t_ab,[0,0,0])
2663  chckerr('')
2664  v_com = 0._dp
2665  do jd = 1,3
2666  do id = 1,3
2667  v_com(:,:,:,jd,1) = v_com(:,:,:,jd,1) + t_ba(:,:,:,&
2668  &c([jd,id],.false.),0,0,0) * v_temp(:,:,:,id,1)
2669  v_com(:,:,:,jd,2) = v_com(:,:,:,jd,2) + t_ab(:,:,:,&
2670  &c([id,jd],.false.),0,0,0) * v_temp(:,:,:,id,2)
2671  end do
2672  end do
2673  if (present(v_mag)) then
2674  v_mag = 0._dp
2675  do id = 1,3
2676  v_mag = v_mag + v_com(:,:,:,id,1)*v_com(:,:,:,id,2)
2677  end do
2678  v_mag = sqrt(v_mag)
2679  end if
2680 
2681  ! save temporary copy, normalized
2682  v_temp = v_com
2683 
2684  ! set up plot variables
2685  if (do_plot) then
2686  description(1) = 'covariant components of the magnetic field in &
2687  &cylindrical coordinates'
2688  description(2) = 'contravariant components of the magnetic field &
2689  &in cylindrical coordinates'
2690  description(3) = 'magnitude of the magnetic field in cylindrical &
2691  &coordinates'
2692  file_names(1) = trim(base_name)//'_C_sub'
2693  file_names(2) = trim(base_name)//'_C_sup'
2694  file_names(3) = trim(base_name)//'_C_mag'
2695  if (compare_tor_pos_loc) then
2696  do id = 1,3
2697  file_names(id) = trim(file_names(id))//'_COMP'
2698  end do
2699  end if
2700  coord_names(1) = 'R'
2701  coord_names(2) = 'phi'
2702  coord_names(3) = 'Z'
2703  var_names = trim(base_name)
2704  do id = 1,3
2705  var_names(id,1) = trim(var_names(id,1))//'_sub_'//&
2706  &trim(coord_names(id))
2707  var_names(id,2) = trim(var_names(id,2))//'_sup_'//&
2708  &trim(coord_names(id))
2709  end do
2710  if (use_normalization) then
2711  !v_com(:,:,:,1,1) = v_com(:,:,:,1,1) ! norm factor for e_R
2712  v_com(:,:,:,2,1) = v_com(:,:,:,2,1) * r_0 ! norm factor for e_phi
2713  !v_com(:,:,:,3,1) = v_com(:,:,:,3,1) ! norm factor for e_Z
2714  !v_com(:,:,:,1,2) = v_com(:,:,:,1,2) ! norm factor for e^R
2715  v_com(:,:,:,2,2) = v_com(:,:,:,2,2) / r_0 ! norm factor for e^phi
2716  !v_com(:,:,:,3,2) = v_com(:,:,:,3,2) ! norm factor for e^Z
2717  end if
2718  if (compare_tor_pos_loc) then
2719  ierr = calc_tor_diff(v_com,theta_geo,norm_disc_prec)
2720  chckerr('')
2721  end if
2722  do id = 1,2
2723  call plot_hdf5(var_names(:,id),trim(file_names(id)),&
2724  &v_com(:,tor_id(1):tor_id(2),norm_id(1):norm_id(2),:,id),&
2725  &tot_dim=plot_dim,loc_offset=plot_offset,&
2726  &x=x(:,tor_id(1):tor_id(2),:,:),&
2727  &y=y(:,tor_id(1):tor_id(2),:,:),&
2728  &z=z(:,tor_id(1):tor_id(2),:,:),&
2729  &cont_plot=cont_plot,descr=description(id))
2730  end do
2731  if (present(v_mag)) then
2732  if (compare_tor_pos_loc) then
2733  ierr = calc_tor_diff(v_mag,theta_geo,norm_disc_prec)
2734  chckerr('')
2735  end if
2736  call plot_hdf5(trim(base_name),trim(file_names(3)),&
2737  &v_mag(:,tor_id(1):tor_id(2),norm_id(1):norm_id(2)),&
2738  &tot_dim=plot_dim(1:3),loc_offset=plot_offset(1:3),&
2739  &x=x(:,tor_id(1):tor_id(2),:,1),&
2740  &y=y(:,tor_id(1):tor_id(2),:,1),&
2741  &z=z(:,tor_id(1):tor_id(2),:,1),&
2742  &cont_plot=cont_plot,descr=description(3))
2743  end if
2744  end if
2745 
2746  call lvl_ud(-1)
2747 
2748  if (max_transf_loc.eq.4) return
2749 
2750  ! 5. Cartesian coordinates (X,Y,Z) by converting using transformation
2751  ! matrices
2752  ! (e_R) ( cos(phi) sin(phi) 0 ) (e_X)
2753  ! (e_phi) = ( -R sin(phi) R cos(phi) 0 ) (e_Y)
2754  ! (e_Z) ( 0 0 1 ) (e_Z)
2755  ! with phi = - zeta_F. For the transformation of contravariant
2756  ! components, the factor R becomes 1/R and the transpose is taken.
2757  call writo('Cartesian coordinates')
2758  call lvl_ud(1)
2759  t_ba = 0._dp
2760  t_ab = 0._dp
2761  t_ab(:,:,:,c([1,1],.false.),0,0,0) = cos(-grid%zeta_F)
2762  t_ab(:,:,:,c([1,2],.false.),0,0,0) = sin(-grid%zeta_F)
2763  t_ab(:,:,:,c([2,1],.false.),0,0,0) = -xyzr(:,:,:,4)/norm_len*&
2764  &sin(-grid%zeta_F)
2765  t_ab(:,:,:,c([2,2],.false.),0,0,0) = xyzr(:,:,:,4)/norm_len*&
2766  &cos(-grid%zeta_F)
2767  t_ab(:,:,:,c([3,3],.false.),0,0,0) = 1._dp
2768  ierr = calc_inv_met(t_ba,t_ab,[0,0,0])
2769  chckerr('')
2770  v_com = 0._dp
2771  do jd = 1,3
2772  do id = 1,3
2773  v_com(:,:,:,jd,1) = v_com(:,:,:,jd,1) + t_ba(:,:,:,&
2774  &c([jd,id],.false.),0,0,0) * v_temp(:,:,:,id,1)
2775  v_com(:,:,:,jd,2) = v_com(:,:,:,jd,2) + t_ab(:,:,:,&
2776  &c([id,jd],.false.),0,0,0) * v_temp(:,:,:,id,2)
2777  end do
2778  end do
2779  if (present(v_mag)) then
2780  v_mag = 0._dp
2781  do id = 1,3
2782  v_mag = v_mag + v_com(:,:,:,id,1)*v_com(:,:,:,id,2)
2783  end do
2784  v_mag = sqrt(v_mag)
2785  end if
2786 
2787  ! set up plot variables
2788  if (do_plot) then
2789  description(1) = 'covariant components of the magnetic field in &
2790  &cartesian coordinates'
2791  description(2) = 'contravariant components of the magnetic field &
2792  &in cartesian coordinates'
2793  description(3) = 'magnitude of the magnetic field in cartesian &
2794  &coordinates'
2795  file_names(1) = trim(base_name)//'_X_sub'
2796  file_names(2) = trim(base_name)//'_X_sup'
2797  file_names(3) = trim(base_name)//'_X_mag'
2798  if (compare_tor_pos_loc) then
2799  do id = 1,3
2800  file_names(id) = trim(file_names(id))//'_COMP'
2801  end do
2802  end if
2803  coord_names(1) = 'X'
2804  coord_names(2) = 'Y'
2805  coord_names(3) = 'Z'
2806  var_names = trim(base_name)
2807  do id = 1,3
2808  var_names(id,1) = trim(var_names(id,1))//'_sub_'//&
2809  &trim(coord_names(id))
2810  var_names(id,2) = trim(var_names(id,2))//'_sup_'//&
2811  &trim(coord_names(id))
2812  end do
2813  if (compare_tor_pos_loc) then
2814  ierr = calc_tor_diff(v_com,theta_geo,norm_disc_prec)
2815  chckerr('')
2816  end if
2817  do id = 1,2
2818  call plot_hdf5(var_names(:,id),trim(file_names(id)),&
2819  &v_com(:,tor_id(1):tor_id(2),norm_id(1):norm_id(2),:,id),&
2820  &tot_dim=plot_dim,loc_offset=plot_offset,&
2821  &x=x(:,tor_id(1):tor_id(2),:,:),&
2822  &y=y(:,tor_id(1):tor_id(2),:,:),&
2823  &z=z(:,tor_id(1):tor_id(2),:,:),&
2824  &cont_plot=cont_plot,descr=description(id))
2825  end do
2826  if (present(v_mag)) then
2827  if (compare_tor_pos_loc) then
2828  ierr = calc_tor_diff(v_mag,theta_geo,norm_disc_prec)
2829  chckerr('')
2830  end if
2831  call plot_hdf5(trim(base_name),trim(file_names(3)),&
2832  &v_mag(:,tor_id(1):tor_id(2),norm_id(1):norm_id(2)),&
2833  &tot_dim=plot_dim(1:3),loc_offset=plot_offset(1:3),&
2834  &x=x(:,tor_id(1):tor_id(2),:,1),&
2835  &y=y(:,tor_id(1):tor_id(2),:,1),&
2836  &z=z(:,tor_id(1):tor_id(2),:,1),&
2837  &cont_plot=cont_plot,descr=description(3))
2838  end if
2839 
2840  ! plot vector as well
2841  call plot_hdf5([trim(base_name)],trim(base_name)//'_vec',&
2842  &v_com(:,tor_id(1):tor_id(2),norm_id(1):norm_id(2),:,1),&
2843  &tot_dim=plot_dim,loc_offset=plot_offset,&
2844  &x=x(:,tor_id(1):tor_id(2),:,:),y=y(:,tor_id(1):tor_id(2),:,:),&
2845  &z=z(:,tor_id(1):tor_id(2),:,:),col=4,cont_plot=cont_plot,&
2846  &descr='magnetic field vector')
2847  end if
2848 
2849  call lvl_ud(-1)
2850 
2851  ! clean up
2852  call grid_trim%dealloc()
2853  end function calc_vec_comp
2854 
2861  integer function calc_n_par_x_rich(n_par_X_rich,only_half_grid) result(ierr)
2862  use rich_vars, only: n_par_x
2863 
2864  character(*), parameter :: rout_name = 'calc_n_par_X_rich'
2865 
2866  ! input / output
2867  integer, intent(inout) :: n_par_x_rich
2868  logical, intent(in), optional :: only_half_grid
2869 
2870  ! local variables
2871  character(len=max_str_ln) :: err_msg ! error message
2872 
2873  ! initialize ierr
2874  ierr = 0
2875 
2876  ! possibly divide n_par_X by 2
2877  n_par_x_rich = n_par_x
2878  if (present(only_half_grid)) then
2879  if (only_half_grid) then
2880  if (mod(n_par_x,2).eq.1) then
2881  n_par_x_rich = (n_par_x-1)/2
2882  else
2883  ierr = 1
2884  err_msg = 'Need odd number of points'
2885  chckerr(err_msg)
2886  end if
2887  end if
2888  end if
2889  end function calc_n_par_x_rich
2890 
2900  integer function nufft(x,f,f_F,plot_name) result(ierr)
2901  use num_utilities, only: order_per_fun, spline
2902 
2903  character(*), parameter :: rout_name = 'nufft'
2904 
2905  ! input / output
2906  real(dp), intent(in) :: x(:)
2907  real(dp), intent(in) :: f(:)
2908  real(dp), intent(inout), allocatable :: f_f(:,:)
2909  character(len=*), intent(in), optional :: plot_name
2910 
2911  ! local variables
2912  integer :: m_f ! nr. of modes
2913  integer :: n_x ! nr. of points
2914  integer :: id ! counter
2915  logical :: print_log = .false. ! print log plot as well
2916  real(dp), allocatable :: x_loc(:) ! local x, extended
2917  real(dp), allocatable :: work(:) ! work array
2918  real(dp), allocatable :: f_loc(:) ! local f, extended
2919  real(dp), allocatable :: f_int(:) ! interpolated f, later Fourier modes
2920  character(len=max_str_ln) :: plot_title(2) ! name of plot
2921 
2922  ! tests
2923  if (size(x).ne.size(f)) then
2924  ierr = 1
2925  chckerr('x and f must have same size')
2926  end if
2927  if (maxval(x)-minval(x).gt.2*pi) then
2928  ierr = 1
2929  chckerr('x is not a fundamental interval')
2930  end if
2931 
2932  ! create local x and f that go from <0 to <2pi, with a royal overlap for
2933  ! interpolation
2934  ierr = order_per_fun(x,f,x_loc,f_loc,10)
2935  chckerr('')
2936 
2937  ! set local f and interpolate
2938  n_x = size(x)
2939  allocate(f_int(n_x))
2940  ierr = spline(x_loc,f_loc,[((id-1._dp)/n_x*2*pi,id=1,n_x)],f_int,&
2941  &ord=3)
2942  chckerr('')
2943 
2944  ! clean up
2945  deallocate(x_loc)
2946  deallocate(f_loc)
2947 
2948  !!do id = 1,n_x
2949  !!f_int(id) = -4._dp+&
2950  !!&3*cos((id-1._dp)/n_x*2*pi*1)+&
2951  !!&0.5*cos((id-1._dp)/n_x*2*pi*100)+&
2952  !!&1.5*cos((id-1._dp)/n_x*2*pi*101)+&
2953  !!&4*sin((id-1._dp)/n_x*2*pi*1)+&
2954  !!&2*sin((id-1._dp)/n_x*2*pi*50)+&
2955  !!&4.5*sin((id-1._dp)/n_x*2*pi*55)
2956  !!end do
2957 
2958  ! set up variables for fft
2959  m_f = (n_x-1)/2 ! nr. of modes
2960  allocate(work(3*n_x+15)) ! from fftpack manual
2961 
2962  ! calculate fft
2963  call dffti(n_x,work)
2964  call dfftf(n_x,f_int,work)
2965  deallocate(work)
2966 
2967  ! rescale
2968  f_int(:) = f_int(:)*2/n_x
2969  f_int(1) = f_int(1)/2
2970 
2971  ! separate cos and sine
2972  if (allocated(f_f)) deallocate(f_f)
2973  allocate(f_f(m_f+1,2))
2974  f_f(1,1) = f_int(1)
2975  f_f(1,2) = 0._dp
2976  f_f(2:m_f+1,1) = f_int(2:2*m_f:2)
2977  f_f(2:m_f+1,2) = -f_int(3:2*m_f+1:2) ! routine returns - sine factors
2978 
2979  ! output in plot if requested
2980  if (present(plot_name)) then
2981  plot_title = ['cos','sin']
2982  call print_ex_2d(plot_title,plot_name,f_f,draw=.false.)
2983  call draw_ex(plot_title,plot_name,2,1,.false.)
2984  if (print_log) then
2985  plot_title = ['cos [log]','sin [log]']
2986  call print_ex_2d(plot_title,trim(plot_name)//'_log',&
2987  &log10(abs(f_f)),draw=.false.)
2988  call draw_ex(plot_title,trim(plot_name)//'_log',2,1,.false.)
2989  end if
2990  end if
2991  end function nufft
2992 
2994  subroutine find_compr_range(r_F,lim_r,lim_id)
2995  ! input / output
2996  real(dp), intent(in) :: r_f(:)
2997  real(dp), intent(in) :: lim_r(2)
2998  integer, intent(inout) :: lim_id(2)
2999 
3000  ! local variables
3001  integer :: id ! counter
3002  integer :: n_r ! number of points in coordinate
3003  real(dp), parameter :: tol = 1.e-6 ! tolerance for grids
3004 
3005  ! set n_r
3006  n_r = size(r_f)
3007 
3008  lim_id = [1,n_r] ! initialize with full range
3009  if (r_f(1).lt.r_f(n_r)) then ! ascending r_F
3010  do id = 1,n_r
3011  if (r_f(id).le.minval(lim_r)-tol) lim_id(1) = id ! move lower limit up
3012  if (r_f(n_r-id+1).ge.maxval(lim_r)+tol) &
3013  &lim_id(2) = n_r-id+1 ! move upper limit down
3014  end do
3015  else ! descending r_F
3016  do id = 1,n_r
3017  if (r_f(id).ge.maxval(lim_r)+tol) lim_id(1) = id ! move lower limit up
3018  if (r_f(n_r-id+1).le.minval(lim_r)-tol) &
3019  &lim_id(2) = n_r-id+1 ! move upper limit down
3020  end do
3021  end if
3022  end subroutine find_compr_range
3023 
3038  integer function calc_arc_angle(grid,eq_1,eq_2,arc,use_E) result(ierr)
3039  use eq_vars, only: eq_1_type, eq_2_type, &
3040  &r_0
3041  use num_vars, only: use_normalization, use_pol_flux_f
3042  use num_utilities, only: c, calc_int
3043 
3044  character(*), parameter :: rout_name = 'calc_arc_angle'
3045 
3046  ! input / output
3047  type(grid_type), intent(in) :: grid
3048  type(eq_1_type), intent(in) :: eq_1
3049  type(eq_2_type), intent(in) :: eq_2
3050  real(dp), intent(inout), allocatable :: arc(:,:,:)
3051  logical, intent(in), optional :: use_e
3052 
3053  ! local variables
3054  integer :: jd, kd ! counters
3055  logical :: use_e_loc ! local use_E
3056 
3057  ! initialize ierr
3058  ierr = 0
3059 
3060  ! set local use_E
3061  use_e_loc = .false.
3062  if (present(use_e)) use_e_loc = use_e
3063 
3064  ! calculate arclength: int dtheta |e_theta|
3065  allocate(arc(grid%n(1),grid%n(2),grid%loc_n_r))
3066  do kd = 1,grid%loc_n_r
3067  do jd = 1,grid%n(2)
3068  if (use_e_loc) then
3069  ierr = calc_int(&
3070  &eq_2%g_E(:,jd,kd,c([2,2],.true.),0,0,0)**(0.5),&
3071  &grid%theta_E(:,jd,kd),arc(:,jd,kd))
3072  chckerr('')
3073  if (use_normalization) arc(:,jd,kd) = arc(:,jd,kd) * r_0 ! not really necessary as we take fractional arc length
3074  arc(:,jd,kd) = grid%theta_E(1,jd,kd) + &
3075  &arc(:,jd,kd) / arc(grid%n(1),jd,kd) * &
3076  &(grid%theta_E(grid%n(1),jd,kd)-grid%theta_E(1,jd,kd))
3077  else
3078  if (use_pol_flux_f) then
3079  ! e_arc = e_theta - q e_alpha
3080  ierr = calc_int((&
3081  &eq_2%g_FD(:,jd,kd,c([3,3],.true.),0,0,0) - &
3082  &eq_2%g_FD(:,jd,kd,c([1,3],.true.),0,0,0) * &
3083  &2*eq_1%q_saf_FD(kd,0) + &
3084  &eq_2%g_FD(:,jd,kd,c([1,1],.true.),0,0,0) * &
3085  &eq_1%q_saf_FD(kd,0)**2)**(0.5),&
3086  &grid%theta_F(:,jd,kd),arc(:,jd,kd))
3087  chckerr('')
3088  else
3089  ! e_arc = -e_alpha
3090  ierr = calc_int(&
3091  &(-eq_2%g_FD(:,jd,kd,c([1,1],.true.),0,0,0))&
3092  &**(0.5),grid%theta_F(:,jd,kd),arc(:,jd,kd))
3093  chckerr('')
3094  end if
3095  if (use_normalization) arc(:,jd,kd) = arc(:,jd,kd) * r_0 ! not really necessary as we take fractional arc length
3096  arc(:,jd,kd) = grid%theta_F(1,jd,kd) + &
3097  &arc(:,jd,kd) / arc(grid%n(1),jd,kd) * &
3098  &(grid%theta_F(grid%n(1),jd,kd)-grid%theta_F(1,jd,kd))
3099  end if
3100  end do
3101  end do
3102  end function calc_arc_angle
3103 end module grid_utilities
num_utilities::c
integer function, public c(ij, sym, n, lim_n)
Convert 2-D coordinates (i,j) to the storage convention used in matrices.
Definition: num_utilities.f90:2556
mpi_utilities::get_ser_var
Gather parallel variable in serial version on group master.
Definition: MPI_utilities.f90:55
grid_utilities::debug_calc_int_vol
logical, public debug_calc_int_vol
plot debug information for calc_int_vol()
Definition: grid_utilities.f90:25
num_utilities::calc_int
Integrates a function using the trapezoidal rule.
Definition: num_utilities.f90:160
num_utilities::round_with_tol
Rounds an arry of values to limits, with a tolerance that can optionally be modified.
Definition: num_utilities.f90:169
num_vars::dp
integer, parameter, public dp
double precision
Definition: num_vars.f90:46
eq_vars
Variables that have to do with equilibrium quantities and the grid used in the calculations:
Definition: eq_vars.f90:27
num_vars::use_normalization
logical, public use_normalization
whether to use normalization or not
Definition: num_vars.f90:115
mpi_utilities
Numerical utilities related to MPI.
Definition: MPI_utilities.f90:20
num_vars
Numerical variables used by most other modules.
Definition: num_vars.f90:4
num_vars::max_str_ln
integer, parameter, public max_str_ln
maximum length of strings
Definition: num_vars.f90:50
mpi_utilities::get_ghost_arr
Fill the ghost regions in an array.
Definition: MPI_utilities.f90:73
rich_vars
Variables concerning Richardson extrapolation.
Definition: rich_vars.f90:4
num_vars::norm_disc_prec_eq
integer, public norm_disc_prec_eq
precision for normal discretization for equilibrium
Definition: num_vars.f90:120
num_ops
Numerical operations.
Definition: num_ops.f90:4
str_utilities::i2str
elemental character(len=max_str_ln) function, public i2str(k)
Convert an integer to string.
Definition: str_utilities.f90:18
grid_utilities::untrim_grid
integer function, public untrim_grid(grid_in, grid_out, size_ghost)
Untrims a trimmed grid by introducing an asymmetric ghost regions at the right.
Definition: grid_utilities.f90:1764
grid_utilities::calc_n_par_x_rich
integer function, public calc_n_par_x_rich(n_par_X_rich, only_half_grid)
Calculates the local number of parallel grid points for this Richardson level, taking into account th...
Definition: grid_utilities.f90:2862
num_vars::n_zeta_plot
integer, public n_zeta_plot
nr. of toroidal points in plot
Definition: num_vars.f90:157
grid_utilities::coord_f2e
Converts Flux coordinates to Equilibrium coordinates .
Definition: grid_utilities.f90:47
eq_utilities::calc_inv_met
Calculate from and where and , according to .
Definition: eq_utilities.f90:61
num_vars::eq_job_nr
integer, public eq_job_nr
nr. of eq job
Definition: num_vars.f90:79
num_vars::iu
complex(dp), parameter, public iu
complex unit
Definition: num_vars.f90:85
vmec_utilities::calc_trigon_factors
integer function, public calc_trigon_factors(theta, zeta, trigon_factors)
Calculate the trigonometric cosine and sine factors.
Definition: VMEC_utilities.f90:275
num_vars::max_theta_plot
real(dp), public max_theta_plot
max. of theta_plot
Definition: num_vars.f90:160
num_utilities::order_per_fun
Order a periodic function to include and an overlap.
Definition: num_utilities.f90:248
num_vars::n_procs
integer, public n_procs
nr. of MPI processes
Definition: num_vars.f90:69
output_ops::print_ex_2d
Print 2-D output on a file.
Definition: output_ops.f90:47
str_utilities
Operations on strings.
Definition: str_utilities.f90:4
grid_vars::grid_type
Type for grids.
Definition: grid_vars.f90:59
grid_utilities::calc_int_vol
integer function, public calc_int_vol(ang_1, ang_2, norm, J, f, f_int)
Calculates volume integral on a 3D grid.
Definition: grid_utilities.f90:1320
str_utilities::r2strt
elemental character(len=max_str_ln) function, public r2strt(k)
Convert a real (double) to string.
Definition: str_utilities.f90:54
eq_vars::r_0
real(dp), public r_0
independent normalization constant for nondimensionalization
Definition: eq_vars.f90:42
eq_vars::eq_1_type
flux equilibrium type
Definition: eq_vars.f90:63
eq_utilities
Numerical utilities related to equilibrium variables.
Definition: eq_utilities.f90:4
grid_utilities::debug_calc_vec_comp
logical, public debug_calc_vec_comp
plot debug information for calc_vec_comp()
Definition: grid_utilities.f90:27
num_vars::min_zeta_plot
real(dp), public min_zeta_plot
min. of zeta_plot
Definition: num_vars.f90:161
num_vars::compare_tor_pos
logical, public compare_tor_pos
compare quantities at toroidal positions (only for POST)
Definition: num_vars.f90:151
helena_vars::r_h
real(dp), dimension(:,:), allocatable, public r_h
major radius (xout)
Definition: HELENA_vars.f90:33
num_vars::min_theta_plot
real(dp), public min_theta_plot
min. of theta_plot
Definition: num_vars.f90:159
helena_vars::z_h
real(dp), dimension(:,:), allocatable, public z_h
height (yout)
Definition: HELENA_vars.f90:34
grid_utilities::extend_grid_f
integer function, public extend_grid_f(grid_in, grid_ext, grid_eq, n_theta_plot, n_zeta_plot, lim_theta_plot, lim_zeta_plot)
Extend a grid angularly.
Definition: grid_utilities.f90:1096
num_utilities::spline
Wrapper to the pspline library, making it easier to use for 1-D applications where speed is not the m...
Definition: num_utilities.f90:276
output_ops::draw_ex
subroutine, public draw_ex(var_names, draw_name, nplt, draw_dim, plot_on_screen, ex_plot_style, data_name, draw_ops, extra_ops, is_animated, ranges, delay, persistent)
Use external program to draw a plot.
Definition: output_ops.f90:1079
num_vars::rz_0
real(dp), dimension(2), public rz_0
origin of geometrical poloidal coordinate
Definition: num_vars.f90:179
output_ops::plot_hdf5
Prints variables vars with names var_names in an HDF5 file with name c file_name and accompanying XDM...
Definition: output_ops.f90:137
vmec_utilities
Numerical utilities related to the output of VMEC.
Definition: VMEC_utilities.f90:4
num_vars::eq_style
integer, public eq_style
either 1 (VMEC) or 2 (HELENA)
Definition: num_vars.f90:89
rich_vars::n_par_x
integer, public n_par_x
nr. of parallel points in field-aligned grid
Definition: rich_vars.f90:20
eq_vars::max_flux_f
real(dp), public max_flux_f
max. flux in Flux coordinates, set in calc_norm_range_PB3D_in
Definition: eq_vars.f90:50
num_vars::tol_zero
real(dp), public tol_zero
tolerance for zeros
Definition: num_vars.f90:133
num_ops::calc_zero_hh
Finds the zero of a function using Householder iteration.
Definition: num_ops.f90:35
grid_utilities::copy_grid
integer function, public copy_grid(grid_A, grid_B, lims_B, i_lim)
Copy a grid A to a new grid B, that was not yet initialized.
Definition: grid_utilities.f90:1189
grid_utilities::calc_eqd_grid
Calculate grid of equidistant points, where optionally the last point can be excluded.
Definition: grid_utilities.f90:75
num_utilities
Numerical utilities.
Definition: num_utilities.f90:4
str_utilities::r2str
elemental character(len=max_str_ln) function, public r2str(k)
Convert a real (double) to string.
Definition: str_utilities.f90:42
messages::writo
subroutine, public writo(input_str, persistent, error, warning, alert)
Write output to file identified by output_i.
Definition: messages.f90:275
messages
Numerical utilities related to giving output.
Definition: messages.f90:4
vmec_vars::z_v_s
real(dp), dimension(:,:,:), allocatable, public z_v_s
Coeff. of in cosine series (FM) and norm. deriv.
Definition: VMEC_vars.f90:42
helena_vars
Variables that have to do with HELENA quantities.
Definition: HELENA_vars.f90:4
num_vars::pi
real(dp), parameter, public pi
Definition: num_vars.f90:83
grid_utilities
Numerical utilities related to the grids and different coordinate systems.
Definition: grid_utilities.f90:4
grid_utilities::coord_e2f_r
integer function coord_e2f_r(grid_eq, r_E, r_F, r_E_array, r_F_array)
version with only r
Definition: grid_utilities.f90:473
vmec_vars::r_v_c
real(dp), dimension(:,:,:), allocatable, public r_v_c
Coeff. of in sine series (FM) and norm. deriv.
Definition: VMEC_vars.f90:39
grid_vars
Variables pertaining to the different grids used.
Definition: grid_vars.f90:4
grid_utilities::calc_vec_comp
integer function, public calc_vec_comp(grid, grid_eq, eq_1, eq_2, v_com, norm_disc_prec, v_mag, base_name, max_transf, v_flux_tor, v_flux_pol, XYZ, compare_tor_pos)
Calculates contra- and covariant components of a vector in multiple coordinate systems.
Definition: grid_utilities.f90:1859
helena_vars::ias
integer, public ias
0 if top-bottom symmetric, 1 if not
Definition: HELENA_vars.f90:36
messages::lvl_ud
subroutine, public lvl_ud(inc)
Increases/decreases lvl of output.
Definition: messages.f90:254
vmec_vars
Variables that concern the output of VMEC.
Definition: VMEC_vars.f90:4
vmec_vars::l_v_c
real(dp), dimension(:,:,:), allocatable, public l_v_c
Coeff. of in sine series (HM) and norm. deriv.
Definition: VMEC_vars.f90:43
num_vars::use_pol_flux_f
logical, public use_pol_flux_f
whether poloidal flux is used in F coords.
Definition: num_vars.f90:114
grid_utilities::find_compr_range
subroutine, public find_compr_range(r_F, lim_r, lim_id)
finds smallest range that comprises a minimum and maximum value.
Definition: grid_utilities.f90:2995
grid_utilities::coord_e2f_rtz
integer function coord_e2f_rtz(grid_eq, r_E, theta_E, zeta_E, r_F, theta_F, zeta_F, r_E_array, r_F_array)
version with r, theta and zeta
Definition: grid_utilities.f90:359
eq_vars::psi_0
real(dp), public psi_0
derived normalization constant for nondimensionalization
Definition: eq_vars.f90:46
vmec_vars::r_v_s
real(dp), dimension(:,:,:), allocatable, public r_v_s
Coeff. of in cosine series (FM) and norm. deriv.
Definition: VMEC_vars.f90:40
num_vars::eq_jobs_lims
integer, dimension(:,:), allocatable, public eq_jobs_lims
data about eq jobs: [ , ] for all jobs
Definition: num_vars.f90:77
helena_vars::chi_h
real(dp), dimension(:), allocatable, public chi_h
poloidal angle
Definition: HELENA_vars.f90:23
output_ops
Operations concerning giving output, on the screen as well as in output files.
Definition: output_ops.f90:5
num_vars::max_zeta_plot
real(dp), public max_zeta_plot
max. of zeta_plot
Definition: num_vars.f90:162
grid_utilities::calc_tor_diff
Calculates the toroidal difference for a magnitude calculated on three toroidal points: two extremiti...
Definition: grid_utilities.f90:109
num_vars::rank
integer, public rank
MPI rank.
Definition: num_vars.f90:68
grid_utilities::nufft
integer function, public nufft(x, f, f_F, plot_name)
calculates the cosine and sine mode numbers of a function defined on a non-regular grid.
Definition: grid_utilities.f90:2901
str_utilities::c2str
elemental character(len=max_str_ln) function, public c2str(k)
Convert a complex (double) to string.
Definition: str_utilities.f90:66
vmec_utilities::fourier2real
Inverse Fourier transformation, from VMEC.
Definition: VMEC_utilities.f90:56
vmec_vars::l_v_s
real(dp), dimension(:,:,:), allocatable, public l_v_s
Coeff. of in cosine series (HM) and norm. deriv.
Definition: VMEC_vars.f90:44
grid_utilities::calc_arc_angle
integer function, public calc_arc_angle(grid, eq_1, eq_2, arc, use_E)
Calculate arclength angle.
Definition: grid_utilities.f90:3039
eq_vars::eq_2_type
metric equilibrium type
Definition: eq_vars.f90:114
helena_vars::nchi
integer, public nchi
nr. of poloidal points
Definition: HELENA_vars.f90:35
vmec_vars::z_v_c
real(dp), dimension(:,:,:), allocatable, public z_v_c
Coeff. of in sine series (FM) and norm. deriv.
Definition: VMEC_vars.f90:41
grid_utilities::coord_e2f
Converts Equilibrium coordinates . to Flux coordinates .
Definition: grid_utilities.f90:66
num_vars::n_theta_plot
integer, public n_theta_plot
nr. of poloidal points in plot
Definition: num_vars.f90:156
grid_utilities::calc_xyz_grid
integer function, public calc_xyz_grid(grid_eq, grid_XYZ, X, Y, Z, L, R)
Calculates , and on a grid grid_XYZ, determined through its E(quilibrium) coordinates.
Definition: grid_utilities.f90:799
grid_utilities::trim_grid
integer function, public trim_grid(grid_in, grid_out, norm_id)
Trim a grid, removing any overlap between the different regions.
Definition: grid_utilities.f90:1636
eq_vars::b_0
real(dp), public b_0
derived normalization constant for nondimensionalization
Definition: eq_vars.f90:45
this documentation was created with Doxygen 1.8.17 by Toon Weyens