helena21_mod.f90

Go to the documentation of this file.

module helena_grid
  use itm_types
  implicit none
  integer             :: nr, np
  real (r8),pointer      :: RR(:,:), ZZ(:,:), PSI(:,:)
  real (r8)              :: R_min, R_max, Z_min, Z_max
endmodule helena_grid

module helena_boundary
  use itm_types
  implicit none
  integer             :: mf, n_bnd
  real (r8), allocatable :: fr(:), R_bnd(:), Z_bnd(:)
  real (r8)              :: Bgeo, Rgeo, Zgeo, amin, eps, ellip, tria_u, tria_l, quad_u, quad_l
  real (r8)              :: Reast, Rwest
end module helena_boundary

module helena_profiles
  use itm_types
  implicit none
  integer             :: n_prof
  real (r8),allocatable  :: dp_dpsi(:),fdf_dpsi(:),p_psi(:),f_psi(:),zjz_psi(:),q(:)
  real (r8)              :: p_bnd
endmodule helena_profiles

module helena_elements
  use itm_types
  implicit none
  integer             :: n_elms, n_nodes
  integer,allocatable :: elm_list(:,:), index_list(:), index_inverse(:)
  real (r8)              :: RS(4,2)
  integer             :: IJ(4,2)
  real (r8)              :: H(4,4,4,4),   HS(4,4,4,4),  HR(4,4,4,4)
  real (r8)              :: HRS(4,4,4,4), HRR(4,4,4,4), HSS(4,4,4,4)
end module helena_elements

module helena_matrix
  use itm_types
  implicit none
  integer             :: lda, n_matrix, n_diag
  real (r8),allocatable  :: AA(:,:), BB(:)
end module helena_matrix

module helena_axis
  use itm_types
  implicit none
  real (r8)              :: R_axis, Z_axis, ps_axis, r_ax, s_ax
  real (r8)              :: q_axis, B_axis, CRR_axis, CZZ_axis
  integer             :: ielm_ax
end module helena_axis

module helena_map
  use itm_types
  implicit none
  integer :: n_psi_map, n_chi_map
  real (r8), allocatable :: gem11(:,:), gem12(:,:), gem22(:,:), gem33(:,:)
  real (r8), allocatable :: q_map(:), dq_map(:), p_map(:), dp_map(:), f_map(:), df_map(:)
  real (r8), allocatable :: RR_map(:,:),ZZ_map(:,:)
end module helena_map

module helena_surfaces
  use itm_types
  implicit none
  integer   :: n_pieces_max
  parameter(n_pieces_max = 1000)

  type type_surface
    integer :: n_pieces
    integer :: elm(n_pieces_max)
    real (r8)    :: r(4,n_pieces_max), s(4,n_pieces_max)   ! 4 variables per line piece of the flux surface
  endtype

  integer             :: n_psi
  real (r8) ,allocatable :: psi_values(:)
  type (type_surface), allocatable :: flux_surfaces(:)

end module helena_surfaces

module helena_phys
  use itm_types
  implicit none
  real (r8) :: beta_p, beta_t, beta_n, current, volume, area
end module helena_phys

module helena_gauss_4
  use itm_types
  implicit none
  real (r8) :: xgauss(4),wgauss(4)
  data xgauss /-0.861136311594053, -0.339981043584856, 0.339981043584856,  0.861136311594053 /
  data wgauss / 0.347854845137454,  0.652145154862546, 0.652145154862546,  0.347854845137454 /
end module helena_gauss_4

module helena_constants
  use itm_types
  implicit none
  real (r8)       :: PI, MU_zero
  parameter( pi      = 3.14159265358979323846264338327950288419716939937510 )
  parameter( mu_zero = pi * 4.d-7)
end module helena_constants

module helena21_mod

contains

function fdfdpsi(psi_n)
use helena_profiles
use itm_types
implicit none
real (r8)    :: fdfdpsi, psi_n, dpsi
integer :: n_int

dpsi = 1. /float(n_prof - 1)

n_int = max(int((n_prof-1)*psi_n) + 1,1)
n_int = min(n_int,n_prof-1)

fdfdpsi = ( fdf_dpsi(n_int) + (psi_n - dpsi * (n_int-1))/dpsi * (fdf_dpsi(n_int+1)-fdf_dpsi(n_int)) )

return
end function fdfdpsi

function dpdpsi(psi_n)
use helena_profiles
use itm_types
implicit none
real (r8)    :: dpdpsi, psi_n, dpsi, tmp
integer :: n_int

dpsi = 1. /float(n_prof - 1)

n_int = max(int((n_prof-1)*psi_n) + 1,1)
n_int = min(n_int,n_prof-1)

dpdpsi = ( dp_dpsi(n_int) + (psi_n - dpsi * (n_int-1))/dpsi * (dp_dpsi(n_int+1)-dp_dpsi(n_int)) )

return
end function dpdpsi

function fdia(psi_n)
use helena_profiles
use helena_axis
use helena_boundary
use itm_types
implicit none
real (r8)    :: f2int, f0, fdia, psi_n, dpsi
integer :: n_int

if (psi_n .lt. 0.) write(*,*) ' FDIA : ',psi_n

dpsi = 1. /float(n_prof - 1)

n_int = max(int((n_prof-1)*psi_n) + 1,1)
n_int = min(n_int,n_prof-1)

f2int = ( f_psi(n_int) + (psi_n - dpsi * (n_int-1))/dpsi * (f_psi(n_int+1)-f_psi(n_int)) )

f0 = rgeo * bgeo

fdia = sqrt( abs(f0*f0 + 2. * f2int) )

return
end function fdia

function pressure(psi_n)
use helena_profiles
use helena_axis
use itm_types
implicit none
real (r8)    :: pressure, psi_n, dpsi
integer :: n_int

dpsi = 1. /float(n_prof - 1)

n_int = max(int((n_prof-1)*psi_n) + 1,1)
n_int = min(n_int,n_prof-1)

pressure = ( p_psi(n_int) + (psi_n - dpsi * (n_int-1))/dpsi * (p_psi(n_int+1)-p_psi(n_int)) )

return
end function pressure

function current_profile(psi_n)
use helena_profiles
use itm_types
implicit none
real (r8)    :: current_profile, psi_n, dpsi
integer :: n_int

dpsi = 1. /float(n_prof - 1)

n_int = max(int((n_prof-1)*psi_n) + 1,1)
n_int = min(n_int,n_prof-1)

current_profile = zjz_psi(n_int) + (psi_n - dpsi * (n_int-1))/dpsi * (zjz_psi(n_int+1)-zjz_psi(n_int))

return
end function current_profile
!****************************************************************************************************************

subroutine helena21(R_bnd_in,Z_bnd_in,n_bnd_in,                                                                  &
                    rgeo_in,zgeo_in, amin_in,                                                                    &
                    ellip_in,tria_u_in,tria_l_in,quad_u_in,quad_l_in,                                            &
                    bgeo_in,                                                                                     &
                    psi_in, pprime_in, ffprime_in,                                                               &
                    pressure_in, fdia_in, current_in, n_prof_in, iopt_p, iopt_f,                                 &
                    nr_grid, np_grid,                                                                            &
                    psi_out, pprime_out, ffprime_out,                                                            &
                    pressure_out, fdia_out, current_out, qprof_out, vprime,                                      &
                    fraction_circ, moments, n_moments,                                                           &
                    surface_powers, surface_integrals, n_var_surfaces, n_int_surfaces,                           &
                    volume_powers, volume_integrals, n_var_volumes, n_int_volumes, n_psi_out, n_tht_out,         &
                    amin_out, rgeo_out, zgeo_out,                                                                &
                    area_out, volume_out, betap_out, xip_out, xli_out, beta_out,                                 &
                    r_axis_out, z_axis_out, b_axis_out, psi_axis_out, psi_bnd_out,                               &
                    rrflux,zzflux,psflux)

!****************************************************************************************************************
!
! Helena version 2.0
!
! - removed all normalisations. Quantities are in S.I. units.
! - improved(?) the tracing of flux surfaces (no iterations, should be even more reliable)
!
!
! author  : Guido Huysmans (Association Euratom-CEA)
! date    : 05/02/2009
! version : 0.8
! status  : slowly progressing ...
!
!
! to be done :
!        - adjust radial integration of dpdpsi and fdfdpsi to changing ps_axis (done?)
!        - change modules to derived types and arguments to subroutines
!        - do the mapping to straight field line coordinates (generalised including Boozer, Hamada)
!        - finding all the new bugs...
!
!    !!!    - fdia and pressure spline if available, otherwise integrate from spline pprime,ffprime
!
!****************************************************************************************************************
use helena_grid
use helena_boundary
use helena_profiles
use helena_matrix
use helena_axis
use helena_elements
use helena_surfaces
use helena_phys
use helena_map
use helena_constants

use itm_types
implicit none

integer          :: n_int_surfaces, n_int_volumes, n_var_surfaces, n_var_volumes, n_moments, n_psi_out, n_tht_out, ioption
real (r8)           :: rrflux(4,*),zzflux(4,*),psflux(4,*)
real (r8)           :: r_bnd_in(*),z_bnd_in(*),pressure_in(*),fdia_in(*),current_in(*),psi_in(*), pprime_in(*), ffprime_in(*)
real (r8)           :: bgeo_in,rgeo_in,zgeo_in,amin_in,ellip_in,tria_u_in,tria_l_in,quad_u_in,quad_l_in
integer          :: n_bnd_in, n_prof_in, nr_grid, np_grid, iopt_p, iopt_f
real (r8)           :: surface_powers(n_var_surfaces,n_int_surfaces),surface_integrals(n_psi_out,n_int_surfaces)
real (r8)           :: volume_powers(n_var_volumes,n_int_volumes),   volume_integrals(n_psi_out,n_int_volumes)
real (r8)           :: psi_n, psi_out(*), pressure_out(*), fdia_out(*), current_out(*), pprime_out(*), ffprime_out(*)
real (r8)           ::  qprof_out(*), vprime(*), fraction_circ(*),moments(n_moments,*)
integer          :: i, j, index, n_iter, ifail, inode, n_chi, ileft, iright
real (r8)           :: s, error_iteration, error_large, error_current, small, ss, psi_plot
real (r8)           :: fdf_error, t1, t2, amix, fmix, ymin, ymax
logical          :: solve_only, rhs_only
!real (r8)           :: pressure, fdia, fdfdpsi, dpdpsi, current_profile
real (r8)           :: amin_out, rgeo_out, zgeo_out, area_out, volume_out, betap_out, xip_out, xli_out, beta_out
real (r8)           :: r_axis_out, z_axis_out, b_axis_out, psi_axis_out, psi_bnd_out
real (r8),allocatable :: rho_out(:), xplot(:), qplot(:), pplot(:)
real (r8)           :: qtmp

call begplt('helena.ps')

!----------------------------------------- initialise
write(*,*) '**********************************************'
write(*,*) '*          HELENA 2.0 (alpha_21)             *'
write(*,*) '**********************************************'

write(*,*) ' nr, np        : ',nr_grid,np_grid
write(*,*) ' n_var_surfaces,n_int_surfaces : ',n_var_surfaces,n_int_surfaces
write(*,*) ' n_var_volumes, n_int_volumes  : ',n_var_volumes, n_int_volumes
write(*,*) ' n_psi, n_tht  : ',n_psi_out,n_tht_out

mf     = 128
nr     = nr_grid
np     = 2*int((np_grid+1)/2)
n_iter = 50
n_chi  = 64
error_iteration = 1.d-6
error_large     = 1.d-4
error_current   = 1.d-4

bgeo = bgeo_in

write(*,*) ' Bgeo : ',bgeo

n_bnd = n_bnd_in
allocate(r_bnd(n_bnd),z_bnd(n_bnd))
r_bnd(1:n_bnd) = r_bnd_in(1:n_bnd)
z_bnd(1:n_bnd) = z_bnd_in(1:n_bnd)

call fshape

n_prof = n_prof_in


call initialise_profiles(n_prof_in,iopt_p,iopt_f,psi_in, pprime_in, ffprime_in,pressure_in, fdia_in, current_in)

call initialise_grid

call initialise_elements

solve_only = .false.
rhs_only   = .false.
amix = 0.e0_r8
fmix = 0.e0_r8

!---------------------------------------- iteration for current profile
if (iopt_f .eq. 3) then

  do i=1, n_iter

    call gs_solve(n_iter,rhs_only,solve_only,error_large,amix,ifail)

    if (ifail .ne. 0 ) then
      amix = (1.e0_r8 + amix) /2.e0_r8
      write(*,'(A,f8.3)') ' increasing amix : ',amix
    endif

    call update_fdf(fmix,fdf_error)

    write(*,'(A,e14.6)') ' fdf_error : ',fdf_error
    if (fdf_error .lt. error_current) exit

  enddo

endif
  
n_iter     = 100
solve_only = .false.
rhs_only   = .false.

call gs_solve(n_iter,rhs_only,solve_only,error_iteration,amix,ifail)

!-------------------------------------- output quantities
call phys_values

!call plot_solution

n_psi = n_psi_out - 1           ! n_psi does not count the axis (only the flux surfaces traced)

if (allocated(psi_values))    deallocate(psi_values)
if (allocated(flux_surfaces)) deallocate(flux_surfaces)

allocate(psi_values(n_psi))

small = 1.d-8
write(*,*) 'main : i ss psi_values'
do i=1,n_psi
  ss = float(i)/float(n_psi)
  psi_values(i) = ps_axis - (1.e0_r8 - small) * ss*ss * ps_axis
  write(*,'(A,i4,2e14.6)') ' main : ',i,ss,psi_values(i)
enddo

call find_flux_surfaces
call plot_flux_surfaces

!write(*,*) ' size surface_integrals : ',size(surface_integrals)
!call fluxsurface_integrals(surface_powers,n_int_surfaces,n_var_surfaces,surface_integrals)
!write(*,*) ' fluxsurface integrals (1) : '
!do i=1,n_psi_out
!  write(*,'(i5,15e12.4)') i,surface_integrals(i,1:n_int_surfaces)
!enddo

write(*,*) 'profiles out : i psi_n psi_out fdia_out'
do i=1,n_psi+1
  if (i.eq.1) then
    psi_n = 0.e0_r8
    psi_out(i)      = - ps_axis + psi_in(n_prof_in) 
  else
    psi_n           = -(psi_values(i-1) - ps_axis) / ps_axis
    psi_out(i)      = - psi_values(i-1) + psi_in(n_prof_in) 
  endif
  pprime_out(i)   = dpdpsi(psi_n)
  ffprime_out(i)  = fdfdpsi(psi_n)
  pressure_out(i) = pressure(psi_n)
  fdia_out(i)     = fdia(psi_n)
  current_out(i)  = current_profile(psi_n)

  write(*,'(A,i4,8e14.6)') 'profiles out : ',i,psi_n,psi_out(i),fdia_out(i)
enddo

!call export_helena

call helena_remesh(n_psi_out,n_tht_out,rrflux,zzflux,psflux)

call helena_flux_surface_integrals(n_psi_out,n_tht_out,rrflux,zzflux,psflux, &
                                   surface_powers,n_int_surfaces,n_var_surfaces,surface_integrals, &
                                   qprof_out, vprime) 

do i=1, n_psi_out
  qprof_out(i)   = fdia_out(i) * qprof_out(i) / (2.e0_r8*pi)
enddo

write(*,*) ' fluxsurface integrals (2) : ',n_int_surfaces
do i=1,n_psi_out
  write(*,'(i5,15e12.4)') i,surface_integrals(i,1:n_int_surfaces),qprof_out(i),vprime(i)
enddo
      
!call helena_mapping(RRflux,ZZflux,PSflux,n_tht_out,n_chi)

call helena_volume_integrals(n_psi_out,n_tht_out,rrflux,zzflux,psflux, &
                             volume_powers,volume_integrals,n_var_volumes,n_int_volumes,n_psi+1)

allocate(rho_out(n_psi_out))

write(*,*) 'i psi_out rho_out volume_integrals(1:3) qprof_out qtmp abs(qtmp-qprof_out)/qprof_out fdia_out'
do i=1, n_psi_out
  
  rho_out(i) = sqrt(volume_integrals(i,2)/ (2.e0_r8*pi) / (pi * bgeo_in))          ! B at magnetic axis?
  
  qtmp = 0.e0_r8
  if ( (i .gt. 1) .and. (i .lt. n_psi_out)) then
    qtmp = -(volume_integrals(i+1,2) - volume_integrals(i-1,2)) / (4.*pi*pi) / (psi_out(i+1)-psi_out(i-1))
  endif
  
  write(*,'(i5,20f12.6)') i,psi_out(i),rho_out(i),volume_integrals(i,1:3),qprof_out(i),qtmp,abs(qtmp-qprof_out(i))/qprof_out(i),&
  fdia_out(i)
  
enddo


!call helena_moments(RRflux,ZZflux,n_psi_out,n_tht_out,moments,n_moments)

call helena_circulating(n_psi_out,n_tht_out,rrflux,zzflux,psflux,fraction_circ)
fraction_circ(1) = 1.e0_r8

area_out   = area
volume_out = volume
betap_out  = beta_p
xip_out    = current
xli_out    = 0.e0_r8
beta_out   = beta_t
r_axis_out = r_axis
z_axis_out = z_axis
b_axis_out = b_axis
psi_axis_out = ps_axis
psi_bnd_out  = psi_out(n_psi_out)
amin_out     = amin
rgeo_out     = rgeo
zgeo_out     = zgeo


allocate(xplot(2*n_psi_out),pplot(2*n_psi_out),qplot(2*n_psi_out))

do i=1,n_psi_out
  ileft  = (i-1)*n_tht_out + n_tht_out /2.e0_r8 + 1
  iright = (i-1)*n_tht_out + 1  
  xplot(n_psi_out+i)   = rrflux(1,ileft)
  xplot(n_psi_out-i+1) = rrflux(1,iright)
  pplot(n_psi_out+i)   = pressure_out(i)
  pplot(n_psi_out-i+1) = pressure_out(i)
  qplot(n_psi_out+i)   = qprof_out(i)
  qplot(n_psi_out-i+1) = qprof_out(i)
enddo

call lincol(0)
call lplot(3,2,1,xplot,pplot,2*n_psi_out,1,'pressure [Pa]',13,'R [m]',5,' ',1)
!DPC infinite loop? call lplot(3,3,1,xplot,qplot,2*n_psi_out,1,'q-profile',9,'R [m]',5,' ',1)

call lblbot('HELENA_2(alpha)',15)

call finplt

deallocate(aa,bb,fr,dp_dpsi,fdf_dpsi,p_psi,f_psi)
deallocate(elm_list, index_list, index_inverse)
deallocate(rr, zz, psi)
!DPC
deallocate(flux_surfaces,r_bnd,z_bnd,psi_values,zjz_psi)

end subroutine helena21

subroutine helena_mapping(RR,ZZ,PSI,n_tht,n_chi)
!-----------------------------------------------------------------------
! SUBROUTINE TO CALCULATE THE METRIC COEFFICIENTS NEEDED FOR CASTOR/MISHKA
! a straight field line flux surface coordinate system with a regular
! toroidal angle.
!-----------------------------------------------------------------------
use helena_constants
use helena_axis
use helena_profiles
use helena_surfaces
use helena_map
use helena_gauss_4
use itm_types
implicit none

real (r8)               :: rr(4,*),zz(4,*),psi(4,*)
real (r8), allocatable  :: cchi(:,:), schi(:,:), xchi(:,:), ychi(:,:), chikn(:)
real (r8), allocatable  :: s_values(:)
integer, allocatable :: ijchi(:,:)
logical              :: chin

real (r8)               :: sumq, sumqr, zsumq, zsumqr, psi_n, psir, psirr, ejac, er, bigr
real (r8)               :: s, maxerr, x, xr, xs, xrs, xrr, xss, y, yr, ys, yrs, yrr, yss
real (r8)               :: dum, dum1, dum2, dum3, dum4, qq, dqq, rb, drb
real (r8)               :: dchi, chir, chis, chirs, chirr, psix, psiy, chix, chiy, zchi
real (r8)               :: errj, serr, cerr, chiss, ejac2, a0, a1, a2, a3, s1, s2, s3
integer              :: n_tht, n_chi, i,j,k, i_elm, node1, node2, node3, node4, no, nom, nop
integer              :: jbase, n1, n2, n3, n4, ierr, jerr, ifail
!real (r8)               :: fdia, fdfdpsi
real (r8) :: psi_tmp,psir_tmp,psis_tmp,psirs_tmp,psirr_tmp,psiss_tmp

write(*,*) '********************************'
write(*,*) '*        helena mapping        *'
write(*,*) '********************************'
write(*,*) ' n_tht, n_chi : ',n_tht, n_chi

!---------------------------------------------------------------------
maxerr = -1.d20


allocate(s_values(n_psi+1),cchi(4,(n_psi+1)*n_tht),schi((n_psi+1),n_chi),xchi((n_psi+1),n_chi),ychi((n_psi+1),n_chi))
allocate(chikn(n_chi),ijchi((n_psi+1),n_chi))
allocate(q_map(n_psi+1),dq_map(n_psi+1))

do i=1,n_psi
  s_values(i) = float(i)/float(n_psi)
enddo

!------------------------------------------ - Q PROFILE ----------------
!                                           - DQ/DS PROFILE
!                                           - CHI VALUES AT THETA NODES
do i = 2, n_psi+1

  sumq  = 0.e0_r8
  sumqr = 0.e0_r8
  psi_n = s_values(i-1)**2

  psir  = - ps_axis * 2.e0_r8*s_values(i-1) / (2.e0_r8*REAL(n_psi))
  psirr = - ps_axis * 2.e0_r8               / (2.e0_r8*REAL(n_psi))**2

  do j = 1, n_tht - 1

    node1 = (i-2)*n_tht + j
    node2 = node1 + 1
    node3 = node2 + n_tht
    node4 = node1 + n_tht
    
    write(*,'(A,i5,4e16.8)') ' node 1 : ',node1,psi(:,node1)
    write(*,'(A,i5,4e16.8)') ' node 2 : ',node2,psi(:,node2)
    write(*,'(A,i5,4e16.8)') ' node 3 : ',node3,psi(:,node3)
    write(*,'(A,i5,4e16.8)') ' node 4 : ',node4,psi(:,node4)

    do k = 1,4

      s = xgauss(k)

      CALL interp(rr(1,node1),rr(1,node2),rr(1,node3),rr(1,node4),+1.e0_r8,s, x,xr,xs,xrs,xrr,xss)
      CALL interp(zz(1,node1),zz(1,node2),zz(1,node3),zz(1,node4),+1.e0_r8,s, y,yr,ys,yrs,yrr,yss)

CALL interp(psi(1,node1),psi(1,node2),psi(1,node3),psi(1,node4),1.e0_r8,s,psi_tmp,psir_tmp,psis_tmp,psirs_tmp,psirr_tmp,psiss_tmp)
  write(*,'(A,i5,12e12.4)') ' 1 : ',i,s_values(i-1),s,ps_axis,psi_tmp,psir_tmp,psir,psirr_tmp,psirr

      ejac  = xr*ys - xs*yr
      er    = xrr*ys + xr*yrs - xrs*yr - xs*yrr
      bigr  = x
      sumq  = sumq  - wgauss(k) * ejac / ( bigr * abs(psir))
      sumqr = sumqr + psirr * ejac / ((psir**2)*bigr)* wgauss(k)
      sumqr = sumqr - er / (bigr*psir)               * wgauss(k)
      sumqr = sumqr + ejac * xr / ((bigr**2)*psir)   * wgauss(k)

    enddo

    cchi(1,(i-1)*n_tht+j+1) = sumq
    cchi(2,(i-1)*n_tht+j+1) = sumqr

    CALL interp(rr(1,node1),rr(1,node2),rr(1,node3),rr(1,node4),+1.e0_r8,1.e0_r8,x,xr,xs,xrs,xrr,xss)
    CALL interp(zz(1,node1),zz(1,node2),zz(1,node3),zz(1,node4),+1.e0_r8,1.e0_r8,y,yr,ys,yrs,yrr,yss)


    CALL interp(psi(1,node1),psi(1,node2),psi(1,node3),psi(1,node4),  &
         1.e0_r8,1.e0_r8,psi_tmp,psir_tmp,psis_tmp,psirs_tmp,psirr_tmp,psiss_tmp)
    write(*,'(A,i5,12e12.4)') ' end : ',i,s_values(i-1),ps_axis,psi_tmp,  &
         psir_tmp,psir,psirr_tmp,psirr

  
    ejac   = xr*ys - xs*yr
    er     = xrr*ys + xr*yrs - xrs*yr - xs*yrr
    bigr   = x

    zsumq  = - ejac / ( bigr * abs(psir))
    zsumqr = + psirr * ejac / (psir**2 *bigr) - er / (bigr*psir) + ejac*xr/((bigr**2) * psir)

    cchi(3,(i-1)*n_tht+j+1) = zsumq
    cchi(4,(i-1)*n_tht+j+1) = zsumqr

  enddo

  cchi(1,(i-1)*n_tht+1) = 0.e0_r8
  cchi(2,(i-1)*n_tht+1) = 0.e0_r8

  node1 = (i-1)*n_tht + 1
  node2 = node1 + 1
  node3 = node2 + n_tht
  node4 = node1 + n_tht

  CALL interp(rr(1,node1),rr(1,node2),rr(1,node3),rr(1,node4),+1.e0_r8,-1.e0_r8,x,xr,xs,xrs,xrr,xss)
  CALL interp(zz(1,node1),zz(1,node2),zz(1,node3),zz(1,node4),+1.e0_r8,-1.e0_r8,y,yr,ys,yrs,yrr,yss)

  ejac   = xr*ys - xs*yr
  er     = xrr*ys + xr*yrs - xrs*yr - xs*yrr
  bigr   = x

  zsumq  = - ejac / (bigr * abs(psir))
  zsumqr = + psirr * ejac / (psir**2 *bigr) - er / (bigr*psir) + ejac*xr/((bigr**2) * psir)

  cchi(3,(i-1)*n_tht+1) =  zsumq
  cchi(4,(i-1)*n_tht+1) =  zsumqr

  q_map(i)  = 0.5e0_r8/pi *   sumq  * fdia(psi_n)
  dq_map(i) = 0.5e0_r8/pi * ( sumqr * fdia(psi_n) + sumq*fdfdpsi(psi_n)/fdia(psi_n) /(2.e0_r8*(n_psi-1)))

  write(*,*) i,q_map(i),dq_map(i)

enddo

q_map(1)  = fdia(0.e0_r8) /(2.e0_r8*sqrt(crr_axis*czz_axis)*r_axis)
dq_map(1) = 0.e0_r8

do i = 2, n_psi

  psi_n = s_values(i)**2

  do j = 1, n_tht

    dum = cchi(1,i*n_tht)

    no  = (i-1)*n_tht + j

    cchi(1,no) = 2.e0_r8*pi*cchi(1,no) / dum
    dum2       = cchi(2,i*n_tht)
    cchi(2,no) = 2.e0_r8*pi*cchi(2,no) / dum
    cchi(3,no) = 2.e0_r8*pi*cchi(3,no) / dum
    cchi(4,no) = 2.e0_r8*pi*cchi(4,no) / dum

    qq  = q_map(i)
    dqq = dq_map(i)
    rb  = fdia(psi_n)
    drb = fdfdpsi(psi_n) / rb  /  (2*(n_psi-1))

    cchi(2,no)= +(dqq/qq - drb/rb) * cchi(1,no) - cchi(2,no)
    cchi(4,no)= +(dqq/qq - drb/rb) * cchi(3,no) - cchi(4,no)

  enddo

enddo

WRITE(*,'(A,f8.4)') '  Q ON AXIS     = ',q_map(1)
WRITE(*,'(A,f8.4)') '  Q AT BOUNDARY = ',q_map(n_psi)

!--------------------------- DETERMINE POSITIONS OF EQUIDISTANT CHI'S --
!                            AND CALCULATE MATRIX ELEMENTS -------------

allocate(gem11(n_psi,n_chi),gem12(n_psi,n_chi),gem22(n_psi,n_chi),gem33(n_psi,n_chi))
allocate(rr_map(n_psi,n_chi), zz_map(n_psi,n_chi))

do j=1,n_chi
  chikn(j) = 2.e0_r8 * pi * REAL(j-1)/REAL(n_chi)
enddo

do i=1,n_psi-1

  psi_n = s_values(i+1)**2
  psir  = - ps_axis * 2.e0_r8 * s_values(i+1) / (2.e0_r8*REAL(n_psi))

!-------------------------- FIRST POINT IS KNOWN -----------------------
  k = 1
  s = -1.e0_r8
  schi(i,1)  = s
  ijchi(i,1) = k

  node1 = (i-1)*n_tht + k
  node2 = node1 + 1
  node3 = node2 + n_tht
  node4 = node1 + n_tht

  CALL interp(rr(1,node1),rr(1,node2),rr(1,node3),rr(1,node4),1.e0_r8,s,xchi(i,1),xr,xs,xrs,xrr,xss)
  CALL interp(zz(1,node1),zz(1,node2),zz(1,node3),zz(1,node4),1.e0_r8,s,ychi(i,1),yr,ys,yrs,yrr,yss)
  CALL interp(cchi(1,node1),cchi(1,node2),cchi(1,node3),cchi(1,node4),1.e0_r8,s,dchi,chir,chis,chirs,chirr,chiss)
  
  CALL interp(psi(1,node1),psi(1,node2),psi(1,node3),psi(1,node4),1.e0_r8,s,psi_tmp,psir_tmp,psis_tmp,psirs_tmp,psirr_tmp,psiss_tmp)

  write(*,'(A,i5,8e12.4)') ' 2 : ',i,s_values(i+1),ps_axis,psi_tmp,psir_tmp

!-----------------------------------------------------------------------
  ejac  = (xr*ys-xs*yr)
  ejac2 = ejac**2

!----------------------------- DATA FOR VECTOR PLOT TO FILE 21 -----
  psix = psir * ys / ejac
  psiy = - psir * xs / ejac
  chix = (chir * ys - chis * yr) / ejac
  chiy = (-chir * xs + chis * xr)/ ejac
!        WRITE(21,61) SQRT(ZPSI),DCHI,XCHI(NO),YCHI(NO),PSIX,PSIY,CHIX,CHIY

!----------------------------------------------------------------------
  gem11(i,1) = psir**2 * (xs**2 + ys**2) / ejac2
  gem12(i,1) = ( chir * psir * (xs**2 + ys**2) - chis * psir * (xr*xs + yr*ys) ) / ejac2
  gem33(i,1) = xchi(i,1)
  rr_map(i,1)= xchi(i,1)
  zz_map(i,1)= ychi(i,1)

!------------------------------------ CHECK JACOBIAN -----------------
  dum1 = fdia(psi_n)**2 / (q_map(i)**2 * gem33(i,1))
  dum2 = gem11(i,1) * (chix**2 + chiy**2)
  dum3 = gem12(i,1)
  dum4 = dum2 - dum3*dum3
  errj = abs(dum4-dum1)

  IF (errj.GT.maxerr) THEN
    maxerr=errj
    ierr = i
    jerr = j
    serr = s_values(i)
    cerr = dchi
  ENDIF

!---------------------------------------------------------------------
  jbase = 2

  do k = 1,n_tht-1

    chin =.false.

    do j = jbase,n_chi

      zchi = chikn(j)

      IF (((cchi(1,(i-1)*n_tht+k).LE.zchi).AND.(cchi(1,(i-1)*n_tht+k+1) .GE.zchi)) ) THEN

        chin = .true.
        nom = (i-1)*n_tht + k
        nop = nom + 1
        a3 = (cchi(1,nom)+cchi(3,nom)-cchi(1,nop)+cchi(3,nop))/4.e0_r8
        a2 = (- cchi(3,nom) + cchi(3,nop))/4.e0_r8
        a1=(-3*cchi(1,nom)-cchi(3,nom)+3*cchi(1,nop)-cchi(3,nop))/4.e0_r8
        a0=( 2*cchi(1,nom)+cchi(3,nom)+2*cchi(1,nop)-cchi(3,nop))/4.e0_r8 - zchi

        CALL solvp3(a0,a1,a2,a3,s,s2,s3,ifail)

        IF (ifail.EQ.0) THEN

          schi(i,j) = s
          ijchi(i,j) = k

          node1 = (i-1)*n_tht + k
          node2 = node1 + 1
          node3 = node2 + n_tht
          node4 = node1 + n_tht

          CALL interp(rr(1,node1),rr(1,node2),rr(1,node3),rr(1,node4),-1.e0_r8,s,xchi(i,j),xr,xs,xrs,xrr,xss)
          CALL interp(zz(1,node1),zz(1,node2),zz(1,node3),zz(1,node4),-1.e0_r8,s,ychi(i,j),yr,ys,yrs,yrr,yss)
          CALL interp(cchi(1,node1),cchi(1,node2),cchi(1,node3),cchi(1,node4),-1.e0_r8,s,dchi,chir,chis,chirs,chirr,chiss)

!-----------------------------------------------------------------
          ejac  = (xr*ys-xs*yr)
          ejac2 = ejac**2

!-------------------------- DATA FOR VECTOR PLOT TO FILE 21 -----
          psix = psir * ys / ejac
          psiy = - psir * xs / ejac
          chix = (chir * ys - chis * yr) / ejac
          chiy = (-chir * xs + chis * xr)/ ejac

!          WRITE(21,61) SQRT(ZPSI),DCHI,XCHI(NO),YCHI(NO),PSIX,PSIY,CHIX,CHIY
!----------------------------------------------------------------

          gem11(i,j) =  psir**2 * (xs**2 + ys**2) / ejac2
          gem12(i,j) = ( chir * psir * (xs**2 + ys**2) - chis * psir * (xr*xs + yr*ys) ) / ejac2
          gem33(i,j) = xchi(i,j)**2
          rr_map(i,j) = xchi(i,j)
          zz_map(i,j) = ychi(i,j)

!------------------------------ CHECK JACOBIAN -----------------
          dum1 = fdia(psi_n)**2 / (q_map(i)**2 * gem33(i,j))    !!!!!! check if psi_n defined
          dum2 = gem11(i,j) * (chix**2 + chiy**2)
          dum3 = gem12(i,j)
          dum4 = dum2 - dum3*dum3
          errj = abs(dum4-dum1)
          IF (errj.GT.maxerr) THEN
            maxerr=errj
            ierr = i
            jerr = j
            serr = s_values(i)
            cerr = dchi
          ENDIF
!---------------------------------------------------------------------
        ELSE
          WRITE(20,*) 'ERROR IN SOLVP3 I,J,K : ',i,j,k,s,s2,s3
          WRITE(*,*) a0,a1,a2,a3,zchi
          WRITE(*,*) cchi(1,(i-1)*n_tht+k),cchi(1,(i-1)*n_tht+k+1)
        ENDIF
      ELSEIF (chin) THEN
        jbase = j
        exit
      ENDIF

    enddo
  enddo
enddo

WRITE(*,*)
WRITE(*,*) '***************************************************'
WRITE(*,'(3e12.4,2i5)') maxerr,serr,cerr,ierr,jerr
WRITE(*,*) '***************************************************'

n_psi_map = n_psi
n_chi_map = n_chi

RETURN
END subroutine helena_mapping


subroutine helena_moments(RRflux,ZZflux,nr_flux,np_flux,moments,n_moments)
!----------------------------------------------------------------------
! subroutine finds the Shafranov shift, ellipticity and triangularity
! as a function of radius (flux).
!
! moments(1,:) : minor radius (R_east - R_west) / 2.
! moments(2,:) : Shafranov shift (R_east + R_west)/2 - R_geo
! moments(3,:) : (Z_east - Z_west)/2. - Z_geo
! moments(4,:) : Z_top/((R_east - R_west)/2.)
! moments(5,:) : Z_low/((R_east - R_west)/2.)
! moments(6,:) : (R_top - (R_east + R_west)/2.)/((R_east - R_west)/2.)
! moments(7,:) : (R_low - (R_east + R_west)/2.)/((R_east - R_west)/2.)
!----------------------------------------------------------------------
use helena_constants
use helena_axis
use itm_types
implicit none
real (r8)            :: rrflux(4,*),zzflux(4,*),moments(n_moments,*)
integer           :: nr_flux, np_flux, n_moments
real (r8) r,sum1,sum2,sumr,b02max,s,ws,rrg1,drrg1_dr,drrg1_ds,drrg1_drs,drrg1_drr,drrg1_dss,  &
     zzg1,dzzg1_dr,dzzg1_ds,dzzg1_drs,dzzg1_drr,r_geo,z_geo,a_minor,  &
     zzm,zzmr,zzp,zzpr,aa,bb,cc,det,r_top,dummy,z_top,  &
     r_east,r_west,z_east,z_west,r_mid,z_mid,radius
integer i,j,n1,n2,n3,n4,ngs



r_geo   =    ( rrflux(1,(nr_flux-1)*np_flux + 1) + rrflux(1,(nr_flux-1)*np_flux + np_flux/2 + 1)) /2.
z_geo   =    ( zzflux(1,(nr_flux-1)*np_flux + 1) + zzflux(1,(nr_flux-1)*np_flux + np_flux/2 + 1)) /2.
a_minor = abs( rrflux(1,(nr_flux-1)*np_flux + 1) - rrflux(1,(nr_flux-1)*np_flux + np_flux/2 + 1)) /2.

write(*,*) ' moments : ',r_geo,z_geo

moments(n_moments,1:nr_flux) = 0.e0_r8

moments(1,1) = 0.e0_r8
moments(2,1) = rrflux(1,1) - r_geo
moments(3,1) = zzflux(1,1) - z_geo
moments(4,1) = sqrt(crr_axis / czz_axis)
moments(5,1) = sqrt(crr_axis / czz_axis)
moments(6,1) = 0.e0_r8
moments(7,1) = 0.e0_r8

do i = 2, nr_flux

  do j = 1, np_flux - 1

    n1 = (i-1)*np_flux + j
    n2 = (i-1)*np_flux + j + 1

    if (zzflux(3,n1)*zzflux(3,n2) .le. 0.e0_r8) then
!------------------------------------- QUAD. EQ FOR S VALUE AT MINIMUM -

      zzm  = zzflux(1,n1)
      zzmr = zzflux(3,n1)
      zzp  = zzflux(1,n2)
      zzpr = zzflux(3,n2)
      aa =  3.e0_r8 * (zzm + zzmr - zzp + zzpr ) / 4.e0_r8
      bb =  ( - zzmr + zzpr ) / 2.e0_r8
      cc =  ( - 3.e0_r8*zzm - zzmr + 3.e0_r8*zzp - zzpr) / 4.e0_r8
      det = bb*bb - 4.e0_r8*aa*cc
      s  = root(aa,bb,cc,det,1.e0_r8)

      if (abs(s) .GT. 1.e0_r8) then
        s = root(aa,bb,cc,det,-1.e0_r8)
      endif

      call cub1d(rrflux(1,n1),rrflux(3,n1),rrflux(1,n2),rrflux(3,n2),s,r_top,dummy)
      call cub1d(zzflux(1,n1),zzflux(3,n1),zzflux(1,n2),zzflux(3,n2),s,z_top,dummy)

      r_east  = rrflux(1,(i-1)*np_flux + 1)
      r_west  = rrflux(1,(i-1)*np_flux + np_flux/2 + 1)
      z_east  = zzflux(1,(i-1)*np_flux + 1)
      z_west  = zzflux(1,(i-1)*np_flux + np_flux/2 + 1)

      r_mid   = (r_east + r_west)/2.e0_r8
      z_mid   = (z_east + z_west)/2.e0_r8
      radius  = abs(r_east - r_west)/2.e0_r8

      write(*,'(A,i3,4f12.6)') 'mom : ',i,r_mid,z_mid,r_top,z_top

      moments(1,i) = radius
      moments(2,i) = (r_mid - r_geo)
      moments(3,i) = (z_mid - z_geo)

      if ( z_top - z_mid .gt. 0.e0_r8) then
        moments(4,i) =  (z_top - z_mid)  / radius
        moments(6,i) = -(r_top - r_mid)  / radius
      else
        moments(5,i) = -(z_top - z_mid)  / radius
        moments(7,i) = -(r_top - r_mid)  / radius
      endif

   endif
  enddo

enddo

!do i=1,nr_flux
!  write(*,'(i5,7f12.6)') i,moments(1:7,i)
!enddo

return
end subroutine helena_moments


subroutine helena_circulating(nr_flux,np_flux,RRflux,ZZflux,PSflux,fraction_circ)
!----------------------------------------------------------------------
! subroutine for the fraction of circulating (non-trapped) particles
!----------------------------------------------------------------------
use helena_gauss_4
use helena_axis
use helena_constants
use itm_types
implicit none
real (r8),allocatable:: b0max(:),b02av(:),rav(:),sumk(:)
real (r8)            :: rrflux(4,*),zzflux(4,*),psflux(4,*),fraction_circ(*)
integer           :: nr_flux, np_flux
!real (r8)            :: fdia
real (r8) r,sum1,sum2,sumr,b02max,s,ws,  &
     rrg1,drrg1_dr,drrg1_ds,drrg1_drs,drrg1_drr,drrg1_dss,  &
     zzg1,dzzg1_dr,dzzg1_ds,dzzg1_drs,dzzg1_drr,dzzg1_dss,  &
     psg1,dpsg1_dr,dpsg1_ds,dpsg1_drs,dpsg1_drr,dpsg1_dss,  &
     rzjac,gradps2,zjdchi,psi_n,bphi,b02,b0,bm,zlam,sum
integer i,j,n1,n2,n3,n4,ngs,nk,k

!---------------------- find Bmax and B^2 average on every surface

allocate(b0max(nr_flux),b02av(nr_flux),rav(nr_flux))

r = +1.e0_r8

do i=1, nr_flux - 1

  sum1 = 0.e0_r8
  sum2 = 0.e0_r8
  sumr = 0.e0_r8
  b02max = -1.d20

  do j=1, np_flux

    n1 = (i-1)*np_flux + j
    n2 = (i-1)*np_flux + j + 1
    n3 = (i  )*np_flux + j + 1
    n4 = (i  )*np_flux + j

    if (j .eq. np_flux) then
      n1 = (i  )*np_flux
      n2 = (i  )*np_flux - np_flux + 1
      n3 = (i  )*np_flux + 1
      n4 = (i  )*np_flux + np_flux
    endif

!------------------------------------- 4 POINT GAUSSIAN INT. IN S -----
    do ngs = 1, 4

      s  = xgauss(ngs)
      ws = wgauss(ngs)

      CALL interp(rrflux(1,n1),rrflux(1,n2),rrflux(1,n3),rrflux(1,n4),r,s,rrg1,drrg1_dr,drrg1_ds,drrg1_drs,drrg1_drr,drrg1_dss)
      CALL interp(zzflux(1,n1),zzflux(1,n2),zzflux(1,n3),zzflux(1,n4),r,s,zzg1,dzzg1_dr,dzzg1_ds,dzzg1_drs,dzzg1_drr,dzzg1_dss)
      CALL interp(psflux(1,n1),psflux(1,n2),psflux(1,n3),psflux(1,n4),r,s,psg1,dpsg1_dr,dpsg1_ds,dpsg1_drs,dpsg1_drr,dpsg1_dss)

      rzjac =  drrg1_dr * dzzg1_ds - drrg1_ds * dzzg1_dr

      gradps2 = dpsg1_dr**2 * (drrg1_ds**2 + dzzg1_ds**2) / rzjac**2

      zjdchi = rrg1 * rzjac / dabs(dpsg1_dr)

      psi_n = - (psg1 - ps_axis) / ps_axis

      bphi  = fdia(psi_n) / rrg1

      b02 = bphi**2 + gradps2/rrg1**2

      sum1  = sum1 - ws * zjdchi
      sum2  = sum2 - ws * zjdchi * b02
      sumr  = sumr - ws * zjdchi * rrg1

      if (b02 .gt. b02max) b02max = b02

    enddo

  enddo

  b02av(i+1) = sum2 / sum1
  b0max(i+1) = dsqrt(dabs(b02max))
  rav(i+1)   = sumr / sum1

enddo


!----------------------------- calculate average term in integral
nk=101
allocate(sumk(nk))

r = +1.e0_r8

write(*,*) 'i+1 sqrt(psi_n) B02av B0max fraction_circ'
do i= 1, nr_flux-1

  sum1 = 0.e0_r8

  sumk = 0.e0_r8

  do j=1, np_flux

    n1 = (i-1)*np_flux + j
    n2 = (i-1)*np_flux + j + 1
    n3 = (i  )*np_flux + j + 1
    n4 = (i  )*np_flux + j

    if (j .eq. np_flux) then
      n1 = (i  )*np_flux
      n2 = (i  )*np_flux - np_flux + 1
      n3 = (i  )*np_flux + 1
      n4 = (i  )*np_flux + np_flux
    endif

!------------------------------------- 4 POINT GAUSSIAN INT. IN S -----
    do ngs=1,4

      s  = xgauss(ngs)
      ws = wgauss(ngs)

      call interp(rrflux(1,n1),rrflux(1,n2),rrflux(1,n3),rrflux(1,n4),r,s,rrg1,drrg1_dr,drrg1_ds,drrg1_drs,drrg1_drr,drrg1_dss)
      call interp(zzflux(1,n1),zzflux(1,n2),zzflux(1,n3),zzflux(1,n4),r,s,zzg1,dzzg1_dr,dzzg1_ds,dzzg1_drs,dzzg1_drr,dzzg1_dss)
      call interp(psflux(1,n1),psflux(1,n2),psflux(1,n3),psflux(1,n4),r,s,psg1,dpsg1_dr,dpsg1_ds,dpsg1_drs,dpsg1_drr,dpsg1_dss)

      rzjac =  drrg1_dr * dzzg1_ds - drrg1_ds * dzzg1_dr

      gradps2 = dpsg1_dr**2 * (drrg1_ds**2 + dzzg1_ds**2) / rzjac**2

      zjdchi = rrg1 * rzjac / dabs(dpsg1_dr)

      psi_n = - (psg1 - ps_axis) / ps_axis

      bphi  = fdia(psi_n) / rrg1

      b02 = bphi**2 + gradps2/rrg1**2

      b0 = dsqrt(dabs(b02))
      bm = b0max(i+1)

      do k = 1, nk

        zlam = dble(k-1)/dble(nk-1)

        sumk(k) = sumk(k) - ws*zjdchi*dsqrt(dabs(1.-zlam*b0/bm))

      enddo

      sum1  = sum1 - ws * zjdchi
    enddo

  enddo

!------------------------------------------ integrate over lambda
  sum = 0.e0_r8

  do k=2,nk-1

    zlam = float(k-1)/float(nk-1)
    sum  = sum + zlam * sum1 / sumk(k)

  enddo

  fraction_circ(i+1) = (sum + 0.5e0_r8 * sum1/sumk(nk)) / float(nk-1)
  fraction_circ(i+1) = 0.75e0_r8 * fraction_circ(i+1) * b02av(i+1) / b0max(i+1)**2


  psi_n = - (psg1 - ps_axis) / ps_axis
  write(*,'(i5,4f12.6)') i+1,sqrt(psi_n),b02av(i+1),b0max(i+1),fraction_circ(i+1)

enddo

return
end subroutine helena_circulating



subroutine helena_volume_integrals(nr_flux,np_flux,RR_flux,ZZ_flux,PS_flux, &
                                   powers,volume_integrals,n_var,n_int,n_psi)
!------------------------------------------------------------------------
! subroutines calculates some quantities which are volume integrals
! within a fluxsurface.
!
! requires a mesh aligned on fluxsurface (helena_remesh)
!
!        1)   R               2)   B_phi
!------------------------------------------------------------------------
use helena_gauss_4
use helena_axis
use helena_constants
use itm_types
implicit none

real (r8)    :: rr_flux(4,*),zz_flux(4,*), ps_flux(4,*)
real (r8)    :: powers(n_var,n_int),volume_integrals(n_psi,n_int)
integer   :: nr_flux, np_flux, n_var, n_int, n_psi
!real (r8)           :: fdia
real (r8) r,wr,s,ws,rzjac,psi_n,bphi,  &
     rrg1,drrg1_dr,drrg1_ds,zzg1,dzzg1_dr,dzzg1_ds,psg1,dpsg1_dr,dpsg1_ds
integer i,j,n1,n2,n3,n4,ngr,ngs,m


write(*,*) ' helena volume integrals : ',nr_flux,np_flux

volume_integrals(1:n_psi,1:n_int) = 0.e0_r8

do i=1,nr_flux-1

  volume_integrals(i+1,1:n_int) = volume_integrals(i,1:n_int)
  
  do j=1,np_flux

    n1 = (i-1)*np_flux + j
    n2 = (i-1)*np_flux + j + 1
    n3 = (i  )*np_flux + j + 1
    n4 = (i  )*np_flux + j

    if (j .eq. np_flux) then
      n1 = (i  )*np_flux
      n2 = (i  )*np_flux - np_flux + 1
      n3 = (i  )*np_flux + 1
      n4 = (i  )*np_flux + np_flux
    endif

    do ngr=1,4

      r  = xgauss(ngr)
      wr = wgauss(ngr)

      do ngs=1,4

        s  = xgauss(ngs)
        ws = wgauss(ngs)

        CALL interp2(rr_flux(1,n1),rr_flux(1,n2),rr_flux(1,n3),rr_flux(1,n4),r,s,rrg1,drrg1_dr,drrg1_ds)
        CALL interp2(zz_flux(1,n1),zz_flux(1,n2),zz_flux(1,n3),zz_flux(1,n4),r,s,zzg1,dzzg1_dr,dzzg1_ds)
        CALL interp2(ps_flux(1,n1),ps_flux(1,n2),ps_flux(1,n3),ps_flux(1,n4),r,s,psg1,dpsg1_dr,dpsg1_ds)

        rzjac =   drrg1_dr * dzzg1_ds - drrg1_ds * dzzg1_dr
        psi_n = - (psg1 - ps_axis) / ps_axis

        bphi  = fdia(psi_n) / rrg1

        do m=1,n_int

         volume_integrals(i+1,m) = volume_integrals(i+1,m)   &
                                 + rrg1**powers(1,m) *  bphi**powers(2,m) * wr * ws * 2.e0_r8 *pi * rrg1 * rzjac     ! 3D volume integral (not area)

        enddo

      enddo

    enddo

  enddo

enddo

return
end subroutine helena_volume_integrals

subroutine helena_remesh(nrnew,npnew,RRnew,ZZnew,PSInew)
!------------------------------------------------------------------------
! subroutine calculates a new flux surface grid
!------------------------------------------------------------------------
use helena_surfaces
use helena_elements
use helena_grid
use helena_axis
use itm_types
implicit none

real (r8)             :: rrnew(4,*),zznew(4,*),psinew(4,*), abltg(3), xtmp, tol_t
real (r8),allocatable :: s_values(:),tht_start(:),tht_end(:)
real (r8),allocatable :: sp1(:),sp2(:),sp3(:),sp4(:)
real (r8) pi,dpsi_ds,tht_min,tht_max,rr1,rr2,ss1,ss2,  &
     rrg1,drrg1_dr,drrg1_ds,zzg1,dzzg1_dr,dzzg1_ds,rrg2,drrg2_dr,  &
     drrg2_ds,zzg2,dzzg2_dr,dzzg2_ds,tht1,tht2,tht_tmp,theta,  &
     drr1,dss1,drr2,dss2,drrg1_dt,dzzg1_dt,drrg2_dt,dzzg2_dt,  &
     rz1,rz2,drz1,drz2,rz0,a3,a2,a1,a0,t,t2,t3,ri,dri,si,  &
     dsi,drrg1_drs,drrg1_drr,drrg1_dss,dzzg1_drs,dzzg1_drr,  &
     dzzg1_dss,psg1,dpsg1_dr,dpsg1_ds,dpsg1_drs,dpsg1_drr,dpsg1_dss,  &
     check,rad2,th_z,th_r,th_rr,th_zz,th_rz,dth_ds,dth_dr,dth_drr,dth_drs,dth_dss,  &
     rzjac,ps_r,ps_z,ejac,ptjac,rt,st,dptjac_dr,dptjac_ds,rpt,spt,  &
     dr_dpt,dz_dpt,dr_dz,dr_dr,ds_dz,ds_dr,dps_drr,dps_dzz,cx,cy,tn,tn2,cn
integer nrnew,npnew,i,k,i_elm,node1,node2,node3,node4,j,ifail,inode

if (nrnew .ne. n_psi+1) write(*,*) ' MAJOR PROBLEM in REMESH! nrnew <> n_psi+1'

write(*,*) '**************************************'
write(*,*) '*    helena remesh : ',nrnew,npnew,'*'
write(*,*) '**************************************'

tol_t = 1.d-8             ! tolerance

allocate(s_values(n_psi+1),tht_start(n_pieces_max),tht_end(n_pieces_max))

rrnew(1:4,1:nrnew*npnew)  = 0.e0_r8
zznew(1:4,1:nrnew*npnew)  = 0.e0_r8
psinew(1:4,1:nrnew*npnew) = 0.e0_r8

pi = 2.e0_r8*asin(1.e0_r8)

do i=1,n_psi+1
  s_values(i) = float(i-1)/float(n_psi)
enddo

allocate(sp1(n_psi+1),sp2(n_psi+1),sp3(n_psi+1),sp4(n_psi+1))

!call spline(n_psi+1,s_values,psi_values,0.,0.,1,sp1,sp2,sp3,sp4)

do i=1,n_psi

!    xtmp = spwert(n_psi,s_values(i),sp1,sp2,sp3,sp4,s_values,abltg)
!    dpsi_ds = abltg(1)

    dpsi_ds = - ps_axis * 2.e0_r8*s_values(i+1)   ! for an equidistant grid in s

    tht_min =  1d20
    tht_max = -1d20

    do k=1, flux_surfaces(i)%n_pieces

      rr1  = flux_surfaces(i)%r(1,k)
      rr2  = flux_surfaces(i)%r(3,k)

      ss1  = flux_surfaces(i)%s(1,k)
      ss2  = flux_surfaces(i)%s(3,k)

      i_elm = flux_surfaces(i)%elm(k)

      node1 = elm_list(1,i_elm)
      node2 = elm_list(2,i_elm)
      node3 = elm_list(3,i_elm)
      node4 = elm_list(4,i_elm)

      call interp2(rr(:,node1),rr(:,node2),rr(:,node3),rr(:,node4),rr1,ss1,rrg1,drrg1_dr,drrg1_ds)
      call interp2(zz(:,node1),zz(:,node2),zz(:,node3),zz(:,node4),rr1,ss1,zzg1,dzzg1_dr,dzzg1_ds)
      call interp2(rr(:,node1),rr(:,node2),rr(:,node3),rr(:,node4),rr2,ss2,rrg2,drrg2_dr,drrg2_ds)
      call interp2(zz(:,node1),zz(:,node2),zz(:,node3),zz(:,node4),rr2,ss2,zzg2,dzzg2_dr,dzzg2_ds)

      tht1 = atan2(zzg1-z_axis,rrg1-r_axis)
      tht2 = atan2(zzg2-z_axis,rrg2-r_axis)

      if (tht1 .lt. 0.e0_r8) tht1 = tht1 + 2.e0_r8*pi
      if (tht2 .lt. 0.e0_r8) tht2 = tht2 + 2.e0_r8*pi

      tht_start(k) = min(tht1,tht2)
      tht_end(k)   = max(tht1,tht2)

      if ((tht_end(k) - tht_start(k)) .gt. 3.e0_r8*pi/4.e0_r8) then
         tht_tmp      = tht_end(k)
         tht_end(k)   = tht_start(k)
         tht_start(k) = tht_tmp - 2.e0_r8*pi
      endif

      tht_min = min(tht_min,tht_start(k))
      tht_max = max(tht_max,tht_end(k))

    enddo

    do j=1, npnew

      theta = 2.e0_r8 * pi * float(j-1)/float(npnew)

      if (theta .gt. tht_max) theta = theta - 2.e0_r8*pi
      if (theta .lt. tht_min) theta = theta + 2.e0_r8*pi

      do k=1, flux_surfaces(i)%n_pieces

        if ( (theta .ge. tht_start(k)) .and. (theta .le. tht_end(k)) ) then

          rr1  = flux_surfaces(i)%r(1,k);   ss1  = flux_surfaces(i)%s(1,k)
          drr1 = flux_surfaces(i)%r(2,k);   dss1 = flux_surfaces(i)%s(2,k)
          rr2  = flux_surfaces(i)%r(3,k);   ss2  = flux_surfaces(i)%s(3,k)
          drr2 = flux_surfaces(i)%r(4,k);   dss2 = flux_surfaces(i)%s(4,k)

          i_elm = flux_surfaces(i)%elm(k)

          node1 = elm_list(1,i_elm)
          node2 = elm_list(2,i_elm)
          node3 = elm_list(3,i_elm)
          node4 = elm_list(4,i_elm)

          call interp2(rr(:,node1),rr(:,node2),rr(:,node3),rr(:,node4),rr1,ss1,rrg1,drrg1_dr,drrg1_ds)
          call interp2(zz(:,node1),zz(:,node2),zz(:,node3),zz(:,node4),rr1,ss1,zzg1,dzzg1_dr,dzzg1_ds)
          call interp2(rr(:,node1),rr(:,node2),rr(:,node3),rr(:,node4),rr2,ss2,rrg2,drrg2_dr,drrg2_ds)
          call interp2(zz(:,node1),zz(:,node2),zz(:,node3),zz(:,node4),rr2,ss2,zzg2,dzzg2_dr,dzzg2_ds)

          drrg1_dt = drrg1_dr * drr1 + drrg1_ds * dss1
          dzzg1_dt = dzzg1_dr * drr1 + dzzg1_ds * dss1
          drrg2_dt = drrg2_dr * drr2 + drrg2_ds * dss2
          dzzg2_dt = dzzg2_dr * drr2 + dzzg2_ds * dss2

          rz1  = rrg1     * tan(theta) - zzg1
          rz2  = rrg2     * tan(theta) - zzg2
          drz1 = drrg1_dt * tan(theta) - dzzg1_dt
          drz2 = drrg2_dt * tan(theta) - dzzg2_dt

          rz0  = r_axis  * tan(theta) - z_axis

          a3 = (   rz1 + drz1 -   rz2 + drz2 )/4.e0_r8
          a2 = (       - drz1         + drz2 )/4.e0_r8
          a1 = (-3.e0_r8*rz1 - drz1 + 3.e0_r8*rz2 - drz2 )/4.e0_r8
          a0 = ( 2.e0_r8*rz1 + drz1 + 2.e0_r8*rz2 - drz2 )/4.e0_r8 - rz0

          call solvp3(a0,a1,a2,a3,t,t2,t3,ifail)

          if (abs(t) .le. 1.e0_r8+tol_t) then

            call cub1d(rr1, drr1, rr2, drr2, t, ri, dri)
            call cub1d(ss1, dss1, ss2, dss2, t, si, dsi)

            call interp(rr(:,node1),rr(:,node2),rr(:,node3),rr(:,node4),ri,si, &
                        rrg1,drrg1_dr,drrg1_ds,drrg1_drs,drrg1_drr,drrg1_dss)
            call interp(zz(:,node1),zz(:,node2),zz(:,node3),zz(:,node4),ri,si, &
                        zzg1,dzzg1_dr,dzzg1_ds,dzzg1_drs,dzzg1_drr,dzzg1_dss)
            call interp(psi(:,node1),psi(:,node2),psi(:,node3),psi(:,node4),ri,si, &
                        psg1,dpsg1_dr,dpsg1_ds,dpsg1_drs,dpsg1_drr,dpsg1_dss)

            check = atan2(zzg1- z_axis,rrg1-r_axis)
            if (check .lt. 0.e0_r8) check = check + 2.e0_r8*pi

            rad2    =   (rrg1 - r_axis)**2 + (zzg1 - z_axis)**2
            th_z    =   (rrg1 - r_axis) / rad2
            th_r    = - (zzg1 - z_axis) / rad2

            th_rr   = 2*(zzg1 - z_axis) * (rrg1 - r_axis) / rad2**2
            th_zz   = - th_rr
            th_rz   = ( (zzg1 - z_axis)**2 - (rrg1 - r_axis)**2 ) / rad2**2

            dth_ds  = th_r * drrg1_ds + th_z * dzzg1_ds
            dth_dr  = th_r * drrg1_dr + th_z * dzzg1_dr

            dth_drr = th_rr * drrg1_dr * drrg1_dr + th_r * drrg1_drr + 2.* th_rz * drrg1_dr * dzzg1_dr &
                    + th_zz * dzzg1_dr * dzzg1_dr + th_z * dzzg1_drr

            dth_drs = th_rr * drrg1_dr * drrg1_ds + th_rz * drrg1_dr * dzzg1_ds + th_r * drrg1_drs &
                    + th_rz * dzzg1_dr * drrg1_ds + th_zz * dzzg1_dr * dzzg1_ds + th_z * dzzg1_drs

            dth_dss = th_rr * drrg1_ds * drrg1_ds + th_r * drrg1_dss + 2.* th_rz * drrg1_ds * dzzg1_ds &
                    + th_zz * dzzg1_ds * dzzg1_ds + th_z * dzzg1_dss

            rzjac   = drrg1_dr * dzzg1_ds - drrg1_ds * dzzg1_dr

            ps_r    = (  dzzg1_ds * dpsg1_dr - dzzg1_dr * dpsg1_ds) / rzjac
            ps_z    = (- drrg1_ds * dpsg1_dr + drrg1_dr * dpsg1_ds) / rzjac

            ejac   =  (ps_r * th_z - ps_z * th_r)
            ptjac   = (dpsg1_dr * dth_ds - dpsg1_ds * dth_dr)

            rt      = - dpsg1_ds / ptjac
            st      =   dpsg1_dr / ptjac

            dptjac_dr = dpsg1_drr * dth_ds + dpsg1_dr * dth_drs - dpsg1_drs * dth_dr - dpsg1_ds * dth_drr

            dptjac_ds = dpsg1_drs * dth_ds + dpsg1_dr * dth_dss - dpsg1_dss * dth_dr - dpsg1_ds * dth_drs

            rpt = (-dptjac_dr * dth_ds / ptjac**2 + dth_drs/ptjac) * rt + (- dptjac_ds * dth_ds / ptjac**2 + dth_dss/ptjac) * st

            spt = ( dptjac_dr * dth_dr / ptjac**2 - dth_drr/ ptjac) * rt + (  dptjac_ds * dth_dr / ptjac**2 - dth_drs/ptjac) * st

            dr_dpt = - drrg1_drr * dth_ds * dpsg1_ds / ptjac**2 + drrg1_drs * (dth_ds * dpsg1_dr + dth_dr * dpsg1_ds) / ptjac**2 &
                     + drrg1_dr  * rpt    + drrg1_ds * spt - drrg1_dss * dth_dr * dpsg1_dr / ptjac**2

            dz_dpt = - dzzg1_drr * dth_ds * dpsg1_ds / ptjac**2 + dzzg1_drs * (dth_ds * dpsg1_dr + dth_dr * dpsg1_ds) / ptjac**2 &
                     + dzzg1_dr  * rpt    + dzzg1_ds * spt - dzzg1_dss * dth_dr * dpsg1_dr /ptjac**2


            rrnew(1,npnew*(i) + j)  = rrg1
            zznew(1,npnew*(i) + j)  = zzg1

            rrnew(2,npnew*(i) + j) = (  th_z / ejac) /  (2.e0_r8*(nrnew-1)) * dpsi_ds
            zznew(2,npnew*(i) + j) = (- th_r / ejac) /  (2.e0_r8*(nrnew-1)) * dpsi_ds

            rrnew(3,npnew*(i) + j) = + (- ps_z / ejac) /  (npnew/pi)
            zznew(3,npnew*(i) + j) = + (  ps_r / ejac) /  (npnew/pi)

            rrnew(4,npnew*(i) + j) =  dr_dpt /(2.e0_r8*(nrnew-1)*(npnew)/pi)  * dpsi_ds
            zznew(4,npnew*(i) + j) =  dz_dpt /(2.e0_r8*(nrnew-1)*(npnew)/pi)  * dpsi_ds

            psinew(1,npnew*(i) + j) = psi_values(i)
            psinew(2,npnew*(i) + j) = dpsi_ds / (2.e0_r8*(nrnew-1))
            psinew(3,npnew*(i) + j) = 0.e0_r8
            psinew(4,npnew*(i) + j) = 0.e0_r8

          else

            write(*,*) ' T TOO BIG : ',t,t2,t3

          endif

        endif

      enddo

    enddo

enddo

!----------------------------------- magnetic axis

node1 = elm_list(1,ielm_ax)
node2 = elm_list(2,ielm_ax)
node3 = elm_list(3,ielm_ax)
node4 = elm_list(4,ielm_ax)

call interp(rr(:,node1),rr(:,node2),rr(:,node3),rr(:,node4),r_ax,s_ax,rrg1,drrg1_dr,drrg1_ds,drrg1_drs,drrg1_drr,drrg1_dss)
call interp(zz(:,node1),zz(:,node2),zz(:,node3),zz(:,node4),r_ax,s_ax,zzg1,dzzg1_dr,dzzg1_ds,dzzg1_drs,dzzg1_drr,dzzg1_dss)
call interp(psi(:,node1),psi(:,node2),psi(:,node3),psi(:,node4),r_ax,s_ax,psg1,dpsg1_dr,dpsg1_ds,dpsg1_drs,dpsg1_drr,dpsg1_dss)

ejac  = drrg1_dr * dzzg1_ds - drrg1_ds * dzzg1_dr
dr_dz = - drrg1_ds / ejac
dr_dr = + dzzg1_ds / ejac
ds_dz = + drrg1_dr / ejac
ds_dr = - dzzg1_dr / ejac

dps_drr = dpsg1_drr * dr_dr * dr_dr + 2.*dpsg1_drs * dr_dr * ds_dr + dpsg1_dss * ds_dr * ds_dr
dps_dzz = dpsg1_drr * dr_dz * dr_dz + 2.*dpsg1_drs * dr_dz * ds_dz + dpsg1_dss * ds_dz * ds_dz

crr_axis = dps_drr / 2.e0_r8
czz_axis = dps_dzz / 2.e0_r8

!write(*,*) ' CRR, CZZ : ',CRR_axis,CZZ_axis

cx = crr_axis
cy = czz_axis

do j=1,npnew

  inode = j
  rrnew(1,inode) = rrg1
  zznew(1,inode) = zzg1

  theta = float(j-1)/float(npnew) * 2.e0_r8*pi

  tn  = tan(theta)
  tn2 = tn**2
  cn  = cos(theta)

  if (theta .eq. pi/2.e0_r8) then
    rrnew(2,inode) = 0.e0_r8
    zznew(2,inode) = +1.e0_r8/(sqrt(abs(cy))*2.e0_r8*float(nrnew-1))
    rrnew(4,inode) = -1.e0_r8/(sqrt(abs(cy))*2.e0_r8*float(nrnew-1)*float(npnew)/pi)
    zznew(4,inode) = 0.e0_r8
  ELSEIF (theta .eq. (3.e0_r8*pi/2)) THEN
    rrnew(2,inode) = 0.e0_r8
    zznew(2,inode) = +1.e0_r8/(sqrt(abs(cy))*2.*float(nrnew-1))
    rrnew(4,inode) = -1.e0_r8/(sqrt(abs(cy))*2.*float(nrnew-1)*float(npnew)/pi)
    zznew(4,inode) = 0.e0_r8
  ELSE
    rrnew(2,inode) = + sign(1.e0_r8,cn)/(sqrt(abs(cx+cy*tn2))*2.e0_r8*float(nrnew-1))
    zznew(2,inode) = + abs(tn)/(sqrt(abs(cx+cy*tn2))*2.e0_r8*float(nrnew-1))
    rrnew(4,inode) = - (cx+cy*tn2)**(-1.5e0_r8) * cy * abs(tn) / (cn**2 * 2.e0_r8*float(nrnew-1)*float(npnew)/pi)
    zznew(4,inode) = + cx * (cx + cy*tn2)**(-1.5e0_r8) / (cn*abs(cn) * 2.e0_r8*float(nrnew-1)*float(npnew-1)/pi)
  ENDIF
  IF (theta .gt. pi) THEN
    zznew(2,inode) = - zznew(2,inode)
    rrnew(4,inode) = - rrnew(4,inode)
  ENDIF
  rrnew(3,inode) = 0.e0_r8
  zznew(3,inode) = 0.e0_r8
  psinew(1,inode) = ps_axis

enddo

!call plot_grid(RRnew,ZZnew,nrnew,npnew)

return
end subroutine helena_remesh


subroutine export_helena
use helena_boundary
use helena_profiles
use helena_phys
use helena_grid
use helena_axis
use helena_constants
use itm_types
implicit none
integer :: i

write(*,*)  ' export to HELENA namelist'

open(20,file='helena.nml')
write(20,*) ' &SHAPE IAS = 1, ISHAPE = 2,'
write(20,'(A,i5,A,i5)') '    MFM= ',mf,', MHARM = ',mf/2
do i=1,mf/2
  write(20,'(A,i4,A,e14.6,A,i4,A,e14.6)') '    FM(',2*i-1,')=',fr(2*i-1)/amin,' FM(',2*i,')=',fr(2*i)/amin
enddo
write(20,*) ' &END'
write(20,*) ' &PROFILE '
write(20,*) '    IPAI=7, EPI=1.0, FPI=1.0'
write(20,*) '    IGAM=7'
write(20,*) '    ICUR=2, ECUR=1.0, FCUR=1.0'
write(20,*) '    NPTS = ',n_prof
do i=1,n_prof
  write(20,'(A,i4,A,e12.4,A,i4,A,e12.4,A,i4,A,e12.4)') &
        '    DPR(',i,') = ',dp_dpsi(i),', DF2(',i,') = ',fdf_dpsi(i),', ZJZ(',i,') = ',zjz_psi(i)
enddo
write(20,*) ' &END'
write(20,*) ' &PHYS '
write(20,'(A,f10.6)') '   EPS   = ',amin/rgeo
write(20,'(A,f10.6)') '   ALFA  = ',abs(amin**2 * bgeo / ps_axis)
write(20,'(A,f10.6)') '   B     = ',mu_zero * rgeo**2 * dp_dpsi(1)/fdf_dpsi(1)
write(20,'(A,f10.6)') '   BETAP = ',beta_p
write(20,'(A,f10.6)') '   XIAB  = ',mu_zero * abs(current) / (amin * bgeo)
write(20,'(A,f10.6)') '   RVAC  = ',rgeo
write(20,'(A,f10.6)') '   BVAC  = ',bgeo
write(20,*) ' &END'
write(20,*) ' &NUM'
write(20,*) '    NR    = 51, NP    = 33, NRMAP  = 101,   NPMAP = 129, NCHI = 128'
write(20,*) '    NRCUR = 51, NPCUR = 33, ERRCUR = 1.e-5, NITER = 100, NMESH = 100'
write(20,*) ' &END'
write(20,*) ' &PRI  NPR1=1 &END '
write(20,*) ' &PLOT NPL1=1 &END '
write(20,*) ' &BALL        &END '
close(20)
return
end subroutine export_helena

subroutine phys_values
use helena_grid
use helena_elements
use helena_gauss_4
use helena_matrix
use helena_constants
use helena_axis
use helena_boundary
use helena_phys

use itm_types
implicit none
integer :: nodes(4)
real (r8)  :: xx(4,4), xr(4,4), xs(4,4), yr(4,4), ys(4,4), xjac(4,4), ps(4,4), psn(4,4)
real (r8)  :: sum_p_area, sum_p_vol, sum_c_area, sum_vol, sum_area, zjz
integer :: ngr, ngs, kf, kv, i, vertex
!real (r8)           :: dpdpsi, fdfdpsi, pressure

sum_area   = 0.e0_r8
sum_vol    = 0.e0_r8
sum_p_vol  = 0.e0_r8
sum_p_area = 0.e0_r8
sum_c_area = 0.e0_r8

do i=1, n_elms

  nodes(1) = elm_list(1,i)
  nodes(2) = elm_list(2,i)
  nodes(3) = elm_list(3,i)
  nodes(4) = elm_list(4,i)

  xr = 0.e0_r8
  xs = 0.e0_r8
  yr = 0.e0_r8
  ys = 0.e0_r8
  xx = 0.e0_r8
  ps = 0.e0_r8

  do ngr = 1, 4          ! 4 Gaussian points
    do ngs = 1, 4        ! 4 Gaussian points
      do kf = 1, 4       ! 4 basis functions
        do kv = 1, 4     ! 4 vertices

          vertex = nodes(kv)

          xx(ngr,ngs) = xx(ngr,ngs) + rr(kf,vertex)  * h(ngr,ngs,kv,kf)
          ps(ngr,ngs) = ps(ngr,ngs) + psi(kf,vertex) * h(ngr,ngs,kv,kf)

          xr(ngr,ngs) = xr(ngr,ngs) + rr(kf,vertex) * hr(ngr,ngs,kv,kf)
          xs(ngr,ngs) = xs(ngr,ngs) + rr(kf,vertex) * hs(ngr,ngs,kv,kf)
          yr(ngr,ngs) = yr(ngr,ngs) + zz(kf,vertex) * hr(ngr,ngs,kv,kf)
          ys(ngr,ngs) = ys(ngr,ngs) + zz(kf,vertex) * hs(ngr,ngs,kv,kf)

        enddo
      enddo

      xjac(ngr,ngs) = xr(ngr,ngs)*ys(ngr,ngs) - xs(ngr,ngs)*yr(ngr,ngs)

    enddo
  enddo

  psn(:,:) = (ps_axis - ps(:,:)) / ps_axis

  do ngr = 1, 4
    do ngs = 1, 4

       sum_vol    = sum_vol    + wgauss(ngr)*wgauss(ngs) * xjac(ngr,ngs) * xx(ngr,ngs)
       sum_area   = sum_area   + wgauss(ngr)*wgauss(ngs) * xjac(ngr,ngs)

       sum_p_area = sum_p_area + wgauss(ngr)*wgauss(ngs) * xjac(ngr,ngs)               * pressure(psn(ngr,ngs))
       sum_p_vol  = sum_p_vol  + wgauss(ngr)*wgauss(ngs) * xjac(ngr,ngs) * xx(ngr,ngs) * pressure(psn(ngr,ngs))

       zjz = fdfdpsi(psn(ngr,ngs)) / (mu_zero * xx(ngr,ngs)) + dpdpsi(psn(ngr,ngs)) * xx(ngr,ngs)

       sum_c_area  = sum_c_area +  wgauss(ngr)*wgauss(ngs) * xjac(ngr,ngs) * zjz

    enddo
  enddo

enddo

volume  = 2.e0_r8 * pi * sum_vol
area    = sum_area
current = sum_c_area
beta_p  = 8.e0_r8 * pi / mu_zero * sum_p_area / (sum_c_area**2 )
beta_t  = 2.e0_r8 *      mu_zero * sum_p_area / (sum_area * bgeo**2)
beta_n  = 100.e0_r8 * (4.e0_r8*pi/10.e0_r8) * beta_t / (mu_zero * abs(current) /  (amin * bgeo))

write(*,'(A,f16.8,A)') ' total volume    : ',volume,' m^3'
write(*,'(A,f16.8,A)') ' total area      : ',area,' m^2'
write(*,'(A,e12.4,A)') ' total current   : ',current,' A'
write(*,'(A,f16.8)')   ' poloidal beta   : ',beta_p
write(*,'(A,f16.8)')   ' toroidal beta   : ',beta_t
write(*,'(A,f16.8)')   ' normalised beta : ',beta_n
write(*,'(A,3f16.8)')  ' magnetic axis   : ',r_axis,z_axis,(r_axis-rgeo)/amin

return
end subroutine phys_values




subroutine update_fdf(fmix,fdf_error)
use helena_profiles
use helena_axis
use helena_constants
use helena_surfaces
use itm_types
implicit none
real (r8), allocatable :: r_av(:),or_av(:)
!real (r8)              :: current_profile, dpdpsi
real (r8)              :: psi_n, psi_n2, small, fdf_error, fdf_old, fmix
integer           :: n_int, i

write(*,*) ' UPDATE FDF'

n_psi = n_prof
n_int = 5

small = 1.d-8

allocate(psi_values(n_psi),r_av(n_psi),or_av(n_psi))

do i=1,n_prof
  psi_n         = float(i-1)/float(n_prof-1)
  psi_values(i) = ps_axis - (1. - small) * psi_n  * ps_axis
enddo

call find_flux_surfaces
call fluxsurface_current(r_av, or_av)

fdf_old   = 0.e0_r8
fdf_error = 0.e0_r8

do i=1,n_prof

  psi_n  = 1. - psi_values(i)/ps_axis
  psi_n2 = float(i-1)/float(n_prof-1)

  fdf_old = fdf_dpsi(i)

  fdf_dpsi(i) = mu_zero *(- current_profile(psi_n) - r_av(i) * dpdpsi(psi_n) ) / or_av(i)

  fdf_dpsi(i) = (1.-fmix)*fdf_dpsi(i) + fmix * fdf_old

!  write(*,'(i4,8e12.4)') i,R_av,OR_av,fdf_dpsi(i),dp_dpsi(i),dpdpsi(psi_n),psi_n,psi_n2

  fdf_error= fdf_error + abs(fdf_old - fdf_dpsi(i))

enddo

fdf_error = fdf_error / float(n_prof)

!call lplot6(2,1,psi_values,fdf_dpsi,n_prof,'fdf')

deallocate(psi_values,r_av,or_av)

return
end subroutine update_fdf

subroutine gs_solve(n_iter,rhs_only,solve_only,error_iteration,amix,ifail)
use helena_axis
use helena_grid
use itm_types
implicit none
integer :: i, n_iter, ifail
logical :: solve_only, rhs_only
real (r8)  :: error, error_iteration, amix, t_start, t_end

ifail = 999

call cpu_time(t_start)

do i = 1, n_iter

  call findaxis
  call matrix_gs(ps_axis,rhs_only)
  call solve_matrix(amix,solve_only,error)

!  write(*,'(A,i5,f12.8,e16.8)') ' GS_solve : ',i,R_axis,error

  if (error .lt. error_iteration) then
    write(*,'(A,i5,3e16.8)') ' GS_solve : ',i,r_axis,ps_axis,error
    ifail = 0
    exit
  endif

  solve_only = .true.
  rhs_only   = .true.

enddo

call cpu_time(t_end)

if (error .gt. error_iteration) then
  ifail = n_iter
  write(*,'(A,i5,3e16.8)') ' GS_solve failed : ',i,r_axis,ps_axis,error
endif

write(*,'(A,f12.8)') ' GS, CPU time : ',t_end-t_start

return

end subroutine gs_solve


subroutine find_flux_surfaces
use helena_grid
use helena_elements
use helena_surfaces
use helena_constants
use itm_types
implicit none
real (r8)  :: psimin, psimax, a0, a1, a2, a3
real (r8)  :: psi_test, dpsi_dr(4),dpsi_ds(4), drr_dr(4), drr_ds(4), dzz_dr(4), dzz_ds(4)
real (r8)  :: rr_psi(4), zz_psi(4), r_psi(4), s_psi(4), tht(4)
real (r8)  :: s, s2, s3, r_tmp, s_tmp, psr_tmp, pss_tmp, ttmp
integer :: i, j, ifound, iv, im, is, n1, n2, n3
integer :: ifail, itht(4), itmp


write(*,*) ' find_flux_surfaces : ',n_psi

if (.not. allocated(flux_surfaces)) allocate(flux_surfaces(n_psi))

do j=1, n_psi
  flux_surfaces(j)%n_pieces = 0
  flux_surfaces(j)%elm      = 0
  flux_surfaces(j)%r        = 0
  flux_surfaces(j)%s        = 0
enddo

do i=1, n_elms

  call psi_minmax(i,psimin,psimax)

  do j=1, n_psi

    ifound = 0

    if ((psi_values(j) .ge. psimin) .and. (psi_values(j) .le. psimax)) then

      do iv=1, 4

        im = mod(iv,4) + 1
        n1 = elm_list(iv,i)
        n2 = elm_list(im,i)

        is = mod(iv,2) + 2

        if (iv .gt. 2) then
          n3 = n2
          n2 = n1
          n1 = n3
        endif

        a3 = (   psi(1,n1) + psi(is,n1) -   psi(1,n2) + psi(is,n2))/4.
        a2 = (             - psi(is,n1)               + psi(is,n2))/4.
        a1 = (-3*psi(1,n1) - psi(is,n1) + 3*psi(1,n2) - psi(is,n2))/4.
        a0 = ( 2*psi(1,n1) + psi(is,n1) + 2*psi(1,n2) - psi(is,n2))/4. - psi_values(j)

        call solvp3(a0,a1,a2,a3,s,s2,s3,ifail)

        if (abs(s) .le. 1.) then

          ifound = ifound + 1

          call flux_surface_add_point(s,i,iv,ifound,r_psi,s_psi,dpsi_dr,dpsi_ds)

        endif

        if (abs(s2) .le. 1.) then

          ifound = ifound + 1

          call flux_surface_add_point(s2,i,iv,ifound,r_psi,s_psi,dpsi_dr,dpsi_ds)

          if (abs(s3) .le. 1.)   write(*,*) ' another solution : ',s3

        endif

      enddo ! end of 4 edges

      if (ifound .eq. 2) then

        call flux_surface_add_line(i,j,r_psi(1:2),s_psi(1:2),dpsi_dr(1:2),dpsi_ds(1:2))

      elseif (ifound .eq. 4) then

! complicated : 2 line pieces but which point belongs to which line piece?

        tht(1:4) = atan2(s_psi(1:4),r_psi(1:4))

        where (tht .lt. 0.) tht = tht + 2.e0_r8*pi

        itht(1)= 1; itht(2) = 2; itht(3) = 3; itht(4) = 4

        if (tht(2) .lt. tht(1)) then
          itmp = itht(1); itht(1) = itht(2) ; itht(2) = itmp;
        endif
        if (tht(4) .lt. tht(3)) then
          itmp = itht(3); itht(3) = itht(4) ; itht(4) = itmp;
        endif
        if (tht(itht(3)) .lt. tht(itht(2))) then
          itmp = itht(2); itht(2) = itht(3) ; itht(3) = itmp;
        endif
        if (tht(itht(2)) .lt. tht(itht(1))) then
          itmp = itht(1); itht(1) = itht(2) ; itht(2) = itmp;
        endif
        if (tht(itht(4)) .lt. tht(itht(3))) then
          itmp = itht(3); itht(3) = itht(4) ; itht(4) = itmp;
        endif
        if (tht(itht(3)) .lt. tht(itht(2))) then
          itmp = itht(2); itht(2) = itht(3) ; itht(3) = itmp;
        endif

!        write(*,*) i,j
!        write(*,'(4f12.8)') tht
!        write(*,'(4i5)') itht
!        write(*,'(4f12.8)') tht(itht)
!        write(*,'(4f12.8)') r_psi
!        write(*,'(4f12.8)') s_psi
!        write(*,'(4f12.8)') dpsi_dr
!        write(*,'(4f12.8)') dpsi_ds

        call flux_surface_add_line(i,j,r_psi(itht(1:2)),s_psi(itht(1:2)),dpsi_dr(itht(1:2)),dpsi_ds(itht(1:2)))
        call flux_surface_add_line(i,j,r_psi(itht(3:4)),s_psi(itht(3:4)),dpsi_dr(itht(3:4)),dpsi_ds(itht(3:4)))

      endif

    endif

  enddo

enddo

return
end subroutine find_flux_surfaces

subroutine flux_surface_add_line(i,j,r_psi,s_psi,dpsi_dr,dpsi_ds)
use helena_grid
use helena_elements
use helena_surfaces
use itm_types
implicit none
integer :: i, j,node1, node2, node3, node4
real (r8)  :: psi_value, r_psi(*), s_psi(*), dpsi_dr(*), dpsi_ds(*)
real (r8)  :: rr1, rr2, ss1, ss2,sgn, drs, xl1 ,xl2, drr1, drr2, dss1, dss2, t
real (r8)  :: ri, si, dri, dsi, delta_ri, delta_si,psi_test, dummy_r, dummy_s, dl1, dl2

rr1  = r_psi(1)
rr2  = r_psi(2)
ss1  = s_psi(1)
ss2  = s_psi(2)

sgn = 1.e0_r8

if ((rr1 .eq. -1.e0_r8).and. (dpsi_ds(1).lt. 0.e0_r8)) then
  sgn = -1.e0_r8
elseif ((rr1 .eq. +1.e0_r8).and. (dpsi_ds(1).gt. 0.e0_r8)) then
  sgn = -1.e0_r8
endif
if ((ss1 .eq. -1.e0_r8).and. (dpsi_dr(1).gt. 0.e0_r8)) then
  sgn = -1.e0_r8
elseif ((ss1 .eq. +1.e0_r8).and. (dpsi_dr(1).lt. 0.e0_r8)) then
  sgn = -1.e0_r8
endif

drs = sgn * sqrt((rr2-rr1)**2 + (ss2-ss1)**2)

xl1 = 0.5e0_r8 * drs / sqrt(dpsi_dr(1)**2 + dpsi_ds(1)**2)  ! temporary fix; waiting for ideas
xl2 = 0.5e0_r8 * drs / sqrt(dpsi_dr(2)**2 + dpsi_ds(2)**2)

drr1 =   xl1 * dpsi_ds(1)
drr2 =   xl2 * dpsi_ds(2)
dss1 = - xl1 * dpsi_dr(1)
dss2 = - xl2 * dpsi_dr(2)

!------------------------------------------- attempt a correction step to correct drr1,drr2,dss1,dss2

t = 0.e0_r8                 ! make the curve fit exactly at t=0.

call cub1d(rr1, drr1, rr2, drr2, t, ri, dri)
call cub1d(ss1, dss1, ss2, dss2, t, si, dsi)

node1 = elm_list(1,i);  node2 = elm_list(2,i);  node3 = elm_list(3,i);  node4 = elm_list(4,i)

call interp2(psi(:,node1),psi(:,node2),psi(:,node3),psi(:,node4),ri,si,psi_test,dummy_r,dummy_s)

delta_ri = 0.5e0_r8 * (psi_values(j) - psi_test) / dummy_r
delta_si = 0.5e0_r8 * (psi_values(j) - psi_test) / dummy_s

call solvem2(0.25e0_r8*drr1,-0.25e0_r8*drr2,0.25e0_r8*dss1,-0.25e0_r8*dss2,delta_ri,delta_si,dl1,dl2)

if ((abs(dl1) .lt. 0.1e0_r8) .and. (abs(dl2).lt.0.1e0_r8)) then

  drr1 = (1.e0_r8+ dl1) * drr1
  dss1 = (1.e0_r8+ dl1) * dss1
  drr2 = (1.e0_r8+ dl2) * drr2
  dss2 = (1.e0_r8+ dl2) * dss2

endif

flux_surfaces(j)%n_pieces                         = flux_surfaces(j)%n_pieces + 1
flux_surfaces(j)%elm(flux_surfaces(j)%n_pieces)   = i
flux_surfaces(j)%r(1:4,flux_surfaces(j)%n_pieces) = (/ rr1, drr1, rr2, drr2 /)
flux_surfaces(j)%s(1:4,flux_surfaces(j)%n_pieces) = (/ ss1, dss1, ss2, dss2 /)

return
end subroutine flux_surface_add_line

subroutine flux_surface_add_point(s,i,iv,ifound,r_psi,s_psi,dpsi_dr,dpsi_ds)

use helena_elements
use helena_grid
use itm_types
implicit none

integer :: i, iv, ifound, node1, node2, node3, node4
real (r8)  :: s, r_psi(*), s_psi(*), dpsi_dr(*), dpsi_ds(*), psi_test

if (iv .eq. 1) then
  r_psi(ifound) = -1.e0_r8 ; s_psi(ifound) = s
elseif (iv .eq. 2) then
  r_psi(ifound) = s   ; s_psi(ifound) = 1.e0_r8
elseif (iv .eq. 3) then
  r_psi(ifound) = +1.e0_r8 ; s_psi(ifound) = s
elseif (iv .eq. 4) then
  r_psi(ifound) = s   ; s_psi(ifound) = -1.e0_r8
endif

node1 = elm_list(1,i);  node2 = elm_list(2,i);  node3 = elm_list(3,i);  node4 = elm_list(4,i)

call interp2(psi(:,node1),psi(:,node2),psi(:,node3),psi(:,node4), &
             r_psi(ifound),s_psi(ifound),psi_test,dpsi_dr(ifound),dpsi_ds(ifound))

return
end subroutine flux_surface_add_point

subroutine plot_flux_surfaces
use helena_surfaces
use helena_elements
use helena_grid
use itm_types
implicit none
integer          :: j, k,ip, nplot, node1, node2, node3, node4, i_elm
real (r8)           :: t, rr1, rr2, drr1, drr2, ss1, ss2, dss1, dss2, ri, si, dri, dsi, dummy1, dummy2
real (r8),allocatable :: rplot(:), zplot(:)
character*13     :: label

nplot=5
allocate(rplot(nplot),zplot(nplot))

label= 'Flux surfaces'

if (((z_max-z_min)/(r_max-r_min)) .lt. 1.5e0_r8 ) then
  CALL nframe(2,21,1,r_min,r_max,z_min,z_max,label,13,'R [m]',5,'Z [m]',5)
else
  CALL nframe(22,1,1,r_min,r_max,z_min,z_max,label,13,'R [m]',5,'Z [m]',5)
endif
  
!call plot_elements

do j=1,n_psi

  do k=1,flux_surfaces(j)%n_pieces

    i_elm = flux_surfaces(j)%elm(k)

    node1 = elm_list(1,i_elm)
    node2 = elm_list(2,i_elm)
    node3 = elm_list(3,i_elm)
    node4 = elm_list(4,i_elm)

    rr1  = flux_surfaces(j)%r(1,k)
    drr1 = flux_surfaces(j)%r(2,k)
    rr2  = flux_surfaces(j)%r(3,k)
    drr2 = flux_surfaces(j)%r(4,k)

    ss1  = flux_surfaces(j)%s(1,k)
    dss1 = flux_surfaces(j)%s(2,k)
    ss2  = flux_surfaces(j)%s(3,k)
    dss2 = flux_surfaces(j)%s(4,k)

    do ip = 1, nplot

      t = -1. + 2.*float(ip-1)/float(nplot-1)

      call cub1d(rr1, drr1, rr2, drr2, t, ri, dri)
      call cub1d(ss1, dss1, ss2, dss2, t, si, dsi)

      call interp2(rr(:,node1),rr(:,node2),rr(:,node3),rr(:,node4),ri,si,rplot(ip),dummy1,dummy2)
      call interp2(zz(:,node1),zz(:,node2),zz(:,node3),zz(:,node4),ri,si,zplot(ip),dummy1,dummy2)

    enddo

    call lincol(1)
    if (k .gt. 1) then
      if (i_elm .eq. flux_surfaces(j)%elm(k-1)) then
        call lincol(2)
        write(51,*) ' 1 setlinewidth '
      elseif (i_elm .eq. flux_surfaces(j)%elm(k+1)) then
        call lincol(3)
        write(51,*) ' 2 setlinewidth '
      else
        call lincol(1)
        write(51,*) ' .5 setlinewidth '
      endif
    endif

    call lplot6(22,1,rplot,zplot,-nplot,' ')

  enddo

enddo

return
end subroutine plot_flux_surfaces

subroutine fluxsurface_current(R_av, OR_av)
use helena_grid
use helena_elements
use helena_axis
use helena_surfaces
use helena_constants
use helena_gauss_4
use itm_types
implicit none
real (r8)  :: r_av(*), or_av(*), sum_dl, sum_rav, sum_orav
integer :: i_elm, j, k, n1, n2, n3, node1, node2, node3, node4
real (r8)  :: t,rr1, rr2, drr1, drr2, ss1, ss2, dss1, dss2, ri, si, dri, dsi, dl
real (r8)  :: rrgi, drrgi_dr, drrgi_ds, zzgi, dzzgi_dr, dzzgi_ds, drrgi_dt, dzzgi_dt
real (r8)  :: psgi, dpsgi_dr, dpsgi_ds, psi_r, psi_z, rzjac, grad_psi
integer :: ig, ip


do j=1, n_psi

  sum_dl   = 0.e0_r8
  sum_rav  = 0.e0_r8
  sum_orav = 0.e0_r8

  do k=1, flux_surfaces(j)%n_pieces

    do ig = 1, 4

      t = xgauss(ig)

      rr1  = flux_surfaces(j)%r(1,k)
      drr1 = flux_surfaces(j)%r(2,k)
      rr2  = flux_surfaces(j)%r(3,k)
      drr2 = flux_surfaces(j)%r(4,k)

      ss1  = flux_surfaces(j)%s(1,k)
      dss1 = flux_surfaces(j)%s(2,k)
      ss2  = flux_surfaces(j)%s(3,k)
      dss2 = flux_surfaces(j)%s(4,k)

      call cub1d(rr1, drr1, rr2, drr2, t, ri, dri)
      call cub1d(ss1, dss1, ss2, dss2, t, si, dsi)

      i_elm = flux_surfaces(j)%elm(k)

      node1 = elm_list(1,i_elm)
      node2 = elm_list(2,i_elm)
      node3 = elm_list(3,i_elm)
      node4 = elm_list(4,i_elm)

!-----------------------------------------------------------
      call interp2(psi(:,node1),psi(:,node2),psi(:,node3),psi(:,node4),ri,si,psgi,dpsgi_dr,dpsgi_ds)
      call interp2(rr(:,node1),rr(:,node2),rr(:,node3),rr(:,node4),ri,si,rrgi,drrgi_dr,drrgi_ds)
      call interp2(zz(:,node1),zz(:,node2),zz(:,node3),zz(:,node4),ri,si,zzgi,dzzgi_dr,dzzgi_ds)

      drrgi_dt = drrgi_dr * dri + drrgi_ds * dsi
      dzzgi_dt = dzzgi_dr * dri + dzzgi_ds * dsi

      dl = sqrt(drrgi_dt**2 + dzzgi_dt**2)

      rzjac  = drrgi_dr * dzzgi_ds - drrgi_ds * dzzgi_dr

      psi_r = (   dpsgi_dr * dzzgi_ds - dpsgi_ds * dzzgi_dr ) / rzjac
      psi_z = ( - dpsgi_dr * drrgi_ds + dpsgi_ds * drrgi_dr ) / rzjac

      grad_psi = sqrt(psi_r * psi_r + psi_z * psi_z)
!----------------------------------------------------------

      sum_dl   = sum_dl   + wgauss(ig)        * dl * rrgi / grad_psi
      sum_rav  = sum_rav  + wgauss(ig) * rrgi * dl * rrgi / grad_psi
      sum_orav = sum_orav + wgauss(ig) / rrgi * dl * rrgi / grad_psi

    enddo

  enddo

  r_av(j)  = sum_rav / sum_dl
  or_av(j) = sum_orav / sum_dl

enddo

r_av(1)  = r_axis
or_av(1) = 1.e0_r8/ r_axis

return
end subroutine fluxsurface_current

subroutine fluxsurface_integrals(powers,n_int,n_var,results)
!----------------------------------------------------------------------
! subroutine calculats arbitrary flux surfcae quantities as defined by
! the matrix powers. The entries in the matrix define the power of
! the variables in the integrals
!
!  powers(n_int,n_var)
!    n_int : the number of flux surface quantities
!    n_var : the number of variables in the integrand
!
!    nvar = 1: R     2: B^2          3: grad(psi)
!---------------------------------------------------------------------
use helena_grid
use helena_elements
use helena_axis
use helena_surfaces
use helena_constants
use helena_gauss_4
use itm_types
implicit none
integer :: n_int,n_var
real (r8)  :: powers(n_var,n_int), results(n_psi+1,n_int)
integer :: i_elm, j, k, n1, n2, n3, node1, node2, node3, node4
real (r8)  :: t,rr1, rr2, drr1, drr2, ss1, ss2, dss1, dss2, ri, si, dri, dsi, dl
real (r8)  :: rrgi, drrgi_dr, drrgi_ds, zzgi, dzzgi_dr, dzzgi_ds, drrgi_dt, dzzgi_dt
real (r8)  :: psgi, dpsgi_dr, dpsgi_ds, psi_r, psi_z, rzjac, grad_psi, psi_n
real (r8)  :: sum_dl, b_tot2
integer :: i,m, ig, ip
!real (r8)           :: fdia

write(*,*) ' flux_surface_integrals : ',n_psi,n_int,n_var

results(1:n_psi+1,1:n_int)    = 0.e0_r8

do i=1, n_psi

  psi_n = 1.e0_r8 - psi_values(i)/ps_axis
  
  sum_dl = 0.e0_r8
  
  do k=1, flux_surfaces(i)%n_pieces

    do ig = 1, 4

      t = xgauss(ig)

      rr1  = flux_surfaces(i)%r(1,k)
      drr1 = flux_surfaces(i)%r(2,k)
      rr2  = flux_surfaces(i)%r(3,k)
      drr2 = flux_surfaces(i)%r(4,k)

      ss1  = flux_surfaces(i)%s(1,k)
      dss1 = flux_surfaces(i)%s(2,k)
      ss2  = flux_surfaces(i)%s(3,k)
      dss2 = flux_surfaces(i)%s(4,k)

      call cub1d(rr1, drr1, rr2, drr2, t, ri, dri)
      call cub1d(ss1, dss1, ss2, dss2, t, si, dsi)

      i_elm = flux_surfaces(i)%elm(k)

      node1 = elm_list(1,i_elm)
      node2 = elm_list(2,i_elm)
      node3 = elm_list(3,i_elm)
      node4 = elm_list(4,i_elm)

      call interp2(psi(:,node1),psi(:,node2),psi(:,node3),psi(:,node4),ri,si,psgi,dpsgi_dr,dpsgi_ds)
      call interp2(rr(:,node1),rr(:,node2),rr(:,node3),rr(:,node4),ri,si,rrgi,drrgi_dr,drrgi_ds)
      call interp2(zz(:,node1),zz(:,node2),zz(:,node3),zz(:,node4),ri,si,zzgi,dzzgi_dr,dzzgi_ds)

      drrgi_dt = drrgi_dr * dri + drrgi_ds * dsi
      dzzgi_dt = dzzgi_dr * dri + dzzgi_ds * dsi

      dl = sqrt(drrgi_dt**2 + dzzgi_dt**2)

      rzjac  = drrgi_dr * dzzgi_ds - drrgi_ds * dzzgi_dr

      psi_r = (   dpsgi_dr * dzzgi_ds - dpsgi_ds * dzzgi_dr ) / rzjac
      psi_z = ( - dpsgi_dr * drrgi_ds + dpsgi_ds * drrgi_dr ) / rzjac

      grad_psi = sqrt(psi_r * psi_r + psi_z * psi_z)

      b_tot2 =  (fdia(psi_n) / rrgi)**2 + (grad_psi/ rrgi)**2

      sum_dl = sum_dl +  wgauss(ig) * dl

      do m = 1, n_int
        results(i+1,m) = results(i+1,m) + wgauss(ig) * dl   &
                       * rrgi**powers(1,m)                  &
                       * b_tot2**powers(2,m)                &
                       * grad_psi**powers(3,m)

      enddo

    enddo

  enddo
  
!  results(i+1,1:n_int) = results(i+1,1:n_int) / sum_dl

enddo

!----------------------------------- values on axis

do m=1, n_int
 
  results(1,m) =  r_axis**powers(1,m) * b_axis**powers(2,m)
  
  if (powers(3,m) .eq. -1.e0_r8) then
    results(1,m) = results(1,m) * pi / sqrt(crr_axis*czz_axis)
  elseif (powers(3,m) .gt. 0.e0_r8) then
    results(1,m) = 0.e0_r8
  endif

enddo

return
end subroutine fluxsurface_integrals

subroutine helena_flux_surface_integrals(nr_flux,np_flux,RR_flux,ZZ_flux,PS_flux, &
                                         powers,n_int,n_var,results,q, vprime)
!----------------------------------------------------------------------
! subroutine calculats arbitrary flux surfcae quantities as defined by
! the matrix powers. The entries in the matrix define the power of
! the variables in the integrals
!
! requires fluxsurface aligned mesh
!
!  powers(n_int,n_var)
!    n_int : the number of flux surface quantities
!    n_var : the number of variables in the integrand
!
!    nvar = 1: R     2: B^2          3: grad(psi)
!---------------------------------------------------------------------
use helena_grid
use helena_elements
use helena_axis
use helena_surfaces
use helena_constants
use helena_gauss_4
use itm_types
implicit none
real (r8)  :: rr_flux(4,*),zz_flux(4,*),ps_flux(4,*)
integer :: n_int,n_var, nr_flux, np_flux
real (r8)  :: powers(n_var,n_int), results(nr_flux,n_int), q(nr_flux), vprime(nr_flux)
integer :: i, j, k, m, n1, n2, n3, n4, ngs
real (r8)  :: rrg1, drrg1_dr, drrg1_ds, zzg1, dzzg1_dr, dzzg1_ds, drrg1_dt, dzzg1_dt
real (r8)  :: psg1, dpsg1_dr, dpsg1_ds, psi_r, psi_z, rzjac, grad_psi, ws
real (r8)  :: r, s, dl, sum_dl, b_tot2, bphi, grad_psi2, zjdchi, psi_n
!real (r8)           :: fdia

write(*,*) ' flux_surface_integrals : ',n_psi,n_int,n_var
write(*,*) ' nr_flux : ',nr_flux

if (nr_flux .ne. n_psi+1) write(*,*) ' MAJOR PORBLEM in helena_fluxsurface_integral',nr_flux,n_psi


results(1:n_psi+1,1:n_int)    = 0.e0_r8
q      = 0.e0_r8
vprime = 0.e0_r8

r = + 1.e0_r8

do i=1, nr_flux-1

  sum_dl = 0.e0_r8
  
  do j=1, np_flux

    n1 = (i-1)*np_flux + j
    n2 = (i-1)*np_flux + j + 1
    n3 = (i  )*np_flux + j + 1
    n4 = (i  )*np_flux + j

    if (j .eq. np_flux) then
      n1 = (i  )*np_flux
      n2 = (i  )*np_flux - np_flux + 1
      n3 = (i  )*np_flux + 1
      n4 = (i  )*np_flux + np_flux
    endif

!------------------------------------- 4 POINT GAUSSIAN INT. IN S -----
    do ngs = 1, 4

      s  = xgauss(ngs)
      ws = wgauss(ngs)

      call interp2(rr_flux(1,n1),rr_flux(1,n2),rr_flux(1,n3),rr_flux(1,n4),r,s,rrg1,drrg1_dr,drrg1_ds)
      call interp2(zz_flux(1,n1),zz_flux(1,n2),zz_flux(1,n3),zz_flux(1,n4),r,s,zzg1,dzzg1_dr,dzzg1_ds)
      call interp2(ps_flux(1,n1),ps_flux(1,n2),ps_flux(1,n3),ps_flux(1,n4),r,s,psg1,dpsg1_dr,dpsg1_ds)

      rzjac =  drrg1_dr * dzzg1_ds - drrg1_ds * dzzg1_dr
      
      dl  = sqrt(drrg1_ds**2 + dzzg1_ds**2)

      if(rzjac .eq. 0.0_r8) then
         write(*,*) 'Warning: RZjac == 0, grad_psi2 set to 0.0'
         grad_psi2= 0.0_r8   ! DPC
      else
         grad_psi2= dpsg1_dr**2 * (drrg1_ds**2 + dzzg1_ds**2) / rzjac**2
      endif

      if(dpsg1_dr .eq. 0.0_r8) then
         write(*,*) 'Warning: dPSg1_dr == 0, ZJDCHI set to 0.0'
         zjdchi = 0.0_r8   ! DPC
      else
         zjdchi = 2.e0_r8 * pi * rrg1 * rzjac / dabs(dpsg1_dr)
      endif
      
      psi_n = - (psg1 - ps_axis) / ps_axis

      bphi  = fdia(psi_n) / rrg1

      b_tot2 = bphi**2 + grad_psi2/rrg1**2

      sum_dl = sum_dl +  ws * dl
            
      q(i+1)      = q(i+1)      + ws / (rrg1 * sqrt(grad_psi2)) * dl
      vprime(i+1) = vprime(i+1) + ws * zjdchi 
      
!      write(*,'(A,6i5,12f12.8)') ' test : ',i,j,n1,n2,n3,n4,r,s,PSg1,RRg1,RR_flux(1,N1),RR_flux(1,N2),RR_flux(1,N3),RR_flux(1,N4)

      do m = 1, n_int

!        results(i+1,m) = results(i+1,m) + ws * dl         &
!                       * RRg1**powers(1,m)                &
!                       * B_tot2**powers(2,m)              &
!                       * sqrt(grad_psi2)**powers(3,m)

        results(i+1,m) = results(i+1,m) + ws * zjdchi     &
                       * rrg1**powers(1,m)                &
                       * b_tot2**powers(2,m)              &
                       * sqrt(grad_psi2)**powers(3,m)

      enddo

    enddo

  enddo
  
  results(i+1,1:n_int) = results(i+1,1:n_int) / vprime(i+1)
  
enddo

!----------------------------------- values on axis
q(1) = pi / (r_axis * sqrt(crr_axis*czz_axis))

do m=1, n_int
 
  results(1,m) =  r_axis**powers(1,m) * b_axis**powers(2,m)
  
  if (powers(3,m) .eq. -1.e0_r8) then
    results(1,m) = results(1,m) * pi / sqrt(crr_axis*czz_axis)
  elseif (powers(3,m) .gt. 0.e0_r8) then
    results(1,m) = 0.e0_r8
  endif
enddo

return
end subroutine helena_flux_surface_integrals

subroutine solvem2(a,b,c,d,e,f,x,y)
!-----------------------------------------------------------------------
! solves a 2x2 set of equations
!-----------------------------------------------------------------------
use itm_types
implicit none
real (r8) :: a,b,c,d,e,f,x,y, det

det = a*d - b*c

if (det .ne. 0.e0_r8) then
  x = (   d * e - b * f ) / det
  y = ( - c * e + a * f ) / det
else
  x = 0.e0_r8
  y = 0.e0_r8
endif

return
end subroutine solvem2

SUBROUTINE solvp3(C0,C1,C2,C3,X1,X2,X3,IFAIL)
!-----------------------------------------------------------------------
! SOLVES A CUBIC EQUATION WITH A SOLUTION WITH -1.< X < 1
! CN : THE COEFFICIENT OF X**N, X : THE REAL SOLUTION WITH -1.< X < 1.
!-----------------------------------------------------------------------
use itm_types
implicit none
real (r8)  :: c0, c1, c2, c3, x1, x2, x3, dum, tol, angle
real (r8)  :: a0, a1, a2, aa, bb, cc, det, pi, p, q, u, v
integer :: ifail

x1    = 99.e0_r8
x2    = 999.e0_r8
x3    = 9999.e0_r8
tol   = 0.
ifail = 0

!------------------------------------- 2nd order poly for small c3
IF (abs(c3)/(abs(c1)+abs(c2)+abs(c3)) .LT. 1.d-9) THEN
  aa = c2
  bb = c1
  cc = c0
  det = bb**2 - 4*aa*cc
  IF (det.GE.0.) THEN
    x1 = root(aa,bb,cc,det,1.e0_r8)
    IF (abs(x1).GT. 1.e0_r8 + tol) THEN
      x1 = root(aa,bb,cc,det,-1.e0_r8)
    ENDIF
  ELSE
    ifail = 1
  ENDIF

ELSE
!------------------------------------- 3rd order poly solution
  pi = 2*asin(1.e0_r8)
  a0 = c0 / c3
  a1 = c1 / c3
  a2 = c2 / c3
  p = - (a2**2)/3.e0_r8 + a1
  q = 2.e0_r8/27.e0_r8*(a2**3) - a2 * a1/3.e0_r8 + a0
  det = (p/3.e0_r8)**3 + (q/2.e0_r8)**2
  IF (det .GE. 0) THEN
    u  = sign(1.e0_r8,-q/2.e0_r8+sqrt(det))*abs(-q/2.e0_r8 + sqrt(det))**(1.e0_r8/3.e0_r8)
    v  = sign(1.e0_r8,-q/2.e0_r8-sqrt(det))*abs(-q/2.e0_r8 - sqrt(det))**(1.e0_r8/3.e0_r8)
    x1 =  u + v - a2/3.e0_r8
    IF (abs(x1) .GE. (1.e0_r8+tol)) ifail = 1
  ELSE
    p = -p
    angle = sign(1.e0_r8,p)*acos((q/2.e0_r8)/sqrt(abs(p)/3.e0_r8)**3)
    x1 = -2.e0_r8*sqrt(abs(p)/3.e0_r8)*cos(angle/3.e0_r8) - a2/3.e0_r8
    x2 = -2.e0_r8*sqrt(abs(p)/3.e0_r8)*cos(2*pi/3.e0_r8 - angle/3.e0_r8) - a2/3.e0_r8
    x3 = -2.e0_r8*sqrt(abs(p)/3.e0_r8)*cos(2*pi/3.e0_r8 + angle/3.e0_r8) - a2/3.e0_r8
  ENDIF
  IF (abs(x1) .GT. abs(x2)) THEN
    dum = x1
    x1 = x2
    x2 = dum
  ENDIF
  IF (abs(x2) .GT. abs(x3)) THEN
    dum = x2
    x2 = x3
    x3 = dum
  ENDIF
  IF (abs(x1) .GT. abs(x2)) THEN
    dum = x1
    x1 = x2
    x2 = dum
  ENDIF
ENDIF
IF (abs(x1) .GT. (1.e0_r8 + tol)) ifail=1

RETURN
END SUBROUTINE solvp3

FUNCTION root(A,B,C,D,SGN)
!---------------------------------------------------------------------
! THIS FUNCTION GIVES BETTER ROOTS OF QUADRATICS BY AVOIDING
! CANCELLATION OF SMALLER ROOT
!---------------------------------------------------------------------
use itm_types
implicit none
real (r8) :: root,a, b, c, d, sgn

IF (b*sgn .GE. 0.e0_r8) THEN
  root = -2.e0_r8*c/(b+sgn*sqrt(d))
ELSE
  root = (-b + sgn*sqrt(d)) / (2.e0_r8 * a)
ENDIF
RETURN
END FUNCTION root


subroutine psi_minmax(n,psimin,psimax)
use helena_grid
use helena_elements
use itm_types
implicit none
real (r8)  :: psimin, psimax, psma, psmi, psmima, psim, psimr, psip, psipr
real (r8)  :: aa, bb, cc, det, r, dummy
integer :: i, n, im, n1, n2

psimin = 1d10
psimax =-1d10

do i=1, 4

  im = mod(i,4) + 1
  n1 = elm_list(i,n)
  n2 = elm_list(im,n)

  IF (i .eq. 1) THEN
    psim  = psi(1,n1)
    psimr = psi(3,n1)
    psip  = psi(1,n2)
    psipr = psi(3,n2)
  ELSEif (i .eq. 2) then
    psim  = psi(1,n1)
    psimr = psi(2,n1)
    psip  = psi(1,n2)
    psipr = psi(2,n2)
  ELSEIF (i .eq. 3) then
    psim  =   psi(1,n1)
    psimr = - psi(3,n1)
    psip  =   psi(1,n2)
    psipr = - psi(3,n2)
  ELSEIF (i .eq. 4) then
    psim  =   psi(1,n1)
    psimr = - psi(2,n1)
    psip  =   psi(1,n2)
    psipr = - psi(2,n2)
  ENDIF

  psma = max(psim,psip)
  psmi = min(psim,psip)
  aa =  3.e0_r8 * (psim + psimr - psip + psipr ) / 4.e0_r8
  bb =  ( - psimr + psipr ) / 2.e0_r8
  cc =  ( - 3*psim - psimr + 3*psip - psipr) / 4.e0_r8
  det = bb**2 - 4*aa*cc
  IF (det .GE. 0.e0_r8) THEN
    r = (-bb + sqrt(bb**2-4*aa*cc) ) / (2*aa)
    IF (abs(r) .GT. 1.e0_r8) THEN
      r = (-bb - sqrt(bb**2-4*aa*cc) ) / (2*aa)
    ENDIF
    IF (abs(r) .LE. 1.e0_r8) THEN
      CALL cub1d(psim,psimr,psip,psipr,r,psmima,dummy)
      psma = max(psma,psmima)
      psmi = min(psmi,psmima)
    ENDIF
  ENDIF
  psimin = min(psimin,psmi)
  psimax = max(psimax,psma)
ENDDO

RETURN
END subroutine psi_minmax

SUBROUTINE tht_minmax(n,thtmin,thtmax)
use helena_grid
use helena_elements
use helena_axis
use helena_constants
use itm_types
implicit none
real (r8)  :: thtmin, thtmax, theta
integer :: i, node, n

thtmin = 1d10
thtmax =-1d10

do i=1,4
  node = elm_list(i,n)
  theta = atan2(zz(1,node)-z_axis,rr(1,node)-r_axis)
  if (theta .lt. 0.) then
    theta = theta + 2.e0_r8*pi
  endif
  thtmin = min(theta,thtmin)
  thtmax = max(theta,thtmax)
enddo

return
end SUBROUTINE tht_minmax

subroutine initialise_profiles(n_prof_in,iopt_p,iopt_f,psi_in, pprime_in, ffprime_in,pressure_in, fdia_in, current_in)
!-------------------------------------------------------------------------------
! initialising the profiles contained in helena_profiles
!
!  iopt_p = 1  : use pressure and psi profile
!  iopt_p = 2  : use pprime profile
!
!  iopt_f = 1  : use f_dia and psi profile
!  iopt_f = 2  : use ffprime profile
!  iopt_f = 3  : use current profile and ffprime profile
!  iopt_f = 4  : use current profile and reinitialise ffprime
!-------------------------------------------------------------------------------
use helena_profiles
use helena_boundary
use helena_axis
use helena_constants
use itm_types
implicit none
real (r8)             :: psi, dpsi, sump, sumf, f0, f_old, f_tmp
real (r8),allocatable :: psi_norm(:), psi_prof(:), sp1(:) ,sp2(:), sp3(:), sp4(:)
real (r8)             :: pprime_in(*), ffprime_in(*), pressure_in(*), fdia_in(*), current_in(*), psi_in(*)
integer          :: i, n_prof_in, iopt_p, iopt_f
real (r8)             :: abltg(3), p_tmp, si, rbnd_av, orbnd_av, r_av, or_av, ps0, ps1
!real (r8)             :: spwert

write(*,'(A)')      ' **************************************'
write(*,'(A,i7,A)') ' * initialising profiles',n_prof_in,'       *'
write(*,'(A)')      ' **************************************'

ps0     = psi_in(1)
ps1     = psi_in(n_prof_in)
ps_axis = abs(ps1 - ps0)

p_bnd   = pressure_in(n_prof_in)

allocate(psi_norm(n_prof_in))

do i=1,n_prof_in
  psi_norm(i) = (psi_in(i)- ps0)/(ps1 - ps0)
!  write(*,'(i5,8e12.4)') i,psi_in(i),psi_norm(i),pprime_in(i),ffprime_in(i),pressure_in(i),fdia_in(i),current_in(i)
enddo
  
allocate(sp1(n_prof_in),sp2(n_prof_in),sp3(n_prof_in),sp4(n_prof_in))

n_prof = 101
allocate(psi_prof(n_prof), dp_dpsi(n_prof), zjz_psi(n_prof), fdf_dpsi(n_prof))

if (iopt_p .eq. 1) then       ! calculate pprime from pressure profile

  write(*,*) ' HELENA : using pressure profile as input'

  call spline(n_prof_in,psi_norm,pressure_in,0.e0_r8,0.e0_r8,2,sp1,sp2,sp3,sp4)

  do i = 1, n_prof

    psi_prof(i) = float(i-1)/float(n_prof-1)

    p_tmp = spwert(n_prof_in,psi_prof(i),sp1,sp2,sp3,sp4,psi_norm,abltg)

    dp_dpsi(i) = abltg(1) / abs(ps_axis)

  enddo

  dp_dpsi(1) = dp_dpsi(3)
  dp_dpsi(2) = dp_dpsi(3)

else

  write(*,*) ' HELENA : using pprime profile as input'
    
  call spline(n_prof_in,psi_norm,pprime_in,0.e0_r8,0.e0_r8,2,sp1,sp2,sp3,sp4)

  do i = 1, n_prof

    psi_prof(i) = float(i-1)/float(n_prof-1)

    dp_dpsi(i) = spwert(n_prof_in,psi_prof(i),sp1,sp2,sp3,sp4,psi_norm,abltg)

  enddo

endif

if (iopt_f .eq. 1) then       ! calculate ffprime from fdia profile

  write(*,*) ' HELENA : using fdia profile as input'

  call spline(n_prof_in,psi_norm,fdia_in,0.0_r8,0.0_r8,2,sp1,sp2,sp3,sp4)

  do i = 1, n_prof

    psi_prof(i) = float(i-1)/float(n_prof-1)

    f_tmp = spwert(n_prof_in,psi_prof(i),sp1,sp2,sp3,sp4,psi_norm,abltg)

    fdf_dpsi(i) = f_tmp * abltg(1) / abs(ps_axis)

  enddo

  fdf_dpsi(1) = fdf_dpsi(3)
  fdf_dpsi(2) = fdf_dpsi(3)

elseif ((iopt_f .eq. 2) .or. (iopt_f .eq. 3)) then

  write(*,*) ' HELENA : using ffprime profile as input'
    
  call spline(n_prof_in,psi_norm,ffprime_in,0.e0_r8,0.e0_r8,2,sp1,sp2,sp3,sp4)

  do i = 1, n_prof

    psi_prof(i) = float(i-1)/float(n_prof-1)

    fdf_dpsi(i) = spwert(n_prof_in,psi_prof(i),sp1,sp2,sp3,sp4,psi_norm,abltg)

  enddo

elseif ((iopt_f .eq. 3) .or. (iopt_f .eq. 4)) then

!DPC this is a guess!  I'm also not sure why both this and the previous option overlap
  write(*,*) ' HELENA : using current profile as input'

  call spline(n_prof_in,psi_norm,current_in,0.e0_r8,0.e0_r8,2,sp1,sp2,sp3,sp4)

  rbnd_av  = 0.e0_r8
  orbnd_av = 0.e0_r8
  do i=1,n_bnd
    rbnd_av  = rbnd_av  + r_bnd(i)
    orbnd_av = orbnd_av + 1.e0_r8/ r_bnd(i)
  enddo
  rbnd_av  = rbnd_av  / float(n_bnd)
  orbnd_av = orbnd_av / float(n_bnd)

  do i = 1, n_prof

    zjz_psi(i) = spwert(n_prof_in,psi_prof(i),sp1,sp2,sp3,sp4,psi_norm,abltg)

    si = sqrt(psi_prof(i))

    r_av  = rgeo    + si * ( rbnd_av - rgeo)
    or_av = 1./rgeo + si * (orbnd_av - 1./rgeo)

    if (iopt_f .eq. 4) fdf_dpsi(i) = ( mu_zero * ( - zjz_psi(i) - r_av * dp_dpsi(i) )/ or_av )

  enddo

endif

allocate(p_psi(n_prof),f_psi(n_prof))

dpsi = 1./float(n_prof)
sump = 0.
sumf = 0.
p_psi = 0.
f_psi = 0.

p_psi(n_prof) = p_bnd

do i = 1, n_prof - 1
  sump = sump - (dp_dpsi(n_prof-i+1)+dp_dpsi(n_prof-i))*dpsi/2.
  p_psi(n_prof-i) = sump * abs(ps_axis) + p_bnd

  sumf = sumf - (fdf_dpsi(n_prof-i+1)+fdf_dpsi(n_prof-i))*dpsi/2.
  f_psi(n_prof-i) = sumf * abs(ps_axis)

!DPC  write(*,*) i, p_psi(n_prof-i),f_psi(n_prof-i)
enddo

call lplot6(2,2,psi_prof,dp_dpsi,n_prof,'dp/dpsi')
call lplot6(2,3,psi_prof,fdf_dpsi,n_prof,'fdf/dpsi')
call lplot6(3,2,psi_prof,p_psi,n_prof,'pressure')
call lplot6(3,3,psi_prof,f_psi,n_prof,'F**2-F0**2')

return
end subroutine initialise_profiles

subroutine plot_solution
use helena_grid
use helena_elements
use itm_types
implicit none
integer :: i, ngr, ngs, kf, kv, vertex, nodes(4), nc
real (r8)    :: psi_max, psi_min, rmin, rmax, zmin, zmax, zc(3)
real (r8)    :: xx(4,4), yy(4,4), ps(4,4)
character*19 :: label

label= 'SOLUTION'

psi_max = maxval(psi(1,:))
psi_min = minval(psi(1,:))

CALL nframe(22,11,1,r_min,r_max,z_min,z_max,label,19,'R [m]',5,'Z [m]',5)

nc = 3
zc(1) = -max(abs(psi_max),abs(psi_min))
zc(2) = 0.
zc(3) = abs(zc(1))

do i=1,n_elms

  nodes(1) = elm_list(1,i)
  nodes(2) = elm_list(2,i)
  nodes(3) = elm_list(3,i)
  nodes(4) = elm_list(4,i)

  xx = 0.
  yy = 0.
  ps = 0.

  do ngr = 1, 4          ! 4 Gaussian points
    do ngs = 1, 4        ! 4 Gaussian points

      do kf = 1, 4       ! 4 basis functions
        do kv = 1, 4     ! 4 vertices

           vertex = nodes(kv)

           xx(ngr,ngs) = xx(ngr,ngs) + rr(kf,vertex)  * h(ngr,ngs,kv,kf)
           yy(ngr,ngs) = yy(ngr,ngs) + zz(kf,vertex)  * h(ngr,ngs,kv,kf)
           ps(ngr,ngs) = ps(ngr,ngs) + psi(kf,vertex) * h(ngr,ngs,kv,kf)

        enddo
      enddo
    enddo
  enddo

  call cplotm(2,1,-2,xx,yy,4,-4,1,1,ps,4,zc,nc,'Solution',8,'R [m]',5,'Z [m]',5)

enddo
call lincol(0)
return
end subroutine plot_solution


subroutine plot_elements
use helena_grid
use helena_elements
use itm_types
implicit none
integer      :: i, n1, n2, n3, n4
real (r8)         :: rmin, rmax, zmin, zmax
character*19 :: label

label= 'FINITE ELEMENT GRID'

rmin = minval(rr(1,1:nr*np)); rmax = maxval(rr(1,1:nr*np))
zmin = minval(zz(1,1:nr*np)); zmax = maxval(zz(1,1:nr*np))

CALL nframe(22,11,1,rmin,rmax,zmin,zmax,label,19,'R [m]',5,'Z [m]',5)

do i=1, n_elms

    n1 = elm_list(1,i)
    n2 = elm_list(2,i)
    n3 = elm_list(3,i)
    n4 = elm_list(4,i)

    call plotcu(rr(1,n1), rr(3,n1), zz(1,n1), zz(3,n1), rr(1,n2), rr(3,n2), zz(1,n2), zz(3,n2))
    call plotcu(rr(1,n2), rr(2,n2), zz(1,n2), zz(2,n2), rr(1,n3), rr(2,n3), zz(1,n3), zz(2,n3))
    call plotcu(rr(1,n4), rr(3,n4), zz(1,n4), zz(3,n4), rr(1,n3), rr(3,n3), zz(1,n3), zz(3,n3))
    call plotcu(rr(1,n1), rr(2,n1), zz(1,n1), zz(2,n1), rr(1,n4), rr(2,n4), zz(1,n4), zz(2,n4))

enddo

return
end subroutine plot_elements


subroutine solve_matrix(amix,solve_only,error_iteration)
use helena_matrix
use helena_grid
use helena_elements
use itm_types
implicit none
logical :: solve_only
integer :: i, j, istart, info
real (r8)    :: error_iteration, amix

!write(*,*) n_matrix,n_diag,lda

if (.not. solve_only ) then
  call dpbtrf('L',n_matrix,n_diag,aa,lda,info) ! factorisation (lapack)
  if (info .ne. 0) write(*,*) ' dpbtrf : ',info
endif

call dpbtrs('L',n_matrix,n_diag,1,aa,lda,bb,n_matrix,info)  ! solve (lapack)
if (info .ne. 0) write(*,*) ' dpbtrs : ',info

error_iteration = 0

do i=1,n_nodes
  do j=1,4
    error_iteration = error_iteration + abs(psi(j,index_inverse(i)) - bb(4*(i-1)+j))
    psi(j,index_inverse(i)) = amix * psi(j,index_inverse(i)) + (1. - amix) * bb(4*(i-1)+j)
  enddo
enddo

error_iteration = error_iteration / float(4*n_nodes)

return
end subroutine solve_matrix


subroutine matrix_gs(ps_axis,rhs_only)
use helena_grid
use helena_elements
use helena_gauss_4
use helena_matrix
use helena_constants
use itm_types
implicit none

logical :: rhs_only
integer :: nodes(4)
real (r8)    :: xx(4,4), xr(4,4), xs(4,4), yr(4,4), ys(4,4), xjac(4,4), ps(4,4), psn(4,4)
real (r8)    :: vx(4,4,4,4), vy(4,4,4,4), rhs(4,4)
real (r8)    :: ps_axis, sumq, sumk
integer :: ngr, ngs, kf, kv, lf ,lv, i, vertex, nrow, nrow2, ncol, ncol2
!real (r8) fdfdpsi,dpdpsi
integer n_dof,noff,jstart,j

lda      = 4*(np+1) + 9
n_dof    = 4*nr*np
n_matrix = n_dof
n_diag   = 4*(np+1) + 7

if (.not. allocated(aa)) allocate(aa(lda,n_dof))
if (.not. allocated(bb)) allocate(bb(n_dof))

if (.not. rhs_only) aa = 0.
bb = 0.

do i=1, n_elms

  nodes(1) = elm_list(1,i)
  nodes(2) = elm_list(2,i)
  nodes(3) = elm_list(3,i)
  nodes(4) = elm_list(4,i)

  xr = 0.
  xs = 0.
  yr = 0.
  ys = 0.
  xx = 0.
  ps = 0.

  do ngr = 1, 4          ! 4 Gaussian points
    do ngs = 1, 4        ! 4 Gaussian points
      do kf = 1, 4       ! 4 basis functions
        do kv = 1, 4     ! 4 vertices

          vertex = nodes(kv)

          xx(ngr,ngs) = xx(ngr,ngs) + rr(kf,vertex)  * h(ngr,ngs,kv,kf)
          ps(ngr,ngs) = ps(ngr,ngs) + psi(kf,vertex) * h(ngr,ngs,kv,kf)

          xr(ngr,ngs) = xr(ngr,ngs) + rr(kf,vertex) * hr(ngr,ngs,kv,kf)
          xs(ngr,ngs) = xs(ngr,ngs) + rr(kf,vertex) * hs(ngr,ngs,kv,kf)
          yr(ngr,ngs) = yr(ngr,ngs) + zz(kf,vertex) * hr(ngr,ngs,kv,kf)
          ys(ngr,ngs) = ys(ngr,ngs) + zz(kf,vertex) * hs(ngr,ngs,kv,kf)

        enddo
      enddo

      xjac(ngr,ngs) = xr(ngr,ngs)*ys(ngr,ngs) - xs(ngr,ngs)*yr(ngr,ngs)

    enddo
  enddo

  psn(:,:) = (ps_axis - ps(:,:)) / ps_axis

  do ngr = 1, 4
    do ngs = 1, 4

      rhs(ngr,ngs) = - ( fdfdpsi(psn(ngr,ngs)) + mu_zero * xx(ngr,ngs)*xx(ngr,ngs)*dpdpsi(psn(ngr,ngs)) )  &
                   *  wgauss(ngr)*wgauss(ngs) * xjac(ngr,ngs) / xx(ngr,ngs)

    enddo
  enddo

  do kf = 1, 4
    do kv = 1, 4
      do ngr = 1, 4
        do ngs = 1, 4

          vx(ngr,ngs,kf,kv) =  ys(ngr,ngs) * hr(ngr,ngs,kv,kf)  - yr(ngr,ngs) * hs(ngr,ngs,kv,kf)
          vy(ngr,ngs,kf,kv) = -xs(ngr,ngs) * hr(ngr,ngs,kv,kf)  + xr(ngr,ngs) * hs(ngr,ngs,kv,kf)

        enddo
      enddo
    enddo
  enddo

  do kv = 1, 4

    nrow = 4*(index_list(nodes(kv)) - 1)

    do kf = 1, 4

      nrow2 = nrow + kf

      sumq = 0.

      do ngr = 1, 4
        do ngs = 1, 4

          sumq = sumq - rhs(ngr,ngs) * h(ngr,ngs,kv,kf)

        enddo
      enddo

      bb(nrow2) = bb(nrow2) + sumq

      if (.not. rhs_only) then

        do lv = 1, 4

          ncol = 4*(index_list(nodes(lv)) - 1)

          do lf = 1, 4

            ncol2 = ncol + lf

            noff = nrow2 - ncol2

            if (noff .ge. 0) then

              sumk = 0.

              do ngr = 1, 4
                do ngs = 1, 4

                  sumk = sumk + wgauss(ngr)*wgauss(ngs) / (xjac(ngr,ngs) * xx(ngr,ngs)) &
                       * (vx(ngr,ngs,lf,lv) * vx(ngr,ngs,kf,kv) + vy(ngr,ngs,lf,lv) * vy(ngr,ngs,kf,kv))

                enddo
              enddo

              aa(noff+1,ncol2) = aa(noff+1,ncol2) + sumk

            endif

          enddo
        enddo
      endif
    enddo
  enddo

enddo

!------------------------------- boundary conditions
if (.not. rhs_only) then
  jstart = 4.*(nr-1)*np
  do j=1,np
    aa(1,jstart+4*j-1) = 1.e20
    aa(1,jstart+4*j-3) = 1.e20
  enddo
endif
return
end subroutine matrix_gs


subroutine initialise_elements
use helena_grid
use helena_elements
use helena_gauss_4
use itm_types
implicit none
integer :: i,j, ij1, ijn, jn, k, l, ngr, ngs, kv, kf, i0, j0
real (r8)    :: r, s, r0, s0

n_elms  = (nr-1)*np
n_nodes = nr*np
allocate(elm_list(4,n_elms), index_list(n_nodes), index_inverse(n_nodes))

k=0
do i=1,nr-1
  do j=1,np-1
    k = k + 1
    elm_list(1,k) = (i-1)*np + j
    elm_list(2,k) = (i-1)*np + j + 1
    elm_list(3,k) = (i  )*np + j + 1
    elm_list(4,k) = (i  )*np + j
  enddo
  k = k + 1
  elm_list(1,k) = (i  )*np
  elm_list(2,k) = (i  )*np - np + 1
  elm_list(3,k) = (i  )*np + 1
  elm_list(4,k) = (i  )*np + np
enddo

write(*,*) ' number of elements : ',k

do i = 1, nr

  ij1 = (i-1)*np  + 1         ! the first poloidal point
  index_list(ij1)    = ij1
  index_inverse(ij1) = ij1

  do j = 1, np/2               ! the upper part of the poloidal plane
    ij1 = (i-1)*np + j   + 1
    ijn = (i-1)*np + 2*j
    index_list(ij1)    = ijn
    index_inverse(ijn) = ij1
  enddo

  do j = 1, np/2 - 1
    ij1 = i*np - j + 1
    ijn = (i-1)*np + 2*j + 1
    index_list(ij1)    = ijn
    index_inverse(ijn) = ij1
  enddo

enddo

rs(1,1) = -1.; rs(1,2) = -1.; rs(2,1) = -1.; rs(2,2) = +1.
rs(3,1) = +1.; rs(3,2) = +1.; rs(4,1) = +1.; rs(4,2) = -1.
ij(1,1) = 0;   ij(1,2) = 0;   ij(2,1) = 1;   ij(2,2) = 0
ij(3,1) = 0;   ij(3,2) = 1;   ij(4,1) = 1;   ij(4,2) = 1

do ngr = 1, 4        ! 4 Gaussian points
  do ngs = 1, 4      ! 4 Gaussian points

    r = xgauss(ngr)
    s = xgauss(ngs)

    do kv=1,4    ! 4 corners

      r0 = rs(kv,1)
      s0 = rs(kv,2)

      do kf=1,4  ! 4 basis functions

        i0 = ij(kf,1)
        j0 = ij(kf,2)
        call cubich(i0,j0,r0,s0,r,s,h(ngr,ngs,kv,kf),  hr(ngr,ngs,kv,kf), hs(ngr,ngs,kv,kf),&
                                    hrs(ngr,ngs,kv,kf),hrr(ngr,ngs,kv,kf),hss(ngr,ngs,kv,kf))

      enddo

    enddo
  enddo
enddo

return
end subroutine initialise_elements



subroutine initialise_grid
!--------------------------------------------------------------------
!
!--------------------------------------------------------------------
use helena_grid
use helena_boundary
use itm_types
implicit none
integer          :: i, j, m, node
real (r8)             :: theta(np), s(nr)
real (r8)             :: pi, dt, dr, rm, drm, drmt, drmtr, radius, thtj

allocate(rr(4,nr*np),zz(4,nr*np),psi(4,nr*np))

pi = 2.*asin(1.)
dt = 2.*pi/real(np)
dr = 1./real(nr-1)

do j=1,np
  theta(j) = dt * real(j-1)
enddo

do i=1,nr
  s(i)   = float(i-1)/float(nr-1)
enddo

do i=1,nr
  do  j=1,np
    node   = np*(i-1) + j
    thtj   = theta(j)
    radius = s(i)

    rr(1,node) = rgeo + radius * fr(1) * cos(thtj) / 2.
    rr(2,node) =                 fr(1) * cos(thtj) / 2.
    rr(3,node) =      - radius * fr(1) * sin(thtj) / 2.
    rr(4,node) =               - fr(1) * sin(thtj) / 2.
    zz(1,node) = zgeo + radius * fr(1) * sin(thtj) / 2.
    zz(2,node) =                 fr(1) * sin(thtj) / 2.
    zz(3,node) =        radius * fr(1) * cos(thtj) / 2.
    zz(4,node) =                 fr(1) * cos(thtj) / 2.

!---------------------------- KEEP ELLIPTICITY ON AXIS -----------
    do m = 2, mf/2
      if (m .eq. 2) then
        rm   = radius * ( fr(2*m-1) * cos((m-1)*thtj)           + fr(2*m) * sin((m-1)*thtj) )
        drm  =          ( fr(2*m-1) * cos((m-1)*thtj)           + fr(2*m) * sin((m-1)*thtj))
        drmt = radius * (-fr(2*m-1) * (m-1)*sin((m-1)*thtj)     + fr(2*m) * (m-1)*cos((m-1)*thtj))
        drmtr=          (-fr(2*m-1) * (m-1)*sin((m-1)*thtj)     + fr(2*m) * (m-1)*cos((m-1)*thtj))
      else
        rm   =      radius**(m-1) * ( fr(2*m-1) * cos((m-1)*thtj)       + fr(2*m) * sin((m-1)*thtj) )
        drm  =(m-1)*radius**(m-2) * ( fr(2*m-1) * cos((m-1)*thtj)       + fr(2*m) * sin((m-1)*thtj))
        drmt =      radius**(m-1) * (-fr(2*m-1) * (m-1)*sin((m-1)*thtj) + fr(2*m) *(m-1)*cos((m-1)*thtj))
        drmtr=(m-1)*radius**(m-2) * (-fr(2*m-1) * (m-1)*sin((m-1)*thtj) + fr(2*m) *(m-1)*cos((m-1)*thtj))
      endif
      rr(1,node) = rr(1,node) + rm  * cos(thtj)
      zz(1,node) = zz(1,node) + rm  * sin(thtj)
      rr(2,node) = rr(2,node) + drm * cos(thtj)
      zz(2,node) = zz(2,node) + drm * sin(thtj)
      rr(3,node) = rr(3,node) - rm  * sin(thtj) + drmt  * cos(thtj)
      zz(3,node) = zz(3,node) + rm  * cos(thtj) + drmt  * sin(thtj)
      rr(4,node) = rr(4,node) - drm * sin(thtj) + drmtr * cos(thtj)
      zz(4,node) = zz(4,node) + drm * cos(thtj) + drmtr * sin(thtj)
    enddo
    rr(2,node) = rr(2,node) * dr/2.
    rr(3,node) = rr(3,node) * dt/2.
    rr(4,node) = rr(4,node) * dr/2. * dt/2.
    zz(2,node) = zz(2,node) * dr/2.
    zz(3,node) = zz(3,node) * dt/2.
    zz(4,node) = zz(4,node) * dr/2. * dt/2.

    psi(1,node) = radius **2  - 1.
    psi(2,node) = 2.* radius * dr / 2.
    psi(3,node) = 0.
    psi(4,node) = 0.

  enddo
enddo

r_min = minval(rr(1,:))
r_max = maxval(rr(1,:))
z_min = minval(zz(1,:))
z_max = maxval(zz(1,:))

return
end subroutine initialise_grid

subroutine plot_grid(R,Z,nr,np)
!-------------------------------------------------------------------
! THE X,Y GRID IS PLOTTED USING THE ISOPARAMETRIC REPRESENTATION
!-------------------------------------------------------------------
use itm_types
implicit none
real (r8)         :: r(4,*),z(4,*), rmin, rmax, zmin, zmax
character*19 :: label
integer      :: nr, np, nbase, i, j

label= 'FINITE ELEMENT GRID'

rmin = minval(r(1,1:nr*np));  rmax = maxval(r(1,1:nr*np))
zmin = minval(z(1,1:nr*np));  zmax = maxval(z(1,1:nr*np))

CALL nframe(22,1,1,rmin,rmax,zmin,zmax,label,19,'R [m]',5,'Z [m]',5)

do i=1, nr
  do j=1, np-1
    nbase = j + (i-1)*np
    call plotcu(r(1,nbase),  r(3,nbase),  z(1,nbase),  z(3,nbase), &
                r(1,nbase+1),r(3,nbase+1),z(1,nbase+1),z(3,nbase+1))
  enddo
enddo
do i=1, nr-1
  do j=1, np-1
    nbase = j + (i-1)*np
    call plotcu(r(1,nbase),   r(2,nbase),    z(1,nbase),   z(2,nbase), &
                r(1,nbase+np),r(2,nbase+np), z(1,nbase+np),z(2,nbase+np))
  enddo
enddo

return
end subroutine plot_grid


subroutine plotcu(X1,X1S,Y1,Y1S,X2,X2S,Y2,Y2S)
!-----------------------------------------------------------------------
! PLOTS A CUBIC LINE FROM X1,Y1 TO X2,Y2 GIVEN BY THE ARRAYS XI(1..4)
! AND YI(1..4) in an existing frame of PPPLIB
!-----------------------------------------------------------------------
use itm_types
implicit none
integer :: nplot, i
parameter(nplot=21)
real (r8)    :: x1,x1s,y1,y1s,x2,x2s,y2,y2s
real (r8)    :: xp(nplot),yp(nplot), s, dummy

do i=1,nplot
  s = -1. + 2.*float(i-1)/float(nplot-1)
  call cub1d(x1,x1s,x2,x2s,s,xp(i),dummy)
  call cub1d(y1,y1s,y2,y2s,s,yp(i),dummy)
enddo
call lplot6(2,1,xp,yp,-nplot,' ')

return
end subroutine plotcu

SUBROUTINE cubich(I,J,R0,S0,R,S,H,HR,HS,HRS,HRR,HSS)
!------------------------------------------------------------------
! SUBROUTINE TO CALCULATE THE VALUE OF THE CUBIC POLYNOMIALS AND
! THE DERIVATIVES OF THE CORNER MARKED BY R0,S0 AT THE POSITION R,S
!------------------------------------------------------------------
use itm_types
implicit none
integer :: i,j
real (r8)    :: h,hr,hs,hrs,hi,hri,hj,hsj,hrr,hss,hrri,hssj, r0, s0, r, s

IF (i.EQ.0) THEN
  hi = - (r+r0)**2 * (r*r0-2.) / 4.
  hri = - (r+r0)*(r*r0-2.)/2. - r0*(r+r0)**2 / 4.
  hrri = - 1.5 * r * r0
ELSE
  hi = + r0 * (r+r0)**2 * (r*r0 - 1.) / 4.
  hri = + r0*(r+r0)*(r*r0-1.)/2. + (r+r0)**2 /4.
  hrri = 1.5*r + .5*r0
ENDIF
IF (j.EQ.0) THEN
  hj = - (s+s0)**2 * (s*s0-2.) / 4.
  hsj = - (s+s0)*(s*s0-2.)/2. - s0*(s+s0)**2 / 4.
  hssj = - 1.5 * s * s0
ELSE
  hj = + s0 * (s+s0)**2 * (s*s0 - 1.) / 4.
  hsj = + s0*(s+s0)*(s*s0-1.)/2. + (s+s0)**2 / 4.
  hssj = 1.5*s + .5*s0
ENDIF
h   = hi   * hj
hr  = hri  * hj
hs  = hi   * hsj
hrs = hri  * hsj
hrr = hrri * hj
hss = hi   * hssj
RETURN
END SUBROUTINE cubich

subroutine cub1d(X1,X1S,X2,X2S,S,X,XS)
!-----------------------------------------------------------------------
! CUBIC HERMITE INTERPOLATION IN ONE DIMENSION
!-----------------------------------------------------------------------
use itm_types
implicit none
real (r8)   :: x1,x1s,x2,x2s,s,x,xs
real (r8)   :: h0m,h0p,h1m,h1p,h0ms,h0ps,h1ms,h1ps

h0m  =  (s-1.e0_r8)**2 *(s+2.e0_r8) * 0.25e0_r8
h0ms =  (s-1.e0_r8)*(s+2.e0_r8)/2.e0_r8 + (s-1.e0_r8)**2 * 0.25e0_r8
h0p  = -(s+1.e0_r8)**2 *(s-2.e0_r8) * 0.25e0_r8
h0ps = -(s+1.e0_r8)*(s-2.e0_r8)/2.e0_r8 - (s+1.e0_r8)**2 * 0.25e0_r8
h1m  =  (s-1.e0_r8)**2 *(s+1.e0_r8) * 0.25e0_r8
h1ms =  (s-1.e0_r8)*(s+1.e0_r8)/2. + (s-1.e0_r8)**2 * 0.25e0_r8
h1p  =  (s+1.e0_r8)**2 *(s-1.e0_r8) * 0.25e0_r8
h1ps =  (s+1.e0_r8)*(s-1.e0_r8)/2.e0_r8 + (s+1.e0_r8)**2 * 0.25e0_r8

x  = x1*h0m  + x1s*h1m +  x2*h0p  + x2s*h1p
xs = x1*h0ms + x1s*h1ms + x2*h0ps + x2s*h1ps

return
end subroutine cub1d

subroutine fshape
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
use helena_boundary
use helena_constants
use itm_types
implicit none
real (r8)               :: xj, yj, ga
real (r8), allocatable  :: theta(:), gamma(:), xv(:),yv(:)
real (r8)               :: angle, error, gamm
real (r8), allocatable  :: tht_tmp(:),  fr_tmp(:), work(:)
real (r8), allocatable  :: tht_sort(:), fr_sort(:), dfr_sort(:)
integer, allocatable :: index_order(:)
real (r8)               :: tht, rbnd_av, orbnd_av, values(4)
integer              :: m, igrinv, i, j, ishape, ieast(1), iwest(1), n_bnd_short
parameter(error = 1.d-8)

allocate(fr(mf+2))
allocate(theta(mf),gamma(mf),xv(mf),yv(mf))

if (n_bnd .gt. 1) then

  write(*,*) ' fshape : (R,Z) set given on ',n_bnd,' points'

  reast = maxval(r_bnd)
  rwest = minval(r_bnd)
  ieast = maxloc(r_bnd)
  iwest = minloc(r_bnd)
  rgeo = (reast + rwest) /2.
  zgeo = (z_bnd(ieast(1))+z_bnd(iwest(1)))/2.
  amin = (reast - rwest)/2.

  write(*,'(A,3f12.8)') ' Rgeo, Zgeo : ',rgeo,zgeo,amin

  allocate(tht_tmp(n_bnd),fr_tmp(n_bnd),work(3*n_bnd+6))
  allocate(tht_sort(n_bnd+2),fr_sort(n_bnd+2),dfr_sort(n_bnd+2))
  allocate(index_order(n_bnd+2))

  rbnd_av  = 0.
  orbnd_av = 0.

  do i=1,n_bnd

    tht_tmp(i) = atan2(z_bnd(i)-zgeo,r_bnd(i)-rgeo)

!    if (i .gt. 1) then
!      if (tht_tmp(i) .lt. tht_tmp(i-1)) tht_tmp(i) = tht_tmp(i) + 2.*PI
!    endif

    fr_tmp(i)  = sqrt((r_bnd(i)-rgeo)**2 + (z_bnd(i)-zgeo)**2)

!    write(*,'(i5,2f12.8)') i,tht_tmp(i),fr_tmp(i)

    rbnd_av  =  rbnd_av + r_bnd(i)
    orbnd_av = orbnd_av +  1. / r_bnd(i)

  enddo

  if (abs(tht_tmp(n_bnd) - tht_tmp(1)) .lt. 1.e-6)  n_bnd = n_bnd - 1

  rbnd_av  = rbnd_av  / float(n_bnd)
  orbnd_av = orbnd_av / float(n_bnd)

!  write(*,*)  ' Rbnd_av  : ',Rbnd_av
!  write(*,*)  ' ORbnd_av : ',ORbnd_av

  call qsort2(index_order,n_bnd,tht_tmp)

  do i=1,n_bnd
    tht_sort(i) = tht_tmp(index_order(i))
    fr_sort(i)  = fr_tmp(index_order(i))
!    write(*,*) i,tht_sort(i),fr_sort(i)
  enddo

  n_bnd_short = n_bnd
  do i=2,n_bnd

    if ((tht_sort(i) - tht_sort(1)) .gt. 2.*pi) then
      n_bnd_short = i - 1
      exit
    endif

  enddo

  if (abs(tht_sort(n_bnd_short)- tht_sort(1)) .lt. 1.e-6) then
    tht_sort(n_bnd_short) = tht_sort(1) + 2.*pi
    fr_sort(n_bnd_short)  = fr_sort(1)
  else
    tht_sort(n_bnd_short+1) = tht_sort(1) + 2.*pi
    fr_sort(n_bnd_short+1)  = fr_sort(1)
    n_bnd_short = n_bnd_short + 1
  endif

!  write(*,*) ' n_bnd, n_bnd_short : ',n_bnd, n_bnd_short

  call tb15a(n_bnd_short,tht_sort,fr_sort,dfr_sort,work,6)

  do i = 1, mf

   tht = 2.*pi * float(i-1)/float(mf)

   if (tht .lt. tht_sort(1))          tht = tht + 2.*pi
   if (tht .gt. tht_sort(n_bnd_short)) tht = tht - 2.*pi

   call tg02a(0,n_bnd_short,tht_sort,fr_sort,dfr_sort,tht,values)

   fr(i) = values(1)
  enddo

else

 write(*,*) ' fshape : using moments'

!------------------------------------ THETA(GAMMA(J)) ------------
  do j=1,mf

    ga = 2*pi*(j-1.)/REAL(mf)

    if (ga .le. pi) then
      xj=         amin * cos(ga + tria_u*sin(ga) + quad_u*sin(2*ga))
      yj= ellip * amin * sin(ga)
    else
      xj=         amin * cos(ga + tria_l*sin(ga) + quad_l*sin(2*ga))
      yj= ellip * amin * sin(ga)
    endif
    theta(j) = atan2(yj,xj)

  enddo

!-------- INVERSION OF THETA(GAMMA(J)) TO GAMMA(THETA(J)) ------  --
  CALL grid2nv(theta,gamma,mf,error,igrinv)

  do j=1,mf
    gamm=gamma(j)

    if (ga .le. pi) then
      xv(j) = amin * cos(gamm + tria_u*sin(gamm) + quad_u*sin(2*gamm))
      yv(j) = ellip * amin * sin(gamm)
    else
      xv(j) = amin * cos(gamm + tria_l*sin(gamm) + quad_l*sin(2*gamm))
      yv(j) = ellip * amin * sin(gamm)
   endif

  enddo

  DO j=1,mf
    fr(j) = sqrt(xv(j)**2 + yv(j)**2)
  enddo

  deallocate(theta, gamma, xv, yv)

endif

!-------------- FOURIER COEFFICIENTS FRFNUL AND FRF(M) OF FR(J).
call rft2(fr,mf,1)

do m=1,mf
  fr(m) = 2. * fr(m) / float(mf)
enddo
do m=2,mf,2
  fr(m) = - fr(m)
enddo
RETURN
END subroutine fshape

subroutine findaxis
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
use helena_grid
use helena_elements
use helena_boundary
use helena_gauss_4
use helena_axis
use itm_types
implicit none
real (r8)  :: ps_tmp(4,4), psi_min, axis_min
integer :: nodes(4), i, ngr, ngs, kf, kv, vertex, elm_min, ij_axis(2), ifail
real (r8)  :: x(2), r, s, xerr, ferr, rs_tolerance, zpsir, zpsis
real (r8)  :: drrg1_dr,drrg1_ds,drrg1_drs,drrg1_drr,drrg1_dss
real (r8)  :: dzzg1_dr,dzzg1_ds,dzzg1_drs,dzzg1_drr,dzzg1_dss
real (r8)  :: dpsg1_dr,dpsg1_ds,dpsg1_drs,dpsg1_drr,dpsg1_dss
real (r8)  :: ejac, dr_dz, dr_dr, ds_dz, ds_dr, dps_drr, dps_dzz
integer :: node1, node2, node3, node4
!real (r8)          :: fdia
logical :: early_exit
parameter(rs_tolerance = 1.d-8)

early_exit = .false.
!---------------------------------- first see if axis is in the same element as before
if ((ielm_ax .ge. 1) .and. (ielm_ax .le. n_elms)) then
  nodes(1) = elm_list(1,ielm_ax)
  nodes(2) = elm_list(2,ielm_ax)
  nodes(3) = elm_list(3,ielm_ax)
  nodes(4) = elm_list(4,ielm_ax)

  call mnewtax(psi(1:4,nodes(1)),psi(1:4,nodes(2)),psi(1:4,nodes(3)),psi(1:4,nodes(4)),r,s,xerr,ferr,ifail)

  if ((ifail .eq. 0) .and. ((abs(r)-1.) .le. rs_tolerance) .and. ((abs(s)-1.) .le. rs_tolerance)) then

    call interp2(psi(1:4,nodes(1)),psi(1:4,nodes(2)),psi(1:4,nodes(3)),psi(1:4,nodes(4)),r,s,ps_axis,zpsir,zpsis)
    call interp1(rr(1:4,nodes(1)),rr(1:4,nodes(2)),rr(1:4,nodes(3)),rr(1:4,nodes(4)),r,s,r_axis)
    call interp1(zz(1:4,nodes(1)),zz(1:4,nodes(2)),zz(1:4,nodes(3)),zz(1:4,nodes(4)),r,s,z_axis)

    r_ax = r
    s_ax = s

    write(*,*) ' EARLY EXIT !'
    write(*,'(A,3f12.8)') ' magnetic axis : ',r_axis,z_axis,ps_axis
    write(*,'(A,f12.8)') ' delta (EARLY EXIT) : ',(r_axis-rgeo)/amin
    early_exit = .true.

  endif
endif

if (.not. early_exit) then

  psi_min = 1.d20

  do i=1,n_elms/2

    nodes(1) = elm_list(1,i)
    nodes(2) = elm_list(2,i)
    nodes(3) = elm_list(3,i)
    nodes(4) = elm_list(4,i)

    ps_tmp = 0.

    do ngr = 1, 4
      do ngs= 1, 4

        do kf = 1, 4       ! 4 basis functions
          do kv = 1, 4     ! 4 vertices

             vertex = nodes(kv)

             ps_tmp(ngr,ngs) = ps_tmp(ngr,ngs) + psi(kf,vertex) * h(ngr,ngs,kv,kf)

          enddo
        enddo
      enddo
    enddo

    axis_min = minval(ps_tmp)

    if (psi_min .gt. axis_min) then
      psi_min = axis_min
      elm_min = i
      ij_axis = minloc(ps_tmp)
    endif

  enddo

  nodes(1) = elm_list(1,elm_min)
  nodes(2) = elm_list(2,elm_min)
  nodes(3) = elm_list(3,elm_min)
  nodes(4) = elm_list(4,elm_min)

  r=xgauss(ij_axis(1)) ; s=xgauss(ij_axis(2))
  
!  write(*,'(A,e14.6,i6,2e14.6)') 'minimum psi estimate : ',psi_min,elm_min,r,s

  call mnewtax(psi(1:4,nodes(1)),psi(1:4,nodes(2)),psi(1:4,nodes(3)),psi(1:4,nodes(4)),r,s,xerr,ferr,ifail)

  if (ifail .ne. 0 ) write(*,*) ' MNEWTAX : ifail = ',ifail

  r_ax = r
  s_ax = s
  ielm_ax = elm_min

!  write(*,*) ' find_axis : ',r_ax,s_ax,ielm_ax

endif

node1 = elm_list(1,ielm_ax)
node2 = elm_list(2,ielm_ax)
node3 = elm_list(3,ielm_ax)
node4 = elm_list(4,ielm_ax)

call interp(rr(:,nodes(1)),rr(:,nodes(2)),rr(:,nodes(3)),rr(:,nodes(4)),r_ax,s_ax, &
            r_axis,drrg1_dr,drrg1_ds,drrg1_drs,drrg1_drr,drrg1_dss)
call interp(zz(:,nodes(1)),zz(:,nodes(2)),zz(:,nodes(3)),zz(:,nodes(4)),r_ax,s_ax, &
            z_axis,dzzg1_dr,dzzg1_ds,dzzg1_drs,dzzg1_drr,dzzg1_dss)
call interp(psi(:,nodes(1)),psi(:,nodes(2)),psi(:,nodes(3)),psi(:,nodes(4)),r_ax,s_ax, &
            ps_axis,dpsg1_dr,dpsg1_ds,dpsg1_drs,dpsg1_drr,dpsg1_dss)

ejac  = drrg1_dr * dzzg1_ds - drrg1_ds * dzzg1_dr

if (ejac .ne. 0.e0_r8) then
  dr_dz = - drrg1_ds / ejac
  dr_dr = + dzzg1_ds / ejac
  ds_dz = + drrg1_dr / ejac
  ds_dr = - dzzg1_dr / ejac

  dps_drr = dpsg1_drr * dr_dr * dr_dr + 2.e0_r8*dpsg1_drs * dr_dr * ds_dr + dpsg1_dss * ds_dr * ds_dr
  dps_dzz = dpsg1_drr * dr_dz * dr_dz + 2.e0_r8*dpsg1_drs * dr_dz * ds_dz + dpsg1_dss * ds_dz * ds_dz

  crr_axis = dps_drr / 2.
  czz_axis = dps_dzz / 2.

  b_axis = fdia(0.e0_r8) / r_axis
  q_axis = b_axis   /(2.e0_r8*sqrt(crr_axis*czz_axis))
endif

write(*,'(A,4f14.8)') ' magnetic axis : ',r_axis,z_axis,ps_axis,q_axis

RETURN
END subroutine findaxis

SUBROUTINE mnewtax(ps1, ps2, ps3, ps4, r, s, errx, errf, ifail)
!-------------------------------------------------------------------------
! ROUTINE TO SOLVE TWO NONLINEAR EQUATIONS USING NEWTONS METHOD FROM
! NUMERICAL RECIPES.
! LU DECOMPOSITION REPLACED BY EXPLICIT SOLUTION OF 2X2 MATRIX.
!-------------------------------------------------------------------------
use itm_types
implicit none
REAL (R8)  :: ps1(4), ps2(4), ps3(4), ps4(4)
REAL (R8)  :: r, s, x(2), fvec(2),fjac(2,2)
REAL (R8)  :: tolf,tolx, errf, errx
INTEGER :: ntrial, i, k, ifail
REAL (R8)  :: p(2)
real (r8) zpsi,zpsir,zpsis,zpsirs,zpsirr,zpsiss,temp,dis

ntrial = 50
tolx = 1.d-8
tolf = 1.d-16

x(1) = r
x(2) = s

ifail = 999

do k=1,ntrial

  call interp(ps1,ps2,ps3,ps4,x(1),x(2),zpsi,zpsir,zpsis,zpsirs,zpsirr,zpsiss)

!  write(*,'(i3,5e16.8)') k,x(1),x(2),zpsi,zpsir,zpsis

  fvec(1)   = zpsir
  fvec(2)   = zpsis
  fjac(1,1) = zpsirr
  fjac(1,2) = zpsirs
  fjac(2,1) = zpsirs
  fjac(2,2) = zpsiss

  errf=abs(fvec(1))+abs(fvec(2))

  if (errf .le. tolf) then
    r = x(1)
    s = x(2)
!    write(*,'(A,2e16.8,i5)') ' newton (1) : ',errf,errx,k
    ifail = 0
    return
  endif
  
  p = -fvec

  temp = p(1)
  dis  = fjac(2,2)*fjac(1,1)-fjac(1,2)*fjac(2,1)
  if (dis .ne. 0.e0_r8) then
    p(1) = (fjac(2,2)*p(1)-fjac(1,2)*p(2))/dis
    p(2) = (fjac(1,1)*p(2)-fjac(2,1)*temp)/dis
  endif

  errx=abs(p(1)) + abs(p(2))

  p = min(p,+0.5e0_r8)
  p = max(p,-0.5e0_r8)

  x = x + p

  x = max(x,-1.e0_r8)
  x = min(x,+1.e0_r8)
    
  if (errx .le. tolx) then
    r = x(1)
    s = x(2)
!    write(*,'(A,2e16.8,i5)') ' newton (2) : ',errf,errx,k
    ifail = 0
    return
  endif

enddo

return
end SUBROUTINE mnewtax

SUBROUTINE interp(XN1,XN2,XN3,XN4,R,S,X,XR,XS,XRS,XRR,XSS)
!----------------------------------------------------------------------
! SUBROUTINE CALCULATES THE INTERPOLATED VALUE OF THE FUNCTION X GIVEN
! BY XI(1..4) AT THE FOUR NODES USING BI-CUBIC HERMITE ELEMENTS
!----------------------------------------------------------------------
use itm_types
implicit none
REAL (R8) :: xn1(4),xn2(4),xn3(4),xn4(4)
REAL (R8) :: r,s,x,xr,xs,xrs,xrr,xss
REAL (R8) :: hi0m,hri0m,hrri0m,hi1m,hri1m,hrri1m,hj0m,hsj0m,hssj0m,hj1m,hsj1m,hssj1m
REAL (R8) :: hi0p,hri0p,hrri0p,hi1p,hri1p,hrri1p,hj0p,hsj0p,hssj0p,hj1p,hsj1p,hssj1p

hi0m   = - (r-1.e0_r8)**2 * (-r-2.e0_r8) * 0.25e0_r8
hri0m  = - (r-1.e0_r8)*(-r-2.e0_r8)*0.5e0_r8 +(r-1.e0_r8)**2 * 0.25e0_r8
hrri0m = + 1.5e0_r8 * r
hi1m   = - (r-1.e0_r8)**2 * (-r-1.e0_r8) * 0.25e0_r8
hri1m  = - (r-1.e0_r8)*(-r-1.e0_r8)*0.5e0_r8 + (r-1.e0_r8)**2 *0.25e0_r8
hrri1m = + 1.5e0_r8 * r - .5e0_r8

hj0m   = - (s-1.e0_r8)**2 * (-s-2.e0_r8) * 0.25e0_r8
hsj0m  = - (s-1.e0_r8)*(-s-2.e0_r8)*0.5 +(s-1.e0_r8)**2 * 0.25e0_r8
hssj0m = + 1.5e0_r8 * s
hj1m   = - (s-1.e0_r8)**2 * (-s-1.e0_r8) * 0.25e0_r8
hsj1m  = - (s-1.e0_r8)*(-s-1.e0_r8)*0.5e0_r8 + (s-1.e0_r8)**2 * 0.25e0_r8
hssj1m = + 1.5e0_r8 * s - .5e0_r8

hi0p   = - (r+1.e0_r8)**2 * (r-2.e0_r8) * 0.25e0_r8
hri0p  = - (r+1.e0_r8)*(r-2.e0_r8)*0.5e0_r8 - (r+1.e0_r8)**2 * 0.25e0_r8
hrri0p = - 1.5e0_r8 * r
hi1p   = + (r+1.e0_r8)**2 * (r-1.e0_r8) * 0.25e0_r8
hri1p  = + (r+1.e0_r8)*(r-1.e0_r8)*0.5e0_r8 + (r+1.e0_r8)**2 * 0.25e0_r8
hrri1p = + 1.5e0_r8 * r + .5e0_r8

hj0p  = - (s+1.e0_r8)**2 * (s-2.e0_r8) * 0.25e0_r8
hsj0p = - (s+1.e0_r8)*(s-2.e0_r8)*0.5e0_r8 - (s+1.e0_r8)**2 * 0.25e0_r8
hssj0p = - 1.5e0_r8 * s
hj1p  = + (s+1.e0_r8)**2 * (s-1.e0_r8) * 0.25e0_r8
hsj1p = + (s+1.e0_r8)*(s-1.e0_r8)*0.5e0_r8 + (s+1.e0_r8)**2 * 0.25e0_r8
hssj1p = + 1.5e0_r8 * s + .5e0_r8

x = hi0m*hj0m * xn1(1) + hi1m*hj0m * xn1(2) + hi0m*hj1m * xn1(3) + hi1m*hj1m * xn1(4) &
  + hi0m*hj0p * xn2(1) + hi1m*hj0p * xn2(2) + hi0m*hj1p * xn2(3) + hi1m*hj1p * xn2(4) &
  + hi0p*hj0m * xn4(1) + hi1p*hj0m * xn4(2) + hi0p*hj1m * xn4(3) + hi1p*hj1m * xn4(4) &
  + hi0p*hj0p * xn3(1) + hi1p*hj0p * xn3(2) + hi0p*hj1p * xn3(3) + hi1p*hj1p * xn3(4)

xr = hri0m*hj0m * xn1(1) + hri1m*hj0m * xn1(2) + hri0m*hj1m * xn1(3) + hri1m*hj1m * xn1(4) &
   + hri0m*hj0p * xn2(1) + hri1m*hj0p * xn2(2) + hri0m*hj1p * xn2(3) + hri1m*hj1p * xn2(4) &
   + hri0p*hj0m * xn4(1) + hri1p*hj0m * xn4(2) + hri0p*hj1m * xn4(3) + hri1p*hj1m * xn4(4) &
   + hri0p*hj0p * xn3(1) + hri1p*hj0p * xn3(2) + hri0p*hj1p * xn3(3) + hri1p*hj1p * xn3(4)

xs = hi0m*hsj0m * xn1(1) + hi1m*hsj0m * xn1(2) + hi0m*hsj1m * xn1(3) + hi1m*hsj1m * xn1(4) &
   + hi0m*hsj0p * xn2(1) + hi1m*hsj0p * xn2(2) + hi0m*hsj1p * xn2(3) + hi1m*hsj1p * xn2(4) &
   + hi0p*hsj0m * xn4(1) + hi1p*hsj0m * xn4(2) + hi0p*hsj1m * xn4(3) + hi1p*hsj1m * xn4(4) &
   + hi0p*hsj0p * xn3(1) + hi1p*hsj0p * xn3(2) + hi0p*hsj1p * xn3(3) + hi1p*hsj1p * xn3(4)

xrr = hrri0m*hj0m * xn1(1) + hrri1m*hj0m * xn1(2) + hrri0m*hj1m * xn1(3) + hrri1m*hj1m * xn1(4) &
    + hrri0m*hj0p * xn2(1) + hrri1m*hj0p * xn2(2) + hrri0m*hj1p * xn2(3) + hrri1m*hj1p * xn2(4) &
    + hrri0p*hj0m * xn4(1) + hrri1p*hj0m * xn4(2) + hrri0p*hj1m * xn4(3) + hrri1p*hj1m * xn4(4) &
    + hrri0p*hj0p * xn3(1) + hrri1p*hj0p * xn3(2) + hrri0p*hj1p * xn3(3) + hrri1p*hj1p * xn3(4)

xss = hi0m*hssj0m * xn1(1) + hi1m*hssj0m * xn1(2) + hi0m*hssj1m * xn1(3) + hi1m*hssj1m * xn1(4) &
    + hi0m*hssj0p * xn2(1) + hi1m*hssj0p * xn2(2) + hi0m*hssj1p * xn2(3) + hi1m*hssj1p * xn2(4) &
    + hi0p*hssj0m * xn4(1) + hi1p*hssj0m * xn4(2) + hi0p*hssj1m * xn4(3) + hi1p*hssj1m * xn4(4) &
    + hi0p*hssj0p * xn3(1) + hi1p*hssj0p * xn3(2) + hi0p*hssj1p * xn3(3) + hi1p*hssj1p * xn3(4)

xrs = hri0m*hsj0m * xn1(1) + hri1m*hsj0m * xn1(2) + hri0m*hsj1m * xn1(3) + hri1m*hsj1m * xn1(4) &
    + hri0m*hsj0p * xn2(1) + hri1m*hsj0p * xn2(2) + hri0m*hsj1p * xn2(3) + hri1m*hsj1p * xn2(4) &
    + hri0p*hsj0m * xn4(1) + hri1p*hsj0m * xn4(2) + hri0p*hsj1m * xn4(3) + hri1p*hsj1m * xn4(4) &
    + hri0p*hsj0p * xn3(1) + hri1p*hsj0p * xn3(2) + hri0p*hsj1p * xn3(3) + hri1p*hsj1p * xn3(4)

RETURN
END SUBROUTINE interp


SUBROUTINE interp1(XN1,XN2,XN3,XN4,R,S,X)
!----------------------------------------------------------------------
! SUBROUTINE CALCULATES THE INTERPOLATED VALUE OF THE FUNCTION X GIVEN
! BY XI(1..4) AT THE FOUR NODES USING BI-CUBIC HERMITE ELEMENTS
!----------------------------------------------------------------------
use itm_types
implicit none
REAL (R8) xn1(4),xn2(4),xn3(4),xn4(4)
real (r8) r,s,x,hi0m,hi1m,hj0m,hj1m,hi0p,hi1p,hj0p,hj1p

hi0m  = - (r-1.)**2 * (-r-2.) * 0.25
hi1m  = - (r-1.)**2 * (-r-1.) * 0.25

hj0m  = - (s-1.)**2 * (-s-2.) * 0.25
hj1m  = - (s-1.)**2 * (-s-1.) * 0.25

hi0p  = - (r+1.)**2 * (r-2.)  * 0.25
hi1p  = + (r+1.)**2 * (r-1.)  * 0.25

hj0p  = - (s+1.)**2 * (s-2.)  * 0.25
hj1p  = + (s+1.)**2 * (s-1.)  * 0.25

x = hi0m*hj0m * xn1(1) + hi1m*hj0m * xn1(2) + hi0m*hj1m * xn1(3) + hi1m*hj1m * xn1(4) &
  + hi0m*hj0p * xn2(1) + hi1m*hj0p * xn2(2) + hi0m*hj1p * xn2(3) + hi1m*hj1p * xn2(4) &
  + hi0p*hj0m * xn4(1) + hi1p*hj0m * xn4(2) + hi0p*hj1m * xn4(3) + hi1p*hj1m * xn4(4) &
  + hi0p*hj0p * xn3(1) + hi1p*hj0p * xn3(2) + hi0p*hj1p * xn3(3) + hi1p*hj1p * xn3(4)

RETURN
END SUBROUTINE interp1

SUBROUTINE interp2(XN1,XN2,XN3,XN4,R,S,X,XR,XS)
!----------------------------------------------------------------------
! SUBROUTINE CALCULATES THE INTERPOLATED VALUE OF THE FUNCTION X GIVEN
! BY XI(1..4) AT THE FOUR NODES USING BI-CUBIC HERMITE ELEMENTS
!----------------------------------------------------------------------
use itm_types
implicit none
REAL (R8) xn1(4),xn2(4),xn3(4),xn4(4)
real (r8) r,s,x,xr,xs,hi0m,hri0m,hi1m,hri1m,hj0m,hsj0m,hj1m,hsj1m,hi0p,hri0p,  &
     hi1p,hri1p,hj0p,hsj0p,hj1p,hsj1p

hi0m  = - (r-1.)**2 * (-r-2.) * 0.25
hri0m = - (r-1.)*(-r-2.)*0.5 +(r-1.)**2 * 0.25
hi1m  = - (r-1.)**2 * (-r-1.) * 0.25
hri1m = - (r-1.)*(-r-1.)*0.5 + (r-1.)**2 * 0.25
hj0m  = - (s-1.)**2 * (-s-2.) * 0.25
hsj0m = - (s-1.)*(-s-2.)*0.5 +(s-1.)**2  * 0.25
hj1m  = - (s-1.)**2 * (-s-1.) * 0.25
hsj1m = - (s-1.)*(-s-1.)*0.5 + (s-1.)**2 * 0.25

hi0p  = - (r+1.)**2 * (r-2.) * 0.25
hri0p = - (r+1.)*(r-2.)*0.5 - (r+1.)**2 * 0.25
hi1p  = + (r+1.)**2 * (r-1.) * 0.25
hri1p = + (r+1.)*(r-1.)*0.5 + (r+1.)**2 * 0.25

hj0p  = - (s+1.)**2 * (s-2.) * 0.25
hsj0p = - (s+1.)*(s-2.)*0.5 - (s+1.)**2 * 0.25
hj1p  = + (s+1.)**2 * (s-1.) * 0.25
hsj1p = + (s+1.)*(s-1.)*0.5 + (s+1.)**2 * 0.25

x = hi0m*hj0m * xn1(1) + hi1m*hj0m * xn1(2) + hi0m*hj1m * xn1(3) + hi1m*hj1m * xn1(4) &
  + hi0m*hj0p * xn2(1) + hi1m*hj0p * xn2(2) + hi0m*hj1p * xn2(3) + hi1m*hj1p * xn2(4) &
  + hi0p*hj0m * xn4(1) + hi1p*hj0m * xn4(2) + hi0p*hj1m * xn4(3) + hi1p*hj1m * xn4(4) &
  + hi0p*hj0p * xn3(1) + hi1p*hj0p * xn3(2) + hi0p*hj1p * xn3(3) + hi1p*hj1p * xn3(4)

xr = hri0m*hj0m * xn1(1) + hri1m*hj0m * xn1(2) + hri0m*hj1m * xn1(3) + hri1m*hj1m * xn1(4) &
   + hri0m*hj0p * xn2(1) + hri1m*hj0p * xn2(2) + hri0m*hj1p * xn2(3) + hri1m*hj1p * xn2(4) &
   + hri0p*hj0m * xn4(1) + hri1p*hj0m * xn4(2) + hri0p*hj1m * xn4(3) + hri1p*hj1m * xn4(4) &
   + hri0p*hj0p * xn3(1) + hri1p*hj0p * xn3(2) + hri0p*hj1p * xn3(3) + hri1p*hj1p * xn3(4)

xs = hi0m*hsj0m * xn1(1) + hi1m*hsj0m * xn1(2) + hi0m*hsj1m * xn1(3) + hi1m*hsj1m * xn1(4) &
   + hi0m*hsj0p * xn2(1) + hi1m*hsj0p * xn2(2) + hi0m*hsj1p * xn2(3) + hi1m*hsj1p * xn2(4) &
   + hi0p*hsj0m * xn4(1) + hi1p*hsj0m * xn4(2) + hi0p*hsj1m * xn4(3) + hi1p*hsj1m * xn4(4) &
   + hi0p*hsj0p * xn3(1) + hi1p*hsj0p * xn3(2) + hi0p*hsj1p * xn3(3) + hi1p*hsj1p * xn3(4)

RETURN
END SUBROUTINE interp2

SUBROUTINE interp3(XN1,XN2,XN3,XN4,YN1,YN2,YN3,YN4,PN1,PN2,PN3,PN4,R,S, &
                   x,xr,xs,yr,ys,ps)
!----------------------------------------------------------------------
! SUBROUTINE CALCULATES THE INTERPOLATED VALUE OF THE FUNCTION X GIVEN
! BY XI(1..4) AT THE FOUR NODES USING BI-CUBIC HERMITE ELEMENTS
!----------------------------------------------------------------------
use itm_types
implicit none
REAL (R8) xn1(4),xn2(4),xn3(4),xn4(4)
REAL (R8) yn1(4),yn2(4),yn3(4),yn4(4)
REAL (R8) pn1(4),pn2(4),pn3(4),pn4(4)
real (r8) r,s,x,xr,xs,yr,ys,ps,hi0m,hri0m,hi1m,hri1m,hj0m,hsj0m,  &
     hj1m,hsj1m,hi0p,hri0p,hi1p,hri1p,hj0p,hsj0p,hj1p,hsj1p

hi0m  = - (r-1.)**2 * (-r-2.) * 0.25
hri0m = - (r-1.)*(-r-2.)*0.5 +(r-1.)**2 * 0.25
hi1m  = - (r-1.)**2 * (-r-1.) * 0.25
hri1m = - (r-1.)*(-r-1.)*0.5 + (r-1.)**2 * 0.25
hj0m  = - (s-1.)**2 * (-s-2.) * 0.25
hsj0m = - (s-1.)*(-s-2.)*0.5 +(s-1.)**2  * 0.25
hj1m  = - (s-1.)**2 * (-s-1.) * 0.25
hsj1m = - (s-1.)*(-s-1.)*0.5 + (s-1.)**2 * 0.25

hi0p  = - (r+1.)**2 * (r-2.) * 0.25
hri0p = - (r+1.)*(r-2.)*0.5 - (r+1.)**2 * 0.25
hi1p  = + (r+1.)**2 * (r-1.) * 0.25
hri1p = + (r+1.)*(r-1.)*0.5 + (r+1.)**2 * 0.25

hj0p  = - (s+1.)**2 * (s-2.) * 0.25
hsj0p = - (s+1.)*(s-2.)*0.5 - (s+1.)**2 * 0.25
hj1p  = + (s+1.)**2 * (s-1.) * 0.25
hsj1p = + (s+1.)*(s-1.)*0.5 + (s+1.)**2 * 0.25

x = hi0m*hj0m * xn1(1) + hi1m*hj0m * xn1(2) + hi0m*hj1m * xn1(3) + hi1m*hj1m * xn1(4) &
  + hi0m*hj0p * xn2(1) + hi1m*hj0p * xn2(2) + hi0m*hj1p * xn2(3) + hi1m*hj1p * xn2(4) &
  + hi0p*hj0m * xn4(1) + hi1p*hj0m * xn4(2) + hi0p*hj1m * xn4(3) + hi1p*hj1m * xn4(4) &
  + hi0p*hj0p * xn3(1) + hi1p*hj0p * xn3(2) + hi0p*hj1p * xn3(3) + hi1p*hj1p * xn3(4)

xr = hri0m*hj0m * xn1(1) + hri1m*hj0m * xn1(2) + hri0m*hj1m * xn1(3) + hri1m*hj1m * xn1(4) &
   + hri0m*hj0p * xn2(1) + hri1m*hj0p * xn2(2) + hri0m*hj1p * xn2(3) + hri1m*hj1p * xn2(4) &
   + hri0p*hj0m * xn4(1) + hri1p*hj0m * xn4(2) + hri0p*hj1m * xn4(3) + hri1p*hj1m * xn4(4) &
   + hri0p*hj0p * xn3(1) + hri1p*hj0p * xn3(2) + hri0p*hj1p * xn3(3) + hri1p*hj1p * xn3(4)

xs = hi0m*hsj0m * xn1(1) + hi1m*hsj0m * xn1(2) + hi0m*hsj1m * xn1(3) + hi1m*hsj1m * xn1(4) &
   + hi0m*hsj0p * xn2(1) + hi1m*hsj0p * xn2(2) + hi0m*hsj1p * xn2(3) + hi1m*hsj1p * xn2(4) &
   + hi0p*hsj0m * xn4(1) + hi1p*hsj0m * xn4(2) + hi0p*hsj1m * xn4(3) + hi1p*hsj1m * xn4(4) &
   + hi0p*hsj0p * xn3(1) + hi1p*hsj0p * xn3(2) + hi0p*hsj1p * xn3(3) + hi1p*hsj1p * xn3(4)

ps = hi0m*hj0m * pn1(1) + hi1m*hj0m * pn1(2) + hi0m*hj1m * pn1(3) + hi1m*hj1m * pn1(4) &
   + hi0m*hj0p * pn2(1) + hi1m*hj0p * pn2(2) + hi0m*hj1p * pn2(3) + hi1m*hj1p * pn2(4) &
   + hi0p*hj0m * pn4(1) + hi1p*hj0m * pn4(2) + hi0p*hj1m * pn4(3) + hi1p*hj1m * pn4(4) &
   + hi0p*hj0p * pn3(1) + hi1p*hj0p * pn3(2) + hi0p*hj1p * pn3(3) + hi1p*hj1p * pn3(4)

yr = hri0m*hj0m * yn1(1) + hri1m*hj0m * yn1(2) + hri0m*hj1m * yn1(3) + hri1m*hj1m * yn1(4) &
   + hri0m*hj0p * yn2(1) + hri1m*hj0p * yn2(2) + hri0m*hj1p * yn2(3) + hri1m*hj1p * yn2(4) &
   + hri0p*hj0m * yn4(1) + hri1p*hj0m * yn4(2) + hri0p*hj1m * yn4(3) + hri1p*hj1m * yn4(4) &
   + hri0p*hj0p * yn3(1) + hri1p*hj0p * yn3(2) + hri0p*hj1p * yn3(3) + hri1p*hj1p * yn3(4)

ys = hi0m*hsj0m * yn1(1) + hi1m*hsj0m * yn1(2) + hi0m*hsj1m * yn1(3) + hi1m*hsj1m * yn1(4) &
   + hi0m*hsj0p * yn2(1) + hi1m*hsj0p * yn2(2) + hi0m*hsj1p * yn2(3) + hi1m*hsj1p * yn2(4) &
   + hi0p*hsj0m * yn4(1) + hi1p*hsj0m * yn4(2) + hi0p*hsj1m * yn4(3) + hi1p*hsj1m * yn4(4) &
   + hi0p*hsj0p * yn3(1) + hi1p*hsj0p * yn3(2) + hi0p*hsj1p * yn3(3) + hi1p*hsj1p * yn3(4)

RETURN
END SUBROUTINE interp3


SUBROUTINE rft2(DATA,NR,KR)
!*****************************************************************
! REAL FOURIER TRANSFORM.                                        *
! INPUT:  NR REAL COEFFICIENTS                                   *
!             DATA(1),DATA(1+KR),....,DATA(1+(NR-1)*KR).         *
! OUTPUT: NR/2+1 COMPLEX COEFFICIENTS                            *
!            (DATA(1),      DATA(1+KR))                          *
!            (DATA(1+2*KR), DATA(1+3*KR))                        *
!             .............................                      *
!            (DATA(1+NR*KR),DATA(1+(NR+1)*KR).                   *
! THE CALLING PROGRAM SHOULD HAVE DATA DIMENSIONED WITH AT LEAST *
! (NR+1)*KR+1 ELEMENTS. (I.E., NR+2 IF INCREMENT KR=1).          *
! LASL ROUTINE MAY 75, CALLING FFT2 AND RTRAN2.                  *
!*****************************************************************
use itm_types
implicit none
real (r8)    :: data(*)
integer :: kr,nr, ktran

CALL fft2(DATA(1),DATA(kr+1),nr/2,-(kr+kr))
CALL rtran2(DATA,nr,kr,1)
RETURN
END SUBROUTINE rft2

SUBROUTINE rtran2(DATA,NR,KR,KTRAN)
!*****************************************************************
! INTERFACE BETWEEN RFT2, RFI2, AND FFT2.                        *
! THE CALLING PROGRAM SHOULD HAVE DATA DIMENSIONED WITH AT LEAST *
! (NR+1)*KR+1 ELEMENTS.                                          *
! LASL ROUTINE MAY 75, CALLED FROM RFT2 AND RFI2.                *
!*****************************************************************
use itm_types
implicit none
real (r8)    :: data(*), theta, dc, ds, ws, wc, sumr, difr, sumi, difi
real (r8)    :: tr, ti, wca
integer :: nr, kr, ktran, ks, n, nmax, kmax, k, nk

ks=2*kr
n=nr/2
nmax=n*ks+2
kmax=nmax/2
theta=1.5707963267949/float(n)
dc=2.*sin(theta)**2
ds=sin(2.*theta)
ws=0.

IF (ktran .LE. 0) THEN
   wc=-1.0
   ds=-ds
ELSE
   wc=1.0
   DATA(nmax-1)=DATA(1)
   DATA(nmax-1+kr)=DATA(kr+1)
ENDIF
DO k=1,kmax,ks
   nk=nmax-k
   sumr=.5*(DATA(k)+DATA(nk))
   difr=.5*(DATA(k)-DATA(nk))
   sumi=.5*(DATA(k+kr)+DATA(nk+kr))
   difi=.5*(DATA(k+kr)-DATA(nk+kr))
   tr=wc*sumi-ws*difr
   ti=ws*sumi+wc*difr
   DATA(k)=sumr+tr
   DATA(k+kr)=difi-ti
   DATA(nk)=sumr-tr
   DATA(nk+kr)=-difi-ti
   wca=wc-dc*wc-ds*ws
   ws=ws+ds*wc-dc*ws
   wc=wca
enddo
return
end SUBROUTINE rtran2

SUBROUTINE fft2 (DATAR,DATAI,N,INC)
!*****************************************************************
! FFT2 FORTRAN VERSION CLAIR NIELSON MAY 75.                     *
!*****************************************************************
use itm_types
implicit none
real (r8)    :: datar(*), datai(*)
integer :: n, ninc
real (r8) tempr,tempi,theta,sinth,wstpr,wstpi,wr,wi
integer inc,ktran,ks,ip0,ip3,irev,i,ibit,ip1,ip2,i1,i3,j0,j1

ktran=isign(-1,inc)
ks=iabs(inc)
ip0=ks
ip3=ip0*n
irev=1

      DO i=1,ip3,ip0
         IF(i.LT.irev) THEN
            tempr=datar(i)
            tempi=datai(i)
            datar(i)=datar(irev)
            datai(i)=datai(irev)
            datar(irev)=tempr
            datai(irev)=tempi
         ENDIF
         ibit=ip3/2
   10    IF(irev.GT.ibit) THEN
            irev=irev-ibit
            ibit=ibit/2
            IF(ibit.GE.ip0) goto 10
         ENDIF
         irev=irev+ibit
      enddo
      ip1=ip0
      theta=REAL(ktran)*3.1415926535898
   30 IF(ip1.GE.ip3) return
      ip2=ip1+ip1
      sinth=sin(.5*theta)
      wstpr=-2.*sinth*sinth
      wstpi=sin(theta)
      wr=1.
      wi=0.
      DO i1=1,ip1,ip0
         DO i3=i1,ip3,ip2
            j0=i3
            j1=j0+ip1
            tempr=wr*datar(j1)-wi*datai(j1)
            tempi=wr*datai(j1)+wi*datar(j1)
            datar(j1)=datar(j0)-tempr
            datai(j1)=datai(j0)-tempi
            datar(j0)=datar(j0)+tempr
            datai(j0)=datai(j0)+tempi
         enddo
         tempr=wr
         wr=wr*wstpr-wi*wstpi+wr
         wi=wi*wstpr+tempr*wstpi+wi
      enddo
      ip1=ip2
      theta=.5*theta
      goto 30
RETURN
END SUBROUTINE fft2


SUBROUTINE fsum2(F,T,FFNUL,FFCOS,FFSIN,MHARM)
!-----------------------------------------------------------------------
! FOURIER SYNTHESIS OF GENERAL  FUNCTION F(T) AT SINGLE POINT T.
!-----------------------------------------------------------------------
use itm_types
implicit none
integer :: mharm, m
real (r8)    :: ffnul, ffcos(*), ffsin(*), f, t, s, c, co, ca, si, sum

co=cos(t)
si=sin(t)
c=1.
s=0.
sum=.5*ffnul
do m=1,mharm
  ca=c*co-s*si
  s=s*co+c*si
  c=ca
  sum=sum+ffcos(m)*c + ffsin(m)*s
enddo
f=sum
RETURN
END SUBROUTINE fsum2

SUBROUTINE grid2nv(TIN,TOUT,JPTS,ACC,IGRD)
!------------------------------------------------------------------------
!  THE FUNCTION TIN(TOUT), GIVEN ON THE GRID TOUT=2*PI*(J-1)/JPTS,
!  IS INVERTED TO GIVE TOUT(TIN) ON THE GRID TIN=2*PI*(I-1)/JPTS.
!  THIS IS DONE BY DETERMINING THE ZEROS OF THE FUNCTION
!     Y(T)=T+SUM(GF(M)*SIN(M*T))-2*PI*(I-1)/JPTS,
!  WHERE GF(M) ARE THE FOURIER COEFFICIENTS OF G(T)=TIN(T)-T.
!-----------------------------------------------------------------------
use itm_types
implicit none
integer jmax,ninv
parameter(jmax=1024,ninv=100)
dimension tin(*),tout(*),t(jmax+1),g(jmax+1),gfcos(jmax/2-1),gfsin(jmax/2-1)
equivalence(t(1),g(1))
real (r8) tin,tout,t,g,gfcos,gfsin,acc,pi,t1,sum1,y1,t0,y0,t2,sum2,y2,gfnul
integer jpts,igrd,mharm,jj,i,j1,igrdnv,ifirst,icirc,j,n

pi=2. * asin(1.)
mharm=jpts/2-1

DO jj=2,jpts
  IF (tin(jj-1).GT.tin(jj))  tin(jj)=tin(jj)+2*pi
ENDDO

DO i=1,jpts
  g(i)=tin(i)-2.*pi*(i-1.)/jpts
ENDDO

CALL rft(g,gfnul,gfcos,gfsin,jpts,mharm)

do i=1,jpts+1
  t(i)=2.*pi*(i-1.)/jpts
enddo
j1=1
igrdnv=1
ifirst=0
icirc= - (int(tin(1)/(2*pi)+10000) -  9999)
IF (abs(tin(1)).LT.1e-12) icirc=0

      DO 80 i=1,jpts
        j=j1
        t1=t(j) + icirc*2*pi
        CALL fsum2(sum1,t1,gfnul,gfcos,gfsin,mharm)
        y1=t1+sum1-t(i)
   30   CONTINUE
          t0=t1
          y0=y1
          IF (abs(y0).LE.acc) THEN
            tout(i)=t0
            goto 80
          ENDIF
          IF (j.NE.jpts+1) goto 31
            IF (ifirst.EQ.0)  THEN
              j=1
              icirc=icirc+1
              ifirst=1
            ELSE
              goto 90
            ENDIF
   31     j=j+1
          t1=t(j) + icirc*2*pi
          CALL fsum2(sum1,t1,gfnul,gfcos,gfsin,mharm)
          y1=t1+sum1-t(i)
          IF(sign(1.0_r8,y0).EQ.sign(1.0_r8,y1)) goto 30
        j1=j-1
        DO 40 n=1,ninv
          t2=t0-(t1-t0)*y0/(y1-y0)
          CALL fsum2(sum2,t2,gfnul,gfcos,gfsin,mharm)
          y2=t2+sum2-t(i)
          IF(abs(y2).LE.acc) goto 50
          IF(sign(1.0_r8,y2).EQ.sign(1.0_r8,y1)) goto 45
          t0=t2
          y0=y2
          goto 40
   45     t1=t2
          y1=y2
   40   CONTINUE
   50   tout(i)=t2
        IF(n.GT.igrdnv) igrdnv=n
   80 CONTINUE
   90 RETURN

END SUBROUTINE grid2nv

SUBROUTINE rft(F,FFNUL,FFCOS,FFSIN,JPTS,MHARM)
!-----------------------------------------------------------------------
!  CALCULATES FOURIER COSINE AND SINE COEFFICIENTS FFCOS AND
!  FFSIN OF THE ARRAY FF CORRESPONDING TO THE  FUNCTION
!  F(T)=.5*FFNUL+SUM(FFCOS(M)*COS(M*T)+FFSIN(M)*SIN(M*T))
!  WHERE MHARM.LE.JPTS/2-1, FFNUL=FF(0) AND T=2*PI*(J-1)/JPTS.
!  THE INPUT ARRAY F(J) IS NOT DESTROYED BY CALLING RFTCOS.
!  TYPICAL USE IS FOR MHARM MUCH SMALLER THAN JPTS/2-1, SO THAT
!  RFT2 CANNOT BE USED DIRECTLY.
!-----------------------------------------------------------------------
use itm_types
implicit none
integer :: jpts, jmax, mharm, j, m
parameter(jmax=1024)
real (r8)    :: f(*),ffnul,ffcos(*),ffsin(*),fstore(jmax+2)
real (r8)    :: fac

DO j=1,jpts
  fstore(j)=f(j)
enddo
CALL rft2(fstore,jpts,1)
fac=2./jpts
ffnul=fstore(1)*fac
DO m=1,mharm
  ffcos(m)=fstore(2*m+1)*fac
  ffsin(m) = - fstore(2*m+2)*fac
ENDDO
RETURN
END SUBROUTINE rft

SUBROUTINE spline(N,X,Y,ALFA,BETA,TYP,A,B,C,D)
!-----------------------------------------------------------------------
!     INPUT:
!
!     N     ANZAHL DER KNOTEN
!     X     ARRAY DER X-WERTE
!     Y     ARRAY DER Y-WERTE
!     ALFA  RANDBEDINGUNG IN X(1)
!     BETA        "       IN X(N)
!     TYP   =  0  NOT-A-KNOT SPLINE
!              1  ALFA, BETA 1. ABLEITUNGEN VORGEGEBEN
!              2    "    "   2.     "           "
!              3    "    "   3.     "           "
!
!     BEMERKUNG: MIT TYP = 2 UND ALFA = BETA = 0 ERHAELT MAN
!           EINEN NATUERLICHEN SPLINE
!
!     OUTPUT:
!
!     A, B, C, D     ARRAYS DER SPLINEKOEFFIZIENTEN
!       S = A(I) + B(I)*(X-X(I)) + C(I)*(X-X(I))**2+ D(I)*(X-X(I))**3
!
!     BEI ANWENDUNGSFEHLERN WIRD DAS PROGRAMM MIT ENTSPRECHENDER
!     FEHLERMELDUNG ABGEBROCHEN
!-----------------------------------------------------------------------
use itm_types
implicit none
      INTEGER  n, typ
      REAL (R8)     x(n), y(n), alfa, beta, a(n), b(n), c(n), d(n)
      INTEGER  i, ierr
      REAL (R8)     h(1001)


      IF((typ.LT.0).OR.(typ.GT.3)) THEN
         WRITE(*,*) ' ERROR IN ROUTINE SPLINE: WRONG TYPE'
         stop
      ENDIF

      IF (n.LT.3) THEN
         WRITE(*,*) ' ERROR IN ROUTINE  SPLINE: N < 3 '
         stop
      ENDIF

!     BERECHNE DIFFERENZ AUFEINENDERFOLGENDER X-WERTE UND
!     UNTERSUCHE MONOTONIE

      DO i = 1, n-1
        h(i) = x(i+1)- x(i)
        IF ( h(i).LE. 0.0) THEN
         WRITE(*,*) ' NON-MONOTONIC COORDINATE IN SPLINE: X(I-1)>=X(I)'
         stop
       ENDIF
      ENDDO

!     AUFSTELLEN DES GLEICHUNGSSYSTEMS

      DO i = 1, n-2
         a(i) = 3.0 * ((y(i+2)-y(i+1)) / h(i+1) - (y(i+1)-y(i)) / h(i))
         b(i) = h(i)
         c(i) = h(i+1)
         d(i) = 2.0 * (h(i) + h(i+1))
      ENDDO

!     BERUECKSICHTIGEN DER RANDBEDINGUNGEN
!
!     NOT-A-KNOT

      IF (typ.EQ.0) THEN
         a(1)   = a(1) * h(2) / (h(1) + h(2))
         a(n-2) = a(n-2) * h(n-2) / (h(n-1) + h(n-2))
         d(1)   = d(1) - h(1)
         d(n-2) = d(n-2) - h(n-1)
         c(1)   = c(1) - h(1)
         b(n-2) = b(n-2) - h(n-1)
      ENDIF

!     1. ABLEITUNG VORGEGEBEN

      IF (typ.EQ.1) THEN
         a(1)   = a(1) - 1.5 * ((y(2)-y(1)) / h(1) - alfa)
         a(n-2) = a(n-2) - 1.5 * (beta - (y(n)-y(n-1)) / h(n-1))
         d(1)   = d(1) - 0.5 * h(1)
         d(n-2) = d(n-2) - 0.5 * h(n-1)
      ENDIF

!     2. ABLEITUNG VORGEGEBEN

      IF (typ.EQ.2) THEN
         a(1)   = a(1) - 0.5 * alfa * h(1)
         a(n-2) = a(n-2) - 0.5 * beta * h(n-1)
      ENDIF

!     3. ABLEITUNG VORGEGEBEN

      IF (typ.EQ.3 ) THEN
         a(1)   = a(1) + 0.5 * alfa * h(1) * h(1)
         a(n-2) = a(n-2) - 0.5 * beta * h(n-1)* h(n-1)
         d(1)   = d(1) + h(1)
         d(n-2) = d(n-2) + h(n-1)
      ENDIF

!     BERECHNUNG DER KOEFFIZIENTEN

      CALL sgtsl(n-2,b,d,c,a,ierr)

      IF(ierr.NE.0) THEN
         WRITE(*,21)
         stop
      ENDIF

!     UEBERSCHREIBEN DES LOESUNGSVEKTORS

      CALL dcopy(n-2,a,1,c(2),1)

!     IN ABHAENGIGKEIT VON DEN RANDBEDINGUNGEN WIRD DER 1. UND
!     DER LETZTE WERT VON C KORRIGIERT

      IF (typ.EQ.0) THEN
         c(1) = c(2) + h(1) * (c(2)-c(3)) / h(2)
         c(n) = c(n-1) + h(n-1) * (c(n-1)-c(n-2)) / h(n-2)
      ENDIF

      IF (typ.EQ.1) THEN
         c(1) = 1.5*((y(2)-y(1)) / h(1) - alfa) / h(1) - 0.5 * c(2)
         c(n) = -1.5*((y(n)-y(n-1)) / h(n-1)-beta) / h(n-1)-0.5*c(n-1)
      ENDIF

      IF (typ.EQ.2) THEN
         c(1) = 0.5 * alfa
         c(n) = 0.5 * beta
      ENDIF

      IF (typ.EQ.3) THEN
         c(1) = c(2) - 0.5 * alfa * h(1)
         c(n) = c(n-1) + 0.5 * beta * h(n-1)
      ENDIF

      CALL dcopy(n,y,1,a,1)

      DO i = 1, n-1
         b(i) = (a(i+1)-a(i)) / h(i) - h(i) * (c(i+1)+2.0 * c(i)) / 3.0
         d(i) = (c(i+1)-c(i)) / (3.0 * h(i))
      ENDDO

      b(n) = (3.0 * d(n-1) * h(n-1) + 2.0 * c(n-1)) * h(n-1) + b(n-1)

      RETURN

   21 FORMAT(1x,'ERROR IN SGTSL: MATRIX SINGULAR')
    END SUBROUTINE spline



FUNCTION spwert(N,XWERT,A,B,C,D,X,ABLTG)
!-----------------------------------------------------------------------
!     INPUT:
!
!     N           ANZAHL DER KNOTENPUNKTE
!     XWERT       STELLE AN DER FUNKTIONSWERTE BERECHNET WERDEN
!     A, B, C, D  ARRAYS DER SPLINEKOEFFIZIENTEN (AUS SPLINE)
!     X           ARRAY DER KNOTENPUNKTE
!
!     OUTPUT:
!
!     SPWERT   FUNKTIONSWERT AN DER STELLE XWERT
!     ABLTG(I) WERT DER I-TEN ABLEITUNG BEI XWERT
!-----------------------------------------------------------------------
use itm_types
implicit none
REAL (R8) spwert
      INTEGER  n
      REAL (R8)     xwert, a(n), b(n), c(n), d(n), x(n), abltg(3)
      INTEGER  i, k, m
      real (r8) xx

!     SUCHE PASSENDES INTERVALL (BINAERE SUCHE)

      i = 1
      k = n

   10 m = (i+k) / 2

      IF(m.NE.i) THEN
         IF(xwert.LT.x(m)) THEN
            k = m
         ELSE
            i = m
         ENDIF
         goto 10
      ENDIF

      xx = xwert - x(i)

      abltg(1) = (3.0 * d(i) * xx + 2.0 * c(i)) * xx + b(i)
      abltg(2) = 6.0 * d(i) * xx + 2.0 * c(i)
      abltg(3) = 6.0 * d(i)

      spwert = ((d(i)*xx + c(i))*xx + b(i))*xx + a(i)

      RETURN
    END FUNCTION spwert

SUBROUTINE sgtsl(N,C,D,E,B,INFO)
use itm_types
implicit none
INTEGER n,info
REAL (R8) c(*),d(*),e(*),b(*)
!
!     SGTSL GIVEN A GENERAL TRIDIAGONAL MATRIX AND A RIGHT HAND
!     SIDE WILL FIND THE SOLUTION.
!
!     ON ENTRY
!
!        N       INTEGER
!                IS THE ORDER OF THE TRIDIAGONAL MATRIX.
!
!        C       REAL(N)
!                IS THE SUBDIAGONAL OF THE TRIDIAGONAL MATRIX.
!                C(2) THROUGH C(N) SHOULD CONTAIN THE SUBDIAGONAL.
!                ON OUTPUT C IS DESTROYED.
!
!        D       REAL(N)
!                IS THE DIAGONAL OF THE TRIDIAGONAL MATRIX.
!                ON OUTPUT D IS DESTROYED.
!
!        E       REAL(N)
!                IS THE SUPERDIAGONAL OF THE TRIDIAGONAL MATRIX.
!                E(1) THROUGH E(N-1) SHOULD CONTAIN THE SUPERDIAGONAL.
!                ON OUTPUT E IS DESTROYED.
!
!        B       REAL(N)
!                IS THE RIGHT HAND SIDE VECTOR.
!
!     ON RETURN
!
!        B       IS THE SOLUTION VECTOR.
!
!        INFO    INTEGER
!                = 0 NORMAL VALUE.
!                = K IF THE K-TH ELEMENT OF THE DIAGONAL BECOMES
!                    EXACTLY ZERO.  THE SUBROUTINE RETURNS WHEN
!                    THIS IS DETECTED.
!
!     LINPACK. THIS VERSION DATED 08/14/78 .
!     JACK DONGARRA, ARGONNE NATIONAL LABORATORY.
!
!     NO EXTERNALS
!     FORTRAN ABS
!
!     INTERNAL VARIABLES
!
      INTEGER k,kb,kp1,nm1,nm2
      REAL (R8) t
!     BEGIN BLOCK PERMITTING ...EXITS TO 100

         info = 0
         c(1) = d(1)
         nm1 = n - 1
         IF (nm1 .LT. 1) go to 40
            d(1) = e(1)
            e(1) = 0.0e0
            e(n) = 0.0e0

            DO 30 k = 1, nm1
               kp1 = k + 1

!              FIND THE LARGEST OF THE TWO ROWS

               IF (abs(c(kp1)) .LT. abs(c(k))) go to 10

!                 INTERCHANGE ROW

                  t = c(kp1)
                  c(kp1) = c(k)
                  c(k) = t
                  t = d(kp1)
                  d(kp1) = d(k)
                  d(k) = t
                  t = e(kp1)
                  e(kp1) = e(k)
                  e(k) = t
                  t = b(kp1)
                  b(kp1) = b(k)
                  b(k) = t
   10          CONTINUE

!              ZERO ELEMENTS

               IF (c(k) .NE. 0.0e0) go to 20
                  info = k
!     ............EXIT
                  go to 100
   20          CONTINUE
               t = -c(kp1)/c(k)
               c(kp1) = d(kp1) + t*d(k)
               d(kp1) = e(kp1) + t*e(k)
               e(kp1) = 0.0e0
               b(kp1) = b(kp1) + t*b(k)
   30       CONTINUE
   40    CONTINUE
         IF (c(n) .NE. 0.0e0) go to 50
            info = n
         go to 90
   50    CONTINUE

!           BACK SOLVE

            nm2 = n - 2
            b(n) = b(n)/c(n)
            IF (n .EQ. 1) go to 80
               b(nm1) = (b(nm1) - d(nm1)*b(n))/c(nm1)
               IF (nm2 .LT. 1) go to 70
               DO 60 kb = 1, nm2
                  k = nm2 - kb + 1
                  b(k) = (b(k) - d(k)*b(k+1) - e(k)*b(k+2))/c(k)
   60          CONTINUE
   70          CONTINUE
   80       CONTINUE
   90    CONTINUE
  100 CONTINUE

RETURN
END SUBROUTINE sgtsl

SUBROUTINE tb15a(N,X,F,D,W,LP)
!------------------------------------------------------------------
! HSL routine for cubic spline with periodic boundary conditions
! first point must be the same as last : f(1)=f(n)
!    N : number of points
!    X : coordinate (input)
!    F : the function values to be splined (input)
!    D : the derivatives at the points (output)
!    W : workspace (dimension 3N)
!   LP : unit number for output
!------------------------------------------------------------------
use itm_types
implicit none
REAL (R8) zero,one,two,three
parameter(zero=0.0e0,one=1.0e0,two=2.0e0,three=3.0e0)
INTEGER lp,n
REAL (R8) d(n),f(n),w(*),x(n)
REAL (R8) a3n1,f1,f2,h1,h2,p
INTEGER i,j,k,n2

WRITE(*,*) f(1),f(n)
IF (n.LT.4) THEN
  WRITE (lp,'(A39)')  'RETURN FROM TB15AD BECAUSE N TOO SMALL'
  w(1) = one
  RETURN
END IF
DO i = 2,n
  IF (x(i).LE.x(i-1)) THEN
    WRITE (lp,'(A29,I3,A13)') ' RETURN FROM TB15AD BECAUSE  ',i,' OUT OF ORDER'
    write(*,*) x
    w(1) = two
    RETURN
  END IF
ENDDO
IF (f(1).NE.f(n)) THEN
  WRITE (lp,'(A40)')  .NE.'RETURN FROM TB15AD BECAUSE F(1)F(N)'
  w(1) = three
  RETURN
END IF
DO i = 2,n
  h1 = one/ (x(i)-x(i-1))
  f1 = f(i-1)
  IF (i.EQ.n) THEN
    h2 = one/ (x(2)-x(1))
    f2 = f(2)
  ELSE
    h2 = one/ (x(i+1)-x(i))
    f2 = f(i+1)
  END IF
  w(3*i-2) = h1
  w(3*i-1) = two* (h1+h2)
  w(3*i) = h2
  d(i) = 3.0* (f2*h2*h2+f(i)* (h1*h1-h2*h2)-f1*h1*h1)
ENDDO
n2 = n - 2
k = 5
a3n1 = w(3*n-1)
DO i = 2,n2
  p = w(k+2)/w(k)
  w(k+3) = w(k+3) - p*w(k+1)
  d(i+1) = d(i+1) - p*d(i)
  w(k+2) = -p*w(k-1)
  p = w(k-1)/w(k)
  a3n1 = -p*w(k-1) + a3n1
  d(n) = d(n) - p*d(i)
  k = k + 3
ENDDO
p = (w(k+2)+w(k-1))/w(k)
a3n1 = a3n1 - p* (w(k+1)+w(k-1))
d(n) = (d(n)-p*d(n-1))/a3n1
DO i = 3,n
  j = n + 2 - i
  d(j) = (d(j)-w(3*j)*d(j+1)-w(3*j-2)*d(n))/w(3*j-1)
ENDDO
d(1) = d(n)
w(1) = zero
RETURN
END SUBROUTINE tb15a



SUBROUTINE tg02a(IX,N,U,S,D,X,V)
!------------------------------------------------------------------
! HSL subroutine to calculate splined values
!    N  : number of points
!    IX : negative 0 -> no initial guess for where xi is
!        positive -> gues for index close to value X
!   U   : the coordinates of the spline points
!   S   : the function values of the spline points
!   D   : the derivatives on the spline points
!   X   : the coordinate where the output is wanted
!   V(1-4) : value and derivatives of the spline interpolation
!------------------------------------------------------------------
use itm_types
implicit none
REAL (R8) x
INTEGER ix,n
REAL (R8) d(*),s(*),u(*),v(*)
REAL (R8) a,b,c,c3,e0_r8,d1,eps,gama,h,hr,hrr,phi,s0,s1,t,theta
INTEGER i,iflg,j,k

eps = 1.e-33
k = 0
iflg = 0
IF (x.LT.u(1)) go to 990
IF (x.GT.u(n)) go to 991
IF (ix.LT.0 .OR. iflg.EQ.0) go to 12
IF (x.GT.u(j+1)) go to 1
IF (x.GE.u(j)) go to 18
go to 2

    1 j = j + 1
   11 IF (x.GT.u(j+1)) go to 1
      go to 7
   12 j = abs(x-u(1))/ (u(n)-u(1))* (n-1) + 1
      j = min(j,n-1)
      iflg = 1
      IF (x.GE.u(j)) go to 11
    2 j = j - 1
      IF (x.LT.u(j)) go to 2
    7 k = j
      h = u(j+1) - u(j)
      hr = 1./h
      hrr = (hr+hr)*hr
      s0 = s(j)
      s1 = s(j+1)
      e0_r8 = d(j)
      d1 = d(j+1)
      a = s1 - s0
      b = a - h*d1
      a = a - h*e0_r8
      c = a + b
      c3 = c*3.
   18 theta = (x-u(j))*hr
      phi = 1. - theta
      t = theta*phi
      gama = theta*b - phi*a
      v(1) = theta*s1 + phi*s0 + t*gama
      v(2) = theta*d1 + phi*e0_r8 + t*c3*hr
      v(3) = (c* (phi-theta)-gama)*hrr
      v(4) = -c3*hrr*hr
      RETURN
  990 IF (x.LE.u(1)-eps*max(abs(u(1)),abs(u(n)))) go to 99
      j = 1
      go to 7
  991 IF (x.GE.u(n)+eps*max(abs(u(1)),abs(u(n)))) go to 995
      j = n - 1
      go to 7
  995 k = n
   99 iflg = 0
      DO i = 1,4
        v(i) = 0.
      ENDDO
RETURN
END SUBROUTINE tg02a

SUBROUTINE qsort2 (ORD,N,A)
use itm_types
implicit none
!
!==============SORTS THE ARRAY A(I),I=1,2,...,N BY PUTTING THE
!   ASCENDING ORDER VECTOR IN ORD.  THAT IS ASCENDING ORDERED A
!   IS A(ORD(I)),I=1,2,...,N; DESCENDING ORDER A IS A(ORD(N-I+1)),
!   I=1,2,...,N .  THIS SORT RUNS IN TIME PROPORTIONAL TO N LOG N .
!
!
!     ACM QUICKSORT - ALGORITHM #402 - IMPLEMENTED IN FORTRAN BY
!                                 WILLIAM H. VERITY
!                                 COMPUTATION CENTER
!                                 PENNSYLVANIA STATE UNIVERSITY
!                                 UNIVERSITY PARK, PA.  16802
!     With correction to that algorithm.
!
!
      dimension ord(n),poplst(2,20)
      integer ord,poplst
!
!     To sort different input types change the following
!     specification statements; FOR EXAMPLE,  REAL A(N) or
!     CHARACTER *(L) A(N)  for REAL or CHARACTER sorting
!     respectively  similarly for X,XX,Z,ZZ,Y. L is the
!     character length of the elements of A.
!
      REAL (R8) a(n)
      REAL (R8) x,xx,z,zz,y
      integer n,ndeep,u1,l1,i,l,u,p,q,yp,ix,iz,ip,iq
!
      ndeep=0
      u1=n
      l1=1
      DO 1  i=1,n
    1 ord(i)=i
    2 IF (u1.GT.l1) go to 3
      RETURN
!
    3 l=l1
      u=u1
!
! PART
!
    4 p=l
      q=u
      x=a(ord(p))
      z=a(ord(q))
      IF (x.LE.z) go to 5
      y=x
      x=z
      z=y
      yp=ord(p)
      ord(p)=ord(q)
      ord(q)=yp
    5 IF (u-l.LE.1) go to 15
      xx=x
      ix=p
      zz=z
      iz=q
!
! LEFT
!
    6 p=p+1
      IF (p.GE.q) go to 7
      x=a(ord(p))
      IF (x.GE.xx) go to 8
      go to 6
    7 p=q-1
      go to 13
!
! RIGHT
!
    8 q=q-1
      IF (q.LE.p) go to 9
      z=a(ord(q))
      IF (z.LE.zz) go to 10
      go to 8
    9 q=p
      p=p-1
      z=x
      x=a(ord(p))
!
! DIST
!
   10 IF (x.LE.z) go to 11
      y=x
      x=z
      z=y
      ip=ord(p)
      ord(p)=ord(q)
      ord(q)=ip
   11 IF (x.LE.xx) go to 12
      xx=x
      ix=p
   12 IF (z.GE.zz) go to 6
      zz=z
      iz=q
      go to 6
!
! OUT
!
   13 CONTINUE
      IF (.NOT.(p.NE.ix.AND.x.NE.xx)) go to 14
      ip=ord(p)
      ord(p)=ord(ix)
      ord(ix)=ip
   14 CONTINUE
      IF (.NOT.(q.NE.iz.AND.z.NE.zz)) go to 15
      iq=ord(q)
      ord(q)=ord(iz)
      ord(iz)=iq
   15 CONTINUE
      IF (u-q.LE.p-l) go to 16
      l1=l
      u1=p-1
      l=q+1
      go to 17
   16 u1=u
      l1=q+1
      u=p-1
   17 CONTINUE
      IF (u1.LE.l1) go to 18
!
! START RECURSIVE CALL
!
      ndeep=ndeep+1
      poplst(1,ndeep)=u
      poplst(2,ndeep)=l
      go to 3
   18 IF (u.GT.l) go to 4
!
! POP BACK UP IN THE RECURSION LIST
!
      IF (ndeep.EQ.0) go to 2
      u=poplst(1,ndeep)
      l=poplst(2,ndeep)
      ndeep=ndeep-1
      go to 18
!
! END QSORT
END SUBROUTINE qsort2

end module helena21_mod