diff --git a/src/coreprofiles.f90 b/src/coreprofiles.f90
index f1dfa89..a5ea467 100644
--- a/src/coreprofiles.f90
+++ b/src/coreprofiles.f90
@@ -1,50 +1,15 @@
 module coreprofiles
-  use const_and_precisions, only : wp_
+  use const_and_precisions, only : wp_,zero,one
   implicit none
 
   INTEGER, SAVE :: npp,nsfd
   REAL(wp_), SAVE :: psdbnd,psnpp,denpp,ddenpp,d2denpp
-  REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: psrad,derad,terad,zfc
-  REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: tfn,cfn
+  REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: tfn,cfn,psrad
   REAL(wp_), DIMENSION(:,:), ALLOCATABLE, SAVE :: ct,cz
 
 contains
 
-  subroutine alloc_profvec(ier)
-    implicit none
-    integer, intent(out) :: ier
-    integer :: npest
-
-    ier=0
-    if(npp.le.0) then
-      ier = -1
-      return
-    end if
-
-    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
-    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)
-
-  end subroutine dealloc_profvec
-
   subroutine density(psin,dens,ddens)
-    use const_and_precisions, only : wp_
     use graydata_flags, only : iprof
     use graydata_anequil, only : dens0,aln1,aln2
     use dierckx, only : splev,splder
@@ -62,16 +27,16 @@ contains
 !
 !     computation of density [10^19 m^-3] and derivative wrt psi
 !
-    dens=0.0_wp_
-    ddens=0.0_wp_
-    if(psin.ge.psdbnd.or.psin.lt.0.0_wp_) return
+    dens=zero
+    ddens=zero
+    if(psin.ge.psdbnd.or.psin.lt.zero) return
 !
     if(iprof.eq.0) then
-      if(psin.gt.1.0_wp_) return
-      profd=(1.0_wp_-psin**aln1)**aln2
+      if(psin.gt.one) return
+      profd=(one-psin**aln1)**aln2
       dens=dens0*profd
-      dprofd=-aln1*aln2*psin**(aln1-1.0_wp_)           &
-              *(1.0_wp_-psin**aln1)**(aln2-1.0_wp_)
+      dprofd=-aln1*aln2*psin**(aln1-one)           &
+              *(one-psin**aln1)**(aln2-one)
       ddens=dens0*dprofd
     else
       if(psin.gt.psnpp) then
@@ -86,8 +51,8 @@ contains
         fp=denpp+dpsib*ddenpp+0.5_wp_*dpsib**2*d2denpp
         dfp=ddenpp+dpsib*d2denpp
         tt=dpsib/(psdbnd-psnpp)
-        fh=(1.0_wp_-tt)**3*(1.0_wp_+3.0_wp_*tt+6.0_wp_*tt**2)
-        dfh=-30.0_wp_*(1.0_wp_-tt)**2*tt**2/(psdbnd-psnpp)
+        fh=(one-tt)**3*(one+3.0_wp_*tt+6.0_wp_*tt**2)
+        dfh=-30.0_wp_*(one-tt)**2*tt**2/(psdbnd-psnpp)
         dens=fp*fh
         ddens=dfp*fh+fp*dfh
       else
@@ -100,9 +65,9 @@ contains
         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_
+        if(abs(dens).lt.1.0e-10_wp_) dens=zero
       end if
-      if(dens.lt.0.0_wp_) print*,'    DENSITY NEGATIVE',dens
+      if(dens.lt.zero) print*,'    DENSITY NEGATIVE',dens
     end if
   end subroutine density
 
@@ -158,5 +123,148 @@ contains
       fzeff=spli(cz,npp,k,dps)
     endif
   end function fzeff
+  
+  subroutine read_profiles(filenm,psin,te,ne,zeff,unit)
+    use utils, only : get_free_unit
+    implicit none
+! arguments
+    character(len=*), intent(in) :: filenm
+    real(wp_), dimension(:), allocatable, intent(out) :: psin,te,ne,zeff
+    integer, optional, intent(in) :: unit
+! local variables
+    integer :: u, i, n
+    
+    if (present(unit)) then
+      u=unit
+    else
+      u=get_free_unit()
+    end if
+    open(file=trim(filenm),status='old',unit=u)
+    read(u,*) n
+    if(allocated(psin)) deallocate(psin)
+    if(allocated(te)) deallocate(te)
+    if(allocated(ne)) deallocate(ne)
+    if(allocated(zeff)) deallocate(zeff)
+    allocate(psin(n),te(n),ne(n),zeff(n))
+    do i=1,n
+      read(u,*) psin(i),te(i),ne(i),zeff(i)
+    end do
+    psin(1)=max(psin(1),zero)
+    close(u)
+  end subroutine read_profiles
+
+  subroutine tene_scal(te,ne,tfact,nfact,bfact,iscal)
+    implicit none
+! arguments
+    real(wp_), dimension(:), intent(inout) :: te,ne
+    real(wp_), intent(in) :: tfact,nfact,bfact
+    integer, intent(in) :: iscal
+! local variables
+    real(wp_) :: aat,aan,ffact
+
+    if (iscal==0) then
+      aat=2.0_wp_/3.0_wp_
+      aan=4.0_wp_/3.0_wp_
+    else
+      aat=1.0_wp_
+      aan=1.0_wp_
+    end if
+    if(iscal==2) then
+      ffact=1.0_wp_
+    else
+      ffact=bfact
+    end if
+    te(:)=te(:)*ffact**aat*tfact
+    ne(:)=ne(:)*ffact**aan*nfact
+  end subroutine tene_scal
+
+  subroutine set_prfspl(psin,te,ne,zeff,ssplne)
+    use simplespline, only : difcs
+    use dierckx, only : curfit, splev, splder
+    implicit none
+! arguments
+    real(wp_), dimension(:), intent(in) :: psin,te,ne,zeff
+    real(wp_), intent(in) :: ssplne
+! local variables
+    integer, parameter :: iopt=0, kspl=3
+    integer :: n, npest, lwrkf, ier
+    real(wp_) :: xb, xe, fp, xnv, ynv
+    real(wp_), dimension(:), allocatable :: wf, wrkf
+    integer, dimension(:), allocatable :: iwrkf
+    real(wp_), dimension(1) :: dedge,ddedge,d2dedge
+
+    n=size(psin)
+    npest=n+4
+    lwrkf=n*4+npest*16
+    allocate(wrkf(lwrkf),iwrkf(npest),wf(n))
+
+!   if necessary, reallocate spline arrays
+    if(.not.allocated(psrad)) then
+      allocate(psrad(n),ct(n,4),cz(n,4))
+    else
+      if(size(psrad)<n) then
+        deallocate(psrad,ct,cz)
+        allocate(psrad(n),ct(n,4),cz(n,4))
+      end if
+    end if
+    if(.not.allocated(cfn)) then
+      allocate(tfn(npest),cfn(npest))
+    else
+      if(size(cfn)<npest) then
+        deallocate(tfn,cfn)
+        allocate(tfn(npest),cfn(npest))
+      end if
+    end if
+
+!   spline approximation of  temperature and Zeff
+    call difcs(psin,te,  n,iopt,ct,ier)
+    call difcs(psin,zeff,n,iopt,cz,ier)
+    psrad=psin
+    npp=n
+
+!   spline approximation of density
+    xb=zero
+    xe=psin(n)
+    wf(:)=one
+    call curfit(iopt,n,psin,ne,wf,xb,xe,kspl,ssplne,npest, &
+       nsfd,tfn,cfn,fp,wrkf,lwrkf,iwrkf,ier)
+
+!   compute polinomial extrapolation matching the spline boundary up to the
+!   2nd order derivative, extending the profile up to psi=psdbnd where
+!   ne=ne'=ne''=0
+!   spline value and derivatives at the edge
+    call splev(tfn,nsfd,cfn,3,psin(n:n),dedge(1:1),1,ier)
+    call splder(tfn,nsfd,cfn,3,1,psin(n:n),ddedge(1:1), 1,wrkf(1:nsfd),ier)
+    call splder(tfn,nsfd,cfn,3,2,psin(n:n),d2dedge(1:1),1,wrkf(1:nsfd),ier)
+!   determination of the boundary
+    psnpp=psin(n)
+    denpp=dedge(1)
+    ddenpp=ddedge(1)
+    d2denpp=d2dedge(1)
+    psdbnd=psnpp
+    if(ddenpp.lt.zero) then
+      xnv=psnpp-ddenpp/d2denpp
+      ynv=denpp-0.5_wp_*ddenpp**2/d2denpp
+      if(d2denpp.gt.zero.and.ynv.ge.zero)  then
+        psdbnd=xnv
+      else
+        psdbnd=xnv+sqrt((ddenpp/d2denpp)**2-2.0_wp_*denpp/d2denpp)
+      end if  
+      print*,psnpp,denpp,ddenpp,d2denpp,xnv,ynv
+      print*,psdbnd
+    end if
+
+    deallocate(iwrkf,wrkf,wf)
+  end subroutine set_prfspl
+
+  subroutine unset_prfspl
+    implicit none
+
+    if(allocated(psrad)) deallocate(psrad)
+    if(allocated(ct)) deallocate(ct)
+    if(allocated(cz)) deallocate(cz)
+    if(allocated(tfn)) deallocate(tfn)
+    if(allocated(cfn)) deallocate(cfn)
+  end subroutine unset_prfspl
 
 end module coreprofiles
\ No newline at end of file
diff --git a/src/gray.f b/src/gray.f
index 10d025c..4f79a44 100644
--- a/src/gray.f
+++ b/src/gray.f
@@ -657,15 +657,16 @@ c
       use graydata_flags
       use graydata_par
       use graydata_anequil
-      use coreprofiles, only : psdbnd
+      use coreprofiles, only : psdbnd,read_profiles,tene_scal,set_prfspl
       use interp_eqprof, only : rmxm,zbmin,zbmax,btrcen,rcen
       use reflections, only : rlim,zlim,nlim,alloc_lim
       use beamdata, only : nrayr,nrayth,nstep
       implicit none
 c local variables
-      integer :: nfil,iox,ierr,leqq,lprf
+      integer :: nfil,iox,ierr,leqq
       real(wp_) :: dummy,bresg,r00,anr0c,anphi0c,fghz,zeff,
      . w0csi,w0eta,d0csi,d0eta,phi0,zrcsi,zreta
+      real(wp_), dimension(:), allocatable :: psrad,terad,derad,zfc
       character(len=8) :: wdat
       character(len=10) :: wtim
 c common/external functions/variables
@@ -908,11 +909,11 @@ c
 
       if (iprof.eq.1) then
         nprof=98
-        lprf=len_trim(filenmprf)
-        open(file=filenmprf(1:lprf)//'.prf',
-     .                    status= 'unknown',unit=nprof)
-        call profdata
-        close(nprof)
+        call read_profiles(trim(filenmprf)//'.prf',psrad,terad,derad,
+     .     zfc,nprof)
+        call tene_scal(terad,derad,factt,factn,factb,iscal)
+        call set_prfspl(psrad,terad,derad,zfc,0.001_wp_)
+        deallocate(psrad,terad,derad,zfc)
       end if
 c
 c     read equilibrium data  from file if iequil=2
@@ -1905,117 +1906,6 @@ c
  113  format(i6,12(1x,e12.5))
       end
 c
-c
-      subroutine profdata
-      use const_and_precisions, only : wp_,zero
-      use graydata_flags, only : iprof,iscal,nprof
-      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
-      integer :: ierr,i,iopt,ier,kspl,npest,lwrkf,nu,nn,nn1,nn2
-      integer, dimension(:), allocatable :: iwrkf
-      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
-      real(wp_), dimension(1) :: densi,ddensi,d2densi
-c
-c     read plasma profiles from file if iprof>0
-c
-      if (iscal.eq.0) then
-        aat=2.0_wp_/3.0_wp_
-        aan=4.0_wp_/3.0_wp_
-      else
-        aat=1.0_wp_
-        aan=1.0_wp_
-      end if
-      ffact=factb
-      if(iscal.eq.2) ffact=1.0_wp_
-c
-      if (iprof.gt.0) then
-        read(nprof,*) npp
-c
-        call alloc_profvec(ierr)
-        if (ierr.ne.0) stop
-c
-        npest=npp+4
-        lwrkf=npp*4+npest*16
-        if(allocated(wrkf)) deallocate(wrkf)
-        if(allocated(iwrkf)) deallocate(iwrkf)
-        if(allocated(wf)) deallocate(wf)
-        if(allocated(wrkfd)) deallocate(wrkfd)
-        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_
-        psrad(1)=psrad0
-        terad(1)=terad0*ffact**aat*factt
-        derad(1)=derad0*ffact**aan*factn
-        zfc(1)=zfc0
-        wf(1)=1.0_wp_
-        do i=2,npp
-          read(nprof,*) psradi,teradi,deradi,zfci
-          psrad(i)=psradi
-          terad(i)=teradi*ffact**aat*factt
-          derad(i)=deradi*ffact**aan*factn
-          zfc(i)=zfci
-          wf(i)=1.0_wp_
-        end do
-
-c     spline approximation of  temperature and Zeff
-
-        iopt=0
-        call difcs(psrad,terad,npp,iopt,ct,ier)
-
-        iopt=0
-        call difcs(psrad,zfc,npp,iopt,cz,ier)
-
-c     spline approximation of density
-
-        iopt=0
-        xb=0.0_wp_
-        xe=psrad(npp)
-        kspl=3
-        sspl=.001_wp_
-
-        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(npp),densi(1),1,ier)
-        nu=1
-        call splder(tfn,nsfd,cfn,3,nu,psrad(npp),ddensi(1),1,wrkfd,ier)
-        nu=2
-        call splder(tfn,nsfd,cfn,3,nu,psrad(npp),d2densi(1),1,wrkfd,ier)
-
-        psnpp=psrad(npp)
-        denpp=densi(1)
-        ddenpp=ddensi(1)
-        d2denpp=d2densi(1)
-        psdbnd=psnpp
-        if(ddenpp.lt.0.0_wp_) then
-          xnv=psnpp-ddenpp/d2denpp
-          ynv=denpp-0.5_wp_*ddenpp**2/d2denpp
-          if(d2denpp.gt.zero.and.ynv.ge.zero)  then
-            psdbnd=xnv
-          else
-            psdbnd=xnv+sqrt((ddenpp/d2denpp)**2-2.0_wp_*denpp/d2denpp)
-          end if  
-      print*,psnpp,denpp,ddenpp,d2denpp,xnv,ynv
-      print*,psdbnd
-        end if
-
-      end if
-
-      deallocate(iwrkf,wrkf,wf,wrkfd)
-      
-      return
-      end
-c
-c
 c
       subroutine rhotor(rhotsx)
       use const_and_precisions, only : wp_
diff --git a/src/utils.f90 b/src/utils.f90
index 84249ca..04a3157 100644
--- a/src/utils.f90
+++ b/src/utils.f90
@@ -246,4 +246,33 @@ contains
     end do
   end subroutine bubble
 
+  function get_free_unit(umin,umax) result(i)
+    implicit none
+    integer :: i
+    integer, intent(in), optional :: umin, umax
+    integer, parameter :: max_allowed = 999
+    integer :: ierr, iend
+    logical :: ex, op
+
+    if (present(umin)) then
+      i = max(0,umin)    ! start searching from unit min
+    else
+      i = 0
+    end if
+    if (present(umax)) then
+      iend = min(max(0,umax),max_allowed)
+    else
+      iend = max_allowed
+    end if
+    do
+      if (i>iend) then
+        i=-1  ! no free units found
+        exit
+      end if
+      inquire(unit=i,exist=ex,opened=op,iostat=ierr)
+      if (ierr==0.and.ex.and..not.op) exit  ! unit i exists and is not open
+      i = i + 1
+    end do
+  end function get_free_unit
+
 end module utils
\ No newline at end of file