matlab_class/gcmfaces_IO/write2nctiles.m

function [dimOut]=write2nctiles(fileOut,fldIn,doCreate,varargin);
%object : write gcmfaces object to netcdf files (tiled)
%inputs : fileOut is the name of the file to be created
%         fldIn is the array to write to disk
%         doCreate is a 0/1 switch; 1 => create file ; 0 => append to file. 
%optional paramaters :
%         can be provided in the form {'name',value}
%         those that are currently active are
%               'descr' is the file description (default '').
%               'rdm' is the extended estimate description (default '').
%               'fldName' is the nc variable name for fld (default : the outside name of fldIn).
%               'longName' is the corresponding long name (default : '').
%               'units' is the unit of fld (default : '(unknown)').
%               'missval' is the missing value (default : NaN).
%               'fillval' is the fill value (default : NaN).
%               'tileNo' is a map of tile indices (default is face number)
%               'coord' is auxilliary coordinates attribute (e.g. 'lon lat dep')
%               'dimsize' is the array size associated with 'coord'
%               'xtype' is 'double' by default 
%netcdf dimensions : array dimensions are simply set to 'i1,i2,...'

gcmfaces_global;
if ~(myenv.useNativeMatlabNetcdf);
  error('only native matlab nectdf is supported in write2nctiles');
end;

doCheck=0;%set to one to print stuff to screen

fldInIsaGcmfaces=isa(fldIn,'gcmfaces');

if fldInIsaGcmfaces;
  fldInIsEmpty=isempty(fldIn{1});
else;
  fldInIsEmpty=isempty(fldIn);
end;

%set more optional paramaters to default values
descr=''; 
rdm='';
fldName=inputname(2); longName=''; 
units='(unknown)'; missval=NaN; fillval=NaN; dimIn=[];
tileNo=mygrid.XC; for ff=1:mygrid.nFaces; tileNo{ff}(:)=ff; end;
coord=''; dimsize=[]; xtype='double';

%set more optional paramaters to user defined values
for ii=1:nargin-3;
    if ~iscell(varargin{ii});
        warning('inputCheck:write2nctiles_1',...
            ['write2nctiles expects \n'...
            '         its optional parameters as cell arrays. \n'...
            '         Argument no. ' num2str(ii+1) ' was ignored \n'...
            '         Type ''help write2nctiles'' for details.']);
    elseif ~ischar(varargin{ii}{1});
        warning('inputCheck:write2nctiles_2',...
            ['write2nctiles expects \n'...
            '         its optional parameters cell arrays \n'...
            '         to start with character string. \n'...
            '         Argument no. ' num2str(ii+1) ' was ignored \n'...
            '         Type ''help write2nctiles'' for details.']);
    else;
        if strcmp(varargin{ii}{1},'descr')|...
                strcmp(varargin{ii}{1},'rdm')|...
                strcmp(varargin{ii}{1},'fldName')|...
                strcmp(varargin{ii}{1},'longName')|...
                strcmp(varargin{ii}{1},'units')|...
                strcmp(varargin{ii}{1},'missval')|...
                strcmp(varargin{ii}{1},'fillval')|...
                strcmp(varargin{ii}{1},'tileNo')|...
                strcmp(varargin{ii}{1},'coord')|...
                strcmp(varargin{ii}{1},'dimsize')|...
                strcmp(varargin{ii}{1},'xtype')|...
                strcmp(varargin{ii}{1},'dimIn');
            eval([varargin{ii}{1} '=varargin{ii}{2};']);
        else;
            warning('inputCheck:write2nctiles_3',...
                ['unknown option ''' varargin{ii}{1} ''' was ignored']);
        end;
    end;
end;

%split fldIn (if isa gcmfaces) into tiles
tileList=unique(convert2vector(tileNo)); 
tileList=tileList(~isnan(tileList));
ntile=length(tileList);
%
if fldInIsaGcmfaces&~fldInIsEmpty;
 for ff=1:ntile; 
   tmp1=[];
   for gg=1:mygrid.nFaces;
     [tmpi,tmpj]=find(tileNo{gg}==ff);
     if ~isempty(tmpi);
       tmpi=[min(tmpi(:)):max(tmpi(:))];
       tmpj=[min(tmpj(:)):max(tmpj(:))];
       tmp1=fldIn{gg}(tmpi,tmpj,:,:);
     end;
   end;
   fldTiles{ff}=tmp1;
 end;
 clear fldIn;
elseif fldInIsaGcmfaces;
 for ff=1:ntile;
   fldTiles{ff}=[];
 end;
end;

%start processing loop
for ff=1:ntile;

if fldInIsaGcmfaces;
  fldTile=fldTiles{ff};
  %reverse order of dimensions
  nn=length(size(fldTile));
  fldTile=permute(fldTile,[nn:-1:1]);
  if ntile==1; clear fldTiles; end;
else;
  fldTile=fldIn;
end;
fileTile=[fileOut sprintf('.%04d.nc',ff)];

%select dimensions of relevance:
if ~fldInIsEmpty;
  nDim=length(size(fldTile));
  dimsize=size(fldTile);
elseif ~isempty(dimsize);
  nDim=length(dimsize);
else;
  error('undertermined array size');
end;

for iDim=1:nDim;
if dimsize(iDim)~=1;
    dimlist{iDim}=['i' num2str(iDim)];
    dimName{iDim}=['array index ' num2str(iDim)];
    eval(['dimvec.i' num2str(iDim) '=[1:dimsize(iDim)];']);
end;
end;

%omit singleton dimensions:
ii=find(dimsize~=1);
dimsize=dimsize(ii);
dimlist={dimlist{ii}};
dimName={dimName{ii}};

%check : 
if doCheck;
whos fldTile
descr
fldName
longName
units
missval
fillval
dimlist
dimName
dimsize
dimvec
keyboard;
end;

%output dimension information
dimOut{ff}=dimlist;

if doCreate;
  %create netcdf file:
  %-------------------
  if prod(dimsize)*4/1e9<1.5;%use (always available) basic netcdf:
    mode='clobber';
  else;%to allow for large file:
    mode='NETCDF4';
  end;
  ncid=nccreate(fileTile,mode);
  nc_global=netcdf.getConstant('NC_GLOBAL');

  if ~isempty(rdm); 
    descr2=[descr ' -- ' rdm{1}];
  else;
    descr2=descr;
  end;
  ncputAtt(ncid,'','description',descr2);
  for pp=2:length(rdm);
  tmp1=char(pp+63);
  netcdf.putAtt(ncid,nc_global,tmp1,rdm{pp});
  end;
  %append readme
  pp=length(rdm)+1;
  netcdf.putAtt(ncid,nc_global,tmp1,'file created using gcmfaces_IO/write2nctiles.m');
  ncputAtt(ncid,'','date',date);
  netcdf.putAtt(ncid,nc_global,'Conventions','CF-1.6')   

  ncputAtt(ncid,'','_FillValue',fillval);
  ncputAtt(ncid,'','missing_value',missval);

  ncdefDim(ncid,'itxt',30);
  for dd=1:length(dimlist); ncdefDim(ncid,dimlist{dd},dimsize(dd)); end;

  for dd=1:length(dimlist);
    ncdefVar(ncid,dimlist{dd},'double',{dimlist{dd}});
    ncputAtt(ncid,dimlist{dd},'long_name',dimName{dd});
    ncputAtt(ncid,dimlist{dd},'units','1');
  end;
  ncclose(ncid);

  %fill in the dimensions dimensions values vectors:
  %-------------------------------------------------
  ncid=ncopen(fileTile,'write');
  for dd=1:length(dimlist);
    ncputvar(ncid,dimlist{dd},getfield(dimvec,dimlist{dd}));
  end;
  ncclose(ncid);
end;

%use dimentsion specified by user
if ~isempty(dimIn); dimlist=dimIn{ff}; end;

%define and fill field:
%----------------------
ncid=ncopen(fileTile,'write');
%
netcdf.reDef(ncid);
ncdefVar(ncid,fldName,xtype,flipdim(dimlist,2));%note the direction flip
if ~isempty(longName); ncputAtt(ncid,fldName,'long_name',longName); end;
if ~isempty(units); ncputAtt(ncid,fldName,'units',units); end;
if ~isempty(coord); ncputAtt(ncid,fldName,'coordinates',coord); end; 
if strcmp(fldName,'lon'); ncputAtt(ncid,fldName,'standard_name','longitude'); end;
if strcmp(fldName,'lat'); ncputAtt(ncid,fldName,'standard_name','latitude'); end;
if strcmp(fldName,'dep'); ncputAtt(ncid,fldName,'standard_name','depth'); end;
if strcmp(fldName,'tim'); ncputAtt(ncid,fldName,'standard_name','time'); end;
if strcmp(fldName,'land'); ncputAtt(ncid,fldName,'standard_name','land_binary_mask'); end;
if strcmp(fldName,'area'); ncputAtt(ncid,fldName,'standard_name','cell_area'); end;
if strcmp(fldName,'thic'); ncputAtt(ncid,fldName,'standard_name','cell_thickness'); end;
netcdf.endDef(ncid);
%
if ~fldInIsEmpty; ncputvar(ncid,fldName,fldTile); end;
%
ncclose(ncid);

end;%for ff=1:ntile;

1	gforget	1.1	function [dimOut]=write2nctiles(fileOut,fldIn,doCreate,varargin);
2			%object : write gcmfaces object to netcdf files (tiled)
3			%inputs : fileOut is the name of the file to be created
4			% fldIn is the array to write to disk
5			% doCreate is a 0/1 switch; 1 => create file ; 0 => append to file.
6			%optional paramaters :
7			% can be provided in the form {'name',value}
8			% those that are currently active are
9			% 'descr' is the file description (default '').
10	gforget	1.2	% 'rdm' is the extended estimate description (default '').
11	gforget	1.1	% 'fldName' is the nc variable name for fld (default : the outside name of fldIn).
12			% 'longName' is the corresponding long name (default : '').
13			% 'units' is the unit of fld (default : '(unknown)').
14			% 'missval' is the missing value (default : NaN).
15			% 'fillval' is the fill value (default : NaN).
16	gforget	1.4	% 'tileNo' is a map of tile indices (default is face number)
17	gforget	1.5	% 'coord' is auxilliary coordinates attribute (e.g. 'lon lat dep')
18	gforget	1.8	% 'dimsize' is the array size associated with 'coord'
19			% 'xtype' is 'double' by default
20	gforget	1.1	%netcdf dimensions : array dimensions are simply set to 'i1,i2,...'
21
22			gcmfaces_global;
23	gforget	1.5	if ~(myenv.useNativeMatlabNetcdf);
24			error('only native matlab nectdf is supported in write2nctiles');
25			end;
26	gforget	1.1
27			doCheck=0;%set to one to print stuff to screen
28
29	gforget	1.7	fldInIsaGcmfaces=isa(fldIn,'gcmfaces');
30
31	gforget	1.8	if fldInIsaGcmfaces;
32			fldInIsEmpty=isempty(fldIn{1});
33			else;
34			fldInIsEmpty=isempty(fldIn);
35			end;
36
37	gforget	1.1	%set more optional paramaters to default values
38			descr='';
39	gforget	1.2	rdm='';
40	gforget	1.1	fldName=inputname(2); longName='';
41	gforget	1.4	units='(unknown)'; missval=NaN; fillval=NaN; dimIn=[];
42			tileNo=mygrid.XC; for ff=1:mygrid.nFaces; tileNo{ff}(:)=ff; end;
43	gforget	1.8	coord=''; dimsize=[]; xtype='double';
44	gforget	1.1
45			%set more optional paramaters to user defined values
46			for ii=1:nargin-3;
47			if ~iscell(varargin{ii});
48			warning('inputCheck:write2nctiles_1',...
49			['write2nctiles expects \n'...
50			' its optional parameters as cell arrays. \n'...
51			' Argument no. ' num2str(ii+1) ' was ignored \n'...
52			' Type ''help write2nctiles'' for details.']);
53			elseif ~ischar(varargin{ii}{1});
54			warning('inputCheck:write2nctiles_2',...
55			['write2nctiles expects \n'...
56			' its optional parameters cell arrays \n'...
57			' to start with character string. \n'...
58			' Argument no. ' num2str(ii+1) ' was ignored \n'...
59			' Type ''help write2nctiles'' for details.']);
60			else;
61			if strcmp(varargin{ii}{1},'descr')\|...
62	gforget	1.2	strcmp(varargin{ii}{1},'rdm')\|...
63	gforget	1.1	strcmp(varargin{ii}{1},'fldName')\|...
64			strcmp(varargin{ii}{1},'longName')\|...
65			strcmp(varargin{ii}{1},'units')\|...
66			strcmp(varargin{ii}{1},'missval')\|...
67			strcmp(varargin{ii}{1},'fillval')\|...
68	gforget	1.4	strcmp(varargin{ii}{1},'tileNo')\|...
69	gforget	1.5	strcmp(varargin{ii}{1},'coord')\|...
70	gforget	1.8	strcmp(varargin{ii}{1},'dimsize')\|...
71			strcmp(varargin{ii}{1},'xtype')\|...
72	gforget	1.1	strcmp(varargin{ii}{1},'dimIn');
73			eval([varargin{ii}{1} '=varargin{ii}{2};']);
74			else;
75			warning('inputCheck:write2nctiles_3',...
76			['unknown option ''' varargin{ii}{1} ''' was ignored']);
77			end;
78			end;
79			end;
80
81	gforget	1.4	%split fldIn (if isa gcmfaces) into tiles
82			tileList=unique(convert2vector(tileNo));
83			tileList=tileList(~isnan(tileList));
84			ntile=length(tileList);
85			%
86	gforget	1.8	if fldInIsaGcmfaces&~fldInIsEmpty;
87	gforget	1.4	for ff=1:ntile;
88			tmp1=[];
89			for gg=1:mygrid.nFaces;
90			[tmpi,tmpj]=find(tileNo{gg}==ff);
91			if ~isempty(tmpi);
92			tmpi=[min(tmpi(:)):max(tmpi(:))];
93			tmpj=[min(tmpj(:)):max(tmpj(:))];
94			tmp1=fldIn{gg}(tmpi,tmpj,:,:);
95			end;
96			end;
97			fldTiles{ff}=tmp1;
98			end;
99	gforget	1.7	clear fldIn;
100	gforget	1.8	elseif fldInIsaGcmfaces;
101			for ff=1:ntile;
102			fldTiles{ff}=[];
103			end;
104	gforget	1.4	end;
105
106			%start processing loop
107			for ff=1:ntile;
108	gforget	1.1
109	gforget	1.7	if fldInIsaGcmfaces;
110	gforget	1.4	fldTile=fldTiles{ff};
111	gforget	1.1	%reverse order of dimensions
112			nn=length(size(fldTile));
113			fldTile=permute(fldTile,[nn:-1:1]);
114	gforget	1.7	if ntile==1; clear fldTiles; end;
115	gforget	1.1	else;
116			fldTile=fldIn;
117			end;
118	gforget	1.4	fileTile=[fileOut sprintf('.%04d.nc',ff)];
119	gforget	1.1
120			%select dimensions of relevance:
121	gforget	1.8	if ~fldInIsEmpty;
122			nDim=length(size(fldTile));
123			dimsize=size(fldTile);
124			elseif ~isempty(dimsize);
125			nDim=length(dimsize);
126			else;
127			error('undertermined array size');
128			end;
129
130	gforget	1.1	for iDim=1:nDim;
131	gforget	1.8	if dimsize(iDim)~=1;
132	gforget	1.1	dimlist{iDim}=['i' num2str(iDim)];
133			dimName{iDim}=['array index ' num2str(iDim)];
134	gforget	1.8	eval(['dimvec.i' num2str(iDim) '=[1:dimsize(iDim)];']);
135	gforget	1.1	end;
136			end;
137
138			%omit singleton dimensions:
139			ii=find(dimsize~=1);
140			dimsize=dimsize(ii);
141			dimlist={dimlist{ii}};
142			dimName={dimName{ii}};
143
144			%check :
145			if doCheck;
146			whos fldTile
147			descr
148			fldName
149			longName
150			units
151			missval
152			fillval
153			dimlist
154			dimName
155			dimsize
156			dimvec
157			keyboard;
158			end;
159
160			%output dimension information
161			dimOut{ff}=dimlist;
162
163			if doCreate;
164			%create netcdf file:
165			%-------------------
166	gforget	1.8	if prod(dimsize)*4/1e9<1.5;%use (always available) basic netcdf:
167	gforget	1.7	mode='clobber';
168			else;%to allow for large file:
169			mode='NETCDF4';
170			end;
171			ncid=nccreate(fileTile,mode);
172	gforget	1.6	nc_global=netcdf.getConstant('NC_GLOBAL');
173	gforget	1.1
174	gforget	1.3	if ~isempty(rdm);
175			descr2=[descr ' -- ' rdm{1}];
176			else;
177			descr2=descr;
178			end;
179	gforget	1.4	ncputAtt(ncid,'','description',descr2);
180	gforget	1.2	for pp=2:length(rdm);
181			tmp1=char(pp+63);
182			netcdf.putAtt(ncid,nc_global,tmp1,rdm{pp});
183			end;
184			%append readme
185			pp=length(rdm)+1;
186			netcdf.putAtt(ncid,nc_global,tmp1,'file created using gcmfaces_IO/write2nctiles.m');
187	gforget	1.1	ncputAtt(ncid,'','date',date);
188	gforget	1.5	netcdf.putAtt(ncid,nc_global,'Conventions','CF-1.6')
189
190			ncputAtt(ncid,'','_FillValue',fillval);
191			ncputAtt(ncid,'','missing_value',missval);
192	gforget	1.1
193			ncdefDim(ncid,'itxt',30);
194			for dd=1:length(dimlist); ncdefDim(ncid,dimlist{dd},dimsize(dd)); end;
195
196			for dd=1:length(dimlist);
197			ncdefVar(ncid,dimlist{dd},'double',{dimlist{dd}});
198			ncputAtt(ncid,dimlist{dd},'long_name',dimName{dd});
199	gforget	1.5	ncputAtt(ncid,dimlist{dd},'units','1');
200	gforget	1.1	end;
201			ncclose(ncid);
202
203			%fill in the dimensions dimensions values vectors:
204			%-------------------------------------------------
205			ncid=ncopen(fileTile,'write');
206			for dd=1:length(dimlist);
207			ncputvar(ncid,dimlist{dd},getfield(dimvec,dimlist{dd}));
208			end;
209			ncclose(ncid);
210			end;
211
212			%use dimentsion specified by user
213			if ~isempty(dimIn); dimlist=dimIn{ff}; end;
214
215			%define and fill field:
216			%----------------------
217			ncid=ncopen(fileTile,'write');
218	gforget	1.5	%
219			netcdf.reDef(ncid);
220	gforget	1.8	ncdefVar(ncid,fldName,xtype,flipdim(dimlist,2));%note the direction flip
221	gforget	1.1	if ~isempty(longName); ncputAtt(ncid,fldName,'long_name',longName); end;
222			if ~isempty(units); ncputAtt(ncid,fldName,'units',units); end;
223	gforget	1.5	if ~isempty(coord); ncputAtt(ncid,fldName,'coordinates',coord); end;
224			if strcmp(fldName,'lon'); ncputAtt(ncid,fldName,'standard_name','longitude'); end;
225			if strcmp(fldName,'lat'); ncputAtt(ncid,fldName,'standard_name','latitude'); end;
226			if strcmp(fldName,'dep'); ncputAtt(ncid,fldName,'standard_name','depth'); end;
227			if strcmp(fldName,'tim'); ncputAtt(ncid,fldName,'standard_name','time'); end;
228			if strcmp(fldName,'land'); ncputAtt(ncid,fldName,'standard_name','land_binary_mask'); end;
229			if strcmp(fldName,'area'); ncputAtt(ncid,fldName,'standard_name','cell_area'); end;
230			if strcmp(fldName,'thic'); ncputAtt(ncid,fldName,'standard_name','cell_thickness'); end;
231			netcdf.endDef(ncid);
232			%
233	gforget	1.8	if ~fldInIsEmpty; ncputvar(ncid,fldName,fldTile); end;
234			%
235	gforget	1.1	ncclose(ncid);
236
237	gforget	1.4	end;%for ff=1:ntile;
238	gforget	1.1