% Read RSI's optimally interpolated TMI SST. This requires the simple
% read_rss_oi_sst.m, which I have put in the directory added to the path:
addpath('/home/sol/deszoeke/data/oi_sst/read_routines/matlab/');

% The daily files September-October 2001 have been downloaded and
% unzipped in:
directory = '/home/sol/deszoeke/data/oi_sst/tmi/2001/';

% Simon de Szoeke
% IPRC
% 17 September 2004
% Modified 8 October 2001 to read in all of 2001, and go from -21 31N lat

prefix = 'tmi.fusion.';
suffix = '.v01';
year = 2001;

% The center of the first cell of the 1440 column and 720 row map is at 0.125 
% E longitude and -89.875 latitude. The center of the second cell is 0.375 E
% longitude, -89.875 latitude.
lonc =   0.125:0.25:359.875; % center of the patches
latc = -89.875:0.25:89.875;
%lonw = lonc-0.125; % corner of the patches
%lats = latc-0.125;
lonrange = 360+[-181 -20];
latrange = [-31 31];

lon = lonc(lonc>=lonrange(1) & lonc<=lonrange(2));
lat = latc(latc>=latrange(1) & latc<=latrange(2));
[Y,X] = meshgrid(lat,lon);

temp = NaN+zeros(length(lonc),length(latc));
% The yeardays for September-October are 244-304.
% yday = 244:304;
% Use all days in 2001
yday = 1:365;

% setup the file for writing
missingv = -1e33;
nc = netcdf('sst_fill_2001_big.nc','clobber'); % don't clobber!
nc.description = 'OI "daily" TMI SST from RSS, interpolated over land';
nc('time') = 0;
nc('lon') = length(lon);
nc('lat') = length(lat);

nc{'time'} = 'time';
nc{'time'}.units = 'days since December 31, 2000';
nc{'lon'} = 'lon';
nc{'lon'}.units = 'degrees east longitude';
nc{'lat'} = 'lat';
nc{'lat'}.units = 'degrees north longitude';

nc{'sst'} = {'time','lon','lat'};
nc{'sst'}.units = 'degrees C';
nc{'sst'}.missing_value = missingv;
nc{'available_mask'} = {'time','lon','lat'};
nc{'available_mask'}.description = '1 for a datum from OI_RSS, 0 otherwise';

nc{'lon'}(:) = lon;
nc{'lat'}(:) = lat;
% loop over and load each day from its file
for di = 1:length(yday)
  filename = [directory,prefix, num2str(year,'%4.4i'),'.',num2str(yday(di),'%3.3i'),suffix];
  fprintf(1,['Reading ' filename '\n']);
  temp = read_rss_oi_sst(filename);
  t = temp(lonc>=lonrange(1) & lonc<=lonrange(2), latc>=latrange(1) & latc<=latrange(2));
  % fill in gaps near continents, islands, etc.
  mask = t<200;
  t(~mask) = griddata(X(mask),Y(mask),t(mask),X(~mask),Y(~mask),'cubic');
  if any(any(isnan(t))); t(isnan(t)) = missingv; end
  nc{'sst'}(di,:,:) = t;
  nc{'available_mask'}(di,:,:) = 1*mask;
end

nc = close(nc);

if 0
% sample read
nc = netcdf('sst_fill.nc');
lon = nc{'lon'}(:);
lat = nc{'lat'}(:);
sst = nc{'sst'}(:);
sst(sst==nc{'sst'}.missing_value(:)) = NaN;
am = nc{'available_mask'}(:); % FYI the same every day

M = moviein(size(sst,1));
for j = 1:size(sst,1)
  [c,h] = contourf(360-lon,lat,squeeze(sst(j,:,:))',15:0.5:33);
  set(h(:),'linestyle','none');
  title(['SST (C), ' datestr(j+datenum('31-aug-2001'))]);
  colorbar
  set(gca,'xdir','reverse');
  xlabel('degrees W longitude');
  ylabel('degrees N latitude');
  M(:,j) = getframe
end
end