diff --git a/src/dispersion.f90 b/src/dispersion.f90
index f199454..3818d54 100644
--- a/src/dispersion.f90
+++ b/src/dispersion.f90
@@ -110,7 +110,7 @@ subroutine warmdisp(xg,yg,mu,npl,nprf,sox,lrm,err,nprr,npri,fast,imx,ex,ey,ez)
   complex(wp_) :: ex,ey,ez,den
 ! local variables
   integer :: i,j,k,imxx,ilrm
-  real(wp_) :: errnpr,npl2,s,cr,ci,ez2,ey2,ex2,enx2
+  real(wp_) :: errnpr,npl2,s,ez2,ey2,ex2,enx2
   complex(wp_) :: cc0,cc2,cc4,discr,npra2,npra,npr,npr2,e330, &
                   e11,e22,e12,e13,e23,a13,a31,a23,a32,a33
   complex(wp_) :: epsl(3,3,lrm),sepsl(3,3)
diff --git a/src/gray.f b/src/gray.f
index 4f79a44..cc31738 100644
--- a/src/gray.f
+++ b/src/gray.f
@@ -658,15 +658,19 @@ c
       use graydata_par
       use graydata_anequil
       use coreprofiles, only : psdbnd,read_profiles,tene_scal,set_prfspl
-      use interp_eqprof, only : rmxm,zbmin,zbmax,btrcen,rcen
+      use interp_eqprof, only : rmxm,zbmin,zbmax,rcen,psia,rmaxis,
+     . zmaxis,rbbbs,zbbbs,read_eqdsk,change_cocos,eq_scal,set_eqspl,
+     . btrcen,psinr=>psiq,qpsi=>q
       use reflections, only : rlim,zlim,nlim,alloc_lim
       use beamdata, only : nrayr,nrayth,nstep
       implicit none
 c local variables
-      integer :: nfil,iox,ierr,leqq
+      integer :: nfil,iox,ierr,isgniphi,isgnbphi
       real(wp_) :: dummy,bresg,r00,anr0c,anphi0c,fghz,zeff,
-     . w0csi,w0eta,d0csi,d0eta,phi0,zrcsi,zreta
-      real(wp_), dimension(:), allocatable :: psrad,terad,derad,zfc
+     . w0csi,w0eta,d0csi,d0eta,phi0,zrcsi,zreta,ssplf
+      real(wp_), dimension(:), allocatable :: psrad,terad,derad,zfc,
+     . rv,zv,fpol
+      real(wp_), dimension(:,:), allocatable :: psin
       character(len=8) :: wdat
       character(len=10) :: wtim
 c common/external functions/variables
@@ -790,7 +794,7 @@ c
 c
 c     signum of toroidal B and I
 c     icocos index for equilibrium from COCOS - O. Sauter Feb 2012
-      read(2,*) sgnbphi,sgniphi,icocos
+      read(2,*) isgnbphi,isgniphi,icocos
 c
 c     read parameters for analytical density and temperature
 c     profiles if iprof=0
@@ -920,16 +924,27 @@ c     read equilibrium data  from file if iequil=2
 c
       if (iequil.eq.2) then
         neqdsk=99
-        leqq=len_trim(filenmeqq)
-         open(file=filenmeqq(1:leqq)//'.eqdsk',
-     .        status= 'unknown', unit=neqdsk)
+        call read_eqdsk(trim(filenmeqq)//'.eqdsk',rv,zv,psin,psia,
+     .     psinr,fpol,qpsi,rcen,rmaxis,zmaxis,rbbbs,zbbbs,rlim,zlim,
+     .     ipsinorm,0,ipsinorm,neqdsk) !,idesc,ifreefmt,neqdsk)
+        call change_cocos(psia,fpol,qpsi,icocos,3)
+        call eq_scal(psia,fpol,isgniphi,isgnbphi,factb)
+        ssplf=0.01_wp_
+        call set_eqspl(rv,zv,psin,psinr,fpol,sspl,ssplf)
+        nlim=size(rlim)
         call equidata
-        close(neqdsk)
+c
+c     locate  btot=bres
+c
+        call bfield_res(rv,zv,size(rv),size(zv))
+        deallocate(rv,zv,psin,fpol)
         
 c     print density, temperature, safecty factor, toroidal current dens
 c     versus psi, rhop, rhot
         call print_prof
       end if
+      sgnbphi=isgnbphi
+      sgnbphi=isgniphi
 
       if (iequil.eq.1)   call bres_anal
 
@@ -1211,239 +1226,23 @@ c
       subroutine equidata
       use const_and_precisions, only : wp_,pi
       use utils, only : vmaxmini
-      use dierckx, only : curfit,splev,regrid,bispev,coeff_parder
-      use graydata_flags, only : ipsinorm,sspl,ixp,icocos,neqdsk
-      use graydata_par, only : sgnbphi,sgniphi,factb
-      use interp_eqprof, only : nsr,nsz,nsf,nr,nz,
-     . psia,psiant,psinop,btrcen,rcen,btaxis,rmaxis,zmaxis,rmnm,
-     . rmxm,zmnm,zmxm,dr,dz,zbmin,zbmax,fpolas,phitedge,rup,zup,
-     . rlw,zlw,tr,tz,tfp,cc=>cceq,cfp,cc01=>cceq01,cc10=>cceq10,
-     . cc20=>cceq20,cc02=>cceq02,cc11=>cceq11,psinr,qpsi,rv,zv,psin,
-     . fpol,rbbbs,zbbbs,nbbbs,alloc_equilvec,alloc_bnd
-      use reflections, only : rlim,zlim,nlim,alloc_lim
+      use graydata_flags, only : ixp
+      use graydata_par, only : factb
+      use interp_eqprof, only : 
+     . psia,psiant,psinop,btrcen,btaxis,rmaxis,zmaxis,rmnm,
+     . rmxm,zmnm,zmxm,dr,dz,zbmin,zbmax,phitedge,rup,zup,
+     . rlw,zlw,rbbbs,zbbbs,nbbbs,qpsi=>q,nq,equinum_fpol
       implicit none
-c local constants
-      integer, parameter :: kspl=3
 c local variables
-      integer :: idum,nrest,nzest,lwrk,liwrk,lw10,lw01,lw20,lw02,lw11,
-     . lwrkf,i,j,ij,iopt,ier,nur,nuz,info,
-     . iqmn,iqmx,icocosmod,ierr
-      integer, dimension(:), allocatable :: iwrk,iwrkf
-      real(wp_) :: xdum,drnr1,dznz1,rleft,zmid,psiaxis,psiedge,
-     . dpsinr,psia0,fp,xb,xe,ssfp,psinxp,rmop,zmop,psinoptmp,r10,z10,
+      integer :: i,info,iqmn,iqmx
+      real(wp_) :: psinxp,rmop,zmop,psinoptmp,r10,z10,
      . rxp,zxp,psinxptmp,r1,z1,psin1,rbmin,rbmax,q2,q15,qmin,qmax,
      . rhot15,psi15,rhot2,psi2,rhotsx
-      real(wp_), dimension(:), allocatable :: fpoli,fvpsi,
-     . wrkf,wf,wrk
-      character(len=48) :: stringa
 c common/external functions/variables
       real(wp_) :: frhotor
 c
-c     read from file eqdsk
-c     (see http://fusion.gat.com/efit/g_eqdsk.html)
-c
-c     spline interpolation of psi and derivatives
-c
-      if(icocos.gt.0) then 
-        read (neqdsk,'(a48,3i4)') stringa,idum,nr,nz
-      else
-        read (neqdsk,*) nr,nz
-      end if
-c
-      call alloc_equilvec(ierr)
-      if (ierr.ne.0) stop
-
-      nrest=nr+4
-      nzest=nz+4
-      lwrk=4*(nr+nz)+11*(nrest+nzest)+nrest*nz+nzest+54
-      liwrk=nz+nr+nrest+nzest+3
-      lw10=nr*3+nz*4+nr*nz
-      lw01=nr*4+nz*3+nr*nz
-      lw20=nr*2+nz*4+nr*nz
-      lw02=nr*4+nz*2+nr*nz
-      lw11=nr*3+nz*3+nr*nz
-      lwrkf=nr*4+nrest*16
-c
-      if(allocated(fvpsi)) deallocate(fvpsi)
-      if(allocated(fpoli)) deallocate(fpoli)
-      if(allocated(wrkf)) deallocate(wrkf)
-      if(allocated(iwrkf)) deallocate(iwrkf)
-      if(allocated(wf)) deallocate(wf)
-      if(allocated(wrk)) deallocate(wrk)
-      if(allocated(iwrk)) deallocate(iwrk)
-c
-      allocate(fvpsi(nr*nz),wf(nr),wrk(lwrk),iwrk(liwrk),
-     . fpoli(nr),wrkf(lwrkf),iwrkf(nrest))
-c
-      if(ipsinorm.eq.0) then
-        read (neqdsk,2020) drnr1,dznz1,rcen,rleft,zmid
-        read (neqdsk,2020) rmaxis,zmaxis,psiaxis,psiedge,xdum
-        read (neqdsk,2020) xdum,xdum,xdum,xdum,xdum
-        read (neqdsk,2020) xdum,xdum,xdum,xdum,xdum
-        read (neqdsk,2020) (fpol(i),i=1,nr)
-        read (neqdsk,2020) (xdum,i=1,nr)
-        read (neqdsk,2020) (xdum,i=1,nr)
-        read (neqdsk,2020) (xdum,i=1,nr)
-        read (neqdsk,2020) ((psin(i,j),i=1,nr),j=1,nz)
-        read (neqdsk,2020) (qpsi(i),i=1,nr)
-      else
-        read (neqdsk,*) drnr1,dznz1,rcen,rleft,zmid
-        read (neqdsk,*) rmaxis,zmaxis,psiaxis,psiedge,xdum
-        read (neqdsk,*) xdum,xdum,xdum,xdum,xdum
-        read (neqdsk,*) xdum,xdum,xdum,xdum,xdum
-        read (neqdsk,*) (fpol(i),i=1,nr)
-        read (neqdsk,*) (xdum,i=1,nr)
-        read (neqdsk,*) (xdum,i=1,nr)
-        read (neqdsk,*) (xdum,i=1,nr)
-        read (neqdsk,*) ((psin(i,j),i=1,nr),j=1,nz)
-        read (neqdsk,*) (qpsi(i),i=1,nr)
-      end if
-2020  format (5e16.9)
-
-      if(ipsinorm.eq.0) psin=(psin-psiaxis)/(psiedge-psiaxis)
-      psia0=psiedge-psiaxis
-c
-c  compensate for different reference systems
-c
-      icocosmod=mod(icocos,10)
-c
-      if (mod(icocos,2).eq.0.and.icocosmod.gt.0) then
-c  icocos mod 2 = 0: toroidal angle phi CW (opposite to gray convention)
-        do i=1,nr
-          fpol(i)=-fpol(i)
-        end do
-      end if
-c
-      if (icocosmod.eq.1 .or. icocosmod.eq.4 .or. 
-     &    icocosmod.eq.5 .or. icocosmod.eq.8) then
-c  icocos mod 10 = 1,4,5,8: psi increasing with CCW Ip
-c  icocos mod 10 = 2,3,6,7: psi decreasing with CCW Ip
-        psia0=-psia0
-      end if
-c
-c     length in m !!!
-c
-      dr=drnr1/dble(nr-1)
-      dz=dznz1/dble(nz-1)
-      rv(1)=rleft
-      zv(1)=zmid-dznz1/2.0_wp_
-      dpsinr=1.0_wp_/dble(nr-1)
-c
-      do i=1,nr
-        psinr(i)=(i-1)*dpsinr
-        rv(i)=rv(1)+(i-1)*dr
-      end do
-c
-      do j=1,nz
-        zv(j)=zv(1)+(j-1)*dz
-      end do
-c
-      rmnm=rv(1)
-      rmxm=rv(nr)
-      zmnm=zv(1)
-      zmxm=zv(nz)
-c
-c     psi function
-c
-c  icocos=0: adapt psi to force Ip sign, otherwise maintain psi
-      if (icocosmod.ne.0) sgniphi=sign(1.0_wp_,-psia0)
-c  icocos>10: input psi is in Wb
-c  icocos<10: input psi is in Wb/rad (gray convention)
-      if (icocos.ge.10) psia0=psia0/(2.0_wp_*pi)
-c
-      psia=sign(psia0,-sgniphi)*factb
-c
-      do j=1,nz
-        do i=1,nr
-          ij=nz*(i-1)+j
-          fvpsi(ij)=psin(i,j)
-        enddo
-      enddo
-c
-c     spline fitting/interpolation of psin(i,j) and derivatives
-c
-      iopt=0
-      call regrid(iopt,nr,rv,nz,zv,fvpsi,rmnm,rmxm,zmnm,zmxm,
-     .              kspl,kspl,sspl,nrest,nzest,nsr,tr,nsz,tz,cc,fp,
-     .              wrk,lwrk,iwrk,liwrk,ier)
-c     if ier=-1 data are fitted using sspl=0    
-      if(ier.eq.-1) then
-      sspl=0.0_wp_
-      call regrid(iopt,nr,rv,nz,zv,fvpsi,rmnm,rmxm,zmnm,zmxm,
-     .              kspl,kspl,sspl,nrest,nzest,nsr,tr,nsz,tz,cc,fp,
-     .              wrk,lwrk,iwrk,liwrk,ier)
-      end if
-c
-      nur=1
-      nuz=0
-      call coeff_parder(tr,nsr,tz,nsz,cc,kspl,kspl,nur,nuz,
-     . cc10,lw10,ier)
-c
-      nur=0
-      nuz=1
-      call coeff_parder(tr,nsr,tz,nsz,cc,kspl,kspl,nur,nuz,
-     . cc01,lw01,ier)
-c
-      nur=2
-      nuz=0
-      call coeff_parder(tr,nsr,tz,nsz,cc,kspl,kspl,nur,nuz,
-     . cc20,lw20,ier)
-c
-      nur=0
-      nuz=2
-      call coeff_parder(tr,nsr,tz,nsz,cc,kspl,kspl,nur,nuz,
-     . cc02,lw02,ier)
-c
-      nur=1
-      nuz=1
-      call coeff_parder(tr,nsr,tz,nsz,cc,kspl,kspl,nur,nuz,
-     . cc11,lw11,ier)
-c
-c     scaling of f_poloidal
-c
-c  icocos=0: adapt fpol to force Ip sign, otherwise maintain fpol
-      if (icocosmod.ne.0) sgnbphi=sign(1.0_wp_,fpol(nr))
-c
-      do i=1,nr
-        fpol(i)=sign(fpol(i),sgnbphi)*factb
-        wf(i)=1.0_wp_
-      end do
-      wf(nr)=1.0e2_wp_
-c
-c     spline interpolation of fpol(i)
-c
-      iopt=0
-      xb=0.0_wp_
-      xe=1.0_wp_
-      ssfp=0.01_wp_
-      call curfit(iopt,nr,psinr,fpol,wf,xb,xe,kspl,ssfp,nrest,nsf,
-     . tfp,cfp,fp,wrkf,lwrkf,iwrkf,ier)
-c
-      call splev(tfp,nsf,cfp,3,psinr,fpoli,nr,ier)
-      fpolas=fpoli(nr)
-      btrcen=fpolas/rcen
-c
-c     read plasma boundaries from eqdsk file
-c
-      read (neqdsk,*) nbbbs,nlim
-c
-      call alloc_bnd(ierr)
-      if (ierr.ne.0) stop
-      call alloc_lim(ierr)
-      if (ierr.ne.0) stop
-c
-      if(nbbbs.gt.0) then
-        if(ipsinorm.eq.1)
-     .   read(neqdsk,*) (rbbbs(i),zbbbs(i),i=1,nbbbs)
-        if(ipsinorm.eq.0)
-     .   read(neqdsk,2020) (rbbbs(i),zbbbs(i),i=1,nbbbs)
-      end if
-      if(nlim.gt.0) then
-        if(ipsinorm.eq.1)
-     .   read(neqdsk,*) (rlim(i),zlim(i),i=1,nlim)
-        if(ipsinorm.eq.0) 
-     .   read(neqdsk,2020) (rlim(i),zlim(i),i=1,nlim)
-      end if
+      nbbbs=size(rbbbs)
+      nq=size(qpsi)
 c
 c     compute max and min z of last closed surface
 c
@@ -1548,7 +1347,8 @@ c
 
 c     compute B_toroidal on axis
 
-      btaxis=fpol(1)/rmaxis
+      call equinum_fpol(0.0_wp_,btaxis)
+      btaxis=btaxis/rmaxis
       print'(a,f8.4)','factb =   ',factb
       print'(a,f8.4)','BT_centr= ',btrcen
       print'(a,f8.4)','BT_axis = ',btaxis
@@ -1581,7 +1381,7 @@ c     locate  psi surface for q=1.5 and q=2
       zlw=(zmaxis+zbmin)/2.0_wp_
       q2=2.0_wp_
       q15=1.5_wp_
-      call vmaxmini(qpsi,nr,qmin,qmax,iqmn,iqmx)
+      call vmaxmini(abs(qpsi),nq,qmin,qmax,iqmn,iqmx)
       if (q15.gt.qmin.and.q15.lt.qmax) then
         call surfq(q15,psi15)
         rhot15=frhotor(sqrt(psi15))
@@ -1595,11 +1395,6 @@ c     locate  psi surface for q=1.5 and q=2
      .   ,sqrt(psi2),'  rhot_2   = ',rhot2
       end if
 c
-c     locate  btot=bres
-c
-      call bfield_res
-c
-      deallocate(iwrk,iwrkf,fpoli,fvpsi,wrkf,wf,wrk)
       return
       end
 c
@@ -1747,7 +1542,7 @@ c
 c
       subroutine print_prof
       use const_and_precisions, only : wp_
-      use interp_eqprof, only : psinr,qpsi,nr
+      use interp_eqprof, only : psinr=>psiq,qpsi=>q,nq
       use coreprofiles, only : density, temp
       use utils, only : intlin
       implicit none
@@ -1771,7 +1566,7 @@ c
 c
       write(55,111) psin,rhot,dens,te,qq,ajphi*1.0e-6_wp_
 c
-      nst=nr
+      nst=nq
       do i=2,nst
         psin=dble(i-1)/dble(nst-1)
         rhop=sqrt(psin)
@@ -1779,12 +1574,12 @@ c
         call density(psin,dens,ddens)
         te=temp(psin)
 c
-        ips=int((nr-1)*psin+1)
+        ips=int((nq-1)*psin+1)
         if(i.lt.nst) then
           call intlin(psinr(ips),qpsi(ips),psinr(ips+1),qpsi(ips+1),
      .            psin,qq)
         else
-          qq=qpsi(nr)
+          qq=qpsi(nq)
         end if
         rhot=frhotor(rhop)
         call tor_curr_psi(psin,ajphi)
@@ -1829,7 +1624,7 @@ c
       subroutine surfq(qval,psival)
       use const_and_precisions, only : wp_
       use magsurf_data, only : npoints
-      use interp_eqprof, only : nr,psinr,qpsi
+      use interp_eqprof, only : psinr=>psiq,qpsi=>q,nq
       use utils, only : locate, intlin
       implicit none
 c arguments
@@ -1842,8 +1637,8 @@ c
 c
 c     locate  psi surface for q=qval
 c
-      call locate(qpsi,nr,qval,i1)
-      call intlin(qpsi(i1),psinr(i1),qpsi(i1+1),psinr(i1+1),
+      call locate(abs(qpsi),nq,qval,i1)
+      call intlin(abs(qpsi(i1)),psinr(i1),abs(qpsi(i1+1)),psinr(i1+1),
      .  qval,psival)
       ipr=1
       call contours_psi(psival,rcn,zcn,ipr)
@@ -1853,10 +1648,13 @@ c
 c
 c
 c
-      subroutine bfield_res
+      subroutine bfield_res(rv,zv,nr,nz)
       use const_and_precisions, only : wp_
-      use interp_eqprof, only : rv,zv,nr,nz,bfield
+      use interp_eqprof, only : bfield
       implicit none
+c arguments
+      integer, intent(in) :: nr, nz
+      real(wp_), intent(in) :: rv(nr), zv(nz)
 c local constants
       integer, parameter :: icmx=2002
 c local variables
@@ -1909,7 +1707,7 @@ c
 c
       subroutine rhotor(rhotsx)
       use const_and_precisions, only : wp_
-      use interp_eqprof, only : nr,psinr,rhopnr,qpsi,crhot,cq
+      use interp_eqprof, only : nq,psinr=>psiq,qpsi=>q,crhot,cq,rhopq
       use simplespline, only : difcs
       implicit none
 c arguments
@@ -1917,34 +1715,34 @@ c arguments
 c local variables
       integer :: iopt,ier,k
       real(wp_) :: dx,drhot
-      real(wp_), dimension(nr) :: rhotnr
+      real(wp_), dimension(nq) :: rhotnr
 c
       if(allocated(cq)) deallocate(cq)
       if(allocated(crhot)) deallocate(crhot)
-      allocate(cq(nr,4),crhot(nr,4))
+      if(allocated(rhopq)) deallocate(rhopq)
+      allocate(cq(nq,4),crhot(nq,4),rhopq(nq))
 c
 c     normalized toroidal rho : ~ Integral q(psi) dpsi
 c
       iopt=0
-      call difcs(psinr,qpsi,nr,iopt,cq,ier)
+      call difcs(psinr,abs(qpsi),nq,iopt,cq,ier)
 c
       rhotnr(1)=0.0_wp_
-      do k=1,nr-1
+      do k=1,nq-1
         dx=psinr(k+1)-psinr(k)
         drhot=dx*(cq(k,1)+dx*(cq(k,2)/2.0_wp_+dx*(cq(k,3)/3.0_wp_
      .        +dx*cq(k,4)/4.0_wp_)))
         rhotnr(k+1)=rhotnr(k)+drhot
       end do
-      rhotsx=rhotnr(nr)
-      do k=1,nr
-        rhotnr(k)=sqrt(rhotnr(k)/rhotnr(nr))
-      end do
+      rhotsx=rhotnr(nq)
+      rhotnr=sqrt(rhotnr/rhotsx)
+      rhotsx=sign(rhotsx,qpsi(nq))
 c
 c     spline interpolation of rhotor
 c
       iopt=0
-      rhopnr=sqrt(psinr)
-      call difcs(rhopnr,rhotnr,nr,iopt,crhot,ier)
+      rhopq=sqrt(psinr)
+      call difcs(rhopq,rhotnr,nq,iopt,crhot,ier)
       return
       end
 c
@@ -1952,7 +1750,7 @@ c
 c
       function fq_eq(psi)
       use const_and_precisions, only : wp_
-      use interp_eqprof, only : psinr,nr,cq
+      use interp_eqprof, only : psinr=>psiq,nq,cq
       use simplespline, only :spli
       implicit none
 c arguments
@@ -1961,11 +1759,11 @@ c local variables
       integer :: irt
       real(wp_) :: dps
 c
-      irt=int((nr-1)*psi+1)
+      irt=int((nq-1)*psi+1)
       if(irt.eq.0) irt=1
-      if(irt.eq.nr) irt=nr-1
+      if(irt.eq.nq) irt=nq-1
       dps=psi-psinr(irt)
-      fq_eq=spli(cq,nr,irt,dps)
+      fq_eq=spli(cq,nq,irt,dps)
 c
       return
       end
@@ -1974,7 +1772,7 @@ c
 c
       function frhotor_eq(rhop)
       use const_and_precisions, only : wp_
-      use interp_eqprof, only : rhopnr,nr,crhot
+      use interp_eqprof, only : rhopq,nq,crhot
       use simplespline, only : spli
       use utils, only : locate
       implicit none
@@ -1984,10 +1782,10 @@ c local variables
       integer :: irt
       real(wp_) :: drh
 c
-      call locate(rhopnr,nr,rhop,irt)
-      irt=min(max(1,irt),nr-1)
-      drh=rhop-rhopnr(irt)
-      frhotor_eq=spli(crhot,nr,irt,drh)
+      call locate(rhopq,nq,rhop,irt)
+      irt=min(max(1,irt),nq-1)
+      drh=rhop-rhopq(irt)
+      frhotor_eq=spli(crhot,nq,irt,drh)
 c
       return
       end
@@ -2304,11 +2102,11 @@ c
       use graydata_par, only : rwallm
       use magsurf_data, only : npoints
       use interp_eqprof, only : psiant,psinop,nsr,nsz,
-     . cc=>cceq,tr,tz,rup,zup,rlw,zlw,nr
+     . cc=>cceq,tr,tz,rup,zup,rlw,zlw,kspl
       use dierckx, only : profil,sproota
       implicit none
 c local constants
-      integer, parameter :: mest=4,kspl=3
+      integer, parameter :: mest=4
 c arguments
       integer :: ipr
       real(wp_) :: h
@@ -2317,7 +2115,7 @@ c local variables
       integer :: np,info,ic,ier,ii,iopt,m
       real(wp_) :: ra,rb,za,zb,th,zc,val
       real(wp_), dimension(mest) :: zeroc
-      real(wp_), dimension(nr+4) :: czc
+      real(wp_), dimension(nsr) :: czc
 c
       np=(npoints-1)/2
 c
@@ -4161,15 +3959,12 @@ c local variables
       real(wp_) :: cst,dpsidrp,d2psidrp,dqq,qq,rn,rpm,sgn,snt,rhop,rhot
 c common/external functions/variables
       real(wp_) :: dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,
-     . dfpv,xg,xgcn,dxgdpsi,dens,ddens
+     . dfpv
       real(wp_) :: frhopol
 c
       common/derip1/dpsidr,dpsidz
       common/derip2/ddpsidrr,ddpsidzz,ddpsidrz
       common/fpol/fpolv,dfpv
-      common/xgxg/xg
-      common/dxgdps/dxgdpsi
-      common/xgcn/xgcn
 c
       dpsidrp=0.0_wp_
       d2psidrp=0.0_wp_
@@ -4263,18 +4058,18 @@ c
       subroutine psi_raxis(h,r1,r2)
       use const_and_precisions, only : wp_
       use interp_eqprof, only : psiant,psinop,zmaxis,nsr,
-     . nsz,cc=>cceq,tr,tz,nr
+     . nsz,cc=>cceq,tr,tz,kspl
       use dierckx, only : profil,sproota
       implicit none
 c local constants
-      integer, parameter :: mest=4,kspl=3
+      integer, parameter :: mest=4
 c arguments
       real(wp_) :: h,r1,r2
 c local variables
       integer :: iopt,ier,m
       real(wp_) :: zc,val
       real(wp_), dimension(mest) :: zeroc
-      real(wp_), dimension(nr+4) :: czc
+      real(wp_), dimension(nsr) :: czc
 c
       iopt=1
       zc=zmaxis
diff --git a/src/interp_eqprof.f90 b/src/interp_eqprof.f90
index b8dbd65..e681cae 100644
--- a/src/interp_eqprof.f90
+++ b/src/interp_eqprof.f90
@@ -3,13 +3,9 @@ module interp_eqprof
   implicit none
   
 ! equidata
-  ! === 2D array psi(R,z) ==========================================
-  INTEGER, SAVE :: nr,nz
-  REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: rv,zv
-  REAL(wp_), DIMENSION(:,:), ALLOCATABLE, SAVE :: psin
   ! === 1D array Fpol(psi), q(psi), rhotor(psi) ====================
   !INTEGER, SAVE :: npsieq
-  REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: psinr,rhopnr,fpol,qpsi
+  !REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: rhopnr
   
   ! === 1D array plasma boundary Rbnd_i, Zbnd_i ====================
   INTEGER, SAVE :: nbbbs
@@ -20,11 +16,11 @@ module interp_eqprof
   REAL(wp_), SAVE :: zbmin,zbmax
   REAL(wp_), SAVE :: phitedge
   REAL(wp_), SAVE :: rup,zup,rlw,zlw
-  REAL(wp_), SAVE :: rcen,btrcen ! rcen unused, btrcen used only for Jcd_ASTRA def.
+  REAL(wp_), SAVE :: rcen,btrcen ! rcen used to compute btrcen from fpol, btrcen used only for Jcd_ASTRA def.
 
+  INTEGER, PARAMETER :: kspl=3,ksplp=kspl+1
   ! === 2D spline psi(R,z), normalization and derivatives ==========
-  INTEGER, SAVE :: nrest, nzest, lw10, lw01, lw20, lw02, lw11
-  INTEGER :: nsr,nsz
+  INTEGER, SAVE :: nsr, nsz
   REAL(wp_), SAVE :: psia, psiant, psinop
   REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: tr,tz
   REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: cceq,   cceq01, cceq10, &
@@ -32,43 +28,314 @@ module interp_eqprof
   ! === 1D spline Fpol(psi) ========================================
   ! INTEGER, SAVE :: npsiest
   INTEGER, SAVE :: nsf
-  REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: tfp,cfp
+  REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: tfp, cfp
   REAL(wp_), SAVE :: fpolas
   ! === 1D spline rhot(rhop), rhop(rhot), q(psi) ===================
   ! computed on psinr,rhopnr [,rhotnr] arrays
+  INTEGER, SAVE :: nq
+  REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: rhopq,psiq,q
   REAL(wp_), DIMENSION(:,:), ALLOCATABLE, SAVE :: crhot,cq
 
 contains
 
-  subroutine alloc_equilvec(ier)
+  subroutine read_eqdsk(filenm,rv,zv,psin,psia,psinr,fpol,q,rvac,rax,zax, &
+                        rbnd,zbnd,rlim,zlim,ipsinorm,idesc,ifreefmt,unit)
+    use const_and_precisions, only : one
+    use utils, only : get_free_unit
     implicit none
-    integer, intent(out) :: ier
-
-    if(nr.le.0.or.nz.le.0) then
-      ier = -1
-      return
+! arguments
+    character(len=*), intent(in) :: filenm
+    real(wp_), intent(out) :: psia,rvac,rax,zax
+    real(wp_), dimension(:), allocatable, intent(out) :: rv,zv,psinr,fpol,q
+    real(wp_), dimension(:), allocatable, intent(out) :: rbnd,zbnd,rlim,zlim
+    real(wp_), dimension(:,:), allocatable, intent(out) :: psin
+    integer, optional, intent(in) :: ipsinorm,idesc,ifreefmt,unit
+! local variables
+    integer, parameter :: indef=0,iddef=1,iffdef=0
+    integer :: in,id,iffmt,u,idum,i,j,nr,nz,nbnd,nlim
+    character(len=48) :: string
+    real(wp_) :: dr,dz,dps,rleft,zmid,zleft,xdum,psiedge,psiaxis
+    
+! set default values if optional arguments are absent
+    in=indef; id=iddef; iffmt=iffdef
+    if(present(ipsinorm)) in=ipsinorm
+    if(present(idesc)) id=idesc
+    if(present(ifreefmt)) iffmt=ifreefmt
+    if (present(unit)) then
+      u=unit
+    else
+      u=get_free_unit()
     end if
 
-    nrest=nr+4
-    nzest=nz+4
-    lw10=nr*3+nz*4+nr*nz
-    lw01=nr*4+nz*3+nr*nz
-    lw20=nr*2+nz*4+nr*nz
-    lw02=nr*4+nz*2+nr*nz
-    lw11=nr*3+nz*3+nr*nz
-    
-    call dealloc_equilvec
-    allocate(rv(nr),zv(nz),tr(nrest),tz(nzest),tfp(nrest),cfp(nrest), &
-     cceq(nr*nz),cceq01(lw01),cceq10(lw10),cceq02(lw02), &
-     cceq20(lw20),cceq11(lw11),psin(nr,nz),psinr(nr),rhopnr(nr),fpol(nr), &
-     qpsi(nr), stat=ier)
-    if (ier/=0) call dealloc_equilvec
-  end subroutine alloc_equilvec
+! open G EQDSK file (see http://fusion.gat.com/efit/g_eqdsk.html)
+    open(file=trim(filenm),status='old',unit=u)
 
-  subroutine dealloc_equilvec
+! get size of main arrays and allocate them
+    if (id==1) then
+      read (u,'(a48,3i4)') string,idum,nr,nz
+    else
+      read (u,*) nr,nz
+    end if
+    if (allocated(rv)) deallocate(rv)
+    if (allocated(zv)) deallocate(zv)
+    if (allocated(psin)) deallocate(psin)
+    if (allocated(psinr)) deallocate(psinr)
+    if (allocated(fpol)) deallocate(fpol)
+    if (allocated(q)) deallocate(q)
+    allocate(rv(nr),zv(nz),psin(nr,nz),psinr(nr),fpol(nr),q(nr))
+
+! store 0D data and main arrays
+    if (iffmt==1) then
+      read (u,*) dr,dz,rvac,rleft,zmid
+      read (u,*) rax,zax,psiaxis,psiedge,xdum
+      read (u,*) xdum,xdum,xdum,xdum,xdum
+      read (u,*) xdum,xdum,xdum,xdum,xdum
+      read (u,*) (fpol(i),i=1,nr)
+      read (u,*) (xdum,i=1,nr)
+      read (u,*) (xdum,i=1,nr)
+      read (u,*) (xdum,i=1,nr)
+      read (u,*) ((psin(i,j),i=1,nr),j=1,nz)
+      read (u,*) (q(i),i=1,nr)
+    else
+      read (u,'(5e16.9)') dr,dz,rvac,rleft,zmid
+      read (u,'(5e16.9)') rax,zax,psiaxis,psiedge,xdum
+      read (u,'(5e16.9)') xdum,xdum,xdum,xdum,xdum
+      read (u,'(5e16.9)') xdum,xdum,xdum,xdum,xdum
+      read (u,'(5e16.9)') (fpol(i),i=1,nr)
+      read (u,'(5e16.9)') (xdum,i=1,nr)
+      read (u,'(5e16.9)') (xdum,i=1,nr)
+      read (u,'(5e16.9)') (xdum,i=1,nr)
+      read (u,'(5e16.9)') ((psin(i,j),i=1,nr),j=1,nz)
+      read (u,'(5e16.9)') (q(i),i=1,nr)
+    end if
+
+! get size of boundary and limiter arrays and allocate them
+    read (u,*) nbnd,nlim
+    if (allocated(rbnd)) deallocate(rbnd)
+    if (allocated(zbnd)) deallocate(zbnd)
+    if (allocated(rlim)) deallocate(rlim)
+    if (allocated(zlim)) deallocate(zlim)
+
+! store boundary and limiter data
+    if(nbnd>0) then
+      allocate(rbnd(nbnd),zbnd(nbnd))
+      if (iffmt==1) then
+        read(u,*) (rbnd(i),zbnd(i),i=1,nbnd)
+      else
+        read(u,'(5e16.9)') (rbnd(i),zbnd(i),i=1,nbnd)
+      end if
+    end if
+    if(nlim>0) then
+      allocate(rlim(nlim),zlim(nlim))
+      if (iffmt==1) then
+        read(u,*) (rlim(i),zlim(i),i=1,nlim)
+      else
+        read(u,'(5e16.9)') (rlim(i),zlim(i),i=1,nlim)
+      end if
+    end if
+
+! reading of G EQDSK file completed    
+    close(u)
+    
+! build rv,zv,psinr arrays and normalize psin
+    zleft=zmid-0.5_wp_*dz
+    dr=dr/(nr-1)
+    dz=dz/(nz-1)
+    dps=one/(nr-1)
+    do i=1,nr
+      psinr(i)=(i-1)*dps
+      rv(i)=rleft+(i-1)*dr
+    end do
+    do i=1,nz
+      zv(i)=zleft+(i-1)*dz
+    end do
+    psia=psiedge-psiaxis
+    if(in==0) psin=(psin-psiaxis)/psia
+  end subroutine read_eqdsk
+  
+  subroutine change_cocos(psia,fpol,q,cocosin,cocosout,ierr)
+    use const_and_precisions, only : zero,one,pi
+    implicit none
+! arguments
+    real(wp_), intent(inout) :: psia
+    real(wp_), dimension(:), intent(inout) :: fpol,q
+    integer, intent(in) :: cocosin, cocosout
+!    real(wp_), intent(out) :: isign,bsign
+    integer, intent(out), optional :: ierr
+! local variables
+    real(wp_) :: isign,bsign
+    integer :: exp2pi,exp2piout
+    logical :: phiccw,psiincr,qpos,phiccwout,psiincrout,qposout
+
+    call decode_cocos(cocosin,exp2pi,phiccw,psiincr,qpos)
+    call decode_cocos(cocosout,exp2piout,phiccwout,psiincrout,qposout)
+
+!   check sign consistency
+    isign=sign(one,psia)
+    if (.not.psiincr) isign=-isign
+    bsign=sign(one,fpol(size(fpol)))
+    if (qpos.neqv.isign*bsign*q(size(q))>zero) then
+!     warning: sign inconsistency found among q, Ipla and Bref
+      q=-q
+      if(present(ierr)) ierr=1
+    else
+      if(present(ierr)) ierr=0
+    end if
+
+! convert cocosin to cocosout
+    if (phiccw.neqv.phiccwout) then
+!     opposite direction of toroidal angle phi in cocosin and cocosout
+!      bsign=-bsign
+!      isign=-isign
+      fpol=-fpol
+    end if
+!   q has opposite sign for given sign of Bphi*Ip
+    if (qpos .neqv. qposout) q=-q
+!   psi and Ip signs don't change accordingly
+    if ((phiccw.eqv.phiccwout) .neqv. (psiincr.eqv.psiincrout)) psia=-psia
+!   convert Wb to Wb/rad or viceversa
+    psia=psia*(2.0_wp_*pi)**(exp2piout-exp2pi)
+  end subroutine change_cocos
+
+  subroutine decode_cocos(cocos,exp2pi,phiccw,psiincr,qpos)
+    implicit none
+    integer, intent(in) :: cocos
+    integer, intent(out) :: exp2pi
+    logical, intent(out) :: phiccw,psiincr,qpos
+    integer :: cmod10,cmod4
+
+    cmod10=mod(cocos,10)
+    cmod4=mod(cmod10,4)
+!  cocos>10 psi in Wb, cocos<10 psi in Wb/rad
+    exp2pi=(cocos-cmod10)/10
+!  cocos mod 10 = 1,3,5,7: toroidal angle phi increasing CCW
+    phiccw=(mod(cmod10,2)==1)
+!  cocos mod 10 = 1,2,5,6: psi increasing with positive Ip
+    psiincr=(cmod4==1 .or. cmod4==2)
+!  cocos mod 10 = 1,2,7,8: q positive for positive Bphi*Ip
+    qpos=(cmod10<3 .or. cmod10>6)
+  end subroutine decode_cocos
+
+  subroutine eq_scal(psia,fpol,isign,bsign,factb)
+    use const_and_precisions, only : one
+    implicit none
+    real(wp_), intent(inout) :: psia
+    integer, intent(inout) :: isign,bsign
+    real(wp_), dimension(:), intent(inout) :: fpol
+    real(wp_), intent(in) :: factb
+    
+    ! isign and bsign ignored on input if equal to 0
+    ! used to assign the direction of Bphi and Ipla BEFORE scaling
+    ! cocos=3 assumed: CCW direction is >0
+    ! Bphi and Ipla scaled by the same factor factb to keep q unchanged
+    ! factb<0 reverses the directions of Bphi and Ipla
+    if(isign/=0) psia=sign(psia,real(-isign,wp_))
+    if(bsign/=0 .and. bsign*fpol(size(fpol))<0) fpol=-fpol
+    psia=psia*factb
+    fpol=fpol*factb
+    isign=int(sign(one,-psia))
+    bsign=int(sign(one,fpol(size(fpol))))
+  end subroutine eq_scal
+
+  subroutine set_eqspl(rv,zv,psin,psinr,fpol,sspl,ssfp)
+    use const_and_precisions, only : zero,one
+    use dierckx, only : regrid,coeff_parder,curfit,splev
+    implicit none
+! arguments
+    real(wp_), dimension(:), intent(in) :: rv,zv,psinr,fpol
+    real(wp_), dimension(:,:), intent(in) :: psin
+    real(wp_), intent(inout) :: sspl,ssfp
+! local constants
+    integer, parameter :: iopt=0
+    integer :: liwrk,lwrk,lw10,lw01,lw20,lw02,lw11,lwrkf
+    integer :: nr,nz,nrest,nzest,npsest,nrz,npsi
+    real(wp_) :: fp
+    real(wp_), dimension(1) :: fpoli
+    real(wp_), dimension(:), allocatable :: fvpsi,wf,wrk
+    integer, dimension(:), allocatable :: iwrk
+! local variables
+    integer :: ier
+!
+!   compute array sizes and prepare working space arrays
+    nr=size(rv)
+    nz=size(zv)
+    nrz=nr*nz
+    nrest=nr+ksplp
+    nzest=nz+ksplp
+    lwrk=4+nrest*nz+(nrest+nzest)*(2*kspl+5)+(nr+nz)*ksplp+max(nz,nrest)
+    liwrk=nz+nr+nrest+nzest+kspl
+!
+    npsi=size(psinr)
+    npsest=npsi+ksplp
+    lwrkf=npsi*ksplp+npsest*(7+3*kspl)
+!
+    allocate(wrk(max(lwrk,lwrkf)),iwrk(max(liwrk,npsest)))
+!
+!   spline fitting/interpolation of psin(i,j) and derivatives
+!
+!   length in m !!!
+!
+    rmnm=rv(1)
+    rmxm=rv(nr)
+    zmnm=zv(1)
+    zmxm=zv(nz)
+!   allocate knots and spline coefficients arrays
+    if (allocated(tr)) deallocate(tr)
+    if (allocated(tz)) deallocate(tz)
+    allocate(tr(nrest),tz(nzest),cceq(nrz))
+!   allocate work arrays
+!   reshape 2D psi array to 1D (transposed) array and compute spline coeffs
+    allocate(fvpsi(nrz))
+    fvpsi=reshape(transpose(psin),(/nrz/))
+    call regrid(iopt,nr,rv,nz,zv,fvpsi,rmnm,rmxm,zmnm,zmxm,         &
+                kspl,kspl,sspl,nrest,nzest,nsr,tr,nsz,tz,cceq,fp,   &
+                wrk(1:lwrk),lwrk,iwrk(1:liwrk),liwrk,ier)
+!   if ier=-1 data are re-fitted using sspl=0    
+    if(ier==-1) then
+      sspl=0.0_wp_
+      call regrid(iopt,nr,rv,nz,zv,fvpsi,rmnm,rmxm,zmnm,zmxm,       &
+                  kspl,kspl,sspl,nrest,nzest,nsr,tr,nsz,tz,cceq,fp, &
+                  wrk(1:lwrk),lwrk,iwrk(1:liwrk),liwrk,ier)
+    end if
+    deallocate(fvpsi)
+!   compute spline coefficients for psi partial derivatives
+    lw10 = nr*(ksplp-1) + nz*ksplp     + nrz
+    lw01 = nr*ksplp     + nz*(ksplp-1) + nrz
+    lw20 = nr*(ksplp-2) + nz*ksplp     + nrz
+    lw02 = nr*ksplp     + nz*(ksplp-2) + nrz
+    lw11 = nr*(ksplp-1) + nz*(ksplp-1) + nrz
+    if (allocated(cceq10)) deallocate(cceq10)
+    if (allocated(cceq01)) deallocate(cceq01)
+    if (allocated(cceq20)) deallocate(cceq20)
+    if (allocated(cceq02)) deallocate(cceq02)
+    if (allocated(cceq11)) deallocate(cceq11)
+    allocate(cceq10(lw10),cceq01(lw01),cceq20(lw20),cceq02(lw02),cceq11(lw11))
+    call coeff_parder(tr,nsr,tz,nsz,cceq,kspl,kspl,1,0,cceq10,lw10,ier)
+    call coeff_parder(tr,nsr,tz,nsz,cceq,kspl,kspl,0,1,cceq01,lw01,ier)
+    call coeff_parder(tr,nsr,tz,nsz,cceq,kspl,kspl,2,0,cceq20,lw20,ier)
+    call coeff_parder(tr,nsr,tz,nsz,cceq,kspl,kspl,0,2,cceq02,lw02,ier)
+    call coeff_parder(tr,nsr,tz,nsz,cceq,kspl,kspl,1,1,cceq11,lw11,ier)
+!
+!   spline interpolation of fpol(i)
+!
+!   allocate knots and spline coefficients arrays
+    if (allocated(tfp)) deallocate(tfp)
+    if (allocated(cfp)) deallocate(cfp)
+    allocate(tfp(npsest),cfp(npsest))
+    allocate(wf(npsi))
+    wf(1:npsi-1)=one
+    wf(npsi)=1.0e2_wp_
+    call curfit(iopt,npsi,psinr,fpol,wf,zero,one,kspl,ssfp,npsest,nsf, &
+                tfp,cfp,fp,wrk(1:lwrkf),lwrkf,iwrk(1:npsest),ier)
+    call splev(tfp,nsf,cfp,kspl,psinr(npsi:npsi),fpoli,1,ier)
+    fpolas=fpoli(1)
+    btrcen=fpolas/rcen
+!   free temporary arrays
+    deallocate(wrk,iwrk,wf)
+  end subroutine set_eqspl
+
+  subroutine unset_eqspl
     implicit none
-    if(allocated(rv)) deallocate(rv)
-    if(allocated(zv)) deallocate(zv)
     if(allocated(tr)) deallocate(tr)
     if(allocated(tz)) deallocate(tz)
     if(allocated(tfp)) deallocate(tfp)
@@ -79,28 +346,7 @@ contains
     if(allocated(cceq02)) deallocate(cceq02)
     if(allocated(cceq20)) deallocate(cceq20)
     if(allocated(cceq11)) deallocate(cceq11)
-    if(allocated(psin)) deallocate(psin)
-    if(allocated(psinr)) deallocate(psinr)
-    if(allocated(fpol)) deallocate(fpol)
-    if(allocated(rhopnr)) deallocate(rhopnr)
-    if(allocated(qpsi)) deallocate(qpsi)
-
-  end subroutine dealloc_equilvec
-
-  subroutine alloc_bnd(ier)
-    implicit none
-    integer, intent(out) :: ier
-    
-    if(nbbbs.lt.0) then
-      ier = -1
-      return
-    end if
-    
-    call dealloc_bnd
-    allocate(rbbbs(nbbbs),zbbbs(nbbbs), &
-     stat=ier)
-    if (ier/=0) call dealloc_bnd
-  end subroutine alloc_bnd
+  end subroutine unset_eqspl
 
   subroutine dealloc_bnd
     implicit none
@@ -123,7 +369,6 @@ contains
     integer, dimension(liwrk) :: iwrk
     real(wp_), dimension(1) :: rrs,zzs,ffspl
     real(wp_), dimension(lwrk) :: wrk
-
 !
 !     here lengths are measured in meters
 !
@@ -138,19 +383,19 @@ contains
         psinv=(ffspl(1)-psinop)/psiant
       end if
       if (present(dpsidr)) then
-        call sub_derpsi(rpsim,zpsim,1,0,dpsidr,cceq10,lw10)
+        call sub_derpsi(rpsim,zpsim,1,0,dpsidr,cceq10)
       end if
       if (present(dpsidz)) then
-        call sub_derpsi(rpsim,zpsim,0,1,dpsidz,cceq01,lw01)
+        call sub_derpsi(rpsim,zpsim,0,1,dpsidz,cceq01)
       end if
       if (present(ddpsidrr)) then
-        call sub_derpsi(rpsim,zpsim,2,0,ddpsidrr,cceq20,lw20)
+        call sub_derpsi(rpsim,zpsim,2,0,ddpsidrr,cceq20)
       end if
       if (present(ddpsidzz)) then
-        call sub_derpsi(rpsim,zpsim,0,2,ddpsidzz,cceq02,lw02)
+        call sub_derpsi(rpsim,zpsim,0,2,ddpsidzz,cceq02)
       end if
       if (present(ddpsidrz)) then
-        call sub_derpsi(rpsim,zpsim,1,1,ddpsidrz,cceq11,lw11)
+        call sub_derpsi(rpsim,zpsim,1,1,ddpsidrz,cceq11)
       end if
     else
       if(present(psinv)) psinv=-1.0_wp_
@@ -162,26 +407,24 @@ contains
     end if
   end subroutine equinum_psi
 
-  subroutine sub_derpsi(rpsim,zpsim,nur,nuz,derpsi,cc,lw)
+  subroutine sub_derpsi(rpsim,zpsim,nur,nuz,derpsi,cc)
     use dierckx, only : fpbisp
     implicit none
-! local constants
-    integer, parameter :: liwrk=2,nrs=1,nzs=1
 ! arguments
-    integer :: nur,nuz,lw
-    real(wp_) :: rpsim,zpsim,derpsi
-    real(wp_), dimension(lw) :: cc
+    real(wp_), intent(in) :: rpsim,zpsim
+    integer, intent(in) :: nur,nuz
+    real(wp_), intent(out) :: derpsi
+    real(wp_), dimension(:) :: cc
 ! local variables
-    integer :: iwr,iwz
-    integer, dimension(liwrk) :: iwrk
+    integer, dimension(1) :: iwrkr,iwrkz
     real(wp_), dimension(1) :: rrs,zzs,ffspl
+    real(wp_), dimension(1,ksplp) :: wrkr
+    real(wp_), dimension(1,ksplp) :: wrkz
 
     rrs(1)=rpsim
     zzs(1)=zpsim
-    iwr=1+(nr-nur-4)*(nz-nuz-4)
-    iwz=iwr+4-nur
-    call fpbisp(tr(nur+1),nsr-2*nur,tz(nuz+1),nsz-2*nuz,cc,3-nur,3-nuz, &
-         rrs,nrs,zzs,nzs,ffspl,cc(iwr),cc(iwz),iwrk(1),iwrk(2))
+    call fpbisp(tr(nur+1),nsr-2*nur,tz(nuz+1),nsz-2*nuz,cc,kspl-nur,kspl-nuz, &
+         rrs,1,zzs,1,ffspl,wrkr,wrkz,iwrkr,iwrkz)
     derpsi=ffspl(1)*psia
   end subroutine sub_derpsi
 
@@ -254,5 +497,5 @@ contains
     if (present(br)) br=-br/rpsim
     if (present(bz)) bz= bz/rpsim
   end subroutine bfield
-  
+
 end module interp_eqprof
diff --git a/src/reflections.f90 b/src/reflections.f90
index 567db67..42f5fa6 100644
--- a/src/reflections.f90
+++ b/src/reflections.f90
@@ -35,7 +35,7 @@ subroutine inters_linewall(xv,kv,rw,zw,nw,sint,normw)
   real(wp_), dimension(nw), intent(in) :: rw,zw
   real(wp_), intent(out) :: sint
   real(wp_), dimension(3), intent(out) :: normw
-  integer :: i,j,ni,iint
+  integer :: i,j,ni,iint,nneg
   real(wp_), dimension(2) :: si,ti
   real(wp_) :: drw,dzw,xint,yint,rint,l,kxy
   real(wp_) :: tol
@@ -46,13 +46,21 @@ subroutine inters_linewall(xv,kv,rw,zw,nw,sint,normw)
   do i=1,nw-1
     !search intersections with i-th wall segment
     call linecone_coord(xv,kv,rw(i:i+1),zw(i:i+1),si,ti,ni)
-    do while (ni>0 .and. si(1)<=tol)
-      !remove solutions with s<=0
-      ni = ni-1
-      si(1) = si(2)
-      ti(1) = ti(2)
-    end do
+    !discard solutions with s<=0
+    nneg=0
     do j=1,ni
+      if (si(j)<=tol) then
+        nneg=j
+      else
+        exit
+      end if
+    end do
+!    do while (ni>0 .and. si(1)<=tol)
+!      ni = ni-1
+!      si(1) = si(2)   ???
+!      ti(1) = ti(2)   ???
+!    end do
+    do j=nneg+1,ni
       if ((si(j)<sint .or. iint==0) .and. ti(j)>=zero .and. ti(j)<=one) then
         !check intersection is in r,z range and keep the closest
         sint = si(j)