diagnostics.f90

Go to the documentation of this file.

subroutine diagnostics(a, xaxis, yaxis, zvol, zvolp, current)
!----------------------------------------------------------------------
! subroutine to evaluate output quantities
!----------------------------------------------------------------------

  use itm_types
  use mod_parameter
  use mod_corners
  use mod_dat
  use mod_gauss
  use mod_mesh
  use mod_output
  use mod_helena_io, only : out_he
  use mod_interpolation
  use mod_output
  use phys_profiles, only : j_phi
  use mod_parameter
  use mod_profiles
  use mod_spline
  use matrix_calculations
  use phys_functions, only : rh_side
  use helena_spline

  implicit none

  real(r8), intent(in) :: a, xaxis, yaxis
  real(r8), dimension(nr), intent(out) :: zvol, zvolp
  real(r8), intent(out) :: current

  real(r8), dimension(nr) :: avc
  real(r8), dimension(nr - 1) :: zjar, zpar, zar, zps, xl
  real(r8), dimension(3) :: ablt
  real(r8) :: factas
  real(r8) :: area, volume
  real(r8) :: parea, carea, parea2
  real(r8) :: grps2, bp2, bp2vol
  real(r8) :: r, wr, s, ws
  real(r8) :: x, xr, xs, y, yr, ys, ps, psr, pss
  real(r8) :: xjac, arhs
  real(r8) :: rav, xli
  real(r8) :: beta, betapl, betastar
  real(r8) :: xlength, avcur, sumc3, sumc5
  real(r8) :: bigr, zjdchi, dl, dlv
  real(r8) :: prs, bussac, err
  integer(itm_i4) :: i, j, ngr, ngs
  integer(itm_i4) :: nelm
  integer(itm_i4) :: n1, n2, n3, n4

  factas = 2._r8
  if (ias == 1) factas = 1._r8

  area = 0._r8
  volume = 0._r8
  bp2vol = 0._r8
  parea = 0._r8
  carea = 0._r8
  parea2 = 0._r8

!-- spline p_int and gam_int with use of known derivatives
  call allocate_spline_coefficients(p_spline, npts)
  call spline(npts, psivec, p_int, dpres(1) * sign(1._r8, cpsurfin), &
   dpres(npts) * sign(1._r8, cpsurfin), 1, p_spline)

!------------------------------------- nelm elements ------------------
  do i = nr - 1, 1, -1
    do j = 1, np - 1
      nelm = (i - 1) * (np - 1) + j
      call set_node_number(nelm, n1, n2, n3, n4)
!------------------------------------- 4 point gaussian int. in r -----
      do ngr = 1, 4
        r = xgauss(ngr)
        wr = wgauss(ngr)
!------------------------------------- 4 point gaussian int. in s -----
        do ngs = 1, 4
          s = xgauss(ngs)
          ws = wgauss(ngs)
          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
          area = area + wr * ws * xjac
          volume = volume + (1. + eps * x) * wr * ws * xjac
          grps2 = psr**2 * (xs**2 + ys**2) / xjac**2
          bp2 = grps2 / (1._r8 + eps * x)**2    
          bp2vol = bp2vol + (1._r8 + eps * x) * bp2 *wr * ws * xjac
          parea = parea - wr * ws * xjac * spwert(npts, ps, p_spline, &
           psivec, ablt, 0)
          parea2 = parea2 + wr * ws * xjac * spwert(npts, ps, p_spline, &
           psivec, ablt, 0)**2
          arhs = rh_side(1._r8, x, 0._r8, 0._r8, 0._r8, 0._r8, ps, 0._r8, &
           0._r8)
          carea = carea + wr * ws * xjac * arhs
        end do
      end do
    end do
    if (hbt) then
      zjar(i) = -factas * a * carea 
      zpar(i) = -factas * 0.5_r8 * a * b * parea 
    else  
      zjar(i) = -factas * a * carea 
      zpar(i) = -factas * eps * a * b * parea
    end if
    zar(i)  = factas * abs(area)
    zvol(nr - i + 1) = factas * abs(volume)
  end do
  area = factas * abs(area)
  volume = factas * abs(volume)
  rav = volume / area
  if (hbt) then
    parea = factas * 0.5_r8 * a * b * parea
    parea2 = factas * (0.5_r8*a*b)**2 * parea2
    carea = factas * a * carea
  else
    parea = factas * eps * a * b * parea
    parea2 = factas * (eps * a * b)**2 * parea2
    carea = factas * a * carea
  end if

  volume = 2._r8 * pi * volume
  xli = abs(4._r8 * factas * pi * bp2vol / (rav * carea**2))
  beta = 2._r8 * (eps / alfa**2) * abs(parea) / area
  betastar = 2._r8 * (eps / alfa**2) * sqrt(abs(parea2) / area)
  betapl = (8._r8 * pi / eps) * abs(parea) / carea**2
  current = abs(eps * carea / alfa)

!-- denormalize current
  current = current * rvac * bvac * eps / (4.e-7_r8 * pi)

  if (standard_output) then
    write(out_he, *)
    write(out_he, *) '***************************************'
    write(out_he, 11) xaxis, yaxis
    write(out_he, 2) betapl
    write(out_he, 3) beta
    write(out_he, 8) betastar
    write(out_he, 12) 1.256637_r8 * beta / current * 1.e+6_r8
    write(out_he, 4) current
    write(out_he, 5) area
    write(out_he, 6) volume
    write(out_he, 7) xli
    write(out_he, 9) alfa
    write(out_he, 13) a, b
    write(out_he, *) '***************************************'
    write(out_he, *)
  end if

  write(txtout(20), 2) betapl
  write(txtout(21), 3) beta
  write(txtout(22), 8) betastar
  write(txtout(23), 4) current
  write(txtout(24), 5) area
  write(txtout(25), 6) volume
  write(txtout(26), 7) xli
  write(txtout(27), 9) alfa

!------------------------------------- nelm elements ------------------
  r = -1._r8
  do i = 1, nr - 1
    xlength = 0._r8
    avcur = 0._r8
    sumc3 = 0._r8
    sumc5 = 0._r8
    do j = 1, np - 1
      nelm = (i - 1) * (np - 1) + j
      call set_node_number(nelm, n1, n2, n3, n4)
!------------------------------------- 4 point gaussian int. in s -----
      do ngs = 1, 4
        s = xgauss(ngs)
        ws = wgauss(ngs)
        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)                                         
        zjdchi = bigr * xjac / abs(psr)  
        dl = sqrt(xs**2 + ys**2)
        xlength = xlength +  dl * ws
        sumc5 = sumc5 + ws * zjdchi
        if (trim(current_averaging) == 'line') then
          dlv = dl
        else
          dlv = zjdchi
        endif
        sumc3 = sumc3 + ws * dlv
        arhs = rh_side(1._r8, x, 0._r8, 0._r8, 0._r8, 0._r8, ps, 0._r8, 0._r8)
        avcur = avcur + ws * arhs * dlv
      end do
    end do
    zps(i) = ps
    xl(i) = factas * xlength
    avc(i) = factas * a * eps * avcur / (alfa * sumc3)
    zvolp(nr - i + 1) = factas * abs(sumc5)
  end do

  if (standard_output) then
    write(out_he, *)
    write(out_he, *) '***********************************************' &
               //'********'
    write(out_he, *) '* i   psi     s      <j>     error   length    ' &
               //'bussac *'
    write(out_he, *) '***********************************************' &
               //'********'
    avc(nr) = avc(nr - 1) - (avc(nr - 2) - avc(nr - 1)) &
     / (zps(nr - 2) - zps(nr - 1)) * zps(nr - 1)
  
    do i = nr - 1, 1, -1
      ps = zps(i)
      if (hbt) then
        prs = 0.5 * a * b * spwert(npts, ps, p_spline, psivec, &
         ablt, 0)
      else
        prs = eps * a * b * spwert(npts, ps, p_spline, psivec, &
         ablt, 0)
      end if
      bussac = -(2._r8 / eps) * (prs - zpar(i) / zar(i)) / (zjar(i) &
       / xl(i))**2
      err = avc(i) / avc(nr) - j_phi(zps(i))
      write(out_he, 72) i, zps(i), sqrt(zps(i)), avc(i) / avc(nr), &
       err, xl(i), bussac
    end do
    write(out_he, *) '***********************************************' &
               //'********'
    write(out_he, *)
  end if


  2 format('  poloidal beta : ', e12.4)
  3 format('  toroidal beta : ', e12.4)
  4 format('  total current : ', e12.4)
  5 format('  total area    : ', e12.4)
  6 format('  total volume  : ', e12.4)
  7 format('  int. inductance : ', e12.4)
  8 format('  beta star     : ', e12.4)
  9 format('  pol. flux     : ', e12.4)
 11 format('  magnetic axis : ', 2f9.5)
 12 format('  norm. beta    : ', f9.5)
 13 format('  a, b          : ', 2e12.4)
 72 format(i3, 3f12.4, es10.2, 2f12.4)   

  call deallocate_spline_coefficients(p_spline)

  return
end subroutine diagnostics