cs_grid/latloncap/quikread_llc.m

function [fld fc ix jx]=quikread_llc(fnam,nx,kx,prec,gdir,minlat,maxlat,minlon,maxlon);

% Function [fld fc ix jx]=quikread_llc(fnam,nx,kx,prec, ...
%                                      gdir,minlat,maxlat,minlon,maxlon);
% Read lat-lon-cap field
% If there is less than 5 input arguments: read the complete field.
% If there is more than 5 input arguments: read a region.
%
% INPUTS
% fnam   input path and file name
% nx     tile dimension (default 270)
% kx     vertical indices to read, e.g., 1:50 (default 1)
% prec   numeric precision (see fread; default 'real*4')
% gdir   directory path name that contains grid files XC.data and YC.data
% minlat minimum latitude of region to extract
% maxlat maximum latitude of region to extract
% minlon minimum longitude of region to extract
% maxlon maximum longitude of region to extract
% (valid longitude range is -180 to 180 or 0 to 360)
%
% OUTPUTS
% fld  output array
% fc   faces that contain requested region
% ix   i-indices for requested region
% jx   j-indices for requested region
% (note that if length(fc)>1, fld, ix, and jx are matlab
%  cell arrays that can be accessed as fld{fc(1)}, etc.)
%
% NOTES
%  For faces 4 and 5, southwest velocity is:
%   u_East (i,j) =   v_Model(i,j)
%   v_North(i,j) = - u_Model(i,j-1)
%
% EXAMPLES
%
% % read and plot complete field
% fld=quikread_llc('Depth.data',270);
% quikplot_llc(fld)
%
% % read and plot the region 120W to 40W and 80S and 60N
% fld=quikread_llc('Depth.data',270,1,'real*4','',-80,60,-120,-40);
% quikpcolor(fld')
%
% SEE ALSO
% read_llc_fkij quilplot_llc quikpcolor

if nargin < 9, maxlon=180; end
if nargin < 8, minlon=-180; end
if nargin < 7, maxlat=90; end
if nargin < 6, minlat=-90; end
if nargin < 5, gdir=''; end
if nargin < 4, prec='real*4'; end
if nargin < 3, kx=1; end
if nargin < 2, nx=270; end
if nargin < 1, error('please specify input file name'); end

switch prec
 case {'integer*1'}
  prec='int8';
 case {'integer*2'}
  prec='int16';
 case {'integer*4'}
  prec='int32';
 case {'real*4','float32'}
  prec='single';
 case {'integer*8'}
  prec='int64';
 case {'real*8','float64'}
  prec='double';
end

switch prec
 case {'int8'}
  preclength=1;
 case {'int16','uint16'}
  preclength=2;
 case {'int32','uint32','single'}
  preclength=4;
 case {'int64','uint64','double'}
  preclength=8;
end

if nargin < 6
    fld=zeros(nx,nx*13,length(kx),prec);
    fid=fopen(fnam,'r','ieee-be');
    for k=1:length(kx)
        if kx(k) > 1
            skip=(kx(k)-1)*nx*nx*13;
            if(fseek(fid,skip*preclength,'bof')<0)
                error('past end of file');
            end
        end
        fld(:,:,k)=reshape(fread(fid,nx*nx*13,prec),nx,nx*13);
    end
    fid=fclose(fid);
else
    if minlon < -180
        error('minlon<-180 not yet implemented: please email menemenlis@jpl.nasa.gov')
    end
    if maxlon > 360
        error('maxlon>360 not yet implemented: please email menemenlis@jpl.nasa.gov')
    end    
    if minlat >= maxlat
        error('maxlat must be greater than minlat')
    end
    if minlon >= maxlon
        error('maxlon must be greater than minlon')
    end    
    fld=[];
    fc =[];
    ix =[];
    jx =[];
    for f=1:5
        yc=read_llc_fkij([gdir 'YC.data'],nx,f);
        xc=read_llc_fkij([gdir 'XC.data'],nx,f);
        if maxlon>180
            in=find(xc<0);
            xc(in)=xc(in)+360;
        end
        [i j]=find(yc>=minlat&yc<=maxlat&xc>=minlon&xc<=maxlon);
        if ~isempty(i)
            fc=[fc f];
            ix{f}=min(i):max(i);
            jx{f}=min(j):max(j);
            fld{f}=zeros(length(ix{f}),length(jx{f}),length(kx),prec);
            for k=1:length(kx)
                fld{f}(:,:,k)=read_llc_fkij(fnam,nx,f,kx(k),ix{f},jx{f},prec);
            end
        end
    end
    if length(fc) == 1
        fld=fld{fc};
        ix = ix{fc};
        jx = jx{fc};
    end
end
1	dimitri	1.8	function [fld fc ix jx]=quikread_llc(fnam,nx,kx,prec,gdir,minlat,maxlat,minlon,maxlon);
2	dimitri	1.1
3	dimitri	1.8	% Function [fld fc ix jx]=quikread_llc(fnam,nx,kx,prec, ...
4			% gdir,minlat,maxlat,minlon,maxlon);
5	dimitri	1.7	% Read lat-lon-cap field
6			% If there is less than 5 input arguments: read the complete field.
7			% If there is more than 5 input arguments: read a region.
8	dimitri	1.1	%
9			% INPUTS
10	dimitri	1.7	% fnam input path and file name
11			% nx tile dimension (default 270)
12			% kx vertical indices to read, e.g., 1:50 (default 1)
13			% prec numeric precision (see fread; default 'real*4')
14	dimitri	1.8	% gdir directory path name that contains grid files XC.data and YC.data
15	dimitri	1.7	% minlat minimum latitude of region to extract
16			% maxlat maximum latitude of region to extract
17			% minlon minimum longitude of region to extract
18			% maxlon maximum longitude of region to extract
19	dimitri	1.15	% (valid longitude range is -180 to 180 or 0 to 360)
20	dimitri	1.1	%
21			% OUTPUTS
22	dimitri	1.7	% fld output array
23	dimitri	1.10	% fc faces that contain requested region
24			% ix i-indices for requested region
25			% jx j-indices for requested region
26	dimitri	1.11	% (note that if length(fc)>1, fld, ix, and jx are matlab
27			% cell arrays that can be accessed as fld{fc(1)}, etc.)
28	dimitri	1.2	%
29	dimitri	1.16	% NOTES
30			% For faces 4 and 5, southwest velocity is:
31			% u_East (i,j) = v_Model(i,j)
32			% v_North(i,j) = - u_Model(i,j-1)
33			%
34	dimitri	1.2	% EXAMPLES
35	dimitri	1.4	%
36	dimitri	1.7	% % read and plot complete field
37			% fld=quikread_llc('Depth.data',270);
38			% quikplot_llc(fld)
39			%
40			% % read and plot the region 120W to 40W and 80S and 60N
41			% fld=quikread_llc('Depth.data',270,1,'real*4','',-80,60,-120,-40);
42			% quikpcolor(fld')
43	dimitri	1.8	%
44			% SEE ALSO
45			% read_llc_fkij quilplot_llc quikpcolor
46	dimitri	1.7
47			if nargin < 9, maxlon=180; end
48			if nargin < 8, minlon=-180; end
49			if nargin < 7, maxlat=90; end
50			if nargin < 6, minlat=-90; end
51	dimitri	1.8	if nargin < 5, gdir=''; end
52	dimitri	1.1	if nargin < 4, prec='real*4'; end
53			if nargin < 3, kx=1; end
54			if nargin < 2, nx=270; end
55			if nargin < 1, error('please specify input file name'); end
56
57			switch prec
58	dimitri	1.17	case {'integer*1'}
59			prec='int8';
60			case {'integer*2'}
61			prec='int16';
62			case {'integer*4'}
63			prec='int32';
64			case {'real*4','float32'}
65			prec='single';
66			case {'integer*8'}
67			prec='int64';
68			case {'real*8','float64'}
69			prec='double';
70			end
71
72			switch prec
73			case {'int8'}
74	dimitri	1.1	preclength=1;
75	dimitri	1.17	case {'int16','uint16'}
76	dimitri	1.1	preclength=2;
77	dimitri	1.17	case {'int32','uint32','single'}
78	dimitri	1.1	preclength=4;
79	dimitri	1.17	case {'int64','uint64','double'}
80	dimitri	1.1	preclength=8;
81			end
82
83	dimitri	1.7	if nargin < 6
84	dimitri	1.17	fld=zeros(nx,nx*13,length(kx),prec);
85	dimitri	1.6	fid=fopen(fnam,'r','ieee-be');
86	dimitri	1.2	for k=1:length(kx)
87	dimitri	1.7	if kx(k) > 1
88			skip=(kx(k)-1)nxnx*13;
89			if(fseek(fid,skip*preclength,'bof')<0)
90			error('past end of file');
91			end
92			end
93	dimitri	1.2	fld(:,:,k)=reshape(fread(fid,nxnx13,prec),nx,nx*13);
94	dimitri	1.7	end
95	dimitri	1.6	fid=fclose(fid);
96	dimitri	1.2	else
97	dimitri	1.13	if minlon < -180
98			error('minlon<-180 not yet implemented: please email menemenlis@jpl.nasa.gov')
99			end
100	dimitri	1.15	if maxlon > 360
101			error('maxlon>360 not yet implemented: please email menemenlis@jpl.nasa.gov')
102			end
103	dimitri	1.13	if minlat >= maxlat
104			error('maxlat must be greater than minlat')
105			end
106			if minlon >= maxlon
107			error('maxlon must be greater than minlon')
108			end
109	dimitri	1.11	fld=[];
110			fc =[];
111			ix =[];
112			jx =[];
113	dimitri	1.7	for f=1:5
114	dimitri	1.8	yc=read_llc_fkij([gdir 'YC.data'],nx,f);
115			xc=read_llc_fkij([gdir 'XC.data'],nx,f);
116	dimitri	1.15	if maxlon>180
117			in=find(xc<0);
118			xc(in)=xc(in)+360;
119			end
120	dimitri	1.7	[i j]=find(yc>=minlat&yc<=maxlat&xc>=minlon&xc<=maxlon);
121			if ~isempty(i)
122	dimitri	1.9	fc=[fc f];
123			ix{f}=min(i):max(i);
124			jx{f}=min(j):max(j);
125	dimitri	1.17	fld{f}=zeros(length(ix{f}),length(jx{f}),length(kx),prec);
126	dimitri	1.12	for k=1:length(kx)
127	dimitri	1.14	fld{f}(:,:,k)=read_llc_fkij(fnam,nx,f,kx(k),ix{f},jx{f},prec);
128	dimitri	1.12	end
129	dimitri	1.2	end
130			end
131	dimitri	1.9	if length(fc) == 1
132	dimitri	1.11	fld=fld{fc};
133			ix = ix{fc};
134			jx = jx{fc};
135	dimitri	1.7	end
136	dimitri	1.1	end