% Compare two TMI SST products from RSS:
% The new (2004) optimally interpolated "daily" data and
% the standard weekly SST data
%
% Simon de Szoeke
% 2004-2005

lonc =   0.125:0.25:359.875; % center of the patches
latc = -89.875:0.25:89.875;
lonrange = [0 360];
latrange = [-40 40];
lon = lonc(lonc>=lonrange(1) & lonc<=lonrange(2));
lat = latc(latc>=latrange(1) & latc<=latrange(2));
scale_factor = 0.15;
add_offset   = -3;

% start and end of averaging interval in "Julian" yearday
start = datenum('1-may-2001')-datenum('0-jan-2001'); % no offset for synch to weeks
endt  = datenum('31-oct-2001')-datenum('0-jan-2001');
dt = 7;
sm   = zeros(1440,320);
daily1  = sm;
weekly1 = sm;
vari1 = sm;

% Start: read the previous and next weekly SST for Wednesdays surrounding
% the first day, start.
mday = start+datenum('0-jan-2001');
mdayp = mday - mod(mday-datenum('3-dec-1997'),dt); % previous Weds
indp = (mdayp-datenum('3-dec-1997'))/dt; % c(0:n-1) index
% read the current week sst (indp)
sstp = loaddods(['http://apdrc.soest.hawaii.edu:9090/dods/public_data/satellite_product/TRMM/trmm_weekly?sst',sprintf('[%i]',indp),'[0:0][0:319][0:1439]'])';
sstp(sstp<0) = NaN;
% read the next week sst (indp+1)
sstn = loaddods(['http://apdrc.soest.hawaii.edu:9090/dods/public_data/satellite_product/TRMM/trmm_weekly?sst',sprintf('[%i]',indp+1),'[0:0][0:319][0:1439]'])';
sstn(sstn<0) = NaN;
% loop over desired days for jday
for jday = start:endt
  mday  = jday+datenum('0-jan-2001');
  mdayp = mday - mod(mday-datenum('3-dec-1997'),dt);
  %disp(sprintf('jday=%3d, days since last Weds=%1d',jday,mod(mday-datenum('3-dec-1997'),dt)));
  if mday==mdayp % mday is a Wednesday, increment weekly data
    indp = (mdayp-datenum('3-dec-1997'))/dt;
    sstp = sstn; % first timestep sstn undefined!
    % read the next week sst (indp+1)
    sstn = loaddods(['http://apdrc.soest.hawaii.edu:9090/dods/public_data/satellite_product/TRMM/trmm_weekly?sst',sprintf('[%i]',indp+1),'[0:0][0:319][0:1439]'])';
    sstn(sstn<0) = NaN;
  end
  % read the daily sst -> t
  fname = ['/home/sol/deszoeke/data/oi_sst/tmi/2001/tmi.fusion.2001.',sprintf('%03d',jday),'.v01'];
  fid = fopen(fname);
  temp = reshape(fread(fid,1440*720,'char'),1440,720);
  fid = fclose(fid);
  temp(temp>250) = NaN;
  t = temp(:,latc>=latrange(1) & latc<=latrange(2))*scale_factor+add_offset;
  
  ti = sstp + (sstn-sstp)*mod(mday-datenum('3-dec-1997'),dt)/dt;
  mask = isfinite(t) & isfinite(ti);
  vari1(mask) = vari1(mask) + (t(mask)-ti(mask)).^2;
  daily1(mask)  =  daily1(mask) + t(mask); % mean
  weekly1(mask) = weekly1(mask) + ti(mask);
  sm = sm + mask;
end

daily  =  daily1./sm;
weekly = weekly1./sm;
bias = daily-weekly;
vari = vari1./sm - bias.^2;

mask = sm>(2*dt);
vari(~mask) = NaN;
bias(~mask) = NaN;

save vari_TMI.mat daily weekly vari sm

% plots
% mean daily and weekly
subplot(2,1,1)
pcolor(lon,lat,daily');shading flat;
colorbar;
title({'TMI' 'daily (C)'})
ylabel('degrees latitude')
subplot(2,1,2)
pcolor(lon,lat,weekly');shading flat;
colorbar;
title('weekly (C)')
xlabel('degrees E longitude')
% don't see much 

% look at the differences
subplot(2,1,1)
pcolor(lon,lat,bias');shading flat;
caxis([-1 1]); colorbar;
title({'TMI' 'SST_{daily} - SST_{weekly} bias (K)'})
ylabel('degrees latitude')
subplot(2,1,2)
pcolor(lon,lat,vari');shading flat;
caxis([0 1.5]); colorbar;
title('RMS difference (K)')
xlabel('degrees E longitude')

orient landscape
print -depsc tmi_daily_vs_weekly_sst.eps
