BoundaryFlux_DGModel3D_t Subroutine

  subroutine BoundaryFlux_DGModel3D_t(this)
    ! this method uses an linear upwind solver for the
    ! advective flux and the bassi-rebay method for the
    ! diffusive fluxes
    implicit none
    class(DGModel3D_t),intent(inout) :: this
    ! Local
    integer :: i,j,k,iel
    real(prec) :: sL(1:this%nvar),sR(1:this%nvar)
    real(prec) :: dsdx(1:this%nvar,1:3)
    real(prec) :: nhat(1:3),nmag

    do concurrent(i=1:this%solution%N+1,j=1:this%solution%N+1, &
                  k=1:6,iel=1:this%mesh%nElem)
      ! Get the boundary normals on cell edges from the mesh geometry
      nhat = this%geometry%nHat%boundary(i,j,k,iEl,1,1:3)
      sL = this%solution%boundary(i,j,k,iel,1:this%nvar) ! interior solution
      sR = this%solution%extboundary(i,j,k,iel,1:this%nvar) ! exterior solution
      dsdx = this%solutiongradient%avgboundary(i,j,k,iel,1:this%nvar,1:3)
      nmag = this%geometry%nScale%boundary(i,j,k,iEl,1)

      this%flux%boundaryNormal(i,j,k,iEl,1:this%nvar) = this%riemannflux3d(sL,sR,dsdx,nhat)*nmag

    enddo

  endsubroutine BoundaryFlux_DGModel3D_t