Report Typos and Errors    
Semi-Lagrangian Library
Modular library for kinetic and gyrokinetic simulations of plasmas in fusion energy devices.
sll_m_gyroaverage_2d_polar.F90
Go to the documentation of this file.
1 !**************************************************************
2 ! Copyright INRIA
3 ! Authors :
4 ! CALVI project team
5 !
6 ! This code SeLaLib (for Semi-Lagrangian-Library)
7 ! is a parallel library for simulating the plasma turbulence
8 ! in a tokamak.
9 !
10 ! This software is governed by the CeCILL-B license
11 ! under French law and abiding by the rules of distribution
12 ! of free software. You can use, modify and redistribute
13 ! the software under the terms of the CeCILL-B license as
14 ! circulated by CEA, CNRS and INRIA at the following URL
15 ! "http://www.cecill.info".
16 !**************************************************************
17 
18 module sll_m_gyroaverage_2d_polar
19 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
20 #include "sll_memory.h"
21 #include "sll_working_precision.h"
22 
23 ! use F77_fftpack, only: &
24 ! dfftb, &
25 ! dfftf, &
26 ! dffti
27 
28  use sll_m_constants, only: &
29  sll_p_pi
30 
31  use sll_m_gyroaverage_utilities, only: &
32  sll_s_compute_shape_circle
33 
34  implicit none
35 
36  public :: &
37  sll_s_compute_gyroaverage_pade_high_order_polar, &
38  sll_s_compute_gyroaverage_pade_polar, &
39  sll_s_compute_gyroaverage_points_polar_hermite, &
40  sll_s_compute_gyroaverage_points_polar_hermite_c1, &
41  sll_s_compute_gyroaverage_points_polar_spl, &
42  sll_s_compute_gyroaverage_points_polar_with_invar_hermite_c1, &
43  sll_s_compute_gyroaverage_points_polar_with_invar_spl, &
44  sll_s_compute_gyroaverage_pre_compute_polar_hermite_c1, &
45  sll_s_compute_gyroaverage_pre_compute_polar_spl, &
46  sll_s_compute_gyroaverage_pre_compute_polar_spl_fft, &
47  sll_f_new_plan_gyroaverage_polar_hermite, &
48  sll_f_new_plan_gyroaverage_polar_pade, &
49  sll_f_new_plan_gyroaverage_polar_splines, &
50  sll_s_penta, &
51  sll_s_pre_compute_gyroaverage_polar_hermite_c1, &
52  sll_s_pre_compute_gyroaverage_polar_spl, &
53  sll_s_pre_compute_gyroaverage_polar_spl_fft, &
54  sll_t_plan_gyroaverage_polar
55 
56  private
57 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
58 
59  type sll_t_plan_gyroaverage_polar
60 
61  sll_real64 :: eta_min(2)
62  sll_real64 :: eta_max(2)
63  sll_int32 :: nc(2)
64 
65  sll_int32 :: n_points
66  sll_int32 :: interp_degree(2)
67 
68  sll_real64, dimension(:, :), pointer :: points
69  sll_real64, dimension(:, :, :), pointer :: deriv
70  sll_int32, dimension(:), pointer :: pre_compute_n
71  sll_real64, dimension(:, :), pointer :: pre_compute_coeff
72  sll_int32, dimension(:, :), pointer :: pre_compute_index
73  sll_real64, dimension(:), pointer :: pre_compute_coeff_spl
74  sll_real64, dimension(:, :), pointer :: lunat
75  sll_real64, dimension(:), pointer :: luper
76  sll_real64, dimension(:, :, :), pointer :: a_fft
77 
78  end type sll_t_plan_gyroaverage_polar
79 
80 contains
81 
82  function sll_f_new_plan_gyroaverage_polar_hermite(eta_min, eta_max, Nc, N_points, interp_degree, deriv_size) result(this)
83 
84  implicit none
85 
86  sll_real64, intent(in) :: eta_min(2)
87  sll_real64, intent(in) :: eta_max(2)
88  sll_int32, intent(in) :: nc(2)
89  sll_int32, intent(in) :: n_points
90  sll_int32, intent(in) :: interp_degree(2)
91  sll_int32, intent(in) :: deriv_size
92  type(sll_t_plan_gyroaverage_polar), pointer :: this
93 
94  sll_int32 :: err
95 
96  sll_allocate(this, err)
97  sll_allocate(this%deriv(deriv_size, nc(1) + 1, nc(2) + 1), err)
98  sll_allocate(this%points(3, n_points), err)
99 
100  call sll_s_compute_shape_circle(this%points, n_points)
101 
102  this%eta_min = eta_min
103  this%eta_max = eta_max
104  this%Nc = nc
105  this%N_points = n_points
106  this%interp_degree = interp_degree
107 
108  end function sll_f_new_plan_gyroaverage_polar_hermite
109 
110  function sll_f_new_plan_gyroaverage_polar_splines(eta_min, eta_max, Nc, N_points) result(this)
111 
112  implicit none
113 
114  sll_real64, intent(in) :: eta_min(2)
115  sll_real64, intent(in) :: eta_max(2)
116  sll_int32, intent(in) :: nc(2)
117  sll_int32, intent(in) :: n_points
118  type(sll_t_plan_gyroaverage_polar), pointer :: this
119 
120  sll_int32 :: err
121 
122  sll_allocate(this, err)
123  sll_allocate(this%points(3, n_points), err)
124 
125  call sll_s_compute_shape_circle(this%points, n_points)
126 
127  this%eta_min = eta_min
128  this%eta_max = eta_max
129  this%Nc = nc
130  this%N_points = n_points
131 
132  end function sll_f_new_plan_gyroaverage_polar_splines
133 
134  function sll_f_new_plan_gyroaverage_polar_pade(eta_min, eta_max, Nc) result(this)
135 
136  implicit none
137 
138  sll_real64, intent(in) :: eta_min(2)
139  sll_real64, intent(in) :: eta_max(2)
140  sll_int32, intent(in) :: nc(2)
141  type(sll_t_plan_gyroaverage_polar), pointer :: this
142  sll_int32 :: err
143 
144  sll_allocate(this, err)
145 
146  this%eta_min = eta_min
147  this%eta_max = eta_max
148  this%Nc = nc
149 
150  end function sll_f_new_plan_gyroaverage_polar_pade
151 
152  subroutine sll_s_compute_gyroaverage_points_polar_hermite(gyro, f, rho)
153  type(sll_t_plan_gyroaverage_polar) :: gyro
154  sll_real64, dimension(:, :), intent(inout) :: f
155  sll_real64, intent(in) :: rho
156  sll_int32 :: i, j, k, ii(2)
157  sll_real64::fval, sum_fval, eta_star(2), eta(2), delta_eta(2), x(2)
158 
159  fval = 0._f64
160  delta_eta(1) = (gyro%eta_max(1) - gyro%eta_min(1))/real(gyro%Nc(1), f64)
161  delta_eta(2) = (gyro%eta_max(2) - gyro%eta_min(2))/real(gyro%Nc(2), f64)
162 
163  call hermite_coef_nat_per(f(1:gyro%Nc(1) + 1, 1:gyro%Nc(2)), gyro%deriv, gyro%Nc, gyro%interp_degree)
164 
165  do j = 1, gyro%Nc(2)
166  eta(2) = gyro%eta_min(2) + real(j - 1, f64)*delta_eta(2)
167  do i = 1, gyro%Nc(1) + 1
168  eta(1) = gyro%eta_min(1) + real(i - 1, f64)*delta_eta(1)
169  sum_fval = 0._f64
170  do k = 1, gyro%N_points
171  x(1) = eta(1)*cos(eta(2)) + rho*gyro%points(1, k)
172  x(2) = eta(1)*sin(eta(2)) + rho*gyro%points(2, k)
173  call localize_polar(x, gyro%eta_min, gyro%eta_max, ii, eta_star, gyro%Nc)
174  call interpolate_hermite(gyro%deriv, ii, eta_star, fval, gyro%Nc)
175  sum_fval = sum_fval + gyro%points(3, k)*fval
176  end do
177  f(i, j) = sum_fval
178  end do
179  end do
180 
181  f(1:gyro%Nc(1) + 1, gyro%Nc(2) + 1) = f(1:gyro%Nc(1) + 1, 1)
182 
183  ! We force zero condition at radius r_min and r_max
184  f(1, 1:gyro%Nc(2) + 1) = 0._f64
185  f(gyro%Nc(1) + 1, 1:gyro%Nc(2) + 1) = 0._f64
186 
187  end subroutine sll_s_compute_gyroaverage_points_polar_hermite
188 
189  subroutine sll_s_compute_gyroaverage_points_polar_hermite_c1(gyro, f, rho)
190  type(sll_t_plan_gyroaverage_polar) :: gyro
191  sll_real64, dimension(:, :), intent(inout) :: f
192  sll_real64, intent(in)::rho
193  sll_int32 ::i, j, k, ii(2)
194  sll_real64::fval, sum_fval, eta_star(2), eta(2), delta_eta(2), x(2)
195 
196  fval = 0._f64
197  delta_eta(1) = (gyro%eta_max(1) - gyro%eta_min(1))/real(gyro%Nc(1), f64)
198  delta_eta(2) = (gyro%eta_max(2) - gyro%eta_min(2))/real(gyro%Nc(2), f64)
199 
200  call hermite_c1_coef_nat_per(f(1:gyro%Nc(1) + 1, 1:gyro%Nc(2)), gyro%deriv, gyro%Nc, gyro%interp_degree)
201 
202  do j = 1, gyro%Nc(2)
203  eta(2) = gyro%eta_min(2) + real(j - 1, f64)*delta_eta(2)
204  !eta(2)=gyro%eta_min(2) ! to uncomment for compatibility with precompute
205  do i = 1, gyro%Nc(1) + 1
206  eta(1) = gyro%eta_min(1) + real(i - 1, f64)*delta_eta(1)
207  sum_fval = 0._f64
208  do k = 1, gyro%N_points
209  x(1) = eta(1)*cos(eta(2)) + rho*gyro%points(1, k)
210  x(2) = eta(1)*sin(eta(2)) + rho*gyro%points(2, k)
211  eta_star(1) = sqrt(x(1)**2 + x(2)**2)
212  call localize_nat(ii(1), eta_star(1), gyro%eta_min(1), gyro%eta_max(1), gyro%Nc(1))
213  eta_star(2) = modulo(atan2(x(2), x(1)), 2._f64*sll_p_pi)
214  !eta_star(2)=modulo(atan2(x(2),x(1))+real(j-1,f64)*delta_eta(2),2._f64*M_PI) &
215  ! to uncomment for compatibility with precompute
216  call localize_per(ii(2), eta_star(2), gyro%eta_min(2), gyro%eta_max(2), gyro%Nc(2))
217  call interpolate_hermite_c1(gyro%deriv, ii, eta_star, fval, gyro%Nc)
218  sum_fval = sum_fval + gyro%points(3, k)*fval
219  end do
220  f(i, j) = sum_fval
221  end do
222  end do
223 
224  f(1:gyro%Nc(1) + 1, gyro%Nc(2) + 1) = f(1:gyro%Nc(1) + 1, 1)
225 
226  ! We force zero condition at radius r_min and r_max
227  f(1, 1:gyro%Nc(2) + 1) = 0._f64
228  f(gyro%Nc(1) + 1, 1:gyro%Nc(2) + 1) = 0._f64
229 
230  end subroutine sll_s_compute_gyroaverage_points_polar_hermite_c1
231 
232  subroutine sll_s_compute_gyroaverage_points_polar_with_invar_hermite_c1(gyro, f, rho)
233  type(sll_t_plan_gyroaverage_polar) :: gyro
234  sll_real64, dimension(:, :), intent(inout) :: f
235  sll_real64, intent(in)::rho
236  sll_int32 ::i, j, k, ii(2)
237  sll_real64::fval, eta_star(2), eta(2), delta_eta(2), x(2), angle
238  sll_real64, dimension(:, :), allocatable::sum_fval
239  sll_int32 ::error
240 
241  sll_allocate(sum_fval(0:gyro%Nc(1), 0:gyro%Nc(2) - 1), error)
242 
243  sum_fval = 0._f64
244  delta_eta(1) = (gyro%eta_max(1) - gyro%eta_min(1))/real(gyro%Nc(1), f64)
245  delta_eta(2) = (gyro%eta_max(2) - gyro%eta_min(2))/real(gyro%Nc(2), f64)
246 
247  call hermite_c1_coef_nat_per(f(1:gyro%Nc(1) + 1, 1:gyro%Nc(2)), gyro%deriv, gyro%Nc, gyro%interp_degree)
248 
249  do i = 0, gyro%Nc(1)
250  eta(1) = gyro%eta_min(1) + real(i, f64)*delta_eta(1) ! rayon du centre
251  eta(2) = gyro%eta_min(2) ! angle quelconque pour le centre
252  do k = 1, gyro%N_points
253  x(1) = eta(1)*cos(eta(2)) + rho*gyro%points(1, k)
254  x(2) = eta(1)*sin(eta(2)) + rho*gyro%points(2, k)
255  eta_star(1) = sqrt(x(1)**2 + x(2)**2)
256  call localize_nat(ii(1), eta_star(1), gyro%eta_min(1), gyro%eta_max(1), gyro%Nc(1))
257  angle = atan2(x(2), x(1))
258  do j = 0, gyro%Nc(2) - 1
259  eta_star(2) = modulo(angle + real(j, f64)*delta_eta(2), 2._f64*sll_p_pi)
260  call localize_per(ii(2), eta_star(2), gyro%eta_min(2), gyro%eta_max(2), gyro%Nc(2))
261  call interpolate_hermite_c1(gyro%deriv, ii, eta_star, fval, gyro%Nc)
262  sum_fval(i, j) = sum_fval(i, j) + gyro%points(3, k)*fval
263  end do
264  end do
265  end do
266 
267  f(1:gyro%Nc(1) + 1, 1:gyro%Nc(2)) = sum_fval(0:gyro%Nc(1), 0:gyro%Nc(2) - 1)
268  f(1:gyro%Nc(1) + 1, gyro%Nc(2) + 1) = f(1:gyro%Nc(1) + 1, 1)
269 
270  ! We force zero condition at radius r_min and r_max
271  f(1, 1:gyro%Nc(2) + 1) = 0._f64
272  f(gyro%Nc(1) + 1, 1:gyro%Nc(2) + 1) = 0._f64
273 
274  sll_deallocate_array(sum_fval, error)
275 
276  end subroutine sll_s_compute_gyroaverage_points_polar_with_invar_hermite_c1
277 
278  subroutine sll_s_pre_compute_gyroaverage_polar_hermite_c1(gyro, rho)
279  type(sll_t_plan_gyroaverage_polar) :: gyro
280  sll_real64, intent(in)::rho
281  sll_int32, dimension(:, :), allocatable :: buf
282  sll_int32 ::i, j, k, ell_1, ell_2, ii(2), s, nb, ind(2)
283  sll_real64::val(4, 0:1, 0:1), eta_star(2), eta(2), delta_eta(2), x(2)
284  sll_int32 ::error
285 
286  val = 0._f64
287  delta_eta(1) = (gyro%eta_max(1) - gyro%eta_min(1))/real(gyro%Nc(1), f64)
288  delta_eta(2) = (gyro%eta_max(2) - gyro%eta_min(2))/real(gyro%Nc(2), f64)
289 
290  sll_allocate(buf(0:gyro%Nc(1), 0:gyro%Nc(2)), error)
291 
292  sll_allocate(gyro%pre_compute_N(gyro%Nc(1) + 1), error)
293 
294  !gyro%pre_compute_N(i)
295  !gyro%pre_compute_index(1:2,s)
296  !gyro%deriv(ell,ii(1),ii(2))*gyro_pre_compute_coeff(s)
297 
298  eta(2) = gyro%eta_min(2)
299  buf = 0
300  nb = 0
301  do i = 1, gyro%Nc(1) + 1
302  eta(1) = gyro%eta_min(1) + real(i - 1, f64)*delta_eta(1)
303  s = 0
304  do k = 1, gyro%N_points
305  x(1) = eta(1)*cos(eta(2)) + rho*gyro%points(1, k)
306  x(2) = eta(1)*sin(eta(2)) + rho*gyro%points(2, k)
307  call localize_polar(x, gyro%eta_min, gyro%eta_max, ii, eta_star, gyro%Nc)
308  do ell_2 = 0, 1
309  ind(2) = ii(2) + ell_2
310  do ell_1 = 0, 1
311  ind(1) = ii(1) + ell_1
312  if (buf(ind(1), ind(2)) .ne. i) then
313  s = s + 1
314  buf(ind(1), ind(2)) = i
315  end if
316  end do
317  end do
318  end do
319  gyro%pre_compute_N(i) = s
320  nb = nb + s
321  end do
322 
323  print *, '#N_points pre_compute=', nb, gyro%Nc(1), nb/gyro%Nc(1)
324 
325  sll_allocate(gyro%pre_compute_index(1:2, nb), error)
326  sll_allocate(gyro%pre_compute_coeff(4, nb), error)
327 
328  buf = 0
329  nb = 0
330  s = 0
331  val = 0._f64
332  eta(2) = gyro%eta_min(2)
333  do i = 1, gyro%Nc(1) + 1
334  eta(1) = gyro%eta_min(1) + real(i - 1, f64)*delta_eta(1)
335  do k = 1, gyro%N_points
336  x(1) = eta(1)*cos(eta(2)) + rho*gyro%points(1, k)
337  x(2) = eta(1)*sin(eta(2)) + rho*gyro%points(2, k)
338  call localize_polar(x, gyro%eta_min, gyro%eta_max, ii, eta_star, gyro%Nc)
339  call contribution_hermite_c1(eta_star, val)
340  !print *,eta_star
341  !print *,i,k,val(1,0,0),val(1,1,0),ii(1),eta_star(1)
342 
343  val = val*gyro%points(3, k)
344  do ell_2 = 0, 1
345  ind(2) = ii(2) + ell_2
346  do ell_1 = 0, 1
347  ind(1) = ii(1) + ell_1
348  j = buf(ind(1), ind(2))
349  if (j <= nb) then
350  s = s + 1
351  buf(ind(1), ind(2)) = s
352  gyro%pre_compute_coeff(1:4, s) = val(1:4, ell_1, ell_2)
353  gyro%pre_compute_index(1, s) = ind(1)
354  gyro%pre_compute_index(2, s) = ind(2)
355  else
356  gyro%pre_compute_coeff(1:4, j) = gyro%pre_compute_coeff(1:4, j) + val(1:4, ell_1, ell_2)
357  end if
358  end do
359  end do
360  end do
361  nb = s
362  end do
363 
364  !print *,gyro%pre_compute_coeff
365  !stop
366 
367  print *, '#min/max=', minval(gyro%pre_compute_coeff(1:4, :)), maxval(gyro%pre_compute_coeff(1:4, :))
368  sll_deallocate_array(buf, error)
369  end subroutine sll_s_pre_compute_gyroaverage_polar_hermite_c1
370 
371 ! example of CSR format see
372 ! http://netlib.org/linalg/html_templates/node91.html#SECTION00931100000000000000
373 !A= 10 0 0 0 -2 0
374 ! 3 9 0 0 0 3
375 ! 0 7 8 7 0 0
376 ! 3 0 8 7 5 0
377 ! 0 8 0 9 9 13
378 !
379 !val = 10 -2 3 9 3 7 ...
380 !col_ind = 1 5 1 2 6 2 ...
381 !
382 !row_ptr = 1 3 6 9 13 17 20
383 !
384 ! here ia stands for row_prt
385 ! ja for col_ind
386 ! and a for val
387 
388  subroutine assemble_csr_from_pre_compute( &
389  pre_compute_coeff, &
390  pre_compute_N, &
391  pre_compute_index, &
392  nb, &
393  Nc, &
394  ia, &
395  ja, &
396  a, &
397  num_rows, &
398  num_nz)
399  sll_real64, dimension(:, :), intent(in) :: pre_compute_coeff
400  sll_int32, dimension(:), intent(in) :: pre_compute_n
401  sll_int32, dimension(:, :), intent(in) :: pre_compute_index
402  sll_int32, intent(in) :: nb
403  sll_int32, intent(in) :: nc(2)
404  sll_int32, dimension(:), pointer :: ia
405  sll_int32, dimension(:), pointer :: ja
406  sll_real64, dimension(:), pointer :: a
407  sll_int32, intent(out) :: num_rows
408  sll_int32, intent(out) :: num_nz
409  sll_int32 :: ierr
410  sll_int32 :: i
411 
412  !SLL_ALLOCATE(gyro%pre_compute_N(gyro%Nc(1)+1),error)
413  !SLL_ALLOCATE(gyro%pre_compute_index(1:2,nb),error)
414  !SLL_ALLOCATE(gyro%pre_compute_coeff(4,nb),error)
415 
416  if ((size(pre_compute_coeff, 1) < 4) .or. (size(pre_compute_coeff, 2) < nb)) then
417  print *, '#bad size for pre_compute_coeff'
418  print *, '#in assemble_csr_from_pre_comput'
419  stop
420  end if
421 
422  if (size(pre_compute_n, 1) < nc(1) + 1) then
423  print *, '#bad size for pre_compute_N'
424  print *, '#in assemble_csr_from_pre_comput'
425  stop
426  end if
427 
428  if ((size(pre_compute_index, 1) < 2) .or. (size(pre_compute_index, 2) < nb)) then
429  print *, '#bad size for pre_compute_index'
430  print *, '#in assemble_csr_from_pre_comput'
431  stop
432  end if
433 
434  num_rows = (nc(1) + 1)*nc(2)
435 
436  num_nz = nc(2)
437  do i = 1, nc(1) + 1
438  num_nz = num_nz*pre_compute_n(i)
439  end do
440 
441  sll_allocate(ia(num_rows + 1), ierr)
442  sll_allocate(ja(num_nz), ierr)
443  sll_allocate(a(num_nz), ierr)
444 
445  ! inum_rows = 0
446 ! s = 0
447 ! do j=1,Nc(2)
448 ! do i=1,Nc(1)+1
449 ! inum_rows = inum_rows+1
450 ! do k=1,gyro%pre_compute_N(i)
451 ! s = s+1
452 ! do ell=1,4
453 ! enddo
454 ! enddo
455 ! enddo
456 ! enddo
457 !
458 
459  end subroutine assemble_csr_from_pre_compute
460 
461  subroutine sll_s_compute_gyroaverage_pre_compute_polar_hermite_c1(gyro, f)
462  type(sll_t_plan_gyroaverage_polar) :: gyro
463  sll_real64, dimension(:, :), intent(inout) :: f
464  sll_int32 ::i, j, k, ell, ii(2), s
465  sll_real64::fval, delta_eta(2)
466 
467  fval = 0._f64
468  delta_eta(1) = (gyro%eta_max(1) - gyro%eta_min(1))/real(gyro%Nc(1), f64)
469  delta_eta(2) = (gyro%eta_max(2) - gyro%eta_min(2))/real(gyro%Nc(2), f64)
470 
471  call hermite_c1_coef_nat_per(f(1:gyro%Nc(1) + 1, 1:gyro%Nc(2)), gyro%deriv, gyro%Nc, gyro%interp_degree)
472 
473  do j = 1, gyro%Nc(2)
474  s = 0
475  do i = 1, gyro%Nc(1) + 1
476  fval = 0._f64
477  do k = 1, gyro%pre_compute_N(i)
478  s = s + 1
479  do ell = 1, 4
480  ii(1) = gyro%pre_compute_index(1, s)
481  ii(2) = modulo(j - 1 + gyro%pre_compute_index(2, s) + gyro%Nc(2), gyro%Nc(2))
482  fval = fval + gyro%deriv(ell, ii(1) + 1, ii(2) + 1)*gyro%pre_compute_coeff(ell, s)
483  !print *,k,ell,s,ii(1),ii(2),gyro%deriv(ell,ii(1)+1,ii(2)+1),gyro%pre_compute_coeff(ell,s)
484  end do
485  end do
486  f(i, j) = fval
487  !print *,'#',i,j,fval,gyro%pre_compute_N(i)
488  !stop
489  end do
490  !stop
491  end do
492 
493  !print *,f
494  !stop
495 
496  f(1:gyro%Nc(1) + 1, gyro%Nc(2) + 1) = f(1:gyro%Nc(1) + 1, 1)
497 
498  ! We force zero condition at radius r_min and r_max
499  f(1, 1:gyro%Nc(2) + 1) = 0._f64
500  f(gyro%Nc(1) + 1, 1:gyro%Nc(2) + 1) = 0._f64
501 
502  end subroutine sll_s_compute_gyroaverage_pre_compute_polar_hermite_c1
503 
504  subroutine sll_s_compute_gyroaverage_points_polar_spl(gyro, f, rho)
505  type(sll_t_plan_gyroaverage_polar) :: gyro
506  sll_real64, dimension(:, :), intent(inout) :: f
507  sll_real64, dimension(:, :), allocatable, target::fbords
508  sll_real64, dimension(:, :), pointer::pointer_f_bords
509  sll_real64, dimension(:), allocatable, target::buf, dnat, lnat, dper, lper, mper
510  sll_real64, dimension(:), pointer::pointer_buf, pointer_dnat, pointer_lnat
511  sll_real64, dimension(:), pointer::pointer_dper, pointer_lper, pointer_mper
512  sll_real64, intent(in)::rho
513  sll_int32 ::i, j, k
514  sll_real64::fval, sum_fval, eta(2), delta_eta(2), x(2), xx(2)
515  sll_int32 ::error
516 
517  fval = 0._f64
518  delta_eta(1) = (gyro%eta_max(1) - gyro%eta_min(1))/real(gyro%Nc(1), f64)
519  delta_eta(2) = (gyro%eta_max(2) - gyro%eta_min(2))/real(gyro%Nc(2), f64)
520 
521  sll_allocate(dnat(0:gyro%Nc(1) + 2), error)
522  sll_allocate(lnat(0:gyro%Nc(1) + 2), error)
523  sll_allocate(dper(0:gyro%Nc(2) - 1), error)
524  sll_allocate(lper(0:gyro%Nc(2) - 1), error)
525  sll_allocate(mper(0:gyro%Nc(2) - 1), error)
526 
527  sll_allocate(fbords(0:gyro%Nc(1) + 2, 0:gyro%Nc(2) - 1), error)
528  sll_allocate(buf(0:max(gyro%Nc(1) + 2, gyro%Nc(2) - 1)), error)
529 
530  call splcoefnat1d0old(dnat, lnat, gyro%Nc(1))
531  call splcoefper1d0old(dper, lper, mper, gyro%Nc(2))
532 
533  fbords(0, 0:gyro%Nc(2) - 1) = 0._f64
534  fbords(1:gyro%Nc(1) + 1, 0:gyro%Nc(2) - 1) = f(1:gyro%Nc(1) + 1, 1:gyro%Nc(2))
535  fbords(gyro%Nc(1) + 2, 0:gyro%Nc(2) - 1) = 0._f64
536 
537  pointer_f_bords => fbords
538  pointer_buf => buf
539  pointer_dnat => dnat
540  pointer_lnat => lnat
541  pointer_dper => dper
542  pointer_lper => lper
543  pointer_mper => mper
544 
545  call splcoefnatper2d(pointer_f_bords, pointer_buf, pointer_dnat, pointer_lnat, &
546  pointer_dper, pointer_lper, pointer_mper, gyro%Nc(1), gyro%Nc(2))
547 
548  do j = 1, gyro%Nc(2)
549  eta(2) = gyro%eta_min(2) + real(j - 1, f64)*delta_eta(2)
550  !eta(2)=gyro%eta_min(2) ! to uncomment for compatibility with precompute
551  do i = 1, gyro%Nc(1) + 1
552  eta(1) = gyro%eta_min(1) + real(i - 1, f64)*delta_eta(1)
553  sum_fval = 0._f64
554  do k = 1, gyro%N_points
555  x(1) = eta(1)*cos(eta(2)) + rho*gyro%points(1, k)
556  x(2) = eta(1)*sin(eta(2)) + rho*gyro%points(2, k)
557  xx(1) = sqrt(x(1)**2 + x(2)**2)
558  xx(2) = modulo(atan2(x(2), x(1)), 2._f64*sll_p_pi)
559  !xx(2)=modulo(atan2(x(2),x(1))+real(j-1,f64)*delta_eta(2),2._f64*M_PI) &
560  ! to uncomment for compatibility with precompute
561  call splnatper2d(pointer_f_bords, xx(1), gyro%eta_min(1), gyro%eta_max(1), &
562  xx(2), gyro%eta_min(2), gyro%eta_max(2), fval, gyro%Nc(1), gyro%Nc(2))
563  sum_fval = sum_fval + gyro%points(3, k)*fval
564  end do
565  f(i, j) = sum_fval
566  end do
567  end do
568  f(1:gyro%Nc(1) + 1, gyro%Nc(2) + 1) = f(1:gyro%Nc(1) + 1, 1)
569 
570  ! We force zero condition at radius r_min and r_max
571  f(1, 1:gyro%Nc(2) + 1) = 0._f64
572  f(gyro%Nc(1) + 1, 1:gyro%Nc(2) + 1) = 0._f64
573 
574  sll_deallocate_array(dnat, error)
575  sll_deallocate_array(lnat, error)
576  sll_deallocate_array(dper, error)
577  sll_deallocate_array(lper, error)
578  sll_deallocate_array(mper, error)
579 
580  sll_deallocate_array(fbords, error)
581  sll_deallocate_array(buf, error)
582 
583  end subroutine sll_s_compute_gyroaverage_points_polar_spl
584 
585  subroutine sll_s_compute_gyroaverage_points_polar_with_invar_spl(gyro, f, rho)
586  type(sll_t_plan_gyroaverage_polar) :: gyro
587  sll_real64, dimension(:, :), intent(inout) :: f
588  sll_real64, intent(in)::rho
589  sll_int32 ::i, j, k
590  sll_real64::fval, eta(2), delta_eta(2), x(2), xx(2), angle
591  sll_real64, dimension(:, :), allocatable::buf, sum_fval
592  sll_real64, dimension(:), allocatable::buf2
593  sll_int32 ::error
594 
595  sll_allocate(buf(0:gyro%Nc(1) + 2, 0:gyro%Nc(2)), error)
596  sll_allocate(buf2(0:gyro%Nc(2) - 1), error)
597  sll_allocate(sum_fval(0:gyro%Nc(1), 0:gyro%Nc(2) - 1), error)
598  sll_allocate(gyro%lunat(0:1, 0:gyro%Nc(1) + 2), error)
599  sll_allocate(gyro%luper(0:3*gyro%Nc(2) - 1), error)
600 
601  delta_eta(1) = (gyro%eta_max(1) - gyro%eta_min(1))/real(gyro%Nc(1), f64)
602  delta_eta(2) = (gyro%eta_max(2) - gyro%eta_min(2))/real(gyro%Nc(2), f64)
603 
604  call splcoefnat1d0(gyro%lunat, gyro%Nc(1))
605  call splcoefper1d0(gyro%luper, gyro%Nc(2))
606 
607  do j = 0, gyro%Nc(2)
608  buf(0, j) = 0._f64
609  do i = 0, gyro%Nc(1)
610  buf(i + 1, j) = f(i + 1, j + 1)
611  end do
612  buf(gyro%Nc(1) + 2, j) = 0._f64
613  end do
614 
615  do j = 0, gyro%Nc(2) - 1
616  call splcoefnat1d(buf(:, j), gyro%lunat, gyro%Nc(1))
617  end do
618  buf(0:gyro%Nc(1) + 2, gyro%Nc(2)) = buf(0:gyro%Nc(1) + 2, 0)
619 
620  !spline decomposition in theta
621  do i = 0, gyro%Nc(1) + 2
622  call splcoefper1d(buf(i, :), gyro%luper, gyro%Nc(2))
623  end do
624 
625  sum_fval = 0._f64
626 
627  do i = 0, gyro%Nc(1)
628  eta(1) = gyro%eta_min(1) + real(i, f64)*delta_eta(1) ! rayon du centre
629  eta(2) = gyro%eta_min(2) ! angle quelconque pour le centre
630  do k = 1, gyro%N_points
631  x(1) = eta(1)*cos(eta(2)) + rho*gyro%points(1, k)
632  x(2) = eta(1)*sin(eta(2)) + rho*gyro%points(2, k)
633  xx(1) = sqrt(x(1)**2 + x(2)**2)
634  xx(2) = atan2(x(2), x(1))
635  do j = 0, gyro%Nc(2) - 1
636  call splnat1d(buf(:, j), xx(1), gyro%eta_min(1), gyro%eta_max(1), buf2(j), gyro%Nc(1))
637  end do
638  !call splcoefper1d(buf2,gyro%luper,gyro%N(2))
639  do j = 0, gyro%Nc(2) - 1
640  angle = modulo(xx(2) + real(j, f64)*delta_eta(2), 2._f64*sll_p_pi)
641  call splper1d(buf2, angle, gyro%eta_min(2), gyro%eta_max(2), fval, gyro%Nc(2))
642  sum_fval(i, j) = sum_fval(i, j) + gyro%points(3, k)*fval
643  end do
644  end do
645  end do
646 
647  f(1:gyro%Nc(1) + 1, 1:gyro%Nc(2)) = sum_fval(0:gyro%Nc(1), 0:gyro%Nc(2) - 1)
648  f(1:gyro%Nc(1) + 1, gyro%Nc(2) + 1) = f(1:gyro%Nc(1) + 1, 1)
649 
650  ! We force zero condition at radius r_min and r_max
651  f(1, 1:gyro%Nc(2) + 1) = 0._f64
652  f(gyro%Nc(1) + 1, 1:gyro%Nc(2) + 1) = 0._f64
653 
654  sll_deallocate_array(gyro%lunat, error)
655  sll_deallocate_array(gyro%luper, error)
656  sll_deallocate_array(buf, error)
657  sll_deallocate_array(buf2, error)
658  sll_deallocate_array(sum_fval, error)
659 
660  end subroutine sll_s_compute_gyroaverage_points_polar_with_invar_spl
661 
662  subroutine sll_s_pre_compute_gyroaverage_polar_spl(gyro, rho)
663  type(sll_t_plan_gyroaverage_polar) :: gyro
664  sll_real64, intent(in)::rho
665  sll_int32, dimension(:, :), allocatable :: buf
666  sll_int32 ::i, j, k, ell_1, ell_2, ii(2), s, nb, ind(2)
667  sll_real64::val(-1:2, -1:2), eta_star(2), eta(2), delta_eta(2), x(2)
668  sll_int32 ::error
669 
670  delta_eta(1) = (gyro%eta_max(1) - gyro%eta_min(1))/real(gyro%Nc(1), f64)
671  delta_eta(2) = (gyro%eta_max(2) - gyro%eta_min(2))/real(gyro%Nc(2), f64)
672 
673  sll_allocate(buf(0:gyro%Nc(1) + 2, 0:gyro%Nc(2) - 1), error)
674 
675  sll_allocate(gyro%pre_compute_N(gyro%Nc(1) + 1), error)
676 
677  eta(2) = gyro%eta_min(2)
678  buf = 0
679  nb = 0
680  do i = 1, gyro%Nc(1) + 1
681  eta(1) = gyro%eta_min(1) + real(i - 1, f64)*delta_eta(1)
682  s = 0
683  do k = 1, gyro%N_points
684  x(1) = eta(1)*cos(eta(2)) + rho*gyro%points(1, k)
685  x(2) = eta(1)*sin(eta(2)) + rho*gyro%points(2, k)
686  call localize_polar(x, gyro%eta_min, gyro%eta_max, ii, eta_star, gyro%Nc)
687  do ell_2 = -1, 2
688  ind(2) = modulo(ii(2) + ell_2, gyro%Nc(2))
689  do ell_1 = -1, 2
690  ind(1) = ii(1) + 1 + ell_1
691  if (buf(ind(1), ind(2)) .ne. i) then
692  s = s + 1
693  buf(ind(1), ind(2)) = i
694  end if
695  end do
696  end do
697  end do
698  gyro%pre_compute_N(i) = s
699  nb = nb + s
700  end do
701 
702  print *, '#N_points pre_compute=', nb, gyro%Nc(1), nb/gyro%Nc(1)
703 
704  sll_allocate(gyro%pre_compute_index(1:2, nb), error)
705  sll_allocate(gyro%pre_compute_coeff_spl(nb), error)
706 
707  buf = 0
708  nb = 0
709  s = 0
710  val = 0._f64
711  eta(2) = gyro%eta_min(2)
712  do i = 1, gyro%Nc(1) + 1
713  eta(1) = gyro%eta_min(1) + real(i - 1, f64)*delta_eta(1)
714  do k = 1, gyro%N_points
715  x(1) = eta(1)*cos(eta(2)) + rho*gyro%points(1, k)
716  x(2) = eta(1)*sin(eta(2)) + rho*gyro%points(2, k)
717  call localize_polar(x, gyro%eta_min, gyro%eta_max, ii, eta_star, gyro%Nc)
718 
719  call contribution_spl(eta_star, val)
720 
721  val = val*gyro%points(3, k)
722 
723  do ell_2 = -1, 2
724  ind(2) = modulo(ii(2) + ell_2, gyro%Nc(2))
725  do ell_1 = -1, 2
726  ind(1) = ii(1) + 1 + ell_1
727  j = buf(ind(1), ind(2))
728  if (j <= nb) then
729  s = s + 1
730  buf(ind(1), ind(2)) = s
731  gyro%pre_compute_coeff_spl(s) = val(ell_1, ell_2)
732  gyro%pre_compute_index(1, s) = ind(1)
733  gyro%pre_compute_index(2, s) = ind(2)
734  else
735  gyro%pre_compute_coeff_spl(j) = &
736  gyro%pre_compute_coeff_spl(j) + val(ell_1, ell_2)
737  end if
738  end do
739  end do
740  end do
741  nb = s
742  end do
743 
744  print *, '#min/max=', minval(gyro%pre_compute_coeff_spl(:)), maxval(gyro%pre_compute_coeff_spl(:))
745  sll_deallocate_array(buf, error)
746  end subroutine sll_s_pre_compute_gyroaverage_polar_spl
747 
748  subroutine sll_s_compute_gyroaverage_pre_compute_polar_spl(gyro, f)
749  type(sll_t_plan_gyroaverage_polar) :: gyro
750  sll_real64, dimension(:, :), intent(inout) :: f
751  sll_real64, dimension(:, :), allocatable, target::fbords
752  sll_real64, dimension(:, :), pointer::pointer_f_bords
753  sll_real64, dimension(:), allocatable, target::buf, dnat, lnat, dper, lper, mper
754  sll_real64, dimension(:), pointer::pointer_buf, pointer_dnat, pointer_lnat
755  sll_real64, dimension(:), pointer::pointer_dper, pointer_lper, pointer_mper
756  sll_real64, dimension(:, :), allocatable ::tmp_f
757  sll_int32 ::i, j, k, ii(2), s
758  sll_real64::fval, delta_eta(2)
759  sll_int32 ::error
760 
761  sll_allocate(dnat(0:gyro%Nc(1) + 2), error)
762  sll_allocate(lnat(0:gyro%Nc(1) + 2), error)
763  sll_allocate(dper(0:gyro%Nc(2) - 1), error)
764  sll_allocate(lper(0:gyro%Nc(2) - 1), error)
765  sll_allocate(mper(0:gyro%Nc(2) - 1), error)
766  sll_allocate(fbords(0:gyro%Nc(1) + 2, 0:gyro%Nc(2) - 1), error)
767  sll_allocate(buf(0:max(gyro%Nc(1) + 2, gyro%Nc(2) - 1)), error)
768  sll_allocate(tmp_f(1:gyro%Nc(1) + 1, 1:gyro%Nc(2)), error)
769 
770  fval = 0._f64
771  delta_eta(1) = (gyro%eta_max(1) - gyro%eta_min(1))/real(gyro%Nc(1), f64)
772  delta_eta(2) = (gyro%eta_max(2) - gyro%eta_min(2))/real(gyro%Nc(2), f64)
773 
774  call splcoefnat1d0old(dnat, lnat, gyro%Nc(1))
775  call splcoefper1d0old(dper, lper, mper, gyro%Nc(2))
776 
777  fbords(0, 0:gyro%Nc(2) - 1) = 0._f64
778  fbords(1:gyro%Nc(1) + 1, 0:gyro%Nc(2) - 1) = f(1:gyro%Nc(1) + 1, 1:gyro%Nc(2))
779  fbords(gyro%Nc(1) + 2, 0:gyro%Nc(2) - 1) = 0._f64
780 
781  pointer_f_bords => fbords
782  pointer_buf => buf
783  pointer_dnat => dnat
784  pointer_lnat => lnat
785  pointer_dper => dper
786  pointer_lper => lper
787  pointer_mper => mper
788 
789  call splcoefnatper2d(pointer_f_bords, pointer_buf, pointer_dnat, pointer_lnat, &
790  pointer_dper, pointer_lper, pointer_mper, gyro%Nc(1), gyro%Nc(2))
791 
792  do j = 1, gyro%Nc(2)
793  s = 0
794  do i = 1, gyro%Nc(1) + 1
795  fval = 0._f64
796  do k = 1, gyro%pre_compute_N(i)
797  s = s + 1
798  ii(1) = gyro%pre_compute_index(1, s)
799  ii(2) = modulo(j - 1 + gyro%pre_compute_index(2, s), gyro%Nc(2))
800  fval = fval + pointer_f_bords(ii(1), ii(2))*gyro%pre_compute_coeff_spl(s)
801  end do
802  tmp_f(i, j) = fval
803  end do
804  end do
805 
806  f(1:gyro%Nc(1) + 1, 1:gyro%Nc(2)) = tmp_f(1:gyro%Nc(1) + 1, 1:gyro%Nc(2))
807  f(1:gyro%Nc(1) + 1, gyro%Nc(2) + 1) = f(1:gyro%Nc(1) + 1, 1)
808 
809  ! We force zero condition at radius r_min and r_max
810  f(1, 1:gyro%Nc(2) + 1) = 0._f64
811  f(gyro%Nc(1) + 1, 1:gyro%Nc(2) + 1) = 0._f64
812 
813  sll_deallocate_array(dnat, error)
814  sll_deallocate_array(lnat, error)
815  sll_deallocate_array(dper, error)
816  sll_deallocate_array(lper, error)
817  sll_deallocate_array(mper, error)
818  sll_deallocate_array(fbords, error)
819  sll_deallocate_array(buf, error)
820  sll_deallocate_array(tmp_f, error)
821 
822  end subroutine sll_s_compute_gyroaverage_pre_compute_polar_spl
823 
824  subroutine sll_s_pre_compute_gyroaverage_polar_spl_fft(gyro, rho)
825  type(sll_t_plan_gyroaverage_polar) :: gyro
826  sll_real64, intent(in)::rho
827  sll_int32, dimension(:, :), allocatable :: buf
828  sll_real64, dimension(:), allocatable::buf_fft
829  sll_int32 ::i, j, k, ell_1, ell_2, ii(2), s, nb, ind(2)
830  sll_real64::val(-1:2, -1:2), eta_star(2), eta(2), delta_eta(2), x(2)
831  sll_int32 ::error
832 
833  delta_eta(1) = (gyro%eta_max(1) - gyro%eta_min(1))/real(gyro%Nc(1), f64)
834  delta_eta(2) = (gyro%eta_max(2) - gyro%eta_min(2))/real(gyro%Nc(2), f64)
835 
836  sll_allocate(buf(1:gyro%Nc(1) + 1, 0:gyro%Nc(1) + 2), error)
837 
838  sll_allocate(gyro%A_fft(1:gyro%Nc(1) + 1, 0:gyro%Nc(1) + 2, 0:gyro%Nc(2) - 1), error)
839  sll_allocate(buf_fft(1:2*gyro%Nc(2) + 15), error)
840 
841  eta(2) = gyro%eta_min(2)
842  buf = 0
843  nb = 0
844  do i = 1, gyro%Nc(1) + 1
845  eta(1) = gyro%eta_min(1) + real(i - 1, f64)*delta_eta(1)
846  s = 0
847  do k = 1, gyro%N_points
848  x(1) = eta(1)*cos(eta(2)) + rho*gyro%points(1, k)
849  x(2) = eta(1)*sin(eta(2)) + rho*gyro%points(2, k)
850  call localize_polar(x, gyro%eta_min, gyro%eta_max, ii, eta_star, gyro%Nc)
851  do ell_2 = -1, 2
852  ind(2) = modulo(ii(2) + ell_2, gyro%Nc(2))
853  do ell_1 = -1, 2
854  ind(1) = ii(1) + 1 + ell_1
855  if (buf(i, ind(1)) .ne. 1) then
856  s = s + 1
857  buf(i, ind(1)) = 1
858  end if
859  end do
860  end do
861  end do
862  nb = nb + s
863  end do
864 
865  sll_allocate(gyro%pre_compute_index(1:2, nb), error)
866 
867  gyro%A_fft = 0._f64
868  buf = 0
869  s = 0
870  val = 0._f64
871  eta(2) = gyro%eta_min(2)
872  do i = 1, gyro%Nc(1) + 1
873  eta(1) = gyro%eta_min(1) + real(i - 1, f64)*delta_eta(1)
874  do k = 1, gyro%N_points
875  x(1) = eta(1)*cos(eta(2)) + rho*gyro%points(1, k)
876  x(2) = eta(1)*sin(eta(2)) + rho*gyro%points(2, k)
877  call localize_polar(x, gyro%eta_min, gyro%eta_max, ii, eta_star, gyro%Nc)
878 
879  call contribution_spl(eta_star, val)
880 
881  val = val*gyro%points(3, k)
882 
883  do ell_2 = -1, 2
884  ind(2) = modulo(ii(2) + ell_2, gyro%Nc(2))
885  do ell_1 = -1, 2
886  ind(1) = ii(1) + 1 + ell_1
887  j = buf(i, ind(1))
888  if (j == 0) then
889  s = s + 1
890  buf(i, ind(1)) = 1
891  gyro%pre_compute_index(1, s) = i
892  gyro%pre_compute_index(2, s) = ind(1)
893  gyro%A_fft(i, ind(1), ind(2)) = val(ell_1, ell_2)
894  else
895  gyro%A_fft(i, ind(1), ind(2)) = &
896  gyro%A_fft(i, ind(1), ind(2)) + val(ell_1, ell_2)
897  end if
898  end do
899  end do
900  end do
901  end do
902 
903  call dffti(gyro%Nc(2), buf_fft)
904  do s = 1, nb
905  call dfftf(gyro%Nc(2), gyro%A_fft(gyro%pre_compute_index(1, s), gyro%pre_compute_index(2, s), :), buf_fft)
906  end do
907 
908  sll_deallocate_array(buf, error)
909  sll_deallocate_array(buf_fft, error)
910 
911  end subroutine sll_s_pre_compute_gyroaverage_polar_spl_fft
912 
913  subroutine sll_s_compute_gyroaverage_pre_compute_polar_spl_fft(gyro, f)
914  type(sll_t_plan_gyroaverage_polar) :: gyro
915  sll_real64, dimension(:, :), intent(inout) :: f
916  sll_real64, dimension(:, :), allocatable, target::fbords
917  sll_real64, dimension(:, :), pointer::pointer_f_bords
918  sll_real64, dimension(:), allocatable, target::buf, dnat, lnat, dper, lper, mper
919  sll_real64, dimension(:), pointer::pointer_buf, pointer_dnat, pointer_lnat
920  sll_real64, dimension(:), pointer::pointer_dper, pointer_lper, pointer_mper
921  sll_real64, dimension(:, :), allocatable ::tmp_f
922  sll_real64, dimension(:), allocatable::buf_fft
923  sll_int32 ::i, j, s
924  sll_real64::fval, delta_eta(2)
925  sll_int32 ::error
926 
927  sll_allocate(dnat(0:gyro%Nc(1) + 2), error)
928  sll_allocate(lnat(0:gyro%Nc(1) + 2), error)
929  sll_allocate(dper(0:gyro%Nc(2) - 1), error)
930  sll_allocate(lper(0:gyro%Nc(2) - 1), error)
931  sll_allocate(mper(0:gyro%Nc(2) - 1), error)
932  sll_allocate(fbords(0:gyro%Nc(1) + 2, 0:gyro%Nc(2) - 1), error)
933  sll_allocate(buf(0:max(gyro%Nc(1) + 2, gyro%Nc(2) - 1)), error)
934  sll_allocate(tmp_f(1:gyro%Nc(1) + 1, 0:gyro%Nc(2) - 1), error)
935  sll_allocate(buf_fft(1:2*gyro%Nc(2) + 15), error)
936 
937  fval = 0._f64
938  tmp_f = 0._f64
939  delta_eta(1) = (gyro%eta_max(1) - gyro%eta_min(1))/real(gyro%Nc(1), f64)
940  delta_eta(2) = (gyro%eta_max(2) - gyro%eta_min(2))/real(gyro%Nc(2), f64)
941 
942  call splcoefnat1d0old(dnat, lnat, gyro%Nc(1))
943  call splcoefper1d0old(dper, lper, mper, gyro%Nc(2))
944 
945  fbords(0, 0:gyro%Nc(2) - 1) = 0._f64
946  fbords(1:gyro%Nc(1) + 1, 0:gyro%Nc(2) - 1) = f(1:gyro%Nc(1) + 1, 1:gyro%Nc(2))
947  fbords(gyro%Nc(1) + 2, 0:gyro%Nc(2) - 1) = 0._f64
948 
949  pointer_f_bords => fbords
950  pointer_buf => buf
951  pointer_dnat => dnat
952  pointer_lnat => lnat
953  pointer_dper => dper
954  pointer_lper => lper
955  pointer_mper => mper
956 
957  call splcoefnatper2d(pointer_f_bords, pointer_buf, pointer_dnat, pointer_lnat, &
958  pointer_dper, pointer_lper, pointer_mper, gyro%Nc(1), gyro%Nc(2))
959 
960  call dffti(gyro%Nc(2), buf_fft)
961  do i = 0, gyro%Nc(1) + 2
962  call dfftf(gyro%Nc(2), pointer_f_bords(i, :), buf_fft)
963  end do
964 
965  do j = 0, gyro%Nc(2) - 1
966  do s = 1, size(gyro%pre_compute_index(1, :))
967  tmp_f(gyro%pre_compute_index(1, s), j) = &
968  tmp_f(gyro%pre_compute_index(1, s), j) + &
969  gyro%A_fft(gyro%pre_compute_index(1, s), &
970  gyro%pre_compute_index(2, s), j)* &
971  pointer_f_bords(gyro%pre_compute_index(2, s), j)
972  end do
973  end do
974 
975  do i = 1, gyro%Nc(1) + 1
976  call dfftb(gyro%Nc(2), tmp_f(i, :), buf_fft)
977  end do
978 
979  tmp_f = tmp_f/real(gyro%Nc(2), f64)
980 
981  f(1:gyro%Nc(1) + 1, 1:gyro%Nc(2)) = tmp_f(1:gyro%Nc(1) + 1, 0:gyro%Nc(2) - 1)
982  f(1:gyro%Nc(1) + 1, gyro%Nc(2) + 1) = f(1:gyro%Nc(1) + 1, 1)
983 
984  ! We force zero condition at radius r_min and r_max
985  f(1, 1:gyro%Nc(2) + 1) = 0._f64
986  f(gyro%Nc(1) + 1, 1:gyro%Nc(2) + 1) = 0._f64
987 
988  sll_deallocate_array(dnat, error)
989  sll_deallocate_array(lnat, error)
990  sll_deallocate_array(dper, error)
991  sll_deallocate_array(lper, error)
992  sll_deallocate_array(mper, error)
993  sll_deallocate_array(fbords, error)
994  sll_deallocate_array(buf, error)
995  sll_deallocate_array(buf_fft, error)
996  sll_deallocate_array(tmp_f, error)
997 
998  end subroutine sll_s_compute_gyroaverage_pre_compute_polar_spl_fft
999 
1000  subroutine sll_s_compute_gyroaverage_pade_polar(gyro, f, rho)
1001  type(sll_t_plan_gyroaverage_polar) :: gyro
1002  sll_real64, dimension(:, :), intent(inout) :: f
1003  sll_real64, dimension(:, :), allocatable :: fcomp
1004  sll_real64, dimension(:), allocatable :: buf, diagm1, diag, diagp1
1005  sll_real64, intent(in)::rho
1006  sll_int32 ::i, k
1007  sll_real64::dr
1008  sll_int32 ::error
1009 
1010  dr = (gyro%eta_max(1) - gyro%eta_min(1))/real(gyro%Nc(1), f64)
1011 
1012  sll_allocate(buf(1:2*gyro%Nc(2) + 15), error)
1013  sll_allocate(fcomp(1:gyro%Nc(1) + 1, 1:gyro%Nc(2)), error)
1014  sll_allocate(diagm1(1:gyro%Nc(1) + 1), error)
1015  sll_allocate(diag(1:gyro%Nc(1) + 1), error)
1016  sll_allocate(diagp1(1:gyro%Nc(1) + 1), error)
1017 
1018  fcomp(1:gyro%Nc(1) + 1, 1:gyro%Nc(2)) = f(1:gyro%Nc(1) + 1, 1:gyro%Nc(2))
1019 
1020  !*** Perform FFT 1D in theta direction of ***
1021  !*** the system solution ***
1022  call dffti(gyro%Nc(2), buf)
1023  do i = 1, gyro%Nc(1) + 1
1024  call dfftf(gyro%Nc(2), fcomp(i, :), buf)
1025  end do
1026  fcomp = fcomp/real(gyro%Nc(2), f64)
1027 
1028  !***POISSON
1029  do k = 1, gyro%Nc(2)
1030  do i = 1, gyro%Nc(1)
1031  diagm1(i + 1) = -(rho**2/4._f64)*(1._f64/dr**2 - 1._f64/(2._f64*dr*(gyro%eta_min(1) + &
1032  (gyro%eta_max(1) - gyro%eta_min(1))*real(i, f64) &
1033  /real(gyro%Nc(1), f64))))
1034  diag(i) = 1._f64 - (rho**2/4._f64)*(-(2._f64/dr**2) - (real(floor(k/2._f64), kind=f64)/ &
1035  (gyro%eta_min(1) + (gyro%eta_max(1) - gyro%eta_min(1))* &
1036  real(i - 1, f64)/real(gyro%Nc(1), f64)))**2)
1037  diagp1(i) = -(rho**2/4._f64)*(1._f64/dr**2 + 1._f64/(2._f64*dr*(gyro%eta_min(1) + &
1038  (gyro%eta_max(1) - gyro%eta_min(1))*real(i - 1, f64) &
1039  /real(gyro%Nc(1), f64))))
1040  end do
1041  diagm1(1) = 0._f64
1042  diagp1(gyro%Nc(1) + 1) = 0._f64
1043  diag(1) = 1._f64
1044  diag(gyro%Nc(1) + 1) = 1._f64
1045  !*** Dirichlet boundary conditions ***
1046  diagp1(1) = 0._f64
1047  diagm1(gyro%Nc(1) + 1) = 0._f64
1048  !*** Neumann boundary conditions ***
1049  ! diagp1(1)=-1._f64
1050  ! diagm1(gyro%Nc(1)+1)=-1._f64
1051  call solve_tridiag(diagm1, diag, diagp1, fcomp(1:gyro%Nc(1) + 1, k), f(1:gyro%Nc(1) + 1, k), gyro%Nc(1) + 1)
1052  end do
1053 
1054  !*** Perform FFT 1D inverse ***
1055  do i = 1, gyro%Nc(1) + 1
1056  call dfftb(gyro%Nc(2), f(i, 1:gyro%Nc(2)), buf)
1057  end do
1058 
1059  !*** duplicate periodic value ***
1060  f(1:gyro%Nc(1) + 1, gyro%Nc(2) + 1) = f(1:gyro%Nc(1) + 1, 1)
1061 
1062  end subroutine sll_s_compute_gyroaverage_pade_polar
1063 
1064  subroutine sll_s_compute_gyroaverage_pade_high_order_polar(gyro, f, rho, order)
1065  type(sll_t_plan_gyroaverage_polar) :: gyro
1066  sll_real64, dimension(:, :), intent(inout) :: f
1067  sll_real64, dimension(:, :), allocatable :: fcomp
1068  sll_real64, dimension(:), allocatable :: buf
1069  sll_real64, dimension(:), allocatable :: diagm1_left, diag_left, diagp1_left, diagm2_left, diagp2_left
1070  sll_real64, dimension(:), allocatable :: diagm1_right, diag_right, diagp1_right, fbuf
1071  sll_int32, intent(in) :: order(2)
1072  sll_real64, intent(in)::rho
1073  sll_int32 ::i, k
1074  sll_real64::dr, a, b, c, ri, rip1, rip2, nfloat
1075  sll_int32 ::error
1076 
1077  if ((order(1) == 0) .and. (order(2) == 2)) then
1078  a = -(rho**2)/4._f64
1079  b = 0._f64
1080  c = 0._f64
1081  elseif ((order(1) == 0) .and. (order(2) == 4)) then
1082  a = -(rho**2)/4._f64
1083  b = 3._f64*(rho**4)/64._f64
1084  c = 0._f64
1085  elseif ((order(1) == 2) .and. (order(2) == 4)) then
1086  a = -2._f64*(rho**2)/27._f64
1087  b = 5._f64*(rho**4)/1728._f64
1088  c = 19._f64*(rho**2)/108._f64
1089  else
1090  print *, '#bad value of pade_order=', order
1091  print *, '#not implemented'
1092  print *, '#in sll_s_compute_gyroaverage_pade_high_order_polar'
1093  stop
1094  end if
1095 
1096  dr = (gyro%eta_max(1) - gyro%eta_min(1))/real(gyro%Nc(1), f64)
1097 
1098  sll_allocate(buf(1:2*gyro%Nc(2) + 15), error)
1099  sll_allocate(fcomp(1:gyro%Nc(1) + 1, 1:gyro%Nc(2)), error)
1100  sll_allocate(diagm2_left(1:gyro%Nc(1) + 1), error)
1101  sll_allocate(diagm1_left(1:gyro%Nc(1) + 1), error)
1102  sll_allocate(diag_left(1:gyro%Nc(1) + 1), error)
1103  sll_allocate(diagp1_left(1:gyro%Nc(1) + 1), error)
1104  sll_allocate(diagp2_left(1:gyro%Nc(1) + 1), error)
1105  sll_allocate(diagm1_right(1:gyro%Nc(1) + 1), error)
1106  sll_allocate(diag_right(1:gyro%Nc(1) + 1), error)
1107  sll_allocate(diagp1_right(1:gyro%Nc(1) + 1), error)
1108  sll_allocate(fbuf(1:gyro%Nc(1) + 1), error)
1109 
1110  fcomp(1:gyro%Nc(1) + 1, 1:gyro%Nc(2)) = f(1:gyro%Nc(1) + 1, 1:gyro%Nc(2))
1111 
1112  !*** Perform FFT 1D in theta direction of ***
1113  !*** the system solution ***
1114  call dffti(gyro%Nc(2), buf)
1115  do i = 1, gyro%Nc(1) + 1
1116  call dfftf(gyro%Nc(2), fcomp(i, :), buf)
1117  end do
1118  fcomp = fcomp/real(gyro%Nc(2), f64)
1119 
1120  !***POISSON
1121  do k = 1, gyro%Nc(2)
1122  nfloat = real(floor(k/2._f64), kind=f64)
1123  do i = 1, gyro%Nc(1) - 1
1124  ri = gyro%eta_min(1) + (gyro%eta_max(1) - gyro%eta_min(1))*real(i - 1, f64)/real(gyro%Nc(1), f64)
1125  rip1 = gyro%eta_min(1) + (gyro%eta_max(1) - gyro%eta_min(1))*real(i, f64)/real(gyro%Nc(1), f64)
1126  rip2 = gyro%eta_min(1) + (gyro%eta_max(1) - gyro%eta_min(1))*real(i + 1, f64)/real(gyro%Nc(1), f64)
1127 
1128  ! gauche
1129 
1130  diagm2_left(i + 2) = b/dr**4 + &
1131  (2._f64*b/rip2)*(-1._f64/(2._f64*dr**3))
1132 
1133  diagm1_left(i + 1) = -4._f64*b/dr**4 + &
1134  (2._f64*b/rip1)*(2._f64/(2._f64*dr**3)) + &
1135  (a - (b*(2._f64*(nfloat**2) + 1._f64)/rip1**2))*(1._f64/dr**2) + &
1136  (a/rip1 + b*(2._f64*(nfloat**2) + 1._f64)/(rip1**3))*(-1._f64/(2._f64*dr))
1137 
1138  diag_left(i) = 6._f64*b/dr**4 + &
1139  (a - (b*(2._f64*(nfloat**2) + 1._f64)/ri**2))*(-2._f64/dr**2) + &
1140  (1._f64 - a*(nfloat**2)/(rip1**2) + b*(nfloat**2)*(nfloat**2 - 4._f64)/(ri**4))
1141 
1142  diagp1_left(i) = -4._f64*b/dr**4 + &
1143  (2._f64*b/ri)*(-2._f64/(2._f64*dr**3)) + &
1144  (a - (b*(2._f64*(nfloat**2) + 1._f64)/ri**2))*(1._f64/dr**2) + &
1145  (a/ri + b*(2._f64*(nfloat**2) + 1._f64)/(ri**3))*(1._f64/(2._f64*dr))
1146 
1147  diagp2_left(i) = b/dr**4 + &
1148  (2._f64*b/ri)*(1._f64/(2._f64*dr**3))
1149 
1150  ! droite
1151 
1152  diagm1_right(i + 1) = c*(1._f64/dr**2) + &
1153  (c/rip1)*(-1._f64/(2._f64*dr))
1154 
1155  diag_right(i) = c*(-2._f64/dr**2) + &
1156  (1._f64 - c*(nfloat**2)/(rip1**2))
1157 
1158  diagp1_right(i) = c*(1._f64/dr**2) + &
1159  (c/ri)*(1._f64/(2._f64*dr))
1160 
1161  end do
1162  diagm1_left(1) = 0._f64
1163  diagm2_left(1) = 0._f64
1164  diagm2_left(2) = 0._f64
1165  diagp1_left(gyro%Nc(1) + 1) = 0._f64
1166  diagp2_left(gyro%Nc(1) + 1) = 0._f64
1167  diagp2_left(gyro%Nc(1)) = 0._f64
1168  !*** Dirichlet boundary conditions ***
1169  diag_left(1) = 1._f64
1170  diagp1_left(1) = 0._f64
1171  diagp2_left(1) = 0._f64
1172  diagm1_left(2) = 0._f64
1173  diag_left(2) = 1._f64
1174  diagp1_left(2) = 0._f64
1175  diagp2_left(2) = 0._f64
1176  diagm1_left(gyro%Nc(1)) = 0._f64
1177  diagm2_left(gyro%Nc(1)) = 0._f64
1178  diag_left(gyro%Nc(1)) = 1._f64
1179  diagp1_left(gyro%Nc(1)) = 0._f64
1180  diagm1_left(gyro%Nc(1) + 1) = 0._f64
1181  diagm2_left(gyro%Nc(1) + 1) = 0._f64
1182  diag_left(gyro%Nc(1) + 1) = 1._f64
1183 
1184  diagm1_right(1) = 0._f64
1185  diagp1_right(gyro%Nc(1) + 1) = 0._f64
1186  !*** Dirichlet boundary conditions ***
1187  diag_right(1) = 1._f64
1188  diagp1_right(1) = 0._f64
1189  diagm1_right(2) = 0._f64
1190  diag_right(2) = 1._f64
1191  diagp1_right(2) = 0._f64
1192  diagm1_right(gyro%Nc(1)) = 0._f64
1193  diag_right(gyro%Nc(1)) = 1._f64
1194  diagp1_right(gyro%Nc(1)) = 0._f64
1195  diagm1_right(gyro%Nc(1) + 1) = 0._f64
1196  diag_right(gyro%Nc(1) + 1) = 1._f64
1197 
1198  ! droite
1199  fbuf(1:gyro%Nc(1) + 1) = fcomp(1:gyro%Nc(1) + 1, k)
1200  do i = 2, gyro%Nc(1)
1201  fcomp(i, k) = diagm1_right(i)*fbuf(i - 1) + diag_right(i)*fbuf(i) + diagp1_right(i)*fbuf(i + 1)
1202  end do
1203 
1204  !*** Solve sll_s_penta ***
1205 call sll_s_penta(gyro%Nc(1)+1,diagm2_left,diagm1_left,diag_left,diagp1_left,diagp2_left,fcomp(1:gyro%Nc(1)+1,k),f(1:gyro%Nc(1)+1,k))
1206  end do
1207 
1208  !*** Perform FFT 1D inverse ***
1209  do i = 1, gyro%Nc(1) + 1
1210  call dfftb(gyro%Nc(2), f(i, 1:gyro%Nc(2)), buf)
1211  end do
1212 
1213  !*** duplicate periodic value ***
1214  f(1:gyro%Nc(1) + 1, gyro%Nc(2) + 1) = f(1:gyro%Nc(1) + 1, 1)
1215 
1216  end subroutine sll_s_compute_gyroaverage_pade_high_order_polar
1217 
1218  subroutine hermite_coef_nat_per(f, buf3d, N, d)
1219  sll_int32, intent(in)::n(2), d(2)
1220  sll_real64, dimension(N(1) + 1, N(2)), intent(in)::f
1221  sll_real64, dimension(9, N(1) + 1, N(2) + 1), intent(out)::buf3d
1222  sll_real64 ::w_left_1(-d(1)/2:(d(1) + 1)/2), w_right_1((-d(1) + 1)/2:d(1)/2 + 1)
1223  sll_real64 ::w_left_2(-d(2)/2:(d(2) + 1)/2), w_right_2((-d(2) + 1)/2:d(2)/2 + 1)
1224  sll_real64 ::tmp
1225  sll_int32 ::i, j, ii, r_left(2), r_right(2), s_left(2), s_right(2), ind
1226  r_left = -d/2
1227  s_left = (d + 1)/2
1228  r_right = (-d + 1)/2
1229  s_right = d/2 + 1
1230 
1231  call compute_w_hermite(w_left_1, r_left(1), s_left(1))
1232  call compute_w_hermite(w_left_2, r_left(2), s_left(2))
1233  if ((2*(d(1)/2) - d(1)) == 0) then
1234  w_right_1(r_right(1):s_right(1)) = w_left_1(r_left(1):s_left(1))
1235  else
1236  w_right_1(r_right(1):s_right(1)) = -w_left_1(s_left(1):r_left(1):-1)
1237  end if
1238 
1239  if ((2*(d(2)/2) - d(2)) == 0) then
1240  w_right_2(r_right(2):s_right(2)) = w_left_2(r_left(2):s_left(2))
1241  else
1242  w_right_2(r_right(2):s_right(2)) = -w_left_2(s_left(2):r_left(2):-1)
1243  end if
1244 
1245  !print *,'w(',r_left(1),':',s_left(1),')=',w_left_1(r_left(1):s_left(1))
1246  !print *,'w(',r_right(1),':',s_right(1),')=',w_right_1(r_right(1):s_right(1))
1247 
1248  do j = 1, n(2)
1249  do i = 1, n(1) + 1
1250  buf3d(1, i, j) = f(i, j) !f(0,0)
1251  tmp = 0._f64
1252  do ii = r_left(1), s_left(1)
1253  ind = i + ii; if (ind > n(1) + 1) ind = 2*(n(1) + 1) - ind; if (ind < 1) ind = 2 - ind
1254  tmp = tmp + w_left_1(ii)*f(ind, j)
1255  end do
1256  buf3d(2, i, j) = tmp !fx(0,0)
1257  tmp = 0._f64
1258  do ii = r_right(1), s_right(1)
1259  ind = i + ii; if (ind > n(1) + 1) ind = 2*(n(1) + 1) - ind; if (ind < 1) ind = 2 - ind
1260  tmp = tmp + w_right_1(ii)*f(ind, j)
1261  end do
1262  buf3d(3, i, j) = tmp !fx(1,0)
1263  end do
1264  end do
1265  do i = 1, n(1) + 1
1266  do j = 1, n(2)
1267  tmp = 0._f64
1268  do ii = r_left(2), s_left(2)
1269  ind = modulo(j + ii - 1 + n(2), n(2)) + 1
1270  tmp = tmp + w_left_2(ii)*f(i, ind)
1271  end do
1272  buf3d(4, i, j) = tmp !fy(0,0)
1273  tmp = 0._f64
1274  do ii = r_right(2), s_right(2)
1275  ind = modulo(j + ii - 1 + n(2), n(2)) + 1
1276  tmp = tmp + w_right_2(ii)*f(i, ind)
1277  end do
1278  buf3d(5, i, j) = tmp !fy(0,1)
1279  end do
1280  end do
1281 
1282  do j = 1, n(2)
1283  do i = 1, n(1) + 1
1284  tmp = 0._f64
1285  do ii = r_left(1), s_left(1)
1286  ind = i + ii; if (ind > n(1) + 1) ind = 2*(n(1) + 1) - ind; if (ind < 1) ind = 2 - ind
1287  tmp = tmp + w_left_1(ii)*buf3d(4, ind, j)
1288  end do
1289  buf3d(6, i, j) = tmp !fxy(0,0)
1290  tmp = 0._f64
1291  do ii = r_right(1), s_right(1)
1292  ind = i + ii; if (ind > n(1) + 1) ind = 2*(n(1) + 1) - ind; if (ind < 1) ind = 2 - ind
1293  tmp = tmp + w_right_1(ii)*buf3d(4, ind, j)
1294  end do
1295  buf3d(7, i, j) = tmp !fxy(1,0)
1296  tmp = 0._f64
1297  do ii = r_left(1), s_left(1)
1298  ind = i + ii; if (ind > n(1) + 1) ind = 2*(n(1) + 1) - ind; if (ind < 1) ind = 2 - ind
1299  tmp = tmp + w_left_1(ii)*buf3d(5, ind, j)
1300  end do
1301  buf3d(8, i, j) = tmp !fxy(0,1)
1302  tmp = 0._f64
1303  do ii = r_right(1), s_right(1)
1304  ind = i + ii; if (ind > n(1) + 1) ind = 2*(n(1) + 1) - ind; if (ind < 1) ind = 2 - ind
1305  tmp = tmp + w_right_1(ii)*buf3d(5, ind, j)
1306  end do
1307  buf3d(9, i, j) = tmp !fxy(1,1)
1308  end do
1309  end do
1310 
1311  buf3d(:, :, n(2) + 1) = buf3d(:, :, 1)
1312 
1313  end subroutine hermite_coef_nat_per
1314 
1315  subroutine hermite_c1_coef_nat_per(f, buf3d, N, d)
1316  integer, intent(in)::N(2), d(2)
1317  sll_real64, dimension(N(1) + 1, N(2)), intent(in)::f
1318  sll_real64, dimension(4, N(1) + 1, N(2) + 1), intent(out)::buf3d
1319  sll_real64, dimension(:), allocatable ::w_left_1, w_left_2
1320  sll_real64 ::tmp
1321  sll_int32::i, j, ii, r_left(2), s_left(2), ind, dd(2)
1322  sll_int32::error
1323  dd(1) = 2*((d(1) + 1)/2)
1324  dd(2) = 2*((d(2) + 1)/2)
1325 
1326  r_left = -dd/2
1327  s_left = (dd + 1)/2
1328 
1329  sll_allocate(w_left_1(-dd(1)/2:dd(1)/2), error)
1330  sll_allocate(w_left_2(-dd(2)/2:dd(2)/2), error)
1331 
1332  call compute_w_hermite(w_left_1, r_left(1), s_left(1))
1333  call compute_w_hermite(w_left_2, r_left(2), s_left(2))
1334 
1335  !print *,'w(',r_left(1),':',s_left(1),')=',w_left_1(r_left(1):s_left(1))
1336  !print *,'w(',r_right(1),':',s_right(1),')=',w_right_1(r_right(1):s_right(1))
1337 
1338  do j = 1, n(2)
1339  do i = 1, n(1) + 1
1340  buf3d(1, i, j) = f(i, j) !f(0,0)
1341  tmp = 0._f64
1342  do ii = r_left(1), s_left(1)
1343  ind = i + ii; if (ind > n(1) + 1) ind = 2*(n(1) + 1) - ind; if (ind < 1) ind = 2 - ind
1344  tmp = tmp + w_left_1(ii)*f(ind, j)
1345  end do
1346  buf3d(2, i, j) = tmp !fx(0,0)
1347  end do
1348  end do
1349  do i = 1, n(1) + 1
1350  do j = 1, n(2)
1351  tmp = 0._f64
1352  do ii = r_left(2), s_left(2)
1353  ind = modulo(j + ii - 1 + n(2), n(2)) + 1
1354  tmp = tmp + w_left_2(ii)*f(i, ind)
1355  end do
1356  buf3d(3, i, j) = tmp !fy(0,0)
1357  end do
1358  end do
1359 
1360  do j = 1, n(2)
1361  do i = 1, n(1) + 1
1362  tmp = 0._f64
1363  do ii = r_left(1), s_left(1)
1364  ind = i + ii; if (ind > n(1) + 1) ind = 2*(n(1) + 1) - ind; if (ind < 1) ind = 2 - ind
1365  tmp = tmp + w_left_1(ii)*buf3d(3, ind, j)
1366  end do
1367  buf3d(4, i, j) = tmp !fxy(0,0)
1368  end do
1369  end do
1370 
1371  buf3d(:, :, n(2) + 1) = buf3d(:, :, 1)
1372 
1373  sll_deallocate_array(w_left_1, error)
1374  sll_deallocate_array(w_left_2, error)
1375 
1376  end subroutine hermite_c1_coef_nat_per
1377 
1378  subroutine compute_w_hermite(w, r, s)
1379  sll_int32, intent(in)::r, s
1380  sll_real64, dimension(r:s), intent(out)::w
1381  sll_int32 ::i, j
1382  sll_real64::tmp
1383 
1384  !maple code for generation of w
1385  !for k from r to -1 do
1386  ! C[k]:=product((k-j),j=r..k-1)*product((k-j),j=k+1..s):
1387  ! C[k]:=1/C[k]*product((-j),j=r..k-1)*product((-j),j=k+1..-1)*product((-j),j=1..s):
1388  !od:
1389  !for k from 1 to s do
1390  ! C[k]:=product((k-j),j=r..k-1)*product((k-j),j=k+1..s):
1391  ! C[k]:=1/C[k]*product((-j),j=r..-1)*product((-j),j=1..k-1)*product((-j),j=k+1..s):
1392  !od:
1393  !C[0]:=-add(C[k],k=r..-1)-add(C[k],k=1..s):
1394 
1395  do i = r, -1
1396  tmp = 1._f64
1397  do j = r, i - 1
1398  tmp = tmp*real(i - j, f64)
1399  end do
1400  do j = i + 1, s
1401  tmp = tmp*real(i - j, f64)
1402  end do
1403  tmp = 1._f64/tmp
1404  do j = r, i - 1
1405  tmp = tmp*real(-j, f64)
1406  end do
1407  do j = i + 1, -1
1408  tmp = tmp*real(-j, f64)
1409  end do
1410  do j = 1, s
1411  tmp = tmp*real(-j, f64)
1412  end do
1413  w(i) = tmp
1414  end do
1415 
1416  do i = 1, s
1417  tmp = 1._f64
1418  do j = r, i - 1
1419  tmp = tmp*real(i - j, f64)
1420  end do
1421  do j = i + 1, s
1422  tmp = tmp*real(i - j, f64)
1423  end do
1424  tmp = 1._f64/tmp
1425  do j = r, -1
1426  tmp = tmp*real(-j, f64)
1427  end do
1428  do j = 1, i - 1
1429  tmp = tmp*real(-j, f64)
1430  end do
1431  do j = i + 1, s
1432  tmp = tmp*real(-j, f64)
1433  end do
1434  w(i) = tmp
1435  end do
1436  tmp = 0._f64
1437  do i = r, -1
1438  tmp = tmp + w(i)
1439  end do
1440  do i = 1, s
1441  tmp = tmp + w(i)
1442  end do
1443  w(0) = -tmp
1444 
1445  !print *,'w',w
1446  !do ii=r,s
1447  ! print *,ii,w(r+s-ii)
1448  !enddo
1449  !
1450 
1451  end subroutine compute_w_hermite
1452 
1453  subroutine localize_polar(x, eta_min, eta_max, ii, eta, N)
1454  sll_real64, intent(in)::x(2), eta_min(2), eta_max(2)
1455  sll_int32, intent(out)::ii(2)
1456  sll_int32, intent(in)::n(2)
1457  sll_real64, intent(out)::eta(2)
1458 
1459  eta(1) = sqrt(x(1)**2 + x(2)**2)
1460  call localize_nat(ii(1), eta(1), eta_min(1), eta_max(1), n(1))
1461  eta(2) = atan2(x(2), x(1))
1462  call localize_per(ii(2), eta(2), eta_min(2), eta_max(2), n(2))
1463  end subroutine localize_polar
1464 
1465  subroutine localize_per(i, x, xmin, xmax, N)
1466  sll_int32, intent(out)::i
1467  sll_real64, intent(inout)::x
1468  sll_real64, intent(in)::xmin, xmax
1469  sll_int32, intent(in)::n
1470  x = (x - xmin)/(xmax - xmin)
1471  x = x - real(floor(x), f64)
1472  x = x*real(n, f64)
1473  i = floor(x)
1474  x = x - real(i, f64)
1475  if (i == n) then
1476  i = 0
1477  x = 0._f64
1478  end if
1479  end subroutine localize_per
1480 
1481  subroutine localize_nat(i, x, xmin, xmax, N)
1482  sll_int32, intent(out)::i
1483  sll_real64, intent(inout)::x
1484  sll_real64, intent(in)::xmin, xmax
1485  sll_int32, intent(in)::n
1486  x = (x - xmin)/(xmax - xmin)
1487  x = x*real(n, f64)
1488  if (x >= real(n, f64)) then
1489  x = real(n, f64)
1490  end if
1491  if (x <= 0._f64) then
1492  x = 0._f64
1493  end if
1494  i = floor(x)
1495  x = x - real(i, f64)
1496  if (i == n) then
1497  i = n - 1
1498  x = 1._f64
1499  end if
1500  end subroutine localize_nat
1501 
1502  subroutine interpolate_hermite(f, i, x, fval, N)
1503  sll_int32, intent(in)::i(2), n(2)
1504  !real(f64),intent(in)::xx(2),xmin(2),xmax(2)
1505  sll_real64, intent(in)::x(2)
1506  sll_real64, intent(out)::fval
1507  sll_real64, dimension(0:8, 0:N(1), 0:N(2))::f
1508  !integer::i(2),i1(2),s
1509  sll_int32::i1(2), s
1510  sll_real64::w(2, 0:3), tmp(0:3)
1511  sll_real64::g(0:3, 0:3)
1512 
1513  !fval =f(0,i(1),i(2))!real(i(1),f64)
1514 
1515  !return
1516 
1517  do s = 1, 2
1518  w(s, 0) = (2._f64*x(s) + 1)*(1._f64 - x(s))*(1._f64 - x(s));
1519  w(s, 1) = x(s)*x(s)*(3._f64 - 2._f64*x(s))
1520  w(s, 2) = x(s)*(1._f64 - x(s))*(1._f64 - x(s))
1521  w(s, 3) = x(s)*x(s)*(x(s) - 1._f64)
1522  i1(s) = i(s) + 1
1523  end do
1524 
1525  g(0, 0) = f(0, i(1), i(2)) !f(0,0)
1526  g(1, 0) = f(0, i1(1), i(2)) !f(1,0)
1527  g(2, 0) = f(1, i(1), i(2)) !fx(0,0)
1528  g(3, 0) = f(2, i(1), i(2)) !fx(1,0)
1529  g(0, 1) = f(0, i(1), i1(2)) !f(0,1)
1530  g(1, 1) = f(0, i1(1), i1(2)) !f(1,1)
1531  g(2, 1) = f(1, i(1), i1(2)) !fx(0,1)
1532  g(3, 1) = f(2, i(1), i1(2)) !fx(1,1)
1533  g(0, 2) = f(3, i(1), i(2)) !fy(0,0)
1534  g(1, 2) = f(3, i1(1), i(2)) !fy(1,0)
1535  g(2, 2) = f(5, i(1), i(2)) !fxy(0,0)
1536  g(3, 2) = f(6, i(1), i(2)) !fxy(1,0)
1537  g(0, 3) = f(4, i(1), i(2)) !fy(0,1)
1538  g(1, 3) = f(4, i1(1), i(2)) !fy(1,1)
1539  g(2, 3) = f(7, i(1), i(2)) !fxy(0,1)
1540  g(3, 3) = f(8, i(1), i(2)) !fxy(1,1)
1541 
1542  do s = 0, 3
1543  tmp(s) = w(1, 0)*g(0, s) + w(1, 1)*g(1, s) + w(1, 2)*g(2, s) + w(1, 3)*g(3, s)
1544  end do
1545 
1546  fval = w(2, 0)*tmp(0) + w(2, 1)*tmp(1) + w(2, 2)*tmp(2) + w(2, 3)*tmp(3)
1547 
1548  !print *,fval,' t',f
1549  end subroutine interpolate_hermite
1550 
1551  subroutine interpolate_hermite_c1(f, i, x, fval, N)
1552  sll_int32, intent(in)::i(2), n(2)
1553  !real(f64),intent(in)::xx(2),xmin(2),xmax(2)
1554  sll_real64, intent(in)::x(2)
1555  sll_real64, intent(out)::fval
1556  sll_real64, dimension(0:3, 0:N(1), 0:N(2))::f
1557  !integer::i(2),i1(2),s
1558  sll_int32::i1(2), s
1559  sll_real64::w(2, 0:3), tmp(0:3)
1560  sll_real64::g(0:3, 0:3)
1561 
1562  !fval =f(0,i(1),i(2))!real(i(1),f64)
1563 
1564  !return
1565 
1566  do s = 1, 2
1567  w(s, 0) = (2._f64*x(s) + 1)*(1._f64 - x(s))*(1._f64 - x(s));
1568  w(s, 1) = x(s)*x(s)*(3._f64 - 2._f64*x(s))
1569  w(s, 2) = x(s)*(1._f64 - x(s))*(1._f64 - x(s))
1570  w(s, 3) = x(s)*x(s)*(x(s) - 1._f64)
1571  i1(s) = i(s) + 1
1572  end do
1573 
1574  g(0, 0) = f(0, i(1), i(2)) !f(0,0)
1575  g(1, 0) = f(0, i1(1), i(2)) !f(1,0)
1576  g(2, 0) = f(1, i(1), i(2)) !fx(0,0)
1577  g(3, 0) = f(1, i1(1), i(2)) !fx(1,0)
1578  g(0, 1) = f(0, i(1), i1(2)) !f(0,1)
1579  g(1, 1) = f(0, i1(1), i1(2)) !f(1,1)
1580  g(2, 1) = f(1, i(1), i1(2)) !fx(0,1)
1581  g(3, 1) = f(1, i1(1), i1(2)) !fx(1,1)
1582  g(0, 2) = f(2, i(1), i(2)) !fy(0,0)
1583  g(1, 2) = f(2, i1(1), i(2)) !fy(1,0)
1584  g(2, 2) = f(3, i(1), i(2)) !fxy(0,0)
1585  g(3, 2) = f(3, i1(1), i(2)) !fxy(1,0)
1586  g(0, 3) = f(2, i(1), i1(2)) !fy(0,1)
1587  g(1, 3) = f(2, i1(1), i1(2)) !fy(1,1)
1588  g(2, 3) = f(3, i(1), i1(2)) !fxy(0,1)
1589  g(3, 3) = f(3, i1(1), i1(2)) !fxy(1,1)
1590 
1591  do s = 0, 3
1592  tmp(s) = w(1, 0)*g(0, s) + w(1, 1)*g(1, s) + w(1, 2)*g(2, s) + w(1, 3)*g(3, s)
1593  end do
1594 
1595  fval = w(2, 0)*tmp(0) + w(2, 1)*tmp(1) + w(2, 2)*tmp(2) + w(2, 3)*tmp(3)
1596  !print *,fval,' t',f
1597  end subroutine interpolate_hermite_c1
1598 
1599  subroutine contribution_hermite_c1(x, val)
1600  sll_real64, intent(in)::x(0:1)
1601  sll_real64, intent(out)::val(4, 0:1, 0:1)
1602  sll_int32::s
1603  sll_real64::w(0:3, 0:1)!,tmp(0:3)
1604  do s = 0, 1
1605  w(0, s) = (2._f64*x(s) + 1)*(1._f64 - x(s))*(1._f64 - x(s));
1606  w(1, s) = x(s)*x(s)*(3._f64 - 2._f64*x(s))
1607  w(2, s) = x(s)*(1._f64 - x(s))*(1._f64 - x(s))
1608  w(3, s) = x(s)*x(s)*(x(s) - 1._f64)
1609  end do
1610 
1611  val(1, 0, 0) = w(0, 0)*w(0, 1)
1612  val(2, 0, 0) = w(2, 0)*w(0, 1)
1613  val(3, 0, 0) = w(0, 0)*w(2, 1)
1614  val(4, 0, 0) = w(2, 0)*w(2, 1)
1615 
1616  val(1, 1, 0) = w(1, 0)*w(0, 1)
1617  val(2, 1, 0) = w(3, 0)*w(0, 1)
1618  val(3, 1, 0) = w(1, 0)*w(2, 1)
1619  val(4, 1, 0) = w(3, 0)*w(2, 1)
1620 
1621  val(1, 0, 1) = w(0, 0)*w(1, 1)
1622  val(2, 0, 1) = w(2, 0)*w(1, 1)
1623  val(3, 0, 1) = w(0, 0)*w(3, 1)
1624  val(4, 0, 1) = w(2, 0)*w(3, 1)
1625 
1626  val(1, 1, 1) = w(1, 0)*w(1, 1)
1627  val(2, 1, 1) = w(3, 0)*w(1, 1)
1628  val(3, 1, 1) = w(1, 0)*w(3, 1)
1629  val(4, 1, 1) = w(3, 0)*w(3, 1)
1630 
1631  !print *,val
1632  !stop
1633 
1634  end subroutine contribution_hermite_c1
1635 
1636  subroutine contribution_spl(x, val)
1637  sll_real64, intent(in)::x(0:1)
1638  sll_real64, intent(out)::val(-1:2, -1:2)
1639  sll_int32::s
1640  sll_real64::w(-1:2, 0:1)
1641  do s = 0, 1
1642  w(-1, s) = (1._f64/6._f64)*(1._f64 - x(s))*(1._f64 - x(s))*(1._f64 - x(s));
1643  w(0, s) = 1._f64/6._f64 + 0.5_f64*(1._f64 - x(s))*(-(1._f64 - x(s))* &
1644  (1._f64 - x(s)) + (1._f64 - x(s)) + 1._f64);
1645  w(1, s) = 1._f64/6._f64 + 0.5_f64*x(s)*(-x(s)*x(s) + x(s) + 1._f64);
1646  w(2, s) = (1._f64/6._f64)*x(s)*x(s)*x(s);
1647  end do
1648 
1649  val(-1, -1) = w(-1, 0)*w(-1, 1)
1650  val(-1, 0) = w(-1, 0)*w(0, 1)
1651  val(-1, 1) = w(-1, 0)*w(1, 1)
1652  val(-1, 2) = w(-1, 0)*w(2, 1)
1653 
1654  val(0, -1) = w(0, 0)*w(-1, 1)
1655  val(0, 0) = w(0, 0)*w(0, 1)
1656  val(0, 1) = w(0, 0)*w(1, 1)
1657  val(0, 2) = w(0, 0)*w(2, 1)
1658 
1659  val(1, -1) = w(1, 0)*w(-1, 1)
1660  val(1, 0) = w(1, 0)*w(0, 1)
1661  val(1, 1) = w(1, 0)*w(1, 1)
1662  val(1, 2) = w(1, 0)*w(2, 1)
1663 
1664  val(2, -1) = w(2, 0)*w(-1, 1)
1665  val(2, 0) = w(2, 0)*w(0, 1)
1666  val(2, 1) = w(2, 0)*w(1, 1)
1667  val(2, 2) = w(2, 0)*w(2, 1)
1668 
1669  end subroutine contribution_spl
1670 
1671  subroutine splcoefnat1d0(lunat, N)
1672  sll_int32, intent(in)::n
1673  sll_real64, dimension(0:1, 0:N + 2), intent(out)::lunat
1674  sll_int32::i
1675  sll_real64::a(0:2)
1676  a(0) = -3._f64; a(1) = 0._f64; a(2) = 3._f64;
1677  lunat(0, 0) = a(0);
1678  lunat(1, 0) = 1._f64/lunat(0, 0);
1679  lunat(0, 1) = 4._f64 - a(1)*lunat(1, 0);
1680  lunat(1, 1) = 1._f64/lunat(0, 1);
1681  lunat(0, 2) = 4._f64 - lunat(1, 1)*(1._f64 - a(2)/a(0));
1682  do i = 2, n
1683  lunat(1, i) = 1._f64/lunat(0, i);
1684  lunat(0, i + 1) = 4._f64 - lunat(1, i);
1685  end do
1686  lunat(1, n + 2) = a(0)/lunat(0, n);
1687  lunat(1, n + 1) = (a(1) - lunat(1, n + 2))/lunat(0, n + 1);
1688  lunat(0, n + 2) = a(2) - lunat(1, n + 1);
1689  end subroutine splcoefnat1d0
1690 
1691  subroutine splcoefnat1d(p, lunat, N)
1692  sll_int32, intent(in)::n
1693  sll_real64, dimension(0:1, 0:N + 2), intent(in)::lunat
1694  sll_real64, dimension(0:N + 2), intent(inout)::p
1695  sll_int32::i
1696  sll_real64::a(0:2)
1697  a(0) = -3._f64; a(1) = 0._f64; a(2) = 3._f64;
1698  !p(0)=0._f64;
1699  !p(N+2)=0._f64;
1700  do i = 0, n + 2;
1701  p(i) = 6._f64*p(i);
1702  end do;
1703  do i = 1, n + 1; p(i) = p(i) - lunat(1, i - 1)*p(i - 1); end do
1704  p(n + 2) = p(n + 2) - (lunat(1, n + 1)*p(n + 1) + lunat(1, n + 2)*p(n));
1705  p(n + 2) = p(n + 2)/lunat(0, n + 2);
1706  do i = n + 1, 2, -1; p(i) = (p(i) - p(i + 1))/lunat(0, i); end do
1707  p(1) = (p(1) - (1._f64 - a(2)/a(0))*p(2))/lunat(0, 1);
1708  p(0) = (p(0) - a(1)*p(1) - a(2)*p(2))/lunat(0, 0);
1709  !p(i-1)+4*p(i)+p(i+1)=6ftab(i-1,j), i=1..N+1
1710  !a(0)*p(0)+a(1)*p(1)+a(2)*p(2)=f'(rmin)=0);
1711  !a(0)*p(Na)+a(1)*p(Na+1)+a(2)*p(Na+2)=f'(rmax)=0);
1712  end subroutine splcoefnat1d
1713 
1714  subroutine splcoefnat1d0old(dnat, lnat, N)
1715  sll_int32, intent(in)::n
1716  sll_real64, dimension(0:N + 2), intent(inout)::dnat, lnat
1717  sll_int32::i
1718  sll_real64::a(0:2)
1719  a(0) = -3._f64; a(1) = 0._f64; a(2) = 3._f64;
1720  dnat(0) = a(0);
1721  lnat(0) = 1._f64/dnat(0);
1722  dnat(1) = 4._f64 - a(1)*lnat(0);
1723  lnat(1) = 1._f64/dnat(1);
1724  dnat(2) = 4._f64 - lnat(1)*(1._f64 - a(2)/a(0));
1725  do i = 2, n
1726  lnat(i) = 1._f64/dnat(i);
1727  dnat(i + 1) = 4._f64 - lnat(i);
1728  end do
1729  lnat(n + 2) = a(0)/dnat(n);
1730  lnat(n + 1) = (a(1) - lnat(n + 2))/dnat(n + 1);
1731  dnat(n + 2) = a(2) - lnat(n + 1);
1732  end subroutine splcoefnat1d0old
1733 
1734  subroutine splcoefnatper2d(f, buf, dnatx, lnatx, dpery, lpery, mpery, Nx, Ny)
1735  sll_int32, intent(in)::nx, ny
1736  sll_real64, dimension(:, :), pointer::f
1737  sll_real64, dimension(:), pointer::buf, dpery, lpery, mpery, dnatx, lnatx
1738  sll_int32::i, j
1739  !natural spline coefficients in x
1740  do j = 0, ny - 1
1741  do i = 0, nx + 2; buf(i) = f(i, j); end do
1742  call splcoefnat1dold(buf, dnatx, lnatx, nx)
1743  do i = 0, nx + 2; f(i, j) = buf(i); end do
1744  end do
1745  !periodic spline coefficients in y
1746  do i = 0, nx + 2
1747  do j = 0, ny - 1; buf(j) = f(i, j); end do
1748  call splcoefper1dold(buf, dpery, lpery, mpery, ny)
1749  do j = 0, ny - 1; f(i, j) = buf(j); end do
1750  end do
1751 
1752  end subroutine splcoefnatper2d
1753 
1754  subroutine splnatper2d(f, xx, xmin, xmax, yy, ymin, ymax, fval, Nx, Ny)
1755  sll_int32, intent(in)::nx, ny
1756  sll_real64, intent(in)::xx, xmin, xmax, yy, ymin, ymax
1757  sll_real64, intent(out)::fval
1758  sll_real64, dimension(:, :), pointer::f
1759  sll_int32::i, j
1760  sll_int32::ix(0:3), iy(0:3)
1761  sll_real64::x, y
1762  sll_real64::wx(0:3), wy(0:3), tmp(0:3)
1763 
1764  x = (xx - xmin)/(xmax - xmin)
1765  !x=x-real(floor(x),f64)
1766  if (x < 0._f64) x = 0._f64;
1767  if (x >= 1.0_f64) then
1768  x = 1.0_f64 - 1.e-12_f64;
1769  end if
1770  x = x*real(nx, f64)
1771  i = floor(x)
1772  x = x - real(i, f64)
1773 
1774  y = (yy - ymin)/(ymax - ymin)
1775  y = y - real(floor(y), f64)
1776  y = y*real(ny, f64)
1777  j = floor(y)
1778  y = y - real(j, f64)
1779 
1780  wx(0) = (1.0_f64/6._f64)*(1.0_f64 - x)*(1.0_f64 - x)*(1.0_f64 - x);
1781  wx(1) = &
1782  1.0_f64/6._f64 + &
1783  0.5_f64*(1.0_f64 - x)*(-(1.0_f64 - x)*(1.0_f64 - x) + (1.0_f64 - x) + 1.0_f64);
1784  wx(2) = 1.0_f64/6._f64 + 0.5_f64*x*(-x*x + x + 1.0_f64);
1785  wx(3) = (1.0_f64/6._f64)*x*x*x;
1786  wy(0) = (1.0_f64/6._f64)*(1.0_f64 - y)*(1.0_f64 - y)*(1.0_f64 - y);
1787  wy(1) = &
1788  1.0_f64/6._f64 + &
1789  0.5_f64*(1.0_f64 - y)*(-(1.0_f64 - y)*(1.0_f64 - y) + (1.0_f64 - y) + 1.0_f64);
1790  wy(2) = 1.0_f64/6._f64 + 0.5_f64*y*(-y*y + y + 1.0_f64);
1791  wy(3) = (1.0_f64/6._f64)*y*y*y;
1792  iy(0) = mod(j + ny - 1, ny)
1793  iy(1) = j
1794  iy(2) = mod(j + 1, ny)
1795  iy(3) = mod(j + 2, ny)
1796 
1797  ix(0) = i
1798  ix(1) = i + 1
1799  ix(2) = i + 2
1800  ix(3) = i + 3
1801 
1802  tmp(0) = wx(0)*f(ix(0), iy(0)) + wx(1)*f(ix(1), iy(0)) &
1803  + wx(2)*f(ix(2), iy(0)) + wx(3)*f(ix(3), iy(0))
1804  tmp(1) = wx(0)*f(ix(0), iy(1)) + wx(1)*f(ix(1), iy(1)) &
1805  + wx(2)*f(ix(2), iy(1)) + wx(3)*f(ix(3), iy(1))
1806  tmp(2) = wx(0)*f(ix(0), iy(2)) + wx(1)*f(ix(1), iy(2)) &
1807  + wx(2)*f(ix(2), iy(2)) + wx(3)*f(ix(3), iy(2))
1808  tmp(3) = wx(0)*f(ix(0), iy(3)) + wx(1)*f(ix(1), iy(3)) &
1809  + wx(2)*f(ix(2), iy(3)) + wx(3)*f(ix(3), iy(3))
1810 
1811  fval = wy(0)*tmp(0) + wy(1)*tmp(1) + wy(2)*tmp(2) + wy(3)*tmp(3)
1812 
1813  end subroutine splnatper2d
1814 
1815  subroutine splper1d(f, xx, xmin, xmax, fval, N)
1816  sll_int32, intent(in)::n
1817  sll_real64, intent(in)::xx, xmin, xmax
1818  sll_real64, intent(out)::fval
1819  sll_real64, dimension(0:N - 1), intent(inout)::f
1820  sll_int32::i
1821  real(f64)::x
1822  real(f64)::w(0:3)
1823  x = (xx - xmin)/(xmax - xmin)
1824  x = x - real(floor(x), f64)
1825  x = x*real(n, f64)
1826  i = floor(x)
1827  x = x - real(i, f64)
1828  w(0) = (1.0_f64/6._f64)*(1.0_f64 - x)*(1.0_f64 - x)*(1.0_f64 - x);
1829  w(1) = &
1830  1.0_f64/6._f64 + &
1831  0.5_f64*(1.0_f64 - x)*(-(1.0_f64 - x)*(1.0_f64 - x) + (1.0_f64 - x) + 1.0_f64);
1832  w(2) = 1.0_f64/6._f64 + 0.5_f64*x*(-x*x + x + 1.0_f64);
1833  w(3) = (1.0_f64/6._f64)*x*x*x;
1834  fval = w(0)*f(mod(i + n - 1, n)) + w(1)*f(i) + w(2)*f(mod(i + 1, n)) + w(3)*f(mod(i + 2, n))
1835  end subroutine splper1d
1836 
1837  subroutine splnat1d(f, xx, xmin, xmax, fval, N)
1838  sll_int32, intent(in)::n
1839  sll_real64, intent(in)::xx, xmin, xmax
1840  sll_real64, intent(out)::fval
1841  sll_real64, dimension(0:N + 2), intent(inout)::f
1842  sll_int32::i
1843  sll_real64::x
1844  sll_real64::w(0:3)
1845  x = (xx - xmin)/(xmax - xmin)
1846  !x=x-real(floor(x),f64)
1847  if (x < 0._f64) x = 0._f64;
1848  if (x > 1.0_f64) then
1849  x = 1.0_f64 - 1.e-12_f64;
1850  end if
1851  x = x*real(n, f64)
1852  i = floor(x)
1853  x = x - real(i, f64)
1854 
1855  w(0) = (1.0_f64/6._f64)*(1.0_f64 - x)*(1.0_f64 - x)*(1.0_f64 - x);
1856  w(1) = &
1857  1.0_f64/6.0_f64 + &
1858  0.5_f64*(1.0_f64 - x)*(-(1.0_f64 - x)*(1.0_f64 - x) + (1.0_f64 - x) + 1.0_f64);
1859  w(2) = 1.0_f64/6._f64 + 0.5_f64*x*(-x*x + x + 1.0_f64);
1860  w(3) = (1.0_f64/6._f64)*x*x*x;
1861  fval = w(0)*f(i) + w(1)*f(i + 1) + w(2)*f(i + 2) + w(3)*f(i + 3)
1862  end subroutine splnat1d
1863 
1864  subroutine splcoefper1d0old(dper, lper, mper, N)
1865  sll_int32, intent(in)::n
1866  sll_real64, dimension(0:N - 1), intent(out)::dper, lper, mper
1867  sll_int32::i
1868 
1869  dper(0) = 4._f64
1870  mper(0) = 0.25_f64
1871  do i = 0, n - 2
1872  lper(i) = 1.0_f64/dper(i)
1873  dper(i + 1) = 4._f64 - lper(i)
1874  mper(i + 1) = -mper(i)/dper(i + 1)
1875  end do
1876  dper(n - 1) = dper(n - 1) - (lper(n - 2) + 2._f64*mper(n - 2))
1877  do i = 0, n - 1
1878  dper(i) = 1.0_f64/dper(i)
1879  end do
1880  end subroutine splcoefper1d0old
1881 
1882  subroutine splcoefper1d0(luper, N)
1883  sll_int32, intent(in)::n
1884  sll_real64, dimension(0:3*N - 1), intent(out)::luper
1885  sll_int32::i
1886 
1887  luper(0 + 3*0) = 4._f64
1888  luper(2 + 3*0) = 0.25_f64
1889  do i = 0, n - 2
1890  luper(1 + 3*i) = 1.0_f64/luper(0 + 3*i)
1891  luper(0 + 3*(i + 1)) = 4._f64 - luper(1 + 3*i)
1892  luper(2 + 3*(i + 1)) = -luper(2 + 3*i)/luper(0 + 3*(i + 1))
1893  end do
1894  luper(0 + 3*(n - 1)) = &
1895  luper(0 + 3*(n - 1)) - (luper(1 + 3*(n - 2)) + 2._f64*luper(2 + 3*(n - 2)))
1896  do i = 0, n - 1
1897  luper(0 + 3*i) = 1.0_f64/luper(0 + 3*i)
1898  end do
1899  end subroutine splcoefper1d0
1900 
1901  subroutine splcoefper1d(f, luper, N)
1902  sll_int32, intent(in)::n
1903  sll_real64, dimension(0:3*N - 1), intent(in)::luper
1904  sll_real64, dimension(0:N - 1), intent(inout)::f
1905  sll_int32::i
1906  do i = 0, n - 1; f(i) = 6._f64*f(i); end do;
1907  do i = 1, n - 1
1908  f(i) = f(i) - f(i - 1)*luper(1 + 3*(i - 1))
1909  end do
1910  do i = 0, n - 2
1911  f(n - 1) = f(n - 1) - luper(2 + 3*i)*f(i)
1912  end do
1913  f(n - 1) = f(n - 1)*luper(0 + 3*(n - 1)); f(n - 2) = &
1914  luper(0 + 3*(n - 2))*(f(n - 2) - (1.0_f64 - luper(2 + 3*(n - 3)))*f(n - 1))
1915  do i = n - 3, 1, -1
1916  f(i) = luper(0 + 3*i)*(f(i) - f(i + 1) + luper(2 + 3*(i - 1))*f(n - 1))
1917  end do
1918  f(0) = luper(0 + 3*0)*(f(0) - f(1) - f(n - 1));
1919  end subroutine splcoefper1d
1920 
1921  subroutine splcoefnat1dold(p, dnat, lnat, N)
1922  sll_int32, intent(in)::n
1923  sll_real64, dimension(0:N + 2), intent(in)::dnat, lnat
1924  sll_real64, dimension(0:N + 2), intent(inout)::p
1925  sll_int32::i
1926  sll_real64::a(0:2)
1927  a(0) = -3._f64; a(1) = 0._f64; a(2) = 3._f64;
1928  !p(0)=0._f64;
1929  !p(N+2)=0._f64;
1930  do i = 0, n + 2;
1931  p(i) = 6._f64*p(i);
1932  end do;
1933  do i = 1, n + 1; p(i) = p(i) - lnat(i - 1)*p(i - 1); end do
1934  p(n + 2) = p(n + 2) - (lnat(n + 1)*p(n + 1) + lnat(n + 2)*p(n));
1935  p(n + 2) = p(n + 2)/dnat(n + 2);
1936  do i = n + 1, 2, -1; p(i) = (p(i) - p(i + 1))/dnat(i); end do
1937  p(1) = (p(1) - (1.0_f64 - a(2)/a(0))*p(2))/dnat(1);
1938  p(0) = (p(0) - a(1)*p(1) - a(2)*p(2))/dnat(0);
1939  !p(i-1)+4*p(i)+p(i+1)=6ftab(i-1,j), i=1..N+1
1940  !a(0)*p(0)+a(1)*p(1)+a(2)*p(2)=f'(rmin)=0);
1941  !a(0)*p(Na)+a(1)*p(Na+1)+a(2)*p(Na+2)=f'(rmax)=0);
1942  end subroutine splcoefnat1dold
1943 
1944  subroutine splcoefper1dold(f, dper, lper, mper, N)
1945  sll_int32, intent(in)::n
1946  sll_real64, dimension(0:N - 1), intent(in)::dper, lper, mper
1947  sll_real64, dimension(0:N - 1), intent(inout)::f
1948  sll_int32::i
1949  do i = 0, n - 1; f(i) = 6._f64*f(i); end do;
1950  do i = 1, n - 1
1951  f(i) = f(i) - f(i - 1)*lper(i - 1)
1952  end do
1953  do i = 0, n - 2
1954  f(n - 1) = f(n - 1) - mper(i)*f(i)
1955  end do
1956  f(n - 1) = f(n - 1)*dper(n - 1); f(n - 2) = dper(n - 2)*(f(n - 2) - (1.0_f64 - mper(n - 3))*f(n - 1))
1957  do i = n - 3, 1, -1
1958  f(i) = dper(i)*(f(i) - f(i + 1) + mper(i - 1)*f(n - 1))
1959  end do
1960  f(0) = dper(0)*(f(0) - f(1) - f(n - 1));
1961  end subroutine splcoefper1dold
1962 
1963  subroutine solve_tridiag(a, b, c, v, x, n)
1964 ! Using Thomas' Algorithm
1965  implicit none
1966  ! a - sub-diagonal (means it is the diagonal below the main diagonal)
1967  ! b - the main diagonal
1968  ! c - sup-diagonal (means it is the diagonal above the main diagonal)
1969  ! v - right part
1970  ! x - the answer
1971  ! n - number of equations
1972 
1973  sll_int32, intent(in) :: n
1974  sll_real64, dimension(n), intent(in) :: a, b, c, v
1975  sll_real64, dimension(n), intent(out) :: x
1976  sll_real64, dimension(n) :: bp, vp
1977  sll_real64 :: m
1978  sll_int32 i
1979 
1980  ! Make copies of the b and v variables so that they are unaltered by this sub
1981  bp(1) = b(1)
1982  vp(1) = v(1)
1983 
1984  !The first pass (setting coefficients):
1985  firstpass: do i = 2, n
1986  m = a(i)/bp(i - 1)
1987  bp(i) = b(i) - m*c(i - 1)
1988  vp(i) = v(i) - m*vp(i - 1)
1989  end do firstpass
1990 
1991  x(n) = vp(n)/bp(n)
1992  !The second pass (back-substition)
1993  backsub: do i = n - 1, 1, -1
1994  x(i) = (vp(i) - c(i)*x(i + 1))/bp(i)
1995  end do backsub
1996 
1997  end subroutine solve_tridiag
1998 
1999  subroutine sll_s_penta(n, e, a, b, c, f, d, x)
2000 !************************** P E N T A ***************************
2001 ! *
2002 ! Solution of a linear system of algebraic equations with *
2003 ! a pentadiagonal matrix of coefficients.(No pivoting) *
2004 ! Equation no. i : *
2005 ! *
2006 ! e(i)*x(i-2)+ a(i)*x(i-1) + b(i)*x(i) + *
2007 ! c(i)*x(i+1) + f(i)*x(i+2) = d(i), i = 1,2,...n *
2008 ! *
2009 ! === Use === *
2010 ! *
2011 ! call sll_s_penta(n,a,b,c,d,x) *
2012 ! or *
2013 ! call sll_s_penta(n,a,b,c,d,d) *
2014 ! *
2015 ! In the last case, vector d contains the solution. *
2016 ! *
2017 ! === Input === *
2018 ! *
2019 ! n ....... integer . Number of equations *
2020 ! e(1:n) .. real vector . Lowest diagonal.Elements e(1) *
2021 ! and e(2) are not used *
2022 ! a(1:n) .. real vector . Lower diagonal.Element a(1) *
2023 ! is not used. *
2024 ! b(1:n) .. real vector . Main diagonal *
2025 ! c(1:n) .. real vector . Upper diagonal.Element c(n) *
2026 ! is not used. *
2027 ! f(1:n) .. real vector . Uppermost diagonal.Elements *
2028 ! f(n-1)and f(n) are not used. *
2029 ! d(1:n) .. real vector . Right hand side of the system. *
2030 ! *
2031 ! === Output === *
2032 ! *
2033 ! x(1:n) .. real vector . The solution vector *
2034 ! *
2035 !********************** fortran 90 ******************************
2036  implicit none
2037  sll_int32, intent(in) :: n
2038  sll_real64, intent(inout), dimension(n) :: e, a, b, c, f, d
2039  sll_real64, intent(out) :: x(n)
2040  ! --- Local variables ---
2041  sll_int32 :: i
2042  sll_real64 :: q
2043  do i = 2, n - 1
2044  q = a(i)/b(i - 1)
2045  b(i) = b(i) - q*c(i - 1)
2046  c(i) = c(i) - q*f(i - 1)
2047  d(i) = d(i) - q*d(i - 1)
2048  q = e(i + 1)/b(i - 1)
2049  a(i + 1) = a(i + 1) - q*c(i - 1)
2050  b(i + 1) = b(i + 1) - q*f(i - 1)
2051  d(i + 1) = d(i + 1) - q*d(i - 1)
2052  end do
2053  q = a(n)/b(n - 1)
2054  b(n) = b(n) - q*c(n - 1)
2055  x(n) = (d(n) - q*d(n - 1))/b(n)
2056  x(n - 1) = (d(n - 1) - c(n - 1)*x(n))/b(n - 1)
2057  do i = n - 2, 1, -1
2058  x(i) = (d(i) - f(i)*x(i + 2) - c(i)*x(i + 1))/b(i)
2059  end do
2060  return
2061  end subroutine sll_s_penta
2062 
2063 !subroutine sll_s_compute_shape_circle(points,N_points)
2064 ! sll_int32,intent(in) :: N_points
2065 ! sll_real64,dimension(:,:) ::points
2066 ! sll_int32 :: i
2067 ! sll_real64 :: x
2068 ! do i=1,N_points
2069 ! x = 2._f64*sll_p_pi*real(i,f64)/(real(N_points,f64))
2070 ! points(1,i) = cos(x)
2071 ! points(2,i) = sin(x)
2072 ! points(3,i) = 1.0_f64/real(N_points,f64)
2073 ! enddo
2074 !
2075 !end subroutine sll_s_compute_shape_circle
2076 !
2077 
2078 ! subroutine sll_s_compute_init_f_polar(f,mode,N,eta_min,eta_max)
2079 ! sll_real64,dimension(:,:),intent(out)::f
2080 ! sll_int32,intent(in)::N(2),mode(2)
2081 ! sll_real64,intent(in)::eta_min(2),eta_max(2)
2082 ! sll_int32::i,j
2083 ! sll_real64::eta(2),delta_eta(2),kmode,val
2084 !
2085 ! call sll_s_zero_bessel_dir_dir(mode,eta_min(1),eta_max(1),val)
2086 ! delta_eta(1)=(eta_max(1)-eta_min(1))/real(N(1),f64)
2087 ! delta_eta(2)=(eta_max(2)-eta_min(2))/real(N(2),f64)
2088 !
2089 ! kmode=real(mode(2),f64)*(2._f64*sll_p_pi)/(eta_max(2)-eta_min(2))
2090 !
2091 ! do j=1,N(2)+1
2092 ! eta(2)=eta_min(2)+real(j-1,f64)*delta_eta(2)
2093 ! do i=1,N(1)+1
2094 ! eta(1)=eta_min(1)+real(i-1,f64)*delta_eta(1)
2095 ! eta(1)=val*eta(1)/eta_max(1)
2096 ! f(i,j)= 0._f64 !temporary, because DBESJ not recognized on helios
2097 ! f(i,j) = (DBESJN(mode(2),val)*DBESYN(mode(2),eta(1))-DBESYN(mode(2),val)*DBESJN(mode(2),eta(1)))*cos(kmode*eta(2))
2098 ! enddo
2099 ! enddo
2100 !
2101 ! end subroutine sll_s_compute_init_f_polar
2102 !
2103 !
2104 ! subroutine sll_s_zero_bessel_dir_dir(mode,eta_min,eta_max,val)
2105 ! sll_real64,intent(in)::eta_min,eta_max
2106 ! sll_int32,intent(in)::mode(2)
2107 ! sll_real64,intent(out)::val
2108 ! sll_real64::alpha,tmp
2109 ! sll_int32::mode_max(2),i,j
2110 ! logical::is_file
2111 ! val=0._f64
2112 ! INQUIRE(FILE="zeros_bessel.txt", EXIST=is_file)
2113 ! if((is_file).eqv.(.false.))then
2114 ! print *,'#file zeros_bessel.txt does not exist'
2115 ! return
2116 ! endif
2117 ! open(27,file='zeros_bessel.txt',action="read")
2118 ! read(27,*) mode_max(1),mode_max(2),alpha
2119 ! close(27)
2120 ! if((mode(1)<1).or.(mode(1)>mode_max(1)))then
2121 ! print *,'#bad value of mode(1) vs mode_max(1)',mode(1),mode_max(1)
2122 ! return
2123 ! endif
2124 ! if((mode(2)<0).or.(mode(2)>mode_max(2)))then
2125 ! print *,'#bad value of mode(2) vs mode_max(2)',mode(2),mode_max(2)
2126 ! return
2127 ! endif
2128 ! if(abs(alpha-eta_min/eta_max)>1.e-12)then
2129 ! print *,'#bad value of rmin/rmax w.r.t zeros_bessel.txt',eta_min/eta_max,alpha
2130 ! return
2131 ! endif
2132 ! open(27,file='zeros_bessel.txt',action="read")
2133 ! read(27,*) mode_max(1),mode_max(2),alpha
2134 ! read(27,*) i,j,tmp
2135 ! do while((i.ne.mode(1)).or.(j.ne.mode(2)))
2136 ! read(27,*) i,j,tmp
2137 ! enddo
2138 ! close(27)
2139 ! val = tmp
2140 !
2141 ! end subroutine sll_s_zero_bessel_dir_dir
2142 !
2143 
2144 end module sll_m_gyroaverage_2d_polar
sll_m_constants
Fortran module where set some physical and mathematical constants.
Definition: sll_m_constants.F90:23
sll_m_constants::sll_p_pi
real(kind=f64), parameter, public sll_p_pi
Definition: sll_m_constants.F90:51
sll_m_gyroaverage_2d_polar
Definition: sll_m_gyroaverage_2d_polar.F90:18
sll_m_gyroaverage_2d_polar::sll_f_new_plan_gyroaverage_polar_splines
type(sll_t_plan_gyroaverage_polar) function, pointer, public sll_f_new_plan_gyroaverage_polar_splines(eta_min, eta_max, Nc, N_points)
Definition: sll_m_gyroaverage_2d_polar.F90:111
sll_m_gyroaverage_2d_polar::splcoefnat1d0
subroutine splcoefnat1d0(lunat, N)
Definition: sll_m_gyroaverage_2d_polar.F90:1672
sll_m_gyroaverage_2d_polar::splcoefper1d
subroutine splcoefper1d(f, luper, N)
Definition: sll_m_gyroaverage_2d_polar.F90:1902
sll_m_gyroaverage_2d_polar::sll_s_compute_gyroaverage_pade_polar
subroutine, public sll_s_compute_gyroaverage_pade_polar(gyro, f, rho)
Definition: sll_m_gyroaverage_2d_polar.F90:1001
sll_m_gyroaverage_2d_polar::sll_s_compute_gyroaverage_points_polar_hermite
subroutine, public sll_s_compute_gyroaverage_points_polar_hermite(gyro, f, rho)
Definition: sll_m_gyroaverage_2d_polar.F90:153
sll_m_gyroaverage_2d_polar::splcoefper1d0old
subroutine splcoefper1d0old(dper, lper, mper, N)
Definition: sll_m_gyroaverage_2d_polar.F90:1865
sll_m_gyroaverage_2d_polar::sll_s_compute_gyroaverage_points_polar_spl
subroutine, public sll_s_compute_gyroaverage_points_polar_spl(gyro, f, rho)
Definition: sll_m_gyroaverage_2d_polar.F90:505
sll_m_gyroaverage_2d_polar::splcoefnat1dold
subroutine splcoefnat1dold(p, dnat, lnat, N)
Definition: sll_m_gyroaverage_2d_polar.F90:1922
sll_m_gyroaverage_2d_polar::localize_polar
subroutine localize_polar(x, eta_min, eta_max, ii, eta, N)
Definition: sll_m_gyroaverage_2d_polar.F90:1454
sll_m_gyroaverage_2d_polar::sll_s_pre_compute_gyroaverage_polar_spl
subroutine, public sll_s_pre_compute_gyroaverage_polar_spl(gyro, rho)
Definition: sll_m_gyroaverage_2d_polar.F90:663
sll_m_gyroaverage_2d_polar::splcoefper1dold
subroutine splcoefper1dold(f, dper, lper, mper, N)
Definition: sll_m_gyroaverage_2d_polar.F90:1945
sll_m_gyroaverage_2d_polar::sll_s_pre_compute_gyroaverage_polar_spl_fft
subroutine, public sll_s_pre_compute_gyroaverage_polar_spl_fft(gyro, rho)
Definition: sll_m_gyroaverage_2d_polar.F90:825
sll_m_gyroaverage_2d_polar::compute_w_hermite
subroutine compute_w_hermite(w, r, s)
Definition: sll_m_gyroaverage_2d_polar.F90:1379
sll_m_gyroaverage_2d_polar::solve_tridiag
subroutine solve_tridiag(a, b, c, v, x, n)
Definition: sll_m_gyroaverage_2d_polar.F90:1964
sll_m_gyroaverage_2d_polar::sll_s_compute_gyroaverage_points_polar_hermite_c1
subroutine, public sll_s_compute_gyroaverage_points_polar_hermite_c1(gyro, f, rho)
Definition: sll_m_gyroaverage_2d_polar.F90:190
sll_m_gyroaverage_2d_polar::sll_s_compute_gyroaverage_pre_compute_polar_spl
subroutine, public sll_s_compute_gyroaverage_pre_compute_polar_spl(gyro, f)
Definition: sll_m_gyroaverage_2d_polar.F90:749
sll_m_gyroaverage_2d_polar::localize_nat
subroutine localize_nat(i, x, xmin, xmax, N)
Definition: sll_m_gyroaverage_2d_polar.F90:1482
sll_m_gyroaverage_2d_polar::splper1d
subroutine splper1d(f, xx, xmin, xmax, fval, N)
Definition: sll_m_gyroaverage_2d_polar.F90:1816
sll_m_gyroaverage_2d_polar::sll_f_new_plan_gyroaverage_polar_pade
type(sll_t_plan_gyroaverage_polar) function, pointer, public sll_f_new_plan_gyroaverage_polar_pade(eta_min, eta_max, Nc)
Definition: sll_m_gyroaverage_2d_polar.F90:135
sll_m_gyroaverage_2d_polar::sll_s_compute_gyroaverage_pade_high_order_polar
subroutine, public sll_s_compute_gyroaverage_pade_high_order_polar(gyro, f, rho, order)
Definition: sll_m_gyroaverage_2d_polar.F90:1065
sll_m_gyroaverage_2d_polar::splcoefnat1d
subroutine splcoefnat1d(p, lunat, N)
Definition: sll_m_gyroaverage_2d_polar.F90:1692
sll_m_gyroaverage_2d_polar::splcoefnat1d0old
subroutine splcoefnat1d0old(dnat, lnat, N)
Definition: sll_m_gyroaverage_2d_polar.F90:1715
sll_m_gyroaverage_2d_polar::interpolate_hermite
subroutine interpolate_hermite(f, i, x, fval, N)
Definition: sll_m_gyroaverage_2d_polar.F90:1503
sll_m_gyroaverage_2d_polar::splcoefnatper2d
subroutine splcoefnatper2d(f, buf, dnatx, lnatx, dpery, lpery, mpery, Nx, Ny)
Definition: sll_m_gyroaverage_2d_polar.F90:1735
sll_m_gyroaverage_2d_polar::sll_f_new_plan_gyroaverage_polar_hermite
type(sll_t_plan_gyroaverage_polar) function, pointer, public sll_f_new_plan_gyroaverage_polar_hermite(eta_min, eta_max, Nc, N_points, interp_degree, deriv_size)
Definition: sll_m_gyroaverage_2d_polar.F90:83
sll_m_gyroaverage_2d_polar::hermite_coef_nat_per
subroutine hermite_coef_nat_per(f, buf3d, N, d)
Definition: sll_m_gyroaverage_2d_polar.F90:1219
sll_m_gyroaverage_2d_polar::hermite_c1_coef_nat_per
subroutine hermite_c1_coef_nat_per(f, buf3d, N, d)
Definition: sll_m_gyroaverage_2d_polar.F90:1316
sll_m_gyroaverage_2d_polar::sll_s_compute_gyroaverage_points_polar_with_invar_spl
subroutine, public sll_s_compute_gyroaverage_points_polar_with_invar_spl(gyro, f, rho)
Definition: sll_m_gyroaverage_2d_polar.F90:586
sll_m_gyroaverage_2d_polar::assemble_csr_from_pre_compute
subroutine assemble_csr_from_pre_compute(pre_compute_coeff, pre_compute_N, pre_compute_index, nb, Nc, ia, ja, a, num_rows, num_nz)
Definition: sll_m_gyroaverage_2d_polar.F90:399
sll_m_gyroaverage_2d_polar::sll_s_compute_gyroaverage_pre_compute_polar_hermite_c1
subroutine, public sll_s_compute_gyroaverage_pre_compute_polar_hermite_c1(gyro, f)
Definition: sll_m_gyroaverage_2d_polar.F90:462
sll_m_gyroaverage_2d_polar::sll_s_pre_compute_gyroaverage_polar_hermite_c1
subroutine, public sll_s_pre_compute_gyroaverage_polar_hermite_c1(gyro, rho)
Definition: sll_m_gyroaverage_2d_polar.F90:279
sll_m_gyroaverage_2d_polar::sll_s_compute_gyroaverage_points_polar_with_invar_hermite_c1
subroutine, public sll_s_compute_gyroaverage_points_polar_with_invar_hermite_c1(gyro, f, rho)
Definition: sll_m_gyroaverage_2d_polar.F90:233
sll_m_gyroaverage_2d_polar::contribution_hermite_c1
subroutine contribution_hermite_c1(x, val)
Definition: sll_m_gyroaverage_2d_polar.F90:1600
sll_m_gyroaverage_2d_polar::splcoefper1d0
subroutine splcoefper1d0(luper, N)
Definition: sll_m_gyroaverage_2d_polar.F90:1883
sll_m_gyroaverage_2d_polar::sll_s_compute_gyroaverage_pre_compute_polar_spl_fft
subroutine, public sll_s_compute_gyroaverage_pre_compute_polar_spl_fft(gyro, f)
Definition: sll_m_gyroaverage_2d_polar.F90:914
sll_m_gyroaverage_2d_polar::interpolate_hermite_c1
subroutine interpolate_hermite_c1(f, i, x, fval, N)
Definition: sll_m_gyroaverage_2d_polar.F90:1552
sll_m_gyroaverage_2d_polar::splnat1d
subroutine splnat1d(f, xx, xmin, xmax, fval, N)
Definition: sll_m_gyroaverage_2d_polar.F90:1838
sll_m_gyroaverage_2d_polar::contribution_spl
subroutine contribution_spl(x, val)
Definition: sll_m_gyroaverage_2d_polar.F90:1637
sll_m_gyroaverage_2d_polar::localize_per
subroutine localize_per(i, x, xmin, xmax, N)
Definition: sll_m_gyroaverage_2d_polar.F90:1466
sll_m_gyroaverage_2d_polar::sll_s_penta
subroutine, public sll_s_penta(n, e, a, b, c, f, d, x)
Definition: sll_m_gyroaverage_2d_polar.F90:2000
sll_m_gyroaverage_2d_polar::splnatper2d
subroutine splnatper2d(f, xx, xmin, xmax, yy, ymin, ymax, fval, Nx, Ny)
Definition: sll_m_gyroaverage_2d_polar.F90:1755
sll_m_gyroaverage_utilities
Definition: sll_m_gyroaverage_utilities.F90:1
sll_m_gyroaverage_utilities::sll_s_compute_shape_circle
subroutine, public sll_s_compute_shape_circle(points, N_points)
Definition: sll_m_gyroaverage_utilities.F90:31
sll_m_gyroaverage_2d_polar::sll_t_plan_gyroaverage_polar
Definition: sll_m_gyroaverage_2d_polar.F90:59
    Report Typos and Errors
Generated on Mon Oct 23 2023 19:15:39 for Semi-Lagrangian Library by   doxygen 1.9.1