!--------------------------------------------------------------------------
   !---------------------------- div ----------------------------------
   !$OMP DO SCHEDULE(STATIC)
   DO ij = 1, primal_num
      DO l = 1, llm
         div_ij=0.
         DO iedge = 1, primal_deg(ij)
            edge = primal_edge(iedge,ij)
            div_ij = div_ij + primal_ne(iedge,ij)*le_de(edge)*u(l,edge)
         END DO
         divu(l,ij) = div_ij / Ai(ij)
      END DO
   END DO
   !$OMP END DO
   !---------------------------- div ----------------------------------
   !--------------------------------------------------------------------------