natl_box_forward/matlab/writegcm.m

function writegcm(fnam,fld,tme,nPx,nPy,prec,oflag)

% Function writegcm(fnam,fld,tme,nPx,nPy,prec,oflag)
% write array fld in MIT GCM UV format
%
% INPUTS
% fnam      output file name
% fld       Nx*Ny*Nz*Nt input array
% tme       model integration times in days has dimension Nt
% nPx, nPy  number of processes in x and y direction (default 1,1)
% prec      numeric precision (default 'real*4')
% oflag     overwrite flag specifies field to be overwritten
%           oflag=0 means create new file (the default)
%           oflag>0 must have dimension Nt
%
% SEE ALSO
% readgcm

if nargin < 2, error('please specify array and output file name'); end
[Nx Ny Nz Nt]=size(fld);
if nargin < 7, oflag=0; end
if oflag > 0
  if length(oflag) ~= Nt, error('length(oflag) must equal Nt'); end
end
if nargin < 6, prec='real*4'; end
if nargin < 5, nPy=1; end
if nargin < 4, nPx=1; end
if nargin < 3, tme=1:Nt; end
if length(tme) ~= Nt, error('length(tme) must equal Nt'); end

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

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

if oflag == 0

  fid=fopen(fnam,'w','ieee-be');
  for t=1:Nt
    for j=1:nPy
      jx=((j-1)*sNy+1):(j*sNy);
      for i=1:nPx
        ix=((i-1)*sNx+1):(i*sNx);
        fwrite(fid,tme(t),prec);
        fwrite(fid,fld(ix,jx,:,t),prec);
      end
    end
  end

else

  fid=fopen(fnam,'r+','ieee-be');
  for t=1:Nt
    if fseek(fid,nPx*nPy*(oflag(t)-1)*rlength,'bof')
      error('unexpected end of file reached')
    end
    for j=1:nPy
      jx=((j-1)*sNy+1):(j*sNy);
      for i=1:nPx
        ix=((i-1)*sNx+1):(i*sNx);
        fwrite(fid,tme(t),prec);
        fwrite(fid,fld(ix,jx,:,t),prec);
      end
    end
  end

end

fid=fclose(fid);
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 NxNyNz*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=(sNxsNyNz+1)*4;
40	case {'float64', 'real*8'}
41	rlength=(sNxsNyNz+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):(jsNy);
50	for i=1:nPx
51	ix=((i-1)sNx+1):(isNx);
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,nPxnPy(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):(jsNy);
67	for i=1:nPx
68	ix=((i-1)sNx+1):(isNx);
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);