gll.f90

Spectral Element Method for wave propagation and rupture dynamics.

   else
      x1 = cos((2.d0*(dble(j)-1.d0)+1.d0)*dth)
      x2 = xlast
      x  = (x1+x2)/2.d0
   endif
   do 30 k=1,kstop
      call jacobf (p,pd,pm1,pdm1,pm2,pdm2,np,alpha,beta,x)
      recsum = 0.d0
      jm = j-1
      do 29 i=1,jm
         recsum = recsum+1.d0/(x-xjac(np-i+1))
 29         continue
      delx = -p/(pd-recsum*p)
      x    = x+delx
      if (abs(delx) < eps) goto 31
 30      continue
 31      continue
   xjac(np-j+1) = x
   xlast        = x
 40   continue
  do 200 i=1,np
   xmin = 2.d0
   do 100 j=i,np
      if (xjac(j) < xmin) then
         xmin = xjac(j)
         jmin = j
      endif
 100     continue
   if (jmin /= 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)
!
!=======================================================================
!
!     J a c o b f : Compute the Jacobi polynomial and its derivative
!     -----------   of degree n at x.
!
!=======================================================================
  

  double precision poly,pder,polym1,pderm1,polym2,pderm2,alp,bet,x
  integer n

  double precision apb,polyl,pderl,dk,a1,a2,b3,a3,a4,polyn,pdern,psave,pdsave
  integer k

  apb  = alp+bet
  poly = 1.d0
  pder = 0.d0
  psave = 0.d0
  pdsave = 0.d0

  if (n  ==  0) return

  polyl = poly
  pderl = pder
  poly  = (alp-bet+(apb+2.d0)*x)/2.d0
  pder  = (apb+2.d0)/2.d0
  if (n  ==  1) return
  do 20 k=2,n
   dk = dble(k)
   a1 = 2.d0*dk*(dk+apb)*(2.d0*dk+apb-2.d0)
   a2 = (2.d0*dk+apb-1.d0)*(alp**2-bet**2)
   b3 = (2.d0*dk+apb-2.d0)
   a3 = b3*(b3+1.d0)*(b3+2.d0)
   a4 = 2.d0*(dk+alp-1.d0)*(dk+bet-1.d0)*(2.d0*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
!
!=====================================================================

  double precision FUNCTION PNDLEG (Z,N)
!-------------------------------------------------------------
!
!     Compute the derivative of the Nth order Legendre polynomial at Z.
!     Based on the recursion formula for the Legendre polynomials.
!
!-------------------------------------------------------------
  

  double precision z
  integer n

  double precision P1,P2,P1D,P2D,P3D,FK,P3
  integer k

  P1   = 1.d0
  P2   = Z
  P1D  = 0.d0
  P2D  = 1.d0
  P3D  = 1.d0
  DO 10 K = 1, N-1
   FK  = dble(K)
   P3  = ((2.d0*FK+1.d0)*Z*P2 - FK*P1)/(FK+1.d0)
   P3D = ((2.d0*FK+1.d0)*P2 + (2.d0*FK+1.d0)*Z*P2D - FK*P1D) &
                          /(FK+1.d0)
   P1  = P2
   P2  = P3
   P1D = P2D
   P2D = P3D
 10   CONTINUE
  PNDLEG = P3D
  RETURN
  end function pndleg
!
!=====================================================================

  double precision FUNCTION PNLEG (Z,N)
!-------------------------------------------------------------
!
!     Compute the value of the Nth order Legendre polynomial at Z.
!     Based on the recursion formula for the Legendre polynomials.
!
!-------------------------------------------------------------
  

  double precision z
  integer n

  double precision P1,P2,P3,FK
  integer k

  P1   = 1.d0
  P2   = Z
  P3   = P2
  DO 10 K = 1, N-1
   FK  = dble(K)
   P3  = ((2.d0*FK+1.d0)*Z*P2 - FK*P1)/(FK+1.d0)
   P1  = P2
   P2  = P3
 10   CONTINUE
  PNLEG = P3
  RETURN
  end function pnleg
!
!=====================================================================

  double precision function pnormj (n,alpha,beta)
!
!=======================================================================
!
!     P n o r m j
!     -----------
!
!=======================================================================
!
  

  double precision alpha,beta
  integer n

  double precision one,two,dn,const,prod,dindx,frac
  integer i

  one   = 1.d0
  two   = 2.d0
  dn    = dble(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 = dble(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 zwgjd(z,w,np,alpha,beta)
!
!=======================================================================
!
!     Z w g j d : Generate np Gauss-Jacobi points and weights
!                 associated with Jacobi polynomial of degree n = np-1
!
!=======================================================================
!
!     Note : Coefficients alpha and beta must be greater than -1.
!     ----
!
!=======================================================================
!
  

  double precision, parameter :: zero=0.d0,one=1.d0,two=2.d0

  integer np
  double precision z(np),w(np)
  double precision alpha,beta

  integer n,np1,np2,i
  double precision p,pd,pm1,pdm1,pm2,pdm2
  double precision apb,dnp1,dnp2,fac1,fac2,fac3,fnorm,rcoef
!
!-----------------------------------------------------------------------
!
  pd = zero
  pm1 = zero
  pm2 = zero
  pdm1 = zero
  pdm2 = zero

  n    = np-1
  apb  = alpha+beta
  p    = zero
  pdm1 = zero

  if (np <= 0) stop 'Minimum number of Gauss points is 1'

  if ((alpha <= -one).or.(beta <= -one)) &
    stop 'Alpha and Beta must be greater than -1'

  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  = dble(np1)
  dnp2  = dble(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 i=1,np
    call jacobf (p,pd,pm1,pdm1,pm2,pdm2,np2,alpha,beta,z(i))
    w(i) = -rcoef/(p*pdm1)
  enddo

  return
  end subroutine zwgjd
!
!=====================================================================

  subroutine zwgljd (z,w,np,alpha,beta)
!
!=======================================================================
!
!     Z w g l j d : Generate np Gauss-Lobatto-Jacobi points and the
!     -----------   weights associated with Jacobi polynomials of degree
!                   n = np-1.
!
!=======================================================================
!
!     Note : alpha and beta coefficients must be greater than -1.
!     ----
!            Legendre polynomials are special case of Jacobi polynomials
!            just by setting alpha and beta to 0.
!
!=======================================================================
!
  

  double precision, parameter :: zero=0.d0,one=1.d0,two=2.d0

  integer np
  double precision alpha,beta
  double precision z(np), w(np)

  integer n,nm1,i
  double precision p,pd,pm1,pdm1,pm2,pdm2
  double precision alpg,betg
!
!-----------------------------------------------------------------------
!
  p = zero
  pm1 = zero
  pm2 = zero
  pdm1 = zero
  pdm2 = zero

  n   = np-1
  nm1 = n-1
  pd  = zero

  if (np <= 1) stop 'Minimum number of Gauss-Lobatto points is 2'

  if ((alpha <= -one).or.(beta <= -one)) &
    stop 'Alpha and Beta must be greater than -1'

  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  110 i=2,np-1
   w(i) = w(i)/(one-z(i)**2)
  110 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

end module gll