From b7c4630657316d9787ab392e8498cd57b64d9617 Mon Sep 17 00:00:00 2001
From: Daniela Farina <daniela.farina@istp.cnr.it>
Date: Fri, 13 Jul 2012 16:24:03 +0000
Subject: [PATCH] GRAY.f bugs fixed to avoid NaN when rays do not enter into
 plasma GRAY.f "ad-hoc" fix in subroutine contour_psi, that MUST be checked

---
 src/gray.f | 76 +++++++++++++++++++++++++++++++++---------------------
 1 file changed, 46 insertions(+), 30 deletions(-)

diff --git a/src/gray.f b/src/gray.f
index 8c5124e..2c8a9a6 100644
--- a/src/gray.f
+++ b/src/gray.f
@@ -208,7 +208,8 @@ c
       common/dsds/dst
       common/index_rt/index_rt
       common/ipass/ipass
-c
+      common/rwallm/rwallm
+
       pabstot=0.0d0
       currtot=0.0d0
       taumn=1d+30
@@ -324,6 +325,9 @@ c              print*,' '
       if(psimin.gt.1.0d0.and.rrm11.gt.rcen.and.index_rt.gt.1)
      .    istop=1
 
+      if(rrm11.lt.rwallm.and.ipass.eq.1) istop=1
+      if(iovmin.eq.2.and.ipass.eq.1) istop=1
+      
       if(iovmin.eq.3) istop=1
      
 c     print ray positions for j=nray in local reference system
@@ -1389,31 +1393,29 @@ c
       end if
       print*,' '
 
-c
 c     compute B_toroidal on axis
-c
+
       btaxis=fpol(1)/rmaxis
       btrcen=abs(btrcen)*factb
       print'(a,f8.4)','factb =   ',factb
       print'(a,f8.4)','|BT_centr|= ',btrcen
       print'(a,f8.4)','|BT_axis| = ',abs(btaxis)
-c
+
 c     compute normalized rho_tor from eqdsk q profile
-c
       call rhotor(nr)
 c      phitedge=deltapsi*rhotsx*2*pi
 c      rrtor=sqrt(phitedge/abs(btrcen)/pi)
-c
+
 c     compute flux surface averaged quantities
-c
+
       rup=rmaxis
       rlw=rmaxis
-      zup=(zbmax+zmaxis)/2.0d0
-      zlw=(zmaxis+zbmin)/2.0d0
+      zup=zmaxis+(zbmax-zmaxis)/10.0d0
+      zlw=zmaxis-(zmaxis-zbmin)/10.0d0
       call flux_average
-c
+
 c     locate  psi surface for q=1.5 and q=2
-c
+
       rup=rmaxis
       rlw=rmaxis
       zup=(zbmax+zmaxis)/2.0d0
@@ -1809,7 +1811,7 @@ c
       
       function frhotor_av(psi)
       implicit real*8(a-h,o-z)
-      parameter(nnintp=101)
+      parameter(nnintp=41)
       dimension rpstab(nnintp),crhotq(nnintp,4)
       common/pstab/rpstab
       common/eqnn/nr,nz,npp,nintp
@@ -2055,18 +2057,17 @@ c
         if(ier.gt.0) print*,'    profil =',ier
         val=h+psiaxis0/psia
         call sproota(val,tr,nsr,czc,zeroc,mest,m,ier)
-       
+        call sort (mest, zeroc)
         if (zeroc(1).gt.rwallm) then
-        rcn(ic)=zeroc(1)
-        zcn(ic)=zc
-        rcn(2*np+2-ic)=zeroc(2)
-        zcn(2*np+2-ic)=zc
+          iic=1
         else
-        rcn(ic)=zeroc(2)
-        zcn(ic)=zc
-        rcn(2*np+2-ic)=zeroc(3)
-        zcn(2*np+2-ic)=zc        
+          iic=2
+          if (zeroc(2).lt.rwallm) iic=3
         end if
+        rcn(ic)=zeroc(iic)
+        zcn(ic)=zc
+        rcn(2*np+2-ic)=zeroc(iic+1)
+        zcn(2*np+2-ic)=zc
       end do
       if (ipr.gt.0) then
       do ii=1,2*np+1
@@ -2077,7 +2078,22 @@ c
       end if
       return
 111   format(i6,12(1x,e12.5))
-99    format(2i6,12(1x,e16.8e3))
+99    format(12(1x,e12.5))
+      end
+
+      subroutine sort(n,a)
+      implicit none
+      integer n,i,j
+      double precision a(n),temp
+      do i = 2, n
+        j = i - 1
+        temp = a(i)
+        do while (j.ge.1.and.a(j).gt.temp)
+          a(j+1) = a(j)
+          j = j - 1
+        end do
+        a(j+1) = temp
+      end do
       end
 c
 c
@@ -2086,7 +2102,7 @@ c
       implicit real*8 (a-h,o-z)
       real*8 lam
 
-      parameter(nnintp=101,ncnt=100,ncntt=2*ncnt+1,nlam=101)
+      parameter(nnintp=41,ncnt=100,ncntt=2*ncnt+1,nlam=41)
       parameter(zero=0.0d0,one=1.0d0)
       parameter(pi=3.14159265358979d0,ccj=1.0d+7/(4.0d0*pi))
       parameter(ksp=3,njest=nnintp+ksp+1,nlest=nlam+ksp+1)
@@ -2179,10 +2195,10 @@ c     computation of flux surface averaged quantities
       ffc(1)=fc
       rhot2q=0.0d0
       
-      rup=rmaxis
-      rlw=rmaxis
-      zup=(zbmax+zmaxis)/2.0d0
-      zlw=(zmaxis+zbmin)/2.0d0
+C      rup=rmaxis
+C      rlw=rmaxis
+C      zup=(zbmax+zmaxis)/2.0d0
+C      zlw=(zmaxis+zbmin)/2.0d0
 
       do jp=2,nintp
         height=dble(jp-1)/dble(nintp-1)
@@ -4208,7 +4224,7 @@ c
       subroutine valpsispl(rpsi,voli,dervoli,areai,rrii,rbavi,
      .                     bmxi,bmni,fci,intp)
       implicit real*8 (a-h,o-z)
-      parameter(nnintp=101)
+      parameter(nnintp=41)
       dimension rpstab(nnintp),cbmx(nnintp,4),cbmn(nnintp,4)
       dimension cvol(nnintp,4),crri(nnintp,4),crbav(nnintp,4)
       dimension carea(nnintp,4),cfc(nnintp,4)
@@ -4242,7 +4258,7 @@ c
 c
       subroutine ratioj(rpsi,ratj1i,ratj2i)
       implicit real*8 (a-h,o-z)
-      parameter(nnintp=101)
+      parameter(nnintp=41)
       dimension rpstab(nnintp),cratj1(nnintp,4),cratj2(nnintp,4)
       common/pstab/rpstab
       common/eqnn/nr,nz,npp,nintp
@@ -5550,7 +5566,7 @@ c     gg=F(u)/u with F(u) as in Cohen paper
 
       subroutine vlambda(alam,psi,fv,dfv)
       implicit real*8 (a-h,o-z)
-      parameter (nnintp=101,nlam=101)
+      parameter (nnintp=41,nlam=41)
       parameter(ksp=3,njest=nnintp+ksp+1,nlest=nlam+ksp+1)
       parameter(lwrk=4*(nnintp+nlam)+11*(njest+nlest)+
      .          njest*nnintp+nlest+54)