matlab_class/gcmfaces_calc/gcmfaces_loc_tile.m

function [varargout]=gcmfaces_loc_tile(ni,nj,varargin);
%object : compute map of tile indices and (optional) associate
%         select location (1D vectors) with tile indices etc.
%input :  ni,nj is the MITgcm tile size (2 numbers total)
%         (optional)
%  AND    XC,YC (vectors)
%  OR     iTile, jTile, tileNo (vectors)
%  OR     iTile, jTile, XC11, YC11, XCNINJ, YCNINJ (vectors)
%where :
%         XC, YC are lon/lat of the grid point
%         tileNo is the grid point s tile number (assuming no blank tiles)
%         iTile, jTile is the grid point tile index
%         XC11, YC11 is lat/lon for the tile SE corner
%         XCNINJ, YCNINJ is lat/lon for the tile NW corner
%
%output : (nargin==2) tileNo is the map of tile numbers
%         (nargin >3) loc_tile incl. XC,YC, tileNo, XC11, YC11, XCNINJ, YCNINJ
%
%notes : - if XC, YC are provided as input, then loc_tile.XC, YC is nearest neighbor
%        - by assumption tile numbers increase from w to e (2nd dim), 
%           then from s to n (1st dim), then with face number 

gcmfaces_global;
global mytiles;
loc_arrays={'XC','YC','XC11','YC11','XCNINJ','YCNINJ','iTile','jTile','tileNo'};

%1) check that tile arrays are up to date
if isempty(mytiles); 
  mytiles.dirGrid=mygrid.dirGrid; mytiles.nFaces=mygrid.nFaces;
  mytiles.fileFormat=mygrid.fileFormat; mytiles.ioSize=mygrid.ioSize;
  mytiles.ni=-1; mytiles.nj=-1;
end;

test1=strcmp(mytiles.dirGrid,mygrid.dirGrid)&(mytiles.nFaces==mygrid.nFaces)&...
      strcmp(mytiles.fileFormat,mygrid.fileFormat)&(sum(mytiles.ioSize~=mygrid.ioSize)==0);
test2=(mytiles.ni==ni)&(mytiles.nj==nj);

%2) update tile arrays if not up to date
if ~test1|~test2;
  %
  mytiles.dirGrid=mygrid.dirGrid; mytiles.nFaces=mygrid.nFaces;
  mytiles.fileFormat=mygrid.fileFormat; mytiles.ioSize=mygrid.ioSize;
  %
  XC=mygrid.XC; YC=mygrid.YC;
  XC11=XC; YC11=XC; XCNINJ=XC; YCNINJ=XC; iTile=XC; jTile=XC; tileNo=XC;
  tileCount=0;
  for iF=1:XC11.nFaces;
    face_XC=XC{iF}; face_YC=YC{iF};
    for ii=1:size(face_XC,1)/ni;
        for jj=1:size(face_XC,2)/nj;
            tileCount=tileCount+1;
            tmp_i=[1:ni]+ni*(ii-1);
            tmp_j=[1:nj]+nj*(jj-1);
            tmp_XC=face_XC(tmp_i,tmp_j);
            tmp_YC=face_YC(tmp_i,tmp_j);
            XC11{iF}(tmp_i,tmp_j)=tmp_XC(1,1);
            YC11{iF}(tmp_i,tmp_j)=tmp_YC(1,1);
            XCNINJ{iF}(tmp_i,tmp_j)=tmp_XC(end,end);
            YCNINJ{iF}(tmp_i,tmp_j)=tmp_YC(end,end);
            iTile{iF}(tmp_i,tmp_j)=[1:ni]'*ones(1,nj);
            jTile{iF}(tmp_i,tmp_j)=ones(ni,1)*[1:nj];
            tileNo{iF}(tmp_i,tmp_j)=tileCount*ones(ni,nj);
        end;
    end;
  end;

  for iF=1:length(loc_arrays);
    eval(['mytiles.' loc_arrays{iF} '=convert2array(' loc_arrays{iF} ');']);
    eval(['clear ' loc_arrays{iF} ';']);
  end;
end;

%3) decide what is needed based on the number or arguments
%   (this needs to happen after mytiles is computed, since 
%    I use the same variable names in that bloc and below) 
if nargin==2;
  %simply return mytiles as loc_tile
  tileNo=convert2array(mytiles.tileNo);
  varargout={tileNo};
  return;
elseif nargin==4;
  XC=varargin{1}; YC=varargin{2}; 
  tmp1=find(XC>180); XC(tmp1)=XC(tmp1)-360;
  iTile=[]; jTile=[]; tileNo=[];
  XC11=[]; YC11=[]; XCNINJ=[]; YCNINJ=[];
elseif nargin==5;
  XC=[]; YC=[];
  iTile=varargin{1}; jTile=varargin{2}; tileNo=varargin{3};
  XC11=[]; YC11=[]; XCNINJ=[]; YCNINJ=[];
elseif nargin==8;
  XC=[]; YC=[];
  iTile=varargin{1}; jTile=varargin{2}; tileNo=[];
  XC11=varargin{3}; YC11=varargin{4}; XCNINJ=varargin{5}; YCNINJ=varargin{6};
else;
  error('wrong argument list');
end;

%3) find the grid points indices in array format
if ~isempty(XC);
  %identify nearest neighbor on the sphere
  XCgrid=convert2array(mygrid.XC);
  YCgrid=convert2array(mygrid.YC);
  x=sin(pi/2-YCgrid*pi/180).*cos(XCgrid*pi/180);
  y=sin(pi/2-YCgrid*pi/180).*sin(XCgrid*pi/180);
  z=cos(pi/2-YCgrid*pi/180);
  xx=sin(pi/2-YC*pi/180).*cos(XC*pi/180);
  yy=sin(pi/2-YC*pi/180).*sin(XC*pi/180);
  zz=cos(pi/2-YC*pi/180);
  kk=find(~isnan(x));
  ik = knnsearch([x(kk) y(kk) z(kk)],[xx yy zz]);
  ik=kk(ik);
  %
  %(old method that uses longitude, latitude directly)
  %tmp1=find(XC>180); XC(tmp1)=XC(tmp1)-360;
  %ik=gcmfaces_bindata(XC,YC);
else;
  %use mytiles as a look up table
  ik=zeros(size(iTile));
  for ikk=1:length(iTile);
    if ~isempty(tileNo);
      tmp1=find(mytiles.iTile==iTile(ikk)&mytiles.jTile==jTile(ikk)&mytiles.tileNo==tileNo(ikk));
    else;
      tmp1=find(mytiles.iTile==iTile(ikk)&mytiles.jTile==jTile(ikk)&...
                abs(mytiles.XC-XC(ikk))<1e-4&...
                abs(mytiles.YC-YC(ikk))<1e-4&...
                abs(mytiles.XC11-XC11(ikk))<1e-4&...  
                abs(mytiles.YC11-YC11(ikk))<1e-4&...    
                abs(mytiles.XCNINJ-XCNINJ(ikk))<1e-4&...
                abs(mytiles.YCNINJ-YCNINJ(ikk))<1e-4);
    end;
    if isempty(tmp1); error('grid point could not be identified from inputs'); end;
    ik(ikk)=tmp1(1);
  end;
end;

%4) complement location arrays and prepare output

for iF=1:length(loc_arrays);
  eval(['loc_tile.' loc_arrays{iF} '=mytiles.' loc_arrays{iF} '(ik);']);
end;

%5) output result
varargout={loc_tile};


1	gforget	1.2	function [varargout]=gcmfaces_loc_tile(ni,nj,varargin);
2			%object : compute map of tile indices and (optional) associate
3			% select location (1D vectors) with tile indices etc.
4	gforget	1.1	%input : ni,nj is the MITgcm tile size (2 numbers total)
5	gforget	1.2	% (optional)
6	gforget	1.1	% AND XC,YC (vectors)
7			% OR iTile, jTile, tileNo (vectors)
8			% OR iTile, jTile, XC11, YC11, XCNINJ, YCNINJ (vectors)
9			%where :
10			% XC, YC are lon/lat of the grid point
11			% tileNo is the grid point s tile number (assuming no blank tiles)
12			% iTile, jTile is the grid point tile index
13			% XC11, YC11 is lat/lon for the tile SE corner
14			% XCNINJ, YCNINJ is lat/lon for the tile NW corner
15			%
16	gforget	1.2	%output : (nargin==2) tileNo is the map of tile numbers
17			% (nargin >3) loc_tile incl. XC,YC, tileNo, XC11, YC11, XCNINJ, YCNINJ
18	gforget	1.1	%
19			%notes : - if XC, YC are provided as input, then loc_tile.XC, YC is nearest neighbor
20	gforget	1.2	% - by assumption tile numbers increase from w to e (2nd dim),
21			% then from s to n (1st dim), then with face number
22	gforget	1.1
23			gcmfaces_global;
24			global mytiles;
25			loc_arrays={'XC','YC','XC11','YC11','XCNINJ','YCNINJ','iTile','jTile','tileNo'};
26
27			%1) check that tile arrays are up to date
28			if isempty(mytiles);
29			mytiles.dirGrid=mygrid.dirGrid; mytiles.nFaces=mygrid.nFaces;
30			mytiles.fileFormat=mygrid.fileFormat; mytiles.ioSize=mygrid.ioSize;
31	gforget	1.2	mytiles.ni=-1; mytiles.nj=-1;
32	gforget	1.1	end;
33
34			test1=strcmp(mytiles.dirGrid,mygrid.dirGrid)&(mytiles.nFaces==mygrid.nFaces)&...
35			strcmp(mytiles.fileFormat,mygrid.fileFormat)&(sum(mytiles.ioSize~=mygrid.ioSize)==0);
36	gforget	1.2	test2=(mytiles.ni==ni)&(mytiles.nj==nj);
37	gforget	1.1
38			%2) update tile arrays if not up to date
39			if ~test1\|~test2;
40			%
41			mytiles.dirGrid=mygrid.dirGrid; mytiles.nFaces=mygrid.nFaces;
42			mytiles.fileFormat=mygrid.fileFormat; mytiles.ioSize=mygrid.ioSize;
43			%
44			XC=mygrid.XC; YC=mygrid.YC;
45			XC11=XC; YC11=XC; XCNINJ=XC; YCNINJ=XC; iTile=XC; jTile=XC; tileNo=XC;
46			tileCount=0;
47			for iF=1:XC11.nFaces;
48			face_XC=XC{iF}; face_YC=YC{iF};
49			for ii=1:size(face_XC,1)/ni;
50			for jj=1:size(face_XC,2)/nj;
51			tileCount=tileCount+1;
52			tmp_i=[1:ni]+ni*(ii-1);
53			tmp_j=[1:nj]+nj*(jj-1);
54			tmp_XC=face_XC(tmp_i,tmp_j);
55			tmp_YC=face_YC(tmp_i,tmp_j);
56			XC11{iF}(tmp_i,tmp_j)=tmp_XC(1,1);
57			YC11{iF}(tmp_i,tmp_j)=tmp_YC(1,1);
58			XCNINJ{iF}(tmp_i,tmp_j)=tmp_XC(end,end);
59			YCNINJ{iF}(tmp_i,tmp_j)=tmp_YC(end,end);
60			iTile{iF}(tmp_i,tmp_j)=[1:ni]'*ones(1,nj);
61			jTile{iF}(tmp_i,tmp_j)=ones(ni,1)*[1:nj];
62			tileNo{iF}(tmp_i,tmp_j)=tileCount*ones(ni,nj);
63			end;
64			end;
65			end;
66
67			for iF=1:length(loc_arrays);
68			eval(['mytiles.' loc_arrays{iF} '=convert2array(' loc_arrays{iF} ');']);
69			eval(['clear ' loc_arrays{iF} ';']);
70			end;
71			end;
72
73			%3) decide what is needed based on the number or arguments
74			% (this needs to happen after mytiles is computed, since
75			% I use the same variable names in that bloc and below)
76	gforget	1.2	if nargin==2;
77			%simply return mytiles as loc_tile
78			tileNo=convert2array(mytiles.tileNo);
79			varargout={tileNo};
80			return;
81			elseif nargin==4;
82	gforget	1.1	XC=varargin{1}; YC=varargin{2};
83			tmp1=find(XC>180); XC(tmp1)=XC(tmp1)-360;
84			iTile=[]; jTile=[]; tileNo=[];
85			XC11=[]; YC11=[]; XCNINJ=[]; YCNINJ=[];
86			elseif nargin==5;
87			XC=[]; YC=[];
88			iTile=varargin{1}; jTile=varargin{2}; tileNo=varargin{3};
89			XC11=[]; YC11=[]; XCNINJ=[]; YCNINJ=[];
90			elseif nargin==8;
91			XC=[]; YC=[];
92			iTile=varargin{1}; jTile=varargin{2}; tileNo=[];
93			XC11=varargin{3}; YC11=varargin{4}; XCNINJ=varargin{5}; YCNINJ=varargin{6};
94			else;
95			error('wrong argument list');
96			end;
97
98			%3) find the grid points indices in array format
99			if ~isempty(XC);
100	gforget	1.3	%identify nearest neighbor on the sphere
101			XCgrid=convert2array(mygrid.XC);
102			YCgrid=convert2array(mygrid.YC);
103			x=sin(pi/2-YCgridpi/180).cos(XCgrid*pi/180);
104			y=sin(pi/2-YCgridpi/180).sin(XCgrid*pi/180);
105			z=cos(pi/2-YCgrid*pi/180);
106			xx=sin(pi/2-YCpi/180).cos(XC*pi/180);
107			yy=sin(pi/2-YCpi/180).sin(XC*pi/180);
108			zz=cos(pi/2-YC*pi/180);
109			kk=find(~isnan(x));
110			ik = knnsearch([x(kk) y(kk) z(kk)],[xx yy zz]);
111			ik=kk(ik);
112			%
113			%(old method that uses longitude, latitude directly)
114			%tmp1=find(XC>180); XC(tmp1)=XC(tmp1)-360;
115			%ik=gcmfaces_bindata(XC,YC);
116	gforget	1.1	else;
117			%use mytiles as a look up table
118			ik=zeros(size(iTile));
119			for ikk=1:length(iTile);
120			if ~isempty(tileNo);
121			tmp1=find(mytiles.iTile==iTile(ikk)&mytiles.jTile==jTile(ikk)&mytiles.tileNo==tileNo(ikk));
122			else;
123			tmp1=find(mytiles.iTile==iTile(ikk)&mytiles.jTile==jTile(ikk)&...
124			abs(mytiles.XC-XC(ikk))<1e-4&...
125			abs(mytiles.YC-YC(ikk))<1e-4&...
126			abs(mytiles.XC11-XC11(ikk))<1e-4&...
127			abs(mytiles.YC11-YC11(ikk))<1e-4&...
128			abs(mytiles.XCNINJ-XCNINJ(ikk))<1e-4&...
129			abs(mytiles.YCNINJ-YCNINJ(ikk))<1e-4);
130			end;
131			if isempty(tmp1); error('grid point could not be identified from inputs'); end;
132			ik(ikk)=tmp1(1);
133			end;
134			end;
135
136			%4) complement location arrays and prepare output
137
138			for iF=1:length(loc_arrays);
139			eval(['loc_tile.' loc_arrays{iF} '=mytiles.' loc_arrays{iF} '(ik);']);
140			end;
141
142	gforget	1.2	%5) output result
143			varargout={loc_tile};
144
145
146