clear

casenum=input('Case number: ');

nlat = 142;
nlon = 362;
nlev = 20; % control file is wrong; 20 not 30 levels.

mon = char('jan','feb','mar','apr','may','jun','jul','aug','sep','oct','nov','dec');
nday = [31 28 31 30 31 30 31 31 30 31 30 31];
year = [1999 2000 2001 2002];


% allocate data array(s)
flatsz = [nlon nlat]; % for variables with no vertical dimension
sz = [nlev nlon nlat];
ur = NaN+zeros(sz);
vr = NaN+zeros(sz);
tr = vr;
kv= vr;
hflux=NaN+zeros(flatsz);
sflux=hflux;

datadir = sprintf('/disk/a2/iroam/DATA%d/',casenum);

nt=0;


for iyear=1:length(year)
for imon=1:12
  filename = [datadir sprintf('oceandaily.%s.%d.dat',mon(imon,:),year(iyear))];
filename
% open the appointed files
fid = fopen(filename,'r','ieee-be');

clear v

for nt=1:nday(imon)

% READ one 3-D volume of data (read prod(sz) floating point reals)
for i = 1:nlev
  recl = fread(fid,1,'uint32');
  ur(i,:) = fread(fid,prod(flatsz),'float');
  recl2 = fread(fid,1,'uint32');
  if recl2~=recl; disp('Unexpected record length.'); break; end
end

for i = 1:nlev
  recl = fread(fid,1,'uint32');
  vr(i,:) = fread(fid,prod(flatsz),'float');
  recl2 = fread(fid,1,'uint32');
  if recl2~=recl; disp('Unexpected record length.'); break; end
end
  v(nt,:,:,:)=vr(:,:,:);
for i = 1:nlev
  recl = fread(fid,1,'uint32');
  tr(i,:) = fread(fid,prod(flatsz),'float');
  recl2 = fread(fid,1,'uint32');
  if recl2~=recl; disp('Unexpected record length.'); break; end
end

  recl = fread(fid,1,'uint32');
  hflux = fread(fid,prod(flatsz),'float');
  recl2 = fread(fid,1,'uint32');
  if recl2~=recl; disp('Unexpected record length.'); break; end
  recl = fread(fid,1,'uint32');
  sflux = fread(fid,prod(flatsz),'float');
  recl2 = fread(fid,1,'uint32');
  if recl2~=recl; disp('Unexpected record length.'); break; end
for i = 1:nlev
  recl = fread(fid,1,'uint32');
  kv(i,:) = fread(fid,prod(flatsz),'float');
  recl2 = fread(fid,1,'uint32');
  if recl2~=recl; disp('Unexpected record length.'); break; end
end
end



fclose(fid);

v(v<-1e10) = NaN;

filename=['/export/t2d2/kelvin/iROAM/daily_output/case' int2str(casenum) '/daily_v_' mon(imon,:) '_' int2str(year(iyear))];
eval(['save ' filename ' v'])

end
end
