From 621c7259484bba2c1c61c245b7ce17caf687413b Mon Sep 17 00:00:00 2001
From: Lorenzo Figini <lorenzo.figini@istp.cnr.it>
Date: Fri, 20 Oct 2023 15:19:56 +0200
Subject: [PATCH] Fix unsafe use of merge and missing igrad override

---
 src/dispersion.f90  | 23 +++++++++++++++++------
 src/gray_params.f90 | 14 +++++++-------
 2 files changed, 24 insertions(+), 13 deletions(-)

diff --git a/src/dispersion.f90 b/src/dispersion.f90
index a4228ad..a2c0eda 100644
--- a/src/dispersion.f90
+++ b/src/dispersion.f90
@@ -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
diff --git a/src/gray_params.f90 b/src/gray_params.f90
index 0e721ff..416ca85 100644
--- a/src/gray_params.f90
+++ b/src/gray_params.f90
@@ -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