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
