Normalized legendre and quadrature basis for discontinous Galerkin method

Question

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.
Solution:

The standard way of computing tensor products is $$ u_{h}(x,y) = \sum_{k = 0}^{N} \sum_{l = 0}^{M} u_{ij}(t) \phi_{k}(x) \phi_{l}(y) (1), $$ The tensor product formula adapted in the last two papers are $$ U_{h}(x,y) = \sum_{k=0}^{\frac{1}{2}(k + 1)(k + 2)} u_{ij}^{k}(t) \phi_{k}(x,y) (2) $$ Is there a reason to use formula (1) over formula (2)? Mathematically, they seem both correct.

Answer 1

Per the regulations of Stack Exchange. I have included the "answer" to my question.

So I managed to fix my problem. First, off I was missing a factor of 0.5 when calculating the surface integral. Additionally, I am using the following normalization $\phi(x) = \sqrt{2m + 1} L(\xi)$.