      subroutine t1(input,rad,ndata)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c     Transformation from DYNAMO/funcfl EAM file format
c     into format used by program BULK
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit double precision (a-h,o-z)
      parameter (np=10000)
      dimension da(np),db(np),dc(np),bb(np)
      dimension a(np),b(np),c(np),x(np),ro(np)
      character trash
      character*80 input
      character*20 embed_out,pair_out,den_out
c-----------------------------------------------------------------------
      embed_out='out.embed' !output files
      pair_out='out.pair'
      den_out='out.den'
c-----------------------------------------------------------------------
      open(1,file=input)
      do i=1,4
         read (1,*) trash
      end do
      read (1,*) ndata, rho_step, ntrash, d_step, rad
      read (1,*) trash
! creates distance and density data
      x(1)=0.0d0
      ro(1)=0.0d0
      dro=rho_step
      dx=d_step
      do i=1,ndata-1
         x(i+1)=real(i)*dx
         ro(i+1)=real(i)*dro
      end do
c-----------------------------------------------------------------------
! reads 1-st potential embed
      n1=ndata/5
      n2=ndata-n1*5
      n=1
      do i=1,n1
         read (1,*) a(n),a(n+1),a(n+2),a(n+3),a(n+4)
         n=n+5
      end do
      if (n2.gt.0) then
         read (1,*) (a(k),k=n,ndata)
      end if
! reads 2-nd potential rho
      n=1
      do i=1,n1
         read (1,*) b(n),b(n+1),b(n+2),b(n+3),b(n+4)
         n=n+5
      end do
      if (n2.gt.0) then
         read (1,*) (b(k),k=n,ndata)
      end if
! reads 3-rd potential r*pair
      n=1
      do i=1,n1
         read (1,*) c(n),c(n+1),c(n+2),c(n+3),c(n+4)
         n=n+5
      end do
      if (n2.gt.0) then
         read (1,*) (c(k),k=n,ndata)
      end if
      close(1)
! transformation to pair
      do i=2,ndata
         c(i)=c(i)/x(i)
      end do
      c(1)=c(2)
c-----------------------------------------------------------------------
! differentiation
      call derivacia(ro,a,ndata,da,dro)
      call derivacia(x,b,ndata,db,dx)
      call derivacia(x,c,ndata,dc,dx)
! output files
      open (2,file=embed_out)
      do i=1,ndata
         write (2,10) ro(i),a(i),da(i)
      end do
      close (2)
      open (2,file=den_out)
      do i=1,ndata
         write (2,10) x(i),b(i),db(i)
      end do
      close (2)
      open (2,file=pair_out)
      do i=1,ndata
         write (2,10) x(i),c(i),dc(i)
      end do
      close (2)
  10  format (F14.8,E24.9,E24.9)
      end
c
c========================= subroutines =================================
c
      subroutine derivacia (x,y,n,dy,dd)
      implicit double precision (a-h,o-z)
      dimension dy(n),x(n),y(n)
      do i=1,n
         dy(i)=0.0d0
      end do
      yp1=(y(2)-y(1))/dd
      ypn=(y(n)-y(n-1))/dd
      call spline (x,y,n,yp1,ypn,dy)
      end
c
c
c============================== downloaded subroutines =================
c
       SUBROUTINE spline(x,y,n,yp1,ypn,y2)
! Given arrays x(1:n) and y(1:n) containing a tabulated function, i.e., yi = f(xi), with
! x1 < x2 < ... < xN, and given values yp1 and ypn for the first derivative of the interpolating
! function at points 1 and n, respectively, this routine returns an array y2(1:n) of
! length n which contains the second derivatives of the interpolating function at the tabulated
! points xi. If yp1 and/or ypn are equal to 1x10^30 or larger, the routine is signaled to set
! the corresponding boundary condition for a natural spline, with zero second derivative on
! that boundary. Parameter: NMAX is the largest anticipated value of n.
!
! Please notice that the program spline is called only once to process an entire tabulated
! function in arrays xi and yi. Once this has been done, the values of the interpolated function
! for any value of x are obtained by calls (as many as desired) to a separate routine splint.
       INTEGER n,NMAX
       REAL*8 yp1,ypn,x(n),y(n),y2(n)
       PARAMETER (NMAX=500)
       INTEGER i,k
       REAL*8 p,qn,sig,un,u(n)
       if (yp1.gt..99e30) then
       y2(1)=0.
       u(1)=0.
       else
       y2(1)=-0.5
       u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
       endif
       do 11 i=2,n-1
       sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
       p=sig*y2(i-1)+2.
       y2(i)=(sig-1.)/p
       u(i)=(6.*((y(i+1)-y(i))/(x(i+
     *1)-x(i))-(y(i)-y(i-1))/(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*
     *u(i-1))/p
   11 continue
       if (ypn.gt..99e30) then
       qn=0.
       un=0.
       else
       qn=0.5
       un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
       endif
       y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.)
       do 12 k=n-1,1,-1
       y2(k)=y2(k)*y2(k+1)+u(k)
   12 continue
       return
       END
c ======================================================================
c
c                     Last modified: 18.03.2013              
c                                 
c                 E. Vitkovska (atavaha@hotmail.com)
c                                
c=======================================================================
