clear all

% cd('/Users/skrueger/ RESEARCH/ MEETINGS/2008 Meetings/POST/NCEP postgfs')
fn = 'postgfs.2008052212'
fn = 'postgfs.2008061412'

% [FID, MESSAGE] = FOPEN(FILENAME,PERMISSION,MACHINEFORMAT)
[fid, MESSAGE] = fopen(fn,'r','ieee-be')    % specify correct machine format

% add 4 extra bytes at start and end of each record
% add 130*4 extra bytes before loop over time & points

    dummy = fread(fid,1, 'integer*4');

hour  = fread(fid,1, 'integer*4') 
month = fread(fid,1, 'integer*4') 
day   = fread(fid,1, 'integer*4')
year  = fread(fid,1, 'integer*4')
nsfc  = fread(fid,1, 'integer*4')
nflx  = fread(fid,1, 'integer*4')
nvar  = fread(fid,1, 'integer*4')
levs = fread(fid,1, 'integer*4')
npoint= fread(fid,1, 'integer*4')
starting= fread(fid,1, 'integer*4')
ending  = fread(fid,1, 'integer*4')
interval= fread(fid,1, 'integer*4')

    dummy = fread(fid,1, 'integer*4');
    dummy = fread(fid,1, 'integer*4');

sigi = fread(fid,levs+1, 'real*4'); 
sigl = fread(fid,levs, 'real*4'); 

    dummy = fread(fid,1, 'integer*4');

extra = fread(fid,130, 'real*4'); % this works; checked results

ntime=(ending-starting)/interval+1

for i=1:ntime
for j=1:npoint

    dummy = fread(fid,1, 'integer*4');
sfc = fread(fid,nsfc, 'real*4');
    dummy = fread(fid,1, 'integer*4');
    if i == 1; disp( [sfc(1) sfc(2)] ); end

    dummy = fread(fid,1, 'integer*4');
flx = fread(fid,nflx, 'real*4');
    dummy = fread(fid,1, 'integer*4');


% ! record 5 ... levs of model zonal wind velocity (m/s)
% ! record 6 ... levs of model meridional wind velocity (m/s)
% ! record 7 ... levs of model temperature (K)
% ! record 8 ... levs of model specific humidity (g/g)
% ! record 9 ... levs of model pressure (mPa)

    dummy = fread(fid,1, 'integer*4');
u_wind = fread(fid,levs, 'real*4');
    dummy = fread(fid,1, 'integer*4');    
    dummy = fread(fid,1, 'integer*4');
v_wind = fread(fid,levs, 'real*4');
    dummy = fread(fid,1, 'integer*4'); 
    dummy = fread(fid,1, 'integer*4');
tl = fread(fid,levs, 'real*4');
    dummy = fread(fid,1, 'integer*4');
    dummy = fread(fid,1, 'integer*4');
ql = fread(fid,levs, 'real*4');
    dummy = fread(fid,1, 'integer*4');
    dummy = fread(fid,1, 'integer*4');
pl = fread(fid,levs, 'real*4');
    dummy = fread(fid,1, 'integer*4');

% ! record 10 ... levs of model derived omega (mPa/sec)
% !             This is an inadvertant unit selection due to making two unit
% !             conversions. mPa is milli-Pascal.
    dummy = fread(fid,1, 'integer*4');
omega = fread(fid,levs, 'real*4');
    dummy = fread(fid,1, 'integer*4');
% ! record 11 ... levs of model derived dtdt (advection) (K/sec)
% ! record 12 ... levs of model derived dqdt (advection) (g/g/sec)
       
    dummy = fread(fid,1, 'integer*4');
dtdt = fread(fid,levs, 'real*4');
    dummy = fread(fid,1, 'integer*4');
    dummy = fread(fid,1, 'integer*4');
dqdt = fread(fid,levs, 'real*4');
    dummy = fread(fid,1, 'integer*4');

% ! record 13 ... levs of model cloud water/ice concentration (g/g)
% ! record 14 ... levs of model derived cloud water tendency (advection) (g/g/sec)
% ! record 15 ... levs of model cloud fraction 

    dummy = fread(fid,1, 'integer*4');
qci = fread(fid,levs, 'real*4');
    dummy = fread(fid,1, 'integer*4');
    dummy = fread(fid,1, 'integer*4');
dqci = fread(fid,levs, 'real*4');
    dummy = fread(fid,1, 'integer*4');
    dummy = fread(fid,1, 'integer*4');
cld = fread(fid,levs, 'real*4');
    dummy = fread(fid,1, 'integer*4');

end
end

dummy130 = fread(fid,inf, 'integer*4');