%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright (c) 1996 Rutgers University.                                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                           %
%  Matlab script to process SCRUM's 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

%  Read in positions at U-, V-, and PSI-points.

[ulon,ulat,vlon,vlat,plon,plat]=read_uvpgrid(grdname);

%----------------------------------------------------------------------------
%  If spherical set-up, get geographical map projection.
%----------------------------------------------------------------------------

if (spherical == 'T' | spherical == 't'),

  disp(' ');
  disp('Available Map Projections: ');
  disp(' ');
  disp('         [ 1] Polar Stereographic');
  disp('         [ 2] Orthographic');
  disp('         [ 3] Azimuthal Equal-Area');
  disp('         [ 4] Azimuthal Equidistant');
  disp('         [ 5] Gnomonic');
  disp('         [ 6] Satellite');
  disp('         [ 7] Albers Equal-Area Conic');
  disp('         [ 8] Lambert Conformal Conic');
  disp('         [ 9] Mercator');
  disp('         [10] Miller Cylindrical');
  disp('         [11] Equidistant Cylindrical');
  disp('         [12] Oblique Mercator');
  disp('         [13] Transverse Mercator');
  disp('         [14] Sinusoidal');
  disp(' ');

  iproj=input('Enter Map Projection Index:  ');

% Set and initialize Map Projection.

  switch iproj,
    case 1,
      m_proj('Stereographic');
    case 2,
      m_proj('Orthographic');
    case 3,
      m_proj('Azimuthal Equal-Area');
    case 4,
      m_proj('Azimuthal Equidistant');
    case 5,
      m_proj('Gnomonic');
    case 6,
      m_proj('Satellite');
    case 7,
      m_proj('Albers Equal-Area Conic');
    case 8,
      m_proj('Lambert Conformal Conic');
    case 9,
      m_proj('Mercator');
    case 10,
      m_proj('Miller Cylindrical');
    case 11,
      m_proj('Equidistant Cylindrical');
    case 12,
      m_proj('Oblique Mercator');
    case 13,
      m_proj('Transverse Mercator');
    case 14,
      m_proj('Sinusoidal');
    otherwise
      disp([' Unknown Map Projection Option: ',num2str(iproj)]);
  end,
end,

%----------------------------------------------------------------------------
%  If spherical set-up, convert grid to map projection coordinates.
%----------------------------------------------------------------------------

if (spherical == 'T' | spherical == 't'),
  [xrgrd,yrgrd]=m_ll2xy(rlon,rlat);
  [xugrd,yugrd]=m_ll2xy(ulon,ulat);
  [xvgrd,yvgrd]=m_ll2xy(vlon,vlat);
  [xpgrd,ypgrd]=m_ll2xy(plon,plat);
else,
  xrgrd=rlon; yrgrd=rlat;
  xugrd=ulon; yugrd=ulat;
  xvgrd=vlon; yvgrd=vlat;
  xpgrd=plon; ypgrd=plat;
end,

%clear rlon rlat ulon ulat vlon vlat plon plat

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

if (spherical == 'T' | spherical == 't'),
  cstname=input('Enter coastline (lat,lon) data file name: ','s');
  [cstlon,cstlat]=rcoastline(cstname);
  spval=900.0;
  [icst]=find(cstlon>=spval);
  [xcst,ycst]=m_ll2xy(cstlon,cstlat);
  xcst(icst)=NaN;
  ycst(icst)=NaN;
  clear spval cstlon cstlat
else,
  xcst=[];
  ycst=[];
end,

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

msktitle=['Land/Sea mask on RHO-points'];
pltmask(mask,xrgrd,yrgrd,spherical,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'];
  pltmask(umask,xugrd,yugrd,spherical,xcst,ycst,msktitle);

%  Plot Land/Sea mask on V-points.

  figure;
  msktitle=['Land/Sea mask on V-points'];
  pltmask(vmask,xvgrd,yvgrd,spherical,xcst,ycst,msktitle);

%  Plot Land/Sea mask on PSI-points.

  figure;
  msktitle=['Land/Sea mask on PSI-points'];
  pltmask(pmask,xpgrd,ypgrd,spherical,xcst,ycst,msktitle);

end,
