161 REAL(KIND=
rp) z(1),w(1), alpha, beta
164 CALL zwgj (z,w,np,alpha,beta)
178 REAL(KIND=
rp) z(1),w(1), alpha, beta
181 CALL zwglj (z,w,np,alpha,beta)
185 SUBROUTINE zwgj (Z,W,NP,ALPHA,BETA)
195 parameter(nzd = nmax)
196 REAL(KIND=
rp) zd(nzd),wd(nzd),alphad,betad
197 REAL(KIND=
rp) z(1),w(1),alpha,beta
200 IF (np.GT.npmax)
THEN
201 WRITE (6,*)
'Too large polynomial degree in ZWGJ'
202 WRITE (6,*)
'Maximum polynomial degree is',nmax
203 WRITE (6,*)
'Here NP=',np
208 CALL zwgjd (zd,wd,np,alphad,betad)
216 SUBROUTINE zwgjd (Z,W,NP,ALPHA,BETA)
225 IMPLICIT REAL(KIND=RP) (a-h,o-z)
226 REAL(KIND=
rp) z(1),w(1),alpha,beta
235 WRITE (6,*)
'ZWGJD: Minimum number of Gauss points is 1',np
238 IF ((alpha.LE.-one).OR.(beta.LE.-one))
THEN
239 WRITE (6,*)
'ZWGJD: Alpha and Beta must be greater than -1'
244 z(1) = (beta-alpha)/(apb+two)
250 CALL jacg (z,np,alpha,beta)
256 fac1 = dnp1+alpha+beta+one
259 fnorm =
pnormj(np1,alpha,beta)
260 rcoef = (fnorm*fac2*fac3)/(two*fac1*dnp2)
262 CALL jacobf (p,pd,pm1,pdm1,pm2,pdm2,np2,alpha,beta,z(i))
263 w(i) = -rcoef/(p*pdm1)
268 SUBROUTINE zwglj (Z,W,NP,ALPHA,BETA)
278 parameter(nzd = nmax)
279 REAL(KIND=
rp) zd(nzd),wd(nzd),alphad,betad
280 REAL(KIND=
rp) z(1),w(1),alpha,beta
283 IF (np.GT.npmax)
THEN
284 WRITE (6,*)
'Too large polynomial degree in ZWGLJ'
285 WRITE (6,*)
'Maximum polynomial degree is',nmax
286 WRITE (6,*)
'Here NP=',np
291 CALL zwgljd (zd,wd,np,alphad,betad)
308 IMPLICIT REAL(KIND=RP) (a-h,o-z)
309 REAL(KIND=
rp) z(np),w(np),alpha,beta
317 WRITE (6,*)
'ZWGLJD: Minimum number of Gauss-Lobatto points is 2'
318 WRITE (6,*)
'ZWGLJD: alpha,beta:',alpha,beta,np
321 IF ((alpha.LE.-one).OR.(beta.LE.-one))
THEN
322 WRITE (6,*)
'ZWGLJD: Alpha and Beta must be greater than -1'
329 CALL zwgjd (z(2),w(2),nm1,alpg,betg)
334 w(i) = w(i)/(one-z(i)**2)
336 CALL jacobf (p,pd,pm1,pdm1,pm2,pdm2,n,alpha,beta,z(1))
337 w(1) =
endw1(n,alpha,beta)/(two*pd)
338 CALL jacobf (p,pd,pm1,pdm1,pm2,pdm2,n,alpha,beta,z(np))
339 w(np) =
endw2(n,alpha,beta)/(two*pd)
345 IMPLICIT REAL(KIND=RP) (a-h,o-z)
346 REAL(kind=
rp) alpha,beta
358 f1 = f1*(apb+two)*two**(apb+two)/two
364 fint1 = fint1*two**(apb+two)
366 fint2 = fint2*two**(apb+three)
367 f2 = (-two*(beta+two)*fint1 + (apb+four)*fint2) &
377 a1 = -(two*(di+alpha)*(di+beta))/(abn*abnn*(abnn+one))
378 a2 = (two*(alpha-beta))/(abnn*(abnn+two))
379 a3 = (two*(abn+one))/((abnn+two)*(abnn+one))
380 f3 = -(a2*f2+a1*f1)/a3
389 IMPLICIT REAL(KIND=RP) (a-h,o-z)
390 REAL(kind=
rp) alpha,beta
402 f1 = f1*(apb+two)*two**(apb+two)/two
408 fint1 = fint1*two**(apb+two)
410 fint2 = fint2*two**(apb+three)
411 f2 = (two*(alpha+two)*fint1 - (apb+four)*fint2) &
421 a1 = -(two*(di+alpha)*(di+beta))/(abn*abnn*(abnn+one))
422 a2 = (two*(alpha-beta))/(abnn*(abnn+two))
423 a3 = (two*(abn+one))/((abnn+two)*(abnn+one))
424 f3 = -(a2*f2+a1*f1)/a3
433 IMPLICIT REAL(KIND=RP) (a-h,o-z)
442 IF (x.EQ.-half)
gammaf = -two*sqrt(pi)
443 IF (x.EQ. half)
gammaf = sqrt(pi)
444 IF (x.EQ. one )
gammaf = one
445 IF (x.EQ. two )
gammaf = one
446 IF (x.EQ. 1.5 )
gammaf = sqrt(pi)/2.
447 IF (x.EQ. 2.5)
gammaf = 1.5*sqrt(pi)/2.
448 IF (x.EQ. 3.5)
gammaf = 0.5*(2.5*(1.5*sqrt(pi)))
449 IF (x.EQ. 3. )
gammaf = 2.
450 IF (x.EQ. 4. )
gammaf = 6.
451 IF (x.EQ. 5. )
gammaf = 24.
452 IF (x.EQ. 6. )
gammaf = 120.
457 IMPLICIT REAL(KIND=RP) (a-h,o-z)
458 REAL(kind=
rp) alpha,beta
462 const = alpha+beta+one
466 pnormj = prod * two**const/(two*dn+const)
470 prod = prod/(two*(one+const)*
gammaf(const+one))
471 prod = prod*(one+alpha)*(two+alpha)
472 prod = prod*(one+beta)*(two+beta)
475 frac = (dindx+alpha)*(dindx+beta)/(dindx*(dindx+alpha+beta))
478 pnormj = prod * two**const/(two*dn+const)
482 SUBROUTINE jacg (XJAC,NP,ALPHA,BETA)
493 IMPLICIT REAL(KIND=RP) (a-h,o-z)
494 REAL(KIND=
rp) xjac(1)
499 dth = 4.*atan(one)/(2.*((n))+2.)
502 x = cos((2.*(((j))-1.)+1.)*dth)
504 x1 = cos((2.*(((j))-1.)+1.)*dth)
509 CALL jacobf (p,pd,pm1,pdm1,pm2,pdm2,np,alpha,beta,x)
513 recsum = recsum+1./(x-xjac(np-i+1))
515 delx = -p/(pd-recsum*p)
517 IF (abs(delx) .LT. eps)
GOTO 31
526 IF (xjac(j).LT.xmin)
THEN
540 SUBROUTINE jacobf (POLY,PDER,POLYM1,PDERM1,POLYM2,PDERM2,N,ALP,BET,X)
547 IMPLICIT REAL(KIND=RP) (a-h,o-z)
554 poly = (alp-bet+(apb+2.)*x)/2.
559 a1 = 2.*dk*(dk+apb)*(2.*dk+apb-2.)
560 a2 = (2.*dk+apb-1.)*(alp**2-bet**2)
562 a3 = b3*(b3+1.)*(b3+2.)
563 a4 = 2.*(dk+alp-1.)*(dk+bet-1.)*(2.*dk+apb)
564 polyn = ((a2+a3*x)*poly-a4*polyl)/a1
565 pdern = ((a2+a3*x)*pder-a4*pderl+a3*poly)/a1
580 REAL(kind=
rp)
FUNCTION hgj (II,Z,ZGJ,NP,ALPHA,BETA)
589 parameter(nzd = nmax)
590 REAL(kind=
rp) zd,zgjd(nzd),alphad,betad
591 REAL(kind=
rp) z,zgj(1),alpha,beta
593 IF (np.GT.npmax)
THEN
594 WRITE (6,*)
'Too large polynomial degree in HGJ'
595 WRITE (6,*)
'Maximum polynomial degree is',nmax
596 WRITE (6,*)
'Here NP=',np
605 hgj =
hgjd(ii,zd,zgjd,np,alphad,betad)
609 REAL(kind=
rp)
FUNCTION hgjd (II,Z,ZGJ,NP,ALPHA,BETA)
617 IMPLICIT REAL(KIND=RP) (a-h,o-z)
618 REAL(kind=
rp) z,zgj(1),alpha,beta
623 IF (abs(dz).LT.eps)
THEN
627 CALL jacobf (pzi,pdzi,pm1,pdm1,pm2,pdm2,np,alpha,beta,zi)
628 CALL jacobf (pz,pdz,pm1,pdm1,pm2,pdm2,np,alpha,beta,z)
629 hgjd = pz/(pdzi*(z-zi))
633 REAL(kind=
rp)
FUNCTION hglj (II,Z,ZGLJ,NP,ALPHA,BETA)
642 parameter(nzd = nmax)
643 REAL(kind=
rp) zd,zgljd(nzd),alphad,betad
644 REAL(kind=
rp) z,zglj(1),alpha,beta
646 IF (np.GT.npmax)
THEN
647 WRITE (6,*)
'Too large polynomial degree in HGLJ'
648 WRITE (6,*)
'Maximum polynomial degree is',nmax
649 WRITE (6,*)
'Here NP=',np
658 hglj =
hgljd(ii,zd,zgljd,np,alphad,betad)
662 REAL(kind=
rp)
FUNCTION hgljd (I,Z,ZGLJ,NP,ALPHA,BETA)
670 IMPLICIT REAL(KIND=RP) (a-h,o-z)
671 REAL(kind=
rp) z,zglj(1),alpha,beta
676 IF (abs(dz).LT.eps)
THEN
682 eigval = -dn*(dn+alpha+beta+one)
683 CALL jacobf (pi,pdi,pm1,pdm1,pm2,pdm2,n,alpha,beta,zi)
684 const = eigval*pi+alpha*(one+zi)*pdi-beta*(one-zi)*pdi
685 CALL jacobf (p,pd,pm1,pdm1,pm2,pdm2,n,alpha,beta,z)
686 hgljd = (one-z**2)*pd/(const*(z-zi))
690 SUBROUTINE dgj (D,DT,Z,NZ,NZD,ALPHA,BETA)
701 parameter(nzdd = nmax)
702 REAL(KIND=
rp) dd(nzdd,nzdd),dtd(nzdd,nzdd),zd(nzdd),alphad,betad
703 REAL(KIND=
rp) d(nzd,nzd),dt(nzd,nzd),z(1),alpha,beta
706 WRITE (6,*)
'DGJ: Minimum number of Gauss points is 1'
709 IF (nz .GT. nmax)
THEN
710 WRITE (6,*)
'Too large polynomial degree in DGJ'
711 WRITE (6,*)
'Maximum polynomial degree is',nmax
712 WRITE (6,*)
'Here Nz=',nz
715 IF ((alpha.LE.-1.).OR.(beta.LE.-1.))
THEN
716 WRITE (6,*)
'DGJ: Alpha and Beta must be greater than -1'
724 CALL dgjd (dd,dtd,zd,nz,nzdd,alphad,betad)
734 SUBROUTINE dgjd (D,DT,Z,NZ,NZD,ALPHA,BETA)
744 IMPLICIT REAL(KIND=RP) (a-h,o-z)
745 REAL(KIND=
rp) d(nzd,nzd),dt(nzd,nzd),z(1),alpha,beta
752 WRITE (6,*)
'DGJD: Minimum number of Gauss-Lobatto points is 2'
755 IF ((alpha.LE.-one).OR.(beta.LE.-one))
THEN
756 WRITE (6,*)
'DGJD: Alpha and Beta must be greater than -1'
762 CALL jacobf (pi,pdi,pm1,pdm1,pm2,pdm2,nz,alpha,beta,z(i))
763 CALL jacobf (pj,pdj,pm1,pdm1,pm2,pdm2,nz,alpha,beta,z(j))
764 IF (i.NE.j) d(i,j) = pdi/(pdj*(z(i)-z(j)))
765 IF (i.EQ.j) d(i,j) = ((alpha+beta+two)*z(i)+alpha-beta)/ &
773 SUBROUTINE dglj (D,DT,Z,NZ,NZD,ALPHA,BETA)
784 parameter(nzdd = nmax)
785 REAL(KIND=
rp) dd(nzdd,nzdd),dtd(nzdd,nzdd),zd(nzdd),alphad,betad
786 REAL(KIND=
rp) d(nzd,nzd),dt(nzd,nzd),z(1),alpha,beta
789 WRITE (6,*)
'DGLJ: Minimum number of Gauss-Lobatto points is 2'
792 IF (nz .GT. nmax)
THEN
793 WRITE (6,*)
'Too large polynomial degree in DGLJ'
794 WRITE (6,*)
'Maximum polynomial degree is',nmax
795 WRITE (6,*)
'Here NZ=',nz
798 IF ((alpha.LE.-1.).OR.(beta.LE.-1.))
THEN
799 WRITE (6,*)
'DGLJ: Alpha and Beta must be greater than -1'
807 CALL dgljd (dd,dtd,zd,nz,nzdd,alphad,betad)
817 SUBROUTINE dgljd (D,DT,Z,NZ,NZD,ALPHA,BETA)
827 IMPLICIT REAL(KIND=RP) (a-h,o-z)
828 REAL(KIND=
rp) d(nzd,nzd),dt(nzd,nzd),z(1),alpha,beta
833 eigval = -dn*(dn+alpha+beta+one)
836 WRITE (6,*)
'DGLJD: Minimum number of Gauss-Lobatto points is 2'
839 IF ((alpha.LE.-one).OR.(beta.LE.-one))
THEN
840 WRITE (6,*)
'DGLJD: Alpha and Beta must be greater than -1'
846 CALL jacobf (pi,pdi,pm1,pdm1,pm2,pdm2,n,alpha,beta,z(i))
847 CALL jacobf (pj,pdj,pm1,pdm1,pm2,pdm2,n,alpha,beta,z(j))
848 ci = eigval*pi-(beta*(one-z(i))-alpha*(one+z(i)))*pdi
849 cj = eigval*pj-(beta*(one-z(j))-alpha*(one+z(j)))*pdj
850 IF (i.NE.j) d(i,j) = ci/(cj*(z(i)-z(j)))
851 IF ((i.EQ.j).AND.(i.NE.1).AND.(i.NE.nz)) &
852 d(i,j) = (alpha*(one+z(i))-beta*(one-z(i)))/ &
854 IF ((i.EQ.j).AND.(i.EQ.1)) &
855 d(i,j) = (eigval+alpha)/(two*(beta+two))
856 IF ((i.EQ.j).AND.(i.EQ.nz)) &
857 d(i,j) = -(eigval+beta)/(two*(alpha+two))
864 SUBROUTINE dgll (D,DT,Z,NZ,NZD)
873 IMPLICIT REAL(KIND=RP) (a-h,o-z)
875 REAL(KIND=
rp) d(nzd,nzd),dt(nzd,nzd),z(1)
877 IF (nz .GT. nmax)
THEN
878 WRITE (6,*)
'Subroutine DGLL'
879 WRITE (6,*)
'Maximum polynomial degree =',nmax
880 WRITE (6,*)
'Polynomial degree =',nz
891 IF (i.NE.j) d(i,j) =
pnleg(z(i),n)/ &
892 (
pnleg(z(j),n)*(z(i)-z(j)))
893 IF ((i.EQ.j).AND.(i.EQ.1)) d(i,j) = -d0
894 IF ((i.EQ.j).AND.(i.EQ.nz)) d(i,j) = d0
901 REAL(kind=
rp)
FUNCTION hgll (I,Z,ZGLL,NZ)
908 IMPLICIT REAL(KIND=RP) (a-h,o-z)
909 REAL(kind=
rp) zgll(1), eps, dz, z
912 IF (abs(dz) .LT. eps)
THEN
919 (alfan*
pnleg(zgll(i),n)*(z-zgll(i)))
923 REAL(kind=
rp)
FUNCTION hgl (I,Z,ZGL,NZ)
930 REAL(kind=
rp) zgl(1), z, eps, dz
933 IF (abs(dz) .LT. eps)
THEN
954 IMPLICIT REAL(KIND=RP) (a-h,o-z)
955 REAL(kind=
rp) z, p1, p2, p3
956 IF(abs(z) .LT. 1.0e-25) z = 0.0
968 p3 = ((2.*fk+1.)*z*p2 - fk*p1)/(fk+1.)
973 if (n.eq.0)
pnleg = 1.
979 real(kind=
rp),
intent(inout):: l(1:n+1)
987 l(j) = ( (2*j-1) * x * l(j-1) - (j-1) * l(j-2) ) / j
999 IMPLICIT REAL(KIND=RP) (a-h,o-z)
1000 REAL(kind=
rp) p1, p2, p1d, p2d, p3d, z
1008 p3 = ((2.*fk+1.)*z*p2 - fk*p1)/(fk+1.)
1009 p3d = ((2.*fk+1.)*p2 + (2.*fk+1.)*z*p2d - fk*p1d)/(fk+1.)
1020 SUBROUTINE dgllgl (D,DT,ZM1,ZM2,IM12,NZM1,NZM2,ND1,ND2)
1032 REAL(KIND=
rp) d(nd2,nd1), dt(nd1,nd2), zm1(nd1), zm2(nd2), im12(nd2,nd1)
1033 REAL(KIND=
rp) eps, zp, zq
1045 IF ((abs(zp) .LT. eps).AND.(abs(zq) .LT. eps))
THEN
1049 -im12(ip,jq))/(zp-zq)
1051 dt(jq,ip) = d(ip,jq)
1057 SUBROUTINE dgljgj (D,DT,ZGL,ZG,IGLG,NPGL,NPG,ND1,ND2,ALPHA,BETA)
1070 REAL(KIND=
rp) d(nd2,nd1), dt(nd1,nd2), zgl(nd1), zg(nd2), iglg(nd2,nd1)
1072 parameter(ndd = nmax)
1073 REAL(KIND=
rp) dd(ndd,ndd), dtd(ndd,ndd)
1074 REAL(KIND=
rp) zgd(ndd), zgld(ndd), iglgd(ndd,ndd)
1075 REAL(KIND=
rp) alphad, betad
1078 WRITE(6,*)
'DGLJGJ: Minimum number of Gauss-Lobatto points is 2'
1081 IF (npgl.GT.nmax)
THEN
1082 WRITE(6,*)
'Polynomial degree too high in DGLJGJ'
1083 WRITE(6,*)
'Maximum polynomial degree is',nmax
1084 WRITE(6,*)
'Here NPGL=',npgl
1087 IF ((alpha.LE.-1.).OR.(beta.LE.-1.))
THEN
1088 WRITE(6,*)
'DGLJGJ: Alpha and Beta must be greater than -1'
1097 iglgd(i,j) = iglg(i,j)
1103 CALL dgljgjd (dd,dtd,zgld,zgd,iglgd,npgl,npg,ndd,ndd,alphad,betad)
1113 SUBROUTINE dgljgjd (D,DT,ZGL,ZG,IGLG,NPGL,NPG,ND1,ND2,ALPHA,BETA)
1126 IMPLICIT REAL(KIND=RP) (a-h,o-z)
1127 REAL(KIND=
rp) d(nd2,nd1), dt(nd1,nd2), zgl(nd1), zg(nd2)
1128 REAL(KIND=
rp) iglg(nd2,nd1), alpha, beta
1131 WRITE(6,*)
'DGLJGJD: Minimum number of Gauss-Lobatto points is 2'
1134 IF ((alpha.LE.-1.).OR.(beta.LE.-1.))
THEN
1135 WRITE(6,*)
'DGLJGJD: Alpha and Beta must be greater than -1'
1144 eigval = -dn*(dn+alpha+beta+one)
1148 dz = abs(zg(i)-zgl(j))
1150 d(i,j) = (alpha*(one+zg(i))-beta*(one-zg(i)))/ &
1151 (two*(one-zg(i)**2))
1153 CALL jacobf (pi,pdi,pm1,pdm1,pm2,pdm2,ngl,alpha,beta,zg(i))
1154 CALL jacobf (pj,pdj,pm1,pdm1,pm2,pdm2,ngl,alpha,beta,zgl(j))
1155 faci = alpha*(one+zg(i))-beta*(one-zg(i))
1156 facj = alpha*(one+zgl(j))-beta*(one-zgl(j))
1157 const = eigval*pj+facj*pdj
1158 d(i,j) = ((eigval*pi+faci*pdi)*(zg(i)-zgl(j)) &
1159 -(one-zg(i)**2)*pdi)/(const*(zg(i)-zgl(j))**2)
1167 SUBROUTINE iglm (I12,IT12,Z1,Z2,NZ1,NZ2,ND1,ND2)
1177 REAL(KIND=
rp) i12(nd2,nd1),it12(nd1,nd2),z1(nd1),z2(nd2), zi
1178 IF (nz1 .EQ. 1)
THEN
1186 i12(i,j) =
hgl(j,zi,z1,nz1)
1187 it12(j,i) = i12(i,j)
1193 SUBROUTINE igllm (I12,IT12,Z1,Z2,NZ1,NZ2,ND1,ND2)
1203 REAL(KIND=
rp) i12(nd2,nd1),it12(nd1,nd2),z1(nd1),z2(nd2),zi
1204 IF (nz1 .EQ. 1)
THEN
1212 i12(i,j) =
hgll(j,zi,z1,nz1)
1213 it12(j,i) = i12(i,j)
1217 end subroutine igllm
1219 SUBROUTINE igjm (I12,IT12,Z1,Z2,NZ1,NZ2,ND1,ND2,ALPHA,BETA)
1230 REAL(KIND=
rp) i12(nd2,nd1),it12(nd1,nd2),z1(nd1),z2(nd2),zi,alpha,beta
1231 IF (nz1 .EQ. 1)
THEN
1239 i12(i,j) =
hgj(j,zi,z1,nz1,alpha,beta)
1240 it12(j,i) = i12(i,j)
1246 SUBROUTINE igljm (I12,IT12,Z1,Z2,NZ1,NZ2,ND1,ND2,ALPHA,BETA)
1257 REAL(KIND=
rp) i12(nd2,nd1),it12(nd1,nd2),z1(nd1),z2(nd2),zi,alpha,beta
1258 IF (nz1 .EQ. 1)
THEN
1266 i12(i,j) =
hglj(j,zi,z1,nz1,alpha,beta)
1267 it12(j,i) = i12(i,j)
1271 end subroutine igljm
integer, parameter, public rp
Global precision used in computations.
LIBRARY ROUTINES FOR SPECTRAL METHODS.
real(kind=rp) function pnormj(N, ALPHA, BETA)
subroutine zwgll(Z, W, NP)
real(kind=rp) function pnleg(Z, N)
real(kind=rp) function hgj(II, Z, ZGJ, NP, ALPHA, BETA)
subroutine dgljgj(D, DT, ZGL, ZG, IGLG, NPGL, NPG, ND1, ND2, ALPHA, BETA)
subroutine zwgj(Z, W, NP, ALPHA, BETA)
subroutine legendre_poly(L, x, N)
subroutine dgj(D, DT, Z, NZ, NZD, ALPHA, BETA)
subroutine jacg(XJAC, NP, ALPHA, BETA)
subroutine zwgjd(Z, W, NP, ALPHA, BETA)
subroutine dgll(D, DT, Z, NZ, NZD)
real(kind=rp) function hglj(II, Z, ZGLJ, NP, ALPHA, BETA)
subroutine zwglj(Z, W, NP, ALPHA, BETA)
subroutine dglj(D, DT, Z, NZ, NZD, ALPHA, BETA)
real(kind=rp) function pndleg(Z, N)
subroutine igjm(I12, IT12, Z1, Z2, NZ1, NZ2, ND1, ND2, ALPHA, BETA)
real(kind=rp) function hgl(I, Z, ZGL, NZ)
real(kind=rp) function gammaf(X)
subroutine dgljgjd(D, DT, ZGL, ZG, IGLG, NPGL, NPG, ND1, ND2, ALPHA, BETA)
subroutine igljm(I12, IT12, Z1, Z2, NZ1, NZ2, ND1, ND2, ALPHA, BETA)
real(kind=rp) function hgjd(II, Z, ZGJ, NP, ALPHA, BETA)
subroutine iglm(I12, IT12, Z1, Z2, NZ1, NZ2, ND1, ND2)
subroutine zwgl(Z, W, NP)
Generate NP Gauss Legendre points Z and weights W associated with Jacobi polynomial ....
real(kind=rp) function hgljd(I, Z, ZGLJ, NP, ALPHA, BETA)
real(kind=rp) function endw1(N, ALPHA, BETA)
subroutine dgljd(D, DT, Z, NZ, NZD, ALPHA, BETA)
real(kind=rp) function hgll(I, Z, ZGLL, NZ)
real(kind=rp) function endw2(N, ALPHA, BETA)
subroutine zwgljd(Z, W, NP, ALPHA, BETA)
subroutine dgjd(D, DT, Z, NZ, NZD, ALPHA, BETA)
subroutine jacobf(POLY, PDER, POLYM1, PDERM1, POLYM2, PDERM2, N, ALP, BET, X)
subroutine dgllgl(D, DT, ZM1, ZM2, IM12, NZM1, NZM2, ND1, ND2)
subroutine igllm(I12, IT12, Z1, Z2, NZ1, NZ2, ND1, ND2)