matlab_class/gcmfaces_IO/process2nctiles.m

function []=process2nctiles(dirModel,fileModel,fldModel);
%process2nctiles(dirModel);
%object : convert MITgcm binary output to netcdf files (tiled) 
%inputs : dirModel is the MITgcm run directory
%           It is expected to contain binaries in 
%           'diags/STATE/', 'diags/TRSP/', etc. as well
%           as the 'available_diagnostics.log' text file.
%         fileModel the file name base e.g. 'state_2d_set1'
%           By default : all variables in e.g. 'state_2d_set1*' 
%           files will be processed, and writen individually to
%           nctiles (tiled netcdf) that will be located in 'nctiles/'
%         fldModel (optional) can be specifiec as e.g. 'ETAN'
%output : (netcdf files)

gcmfaces_global;

fprintf('to do list : embed README in global attributes\n');
fprintf('to do list : create tim variable\n');

%directory names
listDirs={'STATE/','TRSP/'};
filAvailDiag=[dirModel 'available_diagnostics.log'];
dirOut=[dirModel 'nctiles/'];
if ~isdir(dirOut); mkdir(dirOut); end;

%seacrh in subdirectories
subDir=[];
for ff=1:length(listDirs);
tmp1=dir([dirModel 'diags/' listDirs{ff} fileModel '*']);
if ~isempty(tmp1); subDir=listDirs{ff}; end;
end;

if isempty(subDir); 
  error(['file ' fileModel ' was not found']);
else;
  dirIn=[dirModel 'diags/' subDir];
  nn=length(dir([dirIn fileModel '*data']));
  fprintf('%s (%d files) was found in \n %s \n',fileModel,nn,dirIn);
end;

%set list of variables to process
if ~isempty(who('fldModel'));
   if ischar(fldModel); listFlds={fldModel};
   else; listFlds=fldModel;
   end;
else;
   meta=read_meta([dirIn fileModel '*']);
   listFlds=meta.fldList;
end;

%now do the actual processing
for vv=1:length(listFlds);
nameDiag=deblank(listFlds{vv}) 

%get meta information
meta=read_meta([dirIn fileModel '*']);
irec=find(strcmp(deblank(meta.fldList),nameDiag));
if length(irec)~=1; error('field not in file\n'); end;

%read time series
myDiag=rdmds2gcmfaces([dirIn fileModel '*'],NaN,'rec',irec);

%ancilliary fields for netcdf file
tim=[1:size(myDiag{1},3)];
timName='month index';
descr=nameDiag;

%get units and long name from available_diagnostics.log
[avail_diag]=read_avail_diag(filAvailDiag,nameDiag);

%convert to MITgcm format (90x1170 array)
myFile=[dirOut nameDiag];

%get grid params
[grid_diag]=set_grid_diag(avail_diag);

%create netcdf file using write2nctiles (works only with old matlab, thus far ...)
doCreate=1; 
dimOut=write2nctiles(myFile,myDiag,doCreate,...
    {'fldName',nameDiag},{'longName',avail_diag.longNameDiag},...
    {'units',avail_diag.units},{'descr',descr});

%prepare to add fields
doCreate=0;

%determine relevant dimensions (note the reverse order)
for ff=1:mygrid.nFaces;
  dim2d{ff}={dimOut{ff}{end-1:end}};
  if avail_diag.nr~=1; 
    dimmsk{ff}={dimOut{ff}{end-2:end}};
    dimdep{ff}={dimOut{ff}{end-2}};
  else;
    dimmsk{ff}=dim2d{ff};
    dimdep{ff}=[];
  end;
end;

%now add fields
write2nctiles(myFile,grid_diag.lon,doCreate,...
  {'fldName','longitude'},{'units',''},{'dimIn',dim2d});
write2nctiles(myFile,grid_diag.lat,doCreate,...
  {'fldName','latitude'},{'units',''},{'dimIn',dim2d});
write2nctiles(myFile,grid_diag.msk,doCreate,...
  {'fldName','landmask'},{'units',''},...
  {'longName','land mask'},{'dimIn',dimmsk});
write2nctiles(myFile,grid_diag.RAC,doCreate,...
  {'fldName','cellarea'},{'units','m^2'},...
  {'longName','grid cell area'},{'dimIn',dim2d});
if isfield(grid_diag,'dep');
    write2nctiles(myFile,grid_diag.dep,doCreate,...
      {'fldName','depth'},{'units','m'},{'dimIn',dimdep});
    write2nctiles(myFile,grid_diag.dz,doCreate,...
      {'fldName','cellthick'},{'units','m'},{'dimIn',dimdep});
end;

end;%for vv=1:length(listFlds);

function [meta]=read_meta(fileName);

%read meta file
tmp1=dir([fileName '*.meta']); tmp1=tmp1(1).name;
tmp2=strfind(fileName,'/');
if ~isempty(tmp2); tmp2=tmp2(end); else; tmp2=0; end;
tmp1=[fileName(1:tmp2) tmp1]; fid=fopen(tmp1);
while 1;
    tline = fgetl(fid);
    if ~ischar(tline), break, end
    if isempty(whos('tmp3')); tmp3=tline; else; tmp3=[tmp3 ' ' tline]; end;
end
fclose(fid);

%add meta variables to workspace
eval(tmp3);

%reformat to meta structure
meta.dataprec=dataprec;
meta.nDims=nDims;
meta.nFlds=nFlds;
meta.nrecords=nrecords;
meta.fldList=fldList;
meta.dimList=dimList;
if ~isempty(who('timeInterval')); meta.timeInterval=timeInterval; end;
if ~isempty(who('timeStepNumber'));  meta.timeStepNumber=timeStepNumber; end;

%%

function [avail_diag]=read_avail_diag(filAvailDiag,nameDiag);

gcmfaces_global;

avail_diag=[];

fid=fopen(filAvailDiag,'rt');
while ~feof(fid);
    tline = fgetl(fid);
    tmp1=8-length(nameDiag); tmp1=repmat(' ',[1 tmp1]);
    tname = ['|' sprintf('%s',nameDiag) tmp1 '|'];
    if ~isempty(strfind(tline,tname));
        %e.g. tline='   235 |SIatmQnt|  1 |       |SM      U1|W/m^2           |Net atmospheric heat flux, >0 decreases theta';
        %
        tmp1=strfind(tline,'|'); tmp1=tmp1(end-1:end);
        avail_diag.units=strtrim(tline(tmp1(1)+1:tmp1(2)-1));
        avail_diag.longNameDiag=tline(tmp1(2)+1:end);
        %
        tmp1=strfind(tline,'|'); tmp1=tmp1(4:5);
        pars=tline(tmp1(1)+1:tmp1(2)-1);
        %
        if strcmp(pars(2),'M'); avail_diag.loc_h='C';
        elseif strcmp(pars(2),'U'); avail_diag.loc_h='W';
        elseif strcmp(pars(2),'V'); avail_diag.loc_h='S';
        end;
        %
        avail_diag.loc_z=pars(9);
        %
        if strcmp(pars(10),'1'); avail_diag.nr=1;
        else; avail_diag.nr=length(mygrid.RC);
        end;
    end;
end;
fclose(fid);

%%

function [grid_diag]=set_grid_diag(avail_diag);

gcmfaces_global;

%switch for non-tracer point values
if strcmp(avail_diag.loc_h,'C');
    grid_diag.lon=mygrid.XC;
    grid_diag.lat=mygrid.YC;
    grid_diag.msk=mygrid.mskC(:,:,1:avail_diag.nr);
elseif strcmp(avail_diag.loc_h,'W');
    grid_diag.lon=mygrid.XG;
    grid_diag.lat=mygrid.YC;
    grid_diag.msk=mygrid.mskW(:,:,1:avail_diag.nr);
elseif strcmp(avail_diag.loc_h,'S');
    grid_diag.lon=mygrid.XC;
    grid_diag.lat=mygrid.YG;
    grid_diag.msk=mygrid.mskS(:,:,1:avail_diag.nr);
end;
grid_diag.RAC=mygrid.RAC;

%vertical grid
if avail_diag.nr~=1;
    if strcmp(avail_diag.loc_z,'M');
        grid_diag.dep=-mygrid.RC;
        grid_diag.dz=mygrid.DRF;
    elseif strcmp(avail_diag.loc_z,'L');
        grid_diag.dep=-mygrid.RF(2:end);
        grid_diag.dz=[mygrid.DRC(2:end) ; 228.25];%quick fix
    else;
        error('unknown vertical grid');
    end;
    grid_diag.dep=reshape(grid_diag.dep,[1 1 avail_diag.nr]);
    grid_diag.dz=reshape(grid_diag.dz,[1 1 avail_diag.nr]);
end;
1	gforget	1.1	function []=process2nctiles(dirModel,fileModel,fldModel);
2			%process2nctiles(dirModel);
3			%object : convert MITgcm binary output to netcdf files (tiled)
4			%inputs : dirModel is the MITgcm run directory
5			% It is expected to contain binaries in
6			% 'diags/STATE/', 'diags/TRSP/', etc. as well
7			% as the 'available_diagnostics.log' text file.
8			% fileModel the file name base e.g. 'state_2d_set1'
9			% By default : all variables in e.g. 'state_2d_set1*'
10			% files will be processed, and writen individually to
11			% nctiles (tiled netcdf) that will be located in 'nctiles/'
12			% fldModel (optional) can be specifiec as e.g. 'ETAN'
13			%output : (netcdf files)
14
15			gcmfaces_global;
16
17			fprintf('to do list : embed README in global attributes\n');
18			fprintf('to do list : create tim variable\n');
19
20			%directory names
21			listDirs={'STATE/','TRSP/'};
22			filAvailDiag=[dirModel 'available_diagnostics.log'];
23			dirOut=[dirModel 'nctiles/'];
24			if ~isdir(dirOut); mkdir(dirOut); end;
25
26			%seacrh in subdirectories
27			subDir=[];
28			for ff=1:length(listDirs);
29			tmp1=dir([dirModel 'diags/' listDirs{ff} fileModel '*']);
30			if ~isempty(tmp1); subDir=listDirs{ff}; end;
31			end;
32
33			if isempty(subDir);
34			error(['file ' fileModel ' was not found']);
35			else;
36			dirIn=[dirModel 'diags/' subDir];
37			nn=length(dir([dirIn fileModel '*data']));
38			fprintf('%s (%d files) was found in \n %s \n',fileModel,nn,dirIn);
39			end;
40
41			%set list of variables to process
42			if ~isempty(who('fldModel'));
43			if ischar(fldModel); listFlds={fldModel};
44			else; listFlds=fldModel;
45			end;
46			else;
47			meta=read_meta([dirIn fileModel '*']);
48			listFlds=meta.fldList;
49			end;
50
51			%now do the actual processing
52			for vv=1:length(listFlds);
53			nameDiag=deblank(listFlds{vv})
54
55			%get meta information
56			meta=read_meta([dirIn fileModel '*']);
57			irec=find(strcmp(deblank(meta.fldList),nameDiag));
58			if length(irec)~=1; error('field not in file\n'); end;
59
60			%read time series
61			myDiag=rdmds2gcmfaces([dirIn fileModel '*'],NaN,'rec',irec);
62
63			%ancilliary fields for netcdf file
64			tim=[1:size(myDiag{1},3)];
65			timName='month index';
66			descr=nameDiag;
67
68			%get units and long name from available_diagnostics.log
69			[avail_diag]=read_avail_diag(filAvailDiag,nameDiag);
70
71			%convert to MITgcm format (90x1170 array)
72			myFile=[dirOut nameDiag];
73
74			%get grid params
75			[grid_diag]=set_grid_diag(avail_diag);
76
77			%create netcdf file using write2nctiles (works only with old matlab, thus far ...)
78			doCreate=1;
79			dimOut=write2nctiles(myFile,myDiag,doCreate,...
80			{'fldName',nameDiag},{'longName',avail_diag.longNameDiag},...
81			{'units',avail_diag.units},{'descr',descr});
82
83			%prepare to add fields
84			doCreate=0;
85
86			%determine relevant dimensions (note the reverse order)
87			for ff=1:mygrid.nFaces;
88			dim2d{ff}={dimOut{ff}{end-1:end}};
89			if avail_diag.nr~=1;
90			dimmsk{ff}={dimOut{ff}{end-2:end}};
91			dimdep{ff}={dimOut{ff}{end-2}};
92			else;
93			dimmsk{ff}=dim2d{ff};
94			dimdep{ff}=[];
95			end;
96			end;
97
98			%now add fields
99			write2nctiles(myFile,grid_diag.lon,doCreate,...
100			{'fldName','longitude'},{'units',''},{'dimIn',dim2d});
101			write2nctiles(myFile,grid_diag.lat,doCreate,...
102			{'fldName','latitude'},{'units',''},{'dimIn',dim2d});
103			write2nctiles(myFile,grid_diag.msk,doCreate,...
104			{'fldName','landmask'},{'units',''},...
105			{'longName','land mask'},{'dimIn',dimmsk});
106			write2nctiles(myFile,grid_diag.RAC,doCreate,...
107			{'fldName','cellarea'},{'units','m^2'},...
108			{'longName','grid cell area'},{'dimIn',dim2d});
109			if isfield(grid_diag,'dep');
110			write2nctiles(myFile,grid_diag.dep,doCreate,...
111			{'fldName','depth'},{'units','m'},{'dimIn',dimdep});
112			write2nctiles(myFile,grid_diag.dz,doCreate,...
113			{'fldName','cellthick'},{'units','m'},{'dimIn',dimdep});
114			end;
115
116			end;%for vv=1:length(listFlds);
117
118			function [meta]=read_meta(fileName);
119
120			%read meta file
121			tmp1=dir([fileName '*.meta']); tmp1=tmp1(1).name;
122			tmp2=strfind(fileName,'/');
123			if ~isempty(tmp2); tmp2=tmp2(end); else; tmp2=0; end;
124			tmp1=[fileName(1:tmp2) tmp1]; fid=fopen(tmp1);
125			while 1;
126			tline = fgetl(fid);
127			if ~ischar(tline), break, end
128			if isempty(whos('tmp3')); tmp3=tline; else; tmp3=[tmp3 ' ' tline]; end;
129			end
130			fclose(fid);
131
132			%add meta variables to workspace
133			eval(tmp3);
134
135			%reformat to meta structure
136			meta.dataprec=dataprec;
137			meta.nDims=nDims;
138			meta.nFlds=nFlds;
139			meta.nrecords=nrecords;
140			meta.fldList=fldList;
141			meta.dimList=dimList;
142			if ~isempty(who('timeInterval')); meta.timeInterval=timeInterval; end;
143			if ~isempty(who('timeStepNumber')); meta.timeStepNumber=timeStepNumber; end;
144
145			%%
146
147			function [avail_diag]=read_avail_diag(filAvailDiag,nameDiag);
148
149			gcmfaces_global;
150
151			avail_diag=[];
152
153			fid=fopen(filAvailDiag,'rt');
154			while ~feof(fid);
155			tline = fgetl(fid);
156			tmp1=8-length(nameDiag); tmp1=repmat(' ',[1 tmp1]);
157			tname = ['\|' sprintf('%s',nameDiag) tmp1 '\|'];
158			if ~isempty(strfind(tline,tname));
159			%e.g. tline=' 235 \|SIatmQnt\| 1 \| \|SM U1\|W/m^2 \|Net atmospheric heat flux, >0 decreases theta';
160			%
161			tmp1=strfind(tline,'\|'); tmp1=tmp1(end-1:end);
162			avail_diag.units=strtrim(tline(tmp1(1)+1:tmp1(2)-1));
163			avail_diag.longNameDiag=tline(tmp1(2)+1:end);
164			%
165			tmp1=strfind(tline,'\|'); tmp1=tmp1(4:5);
166			pars=tline(tmp1(1)+1:tmp1(2)-1);
167			%
168			if strcmp(pars(2),'M'); avail_diag.loc_h='C';
169			elseif strcmp(pars(2),'U'); avail_diag.loc_h='W';
170			elseif strcmp(pars(2),'V'); avail_diag.loc_h='S';
171			end;
172			%
173			avail_diag.loc_z=pars(9);
174			%
175			if strcmp(pars(10),'1'); avail_diag.nr=1;
176			else; avail_diag.nr=length(mygrid.RC);
177			end;
178			end;
179			end;
180			fclose(fid);
181
182			%%
183
184			function [grid_diag]=set_grid_diag(avail_diag);
185
186			gcmfaces_global;
187
188			%switch for non-tracer point values
189			if strcmp(avail_diag.loc_h,'C');
190			grid_diag.lon=mygrid.XC;
191			grid_diag.lat=mygrid.YC;
192			grid_diag.msk=mygrid.mskC(:,:,1:avail_diag.nr);
193			elseif strcmp(avail_diag.loc_h,'W');
194			grid_diag.lon=mygrid.XG;
195			grid_diag.lat=mygrid.YC;
196			grid_diag.msk=mygrid.mskW(:,:,1:avail_diag.nr);
197			elseif strcmp(avail_diag.loc_h,'S');
198			grid_diag.lon=mygrid.XC;
199			grid_diag.lat=mygrid.YG;
200			grid_diag.msk=mygrid.mskS(:,:,1:avail_diag.nr);
201			end;
202			grid_diag.RAC=mygrid.RAC;
203
204			%vertical grid
205			if avail_diag.nr~=1;
206			if strcmp(avail_diag.loc_z,'M');
207			grid_diag.dep=-mygrid.RC;
208			grid_diag.dz=mygrid.DRF;
209			elseif strcmp(avail_diag.loc_z,'L');
210			grid_diag.dep=-mygrid.RF(2:end);
211			grid_diag.dz=[mygrid.DRC(2:end) ; 228.25];%quick fix
212			else;
213			error('unknown vertical grid');
214			end;
215			grid_diag.dep=reshape(grid_diag.dep,[1 1 avail_diag.nr]);
216			grid_diag.dz=reshape(grid_diag.dz,[1 1 avail_diag.nr]);
217			end;