matlab_class/ecco_v4/v4_basin.m

function [varargout]=v4_basin(nameBasin,varargin);
%object:        obtain the mask of an ocean basin
%
%input:         nameBasin       name of the basin of interest (atl, pac, ind, arctic, etc.)
%optional:      msk0            value for the masked region (0 by default)
%
%output:        mskC            mask for tracer points (msk0=outside basin; 1=inside basin)
%optional:      mskW,mskS       mask for velocity points (0=oustide/land; 1/2=edge; 1=inside)

if nargin==2; msk0=varargin{1}; else; msk0=0; end;

gcmfaces_global;

dir0=which('gcmfaces_demo');
tmp1=strfind(dir0,'/'); 
dir0=[dir0(1:tmp1(end)) 'sample_input/OCCAetcONv4GRID/'];
msk=v4_read_bin([dir0 'basin_masks_eccollc_90x50.bin'],0,1);

%list of available basins:
list0={ 'pac','atl','ind','arct','bering','southChina','mexico','okhotsk','hudson',...
        'med','java','north','japan','timor','eastChina','red','gulf'};
%list of selected basins:
if ischar(nameBasin); nameBasin={nameBasin}; end;
list1={};
for ii=1:length(nameBasin);
  if strcmp(nameBasin{ii},'atlExt');
    list1={list1{:},'atl','mexico','hudson','med','north'};
  elseif strcmp(nameBasin{ii},'pacExt');
    list1={list1{:},'pac','bering','okhotsk','japan','eastChina'};
  elseif strcmp(nameBasin{ii},'indExt');
    list1={list1{:},'ind','southChina','java','timor','red','gulf'};
  else;
    list1={list1{:},nameBasin{ii}};
  end;
end 

%derive tracer points mask:
mskC=0*msk;
for ii=1:length(list1);
  jj=find(strcmp(list1{ii},list0));
  if ~isempty(jj); mskC(find(msk==jj))=1; end;
end;

%determine velocity points masks, if needed:
if nargout>1;
  %flag velocity points according to neighboring pair: 
  fld=3*mskC+1*(~isnan(mygrid.mskC(:,:,1)));
  FLD=exch_T_N(fld);
  fldW=fld; fldS=fld;
  for iF=1:FLD.nFaces;
     tmpA=FLD{iF}(2:end-1,2:end-1);
     tmpB=FLD{iF}(1:end-2,2:end-1);
     fldW{iF}=(tmpA+tmpB)/2;
     tmpA=FLD{iF}(2:end-1,2:end-1);
     tmpB=FLD{iF}(2:end-1,1:end-2);
     fldS{iF}=(tmpA+tmpB)/2;
  end;
  %compute corresponding masks:
  mskW=0*mskC;
  mskW(find(fldW==4))=1;%inside points
  mskW(find(fldW==2.5))=0.5;%basin edge points
  mskS=0*mskC; 
  mskS(find(fldS==4))=1;%inside points
  mskS(find(fldS==2.5))=0.5;%basin edge points
  %for checking:
  if 0;
  mskWout=0*mskC;
  mskWout(find(fldW==1))=1;%outside points
  mskWout(find(fldW==2.5))=0.5;%basin edge points
  mskSout=0*mskC;
  mskSout(find(fldS==1))=1;%outside points
  mskSout(find(fldS==2.5))=0.5;%basin edge points
  end;
end;

%replace 0 with msk0:
mskC(find(mskC==0))=msk0;
if nargout>1; mskW(find(mskW==0))=msk0; mskS(find(mskS==0))=msk0; end;

%output(s):
if nargout==1; varargout={mskC}; else; varargout={mskC,mskW,mskS}; end;

%for checking:
if 0;
figure;
msk0=1*(msk0>0); msk0(find(msk0==0))=NaN;
subplot(2,1,1); imagescnan(convert2array(msk0)'); axis xy; caxis([-1 2]);
subplot(2,1,2); imagescnan(convert2array(mskC.*msk0)'); axis xy; caxis([-1 2]);  
drawnow;
end;


1	gforget	1.1	function [varargout]=v4_basin(nameBasin,varargin);
2			%object: obtain the mask of an ocean basin
3			%
4			%input: nameBasin name of the basin of interest (atl, pac, ind, arctic, etc.)
5			%optional: msk0 value for the masked region (0 by default)
6			%
7			%output: mskC mask for tracer points (msk0=outside basin; 1=inside basin)
8			%optional: mskW,mskS mask for velocity points (0=oustide/land; 1/2=edge; 1=inside)
9
10			if nargin==2; msk0=varargin{1}; else; msk0=0; end;
11
12			gcmfaces_global;
13
14			dir0=which('gcmfaces_demo');
15			tmp1=strfind(dir0,'/');
16			dir0=[dir0(1:tmp1(end)) 'sample_input/OCCAetcONv4GRID/'];
17			msk=v4_read_bin([dir0 'basin_masks_eccollc_90x50.bin'],0,1);
18
19			%list of available basins:
20			list0={ 'pac','atl','ind','arct','bering','southChina','mexico','okhotsk','hudson',...
21			'med','java','north','japan','timor','eastChina','red','gulf'};
22			%list of selected basins:
23			if ischar(nameBasin); nameBasin={nameBasin}; end;
24			list1={};
25			for ii=1:length(nameBasin);
26			if strcmp(nameBasin{ii},'atlExt');
27			list1={list1{:},'atl','mexico','hudson','med','north'};
28			elseif strcmp(nameBasin{ii},'pacExt');
29			list1={list1{:},'pac','bering','okhotsk','japan','eastChina'};
30			elseif strcmp(nameBasin{ii},'indExt');
31			list1={list1{:},'ind','southChina','java','timor','red','gulf'};
32			else;
33			list1={list1{:},nameBasin{ii}};
34			end;
35			end
36
37			%derive tracer points mask:
38			mskC=0*msk;
39			for ii=1:length(list1);
40			jj=find(strcmp(list1{ii},list0));
41			if ~isempty(jj); mskC(find(msk==jj))=1; end;
42			end;
43
44			%determine velocity points masks, if needed:
45			if nargout>1;
46			%flag velocity points according to neighboring pair:
47			fld=3mskC+1(~isnan(mygrid.mskC(:,:,1)));
48			FLD=exch_T_N(fld);
49			fldW=fld; fldS=fld;
50			for iF=1:FLD.nFaces;
51			tmpA=FLD{iF}(2:end-1,2:end-1);
52			tmpB=FLD{iF}(1:end-2,2:end-1);
53			fldW{iF}=(tmpA+tmpB)/2;
54			tmpA=FLD{iF}(2:end-1,2:end-1);
55			tmpB=FLD{iF}(2:end-1,1:end-2);
56			fldS{iF}=(tmpA+tmpB)/2;
57			end;
58			%compute corresponding masks:
59			mskW=0*mskC;
60			mskW(find(fldW==4))=1;%inside points
61			mskW(find(fldW==2.5))=0.5;%basin edge points
62			mskS=0*mskC;
63			mskS(find(fldS==4))=1;%inside points
64			mskS(find(fldS==2.5))=0.5;%basin edge points
65			%for checking:
66			if 0;
67			mskWout=0*mskC;
68			mskWout(find(fldW==1))=1;%outside points
69			mskWout(find(fldW==2.5))=0.5;%basin edge points
70			mskSout=0*mskC;
71			mskSout(find(fldS==1))=1;%outside points
72			mskSout(find(fldS==2.5))=0.5;%basin edge points
73			end;
74			end;
75
76			%replace 0 with msk0:
77			mskC(find(mskC==0))=msk0;
78			if nargout>1; mskW(find(mskW==0))=msk0; mskS(find(mskS==0))=msk0; end;
79
80			%output(s):
81			if nargout==1; varargout={mskC}; else; varargout={mskC,mskW,mskS}; end;
82
83			%for checking:
84			if 0;
85			figure;
86			msk0=1*(msk0>0); msk0(find(msk0==0))=NaN;
87			subplot(2,1,1); imagescnan(convert2array(msk0)'); axis xy; caxis([-1 2]);
88			subplot(2,1,2); imagescnan(convert2array(mskC.*msk0)'); axis xy; caxis([-1 2]);
89			drawnow;
90			end;
91
92