sll_m_qn_2d_polar.F90

Go to the documentation of this file.

!**************************************************************
!  Copyright INRIA
!  Authors :
!     CALVI project team
!
!  This code SeLaLib (for Semi-Lagrangian-Library)
!  is a parallel library for simulating the plasma turbulence
!  in a tokamak.
!
!  This software is governed by the CeCILL-B license
!  under French law and abiding by the rules of distribution
!  of free software.  You can  use, modify and redistribute
!  the software under the terms of the CeCILL-B license as
!  circulated by CEA, CNRS and INRIA at the following URL
!  "http://www.cecill.info".
!**************************************************************
 
module sll_m_qn_2d_polar
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#include "sll_memory.h"
#include "sll_working_precision.h"
 
   use sll_m_constants, only: sll_p_pi
 
   use sll_m_gnuplot, only: sll_o_gnuplot_2d
 
   use sll_m_fft, only: sll_t_fft, &
                        sll_s_fft_init_c2r_1d, &
                        sll_s_fft_exec_c2r_1d, &
                        sll_s_fft_init_r2c_1d, &
                        sll_s_fft_exec_r2c_1d, &
                        sll_s_fft_free
 
   use sll_m_gyroaverage_utilities, only: &
      sll_s_compute_shape_circle, &
      sll_s_zero_bessel_dir_dir
 
   implicit none
 
   public :: &
      sll_s_compute_error, &
      sll_s_compute_gamma0, &
      sll_f_compute_gamma0_quadrature, &
      sll_s_compute_splines_coefs_matrix_nat_1d, &
      sll_s_compute_splines_coefs_matrix_per_1d, &
      sll_s_contribution_spl, &
      sll_s_mu_quadr_for_phi_init, &
      sll_s_localize_polar, &
      sll_s_matrix_product_compf, &
      sll_s_qn_2d_polar_init, &
      sll_s_precompute_gyroaverage_index, &
      sll_s_precompute_inverse_qn_matrix_polar_splines, &
      sll_t_qn_2d_polar, &
      sll_s_qn_2d_polar_solve, &
      sll_s_splcoefnat1d0old, &
      sll_s_splcoefper1d0old, &
      sll_s_qn_2d_polar_test_solve
 
   private
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
   type sll_t_qn_2d_polar
 
      sll_real64          :: eta_min(2)     
      sll_real64          :: eta_max(2)     
      sll_int32           :: nc(2)          
      sll_int32           :: n_points       
 
      sll_real64, dimension(:, :), pointer   :: points
      sll_int32, dimension(:, :), pointer   :: pre_compute_n
      sll_real64, dimension(:, :, :), pointer   :: pre_compute_coeff
      sll_int32, dimension(:, :, :), pointer   :: pre_compute_index
      sll_real64, dimension(:, :), pointer   :: pre_compute_coeff_spl
      sll_int32 :: size_pre_compute
 
      sll_real64, dimension(:, :, :), pointer :: mat_gyro
      sll_real64, dimension(:, :, :), pointer :: mat_double_gyro
      sll_real64, dimension(:, :, :, :), pointer :: mat_gyro_circ
      sll_real64, dimension(:, :, :, :), pointer :: mat_double_gyro_circ
 
      sll_comp64, dimension(:, :, :), pointer :: mat_qn_inverse
 
      sll_real64, dimension(:), allocatable   :: lambda
      sll_real64, dimension(:), allocatable   :: t_i
      ! solve \lambda(r)\phi-\tilde{\phi} = second member
 
      sll_real64, dimension(:), pointer       :: mu_points_for_phi
      sll_real64, dimension(:), pointer       :: mu_weights_for_phi
      sll_int32                               :: n_mu_for_phi
 
   end type sll_t_qn_2d_polar
 
contains
 
   subroutine sll_s_qn_2d_polar_init(this, eta_min, eta_max, Nc, N_points, lambda, T_i)
 
      sll_real64, intent(in) :: eta_min(2)
      sll_real64, intent(in) :: eta_max(2)
      sll_int32, intent(in)  :: nc(2)
      sll_int32, intent(in)  :: n_points
      sll_real64, dimension(:), intent(in)    :: lambda
      sll_real64, dimension(:), intent(in)    :: t_i
      type(sll_t_qn_2d_polar) :: this
 
      sll_int32 :: err
 
      sll_allocate(this%points(3, n_points), err)
      sll_allocate(this%lambda(1:nc(1) + 1), err)
      sll_allocate(this%T_i(1:nc(1) + 1), err)
 
      call sll_s_compute_shape_circle(this%points, n_points)
 
      this%eta_min = eta_min
      this%eta_max = eta_max
      this%Nc = nc
      this%N_points = n_points
 
      this%lambda = lambda
      this%T_i = t_i
 
   end subroutine sll_s_qn_2d_polar_init
 
   subroutine sll_s_compute_splines_coefs_matrix_nat_1d(mat, dnat, lnat, N)
      sll_int32, intent(in)::n
      sll_real64, dimension(0:N + 2, 1:N + 1), intent(inout) :: mat
      sll_real64, dimension(0:N + 2), intent(in)::dnat, lnat
      sll_real64, dimension(:, :), allocatable :: matb
      sll_int32 :: err, i, j
 
      ! Conditions aux bords nulles
 
      sll_allocate(matb(0:n + 2, 0:n + 2), err)
 
      matb = 0._f64
 
      do j = 0, n + 2
         matb(j, j) = 6._f64
      end do
 
      do i = 1, n + 1
         do j = 0, n + 2
            matb(i, j) = matb(i, j) - lnat(i - 1)*matb(i - 1, j)
         end do
      end do
 
      do j = 0, n + 2
         matb(n + 2, j) = (matb(n + 2, j) - lnat(n + 1)*matb(n + 1, j) - lnat(n + 2)*matb(n, j))/dnat(n + 2)
      end do
 
      do i = n + 1, 2, -1
         do j = 0, n + 2
            matb(i, j) = (matb(i, j) - matb(i + 1, j))/dnat(i)
         end do
      end do
 
      do j = 0, n + 2
         matb(1, j) = (matb(1, j) - 2._f64*matb(2, j))/dnat(1)
      end do
 
      do j = 0, n + 2
         matb(0, j) = (matb(0, j) - 3._f64*matb(2, j))/dnat(0)
      end do
 
      mat(0:n + 2, 1:n + 1) = matb(0:n + 2, 1:n + 1)
 
   end subroutine sll_s_compute_splines_coefs_matrix_nat_1d
 
   subroutine sll_s_compute_splines_coefs_matrix_per_1d(mat, dper, lper, mper, N)
      sll_int32, intent(in)::n
      sll_real64, dimension(0:N - 1, 0:N - 1), intent(inout) :: mat
      sll_real64, dimension(0:N + 2), intent(in)::dper, lper, mper
      sll_int32 :: i, j
 
      mat = 0._f64
 
      do j = 0, n - 1
         mat(j, j) = 6._f64
      end do
 
      do i = 1, n - 1
         do j = 0, n - 1
            mat(i, j) = mat(i, j) - lper(i - 1)*mat(i - 1, j)
         end do
      end do
 
      do i = 0, n - 2
         do j = 0, n - 1
            mat(n - 1, j) = mat(n - 1, j) - mper(i)*mat(i, j)
         end do
      end do
 
      do j = 0, n - 1
         mat(n - 1, j) = mat(n - 1, j)*dper(n - 1)
      end do
 
      do j = 0, n - 1
         mat(n - 2, j) = dper(n - 2)*(mat(n - 2, j) - (1._f64 - mper(n - 3))*mat(n - 1, j))
      end do
 
      do i = n - 3, 1, -1
         do j = 0, n - 1
            mat(i, j) = dper(i)*(mat(i, j) - mat(i + 1, j) + mper(i - 1)*mat(n - 1, j))
         end do
      end do
 
      do j = 0, n - 1
         mat(0, j) = dper(0)*(mat(0, j) - mat(1, j) - mat(n - 1, j))
      end do
 
   end subroutine sll_s_compute_splines_coefs_matrix_per_1d
 
   subroutine kronecker_product(A, NxA, NyA, B, NxB, NyB, kronecker)
      sll_int32, intent(in)::nxa, nya, nxb, nyb
      sll_real64, dimension(0:NxA - 1, 0:NyA - 1), intent(in)::a
      sll_real64, dimension(0:NxB - 1, 0:NyB - 1), intent(in)::b
      sll_real64, dimension(0:NxA*NxB - 1, 0:NyA*NyB - 1), intent(inout) :: kronecker
      sll_int32 :: ia, ib, ja, jb
 
      do ia = 0, nxa - 1
         do ib = 0, nxb - 1
            do ja = 0, nya - 1
               do jb = 0, nyb - 1
                  kronecker(ib + nxb*ia, jb + nyb*ja) = a(ia, ja)*b(ib, jb)
               end do
            end do
         end do
      end do
 
   end subroutine kronecker_product
 
   subroutine matrix_product(A, NxA, NyA, B, NxB, NyB, prod)
      sll_int32, intent(in)::nxa, nya, nxb, nyb
      sll_real64, dimension(0:NxA - 1, 0:NyA - 1), intent(in)::a
      sll_real64, dimension(0:NxB - 1, 0:NyB - 1), intent(in)::b
      sll_real64, dimension(0:NxA - 1, 0:NyB - 1), intent(inout) :: prod
      sll_int32 :: i, j, k
      sll_real64 :: result
 
      if (nya /= nxb) then
         print *, '#incompatible sizes in matrix_product'
         stop
      else
         do i = 0, nxa - 1
            do j = 0, nyb - 1
               result = 0._f64
               do k = 0, nya - 1
                  result = result + a(i, k)*b(k, j)
               end do
               prod(i, j) = result
            end do
         end do
      end if
 
   end subroutine matrix_product
 
   subroutine sll_s_matrix_product_compf(A, NxA, NyA, B, NxB, NyB, prod)
      sll_int32, intent(in)::nxa, nya, nxb, nyb
      sll_comp64, dimension(:, :), intent(in)::a
      sll_comp64, dimension(:, :), intent(in)::b
      sll_comp64, dimension(:, :), intent(inout) :: prod
      sll_int32 :: i, j, k
      sll_comp64 :: result
 
      if (nya /= nxb) then
         print *, '#incompatible sizes in matrix_product'
         stop
      else
         do i = 1, nxa
            do j = 1, nyb
               result = (0.0_f64, 0.0_f64)
               do k = 1, nya
                  result = result + a(i, k)*b(k, j)
               end do
               prod(i, j) = result
            end do
         end do
      end if
 
   end subroutine sll_s_matrix_product_compf
 
   subroutine matrix_product_comp(A, NxA, NyA, B, NxB, NyB, prod)
      sll_int32, intent(in)::nxa, nya, nxb, nyb
      sll_comp64, dimension(0:NxA - 1, 0:NyA - 1), intent(in)::a
      sll_comp64, dimension(0:NxB - 1, 0:NyB - 1), intent(in)::b
      sll_comp64, dimension(0:NxA - 1, 0:NyB - 1), intent(inout) :: prod
      sll_int32 :: i, j, k
      sll_comp64 :: result
 
      if (nya /= nxb) then
         print *, '#incompatible sizes in matrix_product'
         stop
      else
         do i = 0, nxa - 1
            do j = 0, nyb - 1
               result = (0._f64, 0.0_f64)
               do k = 0, nya - 1
                  result = result + a(i, k)*b(k, j)
               end do
               prod(i, j) = result
            end do
         end do
      end if
 
   end subroutine matrix_product_comp
 
   subroutine matrix_product_circ(A, NxA, NyA, B, NxB, NyB, prod, Ntheta)
      sll_int32, intent(in)::ntheta, nxa, nya, nxb, nyb
      sll_real64, dimension(0:Ntheta - 1, 0:NxA - 1, 0:NyA - 1), intent(in)::a
      sll_real64, dimension(0:Ntheta - 1, 0:NxB - 1, 0:NyB - 1), intent(in)::b
      sll_real64, dimension(0:Ntheta - 1, 0:NxA - 1, 0:NyB - 1), intent(inout) :: prod
      sll_real64, dimension(:, :), allocatable :: mat_stock, mat_stock_sum
      sll_int32 :: i, j, error
      !sll_real64 :: result
 
      sll_allocate(mat_stock(0:nxa - 1, 0:nyb - 1), error)
      sll_allocate(mat_stock_sum(0:nxa - 1, 0:nyb - 1), error)
 
      if (nya /= nxb) then
         print *, '#incompatible sizes in matrix_product_circ'
         stop
      else
         do i = 0, ntheta - 1
            mat_stock_sum = 0._f64
            do j = 0, ntheta - 1
               call matrix_product(a(modulo(i + j, ntheta), :, :), nxa, nya, b(j, :, :), nxb, nyb, mat_stock)
               mat_stock_sum = mat_stock_sum + mat_stock
            end do
            prod(i, :, :) = mat_stock_sum
         end do
      end if
 
   end subroutine matrix_product_circ
 
   subroutine sll_s_precompute_gyroaverage_index(quasineutral, rho, N_rho)
      type(sll_t_qn_2d_polar)  :: quasineutral
      sll_int32 :: n_rho
      sll_real64, dimension(1:N_rho), intent(in) :: rho
      sll_int32, dimension(:, :), allocatable :: buf
      sll_int32 ::i, j, k, p, ell_1, ell_2, ii(2), s, nb, ind(2)
      sll_real64::val(-1:2, -1:2), eta_star(2), eta(2), delta_eta(2), x(2)
      sll_int32 ::error, max_nb
 
      delta_eta(1) = (quasineutral%eta_max(1) - quasineutral%eta_min(1))/real(quasineutral%Nc(1), f64)
      delta_eta(2) = (quasineutral%eta_max(2) - quasineutral%eta_min(2))/real(quasineutral%Nc(2), f64)
 
      sll_allocate(buf(0:quasineutral%Nc(1) + 2, 0:quasineutral%Nc(2) - 1), error)
      sll_allocate(quasineutral%pre_compute_N(1:n_rho, quasineutral%Nc(1) + 1), error)
 
      eta(2) = quasineutral%eta_min(2)
 
      max_nb = 0
 
      do p = 1, n_rho
 
         buf = 0
         nb = 0
         do i = 1, quasineutral%Nc(1) + 1
            eta(1) = quasineutral%eta_min(1) + real(i - 1, f64)*delta_eta(1)
            s = 0
            do k = 1, quasineutral%N_points
               x(1) = eta(1)*cos(eta(2)) + rho(p)*quasineutral%points(1, k)
               x(2) = eta(1)*sin(eta(2)) + rho(p)*quasineutral%points(2, k)
               call sll_s_localize_polar(x, quasineutral%eta_min, quasineutral%eta_max, ii, eta_star, quasineutral%Nc)
               do ell_2 = -1, 2
                  ind(2) = modulo(ii(2) + ell_2, quasineutral%Nc(2))
                  do ell_1 = -1, 2
                     ind(1) = ii(1) + 1 + ell_1
                     if (buf(ind(1), ind(2)) .ne. i) then
                        s = s + 1
                        buf(ind(1), ind(2)) = i
                     end if
                  end do
               end do
            end do
            quasineutral%pre_compute_N(p, i) = s
            nb = nb + s
         end do
         max_nb = max(max_nb, nb)
 
         quasineutral%size_pre_compute = max_nb
         ! print*,'#N_points pre_compute=',nb,quasineutral%Nc(1),nb/quasineutral%Nc(1)
 
      end do ! p=1,N_rho
 
      sll_allocate(quasineutral%pre_compute_index(1:n_rho, 1:2, 1:max_nb), error)
      sll_allocate(quasineutral%pre_compute_coeff_spl(1:n_rho, 1:max_nb), error)
 
      do p = 1, n_rho
 
         buf = 0
         nb = 0
         s = 0
         val = 0._f64
         eta(2) = quasineutral%eta_min(2)
         do i = 1, quasineutral%Nc(1) + 1
            eta(1) = quasineutral%eta_min(1) + real(i - 1, f64)*delta_eta(1)
            do k = 1, quasineutral%N_points
               x(1) = eta(1)*cos(eta(2)) + rho(p)*quasineutral%points(1, k)
               x(2) = eta(1)*sin(eta(2)) + rho(p)*quasineutral%points(2, k)
               call sll_s_localize_polar(x, quasineutral%eta_min, quasineutral%eta_max, ii, eta_star, quasineutral%Nc)
 
               call sll_s_contribution_spl(eta_star, val)
 
               val = val*quasineutral%points(3, k)
 
               do ell_2 = -1, 2
                  ind(2) = modulo(ii(2) + ell_2, quasineutral%Nc(2))
                  do ell_1 = -1, 2
                     ind(1) = ii(1) + 1 + ell_1
                     j = buf(ind(1), ind(2))
                     if (j <= nb) then
                        s = s + 1
                        buf(ind(1), ind(2)) = s
                        quasineutral%pre_compute_coeff_spl(p, s) = val(ell_1, ell_2)
                        quasineutral%pre_compute_index(p, 1, s) = ind(1)
                        quasineutral%pre_compute_index(p, 2, s) = ind(2)
                     else
                        quasineutral%pre_compute_coeff_spl(p, j) = quasineutral%pre_compute_coeff_spl(p, j) + val(ell_1, ell_2)
                     end if
                  end do
               end do
            end do
            nb = s
         end do
 
      end do ! p=1,N_rho
 
      !   print *,'#min/max=',minval(quasineutral%pre_compute_coeff_spl(:)),maxval(quasineutral%pre_compute_coeff_spl(:))
 
      sll_deallocate_array(buf, error)
 
   end subroutine sll_s_precompute_gyroaverage_index
 
   subroutine sll_s_precompute_inverse_qn_matrix_polar_splines(quasineutral, mu_points, mu_weights, N_mu)
      type(sll_t_qn_2d_polar) :: quasineutral
      sll_int32, intent(in) :: n_mu
      sll_real64, dimension(1:N_mu), intent(in) :: mu_points
      sll_real64, dimension(1:N_mu), intent(in) :: mu_weights
      sll_comp64, dimension(:, :), allocatable :: mat_stock1, mat_stock2
      sll_real64, dimension(:), allocatable, target::dnat, lnat, dper, lper, mper
      sll_real64, dimension(:), pointer::pointer_dnat, pointer_lnat
      sll_real64, dimension(:), pointer::pointer_dper, pointer_lper, pointer_mper
      sll_real64, dimension(:, :), allocatable, target::mat_nat, mat_per
      sll_real64, dimension(:, :), pointer::pointer_mat_nat, pointer_mat_per
      sll_real64, dimension(:, :, :), allocatable, target::mat_spl2d_circ
      sll_real64, dimension(:, :, :), pointer::pointer_mat_spl2d_circ
      sll_real64, dimension(:, :, :, :), allocatable, target::mat_contribution_circ
      sll_real64, dimension(:, :, :, :), pointer::pointer_mat_contribution_circ
      sll_comp64, dimension(:, :, :), allocatable::d_spl2d
      sll_comp64, dimension(:, :, :, :), allocatable::d_contr
      sll_int32, dimension(:), allocatable :: ipiv
      sll_comp64, dimension(:), allocatable :: work
      sll_real64 :: mode
      sll_comp64 :: exp_comp
      sll_int32 :: info
      sll_int32 :: nr, ntheta
      sll_int32 :: error
      sll_int32 :: i, j, k, m, p, s
      sll_int32 :: ii(2)
 
      ! Taille du maillage
      nr = quasineutral%Nc(1)
      ntheta = quasineutral%Nc(2)
 
      ! Allocate
      sll_allocate(mat_stock1(0:nr, 0:nr), error)
      sll_allocate(mat_stock2(0:nr, 0:nr), error)
      sll_allocate(quasineutral%mat_qn_inverse(0:ntheta - 1, 0:nr, 0:nr), error)
      sll_allocate(dnat(0:nr + 2), error)
      sll_allocate(lnat(0:nr + 2), error)
      sll_allocate(dper(0:ntheta - 1), error)
      sll_allocate(lper(0:ntheta - 1), error)
      sll_allocate(mper(0:ntheta - 1), error)
      sll_allocate(mat_nat(0:nr + 2, 0:nr), error)
      sll_allocate(mat_per(0:ntheta - 1, 0:ntheta - 1), error)
      sll_allocate(mat_spl2d_circ(0:ntheta - 1, 0:nr + 2, 0:nr), error)
      sll_allocate(d_spl2d(0:ntheta - 1, 0:nr + 2, 0:nr), error)
      sll_allocate(mat_contribution_circ(0:n_mu, 0:ntheta - 1, 0:nr, 0:nr + 2), error)
      sll_allocate(d_contr(0:n_mu, 0:ntheta - 1, 0:nr, 0:nr + 2), error)
      sll_allocate(ipiv(nr + 1), error)
      sll_allocate(work((nr + 1)**2), error)
 
      sll_allocate(quasineutral%mat_double_gyro_circ(1:n_mu, 0:ntheta - 1, 0:nr, 0:nr), error)
      sll_allocate(quasineutral%mat_gyro_circ(1:n_mu, 0:ntheta - 1, 0:nr, 0:nr), error)
 
      ! Initialise les coefficients de splines
      call sll_s_splcoefnat1d0old(dnat, lnat, nr)
      call sll_s_splcoefper1d0old(dper, lper, mper, ntheta)
 
      ! Pointeurs
      pointer_dnat => dnat
      pointer_lnat => lnat
      pointer_dper => dper
      pointer_lper => lper
      pointer_mper => mper
      pointer_mat_nat => mat_nat
      pointer_mat_per => mat_per
      pointer_mat_spl2d_circ => mat_spl2d_circ
      pointer_mat_contribution_circ => mat_contribution_circ
 
      ! Construction de la matrice des coefficients de splines
      call sll_s_compute_splines_coefs_matrix_nat_1d(pointer_mat_nat, pointer_dnat, pointer_lnat, nr)
      call sll_s_compute_splines_coefs_matrix_per_1d(pointer_mat_per, pointer_dper, pointer_lper, pointer_mper, ntheta)
      do j = 0, ntheta - 1
         pointer_mat_spl2d_circ(j, :, :) = pointer_mat_per(0, j)*pointer_mat_nat
      end do
 
      ! Construction de la matrice de contribution
      pointer_mat_contribution_circ = 0._f64
      do p = 1, n_mu
         s = 0
         do i = 1, nr + 1
            do k = 1, quasineutral%pre_compute_N(p, i)
               s = s + 1
               ii(1) = quasineutral%pre_compute_index(p, 1, s)
               ii(2) = modulo(quasineutral%pre_compute_index(p, 2, s), ntheta)
          pointer_mat_contribution_circ(p,ii(2),i-1,ii(1))=pointer_mat_contribution_circ(p,ii(2),i-1,ii(1))+quasineutral%pre_compute_coeff_spl(p,s) 
            end do
         end do
      end do
 
      ! Matrices D^{spl} et D^{contr}
      d_spl2d = (0.0_f64, 0.0_f64)
      d_contr = (0.0_f64, 0.0_f64)
      do m = 0, ntheta - 1
         do j = 0, ntheta - 1
            mode = real(-2._f64*sll_p_pi*real(j, f64)*real(m, f64)/real(ntheta, f64), f64)
            exp_comp = cmplx(cos(mode), sin(mode), kind=f64)
            d_spl2d(m, :, :) = d_spl2d(m, :, :) + pointer_mat_spl2d_circ(j, :, :)*exp_comp
            do p = 1, n_mu
               d_contr(p, m, :, :) = d_contr(p, m, :, :) + pointer_mat_contribution_circ(p, j, :, :)*exp_comp
            end do
         end do
      end do
 
      ! Construction de la matrice à inverser
      quasineutral%mat_qn_inverse = (0._f64, 0._f64)
      do m = 0, ntheta - 1
         do p = 1, n_mu
            mat_stock1 = (0._f64, 0._f64)
            mat_stock2 = (0._f64, 0._f64)
            call matrix_product_comp(d_contr(p, m, :, :), nr + 1, nr + 3, d_spl2d(m, :, :), nr + 3, nr + 1, mat_stock1)
            call matrix_product_comp(mat_stock1(:, :), nr + 1, nr + 1, mat_stock1(:, :), nr + 1, nr + 1, mat_stock2)
 
            !call matrix_product_comp(transpose(mat_stock1(:,:)),Nr+1,Nr+1,mat_stock1(:,:),Nr+1,Nr+1,mat_stock2)
 
            do i = 0, nr
               mat_stock2(i, i) = mat_stock2(i, i) - quasineutral%lambda(i + 1)*(1._f64, 0._f64)
            end do
            do i = 0, nr
               quasineutral%mat_qn_inverse(m, i, :) = &
                  quasineutral%mat_qn_inverse(m, i, :) &
                  - mu_weights(p)*mat_stock2(i, :)* &
                  cmplx(exp(-mu_points(p)/quasineutral%T_i(i + 1)), 0._f64, f64)
            end do
         end do
      end do
      ! Inversion des blocs
      do m = 0, ntheta - 1
         call zgetrf(nr + 1, nr + 1, quasineutral%mat_qn_inverse(m, :, :), nr + 1, ipiv, info)
         call zgetri(nr + 1, quasineutral%mat_qn_inverse(m, :, :), nr + 1, ipiv, work, (nr + 1)**2, info)
      end do
 
      ! Deallocate
      sll_deallocate_array(mat_stock1, error)
      sll_deallocate_array(mat_stock2, error)
      sll_deallocate_array(dnat, error)
      sll_deallocate_array(lnat, error)
      sll_deallocate_array(dper, error)
      sll_deallocate_array(lper, error)
      sll_deallocate_array(mper, error)
      sll_deallocate_array(mat_nat, error)
      sll_deallocate_array(mat_per, error)
      sll_deallocate_array(mat_spl2d_circ, error)
      sll_deallocate_array(d_spl2d, error)
      sll_deallocate_array(mat_contribution_circ, error)
      sll_deallocate_array(d_contr, error)
      sll_deallocate_array(ipiv, error)
      sll_deallocate_array(work, error)
 
   end subroutine sll_s_precompute_inverse_qn_matrix_polar_splines
 
   subroutine sll_s_qn_2d_polar_solve(quasineutral, phi)
      type(sll_t_qn_2d_polar) :: quasineutral
      sll_real64, dimension(1:quasineutral%Nc(1) + 1, 1:quasineutral%Nc(2)), intent(inout) :: phi
      sll_comp64, dimension(:, :), allocatable :: phi_comp, phi_old
      sll_real64, dimension(:), allocatable::buf_fft
      !type(sll_t_fft) :: fw, bw
      sll_comp64 :: result
      sll_int32 :: nr, ntheta
      sll_int32 :: error
      sll_int32 :: i, j, m
 
      ! Taille du maillage
      nr = quasineutral%Nc(1)
      ntheta = quasineutral%Nc(2)
 
      ! Allocate
      sll_allocate(phi_comp(1:nr + 1, 1:ntheta), error)
      sll_allocate(phi_old(1:nr + 1, 1:ntheta), error)
      sll_clear_allocate(buf_fft(1:4*ntheta + 15), error)
 
      ! FFT(PHI)
      phi_comp = phi*(1._f64, 0._f64)
      call zffti(ntheta, buf_fft)
      !call sll_s_fft_init_c2r_1d(fw, 2*ntheta-2,phi_comp(1,:),buf_fft)
      !call sll_s_fft_init_r2c_1d(bw, 2*ntheta-2,buf_fft,phi_comp(1,:))
      do i = 1, nr + 1
         call zfftf(ntheta, phi_comp(i, :), buf_fft)
         !call sll_s_fft_exec_c2r_1d(fw, phi_comp(i,:),buf_fft)
      end do
 
      ! Produit matrice/vecteur
      phi_old = phi_comp
      do m = 1, ntheta
         do i = 1, nr + 1
            result = (0._f64, 0.0_f64)
            do j = 1, nr + 1
               result = result + quasineutral%mat_qn_inverse(m - 1, i - 1, j - 1)*phi_old(j, m)
            end do
            phi_comp(i, m) = result
         end do
      end do
 
      ! FFT^-1
      do i = 1, nr + 1
         call zfftb(ntheta, phi_comp(i, :), buf_fft)
         !call sll_s_fft_exec_r2c_1d(bw, buf_fft, phi_comp(i,:))
      end do
      phi = real(phi_comp, f64)/real(ntheta, f64)
 
      ! Sorties
      !   print *,"phi_min : ",minval(phi)
      !   print *,"phi_max : ",maxval(phi)
 
      !call sll_s_fft_free(fw)
      !call sll_s_fft_free(bw)
      ! Deallocate
      sll_deallocate_array(phi_comp, error)
      sll_deallocate_array(phi_old, error)
      sll_deallocate_array(buf_fft, error)
 
   end subroutine sll_s_qn_2d_polar_solve
 
   subroutine sll_s_qn_2d_polar_test_solve(Nc, eta_min, eta_max, mu_points, mu_weights, N_mu, mode, lambda, T_i, phi_init, phi_qn)
      sll_int32, intent(in)  :: nc(2)
      sll_real64, intent(in) :: eta_min(2)
      sll_real64, intent(in) :: eta_max(2)
      sll_int32, intent(in) :: n_mu
      sll_int32, intent(in)  :: mode(2)
      sll_real64, dimension(1:N_mu), intent(in) :: mu_points
      sll_real64, dimension(1:N_mu), intent(in) :: mu_weights
      sll_real64, dimension(1:Nc(1) + 1, 1:Nc(2) + 1), intent(in) :: phi_init
      sll_real64, dimension(1:Nc(1) + 1, 1:Nc(2) + 1), intent(in) :: phi_qn
      sll_real64, dimension(1:Nc(1) + 1), intent(in) :: lambda, t_i
      sll_int32  :: n_min(2), n_max(2)
      sll_real64, dimension(:), allocatable :: gamma0
      sll_real64  :: tmp1
      sll_real64 :: eps, rho2d(2), mu2dmax(2), error(3)
      sll_int32 :: i, ierr, p
 
      sll_allocate(gamma0(1:nc(1) + 1), ierr)
 
      n_min = 1
      n_max(1) = nc(1)
      n_max(2) = nc(2)
 
      mu2dmax(1) = maxval(mu_points)
      mu2dmax(2) = maxval(mu_points)
 
      eps = 1.d-10
      gamma0 = 0._f64
 
      call compute_n_bounds_polar_circle(n_min(1), n_max(1), nc(1), 2._f64*sqrt(2._f64*mu2dmax), eta_min(1), eta_max(1))
 
      do p = 1, n_mu
         rho2d(1) = sqrt(2._f64*mu_points(p))
         rho2d(2) = sqrt(2._f64*mu_points(p))
         call solution_polar_circle(rho2d, mode, eta_min, eta_max, tmp1)
         do i = 1, nc(1) + 1
            gamma0(i) = gamma0(i) + mu_weights(p)*exp(-mu_points(p)/t_i(i))*(lambda(i) - tmp1**2)
         end do
      end do
      do i = 1, nc(1) + 1
         gamma0(i) = 1._f64/gamma0(i)
      end do
 
      ! print *,'#gamma0val=',gamma0
 
      print *, '#N_min(1:2) / N_max(1:2) = ', n_min, ' / ', n_max
      call compute_error_1d(phi_qn, phi_init, gamma0, error, n_min, n_max)
      print *, '#error subdomain=', error
      call compute_error_1d(phi_qn, phi_init, gamma0, error, (/1, 1/), nc)
      print *, '#error whole domain=', error
 
      call sll_o_gnuplot_2d( &
         eta_min(1), &
         eta_max(1), &
         nc(1) + 1, &
         eta_min(2), &
         eta_max(2), &
         nc(2) + 1, &
         phi_qn, &
         'fdiff', &
         1, &
         ierr)
 
   end subroutine sll_s_qn_2d_polar_test_solve
 
   function sll_f_compute_gamma0_quadrature( &
      Nc, &
      eta_min, &
      eta_max, &
      mu_points, &
      mu_weights, &
      N_mu, &
      mode) &
      result(gamma0)
      sll_int32, intent(in)  :: nc(2)
      sll_real64, intent(in) :: eta_min(2)
      sll_real64, intent(in) :: eta_max(2)
      sll_int32, intent(in) :: n_mu
      sll_int32, intent(in)  :: mode(2)
      sll_real64, dimension(1:N_mu), intent(in) :: mu_points
      sll_real64, dimension(1:N_mu), intent(in) :: mu_weights
      sll_real64 :: gamma0
      sll_real64 :: tmp1
      sll_real64 :: rho2d(2)
      sll_int32 :: p
 
      gamma0 = 0._f64
      do p = 1, n_mu
         rho2d(1) = sqrt(2._f64*mu_points(p))
         rho2d(2) = sqrt(2._f64*mu_points(p))
         call solution_polar_circle(rho2d, mode, eta_min, eta_max, tmp1)
         !gamma0 = gamma0 + mu_weights(p)*exp(-mu_points(p))*(1._f64-tmp1**2)
         gamma0 = gamma0 + mu_weights(p)*(1._f64 - tmp1**2)
         !gamma0 = gamma0 + mu_weights(p)*(tmp1**2)
      end do
 
      !gamma0 = 1._f64-gamma0
 
      return
      print *, nc
 
   end function sll_f_compute_gamma0_quadrature
 
   subroutine solve_circulant_system(Ntheta, Nr, mat_circ, sol)
      ! Solve mat_circ*X=sol where mat_circ(0:Ntheta-1,0:Nr,0:Nr)
      ! is a circulent matrix of size Ntheta*(Nr+1) Ntheta*(Nr+1)
      ! and sol(1:Nr+1,1:Ntheta)
      sll_int32                                :: ntheta
      sll_int32                                :: nr
      sll_int32                                :: i
      sll_int32                                :: j
      sll_int32                                :: m
      sll_int32                                :: error
      sll_comp64, dimension(:, :, :), allocatable :: dm
      sll_comp64, dimension(:, :), allocatable :: sol_comp
      sll_comp64, dimension(:, :), allocatable :: sol_old
      sll_real64, dimension(:), allocatable :: buf_fft
      sll_real64 :: mode
      sll_comp64 :: exp_comp, result
      sll_int32, dimension(:), allocatable :: ipiv
      sll_comp64, dimension(:), allocatable :: work
      sll_int32 :: info
      sll_real64, dimension(0:Ntheta - 1, 0:Nr, 0:Nr), intent(in) :: mat_circ
      sll_real64, dimension(1:Nr + 1, 1:Ntheta), intent(inout) :: sol
      !type(sll_t_fft) :: fw, bw
 
      sll_allocate(dm(0:ntheta - 1, 0:nr, 0:nr), error)
      sll_allocate(sol_comp(1:nr + 1, 1:ntheta), error)
      sll_allocate(sol_old(1:nr + 1, 1:ntheta), error)
      sll_allocate(buf_fft(1:4*ntheta + 15), error)
      sll_allocate(ipiv(nr + 1), error)
      sll_allocate(work((nr + 1)**2), error)
 
      sol_comp = sol*(1._f64, 0._f64)
 
      ! FFT(PHI)
      call zffti(ntheta, buf_fft)
      !call sll_s_fft_init_c2r_1d(fw,2*ntheta-2,sol_comp(1,:),buf_fft)
      !call sll_s_fft_init_r2c_1d(bw,2*ntheta-2,buf_fft,sol_comp(1,:))
      do i = 1, nr + 1
         call zfftf(ntheta, sol_comp(i, :), buf_fft)
         !call sll_s_fft_exec_c2r_1d(fw,sol_comp(i,:),buf_fft)
      end do
 
      ! Matrices Dm
      dm = (0._f64, 0.0_f64)
      do m = 0, ntheta - 1
         do j = 0, ntheta - 1
            mode = real(-2._f64*sll_p_pi*real(j, f64)*real(m, f64)/real(ntheta, f64), f64)
            exp_comp = cmplx(cos(mode), sin(mode), kind=f64)
            dm(m, :, :) = dm(m, :, :) + mat_circ(j, :, :)*exp_comp
         end do
 
         call zgetrf(nr + 1, nr + 1, dm(m, :, :), nr + 1, ipiv, info)
         call zgetri(nr + 1, dm(m, :, :), nr + 1, ipiv, work, (nr + 1)**2, info)
 
      end do
 
      ! Produit Dm
      sol_old = sol_comp
      do m = 1, ntheta
         do i = 1, nr + 1
            result = (0._f64, 0.0_f64)
            do j = 1, nr + 1
               result = result + dm(m - 1, i - 1, j - 1)*sol_old(j, m)
            end do
            sol_comp(i, m) = result
         end do
      end do
 
      ! FFT^-1
      do i = 1, nr + 1
         call zfftb(ntheta, sol_comp(i, :), buf_fft)
         !call sll_s_fft_exec_r2c_1d(bw,buf_fft,sol_comp(i,:))
      end do
 
      sol = real(sol_comp/cmplx(ntheta, 0._f64, f64), f64)
 
      !call sll_s_fft_free(fw)
      !call sll_s_fft_free(bw)
 
   end subroutine solve_circulant_system
 
   subroutine test_solve_circulant_system(mat, Nr, Ntheta)
      sll_int32, intent(in) :: nr, ntheta
      sll_real64, dimension(0:Ntheta - 1, 0:Nr, 0:Nr), intent(in) :: mat
      sll_real64, dimension(:, :, :), allocatable :: inv, test_id
      sll_real64, dimension(:, :), allocatable :: phi
      sll_int32 :: i, j, m, error
 
      sll_allocate(inv(0:ntheta - 1, 0:nr, 0:nr), error)
      sll_allocate(test_id(0:ntheta - 1, 0:nr, 0:nr), error)
      sll_allocate(phi(1:nr + 1, 1:ntheta), error)
 
      do m = 1, nr + 1
         phi = 0._f64
         phi(m, 1) = 1._f64
 
         call solve_circulant_system(ntheta, nr, mat, phi)
 
         do i = 0, nr
            do j = 0, ntheta - 1
               inv(modulo(ntheta - j, ntheta), i, m - 1) = phi(i + 1, j + 1)
            end do
         end do
      end do
 
      call matrix_product_circ(mat, nr + 1, nr + 1, inv, nr + 1, nr + 1, test_id, ntheta)
 
      do i = 0, nr
         test_id(0, i, i) = test_id(0, i, i) - 1._f64
      end do
 
      print *, "Inversion error = ", maxval(test_id)
 
   end subroutine test_solve_circulant_system
 
   !----------------------------------------
 
   subroutine hermite_coef_nat_per(f, buf3d, N, d)
      sll_int32, intent(in)::n(2), d(2)
      sll_real64, dimension(N(1) + 1, N(2)), intent(in)::f
      sll_real64, dimension(9, N(1) + 1, N(2) + 1), intent(out)::buf3d
      sll_real64 ::w_left_1(-d(1)/2:(d(1) + 1)/2), w_right_1((-d(1) + 1)/2:d(1)/2 + 1)
      sll_real64 ::w_left_2(-d(2)/2:(d(2) + 1)/2), w_right_2((-d(2) + 1)/2:d(2)/2 + 1)
      sll_real64 ::tmp
      sll_int32  ::i, j, ii, r_left(2), r_right(2), s_left(2), s_right(2), ind
      r_left = -d/2
      s_left = (d + 1)/2
      r_right = (-d + 1)/2
      s_right = d/2 + 1
 
      call compute_w_hermite(w_left_1, r_left(1), s_left(1))
      call compute_w_hermite(w_left_2, r_left(2), s_left(2))
      if ((2*(d(1)/2) - d(1)) == 0) then
         w_right_1(r_right(1):s_right(1)) = w_left_1(r_left(1):s_left(1))
      else
         w_right_1(r_right(1):s_right(1)) = -w_left_1(s_left(1):r_left(1):-1)
      end if
 
      if ((2*(d(2)/2) - d(2)) == 0) then
         w_right_2(r_right(2):s_right(2)) = w_left_2(r_left(2):s_left(2))
      else
         w_right_2(r_right(2):s_right(2)) = -w_left_2(s_left(2):r_left(2):-1)
      end if
 
      !print *,'w(',r_left(1),':',s_left(1),')=',w_left_1(r_left(1):s_left(1))
      !print *,'w(',r_right(1),':',s_right(1),')=',w_right_1(r_right(1):s_right(1))
 
      do j = 1, n(2)
         do i = 1, n(1) + 1
            buf3d(1, i, j) = f(i, j) !f(0,0)
            tmp = 0._f64
            do ii = r_left(1), s_left(1)
               ind = i + ii; if (ind > n(1) + 1) ind = 2*(n(1) + 1) - ind; if (ind < 1) ind = 2 - ind
               tmp = tmp + w_left_1(ii)*f(ind, j)
            end do
            buf3d(2, i, j) = tmp !fx(0,0)
            tmp = 0._f64
            do ii = r_right(1), s_right(1)
               ind = i + ii; if (ind > n(1) + 1) ind = 2*(n(1) + 1) - ind; if (ind < 1) ind = 2 - ind
               tmp = tmp + w_right_1(ii)*f(ind, j)
            end do
            buf3d(3, i, j) = tmp !fx(1,0)
         end do
      end do
      do i = 1, n(1) + 1
         do j = 1, n(2)
            tmp = 0._f64
            do ii = r_left(2), s_left(2)
               ind = modulo(j + ii - 1 + n(2), n(2)) + 1
               tmp = tmp + w_left_2(ii)*f(i, ind)
            end do
            buf3d(4, i, j) = tmp !fy(0,0)
            tmp = 0._f64
            do ii = r_right(2), s_right(2)
               ind = modulo(j + ii - 1 + n(2), n(2)) + 1
               tmp = tmp + w_right_2(ii)*f(i, ind)
            end do
            buf3d(5, i, j) = tmp !fy(0,1)
         end do
      end do
 
      do j = 1, n(2)
         do i = 1, n(1) + 1
            tmp = 0._f64
            do ii = r_left(1), s_left(1)
               ind = i + ii; if (ind > n(1) + 1) ind = 2*(n(1) + 1) - ind; if (ind < 1) ind = 2 - ind
               tmp = tmp + w_left_1(ii)*buf3d(4, ind, j)
            end do
            buf3d(6, i, j) = tmp !fxy(0,0)
            tmp = 0._f64
            do ii = r_right(1), s_right(1)
               ind = i + ii; if (ind > n(1) + 1) ind = 2*(n(1) + 1) - ind; if (ind < 1) ind = 2 - ind
               tmp = tmp + w_right_1(ii)*buf3d(4, ind, j)
            end do
            buf3d(7, i, j) = tmp !fxy(1,0)
            tmp = 0._f64
            do ii = r_left(1), s_left(1)
               ind = i + ii; if (ind > n(1) + 1) ind = 2*(n(1) + 1) - ind; if (ind < 1) ind = 2 - ind
               tmp = tmp + w_left_1(ii)*buf3d(5, ind, j)
            end do
            buf3d(8, i, j) = tmp  !fxy(0,1)
            tmp = 0._f64
            do ii = r_right(1), s_right(1)
               ind = i + ii; if (ind > n(1) + 1) ind = 2*(n(1) + 1) - ind; if (ind < 1) ind = 2 - ind
               tmp = tmp + w_right_1(ii)*buf3d(5, ind, j)
            end do
            buf3d(9, i, j) = tmp !fxy(1,1)
         end do
      end do
 
      buf3d(:, :, n(2) + 1) = buf3d(:, :, 1)
 
   end subroutine hermite_coef_nat_per
 
   subroutine hermite_c1_coef_nat_per(f, buf3d, N, d)
      integer, intent(in)::N(2), d(2)
      sll_real64, dimension(N(1) + 1, N(2)), intent(in)::f
      sll_real64, dimension(4, N(1) + 1, N(2) + 1), intent(out)::buf3d
      sll_real64, dimension(:), allocatable ::w_left_1, w_left_2
      sll_real64 ::tmp
      sll_int32::i, j, ii, r_left(2), s_left(2), ind, dd(2)
      sll_int32::error
      dd(1) = 2*((d(1) + 1)/2)
      dd(2) = 2*((d(2) + 1)/2)
 
      r_left = -dd/2
      s_left = (dd + 1)/2
 
      sll_allocate(w_left_1(-dd(1)/2:dd(1)/2), error)
      sll_allocate(w_left_2(-dd(2)/2:dd(2)/2), error)
 
      call compute_w_hermite(w_left_1, r_left(1), s_left(1))
      call compute_w_hermite(w_left_2, r_left(2), s_left(2))
 
      !print *,'w(',r_left(1),':',s_left(1),')=',w_left_1(r_left(1):s_left(1))
      !print *,'w(',r_right(1),':',s_right(1),')=',w_right_1(r_right(1):s_right(1))
 
      do j = 1, n(2)
         do i = 1, n(1) + 1
            buf3d(1, i, j) = f(i, j) !f(0,0)
            tmp = 0._f64
            do ii = r_left(1), s_left(1)
               ind = i + ii; if (ind > n(1) + 1) ind = 2*(n(1) + 1) - ind; if (ind < 1) ind = 2 - ind
               tmp = tmp + w_left_1(ii)*f(ind, j)
            end do
            buf3d(2, i, j) = tmp !fx(0,0)
         end do
      end do
      do i = 1, n(1) + 1
         do j = 1, n(2)
            tmp = 0._f64
            do ii = r_left(2), s_left(2)
               ind = modulo(j + ii - 1 + n(2), n(2)) + 1
               tmp = tmp + w_left_2(ii)*f(i, ind)
            end do
            buf3d(3, i, j) = tmp !fy(0,0)
         end do
      end do
 
      do j = 1, n(2)
         do i = 1, n(1) + 1
            tmp = 0._f64
            do ii = r_left(1), s_left(1)
               ind = i + ii; if (ind > n(1) + 1) ind = 2*(n(1) + 1) - ind; if (ind < 1) ind = 2 - ind
               tmp = tmp + w_left_1(ii)*buf3d(3, ind, j)
            end do
            buf3d(4, i, j) = tmp !fxy(0,0)
         end do
      end do
 
      buf3d(:, :, n(2) + 1) = buf3d(:, :, 1)
 
      sll_deallocate_array(w_left_1, error)
      sll_deallocate_array(w_left_2, error)
 
   end subroutine hermite_c1_coef_nat_per
 
   subroutine compute_w_hermite(w, r, s)
      sll_int32, intent(in)::r, s
      sll_real64, dimension(r:s), intent(out)::w
      sll_int32 ::i, j
      sll_real64::tmp
 
      !maple code for generation of w
      !for k from r to -1 do
      !  C[k]:=product((k-j),j=r..k-1)*product((k-j),j=k+1..s):
      !  C[k]:=1/C[k]*product((-j),j=r..k-1)*product((-j),j=k+1..-1)*product((-j),j=1..s):
      !od:
      !for k from 1 to s do
      !  C[k]:=product((k-j),j=r..k-1)*product((k-j),j=k+1..s):
      !  C[k]:=1/C[k]*product((-j),j=r..-1)*product((-j),j=1..k-1)*product((-j),j=k+1..s):
      !od:
      !C[0]:=-add(C[k],k=r..-1)-add(C[k],k=1..s):
 
      do i = r, -1
         tmp = 1._f64
         do j = r, i - 1
            tmp = tmp*real(i - j, f64)
         end do
         do j = i + 1, s
            tmp = tmp*real(i - j, f64)
         end do
         tmp = 1._f64/tmp
         do j = r, i - 1
            tmp = tmp*real(-j, f64)
         end do
         do j = i + 1, -1
            tmp = tmp*real(-j, f64)
         end do
         do j = 1, s
            tmp = tmp*real(-j, f64)
         end do
         w(i) = tmp
      end do
 
      do i = 1, s
         tmp = 1._f64
         do j = r, i - 1
            tmp = tmp*real(i - j, f64)
         end do
         do j = i + 1, s
            tmp = tmp*real(i - j, f64)
         end do
         tmp = 1._f64/tmp
         do j = r, -1
            tmp = tmp*real(-j, f64)
         end do
         do j = 1, i - 1
            tmp = tmp*real(-j, f64)
         end do
         do j = i + 1, s
            tmp = tmp*real(-j, f64)
         end do
         w(i) = tmp
      end do
      tmp = 0._f64
      do i = r, -1
         tmp = tmp + w(i)
      end do
      do i = 1, s
         tmp = tmp + w(i)
      end do
      w(0) = -tmp
 
      !print *,'w',w
      !do ii=r,s
      !  print *,ii,w(r+s-ii)
      !enddo
      !
 
   end subroutine compute_w_hermite
 
   subroutine sll_s_localize_polar(x, eta_min, eta_max, ii, eta, N)
      sll_real64, intent(in)::x(2), eta_min(2), eta_max(2)
      sll_int32, intent(out)::ii(2)
      sll_int32, intent(in)::n(2)
      sll_real64, intent(out)::eta(2)
 
      eta(1) = sqrt(x(1)**2 + x(2)**2)
      call localize_nat(ii(1), eta(1), eta_min(1), eta_max(1), n(1))
      eta(2) = atan2(x(2), x(1))
      call localize_per(ii(2), eta(2), eta_min(2), eta_max(2), n(2))
   end subroutine sll_s_localize_polar
 
   subroutine localize_per(i, x, xmin, xmax, N)
      sll_int32, intent(out)::i
      sll_real64, intent(inout)::x
      sll_real64, intent(in)::xmin, xmax
      sll_int32, intent(in)::n
      x = (x - xmin)/(xmax - xmin)
      x = x - real(floor(x), f64)
      x = x*real(n, f64)
      i = floor(x)
      x = x - real(i, f64)
      if (i == n) then
         i = 0
         x = 0._f64
      end if
   end subroutine localize_per
 
   subroutine localize_nat(i, x, xmin, xmax, N)
      sll_int32, intent(out)::i
      sll_real64, intent(inout)::x
      sll_real64, intent(in)::xmin, xmax
      sll_int32, intent(in)::n
      x = (x - xmin)/(xmax - xmin)
      x = x*real(n, f64)
      if (x >= real(n, f64)) then
         x = real(n, f64)
      end if
      if (x <= 0._f64) then
         x = 0._f64
      end if
      i = floor(x)
      x = x - real(i, f64)
      if (i == n) then
         i = n - 1
         x = 1._f64
      end if
   end subroutine localize_nat
 
   subroutine interpolate_hermite(f, i, x, fval, N)
      sll_int32, intent(in)::i(2), n(2)
      !real(f64),intent(in)::xx(2),xmin(2),xmax(2)
      sll_real64, intent(in)::x(2)
      sll_real64, intent(out)::fval
      sll_real64, dimension(0:8, 0:N(1), 0:N(2))::f
      !integer::i(2),i1(2),s
      sll_int32::i1(2), s
      sll_real64::w(2, 0:3), tmp(0:3)
      sll_real64::g(0:3, 0:3)
 
      !fval =f(0,i(1),i(2))!real(i(1),f64)
 
      !return
 
      do s = 1, 2
         w(s, 0) = (2._f64*x(s) + 1)*(1._f64 - x(s))*(1._f64 - x(s)); 
         w(s, 1) = x(s)*x(s)*(3._f64 - 2._f64*x(s))
         w(s, 2) = x(s)*(1._f64 - x(s))*(1._f64 - x(s))
         w(s, 3) = x(s)*x(s)*(x(s) - 1._f64)
         i1(s) = i(s) + 1
      end do
 
      g(0, 0) = f(0, i(1), i(2))          !f(0,0)
      g(1, 0) = f(0, i1(1), i(2))         !f(1,0)
      g(2, 0) = f(1, i(1), i(2))          !fx(0,0)
      g(3, 0) = f(2, i(1), i(2))          !fx(1,0)
      g(0, 1) = f(0, i(1), i1(2))         !f(0,1)
      g(1, 1) = f(0, i1(1), i1(2))        !f(1,1)
      g(2, 1) = f(1, i(1), i1(2))         !fx(0,1)
      g(3, 1) = f(2, i(1), i1(2))         !fx(1,1)
      g(0, 2) = f(3, i(1), i(2))          !fy(0,0)
      g(1, 2) = f(3, i1(1), i(2))         !fy(1,0)
      g(2, 2) = f(5, i(1), i(2))          !fxy(0,0)
      g(3, 2) = f(6, i(1), i(2))          !fxy(1,0)
      g(0, 3) = f(4, i(1), i(2))          !fy(0,1)
      g(1, 3) = f(4, i1(1), i(2))         !fy(1,1)
      g(2, 3) = f(7, i(1), i(2))          !fxy(0,1)
      g(3, 3) = f(8, i(1), i(2))          !fxy(1,1)
 
      do s = 0, 3
         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)
      end do
 
      fval = w(2, 0)*tmp(0) + w(2, 1)*tmp(1) + w(2, 2)*tmp(2) + w(2, 3)*tmp(3)
 
      !print *,fval,' t',f
   end subroutine interpolate_hermite
 
   subroutine interpolate_hermite_c1(f, i, x, fval, N)
      sll_int32, intent(in)::i(2), n(2)
      !real(f64),intent(in)::xx(2),xmin(2),xmax(2)
      sll_real64, intent(in)::x(2)
      sll_real64, intent(out)::fval
      sll_real64, dimension(0:3, 0:N(1), 0:N(2))::f
      !integer::i(2),i1(2),s
      sll_int32::i1(2), s
      sll_real64::w(2, 0:3), tmp(0:3)
      sll_real64::g(0:3, 0:3)
 
      !fval =f(0,i(1),i(2))!real(i(1),f64)
 
      !return
 
      do s = 1, 2
         w(s, 0) = (2._f64*x(s) + 1)*(1._f64 - x(s))*(1._f64 - x(s)); 
         w(s, 1) = x(s)*x(s)*(3._f64 - 2._f64*x(s))
         w(s, 2) = x(s)*(1._f64 - x(s))*(1._f64 - x(s))
         w(s, 3) = x(s)*x(s)*(x(s) - 1._f64)
         i1(s) = i(s) + 1
      end do
 
      g(0, 0) = f(0, i(1), i(2))          !f(0,0)
      g(1, 0) = f(0, i1(1), i(2))         !f(1,0)
      g(2, 0) = f(1, i(1), i(2))          !fx(0,0)
      g(3, 0) = f(1, i1(1), i(2))          !fx(1,0)
      g(0, 1) = f(0, i(1), i1(2))         !f(0,1)
      g(1, 1) = f(0, i1(1), i1(2))        !f(1,1)
      g(2, 1) = f(1, i(1), i1(2))         !fx(0,1)
      g(3, 1) = f(1, i1(1), i1(2))         !fx(1,1)
      g(0, 2) = f(2, i(1), i(2))          !fy(0,0)
      g(1, 2) = f(2, i1(1), i(2))         !fy(1,0)
      g(2, 2) = f(3, i(1), i(2))          !fxy(0,0)
      g(3, 2) = f(3, i1(1), i(2))          !fxy(1,0)
      g(0, 3) = f(2, i(1), i1(2))          !fy(0,1)
      g(1, 3) = f(2, i1(1), i1(2))         !fy(1,1)
      g(2, 3) = f(3, i(1), i1(2))          !fxy(0,1)
      g(3, 3) = f(3, i1(1), i1(2))          !fxy(1,1)
 
      do s = 0, 3
         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)
      end do
 
      fval = w(2, 0)*tmp(0) + w(2, 1)*tmp(1) + w(2, 2)*tmp(2) + w(2, 3)*tmp(3)
      !print *,fval,' t',f
   end subroutine interpolate_hermite_c1
 
   subroutine contribution_hermite_c1(x, val)
      sll_real64, intent(in)::x(0:1)
      sll_real64, intent(out)::val(4, 0:1, 0:1)
      sll_int32::s
      sll_real64::w(0:3, 0:1)!,tmp(0:3)
      do s = 0, 1
         w(0, s) = (2._f64*x(s) + 1)*(1._f64 - x(s))*(1._f64 - x(s)); 
         w(1, s) = x(s)*x(s)*(3._f64 - 2._f64*x(s))
         w(2, s) = x(s)*(1._f64 - x(s))*(1._f64 - x(s))
         w(3, s) = x(s)*x(s)*(x(s) - 1._f64)
      end do
 
      val(1, 0, 0) = w(0, 0)*w(0, 1)
      val(2, 0, 0) = w(2, 0)*w(0, 1)
      val(3, 0, 0) = w(0, 0)*w(2, 1)
      val(4, 0, 0) = w(2, 0)*w(2, 1)
 
      val(1, 1, 0) = w(1, 0)*w(0, 1)
      val(2, 1, 0) = w(3, 0)*w(0, 1)
      val(3, 1, 0) = w(1, 0)*w(2, 1)
      val(4, 1, 0) = w(3, 0)*w(2, 1)
 
      val(1, 0, 1) = w(0, 0)*w(1, 1)
      val(2, 0, 1) = w(2, 0)*w(1, 1)
      val(3, 0, 1) = w(0, 0)*w(3, 1)
      val(4, 0, 1) = w(2, 0)*w(3, 1)
 
      val(1, 1, 1) = w(1, 0)*w(1, 1)
      val(2, 1, 1) = w(3, 0)*w(1, 1)
      val(3, 1, 1) = w(1, 0)*w(3, 1)
      val(4, 1, 1) = w(3, 0)*w(3, 1)
 
      !print *,val
      !stop
 
   end subroutine contribution_hermite_c1
 
   subroutine sll_s_contribution_spl(x, val)
      sll_real64, intent(in)::x(0:1)
      sll_real64, intent(out)::val(-1:2, -1:2)
      sll_int32::s
      sll_real64::w(-1:2, 0:1)
      do s = 0, 1
         w(-1, s) = (1._f64/6._f64)*(1._f64 - x(s))*(1._f64 - x(s))*(1._f64 - x(s)); 
         w(0, s) = 1._f64/6._f64 + 0.5_f64*(1._f64 - x(s))*(-(1._f64 - x(s))* &
                                                            (1._f64 - x(s)) + (1._f64 - x(s)) + 1._f64); 
         w(1, s) = 1._f64/6._f64 + 0.5_f64*x(s)*(-x(s)*x(s) + x(s) + 1._f64); 
         w(2, s) = (1._f64/6._f64)*x(s)*x(s)*x(s); 
      end do
 
      val(-1, -1) = w(-1, 0)*w(-1, 1)
      val(-1, 0) = w(-1, 0)*w(0, 1)
      val(-1, 1) = w(-1, 0)*w(1, 1)
      val(-1, 2) = w(-1, 0)*w(2, 1)
 
      val(0, -1) = w(0, 0)*w(-1, 1)
      val(0, 0) = w(0, 0)*w(0, 1)
      val(0, 1) = w(0, 0)*w(1, 1)
      val(0, 2) = w(0, 0)*w(2, 1)
 
      val(1, -1) = w(1, 0)*w(-1, 1)
      val(1, 0) = w(1, 0)*w(0, 1)
      val(1, 1) = w(1, 0)*w(1, 1)
      val(1, 2) = w(1, 0)*w(2, 1)
 
      val(2, -1) = w(2, 0)*w(-1, 1)
      val(2, 0) = w(2, 0)*w(0, 1)
      val(2, 1) = w(2, 0)*w(1, 1)
      val(2, 2) = w(2, 0)*w(2, 1)
 
   end subroutine sll_s_contribution_spl
 
   subroutine splcoefnat1d0(lunat, N)
      sll_int32, intent(in)::n
      sll_real64, dimension(0:1, 0:N + 2), intent(out)::lunat
      sll_int32::i
      sll_real64::a(0:2)
      a(0) = -3._f64; a(1) = 0._f64; a(2) = 3._f64; 
      lunat(0, 0) = a(0); 
      lunat(1, 0) = 1._f64/lunat(0, 0); 
      lunat(0, 1) = 4._f64 - a(1)*lunat(1, 0); 
      lunat(1, 1) = 1._f64/lunat(0, 1); 
      lunat(0, 2) = 4._f64 - lunat(1, 1)*(1._f64 - a(2)/a(0)); 
      do i = 2, n
         lunat(1, i) = 1._f64/lunat(0, i); 
         lunat(0, i + 1) = 4._f64 - lunat(1, i); 
      end do
      lunat(1, n + 2) = a(0)/lunat(0, n); 
      lunat(1, n + 1) = (a(1) - lunat(1, n + 2))/lunat(0, n + 1); 
      lunat(0, n + 2) = a(2) - lunat(1, n + 1); 
   end subroutine splcoefnat1d0
 
   subroutine splcoefnat1d(p, lunat, N)
      sll_int32, intent(in)::n
      sll_real64, dimension(0:1, 0:N + 2), intent(in)::lunat
      sll_real64, dimension(0:N + 2), intent(inout)::p
      sll_int32::i
      sll_real64::a(0:2)
      a(0) = -3._f64; a(1) = 0._f64; a(2) = 3._f64; 
      !p(0)=0._f64;
      !p(N+2)=0._f64;
      do i = 0, n + 2; 
         p(i) = 6._f64*p(i); 
      end do; 
      do i = 1, n + 1; p(i) = p(i) - lunat(1, i - 1)*p(i - 1); end do
      p(n + 2) = p(n + 2) - (lunat(1, n + 1)*p(n + 1) + lunat(1, n + 2)*p(n)); 
      p(n + 2) = p(n + 2)/lunat(0, n + 2); 
      do i = n + 1, 2, -1; p(i) = (p(i) - p(i + 1))/lunat(0, i); end do
      p(1) = (p(1) - (1._f64 - a(2)/a(0))*p(2))/lunat(0, 1); 
      p(0) = (p(0) - a(1)*p(1) - a(2)*p(2))/lunat(0, 0); 
      !p(i-1)+4*p(i)+p(i+1)=6ftab(i-1,j), i=1..N+1
      !a(0)*p(0)+a(1)*p(1)+a(2)*p(2)=f'(rmin)=0);
      !a(0)*p(Na)+a(1)*p(Na+1)+a(2)*p(Na+2)=f'(rmax)=0);
   end subroutine splcoefnat1d
 
   subroutine sll_s_splcoefnat1d0old(dnat, lnat, N)
      sll_int32, intent(in)::n
      sll_real64, dimension(0:N + 2), intent(inout)::dnat, lnat
      sll_int32::i
      sll_real64::a(0:2)
      a(0) = -3._f64; a(1) = 0._f64; a(2) = 3._f64; 
      dnat(0) = a(0); 
      lnat(0) = 1._f64/dnat(0); 
      dnat(1) = 4._f64 - a(1)*lnat(0); 
      lnat(1) = 1._f64/dnat(1); 
      dnat(2) = 4._f64 - lnat(1)*(1._f64 - a(2)/a(0)); 
      do i = 2, n
         lnat(i) = 1._f64/dnat(i); 
         dnat(i + 1) = 4._f64 - lnat(i); 
      end do
      lnat(n + 2) = a(0)/dnat(n); 
      lnat(n + 1) = (a(1) - lnat(n + 2))/dnat(n + 1); 
      dnat(n + 2) = a(2) - lnat(n + 1); 
   end subroutine sll_s_splcoefnat1d0old
 
   subroutine splcoefnatper2d(f, buf, dnatx, lnatx, dpery, lpery, mpery, Nx, Ny)
      sll_int32, intent(in)::nx, ny
      sll_real64, dimension(:, :), pointer::f
      sll_real64, dimension(:), pointer::buf, dpery, lpery, mpery, dnatx, lnatx
      sll_int32::i, j
      !natural spline coefficients in x
      do j = 0, ny - 1
         do i = 0, nx + 2; buf(i) = f(i, j); end do
         call splcoefnat1dold(buf, dnatx, lnatx, nx)
         do i = 0, nx + 2; f(i, j) = buf(i); end do
      end do
      !periodic spline coefficients in y
      do i = 0, nx + 2
         do j = 0, ny - 1; buf(j) = f(i, j); end do
         call splcoefper1dold(buf, dpery, lpery, mpery, ny)
         do j = 0, ny - 1; f(i, j) = buf(j); end do
      end do
 
   end subroutine splcoefnatper2d
 
   subroutine splnatper2d(f, xx, xmin, xmax, yy, ymin, ymax, fval, Nx, Ny)
      sll_int32, intent(in)::nx, ny
      sll_real64, intent(in)::xx, xmin, xmax, yy, ymin, ymax
      sll_real64, intent(out)::fval
      sll_real64, dimension(:, :), pointer::f
      sll_int32::i, j
      sll_int32::ix(0:3), iy(0:3)
      sll_real64::x, y
      sll_real64::wx(0:3), wy(0:3), tmp(0:3)
 
      x = (xx - xmin)/(xmax - xmin)
      !x=x-real(floor(x),f64)
      if (x < 0._f64) x = 0._f64; 
      if (x >= 1._f64) x = 1._f64 - 1.e-12_f64; 
      x = x*real(nx, f64)
      i = floor(x)
      x = x - real(i, f64)
 
      y = (yy - ymin)/(ymax - ymin)
      y = y - real(floor(y), f64)
      y = y*real(ny, f64)
      j = floor(y)
      y = y - real(j, f64)
 
      wx(0) = (1._f64/6._f64)*(1._f64 - x)*(1._f64 - x)*(1._f64 - x); 
      wx(1) = 1._f64/6._f64 + 0.5_f64*(1._f64 - x)*(-(1._f64 - x)* &
                                                    (1._f64 - x) + (1._f64 - x) + 1._f64); 
      wx(2) = 1._f64/6._f64 + 0.5_f64*x*(-x*x + x + 1._f64); 
      wx(3) = (1._f64/6._f64)*x*x*x; 
      wy(0) = (1._f64/6._f64)*(1._f64 - y)*(1._f64 - y)*(1._f64 - y); 
      wy(1) = 1._f64/6._f64 + 0.5_f64*(1._f64 - y)*(-(1._f64 - y)* &
                                                    (1._f64 - y) + (1._f64 - y) + 1._f64); 
      wy(2) = 1._f64/6._f64 + 0.5_f64*y*(-y*y + y + 1._f64); 
      wy(3) = (1._f64/6._f64)*y*y*y; 
      iy(0) = mod(j + ny - 1, ny)
      iy(1) = j
      iy(2) = mod(j + 1, ny)
      iy(3) = mod(j + 2, ny)
 
      ix(0) = i
      ix(1) = i + 1
      ix(2) = i + 2
      ix(3) = i + 3
 
      tmp(0) = wx(0)*f(ix(0), iy(0)) + wx(1)*f(ix(1), iy(0)) &
               + wx(2)*f(ix(2), iy(0)) + wx(3)*f(ix(3), iy(0))
      tmp(1) = wx(0)*f(ix(0), iy(1)) + wx(1)*f(ix(1), iy(1)) &
               + wx(2)*f(ix(2), iy(1)) + wx(3)*f(ix(3), iy(1))
      tmp(2) = wx(0)*f(ix(0), iy(2)) + wx(1)*f(ix(1), iy(2)) &
               + wx(2)*f(ix(2), iy(2)) + wx(3)*f(ix(3), iy(2))
      tmp(3) = wx(0)*f(ix(0), iy(3)) + wx(1)*f(ix(1), iy(3)) &
               + wx(2)*f(ix(2), iy(3)) + wx(3)*f(ix(3), iy(3))
 
      fval = wy(0)*tmp(0) + wy(1)*tmp(1) + wy(2)*tmp(2) + wy(3)*tmp(3)
 
   end subroutine splnatper2d
 
   subroutine splper1d(f, xx, xmin, xmax, fval, N)
      sll_int32, intent(in)::n
      sll_real64, intent(in)::xx, xmin, xmax
      sll_real64, intent(out)::fval
      sll_real64, dimension(0:N - 1), intent(inout)::f
      sll_int32::i
      real(f64)::x
      real(f64)::w(0:3)
      x = (xx - xmin)/(xmax - xmin)
      x = x - real(floor(x), f64)
      x = x*real(n, f64)
      i = floor(x)
      x = x - real(i, f64)
      w(0) = (1._f64/6._f64)*(1._f64 - x)*(1._f64 - x)*(1._f64 - x); 
      w(1) = 1._f64/6._f64 + 0.5_f64*(1._f64 - x)*(-(1._f64 - x)* &
                                                   (1._f64 - x) + (1._f64 - x) + 1._f64); 
      w(2) = 1._f64/6._f64 + 0.5_f64*x*(-x*x + x + 1._f64); 
      w(3) = (1._f64/6._f64)*x*x*x; 
      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))
   end subroutine splper1d
 
   subroutine splnat1d(f, xx, xmin, xmax, fval, N)
      sll_int32, intent(in)::n
      sll_real64, intent(in)::xx, xmin, xmax
      sll_real64, intent(out)::fval
      sll_real64, dimension(0:N + 2), intent(inout)::f
      sll_int32::i
      sll_real64::x
      sll_real64::w(0:3)
      x = (xx - xmin)/(xmax - xmin)
      !x=x-real(floor(x),f64)
      if (x < 0._f64) x = 0._f64; 
      if (x > 1._f64) x = 1._f64 - 1.e-12_f64; 
      x = x*real(n, f64)
      i = floor(x)
      x = x - real(i, f64)
 
      w(0) = (1._f64/6._f64)*(1._f64 - x)*(1._f64 - x)*(1._f64 - x); 
      w(1) = 1._f64/6._f64 + 0.5_f64*(1._f64 - x)*(-(1._f64 - x)* &
                                                   (1._f64 - x) + (1._f64 - x) + 1._f64); 
      w(2) = 1._f64/6._f64 + 0.5_f64*x*(-x*x + x + 1._f64); 
      w(3) = (1._f64/6._f64)*x*x*x; 
      fval = w(0)*f(i) + w(1)*f(i + 1) + w(2)*f(i + 2) + w(3)*f(i + 3)
   end subroutine splnat1d
 
   subroutine sll_s_splcoefper1d0old(dper, lper, mper, N)
      sll_int32, intent(in)::n
      sll_real64, dimension(0:N - 1), intent(out)::dper, lper, mper
      sll_int32::i
 
      dper(0) = 4._f64
      mper(0) = 0.25_f64
      do i = 0, n - 2
         lper(i) = 1._f64/dper(i)
         dper(i + 1) = 4._f64 - lper(i)
         mper(i + 1) = -mper(i)/dper(i + 1)
      end do
      dper(n - 1) = dper(n - 1) - (lper(n - 2) + 2._f64*mper(n - 2))
      do i = 0, n - 1
         dper(i) = 1._f64/dper(i)
      end do
   end subroutine sll_s_splcoefper1d0old
 
   subroutine splcoefper1d0(luper, N)
      sll_int32, intent(in)::n
      sll_real64, dimension(0:3*N - 1), intent(out)::luper
      sll_int32::i
 
      luper(0 + 3*0) = 4._f64
      luper(2 + 3*0) = 0.25_f64
      do i = 0, n - 2
         luper(1 + 3*i) = 1._f64/luper(0 + 3*i)
         luper(0 + 3*(i + 1)) = 4._f64 - luper(1 + 3*i)
         luper(2 + 3*(i + 1)) = -luper(2 + 3*i)/luper(0 + 3*(i + 1))
      end do
      luper(0 + 3*(n - 1)) = luper(0 + 3*(n - 1)) - (luper(1 + 3*(n - 2)) + 2._f64*luper(2 + 3*(n - 2)))
      do i = 0, n - 1
         luper(0 + 3*i) = 1._f64/luper(0 + 3*i)
      end do
   end subroutine splcoefper1d0
 
   subroutine splcoefper1d(f, luper, N)
      sll_int32, intent(in)::n
      sll_real64, dimension(0:3*N - 1), intent(in)::luper
      sll_real64, dimension(0:N - 1), intent(inout)::f
      sll_int32::i
      do i = 0, n - 1; f(i) = 6._f64*f(i); end do; 
      do i = 1, n - 1
         f(i) = f(i) - f(i - 1)*luper(1 + 3*(i - 1))
      end do
      do i = 0, n - 2
         f(n - 1) = f(n - 1) - luper(2 + 3*i)*f(i)
      end do
      f(n - 1) = f(n - 1)*luper(0 + 3*(n - 1)); f(n - 2) = luper(0 + 3*(n - 2))*(f(n - 2) - (1._f64 - luper(2 + 3*(n - 3)))*f(n - 1))
      do i = n - 3, 1, -1
         f(i) = luper(0 + 3*i)*(f(i) - f(i + 1) + luper(2 + 3*(i - 1))*f(n - 1))
      end do
      f(0) = luper(0 + 3*0)*(f(0) - f(1) - f(n - 1)); 
   end subroutine splcoefper1d
 
   subroutine splcoefnat1dold(p, dnat, lnat, N)
      sll_int32, intent(in)::n
      sll_real64, dimension(0:N + 2), intent(in)::dnat, lnat
      sll_real64, dimension(0:N + 2), intent(inout)::p
      sll_int32::i
      sll_real64::a(0:2)
      a(0) = -3._f64; a(1) = 0._f64; a(2) = 3._f64; 
      !p(0)=0._f64;
      !p(N+2)=0._f64;
      do i = 0, n + 2; 
         p(i) = 6._f64*p(i); 
      end do; 
      do i = 1, n + 1; p(i) = p(i) - lnat(i - 1)*p(i - 1); end do
      p(n + 2) = p(n + 2) - (lnat(n + 1)*p(n + 1) + lnat(n + 2)*p(n)); 
      p(n + 2) = p(n + 2)/dnat(n + 2); 
      do i = n + 1, 2, -1; p(i) = (p(i) - p(i + 1))/dnat(i); end do
      p(1) = (p(1) - (1._f64 - a(2)/a(0))*p(2))/dnat(1); 
      p(0) = (p(0) - a(1)*p(1) - a(2)*p(2))/dnat(0); 
      !p(i-1)+4*p(i)+p(i+1)=6ftab(i-1,j), i=1..N+1
      !a(0)*p(0)+a(1)*p(1)+a(2)*p(2)=f'(rmin)=0);
      !a(0)*p(Na)+a(1)*p(Na+1)+a(2)*p(Na+2)=f'(rmax)=0);
   end subroutine splcoefnat1dold
 
   subroutine splcoefper1dold(f, dper, lper, mper, N)
      sll_int32, intent(in)::n
      sll_real64, dimension(0:N - 1), intent(in)::dper, lper, mper
      sll_real64, dimension(0:N - 1), intent(inout)::f
      sll_int32::i
      do i = 0, n - 1; f(i) = 6._f64*f(i); end do; 
      do i = 1, n - 1
         f(i) = f(i) - f(i - 1)*lper(i - 1)
      end do
      do i = 0, n - 2
         f(n - 1) = f(n - 1) - mper(i)*f(i)
      end do
      f(n - 1) = f(n - 1)*dper(n - 1); f(n - 2) = dper(n - 2)*(f(n - 2) - (1._f64 - mper(n - 3))*f(n - 1))
      do i = n - 3, 1, -1
         f(i) = dper(i)*(f(i) - f(i + 1) + mper(i - 1)*f(n - 1))
      end do
      f(0) = dper(0)*(f(0) - f(1) - f(n - 1)); 
   end subroutine splcoefper1dold
 
   subroutine solve_tridiag(a, b, c, v, x, n)
! Using Thomas' Algorithm
      implicit none
!        a - sub-diagonal (means it is the diagonal below the main diagonal)
!        b - the main diagonal
!        c - sup-diagonal (means it is the diagonal above the main diagonal)
!        v - right part
!        x - the answer
!        n - number of equations
 
      sll_int32, intent(in) :: n
      sll_real64, dimension(n), intent(in) :: a, b, c, v
      sll_real64, dimension(n), intent(out) :: x
      sll_real64, dimension(n) :: bp, vp
      sll_real64 :: m
      sll_int32 i
 
      ! Make copies of the b and v variables so that they are unaltered by this sub
      bp(1) = b(1)
      vp(1) = v(1)
 
      !The first pass (setting coefficients):
      firstpass: do i = 2, n
         m = a(i)/bp(i - 1)
         bp(i) = b(i) - m*c(i - 1)
         vp(i) = v(i) - m*vp(i - 1)
      end do firstpass
 
      x(n) = vp(n)/bp(n)
      !The second pass (back-substition)
      backsub: do i = n - 1, 1, -1
         x(i) = (vp(i) - c(i)*x(i + 1))/bp(i)
      end do backsub
 
   end subroutine solve_tridiag
 
   subroutine compute_n_bounds_polar_circle(N_min, N_max, N, rho, eta_min, eta_max)
      ! version modifiée : N -> N+1
      sll_int32, intent(out)::n_min, n_max
      sll_int32, intent(in)::n
      sll_real64, intent(in)::rho(2), eta_min, eta_max
      sll_real64::delta_eta
      sll_int32::n_rho
      delta_eta = (eta_max - eta_min)/real(n, f64)
      n_rho = floor(max(rho(1), rho(2))/delta_eta) + 1
      n_min = 1 + n_rho
      n_max = n + 1 - n_rho
      if ((n_min > n + 1) .or. (n_max < 1) .or. (n_min > n_max)) then
         print *, '#Warning: rho is too big'
         print *, '#Bad computation of N_min and N_max'
         n_min = 1
         n_max = n + 1
      end if
   end subroutine compute_n_bounds_polar_circle
 
   subroutine compute_factor_bounds(f, f_init, eps, bounds, N_min, N_max)
      sll_real64, dimension(:, :), intent(in)::f, f_init
      sll_int32, intent(in)::n_min(2), n_max(2)
      sll_real64, intent(in)::eps
      sll_real64, intent(out)::bounds(2)
      sll_int32::i, j
      sll_real64::tmp
 
      bounds(1) = 1.e10_f64
      bounds(2) = -bounds(1)
      do j = n_min(2), n_max(2) + 1
         do i = n_min(1), n_max(1) + 1
            if (abs(f_init(i, j)) > eps) then
               tmp = f(i, j)/f_init(i, j)
               if (tmp > bounds(2)) bounds(2) = tmp; 
               if (tmp < bounds(1)) bounds(1) = tmp; 
            end if
         end do
      end do
   end subroutine compute_factor_bounds
 
   subroutine sll_s_compute_error(f, f_init, J_factor, err, N_min, N_max)
      sll_real64, dimension(:, :), intent(in)::f, f_init
      sll_int32, intent(in)::n_min(2), n_max(2)
      sll_real64, intent(in)::j_factor
      sll_real64, intent(out)::err(3)
      sll_int32::i, j
      sll_real64::tmp, delta
 
      err = 0._f64
 
      do j = n_min(2), n_max(2) + 1
         do i = n_min(1), n_max(1) + 1
            tmp = f(i, j) - j_factor*f_init(i, j)
            err(1) = err(1) + abs(tmp)
            err(2) = err(2) + abs(tmp)**2
            err(3) = max(err(3), abs(tmp))
         end do
      end do
      delta = 1._f64/(real((n_max(2) - n_min(2)), f64)*real((n_max(1) - n_min(1)), f64))
      err(1) = err(1)*delta
      err(2) = sqrt(err(2)*delta)
 
   end subroutine sll_s_compute_error
 
   subroutine compute_error_1d(f, f_init, J_factor, err, N_min, N_max)
      sll_real64, dimension(:, :), intent(in)::f, f_init
      sll_int32, intent(in)::n_min(2), n_max(2)
      sll_real64, dimension(:), intent(in) :: j_factor
      sll_real64, intent(out)::err(3)
      sll_int32::i, j
      sll_real64::tmp, delta
 
      err = 0._f64
 
      do j = n_min(2), n_max(2) + 1
         do i = n_min(1), n_max(1) + 1
            tmp = f(i, j) - j_factor(i)*f_init(i, j)
            err(1) = err(1) + abs(tmp)
            err(2) = err(2) + abs(tmp)**2
            err(3) = max(err(3), abs(tmp))
         end do
      end do
      delta = 1._f64/(real((n_max(2) - n_min(2)), f64)*real((n_max(1) - n_min(1)), f64))
      err(1) = err(1)*delta
      err(2) = sqrt(err(2)*delta)
 
   end subroutine compute_error_1d
 
   subroutine solution_polar_circle(rho, mode, eta_min, eta_max, val)
      sll_real64, intent(in)::rho(2), eta_min(2), eta_max(2)
      sll_int32, intent(in)::mode(2)
      sll_real64, intent(out)::val
      sll_real64::tmp
      !integer::i
      !logical::is_file
      val = 0._f64
      if (abs(rho(1) - rho(2)) > 1.e-12) then
         print *, '#for the moment rho(1)=rho(2) is needed'
         return
      end if
      call sll_s_zero_bessel_dir_dir(mode, eta_min(1), eta_max(1), tmp)
      val = 0._f64 !temporary because DBESJ not recognized on helios and curie
      !val = DBESJN(0,tmp*rho(1)/eta_max(1))
      !print *,i,j,mode_max,alpha,tmp
   end subroutine solution_polar_circle
 
   subroutine sll_s_compute_gamma0(mode, eta_min, eta_max, val)
      sll_real64, intent(in)::eta_min(2), eta_max(2)
      sll_int32, intent(in)::mode(2)
      sll_real64, intent(out)::val
      sll_real64::tmp!,mu,mu_max
      !integer::N_approx
      !logical::is_file
!    N_approx = 2**15
!    mu_max = 50._f64
!    sum1=0._f64
!    sum2=0._f64
      call sll_s_zero_bessel_dir_dir(mode, eta_min(1), eta_max(1), tmp)
      !print *,"tmp=",tmp
      !print *,"eta_max(1)=",eta_max(1)
      if (abs(tmp - 3.9651944700162098_f64) > 1.e-12) then
         print *, '#error : gamma0 not computed for this mode : ', mode
         return
      end if
      if (abs(eta_max(1) - 18._f64) > 1.e-12) then
         print *, '#error : gamma0 not computed for this eta_max : ', eta_max(1)
         return
      end if
      val = 0.95319247622058476357_f64
!    delta_x=real(real(mu_max,f64)/real(N_approx,f64),f64)
!    do i=1,N_approx/2-1
!      mu = real(2._f64*real(i,f64)*delta_x,f64)
!      sum1 = sum1 + DBESJN(0,tmp*sqrt(2._f64*mu)/eta_max(1))**2*exp(-mu)
!    enddo
!    do i=1,N_approx/2
!      mu = real((2._f64*real(i,f64)-1._f64)*delta_x,f64)
!      sum2 = sum2 + DBESJN(0,tmp*sqrt(2._f64*mu)/eta_max(1))**2*exp(-mu)
!    enddo
!    val = 2._f64*sum1 + 4._f64*sum2 + DBESJN(0,0._f64)**2 + DBESJN(0,tmp*sqrt(2._f64*mu_max)/eta_max(1))**2*exp(-mu_max)
!    val = val*real(delta_x/3._f64,f64)
   end subroutine sll_s_compute_gamma0
 
   subroutine sll_s_mu_quadr_for_phi_init( &
      quasineutral, &
      mu_quadr_for_phi_case, &
      N_mu_for_phi, &
      mu_max_for_phi, &
      mu_points_user_defined, &
      mu_weights_user_defined, &
      N_mu_user_defined)
 
      type(sll_t_qn_2d_polar)                        :: quasineutral
      character(len=256), intent(in)                 :: mu_quadr_for_phi_case
      sll_int32, intent(in)                          :: n_mu_for_phi
      sll_real64, intent(in)                         :: mu_max_for_phi
      sll_real64, dimension(:), intent(in), optional :: mu_points_user_defined
      sll_real64, dimension(:), intent(in), optional :: mu_weights_user_defined
      sll_int32, intent(in)                          :: n_mu_user_defined
      sll_real64 :: h
      sll_int32 :: ierr, i
 
      if (mu_quadr_for_phi_case == "SLL_USER_DEFINED") then
         quasineutral%N_mu_for_phi = n_mu_user_defined
         sll_allocate(quasineutral%mu_points_for_phi(1:n_mu_user_defined), ierr)
         sll_allocate(quasineutral%mu_weights_for_phi(1:n_mu_user_defined), ierr)
         if (.not. ((present(mu_points_user_defined)) &
                    .and. (present(mu_weights_user_defined)))) then
            print *, '#provide mu_points_user_defined, mu_weights_user_defined'
            print *, '#in sll_s_mu_quadr_for_phi_init'
            stop
         end if
!    else
!      if(&
!        (present(mu_points_user_defined)) &
!        .or.(present(mu_weights_user_defined)) &
!        )then
!        print *,'# do not provide mu_points_user_defined, mu_weights_user_defined'
!        print *,'#in sll_s_mu_quadr_for_phi_init'
!        stop
!      endif
      end if
 
      select case (mu_quadr_for_phi_case)
      case ("SLL_USER_DEFINED")
         if (size(mu_points_user_defined) /= n_mu_user_defined) then
            print *, '#incompatible sizes for mu_points_user_defined :', size(mu_points_user_defined)
            print *, '#and N_mu_user_defined :', n_mu_user_defined
            print *, '#in sll_s_mu_quadr_for_phi_init'
            stop
         elseif (size(mu_weights_user_defined) /= n_mu_user_defined) then
            print *, '#incompatible sizes for mu_weights_user_defined :', size(mu_weights_user_defined)
            print *, '#and N_mu_user_defined :', n_mu_user_defined
            print *, '#in sll_s_mu_quadr_for_phi_init'
            stop
         else
            quasineutral%mu_points_for_phi(1:n_mu_user_defined) = mu_points_user_defined
            quasineutral%mu_weights_for_phi(1:n_mu_user_defined) = mu_weights_user_defined
         end if
      case ("SLL_LIN_LEE")
         quasineutral%N_mu_for_phi = n_mu_for_phi
         sll_allocate(quasineutral%mu_points_for_phi(1:n_mu_for_phi), ierr)
         sll_allocate(quasineutral%mu_weights_for_phi(1:n_mu_for_phi), ierr)
         select case (n_mu_for_phi)
         case (1)
            quasineutral%mu_points_for_phi(1) = 1._f64
            quasineutral%mu_weights_for_phi(1) = 1._f64
         case (2)
            quasineutral%mu_points_for_phi(1) = 0.4167845_f64
            quasineutral%mu_points_for_phi(2) = 2.495154605_f64
            quasineutral%mu_weights_for_phi(1) = 0.7194_f64
            quasineutral%mu_weights_for_phi(2) = 0.2806_f64
         case (3)
            quasineutral%mu_points_for_phi(1) = 0.148131245_f64
            quasineutral%mu_points_for_phi(2) = 1._f64
            quasineutral%mu_points_for_phi(3) = 3.15959522_f64
            quasineutral%mu_weights_for_phi(1) = 0.3583_f64
            quasineutral%mu_weights_for_phi(2) = 0.5004_f64
            quasineutral%mu_weights_for_phi(3) = 0.1413_f64
         case default
            print *, '# bad N_mu_user_defined in SLL_LIN_LEE (must be 1, 2 or 3) : ', n_mu_for_phi
            print *, '#in sll_s_mu_quadr_for_phi_init'
            stop
         end select
      case ("SLL_RECTANGLES")
         quasineutral%N_mu_for_phi = n_mu_for_phi
         sll_allocate(quasineutral%mu_points_for_phi(1:n_mu_for_phi), ierr)
         sll_allocate(quasineutral%mu_weights_for_phi(1:n_mu_for_phi), ierr)
         do i = 1, n_mu_for_phi
            quasineutral%mu_points_for_phi(i) = mu_max_for_phi*real(i - 1, f64)/real(n_mu_for_phi, f64)
            quasineutral%mu_weights_for_phi(i) = mu_max_for_phi/real(n_mu_for_phi, f64)
         end do
      case ("SLL_SIMPSON")
         quasineutral%N_mu_for_phi = n_mu_for_phi
         sll_allocate(quasineutral%mu_points_for_phi(1:n_mu_for_phi), ierr)
         sll_allocate(quasineutral%mu_weights_for_phi(1:n_mu_for_phi), ierr)
         do i = 1, n_mu_for_phi
            quasineutral%mu_points_for_phi(i) = mu_max_for_phi*real(i - 1, f64)/real(n_mu_for_phi - 1, f64)
         end do
         h = mu_max_for_phi/(3._f64*real(n_mu_for_phi, f64))
         quasineutral%mu_weights_for_phi(1) = h
         do i = 1, (n_mu_for_phi - 1)/2 - 1
            quasineutral%mu_weights_for_phi(2*i) = 4._f64*h
            quasineutral%mu_weights_for_phi(2*i + 1) = 2._f64*h
         end do
         quasineutral%mu_weights_for_phi(n_mu_for_phi - 1) = 4._f64*h
         quasineutral%mu_weights_for_phi(n_mu_for_phi) = h
      case default
         print *, '#mu_quadr_for_phi_case not defined'
         print *, '#in sll_s_mu_quadr_for_phi_init'
         stop
      end select
   end subroutine sll_s_mu_quadr_for_phi_init
 
end module sll_m_qn_2d_polar