write_eqdsk.f90

Go to the documentation of this file.

subroutine write_eqdsk(a, curr, r0, b0, qprof, path)
!-----------------------------------------------------------------------
! write eqdsk file, but without the flux surfaces, which are set to 0.
!-----------------------------------------------------------------------

  use itm_types
  use mod_parameter
  use mod_dat
  use mod_mesh
  use mod_helena_io
  use mod_profiles
  use mod_spline
  use helena_spline

  implicit none

  character(len = 132) :: path 
  integer(itm_i4), parameter :: out_eqdsk = 11

  real(r8), intent(in) :: a, curr, r0, b0
  real(r8), dimension(nr), intent(in) :: qprof

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

  real(r8), dimension(nr) :: p0tmp, dptmp, ftmp, dftmp
  real(r8), dimension(3) :: ablt
  real(r8) :: flux
  real(r8) :: radius, cpsurf
  real(r8) :: rbscale
  real(r8) :: reast, rwest, rgridl
  real(r8) :: zmid, rmaxis, zmaxis, xip
  integer(itm_i4) :: i, j, np2

  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)

  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)

  do i = 1, nr
    flux = (i - 1) / dble(nr - 1)
    if (hbt) then
      p0tmp(i) = 0.5_r8 * a * b * spwert(npts, flux, p_spline, psivec, ablt, &
       0) * sign(1.0_r8, -alfa)
      dptmp(i) = -0.5_r8 * a * b * dp_dpsi(flux)
      ftmp(i) = p0tmp(i) + eps * a * spwert(npts, flux, rbphi_spline, psivec, &
       ablt, 0) * sign(1.0_r8, -alfa) * sign(1.0_r8, a)
      dftmp(i) = dptmp(i) - eps * a * dgamma_dpsi(flux)
    else
      p0tmp(i) = eps * a * b * spwert(npts, flux, p_spline, psivec, ablt, 0) &
       * sign(1.0_r8, -alfa)
      dptmp(i) = -eps * a * b * dp_dpsi(flux) 
      ftmp(i) = eps * a * spwert(npts, flux, rbphi_spline, psivec, ablt, 0) &
       * sign(1.0_r8, -alfa) * sign(1.0_r8, a)
      dftmp(i) = -eps * a * dgamma_dpsi(flux)
    endif
    ftmp(i) = sqrt(1._r8 + 2._r8 * eps * ftmp(i) / alfa**2)
    dftmp(i) = 1._r8 / (2._r8 * ftmp(i)) * 2._r8 * eps * dftmp(i) / alfa**2
  end do

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

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

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

  np2 = np
  if (ias == 1) np2 = (np - 1) / 2 + 1
  reast = xx(1, 1) * eps * r0
  rwest = xx(1, np2) * eps * r0
  rgridl = reast - rwest

  zmid = yy(1, nr * np) * eps * r0
  rmaxis = (1._r8 + xx(1, nr * np) * eps) * r0
  zmaxis = zmid
  xip = curr * eps * r0 * b0 / zmu0

  open (unit = out_eqdsk, file = trim(adjustl(path)) // file_eqdsk_out, &
   status = 'replace', form = 'formatted', action = 'write', iostat = i_error)

  write(out_eqdsk, *) ' helena produced eqdsk file ', nr, np
  write(out_eqdsk, 11) 0., 0., r0, rgridl, zmid
  write(out_eqdsk, 11) rmaxis, zmaxis, 0., 0., b0
  write(out_eqdsk, 11) xip, 0., 0., rmaxis, 0.
  write(out_eqdsk, 11) zmaxis, 0., 0., 0., 0.

  write(out_eqdsk, 11) (ftmp(i), i = 1, nr)
  write(out_eqdsk, 11) (p0tmp(i), i = 1, nr)
  write(out_eqdsk, 11) (dftmp(i), i = 1, nr)
  write(out_eqdsk, 11) (dptmp(i), i = 1, nr)
  write(out_eqdsk, 11) (0., i = 1, nr * np)
  write(out_eqdsk, 11) (qprof(i), i = 1, nr)

  if (ias == 1) then
    write(out_eqdsk, *) np - 1, 1
    write(out_eqdsk, 11) (r0 * (1._r8 + eps * xx(1, j)), eps * r0 &
     * yy(1, j), j = 1, np - 1)
  else
    write(out_eqdsk, *) 2 * (np - 1), 1
    write(out_eqdsk, 11) (r0 * (1._r8 + eps * xx(1, j)), eps * r0 &
     * yy(1, j), j = 1, np)
    write(out_eqdsk, 11) (r0 * (1._r8 + eps * xx(1, np - j)), -eps * r0 &
     * yy(1, np - j), j = 1, np - 2)
  end if
  close(out_eqdsk)
!---------------------------------------------- end of eqdsk writing
  
 11 format(5e17.9)

  call deallocate_spline_coefficients(p_spline)
  call deallocate_spline_coefficients(rbphi_spline)

  return
end subroutine write_eqdsk