From 321d8704319df88bc5670aeabdf80f9e1d02fe22 Mon Sep 17 00:00:00 2001
From: Lorenzo Figini <lorenzo.figini@istp.cnr.it>
Date: Mon, 13 Jul 2015 14:50:41 +0000
Subject: [PATCH] started reorganization of interp_eqprof (equilibrium) and
 coreprofiles modules

---
 Makefile              |   8 +-
 src/eccd.f90          |   7 -
 src/gray.f            | 688 +++++++++---------------------------------
 src/graydata_par.f90  |   2 +-
 src/interp_eqprof.f90 | 183 ++++++++---
 5 files changed, 279 insertions(+), 609 deletions(-)

diff --git a/Makefile b/Makefile
index d9a88e5..3ff231f 100644
--- a/Makefile
+++ b/Makefile
@@ -3,8 +3,8 @@ EXE=gray
 
 # Objects list
 MAINOBJ=gray.o
-OTHOBJ=conical.o const_and_precisions.o dierckx.o dispersion.o eccd.o eierf.o \
-    graydata_anequil.o graydata_flags.o graydata_par.o \
+OTHOBJ=conical.o const_and_precisions.o coreprofiles.o dierckx.o dispersion.o \
+     eccd.o eierf.o graydata_anequil.o graydata_flags.o graydata_par.o \
     interp_eqprof.o magsurf_data.o math.o minpack.o numint.o quadpack.o \
     reflections.o simplespline.o utils.o beamdata.o 
     
@@ -25,11 +25,13 @@ $(EXE): $(MAINOBJ) $(OTHOBJ)
 	$(FC) $(FFLAGS) -o $@ $^
 
 # Dependencies on modules
-gray.o: const_and_precisions.o dierckx.o dispersion.o eccd.o \
+gray.o: const_and_precisions.o coreprofiles.o dierckx.o dispersion.o eccd.o \
         graydata_anequil.o graydata_flags.o graydata_par.o \
         interp_eqprof.o magsurf_data.o math.o minpack.o numint.o quadpack.o \
         reflections.o simplespline.o utils.o beamdata.o 
 conical.o: const_and_precisions.o
+coreprofiles.o: const_and_precisions.o dierckx.o graydata_anequil.o \
+                graydata_flags.o simplespline.o utils.o
 dierckx.o: const_and_precisions.o
 dispersion.o: const_and_precisions.o eierf.o math.o quadpack.o
 eccd.o: const_and_precisions.o conical.o magsurf_data.o dierckx.o numint.o
diff --git a/src/eccd.f90 b/src/eccd.f90
index ea4b391..597d61c 100644
--- a/src/eccd.f90
+++ b/src/eccd.f90
@@ -99,11 +99,6 @@ module eccd
 
 contains
 
-  subroutine initeccd(ieccd)
-    implicit none
-    integer, intent(in) :: ieccd
-  end subroutine initeccd
-  
   subroutine setcdcoeff_notrap(zeff,cst2,eccdpar)
     implicit none
     real(wp_), intent(in) :: zeff
@@ -151,7 +146,6 @@ contains
     real(wp_), intent(in) :: zeff,rbx,fc,amu,rhop
     real(wp_), intent(out) :: cst2
     real(wp_), dimension(:), allocatable, intent(out) :: eccdpar
-    real(wp_) :: alams,pa,fp0s
     real(wp_), dimension(nlmt) :: chlm
     integer :: nlm,ierr,npar
 
@@ -424,7 +418,6 @@ contains
     real(wp_) :: upl,fjch
     real(wp_), dimension(npar) :: extrapar
 ! local variables
-    integer :: nhn
     real(wp_) :: anpl,anprre,ygn,zeff,rb,alams,pa,fp0s,          &
       upr2,gam,u2,u,z5,xi,xib,xibi,fu2b,fu2,gu,gg,dgg,alam,fp0,  &
       dfp0,fh,dfhl,eta
diff --git a/src/gray.f b/src/gray.f
index 644f106..14bb662 100644
--- a/src/gray.f
+++ b/src/gray.f
@@ -118,7 +118,7 @@ c
       subroutine after_gray_integration
       use const_and_precisions, only : wp_,zero
       use graydata_flags, only : ibeam,iwarm,iequil,iprof,
-     . nnd,ipec,ieccd,filenmeqq,filenmprf,filenmbm
+     . nnd,ipec,filenmeqq,filenmprf,filenmbm
       use graydata_anequil, only : dens0,te0
       use beamdata, only : nrayr
       implicit none
@@ -128,11 +128,10 @@ c local variables
       
       real(wp_), dimension(nnd) :: dpdv,ajphiv,pins,currins
       real(wp_), dimension(nnd) :: dvol,darea,ratjav,ratjbv,ratjplv
-      real(wp_), dimension(nnd) :: ajcdav,ajcdbv,ajplv
       real(wp_), dimension(nnd) :: rt_in,rhop_tab,rhot_tab
       real(wp_), dimension(0:nnd) :: rtabpsi1
       real(wp_) :: rhotpav,drhotpav
-      real(wp_) :: rhotjav,rhotjava,drhotjav,drhotjava
+      real(wp_) :: rhotjava,drhotjava
 c common/external functions/variables
       integer :: index_rt
       real(wp_) :: st,taumn,taumx,pabstot,currtot
@@ -209,7 +208,7 @@ c
       use graydata_par, only : psipol0,chipol0
       use interp_eqprof, only : zbmin,zbmax,rlim,zlim,nlim
       use beamdata, only : nrayr,nrayth,psjki,ppabs,ccci,iiv,tauv,ihcd,
-     . istore,iop,iow,tau1v,yyrfl,ywrk,ypwrk,ext,eyt
+     . istore,iop,iow,tau1v,yyrfl,ywrk,ext,eyt
       implicit none
 c local constants
       real(wp_), parameter :: taucr=12.0_wp_
@@ -488,9 +487,9 @@ c
       subroutine print_output(i,j,k)
       use const_and_precisions, only : wp_,pi
       use graydata_flags, only : istpl0
-      use graydata_par, only : psdbnd
+      use coreprofiles, only : psdbnd
       use beamdata, only : ywrk,psjki,tauv,alphav,pdjki,
-     . currj,didst,q,tau1v,nrayr,nrayth
+     . currj,didst,q,tau1v,nrayr !,nrayth
       implicit none
 c local constants
       integer, parameter :: ndim=6
@@ -658,6 +657,7 @@ c
       use graydata_flags
       use graydata_par
       use graydata_anequil
+      use coreprofiles, only : psdbnd
       use interp_eqprof, only : rmxm,rlim,zlim,nlim,zbmin,zbmax,
      . btrcen,rcen,alloc_lim
       use beamdata, only : nrayr,nrayth,nstep
@@ -1213,26 +1213,26 @@ c
       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 : nsrt,nszt,nsft,rlim,zlim,nlim,nr,nz,
+      use interp_eqprof, only : nsr,nsz,nsf,rlim,zlim,nlim,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,
-     . psi,rbbbs,zbbbs,nbbbs,alloc_equilvec,alloc_bnd!,rrtor
+     . fpol,rbbbs,zbbbs,nbbbs,alloc_equilvec,alloc_bnd
       implicit none
 c local constants
       integer, parameter :: kspl=3
 c local variables
       integer :: idum,nrest,nzest,lwrk,liwrk,lw10,lw01,lw20,lw02,lw11,
-     . lwrkf,lwrkbsp,liwrkbsp,i,j,ij,iopt,ier,nsr,nsz,nur,nuz,info,
+     . lwrkf,i,j,ij,iopt,ier,nur,nuz,info,
      . iqmn,iqmx,icocosmod,ierr
-      integer, dimension(:), allocatable :: iwrk,iwrkf,iwrkbsp
-      real(wp_) :: xdum,drnr1,dznz1,rleft,zmid,psiaxis,psiedge,current,
+      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,
      . rxp,zxp,psinxptmp,r1,z1,psin1,rbmin,rbmax,q2,q15,qmin,qmax,
      . rhot15,psi15,rhot2,psi2,rhotsx
-      real(wp_), dimension(:), allocatable :: fpoli,fvpsi,ffvpsi,fpol,
-     . wrkf,wf,wrk,wrkbsp,pres,ffprim,pprim
+      real(wp_), dimension(:), allocatable :: fpoli,fvpsi,
+     . wrkf,wf,wrk
       character(len=48) :: stringa
 c common/external functions/variables
       real(wp_) :: frhotor
@@ -1261,52 +1261,45 @@ c
       lw02=nr*4+nz*2+nr*nz
       lw11=nr*3+nz*3+nr*nz
       lwrkf=nr*4+nrest*16
-      lwrkbsp=4*(nr+nz)
-      liwrkbsp=nr+nz
 c
       if(allocated(fvpsi)) deallocate(fvpsi)
-      if(allocated(pres)) deallocate(pres)
-      if(allocated(ffprim)) deallocate(ffprim)
-      if(allocated(pprim)) deallocate(pprim)
-      if(allocated(ffvpsi)) deallocate(ffvpsi)
-      if(allocated(fpol)) deallocate(fpol)
       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)
-      if(allocated(wrkbsp)) deallocate(wrkbsp)
-      if(allocated(iwrkbsp)) deallocate(iwrkbsp)
 c
-      allocate(fvpsi(nr*nz),pres(nr),ffprim(nr),pprim(nr),
-     . ffvpsi(nr*nz),fpol(nr),fpoli(nr),wrkf(lwrkf),iwrkf(nrest),
-     . wf(nr),wrk(lwrk),iwrk(liwrk),wrkbsp(lwrkbsp),iwrkbsp(liwrkbsp))
+      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,btrcen
-        read (neqdsk,2020) current,xdum,xdum,xdum,xdum
+        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) (pres(i),i=1,nr)
-        read (neqdsk,2020) (ffprim(i),i=1,nr)
-        read (neqdsk,2020) (pprim(i),i=1,nr)
-        read (neqdsk,2020) ((psi(i,j),i=1,nr),j=1,nz)
+        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,btrcen
-        read (neqdsk,*) current,xdum,xdum,xdum,xdum
+        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,*) (pres(i),i=1,nr)
-        read (neqdsk,*) (ffprim(i),i=1,nr)
-        read (neqdsk,*) (pprim(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
@@ -1314,8 +1307,6 @@ c
 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)
-        btrcen=-btrcen
-        current=-current
         do i=1,nr
           fpol(i)=-fpol(i)
         end do
@@ -1325,22 +1316,9 @@ c
      &    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
-        psiedge=-psiedge
-        psiaxis=-psiaxis
-        if (ipsinorm.eq.0) then
-          do j=1,nz
-            do i=1,nr
-              psi(i,j)=-psi(i,j)
-            end do
-          end do
-        end if
+        psia0=-psia0
       end if
 c
-c  add check for Ip/psi and B0/Fpol sign consistency?
-c
-      current=sign(current,psiaxis-psiedge)
-      btrcen=sign(btrcen,fpol(nr))
-c
 c     length in m !!!
 c
       dr=drnr1/dble(nr-1)
@@ -1365,31 +1343,16 @@ c
 c
 c     psi function
 c
-      psia0=psiedge-psiaxis
 c  icocos=0: adapt psi to force Ip sign, otherwise maintain psi
       if (icocosmod.ne.0) sgniphi=sign(1.0_wp_,-psia0)
-      current=sign(current,sgniphi)
-c
-      psia=-sgniphi*abs(psia0)*factb
 c  icocos>10: input psi is in Wb
 c  icocos<10: input psi is in Wb/rad (gray convention)
-      if (icocos.ge.10) psia=psia/(2.0_wp_*pi)
+      if (icocos.ge.10) psia0=psia0/(2.0_wp_*pi)
 c
-c      do j=1,nz
-c        do i=1,nr
-c          write(80,2021) rv(i),zv(j),psi(i,j)
-c        enddo
-c        write(80,*) ' '
-c      enddo
+      psia=sign(psia0,-sgniphi)*factb
 c
       do j=1,nz
         do i=1,nr
-          if(ipsinorm.eq.0) then
-            psin(i,j)=(psi(i,j)-psiaxis)/psia0
-            psi(i,j)=psin(i,j)*psia
-          else
-            psi(i,j)=psin(i,j)*psia
-          end if
           ij=nz*(i-1)+j
           fvpsi(ij)=psin(i,j)
         enddo
@@ -1408,23 +1371,6 @@ c     if ier=-1 data are fitted using sspl=0
      .              kspl,kspl,sspl,nrest,nzest,nsr,tr,nsz,tz,cc,fp,
      .              wrk,lwrk,iwrk,liwrk,ier)
       end if
-      nsrt=nsr
-      nszt=nsz
-      if (sspl.gt.0.0_wp_) then
-        call bispev(tr,nsr,tz,nsz,cc,kspl,kspl,rv,nr,zv,nz,ffvpsi,
-     .   wrkbsp,lwrkbsp,iwrkbsp,liwrkbsp,ier)
-c
-        do j=1,nz
-          do i=1,nr
-            ij=nz*(i-1)+j
-            psin(i,j)=ffvpsi(ij)
-            psi(i,j)=psin(i,j)*psia
-c            write(79,2021) rv(i),zv(j),psin(i,j)
-          enddo
-c            write(79,*) ' '
-        enddo
-      end if
-c2021  format(5(1x,e16.9))
 c
       nur=1
       nuz=0
@@ -1455,10 +1401,9 @@ 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))
-      btrcen=sign(btrcen,sgnbphi)
 c
       do i=1,nr
-        fpol(i)=sgnbphi*abs(fpol(i))*factb
+        fpol(i)=sign(fpol(i),sgnbphi)*factb
         wf(i)=1.0_wp_
       end do
       wf(nr)=1.0e2_wp_
@@ -1469,11 +1414,12 @@ c
       xb=0.0_wp_
       xe=1.0_wp_
       ssfp=0.01_wp_
-      call curfit(iopt,nr,psinr,fpol,wf,xb,xe,kspl,ssfp,nrest,nsft,
+      call curfit(iopt,nr,psinr,fpol,wf,xb,xe,kspl,ssfp,nrest,nsf,
      . tfp,cfp,fp,wrkf,lwrkf,iwrkf,ier)
 c
-      call splev(tfp,nsft,cfp,3,psinr,fpoli,nr,ier)
+      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
@@ -1487,9 +1433,6 @@ c
      .   read(neqdsk,*) (rbbbs(i),zbbbs(i),i=1,nbbbs)
         if(ipsinorm.eq.0)
      .   read(neqdsk,2020) (rbbbs(i),zbbbs(i),i=1,nbbbs)
-c          do i=1,nbbbs
-c            write(51,*) rbbbs(i),zbbbs(i)
-c          end do
       end if
       if(nlim.gt.0) then
         if(ipsinorm.eq.1)
@@ -1602,7 +1545,6 @@ c
 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 = ',btaxis
@@ -1653,8 +1595,7 @@ c     locate  btot=bres
 c
       call bfield_res
 c
-      deallocate(iwrk,iwrkf,iwrkbsp,fpoli,fvpsi,ffvpsi,fpol,
-     . wrkf,wf,wrk,wrkbsp,pres,ffprim,pprim)
+      deallocate(iwrk,iwrkf,fpoli,fvpsi,wrkf,wf,wrk)
       return
       end
 c
@@ -1662,6 +1603,7 @@ c
 c
       subroutine points_ox(rz,zz,rf,zf,psinvf,info)
       use const_and_precisions, only : wp_,comp_eps
+      use interp_eqprof, only : equinum_psi
       use minpack, only : hybrj1
       implicit none
 c local constants
@@ -1674,12 +1616,9 @@ c local variables
       real(wp_), dimension(n) :: xvec,fvec
       real(wp_), dimension(lwa) :: wa
       real(wp_), dimension(ldfjac,n) :: fjac
-c common/external functions/variables
       real(wp_) :: psinv
 c
       external fcnox
-c
-      common/psival/psinv
 c
       xvec(1)=rz
       xvec(2)=zz
@@ -1691,7 +1630,7 @@ c
       end if
       rf=xvec(1)
       zf=xvec(2)
-      call equinum_psi(rf,zf)
+      call equinum_psi(rf,zf,psinv)
       psinvf=psinv
       return
 c
@@ -1701,25 +1640,23 @@ c
 c
       subroutine fcnox(n,x,fvec,fjac,ldfjac,iflag)
       use const_and_precisions, only : wp_
-      use interp_eqprof, only: psia
+      use interp_eqprof, only: psia,equinum_psi
       implicit none
 c arguments
       integer :: n,iflag,ldfjac
       real(wp_), dimension(n) :: x,fvec
       real(wp_), dimension(ldfjac,n) :: fjac
-c common/external functions/variables
+c local variables
       real(wp_) :: dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz
-c
-      common/derip1/dpsidr,dpsidz
-      common/derip2/ddpsidrr,ddpsidzz,ddpsidrz
 c
       select case(iflag)
       case(1)
-          call equinum_derpsi(x(1),x(2),iflag)
+          call equinum_psi(x(1),x(2),dpsidr=dpsidr,dpsidz=dpsidz)
           fvec(1) = dpsidr/psia
           fvec(2) = dpsidz/psia    
       case(2)
-          call equinum_derpsi(x(1),x(2),iflag)
+          call equinum_psi(x(1),x(2),ddpsidrr=ddpsidrr,
+     .        ddpsidzz=ddpsidzz,ddpsidrz=ddpsidrz)
           fjac(1,1) = ddpsidrr/psia
           fjac(1,2) = ddpsidrz/psia
           fjac(2,1) = ddpsidrz/psia
@@ -1770,32 +1707,27 @@ c
 c
       subroutine fcntgo(n,x,fvec,fjac,ldfjac,iflag)
       use const_and_precisions, only : wp_
-      use interp_eqprof, only: psia
+      use interp_eqprof, only: psia, equinum_psi
       implicit none
 c arguments
       integer :: n,ldfjac,iflag
       real(wp_), dimension(n) :: x,fvec
       real(wp_), dimension(ldfjac,n) :: fjac
 c internal variables
-      integer :: ii
+      real(wp_) :: psinv,dpsidr,dpsidz,ddpsidrr,ddpsidrz
 c common/external functions/variables
-      real(wp_) :: dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,
-     . psinv,h
+      real(wp_) :: h
 c
-      common/derip1/dpsidr,dpsidz
-      common/derip2/ddpsidrr,ddpsidzz,ddpsidrz
-      common/psival/psinv
       common/cnpsi/h
 c
       select case(iflag)
         case(1)
-          call equinum_psi(x(1),x(2))
-          call equinum_derpsi(x(1),x(2),iflag)
+          call equinum_psi(x(1),x(2),psinv,dpsidr)
           fvec(1) = psinv-h
           fvec(2) = dpsidr/psia
         case(2)
-          ii=iflag+1
-          call equinum_derpsi(x(1),x(2),ii)
+          call equinum_psi(x(1),x(2),dpsidr=dpsidr,dpsidz=dpsidz,
+     .        ddpsidrr=ddpsidrr,ddpsidrz=ddpsidrz)
           fjac(1,1) = dpsidr/psia
           fjac(1,2) = dpsidz/psia
           fjac(2,1) = ddpsidrr/psia
@@ -1812,6 +1744,7 @@ c
       subroutine print_prof
       use const_and_precisions, only : wp_
       use interp_eqprof, only : psinr,qpsi,nr
+      use coreprofiles, only : density, temp
       use utils, only : intlin
       implicit none
 c local constants
@@ -1819,17 +1752,15 @@ c local constants
 c local variables
       integer :: i,ips,nst
       real(wp_) :: psin,rhop,rhot,ajphi,te,qq
-c common/external functions/variables
       real(wp_) :: dens,ddens
-      real(wp_) :: frhotor,temp
-c
-      common/dens/dens,ddens
+c common/external functions/variables
+      real(wp_) :: frhotor
 c
       write(55,*) ' #psi rhot ne Te q Jphi'
       psin=0.0_wp_
       rhop=0.0_wp_
       rhot=0.0_wp_
-      call density(psin)
+      call density(psin,dens,ddens)
       call tor_curr_psi(eps,ajphi)
       te=temp(psin)
       qq=qpsi(1)
@@ -1841,7 +1772,7 @@ c
         psin=dble(i-1)/dble(nst-1)
         rhop=sqrt(psin)
 c
-        call density(psin)
+        call density(psin,dens,ddens)
         te=temp(psin)
 c
         ips=int((nr-1)*psin+1)
@@ -1864,24 +1795,23 @@ c
 c
       subroutine print_prof_an
       use const_and_precisions, only : wp_
+      use coreprofiles, only : density, temp
       implicit none
 c local constants
       integer, parameter :: nst=51
 c local variables
       integer :: i
       real(wp_) :: psin,rhop,rhot,te
-c common/external functions/variables
       real(wp_) :: dens,ddens
-      real(wp_) :: temp,frhotor
-c
-      common/dens/dens,ddens
+c common/external functions/variables
+      real(wp_) :: frhotor
 c
       write(55,*) ' #psi rhot ne Te'
       do i=1,nst
         psin=dble(i-1)/dble(nst-1)
         rhop=sqrt(psin)
         rhot=frhotor(rhop)
-        call density(psin)
+        call density(psin,dens,ddens)
         te=temp(psin)
         write(55,111) psin,rhot,dens,te
       end do
@@ -1974,10 +1904,10 @@ c
       subroutine profdata
       use const_and_precisions, only : wp_,zero
       use graydata_flags, only : iprof,iscal,nprof
-      use graydata_par, only : psdbnd,factb,factt,factn
-      use interp_eqprof, only : nsfd,npp,psnpp,denpp,ddenpp,d2denpp,
-     . psrad,derad,terad,zfc,tfn,cfn,ct,cz,alloc_profvec
-      use simplespline, only : difcsn
+      use graydata_par, only : factb,factt,factn
+      use coreprofiles, only : nsfd,npp,psnpp,denpp,ddenpp,d2denpp,
+     . psdbnd,psrad,derad,terad,zfc,tfn,cfn,ct,cz,alloc_profvec
+      use simplespline, only : difcs
       use dierckx, only : curfit,splev,splder
       implicit none
 c local variables
@@ -1986,8 +1916,8 @@ c local variables
       real(wp_) :: aat,aan,ffact,psrad0,terad0,derad0,zfc0,psradi,
      . teradi,deradi,zfci,xb,xe,sspl,dpsb,fp
       real(wp_) :: xnv,ynv
-      real(wp_), dimension(:), allocatable :: wf,wrkf,wrkfd,densi,
-     . ddensi,d2densi
+      real(wp_), dimension(:), allocatable :: wf,wrkf,wrkfd
+      real(wp_), dimension(1) :: densi,ddensi,d2densi
 c
 c     read plasma profiles from file if iprof>0
 c
@@ -2013,11 +1943,7 @@ c
         if(allocated(iwrkf)) deallocate(iwrkf)
         if(allocated(wf)) deallocate(wf)
         if(allocated(wrkfd)) deallocate(wrkfd)
-        if(allocated(densi)) deallocate(densi)
-        if(allocated(ddensi)) deallocate(ddensi)
-        if(allocated(d2densi)) deallocate(d2densi)
-        allocate(wrkf(lwrkf),iwrkf(npest),wf(npp),densi(npest),
-     .   wrkfd(npest),ddensi(npest),d2densi(npest))
+        allocate(wrkf(lwrkf),iwrkf(npest),wf(npp),wrkfd(npest))
 c
         read(nprof,*) psrad0,terad0,derad0,zfc0
         if(psrad0.ne.0.0_wp_) psrad0=0.0_wp_
@@ -2038,10 +1964,10 @@ c
 c     spline approximation of  temperature and Zeff
 
         iopt=0
-        call difcsn(psrad,terad,npp,npp,iopt,ct,ier)
+        call difcs(psrad,terad,npp,iopt,ct,ier)
 
         iopt=0
-        call difcsn(psrad,zfc,npp,npp,iopt,cz,ier)
+        call difcs(psrad,zfc,npp,iopt,cz,ier)
 
 c     spline approximation of density
 
@@ -2054,17 +1980,16 @@ c     spline approximation of density
         call curfit(iopt,npp,psrad,derad,wf,xb,xe,kspl,sspl,npest,nsfd,
      .   tfn,cfn,fp,wrkf,lwrkf,iwrkf,ier)
 
-        call splev(tfn,nsfd,cfn,3,psrad,densi,npp,ier)
+        call splev(tfn,nsfd,cfn,3,psrad(npp),densi(1),1,ier)
         nu=1
-        call splder(tfn,nsfd,cfn,3,nu,psrad,ddensi,npp,wrkfd,ier)
-        
+        call splder(tfn,nsfd,cfn,3,nu,psrad(npp),ddensi(1),1,wrkfd,ier)
         nu=2
-        call splder(tfn,nsfd,cfn,3,nu,psrad,d2densi,npp,wrkfd,ier)
+        call splder(tfn,nsfd,cfn,3,nu,psrad(npp),d2densi(1),1,wrkfd,ier)
 
         psnpp=psrad(npp)
-        denpp=densi(npp)
-        ddenpp=ddensi(npp)
-        d2denpp=d2densi(npp)
+        denpp=densi(1)
+        ddenpp=ddensi(1)
+        d2denpp=d2densi(1)
         psdbnd=psnpp
         if(ddenpp.lt.0.0_wp_) then
           xnv=psnpp-ddenpp/d2denpp
@@ -2080,7 +2005,7 @@ c     spline approximation of density
 
       end if
 
-      deallocate(iwrkf,wrkf,wf,densi,ddensi,d2densi,wrkfd)
+      deallocate(iwrkf,wrkf,wf,wrkfd)
       
       return
       end
@@ -2498,7 +2423,7 @@ c
       use const_and_precisions, only : wp_,pi
       use graydata_par, only : rwallm
       use magsurf_data, only : npoints
-      use interp_eqprof, only : psiant,psinop,nsr=>nsrt,nsz=>nszt,
+      use interp_eqprof, only : psiant,psinop,nsr,nsz,
      . cc=>cceq,tr,tz,rup,zup,rlw,zlw,nr
       use dierckx, only : profil,sproota
       implicit none
@@ -2566,7 +2491,8 @@ c
       subroutine flux_average
       use const_and_precisions, only : wp_,zero,one,pi,ccj=>mu0inv
       use magsurf_data
-      use interp_eqprof, only : btrcen,btaxis,rmaxis,zmaxis,phitedge
+      use interp_eqprof, only : btrcen,btaxis,rmaxis,zmaxis,phitedge,
+     . equinum_fpol
       use simplespline, only : difcs
       use dierckx, only : regrid,coeff_parder
       implicit none
@@ -2593,9 +2519,7 @@ c local variables
       real(wp_), dimension(lwrk) :: wrk
       real(wp_), dimension(:), allocatable :: rctemp,zctemp
 c common/external functions/variables
-      real(wp_) :: fpolv,ffpv
-c
-      common/fpol/fpolv,ffpv
+      real(wp_) :: fpolv
 c
       npsi=nnintp
       ninpr=(npsi-1)/10
@@ -2696,7 +2620,7 @@ c
 
         call tor_curr(rctemp(1),zctemp(1),ajphi0)
         call bpol(rctemp(1),zctemp(1),brr,bzz)
-        call equinum_fpol(0)
+        call equinum_fpol(height2,fpolv)
         bphi=fpolv/rctemp(1)
         btot0=sqrt(bphi**2+brr**2+bzz**2)
         bpoloid0=sqrt(brr**2+bzz**2)
@@ -2727,7 +2651,7 @@ c         compute line integral on the contour psi=height^2
           zpsim=zctemp(inc1)
           call bpol(rpsim,zpsim,brr,bzz)
           call tor_curr(rpsim,zpsim,ajphi)
-          call equinum_fpol(0)
+!          call equinum_fpol(height2,fpolv)
           bphi=fpolv/rpsim
           btot=sqrt(bphi**2+brr**2+bzz**2)
           bpoloid=sqrt(brr**2+bzz**2)
@@ -2995,10 +2919,10 @@ c local variables
       real(wp_), dimension(lwrk) :: wrk
       real(wp_), dimension(:), allocatable :: rctemp,zctemp
 c common/external functions/variables
-      real(wp_) :: dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv
+      real(wp_) :: dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,dfpv
       real(wp_) :: frhotor_an
 c
-      common/fpol/fpolv,ffpv
+      common/fpol/fpolv,dfpv
       common/derip1/dpsidr,dpsidz
       common/derip2/ddpsidrr,ddpsidzz,ddpsidrz
 c
@@ -4169,18 +4093,19 @@ c
       use const_and_precisions, only : wp_,pi,ccj=>mu0inv
       use graydata_flags, only : iequil
       use graydata_par, only : sgnbphi
+      use interp_eqprof, only : equinum_fpol, equinum_psi
       implicit none
 c arguments
       real(wp_) :: xx,yy,zz
 c local variables
       integer :: iv,jv
-      real(wp_) :: b2tot,csphi,dfpolv,drrdx,drrdy,dphidx,dphidy,
+      real(wp_) :: b2tot,csphi,drrdx,drrdy,dphidx,dphidy,
      . rr,rr2,rrm,snphi,zzm
       real(wp_), dimension(3) :: dbtot,bvc
       real(wp_), dimension(3,3) :: dbvcdc,dbvdc,dbv
 c common/external functions/variables
       real(wp_) :: bres,brr,bphi,bzz,ajphi,btot,xg,yg,dxgdpsi,dpsidr,
-     . dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv,psinv
+     . dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,dfpv,psinv
       real(wp_), dimension(3) :: bv,derxg,deryg
       real(wp_), dimension(3,3) :: derbv
 c
@@ -4195,7 +4120,7 @@ c
       common/dxgyg/derxg,deryg
       common/derip1/dpsidr,dpsidz
       common/derip2/ddpsidrr,ddpsidzz,ddpsidrz
-      common/fpol/fpolv,ffpv
+      common/fpol/fpolv,dfpv
       common/psival/psinv
 
       xg=0.0_wp_
@@ -4237,16 +4162,16 @@ c
       end if
 c
       if(iequil.eq.2) then
-        call equinum_psi(rrm,zzm)
-        call equinum_derpsi(rrm,zzm,3)
-        call equinum_fpol(1)
+        call equinum_psi(rrm,zzm,psinv,dpsidr,dpsidz,ddpsidrr,ddpsidzz,
+     .      ddpsidrz)
+        call equinum_fpol(psinv,fpolv,dfpv)
       end if
 
-        call sub_xg_derxg
-        yg=fpolv/(rrm*bres)
-        bphi=fpolv/rrm
-        btot=abs(bphi)
-        if(psinv.lt.0.0_wp_) return
+      call sub_xg_derxg(psinv)
+      yg=fpolv/(rrm*bres)
+      bphi=fpolv/rrm
+      btot=abs(bphi)
+      if(psinv.lt.0.0_wp_) return
 c
 c     B = f(psi)/R e_phi+  grad psi x e_phi/R
 c
@@ -4256,8 +4181,6 @@ c
       bphi=fpolv/rrm
       brr=-dpsidz/rrm
       bzz= dpsidr/rrm
-c
-      dfpolv=ffpv/fpolv
 c
       bvc(1)=brr
       bvc(2)=bphi
@@ -4274,8 +4197,8 @@ c     dbvcdc(iv,jv) = d Bcil(iv) / dxvcil(jv)
 c
       dbvcdc(1,1)=-ddpsidrz/rrm-brr/rrm
       dbvcdc(1,3)=-ddpsidzz/rrm
-      dbvcdc(2,1)= dfpolv*dpsidr/rrm-bphi/rrm
-      dbvcdc(2,3)= dfpolv*dpsidz/rrm
+      dbvcdc(2,1)= dfpv*dpsidr/rrm-bphi/rrm
+      dbvcdc(2,3)= dfpv*dpsidz/rrm
       dbvcdc(3,1)= ddpsidrr/rrm-bzz/rrm
       dbvcdc(3,3)= ddpsidrz/rrm
 c
@@ -4352,21 +4275,19 @@ c
 c arguments
       real(wp_) :: rrm,zzm
 c local variables
-      real(wp_) :: cst,dpsidrp,d2psidrp,dqq,qq,rn,rpm,sgn,snt,dfpolv,
-     . rhop,rhot
+      real(wp_) :: cst,dpsidrp,d2psidrp,dqq,qq,rn,rpm,sgn,snt,rhop,rhot
 c common/external functions/variables
       real(wp_) :: psinv,dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,
-     . ffpv,xg,xgcn,dxgdpsi,dens,ddens
+     . dfpv,xg,xgcn,dxgdpsi,dens,ddens
       real(wp_) :: frhopol
 c
       common/psival/psinv
       common/derip1/dpsidr,dpsidz
       common/derip2/ddpsidrr,ddpsidzz,ddpsidrz
-      common/fpol/fpolv,ffpv
+      common/fpol/fpolv,dfpv
       common/xgxg/xg
       common/dxgdps/dxgdpsi
       common/xgcn/xgcn
-      common/dens/dens,ddens
 c
       dpsidrp=0.0_wp_
       d2psidrp=0.0_wp_
@@ -4405,8 +4326,7 @@ c
       end if
 c
       fpolv=sgnbphi*b0*rr0m
-      dfpolv=0.0_wp_
-      ffpv=0.0_wp_
+      dfpv=0.0_wp_
 c
       dpsidr=dpsidrp*cst
       dpsidz=dpsidrp*snt
@@ -4418,209 +4338,6 @@ c
       end
 c
 c
-c
-      subroutine equinum_psi(rpsim,zpsim)
-      use const_and_precisions, only : wp_
-      use interp_eqprof, only : rmnm,rmxm,zmnm,zmxm,psiant,psinop,
-     . tr,tz,ccspl=>cceq,nsrt,nszt
-      use dierckx, only : fpbisp
-      implicit none
-c local constants
-      integer, parameter :: lwrk=8,liwrk=2,nrs=1,nzs=1
-c arguments
-      real(wp_) :: rpsim,zpsim
-c local variables
-      integer :: nsr,nsz
-      integer, dimension(liwrk) :: iwrk
-      real(wp_), dimension(1) :: rrs,zzs,ffspl
-      real(wp_), dimension(lwrk) :: wrk
-c common/external functions/variables
-      real(wp_) :: psinv
-c
-      common/psival/psinv
-c
-      psinv=-1.0_wp_
-c
-c     here lengths are measured in meters
-c
-      if (rpsim.gt.rmxm.or.rpsim.lt.rmnm) return
-      if (zpsim.gt.zmxm.or.zpsim.lt.zmnm) return
-c
-      rrs(1)=rpsim
-      zzs(1)=zpsim
-      nsr=nsrt
-      nsz=nszt
-      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
-      return
-      end
-c
-c
-c
-      subroutine equinum_derpsi(rpsim,zpsim,iderpsi)
-      use const_and_precisions, only : wp_
-      use interp_eqprof, only : nr,nz,rmnm,rmxm,zmnm,zmxm,cc01=>cceq01,
-     . cc10=>cceq10,cc20=>cceq20,cc02=>cceq02,cc11=>cceq11
-      implicit none
-c arguments
-      integer :: iderpsi
-      real(wp_) :: rpsim,zpsim
-c local variables
-      integer :: lw10,lw01,lw20,lw02,lw11,nur,nuz
-      real(wp_) :: derpsi
-c common/external functions/variables
-      real(wp_) :: dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz
-c
-      common/derip1/dpsidr,dpsidz
-      common/derip2/ddpsidrr,ddpsidzz,ddpsidrz
-c
-      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
-c
-      dpsidr=0.0_wp_
-      dpsidz=0.0_wp_
-      ddpsidrr=0.0_wp_
-      ddpsidzz=0.0_wp_
-      ddpsidrz=0.0_wp_
- 
-c     here lengths are measured in meters
-      if (rpsim.gt.rmxm.or.rpsim.lt.rmnm) return
-      if (zpsim.gt.zmxm.or.zpsim.lt.zmnm) return
-c
-      select case(iderpsi)
-
-      case(1)
-      nur=1
-      nuz=0
-      call sub_derpsi(rpsim,zpsim,nur,nuz,derpsi,cc10,lw10)
-      dpsidr=derpsi     
-      nur=0
-      nuz=1
-      call sub_derpsi(rpsim,zpsim,nur,nuz,derpsi,cc01,lw01)
-      dpsidz=derpsi
-
-      case(2)
-      nur=2
-      nuz=0
-      call sub_derpsi(rpsim,zpsim,nur,nuz,derpsi,cc20,lw20)
-      ddpsidrr=derpsi      
-      nur=0
-      nuz=2
-      call sub_derpsi(rpsim,zpsim,nur,nuz,derpsi,cc02,lw02)
-      ddpsidzz=derpsi      
-      nur=1
-      nuz=1
-      call sub_derpsi(rpsim,zpsim,nur,nuz,derpsi,cc11,lw11)
-      ddpsidrz=derpsi
-
-      case(3)
-      nur=1
-      nuz=0
-      call sub_derpsi(rpsim,zpsim,nur,nuz,derpsi,cc10,lw10)
-      dpsidr=derpsi      
-      nur=0
-      nuz=1
-      call sub_derpsi(rpsim,zpsim,nur,nuz,derpsi,cc01,lw01)
-      dpsidz=derpsi
-      nur=2
-      nuz=0
-      call sub_derpsi(rpsim,zpsim,nur,nuz,derpsi,cc20,lw20)
-      ddpsidrr=derpsi      
-      nur=0
-      nuz=2
-      call sub_derpsi(rpsim,zpsim,nur,nuz,derpsi,cc02,lw02)
-      ddpsidzz=derpsi      
-      nur=1
-      nuz=1
-      call sub_derpsi(rpsim,zpsim,nur,nuz,derpsi,cc11,lw11)
-      ddpsidrz=derpsi
-
-      case default
-      print*,'iderpsi undefined'
-
-      end select
-c
-      return
-      end
-c
-c
-c
-      subroutine sub_derpsi(rpsim,zpsim,nur,nuz,derpsi,cc,lw)
-      use const_and_precisions, only : wp_
-      use dierckx, only : fpbisp
-      use interp_eqprof, only : nr,nz,psia,tr,tz,nsrt,nszt
-      implicit none
-c local constants
-      integer, parameter :: liwrk=2,nrs=1,nzs=1
-c arguments
-      integer :: nur,nuz,lw
-      real(wp_) :: rpsim,zpsim,derpsi
-      real(wp_), dimension(lw) :: cc
-c local variables
-      integer :: iwr,iwz,kkr,kkz,nsr,nsz
-      integer, dimension(liwrk) :: iwrk
-      real(wp_), dimension(1) :: rrs,zzs,ffspl
-c
-      rrs(1)=rpsim
-      zzs(1)=zpsim
-      nsr=nsrt
-      nsz=nszt
-      kkr=3-nur
-      kkz=3-nuz
-      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,kkr,kkz
-     .     ,rrs,nrs,zzs,nzs,ffspl,cc(iwr),cc(iwz),iwrk(1),iwrk(2))
-      derpsi=ffspl(1)*psia
-c
-      return
-      end
-c
-c
-c
-      subroutine equinum_fpol(iderfpol)
-      use const_and_precisions, only : wp_
-      use dierckx, only : splev,splder
-      use interp_eqprof, only : nr,psia,tfp,nsft,cfp,fpolas
-      implicit none
-c arguments
-      integer :: iderfpol
-c local variables
-      integer :: ier
-      real(wp_) :: dfpolv
-      real(wp_), dimension(1) :: rrs,ffspl
-      real(wp_), dimension(nr+4) :: wrkfd
-c common/external functions/variables
-      real(wp_) :: psinv,fpolv,ffpv
-c
-      common/psival/psinv
-      common/fpol/fpolv,ffpv
-c      
-      fpolv=fpolas
-      dfpolv=0.0_wp_
-      ffpv=0.0_wp_
-      if(iderfpol.lt.0.or.iderfpol.gt.1) return
-c
-      if(psinv.le.1.0_wp_.and.psinv.gt.0.0_wp_) then
-        rrs(1)=psinv
-        call splev(tfp,nsft,cfp,3,rrs,ffspl,1,ier)
-        fpolv=ffspl(1)
-        if(iderfpol.eq.1) then       
-          call splder(tfp,nsft,cfp,3,1,rrs,ffspl,1,wrkfd,ier)
-          dfpolv=ffspl(1)
-          ffpv=fpolv*dfpolv/psia
-        end if
-      end if
-c
-      return
-      end
-c
-c
 c
       subroutine bfield(rpsim,zpsim,bphi,brr,bzz)
       use const_and_precisions, only : wp_
@@ -4638,17 +4355,15 @@ c
 c
       subroutine btor(rpsim,zpsim,bphi)
       use const_and_precisions, only : wp_
+      use interp_eqprof, only : equinum_fpol, equinum_psi
       implicit none
 c arguments
       real(wp_) :: rpsim,zpsim,bphi
-c common/external functions/variables
-      real(wp_) :: psinv,fpolv,ffpv
+c local variables
+      real(wp_) :: psinv,fpolv
 c
-      common/psival/psinv
-      common/fpol/fpolv,ffpv
-c
-      call equinum_psi(rpsim,zpsim)
-      call equinum_fpol(0)
+      call equinum_psi(rpsim,zpsim,psinv)
+      call equinum_fpol(psinv,fpolv)
       bphi=fpolv/rpsim
 c
       return
@@ -4658,15 +4373,14 @@ c
 c
       subroutine bpol(rpsim,zpsim,brr,bzz)
       use const_and_precisions, only : wp_
+      use interp_eqprof, only : equinum_psi
       implicit none
 c arguments
       real(wp_) :: rpsim,zpsim,brr,bzz
-c common/external functions/variables
+c local variables
       real(wp_) :: dpsidr,dpsidz
 c
-      common/derip1/dpsidr,dpsidz
-c
-      call equinum_derpsi(rpsim,zpsim,1)
+      call equinum_psi(rpsim,zpsim,dpsidr=dpsidr,dpsidz=dpsidz)
       brr=-dpsidz/rpsim
       bzz= dpsidr/rpsim
 c
@@ -4694,18 +4408,16 @@ c
 c
       subroutine tor_curr(rpsim,zpsim,ajphi)
       use const_and_precisions, only : wp_,ccj=>mu0inv
+      use interp_eqprof, only : equinum_psi
       implicit none
 c arguments
       real(wp_) :: rpsim,zpsim,ajphi
 c local variables
       real(wp_) :: bzz,dbvcdc13,dbvcdc31
-c common/external functions/variables
-      real(wp_) :: dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz
+      real(wp_) :: dpsidr,ddpsidrr,ddpsidzz
 c
-      common/derip1/dpsidr,dpsidz
-      common/derip2/ddpsidrr,ddpsidzz,ddpsidrz
-c
-      call equinum_derpsi(rpsim,zpsim,3)
+      call equinum_psi(rpsim,zpsim,dpsidr=dpsidr,ddpsidrr=ddpsidrr,
+     .    ddpsidzz=ddpsidzz)
       bzz= dpsidr/rpsim
       dbvcdc13=-ddpsidzz/rpsim
       dbvcdc31= ddpsidrr/rpsim-bzz/rpsim
@@ -4718,8 +4430,8 @@ c
 c
       subroutine psi_raxis(h,r1,r2)
       use const_and_precisions, only : wp_
-      use interp_eqprof, only : psiant,psinop,zmaxis,nsr=>nsrt,
-     . nsz=>nszt,cc=>cceq,tr,tz,nr
+      use interp_eqprof, only : psiant,psinop,zmaxis,nsr,
+     . nsz,cc=>cceq,tr,tz,nr
       use dierckx, only : profil,sproota
       implicit none
 c local constants
@@ -4746,14 +4458,17 @@ c
 c
 c
 c
-      subroutine sub_xg_derxg
+      subroutine sub_xg_derxg(psinv)
       use const_and_precisions, only : wp_
       use interp_eqprof, only : psia
+      use coreprofiles, only : density
       implicit none
+c arguments
+      real(wp_) :: psinv
 c common/external functions/variables
-      real(wp_) :: psinv,xg,xgcn,dxgdpsi,dens,ddenspsin
+      real(wp_) :: xg,xgcn,dxgdpsi
+      real(wp_) :: dens,ddenspsin
 c
-      common/psival/psinv
       common/xgxg/xg
       common/dxgdps/dxgdpsi
       common/xgcn/xgcn
@@ -4762,7 +4477,7 @@ c
       xg=0.0_wp_
       dxgdpsi=0.0_wp_
 c      if(psinv.le.psdbnd.and.psinv.ge.0) then
-        call density(psinv)
+        call density(psinv,dens,ddenspsin)
         xg=xgcn*dens
         dxgdpsi=xgcn*ddenspsin/psia
 c      end if
@@ -4771,139 +4486,6 @@ c
       end
 c
 c
-c
-      subroutine density(arg)
-      use const_and_precisions, only : wp_
-      use graydata_flags, only : iprof
-      use graydata_par, only : psdbnd
-      use graydata_anequil, only : dens0,aln1,aln2
-      use interp_eqprof, only : psnpp,denpp,ddenpp,d2denpp,tfn,nsfd,
-     . cfn,npp
-      use dierckx, only : splev,splder
-      implicit none
-c arguments
-      real(wp_) :: arg
-c local variables
-      integer :: ier,nn,nn1,nn2,nu
-      real(wp_) :: profd,dprofd,dpsib,tt,fp,dfp,fh,dfh
-      real(wp_), dimension(1) :: xxs,ffs
-      real(wp_), dimension(npp+4) :: wrkfd
-c common/external functions/variables
-      real(wp_) :: dens,ddens
-c
-      common/dens/dens,ddens
-c
-c     computation of density [10^19 m^-3] and derivative wrt psi
-c
-      dens=0.0_wp_
-      ddens=0.0_wp_
-      if(arg.gt.psdbnd.or.arg.lt.0.0_wp_) return
-c
-      if(iprof.eq.0) then
-        if(arg.gt.1.0_wp_) return
-        profd=(1.0_wp_-arg**aln1)**aln2
-        dens=dens0*profd
-        dprofd=-aln1*aln2*arg**(aln1-1.0_wp_)
-     .         *(1.0_wp_-arg**aln1)**(aln2-1.0_wp_)
-        ddens=dens0*dprofd
-      else
-        if(arg.le.psdbnd.and.arg.gt.psnpp) then
-
-c     smooth interpolation for psnpp < psi < psdbnd
-c     dens = fp * fh
-c     fp: parabola matched at psi=psnpp with given profile density
-c     fh=(1-t)^3(1+3t+6t^2) is a smoothing function: 
-c     fh(0)=1, fh(1)=0 and zero first and second deriv at t=0,1 
-
-          fp=denpp+(arg-psnpp)*ddenpp+0.5_wp_*(arg-psnpp)**2*d2denpp
-          dfp=ddenpp+(arg-psnpp)*d2denpp
-          tt=(arg-psnpp)/(psdbnd-psnpp)
-          fh=(1.0_wp_-tt)**3*(1.0_wp_+3.0_wp_*tt+6.0_wp_*tt*tt)
-          dfh=-30.0_wp_*(1.0_wp_-tt)**2*tt*tt/(psdbnd-psnpp)
-          dens=fp*fh
-          ddens=dfp*fh+fp*dfh
-        else
-          xxs(1)=arg
-          ier=0
-          call splev(tfn,nsfd,cfn,3,xxs,ffs,1,ier)
-          dens=ffs(1)
-          nu=1
-          ier=0
-          call splder(tfn,nsfd,cfn,3,nu,xxs,ffs,1,wrkfd,ier)
-          ddens=ffs(1)
-          if(ier.gt.0) print*,ier
-          if(abs(dens).lt.1.0e-10_wp_) dens=0.0_wp_
-        end if
-        if(dens.lt.0.0_wp_) print*,'    DENSITY NEGATIVE',dens
-c
-      end if
-c
-      return
-      end
-c
-c
-c
-      function temp(arg)
-      use const_and_precisions, only : wp_,zero,one
-      use graydata_flags, only : iprof
-      use graydata_anequil, only : te0,dte0,alt1,alt2
-      use interp_eqprof, only : psrad,ct,npp
-      use utils, only : locate
-      use simplespline, only :spli
-      implicit none
-c arguments
-      real(wp_) :: arg,temp
-c local variables
-      integer :: k
-      real(wp_) :: proft,dps
-c
-      temp=zero
-      if(arg.ge.one.or.arg.lt.zero) return
-      if(iprof.eq.0) then
-        proft=(1.0_wp_-arg**alt1)**alt2
-        temp=(te0-dte0)*proft+dte0
-      else
-        call locate(psrad,npp,arg,k)
-        k=max(1,min(k,npp-1))
-        dps=arg-psrad(k)
-        temp=spli(ct,npp,k,dps)
-      endif
-c
-      return
-      end
-c
-c
-c
-      function fzeff(ps)
-      use const_and_precisions, only : wp_,zero,one
-      use graydata_flags, only : iprof
-      use graydata_anequil, only : zeffan
-      use interp_eqprof, only : psrad,cz,npp
-      use utils, only : locate
-      use simplespline, only :spli
-      implicit none
-c arguments
-      real(wp_), intent(in) :: ps
-      real(wp_) :: fzeff
-c local variables
-      integer :: k
-      real(wp_) :: dps
-c
-      fzeff=one
-      if(ps.gt.one.or.ps.lt.zero) return
-      if(iprof.eq.0) then
-        fzeff=zeffan
-      else
-        call locate(psrad,npp,ps,k)
-        k=max(1,min(k,npp-1))
-        dps=ps-psrad(k)
-        fzeff=spli(cz,npp,k,dps)
-      endif
-c
-      return
-      end
-c
-c
 c
       subroutine ic_gb
 c     beam tracing initial conditions  igrad=1
@@ -5536,7 +5118,7 @@ c
       use beamdata, only : nrayr,nrayth,q
       implicit none
 c local variables
-      integer :: j,k
+      integer :: j
       real(wp_) :: dr,r0,r,w,w0,summ
 c
       q(1)=1.0_wp_
@@ -5618,6 +5200,7 @@ c
       use graydata_anequil, only : rr0m
       use beamdata, only : psjki,tauv,alphav,pdjki,ppabs,currj,didst,
      . ccci,q,tau1v
+      use coreprofiles, only : temp, fzeff
       use dispersion, only : harmnumber, warmdisp
       use eccd, only : setcdcoeff,eccdeff,fjch0,fjch,fjncl
       implicit none
@@ -5626,8 +5209,7 @@ c arguments
 c local constants
       real(wp_), parameter :: taucr=12.0_wp_,xxcr=16.0_wp_,eps=1.e-8_wp_
 c local variables
-      integer :: intp
-      real(wp_) :: dvvl,rbavi,rrii,rhop0,rhop,voli,dervoli,area,
+      real(wp_) :: dvvl,rbavi,rrii,rhop0,rhop,dervoli,
      . tau,effjcdav,dpdst,p0,pow,alpha0,didst0,ccci0
       integer :: lrm,ithn
       real(wp_) :: amu,ratiovgr,rbn,rbx,zeff
@@ -5640,7 +5222,6 @@ c common/external functions/variables
       real(wp_) :: anpl,anpr,vgm,derdnm,sox,anprre,anprim,alpha,effjcd,
      .             akim,tau0
       complex(wp_) :: ex,ey,ez
-      real(wp_), external :: fzeff,temp
 c
       common/p0/p0mw
       common/dersdst/dersdst
@@ -5884,8 +5465,7 @@ c
       integer :: it,ipec
       real(wp_) :: drt,rt,rt1,rhop1
       real(wp_) :: ratjai,ratjbi,ratjpli
-      real(wp_) :: voli0,voli1,dervoli,areai0,areai1
-      real(wp_) :: rrii,rbavi,bmxi,bmni,fci
+      real(wp_) :: voli0,voli1,areai0,areai1
       real(wp_), dimension(nnd), intent(in)  ::  rt_in
       real(wp_), dimension(nnd), intent(out)  :: rhop_tab,rhot_tab
       real(wp_), dimension(0:nnd), intent(out) :: rtabpsi1
@@ -5949,7 +5529,7 @@ c     psi grid at mid points, dimension nnd+1, for use in pec_tab
       subroutine spec(pabs,currt,rtabpsi1,dvol,darea,dpdv,ajphiv,
      .                pins,currins)
       use const_and_precisions, only : wp_,one,zero,izero
-      use graydata_flags, only : nnd,ipec,ieccd
+      use graydata_flags, only : nnd
       use beamdata, only : ppabs,ccci,psjki,nrayr,nrayth,nstep,iiv
       implicit none
 c local constants
@@ -6134,7 +5714,6 @@ c     radial average values over power and current density profile
       real(wp_), dimension(nnd), intent(in) :: dpdv,ajphiv
       real(wp_), external :: frhopol,fdadrhot,fdvdrhot
     
-      integer :: i
       real(wp_) :: sccsa
       real(wp_) :: rhotjav,rhot2pav,rhot2java
     
@@ -6192,7 +5771,6 @@ c     factor sqrt(8) = 2sqrt(2) to match full width of gaussian profile
 
       subroutine profwidth(nd,xx,yy,xpk,ypk,dxxe)
       use const_and_precisions, only : wp_,emn1
-      use graydata_flags, only : ipec
       use utils, only : locatex, locate, intlin, vmaxmini
       implicit none
 c arguments
@@ -6246,7 +5824,11 @@ c
         rte1=0.0_wp_
         yye=ypk*emn1
         call locate(yy,nd,yye,ie)
-        call intlin(yy(ie),xx(ie),yy(ie+1),xx(ie+1),yye,rte2)
+        if(ie.gt.0.and.ie.lt.nd) then
+          call intlin(yy(ie),xx(ie),yy(ie+1),xx(ie+1),yye,rte2)
+        else
+          rte2=0.0_wp_
+        end if
       end if
       dxxe=rte2-rte1
       if(ymx.ne.0.and.ymn.ne.0) dxxe=-dxxe
@@ -6428,8 +6010,8 @@ c
       logical function inside_plasma(rrm,zzm)
       use const_and_precisions, only : wp_
       use graydata_flags, only : iequil
-      use graydata_par, only : psdbnd
-      use interp_eqprof, only : zbmin,zbmax
+      use coreprofiles, only : psdbnd
+      use interp_eqprof, only : zbmin,zbmax,equinum_psi
       implicit none
 c arguments
       real(wp_) :: rrm,zzm
@@ -6441,7 +6023,7 @@ c
       if(iequil.eq.1) then
         call equian(rrm,zzm)
       else
-        call equinum_psi(rrm,zzm)
+        call equinum_psi(rrm,zzm,psinv)
       end if
 c
       if (psinv.ge.0.0_wp_.and.psinv.lt.psdbnd) then
diff --git a/src/graydata_par.f90 b/src/graydata_par.f90
index 0d362a8..4e8a29c 100644
--- a/src/graydata_par.f90
+++ b/src/graydata_par.f90
@@ -4,7 +4,7 @@ module graydata_par
   
   real(wp_), save :: rwmax,rwallm
   real(wp_), save :: psipol0,chipol0
-  real(wp_), save :: factb,factt,factn,psdbnd
+  real(wp_), save :: factb,factt,factn
   real(wp_), save :: sgnbphi,sgniphi
 
 end module graydata_par
diff --git a/src/interp_eqprof.f90 b/src/interp_eqprof.f90
index 2f7bbe7..1da17e1 100644
--- a/src/interp_eqprof.f90
+++ b/src/interp_eqprof.f90
@@ -3,29 +3,53 @@ module interp_eqprof
   implicit none
   
 ! equidata
-  INTEGER, SAVE :: nlim,nr,nz,nbbbs,nsrt,nszt,nsft
-  REAL(wp_), SAVE :: psia,psiant,psinop,btrcen,rcen
-  REAL(wp_), SAVE :: btaxis,rmaxis,zmaxis,rmnm,rmxm,zmnm,zmxm,dr,dz
-  REAL(wp_), SAVE :: zbmin,zbmax,fpolas,phitedge,rrtor,rup,zup,rlw,zlw   !rrtor non usato, solo in equidata
-  REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: tr,tz,tfp
-  REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: cceq,cfp
-  REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: cceq01,cceq10,cceq20,cceq02,cceq11
-  REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: psinr,rhopnr,qpsi,rv,zv
-  REAL(wp_), DIMENSION(:,:), ALLOCATABLE, SAVE :: psin,psi,btotal,crhot,cq
-  REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: rlim,zlim,rbbbs,zbbbs
+  ! === 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
+  
+  ! === 1D array plasma boundary Rbnd_i, Zbnd_i ====================
+  INTEGER, SAVE :: nbbbs
+  REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: rbbbs,zbbbs
+  ! === 1D array limiter Rlim_i, Zlim_i   ==> move in wall/reflections
+  INTEGER, SAVE :: nlim
+  REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: rlim,zlim
 
-! profdata
-  INTEGER, SAVE :: npp,nsfd
-  REAL(wp_), SAVE :: psnpp,denpp,ddenpp,d2denpp
-  REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: psrad,derad,terad,zfc,tfn,cfn
-  REAL(wp_), DIMENSION(:,:), ALLOCATABLE, SAVE :: ct,cz
+  REAL(wp_), SAVE :: btaxis,rmaxis,zmaxis
+  REAL(wp_), SAVE :: rmnm,rmxm,zmnm,zmxm,dr,dz
+  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.
+
+  ! === 2D spline psi(R,z), normalization and derivatives ==========
+  INTEGER, SAVE :: nrest, nzest, lw10, lw01, lw20, lw02, lw11
+  INTEGER :: nsr,nsz
+  REAL(wp_), SAVE :: psia, psiant, psinop
+  REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: tr,tz
+  REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: cceq,   cceq01, cceq10, &
+                                                cceq20, cceq02, cceq11
+  ! === 1D spline Fpol(psi) ========================================
+  ! INTEGER, SAVE :: npsiest
+  INTEGER, SAVE :: nsf
+  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
+  REAL(wp_), DIMENSION(:,:), ALLOCATABLE, SAVE :: crhot,cq
+
+  ! !!! 2D B(R,z) array. Computed in bfield_res,
+  ! !!! used by cniteq to plot resonant field contour lines.
+  REAL(wp_), DIMENSION(:,:), ALLOCATABLE, SAVE :: btotal
 
 contains
 
   subroutine alloc_equilvec(ier)
     implicit none
     integer, intent(out) :: ier
-    integer :: nrest,nzest,lw10,lw01,lw20,lw02,lw11
 
     if(nr.le.0.or.nz.le.0) then
       ier = -1
@@ -43,7 +67,7 @@ contains
     call dealloc_equilvec
     allocate(rv(nr),zv(nz),tr(nrest),tz(nzest),tfp(nrest),cfp(nrest), &
      btotal(nr,nz),cceq(nr*nz),cceq01(lw01),cceq10(lw10),cceq02(lw02), &
-     cceq20(lw20),cceq11(lw11),psi(nr,nz),psin(nr,nz),psinr(nr),rhopnr(nr), &
+     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
@@ -63,9 +87,9 @@ contains
     if(allocated(cceq02)) deallocate(cceq02)
     if(allocated(cceq20)) deallocate(cceq20)
     if(allocated(cceq11)) deallocate(cceq11)
-    if(allocated(psi)) deallocate(psi)
     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)
 
@@ -117,38 +141,107 @@ contains
     if(allocated(zlim))   deallocate(zlim)
   end subroutine dealloc_lim
 
-
-  subroutine alloc_profvec(ier)
+  subroutine equinum_psi(rpsim,zpsim,psinv,dpsidr,dpsidz, &
+        ddpsidrr,ddpsidzz,ddpsidrz)
+    use dierckx, only : fpbisp
     implicit none
-    integer, intent(out) :: ier
-    integer :: npest
+! local constants
+    integer, parameter :: lwrk=8,liwrk=2
+! arguments
+    real(wp_), intent(in) :: rpsim,zpsim
+    real(wp_), intent(out), optional :: psinv,dpsidr,dpsidz, &
+        ddpsidrr,ddpsidzz,ddpsidrz
 
-    ier=0
-    if(npp.le.0) then
-      ier = -1
-      return
+! local variables
+    integer, dimension(liwrk) :: iwrk
+    real(wp_), dimension(1) :: rrs,zzs,ffspl
+    real(wp_), dimension(lwrk) :: wrk
+
+!
+!     here lengths are measured in meters
+!
+    if (rpsim.le.rmxm .and. rpsim.ge.rmnm .and. &
+        zpsim.le.zmxm .and. zpsim.ge.zmnm) then
+
+      if (present(psinv)) then
+        rrs(1)=rpsim
+        zzs(1)=zpsim
+        call fpbisp(tr,nsr,tz,nsz,cceq,3,3,rrs,1,zzs,1,ffspl, &
+                   wrk(1),wrk(5),iwrk(1),iwrk(2))
+        psinv=(ffspl(1)-psinop)/psiant
+      end if
+      if (present(dpsidr)) then
+        call sub_derpsi(rpsim,zpsim,1,0,dpsidr,cceq10,lw10)
+      end if
+      if (present(dpsidz)) then
+        call sub_derpsi(rpsim,zpsim,0,1,dpsidz,cceq01,lw01)
+      end if
+      if (present(ddpsidrr)) then
+        call sub_derpsi(rpsim,zpsim,2,0,ddpsidrr,cceq20,lw20)
+      end if
+      if (present(ddpsidzz)) then
+        call sub_derpsi(rpsim,zpsim,0,2,ddpsidzz,cceq02,lw02)
+      end if
+      if (present(ddpsidrz)) then
+        call sub_derpsi(rpsim,zpsim,1,1,ddpsidrz,cceq11,lw11)
+      end if
+    else
+      if(present(psinv)) psinv=-1.0_wp_
+      if(present(dpsidr)) dpsidr=0.0_wp_
+      if(present(dpsidz)) dpsidz=0.0_wp_
+      if(present(ddpsidrr)) ddpsidrr=0.0_wp_
+      if(present(ddpsidzz)) ddpsidzz=0.0_wp_
+      if(present(ddpsidrz)) ddpsidrz=0.0_wp_
     end if
+  end subroutine equinum_psi
 
-    npest=npp+4
-
-    call dealloc_profvec
-    allocate(psrad(npp),terad(npp),derad(npp),zfc(npp),ct(npp,4), &
-     cz(npp,4),tfn(npest),cfn(npest), &
-     stat=ier)
-    if (ier/=0) call dealloc_profvec
-  end subroutine alloc_profvec
-
-  subroutine dealloc_profvec
+  subroutine sub_derpsi(rpsim,zpsim,nur,nuz,derpsi,cc,lw)
+    use dierckx, only : fpbisp
     implicit none
-    if(allocated(psrad)) deallocate(psrad)
-    if(allocated(terad)) deallocate(terad)
-    if(allocated(derad)) deallocate(derad)
-    if(allocated(zfc)) deallocate(zfc)
-    if(allocated(ct)) deallocate(ct)
-    if(allocated(cz)) deallocate(cz)
-    if(allocated(tfn)) deallocate(tfn)
-    if(allocated(cfn)) deallocate(cfn)
+! 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
+! local variables
+    integer :: iwr,iwz
+    integer, dimension(liwrk) :: iwrk
+    real(wp_), dimension(1) :: rrs,zzs,ffspl
 
-  end subroutine dealloc_profvec
+    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))
+    derpsi=ffspl(1)*psia
+  end subroutine sub_derpsi
 
+  subroutine equinum_fpol(psinv,fpolv,dfpv)
+    use dierckx, only : splev,splder
+    implicit none
+! arguments
+    real(wp_), intent(in) :: psinv
+    real(wp_), intent(out) :: fpolv
+    real(wp_), intent(out), optional :: dfpv
+! local variables
+    integer :: ier
+    real(wp_), dimension(1) :: rrs,ffspl
+    real(wp_), dimension(nsf) :: wrkfd
+!      
+    if(psinv.le.1.0_wp_.and.psinv.gt.0.0_wp_) then
+      rrs(1)=psinv
+      call splev(tfp,nsf,cfp,3,rrs,ffspl,1,ier)
+      fpolv=ffspl(1)
+      if(present(dfpv)) then       
+        call splder(tfp,nsf,cfp,3,1,rrs,ffspl,1,wrkfd,ier)
+        dfpv=ffspl(1)/psia
+      end if
+    else
+      fpolv=fpolas
+      if (present(dfpv)) dfpv=0._wp_
+    end if
+  end subroutine equinum_fpol
+  
 end module interp_eqprof