% 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/sepoct2001/';

% Simon de Szoeke
% IPRC
% 17 September 2004

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;
temp = NaN+zeros(length(lonc),length(latc));

% The yeardays for September-October are 244-304.
yday = 244:304;

% this array is megabig, and runs sol out of memory
%sst = zeros(length(lonc),sum(abs(latc)<40),length(yday));

% 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];
  temp = read_rss_oi_sst(filename);
  % cut out areas north and south of 40degrees for TMI
  t = temp(:,abs(latc)<40);
  lon = lonc;
  lat = latc(abs(latc)<40);
  if length(t) > 101 % do only once
    switch upper('EPIC')
     case 'EPIC'
      t = t(lon>=360-129 & lon<=360-64, lat>=-14 & lat<=21);
      lon = lon(lon>=360-129 & lon<=360-64);
      lat = lat(lat>=-14 & lat<=21);
     case 'HAWAII'
      % zoom in on Hawaiian islands to test
      % gridding & filling on a smaller domain
      t = t(750:850,220:260);
      lon = lon(750:850);
      lat = lat(220:260);
    end
  end
  % fill in gaps near continents, islands, etc.
  mask = t<200;
  [Y,X] = meshgrid(lat,lon);
  %tf = NaN+zeros(size(t));
  %tf(~mask) = griddata(X(mask),Y(mask),t(mask),X(~mask),Y(~mask),'cubic');
  %t(~mask) = tf(~mask)
  t(~mask) = griddata(X(mask),Y(mask),t(mask),X(~mask),Y(~mask),'cubic');
  save sst_fill.mat t -append
end

% plot 1 day
subplot(3,1,1)
pcolor(lon,lat,t');shading flat
axis equal;axis tight
title 'SST'
subplot(3,1,2)
pcolor(lon,lat,tf');shading flat
axis equal;axis tight
ylabel 'latitude'
title 'filled'
subplot(3,1,3)
pcolor(lon,lat,t');shading flat
hold on
pcolor(lon,lat,tf');shading flat
axis equal;axis tight
title 'both'
xlabel 'longitude'
orient tall
print -depsc2 SST_fill.eps