!--------------------------------------------------------------------------
   !---------------------------- caldyn_slow_hydro ----------------------------------
   !$OMP DO SCHEDULE(STATIC)
   DO edge = 1, edge_num
      ij_left = left(edge)
      ij_right = right(edge)
      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
      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 DO
   !$OMP END DO
   IF(zero) THEN
      !$OMP DO SCHEDULE(STATIC)
      DO edge = 1, edge_num
         ij_left = left(edge)
         ij_right = right(edge)
         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)
         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 ----------------------------------
   !--------------------------------------------------------------------------