2274 lines
82 KiB
Fortran
2274 lines
82 KiB
Fortran
module gray_core
|
||
use const_and_precisions, only : wp_
|
||
implicit none
|
||
|
||
contains
|
||
|
||
subroutine gray_main(params, data, results, error, rhout)
|
||
use const_and_precisions, only : zero, one, degree, comp_tiny
|
||
use coreprofiles, only : set_prfan, set_prfspl, temp, fzeff, unset_prfspl
|
||
use dispersion, only : expinit
|
||
use gray_params, only : gray_parameters, gray_data, gray_results, print_parameters, &
|
||
iwarm, ipec, istpr0, igrad, headw, headl, ipass
|
||
use beams, only : xgygcoeff, launchangles2n
|
||
use beamdata, only : pweight, rayi2jk
|
||
use equilibrium, only : unset_eqspl, unset_rhospl, unset_q
|
||
use errcodes, only : check_err, print_errn, print_errhcd
|
||
use magsurf_data, only : flux_average, dealloc_surfvec
|
||
use beamdata, only : init_btr, dealloc_beam, nray, nstep, dst
|
||
use pec, only : pec_init, spec, postproc_profiles, dealloc_pec, &
|
||
rhop_tab, rhot_tab
|
||
use limiter, only : limiter_unset_globals=>unset_globals
|
||
use utils, only : vmaxmin
|
||
use reflections, only : inside
|
||
use multipass, only : alloc_multipass, dealloc_multipass, initbeam, &
|
||
initmultipass, turnoffray, plasma_in, plasma_out, wall_out
|
||
use units, only : ucenr
|
||
use logger, only : log_info, log_debug
|
||
|
||
implicit none
|
||
|
||
! Subroutine arguments
|
||
type(gray_parameters), intent(in) :: 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
|
||
|
||
! Exit code
|
||
integer, intent(out) :: error
|
||
|
||
! local variables
|
||
real(wp_), parameter :: taucr = 12._wp_, etaucr = exp(-taucr)
|
||
character, dimension(2), parameter :: mode=(/'O','X'/)
|
||
|
||
real(wp_) :: sox,ak0,bres,xgcn,xg,yg,rrm,zzm,alpha,didp,anpl,anpr,anprim,anprre
|
||
real(wp_) :: chipol,psipol,btot,psinv,dens,tekev,dersdst,derdnm
|
||
real(wp_) :: tau,pow,dids,ddr,ddi,taumn,taumx
|
||
real(wp_) :: rhotpav,drhotpav,rhotjava,drhotjava,dpdvp,jphip
|
||
real(wp_) :: rhotp,drhotp,rhotj,drhotj,dpdvmx,jphimx,ratjamx,ratjbmx
|
||
real(wp_) :: pabs_beam,icd_beam,cpl_beam1,cpl_beam2,cpl_cbeam1,cpl_cbeam2
|
||
|
||
real(wp_), dimension(2) :: pabs_pass,icd_pass,cpl,cpl0
|
||
real(wp_), dimension(3) :: xv,anv0,anv,bv,derxg
|
||
|
||
! Ray variables
|
||
real(wp_), dimension(:,:), pointer :: yw=>null(),ypw=>null(),gri=>null()
|
||
real(wp_), dimension(:,:,:), pointer :: xc=>null(),du1=>null(),ggri=>null()
|
||
|
||
! i: integration step, jk: global ray index
|
||
integer :: i, jk
|
||
|
||
integer :: iox,nharm,nhf,nnd,iokhawa,istop,ierrn,ierrhcd,index_rt
|
||
integer :: ip,ib,iopmin,ipar,iO
|
||
integer :: igrad_b,istop_pass,nbeam_pass,nlim
|
||
logical :: ins_pl,ins_wl,ent_pl,ext_pl,ent_wl,ext_wl,iboff
|
||
|
||
real(wp_), dimension(:,:,:), pointer :: yynext=>null(),yypnext=>null()
|
||
real(wp_), dimension(:,:), pointer :: psjki=>null(),ppabs=>null(),ccci=>null()
|
||
real(wp_), dimension(:,:), pointer :: taus=>null(),stnext=>null(), &
|
||
yw0=>null(),ypw0=>null(),cpls=>null()
|
||
real(wp_), dimension(:), pointer :: p0ray=>null(),tau0=>null(),alphaabs0=>null(), &
|
||
dids0=>null(),ccci0=>null(),tau1=>null(),etau1=>null(),cpl1=>null(),lgcpl1=>null()
|
||
real(wp_), dimension(:), pointer :: p0jk=>null()
|
||
real(wp_), dimension(:), pointer :: jphi_beam=>null(),pins_beam=>null(), &
|
||
currins_beam=>null(), dpdv_beam=>null(),jcd_beam=>null(),stv=>null(), &
|
||
psipv=>null(),chipv=>null()
|
||
complex(wp_), dimension(:), pointer :: ext=>null(), eyt=>null()
|
||
integer, dimension(:), pointer :: iiv=>null(),iop=>null(),iow=>null()
|
||
logical, dimension(:), pointer :: iwait=>null()
|
||
logical, dimension(:,:), pointer :: iroff=>null()
|
||
|
||
! parameters log in file headers
|
||
character(len=headw), dimension(headl) :: strheader
|
||
! buffer for formatting log messages
|
||
character(256) :: msg
|
||
|
||
! ======== set environment BEGIN ========
|
||
! Number of limiter contourn points
|
||
nlim = size(data%equilibrium%zlim)
|
||
|
||
! Compute X=ω/ω_ce and Y=(ω/ω_pe)² (with B=1)
|
||
call xgygcoeff(params%antenna%fghz, ak0, bres, xgcn)
|
||
|
||
! Compute the initial cartesian wavevector (anv0)
|
||
! from launch angles α,β and the position x₀:
|
||
! NR(α, β, x₀)
|
||
! Nφ(α, β, x₀)
|
||
! Nz(α, β, x₀)
|
||
call launchangles2n(params%antenna, anv0)
|
||
|
||
! Initialise the ray variables (beamtracing)
|
||
call init_btr(params%raytracing, yw, ypw, xc, du1, gri, ggri, psjki, ppabs, ccci, &
|
||
tau0, alphaabs0, dids0, ccci0, p0jk, ext, eyt, iiv)
|
||
|
||
! Initialise the dispersion module
|
||
if(iwarm > 1) call expinit
|
||
|
||
! Initialise the magsurf_data module
|
||
call flux_average ! requires frhotor for dadrhot,dvdrhot
|
||
|
||
! Initialise the output profiles
|
||
call pec_init(ipec, rhout)
|
||
nnd = size(rhop_tab) ! number of radial profile points
|
||
|
||
call alloc_multipass(nnd, iwait, iroff, iop, iow, yynext, yypnext, yw0, ypw0, stnext, &
|
||
stv, p0ray, taus, tau1, etau1, cpls, cpl1, lgcpl1, jphi_beam, &
|
||
pins_beam, currins_beam, dpdv_beam, jcd_beam, psipv, chipv)
|
||
|
||
! 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_parameters(params, strheader)
|
||
call print_headers(strheader)
|
||
|
||
! 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
|
||
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
|
||
call print_bres(bres)
|
||
call print_prof
|
||
call print_maps(bres, xgcn, &
|
||
0.01_wp_*sqrt(params%antenna%pos(1)**2 + params%antenna%pos(2)**2), &
|
||
sin(params%antenna%beta*degree))
|
||
! ========= 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)
|
||
|
||
if(params%raytracing%ipol .eq. 0) then
|
||
if(params%antenna%iox .eq. 2) then ! only X mode on 1st pass
|
||
cpl0 = (/zero,one/)
|
||
else ! only O mode on 1st pass
|
||
cpl0 = (/one,zero/)
|
||
end if
|
||
end if
|
||
|
||
sox=one ! mode inverted for each beam
|
||
iox=2 ! start with O: sox=-1, iox=1
|
||
|
||
psipol = params%antenna%psi
|
||
chipol = params%antenna%chi
|
||
call pweight(params%antenna%power, p0jk)
|
||
|
||
nbeam_pass=1 ! max n of beam per pass
|
||
index_rt=0 ! global beam index: 1,O 2,X 1st pass
|
||
! | | | |
|
||
call log_debug('pass loop start', mod='gray_core', proc='gray_main') ! 3,O 4,X 5,O 6,X 2nd pass
|
||
do ip=1,ipass
|
||
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.gt.1) then
|
||
du1 = zero
|
||
gri = zero
|
||
ggri = zero
|
||
if(ip.eq.ipass) cpl = (/zero,zero/) ! no successive passes
|
||
end if
|
||
|
||
! =========== beam loop BEGIN ===========
|
||
call log_debug('beam loop start', mod='gray_core', proc='gray_main')
|
||
do ib=1,nbeam_pass
|
||
|
||
sox = -sox ! invert mode
|
||
iox = 3-iox ! O-mode at ip=1,ib=1
|
||
|
||
index_rt = index_rt +1
|
||
iO = 2*index_rt +1 ! * index_rt of O-mode derived ray (iX=iO+1)
|
||
|
||
call initbeam(index_rt,iroff,iboff,iwait,stv,jphi_beam, &
|
||
pins_beam,currins_beam,dpdv_beam,jcd_beam)
|
||
|
||
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
|
||
call log_info(" beam is off", mod='gray_core', proc='gray_main')
|
||
cycle
|
||
end if
|
||
|
||
call vectinit(psjki,ppabs,ccci,tau0,alphaabs0,dids0,ccci0,iiv)
|
||
|
||
if(ip.eq.1) then ! 1st pass
|
||
igrad_b = igrad ! * input value, igrad_b=0 from 2nd pass
|
||
|
||
tau1 = zero ! * tau from previous passes
|
||
etau1 = one
|
||
cpl1 = one ! * coupling from previous passes
|
||
lgcpl1 = zero
|
||
p0ray = p0jk ! * initial beam power
|
||
|
||
call ic_gb(params%antenna%pos, anv0, ak0, &
|
||
params%antenna%w(1),params%antenna%w(2), &
|
||
params%antenna%ri(1),params%antenna%ri(2), &
|
||
params%antenna%phi(1),params%antenna%phi(2), &
|
||
yw,ypw,xc,du1,gri,ggri,index_rt) ! * initial conditions
|
||
call set_pol(yw,bres,sox,psipol,chipol,ext,eyt) ! * initial polarization
|
||
|
||
do jk=1,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, &
|
||
nlim, rrm, zzm)) then ! * start propagation in/outside vessel?
|
||
iow(jk) = 1 ! + inside
|
||
else
|
||
iow(jk) = 0 ! + outside
|
||
end if
|
||
end do
|
||
|
||
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(nray>1 .and. all(.not.iwait)) call print_projxyzt(stv,yw,0) ! iproj=0 ==> nfilp=8
|
||
|
||
! ======= propagation loop BEGIN =======
|
||
call log_debug(' propagation loop start', mod='gray_core', proc='gray_main')
|
||
|
||
do i=1,nstep
|
||
! advance one step with "frozen" grad(S_I)
|
||
do jk=1,nray
|
||
if(iwait(jk)) cycle ! jk ray is waiting for next pass
|
||
stv(jk) = stv(jk) + dst ! current ray step
|
||
call rkstep(sox,bres,xgcn,yw(:,jk),ypw(:,jk),gri(:,jk),ggri(:,:,jk),igrad_b)
|
||
end do
|
||
! update position and grad
|
||
if(igrad_b == 1) call gradi_upd(yw,ak0,xc,du1,gri,ggri)
|
||
|
||
error = 0
|
||
istop = 0 ! stop flag for current beam
|
||
iopmin = 10
|
||
|
||
! =========== rays loop BEGIN ===========
|
||
do jk=1,nray
|
||
if(iwait(jk)) cycle ! jk ray is waiting for next pass
|
||
|
||
! compute derivatives with updated gradient and local plasma values
|
||
xv = yw(1:3,jk)
|
||
anv = yw(4:6,jk)
|
||
call ywppla_upd(xv,anv,gri(:,jk),ggri(:,:,jk),sox,bres,xgcn,ypw(:,jk), &
|
||
psinv,dens,btot,bv,xg,yg,derxg,anpl,anpr,ddr,ddi,dersdst,derdnm, &
|
||
ierrn,igrad_b)
|
||
! update global error code and print message
|
||
if(ierrn/=0) then
|
||
error = ior(error,ierrn)
|
||
call print_errn(ierrn,i,anpl)
|
||
end if
|
||
|
||
! check entrance/exit plasma/wall
|
||
zzm = xv(3)*0.01_wp_
|
||
rrm = sqrt(xv(1)*xv(1)+xv(2)*xv(2))*0.01_wp_
|
||
|
||
ins_pl = (psinv>=zero .and. psinv<params%profiles%psnbnd) ! in/out plasma?
|
||
ins_wl = inside(data%equilibrium%rlim, data%equilibrium%zlim, &
|
||
nlim,rrm,zzm) ! in/out vessel?
|
||
ent_pl = (mod(iop(jk),2).eq.0 .and. ins_pl) ! enter plasma
|
||
ext_pl = (mod(iop(jk),2).eq.1 .and. .not.ins_pl) ! exit plasma
|
||
ent_wl = (mod(iow(jk),2).eq.0 .and. ins_wl) ! enter vessel
|
||
ext_wl = (mod(iow(jk),2).eq.1 .and. .not.ins_wl) ! exit vessel
|
||
|
||
if(ent_pl) then ! ray enters plasma
|
||
call plasma_in(jk,xv,anv,bres,sox,cpl,psipol,chipol,iop,ext,eyt)
|
||
|
||
if(iop(jk).eq.1 .and. ip==1) then ! * 1st entrance on 1st pass (ray hasn't entered in plasma yet) => continue current pass
|
||
|
||
if(params%raytracing%ipol .eq. 0) then ! + IF single mode propagation
|
||
cpl = cpl0
|
||
p0ray(jk) = p0ray(jk)*cpl(iox)
|
||
else if(cpl(iox).lt.etaucr) then ! + ELSE IF low coupled power for current mode => de-activate derived rays
|
||
call turnoffray(jk,ip+1,2*ib+2-iox,iroff)
|
||
iwait(jk) = .true. ! . stop advancement and H&CD computation for current ray
|
||
if(cpl(iox).le.comp_tiny) cpl(iox)=etaucr
|
||
else ! + ELSE assign coupled power to current ray
|
||
p0ray(jk) = p0ray(jk)*cpl(iox)
|
||
end if
|
||
cpls(jk,index_rt) = cpl(iox)
|
||
|
||
if(jk.eq.1) then
|
||
write (msg,'(" 1st pass - central ray (",a1,"-mode) c=",g0.4)') &
|
||
mode(iox), cpl(iox)
|
||
call log_info(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).gt.2) then ! * 2nd entrance on 1st pass / entrance on 2nd+ pass => end of current pass for ray jk
|
||
igrad_b = 0 ! + switch to ray-tracing
|
||
iwait(jk) = .true. ! + stop advancement and H&CD computation for current ray
|
||
|
||
if(ip.lt.ipass) then ! + not last pass
|
||
yynext(:,jk,index_rt) = yw0(:,jk) ! . copy starting coordinates
|
||
yypnext(:,jk,index_rt) = ypw0(:,jk) ! for next pass from last step
|
||
stnext(jk,index_rt) = stv(jk) - dst ! . starting step for next pass = last step
|
||
|
||
if(cpl(1).lt.etaucr) then ! . low coupled power for O-mode => de-activate derived rays
|
||
call turnoffray(jk,ip+1,2*ib-1,iroff)
|
||
if(cpl(1).le.comp_tiny) cpl(1)=etaucr
|
||
end if
|
||
if(cpl(2).lt.etaucr) then ! . low coupled power for X-mode => de-activate derived rays
|
||
call turnoffray(jk,ip+1,2*ib,iroff)
|
||
if(cpl(2).le.comp_tiny) cpl(2)=etaucr
|
||
end if
|
||
|
||
taus(jk,iO:iO+1) = tau1(jk) + tau0(jk) ! . starting tau for next O-mode pass
|
||
cpls(jk,iO) = cpl1(jk) * cpl(1) ! . cumulative coupling for next O-mode pass
|
||
cpls(jk,iO+1) = cpl1(jk) * cpl(2) ! . cumulative coupling for next X-mode pass
|
||
|
||
if(jk.eq.1) then ! . polarization angles at plasma boundary for central ray
|
||
psipv(iO:iO+1) = psipol
|
||
chipv(iO:iO+1) = chipol
|
||
end if
|
||
else ! * 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
|
||
call plasma_out(jk,xv,anv,bres,sox,iop,ext,eyt)
|
||
end if
|
||
|
||
if(ent_wl) then ! ray enters vessel
|
||
iow(jk)=iow(jk)+1 ! * 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_errn(ierrn,i,anpl)
|
||
end if
|
||
|
||
if(jk.eq.1 .and. ip.eq.1) then ! * 1st pass, polarization angles at reflection for central ray
|
||
psipv(index_rt) = psipol
|
||
chipv(index_rt) = chipol
|
||
end if
|
||
|
||
if(ins_pl) then ! * plasma-wall overlapping => wall+plasma crossing => end of current pass
|
||
iwait(jk) = .true. ! + stop advancement and H&CD computation for current ray
|
||
|
||
if(ip.lt.ipass) then ! + not last pass
|
||
yynext(:,jk,index_rt) = (/xv,anv/) ! . starting coordinates
|
||
yypnext(:,jk,index_rt) = ypw(:,jk) ! for next pass = reflection point
|
||
stnext(jk,index_rt) = stv(jk) ! . starting step for next pass = step after reflection
|
||
|
||
call plasma_in(jk,xv,anv,bres,sox,cpl,psipol,chipol,iop,ext,eyt) ! . ray re-enters plasma after reflection
|
||
|
||
if(cpl(1).lt.etaucr) then ! . low coupled power for O-mode? => de-activate derived rays
|
||
call turnoffray(jk,ip+1,2*ib-1,iroff)
|
||
if(cpl(1).le.comp_tiny) cpl(1)=etaucr
|
||
end if
|
||
if(cpl(2).lt.etaucr) then ! . low coupled power for X-mode? => de-activate derived rays
|
||
call turnoffray(jk,ip+1,2*ib,iroff)
|
||
if(cpl(2).le.comp_tiny) cpl(2)=etaucr
|
||
end if
|
||
|
||
taus(jk,iO:iO+1) = tau1(jk) + tau0(jk) ! . starting tau for next O-mode pass
|
||
cpls(jk,iO) = cpl1(jk) * cpl(1) ! . cumulative coupling for next O-mode pass
|
||
cpls(jk,iO+1) = cpl1(jk) * cpl(2) ! . cumulative coupling for next X-mode pass
|
||
|
||
if(jk.eq.1) then ! + polarization angles at plasma boundary for central ray
|
||
psipv(iO:iO+1) = psipol
|
||
chipv(iO:iO+1) = chipol
|
||
end if
|
||
end if
|
||
end if
|
||
end if
|
||
|
||
iopmin = min(iopmin,iop(jk))
|
||
if(ip.lt.ipass) then ! not last pass
|
||
yw0(:,jk) = yw(:,jk) ! * store current coordinates in case
|
||
ypw0(:,jk) = ypw(:,jk) ! current pass ends on next step
|
||
end if
|
||
|
||
! compute ECRH&CD if (inside plasma & power available>0 & ray still active)
|
||
if(ierrn==0 .and. iwarm>0 .and. ins_pl .and. &
|
||
(tau1(jk)+tau0(jk)+lgcpl1(jk))<=taucr .and. .not.iwait(jk)) then ! H&CD computation check
|
||
tekev=temp(psinv)
|
||
call alpha_effj(psinv,xg,yg,dens,tekev,ak0,bres,derdnm,anpl,anpr, &
|
||
sox,anprre,anprim,alpha,didp,nharm,nhf,iokhawa,ierrhcd)
|
||
if(ierrhcd/=0) then
|
||
error = ior(error,ierrhcd)
|
||
call print_errhcd(ierrhcd,i,anprre,anprim,alpha)
|
||
end if
|
||
else
|
||
tekev=zero
|
||
alpha=zero
|
||
didp=zero
|
||
anprim=zero
|
||
anprre=anpr
|
||
nharm=0
|
||
nhf=0
|
||
iokhawa=0
|
||
end if
|
||
if(nharm>0) iiv(jk)=i
|
||
|
||
psjki(jk,i) = psinv
|
||
! computation of optical depth tau, dP/ds, P(s), dI/ds, I(s)
|
||
tau=tau0(jk)+0.5_wp_*(alphaabs0(jk)+alpha)*dersdst*dst
|
||
pow=p0ray(jk)*exp(-tau) !*exp(-tau1v(jk))
|
||
ppabs(jk,i)=p0ray(jk)-pow
|
||
dids=didp*pow*alpha
|
||
ccci(jk,i)=ccci0(jk)+0.5_wp_*(dids0(jk)+dids)*dersdst*dst
|
||
tau0(jk)=tau
|
||
alphaabs0(jk)=alpha
|
||
dids0(jk)=dids
|
||
ccci0(jk)=ccci(jk,i)
|
||
|
||
if(iwait(jk)) then ! copy values from last pass for inactive ray
|
||
ppabs(jk,i:nstep) = ppabs(jk,i-1)
|
||
ccci(jk,i:nstep) = ccci(jk,i-1)
|
||
psjki(jk,i:nstep) = psjki(jk,i-1)
|
||
else
|
||
call print_output(i,jk,stv(jk),p0ray(jk),xv,psinv, &
|
||
btot,bv,ak0,anpl,anpr,anv,anprim,dens,tekev,alpha,tau,dids, &
|
||
nharm,nhf,iokhawa,index_rt,ddr,ddi,xg,yg,derxg) ! p0ray/etau1 [dids normalization] = fraction of p0 coupled to this ray (not including absorption from previous passes)
|
||
end if
|
||
|
||
end do
|
||
! ============ rays loop END ============
|
||
|
||
if(i==nstep) then ! step limit reached?
|
||
do jk=1,nray
|
||
if(iop(jk)<3) call turnoffray(jk,ip,ib,iroff) ! * ray hasn't exited+reentered the plasma by last step => stop ray
|
||
end do
|
||
end if
|
||
|
||
! print ray positions for j=nrayr in local reference system
|
||
if(mod(i,istpr0) == 0) then
|
||
if(nray > 1 .and. all(.not.iwait)) call print_projxyzt(stv,yw,0)
|
||
end if
|
||
|
||
! check for any error code and stop if necessary
|
||
call check_err(error,istop)
|
||
! test whether further trajectory integration is unnecessary
|
||
call vmaxmin(tau1+tau0+lgcpl1,nray,taumn,taumx) ! test on tau + coupling
|
||
! if(taumn > taucr .or. all(iroff(:,index_rt))) istop = 1 ! (residual power~0) or (no ray active) => stop beam
|
||
|
||
if(istop == 1) then ! stop propagation for current beam
|
||
istop_pass = istop_pass +1 ! * +1 non propagating beam
|
||
if(ip.lt.ipass) call turnoffray(0,ip,ib,iroff) ! * de-activate derived beams
|
||
exit
|
||
else if(all(iwait)) then ! all rays in current beam are waiting for next pass => do not increase istop_pass
|
||
exit
|
||
end if
|
||
end do
|
||
call log_debug(' propagation loop end', mod='gray_core', proc='gray_main')
|
||
! ======== propagation loop END ========
|
||
|
||
! print all ray positions in local reference system
|
||
if(nray > 1 .and. all(.not.iwait)) call print_projxyzt(stv,yw,1)
|
||
|
||
! =========== post-proc BEGIN ===========
|
||
! compute total absorbed power and driven current for current beam
|
||
if(i>nstep) i=nstep
|
||
pabs_beam = sum(ppabs(:,i))
|
||
icd_beam = sum(ccci(:,i))
|
||
call vmaxmin(tau0,nray,taumn,taumx) ! taumn,taumx for print
|
||
|
||
! compute power and current density profiles for all rays
|
||
call spec(psjki,ppabs,ccci,iiv,pabs_beam,icd_beam,dpdv_beam,jphi_beam,jcd_beam, &
|
||
pins_beam,currins_beam)
|
||
|
||
pabs_pass(iox) = pabs_pass(iox) + pabs_beam ! 0D results for current pass, sum on O/X mode beams
|
||
icd_pass(iox) = icd_pass(iox) + icd_beam
|
||
|
||
if(ip.lt.ipass .and. iopmin.gt.2) then ! not last pass AND at least one ray re-entered plasma
|
||
cpl_beam1 = sum(p0ray * exp(-tau0) * cpls(:,iO)/cpl1, MASK=iop.gt.2) / &
|
||
sum(p0ray * exp(-tau0), MASK=iop.gt.2) ! * average O-mode coupling for next beam (on active rays)
|
||
cpl_beam2 = one - cpl_beam1 ! * average X-mode coupling for next beam
|
||
|
||
if(iop(1).gt.2) then ! * central ray O/X-mode coupling for next beam
|
||
cpl_cbeam1 = cpls(1,iO)/cpl1(1)
|
||
cpl_cbeam2 = one - cpl_cbeam1
|
||
end if
|
||
else ! last pass OR no ray re-entered plasma
|
||
cpl_beam1 = zero
|
||
cpl_beam2 = zero
|
||
end if
|
||
|
||
! print final results for pass on screen
|
||
call log_info(' partial results:', mod='gray_core', proc='gray_main')
|
||
write(msg, '(3x,a,g0.4)') 'final step: (s, ct, Sr)=' ,stv(1)
|
||
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')
|
||
|
||
write(msg, '(3x,a,g0.3," MW")') 'current drive: I=', icd_beam*1.0e3_wp_
|
||
call log_info(msg, mod='gray_core', proc='gray_main')
|
||
|
||
if(ip < ipass) then
|
||
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
|
||
end if
|
||
end if
|
||
|
||
write(ucenr,*) ''
|
||
|
||
call print_pec(rhop_tab,rhot_tab,jphi_beam,jcd_beam,dpdv_beam,currins_beam, &
|
||
pins_beam,ip) ! *print power and current density profiles for current beam
|
||
|
||
call postproc_profiles(pabs_beam,icd_beam,rhot_tab,dpdv_beam,jphi_beam, &
|
||
rhotpav,drhotpav,rhotjava,drhotjava,dpdvp,jphip,rhotp,drhotp,rhotj, &
|
||
drhotj,dpdvmx,jphimx,ratjamx,ratjbmx) ! *compute profiles width for current beam
|
||
|
||
call print_finals(pabs_beam,icd_beam,dpdvp,jphip,rhotpav,rhotjava, &
|
||
drhotpav,drhotjava,dpdvmx,jphimx,rhotp,rhotj,drhotp,drhotj,ratjamx, &
|
||
ratjbmx,stv(1),psipv(index_rt),chipv(index_rt),index_rt,sum(p0ray), &
|
||
cpl_beam1,cpl_beam2) ! *print 0D results for current beam
|
||
! ============ post-proc END ============
|
||
|
||
end do
|
||
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)
|
||
|
||
! 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
|
||
call log_debug('pass loop end', mod='gray_core', proc='gray_main')
|
||
! ============ main loop END ============
|
||
|
||
! ========== free memory BEGIN ==========
|
||
call dealloc_surfvec
|
||
call dealloc_beam(yw,ypw,xc,du1,gri,ggri,psjki,ppabs,ccci,tau0, &
|
||
alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv)
|
||
call dealloc_pec
|
||
call dealloc_multipass(iwait,iroff,iop,iow,yynext,yypnext,yw0,ypw0, &
|
||
stnext,stv,p0ray,taus,tau1,etau1,cpls,cpl1,lgcpl1,jphi_beam, &
|
||
pins_beam,currins_beam,dpdv_beam,jcd_beam,psipv,chipv)
|
||
! =========== free memory END ===========
|
||
end subroutine gray_main
|
||
|
||
|
||
|
||
|
||
subroutine vectinit(psjki,ppabs,ccci,tau0,alphaabs0,dids0,ccci0,iiv)
|
||
use const_and_precisions, only : wp_, zero
|
||
implicit none
|
||
! arguments
|
||
real(wp_), dimension(:,:), intent(out) :: psjki,ppabs,ccci
|
||
real(wp_), dimension(:), intent(out) :: tau0,alphaabs0,dids0,ccci0
|
||
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
|
||
|
||
end subroutine vectinit
|
||
|
||
|
||
|
||
subroutine ic_gb(xv0c,anv0c,ak0,wcsi,weta,rcicsi,rcieta,phiw,phir, &
|
||
ywrk0,ypwrk0,xc0,du10,gri,ggri,index_rt)
|
||
! beam tracing initial conditions igrad=1
|
||
! !!!!!! check ray tracing initial conditions igrad=0 !!!!!!
|
||
use const_and_precisions, only : wp_,zero,one,pi,half,two,degree,ui=>im
|
||
use math, only : catand
|
||
use gray_params, only : idst
|
||
use beamdata, only : nray,nrayr,nrayth,rwmax
|
||
implicit none
|
||
! arguments
|
||
integer, intent(in) :: index_rt
|
||
real(wp_), dimension(3), intent(in) :: xv0c,anv0c
|
||
real(wp_), intent(in) :: ak0
|
||
real(wp_), intent(in) :: wcsi,weta,rcicsi,rcieta,phiw,phir
|
||
real(wp_), dimension(6,nray), intent(out) :: ywrk0,ypwrk0
|
||
real(wp_), dimension(3,nray), intent(out) :: gri
|
||
real(wp_), dimension(3,3,nray), intent(out) :: ggri
|
||
real(wp_), dimension(3,nrayth,nrayr), intent(out) :: xc0,du10
|
||
|
||
! local variables
|
||
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
|
||
|
||
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
|
||
|
||
phiwrad = phiw*degree
|
||
phirrad = phir*degree
|
||
csphiw = cos(phiwrad)
|
||
snphiw = sin(phiwrad)
|
||
! csphir = cos(phirrad)
|
||
! snphir = sin(phirrad)
|
||
|
||
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*catand(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 = dble(qi1)
|
||
rci2 = dble(qi2)
|
||
ww1 = -dimag(qi1)
|
||
ww2 = -dimag(qi2)
|
||
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)
|
||
|
||
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 = -dimag(qqxx)
|
||
wwyy = -dimag(qqyy)
|
||
wwxy = -dimag(qqxy)
|
||
rcixx = dble(qqxx)
|
||
rciyy = dble(qqyy)
|
||
rcixy = dble(qqxy)
|
||
|
||
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)
|
||
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 = -dimag(dqqxx)
|
||
dwwyy = -dimag(dqqyy)
|
||
dwwxy = -dimag(dqqxy)
|
||
d2wwxx = -dimag(d2qqxx)
|
||
d2wwyy = -dimag(d2qqyy)
|
||
d2wwxy = -dimag(d2qqxy)
|
||
drcixx = dble(dqqxx)
|
||
drciyy = dble(dqqyy)
|
||
drcixy = dble(dqqxy)
|
||
|
||
if(nrayr > 1) then
|
||
dr = rwmax/dble(nrayr-1)
|
||
else
|
||
dr = one
|
||
end if
|
||
ddfu = two*dr**2/ak0 ! twodr2 = 2*dr**2 = 2*rwmax/dble(nrayr-1)
|
||
do j = 1, nrayr
|
||
uj(j) = dble(j-1)
|
||
end do
|
||
|
||
da=2*pi/dble(nrayth)
|
||
do k=1,nrayth
|
||
alfak = (k-1)*da
|
||
sna(k) = sin(alfak)
|
||
csa(k) = cos(alfak)
|
||
end do
|
||
|
||
! central ray
|
||
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
|
||
|
||
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
|
||
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
|
||
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 = -(half*(rcixx*x0t**2 + rciyy*y0t**2) + rcixy*x0t*y0t)
|
||
|
||
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
|
||
|
||
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
|
||
|
||
select case(idst)
|
||
case(1)
|
||
! integration variable: c*t
|
||
denom = one
|
||
case(2)
|
||
! integration variable: Sr
|
||
denom = an20
|
||
case default ! idst=0
|
||
! integration variable: s
|
||
denom = an0
|
||
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)
|
||
call print_output(0,jk,zero,one,xc0(:,k,j),-one,zero,(/zero,zero,zero/), &
|
||
ak0,zero,zero,(/zero,zero,zero/),zero,zero,zero,zero,zero,zero, &
|
||
0,0,0,index_rt,ddr,ddi,zero,zero,(/zero,zero,zero/)) ! st=0, index_rt=1, B=0, N=0, psin=-1, Xg=0, Yg=0, gradXg=0
|
||
end do
|
||
end subroutine ic_gb
|
||
|
||
|
||
|
||
subroutine rkstep(sox,bres,xgcn,y,yp,dgr,ddgr,igrad)
|
||
! Runge-Kutta integrator
|
||
use const_and_precisions, only : wp_
|
||
! use gray_params, only : igrad
|
||
use beamdata, only : h,hh,h6
|
||
implicit none
|
||
real(wp_), intent(in) :: sox,bres,xgcn
|
||
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
|
||
|
||
real(wp_), dimension(6) :: yy,fk1,fk2,fk3,fk4
|
||
real(wp_) :: gr2
|
||
real(wp_), dimension(3) :: dgr2
|
||
|
||
! if(igrad.eq.1) then
|
||
gr2 = dgr(1)**2 + dgr(2)**2 + dgr(3)**2
|
||
dgr2 = 2*(dgr(1)*ddgr(:,1) + dgr(2)*ddgr(:,2) + dgr(3)*ddgr(:,3))
|
||
! end if
|
||
fk1 = yp
|
||
|
||
yy = y + fk1*hh
|
||
call rhs(sox,bres,xgcn,yy,gr2,dgr2,dgr,ddgr,fk2,igrad)
|
||
yy = y + fk2*hh
|
||
call rhs(sox,bres,xgcn,yy,gr2,dgr2,dgr,ddgr,fk3,igrad)
|
||
yy = y + fk3*h
|
||
call rhs(sox,bres,xgcn,yy,gr2,dgr2,dgr,ddgr,fk4,igrad)
|
||
|
||
y = y + h6*(fk1 + 2*fk2 + 2*fk3 + fk4)
|
||
end subroutine rkstep
|
||
|
||
|
||
|
||
subroutine rhs(sox,bres,xgcn,y,gr2,dgr2,dgr,ddgr,dery,igrad)
|
||
! Compute right-hand side terms of the ray equations (dery)
|
||
! used in R-K integrator
|
||
use const_and_precisions, only : wp_
|
||
implicit none
|
||
! arguments
|
||
real(wp_), dimension(6), intent(in) :: y
|
||
real(wp_), intent(in) :: sox,bres,xgcn,gr2
|
||
real(wp_), dimension(3), intent(in) :: dgr2,dgr
|
||
real(wp_), dimension(3,3), intent(in) :: ddgr
|
||
real(wp_), dimension(6), intent(out) :: dery
|
||
integer, intent(in) :: igrad
|
||
! local variables
|
||
real(wp_) :: psinv,dens,btot,xg,yg,anpl,anpr,ajphi
|
||
real(wp_) :: ddr,ddi,dersdst,derdnm
|
||
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,psinv,dens,btot,bv,derbv,xg,yg,derxg,deryg, &
|
||
ajphi)
|
||
|
||
anv = y(4:6)
|
||
call disp_deriv(anv,sox,xg,yg,derxg,deryg,bv,derbv,gr2,dgr2,dgr,ddgr, &
|
||
dery,anpl,anpr,ddr,ddi,dersdst,derdnm,igrad)
|
||
end subroutine rhs
|
||
|
||
|
||
|
||
subroutine ywppla_upd(xv,anv,dgr,ddgr,sox,bres,xgcn,dery,psinv,dens,btot, &
|
||
bv,xg,yg,derxg,anpl,anpr,ddr,ddi,dersdst,derdnm,error,igrad)
|
||
! Compute right-hand side terms of the ray equations (dery)
|
||
! used after full R-K step and grad(S_I) update
|
||
use errcodes, only : pnpl
|
||
implicit none
|
||
! arguments
|
||
real(wp_), dimension(3), intent(in) :: xv,anv
|
||
real(wp_), dimension(3), intent(in) :: dgr
|
||
real(wp_), dimension(3,3), intent(in) :: ddgr
|
||
real(wp_), intent(in) :: sox,bres,xgcn
|
||
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
|
||
real(wp_), dimension(3), intent(out) :: derxg
|
||
integer, intent(in) :: igrad
|
||
! local variables
|
||
real(wp_) :: gr2,ajphi
|
||
real(wp_), dimension(3) :: dgr2,deryg
|
||
real(wp_), dimension(3,3) :: derbv
|
||
real(wp_), parameter :: anplth1 = 0.99_wp_, anplth2 = 1.05_wp_
|
||
|
||
gr2 = dgr(1)**2 + dgr(2)**2 + dgr(3)**2
|
||
dgr2 = 2*(dgr(1)*ddgr(:,1) + dgr(2)*ddgr(:,2) + dgr(3)*ddgr(:,3))
|
||
call plas_deriv(xv,bres,xgcn,psinv,dens,btot,bv,derbv,xg,yg,derxg,deryg,ajphi)
|
||
call disp_deriv(anv,sox,xg,yg,derxg,deryg,bv,derbv,gr2,dgr2,dgr,ddgr, &
|
||
dery,anpl,anpr,ddr,ddi,dersdst,derdnm,igrad)
|
||
|
||
error=0
|
||
if( abs(anpl) > anplth1) then
|
||
if(abs(anpl) > anplth2) then
|
||
error=ibset(error,pnpl+1)
|
||
else
|
||
error=ibset(error,pnpl)
|
||
end if
|
||
end if
|
||
end subroutine ywppla_upd
|
||
|
||
|
||
|
||
subroutine gradi_upd(ywrk,ak0,xc,du1,gri,ggri)
|
||
use const_and_precisions, only : wp_,zero,half
|
||
use beamdata, only : nray,nrayr,nrayth,twodr2
|
||
implicit none
|
||
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
|
||
|
||
! 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
|
||
|
||
end subroutine gradi_upd
|
||
|
||
|
||
|
||
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)
|
||
use const_and_precisions, only : wp_
|
||
implicit none
|
||
! 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
|
||
use const_and_precisions, only : wp_
|
||
implicit none
|
||
! 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,psinv,dens,btot,bv,derbv, &
|
||
xg,yg,derxg,deryg,ajphi)
|
||
use const_and_precisions, only : wp_,zero,ccj=>mu0inv
|
||
use gray_params, only : iequil
|
||
use equilibrium, only : psia,equinum_fpol,equinum_psi,equian,sgnbphi
|
||
use coreprofiles, only : density
|
||
implicit none
|
||
! arguments
|
||
real(wp_), dimension(3), intent(in) :: xv
|
||
real(wp_), intent(in) :: xgcn,bres
|
||
real(wp_), intent(out) :: psinv,dens,btot,xg,yg
|
||
real(wp_), dimension(3), intent(out) :: bv,derxg,deryg
|
||
real(wp_), dimension(3,3), intent(out) :: derbv
|
||
! local variables
|
||
integer :: jv
|
||
real(wp_) :: xx,yy,zz
|
||
real(wp_) :: b2tot,csphi,drrdx,drrdy,dphidx,dphidy,rr,rr2,rrm,snphi,zzm
|
||
real(wp_), dimension(3) :: dbtot,bvc
|
||
real(wp_), dimension(3,3) :: dbvcdc,dbvdc,dbv
|
||
real(wp_) :: brr,bphi,bzz,ajphi,dxgdpsi
|
||
real(wp_) :: dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,dfpv,ddenspsin
|
||
|
||
xg = zero
|
||
yg = 99._wp_
|
||
psinv = -1._wp_
|
||
dens = zero
|
||
btot = zero
|
||
ajphi = zero
|
||
derxg = zero
|
||
deryg = zero
|
||
bv = zero
|
||
derbv = zero
|
||
|
||
if(iequil==0) return
|
||
|
||
dbtot = zero
|
||
dbv = zero
|
||
dbvcdc = zero
|
||
dbvcdc = zero
|
||
dbvdc = zero
|
||
|
||
xx = xv(1)
|
||
yy = xv(2)
|
||
zz = xv(3)
|
||
|
||
! cylindrical coordinates
|
||
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
|
||
zzm = 1.0e-2_wp_*zz
|
||
rrm = 1.0e-2_wp_*rr
|
||
|
||
if(iequil==1) then
|
||
call equian(rrm,zzm,psinv,fpolv,dfpv,dpsidr,dpsidz, &
|
||
ddpsidrr,ddpsidzz,ddpsidrz)
|
||
else
|
||
call equinum_psi(rrm,zzm,psinv,dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz)
|
||
call equinum_fpol(psinv,fpolv,dfpv)
|
||
end if
|
||
|
||
! compute yg and derivative
|
||
if(psinv < zero) then
|
||
bphi = fpolv/rrm
|
||
btot = abs(bphi)
|
||
yg = btot/bres
|
||
return
|
||
end if
|
||
|
||
! compute xg and derivative
|
||
call density(psinv,dens,ddenspsin)
|
||
xg = xgcn*dens
|
||
dxgdpsi = xgcn*ddenspsin/psia
|
||
|
||
! B = f(psi)/R e_phi+ grad psi x e_phi/R
|
||
bphi = fpolv/rrm
|
||
brr =-dpsidz/rrm
|
||
bzz = dpsidr/rrm
|
||
|
||
! bvc(i) = B_i in cylindrical coordinates
|
||
bvc(1) = brr
|
||
bvc(2) = bphi
|
||
bvc(3) = bzz
|
||
|
||
! bv(i) = B_i in cartesian coordinates
|
||
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/rrm - brr/rrm
|
||
dbvcdc(2,1) = dfpv*dpsidr/rrm - bphi/rrm
|
||
dbvcdc(3,1) = ddpsidrr/rrm - bzz/rrm
|
||
dbvcdc(1,3) = -ddpsidzz/rrm
|
||
dbvcdc(2,3) = dfpv*dpsidz/rrm
|
||
dbvcdc(3,3) = ddpsidrz/rrm
|
||
|
||
! dbvdc(iv,jv) = d Bcart(iv) / dxvcil(jv)
|
||
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)
|
||
dbvdc(3,2) = dbvcdc(3,2)
|
||
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)
|
||
dbv(:,1) = drrdx*dbvdc(:,1) + dphidx*dbvdc(:,2)
|
||
dbv(:,2) = drrdy*dbvdc(:,1) + dphidy*dbvdc(:,2)
|
||
dbv(:,3) = dbvdc(:,3)
|
||
|
||
! B magnitude and derivatives
|
||
b2tot = bv(1)**2 + bv(2)**2 + bv(3)**2
|
||
btot = sqrt(b2tot)
|
||
|
||
dbtot = (bv(1)*dbv(1,:) + bv(2)*dbv(2,:) + bv(3)*dbv(3,:))/btot
|
||
|
||
yg = btot/Bres
|
||
|
||
! convert spatial derivatives from dummy/m -> dummy/cm
|
||
! to be used in rhs
|
||
|
||
! bv(i) = B_i / B ; derbv(i,j) = d (B_i / B) /d x,y,z
|
||
deryg = 1.0e-2_wp_*dbtot/Bres
|
||
bv = bv/btot
|
||
do jv=1,3
|
||
derbv(:,jv) = 1.0e-2_wp_*(dbv(:,jv) - bv(:)*dbtot(jv))/btot
|
||
end do
|
||
|
||
derxg(1) = 1.0e-2_wp_*drrdx*dpsidr*dxgdpsi
|
||
derxg(2) = 1.0e-2_wp_*drrdy*dpsidr*dxgdpsi
|
||
derxg(3) = 1.0e-2_wp_*dpsidz *dxgdpsi
|
||
|
||
! current density computation in Ampere/m^2, ccj==1/mu_0
|
||
ajphi = ccj*(dbvcdc(1,3) - dbvcdc(3,1))
|
||
! ajr=ccj*(dbvcdc(3,2)/rrm-dbvcdc(2,3))
|
||
! ajz=ccj*(bvc(2)/rrm+dbvcdc(2,1)-dbvcdc(1,2))
|
||
|
||
end subroutine plas_deriv
|
||
|
||
|
||
|
||
subroutine disp_deriv(anv,sox,xg,yg,derxg,deryg,bv,derbv,gr2,dgr2,dgr,ddgr, &
|
||
dery,anpl,anpr,ddr,ddi,dersdst,derdnm,igrad)
|
||
use const_and_precisions, only : wp_,zero,one,half,two
|
||
use gray_params, only : idst
|
||
implicit none
|
||
! arguments
|
||
real(wp_), intent(in) :: xg,yg,gr2,sox
|
||
real(wp_), intent(out) :: anpl,anpr,ddr,ddi,derdnm,dersdst
|
||
real(wp_), dimension(3), intent(in) :: anv,bv,derxg,deryg
|
||
real(wp_), dimension(3), intent(in) :: dgr2,dgr
|
||
real(wp_), dimension(3,3), intent(in) :: ddgr,derbv
|
||
real(wp_), dimension(6), intent(out) :: dery
|
||
integer, intent(in) :: igrad
|
||
! local variables
|
||
integer :: iv
|
||
real(wp_) :: yg2,anpl2,anpr2,del,dnl,duh,dan2sdnpl,an2,an2s
|
||
real(wp_) :: dan2sdxg,dan2sdyg,ddelnpl2,ddelnpl2x,ddelnpl2y,denom,derdel
|
||
real(wp_) :: derdom,dfdiadnpl,dfdiadxg,dfdiadyg,fdia,bdotgr !,vgm
|
||
real(wp_), dimension(3) :: derdxv,danpldxv,derdnv,dbgr !,vgv
|
||
|
||
an2 = anv(1)*anv(1) + anv(2)*anv(2) + anv(3)*anv(3)
|
||
anpl = anv(1)*bv(1) + anv(2)*bv(2) + anv(3)*bv(3)
|
||
|
||
anpl2 = anpl**2
|
||
dnl = one - anpl2
|
||
anpr2 = max(an2-anpl2,zero)
|
||
anpr = sqrt(anpr2)
|
||
yg2 = yg**2
|
||
|
||
an2s = one
|
||
dan2sdxg = zero
|
||
dan2sdyg = zero
|
||
dan2sdnpl = zero
|
||
del = zero
|
||
fdia = zero
|
||
dfdiadnpl = zero
|
||
dfdiadxg = zero
|
||
dfdiadyg = zero
|
||
|
||
duh = one - xg - yg2
|
||
if(xg > zero) then
|
||
del = sqrt(dnl**2 + 4.0_wp_*anpl2*(one - xg)/yg2)
|
||
an2s = one - xg - half*xg*yg2*(one + anpl2 + sox*del)/duh
|
||
|
||
dan2sdxg = - half*yg2*(one - yg2)*(one + anpl2 + sox*del)/duh**2 &
|
||
+ sox*xg*anpl2/(del*duh) - one
|
||
dan2sdyg = - xg*yg*(one - xg)*(one + anpl2 + sox*del)/duh**2 &
|
||
+ two*sox*xg*(one - xg)*anpl2/(yg*del*duh)
|
||
dan2sdnpl = - xg*yg2*anpl/duh &
|
||
- sox*xg*anpl*(two*(one - xg) - yg2*dnl)/(del*duh)
|
||
|
||
if(igrad > 0) then
|
||
ddelnpl2 = two*(two*(one - xg)*(one + 3.0_wp_*anpl2**2) &
|
||
- yg2*dnl**3)/yg2/del**3
|
||
fdia = - xg*yg2*(one + half*sox*ddelnpl2)/duh
|
||
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
|
||
dfdiadnpl = 24.0_wp_*sox*xg*(one - xg)*anpl*(one - anpl2**2) &
|
||
/(yg2*del**5)
|
||
dfdiadxg = - yg2*(one - yg2)/duh**2 - sox*yg2*((one - yg2) &
|
||
*ddelnpl2 + xg*duh*ddelnpl2x)/(two*duh**2)
|
||
dfdiadyg = - two*yg*xg*(one - xg)/duh**2 &
|
||
- sox*xg*yg*(two*(one - xg)*ddelnpl2 &
|
||
+ yg*duh*ddelnpl2y)/(two*duh**2)
|
||
end if
|
||
end if
|
||
|
||
bdotgr = bv(1)*dgr(1) + bv(2)*dgr(2) + bv(3)*dgr(3)
|
||
do iv=1,3
|
||
dbgr(iv) = dgr(1)*derbv(1,iv) + bv(1)*ddgr(1,iv) &
|
||
+ dgr(2)*derbv(2,iv) + bv(2)*ddgr(2,iv) &
|
||
+ dgr(3)*derbv(3,iv) + bv(3)*ddgr(3,iv)
|
||
danpldxv(iv) = anv(1)*derbv(1,iv) + anv(2)*derbv(2,iv) + anv(3)*derbv(3,iv)
|
||
end do
|
||
|
||
derdxv = -(derxg*dan2sdxg + deryg*dan2sdyg + danpldxv*dan2sdnpl + &
|
||
igrad*dgr2) &
|
||
+ fdia*bdotgr*dbgr + half*bdotgr**2 &
|
||
*(derxg*dfdiadxg + deryg*dfdiadyg + danpldxv*dfdiadnpl)
|
||
derdnv = two*anv + (half*bdotgr**2*dfdiadnpl - dan2sdnpl)*bv
|
||
|
||
derdnm = sqrt(derdnv(1)**2 + derdnv(2)**2 + derdnv(3)**2)
|
||
|
||
derdom = -two*an2 + two*xg*dan2sdxg + yg*dan2sdyg + anpl*dan2sdnpl &
|
||
+ two*igrad*gr2 - bdotgr**2*(fdia + xg*dfdiadxg &
|
||
+ half*yg*dfdiadyg &
|
||
+ half*anpl*dfdiadnpl)
|
||
|
||
if (idst == 0) then
|
||
! integration variable: s
|
||
denom = derdnm
|
||
else if (idst == 1) then
|
||
! integration variable: c*t
|
||
denom = -derdom
|
||
else
|
||
! integration variable: Sr
|
||
denom = anv(1)*derdnv(1) + anv(2)*derdnv(2) + anv(3)*derdnv(3)
|
||
end if
|
||
|
||
! coefficient for integration in s
|
||
! ds/dst, where st is the integration variable
|
||
dersdst = derdnm/denom
|
||
|
||
! rhs vector
|
||
dery(1:3) = derdnv(:)/denom
|
||
dery(4:6) = -derdxv(:)/denom
|
||
|
||
! vgv : ~ group velocity
|
||
! vgm=0
|
||
! do iv=1,3
|
||
! vgv(iv)=-derdnv(iv)/derdom
|
||
! vgm=vgm+vgv(iv)**2
|
||
! end do
|
||
! vgm=sqrt(vgm)
|
||
|
||
! ddr : dispersion relation (real part)
|
||
! ddi : dispersion relation (imaginary part)
|
||
ddr = an2 - an2s - igrad*(gr2 - half*bdotgr**2*fdia)
|
||
ddi = derdnv(1)*dgr(1) + derdnv(2)*dgr(2) + derdnv(3)*dgr(3)
|
||
|
||
end subroutine disp_deriv
|
||
|
||
|
||
|
||
subroutine alpha_effj(psinv,xg,yg,dens,tekev,ak0,bres,derdnm,anpl,anpr, &
|
||
sox,anprre,anprim,alpha,didp,nhmin,nhmax,iokhawa,error)
|
||
use const_and_precisions, only : wp_,zero,pi,mc2=>mc2_
|
||
use gray_params, only : iwarm,ilarm,ieccd,imx
|
||
use coreprofiles, only : fzeff
|
||
use equilibrium, only : sgnbphi
|
||
use dispersion, only : harmnumber, warmdisp
|
||
use eccd, only : setcdcoeff,eccdeff,fjch0,fjch,fjncl
|
||
use errcodes, only : palph
|
||
use magsurf_data, only : fluxval
|
||
implicit none
|
||
! arguments
|
||
real(wp_),intent(in) ::psinv,ak0,bres
|
||
real(wp_),intent(in) :: xg,yg,tekev,dens,anpl,anpr,derdnm,sox
|
||
real(wp_),intent(out) :: anprre,anprim,alpha,didp
|
||
integer, intent(out) :: nhmin,nhmax,iokhawa
|
||
integer, intent(out) :: error
|
||
! local variables
|
||
real(wp_) :: rbavi,rrii,rhop
|
||
integer :: lrm,ithn,ierrcd
|
||
real(wp_) :: amu,ratiovgr,rbn,rbx
|
||
real(wp_) :: zeff,cst2,bmxi,bmni,fci
|
||
real(wp_), dimension(:), pointer :: eccdpar=>null()
|
||
real(wp_) :: effjcd,effjcdav,akim,btot
|
||
complex(wp_) :: ex,ey,ez
|
||
|
||
alpha=zero
|
||
anprim=zero
|
||
anprre=zero
|
||
didp=zero
|
||
nhmin=0
|
||
nhmax=0
|
||
iokhawa=0
|
||
error=0
|
||
|
||
if(tekev>zero) then
|
||
! absorption computation
|
||
amu=mc2/tekev
|
||
call harmnumber(yg,amu,anpl,nhmin,nhmax,iwarm)
|
||
if(nhmin.gt.0) then
|
||
lrm=max(ilarm,nhmax)
|
||
call warmdisp(xg,yg,amu,anpl,anpr,sox,lrm,error,anprre,anprim, &
|
||
iwarm,imx,ex,ey,ez)
|
||
akim=ak0*anprim
|
||
ratiovgr=2.0_wp_*anpr/derdnm!*vgm
|
||
alpha=2.0_wp_*akim*ratiovgr
|
||
if(alpha<zero) then
|
||
error=ibset(error,palph)
|
||
return
|
||
end if
|
||
|
||
! calcolo della efficienza <j/p>: effjcdav [A m/W ]
|
||
if(ieccd>0) then
|
||
! current drive computation
|
||
zeff=fzeff(psinv)
|
||
ithn=1
|
||
if(lrm>nhmin) ithn=2
|
||
rhop=sqrt(psinv)
|
||
call fluxval(rhop,rri=rrii,rbav=rbavi,bmn=bmni,bmx=bmxi,fc=fci)
|
||
btot=yg*bres
|
||
rbn=btot/bmni
|
||
rbx=btot/bmxi
|
||
|
||
select case(ieccd)
|
||
case(1)
|
||
! cohen model
|
||
call setcdcoeff(zeff,rbn,rbx,cst2,eccdpar)
|
||
call eccdeff(yg,anpl,anprre,dens,amu,ex,ey,ez,nhmin,nhmax, &
|
||
ithn,cst2,fjch,eccdpar,effjcd,iokhawa,ierrcd)
|
||
case(2)
|
||
! no trapping
|
||
call setcdcoeff(zeff,cst2,eccdpar)
|
||
call eccdeff(yg,anpl,anprre,dens,amu,ex,ey,ez,nhmin,nhmax, &
|
||
ithn,cst2,fjch0,eccdpar,effjcd,iokhawa,ierrcd)
|
||
case default
|
||
! neoclassical model
|
||
call setcdcoeff(zeff,rbx,fci,amu,rhop,cst2,eccdpar)
|
||
call eccdeff(yg,anpl,anprre,dens,amu,ex,ey,ez,nhmin,nhmax, &
|
||
ithn,cst2,fjncl,eccdpar,effjcd,iokhawa,ierrcd)
|
||
end select
|
||
error=error+ierrcd
|
||
if(associated(eccdpar)) deallocate(eccdpar)
|
||
effjcdav=rbavi*effjcd
|
||
didp=sgnbphi*effjcdav/(2.0_wp_*pi*rrii)
|
||
end if
|
||
end if
|
||
end if
|
||
end subroutine alpha_effj
|
||
|
||
|
||
|
||
subroutine set_pol(ywrk0,bres,sox,psipol0,chipol0,ext0,eyt0)
|
||
use const_and_precisions, only : wp_,degree,zero,one,half,im
|
||
use beamdata, only : nray,nrayth
|
||
use equilibrium, only : bfield
|
||
use gray_params, only : ipol
|
||
use polarization, only : pol_limit, polellipse, stokes_ce, stokes_ell
|
||
implicit none
|
||
! arguments
|
||
real(wp_), dimension(6,nray), intent(in) :: ywrk0
|
||
real(wp_), intent(in) :: sox,bres
|
||
real(wp_), intent(inout) :: psipol0, chipol0
|
||
complex(wp_), dimension(nray), intent(out) :: ext0, eyt0
|
||
! local variables
|
||
integer :: j,k,jk
|
||
real(wp_), dimension(3) :: xmv, anv, bv
|
||
real(wp_) :: rm, csphi, snphi, bphi, br, bz, qq, uu, vv, deltapol
|
||
|
||
j=1
|
||
k=0
|
||
do jk=1,nray
|
||
k=k+1
|
||
if(jk == 2 .or. k > nrayth) then
|
||
j=j+1
|
||
k=1
|
||
end if
|
||
|
||
if(ipol == 0) then
|
||
xmv=ywrk0(1:3,jk)*0.01_wp_ ! convert from cm to m
|
||
anv=ywrk0(4:6,jk)
|
||
rm=sqrt(xmv(1)**2+xmv(2)**2)
|
||
csphi=xmv(1)/rm
|
||
snphi=xmv(2)/rm
|
||
call bfield(rm,xmv(3),bphi,br,bz)
|
||
! bv(i) = B_i in cartesian coordinates
|
||
bv(1)=br*csphi-bphi*snphi
|
||
bv(2)=br*snphi+bphi*csphi
|
||
bv(3)=bz
|
||
call pol_limit(anv,bv,bres,sox,ext0(jk),eyt0(jk))
|
||
|
||
if (jk == 1) then
|
||
call stokes_ce(ext0(jk),eyt0(jk),qq,uu,vv)
|
||
call polellipse(qq,uu,vv,psipol0,chipol0)
|
||
psipol0=psipol0/degree ! convert from rad to degree
|
||
chipol0=chipol0/degree
|
||
end if
|
||
else
|
||
call stokes_ell(chipol0*degree,psipol0*degree,qq,uu,vv)
|
||
if(qq**2 < one) then
|
||
deltapol=asin(vv/sqrt(one - qq**2))
|
||
ext0(jk)= sqrt(half*(one + qq))
|
||
eyt0(jk)= sqrt(half*(one - qq))*exp(-im*deltapol)
|
||
else
|
||
ext0(jk)= one
|
||
eyt0(jk)= zero
|
||
end if
|
||
endif
|
||
end do
|
||
end subroutine set_pol
|
||
|
||
|
||
subroutine cniteq(rqgrid,zqgrid,matr2dgrid,nr,nz,h,ncon,npts,icount,rcon,zcon)
|
||
use const_and_precisions, only : wp_
|
||
! v2.01 12/07/95 -- written by d v bartlett, jet joint undertaking.
|
||
! (based on an older code)
|
||
use const_and_precisions, only : wp_
|
||
implicit none
|
||
! 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 .gt. 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 .eq. 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(strheader)
|
||
use units, only : uprj0,uwbm,udisp,ucenr,uoutr,upec,usumm
|
||
implicit none
|
||
! subroutine arguments
|
||
character(len=*), dimension(:), intent(in) :: strheader
|
||
! local variables
|
||
integer :: i,l
|
||
|
||
l=size(strheader)
|
||
do i=1,l
|
||
write(uprj0,'(1x,a)') strheader(i)
|
||
write(uprj0+1,'(1x,a)') strheader(i)
|
||
write(uwbm,'(1x,a)') strheader(i)
|
||
write(udisp,'(1x,a)') strheader(i)
|
||
write(ucenr,'(1x,a)') strheader(i)
|
||
write(uoutr,'(1x,a)') strheader(i)
|
||
write(upec,'(1x,a)') strheader(i)
|
||
write(usumm,'(1x,a)') strheader(i)
|
||
end do
|
||
|
||
write(uprj0,'(1x,a)') '#sst j k xt yt zt rt'
|
||
write(uprj0+1,'(1x,a)') '#sst j k xt yt zt rt'
|
||
write(uwbm,'(1x,a)') '#sst w1 w2'
|
||
write(udisp,'(1x,a)') '#sst Dr_Nr Di_Nr'
|
||
write(ucenr,'(1x,a)') '#sst R z phi psin rhot ne Te Btot Bx By Bz Nperp Npl '// &
|
||
'Nx Ny Nz ki alpha tau Pt dIds nhmin nhmax iohkw index_rt ddr Xg Yg dXgdx dXgdy dXgdz'
|
||
write(uoutr,'(1x,a)') '#i k sst x y R z psin tau Npl alpha index_rt'
|
||
write(upec,'(1x,a)') '#rhop rhot Jphi Jcdb dPdV Icdins Pins index_rt'
|
||
write(usumm,'(1x,a)') '#Icd Pa Jphip dPdVp rhotj rhotjava rhotp rhotpav ' // &
|
||
'drhotjava drhotpav ratjamx ratjbmx stmx psipol chipol index_rt ' // &
|
||
'Jphimx dPdVmx drhotj drhotp P0 cplO cplX'
|
||
end subroutine print_headers
|
||
|
||
|
||
|
||
subroutine print_prof
|
||
use const_and_precisions, only : wp_
|
||
use equilibrium, only : psinr,nq,fq,frhotor,tor_curr_psi
|
||
use coreprofiles, only : density, temp
|
||
use units, only : uprfin
|
||
implicit none
|
||
! local constants
|
||
real(wp_), parameter :: eps=1.e-4_wp_
|
||
! local variables
|
||
integer :: i
|
||
real(wp_) :: psin,rhot,ajphi,dens,ddens
|
||
|
||
write(uprfin,*) ' #psi rhot ne Te q Jphi'
|
||
do i=1,nq
|
||
psin=psinr(i)
|
||
rhot=frhotor(sqrt(psin))
|
||
call density(psin,dens,ddens)
|
||
call tor_curr_psi(max(eps,psin),ajphi)
|
||
|
||
write(uprfin,"(12(1x,e12.5))") psin,rhot,dens,temp(psin),fq(psin),ajphi*1.e-6_wp_
|
||
end do
|
||
end subroutine print_prof
|
||
|
||
|
||
|
||
subroutine print_bres(bres)
|
||
use const_and_precisions, only : wp_
|
||
use equilibrium, only : rmnm, rmxm, zmnm, zmxm, bfield, nq
|
||
use units, only : ubres
|
||
implicit none
|
||
! arguments
|
||
real(wp_) :: 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(nq), zv(nq)
|
||
real(wp_), dimension(nq,nq) :: btotal
|
||
|
||
|
||
dr = (rmxm-rmnm)/(nq-1)
|
||
dz = (zmxm-zmnm)/(nq-1)
|
||
do j=1,nq
|
||
rv(j) = 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,nq
|
||
zzk=zv(k)
|
||
do j=1,nq
|
||
rrj=rv(j)
|
||
call bfield(rrj,zzk,bbphi,bbr,bbz)
|
||
btotal(j,k)=sqrt(bbr**2+bbz**2+bbphi**2)
|
||
if(btotal(j,k).ge.btmx) btmx=btotal(j,k)
|
||
if(btotal(j,k).le.btmn) btmn=btotal(j,k)
|
||
enddo
|
||
enddo
|
||
|
||
! compute Btot=Bres/n with n=1,5
|
||
write(ubres,*)'#i Btot R z'
|
||
do n=1,5
|
||
bbb=bres/dble(n)
|
||
if (bbb.ge.btmn.and.bbb.le.btmx) then
|
||
nconts=size(ncpts)
|
||
nctot=size(rrcb)
|
||
call cniteq(rv,zv,btotal,nq,nq,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(bres,xgcn,r0,anpl0)
|
||
use const_and_precisions, only : wp_
|
||
use gray_params, only : iequil
|
||
use equilibrium, only : rmnm, rmxm, zmnm, zmxm, equian, equinum_psi, &
|
||
equinum_fpol, nq
|
||
use coreprofiles, only : density, temp
|
||
use units, only : umaps
|
||
implicit none
|
||
! arguments
|
||
real(wp_), intent(in) :: bres,xgcn,r0,anpl0
|
||
! local variables
|
||
integer :: j,k
|
||
real(wp_) :: dr,dz,zk,rj,bphi,br,bz,btot,psin,ne,dne,te,xg,yg,anpl
|
||
real(wp_), dimension(nq) :: r, z
|
||
|
||
|
||
dr = (rmxm-rmnm)/(nq-1)
|
||
dz = (zmxm-zmnm)/(nq-1)
|
||
do j=1,nq
|
||
r(j) = 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,nq
|
||
rj=r(j)
|
||
anpl=anpl0*r0/rj
|
||
do k=1,nq
|
||
zk=z(k)
|
||
if (iequil < 2) then
|
||
call equian(rj,zk,psinv=psin,fpolv=bphi,dpsidr=bz,dpsidz=br)
|
||
else
|
||
call equinum_psi(rj,zk,psinv=psin,dpsidr=bz,dpsidz=br)
|
||
call equinum_fpol(psin,fpolv=bphi)
|
||
end if
|
||
br = -br/rj
|
||
bphi = bphi/rj
|
||
bz = bz/rj
|
||
btot = sqrt(br**2+bphi**2+bz**2)
|
||
yg = btot/bres
|
||
te = temp(psin)
|
||
call density(psin,ne,dne)
|
||
xg = xgcn*ne
|
||
write(umaps,'(12(x,e12.5))') rj,zk,psin,br,bphi,bz,btot,ne,te,xg,yg,anpl
|
||
enddo
|
||
write(umaps,*)
|
||
enddo
|
||
|
||
end subroutine print_maps
|
||
|
||
|
||
subroutine print_surfq(qval)
|
||
use equilibrium, only : psinr, nq, fq, frhotor, &
|
||
rmaxis, zmaxis, zbsup, zbinf
|
||
use magsurf_data, only : contours_psi, npoints, print_contour
|
||
use utils, only : locate, intlin
|
||
use logger, only : log_info
|
||
|
||
implicit none
|
||
|
||
! subroutine arguments
|
||
real(wp_), dimension(:), intent(in) :: qval
|
||
|
||
! local variables
|
||
integer :: i1,i
|
||
real(wp_) :: rup,zup,rlw,zlw,rhot,psival
|
||
real(wp_), dimension(npoints) :: rcn,zcn
|
||
real(wp_), dimension(nq) :: qpsi
|
||
character(256) :: msg ! for log messages formatting
|
||
|
||
! build q profile on psin grid
|
||
do i=1,nq
|
||
qpsi(i) = fq(psinr(i))
|
||
end do
|
||
|
||
! locate ψ surface for q=qval
|
||
call log_info('constant ψ surfaces for:', &
|
||
mod='gray_core', proc='print_surfq')
|
||
do i=1,size(qval)
|
||
! FIXME: check for non monotonous q profile
|
||
call locate(abs(qpsi),nq,qval(i),i1)
|
||
if (i1>0.and.i1<nq) then
|
||
call intlin(abs(qpsi(i1)),psinr(i1),abs(qpsi(i1+1)),psinr(i1+1), &
|
||
qval(i),psival)
|
||
rup=rmaxis
|
||
rlw=rmaxis
|
||
zup=(zbsup+zmaxis)/2.0_wp_
|
||
zlw=(zmaxis+zbinf)/2.0_wp_
|
||
call contours_psi(psival,rcn,zcn,rup,zup,rlw,zlw)
|
||
call print_contour(psival,rcn,zcn)
|
||
rhot=frhotor(sqrt(psival))
|
||
write (msg, '(4(x,a,"=",g0.3))') &
|
||
'q', qval(i), 'ψ', psival, 'rhop', sqrt(psival), 'rhot', rhot
|
||
call log_info(msg, mod='gray_core', proc='print_surfq')
|
||
end if
|
||
end do
|
||
|
||
end subroutine print_surfq
|
||
|
||
|
||
subroutine print_projxyzt(stv,ywrk,iproj)
|
||
use const_and_precisions, only : wp_, comp_huge, zero, one
|
||
use beamdata, only : nray, nrayr, nrayth, rayi2jk
|
||
use units, only : uprj0,uwbm
|
||
implicit none
|
||
! arguments
|
||
real(wp_), 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_) :: dirm,rtimn,rtimx,csth1,snth1,csps1,snps1,xti,yti,zti,rti
|
||
! common/external functions/variables
|
||
|
||
uprj = uprj0 + iproj
|
||
|
||
xv1 = ywrk(1:3,1)
|
||
dir = ywrk(4:6,1)
|
||
dirm = sqrt(dir(1)**2 + dir(2)**2 + dir(3)**2)
|
||
dir = dir/dirm
|
||
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 = sqrt(xti**2 + yti**2)
|
||
|
||
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)
|
||
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
|
||
implicit none
|
||
! 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
|
||
real(wp_), dimension(3), intent(in) :: derxg
|
||
|
||
! local variables
|
||
real(wp_) :: stm,xxm,yym,zzm,rrm,phideg,rhot,akim,pt,didsn
|
||
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=sqrt(xxm**2 + yym**2)
|
||
|
||
! print central ray trajectory. dIds in A/m/W, ki in m^-1
|
||
if(jk.eq.1) then
|
||
phideg=atan2(yym,xxm)/degree
|
||
if(psinv>=zero .and. psinv<=one) then
|
||
rhot=frhotor(sqrt(psinv))
|
||
else
|
||
rhot=1.0_wp_
|
||
end if
|
||
akim=anprim*ak0*1.0e2_wp_
|
||
pt=qj*exp(-tau)
|
||
didsn=dids*1.0e2_wp_/qj
|
||
|
||
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
|
||
end if
|
||
|
||
! print conservation of dispersion relation
|
||
if(jk==nray) write(udisp,'(30(1x,e16.8e3))') st,ddr,ddi
|
||
|
||
! print outer trajectories
|
||
if(mod(i,istpl0)==0) then
|
||
k = jk - jkray1 + 1
|
||
if(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
|
||
end if
|
||
end if
|
||
end subroutine print_output
|
||
|
||
|
||
|
||
subroutine print_pec(rhop_tab,rhot_tab,jphi,jcd,dpdv,currins,pins,index_rt)
|
||
use const_and_precisions, only : wp_
|
||
use units, only : upec
|
||
implicit none
|
||
! arguments
|
||
real(wp_), dimension(:), intent(in) :: rhop_tab,rhot_tab,jphi,jcd,dpdv, &
|
||
currins,pins
|
||
integer, intent(in) :: index_rt
|
||
! local variables
|
||
integer :: i
|
||
|
||
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)
|
||
use const_and_precisions, only : wp_
|
||
use units, only : usumm
|
||
implicit none
|
||
real(wp_), intent(in) :: pabs,icd,dpdvp,jphip,rhotpav,rhotjava,drhotpav, &
|
||
drhotjava,dpdvmx,jphimx,rhotp,rhotj,drhotp,drhotj,ratjamx,ratjbmx, &
|
||
stmx,psipol,chipol,p0,cpl1,cpl2
|
||
integer, intent(in) :: index_rt
|
||
|
||
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
|
||
! write(usumm,*) ''
|
||
end subroutine print_finals
|
||
|
||
end module gray_core
|