diff --git a/src/gray.f b/src/gray.f
index 95d8f50..d4804e3 100644
--- a/src/gray.f
+++ b/src/gray.f
@@ -167,9 +167,10 @@ c
 c
 c
       subroutine after_onestep(i,istop)
+      use reflections, only : inside
       implicit real*8 (a-h,o-z)
       complex*16 ext,eyt,extr,eytr,exin2,eyin2
-      parameter(jmx=31,kmx=36,nmx=8000)
+      parameter(jmx=31,kmx=36,nmx=8000,nbb=5000)
       parameter(taucr=12.0d0,pi=3.14159265358979d0,cvdr=pi/180.0d0)
       dimension psjki(jmx,kmx,nmx),ppabs(jmx,kmx,nmx),ccci(jmx,kmx,nmx)
       dimension iiv(jmx,kmx),tauv(jmx,kmx,nmx),alphav(jmx,kmx,nmx)
@@ -177,8 +178,11 @@ c
       dimension tau1v(jmx,kmx),yyrfl(jmx,kmx,6),y(6),dery(6)
       dimension ext(jmx,kmx,0:3),eyt(jmx,kmx,0:3)
       dimension xv(3),anv(3),xvrfl(3),anvrfl(3),xvvac(3),anw(3),anwcl(3)
-      dimension xvjk(jmx,kmx),anvjk(jmx,kmx)
+      dimension xvjk(3,jmx,kmx),anvjk(3,jmx,kmx)
       dimension ywrk(6,jmx,kmx),ypwrk(6,jmx,kmx)
+      real*8 rlim(nbb),zlim(nbb)
+      logical inside_plasma
+      external inside_plasma
 c
       common/pcjki/ppabs,ccci
       common/atjki/tauv,alphav
@@ -278,7 +282,8 @@ cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
               yyrfl(j,k,4:6)=anv
               ihcd(j,k)=0
             end if
-          else if (mod(iop(j,k),2).eq.1.and..not.inside_plasma(rrm,zzm))
+          else if (mod(iop(j,k),2).eq.1.and.
+     .             .not.inside_plasma(rrm,zzm)) then
             iop(j,k)=iop(j,k)+1
             call pol_limit(ext(j,k,iop(j,k)),eyt(j,k,iop(j,k)))
           end if
@@ -287,7 +292,7 @@ cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
             if (iow(j,k).eq.0 .and. inside(rlim,zlim,nlim,rrm,zzm)) then
               iow(j,k)=1
             else if (iow(j,k).eq.1 .and. 
-     .               .not.inside(rlim,zlim,nlim,rrm,zzm))
+     .               .not.inside(rlim,zlim,nlim,rrm,zzm)) then
               iow(j,k)=2
               if (inside_plasma(rrm,zzm)) then
                 iop(j,k)=iop(j,k)+1
@@ -317,7 +322,6 @@ cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
                 iop(j,k)=iop(j,k)+1
                 call pol_limit(ext(j,k,iop(j,k)),eyt(j,k,iop(j,k)))
                 if (index_rt.eq.1) ihcd(j,k)=0
-                end if
               end if
             end if
           end if
@@ -327,8 +331,8 @@ cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
           if(iow(j,k).lt.iowmin) iowmin=iow(j,k)
           if(iow(j,k).gt.iowmax) iowmax=iow(j,k)
 
-          xvjk(j,k)=xv
-          anvjk(j,k)=anv
+          xvjk(:,j,k)=xv
+          anvjk(:,j,k)=anv
 
         end do
       end do
@@ -339,7 +343,7 @@ cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
           kkk=nrayth
           if(j.eq.1) kkk=1
           do k=1,kkk
-            akdotn=dot_product(anvjk(j,k),anwcl)
+            akdotn=dot_product(anvjk(:,j,k),anwcl)
             if(akdotn.lt.aknmin) aknmin=akdotn
           end do
         end do
@@ -383,7 +387,7 @@ c         central ray re-enters the plasma
         istop=1
       else if(ipass.gt.1 .and. index_rt.eq.1 .and.
      .   ((iowmin.gt.1 .and. aknmin.gt.0) .or.
-     .    (iowmax.le.1 .and. iop(1,1).gt.2)) then
+     .    (iowmax.le.1 .and. iop(1,1).gt.2))) then
 c     flag second pass mode coupling as unset
         powrfl=-1.0d0
         qqout=0.0d0
@@ -396,18 +400,18 @@ c     flag second pass mode coupling as unset
 c     store missing initial conditions for the second pass
             if (istore(j,k).eq.0) then
               istore(j,k)=istore(j,k)+1
-              yyrfl(j,k,1:3)=xvjk(j,k)
-              yyrfl(j,k,4:6)=anvjk(j,k)
+              yyrfl(j,k,1:3)=xvjk(:,j,k)
+              yyrfl(j,k,4:6)=anvjk(:,j,k)
               tau1v(j,k)=tauv(j,k,iiv(j,k))
             end if
 c     determine mode coupling at the plasma boundary
             if (powrfl.lt.0.0d0) then
-              call vacuum_rt(xvjk(j,k),anvjk(j,k),xvvac,ivac)
+              call vacuum_rt(xvjk(:,j,k),anvjk(:,j,k),xvvac,ivac)
 c     look for first ray hitting the plasma, starting from the central
 c     and evaluate polarization
               if (ivac.eq.1) then
-                y(1:3)=xvjk(j,k)
-                y(4:6)=anvjk(j,k)
+                y(1:3)=xvjk(:,j,k)
+                y(4:6)=anvjk(:,j,k)
                 call fwork(y,dery)
                 call pol_limit(exin2,eyin2)
                 call stokes(exin2,eyin2,qqin2,uuin2,vvin2)
@@ -423,7 +427,7 @@ c     if found, use its polarization state to compute mode coupling
                       exit powloop
                     end if
                   end do
-                end do
+                end do powloop
 c     if no ray completed a first pass in the plasma, use central ray
 c     initial polarization (possibly reflected)
                 if (qqout.le.0.0d0) then
@@ -6784,11 +6788,12 @@ c      ell=bb/aa
 
       subroutine wall_refl(xv,anv,ext,eyt,xvrfl,anvrfl,extr,eytr,walln,
      .                     irfl)
+      use reflections, only : inters_linewall,inside
       implicit none
       integer*4 irfl
       real*8 anv(3),anv0(3),xv(3),xvrfl(3)
       real*8 walln(3),anvrfl(3),vv1(3),vv2(3),vv3(3)
-      real*8 pi,smax
+      real*8 pi,smax,rrm,zzm
       complex*16 ui,extr,eytr,eztr,ext,eyt
       complex*16 evin(3),evrfl(3) 
       parameter(ui=(0.0d0,1.0d0),pi=3.14159265358979d0)
@@ -6873,7 +6878,8 @@ c     wave vector and electric field after reflection in lab frame
       real*8 anv(3),xvstart(3),xvend(3),walln(3),y(6),dery(6)
       real*8 xv0(3),anv0(3)
       real*8 rlim(nbb),zlim(nbb)
-      logical plfound
+      logical plfound,inside_plasma
+      external inside_plasma
       
       common/limiter/rlim,zlim,nlim
       common/dsds/dst