function patchsurf3D(gridX1,gridY1,gridZ1,RG,opt)
if nargin<5
    opt=1;
end
if gridZ1==0
    gridZ1=zeros(size(gridX1));
    dimgrid=2;
else
    dimgrid=3;
end
hold on
for r=1:size(gridX1,1)-1
    for c=1:size(gridX1,2)-1
        xcoord1=[gridX1(r,c) gridX1(r,c+1) gridX1(r+1,c+1) gridX1(r+1,c)]';
        ycoord1=[gridY1(r,c) gridY1(r,c+1) gridY1(r+1,c+1) gridY1(r+1,c)]';
        zcoord1=[gridZ1(r,c) gridZ1(r,c+1) gridZ1(r+1,c+1) gridZ1(r+1,c)]';
        %area1=polyarea(xcoord1,ycoord1);
        %decompose quadrilateral (which may be non planar) into 2 triangles
        %and add areas of the 2 triangles
        
        if ~isnan(RG(r,c)) && dimgrid==3
            if opt==2
                h=patch(xcoord1,ycoord1,zcoord1,[1 1 1]);
            else
                h=patch(xcoord1,ycoord1,zcoord1,RG(r,c));
            end
            if opt==0
                set(h,'LineStyle','none');
            end
        elseif dimgrid==2 && ~isnan(RG(r,c))
            if opt==2
                h=patch(xcoord1,ycoord1,[1 1 1]);
            else
                h=patch(xcoord1,ycoord1,RG(r,c));
            end
            if opt==0
                set(h,'LineStyle','none');
            end
        end
        
    end
end

disp('')