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,filesep); 
dir0=[dir0(1:tmp1(end)) 'sample_input/OCCAetcONv4GRID/'];
fil0='basin_masks_eccollc_90x50.bin';

if ~isdir([dir0 fil0]);
    warning(['v4_basin requires ' fil0 ' that was not found ---> skip v4_basin!']);
    varargout={[]};
    return;
end;

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