matlab_class/gcmfaces_IO/process2nctiles.m

function []=process2nctiles(dirModel,fileModel,fldModel,tileSize);
%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 (by default []) can be specified (as e.g. 'ETAN')
%            when fldModel is empty, all fields are processed
%         tileSize (optional) is e.g. [90 90] (by default tiles=faces) 
%output : (netcdf files)

gcmfaces_global;

%listFiles={'state_2d_set1','state_2d_set2','state_3d_set1','state_3d_set2'};
%listFiles={'trsp_3d_set1','trsp_3d_set2','trsp_3d_set3'};
%for ff=1:length(listFiles); process2nctiles('iter12/',listFiles{ff},[],[90 90]); end;

%replace time series with monthly climatology?
doClim=0;

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

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

if isempty(subDir);
tmp1=dir([dirModel diagsDir fileModel '/' fileModel '*']);
if ~isempty(tmp1); subDir=[fileModel '/']; end;
end;

if isempty(subDir); 
  error(['file ' fileModel ' was not found']);
else;
  dirIn=[dirModel diagsDir 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(fldModel);
   if ischar(fldModel); listFlds={fldModel};
   else; listFlds=fldModel;
   end;
else;
   meta=read_meta([dirIn fileModel '*']);
   listFlds=meta.fldList;
end;

%determine map of tile indices (by default tiles=faces)
if isempty(whos('tileSize'));
  tileNo=mygrid.XC; 
  for ff=1:mygrid.nFaces; tileNo{ff}(:)=ff; end;
else;
  tileNo=gcmfaces_loc_tile(tileSize(1),tileSize(2));
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);

%replace time series with monthly climatology
if doClim; myDiag=compClim(myDiag); end;

%set ancilliary time variable
nn=length(size(myDiag{1}));
nn=size(myDiag{1},nn);
%tim=[1:nn];
tim=[1992*ones(nn,1) [1:nn]' 15*ones(nn,1)];
tim=datenum(tim)-datenum([1992 1 0]);
timUnits='days since 1992-1-1 0:0:0';

%get time step axis
[listTimes]=diags_list_times({dirIn},{fileModel});

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

%get description of estimate from README
[rdm]=read_readme(filReadme);
disp(rdm');

%set output directory/file name
myFile=[dirOut nameDiag];%first instance is for subdirectory name
if ~isdir(myFile); mkdir(myFile); end;
myFile=[myFile filesep nameDiag];%second instance is for file name base

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

%apply mask, and convert to land mask
if ~isempty(mygrid.RAC);
  msk=grid_diag.msk;
  if length(size(myDiag{1}))==3;
    msk=repmat(msk(:,:,1),[1 1 size(myDiag{1},3)]);
  else;
    msk=repmat(msk,[1 1 1 size(myDiag{1},4)]);
  end;
  myDiag=myDiag.*msk;
  clear msk;
  %
  land=isnan(grid_diag.msk);
end;

%set 'coord' attribute
if avail_diag.nr~=1;
  coord='lon lat dep tim';
else;
  coord='lon lat tim';
end;

%replace time series with monthly climatology
if doClim;
  listTimes=listTimes(1:12);
  timUnits='days since year-1-1 0:0:0';
  avail_diag.longNameDiag=[avail_diag.longNameDiag ' (climatology) '];
end;

%create netcdf file using write2nctiles
doCreate=1; 
dimlist=write2nctiles(myFile,myDiag,doCreate,{'tileNo',tileNo},...
    {'fldName',nameDiag},{'longName',avail_diag.longNameDiag},...
    {'units',avail_diag.units},{'descr',nameDiag},{'coord',coord},{'rdm',rdm});

%determine relevant dimensions
for ff=1:length(dimlist);
  dim.tim{ff}={dimlist{ff}{1}};
  dim.twoD{ff}={dimlist{ff}{end-1:end}};
  if avail_diag.nr~=1;
    dim.threeD{ff}={dimlist{ff}{end-2:end}};
    dim.dep{ff}={dimlist{ff}{end-2}};
  else;
    dim.threeD{ff}=dim.twoD{ff};
    dim.dep{ff}=[];
  end;
end;

%prepare to add fields
doCreate=0;

%now add fields
write2nctiles(myFile,grid_diag.lon,doCreate,{'tileNo',tileNo},...
  {'fldName','lon'},{'units','degrees_east'},{'dimIn',dim.twoD});
write2nctiles(myFile,grid_diag.lat,doCreate,{'tileNo',tileNo},...
  {'fldName','lat'},{'units','degrees_north'},{'dimIn',dim.twoD});
if isfield(grid_diag,'dep');
    write2nctiles(myFile,grid_diag.dep,doCreate,{'tileNo',tileNo},...
      {'fldName','dep'},{'units','m'},{'dimIn',dim.dep});
end;
write2nctiles(myFile,tim,doCreate,{'tileNo',tileNo},...
  {'fldName','tim'},{'longName','time'},...
  {'units',timUnits},{'dimIn',dim.tim});
if ~isempty(mygrid.RAC);
  write2nctiles(myFile,listTimes,doCreate,{'tileNo',tileNo},...
    {'fldName','timstep'},{'longName','final time step number'},...
    {'units','1'},{'dimIn',dim.tim});
  write2nctiles(myFile,grid_diag.msk,doCreate,{'tileNo',tileNo},...
    {'fldName','land'},{'units','1'},{'longName','land mask'},{'dimIn',dim.threeD});
  write2nctiles(myFile,grid_diag.RAC,doCreate,{'tileNo',tileNo},...
    {'fldName','area'},{'units','m^2'},{'longName','grid cell area'},{'dimIn',dim.twoD});
  if isfield(grid_diag,'dz');
    write2nctiles(myFile,grid_diag.dz,doCreate,{'tileNo',tileNo},...
      {'fldName','thic'},{'units','m'},{'dimIn',dim.dep});
  end;
end;

clear myDiag;

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,filesep);
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 [rdm]=read_readme(filReadme);

gcmfaces_global;

rdm=[];

fid=fopen(filReadme,'rt');
while ~feof(fid);
    nn=length(rdm);
    rdm{nn+1} = fgetl(fid);
end;
fclose(fid);

%%

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;

%%replace time series with monthly climatology
function [FLD]=compClim(fld);

gcmfaces_global;

ndim=length(size(fld{1}));
nyear=size(fld{1},ndim)/12;

if ndim==3; FLD=NaN*fld(:,:,1:12); end;
if ndim==4; FLD=NaN*fld(:,:,:,1:12); end;

for mm=1:12;
if ndim==3; FLD(:,:,mm)=mean(fld(:,:,mm:12:12*nyear),ndim); end;
if ndim==4; FLD(:,:,:,mm)=mean(fld(:,:,:,mm:12:12*nyear),ndim); end;
end;

1	gforget	1.4	function []=process2nctiles(dirModel,fileModel,fldModel,tileSize);
2	gforget	1.1	%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	gforget	1.4	% fldModel (by default []) can be specified (as e.g. 'ETAN')
13			% when fldModel is empty, all fields are processed
14			% tileSize (optional) is e.g. [90 90] (by default tiles=faces)
15	gforget	1.1	%output : (netcdf files)
16
17			gcmfaces_global;
18
19	gforget	1.5	%listFiles={'state_2d_set1','state_2d_set2','state_3d_set1','state_3d_set2'};
20			%listFiles={'trsp_3d_set1','trsp_3d_set2','trsp_3d_set3'};
21			%for ff=1:length(listFiles); process2nctiles('iter12/',listFiles{ff},[],[90 90]); end;
22
23	gforget	1.9	%replace time series with monthly climatology?
24			doClim=0;
25
26	gforget	1.1	%directory names
27	gforget	1.10	listDirs={'STATE/','TRSP/'};%BUDG?
28	gforget	1.1	filAvailDiag=[dirModel 'available_diagnostics.log'];
29	gforget	1.2	filReadme=[dirModel 'README'];
30	gforget	1.9	dirOut=[dirModel 'nctiles_tmp/'];
31	gforget	1.13	%dirOut=[dirModel 'nctiles_post_tmp/'];
32	gforget	1.1	if ~isdir(dirOut); mkdir(dirOut); end;
33
34	gforget	1.9	%search in subdirectories
35	gforget	1.1	subDir=[];
36	gforget	1.13	diagsDir='diags/';
37			%diagsDir='diags_post/';
38			%diagsDir='diags_interp/';
39	gforget	1.1	for ff=1:length(listDirs);
40	gforget	1.13	tmp1=dir([dirModel diagsDir listDirs{ff} fileModel '*']);
41	gforget	1.1	if ~isempty(tmp1); subDir=listDirs{ff}; end;
42			end;
43
44	gforget	1.13	if isempty(subDir);
45			tmp1=dir([dirModel diagsDir fileModel '/' fileModel '*']);
46			if ~isempty(tmp1); subDir=[fileModel '/']; end;
47			end;
48
49	gforget	1.1	if isempty(subDir);
50			error(['file ' fileModel ' was not found']);
51			else;
52	gforget	1.13	dirIn=[dirModel diagsDir subDir];
53	gforget	1.1	nn=length(dir([dirIn fileModel '*data']));
54			fprintf('%s (%d files) was found in \n %s \n',fileModel,nn,dirIn);
55			end;
56
57			%set list of variables to process
58	gforget	1.4	if ~isempty(fldModel);
59	gforget	1.1	if ischar(fldModel); listFlds={fldModel};
60			else; listFlds=fldModel;
61			end;
62			else;
63			meta=read_meta([dirIn fileModel '*']);
64			listFlds=meta.fldList;
65			end;
66
67	gforget	1.4	%determine map of tile indices (by default tiles=faces)
68			if isempty(whos('tileSize'));
69			tileNo=mygrid.XC;
70			for ff=1:mygrid.nFaces; tileNo{ff}(:)=ff; end;
71			else;
72			tileNo=gcmfaces_loc_tile(tileSize(1),tileSize(2));
73			end;
74
75	gforget	1.1	%now do the actual processing
76			for vv=1:length(listFlds);
77			nameDiag=deblank(listFlds{vv})
78
79			%get meta information
80			meta=read_meta([dirIn fileModel '*']);
81			irec=find(strcmp(deblank(meta.fldList),nameDiag));
82			if length(irec)~=1; error('field not in file\n'); end;
83
84			%read time series
85			myDiag=rdmds2gcmfaces([dirIn fileModel '*'],NaN,'rec',irec);
86
87	gforget	1.9	%replace time series with monthly climatology
88			if doClim; myDiag=compClim(myDiag); end;
89
90	gforget	1.5	%set ancilliary time variable
91	gforget	1.3	nn=length(size(myDiag{1}));
92	gforget	1.5	nn=size(myDiag{1},nn);
93			%tim=[1:nn];
94			tim=[1992ones(nn,1) [1:nn]' 15ones(nn,1)];
95			tim=datenum(tim)-datenum([1992 1 0]);
96			timUnits='days since 1992-1-1 0:0:0';
97	gforget	1.1
98	gforget	1.7	%get time step axis
99	gforget	1.10	[listTimes]=diags_list_times({dirIn},{fileModel});
100	gforget	1.7
101	gforget	1.1	%get units and long name from available_diagnostics.log
102			[avail_diag]=read_avail_diag(filAvailDiag,nameDiag);
103
104	gforget	1.2	%get description of estimate from README
105			[rdm]=read_readme(filReadme);
106			disp(rdm');
107
108	gforget	1.12	%set output directory/file name
109	gforget	1.4	myFile=[dirOut nameDiag];%first instance is for subdirectory name
110			if ~isdir(myFile); mkdir(myFile); end;
111	gforget	1.12	myFile=[myFile filesep nameDiag];%second instance is for file name base
112	gforget	1.1
113			%get grid params
114			[grid_diag]=set_grid_diag(avail_diag);
115
116	gforget	1.5	%apply mask, and convert to land mask
117	gforget	1.13	if ~isempty(mygrid.RAC);
118			msk=grid_diag.msk;
119			if length(size(myDiag{1}))==3;
120			msk=repmat(msk(:,:,1),[1 1 size(myDiag{1},3)]);
121			else;
122			msk=repmat(msk,[1 1 1 size(myDiag{1},4)]);
123			end;
124			myDiag=myDiag.*msk;
125			clear msk;
126			%
127			land=isnan(grid_diag.msk);
128			end;
129	gforget	1.5
130			%set 'coord' attribute
131			if avail_diag.nr~=1;
132			coord='lon lat dep tim';
133			else;
134			coord='lon lat tim';
135			end;
136
137	gforget	1.9	%replace time series with monthly climatology
138			if doClim;
139			listTimes=listTimes(1:12);
140			timUnits='days since year-1-1 0:0:0';
141			avail_diag.longNameDiag=[avail_diag.longNameDiag ' (climatology) '];
142			end;
143
144	gforget	1.11	%create netcdf file using write2nctiles
145	gforget	1.1	doCreate=1;
146	gforget	1.4	dimlist=write2nctiles(myFile,myDiag,doCreate,{'tileNo',tileNo},...
147	gforget	1.1	{'fldName',nameDiag},{'longName',avail_diag.longNameDiag},...
148	gforget	1.5	{'units',avail_diag.units},{'descr',nameDiag},{'coord',coord},{'rdm',rdm});
149	gforget	1.1
150	gforget	1.4	%determine relevant dimensions
151			for ff=1:length(dimlist);
152			dim.tim{ff}={dimlist{ff}{1}};
153			dim.twoD{ff}={dimlist{ff}{end-1:end}};
154			if avail_diag.nr~=1;
155			dim.threeD{ff}={dimlist{ff}{end-2:end}};
156			dim.dep{ff}={dimlist{ff}{end-2}};
157	gforget	1.1	else;
158	gforget	1.4	dim.threeD{ff}=dim.twoD{ff};
159			dim.dep{ff}=[];
160	gforget	1.1	end;
161			end;
162
163	gforget	1.4	%prepare to add fields
164			doCreate=0;
165
166	gforget	1.1	%now add fields
167	gforget	1.4	write2nctiles(myFile,grid_diag.lon,doCreate,{'tileNo',tileNo},...
168	gforget	1.5	{'fldName','lon'},{'units','degrees_east'},{'dimIn',dim.twoD});
169	gforget	1.4	write2nctiles(myFile,grid_diag.lat,doCreate,{'tileNo',tileNo},...
170	gforget	1.5	{'fldName','lat'},{'units','degrees_north'},{'dimIn',dim.twoD});
171	gforget	1.1	if isfield(grid_diag,'dep');
172	gforget	1.4	write2nctiles(myFile,grid_diag.dep,doCreate,{'tileNo',tileNo},...
173	gforget	1.5	{'fldName','dep'},{'units','m'},{'dimIn',dim.dep});
174	gforget	1.13	end;
175			write2nctiles(myFile,tim,doCreate,{'tileNo',tileNo},...
176			{'fldName','tim'},{'longName','time'},...
177			{'units',timUnits},{'dimIn',dim.tim});
178			if ~isempty(mygrid.RAC);
179			write2nctiles(myFile,listTimes,doCreate,{'tileNo',tileNo},...
180			{'fldName','timstep'},{'longName','final time step number'},...
181			{'units','1'},{'dimIn',dim.tim});
182			write2nctiles(myFile,grid_diag.msk,doCreate,{'tileNo',tileNo},...
183			{'fldName','land'},{'units','1'},{'longName','land mask'},{'dimIn',dim.threeD});
184			write2nctiles(myFile,grid_diag.RAC,doCreate,{'tileNo',tileNo},...
185			{'fldName','area'},{'units','m^2'},{'longName','grid cell area'},{'dimIn',dim.twoD});
186			if isfield(grid_diag,'dz');
187	gforget	1.4	write2nctiles(myFile,grid_diag.dz,doCreate,{'tileNo',tileNo},...
188	gforget	1.5	{'fldName','thic'},{'units','m'},{'dimIn',dim.dep});
189	gforget	1.13	end;
190	gforget	1.1	end;
191
192	gforget	1.13	clear myDiag;
193
194	gforget	1.1	end;%for vv=1:length(listFlds);
195
196			function [meta]=read_meta(fileName);
197
198			%read meta file
199			tmp1=dir([fileName '*.meta']); tmp1=tmp1(1).name;
200	gforget	1.8	tmp2=strfind(fileName,filesep);
201	gforget	1.1	if ~isempty(tmp2); tmp2=tmp2(end); else; tmp2=0; end;
202			tmp1=[fileName(1:tmp2) tmp1]; fid=fopen(tmp1);
203			while 1;
204			tline = fgetl(fid);
205			if ~ischar(tline), break, end
206			if isempty(whos('tmp3')); tmp3=tline; else; tmp3=[tmp3 ' ' tline]; end;
207			end
208			fclose(fid);
209
210			%add meta variables to workspace
211			eval(tmp3);
212
213			%reformat to meta structure
214			meta.dataprec=dataprec;
215			meta.nDims=nDims;
216			meta.nFlds=nFlds;
217			meta.nrecords=nrecords;
218			meta.fldList=fldList;
219			meta.dimList=dimList;
220			if ~isempty(who('timeInterval')); meta.timeInterval=timeInterval; end;
221			if ~isempty(who('timeStepNumber')); meta.timeStepNumber=timeStepNumber; end;
222
223			%%
224
225	gforget	1.2	function [rdm]=read_readme(filReadme);
226
227			gcmfaces_global;
228
229			rdm=[];
230
231			fid=fopen(filReadme,'rt');
232			while ~feof(fid);
233			nn=length(rdm);
234			rdm{nn+1} = fgetl(fid);
235			end;
236			fclose(fid);
237
238			%%
239
240	gforget	1.1	function [avail_diag]=read_avail_diag(filAvailDiag,nameDiag);
241
242			gcmfaces_global;
243
244			avail_diag=[];
245
246			fid=fopen(filAvailDiag,'rt');
247			while ~feof(fid);
248			tline = fgetl(fid);
249			tmp1=8-length(nameDiag); tmp1=repmat(' ',[1 tmp1]);
250			tname = ['\|' sprintf('%s',nameDiag) tmp1 '\|'];
251			if ~isempty(strfind(tline,tname));
252			%e.g. tline=' 235 \|SIatmQnt\| 1 \| \|SM U1\|W/m^2 \|Net atmospheric heat flux, >0 decreases theta';
253			%
254			tmp1=strfind(tline,'\|'); tmp1=tmp1(end-1:end);
255			avail_diag.units=strtrim(tline(tmp1(1)+1:tmp1(2)-1));
256			avail_diag.longNameDiag=tline(tmp1(2)+1:end);
257			%
258			tmp1=strfind(tline,'\|'); tmp1=tmp1(4:5);
259			pars=tline(tmp1(1)+1:tmp1(2)-1);
260			%
261			if strcmp(pars(2),'M'); avail_diag.loc_h='C';
262			elseif strcmp(pars(2),'U'); avail_diag.loc_h='W';
263			elseif strcmp(pars(2),'V'); avail_diag.loc_h='S';
264			end;
265			%
266			avail_diag.loc_z=pars(9);
267			%
268			if strcmp(pars(10),'1'); avail_diag.nr=1;
269			else; avail_diag.nr=length(mygrid.RC);
270			end;
271			end;
272			end;
273			fclose(fid);
274
275			%%
276
277			function [grid_diag]=set_grid_diag(avail_diag);
278
279			gcmfaces_global;
280
281			%switch for non-tracer point values
282			if strcmp(avail_diag.loc_h,'C');
283			grid_diag.lon=mygrid.XC;
284			grid_diag.lat=mygrid.YC;
285			grid_diag.msk=mygrid.mskC(:,:,1:avail_diag.nr);
286			elseif strcmp(avail_diag.loc_h,'W');
287			grid_diag.lon=mygrid.XG;
288			grid_diag.lat=mygrid.YC;
289			grid_diag.msk=mygrid.mskW(:,:,1:avail_diag.nr);
290			elseif strcmp(avail_diag.loc_h,'S');
291			grid_diag.lon=mygrid.XC;
292			grid_diag.lat=mygrid.YG;
293			grid_diag.msk=mygrid.mskS(:,:,1:avail_diag.nr);
294			end;
295			grid_diag.RAC=mygrid.RAC;
296
297			%vertical grid
298			if avail_diag.nr~=1;
299			if strcmp(avail_diag.loc_z,'M');
300			grid_diag.dep=-mygrid.RC;
301			grid_diag.dz=mygrid.DRF;
302			elseif strcmp(avail_diag.loc_z,'L');
303			grid_diag.dep=-mygrid.RF(2:end);
304			grid_diag.dz=[mygrid.DRC(2:end) ; 228.25];%quick fix
305			else;
306			error('unknown vertical grid');
307			end;
308			grid_diag.dep=reshape(grid_diag.dep,[1 1 avail_diag.nr]);
309			grid_diag.dz=reshape(grid_diag.dz,[1 1 avail_diag.nr]);
310			end;
311	gforget	1.9
312			%%replace time series with monthly climatology
313			function [FLD]=compClim(fld);
314
315			gcmfaces_global;
316
317			ndim=length(size(fld{1}));
318			nyear=size(fld{1},ndim)/12;
319
320			if ndim==3; FLD=NaN*fld(:,:,1:12); end;
321			if ndim==4; FLD=NaN*fld(:,:,:,1:12); end;
322
323			for mm=1:12;
324			if ndim==3; FLD(:,:,mm)=mean(fld(:,:,mm:12:12*nyear),ndim); end;
325			if ndim==4; FLD(:,:,:,mm)=mean(fld(:,:,:,mm:12:12*nyear),ndim); end;
326			end;
327