natl_box/matlab/readgcm.m

function [fld,tme,nPx,nPy]=readgcm(fnam,t,Nx,Ny,Nz,Ix,Iy,Iz,prec,nPx,nPy)

% Function [fld,tme,nPx,nPy]=readgcm(fnam,t,Nx,Ny,Nz,Ix,Iy,Iz,prec,nPx,nPy)
% read and output model field from file fnam time step t
%
% INPUTS
% fnam      input path and file name
% t         indices of time steps to read (0 means read all, the default)
% Nx*Ny*Nz  grid dimension (default 360*224*46)
% Ix,Iy,Iz  subsample indices (default Ix=1:Nx, Iy=1:Ny, and Iz=1:Nz)
% prec      numeric precision (default 'real*4')
%
% INPUT/OUTPUT
% nPx, nPy  number of processes in x and y direction;
%           when not specified as INPUTS, nPx and nPy are determined based
%           on the following file naming convention:
%             U_06_04.00060_00120_10 is
%             U for nPx=6, nPy=4, hour 60 to 120 at 10-hour intervals
%           (defaults to 1 when filename contains no underscores)
%
% OUTPUTS
% fld   Nx*Ny*Nz*Nt output array
%          e.g., to contour level k, use "contourf(fld(:,:,k,1)'), colorbar"
% tme   model integration times in days
%
% The following are valid calling forms:
%    readgcm(fnam,t,Nx,Ny,Nz,Ix,Iy,Iz,prec,nPx,nPy)
%    readgcm(fnam, ...                            )
%    readgcm(fnam, ...               ,prec, ...   )
% where "..." represents 0 or more ordered arguments.
%
% SEE ALSO
% writegcm

% "prec" is set to value of second string argument if it exists
nprec=9;
if exist('t') , if ischar(t)
    prec=t;  nprec=2;
    if exist('Nx'), nPx=Nx; end
    if exist('Ny'), nPy=Ny; end
end, end
if exist('Nx'), if ischar(Nx)
    prec=Nx; nprec=3;
    if exist('Ny'), nPx=Ny; end
    if exist('Nz'), nPy=Nz; end
end, end
if exist('Ny'), if ischar(Ny)
    prec=Ny; nprec=4;
    if exist('Nz'), nPx=Nz; end
    if exist('Ix'), nPy=Ix; end
end, end
if exist('Nz'), if ischar(Nz)
    prec=Nz; nprec=5;
    if exist('Ix'), nPx=Ix; end
    if exist('Iy'), nPy=Iy; end
end, end
if exist('Ix'), if ischar(Ix)
    prec=Ix; nprec=6;
    if exist('Iy'), nPx=Iy; end
    if exist('Iz'), nPy=Iz; end
end, end
if exist('Iy'), if ischar(Iy)
    prec=Iy; nprec=7;
    if exist('Iz'), nPx=Iz; end
    if exist('prec'), nPy=prec; end
end, end
if exist('Iz'), if ischar(Iz)
    prec=Iz; nprec=8;
    if exist('prec'), nPx=prec; end
    if exist('nPx'), nPy=nPx; end
end, end

% set default argument values
if nargin < 9 & nprec==9, prec='real*4'; end
if nargin < 5 | nprec <= 5
  if isempty(findstr(fnam,'KPPhbl')) & ...
        isempty(findstr(fnam,'H')) & ...
        isempty(findstr(fnam,'SF')) & ...
        isempty(findstr(fnam,'Pbottom')) & ...
        isempty(findstr(fnam,'Pbottom')) & ...
        isempty(findstr(fnam,'BU')) & ...
        isempty(findstr(fnam,'BV'))
    Nz=46;
  else 
    Nz=1;
  end
end
if nargin < 4 | nprec <= 4, Ny=224; end
if nargin < 3 | nprec <= 3, Nx=360; end
if nargin < 8 | nprec <= 8, Iz=1:Nz; end
if nargin < 7 | nprec <= 7, Iy=1:Ny; end
if nargin < 6 | nprec <= 6, Ix=1:Nx; end
if nargin < 2 | nprec <= 2, t=0; end
if nargin < 1, error('please specify input file name'); end
if nargin < nprec+1
  % locate beginning filename character in path
  slash_loc=findstr('/',fnam);
  if isempty(slash_loc), slash_loc=0; end
  slash_loc=max(slash_loc);
  bar_loc=findstr('_',fnam((slash_loc+1):length(fnam)));
  bar_loc=min(bar_loc);
  if isempty(bar_loc)
    nPx=1; nPy=1;
  elseif length(fnam)<slash_loc+bar_loc+5
    nPx=1; nPy=1;
  else
    nPx=str2num(fnam(slash_loc+bar_loc+(1:2)));  % number of processes in X
    nPy=str2num(fnam(slash_loc+bar_loc+(4:5)));  % number of processes in Y
    if isempty(nPx) | isempty(nPy), nPx=1; nPy=1; end
  end
elseif nargin < nprec+2
  nPy=nPx;                                  % number of processes in X and Y
end

% compute number of sub-grids
sNx=Nx/nPx;                                 % number of X points in sub-grid
sNy=Ny/nPy;                                 % number of Y points in sub-grid
if (sNx-floor(sNx))~=0 | (sNy-floor(sNy))~=0
  error('Nx/nPx and Ny/nPy must be integer')
end

% compute file record length
switch prec
  case {'float32', 'real*4'}
    rlength=(sNx*sNy*Nz+1)*4;
  case {'float64', 'real*8'}
    rlength=(sNx*sNy*Nz+1)*8;
end

% compute total number of time steps
fid=fopen(fnam,'r','ieee-be');
tmp=fseek(fid,0,'eof');
tmp=ftell(fid);
tmp=tmp/(rlength*nPx*nPy);
if t==0 | tmp<max(t)
  t=1:tmp;
end

% read model file
fld=zeros(length(Ix),length(Iy),length(Iz),length(t));
tme=zeros(1,length(t));
disp(['Reading ' int2str(length(t)) ' time steps:'])

if length(Ix)==Nx & length(Iy)==Ny & length(Iz)==Nz

  for k=1:length(t)
    fprintf(1,'\b\b\b%0.0f',k)
    tmp=fseek(fid,nPx*nPy*(t(k)-1)*rlength,'bof');
    for j=1:nPy
      jx=((j-1)*sNy+1):(j*sNy);
      for i=1:nPx
        ix=((i-1)*sNx+1):(i*sNx);
        tm=fread(fid,1,prec);
        [tmp count]=fread(fid,[sNx,sNy*Nz],prec);
        tme(k)=tm;
        fld(ix,jx,:,k)=reshape(tmp,sNx,sNy,Nz);
      end
    end
  end
  
elseif length(Ix)<Nx | length(Iy)<Ny

  % If only part of the file needs to be read, determine
  % which tiles contain the desired data.
  nPx_id=zeros(Nx,Ny);
  nPy_id=zeros(Nx,Ny);
  for j=1:nPy
    jx=((j-1)*sNy+1):(j*sNy);
    for i=1:nPx
      ix=((i-1)*sNx+1):(i*sNx);
      nPx_id(ix,jx)=i;
      nPy_id(ix,jx)=j;
    end
  end
  nPx_id=nPx_id(Ix,Iy);
  nPy_id=nPy_id(Ix,Iy);
  nPxy_id=[nPx_id(:) nPy_id(:)];
  if length(nPx_id)>1
    nPxy_id=sortrows(nPxy_id);
    nPxy_id(find(sum(abs(diff(nPxy_id)'))'==0)+1,:)=[];
  end

  fld_tmp=zeros(Nx,Ny,Nz);
  for k=1:length(t)
    fprintf(1,'\b\b\b%0.0f',k)
    for l=1:size(nPxy_id,1)
      i=nPxy_id(l,1);
      j=nPxy_id(l,2);
      rnumber=nPx*nPy*(t(k)-1)+(j-1)*nPx+i-1;
      tmp=fseek(fid,rnumber*rlength,'bof');
      ix=((i-1)*sNx+1):(i*sNx);
      jx=((j-1)*sNy+1):(j*sNy);
      tm=fread(fid,1,prec);
      [tmp count]=fread(fid,[sNx,sNy*Nz],prec);
      tme(k)=tm;
      fld_tmp(ix,jx,:)=reshape(tmp,sNx,sNy,Nz);
    end
    fld(:,:,:,k)=fld_tmp(Ix,Iy,Iz);
  end
  
elseif length(Iz)<Nz

  fld_tmp=zeros(Nx,Ny,Nz);
  for k=1:length(t)
    fprintf(1,'\b\b\b%0.0f',k)
    tmp=fseek(fid,nPx*nPy*(t(k)-1)*rlength,'bof');
    for j=1:nPy
      jx=((j-1)*sNy+1):(j*sNy);
      for i=1:nPx
        ix=((i-1)*sNx+1):(i*sNx);
        tm=fread(fid,1,prec);
        [tmp count]=fread(fid,[sNx,sNy*Nz],prec);
        tme(k)=tm;
        fld_tmp(ix,jx,:)=reshape(tmp,sNx,sNy,Nz);
      end
    end
    fld(:,:,:,k)=fld_tmp(:,:,Iz);
  end

end

fid=fclose(fid);
fprintf(1,'\b\b\b',k)

1	function [fld,tme,nPx,nPy]=readgcm(fnam,t,Nx,Ny,Nz,Ix,Iy,Iz,prec,nPx,nPy)
2
3	% Function [fld,tme,nPx,nPy]=readgcm(fnam,t,Nx,Ny,Nz,Ix,Iy,Iz,prec,nPx,nPy)
4	% read and output model field from file fnam time step t
5	%
6	% INPUTS
7	% fnam input path and file name
8	% t indices of time steps to read (0 means read all, the default)
9	% NxNyNz grid dimension (default 36022446)
10	% Ix,Iy,Iz subsample indices (default Ix=1:Nx, Iy=1:Ny, and Iz=1:Nz)
11	% prec numeric precision (default 'real*4')
12	%
13	% INPUT/OUTPUT
14	% nPx, nPy number of processes in x and y direction;
15	% when not specified as INPUTS, nPx and nPy are determined based
16	% on the following file naming convention:
17	% U_06_04.00060_00120_10 is
18	% U for nPx=6, nPy=4, hour 60 to 120 at 10-hour intervals
19	% (defaults to 1 when filename contains no underscores)
20	%
21	% OUTPUTS
22	% fld NxNyNz*Nt output array
23	% e.g., to contour level k, use "contourf(fld(:,:,k,1)'), colorbar"
24	% tme model integration times in days
25	%
26	% The following are valid calling forms:
27	% readgcm(fnam,t,Nx,Ny,Nz,Ix,Iy,Iz,prec,nPx,nPy)
28	% readgcm(fnam, ... )
29	% readgcm(fnam, ... ,prec, ... )
30	% where "..." represents 0 or more ordered arguments.
31	%
32	% SEE ALSO
33	% writegcm
34
35	% "prec" is set to value of second string argument if it exists
36	nprec=9;
37	if exist('t') , if ischar(t)
38	prec=t; nprec=2;
39	if exist('Nx'), nPx=Nx; end
40	if exist('Ny'), nPy=Ny; end
41	end, end
42	if exist('Nx'), if ischar(Nx)
43	prec=Nx; nprec=3;
44	if exist('Ny'), nPx=Ny; end
45	if exist('Nz'), nPy=Nz; end
46	end, end
47	if exist('Ny'), if ischar(Ny)
48	prec=Ny; nprec=4;
49	if exist('Nz'), nPx=Nz; end
50	if exist('Ix'), nPy=Ix; end
51	end, end
52	if exist('Nz'), if ischar(Nz)
53	prec=Nz; nprec=5;
54	if exist('Ix'), nPx=Ix; end
55	if exist('Iy'), nPy=Iy; end
56	end, end
57	if exist('Ix'), if ischar(Ix)
58	prec=Ix; nprec=6;
59	if exist('Iy'), nPx=Iy; end
60	if exist('Iz'), nPy=Iz; end
61	end, end
62	if exist('Iy'), if ischar(Iy)
63	prec=Iy; nprec=7;
64	if exist('Iz'), nPx=Iz; end
65	if exist('prec'), nPy=prec; end
66	end, end
67	if exist('Iz'), if ischar(Iz)
68	prec=Iz; nprec=8;
69	if exist('prec'), nPx=prec; end
70	if exist('nPx'), nPy=nPx; end
71	end, end
72
73	% set default argument values
74	if nargin < 9 & nprec==9, prec='real*4'; end
75	if nargin < 5 \| nprec <= 5
76	if isempty(findstr(fnam,'KPPhbl')) & ...
77	isempty(findstr(fnam,'H')) & ...
78	isempty(findstr(fnam,'SF')) & ...
79	isempty(findstr(fnam,'Pbottom')) & ...
80	isempty(findstr(fnam,'Pbottom')) & ...
81	isempty(findstr(fnam,'BU')) & ...
82	isempty(findstr(fnam,'BV'))
83	Nz=46;
84	else
85	Nz=1;
86	end
87	end
88	if nargin < 4 \| nprec <= 4, Ny=224; end
89	if nargin < 3 \| nprec <= 3, Nx=360; end
90	if nargin < 8 \| nprec <= 8, Iz=1:Nz; end
91	if nargin < 7 \| nprec <= 7, Iy=1:Ny; end
92	if nargin < 6 \| nprec <= 6, Ix=1:Nx; end
93	if nargin < 2 \| nprec <= 2, t=0; end
94	if nargin < 1, error('please specify input file name'); end
95	if nargin < nprec+1
96	% locate beginning filename character in path
97	slash_loc=findstr('/',fnam);
98	if isempty(slash_loc), slash_loc=0; end
99	slash_loc=max(slash_loc);
100	bar_loc=findstr('_',fnam((slash_loc+1):length(fnam)));
101	bar_loc=min(bar_loc);
102	if isempty(bar_loc)
103	nPx=1; nPy=1;
104	elseif length(fnam)<slash_loc+bar_loc+5
105	nPx=1; nPy=1;
106	else
107	nPx=str2num(fnam(slash_loc+bar_loc+(1:2))); % number of processes in X
108	nPy=str2num(fnam(slash_loc+bar_loc+(4:5))); % number of processes in Y
109	if isempty(nPx) \| isempty(nPy), nPx=1; nPy=1; end
110	end
111	elseif nargin < nprec+2
112	nPy=nPx; % number of processes in X and Y
113	end
114
115	% compute number of sub-grids
116	sNx=Nx/nPx; % number of X points in sub-grid
117	sNy=Ny/nPy; % number of Y points in sub-grid
118	if (sNx-floor(sNx))~=0 \| (sNy-floor(sNy))~=0
119	error('Nx/nPx and Ny/nPy must be integer')
120	end
121
122	% compute file record length
123	switch prec
124	case {'float32', 'real*4'}
125	rlength=(sNxsNyNz+1)*4;
126	case {'float64', 'real*8'}
127	rlength=(sNxsNyNz+1)*8;
128	end
129
130	% compute total number of time steps
131	fid=fopen(fnam,'r','ieee-be');
132	tmp=fseek(fid,0,'eof');
133	tmp=ftell(fid);
134	tmp=tmp/(rlengthnPxnPy);
135	if t==0 \| tmp<max(t)
136	t=1:tmp;
137	end
138
139	% read model file
140	fld=zeros(length(Ix),length(Iy),length(Iz),length(t));
141	tme=zeros(1,length(t));
142	disp(['Reading ' int2str(length(t)) ' time steps:'])
143
144	if length(Ix)==Nx & length(Iy)==Ny & length(Iz)==Nz
145
146	for k=1:length(t)
147	fprintf(1,'\b\b\b%0.0f',k)
148	tmp=fseek(fid,nPxnPy(t(k)-1)*rlength,'bof');
149	for j=1:nPy
150	jx=((j-1)sNy+1):(jsNy);
151	for i=1:nPx
152	ix=((i-1)sNx+1):(isNx);
153	tm=fread(fid,1,prec);
154	[tmp count]=fread(fid,[sNx,sNy*Nz],prec);
155	tme(k)=tm;
156	fld(ix,jx,:,k)=reshape(tmp,sNx,sNy,Nz);
157	end
158	end
159	end
160
161	elseif length(Ix)<Nx \| length(Iy)<Ny
162
163	% If only part of the file needs to be read, determine
164	% which tiles contain the desired data.
165	nPx_id=zeros(Nx,Ny);
166	nPy_id=zeros(Nx,Ny);
167	for j=1:nPy
168	jx=((j-1)sNy+1):(jsNy);
169	for i=1:nPx
170	ix=((i-1)sNx+1):(isNx);
171	nPx_id(ix,jx)=i;
172	nPy_id(ix,jx)=j;
173	end
174	end
175	nPx_id=nPx_id(Ix,Iy);
176	nPy_id=nPy_id(Ix,Iy);
177	nPxy_id=[nPx_id(:) nPy_id(:)];
178	if length(nPx_id)>1
179	nPxy_id=sortrows(nPxy_id);
180	nPxy_id(find(sum(abs(diff(nPxy_id)'))'==0)+1,:)=[];
181	end
182
183	fld_tmp=zeros(Nx,Ny,Nz);
184	for k=1:length(t)
185	fprintf(1,'\b\b\b%0.0f',k)
186	for l=1:size(nPxy_id,1)
187	i=nPxy_id(l,1);
188	j=nPxy_id(l,2);
189	rnumber=nPxnPy(t(k)-1)+(j-1)*nPx+i-1;
190	tmp=fseek(fid,rnumber*rlength,'bof');
191	ix=((i-1)sNx+1):(isNx);
192	jx=((j-1)sNy+1):(jsNy);
193	tm=fread(fid,1,prec);
194	[tmp count]=fread(fid,[sNx,sNy*Nz],prec);
195	tme(k)=tm;
196	fld_tmp(ix,jx,:)=reshape(tmp,sNx,sNy,Nz);
197	end
198	fld(:,:,:,k)=fld_tmp(Ix,Iy,Iz);
199	end
200
201	elseif length(Iz)<Nz
202
203	fld_tmp=zeros(Nx,Ny,Nz);
204	for k=1:length(t)
205	fprintf(1,'\b\b\b%0.0f',k)
206	tmp=fseek(fid,nPxnPy(t(k)-1)*rlength,'bof');
207	for j=1:nPy
208	jx=((j-1)sNy+1):(jsNy);
209	for i=1:nPx
210	ix=((i-1)sNx+1):(isNx);
211	tm=fread(fid,1,prec);
212	[tmp count]=fread(fid,[sNx,sNy*Nz],prec);
213	tme(k)=tm;
214	fld_tmp(ix,jx,:)=reshape(tmp,sNx,sNy,Nz);
215	end
216	end
217	fld(:,:,:,k)=fld_tmp(:,:,Iz);
218	end
219
220	end
221
222	fid=fclose(fid);
223	fprintf(1,'\b\b\b',k)
224