speclib.f90

Go to the documentation of this file.

! Copyright (c) 2008-2020, UCHICAGO ARGONNE, LLC.
!
! The UChicago Argonne, LLC as Operator of Argonne National
! Laboratory holds copyright in the Software. The copyright holder
! reserves all rights except those expressly granted to licensees,
! and U.S. Government license rights.
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
!
! 1. Redistributions of source code must retain the above copyright
! notice, this list of conditions and the disclaimer below.
!
! 2. Redistributions in binary form must reproduce the above copyright
! notice, this list of conditions and the disclaimer (as noted below)
! in the documentation and/or other materials provided with the
! distribution.
!
! 3. Neither the name of ANL nor the names of its contributors
! may be used to endorse or promote products derived from this software
! without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
! FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL
! UCHICAGO ARGONNE, LLC, THE U.S. DEPARTMENT OF
! ENERGY OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
! SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
! TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
! DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
!
! Additional BSD Notice
! ---------------------
! 1. This notice is required to be provided under our contract with
! the U.S. Department of Energy (DOE). This work was produced at
! Argonne National Laboratory under Contract
! No. DE-AC02-06CH11357 with the DOE.
!
! 2. Neither the United States Government nor UCHICAGO ARGONNE,
! LLC nor any of their employees, makes any warranty,
! express or implied, or assumes any liability or responsibility for the
! accuracy, completeness, or usefulness of any information, apparatus,
! product, or process disclosed, or represents that its use would not
! infringe privately-owned rights.
!
! 3. Also, reference herein to any specific commercial products, process,
! or services by trade name, trademark, manufacturer or otherwise does
! not necessarily constitute or imply its endorsement, recommendation,
! or favoring by the United States Government or UCHICAGO ARGONNE LLC.
! The views and opinions of authors expressed
! herein do not necessarily state or reflect those of the United States
! Government or UCHICAGO ARGONNE, LLC, and shall
! not be used for advertising or product endorsement purposes.
!
!==============================================================================
!
!
!------------------------------------------------------------------------------
!
!     ABBRIVIATIONS:
!
!     M   - Set of mesh points
!     Z   - Set of collocation/quadrature points
!     W   - Set of quadrature weights
!     H   - Lagrangian interpolant
!     D   - Derivative operator
!     I   - Interpolation operator
!     GL  - Gauss Legendre
!     GLL - Gauss-Lobatto Legendre
!     GJ  - Gauss Jacobi
!     GLJ - Gauss-Lobatto Jacobi
!
!
!     MAIN ROUTINES:
!
!     Points and weights:
!
!     ZWGL      Compute Gauss Legendre points and weights
!     ZWGLL     Compute Gauss-Lobatto Legendre points and weights
!     ZWGJ      Compute Gauss Jacobi points and weights (general)
!     ZWGLJ     Compute Gauss-Lobatto Jacobi points and weights (general)
!
!     Lagrangian interpolants:
!
!     HGL       Compute Gauss Legendre Lagrangian interpolant
!     HGLL      Compute Gauss-Lobatto Legendre Lagrangian interpolant
!     HGJ       Compute Gauss Jacobi Lagrangian interpolant (general)
!     HGLJ      Compute Gauss-Lobatto Jacobi Lagrangian interpolant (general)
!
!     Derivative operators:
!
!     DGLL      Compute Gauss-Lobatto Legendre derivative matrix
!     DGLLGL    Compute derivative matrix for a staggered mesh (GLL->GL)
!     DGJ       Compute Gauss Jacobi derivative matrix (general)
!     DGLJ      Compute Gauss-Lobatto Jacobi derivative matrix (general)
!     DGLJGJ    Compute derivative matrix for a staggered mesh (GLJ->GJ) (general)
!
!     Interpolation operators:
!
!     IGLM      Compute interpolation operator GL  -> M
!     IGLLM     Compute interpolation operator GLL -> M
!     IGJM      Compute interpolation operator GJ  -> M  (general)
!     IGLJM     Compute interpolation operator GLJ -> M  (general)
!
!     Other:
!
!     PNLEG     Compute Legendre polynomial of degree N
!     PNDLEG    Compute derivative of Legendre polynomial of degree N
!
!     Comments:
!
!     Note that many of the above routines exist in both single and
!     double precision. If the name of the single precision routine is
!     SUB, the double precision version is called SUBD. In most cases
!     all the "low-level" arithmetic is done in double precision, even
!     for the single precsion versions.
!
!     Useful references:
!
! [1] Gabor Szego: Orthogonal Polynomials, American Mathematical Society,
!     Providence, Rhode Island, 1939.
! [2] Abramowitz & Stegun: Handbook of Mathematical Functions,
!     Dover, New York, 1972.
! [3] Canuto, Hussaini, Quarteroni & Zang: Spectral Methods in Fluid
!     Dynamics, Springer-Verlag, 1988.
!
!
!==============================================================================
module speclib
  use num_types, only : rp
  use utils, only: neko_error
 
 
contains
  SUBROUTINE zwgl (Z,W,NP)
    REAL(KIND=rp) z(1),w(1), alpha, beta
    alpha = 0.
    beta  = 0.
    CALL zwgj (z,w,np,alpha,beta)
    RETURN
  end subroutine zwgl
 
  SUBROUTINE zwgll (Z,W,NP)
!--------------------------------------------------------------------
!
!     Generate NP Gauss-Lobatto Legendre points (Z) and weights (W)
!     associated with Jacobi polynomial P(N)(alpha=0,beta=0).
!     The polynomial degree N=NP-1.
!     Z and W are in single precision, but all the arithmetic
!     operations are done in double precision.
!
!--------------------------------------------------------------------
    REAL(KIND=rp) z(1),w(1), alpha, beta
    alpha = 0.
    beta  = 0.
    CALL zwglj (z,w,np,alpha,beta)
    RETURN
  end subroutine zwgll
 
  SUBROUTINE zwgj (Z,W,NP,ALPHA,BETA)
!--------------------------------------------------------------------
!
!     Generate NP GAUSS JACOBI points (Z) and weights (W)
!     associated with Jacobi polynomial P(N)(alpha>-1,beta>-1).
!     The polynomial degree N=NP-1.
!     Single precision version.
!
!--------------------------------------------------------------------
    parameter(nmax=84)
    parameter(nzd = nmax)
    REAL(KIND=rp)  zd(nzd),wd(nzd),alphad,betad
    REAL(KIND=rp) z(1),w(1),alpha,beta
 
    npmax = nzd
    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 neko_error
    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 subroutine zwgj
 
  SUBROUTINE zwgjd (Z,W,NP,ALPHA,BETA)
!--------------------------------------------------------------------
!
!     Generate NP GAUSS JACOBI points (Z) and weights (W)
!     associated with Jacobi polynomial P(N)(alpha>-1,beta>-1).
!     The polynomial degree N=NP-1.
!     Double precision version.
!
!--------------------------------------------------------------------
    IMPLICIT REAL(KIND=RP)  (a-h,o-z)
    REAL(KIND=rp)  z(1),w(1),alpha,beta
 
    n     = np-1
    dn    = ((n))
    one   = 1.
    two   = 2.
    apb   = alpha+beta
 
    IF (np.LE.0) THEN
       WRITE (6,*) 'ZWGJD: Minimum number of Gauss points is 1',np
       call neko_error
    ENDIF
    IF ((alpha.LE.-one).OR.(beta.LE.-one)) THEN
       WRITE (6,*) 'ZWGJD: Alpha and Beta must be greater than -1'
       call neko_error
    ENDIF
 
    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
 
    CALL jacg (z,np,alpha,beta)
 
    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 subroutine zwgjd
 
  SUBROUTINE zwglj (Z,W,NP,ALPHA,BETA)
!--------------------------------------------------------------------
!
!     Generate NP GAUSS LOBATTO JACOBI points (Z) and weights (W)
!     associated with Jacobi polynomial P(N)(alpha>-1,beta>-1).
!     The polynomial degree N=NP-1.
!     Single precision version.
!
!--------------------------------------------------------------------
    parameter(nmax=84)
    parameter(nzd = nmax)
    REAL(KIND=rp)  zd(nzd),wd(nzd),alphad,betad
    REAL(KIND=rp) z(1),w(1),alpha,beta
 
    npmax = nzd
    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 neko_error
    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 subroutine zwglj
 
  SUBROUTINE zwgljd (Z,W,NP,ALPHA,BETA)
!--------------------------------------------------------------------
!
!     Generate NP GAUSS LOBATTO JACOBI points (Z) and weights (W)
!     associated with Jacobi polynomial P(N)(alpha>-1,beta>-1).
!     The polynomial degree N=NP-1.
!     Double precision version.
!
!--------------------------------------------------------------------
    IMPLICIT REAL(KIND=RP)  (a-h,o-z)
    REAL(KIND=rp)  z(np),w(np),alpha,beta
 
    n     = np-1
    nm1   = n-1
    one   = 1.
    two   = 2.
 
    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 neko_error
    ENDIF
    IF ((alpha.LE.-one).OR.(beta.LE.-one)) THEN
       WRITE (6,*) 'ZWGLJD: Alpha and Beta must be greater than -1'
       call neko_error
    ENDIF
 
    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)
 
!      RETURN
  end subroutine zwgljd
 
  REAL(kind=rp)  FUNCTION endw1 (N,ALPHA,BETA)
    IMPLICIT REAL(KIND=RP)  (a-h,o-z)
    REAL(kind=rp)  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 function endw1
 
  REAL(kind=rp)  FUNCTION endw2 (N,ALPHA,BETA)
    IMPLICIT REAL(KIND=RP)  (a-h,o-z)
    REAL(kind=rp)  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 function endw2
 
  REAL(kind=rp)  FUNCTION gammaf (X)
    IMPLICIT REAL(KIND=RP)  (a-h,o-z)
    REAL(kind=rp)  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 function gammaf
 
  REAL(kind=rp)  FUNCTION pnormj (N,ALPHA,BETA)
    IMPLICIT REAL(KIND=RP)  (a-h,o-z)
    REAL(kind=rp)  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 function pnormj
 
  SUBROUTINE jacg (XJAC,NP,ALPHA,BETA)
!--------------------------------------------------------------------
!
!     Compute NP Gauss points XJAC, which are the zeros of the
!     Jacobi polynomial J(NP) with parameters ALPHA and BETA.
!     ALPHA and BETA determines the specific type of Gauss points.
!     Examples:
!     ALPHA = BETA =  0.0  ->  Legendre points
!     ALPHA = BETA = -0.5  ->  Chebyshev points
!
!--------------------------------------------------------------------
    IMPLICIT REAL(KIND=RP)  (a-h,o-z)
    REAL(KIND=rp)  xjac(1)
    DATA kstop /10/
    DATA eps/1.0e-12_rp/
    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 subroutine jacg
 
  SUBROUTINE jacobf (POLY,PDER,POLYM1,PDERM1,POLYM2,PDERM2,N,ALP,BET,X)
!--------------------------------------------------------------------
!
!     Computes the Jacobi polynomial (POLY) and its derivative (PDER)
!     of degree N at X.
!
!--------------------------------------------------------------------
    IMPLICIT REAL(KIND=RP)  (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 subroutine jacobf
 
  REAL(kind=rp) FUNCTION hgj (II,Z,ZGJ,NP,ALPHA,BETA)
!---------------------------------------------------------------------
!
!     Compute the value of the Lagrangian interpolant HGJ through
!     the NP Gauss Jacobi points ZGJ at the point Z.
!     Single precision version.
!
!---------------------------------------------------------------------
    parameter(nmax=84)
    parameter(nzd = nmax)
    REAL(kind=rp)  zd,zgjd(nzd),alphad,betad
    REAL(kind=rp)  z,zgj(1),alpha,beta
    npmax = nzd
    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 neko_error
    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 function hgj
 
  REAL(kind=rp)  FUNCTION hgjd (II,Z,ZGJ,NP,ALPHA,BETA)
!---------------------------------------------------------------------
!
!     Compute the value of the Lagrangian interpolant HGJD through
!     the NZ Gauss-Lobatto Jacobi points ZGJ at the point Z.
!     Double precision version.
!
!---------------------------------------------------------------------
    IMPLICIT REAL(KIND=RP)  (a-h,o-z)
    REAL(kind=rp)  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 function hgjd
 
  REAL(kind=rp) FUNCTION hglj (II,Z,ZGLJ,NP,ALPHA,BETA)
!---------------------------------------------------------------------
!
!     Compute the value of the Lagrangian interpolant HGLJ through
!     the NZ Gauss-Lobatto Jacobi points ZGLJ at the point Z.
!     Single precision version.
!
!---------------------------------------------------------------------
    parameter(nmax=84)
    parameter(nzd = nmax)
    REAL(kind=rp)  zd,zgljd(nzd),alphad,betad
    REAL(kind=rp)  z,zglj(1),alpha,beta
    npmax = nzd
    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 neko_error
    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 function hglj
 
  REAL(kind=rp)  FUNCTION hgljd (I,Z,ZGLJ,NP,ALPHA,BETA)
!---------------------------------------------------------------------
!
!     Compute the value of the Lagrangian interpolant HGLJD through
!     the NZ Gauss-Lobatto Jacobi points ZJACL at the point Z.
!     Double precision version.
!
!---------------------------------------------------------------------
    IMPLICIT REAL(KIND=RP)  (a-h,o-z)
    REAL(kind=rp)  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 function hgljd
 
  SUBROUTINE dgj (D,DT,Z,NZ,NZD,ALPHA,BETA)
!-----------------------------------------------------------------
!
!     Compute the derivative matrix D and its transpose DT
!     associated with the Nth order Lagrangian interpolants
!     through the NZ Gauss Jacobi points Z.
!     Note: D and DT are square matrices.
!     Single precision version.
!
!-----------------------------------------------------------------
    parameter(nmax=84)
    parameter(nzdd = nmax)
    REAL(KIND=rp)  dd(nzdd,nzdd),dtd(nzdd,nzdd),zd(nzdd),alphad,betad
    REAL(KIND=rp) d(nzd,nzd),dt(nzd,nzd),z(1),alpha,beta
 
    IF (nz.LE.0) THEN
       WRITE (6,*) 'DGJ: Minimum number of Gauss points is 1'
       call neko_error
    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 neko_error
    ENDIF
    IF ((alpha.LE.-1.).OR.(beta.LE.-1.)) THEN
       WRITE (6,*) 'DGJ: Alpha and Beta must be greater than -1'
       call neko_error
    ENDIF
    alphad = alpha
    betad  = beta
    DO 100 i=1,nz
       zd(i) = z(i)
100 CONTINUE
    CALL dgjd (dd,dtd,zd,nz,nzdd,alphad,betad)
    DO i=1,nz
       DO j=1,nz
          d(i,j)  = dd(i,j)
          dt(i,j) = dtd(i,j)
       END DO
    END DO
    RETURN
  end subroutine dgj
 
  SUBROUTINE dgjd (D,DT,Z,NZ,NZD,ALPHA,BETA)
!-----------------------------------------------------------------
!
!     Compute the derivative matrix D and its transpose DT
!     associated with the Nth order Lagrangian interpolants
!     through the NZ Gauss Jacobi points Z.
!     Note: D and DT are square matrices.
!     Double precision version.
!
!-----------------------------------------------------------------
    IMPLICIT REAL(KIND=RP)  (a-h,o-z)
    REAL(KIND=rp)  d(nzd,nzd),dt(nzd,nzd),z(1),alpha,beta
    n    = nz-1
    dn   = ((n))
    one  = 1.
    two  = 2.
 
    IF (nz.LE.1) THEN
       WRITE (6,*) 'DGJD: Minimum number of Gauss-Lobatto points is 2'
       call neko_error
    ENDIF
    IF ((alpha.LE.-one).OR.(beta.LE.-one)) THEN
       WRITE (6,*) 'DGJD: Alpha and Beta must be greater than -1'
       call neko_error
    ENDIF
 
    DO i=1,nz
       DO 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)
       END DO
    END DO
    RETURN
  end subroutine dgjd
 
  SUBROUTINE dglj (D,DT,Z,NZ,NZD,ALPHA,BETA)
!-----------------------------------------------------------------
!
!     Compute the derivative matrix D and its transpose DT
!     associated with the Nth order Lagrangian interpolants
!     through the NZ Gauss-Lobatto Jacobi points Z.
!     Note: D and DT are square matrices.
!     Single precision version.
!
!-----------------------------------------------------------------
    parameter(nmax=84)
    parameter(nzdd = nmax)
    REAL(KIND=rp)  dd(nzdd,nzdd),dtd(nzdd,nzdd),zd(nzdd),alphad,betad
    REAL(KIND=rp) d(nzd,nzd),dt(nzd,nzd),z(1),alpha,beta
 
    IF (nz.LE.1) THEN
       WRITE (6,*) 'DGLJ: Minimum number of Gauss-Lobatto points is 2'
       call neko_error
    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 neko_error
    ENDIF
    IF ((alpha.LE.-1.).OR.(beta.LE.-1.)) THEN
       WRITE (6,*) 'DGLJ: Alpha and Beta must be greater than -1'
       call neko_error
    ENDIF
    alphad = alpha
    betad  = beta
    DO 100 i=1,nz
       zd(i) = z(i)
100 CONTINUE
    CALL dgljd (dd,dtd,zd,nz,nzdd,alphad,betad)
    DO i=1,nz
       DO j=1,nz
          d(i,j)  = dd(i,j)
          dt(i,j) = dtd(i,j)
       END DO
    END DO
    RETURN
  end subroutine dglj
 
  SUBROUTINE dgljd (D,DT,Z,NZ,NZD,ALPHA,BETA)
!-----------------------------------------------------------------
!
!     Compute the derivative matrix D and its transpose DT
!     associated with the Nth order Lagrangian interpolants
!     through the NZ Gauss-Lobatto Jacobi points Z.
!     Note: D and DT are square matrices.
!     Double precision version.
!
!-----------------------------------------------------------------
    IMPLICIT REAL(KIND=RP)  (a-h,o-z)
    REAL(KIND=rp)  d(nzd,nzd),dt(nzd,nzd),z(1),alpha,beta
    n    = nz-1
    dn   = ((n))
    one  = 1.
    two  = 2.
    eigval = -dn*(dn+alpha+beta+one)
 
    IF (nz.LE.1) THEN
       WRITE (6,*) 'DGLJD: Minimum number of Gauss-Lobatto points is 2'
       call neko_error
    ENDIF
    IF ((alpha.LE.-one).OR.(beta.LE.-one)) THEN
       WRITE (6,*) 'DGLJD: Alpha and Beta must be greater than -1'
       call neko_error
    ENDIF
 
    DO i=1,nz
       DO 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)
       END DO
    END DO
    RETURN
  end subroutine dgljd
 
  SUBROUTINE dgll (D,DT,Z,NZ,NZD)
!-----------------------------------------------------------------
!
!     Compute the derivative matrix D and its transpose DT
!     associated with the Nth order Lagrangian interpolants
!     through the NZ Gauss-Lobatto Legendre points Z.
!     Note: D and DT are square matrices.
!
!-----------------------------------------------------------------
    IMPLICIT REAL(KIND=RP)  (a-h,o-z)
    parameter(nmax=84)
    REAL(KIND=rp) d(nzd,nzd),dt(nzd,nzd),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 i=1,nz
       DO 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)
       END DO
    END DO
    RETURN
  end subroutine dgll
 
  REAL(kind=rp) FUNCTION hgll (I,Z,ZGLL,NZ)
!---------------------------------------------------------------------
!
!     Compute the value of the Lagrangian interpolant L through
!     the NZ Gauss-Lobatto Legendre points ZGLL at the point Z.
!
!---------------------------------------------------------------------
    IMPLICIT REAL(KIND=RP)  (a-h,o-z)
    REAL(kind=rp) zgll(1), eps, dz, z
    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 function hgll
 
  REAL(kind=rp) FUNCTION hgl (I,Z,ZGL,NZ)
!---------------------------------------------------------------------
!
!     Compute the value of the Lagrangian interpolant HGL through
!     the NZ Gauss Legendre points ZGL at the point Z.
!
!---------------------------------------------------------------------
    REAL(kind=rp) zgl(1), z, eps, dz
    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 function hgl
 
  REAL(kind=rp) FUNCTION pnleg (Z,N)
!---------------------------------------------------------------------
!
!     Compute the value of the Nth order Legendre polynomial at Z.
!     (Simpler than JACOBF)
!     Based on the recursion formula for the Legendre polynomials.
!
!---------------------------------------------------------------------
!
!     This next statement is to overcome the underflow bug in the i860.
!     It can be removed at a later date.  11 Aug 1990   pff.
!
    IMPLICIT REAL(KIND=RP)  (a-h,o-z)
    REAL(kind=rp) z, p1, p2, p3
    IF(abs(z) .LT. 1.0e-25) z = 0.0
 
 
    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 function pnleg
 
  REAL(kind=rp) FUNCTION pndleg (Z,N)
!----------------------------------------------------------------------
!
!     Compute the derivative of the Nth order Legendre polynomial at Z.
!     (Simpler than JACOBF)
!     Based on the recursion formula for the Legendre polynomials.
!
!----------------------------------------------------------------------
    IMPLICIT REAL(KIND=RP)  (a-h,o-z)
    REAL(kind=rp) p1, p2, p1d, p2d, p3d, z
    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 function pndleg
 
  SUBROUTINE dgllgl (D,DT,ZM1,ZM2,IM12,NZM1,NZM2,ND1,ND2)
!-----------------------------------------------------------------------
!
!     Compute the (one-dimensional) derivative matrix D and its
!     transpose DT associated with taking the derivative of a variable
!     expanded on a Gauss-Lobatto Legendre mesh (M1), and evaluate its
!     derivative on a Guass Legendre mesh (M2).
!     Need the one-dimensional interpolation operator IM12
!     (see subroutine IGLLGL).
!     Note: D and DT are rectangular matrices.
!
!-----------------------------------------------------------------------
    REAL(KIND=rp) d(nd2,nd1), dt(nd1,nd2), zm1(nd1), zm2(nd2), im12(nd2,nd1)
    REAL(KIND=rp) eps, zp, zq
    IF (nzm1.EQ.1) THEN
       d(1,1) = 0.
       dt(1,1) = 0.
       RETURN
    ENDIF
    eps = 1.e-6
    nm1 = nzm1-1
    DO ip = 1, nzm2
       DO 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)
       END DO
    END DO
    RETURN
  end subroutine dgllgl
 
  SUBROUTINE dgljgj (D,DT,ZGL,ZG,IGLG,NPGL,NPG,ND1,ND2,ALPHA,BETA)
!-----------------------------------------------------------------------
!
!     Compute the (one-dimensional) derivative matrix D and its
!     transpose DT associated with taking the derivative of a variable
!     expanded on a Gauss-Lobatto Jacobi mesh (M1), and evaluate its
!     derivative on a Guass Jacobi mesh (M2).
!     Need the one-dimensional interpolation operator IM12
!     (see subroutine IGLJGJ).
!     Note: D and DT are rectangular matrices.
!     Single precision version.
!
!-----------------------------------------------------------------------
    REAL(KIND=rp) d(nd2,nd1), dt(nd1,nd2), zgl(nd1), zg(nd2), iglg(nd2,nd1)
    parameter(nmax=84)
    parameter(ndd = nmax)
    REAL(KIND=rp)  dd(ndd,ndd), dtd(ndd,ndd)
    REAL(KIND=rp)  zgd(ndd), zgld(ndd), iglgd(ndd,ndd)
    REAL(KIND=rp)  alphad, betad
 
    IF (npgl.LE.1) THEN
       WRITE(6,*) 'DGLJGJ: Minimum number of Gauss-Lobatto points is 2'
       call neko_error
    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 neko_error
    ENDIF
    IF ((alpha.LE.-1.).OR.(beta.LE.-1.)) THEN
       WRITE(6,*) 'DGLJGJ: Alpha and Beta must be greater than -1'
       call neko_error
    ENDIF
 
    alphad = alpha
    betad  = beta
    DO i=1,npg
       zgd(i) = zg(i)
       DO j=1,npgl
          iglgd(i,j) = iglg(i,j)
       END DO
    END DO
    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 i=1,npg
       DO j=1,npgl
          d(i,j)  = dd(i,j)
          dt(j,i) = dtd(j,i)
       END DO
    END DO
    RETURN
  end subroutine dgljgj
 
  SUBROUTINE dgljgjd (D,DT,ZGL,ZG,IGLG,NPGL,NPG,ND1,ND2,ALPHA,BETA)
!-----------------------------------------------------------------------
!
!     Compute the (one-dimensional) derivative matrix D and its
!     transpose DT associated with taking the derivative of a variable
!     expanded on a Gauss-Lobatto Jacobi mesh (M1), and evaluate its
!     derivative on a Guass Jacobi mesh (M2).
!     Need the one-dimensional interpolation operator IM12
!     (see subroutine IGLJGJ).
!     Note: D and DT are rectangular matrices.
!     Double precision version.
!
!-----------------------------------------------------------------------
    IMPLICIT REAL(KIND=RP)  (a-h,o-z)
    REAL(KIND=rp)  d(nd2,nd1), dt(nd1,nd2), zgl(nd1), zg(nd2)
    REAL(KIND=rp)  iglg(nd2,nd1), alpha, beta
 
    IF (npgl.LE.1) THEN
       WRITE(6,*) 'DGLJGJD: Minimum number of Gauss-Lobatto points is 2'
       call neko_error
    ENDIF
    IF ((alpha.LE.-1.).OR.(beta.LE.-1.)) THEN
       WRITE(6,*) 'DGLJGJD: Alpha and Beta must be greater than -1'
       call neko_error
    ENDIF
 
    eps    = 1.e-6
    one    = 1.
    two    = 2.
    ngl    = npgl-1
    dn     = ((ngl))
    eigval = -dn*(dn+alpha+beta+one)
 
    DO i=1,npg
       DO 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)
       END DO
    END DO
    RETURN
  end subroutine dgljgjd
 
  SUBROUTINE iglm (I12,IT12,Z1,Z2,NZ1,NZ2,ND1,ND2)
!----------------------------------------------------------------------
!
!     Compute the one-dimensional interpolation operator (matrix) I12
!     ands its transpose IT12 for interpolating a variable from a
!     Gauss Legendre mesh (1) to a another mesh M (2).
!     Z1 : NZ1 Gauss Legendre points.
!     Z2 : NZ2 points on mesh M.
!
!--------------------------------------------------------------------
    REAL(KIND=rp) i12(nd2,nd1),it12(nd1,nd2),z1(nd1),z2(nd2), zi
    IF (nz1 .EQ. 1) THEN
       i12(1,1) = 1.
       it12(1,1) = 1.
       RETURN
    ENDIF
    DO i=1,nz2
       zi = z2(i)
       DO j=1,nz1
          i12(i,j) = hgl(j,zi,z1,nz1)
          it12(j,i) = i12(i,j)
       END DO
    END DO
    RETURN
  end subroutine iglm
 
  SUBROUTINE igllm (I12,IT12,Z1,Z2,NZ1,NZ2,ND1,ND2)
!----------------------------------------------------------------------
!
!     Compute the one-dimensional interpolation operator (matrix) I12
!     ands its transpose IT12 for interpolating a variable from a
!     Gauss-Lobatto Legendre mesh (1) to a another mesh M (2).
!     Z1 : NZ1 Gauss-Lobatto Legendre points.
!     Z2 : NZ2 points on mesh M.
!
!--------------------------------------------------------------------
    REAL(KIND=rp) i12(nd2,nd1),it12(nd1,nd2),z1(nd1),z2(nd2),zi
    IF (nz1 .EQ. 1) THEN
       i12(1,1) = 1.
       it12(1,1) = 1.
       RETURN
    ENDIF
    DO i=1,nz2
       zi = z2(i)
       DO j=1,nz1
          i12(i,j) = hgll(j,zi,z1,nz1)
          it12(j,i) = i12(i,j)
       END DO
    END DO
    RETURN
  end subroutine igllm
 
  SUBROUTINE igjm (I12,IT12,Z1,Z2,NZ1,NZ2,ND1,ND2,ALPHA,BETA)
!----------------------------------------------------------------------
!
!     Compute the one-dimensional interpolation operator (matrix) I12
!     ands its transpose IT12 for interpolating a variable from a
!     Gauss Jacobi mesh (1) to a another mesh M (2).
!     Z1 : NZ1 Gauss Jacobi points.
!     Z2 : NZ2 points on mesh M.
!     Single precision version.
!
!--------------------------------------------------------------------
    REAL(KIND=rp) i12(nd2,nd1),it12(nd1,nd2),z1(nd1),z2(nd2),zi,alpha,beta
    IF (nz1 .EQ. 1) THEN
       i12(1,1) = 1.
       it12(1,1) = 1.
       RETURN
    ENDIF
    DO i=1,nz2
       zi = z2(i)
       DO j=1,nz1
          i12(i,j) = hgj(j,zi,z1,nz1,alpha,beta)
          it12(j,i) = i12(i,j)
       END DO
    END DO
    RETURN
  end subroutine igjm
 
  SUBROUTINE igljm (I12,IT12,Z1,Z2,NZ1,NZ2,ND1,ND2,ALPHA,BETA)
!----------------------------------------------------------------------
!
!     Compute the one-dimensional interpolation operator (matrix) I12
!     ands its transpose IT12 for interpolating a variable from a
!     Gauss-Lobatto Jacobi mesh (1) to a another mesh M (2).
!     Z1 : NZ1 Gauss-Lobatto Jacobi points.
!     Z2 : NZ2 points on mesh M.
!     Single precision version.
!
!--------------------------------------------------------------------
    REAL(KIND=rp) i12(nd2,nd1),it12(nd1,nd2),z1(nd1),z2(nd2),zi,alpha,beta
    IF (nz1 .EQ. 1) THEN
       i12(1,1) = 1.
       it12(1,1) = 1.
       RETURN
    ENDIF
    DO i=1,nz2
       zi = z2(i)
       DO j=1,nz1
          i12(i,j) = hglj(j,zi,z1,nz1,alpha,beta)
          it12(j,i) = i12(i,j)
       END DO
    END DO
    RETURN
  end subroutine igljm
end module speclib