clear

load basin_mask60S60N

load VertShear_era_clim
load AbsVort_era_clim
load RelHum_era_clim

absvcl_era = absvcl;
vshearcl_era = vshearcl;
rhcl_era1 = rhcl;

load FP_dens_clim
fpdens_cl = fpdens_cl;

load AbsVort_clim
load RelHum_clim
load VertShear_clim

rhcl1 = rhcl;

load SST_minus_trop_clim

%new RH
load RH_clim_SSMI_ERA40_mod.mat
rhcl_era = rh_clim*100;

%new RH
load RH_clim_SSMI_NCEP_mod.mat
rhcl = rh_clim*100;

load ERA40_Vmax_clim_61_00.mat
Vmax_cl_era=Vmax_cl;

load NCEP_Vmax_clim_61_00.mat

[ny,nx,nt]=size(fpdens_cl);

YY = sst;
for iy=1:length(Y)
  YY(iy,:,:) = Y(iy);
end
area=cos(YY/180*pi);

mask = ~isnan(rhcl_era+rh_clim+sst);
iok=~isnan(rhcl_era(:)+rh_clim(:)+sst(:));

y = fpdens_cl(iok);

n=length(y);
avmax=3.7

x=[min(abs(1e5*absvcl(iok)),avmax),rhcl(iok),sst(iok),vshearcl(iok)];
x_era=[min(abs(1e5*absvcl_era(iok)),avmax),rhcl_era(iok),sst(iok),vshearcl_era(iok)];

x2=[x;x_era];
y2=[y;y];

[b,dev0,stats] = glmfit(x2,y2,'poisson','offset',log([area(iok);area(iok)]));

wrk = glmval(b,x,'log','offset',log(area(iok)));
ypoisson=zeros(ny,nx,nt)+NaN;
ypoisson(mask) = wrk;
ypoisson_ncep=ypoisson;

wrk = glmval(b,x_era,'log','offset',log(area(iok)));
ypoisson=zeros(ny,nx,nt)+NaN;
ypoisson(mask) = wrk;
ypoisson_era=ypoisson;

data(:,:,:,1) =  fpdens_cl+0*sst;
data(:,:,:,2) =  ypoisson_ncep;
data(:,:,:,3) =  ypoisson_era;

leg_titles={'Obs','NCEP','ERA'};

dtas = {[1,2,3]};

for id = 1:length(dtas)
  c1 = [1,1,1]*0;
  c3 = [176 226 255]/255;
  c2 = [1,1,1];
  figure(1)
  dta = dtas{id}

  subplot(5,2,1)
  indy=Y>=0;
  h = bar(squeeze(nansum(nansum(data(indy,:,:,dta)))));
  set(h(1),'facecolor',c1)
  set(h(2),'facecolor',c2)
  set(h(3),'facecolor',c3)
  set(gca,'xtick',[1:12],'xticklabel',['J','F','M','A','M','J','J','A','S','O','N','D']')
  title('NH')
  axis('tight')
  
  subplot(5,2,2)
  indy=Y<=0;
  h = bar(squeeze(nansum(nansum(data(indy,:,:,dta)))));
  set(h(1),'facecolor',c1)
  set(h(2),'facecolor',c2)
  set(h(3),'facecolor',c3)
  set(gca,'xtick',[1:12],'xticklabel',['J','F','M','A','M','J','J','A','S','O','N','D']')
  title('SH')
  axis('tight')
  
  titles={'South Indian','Australian','South Pacific','North Indian','Western North Pacific','Central North Pacific','Eastern North Pacific','Atlantic'};
  
  for ii=1:8
    subplot(5,2,ii+2)
    wrk = basin_mask==ii;
    wrk = area.*wrk(:,:,ones(1,nt));
    wrk = wrk(:,:,:,ones(1,length(dta)));
    h = bar(squeeze(nansum(nansum(wrk.*data(:,:,:,dta)))));
    set(h(1),'facecolor',c1)
    set(h(2),'facecolor',c2)
    set(h(3),'facecolor',c3);
    set(gca,'xtick',[1:12],'xticklabel',['J','F','M','A','M','J','J','A','S','O','N','D']')
    
    title(titles{ii})
    axis('tight')
  end
  subplot(5,2,1)
  h=legend(leg_titles{dta},'location','Northwest','Orientation','vertical');
  set(h,'fontsize',6)
  
  %fn = ['basin_climatology',int2str(id)]
  %exportfig(gcf,fn,'width',6,'height',7,'color','cmyk','FontMode','scaled','fontsize',0.8,'linemode','scaled')
  pause
end

letters = {'a','b','c','d','e','f','g','h','i','j','k','l'};
figure(2)
colormap('default');
c=colormap;
c(1,:)=[1,1,1]*0.65;
colormap(c)

dtas={[1,2],[1,3]}
for ip = 1:2

ilet =1;
  dta = dtas{ip};
  for id =1:length(dta)
    
    subplot(length(dta),2,2*id-1)
    my_pcolor(X,Y,sum(data(:,:,[1:3],dta(id)),3))
    title(['(',letters{ilet},')',' JFM ',leg_titles{dta(id)}])
    ilet = ilet + 1;
    caxis([0,0.25])
    pacific_noresize_at_all
    set(gca,'ylim',[-30,40])
    
    subplot(length(dta),2,2*id)
    my_pcolor(X,Y,sum(data(:,:,[8:10],dta(id)),3))
    title(['(',letters{ilet},')',' ASO ',leg_titles{dta(id)}])
    ilet = ilet + 1;
    caxis([0,0.25])
    pacific_noresize_at_all
    set(gca,'ylim',[-30,40])
    
    
  end
  
  %fn = ['climatology_map',int2str(ip)];;
  %exportfig(gcf,fn,'width',5.5,'height',7,'color','cmyk','FontMode','fixed','fontsize',8,'linemode','scaled','renderer','painters')
  
end
