function plot_gridded(md,data,options,plotlines,plotcols,i)
%PLOT_OVERLAY - superimpose radar image to a given field
%
%   Usage:
%      plot_gridded(md,options,plotlines,plotcols,i)
%
%   See also: PLOTMODEL

whitepos=getfieldvalue(options,'whitepos',2); %1: up, 2: down, else: none

%process mesh and data
[x y z elements is2d isplanet]=processmesh(md,[],options);
[data datatype]=processdata(md,data,options);

%check is2d
if ~is2d, 
	error('buildgridded error message: gridded not supported for 3d meshes, project on a layer');
end

%Get xlim and ylim (used to extract radar image)
xlim=getfieldvalue(options,'xlim',[min(x) max(x)]);
ylim=getfieldvalue(options,'ylim',[min(y) max(y)]);
post=getfieldvalue(options,'posting',diff(xlim)/1000);

%Interpolating data on grid
[x_m y_m data_grid]=InterpFromMeshToGrid(elements,x,y,data,xlim(1),ylim(2),post,post,round(diff(ylim)/post),round(diff(xlim)/post),NaN);
if size(data_grid,1)<3 | size(data_grid,2)<3,
	error('data_grid size too small in plot_gridded, check posting and uni');
end

%Process data_grid: add white in NaN and correct caxis accordingly
if exist(options,'caxis'),
	caxis_opt=getfieldvalue(options,'caxis');
	data_grid(find(data_grid<caxis_opt(1)))=caxis_opt(1);
	data_grid(find(data_grid>caxis_opt(2)))=caxis_opt(2);
	data_min=caxis_opt(1);
	data_max=caxis_opt(2);
else
	data_min=min(data_grid(:));
	data_max=max(data_grid(:));
end
options=changefieldvalue(options,'cbYLim',[data_min data_max]);
if whitepos==1,
	white  =data_max + (data_max-data_min)/55;
	options=changefieldvalue(options,'caxis',[data_min white]);
	data_grid(isnan(data_grid))=white;
elseif whitepos==2,
	white  =data_min - (data_max-data_min)/55;
	options=changefieldvalue(options,'caxis',[white data_max]);
	data_grid(isnan(data_grid))=white;
end

%Select plot area 
subplot(plotlines,plotcols,i);

%shading interp;
h=imagesc(xlim,ylim,data_grid);set(gca,'YDir','normal');
map=getfieldvalue(options,'colormap',jet);
if whitepos==1,
	map(end,:)=[1 1 1];
elseif whitepos==2,
	map(1,:)=[1 1 1];
end
options=changefieldvalue(options,'colormap',map);

%last step: mesh gridded?
if exist(options,'edgecolor'),
	A=elements(:,1); B=elements(:,2); C=elements(:,3); 
	patch('Faces',[A B C],'Vertices', [x y z],'FaceVertexCData',data_grid(1)*ones(size(x)),'FaceColor','none','EdgeColor',getfieldvalue(options,'edgecolor'));
end

%Apply options
applyoptions(md,data,options);
