From 608d63acfe5b28082db4c676b76a9cf70ebc4f18 Mon Sep 17 00:00:00 2001
From: Lorenzo Figini <lorenzo.figini@istp.cnr.it>
Date: Fri, 23 Jan 2015 15:08:41 +0000
Subject: [PATCH] fixed few out of bounds checks. added imx (dispersion) read
 from gray.data

---
 Makefile   |  3 +--
 src/gray.f | 18 +++++++++++++-----
 2 files changed, 14 insertions(+), 7 deletions(-)

diff --git a/Makefile b/Makefile
index 4b4c807..0c94f4a 100644
--- a/Makefile
+++ b/Makefile
@@ -11,8 +11,7 @@ vpath %.f src
 
 # Fortran compiler name and flags
 FC=gfortran
-#FFLAGS=-O3 
-FFLAGS=-O3 -Wall -fcheck=all
+FFLAGS=-Wall -fcheck=all
 
 DIRECTIVES = -DREVISION="'$(shell svnversion src)'"
 
diff --git a/src/gray.f b/src/gray.f
index c59513b..f05ecff 100644
--- a/src/gray.f
+++ b/src/gray.f
@@ -663,6 +663,7 @@ c
       common/warm/iwarm,ilarm
       common/ieccd/ieccd
       common/idst/idst
+      common/imx/imx
 c
       common/filesn/filenmeqq,filenmprf,filenmbm
 c
@@ -758,8 +759,10 @@ c     iwarm=1  :weakly relativistic absorption
 c     iwarm=2  :relativistic absorption, n<1 asymptotic expansion
 c     iwarm=3  :relativistic absorption, numerical integration
 c     ilarm    :order of larmor expansion
+c     imx      :max n of iterations in dispersion, imx<0 uses 1st
+c               iteration in case of failure after |imx| iterations
 
-      read(2,*) iwarm,ilarm
+      read(2,*) iwarm,ilarm,imx
 c      if(iwarm.gt.2) iwarm=2
 c
 c     ieccd 0/1 NO/YES  ECCD calculation ieccd>0 different CD models
@@ -1987,8 +1990,9 @@ c
       common/crhot/crhot
 c
       irt=int((nr-1)*psi+1)
-      if(irt.eq.0) irt=1
-      if(irt.eq.nr) irt=nr-1
+c      if(irt.eq.0) irt=1
+c      if(irt.eq.nr) irt=nr-1
+      irt=min(max(1,irt),nr-1)
       dps=psi-psinr(irt)
       frhotor_eq=spli(crhot,nr,irt,dps)
       return
@@ -4852,7 +4856,6 @@ c
 c
       subroutine dispersion(lrm)
       implicit real*8(a-h,o-z)
-      parameter(imx=20)
       complex*16 cc0,cc2,cc4,rr
       complex*16 anpr2a,anpra
       complex*16 anpr,anpr2,ex,ey,ez,den
@@ -4869,6 +4872,7 @@ c
       common/nprw/anprr,anpri
       common/epolar/ex,ey,ez
       common/amut/amu
+      common/imx/imx
 
       errnpr=1.0d0
       anpr2a=anprf**2
@@ -6681,7 +6685,11 @@ c
           rte1=0.0d0
         end if
         call locatex(yy,nd,ipk,nd,yye,ie2)
-        call intlin(yy(ie2),xx(ie2),yy(ie2+1),xx(ie2+1),yye,rte2)
+        if(ie2.gt.0.and.ie2.lt.nd) then
+          call intlin(yy(ie2),xx(ie2),yy(ie2+1),xx(ie2+1),yye,rte2)
+        else
+          rte2=0.0d0
+        end if
       else
         ipk=2
         xpk=xx(2)