! An arbitrary polynomial of degree NP-1, using the coefficients
      ! in COEF and the point X to evaluate
      !  RES = COEFS(1) + COEFS(2)*X + COEFS(3)*X**2 + COEFS(4)*X**3 + ...
      REAL FUNCTION POLY(np, coefs, x) !RESULT (RES)
        IMPLICIT NONE
        INTEGER np, i
        REAL coefs(np)
        REAL x, xe, res, poly
        res = 0
        xe = 1
        DO i=1,np
          res = res + coefs(i)*xe
          xe = xe*x
        END DO
        poly = res
      END FUNCTION POLY

      !        Generated by TAPENADE     (INRIA, Ecuador team)
      !  Tapenade 3.16 (develop) - 17 Dec 2020 10:03
      !
      !  Differentiation of poly in forward (tangent) mode:
      !   variations   of useful results: res
      !   with respect to varying inputs: x
      !   RW status of diff variables: res:out x:in
      REAL FUNCTION POLY_D(np, coefs, x, xd, res) !RESULT (RESD)
        IMPLICIT NONE
        INTEGER  np, i
        REAL  coefs(np)
        REAL  x, xe, res
        REAL  xd, xed, resd, poly_d
        res = 0
        xe = 1
        resd = 0.0
        xed = 0.0
        DO i=1,np
          resd = resd + coefs(i)*xed
          res = res + coefs(i)*xe
          xed = x*xed + xe*xd
          xe = xe*x
        END DO
        poly_d = resd
      END FUNCTION POLY_D

      ! This is the function that we want to integrate. It is actually
      ! defined as the derivative of POLY, which itself is an arbitrary
      ! polynomial of degree NP-1. POLYD is thus a polynomial of degree
      ! NP-2. The function is just a thin wrapper around POLY_D, which
      ! is generated by applying the Automatic Differentiation tool
      ! Tapenade to the polynomial function.
      REAL FUNCTION POLYD(np, coefs, x) !RESULT(RESD)
      IMPLICIT NONE
      INTEGER  np
      REAL  coefs(np)
      REAL  x, res
      REAL  xd, resd, polyd
      xd = 1.0
      resd = POLY_D(np, coefs, x, xd, res)
      polyd = resd
      END FUNCTION

      SUBROUTINE TEST_ZWGLL(NP,DEG)
      IMPLICIT NONE
      INTEGER NP, I, DEG
      REAL RX, RY, RZ, COEF(DEG)
      REAL Z(NP), W(NP), COEFS(DEG), COEFS2(DEG), RES, RESD
      REAL QUADA, QUADB, ONE, NEGONE, DIFF, MINDIFF, TMP
      
      DO I=1, 2*NP, 2
          TMP = 2.0*I*ATAN(1)/NP
!$CVL $ASSUME(COS(TMP) + COS(4*ATAN(1) - TMP) == 0.0);
      END DO
!$CVL $ASSUME(2.0*COS(ATAN(1))*COS(ATAN(1)) == 1.0);
!$CVL $ASSUME(4.0*COS(10.0/3.0*ATAN(1))*COS(10.0/3.0*ATAN(1)) == 3.0)
!$CVL $ASSUME(SIN(4) .NE. 0.0)
!$CVL $ASSUME(SIN(6) .NE. 0.0)
      
      ! COEFS should be a symbolic input. It is currently
      ! set to some arbitrary values.
      DO i=1,DEG
        COEFS(i) = sin(i*1.0)
        COEFS2(i) = COEFS(i)
      END DO
      ! Compute the reference solution by using analytical integration.
      ONE = 1.0
      NEGONE = -1.0
      QUADA = POLY(DEG,COEFS,ONE)
      QUADA = QUADA  - POLY(DEG,COEFS2,NEGONE)
      ! Compute Gauss quadrature points/weights
      CALL ZWGLL(Z,W,NP)
      ! Compute the Integral using Gauss quadrature points/weights
      QUADB = 0.0
      DO i=1,NP
        QUADB = QUADB + W(i) * POLYD(DEG,COEFS,Z(i))
      END DO
      ! Compare results and print some diagnostics.
      DIFF = QUADA - QUADB
      MINDIFF = 0.0
      PRINT *, "N Points ", NP, "Degree ", DEG
      PRINT *, "Ref soln ", QUADA
      PRINT *, "Quadrature ", QUADB
      IF(DEG .le. 2*NP-1) THEN
        PRINT *, "Expected error: ZERO"
!$CVL $ASSERT(DIFF .EQ. MINDIFF)
      ELSE
        PRINT *, "Expected error: NON-ZERO"
!$CVL $ASSERT(DIFF .NE. MINDIFF)
      END IF
      !WRITE(*,*) "Points ", Z(1)
      !WRITE(*,*) "Weights", W(1)
      PRINT *, "\n"
      END SUBROUTINE

      PROGRAM DRIVER
      IMPLICIT NONE
      INTEGER NP, DEG
      ! The printed error should be zero as long as DEG <= (2*NP)-1.
      ! In the below examples, we have:
      NP = 2
      DEG = 2
      CALL TEST_ZWGLL(NP,DEG) ! Should be correct
      DEG = 3
      CALL TEST_ZWGLL(NP,DEG) ! Should be correct
      DEG = 4
      CALL TEST_ZWGLL(NP,DEG) ! Should be incorrect
      NP = 3
      CALL TEST_ZWGLL(NP,DEG) ! Should be correct
      DEG = 5
      CALL TEST_ZWGLL(NP,DEG) ! Should be correct
      DEG = 6
      CALL TEST_ZWGLL(NP,DEG) ! Should be incorrect
      END