I've successfully implemented a 1D DG code with non-normalized Legendre basis and I've now moved onto developing a 2D code using tensor products.  For my 2D code I've chosen to have normalized Legendre polynomials,i.e.

$$
L_{0}(x) = \sqrt{\frac{1}{2}}\\
L_{1}(x) = \sqrt{\frac{3}{2}} x.
$$

My question is:  If we normalize the Legendre polynomials, do I also have to normalized my Gauss-Legendre weights as well?  Or do we only normalized the polynomials in terms of the basis?

    subroutine GaussQuad (xq,wq,n)
    use parameters 
    implicit none 
    integer :: n
    real(kind=8),dimension(n) :: xq,wq

    integer :: i,iter
    real(kind=8) :: xx
    real(kind=8) :: legendre,dlegendre

    do i=1,n
     xx = cos(dpi*(i - 0.25d0)/(n + 0.5d0))

     do iter=1,500
        xx = xx - legendre(xx,n)/dlegendre(xx,n)
     end do

     xq(i) = xx
     wq(i) = (2.0*dble(n) + 1.0)*2.0d0/((1.0d0-xx**2.0)*dlegendre(xx,n)**2.0)
    end do
    end subroutine GaussQuad


In order for my 2D code to work I need to use the normalized Legendre polynomials in the Gauss quadrature routine along with the $2n + 1$ normalization on the weights (see how its been added to wq(i)).  However, I only got this to work due to an ad hoc guess.  I would like to avoid this as I don't personally understand why this is required for my solver to work.

EDIT 4/1/2020:  

1. In my code, I do all my operations in the reference element.  In fact, double checked and my Gauss-Legendre points values are within $[-1,+1]$ and my weights also sum to 2.  
2. To map from physical space to the reference space, I use the following
$$
X(\eta) = \frac{x_{i} - x_{i-1}}{2} \eta + \frac{x_{i} + x_{i-1}}{2}\\
Y(\zeta) = \frac{y_{j} - y_{j-1}}{2} \zeta + \frac{y_{j} + y_{j-1}}{2}
$$
3.  To integrate the volume fluxes I use a tensor product Gauss-Legendre scheme with $M^{2}$ points

        do ix=1,nx
         do iy=1,ny
          do i=1,mx
           do j=1,my
            do inode=1,mx
              do jnode=1,my
               call Flux(un(:,ix,iy,inode,jnode),FFlux,GFlux)

               flux_vol1(:,ix,iy,i,j) = flux_vol1(:,ix,iy,i,j) + &
               & 0.5*FFlux(:)* &
               & dlegendre(xg(inode),i-1)*& 
               & legendre(xg(jnode),j-1)*&
               & wg(inode)*&
               & wg(jnode)
            
               flux_vol2(:,ix,iy,i,j) = flux_vol2(:,ix,iy,i,j) + &
               & 0.5*GFlux(:)* &
               & dlegendre(xg(jnode),j-1)*&
               & legendre(xg(inode),i-1)*&
               & wg(inode)* &
               & wg(jnode)
             end do
            end do
           end do
          end do
         end do
        end do

Above is a snippet in how I perform the tensor product quadrature.  I'm a little confused on why removing the normalizing of the Gauss quadrature is causing such severe results.  Below are images of when I remove normalization of the weights and Legendre polynomials (but normalize the basis).
[![No normalizaion of wieghts][1]][1]
[![Normalized weights][2]][2]


  [1]: https://i.sstatic.net/Qc9k5.png
  [2]: https://i.sstatic.net/ISyut.png