mesh_accumulation.f90

Go to the documentation of this file.

subroutine mesh_accumulation(nrtmp, nv, s_acc, sig, weights, &
 equidistant, sg, dsg, ddsg, in_psi)
!-----------------------------------------------------------------------
! subroutine to construct non-equidistant radial mesh in helena
! for in_psi == .true. accumulation in psi, else in s
!-----------------------------------------------------------------------

  use itm_types
  use mod_helena_io, only : iu6
  use mod_profiles
  use helena_spline

  implicit none

  integer(itm_i4), intent(in) :: nrtmp, nv
  real(r8), dimension(nv), intent(in) :: s_acc, sig, weights
  real(r8), intent(in) :: equidistant
  logical, intent(in) :: in_psi  

  real(r8), dimension(nrtmp), intent(out) :: sg, dsg, ddsg
  
  integer(itm_i4), parameter :: nmax = 1001

  type(spline_coefficients) :: f_spline
  real(r8), allocatable :: s1(:), fsum(:)
  real(r8), dimension(3) :: abltg
  real(r8) :: alfa, beta
  real(r8) :: fi
  real(r8) :: dx
  real(r8) :: weight, gauss_integral
  integer(itm_i4) :: typ
  integer(itm_i4) :: i, n_int
  
  real(r8) :: fgauss

  if (any(s_acc(1 : nv) < 0._r8) .or. any(s_acc(1 : nv) > 1._r8)) then
    write(iu6, 2)
    stop
  else if (any(sig(1 : nv) <= 0._r8)) then
    write(iu6, 3)
    stop
  else if (any(weights(1 : nv) < 0._r8) .or. any(weights(1 : nv) > 1)) then
    write(iu6, 4)
    stop
  end if

  weight = 0._r8
  do i = 1, nv
    weight = weight + weights(i)
  end do
  if (weight <= 0._r8) then
    write(iu6, 5)
    stop
  end if

! --- first integrate the complete Gauss function for normalization
  n_int = max(nint(1._r8 / minval(sig(1 : nv))), nmax)

  allocate(s1(n_int))
  call allocate_spline_coefficients(f_spline, n_int)
  allocate(fsum(n_int))

  dx = 1._r8 / dble(n_int - 1)
  gauss_integral = 0._r8
  fsum(1) = 0._r8
  do i = 1, n_int - 1
    gauss_integral = gauss_integral + fgauss(dble(i - 0.5_r8) * dx, &
     nv, equidistant, s_acc(1 : nv), weights(1 : nv), sig(1 : nv)) * dx
    fsum(i + 1) = gauss_integral
  end do

  do i = 1, n_int
    s1(i) = dble(i - 1) / dble(n_int - 1)
    fsum(i) = fsum(i) / fsum(n_int)
  end do

  if (in_psi) fsum = fsum**2

  alfa = 0._r8
  beta = 0._r8
  typ = 3
  call spline(n_int, fsum, s1, alfa, beta, typ, f_spline)
  
  do i = 1, nrtmp
    fi = dble(i - 1) / dble(nrtmp - 1)
    sg(i) = spwert(n_int, fi, f_spline, fsum, abltg, 2)
    dsg(i) = abltg(1)
    ddsg(i) = abltg(2)
  end do
  
  deallocate(s1)
  call deallocate_spline_coefficients(f_spline)
  deallocate(fsum)

 2  format(//, 1x, 'ERROR:', /, 1x, &
     'accumulation points s_acc must be >= 0 and <= 1', /, 1x, &
     'stop in subroutine mesh_accumulation')
 3  format(//, 1x, 'ERROR:', /, 1x, &
     'Gaussian widths sig must be > 0', /, 1x, &
     'stop in subroutine mesh_accumulation')
 4  format(//, 1x, 'ERROR:', /, 1x, &
     'weighting factors weights must be >= 0 and <= 1', /, 1x, &
     'stop in subroutine mesh_accumulation')
 5  format(//, 1x, 'ERROR:', /, 1x, &
     'at least one weighting factor must be > 0', /, 1x, &
     'stop in subroutine mesh_accumulation')
  return

end subroutine mesh_accumulation