diff --git a/src/coreprofiles.f90 b/src/coreprofiles.f90
new file mode 100644
index 0000000..f1dfa89
--- /dev/null
+++ b/src/coreprofiles.f90
@@ -0,0 +1,162 @@
+module coreprofiles
+  use const_and_precisions, only : wp_
+  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 :: 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
+    implicit none
+! arguments
+    real(wp_), intent(in) :: psin
+    real(wp_), intent(out) :: dens,ddens
+! local variables
+    integer, parameter :: nn=3, nn1=nn+1, nn2=nn+2
+    integer :: ier,nu
+    real(wp_) :: profd,dprofd,dpsib,tt,fp,dfp,fh,dfh
+    real(wp_), dimension(1) :: xxs,ffs
+    real(wp_), dimension(npp+4) :: wrkfd
+
+!
+!     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
+!
+    if(iprof.eq.0) then
+      if(psin.gt.1.0_wp_) return
+      profd=(1.0_wp_-psin**aln1)**aln2
+      dens=dens0*profd
+      dprofd=-aln1*aln2*psin**(aln1-1.0_wp_)           &
+              *(1.0_wp_-psin**aln1)**(aln2-1.0_wp_)
+      ddens=dens0*dprofd
+    else
+      if(psin.gt.psnpp) then
+!
+!     smooth interpolation for psnpp < psi < psdbnd
+!     dens = fp * fh
+!     fp: parabola matched at psi=psnpp with given profile density
+!     fh=(1-t)^3(1+3t+6t^2) is a smoothing function: 
+!     fh(0)=1, fh(1)=0 and zero first and second deriv at t=0,1 
+!
+        dpsib=psin-psnpp
+        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)
+        dens=fp*fh
+        ddens=dfp*fh+fp*dfh
+      else
+        xxs(1)=psin
+        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
+    end if
+  end subroutine density
+
+  function temp(psin)
+    use const_and_precisions, only : wp_,zero,one
+    use graydata_flags, only : iprof
+    use graydata_anequil, only : te0,dte0,alt1,alt2
+    use utils, only : locate
+    use simplespline, only :spli
+    implicit none
+! arguments
+    real(wp_), intent(in) :: psin
+    real(wp_) :: temp
+! local variables
+    integer :: k
+    real(wp_) :: proft,dps
+
+    temp=zero
+    if(psin.ge.one.or.psin.lt.zero) return
+    if(iprof.eq.0) then
+      proft=(1.0_wp_-psin**alt1)**alt2
+      temp=(te0-dte0)*proft+dte0
+    else
+      call locate(psrad,npp,psin,k)
+      k=max(1,min(k,npp-1))
+      dps=psin-psrad(k)
+      temp=spli(ct,npp,k,dps)
+    endif
+  end function temp
+
+  function fzeff(psin)
+    use const_and_precisions, only : wp_,zero,one
+    use graydata_flags, only : iprof
+    use graydata_anequil, only : zeffan
+    use utils, only : locate
+    use simplespline, only :spli
+    implicit none
+! arguments
+    real(wp_), intent(in) :: psin
+    real(wp_) :: fzeff
+! local variables
+    integer :: k
+    real(wp_) :: dps
+
+    fzeff=one
+    if(psin.gt.one.or.psin.lt.zero) return
+    if(iprof.eq.0) then
+      fzeff=zeffan
+    else
+      call locate(psrad,npp,psin,k)
+      k=max(1,min(k,npp-1))
+      dps=psin-psrad(k)
+      fzeff=spli(cz,npp,k,dps)
+    endif
+  end function fzeff
+
+end module coreprofiles
\ No newline at end of file