src/gray_core.f90: cleanup
- add some comments - annotate loops - indent comments - remove trailing whitespace - reduce usage of opaque global variables - use Fortran 90 logical operators - use Fortran 2003 array syntax
This commit is contained in:
Michele Guerini Rocco
parent
6010a9361b
commit
f5ab40f54f
GPG Key ID:
BFBAF4C975F76450
@ -12,7 +12,7 @@ contains
|
||||
use gray_params, only : raytracing_parameters
|
||||
use const_and_precisions, only : half,two
|
||||
implicit none
|
||||
type(raytracing_parameters), intent(in) :: rtrparam
|
||||
type(raytracing_parameters), intent(inout) :: rtrparam
|
||||
real(wp_), dimension(:,:), intent(out), pointer :: ywork,ypwork, &
|
||||
gri,psjki,ppabs,ccci
|
||||
real(wp_), dimension(:,:,:), intent(out), pointer :: xc,du1,ggri
|
||||
@ -41,6 +41,7 @@ contains
|
||||
end if
|
||||
nray = (nrayr-1)*nrayth + 1
|
||||
jkray1 = (jray1-2)*nrayth + 2
|
||||
rtrparam%nray = nray
|
||||
|
||||
if(nrayr>1) then
|
||||
twodr2 = two*(rwmax/(nrayr-1))**2
|
||||
|
||||
@ -12,13 +12,12 @@ contains
|
||||
use coreprofiles, only : temp, fzeff
|
||||
use dispersion, only : expinit
|
||||
use gray_params, only : gray_parameters, gray_data, gray_results, &
|
||||
print_parameters, &
|
||||
iwarm, ipec, istpr0, igrad, headw, headl, ipass
|
||||
print_parameters
|
||||
use beams, only : xgygcoeff, launchangles2n
|
||||
use beamdata, only : pweight, rayi2jk
|
||||
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
|
||||
use beamdata, only : init_btr, dealloc_beam
|
||||
use pec, only : pec_init, spec, postproc_profiles, dealloc_pec, &
|
||||
rhop_tab, rhot_tab
|
||||
use utils, only : vmaxmin
|
||||
@ -31,7 +30,7 @@ contains
|
||||
implicit none
|
||||
|
||||
! subroutine arguments
|
||||
type(gray_parameters), intent(in) :: params
|
||||
type(gray_parameters), intent(inout) :: params
|
||||
type(gray_data), intent(in) :: data
|
||||
type(gray_results), intent(out) :: results
|
||||
|
||||
@ -43,7 +42,7 @@ contains
|
||||
|
||||
! local variables
|
||||
real(wp_), parameter :: taucr = 12._wp_, etaucr = exp(-taucr)
|
||||
character, dimension(2), parameter :: mode=(/'O','X'/)
|
||||
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
|
||||
@ -104,13 +103,13 @@ contains
|
||||
tau0, alphaabs0, dids0, ccci0, p0jk, ext, eyt, iiv)
|
||||
|
||||
! Initialise the dispersion module
|
||||
if(iwarm > 1) call expinit
|
||||
if(params%ecrh_cd%iwarm > 1) call expinit
|
||||
|
||||
! Initialise the magsurf_data module
|
||||
call flux_average ! requires frhotor for dadrhot,dvdrhot
|
||||
|
||||
! Initialise the output profiles
|
||||
call pec_init(ipec, rhout)
|
||||
call pec_init(params%output%ipec, rhout)
|
||||
nnd = size(rhop_tab) ! number of radial profile points
|
||||
|
||||
call alloc_multipass(nnd, iwait, iroff, iop, iow, yynext, yypnext, yw0, ypw0, stnext, &
|
||||
@ -156,11 +155,11 @@ contains
|
||||
iroff,yynext,yypnext,yw0,ypw0, &
|
||||
stnext,p0ray,taus,tau1,etau1,cpls,cpl1,lgcpl1,psipv,chipv)
|
||||
|
||||
if(params%raytracing%ipol .eq. 0) then
|
||||
if(params%antenna%iox .eq. 2) then ! only X mode on 1st pass
|
||||
cpl0 = (/zero,one/)
|
||||
if(params%raytracing%ipol == 0) then
|
||||
if(params%antenna%iox == 2) then ! only X mode on 1st pass
|
||||
cpl0 = [zero, one]
|
||||
else ! only O mode on 1st pass
|
||||
cpl0 = (/one,zero/)
|
||||
cpl0 = [one, zero]
|
||||
end if
|
||||
end if
|
||||
|
||||
@ -175,7 +174,7 @@ contains
|
||||
index_rt=0 ! global beam index: 1,O 2,X 1st pass
|
||||
! | | | |
|
||||
call log_debug('pass loop start', mod='gray_core', proc='gray_main') ! 3,O 4,X 5,O 6,X 2nd pass
|
||||
do ip=1,ipass
|
||||
main_loop: do ip=1,params%raytracing%ipass
|
||||
write (msg, '("pass: ",g0)') ip
|
||||
call log_info(msg, mod='gray_core', proc='gray_main')
|
||||
|
||||
@ -184,16 +183,17 @@ contains
|
||||
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
|
||||
if(ip > 1) then
|
||||
du1 = zero
|
||||
gri = zero
|
||||
ggri = zero
|
||||
if(ip.eq.ipass) cpl = (/zero,zero/) ! no successive passes
|
||||
if(ip == params%raytracing%ipass) cpl = [zero, zero] ! no successive passes
|
||||
end if
|
||||
|
||||
! =========== beam loop BEGIN ===========
|
||||
call log_debug('beam loop start', mod='gray_core', proc='gray_main')
|
||||
do ib=1,nbeam_pass
|
||||
|
||||
beam_loop: do ib=1,nbeam_pass
|
||||
|
||||
sox = -sox ! invert mode
|
||||
iox = 3-iox ! O-mode at ip=1,ib=1
|
||||
@ -215,8 +215,8 @@ contains
|
||||
|
||||
call vectinit(psjki,ppabs,ccci,tau0,alphaabs0,dids0,ccci0,iiv)
|
||||
|
||||
if(ip.eq.1) then ! 1st pass
|
||||
igrad_b = igrad ! * input value, igrad_b=0 from 2nd pass
|
||||
if(ip == 1) then ! 1st pass
|
||||
igrad_b = params%raytracing%igrad ! * input value, igrad_b=0 from 2nd pass
|
||||
|
||||
tau1 = zero ! * tau from previous passes
|
||||
etau1 = one
|
||||
@ -224,14 +224,12 @@ contains
|
||||
lgcpl1 = zero
|
||||
p0ray = p0jk ! * initial beam power
|
||||
|
||||
call ic_gb(params%antenna%pos, anv0, ak0, &
|
||||
params%antenna%w(1),params%antenna%w(2), &
|
||||
params%antenna%ri(1),params%antenna%ri(2), &
|
||||
params%antenna%phi(1),params%antenna%phi(2), &
|
||||
yw,ypw,xc,du1,gri,ggri,index_rt) ! * initial conditions
|
||||
call set_pol(yw,bres,sox,psipol,chipol,ext,eyt) ! * initial polarization
|
||||
call ic_gb(params, anv0, ak0, yw, ypw, & ! * initial conditions
|
||||
xc, du1, gri, ggri, index_rt)
|
||||
call set_pol(yw, bres, sox, params%raytracing%ipol, & ! * initial polarization
|
||||
psipol, chipol, ext, eyt)
|
||||
|
||||
do jk=1,nray
|
||||
do jk=1,params%raytracing%nray
|
||||
zzm = yw(3,jk)*0.01_wp_
|
||||
rrm = sqrt(yw(1,jk)*yw(1,jk)+yw(2,jk)*yw(2,jk))*0.01_wp_
|
||||
|
||||
@ -258,16 +256,17 @@ contains
|
||||
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
|
||||
if(params%raytracing%nray>1 .and. all(.not.iwait)) & ! iproj=0 ==> nfilp=8
|
||||
call print_projxyzt(stv,yw,0)
|
||||
|
||||
! ======= propagation loop BEGIN =======
|
||||
call log_debug(' propagation loop start', mod='gray_core', proc='gray_main')
|
||||
|
||||
do i=1,nstep
|
||||
propagation_loop: do i=1,params%raytracing%nstep
|
||||
! advance one step with "frozen" grad(S_I)
|
||||
do jk=1,nray
|
||||
do jk=1,params%raytracing%nray
|
||||
if(iwait(jk)) cycle ! jk ray is waiting for next pass
|
||||
stv(jk) = stv(jk) + dst ! current ray step
|
||||
stv(jk) = stv(jk) + params%raytracing%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
|
||||
@ -278,7 +277,7 @@ contains
|
||||
iopmin = 10
|
||||
|
||||
! =========== rays loop BEGIN ===========
|
||||
do jk=1,nray
|
||||
rays_loop: do jk=1,params%raytracing%nray
|
||||
if(iwait(jk)) cycle ! jk ray is waiting for next pass
|
||||
|
||||
! compute derivatives with updated gradient and local plasma values
|
||||
@ -300,20 +299,20 @@ contains
|
||||
ins_pl = (psinv>=zero .and. psinv<params%profiles%psnbnd) ! in/out plasma?
|
||||
ins_wl = inside(data%equilibrium%rlim, data%equilibrium%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
|
||||
ent_pl = (mod(iop(jk),2) == 0 .and. ins_pl) ! enter plasma
|
||||
ext_pl = (mod(iop(jk),2) == 1 .and. .not.ins_pl) ! exit plasma
|
||||
ent_wl = (mod(iow(jk),2) == 0 .and. ins_wl) ! enter vessel
|
||||
ext_wl = (mod(iow(jk),2) == 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 hasn't entered in plasma yet) => continue current pass
|
||||
if(iop(jk) == 1 .and. ip==1) then ! * 1st entrance on 1st pass (ray hasn't entered in plasma yet) => continue current pass
|
||||
|
||||
if(params%raytracing%ipol .eq. 0) then ! + IF single mode propagation
|
||||
if(params%raytracing%ipol == 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
|
||||
else if(cpl(iox) < 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
|
||||
@ -322,7 +321,7 @@ contains
|
||||
end if
|
||||
cpls(jk,index_rt) = cpl(iox)
|
||||
|
||||
if(jk.eq.1) then
|
||||
if(jk == 1) then
|
||||
write (msg,'(" 1st pass - central ray (",a1,"-mode) c=",g0.4)') &
|
||||
mode(iox), cpl(iox)
|
||||
call log_info(msg, mod='gray_core', proc='gray_main')
|
||||
@ -330,20 +329,20 @@ contains
|
||||
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
|
||||
else if(iop(jk) > 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
|
||||
if(ip < params%raytracing%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
|
||||
stnext(jk,index_rt) = stv(jk) - params%raytracing%dst ! . starting step for next pass = last step
|
||||
|
||||
if(cpl(1).lt.etaucr) then ! . low coupled power for O-mode => de-activate derived rays
|
||||
if(cpl(1) < 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
|
||||
if(cpl(2) < 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
|
||||
@ -352,12 +351,12 @@ contains
|
||||
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
|
||||
if(jk == 1) then ! . polarization angles at plasma boundary for central ray
|
||||
psipv(iO:iO+1) = psipol
|
||||
chipv(iO:iO+1) = chipol
|
||||
end if
|
||||
else ! * 1st entrance on 2nd+ pass (ray hasn't entered in plasma since end of previous pass) => continue current pass
|
||||
cpl = (/zero,zero/)
|
||||
cpl = [zero, zero]
|
||||
end if
|
||||
end if
|
||||
|
||||
@ -370,7 +369,7 @@ contains
|
||||
|
||||
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)
|
||||
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), &
|
||||
@ -381,7 +380,7 @@ contains
|
||||
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
|
||||
if(jk == 1 .and. ip == 1) then ! * 1st pass, polarization angles at reflection for central ray
|
||||
psipv(index_rt) = psipol
|
||||
chipv(index_rt) = chipol
|
||||
end if
|
||||
@ -389,18 +388,18 @@ contains
|
||||
if(ins_pl) then ! * plasma-wall overlapping => wall+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
|
||||
if(ip < params%raytracing%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
|
||||
if(cpl(1) < 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
|
||||
if(cpl(2) < 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
|
||||
@ -409,7 +408,7 @@ contains
|
||||
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
|
||||
if(jk == 1) then ! + polarization angles at plasma boundary for central ray
|
||||
psipv(iO:iO+1) = psipol
|
||||
chipv(iO:iO+1) = chipol
|
||||
end if
|
||||
@ -418,17 +417,20 @@ contains
|
||||
end if
|
||||
|
||||
iopmin = min(iopmin,iop(jk))
|
||||
if(ip.lt.ipass) then ! not last pass
|
||||
if(ip < params%raytracing%ipass) then ! not last pass
|
||||
yw0(:,jk) = yw(:,jk) ! * store current coordinates in case
|
||||
ypw0(:,jk) = ypw(:,jk) ! current pass ends on next step
|
||||
end if
|
||||
|
||||
! compute ECRH&CD if (inside plasma & power available>0 & ray still active)
|
||||
if(ierrn==0 .and. iwarm>0 .and. ins_pl .and. &
|
||||
! Compute ECRH&CD if (inside plasma & power available>0 & ray still active)
|
||||
! Note: power check is τ + τ₀ + log(coupling O/X) < τ_critic
|
||||
|
||||
if(ierrn==0 .and. params%ecrh_cd%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)
|
||||
call alpha_effj(params%ecrh_cd, psinv, xg, yg, dens, tekev, &
|
||||
ak0, bres, derdnm, anpl, anpr, sox, anprre, &
|
||||
anprim, alpha, didp, nharm, nhf, iokhawa, ierrhcd)
|
||||
if(ierrhcd/=0) then
|
||||
error = ior(error,ierrhcd)
|
||||
call print_errhcd(ierrhcd,i,anprre,anprim,alpha)
|
||||
@ -446,67 +448,74 @@ contains
|
||||
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
|
||||
|
||||
! Computation of the ray τ, dP/ds, P(s), dI/ds, I(s)
|
||||
|
||||
! optical depth: τ = ∫α(s)ds using the trapezoid rule
|
||||
tau = tau0(jk) + 0.5_wp_*(alphaabs0(jk) + alpha) * dersdst * params%raytracing%dst
|
||||
|
||||
pow = p0ray(jk) * exp(-tau) ! residual power: P = P₀exp(-τ)
|
||||
ppabs(jk,i) = p0ray(jk) - pow ! absorbed power: P_abs = P₀ - P
|
||||
dids = didp * pow * alpha ! current drive: dI/ds = dI/dP⋅P⋅α
|
||||
|
||||
! current: I = ∫dI/ds⋅ds using the trapezoid rule
|
||||
ccci(jk,i) = ccci0(jk) + 0.5_wp_*(dids0(jk) + dids) * dersdst * params%raytracing%dst
|
||||
|
||||
tau0(jk) = tau
|
||||
alphaabs0(jk) = alpha
|
||||
dids0(jk) = dids
|
||||
ccci0(jk) = ccci(jk,i)
|
||||
|
||||
if(iwait(jk)) then ! copy values from last pass for inactive ray
|
||||
ppabs(jk,i:nstep) = ppabs(jk,i-1)
|
||||
ccci(jk,i:nstep) = ccci(jk,i-1)
|
||||
psjki(jk,i:nstep) = psjki(jk,i-1)
|
||||
ppabs(jk,i:params%raytracing%nstep) = ppabs(jk,i-1)
|
||||
ccci(jk,i:params%raytracing%nstep) = ccci(jk,i-1)
|
||||
psjki(jk,i:params%raytracing%nstep) = psjki(jk,i-1)
|
||||
else
|
||||
call print_output(i,jk,stv(jk),p0ray(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/etau1 [dids normalization] = fraction of p0 coupled to this ray (not including absorption from previous passes)
|
||||
end if
|
||||
|
||||
end do
|
||||
end do rays_loop
|
||||
! ============ rays loop END ============
|
||||
|
||||
if(i==nstep) then ! step limit reached?
|
||||
do jk=1,nray
|
||||
if(i==params%raytracing%nstep) then ! step limit reached?
|
||||
do jk=1,params%raytracing%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)
|
||||
if(mod(i,params%output%istpr) == 0) then
|
||||
if(params%raytracing%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(error,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
|
||||
call vmaxmin(tau1+tau0+lgcpl1,params%raytracing%nray,taumn,taumx) ! test on tau + coupling
|
||||
|
||||
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
|
||||
if(ip < params%raytracing%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 => do not increase istop_pass
|
||||
exit
|
||||
end if
|
||||
end do
|
||||
end do propagation_loop
|
||||
call log_debug(' propagation loop end', mod='gray_core', proc='gray_main')
|
||||
! ======== propagation loop END ========
|
||||
|
||||
! print all ray positions in local reference system
|
||||
if(nray > 1 .and. all(.not.iwait)) call print_projxyzt(stv,yw,1)
|
||||
if(params%raytracing%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
|
||||
if(i>params%raytracing%nstep) i=params%raytracing%nstep
|
||||
pabs_beam = sum(ppabs(:,i))
|
||||
icd_beam = sum(ccci(:,i))
|
||||
call vmaxmin(tau0,nray,taumn,taumx) ! taumn,taumx for print
|
||||
call vmaxmin(tau0,params%raytracing%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, &
|
||||
@ -515,12 +524,12 @@ contains
|
||||
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 .and. iopmin.gt.2) then ! not last pass AND at least one ray re-entered plasma
|
||||
cpl_beam1 = sum(p0ray * exp(-tau0) * cpls(:,iO)/cpl1, MASK=iop.gt.2) / &
|
||||
sum(p0ray * exp(-tau0), MASK=iop.gt.2) ! * average O-mode coupling for next beam (on active rays)
|
||||
if(ip < params%raytracing%ipass .and. iopmin > 2) then ! not last pass AND at least one ray re-entered plasma
|
||||
cpl_beam1 = sum(p0ray * exp(-tau0) * cpls(:,iO)/cpl1, MASK=iop > 2) / &
|
||||
sum(p0ray * exp(-tau0), MASK=iop > 2) ! * average O-mode coupling for next beam (on active rays)
|
||||
cpl_beam2 = one - cpl_beam1 ! * average X-mode coupling for next beam
|
||||
|
||||
if(iop(1).gt.2) then ! * central ray O/X-mode coupling for next beam
|
||||
if(iop(1) > 2) then ! * central ray O/X-mode coupling for next beam
|
||||
cpl_cbeam1 = cpls(1,iO)/cpl1(1)
|
||||
cpl_cbeam2 = one - cpl_cbeam1
|
||||
end if
|
||||
@ -543,7 +552,7 @@ contains
|
||||
write(msg, '(3x,a,g0.3," MW")') 'current drive: I=', icd_beam*1.0e3_wp_
|
||||
call log_info(msg, mod='gray_core', proc='gray_main')
|
||||
|
||||
if(ip < ipass) then
|
||||
if(ip < params%raytracing%ipass) then
|
||||
write (msg,'(3x,a,(g0.4,", ",g0.4))') & ! average coupling for next O/X beams (=0 if no ray re-entered plasma)
|
||||
'next couplings [O,X mode]: c=', cpl_beam1, cpl_beam2
|
||||
call log_info(msg, mod='gray_core', proc='gray_main')
|
||||
@ -566,7 +575,7 @@ contains
|
||||
cpl_beam1,cpl_beam2) ! *print 0D results for current beam
|
||||
! ============ post-proc END ============
|
||||
|
||||
end do
|
||||
end do beam_loop
|
||||
call log_debug('beam loop end', mod='gray_core', proc='gray_main')
|
||||
! ============ beam loop END ============
|
||||
|
||||
@ -587,7 +596,7 @@ contains
|
||||
! ======== cumulative prints END ========
|
||||
|
||||
if(istop_pass == nbeam_pass) exit ! no active beams
|
||||
end do
|
||||
end do main_loop
|
||||
call log_debug('pass loop end', mod='gray_core', proc='gray_main')
|
||||
! ============ main loop END ============
|
||||
|
||||
@ -635,26 +644,30 @@ contains
|
||||
|
||||
|
||||
|
||||
subroutine ic_gb(xv0c,anv0c,ak0,wcsi,weta,rcicsi,rcieta,phiw,phir, &
|
||||
ywrk0,ypwrk0,xc0,du10,gri,ggri,index_rt)
|
||||
subroutine ic_gb(params, anv0c, ak0, 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 : zero,one,pi,half,two,degree,ui=>im
|
||||
use math, only : catand
|
||||
use gray_params, only : idst
|
||||
use gray_params, only : gray_parameters
|
||||
use beamdata, only : nray,nrayr,nrayth,rwmax
|
||||
|
||||
implicit none
|
||||
! arguments
|
||||
integer, intent(in) :: index_rt
|
||||
real(wp_), dimension(3), intent(in) :: xv0c,anv0c
|
||||
|
||||
! subroutine arguments
|
||||
type(gray_parameters), intent(in) :: params
|
||||
real(wp_), dimension(3), intent(in) :: anv0c
|
||||
real(wp_), intent(in) :: ak0
|
||||
real(wp_), intent(in) :: wcsi,weta,rcicsi,rcieta,phiw,phir
|
||||
real(wp_), dimension(6, nray), intent(out) :: ywrk0, ypwrk0
|
||||
real(wp_), dimension(3, nrayth, nrayr), intent(out) :: xc0, du10
|
||||
real(wp_), dimension(3, nray), intent(out) :: gri
|
||||
real(wp_), dimension(3, 3, nray), intent(out) :: ggri
|
||||
real(wp_), dimension(3,nrayth,nrayr), intent(out) :: xc0,du10
|
||||
integer, intent(in) :: index_rt
|
||||
|
||||
! local variables
|
||||
real(wp_), dimension(3) :: xv0c
|
||||
real(wp_) :: wcsi,weta,rcicsi,rcieta,phiw,phir
|
||||
integer :: j,k,jk
|
||||
real(wp_) :: csth,snth,csps,snps,phiwrad,phirrad,csphiw,snphiw,alfak
|
||||
real(wp_) :: wwcsi,wweta,sk,sw,dk,dw,rci1,ww1,rci2,ww2,wwxx,wwyy,wwxy
|
||||
@ -670,6 +683,15 @@ contains
|
||||
complex(wp_) :: sss,ddd,phic,qi1,qi2,tc,ts,qqxx,qqxy,qqyy,dqi1,dqi2
|
||||
complex(wp_) :: dqqxx,dqqyy,dqqxy,d2qi1,d2qi2,d2qqxx,d2qqyy,d2qqxy
|
||||
|
||||
! initial beam parameters
|
||||
xv0c = params%antenna%pos
|
||||
wcsi = params%antenna%w(1)
|
||||
weta = params%antenna%w(2)
|
||||
rcicsi = params%antenna%ri(1)
|
||||
rcieta = params%antenna%ri(2)
|
||||
phiw = params%antenna%phi(1)
|
||||
phir = params%antenna%phi(2)
|
||||
|
||||
csth=anv0c(3)
|
||||
snth=sqrt(one-csth**2)
|
||||
if(snth > zero) then
|
||||
@ -911,16 +933,17 @@ contains
|
||||
ywrk0(5,jk) = any
|
||||
ywrk0(6,jk) = anz
|
||||
|
||||
select case(idst)
|
||||
! integration variable
|
||||
select case(params%raytracing%idst)
|
||||
case(0)
|
||||
! optical path: s
|
||||
denom = an0
|
||||
case(1)
|
||||
! integration variable: c*t
|
||||
! "time": c⋅t
|
||||
denom = one
|
||||
case(2)
|
||||
! integration variable: Sr
|
||||
! real eikonal: S_R
|
||||
denom = an20
|
||||
case default ! idst=0
|
||||
! integration variable: s
|
||||
denom = an0
|
||||
end select
|
||||
ypwrk0(1,jk) = anx/denom
|
||||
ypwrk0(2,jk) = any/denom
|
||||
@ -931,9 +954,10 @@ contains
|
||||
|
||||
ddr = anx**2 + any**2 + anz**2 - an20
|
||||
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,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
|
||||
! st=0, index_rt=1, B=0, N=0, psin=-1, Xg=0, Yg=0, gradXg=0
|
||||
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,index_rt,ddr,ddi,zero,zero,[zero,zero,zero])
|
||||
end do
|
||||
end subroutine ic_gb
|
||||
|
||||
@ -1247,13 +1271,16 @@ contains
|
||||
use gray_params, only : iequil
|
||||
use equilibrium, only : psia, equinum_fpol, equinum_psi, equian, sgnbphi
|
||||
use coreprofiles, only : density
|
||||
|
||||
implicit none
|
||||
! arguments
|
||||
|
||||
! subroutine arguments
|
||||
real(wp_), dimension(3), intent(in) :: xv
|
||||
real(wp_), intent(in) :: xgcn,bres
|
||||
real(wp_), intent(out) :: psinv,dens,btot,xg,yg
|
||||
real(wp_), dimension(3), intent(out) :: bv,derxg,deryg
|
||||
real(wp_), dimension(3,3), intent(out) :: derbv
|
||||
|
||||
! local variables
|
||||
integer :: jv
|
||||
real(wp_) :: xx,yy,zz
|
||||
@ -1326,9 +1353,7 @@ contains
|
||||
bzz = dpsidr/rrm
|
||||
|
||||
! bvc(i) = B_i in cylindrical coordinates
|
||||
bvc(1) = brr
|
||||
bvc(2) = bphi
|
||||
bvc(3) = bzz
|
||||
bvc = [brr, bphi, bzz]
|
||||
|
||||
! bv(i) = B_i in cartesian coordinates
|
||||
bv(1)=bvc(1)*csphi - bvc(2)*snphi
|
||||
@ -1394,12 +1419,15 @@ contains
|
||||
|
||||
|
||||
|
||||
subroutine disp_deriv(anv,sox,xg,yg,derxg,deryg,bv,derbv,gr2,dgr2,dgr,ddgr, &
|
||||
dery,anpl,anpr,ddr,ddi,dersdst,derdnm,igrad)
|
||||
subroutine disp_deriv(anv, sox, xg, yg, derxg, deryg, bv, derbv, &
|
||||
gr2, dgr2, dgr, ddgr, dery, anpl, anpr, &
|
||||
ddr, ddi, dersdst, derdnm, igrad)
|
||||
use const_and_precisions, only : zero, one, half, two
|
||||
use gray_params, only : idst
|
||||
|
||||
implicit none
|
||||
! arguments
|
||||
|
||||
! subroutine arguments
|
||||
real(wp_), intent(in) :: xg,yg,gr2,sox
|
||||
real(wp_), intent(out) :: anpl,anpr,ddr,ddi,derdnm,dersdst
|
||||
real(wp_), dimension(3), intent(in) :: anv,bv,derxg,deryg
|
||||
@ -1407,12 +1435,12 @@ contains
|
||||
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
|
||||
real(wp_) :: dan2sdxg, dan2sdyg, ddelnpl2, ddelnpl2x, ddelnpl2y, denom, derdel
|
||||
real(wp_) :: derdom,dfdiadnpl,dfdiadxg,dfdiadyg,fdia,bdotgr !,vgm
|
||||
real(wp_), dimension(3) :: derdxv,danpldxv,derdnv,dbgr !,vgv
|
||||
real(wp_) :: derdom, dfdiadnpl, dfdiadxg, dfdiadyg, fdia, bdotgr
|
||||
real(wp_), dimension(3) :: derdxv, danpldxv, derdnv, dbgr
|
||||
|
||||
an2 = dot_product(anv, anv)
|
||||
anpl = dot_product(anv, bv)
|
||||
@ -1435,17 +1463,28 @@ contains
|
||||
|
||||
duh = one - xg - yg2
|
||||
if(xg > zero) then
|
||||
! Derivatives of the cold plasma dispersion relation
|
||||
!
|
||||
! Λ = N²⊥ - N²s(X,Y,N∥) = 0
|
||||
!
|
||||
! where N²s = 1 - X - XY²⋅(1 + N∥ ∓ √Δ)/[2(1 - X - Y²)]
|
||||
! Δ = (1 - N∥)² + 4N∥⋅(1 - X)/Y²
|
||||
!
|
||||
del = sqrt(dnl**2 + 4.0_wp_*anpl2*(one - xg)/yg2)
|
||||
an2s = one - xg - half*xg*yg2*(one + anpl2 + sox*del)/duh
|
||||
|
||||
! ∂(N²s)/∂X
|
||||
dan2sdxg = - half*yg2*(one - yg2)*(one + anpl2 + sox*del)/duh**2 &
|
||||
+ sox*xg*anpl2/(del*duh) - one
|
||||
! ∂(N²s)/∂Y
|
||||
dan2sdyg = - xg*yg*(one - xg)*(one + anpl2 + sox*del)/duh**2 &
|
||||
+ two*sox*xg*(one - xg)*anpl2/(yg*del*duh)
|
||||
! ∂(N²s)/∂N∥
|
||||
dan2sdnpl = - xg*yg2*anpl/duh &
|
||||
- sox*xg*anpl*(two*(one - xg) - yg2*dnl)/(del*duh)
|
||||
|
||||
if(igrad > 0) then
|
||||
! Derivatives of the eikonal (beamtracing only)
|
||||
ddelnpl2 = two*(two*(one - xg)*(one + 3.0_wp_*anpl2**2) &
|
||||
- yg2*dnl**3)/yg2/del**3
|
||||
fdia = - xg*yg2*(one + half*sox*ddelnpl2)/duh
|
||||
@ -1465,12 +1504,8 @@ contains
|
||||
end if
|
||||
|
||||
bdotgr = dot_product(bv, dgr)
|
||||
do iv=1,3
|
||||
dbgr(iv) = dgr(1)*derbv(1,iv) + bv(1)*ddgr(1,iv) &
|
||||
+ dgr(2)*derbv(2,iv) + bv(2)*ddgr(2,iv) &
|
||||
+ dgr(3)*derbv(3,iv) + bv(3)*ddgr(3,iv)
|
||||
danpldxv(iv) = anv(1)*derbv(1,iv) + anv(2)*derbv(2,iv) + anv(3)*derbv(3,iv)
|
||||
end do
|
||||
dbgr = matmul(dgr, derbv) + matmul(bv, ddgr)
|
||||
danpldxv = matmul(anv, derbv)
|
||||
|
||||
derdxv = -(derxg*dan2sdxg + deryg*dan2sdyg + danpldxv*dan2sdnpl + &
|
||||
igrad*dgr2) &
|
||||
@ -1485,67 +1520,91 @@ contains
|
||||
+ half*yg*dfdiadyg &
|
||||
+ half*anpl*dfdiadnpl)
|
||||
|
||||
! integration variable
|
||||
if (idst == 0) then
|
||||
! integration variable σ
|
||||
select case(idst)
|
||||
case(0)
|
||||
! optical path: s
|
||||
denom = derdnm
|
||||
else if (idst == 1) then
|
||||
! "time": c*t
|
||||
case(1)
|
||||
! "time": c⋅t
|
||||
denom = -derdom
|
||||
else
|
||||
case(2)
|
||||
! real eikonal: S_R
|
||||
denom = dot_product(anv, derdnv)
|
||||
end if
|
||||
end select
|
||||
|
||||
! coefficient for integration in s
|
||||
! ds/dst, where st is the integration variable
|
||||
dersdst = derdnm/denom
|
||||
|
||||
! rhs vector
|
||||
dery(1:3) = derdnv(:)/denom
|
||||
dery(4:6) = -derdxv(:)/denom
|
||||
! r.h.s. vector
|
||||
dery(1:3) = derdnv(:)/denom ! ∂Λ/∂N̅
|
||||
dery(4:6) = -derdxv(:)/denom ! -∂Λ/∂x̅
|
||||
|
||||
! vgv : ~ group velocity
|
||||
! vgm=0
|
||||
! do iv=1,3
|
||||
! vgv(iv)=-derdnv(iv)/derdom
|
||||
! vgm=vgm+vgv(iv)**2
|
||||
! end do
|
||||
! vgm=sqrt(vgm)
|
||||
|
||||
! ddr : dispersion relation (real part)
|
||||
! ddi : dispersion relation (imaginary part)
|
||||
ddr = an2 - an2s - igrad*(gr2 - half*bdotgr**2*fdia)
|
||||
ddi = derdnv(1)*dgr(1) + derdnv(2)*dgr(2) + derdnv(3)*dgr(3)
|
||||
! dispersion relation
|
||||
ddr = an2 - an2s - igrad*(gr2 - half*bdotgr**2*fdia) ! real part
|
||||
ddi = dot_product(derdnv, dgr) ! imaginary part
|
||||
|
||||
end subroutine disp_deriv
|
||||
|
||||
|
||||
|
||||
subroutine alpha_effj(psinv,xg,yg,dens,tekev,ak0,bres,derdnm,anpl,anpr, &
|
||||
sox,anprre,anprim,alpha,didp,nhmin,nhmax,iokhawa,error)
|
||||
subroutine alpha_effj(params, psinv, xg, yg, dens, tekev, ak0, bres, &
|
||||
derdnm, anpl, anpr, sox, anprre, anprim, alpha, &
|
||||
didp, nhmin, nhmax, iokhawa, error)
|
||||
! Computes the absorption coefficient and effective current
|
||||
|
||||
use const_and_precisions, only : zero, pi, mc2=>mc2_
|
||||
use gray_params, only : iwarm,ilarm,ieccd,imx
|
||||
use gray_params, only : ecrh_cd_parameters
|
||||
use coreprofiles, only : fzeff
|
||||
use equilibrium, only : sgnbphi
|
||||
use dispersion, only : harmnumber, warmdisp
|
||||
use eccd, only : setcdcoeff, eccdeff, fjch0, fjch, fjncl
|
||||
use errcodes, only : palph
|
||||
use magsurf_data, only : fluxval
|
||||
|
||||
implicit none
|
||||
! arguments
|
||||
real(wp_),intent(in) ::psinv,ak0,bres
|
||||
real(wp_),intent(in) :: xg,yg,tekev,dens,anpl,anpr,derdnm,sox
|
||||
real(wp_),intent(out) :: anprre,anprim,alpha,didp
|
||||
integer, intent(out) :: nhmin,nhmax,iokhawa
|
||||
integer, intent(out) :: error
|
||||
|
||||
! subroutine arguments
|
||||
|
||||
! Inputs
|
||||
|
||||
! ECRH & CD parameters
|
||||
type(ecrh_cd_parameters) :: params
|
||||
! poloidal flux ψ
|
||||
real(wp_), intent(in) :: psinv
|
||||
! CMA diagram variables: X=(ω_pe/ω)², Y=ω_ce/ω
|
||||
real(wp_), intent(in) :: xg, yg
|
||||
! density [10¹⁹ m⁻³], temperature [keV]
|
||||
real(wp_), intent(in) :: dens, tekev
|
||||
! vacuum wavenumber k₀=ω/c, resonant B field
|
||||
real(wp_), intent(in) :: ak0, bres
|
||||
! group velocity: |∂Λ/∂N̅| where Λ=c²/ω²⋅D_R
|
||||
real(wp_), intent(in) :: derdnm
|
||||
! refractive index: N∥, N⊥ (cold dispersion)
|
||||
real(wp_), intent(in) :: anpl, anpr
|
||||
! sign of polarisation mode: -1 ⇒ O, +1 ⇒ X
|
||||
real(wp_), intent(in) :: sox
|
||||
|
||||
! Outputs
|
||||
|
||||
! real/imaginary parts of N⊥ (warm dispersion)
|
||||
real(wp_), intent(out) :: anprre, anprim
|
||||
! absorption coefficient, current density
|
||||
real(wp_), intent(out) :: alpha, didp
|
||||
! lowest/highest resonant harmonic numbers
|
||||
integer, intent(out) :: nhmin, nhmax
|
||||
! ECCD/overall error codes
|
||||
integer, intent(out) :: iokhawa, error
|
||||
|
||||
! local variables
|
||||
real(wp_) :: rbavi, rrii, rhop
|
||||
integer :: lrm, ithn, ierrcd
|
||||
real(wp_) :: amu,ratiovgr,rbn,rbx
|
||||
real(wp_) :: amu, rbn, rbx
|
||||
real(wp_) :: zeff, cst2, bmxi, bmni, fci
|
||||
real(wp_), dimension(:), pointer :: eccdpar=>null()
|
||||
real(wp_) :: effjcd,effjcdav,akim,btot
|
||||
real(wp_) :: effjcd, effjcdav, btot
|
||||
real(wp_) :: akim
|
||||
complex(wp_) :: ex, ey, ez
|
||||
|
||||
alpha = zero
|
||||
@ -1557,25 +1616,47 @@ contains
|
||||
iokhawa = 0
|
||||
error = 0
|
||||
|
||||
if(tekev>zero) then
|
||||
! absorption computation
|
||||
amu=mc2/tekev
|
||||
call harmnumber(yg,amu,anpl,nhmin,nhmax,iwarm)
|
||||
if(nhmin.gt.0) then
|
||||
lrm=max(ilarm,nhmax)
|
||||
call warmdisp(xg,yg,amu,anpl,anpr,sox,lrm,error,anprre,anprim, &
|
||||
iwarm,imx,ex,ey,ez)
|
||||
akim=ak0*anprim
|
||||
ratiovgr=2.0_wp_*anpr/derdnm!*vgm
|
||||
alpha=2.0_wp_*akim*ratiovgr
|
||||
! Absorption computation
|
||||
|
||||
! Skip when the temperature is zero (no plasma)
|
||||
if (tekev <= zero) return
|
||||
|
||||
! Skip when the lowest resonant harmonic number is zero
|
||||
amu = mc2/tekev ! μ=(mc²/kT)
|
||||
call harmnumber(yg, amu, anpl, nhmin, nhmax, params%iwarm)
|
||||
if (nhmin <= 0) return
|
||||
|
||||
! Solve the dispersion relation
|
||||
lrm = max(params%ilarm, nhmax)
|
||||
call warmdisp(xg, yg, amu, anpl, anpr, sox, lrm, error, &
|
||||
anprre, anprim, params%iwarm, params%imx, &
|
||||
ex, ey, ez)
|
||||
|
||||
! The absoption coefficient is defined as
|
||||
!
|
||||
! α = 2 Im(k̅)⋅s̅
|
||||
!
|
||||
! where s̅ = v̅_g/|v_g|, the direction of the energy flow.
|
||||
! Since v̅_g = - ∂Λ/∂N̅ / ∂Λ/∂ω we have:
|
||||
!
|
||||
! s̅ = ∂Λ/∂N̅ / |∂Λ/∂N̅|
|
||||
! = -2N̅ / |∂Λ/∂N̅| (using the cold dispersion)
|
||||
!
|
||||
! Assuming Im(k∥)=0:
|
||||
!
|
||||
! α = 4 Im(k⊥)⋅N⊥ / |∂Λ/∂N̅|
|
||||
!
|
||||
akim = ak0 * anprim ! imaginary part of k⊥ = k₀N⊥
|
||||
alpha = 4 * akim * anpr/derdnm
|
||||
|
||||
if (alpha < zero) then
|
||||
error = ibset(error, palph)
|
||||
return
|
||||
end if
|
||||
|
||||
! calcolo della efficienza <j/p>: effjcdav [A m/W ]
|
||||
if(ieccd>0) then
|
||||
! current drive computation
|
||||
! Current drive computation
|
||||
if (params%ieccd <= 0) return
|
||||
|
||||
zeff = fzeff(psinv)
|
||||
ithn = 1
|
||||
if (lrm > nhmin) ithn = 2
|
||||
@ -1585,46 +1666,50 @@ contains
|
||||
rbn = btot/bmni
|
||||
rbx = btot/bmxi
|
||||
|
||||
select case(ieccd)
|
||||
select case(params%ieccd)
|
||||
case(1)
|
||||
! cohen model
|
||||
! Cohen model
|
||||
call setcdcoeff(zeff,rbn,rbx,cst2,eccdpar)
|
||||
call eccdeff(yg,anpl,anprre,dens,amu,ex,ey,ez,nhmin,nhmax, &
|
||||
ithn,cst2,fjch,eccdpar,effjcd,iokhawa,ierrcd)
|
||||
case(2)
|
||||
! no trapping
|
||||
! No trapping
|
||||
call setcdcoeff(zeff,cst2,eccdpar)
|
||||
call eccdeff(yg,anpl,anprre,dens,amu,ex,ey,ez,nhmin,nhmax, &
|
||||
ithn,cst2,fjch0,eccdpar,effjcd,iokhawa,ierrcd)
|
||||
case default
|
||||
! neoclassical model
|
||||
! Neoclassical model
|
||||
call setcdcoeff(zeff,rbx,fci,amu,rhop,cst2,eccdpar)
|
||||
call eccdeff(yg,anpl,anprre,dens,amu,ex,ey,ez,nhmin,nhmax, &
|
||||
ithn,cst2,fjncl,eccdpar,effjcd,iokhawa,ierrcd)
|
||||
end select
|
||||
error = error + ierrcd
|
||||
if (associated(eccdpar)) deallocate(eccdpar)
|
||||
|
||||
! current drive efficiency <j/p> [A⋅m/W]
|
||||
effjcdav = rbavi*effjcd
|
||||
didp = sgnbphi * effjcdav / (2.0_wp_*pi*rrii)
|
||||
end if
|
||||
end if
|
||||
end if
|
||||
|
||||
end subroutine alpha_effj
|
||||
|
||||
|
||||
|
||||
subroutine set_pol(ywrk0,bres,sox,psipol0,chipol0,ext0,eyt0)
|
||||
subroutine set_pol(ywrk0, bres, sox, ipol, psipol0, chipol0, ext0, eyt0)
|
||||
use const_and_precisions, only : degree, zero, one, half, im
|
||||
use beamdata, only : nray,nrayth
|
||||
use equilibrium, only : bfield
|
||||
use gray_params, only : ipol
|
||||
use polarization, only : pol_limit, polellipse, stokes_ce, stokes_ell
|
||||
use polarization, only : pol_limit, polellipse, &
|
||||
stokes_ce, stokes_ell
|
||||
|
||||
implicit none
|
||||
! arguments
|
||||
|
||||
! subroutine arguments
|
||||
real(wp_), dimension(6, nray), intent(in) :: ywrk0
|
||||
real(wp_), intent(in) :: sox,bres
|
||||
real(wp_), intent(in) :: bres, sox
|
||||
integer, intent(in) :: ipol
|
||||
real(wp_), intent(inout) :: psipol0, chipol0
|
||||
complex(wp_), dimension(nray), intent(out) :: ext0, eyt0
|
||||
|
||||
! local variables
|
||||
integer :: j,k,jk
|
||||
real(wp_), dimension(3) :: xmv, anv, bv
|
||||
@ -1697,7 +1782,7 @@ contains
|
||||
|
||||
px = 0.5_wp_
|
||||
|
||||
a = reshape(matr2dgrid,(/nr*nz/))
|
||||
a = reshape(matr2dgrid, [nr*nz])
|
||||
|
||||
rcon = 0.0_wp_
|
||||
zcon = 0.0_wp_
|
||||
@ -1842,7 +1927,7 @@ bb: do
|
||||
itm=itm+1
|
||||
inext= i
|
||||
if(jfor/=0) exit aa
|
||||
if(itm .gt. 3) then
|
||||
if(itm > 3) then
|
||||
flag1=.true.
|
||||
exit aa
|
||||
end if
|
||||
@ -1853,7 +1938,7 @@ bb: do
|
||||
|
||||
if(.not.flag1) then
|
||||
lx(in)=0
|
||||
if(itm .eq. 1) exit
|
||||
if(itm == 1) exit
|
||||
end if
|
||||
|
||||
jfor=jx
|
||||
|
||||
@ -51,7 +51,7 @@ module gray_params
|
||||
type raytracing_parameters
|
||||
real(wp_) :: dst ! Integration step size
|
||||
real(wp_) :: rwmax ! Normalized maximum radius of beam power
|
||||
integer :: nrayr, nrayth ! Radial/angular number of rays
|
||||
integer :: nray, nrayr, nrayth ! Total/radial/angular number of rays
|
||||
integer :: igrad ! Complex eikonal switch
|
||||
integer :: nstep ! Max number of integration steps
|
||||
integer :: idst ! Choice of the integration variable
|
||||
|
||||
@ -408,8 +408,7 @@ contains
|
||||
use const_and_precisions, only : zero, degree
|
||||
use coreprofiles, only : set_prfan, set_prfspl, temp, fzeff
|
||||
use dispersion, only : expinit
|
||||
use gray_params, only : gray_parameters, print_parameters, &
|
||||
headw, headl
|
||||
use gray_params, only : gray_parameters, print_parameters
|
||||
use beams, only : launchangles2n, xgygcoeff
|
||||
use magsurf_data, only : flux_average, dealloc_surfvec
|
||||
use beamdata, only : init_btr, dealloc_beam
|
||||
@ -421,7 +420,7 @@ contains
|
||||
implicit none
|
||||
|
||||
! subroutine arguments
|
||||
type(gray_parameters), intent(in) :: params
|
||||
type(gray_parameters), intent(inout) :: params
|
||||
real(wp_), intent(in) :: pabs, icd
|
||||
real(wp_), dimension(:), intent(in) :: jphi, jcd, dpdv, currins, pins
|
||||
real(wp_), intent(out) :: jphip, dpdvp, rhotj, rhotjava, &
|
||||
|
||||
Loading…
Reference in New Issue
Block a user