%% Copyright notice
% Cite following references when you use or modify this function.
%
% Odériz, I.,  Mori, N., Shimura, T., Webb, A., Silva, R., Mortlock, T. R., (2022). 
% Transitional Wave Climate Regions on Continental and Polar Coasts in a Warming World. Nature
% Climate Change

% Odériz, I., Mori, N., Shimura, T., Webb, A., Silva, R., Mortlock, T. R., Dataset for Transitional
% Wave Climate Regions on Continental and Polar Coasts in a Warming World (Version 1). Zenodo. DOI: 10.5281/zenodo.6482839 (2022).

% --------------------------------------------------------------------
%       Itxaso Odériz
%       itxaso.oderiz@gmail.com	
%%   --------------------------------------------------------------------
clear all;close all;clc

%% Change filename and paths
nfile=1
FILE={'ENSAMBLE_'};
PathHIST=strcat('*\01.WaveClimates\CMIP5_HIST/');
PathRCP=strcat('*\01.WaveClimates\CMIP5_RCP8.5/');
PathOut=strcat('*\01.WaveClimates\');

L1=-30;
L2=30;

%% load data
        lon=ncread([PathHIST 'ENSEMBLE-HIST_WC.nc'],'lon');
        lat=ncread([PathHIST 'ENSEMBLE-HIST_WC.nc'],'lat');
        Lat=repmat(flip(lat),[1 length(lon)]);
        Lon=repmat(lon',[  length(lat) 1]);
        monthHIST=ncread([PathHIST 'ENSEMBLE-HIST_WC.nc'],'time');monthHIST=cell2mat(monthHIST);monthHIST=str2num(monthHIST(:,6:7));
        Pw_HIST=ncread([PathHIST 'ENSEMBLE-HIST_WC.nc'],'dirm');
        idxHIST=ncread([PathHIST 'ENSEMBLE-HIST_WC.nc'],'WC');
        monthRCP=ncread([PathRCP 'ENSEMBLE-RCP8.5_WC.nc'],'time');monthRCP=cell2mat(monthRCP);monthRCP=str2num(monthRCP(:,6:7));
        Pw_RCP=ncread([PathRCP 'ENSEMBLE-RCP8.5_WC.nc'],'dirm');
        idxRCP=ncread([PathRCP 'ENSEMBLE-RCP8.5_WC.nc'],'WC');
        
        [nr1 nc1 nt]=size(idxRCP);


%% calculo la diferencia y lo guardo en un archivo

for WCT=1:3 % wave climates
    
    varHIST=Pw_HIST;
    varHIST(idxHIST~=WCT)=NaN;

    varRCP=Pw_RCP;
    varRCP(idxRCP~=WCT)=NaN;

    
for i=1:5 % seasons
         aa=[];   dataHIST=[];seasonHIST=[];  dataRCP=[];H=[];seasonRCP=[];Estd=[];N=[];t_PwT=[];
         if i==1
             sea='year';
                dataHIST=varHIST;
                dataRCP=varRCP;
                seasonHIST=find(monthHIST==1 | monthHIST==2 | monthHIST==3 | monthHIST==4 | monthHIST==5 | monthHIST==6 | monthHIST==7 | monthHIST==8 | monthHIST==9 | monthHIST==10 | monthHIST==11 | monthHIST==12);
                seasonRCP=find(monthRCP==1 | monthRCP==2 | monthRCP==3 | monthRCP==4 | monthRCP==5 | monthRCP==6 | monthRCP==7 | monthRCP==8 | monthRCP==9 | monthRCP==10 | monthRCP==11 | monthRCP==12);
        elseif i==2
                sea='DJF';
                seasonHIST=find(monthHIST==1 | monthHIST==2 | monthHIST==12);
                seasonRCP=find(monthRCP==1 | monthRCP==2 | monthRCP==12);
                dataHIST=varHIST(:,:,seasonHIST);
                dataRCP=varRCP(:,:,seasonRCP);

        elseif i==3
                sea='MAM';
                seasonHIST=find(monthHIST==3 | monthHIST==4 | monthHIST==5);
                seasonRCP=find(monthRCP==3 | monthRCP==4 | monthRCP==5);
                dataHIST=varHIST(:,:,seasonHIST);
                dataRCP=varRCP(:,:,seasonRCP);

        elseif i==4
                sea='JJA';
                seasonHIST=find(monthHIST==6 | monthHIST==7 | monthHIST==8);
                seasonRCP=find(monthRCP==6 | monthRCP==7 | monthRCP==8);
                dataHIST=varHIST(:,:,seasonHIST);
                dataRCP=varRCP(:,:,seasonRCP);

        elseif i==5
                sea='SON';
                seasonHIST=find(monthHIST==9 | monthHIST==10 | monthHIST==11);
                seasonRCP=find(monthRCP==9 | monthRCP==10 | monthRCP==11);
                dataHIST=varHIST(:,:,seasonHIST);
                dataRCP=varRCP(:,:,seasonRCP);
         end
        
                %% significance test of the 40% of the sample
                N=[];Esdt=[];Neff=[]; T_PwT=[];  t_PwT=[]; aa=[];
                NN_RCP=nansum(dataRCP-dataRCP+1,3);
                NN_HIST=nansum(dataHIST-dataHIST+1,3);
                FRE_RCP=100*NN_RCP./length(seasonRCP);
                FRE_HIST=100*NN_HIST./length(seasonHIST);
               
                FRE_RCP(FRE_RCP<40)=NaN;FRE_RCP=FRE_RCP-FRE_RCP+1;
                FRE_HIST(FRE_HIST<40)=NaN;FRE_HIST=FRE_HIST-FRE_HIST+1;
                
                N=nansum(dataRCP-dataRCP+1,3);

                %% calculo la frecuencia 
                Estd=nanstd(dataRCP,[],3)./sqrt(N);
                
                %% calculo los grados de libertad son el número de la muestra más pequña menos 1.
                Neff=N-1;

                %% busco le valor crítico de la tabla T-STUDENT CON LOS GRADOS DE LIBERTAD Y P
                p=0.05;
                T_PwT=abs(tinv(1-p/2,Neff));
                
                %% utlizimos el estadistico t (variable aleatoria de una Normal estándar)
                t_PwT=abs((nanmean(dataRCP,3).*FRE_RCP-nanmean(dataHIST,3).*FRE_HIST)./Estd);
                
                %% compruebo si es estadisticamente significante
                H= abs(T_PwT<t_PwT);%PwTmean
                aa= nanmean(dataRCP,3).*FRE_RCP-nanmean(dataHIST,3).*FRE_HIST;
                R(:,2*i+1)=aa(:);
                R(:,2*i+2)=H(:);
                 
   %% Dibujo los resultados
    % factor to visualize dots
     nn=1;

   close all
   figura=figure(i)
   cmap= load ('ColorScale_BlueRed.txt'); 
   fileOut= strcat(FILE{nfile},'-WCT-',num2str(WCT),'-difRCP85-HIST_DIR-',sea,'.png');

   axesm('robinson','MapLatLimit',[-90 90],'MapLonLimit',[0 360])
%  axesm('mapprojection','ortho','origin',[-90 180 0])
%  axesm('mapprojection','ortho','origin',[ 90 180 0])

     geoshow(reshape(Lat,[nr1,nc1]),reshape(Lon,[nr1,nc1]),aa/nn,'DisplayType','surface','FaceColor','flat');
     cb=colorbar('horiz');
%     set(get(cb,'xlabel'),'string',{'Dir(°)',FILE{nfile}});
     set(cb, 'ylim', [L1/nn L2/nn]);
     colormap(cmap);
     caxis([L1/nn L2/nn]);

     load coast;
     geoshow(lat, long,'DisplayType', 'polygon', 'FaceColor', '[.6 .6 .6]','edgecolor','[.6 .6 .6]','LineWidth',1.5); 
     LatH=Lat;LonH=Lon;
     a=LatH(:,1)-LatH(:,1)+NaN;p1=1:4:length(a);a(p1)=1;
     b=LatH(1,:)-LatH(1,:)+NaN;p2=1:1:length(b);b(p2)=NaN;
     A=repmat (a,[1,length(b)]);
     LatH=LatH.*A;LonH=LonH.*A;
     LatH(H==0)=NaN;LonH(H==0)=NaN;
       
     geoshow(LatH(:), LonH(:), 'DisplayType', 'Point', 'Marker', '.', 'MarkerEdgeColor', 'k','MarkerSize', 0.001);
     setm(gca,'LabelFormat','signed');
     setm(gca,'MLabelParallel',0,'PLabelMeridian',10,'PLabelMeridian',-45);
     
     gridm on;
     title(strcat(sea,' difference RCP8.5-HIST'));
     saveas(figura,[PathOut fileOut]);

        
end


end



