% This script loads a grid object already created by a previous
% gridgen session, but creates the model bathymetry and land mask
% by reading in a high-resolution lat/lon topographic grid that contains
% both land and ocean values (such as etopo5 or the new NGDC coastal
% relief model) and using Matlab's INTERP2D to determine the model
% grid values.

% 1. Load the grid object.  Here we previously ran "gridgen"
% up to the point where we were happy with the grid resolution, but instead
% of clicking "commit" and proceeding to grid masking, we typed
% "global grid_obj" in the Matlab command window and then 
% "save casco_grid_obj grid_obj" to save the  grid object with a unique name.

load casco_grid_obj
global grid_obj

% specify the projection to be what was used in grid generation
m_proj(grid_obj.projection);

% call the masking command, but only execute "special" parts of it
% (the parts that compute grid parameters, but don't actually mask)

gridgen_masking_command('special');
x=grid_obj.x_rho;
y=grid_obj.y_rho;

% read in a hi-res geographic topography grid
[xb,yb,zb]=read_gmt('sharman.grd');

% convert the geographic grid (lon/lat) to mercator
[xb,yb]=m_ll2xy(xb,yb);

% determine the topography at the model grid cell centers
zi=interp2(xb,yb,zb,x,y);

% 
zi=-zi;   %make depths positive 
land=-99999.;
%land=nan;
zi(zi<=0)=land;  % make land values -99999.
zi(zi~=land & zi<3)=3;   % 3 m minimum depth

[m,n]=size(zi);

% one dummy row around whole grid
zi(:,[1 n])=land;
zi([1 m],:)=land;
zi([2 m-1],[2 n-1])=land;   %

grid_obj.grid_bathymetry=zi;

grid_obj.ecomfile='model_grid';

% output the grid information to ECOM "model_grid" format
gridgen_output_command('special');
