gray/src/equilibrium.f90

1280 lines
39 KiB
Fortran
Raw Normal View History

! This module handles the loading, interpolation and evaluation of the
! MHD equilibrium data (poloidal current function, safety factor,
! poloidal flux, magnetic field, etc.)
!
! Two kinds of data are supported: analytical (suffix `_an` in the
! subroutine names) or numerical (suffix `_spline`). For the latter, the
! the data is interpolated using splines.
2015-11-18 17:34:33 +01:00
module equilibrium
use const_and_precisions, only : wp_
use splines, only : spline_simple, spline_1d, spline_2d, linear_1d
2015-11-18 17:34:33 +01:00
implicit none
! Parameters of the analytical equilibrium model
type analytic_model
real(wp_) :: q0 ! Safety factor at the magnetic axis
real(wp_) :: q1 ! Safety factor at the edge
real(wp_) :: lq ! Exponent for the q(ρ) power law
end type
! Order of the splines
integer, parameter :: kspl=3, ksplp=kspl + 1
! Global variable storing the state of the module
! Splines
type(spline_1d), save :: fpol_spline ! Poloidal current function F(ψ)
type(spline_2d), save :: psi_spline ! Poloidal flux ψ(R,z)
type(spline_simple), save :: q_spline ! Safey factor q(ψ)
type(linear_1d), save :: rhop_spline, rhot_spline ! Normalised radii ρ_p(ρ_t), ρ_t(ρ_p)
! Analytic model
type(analytic_model), save :: model
! More parameters
real(wp_), save :: fpolas
real(wp_), save :: psia, psiant, psinop
real(wp_), save :: phitedge, aminor
real(wp_), save :: btaxis,rmaxis,zmaxis,sgnbphi
real(wp_), save :: btrcen ! used only for jcd_astra def.
real(wp_), save :: rcen ! computed as fpol(a)/btrcen
real(wp_), save :: rmnm, rmxm, zmnm, zmxm
real(wp_), save :: zbinf, zbsup
private
public read_eqdsk, read_equil_an ! Reading data files
public scale_equil, change_cocos ! Transforming data
public equian ! Analytical model
public equinum_psi, equinum_fpol ! Interpolated data
public fq, bfield, tor_curr, tor_curr_psi ! Accessing local quantities
public frhopol, frhopolv, frhotor ! Converting between poloidal/toroidal flux
public set_equil_spline, set_equil_an ! Initialising internal state
public unset_equil_spline ! Deinitialising internal state
! Members exposed for magsurf_data
public kspl, psi_spline, q_spline, points_tgo
! Members exposed to gray_core and more
public psia, psiant, psinop
public phitedge, aminor
public btaxis,rmaxis,zmaxis,sgnbphi
public btrcen, rcen
public rmnm, rmxm, zmnm, zmxm
public zbinf, zbsup
2015-11-18 17:34:33 +01:00
contains
subroutine read_eqdsk(params, data, unit)
! Reads the MHD equilibrium `data` from a G-EQDSK file (params%filenm).
! If given, the file is opened in the `unit` number.
! For a description of the G-EQDSK, see the GRAY user manual.
2015-11-18 17:34:33 +01:00
use const_and_precisions, only : one
use gray_params, only : equilibrium_parameters, equilibrium_data
use utils, only : get_free_unit
2021-12-18 18:57:38 +01:00
use logger, only : log_error
2015-11-18 17:34:33 +01:00
implicit none
! subroutine arguments
type(equilibrium_parameters), intent(in) :: params
type(equilibrium_data), intent(out) :: data
integer, optional, intent(in) :: unit
! local variables
integer :: u, idum, i, j, nr, nz, nbnd, nlim
2015-11-18 17:34:33 +01:00
character(len=48) :: string
real(wp_) :: dr, dz, dps, rleft, zmid, zleft, psiedge, psiaxis
real(wp_) :: xdum ! dummy variable, used to discard data
2021-12-18 18:57:38 +01:00
integer :: err
u = get_free_unit(unit)
2015-11-18 17:34:33 +01:00
! Open the G-EQDSK file
2021-12-18 18:57:38 +01:00
open(file=params%filenm, status='old', action='read', unit=u, iostat=err)
if (err /= 0) then
call log_error('opening eqdsk file ('//trim(params%filenm)//') failed!', &
mod='equilibrium', proc='read_eqdsk')
call exit(1)
end if
2015-11-18 17:34:33 +01:00
! get size of main arrays and allocate them
if (params%idesc == 1) then
2015-11-18 17:34:33 +01:00
read (u,'(a48,3i4)') string,idum,nr,nz
else
read (u,*) nr, nz
2015-11-18 17:34:33 +01:00
end if
if (allocated(data%rv)) deallocate(data%rv)
if (allocated(data%zv)) deallocate(data%zv)
if (allocated(data%psin)) deallocate(data%psin)
if (allocated(data%psinr)) deallocate(data%psinr)
if (allocated(data%fpol)) deallocate(data%fpol)
if (allocated(data%qpsi)) deallocate(data%qpsi)
allocate(data%rv(nr), data%zv(nz), &
data%psin(nr, nz), &
data%psinr(nr), &
data%fpol(nr), &
data%qpsi(nr))
! Store 0D data and main arrays
if (params%ifreefmt==1) then
read (u, *) dr, dz, data%rvac, rleft, zmid
read (u, *) data%rax, data%zax, psiaxis, psiedge, xdum
read (u, *) xdum, xdum, xdum, xdum, xdum
read (u, *) xdum, xdum, xdum, xdum, xdum
read (u, *) (data%fpol(i), i=1,nr)
read (u, *) (xdum,i=1, nr)
read (u, *) (xdum,i=1, nr)
read (u, *) (xdum,i=1, nr)
read (u, *) ((data%psin(i,j), i=1,nr), j=1,nz)
read (u, *) (data%qpsi(i), i=1,nr)
2015-11-18 17:34:33 +01:00
else
read (u, '(5e16.9)') dr,dz,data%rvac,rleft,zmid
read (u, '(5e16.9)') data%rax,data%zax,psiaxis,psiedge,xdum
read (u, '(5e16.9)') xdum,xdum,xdum,xdum,xdum
read (u, '(5e16.9)') xdum,xdum,xdum,xdum,xdum
read (u, '(5e16.9)') (data%fpol(i),i=1,nr)
read (u, '(5e16.9)') (xdum,i=1,nr)
read (u, '(5e16.9)') (xdum,i=1,nr)
read (u, '(5e16.9)') (xdum,i=1,nr)
read (u, '(5e16.9)') ((data%psin(i,j),i=1,nr),j=1,nz)
read (u, '(5e16.9)') (data%qpsi(i),i=1,nr)
2015-11-18 17:34:33 +01:00
end if
! Get size of boundary and limiter arrays and allocate them
read (u,*) nbnd, nlim
if (allocated(data%rbnd)) deallocate(data%rbnd)
if (allocated(data%zbnd)) deallocate(data%zbnd)
if (allocated(data%rlim)) deallocate(data%rlim)
if (allocated(data%zlim)) deallocate(data%zlim)
! Load plasma boundary data
if(nbnd > 0) then
allocate(data%rbnd(nbnd), data%zbnd(nbnd))
if (params%ifreefmt == 1) then
read(u, *) (data%rbnd(i), data%zbnd(i), i=1,nbnd)
2015-11-18 17:34:33 +01:00
else
read(u, '(5e16.9)') (data%rbnd(i), data%zbnd(i), i=1,nbnd)
2015-11-18 17:34:33 +01:00
end if
end if
! Load limiter data
if(nlim > 0) then
allocate(data%rlim(nlim), data%zlim(nlim))
if (params%ifreefmt == 1) then
read(u, *) (data%rlim(i), data%zlim(i), i=1,nlim)
2015-11-18 17:34:33 +01:00
else
read(u, '(5e16.9)') (data%rlim(i), data%zlim(i), i=1,nlim)
2015-11-18 17:34:33 +01:00
end if
end if
! End of G-EQDSK file
2015-11-18 17:34:33 +01:00
close(u)
! Build rv,zv,psinr arrays
zleft = zmid-0.5_wp_*dz
dr = dr/(nr-1)
dz = dz/(nz-1)
dps = one/(nr-1)
2015-11-18 17:34:33 +01:00
do i=1,nr
data%psinr(i) = (i-1)*dps
data%rv(i) = rleft + (i-1)*dr
2015-11-18 17:34:33 +01:00
end do
do i=1,nz
data%zv(i) = zleft + (i-1)*dz
2015-11-18 17:34:33 +01:00
end do
! Normalize psin
data%psia = psiedge - psiaxis
if(params%ipsinorm == 0) data%psin = (data%psin - psiaxis)/data%psia
end subroutine read_eqdsk
2015-11-18 17:34:33 +01:00
subroutine read_equil_an(filenm, ipass, data, unit)
! Reads the MHD equilibrium `data` in the analytical format
! from params%filenm.
! If given, the file is opened in the `unit` number.
!
! TODO: add format description
use gray_params, only : equilibrium_data
use utils, only : get_free_unit
use logger, only : log_error
2021-12-18 18:57:38 +01:00
2015-11-18 17:34:33 +01:00
implicit none
! subroutine arguments
character(len=*), intent(in) :: filenm
integer, intent(in) :: ipass
type(equilibrium_data), intent(out) :: data
integer, optional, intent(in) :: unit
! local variables
integer :: i, u, nlim
2021-12-18 18:57:38 +01:00
integer :: err
real(wp_) :: rr0m, zr0m, rpam, b0
2015-11-18 17:34:33 +01:00
u = get_free_unit(unit)
2021-12-18 18:57:38 +01:00
open(file=filenm, status='old', action='read', unit=u, iostat=err)
if (err /= 0) then
call log_error('opening equilibrium file ('//trim(filenm)//') failed!', &
mod='equilibrium', proc='read_equil_an')
call exit(1)
end if
read(u, *) rr0m, zr0m, rpam
read(u, *) b0
read(u, *) model%q0, model%q1, model%lq
if(allocated(data%rv)) deallocate(data%rv)
if(allocated(data%zv)) deallocate(data%zv)
if(allocated(data%fpol)) deallocate(data%fpol)
if(allocated(data%qpsi)) deallocate(data%qpsi)
allocate(data%rv(2), data%zv(1), data%fpol(1), data%qpsi(3))
data%rv = [rr0m, rpam]
data%zv = [zr0m]
data%fpol = [b0*rr0m]
data%qpsi = [model%q0, model%q1, model%lq]
if(ipass >= 2) then
! get size of boundary and limiter arrays and allocate them
read (u,*) nlim
if (allocated(data%rlim)) deallocate(data%rlim)
if (allocated(data%zlim)) deallocate(data%zlim)
! store boundary and limiter data
if(nlim > 0) then
allocate(data%rlim(nlim), data%zlim(nlim))
read(u,*) (data%rlim(i), data%zlim(i), i = 1, nlim)
end if
end if
2015-11-18 17:34:33 +01:00
close(u)
end subroutine read_equil_an
subroutine change_cocos(data, cocosin, cocosout, error)
! Convert the MHD equilibrium data from one coordinate convention
! (COCOS) to another. These are specified by `cocosin` and
! `cocosout`, respectively.
!
! For more information, see: https://doi.org/10.1016/j.cpc.2012.09.010
use const_and_precisions, only : zero, one, pi
use gray_params, only : equilibrium_data
2015-11-18 17:34:33 +01:00
implicit none
! subroutine arguments
type(equilibrium_data), intent(inout) :: data
2015-11-18 17:34:33 +01:00
integer, intent(in) :: cocosin, cocosout
integer, intent(out), optional :: error
! local variables
real(wp_) :: isign, bsign
integer :: exp2pi, exp2piout
logical :: phiccw, psiincr, qpos, phiccwout, psiincrout, qposout
call decode_cocos(cocosin, exp2pi, phiccw, psiincr, qpos)
call decode_cocos(cocosout, exp2piout, phiccwout, psiincrout, qposout)
! Check sign consistency
isign = sign(one, data%psia)
if (.not.psiincr) isign = -isign
bsign = sign(one, data%fpol(size(data%fpol)))
if (qpos .neqv. isign * bsign * data%qpsi(size(data%qpsi)) > zero) then
! Warning: sign inconsistency found among q, Ipla and Bref
data%qpsi = -data%qpsi
if (present(error)) error = 1
2015-11-18 17:34:33 +01:00
else
if (present(error)) error = 0
2015-11-18 17:34:33 +01:00
end if
! Convert cocosin to cocosout
! Opposite direction of toroidal angle phi in cocosin and cocosout
if (phiccw .neqv. phiccwout) data%fpol = -data%fpol
! q has opposite sign for given sign of Bphi*Ip
if (qpos .neqv. qposout) data%qpsi = -data%qpsi
! psi and Ip signs don't change accordingly
if ((phiccw .eqv. phiccwout) .neqv. (psiincr .eqv. psiincrout)) &
data%psia = -data%psia
2015-11-18 17:34:33 +01:00
! Convert Wb to Wb/rad or viceversa
data%psia = data%psia * (2.0_wp_*pi)**(exp2piout - exp2pi)
2015-11-18 17:34:33 +01:00
end subroutine change_cocos
subroutine decode_cocos(cocos, exp2pi, phiccw, psiincr, qpos)
! Extracts the sign and units conventions from a COCOS index
2015-11-18 17:34:33 +01:00
implicit none
! subroutine arguments
integer, intent(in) :: cocos
2015-11-18 17:34:33 +01:00
integer, intent(out) :: exp2pi
logical, intent(out) :: phiccw, psiincr, qpos
! local variables
integer :: cmod10, cmod4
cmod10 = mod(cocos, 10)
cmod4 = mod(cmod10, 4)
! cocos>10 ψ in Wb, cocos<10 ψ in Wb/rad
exp2pi = (cocos - cmod10)/10
! cocos mod 10 = 1,3,5,7: toroidal angle φ increasing CCW
phiccw = (mod(cmod10, 2)== 1)
! cocos mod 10 = 1,2,5,6: ψ increasing with positive Ip
psiincr = (cmod4==1 .or. cmod4==2)
! cocos mod 10 = 1,2,7,8: q positive for positive Bφ*Ip
qpos = (cmod10<3 .or. cmod10>6)
2015-11-18 17:34:33 +01:00
end subroutine decode_cocos
subroutine scale_equil(params, data)
! Rescale the magnetic field (B) and the plasma current (I)
! and/or force their signs.
!
! Notes:
! 1. signi and signb are ignored on input if equal to 0.
! They are used to assign the direction of Bphi and Ipla BEFORE scaling.
! 2. cocos=3 assumed: CCW direction is >0
! 3. Bphi and Ipla scaled by the same factor factb to keep q unchanged
! 4. factb<0 reverses the directions of Bphi and Ipla
2015-11-18 17:34:33 +01:00
use const_and_precisions, only : one
use gray_params, only : equilibrium_parameters, equilibrium_data
2015-11-18 17:34:33 +01:00
implicit none
! subroutine arguments
type(equilibrium_parameters), intent(inout) :: params
type(equilibrium_data), intent(inout) :: data
! local variables
real(wp_) :: last_fpol
last_fpol = data%fpol(size(data%fpol))
2015-11-18 17:34:33 +01:00
if (params%sgni /=0) &
data%psia = sign(data%psia, real(-params%sgni, wp_))
if (params%sgnb /=0 .and. params%sgnb * last_fpol < 0) &
data%fpol = -data%fpol
data%psia = data%psia * params%factb
data%fpol = data%fpol * params%factb
params%sgni = int(sign(one, -data%psia))
params%sgnb = int(sign(one, last_fpol))
end subroutine scale_equil
2015-11-18 17:34:33 +01:00
subroutine set_equil_spline(params, data)
! Computes splines for the MHD equilibrium data and stores them
! in their respective global variables, see the top of this file.
use const_and_precisions, only : zero, one
use gray_params, only : equilibrium_parameters, equilibrium_data
use gray_params, only : iequil
use reflections, only : inside
use utils, only : vmaxmin, vmaxmini
2021-12-18 18:57:38 +01:00
use logger, only : log_info
2015-11-18 17:34:33 +01:00
implicit none
! subroutine arguments
type(equilibrium_parameters), intent(in) :: params
type(equilibrium_data), intent(in) :: data
! local variables
integer :: nr, nz, nrest, nzest, npsest, nrz, npsi, nbnd, ibinf, ibsup
real(wp_) :: tension, rax0, zax0, psinoptmp, psinxptmp
real(wp_) :: rbmin, rbmax, rbinf, rbsup, r1, z1
real(wp_), dimension(size(data%psinr)) :: rhotn
real(wp_), dimension(:), allocatable :: rv1d, zv1d, fvpsi, wf
integer :: ier, ixploc, info, i, j, ij
2021-12-18 18:57:38 +01:00
character(256) :: msg ! for log messages formatting
2015-11-18 17:34:33 +01:00
! compute array sizes
nr = size(data%rv)
nz = size(data%zv)
nrz = nr*nz
npsi = size(data%psinr)
nrest = nr + ksplp
nzest = nz + ksplp
npsest = npsi + ksplp
! length in m !!!
rmnm = data%rv(1)
rmxm = data%rv(nr)
zmnm = data%zv(1)
zmxm = data%zv(nz)
! Spline interpolation of ψ(R, z)
if (iequil>2) then
! data valid only inside boundary (data%psin=0 outside), e.g. source==ESCO
! presence of boundary anticipated here to filter invalid data
nbnd = min(size(data%rbnd), size(data%zbnd))
! allocate knots and spline coefficients arrays
if (allocated(psi_spline%knots_x)) deallocate(psi_spline%knots_x)
if (allocated(psi_spline%knots_y)) deallocate(psi_spline%knots_y)
if (allocated(psi_spline%coeffs)) deallocate(psi_spline%coeffs)
allocate(psi_spline%knots_x(nrest), psi_spline%knots_y(nzest))
allocate(psi_spline%coeffs(nrz))
! determine number of valid grid points
nrz=0
do j=1,nz
do i=1,nr
if (nbnd.gt.0) then
if(.not.inside(data%rbnd,data%zbnd,nbnd,data%rv(i),data%zv(j))) cycle
else
if(data%psin(i,j).le.0.0d0) cycle
end if
nrz=nrz+1
end do
end do
! store valid data
allocate(rv1d(nrz),zv1d(nrz),fvpsi(nrz),wf(nrz))
ij=0
do j=1,nz
do i=1,nr
if (nbnd.gt.0) then
if(.not.inside(data%rbnd,data%zbnd,nbnd,data%rv(i),data%zv(j))) cycle
else
if(data%psin(i,j).le.0.0d0) cycle
end if
ij=ij+1
rv1d(ij)=data%rv(i)
zv1d(ij)=data%zv(j)
fvpsi(ij)=data%psin(i,j)
wf(ij)=1.0d0
end do
end do
! Fit as a scattered set of points
! use reduced number of knots to limit memory comsumption ?
psi_spline%nknots_x=nr/4+4
psi_spline%nknots_y=nz/4+4
tension = params%ssplps
call scatterspl(rv1d, zv1d, fvpsi, wf, nrz, kspl, tension, &
rmnm, rmxm, zmnm, zmxm, &
psi_spline%knots_x, psi_spline%nknots_x, &
psi_spline%knots_y, psi_spline%nknots_y, &
psi_spline%coeffs, ier)
! if failed, re-fit with an interpolating spline (zero tension)
if(ier == -1) then
tension = 0
psi_spline%nknots_x=nr/4+4
psi_spline%nknots_y=nz/4+4
call scatterspl(rv1d, zv1d, fvpsi, wf, nrz, kspl, tension, &
rmnm, rmxm, zmnm, zmxm, &
psi_spline%knots_x, psi_spline%nknots_x, &
psi_spline%knots_y, psi_spline%nknots_y, &
psi_spline%coeffs, ier)
end if
deallocate(rv1d, zv1d, wf, fvpsi)
! reset nrz to the total number of grid points for next allocations
nrz = nr*nz
else
! iequil==2: data are valid on the full R,z grid
! reshape 2D ψ array to 1D (transposed)
allocate(fvpsi(nrz))
fvpsi = reshape(transpose(data%psin), [nrz])
! compute spline coefficients
call psi_spline%init(data%rv, data%zv, fvpsi, nr, nz, &
range=[rmnm, rmxm, zmnm, zmxm], &
tension=params%ssplps, err=ier)
! if failed, re-fit with an interpolating spline (zero tension)
if(ier == -1) then
call psi_spline%init(data%rv, data%zv, fvpsi, nr, nz, &
range=[rmnm, rmxm, zmnm, zmxm], &
tension=zero)
end if
deallocate(fvpsi)
2015-11-18 17:34:33 +01:00
end if
! compute spline coefficients for ψ(R,z) partial derivatives
call psi_spline%init_deriv(nr, nz, 1, 0) ! ∂ψ/∂R
call psi_spline%init_deriv(nr, nz, 0, 1) ! ∂ψ/∂z
call psi_spline%init_deriv(nr, nz, 1, 1) ! ∂²ψ/∂R∂z
call psi_spline%init_deriv(nr, nz, 2, 0) ! ∂²ψ/∂R²
call psi_spline%init_deriv(nr, nz, 0, 2) ! ∂²ψ/∂z²
! Spline interpolation of F(ψ)
! give a small weight to the last point
2015-11-18 17:34:33 +01:00
allocate(wf(npsi))
wf(1:npsi-1) = 1
wf(npsi) = 1.0e2_wp_
call fpol_spline%init(data%psinr, data%fpol, npsi, range=[zero, one], &
weights=wf, tension=params%ssplf)
deallocate(wf)
! set vacuum value used outside 0 ≤ ψ ≤ 1 range
fpolas = fpol_spline%eval(data%psinr(npsi))
sgnbphi = sign(one ,fpolas)
! Re-normalize ψ after spline computation
2015-11-18 17:34:33 +01:00
! Start with un-corrected psi
psia = data%psia
psinop = 0
psiant = 1
2015-11-18 17:34:33 +01:00
! Use provided boundary to set an initial guess
! for the search of O/X points
nbnd=min(size(data%rbnd), size(data%zbnd))
2015-11-18 17:34:33 +01:00
if (nbnd>0) then
call vmaxmini(data%zbnd,nbnd,zbinf,zbsup,ibinf,ibsup)
rbinf=data%rbnd(ibinf)
rbsup=data%rbnd(ibsup)
call vmaxmin(data%rbnd,nbnd,rbmin,rbmax)
2015-11-18 17:34:33 +01:00
else
zbinf=data%zv(2)
zbsup=data%zv(nz-1)
rbinf=data%rv((nr+1)/2)
2015-11-18 17:34:33 +01:00
rbsup=rbinf
rbmin=data%rv(2)
rbmax=data%rv(nr-1)
2015-11-18 17:34:33 +01:00
end if
! Search for exact location of the magnetic axis
rax0=data%rax
zax0=data%zax
2015-11-18 17:34:33 +01:00
call points_ox(rax0,zax0,rmaxis,zmaxis,psinoptmp,info)
2021-12-18 18:57:38 +01:00
write (msg, '("O-point found:", 3(x,a,"=",g0.3))') &
'r', rmaxis, 'z', zmaxis, 'ψ', psinoptmp
call log_info(msg, mod='equilibrium', proc='set_equil_spline')
2015-11-18 17:34:33 +01:00
! search for X-point if params%ixp /= 0
2015-11-18 17:34:33 +01:00
ixploc = params%ixp
2015-11-18 17:34:33 +01:00
if(ixploc/=0) then
if(ixploc<0) then
call points_ox(rbinf,zbinf,r1,z1,psinxptmp,info)
if(psinxptmp/=-1.0_wp_) then
2021-12-18 18:57:38 +01:00
write (msg, '("X-point found:", 3(x,a,"=",g0.3))') &
'r', r1, 'z', z1, 'ψ', psinxptmp
call log_info(msg, mod='equilibrium', proc='set_equil_spline')
2021-12-18 18:57:38 +01:00
2015-11-18 17:34:33 +01:00
zbinf=z1
psinop=psinoptmp
psiant=psinxptmp-psinop
call points_tgo(rmaxis,0.5_wp_*(zmaxis+zbsup),r1,z1,one,info)
zbsup=z1
else
ixploc=0
end if
else
call points_ox(rbsup,zbsup,r1,z1,psinxptmp,info)
if(psinxptmp.ne.-1.0_wp_) then
2021-12-18 18:57:38 +01:00
write (msg, '("X-point found:", 3(x,a,"=",g0.3))') &
'r', r1, 'z', z1, 'ψ', psinxptmp
call log_info(msg, mod='equilibrium', proc='set_equil_spline')
2021-12-18 18:57:38 +01:00
2015-11-18 17:34:33 +01:00
zbsup=z1
psinop=psinoptmp
psiant=psinxptmp-psinop
call points_tgo(rmaxis,0.5_wp_*(zmaxis+zbinf),r1,z1,one,info)
zbinf=z1
else
ixploc=0
end if
end if
end if
if (ixploc==0) then
psinop=psinoptmp
psiant=one-psinop
! Find upper horizontal tangent point
2015-11-18 17:34:33 +01:00
call points_tgo(rmaxis,0.5_wp_*(zmaxis+zbsup),r1,z1,one,info)
zbsup=z1
rbsup=r1
! Find lower horizontal tangent point
2015-11-18 17:34:33 +01:00
call points_tgo(rmaxis,0.5_wp_*(zmaxis+zbinf),r1,z1,one,info)
zbinf=z1
rbinf=r1
2021-12-18 18:57:38 +01:00
write (msg, '("X-point not found in", 2(x,a,"∈[",g0.3,",",g0.3,"]"))') &
'r', rbinf, rbsup, 'z', zbinf, zbsup
call log_info(msg, mod='equilibrium', proc='set_equil_spline')
2015-11-18 17:34:33 +01:00
end if
! Save Bt value on axis (required in flux_average and used in Jcd def)
! and vacuum value B0 at ref. radius data%rvac (used in Jcd_astra def)
2015-11-18 17:34:33 +01:00
call equinum_fpol(zero, btaxis)
btaxis = btaxis/rmaxis
btrcen = fpolas/data%rvac
rcen = data%rvac
2021-12-18 18:57:38 +01:00
write (msg, '(2(a,g0.3))') 'Bt_center=', btrcen, ' Bt_axis=', btaxis
call log_info(msg, mod='equilibrium', proc='set_equil_spline')
! Compute ρ_p/ρ_t mapping based on the input q profile
call setqphi_num(data%psinr, abs(data%qpsi), abs(psia), rhotn)
call set_rho_spline(sqrt(data%psinr), rhotn)
end subroutine set_equil_spline
2015-11-18 17:34:33 +01:00
subroutine scatterspl(x,y,z,w,n,kspl,sspl,xmin,xmax,ymin,ymax, &
tx,nknt_x,ty,nknt_y,coeff,ierr)
! Computes the spline interpolation of a surface when
! the data points are irregular, i.e. not on a uniform grid
use const_and_precisions, only : comp_eps
use dierckx, only : surfit
2015-11-18 17:34:33 +01:00
implicit none
! subroutine arguments
integer, intent(in) :: n
real(wp_), dimension(n), intent(in) :: x, y, z
real(wp_), dimension(n), intent(in) :: w
integer, intent(in) :: kspl
real(wp_), intent(in) :: sspl
real(wp_), intent(in) :: xmin, xmax, ymin, ymax
real(wp_), dimension(nknt_x), intent(inout) :: tx
real(wp_), dimension(nknt_y), intent(inout) :: ty
integer, intent(inout) :: nknt_x, nknt_y
real(wp_), dimension(nknt_x*nknt_y), intent(out) :: coeff
integer, intent(out) :: ierr
! local variables
integer :: iopt
real(wp_) :: resid
integer :: u,v,km,ne,b1,b2,lwrk1,lwrk2,kwrk,nxest,nyest
real(wp_), dimension(:), allocatable :: wrk1, wrk2
integer, dimension(:), allocatable :: iwrk
nxest=nknt_x
nyest=nknt_y
ne = max(nxest,nyest)
km = kspl+1
u = nxest-km
v = nyest-km
b1 = kspl*min(u,v)+kspl+1
b2 = (kspl+1)*min(u,v)+1
lwrk1 = u*v*(2+b1+b2)+2*(u+v+km*(n+ne)+ne-2*kspl)+b2+1
lwrk2 = u*v*(b2+1)+b2
kwrk = n+(nknt_x-2*kspl-1)*(nknt_y-2*kspl-1)
allocate(wrk1(lwrk1),wrk2(lwrk2),iwrk(kwrk))
iopt=0
call surfit(iopt,n,x,y,z,w,xmin,xmax,ymin,ymax,kspl,kspl, &
sspl,nxest,nyest,ne,comp_eps,nknt_x,tx,nknt_y,ty, &
coeff,resid,wrk1,lwrk1,wrk2,lwrk2,iwrk,kwrk,ierr)
deallocate(wrk1,wrk2,iwrk)
end subroutine scatterspl
subroutine setqphi_num(psinq,q,psia,rhotn)
! Computes the spline of the safety factor q(ψ)
use const_and_precisions, only : pi
implicit none
! subroutine arguments
real(wp_), dimension(:), intent(in) :: psinq,q
real(wp_), intent(in) :: psia
real(wp_), dimension(:), intent(out), optional :: rhotn
! local variables
real(wp_), dimension(size(q)) :: phit
real(wp_) :: dx
integer :: k
call q_spline%init(psinq, q, size(q))
! Toroidal flux φ = 2π ∫q(ψ)dψ
phit(1)=0
do k=1,q_spline%ndata-1
dx=q_spline%data(k+1)-q_spline%data(k)
phit(k+1)=phit(k) + dx*(q_spline%coeffs(k,1) + dx*(q_spline%coeffs(k,2)/2 + &
dx*(q_spline%coeffs(k,3)/3 + dx* q_spline%coeffs(k,4)/4) ) )
end do
phitedge=phit(q_spline%ndata)
if(present(rhotn)) rhotn(1:q_spline%ndata)=sqrt(phit/phitedge)
phitedge=2*pi*psia*phitedge
end subroutine setqphi_num
subroutine set_equil_an(data, n)
! Computes the analytical equilibrium data and stores them
! in their respective global variables, see the top of this file.
use const_and_precisions, only : pi, zero, one
use gray_params, only : equilibrium_data
implicit none
! subroutine arguments
type(equilibrium_data), intent(in) :: data
integer, optional, intent(in) :: n
! local variables
integer, parameter :: nqdef=101
integer :: i
real(wp_) :: dr, fq0, fq1, qq, res, rn
real(wp_) :: rax, zax, a, bax, qax, q1, qexp
real(wp_), dimension(:), allocatable :: rhotn,rhopn
rax = data%rv(1)
zax = data%zv(1)
a = data%rv(2)
bax = data%fpol(1) / rax
qax = data%qpsi(1)
q1 = data%qpsi(2)
qexp = data%qpsi(3)
btaxis = bax
rmaxis = rax
zmaxis = zax
btrcen = bax
rcen = rax
aminor = a
zbinf = zmaxis-a
zbsup = zmaxis+a
model%q0 = qax
model%q1 = q1
model%lq = qexp
sgnbphi = sign(one, bax)
rmxm = rmaxis+aminor
rmnm = rmaxis-aminor
zmxm = zbsup
zmnm = zbinf
if (present(n)) then
q_spline%ndata = n
else
q_spline%ndata = nqdef
end if
if (allocated(q_spline%data)) deallocate(q_spline%data)
allocate(q_spline%data(q_spline%ndata))
allocate(rhotn(q_spline%ndata), rhopn(q_spline%ndata))
dr = one/(q_spline%ndata - 1)
rhotn(1) = zero
q_spline%data(1) = zero
res = zero
fq0 = zero
do i = 2, q_spline%ndata
rn = (i - 1)*dr
qq = model%q0 + (model%q1 - model%q0) * rn**model%lq
fq1 = rn / qq
res = res + (fq1 + fq0)/2 * dr
fq0 = fq1
rhotn(i) = rn
q_spline%data(i) = res
end do
phitedge = btaxis*aminor**2 ! temporary
psia = res*phitedge
phitedge = pi*phitedge ! final
q_spline%data = q_spline%data/res
rhopn = sqrt(q_spline%data)
call set_rho_spline(rhopn, rhotn)
end subroutine set_equil_an
subroutine set_rho_spline(rhop, rhot)
! Computes the splines for converting between the poloidal (ρ_p)
! and toroidal (ρ_t) normalised radii
implicit none
! subroutine arguments
real(wp_), dimension(:), intent(in) :: rhop, rhot
call rhop_spline%init(rhot, rhop, size(rhop))
call rhot_spline%init(rhop, rhot, size(rhot))
end subroutine set_rho_spline
subroutine equinum_psi(R, z, psi, dpsidr, dpsidz, &
ddpsidrr, ddpsidzz, ddpsidrz)
! Computes the normalised poloidal flux ψ (and its derivatives)
! as a function of (R, z) using the numerical equilibrium
2015-11-18 17:34:33 +01:00
implicit none
! subroutine arguments
real(wp_), intent(in) :: R, z
real(wp_), intent(out), optional :: psi, dpsidr, dpsidz
real(wp_), intent(out), optional :: ddpsidrr, ddpsidzz, ddpsidrz
! Note: here lengths are measured in meters
if (R <= rmxm .and. R >= rmnm .and. &
z <= zmxm .and. z >= zmnm) then
if (present(psi)) psi = (psi_spline%eval(R, z) - psinop) / psiant
if (present(dpsidr)) dpsidr = psi_spline%deriv(R, z, 1, 0) * psia
if (present(dpsidz)) dpsidz = psi_spline%deriv(R, z, 0, 1) * psia
if (present(ddpsidrr)) ddpsidrr = psi_spline%deriv(R, z, 2, 0) * psia
if (present(ddpsidzz)) ddpsidzz = psi_spline%deriv(R, z, 0, 2) * psia
if (present(ddpsidrz)) ddpsidrz = psi_spline%deriv(R, z, 1, 1) * psia
2015-11-18 17:34:33 +01:00
else
if (present(psi)) psi = -1
if (present(dpsidr)) dpsidr = 0
if (present(dpsidz)) dpsidz = 0
if (present(ddpsidrr)) ddpsidrr = 0
if (present(ddpsidzz)) ddpsidzz = 0
if (present(ddpsidrz)) ddpsidrz = 0
2015-11-18 17:34:33 +01:00
end if
end subroutine equinum_psi
subroutine equinum_fpol(psi, fpol, dfpol)
! Computes the poloidal current function F(ψ)
! (and its derivative) using the numerical equilibrium
2015-11-18 17:34:33 +01:00
implicit none
! subroutine arguments
real(wp_), intent(in) :: psi
real(wp_), intent(out) :: fpol
real(wp_), intent(out), optional :: dfpol
if(psi <= 1 .and. psi >= 0) then
fpol = fpol_spline%eval(psi)
if (present(dfpol)) dfpol = fpol_spline%deriv(psi) / psia
2015-11-18 17:34:33 +01:00
else
fpol = fpolas
if (present(dfpol)) dfpol = 0
2015-11-18 17:34:33 +01:00
end if
end subroutine equinum_fpol
subroutine equian(R, z, psi, fpolv, dfpv, dpsidr, dpsidz, &
ddpsidrr, ddpsidzz, ddpsidrz)
! Computes all MHD equilibrium parameters from a simple analytical model
2015-11-18 17:34:33 +01:00
implicit none
! subroutine arguments
real(wp_), intent(in) :: R,z
real(wp_), intent(out), optional :: psi, fpolv, dfpv, dpsidr, dpsidz, &
ddpsidrr, ddpsidzz, ddpsidrz
! local variables
real(wp_) :: cst, dpsidrp, d2psidrp, dqq, qq, &
rn, rpm, snt, rhop, rhot, btaxqq
! simple model for equilibrium: large aspect ratio
! outside plasma: analytical continuation, not solution Maxwell eqs
! Note: rpm is ρ_t in meters, rn is ρ_t normalised
rpm = hypot(R - rmaxis, z - zmaxis)
rn = rpm/aminor
snt = 0
cst = 1
if (rpm > 0) then
snt = (z - zmaxis)/rpm
cst = (R - rmaxis)/rpm
end if
if (present(psi)) then
rhot=rn
if(rn <= 1) then
rhop = frhopol(rhot)
psi = rhop**2
else
psi = 1 + btaxis / (2*psia*model%q1) * (rpm**2 - aminor**2)
rhop = sqrt(psi)
2015-11-18 17:34:33 +01:00
end if
end if
if(rn <= 1) then
qq = model%q0 + (model%q1 - model%q0) * rn**model%lq
btaxqq = btaxis/qq
dpsidrp = btaxqq*rpm
dqq = model%lq * (model%q1 - model%q0) * rn**(model%lq - 1)
d2psidrp=btaxqq*(1 - rn*dqq/qq)
else
btaxqq = btaxis / model%q1
dpsidrp = btaxqq * rpm
d2psidrp = btaxqq
end if
if(present(fpolv)) fpolv = btaxis * rmaxis
if(present(dfpv)) dfpv = 0
if(present(dpsidr)) dpsidr = dpsidrp*cst
if(present(dpsidz)) dpsidz = dpsidrp*snt
if(present(ddpsidrr)) ddpsidrr = btaxqq*snt**2 + cst**2*d2psidrp
if(present(ddpsidrz)) ddpsidrz = cst * snt * (d2psidrp - btaxqq)
if(present(ddpsidzz)) ddpsidzz = btaxqq * cst**2 + snt**2 * d2psidrp
end subroutine equian
function frhopol(rhot)
! Converts from toroidal (ρ_t) to poloidal (ρ_p) normalised radius
2015-11-18 17:34:33 +01:00
implicit none
! function arguments
real(wp_), intent(in) :: rhot
real(wp_) :: frhopol
2015-11-18 17:34:33 +01:00
frhopol = rhop_spline%eval(rhot)
end function frhopol
2015-11-18 17:34:33 +01:00
function frhopolv(rhot)
! Vector variant of `frhopol`
use utils, only : locate
implicit none
! function arguments
real(wp_), intent(in) :: rhot(:)
real(wp_) :: frhopolv(size(rhot))
! local variables
integer :: i, i0, j
i0 = 1
do j = 1, size(rhot)
call locate(rhop_spline%xdata(i0:), rhop_spline%ndata-i0+1, rhot(j), i)
i = min(max(1,i), rhop_spline%ndata - i0) + i0 - 1
frhopolv(j) = rhop_spline%raw_eval(i, rhot(j))
i0 = i
2015-11-18 17:34:33 +01:00
end do
end function frhopolv
2015-11-18 17:34:33 +01:00
function frhotor(rhop)
! Converts from poloidal (ρ_p) to toroidal (ρ_t) normalised radius
2015-11-18 17:34:33 +01:00
implicit none
! function arguments
real(wp_), intent(in) :: rhop
real(wp_) :: frhotor
frhotor = rhot_spline%eval(rhop)
end function frhotor
2015-11-18 17:34:33 +01:00
function fq(psin)
! Computes the safety factor q as a function of the
! normalised poloidal flux ψ
2015-11-18 17:34:33 +01:00
use gray_params, only : iequil
2015-11-18 17:34:33 +01:00
implicit none
! function arguments
2015-11-18 17:34:33 +01:00
real(wp_), intent(in) :: psin
real(wp_) :: fq
! local variables
real(wp_) :: rn
if (iequil < 2) then
! q = q₀ + (q₁ - q₀)ρ^l
rn = frhotor(sqrt(psin))
fq = model%q0 + (model%q1 - model%q0) * rn**model%lq
2015-11-18 17:34:33 +01:00
else
fq = q_spline%eval(psin)
2015-11-18 17:34:33 +01:00
end if
end function fq
subroutine bfield(rpsim, zpsim, bphi, br, bz)
! Computes the magnetic field as a function of
! (R, z) in cylindrical coordinates
use gray_params, only : iequil
2015-11-18 17:34:33 +01:00
implicit none
! subroutine arguments
real(wp_), intent(in) :: rpsim, zpsim
real(wp_), intent(out), optional :: bphi,br,bz
! local variables
real(wp_) :: psin,fpol
2015-11-18 17:34:33 +01:00
if (iequil < 2) then
call equian(rpsim,zpsim,fpolv=bphi,dpsidr=bz,dpsidz=br)
if (present(bphi)) bphi=bphi/rpsim
else
call equinum_psi(rpsim,zpsim,psi=bphi,dpsidr=bz,dpsidz=br)
if (present(bphi)) then
psin=bphi
call equinum_fpol(psin,fpol)
bphi=fpol/rpsim
end if
end if
if (present(br)) br=-br/rpsim
if (present(bz)) bz= bz/rpsim
end subroutine bfield
2015-11-18 17:34:33 +01:00
function tor_curr(rpsim,zpsim) result(jphi)
! Computes the toroidal current Jφ as a function of (R, z)
use const_and_precisions, only : ccj=>mu0inv
use gray_params, only : iequil
2015-11-18 17:34:33 +01:00
implicit none
! function arguments
real(wp_), intent(in) :: rpsim, zpsim
real(wp_) :: jphi
! local variables
real(wp_) :: bzz,dbvcdc13,dbvcdc31
real(wp_) :: dpsidr,ddpsidrr,ddpsidzz
if(iequil < 2) then
call equian(rpsim,zpsim,dpsidr=dpsidr, &
ddpsidrr=ddpsidrr,ddpsidzz=ddpsidzz)
else
call equinum_psi(rpsim,zpsim,dpsidr=dpsidr, &
ddpsidrr=ddpsidrr,ddpsidzz=ddpsidzz)
end if
bzz= dpsidr/rpsim
dbvcdc13=-ddpsidzz/rpsim
dbvcdc31= ddpsidrr/rpsim-bzz/rpsim
jphi=ccj*(dbvcdc13-dbvcdc31)
end function tor_curr
function tor_curr_psi(psin) result(jphi)
! Computes the toroidal current Jφ as a function of ψ
implicit none
2015-11-18 17:34:33 +01:00
! function arguments
real(wp_), intent(in) :: psin
real(wp_) :: jphi
! local variables
real(wp_) :: r1, r2
call psi_raxis(psin, r1, r2)
jphi = tor_curr(r2, zmaxis)
end function tor_curr_psi
2015-11-18 17:34:33 +01:00
subroutine psi_raxis(psin,r1,r2)
2021-12-18 18:57:38 +01:00
use gray_params, only : iequil
use dierckx, only : profil, sproota
use logger, only : log_error
2015-11-18 17:34:33 +01:00
implicit none
2021-12-18 18:57:38 +01:00
! subroutine arguments
real(wp_) :: psin,r1,r2
2021-12-18 18:57:38 +01:00
! local constants
integer, parameter :: mest=4
! local variables
integer :: iopt,ier,m
real(wp_) :: zc,val
real(wp_), dimension(mest) :: zeroc
real(wp_), dimension(psi_spline%nknots_x) :: czc
2021-12-18 18:57:38 +01:00
character(64) :: msg
2015-11-18 17:34:33 +01:00
if (iequil < 2) then
val=frhotor(sqrt(psin))
r1=rmaxis-val*aminor
r2=rmaxis+val*aminor
else
iopt=1
zc=zmaxis
call profil(iopt, psi_spline%knots_x, psi_spline%nknots_x, &
psi_spline%knots_y, psi_spline%nknots_y, &
psi_spline%coeffs, kspl, kspl, zc, &
psi_spline%nknots_x, czc, ier)
2021-12-18 18:57:38 +01:00
if (ier > 0) then
write (msg, '("profil failed with error ",g0)') ier
call log_error(msg, mod='equilibrium', proc='psi_raxis')
end if
val=psin*psiant+psinop
call sproota(val, psi_spline%knots_x, psi_spline%nknots_x, &
czc, zeroc, mest, m, ier)
r1=zeroc(1)
r2=zeroc(2)
end if
end subroutine psi_raxis
2015-11-18 17:34:33 +01:00
2015-11-18 17:34:33 +01:00
subroutine points_ox(rz,zz,rf,zf,psinvf,info)
! Finds the location of the O,X points
2015-11-18 17:34:33 +01:00
use const_and_precisions, only : comp_eps
2021-12-18 18:57:38 +01:00
use minpack, only : hybrj1
use logger, only : log_error, log_debug
2015-11-18 17:34:33 +01:00
implicit none
2021-12-18 18:57:38 +01:00
! local constants
2015-11-18 17:34:33 +01:00
integer, parameter :: n=2,ldfjac=n,lwa=(n*(n+13))/2
2021-12-18 18:57:38 +01:00
! arguments
2015-11-18 17:34:33 +01:00
real(wp_), intent(in) :: rz,zz
real(wp_), intent(out) :: rf,zf,psinvf
integer, intent(out) :: info
2021-12-18 18:57:38 +01:00
! local variables
2015-11-18 17:34:33 +01:00
real(wp_) :: tol
real(wp_), dimension(n) :: xvec,fvec
real(wp_), dimension(lwa) :: wa
real(wp_), dimension(ldfjac,n) :: fjac
2021-12-18 18:57:38 +01:00
character(256) :: msg
2015-11-18 17:34:33 +01:00
xvec(1)=rz
xvec(2)=zz
tol = sqrt(comp_eps)
call hybrj1(fcnox,n,xvec,fvec,fjac,ldfjac,tol,info,wa,lwa)
2021-12-18 18:57:38 +01:00
if(info /= 1) then
write (msg, '("O,X coordinates:",2(x,", ",g0.3))') xvec
call log_debug(msg, mod='equilibrium', proc='points_ox')
write (msg, '("hybrj1 failed with error ",g0)') info
call log_error(msg, mod='equilibrium', proc='points_ox')
end if
rf=xvec(1)
zf=xvec(2)
call equinum_psi(rf,zf,psinvf)
2015-11-18 17:34:33 +01:00
end subroutine points_ox
2015-11-18 17:34:33 +01:00
subroutine fcnox(n,x,fvec,fjac,ldfjac,iflag)
2021-12-18 18:57:38 +01:00
use logger, only : log_error
2015-11-18 17:34:33 +01:00
implicit none
2021-12-18 18:57:38 +01:00
! subroutine arguments
2015-11-18 17:34:33 +01:00
integer, intent(in) :: n,iflag,ldfjac
real(wp_), dimension(n), intent(in) :: x
real(wp_), dimension(n), intent(inout) :: fvec
real(wp_), dimension(ldfjac,n), intent(inout) :: fjac
2021-12-18 18:57:38 +01:00
! local variables
2015-11-18 17:34:33 +01:00
real(wp_) :: dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz
2021-12-18 18:57:38 +01:00
character(64) :: msg
2015-11-18 17:34:33 +01:00
select case(iflag)
case(1)
call equinum_psi(x(1),x(2),dpsidr=dpsidr,dpsidz=dpsidz)
fvec(1) = dpsidr/psia
fvec(2) = dpsidz/psia
case(2)
call equinum_psi(x(1),x(2),ddpsidrr=ddpsidrr,ddpsidzz=ddpsidzz, &
ddpsidrz=ddpsidrz)
fjac(1,1) = ddpsidrr/psia
fjac(1,2) = ddpsidrz/psia
fjac(2,1) = ddpsidrz/psia
fjac(2,2) = ddpsidzz/psia
case default
2021-12-18 18:57:38 +01:00
write (msg, '("invalid iflag: ",g0)')
call log_error(msg, mod='equilibrium', proc='fcnox')
2015-11-18 17:34:33 +01:00
end select
end subroutine fcnox
2015-11-18 17:34:33 +01:00
subroutine points_tgo(rz,zz,rf,zf,psin0,info)
use const_and_precisions, only : comp_eps
2021-12-18 18:57:38 +01:00
use minpack, only : hybrj1mv
use logger, only : log_error, log_debug
2015-11-18 17:34:33 +01:00
implicit none
2021-12-18 18:57:38 +01:00
! local constants
2015-11-18 17:34:33 +01:00
integer, parameter :: n=2,ldfjac=n,lwa=(n*(n+13))/2
2021-12-18 18:57:38 +01:00
! arguments
2015-11-18 17:34:33 +01:00
real(wp_), intent(in) :: rz,zz,psin0
real(wp_), intent(out) :: rf,zf
integer, intent(out) :: info
2021-12-18 18:57:38 +01:00
character(256) :: msg
! local variables
2015-11-18 17:34:33 +01:00
real(wp_) :: tol
real(wp_), dimension(n) :: xvec,fvec,f0
real(wp_), dimension(lwa) :: wa
real(wp_), dimension(ldfjac,n) :: fjac
xvec(1)=rz
xvec(2)=zz
f0(1)=psin0
f0(2)=0.0_wp_
tol = sqrt(comp_eps)
call hybrj1mv(fcntgo,n,xvec,f0,fvec,fjac,ldfjac,tol,info,wa,lwa)
2021-12-18 18:57:38 +01:00
if(info /= 1) then
write (msg, '("R,z coordinates:",5(x,g0.3))') xvec, rz, zz, psin0
call log_debug(msg, mod='equilibrium', proc='points_tgo')
write (msg, '("hybrj1mv failed with error ",g0)') info
call log_error(msg, mod='equilibrium', proc='points_tgo')
2015-11-18 17:34:33 +01:00
end if
rf=xvec(1)
zf=xvec(2)
end
2015-11-18 17:34:33 +01:00
subroutine fcntgo(n,x,f0,fvec,fjac,ldfjac,iflag)
use logger, only : log_error
2021-12-18 18:57:38 +01:00
2015-11-18 17:34:33 +01:00
implicit none
2021-12-18 18:57:38 +01:00
! subroutine arguments
2015-11-18 17:34:33 +01:00
integer, intent(in) :: n,ldfjac,iflag
real(wp_), dimension(n), intent(in) :: x,f0
real(wp_), dimension(n), intent(inout) :: fvec
real(wp_), dimension(ldfjac,n), intent(inout) :: fjac
2021-12-18 18:57:38 +01:00
! local variables
2015-11-18 17:34:33 +01:00
real(wp_) :: psinv,dpsidr,dpsidz,ddpsidrr,ddpsidrz
2021-12-18 18:57:38 +01:00
character(64) :: msg
2015-11-18 17:34:33 +01:00
select case(iflag)
case(1)
call equinum_psi(x(1),x(2),psinv,dpsidr)
fvec(1) = psinv-f0(1)
fvec(2) = dpsidr/psia-f0(2)
case(2)
call equinum_psi(x(1),x(2),dpsidr=dpsidr,dpsidz=dpsidz, &
ddpsidrr=ddpsidrr,ddpsidrz=ddpsidrz)
fjac(1,1) = dpsidr/psia
fjac(1,2) = dpsidz/psia
fjac(2,1) = ddpsidrr/psia
fjac(2,2) = ddpsidrz/psia
case default
2021-12-18 18:57:38 +01:00
write (msg, '("invalid iflag: ",g0)')
call log_error(msg, mod='equilibrium', proc='fcntgo')
2015-11-18 17:34:33 +01:00
end select
end subroutine fcntgo
subroutine unset_equil_spline
! Unsets the splines global variables, see the top of this file.
2015-11-18 17:34:33 +01:00
implicit none
call fpol_spline%deinit
call psi_spline%deinit
call q_spline%deinit
call rhop_spline%deinit
call rhot_spline%deinit
end subroutine unset_equil_spline
2015-11-18 17:34:33 +01:00
end module equilibrium