/[MITgcm]/MITgcm_contrib/gael/matlab_class/ecco_v4/v4_basin.m
ViewVC logotype

Annotation of /MITgcm_contrib/gael/matlab_class/ecco_v4/v4_basin.m

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.8 - (hide annotations) (download)
Wed Jan 13 15:54:55 2016 UTC (9 years, 6 months ago) by gforget
Branch: MAIN
CVS Tags: checkpoint65x, checkpoint65v, checkpoint65w, checkpoint65t, checkpoint65u, checkpoint66f, checkpoint66e, checkpoint66d, checkpoint66c, checkpoint66b, checkpoint66a, checkpoint66o, HEAD
Changes since 1.7: +10 -8 lines
- add v4_basin.bin in ecco_v4/
- fix detection of basin_masks_eccollc_90x50.bin in old location

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

  ViewVC Help
Powered by ViewVC 1.1.22