PB3D  [2.45]
Ideal linear high-n MHD stability in 3-D
X_ops.f90
Go to the documentation of this file.
1 !------------------------------------------------------------------------------!
3 !------------------------------------------------------------------------------!
4 module x_ops
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, iu, max_str_ln, max_name_ln, pi
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, modes_type
14 
15  implicit none
16  private
17  public calc_x, calc_magn_ints, print_output_x, resonance_plot, &
18  &calc_res_surf, check_x_modes, init_modes, setup_modes, divide_x_jobs, &
19  &redistribute_output_x
20 #if ldebug
21  public debug_check_x_modes_2, print_debug_x_1, print_debug_x_2
22 #endif
23 
24  ! global variables
25 #if ldebug
26  logical :: debug_check_x_modes_2 = .false.
27  logical :: debug_setup_modes = .false.
28 #endif
29 
30  ! interfaces
31 
39  interface calc_x
41  module procedure calc_x_1
43  module procedure calc_x_2
44  end interface
45 
51  interface redistribute_output_x
53  module procedure redistribute_output_x_1
55  module procedure redistribute_output_x_2
56  end interface
57 
85  interface print_output_x
87  module procedure print_output_x_1
89  module procedure print_output_x_2
90  end interface
91 
92 contains
94  integer function calc_x_1(mds,grid_eq,grid_X,eq_1,eq_2,X,lim_sec_X) &
95  &result(ierr)
96 
97  character(*), parameter :: rout_name = 'calc_X_1'
98 
99  ! input / output
100  type(modes_type), intent(in) :: mds
101  type(grid_type), intent(in) :: grid_eq
102  type(grid_type), intent(in) :: grid_x
103  type(eq_1_type), intent(in) :: eq_1
104  type(eq_2_type), intent(in) :: eq_2
105  type(x_1_type), intent(inout) :: x
106  integer, intent(in), optional :: lim_sec_x(2)
107 
108  ! initialize ierr
109  ierr = 0
110 
111  ! user output
112  call writo('Calculate vectorial perturbation variables')
113  call lvl_ud(1)
114 
115  ! create perturbation with modes of current X job
116  call x%init(mds,grid_x,lim_sec_x)
117 
118  ! calculate U and DU
119  call writo('Calculating U and DU...')
120  call lvl_ud(1)
121  ierr = calc_u(grid_eq,grid_x,eq_1,eq_2,x)
122  chckerr('')
123  call lvl_ud(-1)
124 
125  ! user output
126  call lvl_ud(-1)
127  end function calc_x_1
129  integer function calc_x_2(mds,grid_eq,grid_X,eq_1,eq_2,X_a,X_b,X,&
130  &lim_sec_X) result(ierr)
131 
132  character(*), parameter :: rout_name = 'calc_X_2'
133 
134  ! input / output
135  type(modes_type), intent(in) :: mds
136  type(grid_type), intent(in) :: grid_eq
137  type(grid_type), intent(in) :: grid_x
138  type(eq_1_type), intent(in) :: eq_1
139  type(eq_2_type), intent(in) :: eq_2
140  type(x_1_type), intent(inout) :: x_a
141  type(x_1_type), intent(inout) :: x_b
142  type(x_2_type), intent(inout) :: x
143  integer, intent(in), optional :: lim_sec_x(2,2)
144 
145  ! initialize ierr
146  ierr = 0
147 
148  ! user output
149  call writo('Calculate tensorial perturbation variables')
150  call lvl_ud(1)
151 
152  ! create perturbation with modes of current X job
153  call x%init(mds,grid_x,lim_sec_x)
154 
155  ! Calculate PV_i for all (k,m) pairs and n_r (equilibrium)
156  ! values of the normal coordinate
157  call writo('Calculating PV...')
158  call lvl_ud(1)
159  ierr = calc_pv(grid_eq,grid_x,eq_1,eq_2,x_a,x_b,x,lim_sec_x)
160  chckerr('')
161  call lvl_ud(-1)
162 
163  ! Calculate KV_i for all (k,m) pairs and n_r (equilibrium)
164  ! values of the normal coordinate
165  call writo('Calculating KV...')
166  call lvl_ud(1)
167  ierr = calc_kv(grid_eq,grid_x,eq_1,eq_2,x_a,x_b,x,lim_sec_x)
168  chckerr('')
169  call lvl_ud(-1)
170 
171  ! user output
172  call lvl_ud(-1)
173  end function calc_x_2
174 
176  integer function redistribute_output_x_1(mds,grid,grid_out,X,X_out) &
177  &result(ierr)
178  use mpi_utilities, only: redistribute_var
179 
180  character(*), parameter :: rout_name = 'redistribute_output_X_1'
181 
182  ! input / output
183  type(modes_type), intent(in) :: mds
184  type(grid_type), intent(in) :: grid
185  type(grid_type), intent(in) :: grid_out
186  type(x_1_type), intent(in) :: x
187  type(x_1_type), intent(inout) :: x_out
188 
189  ! local variables
190  integer :: ld ! counter
191  integer :: lims(2), lims_dis(2) ! limits and distributed limits, taking into account the angular extent
192  integer :: siz(3), siz_dis(3) ! size for geometric part of variable
193  real(dp), allocatable :: temp_var(:) ! temporary variable
194 
195  ! initialize ierr
196  ierr = 0
197 
198  ! user output
199  call writo('Redistribute vectorial perturbation variables')
200  call lvl_ud(1)
201 
202  ! test
203  if (grid%n(1).ne.grid_out%n(1) .or. grid%n(2).ne.grid_out%n(2)) then
204  ierr = 1
205  chckerr('grid and grid_out are not compatible')
206  end if
207 
208  ! create redistributed vectorial perturbation variables
209  call x_out%init(mds,grid_out,lim_sec_x=x%lim_sec_X)
210 
211  ! set up limits taking into account angular extent and temporary var
212  lims(1) = product(grid%n(1:2))*(grid%i_min-1)+1
213  lims(2) = product(grid%n(1:2))*grid%i_max
214  lims_dis(1) = product(grid%n(1:2))*(grid_out%i_min-1)+1
215  lims_dis(2) = product(grid%n(1:2))*grid_out%i_max
216  siz = [grid%n(1:2),grid%loc_n_r]
217  siz_dis = [grid%n(1:2),grid_out%loc_n_r]
218  allocate(temp_var(product(siz_dis)))
219 
220  ! for all mode numbers
221  do ld = 1,size(x%U_0,4)
222  ! U_0
223  ierr = redistribute_var(reshape(rp(x%U_0(:,:,:,ld)),&
224  &[product(siz)]),temp_var,lims,lims_dis)
225  chckerr('')
226  x_out%U_0(:,:,:,ld) = &
227  &reshape(temp_var(1:product(siz_dis)),siz_dis)
228  ierr = redistribute_var(reshape(ip(x%U_0(:,:,:,ld)),&
229  &[product(siz)]),temp_var,lims,lims_dis)
230  chckerr('')
231  x_out%U_0(:,:,:,ld) = x_out%U_0(:,:,:,ld) + iu*&
232  &reshape(temp_var(1:product(siz_dis)),siz_dis)
233 
234  ! U_1
235  ierr = redistribute_var(reshape(rp(x%U_1(:,:,:,ld)),&
236  &[product(siz)]),temp_var,lims,lims_dis)
237  chckerr('')
238  x_out%U_1(:,:,:,ld) = &
239  &reshape(temp_var(1:product(siz_dis)),siz_dis)
240  ierr = redistribute_var(reshape(ip(x%U_1(:,:,:,ld)),&
241  &[product(siz)]),temp_var,lims,lims_dis)
242  chckerr('')
243  x_out%U_1(:,:,:,ld) = x_out%U_1(:,:,:,ld) + iu*&
244  &reshape(temp_var(1:product(siz_dis)),siz_dis)
245 
246  ! DU_0
247  ierr = redistribute_var(reshape(rp(x%DU_0(:,:,:,ld)),&
248  &[product(siz)]),temp_var,lims,lims_dis)
249  chckerr('')
250  x_out%DU_0(:,:,:,ld) = &
251  &reshape(temp_var(1:product(siz_dis)),siz_dis)
252  ierr = redistribute_var(reshape(ip(x%DU_0(:,:,:,ld)),&
253  &[product(siz)]),temp_var,lims,lims_dis)
254  chckerr('')
255  x_out%DU_0(:,:,:,ld) = x_out%DU_0(:,:,:,ld) + iu*&
256  &reshape(temp_var(1:product(siz_dis)),siz_dis)
257 
258  ! DU_1
259  ierr = redistribute_var(reshape(rp(x%DU_1(:,:,:,ld)),&
260  &[product(siz)]),temp_var,lims,lims_dis)
261  chckerr('')
262  x_out%DU_1(:,:,:,ld) = &
263  &reshape(temp_var(1:product(siz_dis)),siz_dis)
264  ierr = redistribute_var(reshape(ip(x%DU_1(:,:,:,ld)),&
265  &[product(siz)]),temp_var,lims,lims_dis)
266  chckerr('')
267  x_out%DU_1(:,:,:,ld) = x_out%DU_1(:,:,:,ld) + iu*&
268  &reshape(temp_var(1:product(siz_dis)),siz_dis)
269  end do
270 
271  ! user output
272  call lvl_ud(-1)
273  end function redistribute_output_x_1
275  integer function redistribute_output_x_2(mds,grid,grid_out,X,X_out) &
276  &result(ierr)
277  use mpi_utilities, only: redistribute_var
278 
279  character(*), parameter :: rout_name = 'redistribute_output_X_2'
280 
281  ! input / output
282  type(modes_type), intent(in) :: mds
283  type(grid_type), intent(in) :: grid
284  type(grid_type), intent(in) :: grid_out
285  type(x_2_type), intent(in) :: x
286  type(x_2_type), intent(inout) :: x_out
287 
288  ! local variables
289  integer :: ld ! counters
290  integer :: n_loc(2) ! local n
291  integer :: lims(2), lims_dis(2) ! limits and distributed limits, taking into account the angular extent
292  integer :: siz(3), siz_dis(3) ! size for geometric part of variable
293  logical :: is_field_averaged ! whether field-averaged, so only one dimension for first index
294  real(dp), allocatable :: temp_var(:) ! temporary variable
295 
296  ! initialize ierr
297  ierr = 0
298 
299  ! user output
300  call writo('Redistribute tensorial perturbation variables')
301  call lvl_ud(1)
302 
303  ! test
304  is_field_averaged = size(x%PV_0,1).eq.1
305  if ((grid%n(1).ne.grid_out%n(1) .or. grid%n(2).ne.grid_out%n(2)) .and. &
306  &.not.is_field_averaged) then ! for field-averaged grids, don't do tests
307  ierr = 1
308  chckerr('grid and grid_out are not compatible')
309  end if
310 
311  ! create redistributed tensorial perturbation variables
312  call x_out%init(mds,grid_out,lim_sec_x=x%lim_sec_X,&
313  &is_field_averaged=is_field_averaged)
314 
315  ! set up limits taking into account angular extent and temporary var
316  n_loc = grid%n(1:2)
317  if (is_field_averaged) n_loc(1) = 1
318  lims(1) = product(n_loc)*(grid%i_min-1)+1
319  lims(2) = product(n_loc)*grid%i_max
320  lims_dis(1) = product(n_loc)*(grid_out%i_min-1)+1
321  lims_dis(2) = product(n_loc)*grid_out%i_max
322  siz = [n_loc,grid%loc_n_r]
323  siz_dis = [n_loc,grid_out%loc_n_r]
324  allocate(temp_var(product(siz_dis)))
325 
326  ! for all mode number combinations, symmetric
327  do ld = 1,size(x%PV_0,4)
328  ! PV_0
329  ierr = redistribute_var(reshape(rp(x%PV_0(1:n_loc(1),:,:,ld)),&
330  &[product(siz)]),temp_var,lims,lims_dis)
331  chckerr('')
332  x_out%PV_0(1:n_loc(1),:,:,ld) = &
333  &reshape(temp_var(1:product(siz_dis)),siz_dis)
334  ierr = redistribute_var(reshape(ip(x%PV_0(1:n_loc(1),:,:,ld)),&
335  &[product(siz)]),temp_var,lims,lims_dis)
336  chckerr('')
337  x_out%PV_0(1:n_loc(1),:,:,ld) = x_out%PV_0(1:n_loc(1),:,:,ld) + iu*&
338  &reshape(temp_var(1:product(siz_dis)),siz_dis)
339 
340  ! PV_2
341  ierr = redistribute_var(reshape(rp(x%PV_2(1:n_loc(1),:,:,ld)),&
342  &[product(siz)]),temp_var,lims,lims_dis)
343  chckerr('')
344  x_out%PV_2(1:n_loc(1),:,:,ld) = &
345  &reshape(temp_var(1:product(siz_dis)),siz_dis)
346  ierr = redistribute_var(reshape(ip(x%PV_2(1:n_loc(1),:,:,ld)),&
347  &[product(siz)]),temp_var,lims,lims_dis)
348  chckerr('')
349  x_out%PV_2(1:n_loc(1),:,:,ld) = x_out%PV_2(1:n_loc(1),:,:,ld) + iu*&
350  &reshape(temp_var(1:product(siz_dis)),siz_dis)
351 
352  ! KV_0
353  ierr = redistribute_var(reshape(rp(x%KV_0(1:n_loc(1),:,:,ld)),&
354  &[product(siz)]),temp_var,lims,lims_dis)
355  chckerr('')
356  x_out%KV_0(1:n_loc(1),:,:,ld) = &
357  &reshape(temp_var(1:product(siz_dis)),siz_dis)
358  ierr = redistribute_var(reshape(ip(x%KV_0(1:n_loc(1),:,:,ld)),&
359  &[product(siz)]),temp_var,lims,lims_dis)
360  chckerr('')
361  x_out%KV_0(1:n_loc(1),:,:,ld) = x_out%KV_0(1:n_loc(1),:,:,ld) + iu*&
362  &reshape(temp_var(1:product(siz_dis)),siz_dis)
363 
364  ! KV_2
365  ierr = redistribute_var(reshape(rp(x%KV_2(1:n_loc(1),:,:,ld)),&
366  &[product(siz)]),temp_var,lims,lims_dis)
367  chckerr('')
368  x_out%KV_2(1:n_loc(1),:,:,ld) = &
369  &reshape(temp_var(1:product(siz_dis)),siz_dis)
370  ierr = redistribute_var(reshape(ip(x%KV_2(1:n_loc(1),:,:,ld)),&
371  &[product(siz)]),temp_var,lims,lims_dis)
372  chckerr('')
373  x_out%KV_2(1:n_loc(1),:,:,ld) = x_out%KV_2(1:n_loc(1),:,:,ld) + iu*&
374  &reshape(temp_var(1:product(siz_dis)),siz_dis)
375  end do
376 
377  ! for all mode number combinations, asymmetric
378  do ld = 1,size(x%PV_1,4)
379  ! PV_1
380  ierr = redistribute_var(reshape(rp(x%PV_1(1:n_loc(1),:,:,ld)),&
381  &[product(siz)]),temp_var,lims,lims_dis)
382  chckerr('')
383  x_out%PV_1(1:n_loc(1),:,:,ld) = &
384  &reshape(temp_var(1:product(siz_dis)),siz_dis)
385  ierr = redistribute_var(reshape(ip(x%PV_1(1:n_loc(1),:,:,ld)),&
386  &[product(siz)]),temp_var,lims,lims_dis)
387  chckerr('')
388  x_out%PV_1(1:n_loc(1),:,:,ld) = x_out%PV_1(1:n_loc(1),:,:,ld) + iu*&
389  &reshape(temp_var(1:product(siz_dis)),siz_dis)
390 
391  ! KV_1
392  ierr = redistribute_var(reshape(rp(x%KV_1(1:n_loc(1),:,:,ld)),&
393  &[product(siz)]),temp_var,lims,lims_dis)
394  chckerr('')
395  x_out%KV_1(1:n_loc(1),:,:,ld) = &
396  &reshape(temp_var(1:product(siz_dis)),siz_dis)
397  ierr = redistribute_var(reshape(ip(x%KV_1(1:n_loc(1),:,:,ld)),&
398  &[product(siz)]),temp_var,lims,lims_dis)
399  chckerr('')
400  x_out%KV_1(1:n_loc(1),:,:,ld) = x_out%KV_1(1:n_loc(1),:,:,ld) + iu*&
401  &reshape(temp_var(1:product(siz_dis)),siz_dis)
402  end do
403 
404  ! user output
405  call lvl_ud(-1)
406  end function redistribute_output_x_2
407 
409  integer function print_output_x_1(grid_X,X,data_name,rich_lvl,par_div,&
410  &lim_sec_X) result(ierr)
411 
412  use num_vars, only: pb3d_name, eq_jobs_lims, eq_job_nr
413  use grid_utilities, only: trim_grid
414  use hdf5_ops, only: print_hdf5_arrs
415  use hdf5_vars, only: dealloc_var_1d, var_1d_type, &
416  &max_dim_var_1d
417  use x_vars, only: n_mod_x
418 
419  character(*), parameter :: rout_name = 'print_output_X_1'
420 
421  ! input / output
422  type(grid_type), intent(in) :: grid_x
423  type(x_1_type), intent(in) :: x
424  character(len=*), intent(in) :: data_name
425  integer, intent(in), optional :: rich_lvl
426  logical, intent(in), optional :: par_div
427  integer, intent(in), optional :: lim_sec_x(2)
428 
429  ! local variables
430  type(grid_type) :: grid_trim ! trimmed grid
431  integer :: n_tot(3) ! total n
432  integer :: par_id(2) ! parallel interval
433  integer :: norm_id(2) ! untrimmed normal indices for trimmed grids
434  type(var_1d_type), allocatable, target :: x_1d(:) ! 1D equivalent of X variables
435  type(var_1d_type), pointer :: x_1d_loc => null() ! local element in X_1D
436  integer :: id ! counters
437  logical :: par_div_loc ! local par_div
438  integer :: lim_sec_x_loc(2) ! local lim_sec_X
439  integer :: loc_size ! local size
440 
441  ! initialize ierr
442  ierr = 0
443 
444  ! user output
445  call writo('Write vectorial perturbation variables to output file')
446  call lvl_ud(1)
447 
448  ! trim grid
449  ierr = trim_grid(grid_x,grid_trim,norm_id)
450  chckerr('')
451 
452  ! set local par_div
453  par_div_loc = .false.
454  if (present(par_div)) par_div_loc = par_div
455 
456  ! set total n and parallel interval
457  n_tot = grid_trim%n
458  par_id = [1,n_tot(1)]
459  if (grid_trim%n(1).gt.0 .and. par_div_loc) then ! total grid includes all equilibrium jobs
460  n_tot(1) = maxval(eq_jobs_lims)-minval(eq_jobs_lims)+1
461  par_id = eq_jobs_lims(:,eq_job_nr)
462  end if
463 
464  ! set local size and lim_sec_X
465  loc_size = size(x%U_0(:,:,norm_id(1):norm_id(2),:))
466  lim_sec_x_loc = [1,n_mod_x]
467  if (present(lim_sec_x)) lim_sec_x_loc = lim_sec_x
468 
469  ! Set up the 1D equivalents of the perturbation variables
470  allocate(x_1d(max_dim_var_1d))
471 
472  ! set up variables X_1D
473  id = 1
474 
475  ! RE_U_0
476  x_1d_loc => x_1d(id); id = id+1
477  x_1d_loc%var_name = 'RE_U_0'
478  allocate(x_1d_loc%tot_i_min(4),x_1d_loc%tot_i_max(4))
479  allocate(x_1d_loc%loc_i_min(4),x_1d_loc%loc_i_max(4))
480  x_1d_loc%tot_i_min = [1,1,1,1]
481  x_1d_loc%tot_i_max = [n_tot,n_mod_x]
482  x_1d_loc%loc_i_min = [par_id(1),1,grid_trim%i_min,lim_sec_x_loc(1)]
483  x_1d_loc%loc_i_max = [par_id(2),n_tot(2),grid_trim%i_max,&
484  &lim_sec_x_loc(2)]
485  allocate(x_1d_loc%p(loc_size))
486  x_1d_loc%p = reshape(rp(x%U_0(:,:,norm_id(1):norm_id(2),:)),&
487  &[loc_size])
488 
489  ! IM_U_0
490  x_1d_loc => x_1d(id); id = id+1
491  x_1d_loc%var_name = 'IM_U_0'
492  allocate(x_1d_loc%tot_i_min(4),x_1d_loc%tot_i_max(4))
493  allocate(x_1d_loc%loc_i_min(4),x_1d_loc%loc_i_max(4))
494  x_1d_loc%tot_i_min = [1,1,1,1]
495  x_1d_loc%tot_i_max = [n_tot,n_mod_x]
496  x_1d_loc%loc_i_min = [par_id(1),1,grid_trim%i_min,lim_sec_x_loc(1)]
497  x_1d_loc%loc_i_max = [par_id(2),n_tot(2),grid_trim%i_max,&
498  &lim_sec_x_loc(2)]
499  allocate(x_1d_loc%p(loc_size))
500  x_1d_loc%p = reshape(ip(x%U_0(:,:,norm_id(1):norm_id(2),:)),&
501  &[loc_size])
502 
503  ! RE_U_1
504  x_1d_loc => x_1d(id); id = id+1
505  x_1d_loc%var_name = 'RE_U_1'
506  allocate(x_1d_loc%tot_i_min(4),x_1d_loc%tot_i_max(4))
507  allocate(x_1d_loc%loc_i_min(4),x_1d_loc%loc_i_max(4))
508  x_1d_loc%tot_i_min = [1,1,1,1]
509  x_1d_loc%tot_i_max = [n_tot,n_mod_x]
510  x_1d_loc%loc_i_min = [par_id(1),1,grid_trim%i_min,lim_sec_x_loc(1)]
511  x_1d_loc%loc_i_max = [par_id(2),n_tot(2),grid_trim%i_max,&
512  &lim_sec_x_loc(2)]
513  allocate(x_1d_loc%p(loc_size))
514  x_1d_loc%p = reshape(rp(x%U_1(:,:,norm_id(1):norm_id(2),:)),&
515  &[loc_size])
516 
517  ! IM_U_1
518  x_1d_loc => x_1d(id); id = id+1
519  x_1d_loc%var_name = 'IM_U_1'
520  allocate(x_1d_loc%tot_i_min(4),x_1d_loc%tot_i_max(4))
521  allocate(x_1d_loc%loc_i_min(4),x_1d_loc%loc_i_max(4))
522  x_1d_loc%tot_i_min = [1,1,1,1]
523  x_1d_loc%tot_i_max = [n_tot,n_mod_x]
524  x_1d_loc%loc_i_min = [par_id(1),1,grid_trim%i_min,lim_sec_x_loc(1)]
525  x_1d_loc%loc_i_max = [par_id(2),n_tot(2),grid_trim%i_max,&
526  &lim_sec_x_loc(2)]
527  allocate(x_1d_loc%p(loc_size))
528  x_1d_loc%p = reshape(ip(x%U_1(:,:,norm_id(1):norm_id(2),:)),&
529  &[loc_size])
530 
531  ! RE_DU_0
532  x_1d_loc => x_1d(id); id = id+1
533  x_1d_loc%var_name = 'RE_DU_0'
534  allocate(x_1d_loc%tot_i_min(4),x_1d_loc%tot_i_max(4))
535  allocate(x_1d_loc%loc_i_min(4),x_1d_loc%loc_i_max(4))
536  x_1d_loc%tot_i_min = [1,1,1,1]
537  x_1d_loc%tot_i_max = [n_tot,n_mod_x]
538  x_1d_loc%loc_i_min = [par_id(1),1,grid_trim%i_min,lim_sec_x_loc(1)]
539  x_1d_loc%loc_i_max = [par_id(2),n_tot(2),grid_trim%i_max,&
540  &lim_sec_x_loc(2)]
541  allocate(x_1d_loc%p(loc_size))
542  x_1d_loc%p = reshape(rp(x%DU_0(:,:,norm_id(1):norm_id(2),:)),&
543  &[loc_size])
544 
545  ! IM_DU_0
546  x_1d_loc => x_1d(id); id = id+1
547  x_1d_loc%var_name = 'IM_DU_0'
548  allocate(x_1d_loc%tot_i_min(4),x_1d_loc%tot_i_max(4))
549  allocate(x_1d_loc%loc_i_min(4),x_1d_loc%loc_i_max(4))
550  x_1d_loc%tot_i_min = [1,1,1,1]
551  x_1d_loc%tot_i_max = [n_tot,n_mod_x]
552  x_1d_loc%loc_i_min = [par_id(1),1,grid_trim%i_min,lim_sec_x_loc(1)]
553  x_1d_loc%loc_i_max = [par_id(2),n_tot(2),grid_trim%i_max,&
554  &lim_sec_x_loc(2)]
555  allocate(x_1d_loc%p(loc_size))
556  x_1d_loc%p = reshape(ip(x%DU_0(:,:,norm_id(1):norm_id(2),:)),&
557  &[loc_size])
558 
559  ! RE_DU_1
560  x_1d_loc => x_1d(id); id = id+1
561  x_1d_loc%var_name = 'RE_DU_1'
562  allocate(x_1d_loc%tot_i_min(4),x_1d_loc%tot_i_max(4))
563  allocate(x_1d_loc%loc_i_min(4),x_1d_loc%loc_i_max(4))
564  x_1d_loc%tot_i_min = [1,1,1,1]
565  x_1d_loc%tot_i_max = [n_tot,n_mod_x]
566  x_1d_loc%loc_i_min = [par_id(1),1,grid_trim%i_min,lim_sec_x_loc(1)]
567  x_1d_loc%loc_i_max = [par_id(2),n_tot(2),grid_trim%i_max,&
568  &lim_sec_x_loc(2)]
569  allocate(x_1d_loc%p(loc_size))
570  x_1d_loc%p = reshape(rp(x%DU_1(:,:,norm_id(1):norm_id(2),:)),&
571  &[loc_size])
572 
573  ! IM_DU_1
574  x_1d_loc => x_1d(id); id = id+1
575  x_1d_loc%var_name = 'IM_DU_1'
576  allocate(x_1d_loc%tot_i_min(4),x_1d_loc%tot_i_max(4))
577  allocate(x_1d_loc%loc_i_min(4),x_1d_loc%loc_i_max(4))
578  x_1d_loc%tot_i_min = [1,1,1,1]
579  x_1d_loc%tot_i_max = [n_tot,n_mod_x]
580  x_1d_loc%loc_i_min = [par_id(1),1,grid_trim%i_min,lim_sec_x_loc(1)]
581  x_1d_loc%loc_i_max = [par_id(2),n_tot(2),grid_trim%i_max,&
582  &lim_sec_x_loc(2)]
583  allocate(x_1d_loc%p(loc_size))
584  x_1d_loc%p = reshape(ip(x%DU_1(:,:,norm_id(1):norm_id(2),:)),&
585  &[loc_size])
586 
587  ! write
588  ierr = print_hdf5_arrs(x_1d(1:id-1),pb3d_name,trim(data_name),&
589  &rich_lvl=rich_lvl,ind_print=.not.grid_trim%divided)
590  chckerr('')
591 
592  ! clean up
593  call grid_trim%dealloc()
594  call dealloc_var_1d(x_1d)
595  nullify(x_1d_loc)
596 
597  ! user output
598  call lvl_ud(-1)
599  end function print_output_x_1
601  integer function print_output_x_2(grid_X,X,data_name,rich_lvl,par_div,&
602  &lim_sec_X,is_field_averaged) result(ierr)
603  use num_vars, only: pb3d_name, eq_jobs_lims, eq_job_nr
604  use grid_utilities, only: trim_grid
605  use hdf5_ops, only: print_hdf5_arrs
606  use hdf5_vars, only: dealloc_var_1d, var_1d_type, &
607  &max_dim_var_1d
608  use x_utilities, only: is_necessary_x, get_sec_x_range
609  use x_vars, only: set_nn_mod
610  use num_utilities, only: c
611 
612  character(*), parameter :: rout_name = 'print_output_X_2'
613 
614  ! input / output
615  type(grid_type), intent(in) :: grid_x
616  type(x_2_type), intent(in) :: x
617  character(len=*), intent(in) :: data_name
618  integer, intent(in), optional :: rich_lvl
619  logical, intent(in), optional :: par_div
620  integer, intent(in), optional :: lim_sec_x(2,2)
621  logical, intent(in), optional :: is_field_averaged
622 
623  ! local variables
624  type(grid_type) :: grid_trim ! trimmed grid
625  integer :: n_tot(3) ! total n
626  integer :: par_id(2) ! parallel interval
627  integer :: par_id_loc(2) ! local parallel interval
628  integer :: norm_id(2) ! untrimmed normal indices for trimmed grids
629  type(var_1d_type), allocatable, target :: x_1d(:) ! 1D equivalent of X variables
630  type(var_1d_type), pointer :: x_1d_loc => null() ! local element in X_1D
631  logical :: print_this(2) ! whether symmetric and asymmetric variables need to be printed
632  integer :: nn_mod_tot(2) ! total nr. of modes for symmetric and asymmetric variables
633  integer :: nn_mod_loc(2) ! local nr. of modes for symmetric and asymmetric variables
634  integer :: id ! counter
635  logical :: par_div_loc ! local par_div
636  integer :: loc_size(2) ! local size for symmetric and asymmetric variables
637  integer :: m, k ! counters
638  integer :: sxr_loc(2,2) ! local secondary X limits for symmetric and asymmetric variables
639  integer :: sxr_tot(2,2) ! total secondary X limits for symmetric and asymmetric variables
640 
641  ! initialize ierr
642  ierr = 0
643 
644  ! user output
645  call writo('Write tensorial perturbation variables to output file')
646  call lvl_ud(1)
647 
648  ! trim grid
649  ierr = trim_grid(grid_x,grid_trim,norm_id)
650  chckerr('')
651 
652  ! set local par_div
653  par_div_loc = .false.
654  if (present(par_div)) par_div_loc = par_div
655 
656  ! set total n, parallel interval and total n_mod_X
657  n_tot = grid_trim%n
658  par_id = [1,n_tot(1)]
659  par_id_loc = [1,n_tot(1)]
660  if (grid_trim%n(1).gt.0 .and. par_div_loc) then ! total grid includes all equilibrium jobs
661  n_tot(1) = maxval(eq_jobs_lims)-minval(eq_jobs_lims)+1
662  par_id = eq_jobs_lims(:,eq_job_nr)
663  par_id_loc = par_id - par_id(1)+1
664  else if (present(is_field_averaged)) then
665  if (is_field_averaged) then ! only first point
666  par_id = [1,1]
667  par_id_loc = [1,1]
668  n_tot(1) = 1
669  end if
670  end if
671  nn_mod_tot(1) = set_nn_mod(.true.)
672  nn_mod_tot(2) = set_nn_mod(.false.)
673 
674  ! set dimensions, parallel limits and local and total nn_mod_X
675 
676  ! Set up the 1D equivalents of the perturbation variables
677  allocate(x_1d(max_dim_var_1d))
678 
679  ! set up variables X_1D
680  id = 1
681 
682  ! loop over modes of dimension 2
683  do m = 1,x%n_mod(2)
684  ! get contiguous range of modes of this m
685  call get_sec_x_range(sxr_loc(:,1),sxr_tot(:,1),m,.true.,lim_sec_x)
686  call get_sec_x_range(sxr_loc(:,2),sxr_tot(:,2),m,.false.,lim_sec_x)
687  nn_mod_loc = sxr_loc(2,:)-sxr_loc(1,:)+1
688  print_this = .false.
689  do k = 1,2
690  if (sxr_loc(1,k).le.sxr_loc(2,k)) print_this(k) = .true. ! a bound is found
691  end do
692 
693  ! set local size
694  loc_size(1) = (par_id(2)-par_id(1)+1)*n_tot(2)*&
695  &(norm_id(2)-norm_id(1)+1)*nn_mod_loc(1)
696  loc_size(2) = (par_id(2)-par_id(1)+1)*n_tot(2)*&
697  &(norm_id(2)-norm_id(1)+1)*nn_mod_loc(2)
698 
699  if (print_this(1)) then
700  ! RE_PV_0
701  x_1d_loc => x_1d(id); id = id+1
702  x_1d_loc%var_name = 'RE_PV_0'
703  allocate(x_1d_loc%tot_i_min(4),x_1d_loc%tot_i_max(4))
704  allocate(x_1d_loc%loc_i_min(4),x_1d_loc%loc_i_max(4))
705  x_1d_loc%tot_i_min = [1,1,1,1]
706  x_1d_loc%tot_i_max = [n_tot,nn_mod_tot(1)]
707  x_1d_loc%loc_i_min = [par_id(1),1,grid_trim%i_min,sxr_tot(1,1)]
708  x_1d_loc%loc_i_max = [par_id(2),n_tot(2),grid_trim%i_max,&
709  &sxr_tot(2,1)]
710  allocate(x_1d_loc%p(loc_size(1)))
711  x_1d_loc%p = reshape(rp(x%PV_0(par_id_loc(1):par_id_loc(2),:,&
712  &norm_id(1):norm_id(2),sxr_loc(1,1):sxr_loc(2,1))),&
713  &[loc_size(1)])
714 
715  ! IM_PV_0
716  x_1d_loc => x_1d(id); id = id+1
717  x_1d_loc%var_name = 'IM_PV_0'
718  allocate(x_1d_loc%tot_i_min(4),x_1d_loc%tot_i_max(4))
719  allocate(x_1d_loc%loc_i_min(4),x_1d_loc%loc_i_max(4))
720  x_1d_loc%tot_i_min = [1,1,1,1]
721  x_1d_loc%tot_i_max = [n_tot,nn_mod_tot(1)]
722  x_1d_loc%loc_i_min = [par_id(1),1,grid_trim%i_min,sxr_tot(1,1)]
723  x_1d_loc%loc_i_max = [par_id(2),n_tot(2),grid_trim%i_max,&
724  &sxr_tot(2,1)]
725  allocate(x_1d_loc%p(loc_size(1)))
726  x_1d_loc%p = reshape(ip(x%PV_0(par_id_loc(1):par_id_loc(2),:,&
727  &norm_id(1):norm_id(2),sxr_loc(1,1):sxr_loc(2,1))),&
728  &[loc_size(1)])
729 
730  ! RE_PV_2
731  x_1d_loc => x_1d(id); id = id+1
732  x_1d_loc%var_name = 'RE_PV_2'
733  allocate(x_1d_loc%tot_i_min(4),x_1d_loc%tot_i_max(4))
734  allocate(x_1d_loc%loc_i_min(4),x_1d_loc%loc_i_max(4))
735  x_1d_loc%tot_i_min = [1,1,1,1]
736  x_1d_loc%tot_i_max = [n_tot,nn_mod_tot(1)]
737  x_1d_loc%loc_i_min = [par_id(1),1,grid_trim%i_min,sxr_tot(1,1)]
738  x_1d_loc%loc_i_max = [par_id(2),n_tot(2),grid_trim%i_max,&
739  &sxr_tot(2,1)]
740  allocate(x_1d_loc%p(loc_size(1)))
741  x_1d_loc%p = reshape(rp(x%PV_2(par_id_loc(1):par_id_loc(2),:,&
742  &norm_id(1):norm_id(2),sxr_loc(1,1):sxr_loc(2,1))),&
743  &[loc_size(1)])
744 
745  ! IM_PV_2
746  x_1d_loc => x_1d(id); id = id+1
747  x_1d_loc%var_name = 'IM_PV_2'
748  allocate(x_1d_loc%tot_i_min(4),x_1d_loc%tot_i_max(4))
749  allocate(x_1d_loc%loc_i_min(4),x_1d_loc%loc_i_max(4))
750  x_1d_loc%tot_i_min = [1,1,1,1]
751  x_1d_loc%tot_i_max = [n_tot,nn_mod_tot(1)]
752  x_1d_loc%loc_i_min = [par_id(1),1,grid_trim%i_min,sxr_tot(1,1)]
753  x_1d_loc%loc_i_max = [par_id(2),n_tot(2),grid_trim%i_max,&
754  &sxr_tot(2,1)]
755  allocate(x_1d_loc%p(loc_size(1)))
756  x_1d_loc%p = reshape(ip(x%PV_2(par_id_loc(1):par_id_loc(2),:,&
757  &norm_id(1):norm_id(2),sxr_loc(1,1):sxr_loc(2,1))),&
758  &[loc_size(1)])
759 
760  ! RE_KV_0
761  x_1d_loc => x_1d(id); id = id+1
762  x_1d_loc%var_name = 'RE_KV_0'
763  allocate(x_1d_loc%tot_i_min(4),x_1d_loc%tot_i_max(4))
764  allocate(x_1d_loc%loc_i_min(4),x_1d_loc%loc_i_max(4))
765  x_1d_loc%tot_i_min = [1,1,1,1]
766  x_1d_loc%tot_i_max = [n_tot,nn_mod_tot(1)]
767  x_1d_loc%loc_i_min = [par_id(1),1,grid_trim%i_min,sxr_tot(1,1)]
768  x_1d_loc%loc_i_max = [par_id(2),n_tot(2),grid_trim%i_max,&
769  &sxr_tot(2,1)]
770  allocate(x_1d_loc%p(loc_size(1)))
771  x_1d_loc%p = reshape(rp(x%KV_0(par_id_loc(1):par_id_loc(2),:,&
772  &norm_id(1):norm_id(2),sxr_loc(1,1):sxr_loc(2,1))),&
773  &[loc_size(1)])
774 
775  ! IM_KV_0
776  x_1d_loc => x_1d(id); id = id+1
777  x_1d_loc%var_name = 'IM_KV_0'
778  allocate(x_1d_loc%tot_i_min(4),x_1d_loc%tot_i_max(4))
779  allocate(x_1d_loc%loc_i_min(4),x_1d_loc%loc_i_max(4))
780  x_1d_loc%tot_i_min = [1,1,1,1]
781  x_1d_loc%tot_i_max = [n_tot,nn_mod_tot(1)]
782  x_1d_loc%loc_i_min = [par_id(1),1,grid_trim%i_min,sxr_tot(1,1)]
783  x_1d_loc%loc_i_max = [par_id(2),n_tot(2),grid_trim%i_max,&
784  &sxr_tot(2,1)]
785  allocate(x_1d_loc%p(loc_size(1)))
786  x_1d_loc%p = reshape(ip(x%KV_0(par_id_loc(1):par_id_loc(2),:,&
787  &norm_id(1):norm_id(2),sxr_loc(1,1):sxr_loc(2,1))),&
788  &[loc_size(1)])
789 
790  ! RE_KV_2
791  x_1d_loc => x_1d(id); id = id+1
792  x_1d_loc%var_name = 'RE_KV_2'
793  allocate(x_1d_loc%tot_i_min(4),x_1d_loc%tot_i_max(4))
794  allocate(x_1d_loc%loc_i_min(4),x_1d_loc%loc_i_max(4))
795  x_1d_loc%tot_i_min = [1,1,1,1]
796  x_1d_loc%tot_i_max = [n_tot,nn_mod_tot(1)]
797  x_1d_loc%loc_i_min = [par_id(1),1,grid_trim%i_min,sxr_tot(1,1)]
798  x_1d_loc%loc_i_max = [par_id(2),n_tot(2),grid_trim%i_max,&
799  &sxr_tot(2,1)]
800  allocate(x_1d_loc%p(loc_size(1)))
801  x_1d_loc%p = reshape(rp(x%KV_2(par_id_loc(1):par_id_loc(2),:,&
802  &norm_id(1):norm_id(2),sxr_loc(1,1):sxr_loc(2,1))),&
803  &[loc_size(1)])
804 
805  ! IM_KV_2
806  x_1d_loc => x_1d(id); id = id+1
807  x_1d_loc%var_name = 'IM_KV_2'
808  allocate(x_1d_loc%tot_i_min(4),x_1d_loc%tot_i_max(4))
809  allocate(x_1d_loc%loc_i_min(4),x_1d_loc%loc_i_max(4))
810  x_1d_loc%tot_i_min = [1,1,1,1]
811  x_1d_loc%tot_i_max = [n_tot,nn_mod_tot(1)]
812  x_1d_loc%loc_i_min = [par_id(1),1,grid_trim%i_min,sxr_tot(1,1)]
813  x_1d_loc%loc_i_max = [par_id(2),n_tot(2),grid_trim%i_max,&
814  &sxr_tot(2,1)]
815  allocate(x_1d_loc%p(loc_size(1)))
816  x_1d_loc%p = reshape(ip(x%KV_2(par_id_loc(1):par_id_loc(2),:,&
817  &norm_id(1):norm_id(2),sxr_loc(1,1):sxr_loc(2,1))),&
818  &[loc_size(1)])
819  end if
820 
821  if (print_this(2)) then
822  ! RE_PV_1
823  x_1d_loc => x_1d(id); id = id+1
824  x_1d_loc%var_name = 'RE_PV_1'
825  allocate(x_1d_loc%tot_i_min(4),x_1d_loc%tot_i_max(4))
826  allocate(x_1d_loc%loc_i_min(4),x_1d_loc%loc_i_max(4))
827  x_1d_loc%tot_i_min = [1,1,1,1]
828  x_1d_loc%tot_i_max = [n_tot,nn_mod_tot(2)]
829  x_1d_loc%loc_i_min = [par_id(1),1,grid_trim%i_min,sxr_tot(1,2)]
830  x_1d_loc%loc_i_max = [par_id(2),n_tot(2),grid_trim%i_max,&
831  &sxr_tot(2,2)]
832  allocate(x_1d_loc%p(loc_size(2)))
833  x_1d_loc%p = reshape(rp(x%PV_1(par_id_loc(1):par_id_loc(2),:,&
834  &norm_id(1):norm_id(2),sxr_loc(1,2):sxr_loc(2,2))),&
835  &[loc_size(2)])
836 
837  ! IM_PV_1
838  x_1d_loc => x_1d(id); id = id+1
839  x_1d_loc%var_name = 'IM_PV_1'
840  allocate(x_1d_loc%tot_i_min(4),x_1d_loc%tot_i_max(4))
841  allocate(x_1d_loc%loc_i_min(4),x_1d_loc%loc_i_max(4))
842  x_1d_loc%tot_i_min = [1,1,1,1]
843  x_1d_loc%tot_i_max = [n_tot,nn_mod_tot(2)]
844  x_1d_loc%loc_i_min = [par_id(1),1,grid_trim%i_min,sxr_tot(1,2)]
845  x_1d_loc%loc_i_max = [par_id(2),n_tot(2),grid_trim%i_max,&
846  &sxr_tot(2,2)]
847  allocate(x_1d_loc%p(loc_size(2)))
848  x_1d_loc%p = reshape(ip(x%PV_1(par_id_loc(1):par_id_loc(2),:,&
849  &norm_id(1):norm_id(2),sxr_loc(1,2):sxr_loc(2,2))),&
850  &[loc_size(2)])
851 
852  ! RE_KV_1
853  x_1d_loc => x_1d(id); id = id+1
854  x_1d_loc%var_name = 'RE_KV_1'
855  allocate(x_1d_loc%tot_i_min(4),x_1d_loc%tot_i_max(4))
856  allocate(x_1d_loc%loc_i_min(4),x_1d_loc%loc_i_max(4))
857  x_1d_loc%tot_i_min = [1,1,1,1]
858  x_1d_loc%tot_i_max = [n_tot,nn_mod_tot(2)]
859  x_1d_loc%loc_i_min = [par_id(1),1,grid_trim%i_min,sxr_tot(1,2)]
860  x_1d_loc%loc_i_max = [par_id(2),n_tot(2),grid_trim%i_max,&
861  &sxr_tot(2,2)]
862  allocate(x_1d_loc%p(loc_size(2)))
863  x_1d_loc%p = reshape(rp(x%KV_1(par_id_loc(1):par_id_loc(2),:,&
864  &norm_id(1):norm_id(2),sxr_loc(1,2):sxr_loc(2,2))),&
865  &[loc_size(2)])
866 
867  ! IM_KV_1
868  x_1d_loc => x_1d(id); id = id+1
869  x_1d_loc%var_name = 'IM_KV_1'
870  allocate(x_1d_loc%tot_i_min(4),x_1d_loc%tot_i_max(4))
871  allocate(x_1d_loc%loc_i_min(4),x_1d_loc%loc_i_max(4))
872  x_1d_loc%tot_i_min = [1,1,1,1]
873  x_1d_loc%tot_i_max = [n_tot,nn_mod_tot(2)]
874  x_1d_loc%loc_i_min = [par_id(1),1,grid_trim%i_min,sxr_tot(1,2)]
875  x_1d_loc%loc_i_max = [par_id(2),n_tot(2),grid_trim%i_max,&
876  &sxr_tot(2,2)]
877  allocate(x_1d_loc%p(loc_size(2)))
878  x_1d_loc%p = reshape(ip(x%KV_1(par_id_loc(1):par_id_loc(2),:,&
879  &norm_id(1):norm_id(2),sxr_loc(1,2):sxr_loc(2,2))),&
880  &[loc_size(2)])
881  end if
882  end do
883 
884  ! write
885  ierr = print_hdf5_arrs(x_1d(1:id-1),pb3d_name,trim(data_name),&
886  &rich_lvl=rich_lvl,ind_print=.not.grid_trim%divided)
887  chckerr('')
888 
889  ! clean up
890  call grid_trim%dealloc()
891  call dealloc_var_1d(x_1d)
892  nullify(x_1d_loc)
893 
894  ! user output
895  call lvl_ud(-1)
896  end function print_output_x_2
897 
918  integer function init_modes(grid_eq,eq) result(ierr)
919  use num_vars, only: use_pol_flux_f, x_style
920  use x_vars, only: prim_x, min_sec_x, max_sec_x, n_mod_x, min_n_x, &
921  &max_n_x, min_m_x, max_m_x, min_nm_x
922  use mpi_utilities, only: get_ser_var
923  use grid_utilities, only: trim_grid
924 
925  character(*), parameter :: rout_name = 'init_modes'
926 
927  ! input / output
928  type(grid_type), intent(in) :: grid_eq
929  type(eq_1_type), intent(in) :: eq
930 
931  ! local variables
932  type(grid_type) :: grid_eq_trim ! trimmed version of equilibrium
933  real(dp), allocatable :: jq_tot(:) ! saf. fac. or rot. transf. and derivs. in Flux coords.
934  integer :: norm_eq_id(2) ! untrimmed normal indices for trimmed grids
935 
936  ! initialize ierr
937  ierr = 0
938 
939  ! user output
940  call writo('Initializing perturbation modes variables')
941  call lvl_ud(1)
942 
943  ! get trimmed grids
944  ierr = trim_grid(grid_eq,grid_eq_trim,norm_eq_id)
945  chckerr('')
946 
947  ! (re)allocate variables
948  if (allocated(min_n_x)) deallocate(min_n_x)
949  if (allocated(max_n_x)) deallocate(max_n_x)
950  if (allocated(min_m_x)) deallocate(min_m_x)
951  if (allocated(max_m_x)) deallocate(max_m_x)
952  allocate(min_n_x(grid_eq_trim%n(3)),max_n_x(grid_eq_trim%n(3)))
953  allocate(min_m_x(grid_eq_trim%n(3)),max_m_x(grid_eq_trim%n(3)))
954 
955  ! get serial version of safety factor or rot. transform
956  allocate(jq_tot(grid_eq_trim%n(3)))
957  if (use_pol_flux_f) then
958  if (grid_eq%divided) then
959  ierr = get_ser_var(eq%q_saf_FD(norm_eq_id(1):norm_eq_id(2),0),&
960  &jq_tot,scatter=.true.)
961  chckerr('')
962  else
963  jq_tot = eq%q_saf_FD(:,0)
964  end if
965  else
966  if (grid_eq%divided) then
967  ierr = get_ser_var(eq%rot_t_FD(norm_eq_id(1):norm_eq_id(2),0),&
968  &jq_tot,scatter=.true.)
969  chckerr('')
970  else
971  jq_tot = eq%rot_t_FD(:,0)
972  end if
973  end if
974 
975  ! set the limits depending on the X style
976  select case (x_style)
977  case (1) ! prescribed
978  if (use_pol_flux_f) then
979  min_n_x = prim_x
980  max_n_x = prim_x
981  min_m_x = min_sec_x
982  max_m_x = max_sec_x
983  else
984  min_n_x = min_sec_x
985  max_n_x = max_sec_x
986  min_m_x = prim_x
987  max_m_x = prim_x
988  end if
989  case (2) ! fast
990  if (use_pol_flux_f) then
991  min_n_x = prim_x
992  max_n_x = prim_x
993  if (prim_x*jq_tot(1).gt.0) then ! nq > 0: m > 0
994  min_m_x = nint(prim_x*jq_tot-(n_mod_x-1)*0.5)
995  min_m_x = max(min_nm_x,min_m_x)
996  max_m_x = min_m_x + n_mod_x - 1
997  else ! nq < 0: m < 0
998  max_m_x = nint(prim_x*jq_tot+(n_mod_x-1)*0.5)
999  max_m_x = min(-min_nm_x,max_m_x)
1000  min_m_x = max_m_x - n_mod_x + 1
1001  end if
1002  else
1003  if (prim_x*jq_tot(1).gt.0) then ! m iota > 0: n > 0
1004  min_n_x = nint(prim_x*jq_tot-(n_mod_x-1)*0.5)
1005  min_n_x = max(min_nm_x,min_n_x)
1006  max_n_x = min_n_x + n_mod_x - 1
1007  else ! m iota < 0: n < 0
1008  max_n_x = nint(prim_x*jq_tot+(n_mod_x-1)*0.5)
1009  max_n_x = min(-min_nm_x,max_n_x)
1010  min_n_x = max_n_x - n_mod_x + 1
1011  end if
1012  min_m_x = prim_x
1013  max_m_x = prim_x
1014  end if
1015  end select
1016 
1017  ! clean up
1018  call grid_eq_trim%dealloc()
1019 
1020  call lvl_ud(-1)
1021  call writo('Perturbation modes variables initialized')
1022  end function init_modes
1023 
1059  integer function setup_modes(mds,grid_eq,grid,plot_name) result(ierr)
1060  use num_vars, only: use_pol_flux_f, rank
1061  use x_vars, only: n_mod_x, min_n_x, min_m_x, max_n_x, max_m_x
1062  use grid_utilities, only: trim_grid
1063  use eq_vars, only: max_flux_f
1064  use num_utilities, only: spline
1065 
1066  character(*), parameter :: rout_name = 'setup_modes'
1067 
1068  ! input / output
1069  type(modes_type), intent(inout), target :: mds
1070  type(grid_type), intent(in) :: grid_eq
1071  type(grid_type), intent(in) :: grid
1072  character(len=*), intent(in), optional :: plot_name
1073 
1074  ! local variables
1075  real(dp), allocatable :: min_nm_x(:) ! interpolated minimum mode number
1076  real(dp), allocatable :: x_plot(:,:) ! abscissa of plot
1077  integer, pointer :: nm_x(:,:) ! either n or m
1078  integer :: id, ld, kd ! counters
1079  integer :: ld_loc ! shfited ld
1080  integer :: n_mod_tot ! total number of modes
1081  integer :: ind_id ! current size of ind_tot
1082  integer :: ld_shift ! shift in table indices
1083  integer :: mod_x_range ! total mode range
1084  integer :: delta_ld ! change in total mode numbers
1085  integer, allocatable :: ind_cur(:) ! current indices
1086  integer, allocatable :: ind_tot(:,:) ! total index information, will be later cut to sec
1087  character(len=max_str_ln), allocatable :: plot_titles(:) ! title for plots
1088  character(len=max_str_ln) :: plot_name_tot ! file name for plots
1089  character(len=max_str_ln) :: err_msg ! error message
1090 #if ldebug
1091  real(dp), allocatable :: y_plot(:,:) ! ordinate of plot
1092 #endif
1093 
1094  ! initialize ierr
1095  ierr = 0
1096 
1097  ! user output
1098  call writo('Setting up general modes variables')
1099  call lvl_ud(1)
1100 
1101  ! set up normal tabulation values
1102  allocate(ind_tot(-grid%n(3)*n_mod_x:grid%n(3)*n_mod_x,4)) ! not quite absolute maximum but should never be reached
1103  allocate(ind_cur(n_mod_x)) ! indices of total modes currently being treated
1104 
1105  ! calculate n and m
1106  allocate(mds%n(grid%n(3),n_mod_x))
1107  allocate(mds%m(grid%n(3),n_mod_x))
1108 
1109  ! auxiliary variable
1110  allocate(min_nm_x(grid%n(3)))
1111 
1112  ! iterate over n (id=1) and m (id=2)
1113  do id = 1,2
1114  ! point nm_X, and interpolate mode minimum (linear, to preserve
1115  ! form)
1116  if (id.eq.1) then
1117  nm_x => mds%n
1118  mod_x_range = max_n_x(1)-min_n_x(1) ! assumes that this is true for all normal points
1119  ierr = spline(grid_eq%r_F,min_n_x*1._dp,grid%r_F,min_nm_x,ord=1)
1120  chckerr('')
1121  else
1122  nm_x => mds%m
1123  mod_x_range = max_m_x(1)-min_m_x(1) ! assumes that this is true for all normal points
1124  ierr = spline(grid_eq%r_F,min_m_x*1._dp,grid%r_F,min_nm_x,ord=1)
1125  chckerr('')
1126  end if
1127 
1128  !!! For non-monotonous tests:
1129  !!!if (id .eq. 2) then
1130  !!!min_nm_X(nint(grid%n(3)*0.7):grid%n(3)) = &
1131  !!!&2*min_nm_X(nint(grid%n(3)*0.7)) - &
1132  !!!&min_nm_X(nint(grid%n(3)*0.7):grid%n(3))
1133  !!!end if
1134 
1135  ! initialize variables for first normal grid point
1136  ld_shift = 0
1137  ind_id = 0
1138  do ld = 1,n_mod_x
1139  nm_x(1,ld) = nint(min_nm_x(1)) + mod_x_range*(ld-1)/(n_mod_x-1)
1140  ind_tot(ld,:) = [nm_x(1,ld),1,0,ld]
1141  ind_id = ind_id + 1
1142 #if ldebug
1143  if (debug_setup_modes) write(*,*) 'starting with', ld, ':', &
1144  &ind_tot(ld,:)
1145 #endif
1146  end do
1147  ind_cur = [(ld, ld=1,n_mod_x)]
1148 
1149  ! iterate over all next normal grid points
1150  do kd = 2,grid%n(3)
1151  ! update ld delta and shift
1152  delta_ld = nint(min_nm_x(kd)) - nint(min_nm_x(kd-1))
1153  ld_shift = ld_shift + delta_ld
1154 
1155  ! set n or m
1156  do ld = 1,n_mod_x
1157  ! shift ld and wrap back to [1..n_mod_X]
1158  ld_loc = modulo(ld+ld_shift-1,n_mod_x)+1
1159  nm_x(kd,ld_loc) = nint(min_nm_x(kd)) + &
1160  &mod_x_range*(ld-1)/(n_mod_x-1)
1161  end do
1162 
1163  ! if there was a shift, set new total index information
1164  if (abs(delta_ld).gt.0) then
1165 #if ldebug
1166  if (debug_setup_modes) write(*,*) 'kd', kd, 'delta', &
1167  &delta_ld
1168 #endif
1169  do ld = 1,abs(delta_ld)
1170  if (delta_ld.gt.0) then ! mode number has increased
1171  ind_tot(ind_cur(1),3) = kd - 1 ! upper limit in normal range
1172  ind_tot(ind_id+ld,:) = [ &
1173  &ind_tot(ind_cur(n_mod_x),1) + 1, & ! total mode number
1174  &kd, & ! lower limit in normal range
1175  &0, & ! initalize upper limit normal range
1176  &modulo(ind_tot(ind_id,4)+ld-1,n_mod_x) + 1] ! index in tables, shifted and wrapped around
1177 #if ldebug
1178  if (debug_setup_modes) then
1179  write(*,*) ' FINISHES', ind_cur(1), &
1180  &':', ind_tot(ind_cur(1),:)
1181  write(*,*) ' STARTS', ind_id+ld, ':', &
1182  &ind_tot(ind_id+ld,:)
1183  end if
1184 #endif
1185  ind_cur = [ind_cur(2:n_mod_x),ind_id+ld] ! add ind_id+ld at top
1186  else ! mode number has decreased
1187  ind_tot(ind_cur(n_mod_x),3) = kd - 1 ! upper limit normal range
1188  ind_tot(ind_id+ld,:) = [ &
1189  &ind_tot(ind_cur(1),1) - 1, & ! total mode number
1190  &kd, & ! lower limit in normal range
1191  &0, & ! initalize upper limit normal range
1192  &modulo(ind_tot(ind_id,4)-ld,n_mod_x) + 1] ! index in tables, shifted and wrapped around
1193 #if ldebug
1194  if (debug_setup_modes) then
1195  write(*,*) ' FINISHES', ind_cur(n_mod_x), &
1196  &':', ind_tot(ind_cur(n_mod_x),:)
1197  write(*,*) ' STARTS', ind_id+ld, ':', &
1198  &ind_tot(ind_id+ld,:)
1199  end if
1200 #endif
1201  ind_cur = [ind_id+ld,ind_cur(1:n_mod_x-1)] ! add ind_id+ld at bottom
1202  end if
1203  end do
1204  ind_id = ind_id+abs(delta_ld)
1205 #if ldebug
1206  if (debug_setup_modes) write(*,*) ' current indices:', &
1207  &ind_cur
1208 #endif
1209  end if
1210  end do
1211 
1212  ! close last total modes
1213  do ld = 1,size(ind_cur)
1214  ind_tot(ind_cur(ld),3) = grid%n(3)
1215  end do
1216 
1217  ! save in index information
1218  if ((use_pol_flux_f .and. id.eq.2) .or. &
1219  &(.not.use_pol_flux_f .and. id.eq.1)) then
1220  allocate(mds%sec(ind_id,4))
1221  mds%sec = ind_tot(1:ind_id,:)
1222  end if
1223  end do
1224 
1225  ! master plots output if requested
1226  if (rank.eq.0 .and. present(plot_name)) then
1227  call writo('Plotting mode numbers')
1228  call lvl_ud(1)
1229 
1230  ! total number of modes
1231  n_mod_tot = size(mds%sec,1)
1232 
1233  ! set up x values
1234  allocate(x_plot(grid%n(3),n_mod_x))
1235  do ld = 1,n_mod_x
1236  x_plot(:,ld) = grid%r_F
1237  end do
1238  x_plot = x_plot*2*pi/max_flux_f
1239  allocate(plot_titles(1))
1240 
1241  ! plot poloidal modes
1242  plot_titles(1) = 'poloidal mode numbers'
1243  plot_name_tot = 'modes_m_'//trim(plot_name)
1244  call print_ex_2d(plot_titles(1:1),plot_name_tot,mds%m*1._dp,&
1245  &x=x_plot,draw=.false.)
1246  call draw_ex(plot_titles(1:1),plot_name_tot,n_mod_x,1,.false.)
1247 
1248  ! plot toroidal modes
1249  plot_titles(1) = 'toroidal mode numbers'
1250  plot_name_tot = 'modes_n_'//trim(plot_name)
1251  call print_ex_2d(plot_titles(1:1),plot_name_tot,mds%n*1._dp,&
1252  &x=x_plot,draw=.false.)
1253  call draw_ex(plot_titles(1:1),plot_name_tot,n_mod_x,1,.false.)
1254 
1255  ! clean up
1256  deallocate(x_plot)
1257  deallocate(plot_titles)
1258 
1259 #if ldebug
1260  ! plot secondary limits
1261  allocate(x_plot(n_mod_tot,n_mod_tot+2))
1262  allocate(y_plot(n_mod_tot,n_mod_tot+2))
1263  allocate(plot_titles(n_mod_tot+2))
1264  do ld = 1,n_mod_tot
1265  x_plot(:,ld) = ld
1266  y_plot(:,ld) = mds%sec(ld,2) ! lower boundary
1267  y_plot(n_mod_tot,ld) = mds%sec(ld,3) ! upper boundary
1268  plot_titles(ld) = 'mode '//trim(i2str(mds%sec(ld,1)))
1269  end do
1270  x_plot(:,n_mod_tot+1) = x_plot(1,1:n_mod_tot)
1271  y_plot(:,n_mod_tot+1) = mds%sec(:,1) ! mode number
1272  plot_titles(n_mod_tot+1) = 'mode number'
1273  x_plot(:,n_mod_tot+2) = x_plot(1,1:n_mod_tot)
1274  y_plot(:,n_mod_tot+2) = mds%sec(:,4) ! index
1275  plot_titles(n_mod_tot+2) = 'index'
1276  plot_name_tot = 'modes_sec_'//trim(plot_name)
1277  call print_ex_2d(plot_titles,plot_name_tot,y_plot,x=x_plot,&
1278  &draw=.false.)
1279  call draw_ex(plot_titles,plot_name_tot,n_mod_tot+2,1,.false.)
1280  plot_name_tot = 'modes_sec_flux_'//trim(plot_name)
1281  do ld = 1,n_mod_tot
1282  y_plot(:,ld) = grid%r_F(mds%sec(ld,2))*2*pi/max_flux_f ! lower boundary
1283  y_plot(n_mod_tot,ld) = grid%r_F(mds%sec(ld,3))*2*pi/max_flux_f ! upper boundary
1284  end do
1285  call print_ex_2d(plot_titles(1:n_mod_tot),plot_name_tot,&
1286  &y_plot(:,1:n_mod_tot),x=x_plot(:,1:n_mod_tot),draw=.false.)
1287  call draw_ex(plot_titles(1:n_mod_tot),plot_name_tot,n_mod_tot,1,&
1288  &.false.)
1289  deallocate(x_plot)
1290  deallocate(y_plot)
1291  deallocate(plot_titles)
1292 #endif
1293 
1294  call lvl_ud(-1)
1295  end if
1296 
1297  ! test whether all modes have at least one normal position
1298  if (minval(mds%sec(:,3)-mds%sec(:,2)+1).lt.1) then
1299  ierr = 1
1300  call writo('Mode '//&
1301  &trim(i2str(minloc(mds%sec(:,3)-mds%sec(:,2),1)))//&
1302  &' is not present in any flux surface')
1303  err_msg = 'Aument number of modes or choose finer equilibrium'
1304  chckerr(err_msg)
1305  end if
1306 
1307  ! clean up
1308  nullify(nm_x)
1309 
1310  call lvl_ud(-1)
1311  call writo('General mode variables set up')
1312  end function setup_modes
1313 
1349  integer function check_x_modes(grid_eq,eq) result(ierr)
1350  use num_vars, only: x_style, rank, use_pol_flux_f, n_procs
1351  use mpi_utilities, only: get_ser_var
1352 
1353  character(*), parameter :: rout_name = 'check_X_modes'
1354 
1355  ! input / output
1356  type(grid_type), intent(in) :: grid_eq
1357  type(eq_1_type), intent(in) :: eq
1358 
1359  ! local variables
1360  character(len=max_str_ln) :: err_msg ! error message
1361 
1362  ! initialize ierr
1363  ierr = 0
1364 
1365  call writo('Checking mode numbers')
1366  call lvl_ud(1)
1367 
1368  ! check the modes depending on the X style
1369  select case (x_style)
1370  case (1) ! prescribed
1371  ierr = check_x_modes_1(eq)
1372  chckerr('')
1373  case (2) ! fast
1374  ierr = check_x_modes_2(grid_eq,eq)
1375  chckerr('')
1376  end select
1377 
1378  call lvl_ud(-1)
1379  call writo('Mode numbers checked')
1380  contains
1381  ! Version for X style 1: Checks whether there exists a range in which
1382  ! each of the modes resonates, with a certain tolerance tol_norm.
1384  integer function check_x_modes_1(eq) result(ierr)
1385  use x_vars, only: prim_x, min_sec_x, n_mod_x
1386  use num_vars, only: tol_norm
1387 
1388  character(*), parameter :: rout_name = 'check_X_modes_1'
1389 
1390  ! input / output
1391  type(eq_1_type), intent(in) :: eq ! flux equilibrium
1392 
1393  ! local variables
1394  integer :: id ! counter
1395  integer :: pmone ! plus or minus one
1396  real(dp) :: min_jq, max_jq ! min. and max. values of q or iota
1397  real(dp) :: lim_lo, lim_hi ! lower and upper limit on n/m (or m/n)
1398  real(dp), allocatable :: temp_var(:) ! temporary variable
1399  character(len=max_str_ln) :: jq_name ! either safety factor or rotational transform
1400  character(len=1) :: mode_name ! either n or m
1401  logical :: all_modes_fine ! all modes are fine
1402 
1403  ! initialize ierr
1404  ierr = 0
1405 
1406  ! user output
1407  call writo('The tolerance used is '//trim(r2strt(tol_norm))//'...')
1408 
1409  ! set min_jq and max_jq in Flux coordinate system
1410  if (use_pol_flux_f) then
1411  min_jq = minval(eq%q_saf_FD(:,0))
1412  max_jq = maxval(eq%q_saf_FD(:,0))
1413  else
1414  min_jq = minval(eq%rot_t_FD(:,0))
1415  max_jq = maxval(eq%rot_t_FD(:,0))
1416  end if
1417 
1418  ! combine all processes
1419  ierr = get_ser_var([min_jq],temp_var)
1420  chckerr('')
1421  if (rank.eq.0) min_jq = minval(temp_var)
1422  ierr = get_ser_var([max_jq],temp_var)
1423  chckerr('')
1424  if (rank.eq.0) max_jq = maxval(temp_var)
1425 
1426  ! output if master
1427  if (rank.eq.0) then
1428  ! set up jq name and all_modes_fine
1429  if (use_pol_flux_f) then
1430  jq_name = 'safety factor'
1431  mode_name = 'm'
1432  else
1433  jq_name = 'rotational transform'
1434  mode_name = 'n'
1435  end if
1436  all_modes_fine = .true.
1437 
1438  ! set up plus minus one, according to the sign of jq
1439  if (min_jq.lt.0 .and. max_jq.lt.0) then
1440  pmone = -1
1441  else if (min_jq.gt.0 .and. max_jq.gt.0) then
1442  pmone = 1
1443  else
1444  err_msg = trim(jq_name)//' cannot change sign'
1445  ierr = 1
1446  chckerr(err_msg)
1447  end if
1448 
1449  ! calculate upper and lower limits
1450  lim_lo = max(min_jq-tol_norm,min_jq/(1+pmone*tol_norm))
1451  lim_hi = min(max_jq+tol_norm,max_jq/(1-pmone*tol_norm))
1452 
1453  ! for every mode (n,m) check whether m/n is inside the range of
1454  ! q values or n/m inside the range of iota values
1455  do id = 1,n_mod_x
1456  ! check if limits are met
1457  if (use_pol_flux_f) then
1458  if ((min_sec_x+id-1._dp)/prim_x.lt.lim_lo .or. &
1459  &(min_sec_x+id-1._dp)/prim_x.gt.lim_hi) then
1460  call writo('for (n,m) = ('//trim(i2str(prim_x))//&
1461  &','//trim(i2str(min_sec_x+id-1))//&
1462  &'), there is no range in the plasma where the &
1463  &ratio |n q - m| << |n|,|m| is met',alert=.true.)
1464  all_modes_fine = .false.
1465  end if
1466  else
1467  if ((min_sec_x+id-1._dp)/prim_x.lt.lim_lo .or. &
1468  &(min_sec_x+id-1._dp)/prim_x.gt.lim_hi) then
1469  call writo('for (n,m) = ('//&
1470  &trim(i2str(min_sec_x+id-1))//','//&
1471  &trim(i2str(prim_x))//'), there is no range in &
1472  &the plasma where the ratio |n - iota m| << &
1473  &|m|,|n| is met',alert=.true.)
1474  all_modes_fine = .false.
1475  end if
1476  end if
1477  end do
1478 
1479  ! output message
1480  if (all_modes_fine) then
1481  call writo('The modes are all within the allowed range &
1482  &of '//trim(i2str(ceiling(prim_x*lim_lo)))//' < '//&
1483  &mode_name//' < '//trim(i2str(floor(prim_x*lim_hi)))//&
1484  &'...')
1485  else
1486  call writo('Not all modes are within the allowed range &
1487  &of '//trim(i2str(ceiling(prim_x*lim_lo)))//' < '//&
1488  &mode_name//' < '//trim(i2str(floor(prim_x*lim_hi)))//&
1489  &'...')
1490  end if
1491  end if
1492  end function check_x_modes_1
1493 
1494  ! version for X style 2: Check how efficient the chosen number of modes
1495  ! is.
1497  integer function check_x_modes_2(grid_eq,eq) result(ierr)
1498  use x_vars, only: min_n_x, min_m_x, n_mod_x
1499  use grid_utilities, only: trim_grid
1500 
1501  character(*), parameter :: rout_name = 'check_X_modes_2'
1502 
1503  ! input / output
1504  type(grid_type), intent(in) :: grid_eq ! equilibrium grid
1505  type(eq_1_type), intent(in) :: eq ! flux equilibrium
1506 
1507  ! local variables
1508  type(grid_type) :: grid_eq_trim ! trimmed equilibrium grid
1509  integer :: norm_id(2) ! untrimmed normal indices for trimmed grid
1510  integer :: jd, kd, ld ! counters
1511  real(dp), allocatable :: max_frac(:,:) ! maximum fraction
1512  real(dp), allocatable :: max_frac_sum(:) ! sum of maximum fraction for all processes
1513  real(dp), allocatable :: max_frac_max(:) ! maximum maximum fraction for all processes
1514  real(dp), allocatable :: max_frac_temp(:) ! auxilliary variable
1515  real(dp), allocatable :: n(:,:), m(:,:) ! n and m
1516  real(dp), allocatable :: fac_n(:), fac_m(:) ! factors to be multiplied with n m m
1517  character(len=max_str_ln) :: frac_name ! name of fraction
1518 #if ldebug
1519  real(dp), allocatable :: x_plot(:,:) ! for plotting
1520  character(len=max_str_ln) :: plot_title ! title for plots
1521  character(len=max_str_ln) :: plot_name ! file name for plots
1522 #endif
1523 
1524  ! initialize ierr
1525  ierr = 0
1526 
1527  ! set up helper variables
1528  allocate(max_frac(grid_eq%loc_n_r,n_mod_x))
1529  allocate(fac_n(grid_eq%loc_n_r))
1530  allocate(fac_m(grid_eq%loc_n_r))
1531  allocate(n(grid_eq%loc_n_r,n_mod_x))
1532  allocate(m(grid_eq%loc_n_r,n_mod_x))
1533  if (use_pol_flux_f) then
1534  do jd = 1,n_mod_x
1535  n(:,jd) = min_n_x(grid_eq%i_min:grid_eq%i_max)
1536  m(:,jd) = min_m_x(grid_eq%i_min:grid_eq%i_max)+jd-1
1537  end do
1538  fac_n = eq%q_saf_FD(:,0)
1539  fac_m = 1._dp
1540  else
1541  do jd = 1,n_mod_x
1542  n(:,jd) = min_n_x(grid_eq%i_min:grid_eq%i_max)+jd-1
1543  m(:,jd) = min_m_x(grid_eq%i_min:grid_eq%i_max)
1544  end do
1545  fac_n = 1._dp
1546  fac_m = eq%rot_t_FD(:,0)
1547  end if
1548 
1549  ! loop over all modes
1550  do ld = 1,n_mod_x
1551  do kd = 1,grid_eq%loc_n_r
1552  max_frac(kd,ld) = &
1553  &abs(n(kd,ld)*fac_n(kd)-m(kd,ld)*fac_m(kd))/&
1554  &min(abs(n(kd,ld)),abs(m(kd,ld)))
1555  end do
1556  end do
1557 
1558  ! trim grid
1559  ierr = trim_grid(grid_eq,grid_eq_trim,norm_id=norm_id)
1560  chckerr('')
1561 
1562  ! gather all fractions
1563  if (rank.eq.0) allocate(max_frac_temp(n_procs))
1564  if (rank.eq.0) allocate(max_frac_sum(n_mod_x))
1565  if (rank.eq.0) allocate(max_frac_max(n_mod_x))
1566  do ld = 1,n_mod_x
1567  ierr = get_ser_var([sum(max_frac(norm_id(1):norm_id(2),ld))],&
1568  &max_frac_temp)
1569  chckerr('')
1570  if (rank.eq.0) max_frac_sum(ld) = sum(max_frac_temp)
1571  ierr = get_ser_var(&
1572  &[maxval(max_frac(norm_id(1):norm_id(2),ld))],max_frac_temp)
1573  chckerr('')
1574  if (rank.eq.0) max_frac_max(ld) = maxval(max_frac_temp)
1575  end do
1576 
1577  ! output if master
1578  if (rank.eq.0) then
1579  if (use_pol_flux_f) then
1580  frac_name = '|nq-m|/|n| or |nq-m|/|m|'
1581  else
1582  frac_name = '|n-iota m|/|n| or |n-iota m|/|m|'
1583  end if
1584  call writo('The fraction '//trim(frac_name)//' is')
1585  call lvl_ud(1)
1586  do ld = 1,n_mod_x
1587  call writo('for mode '//trim(i2str(ld))//' maximally '//&
1588  &trim(r2strt(max_frac_max(ld)))//&
1589  &', average '//&
1590  &trim(r2strt(max_frac_sum(ld)/grid_eq%n(3))))
1591  end do
1592  if (n_mod_x.gt.1) call writo('so for all modes maximally '//&
1593  &trim(r2strt(maxval(max_frac_max)))//', average '//&
1594  &trim(r2strt(sum(max_frac_sum)/(grid_eq%n(3)*n_mod_x))))
1595  call lvl_ud(-1)
1596  end if
1597 #if ldebug
1598  if (debug_check_x_modes_2) then
1599  call writo('Plotting the fraction for all modes')
1600  allocate(x_plot(grid_eq%n(3),n_mod_x))
1601  do ld = 1,n_mod_x
1602  x_plot(:,ld) = grid_eq%r_F
1603  end do
1604  plot_name = 'TEST_max_frac'
1605  plot_title = 'maximum fraction'
1606  call print_ex_2d([plot_title],plot_name,max_frac,x=x_plot,&
1607  &draw=.false.)
1608  call draw_ex([plot_title],plot_name,n_mod_x,1,.false.)
1609  end if
1610 #endif
1611 
1612  ! clean up
1613  call grid_eq_trim%dealloc()
1614  end function check_x_modes_2
1615  end function check_x_modes
1616 
1637  integer function calc_res_surf(mds,grid_eq,eq,res_surf,info,jq) result(ierr)
1638  use x_vars, only: prim_x
1639  use num_vars, only: use_pol_flux_f, norm_disc_prec_eq
1640  use eq_vars, only: max_flux_f, max_flux_f
1641  use num_ops, only: calc_zero_zhang
1642  use grid_utilities, only: trim_grid
1643  use mpi_utilities, only: get_ser_var
1644 
1645  character(*), parameter :: rout_name = 'calc_res_surf'
1646 
1647  ! input / output
1648  type(modes_type), intent(in) :: mds
1649  type(grid_type), intent(in) :: grid_eq
1650  type(eq_1_type), intent(in) :: eq
1651  real(dp), intent(inout), allocatable :: res_surf(:,:)
1652  logical, intent(in), optional :: info
1653  real(dp), intent(inout), optional, allocatable :: jq(:)
1654 
1655  ! local variables (also used in child functions)
1656  real(dp) :: nmfrac_fun ! fraction m/n or n/m to determine resonant flux surface
1657 
1658  ! local variables (not to be used in child functions)
1659  integer :: norm_id(2) ! untrimmed normal indices for trimmed grid
1660  integer :: ld ! counter
1661  integer :: ld_loc ! local ld
1662  logical :: info_loc ! local version of info
1663  real(dp), allocatable :: res_surf_loc(:,:) ! local copy of res_surf
1664  real(dp), allocatable :: jq_tot(:) ! saf. fac. or rot. transf. and derivs. in Flux coords.
1665  real(dp), allocatable :: jq_loc(:) ! local version of jq
1666  real(dp) :: norm_factor ! normalization factor to make normal coordinate 0..1
1667  type(grid_type) :: grid_eq_trim ! trimmed version of equilibrium grid
1668  integer :: n_loc, m_loc ! local n and m
1669  character(len=max_str_ln) :: err_msg ! error message
1670 
1671  ! initialize ierr
1672  ierr = 0
1673 
1674  ! set local info
1675  info_loc = .false.
1676  if (present(info)) info_loc = info
1677 
1678  ! get trimmed grid
1679  ierr = trim_grid(grid_eq,grid_eq_trim,norm_id)
1680  chckerr('')
1681 
1682  ! get serial version of safety factor or rot. transform
1683  allocate(jq_tot(grid_eq_trim%n(3)))
1684  if (use_pol_flux_f) then
1685  if (grid_eq%divided) then
1686  ierr = get_ser_var(eq%q_saf_FD(norm_id(1):norm_id(2),0),&
1687  &jq_loc,scatter=.true.)
1688  chckerr('')
1689  jq_tot = jq_loc
1690  else
1691  jq_tot = eq%q_saf_FD(:,0)
1692  end if
1693  else
1694  if (grid_eq%divided) then
1695  ierr = get_ser_var(eq%rot_t_FD(norm_id(1):norm_id(2),0),&
1696  &jq_loc,scatter=.true.)
1697  chckerr('')
1698  jq_tot = jq_loc
1699  else
1700  jq_tot = eq%rot_t_FD(:,0)
1701  end if
1702  end if
1703  if (present(jq)) then
1704  allocate(jq(grid_eq_trim%n(3)))
1705  jq = jq_tot
1706  end if
1707 
1708  ! initialize local res_surf
1709  allocate(res_surf_loc(1:size(mds%sec,1),3))
1710 
1711  ! calculate normalization factor max_flux / 2pi
1712  norm_factor = max_flux_f/(2*pi)
1713 
1714  ! loop over all modes (and shift the index in m_loc and n_loc by 1)
1715  ld_loc = 1
1716  do ld = 1,size(mds%sec,1)
1717  ! Find place where q = m/n or iota = n/m in Flux coordinates by
1718  ! solving q-m/n = 0 or iota-n/m=0, using the functin jq_fun.
1719 
1720  ! set up local m, n and nmfrac for function
1721  if (use_pol_flux_f) then
1722  n_loc = prim_x
1723  m_loc = mds%sec(ld,1)
1724  nmfrac_fun = 1.0_dp*m_loc/n_loc
1725  else
1726  n_loc = mds%sec(ld,1)
1727  m_loc = prim_x
1728  nmfrac_fun = 1.0_dp*n_loc/m_loc
1729  end if
1730 
1731  ! calculate zero using Zhang
1732  if (use_pol_flux_f) then
1733  res_surf_loc(ld_loc,1) = m_loc
1734  else
1735  res_surf_loc(ld_loc,1) = n_loc
1736  end if
1737  res_surf_loc(ld_loc,3) = nmfrac_fun
1738  err_msg = calc_zero_zhang(res_surf_loc(ld_loc,2),jq_fun,&
1739  &[minval(grid_eq_trim%r_F),maxval(grid_eq_trim%r_F)])
1740 
1741  ! intercept error
1742  if (err_msg.ne.'') then
1743  if (info_loc) then
1744  call writo('Mode (n,m) = ('//trim(i2str(n_loc))//','//&
1745  &trim(i2str(m_loc))//') is not found')
1746  call lvl_ud(1)
1747  call writo(trim(err_msg))
1748  call lvl_ud(-1)
1749  end if
1750  else if (res_surf_loc(ld_loc,2).lt.minval(grid_eq_trim%r_F) .or. &
1751  &res_surf_loc(ld_loc,2).gt.maxval(grid_eq_trim%r_F)) then
1752  if (info_loc) call writo('Mode (n,m) = ('//&
1753  &trim(i2str(n_loc))//','//trim(i2str(m_loc))//&
1754  &') does not resonate in plasma')
1755  else
1756  if (info_loc) call writo('Mode (n,m) = ('//&
1757  &trim(i2str(n_loc))//','//trim(i2str(m_loc))//&
1758  &') resonates in plasma at normalized flux surface '//&
1759  &trim(r2str(res_surf_loc(ld_loc,2)/norm_factor)))
1760  ld_loc = ld_loc + 1 ! advance ld_loc
1761  end if
1762  end do
1763 
1764  ! set res_surf from local copy
1765  allocate(res_surf(ld_loc-1,3))
1766  res_surf = res_surf_loc(1:ld_loc-1,:)
1767 
1768  ! clean up
1769  call grid_eq_trim%dealloc()
1770  contains
1771  ! Returns q-m/n or iota-n/m in Flux coordinates, used to solve for q =
1772  ! m/n or iota = n/m.
1774  real(dp) function jq_fun(pt) result(res)
1775  use num_utilities, only: spline
1776 
1777  ! input / output
1778  real(dp), intent(in) :: pt ! normal position at which to evaluate
1779 
1780  ! local variables
1781  integer :: i_min, i_max ! index of minimum and maximum value of x
1782  real(dp) :: x_min, x_max ! minimum and maximum value of x
1783  real(dp) :: res_loc(1) ! local copy of res
1784 
1785  ! initialize res
1786  res = 0
1787 
1788  ! set up min. and max index and value
1789  x_min = minval(grid_eq_trim%r_F)
1790  x_max = maxval(grid_eq_trim%r_F)
1791  i_min = minloc(grid_eq_trim%r_F,1)
1792  i_max = maxloc(grid_eq_trim%r_F,1)
1793 
1794  ! check whether to interpolate or extrapolate
1795  if (pt.ge.x_min .and. pt.le.x_max) then ! point requested within range
1796  ! interpolate
1797  ierr = spline(grid_eq_trim%r_F,jq_tot-nmfrac_fun,[pt],res_loc,&
1798  &ord=norm_disc_prec_eq)
1799  chckerr('')
1800  ! copy
1801  res = res_loc(1)
1802  end if
1803  end function jq_fun
1804  end function calc_res_surf
1805 
1813  integer function resonance_plot(mds,grid_eq,eq) result(ierr)
1814  use num_vars, only: use_pol_flux_f, no_plots, n_theta_plot, &
1815  &n_zeta_plot, rank, eq_style, min_theta_plot, max_theta_plot, &
1816  &min_zeta_plot, max_zeta_plot
1817  use eq_vars, only: max_flux_f
1818  use grid_utilities, only: calc_xyz_grid, calc_eqd_grid, coord_f2e, &
1819  &trim_grid
1820  use x_vars, only: prim_x
1821  use vmec_utilities, only: calc_trigon_factors
1822 
1823  character(*), parameter :: rout_name = 'resonance_plot'
1824 
1825  ! input / output
1826  type(modes_type), intent(in) :: mds
1827  type(grid_type), intent(in) :: grid_eq
1828  type(eq_1_type), intent(in) :: eq
1829 
1830  ! local variables (not to be used in child functions)
1831  integer :: ld ! counter
1832  integer :: n_mod_loc ! local n_mod
1833  real(dp), allocatable :: res_surf(:,:) ! resonant surfaces
1834  real(dp), allocatable :: x_plot_loc(:,:) ! for plotting
1835  real(dp), allocatable :: y_plot_loc(:,:) ! for plotting
1836  character(len=max_str_ln) :: plot_title, file_name ! name of plot, of file
1837  real(dp), allocatable :: jq(:) ! saf. fac. or rot. transf. in Flux coords.
1838  integer :: n_r ! total number of normal points
1839  integer :: plot_dim(4) ! plot dimensions (total = local because only masters)
1840  type(grid_type) :: grid_trim ! trimmed version of grid
1841  type(grid_type) :: grid_plot ! grid for plotting
1842  real(dp), allocatable :: r_plot_e(:) ! normal E coordinates of plot (needed to calculate X, Y and Z)
1843  real(dp), allocatable :: theta_plot(:,:,:), zeta_plot(:,:,:) ! pol. and tor. angle of plot
1844  real(dp), allocatable :: x_plot(:,:,:,:), y_plot(:,:,:,:), &
1845  &Z_plot(:,:,:,:) ! X, Y and Z of plot of all surfaces
1846  real(dp), allocatable :: vars(:,:,:,:) ! variable to plot
1847  character(len=max_str_ln), allocatable :: plot_titles(:) ! name of plots
1848 
1849  ! initialize ierr
1850  ierr = 0
1851 
1852  ! bypass plots if no_plots
1853  if (no_plots) return
1854 
1855  ! get trimmed grid
1856  ierr = trim_grid(grid_eq,grid_trim)
1857  chckerr('')
1858 
1859  ! initialize variables
1860  n_r = grid_trim%n(3)
1861 
1862  ! print user output
1863  if (use_pol_flux_f) then
1864  call writo('Plotting safety factor q and resonant surfaces &
1865  &q = m/n...')
1866  plot_title = 'safety factor q'
1867  file_name = 'q_saf'
1868  else
1869  call writo('Plotting rotational transform iota and resonant &
1870  &surfaces iota = n/m...')
1871  plot_title = 'rotational transform iota'
1872  file_name = 'rot_t'
1873  end if
1874 
1875  call lvl_ud(1)
1876 
1877  call writo('Calculate resonant surfaces')
1878  call lvl_ud(1)
1879 
1880  ! find resonating surfaces
1881  ierr = calc_res_surf(mds,grid_eq,eq,res_surf,info=.true.,jq=jq)
1882  chckerr('')
1883 
1884  call lvl_ud(-1)
1885 
1886  ! only master
1887  if (rank.eq.0) then
1888  ! set local n_mod
1889  n_mod_loc = size(res_surf,1)
1890 
1891  ! set up plot titles
1892  allocate(plot_titles(n_mod_loc))
1893  if (use_pol_flux_f) then ! n is fixed and m = m/n n
1894  do ld = 1,n_mod_loc
1895  plot_titles(ld) = 'resonance for m,n = '//&
1896  &trim(i2str(nint(res_surf(ld,3)*prim_x)))//','//&
1897  &trim(i2str(prim_x))
1898  end do
1899  else ! m is fixed and n = n/m m
1900  do ld = 1,n_mod_loc
1901  plot_titles(ld) = 'resonance for m,n = '//&
1902  &trim(i2str(prim_x))//','//&
1903  &trim(i2str(nint(res_surf(ld,3)*prim_x)))
1904  end do
1905  end if
1906 
1907  ! initialize x_plot_loc and y_plot_loc
1908  allocate(x_plot_loc(n_r,n_mod_loc+1)); x_plot_loc = 0
1909  allocate(y_plot_loc(n_r,n_mod_loc+1)); y_plot_loc = 0
1910 
1911  ! set x_plot_loc and y_plot_loc for first column
1912  x_plot_loc(:,1) = grid_trim%r_F
1913  y_plot_loc(:,1) = jq(:)
1914 
1915  ! set x_plot_loc and y_plot_loc for other columns
1916  do ld = 1,n_mod_loc
1917  x_plot_loc(:,ld+1) = res_surf(ld,2)
1918  y_plot_loc(n_r,ld+1) = res_surf(ld,3)
1919  end do
1920 
1921  ! only if resonant surfaces
1922  if (size(res_surf,1).gt.0) then
1923  ! user message
1924  call writo('Plot resonant flux surfaces using HDF5')
1925 
1926  call lvl_ud(1)
1927 
1928  ! set up pol. and tor. angle for plot
1929  allocate(theta_plot(n_theta_plot,n_zeta_plot,1))
1930  allocate(zeta_plot(n_theta_plot,n_zeta_plot,1))
1931  ierr = calc_eqd_grid(theta_plot,min_theta_plot*pi,&
1932  &max_theta_plot*pi,1)
1933  chckerr('')
1934  ierr = calc_eqd_grid(zeta_plot,min_zeta_plot*pi,&
1935  &max_zeta_plot*pi,2)
1936  chckerr('')
1937 
1938  ! set up vars
1939  allocate(vars(n_theta_plot,n_zeta_plot,1,n_mod_loc))
1940  do ld = 1,n_mod_loc
1941  vars(:,:,:,ld) = y_plot_loc(n_r,ld+1)
1942  end do
1943 
1944  ! set dimensions for single flux surface
1945  plot_dim = [n_theta_plot,n_zeta_plot,1,n_mod_loc]
1946 
1947  ! calculate normal vars in Equilibrium coords.
1948  allocate(r_plot_e(n_mod_loc))
1949  ierr = coord_f2e(grid_eq,x_plot_loc(n_r,2:n_mod_loc+1),&
1950  &r_plot_e,r_f_array=grid_eq%r_F,r_e_array=grid_eq%r_E)
1951  chckerr('')
1952 
1953  ! create plot grid for single flux surface
1954  ierr = grid_plot%init(plot_dim(1:3))
1955  chckerr('')
1956  grid_plot%theta_E = theta_plot
1957  grid_plot%zeta_E = zeta_plot
1958 
1959  ! if VMEC, calculate trigonometric factors of plot grid
1960  if (eq_style.eq.1) then
1961  ierr = calc_trigon_factors(grid_plot%theta_E,&
1962  &grid_plot%zeta_E,grid_plot%trigon_factors)
1963  chckerr('')
1964  end if
1965 
1966  ! set up plot X, Y and Z
1967  allocate(x_plot(n_theta_plot,n_zeta_plot,1,n_mod_loc))
1968  allocate(y_plot(n_theta_plot,n_zeta_plot,1,n_mod_loc))
1969  allocate(z_plot(n_theta_plot,n_zeta_plot,1,n_mod_loc))
1970 
1971  ! loop over all resonant surfaces to calculate X, Y and Z values
1972  do ld = 1,n_mod_loc
1973  ! set loc_r_E of plot grid
1974  grid_plot%loc_r_E = r_plot_e(ld)
1975 
1976  ! calculate X, Y and Z
1977  ierr = calc_xyz_grid(grid_eq,grid_plot,x_plot(:,:,:,ld),&
1978  &y_plot(:,:,:,ld),z_plot(:,:,:,ld))
1979  chckerr('')
1980  end do
1981 
1982  ! plot using HDF5
1983  call plot_hdf5(plot_titles,file_name,vars,x=x_plot,y=y_plot,&
1984  &z=z_plot,col=1,descr='resonant surfaces')
1985 
1986  ! deallocate local variables
1987  deallocate(vars)
1988  deallocate(theta_plot,zeta_plot,r_plot_e)
1989  deallocate(x_plot,y_plot,z_plot)
1990  call grid_plot%dealloc()
1991 
1992  call lvl_ud(-1)
1993  end if
1994 
1995  call writo('Plot safety factor with resonant flux surfaces using &
1996  &external program')
1997 
1998  call lvl_ud(1)
1999 
2000  ! rescale x_plot_loc by max_flux_F/2*pi
2001  x_plot_loc = x_plot_loc*2*pi/max_flux_f
2002 
2003  ! plot using external program
2004  call print_ex_2d([plot_title,plot_titles],file_name,y_plot_loc,&
2005  &x=x_plot_loc,draw=.false.)
2006  call draw_ex([plot_title,plot_titles],file_name,n_mod_loc+1,1,&
2007  &.false.)
2008 
2009  call lvl_ud(-1)
2010 
2011  ! only if resonant surfaces
2012  if (size(res_surf,1).gt.0) then
2013  call writo('Plot 2D map of resonances using HDF5')
2014 
2015  call lvl_ud(1)
2016 
2017  ! initialize plot dimensions
2018  plot_dim = [1,n_mod_loc,1,0]
2019 
2020  ! set up X, Y and Z
2021  allocate(x_plot(plot_dim(1),plot_dim(2),plot_dim(3),1))
2022  allocate(y_plot(plot_dim(1),plot_dim(2),plot_dim(3),1))
2023  allocate(z_plot(plot_dim(1),plot_dim(2),plot_dim(3),1))
2024 
2025  do ld = 1,n_mod_loc
2026  x_plot(1,ld,1,1) = x_plot_loc(1,ld+1)
2027  y_plot(1,ld,1,1) = res_surf(ld,1)
2028  z_plot(1,ld,1,1) = 1._dp
2029  end do
2030 
2031  ! plot 2D (to be used with plots of harmonics in sol_ops)
2032  call plot_hdf5('resonating surfaces','res_surf',&
2033  &z_plot(:,:,:,1),x=x_plot(:,:,:,1),y=y_plot(:,:,:,1))
2034 
2035  ! deallocate local variables
2036  deallocate(x_plot,y_plot,z_plot)
2037 
2038  call lvl_ud(-1)
2039  end if
2040  end if
2041 
2042  ! clean up
2043  call grid_trim%dealloc()
2044 
2045  call lvl_ud(-1)
2046  end function resonance_plot
2047 
2096  integer function calc_u(grid_eq,grid_X,eq_1,eq_2,X) result(ierr)
2097  use num_vars, only: use_pol_flux_f, eq_style, u_style, &
2098  &norm_disc_prec_x, x_grid_style
2099  use num_utilities, only: c, spline
2100  use input_utilities, only: get_log, pause_prog
2101  use eq_vars, only: vac_perm
2102 #if ldebug
2103  use num_vars, only: ltest
2104 #endif
2105 
2106  character(*), parameter :: rout_name = 'calc_U'
2107 
2108  ! input / output
2109  type(grid_type), intent(in) :: grid_eq
2110  type(grid_type), intent(in) :: grid_x
2111  type(eq_1_type), intent(in), target :: eq_1
2112  type(eq_2_type), intent(in), target :: eq_2
2113  type(x_1_type), intent(inout) :: x
2114 
2115  ! local variables
2116  integer :: id, jd, kd, ld ! counters
2117  integer :: t_size ! 2 for VMEC and 1 for HELENA
2118 
2119  ! Jacobian
2120  real(dp), pointer :: j(:,:,:) ! Jacobian
2121  real(dp), pointer :: d3j(:,:,:) ! D_theta Jacobian
2122  ! lower metric factors
2123  real(dp), pointer :: g13(:,:,:) ! g_alpha,theta
2124  real(dp), pointer :: d3g13(:,:,:) ! D_theta g_alpha,theta
2125  real(dp), pointer :: g23(:,:,:) ! g_psi,theta
2126  real(dp), pointer :: d3g23(:,:,:) ! D_theta g_psi,theta
2127  real(dp), pointer :: g33(:,:,:) ! g_theta,theta
2128  real(dp), pointer :: d3g33(:,:,:) ! D_theta g_theta,theta
2129  ! upper metric factors
2130  real(dp), pointer :: h12(:,:,:) ! h^alpha,psi
2131  real(dp), pointer :: d3h12(:,:,:) ! D_theta h^alpha,psi
2132  real(dp), pointer :: h22(:,:,:) ! h^psi,psi
2133  real(dp), pointer :: d1h22(:,:,:) ! D_alpha h^psi,psi
2134  real(dp), pointer :: d3h22(:,:,:) ! D_theta h^psi,psi
2135  real(dp), pointer :: d13h22(:,:,:) ! D_alpha,theta h^psi,psi
2136  real(dp), pointer :: d33h22(:,:,:) ! D_theta,theta h^psi,psi
2137  real(dp), pointer :: h23(:,:,:) ! h^psi,theta
2138  real(dp), pointer :: d1h23(:,:,:) ! D_alpha h^psi,theta
2139  real(dp), pointer :: d3h23(:,:,:) ! D_theta h^psi,theta
2140  real(dp), pointer :: d13h23(:,:,:) ! D_alpha,theta h^psi,theta
2141  real(dp), pointer :: d33h23(:,:,:) ! D_theta,theta h^psi,theta
2142  ! helper variables
2143  real(dp), pointer :: ang_par_f(:,:,:) ! equilibrium parallel angle in flux coordinates
2144  real(dp), pointer :: ang_geo_f(:,:,:) ! equilibrium geodesical angle in flux coordinates
2145  real(dp), pointer :: q_saf(:), rot_t(:) ! safety factor and rotational transform in X grid
2146  real(dp), allocatable :: djq(:) ! either q' (pol. flux) or -iota' (tor. flux)
2147  real(dp), allocatable :: theta_3(:,:,:), d1theta_3(:,:,:), &
2148  &D3Theta_3(:,:,:) ! Theta^theta and derivatives
2149  real(dp), allocatable :: d13theta_3(:,:,:), d33theta_3(:,:,:) ! Theta^theta and derivatives
2150  real(dp), allocatable :: n_frac(:,:) ! nq-m (pol. flux) or n-iotam (tor. flux) for all modes
2151  real(dp), allocatable :: nm(:,:) ! either n*A_0 (pol. flux) or m (tor.flux)
2152  complex(dp), allocatable :: u_corr(:,:,:,:) ! correction to U_0 and U_1 for a certain (n,m)
2153  complex(dp), allocatable :: d3u_corr(:,:,:,:) ! D_theta U_corr
2154  ! U factors
2155  real(dp), allocatable, target :: t1(:,:,:,:) ! B_alpha/B_theta and par. deriv.
2156  real(dp), allocatable, target :: t2(:,:,:,:) ! Theta^alpha + q' theta and par. deriv.
2157  real(dp), allocatable, target :: t3(:,:,:,:) ! B_alpha/B_theta q' + J^2/B_theta mu_0 p' and par. deriv.
2158  real(dp), allocatable, target :: t4(:,:,:,:) ! B_alpha/B_theta q' theta - B_psi/B_theta and par. deriv.
2159  real(dp), allocatable, target :: t5(:,:,:,:) ! B_alpha/B_theta D3Theta^theta - D1Theta^theta and par. deriv.
2160  real(dp), allocatable, target :: t6(:,:,:,:) ! B_alpha/B_theta Theta^theta and par. deriv.
2161  real(dp), pointer :: t1_x(:,:,:,:) ! T1 in X grid and par. deriv.
2162  real(dp), pointer :: t2_x(:,:,:,:) ! T2 in X grid and par. deriv.
2163  real(dp), pointer :: t3_x(:,:,:,:) ! T3 in X grid and par. deriv.
2164  real(dp), pointer :: t4_x(:,:,:,:) ! T4 in X grid and par. deriv.
2165  real(dp), pointer :: t5_x(:,:,:,:) ! T5 in X grid and par. deriv.
2166  real(dp), pointer :: t6_x(:,:,:,:) ! T6 in X grid and par. deriv.
2167 
2168  ! initialize ierr
2169  ierr = 0
2170 
2171  ! allocate variables
2172  ! helper variables
2173  allocate(nm(grid_x%loc_n_r,x%n_mod))
2174  allocate(q_saf(grid_x%loc_n_r))
2175  allocate(rot_t(grid_x%loc_n_r))
2176  allocate(djq(grid_eq%loc_n_r))
2177  allocate(theta_3(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r))
2178  allocate(d1theta_3(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r))
2179  allocate(d3theta_3(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r))
2180  allocate(d13theta_3(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r))
2181  allocate(d33theta_3(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r))
2182  allocate(n_frac(grid_x%loc_n_r,x%n_mod))
2183  allocate(u_corr(grid_x%n(1),grid_x%n(2),grid_x%loc_n_r,2))
2184  allocate(d3u_corr(grid_x%n(1),grid_x%n(2),grid_x%loc_n_r,2))
2185  ! U factors
2186  select case (eq_style)
2187  case (1) ! VMEC
2188  t_size = 2
2189  case (2) ! HELENA
2190  t_size = 1
2191  end select
2192  allocate(t1(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r,t_size))
2193  allocate(t2(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r,t_size))
2194  allocate(t3(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r,t_size))
2195  allocate(t4(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r,t_size))
2196  allocate(t5(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r,t_size))
2197  allocate(t6(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r,t_size))
2198  select case (x_grid_style)
2199  case (1) ! equilibrium
2200  ! do nothing
2201  case (2,3) ! solution or enriched
2202  allocate(q_saf(grid_x%loc_n_r))
2203  allocate(rot_t(grid_x%loc_n_r))
2204  allocate(t1_x(grid_x%n(1),grid_x%n(2),grid_x%loc_n_r,t_size))
2205  allocate(t2_x(grid_x%n(1),grid_x%n(2),grid_x%loc_n_r,t_size))
2206  allocate(t3_x(grid_x%n(1),grid_x%n(2),grid_x%loc_n_r,t_size))
2207  allocate(t4_x(grid_x%n(1),grid_x%n(2),grid_x%loc_n_r,t_size))
2208  allocate(t5_x(grid_x%n(1),grid_x%n(2),grid_x%loc_n_r,t_size))
2209  allocate(t6_x(grid_x%n(1),grid_x%n(2),grid_x%loc_n_r,t_size))
2210  end select
2211 
2212  ! set pointers
2213  if (use_pol_flux_f) then
2214  ang_par_f => grid_eq%theta_F
2215  ang_geo_f => grid_eq%zeta_F
2216  else
2217  ang_par_f => grid_eq%zeta_F
2218  ang_geo_f => grid_eq%theta_F
2219  end if
2220  j => eq_2%jac_FD(:,:,:,0,0,0)
2221  d3j => eq_2%jac_FD(:,:,:,0,0,1)
2222  g13 => eq_2%g_FD(:,:,:,c([1,3],.true.),0,0,0)
2223  d3g13 => eq_2%g_FD(:,:,:,c([1,3],.true.),0,0,1)
2224  g23 => eq_2%g_FD(:,:,:,c([2,3],.true.),0,0,0)
2225  d3g23 => eq_2%g_FD(:,:,:,c([2,3],.true.),0,0,1)
2226  g33 => eq_2%g_FD(:,:,:,c([3,3],.true.),0,0,0)
2227  d3g33 => eq_2%g_FD(:,:,:,c([3,3],.true.),0,0,1)
2228  h12 => eq_2%h_FD(:,:,:,c([1,2],.true.),0,0,0)
2229  d3h12 => eq_2%h_FD(:,:,:,c([1,2],.true.),0,0,1)
2230  h22 => eq_2%h_FD(:,:,:,c([2,2],.true.),0,0,0)
2231  d1h22 => eq_2%h_FD(:,:,:,c([2,2],.true.),1,0,0)
2232  d3h22 => eq_2%h_FD(:,:,:,c([2,2],.true.),0,0,1)
2233  d13h22 => eq_2%h_FD(:,:,:,c([2,2],.true.),1,0,1)
2234  d33h22 => eq_2%h_FD(:,:,:,c([2,2],.true.),0,0,2)
2235  h23 => eq_2%h_FD(:,:,:,c([2,3],.true.),0,0,0)
2236  d1h23 => eq_2%h_FD(:,:,:,c([2,3],.true.),1,0,0)
2237  d3h23 => eq_2%h_FD(:,:,:,c([2,3],.true.),0,0,1)
2238  d13h23 => eq_2%h_FD(:,:,:,c([2,3],.true.),1,0,1)
2239  d33h23 => eq_2%h_FD(:,:,:,c([2,3],.true.),0,0,2)
2240 
2241  ! set up helper variables in eq grid
2242  if (use_pol_flux_f) then
2243  nm = x%n
2244  djq = eq_1%q_saf_FD(:,1)
2245  else
2246  djq = -eq_1%rot_t_FD(:,1)
2247  nm = x%m
2248  end if
2249  theta_3 = h23/h22
2250  select case (eq_style)
2251  case (1) ! VMEC
2252  d1theta_3 = d1h23/h22 - h23*d1h22/(h22**2)
2253  d3theta_3 = d3h23/h22 - h23*d3h22/(h22**2)
2254  d13theta_3 = d13h23/h22 - (d3h23*d1h22+d1h23*d3h22)/(h22**2) &
2255  &- h23*d13h22/(h22**2) + 2*h23*d3h22*d1h22/(h22**3)
2256  d33theta_3 = d33h23/h22 - 2*d3h23*d3h22/(h22**2) &
2257  &- h23*d33h22/(h22**2) + 2*h23*d3h22**2/(h22**3)
2258  case (2) ! HELENA
2259  ! geodesic derivative
2260  if (grid_eq%n(2).gt.norm_disc_prec_x+3) then ! need enough terms
2261  do kd = 1,grid_eq%loc_n_r
2262  do id = 1,grid_eq%n(1)
2263  ierr = spline(ang_geo_f(id,:,kd),theta_3(id,:,kd),&
2264  &ang_geo_f(id,:,kd),d1theta_3(id,:,kd),&
2265  &ord=norm_disc_prec_x,deriv=1)
2266  chckerr('')
2267  end do
2268  end do
2269  else
2270  d1theta_3 = 0._dp
2271  end if
2272  ! parallel derivative
2273  if (grid_eq%n(1).gt.norm_disc_prec_x+3) then ! need enough terms
2274  do kd = 1,grid_eq%loc_n_r
2275  do jd = 1,grid_eq%n(2)
2276  ierr = spline(ang_par_f(:,jd,kd),theta_3(:,jd,kd),&
2277  &ang_par_f(:,jd,kd),d3theta_3(:,jd,kd),&
2278  &ord=norm_disc_prec_x,deriv=1)
2279  chckerr('')
2280  end do
2281  end do
2282  else
2283  d3theta_3 = 0._dp
2284  end if
2285  end select
2286 
2287  ! set up U factors in eq grid
2288  t1(:,:,:,1) = g13/g33
2289  do kd = 1,grid_eq%loc_n_r
2290  t2(:,:,kd,1) = h12(:,:,kd)/h22(:,:,kd) + djq(kd)*ang_par_f(:,:,kd)
2291  t3(:,:,kd,1) = t1(:,:,kd,1)*djq(kd) + &
2292  &j(:,:,kd)**2*vac_perm*eq_1%pres_FD(kd,1)/g33(:,:,kd)
2293  t4(:,:,kd,1) = t1(:,:,kd,1)*djq(kd)*ang_par_f(:,:,kd) - &
2294  &g23(:,:,kd)/g33(:,:,kd)
2295  end do
2296  t5(:,:,:,1) = t1(:,:,:,1)*d3theta_3 - d1theta_3
2297  t6(:,:,:,1) = t1(:,:,:,1)*theta_3
2298 
2299  ! set up parallel derivatives of U in eq grid if VMEC is used
2300  if (eq_style.eq.1) then
2301  t1(:,:,:,2) = d3g13/g33 - g13*d3g33/g33**2
2302  do kd = 1,grid_eq%loc_n_r
2303  t2(:,:,kd,2) = d3h12(:,:,kd)/h22(:,:,kd) - &
2304  &h12(:,:,kd)*d3h22(:,:,kd)/h22(:,:,kd)**2 + djq(kd)
2305  t3(:,:,kd,2) = t1(:,:,kd,2)*djq(kd) + &
2306  &j(:,:,kd)*vac_perm*eq_1%pres_FD(kd,1)/g33(:,:,kd) * &
2307  &(2*d3j(:,:,kd)-d3g33(:,:,kd)*j(:,:,kd)/g33(:,:,kd))
2308  t4(:,:,kd,2) = (t1(:,:,kd,2)*ang_par_f(:,:,kd)+t1(:,:,kd,1))*&
2309  &djq(kd) - d3g23(:,:,kd)/g33(:,:,kd) + &
2310  &g23(:,:,kd)*d3g33(:,:,kd)/g33(:,:,kd)**2
2311  end do
2312  t5(:,:,:,2) = t1(:,:,:,2)*d3theta_3 + t1(:,:,:,1)*d33theta_3 - &
2313  &d13theta_3
2314  t6(:,:,:,2) = t1(:,:,:,2)*theta_3 + t1(:,:,:,1)*d3theta_3
2315  end if
2316 
2317  select case (x_grid_style)
2318  case (1) ! equilibrium
2319  ! point
2320  q_saf => eq_1%q_saf_FD(:,0)
2321  rot_t => eq_1%rot_t_FD(:,0)
2322  t1_x => t1
2323  t2_x => t2
2324  t3_x => t3
2325  t4_x => t4
2326  t5_x => t5
2327  t6_x => t6
2328  case (2,3) ! solution or enriched
2329  ! interpolate
2330  ierr = spline(grid_eq%loc_r_F,eq_1%q_saf_FD(:,0),&
2331  &grid_x%loc_r_F,q_saf,ord=norm_disc_prec_x)
2332  chckerr('')
2333  ierr = spline(grid_eq%loc_r_F,eq_1%rot_t_FD(:,0),&
2334  &grid_x%loc_r_F,rot_t,ord=norm_disc_prec_x)
2335  chckerr('')
2336  do ld = 1,t_size
2337  do jd = 1,grid_x%n(2)
2338  do id = 1,grid_x%n(1)
2339  ierr = spline(grid_eq%loc_r_F,t1(id,jd,:,ld),&
2340  &grid_x%loc_r_F,t1_x(id,jd,:,ld),&
2341  &ord=norm_disc_prec_x)
2342  chckerr('')
2343  ierr = spline(grid_eq%loc_r_F,t2(id,jd,:,ld),&
2344  &grid_x%loc_r_F,t2_x(id,jd,:,ld),&
2345  &ord=norm_disc_prec_x)
2346  chckerr('')
2347  ierr = spline(grid_eq%loc_r_F,t3(id,jd,:,ld),&
2348  &grid_x%loc_r_F,t3_x(id,jd,:,ld),&
2349  &ord=norm_disc_prec_x)
2350  chckerr('')
2351  ierr = spline(grid_eq%loc_r_F,t4(id,jd,:,ld),&
2352  &grid_x%loc_r_F,t4_x(id,jd,:,ld),&
2353  &ord=norm_disc_prec_x)
2354  chckerr('')
2355  ierr = spline(grid_eq%loc_r_F,t5(id,jd,:,ld),&
2356  &grid_x%loc_r_F,t5_x(id,jd,:,ld),&
2357  &ord=norm_disc_prec_x)
2358  chckerr('')
2359  ierr = spline(grid_eq%loc_r_F,t6(id,jd,:,ld),&
2360  &grid_x%loc_r_F,t6_x(id,jd,:,ld),&
2361  &ord=norm_disc_prec_x)
2362  chckerr('')
2363  end do
2364  end do
2365  end do
2366  end select
2367 
2368  ! set up n_frac
2369  do kd = 1,grid_x%loc_n_r
2370  if (use_pol_flux_f) then
2371  n_frac(kd,:) = x%n(kd,:)*q_saf(kd) - x%m(kd,:)
2372  else
2373  n_frac(kd,:) = -x%m(kd,:)*rot_t(kd) + x%n(kd,:)
2374  end if
2375  end do
2376 
2377  ! calculate U_0 and U_1 for all modes
2378  do ld = 1,x%n_mod
2379  ! calculate order 1 of U_0 and U_1
2380  if (u_style.ge.1) then
2381  x%U_0(:,:,:,ld) = -t2_x(:,:,:,1)
2382  do kd = 1,grid_x%loc_n_r
2383  x%U_1(:,:,kd,ld) = iu/nm(kd,ld)
2384  end do
2385  end if
2386  ! include order 2
2387  if (u_style.ge.2) then
2388  do kd = 1,grid_x%loc_n_r
2389  u_corr(:,:,kd,1) = t3_x(:,:,kd,1) + &
2390  &iu*n_frac(kd,ld)*t4_x(:,:,kd,1)
2391  u_corr(:,:,kd,2) = n_frac(kd,ld)/nm(kd,ld)*t1_x(:,:,kd,1)
2392  x%U_0(:,:,kd,ld) = x%U_0(:,:,kd,ld) + iu/nm(kd,ld)*&
2393  &u_corr(:,:,kd,1)
2394  x%U_1(:,:,kd,ld) = x%U_1(:,:,kd,ld) + iu/nm(kd,ld)*&
2395  &u_corr(:,:,kd,2)
2396  end do
2397  end if
2398  ! include order 3
2399  if (u_style.ge.3) then
2400  do kd = 1,grid_x%loc_n_r
2401  u_corr(:,:,kd,1) = iu*n_frac(kd,ld)*&
2402  &(-t5_x(:,:,kd,1)-iu*n_frac(kd,ld)*t6_x(:,:,kd,1))
2403  x%U_0(:,:,kd,ld) = x%U_0(:,:,kd,ld) + (iu/nm(kd,ld))**2*&
2404  &u_corr(:,:,kd,1)
2405  end do
2406  end if
2407  if (u_style.ge.4) then
2408  call writo('The geodesic perturbation U is implemented up to &
2409  &order 3',warning=.true.)
2410  end if
2411  end do
2412 
2413  ! calculate DU_0 and DU_1, starting with first part
2414  select case (eq_style)
2415  case (1) ! VMEC
2416  do ld = 1,x%n_mod
2417  ! calculate order 1 of DU_0 and DU_1
2418  if (u_style.ge.1) then
2419  x%DU_0(:,:,:,ld) = -t2_x(:,:,:,2)
2420  x%DU_1(:,:,:,ld) = 0._dp
2421  end if
2422  ! include order 2
2423  if (u_style.ge.2) then
2424  do kd = 1,grid_x%loc_n_r
2425  u_corr(:,:,kd,1) = t3_x(:,:,kd,2) + &
2426  &iu*n_frac(kd,ld)*t4_x(:,:,kd,2)
2427  u_corr(:,:,kd,2) = n_frac(kd,ld)/nm(kd,ld)*&
2428  &t1_x(:,:,kd,2)
2429  x%DU_0(:,:,kd,ld) = x%DU_0(:,:,kd,ld) + &
2430  &iu/nm(kd,ld)*u_corr(:,:,kd,1)
2431  x%DU_1(:,:,kd,ld) = x%DU_1(:,:,kd,ld) + &
2432  &iu/nm(kd,ld)*u_corr(:,:,kd,2)
2433  end do
2434  end if
2435  ! include order 3
2436  if (u_style.ge.3) then
2437  do kd = 1,grid_x%loc_n_r
2438  u_corr(:,:,kd,1) = iu*n_frac(kd,ld)*&
2439  &(-t5_x(:,:,kd,2)-iu*n_frac(kd,ld)*&
2440  &t6_x(:,:,kd,2))
2441  x%DU_0(:,:,kd,ld) = x%DU_0(:,:,kd,ld) + &
2442  &(iu/nm(kd,ld))**2*u_corr(:,:,kd,1)
2443  end do
2444  end if
2445  if (u_style.ge.4) then
2446  call writo('The geodesic perturbation U is implemented &
2447  &only up to order 3',warning=.true.)
2448  end if
2449  end do
2450 #if ldebug
2451  if (ltest) then
2452  call writo('Test calculation of DU')
2453  if(get_log(.false.,ind=.true.)) then
2454  ierr = test_du()
2455  chckerr('')
2456  call pause_prog(ind=.true.)
2457  end if
2458  end if
2459 #endif
2460  case (2) ! HELENA
2461  ! reset pointers
2462  if (use_pol_flux_f) then
2463  ang_par_f => grid_x%theta_F
2464  else
2465  ang_par_f => grid_x%zeta_F
2466  end if
2467  ! set up helper variable for derivative
2468  ! numerically derive U_0 and U_1
2469  do ld = 1,x%n_mod
2470  do kd = 1,grid_x%loc_n_r
2471  do jd = 1,grid_x%n(2)
2472  ierr = spline(ang_par_f(:,jd,kd),x%U_0(:,jd,kd,ld),&
2473  &ang_par_f(:,jd,kd),x%DU_0(:,jd,kd,ld),&
2474  &ord=norm_disc_prec_x,deriv=1)
2475  chckerr('')
2476  ierr = spline(ang_par_f(:,jd,kd),x%U_1(:,jd,kd,ld),&
2477  &ang_par_f(:,jd,kd),x%DU_1(:,jd,kd,ld),&
2478  &ord=norm_disc_prec_x,deriv=1)
2479  chckerr('')
2480  end do
2481  end do
2482  end do
2483  end select
2484 
2485  ! add second part of derivatives ~ n_frac
2486  do kd = 1,grid_x%loc_n_r
2487  do ld = 1,x%n_mod
2488  x%DU_0(:,:,kd,ld) = x%DU_0(:,:,kd,ld) + &
2489  &iu*n_frac(kd,ld)*x%U_0(:,:,kd,ld)
2490  x%DU_1(:,:,kd,ld) = x%DU_1(:,:,kd,ld) + &
2491  &iu*n_frac(kd,ld)*x%U_1(:,:,kd,ld)
2492  end do
2493  end do
2494 
2495  ! clean up
2496  nullify(ang_par_f,ang_geo_f)
2497  nullify(j,d3j)
2498  nullify(g13,d3g13)
2499  nullify(g23,d3g23)
2500  nullify(g33,d3g33)
2501  nullify(h12,d3h12)
2502  nullify(h22,d1h22,d3h22,d13h22,d33h22)
2503  nullify(h23,d1h23,d3h23,d13h23,d33h23)
2504  select case (x_grid_style)
2505  case (1) ! equilibrium
2506  ! do nothing
2507  case (2,3) ! solution or enriched
2508  deallocate(q_saf,rot_t,t1_x,t2_x,t3_x,t4_x,t5_x,t6_x)
2509  end select
2510  nullify(q_saf,rot_t,t1_x,t2_x,t3_x,t4_x,t5_x,t6_x)
2511 #if ldebug
2512  contains
2513  ! Test calculation of DU by deriving U numerically.
2514  integer function test_du() result(ierr)
2515  use num_vars, only: norm_disc_prec_x
2516 
2517  character(*), parameter :: rout_name = 'test_DU'
2518 
2519  ! local variables
2520  integer :: jd, kd, ld ! counters
2521  complex(dp), allocatable :: du_0(:,:,:) ! alternative calculation for DU_0
2522  complex(dp), allocatable :: du_1(:,:,:) ! alternative calculation for DU_1
2523  character(len=max_str_ln) :: file_name ! name of plot file
2524  character(len=max_str_ln) :: description ! description of plot
2525 
2526  ! initialize ierr
2527  ierr = 0
2528 
2529  ! warning
2530  call writo('This test is representable only if there are enough &
2531  &parallel points in the grid!',warning=.true.)
2532 
2533  ! output
2534  call writo('Going to test whether DU is consistent with U')
2535  call lvl_ud(1)
2536 
2537  ! set up DU_0 and DU_1
2538  allocate(du_0(grid_x%n(1),grid_x%n(2),grid_x%loc_n_r))
2539  allocate(du_1(grid_x%n(1),grid_x%n(2),grid_x%loc_n_r))
2540 
2541  ! reset pointers
2542  if (use_pol_flux_f) then
2543  ang_par_f => grid_x%theta_F
2544  else
2545  ang_par_f => grid_x%zeta_F
2546  end if
2547 
2548  ! loop over all modes
2549  do ld = 1,x%n_mod
2550  ! derivate numerically
2551  do kd = 1,grid_x%loc_n_r
2552  do jd = 1,grid_x%n(2)
2553  ierr = spline(ang_par_f(:,jd,kd),x%U_0(:,jd,kd,ld),&
2554  &ang_par_f(:,jd,kd),du_0(:,jd,kd),&
2555  &ord=norm_disc_prec_x,deriv=1)
2556  chckerr('')
2557  ierr = spline(ang_par_f(:,jd,kd),x%U_1(:,jd,kd,ld),&
2558  &ang_par_f(:,jd,kd),du_1(:,jd,kd),&
2559  &ord=norm_disc_prec_x,deriv=1)
2560  chckerr('')
2561  end do
2562  end do
2563 
2564  ! set some variables
2565  file_name = 'TEST_RE_DU_0_'//trim(i2str(ld))
2566  description = 'Verifying DU_0 by deriving U_0 numerically for &
2567  &mode '//trim(i2str(ld)//', real part')
2568 
2569  ! plot difference for RE DU_0
2570  call plot_diff_hdf5(rp(du_0),rp(x%DU_0(:,:,:,ld)),&
2571  &file_name,descr=description,output_message=.true.)
2572 
2573  ! set some variables
2574  file_name = 'TEST_IM_DU_0_'//trim(i2str(ld))
2575  description = 'Verifying DU_0 by deriving U_0 numerically for &
2576  &mode '//trim(i2str(ld)//', imaginary part')
2577 
2578  ! plot difference for IM DU_0
2579  call plot_diff_hdf5(ip(du_0),ip(x%DU_0(:,:,:,ld)),&
2580  &file_name,descr=description,output_message=.true.)
2581 
2582  ! set some variables
2583  file_name = 'TEST_RE_DU_1_'//trim(i2str(ld))
2584  description = 'Verifying DU_1 by deriving U_1 numerically for &
2585  &mode '//trim(i2str(ld)//', real part')
2586 
2587  ! plot difference for RE DU_1
2588  call plot_diff_hdf5(rp(du_1),rp(x%DU_1(:,:,:,ld)),&
2589  &file_name,descr=description,output_message=.true.)
2590 
2591  ! set some variables
2592  file_name = 'TEST_IM_DU_1_'//trim(i2str(ld))
2593  description = 'Verifying DU_1 by deriving U_1 numerically for &
2594  &mode '//trim(i2str(ld)//', imaginary part')
2595 
2596  ! plot difference for IM DU_1
2597  call plot_diff_hdf5(ip(du_1),ip(x%DU_1(:,:,:,ld)),&
2598  &file_name,descr=description,output_message=.true.)
2599  end do
2600 
2601  ! user output
2602  call lvl_ud(-1)
2603  call writo('Test complete')
2604  end function test_du
2605 #endif
2606  end function calc_u
2607 
2644  integer function calc_pv(grid_eq,grid_X,eq_1,eq_2,X_a,X_b,X,lim_sec_X) &
2645  &result(ierr)
2646  use num_vars, only: use_pol_flux_f, norm_disc_prec_x, x_grid_style
2647  use eq_vars, only: vac_perm
2648  use num_utilities, only: c, spline
2649  use x_utilities, only: is_necessary_x
2650  use x_vars, only: n_mod_x
2651 
2652  character(*), parameter :: rout_name = 'calc_PV'
2653 
2654  ! use input / output
2655  type(grid_type), intent(in) :: grid_eq
2656  type(grid_type), intent(in) :: grid_x
2657  type(eq_1_type), intent(in), target :: eq_1
2658  type(eq_2_type), intent(in), target :: eq_2
2659  type(x_1_type), intent(in) :: x_a
2660  type(x_1_type), intent(in) :: x_b
2661  type(x_2_type), intent(inout) :: x
2662  integer, intent(in), optional :: lim_sec_x(2,2)
2663 
2664  ! local variables
2665  integer :: m, k ! counters
2666  integer :: id, jd, kd ! counters
2667  integer :: c_loc(2) ! local c for symmetric and asymmetric variables
2668  logical :: calc_this(2) ! whether this combination needs to be calculated
2669 
2670  ! jacobian
2671  real(dp), pointer :: j(:,:,:) ! jac
2672  ! lower metric factors
2673  real(dp), pointer :: g33(:,:,:) ! h_theta,theta or h_zeta,zeta
2674  ! upper metric factors
2675  real(dp), pointer :: h22(:,:,:) ! h^psi,psi
2676  ! helper variables
2677  real(dp), pointer :: q_saf(:), rot_t(:) ! safety factor and rotational transform in X grid
2678  real(dp), allocatable :: fac_n(:), fac_m(:) ! multiplicative factors for n and m
2679  ! PV factors
2680  real(dp), allocatable, target :: t1(:,:,:) ! h22/mu_0 g33
2681  real(dp), allocatable, target :: t2(:,:,:) ! JS + mu_0 sigma g33/J h22
2682  real(dp), allocatable, target :: t3(:,:,:) ! sigma/J T2
2683  real(dp), allocatable, target :: t4(:,:,:) ! 1/mu_0 J^2 h22
2684  real(dp), allocatable, target :: t5(:,:,:) ! 2 p' kappa_n
2685  real(dp), pointer :: t1_x(:,:,:) ! T1 in X grid and par. deriv.
2686  real(dp), pointer :: t2_x(:,:,:) ! T2 in X grid and par. deriv.
2687  real(dp), pointer :: t3_x(:,:,:) ! T3 in X grid and par. deriv.
2688  real(dp), pointer :: t4_x(:,:,:) ! T4 in X grid and par. deriv.
2689  real(dp), pointer :: t5_x(:,:,:) ! T5 in X grid and par. deriv.
2690 
2691  ! initialize ierr
2692  ierr = 0
2693 
2694  ! allocate variables
2695  ! helper variables
2696  allocate(q_saf(grid_x%loc_n_r))
2697  allocate(rot_t(grid_x%loc_n_r))
2698  allocate(fac_n(grid_x%loc_n_r),fac_m(grid_x%loc_n_r))
2699  ! PV factors
2700  allocate(t1(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r))
2701  allocate(t2(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r))
2702  allocate(t3(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r))
2703  allocate(t4(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r))
2704  allocate(t5(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r))
2705  select case (x_grid_style)
2706  case (1) ! equilibrium
2707  ! do nothing
2708  case (2,3) ! solution or enriched
2709  allocate(q_saf(grid_x%loc_n_r))
2710  allocate(rot_t(grid_x%loc_n_r))
2711  allocate(t1_x(grid_x%n(1),grid_x%n(2),grid_x%loc_n_r))
2712  allocate(t2_x(grid_x%n(1),grid_x%n(2),grid_x%loc_n_r))
2713  allocate(t3_x(grid_x%n(1),grid_x%n(2),grid_x%loc_n_r))
2714  allocate(t4_x(grid_x%n(1),grid_x%n(2),grid_x%loc_n_r))
2715  allocate(t5_x(grid_x%n(1),grid_x%n(2),grid_x%loc_n_r))
2716  end select
2717 
2718  ! set pointers
2719  j => eq_2%jac_FD(:,:,:,0,0,0)
2720  g33 => eq_2%g_FD(:,:,:,c([3,3],.true.),0,0,0)
2721  h22 => eq_2%h_FD(:,:,:,c([2,2],.true.),0,0,0)
2722 
2723  ! set up PV factors in eq grid
2724  t1 = h22/(vac_perm*g33)
2725  t2 = j*eq_2%S + vac_perm*eq_2%sigma*g33/(j*h22)
2726  t3 = eq_2%sigma/j*t2
2727  t4 = 1._dp/(vac_perm*j**2*h22)
2728  do kd = 1,grid_eq%loc_n_r
2729  t5(:,:,kd) = 2*eq_1%pres_FD(kd,1)*eq_2%kappa_n(:,:,kd)
2730  end do
2731 
2732  select case (x_grid_style)
2733  case (1) ! equilibrium
2734  ! point
2735  q_saf => eq_1%q_saf_FD(:,0)
2736  rot_t => eq_1%rot_t_FD(:,0)
2737  t1_x => t1
2738  t2_x => t2
2739  t3_x => t3
2740  t4_x => t4
2741  t5_x => t5
2742  case (2,3) ! solution or enriched
2743  ! interpolate
2744  ierr = spline(grid_eq%loc_r_F,eq_1%q_saf_FD(:,0),&
2745  &grid_x%loc_r_F,q_saf,ord=norm_disc_prec_x)
2746  chckerr('')
2747  ierr = spline(grid_eq%loc_r_F,eq_1%rot_t_FD(:,0),&
2748  &grid_x%loc_r_F,rot_t,ord=norm_disc_prec_x)
2749  chckerr('')
2750  do jd = 1,grid_x%n(2)
2751  do id = 1,grid_x%n(1)
2752  ierr = spline(grid_eq%loc_r_F,t1(id,jd,:),&
2753  &grid_x%loc_r_F,t1_x(id,jd,:),&
2754  &ord=norm_disc_prec_x)
2755  chckerr('')
2756  ierr = spline(grid_eq%loc_r_F,t2(id,jd,:),&
2757  &grid_x%loc_r_F,t2_x(id,jd,:),&
2758  &ord=norm_disc_prec_x)
2759  chckerr('')
2760  ierr = spline(grid_eq%loc_r_F,t3(id,jd,:),&
2761  &grid_x%loc_r_F,t3_x(id,jd,:),&
2762  &ord=norm_disc_prec_x)
2763  chckerr('')
2764  ierr = spline(grid_eq%loc_r_F,t4(id,jd,:),&
2765  &grid_x%loc_r_F,t4_x(id,jd,:),&
2766  &ord=norm_disc_prec_x)
2767  chckerr('')
2768  ierr = spline(grid_eq%loc_r_F,t5(id,jd,:),&
2769  &grid_x%loc_r_F,t5_x(id,jd,:),&
2770  &ord=norm_disc_prec_x)
2771  chckerr('')
2772  end do
2773  end do
2774  end select
2775 
2776  ! set up fac_n and fac_m
2777  if (use_pol_flux_f) then
2778  fac_n = q_saf
2779  fac_m = 1.0_dp
2780  else
2781  fac_n = 1.0_dp
2782  fac_m = rot_t
2783  end if
2784 
2785  ! loop over all modes
2786  do m = 1,x_b%n_mod
2787  do k = 1,x_a%n_mod
2788  ! check whether modes are correct
2789  do kd = 1,grid_x%loc_n_r
2790  if (x%n_1(kd,k).ne.x_a%n(kd,k) .or. &
2791  &x%n_2(kd,m).ne.x_b%n(kd,m) .or. &
2792  &x%m_1(kd,k).ne.x_a%m(kd,k) .or. &
2793  &x%m_2(kd,m).ne.x_b%m(kd,m)) then
2794  ierr = 1
2795  chckerr('Modes do not match')
2796  end if
2797  end do
2798 
2799  ! check whether mode combination needs to be calculated
2800  calc_this(1) = is_necessary_x(.true.,[k,m],lim_sec_x)
2801  calc_this(2) = is_necessary_x(.false.,[k,m],lim_sec_x)
2802 
2803  ! set up c_loc
2804  c_loc(1) = c([k,m],.true.,n_mod_x,lim_sec_x)
2805  c_loc(2) = c([k,m],.false.,n_mod_x,lim_sec_x)
2806 
2807  ! calculate PV_0
2808  if (calc_this(1)) then
2809  do kd = 1,grid_x%loc_n_r
2810  x%PV_0(:,:,kd,c_loc(1)) = t1_x(:,:,kd)*&
2811  &(x_b%DU_0(:,:,kd,m) - t2_x(:,:,kd)) * &
2812  &(conjg(x_a%DU_0(:,:,kd,k)) - t2_x(:,:,kd)) - &
2813  &t3_x(:,:,kd) - t5_x(:,:,kd) + t4_x(:,:,kd)*&
2814  &(x_b%n(kd,m)*fac_n(kd)-x_b%m(kd,m)*fac_m(kd))*&
2815  &(x_a%n(kd,k)*fac_n(kd)-x_a%m(kd,k)*fac_m(kd))
2816  end do
2817  end if
2818 
2819  ! calculate PV_1
2820  if (calc_this(2)) then
2821  x%PV_1(:,:,:,c_loc(2)) = t1_x * x_b%DU_1(:,:,:,m) * &
2822  &(conjg(x_a%DU_0(:,:,:,k))-t2_x)
2823  end if
2824 
2825  ! calculate PV_2
2826  if (calc_this(1)) then
2827  x%PV_2(:,:,:,c_loc(1)) = t1_x * x_b%DU_1(:,:,:,m) * &
2828  &conjg(x_a%DU_1(:,:,:,k))
2829  end if
2830  end do
2831  end do
2832 
2833  ! clean up
2834  nullify(j)
2835  nullify(g33)
2836  nullify(h22)
2837  select case (x_grid_style)
2838  case (1) ! equilibrium
2839  ! do nothing
2840  case (2,3) ! solution or enriched
2841  deallocate(q_saf,rot_t,t1_x,t2_x,t3_x,t4_x,t5_x)
2842  end select
2843  nullify(q_saf,rot_t,t1_x,t2_x,t3_x,t4_x,t5_x)
2844  end function calc_pv
2845 
2877  integer function calc_kv(grid_eq,grid_X,eq_1,eq_2,X_a,X_b,X,lim_sec_X) &
2878  &result(ierr)
2879  use num_vars, only: k_style, norm_disc_prec_x, x_grid_style
2880  use num_utilities, only: c, spline
2881  use x_utilities, only: is_necessary_x
2882  use x_vars, only: n_mod_x
2883 
2884  character(*), parameter :: rout_name = 'calc_KV'
2885 
2886  ! use input / output
2887  type(grid_type), intent(in) :: grid_eq
2888  type(grid_type), intent(in) :: grid_x
2889  type(eq_1_type), intent(in) :: eq_1
2890  type(eq_2_type), intent(in), target :: eq_2
2891  type(x_1_type), intent(in) :: x_a
2892  type(x_1_type), intent(in) :: x_b
2893  type(x_2_type), intent(inout) :: x
2894  integer, intent(in), optional :: lim_sec_x(2,2)
2895 
2896  ! local variables
2897  integer :: m, k ! counters
2898  integer :: id, jd, kd ! counters
2899  integer :: c_loc(2) ! local c for symmetric and asymmetric variables
2900  logical :: calc_this(2) ! whether this combination needs to be calculated
2901 
2902  ! jacobian
2903  real(dp), pointer :: j(:,:,:) ! jac
2904  ! lower metric factors
2905  real(dp), pointer :: g33(:,:,:) ! h_theta,theta or h_zeta,zeta
2906  ! upper metric factors
2907  real(dp), pointer :: h22(:,:,:) ! h^psi,psi
2908  ! KV factors
2909  real(dp), allocatable, target :: t1(:,:,:) ! h22/mu_0 J g33
2910  real(dp), allocatable, target :: t2(:,:,:) ! JS + mu_0 sigma g33/h22
2911  real(dp), pointer :: t1_x(:,:,:) ! T1 in X grid and par. deriv.
2912  real(dp), pointer :: t2_x(:,:,:) ! T2 in X grid and par. deriv.
2913 
2914  ! initialize ierr
2915  ierr = 0
2916 
2917  ! allocate variables
2918  ! KV factors
2919  allocate(t1(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r))
2920  allocate(t2(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r))
2921  select case (x_grid_style)
2922  case (1) ! equilibrium
2923  ! do nothing
2924  case (2,3) ! solution or enriched
2925  allocate(t1_x(grid_x%n(1),grid_x%n(2),grid_x%loc_n_r))
2926  allocate(t2_x(grid_x%n(1),grid_x%n(2),grid_x%loc_n_r))
2927  end select
2928 
2929  ! set pointers
2930  j => eq_2%jac_FD(:,:,:,0,0,0)
2931  g33 => eq_2%g_FD(:,:,:,c([3,3],.true.),0,0,0)
2932  h22 => eq_2%h_FD(:,:,:,c([2,2],.true.),0,0,0)
2933 
2934  ! set up KV factors in eq grid
2935  do kd = 1,grid_eq%loc_n_r
2936  t1(:,:,kd) = eq_1%rho(kd)*j(:,:,kd)**2*h22(:,:,kd)/g33(:,:,kd)
2937  t2(:,:,kd) = eq_1%rho(kd)/h22(:,:,kd)
2938  end do
2939 
2940  select case (x_grid_style)
2941  case (1) ! equilibrium
2942  ! point
2943  t1_x => t1
2944  t2_x => t2
2945  case (2,3) ! solution or enriched
2946  ! interpolate
2947  do jd = 1,grid_x%n(2)
2948  do id = 1,grid_x%n(1)
2949  ierr = spline(grid_eq%loc_r_F,t1(id,jd,:),&
2950  &grid_x%loc_r_F,t1_x(id,jd,:),&
2951  &ord=norm_disc_prec_x)
2952  chckerr('')
2953  ierr = spline(grid_eq%loc_r_F,t2(id,jd,:),&
2954  &grid_x%loc_r_F,t2_x(id,jd,:),&
2955  &ord=norm_disc_prec_x)
2956  chckerr('')
2957  end do
2958  end do
2959  end select
2960 
2961  ! loop over all modes
2962  do m = 1,x_b%n_mod
2963  do k = 1,x_a%n_mod
2964  ! check whether modes are correct
2965  do kd = 1,grid_x%loc_n_r
2966  if (x%n_1(kd,k).ne.x_a%n(kd,k) .or. &
2967  &x%n_2(kd,m).ne.x_b%n(kd,m) .or. &
2968  &x%m_1(kd,k).ne.x_a%m(kd,k) .or. &
2969  &x%m_2(kd,m).ne.x_b%m(kd,m)) then
2970  ierr = 1
2971  chckerr('Modes do not match')
2972  end if
2973  end do
2974 
2975  ! check whether mode combination needs to be calculated
2976  calc_this(1) = is_necessary_x(.true.,[k,m],lim_sec_x)
2977  calc_this(2) = is_necessary_x(.false.,[k,m],lim_sec_x)
2978 
2979  ! set up c_loc
2980  c_loc(1) = c([k,m],.true.,n_mod_x,lim_sec_x)
2981  c_loc(2) = c([k,m],.false.,n_mod_x,lim_sec_x)
2982 
2983  select case (k_style)
2984  case (1) ! normalization of full perpendicular component
2985  ! calculate KV_0
2986  if (calc_this(1)) then
2987  x%KV_0(:,:,:,c_loc(1)) = t2_x + t1_x * &
2988  &x_b%U_0(:,:,:,m) * conjg(x_a%U_0(:,:,:,k))
2989  end if
2990 
2991  ! calculate KV_1
2992  if (calc_this(2)) then
2993  x%KV_1(:,:,:,c_loc(2)) = t1_x * &
2994  &x_b%U_1(:,:,:,m) * conjg(x_a%U_0(:,:,:,k))
2995  end if
2996 
2997  ! calculate KV_2
2998  if (calc_this(1)) then
2999  x%KV_2(:,:,:,c_loc(1)) = t1_x * &
3000  &x_b%U_1(:,:,:,m) * conjg(x_a%U_1(:,:,:,k))
3001  end if
3002  case (2) ! normalization of only normal component
3003  ! calculate KV_0
3004  if (calc_this(1)) then
3005  x%KV_0(:,:,:,c_loc(1)) = t2_x
3006  end if
3007 
3008  ! calculate KV_1
3009  if (calc_this(2)) then
3010  x%KV_1(:,:,:,c_loc(2)) = 0._dp
3011  end if
3012 
3013  ! calculate KV_2
3014  if (calc_this(1)) then
3015  x%KV_2(:,:,:,c_loc(1)) = 0._dp
3016  end if
3017  end select
3018  end do
3019  end do
3020 
3021  ! clean up
3022  nullify(j)
3023  nullify(g33)
3024  nullify(h22)
3025  select case (x_grid_style)
3026  case (1) ! equilibrium
3027  ! do nothing
3028  case (2,3) ! solution or enriched
3029  deallocate(t1_x,t2_x)
3030  end select
3031  nullify(t1_x,t2_x)
3032  end function calc_kv
3033 
3079  integer function calc_magn_ints(grid_eq,grid_X,eq,X,X_int,prev_style,&
3080  &lim_sec_X) result(ierr)
3081  use num_vars, only: use_pol_flux_f, norm_disc_prec_x, magn_int_style, &
3082  &alpha_style, x_grid_style
3083  use x_utilities, only: is_necessary_x
3084  use x_vars, only: n_mod_x
3085  use num_utilities, only: c, spline
3086  use grid_vars, only: min_par_x, max_par_x, n_alpha
3087  use rich_vars, only: rich_lvl
3088 
3089  character(*), parameter :: rout_name = 'calc_magn_ints'
3090 
3091  ! input / output
3092  type(grid_type), intent(in) :: grid_eq
3093  type(grid_type), intent(in) :: grid_x
3094  type(eq_2_type), intent(in), target :: eq
3095  type(x_2_type), intent(in) :: x
3096  type(x_2_type), intent(inout) :: x_int
3097  integer, intent(in), optional :: prev_style
3098  integer, intent(in), optional :: lim_sec_x(2,2)
3099 
3100  ! local variables, not used in subroutines
3101  logical :: calc_this(2) ! whether this combination needs to be calculated
3102  integer :: k, m ! counters
3103  integer :: id, jd, kd ! counters
3104  integer :: c_loc(2) ! local c for symmetric and asymmetric variables
3105  integer :: c_tot(2) ! total c for symmetric and asymmetric variables
3106  integer :: prev_style_loc ! local prev_style
3107  real(dp) :: alpha_fac ! factor for alpha (style 1: Weyl factor, style 2: step size)
3108  real(dp) :: prev_mult_fac ! multiplicative factor for previous style
3109  real(dp), allocatable :: step_size(:,:) ! step size for every field line and normal point
3110  real(dp), pointer :: ang_par_f(:,:,:) => null() ! parallel angle in flux coordinates
3111  complex(dp), allocatable :: v_int_work(:,:) ! work V_int
3112 
3113  ! local variables, also used in subroutines
3114  integer :: nr_int_regions ! nr. of integration regions
3115  integer, allocatable :: int_dims(:,:) ! dimension information of integration regions
3116  real(dp), allocatable :: int_facs(:) ! integration factors for each of the regions
3117  complex(dp), allocatable :: j_exp_ang(:,:,:) ! J * exponential of Flux parallel angle
3118 
3119  ! Makes use of nr_int_regions, int_dims, int_facs and J_exp_ang
3120 
3121  ! jacobian
3122  real(dp), pointer :: j(:,:,:) ! jac
3123 
3124  ! initialize ierr
3125  ierr = 0
3126 
3127  ! user output
3128  call writo('Calculate field-line averages')
3129  call lvl_ud(1)
3130 
3131  ! set up local prev_style
3132  prev_style_loc = 0
3133  if (present(prev_style)) prev_style_loc = prev_style
3134 
3135  ! allocate variables
3136  allocate(j_exp_ang(grid_x%n(1),grid_x%n(2),grid_x%loc_n_r))
3137  select case (x_grid_style)
3138  case (1) ! equilibrium
3139  ! point
3140  j => eq%jac_FD(:,:,:,0,0,0)
3141  case (2,3) ! solution or enriched
3142  ! allocate
3143  allocate(j(grid_eq%n(1),grid_eq%n(2),grid_x%loc_n_r))
3144 
3145  ! interpolate
3146  do jd = 1,grid_x%n(2)
3147  do id = 1,grid_x%n(1)
3148  ierr = spline(grid_eq%loc_r_F,eq%jac_FD(id,jd,:,0,0,0),&
3149  &grid_x%loc_r_F,j(id,jd,:),ord=norm_disc_prec_x)
3150  chckerr('')
3151  end do
3152  end do
3153  end select
3154 
3155  ! set up parallel angle in flux coordinates
3156  if (use_pol_flux_f) then
3157  ang_par_f => grid_x%theta_F
3158  else
3159  ang_par_f => grid_x%zeta_F
3160  end if
3161 
3162  ! set up alpha factor
3163  select case (alpha_style)
3164  case (1) ! single field line, multiple turns
3165  alpha_fac = (2*pi)**2/((max_par_x-min_par_x)*pi) ! Weyl factor
3166  case (2) ! multiple field lines, single turns
3167  alpha_fac = (2*pi)/n_alpha ! grid size
3168  end select
3169 
3170  ! set up integration variables
3171 
3172  ! set nr_int_regions
3173  select case (magn_int_style)
3174  case (1) ! Trapezoidal rule
3175  select case (prev_style_loc)
3176  case (2,3) ! change indices
3177  nr_int_regions = 1
3178  case default ! don't change indices
3179  nr_int_regions = 2
3180  end select
3181  case (2) ! Simpson's 3/8 rule
3182  select case (prev_style_loc)
3183  case (2,3) ! change indices
3184  nr_int_regions = 3
3185  case default ! don't change indices
3186  nr_int_regions = 4
3187  end select
3188  end select
3189 
3190  ! set up integration factors and dimensions
3191  allocate(int_facs(nr_int_regions))
3192  allocate(int_dims(nr_int_regions,3))
3193  select case (magn_int_style)
3194  case (1) ! Trapezoidal rule
3195  select case (prev_style_loc)
3196  case (2,3) ! change indices
3197  int_facs = 0.5_dp*[2]
3198  int_dims(1,:) = [1,grid_x%n(1),1] ! region 1: all points
3199  case default ! don't change indices
3200  int_facs = 0.5_dp*[1,2]
3201  int_dims(1,:) = [1,grid_x%n(1),grid_x%n(1)-1] ! region 1: first and last points
3202  int_dims(2,:) = [2,grid_x%n(1)-1,1] ! region 2: intermediate points
3203  end select
3204  case (2) ! Simpson's 3/8 rule
3205  select case (prev_style_loc)
3206  case (2,3) ! change indices
3207  int_facs = 0.375_dp*[3,2,3]
3208  int_dims(1,:) = [1,grid_x%n(1)-2,3] ! region 1: intermediate points ~ 3
3209  int_dims(2,:) = [2,grid_x%n(1)-1,3] ! region 2: intermediate points ~ 2
3210  int_dims(3,:) = [3,grid_x%n(1),3] ! region 3: intermediate points ~ 3
3211  case default ! don't change indices
3212  int_facs = 0.375_dp*[1,3,3,2]
3213  int_dims(1,:) = [1,grid_x%n(1),grid_x%n(1)-1] ! region 1: first and last points
3214  int_dims(2,:) = [2,grid_x%n(1)-2,3] ! region 2: intermediate points ~ 3
3215  int_dims(3,:) = [3,grid_x%n(1)-1,3] ! region 3: intermediate points ~ 3
3216  int_dims(4,:) = [4,grid_x%n(1)-3,3] ! region 4: intermediate points ~ 2
3217  end select
3218  end select
3219 
3220  ! set up multiplicative factor for previous level
3221  select case (prev_style_loc)
3222  case (1,3) ! add to integral of previous equilibrium job
3223  prev_mult_fac = 1.0_dp
3224  case (2) ! add to integral of previous Richardson level / 2
3225  prev_mult_fac = 0.5_dp
3226  case default ! don't add, overwrite
3227  prev_mult_fac = 0.0_dp
3228  end select
3229 
3230  ! set up step size
3231  allocate(step_size(grid_x%n(2),grid_x%loc_n_r))
3232  step_size = ang_par_f(2,:,:)-ang_par_f(1,:,:)
3233  if (rich_lvl.gt.1) step_size = step_size*0.5_dp ! only half of points present: step size is actually half
3234 
3235  ! initialize local V_int
3236  allocate(v_int_work(grid_x%n(2),grid_x%loc_n_r))
3237 
3238  ! loop over all modes
3239  do m = 1,x%n_mod(2)
3240  do k = 1,x%n_mod(1)
3241  ! check whether mode combination needs to be calculated
3242  calc_this(1) = is_necessary_x(.true.,[k,m],lim_sec_x)
3243  calc_this(2) = is_necessary_x(.false.,[k,m],lim_sec_x)
3244 
3245  ! set up local and total indices
3246  c_loc(1) = c([k,m],.true.,n_mod_x,lim_sec_x)
3247  c_loc(2) = c([k,m],.false.,n_mod_x,lim_sec_x)
3248  c_tot(1) = c([lim_sec_x(1,1)-1+k,lim_sec_x(1,2)-1+m],.true.,&
3249  &n_mod_x)
3250  c_tot(2) = c([lim_sec_x(1,1)-1+k,lim_sec_x(1,2)-1+m],.false.,&
3251  &n_mod_x)
3252 
3253  ! calculate J_exp_ang
3254  do kd = 1,grid_x%loc_n_r
3255  if (use_pol_flux_f) then
3256  j_exp_ang(:,:,kd) = j(:,:,kd)*&
3257  &exp(iu*(x%m_1(kd,k)-x%m_2(kd,m))*&
3258  &ang_par_f(:,:,kd))
3259  else
3260  j_exp_ang(:,:,kd) = j(:,:,kd)*&
3261  &exp(-iu*(x%n_1(kd,k)-x%n_2(kd,m))*&
3262  &ang_par_f(:,:,kd))
3263  end if
3264  end do
3265 
3266  ! integrate PV_0 and KV_0
3267  if (calc_this(1)) then
3268  x_int%PV_0(1,:,:,c_tot(1)) = &
3269  &prev_mult_fac*x_int%PV_0(1,:,:,c_tot(1))
3270  x_int%KV_0(1,:,:,c_tot(1)) = &
3271  &prev_mult_fac*x_int%KV_0(1,:,:,c_tot(1))
3272  call calc_magn_int_loc(x%PV_0(:,:,:,c_loc(1))*alpha_fac,&
3273  &x_int%PV_0(1,:,:,c_tot(1)),v_int_work,step_size)
3274  call calc_magn_int_loc(x%KV_0(:,:,:,c_loc(1))*alpha_fac,&
3275  &x_int%KV_0(1,:,:,c_tot(1)),v_int_work,step_size)
3276  end if
3277 
3278  ! integrate PV_1 and KV_1
3279  if (calc_this(2)) then
3280  x_int%PV_1(1,:,:,c_tot(2)) = &
3281  &prev_mult_fac*x_int%PV_1(1,:,:,c_tot(2))
3282  x_int%KV_1(1,:,:,c_tot(2)) = &
3283  &prev_mult_fac*x_int%KV_1(1,:,:,c_tot(2))
3284  call calc_magn_int_loc(x%PV_1(:,:,:,c_loc(2))*alpha_fac,&
3285  &x_int%PV_1(1,:,:,c_tot(2)),v_int_work,step_size)
3286  call calc_magn_int_loc(x%KV_1(:,:,:,c_loc(2))*alpha_fac,&
3287  &x_int%KV_1(1,:,:,c_tot(2)),v_int_work,step_size)
3288  end if
3289 
3290  ! integrate PV_2 and KV_2
3291  if (calc_this(1)) then
3292  x_int%PV_2(1,:,:,c_tot(1)) = &
3293  &prev_mult_fac*x_int%PV_2(1,:,:,c_tot(1))
3294  x_int%KV_2(1,:,:,c_tot(1)) = &
3295  &prev_mult_fac*x_int%KV_2(1,:,:,c_tot(1))
3296  call calc_magn_int_loc(x%PV_2(:,:,:,c_loc(1))*alpha_fac,&
3297  &x_int%PV_2(1,:,:,c_tot(1)),v_int_work,step_size)
3298  call calc_magn_int_loc(x%KV_2(:,:,:,c_loc(1))*alpha_fac,&
3299  &x_int%KV_2(1,:,:,c_tot(1)),v_int_work,step_size)
3300  end if
3301  end do
3302  end do
3303 
3304  ! clean up
3305  nullify(ang_par_f)
3306  select case (x_grid_style)
3307  case (1) ! equilibrium
3308  ! do nothing
3309  case (2,3) ! solution or enriched
3310  deallocate(j)
3311  end select
3312  nullify(j)
3313 
3314  ! user output
3315  call lvl_ud(-1)
3316  contains
3317  ! Integrate local magnetic integral, adding to the previous result.
3318  !
3319  ! Makes use of nr_int_regions, int_dims, int_facs and J_exp_ang.
3321  subroutine calc_magn_int_loc(V,V_int,V_int_work,step_size)
3322  ! input / output
3323  complex(dp), intent(in) :: v(:,:,:) ! V
3324  complex(dp), intent(inout) :: v_int(:,:) ! integrated V
3325  complex(dp), intent(inout) :: v_int_work(:,:) ! workint variable
3326  real(dp) :: step_size(:,:) ! step size for every normal point
3327 
3328  ! local variables
3329  integer :: id, kd, ld ! counters
3330 
3331  ! loop over regions
3332  do ld = 1,nr_int_regions
3333  ! initialize
3334  v_int_work = 0
3335 
3336  ! integrate
3337  do kd = 1,size(step_size,2)
3338  do id = int_dims(ld,1),int_dims(ld,2),int_dims(ld,3)
3339  v_int_work(:,kd) = v_int_work(:,kd) + &
3340  &j_exp_ang(id,:,kd)*v(id,:,kd)*step_size(:,kd)
3341  end do
3342  end do
3343 
3344  ! scale by integration factor and update V_int
3345  v_int = v_int + v_int_work*int_facs(ld)
3346  end do
3347  end subroutine calc_magn_int_loc
3348  end function calc_magn_ints
3349 
3363  integer function divide_x_jobs(arr_size) result(ierr)
3364  use num_vars, only: max_x_mem, n_procs, x_jobs_lims
3365  use x_vars, only: n_mod_x
3366  use x_utilities, only: calc_memory_x
3367 
3368  character(*), parameter :: rout_name = 'divide_X_jobs'
3369 
3370  ! input / output
3371  integer, intent(in) :: arr_size
3372 
3373  ! local variables
3374  integer :: n_div ! factor by which to divide the total size
3375  real(dp) :: mem_size ! approximation of memory required for X variables
3376  integer :: n_mod_block ! nr. of modes in block
3377  integer, allocatable :: n_mod_loc(:) ! number of modes per block
3378  character(len=max_str_ln) :: block_message ! message about how many blocks
3379  character(len=max_str_ln) :: err_msg ! error message
3380 
3381  ! initialize ierr
3382  ierr = 0
3383 
3384  ! user output
3385  call writo('Dividing the perturbation jobs')
3386  call lvl_ud(1)
3387 
3388  ! calculate largest possible block of (k,m) values
3389  n_div = 0
3390  mem_size = huge(1._dp)
3391  do while (mem_size.gt.max_x_mem)
3392  n_div = n_div + 1
3393  n_mod_block = ceiling(n_mod_x*1._dp/n_div)
3394  ierr = calc_memory_x(2,arr_size,n_mod_block,mem_size)
3395  chckerr('')
3396  if (n_div.gt.n_mod_x) then
3397  ierr = 1
3398  err_msg = 'The memory limit is too low'
3399  chckerr(err_msg)
3400  end if
3401  end do
3402  if (n_div.gt.1) then
3403  block_message = 'The '//trim(i2str(n_mod_x))//&
3404  &' Fourier modes are split into '//trim(i2str(n_div))//&
3405  &' and '//trim(i2str(n_div**2))//' jobs are done separately'
3406  if (n_procs.lt.n_div**2) then
3407  block_message = trim(block_message)//', '//&
3408  &trim(i2str(n_procs))//' at a time'
3409  else
3410  block_message = trim(block_message)//', but simultaneously'
3411  end if
3412  else
3413  block_message = 'The '//trim(i2str(n_mod_x))//' Fourier modes &
3414  &can be done without splitting them'
3415  end if
3416  call writo(block_message)
3417  call writo('The total memory for all processes together is estimated &
3418  &to be about '//trim(i2str(ceiling(mem_size)))//'MB')
3419  call lvl_ud(1)
3420  call writo('(maximum: '//trim(i2str(ceiling(max_x_mem)))//&
3421  &'MB left over from equilibrium variables)')
3422  call lvl_ud(-1)
3423 
3424  ! set up jobs data as illustrated below for 3 divisions, order 1:
3425  ! [1,2,3]
3426  ! or order 2:
3427  ! [1,4,7]
3428  ! [2,5,8]
3429  ! [3,6,9]
3430  ! etc.
3431  ! Also initialize the jobs taken to false.
3432  allocate(n_mod_loc(n_div))
3433  n_mod_loc = n_mod_x/n_div ! number of radial points on this processor
3434  n_mod_loc(1:mod(n_mod_x,n_div)) = n_mod_loc(1:mod(n_mod_x,n_div)) + 1 ! add a mode to if there is a remainder
3435  call calc_x_jobs_lims(x_jobs_lims,n_mod_loc,2)
3436 
3437  ! user output
3438  call lvl_ud(-1)
3439  call writo('Perturbation jobs divided')
3440  contains
3441  ! Calculate X_jobs_lims.
3442  !
3443  ! These are stored as follows: The job index is the second one, while
3444  ! the first index ranges from 1..2. For every job in the range of job
3445  ! indices, thee first index shows the limits of
3446  ! the different orders sequentially. For example:
3447  ! - for order 1: [3,4],
3448  ! which corresponds to the 1-D range 3..4,
3449  ! - for order 2: [3,4,1,5],
3450  ! which corresponds to the 2-D range 3..4 x 1..5,
3451  !
3452  ! etc.
3454  recursive subroutine calc_x_jobs_lims(res,n_mod,ord)
3455  ! input / output
3456  integer, intent(inout), allocatable :: res(:,:) ! result
3457  integer, intent(in) :: n_mod(:) ! X jobs data
3458  integer, intent(in) :: ord ! order of data
3459 
3460  ! local variables
3461  integer, allocatable :: res_loc(:,:) ! local result
3462  integer :: n_div ! nr. of divisions of modes
3463  integer :: id ! counter
3464 
3465  ! set up nr. of divisions
3466  n_div = size(n_mod)
3467 
3468  ! (re)allocate result
3469  if (allocated(res)) deallocate(res)
3470  allocate(res(2*ord,n_div**ord))
3471 
3472  ! loop over divisions
3473  do id = 1,n_div
3474  if (ord.gt.1) then ! call lower order
3475  call calc_x_jobs_lims(res_loc,n_mod,ord-1)
3476  res(1:2*ord-2,(id-1)*n_div**(ord-1)+1:id*n_div**(ord-1)) = &
3477  &res_loc
3478  end if ! set for own order
3479  res(2*ord-1,(id-1)*n_div**(ord-1)+1:id*n_div**(ord-1)) = &
3480  &sum(n_mod(1:id-1))+1
3481  res(2*ord,(id-1)*n_div**(ord-1)+1:id*n_div**(ord-1)) = &
3482  &sum(n_mod(1:id))
3483  end do
3484  end subroutine calc_x_jobs_lims
3485  end function divide_x_jobs
3486 
3487 #if ldebug
3488 
3489  integer function print_debug_x_1(mds,grid_X,X_1) result(ierr)
3490  use num_vars, only: use_pol_flux_f
3491  use grid_utilities, only: trim_grid
3492  use grid_vars, only: alpha, n_alpha
3493  use eq_vars, only: max_flux_f
3494  use rich_vars, only: rich_lvl
3495 
3496  character(*), parameter :: rout_name = 'print_debug_X_1'
3497 
3498  ! input / output
3499  type(modes_type), intent(in) :: mds
3500  type(grid_type), intent(in) :: grid_x
3501  type(x_1_type), intent(in) :: x_1
3502 
3503  ! local variabbles
3504  type(grid_type) :: grid_trim ! trimmed grid
3505  character(len=max_str_ln), allocatable :: var_names(:) ! names of variables
3506  character(len=max_str_ln) :: file_name ! name of file
3507  integer :: k ! local row index in flux surface
3508  integer :: ld, kd, jd, rd, id ! counters
3509  integer :: min_nm_x ! minimal n (tor. flux) or m (pol. flux)
3510  integer :: n_mod_tot ! total number of modes
3511  integer :: kdl_tot(2) ! limits on total normal index for a mode
3512  integer :: rdl ! rd in local tables
3513  integer :: kd_loc ! local kd
3514  integer :: kd_loc_trim ! local kd on trimmed grid
3515  integer :: plot_dim(4) ! dimensions of plot
3516  integer :: plot_offset(4) ! local offset of plot
3517  complex(dp), allocatable :: u(:,:,:,:,:) ! U_i for all modes
3518  complex(dp), allocatable :: du(:,:,:,:,:) ! DU_i for all modes
3519  real(dp), allocatable :: x_plot(:,:,:,:) ! X of plot
3520  real(dp), allocatable :: y_plot(:,:,:,:) ! Y of plot
3521  real(dp), allocatable :: z_plot(:,:,:,:) ! Y of plot
3522 
3523  ! initiaiize ierr
3524  ierr = 0
3525 
3526  ! trim grid
3527  ierr = trim_grid(grid_x,grid_trim)
3528  chckerr('')
3529 
3530  ! set local n_mod and allocate integrated quantities
3531  min_nm_x = minval(mds%sec(:,1),1)
3532  n_mod_tot = size(mds%sec,1)
3533  allocate(u(n_mod_tot,grid_trim%n(1),grid_trim%loc_n_r,grid_trim%n(2),&
3534  &0:1))
3535  allocate(du(n_mod_tot,grid_trim%n(1),grid_trim%loc_n_r,grid_trim%n(2),&
3536  &0:1))
3537  allocate(x_plot(n_mod_tot,grid_trim%n(1),grid_trim%loc_n_r,&
3538  &grid_trim%n(2)))
3539  allocate(y_plot(n_mod_tot,grid_trim%n(1),grid_trim%loc_n_r,&
3540  &grid_trim%n(2)))
3541  allocate(z_plot(n_mod_tot,grid_trim%n(1),grid_trim%loc_n_r,&
3542  &grid_trim%n(2)))
3543  u = 0._dp
3544  du = 0._dp
3545 
3546  ! setup X, Y and Z of plot
3547  do kd = 1,grid_trim%loc_n_r
3548  z_plot(:,:,kd,1) = 1000*grid_trim%loc_r_F(kd)/max_flux_f*2*pi
3549  do rd = 1,n_mod_tot
3550  do id = 1,grid_trim%n(1)
3551  x_plot(rd,id,kd,:) = min_nm_x+rd-1
3552  do jd = 1,grid_trim%n(2)
3553  if (use_pol_flux_f) then
3554  y_plot(rd,id,kd,jd) = grid_trim%theta_F(id,jd,kd)*10
3555  else
3556  y_plot(rd,id,kd,jd) = grid_trim%zeta_F(id,jd,kd)*10
3557  end if
3558  end do
3559  end do
3560  end do
3561  end do
3562 
3563  ! select all modes combinations
3564  do k = 1,size(mds%sec,1)
3565  ! set normal limits for mode k
3566  kdl_tot(1) = mds%sec(k,2)
3567  kdl_tot(2) = mds%sec(k,3)
3568 
3569  ! limit to trimmed grid range
3570  kdl_tot(1) = max(kdl_tot(1),grid_trim%i_min)
3571  kdl_tot(2) = min(kdl_tot(2),grid_trim%i_max)
3572 
3573  ! skip if out of normal range
3574  if (kdl_tot(1).gt.kdl_tot(2)) cycle
3575 
3576  ! total mode index (starting from 1)
3577  rd = mds%sec(k,1)-min_nm_x+1
3578 
3579  ! indices in local tables
3580  rdl = mds%sec(k,4)
3581 
3582  ! loop over all normal points
3583  do kd = kdl_tot(1),kdl_tot(2)
3584  ! set indices
3585  kd_loc = kd - grid_x%i_min + 1
3586  kd_loc_trim = kd - grid_trim%i_min + 1
3587 
3588  u(rd,:,kd_loc_trim,:,0) = x_1%U_0(:,:,kd_loc,rdl)
3589  u(rd,:,kd_loc_trim,:,1) = x_1%U_1(:,:,kd_loc,rdl)
3590  du(rd,:,kd_loc_trim,:,0) = x_1%DU_0(:,:,kd_loc,rdl)
3591  du(rd,:,kd_loc_trim,:,1) = x_1%DU_1(:,:,kd_loc,rdl)
3592  end do
3593  end do
3594 
3595  ! plot
3596  allocate(var_names(n_alpha))
3597  do ld = 1,size(var_names)
3598  var_names(ld) = 'alpha = '//trim(r2strt(alpha(ld)))
3599  end do
3600  plot_dim = [n_mod_tot,grid_trim%n(1),grid_trim%n(3),n_alpha]
3601  plot_offset = [0,0,grid_trim%i_min-1,0]
3602  do id = 0,1
3603  file_name = 'U_'//trim(i2str(id))//'_R'//&
3604  &trim(i2str(rich_lvl))
3605  call plot_hdf5(var_names,'RE_'//trim(file_name),&
3606  &rp(u(:,:,:,:,id)),x=x_plot,y=y_plot,z=z_plot,&
3607  &tot_dim=plot_dim, loc_offset=plot_offset,&
3608  &col_id=4,col=1)
3609  call plot_hdf5(var_names,'IM_'//trim(file_name),&
3610  &ip(u(:,:,:,:,id)),x=x_plot,y=y_plot,z=z_plot,&
3611  &tot_dim=plot_dim, loc_offset=plot_offset,&
3612  &col_id=4,col=1)
3613 
3614  file_name = 'DU_'//trim(i2str(id))//'_R'//&
3615  &trim(i2str(rich_lvl))
3616  call plot_hdf5(var_names,'RE_'//trim(file_name),&
3617  &rp(du(:,:,:,:,id)),x=x_plot,y=y_plot,z=z_plot,&
3618  &tot_dim=plot_dim, loc_offset=plot_offset,&
3619  &col_id=4,col=1)
3620  call plot_hdf5(var_names,'IM_'//trim(file_name),&
3621  &ip(du(:,:,:,:,id)),x=x_plot,y=y_plot,z=z_plot,&
3622  &tot_dim=plot_dim, loc_offset=plot_offset,&
3623  &col_id=4,col=1)
3624  end do
3625  deallocate(var_names)
3626 
3627  ! clean up
3628  call grid_trim%dealloc()
3629  end function print_debug_x_1
3630 
3632  integer function print_debug_x_2(mds,grid_X,X_2_int) result(ierr)
3633  use num_utilities, only: c, con
3634  use x_vars, only: n_mod_x
3635  use grid_vars, only: alpha, n_alpha
3636  use grid_utilities, only: trim_grid
3637  use eq_vars, only: max_flux_f
3638  use rich_vars, only: rich_lvl
3639 
3640  character(*), parameter :: rout_name = 'print_debug_X_2'
3641 
3642  ! input / output
3643  type(modes_type), intent(in) :: mds
3644  type(grid_type), intent(in) :: grid_x
3645  type(x_2_type), intent(in) :: x_2_int
3646 
3647  ! local variables
3648  type(grid_type) :: grid_trim ! trimmed grid
3649  character(len=max_str_ln), allocatable :: var_names(:) ! names of variables
3650  character(len=max_str_ln) :: file_name ! name of file
3651  integer :: k, m ! local row and column index in flux surface
3652  integer :: id, ld, kd, rd, cd ! counters
3653  integer :: min_nm_x ! minimal n (tor. flux) or m (pol. flux)
3654  integer :: n_mod_tot ! total number of modes
3655  integer :: kdl_tot(2) ! limits on total normal index for a mode combination
3656  integer :: kd_loc ! local kd
3657  integer :: kd_loc_trim ! local kd on trimmed grid
3658  integer :: plot_dim(4) ! dimensions of plot
3659  integer :: plot_offset(4) ! local offset of plot
3660  integer :: c_loc(2) ! table index for symmetric and asymmetric quantities
3661  integer :: rdl, cdl ! rd and cd in local tables
3662  complex(dp), allocatable :: pv_int(:,:,:,:,:) ! integrated PV_i for all mode combinations
3663  complex(dp), allocatable :: kv_int(:,:,:,:,:) ! integrated KV_i for all mode combinations
3664  real(dp), allocatable :: x_plot(:,:,:,:) ! X of plot
3665  real(dp), allocatable :: y_plot(:,:,:,:) ! Y of plot
3666  real(dp), allocatable :: z_plot(:,:,:,:) ! Y of plot
3667 
3668  ! initiaiize ierr
3669  ierr = 0
3670 
3671  ! trim grid
3672  ierr = trim_grid(grid_x,grid_trim)
3673  chckerr('')
3674 
3675  ! set local n_mod and allocate integrated quantities
3676  min_nm_x = minval(mds%sec(:,1),1)
3677  n_mod_tot = size(mds%sec,1)
3678  allocate(pv_int(n_mod_tot,n_mod_tot,grid_trim%loc_n_r,&
3679  &grid_trim%n(2),0:2))
3680  allocate(kv_int(n_mod_tot,n_mod_tot,grid_trim%loc_n_r,&
3681  &grid_trim%n(2),0:2))
3682  allocate(x_plot(n_mod_tot,n_mod_tot,grid_trim%loc_n_r,1))
3683  allocate(y_plot(n_mod_tot,n_mod_tot,grid_trim%loc_n_r,1))
3684  allocate(z_plot(n_mod_tot,n_mod_tot,grid_trim%loc_n_r,1))
3685  pv_int = 0._dp
3686  kv_int = 0._dp
3687 
3688  ! setup X, Y and Z of plot
3689  do kd = 1,grid_trim%loc_n_r
3690  z_plot(:,:,kd,1) = 1000*grid_trim%loc_r_F(kd)/max_flux_f*2*pi
3691  do cd = 1,n_mod_tot
3692  do rd = 1,n_mod_tot
3693  x_plot(rd,cd,kd,1) = min_nm_x+rd-1
3694  y_plot(rd,cd,kd,1) = min_nm_x+cd-1
3695  end do
3696  end do
3697  end do
3698 
3699  ! select all mode combinations
3700  do m = 1,size(mds%sec,1)
3701  do k = 1,size(mds%sec,1)
3702  ! set normal limits for mode pair (k,m)
3703  kdl_tot(1) = max(mds%sec(k,2),mds%sec(m,2))
3704  kdl_tot(2) = min(mds%sec(k,3),mds%sec(m,3))
3705 
3706  ! limit to trimmed grid range
3707  kdl_tot(1) = max(kdl_tot(1),grid_trim%i_min)
3708  kdl_tot(2) = min(kdl_tot(2),grid_trim%i_max)
3709 
3710  ! skip if out of normal range
3711  if (kdl_tot(1).gt.kdl_tot(2)) cycle
3712 
3713  ! total mode combinations indices (starting from 1)
3714  rd = mds%sec(k,1)-min_nm_x+1
3715  cd = mds%sec(m,1)-min_nm_x+1
3716 
3717  ! indices in local tables
3718  rdl = mds%sec(k,4)
3719  cdl = mds%sec(m,4)
3720 
3721  ! index in tables
3722  c_loc(1) = c([rdl,cdl],.true.,n_mod_x)
3723  c_loc(2) = c([rdl,cdl],.false.,n_mod_x)
3724 
3725  ! loop over all normal points
3726  do kd = kdl_tot(1),kdl_tot(2)
3727  ! set indices
3728  kd_loc = kd - grid_x%i_min + 1
3729  kd_loc_trim = kd - grid_trim%i_min + 1
3730 
3731  ! symmetric quantities
3732  pv_int(rd,cd,kd_loc_trim,:,0) = &
3733  &con(x_2_int%PV_0(1,:,kd_loc,c_loc(1)),&
3734  &[rdl,cdl],.true.,[n_alpha])
3735  pv_int(rd,cd,kd_loc_trim,:,2) = &
3736  &con(x_2_int%PV_2(1,:,kd_loc,c_loc(1)),&
3737  &[rdl,cdl],.true.,[n_alpha])
3738  kv_int(rd,cd,kd_loc_trim,:,0) = &
3739  &con(x_2_int%KV_0(1,:,kd_loc,c_loc(1)),&
3740  &[rdl,cdl],.true.,[n_alpha])
3741  kv_int(rd,cd,kd_loc_trim,:,2) = &
3742  &con(x_2_int%KV_2(1,:,kd_loc,c_loc(1)),&
3743  &[rdl,cdl],.true.,[n_alpha])
3744 
3745  ! asymmetric quantities
3746  pv_int(rd,cd,kd_loc_trim,:,1) = &
3747  &x_2_int%PV_1(1,:,kd_loc,c_loc(2))
3748  kv_int(rd,cd,kd_loc_trim,:,1) = &
3749  &x_2_int%KV_1(1,:,kd_loc,c_loc(2))
3750  end do
3751  end do
3752  end do
3753 
3754  ! plot
3755  allocate(var_names(n_alpha))
3756  do ld = 1,size(var_names)
3757  var_names(ld) = 'alpha = '//trim(r2strt(alpha(ld)))
3758  end do
3759  plot_dim = [n_mod_tot,n_mod_tot,grid_trim%n(3),n_alpha]
3760  plot_offset = [0,0,grid_trim%i_min-1,0]
3761  do id = 0,2
3762  file_name = 'PV_'//trim(i2str(id))//'_int_R'//&
3763  &trim(i2str(rich_lvl))
3764  call plot_hdf5(var_names,'RE_'//trim(file_name),&
3765  &rp(pv_int(:,:,:,:,id)),x=x_plot,y=y_plot,z=z_plot,&
3766  &tot_dim=plot_dim, loc_offset=plot_offset,&
3767  &col_id=4,col=1)
3768  call plot_hdf5(var_names,'IM_'//trim(file_name),&
3769  &ip(pv_int(:,:,:,:,id)),x=x_plot,y=y_plot,z=z_plot,&
3770  &tot_dim=plot_dim, loc_offset=plot_offset,&
3771  &col_id=4,col=1)
3772 
3773  file_name = 'KV_'//trim(i2str(id))//'_int_R'//&
3774  &trim(i2str(rich_lvl))
3775  call plot_hdf5(var_names,'RE_'//trim(file_name),&
3776  &rp(kv_int(:,:,:,:,id)),x=x_plot,y=y_plot,z=z_plot,&
3777  &tot_dim=plot_dim, loc_offset=plot_offset,&
3778  &col_id=4,col=1)
3779  call plot_hdf5(var_names,'IM_'//trim(file_name),&
3780  &ip(kv_int(:,:,:,:,id)),x=x_plot,y=y_plot,z=z_plot,&
3781  &tot_dim=plot_dim, loc_offset=plot_offset,&
3782  &col_id=4,col=1)
3783  end do
3784  deallocate(var_names)
3785 
3786  ! clean up
3787  call grid_trim%dealloc()
3788  end function print_debug_x_2
3789 #endif
3790 end module
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
num_vars::x_grid_style
integer, public x_grid_style
style for normal component of X grid (1: eq, 2: sol, 3: enriched)
Definition: num_vars.f90:99
x_ops::calc_pv
integer function calc_pv(grid_eq, grid_X, eq_1, eq_2, X_a, X_b, X, lim_sec_X)
calculate (pol. flux) or (tor. flux) at all equilibrium loc_n_r values.
Definition: X_ops.f90:2646
num_vars::u_style
integer, public u_style
style for calculation of U (1: ord.2, 2: ord.1, 1: ord.0)
Definition: num_vars.f90:91
num_vars::max_name_ln
integer, parameter, public max_name_ln
maximum length of filenames
Definition: num_vars.f90:51
hdf5_vars::dealloc_var_1d
Deallocates 1D variable.
Definition: HDF5_vars.f90:68
x_utilities::get_sec_x_range
subroutine, public get_sec_x_range(sec_X_range_loc, sec_X_range_tot, m, sym, lim_sec_X)
Gets one of the the local ranges of contiguous tensorial perturbation variables to be printed or read...
Definition: X_utilities.f90:128
x_ops
Operations considering perturbation quantities.
Definition: X_ops.f90:4
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_utilities::con
Either takes the complex conjugate of a square matrix element A, defined on a 3-D grid,...
Definition: num_utilities.f90:221
eq_vars::vac_perm
real(dp), public vac_perm
either usual mu_0 (default) or normalized
Definition: eq_vars.f90:48
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
input_utilities::get_log
logical function, public get_log(yes, ind)
Queries for a logical value yes or no, where the default answer is also to be provided.
Definition: input_utilities.f90:22
num_vars::max_str_ln
integer, parameter, public max_str_ln
maximum length of strings
Definition: num_vars.f90:50
x_ops::calc_res_surf
integer function, public calc_res_surf(mds, grid_eq, eq, res_surf, info, jq)
Calculates resonating flux surfaces for the perturbation modes.
Definition: X_ops.f90:1638
x_vars::min_m_x
integer, dimension(:), allocatable, public min_m_x
lowest poloidal mode number m_X, in total eq grid
Definition: X_vars.f90:133
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
hdf5_ops
Operations on HDF5 and XDMF variables.
Definition: HDF5_ops.f90:27
grid_vars::max_par_x
real(dp), public max_par_x
max. of parallel coordinate [ ] in field-aligned grid
Definition: grid_vars.f90:25
x_ops::print_debug_x_2
integer function, public print_debug_x_2(mds, grid_X, X_2_int)
Prints debug information for X_2 driver.
Definition: X_ops.f90:3633
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
x_ops::init_modes
integer function, public init_modes(grid_eq, eq)
Initializes some variables concerning the mode numbers.
Definition: X_ops.f90:919
num_vars::eq_job_nr
integer, public eq_job_nr
nr. of eq job
Definition: num_vars.f90:79
x_utilities::is_necessary_x
logical function, public is_necessary_x(sym, sec_X_id, lim_sec_X)
Determines whether a variable needs to be considered:
Definition: X_utilities.f90:197
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_vars::n_procs
integer, public n_procs
nr. of MPI processes
Definition: num_vars.f90:69
hdf5_vars
Variables pertaining to HDF5 and XDMF.
Definition: HDF5_vars.f90:4
output_ops::print_ex_2d
Print 2-D output on a file.
Definition: output_ops.f90:47
num_vars::x_style
integer, public x_style
style for secondary mode numbers (1: prescribed, 2: fast)
Definition: num_vars.f90:95
str_utilities
Operations on strings.
Definition: str_utilities.f90:4
x_vars::min_sec_x
integer, public min_sec_x
m_X (pol. flux) or n_X (tor. flux) (only for X style 1)
Definition: X_vars.f90:127
grid_vars::grid_type
Type for grids.
Definition: grid_vars.f90:59
x_vars::set_nn_mod
integer function, public set_nn_mod(sym, lim_sec_X)
Sets number of entries for tensorial perturbation variables.
Definition: X_vars.f90:383
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::eq_1_type
flux equilibrium type
Definition: eq_vars.f90:63
num_vars::ltest
logical, public ltest
whether or not to call the testing routines
Definition: num_vars.f90:112
test_du
integer function test_du()
Definition: X_ops.f90:2515
num_vars::min_zeta_plot
real(dp), public min_zeta_plot
min. of zeta_plot
Definition: num_vars.f90:161
grid_vars::alpha
real(dp), dimension(:), allocatable, public alpha
field line label alpha
Definition: grid_vars.f90:28
hdf5_vars::max_dim_var_1d
integer, parameter, public max_dim_var_1d
maximum dimension of var_1D
Definition: HDF5_vars.f90:21
num_vars::min_theta_plot
real(dp), public min_theta_plot
min. of theta_plot
Definition: num_vars.f90:159
grid_vars::min_par_x
real(dp), public min_par_x
min. of parallel coordinate [ ] in field-aligned grid
Definition: grid_vars.f90:24
x_vars::prim_x
integer, public prim_x
n_X (pol. flux) or m_X (tor. flux)
Definition: X_vars.f90:126
x_ops::debug_check_x_modes_2
logical, public debug_check_x_modes_2
plot debug information for check_x_modes_2()
Definition: X_ops.f90:26
x_ops::calc_kv
integer function calc_kv(grid_eq, grid_X, eq_1, eq_2, X_a, X_b, X, lim_sec_X)
calculate (pol. flux) or (tor. flux) at all equilibrium loc_n_r values.
Definition: X_ops.f90:2879
hdf5_vars::var_1d_type
1D equivalent of multidimensional variables, used for internal HDF5 storage.
Definition: HDF5_vars.f90:48
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
x_vars::max_n_x
integer, dimension(:), allocatable, public max_n_x
highest poloidal mode number m_X, in total eq grid
Definition: X_vars.f90:132
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
x_vars::x_2_type
tensorial perturbation type
Definition: X_vars.f90:81
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
hdf5_ops::print_hdf5_arrs
integer function, public print_hdf5_arrs(vars, PB3D_name, head_name, rich_lvl, disp_info, ind_print, remove_previous_arrs)
Prints a series of arrays, in the form of an array of pointers, to an HDF5 file.
Definition: HDF5_ops.f90:1132
x_ops::setup_modes
integer function, public setup_modes(mds, grid_eq, grid, plot_name)
Sets up some variables concerning the mode numbers.
Definition: X_ops.f90:1060
x_vars
Variables pertaining to the perturbation quantities.
Definition: X_vars.f90:4
x_ops::check_x_modes
integer function, public check_x_modes(grid_eq, eq)
Checks whether the high-n approximation is valid:
Definition: X_ops.f90:1350
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
x_utilities
Numerical utilities related to perturbation operations.
Definition: X_utilities.f90:4
x_vars::modes_type
mode number type
Definition: X_vars.f90:36
x_ops::print_output_x
Print either vectorial or tensorial perturbation quantities of a certain order to an output file.
Definition: X_ops.f90:85
x_ops::divide_x_jobs
integer function, public divide_x_jobs(arr_size)
Divides the perturbation jobs.
Definition: X_ops.f90:3364
grid_utilities::calc_eqd_grid
Calculate grid of equidistant points, where optionally the last point can be excluded.
Definition: grid_utilities.f90:75
num_vars::magn_int_style
integer, public magn_int_style
style for magnetic integrals (1: trapezoidal, 2: Simpson 3/8)
Definition: num_vars.f90:124
num_utilities
Numerical utilities.
Definition: num_utilities.f90:4
num_vars::alpha_style
integer, public alpha_style
style for alpha (1: one field line, many turns, 2: many field lines, one turn)
Definition: num_vars.f90:100
str_utilities::r2str
elemental character(len=max_str_ln) function, public r2str(k)
Convert a real (double) to string.
Definition: str_utilities.f90:42
num_vars::k_style
integer, public k_style
style for kinetic energy
Definition: num_vars.f90:93
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
x_ops::resonance_plot
integer function, public resonance_plot(mds, grid_eq, eq)
plot -profile or -profile in flux coordinates with or indicated if requested.
Definition: X_ops.f90:1814
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
x_ops::calc_magn_ints
integer function, public calc_magn_ints(grid_eq, grid_X, eq, X, X_int, prev_style, lim_sec_X)
Calculate the magnetic integrals from and in an equidistant grid where the step size can vary depen...
Definition: X_ops.f90:3081
num_vars::pb3d_name
character(len=max_str_ln), public pb3d_name
name of PB3D output file
Definition: num_vars.f90:139
x_vars::max_m_x
integer, dimension(:), allocatable, public max_m_x
highest poloidal mode number m_X, in total eq grid
Definition: X_vars.f90:134
x_vars::min_nm_x
integer, public min_nm_x
minimum for the high-n theory (debable)
Definition: X_vars.f90:130
grid_vars
Variables pertaining to the different grids used.
Definition: grid_vars.f90:4
x_vars::max_sec_x
integer, public max_sec_x
m_X (pol. flux) or n_X (tor. flux) (only for\ c X style 1)
Definition: X_vars.f90:128
messages::lvl_ud
subroutine, public lvl_ud(inc)
Increases/decreases lvl of output.
Definition: messages.f90:254
input_utilities::pause_prog
subroutine, public pause_prog(ind)
Pauses the running of the program.
Definition: input_utilities.f90:226
x_ops::print_debug_x_1
integer function, public print_debug_x_1(mds, grid_X, X_1)
Prints debug information for X_1 driver.
Definition: X_ops.f90:3490
num_vars::no_plots
logical, public no_plots
no plots made
Definition: num_vars.f90:140
x_ops::calc_u
integer function calc_u(grid_eq, grid_X, eq_1, eq_2, X)
Calculate , or , .
Definition: X_ops.f90:2097
x_vars::n_mod_x
integer, public n_mod_x
size of m_X (pol. flux) or n_X (tor. flux)
Definition: X_vars.f90:129
x_vars::min_n_x
integer, dimension(:), allocatable, public min_n_x
lowest poloidal mode number m_X, in total eq grid
Definition: X_vars.f90:131
num_ops::calc_zero_zhang
character(len=max_str_ln) function, public calc_zero_zhang(zero, fun, x_int_in)
Finds the zero of a function using Zhang's method, which is simpler than Brent's method.
Definition: num_ops.f90:462
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
x_ops::calc_x
Calculates either vectorial or tensorial perturbation variables.
Definition: X_ops.f90:39
x_utilities::calc_memory_x
integer function, public calc_memory_x(ord, arr_size, n_mod, mem_size)
Calculate memory in MB necessary for X variables.
Definition: X_utilities.f90:236
rich_vars::rich_lvl
integer, public rich_lvl
current level of Richardson extrapolation
Definition: rich_vars.f90:19
x_vars::x_1_type
vectorial perturbation type
Definition: X_vars.f90:51
num_vars::eq_jobs_lims
integer, dimension(:,:), allocatable, public eq_jobs_lims
data about eq jobs: [ , ] for all jobs
Definition: num_vars.f90:77
x_ops::redistribute_output_x
Redistribute the perturbation variables.
Definition: X_ops.f90:51
grid_vars::n_alpha
integer, public n_alpha
nr. of field-lines
Definition: grid_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
num_vars::rank
integer, public rank
MPI rank.
Definition: num_vars.f90:68
input_utilities
Numerical utilities related to input.
Definition: input_utilities.f90:4
num_vars::tol_norm
real(dp), public tol_norm
tolerance for normal range (normalized to 0..1)
Definition: num_vars.f90:134
num_vars::max_x_mem
real(dp), public max_x_mem
maximum memory for perturbation calculations for all processes [MB]
Definition: num_vars.f90:75
eq_vars::eq_2_type
metric equilibrium type
Definition: eq_vars.f90:114
num_vars::n_theta_plot
integer, public n_theta_plot
nr. of poloidal points in plot
Definition: num_vars.f90:156
num_vars::x_jobs_lims
integer, dimension(:,:), allocatable, public x_jobs_lims
data about X jobs: [ , , , ] for all jobs
Definition: num_vars.f90:76
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
mpi_utilities::redistribute_var
integer function, public redistribute_var(var, dis_var, lims, lims_dis)
Redistribute variables according to new limits.
Definition: MPI_utilities.f90:330
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
num_vars::norm_disc_prec_x
integer, public norm_disc_prec_x
precision for normal discretization for perturbation
Definition: num_vars.f90:121
this documentation was created with Doxygen 1.8.17 by Toon Weyens