DECOMPOSEREP Subroutine

    SUBROUTINE DECOMPOSEREP(CHARSIN, IDECOMP)
        IMPLICIT NONE
        TYPE(Symmetry) IDECOMP
        complex(dp) CHARS(NROT), CHARSIN(NROT), TOT
        real(dp) CNORM
        INTEGER I, J
        real(dp) NORM, DIFF
        character(*), parameter :: this_routine = 'DECOMPOSEREP'
        if (TAbelian) then
            ! We shouldn't be here!  Using symmetry "quantum" numbers
            ! rather than irreps.
            call stop_all(this_routine, "Should not be decomposing irreps with Abelian sym")
        end if
        IDECOMP%s = 0
        CALL DCOPY(NROT * 2, CHARSIN, 1, CHARS, 1)
!         write(stdout,*) "Decompose Rep"
!         CALL writechars(stdout,CHARS,NROT,"REP   ")
!,. First check norm of this state
        CNORM = 0
        DO J = 1, NROT
            CNORM = CNORM + real(CHARS(J) * CHARS(J), dp)
        end do
        DO I = 1, NSYM
            TOT = 0
!            CALL writechars(stdout,IRREPCHARS(1,I),NROT,"IR")
!            CALL writechars(stdout,CHARS(1),NROT,"CH")
            DO J = 1, NROT
                TOT = TOT + CONJG(IRREPCHARS(J, I)) * CHARS(J)
            end do
!            write(stdout,*) I,TOT
            IF (abs(TOT) > 1.0e-12_dp) THEN
!   Calculate the normalization of the state I which matches (if it's an irrep, this will be 1)
                NORM = 0
                DO J = 1, NROT
                    NORM = NORM + real(CONJG(IRREPCHARS(J, I)) * IRREPCHARS(J, I), dp)
                end do
!               write(stdout,*) "IRREP ",I,(TOT+0.0_dp)/NORM
                DIFF = ABS(TOT - NINT(ABS(TOT / NORM)) * NORM)
                IF (DIFF > 1.0e-2_dp) THEN
                    write(stdout, *) 'Symmetry decomposition not complete'
                    CALL writechars(stdout, IRREPCHARS(1, I), NROT, "IRREP ")
                    CALL writechars(stdout, CHARS, NROT, "CHARS ")
                    write(stdout, *) "Dot product: ", (TOT + 0.0_dp) / NORM, TOT, NORM
                    call stop_all(this_routine, 'Incomplete symmetry decomposition')
!   The given representation CHARS has fewer irreps in it than the one in IRREPCHARS, and is an irrep
!   Hurrah!  Remove it from the one in IRREPCHARS, and keep on going)
                else if (ABS(TOT) > 1.0e-2_dp) THEN
!   We've found an (ir)rep which is wholly in CHARS
                    IDECOMP%s = IBSET(IDECOMP%s, I - 1)
                    CNORM = 0
!                  write(stdout,*) I,DIFF,TOT,TOT/NORM
                    DO J = 1, NROT
                        CHARS(J) = CHARS(J) - (IRREPCHARS(J, I) * TOT) / NORM
                        CNORM = CNORM + real(CONJG(CHARS(J)) * CHARS(J), dp)
                    end do
!                  CALL writechars(stdout,IRREPCHARS(1,I),NROT,"DIRREP")
!                  CALL writechars(stdout,CHARS,NROT,"DCHARS")
                end if
            end if
        end do
    END SUBROUTINE DECOMPOSEREP