helena.f90

Go to the documentation of this file.

subroutine helena(equilibrium_in, equilibrium_out, in_path, code_parameters)
!-----------------------------------------------------------------------
!
! main program helena :
! ---------------------
!
!      - solves the 2d grad-shafranov equation for arbitrary
!        continuous plasma boundaries and equilibrium pressure
!        and gamma profiles.
!      - 2d cubic isoparametric finite elements are used for an accurate
!        representation of the solution.
!      - the final solution is obtained on a flux surface grid.
!
!
! program organization :
! ----------------------
!
!        blockdata                       : initialize namelist var.
!        helena                          : main program
!          initialization                : initialize values
!          write_parameters              : print output unit 20
!          gaussian_points               : initialize gaussian points
!            cubic_polynomials           : def. cubic elements
!          shape_soloviev                : boundary shape soloviev eq.
!            (zero)                      : hgolib routine
!            fr_soloviev
!              (rft2)                    : hgolib routine
!          shape_analytic
!            (gridinv)                   : hgolib routine
!            (rft2)
!          node_numbers                  : initialize element numbering
!          initial_grid                  : def. initial grid
!            radb                        : radial profile of elements
!          initkq                        : init. matrix kk and vector q
!          build_matrix                  : calc. matix kk and vector q
!            dp_dpsi                     : derivative of pressure
!            dgamma_dpsi                 : derivative of gamma
!            interpolation               : cubic elements interpolation
!          bounkq                        : insert boundary cond.
!            interpolation
!            delrc                       : delete rows and columns
!          solve                         : solve matrix equations
!            conjgr                      : conjugate gradients
!              shrink                    : remove rows/columns b.c.
!              scale                     : scale matrix problem
!              asub                      : matrix vector inproduct
!                expand                  : insert rows/columns b.c.
!                shrink
!           restore                      : insert b.c. in solution
!           find_axis                    : find magnetic axis
!             root                       : solve quadratic equation
!             cubic_interp_1d
!           normalize_psi                : normalize psi solution
!           remesh                       : calc. new isoparametric mesh
!             radial_mesh
!               rpack
!               drpack
!               ddrpack
!             thtmima                    : find minimum theta in elm.
!             position                   : find psi/theta point in elm.
!               (c05zaf)                 : nag routine
!               fzero
!                 interpolation
!               (c05pbf)                 : nag routine
!             posbnd                     : find psi/theta point boundary
!               solve_cubic              : solve cubic equation
!             interpolation
!             radial_mesh
!           mapping                      : calc. metric flux coord.
!             profiles                   : eq. profiles
!               pressure                 : pressure
!               gamma                    : gamma
!             radial_mesh
!             interpolation
!             solve_cubic
!           diagnostics                  : calc. beta, betapl
!
!
! the numbering used for the interpolating functions h :
!
!                r0,s0
!         h    h
!         h    h        with i,j = 0 or 1 indicates a derivative
!         hhhhhh        with respect to r,s. r0 and s0 refer to the
!         h    h        4 nodes of one element and can be +1. or -1.
!         h    h
!                i,j
!                                s
!                                .
!                             +1 .
!                    n2----------.-----------n3
!                     |          .           |
!                     |          .           |
!                     |          .           |        --> r
!           .......-1.|......................|.+1.......
!                     |          .           |
!                     |          .           |
!                     |          .           |
!                     |          .           |
!                    n1----------.-----------n4
!                                . -1
!                                .
!-----------------------------------------------------------------------

  use itm_types
  use mod_parameter

  use mod_corners
  use mod_dat
  use mod_dete
  use mod_helena_io
  use mod_map
  use mod_mesh
  use mod_meshacc
  use mod_nodes
  use mod_output
  use mod_profiles
  use mod_solver

  use mod_equilibrium
  use current_calculations
  use global_plasma_parameters
  use mod_interpolation, only : transfer_1d_profile
  use math_functions, only : equidistant_profile, integrate_1d, integrate
  use matrix_calculations, only : set_node_number
  use phys_profiles
  use plot_data
  use mapping_mod
  use solver, only : allocate_solver, deallocate_solver
  use helena_spline

  use mod_itm
  use euitm_schemas
  use euitm_xml_parser

  implicit none

!-- call parameters
  type (type_equilibrium), pointer :: equilibrium_in(:)
  type (type_equilibrium), pointer :: equilibrium_out(:)
  character(len = 132), optional :: in_path
  type (type_param) :: code_parameters

  type (type_equilibrium) :: equilibrium

  type(spline_coefficients) :: p_spline

  real(r8), parameter :: rel_change_w_mhd = 1.e-2_r8
  integer(itm_i4), parameter :: max_points_boundary = 1025

  real(r8), allocatable :: zvol(:), zvolp(:)
  real(r8), allocatable :: drmerc(:), dimerc(:), hh(:)
  real(r8), allocatable :: qprof(:), dqprof(:), geonc(:), zjpar(:)
  real(r8), allocatable :: fcirc(:), b02av(:), b0max(:), rav(:)
  real(r8), allocatable :: x_work(:), y_work(:), work(:), theta(:)
  real(r8), allocatable :: fr(:)
  real(r8), dimension(:), pointer :: frmajor, fz
  real(r8), allocatable :: dummy1(:)
  real(r8), allocatable :: pskn1(:)    ! Psi in increasing order
  real(r8), allocatable :: xx_center(:), yy_center(:)

  real(r8), dimension(3) :: abltg

  real(r8) :: a, a_input
  real(r8) :: rminor
  real(r8) :: rgeo, amin
  real(r8) :: b0
  real(r8) :: fscale
  real(r8) :: sax, rax
  real(r8) :: xaxis, yaxis, psaxis
  real(r8) :: cx, cy
  real(r8) :: sumdq
  real(r8) :: q95out, q1out
  real(r8) :: current, mhd_energy, volume
  real(r8) :: w_mhd_ref = 0._r8, b_ref
  real(r8) :: pressure
  integer(itm_i4) :: i, j, m
  integer(itm_i4) :: nax
  integer(itm_i4) :: n_bnd
  integer(itm_i4) :: meshno = 0
  integer(itm_i4) :: nelm, no, n1, n2, n3, n4
  integer(itm_i4) :: end_of_path

  logical :: converged_current = .false.
  logical :: npts_set = .false.
  logical :: first

  integer(itm_i4) :: return_status

  character(len = 132) :: codename(1) = 'HELENA'
  character(len = 132) :: codeversion(1) = '325'

  character(len = 2) :: file_no

  if (.not. present(in_path)) then
    path = './'  ! run locally
  else
    path = in_path
!-- determine end of path, defined by zero byte, required because of Kepler
    end_of_path = index(path, achar(0)) - 1
    if (end_of_path <= 0) then
      end_of_path = len_trim(adjustl(path))
      if (end_of_path == 0) then
        path = ''
      else
        path = trim(adjustl(path))
      end if
    else
      path = trim(adjustl(path(1 : end_of_path)))
    end if
  end if

  write(*, *) 'pre: helena_initialization'
  !call system('date +%s%N')

!-- initialize parameters
  call helena_initialization
  cpsurfin = 1._r8
  rminor = 1._r8
  zgeo = 0._r8

  write(*, *) 'pre: gaussian_points'
  !call system('date +%s%N')

!-- initialize interpolating functions
  call gaussian_points

  write(*, *) 'post: gaussian_points'
  !call system('date +%s%N')

!-- plasma boundary
  if (associated(equilibrium_in(1)%eqgeometry%boundary(1)%r) &
   .and. associated(equilibrium_in(1)%eqgeometry%boundary(1)%r)) then
    frmajor => equilibrium_in(1)%eqgeometry%boundary(1)%r
    fz => equilibrium_in(1)%eqgeometry%boundary(1)%z
!-- find last point in eqgeometry%boundary (may have filling zeros at the end)
    do n_bnd = size(frmajor), 1, -1
      if (frmajor(n_bnd) /= 0._r8 .or. fz(n_bnd) /= 0._r8) exit
    end do
!-- remove duplicated point(s) of closed boundary if necessary
    do
      if (n_bnd == 1) then
        write(iu6, *) 'ERROR: boundary has only 1 point'
        exit
      end if
      if (frmajor(1) == frmajor(n_bnd) .and. fz(1) == fz(n_bnd)) then
        n_bnd = n_bnd - 1
        write(iu6, *) 'removed duplicated point from boundary'
      else
        exit
      end if
    end do
    write(iu6, *) 'n_bnd = ', n_bnd
    rvac = (maxval(frmajor(1 : n_bnd)) + minval(frmajor(1 : n_bnd))) &
     / 2._r8
    rminor = (maxval(frmajor(1 : n_bnd)) - minval(frmajor(1 : n_bnd))) &
     / 2._r8
    write(*,*) 'Guessed rvac, rminor = ', rvac, rminor
  else
    write(iu6, *) 'ERROR: no plasma boundary specified'
    stop
  end if

!TODO: The following is a quick-fix for the CEA ETS workflow which needs
!      to be checked again
  if (associated(equilibrium_in(1)%profiles_1d%psi)) then
    cpsurfin = equilibrium_in(1)%profiles_1d%psi( &
     size(equilibrium_in(1)%profiles_1d%psi)) &
     - equilibrium_in(1)%profiles_1d%psi(1)
  end if

!-- save psi on axis or initialize it to zero
  if (exists_entry(equilibrium_in(1)%global_param%psi_ax)) then
    psi_axis = equilibrium_in(1)%global_param%psi_ax
  else
    psi_axis = 0._r8
  end if

!-- fields which may be overruled by helena.xml
  if (exists_entry(equilibrium_in(1)%global_param%psi_bound) &
   .and. exists_entry(equilibrium_in(1)%global_param%psi_ax)) then
    cpsurfin = (equilibrium_in(1)%global_param%psi_bound &
     - equilibrium_in(1)%global_param%psi_ax)
  else
    if (associated(equilibrium_in(1)%profiles_1d%psi)) then
      equilibrium_in(1)%global_param%psi_ax &
       = equilibrium_in(1)%profiles_1d%psi(1)
      equilibrium_in(1)%global_param%psi_bound &
       = equilibrium_in(1)%profiles_1d%psi(size( &
       equilibrium_in(1)%profiles_1d%psi))
    end if
  end if
  if (exists_entry(equilibrium_in(1)%eqgeometry%geom_axis%r)) &
   rvac = equilibrium_in(1)%eqgeometry%geom_axis%r
  if (exists_entry(equilibrium_in(1)%eqgeometry%a_minor)) &
   rminor = equilibrium_in(1)%eqgeometry%a_minor
  eps = rminor / rvac
  if (exists_entry(equilibrium_in(1)%global_param%mag_axis%bphi)) &
   bvac = equilibrium_in(1)%global_param%mag_axis%bphi
  if (exists_entry(equilibrium_in(1)%global_param%toroid_field%r0) &
   .and. exists_entry(equilibrium_in(1)%global_param%toroid_field%b0)) &
   bvac = equilibrium_in(1)%global_param%toroid_field%r0 &
    * equilibrium_in(1)%global_param%toroid_field%b0 / rvac
  if (exists_entry(equilibrium_in(1)%global_param%i_plasma)) &
   ip = equilibrium_in(1)%global_param%i_plasma

  allocate(equilibrium_out(1))
  allocate(equilibrium_out(1)%codeparam%codename(1))
  allocate(equilibrium_out(1)%codeparam%codeversion(1))
  if (.not. associated(code_parameters%parameters)) then
    write(iu6, *) 'ERROR: code parameters not associated!'
    stop
  else
    allocate(equilibrium_out(1)%codeparam%parameters(size( &
     code_parameters%parameters)))
  end if

!-- add to equilibrium_out
  equilibrium_out(1)%datainfo%cocos = 13
  equilibrium_out(1)%codeparam%codename = codename
  equilibrium_out(1)%codeparam%codeversion = codeversion
  equilibrium_out(1)%codeparam%parameters = code_parameters%parameters

  write(*, *) 'pre: helena_assign_code_parameters'
  !call system('date +%s%N')

!-- assign code parameters to internal variables
  call helena_assign_code_parameters(code_parameters, return_status)

  if (return_status /= 0) then
    write(iu6, *) 'ERROR: Could not assign code parameters.'
    return
  end if

!-- define r0 if not defined
  if (exists_entry(equilibrium_in(1)%global_param%toroid_field%r0)) then
    r_vessel = equilibrium_in(1)%global_param%toroid_field%r0
  else
    r_vessel = rvac
  end if

  if (trim(match) == "Ip" .and. ip == 0._r8) then
    write(iu6, *) "Error: for match == 'Ip' the total current must be" &
     // "specified"
    stop
  end if

  if (verbosity > 0) then
    write(iu6, *) '******************************************************'
    write(iu6, *) '*                     HELENA                         *'
    write(iu6, *) '******************************************************'
  end if

  if (verbosity > 0) write(iu6, *) 'done assigning code parameters'

!-- create run directory
  if (output /= 'none' .and. (standard_output .or. elite_output &
   .or. profiles_output .or. additional_output .or. xmgrace_output)) then
    call system('mkdir -p ' // trim(adjustl(path)))
    write(iu6, *) 'created run directory'
  end if

!-- create plot directory
  if (xmgrace_output) then
    call system('mkdir -p ' // trim(adjustl(path)) // 'plots')
    write(iu6, *) 'created plots directory'
  end if

!-- consistency checks 

  if (ias == 1 .and. mod(nchi, 2) /= 0) then
    write(iu6, *) 'Error: nchi must be 2^n'
    stop
  end if
  if (ias == 0 .and. mod(nchi, 2) /= 1) then 
    write(iu6, *) 'Error: nchi must be 2^n + 1'
    stop
  end if

  if (n_acc_points < 0 .or. n_acc_points > 10) then
    write(iu6, *) 'Error: condition 0 <= n_acc_points <= 10 not fulfilled!'
    stop
  end if
  if (equidistant < 0 .or. equidistant > 1) then
    write(iu6, *) 'Error: condition 0 <= equidistant <= 1 not fulfilled!'
    stop
  end if

  if (te%shape > 2) then
    write(iu6, *) 'Error: 0 <= te%shape <= 2 !'
    stop
  end if
  if (ne%shape > 2) then
    write(iu6, *) 'Error: 0 <= ne%shape <= 2 !'
    stop
  end if

  if (wmhd < 0._r8) then
    write(iu6, *) 'Error: wmhd must be positive !'
    stop
  end if

  if (verbosity > 1) then
    write(iu6, *) 'cpsurfin ', cpsurfin
    write(iu6, *) 'rvac ', rvac
    write(iu6, *) 'rminor ', rminor
    write(iu6, *) 'eps ', eps
    write(iu6, *) 'bvac ', bvac
  end if

  if (verbosity > 2) then
    write(iu6, *) 'W_MHD = ', wmhd
    write(iu6, *) 'mfm   = ', mfm
  end if

  i = 2
  do
    if (mfm < i) then
      write(iu6, *) 'ERROR: mfm not a power of 2'
      return
    end if
    if (mfm == i) exit
    i = i * 2
  end do

  write(*, *) 'pre: shape'
  !call system('date +%s%N')

!-- mfm is the number of points specifying the boundary and at the same
!   time the number of Fourier components
!   i.e. mharm = mfm / 2 is the number of Fourier harmonics
  mharm = mfm / 2

  allocate(fr(mfm + 2))
  fr = 0._r8
!-- calculate radii from boundary
  if (verbosity > 1) write(iu6, *) 'shape from discrete points'
  call shape_from_points(frmajor(1: n_bnd), fz(1 : n_bnd), n_bnd, mfm, rgeo, &
   amin, fr)
  fr = fr / amin  ! normalize to minor radius of input boundary
!-- Rgeo is discarded and replaced by rvac, while Zgeo is kept since
!   vertical shifts do not affect the equilibrium

  allocate(theta(mfm))

  do m = 1, mfm
    theta(m) = 2._r8 * pi * (m - 1) / dble(mfm - 1)
  end do

  if (xmgrace_output) then
    call plot_data_1d('line', theta, fr, mfm, 'xmgr_boundary_orig')
  end if

  deallocate(theta)

!-- Fourier transform of fr(1:mfm) yields mfm Fourier coefficients fr(1:mfm)
!   rft2 requires dimension mfm + 2 for fr !
  call rft2(fr, mfm, 1)

  do m = 1, mfm
    fr(m) = 2._r8 * fr(m) / dble(mfm)
  end do
  do m = 2, mfm, 2
    fr(m) = -fr(m)    ! invert sign of imaginary coefficients
  end do

  call fourier_backtransform(fr, mfm)

  write(*, *) 'post: shape'
  !call system('date +%s%N')

!-- calculate constants
  rminor = eps * rvac
  if (verbosity > 2) write(iu6, *) 'minor radius a = ', rminor

  if (standard_output) &
   open (unit = out_he, file = trim(adjustl(path)) // file_he_out, &
   status = 'replace', form = 'formatted', action = 'write', iostat = i_error)

!-- soloviev settings
  if (isol == 1) then
    a = 2._r8 * (1._r8 + (1._r8 - eps**2 / 4._r8) / ellip**2)
    b = 2._r8 * eps + tria / ellip**2 / (1._r8 + (1._r8 - eps**2 &
     / 4._r8) / ellip**2)
    cpsurf = 0.5_r8 * eps**2 / (1._r8 + eps**2)**2 * ellip &
     / sqrt(1._r8 - eps**2 / 4._r8)
    radius = sqrt(eps**2 / (1._r8 + eps**2))
    b0 = sqrt(1._r8 + eps**2)
    alfa = radius**2 * b0 / cpsurf
    niter = 1
  end if

  txtout(1) = ' '

!-- extract fixed fields from equilibrium_in

!-- 1D profiles
  if (associated(equilibrium_in(1)%profiles_1d%pprime)) then
    if (verbosity > 1) write(iu6, *) 'pprime specified on input'
    npts_set = .true.
    npts = size(equilibrium_in(1)%profiles_1d%pprime)
    allocate(dp(npts))
    dp = equilibrium_in(1)%profiles_1d%pprime
  end if 
  if (associated(equilibrium_in(1)%profiles_1d%pressure)) then
    if (npts_set .and. size(equilibrium_in(1)%profiles_1d%pressure) &
     /= npts) then
      write(iu6, *) 'ERROR: input profiles have inconsistent sizes'
      stop
    else
      if (verbosity > 1) write(iu6, *) 'pressure specified on input'
      npts_set = .true.
      npts = size(equilibrium_in(1)%profiles_1d%pressure)
      p_sep = equilibrium_in(1)%profiles_1d%pressure(npts)
    end if
  else
    p_sep = 0._r8
  end if
  if (associated(equilibrium_in(1)%profiles_1d%ffprime)) then
    if (npts_set .and. size(equilibrium_in(1)%profiles_1d%ffprime) &
     /= npts) then
      write(iu6, *) 'ERROR: input profiles have inconsistent sizes'
      stop
    else
      if (verbosity > 1) write(iu6, *) 'ffprime specified on input'
      npts_set = .true.
      npts = size(equilibrium_in(1)%profiles_1d%ffprime)
      allocate(fdf(npts))
      fdf = equilibrium_in(1)%profiles_1d%ffprime
    end if
  end if
  if (associated(equilibrium_in(1)%profiles_1d%jphi)) then  
    if (npts_set .and. size(equilibrium_in(1)%profiles_1d%jphi) &
     /= npts) then
      write(iu6, *) 'ERROR: input profiles have inconsistent sizes'
      stop
    else
      if (verbosity > 1) write(iu6, *) 'j_tor specified on input'
      npts_set = .true.
      npts = size(equilibrium_in(1)%profiles_1d%jphi)
      allocate(j_tor(npts))
      j_tor = equilibrium_in(1)%profiles_1d%jphi
    end if
  end if
  if (associated(equilibrium_in(1)%profiles_1d%q)) then
    if (npts_set .and. size(equilibrium_in(1)%profiles_1d%q) &
     /= npts) then
      write(iu6, *) 'ERROR: input profiles have inconsistent sizes'
      stop
    else
      if (verbosity > 1) write(iu6, *) 'q specified on input'
      npts_set = .true.
      npts = size(equilibrium_in(1)%profiles_1d%q)
      allocate(q_in(npts))
      q_in = equilibrium_in(1)%profiles_1d%q
    end if
  end if

  if (.not. npts_set) then
    write(iu6, *) 'fatal error: no input profiles specified'
    stop
  end if

  if (associated(equilibrium_in(1)%profiles_1d%psi)) then
    if (npts_set .and. size(equilibrium_in(1)%profiles_1d%psi) &
     /= npts) then
      write(iu6, *) 'ERROR: input profiles have inconsistent sizes'
      stop
    else
      if (verbosity > 1) write(iu6, *) 'psi specified on input'
      npts_set = .true.
      npts = size(equilibrium_in(1)%profiles_1d%psi)
      allocate(psi_in(npts))
      if (equilibrium_in(1)%profiles_1d%psi(1) /= 0._r8) then
        if (verbosity > 2) write(iu6, *) 'psi shifted to be 0 on axis'
      end if
      psi_in = (equilibrium_in(1)%profiles_1d%psi &
       - equilibrium_in(1)%profiles_1d%psi(1))
    end if
  else
!TODO: improve this initial assumption for cases with rho_tor or
!      rho_vol specifications
    if (verbosity >1) then
      write(iu6, *) 'no psi specified on input'
      write(iu6, *) 'assuming equidistant psi'
    end if
    allocate(psi_in(npts))
    do i = 1, npts
      psi_in(i) = dble(i - 1) / dble(npts - 1) * cpsurfin
    end do
  end if

  if (.not. associated(equilibrium_in(1)%profiles_1d%psi)) &
   allocate(equilibrium_in(1)%profiles_1d%psi(npts))

!-- re-set input psi
  equilibrium_in(1)%profiles_1d%psi = psi_in

!-- interpolate input p on input psi to determine derivative dp/dpsi if
!   input_type == "p and j_tor"
  if (trim(input_type) == "p and j_tor") then
    call allocate_spline_coefficients(p_spline, npts)
    call spline(npts, psi_in, equilibrium_in(1)%profiles_1d%pressure, &
     0._r8, 0._r8, 2, p_spline)
    if (.not. associated(equilibrium_in(1)%profiles_1d%pprime)) &
     allocate(equilibrium_in(1)%profiles_1d%pprime(npts))
    do i = 1, npts
      pressure = spwert(npts, psi_in(i), p_spline, psi_in, abltg, 1)
      equilibrium_in(1)%profiles_1d%pprime(i) = abltg(1)
    end do
    call deallocate_spline_coefficients(p_spline)
    if (verbosity > 2) &
     write(iu6, *) "interpolated p to derive p'"
  end if

!-- xmgr output of input fields
  if (xmgrace_output) then
    if (associated(equilibrium_in(1)%profiles_1d%pprime)) &
     call plot_data_1d('line', psi_in, equilibrium_in(1)%profiles_1d%pprime, npts, &
     'xmgr_input_pprime')
    if (associated(equilibrium_in(1)%profiles_1d%pressure)) &
     call plot_data_1d('line', psi_in, equilibrium_in(1)%profiles_1d%pressure, npts, &
     'xmgr_input_pressure')
    if (associated(equilibrium_in(1)%profiles_1d%ffprime)) &
     call plot_data_1d('line', psi_in, equilibrium_in(1)%profiles_1d%ffprime, npts, &
     'xmgr_input_ffprime')
    if (associated(equilibrium_in(1)%profiles_1d%jphi)) &
     call plot_data_1d('line', psi_in, equilibrium_in(1)%profiles_1d%jphi, npts, &
     'xmgr_input_jphi')
    if (associated(equilibrium_in(1)%profiles_1d%q)) &
     call plot_data_1d('line', psi_in, equilibrium_in(1)%profiles_1d%q, npts, &
     'xmgr_input_q')
  end if

  write(*, *) 'pre: profiles'
  !call system('date +%s%N')

!-- allocate missing fields
  if (.not. allocated(fdf)) allocate(fdf(npts))

!-- allocate arrays
  allocate(zvol(nrmap))
  allocate(zvolp(nrmap))

  allocate(psivec(npts))
  allocate(dpres(npts))
  allocate(p_int(npts))
  allocate(dgam(npts))
  allocate(gam_int(npts))

  call allocate_spline_coefficients(dp_spline, npts)
  call allocate_spline_coefficients(gam_spline, npts)
  call allocate_spline_coefficients(j_spline, npts)

!-- define psivec as normalized psi_in
  if (psi_in(1) /= 0._r8) then
    write(iu6, *) 'fatal error: psi on axis must be zero'
    stop
  end if
  psivec = psi_in / psi_in(npts)

!-- set profiles dpres and dgam
  dpres = equilibrium_in(1)%profiles_1d%pprime
  if (trim(input_type) == "p' and FF'") then
    dgam = equilibrium_in(1)%profiles_1d%ffprime
  end if

!-- spline profiles
  call spline(npts, psivec, dpres, 0._r8, 0._r8, 2, dp_spline)
  if (trim(input_type) == "p' and FF'") then
    call spline(npts, psivec, dgam, 0._r8, 0._r8, 2, gam_spline)
  end if
  if (allocated(j_tor)) &
   call spline(npts, psivec, j_tor, 0._r8, 0._r8, 2, j_spline)

!-- integrate profiles
  p_int = integrate(dp_spline, psivec, npts, .true.) &
   * sign(1.0_r8, cpsurfin)
  if (trim(input_type) == "p' and FF'") then
    gam_int = integrate(gam_spline, psivec, npts, .true.) &
     * sign(1.0_r8, cpsurfin)
  end if

!-- initialize a based on gamma on axis
  if (trim(input_type) == "p' and FF'") then
    a = -dgam(1) * (twopi * eps * rvac)**2 / cpsurfin
  else
    a = 4
  end if
  if (verbosity > 2) write(iu6, *) 'initial a = ', a
  a_input = a

!-- initialize alfa if it hasn't been specified on input
  if (alfa == 0._r8) &
   alfa = twopi * (eps * rvac)**2 * bvac / cpsurfin
  if (verbosity > 2) write(iu6, *) 'initial alfa = ', alfa

!-- initial pressure profile, ffprime, and f^2 profiles
  if (xmgrace_output) then
    call plot_data_1d('line', psivec, dpres, npts, 'xmgr_initial_dp_psipol')
    call plot_data_1d('line', psivec, p_int, npts, &
     'xmgr_initial_p_psipol')
    if (trim(input_type) == "p' and FF'") then
      allocate(work(npts))
      call plot_data_1d('line', psivec, dgam, npts, &
       'xmgr_initial_ffprime_psipol')
      work = 2._r8 * gam_int
      call plot_data_1d('line', psivec, work, npts, 'xmgr_initial_f2_psipol')
!-- attention: since gam_int has wrong magnitude due to normalized psivec
!              (rvac * bvac)**2 + work may be negative => use abs
      work = sqrt(abs((rvac * bvac)**2 + work)) * sign(1.0_r8, bvac)
      call plot_data_1d('line', psivec, work, npts, &
       'xmgr_initial_rbphi_psipol')
      deallocate(work)
    end if
  end if

  write(*, *) 'pre: equilibrium constructor'
  !call system('date +%s%N')

!-- construct equilibrium
  call equilibrium_constructor(equilibrium, nr, nchi)

!-- allocate large matrix
  kklda = 4 * np + 9

  allocate(kkbig(kklda, 4 * (nr + 2) * np))

  if (trim(input_type) == "p' and j_tor" &
   .or. trim(input_type) == "p and j_tor" &
   .or. trim(input_type) == "p' and q") then
!-- no equidistant theta points for initial grid with current iteration
    iarc = 0
  end if

!-- normalize all profiles

  if (trim(input_type) == "p' and FF'") then
    fscale = dgam(1)
    dgam = dgam / fscale
    gam_int = gam_int / fscale
    gam_spline%sp1 = gam_spline%sp1 / fscale
    gam_spline%sp2 = gam_spline%sp2 / fscale
    gam_spline%sp3 = gam_spline%sp3 / fscale
    gam_spline%sp4 = gam_spline%sp4 / fscale
  end if

  pscale = dpres(1)
  dpres = dpres / pscale
  p_int = p_int / pscale
  dp_spline%sp1 = dp_spline%sp1 / pscale
  dp_spline%sp2 = dp_spline%sp2 / pscale
  dp_spline%sp3 = dp_spline%sp3 / pscale
  dp_spline%sp4 = dp_spline%sp4 / pscale
 
!-- normalize j_tor to unity on axis (required!)
  if (trim(input_type) == "p' and j_tor" &
   .or. trim(input_type) == "p and j_tor" &
   .or. trim(input_type) == "p' and q") then
    cscale = j_tor(1)
    j_tor = j_tor / cscale
    j_spline%sp1 = j_spline%sp1 / cscale
    j_spline%sp2 = j_spline%sp2 / cscale
    j_spline%sp3 = j_spline%sp3 / cscale
    j_spline%sp4 = j_spline%sp4 / cscale
  end if

!-- if no B specified on input calculate B from pscale and fscale

  if (b == 0._r8 .and. trim(input_type) == "p' and FF'") then
    if (verbosity > 1) then
      write(iu6, *) 'no B specified on input'
      write(iu6, *) 'automatic calculation of B'
    end if
    if (hbt) then
      b = 2._r8 * eps * mu0 * rvac**2 * pscale / fscale
    else
      b = mu0 * rvac**2 * pscale / fscale
    end if
  end if

  if (trim(input_type) == "p' and FF'") then
    if (verbosity > 2) then
      write(iu6, *) "================================================="
      write(iu6, *) "consistency check for input parameter B: "
      write(iu6, *) "input p'(psi) = dp/dpsi and"
      write(iu6, *) "F(psi)F'(psi) = 0.5 * dF^2/dpsi"
      write(iu6, *) "should be in SI units" 
      write(iu6, *) "p'(0) (on input) = ", pscale
      write(iu6, *) "F(0)F'(0) (on input) = ", fscale
      if (hbt) then
        write(iu6, *) "B = ", b, " ?= ", 2._r8 * eps * mu0 * rvac**2 &
         * pscale / fscale
      else
        write(iu6, *) "B = ", b, " ?= ", mu0 * rvac**2 * pscale / fscale
      end if
      write(iu6, *) "================================================="
    end if
  end if

  if (standard_output) call write_parameters
  if (isol == 1) then
    if (standard_output) then
      write(out_he, *)
      write(out_he, *) '               soloviev equilibrium '
      write(out_he, 2) a
      write(out_he, 3) b
      write(out_he, 4) alfa
      write(out_he, 5) radius
      write(out_he, 6) b0
      write(out_he, 7) cpsurf
    end if
    txtout(1) = 'soloviev equilibrium'
  end if
!TODO: re-organize verbose output
  if (standard_output) then
    write(out_he, *) '******************************************'
    write(out_he, *) '*        input profiles :                *'
    write(out_he, *) '******************************************'
    if (trim(input_type) == "p' and FF'" .and. .not. hbt) then
      write(out_he, *) &
       '*  psi,       dp/dpsi,    FdF/dpsi                   *'
    else if (trim(input_type) == "p' and FF'" .and. hbt) then
      write(out_he, *) &
       '*  psi,       dp/dpsi,     gamma                     *'
    else if (trim(input_type) == "p' and q" .and. hbt) then
      write(out_he, *) &
       '*  psi,       dp/dpsi,     gamma,         q          *'
    else if (trim(input_type) == "p' and q" .and. .not. hbt) then
      write(out_he, *) &
       '*  psi,       dp/dpsi,     FdF/dpsi,      q          *'
    else if (trim(input_type) == "p' and j_tor" .and. .not. hbt) then
      write(out_he, *) &
       '*  psi,       dp/dpsi,    FdF/dpsi,    <j_tor>       *'
    else if (trim(input_type) == "p' and j_tor" .and. hbt) then
      write(out_he, *) &
       '*  psi,       dp/dpsi,     gamma,      <j_tor>       *'
    else if (trim(input_type) == "p and j_tor" .and. .not. hbt) then
      write(out_he, *) &
       '*  psi,         p,        FdF/dpsi,    <j_tor>       *'
    else if (trim(input_type) == "p and j_tor" .and. hbt) then
      write(out_he, *) &
       '*  psi,         p,         gamma,      <j_tor>       *'
    end if      
    write(out_he, *) &
     '******************************************************'  
  end if

  write(*, *) 'pre: allocate'
  !call system('date +%s%N')

! allocate arrays to initial grid size
  allocate(psi(4 * nr * np))
  allocate(xx(4, nr * np))
  allocate(yy(4, nr * np))
  allocate(xxold(4, nr * np))
  allocate(yyold(4, nr * np))
  allocate(psiold(4 * nr * np))
  allocate(nodeno((nr - 1) * (np - 1), 4))
             
  psi = 0._r8
  xx = 0._r8
  yy = 0._r8
  psiold = 0._r8
  xxold = 0._r8
  yyold = 0._r8

  allocate(x_work(nr * np))
  allocate(y_work(nr * np))

!-- allocate solver variables
  call allocate_solver(nr, np)

  write(*, *) 'pre: initial grid'
  !call system('date +%s%N')

!-- define initial grid x,y
  allocate(sg(nr))
  allocate(dsg(nr))
  allocate(ddsg(nr))
  call initial_grid(fr)
  xaxis = 0._r8   ! initial axis

  if (xmgrace_output) then
    x_work = xx(1, 1 : nr * np) * rminor + rvac
    y_work = yy(1, 1 : nr * np) * rminor
    call plot_data_1d('2d', x_work, y_work, np, &
     'xmgr_initial_flux_surfaces_real')
    call plot_data_1d('2d_t', x_work, y_work, np, &
     'xmgr_initial_radial_real')
    call plot_data_1d('grid', x_work, y_work, np, &
     'xmgr_initial_grid_x_y_real')
  end if

  write(*, *) 'pre: node_numbers'
  !call system('date +%s%N')

!-- define numbering of elements
  call node_numbers

  write(*, *) 'pre: flxint'
  !call system('date +%s%N')

!-- initialize fdf if input_type /= "p' and FF'"
  if (trim(input_type) == "p' and j_tor" &
   .or. trim(input_type) == "p and j_tor" &
   .or. trim(input_type) == "p' and q") then
!-- initialize fdf based on current profile
    if (trim(input_type) == "p' and j_tor" &
     .or. trim(input_type) == "p and j_tor") then
      call flxint(xaxis, a, 0, fscale)
    else
!-- initialize fdf based on q-profile 
      call flxint2(xaxis, a, cx, cy, sumdq)
    end if

!-- re-calculate gam_spline and gam_int after change of fdf (already
!   normalized)
    dgam = fdf
!-- spline profiles
    call spline(npts, psivec, dgam, 0._r8, 0._r8, 2, gam_spline)
!-- integrate profiles
    gam_int = integrate(gam_spline, psivec, npts, .true.)
  end if

!-- allocate and set psikn for the initial grid
  allocate(psikn(nr))
  allocate(dpsikn(nr))
  allocate(ddpsikn(nr))
  allocate(radpsi(nr))
  do i = 1, nr
    psikn(nr - i + 1) = sg(i)**2
    dpsikn(nr - i + 1) = 2._r8 * sg(i) * dsg(i)
    ddpsikn(nr - i + 1) = 2._r8 * sg(i) * ddsg(i) + 2._r8 * dsg(i)**2
    radpsi(nr - i + 1) = dble(i - 1) / dble(nr - 1)
  end do

!-- fill equilibrium for zero order guesses of W_MHD and Ip
!   (will differ from the final values since no real equilibrium)
!  call fill_equilibrium(a, xaxis, yaxis, current_averaging, equilibrium)

!-- needed for update_fdf
  equilibrium%eqgeometry%geom_axis%r = rvac

!-- initial guess for psi for first iteration
  do i = 1, nr
    equilibrium%profiles_1d%psi(i) = cpsurfin * psikn(nr + 1 - i)
  end do
  equilibrium%profiles_1d%psi(1) = 0._r8

  write(*, *) 'pre: transfer'
  !call system('date +%s%N')

!-- obsolete ----
  if (associated(equilibrium_in(1)%profiles_1d%jphi)) then
    call transfer_1d_profile(psi_in, &
     equilibrium_in(1)%profiles_1d%jphi, &
     equilibrium%profiles_1d%psi, equilibrium%profiles_1d%jphi, &
     spline_type = 1)
  end if  
  if (associated(equilibrium_in(1)%profiles_1d%pprime)) then
    call transfer_1d_profile(psi_in, &
     equilibrium_in(1)%profiles_1d%pprime, &
     equilibrium%profiles_1d%psi, equilibrium%profiles_1d%pprime, &
     spline_type = 1)
  end if
  if (associated(equilibrium_in(1)%profiles_1d%pressure)) then
    call transfer_1d_profile(psi_in, &
     equilibrium_in(1)%profiles_1d%pressure, &
     equilibrium%profiles_1d%psi, equilibrium%profiles_1d%pressure, &
     spline_type = 1)
  end if
!----------------

!-- calculate W_MHD
!TODO: mhd_energy not correct since equilibrium%profiles_1d%pressure not
!      defined
  write(*, *) 'pre: w_mhd'
  !call system('date +%s%N')

  if (verbosity > 3) then
    mhd_energy = w_mhd(equilibrium, rminor, eps)
    volume = plasma_volume(equilibrium, rminor, eps)
    write(iu6, *) 'initial W_MHD = ', mhd_energy, ' J'
    write(iu6, *) 'initial plasma volume = ', volume, ' m^3'
  end if

  first = .true.

!-- allocated quantities for calls to profiles
  allocate(p0(nr))
  allocate(rbphi(nr))
  allocate(dp0(nr))
  allocate(drbphi(nr))

!---------------------------------------------------------------------------
!-- beginning of reference part
!---------------------------------------------------------------------------


!-- initialize old grid before call to p_j_iteration (alternatively call
!   remesh)
  xxold = xx
  yyold = yy
  psiold = psi

  write(*, *) 'pre: p_j'
  !call system('date +%s%N')

!-- loop to adapt <j_tor> and <dp/dpsi>
  call p_j_iteration(equilibrium_in(1), equilibrium, first, a, cx, cy, &
   psaxis, xaxis, yaxis, nax, rax, sax)

!TODO: in p_j_iteration: add rescaling of profiles_1d%pprime and
!      profiles_1d%pressure to match plasma beta (if given) or W_MHD
!      (if given)
!-- adapt plasma beta through change of B
!          call betapol(a, bpl)
!          if (betap >= 0._r8 .and. abs(bpl - betap) > 1e-3_r8) then
!            if (nmg == nbb) then
!              b1 = b
!              bp1 = bpl
!              if (b < 0._r8) then
!                if (bpl > betap) then
!                  b = 1.02_r8 * b1
!                else
!                  b = 0.98_r8 * b1
!                end if
!!                b =  bpl / (bbb * betap) * b1
!              else
!                b =   b1 / bpl * betap
!              end if
!              if (verbosity > 2) write(iu6, '(a18, 2e12.4, a4, e12.4)') &
!               'first change of b ', b1, bpl, ' >> ', b
!            else
!              if (mod(nmg, nbb) == 0) then
!                bold = b
!                deltab = (b1 - bold) / (bp1 - bpl) * (betap - bpl)
!                if (abs(deltab / bold) > bbb) then
!                  b = bold + bbb * deltab / abs(deltab) * abs(bold)
!                else
!                  b = bold + abb * (b1 - bold) / (bp1 - bpl) * (betap - bpl)
!                end if
!                if (verbosity > 2) write(iu6, '(a17, 4f12.4, a4, f12.4)') &
!                 'next change of b ', bold, bpl, b1, bp1, ' >> ', b
!                b1 = bold
!                bp1 = bpl
!              end if
!            end if
!            if (verbosity > 1) write(iu6, *) 'betapol'
!          end if


!-- re-install old mesh for further iterations (new mesh exists on exit from
!   p_j_iteration)
  xx = xxold
  yy = yyold
  psi = psiold

!TODO: implement this into p_j_iteration
!-- adapt B after each current iteration if wmhd and j_tor specified,
!   except for last iteration
!      if (verbosity > 1) then
!        write(iu6, *) 'w_mhd_ref = ', w_mhd_ref
!        write(iu6, *) 'mhd_energy = ', mhd_energy
!      end if
!      if (wmhd /= 0._r8 .and. abs(w_mhd_ref - mhd_energy) / w_mhd_ref &
!       > rel_change_W_MHD) then
!        converged_current = .false.
!        if (verbosity > 1) write(iu6, *) 'converged_current = .false.'
!        if (j /= nouter .and. nmg /= nmesh) then
!          B_ref = B
!          if (verbosity > 1) write(iu6, *) 'old B = ', B_ref
!          B = w_mhd_ref / mhd_energy * B
!          if (verbosity > 1) write(iu6, *) 'new B = ', B
!        else
!          if (verbosity > 1) &
!           write(iu6, *) 'warning: current iteration not converged'
!        end if
!      end if

!-- deallocate solver variables
  call deallocate_solver

  first = .true.

!-- dealllocate arrays
  deallocate(kkbig)
  deallocate(x_work, y_work)
  deallocate(psikn, dpsikn, ddpsikn, radpsi)

  if (verbosity > 2) then
    write(iu6, *) "================================================"
    write(iu6, *) "consistency check for constrained parameter A:  "
    write(iu6, *) "A on input = ", a_input
    write(iu6, *) "A after final iteration = ", a
    write(iu6, *) "================================================"
  end if

!-- remesh on flux coordinates
  if (verbosity > 1) write(iu6, *) ' final remesh'

  allocate(psikn(nrmap))
  allocate(dpsikn(nrmap))
  allocate(ddpsikn(nrmap))
  allocate(radpsi(nrmap))

  write(*, *) 'pre: remesh'
  !call system('date +%s%N')

  call remesh(a, nrmap, npmap, meshno, cx, cy, xaxis, yaxis, nax, rax, sax)
  if (verbosity > 1) write(iu6, *) ' done final remesh'

!-- fill equilibrium for calculation of W_MHD and Ip
  call equilibrium_destructor(equilibrium)
  write(*, *) 'pre: fill_equilibrium'
  !call system('date +%s%N')
  call equilibrium_constructor(equilibrium, nr, nchi)
  call fill_equilibrium(a, xaxis, yaxis, current_averaging, equilibrium)

  write(*, *) 'pre: w_mhd'
  !call system('date +%s%N')
!-- calculate W_MHD
  if (verbosity > 2) then
    mhd_energy = w_mhd(equilibrium, rminor, eps)
    volume = plasma_volume(equilibrium, rminor, eps)
    write(iu6, *) 'W_MHD = ', mhd_energy, ' J'
    write(iu6, *) 'plasma volume = ', volume, ' m^3'
  end if

  if (xmgrace_output) then
    call plot_data_1d('line', equilibrium%profiles_1d%volume, &
     equilibrium%profiles_1d%pressure, nr, 'xmgr_p_final')
    call plot_data_1d('line', equilibrium%profiles_1d%volume, &
     equilibrium%profiles_1d%pprime, nr, 'xmgr_pprime_final')
    call plot_data_1d('line', equilibrium%profiles_1d%volume, &
     equilibrium%profiles_1d%jphi, nr, 'xmgr_jphi_final')
  end if

  deallocate(xxold, yyold, psiold)
  deallocate(sg, dsg, ddsg)

  allocate(fcirc(nr))
  allocate(b02av(nr))
  allocate(b0max(nr))
  allocate(rav(nr))

!TODO: from former current scaling
!  call passing_particles(a, fcirc, b02av, b0max, rav)

!-- allocate j_tor if not allocated (if j_phi was analytically given)
  if (.not. allocated(j_tor)) allocate(j_tor(npts))

  write(*, *) 'pre: rescale alfa'
  !call system('date +%s%N')
!-- rescale alfa to input current or input q95
  call current_calculations_constructor(nr, np)
  if (trim(match) == "q95") then
    if (verbosity > 1) write(iu6, *) '  safety factor', q95
    allocate(dummy1(npts))
    call safety_factor(xaxis, a, cx, cy, dummy1, q95out, q1out)
    deallocate(dummy1)
    if (verbosity > 1) write(iu6, *) q1out, q95, alfa
    alfa = alfa * q95 / q1out
    call total_current(a, xaxis, current)
    ip = current
  else if (trim(match) == "Ip") then      
    if (verbosity > 1) write(iu6, *) ' current Ip = ', ip      
    call total_current(a, xaxis, current)
    if (verbosity > 1) write(iu6, *) ' calculated current = ', current
    alfa = current / ip
    if (verbosity > 1) write(iu6, *) '   alfa : ', alfa
  else    ! match = "psi_bound"
    call total_current(a, xaxis, current)
  end if

  write(*, *) 'pre: si_currents'
  !call system('date +%s%N')
!-- denormalize currents
  call si_currents

  allocate(pskn1(nr))
  do i = 1, nr
    pskn1(i) = psikn(nr + 1 - i)
  end do 

!-- Cartesian coordinates of the geometric centers of each element
  allocate(xx_center((nr - 1) * (np - 1)))
  allocate(yy_center((nr - 1) * (np - 1)))
!-- because of contouring issues in IDL, the magnetic axis is excluded from
!   the output
  no = 0
  do i = 2, nr
    do j = 1, np -1
      no = no + 1
      if (i == 1) then
        xx_center(no) = xaxis
        yy_center(no) = yaxis
      else
!-- elements in HELENA are ordered in inverse radial direction, i.e.,
!   psikn(nr) = 0, psikn(1) = 1
        nelm = (nr - i) * (np - 1) + j
        call set_node_number(nelm, n1, n2, n3, n4)
        xx_center(no) = (xx(1, n1) + xx(1, n2) + xx(1, n3) + xx(1, n4)) &
         / 4._r8
        yy_center(no) = (yy(1, n1) + yy(1, n2) + yy(1, n3) + yy(1, n4)) &
         / 4._r8
      end if
    end do
  end do
!-- output currents, averaged over area

  if (xmgrace_output) then
    call plot_data_1d('line', pskn1, j_tor_av, nr, 'xmgr_j_tor_av_psi')
    call plot_data_1d('line', pskn1, i_tor, nr, 'xmgr_I_tor_psi')
    call plot_data_1d('3d', xx_center, yy_center, nr - 1, 'idl_j_tor_loc', &
     j_tor_loc)
  end if

  deallocate(xx_center, yy_center)
  call current_calculations_destructor

  call deallocate_spline_coefficients(j_spline)
  deallocate(j_tor)

  if (xmgrace_output) then
     call plot_data_1d('line', abs(equilibrium%profiles_1d%psi), equilibrium%profiles_1d%ffprime, nr, &
         'xmgr_resulting_ffprime')
  end if

  call write_equidistant_profiles(equilibrium, pskn1, nr)

!-- j_spline needed in diagnostics only
  call allocate_spline_coefficients(j_spline, npts)
  allocate(j_tor(npts))
  call transfer_1d_profile(pskn1, equilibrium%profiles_1d%jphi, &
   psivec, j_tor, spline_type = 1)
  call spline(npts, psivec, j_tor, 0._r8, 0._r8, 2, j_spline)

!-- normalize j_tor
  cscale = j_tor(1)
  j_tor = j_tor / cscale
  j_spline%sp1 = j_spline%sp1 / cscale
  j_spline%sp2 = j_spline%sp2 / cscale
  j_spline%sp3 = j_spline%sp3 / cscale
  j_spline%sp4 = j_spline%sp4 / cscale

  write(*, *) 'pre: diagnostics'
  !call system('date +%s%N')
!-- diagnostics
  if (verbosity > 3) write(txtout(13), 13) xaxis, yaxis
  call diagnostics(a, xaxis, yaxis, zvol, zvolp, current)
!attention: diagnostics recalculates current !!!

  if (verbosity > 2) then
    write(iu6, *) &
     '-----------------------------------------------------'
    write(iu6, *) &
     'consistency check for Ip and computed alfa:'
    write(iu6, *) &
     'current = ',current,', unnormalized current (from Ip) = ', &
     ip / 1.e+3_r8, ' kA,  alfa = ', alfa
    write(iu6, *) &
     'current = ',current,', unnormalized current (calculated) = ', &
     current / 1.e+3_r8, ' kA, alfa = ',alfa  
    write(iu6, *) &
     '-----------------------------------------------------'
  end if

!TODO: from former current scaling
!  call scale_current(a, rvac, bvac, rav, fcirc, kkk)

  deallocate(fr)

!-- deallocate profiles
  deallocate(p0, rbphi, dp0, drbphi)

!-- mapping of flux coordinates
  allocate(g11_hel(nr, nchi))
  allocate(g12_hel(nr, nchi))
  allocate(g33_hel(nr, nchi))
  allocate(chi(nchi))
  allocate(chikn(nchi))
  allocate(cs(nr))
  allocate(qs(nr))
  allocate(dqs(nr))
  allocate(p0(nr))
  allocate(rbphi(nr))
  allocate(dp0(nr))
  allocate(drbphi(nr))

  write(*, *) 'pre: mapping'
  !call system('date +%s%N')
  call equilibrium_constructor(equilibrium_out(1), nr, nchi)
  call mapping(cx, cy, xaxis, yaxis, a, equilibrium_out(1))
  write(*, *) 'post: mapping'
  !call system('date +%s%N')
  call current_at_nodes(a, equilibrium_out(1))
  if (xmgrace_output) then
    call plot_data_2d('3d', equilibrium_out(1)%coord_sys%position%r, &
     equilibrium_out(1)%coord_sys%position%z, 'idl_jphi_x_y', &
     equilibrium_out(1)%profiles_2d(1)%jphi)
    call plot_data_2d('3d', equilibrium_out(1)%coord_sys%position%r, &
     equilibrium_out(1)%coord_sys%position%z, 'idl_jpar_x_y', &
     equilibrium_out(1)%profiles_2d(1)%jpar)
  end if

!-- deallocate pskn1
  deallocate(pskn1)

  write(*, *) 'pre: w_mhd'
  !call system('date +%s%N')
!-- calculate W_MHD
  if (verbosity > 2) then
    mhd_energy = w_mhd(equilibrium_out(1), rminor, eps)
    volume = plasma_volume(equilibrium_out(1), rminor, eps)
    write(iu6, *) 'W_MHD = ', mhd_energy, ' J'
    write(iu6, *) 'plasma volume = ', volume, ' m^3'
  end if

!-- calculate ballooning stability

  allocate(dimerc(nr))
  allocate(drmerc(nr))
  allocate(hh(nr))
  allocate(qprof(nr))
  allocate(dqprof(nr))
  allocate(geonc(nr))
  allocate(zjpar(nr))

  zjpar = 0._r8  ! initialize array

  if (diagnostics_on) then
    call mercier(a, dimerc, drmerc, hh, qprof, dqprof, geonc, zjpar)
    call suydam_ballooning(equilibrium_out(1), zvol, zvolp, xaxis)
  end if

  deallocate(g11_hel, g12_hel, g33_hel)
  deallocate(chi, chikn)
  deallocate(cs)
  deallocate(qs, dqs)
  deallocate(p0, dp0, rbphi, drbphi)
  deallocate(zvol, zvolp)

  call passing_particles(a, fcirc, b02av, b0max, rav, equilibrium_out(1))

  if (diagnostics_on) then
    call si_output(a, rvac, bvac, fcirc, qprof, dqprof, rav, geonc, &
     zjpar, dimerc, drmerc, hh, current)
  end if

  deallocate(psikn, dpsikn, ddpsikn, radpsi)
  deallocate(dimerc, drmerc)
  deallocate(hh)
  deallocate(dqprof)
  deallocate(geonc)
  deallocate(zjpar)
  deallocate(fcirc)
  deallocate(b02av)
  deallocate(b0max)
  deallocate(rav)

!-- write eqdsk file
  if ((trim(output) == 'equilibrium' .or. trim(output) == 'full') &
   .and. eqdsk_file) then
    call write_eqdsk(a, ip, rvac, bvac, qprof, path)
  end if

!TODO
!  deallocate(dp)
!  deallocate(fdf)
!  deallocate(j_tor)
!  deallocate(q_in)

  deallocate(qprof)

  deallocate(psivec)
  deallocate(dpres, p_int, dgam, gam_int)

  call deallocate_spline_coefficients(gam_spline)
  call deallocate_spline_coefficients(dp_spline)
  call deallocate_spline_coefficients(j_spline)

  deallocate(psi, xx, yy)
  deallocate(nodeno)

  call equilibrium_destructor(equilibrium)
  if (allocated(dp)) deallocate(dp)
  if (allocated(fdf)) deallocate(fdf)
  if (allocated(j_tor)) deallocate(j_tor)
  if (allocated(q_in)) deallocate(q_in)
  if (allocated(psi_in)) deallocate(psi_in)
  deallocate(s_acc, sig, weights)

!-- format statements
  2 format('a   : ', e12.4)
  3 format('b   : ', e14.6)
  4 format('alfa   : ', e14.6)
  5 format('radius : ', e14.6)
  6 format('b0     : ', e14.6)
  7 format('cpsurf : ', e14.6)
 13 format('magnetic axis : x   = ', e14.6, '  y   = ', e14.6)

  write(*, *) 'helena end'
  !call system('date +%s%N')

  if (standard_output) close(out_he)

  return

end subroutine helena 

!***********************************************************************