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