From 4226416c4abc82fed73346fcdadff750e7f09537 Mon Sep 17 00:00:00 2001
From: Lorenzo Figini <lorenzo.figini@istp.cnr.it>
Date: Thu, 12 Nov 2015 17:43:12 +0000
Subject: [PATCH] Optimized tests for call to absorption routine. Fixed out of
 bounds error at 1st step

---
 src/graycore.f90 | 23 +++++++++++++++++++----
 1 file changed, 19 insertions(+), 4 deletions(-)

diff --git a/src/graycore.f90 b/src/graycore.f90
index be38ea8..1076170 100644
--- a/src/graycore.f90
+++ b/src/graycore.f90
@@ -54,6 +54,7 @@ contains
   
     real(wp_) :: sox,ak0,bres,xgcn,xg,yg,zzm,alpha,didp,anpl,anpr,anprim,anprre
     real(wp_) :: chipol,psipol,btot,psinv,dens,tekev,zeff,dersdst,derdnm,st,st0
+    real(wp_) :: tau0,alphaabs0,didst0,ccci0
     real(wp_) :: tau,pow,dpdst,ddr,ddi,taumn,taumx
     real(wp_) :: rhotpav,drhotpav,rhotjava,drhotjava
     real(wp_), dimension(3) :: xv,anv0,anv
@@ -175,13 +176,24 @@ contains
         end if
 
         tekev=zero
+        if(i==1) then
+          tau0=zero
+          alphaabs0=zero
+          didst0=zero
+          ccci0=zero
+        else
+          tau0=tauv(jk,i-1)
+          alphaabs0=alphav(jk,i-1)
+          didst0=didst(jk,i-1)
+          ccci0=ccci(jk,i-1)
+        end if
         zzm = xv(3)*0.01_wp_
         ins_pl = (psinv>=zero .and. psinv<one .and. zzm>=zbinf .and. zzm<=zbsup)
         allout = allout .and. .not.ins_pl
         somein = somein .or. ins_pl
 
 ! compute ECRH&CD
-        if(ierr==0 .and. iwarm>0 .and. ins_pl) then
+        if(ierr==0 .and. iwarm>0 .and. ins_pl .and. tau0<=taucr) then
 !          print*,i,jk,rayi2jk(jk),psinv,zzm,anpl
           tekev=temp(psinv)
           if(tekev>zero) then
@@ -192,15 +204,18 @@ contains
           end if
         else
           alpha=zero
+          didp=zero
           anprim=zero
           anprre=anpr
-          didp=zero
+          nharm=0
+          nhf=0
+          iokhawa=0
         end if
 
         psjki(jk,i) = psinv
 
 ! computation of optical depth tau, dP/ds, P(s), dI/ds, I(s)
-        tau=tauv(jk,i-1)+0.5_wp_*(alpha+alphav(jk,i-1))*dersdst*dst
+        tau=tau0+0.5_wp_*(alpha+alphaabs0)*dersdst*dst
         tauv(jk,i)=tau
         alphav(jk,i)=alpha
         pow=p0jk(jk)*exp(-tau) !*exp(-tau1v(jk))
@@ -208,7 +223,7 @@ contains
 
         dpdst=pow*alpha*dersdst
         didst(jk,i)=didp*dpdst
-        ccci(jk,i)=ccci(jk,i-1)+0.5_wp_*(didst(jk,i)+didst(jk,i-1))*dst
+        ccci(jk,i)=ccci0+0.5_wp_*(didst0+didst(jk,i))*dst
 
         call print_output(i,jk,st,p0jk(jk)/p0,xv,psinv,btot,ak0,anpl,anpr,     &
            anprim,dens,tekev,alphav(jk,i),tauv(jk,i),didst(jk,i),nhf,iokhawa,  &