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