!--------------------------------------------------------------------------
   !---------------------------- caldyn_slow_hydro ----------------------------------
   !$OMP DO SCHEDULE(STATIC)
   DO edge = 1, edge_num
      ij_left = left(edge)
      ij_right = right(edge)
      !DIR$ SIMD
      DO l = 1, llm
         uu = .5*(rhodz(l,ij_left)+rhodz(l,ij_right))*u(l,edge)
         hflux(l,edge) = uu*le_de(edge)
      END DO
   END DO
   !$OMP END DO
   !$OMP DO SCHEDULE(STATIC)
   DO ij = 1, primal_num
      ! this VLOOP iterates over primal cell edges
      SELECT CASE(primal_deg(ij))
      CASE(4)
         edge1 = primal_edge(1,ij)
         edge2 = primal_edge(2,ij)
         edge3 = primal_edge(3,ij)
         edge4 = primal_edge(4,ij)
         le_de1 = le_de(edge1)
         le_de2 = le_de(edge2)
         le_de3 = le_de(edge3)
         le_de4 = le_de(edge4)
         !DIR$ SIMD
         DO l = 1, llm
            ke=0.d0
            ke = ke + le_de1*u(l,edge1)**2
            ke = ke + le_de2*u(l,edge2)**2
            ke = ke + le_de3*u(l,edge3)**2
            ke = ke + le_de4*u(l,edge4)**2
            BERNI(l,ij)=ke*(.25/Ai(ij))
         END DO
      CASE(5)
         edge1 = primal_edge(1,ij)
         edge2 = primal_edge(2,ij)
         edge3 = primal_edge(3,ij)
         edge4 = primal_edge(4,ij)
         edge5 = primal_edge(5,ij)
         le_de1 = le_de(edge1)
         le_de2 = le_de(edge2)
         le_de3 = le_de(edge3)
         le_de4 = le_de(edge4)
         le_de5 = le_de(edge5)
         !DIR$ SIMD
         DO l = 1, llm
            ke=0.d0
            ke = ke + le_de1*u(l,edge1)**2
            ke = ke + le_de2*u(l,edge2)**2
            ke = ke + le_de3*u(l,edge3)**2
            ke = ke + le_de4*u(l,edge4)**2
            ke = ke + le_de5*u(l,edge5)**2
            BERNI(l,ij)=ke*(.25/Ai(ij))
         END DO
      CASE(6)
         edge1 = primal_edge(1,ij)
         edge2 = primal_edge(2,ij)
         edge3 = primal_edge(3,ij)
         edge4 = primal_edge(4,ij)
         edge5 = primal_edge(5,ij)
         edge6 = primal_edge(6,ij)
         le_de1 = le_de(edge1)
         le_de2 = le_de(edge2)
         le_de3 = le_de(edge3)
         le_de4 = le_de(edge4)
         le_de5 = le_de(edge5)
         le_de6 = le_de(edge6)
         !DIR$ SIMD
         DO l = 1, llm
            ke=0.d0
            ke = ke + le_de1*u(l,edge1)**2
            ke = ke + le_de2*u(l,edge2)**2
            ke = ke + le_de3*u(l,edge3)**2
            ke = ke + le_de4*u(l,edge4)**2
            ke = ke + le_de5*u(l,edge5)**2
            ke = ke + le_de6*u(l,edge6)**2
            BERNI(l,ij)=ke*(.25/Ai(ij))
         END DO
      CASE DEFAULT
         !DIR$ SIMD
         DO l = 1, llm
            ke=0.d0
            DO iedge = 1, primal_deg(ij)
               edge = primal_edge(iedge,ij)
               ke = ke + le_de(edge)*u(l,edge)**2
            END DO
            BERNI(l,ij)=ke*(.25/Ai(ij))
         END DO
      END SELECT
   END DO
   !$OMP END DO
   IF(zero) THEN
      !$OMP DO SCHEDULE(STATIC)
      DO edge = 1, edge_num
         ij_left = left(edge)
         ij_right = right(edge)
         !DIR$ SIMD
         DO l = 1, llm
            du(l,edge) = 1.*(berni(l,ij_left)-berni(l,ij_right)) ! minus gradient
         END DO
      END DO
      !$OMP END DO
   ELSE
      !$OMP DO SCHEDULE(STATIC)
      DO edge = 1, edge_num
         ij_left = left(edge)
         ij_right = right(edge)
         !DIR$ SIMD
         DO l = 1, llm
            du(l,edge) = du(l,edge) + 1.*(berni(l,ij_left)-berni(l,ij_right)) ! minus gradient
         END DO
      END DO
      !$OMP END DO
   END IF
   !---------------------------- caldyn_slow_hydro ----------------------------------
   !--------------------------------------------------------------------------