function W=BoundaryWOutput(U,CG,Param)
OP=CG.OrdPoly+1;
NF=Param.Grid.NumFaces;

v1CG=reshape(U(reshape(CG.Glob,OP*OP*NF,1),1,Param.uPos)...
  ,OP,OP,NF);
v2CG=reshape(U(reshape(CG.Glob,OP*OP*NF,1),1,Param.vPos)...
  ,OP,OP,NF,1);

wCG=...
  -(Param.dXdxIC(:,:,:,1,3,1).*v1CG(:,:,:,1)...
  +Param.dXdxIC(:,:,:,1,3,2).*v2CG(:,:,:,1))...
  ./Param.dXdxIC(:,:,:,1,3,3).*Param.JC(:,:,:,1);
W=zeros(size(U,1),1);
for iM=1:size(CG.FaceGlob,1)
  W(reshape(CG.Glob(:,CG.FaceGlob(iM).Ind,:)...
    ,OP*OP*size(CG.FaceGlob(iM).Ind,1),1),:)...
    =W(reshape(CG.Glob(:,CG.FaceGlob(iM).Ind,:)...
    ,OP*OP*size(CG.FaceGlob(iM).Ind,1),1),:)...
    +reshape(wCG(:,:,CG.FaceGlob(iM).Ind,:)...
    ,OP*OP*size(CG.FaceGlob(iM).Ind,1),1);
end
W=W./CG.M(:,1);
end
