1 |
function writegcm(fnam,fld,tme,nPx,nPy,prec,oflag) |
2 |
|
3 |
% Function writegcm(fnam,fld,tme,nPx,nPy,prec,oflag) |
4 |
% write array fld in MIT GCM UV format |
5 |
% |
6 |
% INPUTS |
7 |
% fnam output file name |
8 |
% fld Nx*Ny*Nz*Nt input array |
9 |
% tme model integration times in days has dimension Nt |
10 |
% nPx, nPy number of processes in x and y direction (default 1,1) |
11 |
% prec numeric precision (default 'real*4') |
12 |
% oflag overwrite flag specifies field to be overwritten |
13 |
% oflag=0 means create new file (the default) |
14 |
% oflag>0 must have dimension Nt |
15 |
% |
16 |
% SEE ALSO |
17 |
% readgcm |
18 |
|
19 |
if nargin < 2, error('please specify array and output file name'); end |
20 |
[Nx Ny Nz Nt]=size(fld); |
21 |
if nargin < 7, oflag=0; end |
22 |
if oflag > 0 |
23 |
if length(oflag) ~= Nt, error('length(oflag) must equal Nt'); end |
24 |
end |
25 |
if nargin < 6, prec='real*4'; end |
26 |
if nargin < 5, nPy=1; end |
27 |
if nargin < 4, nPx=1; end |
28 |
if nargin < 3, tme=1:Nt; end |
29 |
if length(tme) ~= Nt, error('length(tme) must equal Nt'); end |
30 |
|
31 |
sNx=Nx/nPx; % number of X points in sub-grid |
32 |
sNy=Ny/nPy; % number of Y points in sub-grid |
33 |
if (sNx-floor(sNx))~=0 | (sNy-floor(sNy))~=0 |
34 |
error('Nx/nPx and Ny/nPy must be integer') |
35 |
end |
36 |
|
37 |
switch prec |
38 |
case {'float32', 'real*4'} |
39 |
rlength=(sNx*sNy*Nz+1)*4; |
40 |
case {'float64', 'real*8'} |
41 |
rlength=(sNx*sNy*Nz+1)*8; |
42 |
end |
43 |
|
44 |
if oflag == 0 |
45 |
|
46 |
fid=fopen(fnam,'w','ieee-be'); |
47 |
for t=1:Nt |
48 |
for j=1:nPy |
49 |
jx=((j-1)*sNy+1):(j*sNy); |
50 |
for i=1:nPx |
51 |
ix=((i-1)*sNx+1):(i*sNx); |
52 |
fwrite(fid,tme(t),prec); |
53 |
fwrite(fid,fld(ix,jx,:,t),prec); |
54 |
end |
55 |
end |
56 |
end |
57 |
|
58 |
else |
59 |
|
60 |
fid=fopen(fnam,'r+','ieee-be'); |
61 |
for t=1:Nt |
62 |
if fseek(fid,nPx*nPy*(oflag(t)-1)*rlength,'bof') |
63 |
error('unexpected end of file reached') |
64 |
end |
65 |
for j=1:nPy |
66 |
jx=((j-1)*sNy+1):(j*sNy); |
67 |
for i=1:nPx |
68 |
ix=((i-1)*sNx+1):(i*sNx); |
69 |
fwrite(fid,tme(t),prec); |
70 |
fwrite(fid,fld(ix,jx,:,t),prec); |
71 |
end |
72 |
end |
73 |
end |
74 |
|
75 |
end |
76 |
|
77 |
fid=fclose(fid); |