diff --git a/src/gray.f b/src/gray.f
index f05ecff..4d73c90 100644
--- a/src/gray.f
+++ b/src/gray.f
@@ -69,7 +69,12 @@ c     second pass into plasma
 
 
       subroutine gray_integration
-      implicit real*8 (a-h,o-z)
+      implicit none
+      integer istep,istop,index_rt
+      integer nstep
+      real*8 st,dst,strfl11
+      integer i
+      real*8 st0
 
       common/ss/st
       common/dsds/dst
@@ -137,24 +142,24 @@ c
 c
 c     print final results on screen
 c
-      print*,' '
-      print'(a,f9.4)','final step (s, ct, Sr)  = ',st
+      write(*,*)
+      write(*,'(a,f9.4)') 'final step (s, ct, Sr)  = ',st
       if(iwarm.gt.0) then
-      print '(a,2e12.5)','taumn, taumx  = ', taumn,taumx
+        write(*,'(a,2e12.5)') 'taumn, taumx  = ', taumn,taumx
       else
-      print '(a,2f9.4)','taumn, taumx  = ', zero,zero
+        write(*,'(a,2f9.4)') 'taumn, taumx  = ', zero,zero
       end if      
 c
-      print'(a,f9.4)','Pabs_tot (MW) = ',pabstot
+      write(*,'(a,f9.4)') 'Pabs_tot (MW) = ',pabstot
       currtka =currtot*1.0d3
-      print'(a,f9.4)','I_tot (kA)    = ',currtka
+      write(*,'(a,f9.4)') 'I_tot (kA)    = ',currtka
 c
       if (index_rt.eq.1) then
-      if(iequil.eq.2) write(6,*) 'EQUILIBRIUM CASE : ',filenmeqq
-      if(iequil.eq.1) write(6,*) 'ANALTYCAL EQUILIBRIUM'
-      if(iprof.eq.1) write(6,*)  'PROFILE file : ',filenmprf
-      if(iprof.eq.0) write(6,*)  'ANALTYCAL PROFILES ne0,Te0',dens0,te0
-      if(ibeam.ge.1) write(6,*)  'LAUNCHER CASE : ',filenmbm
+      if(iequil.eq.2) write(*,*) 'EQUILIBRIUM CASE : ',filenmeqq
+      if(iequil.eq.1) write(*,*) 'ANALTYCAL EQUILIBRIUM'
+      if(iprof.eq.1) write(*,*)  'PROFILE file : ',filenmprf
+      if(iprof.eq.0) write(*,*)  'ANALTYCAL PROFILES ne0,Te0',dens0,te0
+      if(ibeam.ge.1) write(*,*)  'LAUNCHER CASE : ',filenmbm
 
       end if
 c
@@ -248,7 +253,7 @@ c
           call gwork(j,k)
 c 
           if(ierr.gt.0) then
-            print*,'    IERR = ', ierr
+            write(*,*) '    IERR = ', ierr
             if(ierr.eq.97) then
 c              igrad=0
 c              ierr=0 
@@ -959,9 +964,9 @@ c     set simple limiter as two cylindrical walls at rwallm and r00
         rlim(3)=rlim(2)
         rlim(4)=rlim(1)
         rlim(5)=rlim(1)
-        zlim(1)=-dst*nmx*1.d-2
+        zlim(1)=(z00-dst*nmx)*1.d-2
         zlim(2)=zlim(1)
-        zlim(3)=dst*nmx*1.d-2
+        zlim(3)=(z00+dst*nmx)*1.d-2
         zlim(4)=zlim(3)
         zlim(5)=zlim(1)
         ipass=abs(ipass)
@@ -1568,8 +1573,8 @@ c     compute B_toroidal on axis
       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)
+      print'(a,f8.4)','BT_centr= ',btrcen
+      print'(a,f8.4)','BT_axis = ',btaxis
 
 c     compute normalized rho_tor from eqdsk q profile
       call rhotor(nr)
@@ -3693,8 +3698,8 @@ c
       call fpbisp(tr,nsr,tz,nsz,ccspl,3,3,
      .     rrs,nrs,zzs,nzs,ffspl,wrk(1),wrk(5),iwrk(1),iwrk(2))
       psinv=(ffspl(1)-psinop)/psiant
-c      if(psinv.lt.0.0d0)
-c     .   print'(a,3e12.4)', '    psin < 0 , R, z  ',psinv,rpsim,zpsim
+      if(psinv.lt.0.0d0)
+     .   print'(a,3e12.4)', '    psin < 0 , R, z  ',psinv,rpsim,zpsim
 c
       nur=1
       nuz=0
@@ -3945,8 +3950,7 @@ c
         fzeff=zeff
       else
         call locate(psrad,npp,ps,k)
-        if(k.eq.0) k=1
-        if(k.eq.npp) k=npp-1
+        k=max(1,min(k,npp-1))
         dps=ps-psrad(k)
         fzeff=spli(cz,npmx,k,dps)
       endif
@@ -6420,7 +6424,8 @@ c
               else
                 if(xxi(i).gt.rtbc) exit
                 if(xxi(i).lt.xxi(i-1)) then
-                  isev(is)=i
+!                  isev(is)=i    !!!!!!!!!! it should be isev(is)=i-1
+                  isev(is)=i-1
                   is=is+1
                   idecr=-1
                 end if
@@ -6446,12 +6451,13 @@ c
               ind2=nd
               iind=1
             end if
-            do ind=ind1,ind2,iind
-              iind=ind
-              if (idecr.eq.-1) iind=ind-1
-              rt1=rtab1(iind)
+            do ind=ind1,ind2,iind     !!!!!!!!!! is it safe?
+!              iind=ind                !!!!!!!!!! iind reused in the loop!
+              indi=ind                !!!!!!!!!! iind reused in the loop!
+              if (idecr.eq.-1) indi=ind-1
+              rt1=rtab1(indi)
               call locatex(xxi,iise,iise0,iise,rt1,itb1)
-              if(itb1.gt.iise0.and.itb1.lt.iise) then
+              if(itb1.ge.iise0.and.itb1.lt.iise) then
                 call intlin(xxi(itb1),ypt(itb1),xxi(itb1+1),
      .                      ypt(itb1+1),rt1,ppa2)
                 call intlin(xxi(itb1),yamp(itb1),xxi(itb1+1),
@@ -6534,12 +6540,12 @@ c     of gaussian profile
       rhpp=frhopol(rhotpav)
       rhpj=frhopol(rhotjava)
       dpdvp=pabs*2.0d0/(sqrt(pi)*drhotpav*fdvdrhot(rhpp))
+      ajphip=currt*2.0d0/(sqrt(pi)*drhotjava*fdadrhot(rhpj))
+      call ratioj(rhpj,ratjamx,ratjbmx,ratjplmx)
+      
       call profwidth(nd,rtab,dpdv,rhotp,rhopp,dpdvmx,
      .                 drhotp,drhopp)
-
       if(ieccd.ne.0) then
-        ajphip=currt*2.0d0/(sqrt(pi)*drhotjava*fdadrhot(rhpj))
-        call ratioj(rhpj,ratjamx,ratjbmx,ratjplmx)
         call profwidth(nd,rtab,ajphiv,rhotjfi,rhopfi,ajmxfi,
      .                 drhotjfi,drhopfi)
         xps=rhopfi