clear all;
close all;



nemodir1='/home/datawork-lops-siam-moawi/SIAMTEAM/xcouvela/NEMOWAVE/ORCA025/ORCA025_ERAIM/2014/';
%nemodir1='/home1/scratch/xcouvela/ORCA025_WW3_CPL/NCPL/';
%nemodir2='/home/datawork-lops-siam-moawi/SIAMTEAM/xcouvela/NEMOWAVE/ORCA025/ORCA025_ERAIM_CPL_NEW_BREIV16/AXELL03/';
nemodir2='/home/datawork-lops-siam-moawi/SIAMTEAM/xcouvela/NEMOWAVE/ORCA025/ORCA025_ERAIM_CPL_NEW_TAW/TAW_COR/';


nav_lat=ncread([nemodir1,'ORCA025_1d_20140101_20141231_SBC.nc'],'nav_lat');
nav_lon=ncread([nemodir1,'ORCA025_1d_20140101_20141231_SBC.nc'],'nav_lon');

%
%Each days: too long.... 
%
%wind=ncread([nemodir1,'ORCA025_3h_20140101_20141230_SBC.nc'],'windsp',[1 1 1],[inf inf inf],[1 1 8]);
%siconc=ncread([nemodir1,'ORCA025_3h_20140101_20141230_icemod.nc'],'siconc',[1 1 1],[inf inf inf],[1 1 8]);
%Cdncpl=ncread([nemodir1,'ORCA025_3h_20140101_20141230_SBC.nc'],'Cd',[1 1 1],[inf inf inf],[1 1 8]);
%Cdcpl=ncread([nemodir2,'ORCA025_3h_20140101_20141230_SBC.nc'],'Cd',[1 1 1],[inf inf inf],[1 1 8]);

%
%read and avg (daily)
%
all=1;
if all

t=1;
%for i=1:90
%disp(num2str(i))
%wind(:,:,i)  = squeeze(mean(ncread([nemodir1,'ORCA025_1d_20140101_20141231_SBC.nc'],...
%                                                      'windsp',[1 1 i],[inf inf 1]),3));
%siconc(:,:,i)= squeeze(mean(ncread([nemodir2,'ORCA025_1d_20140101_20141231_icemod.nc'],...
%                                                         'siconc',[1 1 i],[inf inf 1]),3));
%Cdncpl(:,:,i)= squeeze(mean(ncread([nemodir1,'ORCA025_1d_20140101_20141231_SBC.nc'],...
%                                                          'Cd',[1 1 i],[inf inf 1]),3));
%Cdcpl(:,:,i) = squeeze(mean(ncread([nemodir2,'ORCA025_1d_20140101_20141231_SBC.nc'],...
%                                                          'Cd',[1 1 i],[inf inf 1]),3));
%end

%for i=1:90
jj=1;
%for i=182:275
for i=1:365
disp(num2str(i))
wind(:,:,jj)  = squeeze(mean(ncread([nemodir1,'ORCA025_1d_20140101_20141231_SBC.nc'],...
                                                      'windsp',[1 224 i],[inf 201 1]),3));
siconc(:,:,jj)= squeeze(mean(ncread([nemodir2,'ORCA025_1d_20140101_20141231_icemod.nc'],...
                                                         'siconc',[1 224 i],[inf 201 1]),3));
Cdncpl(:,:,jj)= squeeze(mean(ncread([nemodir1,'ORCA025_1d_20140101_20141231_SBC.nc'],...
                                                          'taum',[1 224 i],[inf 201 1]),3));
Cdcpl(:,:,jj) = squeeze(mean(ncread([nemodir2,'ORCA025_1d_20140101_20141231_SBC.nc'],...
                                                          'taum',[1 224 i],[inf 201 1]),3));
jj=jj+1;
end

%read by box
%
else
lonmin=-70;lonmax=-50;
latmin=30;latmax=45;

imin=near(nav_lon(:,1),lonmin);
imax=near(nav_lon(:,1),lonmax);
jmin=near(nav_lat(1,:),latmin);
jmax=near(nav_lat(1,:),latmax);

wind  = ncread([nemodir1,'ORCA025_3h_20140101_20141231_SBC.nc'],...
                                                      'windsp',[imin jmin 1],[imax-imin+1 jmax-jmin+1 inf]);
siconc= ncread([nemodir1,'ORCA025_3h_20140101_20141231_icemod.nc'],...
                                                         'siconc',[imin jmin 1],[imax-imin+1 jmax-jmin+1 inf]);
Cdncpl= ncread([nemodir1,'ORCA025_3h_20140101_20141231_SBC.nc'],...
                                                          'Cd',[imin jmin 1],[imax-imin+1 jmax-jmin+1 inf]);
Cdcpl = ncread([nemodir2,'ORCA025_3h_20140101_20141231_SBC.nc'],...
                                                          'Cd',[imin jmin 1],[imax-imin+1 jmax-jmin+1 inf]);
end

maskice=siconc;
clear siconc
maskice(maskice~=0)=NaN;
maskice(~isnan(maskice))=1;

mask=wind;
mask(mask==0)=NaN;
mask(~isnan(mask))=1;

mask=mask.*maskice;

clear maskice

Cdncpl=Cdncpl.*mask;
Cdcpl=Cdcpl.*mask;
wind=wind.*mask;
Cdncpl(find(isnan(Cdncpl)))=[];
Cdcpl(find(isnan(Cdcpl)))=[];
wind(find(isnan(wind)))=[];

wind=reshape(wind,1,size(wind,1)*size(wind,2)*size(wind,3));
Cdncpl=reshape(Cdncpl,1,size(Cdncpl,1)*size(Cdncpl,2)*size(Cdncpl,3));
Cdcpl=reshape(Cdcpl,1,size(Cdcpl,1)*size(Cdcpl,2)*size(Cdcpl,3));

nbins=100;
h1=histogram2(wind,Cdncpl,nbins,'DisplayStyle','tile'); 
caxis([0 6e5])
colorbar;
figure;
h2=histogram2(wind,Cdcpl,nbins,'DisplayStyle','tile');
caxis([0 6e5])
colorbar

hncpl=h1.Values;
hcpl=h2.Values;

xncpl=h1.XBinEdges;
yncpl=h1.YBinEdges;

xcpl=h2.XBinEdges;
ycpl=h2.YBinEdges;

hncpl(hncpl<100)=NaN;
hcpl(hcpl<100)=NaN;
hncpl=hncpl';
hcpl=hcpl';

%
%NEW FIGURE GET THE ENVELOPES
%
%for i=1:nbins
%   if ~isempty(find(~isnan(hcpl(:,i))))
%      top(i)=max(find(~isnan(hcpl(:,i))));
%      bot(i)=min(find(~isnan(hcpl(:,i))));
%   else
%      top(i)=NaN;
%      bot(i)=NaN;
%   end
%end
%pcolor(xncpl(2:end),yncpl(2:end),hcpl);
%hold on 
%plot(ycpl(top),'linewidth',4)
%plot(ycpl(bot),'linewidth',4)
%
%X=[1:size(hncpl,2),size(hncpl,2):-1:1];
%Y=[bot,top(end:-1:1)];
%h=fill(X,Y,[.9,.9,.9]);
%set(h,'edgecolor','none','facealpha',.9)
%%set(h,'edgecolor','b','facealpha',.9)
%hold on ;
%plot(mean(MLD),'linewidth',2,'color','b')
%

c=[0,10,100,1e3,1e4,1e5,1e6];

figure(1);

pcolor(xncpl(2:end),yncpl(2:end),log(hncpl));

caxis(log([0 6e5]));
ylim([0.0002 0.005])
colorbar('FontSize',11,'YTick',log(c),'YTickLabel',c);
shading flat; grid
figure(2);


pcolor(xcpl(2:end),ycpl(2:end),log(hcpl));

caxis(log([0 6e5]));
ylim([0.0002 0.005])
colorbar('FontSize',11,'YTick',log(c),'YTickLabel',c);
shading flat; grid

mask=hcpl;mask(~isnan(mask))=1;
cdcplmean2d=squeeze(repmat(ycpl(2:end),[1,size(xcpl(2:end))])).*mask;
%cdcplmean=nanmean(cdcplmean2d);

mask=hncpl;mask(~isnan(mask))=1;
cdncplmean2d=squeeze(repmat(yncpl(2:end),[1,size(xncpl(2:end))])).*mask;
%cdncplmean=nanmean(cdncplmean2d);


cdncplmean=nansum(cdncplmean2d.*hncpl)./nansum(hncpl);
cdcplmean=nansum(cdcplmean2d.*hcpl)./nansum(hcpl);

stdncpl=std(cdncplmean2d,'omitnan');
stdcpl=std(cdcplmean2d,'omitnan');

stdncpltmp=zeros(size(stdncpl));
stdncpltmp(1:5:end)=stdncpl(1:5:end);
stdncpl=stdncpltmp;

stdcpltmp=zeros(size(stdcpl));
stdcpltmp(1:5:end)=stdcpl(1:5:end);
stdcpl=stdcpltmp;



figure(1)

hold on
%plot(xncpl(2:end),cdncplmean,'color','k','linewidth',2);
errorbar(xncpl(2:end),cdncplmean,stdncpl,'color','k','linewidth',2);
set(gca,'fontsize',13,'fontweight','bold')
xlabel('wind speed (m.s^{-1})')
ylabel('Cd')
title('uncoupled')
figure(2)

hold on 
%plot(xcpl(2:end),cdcplmean,'color','k','linewidth',2);
errorbar(xcpl(2:end),cdcplmean,stdcpl,'color','k','linewidth',2);
set(gca,'fontsize',13,'fontweight','bold')
xlabel('wind speed (m.s^{-1})')
ylabel('Cd')
title('coupled')
figure;

%plot(xncpl(2:end),cdncplmean,'linewidth',2);
errorbar(xncpl(2:end),cdncplmean,stdncpl,'linewidth',2);
set(gca,'fontsize',13,'fontweight','bold')
xlabel('wind speed (m.s^{-1})')
ylabel('N.m^{-2}')

hold on;
%plot(xcpl(2:end),cdcplmean,'linewidth',2);
errorbar(xcpl(2:end),cdcplmean,stdcpl,'linewidth',2);
legend('No\_CPL','WS\_CPL','location','southeast')

grid




