gray/src/gray_core.f90

2548 lines
93 KiB
Fortran
Raw Normal View History

! This modules contains the core GRAY routines
module gray_core
2015-11-18 17:34:33 +01:00
use const_and_precisions, only : wp_
2015-11-18 17:34:33 +01:00
implicit none
contains
subroutine gray_main(params, data, results, error, rhout)
use const_and_precisions, only : zero, one, degree, comp_tiny
use coreprofiles, only : temp, fzeff
use dispersion, only : expinit
use polarization, only : ellipse_to_field
use gray_params, only : gray_parameters, gray_data, gray_results, &
2024-01-30 10:14:24 +01:00
print_parameters, EQ_VACUUM
use beams, only : xgygcoeff, launchangles2n
use beamdata, only : pweight, rayi2jk
use gray_errors, only : is_critical, print_err_raytracing, print_err_ecrh_cd
use magsurf_data, only : flux_average, dealloc_surfvec
use beamdata, only : init_btr
use pec, only : pec_init, spec, postproc_profiles, dealloc_pec, &
rhop_tab, rhot_tab
use utils, only : vmaxmin, inside
use multipass, only : initbeam, initmultipass, turnoffray, &
plasma_in, plasma_out, wall_out
2021-12-18 18:57:38 +01:00
use logger, only : log_info, log_debug
! subroutine arguments
type(gray_parameters), intent(inout) :: params
type(gray_data), intent(in) :: data
type(gray_results), intent(out) :: results
! Predefined grid for the output profiles (optional)
real(wp_), dimension(:), intent(in), optional :: rhout
2015-11-18 17:34:33 +01:00
! Exit code
integer, intent(out) :: error
! local variables
real(wp_), parameter :: taucr = 12._wp_, etaucr = exp(-taucr)
character, dimension(2), parameter :: mode=['O','X']
2021-12-18 18:57:38 +01:00
integer :: sox
real(wp_) :: 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
fix coupling for subsequent beams In situations when multiple beams are traced, either when allowing multiple plasma crossings (raytracing.ipass > 0) or the initial polarisation is mixed (raytracing.ipol == .true.), the couplings of all but the first beam (with least index_rt) were invalid. The bug is due to the re-use of the psipol,chipol variables as the beams are traced sequentially over the beam_loop. For the first beam being traced the psipol,chipol are correctly initialised to the user-defined value and the resulting coupling is correct. However, in each subsequent beam the values were not set to those of the parent beam (or to the user-defined value in the case of the first X mode beam), but to those of the previous beams (current index_rt - 1). This change repurposes the psipv,chipv arrays to store the polarisation of the parent beams, including the initial user-defined value and makes plasma_in always use these to compute the coupling. In addition, in the case the polarisation is not immediately known (i.e. if raytracing.ipol == .false.), this change postpones the computation of the Jones vector (ext, eyt) from the launch point, if the magnetic equilibrium is available, to when the ray actually crosses the plasma boundary. The original code, besides being strictly incorrect, can lead to non-negligible alterations to the coupling. This change also mean: 1. most of the functionality of `set_pol` has been merged with `plasma_in` 2. the polarisation is undefined and the Jones vector is set to the placeholder value [1, 0] till `plasma_im` is called Finally, `polarcold` is removed because it's unused.
2024-03-02 17:32:03 +01:00
integer :: iox,nharm,nhf,nnd,iokhawa,ierrn,ierrhcd, index_rt, parent_index_rt
integer :: ip,ib,iopmin,ipar,child_index_rt
integer :: igrad_b,istop_pass,nbeam_pass,nlim
logical :: ins_pl,ins_wl,ent_pl,ext_pl,ent_wl,ext_wl,iboff
! i: integration step, jk: global ray index
integer :: i, jk
2015-11-18 17:34:33 +01:00
2021-12-18 18:57:38 +01:00
! buffer for formatting log messages
character(256) :: msg
real(wp_), dimension(2) :: pabs_pass,icd_pass,cpl
real(wp_), dimension(3) :: xv,anv0,anv,bv,derxg
associate ( nray => params%raytracing%nray, &
nrayr => params%raytracing%nrayr, &
nrayth => params%raytracing%nrayth, &
nstep => params%raytracing%nstep, &
nbeam_tot => 2**(params%raytracing%ipass+1)-2, &
nbeam_max => 2**(params%raytracing%ipass))
block
! ray variables
real(wp_), dimension(6, nray) :: yw, ypw
real(wp_), dimension(3, nray) :: gri
real(wp_), dimension(3, 3, nray) :: ggri
real(wp_), dimension(3, nrayth, nrayr) :: xc, du1
real(wp_), dimension(nray) :: ccci0, p0jk
real(wp_), dimension(nray) :: tau0, alphaabs0
real(wp_), dimension(nray) :: dids0
integer, dimension(nray) :: iiv
real(wp_), dimension(nray, nstep) :: psjki, ppabs, ccci
complex(wp_), dimension(nray) :: ext, eyt
! multipass variables
logical, dimension(nray) :: iwait
integer, dimension(nray) :: iow, iop
logical, dimension(nray, nbeam_tot) :: iroff
real(wp_), dimension(6, nray, nbeam_max-2) :: yynext, yypnext
real(wp_), dimension(6, nray) :: yw0, ypw0
real(wp_), dimension(nray, nbeam_tot) :: stnext
real(wp_), dimension(nray) :: stv, p0ray
real(wp_), dimension(nray, nbeam_tot) :: taus, cpls
real(wp_), dimension(nray) :: tau1, etau1, cpl1, lgcpl1
real(wp_), dimension(0:nbeam_tot) :: psipv, chipv
! beam variables
real(wp_), dimension(params%output%nrho) :: jphi_beam, pins_beam
real(wp_), dimension(params%output%nrho) :: currins_beam, dpdv_beam, jcd_beam
! ======== set environment BEGIN ========
! Number of limiter contourn points
nlim = size(data%equilibrium%zlim)
2015-11-18 17:34:33 +01:00
! Compute X=ω/ω_ce and Y=(ω/ω_pe)² (with B=1)
call xgygcoeff(params%antenna%fghz, ak0, bres, xgcn)
2015-11-18 17:34:33 +01:00
! Compute the initial cartesian wavevector (anv0)
! from launch angles α,β and the position
call launchangles2n(params%antenna, anv0)
2015-11-18 17:34:33 +01:00
! Initialise the ray variables (beamtracing)
call init_btr(params%raytracing)
2015-11-18 17:34:33 +01:00
! Initialise the dispersion module
if(params%ecrh_cd%iwarm > 1) call expinit
2024-01-30 10:14:24 +01:00
if (params%equilibrium%iequil /= EQ_VACUUM) then
! Initialise the magsurf_data module
call flux_average ! requires frhotor for dadrhot,dvdrhot
2024-01-30 10:14:24 +01:00
! Initialise the output profiles
call pec_init(params%output%ipec, rhout)
nnd = size(rhop_tab) ! number of radial profile points
end if
! Allocate memory for the results...
allocate(results%dpdv(params%output%nrho))
allocate(results%jcd(params%output%nrho))
! ...and initialise them
results%pabs = zero
results%icd = zero
results%dpdv = zero
results%jcd = zero
! ========= set environment END =========
! ======== pre-proc prints BEGIN ========
call print_headers(params)
2021-12-18 18:57:38 +01:00
! print ψ surface for q=1.5 and q=2 on file and log psi,rhot,rhop
call print_surfq([1.5_wp_, 2.0_wp_])
! print initial position
2021-12-18 18:57:38 +01:00
write (msg, '("initial position:",3(x,g0.3))') params%antenna%pos
call log_info(msg, mod='gray_core', proc='gray_main')
write (msg, '("initial direction:",2(x,a,"=",g0.2))') &
'α', params%antenna%alpha, 'β', params%antenna%beta
call log_info(msg, mod='gray_core', proc='gray_main')
! print Btot=Bres
! print ne, Te, q, Jphi versus psi, rhop, rhot
2024-01-30 10:14:24 +01:00
if (params%equilibrium%iequil /= EQ_VACUUM) then
call print_bres(bres)
call print_prof(params%profiles)
call print_maps(bres, xgcn, &
norm2(params%antenna%pos(1:2)) * 0.01_wp_, &
sin(params%antenna%beta*degree))
end if
! ========= pre-proc prints END =========
! =========== main loop BEGIN ===========
call initmultipass(params%raytracing%ipol, params%antenna%iox, &
iroff,yynext,yypnext,yw0,ypw0, &
stnext,p0ray,taus,tau1,etau1,cpls,cpl1,lgcpl1,psipv,chipv)
sox=1 ! mode inverted for each beam
iox=2 ! start with O: sox=-1, iox=1
call pweight(params%antenna%power, p0jk)
2015-11-18 17:34:33 +01:00
fix coupling for subsequent beams In situations when multiple beams are traced, either when allowing multiple plasma crossings (raytracing.ipass > 0) or the initial polarisation is mixed (raytracing.ipol == .true.), the couplings of all but the first beam (with least index_rt) were invalid. The bug is due to the re-use of the psipol,chipol variables as the beams are traced sequentially over the beam_loop. For the first beam being traced the psipol,chipol are correctly initialised to the user-defined value and the resulting coupling is correct. However, in each subsequent beam the values were not set to those of the parent beam (or to the user-defined value in the case of the first X mode beam), but to those of the previous beams (current index_rt - 1). This change repurposes the psipv,chipv arrays to store the polarisation of the parent beams, including the initial user-defined value and makes plasma_in always use these to compute the coupling. In addition, in the case the polarisation is not immediately known (i.e. if raytracing.ipol == .false.), this change postpones the computation of the Jones vector (ext, eyt) from the launch point, if the magnetic equilibrium is available, to when the ray actually crosses the plasma boundary. The original code, besides being strictly incorrect, can lead to non-negligible alterations to the coupling. This change also mean: 1. most of the functionality of `set_pol` has been merged with `plasma_in` 2. the polarisation is undefined and the Jones vector is set to the placeholder value [1, 0] till `plasma_im` is called Finally, `polarcold` is removed because it's unused.
2024-03-02 17:32:03 +01:00
! Set original polarisation
psipv(0) = params%antenna%psi
chipv(0) = params%antenna%chi
nbeam_pass=1 ! max n of beam per pass
index_rt=0 ! global beam index: 1,O 2,X 1st pass
2019-03-28 10:50:28 +01:00
! | | | |
2021-12-18 18:57:38 +01:00
call log_debug('pass loop start', mod='gray_core', proc='gray_main') ! 3,O 4,X 5,O 6,X 2nd pass
main_loop: do ip=1,params%raytracing%ipass
2021-12-18 18:57:38 +01:00
write (msg, '("pass: ",g0)') ip
call log_info(msg, mod='gray_core', proc='gray_main')
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 > 1) then
2019-03-28 10:50:28 +01:00
du1 = zero
gri = zero
ggri = zero
if(ip == params%raytracing%ipass) cpl = [zero, zero] ! no successive passes
end if
! =========== beam loop BEGIN ===========
2021-12-18 18:57:38 +01:00
call log_debug('beam loop start', mod='gray_core', proc='gray_main')
beam_loop: do ib=1,nbeam_pass
sox = -sox ! invert mode
iox = 3-iox ! O-mode at ip=1,ib=1
index_rt = index_rt +1
fix coupling for subsequent beams In situations when multiple beams are traced, either when allowing multiple plasma crossings (raytracing.ipass > 0) or the initial polarisation is mixed (raytracing.ipol == .true.), the couplings of all but the first beam (with least index_rt) were invalid. The bug is due to the re-use of the psipol,chipol variables as the beams are traced sequentially over the beam_loop. For the first beam being traced the psipol,chipol are correctly initialised to the user-defined value and the resulting coupling is correct. However, in each subsequent beam the values were not set to those of the parent beam (or to the user-defined value in the case of the first X mode beam), but to those of the previous beams (current index_rt - 1). This change repurposes the psipv,chipv arrays to store the polarisation of the parent beams, including the initial user-defined value and makes plasma_in always use these to compute the coupling. In addition, in the case the polarisation is not immediately known (i.e. if raytracing.ipol == .false.), this change postpones the computation of the Jones vector (ext, eyt) from the launch point, if the magnetic equilibrium is available, to when the ray actually crosses the plasma boundary. The original code, besides being strictly incorrect, can lead to non-negligible alterations to the coupling. This change also mean: 1. most of the functionality of `set_pol` has been merged with `plasma_in` 2. the polarisation is undefined and the Jones vector is set to the placeholder value [1, 0] till `plasma_im` is called Finally, `polarcold` is removed because it's unused.
2024-03-02 17:32:03 +01:00
child_index_rt = 2*index_rt + 1 ! * index_rt of O-mode child beam
parent_index_rt = ceiling(index_rt / 2.0_wp_) - 1 ! * index_rt of parent beam
2015-11-18 17:34:33 +01:00
call initbeam(index_rt,iroff,iboff,iwait,stv,jphi_beam, &
pins_beam,currins_beam,dpdv_beam,jcd_beam)
2021-12-18 18:57:38 +01:00
write(msg, '(" beam: ",g0," (",a1," mode)")') index_rt, mode(iox)
call log_info(msg, mod='gray_core', proc='gray_main')
if(iboff) then ! no propagation for current beam
istop_pass = istop_pass +1 ! * +1 non propagating beam
2021-12-18 18:57:38 +01:00
call log_info(" beam is off", mod='gray_core', proc='gray_main')
cycle
2015-11-18 17:34:33 +01:00
end if
call vectinit(psjki,ppabs,ccci,tau0,alphaabs0,dids0,ccci0,iiv)
if(ip == 1) then ! 1st pass
igrad_b = params%raytracing%igrad ! * input value, igrad_b=0 from 2nd pass
2019-03-28 10:50:28 +01:00
tau1 = zero ! * tau from previous passes
etau1 = one
2019-03-28 10:50:28 +01:00
cpl1 = one ! * coupling from previous passes
lgcpl1 = zero
p0ray = p0jk ! * initial beam power
fix coupling for subsequent beams In situations when multiple beams are traced, either when allowing multiple plasma crossings (raytracing.ipass > 0) or the initial polarisation is mixed (raytracing.ipol == .true.), the couplings of all but the first beam (with least index_rt) were invalid. The bug is due to the re-use of the psipol,chipol variables as the beams are traced sequentially over the beam_loop. For the first beam being traced the psipol,chipol are correctly initialised to the user-defined value and the resulting coupling is correct. However, in each subsequent beam the values were not set to those of the parent beam (or to the user-defined value in the case of the first X mode beam), but to those of the previous beams (current index_rt - 1). This change repurposes the psipv,chipv arrays to store the polarisation of the parent beams, including the initial user-defined value and makes plasma_in always use these to compute the coupling. In addition, in the case the polarisation is not immediately known (i.e. if raytracing.ipol == .false.), this change postpones the computation of the Jones vector (ext, eyt) from the launch point, if the magnetic equilibrium is available, to when the ray actually crosses the plasma boundary. The original code, besides being strictly incorrect, can lead to non-negligible alterations to the coupling. This change also mean: 1. most of the functionality of `set_pol` has been merged with `plasma_in` 2. the polarisation is undefined and the Jones vector is set to the placeholder value [1, 0] till `plasma_im` is called Finally, `polarcold` is removed because it's unused.
2024-03-02 17:32:03 +01:00
call ic_gb(params, anv0, ak0, yw, ypw, stv, xc, du1, gri, ggri, index_rt) ! * initial conditions
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_
if(inside(data%equilibrium%rlim, data%equilibrium%zlim, &
rrm, zzm)) then ! * start propagation in/outside vessel?
iow(jk) = 1 ! + inside
else
iow(jk) = 0 ! + outside
end if
end do
2019-03-28 10:50:28 +01:00
else ! 2nd+ 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(params%raytracing%nray>1 .and. all(.not.iwait)) & ! iproj=0 ==> nfilp=8
call print_projxyzt(stv,yw,0)
! ======= propagation loop BEGIN =======
2021-12-18 18:57:38 +01:00
call log_debug(' propagation loop start', mod='gray_core', proc='gray_main')
propagation_loop: do i=1,params%raytracing%nstep
! advance one step with "frozen" grad(S_I)
do jk=1,params%raytracing%nray
if(iwait(jk)) cycle ! jk ray is waiting for next pass
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
if(igrad_b == 1) call gradi_upd(yw,ak0,xc,du1,gri,ggri)
error = 0
iopmin = 10
! =========== rays loop BEGIN ===========
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
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
error = ior(error,ierrn)
call print_err_raytracing(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<params%profiles%psnbnd) ! in/out plasma?
ins_wl = inside(data%equilibrium%rlim, data%equilibrium%zlim, &
rrm, zzm) ! in/out 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
2022-05-26 02:33:55 +02:00
write (msg, '(" ray ",g0," entered plasma (",g0," steps)")') jk, i
call log_debug(msg, mod='gray_core', proc='gray_main')
call ellipse_to_field(psipv(parent_index_rt), chipv(parent_index_rt), & ! compute polarisation and couplings
ext, eyt)
fix coupling for subsequent beams In situations when multiple beams are traced, either when allowing multiple plasma crossings (raytracing.ipass > 0) or the initial polarisation is mixed (raytracing.ipol == .true.), the couplings of all but the first beam (with least index_rt) were invalid. The bug is due to the re-use of the psipol,chipol variables as the beams are traced sequentially over the beam_loop. For the first beam being traced the psipol,chipol are correctly initialised to the user-defined value and the resulting coupling is correct. However, in each subsequent beam the values were not set to those of the parent beam (or to the user-defined value in the case of the first X mode beam), but to those of the previous beams (current index_rt - 1). This change repurposes the psipv,chipv arrays to store the polarisation of the parent beams, including the initial user-defined value and makes plasma_in always use these to compute the coupling. In addition, in the case the polarisation is not immediately known (i.e. if raytracing.ipol == .false.), this change postpones the computation of the Jones vector (ext, eyt) from the launch point, if the magnetic equilibrium is available, to when the ray actually crosses the plasma boundary. The original code, besides being strictly incorrect, can lead to non-negligible alterations to the coupling. This change also mean: 1. most of the functionality of `set_pol` has been merged with `plasma_in` 2. the polarisation is undefined and the Jones vector is set to the placeholder value [1, 0] till `plasma_im` is called Finally, `polarcold` is removed because it's unused.
2024-03-02 17:32:03 +01:00
call plasma_in(jk, xv, anv, bres, sox, cpl, psipol, chipol, iop, ext, eyt, &
perfect=.not. params%raytracing%ipol &
.and. params%antenna%iox == iox &
.and. iop(jk) == 0 .and. ip==1)
if(iop(jk) == 1 .and. ip==1) then ! * 1st entrance on 1st pass (ray hasn't entered in plasma yet) => continue current pass
fix coupling for subsequent beams In situations when multiple beams are traced, either when allowing multiple plasma crossings (raytracing.ipass > 0) or the initial polarisation is mixed (raytracing.ipol == .true.), the couplings of all but the first beam (with least index_rt) were invalid. The bug is due to the re-use of the psipol,chipol variables as the beams are traced sequentially over the beam_loop. For the first beam being traced the psipol,chipol are correctly initialised to the user-defined value and the resulting coupling is correct. However, in each subsequent beam the values were not set to those of the parent beam (or to the user-defined value in the case of the first X mode beam), but to those of the previous beams (current index_rt - 1). This change repurposes the psipv,chipv arrays to store the polarisation of the parent beams, including the initial user-defined value and makes plasma_in always use these to compute the coupling. In addition, in the case the polarisation is not immediately known (i.e. if raytracing.ipol == .false.), this change postpones the computation of the Jones vector (ext, eyt) from the launch point, if the magnetic equilibrium is available, to when the ray actually crosses the plasma boundary. The original code, besides being strictly incorrect, can lead to non-negligible alterations to the coupling. This change also mean: 1. most of the functionality of `set_pol` has been merged with `plasma_in` 2. the polarisation is undefined and the Jones vector is set to the placeholder value [1, 0] till `plasma_im` is called Finally, `polarcold` is removed because it's unused.
2024-03-02 17:32:03 +01:00
if(cpl(iox) < etaucr) then ! + 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 == 1) then
2021-12-18 18:57:38 +01:00
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')
fix coupling for subsequent beams In situations when multiple beams are traced, either when allowing multiple plasma crossings (raytracing.ipass > 0) or the initial polarisation is mixed (raytracing.ipol == .true.), the couplings of all but the first beam (with least index_rt) were invalid. The bug is due to the re-use of the psipol,chipol variables as the beams are traced sequentially over the beam_loop. For the first beam being traced the psipol,chipol are correctly initialised to the user-defined value and the resulting coupling is correct. However, in each subsequent beam the values were not set to those of the parent beam (or to the user-defined value in the case of the first X mode beam), but to those of the previous beams (current index_rt - 1). This change repurposes the psipv,chipv arrays to store the polarisation of the parent beams, including the initial user-defined value and makes plasma_in always use these to compute the coupling. In addition, in the case the polarisation is not immediately known (i.e. if raytracing.ipol == .false.), this change postpones the computation of the Jones vector (ext, eyt) from the launch point, if the magnetic equilibrium is available, to when the ray actually crosses the plasma boundary. The original code, besides being strictly incorrect, can lead to non-negligible alterations to the coupling. This change also mean: 1. most of the functionality of `set_pol` has been merged with `plasma_in` 2. the polarisation is undefined and the Jones vector is set to the placeholder value [1, 0] till `plasma_im` is called Finally, `polarcold` is removed because it's unused.
2024-03-02 17:32:03 +01:00
write (msg,'(" polarisation: ψ=",g0.5,"°, χ=",g0.5,"°")') &
psipol, chipol
call log_debug(msg, mod='gray_core', proc='gray_main')
psipv(index_rt) = psipol ! + polarization angles at plasma boundary for central ray
chipv(index_rt) = chipol
end if
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 < params%raytracing%ipass) then ! + not last pass
2019-03-28 10:50:28 +01:00
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) - params%raytracing%dst ! . starting step for next pass = last step
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) < 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
fix coupling for subsequent beams In situations when multiple beams are traced, either when allowing multiple plasma crossings (raytracing.ipass > 0) or the initial polarisation is mixed (raytracing.ipol == .true.), the couplings of all but the first beam (with least index_rt) were invalid. The bug is due to the re-use of the psipol,chipol variables as the beams are traced sequentially over the beam_loop. For the first beam being traced the psipol,chipol are correctly initialised to the user-defined value and the resulting coupling is correct. However, in each subsequent beam the values were not set to those of the parent beam (or to the user-defined value in the case of the first X mode beam), but to those of the previous beams (current index_rt - 1). This change repurposes the psipv,chipv arrays to store the polarisation of the parent beams, including the initial user-defined value and makes plasma_in always use these to compute the coupling. In addition, in the case the polarisation is not immediately known (i.e. if raytracing.ipol == .false.), this change postpones the computation of the Jones vector (ext, eyt) from the launch point, if the magnetic equilibrium is available, to when the ray actually crosses the plasma boundary. The original code, besides being strictly incorrect, can lead to non-negligible alterations to the coupling. This change also mean: 1. most of the functionality of `set_pol` has been merged with `plasma_in` 2. the polarisation is undefined and the Jones vector is set to the placeholder value [1, 0] till `plasma_im` is called Finally, `polarcold` is removed because it's unused.
2024-03-02 17:32:03 +01:00
taus(jk,child_index_rt:child_index_rt+1) = tau1(jk) + tau0(jk) ! . starting tau for next O-mode pass
cpls(jk,child_index_rt) = cpl1(jk) * cpl(1) ! . cumulative coupling for next O-mode pass
cpls(jk,child_index_rt+1) = cpl1(jk) * cpl(2) ! . cumulative coupling for next X-mode pass
if(jk == 1) then ! . polarization angles at plasma boundary for central ray
fix coupling for subsequent beams In situations when multiple beams are traced, either when allowing multiple plasma crossings (raytracing.ipass > 0) or the initial polarisation is mixed (raytracing.ipol == .true.), the couplings of all but the first beam (with least index_rt) were invalid. The bug is due to the re-use of the psipol,chipol variables as the beams are traced sequentially over the beam_loop. For the first beam being traced the psipol,chipol are correctly initialised to the user-defined value and the resulting coupling is correct. However, in each subsequent beam the values were not set to those of the parent beam (or to the user-defined value in the case of the first X mode beam), but to those of the previous beams (current index_rt - 1). This change repurposes the psipv,chipv arrays to store the polarisation of the parent beams, including the initial user-defined value and makes plasma_in always use these to compute the coupling. In addition, in the case the polarisation is not immediately known (i.e. if raytracing.ipol == .false.), this change postpones the computation of the Jones vector (ext, eyt) from the launch point, if the magnetic equilibrium is available, to when the ray actually crosses the plasma boundary. The original code, besides being strictly incorrect, can lead to non-negligible alterations to the coupling. This change also mean: 1. most of the functionality of `set_pol` has been merged with `plasma_in` 2. the polarisation is undefined and the Jones vector is set to the placeholder value [1, 0] till `plasma_im` is called Finally, `polarcold` is removed because it's unused.
2024-03-02 17:32:03 +01:00
psipv(child_index_rt:child_index_rt+1) = psipol
chipv(child_index_rt:child_index_rt+1) = chipol
end if
2019-03-28 10:50:28 +01:00
else ! * 1st entrance on 2nd+ pass (ray hasn't entered in plasma since end of previous pass) => continue current pass
cpl = [zero, zero]
end if
end if
else if(ext_pl) then ! ray exits plasma
2022-05-26 02:33:55 +02:00
write (msg, '(" ray ", g0, " left plasma")') jk
call log_debug(msg, mod='gray_core', proc='gray_main')
call plasma_out(jk,xv,anv,bres,sox,iop,ext,eyt)
end if
if(ent_wl) then ! ray enters vessel
2019-03-28 10:50:28 +01:00
iow(jk)=iow(jk)+1 ! * out->in
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
error = ior(error,ierrn)
call print_err_raytracing(ierrn, i, anpl)
end if
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
2019-03-28 10:50:28 +01:00
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 < 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
fix coupling for subsequent beams In situations when multiple beams are traced, either when allowing multiple plasma crossings (raytracing.ipass > 0) or the initial polarisation is mixed (raytracing.ipol == .true.), the couplings of all but the first beam (with least index_rt) were invalid. The bug is due to the re-use of the psipol,chipol variables as the beams are traced sequentially over the beam_loop. For the first beam being traced the psipol,chipol are correctly initialised to the user-defined value and the resulting coupling is correct. However, in each subsequent beam the values were not set to those of the parent beam (or to the user-defined value in the case of the first X mode beam), but to those of the previous beams (current index_rt - 1). This change repurposes the psipv,chipv arrays to store the polarisation of the parent beams, including the initial user-defined value and makes plasma_in always use these to compute the coupling. In addition, in the case the polarisation is not immediately known (i.e. if raytracing.ipol == .false.), this change postpones the computation of the Jones vector (ext, eyt) from the launch point, if the magnetic equilibrium is available, to when the ray actually crosses the plasma boundary. The original code, besides being strictly incorrect, can lead to non-negligible alterations to the coupling. This change also mean: 1. most of the functionality of `set_pol` has been merged with `plasma_in` 2. the polarisation is undefined and the Jones vector is set to the placeholder value [1, 0] till `plasma_im` is called Finally, `polarcold` is removed because it's unused.
2024-03-02 17:32:03 +01:00
call plasma_in(jk, xv, anv, bres, sox, cpl, psipol, chipol, & ! . ray re-enters plasma after reflection
iop, ext, eyt, perfect=.false.)
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) < 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
fix coupling for subsequent beams In situations when multiple beams are traced, either when allowing multiple plasma crossings (raytracing.ipass > 0) or the initial polarisation is mixed (raytracing.ipol == .true.), the couplings of all but the first beam (with least index_rt) were invalid. The bug is due to the re-use of the psipol,chipol variables as the beams are traced sequentially over the beam_loop. For the first beam being traced the psipol,chipol are correctly initialised to the user-defined value and the resulting coupling is correct. However, in each subsequent beam the values were not set to those of the parent beam (or to the user-defined value in the case of the first X mode beam), but to those of the previous beams (current index_rt - 1). This change repurposes the psipv,chipv arrays to store the polarisation of the parent beams, including the initial user-defined value and makes plasma_in always use these to compute the coupling. In addition, in the case the polarisation is not immediately known (i.e. if raytracing.ipol == .false.), this change postpones the computation of the Jones vector (ext, eyt) from the launch point, if the magnetic equilibrium is available, to when the ray actually crosses the plasma boundary. The original code, besides being strictly incorrect, can lead to non-negligible alterations to the coupling. This change also mean: 1. most of the functionality of `set_pol` has been merged with `plasma_in` 2. the polarisation is undefined and the Jones vector is set to the placeholder value [1, 0] till `plasma_im` is called Finally, `polarcold` is removed because it's unused.
2024-03-02 17:32:03 +01:00
taus(jk,child_index_rt:child_index_rt+1) = tau1(jk) + tau0(jk) ! . starting tau for next O-mode pass
cpls(jk,child_index_rt) = cpl1(jk) * cpl(1) ! . cumulative coupling for next O-mode pass
cpls(jk,child_index_rt+1) = cpl1(jk) * cpl(2) ! . cumulative coupling for next X-mode pass
if(jk == 1) then ! + polarization angles at plasma boundary for central ray
fix coupling for subsequent beams In situations when multiple beams are traced, either when allowing multiple plasma crossings (raytracing.ipass > 0) or the initial polarisation is mixed (raytracing.ipol == .true.), the couplings of all but the first beam (with least index_rt) were invalid. The bug is due to the re-use of the psipol,chipol variables as the beams are traced sequentially over the beam_loop. For the first beam being traced the psipol,chipol are correctly initialised to the user-defined value and the resulting coupling is correct. However, in each subsequent beam the values were not set to those of the parent beam (or to the user-defined value in the case of the first X mode beam), but to those of the previous beams (current index_rt - 1). This change repurposes the psipv,chipv arrays to store the polarisation of the parent beams, including the initial user-defined value and makes plasma_in always use these to compute the coupling. In addition, in the case the polarisation is not immediately known (i.e. if raytracing.ipol == .false.), this change postpones the computation of the Jones vector (ext, eyt) from the launch point, if the magnetic equilibrium is available, to when the ray actually crosses the plasma boundary. The original code, besides being strictly incorrect, can lead to non-negligible alterations to the coupling. This change also mean: 1. most of the functionality of `set_pol` has been merged with `plasma_in` 2. the polarisation is undefined and the Jones vector is set to the placeholder value [1, 0] till `plasma_im` is called Finally, `polarcold` is removed because it's unused.
2024-03-02 17:32:03 +01:00
psipv(child_index_rt:child_index_rt+1) = psipol
chipv(child_index_rt:child_index_rt+1) = chipol
end if
end if
end if
end if
iopmin = min(iopmin,iop(jk))
if(ip < params%raytracing%ipass) then ! not last pass
2019-03-28 10:50:28 +01:00
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)
! 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)
block
complex(wp_) :: Npr_warm
call alpha_effj(params%ecrh_cd, psinv, xg, yg, dens, tekev, &
ak0, bres, derdnm, anpl, anpr, sox, Npr_warm, &
alpha, didp, nharm, nhf, iokhawa, ierrhcd)
anprre = Npr_warm%re
anprim = Npr_warm%im
if(ierrhcd /= 0) then
error = ior(error, ierrhcd)
call print_err_ecrh_cd(ierrhcd, i, Npr_warm, alpha)
end if
end block
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 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
2023-03-29 11:19:00 +02:00
dids = didp * pow * alpha ! current driven: dI/ds = dI/dP⋅dP/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)
2019-03-28 10:50:28 +01:00
if(iwait(jk)) then ! copy values from last pass for inactive ray
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, &
2019-03-28 10:50:28 +01:00
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 rays_loop
! ============ rays loop END ============
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,params%output%istpr) == 0) then
if(params%raytracing%nray > 1 .and. all(.not.iwait)) &
call print_projxyzt(stv,yw,0)
end if
! test whether further trajectory integration is unnecessary
call vmaxmin(tau1+tau0+lgcpl1,params%raytracing%nray,taumn,taumx) ! test on tau + coupling
if(is_critical(error)) then ! stop propagation for current beam
istop_pass = istop_pass +1 ! * +1 non propagating beam
if(ip < params%raytracing%ipass) call turnoffray(0,ip,ib,iroff) ! * de-activate derived beams
exit
2019-03-28 10:50:28 +01:00
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 propagation_loop
2021-12-18 18:57:38 +01:00
call log_debug(' propagation loop end', mod='gray_core', proc='gray_main')
! ======== propagation loop END ========
! print all ray positions in local reference system
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>params%raytracing%nstep) i=params%raytracing%nstep
pabs_beam = sum(ppabs(:,i))
icd_beam = sum(ccci(:,i))
call vmaxmin(tau0,params%raytracing%nray,taumn,taumx) ! taumn,taumx for print
2024-01-30 10:14:24 +01:00
if (params%equilibrium%iequil /= EQ_VACUUM) then
! 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)
end if
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 < params%raytracing%ipass .and. iopmin > 2) then ! not last pass AND at least one ray re-entered plasma
fix coupling for subsequent beams In situations when multiple beams are traced, either when allowing multiple plasma crossings (raytracing.ipass > 0) or the initial polarisation is mixed (raytracing.ipol == .true.), the couplings of all but the first beam (with least index_rt) were invalid. The bug is due to the re-use of the psipol,chipol variables as the beams are traced sequentially over the beam_loop. For the first beam being traced the psipol,chipol are correctly initialised to the user-defined value and the resulting coupling is correct. However, in each subsequent beam the values were not set to those of the parent beam (or to the user-defined value in the case of the first X mode beam), but to those of the previous beams (current index_rt - 1). This change repurposes the psipv,chipv arrays to store the polarisation of the parent beams, including the initial user-defined value and makes plasma_in always use these to compute the coupling. In addition, in the case the polarisation is not immediately known (i.e. if raytracing.ipol == .false.), this change postpones the computation of the Jones vector (ext, eyt) from the launch point, if the magnetic equilibrium is available, to when the ray actually crosses the plasma boundary. The original code, besides being strictly incorrect, can lead to non-negligible alterations to the coupling. This change also mean: 1. most of the functionality of `set_pol` has been merged with `plasma_in` 2. the polarisation is undefined and the Jones vector is set to the placeholder value [1, 0] till `plasma_im` is called Finally, `polarcold` is removed because it's unused.
2024-03-02 17:32:03 +01:00
cpl_beam1 = sum(&
p0ray * exp(-tau0) * cpls(:,child_index_rt)/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) > 2) then ! * central ray O/X-mode coupling for next beam
fix coupling for subsequent beams In situations when multiple beams are traced, either when allowing multiple plasma crossings (raytracing.ipass > 0) or the initial polarisation is mixed (raytracing.ipol == .true.), the couplings of all but the first beam (with least index_rt) were invalid. The bug is due to the re-use of the psipol,chipol variables as the beams are traced sequentially over the beam_loop. For the first beam being traced the psipol,chipol are correctly initialised to the user-defined value and the resulting coupling is correct. However, in each subsequent beam the values were not set to those of the parent beam (or to the user-defined value in the case of the first X mode beam), but to those of the previous beams (current index_rt - 1). This change repurposes the psipv,chipv arrays to store the polarisation of the parent beams, including the initial user-defined value and makes plasma_in always use these to compute the coupling. In addition, in the case the polarisation is not immediately known (i.e. if raytracing.ipol == .false.), this change postpones the computation of the Jones vector (ext, eyt) from the launch point, if the magnetic equilibrium is available, to when the ray actually crosses the plasma boundary. The original code, besides being strictly incorrect, can lead to non-negligible alterations to the coupling. This change also mean: 1. most of the functionality of `set_pol` has been merged with `plasma_in` 2. the polarisation is undefined and the Jones vector is set to the placeholder value [1, 0] till `plasma_im` is called Finally, `polarcold` is removed because it's unused.
2024-03-02 17:32:03 +01:00
cpl_cbeam1 = cpls(1,child_index_rt)/cpl1(1)
cpl_cbeam2 = one - cpl_cbeam1
end if
else ! last pass OR no ray re-entered plasma
cpl_beam1 = zero
cpl_beam2 = zero
2015-11-18 17:34:33 +01:00
end if
2021-12-18 18:57:38 +01:00
! print final results for pass on screen
call log_info(' partial results:', mod='gray_core', proc='gray_main')
2022-05-26 02:33:55 +02:00
write(msg, '(" final step:", x,a,g0, x,a,g0.4)') 'n=', i, '(s, ct, Sr)=', stv(1)
2021-12-18 18:57:38 +01:00
call log_info(msg, mod='gray_core', proc='gray_main')
write(msg, '(3x,a,2(x,a,"=",g0.4))') 'optical depth:', 'τ_min', taumn, 'τ_max', taumx
call log_info(msg, mod='gray_core', proc='gray_main')
write(msg, '(3x,a,g0.3," MW")') 'absoption: P=', pabs_beam
call log_info(msg, mod='gray_core', proc='gray_main')
2023-10-04 16:56:05 +02:00
write(msg, '(3x,a,g0.3," kA")') 'current drive: I=', icd_beam*1.0e3_wp_
2021-12-18 18:57:38 +01:00
call log_info(msg, mod='gray_core', proc='gray_main')
if(ip < params%raytracing%ipass) then
2021-12-18 18:57:38 +01:00
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')
if(iop(1) > 2) then
write(msg, '(3x,a,(g0.4,", ",g0.4))') &
'coupling [ctr ray, O/X]:', cpl_cbeam1, cpl_cbeam2 ! central ray coupling for next O/X beams
fix coupling for subsequent beams In situations when multiple beams are traced, either when allowing multiple plasma crossings (raytracing.ipass > 0) or the initial polarisation is mixed (raytracing.ipol == .true.), the couplings of all but the first beam (with least index_rt) were invalid. The bug is due to the re-use of the psipol,chipol variables as the beams are traced sequentially over the beam_loop. For the first beam being traced the psipol,chipol are correctly initialised to the user-defined value and the resulting coupling is correct. However, in each subsequent beam the values were not set to those of the parent beam (or to the user-defined value in the case of the first X mode beam), but to those of the previous beams (current index_rt - 1). This change repurposes the psipv,chipv arrays to store the polarisation of the parent beams, including the initial user-defined value and makes plasma_in always use these to compute the coupling. In addition, in the case the polarisation is not immediately known (i.e. if raytracing.ipol == .false.), this change postpones the computation of the Jones vector (ext, eyt) from the launch point, if the magnetic equilibrium is available, to when the ray actually crosses the plasma boundary. The original code, besides being strictly incorrect, can lead to non-negligible alterations to the coupling. This change also mean: 1. most of the functionality of `set_pol` has been merged with `plasma_in` 2. the polarisation is undefined and the Jones vector is set to the placeholder value [1, 0] till `plasma_im` is called Finally, `polarcold` is removed because it's unused.
2024-03-02 17:32:03 +01:00
call log_info(msg, mod='gray_core', proc='gray_main')
2021-12-18 18:57:38 +01:00
end if
end if
2021-12-18 18:57:38 +01:00
2024-01-30 10:14:24 +01:00
if (params%equilibrium%iequil /= EQ_VACUUM) then
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
2024-01-30 10:14:24 +01:00
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
end if
! ============ post-proc END ============
end do beam_loop
2021-12-18 18:57:38 +01:00
call log_debug('beam loop end', mod='gray_core', proc='gray_main')
! ============ beam loop END ============
! ======= cumulative prints BEGIN =======
results%pabs = results%pabs + sum(pabs_pass) ! *final results (O+X) [gray_main output]
results%icd = results%icd + sum(icd_pass)
2021-12-18 18:57:38 +01:00
! print final results for pass on screen
call log_info(' comulative results:', mod='gray_core', proc='gray_main')
write(msg, '(" absoption [O,X mode] P=",g0.4,", ",g0.4," MW")') &
pabs_pass(1), pabs_pass(2)
call log_info(msg, mod='gray_core', proc='gray_main')
write(msg, '(" current drive [O,X mode] I=",g0.4,", ",g0.4," kA")') &
icd_pass(1)*1.0e3_wp_, icd_pass(2)*1.0e3_wp_
call log_info(msg, mod='gray_core', proc='gray_main')
! ======== cumulative prints END ========
if(istop_pass == nbeam_pass) exit ! no active beams
end do main_loop
2021-12-18 18:57:38 +01:00
call log_debug('pass loop end', mod='gray_core', proc='gray_main')
! ============ main loop END ============
! ========== free memory BEGIN ==========
call dealloc_surfvec
call dealloc_pec
! =========== free memory END ===========
end block
end associate
end subroutine gray_main
2015-11-18 17:34:33 +01:00
subroutine vectinit(psjki,ppabs,ccci,tau0,alphaabs0,dids0,ccci0,iiv)
use const_and_precisions, only : zero
2015-11-18 17:34:33 +01:00
! arguments
real(wp_), dimension(:,:), intent(out) :: psjki,ppabs,ccci
real(wp_), dimension(:), intent(out) :: tau0,alphaabs0,dids0,ccci0
2015-11-18 17:34:33 +01:00
integer, dimension(:), intent(out) :: iiv
!! common/external functions/variables
! integer :: jclosest
! real(wp_), dimension(3) :: anwcl,xwcl
!
! common/refln/anwcl,xwcl,jclosest
!
! jclosest=nrayr+1
! anwcl(1:3)=0.0_wp_
! xwcl(1:3)=0.0_wp_
psjki = zero
ppabs = zero
ccci = zero
tau0 = zero
alphaabs0 = zero
dids0 = zero
ccci0 = zero
iiv = 1
2015-11-18 17:34:33 +01:00
end subroutine vectinit
subroutine ic_gb(params, anv0c, ak0, ywrk0, ypwrk0, &
stv, 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 gray_params, only : gray_parameters
use beamdata, only : nray,nrayr,nrayth,rwmax
! subroutine arguments
type(gray_parameters), intent(in) :: params
real(wp_), dimension(3), intent(in) :: anv0c
real(wp_), intent(in) :: ak0
real(wp_), dimension(6, nray), intent(out) :: ywrk0, ypwrk0
real(wp_), dimension(nrayth * nrayr), intent(out) :: stv
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
integer, intent(in) :: index_rt
! local variables
real(wp_), dimension(3) :: xv0c
real(wp_) :: wcsi,weta,rcicsi,rcieta,phiw,phir
2015-11-18 17:34:33 +01:00
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
real(wp_) :: rcixx,rciyy,rcixy,dwwxx,dwwyy,dwwxy,d2wwxx,d2wwyy,d2wwxy
real(wp_) :: drcixx,drciyy,drcixy,dr,da,ddfu,dcsiw,detaw,dx0t,dy0t
real(wp_) :: x0t,y0t,z0t,dx0,dy0,dz0,x0,y0,z0,gxt,gyt,gzt,gr2
real(wp_) :: gxxt,gyyt,gzzt,gxyt,gxzt,gyzt,dgr2xt,dgr2yt,dgr2zt
real(wp_) :: dgr2x,dgr2y,dgr2z,pppx,pppy,denpp,ppx,ppy
real(wp_) :: anzt,anxt,anyt,anx,any,anz,an20,an0
real(wp_) :: du1tx,du1ty,du1tz,denom,ddr,ddi
real(wp_), dimension(nrayr) :: uj
real(wp_), dimension(nrayth) :: sna,csa
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)
2015-11-18 17:34:33 +01:00
csth=anv0c(3)
snth=sqrt(one-csth**2)
if(snth > zero) then
csps=anv0c(2)/snth
snps=anv0c(1)/snth
else
csps=one
snps=zero
end if
! Gaussian beam: exp[-ik0 zt] exp[-i k0/2 S(xt,yt,zt)]
! xt,yt,zt, cartesian coordinate system with zt along the beamline and xt in the z = 0 plane
! S(xt,yt,zt) = S_real +i S_imag = Qxx(zt) xt^2 + Qyy(zt) yt^2 + 2 Qxy(zt) xt yt
! (csiw, etaw) and (csiR, etaR) intensity and phase ellipse, rotated by angle phiw and phiR
! S(xt,yt,zt) = csiR^2 / Rccsi +etaR^2 /Rceta - i (csiw^2 Wcsi +etaw^2 Weta)
! Rccsi,eta curvature radius at the launching point
! Wcsi,eta =2/(k0 wcsi,eta^2) with wcsi,eta^2 beam size at the launching point
2015-11-18 17:34:33 +01:00
phiwrad = phiw*degree
phirrad = phir*degree
csphiw = cos(phiwrad)
snphiw = sin(phiwrad)
!csphir = cos(phirrad)
!snphir = sin(phirrad)
2015-11-18 17:34:33 +01:00
wwcsi = two/(ak0*wcsi**2)
wweta = two/(ak0*weta**2)
if(phir/=phiw) then
sk = rcicsi + rcieta
sw = wwcsi + wweta
dk = rcicsi - rcieta
dw = wwcsi - wweta
ts = -(dk*sin(2*phirrad) - ui*dw*sin(2*phiwrad))
tc = (dk*cos(2*phirrad) - ui*dw*cos(2*phiwrad))
phic = half*atan(ts/tc)
ddd = dk*cos(2*(phirrad+phic)) - ui*dw*cos(2*(phiwrad+phic))
sss = sk - ui*sw
qi1 = half*(sss + ddd)
qi2 = half*(sss - ddd)
rci1 = real(qi1)
rci2 = real(qi2)
ww1 = -imag(qi1)
ww2 = -imag(qi2)
2015-11-18 17:34:33 +01:00
else
rci1 = rcicsi
rci2 = rcieta
ww1 = wwcsi
ww2 = wweta
phic = -phiwrad
qi1 = rci1 - ui*ww1
qi2 = rci2 - ui*ww2
end if
!w01=sqrt(2.0_wp_/(ak0*ww1))
!d01=-rci1/(rci1**2+ww1**2)
!w02=sqrt(2.0_wp_/(ak0*ww2))
!d02=-rci2/(rci2**2+ww2**2)
2015-11-18 17:34:33 +01:00
qqxx = qi1*cos(phic)**2 + qi2*sin(phic)**2
qqyy = qi1*sin(phic)**2 + qi2*cos(phic)**2
qqxy = -(qi1 - qi2)*sin(phic)*cos(phic)
wwxx = -imag(qqxx)
wwyy = -imag(qqyy)
wwxy = -imag(qqxy)
rcixx = real(qqxx)
rciyy = real(qqyy)
rcixy = real(qqxy)
2015-11-18 17:34:33 +01:00
dqi1 = -qi1**2
dqi2 = -qi2**2
d2qi1 = 2*qi1**3
d2qi2 = 2*qi2**3
dqqxx = dqi1*cos(phic)**2 + dqi2*sin(phic)**2
dqqyy = dqi1*sin(phic)**2 + dqi2*cos(phic)**2
dqqxy = -(dqi1 - dqi2)*sin(phic)*cos(phic)
2015-11-18 17:34:33 +01:00
d2qqxx = d2qi1*cos(phic)**2 + d2qi2*sin(phic)**2
d2qqyy = d2qi1*sin(phic)**2 + d2qi2*cos(phic)**2
d2qqxy = -(d2qi1 - d2qi2)*sin(phic)*cos(phic)
dwwxx = -imag(dqqxx)
dwwyy = -imag(dqqyy)
dwwxy = -imag(dqqxy)
d2wwxx = -imag(d2qqxx)
d2wwyy = -imag(d2qqyy)
d2wwxy = -imag(d2qqxy)
drcixx = real(dqqxx)
drciyy = real(dqqyy)
drcixy = real(dqqxy)
2015-11-18 17:34:33 +01:00
if(nrayr > 1) then
dr = rwmax/(nrayr-1)
2015-11-18 17:34:33 +01:00
else
dr = one
end if
ddfu = two*dr**2/ak0 ! twodr2 = 2*dr**2 = 2*rwmax/real(nrayr-1)
2015-11-18 17:34:33 +01:00
do j = 1, nrayr
uj(j) = real(j-1, wp_)
end do
2015-11-18 17:34:33 +01:00
da=2*pi/nrayth
2015-11-18 17:34:33 +01:00
do k=1,nrayth
alfak = (k-1)*da
sna(k) = sin(alfak)
csa(k) = cos(alfak)
end do
! central ray
2015-11-18 17:34:33 +01:00
jk=1
gri(:,1) = zero
ggri(:,:,1) = zero
ywrk0(1:3,1) = xv0c
ywrk0(4:6,1) = anv0c
ypwrk0(1:3,1) = anv0c
ypwrk0(4:6,1) = zero
2015-11-18 17:34:33 +01:00
do k=1,nrayth
dcsiw = dr*csa(k)*wcsi
detaw = dr*sna(k)*weta
dx0t = dcsiw*csphiw - detaw*snphiw
dy0t = dcsiw*snphiw + detaw*csphiw
du1tx = (dx0t*wwxx + dy0t*wwxy)/ddfu
du1ty = (dx0t*wwxy + dy0t*wwyy)/ddfu
xc0(:,k,1) = xv0c
du10(1,k,1) = du1tx*csps + snps*du1ty*csth
du10(2,k,1) = -du1tx*snps + csps*du1ty*csth
du10(3,k,1) = -du1ty*snth
end do
ddr = zero
ddi = zero
! loop on rays jk>1
2015-11-18 17:34:33 +01:00
j=2
k=0
do jk=2,nray
k=k+1
if(k > nrayth) then
j=j+1
k=1
end if
!csiw=u*dcsiw
!etaw=u*detaw
!csir=x0t*csphir+y0t*snphir
!etar=-x0t*snphir+y0t*csphir
2015-11-18 17:34:33 +01:00
dcsiw = dr*csa(k)*wcsi
detaw = dr*sna(k)*weta
dx0t = dcsiw*csphiw - detaw*snphiw
dy0t = dcsiw*snphiw + detaw*csphiw
x0t = uj(j)*dx0t
y0t = uj(j)*dy0t
z0t = 0
2015-11-18 17:34:33 +01:00
dx0 = x0t*csps + snps*(y0t*csth + z0t*snth)
dy0 = -x0t*snps + csps*(y0t*csth + z0t*snth)
dz0 = z0t*csth - y0t*snth
x0 = xv0c(1) + dx0
y0 = xv0c(2) + dy0
z0 = xv0c(3) + dz0
gxt = x0t*wwxx + y0t*wwxy
gyt = x0t*wwxy + y0t*wwyy
gzt = half*(x0t**2*dwwxx + y0t**2*dwwyy ) + x0t*y0t*dwwxy
gr2 = gxt*gxt + gyt*gyt + gzt*gzt
gxxt = wwxx
gyyt = wwyy
gzzt = half*(x0t**2*d2wwxx + y0t**2*d2wwyy) + x0t*y0t*d2wwxy
gxyt = wwxy
gxzt = x0t*dwwxx + y0t*dwwxy
gyzt = x0t*dwwxy + y0t*dwwyy
dgr2xt = 2*(gxt*gxxt + gyt*gxyt + gzt*gxzt)
dgr2yt = 2*(gxt*gxyt + gyt*gyyt + gzt*gyzt)
dgr2zt = 2*(gxt*gxzt + gyt*gyzt + gzt*gzzt)
dgr2x = dgr2xt*csps + snps*(dgr2yt*csth + dgr2zt*snth)
dgr2y = -dgr2xt*snps + csps*(dgr2yt*csth + dgr2zt*snth)
dgr2z = dgr2zt*csth - dgr2yt*snth
2015-11-18 17:34:33 +01:00
gri(1,jk) = gxt*csps + snps*(gyt*csth + gzt*snth)
gri(2,jk) = -gxt*snps + csps*(gyt*csth + gzt*snth)
gri(3,jk) = gzt*csth - gyt*snth
ggri(1,1,jk) = gxxt*csps**2 &
+ snps**2 *(gyyt*csth**2 + gzzt*snth**2 + 2*snth*csth*gyzt) &
+2*snps*csps*(gxyt*csth + gxzt*snth)
ggri(2,1,jk) = csps*snps &
*(-gxxt+csth**2*gyyt + snth**2*gzzt + 2*csth*snth*gyzt) &
+(csps**2 - snps**2)*(snth*gxzt + csth*gxyt)
ggri(3,1,jk) = csth*snth*snps*(gzzt - gyyt) + (csth**2 - snth**2) &
*snps*gyzt + csps*(csth*gxzt - snth*gxyt)
ggri(1,2,jk) = ggri(2,1,jk)
ggri(2,2,jk) = gxxt*snps**2 &
+ csps**2 *(gyyt*csth**2 + gzzt*snth**2 + 2*snth*csth*gyzt) &
-2*snps*csps*(gxyt*csth + gxzt*snth)
ggri(3,2,jk) = csth*snth*csps*(gzzt - gyyt) + (csth**2-snth**2) &
*csps*gyzt + snps*(snth*gxyt - csth*gxzt)
ggri(1,3,jk) = ggri(3,1,jk)
ggri(2,3,jk) = ggri(3,2,jk)
ggri(3,3,jk) = gzzt*csth**2 + gyyt*snth**2 - 2*csth*snth*gyzt
du1tx = (dx0t*wwxx + dy0t*wwxy)/ddfu
du1ty = (dx0t*wwxy + dy0t*wwyy)/ddfu
du1tz = half*uj(j)*(dx0t**2*dwwxx + dy0t**2*dwwyy + 2*dx0t*dy0t*dwwxy)/ddfu
du10(1,k,j) = du1tx*csps + snps*(du1ty*csth + du1tz*snth)
du10(2,k,j) = -du1tx*snps + csps*(du1ty*csth + du1tz*snth)
du10(3,k,j) = du1tz*csth - du1ty*snth
pppx = x0t*rcixx + y0t*rcixy
pppy = x0t*rcixy + y0t*rciyy
denpp = pppx*gxt + pppy*gyt
if (denpp/=zero) then
ppx = -pppx*gzt/denpp
ppy = -pppy*gzt/denpp
else
ppx = zero
ppy = zero
end if
anzt = sqrt((one + gr2)/(one + ppx**2 + ppy**2))
anxt = ppx*anzt
anyt = ppy*anzt
anx = anxt*csps + snps*(anyt*csth + anzt*snth)
any =-anxt*snps + csps*(anyt*csth + anzt*snth)
anz = anzt*csth - anyt*snth
an20 = one + gr2
an0 = sqrt(an20)
xc0(1,k,j) = x0
xc0(2,k,j) = y0
xc0(3,k,j) = z0
ywrk0(1,jk) = x0
ywrk0(2,jk) = y0
ywrk0(3,jk) = z0
ywrk0(4,jk) = anx
ywrk0(5,jk) = any
ywrk0(6,jk) = anz
! integration variable
select case(params%raytracing%idst)
case(0)
! optical path: s
stv(jk) = 0
denom = an0
case(1)
! "time": c⋅t
stv(jk) = 0
denom = one
case(2)
! real eikonal: S_R
stv(jk) = half*(rcixx*x0t**2 + rciyy*y0t**2) + rcixy*x0t*y0t
denom = an20
2015-11-18 17:34:33 +01:00
end select
ypwrk0(1,jk) = anx/denom
ypwrk0(2,jk) = any/denom
ypwrk0(3,jk) = anz/denom
ypwrk0(4,jk) = dgr2x/(2*denom)
ypwrk0(5,jk) = dgr2y/(2*denom)
ypwrk0(6,jk) = dgr2z/(2*denom)
ddr = anx**2 + any**2 + anz**2 - an20
ddi = 2*(anxt*gxt + anyt*gyt + anzt*gzt)
! i=0, index_rt=1, B=0, N=0, psin=-1, Xg=0, Yg=0, gradXg=0
call print_output(0,jk,stv(jk),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])
2015-11-18 17:34:33 +01:00
end do
end subroutine ic_gb
subroutine rkstep(sox,bres,xgcn,y,yp,dgr,ddgr,igrad)
2015-11-18 17:34:33 +01:00
! Runge-Kutta integrator
! use gray_params, only : igrad
use beamdata, only : h,hh,h6
real(wp_), intent(in) :: bres,xgcn
2015-11-18 17:34:33 +01:00
real(wp_), dimension(6), intent(inout) :: y
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,sox
2015-11-18 17:34:33 +01:00
real(wp_), dimension(6) :: yy,fk1,fk2,fk3,fk4
fk1 = yp
2015-11-18 17:34:33 +01:00
yy = y + fk1*hh
call rhs(sox,bres,xgcn,yy,dgr,ddgr,fk2,igrad)
2015-11-18 17:34:33 +01:00
yy = y + fk2*hh
call rhs(sox,bres,xgcn,yy,dgr,ddgr,fk3,igrad)
2015-11-18 17:34:33 +01:00
yy = y + fk3*h
call rhs(sox,bres,xgcn,yy,dgr,ddgr,fk4,igrad)
2015-11-18 17:34:33 +01:00
y = y + h6*(fk1 + 2*fk2 + 2*fk3 + fk4)
end subroutine rkstep
subroutine rhs(sox,bres,xgcn,y,dgr,ddgr,dery,igrad)
2015-11-18 17:34:33 +01:00
! Compute right-hand side terms of the ray equations (dery)
! used in R-K integrator
! arguments
real(wp_), dimension(6), intent(in) :: y
real(wp_), intent(in) :: bres,xgcn
real(wp_), dimension(3), intent(in) :: dgr
2015-11-18 17:34:33 +01:00
real(wp_), dimension(3,3), intent(in) :: ddgr
real(wp_), dimension(6), intent(out) :: dery
integer, intent(in) :: igrad,sox
2015-11-18 17:34:33 +01:00
! local variables
real(wp_) :: dens,btot,xg,yg
2015-11-18 17:34:33 +01:00
real(wp_), dimension(3) :: xv,anv,bv,derxg,deryg
real(wp_), dimension(3,3) :: derbv
xv = y(1:3)
call plas_deriv(xv,bres,xgcn,dens,btot,bv,derbv,xg,yg,derxg,deryg)
2015-11-18 17:34:33 +01:00
anv = y(4:6)
call disp_deriv(anv,sox,xg,yg,derxg,deryg,bv,derbv,dgr,ddgr,igrad,dery)
2015-11-18 17:34:33 +01:00
end subroutine rhs
2015-11-18 17:34:33 +01:00
subroutine ywppla_upd(xv,anv,dgr,ddgr,sox,bres,xgcn,dery,psinv,dens,btot, &
bv,xg,yg,derxg,anpl,anpr,ddr,ddi,dersdst,derdnm,error,igrad)
2015-11-18 17:34:33 +01:00
! Compute right-hand side terms of the ray equations (dery)
! used after full R-K step and grad(S_I) update
use gray_errors, only : raise_error, large_npl
2015-11-18 17:34:33 +01:00
! arguments
real(wp_), dimension(3), intent(in) :: xv,anv
real(wp_), dimension(3), intent(in) :: dgr
real(wp_), dimension(3,3), intent(in) :: ddgr
integer, intent(in) :: sox
real(wp_), intent(in) :: bres,xgcn
2015-11-18 17:34:33 +01:00
real(wp_), dimension(6), intent(out) :: dery
real(wp_), intent(out) :: psinv,dens,btot,xg,yg,anpl,anpr
real(wp_), intent(out) :: ddr,ddi,dersdst,derdnm
real(wp_), dimension(3), intent(out) :: bv
integer, intent(out) :: error
2017-09-12 21:37:06 +02:00
real(wp_), dimension(3), intent(out) :: derxg
integer, intent(in) :: igrad
2015-11-18 17:34:33 +01:00
! local variables
real(wp_), dimension(3) :: deryg
2015-11-18 17:34:33 +01:00
real(wp_), dimension(3,3) :: derbv
real(wp_), parameter :: anplth1 = 0.99_wp_, anplth2 = 1.05_wp_
2015-11-18 17:34:33 +01:00
call plas_deriv(xv,bres,xgcn,dens,btot,bv,derbv,xg,yg,derxg,deryg,psinv)
call disp_deriv(anv,sox,xg,yg,derxg,deryg,bv,derbv,dgr,ddgr,igrad, &
dery,anpl,anpr,ddr,ddi,dersdst,derdnm)
if (abs(anpl) > anplth2) then
error = raise_error(error, large_npl, 1)
else if (abs(anpl) > anplth1) then
error = raise_error(error, large_npl, 0)
else
error = 0
end if
2015-11-18 17:34:33 +01:00
end subroutine ywppla_upd
2015-11-18 17:34:33 +01:00
subroutine gradi_upd(ywrk,ak0,xc,du1,gri,ggri)
use const_and_precisions, only : zero,half
2015-11-18 17:34:33 +01:00
use beamdata, only : nray,nrayr,nrayth,twodr2
real(wp_), intent(in) :: ak0
real(wp_), dimension(6,nray), intent(in) :: ywrk
real(wp_), dimension(3,nrayth,nrayr), intent(inout) :: xc,du1
real(wp_), dimension(3,nray), intent(out) :: gri
real(wp_), dimension(3,3,nray), intent(out) :: ggri
! local variables
real(wp_), dimension(3,nrayth,nrayr) :: xco,du1o
integer :: jk,j,jm,jp,k,km,kp
real(wp_) :: ux,uxx,uxy,uxz,uy,uyy,uyz,uz,uzz
real(wp_) :: dfuu,dffiu,gx,gxx,gxy,gxz,gy,gyy,gyz,gz,gzz
real(wp_), dimension(3) :: dxv1,dxv2,dxv3,dgu
real(wp_), dimension(3,3) :: dgg,dff
! update position and du1 vectors
xco = xc
du1o = du1
jk = 1
do j=1,nrayr
do k=1,nrayth
if(j>1) jk=jk+1
xc(1:3,k,j)=ywrk(1:3,jk)
end do
end do
! compute grad u1 for central ray
j = 1
jp = 2
do k=1,nrayth
if(k == 1) then
km = nrayth
else
km = k-1
end if
if(k == nrayth) then
kp = 1
else
kp = k+1
end if
dxv1 = xc(:,k ,jp) - xc(:,k ,j)
dxv2 = xc(:,kp,jp) - xc(:,km,jp)
dxv3 = xc(:,k ,j) - xco(:,k ,j)
call solg0(dxv1,dxv2,dxv3,dgu)
du1(:,k,j) = dgu
end do
gri(:,1) = zero
2015-11-18 17:34:33 +01:00
! compute grad u1 and grad(S_I) for all the other rays
dfuu=twodr2/ak0 ! twodr2 = 2*dr**2 = 2*(rwmax/(nrayr-1))**2
jm=1
j=2
k=0
dffiu = dfuu
do jk=2,nray
k=k+1
if(k > nrayth) then
jm = j
j = j+1
k = 1
dffiu = dfuu*jm
end if
kp = k+1
km = k-1
if (k == 1) then
km=nrayth
else if (k == nrayth) then
kp=1
end if
dxv1 = xc(:,k ,j) - xc(:,k ,jm)
dxv2 = xc(:,kp,j) - xc(:,km,j)
dxv3 = xc(:,k ,j) - xco(:,k ,j)
call solg0(dxv1,dxv2,dxv3,dgu)
du1(:,k,j) = dgu
gri(:,jk) = dgu(:)*dffiu
end do
! compute derivatives of grad u and grad(S_I) for rays jk>1
ggri(:,:,1) = zero
jm=1
j=2
k=0
dffiu = dfuu
do jk=2,nray
k=k+1
if(k > nrayth) then
jm=j
j=j+1
k=1
dffiu = dfuu*jm
end if
kp=k+1
km=k-1
if (k == 1) then
km=nrayth
else if (k == nrayth) then
kp=1
end if
dxv1 = xc(:,k ,j) - xc(:,k ,jm)
dxv2 = xc(:,kp,j) - xc(:,km,j)
dxv3 = xc(:,k ,j) - xco(:,k ,j)
dff(:,1) = du1(:,k ,j) - du1(:,k ,jm)
dff(:,2) = du1(:,kp,j) - du1(:,km,j)
dff(:,3) = du1(:,k ,j) - du1o(:,k ,j)
call solg3(dxv1,dxv2,dxv3,dff,dgg)
! derivatives of u
ux = du1(1,k,j)
uy = du1(2,k,j)
uz = du1(3,k,j)
uxx = dgg(1,1)
uyy = dgg(2,2)
uzz = dgg(3,3)
uxy = (dgg(1,2) + dgg(2,1))*half
uxz = (dgg(1,3) + dgg(3,1))*half
uyz = (dgg(2,3) + dgg(3,2))*half
! derivatives of S_I and Grad(S_I)
gx = ux*dffiu
gy = uy*dffiu
gz = uz*dffiu
gxx = dfuu*ux*ux + dffiu*uxx
gyy = dfuu*uy*uy + dffiu*uyy
gzz = dfuu*uz*uz + dffiu*uzz
gxy = dfuu*ux*uy + dffiu*uxy
gxz = dfuu*ux*uz + dffiu*uxz
gyz = dfuu*uy*uz + dffiu*uyz
ggri(1,1,jk)=gxx
ggri(2,1,jk)=gxy
ggri(3,1,jk)=gxz
ggri(1,2,jk)=gxy
ggri(2,2,jk)=gyy
ggri(3,2,jk)=gyz
ggri(1,3,jk)=gxz
ggri(2,3,jk)=gyz
ggri(3,3,jk)=gzz
end do
2015-11-18 17:34:33 +01:00
end subroutine gradi_upd
2015-11-18 17:34:33 +01:00
subroutine solg0(dxv1,dxv2,dxv3,dgg)
! solution of the linear system of 3 eqs : dgg . dxv = dff
! input vectors : dxv1, dxv2, dxv3, dff
! output vector : dgg
! dff=(1,0,0)
! arguments
real(wp_), dimension(3), intent(in) :: dxv1,dxv2,dxv3
real(wp_), dimension(3), intent(out) :: dgg
! local variables
real(wp_) :: denom,aa1,aa2,aa3
aa1 = (dxv2(2)*dxv3(3) - dxv3(2)*dxv2(3))
aa2 = (dxv1(2)*dxv3(3) - dxv3(2)*dxv1(3))
aa3 = (dxv1(2)*dxv2(3) - dxv2(2)*dxv1(3))
denom = dxv1(1)*aa1 - dxv2(1)*aa2 + dxv3(1)*aa3
dgg(1) = aa1/denom
dgg(2) = -(dxv2(1)*dxv3(3) - dxv3(1)*dxv2(3))/denom
dgg(3) = (dxv2(1)*dxv3(2) - dxv3(1)*dxv2(2))/denom
end subroutine solg0
subroutine solg3(dxv1,dxv2,dxv3,dff,dgg)
! rhs "matrix" dff, result in dgg
! arguments
real(wp_), dimension(3), intent(in) :: dxv1,dxv2,dxv3
real(wp_), dimension(3,3), intent(in) :: dff
real(wp_), dimension(3,3), intent(out) :: dgg
! local variables
real(wp_) denom,a11,a21,a31,a12,a22,a32,a13,a23,a33
a11 = (dxv2(2)*dxv3(3) - dxv3(2)*dxv2(3))
a21 = (dxv1(2)*dxv3(3) - dxv3(2)*dxv1(3))
a31 = (dxv1(2)*dxv2(3) - dxv2(2)*dxv1(3))
a12 = (dxv2(1)*dxv3(3) - dxv3(1)*dxv2(3))
a22 = (dxv1(1)*dxv3(3) - dxv3(1)*dxv1(3))
a32 = (dxv1(1)*dxv2(3) - dxv2(1)*dxv1(3))
a13 = (dxv2(1)*dxv3(2) - dxv3(1)*dxv2(2))
a23 = (dxv1(1)*dxv3(2) - dxv3(1)*dxv1(2))
a33 = (dxv1(1)*dxv2(2) - dxv2(1)*dxv1(2))
denom = dxv1(1)*a11 - dxv2(1)*a21 + dxv3(1)*a31
dgg(:,1) = ( dff(:,1)*a11 - dff(:,2)*a21 + dff(:,3)*a31)/denom
dgg(:,2) = (-dff(:,1)*a12 + dff(:,2)*a22 - dff(:,3)*a32)/denom
dgg(:,3) = ( dff(:,1)*a13 - dff(:,2)*a23 + dff(:,3)*a33)/denom
end subroutine solg3
subroutine plas_deriv(xv, bres, xgcn, dens, btot, bv, derbv, &
xg, yg, derxg, deryg, psinv_opt)
use const_and_precisions, only : zero, cm
use gray_params, only : iequil
use equilibrium, only : psia, pol_flux, pol_curr, sgnbphi
use coreprofiles, only : density
! subroutine arguments
real(wp_), dimension(3), intent(in) :: xv
real(wp_), intent(in) :: xgcn, bres
real(wp_), intent(out) :: dens, btot, xg, yg
real(wp_), dimension(3), intent(out) :: bv, derxg, deryg
2015-11-18 17:34:33 +01:00
real(wp_), dimension(3,3), intent(out) :: derbv
real(wp_), optional, intent(out) :: psinv_opt
! local variables
2015-11-18 17:34:33 +01:00
integer :: jv
real(wp_) :: xx,yy,zz
real(wp_) :: psinv
real(wp_) :: csphi,drrdx,drrdy,dphidx,dphidy,rr,rr2,rrm,snphi,zzm
2015-11-18 17:34:33 +01:00
real(wp_), dimension(3) :: dbtot,bvc
real(wp_), dimension(3,3) :: dbvcdc,dbvdc,dbv
real(wp_) :: brr,bphi,bzz,dxgdpsi
real(wp_) :: dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,dfpv,ddensdpsi
2015-11-18 17:34:33 +01:00
xg = zero
yg = 99._wp_
psinv = -1._wp_
dens = zero
btot = zero
derxg = zero
deryg = zero
bv = zero
derbv = zero
if(iequil==0) then
! copy optional output
if (present(psinv_opt)) psinv_opt = psinv
return
end if
2015-11-18 17:34:33 +01:00
dbtot = zero
dbv = zero
dbvcdc = zero
dbvcdc = zero
dbvdc = zero
xx = xv(1)
yy = xv(2)
zz = xv(3)
! cylindrical coordinates
2015-11-18 17:34:33 +01:00
rr2 = xx**2 + yy**2
rr = sqrt(rr2)
csphi = xx/rr
snphi = yy/rr
bv(1) = -snphi*sgnbphi
bv(2) = csphi*sgnbphi
! convert from cm to meters
2015-11-18 17:34:33 +01:00
zzm = 1.0e-2_wp_*zz
rrm = 1.0e-2_wp_*rr
call pol_flux(rrm, zzm, psinv, dpsidr, dpsidz, ddpsidrr, ddpsidzz, ddpsidrz)
call pol_curr(psinv, fpolv, dfpv)
2015-11-18 17:34:33 +01:00
! copy optional output
if (present(psinv_opt)) psinv_opt = psinv
! compute yg and derivative
2015-11-18 17:34:33 +01:00
if(psinv < zero) then
bphi = fpolv/rrm
btot = abs(bphi)
yg = btot/bres
return
end if
! compute xg and derivative
call density(psinv, dens, ddensdpsi)
2015-11-18 17:34:33 +01:00
xg = xgcn*dens
dxgdpsi = xgcn*ddensdpsi
! B = f(psi)/R e_phi+ grad psi x e_phi/R
bphi = fpolv/rrm
brr = -dpsidz*psia/rrm
bzz = +dpsidr*psia/rrm
2015-11-18 17:34:33 +01:00
! bvc(i) = B_i in cylindrical coordinates
bvc = [brr, bphi, bzz]
2015-11-18 17:34:33 +01:00
! bv(i) = B_i in cartesian coordinates
2015-11-18 17:34:33 +01:00
bv(1)=bvc(1)*csphi - bvc(2)*snphi
bv(2)=bvc(1)*snphi + bvc(2)*csphi
bv(3)=bvc(3)
! dbvcdc(iv,jv) = d Bcil(iv) / dxvcil(jv)
dbvcdc(1,1) = -ddpsidrz*psia/rrm - brr/rrm
dbvcdc(1,3) = -ddpsidzz*psia/rrm
dbvcdc(2,1) = dfpv*dpsidr/rrm - bphi/rrm
dbvcdc(2,3) = dfpv*dpsidz/rrm
dbvcdc(3,1) = +ddpsidrr*psia/rrm - bzz/rrm
dbvcdc(3,3) = +ddpsidrz*psia/rrm
2015-11-18 17:34:33 +01:00
! dbvdc(iv,jv) = d Bcart(iv) / dxvcil(jv)
2015-11-18 17:34:33 +01:00
dbvdc(1,1) = dbvcdc(1,1)*csphi - dbvcdc(2,1)*snphi
dbvdc(2,1) = dbvcdc(1,1)*snphi + dbvcdc(2,1)*csphi
dbvdc(3,1) = dbvcdc(3,1)
dbvdc(1,2) = -bv(2)
dbvdc(2,2) = bv(1)
2023-10-04 16:56:05 +02:00
dbvdc(3,2) = dbvcdc(3,2) ! = 0
2015-11-18 17:34:33 +01:00
dbvdc(1,3) = dbvcdc(1,3)*csphi - dbvcdc(2,3)*snphi
dbvdc(2,3) = dbvcdc(1,3)*snphi + dbvcdc(2,3)*csphi
dbvdc(3,3) = dbvcdc(3,3)
drrdx = csphi
drrdy = snphi
dphidx = -snphi/rrm
dphidy = csphi/rrm
! dbv(iv,jv) = d Bcart(iv) / dxvcart(jv)
2015-11-18 17:34:33 +01:00
dbv(:,1) = drrdx*dbvdc(:,1) + dphidx*dbvdc(:,2)
dbv(:,2) = drrdy*dbvdc(:,1) + dphidy*dbvdc(:,2)
dbv(:,3) = dbvdc(:,3)
! B magnitude and derivatives
btot = norm2(bv)
2015-11-18 17:34:33 +01:00
2023-10-04 16:56:05 +02:00
! dbtot(i) = d |B| / dxvcart(i)
dbtot = matmul(bv, dbv)/btot
2015-11-18 17:34:33 +01:00
yg = btot/Bres
! convert spatial derivatives from dummy/m -> dummy/cm
! to be used in rhs
2015-11-18 17:34:33 +01:00
! bv(i) = B_i / B ; derbv(i,j) = d (B_i / B) /d x,y,z
2015-11-18 17:34:33 +01:00
deryg = 1.0e-2_wp_*dbtot/Bres
bv = bv/btot
do jv=1,3
derbv(:,jv) = cm * (dbv(:,jv) - bv(:)*dbtot(jv))/btot
2015-11-18 17:34:33 +01:00
end do
derxg(1) = cm * drrdx*dpsidr*dxgdpsi
derxg(2) = cm * drrdy*dpsidr*dxgdpsi
derxg(3) = cm * dpsidz *dxgdpsi
2015-11-18 17:34:33 +01:00
end subroutine plas_deriv
subroutine disp_deriv(anv, sox, xg, yg, derxg, deryg, bv, derbv, &
dgr, ddgr, igrad, dery, anpl_, anpr, ddr, ddi, &
dersdst, derdnm_)
! Computes the dispersion relation, derivatives and other
! related quantities
!
! The mandatory outputs are used for both integrating the ray
! trajectory (`rhs` subroutine) and updating the ray state
! (`ywppla_upd` suborutine); while the optional ones are used for
! computing the absoprtion and current drive.
use const_and_precisions, only : zero, one, half, two
use gray_params, only : idst
! subroutine arguments
! Inputs
! refractive index N̅ vector, b̅ = B̅/B magnetic field direction
real(wp_), dimension(3), intent(in) :: anv, bv
! sign of polarisation mode: -1 ⇒ O, +1 ⇒ X
integer, intent(in) :: sox
! CMA diagram variables: X=(ω_pe/ω)², Y=ω_ce/ω
real(wp_), intent(in) :: xg, yg
! gradients of X, Y
real(wp_), dimension(3), intent(in) :: derxg, deryg
! gradients of the complex eikonal
real(wp_), dimension(3), intent(in) :: dgr ! ∇S_I
real(wp_), dimension(3, 3), intent(in) :: ddgr ! ∇∇S_I
! gradient of the magnetic field direction
real(wp_), dimension(3, 3), intent(in) :: derbv ! ∇b̅
! raytracing/beamtracing switch
integer, intent(in) :: igrad
! Outputs
! the actual derivatives: (∂Λ/∂x̅, -∂Λ/∂N̅)
real(wp_), dimension(6), intent(out) :: dery
! additional quantities:
! refractive index
real(wp_), optional, intent(out) :: anpl_, anpr ! N∥, N⊥
! real and imaginary part of the dispersion
real(wp_), optional, intent(out) :: ddr, ddi ! Λ, ∂Λ/∂N̅⋅∇S_I
! Jacobian ds/dσ, where s arclength, σ integration variable
real(wp_), optional, intent(out) :: dersdst
! |∂Λ/∂N̅|
real(wp_), optional, intent(out) :: derdnm_
! Note: assign values to missing optional arguments results in a segfault.
! Since some are always needed anyway, we store them here and copy
! them later, if needed:
real(wp_) :: anpl, derdnm
! local variables
real(wp_) :: gr2, yg2, anpl2, del, dnl, duh, dan2sdnpl, an2, an2s
real(wp_) :: dan2sdxg, dan2sdyg, denom
real(wp_) :: derdom, dfdiadnpl, dfdiadxg, dfdiadyg, fdia, bdotgr
real(wp_), dimension(3) :: derdxv, danpldxv, derdnv, dbgr
2015-11-18 17:34:33 +01:00
an2 = dot_product(anv, anv) ! N²
anpl = dot_product(anv, bv) ! N∥ = N̅⋅B̅
2015-11-18 17:34:33 +01:00
! Shorthands used in the expressions below
yg2 = yg**2 ! Y²
anpl2 = anpl**2 ! N∥²
dnl = one - anpl2 ! 1 - N∥²
duh = one - xg - yg2 ! UH denom (duh=0 on the upper-hybrid resonance)
! Compute/copy optional outputs
if (present(anpr)) anpr = sqrt(max(an2 - anpl2, zero)) ! N⊥
if (present(anpl_)) anpl_ = anpl ! N∥
2015-11-18 17:34:33 +01:00
an2s = one
dan2sdxg = zero
dan2sdyg = zero
dan2sdnpl = zero
del = zero
fdia = zero
dfdiadnpl = zero
dfdiadxg = zero
dfdiadyg = zero
if(xg > zero) then
! Derivatives of the cold plasma refractive index
!
! N²s = 1 - X - XY²⋅(1 + N∥² ± √Δ)/[2(1 - X - Y²)]
!
! where Δ = (1 - N∥²)² + 4N∥²⋅(1 - X)/Y²
! + for the X mode, - for the O mode
! √Δ
del = sqrt(dnl**2 + 4.0_wp_*anpl2*(one - xg)/yg2)
2015-11-18 17:34:33 +01:00
! ∂(N²s)/∂X
2023-03-29 11:19:00 +02:00
! Note: this term is nonzero for X=0, but it multiplies terms
! proportional to X or ∂X/∂ψ which are zero outside the plasma.
2015-11-18 17:34:33 +01:00
dan2sdxg = - half*yg2*(one - yg2)*(one + anpl2 + sox*del)/duh**2 &
+ sox*xg*anpl2/(del*duh) - one
! ∂(N²s)/∂Y
2015-11-18 17:34:33 +01:00
dan2sdyg = - xg*yg*(one - xg)*(one + anpl2 + sox*del)/duh**2 &
+ two*sox*xg*(one - xg)*anpl2/(yg*del*duh)
! ∂(N²s)/∂N∥
2015-11-18 17:34:33 +01:00
dan2sdnpl = - xg*yg2*anpl/duh &
- sox*xg*anpl*(two*(one - xg) - yg2*dnl)/(del*duh)
if(igrad > 0) then
! Derivatives used in the complex eikonal terms (beamtracing only)
block
real(wp_) :: ddelnpl2, ddelnpl2x, ddelnpl2y, derdel
! ∂²Δ/∂N∥²
ddelnpl2 = two*(two*(one - xg)*(one + 3.0_wp_*anpl2**2) &
- yg2*dnl**3)/yg2/del**3
! ∂²(N²s)/∂N∥²
fdia = - xg*yg2*(one + half*sox*ddelnpl2)/duh
! Intermediates results used right below
derdel = two*(one - xg)*anpl2*(one + 3.0_wp_*anpl2**2) &
- dnl**2*(one + 3.0_wp_*anpl2)*yg2
derdel = 4.0_wp_*derdel/(yg*del)**5
ddelnpl2y = two*(one - xg)*derdel
ddelnpl2x = yg*derdel
! ∂³(N²s)/∂N∥³
dfdiadnpl = 24.0_wp_*sox*xg*(one - xg)*anpl*(one - anpl2**2) &
/(yg2*del**5)
! ∂³(N²s)/∂N∥²∂X
dfdiadxg = - yg2*(one - yg2)/duh**2 - sox*yg2*((one - yg2) &
*ddelnpl2 + xg*duh*ddelnpl2x)/(two*duh**2)
! ∂³(N²s)/∂N∥²∂Y
dfdiadyg = - two*yg*xg*(one - xg)/duh**2 &
- sox*xg*yg*(two*(one - xg)*ddelnpl2 &
+ yg*duh*ddelnpl2y)/(two*duh**2)
end block
2015-11-18 17:34:33 +01:00
end if
end if
! ∇(N∥) = ∇(N̅⋅b̅) = N̅⋅∇b̅
danpldxv = matmul(anv, derbv)
2015-11-18 17:34:33 +01:00
! ∂Λ/∂x̅ = - ∇(N²s) = -∂Λ/∂X ∇X - ∂Λ/∂Y ∇Y - ∂Λ/∂N∥ ∇(N∥)
derdxv = -(derxg*dan2sdxg + deryg*dan2sdyg + danpldxv*dan2sdnpl)
! ∂Λ/∂N̅ = 2N̅ - ∂(N²s)/∂N̅ = 2N̅ - ∂(N²s)/∂N∥ b̅
! Note: we used the identity ∇f(v̅⋅b̅) = f' ∇(v̅⋅b̅) = f'b̅.
derdnv = two*anv - dan2sdnpl*bv
! ∂Λ/∂ω = ∂N²/∂ω - ∂N²s/∂X⋅∂X/∂ω - ∂N²s/∂Y⋅∂Y/∂ω - ∂N²s/∂N∥⋅∂N∥/∂ω
! Notes: 1. N depends on ω: N²=c²k²/ω² ⇒ ∂N²/∂ω = -2N²/ω
! N∥=ck∥/ω ⇒ ∂N∥/∂ω = -N∥/ω
! 2. derdom is actually ω⋅∂Λ/∂ω, see below for the reason.
! 3. N gains a dependency on ω because Λ(∇S, ω) is computed
! on the constrains Λ=0.
derdom = -two*an2 + two*xg*dan2sdxg + yg*dan2sdyg + anpl*dan2sdnpl
if (igrad > 0) then
! Complex eikonal terms added to the above expressions
block
real(wp_) :: dgr2(3)
bdotgr = dot_product(bv, dgr) ! b̅⋅∇S_I
gr2 = dot_product(dgr, dgr) ! |∇S_I|²
dgr2 = 2 * matmul(ddgr, dgr) ! ∇|∇S_I|² = 2 ∇∇S_I⋅∇S_I
! ∇(b̅⋅∇S_I) = ∇b̅ᵗ ∇S_I + b̅ᵗ ∇∇S_I
! Notes: 1. We are using the convention ∇b̅_ij = ∂b_i/∂x_j.
! 2. Fortran doesn't distinguish between column/row vectors,
! so matmul(A, v̅) ≅ Av̅ and matmul(v̅, A) ≅ v̅ᵗA ≅ Aᵗv̅.
dbgr = matmul(dgr, derbv) + matmul(bv, ddgr)
! ∂Λ/∂x̅ += - ∇(|∇S_I|²) + ½ ∇[(b̅⋅∇S_I)² ∂²N²s/∂N∥²]
derdxv = derdxv - dgr2 + fdia*bdotgr*dbgr + half*bdotgr**2 &
* (derxg*dfdiadxg + deryg*dfdiadyg + danpldxv*dfdiadnpl)
! ∂Λ/∂N̅ += ½(b̅⋅∇S_I)² ∂³(N²s)/∂N∥³ b̅
! Note: we used again the identity ∇f(v̅⋅b̅) = f'b̅.
derdnv = derdnv + half*bdotgr**2*dfdiadnpl*bv
! ∂Λ/∂ω += ∂|∇S_I|²/∂ω + ½∂(b⋅∇S_I)²/∂ω + ½(b̅⋅∇S_I)² ∂/∂ω (∂²N²s/∂N∥²)
! Note: as above ∇S_I gains a dependency on ω
derdom = derdom + two*gr2 - bdotgr**2 &
* (fdia + xg*dfdiadxg + half*yg*dfdiadyg &
+ half*anpl*dfdiadnpl)
end block
end if
2015-11-18 17:34:33 +01:00
derdnm = norm2(derdnv)
! Denominator of the r.h.s, depending on the
! choice of the integration variable σ:
!
! dx̅/dσ = + ∂Λ/∂N̅ / denom
! dN̅/dσ = - ∂Λ/∂x̅ / denom
!
select case (idst)
case (0) ! σ=s, arclength
! denom = |∂Λ/∂N̅|
! Proof: Normalising ∂Λ/∂N̅ (∥ to the group velocity)
! simply makes σ the arclength parameter s.
denom = derdnm
case (1) ! σ=ct, time
! denom = -ω⋅∂Λ/∂ω
! Proof: The Hamilton equations are
!
! dx̅/dt = + ∂H/∂k̅ (H=ℏΩ, p=ℏk̅)
! dp̅/dt = - ∂H/∂x̅
!
! where Ω(x̅, k̅) is implicitely defined by the dispersion
! D(k̅, ω, x̅) = 0. By differentiating the latter we get:
!
! ∂Ω/∂x̅ = - ∂D/∂x̅ / ∂D/∂ω = - ∂Λ/∂x̅ / ∂Λ/∂ω
! ∂Ω/∂k̅ = - ∂D/∂k̅ / ∂D/∂ω = - ∂Λ/∂k̅ / ∂Λ/∂ω
!
! with Λ=D/k₀, k₀=ω/c. Finally, substituting k̅=k₀N̅:
!
! dx̅/d(ct) = + ∂Λ/∂N̅ / (-ω⋅∂Λ/∂ω)
! dN̅/d(ct) = - ∂Λ/∂x̅ / (-ω⋅∂Λ/∂ω)
!
denom = -derdom
case (2) ! s=S_R, phase
! denom = N̅⋅∂Λ/∂N̅
! Note: This parametrises the curve by the value
! of the (real) phase, with the wave frozen in time.
! Proof: By definition N̅ = k̅/k₀ = ∇S. Using the gradient
! theorem on the ray curve parametrised by the arclength
!
! S(s) = ∫₀^x̅(s) N̅⋅dl̅ = ∫₀^s N̅(σ)⋅dx̅/ds ds
!
! Differentiating gives dS = N̅⋅dx̅/ds ds, so
!
! dS = N̅⋅∂Λ/∂N̅ / |∂Λ/∂N̅| ds ⇒
!
! dx̅/dS = + ∂Λ/∂N̅ / N̅⋅∂Λ/∂N̅
! dN̅/dS = - ∂Λ/∂x̅ / N̅⋅∂Λ/∂N̅
!
denom = dot_product(anv, derdnv)
end select
2015-11-18 17:34:33 +01:00
! Jacobian ds/dσ, where s is the arclength,
! for computing integrals in dσ, like:
!
! τ = ∫α(s)ds = ∫α(σ)(ds/dσ)dσ
!
if (present(dersdst)) dersdst = derdnm/denom
! |∂Λ/∂N̅|: used in the α computation (see alpha_effj)
if (present(derdnm_)) derdnm_ = derdnm
2015-11-18 17:34:33 +01:00
! r.h.s. vector
dery(1:3) = derdnv(:)/denom ! +∂Λ/∂N̅
dery(4:6) = -derdxv(:)/denom ! -∂Λ/∂x̅
2015-11-18 17:34:33 +01:00
if (present(ddr) .or. present(ddi)) then
! Dispersion relation (geometric optics)
! ddr ← Λ = N² - N²s(X,Y,N∥) = 0
an2s = one - xg - half*xg*yg2*(one + anpl2 + sox*del)/duh
ddr = an2 - an2s
end if
if (present(ddr) .and. igrad > 0) then
! Dispersion relation (complex eikonal)
! ddr ← Λ = N² - N²s - |∇S_I|² + ½ (b̅⋅∇S_I)² ∂²N²s/∂N∥²
ddr = ddr - gr2 - half*bdotgr**2*fdia ! real part
end if
! Note: we have to return ddi even for igrad=0 because
! it's printed to udisp unconditionally
if (present(ddi)) then
! Dispersion relation (complex eikonal)
! ddi ← ∂Λ/∂N̅⋅∇S_I
ddi = merge(dot_product(derdnv, dgr), zero, igrad > 0) ! imaginary part
end if
2015-11-18 17:34:33 +01:00
end subroutine disp_deriv
subroutine alpha_effj(params, psinv, X, Y, density, temperature, &
k0, Bres, derdnm, Npl, Npr_cold, sox, &
Npr, alpha, dIdp, nhmin, nhmax, iokhawa, error)
! Computes the absorption coefficient and effective current
use const_and_precisions, only : pi, mc2=>mc2_
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 gray_errors, only : negative_absorption, raise_error
use magsurf_data, only : fluxval
2015-11-18 17:34:33 +01:00
! 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) :: X, Y
! densityity [10¹⁹ m⁻³], temperature [keV]
real(wp_), intent(in) :: density, temperature
! vacuum wavenumber k₀=ω/c, resonant B field
real(wp_), intent(in) :: k0, Bres
! group velocity: |∂Λ/∂N̅| where Λ=c²/ω²⋅D_R
real(wp_), intent(in) :: derdnm
! refractive index: N∥, N⊥ (cold dispersion)
real(wp_), intent(in) :: Npl, Npr_cold
! sign of polarisation mode: -1 ⇒ O, +1 ⇒ X
integer, intent(in) :: sox
! Outputs
! orthogonal refractive index N⊥ (solution of the warm dispersion)
complex(wp_), intent(out) :: Npr
! 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 :: nlarmor, ithn, ierrcd
real(wp_) :: mu, rbn, rbx
real(wp_) :: zeff, cst2, bmxi, bmni, fci
real(wp_), dimension(:), allocatable :: eccdpar
real(wp_) :: effjcd, effjcdav, Btot
complex(wp_) :: e(3)
alpha = 0
Npr = 0
dIdp = 0
nhmin = 0
nhmax = 0
iokhawa = 0
error = 0
! Absorption computation
! Skip when the temperature is zero (no plasma)
if (temperature <= 0) return
! Skip when the lowest resonant harmonic number is zero
mu = mc2/temperature ! μ=(mc²/kT)
call harmnumber(Y, mu, Npl**2, params%iwarm == 1, nhmin, nhmax)
if (nhmin <= 0) return
! Solve the full dispersion only when needed
if (params%iwarm /= 4 .or. params%ieccd /= 0) then
nlarmor = max(params%ilarm, nhmax)
if (params%ieccd /= 0) then
! Compute the polarisation vector only when current drive is on
call warmdisp(X, Y, mu, Npl, Npr_cold, sox, error, Npr, e, &
model=params%iwarm, nlarmor=nlarmor, &
max_iters=abs(params%imx), fallback=params%imx < 0)
else
call warmdisp(X, Y, mu, Npl, Npr_cold, sox, error, Npr, &
model=params%iwarm, nlarmor=nlarmor, &
max_iters=abs(params%imx), fallback=params%imx < 0)
end if
end if
! Compute α from the solution of the dispersion relation
! 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̅ / ∂Λ/∂ω, using the cold dispersion
! relation, we have that
!
! s̅ = ∂Λ/∂N̅ / |∂Λ/∂N̅|
! = [2N̅ - ∂(N²s)/∂N∥ b̅
! + ½(b̅⋅∇S_I)² ∂³(N²s)/∂N∥³ b̅] / |∂Λ/∂N̅|
!
! Assuming Im(k∥)=0:
!
! α = 4 Im(k⊥)⋅N⊥ / |∂Λ/∂N̅|
!
block
real(wp_) :: k_im
k_im = k0 * Npr%im ! imaginary part of k⊥
alpha = 4 * k_im*Npr_cold / derdnm
end block
if (alpha < 0) then
error = raise_error(error, negative_absorption)
return
2015-11-18 17:34:33 +01:00
end if
! Current drive computation
if (params%ieccd <= 0) return
zeff = fzeff(psinv)
ithn = 1
if (nlarmor > nhmin) ithn = 2
rhop = sqrt(psinv)
call fluxval(rhop, rri=rrii, rbav=rbavi, bmn=bmni, bmx=bmxi, fc=fci)
Btot = Y*Bres
rbn = Btot/bmni
rbx = Btot/bmxi
select case(params%ieccd)
case(1)
! Cohen model
call setcdcoeff(zeff, rbn, rbx, cst2, eccdpar)
call eccdeff(Y, Npl, Npr%re, density, mu, e, nhmin, nhmax, &
ithn, cst2, fjch, eccdpar, effjcd, iokhawa, ierrcd)
case(2)
! No trapping
call setcdcoeff(zeff, cst2, eccdpar)
call eccdeff(Y, Npl, Npr%re, density, mu, e, nhmin, nhmax, &
ithn, cst2, fjch0, eccdpar, effjcd, iokhawa, ierrcd)
case default
! Neoclassical model
call setcdcoeff(zeff, rbx, fci, mu, rhop, cst2, eccdpar)
call eccdeff(Y, Npl, Npr%re, density, mu, e, nhmin, nhmax, &
ithn, cst2, fjncl, eccdpar, effjcd, iokhawa, ierrcd)
end select
error = error + ierrcd
if (allocated(eccdpar)) deallocate(eccdpar)
! current drive efficiency <j/p> [A⋅m/W]
effjcdav = rbavi*effjcd
dIdp = sgnbphi * effjcdav / (2*pi*rrii)
2015-11-18 17:34:33 +01:00
end subroutine alpha_effj
subroutine cniteq(rqgrid,zqgrid,matr2dgrid,nr,nz,h,ncon,npts,icount,rcon,zcon)
! v2.01 12/07/95 -- written by d v bartlett, jet joint undertaking.
! (based on an older code)
! arguments
integer, intent(in) :: nr,nz
real(wp_), dimension(nr), intent(in) :: rqgrid
real(wp_), dimension(nz), intent(in) :: zqgrid
real(wp_), dimension(nr,nz), intent(in) :: matr2dgrid
real(wp_), intent(in) :: h
integer, intent(inout) :: ncon, icount
integer, dimension(ncon), intent(out) :: npts
real(wp_), dimension(icount), intent(out) :: rcon,zcon
! local variables
integer :: i,j,k,l,nrqmax,iclast,mpl,ix,jx,mxr,n1,jm,jfor,lda,ldb
integer :: jabs,jnb,kx,ikx,itm,inext,in
integer, dimension(3,2) :: ja
integer, dimension(icount/2-1) :: lx
real(wp_) :: drgrd,dzgrd,ah,adn,px,x,y
real(wp_), dimension(nr*nz) :: a
logical :: flag1
px = 0.5_wp_
a = reshape(matr2dgrid, [nr*nz])
rcon = 0.0_wp_
zcon = 0.0_wp_
nrqmax = nr
drgrd = rqgrid(2) - rqgrid(1)
dzgrd = zqgrid(2) - zqgrid(1)
ncon = 0
npts = 0
iclast = 0
icount = 0
mpl = 0
ix = 0
mxr = nrqmax * (nz - 1)
n1 = nr - 1
do jx=2,n1
do jm=jx,mxr,nrqmax
j = jm + nrqmax
ah=a(j)-h
if (ah <= 0.0_wp_ .and. a(jm) > h .or. &
ah > 0.0_wp_ .and. a(jm) <= h) then
ix=ix+1
lx(ix)=-j
end if
if (ah <= 0.0_wp_ .and. a(j-1) > h .or. &
ah > 0.0_wp_ .and. a(j-1) <= h) then
ix=ix+1
lx(ix)=j
end if
end do
end do
do jm=nr,mxr,nrqmax
j = jm + nrqmax
ah=a(j)-h
if (ah <= 0.0_wp_ .and. a(j-1) > h .or. &
ah > 0.0_wp_ .and. a(j-1) <= h) then
ix=ix+1
lx(ix)=j
end if
if (ah <= 0.0_wp_ .and. a(jm) > h .or. &
ah > 0.0_wp_ .and. a(jm) <= h) then
ix=ix+1
lx(ix)=-j
end if
end do
do jm=1,mxr,nrqmax
j = jm + nrqmax
if (a(j) <= h .and. a(jm) > h .or. &
a(j) > h .and. a(jm) <= h) then
ix=ix+1
lx(ix) =-j
end if
end do
do j=2,nr
if (a(j) <= h .and. a(j-1) > h .or. &
a(j) > h .and. a(j-1) <= h) then
ix=ix+1
lx(ix)=j
end if
end do
if(ix<=0) return
bb: do
in=ix
jx=lx(in)
jfor=0
lda=1
ldb=2
do
if(jx<0) then
jabs=-jx
jnb = jabs - nrqmax
else
jabs=jx
jnb=jabs-1
end if
adn=a(jabs)-a(jnb)
if(adn/=0) px=(a(jabs)-h)/adn
kx = (jabs - 1) / nrqmax
ikx = jabs - nrqmax * kx - 1
if(jx<0) then
x = drgrd * ikx
y = dzgrd * (kx - px)
else
x = drgrd * (ikx - px)
y = dzgrd * kx
end if
icount = icount + 1
rcon(icount) = x + rqgrid(1)
zcon(icount) = y + zqgrid(1)
mpl= icount
itm = 1
ja(1,1) = jabs + nrqmax
j=1
if(jx<=0) then
ja(1,1) = -jabs-1
j=2
end if
ja(2,1) = -ja(1,1)
ja(3,1) = -jx + 1 - nrqmax
ja(3,2) = -jx
ja(j,2) = jabs - nrqmax
k= 3-j
ja(k,2) = 1-jabs
if (kx<=0 .or. ikx<=0) then
lda=1
ldb=lda
else if (ikx + 1 - nr >= 0 .and. jx <= 0) then
lda=2
ldb=lda
else if(jfor/=0) then
lda=2
do i=1,3
if(jfor==ja(i,2)) then
lda=1
exit
end if
end do
ldb=lda
end if
flag1=.false.
aa: do k=1,3
do l=lda,ldb
do i=1,ix
if(lx(i)==ja(k,l)) then
itm=itm+1
inext= i
if(jfor/=0) exit aa
if(itm > 3) then
flag1=.true.
exit aa
end if
end if
end do
end do
end do aa
if(.not.flag1) then
lx(in)=0
if(itm == 1) exit
end if
jfor=jx
jx=lx(inext)
in = inext
end do
do
if(lx(ix)/=0) then
if(mpl>=4) then
ncon = ncon + 1
npts(ncon) = icount - iclast
iclast = icount
end if
exit
end if
ix= ix-1
if(ix<=0) exit bb
end do
end do bb
if(mpl >= 4) then
ncon = ncon + 1
npts(ncon) = icount - iclast
iclast = icount
end if
end subroutine cniteq
subroutine print_headers(params)
! Prints the headers of all gray_main output tables
use units, only : uprj0, uwbm, udisp, ucenr, uoutr, upec, usumm, &
unit_active, active_units
use gray_params, only : gray_parameters, headw, headl, print_parameters
! subroutine arguments
type(gray_parameters), intent(in) :: params
! local variables
integer :: i, j
integer :: main_units(8)
character(256) :: main_headers(8)
character(len=headw), dimension(headl) :: header
main_units = [uprj0, uprj0+1, uwbm, udisp, ucenr, uoutr, upec, usumm]
main_headers(1) = '#sst j k xt yt zt rt'
main_headers(2) = '#sst j k xt yt zt rt'
main_headers(3) = '#sst w1 w2'
main_headers(4) = '#sst Dr_Nr Di_Nr'
main_headers(5) = '#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'
main_headers(6) = '#i k sst x y R z psin tau Npl alpha index_rt'
main_headers(7) = '#rhop rhot Jphi Jcdb dPdV Icdins Pins index_rt'
main_headers(8) = '#Icd Pa Jphip dPdVp rhotj rhotjava rhotp rhotpav ' &
// 'drhotjava drhotpav ratjamx ratjbmx stmx psipol ' &
// 'chipol index_rt Jphimx dPdVmx drhotj drhotp P0 ' &
// 'cplO cplX'
if (all(active_units == 0)) return
call print_parameters(params, header)
do i=1,size(main_units)
if (unit_active(main_units(i))) then
do j=1,size(header)
write (main_units(i), '(1x,a)') header(j)
end do
write (main_units(i), '(1x,a)') trim(main_headers(i))
end if
end do
end subroutine print_headers
subroutine print_prof(params)
! Prints the (input) plasma kinetic profiles (unit 55)
use gray_params, only : profiles_parameters
use equilibrium, only : fq, frhotor
use coreprofiles, only : density, temp
use units, only : uprfin, unit_active
use magsurf_data, only : npsi, vajphiav
! suborutine arguments
type(profiles_parameters), intent(in) :: params
! local variables
integer :: i, ntail
real(wp_) :: rho_t, J_phi, dens, ddens
real(wp_) :: psi_n, rho_p, drho_p
if (.not. unit_active(uprfin)) return
! Δρ_p for the uniform ρ_p grid
! Note: we don't use ψ_n because J_phi has been
! sampled uniformly in ρ_p (see flux_average)
drho_p = 1.0_wp_ / (npsi - 1)
! extra points to reach ψ=ψ_bnd
ntail = ceiling((sqrt(params%psnbnd) - 1) / drho_p)
write (uprfin, *) '#psi rhot ne Te q Jphi'
do i = 0, npsi + ntail
rho_p = i * drho_p
rho_t = frhotor(rho_p)
psi_n = rho_p**2
if (psi_n < 1) then
J_phi = vajphiav(i+1) * 1.e-6_wp_
else
J_phi = 0
end if
call density(psi_n, dens, ddens)
write (uprfin, '(12(1x,e12.5))') &
psi_n, rho_t, dens, temp(psi_n), fq(psi_n), J_phi
end do
end subroutine print_prof
subroutine print_bres(bres)
! Prints the EC resonance surface table (unit 70)
use const_and_precisions, only : comp_eps
use equilibrium, only : rmnm, rmxm, zmnm, zmxm, bfield
use units, only : ubres, unit_active
use magsurf_data, only : npsi
! subroutine arguments
real(wp_), intent(in) :: bres
! local constants
integer, parameter :: icmx = 2002
! local variables
integer :: j,k,n,nconts,nctot
integer, dimension(10) :: ncpts
real(wp_) :: dr,dz,btmx,btmn,zzk,rrj,bbphi,bbr,bbz,bbb
real(wp_), dimension(icmx) :: rrcb,zzcb
real(wp_) :: rv(npsi), zv(npsi)
real(wp_), dimension(npsi,npsi) :: btotal
if (.not. unit_active(ubres)) return
! Build a regular (R, z) grid
dr = (rmxm - rmnm - comp_eps)/(npsi - 1)
dz = (zmxm - zmnm)/(npsi - 1)
do j=1,npsi
rv(j) = comp_eps + rmnm + dr*(j - 1)
zv(j) = zmnm + dz*(j - 1)
end do
! Btotal on psi grid
btmx = -1.0e30_wp_
btmn = 1.0e30_wp_
do k = 1, npsi
zzk = zv(k)
do j = 1, npsi
rrj = rv(j)
call bfield(rrj, zzk, bbr, bbz, bbphi)
btotal(j,k) = sqrt(bbr**2 + bbz**2 + bbphi**2)
if(btotal(j,k) >= btmx) btmx = btotal(j,k)
if(btotal(j,k) <= btmn) btmn = btotal(j,k)
end do
end do
! compute Btot=Bres/n with n=1,5
write (ubres, *) '#i Btot R z'
do n = 1, 5
bbb = bres/n
if (bbb >= btmn .and. bbb <= btmx) then
nconts = size(ncpts)
nctot = size(rrcb)
call cniteq(rv, zv, btotal, npsi, npsi, bbb, &
nconts, ncpts, nctot, rrcb, zzcb)
do j = 1, nctot
write (ubres, '(i6,12(1x,e12.5))') j, bbb, rrcb(j), zzcb(j)
end do
end if
write (ubres, *)
end do
end subroutine print_bres
subroutine print_maps(B_res, xgcn, R0, Npl0)
! Prints several input quantities on a regular grid (unit 72)
use const_and_precisions, only : comp_eps
use equilibrium, only : rmnm, rmxm, zmnm, zmxm, pol_flux, bfield
use coreprofiles, only : density, temp
use units, only : umaps, unit_active
use magsurf_data, only : npsi
! subroutine arguments
real(wp_), intent(in) :: B_res ! resonant magnetic field, e/m_eω
real(wp_), intent(in) :: xgcn ! X normalisation, e²/ε₀m_eω²
real(wp_), intent(in) :: R0 ! initial value of R
real(wp_), intent(in) :: Npl0 ! initial value of N∥
! local variables
integer :: j, k
real(wp_) :: dR, dz, B_R, B_phi, B_z, B
real(wp_) :: psi_n, ne, dne, Te, X, Y, Npl
real(wp_), dimension(npsi) :: R, z
if (.not. unit_active(umaps)) return
! Build a regular (R, z) grid
dR = (rmxm - rmnm - comp_eps)/(npsi - 1)
dz = (zmxm - zmnm)/(npsi - 1)
do j = 1, npsi
R(j) = comp_eps + rmnm + dR*(j - 1)
z(j) = zmnm + dz*(j - 1)
end do
write (umaps, *) '#R z psin Br Bphi Bz Btot ne Te X Y Npl'
do j = 1, npsi
Npl = Npl0 * R0/r(j)
do k = 1, npsi
call pol_flux(r(j), z(k), psi_n)
call bfield(r(j), z(k), B_R, B_z, B_phi)
call density(psi_n, ne, dne)
B = sqrt(B_R**2 + B_phi**2 + B_z**2)
X = xgcn*ne
Y = B/B_res
Te = temp(psi_n)
write (umaps,'(12(x,e12.5))') &
R(j), z(k), psi_n, B_R, B_phi, B_z, B, ne, Te, X, Y, Npl
end do
write (umaps,*)
end do
end subroutine print_maps
subroutine print_surfq(qvals)
! Prints the constant ψ surfaces for a set of values of q
use equilibrium, only : fq, frhotor, rmaxis, zmaxis, zbsup, zbinf
2021-12-18 18:57:38 +01:00
use magsurf_data, only : contours_psi, npoints, print_contour
use logger, only : log_info
use minpack, only : hybrj1
2021-12-18 18:57:38 +01:00
! subroutine arguments
real(wp_), intent(in) :: qvals(:)
2021-12-18 18:57:38 +01:00
! local variables
integer :: i
real(wp_) :: rho_t, rho_p, psi_n
2021-12-18 18:57:38 +01:00
character(256) :: msg ! for log messages formatting
2021-12-18 18:57:38 +01:00
call log_info('constant ψ surfaces for:', &
mod='gray_core', proc='print_surfq')
do i = 1, size(qvals)
! Find value of ψ_n for the current q
block
real(wp_) :: sol(1), fvec(1), fjac(1,1), wa(7)
integer :: info
! Note: since we are mostly interested in q=3/2,2 we start
! searching near the boundary in case q is not monotonic.
sol = [0.8_wp_] ! first guess
! Solve fq(ψ_n) = qvals(i) for ψ_n
call hybrj1(equation, n=1, x=sol, fvec=fvec, fjac=fjac, &
ldfjac=1, tol=1e-3_wp_, info=info, wa=wa, lwa=7)
! Solution
psi_n = sol(1)
! Handle spurious or no solutions
if (info /= 1 .or. psi_n < 0 .or. psi_n > 1) cycle
end block
! Compute and print the ψ_n(R,z) = ψ_n₀ contour
block
real(wp_), dimension(npoints) :: R_cont, z_cont
real(wp_) :: R_hi, R_lo, z_hi, z_lo
! Guesses for the high,low horzizontal-tangent points
R_hi = rmaxis;
z_hi = (zbsup + zmaxis)/2
R_lo = rmaxis
z_lo = (zbinf + zmaxis)/2
call contours_psi(psi_n, R_cont, z_cont, R_hi, z_hi, R_lo, z_lo)
call print_contour(psi_n, R_cont, z_cont)
end block
! Log the values found for ψ_n, ρ_p, ρ_t
rho_p = sqrt(psi_n)
rho_t = frhotor(rho_p)
write (msg, '(4(x,a,"=",g0.3))') &
'q', qvals(i), 'ψ_n', psi_n, 'ρ_p', rho_p, 'ρ_t', rho_t
call log_info(msg, mod='gray_core', proc='print_surfq')
end do
contains
subroutine equation(n, x, f, df, ldf, flag)
! The equation to solve: f(x) = q(x) - q₀ = 0
use const_and_precisions, only : comp_eps
! subroutine arguments
integer, intent(in) :: n, ldf, flag
real(wp_), intent(in) :: x(n)
real(wp_), intent(inout) :: f(n), df(ldf,n)
! optimal step size
real(wp_), parameter :: e = comp_eps**(1/3.0_wp_)
if (flag == 1) then
! return f(x)
f(1) = fq(x(1)) - qvals(i)
else
! return f'(x), computed numerically
df(1,1) = (fq(x(1) + e) - fq(x(1) - e)) / (2*e)
end if
end subroutine
2021-12-18 18:57:38 +01:00
end subroutine print_surfq
subroutine print_projxyzt(stv, ywrk, iproj)
! Prints the beam shape tables (unit 8, 12)
use const_and_precisions, only : comp_huge, zero, one
use beamdata, only : nray, nrayr, nrayth, rayi2jk
use units, only : uprj0, uwbm, unit_active
! subroutine arguments
real(wp_), dimension(:), intent(in) :: stv
real(wp_), dimension(:,:), intent(in) :: ywrk
integer, intent(in) :: iproj
! local variables
integer :: jk, jkz, uprj
integer, dimension(2) ::jkv
real(wp_), dimension(3) :: xv1, dir, dxv
real(wp_) :: rtimn, rtimx, csth1, snth1, &
csps1, snps1, xti, yti, zti, rti
if (.not.(unit_active(uprj0) .or. unit_active(uwbm))) return
! common/external functions/variables
uprj = uprj0 + iproj
xv1 = ywrk(1:3, 1)
dir = ywrk(4:6, 1)
dir = dir/norm2(dir)
csth1 = dir(3)
snth1 = sqrt(one - csth1**2)
if(snth1 > zero) then
csps1 = dir(2)/snth1
snps1 = dir(1)/snth1
else
csps1 = one
snps1 = zero
end if
if(iproj==0) then
jkz = nray - nrayth + 1
else
jkz = 1
end if
rtimn = comp_huge
rtimx = zero
do jk = jkz, nray
dxv = ywrk(1:3,jk) - xv1
xti = dxv(1)*csps1 - dxv(2)*snps1
yti =(dxv(1)*snps1 + dxv(2)*csps1)*csth1 - dxv(3)*snth1
zti =(dxv(1)*snps1 + dxv(2)*csps1)*snth1 + dxv(3)*csth1
rti = hypot(xti, yti)
jkv=rayi2jk(jk)
if(.not.(iproj==0 .and. jk==1)) &
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))') stv(1), rtimn, rtimx
end subroutine print_projxyzt
subroutine print_output(i, jk, st, qj, xv, psinv, btot, bv, ak0, &
anpl, anpr, anv, anprim, dens, tekev, &
alpha, tau, dids, nhm, nhf, iokhawa, &
index_rt, ddr, ddi, xg, yg, derxg)
! Prints the ray variables and Hamiltonian error (units 4, 33, 17)
2015-11-18 17:34:33 +01:00
use const_and_precisions, only : degree,zero,one
use equilibrium, only : frhotor
use gray_params, only : istpl0
use beamdata, only : nray,nrayth,jkray1
use units, only : ucenr, uoutr, udisp, unit_active
! subroutine arguments
integer, intent(in) :: i, jk, nhm, nhf, iokhawa, index_rt
real(wp_), dimension(3), intent(in) :: xv, bv, anv
real(wp_), intent(in) :: st, qj, psinv, btot, ak0, anpl, anpr, anprim
real(wp_), intent(in) :: dens, tekev, alpha, tau, dids, ddr, ddi
real(wp_), intent(in) :: xg, yg
2017-09-12 21:37:06 +02:00
real(wp_), dimension(3), intent(in) :: derxg
! local variables
real(wp_) :: stm, xxm, yym, zzm, rrm, phideg, rhot, akim, pt, didsn
2015-11-18 17:34:33 +01:00
integer :: k
stm = st * 1.0e-2_wp_
xxm = xv(1) * 1.0e-2_wp_
yym = xv(2) * 1.0e-2_wp_
zzm = xv(3) * 1.0e-2_wp_
rrm = hypot(xxm, yym)
2015-11-18 17:34:33 +01:00
! print central ray trajectory. dIds in A/m/W, ki in m⁻¹
if(jk == 1) then
phideg = atan2(yym, xxm)/degree
2015-11-18 17:34:33 +01:00
if(psinv>=zero .and. psinv<=one) then
rhot = frhotor(sqrt(psinv))
2015-11-18 17:34:33 +01:00
else
rhot = 1.0_wp_
end if
akim = anprim * ak0 * 1.0e2_wp_
pt = qj*exp(-tau)
didsn = dids * 1.0e2_wp_/qj
! print central ray data
if (unit_active(ucenr)) then
write (ucenr,'(22(1x,e16.8e3),4i5,6(1x,e16.8e3))') &
stm, rrm, zzm, phideg, psinv, rhot, dens, tekev, &
btot, bv, anpr, anpl, anv, akim, alpha, tau, pt, &
didsn, nhm, nhf, iokhawa, index_rt, ddr, xg, yg, &
derxg
2015-11-18 17:34:33 +01:00
end if
end if
! print conservation of dispersion relation (Hamiltonian error)
if(unit_active(udisp) .and. jk==nray) then
write(udisp,'(30(1x,e16.8e3))') st,ddr,ddi
end if
2015-11-18 17:34:33 +01:00
! print outer ray trajectories
if(mod(i, istpl0) == 0) then
k = jk - jkray1 + 1
if(unit_active(uoutr) .and. k>0 .and. k<=nrayth) then
write(uoutr, '(2i5,9(1x,e16.8e3),i5)') &
i, k, stm, xxm, yym, rrm, zzm, &
psinv, tau, anpl, alpha, index_rt
2015-11-18 17:34:33 +01:00
end if
end if
end subroutine print_output
subroutine print_pec(rhop_tab, rhot_tab, jphi, jcd, &
dpdv, currins, pins, index_rt)
! Prints the ECRH and CD profiles (unit 48)
use units, only : upec, unit_active
! subroutine arguments
real(wp_), dimension(:), intent(in) :: rhop_tab, rhot_tab, jphi, jcd, &
dpdv, currins, pins
integer, intent(in) :: index_rt
! local variables
integer :: i
if (.not. unit_active(upec)) return
do i=1, size(rhop_tab)
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, p0, cpl1, cpl2)
! Prints global data and final results (unit 7)
use units, only : ucenr, usumm, unit_active
! subroutine arguments
real(wp_), intent(in) :: pabs, icd, dpdvp, jphip, rhotpav, rhotjava, &
drhotpav, drhotjava, dpdvmx, jphimx, rhotp, &
rhotj, drhotp, drhotj, ratjamx, ratjbmx, &
stmx, psipol, chipol, p0, cpl1, cpl2
integer, intent(in) :: index_rt
if (unit_active(ucenr)) then
write(ucenr, *) ''
end if
if (unit_active(usumm)) then
write(usumm,'(15(1x,e12.5),i5,7(1x,e12.5))') &
icd, pabs, jphip, dpdvp, rhotj, rhotjava, rhotp, rhotpav, &
drhotjava, drhotpav, ratjamx, ratjbmx, stmx, psipol, &
chipol, index_rt, jphimx, dpdvmx, drhotj, drhotp, &
p0, cpl1, cpl2
end if
end subroutine print_finals
end module gray_core