396 lines
11 KiB
Fortran
396 lines
11 KiB
Fortran
module pec
|
|
use const_and_precisions, only : wp_,zero,one
|
|
implicit none
|
|
real(wp_), dimension(:), allocatable, save :: rhop_tab,rhot_tab
|
|
real(wp_), dimension(:), allocatable, save :: rtabpsi1
|
|
real(wp_), dimension(:), allocatable, save :: dvol,darea
|
|
real(wp_), dimension(:), allocatable, save :: ratjav,ratjbv,ratjplv
|
|
|
|
contains
|
|
|
|
subroutine pec_init(ipec,rt_in)
|
|
use equilibrium, only : frhotor,frhopol
|
|
use gray_params, only : nnd
|
|
use magsurf_data, only : fluxval
|
|
implicit none
|
|
! arguments
|
|
integer, intent(in) :: ipec
|
|
real(wp_), dimension(nnd), intent(in), optional :: rt_in
|
|
! local variables
|
|
integer :: it
|
|
real(wp_) :: drt,rt,rt1,rhop1
|
|
real(wp_) :: ratjai,ratjbi,ratjpli
|
|
real(wp_) :: voli0,voli1,areai0,areai1
|
|
|
|
! rt_in present: read input grid
|
|
! else: build equidistant grid dimension nnd
|
|
|
|
! ipec=1 rho_pol grid
|
|
! ipec=2 rho_tor grid
|
|
call dealloc_pec
|
|
allocate(rhop_tab(nnd),rhot_tab(nnd),rtabpsi1(0:nnd),dvol(nnd),darea(nnd), &
|
|
ratjav(nnd),ratjbv(nnd),ratjplv(nnd))
|
|
|
|
voli0 = zero
|
|
areai0 = zero
|
|
rtabpsi1(0) = zero
|
|
|
|
do it=1,nnd
|
|
if(present(rt_in)) then
|
|
! read radial grid from input
|
|
rt = rt_in(it)
|
|
if(it<nnd) then
|
|
drt = rt_in(it+1)-rt
|
|
end if
|
|
else
|
|
! build equidistant radial grid
|
|
drt = one/dble(nnd-1)
|
|
rt = dble(it-1)*drt
|
|
end if
|
|
! radial coordinate of i-(i+1) interval mid point
|
|
if(it < nnd) then
|
|
rt1 = rt + drt/2.0_wp_
|
|
else
|
|
rt1 = one
|
|
end if
|
|
if (ipec == 1) then
|
|
rhop_tab(it) = rt
|
|
rhot_tab(it) = frhotor(rt)
|
|
rhop1 = rt1
|
|
else
|
|
rhot_tab(it) = rt
|
|
rhop_tab(it) = frhopol(rt)
|
|
rhop1 = frhopol(rt1)
|
|
end if
|
|
! psi grid at mid points, dimension nnd+1, for use in pec_tab
|
|
rtabpsi1(it) = rhop1**2
|
|
|
|
call fluxval(rhop1,area=areai1,vol=voli1)
|
|
dvol(it) = abs(voli1 - voli0)
|
|
darea(it) = abs(areai1 - areai0)
|
|
voli0 = voli1
|
|
areai0 = areai1
|
|
|
|
call fluxval(rhop_tab(it),ratja=ratjai,ratjb=ratjbi,ratjpl=ratjpli)
|
|
ratjav(it) = ratjai
|
|
ratjbv(it) = ratjbi
|
|
ratjplv(it) = ratjpli
|
|
end do
|
|
end subroutine pec_init
|
|
|
|
|
|
subroutine spec(psjki,ppabs,ccci,iiv,pabs,currt,dpdv,ajphiv,ajcd,pins,currins)
|
|
use gray_params, only : nnd
|
|
use beamdata, only : nray,nstep
|
|
implicit none
|
|
! local constants
|
|
real(wp_), parameter :: rtbc=one
|
|
! arguments
|
|
real(wp_), dimension(nray,nstep), intent(in) :: psjki,ppabs,ccci
|
|
integer, dimension(nray), intent(in) :: iiv
|
|
real(wp_), intent(in) :: pabs,currt
|
|
real(wp_), dimension(nnd), intent(out) :: dpdv,ajphiv,ajcd,pins,currins
|
|
! local variables
|
|
integer :: i,ii,jk
|
|
real(wp_) :: spds,sccs,facpds,facjs
|
|
real(wp_), dimension(nstep):: xxi,ypt,yamp
|
|
real(wp_), dimension(nnd) :: wdpdv,wajphiv
|
|
|
|
! calculation of dP and dI over radial grid
|
|
dpdv=zero
|
|
ajphiv=zero
|
|
do jk=1,nray
|
|
ii=iiv(jk)
|
|
if (ii < nstep ) then
|
|
if(psjki(jk,ii+1) /= zero) ii=ii+1 !!! CHECK
|
|
end if
|
|
xxi=zero
|
|
ypt=zero
|
|
yamp=zero
|
|
do i=1,ii
|
|
xxi(i)=abs(psjki(jk,i))
|
|
if(xxi(i) <= one) then
|
|
ypt(i)=ppabs(jk,i)
|
|
yamp(i)=ccci(jk,i)
|
|
end if
|
|
end do
|
|
call pec_tab(xxi,ypt,yamp,ii,rtabpsi1,wdpdv,wajphiv)
|
|
dpdv = dpdv + wdpdv
|
|
ajphiv = ajphiv + wajphiv
|
|
end do
|
|
|
|
! here dpdv is still dP (not divided yet by dV)
|
|
! here ajphiv is still dI (not divided yet by dA)
|
|
spds=zero
|
|
sccs=zero
|
|
do i=1,nnd
|
|
spds=spds+dpdv(i)
|
|
sccs=sccs+ajphiv(i)
|
|
pins(i)=spds
|
|
currins(i)=sccs
|
|
end do
|
|
|
|
facpds=one
|
|
facjs=one
|
|
if(spds > zero) facpds=pabs/spds
|
|
if(sccs /= zero) facjs=currt/sccs
|
|
|
|
dpdv=facpds*(dpdv/dvol)
|
|
ajphiv=facjs*(ajphiv/darea)
|
|
ajcd=ajphiv*ratjbv
|
|
pins=facpds*pins
|
|
currins=facjs*currins
|
|
|
|
! now dpdv is dP/dV [MW/m^3]
|
|
! now ajphiv is J_phi=dI/dA [MA/m^2]
|
|
end subroutine spec
|
|
|
|
|
|
|
|
subroutine pec_tab(xxi,ypt,yamp,ii,xtab1,wdpdv,wajphiv)
|
|
! Power and current projected on psi grid - mid points
|
|
use const_and_precisions, only : wp_,one,zero
|
|
use gray_params, only : nnd
|
|
use utils, only : locatex,intlin
|
|
! arguments
|
|
integer, intent(in) :: ii
|
|
real(wp_), dimension(ii), intent(in) :: xxi,ypt,yamp
|
|
real(wp_), dimension(0:nnd), intent(in) :: xtab1
|
|
real(wp_), dimension(nnd), intent(out) :: wdpdv,wajphiv
|
|
! local variables
|
|
integer, parameter :: llmx = 21
|
|
integer, dimension(llmx) ::isev
|
|
real(wp_) :: ppa1,ppa2,cci1,cci2,dppa,didst,rt1
|
|
integer :: i,is,ise0,idecr,iise0,iise,iis,iis1
|
|
integer :: ind1,ind2,iind,ind,indi,itb1
|
|
|
|
isev = 0
|
|
ise0 = 0
|
|
idecr = -1
|
|
is = 1
|
|
wdpdv = zero
|
|
wajphiv = zero
|
|
do i=1,ii
|
|
if(ise0 == 0) then
|
|
if(xxi(i) < one) then
|
|
ise0 = i
|
|
isev(is) = max(i - 1, 1)
|
|
is = is + 1
|
|
end if
|
|
else
|
|
if (idecr == -1) then
|
|
if(xxi(i) > xxi(i-1)) then
|
|
isev(is) = i - 1
|
|
is = is + 1
|
|
idecr = 1
|
|
end if
|
|
else
|
|
if(xxi(i) > one) exit
|
|
if(xxi(i) < xxi(i-1)) then
|
|
isev(is) = i - 1
|
|
is = is + 1
|
|
idecr = -1
|
|
end if
|
|
end if
|
|
end if
|
|
end do
|
|
|
|
isev(is) = i-1
|
|
ppa1 = zero
|
|
cci1 = zero
|
|
|
|
do iis=1,is-1
|
|
iis1 = iis + 1
|
|
iise0 = isev(iis)
|
|
iise = isev(iis1)
|
|
if (mod(iis,2) /= 0) then
|
|
idecr = -1
|
|
ind1 = nnd
|
|
ind2 = 2
|
|
iind = -1
|
|
else
|
|
idecr = 1
|
|
ind1 = 1
|
|
ind2 = nnd
|
|
iind = 1
|
|
end if
|
|
do ind=ind1,ind2,iind
|
|
indi = ind
|
|
if (idecr == -1) indi = ind - 1
|
|
rt1 = xtab1(indi)
|
|
call locatex(xxi,iise,iise0,iise,rt1,itb1)
|
|
if(itb1 >= iise0 .and. itb1 <= iise) then
|
|
if(itb1 == iise) then
|
|
ppa2=ypt(itb1)
|
|
cci2=yamp(itb1)
|
|
else
|
|
call intlin(xxi(itb1), ypt(itb1),xxi(itb1+1), ypt(itb1+1),rt1,ppa2)
|
|
call intlin(xxi(itb1),yamp(itb1),xxi(itb1+1),yamp(itb1+1),rt1,cci2)
|
|
end if
|
|
dppa = ppa2 - ppa1
|
|
didst = cci2 - cci1
|
|
wdpdv(ind) = wdpdv(ind) + dppa
|
|
wajphiv(ind) = wajphiv(ind) + didst
|
|
ppa1 = ppa2
|
|
cci1 = cci2
|
|
end if
|
|
end do
|
|
end do
|
|
end subroutine pec_tab
|
|
|
|
|
|
subroutine postproc_profiles(pabs,currt,rhot_tab,dpdv,ajphiv, &
|
|
rhotpav,drhotpav,rhotjava,drhotjava,dpdvp, ajphip, &
|
|
rhotp, drhotp, rhotjfi, drhotjfi, dpdvmx,ajmxfi, ratjamx, ratjbmx)
|
|
! radial average values over power and current density profile
|
|
use const_and_precisions, only : pi, comp_eps
|
|
use gray_params, only : nnd
|
|
use equilibrium, only : frhopol
|
|
use magsurf_data, only : fluxval
|
|
implicit none
|
|
real(wp_), intent(in) :: pabs,currt
|
|
real(wp_), dimension(nnd), intent(in) :: rhot_tab
|
|
real(wp_), dimension(nnd), intent(in) :: dpdv,ajphiv
|
|
real(wp_), intent(out) :: rhotpav,drhotpav,dpdvp
|
|
real(wp_), intent(out) :: rhotjava,drhotjava,ajphip
|
|
|
|
real(wp_), intent(out) :: rhotp,drhotp,dpdvmx
|
|
real(wp_), intent(out) :: rhotjfi,drhotjfi,ajmxfi
|
|
real(wp_), intent(out) :: ratjamx,ratjbmx
|
|
|
|
real(wp_) :: sccsa,ratjplmx,rhopjava,rhoppav
|
|
real(wp_) :: rhotjav,rhot2pav,rhot2java,dvdrhotav,dadrhotava
|
|
|
|
rhotpav=zero
|
|
rhot2pav=zero
|
|
rhotjav=zero
|
|
rhotjava=zero
|
|
rhot2java=zero
|
|
|
|
if (pabs > zero) then
|
|
rhotpav = sum(rhot_tab *dpdv*dvol)/pabs
|
|
rhot2pav = sum(rhot_tab**2*dpdv*dvol)/pabs
|
|
end if
|
|
|
|
if (abs(currt) > zero) then
|
|
rhotjav = sum(rhot_tab*ajphiv*darea)/currt
|
|
end if
|
|
sccsa = sum(abs(ajphiv)*darea)
|
|
if (sccsa > zero) then
|
|
rhotjava = sum(rhot_tab *abs(ajphiv)*darea)/sccsa
|
|
rhot2java = sum(rhot_tab**2*abs(ajphiv)*darea)/sccsa
|
|
end if
|
|
|
|
! factor sqrt(8) = 2sqrt(2) to match full width of gaussian profile
|
|
drhotpav = sqrt(8._wp_*(max(rhot2pav -rhotpav**2,comp_eps)))
|
|
drhotjava = sqrt(8._wp_*(max(rhot2java-rhotjava**2,comp_eps)))
|
|
|
|
rhoppav = frhopol(rhotpav)
|
|
rhopjava = frhopol(rhotjava)
|
|
|
|
if (pabs > zero) then
|
|
call fluxval(rhoppav,dvdrhot=dvdrhotav)
|
|
dpdvp = pabs*2.0_wp_/(sqrt(pi)*drhotpav*dvdrhotav)
|
|
call profwidth(nnd,rhot_tab,dpdv,rhotp,dpdvmx,drhotp)
|
|
else
|
|
dpdvp = zero
|
|
rhotp = zero
|
|
dpdvmx = zero
|
|
drhotp = zero
|
|
end if
|
|
|
|
if (sccsa > zero) then
|
|
call fluxval(rhopjava,dadrhot=dadrhotava,ratja=ratjamx,ratjb=ratjbmx, &
|
|
ratjpl=ratjplmx)
|
|
ajphip = currt*2.0_wp_/(sqrt(pi)*drhotjava*dadrhotava)
|
|
call profwidth(nnd,rhot_tab,ajphiv,rhotjfi,ajmxfi,drhotjfi)
|
|
else
|
|
ajphip = zero
|
|
rhotjfi = zero
|
|
ajmxfi = zero
|
|
drhotjfi = zero
|
|
end if
|
|
end subroutine postproc_profiles
|
|
|
|
|
|
|
|
subroutine profwidth(nd,xx,yy,xpk,ypk,dxxe)
|
|
use const_and_precisions, only : wp_,emn1
|
|
use utils, only : locatex, locate, intlin, vmaxmini
|
|
implicit none
|
|
! arguments
|
|
integer :: nd
|
|
real(wp_), dimension(nd) :: xx,yy
|
|
real(wp_), intent(out) :: xpk,ypk,dxxe
|
|
! local variables
|
|
integer :: imn,imx,ipk,ie
|
|
real(wp_) :: xmn,xmx,ymn,ymx,xpkp,xpkm,yye,rte1,rte2
|
|
real(wp_) :: ypkp,ypkm
|
|
|
|
call vmaxmini(yy,nd,ymn,ymx,imn,imx)
|
|
ypk = zero
|
|
xmx = xx(imx)
|
|
xmn = xx(imn)
|
|
if (abs(ymx) > abs(ymn)) then
|
|
ipk = imx
|
|
ypkp = ymx
|
|
xpkp = xmx
|
|
if(abs(ymn/ymx) < 1.0e-2_wp_) ymn = 0.0_wp_
|
|
ypkm = ymn
|
|
xpkm = xmn
|
|
else
|
|
ipk = imn
|
|
ypkp = ymn
|
|
xpkp = xmn
|
|
if(abs(ymx/ymn) < 1.0e-2_wp_) ymx = 0.0_wp_
|
|
ypkm = ymx
|
|
xpkm = xmx
|
|
end if
|
|
if(xpkp > zero) then
|
|
xpk = xpkp
|
|
ypk = ypkp
|
|
yye = ypk*emn1
|
|
call locatex(yy,nd,1,ipk,yye,ie)
|
|
if(ie > 0 .and. ie < nd) then
|
|
call intlin(yy(ie),xx(ie),yy(ie+1),xx(ie+1),yye,rte1)
|
|
else
|
|
rte1 = zero
|
|
end if
|
|
call locatex(yy,nd,ipk,nd,yye,ie)
|
|
if(ie > 0 .and. ie < nd) then
|
|
call intlin(yy(ie),xx(ie),yy(ie+1),xx(ie+1),yye,rte2)
|
|
else
|
|
rte2 = zero
|
|
end if
|
|
else
|
|
ipk=2
|
|
xpk=xx(2)
|
|
ypk=yy(2)
|
|
rte1=0.0_wp_
|
|
yye=ypk*emn1
|
|
call locate(yy,nd,yye,ie)
|
|
if(ie > 0 .and. ie < nd) then
|
|
call intlin(yy(ie),xx(ie),yy(ie+1),xx(ie+1),yye,rte2)
|
|
else
|
|
rte2 = zero
|
|
end if
|
|
end if
|
|
dxxe = rte2 - rte1
|
|
if(ymx /= zero .and. ymn /= zero) dxxe = -dxxe
|
|
end subroutine profwidth
|
|
|
|
subroutine dealloc_pec
|
|
implicit none
|
|
|
|
if (allocated(rhop_tab)) deallocate(rhop_tab)
|
|
if (allocated(rhot_tab)) deallocate(rhot_tab)
|
|
if (allocated(rtabpsi1)) deallocate(rtabpsi1)
|
|
if (allocated(dvol)) deallocate(dvol)
|
|
if (allocated(darea)) deallocate(darea)
|
|
if (allocated(ratjav)) deallocate(ratjav)
|
|
if (allocated(ratjbv)) deallocate(ratjbv)
|
|
if (allocated(ratjplv)) deallocate(ratjplv)
|
|
end subroutine dealloc_pec
|
|
|
|
end module pec
|