% Read sequential GRIB-format ISCCP FD into matlab.
% Makes use of parse_ctl.m read_subvol_sequential.m.
%
% Adapted from QuikSCAT converter to handle
%  -template time conversions %y4%m2
%  -sequential reads handled by read_subvol_sequential.m
%
% SPdeS 10/23/06

fname='/home/alpha2/deszoeke/ISCCP/FD/MPF/flux.ctl'

[nlon,lon,nlat,lat,nlev,lev,ntime,time,dset,dtitle,...
 machineformat,undef,nvars,name,vlev,units,sequential]...
    =parse_ctl(fname);

% open monthly files for reading in the correct architecture
fid = fopen(dset,'r',machineformat);
nbyte=4;   % GrADS is single-precision by default
recbyte=4; % record length indicators are 4-byte integers

fivar=1:length(name); % read all

% initialize arrays
for m=fivar
  eval([name{m} '=NaN+zeros(nlon,nlat,ntime);'])
end

% read variables in order:
% (fastest varying) lon,lat,lev,variable,time (slowest varying)
for l=1:ntime
  % open monthly files for reading in the correct architecture
  [yr,mo] = datevec(time(l));
  datafile=[dset(1:end-6) sprintf('%04i%02i',[yr,mo])]; % handle particular date template
  fid = fopen(datafile,'r',machineformat);
  for m=1:length(name)
    vol=NaN+zeros(nlon,nlat,vlev(m));
    if ismember(m,fivar)
      vol(:)=read_subvol_sequential(fid,1,nlon,1,nlat,1,vlev(m),nlon,nlat,vlev(m),'single');
      vol(vol==undef)=NaN;
      eval([name{m} '(:,:,l)=squeeze(vol);']);
    else
      status=fseek(fid,(recbyte+nbyte*nlat*nlon+recbyte)*vlev(m),0); % skip over the 3D volume with reclen indicators
    end
  end
  fclose(fid);
end

% make a climatology
for mi=1:12
  for m=1:length(name)
    eval([name{m} 'm(:,:,mi)=nanmean(' name{m} '(:,:,mi:12:end),3);']);
  end
end
% save the climatology to disk
%save /home/alpha2/deszoeke/ISCCP/FD/isccpfd_clm1996-2003.mat lat lon srd1crm sru1crm sru5crm sxd1flm sxd5flm sxu1flm sxu5flm trd1crm tru1crm tru5crm txd1flm txu1flm txu5flm
