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

Fix unsafe use of merge and missing igrad override

Browse Source
This commit is contained in:
Lorenzo Figini 2023-10-20 15:19:56 +02:00
parent e3656e8bdd
commit 621c725948
2 changed files with 24 additions and 13 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

23
src/dispersion.f90
View File

@ -597,7 +597,11 @@ subroutine dielectric_tensor(X, Y, mu, Npl, model, nlarmor, e330, epsl, error)
! Compute the next factorial values
fact_lpn = fact_lpn * (l + n+1)
fact_lmn = merge(fact_lmn / (l - n), 1, n < l)
if (n<l) then
fact_lmn = fact_lmn / (l - n)
else
fact_lmn = 1
end if
! Functions Q⁺_hl(n) Q⁻_hl(n), eq. 7
if (model == 1) then
@ -609,11 +613,18 @@ subroutine dielectric_tensor(X, Y, mu, Npl, model, nlarmor, e330, epsl, error)
Qp_2l = Fp(n,1) + mu * Npl2 * (Fp(n,2) + Fp(n,0) - 2*Fp(n,1))
else
! Fully relativistic computation
Qp_0l = rr(n,0,l) + merge(im*ri(n,0,l) + rr(-n,0,l), czero, n > 0)
Qm_0l = rr(n,0,l) + merge(im*ri(n,0,l) - rr(-n,0,l), czero, n > 0)
Qp_1l = rr(n,1,l) + merge(im*ri(n,1,l) + rr(-n,1,l), czero, n > 0)
Qm_1l = rr(n,1,l) + merge(im*ri(n,1,l) - rr(-n,1,l), czero, n > 0)
Qp_2l = rr(n,2,l) + merge(im*ri(n,2,l) + rr(-n,2,l), czero, n > 0)
Qp_0l = rr(n,0,l)
Qm_0l = rr(n,0,l)
Qp_1l = rr(n,1,l)
Qm_1l = rr(n,1,l)
Qp_2l = rr(n,2,l)
if (n>0) then
Qp_0l = Qp_0l + im*ri(n,0,l) + rr(-n,0,l)
Qm_0l = Qm_0l + im*ri(n,0,l) - rr(-n,0,l)
Qp_1l = Qp_1l + im*ri(n,1,l) + rr(-n,1,l)
Qm_1l = Qm_1l + im*ri(n,1,l) - rr(-n,1,l)
Qp_2l = Qp_2l + im*ri(n,2,l) + rr(-n,2,l)
end if
end if
! Components of the ε̅^(l) tensors, eq. 11

14
src/gray_params.f90
View File

@ -376,7 +376,7 @@ contains
implicit none
! subroutine arguments
type(gray_parameters), intent(in) :: params
type(gray_parameters), intent(inout) :: params
iequil = params%equilibrium%iequil
iprof = params%profiles%iprof
@ -386,17 +386,17 @@ contains
istpr0 = params%output%istpr
istpl0 = params%output%istpl
if (params%raytracing%nrayr < 5) then
params%raytracing%igrad = 0
call log_warning('nrayr < 5 ⇒ optical case only', &
mod="gray_params", proc="set_globals")
end if
ipol = params%raytracing%ipol
igrad = params%raytracing%igrad
idst = params%raytracing%idst
ipass = params%raytracing%ipass
if (params%raytracing%nrayr < 5) then
igrad = 0
call log_warning('nrayr < 5 ⇒ optical case only', &
mod="gray_params", proc="set_globals")
end if
iwarm = params%ecrh_cd%iwarm
ilarm = params%ecrh_cd%ilarm
imx = params%ecrh_cd%imx