Back to home page

MITgcm

 
 

    


Warning, /verification/natl_box/matlab/readgcm.m is written in an unsupported language. File is not indexed.

view on githubraw file Latest commit cdf4d124 on 2000-11-13 16:02:33 UTC
cdf4d1244c Patr*0001 function [fld,tme,nPx,nPy]=readgcm(fnam,t,Nx,Ny,Nz,Ix,Iy,Iz,prec,nPx,nPy)
                0002 
                0003 % Function [fld,tme,nPx,nPy]=readgcm(fnam,t,Nx,Ny,Nz,Ix,Iy,Iz,prec,nPx,nPy)
                0004 % read and output model field from file fnam time step t
                0005 %
                0006 % INPUTS
                0007 % fnam      input path and file name
                0008 % t         indices of time steps to read (0 means read all, the default)
                0009 % Nx*Ny*Nz  grid dimension (default 360*224*46)
                0010 % Ix,Iy,Iz  subsample indices (default Ix=1:Nx, Iy=1:Ny, and Iz=1:Nz)
                0011 % prec      numeric precision (default 'real*4')
                0012 %
                0013 % INPUT/OUTPUT
                0014 % nPx, nPy  number of processes in x and y direction;
                0015 %           when not specified as INPUTS, nPx and nPy are determined based
                0016 %           on the following file naming convention:
                0017 %             U_06_04.00060_00120_10 is
                0018 %             U for nPx=6, nPy=4, hour 60 to 120 at 10-hour intervals
                0019 %           (defaults to 1 when filename contains no underscores)
                0020 %
                0021 % OUTPUTS
                0022 % fld   Nx*Ny*Nz*Nt output array
                0023 %          e.g., to contour level k, use "contourf(fld(:,:,k,1)'), colorbar"
                0024 % tme   model integration times in days
                0025 %
                0026 % The following are valid calling forms:
                0027 %    readgcm(fnam,t,Nx,Ny,Nz,Ix,Iy,Iz,prec,nPx,nPy)
                0028 %    readgcm(fnam, ...                            )
                0029 %    readgcm(fnam, ...               ,prec, ...   )
                0030 % where "..." represents 0 or more ordered arguments.
                0031 %
                0032 % SEE ALSO
                0033 % writegcm
                0034 
                0035 % "prec" is set to value of second string argument if it exists
                0036 nprec=9;
                0037 if exist('t') , if ischar(t)
                0038     prec=t;  nprec=2;
                0039     if exist('Nx'), nPx=Nx; end
                0040     if exist('Ny'), nPy=Ny; end
                0041 end, end
                0042 if exist('Nx'), if ischar(Nx)
                0043     prec=Nx; nprec=3;
                0044     if exist('Ny'), nPx=Ny; end
                0045     if exist('Nz'), nPy=Nz; end
                0046 end, end
                0047 if exist('Ny'), if ischar(Ny)
                0048     prec=Ny; nprec=4;
                0049     if exist('Nz'), nPx=Nz; end
                0050     if exist('Ix'), nPy=Ix; end
                0051 end, end
                0052 if exist('Nz'), if ischar(Nz)
                0053     prec=Nz; nprec=5;
                0054     if exist('Ix'), nPx=Ix; end
                0055     if exist('Iy'), nPy=Iy; end
                0056 end, end
                0057 if exist('Ix'), if ischar(Ix)
                0058     prec=Ix; nprec=6;
                0059     if exist('Iy'), nPx=Iy; end
                0060     if exist('Iz'), nPy=Iz; end
                0061 end, end
                0062 if exist('Iy'), if ischar(Iy)
                0063     prec=Iy; nprec=7;
                0064     if exist('Iz'), nPx=Iz; end
                0065     if exist('prec'), nPy=prec; end
                0066 end, end
                0067 if exist('Iz'), if ischar(Iz)
                0068     prec=Iz; nprec=8;
                0069     if exist('prec'), nPx=prec; end
                0070     if exist('nPx'), nPy=nPx; end
                0071 end, end
                0072 
                0073 % set default argument values
                0074 if nargin < 9 & nprec==9, prec='real*4'; end
                0075 if nargin < 5 | nprec <= 5
                0076   if isempty(findstr(fnam,'KPPhbl')) & ...
                0077         isempty(findstr(fnam,'H')) & ...
                0078         isempty(findstr(fnam,'SF')) & ...
                0079         isempty(findstr(fnam,'Pbottom')) & ...
                0080         isempty(findstr(fnam,'Pbottom')) & ...
                0081         isempty(findstr(fnam,'BU')) & ...
                0082         isempty(findstr(fnam,'BV'))
                0083     Nz=46;
                0084   else 
                0085     Nz=1;
                0086   end
                0087 end
                0088 if nargin < 4 | nprec <= 4, Ny=224; end
                0089 if nargin < 3 | nprec <= 3, Nx=360; end
                0090 if nargin < 8 | nprec <= 8, Iz=1:Nz; end
                0091 if nargin < 7 | nprec <= 7, Iy=1:Ny; end
                0092 if nargin < 6 | nprec <= 6, Ix=1:Nx; end
                0093 if nargin < 2 | nprec <= 2, t=0; end
                0094 if nargin < 1, error('please specify input file name'); end
                0095 if nargin < nprec+1
                0096   % locate beginning filename character in path
                0097   slash_loc=findstr('/',fnam);
                0098   if isempty(slash_loc), slash_loc=0; end
                0099   slash_loc=max(slash_loc);
                0100   bar_loc=findstr('_',fnam((slash_loc+1):length(fnam)));
                0101   bar_loc=min(bar_loc);
                0102   if isempty(bar_loc)
                0103     nPx=1; nPy=1;
                0104   elseif length(fnam)<slash_loc+bar_loc+5
                0105     nPx=1; nPy=1;
                0106   else
                0107     nPx=str2num(fnam(slash_loc+bar_loc+(1:2)));  % number of processes in X
                0108     nPy=str2num(fnam(slash_loc+bar_loc+(4:5)));  % number of processes in Y
                0109     if isempty(nPx) | isempty(nPy), nPx=1; nPy=1; end
                0110   end
                0111 elseif nargin < nprec+2
                0112   nPy=nPx;                                  % number of processes in X and Y
                0113 end
                0114 
                0115 % compute number of sub-grids
                0116 sNx=Nx/nPx;                                 % number of X points in sub-grid
                0117 sNy=Ny/nPy;                                 % number of Y points in sub-grid
                0118 if (sNx-floor(sNx))~=0 | (sNy-floor(sNy))~=0
                0119   error('Nx/nPx and Ny/nPy must be integer')
                0120 end
                0121 
                0122 % compute file record length
                0123 switch prec
                0124   case {'float32', 'real*4'}
                0125     rlength=(sNx*sNy*Nz+1)*4;
                0126   case {'float64', 'real*8'}
                0127     rlength=(sNx*sNy*Nz+1)*8;
                0128 end
                0129 
                0130 % compute total number of time steps
                0131 fid=fopen(fnam,'r','ieee-be');
                0132 tmp=fseek(fid,0,'eof');
                0133 tmp=ftell(fid);
                0134 tmp=tmp/(rlength*nPx*nPy);
                0135 if t==0 | tmp<max(t)
                0136   t=1:tmp;
                0137 end
                0138 
                0139 % read model file
                0140 fld=zeros(length(Ix),length(Iy),length(Iz),length(t));
                0141 tme=zeros(1,length(t));
                0142 disp(['Reading ' int2str(length(t)) ' time steps:'])
                0143 
                0144 if length(Ix)==Nx & length(Iy)==Ny & length(Iz)==Nz
                0145 
                0146   for k=1:length(t)
                0147     fprintf(1,'\b\b\b%0.0f',k)
                0148     tmp=fseek(fid,nPx*nPy*(t(k)-1)*rlength,'bof');
                0149     for j=1:nPy
                0150       jx=((j-1)*sNy+1):(j*sNy);
                0151       for i=1:nPx
                0152         ix=((i-1)*sNx+1):(i*sNx);
                0153         tm=fread(fid,1,prec);
                0154         [tmp count]=fread(fid,[sNx,sNy*Nz],prec);
                0155         tme(k)=tm;
                0156         fld(ix,jx,:,k)=reshape(tmp,sNx,sNy,Nz);
                0157       end
                0158     end
                0159   end
                0160   
                0161 elseif length(Ix)<Nx | length(Iy)<Ny
                0162 
                0163   % If only part of the file needs to be read, determine
                0164   % which tiles contain the desired data.
                0165   nPx_id=zeros(Nx,Ny);
                0166   nPy_id=zeros(Nx,Ny);
                0167   for j=1:nPy
                0168     jx=((j-1)*sNy+1):(j*sNy);
                0169     for i=1:nPx
                0170       ix=((i-1)*sNx+1):(i*sNx);
                0171       nPx_id(ix,jx)=i;
                0172       nPy_id(ix,jx)=j;
                0173     end
                0174   end
                0175   nPx_id=nPx_id(Ix,Iy);
                0176   nPy_id=nPy_id(Ix,Iy);
                0177   nPxy_id=[nPx_id(:) nPy_id(:)];
                0178   if length(nPx_id)>1
                0179     nPxy_id=sortrows(nPxy_id);
                0180     nPxy_id(find(sum(abs(diff(nPxy_id)'))'==0)+1,:)=[];
                0181   end
                0182 
                0183   fld_tmp=zeros(Nx,Ny,Nz);
                0184   for k=1:length(t)
                0185     fprintf(1,'\b\b\b%0.0f',k)
                0186     for l=1:size(nPxy_id,1)
                0187       i=nPxy_id(l,1);
                0188       j=nPxy_id(l,2);
                0189       rnumber=nPx*nPy*(t(k)-1)+(j-1)*nPx+i-1;
                0190       tmp=fseek(fid,rnumber*rlength,'bof');
                0191       ix=((i-1)*sNx+1):(i*sNx);
                0192       jx=((j-1)*sNy+1):(j*sNy);
                0193       tm=fread(fid,1,prec);
                0194       [tmp count]=fread(fid,[sNx,sNy*Nz],prec);
                0195       tme(k)=tm;
                0196       fld_tmp(ix,jx,:)=reshape(tmp,sNx,sNy,Nz);
                0197     end
                0198     fld(:,:,:,k)=fld_tmp(Ix,Iy,Iz);
                0199   end
                0200   
                0201 elseif length(Iz)<Nz
                0202 
                0203   fld_tmp=zeros(Nx,Ny,Nz);
                0204   for k=1:length(t)
                0205     fprintf(1,'\b\b\b%0.0f',k)
                0206     tmp=fseek(fid,nPx*nPy*(t(k)-1)*rlength,'bof');
                0207     for j=1:nPy
                0208       jx=((j-1)*sNy+1):(j*sNy);
                0209       for i=1:nPx
                0210         ix=((i-1)*sNx+1):(i*sNx);
                0211         tm=fread(fid,1,prec);
                0212         [tmp count]=fread(fid,[sNx,sNy*Nz],prec);
                0213         tme(k)=tm;
                0214         fld_tmp(ix,jx,:)=reshape(tmp,sNx,sNy,Nz);
                0215       end
                0216     end
                0217     fld(:,:,:,k)=fld_tmp(:,:,Iz);
                0218   end
                0219 
                0220 end
                0221 
                0222 fid=fclose(fid);
                0223 fprintf(1,'\b\b\b',k)
                0224