/[MITgcm]/MITgcm_contrib/gael/matlab_class/gcmfaces_IO/grid_load.m
ViewVC logotype

Annotation of /MITgcm_contrib/gael/matlab_class/gcmfaces_IO/grid_load.m

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


Revision 1.19 - (hide annotations) (download)
Wed Jan 13 16:24:07 2016 UTC (9 years, 6 months ago) by gforget
Branch: MAIN
Changes since 1.18: +1 -1 lines
- replace sample_analysis/line_greatC_TUV_MASKS_v4.m
     with gcmfaces_calc/gcmfaces_lines_pairs.m

1 gforget 1.13 function []=grid_load(dirGrid,nFaces,fileFormat,memoryLimit,omitNativeGrid);
2 gforget 1.5 %object: load grid information, convert it to gcmfaces format
3     % and encapsulate it in the global mygrid structure.
4     %inputs: dirGrid is the directory where the grid files (gcm output) can be found.
5     % nFaces is the number of faces in this gcm set-up of current interest.
6     % fileFormat is the file format ('straight','cube','compact')
7 gforget 1.9 %optional: memoryLimit is a flag that allows the user to omit secondary
8     % grid fields in case memory/storage become an issue. It
9     % takes 3 values : (0; the default) includes everything
10     % (1) omits all 3D fields but hFacC (2) only loads XC & YC.
11 gforget 1.13 % omitNativeGrid is a flag (0 by default) to bypass (when flag 1)
12     % grid_load_native and grid_load_native_RAZ calls (that can
13     % complement mygrid based upon e.g. tile*.mitgrid files)
14 gforget 1.1
15 gforget 1.9 if isempty(whos('memoryLimit')); memoryLimit=0; end;
16 gforget 1.13 if isempty(whos('omitNativeGrid')); omitNativeGrid=0; end;
17 gforget 1.9
18 gforget 1.7 gcmfaces_global; mygrid=[];
19    
20 gforget 1.15 if nargin==0;%assume that LLC90 grid is used
21     if ~isempty(dir('nctiles_grid'));
22 gforget 1.18 mygrid.dirGrid=['nctiles_grid' filesep];
23 gforget 1.15 mygrid.fileFormat='nctiles';
24     elseif ~isempty(dir('GRID'));
25 gforget 1.18 mygrid.dirGrid=['GRID' filesep];
26 gforget 1.15 mygrid.fileFormat='compact';
27     else;
28 gforget 1.16 fprintf('\n The ECCO v4 grid can be obtained as follows: \n');
29     fprintf([' wget --recursive ftp://mit.ecco-group.org/ecco_for_las/' ...
30     'version_4/release1/nctiles_grid \n']);
31     fprintf(' mv mit.ecco-group.org/ecco_for_las/version_4/release1/nctiles_grid . \n\n');
32 gforget 1.15 error('could not find either ''nctiles_grid/'' or ''GRID/''');
33     end;
34     mygrid.nFaces=5;
35     mygrid.gcm2facesFast=false;
36     mygrid.memoryLimit=0;
37     omitNativeGrid=isempty(dir([mygrid.dirGrid 'tile001.mitgrid']));
38     else;%user specified settings
39     mygrid.dirGrid=dirGrid;
40     mygrid.nFaces=nFaces;
41     mygrid.fileFormat=fileFormat;
42     mygrid.gcm2facesFast=false;
43     mygrid.memoryLimit=memoryLimit;
44     end;
45 gforget 1.9
46     if mygrid.memoryLimit>0;
47     gcmfaces_msg(['* Warning from grid_load : memoryLimit>0 ' ...
48     'may precludes advanced gmcfaces functions.'],'');
49     end;
50     if mygrid.memoryLimit>1;
51     gcmfaces_msg(['* Warning from grid_load : memoryLimit>1 ' ...
52     'may only allow for basic fields displays.'],'');
53     end;
54    
55 gforget 1.15 if strcmp(mygrid.fileFormat,'nctiles');
56     mygrid.ioSize=[90 1170];
57     mygrid.facesSize=[[90 270];[90 270];[90 90];[270 90];[270 90];[90 90]];
58     mygrid.facesExpand=[];
59     elseif ~isempty(dir([mygrid.dirGrid 'grid.specs.mat']));
60     specs=open([mygrid.dirGrid 'grid.specs.mat']);
61 gforget 1.9 mygrid.ioSize=specs.ioSize;
62     mygrid.facesSize=specs.facesSize;
63     mygrid.facesExpand=specs.facesExpand;
64     %example for creating grid.specs.mat, to put in dirGrid :
65     %ioSize=[364500 1];
66     %facesSize=[[270 450];[0 0];[270 270];[180 270];[450 270]];
67     %facesExpand=[270 450];
68     %save grid.specs.mat ioSize facesSize facesExpand;
69 gforget 1.15 elseif strcmp(mygrid.fileFormat,'compact');
70     v0=rdmds([mygrid.dirGrid 'XC']);
71 gforget 1.9 mygrid.ioSize=size(v0);
72     nn=size(v0,1); pp=size(v0,2)/nn;
73     mm=(pp+4-mygrid.nFaces)/4*nn;
74 gforget 1.10 mygrid.facesSize=[[nn mm];[nn mm];[nn nn];[mm nn];[mm nn];[nn nn]];
75     mygrid.facesExpand=[];
76 gforget 1.15 elseif strcmp(mygrid.fileFormat,'cube');
77     v0=rdmds([mygrid.dirGrid 'XC']);
78 gforget 1.10 mygrid.ioSize=size(v0);
79     nn=size(v0,2);
80     mygrid.facesSize=[[nn nn];[nn nn];[nn nn];[nn nn];[nn nn];[nn nn]];
81     mygrid.facesExpand=[];
82 gforget 1.15 elseif strcmp(mygrid.fileFormat,'straight');
83     v0=rdmds([mygrid.dirGrid 'XC']);
84 gforget 1.10 mygrid.ioSize=size(v0);
85     mygrid.facesSize=mygrid.ioSize;
86 gforget 1.9 mygrid.facesExpand=[];
87     end;
88 gforget 1.15
89 gforget 1.12 mygrid.missVal=NaN;%will be set to 0 once the grid has been loaded.
90 gforget 1.1
91 gforget 1.15 if ~(mygrid.nFaces==1&strcmp(mygrid.fileFormat,'straight'))&...
92     ~(mygrid.nFaces==6&strcmp(mygrid.fileFormat,'cube'))&...
93     ~(mygrid.nFaces==6&strcmp(mygrid.fileFormat,'compact'))&...
94     ~(mygrid.nFaces==5&strcmp(mygrid.fileFormat,'compact'))&...
95     ~(mygrid.nFaces==5&strcmp(mygrid.fileFormat,'nctiles'));
96     error('non-supported grid topology');
97 gforget 1.1 end;
98    
99 gforget 1.17 if strcmp(mygrid.fileFormat,'nctiles');
100     %place holders (needed for read_nctiles)
101     tmp1=NaN*zeros(90,270);
102     tmp2=NaN*zeros(90,90);
103     tmp3=NaN*zeros(270,90);
104     mygrid.XC=gcmfaces({tmp1,tmp1,tmp2,tmp3,tmp3});
105     mygrid.YC=gcmfaces({tmp1,tmp1,tmp2,tmp3,tmp3});
106     mygrid.RC=NaN*zeros(50,1);
107     clear tmp?;
108     end;
109    
110 gforget 1.6 %the various grid fields:
111 gforget 1.9 if mygrid.memoryLimit==0;
112     list0={'XC','XG','YC','YG','RAC','RAZ','DXC','DXG','DYC','DYG','hFacC','hFacS','hFacW','Depth'};
113     elseif mygrid.memoryLimit==1;
114     list0={'XC','XG','YC','YG','RAC','RAZ','DXC','DXG','DYC','DYG','hFacC','Depth'};
115     elseif mygrid.memoryLimit==2;
116     list0={'XC','YC'};
117     end;
118    
119 gforget 1.1 for iFld=1:length(list0);
120 gforget 1.15 if ~strcmp(mygrid.fileFormat,'nctiles');
121     eval(['mygrid.' list0{iFld} '=rdmds2gcmfaces([mygrid.dirGrid ''' list0{iFld} '*'']);']);
122     else;
123     eval(['mygrid.' list0{iFld} '=read_nctiles([mygrid.dirGrid ''GRID''],''' list0{iFld} ''');']);
124     end;
125 gforget 1.1 end;
126    
127 gforget 1.9 %the vertical grid
128     list0={'RC','RF','DRC','DRF'};
129     for iFld=1:length(list0);
130 gforget 1.15 if ~strcmp(mygrid.fileFormat,'nctiles');
131     eval(['mygrid.' list0{iFld} '=squeeze(rdmds([mygrid.dirGrid ''' list0{iFld} '*'']));']);
132     else;
133     eval(['ncload ' mygrid.dirGrid 'GRID.0001.nc ' list0{iFld} ';']);
134     eval(['mygrid.' list0{iFld} '=' list0{iFld} ''';']);
135     if strcmp(list0{iFld},'RF'); mygrid.RF=[mygrid.RF;NaN]; end;
136     end;
137 gforget 1.9 end;
138    
139     %grid orientation
140     if mygrid.memoryLimit<2;
141     list0={'AngleCS','AngleSN'};
142 gforget 1.15 test0=~isempty(dir([mygrid.dirGrid 'AngleCS*']));
143     if strcmp(mygrid.fileFormat,'nctiles');
144     for iFld=1:length(list0);
145     eval(['mygrid.' list0{iFld} '=read_nctiles([mygrid.dirGrid ''GRID''],''' list0{iFld} ''');']);
146     end;
147     elseif test0;
148 gforget 1.9 for iFld=1:length(list0);
149 gforget 1.15 eval(['mygrid.' list0{iFld} '=rdmds2gcmfaces([mygrid.dirGrid ''' list0{iFld} '*'']);']);
150 gforget 1.9 end;
151     else;
152     warning('\n AngleCS/AngleSN not found; set to 1/0 assuming lat/lon grid.\n');
153     mygrid.AngleCS=mygrid.XC; mygrid.AngleCS(:)=1;
154     mygrid.AngleSN=mygrid.XC; mygrid.AngleSN(:)=0;
155 gforget 1.8 end;
156 gforget 1.4 end;
157    
158 gforget 1.12 %if grid is incomplete (blank tiles) then try to get
159     %additional info from native grid, or apply missing value.
160     test0=sum(isnan(mygrid.XC))>0;
161     test1=prod(mygrid.ioSize)~=sum(isnan(NaN*mygrid.XC(:)));
162 gforget 1.13 if (test0|test1)&~omitNativeGrid;
163 gforget 1.12 %treat fields that are part of the native grid
164     mygrid1=mygrid; mygrid=[];
165 gforget 1.15 grid_load_native(mygrid.dirGrid,nFaces,0);
166 gforget 1.12 mygrid1.XC=mygrid.XC; mygrid1.YC=mygrid.YC;
167 gforget 1.14 mygrid1.XG=mygrid.XG; mygrid1.YG=mygrid.YG;
168 gforget 1.12 mygrid1.RAC=mygrid.RAC;
169     mygrid=mygrid1;
170     %apply missing value for fields that aren't
171     list0={'hFacC','hFacS','hFacW','Depth','AngleCS','AngleSN'};
172     for ii=1:length(list0);
173     eval(['tmp1=mygrid.' list0{ii} ';']);
174     tmp1(isnan(tmp1))=0;
175     eval(['mygrid.' list0{ii} '=tmp1;']);
176     end;
177     %and fix angles if needed
178     tmp1=mygrid.AngleCS.^2+mygrid.AngleSN.^2;
179     tmp1=1*(tmp1>0.999&tmp1<1.001);
180     mygrid.AngleCS(tmp1==0)=1;
181     mygrid.AngleSN(tmp1==0)=0;
182     end;
183    
184 gforget 1.11 %get full RAZ (incl. 'extra line and column') needed for e.g. rotational computations
185 gforget 1.13 if mygrid.memoryLimit<2&~omitNativeGrid;
186 gforget 1.11 grid_load_native_RAZ;
187     end;
188    
189 gforget 1.9 %grid masks
190     if mygrid.memoryLimit<1;
191     mygrid.hFacCsurf=mygrid.hFacC;
192     for ff=1:mygrid.hFacC.nFaces; mygrid.hFacCsurf{ff}=mygrid.hFacC{ff}(:,:,1); end;
193    
194     mskC=mygrid.hFacC; mskC(mskC==0)=NaN; mskC(mskC>0)=1; mygrid.mskC=mskC;
195     mskW=mygrid.hFacW; mskW(mskW==0)=NaN; mskW(mskW>0)=1; mygrid.mskW=mskW;
196     mskS=mygrid.hFacS; mskS(mskS==0)=NaN; mskS(mskS>0)=1; mygrid.mskS=mskS;
197 gforget 1.1 end;
198    
199 gforget 1.9 %zonal mean and sections needed for transport computations
200     if mygrid.memoryLimit<1;
201     if ~isfield(mygrid,'mygrid.LATS_MASKS');
202     gcmfaces_lines_zonal;
203     mygrid.LATS=[mygrid.LATS_MASKS.lat]';
204     end;
205     if ~isfield(mygrid,'LINES_MASKS');
206 gforget 1.19 [lonPairs,latPairs,names]=gcmfaces_lines_pairs;
207 gforget 1.9 gcmfaces_lines_transp(lonPairs,latPairs,names);
208     end;
209     end;
210 gforget 1.2
211 gforget 1.8 %to allow convert2gcmfaces/doFast:
212 gforget 1.9 if isempty(mygrid.facesExpand)&mygrid.memoryLimit<2;
213     tmp1=convert2gcmfaces(mygrid.XC);
214     tmp1(:)=[1:length(tmp1(:))];
215     nn=length(tmp1(:));
216     mygrid.gcm2faces=convert2gcmfaces(tmp1);
217     mygrid.faces2gcmSize=size(tmp1);
218     mygrid.faces2gcm=convert2gcmfaces(tmp1);
219     for iFace=1:mygrid.nFaces;
220     n=length(mygrid.gcm2faces{iFace}(:));
221     mygrid.faces2gcm{iFace}=mygrid.gcm2faces{iFace}(:);
222     mygrid.gcm2faces{iFace}=sparse([1:n],mygrid.gcm2faces{iFace}(:),ones(1,n),n,nn);
223     end;
224     mygrid.gcm2facesFast=true;
225 gforget 1.6 end;
226    
227 gforget 1.12 %reset missVal parameter to 0.
228     %Note : this is only used by convert2widefaces, for runs with cropped grids.
229     %Note : 0 should not be used as a fill for the grid itself (NaN was used).
230     mygrid.missVal=0;
231 gforget 1.1
232 gforget 1.18 mygrid.dirGrid=[pwd filesep mygrid.dirGrid];
233    

  ViewVC Help
Powered by ViewVC 1.1.22