speclib.F90

Go to the documentation of this file.

!==============================================================================

module speclib
  use kinds, only : dp

  integer, PARAMETER :: nmax = 84
  public zwgl, zwgll, igllm, dgll, dgllgl, zwglj
  private

contains

!--------------------------------------------------------------------
    SUBROUTINE zwgl (Z,W,NP)
!--------------------------------------------------------------------
!--------------------------------------------------------------------
    implicit none
    integer, intent(in) :: NP
    REAL(DP), intent(out) :: Z(np), W(np)
    real(DP) :: alpha, beta
    alpha = 0._dp
    beta  = 0._dp
    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.
!--------------------------------------------------------------------
    implicit none
    integer, intent(in) :: NP
    REAL(DP), intent(out) :: Z(np),W(np)
    real(DP) :: alpha, beta
    alpha = 0._dp
    beta  = 0._dp
    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.

!--------------------------------------------------------------------
    implicit none
    integer, parameter :: NZD = NMAX

    integer,  intent(in)  :: NP
    REAL(DP), intent(out) :: Z(np), W(np)
    real(DP), intent(in)  :: ALPHA, BETA
    REAL(DP) :: ZD(nzd),WD(nzd),ALPHAD,BETAD
    integer :: npmax, i

    npmax = nzd
    IF (np > 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 END DO
    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 none
    integer, intent(in)   :: NP
    REAL(DP), intent(out) :: Z(np),W(np)
    real(DP), intent(in)  :: ALPHA,BETA

    integer :: n, np1, np2, i
    real(DP) :: dn, one, two, apb, dnp1, dnp2, fac1, fac2, fac3, fnorm
    real(DP) :: rcoef, P, PD, PM1, PDM1, PM2, PDM2

    n     = np-1
    dn    = ((n))
    one   = 1._dp
    two   = 2._dp
    apb   = alpha+beta

    IF (np <= 0) THEN
        WRITE (6,*) 'ZWGJD: Minimum number of Gauss points is 1',np
        call exitt
    ENDIF
    IF ((alpha <= -one) .OR. (beta <= -one)) THEN
        WRITE (6,*) 'ZWGJD: Alpha and Beta must be greater than -1'
        call exitt
    ENDIF

    IF (np == 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 END DO
    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.

!--------------------------------------------------------------------
    implicit none
    integer, PARAMETER :: NZD = NMAX
    integer,  intent(in)  :: NP
    REAL(DP), intent(out) :: Z(np),W(np)
    real(DP), intent(in)  :: ALPHA,BETA
    REAL(DP) ::  ZD(nzd),WD(nzd),ALPHAD,BETAD

    integer :: npmax, i

    npmax = nzd
    IF (np > 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 END DO
    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 none
    integer, intent(in) :: np
    REAL(DP), intent(out) ::  Z(np),W(np)
    REAL(DP), intent(in) ::  ALPHA,BETA

    integer :: n, nm1, i
    real(DP) :: one, two, alpg, betg, p, pd, pm1, pdm1, pm2, pdm2
    n     = np-1
    nm1   = n-1
    one   = 1._dp
    two   = 2._dp

    IF (np <= 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 <= -one) .OR. (beta <= -one)) THEN
        WRITE (6,*) 'ZWGLJD: Alpha and Beta must be greater than -1'
        call exitt
    ENDIF

    IF (nm1 > 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 END DO
    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(DP) FUNCTION endw1 (N,ALPHA,BETA)
    implicit none
    integer, intent(in) :: n
    real(DP), intent(in) :: alpha, beta

    integer :: i
    real(DP) :: zero, one, two, three, four, apb
    real(DP) :: f1, f2, fint1, fint2, f3, a1, a2, a3, abn, abnn, di
    zero  = 0._dp
    one   = 1._dp
    two   = 2._dp
    three = 3._dp
    four  = 4._dp
    f3    = -1._dp
    apb   = alpha+beta
    IF (n == 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 == 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 == 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 END DO
    endw1  = f3
    RETURN
    END FUNCTION

    REAL(DP) FUNCTION endw2 (N,ALPHA,BETA)
    implicit none
    integer, intent(in) :: n
    real(DP), intent(in) :: alpha, beta

    integer :: i
    real(DP) :: zero, one, two, three, four, apb
    real(DP) :: f1, f2, fint1, fint2, f3, a1, a2, a3, abn, abnn, di

    zero  = 0._dp
    one   = 1._dp
    two   = 2._dp
    three = 3._dp
    four  = 4._dp
    apb   = alpha+beta
    f3 = -1._dp
    IF (n == 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 == 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 == 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 END DO
    endw2  = f3
    RETURN
    END FUNCTION

    REAL(DP) FUNCTION gammaf (X)
    implicit none
    REAL(DP), intent(in) :: X

    real(DP) :: zero, half, one, two, four, pi
    zero = 0.0_dp
    half = 0.5_dp
    one  = 1.0_dp
    two  = 2.0_dp
    four = 4.0_dp
    pi   = four*atan(one)
    gammaf = one
    IF (x == -half) gammaf = -two*sqrt(pi)
    IF (x == half) gammaf =  sqrt(pi)
    IF (x == one ) gammaf =  one
    IF (x == two ) gammaf =  one
    IF (x == 1.5_dp  ) gammaf =  sqrt(pi)/2._dp
    IF (x == 2.5_dp) gammaf =  1.5_dp*sqrt(pi)/2._dp
    IF (x == 3.5_dp) gammaf =  0.5_dp*(2.5_dp*(1.5_dp*sqrt(pi)))
    IF (x == 3._dp ) gammaf =  2._dp
    IF (x == 4._dp ) gammaf = 6._dp
    IF (x == 5._dp ) gammaf = 24._dp
    IF (x == 6._dp ) gammaf = 120._dp
    RETURN
    END FUNCTION

    REAL(DP)  FUNCTION pnormj (N,ALPHA,BETA)
    implicit none
    integer, intent(in) :: n
    REAL(DP), intent(in) ::  ALPHA,BETA

    real(DP) :: one, two, dn, const, prod, dindx, frac
    integer :: i

    one   = 1._dp
    two   = 2._dp
    dn    = ((n))
    const = alpha+beta+one
    IF (n <= 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 END DO
    pnormj = prod * two**const/(two*dn+const)
    RETURN
    END FUNCTION


    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 none
    integer, intent(in) :: np
    REAL(DP), intent(out) ::  XJAC(np)
    real(DP), intent(in) :: alpha, beta

    real(DP) :: eps, one, dth, x, xlast, x1, x2, swap, xmin, delx
    real(DP) :: p, pd, pm1, pdm1, pm2, pdm2, recsum
    integer :: n, kstop, j, k, i, jmin, jm
    DATA kstop /10/

    kstop = 10
    eps = 1.e-12_dp

    n   = np-1
    one = 1._dp
    dth = 4._dp*atan(one)/(2._dp*((n))+2._dp)
    jmin = -1
    DO 40 j=1,np
        IF (j == 1) THEN
            x = cos((2._dp*(((j))-1._dp)+1._dp)*dth)
        ELSE
            x1 = cos((2._dp*(((j))-1._dp)+1._dp)*dth)
            x2 = xlast
            x  = (x1+x2)/2._dp
        ENDIF
        DO 30 k=1,kstop
            CALL jacobf(p,pd,pm1,pdm1,pm2,pdm2,np,alpha,beta,x)
            recsum = 0._dp
            jm = j-1
            DO 29 i=1,jm
                recsum = recsum+1./(x-xjac(np-i+1))
            29 END DO
            delx = -p/(pd-recsum*p)
            x    = x+delx
            IF (abs(delx) < eps) goto 31
        30 END DO
        31 CONTINUE
        xjac(np-j+1) = x
        xlast        = x
    40 END DO
    DO 200 i=1,np
        xmin = 2._dp
        DO 100 j=i,np
            IF (xjac(j) < xmin) THEN
                xmin = xjac(j)
                jmin = j
            ENDIF
        100 END DO
        IF (jmin /= i) THEN
            swap = xjac(i)
            xjac(i) = xjac(jmin)
            xjac(jmin) = swap
        ENDIF
    200 END DO
    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 none
    real(DP), intent(out) :: poly, pder, polym1, pderm1, polym2, pderm2
    real(DP), intent(in) :: alp, bet, x
    integer, intent(in) :: n

    real(DP) :: apb, polyl, pderl, dk, a1, a2, a3, b3, a4, polyn, pdern
    real(DP) :: psave, pdsave
    integer :: k

    apb  = alp+bet
    poly = 1._dp
    pder = 0._dp
    IF (n == 0) RETURN
    polyl = poly
    pderl = pder
    poly  = (alp-bet+(apb+2._dp)*x)/2._dp
    pder  = (apb+2._dp)/2._dp
    psave = 0._dp
    pdsave = 0._dp
    IF (n == 1) RETURN
    DO 20 k=2,n
        dk = ((k))
        a1 = 2._dp*dk*(dk+apb)*(2._dp*dk+apb-2._dp)
        a2 = (2._dp*dk+apb-1._dp)*(alp**2-bet**2)
        b3 = (2._dp*dk+apb-2._dp)
        a3 = b3*(b3+1._dp)*(b3+2._dp)
        a4 = 2._dp*(dk+alp-1._dp)*(dk+bet-1._dp)*(2._dp*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 END DO
    polym1 = polyl
    pderm1 = pderl
    polym2 = psave
    pderm2 = pdsave
    RETURN
    END SUBROUTINE jacobf

    REAL(DP) 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.
!---------------------------------------------------------------------
    implicit none
    integer, PARAMETER :: NZD = NMAX
    integer, intent(in) :: np, ii
    REAL(DP), intent(in) :: Z,ZGJ(np),ALPHA,BETA

    REAL(DP) ::  ZD,ZGJD(nzd),ALPHAD,BETAD
    integer :: i, npmax

    npmax = nzd
    IF (np > 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 END DO
    alphad = alpha
    betad  = beta
    hgj    = hgjd(ii,zd,zgjd,np,alphad,betad)
    RETURN
    END FUNCTION hgj

    REAL(DP) 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 none
    integer, intent(in) :: ii, np
    REAL(DP), intent(in) :: Z,ZGJ(np),ALPHA,BETA

    real(DP) :: eps, one, zi, dz
    real(DP) :: pzi, pdzi, pm1, pdm1, pm2, pdm2, pz, pdz

    eps = 1.e-5_dp
    one = 1._dp
    zi  = zgj(ii)
    dz  = z-zi
    IF (abs(dz) < 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

    REAL(DP) 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.
!---------------------------------------------------------------------
    implicit none
    integer, PARAMETER :: NZD = NMAX
    integer,  intent(in) :: ii, NP
    REAL(DP), intent(in) :: Z,ZGLJ(np)
    REAL(DP), intent(in) :: ALPHA,BETA

    REAL(DP) :: ZD,ZGLJD(nzd),ALPHAD,BETAD
    integer :: npmax, i

    npmax = nzd
    IF (np > 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 END DO
    alphad = alpha
    betad  = beta
    hglj   = hgljd(ii,zd,zgljd,np,alphad,betad)
    RETURN
    END FUNCTION hglj

    REAL(DP) 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 none
    integer,  intent(in) :: i, np
    REAL(DP), intent(in) ::  Z,ZGLJ(np),ALPHA,BETA

    real(DP) :: eps, one, dz, zi, dn, eigval, const
    real(DP) :: p, pd, pi, pdi, pm1, pdm1, pm2, pdm2
    integer :: n

    eps = 1.e-5_dp
    one = 1._dp
    zi  = zglj(i)
    dz  = z-zi
    IF (abs(dz) < 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

    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.
!-----------------------------------------------------------------
    implicit none
    integer, PARAMETER :: NZDD = NMAX
    integer :: nz, nzd
    REAL(DP), intent(out) :: D(nzd,nzd),DT(nzd,nzd)
    real(DP), intent(in) :: Z(nz),ALPHA,BETA

    REAL(DP) ::  DD(nzdd,nzdd),DTD(nzdd,nzdd),ZD(nzdd),ALPHAD,BETAD
    integer :: i, j

    IF (nz <= 0) THEN
        WRITE (6,*) 'DGJ: Minimum number of Gauss points is 1'
        call exitt
    ENDIF
    IF (nz > 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 <= -1.) .OR. (beta <= -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 END DO
    CALL dgjd(dd,dtd,zd,nz,nzdd,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 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 none
    integer,  intent(in)  :: nz, nzd
    REAL(DP), intent(out) :: D(nzd,nzd),DT(nzd,nzd)
    real(DP), intent(in)  :: Z(nz),ALPHA,BETA

    integer :: n, i, j
    real(DP) :: dn, one, two
    real(DP) :: pi, pdi, pm1, pdm1, pm2, pdm2, pj, pdj

    n    = nz-1
    dn   = ((n))
    one  = 1._dp
    two  = 2._dp

    IF (nz <= 1) THEN
        WRITE (6,*) 'DGJD: Minimum number of Gauss-Lobatto points is 2'
        call exitt
    ENDIF
    IF ((alpha <= -one) .OR. (beta <= -one)) THEN
        WRITE (6,*) 'DGJD: Alpha and Beta must be greater than -1'
        call exitt
    ENDIF

    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 /= j) d(i,j) = pdi/(pdj*(z(i)-z(j)))
            IF (i == j) d(i,j) = ((alpha+beta+two)*z(i)+alpha-beta)/ &
            (two*(one-z(i)**2))
            dt(j,i) = d(i,j)
    200 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.
!-----------------------------------------------------------------
    implicit none
    integer, PARAMETER :: NZDD = NMAX
    integer, intent(in) :: nz, nzd
    REAL(DP), intent(out) :: D(nzd,nzd),DT(nzd,nzd)
    REAL(DP), intent(in) :: Z(nz),ALPHA,BETA

    REAL(DP) :: DD(nzdd,nzdd),DTD(nzdd,nzdd),ZD(nzdd),ALPHAD,BETAD
    integer :: i, j

    IF (nz <= 1) THEN
        WRITE (6,*) 'DGLJ: Minimum number of Gauss-Lobatto points is 2'
        call exitt
    ENDIF
    IF (nz > 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 <= -1.) .OR. (beta <= -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 END DO
    CALL dgljd(dd,dtd,zd,nz,nzdd,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 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 none
    integer, intent(in) :: nz, nzd
    REAL(DP), intent(out) ::  D(nzd,nzd),DT(nzd,nzd)
    REAL(DP), intent(in) ::  Z(nz),ALPHA,BETA

    integer :: n, i, j
    real(DP) :: dn, one, two, eigval
    real(DP) :: pi, pdi, pm1, pdm1, pm2, pdm2, pj, pdj, ci, cj
    n    = nz-1
    dn   = ((n))
    one  = 1._dp
    two  = 2._dp
    eigval = -dn*(dn+alpha+beta+one)

    IF (nz <= 1) THEN
        WRITE (6,*) 'DGLJD: Minimum number of Gauss-Lobatto points is 2'
        call exitt
    ENDIF
    IF ((alpha <= -one) .OR. (beta <= -one)) THEN
        WRITE (6,*) 'DGLJD: Alpha and Beta must be greater than -1'
        call exitt
    ENDIF

    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 /= j) d(i,j) = ci/(cj*(z(i)-z(j)))
            IF ((i == j) .AND. (i /= 1) .AND. (i /= nz)) &
            d(i,j) = (alpha*(one+z(i))-beta*(one-z(i)))/ &
            (two*(one-z(i)**2))
            IF ((i == j) .AND. (i == 1)) &
            d(i,j) =  (eigval+alpha)/(two*(beta+two))
            IF ((i == j) .AND. (i == nz)) &
            d(i,j) = -(eigval+beta)/(two*(alpha+two))
            dt(j,i) = d(i,j)
    200 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 none
    integer, intent(in) :: nz, nzd 
    REAL(DP), intent(out) :: D(nzd,nzd),DT(nzd,nzd)
    REAL(DP), intent(in) :: Z(nz)

    integer :: n, i, j
    real(DP) :: fn, d0
    n  = nz-1
    IF (nz > nmax) THEN
        WRITE (6,*) 'Subroutine DGLL'
        WRITE (6,*) 'Maximum polynomial degree =',nmax
        WRITE (6,*) 'Polynomial degree         =',nz
    ENDIF
    IF (nz == 1) THEN
        d(1,1) = 0._dp
        RETURN
    ENDIF
    fn = (n)
    d0 = fn*(fn+1._dp)/4._dp
    DO 200 i=1,nz
        DO 200 j=1,nz
            d(i,j) = 0._dp
            IF  (i /= j) d(i,j) = pnleg(z(i),n)/ &
            (pnleg(z(j),n)*(z(i)-z(j)))
            IF ((i == j) .AND. (i == 1))  d(i,j) = -d0
            IF ((i == j) .AND. (i == nz)) d(i,j) =  d0
            dt(j,i) = d(i,j)
    200 END DO
    RETURN
    END SUBROUTINE dgll

    REAL(DP) 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 none
    integer, intent(in) :: i, nz
    REAL(DP), intent(in) :: z, ZGLL(nz)

    integer :: n
    real(DP) :: eps, dz, alfan

    eps = 1.e-5_dp
    dz = z - zgll(i)
    IF (abs(dz) < eps) THEN
        hgll = 1._dp
        RETURN
    ENDIF
    n = nz - 1
    alfan = (n)*((n)+1._dp)
    hgll = - (1._dp-z*z)*pndleg(z,n)/ &
    (alfan*pnleg(zgll(i),n)*(z-zgll(i)))
    RETURN
    END FUNCTION hgll

    REAL(DP) 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.
!---------------------------------------------------------------------
    implicit none
    integer, intent(in) :: i, nz
    REAL(DP), intent(in) :: z, ZGL(nz)

    real(DP) :: eps, dz
    integer :: n

    eps = 1.e-5_dp
    dz = z - zgl(i)
    IF (abs(dz) < eps) THEN
        hgl = 1._dp
        RETURN
    ENDIF
    n = nz-1
    hgl = pnleg(z,nz)/(pndleg(zgl(i),nz)*(z-zgl(i)))
    RETURN
    END FUNCTION hgl

    REAL(DP) 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.
!---------------------------------------------------------------------
    implicit none
    real(DP), intent(in) :: z
    integer,  intent(in) :: n

    real(DP) :: p1, p2, p3, fk
    integer :: k

    p1   = 1._dp
    IF (n == 0) THEN
        pnleg = p1
        RETURN
    ENDIF
    p2   = z
    p3   = p2
    DO 10 k = 1, n-1
        fk  = (k)
        p3  = ((2._dp*fk+1._dp)*z*p2 - fk*p1)/(fk+1._dp)
        p1  = p2
        p2  = p3
    10 END DO
    pnleg = p3
    if (n == 0) pnleg = 1._dp
    RETURN
    END FUNCTION pnleg

    REAL(DP) 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 none
    real(DP), intent(in) :: z
    integer,  intent(in) :: n

    real(DP) :: p1, p2, p3, p1d, p2d, p3d
    integer :: k, fk
    p1   = 1._dp
    p2   = z
    p1d  = 0._dp
    p2d  = 1._dp
    p3d  = 1._dp
    DO 10 k = 1, n-1
        fk  = (k)
        p3  = ((2._dp*fk+1._dp)*z*p2 - fk*p1)/(fk+1._dp)
        p3d = ((2._dp*fk+1._dp)*p2 + (2._dp*fk+1._dp)*z*p2d - fk*p1d)/(fk+1._dp)
        p1  = p2
        p2  = p3
        p1d = p2d
        p2d = p3d
    10 END DO
    pndleg = p3d
    IF (n == 0) pndleg = 0._dp
    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.
!-----------------------------------------------------------------------
    implicit none
    integer, intent(in) :: nd1, nd2, nzm1, nzm2
    REAL(DP), intent(out) :: D(nd2,nd1), DT(nd1,nd2)
    REAL(DP), intent(in)  :: ZM1(nd1), ZM2(nd2), IM12(nd2,nd1)

    integer :: ip, jq, nm1
    real(DP) :: zp, zq, eps

    IF (nzm1 == 1) THEN
        d(1,1) = 0._dp
        dt(1,1) = 0._dp
        RETURN
    ENDIF
    eps = 1.e-6_dp
    nm1 = nzm1-1
    DO 10 ip = 1, nzm2
        DO 10 jq = 1, nzm1
            zp = zm2(ip)
            zq = zm1(jq)
            IF ((abs(zp) < eps) .AND. (abs(zq) < eps)) THEN
                d(ip,jq) = 0._dp
            ELSE
                d(ip,jq) = (pnleg(zp,nm1)/pnleg(zq,nm1) &
                -im12(ip,jq))/(zp-zq)
            ENDIF
            dt(jq,ip) = d(ip,jq)
    10 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.
!-----------------------------------------------------------------------
    implicit none
    integer, intent(in) :: nd1, nd2, npgl, npg
    REAL(DP), intent(out) :: D(nd2,nd1), DT(nd1,nd2)
    REAL(DP), intent(in)  :: ZGL(nd1), ZG(nd2), IGLG(nd2,nd1)
    real(DP), intent(in)  :: alpha, beta

    integer, PARAMETER :: NDD = NMAX
    REAL(DP) ::  DD(ndd,ndd), DTD(ndd,ndd)
    REAL(DP) ::  ZGD(ndd), ZGLD(ndd), IGLGD(ndd,ndd)
    REAL(DP) ::  ALPHAD, BETAD
    integer :: i, j

    IF (npgl <= 1) THEN
        WRITE(6,*) 'DGLJGJ: Minimum number of Gauss-Lobatto points is 2'
        call exitt
    ENDIF
    IF (npgl > 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 <= -1._dp) .OR. (beta <= -1._dp)) THEN
        WRITE(6,*) 'DGLJGJ: Alpha and Beta must be greater than -1'
        call exitt
    ENDIF

    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 END DO
    DO 200 i=1,npgl
        zgld(i) = zgl(i)
    200 END DO
    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 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 none
    integer, intent(in) :: nd1, nd2, npgl, npg
    REAL(DP), intent(out) ::  D(nd2,nd1), DT(nd1,nd2)
    REAL(DP), intent(in) ::  ZGL(nd1), ZG(nd2), IGLG(nd2,nd1), ALPHA, BETA

    real(DP) :: eps, one, two, faci, facj, const, dz, dn, eigval
    real(DP) :: pj, pdj, pm1, pdm1, pm2, pdm2, pi, pdi
    integer :: i, j, ngl
    IF (npgl <= 1) THEN
        WRITE(6,*) 'DGLJGJD: Minimum number of Gauss-Lobatto points is 2'
        call exitt
    ENDIF
    IF ((alpha <= -1.) .OR. (beta <= -1.)) THEN
        WRITE(6,*) 'DGLJGJD: Alpha and Beta must be greater than -1'
        call exitt
    ENDIF

    eps    = 1.e-6_dp
    one    = 1._dp
    two    = 2._dp
    ngl    = npgl-1
    dn     = ((ngl))
    eigval = -dn*(dn+alpha+beta+one)

    DO 100 i=1,npg
        DO 100 j=1,npgl
            dz = abs(zg(i)-zgl(j))
            IF (dz < 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 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.
!--------------------------------------------------------------------
    implicit none
    integer, intent(in) :: nz1, nz2, nd1, nd2
    REAL(DP), intent(out) :: I12(nd2,nd1),IT12(nd1,nd2)
    REAL(DP), intent(in) :: Z1(nd1),Z2(nd2)

    integer :: i, j
    real(DP) :: zi

    IF (nz1 == 1) THEN
        i12(1,1) = 1._dp
        it12(1,1) = 1._dp
        RETURN
    ENDIF
    DO 10 i=1,nz2
        zi = z2(i)
        DO 10 j=1,nz1
            i12(i,j) = hgl(j,zi,z1,nz1)
            it12(j,i) = i12(i,j)
    10 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.
!--------------------------------------------------------------------
    implicit none
    integer, intent(in) :: nz1, nz2, nd1, nd2
    REAL(DP), intent(out) :: I12(nd2,nd1),IT12(nd1,nd2)
    REAL(DP), intent(in) :: Z1(nd1),Z2(nd2)

    integer :: i, j
    real(DP) :: zi

    IF (nz1 == 1) THEN
        i12(1,1) = 1._dp
        it12(1,1) = 1._dp
        RETURN
    ENDIF
    DO 10 i=1,nz2
        zi = z2(i)
        DO 10 j=1,nz1
            i12(i,j) = hgll(j,zi,z1,nz1)
            it12(j,i) = i12(i,j)
    10 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.
!--------------------------------------------------------------------
    implicit none
    integer, intent(in) :: nz1, nz2, nd1, nd2
    REAL(DP), intent(out) :: I12(nd2,nd1),IT12(nd1,nd2)
    REAL(DP), intent(in) :: Z1(nd1),Z2(nd2), ALPHA, BETA

    integer :: i, j
    real(DP) :: zi
    IF (nz1 == 1) THEN
        i12(1,1) = 1._dp
        it12(1,1) = 1._dp
        RETURN
    ENDIF
    DO 10 i=1,nz2
        zi = z2(i)
        DO 10 j=1,nz1
            i12(i,j) = hgj(j,zi,z1,nz1,alpha,beta)
            it12(j,i) = i12(i,j)
    10 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.
!--------------------------------------------------------------------
    implicit none
    integer, intent(in) :: nz1, nz2, nd1, nd2
    REAL(DP), intent(out) :: I12(nd2,nd1),IT12(nd1,nd2)
    REAL(DP), intent(in) :: Z1(nd1),Z2(nd2), alpha, beta

    integer :: i, j
    real(DP) :: zi
    IF (nz1 == 1) THEN
        i12(1,1) = 1._dp
        it12(1,1) = 1._dp
        RETURN
    ENDIF
    DO 10 i=1,nz2
        zi = z2(i)
        DO 10 j=1,nz1
            i12(i,j) = hglj(j,zi,z1,nz1,alpha,beta)
            it12(j,i) = i12(i,j)
    10 END DO
    RETURN
    END SUBROUTINE igljm

end module speclib