Warning, /utils/matlab/cs_grid/latloncap/quikread_llc.m is written in an unsupported language. File is not indexed.
view on githubraw file Latest commit f91ed0bf on 2016-07-07 16:26:51 UTC
a61e411d79 Dimi*0001 function [fld fc ix jx]=quikread_llc(fnam,nx,kx,prec,gdir,minlat,maxlat,minlon,maxlon);
cb1c045303 Dimi*0002
a61e411d79 Dimi*0003 % Function [fld fc ix jx]=quikread_llc(fnam,nx,kx,prec, ...
0004 % gdir,minlat,maxlat,minlon,maxlon);
c025ca051f Dimi*0005 % Read lat-lon-cap field
0006 % If there is less than 5 input arguments: read the complete field.
0007 % If there is more than 5 input arguments: read a region.
cb1c045303 Dimi*0008 %
0009 % INPUTS
c025ca051f Dimi*0010 % fnam input path and file name
0011 % nx tile dimension (default 270)
0012 % kx vertical indices to read, e.g., 1:50 (default 1)
0013 % prec numeric precision (see fread; default 'real*4')
a61e411d79 Dimi*0014 % gdir directory path name that contains grid files XC.data and YC.data
c025ca051f Dimi*0015 % minlat minimum latitude of region to extract
0016 % maxlat maximum latitude of region to extract
0017 % minlon minimum longitude of region to extract
0018 % maxlon maximum longitude of region to extract
a8cdeb39a4 Dimi*0019 % (valid longitude range is -180 to 180 or 0 to 360)
cb1c045303 Dimi*0020 %
0021 % OUTPUTS
c025ca051f Dimi*0022 % fld output array
0f1ad8376d Dimi*0023 % fc faces that contain requested region
0024 % ix i-indices for requested region
0025 % jx j-indices for requested region
9fbbb66736 Dimi*0026 % (note that if length(fc)>1, fld, ix, and jx are matlab
0027 % cell arrays that can be accessed as fld{fc(1)}, etc.)
c025ca051f Dimi*0028 %
e6b73c8b78 Dimi*0029 % NOTES
0030 % For faces 4 and 5, southwest velocity is:
0031 % u_East (i,j) = v_Model(i,j)
0032 % v_North(i,j) = - u_Model(i,j-1)
0033 %
c025ca051f Dimi*0034 % EXAMPLES
0035 %
0036 % % read and plot complete field
0037 % fld=quikread_llc('Depth.data',270);
0038 % quikplot_llc(fld)
0039 %
0040 % % read and plot the region 120W to 40W and 80S and 60N
0041 % fld=quikread_llc('Depth.data',270,1,'real*4','',-80,60,-120,-40);
0042 % quikpcolor(fld')
a61e411d79 Dimi*0043 %
0044 % SEE ALSO
0045 % read_llc_fkij quilplot_llc quikpcolor
cb1c045303 Dimi*0046
c025ca051f Dimi*0047 if nargin < 9, maxlon=180; end
0048 if nargin < 8, minlon=-180; end
0049 if nargin < 7, maxlat=90; end
0050 if nargin < 6, minlat=-90; end
a61e411d79 Dimi*0051 if nargin < 5, gdir=''; end
cb1c045303 Dimi*0052 if nargin < 4, prec='real*4'; end
0053 if nargin < 3, kx=1; end
0054 if nargin < 2, nx=270; end
0055 if nargin < 1, error('please specify input file name'); end
0056
0057 switch prec
f91ed0bfc5 Dimi*0058 case {'integer*1'}
0059 prec='int8';
0060 case {'integer*2'}
0061 prec='int16';
0062 case {'integer*4'}
0063 prec='int32';
0064 case {'real*4','float32'}
0065 prec='single';
0066 case {'integer*8'}
0067 prec='int64';
0068 case {'real*8','float64'}
0069 prec='double';
0070 end
0071
0072 switch prec
0073 case {'int8'}
cb1c045303 Dimi*0074 preclength=1;
f91ed0bfc5 Dimi*0075 case {'int16','uint16'}
cb1c045303 Dimi*0076 preclength=2;
f91ed0bfc5 Dimi*0077 case {'int32','uint32','single'}
cb1c045303 Dimi*0078 preclength=4;
f91ed0bfc5 Dimi*0079 case {'int64','uint64','double'}
cb1c045303 Dimi*0080 preclength=8;
0081 end
0082
c025ca051f Dimi*0083 if nargin < 6
f91ed0bfc5 Dimi*0084 fld=zeros(nx,nx*13,length(kx),prec);
c4362e6083 Dimi*0085 fid=fopen(fnam,'r','ieee-be');
ebf5c557d5 Dimi*0086 for k=1:length(kx)
c025ca051f Dimi*0087 if kx(k) > 1
0088 skip=(kx(k)-1)*nx*nx*13;
0089 if(fseek(fid,skip*preclength,'bof')<0)
0090 error('past end of file');
0091 end
0092 end
ebf5c557d5 Dimi*0093 fld(:,:,k)=reshape(fread(fid,nx*nx*13,prec),nx,nx*13);
c025ca051f Dimi*0094 end
c4362e6083 Dimi*0095 fid=fclose(fid);
ebf5c557d5 Dimi*0096 else
1c5d463102 Dimi*0097 if minlon < -180
0098 error('minlon<-180 not yet implemented: please email menemenlis@jpl.nasa.gov')
0099 end
a8cdeb39a4 Dimi*0100 if maxlon > 360
0101 error('maxlon>360 not yet implemented: please email menemenlis@jpl.nasa.gov')
0102 end
1c5d463102 Dimi*0103 if minlat >= maxlat
0104 error('maxlat must be greater than minlat')
0105 end
0106 if minlon >= maxlon
0107 error('maxlon must be greater than minlon')
0108 end
9fbbb66736 Dimi*0109 fld=[];
0110 fc =[];
0111 ix =[];
0112 jx =[];
c025ca051f Dimi*0113 for f=1:5
a61e411d79 Dimi*0114 yc=read_llc_fkij([gdir 'YC.data'],nx,f);
0115 xc=read_llc_fkij([gdir 'XC.data'],nx,f);
a8cdeb39a4 Dimi*0116 if maxlon>180
0117 in=find(xc<0);
0118 xc(in)=xc(in)+360;
0119 end
c025ca051f Dimi*0120 [i j]=find(yc>=minlat&yc<=maxlat&xc>=minlon&xc<=maxlon);
0121 if ~isempty(i)
41dfd48a39 Dimi*0122 fc=[fc f];
0123 ix{f}=min(i):max(i);
0124 jx{f}=min(j):max(j);
f91ed0bfc5 Dimi*0125 fld{f}=zeros(length(ix{f}),length(jx{f}),length(kx),prec);
7985d5082e Dimi*0126 for k=1:length(kx)
b2eacc24c6 Dimi*0127 fld{f}(:,:,k)=read_llc_fkij(fnam,nx,f,kx(k),ix{f},jx{f},prec);
7985d5082e Dimi*0128 end
ebf5c557d5 Dimi*0129 end
0130 end
41dfd48a39 Dimi*0131 if length(fc) == 1
9fbbb66736 Dimi*0132 fld=fld{fc};
0133 ix = ix{fc};
0134 jx = jx{fc};
c025ca051f Dimi*0135 end
cb1c045303 Dimi*0136 end
