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

trunk: added multiple passes calculation; added multipass module; graycore/gray_main: subroutine rewritten; equilibrium/read_equil_an: limiter coordinates read from equilibrium file for ipass>1; err_codes/check_err: istop=1 only for alpha<0; pec/pec_tab: fixed index assignements in loops; polarization/pol_limit: fixed ext,eyt normalization

Browse Source
This commit is contained in:
Daniele Micheletti 2019-03-26 14:21:22 +00:00
parent 7a8b66d726
commit 27f1793f14
9 changed files with 832 additions and 205 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
Makefile
View File

@ -5,8 +5,8 @@ EXE=gray
MAINOBJ=main.o
OTHOBJ= beamdata.o beams.o conical.o const_and_precisions.o coreprofiles.o \
dierckx.o dispersion.o eccd.o eierf.o errcodes.o graycore.o \
gray_params.o equilibrium.o limiter.o magsurf_data.o math.o minpack.o numint.o \
pec.o polarization.o quadpack.o reflections.o simplespline.o units.o utils.o
gray_params.o equilibrium.o limiter.o magsurf_data.o math.o minpack.o multipass.o \
numint.o pec.o polarization.o quadpack.o reflections.o simplespline.o units.o utils.o
# Alternative search paths
vpath %.f90 src
@ -15,7 +15,7 @@ vpath %.f src
# Fortran compiler name and flags
FC=gfortran
FFLAGS=-O3
#FFLAGS=-Wall -g -finit-real=nan -ffpe-trap=invalid -fcheck=all -fbounds-check
# FFLAGS=-O0 -Wall -g -finit-real=nan -ffpe-trap=invalid -fcheck=all -fbounds-check
DIRECTIVES = -DREVISION="'$(shell svnversion src)'"
@ -30,7 +30,7 @@ main.o: const_and_precisions.o beams.o coreprofiles.o equilibrium.o \
graycore.o gray_params.o reflections.o
graycore.o: const_and_precisions.o beamdata.o beams.o coreprofiles.o \
dispersion.o eccd.o equilibrium.o errcodes.o gray_params.o \
pec.o polarization.o limiter.o units.o utils.o
pec.o polarization.o limiter.o units.o utils.o reflections.o multipass.o
beams.o: const_and_precisions.o dierckx.o reflections.o simplespline.o utils.o
beamdata.o: const_and_precisions.o gray_params.o
conical.o: const_and_precisions.o
@ -49,6 +49,8 @@ magsurf_data.o: const_and_precisions.o gray_params.o equilibrium.o dierckx.o \
reflections.o simplespline.o units.o utils.o
math.o: const_and_precisions.o
minpack.o: const_and_precisions.o
multipass.o: const_and_precisions.o beamdata.o equilibrium.o gray_params.o \
polarization.o reflections.o
numint.o: const_and_precisions.o
pec.o: const_and_precisions.o beamdata.o equilibrium.o gray_params.o \
magsurf_data.o utils.o

25
src/equilibrium.f90
View File

@ -13,7 +13,7 @@ module equilibrium
! === 2d spline psi(r,z), normalization and derivatives ==========
integer, save :: nsr, nsz
real(wp_), save :: psia, psiant, psinop
real(wp_), save :: psia, psiant, psinop, psnbd
real(wp_), dimension(:), allocatable, save :: tr,tz
real(wp_), dimension(:), allocatable, save :: cceq, cceq01, cceq10, &
cceq20, cceq02, cceq11
@ -151,15 +151,16 @@ contains
subroutine read_equil_an(filenm,rv,zv,fpol,q,unit)
subroutine read_equil_an(filenm,ipass,rv,zv,fpol,q,rlim,zlim,unit)
use utils, only : get_free_unit
implicit none
! arguments
character(len=*), intent(in) :: filenm
integer, intent(in) :: ipass
integer, optional, intent(in) :: unit
real(wp_), dimension(:), allocatable, intent(out) :: rv,zv,fpol,q
real(wp_), dimension(:), allocatable, intent(out) :: rv,zv,fpol,q,rlim,zlim
! local variables
integer :: u
integer :: i, u, nlim
real(wp_) :: rr0m,zr0m,rpam,b0,q0,qa,alq !,rcen,btrcen
if (present(unit)) then
@ -191,6 +192,22 @@ contains
q(1)=q0
q(2)=qa
q(3)=alq
if(ipass.ge.2) then
! get size of boundary and limiter arrays and allocate them
read (u,*) nlim
if (allocated(rlim)) deallocate(rlim)
if (allocated(zlim)) deallocate(zlim)
! store boundary and limiter data
if(nlim>0) then
allocate(rlim(nlim),zlim(nlim))
read(u,*) (rlim(i),zlim(i),i=1,nlim)
end if
end if
close(u)
end subroutine read_equil_an

12
src/errcodes.f90
View File

@ -13,12 +13,16 @@ contains
! arguments
integer, intent(in) :: ierr
integer, intent(inout) :: istop
if(ibits(ierr,pnpl, lnpl )>1 .or. & ! N// too large
ibits(ierr,palph,lalph)==1) then ! alpha < 0
istop = 1
! if(ibits(ierr,pnpl, lnpl )>1 .or. & ! N// too large
! ibits(ierr,palph,lalph)==1) then ! alpha < 0
! istop = 1
! else
! istop = 0
end if
! end if
! if(ibits(ierr,pnpl, lnpl )>1) istop = 1 ! N// too large
if(ibits(ierr,palph,lalph)==1) istop = 1 ! alpha < 0
end subroutine check_err

680
src/graycore.f90
View File

@ -8,16 +8,16 @@ contains
eqp,psrad,terad,derad,zfc,prfp, rlim,zlim, &
p0,fghz,alpha0,beta0,xv0,w1,w2,ri1,ri2,phiw,phir,iox0, &
psipol0,chipol0, dpdv,jcd,pabs,icd, outp,rtrp,hcdp,ierr, rhout)
use const_and_precisions, only : zero, one, degree
use const_and_precisions, only : zero, one, degree, comp_tiny
use coreprofiles, only : set_prfan, set_prfspl, temp, fzeff, unset_prfspl
use dispersion, only : expinit
use gray_params, only : eqparam_type, prfparam_type, outparam_type, &
rtrparam_type, hcdparam_type, antctrl_type, set_codepar, print_params, &
iequil, iprof, iwarm, ipec, istpr0, igrad, headw, headl
iequil, iprof, iwarm, ipec, istpr0, igrad, headw, headl, ipass
use beams, only : read_beam0, read_beam1, launchangles2n, xgygcoeff
use beamdata, only : pweight, rayi2jk
use equilibrium, only : set_equian, set_eqspl, setqphi_num, set_rhospl, &
zbinf, zbsup, unset_eqspl, unset_rhospl, unset_q
zbinf, zbsup, unset_eqspl, unset_rhospl, unset_q, psnbd
use errcodes, only : check_err, print_errn, print_errhcd
use magsurf_data, only : flux_average, dealloc_surfvec
use beamdata, only : init_btr, dealloc_beam, nray, nstep, dst
@ -25,6 +25,10 @@ contains
rhop_tab, rhot_tab
use limiter, only : set_lim, unset_lim
use utils, only : vmaxmin
use reflections, only : inside
use multipass, only : alloc_multipass, dealloc_multipass, initbeam, &
initmultipass, turnoffray, plasma_in, plasma_out, wall_out
implicit none
! arguments
type(eqparam_type), intent(in) :: eqp
@ -49,38 +53,50 @@ contains
integer, intent(out) :: ierr
! local variables
real(wp_), parameter :: taucr = 12._wp_
real(wp_) :: sox,ak0,bres,xgcn,xg,yg,zzm,alpha,didp,anpl,anpr,anprim,anprre
real(wp_) :: chipol,psipol,btot,psinv,dens,tekev,dersdst,derdnm,st,st0
real(wp_), parameter :: taucr = 12._wp_, etaucr = exp(-taucr)
character, dimension(2), parameter :: mode=(/'O','X'/)
real(wp_) :: sox,ak0,bres,xgcn,xg,yg,rrm,zzm,alpha,didp,anpl,anpr,anprim,anprre
real(wp_) :: chipol,psipol,btot,psinv,dens,tekev,dersdst,derdnm
real(wp_) :: tau,pow,dids,ddr,ddi,taumn,taumx
real(wp_) :: rhotpav,drhotpav,rhotjava,drhotjava,dpdvp,jphip
real(wp_) :: rhotp,drhotp,rhotj,drhotj,dpdvmx,jphimx,ratjamx,ratjbmx
real(wp_) :: pabs_beam,icd_beam,cpl_beam1,cpl_beam2,cpl_cbeam1,cpl_cbeam2
real(wp_), dimension(2) :: pabs_pass,icd_pass,cpl,cpl0
real(wp_), dimension(3) :: xv,anv0,anv,bv,derxg
real(wp_), dimension(:,:), pointer :: yw=>null(),ypw=>null(),gri=>null()
real(wp_), dimension(:,:,:), pointer :: xc=>null(),du1=>null(),ggri=>null()
integer :: i,jk,iox,nharm,nhf,nnd,iokhawa,istop,ierrn,ierrhcd,index_rt=1
logical :: ins_pl, somein, allout
integer :: i,j,jk,iox,nharm,nhf,nnd,iokhawa,istop,ierrn,ierrhcd,index_rt
integer :: ip,ib,iopmin,ipar,iO
integer :: igrad_b,ipol,istop_pass,nbeam_pass,nlim
logical :: ins_pl,ins_wl,ent_pl,ext_pl,ent_wl,ext_wl,iboff
real(wp_), dimension(:,:), pointer :: psjki=>null(),ppabs=>null(),ccci=>null()
real(wp_), dimension(:), pointer :: tau0=>null(),alphaabs0=>null(),dids0=>null(), &
ccci0=>null()
real(wp_), dimension(:), pointer :: p0jk=>null()
complex(wp_), dimension(:), pointer :: ext=>null(), eyt=>null()
integer, dimension(:), pointer :: iiv=>null()
real(wp_), dimension(:), allocatable :: jphi,pins,currins
real(wp_), dimension(:,:,:), pointer :: yynext=>null(),yypnext=>null()
real(wp_), dimension(:,:), pointer :: psjki=>null(),ppabs=>null(),ccci=>null()
real(wp_), dimension(:,:), pointer :: taus=>null(),stnext=>null(), &
yw0=>null(),ypw0=>null(),cpls=>null()
real(wp_), dimension(:), pointer :: p0ray=>null(),tau0=>null(),alphaabs0=>null(), &
dids0=>null(),ccci0=>null(),tau1=>null(),etau1=>null(),cpl1=>null(),lgcpl1=>null()
real(wp_), dimension(:), pointer :: p0jk=>null()
real(wp_), dimension(:), pointer :: jphi_beam=>null(),pins_beam=>null(), &
currins_beam=>null(), dpdv_beam=>null(),jcd_beam=>null(),stv=>null(), &
psipv=>null(),chipv=>null()
complex(wp_), dimension(:), pointer :: ext=>null(), eyt=>null()
integer, dimension(:), pointer :: iiv=>null(),iop=>null(),iow=>null()
logical, dimension(:), pointer :: iwait=>null()
logical, dimension(:,:), pointer :: iroff=>null()
! parameters log in file headers
character(len=headw), dimension(headl) :: strheader
type(antctrl_type) :: antp
real(wp_) :: rwall
! ======= set environment BEGIN ======
! ======== set environment BEGIN ========
call set_codepar(eqp,prfp,outp,rtrp,hcdp)
call set_lim(rlim,zlim)
nlim = size(zlim)
if(iequil<2) then
call set_equian(rv(1),zv(1),rv(2), fpol(1)/rv(1), qpsi(1),qpsi(2),qpsi(3))
@ -104,10 +120,16 @@ contains
tau0,alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv)
if(iwarm > 1) call expinit
call pec_init(ipec,rhout)
nnd=size(rhop_tab)
call alloc_multipass(nnd,iwait,iroff,iop,iow,yynext,yypnext,yw0,ypw0,stnext, &
stv,p0ray,taus,tau1,etau1,cpls,cpl1,lgcpl1,jphi_beam, &
pins_beam,currins_beam,dpdv_beam,jcd_beam,psipv,chipv)
! ======= set environment END ======
! ========= set environment END =========
! ======= pre-proc prints BEGIN ======
! ======== pre-proc prints BEGIN ========
antp%alpha=alpha0
antp%beta=beta0
antp%power=p0
@ -118,8 +140,16 @@ contains
antp%ibeam=0
antp%filenm=''
rwall=0._wp_
psnbd=prfp%psnbnd ! plasma boundary
ipol=rtrp%ipol
pabs=zero ! gray_main output
icd=zero
dpdv=zero
jcd=zero
call print_params(rtrp,hcdp,antp,eqp,rwall,prfp,outp,strheader)
call print_headers(strheader)
call print_headers(strheader,0)
! print psi surface for q=1.5 and q=2 on file and psi,rhot,rhop on stdout
call print_surfq((/1.5_wp_,2.0_wp_/))
! print
@ -130,161 +160,435 @@ contains
! print ne, Te, q, Jphi versus psi, rhop, rhot
call print_bres(bres)
call print_prof
call print_maps(bres,xgcn,0.01_wp_*sqrt(xv0(1)**2+xv0(2)**2), &
sin(beta0*degree))
! ======= pre-proc prints END ======
call print_maps(bres,xgcn,0.01_wp_*sqrt(xv0(1)**2+xv0(2)**2),sin(beta0*degree))
! ======= main loop BEGIN ======
iox=iox0
sox=-1.0_wp_
if(iox==2) sox=1.0_wp_
call vectinit(psjki,ppabs,ccci,tau0,alphaabs0,dids0,ccci0,iiv)
call ic_gb(xv0,anv0,ak0,w1,w2,ri1,ri2,phiw,phir,yw,ypw,xc,du1,gri,ggri)
! ========= pre-proc prints END =========
! =========== main loop BEGIN ===========
call initmultipass(ipol,iox0,iroff,yynext,yypnext,yw0,ypw0, &
stnext,p0ray,taus,tau1,etau1,cpls,cpl1,lgcpl1,psipv,chipv)
if(ipol.eq.0) then
if(iox0.eq.2) then ! only X mode on 1st pass
cpl0 = (/zero,one/)
else ! only O mode on 1st pass
cpl0 = (/one,zero/)
end if
end if
sox=one ! mode inverted for each beam
iox=2 ! start with O: sox=-1, iox=1
psipol=psipol0
chipol=chipol0
call set_pol(yw,bres,sox,psipol,chipol,ext,eyt)
call pweight(p0,p0jk)
st0 = zero
if(nray>1) call print_projxyzt(st0,yw,0) ! iproj=0 ==> nfilp=8
somein = .false. ! becomes true if at least part of the beam enters the plasma
! beam/ray propagation
do i=1,nstep
! advance one step with "frozen" grad(S_I)
st=i*dst+st0
do jk=1,nray
call rkstep(sox,bres,xgcn,yw(:,jk),ypw(:,jk),gri(:,jk),ggri(:,:,jk))
end do
! update position and grad
if(igrad == 1) call gradi_upd(yw,ak0,xc,du1,gri,ggri)
allout = .true. ! becomes false if at least part of the beam is inside the plsama
ierr = 0
istop = 0
do jk=1,nray
! compute derivatives with updated gradient and local plasma values
xv = yw(1:3,jk)
anv = yw(4:6,jk)
call ywppla_upd(xv,anv,gri(:,jk),ggri(:,:,jk),sox,bres,xgcn,ypw(:,jk), &
psinv,dens,btot,bv,xg,yg,derxg,anpl,anpr,ddr,ddi,dersdst,derdnm,ierrn)
! update global error code and print message
if (ierrn/=0) then
ierr = ior(ierr,ierrn)
call print_errn(ierrn,i,anpl)
end if
zzm = xv(3)*0.01_wp_
ins_pl = (psinv>=zero .and. psinv<one .and. zzm>=zbinf .and. zzm<=zbsup)
! test if the beam is completely out of the plsama
allout = allout .and. .not.ins_pl
! test if at least part of the beam has entered the plsama
somein = somein .or. ins_pl
! compute ECRH&CD
if(ierrn==0 .and. iwarm>0 .and. ins_pl .and. tau0(jk)<=taucr) then
tekev=temp(psinv)
call alpha_effj(psinv,xg,yg,dens,tekev,ak0,bres,derdnm,anpl,anpr, &
sox,anprre,anprim,alpha,didp,nharm,nhf,iokhawa,ierrhcd)
if (ierrhcd/=0) then
ierr = ior(ierr,ierrhcd)
call print_errhcd(ierrhcd,i,anprre,anprim,alpha)
end if
else
tekev=zero
alpha=zero
didp=zero
anprim=zero
anprre=anpr
nharm=0
nhf=0
iokhawa=0
end if
if(nharm>0) iiv(jk)=i
psjki(jk,i) = psinv
! computation of optical depth tau, dP/ds, P(s), dI/ds, I(s)
tau=tau0(jk)+0.5_wp_*(alphaabs0(jk)+alpha)*dersdst*dst
pow=p0jk(jk)*exp(-tau) !*exp(-tau1v(jk))
ppabs(jk,i)=p0jk(jk)-pow
dids=didp*pow*alpha
ccci(jk,i)=ccci0(jk)+0.5_wp_*(dids0(jk)+dids)*dersdst*dst
tau0(jk)=tau
alphaabs0(jk)=alpha
dids0(jk)=dids
ccci0(jk)=ccci(jk,i)
call print_output(i,jk,st,p0jk(jk)/p0,xv,psinv,btot,bv,ak0,anpl,anpr, &
anv,anprim,dens,tekev,alpha,tau,dids,nharm,nhf,iokhawa,index_rt,ddr,ddi,xg,yg,derxg)
end do
! print ray positions for j=nrayr in local reference system
if (mod(i,istpr0) == 0) then
if(nray > 1) call print_projxyzt(st,yw,0)
nbeam_pass=1 ! max n of beam per pass
index_rt=0 ! global beam index: 1,O 2,X 1st pass
! / \ / \
do ip=1,ipass ! 3,O 4,X 5,O 6,X 2nd pass
pabs_pass = zero
icd_pass = zero
istop_pass = 0 ! stop flag for current pass
nbeam_pass = 2*nbeam_pass ! max n of beams in current pass
if(ip.gt.1) then
du1 = zero ! igrad=0 from 2nd pass
gri = zero
ggri = zero
if(ip.eq.ipass) cpl = (/zero,zero/) ! no successive passes
end if
! =========== beam loop BEGIN ===========
do ib=1,nbeam_pass
sox = -sox ! invert mode
iox = 3-iox ! O-mode at ip=1,ib=1
index_rt = index_rt +1
iO = 2*index_rt +1 ! * index_rt of O-mode derived ray (iX=iO+1)
! check for any error code and stop if necessary
call check_err(ierr,istop)
! test whether further trajectory integration is unnecessary
call vmaxmin(tau0,nray,taumn,taumx)
if ((taumn > taucr) .or. (somein .and. allout)) istop = 1
call initbeam(index_rt,iroff,iboff,iwait,stv,jphi_beam, &
pins_beam,currins_beam,dpdv_beam,jcd_beam)
if(iboff) then ! no propagation for current beam
istop_pass = istop_pass +1 ! * +1 non propagating beam
cycle
end if
call vectinit(psjki,ppabs,ccci,tau0,alphaabs0,dids0,ccci0,iiv)
call print_headers((/' '/),index_rt)
if(ip.eq.1) then ! first pass
igrad_b = igrad
tau1 = zero
etau1 = one
cpl1 = one
lgcpl1 = zero
p0ray = p0jk ! * initial beam power
if(istop == 1) exit
call ic_gb(xv0,anv0,ak0,w1,w2,ri1,ri2,phiw,phir,yw,ypw,xc,du1,gri,ggri,index_rt) ! * initial conditions
call set_pol(yw,bres,sox,psipol,chipol,ext,eyt) ! * initial polarization
do jk=1,nray
if(iwait(jk)) cycle
zzm = yw(3,jk)*0.01_wp_
rrm = sqrt(yw(1,jk)*yw(1,jk)+yw(2,jk)*yw(2,jk))*0.01_wp_
if(inside(rlim,zlim,nlim,rrm,zzm)) then ! * start propagation in/outside vessel?
iow(jk) = 1 ! + inside
else
iow(jk) = 0 ! + outside
end if
end do
else ! successive passes
ipar = (index_rt+1)/2-1 ! * parent beam index
yw = yynext(:,:,ipar) ! * starting coordinates from
ypw = yypnext(:,:,ipar) ! parent beam last step
stv = stnext(:,ipar) ! * starting step from parent beam last step
iow = 1 ! * start propagation inside vessel
tau1 = taus(:,index_rt) ! * tau from previous passes
etau1 = exp(-tau1)
cpl1 = cpls(:,index_rt) ! * coupling from previous passes
lgcpl1 = -log(cpl1)
p0ray = p0jk * etau1 * cpl1 ! * initial beam power
end if
iop = 0 ! start propagation outside plasma
if(nray>1 .and. all(.not.iwait)) call print_projxyzt(stv,yw,0) ! iproj=0 ==> nfilp=8
! ======= propagation loop BEGIN =======
do i=1,nstep
! advance one step with "frozen" grad(S_I)
do jk=1,nray
if(iwait(jk)) cycle ! jk ray is waiting for next pass
stv(jk) = stv(jk) + dst ! current ray step
call rkstep(sox,bres,xgcn,yw(:,jk),ypw(:,jk),gri(:,jk),ggri(:,:,jk),igrad_b)
end do
! update position and grad
if(igrad_b == 1) call gradi_upd(yw,ak0,xc,du1,gri,ggri)
ierr = 0
istop = 0 ! stop flag for current beam
iopmin = 10
! =========== rays loop BEGIN ===========
do jk=1,nray
if(iwait(jk)) cycle ! jk ray is waiting for next pass
! compute derivatives with updated gradient and local plasma values
xv = yw(1:3,jk)
anv = yw(4:6,jk)
call ywppla_upd(xv,anv,gri(:,jk),ggri(:,:,jk),sox,bres,xgcn,ypw(:,jk), &
psinv,dens,btot,bv,xg,yg,derxg,anpl,anpr,ddr,ddi,dersdst,derdnm, &
ierrn,igrad_b)
! update global error code and print message
if(ierrn/=0) then
ierr = ior(ierr,ierrn)
call print_errn(ierrn,i,anpl)
end if
! check entrance/exit plasma/wall
zzm = xv(3)*0.01_wp_
rrm = sqrt(xv(1)*xv(1)+xv(2)*xv(2))*0.01_wp_
ins_pl = (psinv>=zero .and. psinv<psnbd) ! .and. zzm>=zbinf .and. zzm<=zbsup ! in/out plasma?
ins_wl = inside(rlim,zlim,nlim,rrm,zzm) ! in/out vessel?
ent_pl = (mod(iop(jk),2).eq.0 .and. ins_pl) ! enter plasma
ext_pl = (mod(iop(jk),2).eq.1 .and. .not.ins_pl) ! exit plasma
ent_wl = (mod(iow(jk),2).eq.0 .and. ins_wl) ! enter vessel
ext_wl = (mod(iow(jk),2).eq.1 .and. .not.ins_wl) ! exit vessel
if(ent_pl) then ! ray enters plasma
call plasma_in(jk,xv,anv,bres,sox,cpl,psipol,chipol,iop,ext,eyt)
if(iop(jk).eq.1 .and. ip==1) then ! * 1st entrance on 1st pass (ray never entered in plasma) => continue current pass
if(ipol.eq.0) then ! + IF single mode propagation
cpl = cpl0
p0ray(jk) = p0ray(jk)*cpl(iox)
else if(cpl(iox).lt.etaucr) then ! + ELSE IF low coupled power for current mode => de-activate derived rays
call turnoffray(jk,ip+1,2*ib+2-iox,iroff)
iwait(jk) = .true. ! . stop advancement and H&CD computation for current ray
if(cpl(iox).le.comp_tiny) cpl(iox)=etaucr
else ! + ELSE assign coupled power to current ray
p0ray(jk) = p0ray(jk)*cpl(iox)
end if
cpls(jk,index_rt) = cpl(iox)
if(jk.eq.1) then
write(*,*)
write(*,'("1st pass coupling (central ray, ",a1,"-mode)",f9.4)'), &
mode(iox),cpl(iox)
psipv(index_rt) = psipol ! + polarization angles at plasma boundary for central ray
chipv(index_rt) = chipol
end if
else if(iop(jk).gt.2) then ! * 2nd entrance on 1st pass / entrance on 2nd+ pass => end of current pass for ray jk
igrad_b = 0 ! + switch to ray-tracing
iwait(jk) = .true. ! + stop advancement and H&CD computation for current ray
if(ip.lt.ipass) then ! not last pass
yynext(:,jk,index_rt) = yw0(:,jk) ! + copy starting coordinates
yypnext(:,jk,index_rt) = ypw0(:,jk) ! for next pass from last step
stnext(jk,index_rt) = stv(jk) - dst ! + starting step for next pass = last step
if(cpl(1).lt.etaucr) then ! + low coupled power for O-mode => de-activate derived rays
call turnoffray(jk,ip+1,2*ib-1,iroff)
if(cpl(1).le.comp_tiny) cpl(1)=etaucr
end if
if(cpl(2).lt.etaucr) then ! + low coupled power for X-mode => de-activate derived rays
call turnoffray(jk,ip+1,2*ib,iroff)
if(cpl(2).le.comp_tiny) cpl(2)=etaucr
end if
taus(jk,iO:iO+1) = tau1(jk) + tau0(jk) ! + starting tau for next O-mode pass
cpls(jk,iO) = cpl1(jk) * cpl(1) ! + cumulative coupling for next O-mode pass
cpls(jk,iO+1) = cpl1(jk) * cpl(2) ! + cumulative coupling for next X-mode pass
if(jk.eq.1) then ! + polarization angles at plasma boundary for central ray
psipv(iO:iO+1) = psipol
chipv(iO:iO+1) = chipol
end if
else
cpl = (/zero,zero/)
end if
end if
else if(ext_pl) then ! ray exits plasma
call plasma_out(jk,xv,anv,bres,sox,iop,ext,eyt)
end if
if(ent_wl) then ! ray enters vessel
iow(jk)=iow(jk)+1
else if(ext_wl) then ! ray exit vessel
call wall_out(jk,ins_pl,xv,anv,bres,sox,psipol,chipol,iow,iop,ext,eyt)
yw(:,jk) = (/xv,anv/) ! * updated coordinates (reflected)
igrad_b = 0 ! * switch to ray-tracing
call ywppla_upd(xv,anv,gri(:,jk),ggri(:,:,jk),sox,bres,xgcn,ypw(:,jk), &
psinv,dens,btot,bv,xg,yg,derxg,anpl,anpr,ddr,ddi,dersdst,derdnm, &
ierrn,igrad_b) ! * update derivatives after reflection
if(ierrn/=0) then ! * update global error code and print message
ierr = ior(ierr,ierrn)
call print_errn(ierrn,i,anpl)
end if
if(jk.eq.1 .and. ip.eq.1) then ! * 1st pass, polarization angles at reflection for central ray
psipv(index_rt) = psipol
chipv(index_rt) = chipol
end if
if(ins_pl) then ! * plasma-wall overlapping => wall crossing=plasma crossing => end of current pass
iwait(jk) = .true. ! + stop advancement and H&CD computation for current ray
if(ip.lt.ipass) then ! + not last pass
yynext(:,jk,index_rt) = (/xv,anv/) ! . starting coordinates
yypnext(:,jk,index_rt) = ypw(:,jk) ! for next pass = reflection point
stnext(jk,index_rt) = stv(jk) ! . starting step for next pass = step after reflection
call plasma_in(jk,xv,anv,bres,sox,cpl,psipol,chipol,iop,ext,eyt) ! . ray re-enters plasma after reflection
if(cpl(1).lt.etaucr) then ! . low coupled power for O-mode? => de-activate derived rays
call turnoffray(jk,ip+1,2*ib-1,iroff)
if(cpl(1).le.comp_tiny) cpl(1)=etaucr
end if
if(cpl(2).lt.etaucr) then ! . low coupled power for X-mode? => de-activate derived rays
call turnoffray(jk,ip+1,2*ib,iroff)
if(cpl(2).le.comp_tiny) cpl(2)=etaucr
end if
taus(jk,iO:iO+1) = tau1(jk) + tau0(jk) ! . starting tau for next O-mode pass
cpls(jk,iO) = cpl1(jk) * cpl(1) ! . cumulative coupling for next O-mode pass
cpls(jk,iO+1) = cpl1(jk) * cpl(2) ! . cumulative coupling for next X-mode pass
if(jk.eq.1) then ! + polarization angles at plasma boundary for central ray
psipv(iO:iO+1) = psipol
chipv(iO:iO+1) = chipol
end if
end if
end if
end if
iopmin = min(iopmin,iop(jk))
if(ip.lt.ipass) then ! not last pass
yw0(:,jk) = yw(:,jk) ! * store current coordinates
ypw0(:,jk) = ypw(:,jk) ! in case pass ends on next step
end if
! compute ECRH&CD
if(ierrn==0 .and. iwarm>0 .and. ins_pl .and. &
(tau1(jk)+tau0(jk)+lgcpl1(jk))<=taucr .and. .not.iwait(jk)) then ! H&CD computation check
tekev=temp(psinv)
call alpha_effj(psinv,xg,yg,dens,tekev,ak0,bres,derdnm,anpl,anpr, &
sox,anprre,anprim,alpha,didp,nharm,nhf,iokhawa,ierrhcd)
if(ierrhcd/=0) then
ierr = ior(ierr,ierrhcd)
call print_errhcd(ierrhcd,i,anprre,anprim,alpha)
end if
else
tekev=zero
alpha=zero
didp=zero
anprim=zero
anprre=anpr
nharm=0
nhf=0
iokhawa=0
end if
if(nharm>0) iiv(jk)=i
psjki(jk,i) = psinv
! computation of optical depth tau, dP/ds, P(s), dI/ds, I(s)
tau=tau0(jk)+0.5_wp_*(alphaabs0(jk)+alpha)*dersdst*dst
pow=p0ray(jk)*exp(-tau) !*exp(-tau1v(jk))
ppabs(jk,i)=p0ray(jk)-pow
dids=didp*pow*alpha
ccci(jk,i)=ccci0(jk)+0.5_wp_*(dids0(jk)+dids)*dersdst*dst
tau0(jk)=tau
alphaabs0(jk)=alpha
dids0(jk)=dids
ccci0(jk)=ccci(jk,i)
if(iwait(jk)) then
ppabs(jk,i:nstep) = ppabs(jk,i-1)
ccci(jk,i:nstep) = ccci(jk,i-1)
psjki(jk,i:nstep) = psjki(jk,i-1)
else
call print_output(i,jk,stv(jk),p0ray(jk)/etau1(jk),xv,psinv, &
btot,bv,ak0,anpl,anpr,anv,anprim,dens,tekev,alpha,tau,dids, &
nharm,nhf,iokhawa,index_rt,ddr,ddi,xg,yg,derxg) !****** p0ray(jk)/etau1(jk) fattore normalizzazione per dids da controllare
end if
end do
! ============ rays loop END ============
if(i==nstep) then ! step limit reached?
do jk=1,nray
if(iop(jk)<3) call turnoffray(jk,ip,ib,iroff) ! * ray hasn't exited+reentered the plasma by last step => stop ray
end do
end if
! print ray positions for j=nrayr in local reference system
if(mod(i,istpr0) == 0) then
if(nray > 1 .and. all(.not.iwait)) call print_projxyzt(stv,yw,0)
end if
! check for any error code and stop if necessary
call check_err(ierr,istop)
! test whether further trajectory integration is unnecessary
call vmaxmin(tau1+tau0+lgcpl1,nray,taumn,taumx) ! test on tau + coupling
if(taumn > taucr .or. all(iroff(:,index_rt))) istop = 1 ! (residual power~0) or (no ray active) => stop beam
if(istop == 1) then ! stop propagation for current beam
istop_pass = istop_pass +1 ! * +1 non propagating beam
if(ip.lt.ipass) call turnoffray(0,ip,ib,iroff) ! * de-activate derived beams
exit
else if(all(iwait)) then ! all rays in current beam are waiting for next pass
exit
end if
end do
! ======== propagation loop END ========
! print all ray positions in local reference system
if(nray > 1 .and. all(.not.iwait)) call print_projxyzt(stv,yw,1)
! =========== post-proc BEGIN ===========
! compute total absorbed power and driven current for current beam
if(i>nstep) i=nstep
pabs_beam = sum(ppabs(:,i))
icd_beam = sum(ccci(:,i))
call vmaxmin(tau0,nray,taumn,taumx) ! taumn,taumx for print
! compute power and current density profiles for all rays
call spec(psjki,ppabs,ccci,iiv,pabs_beam,icd_beam,dpdv_beam,jphi_beam,jcd_beam, &
pins_beam,currins_beam)
pabs_pass(iox) = pabs_pass(iox) + pabs_beam ! 0D results for current pass, sum on O/X mode beams
icd_pass(iox) = icd_pass(iox) + icd_beam
if(ip.lt.ipass) then
cpl_beam1 = sum(p0ray * exp(-tau0) * cpls(:,iO)/cpl1) / sum(p0ray * exp(-tau0)) ! average O-mode coupling
cpl_beam2 = one - cpl_beam1
cpl_cbeam1 = cpls(1,iO)/cpl1(1)
cpl_cbeam2 = one - cpl_cbeam1
else
cpl_beam1 = zero
cpl_beam2 = zero
end if
! print final results for pass on screen
write(*,*)
write(*,'("End of propagation for beam ",i5," (pass ",i3,", ",a1," mode)")'), &
index_rt,ip,mode(iox)
write(*,'(a,f9.4)') 'final step (s, ct, Sr) = ',stv(1)
write(*,'(a,2e12.5)') 'taumn, taumx = ', taumn,taumx
write(*,'(a,f9.4)') 'Pabs_tot (MW) = ',pabs_beam
write(*,'(a,f9.4)') 'I_tot (kA) = ',icd_beam*1.0e3_wp_
if(ip.lt.ipass) then
write(*,'(a,2(f9.4,1x))') 'Coupling (average, O/X):',cpl_beam1,cpl_beam2 ! average coupling for next O/X beams
write(*,'(a,2(f9.4,1x))') 'Coupling (c. beam, O/X):',cpl_cbeam1,cpl_cbeam2 ! central beam coupling for next O/X beams
end if
call print_pec(rhop_tab,rhot_tab,jphi_beam,jcd_beam,dpdv_beam,currins_beam, &
pins_beam,ip) ! *print power and current density profiles for current beam
call postproc_profiles(pabs_beam,icd_beam,rhot_tab,dpdv_beam,jphi_beam, &
rhotpav,drhotpav,rhotjava,drhotjava,dpdvp,jphip,rhotp,drhotp,rhotj, &
drhotj,dpdvmx,jphimx,ratjamx,ratjbmx) ! *compute profiles width for current beam
call print_finals(pabs_beam,icd_beam,dpdvp,jphip,rhotpav,rhotjava, &
drhotpav,drhotjava,dpdvmx,jphimx,rhotp,rhotj,drhotp,drhotj,ratjamx, &
ratjbmx,stv(1),psipv(index_rt),chipv(index_rt),index_rt,sum(p0ray), &
cpl_beam1,cpl_beam2) ! *print 0D results for current beam
! ============ post-proc END ============
end do
! ============ beam loop END ============
! ======= cumulative prints BEGIN =======
pabs = pabs + sum(pabs_pass) ! *final results (O+X) [gray_main output]
icd = icd + sum(icd_pass)
! print final results for pass on screen
write(*,*)
write(*,'("# End of pass ",i3)'),ip
write(*,'(a,f9.4,f9.4)') '# Pabs_tot (MW) [O,X mode] = ',pabs_pass(1),pabs_pass(2)
write(*,'(a,f9.4,f9.4)') '# I_tot (kA) [O,X mode] = ', &
icd_pass(1)*1.0e3_wp_,icd_pass(2)*1.0e3_wp_
! ======== cumulative prints END ========
if(istop_pass == nbeam_pass) exit ! no active beams
end do
! ============ main loop END ============
! compute total absorbed power and driven current
if (i>nstep) i=nstep
pabs = sum(ppabs(:,i))
icd = sum(ccci(:,i))
! ======= main loop END ======
! ======= post-proc BEGIN ======
! print final results on screen
write(*,*)
write(*,'(a,f9.4)') 'final step (s, ct, Sr) = ',st
write(*,'(a,2e12.5)') 'taumn, taumx = ', taumn,taumx
write(*,'(a,f9.4)') 'Pabs_tot (MW) = ',pabs
write(*,'(a,f9.4)') 'I_tot (kA) = ',icd*1.0e3_wp_
! print all ray positions in local reference system
if(nray > 1) call print_projxyzt(st,yw,1)
! compute power and current density profiles for all rays
call pec_init(ipec,rhout)
nnd=size(rhop_tab)
allocate(jphi(nnd),pins(nnd),currins(nnd))
call spec(psjki,ppabs,ccci,iiv,pabs,icd,dpdv,jphi,jcd,pins,currins)
! print power and current density profiles
call print_pec(rhop_tab,rhot_tab,jphi,jcd,dpdv,currins,pins,index_rt)
! compute profiles width
call postproc_profiles(pabs,icd,rhot_tab,dpdv,jphi, &
rhotpav,drhotpav,rhotjava,drhotjava,dpdvp,jphip, &
rhotp,drhotp,rhotj,drhotj,dpdvmx,jphimx,ratjamx,ratjbmx)
! print 0D results
call print_finals(pabs,icd,dpdvp,jphip,rhotpav,rhotjava,drhotpav, &
drhotjava,dpdvmx,jphimx,rhotp,rhotj,drhotp,drhotj,ratjamx,ratjbmx, &
st,psipol,chipol,index_rt)
! ======= post-proc END ======
! ======= free memory BEGIN ======
! print final results on screen
write(*,*)
write(*,'(a)') '## Final results:'
write(*,'(a,f9.4)') '## Pabs_tot (MW) = ',pabs
write(*,'(a,f9.4)') '## I_tot (kA) = ',icd*1.0e3_wp_
! ========== free memory BEGIN ==========
call unset_eqspl
call unset_q
call unset_rhospl
call unset_prfspl
call unset_lim
call dealloc_surfvec
call dealloc_beam(yw,ypw,xc,du1,gri,ggri,psjki,ppabs,ccci, &
tau0,alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv)
call dealloc_beam(yw,ypw,xc,du1,gri,ggri,psjki,ppabs,ccci,tau0, &
alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv)
call dealloc_pec
deallocate(jphi,pins,currins)
! ======= free memory END ======
call dealloc_multipass(iwait,iroff,iop,iow,yynext,yypnext,yw0,ypw0, &
stnext,stv,p0ray,taus,tau1,etau1,cpls,cpl1,lgcpl1,jphi_beam, &
pins_beam,currins_beam,dpdv_beam,jcd_beam,psipv,chipv)
! =========== free memory END ===========
end subroutine gray_main
@ -320,7 +624,7 @@ contains
subroutine ic_gb(xv0c,anv0c,ak0,wcsi,weta,rcicsi,rcieta,phiw,phir, &
ywrk0,ypwrk0,xc0,du10,gri,ggri)
ywrk0,ypwrk0,xc0,du10,gri,ggri,index_rt)
! beam tracing initial conditions igrad=1
! !!!!!! check ray tracing initial conditions igrad=0 !!!!!!
use const_and_precisions, only : wp_,izero,zero,one,pi,half,two,degree,ui=>im
@ -329,6 +633,7 @@ contains
use beamdata, only : nray,nrayr,nrayth,rwmax
implicit none
! arguments
integer, intent(in) :: index_rt
real(wp_), dimension(3), intent(in) :: xv0c,anv0c
real(wp_), intent(in) :: ak0
real(wp_), intent(in) :: wcsi,weta,rcicsi,rcieta,phiw,phir
@ -616,13 +921,13 @@ contains
ddi = 2*(anxt*gxt + anyt*gyt + anzt*gzt)
call print_output(0,jk,zero,one,xc0(:,k,j),-one,zero,(/zero,zero,zero/), &
ak0,zero,zero,(/zero,zero,zero/),zero,zero,zero,zero,zero,zero, &
0,0,0,1,ddr,ddi,zero,zero,(/zero,zero,zero/)) ! st=0, index_rt=1, B=0, N=0, psin=-1, Xg=0, Yg=0, gradXg=0
0,0,0,index_rt,ddr,ddi,zero,zero,(/zero,zero,zero/)) ! st=0, index_rt=1, B=0, N=0, psin=-1, Xg=0, Yg=0, gradXg=0
end do
end subroutine ic_gb
subroutine rkstep(sox,bres,xgcn,y,yp,dgr,ddgr)
subroutine rkstep(sox,bres,xgcn,y,yp,dgr,ddgr,igrad)
! Runge-Kutta integrator
use const_and_precisions, only : wp_
! use gray_params, only : igrad
@ -633,6 +938,7 @@ contains
real(wp_), dimension(6), intent(in) :: yp
real(wp_), dimension(3), intent(in) :: dgr
real(wp_), dimension(3,3), intent(in) :: ddgr
integer, intent(in) :: igrad
real(wp_), dimension(6) :: yy,fk1,fk2,fk3,fk4
real(wp_) :: gr2
@ -645,18 +951,18 @@ contains
fk1 = yp
yy = y + fk1*hh
call rhs(sox,bres,xgcn,yy,gr2,dgr2,dgr,ddgr,fk2)
call rhs(sox,bres,xgcn,yy,gr2,dgr2,dgr,ddgr,fk2,igrad)
yy = y + fk2*hh
call rhs(sox,bres,xgcn,yy,gr2,dgr2,dgr,ddgr,fk3)
call rhs(sox,bres,xgcn,yy,gr2,dgr2,dgr,ddgr,fk3,igrad)
yy = y + fk3*h
call rhs(sox,bres,xgcn,yy,gr2,dgr2,dgr,ddgr,fk4)
call rhs(sox,bres,xgcn,yy,gr2,dgr2,dgr,ddgr,fk4,igrad)
y = y + h6*(fk1 + 2*fk2 + 2*fk3 + fk4)
end subroutine rkstep
subroutine rhs(sox,bres,xgcn,y,gr2,dgr2,dgr,ddgr,dery)
subroutine rhs(sox,bres,xgcn,y,gr2,dgr2,dgr,ddgr,dery,igrad)
! Compute right-hand side terms of the ray equations (dery)
! used in R-K integrator
use const_and_precisions, only : wp_
@ -667,6 +973,7 @@ contains
real(wp_), dimension(3), intent(in) :: dgr2,dgr
real(wp_), dimension(3,3), intent(in) :: ddgr
real(wp_), dimension(6), intent(out) :: dery
integer, intent(in) :: igrad
! local variables
real(wp_) :: psinv,dens,btot,xg,yg,anpl,anpr,ajphi
real(wp_) :: ddr,ddi,dersdst,derdnm
@ -679,13 +986,13 @@ contains
anv = y(4:6)
call disp_deriv(anv,sox,xg,yg,derxg,deryg,bv,derbv,gr2,dgr2,dgr,ddgr, &
dery,anpl,anpr,ddr,ddi,dersdst,derdnm)
dery,anpl,anpr,ddr,ddi,dersdst,derdnm,igrad)
end subroutine rhs
subroutine ywppla_upd(xv,anv,dgr,ddgr,sox,bres,xgcn,dery,psinv,dens,btot, &
bv,xg,yg,derxg,anpl,anpr,ddr,ddi,dersdst,derdnm,ierr)
bv,xg,yg,derxg,anpl,anpr,ddr,ddi,dersdst,derdnm,ierr,igrad)
! Compute right-hand side terms of the ray equations (dery)
! used after full R-K step and grad(S_I) update
use errcodes, only : pnpl
@ -701,6 +1008,7 @@ contains
real(wp_), dimension(3), intent(out) :: bv
integer, intent(out) :: ierr
real(wp_), dimension(3), intent(out) :: derxg
integer, intent(in) :: igrad
! local variables
real(wp_) :: gr2,ajphi
real(wp_), dimension(3) :: dgr2,deryg
@ -711,7 +1019,7 @@ contains
dgr2 = 2*(dgr(1)*ddgr(:,1) + dgr(2)*ddgr(:,2) + dgr(3)*ddgr(:,3))
call plas_deriv(xv,bres,xgcn,psinv,dens,btot,bv,derbv,xg,yg,derxg,deryg,ajphi)
call disp_deriv(anv,sox,xg,yg,derxg,deryg,bv,derbv,gr2,dgr2,dgr,ddgr, &
dery,anpl,anpr,ddr,ddi,dersdst,derdnm)
dery,anpl,anpr,ddr,ddi,dersdst,derdnm,igrad)
ierr=0
if( abs(anpl) > anplth1) then
@ -1082,9 +1390,9 @@ contains
subroutine disp_deriv(anv,sox,xg,yg,derxg,deryg,bv,derbv,gr2,dgr2,dgr,ddgr, &
dery,anpl,anpr,ddr,ddi,dersdst,derdnm)
dery,anpl,anpr,ddr,ddi,dersdst,derdnm,igrad)
use const_and_precisions, only : wp_,zero,one,half,two
use gray_params, only : idst,igrad
use gray_params, only : idst
implicit none
! arguments
real(wp_), intent(in) :: xg,yg,gr2,sox
@ -1093,6 +1401,7 @@ contains
real(wp_), dimension(3), intent(in) :: dgr2,dgr
real(wp_), dimension(3,3), intent(in) :: ddgr,derbv
real(wp_), dimension(6), intent(out) :: dery
integer, intent(in) :: igrad
! local variables
integer :: iv
real(wp_) :: yg2,anpl2,anpr2,del,dnl,duh,dan2sdnpl,an2,an2s
@ -1662,11 +1971,12 @@ bb: do
subroutine print_headers(strheader)
subroutine print_headers(strheader,index_rt)
use units, only : uprj0,uwbm,udisp,ucenr,uoutr,upec,usumm
implicit none
! arguments
character(len=*), dimension(:), intent(in) :: strheader
integer, intent(in) :: index_rt
! local variables
integer :: i,l
@ -1681,17 +1991,20 @@ bb: do
write(upec,'(1x,a)') strheader(i)
write(usumm,'(1x,a)') strheader(i)
end do
if(index_rt.gt.0) return
write(uprj0,'(1x,a)') '#sst j k xt yt zt rt'
write(uprj0+1,'(1x,a)') '#sst j k xt yt zt rt'
write(uwbm,'(1x,a)') '#sst w1 w2'
write(udisp,'(1x,a)') '#sst Dr_Nr Di_Nr'
write(ucenr,'(1x,a)') '#sst R z phi psin rhot ne Te Btot Bx By Bx Nperp Npl '// &
write(ucenr,'(1x,a)') '#sst R z phi psin rhot ne Te Btot Bx By Bz Nperp Npl '// &
'Nx Ny Nz ki alpha tau Pt dIds nhmin nhmax iohkw index_rt ddr Xg Yg dXgdx dXgdy dXgdz'
write(uoutr,'(1x,a)') '#i k sst x y R z psin tau Npl alpha index_rt'
write(upec,'(1x,a)') '#rhop rhot Jphi Jcdb dPdV Icdins Pins index_rt'
write(usumm,'(1x,a)') '#Icd Pa Jphip dPdVp rhotj rhotjava rhotp rhotpav ' // &
'drhotjava drhotpav ratjamx ratjbmx stmx psipol chipol index_rt ' // &
'Jphimx dPdVmx drhotj drhotp'
'drhotjava drhotpav ratjamx ratjbmx stmx psipol chipol index_rt' // &
'Jphimx dPdVmx drhotj drhotp P0 cplO cplX'
end subroutine print_headers
@ -1872,13 +2185,13 @@ bb: do
subroutine print_projxyzt(st,ywrk,iproj)
subroutine print_projxyzt(stv,ywrk,iproj)
use const_and_precisions, only : wp_, comp_huge, zero, one
use beamdata, only : nray, nrayr, nrayth, rayi2jk
use units, only : uprj0,uwbm
implicit none
! arguments
real(wp_), intent(in) :: st
real(wp_), dimension(:), intent(in) :: stv
real(wp_), dimension(:,:), intent(in) :: ywrk
integer, intent(in) :: iproj
! local variables
@ -1921,14 +2234,14 @@ bb: do
jkv=rayi2jk(jk)
if(.not.(iproj==0 .and. jk==1)) &
write(uprj,'(1x,e16.8e3,2i5,4(1x,e16.8e3))') st,jkv,xti,yti,zti,rti
write(uprj,'(1x,e16.8e3,2i5,4(1x,e16.8e3))') stv(jk),jkv,xti,yti,zti,rti
if(iproj==1 .and. jkv(2)==nrayth) write(uprj,*)
if(rti>=rtimx .and. jkv(1)==nrayr) rtimx = rti
if(rti<=rtimn .and. jkv(1)==nrayr) rtimn = rti
end do
write(uprj,*)
write(uwbm,'(3(1x,e16.8e3))') st,rtimn,rtimx
write(uwbm,'(3(1x,e16.8e3))') stv(1),rtimn,rtimx
end subroutine print_projxyzt
@ -2006,24 +2319,27 @@ bb: do
write(upec,'(7(1x,e16.8e3),i5)') rhop_tab(i),rhot_tab(i), &
jphi(i),jcd(i),dpdv(i),currins(i),pins(i),index_rt
end do
write(upec,*) ''
end subroutine print_pec
subroutine print_finals(pabs,icd,dpdvp,jphip,rhotpav,rhotjava,drhotpav, &
drhotjava,dpdvmx,jphimx,rhotp,rhotj,drhotp,drhotj,ratjamx,ratjbmx, &
stmx,psipol,chipol,index_rt)
stmx,psipol,chipol,index_rt,p0,cpl1,cpl2)
use const_and_precisions, only : wp_
use units, only : usumm
implicit none
real(wp_), intent(in) :: pabs,icd,dpdvp,jphip,rhotpav,rhotjava,drhotpav, &
drhotjava,dpdvmx,jphimx,rhotp,rhotj,drhotp,drhotj,ratjamx,ratjbmx, &
stmx,psipol,chipol
stmx,psipol,chipol,p0,cpl1,cpl2
integer, intent(in) :: index_rt
write(usumm,'(15(1x,e12.5),i5,4(1x,e12.5))') icd,pabs,jphip,dpdvp, &
write(usumm,'(15(1x,e12.5),i5,4(1x,e12.5))') icd,pabs,jphip,dpdvp, &
rhotj,rhotjava,rhotp,rhotpav,drhotjava,drhotpav,ratjamx,ratjbmx, &
stmx,psipol,chipol,index_rt,jphimx,dpdvmx,drhotj,drhotp
stmx,psipol,chipol,index_rt,jphimx,dpdvmx,drhotj,drhotp,p0, &
cpl1,cpl2
write(usumm,*) ''
end subroutine print_finals
end module graycore

4
src/main-sum.f90
View File

@ -46,7 +46,7 @@ program main_sum
! ======= read input data BEGIN =======
!------------ equilibrium ------------
if(eqp%iequil<2) then
call read_equil_an(eqp%filenm, rv, zv, fpol, qpsi)
call read_equil_an(eqp%filenm, rtrp%ipass, rv, zv, fpol, qpsi, rlim, zlim)
! psia sign set to give the correct sign to Iphi (COCOS=3: psia<0 for Iphi>0)
psia = sign(one,qpsi(2)*fpol(1))
else
@ -104,7 +104,7 @@ program main_sum
! set simple limiter if not read from EQDSK
! need to clean up...
r0m=sqrt(x0**2+y0**2)*0.01_wp_
dzmx=rtrp%dst*rtrp%nstep*0.01_wp_
dzmx=abs(rtrp%ipass)*rtrp%dst*rtrp%nstep*0.01_wp_
z0m=z0*0.01_wp_
if (.not.allocated(rlim).or.rtrp%ipass<0) then
if (allocated(rlim)) deallocate(rlim)

4
src/main.f90
View File

@ -38,7 +38,7 @@ program main_std
! ======= read input data BEGIN =======
!------------ equilibrium ------------
if(eqp%iequil<2) then
call read_equil_an(eqp%filenm, rv, zv, fpol, qpsi)
call read_equil_an(eqp%filenm, rtrp%ipass, rv, zv, fpol, qpsi, rlim, zlim)
! psia sign set to give the correct sign to Iphi (COCOS=3: psia<0 for Iphi>0)
psia = sign(one,qpsi(2)*fpol(1))
else
@ -96,7 +96,7 @@ program main_std
! set simple limiter if not read from EQDSK
! need to clean up...
r0m=sqrt(x0**2+y0**2)*0.01_wp_
dzmx=rtrp%dst*rtrp%nstep*0.01_wp_
dzmx=abs(rtrp%ipass)*rtrp%dst*rtrp%nstep*0.01_wp_
z0m=z0*0.01_wp_
if (.not.allocated(rlim).or.rtrp%ipass<0) then
if (allocated(rlim)) deallocate(rlim)

280
src/multipass.f90 Executable file
View File

@ -0,0 +1,280 @@
module multipass
use const_and_precisions, only : wp_, zero, half, one, degree, czero
use beamdata, only : dst, nray
use gray_params, only : ipass
use polarization, only : pol_limit, stokes_ce, polellipse
use reflections, only : wall_refl
use equilibrium, only : bfield
implicit none
integer, save :: nbeam_max
integer, save :: nbeam_tot
contains
! ------------------------------
subroutine plasma_in(i,xv,anv,bres,sox,cpl,psipol1,chipol1,iop,ext,eyt) ! ray enters plasma
implicit none
integer, intent(in) :: i ! ray index
real(wp_), dimension(3), intent(in) :: xv,anv
real(wp_), intent(in) :: bres,sox
real(wp_), dimension(2), intent(out) :: cpl ! coupling (O/X)
real(wp_), intent(out) :: psipol1,chipol1
integer, dimension(:), intent(inout), pointer :: iop
complex(wp_), dimension(:), intent(inout), pointer :: ext,eyt
!
real(wp_) :: rm,csphi,snphi,bphi,br,bz
real(wp_) :: qq1,uu1,vv1,qq,uu,vv,powcpl
real(wp_), dimension(3) :: bv,xmv
complex(wp_) :: ext1,eyt1
!
iop(i)=iop(i)+1
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 bfield(rm,xmv(3),bphi,br,bz)
bv(1)=br*csphi-bphi*snphi
bv(2)=br*snphi+bphi*csphi
bv(3)=bz
call pol_limit(anv,bv,bres,sox,ext1,eyt1)
call stokes_ce(ext1,eyt1,qq1,uu1,vv1) ! stokes parameter at plasma entrance
call stokes_ce(ext(i),eyt(i),qq,uu,vv) ! stokes parameter at plasma exit/wall reflection
powcpl = half*(one + vv*vv1+uu*uu1+qq*qq1)
if(sox.eq.-one) then ! incoming mode = O
cpl=(/powcpl,one-powcpl/)
else ! incoming mode = X
cpl=(/one-powcpl,powcpl/)
end if
if(i.eq.1) then
call polellipse(qq1,uu1,vv1,psipol1,chipol1)
psipol1=psipol1/degree ! convert from rad to degree
chipol1=chipol1/degree
else
psipol1 = zero
chipol1 = zero
end if
end subroutine plasma_in
! ------------------------------
subroutine plasma_out(i,xv,anv,bres,sox,iop,ext,eyt) ! ray exits plasma
implicit none
integer, intent(in) :: i ! ray index
real(wp_), dimension(3), intent(in) :: xv,anv
real(wp_), intent(in) :: bres,sox
integer, dimension(:), intent(inout), pointer :: iop
complex(wp_), dimension(:), intent(out), pointer :: ext,eyt
!
real(wp_) :: rm,csphi,snphi,bphi,br,bz
real(wp_), dimension(3) :: bv,xmv
!
iop(i)=iop(i)+1
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 bfield(rm,xmv(3),bphi,br,bz)
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
end subroutine plasma_out
! ------------------------------
! subroutine wall_in(i) ! ray enters vessel
! implicit none
! integer, intent(in) :: i ! ray index
!
! iow(i)=iow(i)+1
! end subroutine wall_in
! ------------------------------
subroutine wall_out(i,ins,xv,anv,bres,sox,psipol1,chipol1,iow,iop,ext,eyt) ! ray exits vessel
implicit none
integer, intent(in) :: i ! ray index
logical, intent(in) :: ins ! inside plasma? (ins=1 plasma/wall overlap)
real(wp_), dimension(3), intent(inout) :: xv,anv
real(wp_), intent(in) :: bres,sox
real(wp_), intent(out) :: psipol1,chipol1
integer, dimension(:), intent(inout), pointer :: iow,iop
complex(wp_), dimension(:), intent(inout), pointer :: ext,eyt
!
integer :: irfl
real(wp_) :: qq1,uu1,vv1
real(wp_), dimension(3) :: xvrfl,anvrfl,walln
complex(wp_) :: ext1,eyt1
!
iow(i)=iow(i)+1
if(ins) call plasma_out(i,xv,anv,bres,sox,iop,ext,eyt)
call wall_refl(xv-dst*anv,anv,ext(i),eyt(i),xvrfl,anvrfl,ext1,eyt1,walln,irfl)
ext(i) = ext1
eyt(i) = eyt1
xv = xvrfl
anv = anvrfl
if(i.eq.1) then ! polarization angles at wall reflection for central ray
call stokes_ce(ext1,eyt1,qq1,uu1,vv1)
call polellipse(qq1,uu1,vv1,psipol1,chipol1)
psipol1=psipol1/degree ! convert from rad to degree
chipol1=chipol1/degree
else
psipol1 = zero
chipol1 = zero
end if
end subroutine wall_out
! ------------------------------
subroutine initbeam(i,iroff,iboff,iwait,stv,jphi_beam,pins_beam,currins_beam, &
dpdv_beam,jcd_beam) ! initialization at beam propagation start
implicit none
integer, intent(in) :: i ! beam index
logical, dimension(:,:), intent(in), pointer :: iroff
logical, intent(out) :: iboff
logical, dimension(:), intent(out), pointer :: iwait
real(wp_), dimension(:), intent(out), pointer :: jphi_beam,pins_beam, &
currins_beam,dpdv_beam,jcd_beam,stv
iboff = .false. ! beam status (F = active, T = inactive)
iwait = iroff(:,i) ! copy ray status for current beam from global ray status
if(all(iwait)) then ! no rays active => stop beam
iboff = .true.
write(*,*)
write(*,'("Beam ",i5," inactive")'),i
else if(.not.all(.not.iwait)) then ! no all rays active
write(*,*)
write(*,'("WARNING: not all rays in beam ",i5," are active")'),i
end if
stv = zero ! starting step
jphi_beam = zero ! 1D beam profiles
pins_beam = zero
currins_beam = zero
dpdv_beam = zero
jcd_beam = zero
end subroutine initbeam
! ------------------------------
subroutine initmultipass(i,iox,iroff,yynext,yypnext,yw0,ypw0,stnext,p0ray, &
taus,tau1,etau1,cpls,cpl1,lgcpl1,psipv,chipv) ! initialization before pass loop
implicit none
integer, intent(in) :: i, iox ! ipol, mode active on 1st pass
logical, dimension(:,:), intent(out), pointer :: iroff
real(wp_), dimension(:), intent(out), pointer :: p0ray,tau1,etau1,cpl1,lgcpl1, &
psipv,chipv
real(wp_), dimension(:,:), intent(out), pointer :: yw0,ypw0,stnext,taus,cpls
real(wp_), dimension(:,:,:), intent(out), pointer :: yynext,yypnext
!
iroff = .false. ! global ray status (F = active, T = inactive)
if(i.eq.0) call turnoffray(0,1,3-iox,iroff) ! ipol = 0 => stop other mode (iox=1/2 -> stop ib=2/1)
yynext = zero ! starting beam coordinates (1)
yypnext = zero ! starting beam coordinates (2)
yw0 = zero ! temporary beam coordinates (1)
ypw0 = zero ! temporary beam coordinates (2)
stnext = zero ! starting beam step
p0ray = zero ! starting beam power
taus = zero ! starting beam tau
tau1 = zero ! total beam tau
etau1 = one
cpls = one
cpl1 = one
lgcpl1 = zero
psipv = zero ! psi polarization angle at vacuum-plasma boundary
chipv = zero ! chi polarization angle at vacuum-plasma boundary
end subroutine initmultipass
! ------------------------------
subroutine turnoffray(jk,ip,ib,iroff) ! turn off ray propagation
implicit none
integer, intent(in) :: jk, ip, ib ! ray (0=all rays), pass, beam indexes
logical, dimension(:,:), intent(out), pointer :: iroff
!
integer :: ipx, i1, i2
!
if(jk==0) then ! stop all rays
do ipx=ip,ipass
i1 = 2**ipx-2+2**(ipx-ip)*(ib-1)+1
i2 = 2**ipx-2+2**(ipx-ip)*ib
iroff(:,i1:i2) = .true.
end do
else ! only stop ray jk
do ipx=ip,ipass
i1 = 2**ipx-2+2**(ipx-ip)*(ib-1)+1
i2 = 2**ipx-2+2**(ipx-ip)*ib
iroff(jk,i1:i2) = .true.
end do
end if
end subroutine turnoffray
! ------------------------------
subroutine alloc_multipass(dim,iwait,iroff,iop,iow,yynext,yypnext,yw0,ypw0,stnext, &
stv,p0ray,taus,tau1,etau1,cpls,cpl1,lgcpl1,jphi_beam, &
pins_beam,currins_beam,dpdv_beam,jcd_beam,psipv,chipv)
implicit none
integer :: dim
logical, dimension(:), intent(out), pointer :: iwait
logical, dimension(:,:), intent(out), pointer :: iroff
integer, dimension(:), intent(out), pointer :: iop,iow
real(wp_), dimension(:), intent(out), pointer :: jphi_beam,pins_beam,currins_beam, &
dpdv_beam,jcd_beam,stv,tau1,etau1,cpl1,lgcpl1,p0ray,psipv,chipv
real(wp_), dimension(:,:), intent(out), pointer :: taus,cpls,stnext,yw0,ypw0
real(wp_), dimension(:,:,:), intent(out), pointer :: yynext,yypnext
call dealloc_multipass(iwait,iroff,iop,iow,yynext,yypnext,yw0,ypw0,stnext,stv, &
p0ray,taus,tau1,etau1,cpls,cpl1,lgcpl1,jphi_beam,pins_beam,currins_beam, &
dpdv_beam,jcd_beam,psipv,chipv)
nbeam_max = 2**ipass ! max n of beams active at a time
nbeam_tot = 2**(ipass+1)-2 ! total n of beams
allocate(iwait(nray),iroff(nray,nbeam_tot),iop(nray),iow(nray), &
yynext(6,nray,nbeam_max-2),yypnext(6,nray,nbeam_max-2), &
yw0(6,nray),ypw0(6,nray),stnext(nray,nbeam_tot),stv(nray), &
p0ray(nray),taus(nray,nbeam_tot),tau1(nray),etau1(nray), &
cpls(nray,nbeam_tot),cpl1(nray),lgcpl1(nray),jphi_beam(dim), &
pins_beam(dim),currins_beam(dim),dpdv_beam(dim),jcd_beam(dim), &
psipv(nbeam_tot),chipv(nbeam_tot))
end subroutine alloc_multipass
! ------------------------------
subroutine dealloc_multipass(iwait,iroff,iop,iow,yynext,yypnext,yw0,ypw0,stnext, &
stv,p0ray,taus,tau1,etau1,cpls,cpl1,lgcpl1,jphi_beam, &
pins_beam,currins_beam,dpdv_beam,jcd_beam,psipv,chipv)
implicit none
logical, dimension(:), intent(out), pointer :: iwait
logical, dimension(:,:), intent(out), pointer :: iroff
integer, dimension(:), intent(out), pointer :: iop,iow
real(wp_), dimension(:), intent(out), pointer :: stv,p0ray,tau1,etau1,cpl1,lgcpl1, &
jphi_beam,pins_beam,currins_beam,dpdv_beam,jcd_beam,psipv,chipv
real(wp_), dimension(:,:), intent(out), pointer :: yw0,ypw0,stnext,taus,cpls
real(wp_), dimension(:,:,:), intent(out), pointer :: yynext,yypnext
if (associated(iwait)) deallocate(iwait)
if (associated(iroff)) deallocate(iroff)
if (associated(iop)) deallocate(iop)
if (associated(iow)) deallocate(iow)
if (associated(yynext)) deallocate(yynext)
if (associated(yypnext)) deallocate(yypnext)
if (associated(yw0)) deallocate(yw0)
if (associated(ypw0)) deallocate(ypw0)
if (associated(stnext)) deallocate(stnext)
if (associated(stv)) deallocate(stv)
if (associated(p0ray)) deallocate(p0ray)
if (associated(taus)) deallocate(taus)
if (associated(tau1)) deallocate(tau1)
if (associated(etau1)) deallocate(etau1)
if (associated(cpls)) deallocate(cpls)
if (associated(cpl1)) deallocate(cpl1)
if (associated(lgcpl1)) deallocate(lgcpl1)
if (associated(jphi_beam)) deallocate(jphi_beam)
if (associated(pins_beam)) deallocate(pins_beam)
if (associated(currins_beam)) deallocate(currins_beam)
if (associated(dpdv_beam)) deallocate(dpdv_beam)
if (associated(jcd_beam)) deallocate(jcd_beam)
if (associated(psipv)) deallocate(psipv)
if (associated(chipv)) deallocate(chipv)
end subroutine dealloc_multipass
end module multipass

19
src/pec.f90
View File

@ -174,7 +174,7 @@ contains
if(ise0 == 0) then
if(xxi(i) < one) then
ise0 = i
isev(is) = i - 1
isev(is) = max(i - 1, 1)
is = is + 1
end if
else
@ -219,9 +219,14 @@ contains
if (idecr == -1) indi = ind - 1
rt1 = xtab1(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)
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
@ -238,7 +243,7 @@ contains
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
use const_and_precisions, only : pi, comp_eps
use gray_params, only : nnd
use equilibrium, only : frhopol
use magsurf_data, only : fluxval
@ -277,8 +282,8 @@ contains
end if
! 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))
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)

3
src/polarization.f90
View File

@ -86,6 +86,9 @@ contains
if (den>zero) then
ext = (ff*csgam - ui*anpl*sngam)/sqrt(den)
eyt = (-ff*sngam - ui*anpl*csgam)/sqrt(den)
den = sqrt(abs(ext)**2+abs(eyt)**2)
ext = ext/den
eyt = eyt/den
else ! only for XM (sox=+1) when N//=0
ext = -ui*sngam
eyt = -ui*csgam