c#######################################################################
c#
c#                Cadez Integrals
c#
c#######################################################################
c
c#######################################################################
      function admsurfint(Psi,Mask,Xi,dphi,JXzc,JYzf,RHOyf,CJACyf,
     .                            ext,int,XIsing,alpha)
      implicit none
c-----------------------------------------------------------------------
      integer ext(6),int(6)
      real*8 admsurfint
      real*8 Psi(ext(1):ext(4),ext(2):ext(5),ext(3):ext(6))
      real*8 Mask(ext(1):ext(4),ext(2):ext(5),ext(3):ext(6))
      real*8 Xi(ext(2):ext(5))
      real*8 dphi
      real*8 JXzc(ext(2):ext(5))
      real*8 JYzf(ext(3):ext(6))
      real*8 RHOyf(ext(2):ext(5),ext(3):ext(6))
      real*8 CJACyf(ext(2):ext(5),ext(3):ext(6))
      real*8 XIsing,alpha
c-----------------------------------------------------------------------
      real*8 radfac,tmp,tmp1
      integer k,i,j,jf
      real*8 ZERO,ONE,TWO,FOUR,PI,TWOPI
      parameter(ZERO = 0.D0,ONE = 1.D0,TWO = 2.D0,FOUR=4.D0)
      parameter(PI = 3.1415926535897931160D0, TWOPI = 2.D0*PI)
c
c-----------------------------------------------------------------------
      if(Xi(int(2)).lt.XIsing) then
        radfac = ONE
      else
        radfac = alpha
      end if
c
c-----------------------------------------------------------------------
      admsurfint = ZERO
      j = int(3)
      jf = j-1
      do i=int(2),int(5)
        tmp = sqrt(CJACyf(i,jf))
        tmp = (FOUR*tmp*JYzf(jf) - radfac)*JXzc(i)
        tmp = dphi*RHOyf(i,jf)*TWO*JYzf(jf)*radfac/tmp
        do k=int(1),int(4)
          admsurfint = admsurfint + Mask(k,i,j)*Psi(k,i,j)*tmp
        end do
      end do
      admsurfint = admsurfint/TWOPI
c
      return
      end
c
c#######################################################################
c
c#######################################################################
      function admvolint(Psi,KK,
     .                   Mask,Volume,ext,int)
      implicit none
c-----------------------------------------------------------------------
      integer ext(6),int(6)
      real*8 admvolint
      real*8 Psi(ext(1):ext(4),ext(2):ext(5),ext(3):ext(6))
      real*8 KK(ext(1):ext(4),ext(2):ext(5),ext(3):ext(6))
      real*8 Mask(ext(1):ext(4),ext(2):ext(5),ext(3):ext(6))
      real*8 Volume(ext(2):ext(5),ext(3):ext(6))
c-----------------------------------------------------------------------
      integer k,i,j
      real*8 ZERO,F1o8,THREE,PI,TWOPI
      parameter(ZERO = 0.D0,F1o8 = 1.D0/8.D0,THREE = 3.D0)
      parameter(PI = 3.1415926535897931160D0, TWOPI = 2.D0*PI)
c
c-----------------------------------------------------------------------
      admvolint = ZERO
      do j=int(3),int(6)
        do i=int(2),int(5)
          do k=int(1),int(4)
            admvolint = admvolint +
     .        Mask(k,i,j)*Volume(i,j)*(Psi(k,i,j)**(-7))*KK(k,i,j)
          end do
        end do
      end do
      admvolint = F1o8*admvolint/TWOPI
c
      return
      end
c
c#######################################################################
c
c#######################################################################
      function xdipsurfint(Psi,Mask,Phi,Xi,dphi,JXzc,JYzf,
     .                     RHOyf,Zyf,CJACyf,ext,int,XIsing,alpha)
      implicit none
c-----------------------------------------------------------------------
      integer ext(6),int(6)
      real*8 xdipsurfint
      real*8 Psi(ext(1):ext(4),ext(2):ext(5),ext(3):ext(6))
      real*8 Mask(ext(1):ext(4),ext(2):ext(5),ext(3):ext(6))
      real*8 Phi(ext(1):ext(4))
      real*8 Xi(ext(2):ext(5))
      real*8 dphi
      real*8 JXzc(ext(2):ext(5))
      real*8 JYzf(ext(3):ext(6))
      real*8 RHOyf(ext(2):ext(5),ext(3):ext(6))
      real*8 Zyf(ext(2):ext(5),ext(3):ext(6))
      real*8 CJACyf(ext(2):ext(5),ext(3):ext(6))
      real*8 XIsing,alpha
c-----------------------------------------------------------------------
      real*8 radfac,tmpa,tmpc,tmp,zc,rhoc,dedr,dedz
      integer k,i,j,jf
      real*8 ZERO,ONE,TWO,FOUR,PI,TWOPI
      parameter(ZERO = 0.D0,ONE = 1.D0,TWO = 2.D0,FOUR=4.D0)
      parameter(PI = 3.1415926535897931160D0, TWOPI = 2.D0*PI)
c
c-----------------------------------------------------------------------
      if(Xi(int(2)).lt.XIsing) then
        radfac = ONE
      else
        radfac = alpha
      end if
c
c-----------------------------------------------------------------------
      xdipsurfint = ZERO
      j = int(3)
      jf = j-1
      do i=int(2),int(5)
        zc = Zyf(i,jf)
        rhoc = RHOyf(i,jf)
        call diffeta1(zc,rhoc,dedr,dedz)
        tmpc= CJACyf(i,jf)
        tmp = sqrt(tmpc)
        tmpc = ONE/tmp
        tmp = FOUR*tmp*JYzf(jf) - radfac
        tmpa = dphi*rhoc/JXzc(i)
        tmp = tmpa*(radfac*rhoc + TWO*dedr*tmpc)*TWO*JYzf(jf)/tmp
        tmpc = tmpa*dedr*tmpc*tmpc
        do k = int(1),int(4)
          xdipsurfint = xdipsurfint +
     .                        Mask(k,i,j)*(tmp*Psi(k,i,j) - tmpc)*
     .                        cos(Phi(k))
        end do
      end do
      xdipsurfint = xdipsurfint/TWOPI
c
      return
      end
c
c#######################################################################
c
c#######################################################################
      function ydipsurfint(Psi,Mask,Phi,Xi,dphi,JXzc,JYzf,
     .                     RHOyf,Zyf,CJACyf,ext,int,XIsing,alpha)
      implicit none
c-----------------------------------------------------------------------
      integer ext(6),int(6)
      real*8 ydipsurfint
      real*8 Psi(ext(1):ext(4),ext(2):ext(5),ext(3):ext(6))
      real*8 Mask(ext(1):ext(4),ext(2):ext(5),ext(3):ext(6))
      real*8 Phi(ext(1):ext(4))
      real*8 Xi(ext(2):ext(5))
      real*8 dphi
      real*8 JXzc(ext(2):ext(5))
      real*8 JYzf(ext(3):ext(6))
      real*8 RHOyf(ext(2):ext(5),ext(3):ext(6))
      real*8 Zyf(ext(2):ext(5),ext(3):ext(6))
      real*8 CJACyf(ext(2):ext(5),ext(3):ext(6))
      real*8 XIsing,alpha
c-----------------------------------------------------------------------
      real*8 radfac,tmpa,tmpc,tmp,zc,rhoc,dedr,dedz
      integer k,i,j,jf
      real*8 ZERO,ONE,TWO,FOUR,PI,TWOPI
      parameter(ZERO = 0.D0,ONE = 1.D0,TWO = 2.D0,FOUR=4.D0)
      parameter(PI = 3.1415926535897931160D0, TWOPI = 2.D0*PI)
c
c-----------------------------------------------------------------------
      if(Xi(int(2)).lt.XIsing) then
        radfac = ONE
      else
        radfac = alpha
      end if
c
c-----------------------------------------------------------------------
      ydipsurfint = ZERO
      j = int(3)
      jf = j-1
      do i=int(2),int(5)
        zc = Zyf(i,jf)
        rhoc = RHOyf(i,jf)
        call diffeta1(zc,rhoc,dedr,dedz)
        tmpc= CJACyf(i,jf)
        tmp = sqrt(tmpc)
        tmpc = ONE/tmp
        tmp = FOUR*tmp*JYzf(jf) - radfac
        tmpa = dphi*rhoc/JXzc(i)
        tmp = tmpa*(radfac*rhoc + TWO*dedr*tmpc)*TWO*JYzf(jf)/tmp
        tmpc = tmpa*dedr*tmpc*tmpc
        do k = int(1),int(4)
          ydipsurfint = ydipsurfint +
     .                        Mask(k,i,j)*(tmp*Psi(k,i,j) - tmpc)*
     .                        sin(Phi(k))
        end do
      end do
      ydipsurfint = ydipsurfint/TWOPI
c
      return
      end
c
c#######################################################################
c
c#######################################################################
      function zdipsurfint(Psi,Mask,Xi,dphi,JXzc,JYzf,RHOyf,Zyf,CJACyf,
     .                            ext,int,XIsing,alpha)
      implicit none
c-----------------------------------------------------------------------
      integer ext(6),int(6)
      real*8 zdipsurfint
      real*8 Psi(ext(1):ext(4),ext(2):ext(5),ext(3):ext(6))
      real*8 Mask(ext(1):ext(4),ext(2):ext(5),ext(3):ext(6))
      real*8 Xi(ext(2):ext(5))
      real*8 dphi
      real*8 JXzc(ext(2):ext(5))
      real*8 JYzf(ext(3):ext(6))
      real*8 RHOyf(ext(2):ext(5),ext(3):ext(6))
      real*8 Zyf(ext(2):ext(5),ext(3):ext(6))
      real*8 CJACyf(ext(2):ext(5),ext(3):ext(6))
      real*8 XIsing,alpha
c-----------------------------------------------------------------------
      real*8 radfac,tmpa,tmpc,tmp,zc,rhoc,dedr,dedz
      integer k,i,j,jf
      real*8 ZERO,ONE,TWO,FOUR,PI,TWOPI
      parameter(ZERO = 0.D0,ONE = 1.D0,TWO = 2.D0,FOUR=4.D0)
      parameter(PI = 3.1415926535897931160D0, TWOPI = 2.D0*PI)
c
c-----------------------------------------------------------------------
      if(Xi(int(2)).lt.XIsing) then
        radfac = ONE
      else
        radfac = alpha
      end if
c
c-----------------------------------------------------------------------
      zdipsurfint = ZERO
      j = int(3)
      jf = j-1
      do i=int(2),int(5)
        zc = Zyf(i,jf)
        rhoc = RHOyf(i,jf)
        call diffeta1(zc,rhoc,dedr,dedz)
        tmpc= CJACyf(i,jf)
        tmp = sqrt(tmpc)
        tmpc = ONE/tmp
        tmp = FOUR*tmp*JYzf(jf) - radfac
        tmpa = dphi*rhoc/JXzc(i)
        tmp = tmpa*(radfac*zc + TWO*dedz*tmpc)*TWO*JYzf(jf)/tmp
        tmpc = tmpa*dedz*tmpc*tmpc
        do k = int(1),int(4)
          zdipsurfint = zdipsurfint +
     .                        Mask(k,i,j)*(tmp*Psi(k,i,j) - tmpc)
        end do
      end do
      zdipsurfint = zdipsurfint/TWOPI
c
      return
      end
c
c#######################################################################
c
c#######################################################################
      function xdipvolint(Psi,KK,
     .                    Mask,Phi,RHOzc,Volume,ext,int)
      implicit none
c-----------------------------------------------------------------------
      integer ext(6),int(6)
      real*8 xdipvolint
      real*8 Psi(ext(1):ext(4),ext(2):ext(5),ext(3):ext(6))
      real*8 KK(ext(1):ext(4),ext(2):ext(5),ext(3):ext(6))
      real*8 Mask(ext(1):ext(4),ext(2):ext(5),ext(3):ext(6))
      real*8 Phi(ext(1):ext(4))
      real*8 RHOzc(ext(2):ext(5),ext(3):ext(6))
      real*8 Volume(ext(2):ext(5),ext(3):ext(6))
c-----------------------------------------------------------------------
      real*8 tmp
      integer k,i,j
      real*8 ZERO,F1o8,THREE,PI,TWOPI
      parameter(ZERO = 0.D0,F1o8 = 1.D0/8.D0,THREE = 3.D0)
      parameter(PI = 3.1415926535897931160D0, TWOPI = 2.D0*PI)
c
c-----------------------------------------------------------------------
      xdipvolint = ZERO
      do j=int(3),int(6)
        do i=int(2),int(5)
          tmp = RHOzc(i,j)*Volume(i,j)
          do k=int(1),int(4)
            xdipvolint = xdipvolint +
     .                     Mask(k,i,j)*tmp*(Psi(k,i,j)**(-7))*KK(k,i,j)*
     .                     cos(Phi(k))
          end do
        end do
      end do
      xdipvolint = F1o8*xdipvolint/TWOPI
c
      return
      end
c
c#######################################################################
c
c#######################################################################
      function ydipvolint(Psi,KK,
     .                    Mask,Phi,RHOzc,Volume,ext,int)
      implicit none
c-----------------------------------------------------------------------
      integer ext(6),int(6)
      real*8 ydipvolint
      real*8 Psi(ext(1):ext(4),ext(2):ext(5),ext(3):ext(6))
      real*8 KK(ext(1):ext(4),ext(2):ext(5),ext(3):ext(6))
      real*8 Mask(ext(1):ext(4),ext(2):ext(5),ext(3):ext(6))
      real*8 Phi(ext(1):ext(4))
      real*8 RHOzc(ext(2):ext(5),ext(3):ext(6))
      real*8 Volume(ext(2):ext(5),ext(3):ext(6))
c-----------------------------------------------------------------------
      real*8 tmp
      integer k,i,j
      real*8 ZERO,F1o8,THREE,PI,TWOPI
      parameter(ZERO = 0.D0,F1o8 = 1.D0/8.D0,THREE = 3.D0)
      parameter(PI = 3.1415926535897931160D0, TWOPI = 2.D0*PI)
c
c-----------------------------------------------------------------------
      ydipvolint = ZERO
      do j=int(3),int(6)
        do i=int(2),int(5)
          tmp = RHOzc(i,j)*Volume(i,j)
          do k=int(1),int(4)
            ydipvolint = ydipvolint +
     .                     Mask(k,i,j)*tmp*(Psi(k,i,j)**(-7))*KK(k,i,j)*
     .                     sin(Phi(k))
          end do
        end do
      end do
      ydipvolint = F1o8*ydipvolint/TWOPI
c
      return
      end
c
c#######################################################################
c
c#######################################################################
      function zdipvolint(Psi,KK,
     .                    Mask,Zzc,Volume,ext,int)
      implicit none
c-----------------------------------------------------------------------
      integer ext(6),int(6)
      real*8 zdipvolint
      real*8 Psi(ext(1):ext(4),ext(2):ext(5),ext(3):ext(6))
      real*8 KK(ext(1):ext(4),ext(2):ext(5),ext(3):ext(6))
      real*8 Mask(ext(1):ext(4),ext(2):ext(5),ext(3):ext(6))
      real*8 Zzc(ext(2):ext(5),ext(3):ext(6))
      real*8 Volume(ext(2):ext(5),ext(3):ext(6))
c-----------------------------------------------------------------------
      real*8 tmp
      integer k,i,j
      real*8 ZERO,F1o8,THREE,PI,TWOPI
      parameter(ZERO = 0.D0,F1o8 = 1.D0/8.D0,THREE = 3.D0)
      parameter(PI = 3.1415926535897931160D0, TWOPI = 2.D0*PI)
c
c-----------------------------------------------------------------------
      zdipvolint = ZERO
      do j=int(3),int(6)
        do i=int(2),int(5)
          tmp = Zzc(i,j)*Volume(i,j)
          do k=int(1),int(4)
            zdipvolint = zdipvolint +
     .                     Mask(k,i,j)*tmp*(Psi(k,i,j)**(-7))*KK(k,i,j)
          end do
        end do
      end do
      zdipvolint = F1o8*zdipvolint/TWOPI
c
      return
      end
c
c#######################################################################
c
c#######################################################################
      function surfarea(Psi,Mask,Xi,dphi,JXzc,JYzf,RHOyf,CJACyf,
     .                            ext,int,XIsing,alpha)
      implicit none
c-----------------------------------------------------------------------
      integer ext(6),int(6)
      real*8 surfarea
      real*8 Psi(ext(1):ext(4),ext(2):ext(5),ext(3):ext(6))
      real*8 Mask(ext(1):ext(4),ext(2):ext(5),ext(3):ext(6))
      real*8 Xi(ext(2):ext(5))
      real*8 dphi
      real*8 JXzc(ext(2):ext(5))
      real*8 JYzf(ext(3):ext(6))
      real*8 RHOyf(ext(2):ext(5),ext(3):ext(6))
      real*8 CJACyf(ext(2):ext(5),ext(3):ext(6))
      real*8 XIsing,alpha
c-----------------------------------------------------------------------
      real*8 radfac,tmp,tmp1,tmppsi4
      integer k,i,j,jf
      real*8 ZERO,ONE,TWO,FOUR
      parameter(ZERO = 0.D0,ONE = 1.D0,TWO = 2.D0,FOUR=4.D0)
c
c-----------------------------------------------------------------------
      if(Xi(int(2)).lt.XIsing) then
        radfac = ONE
      else
        radfac = alpha
      end if
c
c-----------------------------------------------------------------------
      surfarea = ZERO
      j = int(3)
      jf = j-1
      do i=int(2),int(5)
        tmppsi4 = ZERO
        do k=int(1),int(4)
          tmppsi4 = tmppsi4 + Mask(k,i,j)*(Psi(k,i,j)**4)
        end do
        tmp1 = sqrt(CJACyf(i,jf))
        tmp = FOUR*tmp1*JYzf(jf)
        tmp = tmp/(tmp - radfac)
        surfarea = surfarea +
     .      tmppsi4*(tmp**4)*RHOyf(i,jf)/(tmp1*JXzc(i))
      end do
      surfarea = surfarea*dphi
c
      return
      end
c
c#######################################################################
c
c#######################################################################
      function propsep(psiaxis,zaxis,cnt,n,Psi,Mask,Eta,Xi,JYzf,CJACyf,
     .                 ext,fil,int,Zetap,Zetam,XIsing,alpha)
      implicit none
c-----------------------------------------------------------------------
      integer n,ext(6),fil(6),int(6)
      integer propsep
      real*8 psiaxis(n),zaxis(n)
      integer cnt(n)
      real*8 Psi(ext(1):ext(4),ext(2):ext(5),ext(3):ext(6))
      real*8 Mask(ext(1):ext(4),ext(2):ext(5),ext(3):ext(6))
      real*8 Eta(ext(3):ext(6))
      real*8 Xi(ext(2):ext(4))
      real*8 JYzf(ext(3):ext(6))
      real*8 CJACyf(ext(2):ext(5),ext(3):ext(6))
      real*8 Zetap,Zetam,XIsing,alpha
c-----------------------------------------------------------------------
      integer k,i1,i2,j,region,err,chi2cyl,ptmp
      real*8 radfac,zsurf,rhodumb,tmp1,tmp2,ztmp
      real*8 ZERO,ONE,FOUR,F1o8,F9o8,CADACC
      parameter(ZERO = 0.D0,ONE = 1.D0,F1o8 = 1.D0/8.D0,
     .          F9o8 = 9.D0*F1o8,FOUR = 4.D0)
      parameter(CADACC = 1.D-10)
c
c-----------------------------------------------------------------------
      propsep = 0
      i1 = 0
      i2 = 0
      if(Xi(int(2)).lt.XIsing) then
        region = 0
        radfac = ONE
        zsurf = Zetap - ONE
        if(int(5).eq.fil(5)) then
          i1 = int(5)
          i2 = i1 - 1
        end if
      else
        region = 1
        radfac = alpha
        zsurf = Zetam + ONE/alpha
        if(int(2).eq.fil(2)) then
          i1 = int(2)
          i2 = i1 + 1
        end if
      end if
      if(i1.ne.i2) then
c NOTE!!!  Assumes INVERSION SYMMETRIC boundary conditions on throats
        j = int(3)
        tmp1 = FOUR*JYzf(j-1)
        tmp2 = tmp1*sqrt(CJACyf(i2,j-1))
        tmp1 = tmp1*sqrt(CJACyf(i1,j-1))
        tmp1 = tmp1/(tmp1 - radfac)
        tmp2 = tmp2/(tmp2 - radfac)
        psiaxis(1) = ZERO
        cnt(1) = 0
        do k=int(1),int(4)
          if((int(3).eq.fil(3)).and.(Mask(k,i1,j).eq.ONE)) then
            propsep = 1
            psiaxis(1) = psiaxis(1) +
     .                  (F9o8*Psi(k,i1,j)*tmp1 - F1o8*Psi(k,i2,j)*tmp2)
            zaxis(1) = zsurf
            cnt(1) = cnt(1) + 1
          end if
        end do
        do j = int(3),int(6)
          err = chi2cyl(Eta(j),XIsing,ztmp,
     .                  rhodumb,CADACC,region)
          ptmp = propsep
          psiaxis(ptmp+1) = ZERO
          cnt(ptmp+1) = 0
          do k = int(1),int(4)
            if(Mask(k,i1,j).eq.ONE) then
              ptmp = propsep + 1
              psiaxis(ptmp) = psiaxis(ptmp) +
     .               F9o8*Psi(k,i1,j) - F1o8*Psi(k,i2,j)
              zaxis(ptmp) = ztmp
              cnt(ptmp) = cnt(ptmp) + 1
            end if
          end do
          propsep = ptmp
        end do
      end if
c
      return
      end
c
c#######################################################################
c
c#######################################################################
      SUBROUTINE spline(X,Y,N,YP1,YPN,Y2)
      implicit none
c-----------------------------------------------------------------------
      integer N
      real*8 X(N),Y(N),YP1,YPN,Y2(N)
c-----------------------------------------------------------------------
      integer I,K
      real*8 SIG,P,QN,UN
      real*8 ZERO,HALF,ONE,TWO,THREE,SIX,BIG
      integer NMAX
      PARAMETER (ZERO=0.D0,HALF=0.5D0,
     .           ONE=1.D0,TWO=2.D0,
     .           THREE=3.D0,SIX=6.D0,
     .           BIG=0.99D30,NMAX=500)
      real*8 U(NMAX)
c
c-----------------------------------------------------------------------
      IF (YP1.GT.BIG) THEN
        Y2(1)=ZERO
        U(1)=ZERO
      ELSE
        Y2(1)=-HALF
        U(1)=(THREE/(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)+TWO
        Y2(I)=(SIG-ONE)/P
        U(I)=(SIX*((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.BIG) THEN
        QN=ZERO
        UN=ZERO
      ELSE
        QN=HALF
        UN=(THREE/(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)+ONE)
      DO 12 K=N-1,1,-1
        Y2(K)=Y2(K)*Y2(K+1)+U(K)
12    CONTINUE
      RETURN
      END
c
c#######################################################################
c
c#######################################################################
      function splint(x,y,y2,n)
      implicit none
c-----------------------------------------------------------------------
      integer n
      real*8 splint
      real*8 x(n),y(n),y2(n)
c-----------------------------------------------------------------------
      integer l
      real*8 dx,t1,t2,t3,t4,t5
      real*8 ZERO,F1o7560,F16,F31,F294,F336,F2520
      parameter(ZERO = 0.D0,F1o7560 = 1.D0/7560.D0,F16 = 16.D0,
     .          F31 = 31.D0,F294 = 294.D0,F336 = 336.D0,F2520 = 2520.D0)
c
c-----------------------------------------------------------------------
      splint = ZERO
      do 70 l=1,n-1
        dx = x(l+1) - x(l)
        t1 = y(l)**2 + y(l)*y(l+1) + y(l+1)**2
        t2 = y(l)*y2(l) + y(l+1)*y2(l+1)
        t3 = y(l)*y2(l+1) + y(l+1)*y2(l)
        t4 = y2(l)*y2(l+1)
        t5 = y2(l)**2 + y2(l+1)**2
        splint = splint + (F2520*t1 - dx*dx*(F336*t2 + F294*t3) +
     .                                    (dx**4)*(F31*t4 + F16*t5))*dx
70    continue
      splint = splint*F1o7560
c
      return
      end
c
c#######################################################################
