!--------------------------------------------------------------------------
   !---------------------------- pvort_only ----------------------------------
   !$OMP DO SCHEDULE(STATIC)
   DO ij = 1, dual_num
      ! this VLOOP iterates over dual cell edges
      SELECT CASE(dual_deg(ij))
      CASE(3)
         edge1 = dual_edge(1,ij)
         edge2 = dual_edge(2,ij)
         edge3 = dual_edge(3,ij)
         sign1 = dual_ne(1,ij)
         sign2 = dual_ne(2,ij)
         sign3 = dual_ne(3,ij)
         vertex1 = dual_vertex(1,ij)
         vertex2 = dual_vertex(2,ij)
         vertex3 = dual_vertex(3,ij)
         !DIR$ SIMD
         DO l = 1, llm
            etav = 0.d0
            etav = etav + sign1*u(l,edge1)
            etav = etav + sign2*u(l,edge2)
            etav = etav + sign3*u(l,edge3)
            hv=0.
            hv = hv + Riv2(1,ij)*rhodz(l,vertex1)
            hv = hv + Riv2(2,ij)*rhodz(l,vertex2)
            hv = hv + Riv2(3,ij)*rhodz(l,vertex3)
            qv(l,ij) = (etav + fv(ij)*Av(ij) )/(hv*Av(ij))
         END DO
      CASE(4)
         edge1 = dual_edge(1,ij)
         edge2 = dual_edge(2,ij)
         edge3 = dual_edge(3,ij)
         edge4 = dual_edge(4,ij)
         sign1 = dual_ne(1,ij)
         sign2 = dual_ne(2,ij)
         sign3 = dual_ne(3,ij)
         sign4 = dual_ne(4,ij)
         vertex1 = dual_vertex(1,ij)
         vertex2 = dual_vertex(2,ij)
         vertex3 = dual_vertex(3,ij)
         vertex4 = dual_vertex(4,ij)
         !DIR$ SIMD
         DO l = 1, llm
            etav = 0.d0
            etav = etav + sign1*u(l,edge1)
            etav = etav + sign2*u(l,edge2)
            etav = etav + sign3*u(l,edge3)
            etav = etav + sign4*u(l,edge4)
            hv=0.
            hv = hv + Riv2(1,ij)*rhodz(l,vertex1)
            hv = hv + Riv2(2,ij)*rhodz(l,vertex2)
            hv = hv + Riv2(3,ij)*rhodz(l,vertex3)
            hv = hv + Riv2(4,ij)*rhodz(l,vertex4)
            qv(l,ij) = (etav + fv(ij)*Av(ij) )/(hv*Av(ij))
         END DO
      CASE DEFAULT
         !DIR$ SIMD
         DO l = 1, llm
            etav = 0.d0
            DO iedge = 1, dual_deg(ij)
               edge = dual_edge(iedge,ij)
               etav = etav + dual_ne(iedge,ij)*u(l,edge)
            END DO
            hv=0.
            DO ivertex = 1, dual_deg(ij)
               vertex = dual_vertex(ivertex,ij)
               hv = hv + Riv2(ivertex,ij)*rhodz(l,vertex)
            END DO
            qv(l,ij) = (etav + fv(ij)*Av(ij) )/(hv*Av(ij))
         END DO
      END SELECT
   END DO
   !$OMP END DO
   !$OMP DO SCHEDULE(STATIC)
   DO edge = 1, edge_num
      ij_up = up(edge)
      ij_down = down(edge)
      !DIR$ SIMD
      DO l = 1, llm
         qu(l,edge)=0.5d0*(qv(l,ij_down)+qv(l,ij_up))
      END DO
   END DO
   !$OMP END DO
   !---------------------------- pvort_only ----------------------------------
   !--------------------------------------------------------------------------