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	heimbach	1.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