%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% generate 1/6th, MODIS, OSTIA figure
clf reset, orient tall, wysiwyg, colormap(cmap)
lat1=21; lat2=39; lon1=135; lon2=245; cx=[10 25];

axes('position',[.05 .78 .92 .18])
tmp=readbin('modisbin.dat',8515006); tmp(6)=1300;
lon=tmp(1):(tmp(2)-tmp(1))/(tmp(5)-1):tmp(2);
lat=tmp(3):(tmp(4)-tmp(3))/(tmp(6)-1):tmp(4);
sst=reshape(tmp(7:length(tmp)),length(lon),length(lat))-273.15;
ix=find(lon>=lon1&lon<=lon2); iy=find(lat>=lat1&lat<=lat2);
mypcolor(lon(ix),lat(iy),sst(ix,iy)'); caxis(cx), thincolorbar
title('a) 8-day-averaged SST from MODIS data, March 3 to March 10, 2007')
set(gca,'xtick',140:20:240,'ytick',15:10:35)
set(gca,'yticklabel',{'15N','25N','35N'},'xticklabel',[])

axes('position',[.05 .56 .92 .18])
lon=135:(245-135)/2200:245;
lat=15:(45-15)/600:45;
tmp=readbin('ostiabin.dat',1322807);
lon=tmp(1):(tmp(2)-tmp(1))/(tmp(5)-1):tmp(2);
lat=tmp(3):(tmp(4)-tmp(3))/(tmp(6)-1):tmp(4);
sst=reshape(tmp(7:length(tmp)),length(lon),length(lat))-273.15;
ix=find(lon>=lon1&lon<=lon2); iy=find(lat>=lat1&lat<=lat2);
mypcolor(lon(ix),lat(iy),sst(ix,iy)'); caxis(cx), thincolorbar
title('b) SST on March 3, 2007 from the OSTIA SST analysis')
set(gca,'xtick',140:20:240,'ytick',15:10:35)
set(gca,'yticklabel',{'15N','25N','35N'},'xticklabel',[])

axes('position',[.05 .34 .92 .18])
GRID_0625, lat=lat0625; lon=lon0625; nx=length(lon); ny=length(lat);
T=readbin('SST.0000307440.data',[nx ny]); T(find(~T))=nan;
ix=find(lon>=lon1&lon<=lon2); iy=find(lat>=lat1&lat<=lat2);
mypcolor(lon(ix),lat(iy),T(ix,iy)'); caxis(cx), thincolorbar
title('c) SST snapshot on March 3, 1993 from ECCO2 1/16^o global-ocean simulation')
set(gca,'xtick',140:20:240,'ytick',15:10:35)
set(gca,'yticklabel',{'15N','25N','35N'})
set(gca,'xticklabel',{'140E','160E','180E','160W','140W','120W'})

print -depsc sixteen_modis_ostia

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% generate 1/6th, AVHRR, MODIS, OSTIA figure
clf reset, orient tall, wysiwyg, colormap(cmap)
lat1=17; lat2=42; lon1=135; lon2=245; cx=[4 26];

axes('position',[.05 .78 .92 .18])
GRID_0625, lat=lat0625; lon=lon0625; nx=length(lon); ny=length(lat);
T=readbin('SST.0000307440.data',[nx ny]); T(find(~T))=nan;
ix=find(lon>=lon1&lon<=lon2); iy=find(lat>=lat1&lat<=lat2);
mypcolor(lon(ix),lat(iy),T(ix,iy)'); caxis(cx), thincolorbar
title('a) SST snapshot on March 3, 1993 from ECCO2 1/16^o global-ocean simulation')
set(gca,'xtick',140:20:240,'ytick',20:10:40)
set(gca,'yticklabel',{'20N','30N','40N'},'xticklabel',[])

axes('position',[.05 .56 .92 .18])
dx=360/2048; lon=(-180+dx/2):dx:180; lat=(-90+dx/2):dx:90;
lon(find(lon<0))=lon(find(lon<0))+360; [lon il]=sort(lon); 
land=flipud(readbin('landmask2048x1024.dat',[2048 1024],1,'uchar')');
sst=flipud(hdfread('199308h18ed-gdm.hdf','8-bit Raster Image #3'));
sst=sst(:,il); sst=0.15*double(sst)-3.0;
ix=find(lon>=lon1&lon<=lon2); iy=find(lat>=lat1&lat<=lat2);
mypcolor(lon(ix),lat(iy),sst(iy,ix)); caxis(cx), thincolorbar
title('b) 8-day-averaged SST from AVHRR data, February 26 to March 8, 1993')
set(gca,'xtick',140:20:240,'ytick',20:10:40)
set(gca,'yticklabel',{'20N','30N','40N'},'xticklabel',[])

axes('position',[.05 .34 .92 .18])
tmp=readbin('modisbin.dat',8515006); tmp(6)=1300;
lon=tmp(1):(tmp(2)-tmp(1))/(tmp(5)-1):tmp(2);
lat=tmp(3):(tmp(4)-tmp(3))/(tmp(6)-1):tmp(4);
sst=reshape(tmp(7:length(tmp)),length(lon),length(lat))-273.15;
ix=find(lon>=lon1&lon<=lon2); iy=find(lat>=lat1&lat<=lat2);
mypcolor(lon(ix),lat(iy),sst(ix,iy)'); caxis(cx), thincolorbar
title('c) 8-day-averaged SST from MODIS data, March 3 to March 10, 2007')
set(gca,'xtick',140:20:240,'ytick',20:10:40)
set(gca,'yticklabel',{'20N','30N','40N'},'xticklabel',[])

axes('position',[.05 .12 .92 .18])
lon=135:(245-135)/2200:245;
lat=15:(45-15)/600:45;
tmp=readbin('ostiabin.dat',1322807);
lon=tmp(1):(tmp(2)-tmp(1))/(tmp(5)-1):tmp(2);
lat=tmp(3):(tmp(4)-tmp(3))/(tmp(6)-1):tmp(4);
sst=reshape(tmp(7:length(tmp)),length(lon),length(lat))-273.15;
ix=find(lon>=lon1&lon<=lon2); iy=find(lat>=lat1&lat<=lat2);
mypcolor(lon(ix),lat(iy),sst(ix,iy)'); caxis(cx), thincolorbar
title('d) SST on March 3, 2007 from the OSTIA SST analysis')
set(gca,'xtick',140:20:240,'ytick',20:10:40)
set(gca,'yticklabel',{'20N','30N','40N'})
set(gca,'xticklabel',{'140E','160E','180E','160W','140W','120W'})

print -depsc modis

%%%%%%%%%%%%%%%%%%%%%
% read AVHRR data
lat1=17; lat2=42; lon1=135; lon2=245; cx=[4 26];

dx=360/2048; lon=(-180+dx/2):dx:180; lat=(-90+dx/2):dx:90;
lon(find(lon<0))=lon(find(lon<0))+360; [lon il]=sort(lon); 
land=flipud(readbin('landmask2048x1024.dat',[2048 1024],1,'uchar')');
sst=flipud(hdfread('199308h18ed-gdm.hdf','8-bit Raster Image #3'));
sst=sst(:,il); sst=0.15*double(sst)-3.0;
clf reset, colormap(cmap); axes('position',[.05 .53 .92 .42])
ix=find(lon>=lon1&lon<=lon2); iy=find(lat>=lat1&lat<=lat2);
mypcolor(lon(ix),lat(iy),sst(iy,ix)); caxis(cx), thincolorbar
title('SST on March 2, 1993 from AVHRR data')
set(gca,'xtick',[],'ytick',20:10:40,'yticklabel',{'20N','30N','40N'})

GRID_0625, lat=lat0625; lon=lon0625; nx=length(lon); ny=length(lat);
T=readbin('SST.0000307440.data',[nx ny]); T(find(~T))=nan;
ix=find(lon>=lon1&lon<=lon2); iy=find(lat>=lat1&lat<=lat2);
axes('position',[.05 .05 .92 .42])
mypcolor(lon(ix),lat(iy),T(ix,iy)'); caxis(cx), thincolorbar
title('SST on March 2, 1993 from ECCO2 1/16^o global-ocean simulation')
set(gca,'xtick',140:20:240,'ytick',20:10:40)
set(gca,'yticklabel',{'20N','30N','40N'})
set(gca,'xticklabel',{'140E','160E','180E','160W','140W','120W'})

print -depsc AVHRR_ECCO2_SST

%%%%%%%%%%%%%%%%%%%%%
% Liming's figure
load Std_data_new
clf reset, colormap(cmap)
cx1=[0 6]; cx2=[0 .8]; iy=12:169
longitude(find(longitude<0))=longitude(find(longitude<0))+360;
[longitude ix]=sort(longitude);

subplot(221)
mypcolor(longitude,latitude(iy),std_ori_AMSR(iy,ix));
caxis(cx1), thincolorbar
title('a) Std of Reynolds hi-res SST')
set(gca,'xtick',[0 180 360],'xticklabel',{'0E','180E','360E'})
set(gca,'ytick',[-60 0 60],'yticklabel',{'60S','0N','60N'})

subplot(222)
mypcolor(longitude,latitude(iy),std_ori_CUBE(iy,ix));
caxis(cx1), thincolorbar
title('b) Std of ECCO2 solution')
set(gca,'xtick',[0 180 360],'xticklabel',{'0E','180E','360E'})
set(gca,'ytick',[-60 0 60],'yticklabel',{'60S','0N','60N'})

subplot(223)
mypcolor(longitude,latitude(iy),std_less1mon_AMSR(iy,ix));
caxis(cx2), thincolorbar
title('c) Std of Reynolds, periods < 1 mo')
set(gca,'xtick',[0 180 360],'xticklabel',{'0E','180E','360E'})
set(gca,'ytick',[-60 0 60],'yticklabel',{'60S','0N','60N'})

subplot(224)
mypcolor(longitude,latitude(iy),std_less1mon_CUBE(iy,ix));
caxis(cx2), thincolorbar
title('d) Std of ECCO2, periods < 1 mo')
set(gca,'xtick',[0 180 360],'xticklabel',{'0E','180E','360E'})
set(gca,'ytick',[-60 0 60],'yticklabel',{'60S','0N','60N'})

print -depsc stdsst