safety_factor.f90

Go to the documentation of this file.

subroutine safety_factor(xaxis, a, cx, cy, qprf, q95out, q1out)
!----------------------------------------------------------------------
! subroutine to evaluate the q-profile
!----------------------------------------------------------------------
  
  use itm_types
  use mod_parameter
  use mod_dat
  use mod_gauss
  use mod_mesh
  use mod_profiles
  use mod_interpolation
  use phys_functions
  use mod_helena_io, only : iu6
  use helena_spline

  implicit none
  
  type(spline_coefficients) :: dq_spline
  real(r8), intent(in) :: xaxis, a, cx, cy
  real(r8), dimension(npts), intent(out) :: qprf
  real(r8), intent(out) :: q95out, q1out
  
  real(r8), dimension(nr) :: p0tmp, f2tmp, dptmp, dftmp
  real(r8), dimension(nr) :: zps, qq, avc
  real(r8), dimension(3) :: abltg
  
  integer(itm_i4) :: i, j, k
  integer(itm_i4) :: n1, n2, n3, n4
  real(r8) :: factas, sumq, avcur, xlength
  real(r8) :: r, s, ws, x, xr, xs, y, yr, ys, ps, psr, pss
  real(r8) :: xjac, bigr
  real(r8) :: arhs
  real(r8) :: ss
  
  if (ias == 1) then
    factas = 1._r8
  else
    factas = 2._r8
  endif
  if (verbosity > 2) write(iu6, *) 'safety factor' 
  call profiles(p0tmp, f2tmp, dptmp, dftmp, a)
  
!------------------------------------- nelm elements ------------------
  r = -1._r8
  do i = 1, nr - 1
    sumq = 0._r8
    avcur = 0._r8
    xlength = 0._r8
    do j = 1, np - 1
      n1 = (i - 1) * np + j
      n2 = n1 + 1
      n3 = n2 + np
      n4 = n1 + np
!------------------------------------- 4 point gaussian int. in s -----
      do k = 1, 4
        s = xgauss(k)
        ws = wgauss(k)
        call interpolation(1, xx(1, n1), xx(1, n2), xx(1, n3), xx(1, n4), &
         r, s, x, xr, xs)
        call interpolation(1, yy(1, n1), yy(1, n2), yy(1, n3), yy(1, n4), &
         r, s, y, yr, ys)
        call interpolation(1, 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)
        xjac = xr * ys - xs * yr
        bigr = 1._r8 + eps * x
        sumq = sumq - ws * xjac / (bigr * abs(psr))
        xlength = xlength + sqrt(xs**2 + ys**2) * ws
        arhs = rh_side(1._r8, x, 0._r8, 0._r8, 0._r8, 0._r8, ps, 0._r8, 0._r8)
        avcur = avcur + ws * arhs * sqrt(xs**2 + ys**2)
      end do
    end do
    zps(nr - i + 1) = ps
    qq(nr - i + 1) = factas / 2._r8 * f2tmp(nr - i + 1) * sumq * alfa / pi
    avc(nr - i + 1) = factas * a * eps * avcur / (alfa * xlength)
  end do
!------------------------------------- calc equidistant q profile   
  qq(1) = f2tmp(1) * alfa / (2._r8 * sqrt(cx * cy) * (1._r8 + eps &
   * xaxis))
  zps(1) = 0._r8
  call allocate_spline_coefficients(dq_spline, nr)
  call spline(nr, zps, qq, 0._r8, 0._r8, 2, dq_spline)
  do i = 1, npts
    ss = dble(i - 1) / dble(npts - 1)
    ps = ss * ss
    qprf(i) = spwert(nr, ps, dq_spline, zps, abltg, 0)
  end do 
  ps = 0.95_r8
  q95out = spwert(nr, ps, dq_spline, zps, abltg, 0)
  q1out = qq(nr)
  if (verbosity > 2) write(iu6, *) ' q95 : ', q95out
!------------------------------------- calc equidistant current profile   
  avc(1) = avc(2) - (avc(3) - avc(2)) / (zps(3) - zps(2)) * zps(2)
  call spline(nr, zps, avc, 0._r8, 0._r8, 2, dq_spline)
  do i = 1, npts
    ss = dble(i - 1) / dble(npts - 1)
    ps = ss * ss
    j_tor(i) = spwert(nr, ps, dq_spline, zps, abltg, 0)
  end do

  call deallocate_spline_coefficients(dq_spline)

  return
end subroutine safety_factor