!--------------------------------------------------------------------------
   !---------------------------- wind_centered ----------------------------------
   ! Perot reconstruction based on Gauss theorem
   ! u = sum( u.edge_normal * edge_length * (edge_midpoint-cell_centroid) ) /cell_area
   DO l = ll_begin, ll_end
      !DIR$ SIMD
      DO ij=ij_begin, ij_end
         cx=centroid(ij,1)
         cy=centroid(ij,2)
         cz=centroid(ij,3)
         ux=0. ; uy=0. ; uz=0.
         ue_le = ne_rup*ue(ij+u_rup,l)*le(ij+u_rup)
         ux = ux + ue_le*(.5*(xyz_v(ij+z_rup,1)+xyz_v(ij+z_up,1))-cx)
         uy = uy + ue_le*(.5*(xyz_v(ij+z_rup,2)+xyz_v(ij+z_up,2))-cy)
         uz = uz + ue_le*(.5*(xyz_v(ij+z_rup,3)+xyz_v(ij+z_up,3))-cz)
         ue_le = ne_lup*ue(ij+u_lup,l)*le(ij+u_lup)
         ux = ux + ue_le*(.5*(xyz_v(ij+z_lup,1)+xyz_v(ij+z_up,1))-cx)
         uy = uy + ue_le*(.5*(xyz_v(ij+z_lup,2)+xyz_v(ij+z_up,2))-cy)
         uz = uz + ue_le*(.5*(xyz_v(ij+z_lup,3)+xyz_v(ij+z_up,3))-cz)
         ue_le = ne_left*ue(ij+u_left,l)*le(ij+u_left)
         ux = ux + ue_le*(.5*(xyz_v(ij+z_lup,1)+xyz_v(ij+z_ldown,1))-cx)
         uy = uy + ue_le*(.5*(xyz_v(ij+z_lup,2)+xyz_v(ij+z_ldown,2))-cy)
         uz = uz + ue_le*(.5*(xyz_v(ij+z_lup,3)+xyz_v(ij+z_ldown,3))-cz)
         ue_le = ne_ldown*ue(ij+u_ldown,l)*le(ij+u_ldown)
         ux = ux + ue_le*(.5*(xyz_v(ij+z_ldown,1)+xyz_v(ij+z_down,1))-cx)
         uy = uy + ue_le*(.5*(xyz_v(ij+z_ldown,2)+xyz_v(ij+z_down,2))-cy)
         uz = uz + ue_le*(.5*(xyz_v(ij+z_ldown,3)+xyz_v(ij+z_down,3))-cz)
         ue_le = ne_rdown*ue(ij+u_rdown,l)*le(ij+u_rdown)
         ux = ux + ue_le*(.5*(xyz_v(ij+z_rdown,1)+xyz_v(ij+z_down,1))-cx)
         uy = uy + ue_le*(.5*(xyz_v(ij+z_rdown,2)+xyz_v(ij+z_down,2))-cy)
         uz = uz + ue_le*(.5*(xyz_v(ij+z_rdown,3)+xyz_v(ij+z_down,3))-cz)
         ue_le = ne_right*ue(ij+u_right,l)*le(ij+u_right)
         ux = ux + ue_le*(.5*(xyz_v(ij+z_rup,1)+xyz_v(ij+z_rdown,1))-cx)
         uy = uy + ue_le*(.5*(xyz_v(ij+z_rup,2)+xyz_v(ij+z_rdown,2))-cy)
         uz = uz + ue_le*(.5*(xyz_v(ij+z_rup,3)+xyz_v(ij+z_rdown,3))-cz)
         fac = scale*(1./Ai(ij))
         ucenter(ij,l,1)=ux*fac
         ucenter(ij,l,2)=uy*fac
         ucenter(ij,l,3)=uz*fac
      END DO
   END DO
   !---------------------------- wind_centered ----------------------------------
   !--------------------------------------------------------------------------