C****************************************************

          SUBROUTINE CPU(TIME)

C****************************************************

          real tarray(2)
          double precision time

C     RETURNS THE ELAPSED TOTAL CPU TIME IN MSECS

          time=etime(tarray)
          RETURN
          END

C************************************************

      SUBROUTINE INDEXR(N,ARRIN,INDX)

C************************************************

      DIMENSION INDX(*)
      REAL ARRIN(*)
      DO 11 J=1,N
        INDX(J)=J
11    CONTINUE
      L=N/2+1
      IR=N
10    CONTINUE
        IF(L.GT.1)THEN
          L=L-1
          INDXT=INDX(L)
          Q=ARRIN(INDXT)
        ELSE
          INDXT=INDX(IR)
          Q=ARRIN(INDXT)
          INDX(IR)=INDX(1)
          IR=IR-1
          IF(IR.EQ.1)THEN
            INDX(1)=INDXT
            RETURN
          ENDIF
        ENDIF
        I=L
        J=L+L
20      IF(J.LE.IR)THEN
          IF(J.LT.IR)THEN
            IF(ARRIN(INDX(J)).LT.ARRIN(INDX(J+1)))J=J+1
          ENDIF
          IF(Q.LT.ARRIN(INDX(J)))THEN
            INDX(I)=INDX(J)
            I=J
            J=J+J
          ELSE
            J=IR+1
          ENDIF
        GO TO 20
        ENDIF
        INDX(I)=INDXT
      GO TO 10
      END


C*****************************************************
      FUNCTION RAN1(IDUM)
      DIMENSION R(97)
      PARAMETER (M1=259200,IA1=7141,IC1=54773,RM1=3.8580247E-6)
      PARAMETER (M2=134456,IA2=8121,IC2=28411,RM2=7.4373773E-6)
      PARAMETER (M3=243000,IA3=4561,IC3=51349)
c     DATA IFF /0/
      save
c     IF (IDUM.LT.0.OR.IFF.EQ.0) THEN
      IF (IDUM.LT.0) THEN
       write(*,*) 'idum-seed',idum
        IFF=1
        IX1=MOD(IC1-IDUM,M1)
        IX1=MOD(IA1*IX1+IC1,M1)
        IX2=MOD(IX1,M2)
        IX1=MOD(IA1*IX1+IC1,M1)
        IX3=MOD(IX1,M3)
        DO 11 J=1,97
          IX1=MOD(IA1*IX1+IC1,M1)
          IX2=MOD(IA2*IX2+IC2,M2)
          R(J)=(FLOAT(IX1)+FLOAT(IX2)*RM2)*RM1
11      CONTINUE
        IDUM=1
      ENDIF
      IX1=MOD(IA1*IX1+IC1,M1)
      IX2=MOD(IA2*IX2+IC2,M2)
      IX3=MOD(IA3*IX3+IC3,M3)
      J=1+(97*IX3)/M3
c      IF(J.GT.97.OR.J.LT.1) PAUSE
      RAN1=R(J)
      R(J)=(FLOAT(IX1)+FLOAT(IX2)*RM2)*RM1

c      write(99,*) 'rnd1: ',idum,ran1
      RETURN
      END

C spline-programs from numerical recipes

c*******************************************************************
      SUBROUTINE spline(x,y,n,yp1,ypn,y2)
c*******************************************************************

c     tabulated function y=f(x) with n elements 
c     yp1,ypn = derivatives at the end points
c     y2 = second derivative of f
c	call spline(xxtab,atab1,6,1d30,1d30,datab1)
c              call splint(xxtab,atab1,datab1,6,alpha,aa1)

      INTEGER n,NMAX
      double precision yp1,ypn,x(n),y(n),y2(n)
      PARAMETER (NMAX=1500)
      INTEGER i,k
      double precision p,qn,sig,un,u(NMAX)
      
      
      write(6,*) 'spline: n=',n      
      write(6,*) 'x:'
      write(6,*) (x(i),i=1,n)
      write(6,*) 'y:'
      write(6,*) (y(i),i=1,n)
      write(6,*) 'yp1,yp2',yp1,ypn
      
      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) Copr. 1986-92 Numerical Recipes Software 5.

c*******************************************************************
      SUBROUTINE splint(xa,ya,y2a,n,x,y)
c*******************************************************************

c     tabulated function ya=f(xa) with n elements 
c     y2a = second derivative of f
c     y = interpolated value of f(x)

      INTEGER n
      double precision x,y,xa(n),y2a(n),ya(n)
      INTEGER k,khi,klo
      double precision a,b,h
      klo=1
      khi=n
1     if (khi-klo.gt.1) then
         k=(khi+klo)/2
         if(xa(k).gt.x)then
            khi=k
         else
            klo=k
         endif
         goto 1
      endif
      h=xa(khi)-xa(klo)
c     if (h.eq.0.) pause 'bad xa input in splint'
      a=(xa(khi)-x)/h
      b=(x-xa(klo))/h
      y=a*ya(klo)+b*ya(khi)+((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*
     *     (h**2)/6.
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software 5.