function [D]=diagdiff(FnameS,FnameP,Vname,Trec,Level);

VS=nc_read(FnameS,Vname,Trec);
VP=nc_read(FnameP,Vname,Trec);

if (nargin < 5),
  D=VS-VP;
  Smin=min(min(VS));
  Smax=max(max(VS));
else,
  D=squeeze(VS(:,:,Level)-VP(:,:,Level));
  Smin=min(min(squeeze(VS(:,:,Level))));
  Smax=max(max(squeeze(VS(:,:,Level))));
end,
Dmin=min(min(D));
Dmax=max(max(D));
[Imax,Jmax]=size(D);

[dnames,dsizes,igrid]=nc_vinfo(FnameS,Vname);
switch (igrid),
  case 1
    [x,y]=ndgrid(0:Imax-1,0:Jmax-1);
  case 2
    [x,y]=ndgrid(1:Imax  ,1:Jmax  );
  case 3
    [x,y]=ndgrid(1:Imax  ,0:Jmax-1);
  case 4
    [x,y]=ndgrid(0:Imax-1,1:Jmax  );
  case 5
    [x,y]=ndgrid(0:Imax-1,0:Jmax-1);
end,

ind=find(D == 0);
if ~isempty(D),
  D(ind)=NaN;
end

pcolorjw(x,y,D); shading faceted; colorbar;
title(untexlabel([Vname,', Difference: Min = ', num2str(Dmin), ...
                                     ' Max = ', num2str(Dmax)]));
xlabel(untexlabel(['Serial Run: Min = ', num2str(Smin), ' Max = ', num2str(Smax)]));

return

