si_output.f90

Go to the documentation of this file.

subroutine si_output(a, r0, b0, fcirc, qprof, dqprof, rav, geonc, &
 zjpar, dimerc, drmerc, hh, current)
!----------------------------------------------------------------------
! subroutine to translate back to real world values using the major
! radius (m) , vacuum magnetic field (T) at the geometric centre,
! and the central density (10^19 m-3), rpe = te/(te+ti)
! calculates bootstrap current profile.
!----------------------------------------------------------------------

  use itm_constants
  use mod_parameter
  use mod_corners
  use mod_dat
  use mod_gauss
  use mod_mesh
  use mod_dete
  use mod_helena_io, only : out_he
  use mod_interpolation
  use mod_output
  use plot_data
  use matrix_calculations
  use phys_functions
  use mod_ne_te_profile

  implicit none

  real(r8), parameter :: zmu0 = 4e-7_r8 * pi
  real(r8), parameter :: ze = itm_qe

  real(r8), dimension(nr), intent(in) :: fcirc, qprof, dqprof, rav, &
   geonc, dimerc, drmerc, hh
  real(r8), dimension(nr), intent(inout) :: zjpar
  real(r8), intent(in) :: a, r0, b0, current

  real(r8), dimension(nr) :: p0tmp, ftmp, dptmp, dftmp
  real(r8), dimension(nr) :: sr, avc, dqdpsi, dnc, stab, sp
  real(r8), dimension(nr) :: zni, tte, ti
  real(r8), dimension(nr) :: zjbt, hjbt, signeo, sigspitz, znue, znui
  real(r8), dimension(3) :: qp
  real(r8) :: factas
  real(r8) :: radius, cpsurf
  real(r8) :: rbscale
  real(r8) :: dndpsico, dtdpsico
  real(r8) :: defte, tefte
  real(r8) :: conste, consde
  real(r8) :: ss
  real(r8) :: pwr
  real(r8) :: dni, dtte, dti
  real(r8) :: rr
  real(r8) :: zznue, zznui, zsigsp, zsigneo, zzjbt, zhjbt
  real(r8) :: r, s, ws
  real(r8) :: xlength, avcur, xjac
  real(r8) :: x, xr, xs, y, yr, ys, ps
  real(r8) :: arhs
  real(r8) :: alfas, alfal, ef
  real(r8) :: dlt, ssp, dqp, dnp, dip, drp, hp, zjp, efp, ef2, shrp, dsp
  integer(itm_i4) :: i, j
  integer(itm_i4) :: nelm
  integer(itm_i4) :: n1, n2, n3, n4
  integer(itm_i4) :: ngs, nqp
  integer(itm_i4), parameter :: out_plot_te = 33, out_plot_ne = 34

!---- initializations
  zjbt = 0._r8
  hjbt = 0._r8
  signeo = 0._r8
  sigspitz = 0._r8
  znue = 0._r8
  znui = 0._r8

  factas = 2._r8
  sr(1) = 0._r8
  sr(nr) = 1._r8
  if (ias == 1) factas = 1._r8

  radius = eps * r0
  cpsurf = radius**2 * b0 / alfa

  call profiles(p0tmp, ftmp, dptmp, dftmp, a)

  pscale = b0**2 * eps / (zmu0 * alfa**2)
  rbscale = b0 * r0

  do i = 1, nr
    p0tmp(i) = p0tmp(i) * pscale
    dptmp(i) = dptmp(i) * pscale
    ftmp(i) = ftmp(i) * rbscale
    dftmp(i) = dftmp(i) * rbscale
  end do

  if (standard_output) then
    write(out_he, *) 
    write(out_he, *) '*************************************************'
    write(out_he, *) ' real world output :'
    write(out_he, *) '*************************************************'
    write(out_he, 1)  r0
    write(out_he, 2)  b0
    write(out_he, 10) current / 1000._r8
    write(out_he, 9)  radius
    write(out_he, 8)  cpsurf
    write(out_he, 3)  zn0
    write(out_he, 4)  zeff
    write(out_he, 6)  rpe 
    write(out_he, *) '*************************************************'

    write(out_he, *) '**************************************************'
    write(out_he, *) '   s,   p [pa], ne [10^19m^-3], te [eV],  ti [eV],'
    write(out_he, *) '**************************************************'
  end if

  if (additional_output) then
    write(out_plot_te, *) '\psi   T_e  ', nr - 2 
    write(out_plot_ne, *) '\psi   n_i  ', nr - 2
  end if

  if (te%shape > 0) then
    call ne_te_profile(ne, ne%psi_0, defte, dndpsico, consde)
    call ne_te_profile(te, te%psi_0, tefte, dtdpsico, conste)
  end if
  do i = 2, nr - 1
    ps = psi(4 * (nr - i) * np + 1)
    ss = sqrt(ps)

    if (te%shape == 0) then
      pwr = 1._r8 / (etaei + 1._r8)
      zni(i) = zn0 * (p0tmp(i) / p0tmp(1))**pwr
      dni = zn0 * pwr * abs((p0tmp(i) / p0tmp(1)))**(pwr - 1._r8) &
       * dptmp(i) / p0tmp(1) / (2._r8 * ss)
      tte(i) = rpe * p0tmp(i) / (ze * zni(i) * 1e19_r8)
      ti(i) = (1._r8 - rpe) * p0tmp(i) / (ze * zni(i) * 1e19_r8)
      dtte = rpe / (ze * 1e19_r8) * (dptmp(i) / (2._r8 * ss * zni(i)) &
       - p0tmp(i) * dni / zni(i)**2)
      dti = (1._r8 - rpe) / (ze * 1e19_r8) * (dptmp(i) / (2._r8 &
       * ss * zni(i)) - p0tmp(i) * dni / zni(i)**2)  
    else
      if (ps >= ne%psi_0) then
        call ne_te_profile(ne, ps, zni(i), dni)
      else
        call ne_te_profile(ne, ps, zni(i), dni, consde, dndpsico)
      end if
      dni = -dni
!     te given in eV (profiles defined in keV)
      if (ps >= te%psi_0) then
        call ne_te_profile(te, ps, tte(i), dtte)
      else
        call ne_te_profile(te, ps, tte(i), dtte, conste, dtdpsico)
      end if              
      dtte = -dtte * 1000.0_r8
      ti(i) = tte(i)
      dti = dtte
    end if
    if (additional_output) then
      write (out_plot_te, *) ps, tte(i), dtte
      write (out_plot_ne, *) ps, zni(i), dni
    end if
    dtte = dtte / cpsurf
    dti = dti / cpsurf
    dni = dni / cpsurf

    rr = r0 * rav(i)

! TODO: te, ti in eV but dte, dti in keV ???

    call neo_helena(fcirc(i), tte(i), ti(i), zni(i), zeff, ftmp(i), qprof(i), &
     rr, eps * ss, dni, dtte, dti, zznue, zznui, zsigsp, zsigneo, &
     zzjbt, zhjbt)

    zjbt(i) = zzjbt
    hjbt(i) = zhjbt
    signeo(i) = zsigneo
    sigspitz(i) = zsigsp
    znue(i) = zznue
    znui(i) = zznui
    sr(i) = ss

    zjpar(i) = zjpar(i) / (zmu0 * r0 / b0**2)

    if (standard_output) &
     write (out_he, 7) sqrt(ps), p0tmp(i), zni(i), tte(i), ti(i), &
     zjbt(i), hjbt(i), hjbt(i) * zmu0 * r0 / b0**2, zjpar(i)
    if (additional_output) then
      write (out_plot_ne+10,*) sqrt(ps), zni(i), dni
      write (out_plot_te+10,*) sqrt(ps), tte(i), dtte
    end if
  end do  
!------------------------------- plot neoclassical quantities
  if (xmgrace_output) then
    call plot_data_1d('line', sr, zjbt, nr, 'xmgr_zjbt_rhopol')
    call plot_data_1d('line', sr, hjbt, nr, 'xmgr_hjbt_rhopol')
    call plot_data_1d('line', sr, signeo, nr, 'xmgr_signeo_rhopol')
    call plot_data_1d('line', sr, sigspitz, nr, 'xmgr_sigspitz_rhopol')
    call plot_data_1d('line', sr, znue, nr, 'xmgr_znue_rhopol')
    call plot_data_1d('line', sr, znue, nr, 'xmgr_znui_rhopol')
    call plot_data_1d('line', sr, zjpar, nr, 'xmgr_zjpar_rhopol')
  end if
!------------------------------- flux averaged current profile      
  r = -1._r8
  do i = 2, nr - 1
    xlength = 0._r8
    avcur = 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(0, psi(4 * (n1 - 1) + 1), psi(4 * (n2 - 1) + 1), &
         psi(4 * (n3 - 1) + 1), psi(4 * (n4 - 1) + 1), r, s, ps)
        xjac = xr * ys - xs * yr
        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
    xlength = factas * xlength
    avc(nr - i + 1) = factas * a * eps * avcur/ (alfa * xlength)
    avc(nr - i + 1) = avc(nr - i + 1) * b0 / (zmu0 * r0)
  end do
  if (standard_output) then
    write(out_he, *)
    write(out_he, *) '************************************************', &
     '***********************'
    write(out_he, *) '*   s,       q,      fcirc,     nu_e,     nu_i, ',  &
     '  sig(spitz),sig(neo) *'
    write(out_he, *) '************************************************', &
     '***********************'
    do i = 2, nr - 1
      write(out_he, 7) sr(i), qprof(i), fcirc(i), znue(i), znui(i), &
       sigspitz(i), signeo(i)
    end do
  end if

! here previously plotted: p0temp, te, ti, zni, znue, znui, sigspitz,
!  zjbt, hjbt, zjpar = pressure, temperatures, density, collisionalities,
!  conductivities, bootstrap current, current profile

!----------------------- stability condition Hegna (Phys. Plas. 1999)
  if (standard_output) then
    write(out_he, *)
    write(out_he, *) '*****************************************'
    write(out_he, *) '*  neoclassical stability condition     *'
    write(out_he, *) '*****************************************'
  end if

  j = 0

!-- initialize variables
  dnc = 0._r8

  do i = 2, nr - 1
!----------------------- undo scaling, i.e. back to dimensionless
    p0tmp(i) = p0tmp(i) * zmu0 / b0**2
    dptmp(i) = dptmp(i) * zmu0 / b0**2
    ftmp(i) = ftmp(i) / rbscale
    dftmp(i) = dftmp(i) / rbscale

    zjbt(i) = zjbt(i) * zmu0 * r0 / b0**2
    dqdpsi(i) = dqprof(i) / (2._r8 * sr(i) * cpsurf / (r0**2 * b0) )
    if (dqdpsi(i) /= 0._r8) then
      dnc(i) =  zjbt(i) / dqdpsi(i) * geonc(i)
    else
      dnc(i)  =  0._r8
    end if
    if (dimerc(i) > 0._r8) then
      j = j + 1
      alfas = 0.5_r8 + sqrt(dimerc(i))
      alfal = 0.5_r8 - sqrt(dimerc(i))
      ef = -drmerc(i) - hh(i)**2
      stab(j) = dnc(i) + (-drmerc(i)) / (alfas - hh(i))
      sp(j) = sr(i)
      if (standard_output) &
       write(out_he, 101) sr(i), qprof(i), zjbt(i), dimerc(i), dnc(i), &
       stab(j)
    else
      if (standard_output) &
       write(out_he, 102) sr(i), qprof(i), zjbt(i), dimerc(i), dnc(i)
    end if
  end do
!
!----------------------- write data for analysis of Hegna's paper to file
!
  if (additional_output) then
    write(41, *)  ' r0, b0, cpsurf, eps :'
    write(41, 41) r0, b0, cpsurf, eps
    write(41, *)  ' i, s, p0, dp0, f, df, q, dqdpsi, j.b, d_nc : '
    do i = 2, nr - 1
      write(41, 42) i, sr(i), p0tmp(i), dptmp(i), ftmp(i), dftmp(i), &
       qprof(i), dqdpsi(i), zjbt(i), dnc(i)
    end do
  end if
  nqp = 3
  qp(1) = 1.5_r8
  qp(2) = 2.0_r8
  qp(3) = 3.0_r8
  if (standard_output) then
    write(out_he, *)
    write(out_he, *) '     sp,     q,     shear,      d_i,      d_r,', &     
     '    e+f,     h,      d_nc,       dsp,      zjbt'
  end if
  do i = 1, nqp
    do j = 2, nr - 1
      if (qp(i) > qprof(j) .and. qp(i) < qprof(j+1)) then
        dlt = (qp(i) - qprof(j)) / (qprof(j + 1) - qprof(j))
        ssp = sr(j) + dlt * (sr(j + 1) - sr(j))
        dqp = dqprof(j) + dlt * (dqprof(j + 1) - dqprof(j))
        dnp = dnc(j) + dlt * (dnc(j + 1) - dnc(j))
        dip = dimerc(j) + dlt * (dimerc(j + 1) - dimerc(j))
        drp = drmerc(j) + dlt * (drmerc(j + 1) - drmerc(j))
        hp  = hh(j) + dlt * (hh(j + 1) - hh(j))
        zjp = zjbt(j) + dlt * (zjbt(j + 1) - zjbt(j))
        efp = -drp - hp**2
        ef2 = -dip - hp + 0.25_r8
        shrp = dqp / qp(i) * ssp
        dsp = 0._r8
        if (dip >= 0._r8) dsp = 0.5_r8 - sqrt(dip) - hp  
        if (standard_output) &
         write(out_he, 43) ssp, qp(i), shrp, -dip, -drp, efp, ef2, &
         hp, dnp, dsp, zjp
       end if
     end do
   end do

 1 format(' major radius    : ', f9.4,' [m]')
 2 format(' magnetic field  : ', f9.4,' [T]')
 3 format(' central density : ', f9.4,' 10^19 [m^-3]')
 4 format(' zeff            : ', f9.4)
 6 format(' te/(te+ti)      : ', f9.4)
 7 format(12es11.4)
 8 format(' psi on boundary : ', f9.4,' [wb/rad]')
 9 format(' radius (a)      : ', f9.4,' [m]')
 10 format(' total current   : ', f10.3,' [kA]')
 41 format(12es12.4)
 42 format(i4, 12es12.4)
 43 format(12f8.4)
 102 format(' s =', f7.3, ' q =', f7.3, ' j.b=', f7.3, ' d_i =', &
   es10.3, ' d_nc =', e10.3, ' mercier unstable')
 101 format(' s =', f7.3, ' q =', f7.3, ' j.b=',f7.3, ' d_i =', &
   es10.3, ' d_nc =', es10.3, ' d_nc+d_r/(a-h) =', es10.3)

  return
end subroutine si_output