/[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.20 - (hide annotations) (download)
Sun Jan 17 14:34:22 2016 UTC (9 years, 6 months ago) by gforget
Branch: MAIN
CVS Tags: checkpoint65x, checkpoint65v, checkpoint65w, checkpoint65t, checkpoint65u
Changes since 1.19: +13 -4 lines
- revise treatment of case when neither ''nctiles_grid/''
  nor ''GRID/'' was readily found

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

  ViewVC Help
Powered by ViewVC 1.1.22