suydam_read_helena.f90

Go to the documentation of this file.

subroutine suydam_read_helena(equilibrium)
!-----------------------------------------------------------------------
!    - reads helena mapping file
!    - scales quantities with q on axis
!-----------------------------------------------------------------------

  use itm_types
  use mod_map
  use mod_mesh, only : nr
  use mod_parameter
  use mod_profiles
  use mod_spline
  use mod_suydam
  use mod_helena_io, only : out_he
  use mod_output
  use helena_spline

  use euitm_schemas

  implicit none

  type(type_equilibrium), intent(in) :: equilibrium

  integer(itm_i4) :: i, j
  real(r8) :: c1(nr - 1), dummy(3)
  real(r8) :: scaleq, rbphi02, wurzel
  real(r8) :: dq0, dq1
  real(r8) :: dp_axis, dp_boundary
  real(r8) :: drbphi_axis, drbphi_boundary
  real(r8) :: raxis

  npsi = nr

!-- raxis = Rmag/Rvac
  raxis = equilibrium%global_param%mag_axis%position%r &
   / equilibrium%eqgeometry%geom_axis%r
  g11_hel(1, :) = 0._r8
  g33_hel(1, :) = 1._r8 / raxis**2

  call allocate_spline_coefficients(q_spline, npsi - 1)

  do j = 1, nchi
    c1(:) = g12_hel(2:npsi - 1, j)
    call spline(npsi - 1, cs(2), c1, 0._r8, 0._r8, 2, q_spline)
    g12_hel(1, j) = spwert(npsi - 1, 0._r8, q_spline, cs(2), dummy, 0)
  end do

  call deallocate_spline_coefficients(q_spline)

!------------------------------------------------------------------
!     scale quantities with value of q on axis (total current)
!------------------------------------------------------------------
  scaleq = qs(1) / qaxis

  cpsurf = cpsurf * scaleq
  p0 = p0 * scaleq**2
  g12_hel = -g12_hel * scaleq
  g11_hel = g11_hel * scaleq**2

  rbphi02 = rbphi(1)**2

!-- derivatives at plasma boundaries (all zero except for dp_boundary)
  dp_axis = 0._r8
  drbphi_axis = 0._r8
  drbphi_boundary = 0._r8
  dp_boundary = equilibrium%profiles_1d%pprime(nr) &
   * 2._r8 * equilibrium%profiles_1d%psi(nr) * mu0 &
   / equilibrium%global_param%mag_axis%bphi**2

!-- scaling
  dp_axis = dp_axis * scaleq**2
  dp_boundary = dp_boundary * scaleq**2
  drbphi_axis = drbphi_axis * scaleq**2
  drbphi_boundary = drbphi_boundary * scaleq / sqrt(1._r8 + rbphi02 &
   / rbphi(npsi)**2 * (1._r8 / scaleq**2 - 1._r8))

  do i = 1, npsi
    wurzel = sqrt(1._r8 + rbphi02 / rbphi(i)**2 * (1._r8 / scaleq**2 &
      - 1._r8))
     qs(i) = wurzel * qs(i)
     rbphi(i) = wurzel * scaleq * rbphi(i)
  end do

  if (standard_output) then
    write(out_he, *)  npsi, nchi
    write(out_he, 51) scaleq
    write(out_he, 52) cpsurf
    write(out_he, 53) (qs(i), i = 1, nr)
    write(out_he, 54) (p0(i), i = 1, npsi)
    write(out_he, 55) (rbphi(i), i = 1, npsi)
  end if

!     splines

  dq1 = (qs(npsi) - qs(npsi - 1)) / (cs(npsi) - cs(npsi - 1))
  dq0 = (qs(2) - qs(1)) / (cs(2) - cs(1))

  call allocate_spline_coefficients(q_spline, npsi)
  call allocate_spline_coefficients(p_spline, npsi)
  call allocate_spline_coefficients(rbphi_spline, npsi)

  call spline(npsi, cs, qs, dq0, dq1, 1, q_spline)
  call spline(npsi, cs, p0, dp_axis, dp_boundary, 1, p_spline)
  call spline(npsi, cs, rbphi, drbphi_axis, drbphi_boundary, 1, rbphi_spline)

  dqs = q_spline%sp2

  return

 51 format(/' after scale: scaleq=', es12.4)
 52 format(/' cpsurf = ', es12.4)
 53 format(/' qs'/(1x, 5es16.8))
 54 format(/' p0'/(1x, 5es16.8))
 55 format(/' rbphi'/(1x, 5es16.8))

end subroutine suydam_read_helena