PB3D  [2.45]
Ideal linear high-n MHD stability in 3-D
HELENA_ops.f90
Go to the documentation of this file.
1 !------------------------------------------------------------------------------!
3 !------------------------------------------------------------------------------!
4 module helena_ops
5 #include <PB3D_macros.h>
6  use num_vars, only: pi
7  use str_utilities
8  use output_ops
9  use messages
10  use num_vars, only: dp, max_str_ln
11  use grid_vars, only: grid_type
12  use eq_vars, only: eq_1_type, eq_2_type
13  use x_vars, only: x_1_type, x_2_type
14  use helena_vars
15 
16  implicit none
17  private
18  public read_hel, interp_hel_on_grid
19 #if ldebug
20  public test_metrics_h, test_harm_cont_h
21 #endif
22 
23 contains
83  integer function read_hel(n_r_in,use_pol_flux_H) result(ierr)
84  use num_vars, only: eq_name, eq_i, max_deriv, tol_zero, &
85  &invert_top_bottom_h
86  use num_utilities, only: calc_int, spline
87  use helena_vars, only: pres_h, q_saf_h, rot_t_h, flux_p_h, flux_t_h, &
88  &nchi, chi_h, ias, rbphi_h, r_h, z_h, rmtog_h, bmtog_h
89 
90  character(*), parameter :: rout_name = 'read_HEL'
91 
92  ! input / output
93  integer, intent(inout) :: n_r_in
94  logical, intent(inout) :: use_pol_flux_h
95 
96  ! local variables
97  character(len=max_str_ln) :: err_msg ! error message
98  integer :: id, kd ! counters
99  integer :: nchi_loc ! local nchi
100  integer :: max_loc_z(2) ! location of maximum Z
101  real(dp), allocatable :: s_h(:) ! flux coordinate s
102  real(dp), allocatable :: curj(:) ! toroidal current
103  real(dp), allocatable :: vx(:), vy(:) ! R and Z of surface
104  real(dp) :: diff_dp ! Dp difference
105  real(dp) :: tol_diff_dp ! tolerance on diff_Dp
106  real(dp) :: dp0, dpe ! derivative of pressure on axis and surface
107  real(dp) :: dq0, dqe ! derivative of safety factor on axis and surface
108  real(dp) :: dj0, dje ! derivative of toroidal current on axis and surface
109  real(dp) :: df0, dfe ! derivative of RBphi on axis and surface
110  real(dp) :: radius, raxis ! minor radius, major radius normalized w.r.t. magnetic axis (i.e. R_m)
111  real(dp) :: b0 ! B_geo/B_m
112  real(dp) :: eps ! inverse aspect ratio
113  real(dp) :: cpsurf ! poloidal flux on surface
114  real(dp) :: ellip ! ellipticity
115  real(dp) :: tria ! triangularity
116  real(dp) :: delta_rz(2) ! R and Z range
117 #if ldebug
118  real(dp), allocatable :: inv_loc(:) ! used to invert top-bottom
119 #endif
120 
121  ! initialize ierr
122  ierr = 0
123 
124  ! user output
125  call writo('Reading data from HELENA output "' &
126  &// trim(eq_name) // '"')
127  call lvl_ud(1)
128 
129  ! set error messages
130  err_msg = 'Failed to read HELENA output file. &
131  &Does it output IAS and B0? See tutorial!'
132 
133  ! rewind input file
134  rewind(eq_i)
135 
136  ! read mapped equilibrium from disk
137  read(eq_i,*,iostat=ierr) n_r_in,ias ! nr. normal points (JS0), top-bottom symmetry
138  chckerr(err_msg)
139  n_r_in = n_r_in + 1 ! HELENA outputs nr. of normal points - 1
140 
141  allocate(s_h(n_r_in)) ! flux coordinate
142  read(eq_i,*,iostat=ierr) (s_h(kd),kd=1,n_r_in) ! it is squared below, after reading cpsurf
143  chckerr(err_msg)
144 
145  allocate(q_saf_h(n_r_in,0:max_deriv+1)) ! safety factor
146  read(eq_i,*,iostat=ierr) (q_saf_h(kd,0),kd=1,n_r_in)
147  chckerr(err_msg)
148  read(eq_i,*,iostat=ierr) dq0, dqe ! first point, last point
149  chckerr(err_msg)
150  ! Note: Dq0 an Dqe are in the local coordinate of the finite elements,
151  ! and are therefore not used here: q_saf_H(:,1) will be overwritten
152  q_saf_h(1,1) = dq0
153  read(eq_i,*,iostat=ierr) (q_saf_h(kd,1),kd=2,n_r_in) ! second to last point (again)
154  chckerr(err_msg)
155 
156  allocate(curj(n_r_in)) ! toroidal current
157  read(eq_i,*,iostat=ierr) (curj(kd),kd=1,n_r_in)
158  chckerr(err_msg)
159 
160  read(eq_i,*,iostat=ierr) dj0,dje ! derivative of toroidal current at first and last point
161  chckerr(err_msg)
162 
163  read(eq_i,*,iostat=ierr) nchi ! nr. poloidal points
164  chckerr(err_msg)
165  nchi_loc = nchi
166  if (ias.ne.0) nchi = nchi + 1 ! extend grid to 2pi if asymmetric (so doubling a point!)
167 
168  allocate(chi_h(nchi)) ! poloidal points
169  read(eq_i,*,iostat=ierr) (chi_h(id),id=1,nchi_loc)
170  chckerr(err_msg)
171  if (ias.ne.0) chi_h(nchi) = 2*pi
172 
173  allocate(h_h_11(nchi,n_r_in)) ! upper metric factor 1,1
174  read(eq_i,*,iostat=ierr) &
175  &(h_h_11(mod(id-1,nchi_loc)+1,(id-1)/nchi_loc+1),&
176  &id=nchi_loc+1,n_r_in*nchi_loc) ! (gem11)
177  chckerr(err_msg)
178  h_h_11(:,1) = tol_zero ! first normal point is not given, so set to almost zero
179  if (ias.ne.0) h_h_11(nchi,:) = h_h_11(1,:)
180 
181  allocate(h_h_12(nchi,n_r_in)) ! upper metric factor 1,2
182  read(eq_i,*,iostat=ierr) &
183  &(h_h_12(mod(id-1,nchi_loc)+1,(id-1)/nchi_loc+1),&
184  &id=nchi_loc+1,n_r_in*nchi_loc) ! (gem12)
185  chckerr(err_msg)
186  do id = 1,nchi-1
187  ierr = spline(s_h(2:n_r_in),h_h_12(id,2:n_r_in),s_h(1:1),&
188  &h_h_12(id,1:1),deriv=0,extrap=.true.)
189  chckerr('')
190  end do
191  if (ias.ne.0) h_h_12(nchi,:) = h_h_12(1,:)
192 
193  read(eq_i,*,iostat=ierr) cpsurf, radius ! poloidal flux on surface, minor radius
194  chckerr(err_msg)
195 
196  allocate(h_h_33(nchi,n_r_in)) ! lower metric factor 3,3
197  read(eq_i,*,iostat=ierr) &
198  &(h_h_33(mod(id-1,nchi_loc)+1,(id-1)/nchi_loc+1),id=&
199  &nchi_loc+1,n_r_in*nchi_loc) ! (gem33)
200 
201  read(eq_i,*,iostat=ierr) raxis, b0 ! major radius
202  chckerr(err_msg)
203  h_h_33(:,1) = raxis**2 ! first normal point is degenerate, major radius
204  if (ias.ne.0) h_h_33(nchi,:) = h_h_33(1,:)
205  h_h_33 = 1._dp/h_h_33 ! inverse is given
206 
207  allocate(pres_h(n_r_in,0:max_deriv+1)) ! pressure profile
208  read(eq_i,*,iostat=ierr) (pres_h(kd,0),kd=1,n_r_in)
209  chckerr(err_msg)
210  read(eq_i,*,iostat=ierr) dp0, dpe ! first point, last point
211  chckerr(err_msg)
212 
213  allocate(rbphi_h(n_r_in,0:max_deriv+1)) ! R B_phi (= F)
214  read(eq_i,*,iostat=ierr) (rbphi_h(kd,0),kd=1,n_r_in)
215  chckerr(err_msg)
216  ! Note: DF0 and DFe contain rubbish and are not used.
217  read(eq_i,*,iostat=ierr) df0, dfe ! derivatives of R B_phi on axis and surface
218  chckerr(err_msg)
219 
220  allocate(vx(nchi))
221  read(eq_i,*,iostat=ierr) (vx(id),id=1,nchi_loc) ! R on surface
222  chckerr(err_msg)
223  if (ias.ne.0) vx(nchi) = vx(1)
224 
225  allocate(vy(nchi))
226  read(eq_i,*,iostat=ierr) (vy(id),id=1,nchi_loc) ! Z on surface
227  chckerr(err_msg)
228  if (ias.ne.0) vy(nchi) = vy(1)
229 
230  read(eq_i,*,iostat=ierr) eps ! inverse aspect ratio
231  chckerr(err_msg)
232 
233  allocate(r_h(nchi,n_r_in)) ! major radius R
234  read(eq_i,*,iostat=ierr) (r_h(mod(id-1,nchi_loc)+1,(id-1)/nchi_loc+1),&
235  &id=nchi_loc+1,n_r_in*nchi_loc) ! (xout)
236  chckerr(err_msg)
237 
238  allocate(z_h(nchi,n_r_in)) ! height Z
239  read(eq_i,*,iostat=ierr) (z_h(mod(id-1,nchi_loc)+1,(id-1)/nchi_loc+1),&
240  &id=nchi_loc+1,n_r_in*nchi_loc) ! (yout)
241  chckerr(err_msg)
242 
243  ! transform to MISHKA normalization
244  allocate(flux_p_h(n_r_in,0:max_deriv+1))
245  flux_p_h(:,0) = 2*pi*s_h**2 * cpsurf ! rescale flux coordinate (HELENA uses psi_pol/2pi as flux)
246  radius = radius * raxis ! global length normalization with R_m
247  r_h = radius*(1._dp/eps + r_h) ! local normalization with a
248  z_h = radius*z_h ! local normalization with a
249  r_h(:,1) = raxis ! first normal point is degenerate, major radius
250  z_h(:,1) = 0._dp ! first normal point is degenerate, 0
251  if (ias.ne.0) r_h(nchi,:) = r_h(1,:)
252  if (ias.ne.0) z_h(nchi,:) = z_h(1,:)
253 #if ldebug
254  if (invert_top_bottom_h .and. ias.ne.0) then
255  call print_ex_2d('boundary','RZ',z_h(:,n_r_in),&
256  &x=r_h(:,n_r_in),draw=.false.)
257  call draw_ex(['boundary'],'RZ',1,1,.false.)
258  call writo('Inverting top and bottom to debug HELENA',&
259  &warning=.true.)
260  allocate(inv_loc(nchi))
261  ! chi_H
262  inv_loc = chi_h
263  chi_h = 2*pi-inv_loc(nchi:1:-1)
264  ! Z_H, R_H, h_H_11, h_H12, h_H_33
265  do kd = 1,n_r_in
266  inv_loc = z_h(:,kd)
267  z_h(:,kd) = -inv_loc(nchi:1:-1)
268  inv_loc = r_h(:,kd)
269  r_h(:,kd) = inv_loc(nchi:1:-1)
270  inv_loc = h_h_11(:,kd)
271  h_h_11(:,kd) = inv_loc(nchi:1:-1)
272  inv_loc = h_h_12(:,kd)
273  h_h_12(:,kd) = -inv_loc(nchi:1:-1)
274  inv_loc = h_h_33(:,kd)
275  h_h_33(:,kd) = inv_loc(nchi:1:-1)
276  end do
277  deallocate(inv_loc)
278  call print_ex_2d('boundary','RZ_inv_TB',z_h(:,n_r_in),&
279  &x=r_h(:,n_r_in),draw=.false.)
280  call draw_ex(['boundary'],'RZ_inv_TB',1,1,.false.)
281  end if
282 #endif
283 
284  ! set conversion factors
285  rmtog_h = radius/eps ! R_geo / R_mag
286  bmtog_h = b0 ! B_geo / B_mag
287 
288  ! calculate ellipticity and triangularity
289  max_loc_z = maxloc(z_h)
290  delta_rz(1) = maxval(r_h)-minval(r_h)
291  if (ias.eq.0) then
292  delta_rz(2) = 2*maxval(z_h)
293  else
294  delta_rz(2) = maxval(z_h)-minval(z_h)
295  end if
296  ellip = delta_rz(2)/delta_rz(1)
297  tria = ((maxval(r_h)+minval(r_h))*0.5_dp-&
298  &r_h(max_loc_z(1),max_loc_z(2)))/&
299  &(delta_rz(2)*0.5_dp)
300  call writo('Calculated ellipticity: '//trim(r2str(ellip)))
301  call writo('Calculated triangularity: '//trim(r2str(tria)))
302 
303  ! calculate derivatives of fluxes
304  allocate(flux_t_h(n_r_in,0:max_deriv+1))
305  flux_p_h(:,1) = 2*pi
306  flux_p_h(:,2:max_deriv+1) = 0._dp
307  flux_t_h(:,1) = q_saf_h(:,0)*flux_p_h(:,1)
308  ierr = calc_int(flux_t_h(:,1),flux_p_h(:,0)/(2*pi),flux_t_h(:,0))
309  chckerr('')
310  do kd = 2,max_deriv+1
311  ierr = spline(flux_p_h(:,0)/(2*pi),flux_t_h(:,1),&
312  &flux_p_h(:,0)/(2*pi),flux_t_h(:,kd),deriv=kd-1)
313  chckerr('')
314  end do
315 
316  ! calculate derivatives of other variables
317  allocate(rot_t_h(n_r_in,0:max_deriv+1)) ! rotational transform
318  rot_t_h(:,0) = 1._dp/q_saf_h(:,0)
319  do kd = 1,max_deriv+1
320  ierr = spline(flux_p_h(:,0)/(2*pi),pres_h(:,0),&
321  &flux_p_h(:,0)/(2*pi),pres_h(:,kd),deriv=kd,bcs=[0,1],&
322  &bcs_val=[0._dp,dpe*pi/(flux_p_h(n_r_in,0)*s_h(n_r_in))])
323  chckerr('')
324  ierr = spline(flux_p_h(:,0)/(2*pi),rbphi_h(:,0),&
325  &flux_p_h(:,0)/(2*pi),rbphi_h(:,kd),deriv=kd)
326  chckerr('')
327  ierr = spline(flux_p_h(:,0)/(2*pi),rot_t_h(:,0),&
328  &flux_p_h(:,0)/(2*pi),rot_t_h(:,kd),deriv=kd)
329  chckerr('')
330  ierr = spline(flux_p_h(:,0)/(2*pi),q_saf_h(:,0),&
331  &flux_p_h(:,0)/(2*pi),q_saf_h(:,kd),deriv=kd)
332  chckerr('')
333  end do
334 
335  ! check for consistency
336  tol_diff_dp = abs(pres_h(1,0)-pres_h(n_r_in,0))/n_r_in
337  diff_dp = s_h(1)*flux_p_h(n_r_in,0)/pi*pres_h(1,1)-dp0
338  if (abs(diff_dp).gt.tol_diff_dp) then
339  call writo('Difference between pressure derivative on axis is &
340  &large ('//trim(r2str(diff_dp))//')')
341  call writo('Maybe increase precision or number of normal points',&
342  &warning=.true.)
343  end if
344  diff_dp = s_h(n_r_in)*flux_p_h(n_r_in,0)/pi*pres_h(n_r_in,1)-dpe
345  if (abs(diff_dp).gt.tol_diff_dp) then
346  call writo('Difference between pressure derivative at edge is &
347  &large ('//trim(r2str(diff_dp))//')')
348  call writo('Maybe increase precision or number of normal points',&
349  &warning=.true.)
350  end if
351 
352  !!! To plot the cross-section
353  !!call print_ex_2D(['cross_section'],'cross_section_H',Z_H,x=R_H,&
354  !!&draw=.false.)
355  !!call draw_ex(['cross_section'],'cross_section_H',n_r_in,1,.false.)
356  !!call print_ex_2D(['cross_section_inv'],'cross_section_inv_H',&
357  !!&transpose(Z_H),x=transpose(R_H),draw=.false.)
358  !!call draw_ex(['cross_section_inv'],'cross_section_inv_H',nchi,1,.false.)
359 
360  ! HELENA always uses the poloidal flux
361  use_pol_flux_h = .true.
362 
363  ! user output
364  call writo('HELENA output given on '//trim(i2str(nchi))//&
365  &' poloidal and '//trim(i2str(n_r_in))//' normal points')
366  if (ias.eq.0) call writo('The equilibrium is top-bottom symmetric')
367  call lvl_ud(-1)
368  call writo('Data from HELENA output succesfully read')
369  end function read_hel
370 
411  integer function get_ang_interp_data_hel(grid_in,grid_out,theta_i,&
412  &fund_theta_int_displ,tb_sym,use_E) result(ierr)
413  use num_utilities, only: con2dis
414 
415  character(*), parameter :: rout_name = 'get_ang_interp_data_HEL'
416 
417  ! input / output
418  type(grid_type), intent(in) :: grid_in
419  type(grid_type), intent(in) :: grid_out
420  real(dp), allocatable, intent(inout) :: theta_i(:,:,:)
421  integer, allocatable, intent(inout) :: fund_theta_int_displ(:,:,:)
422  logical, intent(in), optional :: tb_sym
423  logical, intent(in), optional :: use_e
424 
425  ! local variables
426  integer :: id, jd, kd ! counters
427  logical :: tb_sym_loc ! local version of tb_sym
428  logical :: use_e_loc ! local version of use_E
429  real(dp) :: theta_loc ! local theta of output grid
430  real(dp), pointer :: theta_in(:) => null() ! input theta
431  character(len=max_str_ln) :: err_msg ! error message
432  integer :: fund_theta_int_lo, fund_theta_int_hi ! lower and upper limit of fundamental theta interval
433 
434  ! initialize ierr
435  ierr = 0
436 
437  ! test whether axisymmetric grid
438  if (grid_in%n(2).ne.1) then
439  ierr = 1
440  err_msg = 'Not an axisymmetric grid'
441  chckerr(err_msg)
442  end if
443  ! test whether normal sizes compatible
444  if (grid_in%loc_n_r.ne.grid_out%loc_n_r) then
445  ierr = 1
446  write(*,*) grid_in%loc_n_r, grid_out%loc_n_r
447  err_msg = 'Grids are not compatible in normal direction'
448  chckerr(err_msg)
449  end if
450 
451  ! set up local use_E and tb_sym
452  use_e_loc = .false.
453  if (present(use_e)) use_e_loc = use_e
454  tb_sym_loc = .false.
455  if (present(tb_sym)) tb_sym_loc = tb_sym
456 
457  ! allocate theta_i and fundamental theta interval displacement
458  allocate(theta_i(grid_out%n(1),grid_out%n(2),grid_out%loc_n_r))
459  allocate(fund_theta_int_displ(grid_out%n(1),grid_out%n(2),&
460  &grid_out%loc_n_r))
461 
462  ! initialize fundamental theta interval displacement
463  fund_theta_int_displ = 0
464 
465  ! For every point on output grid, check to which half poloidal circle on
466  ! the axisymmetric input grid it belongs to. If there is top-down
467  ! symmetry and the angle theta lies in the bottom part, the quantities
468  ! have to be taken from their symmetric counterpart (2pi-theta).
469  ! loop over all normal points of this rank
470  do kd = 1,grid_out%loc_n_r
471  ! point theta_in
472  if (use_e) then
473  theta_in => grid_in%theta_E(:,1,kd) ! axisymmetric grid should have only one toroidal point
474  else
475  theta_in => grid_in%theta_F(:,1,kd) ! axisymmetric grid should have only one toroidal point
476  end if
477 
478  ! loop over all angles of ang_2
479  do jd = 1,grid_out%n(2)
480  ! loop over all angles of ang_1
481  do id = 1,grid_out%n(1)
482  ! set the local poloidal point from theta
483  if (use_e) then
484  theta_loc = grid_out%theta_E(id,jd,kd)
485  else
486  theta_loc = grid_out%theta_F(id,jd,kd)
487  end if
488 
489  ! set min
490  if (tb_sym) then
491  fund_theta_int_lo = -1
492  fund_theta_int_hi = 1
493  else
494  fund_theta_int_lo = 0
495  fund_theta_int_hi = 2
496  end if
497 
498  ! add or subtract 2pi to the poloidal angle until it is in
499  ! the fundamental range indicated by fund_theta_int_lo and
500  ! fund_theta_int_hi
501  if (theta_loc.lt.fund_theta_int_lo*pi) then
502  do while (theta_loc.lt.fund_theta_int_lo*pi)
503  theta_loc = theta_loc + 2*pi
504  fund_theta_int_displ(id,jd,kd) = &
505  &fund_theta_int_displ(id,jd,kd) - 1 ! theta interval lies to the left of fund. theta interval
506  end do
507  else if (theta_loc.gt.fund_theta_int_hi*pi) then
508  do while (theta_loc.gt.fund_theta_int_hi*pi)
509  theta_loc = theta_loc - 2*pi
510  fund_theta_int_displ(id,jd,kd) = &
511  &fund_theta_int_displ(id,jd,kd) + 1 ! theta interval lies to the right of fund. theta interval
512  end do
513  end if
514  ierr = con2dis(abs(theta_loc),theta_i(id,jd,kd),theta_in)
515  chckerr('')
516  if (theta_loc.lt.0) theta_i(id,jd,kd) = - theta_i(id,jd,kd)
517  end do
518  end do
519  end do
520 
521  ! clean up
522  nullify(theta_in)
523  end function get_ang_interp_data_hel
524 
548  integer function interp_hel_on_grid(grid_in,grid_out,eq_2,X_1,X_2,&
549  &eq_2_out,X_1_out,X_2_out,eq_1,grid_name) result(ierr)
550 
551  use num_vars, only: use_pol_flux_f
552 
553  character(*), parameter :: rout_name = 'interp_HEL_on_grid'
554 
555  ! input / output
556  type(grid_type), intent(in) :: grid_in
557  type(grid_type), intent(in) :: grid_out
558  type(eq_2_type), intent(inout), optional :: eq_2
559  type(x_1_type), intent(inout), optional :: x_1
560  type(x_2_type), intent(inout), optional :: x_2
561  type(eq_2_type), intent(inout), optional :: eq_2_out
562  type(x_1_type), intent(inout), optional :: x_1_out
563  type(x_2_type), intent(inout), optional :: x_2_out
564  type(eq_1_type), intent(in), optional :: eq_1
565  character(len=*), intent(in), optional :: grid_name
566 
567  ! local variables
568  real(dp), allocatable :: theta_i(:,:,:) ! interpolation index
569  integer, allocatable :: fund_theta_int_displ(:,:,:) ! displacement of fundamental theta interval
570  logical :: tb_sym ! whether there is top-bottom symmetry (relevant only for HELENA)
571  character(len=max_str_ln) :: err_msg ! error message
572 
573  ! initialize ierr
574  ierr = 0
575 
576  ! user output
577  if (present(grid_name)) then
578  call writo('Adapt quantities to '//trim(grid_name))
579  call lvl_ud(1)
580  end if
581 
582  ! tests
583  if (present(eq_2).and..not.present(eq_1)) then
584  ierr =1
585  err_msg = 'To interpolate metric equilibrium variables, flux &
586  &equilibrium has to be provided as well through eq_1'
587  chckerr(err_msg)
588  end if
589  if (present(eq_2) .and. (present(x_1).or.present(x_2))) then
590  ierr = 1
591  err_msg = 'Only the quantities on one type of grid can be &
592  &interpolated in a single call'
593  chckerr(err_msg)
594  end if
595 
596  ! set up tb_sym
597  if (ias.eq.0) then
598  tb_sym = .true.
599  else
600  tb_sym = .false.
601  end if
602 
603  ! get angular interpolation factors
604  ierr = get_ang_interp_data_hel(grid_in,grid_out,theta_i,&
605  &fund_theta_int_displ,tb_sym=tb_sym,use_e=.false.)
606  chckerr('')
607 
608  ! adapt equilibrium quantities
609  if (present(eq_2)) then
610  if (present(eq_2_out)) then
611  ierr = interp_var_6d_real(eq_2%jac_FD,eq_2_out%jac_FD)
612  chckerr('')
613  ierr = interp_var_7d_real(eq_2%g_FD,eq_2_out%g_FD,sym_type=3)
614  chckerr('')
615  ierr = interp_var_7d_real(eq_2%h_FD,eq_2_out%h_FD,sym_type=4)
616  chckerr('')
617  ierr = interp_var_3d_real(eq_2%S,eq_2_out%S)
618  chckerr('')
619  ierr = interp_var_3d_real(eq_2%sigma,eq_2_out%sigma)
620  chckerr('')
621  ierr = interp_var_3d_real(eq_2%kappa_n,eq_2_out%kappa_n)
622  chckerr('')
623  ierr = interp_var_3d_real(eq_2%kappa_g,eq_2_out%kappa_g,&
624  &sym_type=2)
625  chckerr('')
626  else
627  ierr = interp_var_6d_real_ow(eq_2%jac_FD)
628  chckerr('')
629  ierr = interp_var_7d_real_ow(eq_2%g_FD,sym_type=3)
630  chckerr('')
631  ierr = interp_var_7d_real_ow(eq_2%h_FD,sym_type=4)
632  chckerr('')
633  ierr = interp_var_3d_real_ow(eq_2%S)
634  chckerr('')
635  ierr = interp_var_3d_real_ow(eq_2%sigma)
636  chckerr('')
637  ierr = interp_var_3d_real_ow(eq_2%kappa_n)
638  chckerr('')
639  ierr = interp_var_3d_real_ow(eq_2%kappa_g,sym_type=2)
640  chckerr('')
641  end if
642  end if
643 
644  ! adapt vectorial perturbation quantities
645  if (present(x_1)) then
646  if (present(x_1_out)) then
647  ierr = interp_var_4d_complex(x_1%U_0,x_1_out%U_0,sym_type=2)
648  chckerr('')
649  ierr = interp_var_4d_complex(x_1%U_1,x_1_out%U_1,sym_type=2)
650  chckerr('')
651  ierr = interp_var_4d_complex(x_1%DU_0,x_1_out%DU_0,sym_type=1)
652  chckerr('')
653  ierr = interp_var_4d_complex(x_1%DU_1,x_1_out%DU_1,sym_type=1)
654  chckerr('')
655  else
656  ierr = interp_var_4d_complex_ow(x_1%U_0,sym_type=2)
657  chckerr('')
658  ierr = interp_var_4d_complex_ow(x_1%U_1,sym_type=2)
659  chckerr('')
660  ierr = interp_var_4d_complex_ow(x_1%DU_0,sym_type=1)
661  chckerr('')
662  ierr = interp_var_4d_complex_ow(x_1%DU_1,sym_type=1)
663  chckerr('')
664  end if
665  end if
666 
667  ! adapt tensorial perturbation quantities
668  if (present(x_2)) then
669  if (present(x_2_out)) then
670  ierr = interp_var_4d_complex(x_2%PV_0,x_2_out%PV_0,sym_type=1)
671  chckerr('')
672  ierr = interp_var_4d_complex(x_2%PV_1,x_2_out%PV_1,sym_type=1)
673  chckerr('')
674  ierr = interp_var_4d_complex(x_2%PV_2,x_2_out%PV_2,sym_type=1)
675  chckerr('')
676  ierr = interp_var_4d_complex(x_2%KV_0,x_2_out%KV_0,sym_type=1)
677  chckerr('')
678  ierr = interp_var_4d_complex(x_2%KV_1,x_2_out%KV_1,sym_type=1)
679  chckerr('')
680  ierr = interp_var_4d_complex(x_2%KV_2,x_2_out%KV_2,sym_type=1)
681  chckerr('')
682  else
683  ierr = interp_var_4d_complex_ow(x_2%PV_0,sym_type=1)
684  chckerr('')
685  ierr = interp_var_4d_complex_ow(x_2%PV_1,sym_type=1)
686  chckerr('')
687  ierr = interp_var_4d_complex_ow(x_2%PV_2,sym_type=1)
688  chckerr('')
689  ierr = interp_var_4d_complex_ow(x_2%KV_0,sym_type=1)
690  chckerr('')
691  ierr = interp_var_4d_complex_ow(x_2%KV_1,sym_type=1)
692  chckerr('')
693  ierr = interp_var_4d_complex_ow(x_2%KV_2,sym_type=1)
694  chckerr('')
695  end if
696  end if
697 
698  ! user output
699  if (present(grid_name)) then
700  call lvl_ud(-1)
701  end if
702  contains
703  ! Interpolate a variable defined on an axisymmetric grid at poloidal
704  ! angles indicated by the interpolation factors theta_i.
705  ! There is an optional variable sym_type that allows for extra
706  ! operations to be done on the variable if top-down symmetry is applied:
707  ! - sym_type = 0: var(-theta) = var(theta)
708  ! - sym_type = 1: var(-theta) = var(theta)*
709  ! - sym_type = 2: var(-theta) = -var(theta)*
710  ! - sym_type = 3: var(-theta) = +/-var(theta)* for g_FD (See note)
711  ! - sym_type = 4: var(-theta) = +/-var(theta)* for h_FD (See note)
712  ! where symmetry type 5 is only valid for tensorial quantities (7D). The
713  ! transformation matrix, a function of only the flux coordinate, has to
714  ! be passed as well. This is useful for the transformations of the
715  ! metric quantities as seen in [ADD REF].
716  ! When top-down symmetry has been used to calculate the interpolation
717  ! factors, this is indicated by a negative factor instead of a positive
718  ! one.
719  ! Note: This is also correct if i_hi = i_lo.
720  ! Note: For the overwrite versions ('ow'), the lower limits of
721  ! indices are 0 for the dimensions corresponding to derivatives. These
722  ! allocations are done here and persist into the calling functions.
723  ! Note: The 6D and 7D overwrite versions take a pointer instead of a
724  ! allocatable variable because they are to be used for jacobians and
725  ! metric coefficients, which are defined as pointers.
726  ! Note: symmetry type 3 for 7D variables (i.e. metric coefficients), the
727  ! relation var(-theta) = (-1)^(i+j) var(theta)* where i and j are the
728  ! indices in g_ij and h^ij. Furthermore, since the metric coefficients
729  ! are for the Flux coordinate system, they need to be adapted for the
730  ! secular coordinate theta or phi.
732  integer function interp_var_3d_real(var,var_int,sym_type) result(ierr) ! 3D_real version, separate output variable
733  ! input / output
734  real(dp), intent(in) :: var(1:,1:,1:) ! variable to be interpolated
735  real(dp), intent(inout) :: var_int(1:,1:,1:) ! interpolated var
736  integer, intent(in), optional :: sym_type ! optionally another type of symmetry
737 
738  ! local variables
739  integer :: i_lo, i_hi ! upper and lower index
740  integer :: id, jd, kd ! counters
741  integer :: sym_type_loc ! local version of symmetry type
742 
743  ! initialize ierr
744  ierr = 0
745 
746  ! set up local symmetry type
747  sym_type_loc = 0
748  if (present(sym_type)) sym_type_loc = sym_type
749 
750  ! test symmetry types
751  if (sym_type_loc.lt.0 .or. sym_type_loc.gt.2) then
752  ierr = 1
753  err_msg = 'Nonsensible symmetry type '//&
754  &trim(i2str(sym_type_loc))
755  chckerr(err_msg)
756  end if
757 
758  ! iterate over all normal points
759  do kd = 1,size(var_int,3)
760  ! iterate over all geodesical points
761  do jd = 1,size(var_int,2)
762  ! iterate over all parallel points
763  do id = 1,size(var_int,1)
764  ! set up i_lo and i_hi
765  i_lo = floor(abs(theta_i(id,jd,kd)))
766  i_hi = ceiling(abs(theta_i(id,jd,kd)))
767 
768  var_int(id,jd,kd) = var(i_lo,1,kd) + &
769  &(var(i_hi,1,kd)-var(i_lo,1,kd))*&
770  &(abs(theta_i(id,jd,kd))-i_lo) ! because i_hi - i_lo = 1
771  if (theta_i(id,jd,kd).lt.0) then
772  if (sym_type_loc.eq.2) var_int(id,jd,kd) = &
773  &- var_int(id,jd,kd)
774  end if
775  end do
776  end do
777  end do
778  end function interp_var_3d_real
780  integer function interp_var_3d_real_ow(var,sym_type) result(ierr) ! 3D_real version, overwriting variable
781  ! input / output
782  real(dp), intent(inout), allocatable :: var(:,:,:) ! variable to be interpolated
783  integer, intent(in), optional :: sym_type ! optionally another type of symmetry
784 
785  ! local variables
786  real(dp), allocatable :: var_3d_loc(:,:,:) ! local 3D variable
787  integer :: var_dim(3,2) ! dimensions of variable
788 
789  ! initialize ierr
790  ierr = 0
791 
792  ! set up variable dimensions
793  var_dim(:,1) = [1,1,1]
794  var_dim(:,2) = [shape(theta_i)]
795 
796  ! set up local 3D variable
797  allocate(var_3d_loc(var_dim(1,1):var_dim(1,2),&
798  &var_dim(2,1):var_dim(2,2),var_dim(3,1):var_dim(3,2)))
799 
800  ! call normal routine
801  ierr = interp_var_3d_real(var,var_3d_loc,sym_type)
802  chckerr('')
803 
804  ! reallocate variable
805  deallocate(var)
806  allocate(var(var_dim(1,1):var_dim(1,2),&
807  &var_dim(2,1):var_dim(2,2),var_dim(3,1):var_dim(3,2)))
808 
809  ! overwrite variable
810  var = var_3d_loc
811  end function interp_var_3d_real_ow
813  integer function interp_var_4d_complex(var,var_int,sym_type) &
814  &result(ierr) ! 4D_complex version, separate output variable
815 
816  ! input / output
817  complex(dp), intent(in) :: var(1:,1:,1:,1:) ! variable to be interpolated
818  complex(dp), intent(inout) :: var_int(1:,1:,1:,1:) ! interpolated var
819  integer, intent(in), optional :: sym_type ! optionally another type of symmetry
820 
821  ! local variables
822  integer :: i_lo, i_hi ! upper and lower index
823  integer :: id, jd, kd ! counters
824  integer :: sym_type_loc ! local version of symmetry type
825 
826  ! initialize ierr
827  ierr = 0
828 
829  ! set up local symmetry type
830  sym_type_loc = 0
831  if (present(sym_type)) sym_type_loc = sym_type
832 
833  ! test symmetry types
834  if (sym_type_loc.lt.0 .or. sym_type_loc.gt.2) then
835  ierr = 1
836  err_msg = 'Nonsensible symmetry type '//&
837  &trim(i2str(sym_type_loc))
838  chckerr(err_msg)
839  end if
840 
841  ! iterate over all normal points
842  do kd = 1,size(var_int,3)
843  ! iterate over all geodesical points
844  do jd = 1,size(var_int,2)
845  ! iterate over all parallel points
846  do id = 1,size(var_int,1)
847  ! set up i_lo and i_hi
848  i_lo = floor(abs(theta_i(id,jd,kd)))
849  i_hi = ceiling(abs(theta_i(id,jd,kd)))
850 
851  var_int(id,jd,kd,:) = var(i_lo,1,kd,:) + &
852  &(var(i_hi,1,kd,:)-var(i_lo,1,kd,:))*&
853  &(abs(theta_i(id,jd,kd))-i_lo) ! because i_hi - i_lo = 1
854  if (theta_i(id,jd,kd).lt.0) then
855  if (sym_type_loc.eq.1) var_int(id,jd,kd,:) = &
856  &conjg(var_int(id,jd,kd,:))
857  if (sym_type_loc.eq.2) var_int(id,jd,kd,:) = &
858  &- conjg(var_int(id,jd,kd,:))
859  end if
860  end do
861  end do
862  end do
863  end function interp_var_4d_complex
865  integer function interp_var_4d_complex_ow(var,sym_type) result(ierr) ! 4D_complex version, overwriting variable
866  ! input / output
867  complex(dp), intent(inout), allocatable :: var(:,:,:,:) ! variable to be interpolated
868  integer, intent(in), optional :: sym_type ! optionally another type of symmetry
869 
870  ! local variables
871  complex(dp), allocatable :: var_4d_loc(:,:,:,:) ! local 4D variable
872  integer :: var_dim(4,2) ! dimensions of variable
873 
874  ! initialize ierr
875  ierr = 0
876 
877  ! set up variable dimensions
878  var_dim(:,1) = [1,1,1,1]
879  var_dim(:,2) = [shape(theta_i),size(var,4)]
880 
881  ! set up local 4D variable
882  allocate(var_4d_loc(var_dim(1,1):var_dim(1,2),&
883  &var_dim(2,1):var_dim(2,2),var_dim(3,1):var_dim(3,2),&
884  &var_dim(4,1):var_dim(4,2)))
885 
886  ! call normal routine
887  ierr = interp_var_4d_complex(var,var_4d_loc,sym_type)
888  chckerr('')
889 
890  ! reallocate variable
891  deallocate(var)
892  allocate(var(var_dim(1,1):var_dim(1,2),&
893  &var_dim(2,1):var_dim(2,2),var_dim(3,1):var_dim(3,2),&
894  &var_dim(4,1):var_dim(4,2)))
895 
896  ! overwrite variable
897  var = var_4d_loc
898  end function interp_var_4d_complex_ow
900  integer function interp_var_6d_real(var,var_int,sym_type) result(ierr) ! 6D_real version, separate output variable
901  ! input / output
902  real(dp), intent(in) :: var(1:,1:,1:,0:,0:,0:) ! variable to be interpolated
903  real(dp), intent(inout) :: var_int(1:,1:,1:,0:,0:,0:) ! interpolated var
904  integer, intent(in), optional :: sym_type ! optionally another type of symmetry
905 
906  ! local variables
907  integer :: i_lo, i_hi ! upper and lower index
908  integer :: id, jd, kd, ld ! counters
909  integer :: sym_type_loc ! local version of symmetry type
910 
911  ! initialize ierr
912  ierr = 0
913 
914  ! set up local symmetry type
915  sym_type_loc = 0
916  if (present(sym_type)) sym_type_loc = sym_type
917 
918  ! test symmetry types
919  if (sym_type_loc.lt.0 .or. sym_type_loc.gt.2) then
920  ierr = 1
921  err_msg = 'Nonsensible symmetry type '//&
922  &trim(i2str(sym_type_loc))
923  chckerr(err_msg)
924  end if
925 
926  ! iterate over all normal points
927  do kd = 1,size(var_int,3)
928  ! iterate over all geodesical points
929  do jd = 1,size(var_int,2)
930  ! iterate over all parallel points
931  do id = 1,size(var_int,1)
932  ! set up i_lo and i_hi
933  i_lo = floor(abs(theta_i(id,jd,kd)))
934  i_hi = ceiling(abs(theta_i(id,jd,kd)))
935 
936  var_int(id,jd,kd,:,:,:) = var(i_lo,1,kd,:,:,:) + &
937  &(var(i_hi,1,kd,:,:,:)-var(i_lo,1,kd,:,:,:))*&
938  &(abs(theta_i(id,jd,kd))-i_lo) ! because i_hi - i_lo = 1
939  if (theta_i(id,jd,kd).lt.0) then
940  if (sym_type_loc.eq.2) var_int(id,jd,kd,:,:,:) = &
941  &- var_int(id,jd,kd,:,:,:)
942  ! adapt the derivatives
943  if (use_pol_flux_f) then
944  do ld = 0,size(var_int,6)-1
945  var_int(id,jd,kd,:,:,ld) = &
946  &(-1)**ld*var_int(id,jd,kd,:,:,ld)
947  end do
948  else
949  ierr = 1
950  err_msg = '!!! INTERP_VAR_6D_REAL NOT YET &
951  &IMPLEMENTED FOR TOROIDAL FLUX !!!'
952  chckerr(err_msg)
953  end if
954  end if
955  end do
956  end do
957  end do
958  end function interp_var_6d_real
960  integer function interp_var_6d_real_ow(var,sym_type) result(ierr) ! 6D_real version, overwriting variable
961  ! input / output
962  real(dp), intent(inout), allocatable :: var(:,:,:,:,:,:) ! variable to be interpolated
963  integer, intent(in), optional :: sym_type ! optionally another type of symmetry
964 
965  ! local variables
966  real(dp), allocatable :: var_6d_loc(:,:,:,:,:,:) ! local 6D variable
967  integer :: var_dim(6,2) ! dimensions of variable
968 
969  ! initialize ierr
970  ierr = 0
971 
972  ! set up variable dimensions
973  var_dim(:,1) = [1,1,1,0,0,0]
974  var_dim(:,2) = [shape(theta_i),size(var,4)-1,size(var,5)-1,&
975  &size(var,6)-1]
976 
977  ! set up local 6D variable
978  allocate(var_6d_loc(var_dim(1,1):var_dim(1,2),&
979  &var_dim(2,1):var_dim(2,2),var_dim(3,1):var_dim(3,2),&
980  &var_dim(4,1):var_dim(4,2),var_dim(5,1):var_dim(5,2),&
981  &var_dim(6,1):var_dim(6,2)))
982 
983  ! call normal routine
984  ierr = interp_var_6d_real(var,var_6d_loc,sym_type)
985  chckerr('')
986 
987  ! reallocate variable
988  deallocate(var)
989  allocate(var(var_dim(1,1):var_dim(1,2),&
990  &var_dim(2,1):var_dim(2,2),var_dim(3,1):var_dim(3,2),&
991  &var_dim(4,1):var_dim(4,2),var_dim(5,1):var_dim(5,2),&
992  &var_dim(6,1):var_dim(6,2)))
993 
994  ! overwrite variable
995  var = var_6d_loc
996  end function interp_var_6d_real_ow
998  integer function interp_var_7d_real(var,var_int,sym_type) result(ierr) ! 7D_real version, separate output variable
999  use num_utilities, only: calc_mult, c, add_arr_mult, derivs
1000  use num_vars, only: max_deriv
1001 
1002  ! input / output
1003  real(dp), intent(in) :: var(1:,1:,1:,1:,0:,0:,0:) ! variable to be interpolated
1004  real(dp), intent(inout) :: var_int(1:,1:,1:,1:,0:,0:,0:) ! interpolated var
1005  integer, intent(in), optional :: sym_type ! optionally another type of symmetry
1006 
1007  ! local variables
1008  integer :: i_lo, i_hi ! upper and lower index
1009  integer :: id, jd, kd, ld ! counters
1010  integer :: sym_type_loc ! local version of symmetry type
1011  real(dp), allocatable :: t_met(:,:,:,:,:,:) ! transformation matrix for metric elements
1012  integer, allocatable :: derivs_loc(:,:) ! array of all unique derivatives
1013  integer :: c_loc(2) ! local c
1014 
1015  ! initialize ierr
1016  ierr = 0
1017 
1018  ! set up local symmetry type
1019  sym_type_loc = 0
1020  if (present(sym_type)) sym_type_loc = sym_type
1021 
1022  ! test symmetry types
1023  if (sym_type_loc.lt.0 .or. sym_type_loc.gt.4) then
1024  ierr = 1
1025  err_msg = 'Nonsensible symmetry type '//&
1026  &trim(i2str(sym_type_loc))
1027  chckerr(err_msg)
1028  end if
1029 
1030  ! interpolate
1031  do kd = 1,size(var_int,3)
1032  do jd = 1,size(var_int,2)
1033  do id = 1,size(var_int,1)
1034  i_lo = floor(abs(theta_i(id,jd,kd)))
1035  i_hi = ceiling(abs(theta_i(id,jd,kd)))
1036 
1037  var_int(id,jd,kd,:,:,:,:) = var(i_lo,1,kd,:,:,:,:) + &
1038  &(var(i_hi,1,kd,:,:,:,:)-var(i_lo,1,kd,:,:,:,:))*&
1039  &(abs(theta_i(id,jd,kd))-i_lo) ! because i_hi - i_lo = 1
1040  end do
1041  end do
1042  end do
1043 
1044  ! inversion of poloidal coordinate if TBS and symmetry type 2, 3, 4
1045  if (sym_type_loc.eq.2 .or. sym_type_loc.eq.3 .or. &
1046  &sym_type_loc.eq.4) then
1047  do kd = 1,size(var_int,3)
1048  do jd = 1,size(var_int,2)
1049  do id = 1,size(var_int,1)
1050  if (theta_i(id,jd,kd).lt.0) then
1051  ! invert value
1052  if (sym_type_loc.eq.2) then
1053  var_int(id,jd,kd,:,:,:,:) = &
1054  &- var_int(id,jd,kd,:,:,:,:)
1055  else if (sym_type_loc.eq.3 .or. &
1056  &sym_type_loc.eq.4) then ! metric coefficients: assuming 3D and symmetric
1057  c_loc(1) = c([2,1],.true.)
1058  c_loc(2) = c([3,2],.true.)
1059  var_int(id,jd,kd,c_loc(1),:,:,:) = &
1060  &- var_int(id,jd,kd,c_loc(1),:,:,:)
1061  var_int(id,jd,kd,c_loc(2),:,:,:) = &
1062  &- var_int(id,jd,kd,c_loc(2),:,:,:)
1063  end if
1064  ! invert derivatives
1065  if (use_pol_flux_f) then ! invert parallel derivative
1066  do ld = 0,size(var_int,7)-1
1067  var_int(id,jd,kd,:,:,:,ld) = (-1)**ld*&
1068  &var_int(id,jd,kd,:,:,:,ld)
1069  end do
1070  else
1071  ierr = 1
1072  err_msg = '!!! INTERP_VAR_7D_REAL NOT YET &
1073  &IMPLEMENTED FOR TOROIDAL FLUX !!!'
1074  chckerr(err_msg)
1075  end if
1076  end if
1077  end do
1078  end do
1079  end do
1080  end if
1081 
1082  ! set up displacement matrix if symmetry type 3 or 4
1083  if (sym_type_loc.eq.3 .or. sym_type_loc.eq.4) then
1084  allocate(t_met(size(var_int,1),size(var_int,2),&
1085  &size(var_int,3),0:size(var_int,5)-1,&
1086  &0:size(var_int,6)-1,0:size(var_int,7)-1))
1087  t_met = 0._dp
1088  if (use_pol_flux_f) then
1089  do ld = 0,max_deriv
1090  do kd = 1,size(theta_i,3)
1091  t_met(:,:,kd,0,ld,0) = &
1092  &2*pi*eq_1%q_saf_FD(kd,ld+1)
1093  end do
1094  end do
1095  else
1096  do ld = 0,max_deriv
1097  do kd = 1,size(theta_i,3)
1098  t_met(:,:,kd,0,ld,0) = &
1099  &-2*pi*eq_1%rot_t_FD(kd,ld+1)
1100  end do
1101  end do
1102  end if
1103  do kd = 1,size(theta_i,3)
1104  do jd = 1,size(theta_i,2)
1105  do id = 1,size(theta_i,1)
1106  t_met(id,jd,kd,:,:,:) = t_met(id,jd,kd,:,:,:)*&
1107  &fund_theta_int_displ(id,jd,kd)
1108  end do
1109  end do
1110  end do
1111  end if
1112 
1113  ! correct g_2i
1114  if (sym_type_loc.eq.3) then ! lower metric elements h_FD
1115  if (use_pol_flux_f) then ! poloidal flux coordinates
1116  do id = 0,max_deriv
1117  derivs_loc = derivs(id)
1118  do jd = 1,size(derivs_loc,2)
1119  ! g_22 -> g_22 + T_met * g_32
1120  c_loc(1) = c([3,2],.true.)
1121  c_loc(2) = c([2,2],.true.)
1122  ierr = add_arr_mult(&
1123  &var_int(:,:,:,c_loc(1),:,:,:),t_met,&
1124  &var_int(:,:,:,c_loc(2),&
1125  &derivs_loc(1,jd),derivs_loc(2,jd),&
1126  &derivs_loc(3,jd)),derivs_loc(:,jd))
1127  chckerr('')
1128  ! g_32 -> g_32 + T_met * g_31
1129  c_loc(1) = c([3,1],.true.)
1130  c_loc(2) = c([3,2],.true.)
1131  ierr = add_arr_mult(&
1132  &var_int(:,:,:,c_loc(1),:,:,:),t_met,&
1133  &var_int(:,:,:,c_loc(2),&
1134  &derivs_loc(1,jd),derivs_loc(2,jd),&
1135  &derivs_loc(3,jd)),derivs_loc(:,jd))
1136  chckerr('')
1137  ! g_22 -> g_22 + T_met * g_32
1138  c_loc(1) = c([3,2],.true.)
1139  c_loc(2) = c([2,2],.true.)
1140  ierr = add_arr_mult(&
1141  &var_int(:,:,:,c_loc(1),:,:,:),t_met,&
1142  &var_int(:,:,:,c_loc(2),&
1143  &derivs_loc(1,jd),derivs_loc(2,jd),&
1144  &derivs_loc(3,jd)),derivs_loc(:,jd))
1145  chckerr('')
1146  ! g_12 -> g_12 + T_met * g_11
1147  c_loc(1) = c([1,1],.true.)
1148  c_loc(2) = c([1,2],.true.)
1149  ierr = add_arr_mult(&
1150  &var_int(:,:,:,c_loc(1),:,:,:),t_met,&
1151  &var_int(:,:,:,c_loc(2),&
1152  &derivs_loc(1,jd),derivs_loc(2,jd),&
1153  &derivs_loc(3,jd)),derivs_loc(:,jd))
1154  chckerr('')
1155  end do
1156  end do
1157  else ! toroidal flux coordinates
1158  ierr = 1
1159  err_msg = '!!! INTERP_VAR_7D_REAL NOT YET &
1160  &IMPLEMENTED FOR TOROIDAL FLUX !!!'
1161  chckerr(err_msg)
1162  end if
1163  end if
1164 
1165  ! correct h_1j
1166  if (sym_type_loc.eq.4) then ! upper metric elements h_FD
1167  if (use_pol_flux_f) then ! poloidal flux coordinates
1168  do id = 0,max_deriv
1169  derivs_loc = derivs(id)
1170  do jd = 1,size(derivs_loc,2)
1171  ! h^11 -> h^11 - T_met * h^12
1172  c_loc(1) = c([1,2],.true.)
1173  c_loc(2) = c([1,1],.true.)
1174  ierr = add_arr_mult(&
1175  &var_int(:,:,:,c_loc(1),:,:,:),-t_met,&
1176  &var_int(:,:,:,c_loc(2),&
1177  &derivs_loc(1,jd),derivs_loc(2,jd),&
1178  &derivs_loc(3,jd)),derivs_loc(:,jd))
1179  chckerr('')
1180  ! h^12 -> h^12 - T_met * h^22
1181  c_loc(1) = c([2,2],.true.)
1182  c_loc(2) = c([1,2],.true.)
1183  ierr = add_arr_mult(&
1184  &var_int(:,:,:,c_loc(1),:,:,:),-t_met,&
1185  &var_int(:,:,:,c_loc(2),&
1186  &derivs_loc(1,jd),derivs_loc(2,jd),&
1187  &derivs_loc(3,jd)),derivs_loc(:,jd))
1188  chckerr('')
1189  ! h^11 -> h^11 - T_met * h^12
1190  c_loc(1) = c([1,2],.true.)
1191  c_loc(2) = c([1,1],.true.)
1192  ierr = add_arr_mult(&
1193  &var_int(:,:,:,c_loc(1),:,:,:),-t_met,&
1194  &var_int(:,:,:,c_loc(2),&
1195  &derivs_loc(1,jd),derivs_loc(2,jd),&
1196  &derivs_loc(3,jd)),derivs_loc(:,jd))
1197  chckerr('')
1198  ! h^13 -> h^13 - T_met * h^23
1199  c_loc(1) = c([2,3],.true.)
1200  c_loc(2) = c([1,3],.true.)
1201  ierr = add_arr_mult(&
1202  &var_int(:,:,:,c_loc(1),:,:,:),-t_met,&
1203  &var_int(:,:,:,c_loc(2),&
1204  &derivs_loc(1,jd),derivs_loc(2,jd),&
1205  &derivs_loc(3,jd)),derivs_loc(:,jd))
1206  chckerr('')
1207  end do
1208  end do
1209  else ! toroidal flux coordinates
1210  ierr = 1
1211  err_msg = '!!! INTERP_VAR_7D_REAL NOT YET &
1212  &IMPLEMENTED FOR TOROIDAL FLUX !!!'
1213  chckerr(err_msg)
1214  end if
1215  end if
1216  end function interp_var_7d_real
1218  integer function interp_var_7d_real_ow(var,sym_type) result(ierr) ! 7D_real version, overwriting variable
1219 
1220  ! input / output
1221  real(dp), intent(inout), allocatable :: var(:,:,:,:,:,:,:) ! variable to be interpolated
1222  integer, intent(in), optional :: sym_type ! optionally another type of symmetry
1223 
1224  ! local variables
1225  real(dp), allocatable :: var_7d_loc(:,:,:,:,:,:,:) ! local 7D variable
1226  integer :: var_dim(7,2) ! dimensions of variable
1227 
1228  ! initialize ierr
1229  ierr = 0
1230 
1231  ! set up variable dimensions
1232  var_dim(:,1) = [1,1,1,1,0,0,0]
1233  var_dim(:,2) = [shape(theta_i),size(var,4),size(var,5)-1,&
1234  &size(var,6)-1,size(var,7)-1]
1235 
1236  ! set up local 7D variable
1237  allocate(var_7d_loc(var_dim(1,1):var_dim(1,2),&
1238  &var_dim(2,1):var_dim(2,2),var_dim(3,1):var_dim(3,2),&
1239  &var_dim(4,1):var_dim(4,2),var_dim(5,1):var_dim(5,2),&
1240  &var_dim(6,1):var_dim(6,2),var_dim(7,1):var_dim(7,2)))
1241 
1242  ! call normal routine
1243  ierr = interp_var_7d_real(var,var_7d_loc,sym_type)
1244  chckerr('')
1245 
1246  ! reallocate variable
1247  deallocate(var)
1248  allocate(var(var_dim(1,1):var_dim(1,2),&
1249  &var_dim(2,1):var_dim(2,2),var_dim(3,1):var_dim(3,2),&
1250  &var_dim(4,1):var_dim(4,2),var_dim(5,1):var_dim(5,2),&
1251  &var_dim(6,1):var_dim(6,2),var_dim(7,1):var_dim(7,2)))
1252 
1253  ! overwrite variable
1254  var = var_7d_loc
1255  end function interp_var_7d_real_ow
1256  end function interp_hel_on_grid
1257 
1258 #if ldebug
1259 
1288  integer function test_metrics_h() result(ierr)
1289  use num_vars, only: rank, norm_disc_prec_eq
1290  use output_ops, only: plot_diff_hdf5
1291  use input_utilities, only: get_int
1292  use grid_vars, only: n_r_eq
1293  use num_utilities, only: spline
1294 
1295  character(*), parameter :: rout_name = 'test_metrics_H'
1296 
1297  ! local variables
1298  integer :: id, kd ! counters
1299  integer :: bcs(2,3) ! boundary conditions (theta(even), theta(odd), r)
1300  real(dp) :: bcs_val(2,3) ! values for boundary conditions
1301  real(dp), allocatable :: rchi(:,:), rpsi(:,:) ! chi and psi derivatives of R
1302  real(dp), allocatable :: zchi(:,:), zpsi(:,:) ! chi and psi derivatives of Z
1303  real(dp), allocatable :: jac(:,:) ! jac as defined above
1304  real(dp), allocatable :: jac_alt(:,:) ! alternative calculation for jac
1305  real(dp), allocatable :: h_h_11_alt(:,:,:) ! alternative calculation for upper metric factor 11
1306  real(dp), allocatable :: h_h_12_alt(:,:,:) ! alternative calculation for upper metric factor 12
1307  real(dp), allocatable :: h_h_33_alt(:,:,:) ! alternative calculation for upper metric factor 33
1308  character(len=max_str_ln) :: file_name ! name of plot file
1309  character(len=max_str_ln) :: description ! description of plot
1310  integer :: r_min = 4 ! first normal index that has meaning
1311  real(dp), allocatable :: tempvar(:,:,:,:) ! temporary variable
1312 
1313  ! initialize ierr
1314  ierr = 0
1315 
1316  if (rank.eq.0) then
1317  ! user output
1318  call writo('Checking consistency of metric factors')
1319  call lvl_ud(1)
1320 
1321  ! calculate the auxiliary quantities Zchi, zpsi, Rchi and Rpsi
1322  ! containing the derivatives as well as jac
1323  allocate(rchi(nchi,n_r_eq),rpsi(nchi,n_r_eq))
1324  allocate(zchi(nchi,n_r_eq),zpsi(nchi,n_r_eq))
1325  allocate(jac(nchi,n_r_eq))
1326 
1327  ! set up boundary conditions
1328  if (ias.eq.0) then ! top-bottom symmetric
1329  bcs(:,1) = [1,1] ! theta(even): zero first derivative
1330  bcs(:,2) = [2,2] ! theta(odd): zero first derivative
1331  else
1332  bcs(:,1) = [-1,-1] ! theta(even): periodic
1333  bcs(:,2) = [-1,-1] ! theta(odd): periodic
1334  end if
1335  bcs(:,3) = [0,0] ! r: not-a-knot
1336  bcs_val = 0._dp
1337 
1338  ! calculate derivatives
1339  do id = 1,nchi
1340  ierr = spline(flux_p_h(:,0)/(2*pi),r_h(id,:),&
1341  &flux_p_h(:,0)/(2*pi),rpsi(id,:),ord=norm_disc_prec_eq,&
1342  &deriv=1,bcs=bcs(:,3),bcs_val=bcs_val(:,3))
1343  chckerr('')
1344  ierr = spline(flux_p_h(:,0)/(2*pi),z_h(id,:),&
1345  &flux_p_h(:,0)/(2*pi),zpsi(id,:),ord=norm_disc_prec_eq,&
1346  &deriv=1,bcs=bcs(:,3),bcs_val=bcs_val(:,3))
1347  chckerr('')
1348  end do
1349  do kd = 1,n_r_eq
1350  ierr = spline(chi_h,r_h(:,kd),chi_h,rchi(:,kd),&
1351  &ord=norm_disc_prec_eq,deriv=1,bcs=bcs(:,1),&
1352  &bcs_val=bcs_val(:,1)) ! even
1353  chckerr('')
1354  ierr = spline(chi_h,z_h(:,kd),chi_h,zchi(:,kd),&
1355  &ord=norm_disc_prec_eq,deriv=1,bcs=bcs(:,2),&
1356  &bcs_val=bcs_val(:,2)) ! odd
1357  chckerr('')
1358  end do
1359 
1360  ! calculate Jacobian
1361  do kd = 1,n_r_eq
1362  jac(:,kd) = q_saf_h(kd,0)/(h_h_33(:,kd)*rbphi_h(kd,0))
1363  end do
1364 
1365  ! calculate Jacobian differently
1366  allocate(jac_alt(nchi,n_r_eq))
1367  jac_alt = r_h*(zchi*rpsi - zpsi*rchi)
1368 
1369  ! output jacobian
1370  ! set some variables
1371  file_name = 'TEST_jac_H'
1372  description = 'Testing whether the HELENA Jacobian is consistent.'
1373 
1374  ! plot difference
1375  call plot_diff_hdf5(&
1376  &reshape(jac(:,r_min:n_r_eq),&
1377  &[nchi,1,n_r_eq-r_min+1]),reshape(jac_alt(:,r_min:n_r_eq),&
1378  &[nchi,1,n_r_eq-r_min+1]),file_name,descr=description,&
1379  &output_message=.true.)
1380 
1381  ! calculate the metric factors directly
1382  allocate(h_h_11_alt(nchi,1,n_r_eq))
1383  allocate(h_h_12_alt(nchi,1,n_r_eq))
1384  allocate(h_h_33_alt(nchi,1,n_r_eq))
1385 
1386  h_h_11_alt(:,1,:) = (r_h/jac)**2 * (zchi**2 + rchi**2)
1387  h_h_12_alt(:,1,:) = -(r_h/jac)**2 * (zchi*zpsi + rchi*rpsi)
1388  h_h_33_alt(:,1,:) = 1._dp/(r_h**2)
1389 
1390  ! output h_H_11
1391  ! set some variables
1392  file_name = 'TEST_h_H_1_1'
1393  description = 'Testing whether the HELENA metric factor h_1_1 is &
1394  &consistent.'
1395 
1396  ! plot difference
1397  call plot_diff_hdf5(h_h_11_alt(:,:,r_min:n_r_eq),&
1398  &reshape(h_h_11(:,r_min:n_r_eq),[nchi,1,n_r_eq-r_min+1]),&
1399  &file_name,descr=description,output_message=.true.)
1400 
1401  ! output h_H_12
1402  ! set some variables
1403  file_name = 'TEST_h_H_1_2'
1404  description = 'Testing whether the HELENA metric factor h_1_2 is &
1405  &consistent.'
1406 
1407  ! plot difference
1408  call plot_diff_hdf5(h_h_12_alt(:,:,r_min:n_r_eq),&
1409  &reshape(h_h_12(:,r_min:n_r_eq),[nchi,1,n_r_eq-r_min+1]),&
1410  &file_name,descr=description,output_message=.true.)
1411 
1412  ! output h_H_33
1413  ! set some variables
1414  file_name = 'TEST_h_H_3_3'
1415  description = 'Testing whether the HELENA metric factor h_3_3 is &
1416  &consistent.'
1417 
1418  ! plot difference
1419  call plot_diff_hdf5(h_h_33_alt(:,:,r_min:n_r_eq),&
1420  &reshape(h_h_33(:,r_min:n_r_eq),[nchi,1,n_r_eq-r_min+1]),&
1421  &file_name,descr=description,output_message=.true.)
1422 
1423  ! user output
1424  call lvl_ud(-1)
1425  call writo('Test complete')
1426 
1427  ! user output
1428  call writo('Checking pressure balance')
1429  call lvl_ud(1)
1430 
1431  ! calculate auxiliary quantities:
1432  ! 1: D1 F ,
1433  ! 2: D1 p ,
1434  ! 3: D1 (q/F h_11) ,
1435  ! 4: D2 (q/F h_12) .
1436  allocate(tempvar(nchi,1,n_r_eq,4))
1437  do kd = 1,n_r_eq
1438  tempvar(:,1,kd,1) = rbphi_h(kd,1)
1439  tempvar(:,1,kd,2) = pres_h(kd,1)
1440  end do
1441  do id = 1,nchi
1442  ierr = spline(flux_p_h(:,0)/(2*pi),&
1443  &q_saf_h(:,0)/rbphi_h(:,0)*h_h_11(id,:),&
1444  &flux_p_h(:,0)/(2*pi),tempvar(id,1,:,3),&
1445  &ord=norm_disc_prec_eq,deriv=1,bcs=bcs(:,3),&
1446  &bcs_val=bcs_val(:,3))
1447  chckerr('')
1448  end do
1449  do kd = 1,n_r_eq
1450  ierr = spline(chi_h,q_saf_h(kd,0)/rbphi_h(kd,0)*h_h_12(:,kd),&
1451  &chi_h,tempvar(:,1,kd,4),ord=norm_disc_prec_eq,deriv=1,&
1452  &bcs=bcs(:,2),bcs_val=bcs_val(:,2)) ! odd
1453  chckerr('')
1454  end do
1455 
1456  ! calculate pressure balance in tempvar(1)
1457  ! mu_0 p' = F/(qR^2) (d/d1 (h_11 q/F) + d/d2 (h_12 q/F) + q F')
1458  do kd = 1,n_r_eq
1459  tempvar(:,1,kd,1) = &
1460  &-rbphi_h(kd,0)*h_h_33(:,kd)/q_saf_h(kd,0) * &
1461  &(tempvar(:,1,kd,1)*q_saf_h(kd,0) + tempvar(:,1,kd,3) + &
1462  &tempvar(:,1,kd,4))
1463  end do
1464 
1465  ! output difference with p'
1466  ! set some variables
1467  file_name = 'TEST_p_H'
1468  description = 'Testing whether the HELENA pressure balance is &
1469  &consistent'
1470 
1471  ! plot difference
1472  call plot_diff_hdf5(tempvar(:,:,:,1),tempvar(:,:,:,2),&
1473  &file_name,descr=description,output_message=.true.)
1474 
1475  ! user output
1476  call lvl_ud(-1)
1477  call writo('Test complete')
1478  end if
1479  end function test_metrics_h
1480 
1484  integer function test_harm_cont_h() result(ierr)
1485  use grid_vars, only: n_r_eq
1486  use grid_utilities, only: nufft
1487  use helena_vars, only: nchi, chi_h, ias, r_h, z_h
1488 
1489  character(*), parameter :: rout_name = 'test_metrics_H'
1490 
1491  ! local variables
1492  integer :: kd ! counter
1493  integer :: nchi_full ! nr. of pol. points in full grid, without doubling
1494  real(dp), allocatable :: f_loc(:,:) ! local Fourier modes
1495  real(dp), allocatable :: chi_full(:) ! chi_H in full poloidal grid
1496  real(dp), allocatable :: r_h_full(:,:) ! R_H in full poloidal grid
1497  real(dp), allocatable :: z_h_full(:,:) ! Z_H in full poloidal grid
1498  real(dp), allocatable :: r_h_f(:,:,:,:) ! Fourier modes of R_H
1499  real(dp), allocatable :: z_h_f(:,:,:,:) ! Fourier modes of Z_H
1500  character(len=max_str_ln) :: plot_name ! name of plot
1501  character(len=2) :: loc_str(2) ! local string
1502 
1503  ! initialize ierr
1504  ierr = 0
1505 
1506  ! set up full R and Z
1507  if (ias.eq.0) then
1508  nchi_full = 2*(nchi-1)
1509  allocate(chi_full(nchi_full))
1510  allocate(r_h_full(nchi_full,n_r_eq))
1511  allocate(z_h_full(nchi_full,n_r_eq))
1512  chi_full(1:nchi-1) = chi_h(1:nchi-1)
1513  chi_full(2*(nchi-1):nchi:-1) = 2*pi - chi_h(2:nchi)
1514  r_h_full(1:nchi-1,:) = r_h(1:nchi-1,:)
1515  r_h_full(2*(nchi-1):nchi:-1,:) = r_h(2:nchi,:)
1516  z_h_full(1:nchi-1,:) = z_h(1:nchi-1,:)
1517  z_h_full(2*(nchi-1):nchi:-1,:) = -z_h(2:nchi,:)
1518  else
1519  nchi_full = nchi-1
1520  allocate(chi_full(nchi_full))
1521  allocate(r_h_full(nchi_full,n_r_eq))
1522  allocate(z_h_full(nchi_full,n_r_eq))
1523  chi_full(1:nchi_full) = chi_h(1:nchi_full)
1524  r_h_full(1:nchi_full,:) = r_h(1:nchi_full,:)
1525  z_h_full(1:nchi_full,:) = z_h(1:nchi_full,:)
1526  end if
1527 
1528  ! calculate NUFFT
1529  do kd = 1,n_r_eq
1530  ierr = nufft(chi_full,r_h_full(:,kd),f_loc)
1531  chckerr('')
1532  if (kd.eq.1) allocate(r_h_f(size(f_loc,1),1,n_r_eq,2))
1533  r_h_f(:,1,kd,:) = f_loc
1534  ierr = nufft(chi_full,z_h_full(:,kd),f_loc)
1535  chckerr('')
1536  if (kd.eq.1) allocate(z_h_f(size(f_loc,1),1,n_r_eq,2))
1537  z_h_f(:,1,kd,:) = f_loc
1538  end do
1539 
1540  ! plot
1541  loc_str(1) = 'Re'
1542  loc_str(2) = 'Im'
1543  do kd = 1,2
1544  plot_name = loc_str(kd)//'_R_H_F'
1545  call print_ex_2d(['1'],plot_name,r_h_f(:,1,:,kd),draw=.false.)
1546  call draw_ex(['1'],plot_name,n_r_eq,1,.false.)
1547  plot_name = loc_str(kd)//'_Z_H_F'
1548  call print_ex_2d(['1'],plot_name,z_h_f(:,1,:,kd),draw=.false.)
1549  call draw_ex(['1'],plot_name,n_r_eq,1,.false.)
1550 
1551  plot_name = loc_str(kd)//'_R_H_F_log'
1552  call print_ex_2d(['1'],plot_name,&
1553  &log10(max(1.e-10_dp,abs(r_h_f(:,1,:,kd)))),&
1554  &draw=.false.)
1555  call draw_ex(['1'],plot_name,n_r_eq,1,.false.)
1556  plot_name = loc_str(kd)//'_Z_H_F_log'
1557  call print_ex_2d(['1'],plot_name,&
1558  &log10(max(1.e-10_dp,abs(z_h_f(:,1,:,kd)))),&
1559  &draw=.false.)
1560  call draw_ex(['1'],plot_name,n_r_eq,1,.false.)
1561  end do
1562  call plot_hdf5(['real','imag'],'R_H_F',r_h_f)
1563  call plot_hdf5(['real','imag'],'R_H_F_log',&
1564  &log10(max(1.e-10_dp,abs(r_h_f))))
1565  call plot_hdf5(['real','imag'],'Z_H_F',z_h_f)
1566  call plot_hdf5(['real','imag'],'Z_H_F_log',&
1567  &log10(max(1.e-10_dp,abs(z_h_f))))
1568  end function
1569 #endif
1570 end module helena_ops
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
output_ops::plot_diff_hdf5
subroutine, public plot_diff_hdf5(A, B, file_name, tot_dim, loc_offset, descr, output_message)
Takes two input vectors and plots these as well as the relative and absolute difference in a HDF5 fil...
Definition: output_ops.f90:1765
num_utilities::calc_int
Integrates a function using the trapezoidal rule.
Definition: num_utilities.f90:160
num_vars::invert_top_bottom_h
logical, public invert_top_bottom_h
invert top and bottom for HELENA equilibria
Definition: num_vars.f90:144
helena_vars::rbphi_h
real(dp), dimension(:,:), allocatable, public rbphi_h
Definition: HELENA_vars.f90:29
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
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
num_utilities::con2dis
Convert between points from a continuous grid to a discrete grid.
Definition: num_utilities.f90:188
helena_vars::rot_t_h
real(dp), dimension(:,:), allocatable, public rot_t_h
rotational transform
Definition: HELENA_vars.f90:28
num_vars::norm_disc_prec_eq
integer, public norm_disc_prec_eq
precision for normal discretization for equilibrium
Definition: num_vars.f90:120
str_utilities::i2str
elemental character(len=max_str_ln) function, public i2str(k)
Convert an integer to string.
Definition: str_utilities.f90:18
helena_vars::q_saf_h
real(dp), dimension(:,:), allocatable, public q_saf_h
safety factor
Definition: HELENA_vars.f90:27
helena_vars::flux_p_h
real(dp), dimension(:,:), allocatable, public flux_p_h
poloidal flux
Definition: HELENA_vars.f90:24
grid_vars::n_r_eq
integer, public n_r_eq
nr. of normal points in equilibrium grid
Definition: grid_vars.f90:20
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
num_utilities::calc_mult
Calculate matrix multiplication of two square matrices .
Definition: num_utilities.f90:95
helena_ops::get_ang_interp_data_hel
integer function get_ang_interp_data_hel(grid_in, grid_out, theta_i, fund_theta_int_displ, tb_sym, use_E)
Calculate interpolation factors for angular interpolation in grid_out of quantities defined on grid_i...
Definition: HELENA_ops.f90:413
eq_vars::eq_1_type
flux equilibrium type
Definition: eq_vars.f90:63
helena_ops::read_hel
integer function, public read_hel(n_r_in, use_pol_flux_H)
Reads the HELENA equilibrium data.
Definition: HELENA_ops.f90:84
helena_vars::h_h_33
real(dp), dimension(:,:), allocatable, public h_h_33
upper metric factor (1 / gem12)
Definition: HELENA_vars.f90:32
helena_vars::r_h
real(dp), dimension(:,:), allocatable, public r_h
major radius (xout)
Definition: HELENA_vars.f90:33
helena_vars::flux_t_h
real(dp), dimension(:,:), allocatable, public flux_t_h
toroidal flux
Definition: HELENA_vars.f90:25
helena_vars::h_h_12
real(dp), dimension(:,:), allocatable, public h_h_12
upper metric factor (gem12)
Definition: HELENA_vars.f90:31
helena_ops::test_harm_cont_h
integer function, public test_harm_cont_h()
Investaige harmonic content of the HELENA variables.
Definition: HELENA_ops.f90:1485
helena_vars::z_h
real(dp), dimension(:,:), allocatable, public z_h
height (yout)
Definition: HELENA_vars.f90:34
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_utilities::derivs
integer function, dimension(:,:), allocatable, public derivs(order, dims)
Returns derivatives of certain order.
Definition: num_utilities.f90:2246
helena_vars::pres_h
real(dp), dimension(:,:), allocatable, public pres_h
pressure profile
Definition: HELENA_vars.f90:26
x_vars::x_2_type
tensorial perturbation type
Definition: X_vars.f90:81
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
x_vars
Variables pertaining to the perturbation quantities.
Definition: X_vars.f90:4
helena_ops::interp_hel_on_grid
integer function, public interp_hel_on_grid(grid_in, grid_out, eq_2, X_1, X_2, eq_2_out, X_1_out, X_2_out, eq_1, grid_name)
Interpolate variables resulting from HELENA equilibria to another grid (angularly).
Definition: HELENA_ops.f90:550
helena_ops
Operations on HELENA variables.
Definition: HELENA_ops.f90:4
num_vars::tol_zero
real(dp), public tol_zero
tolerance for zeros
Definition: num_vars.f90:133
helena_vars::h_h_11
real(dp), dimension(:,:), allocatable, public h_h_11
upper metric factor (gem11)
Definition: HELENA_vars.f90:30
num_vars::max_deriv
integer, parameter, public max_deriv
highest derivatives for metric factors in Flux coords.
Definition: num_vars.f90:52
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
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
helena_vars::bmtog_h
real(dp), public bmtog_h
B_geo/B_mag.
Definition: HELENA_vars.f90:22
grid_vars
Variables pertaining to the different grids used.
Definition: grid_vars.f90:4
helena_vars::ias
integer, public ias
0 if top-bottom symmetric, 1 if not
Definition: HELENA_vars.f90:36
num_vars::eq_i
integer, parameter, public eq_i
file number of equilibrium file from VMEC or HELENA
Definition: num_vars.f90:183
messages::lvl_ud
subroutine, public lvl_ud(inc)
Increases/decreases lvl of output.
Definition: messages.f90:254
input_utilities::get_int
integer function, public get_int(lim_lo, lim_hi, ind)
Queries for user input for an integer value, where allowable range can be provided as well.
Definition: input_utilities.f90:152
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
helena_ops::test_metrics_h
integer function, public test_metrics_h()
Checks whether the metric elements provided by HELENA are consistent with a direct calculation using ...
Definition: HELENA_ops.f90:1289
x_vars::x_1_type
vectorial perturbation type
Definition: X_vars.f90:51
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::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
input_utilities
Numerical utilities related to input.
Definition: input_utilities.f90:4
num_vars::eq_name
character(len=max_str_ln), public eq_name
name of equilibrium file from VMEC or HELENA
Definition: num_vars.f90:138
num_utilities::add_arr_mult
Add to an array (3) the product of arrays (1) and (2).
Definition: num_utilities.f90:39
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
helena_vars::rmtog_h
real(dp), public rmtog_h
R_geo/R_mag.
Definition: HELENA_vars.f90:21
this documentation was created with Doxygen 1.8.17 by Toon Weyens