/[MITgcm]/MITgcm_contrib/gael/matlab_class/gcmfaces_calc/gcmfaces_loc_tile.m
ViewVC logotype

Annotation of /MITgcm_contrib/gael/matlab_class/gcmfaces_calc/gcmfaces_loc_tile.m

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


Revision 1.1 - (hide annotations) (download)
Mon Oct 1 23:35:51 2012 UTC (12 years, 9 months ago) by gforget
Branch: MAIN
- gcmfaces_bindata.m
  implement consistency test between mytri and mygrid; then
  recompute mytri if inconsistency, or if it is empty.
- gcmfaces_loc_tile.m (new)
  complement the location of data points on the grid. Starting
  with either lon/lat location or grid tile indices, returns both.

1 gforget 1.1 function [loc_tile]=gcmfaces_loc_tile(ni,nj,varargin);
2     %object : complement the location of data points on the grid. Starting
3     % with either lon/lat location or grid tile indices, returns both.
4     % Inputs are assumed to be 1D arrays.
5     %input : ni,nj is the MITgcm tile size (2 numbers total)
6     % 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     %output : XC,YC, tileNo, XC11, YC11, XCNINJ, YCNINJ
17     %
18     %pre-requisite : the grid of interest must have been loaded with grid_load
19     %
20     %notes : - if XC, YC are provided as input, then loc_tile.XC, YC is nearest neighbor
21     % - by assumption the tile number increase from w to e, then
22     % from s to n, the from with the face number
23    
24     gcmfaces_global;
25     global mytiles;
26     loc_arrays={'XC','YC','XC11','YC11','XCNINJ','YCNINJ','iTile','jTile','tileNo'};
27    
28     %1) check that tile arrays are up to date
29     if isempty(mytiles);
30     mytiles.dirGrid=mygrid.dirGrid; mytiles.nFaces=mygrid.nFaces;
31     mytiles.fileFormat=mygrid.fileFormat; mytiles.ioSize=mygrid.ioSize;
32     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     test2=isfield(mytiles,'XC11');
37    
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     if nargin==4;
77     XC=varargin{1}; YC=varargin{2};
78     tmp1=find(XC>180); XC(tmp1)=XC(tmp1)-360;
79     iTile=[]; jTile=[]; tileNo=[];
80     XC11=[]; YC11=[]; XCNINJ=[]; YCNINJ=[];
81     elseif nargin==5;
82     XC=[]; YC=[];
83     iTile=varargin{1}; jTile=varargin{2}; tileNo=varargin{3};
84     XC11=[]; YC11=[]; XCNINJ=[]; YCNINJ=[];
85     elseif nargin==8;
86     XC=[]; YC=[];
87     iTile=varargin{1}; jTile=varargin{2}; tileNo=[];
88     XC11=varargin{3}; YC11=varargin{4}; XCNINJ=varargin{5}; YCNINJ=varargin{6};
89     else;
90     error('wrong argument list');
91     end;
92    
93     %3) find the grid points indices in array format
94     if ~isempty(XC);
95     %triangulate nearest neighbor
96     tmp1=find(XC>180); XC(tmp1)=XC(tmp1)-360;
97     ik=gcmfaces_bindata(XC,YC);
98     else;
99     %use mytiles as a look up table
100     ik=zeros(size(iTile));
101     for ikk=1:length(iTile);
102     if ~isempty(tileNo);
103     tmp1=find(mytiles.iTile==iTile(ikk)&mytiles.jTile==jTile(ikk)&mytiles.tileNo==tileNo(ikk));
104     else;
105     tmp1=find(mytiles.iTile==iTile(ikk)&mytiles.jTile==jTile(ikk)&...
106     abs(mytiles.XC-XC(ikk))<1e-4&...
107     abs(mytiles.YC-YC(ikk))<1e-4&...
108     abs(mytiles.XC11-XC11(ikk))<1e-4&...
109     abs(mytiles.YC11-YC11(ikk))<1e-4&...
110     abs(mytiles.XCNINJ-XCNINJ(ikk))<1e-4&...
111     abs(mytiles.YCNINJ-YCNINJ(ikk))<1e-4);
112     end;
113     if isempty(tmp1); error('grid point could not be identified from inputs'); end;
114     ik(ikk)=tmp1(1);
115     end;
116     end;
117    
118     %4) complement location arrays and prepare output
119    
120     for iF=1:length(loc_arrays);
121     eval(['loc_tile.' loc_arrays{iF} '=mytiles.' loc_arrays{iF} '(ik);']);
122     end;
123    

  ViewVC Help
Powered by ViewVC 1.1.22