profiles.f90

Go to the documentation of this file.

subroutine profiles(p0, rbphi, dp0, drbphi, a)
!-----------------------------------------------------------------------
! subroutine to integrate the equilibrium profiles
! using the hbt definitions
!-----------------------------------------------------------------------

  use itm_types
  use mod_dat
  use mod_mesh
  use mod_nodes
  use mod_helena_io, only : iu6
  use mod_profiles
  use mod_spline
  use helena_spline

  implicit none

  integer(itm_i4) :: j

  real(r8), dimension(nr), intent(out) :: p0, rbphi, dp0, drbphi
  real(r8), intent(in) :: a
  real(r8) :: flux

  real(r8), dimension(3) :: ablt

  real(r8) :: dp_dpsi, dgamma_dpsi

!-- spline p_int and gam_int with use of known derivatives
  call allocate_spline_coefficients(p_spline, npts)
  call allocate_spline_coefficients(rbphi_spline, npts)

!-- sign of cpsurfin needed because psivec always positive
  call spline(npts, psivec, p_int, dpres(1) * sign(1._r8, cpsurfin), &
   dpres(npts) * sign(1._r8, cpsurfin), 1, p_spline)
  call spline(npts, psivec, gam_int, dgam(1) * sign(1._r8, cpsurfin), &
   dgam(npts) * sign(1._r8, cpsurfin), 1, rbphi_spline)

!---------------------------- derivatives to r = sqrt(psi) !!!! --------
!  multiplication by sign(-alfa) * sign(bvac) compensates for loss of sign
!  due to normalization with pscale and fscale
  do j = 1, nr
    flux = psikn(nr - j + 1)
    if (hbt) then
      p0(j) = 0.5_r8 * a * b * spwert(npts, flux, p_spline, psivec, ablt, 0) &
       * sign(1.0_r8, -alfa) * sign(1.0_r8, bvac)
      dp0(j) = -a * b * dp_dpsi(flux) * sqrt(flux)
      rbphi(j) = p0(j) + eps * a * spwert(npts, flux, rbphi_spline, psivec, &
       ablt, 0) * sign(1.0_r8, -alfa) * sign(1.0_r8, bvac)
      drbphi(j) = dp0(j) - eps *  2._r8 * a * dgamma_dpsi(flux) &
       * sqrt(flux)
    else
      p0(j) = eps * a * b * spwert(npts, flux, p_spline, psivec, ablt, 0) &
       * sign(1.0_r8, -alfa) * sign(1.0_r8, bvac)
      dp0(j) = -2._r8 * eps * a * b * dp_dpsi(flux) * sqrt(flux)
      rbphi(j) = eps * a * spwert(npts, flux, rbphi_spline, psivec, ablt, 0) &
       * sign(1.0_r8, -alfa) * sign(1.0_r8, bvac)
      drbphi(j) = -eps * 2._r8 * a * dgamma_dpsi(flux) * sqrt(flux)
    end if
!-- rbphi still missing sign(bvac), added in fill_equilibrium
    rbphi(j) = sqrt(1._r8 + 2._r8 * eps * rbphi(j) / alfa**2)
    drbphi(j) = 1._r8 / (2._r8 * rbphi(j)) * 2._r8 * eps * drbphi(j) / alfa**2
  end do

  if (verbosity > 2) write(iu6, *) 'profiles computed with A = ', a

  call deallocate_spline_coefficients(p_spline)
  call deallocate_spline_coefficients(rbphi_spline)

  return
end subroutine profiles