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

fixed few out of bounds checks. added imx (dispersion) read from gray.data

Browse Source
This commit is contained in:
Lorenzo Figini 2015-01-23 15:08:41 +00:00
parent 94212dba99
commit 608d63acfe
2 changed files with 14 additions and 7 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
Makefile
View File

@ -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)'"

18
src/gray.f
View File

@ -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)