!--------------------------------------------------------------------------
   !---------------------------- wind_centered ----------------------------------
   ! Perot reconstruction based on Gauss theorem
   ! u = sum( u.edge_normal * edge_length * (edge_midpoint-cell_centroid) ) /cell_area
   !$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)
         sign1 = primal_ne(1,ij)
         sign2 = primal_ne(2,ij)
         sign3 = primal_ne(3,ij)
         sign4 = primal_ne(4,ij)
         ij_up1 = up(edge1)
         ij_up2 = up(edge2)
         ij_up3 = up(edge3)
         ij_up4 = up(edge4)
         ij_down1 = down(edge1)
         ij_down2 = down(edge2)
         ij_down3 = down(edge3)
         ij_down4 = down(edge4)
         !DIR$ SIMD
         DO l = 1, llm
            cx=centroid(ij,1)
            cy=centroid(ij,2)
            cz=centroid(ij,3)
            ux=0. ; uy=0. ; uz=0.
            ue_le = sign1*ue(l,edge1)*le(edge1)
            ux = ux + ue_le*(.5*(xyz_v(ij_up1,1)+xyz_v(ij_down1,1))-cx)
            uy = uy + ue_le*(.5*(xyz_v(ij_up1,2)+xyz_v(ij_down1,2))-cy)
            uz = uz + ue_le*(.5*(xyz_v(ij_up1,3)+xyz_v(ij_down1,3))-cz)
            ue_le = sign2*ue(l,edge2)*le(edge2)
            ux = ux + ue_le*(.5*(xyz_v(ij_up2,1)+xyz_v(ij_down2,1))-cx)
            uy = uy + ue_le*(.5*(xyz_v(ij_up2,2)+xyz_v(ij_down2,2))-cy)
            uz = uz + ue_le*(.5*(xyz_v(ij_up2,3)+xyz_v(ij_down2,3))-cz)
            ue_le = sign3*ue(l,edge3)*le(edge3)
            ux = ux + ue_le*(.5*(xyz_v(ij_up3,1)+xyz_v(ij_down3,1))-cx)
            uy = uy + ue_le*(.5*(xyz_v(ij_up3,2)+xyz_v(ij_down3,2))-cy)
            uz = uz + ue_le*(.5*(xyz_v(ij_up3,3)+xyz_v(ij_down3,3))-cz)
            ue_le = sign4*ue(l,edge4)*le(edge4)
            ux = ux + ue_le*(.5*(xyz_v(ij_up4,1)+xyz_v(ij_down4,1))-cx)
            uy = uy + ue_le*(.5*(xyz_v(ij_up4,2)+xyz_v(ij_down4,2))-cy)
            uz = uz + ue_le*(.5*(xyz_v(ij_up4,3)+xyz_v(ij_down4,3))-cz)
            fac = scale*(1./Ai(ij))
            ucenter(l,ij,1)=ux*fac
            ucenter(l,ij,2)=uy*fac
            ucenter(l,ij,3)=uz*fac
         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)
         sign1 = primal_ne(1,ij)
         sign2 = primal_ne(2,ij)
         sign3 = primal_ne(3,ij)
         sign4 = primal_ne(4,ij)
         sign5 = primal_ne(5,ij)
         ij_up1 = up(edge1)
         ij_up2 = up(edge2)
         ij_up3 = up(edge3)
         ij_up4 = up(edge4)
         ij_up5 = up(edge5)
         ij_down1 = down(edge1)
         ij_down2 = down(edge2)
         ij_down3 = down(edge3)
         ij_down4 = down(edge4)
         ij_down5 = down(edge5)
         !DIR$ SIMD
         DO l = 1, llm
            cx=centroid(ij,1)
            cy=centroid(ij,2)
            cz=centroid(ij,3)
            ux=0. ; uy=0. ; uz=0.
            ue_le = sign1*ue(l,edge1)*le(edge1)
            ux = ux + ue_le*(.5*(xyz_v(ij_up1,1)+xyz_v(ij_down1,1))-cx)
            uy = uy + ue_le*(.5*(xyz_v(ij_up1,2)+xyz_v(ij_down1,2))-cy)
            uz = uz + ue_le*(.5*(xyz_v(ij_up1,3)+xyz_v(ij_down1,3))-cz)
            ue_le = sign2*ue(l,edge2)*le(edge2)
            ux = ux + ue_le*(.5*(xyz_v(ij_up2,1)+xyz_v(ij_down2,1))-cx)
            uy = uy + ue_le*(.5*(xyz_v(ij_up2,2)+xyz_v(ij_down2,2))-cy)
            uz = uz + ue_le*(.5*(xyz_v(ij_up2,3)+xyz_v(ij_down2,3))-cz)
            ue_le = sign3*ue(l,edge3)*le(edge3)
            ux = ux + ue_le*(.5*(xyz_v(ij_up3,1)+xyz_v(ij_down3,1))-cx)
            uy = uy + ue_le*(.5*(xyz_v(ij_up3,2)+xyz_v(ij_down3,2))-cy)
            uz = uz + ue_le*(.5*(xyz_v(ij_up3,3)+xyz_v(ij_down3,3))-cz)
            ue_le = sign4*ue(l,edge4)*le(edge4)
            ux = ux + ue_le*(.5*(xyz_v(ij_up4,1)+xyz_v(ij_down4,1))-cx)
            uy = uy + ue_le*(.5*(xyz_v(ij_up4,2)+xyz_v(ij_down4,2))-cy)
            uz = uz + ue_le*(.5*(xyz_v(ij_up4,3)+xyz_v(ij_down4,3))-cz)
            ue_le = sign5*ue(l,edge5)*le(edge5)
            ux = ux + ue_le*(.5*(xyz_v(ij_up5,1)+xyz_v(ij_down5,1))-cx)
            uy = uy + ue_le*(.5*(xyz_v(ij_up5,2)+xyz_v(ij_down5,2))-cy)
            uz = uz + ue_le*(.5*(xyz_v(ij_up5,3)+xyz_v(ij_down5,3))-cz)
            fac = scale*(1./Ai(ij))
            ucenter(l,ij,1)=ux*fac
            ucenter(l,ij,2)=uy*fac
            ucenter(l,ij,3)=uz*fac
         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)
         sign1 = primal_ne(1,ij)
         sign2 = primal_ne(2,ij)
         sign3 = primal_ne(3,ij)
         sign4 = primal_ne(4,ij)
         sign5 = primal_ne(5,ij)
         sign6 = primal_ne(6,ij)
         ij_up1 = up(edge1)
         ij_up2 = up(edge2)
         ij_up3 = up(edge3)
         ij_up4 = up(edge4)
         ij_up5 = up(edge5)
         ij_up6 = up(edge6)
         ij_down1 = down(edge1)
         ij_down2 = down(edge2)
         ij_down3 = down(edge3)
         ij_down4 = down(edge4)
         ij_down5 = down(edge5)
         ij_down6 = down(edge6)
         !DIR$ SIMD
         DO l = 1, llm
            cx=centroid(ij,1)
            cy=centroid(ij,2)
            cz=centroid(ij,3)
            ux=0. ; uy=0. ; uz=0.
            ue_le = sign1*ue(l,edge1)*le(edge1)
            ux = ux + ue_le*(.5*(xyz_v(ij_up1,1)+xyz_v(ij_down1,1))-cx)
            uy = uy + ue_le*(.5*(xyz_v(ij_up1,2)+xyz_v(ij_down1,2))-cy)
            uz = uz + ue_le*(.5*(xyz_v(ij_up1,3)+xyz_v(ij_down1,3))-cz)
            ue_le = sign2*ue(l,edge2)*le(edge2)
            ux = ux + ue_le*(.5*(xyz_v(ij_up2,1)+xyz_v(ij_down2,1))-cx)
            uy = uy + ue_le*(.5*(xyz_v(ij_up2,2)+xyz_v(ij_down2,2))-cy)
            uz = uz + ue_le*(.5*(xyz_v(ij_up2,3)+xyz_v(ij_down2,3))-cz)
            ue_le = sign3*ue(l,edge3)*le(edge3)
            ux = ux + ue_le*(.5*(xyz_v(ij_up3,1)+xyz_v(ij_down3,1))-cx)
            uy = uy + ue_le*(.5*(xyz_v(ij_up3,2)+xyz_v(ij_down3,2))-cy)
            uz = uz + ue_le*(.5*(xyz_v(ij_up3,3)+xyz_v(ij_down3,3))-cz)
            ue_le = sign4*ue(l,edge4)*le(edge4)
            ux = ux + ue_le*(.5*(xyz_v(ij_up4,1)+xyz_v(ij_down4,1))-cx)
            uy = uy + ue_le*(.5*(xyz_v(ij_up4,2)+xyz_v(ij_down4,2))-cy)
            uz = uz + ue_le*(.5*(xyz_v(ij_up4,3)+xyz_v(ij_down4,3))-cz)
            ue_le = sign5*ue(l,edge5)*le(edge5)
            ux = ux + ue_le*(.5*(xyz_v(ij_up5,1)+xyz_v(ij_down5,1))-cx)
            uy = uy + ue_le*(.5*(xyz_v(ij_up5,2)+xyz_v(ij_down5,2))-cy)
            uz = uz + ue_le*(.5*(xyz_v(ij_up5,3)+xyz_v(ij_down5,3))-cz)
            ue_le = sign6*ue(l,edge6)*le(edge6)
            ux = ux + ue_le*(.5*(xyz_v(ij_up6,1)+xyz_v(ij_down6,1))-cx)
            uy = uy + ue_le*(.5*(xyz_v(ij_up6,2)+xyz_v(ij_down6,2))-cy)
            uz = uz + ue_le*(.5*(xyz_v(ij_up6,3)+xyz_v(ij_down6,3))-cz)
            fac = scale*(1./Ai(ij))
            ucenter(l,ij,1)=ux*fac
            ucenter(l,ij,2)=uy*fac
            ucenter(l,ij,3)=uz*fac
         END DO
      CASE DEFAULT
         !DIR$ SIMD
         DO l = 1, llm
            cx=centroid(ij,1)
            cy=centroid(ij,2)
            cz=centroid(ij,3)
            ux=0. ; uy=0. ; uz=0.
            DO iedge = 1, primal_deg(ij)
               edge = primal_edge(iedge,ij)
               ij_up = up(edge)
               ij_down = down(edge)
               ue_le = primal_ne(iedge,ij)*ue(l,edge)*le(edge)
               ux = ux + ue_le*(.5*(xyz_v(ij_up,1)+xyz_v(ij_down,1))-cx)
               uy = uy + ue_le*(.5*(xyz_v(ij_up,2)+xyz_v(ij_down,2))-cy)
               uz = uz + ue_le*(.5*(xyz_v(ij_up,3)+xyz_v(ij_down,3))-cz)
            END DO
            fac = scale*(1./Ai(ij))
            ucenter(l,ij,1)=ux*fac
            ucenter(l,ij,2)=uy*fac
            ucenter(l,ij,3)=uz*fac
         END DO
      END SELECT
   END DO
   !$OMP END DO
   !---------------------------- wind_centered ----------------------------------
   !--------------------------------------------------------------------------