update_fdf.f90

Go to the documentation of this file.

subroutine update_fdf(equilibrium, xaxis)
!----------------------------------------------------------------------
! subroutine to calculate FF' from <j_tor> and p' (in SI)
!----------------------------------------------------------------------
  
  use itm_types
  use mod_parameter, only : mu0, twopi
  use mod_dat
  use mod_helena_io, only : iu6
  use mod_mesh, only : nr
  use mod_profiles
  use helena_spline

  use phys_functions
  use math_functions, only : flux_surface_average
  use euitm_schemas

  implicit none
  
  type (type_equilibrium), intent(inout) :: equilibrium

  real(r8), intent(in) :: xaxis

  type(spline_coefficients) :: fdf_spline

  real(r8), dimension(nr) :: psi_new, fdf_new
  real(r8), dimension(3) :: abltg
  integer(itm_i4) :: i

  if (verbosity > 1) write(iu6, *) "updating FF' based on <j_tor> and p'"

!-- calculate new FF' from flux surface averaged Grad-Shafranov equation
!-- Attention: FF' and p' in HELENA are with respect to Psi in [Wb].
!   In Grad-Shafranov equation, however, with respect to Psi [Wb/rad].
  do i = 1, nr
    fdf_new(i) = mu0 * (equilibrium%profiles_1d%jphi(i) &
     - flux_surface_average(i, xaxis, 1._r8, 0._r8, current_averaging, &
     r_major) * equilibrium%eqgeometry%geom_axis%r * twopi &
     * equilibrium%profiles_1d%pprime(i)) &
     / flux_surface_average(i, xaxis, 1._r8, 0._r8, current_averaging, &
     one_over_r) * equilibrium%eqgeometry%geom_axis%r / twopi
  end do

!-- store as new FF'
  equilibrium%profiles_1d%ffprime = fdf_new

!-- calculate new gamma profile from p' and FF'
!   (relative normalization only, not SI)
  if (hbt) then
    fdf_new = mu0 * equilibrium%profiles_1d%pprime &
     + equilibrium%profiles_1d%ffprime &
     / equilibrium%eqgeometry%geom_axis%r**2
  else
    fdf_new = equilibrium%profiles_1d%ffprime
  end if

!-- renomalize fdf_new to 1 on axis
  fdf_new = fdf_new / fdf_new(1)

!-- normalized psi
  psi_new = equilibrium%profiles_1d%psi / equilibrium%profiles_1d%psi(nr)

!-- 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)
  do i = 1, npts
    fdf(i) = spwert(nr, psivec(i), fdf_spline, psi_new, abltg, 0)
  end do

  call deallocate_spline_coefficients(fdf_spline)

  return
end subroutine update_fdf