neo_helena.f90

Go to the documentation of this file.

subroutine neo_helena(fcirc, te, ti, zzne, z, f, q, r, eps, zdne, dte, &
 dti, znue, znui, sigspitz, signeo, zjbt, hjbt)
!-----------------------------------------------------------------------
! formula for neoclassical quantities from Phys. Plasmas 6 (1999) 2834.
! te, ti given in electron Volt, density im m^-3
!-----------------------------------------------------------------------
!TODO: check validity of this subroutine (results differ from Samuli's
!      original but may actually be correct)
  
  use itm_types

  implicit none
  
  real(r8), intent(in) :: fcirc
  real(r8), intent(in) :: zzne, zdne
  real(r8), intent(in) :: z, q, r, f, eps
  real(r8), intent(in) :: te, ti, dte, dti
  
  real(r8), intent(out) :: znue, znui
  real(r8), intent(out) :: sigspitz, signeo
  real(r8), intent(out) :: zjbt, hjbt
  
  real(r8) :: ft, rpe, znz
  real(r8) :: zne, zni, zdni, dni
  real(r8) :: zle, zli
  real(r8) :: f33tef, f31tef, f32tefe, f32tefi, f32ee, f32ei, f34tef
  real(r8) :: x, y
  real(r8) :: zl31, zl32, zl34
  real(r8) :: a0, anu, p
  real(r8) :: xx, dx, alfai, etai
  real(r8) :: a1, a2e, a2i
  real(r8) :: hl31, hl32
  
  ft = 1._r8 - fcirc
  
  zne = zzne * 1e19_r8
  zni = zne / z
  zdni = zdne / z
  dni = zdni * 1e19_r8
  
  rpe = te / (ti + te)
  
  zle = 31.3_r8 - log(sqrt(zne) / te)
  zli = 30.0_r8 - log(z**3 * sqrt(zni) / ti**1.5_r8)      
  
  
  znue = 6.921e-18_r8 * q * r * zne * z * zle / (te**2 * eps**1.5_r8)
  znui = 4.90e-18_r8 * q * r * zni * z**4 * zli / (ti**2 * eps**1.5_r8)

  f33tef = ft / (1._r8 + (0.55_r8 - 0.1_r8 * ft) * sqrt(znue) &
   + 0.45_r8 * (1._r8 - ft) * znue / z**1.5_r8)

  znz = 0.58_r8 + 0.74_r8 / (0.76_r8 + z)

  sigspitz = 1.9012e4_r8 * te**1.5_r8 / (z * znz * zle)

  x = f33tef
  signeo = 1._r8 - (1._r8 + 0.36_r8 / z) * x + 0.59_r8 / z * x &
   * x - 0.23_r8 / z * x**3
  signeo = signeo * sigspitz

  f31tef = ft / (1._r8 + (1._r8 - 0.1_r8 * ft) * sqrt(znue) &
   + 0.5_r8 * (1._r8 - ft) * znue / z)

  x = f31tef
  zl31 = (1._r8 + 1.4_r8 / (z + 1._r8)) * x - 1.9_r8 / (z &
   + 1._r8) * x * x + 0.3_r8 / (z + 1._r8) * x**3 + 0.2_r8 &
   / (z + 1._r8) * x**4

  f32tefe = ft / (1._r8 + 0.26_r8 * (1._r8 - ft) * sqrt(znue) &
   + 0.18_r8 * (1._r8 - 0.37_r8 * ft) * znue / sqrt(z))

  f32tefi = ft / (1._r8 + (1._r8 + 0.6_r8 * ft) * sqrt(znue) &
   + 0.85_r8 * (1._r8 - 0.37_r8 * ft) * znue * (1._r8 + z))

  x = f32tefe
  f32ee = (0.05_r8 + 0.62_r8 * z) / (z * (1._r8 + 0.44_r8 * z)) &
   * (x - x**4) + 1._r8 / (1._r8 + 0.22_r8 * z) * (x**2 - x**4 &
   - 1.2_r8 * (x**3 - x**4)) + 1.2_r8 / (1._r8 + 0.5_r8 * z) &
   * x**4

  y = f32tefi
  f32ei = -(0.56_r8 + 1.93_r8 * z) / (z * (1._r8 + 0.44_r8 * z)) &
   * (y - y**4) + 4.95_r8 / (1._r8 + 2.48_r8 * z) * (y**2 - y**4 &
   - 0.55_r8 * (y**3 - y**4)) - 1.2_r8 / (1._r8 + 0.5_r8 * z) &
   * y**4

  zl32 = f32ee + f32ei

  f34tef = ft / (1._r8 + (1._r8 - 0.1_r8 * ft) * sqrt(znue) &
   + 0.5_r8 * (1._r8 - 0.5_r8 * ft) * znue / z) 

  a0 = -1.17_r8 * (1._r8 - ft) / (1._r8 - 0.22_r8 * ft &
   - 0.19_r8 * ft * ft)

!----------------------------- sign in front of 0.315 changed !!!

  anu = ((a0 + 0.25_r8 * (1._r8 - ft * ft) * sqrt(znui)) / (1._r8 &
   + 0.5_r8 * sqrt(znui)) + 0.315_r8 * znui**2 * ft**6) / (1._r8 &
   + 0.15_r8 * znui**2 * ft**6)

  x = f34tef
  zl34 = (1._r8 + 1.4_r8 / (z + 1._r8)) * x - 1.9_r8 / (z &
   + 1._r8) * x * x + 0.3_r8 / (z + 1._r8) * x**3 + 0.2_r8 &
   / (z + 1._r8) * x**4

  p = (zni * ti + zne * te) * 1.602e-19_r8
  zjbt = -f * p * (zl31 * dni / zni + rpe * (zl31 + zl32) * dte &
   / te + (1._r8 - rpe) * (1._r8 + zl34 / zl31 * anu) * zl31 * dti / ti)

!------------------------------ Hirshman formula at low collisionality
  xx = (1._r8 - fcirc) / fcirc

  dx = 1.414_r8 * z + z * z + xx * (0.754_r8 + 2.657_r8 * z &
   + 2._r8 * z * z) + xx * xx * (0.348_r8 + 1.243_r8 * z + z * z)

  alfai = -1.172_r8 / (1._r8 + 0.462_r8 * xx)

  if (dni /= 0._r8) then
     etai = (dte / te) / (dni / zni)
  else
     etai = 1.e10_r8
  end if

  a1 = (dni * (te + ti) + zni * (dte + dti)) * 1.6022e-19_r8
  
  a2e = etai / (1._r8 + etai) * z / (z + 1._r8) * a1
  a2i = a2e

  hl31 = f * xx * (0.754_r8 + 2.210_r8 * z + z * z + xx &
   * (0.348_r8 + 1.243_r8 * z + z * z)) / dx
       
  hl32 = -f * xx * (0.884_r8 + 2.074_r8 * z) / dx
  
  hjbt = -hl31 * (a1 + alfai / z * a2i) - hl32 * a2e
  
  return
end subroutine neo_helena