speclib.f

Go to the documentation of this file.

C==============================================================================
C
C     LIBRARY ROUTINES FOR SPECTRAL METHODS
C
C     March 1989
C
C     For questions, comments or suggestions, please contact:
C
C     Einar Malvin Ronquist
C     Room 3-243
C     Department of Mechanical Engineering
C     Massachusetts Institute of Technology
C     77 Massachusetts Avenue
C     Cambridge, MA 0299
C     U.S.A.
C
C------------------------------------------------------------------------------
C
C     ABBRIVIATIONS:
C
C     M   - Set of mesh points
C     Z   - Set of collocation/quadrature points
C     W   - Set of quadrature weights
C     H   - Lagrangian interpolant
C     D   - Derivative operator
C     I   - Interpolation operator
C     GL  - Gauss Legendre
C     GLL - Gauss-Lobatto Legendre
C     GJ  - Gauss Jacobi
C     GLJ - Gauss-Lobatto Jacobi
C
C
C     MAIN ROUTINES:
C
C     Points and weights:
C
C     ZWGL      Compute Gauss Legendre points and weights
C     ZWGLL     Compute Gauss-Lobatto Legendre points and weights
C     ZWGJ      Compute Gauss Jacobi points and weights (general)
C     ZWGLJ     Compute Gauss-Lobatto Jacobi points and weights (general)
C
C     Lagrangian interpolants:
C
C     HGL       Compute Gauss Legendre Lagrangian interpolant
C     HGLL      Compute Gauss-Lobatto Legendre Lagrangian interpolant
C     HGJ       Compute Gauss Jacobi Lagrangian interpolant (general)
C     HGLJ      Compute Gauss-Lobatto Jacobi Lagrangian interpolant (general)
C
C     Derivative operators:
C
C     DGLL      Compute Gauss-Lobatto Legendre derivative matrix
C     DGLLGL    Compute derivative matrix for a staggered mesh (GLL->GL)
C     DGJ       Compute Gauss Jacobi derivative matrix (general)
C     DGLJ      Compute Gauss-Lobatto Jacobi derivative matrix (general)
C     DGLJGJ    Compute derivative matrix for a staggered mesh (GLJ->GJ) (general)
C
C     Interpolation operators:
C
C     IGLM      Compute interpolation operator GL  -> M
C     IGLLM     Compute interpolation operator GLL -> M
C     IGJM      Compute interpolation operator GJ  -> M  (general)
C     IGLJM     Compute interpolation operator GLJ -> M  (general)
C
C     Other:
C
C     PNLEG     Compute Legendre polynomial of degree N
C     PNDLEG    Compute derivative of Legendre polynomial of degree N
C
C     Comments:
C
C     Note that many of the above routines exist in both single and
C     double precision. If the name of the single precision routine is
C     SUB, the double precision version is called SUBD. In most cases
C     all the "low-level" arithmetic is done in double precision, even
C     for the single precsion versions.
C
C     Useful references:
C
C [1] Gabor Szego: Orthogonal Polynomials, American Mathematical Society,
C     Providence, Rhode Island, 1939.
C [2] Abramowitz & Stegun: Handbook of Mathematical Functions,
C     Dover, New York, 1972.
C [3] Canuto, Hussaini, Quarteroni & Zang: Spectral Methods in Fluid
C     Dynamics, Springer-Verlag, 1988.
C
C
C==============================================================================
C
C--------------------------------------------------------------------
      SUBROUTINE zwgl (Z,W,NP)
C--------------------------------------------------------------------
C
C     Generate NP Gauss Legendre points (Z) and weights (W)
C     associated with Jacobi polynomial P(N)(alpha=0,beta=0).
C     The polynomial degree N=NP-1.
C     Z and W are in single precision, but all the arithmetic
C     operations are done in double precision.
C
C--------------------------------------------------------------------
      REAL Z(1),W(1)
      alpha = 0.
      beta  = 0.
      CALL zwgj (z,w,np,alpha,beta)
      RETURN
      END
C
      SUBROUTINE zwgll (Z,W,NP)
C--------------------------------------------------------------------
C
C     Generate NP Gauss-Lobatto Legendre points (Z) and weights (W)
C     associated with Jacobi polynomial P(N)(alpha=0,beta=0).
C     The polynomial degree N=NP-1.
C     Z and W are in single precision, but all the arithmetic
C     operations are done in double precision.
C
C--------------------------------------------------------------------
      REAL Z(1),W(1)
      alpha = 0.
      beta  = 0.
      CALL zwglj (z,w,np,alpha,beta)
      RETURN
      END
C
      SUBROUTINE zwgj (Z,W,NP,ALPHA,BETA)
C--------------------------------------------------------------------
C
C     Generate NP GAUSS JACOBI points (Z) and weights (W)
C     associated with Jacobi polynomial P(N)(alpha>-1,beta>-1).
C     The polynomial degree N=NP-1.
C     Single precision version.
C
C--------------------------------------------------------------------
      parameter(nmax=84)
      parameter(lzd = nmax)
      real*8  zd(lzd),wd(lzd),aphad,betad
      REAL Z(1),W(1),ALPHA,BETA
C
      npmax = lzd
      IF (np.GT.npmax) THEN
         WRITE (6,*) 'Too large polynomial degree in ZWGJ'
         WRITE (6,*) 'Maximum polynomial degree is',nmax
         WRITE (6,*) 'Here NP=',np
         call exitt
      ENDIF
      alphad = alpha
      betad  = beta
      CALL zwgjd (zd,wd,np,alphad,betad)
      DO 100 i=1,np
         z(i) = zd(i)
         w(i) = wd(i)
 100  CONTINUE
      RETURN
      END
C
      SUBROUTINE zwgjd (Z,W,NP,ALPHA,BETA)
C--------------------------------------------------------------------
C
C     Generate NP GAUSS JACOBI points (Z) and weights (W)
C     associated with Jacobi polynomial P(N)(alpha>-1,beta>-1).
C     The polynomial degree N=NP-1.
C     Double precision version.
C
C--------------------------------------------------------------------
      IMPLICIT REAL*8  (a-h,o-z)
      real*8  z(1),w(1),alpha,beta
C
      n     = np-1
      dn    = ((n))
      one   = 1.
      two   = 2.
      apb   = alpha+beta
C
      IF (np.LE.0) THEN
         WRITE (6,*) 'ZWGJD: Minimum number of Gauss points is 1',np
         call exitt
      ENDIF
      IF ((alpha.LE.-one).OR.(beta.LE.-one)) THEN
         WRITE (6,*) 'ZWGJD: Alpha and Beta must be greater than -1'
         call exitt
      ENDIF
C
      IF (np.EQ.1) THEN
         z(1) = (beta-alpha)/(apb+two)
         w(1) = gammaf(alpha+one)*gammaf(beta+one)/gammaf(apb+two)
     $          * two**(apb+one)
         RETURN
      ENDIF
C
      CALL jacg (z,np,alpha,beta)
C
      np1   = n+1
      np2   = n+2
      dnp1  = ((np1))
      dnp2  = ((np2))
      fac1  = dnp1+alpha+beta+one
      fac2  = fac1+dnp1
      fac3  = fac2+one
      fnorm = pnormj(np1,alpha,beta)
      rcoef = (fnorm*fac2*fac3)/(two*fac1*dnp2)
      DO 100 i=1,np
         CALL jacobf (p,pd,pm1,pdm1,pm2,pdm2,np2,alpha,beta,z(i))
         w(i) = -rcoef/(p*pdm1)
 100  CONTINUE
      RETURN
      END
C
      SUBROUTINE zwglj (Z,W,NP,ALPHA,BETA)
C--------------------------------------------------------------------
C
C     Generate NP GAUSS LOBATTO JACOBI points (Z) and weights (W)
C     associated with Jacobi polynomial P(N)(alpha>-1,beta>-1).
C     The polynomial degree N=NP-1.
C     Single precision version.
C
C--------------------------------------------------------------------
      parameter(nmax=84)
      parameter(lzd = nmax)
      real*8  zd(lzd),wd(lzd),alphad,betad
      REAL Z(1),W(1),ALPHA,BETA
C
      npmax = lzd
      IF (np.GT.npmax) THEN
         WRITE (6,*) 'Too large polynomial degree in ZWGLJ'
         WRITE (6,*) 'Maximum polynomial degree is',nmax
         WRITE (6,*) 'Here NP=',np
         call exitt
      ENDIF
      alphad = alpha
      betad  = beta
      CALL zwgljd (zd,wd,np,alphad,betad)
      DO 100 i=1,np
         z(i) = zd(i)
         w(i) = wd(i)
 100  CONTINUE
      RETURN
      END
C
      SUBROUTINE zwgljd (Z,W,NP,ALPHA,BETA)
C--------------------------------------------------------------------
C
C     Generate NP GAUSS LOBATTO JACOBI points (Z) and weights (W)
C     associated with Jacobi polynomial P(N)(alpha>-1,beta>-1).
C     The polynomial degree N=NP-1.
C     Double precision version.
C
C--------------------------------------------------------------------
      IMPLICIT REAL*8  (a-h,o-z)
      real*8  z(np),w(np),alpha,beta
C
      n     = np-1
      nm1   = n-1
      one   = 1.
      two   = 2.
C
      IF (np.LE.1) THEN
       WRITE (6,*) 'ZWGLJD: Minimum number of Gauss-Lobatto points is 2'
       WRITE (6,*) 'ZWGLJD: alpha,beta:',alpha,beta,np
       call exitt
      ENDIF
      IF ((alpha.LE.-one).OR.(beta.LE.-one)) THEN
         WRITE (6,*) 'ZWGLJD: Alpha and Beta must be greater than -1'
         call exitt
      ENDIF
C
      IF (nm1.GT.0) THEN
         alpg  = alpha+one
         betg  = beta+one
         CALL zwgjd (z(2),w(2),nm1,alpg,betg)
      ENDIF
      z(1)  = -one
      z(np) =  one
      DO 100  i=2,np-1
         w(i) = w(i)/(one-z(i)**2)
 100  CONTINUE
      CALL jacobf (p,pd,pm1,pdm1,pm2,pdm2,n,alpha,beta,z(1))
      w(1)  = endw1(n,alpha,beta)/(two*pd)
      CALL jacobf (p,pd,pm1,pdm1,pm2,pdm2,n,alpha,beta,z(np))
      w(np) = endw2(n,alpha,beta)/(two*pd)
C
      RETURN
      END
C
      REAL*8  FUNCTION endw1 (N,ALPHA,BETA)
      IMPLICIT REAL*8  (a-h,o-z)
      real*8  alpha,beta
      zero  = 0.
      one   = 1.
      two   = 2.
      three = 3.
      four  = 4.
      apb   = alpha+beta
      IF (n.EQ.0) THEN
         endw1 = zero
         RETURN
      ENDIF
      f1   = gammaf(alpha+two)*gammaf(beta+one)/gammaf(apb+three)
      f1   = f1*(apb+two)*two**(apb+two)/two
      IF (n.EQ.1) THEN
         endw1 = f1
         RETURN
      ENDIF
      fint1 = gammaf(alpha+two)*gammaf(beta+one)/gammaf(apb+three)
      fint1 = fint1*two**(apb+two)
      fint2 = gammaf(alpha+two)*gammaf(beta+two)/gammaf(apb+four)
      fint2 = fint2*two**(apb+three)
      f2    = (-two*(beta+two)*fint1 + (apb+four)*fint2)
     $        * (apb+three)/four
      IF (n.EQ.2) THEN
         endw1 = f2
         RETURN
      ENDIF
      DO 100 i=3,n
         di   = ((i-1))
         abn  = alpha+beta+di
         abnn = abn+di
         a1   = -(two*(di+alpha)*(di+beta))/(abn*abnn*(abnn+one))
         a2   =  (two*(alpha-beta))/(abnn*(abnn+two))
         a3   =  (two*(abn+one))/((abnn+two)*(abnn+one))
         f3   =  -(a2*f2+a1*f1)/a3
         f1   = f2
         f2   = f3
 100  CONTINUE
      endw1  = f3
      RETURN
      END
C
      REAL*8  FUNCTION endw2 (N,ALPHA,BETA)
      IMPLICIT REAL*8  (a-h,o-z)
      real*8  alpha,beta
      zero  = 0.
      one   = 1.
      two   = 2.
      three = 3.
      four  = 4.
      apb   = alpha+beta
      IF (n.EQ.0) THEN
         endw2 = zero
         RETURN
      ENDIF
      f1   = gammaf(alpha+one)*gammaf(beta+two)/gammaf(apb+three)
      f1   = f1*(apb+two)*two**(apb+two)/two
      IF (n.EQ.1) THEN
         endw2 = f1
         RETURN
      ENDIF
      fint1 = gammaf(alpha+one)*gammaf(beta+two)/gammaf(apb+three)
      fint1 = fint1*two**(apb+two)
      fint2 = gammaf(alpha+two)*gammaf(beta+two)/gammaf(apb+four)
      fint2 = fint2*two**(apb+three)
      f2    = (two*(alpha+two)*fint1 - (apb+four)*fint2)
     $        * (apb+three)/four
      IF (n.EQ.2) THEN
         endw2 = f2
         RETURN
      ENDIF
      DO 100 i=3,n
         di   = ((i-1))
         abn  = alpha+beta+di
         abnn = abn+di
         a1   =  -(two*(di+alpha)*(di+beta))/(abn*abnn*(abnn+one))
         a2   =  (two*(alpha-beta))/(abnn*(abnn+two))
         a3   =  (two*(abn+one))/((abnn+two)*(abnn+one))
         f3   =  -(a2*f2+a1*f1)/a3
         f1   = f2
         f2   = f3
 100  CONTINUE
      endw2  = f3
      RETURN
      END
C
      REAL*8  FUNCTION gammaf (X)
      IMPLICIT REAL*8  (a-h,o-z)
      real*8  x
      zero = 0.0
      half = 0.5
      one  = 1.0
      two  = 2.0
      four = 4.0
      pi   = four*atan(one)
      gammaf = one
      IF (x.EQ.-half) gammaf = -two*sqrt(pi)
      IF (x.EQ. half) gammaf =  sqrt(pi)
      IF (x.EQ. one ) gammaf =  one
      IF (x.EQ. two ) gammaf =  one
      IF (x.EQ. 1.5  ) gammaf =  sqrt(pi)/2.
      IF (x.EQ. 2.5) gammaf =  1.5*sqrt(pi)/2.
      IF (x.EQ. 3.5) gammaf =  0.5*(2.5*(1.5*sqrt(pi)))
      IF (x.EQ. 3. ) gammaf =  2.
      IF (x.EQ. 4. ) gammaf = 6.
      IF (x.EQ. 5. ) gammaf = 24.
      IF (x.EQ. 6. ) gammaf = 120.
      RETURN
      END
C
      REAL*8  FUNCTION pnormj (N,ALPHA,BETA)
      IMPLICIT REAL*8  (a-h,o-z)
      real*8  alpha,beta
      one   = 1.
      two   = 2.
      dn    = ((n))
      const = alpha+beta+one
      IF (n.LE.1) THEN
         prod   = gammaf(dn+alpha)*gammaf(dn+beta)
         prod   = prod/(gammaf(dn)*gammaf(dn+alpha+beta))
         pnormj = prod * two**const/(two*dn+const)
         RETURN
      ENDIF
      prod  = gammaf(alpha+one)*gammaf(beta+one)
      prod  = prod/(two*(one+const)*gammaf(const+one))
      prod  = prod*(one+alpha)*(two+alpha)
      prod  = prod*(one+beta)*(two+beta)
      DO 100 i=3,n
         dindx = ((i))
         frac  = (dindx+alpha)*(dindx+beta)/(dindx*(dindx+alpha+beta))
         prod  = prod*frac
 100  CONTINUE
      pnormj = prod * two**const/(two*dn+const)
      RETURN
      END
C
      SUBROUTINE jacg (XJAC,NP,ALPHA,BETA)
C--------------------------------------------------------------------
C
C     Compute NP Gauss points XJAC, which are the zeros of the
C     Jacobi polynomial J(NP) with parameters ALPHA and BETA.
C     ALPHA and BETA determines the specific type of Gauss points.
C     Examples:
C     ALPHA = BETA =  0.0  ->  Legendre points
C     ALPHA = BETA = -0.5  ->  Chebyshev points
C
C--------------------------------------------------------------------
      IMPLICIT REAL*8  (a-h,o-z)
      real*8  xjac(1)
      DATA kstop /10/
      DATA eps/1.0e-12/
      n   = np-1
      one = 1.
      dth = 4.*atan(one)/(2.*((n))+2.)
      DO 40 j=1,np
         IF (j.EQ.1) THEN
            x = cos((2.*(((j))-1.)+1.)*dth)
         ELSE
            x1 = cos((2.*(((j))-1.)+1.)*dth)
            x2 = xlast
            x  = (x1+x2)/2.
         ENDIF
         DO 30 k=1,kstop
            CALL jacobf (p,pd,pm1,pdm1,pm2,pdm2,np,alpha,beta,x)
            recsum = 0.
            jm = j-1
            DO 29 i=1,jm
               recsum = recsum+1./(x-xjac(np-i+1))
 29         CONTINUE
            delx = -p/(pd-recsum*p)
            x    = x+delx
            IF (abs(delx) .LT. eps) GOTO 31
 30      CONTINUE
 31      CONTINUE
         xjac(np-j+1) = x
         xlast        = x
 40   CONTINUE
      DO 200 i=1,np
         xmin = 2.
         DO 100 j=i,np
            IF (xjac(j).LT.xmin) THEN
               xmin = xjac(j)
               jmin = j
            ENDIF
 100     CONTINUE
         IF (jmin.NE.i) THEN
            swap = xjac(i)
            xjac(i) = xjac(jmin)
            xjac(jmin) = swap
         ENDIF
 200  CONTINUE
      RETURN
      END
C
      SUBROUTINE jacobf (POLY,PDER,POLYM1,PDERM1,POLYM2,PDERM2,
     $                   N,ALP,BET,X)
C--------------------------------------------------------------------
C
C     Computes the Jacobi polynomial (POLY) and its derivative (PDER)
C     of degree N at X.
C
C--------------------------------------------------------------------
      IMPLICIT REAL*8  (a-h,o-z)
      apb  = alp+bet
      poly = 1.
      pder = 0.
      IF (n .EQ. 0) RETURN
      polyl = poly
      pderl = pder
      poly  = (alp-bet+(apb+2.)*x)/2.
      pder  = (apb+2.)/2.
      IF (n .EQ. 1) RETURN
      DO 20 k=2,n
         dk = ((k))
         a1 = 2.*dk*(dk+apb)*(2.*dk+apb-2.)
         a2 = (2.*dk+apb-1.)*(alp**2-bet**2)
         b3 = (2.*dk+apb-2.)
         a3 = b3*(b3+1.)*(b3+2.)
         a4 = 2.*(dk+alp-1.)*(dk+bet-1.)*(2.*dk+apb)
         polyn  = ((a2+a3*x)*poly-a4*polyl)/a1
         pdern  = ((a2+a3*x)*pder-a4*pderl+a3*poly)/a1
         psave  = polyl
         pdsave = pderl
         polyl  = poly
         poly   = polyn
         pderl  = pder
         pder   = pdern
 20   CONTINUE
      polym1 = polyl
      pderm1 = pderl
      polym2 = psave
      pderm2 = pdsave
      RETURN
      END
C
      REAL function hgj (ii,z,zgj,np,alpha,beta)
C---------------------------------------------------------------------
C
C     Compute the value of the Lagrangian interpolant HGJ through
C     the NP Gauss Jacobi points ZGJ at the point Z.
C     Single precision version.
C
C---------------------------------------------------------------------
      parameter(nmax=84)
      parameter(lzd = nmax)
      real*8  zd,zgjd(lzd),hgjd,alphad,betad
      REAL z,zgj(1),alpha,beta
      npmax = lzd
      IF (np.GT.npmax) THEN
         WRITE (6,*) 'Too large polynomial degree in HGJ'
         WRITE (6,*) 'Maximum polynomial degree is',nmax
         WRITE (6,*) 'Here NP=',np
         call exitt
      ENDIF
      zd = z
      DO 100 i=1,np
         zgjd(i) = zgj(i)
 100  CONTINUE
      alphad = alpha
      betad  = beta
      hgj    = hgjd(ii,zd,zgjd,np,alphad,betad)
      RETURN
      END
C
      REAL*8  FUNCTION hgjd (II,Z,ZGJ,NP,ALPHA,BETA)
C---------------------------------------------------------------------
C
C     Compute the value of the Lagrangian interpolant HGJD through
C     the NZ Gauss-Lobatto Jacobi points ZGJ at the point Z.
C     Double precision version.
C
C---------------------------------------------------------------------
      IMPLICIT REAL*8  (a-h,o-z)
      real*8  z,zgj(1),alpha,beta
      eps = 1.e-5
      one = 1.
      zi  = zgj(ii)
      dz  = z-zi
      IF (abs(dz).LT.eps) THEN
         hgjd = one
         RETURN
      ENDIF
      CALL jacobf (pzi,pdzi,pm1,pdm1,pm2,pdm2,np,alpha,beta,zi)
      CALL jacobf (pz,pdz,pm1,pdm1,pm2,pdm2,np,alpha,beta,z)
      hgjd  = pz/(pdzi*(z-zi))
      RETURN
      END
C
      REAL function hglj (ii,z,zglj,np,alpha,beta)
C---------------------------------------------------------------------
C
C     Compute the value of the Lagrangian interpolant HGLJ through
C     the NZ Gauss-Lobatto Jacobi points ZGLJ at the point Z.
C     Single precision version.
C
C---------------------------------------------------------------------
      parameter(nmax=84)
      parameter(lzd = nmax)
      real*8  zd,zgljd(lzd),hgljd,alphad,betad
      REAL z,zglj(1),alpha,beta
      npmax = lzd
      IF (np.GT.npmax) THEN
         WRITE (6,*) 'Too large polynomial degree in HGLJ'
         WRITE (6,*) 'Maximum polynomial degree is',nmax
         WRITE (6,*) 'Here NP=',np
         call exitt
      ENDIF
      zd = z
      DO 100 i=1,np
         zgljd(i) = zglj(i)
 100  CONTINUE
      alphad = alpha
      betad  = beta
      hglj   = hgljd(ii,zd,zgljd,np,alphad,betad)
      RETURN
      END
C
      REAL*8  FUNCTION hgljd (I,Z,ZGLJ,NP,ALPHA,BETA)
C---------------------------------------------------------------------
C
C     Compute the value of the Lagrangian interpolant HGLJD through
C     the NZ Gauss-Lobatto Jacobi points ZJACL at the point Z.
C     Double precision version.
C
C---------------------------------------------------------------------
      IMPLICIT REAL*8  (a-h,o-z)
      real*8  z,zglj(1),alpha,beta
      eps = 1.e-5
      one = 1.
      zi  = zglj(i)
      dz  = z-zi
      IF (abs(dz).LT.eps) THEN
         hgljd = one
         RETURN
      ENDIF
      n      = np-1
      dn     = ((n))
      eigval = -dn*(dn+alpha+beta+one)
      CALL jacobf (pi,pdi,pm1,pdm1,pm2,pdm2,n,alpha,beta,zi)
      const  = eigval*pi+alpha*(one+zi)*pdi-beta*(one-zi)*pdi
      CALL jacobf (p,pd,pm1,pdm1,pm2,pdm2,n,alpha,beta,z)
      hgljd  = (one-z**2)*pd/(const*(z-zi))
      RETURN
      END
C
      SUBROUTINE dgj (D,DT,Z,NZ,lzd,ALPHA,BETA)
C-----------------------------------------------------------------
C
C     Compute the derivative matrix D and its transpose DT
C     associated with the Nth order Lagrangian interpolants
C     through the NZ Gauss Jacobi points Z.
C     Note: D and DT are square matrices.
C     Single precision version.
C
C-----------------------------------------------------------------
      parameter(nmax=84)
      parameter(lzdd = nmax)
      real*8  dd(lzdd,lzdd),dtd(lzdd,lzdd),zd(lzdd),alphad,betad
      REAL D(lzd,lzd),DT(lzd,lzd),Z(1),ALPHA,BETA
C
      IF (nz.LE.0) THEN
         WRITE (6,*) 'DGJ: Minimum number of Gauss points is 1'
         call exitt
      ENDIF
      IF (nz .GT. nmax) THEN
         WRITE (6,*) 'Too large polynomial degree in DGJ'
         WRITE (6,*) 'Maximum polynomial degree is',nmax
         WRITE (6,*) 'Here Nz=',nz
         call exitt
      ENDIF
      IF ((alpha.LE.-1.).OR.(beta.LE.-1.)) THEN
         WRITE (6,*) 'DGJ: Alpha and Beta must be greater than -1'
         call exitt
      ENDIF
      alphad = alpha
      betad  = beta
      DO 100 i=1,nz
         zd(i) = z(i)
 100  CONTINUE
      CALL dgjd (dd,dtd,zd,nz,lzdd,alphad,betad)
      DO 200 i=1,nz
      DO 200 j=1,nz
         d(i,j)  = dd(i,j)
         dt(i,j) = dtd(i,j)
 200  CONTINUE
      RETURN
      END
C
      SUBROUTINE dgjd (D,DT,Z,NZ,lzd,ALPHA,BETA)
C-----------------------------------------------------------------
C
C     Compute the derivative matrix D and its transpose DT
C     associated with the Nth order Lagrangian interpolants
C     through the NZ Gauss Jacobi points Z.
C     Note: D and DT are square matrices.
C     Double precision version.
C
C-----------------------------------------------------------------
      IMPLICIT REAL*8  (a-h,o-z)
      real*8  d(lzd,lzd),dt(lzd,lzd),z(1),alpha,beta
      n    = nz-1
      dn   = ((n))
      one  = 1.
      two  = 2.
C
      IF (nz.LE.1) THEN
       WRITE (6,*) 'DGJD: Minimum number of Gauss-Lobatto points is 2'
       call exitt
      ENDIF
      IF ((alpha.LE.-one).OR.(beta.LE.-one)) THEN
         WRITE (6,*) 'DGJD: Alpha and Beta must be greater than -1'
         call exitt
      ENDIF
C
      DO 200 i=1,nz
      DO 200 j=1,nz
         CALL jacobf (pi,pdi,pm1,pdm1,pm2,pdm2,nz,alpha,beta,z(i))
         CALL jacobf (pj,pdj,pm1,pdm1,pm2,pdm2,nz,alpha,beta,z(j))
         IF (i.NE.j) d(i,j) = pdi/(pdj*(z(i)-z(j)))
         IF (i.EQ.j) d(i,j) = ((alpha+beta+two)*z(i)+alpha-beta)/
     $                        (two*(one-z(i)**2))
         dt(j,i) = d(i,j)
 200  CONTINUE
      RETURN
      END
C
      SUBROUTINE dglj (D,DT,Z,NZ,lzd,ALPHA,BETA)
C-----------------------------------------------------------------
C
C     Compute the derivative matrix D and its transpose DT
C     associated with the Nth order Lagrangian interpolants
C     through the NZ Gauss-Lobatto Jacobi points Z.
C     Note: D and DT are square matrices.
C     Single precision version.
C
C-----------------------------------------------------------------
      parameter(nmax=84)
      parameter(lzdd = nmax)
      real*8  dd(lzdd,lzdd),dtd(lzdd,lzdd),zd(lzdd),alphad,betad
      REAL D(lzd,lzd),DT(lzd,lzd),Z(1),ALPHA,BETA
C
      IF (nz.LE.1) THEN
       WRITE (6,*) 'DGLJ: Minimum number of Gauss-Lobatto points is 2'
       call exitt
      ENDIF
      IF (nz .GT. nmax) THEN
         WRITE (6,*) 'Too large polynomial degree in DGLJ'
         WRITE (6,*) 'Maximum polynomial degree is',nmax
         WRITE (6,*) 'Here NZ=',nz
         call exitt
      ENDIF
      IF ((alpha.LE.-1.).OR.(beta.LE.-1.)) THEN
         WRITE (6,*) 'DGLJ: Alpha and Beta must be greater than -1'
         call exitt
      ENDIF
      alphad = alpha
      betad  = beta
      DO 100 i=1,nz
         zd(i) = z(i)
 100  CONTINUE
      CALL dgljd (dd,dtd,zd,nz,lzdd,alphad,betad)
      DO 200 i=1,nz
      DO 200 j=1,nz
         d(i,j)  = dd(i,j)
         dt(i,j) = dtd(i,j)
 200  CONTINUE
      RETURN
      END
C
      SUBROUTINE dgljd (D,DT,Z,NZ,lzd,ALPHA,BETA)
C-----------------------------------------------------------------
C
C     Compute the derivative matrix D and its transpose DT
C     associated with the Nth order Lagrangian interpolants
C     through the NZ Gauss-Lobatto Jacobi points Z.
C     Note: D and DT are square matrices.
C     Double precision version.
C
C-----------------------------------------------------------------
      IMPLICIT REAL*8  (a-h,o-z)
      real*8  d(lzd,lzd),dt(lzd,lzd),z(1),alpha,beta
      n    = nz-1
      dn   = ((n))
      one  = 1.
      two  = 2.
      eigval = -dn*(dn+alpha+beta+one)
C
      IF (nz.LE.1) THEN
       WRITE (6,*) 'DGLJD: Minimum number of Gauss-Lobatto points is 2'
       call exitt
      ENDIF
      IF ((alpha.LE.-one).OR.(beta.LE.-one)) THEN
         WRITE (6,*) 'DGLJD: Alpha and Beta must be greater than -1'
         call exitt
      ENDIF
C
      DO 200 i=1,nz
      DO 200 j=1,nz
         CALL jacobf (pi,pdi,pm1,pdm1,pm2,pdm2,n,alpha,beta,z(i))
         CALL jacobf (pj,pdj,pm1,pdm1,pm2,pdm2,n,alpha,beta,z(j))
         ci = eigval*pi-(beta*(one-z(i))-alpha*(one+z(i)))*pdi
         cj = eigval*pj-(beta*(one-z(j))-alpha*(one+z(j)))*pdj
         IF (i.NE.j) d(i,j) = ci/(cj*(z(i)-z(j)))
         IF ((i.EQ.j).AND.(i.NE.1).AND.(i.NE.nz))
     $   d(i,j) = (alpha*(one+z(i))-beta*(one-z(i)))/
     $            (two*(one-z(i)**2))
         IF ((i.EQ.j).AND.(i.EQ.1))
     $   d(i,j) =  (eigval+alpha)/(two*(beta+two))
         IF ((i.EQ.j).AND.(i.EQ.nz))
     $   d(i,j) = -(eigval+beta)/(two*(alpha+two))
         dt(j,i) = d(i,j)
 200  CONTINUE
      RETURN
      END
C
      SUBROUTINE dgll (D,DT,Z,NZ,lzd)
C-----------------------------------------------------------------
C
C     Compute the derivative matrix D and its transpose DT
C     associated with the Nth order Lagrangian interpolants
C     through the NZ Gauss-Lobatto Legendre points Z.
C     Note: D and DT are square matrices.
C
C-----------------------------------------------------------------
      parameter(nmax=84)
      REAL D(lzd,lzd),DT(lzd,lzd),Z(1)
      n  = nz-1
      IF (nz .GT. nmax) THEN
         WRITE (6,*) 'Subroutine DGLL'
         WRITE (6,*) 'Maximum polynomial degree =',nmax
         WRITE (6,*) 'Polynomial degree         =',nz
      ENDIF
      IF (nz .EQ. 1) THEN
         d(1,1) = 0.
         RETURN
      ENDIF
      fn = (n)
      d0 = fn*(fn+1.)/4.
      DO 200 i=1,nz
      DO 200 j=1,nz
         d(i,j) = 0.
         IF  (i.NE.j) d(i,j) = pnleg(z(i),n)/
     $                        (pnleg(z(j),n)*(z(i)-z(j)))
         IF ((i.EQ.j).AND.(i.EQ.1))  d(i,j) = -d0
         IF ((i.EQ.j).AND.(i.EQ.nz)) d(i,j) =  d0
         dt(j,i) = d(i,j)
 200  CONTINUE
      RETURN
      END
C
      REAL function hgll (i,z,zgll,nz)
C---------------------------------------------------------------------
C
C     Compute the value of the Lagrangian interpolant L through
C     the NZ Gauss-Lobatto Legendre points ZGLL at the point Z.
C
C---------------------------------------------------------------------
      REAL zgll(1)
      eps = 1.e-5
      dz = z - zgll(i)
      IF (abs(dz) .LT. eps) THEN
         hgll = 1.
         RETURN
      ENDIF
      n = nz - 1
      alfan = (n)*((n)+1.)
      hgll = - (1.-z*z)*pndleg(z,n)/
     $         (alfan*pnleg(zgll(i),n)*(z-zgll(i)))
      RETURN
      END
C
      REAL function hgl (i,z,zgl,nz)
C---------------------------------------------------------------------
C
C     Compute the value of the Lagrangian interpolant HGL through
C     the NZ Gauss Legendre points ZGL at the point Z.
C
C---------------------------------------------------------------------
      REAL zgl(1)
      eps = 1.e-5
      dz = z - zgl(i)
      IF (abs(dz) .LT. eps) THEN
         hgl = 1.
         RETURN
      ENDIF
      n = nz-1
      hgl = pnleg(z,nz)/(pndleg(zgl(i),nz)*(z-zgl(i)))
      RETURN
      END
C
      REAL function pnleg (z,n)
C---------------------------------------------------------------------
C
C     Compute the value of the Nth order Legendre polynomial at Z.
C     (Simpler than JACOBF)
C     Based on the recursion formula for the Legendre polynomials.
C
C---------------------------------------------------------------------
C
C     This next statement is to overcome the underflow bug in the i860.  
C     It can be removed at a later date.  11 Aug 1990   pff.
C
      IF(abs(z) .LT. 1.0e-25) z = 0.0
C
      p1   = 1.
      IF (n.EQ.0) THEN
         pnleg = p1
         RETURN
      ENDIF
      p2   = z
      p3   = p2
      DO 10 k = 1, n-1
         fk  = (k)
         p3  = ((2.*fk+1.)*z*p2 - fk*p1)/(fk+1.)
         p1  = p2
         p2  = p3
 10   CONTINUE
      pnleg = p3
      if (n.eq.0) pnleg = 1.
      RETURN
      END
C
      REAL function pndleg (z,n)
C----------------------------------------------------------------------
C
C     Compute the derivative of the Nth order Legendre polynomial at Z.
C     (Simpler than JACOBF)
C     Based on the recursion formula for the Legendre polynomials.
C
C----------------------------------------------------------------------
      p1   = 1.
      p2   = z
      p1d  = 0.
      p2d  = 1.
      p3d  = 1.
      DO 10 k = 1, n-1
         fk  = (k)
         p3  = ((2.*fk+1.)*z*p2 - fk*p1)/(fk+1.)
         p3d = ((2.*fk+1.)*p2 + (2.*fk+1.)*z*p2d - fk*p1d)/(fk+1.)
         p1  = p2
         p2  = p3
         p1d = p2d
         p2d = p3d
 10   CONTINUE
      pndleg = p3d
      IF (n.eq.0) pndleg = 0.
      RETURN
      END
C
      SUBROUTINE dgllgl (D,DT,ZM1,ZM2,IM12,NZM1,NZM2,ND1,ND2)
C-----------------------------------------------------------------------
C
C     Compute the (one-dimensional) derivative matrix D and its
C     transpose DT associated with taking the derivative of a variable
C     expanded on a Gauss-Lobatto Legendre mesh (M1), and evaluate its
C     derivative on a Guass Legendre mesh (M2).
C     Need the one-dimensional interpolation operator IM12
C     (see subroutine IGLLGL).
C     Note: D and DT are rectangular matrices.
C
C-----------------------------------------------------------------------
      REAL D(ND2,ND1), DT(ND1,ND2), ZM1(ND1), ZM2(ND2), IM12(ND2,ND1)
      IF (nzm1.EQ.1) THEN
        d(1,1) = 0.
        dt(1,1) = 0.
        RETURN
      ENDIF
      eps = 1.e-6
      nm1 = nzm1-1
      DO 10 ip = 1, nzm2
         DO 10 jq = 1, nzm1
            zp = zm2(ip)
            zq = zm1(jq)
            IF ((abs(zp) .LT. eps).AND.(abs(zq) .LT. eps)) THEN
                d(ip,jq) = 0.
            ELSE
                d(ip,jq) = (pnleg(zp,nm1)/pnleg(zq,nm1)
     $                     -im12(ip,jq))/(zp-zq)
            ENDIF
            dt(jq,ip) = d(ip,jq)
 10   CONTINUE
      RETURN
      END
C
      SUBROUTINE dgljgj (D,DT,ZGL,ZG,IGLG,NPGL,NPG,ND1,ND2,ALPHA,BETA)
C-----------------------------------------------------------------------
C
C     Compute the (one-dimensional) derivative matrix D and its
C     transpose DT associated with taking the derivative of a variable
C     expanded on a Gauss-Lobatto Jacobi mesh (M1), and evaluate its
C     derivative on a Guass Jacobi mesh (M2).
C     Need the one-dimensional interpolation operator IM12
C     (see subroutine IGLJGJ).
C     Note: D and DT are rectangular matrices.
C     Single precision version.
C
C-----------------------------------------------------------------------
      REAL D(ND2,ND1), DT(ND1,ND2), ZGL(ND1), ZG(ND2), IGLG(ND2,ND1)
      parameter(nmax=84)
      parameter(ndd = nmax)
      real*8  dd(ndd,ndd), dtd(ndd,ndd)
      real*8  zgd(ndd), zgld(ndd), iglgd(ndd,ndd)
      real*8  alphad, betad
C
      IF (npgl.LE.1) THEN
       WRITE(6,*) 'DGLJGJ: Minimum number of Gauss-Lobatto points is 2'
       call exitt
      ENDIF
      IF (npgl.GT.nmax) THEN
         WRITE(6,*) 'Polynomial degree too high in DGLJGJ'
         WRITE(6,*) 'Maximum polynomial degree is',nmax
         WRITE(6,*) 'Here NPGL=',npgl
         call exitt
      ENDIF
      IF ((alpha.LE.-1.).OR.(beta.LE.-1.)) THEN
         WRITE(6,*) 'DGLJGJ: Alpha and Beta must be greater than -1'
         call exitt
      ENDIF
C
      alphad = alpha
      betad  = beta
      DO 100 i=1,npg
         zgd(i) = zg(i)
         DO 100 j=1,npgl
            iglgd(i,j) = iglg(i,j)
 100  CONTINUE
      DO 200 i=1,npgl
         zgld(i) = zgl(i)
 200  CONTINUE
      CALL dgljgjd (dd,dtd,zgld,zgd,iglgd,npgl,npg,ndd,ndd,alphad,betad)
      DO 300 i=1,npg
      DO 300 j=1,npgl
         d(i,j)  = dd(i,j)
         dt(j,i) = dtd(j,i)
 300  CONTINUE
      RETURN
      END
C
      SUBROUTINE dgljgjd (D,DT,ZGL,ZG,IGLG,NPGL,NPG,ND1,ND2,ALPHA,BETA)
C-----------------------------------------------------------------------
C
C     Compute the (one-dimensional) derivative matrix D and its
C     transpose DT associated with taking the derivative of a variable
C     expanded on a Gauss-Lobatto Jacobi mesh (M1), and evaluate its
C     derivative on a Guass Jacobi mesh (M2).
C     Need the one-dimensional interpolation operator IM12
C     (see subroutine IGLJGJ).
C     Note: D and DT are rectangular matrices.
C     Double precision version.
C
C-----------------------------------------------------------------------
      IMPLICIT REAL*8  (a-h,o-z)
      real*8  d(nd2,nd1), dt(nd1,nd2), zgl(nd1), zg(nd2)
      real*8  iglg(nd2,nd1), alpha, beta
C
      IF (npgl.LE.1) THEN
       WRITE(6,*) 'DGLJGJD: Minimum number of Gauss-Lobatto points is 2'
       call exitt
      ENDIF
      IF ((alpha.LE.-1.).OR.(beta.LE.-1.)) THEN
         WRITE(6,*) 'DGLJGJD: Alpha and Beta must be greater than -1'
         call exitt
      ENDIF
C
      eps    = 1.e-6
      one    = 1.
      two    = 2.
      ngl    = npgl-1
      dn     = ((ngl))
      eigval = -dn*(dn+alpha+beta+one)
C
      DO 100 i=1,npg
      DO 100 j=1,npgl
         dz = abs(zg(i)-zgl(j))
         IF (dz.LT.eps) THEN
            d(i,j) = (alpha*(one+zg(i))-beta*(one-zg(i)))/
     $               (two*(one-zg(i)**2))
         ELSE
            CALL jacobf (pi,pdi,pm1,pdm1,pm2,pdm2,ngl,alpha,beta,zg(i))
            CALL jacobf (pj,pdj,pm1,pdm1,pm2,pdm2,ngl,alpha,beta,zgl(j))
            faci   = alpha*(one+zg(i))-beta*(one-zg(i))
            facj   = alpha*(one+zgl(j))-beta*(one-zgl(j))
            const  = eigval*pj+facj*pdj
            d(i,j) = ((eigval*pi+faci*pdi)*(zg(i)-zgl(j))
     $               -(one-zg(i)**2)*pdi)/(const*(zg(i)-zgl(j))**2)
         ENDIF
         dt(j,i) = d(i,j)
 100  CONTINUE
      RETURN
      END
C
      SUBROUTINE iglm (I12,IT12,Z1,Z2,lz1,lz2,ND1,ND2)
C----------------------------------------------------------------------
C
C     Compute the one-dimensional interpolation operator (matrix) I12
C     ands its transpose IT12 for interpolating a variable from a
C     Gauss Legendre mesh (1) to a another mesh M (2).
C     Z1 : lz1 Gauss Legendre points.
C     Z2 : lz2 points on mesh M.
C
C--------------------------------------------------------------------
      REAL I12(ND2,ND1),IT12(ND1,ND2),Z1(ND1),Z2(ND2)
      IF (lz1 .EQ. 1) THEN
         i12(1,1) = 1.
         it12(1,1) = 1.
         RETURN
      ENDIF
      DO 10 i=1,lz2
         zi = z2(i)
         DO 10 j=1,lz1
            i12(i,j) = hgl(j,zi,z1,lz1)
            it12(j,i) = i12(i,j)
 10   CONTINUE
      RETURN
      END
c
      SUBROUTINE igllm (I12,IT12,Z1,Z2,lz1,lz2,ND1,ND2)
C----------------------------------------------------------------------
C
C     Compute the one-dimensional interpolation operator (matrix) I12
C     ands its transpose IT12 for interpolating a variable from a
C     Gauss-Lobatto Legendre mesh (1) to a another mesh M (2).
C     Z1 : lz1 Gauss-Lobatto Legendre points.
C     Z2 : lz2 points on mesh M.
C
C--------------------------------------------------------------------
      REAL I12(ND2,ND1),IT12(ND1,ND2),Z1(ND1),Z2(ND2)
      IF (lz1 .EQ. 1) THEN
         i12(1,1) = 1.
         it12(1,1) = 1.
         RETURN
      ENDIF
      DO 10 i=1,lz2
         zi = z2(i)
         DO 10 j=1,lz1
            i12(i,j) = hgll(j,zi,z1,lz1)
            it12(j,i) = i12(i,j)
 10   CONTINUE
      RETURN
      END
C
      SUBROUTINE igjm (I12,IT12,Z1,Z2,lz1,lz2,ND1,ND2,ALPHA,BETA)
C----------------------------------------------------------------------
C
C     Compute the one-dimensional interpolation operator (matrix) I12
C     ands its transpose IT12 for interpolating a variable from a
C     Gauss Jacobi mesh (1) to a another mesh M (2).
C     Z1 : lz1 Gauss Jacobi points.
C     Z2 : lz2 points on mesh M.
C     Single precision version.
C
C--------------------------------------------------------------------
      REAL I12(ND2,ND1),IT12(ND1,ND2),Z1(ND1),Z2(ND2)
      IF (lz1 .EQ. 1) THEN
         i12(1,1) = 1.
         it12(1,1) = 1.
         RETURN
      ENDIF
      DO 10 i=1,lz2
         zi = z2(i)
         DO 10 j=1,lz1
            i12(i,j) = hgj(j,zi,z1,lz1,alpha,beta)
            it12(j,i) = i12(i,j)
 10   CONTINUE
      RETURN
      END
c
      SUBROUTINE igljm (I12,IT12,Z1,Z2,lz1,lz2,ND1,ND2,ALPHA,BETA)
C----------------------------------------------------------------------
C
C     Compute the one-dimensional interpolation operator (matrix) I12
C     ands its transpose IT12 for interpolating a variable from a
C     Gauss-Lobatto Jacobi mesh (1) to a another mesh M (2).
C     Z1 : lz1 Gauss-Lobatto Jacobi points.
C     Z2 : lz2 points on mesh M.
C     Single precision version.
C
C--------------------------------------------------------------------
      REAL I12(ND2,ND1),IT12(ND1,ND2),Z1(ND1),Z2(ND2)
      IF (lz1 .EQ. 1) THEN
         i12(1,1) = 1.
         it12(1,1) = 1.
         RETURN
      ENDIF
      DO 10 i=1,lz2
         zi = z2(i)
         DO 10 j=1,lz1
            i12(i,j) = hglj(j,zi,z1,lz1,alpha,beta)
            it12(j,i) = i12(i,j)
 10   CONTINUE
      RETURN
      END