matlab_class/gcmfaces_IO/grid_load.m

function []=grid_load(dirGrid,nFaces,fileFormat,memoryLimit);
%object:    load grid information, convert it to gcmfaces format
%           and encapsulate it in the global mygrid structure.
%inputs:    dirGrid is the directory where the grid files (gcm output) can be found.
%           nFaces is the number of faces in this gcm set-up of current interest.
%           fileFormat is the file format ('straight','cube','compact')
%optional:  memoryLimit is a flag that allows the user to omit secondary
%               grid fields in case memory/storage become an issue. It
%               takes 3 values : (0; the default) includes everything
%               (1) omits all 3D fields but hFacC (2) only loads XC & YC.

input_list_check('grid_load',nargin);

if isempty(whos('memoryLimit')); memoryLimit=0; end;

gcmfaces_global; mygrid=[];

mygrid.dirGrid=dirGrid;
mygrid.nFaces=nFaces;
mygrid.fileFormat=fileFormat;
mygrid.gcm2facesFast=false;
mygrid.memoryLimit=memoryLimit;

if mygrid.memoryLimit>0;
    gcmfaces_msg(['* Warning from grid_load : memoryLimit>0 ' ...
        'may precludes advanced gmcfaces functions.'],'');
end;
if mygrid.memoryLimit>1;
    gcmfaces_msg(['* Warning from grid_load : memoryLimit>1 ' ...
        'may only allow for basic fields displays.'],'');
end;

if ~isempty(dir([dirGrid 'grid.specs.mat']));
    specs=open([dirGrid 'grid.specs.mat']);
    mygrid.ioSize=specs.ioSize;
    mygrid.facesSize=specs.facesSize;
    mygrid.facesExpand=specs.facesExpand;
    %example for creating grid.specs.mat, to put in dirGrid :
    %ioSize=[364500 1];
    %facesSize=[[270 450];[0 0];[270 270];[180 270];[450 270]];
    %facesExpand=[270 450];
    %save grid.specs.mat ioSize facesSize facesExpand;
elseif strcmp(fileFormat,'compact');
    v0=rdmds([dirGrid 'XC']);
    mygrid.ioSize=size(v0);
    nn=size(v0,1); pp=size(v0,2)/nn;
    mm=(pp+4-mygrid.nFaces)/4*nn;
    mygrid.facesSize=[[nn mm];[nn mm];[nn nn];[mm nn];[mm nn];[nn nn]];
    mygrid.facesExpand=[];
elseif strcmp(fileFormat,'cube');
    v0=rdmds([dirGrid 'XC']);
    mygrid.ioSize=size(v0);
    nn=size(v0,2);
    mygrid.facesSize=[[nn nn];[nn nn];[nn nn];[nn nn];[nn nn];[nn nn]];
    mygrid.facesExpand=[];
elseif strcmp(fileFormat,'straight');
    v0=rdmds([dirGrid 'XC']);
    mygrid.ioSize=size(v0);
    mygrid.facesSize=mygrid.ioSize;
    mygrid.facesExpand=[];
end;
mygrid.missVal=NaN;%will be set to 0 once the grid has been loaded.

if  ~(nFaces==1&strcmp(fileFormat,'straight'))&...
        ~(nFaces==6&strcmp(fileFormat,'cube'))&...
        ~(nFaces==6&strcmp(fileFormat,'compact'))&...
        ~(nFaces==5&strcmp(fileFormat,'compact'));
    %     if myenv.verbose;
    %         fprintf('\nconvert2gcmfaces.m init: there are several supported file conventions. \n');
    %         fprintf('  By default gcmfaces assumes MITgcm type binary formats as follows: \n')
    %         fprintf('  (1 face) straight global format; (4 or 5 faces) compact global format\n');
    %         fprintf('  (6 faces) cube format with one face after the other. \n');
    %         fprintf('  If this is inadequate, you can change the format below.\n\n');
    %     end;
    error('non-tested topology/fileFormat');
end;

%the various grid fields:
if mygrid.memoryLimit==0;
    list0={'XC','XG','YC','YG','RAC','RAZ','DXC','DXG','DYC','DYG','hFacC','hFacS','hFacW','Depth'};
elseif mygrid.memoryLimit==1;
    list0={'XC','XG','YC','YG','RAC','RAZ','DXC','DXG','DYC','DYG','hFacC','Depth'};
elseif mygrid.memoryLimit==2;
    list0={'XC','YC'};
end;

for iFld=1:length(list0);
    eval(['mygrid.' list0{iFld} '=rdmds2gcmfaces([dirGrid ''' list0{iFld} '*'']);']);
end;

%the vertical grid
list0={'RC','RF','DRC','DRF'};
for iFld=1:length(list0);
    eval(['mygrid.' list0{iFld} '=squeeze(rdmds([dirGrid ''' list0{iFld} '*'']));']);
end;

%grid orientation
if mygrid.memoryLimit<2;
    list0={'AngleCS','AngleSN'};
    test0=~isempty(dir([dirGrid 'AngleCS*']));
    if test0;
        for iFld=1:length(list0);
            eval(['mygrid.' list0{iFld} '=rdmds2gcmfaces([dirGrid ''' list0{iFld} '*'']);']);
        end;
    else;
        warning('\n AngleCS/AngleSN not found; set to 1/0 assuming lat/lon grid.\n');
        mygrid.AngleCS=mygrid.XC; mygrid.AngleCS(:)=1;
        mygrid.AngleSN=mygrid.XC; mygrid.AngleSN(:)=0;
    end;
end;

%if grid is incomplete (blank tiles) then try to get 
%additional info from native grid, or apply missing value.
test0=sum(isnan(mygrid.XC))>0;
test1=prod(mygrid.ioSize)~=sum(isnan(NaN*mygrid.XC(:)));
if test0|test1;
  %treat fields that are part of the native grid
  mygrid1=mygrid; mygrid=[];
  grid_load_native(dirGrid,nFaces,0);
  mygrid1.XC=mygrid.XC; mygrid1.YC=mygrid.YC;
  mygrid1.XC=mygrid.XG; mygrid1.YG=mygrid.YG;
  mygrid1.RAC=mygrid.RAC;
  mygrid=mygrid1;
  %apply missing value for fields that aren't
  list0={'hFacC','hFacS','hFacW','Depth','AngleCS','AngleSN'};
  for ii=1:length(list0); 
    eval(['tmp1=mygrid.' list0{ii} ';']);
    tmp1(isnan(tmp1))=0;
    eval(['mygrid.' list0{ii} '=tmp1;']);
  end;
  %and fix angles if needed
  tmp1=mygrid.AngleCS.^2+mygrid.AngleSN.^2;
  tmp1=1*(tmp1>0.999&tmp1<1.001);
  mygrid.AngleCS(tmp1==0)=1;
  mygrid.AngleSN(tmp1==0)=0;
end;

%get full RAZ (incl. 'extra line and column') needed for e.g. rotational computations
if mygrid.memoryLimit<2;
    grid_load_native_RAZ;
end;

%grid masks
if mygrid.memoryLimit<1;
    mygrid.hFacCsurf=mygrid.hFacC;
    for ff=1:mygrid.hFacC.nFaces; mygrid.hFacCsurf{ff}=mygrid.hFacC{ff}(:,:,1); end;
    
    mskC=mygrid.hFacC; mskC(mskC==0)=NaN; mskC(mskC>0)=1; mygrid.mskC=mskC;
    mskW=mygrid.hFacW; mskW(mskW==0)=NaN; mskW(mskW>0)=1; mygrid.mskW=mskW;
    mskS=mygrid.hFacS; mskS(mskS==0)=NaN; mskS(mskS>0)=1; mygrid.mskS=mskS;
end;

%zonal mean and sections needed for transport computations
if mygrid.memoryLimit<1;
    if ~isfield(mygrid,'mygrid.LATS_MASKS');
        gcmfaces_lines_zonal;
        mygrid.LATS=[mygrid.LATS_MASKS.lat]';
    end;
    if ~isfield(mygrid,'LINES_MASKS');
        [lonPairs,latPairs,names]=line_greatC_TUV_MASKS_v4;
        gcmfaces_lines_transp(lonPairs,latPairs,names);
    end;
end;

%to allow convert2gcmfaces/doFast:
if isempty(mygrid.facesExpand)&mygrid.memoryLimit<2;
    tmp1=convert2gcmfaces(mygrid.XC);
    tmp1(:)=[1:length(tmp1(:))];
    nn=length(tmp1(:));
    mygrid.gcm2faces=convert2gcmfaces(tmp1);
    mygrid.faces2gcmSize=size(tmp1);
    mygrid.faces2gcm=convert2gcmfaces(tmp1);
    for iFace=1:mygrid.nFaces;
        n=length(mygrid.gcm2faces{iFace}(:));
        mygrid.faces2gcm{iFace}=mygrid.gcm2faces{iFace}(:);
        mygrid.gcm2faces{iFace}=sparse([1:n],mygrid.gcm2faces{iFace}(:),ones(1,n),n,nn);
    end;
    mygrid.gcm2facesFast=true;
end;

%reset missVal parameter to 0. 
%Note : this is only used by convert2widefaces, for runs with cropped grids.
%Note : 0 should not be used as a fill for the grid itself (NaN was used).
mygrid.missVal=0;

1	gforget	1.9	function []=grid_load(dirGrid,nFaces,fileFormat,memoryLimit);
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.5
12			input_list_check('grid_load',nargin);
13	gforget	1.1
14	gforget	1.9	if isempty(whos('memoryLimit')); memoryLimit=0; end;
15
16	gforget	1.7	gcmfaces_global; mygrid=[];
17
18	gforget	1.5	mygrid.dirGrid=dirGrid;
19			mygrid.nFaces=nFaces;
20			mygrid.fileFormat=fileFormat;
21	gforget	1.6	mygrid.gcm2facesFast=false;
22	gforget	1.9	mygrid.memoryLimit=memoryLimit;
23
24			if mygrid.memoryLimit>0;
25			gcmfaces_msg(['* Warning from grid_load : memoryLimit>0 ' ...
26			'may precludes advanced gmcfaces functions.'],'');
27			end;
28			if mygrid.memoryLimit>1;
29			gcmfaces_msg(['* Warning from grid_load : memoryLimit>1 ' ...
30			'may only allow for basic fields displays.'],'');
31			end;
32
33			if ~isempty(dir([dirGrid 'grid.specs.mat']));
34			specs=open([dirGrid 'grid.specs.mat']);
35			mygrid.ioSize=specs.ioSize;
36			mygrid.facesSize=specs.facesSize;
37			mygrid.facesExpand=specs.facesExpand;
38			%example for creating grid.specs.mat, to put in dirGrid :
39			%ioSize=[364500 1];
40			%facesSize=[[270 450];[0 0];[270 270];[180 270];[450 270]];
41			%facesExpand=[270 450];
42			%save grid.specs.mat ioSize facesSize facesExpand;
43	gforget	1.10	elseif strcmp(fileFormat,'compact');
44	gforget	1.9	v0=rdmds([dirGrid 'XC']);
45			mygrid.ioSize=size(v0);
46			nn=size(v0,1); pp=size(v0,2)/nn;
47			mm=(pp+4-mygrid.nFaces)/4*nn;
48	gforget	1.10	mygrid.facesSize=[[nn mm];[nn mm];[nn nn];[mm nn];[mm nn];[nn nn]];
49			mygrid.facesExpand=[];
50			elseif strcmp(fileFormat,'cube');
51			v0=rdmds([dirGrid 'XC']);
52			mygrid.ioSize=size(v0);
53			nn=size(v0,2);
54			mygrid.facesSize=[[nn nn];[nn nn];[nn nn];[nn nn];[nn nn];[nn nn]];
55			mygrid.facesExpand=[];
56			elseif strcmp(fileFormat,'straight');
57			v0=rdmds([dirGrid 'XC']);
58			mygrid.ioSize=size(v0);
59			mygrid.facesSize=mygrid.ioSize;
60	gforget	1.9	mygrid.facesExpand=[];
61			end;
62	gforget	1.12	mygrid.missVal=NaN;%will be set to 0 once the grid has been loaded.
63	gforget	1.1
64	gforget	1.5	if ~(nFaces==1&strcmp(fileFormat,'straight'))&...
65			~(nFaces==6&strcmp(fileFormat,'cube'))&...
66			~(nFaces==6&strcmp(fileFormat,'compact'))&...
67			~(nFaces==5&strcmp(fileFormat,'compact'));
68	gforget	1.8	% if myenv.verbose;
69			% fprintf('\nconvert2gcmfaces.m init: there are several supported file conventions. \n');
70			% fprintf(' By default gcmfaces assumes MITgcm type binary formats as follows: \n')
71			% fprintf(' (1 face) straight global format; (4 or 5 faces) compact global format\n');
72			% fprintf(' (6 faces) cube format with one face after the other. \n');
73			% fprintf(' If this is inadequate, you can change the format below.\n\n');
74			% end;
75	gforget	1.5	error('non-tested topology/fileFormat');
76	gforget	1.1	end;
77
78	gforget	1.6	%the various grid fields:
79	gforget	1.9	if mygrid.memoryLimit==0;
80			list0={'XC','XG','YC','YG','RAC','RAZ','DXC','DXG','DYC','DYG','hFacC','hFacS','hFacW','Depth'};
81			elseif mygrid.memoryLimit==1;
82			list0={'XC','XG','YC','YG','RAC','RAZ','DXC','DXG','DYC','DYG','hFacC','Depth'};
83			elseif mygrid.memoryLimit==2;
84			list0={'XC','YC'};
85			end;
86
87	gforget	1.1	for iFld=1:length(list0);
88	gforget	1.8	eval(['mygrid.' list0{iFld} '=rdmds2gcmfaces([dirGrid ''' list0{iFld} '*'']);']);
89	gforget	1.1	end;
90
91	gforget	1.9	%the vertical grid
92			list0={'RC','RF','DRC','DRF'};
93			for iFld=1:length(list0);
94			eval(['mygrid.' list0{iFld} '=squeeze(rdmds([dirGrid ''' list0{iFld} '*'']));']);
95			end;
96
97			%grid orientation
98			if mygrid.memoryLimit<2;
99			list0={'AngleCS','AngleSN'};
100			test0=~isempty(dir([dirGrid 'AngleCS*']));
101			if test0;
102			for iFld=1:length(list0);
103			eval(['mygrid.' list0{iFld} '=rdmds2gcmfaces([dirGrid ''' list0{iFld} '*'']);']);
104			end;
105			else;
106			warning('\n AngleCS/AngleSN not found; set to 1/0 assuming lat/lon grid.\n');
107			mygrid.AngleCS=mygrid.XC; mygrid.AngleCS(:)=1;
108			mygrid.AngleSN=mygrid.XC; mygrid.AngleSN(:)=0;
109	gforget	1.8	end;
110	gforget	1.4	end;
111
112	gforget	1.12	%if grid is incomplete (blank tiles) then try to get
113			%additional info from native grid, or apply missing value.
114			test0=sum(isnan(mygrid.XC))>0;
115			test1=prod(mygrid.ioSize)~=sum(isnan(NaN*mygrid.XC(:)));
116			if test0\|test1;
117			%treat fields that are part of the native grid
118			mygrid1=mygrid; mygrid=[];
119			grid_load_native(dirGrid,nFaces,0);
120			mygrid1.XC=mygrid.XC; mygrid1.YC=mygrid.YC;
121			mygrid1.XC=mygrid.XG; mygrid1.YG=mygrid.YG;
122			mygrid1.RAC=mygrid.RAC;
123			mygrid=mygrid1;
124			%apply missing value for fields that aren't
125			list0={'hFacC','hFacS','hFacW','Depth','AngleCS','AngleSN'};
126			for ii=1:length(list0);
127			eval(['tmp1=mygrid.' list0{ii} ';']);
128			tmp1(isnan(tmp1))=0;
129			eval(['mygrid.' list0{ii} '=tmp1;']);
130			end;
131			%and fix angles if needed
132			tmp1=mygrid.AngleCS.^2+mygrid.AngleSN.^2;
133			tmp1=1*(tmp1>0.999&tmp1<1.001);
134			mygrid.AngleCS(tmp1==0)=1;
135			mygrid.AngleSN(tmp1==0)=0;
136			end;
137
138	gforget	1.11	%get full RAZ (incl. 'extra line and column') needed for e.g. rotational computations
139			if mygrid.memoryLimit<2;
140			grid_load_native_RAZ;
141			end;
142
143	gforget	1.9	%grid masks
144			if mygrid.memoryLimit<1;
145			mygrid.hFacCsurf=mygrid.hFacC;
146			for ff=1:mygrid.hFacC.nFaces; mygrid.hFacCsurf{ff}=mygrid.hFacC{ff}(:,:,1); end;
147
148			mskC=mygrid.hFacC; mskC(mskC==0)=NaN; mskC(mskC>0)=1; mygrid.mskC=mskC;
149			mskW=mygrid.hFacW; mskW(mskW==0)=NaN; mskW(mskW>0)=1; mygrid.mskW=mskW;
150			mskS=mygrid.hFacS; mskS(mskS==0)=NaN; mskS(mskS>0)=1; mygrid.mskS=mskS;
151	gforget	1.1	end;
152
153	gforget	1.9	%zonal mean and sections needed for transport computations
154			if mygrid.memoryLimit<1;
155			if ~isfield(mygrid,'mygrid.LATS_MASKS');
156			gcmfaces_lines_zonal;
157			mygrid.LATS=[mygrid.LATS_MASKS.lat]';
158			end;
159			if ~isfield(mygrid,'LINES_MASKS');
160			[lonPairs,latPairs,names]=line_greatC_TUV_MASKS_v4;
161			gcmfaces_lines_transp(lonPairs,latPairs,names);
162			end;
163			end;
164	gforget	1.2
165	gforget	1.8	%to allow convert2gcmfaces/doFast:
166	gforget	1.9	if isempty(mygrid.facesExpand)&mygrid.memoryLimit<2;
167			tmp1=convert2gcmfaces(mygrid.XC);
168			tmp1(:)=[1:length(tmp1(:))];
169			nn=length(tmp1(:));
170			mygrid.gcm2faces=convert2gcmfaces(tmp1);
171			mygrid.faces2gcmSize=size(tmp1);
172			mygrid.faces2gcm=convert2gcmfaces(tmp1);
173			for iFace=1:mygrid.nFaces;
174			n=length(mygrid.gcm2faces{iFace}(:));
175			mygrid.faces2gcm{iFace}=mygrid.gcm2faces{iFace}(:);
176			mygrid.gcm2faces{iFace}=sparse([1:n],mygrid.gcm2faces{iFace}(:),ones(1,n),n,nn);
177			end;
178			mygrid.gcm2facesFast=true;
179	gforget	1.6	end;
180
181	gforget	1.12	%reset missVal parameter to 0.
182			%Note : this is only used by convert2widefaces, for runs with cropped grids.
183			%Note : 0 should not be used as a fill for the grid itself (NaN was used).
184			mygrid.missVal=0;
185	gforget	1.1