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 |
% Nx*Ny*Nz grid dimension (default 360*224*46) |
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 Nx*Ny*Nz*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=(sNx*sNy*Nz+1)*4; |
126 |
case {'float64', 'real*8'} |
127 |
rlength=(sNx*sNy*Nz+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/(rlength*nPx*nPy); |
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,nPx*nPy*(t(k)-1)*rlength,'bof'); |
149 |
for j=1:nPy |
150 |
jx=((j-1)*sNy+1):(j*sNy); |
151 |
for i=1:nPx |
152 |
ix=((i-1)*sNx+1):(i*sNx); |
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):(j*sNy); |
169 |
for i=1:nPx |
170 |
ix=((i-1)*sNx+1):(i*sNx); |
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=nPx*nPy*(t(k)-1)+(j-1)*nPx+i-1; |
190 |
tmp=fseek(fid,rnumber*rlength,'bof'); |
191 |
ix=((i-1)*sNx+1):(i*sNx); |
192 |
jx=((j-1)*sNy+1):(j*sNy); |
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,nPx*nPy*(t(k)-1)*rlength,'bof'); |
207 |
for j=1:nPy |
208 |
jx=((j-1)*sNy+1):(j*sNy); |
209 |
for i=1:nPx |
210 |
ix=((i-1)*sNx+1):(i*sNx); |
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 |
|