istp/gray
2
0
Fork 0
Code Issues Pull Requests Packages Projects Releases 1 Wiki Activity

src/magsurf_data.f90: cleanup

Browse Source
This commit is contained in:
Michele Guerini Rocco 2024-10-17 00:56:02 +02:00 committed by rnhmjoj
parent 2e2ab16273
commit fdf5ef72fe
Signed by: rnhmjoj
GPG Key ID: BFBAF4C975F76450
7 changed files with 408 additions and 469 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

17
src/eccd.f90
View File

@ -139,21 +139,24 @@ contains
end subroutine setcdcoeff_cohen
subroutine setcdcoeff_ncl(zeff,rbx,fc,amu,rhop,cst2,eccdpar)
use magsurf_data, only : ch,tjp,tlm,njpt,nlmt
use magsurf_data, only : cd_eff
use dierckx, only : profil
use logger, only : log_warning
integer, parameter :: ksp=3
real(wp_), intent(in) :: zeff,rbx,fc,amu,rhop
real(wp_), intent(out) :: cst2
real(wp_), dimension(:), allocatable, intent(out) :: eccdpar
real(wp_), dimension(nlmt) :: chlm
integer :: nlm,ierr,npar
real(wp_), dimension(cd_eff%nknots_y) :: chlm
integer :: nlm, ierr, npar
cst2=1.0_wp_-rbx
if(cst2<cst2min) cst2=0.0_wp_
call SpitzFuncCoeff(amu,Zeff,fc)
nlm=nlmt
call profil(0,tjp,njpt,tlm,nlmt,ch,ksp,ksp,rhop,nlm,chlm,ierr)
nlm = cd_eff%nknots_y
call profil(0, cd_eff%knots_x, cd_eff%nknots_x, &
cd_eff%knots_y, cd_eff%nknots_y, &
cd_eff%coeffs, 3, 3, rhop, &
cd_eff%nknots_y, chlm, ierr)
if (ierr > 0) then
block
character(256) :: msg
@ -167,7 +170,7 @@ contains
eccdpar(1)=zeff
eccdpar(2) = fc
eccdpar(3) = rbx
eccdpar(4:3+nlm) = tlm
eccdpar(4:3+nlm) = cd_eff%knots_y
eccdpar(4+nlm:npar) = chlm
end subroutine setcdcoeff_ncl

5
src/gray_core.f90
View File

@ -21,7 +21,7 @@ contains
print_err_raytracing, print_err_ecrh_cd
use beams, only : xgygcoeff, launchangles2n
use beamdata, only : pweight
use magsurf_data, only : compute_flux_averages, dealloc_surfvec
use magsurf_data, only : compute_flux_averages
use pec, only : pec_init, spec, postproc_profiles, dealloc_pec, &
rhop_tab, rhot_tab
use utils, only : vmaxmin
@ -660,7 +660,6 @@ contains
! ============ main loop END ============
! Free memory
call dealloc_surfvec
call dealloc_pec
end block
@ -1799,7 +1798,7 @@ contains
error = error + ierrcd
if (allocated(eccdpar)) deallocate(eccdpar)
! current drive efficiency <j/p> [A⋅m/W]
! current drive efficiency R* = <J_∥>/<dP/dV> [A⋅m/W]
effjcdav = rbavi*effjcd
dIdp = equil%sgn_bphi * effjcdav / (2*pi*rrii)

41
src/gray_equil.f90
View File

@ -176,19 +176,22 @@ module gray_equil
contains
pure subroutine b_field(self, R, z, B_R, B_z, B_phi)
! Computes the magnetic field as a function of
! (R, z) in cylindrical coordinates
pure subroutine b_field(self, R, z, phi, B_R, B_z, B_phi, B, B_tot)
! Computes the magnetic field as a function of (R, z)
! The outputs include cylindrical coordinates, cartesian vector and norm.
!
! Note: all output arguments are optional.
! subroutine arguments
class(abstract_equil), intent(in) :: self
real(wp_), intent(in) :: R, z
real(wp_), intent(in), optional :: phi
real(wp_), intent(out), optional :: B_R, B_z, B_phi
real(wp_), intent(out), optional :: B(3), B_tot
! local variables
real(wp_) :: psi_n, fpol, dpsidr, dpsidz
real(wp_) :: B_R_, B_z_, B_phi_
call self%pol_flux(R, z, psi_n, dpsidr, dpsidz)
call self%pol_curr(psi_n, fpol)
@ -204,10 +207,34 @@ contains
! and carrying out the cross products:
!
! B = F(ψ)∇φ - ∂ψ/∂z ∇R/R + ∂ψ/∂R ∇z/R
! = F(ψ)φ^/R - ∂ψ/∂z R^/R + ∂ψ/∂R z^/R
!
if (present(B_R)) B_R = - 1/R * dpsidz * self%psi_a
if (present(B_z)) B_z = + 1/R * dpsidr * self%psi_a
if (present(B_phi)) B_phi = fpol / R
B_R_ = - 1/R * dpsidz * self%psi_a
B_z_ = + 1/R * dpsidr * self%psi_a
B_phi_ = fpol / R
if (present(B_R)) B_R = B_R_
if (present(B_z)) B_z = B_z_
if (present(B_phi)) B_phi = B_phi_
! Convert to Cartesian coordinates
!
! Since the unit vectors are given by:
!
! { R^ = x^ cosφ + y^ sinφ
! { φ^ = -x^ sinφ + y^ cosφ,
!
! B = φ^ B_φ + R^ B_R + z^ B_z
! = (B_R cosφ - B_φ sinφ) x^ + (B_R sinφ + B_φ cosφ) y^ + z^ B_Z
! = x^ B_x + y^ B_y + B_Z z^
!
if (present(B)) then
B(1) = B_R_*cos(phi) - B_phi_*sin(phi)
B(2) = B_R_*sin(phi) + B_phi_*cos(phi)
B(3) = B_z_
end if
if (present(B_tot)) B_tot = norm2([B_R_, B_z_, B_phi_])
end subroutine b_field
@ -571,7 +598,7 @@ contains
if (present(dfpol)) dfpol = 0
end if
end subroutine numeric_pol_curr
pure function numeric_safety(self, psi_n) result(q)
! function arguments

72
src/gray_tables.f90
View File

@ -116,7 +116,6 @@ contains
use gray_params, only : gray_parameters, EQ_VACUUM
use gray_equil, only : abstract_equil
use gray_plasma, only : abstract_plasma
use magsurf_data, only : npsi, vajphiav
! function arguments
type(gray_parameters), intent(in) :: params
@ -125,7 +124,8 @@ contains
type(table) :: tbl
! local variables
integer :: i, ntail
integer, parameter :: n_grid = 100
integer :: i, n_tail
real(wp_) :: rho_t, J_phi, dens, ddens
real(wp_) :: psi_n, rho_p, drho_p
@ -137,22 +137,17 @@ contains
if (params%equilibrium%iequil == EQ_VACUUM) return
! Δρ_p for the uniform ρ_p grid
! Note: we don't use ψ_n because J_phi has been
! sampled uniformly in ρ_p (see flux_average)
drho_p = 1.0_wp_ / (npsi - 1)
drho_p = 1.0_wp_ / (n_grid - 1)
! extra points to reach ψ=ψ_bnd
ntail = ceiling((sqrt(params%profiles%psnbnd) - 1) / drho_p)
n_tail = ceiling((sqrt(params%profiles%psnbnd) - 1) / drho_p)
do i = 0, npsi + ntail
do i = 0, n_grid + n_tail
rho_p = i * drho_p
rho_t = equil%pol2tor(rho_p)
psi_n = rho_p**2
if (psi_n < 1) then
J_phi = vajphiav(i+1) * 1.e-6_wp_
else
J_phi = 0
end if
!call equil%flux_average(rho_p, J_phi=J_phi)
J_phi = 0
call plasma%density(psi_n, dens, ddens)
call tbl%append([psi_n, rho_t, dens, plasma%temp(psi_n), &
equil%safety(psi_n), J_phi])
@ -164,10 +159,9 @@ contains
function ec_resonance(params, equil, B_res) result(tbl)
! Generates the EC resonance surface table
use const_and_precisions, only : comp_eps
use const_and_precisions, only : comp_eps, comp_huge
use gray_params, only : gray_parameters, EQ_VACUUM
use gray_equil, only : abstract_equil
use magsurf_data, only : npsi
use types, only : item
use cniteq, only : cniteq_f
@ -185,30 +179,30 @@ contains
block
! local variables
integer, parameter :: n_grid=101, n_cont=2002
integer :: j, k, n
integer :: n_conts, n_points(10), n_tot
real(wp_) :: dr, dz, B_min, B_max
real(wp_) :: B, B_R, B_z, B_phi, B_tot(npsi, npsi)
real(wp_), dimension(2002) :: R_cont, z_cont
real(wp_), dimension(npsi) :: R, z
real(wp_) :: B, B_R, B_z, B_phi, B_tot(n_grid, n_grid)
real(wp_), dimension(n_cont) :: R_cont, z_cont
real(wp_), dimension(n_grid) :: R, z
! Build a regular (R, z) grid
dr = (equil%r_range(2) - equil%r_range(1) - comp_eps)/(npsi - 1)
dz = (equil%z_range(2) - equil%z_range(1))/(npsi - 1)
do j=1,npsi
dR = (equil%r_range(2) - equil%r_range(1) - comp_eps)/(n_grid - 1)
dz = (equil%z_range(2) - equil%z_range(1))/(n_grid - 1)
do j = 1, n_grid
R(j) = comp_eps + equil%r_range(1) + dr*(j - 1)
z(j) = equil%z_range(1) + dz*(j - 1)
end do
! B_tot on psi grid
B_max = -1.0e30_wp_
B_min = +1.0e30_wp_
do k = 1, npsi
do j = 1, npsi
call equil%b_field(R(j), z(k), B_R, B_z, B_phi)
B_tot(j,k) = sqrt(B_R**2 + B_z**2 + B_phi**2)
if(B_tot(j,k) >= B_max) B_max = B_tot(j,k)
if(B_tot(j,k) <= B_min) B_min = B_tot(j,k)
! Compute magnetic field on the (R, z) grid
B_max = -comp_huge
B_min = +comp_huge
do k = 1, n_grid
do j = 1, n_grid
call equil%b_field(R(j), z(k), B_tot=B_tot(j,k))
if (B_tot(j,k) >= B_max) B_max = B_tot(j,k)
if (B_tot(j,k) <= B_min) B_min = B_tot(j,k)
end do
end do
@ -237,7 +231,6 @@ contains
use gray_params, only : gray_parameters, EQ_VACUUM
use gray_equil, only : abstract_equil
use gray_plasma, only : abstract_plasma
use magsurf_data, only : npsi
! function arguments
type(gray_parameters), intent(in) :: params
@ -248,11 +241,12 @@ contains
type(table) :: tbl
! local variables
integer, parameter :: n_grid = 100
integer :: j, k
real(wp_) :: R0, Npl0
real(wp_) :: dR, dz, B_R, B_phi, B_z, B
real(wp_) :: psi_n, ne, dne, Te, X, Y, Npl
real(wp_), dimension(npsi) :: R, z
real(wp_), dimension(n_grid) :: R, z
call tbl%init(title='inputs-maps', id=72, active=is_active(params, 72), &
labels=[character(64) :: 'R', 'z', 'ψ_n', 'B_R', 'B_z', 'B', &
@ -265,20 +259,19 @@ contains
Npl0 = sin(params%antenna%beta*degree) ! initial value of N∥
! Build a regular (R, z) grid
dR = (equil%r_range(2) - equil%r_range(1) - comp_eps)/(npsi - 1)
dz = (equil%z_range(2) - equil%z_range(1))/(npsi - 1)
do concurrent (j = 1:npsi)
dR = (equil%r_range(2) - equil%r_range(1) - comp_eps)/(n_grid - 1)
dz = (equil%z_range(2) - equil%z_range(1))/(n_grid - 1)
do concurrent (j = 1:n_grid)
R(j) = comp_eps + equil%r_range(1) + dR*(j - 1)
z(j) = equil%z_range(1) + dz*(j - 1)
end do
do j = 1, npsi
do j = 1, n_grid
Npl = Npl0 * R0/r(j)
do k = 1, npsi
do k = 1, n_grid
call equil%pol_flux(r(j), z(k), psi_n)
call equil%b_field(r(j), z(k), B_R, B_z, B_phi)
call equil%b_field(r(j), z(k), B_tot=B)
call plasma%density(psi_n, ne, dne)
B = sqrt(B_R**2 + B_phi**2 + B_z**2)
X = X_norm*ne
Y = B/B_res
Te = plasma%temp(psi_n)
@ -296,7 +289,6 @@ contains
use gray_params, only : gray_parameters
use gray_equil, only : abstract_equil
use magsurf_data, only : npoints
use logger, only : log_info
use minpack, only : hybrj1
use types, only : table
@ -336,7 +328,7 @@ contains
! Compute and print the ψ_n(R,z) = ψ_n₀ contour
block
real(wp_), dimension(npoints) :: R_cont, z_cont
real(wp_), dimension(101) :: R_cont, z_cont
real(wp_) :: R_hi, R_lo, z_hi, z_lo
! Guesses for the high,low horizontal-tangent points

698
src/magsurf_data.f90
View File

@ -1,448 +1,372 @@
module magsurf_data
use const_and_precisions, only : wp_
use splines, only : spline_simple
use splines, only : spline_simple, spline_2d
implicit none
integer, save :: npsi, npoints !# sup mag, # punti per sup
integer, save :: njpt, nlmt
real(wp_), save :: rarea
real(wp_), dimension(:), allocatable, save :: psicon,pstab,rhot_eq, &
rhotqv,bav,varea,vcurrp,vajphiav,qqv,ffc,vratja,vratjb
real(wp_), dimension(:), allocatable, save :: rpstab
real(wp_), dimension(:), allocatable, save :: vvol,rri,rbav,bmxpsi,bmnpsi
real(wp_), dimension(:), allocatable, save :: tjp,tlm,ch,ch01
real(wp_), dimension(:,:), allocatable, save :: rcon,zcon
type(spline_simple), save :: cvol, crri, crbav, cbmx, cbmn, carea, cfc
type(spline_simple), save :: crhotq
type(spline_simple), save :: cratja, cratjb, cratjpl
type(spline_simple), save :: cdadrhot, cdvdrhot
type(spline_simple), save :: spline_volume, spline_R_J, spline_B_avg_min
type(spline_simple), save :: spline_B_max, spline_B_min, spline_area, spline_f_c
type(spline_simple), save :: spline_J_astra, spline_J_jintrac, spline_J_paral
type(spline_simple), save :: spline_dAdrho_t, spline_dVdrho_t
type(spline_2d), save :: cd_eff
contains
subroutine alloc_cnt(ierr)
integer, intent(out) :: ierr
if(npsi.le.0.or.npoints.le.0) then
ierr = -1
return
end if
call dealloc_cnt
allocate(psicon(npsi),rcon(npoints,npsi),zcon(npoints,npsi))
end subroutine alloc_cnt
subroutine dealloc_cnt
if(allocated(psicon)) deallocate(psicon)
if(allocated(rcon)) deallocate(rcon)
if(allocated(zcon)) deallocate(zcon)
end subroutine dealloc_cnt
subroutine alloc_surfvec(ierr)
integer, intent(out) :: ierr
if(npsi.le.0.or.npoints.le.0) then
ierr = -1
return
end if
call dealloc_surfvec
call alloc_cnt(ierr)
allocate(pstab(npsi), &
rhot_eq(npsi),rhotqv(npsi),bav(npsi),bmxpsi(npsi),bmnpsi(npsi),varea(npsi), &
vvol(npsi),vcurrp(npsi),vajphiav(npsi),qqv(npsi),ffc(npsi),vratja(npsi), &
vratjb(npsi),rpstab(npsi),rri(npsi),rbav(npsi))
end subroutine alloc_surfvec
subroutine dealloc_surfvec
call dealloc_cnt
if(allocated(pstab)) deallocate(pstab)
if(allocated(rhot_eq)) deallocate(rhot_eq)
if(allocated(rhotqv)) deallocate(rhotqv)
if(allocated(bav)) deallocate(bav)
if(allocated(bmxpsi)) deallocate(bmxpsi)
if(allocated(bmnpsi)) deallocate(bmnpsi)
if(allocated(varea)) deallocate(varea)
if(allocated(vvol)) deallocate(vvol)
if(allocated(vcurrp)) deallocate(vcurrp)
if(allocated(vajphiav)) deallocate(vajphiav)
if(allocated(qqv)) deallocate(qqv)
if(allocated(ffc)) deallocate(ffc)
if(allocated(vratja)) deallocate(vratja)
if(allocated(vratjb)) deallocate(vratjb)
if(allocated(rpstab)) deallocate(rpstab)
if(allocated(rri)) deallocate(rri)
if(allocated(rbav)) deallocate(rbav)
if(allocated(tjp)) deallocate(tjp,tlm,ch)
call cvol%deinit
call crbav%deinit
call crri%deinit
call cbmx%deinit
call cbmn%deinit
call cratja%deinit
call cratjb%deinit
call cratjpl%deinit
call carea%deinit
call cfc%deinit
call cdadrhot%deinit
call cdvdrhot%deinit
end subroutine dealloc_surfvec
pure function cross(v, w)
! z-component of the cross product of 2D vectors
real(wp_), intent(in) :: v(2), w(2)
real(wp_) :: cross
cross = v(1)*w(2) - v(2)*w(1)
end function cross
subroutine compute_flux_averages(params, equil, tables)
use const_and_precisions, only : wp_,zero,one,pi,ccj=>mu0inv
use dierckx, only : regrid, coeff_parder
use types, only : table, wrap
use gray_params, only : gray_parameters, gray_tables
use gray_equil, only : abstract_equil
! Computes the splines of the flux surface averages
! subroutine arguments
use const_and_precisions, only : zero, one, comp_huge, pi, mu0inv
use types, only : table, wrap
use gray_params, only : gray_parameters, gray_tables
use gray_tables, only : store_flux_surface
use gray_equil, only : abstract_equil
! subroutine arguments
type(gray_parameters), intent(in) :: params
class(abstract_equil), intent(in) :: equil
type(gray_tables), intent(inout), optional :: tables
! local constants
integer, parameter :: nnintp=101,ncnt=100,nlam=101,ksp=3, &
njest=nnintp+ksp+1,nlest=nlam+ksp+1, &
lwrk=4*(nnintp+nlam)+11*(njest+nlest)+njest*nnintp+nlest+54, &
kwrk=nnintp+nlam+njest+nlest+3
integer, parameter :: nsurf=101, npoints=101, nlam=101
! local variables
integer :: ier,ierr,l,jp,inc,inc1,iopt,njp,nlm,ninpr
integer, dimension(kwrk) :: iwrk
real(wp_) :: lam,dlam,anorm,b2av,dvdpsi,dadpsi,ratio_cdator, &
ratio_cdbtor,ratio_pltor,fc,height,r2iav,currp, &
area,volume,ajphiav,bbav,bmmx,bmmn,btot0,bpoloid0,rpsim0,dla,dlb, &
dlp,drc,ph,area2,rzp,rz,rpsim,zpsim,btot,bpoloid,dlph,ajphi0, &
shlam,srl,rl2,rl0,rl,dhlam,dhlam0,ccfh,s,ajphi, &
bphi,brr,bzz,riav,fp,psinjp,rhopjp,rhotjp,qq,rup,rlw,zup,zlw
real(wp_), dimension(nnintp) :: dadrhotv,dvdrhotv,vratjpl
real(wp_), dimension(2*ncnt) :: dlpv
real(wp_), dimension(2*ncnt+1) :: bv,bpv
real(wp_), dimension(nlam) :: alam,weights
real(wp_), dimension(nnintp,nlam) :: fhlam
real(wp_), dimension(nnintp*nlam) :: ffhlam,dffhlam
real(wp_), dimension(lwrk) :: wrk
real(wp_), dimension(:), allocatable :: rctemp,zctemp
real(wp_) :: fpolv, ddpsidrr, ddpsidzz
integer :: i, j, l
npsi=nnintp
npoints = 2*ncnt+1
! Arrays of quantities evaluated on a uniform ρ_p grid
real(wp_) :: rho_p(nsurf) ! ρ_p ∈ [0, 1)
real(wp_) :: psi_n(nsurf) ! ψ_n = ρ_p²
real(wp_) :: B_avg(nsurf) ! ⟨B⟩ = 1/N ∮ B dl/B_p, average magnetic field
real(wp_) :: J_phi_avg(nsurf) ! ⟨J_φ⟩ = 1/N ∮ J_φ dl/B_p, average toroidal current
real(wp_) :: I_p(nsurf) ! I_p = 1/μ₀ ∫ B_p⋅dS_φ, plasma current
real(wp_) :: q(nsurf) ! q = 1/2π ∮ B_φ/B_p dl/R, safety factor
real(wp_) :: R_J(nsurf) ! R_J = ⟨B⟩ / (F(ψ) ⋅ ⟨1/R²⟩), effective J_cd radius
real(wp_) :: area(nsurf) ! poloidal area enclosed by the flux surface
real(wp_) :: volume(nsurf) ! volume V enclosed by the flux surface
real(wp_) :: dAdrho_t(nsurf) ! dA/dρ_t, where A is the poloidal area
real(wp_) :: dVdrho_t(nsurf) ! dV/dρ_t, where V is the enclosed volume
real(wp_) :: B_min(nsurf) ! min |B| on the flux surface
real(wp_) :: B_max(nsurf) ! max |B| on the flux surface
real(wp_) :: f_c(nsurf) ! fraction of circulating particles
call alloc_surfvec(ierr)
if(allocated(tjp)) deallocate(tjp)
if(allocated(tlm)) deallocate(tlm)
if(allocated(ch)) deallocate(ch)
allocate(tjp(njest),tlm(nlest),ch((njest-4)*(nlest-4)), &
rctemp(npoints),zctemp(npoints),stat=ierr)
if (ierr.ne.0) return
! Ratios of different J_cd definitions w.r.t. ⟨J_φ⟩
real(wp_) :: ratio_astra_tor(nsurf) ! J_cd_astra = ⟨J_cd⋅B⟩/B₀
real(wp_) :: ratio_jintrac_tor(nsurf) ! J_cd_jintrac = ⟨J_cd⋅B⟩/⟨B⟩
real(wp_) :: ratio_paral_tor(nsurf) ! Jcd_∥ = ⟨J_cd⋅B/|B|⟩
! computation of flux surface averaged quantities
! Arrays of quantities evaluated on the poloidal flux contour
real(wp_) :: R(npoints), z(npoints) ! R, z coordinates
real(wp_) :: B_tots(npoints) ! magnetic field
real(wp_) :: B_pols(npoints) ! poloidal field
real(wp_) :: dls(npoints) ! line element
dlam=1.0_wp_/dble(nlam-1)
do l=1,nlam-1
alam(l)=dble(l-1)*dlam
fhlam(1,l)=sqrt(1.0_wp_-alam(l))
ffhlam(l)=fhlam(1,l)
dffhlam(l)=-0.5_wp_/sqrt(1.0_wp_-alam(l))
weights(l)=1.0_wp_
! Intermediate results
real(wp_) :: norm ! N = ∮ dl/B_p, normalisation of the ⟨⋅⟩ integral
real(wp_) :: B2_avg ! ⟨B²⟩ = 1/N ∮ B² dl/B_p, average squared field
real(wp_) :: dAdpsi ! dA/dψ, where A is the poloidal area
real(wp_) :: dVdpsi ! dV/dψ, where V is the enclosed volume
real(wp_) :: R2_inv_avg ! ⟨R⁻²⟩ = 1/N ∮ R⁻² dl/B_p
real(wp_) :: R_inv_avg ! ⟨R⁻¹⟩ = 1/N ∮ R⁻¹ dl/B_p
real(wp_) :: fpol ! F(ψ), poloidal current function
! Neoclassical CD efficiency variables
real(wp_) :: dlam, lambda(nlam) ! uniform grid λ ∈ [0,1]
real(wp_) :: H(nlam, nsurf) ! H(ρ_p, λ) sampled on a uniform ρ_p×λ grid
! Miscellanea
real(wp_) :: B_R, B_z, B_phi, R_hi, R_lo, z_hi, z_lo
! On the magnetic axis, ρ_p = 0
psi_n(1) = 0
rho_p(1) = 0
area(1) = 0
volume(1) = 0
I_p(1) = 0
J_phi_avg(1) = equil%tor_curr(equil%axis(1), equil%axis(2))
B_max(1) = abs(equil%b_axis)
B_min(1) = abs(equil%b_axis)
B_avg(1) = abs(equil%b_axis)
R_J(1) = equil%axis(1)
f_c(1) = 1
dAdrho_t(1) = 0
dVdrho_t(1) = 0
ratio_paral_tor(1) = 1
ratio_astra_tor(1) = abs(equil%b_axis / equil%b_centre)
ratio_jintrac_tor(1) = 1
! The safety factor at ρ=0 is given by
! q₀ = B₀ / √(∂²ψ/∂R² ⋅ ∂²ψ/∂z²)
! Source: Spheromaks, Bellan
block
real(wp_) :: a, b
call equil%pol_flux(equil%axis(1), equil%axis(2), ddpsidrr=a, ddpsidzz=b)
q(1) = abs(equil%b_axis) / sqrt(a*b)/equil%psi_a
end block
! Compute the uniform λ grid and H(0, λ)
H = 0
dlam = one/(nlam - 1)
do l = 1, nlam
lambda(l) = (l-1) * dlam ! λ
H(l,1) = sqrt(1 - lambda(l)) ! H(0, λ) = √(1-λ)
end do
weights(1)=0.5_wp_
weights(nlam)=0.5_wp_
alam(nlam)=1.0_wp_
fhlam(1,nlam)=0.0_wp_
ffhlam(nlam)=0.0_wp_
dffhlam(nlam)=-99999.0_wp_
jp=1
anorm=2.0_wp_*pi*equil%axis(1)/abs(equil%b_axis)
dvdpsi=2.0_wp_*pi*anorm
dadpsi=2.0_wp_*pi/abs(equil%b_axis)
b2av=equil%b_axis**2
ratio_cdator=abs(equil%b_axis/equil%b_centre)
ratio_cdbtor=1.0_wp_
ratio_pltor=1.0_wp_
fc=1.0_wp_
call equil%pol_flux(equil%axis(1), equil%axis(2), &
ddpsidrr=ddpsidrr, ddpsidzz=ddpsidzz)
qq=abs(equil%b_axis)/sqrt(ddpsidrr*ddpsidzz)/equil%psi_a
ajphiav=-ccj*(ddpsidrr+ddpsidzz)*equil%psi_a/equil%axis(1)
psicon(1)=0.0_wp_
rcon(:,1)=equil%axis(1)
zcon(:,1)=equil%axis(2)
pstab(1)=0.0_wp_
rpstab(1)=0.0_wp_
vcurrp(1)=0.0_wp_
vajphiav(1)=ajphiav
bmxpsi(1)=abs(equil%b_axis)
bmnpsi(1)=abs(equil%b_axis)
bav(1)=abs(equil%b_axis)
rbav(1)=1.0_wp_
rri(1)=equil%axis(1)
varea(1)=0.0_wp_
vvol(1)=0.0_wp_
vratjpl(1)=ratio_pltor
vratja(1)=ratio_cdator
vratjb(1)=ratio_cdbtor
ffc(1)=fc
qqv(1)=qq
dadrhotv(1)=0.0_wp_
dvdrhotv(1)=0.0_wp_
! Guess for first contour reconstruction
R_hi = equil%axis(1)
R_lo = equil%axis(1)
z_hi = equil%axis(2) + (equil%z_boundary(2) - equil%axis(2))/10
z_lo = equil%axis(2) - (equil%axis(2) - equil%z_boundary(1))/10
rup=equil%axis(1)
rlw=equil%axis(1)
zup=equil%axis(2)+(equil%z_boundary(2)-equil%axis(2))/10.0_wp_
zlw=equil%axis(2)-(equil%axis(2)-equil%z_boundary(1))/10.0_wp_
! For each surface with ρ_p > 0
surfaces_loop: do i = 2, nsurf
! Note: keep ρ_p < 1 to avoid the X point
rho_p(i) = merge((i - 1) * one/(nsurf - 1), 0.9999_wp_, i < nsurf)
psi_n(i) = rho_p(i)**2
do jp=2,npsi
height=dble(jp-1)/dble(npsi-1)
if(jp.eq.npsi) height=0.9999_wp_
rhopjp=height
psinjp=height*height
rhotjp=equil%pol2tor(rhopjp)
if (rhotjp /= rhotjp) print *, 'Nan!!!!!!!!!!'
psicon(jp)=height
call equil%flux_contour(psi_n(i), params%misc%rwall, &
R, z, R_hi, z_hi, R_lo, z_lo)
call equil%flux_contour(psinjp, params%misc%rwall, &
rctemp, zctemp, rup, zup, rlw, zlw)
rcon(:,jp) = rctemp
zcon(:,jp) = zctemp
if ((mod(i, 10) == 0 .or. i == nsurf) .and. present(tables)) then
call store_flux_surface(tables%flux_surfaces, psi_n(i), R, z)
end if
r2iav=0.0_wp_
anorm=0.0_wp_
dadpsi=0.0_wp_
currp=0.0_wp_
b2av=0.0_wp_
area=0.0_wp_
volume=0.0_wp_
ajphiav=0.0_wp_
bbav=0.0_wp_
bmmx=-1.0e+30_wp_
bmmn=1.0e+30_wp_
block
! Variables for the previous and current contour point
real(wp_) :: B_tot0, B_pol0, B_tot1, B_pol1, J_phi0, J_phi1
real(wp_) :: dl ! line element
ajphi0 = equil%tor_curr(rctemp(1), zctemp(1))
call equil%b_field(rctemp(1), zctemp(1), brr, bzz, bphi)
fpolv=bphi*rctemp(1)
btot0=sqrt(bphi**2+brr**2+bzz**2)
bpoloid0=sqrt(brr**2+bzz**2)
bv(1)=btot0
bpv(1)=bpoloid0
rpsim0=rctemp(1)
! Initialise all integrals
area(i) = 0
volume(i) = 0
norm = 0
dAdpsi = 0
I_p(i) = 0
B_avg(i) = 0
B2_avg = 0
R2_inv_avg = 0
J_phi_avg(i) = 0
B_max(i) = -comp_huge
B_min(i) = +comp_huge
do inc=1,npoints-1
inc1=inc+1
dla=sqrt((rctemp(inc)-equil%axis(1))**2+(zctemp(inc)-equil%axis(2))**2)
dlb=sqrt((rctemp(inc1)-equil%axis(1))**2+(zctemp(inc1)-equil%axis(2))**2)
dlp=sqrt((rctemp(inc1)-rctemp(inc))**2+(zctemp(inc1)-zctemp(inc))**2)
drc=(rctemp(inc1)-rctemp(inc))
! Evaluate integrands at the first "previous" point
J_phi0 = equil%tor_curr(R(1), z(1))
call equil%b_field(R(1), z(1), B_R=B_R, B_z=B_z, &
B_phi=B_phi, B_tot=B_tot0)
fpol = B_phi*R(1)
B_pol0 = sqrt(B_R**2 + B_z**2)
B_tots(1) = B_tot0
B_pols(1) = B_pol0
! compute length, area and volume defined by psi=psinjp=height^2
ph=0.5_wp_*(dla+dlb+dlp)
area2=ph*(ph-dla)*(ph-dlb)*(ph-dlp)
area=area+sqrt(area2)
rzp=rctemp(inc1)*zctemp(inc1)
rz=rctemp(inc)*zctemp(inc)
volume=pi*(rzp+rz)*drc+volume
contour_loop: do j = 1, npoints-1
block
! Compute the area by decomposing the contour into triangles
! formed by the magnetic axis and two adjacent points, whose
! area is given by the cross-product of the vertices, ½|v×w|.
real(wp_) :: O(2), P(2), Q(2), tri_area
O = equil%axis
P = [R(j), z(j)]
Q = [R(j+1), z(j+1)]
tri_area = abs(cross(P-O, Q-O))/2
area(i) = area(i) + tri_area
! compute line integrals on the contour psi=psinjp=height^2
rpsim=rctemp(inc1)
zpsim=zctemp(inc1)
call equil%b_field(rpsim, zpsim, brr, bzz)
ajphi = equil%tor_curr(rpsim, zpsim)
bphi=fpolv/rpsim
btot=sqrt(bphi**2+brr**2+bzz**2)
bpoloid=sqrt(brr**2+bzz**2)
dlpv(inc)=dlp
bv(inc1)=btot
bpv(inc1)=bpoloid
dlph=0.5_wp_*dlp
anorm=anorm+dlph*(1.0_wp_/bpoloid+1.0_wp_/bpoloid0)
dadpsi=dadpsi+dlph*(1.0_wp_/(bpoloid*rpsim)+1.0_wp_/(bpoloid0*rpsim0))
currp=currp+dlph*(bpoloid+bpoloid0)
b2av=b2av+dlph*(btot0**2/bpoloid0+btot**2/bpoloid)
bbav=bbav+dlph*(btot/bpoloid+btot0/bpoloid0)
r2iav=r2iav+dlph*(1.0_wp_/(bpoloid*rpsim**2)+1.0_wp_/(bpoloid0*rpsim0**2))
ajphiav=ajphiav+dlph*(ajphi0/(bpoloid0*rpsim0)+ajphi/(bpoloid*rpsim))
ajphi0=ajphi
rpsim0=rpsim
bpoloid0=bpoloid
btot0=btot
! Compute the volume using the centroid theorem (V = 2πR⋅A) with
! the centroid R of the contour as the mean of the triangles
! centroids weighted by their area. Note: the total area cancels.
volume(i) = volume(i) + 2*pi * (O(1) + P(1)+ Q(1))/3 * tri_area
! computation maximum/minimum B values on given flux surface
if(btot.le.bmmn) bmmn=btot
if(btot.ge.bmmx) bmmx=btot
end do
! The line differential is difference of the two vertices
dl = norm2(P - Q)
end block
! bav=<B> [T] , b2av=<B^2> [T^2] , rbav=<B>/b_min
! anorm = int d l_p/B_p = dV/dpsi/(2pi)
! r2iav=<1/R^2> [m^-2] ,
! riav=<1/R> [m^-1] = dA/dpsi/(dV/dpsi/(2pi)),
! rri = <B>/(|R B_tor|<1/R^2>) , used to compute I_tor [m^-1]
! currp = plasma current within psi=const
! Evaluate integrands at the current point for line integrals ∮ dl
call equil%b_field(R(j+1), z(j+1), B_R=B_R, B_z=B_z, B_tot=B_tot1)
J_phi1 = equil%tor_curr(R(j+1), z(j+1))
B_pol1 = sqrt(B_R**2 + B_z**2)
dls(j) = dl
B_tots(j+1) = B_tot1
B_pols(j+1) = B_pol1
bbav=bbav/anorm
r2iav=r2iav/anorm
dvdpsi=2.0_wp_*pi*anorm
riav=dadpsi/anorm
b2av=b2av/anorm
vcurrp(jp)=ccj*currp
vajphiav(jp)=ajphiav/dadpsi
! Do one step of the trapezoid method
! N = ∮ dl/B_p
! dA/dψ = ∮ (1/R) dl/B_p
! I_p = 1/μ₀ ∫ B_p⋅dS_φ = 1/μ₀ ∮ B_p dl
! ⟨B⟩ = 1/N ∮ B dl/B_p
! ⟨B²⟩ = 1/N ∮ B² dl/B_p
! ⟨R⁻²⟩ = 1/N ∮ R⁻² dl/B_p
! ⟨J_φ⟩ = 1/N ∮ J_φ dl/B_p
norm = norm + dl*(1/B_pol1 + 1/B_pol0)/2
dAdpsi = dAdpsi + dl*(1/(B_pol1*R(j+1)) + 1/(B_pol0*R(j)))/2
I_p(i) = I_p(i) + dl*(B_pol1 + B_pol0)/2 * mu0inv
B_avg(i) = B_avg(i) + dl*(B_tot1/B_pol1 + B_tot0/B_pol0)/2
B2_avg = B2_avg + dl*(B_tot1**2/B_pol1 + B_tot0**2/B_pol0)/2
R2_inv_avg = R2_inv_avg + dl*(1/(B_pol1*R(j+1)**2) + 1/(B_pol0*R(j)**2))/2
J_phi_avg(i) = J_phi_avg(i) + dl*(J_phi0/(B_pol0*R(j)) + J_phi1/(B_pol1*R(j+1)))/2
! area == varea, volume == vvol
! flux surface minor radius == (area/pi)^1/2
! ratio_cdator = Jcd_astra/J_phi Jcd_astra = <Jcd.B>/B0
! ratio_cdbtor = Jcd_jintrac/J_phi Jcd_jintrac = <Jcd.B>/<B>
! ratio_pltor = Jcd_||/J_phi Jcd_|| = <Jcd.b>
pstab(jp)=psinjp
rpstab(jp)=rhopjp
vvol(jp)=abs(volume)
varea(jp)=area
bav(jp)=bbav
rbav(jp)=bbav/bmmn
bmxpsi(jp)=bmmx
bmnpsi(jp)=bmmn
rri(jp)=bav(jp)/abs(fpolv*r2iav)
ratio_cdator=abs(b2av*riav/(fpolv*r2iav*equil%b_centre))
ratio_cdbtor=abs(b2av*riav/(fpolv*r2iav*bbav))
ratio_pltor=abs(bbav*riav/(fpolv*r2iav))
vratjpl(jp)=ratio_pltor
vratja(jp)=ratio_cdator
vratjb(jp)=ratio_cdbtor
qq=abs(dvdpsi*fpolv*r2iav/(4.0_wp_*pi*pi))
qqv(jp)=qq
dadrhotv(jp)=equil%phi_a*rhotjp/equil%safety(psinjp)*dadpsi/pi
dvdrhotv(jp)=equil%phi_a*rhotjp/equil%safety(psinjp)*dvdpsi/pi
! computation of fraction of circulating/trapped fraction fc, ft
! and of function H(lambda,rhop)
! ffhlam = Bmn/Bmx/fc integral_lambda^1 dlam/<sqrt(1-lam*B(rhop)/Bmx)>
fc=0.0_wp_
shlam=0.0_wp_
do l=nlam,1,-1
lam=alam(l)
srl=0.0_wp_
rl2=1.0_wp_-lam*bv(1)/bmmx
rl0=0.0_wp_
if(rl2.gt.0) rl0=sqrt(rl2)
do inc=1,npoints-1
rl2=1.0_wp_-lam*bv(inc+1)/bmmx
rl=0.0_wp_
if(rl2.gt.0) rl=sqrt(rl2)
srl=srl+0.5_wp_*dlpv(inc)*(rl/bpv(inc+1)+rl0/bpv(inc))
rl0=rl
! Update the previous values
J_phi0 = J_phi1
B_pol0 = B_pol1
B_tot0 = B_tot1
! Compute max,min magnetic field
if(B_tot1 <= B_min(i)) B_min(i) = B_tot1
if(B_tot1 >= B_max(i)) B_max(i) = B_tot1
end do contour_loop
end block
! Normalise integrals to obtain average value
B_avg(i) = B_avg(i) / norm
B2_avg = B2_avg / norm
R2_inv_avg = R2_inv_avg / norm
R_inv_avg = dAdpsi / norm
J_phi_avg(i) = J_phi_avg(i)/dAdpsi
dVdpsi = 2*pi * norm
! J_cd_astra/⟨J_φ⟩ where J_cd_astra = ⟨J_cd⋅B⟩/B₀
! J_cd_jintrac/⟨J_φ⟩ where J_cd_jintrac = ⟨J_cd⋅B⟩/⟨B⟩
! J_cd_∥/⟨J_φ⟩ where Jcd_∥ = ⟨J_cd⋅B/|B|⟩
ratio_astra_tor(i) = abs(B2_avg *R_inv_avg/(fpol*R2_inv_avg*equil%b_centre))
ratio_jintrac_tor(i) = abs(B2_avg *R_inv_avg/(fpol*R2_inv_avg*B_avg(i)))
ratio_paral_tor(i) = abs(B_avg(i)*R_inv_avg/(fpol*R2_inv_avg))
R_J(i) = B_avg(i)/abs(fpol*R2_inv_avg)
! By definition q(γ) = 1/2π ∮_γ dφ, where γ is a field line:
!
! { γ̅(0) = γ̅₀
! { dγ̅/dt = B̅(γ(t))
!
! Using the poloidal angle θ, we rewrite it as:
!
! q(γ) = 1/2π ∮_γ dφ/dθ dθ = 1/2π ∮_γ dφ/dt / dθ/dt dθ
!
! By definition of directional derivative, d/dt = B̅⋅∇, so
!
! q(γ) = 1/2π ∮_γ B̅⋅∇φ / B̅⋅∇θ dθ = 1/2π ∮_γ B_φ/B_p ρ/R dθ
!
! Finally, since ρdθ=dl in the poloidal plane and B_φ=F/R:
!
! q(ρ) = F/2π ∮ 1/R² dl/B_p = F/2π N ⟨R⁻²⟩
!
q(i) = abs(fpol/(2*pi) * norm * R2_inv_avg)
! Using Φ=Φ_a⋅ρ_t² and q=1/2π⋅∂Φ/∂ψ (see gray_equil.f90), we have
!
! ∂ψ/∂ρ_t = ∂ψ/∂Φ ∂Φ/∂ρ_t = 1/(2πq) Φ_a 2ρ_t
! = Φ_a ρ_t / (πq)
block
real(wp_) :: rho_t, dpsidrho_t
rho_t = equil%pol2tor(rho_p(i))
dpsidrho_t = equil%phi_a * rho_t / (equil%safety(psi_n(i)) * pi)
dAdrho_t(i) = dAdpsi * dpsidrho_t
dVdrho_t(i) = dVdpsi * dpsidrho_t
end block
! Compute the fraction of circulating particles:
!
! f_c(ρ_p) = ¾ <B²/B_max²> ∫₀¹ dλ λ/⟨√[1-λ'⋅B(ρ_p)/B_max]⟩
!
! and the function:
!
! H(ρ_p, λ) = ½ B_max/B_min/f_c ∫_λ^1 dλ'/⟨√[1-λ'⋅B(ρ_p)/B_max]⟩
!
f_c(i) = 0
block
real(wp_) :: srl, rl0, rl1, dHdlam1, dHdlam0
do l = nlam, 1, -1
! Compute srl = ⟨√[1-λ B(ρ_p)/B_max]⟩ with the trapezoid method
srl = 0
rl0 = sqrt(max(1 - lambda(l)*B_tots(1)/B_max(i), zero))
do j = 1, npoints-1
rl1 = sqrt(max(1 - lambda(l)*B_tots(j+1)/B_max(i), zero))
srl = srl + dls(j) * (rl1/B_pols(j+1) + rl0/B_pols(j))/2
rl0 = rl1
end do
srl = srl / norm
! Do one step of the Riemann sum of f_c
f_c(i) = f_c(i) + 3/4.0_wp_ * B2_avg/B_max(i)**2 &
* dlam * lambda(l)/srl * merge(1.0_wp_, 0.5_wp_, l > 1 .and. l < nlam)
! Do one step of the trapezoid method for ∫_λ^1 dλ'/srl
dHdlam1 = 1 / srl
if (l /= nlam) H(l, i) = H(l+1, i) + dlam*(dHdlam1 + dHdlam0)/2
dHdlam0 = dHdlam1
end do
srl=srl/anorm
dhlam=0.5_wp_/srl
fc=fc+lam/srl*weights(l)
if(l.eq.nlam) then
fhlam(jp,l)=0.0_wp_
ffhlam(nlam*(jp-1)+l)=0.0_wp_
dffhlam(nlam*(jp-1)+l)=-dhlam
dhlam0=dhlam
else
shlam=shlam+0.5_wp_*(dhlam+dhlam0)*dlam
fhlam(jp,l)=shlam
dffhlam(nlam*(jp-1)+l)=-dhlam
dhlam0=dhlam
end if
end do
fc=0.75_wp_*b2av/bmmx**2*fc*dlam
ffc(jp)=fc
ccfh=bmmn/bmmx/fc
do l=1,nlam
ffhlam(nlam*(jp-1)+l)=ccfh*fhlam(jp,l)
dffhlam(nlam*(jp-1)+l)=ccfh*dffhlam(nlam*(jp-1)+l)
end do
end do
rpstab(npsi)=1.0_wp_
pstab(npsi)=1.0_wp_
end block
! spline coefficients of area,vol,rbav,rri,bmxpsi,bmnpsi,fc,dadrhot,dvdrhot,ratioJs
! used for computations of dP/dV and J_cd
call cvol%init(rpstab, vvol, npsi)
call crbav%init(rpstab, rbav, npsi)
call crri%init(rpstab, rri, npsi)
call cbmx%init(rpstab, bmxpsi, npsi)
call cbmn%init(rpstab, bmnpsi, npsi)
call cratja%init(rpstab, vratja, npsi)
call cratjb%init(rpstab, vratjb, npsi)
call cratjpl%init(rpstab, vratjpl, npsi)
call carea%init(rpstab, varea, npsi)
call cfc%init(rpstab, ffc, npsi)
call cdadrhot%init(rpstab, dadrhotv, npsi)
call cdvdrhot%init(rpstab, dvdrhotv, npsi)
! Add leading factors
H(:, i) = H(:, i)/2 * B_min(i)/B_max(i)/f_c(i)
! spline interpolation of H(lambda,rhop) and dH/dlambda
iopt=0
s=0.0_wp_
call regrid(iopt,npsi,rpstab,nlam,alam,ffhlam,zero,one,zero,one, &
ksp,ksp,s,njest,nlest,njp,tjp,nlm,tlm,ch,fp, &
wrk,lwrk,iwrk,kwrk,ier)
njpt=njp
nlmt=nlm
end do surfaces_loop
deallocate(rctemp, zctemp)
! Fake ρ_p = 1 on the last surface (actually slightly < 1)
rho_p(nsurf) = 1
psi_n(nsurf) = 1
! Compute splines
call spline_volume%init(rho_p, volume, nsurf)
call spline_B_avg_min%init(rho_p, B_avg/B_min, nsurf)
call spline_R_J%init(rho_p, R_J, nsurf)
call spline_B_max%init(rho_p, B_max, nsurf)
call spline_B_min%init(rho_p, B_min, nsurf)
call spline_J_astra%init(rho_p, ratio_astra_tor, nsurf)
call spline_J_jintrac%init(rho_p, ratio_jintrac_tor, nsurf)
call spline_J_paral%init(rho_p, ratio_paral_tor, nsurf)
call spline_area%init(rho_p, area, nsurf)
call spline_f_c%init(rho_p, f_c, nsurf)
call spline_dAdrho_t%init(rho_p, dAdrho_t, nsurf)
call spline_dVdrho_t%init(rho_p, dVdrho_t, nsurf)
call cd_eff%init(rho_p, lambda, reshape(H, [nsurf*nlam]), &
nsurf, nlam, range=[zero, one, zero, one])
if (.not. present(tables)) return
do jp=1,npsi
if (.not. tables%flux_averages%active) exit
if (.not. tables%flux_averages%active) return
do i = 1, nsurf
call tables%flux_averages%append([ &
rpstab(jp), equil%pol2tor(rpstab(jp)), bav(jp), bmxpsi(jp), &
bmnpsi(jp), varea(jp), vvol(jp), vcurrp(jp), vajphiav(jp), &
ffc(jp), vratja(jp), vratjb(jp), qqv(jp)])
end do
ninpr = (npsi-1)/10
do jp = ninpr+1, npsi, ninpr
! Note: can't use store_flux_surface to avoid a circular dependency
if (.not. tables%flux_surfaces%active) exit
do l = 1, npoints
call tables%flux_surfaces%append( &
[wrap(l), wrap(psicon(jp)), wrap(rcon(l,jp)), wrap(zcon(l,jp))])
end do
rho_p(i), equil%pol2tor(rho_p(i)), B_avg(i), B_max(i), &
B_min(i), area(i), volume(i), I_p(i), J_phi_avg(i), &
f_c(i), ratio_astra_tor(i), ratio_jintrac_tor(i), q(i)])
end do
end subroutine compute_flux_averages
subroutine fluxval(rhop,area,vol,dervol,dadrhot,dvdrhot, &
rri,rbav,bmn,bmx,fc,ratja,ratjb,ratjpl)
subroutine fluxval(rho_p, area, vol, dervol, dadrhot, dvdrhot, &
rri, rbav, bmn, bmx, fc, ratja, ratjb, ratjpl)
use const_and_precisions, only : wp_
! subroutine arguments
real(wp_), intent(in) :: rhop
real(wp_), intent(in) :: rho_p
real(wp_), intent(out), optional :: &
vol, area, rri, rbav, dervol, bmn, bmx, fc, &
ratja, ratjb, ratjpl, dadrhot, dvdrhot
if (present(area)) area = carea%eval(rhop)
if (present(vol)) vol = cvol%eval(rhop)
if (present(area)) area = spline_area%eval(rho_p)
if (present(vol)) vol = spline_volume%eval(rho_p)
if (present(dervol)) dervol = cvol%deriv(rhop)
if (present(dadrhot)) dadrhot = cdadrhot%eval(rhop)
if (present(dvdrhot)) dvdrhot = cdvdrhot%eval(rhop)
if (present(dervol)) dervol = spline_volume%deriv(rho_p)
if (present(dadrhot)) dAdrhot = spline_dAdrho_t%eval(rho_p)
if (present(dvdrhot)) dVdrhot = spline_dVdrho_t%eval(rho_p)
if (present(rri)) rri = crri%eval(rhop)
if (present(rbav)) rbav = crbav%eval(rhop)
if (present(bmn)) bmn = cbmn%eval(rhop)
if (present(bmx)) bmx = cbmx%eval(rhop)
if (present(fc)) fc = cfc%eval(rhop)
if (present(ratja)) ratja = cratja%eval(rhop)
if (present(ratjb)) ratjb = cratjb%eval(rhop)
if (present(ratjpl)) ratjpl = cratjpl%eval(rhop)
if (present(rri)) rri = spline_R_J%eval(rho_p)
if (present(rbav)) rbav = spline_B_avg_min%eval(rho_p)
if (present(bmn)) bmn = spline_B_min%eval(rho_p)
if (present(bmx)) bmx = spline_B_max%eval(rho_p)
if (present(fc)) fc = spline_f_c%eval(rho_p)
if (present(ratja)) ratja = spline_J_astra%eval(rho_p)
if (present(ratjb)) ratjb = spline_J_jintrac%eval(rho_p)
if (present(ratjpl)) ratjpl = spline_J_paral%eval(rho_p)
end subroutine fluxval
end module magsurf_data

3
src/main.f90
View File

@ -238,7 +238,7 @@ contains
! (of different beams) and recomputes the summary table
use gray_params, only : gray_parameters
use gray_equil, only : abstract_equil
use magsurf_data, only : compute_flux_averages, dealloc_surfvec
use magsurf_data, only : compute_flux_averages
use pec, only : pec_init, postproc_profiles, dealloc_pec, &
rhot_tab
@ -383,7 +383,6 @@ contains
end do
! Free memory
call dealloc_surfvec
call dealloc_pec
end subroutine sum_profiles

41
src/multipass.f90
View File

@ -26,7 +26,7 @@ contains
logical, intent(in) :: perfect ! whether to assume perfect coupling
! local variables
real(wp_) :: R, z, cosphi, sinphi, B_phi, B_R, B_z
real(wp_) :: R, z, phi
real(wp_) :: B(3)
real(wp_) :: c
complex(wp_) :: e_mode(2), e_ray(2)
@ -34,15 +34,11 @@ contains
! Update the inside/outside flag
iop(i) = iop(i) + 1
! Compute B in cartesian coordinates
! Compute magnetic field
R = norm2(x(1:2)) * cm
z = x(3) * cm
cosphi = x(1)/R * cm
sinphi = x(2)/R * cm
call equil%b_field(R, z, B_R, B_z, B_phi)
B(1) = B_R*cosphi - B_phi*sinphi
B(2) = B_R*sinphi + B_phi*cosphi
B(3) = B_z
phi = atan2(x(2), x(1))
call equil%b_field(R, z, phi, B=B)
! Get the polarisation vector of the given mode
call pol_limit(N, B, Bres, sox, e_mode(1), e_mode(2))
@ -72,33 +68,32 @@ contains
end subroutine plasma_in
subroutine plasma_out(i, equil, xv, anv, bres, sox, iop, ext, eyt)
subroutine plasma_out(i, equil, x, N, Bres, sox, iop, ext, eyt)
! Ray exits plasma
use const_and_precisions, only : cm
! subroutine arguments
integer, intent(in) :: i ! ray index
class(abstract_equil), intent(in) :: equil ! MHD equilibrium
real(wp_), intent(in) :: xv(3), anv(3)
real(wp_), intent(in) :: bres
real(wp_), intent(in) :: x(3), N(3)
real(wp_), intent(in) :: Bres
integer, intent(in) :: sox
integer, intent(inout) :: iop(:) ! in/out plasma flag
complex(wp_), intent(out) :: ext(:), eyt(:)
! local variables
real(wp_) :: rm, csphi, snphi, bphi, br, bz
real(wp_), dimension(3) :: bv, xmv
real(wp_) :: R, z, phi, B(3)
iop(i)=iop(i)+1 ! in->out
iop(i) = iop(i)+1 ! in->out
xmv=xv*0.01_wp_ ! convert from cm to m
rm=sqrt(xmv(1)**2+xmv(2)**2)
csphi=xmv(1)/rm
snphi=xmv(2)/rm
call equil%b_field(rm,xmv(3),br,bz,bphi)
bv(1)=br*csphi-bphi*snphi
bv(2)=br*snphi+bphi*csphi
bv(3)=bz
call pol_limit(anv,bv,bres,sox,ext(i),eyt(i)) ! polarization at plasma exit
! Compute magnetic field
R = norm2(x(1:2)) * cm
z = x(3) * cm
phi = atan2(x(2), x(1))
call equil%b_field(R, z, phi, B=B)
! Compute polarisation on the boundary
call pol_limit(N, B, Bres, sox, ext(i), eyt(i))
end subroutine plasma_out