% Compare two TMI products from RSS:
% The weekly CLW and SST during May-Oct 2001
%
% 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 matlab day
start = datenum('1-may-2001');
endt  = datenum('31-oct-2001');
dt = 7;
mdayp = start - mod(start-datenum('3-dec-1997'),dt); % previous Weds
indp = (mdayp-datenum('3-dec-1997'))/dt + 1; % c(0:n-1) index of first Weds in time interval
mdayn = endt - mod(endt-datenum('3-dec-1997'),dt);
indn = (mdayn-datenum('3-dec-1997'))/dt; % index of the last Weds in interval
sm   = zeros(1440,320);
sst1  = sm;
clw1  = sm;
sst2  = sm;
clw2  = sm;
vari1 = sm;

% loop over desired weeks (reading is slow)
% sum one file at a time to save memory
for ind = indp:indn
  sst = loaddods(['http://apdrc.soest.hawaii.edu:9090/dods/public_data/satellite_product/TRMM/trmm_weekly?sst',sprintf('[%i]',ind),'[0:0][0:319][0:1439]'])';
  sst(sst<0) = NaN;
  clw = loaddods(['http://apdrc.soest.hawaii.edu:9090/dods/public_data/satellite_product/TRMM/trmm_weekly?clw',sprintf('[%i]',ind),'[0:0][0:319][0:1439]'])';
  clw(clw<0) = NaN;
  
  mask = isfinite(sst) & isfinite(clw);
  sst1(mask)  = sst1(mask) + sst(mask);
  clw1(mask)  = clw1(mask) + clw(mask);
  sst2(mask)  = sst2(mask) + sst(mask).^2;
  clw2(mask)  = clw2(mask) + clw(mask).^2;
  vari1(mask) = vari1(mask) + sst(mask).*clw(mask);
  sm = sm + mask;
end

sst  = sst1./sm;
clw  = clw1./sm;
sstv = sst2./sm - sst.^2;
clwv = clw2./sm - clw.^2;
vari = vari1./sm - sst.*clw;
corr = vari./sqrt(sstv.*clwv);

mask = sm>=10; % must have 10 weeks of data
vari(~mask) = NaN;
corr(~mask) = NaN;

save TMI_sst+clw.mat sst clw vari corr sm

% plots
% mean daily and weekly
subplot(3,1,1)
pcolor(lon,lat,sst');shading flat;
colorbar;
title({'weekly TMI May-October 2001 mean' 'SST (C)'})
subplot(3,1,2)
pcolor(lon,lat,clw');shading flat;
colorbar;
title('cloud liquid water (mm)')
ylabel('degrees latitude')
subplot(3,1,3)
pcolor(lon,lat,corr');shading flat;
colorbar;
title('correlation')
xlabel('degrees E longitude')
orient tall

% zoom in on E. Pacific
flonEP = find(lon>=230 & lon<=280);
flatEP = find(lat>=-15 & lat<=15);
subplot(3,1,1)
pcolor(lon(flonEP),lat(flatEP),sst(flonEP,flatEP)');shading flat;
colorbar;
title({'weekly TMI May-October 2001 mean' 'SST (C)'})
subplot(3,1,2)
pcolor(lon(flonEP),lat(flatEP),clw(flonEP,flatEP)');shading flat;
colorbar;
title('cloud liquid water (CLW, mm)')
ylabel('degrees latitude')
subplot(3,1,3)
pcolor(lon(flonEP),lat(flatEP),sqrt(sstv(flonEP,flatEP))'.*corr(flonEP,flatEP)');shading flat;
colorbar;
title('SST regressed on CLW (K)')
xlabel('degrees E longitude')
orient tall

print -depsc TMI_SST+CLW.eps