Calculates density and temperature profile quantities author: P J Knight, CCFE, Culham Science Centre None This subroutine initialises the density and temperature profile averages and peak values, given the main parameters describing these profiles. T&M/PKNIGHT/LOGBOOK24, pp.4-7
subroutine plasma_profiles
!! Calculates density and temperature profile quantities
!! author: P J Knight, CCFE, Culham Science Centre
!! None
!! This subroutine initialises the density and temperature
!! profile averages and peak values, given the main
!! parameters describing these profiles.
!! T&M/PKNIGHT/LOGBOOK24, pp.4-7
!
! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
use constants, only: echarge
use divertor_variables, only: prn1
use maths_library, only: gamfun, sumup3
use physics_variables, only: rhopedt, ten, tin, alphap, tbeta, te0, p0, &
nesep, tesep, pcoef, ipedestal, ni0, ne0, ti0, tratio, dnla, alphat, &
dnitot, neped, ti, rhopedn, dene, teped, alphan, te, rho_ne_max, &
rho_te_max, gradient_length_ne, gradient_length_te, rminor
implicit none
! Arguments
! Local variables
integer, parameter :: nrho = 501
integer :: irho
real(dp) :: drho, rho, integ1, integ2, dens, temp, te_max, ne_max, dndrho_max, dtdrho_max
real(dp), dimension(nrho) :: arg1, arg2, arg3
! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Volume-averaged ion temperature
! (input value used directly if tratio=0.0)
if (tratio > 0.0D0) ti = tratio * te
if (ipedestal == 0) then
! Reset pedestal values to agree with original parabolic profiles
rhopedt = 1.0D0 ; rhopedn = 1.0D0
teped = 0.0D0 ; tesep = 0.0D0
neped = 0.0D0 ; nesep = 0.0D0
tbeta = 2.0D0
! Profile factor; ratio of density-weighted to volume-averaged
! temperature
pcoef = (1.0D0 + alphan)*(1.0D0 + alphat)/(1.0D0+alphan+alphat)
! Line averaged electron density (IPDG89)
! 0.5*gamfun(0.5) = 0.5*sqrt(pi) = 0.886227
dnla = dene*(1.0D0+alphan) * 0.886227D0 * gamfun(alphan+1.0D0) / &
gamfun(alphan+1.5D0)
! Density-weighted temperatures
ten = te * pcoef
tin = ti * pcoef
! Central values for temperature (keV) and density (m**-3)
te0 = te * (1.0D0+alphat)
ti0 = ti * (1.0D0+alphat)
ne0 = dene * (1.0D0+alphan)
ni0 = dnitot * (1.0D0+alphan)
else
! The following reproduces the above results within sensible
! tolerances if rhopedt = rhopedn = 1.0, teped = tesep = neped
! = nesep = 0.0, and tbeta = 2.0
! Central values for temperature (keV) and density (m**-3)
te0 = tcore(rhopedt,teped,tesep,te,alphat,tbeta)
ti0 = ti/te * te0
ne0 = ncore(rhopedn,neped,nesep,dene,alphan)
ni0 = dnitot/dene * ne0
! Perform integrations to calculate ratio of density-weighted
! to volume-averaged temperature, etc.
! Density-weighted temperature = integral(n.T dV) / integral(n dV)
! which is approximately equal to the ratio
! integral(rho.n(rho).T(rho) drho) / integral(rho.n(rho) drho)
drho = 1.0D0/(nrho-1)
do irho = 1,nrho
rho = (irho-1.0D0)/(nrho-1)
dens = nprofile(rho,rhopedn,ne0,neped,nesep,alphan)
temp = tprofile(rho,rhopedt,te0,teped,tesep,alphat,tbeta)
arg1(irho) = rho*dens*temp
arg2(irho) = rho*dens
arg3(irho) = dens
end do
call sumup3(drho,arg1,integ1,nrho)
call sumup3(drho,arg2,integ2,nrho)
! Density-weighted temperatures
ten = integ1/integ2
tin = ti/te * ten
! Profile factor; ratio of density-weighted to volume-averaged
! temperature
pcoef = ten / te
! Line-averaged electron density
! = integral(n(rho).drho)
call sumup3(drho,arg3,dnla,nrho)
! Scrape-off density / volume averaged density
! (Input value is used if ipedestal = 0)
prn1 = max(0.01D0, nesep/dene) ! preventing division by zero later
end if
! Central pressure (Pa), from ideal gas law : p = nkT
p0 = (ne0*te0 + ni0*ti0) * 1.0D3 * echarge
! Pressure profile index (N.B. no pedestal effects included here)
! N.B. p0 is NOT equal to <p> * (1 + alphap), but p(rho) = n(rho)*T(rho)
! and <p> = <n>.T_n where <...> denotes volume-averages and T_n is the
! density-weighted temperature
alphap = alphan + alphat
! The gradient information for ipedestal = 0:
! All formulas can be obtained from the analytical parametric form of the ipedestal profiles
! rho_max is obtained by equalling the second derivative to zero e.g.
if (ipedestal == 0) then
if(alphat > 1.0d0) then
! Rho (normalized radius), where temperature derivative is largest
rho_te_max = 1.0d0/sqrt(-1.0d0 +2.0d0 * alphat)
dtdrho_max = -2.0d0**alphat*(-1.0d0 + alphat)**(-1.0d0 + alphat)*alphat*(-1.0d0 + &
2.0d0 * alphat)**(0.5d0 - alphat)*te0
te_max = te0*(1d0 - rho_te_max**2)**alphat
elseif (alphat .le. 1.0d0 .and. alphaT > 0.0d0) then
! This makes the profiles very 'boxy'
! The gradient diverges here at the edge so define some 'wrong' value of 0.9
! to approximate the gradient
rho_te_max = 0.9d0
dtdrho_max = -2.0d0 * alphat * rho_te_max*(1.d0-rho_te_max**2)**(-1.0d0+alphat)*te0
te_max = te0*(1d0 - rho_te_max**2)**alphat
else
print *, "ERROR: alphat is negative!"
call exit(1)
end if
! Same for density
if(alphan > 1.0d0) then
rho_ne_max = 1.0d0/sqrt(-1.0d0 +2.0d0 * alphan)
dndrho_max = -2.0d0**alphan*(-1.0d0 + alphan)**(-1.0d0 + alphan)*alphan*(-1.0d0 + &
2.0d0 * alphan)**(0.5d0 - alphan)*ne0
ne_max = ne0*(1d0 - rho_ne_max**2)**alphan
elseif (alphan .le. 1.0d0 .and. alphan > 0.0d0) then
! This makes the profiles very 'boxy'
! The gradient diverges here at the edge so define some 'wrong' value of 0.9
! to approximate the gradient
rho_ne_max = 0.9d0
dndrho_max = -2.0d0 * alphan * rho_ne_max*(1.d0-rho_ne_max**2)**(-1.0d0+alphan)*ne0
ne_max = ne0*(1d0 - rho_ne_max**2)**alphan
else
print *, "ERROR: alphan is negative!"
call exit(1)
end if
! set normalized gradient length
! te at rho_te_max
gradient_length_te = -dtdrho_max * rminor*rho_te_max/te_max
! same for density:
gradient_length_ne = -dndrho_max * rminor*rho_ne_max/ne_max
end if
end subroutine plasma_profiles