function [meantime, gdata] = readl2std(fn);

% function [meantime, gdata] = readl2std(fn);
%
% Reads an AIRS level 2 Standard retrieval granule file and returns
% an RTP-like structure of retreived profiles.
%
% Input:
%    fn = (string) Name of an AIRS l1b granule file, something like
%          'AIRS.2003.08.31.067.L2.RetStd.v3.1.9.0.A03261143158.hdf'
%
% Output:
%    meantime = (1x 1) mean time of observations, in hours
%    gdata = (structure) RTP "prof" like structure
%
% Note: if the granule contains no good data, the output variables
% are returned empty.
%

% Created: 01 October 2003, Scott Hannon - based on readl1b_all.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Granule dimensions
nxtrack=30;
natrack=45;
nobs=nxtrack*natrack;
nlevs=28;

% Open granule file
file_name = fn;
file_id   = hdfsw('open',file_name,'read');
%%%
% Uncomment the line below to see what swath names are found in file
% [NSWATH,SWATHLIST] = hdfsw('inqswath',fn)
%%%
swath_id  = hdfsw('attach',file_id,'L2_Standard_atmospheric&surface_product');


% Read "RetQAFlag" and find good FOVs
[junk,s]=hdfsw('readfield',swath_id,'RetQAFlag',[],[],[]);
if s == -1; disp('Error reading RetQAFlag');end;
state = reshape( double(junk), 1,nobs);
i0=find( state == 0);  % Indices of "good" FOVs
n0=length(i0);
%
clear state


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if (n0 > 0)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%
% Uncomment the line below to inquire about attribute names
% [NATTR,ATTRLIST] = hdfsw('inqattrs',swath_id)
%%%

% Read the std pressure levels
[junk,s]=hdfsw('readattr',swath_id,'pressStd');
if (s == -1)
   % pressStd is missing from file!
   disp('WARNING! HDF file is missing pressStd; setting to zeroes')
   plevs=zeros(nlevs,1);
else
   plevs = double(junk);
end


% Read the date/time fields
[junk,s]=hdfsw('readattr',swath_id,'start_year');
l2.start_year = double(junk(1));
%
[junk,s]=hdfsw('readattr',swath_id,'start_month');
l2.start_month = double(junk(1));
%
[junk,s]=hdfsw('readattr',swath_id,'start_day');
l2.start_day = double(junk(1));
%
[junk,s]=hdfsw('readattr',swath_id,'start_hour');
l2.start_hour = double(junk(1));
%
[junk,s]=hdfsw('readattr',swath_id,'start_minute');
l2.start_minute = double(junk(1));
%
[junk,s]=hdfsw('readattr',swath_id,'start_sec');
l2.start_sec = double(junk(1));
%
[junk,s]=hdfsw('readattr',swath_id,'start_Time');
l2.start_Time = double(junk(1));
%
[junk,s]=hdfsw('readattr',swath_id,'end_Time');
l2.end_Time = double(junk(1));


% Compute approximate time for each FOV
deltasec=(l2.end_Time - l2.start_Time);
dtime=deltasec/(nobs - 1);
time=zeros(nxtrack,natrack);
indx=(1:nxtrack)';
for ia=1:natrack
   time(:,ia)=l2.start_Time + ((ia-1)*nxtrack + indx)*dtime; % approximate
end
meantime=(0.5*deltasec + l2.start_sec)/3600 + l2.start_minute/60 + ...
   l2.start_hour;
%
clear l2 deltasec dtime indx


% Read per scanline fields; expand to per FOV later
%
% satheight (1 x natrack)
[junk,s]=hdfsw('readfield',swath_id,'satheight',[],[],[]);
if s == -1; disp('Error reading satheight'); end;
satheight = double(junk');


% Declare temporary variables for expansion
tmp_atrack=zeros(1,nobs);
tmp_xtrack=zeros(1,nobs);
tmp_rtime=zeros(1,nobs);
tmp_zobs=zeros(1,nobs);


% Loop over along-track and fill in temporary variables
ix=1:nxtrack;
for ia=1:natrack
   iobs=nxtrack*(ia-1) + ix;
   %
   % Fill in cross-track
   tmp_atrack(iobs)=ia*3 - 1;  % convert AMSU index to AIRS center FOV
   tmp_xtrack(iobs)=ix*3 - 1;  % convert AMSU index to AIRS center FOV
   tmp_rtime(iobs)=time(:,ia);
   tmp_zobs(iobs)=satheight(ia)*1000;  % convert km to meters
end
%
clear ix ia iobs time satheight


% Subset temporary variables for state and re-assign to gdata
gdata.atrack = tmp_atrack(i0);
gdata.xtrack = tmp_xtrack(i0);
gdata.rtime  = tmp_rtime(i0);
gdata.zobs   = tmp_zobs(i0);
%
clear tmp_atrack tmp_xtrack tmp_rtime tmp_zobs

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Read the per FOV instrument & geolocation data
%
[junk,s]=hdfsw('readfield',swath_id,'latAIRS',[],[],[]);
if s == -1; disp('Error reading latitude');end;
junk2 = reshape( double(junk), 9,nobs);  % All 9 AIRS FOVS
gdata.rlat = junk2(5,i0); % Center AIRS FOV
%
[junk,s]=hdfsw('readfield',swath_id,'lonAIRS',[],[],[]);
if s == -1; disp('Error reading longitude');end;
junk2 = reshape( double(junk), 9,nobs); % All 9 AIRS FOVs
gdata.rlon = junk2(5,i0); % Center AIRS FOV
%
[junk,s]=hdfsw('readfield',swath_id,'satzen',[],[],[]);
if s == -1; disp('Error reading satzen');end;
junk2 = reshape( double(junk), 1,nobs);
gdata.satzen = junk2(i0);
%
[junk,s]=hdfsw('readfield',swath_id,'satazi',[],[],[]);
if s == -1; disp('Error reading satazi');end;
junk2 = reshape( double(junk), 1,nobs);
gdata.satazi = junk2(i0);
%
[junk,s]=hdfsw('readfield',swath_id,'solzen',[],[],[]);
if s == -1; disp('Error reading solzen');end;
junk2 = reshape( double(junk), 1,nobs);
gdata.solzen = junk2(i0);
%
[junk,s]=hdfsw('readfield',swath_id,'solazi',[],[],[]);
if s == -1; disp('Error reading solazi');end;
junk2 = reshape( double(junk), 1,nobs);
gdata.solazi = junk2(i0);
%
[junk,s]=hdfsw('readfield',swath_id,'topog',[],[],[]);
if s == -1; disp('Error reading topog');end;
junk2 = reshape( double(junk), 1,nobs);
gdata.salti =junk2(i0);
%
[junk,s]=hdfsw('readfield',swath_id,'landFrac',[],[],[]);
if s == -1; disp('Error reading landFrac');end;
junk2 = reshape( double(junk), 1,nobs);
gdata.landfrac = junk2(i0);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Read retrieved profile data
%
gdata.plevs = plevs*ones(1,length(i0));
%
[junk,s]=hdfsw('readfield',swath_id,'PSurfStd',[],[],[]);
if s == -1; disp('Error reading PSurfStd');end;
junk2 = reshape( double(junk), 1,nobs);
gdata.spres = junk2(i0);
%
[junk,s]=hdfsw('readfield',swath_id,'TSurfStd',[],[],[]);
if s == -1; disp('Error reading TSurfStd');end;
junk2 = reshape( double(junk), 1,nobs);
gdata.stemp = junk2(i0);
%
[junk,s]=hdfsw('readfield',swath_id,'nSurfStd',[],[],[]);
if s == -1; disp('Error reading nSurfStd');end;
junk2 = reshape( double(junk), 1,nobs);
%%%
%gdata.nlevs = junk2(i0);
% Since pressStd is broken, so is nlevs.  Set to 28
gdata.nlevs = 28*ones(1,length(i0));
%%%

%
[junk,s]=hdfsw('readfield',swath_id,'TAirStd',[],[],[]);
if s == -1; disp('Error reading TAirStd');end;
junk2 = reshape( double(junk), nlevs,nobs);
gdata.ptemp = junk2(:,i0);
%
[junk,s]=hdfsw('readfield',swath_id,'H2OMMRStd',[],[],[]);
if s == -1; disp('Error reading H2OMMRStd');end;
junk2 = reshape( double(junk), nlevs,nobs);
gdata.gas_1 = junk2(:,i0);
%
[junk,s]=hdfsw('readfield',swath_id,'O3VMRStd',[],[],[]);
if s == -1; disp('Error reading O3VMRStd');end;
junk2 = reshape( double(junk), nlevs,nobs);
gdata.gas_3 = junk2(:,i0);
%
clear junk junk2 i0


% Close L2 Std Ret granule file
s = hdfsw('detach',swath_id);
if s == -1; disp('Swatch detach error: L2');end;   
s = hdfsw('close',file_id);
if s == -1; disp('File close error: L2');end;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
else
   disp('No good FOVs in L2 granule file:')
   disp(fn)

   meantime=[];
   f=[];
   gdata=[];

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% end of function %%%
