clear all 
close all

dirOut='jplroms/';
gridFile='jplroms/mb_das_2012012403.nc';

%--------------------
%load and compute grid
%--------------------

%----horizontal grid
wc=load('wdb_west_coast.mat');
lat=ncread(gridFile,'lat');
lon=ncread(gridFile,'lon')-360;
[Y,X]=meshgrid(lat,lon);

%----vertical grid
z=-ncread(gridFile,'depth');
u=ncread(gridFile,'u');
%pcolorjw(X,Y,squeeze(u(:,:,1,1)))

%bottom idx and bathymetry
belowBottomIdx=false(size(u(:,:,:,1)));
Z=min(z)*ones(size(X));
for i=1:size(u,1)
  for j=1:size(u,2)
%    idx_=find(isnan(u(i,j,:)),1);
    idx_=squeeze(isnan(u(i,j,:,1)));
    if sum(idx_)==0%water column extedns to all levels
      idx=size(u,3);
    else
      if sum(idx_)>1%for wet nodes
        idx=[max(1,find(idx_)-1):1:size(u,3)];
      else
        idx=size(u,3);
      end
    end
    belowBottomIdx(i,j,idx)=1;
    Z(i,j)=z(min(idx));
  end
end

zcoord=z;
save(fullfile(dirOut,'grid.mat'),'X','Y','Z','zcoord','wc','belowBottomIdx')

