flxint.f90

Go to the documentation of this file.

subroutine flxint(xaxis, a, nmg, fscale)
!----------------------------------------------------------------------
! subroutine to evaluate new fdf profile
!----------------------------------------------------------------------
  
  use itm_types
  use mod_parameter
  use mod_corners
  use mod_dat
  use mod_gauss
  use mod_mesh
  use mod_profiles
  use mod_interpolation
  use mod_spline
  use current_calculations, only : j0, cscale
  use phys_profiles, only : j_phi
  use phys_functions
  use math_functions, only : flux_surface_average
  use matrix_calculations
  use helena_spline

  implicit none
  
  real(r8), intent(in) :: a, xaxis
  real(r8), intent(inout) :: fscale
  integer(itm_i4), intent(in) :: nmg

  type(spline_coefficients) :: fdf_spline

  real(r8), dimension(nr) :: psi_new, fdf_new
  real(r8), dimension(3) :: abltg
  real(r8) :: factas, cur0
  real(r8) :: r, s, ws
  real(r8) :: ps, psr, pss, psrs, psrr, psss
  integer(itm_i4) :: i, nelm
  integer(itm_i4) :: n1, n2, n3, n4
  integer(itm_i4) :: error
  
  real(r8) :: dp_dpsi

  real(r8) :: pnorm
 
  if (ias == 1) then
    factas = 1._r8
  else
    factas = 2._r8
  end if

!TODO: adapt B such that cur0 = - j_tor(1) / (eps * A * j0) ???
  if (nmg > 0) then
    j0 = bvac / (mu0 * eps * rvac * alfa)
    pnorm = eps * bvac**2 / (mu0 * alfa**2)
    cur0 = -cscale / (eps * a * j0)
!    write(*, *) 'cur0 = ', cur0
!    write(*, *) 'old b = ', b
!    b = -pscale * sign(1._r8, cpsurfin) / (eps * pnorm * a)
!    if (hbt) then
!      b = (cur0 * (1._r8 + eps * xaxis) - 1._r8) / (xaxis * (1._r8 &
!       + eps * xaxis / 2._r8))
!    else
!      b = (cur0 * (1._r8 + eps * xaxis) - 1._r8) / (1._r8 + eps &
!       * xaxis)**2
!    end if
  end if
!  write(*, *) 'new b = ', b

  if (hbt) then
    cur0 = (1._r8 + b * xaxis * (1._r8 + eps * xaxis / 2._r8)) &
     / (1._r8 + eps * xaxis)
  else
    cur0 = (1._r8 + b * (1._r8 + eps * xaxis)**2) / (1._r8 + eps &
     * xaxis)
  end if

  r = -1._r8
  s = xgauss(1)
  ws = wgauss(1)

  do i = 1, nr - 1
!-- first element along flux surface sufficient to determine ps
    nelm = (i - 1) * (np - 1) + 1
    call set_node_number(nelm, n1, n2, n3, n4)
    call interpolation(2, psi(4 * (n1 - 1) + 1), psi(4 * (n2 - 1) + 1), &
     psi(4 * (n3 - 1) + 1), psi(4 * (n4 - 1) + 1), r, s, ps, psr, pss, &
     psrs, psrr, psss)
!-- elements in HELENA are ordered in inverse radial direction, i.e.,
!   psikn(nr) = 0, psikn(1) = 1
    psi_new(nr + 1 - i) = ps
    if (hbt) then
      fdf_new(nr + 1 - i) = (cur0 * j_phi(ps)  - b * dp_dpsi(ps) &
       * flux_surface_average(nr + 1 - i, xaxis, 1._r8, 0._r8, &
       current_averaging, rquad_minus_1_over_2_r)) &
       / flux_surface_average(nr + 1 - i, xaxis, 1._r8, 0._r8, &
       current_averaging, one_over_r)
    else
      fdf_new(nr + 1 - i) = (cur0 * j_phi(ps)  - b * dp_dpsi(ps) &
       * flux_surface_average(nr + 1 - i, xaxis, 1._r8, 0._r8, &
       current_averaging, r_major)) &
       / flux_surface_average(nr + 1 - i, xaxis, 1._r8, 0._r8, &
       current_averaging, one_over_r)
    end if
  end do

!-- extrapolate fdf_new on axis
  psi_new(1) = 0_r8
  fdf_new(1) = (psi_new(3) - psi_new(1)) / (psi_new(3) - psi_new(2)) &
   * fdf_new(2) - (psi_new(2) - psi_new(1)) / (psi_new(3) - psi_new(2)) &
   * fdf_new(3)
!-- renomalize fdf_new to 1 on axis
  fdf_new = fdf_new / fdf_new(1)

!-- calculate new equidistant gamma profile fdf based on fdf_new
  call allocate_spline_coefficients(fdf_spline, nr)
  call spline(nr, psi_new, fdf_new, 0._r8, 0._r8, 2, fdf_spline)
  write(*, *) 'before evaluate'
  call evaluate_spline(fdf_spline, psi_new, psivec, fdf, error)
  write(*, *) 'after evaluate'
  call deallocate_spline_coefficients(fdf_spline)

  return
end subroutine flxint