natl_box/matlab/writegcm.m

cdf4d1244c Patr*0001 function writegcm(fnam,fld,tme,nPx,nPy,prec,oflag)
                0002
                0003 % Function writegcm(fnam,fld,tme,nPx,nPy,prec,oflag)
                0004 % write array fld in MIT GCM UV format
                0005 %
                0006 % INPUTS
                0007 % fnam      output file name
                0008 % fld       Nx*Ny*Nz*Nt input array
                0009 % tme       model integration times in days has dimension Nt
                0010 % nPx, nPy  number of processes in x and y direction (default 1,1)
                0011 % prec      numeric precision (default 'real*4')
                0012 % oflag     overwrite flag specifies field to be overwritten
                0013 %           oflag=0 means create new file (the default)
                0014 %           oflag>0 must have dimension Nt
                0015 %
                0016 % SEE ALSO
                0017 % readgcm
                0018
                0019 if nargin < 2, error('please specify array and output file name'); end
                0020 [Nx Ny Nz Nt]=size(fld);
                0021 if nargin < 7, oflag=0; end
                0022 if oflag > 0
                0023   if length(oflag) ~= Nt, error('length(oflag) must equal Nt'); end
                0024 end
                0025 if nargin < 6, prec='real*4'; end
                0026 if nargin < 5, nPy=1; end
                0027 if nargin < 4, nPx=1; end
                0028 if nargin < 3, tme=1:Nt; end
                0029 if length(tme) ~= Nt, error('length(tme) must equal Nt'); end
                0030
                0031 sNx=Nx/nPx;                                 % number of X points in sub-grid
                0032 sNy=Ny/nPy;                                 % number of Y points in sub-grid
                0033 if (sNx-floor(sNx))~=0 | (sNy-floor(sNy))~=0
                0034   error('Nx/nPx and Ny/nPy must be integer')
                0035 end
                0036
                0037 switch prec
                0038   case {'float32', 'real*4'}
                0039     rlength=(sNx*sNy*Nz+1)*4;
                0040   case {'float64', 'real*8'}
                0041     rlength=(sNx*sNy*Nz+1)*8;
                0042 end
                0043
                0044 if oflag == 0
                0045
                0046   fid=fopen(fnam,'w','ieee-be');
                0047   for t=1:Nt
                0048     for j=1:nPy
                0049       jx=((j-1)*sNy+1):(j*sNy);
                0050       for i=1:nPx
                0051         ix=((i-1)*sNx+1):(i*sNx);
                0052         fwrite(fid,tme(t),prec);
                0053         fwrite(fid,fld(ix,jx,:,t),prec);
                0054       end
                0055     end
                0056   end
                0057
                0058 else
                0059
                0060   fid=fopen(fnam,'r+','ieee-be');
                0061   for t=1:Nt
                0062     if fseek(fid,nPx*nPy*(oflag(t)-1)*rlength,'bof')
                0063       error('unexpected end of file reached')
                0064     end
                0065     for j=1:nPy
                0066       jx=((j-1)*sNy+1):(j*sNy);
                0067       for i=1:nPx
                0068         ix=((i-1)*sNx+1):(i*sNx);
                0069         fwrite(fid,tme(t),prec);
                0070         fwrite(fid,fld(ix,jx,:,t),prec);
                0071       end
                0072     end
                0073   end
                0074
                0075 end
                0076
                0077 fid=fclose(fid);