%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright (c) 2001 Rutgers University.                                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                           %
%  Matlab script to process Land/Sea Mask Data.                             %
%                                                                           %
%  It reads the GRID NetCDF file generated by the grid package and creates  %
%  or updates the Land/Sea masking fields:                                  %
%                                                                           %
%     mask_rho    mask on RHO-points (real matrix).                         %
%     mask_u      mask on U-points (real matrix).                           %
%     mask_v      mask on V-points (real matrix).                           %
%     mask_psi    mask on PSI-points (real matrix).                         %
%                                                                           %
%  Grid points on land have a value of zero whereas grid points on water    %
%  have a value of one.                                                     %
%                                                                           %
%  This script runs interactevely.                                          %
%                                                                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global IPRINT
IPRINT=0;

%----------------------------------------------------------------------------
%  Inquire to the USER about input GRID file name.
%----------------------------------------------------------------------------

grdname=input('Enter GRID NetCDF file name:  ','s');

%----------------------------------------------------------------------------
%  Read in raw bathymetry and Land/Sea mask (if any) at RHO-points from GRID
%  NetCDF file.
%----------------------------------------------------------------------------

[spherical,rlon,rlat,bath,mask]=read_mask(grdname);

%  If MASK is flagged out, inquire the value of the transition isobath
%  and then set the mask data.

maxmask=max(max(mask));
minmask=min(min(mask));
if ((minmask == maxmask) | isempty(mask)),
  zmsk_def=0;
  zmsk=input(['Enter transition isobath (+/- meters) between land and sea',...
             ' points [',num2str(zmsk_def),'] : ']);
  if (isempty(zmsk)),
    zmsk=zmsk_def;
  end,
  mask=ones(size(bath));
  ind=find(bath < zmsk);
  mask(ind)=zeros(size(ind));
  clear ind
end
clear maxmask minmask

%----------------------------------------------------------------------------
%  Set (I,J) mask coordinates on RHO-points.
%----------------------------------------------------------------------------

[Lp,Mp]=size(rlon);
L=Lp-1;
M=Mp-1;

[y,x]=meshgrid(1:Mp,1:Lp);

%----------------------------------------------------------------------------
%  If spherical set-up, read in coastline file and convert to map projection
%  coordinates.
%----------------------------------------------------------------------------

if (spherical),
  cstname=input('Enter coastline (Icst,Jcst) data file name: ','s');
  C = load(cstname);
  xcst = C.lon;
  ycst = C.lat;
  clear C
  %load(cstname);
  %xcst=C.Icst;
  %ycst=C.Jcst;
  %clear C
else,
  xcst=[];
  ycst=[];
end,

%----------------------------------------------------------------------------
%  Plot Land/Sea mask on RHO-points and overlay Map coastlines.
%----------------------------------------------------------------------------

msktitle=['Land/Sea mask on RHO-points'];
h=pltmask(mask,x,y,0,xcst,ycst,msktitle);

%----------------------------------------------------------------------------
%  Process land/sea mask: create a new mask; mask/unmask points.
%----------------------------------------------------------------------------

figure;
set_mask

%----------------------------------------------------------------------------
%  Write out mask data into GRID NetCDF file.
%----------------------------------------------------------------------------

if (~abort),
  status=write_mask(grdname,rmask,umask,vmask,pmask);
end

%----------------------------------------------------------------------------
%  Plot Land/Sea mask on U-, V-, and PSI-points and overlay Map coastlines.
%----------------------------------------------------------------------------

if (~abort),

%  Plot Land/Sea mask on U-points.

  figure;
  msktitle=['Land/Sea mask on U-points'];
  [yugrd,xugrd]=meshgrid(1:Mp,1:L);
  pltmask(umask,xugrd,yugrd,0,xcst,ycst,msktitle);
  clear yugrd xugrd

%  Plot Land/Sea mask on V-points.

  figure;
  msktitle=['Land/Sea mask on V-points'];
  [yvgrd,xvgrd]=meshgrid(1:M,1:Lp);
  pltmask(vmask,xvgrd,yvgrd,0,xcst,ycst,msktitle);
  clear yvgrd xvgrd

%  Plot Land/Sea mask on PSI-points.

  figure;
  msktitle=['Land/Sea mask on PSI-points'];
  [ypgrd,xpgrd]=meshgrid(1:M,1:L);
  pltmask(pmask,xpgrd,ypgrd,0,xcst,ycst,msktitle);
  clear ypgrd xpgrd

end,
