pec subroutine(s) updated
This commit is contained in:
Daniela Farina
parent
771cdb3822
commit
e92ff7cee1
4
Makefile
4
Makefile
@ -6,7 +6,7 @@ MAINOBJ=gray.o
|
|||||||
OTHOBJ=conical.o const_and_precisions.o dierckx.o dispersion.o eierf.o \
|
OTHOBJ=conical.o const_and_precisions.o dierckx.o dispersion.o eierf.o \
|
||||||
graydata_anequil.o graydata_flags.o graydata_par.o green_func_p.o \
|
graydata_anequil.o graydata_flags.o graydata_par.o green_func_p.o \
|
||||||
interp_eqprof.o magsurf_data.o math.o minpack.o numint.o quadpack.o \
|
interp_eqprof.o magsurf_data.o math.o minpack.o numint.o quadpack.o \
|
||||||
reflections.o simplespline.o utils.o beamdata.o
|
reflections.o simplespline.o utils.o beamdata.o
|
||||||
|
|
||||||
# Alternative search paths
|
# Alternative search paths
|
||||||
vpath %.f90 src
|
vpath %.f90 src
|
||||||
@ -28,7 +28,7 @@ $(EXE): $(MAINOBJ) $(OTHOBJ)
|
|||||||
gray.o: const_and_precisions.o conical.o dierckx.o dispersion.o \
|
gray.o: const_and_precisions.o conical.o dierckx.o dispersion.o \
|
||||||
graydata_anequil.o graydata_flags.o graydata_par.o green_func_p.o \
|
graydata_anequil.o graydata_flags.o graydata_par.o green_func_p.o \
|
||||||
interp_eqprof.o magsurf_data.o math.o minpack.o numint.o quadpack.o \
|
interp_eqprof.o magsurf_data.o math.o minpack.o numint.o quadpack.o \
|
||||||
reflections.o simplespline.o utils.o beamdata.o
|
reflections.o simplespline.o utils.o beamdata.o
|
||||||
conical.o: const_and_precisions.o
|
conical.o: const_and_precisions.o
|
||||||
dierckx.o: const_and_precisions.o
|
dierckx.o: const_and_precisions.o
|
||||||
dispersion.o: const_and_precisions.o eierf.o math.o quadpack.o
|
dispersion.o: const_and_precisions.o eierf.o math.o quadpack.o
|
||||||
|
|||||||
702
src/gray.f
702
src/gray.f
@ -118,16 +118,26 @@ c
|
|||||||
subroutine after_gray_integration
|
subroutine after_gray_integration
|
||||||
use const_and_precisions, only : wp_,zero
|
use const_and_precisions, only : wp_,zero
|
||||||
use graydata_flags, only : ibeam,iwarm,iequil,iprof,
|
use graydata_flags, only : ibeam,iwarm,iequil,iprof,
|
||||||
. filenmeqq,filenmprf,filenmbm
|
. nnd,ipec,ieccd,filenmeqq,filenmprf,filenmbm
|
||||||
use graydata_anequil, only : dens0,te0
|
use graydata_anequil, only : dens0,te0
|
||||||
use beamdata, only : nrayr
|
use beamdata, only : nrayr
|
||||||
implicit none
|
implicit none
|
||||||
c local variables
|
c local variables
|
||||||
integer :: iproj,nfilp
|
integer :: iproj,nfilp,i
|
||||||
real(wp_) :: currtka,pabs,currt
|
real(wp_) :: currtka,pabs,currt
|
||||||
|
|
||||||
|
real(wp_), dimension(nnd) :: dpdv,ajphiv,pins,currins
|
||||||
|
real(wp_), dimension(nnd) :: dvol,darea,ratjav,ratjbv,ratjplv
|
||||||
|
real(wp_), dimension(nnd) :: ajcdav,ajcdbv,ajplv
|
||||||
|
real(wp_), dimension(nnd) :: rt_in,rhop_tab,rhot_tab
|
||||||
|
real(wp_), dimension(0:nnd) :: rtabpsi1
|
||||||
|
real(wp_) :: rhotpav,drhotpav
|
||||||
|
real(wp_) :: rhotjav,rhotjava,drhotjav,drhotjava
|
||||||
c common/external functions/variables
|
c common/external functions/variables
|
||||||
integer :: index_rt
|
integer :: index_rt
|
||||||
real(wp_) :: st,taumn,taumx,pabstot,currtot
|
real(wp_) :: st,taumn,taumx,pabstot,currtot
|
||||||
|
|
||||||
|
! arguments
|
||||||
c
|
c
|
||||||
common/ss/st
|
common/ss/st
|
||||||
common/index_rt/index_rt
|
common/index_rt/index_rt
|
||||||
@ -164,13 +174,30 @@ c
|
|||||||
|
|
||||||
end if
|
end if
|
||||||
c
|
c
|
||||||
c compute power and current density profiles for all rays
|
c compute power and current density profiles for all rays
|
||||||
c
|
c
|
||||||
pabs=pabstot
|
pabs=pabstot
|
||||||
currt=currtot
|
currt=currtot
|
||||||
call pec(pabs,currt)
|
|
||||||
c
|
call pec_init(ipec,rt_in,rhop_tab,rhot_tab,rtabpsi1,
|
||||||
|
. dvol,darea,ratjav,ratjbv,ratjplv)
|
||||||
|
|
||||||
|
call spec(pabs,currt,rtabpsi1,dvol,darea,dpdv,ajphiv,pins,currins)
|
||||||
|
|
||||||
|
call postproc_profiles(pabs,currt,rhot_tab,dvol,darea,dpdv,ajphiv,
|
||||||
|
. rhotpav,drhotpav,rhotjava,drhotjava)
|
||||||
|
|
||||||
|
write(*,*) 'rhotpav,drhotpav,rhotjava,drhotjava'
|
||||||
|
write(*,99) rhotpav,drhotpav,rhotjava,drhotjava
|
||||||
|
|
||||||
|
|
||||||
|
do i=1,nnd
|
||||||
|
write(48,99) rhop_tab(i),rhot_tab(i),ajphiv(i),
|
||||||
|
. ajphiv(i)*ratjbv(i),dpdv(i),currins(i),pins(i)
|
||||||
|
end do
|
||||||
|
|
||||||
return
|
return
|
||||||
|
99 format(16(1x,e16.8e3))
|
||||||
end
|
end
|
||||||
c
|
c
|
||||||
c
|
c
|
||||||
@ -605,7 +632,7 @@ c
|
|||||||
. 'rhotj rhotjava rhotp rhotpav drhotjava drhotpav ratjamx '//
|
. 'rhotj rhotjava rhotp rhotpav drhotjava drhotpav ratjamx '//
|
||||||
. 'ratjbmx stmx psipol chipol index_rt Jphimx dPdVmx drhotj '//
|
. 'ratjbmx stmx psipol chipol index_rt Jphimx dPdVmx drhotj '//
|
||||||
. 'drhotp'
|
. 'drhotp'
|
||||||
write(48,*) '#psi rhot Jphi Jcdb dPdV Icdins Pins P% index_rt'
|
write(48,*) '#rhop rhot Jphi Jcdb dPdV Icdins Pins'
|
||||||
write(66,*) "# psipol0 chipol0 powrfl"
|
write(66,*) "# psipol0 chipol0 powrfl"
|
||||||
c
|
c
|
||||||
else
|
else
|
||||||
@ -6374,376 +6401,335 @@ c
|
|||||||
99 format(i5,12(1x,e16.8e3))
|
99 format(i5,12(1x,e16.8e3))
|
||||||
111 format(3i5,12(1x,e16.8e3))
|
111 format(3i5,12(1x,e16.8e3))
|
||||||
end
|
end
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
subroutine pec_init(ipec,rt_in,rhop_tab,rhot_tab,rtabpsi1,
|
||||||
|
. dvol,darea,ratjav,ratjbv,ratjplv)
|
||||||
|
use const_and_precisions, only : wp_,zero,one,izero
|
||||||
|
use graydata_flags, only : nnd
|
||||||
|
implicit none
|
||||||
|
integer :: it,ipec
|
||||||
|
real(wp_) :: drt,rt,rt1,rhop1
|
||||||
|
real(wp_) :: ratjai,ratjbi,ratjpli
|
||||||
|
real(wp_) :: voli0,voli1,dervoli,areai0,areai1
|
||||||
|
real(wp_) :: rrii,rbavi,bmxi,bmni,fci
|
||||||
|
real(wp_), dimension(nnd), intent(in) :: rt_in
|
||||||
|
real(wp_), dimension(nnd), intent(out) :: rhop_tab,rhot_tab
|
||||||
|
real(wp_), dimension(0:nnd), intent(out) :: rtabpsi1
|
||||||
|
real(wp_), dimension(nnd), intent(out) :: dvol,darea
|
||||||
|
real(wp_), dimension(nnd), intent(out) :: ratjav,ratjbv,ratjplv
|
||||||
|
real(wp_), external :: frhotor,frhopol
|
||||||
|
|
||||||
|
c ipec positive build equidistant grid dimension nnd
|
||||||
|
c ipec negative read input grid
|
||||||
c
|
c
|
||||||
c
|
c ipec=+/-1 rho_pol grid
|
||||||
c
|
c ipec=+/-2 rho_tor grid
|
||||||
subroutine pec(pabs,currt)
|
|
||||||
use const_and_precisions, only : wp_,pi,one
|
voli0=zero
|
||||||
use numint, only : simpson
|
areai0=zero
|
||||||
use graydata_flags, only : ipec,ieccd,nnd,dst
|
rtabpsi1(0)=zero
|
||||||
use utils, only : locatex, locate, intlin
|
|
||||||
use beamdata, only : ppabs,ccci,psjki,nrayr,nrayth,nstep,iiv,pdjki
|
do it=1,nnd
|
||||||
|
if(ipec > 0) then
|
||||||
|
c build equidistant radial grid
|
||||||
|
drt=one/dble(nnd-1)
|
||||||
|
rt=dble(it-1)*drt
|
||||||
|
else
|
||||||
|
c read radial grid from input
|
||||||
|
rt=rt_in(it)
|
||||||
|
drt=(rt_in(it+1)-rt)/2.0_wp_
|
||||||
|
end if
|
||||||
|
c 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 (abs(ipec) == 1) then
|
||||||
|
rhop_tab(it)=rt
|
||||||
|
rhot_tab(it)=frhotor(rt*rt)
|
||||||
|
rhop1=rt1
|
||||||
|
else
|
||||||
|
rhot_tab(it)=rt
|
||||||
|
rhop_tab(it)=frhopol(rt)
|
||||||
|
rhop1=frhopol(rt1)
|
||||||
|
end if
|
||||||
|
c psi grid at mid points, dimension nnd+1, for use in pec_tab
|
||||||
|
rtabpsi1(it)=rhop1**2
|
||||||
|
|
||||||
|
call valpsispl(rhop1,voli1,dervoli,areai1,rrii,
|
||||||
|
. rbavi,bmxi,bmni,fci,izero)
|
||||||
|
dvol(it)=abs(voli1-voli0)
|
||||||
|
darea(it)=abs(areai1-areai0)
|
||||||
|
voli0=voli1
|
||||||
|
areai0=areai1
|
||||||
|
|
||||||
|
call ratioj(rhop_tab(it),ratjai,ratjbi,ratjpli)
|
||||||
|
ratjav(it)=ratjai
|
||||||
|
ratjbv(it)=ratjbi
|
||||||
|
ratjplv(it)=ratjpli
|
||||||
|
end do
|
||||||
|
end subroutine pec_init
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
subroutine spec(pabs,currt,rtabpsi1,dvol,darea,dpdv,ajphiv,
|
||||||
|
. pins,currins)
|
||||||
|
use const_and_precisions, only : wp_,one,zero,izero
|
||||||
|
use graydata_flags, only : nnd,ipec,ieccd
|
||||||
|
use beamdata, only : ppabs,ccci,psjki,nrayr,nrayth,nstep,iiv
|
||||||
implicit none
|
implicit none
|
||||||
c local constants
|
c local constants
|
||||||
integer, parameter :: nndmx=5001,llmx=21
|
|
||||||
real(wp_), parameter :: rtbc=one
|
real(wp_), parameter :: rtbc=one
|
||||||
c arguments
|
c arguments
|
||||||
real(wp_) :: pabs,currt
|
real(wp_), intent(in) :: pabs,currt
|
||||||
|
real(wp_), dimension(nstep):: xxi,ypt,yamp
|
||||||
|
real(wp_), dimension(nnd), intent(out) :: dpdv,ajphiv
|
||||||
|
real(wp_), dimension(nnd), intent(out) :: pins,currins
|
||||||
|
real(wp_), dimension(nnd), intent(in) :: dvol,darea
|
||||||
|
real(wp_), dimension(0:nnd), intent(in) :: rtabpsi1
|
||||||
c local variables
|
c local variables
|
||||||
integer :: ll,it,nd,intp,kkk,i,j,k,ii,is,idecr,ise0,iis,iis1,
|
integer :: i,ii,j,k,kkk
|
||||||
. iise0,iise,ind1,ind2,iind,indi,iif1,istmx,i1,ind,itb1
|
real(wp_) :: spds,sccs,facpds,facjs
|
||||||
integer, dimension(llmx) :: isev
|
real(wp_), dimension(nnd) :: wdpdv,wajphiv
|
||||||
real(wp_) :: stf,voli0,areai0,drt,rt,rt1,psit,psit1,voli1,
|
|
||||||
. dervoli,areai1,rrii,rbavi,bmxi,bmni,fci,ppa1,cci1,ppa2,cci2,
|
c calculation of dP and dI over radial grid
|
||||||
. dppa,didst,rhotpav,rhot2pav,rhotjav,rhotjava,rhot2java,spds,
|
dpdv=zero
|
||||||
. sccs,sccsa,drhot2pav,drhotpav,drhot2java,drhotjava,facpds,
|
ajphiv=zero
|
||||||
. facjs,rhpp,rhpj,dpdvp,ajphip,ratjamx,ratjbmx,ratjplmx,rhopp,
|
|
||||||
. drhopp,rhotjfi,rhopfi,ajmxfi,drhotjfi,drhopfi,xps,dpdvmx,
|
|
||||||
. rhotp,drhotp,stmx,pins_02,pins_05,pins_085,xrhot,currtka,rhop,
|
|
||||||
. pinsr,ratjpli,ratjai,ratjbi,frhotor,frhopol,fdadrhot,
|
|
||||||
. fdvdrhot,h,psin
|
|
||||||
real(wp_), dimension(nstep) :: xxi,ypt,yamp
|
|
||||||
real(wp_), dimension(nndmx) :: ajphiv,dpdv,dvolt,darea,ratjav,
|
|
||||||
. ratjbv,ratjplv,ajplv,ajcdav,ajcdbv,pins,currins,fi,rtab,rhotv
|
|
||||||
real(wp_), dimension(0:nndmx) :: rtab1
|
|
||||||
c common/external functions/variables
|
|
||||||
integer :: istep,index_rt
|
|
||||||
real(wp_) :: alpha0,beta0,taumn,taumx,pabstot,currtot,psipol,
|
|
||||||
. chipol
|
|
||||||
c
|
|
||||||
common/istep/istep
|
|
||||||
common/index_rt/index_rt
|
|
||||||
common/angles/alpha0,beta0
|
|
||||||
common/taumnx/taumn,taumx,pabstot,currtot
|
|
||||||
common/polcof/psipol,chipol
|
|
||||||
c
|
|
||||||
stf=istep*dst
|
|
||||||
nd=nnd
|
|
||||||
|
|
||||||
if(pabstot.gt.0.0_wp_) then
|
|
||||||
do ll=1,llmx
|
|
||||||
isev(ll)=0
|
|
||||||
end do
|
|
||||||
|
|
||||||
intp=0
|
|
||||||
voli0=0.0_wp_
|
|
||||||
areai0=0.0_wp_
|
|
||||||
rtab1(0)=0.0_wp_
|
|
||||||
do it=1,nd
|
|
||||||
drt=1.0_wp_/dble(nd-1)
|
|
||||||
rt=dble(it-1)*drt
|
|
||||||
if(it.lt.nd) then
|
|
||||||
rt1=rt+drt/2.0_wp_
|
|
||||||
else
|
|
||||||
rt1=rt
|
|
||||||
end if
|
|
||||||
rtab(it)=rt
|
|
||||||
rtab1(it)=rt1
|
|
||||||
dpdv(it)=0.0_wp_
|
|
||||||
ajphiv(it)=0.0_wp_
|
|
||||||
if (ipec.eq.0) then
|
|
||||||
psit=rt
|
|
||||||
psit1=rt1
|
|
||||||
else
|
|
||||||
psit=rt**2
|
|
||||||
psit1=rt1**2
|
|
||||||
end if
|
|
||||||
rhotv(it)=frhotor(psit)
|
|
||||||
call valpsispl(sqrt(psit1),voli1,dervoli,areai1,rrii,
|
|
||||||
. rbavi,bmxi,bmni,fci,intp)
|
|
||||||
dvolt(it)=abs(voli1-voli0)
|
|
||||||
darea(it)=abs(areai1-areai0)
|
|
||||||
voli0=voli1
|
|
||||||
areai0=areai1
|
|
||||||
call ratioj(sqrt(psit),ratjai,ratjbi,ratjpli)
|
|
||||||
ratjav(it)=ratjai
|
|
||||||
ratjbv(it)=ratjbi
|
|
||||||
ratjplv(it)=ratjpli
|
|
||||||
end do
|
|
||||||
|
|
||||||
kkk=1
|
kkk=1
|
||||||
do j=1,nrayr
|
do j=1,nrayr
|
||||||
if(j.gt.1) kkk=nrayth
|
if(j > 1) kkk=nrayth
|
||||||
do k=1,kkk
|
do k=1,kkk
|
||||||
ise0=0
|
|
||||||
ii=iiv(j,k)
|
ii=iiv(j,k)
|
||||||
if (ii.lt.nstep) then
|
if (ii < nstep .and. psjki(j,k,ii+1) /= zero) ii=ii+1
|
||||||
if(psjki(j,k,ii+1).ne.0.0_wp_) ii=ii+1
|
xxi=zero
|
||||||
end if
|
ypt=zero
|
||||||
idecr=-1
|
yamp=zero
|
||||||
is=1
|
|
||||||
do i=1,ii
|
do i=1,ii
|
||||||
if(ipec.eq.0) xxi(i)=abs(psjki(j,k,i))
|
xxi(i)=abs(psjki(j,k,i))
|
||||||
if(ipec.eq.1) xxi(i)=sqrt(abs(psjki(j,k,i)))
|
if(xxi(i) <= one) then
|
||||||
if(psjki(j,k,i).ge.0.and.psjki(j,k,i).le.rtbc) then
|
|
||||||
ypt(i)=ppabs(j,k,i)
|
ypt(i)=ppabs(j,k,i)
|
||||||
yamp(i)=ccci(j,k,i)
|
yamp(i)=ccci(j,k,i)
|
||||||
else
|
|
||||||
ypt(i)=0.0_wp_
|
|
||||||
yamp(i)=0.0_wp_
|
|
||||||
end if
|
|
||||||
if(ise0.eq.0) then
|
|
||||||
if(xxi(i).lt.rtbc) then
|
|
||||||
ise0=i
|
|
||||||
isev(is)=i-1
|
|
||||||
is=is+1
|
|
||||||
end if
|
|
||||||
else
|
|
||||||
if (idecr.eq.-1) then
|
|
||||||
if(xxi(i).gt.xxi(i-1)) then
|
|
||||||
isev(is)=i-1
|
|
||||||
is=is+1
|
|
||||||
idecr=1
|
|
||||||
end if
|
|
||||||
else
|
|
||||||
if(xxi(i).gt.rtbc) exit
|
|
||||||
if(xxi(i).lt.xxi(i-1)) then
|
|
||||||
c isev(is)=i !!!!!!!!!! it should be isev(is)=i-1
|
|
||||||
isev(is)=i-1
|
|
||||||
is=is+1
|
|
||||||
idecr=-1
|
|
||||||
end if
|
|
||||||
end if
|
|
||||||
end if
|
end if
|
||||||
end do
|
end do
|
||||||
c
|
call pec_tab(xxi,ypt,yamp,ii,rtabpsi1,
|
||||||
isev(is)=i-1
|
. wdpdv,wajphiv)
|
||||||
ppa1=0.0_wp_
|
dpdv=dpdv+wdpdv
|
||||||
cci1=0.0_wp_
|
ajphiv=ajphiv+wajphiv
|
||||||
do iis=1,is-1
|
|
||||||
iis1=iis+1
|
|
||||||
iise0=isev(iis)
|
|
||||||
iise=isev(iis1)
|
|
||||||
if (mod(iis,2).ne.0) then
|
|
||||||
idecr=-1
|
|
||||||
ind1=nd
|
|
||||||
ind2=2
|
|
||||||
iind=-1
|
|
||||||
else
|
|
||||||
idecr=1
|
|
||||||
ind1=1
|
|
||||||
ind2=nd
|
|
||||||
iind=1
|
|
||||||
end if
|
|
||||||
do ind=ind1,ind2,iind !!!!!!!!!! is it safe?
|
|
||||||
c iind=ind !!!!!!!!!! iind reused in the loop!
|
|
||||||
indi=ind !!!!!!!!!! iind reused in the loop!
|
|
||||||
if (idecr.eq.-1) indi=ind-1
|
|
||||||
rt1=rtab1(indi)
|
|
||||||
call locatex(xxi,iise,iise0,iise,rt1,itb1)
|
|
||||||
if(itb1.ge.iise0.and.itb1.lt.iise) then
|
|
||||||
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)
|
|
||||||
dppa=ppa2-ppa1
|
|
||||||
dpdv(ind)=dpdv(ind)+dppa
|
|
||||||
didst=(cci2-cci1)
|
|
||||||
ajphiv(ind)=ajphiv(ind)+didst
|
|
||||||
ppa1=ppa2
|
|
||||||
cci1=cci2
|
|
||||||
end if
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
end do
|
end do
|
||||||
end do
|
end do
|
||||||
|
|
||||||
h=1.0_wp_/dble(nd-1)
|
c here dpdv is still dP (not divided yet by dV)
|
||||||
rhotpav=0.0_wp_
|
c here ajphiv is still dI (not divided yet by dA)
|
||||||
drhotpav=0.0_wp_
|
|
||||||
rhotjav=0.0_wp_
|
spds=zero
|
||||||
rhotjava=0.0_wp_
|
sccs=zero
|
||||||
rhot2java=0.0_wp_
|
do i=1,nnd
|
||||||
|
|
||||||
fi=dpdv/h
|
|
||||||
call simpson (nd,h,fi,spds)
|
|
||||||
fi=rhotv*dpdv/h
|
|
||||||
call simpson (nd,h,fi,rhotpav)
|
|
||||||
fi=rhotv*rhotv*dpdv/h
|
|
||||||
call simpson (nd,h,fi,rhot2pav)
|
|
||||||
rhotpav=rhotpav/spds
|
|
||||||
rhot2pav=rhot2pav/spds
|
|
||||||
|
|
||||||
if (ieccd.ne.0) then
|
|
||||||
fi=ajphiv/h
|
|
||||||
call simpson (nd,h,fi,sccs)
|
|
||||||
fi=rhotv*ajphiv/h
|
|
||||||
call simpson (nd,h,fi,rhotjav)
|
|
||||||
fi=abs(ajphiv)/h
|
|
||||||
call simpson (nd,h,fi,sccsa)
|
|
||||||
fi=rhotv*abs(ajphiv)/h
|
|
||||||
call simpson (nd,h,fi,rhotjava)
|
|
||||||
fi=rhotv*rhotv*abs(ajphiv)/h
|
|
||||||
call simpson (nd,h,fi,rhot2java)
|
|
||||||
rhotjav=rhotjav/sccs
|
|
||||||
rhotjava=rhotjava/sccsa
|
|
||||||
rhot2java=rhot2java/sccsa
|
|
||||||
end if
|
|
||||||
|
|
||||||
c factor sqrt(8)=2 sqrt(2) to match with full width
|
|
||||||
c of gaussian profile
|
|
||||||
drhot2pav=rhot2pav-rhotpav**2
|
|
||||||
drhotpav=sqrt(8.0_wp_*drhot2pav)
|
|
||||||
drhot2java=rhot2java-rhotjava**2
|
|
||||||
drhotjava=sqrt(8.0_wp_*drhot2java)
|
|
||||||
|
|
||||||
spds=0.0_wp_
|
|
||||||
sccs=0.0_wp_
|
|
||||||
do i=1,nd
|
|
||||||
spds=spds+dpdv(i)
|
spds=spds+dpdv(i)
|
||||||
sccs=sccs+ajphiv(i)
|
sccs=sccs+ajphiv(i)
|
||||||
pins(i)=spds
|
pins(i)=spds
|
||||||
currins(i)=sccs
|
currins(i)=sccs
|
||||||
dpdv(i)=dpdv(i)/dvolt(i)
|
dpdv(i)=dpdv(i)/dvol(i)
|
||||||
ajphiv(i)=ajphiv(i)/darea(i)
|
ajphiv(i)=ajphiv(i)/darea(i)
|
||||||
end do
|
end do
|
||||||
|
|
||||||
facpds=1.0_wp_
|
facpds=one
|
||||||
facjs=1.0_wp_
|
facjs=one
|
||||||
if(spds.gt.0.0_wp_) facpds=pabs/spds
|
if(spds > zero) facpds=pabs/spds
|
||||||
if(sccs.ne.0.0_wp_) facjs=currt/sccs
|
if(sccs /= zero) facjs=currt/sccs
|
||||||
|
|
||||||
do i=1,nd
|
do i=1,nnd
|
||||||
dpdv(i)=facpds*dpdv(i)
|
dpdv(i)=facpds*dpdv(i)
|
||||||
ajphiv(i)=facjs*ajphiv(i)
|
ajphiv(i)=facjs*ajphiv(i)
|
||||||
ajcdav(i)=ajphiv(i)*ratjav(i)
|
end do
|
||||||
ajcdbv(i)=ajphiv(i)*ratjbv(i)
|
|
||||||
ajplv(i)=ajphiv(i)*ratjplv(i)
|
|
||||||
end do
|
|
||||||
|
|
||||||
rhpp=frhopol(rhotpav)
|
|
||||||
rhpj=frhopol(rhotjava)
|
|
||||||
dpdvp=pabs*2.0_wp_/(sqrt(pi)*drhotpav*fdvdrhot(rhpp))
|
|
||||||
ajphip=currt*2.0_wp_/(sqrt(pi)*drhotjava*fdadrhot(rhpj))
|
|
||||||
call ratioj(rhpj,ratjamx,ratjbmx,ratjplmx)
|
|
||||||
|
|
||||||
call profwidth(nd,rtab,dpdv,rhotp,rhopp,dpdvmx,
|
c now dpdv is dP/dV [MW/m^3]
|
||||||
. drhotp,drhopp)
|
c now ajphiv is J_phi=dI/dA [MA/m^2]
|
||||||
if(ieccd.ne.0) then
|
|
||||||
call profwidth(nd,rtab,ajphiv,rhotjfi,rhopfi,ajmxfi,
|
end subroutine spec
|
||||||
. drhotjfi,drhopfi)
|
|
||||||
xps=rhopfi
|
|
||||||
else
|
|
||||||
rhotjfi=0.0_wp_
|
subroutine pec_tab(xxi,ypt,yamp,ii,rtabpsi1,wdpdv,wajphiv)
|
||||||
rhopfi=0.0_wp_
|
c Power and current projected on psi grid - mid points
|
||||||
ajmxfi=0.0_wp_
|
use const_and_precisions, only : wp_,one,zero
|
||||||
ajphip=0.0_wp_
|
use graydata_flags, only : nnd
|
||||||
drhotjfi=0.0_wp_
|
use beamdata, only : nstep
|
||||||
drhopfi=0.0_wp_
|
use utils, only : locatex,intlin
|
||||||
xps=rhopp
|
c arguments
|
||||||
ratjamx=1.0_wp_
|
integer, intent(in) :: ii
|
||||||
ratjbmx=1.0_wp_
|
integer, parameter :: llmx = 21
|
||||||
ratjplmx=1.0_wp_
|
integer, dimension(llmx) ::isev
|
||||||
end if
|
real(wp_), dimension(nstep), intent(in) :: xxi,ypt,yamp
|
||||||
|
real(wp_), dimension(0:nnd), intent(in) :: rtabpsi1
|
||||||
iif1=iiv(1,1)
|
real(wp_), dimension(nnd), intent(out) :: wdpdv,wajphiv
|
||||||
istmx=1
|
c local variables
|
||||||
do i=2,iif1
|
real(wp_) :: ppa1,ppa2,cci1,cci2,dppa,didst,rt1
|
||||||
if(psjki(1,1,i).ge.0.0_wp_) then
|
integer :: i,is,ise0,idecr,iise0,iise,iis,iis1
|
||||||
if(pdjki(1,1,i).gt.pdjki(1,1,i-1)) istmx=i
|
integer :: ind1,ind2,iind,ind,indi,itb1
|
||||||
end if
|
|
||||||
end do
|
isev=0
|
||||||
stmx=istmx*dst
|
ise0=0
|
||||||
|
idecr=-1
|
||||||
pins_02=0.0_wp_
|
is=1
|
||||||
pins_05=0.0_wp_
|
wdpdv=zero
|
||||||
pins_085=0.0_wp_
|
wajphiv=zero
|
||||||
|
do i=1,ii
|
||||||
xrhot=0.2_wp_
|
if(ise0 == 0) then
|
||||||
call locate(rhotv,nd,xrhot,i1)
|
if(xxi(i) < one) then
|
||||||
call intlin(rhotv(i1),pins(i1),rhotv(i1+1),pins(i1+1),
|
ise0=i
|
||||||
. xrhot,pins_02)
|
isev(is)=i-1
|
||||||
xrhot=0.5_wp_
|
is=is+1
|
||||||
call locate(rhotv,nd,xrhot,i1)
|
end if
|
||||||
call intlin(rhotv(i1),pins(i1),rhotv(i1+1),pins(i1+1),
|
|
||||||
. xrhot,pins_05)
|
|
||||||
xrhot=0.85_wp_
|
|
||||||
call locate(rhotv,nd,xrhot,i1)
|
|
||||||
call intlin(rhotv(i1),pins(i1),rhotv(i1+1),pins(i1+1),
|
|
||||||
. xrhot,pins_085)
|
|
||||||
|
|
||||||
else
|
|
||||||
ajmxfi=0.0_wp_
|
|
||||||
ajphip=0.0_wp_
|
|
||||||
dpdvp=0.0_wp_
|
|
||||||
dpdvmx=0.0_wp_
|
|
||||||
rhotjfi=1.0_wp_
|
|
||||||
rhotjav=1.0_wp_
|
|
||||||
rhotjava=1.0_wp_
|
|
||||||
rhotp=1.0_wp_
|
|
||||||
rhotpav=1.0_wp_
|
|
||||||
drhotjfi=0.0_wp_
|
|
||||||
drhotjava=0.0_wp_
|
|
||||||
drhotp=0.0_wp_
|
|
||||||
drhotpav=0.0_wp_
|
|
||||||
ratjamx=1.0_wp_
|
|
||||||
ratjbmx=1.0_wp_
|
|
||||||
taumn=0
|
|
||||||
taumx=0
|
|
||||||
stmx=stf
|
|
||||||
pins_02=0.0_wp_
|
|
||||||
pins_05=0.0_wp_
|
|
||||||
pins_085=0.0_wp_
|
|
||||||
c end of pabstot > 0
|
|
||||||
end if
|
|
||||||
|
|
||||||
c dPdV [MW/m^3], Jcd [MA/m^2]
|
|
||||||
|
|
||||||
if(ieccd.eq.0) currt=0.0_wp_
|
|
||||||
currtka=currt*1.0e3_wp_
|
|
||||||
|
|
||||||
write(6,*)' '
|
|
||||||
write(6,*)'#beta0 alpha0 Icd Pa Jphip dPdVp '//
|
|
||||||
.'rhotj rhotjava rhotp rhotpav drhotjava drhotpav ratjamx '//
|
|
||||||
.'ratjbmx stmx psipol chipol index_rt Jphimx dPdVmx drhotj drhotp'
|
|
||||||
write(6,99) beta0,alpha0,currtka,pabstot,ajphip,dpdvp,
|
|
||||||
. rhotjfi,rhotjava,rhotp,rhotpav,
|
|
||||||
. drhotjava,drhotpav,ratjamx,ratjbmx,stmx,psipol,chipol,
|
|
||||||
. real(index_rt),ajmxfi,dpdvmx,drhotjfi,drhotp
|
|
||||||
|
|
||||||
write(7,99) currtka,pabstot,ajphip,dpdvp,
|
|
||||||
. rhotjfi,rhotjava,rhotp,rhotpav,
|
|
||||||
. drhotjava,drhotpav,ratjamx,ratjbmx,stmx,psipol,chipol,
|
|
||||||
. real(index_rt),ajmxfi,dpdvmx,drhotjfi,drhotp
|
|
||||||
|
|
||||||
do i=1,nd
|
|
||||||
if (ipec.eq.0) then
|
|
||||||
psin=rtab(i)
|
|
||||||
rhop=sqrt(rtab(i))
|
|
||||||
else
|
else
|
||||||
psin=rtab(i)**2
|
if (idecr == -1) then
|
||||||
rhop=rtab(i)
|
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 if
|
||||||
pinsr=0.0_wp_
|
|
||||||
if(pabstot.gt.0) pinsr=pins(i)/pabstot
|
|
||||||
write(48,99) psin,rhotv(i),ajphiv(i),ajcdbv(i),dpdv(i),
|
|
||||||
. currins(i),pins(i),pinsr,real(index_rt)
|
|
||||||
end do
|
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=rtabpsi1(indi)
|
||||||
|
call locatex(xxi,iise,iise0,iise,rt1,itb1)
|
||||||
|
if(itb1 >= iise0 .and. itb1 < iise) then
|
||||||
|
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)
|
||||||
|
dppa=ppa2-ppa1
|
||||||
|
wdpdv(ind)=wdpdv(ind)+dppa
|
||||||
|
didst=cci2-cci1
|
||||||
|
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,dvol,darea,
|
||||||
|
. dpdv,ajphiv,rhotpav,drhotpav,rhotjava,drhotjava)
|
||||||
|
c radial average values over power and current density profile
|
||||||
|
use const_and_precisions, only : wp_,one,zero,izero,pi
|
||||||
|
use graydata_flags, only : nnd,ieccd
|
||||||
|
implicit none
|
||||||
|
real(wp_), intent(in) :: pabs,currt
|
||||||
|
real(wp_), intent(out) :: rhotpav,rhotjava
|
||||||
|
real(wp_), intent(out) :: drhotpav,drhotjava
|
||||||
|
real(wp_) :: rhopjava,rhoppav
|
||||||
|
real(wp_) :: dpdvp,dpdvmx,rhotp,drhotp
|
||||||
|
real(wp_) :: ajphip,ajmxfi,rhotjfi,drhotjfi
|
||||||
|
real(wp_) :: ratjamx,ratjbmx,ratjplmx
|
||||||
|
real(wp_), dimension(nnd), intent(in) :: rhot_tab,dvol,darea
|
||||||
|
real(wp_), dimension(nnd), intent(in) :: dpdv,ajphiv
|
||||||
|
real(wp_), external :: frhopol,fdadrhot,fdvdrhot
|
||||||
|
|
||||||
|
integer :: i
|
||||||
|
real(wp_) :: sccsa
|
||||||
|
real(wp_) :: rhotjav,rhot2pav,rhot2java
|
||||||
|
|
||||||
|
rhotpav=zero
|
||||||
|
rhot2pav=zero
|
||||||
|
rhotjav=zero
|
||||||
|
rhotjava=zero
|
||||||
|
rhot2java=zero
|
||||||
|
|
||||||
|
if (pabs > zero) then
|
||||||
|
rhotpav=sum(rhot_tab(1:nnd)*dpdv(1:nnd)*dvol(1:nnd))/pabs
|
||||||
|
rhot2pav=sum(rhot_tab(1:nnd)*rhot_tab(1:nnd)*
|
||||||
|
. dpdv(1:nnd)*dvol(1:nnd))/pabs
|
||||||
|
end if
|
||||||
|
|
||||||
|
if (ieccd /= 0) then
|
||||||
|
if (abs(currt) > zero) then
|
||||||
|
rhotjav=sum(rhot_tab(1:nnd)*ajphiv(1:nnd)*darea(1:nnd))/currt
|
||||||
|
end if
|
||||||
|
sccsa=sum(abs(ajphiv(1:nnd))*darea(1:nnd))
|
||||||
|
if (sccsa > zero) then
|
||||||
|
rhotjava=sum(rhot_tab(1:nnd)*abs(ajphiv(1:nnd))
|
||||||
|
. *darea(1:nnd))/sccsa
|
||||||
|
rhot2java=sum(rhot_tab(1:nnd)*rhot_tab(1:nnd)*
|
||||||
|
. abs(ajphiv(1:nnd))*darea(1:nnd))/sccsa
|
||||||
|
end if
|
||||||
|
end if
|
||||||
|
|
||||||
|
c factor sqrt(8) = 2sqrt(2) to match full width of gaussian profile
|
||||||
|
drhotpav=sqrt(8._wp_*(rhot2pav-rhotpav**2))
|
||||||
|
drhotjava=sqrt(8._wp_*(rhot2java-rhotjava**2))
|
||||||
|
|
||||||
return
|
rhoppav=frhopol(rhotpav)
|
||||||
99 format(30(1x,e12.5))
|
rhopjava=frhopol(rhotjava)
|
||||||
end
|
dpdvp=pabs*2.0_wp_/(sqrt(pi)*drhotpav*fdvdrhot(rhoppav))
|
||||||
c
|
ajphip=currt*2.0_wp_/(sqrt(pi)*drhotjava*fdadrhot(rhopjava))
|
||||||
c
|
call ratioj(rhopjava,ratjamx,ratjbmx,ratjplmx)
|
||||||
c
|
|
||||||
subroutine profwidth(nd,xx,yy,rhotmx,rhopmx,ypk,drhot,drhop)
|
call profwidth(nnd,rhot_tab,dpdv,rhotp,dpdvmx,drhotp)
|
||||||
|
if(ieccd.ne.0) then
|
||||||
|
call profwidth(nnd,rhot_tab,ajphiv,rhotjfi,ajmxfi,drhotjfi)
|
||||||
|
else
|
||||||
|
rhotjfi=0.0d0
|
||||||
|
ajmxfi=0.0d0
|
||||||
|
ajphip=0.0d0
|
||||||
|
drhotjfi=0.0d0
|
||||||
|
ratjamx=1.0d0
|
||||||
|
ratjbmx=1.0d0
|
||||||
|
ratjplmx=1.0d0
|
||||||
|
end if
|
||||||
|
|
||||||
|
end subroutine postproc_profiles
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
subroutine profwidth(nd,xx,yy,xpk,ypk,dxxe)
|
||||||
use const_and_precisions, only : wp_,emn1
|
use const_and_precisions, only : wp_,emn1
|
||||||
use graydata_flags, only : ipec
|
use graydata_flags, only : ipec
|
||||||
use utils, only : locatex, locate, intlin, vmaxmini
|
use utils, only : locatex, locate, intlin, vmaxmini
|
||||||
implicit none
|
implicit none
|
||||||
c arguments
|
c arguments
|
||||||
integer :: nd
|
integer :: nd
|
||||||
real(wp_) :: rhotmx,rhopmx,ypk,drhot,drhop
|
|
||||||
real(wp_), dimension(nd) :: xx,yy
|
real(wp_), dimension(nd) :: xx,yy
|
||||||
|
real(wp_), intent(out) :: xpk,ypk,dxxe
|
||||||
c local variables
|
c local variables
|
||||||
integer :: imn,imx,ipk,ie1,ie2
|
integer :: imn,imx,ipk,ie1,ie2
|
||||||
real(wp_) :: xmn,xmx,ymn,ymx,xpkp,xpkm,xpk,yye,rte1,rte2,
|
real(wp_) :: xmn,xmx,ymn,ymx,xpkp,xpkm,yye,rte1,rte2
|
||||||
. psie1,psie2,rhopmn,rhote2,rhote1,rhotmn
|
real(wp_) :: ypkp,ypkm
|
||||||
c common/external functions/variables
|
|
||||||
real(wp_) :: rhotp,rhotm,ypkp,ypkm
|
|
||||||
real(wp_) :: frhotor
|
|
||||||
c
|
|
||||||
common/jmxmn/rhotp,rhotm,ypkp,ypkm
|
|
||||||
c
|
c
|
||||||
call vmaxmini(yy,nd,ymn,ymx,imn,imx)
|
call vmaxmini(yy,nd,ymn,ymx,imn,imx)
|
||||||
ypk=0.0_wp_
|
ypk=0.0_wp_
|
||||||
@ -6789,47 +6775,9 @@ c
|
|||||||
call locate(yy,nd,yye,ie2)
|
call locate(yy,nd,yye,ie2)
|
||||||
call intlin(yy(ie2),xx(ie2),yy(ie1+2),xx(ie1+2),yye,rte2)
|
call intlin(yy(ie2),xx(ie2),yy(ie1+2),xx(ie1+2),yye,rte2)
|
||||||
end if
|
end if
|
||||||
c
|
dxxe=rte2-rte1
|
||||||
ipk=1
|
if(ymx.ne.0.and.ymn.ne.0) dxxe=-dxxe
|
||||||
if(ymx.ne.0.and.ymn.ne.0) ipk=2
|
|
||||||
c
|
|
||||||
drhop=0.0_wp_
|
|
||||||
drhot=0.0_wp_
|
|
||||||
psie1=0.0_wp_
|
|
||||||
psie2=1.0_wp_
|
|
||||||
rhopmx=1.0_wp_
|
|
||||||
rhopmn=0.0_wp_
|
|
||||||
if (ie1.gt.0.or.ie2.gt.0) then
|
|
||||||
if(ipec.eq.0) then
|
|
||||||
rhopmx=sqrt(xpkp)
|
|
||||||
rhopmn=sqrt(xpkm)
|
|
||||||
psie2=rte2
|
|
||||||
psie1=rte1
|
|
||||||
drhop=sqrt(rte2)-sqrt(rte1)
|
|
||||||
else
|
|
||||||
rhopmx=xpkp
|
|
||||||
rhopmn=xpkm
|
|
||||||
drhop=rte2-rte1
|
|
||||||
psie2=rte2**2
|
|
||||||
psie1=rte1**2
|
|
||||||
end if
|
|
||||||
end if
|
|
||||||
c
|
|
||||||
rhotmx=frhotor(rhopmx**2)
|
|
||||||
rhotmn=frhotor(rhopmn**2)
|
|
||||||
rhote2=frhotor(psie2)
|
|
||||||
rhote1=frhotor(psie1)
|
|
||||||
drhot=rhote2-rhote1
|
|
||||||
c
|
|
||||||
if(ipk.eq.2) then
|
|
||||||
drhop=-drhop
|
|
||||||
drhot=-drhot
|
|
||||||
end if
|
|
||||||
|
|
||||||
ypk=ypkp
|
|
||||||
rhotp=rhotmx
|
|
||||||
rhotm=rhotmn
|
|
||||||
c
|
|
||||||
return
|
return
|
||||||
end
|
end
|
||||||
c
|
c
|
||||||
|
|||||||
Loading…
Reference in New Issue
Block a user