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);
%ph:
%dir0='/raid3/gforget/ecco_v4/input_files/';
%dir0='/net/nares/raid11/ecco-shared/ecco-version-4/input/input_rest/';
%msk=v4_read_bin([dir0 'basin_masks_ext_ecco_v4.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',...
        'baffin','gin','barents'};
%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','baffin','gin'};
  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	gforget	1.3	msk=v4_read_bin([dir0 'basin_masks_eccollc_90x50.bin'],0,1);
18			%ph:
19			%dir0='/raid3/gforget/ecco_v4/input_files/';
20			%dir0='/net/nares/raid11/ecco-shared/ecco-version-4/input/input_rest/';
21			%msk=v4_read_bin([dir0 'basin_masks_ext_ecco_v4.bin'],0,1);
22	gforget	1.1
23			%list of available basins:
24	heimbach	1.2	list0={ 'pac','atl','ind','arct','bering',...
25			'southChina','mexico','okhotsk','hudson','med',...
26			'java','north','japan','timor','eastChina','red','gulf',...
27	gforget	1.3	'baffin','gin','barents'};
28	gforget	1.1	%list of selected basins:
29			if ischar(nameBasin); nameBasin={nameBasin}; end;
30			list1={};
31			for ii=1:length(nameBasin);
32			if strcmp(nameBasin{ii},'atlExt');
33	gforget	1.3	list1={list1{:},'atl','mexico','hudson','med','north','baffin','gin'};
34	gforget	1.1	elseif strcmp(nameBasin{ii},'pacExt');
35	gforget	1.3	list1={list1{:},'pac','bering','okhotsk','japan','eastChina'};
36	gforget	1.1	elseif strcmp(nameBasin{ii},'indExt');
37			list1={list1{:},'ind','southChina','java','timor','red','gulf'};
38			else;
39			list1={list1{:},nameBasin{ii}};
40			end;
41			end
42
43			%derive tracer points mask:
44			mskC=0*msk;
45			for ii=1:length(list1);
46			jj=find(strcmp(list1{ii},list0));
47			if ~isempty(jj); mskC(find(msk==jj))=1; end;
48			end;
49
50			%determine velocity points masks, if needed:
51			if nargout>1;
52			%flag velocity points according to neighboring pair:
53			fld=3mskC+1(~isnan(mygrid.mskC(:,:,1)));
54			FLD=exch_T_N(fld);
55			fldW=fld; fldS=fld;
56			for iF=1:FLD.nFaces;
57			tmpA=FLD{iF}(2:end-1,2:end-1);
58			tmpB=FLD{iF}(1:end-2,2:end-1);
59			fldW{iF}=(tmpA+tmpB)/2;
60			tmpA=FLD{iF}(2:end-1,2:end-1);
61			tmpB=FLD{iF}(2:end-1,1:end-2);
62			fldS{iF}=(tmpA+tmpB)/2;
63			end;
64			%compute corresponding masks:
65			mskW=0*mskC;
66			mskW(find(fldW==4))=1;%inside points
67			mskW(find(fldW==2.5))=0.5;%basin edge points
68			mskS=0*mskC;
69			mskS(find(fldS==4))=1;%inside points
70			mskS(find(fldS==2.5))=0.5;%basin edge points
71			%for checking:
72			if 0;
73			mskWout=0*mskC;
74			mskWout(find(fldW==1))=1;%outside points
75			mskWout(find(fldW==2.5))=0.5;%basin edge points
76			mskSout=0*mskC;
77			mskSout(find(fldS==1))=1;%outside points
78			mskSout(find(fldS==2.5))=0.5;%basin edge points
79			end;
80			end;
81
82			%replace 0 with msk0:
83			mskC(find(mskC==0))=msk0;
84			if nargout>1; mskW(find(mskW==0))=msk0; mskS(find(mskS==0))=msk0; end;
85
86			%output(s):
87			if nargout==1; varargout={mskC}; else; varargout={mskC,mskW,mskS}; end;
88
89			%for checking:
90			if 0;
91			figure;
92			msk0=1*(msk0>0); msk0(find(msk0==0))=NaN;
93			subplot(2,1,1); imagescnan(convert2array(msk0)'); axis xy; caxis([-1 2]);
94			subplot(2,1,2); imagescnan(convert2array(mskC.*msk0)'); axis xy; caxis([-1 2]);
95			drawnow;
96			end;
97
98