% % ERA5 output process every grid cell for U50
% % file output is PB_U50estimates.mat
clearvars
close all
delete(gcp('nocreate'))
parpool('local',24)
cd('/vol_b/data/wrf5/ERA5/hourly');
R50=NaN(1440,681,21);
U50=NaN(1440,681,8);
%this is used in Gumbel fit
for k=1:40;
    pU(k)=-1*log(-1*log((k-0.44)/(40+0.12)));
end
n=0;
for year=1979:2018;
    for month=1:12
        n=n+1;
        file{n,:}=char(strcat('ERA5_hourly_',num2str(year),num2str(month),'.nc'));
    end
end
for y1=21:701;
    y1
    clearvars -except y1 file pU U50 R50 
    WS=NaN(1,1,1440);
    A=NaN(480,1440);
    %read the full WS array to do Weibull fitting.
    %read also maximum values to compute Annual maxima
    parfor n=1:480;
        WS=[WS;....
            permute(sqrt(double(ncread(char(file{n,:}),'v100',[1 y1 1],[Inf 1 Inf])).^2+ ....
            +double(ncread(char(file{n,:}),'u100',[1 y1 1],[Inf 1 Inf])).^2),[3,2,1])];
        A(n,:)=max(sqrt(double(ncread(char(file{n,:}),'v100',[1 y1 1],[1440 1 Inf])).^2+ ....
            +double(ncread(char(file{n,:}),'u100',[1 y1 1],[1440 1 Inf])).^2),[],3);
    end
    Annmax=NaN(40,1440);
    for year=1:40;
        Annmax(year,:)=max(A((((year-1)*12+1):year*12),:));
    end
    clear A
    for x1=1:1440;
        %approach 1; 5 times the mean annual mean;
        U50(x1,y1-20,1)=5*nanmean(WS(2:350641,1,x1));
        %approach 2: Based on Weibull distribution paramaters. Weibull paramaters
        %are generated using maximum likelihood approaches.
        pd=fitdist(WS(2:350641,1,x1),'Weibull');
        %m=N*(1-r1)/(1+r1);
        %N is the total number of observations
        %r1 is the lag-1 correlation coefficient.
        m=length(WS(:,1,x1))*(1-(corr2(WS(2:350640,1,x1),WS(3:350641,1,x1))))/(1+(corr2(WS(2:350640,1,x1),WS(3:350641,1,x1))));
        alpha=1./(pd.ParameterValues(2)/pd.ParameterValues(1)*(log(m)^(1-1/pd.ParameterValues(2))));
        beta=pd.ParameterValues(1)*log(m)^(1/pd.ParameterValues(2));
        U50(x1,y1-20,2)=-1.*alpha'.*log(log((50)/(50-1)))+beta';
        %approach 3: GEV-Gumbel parameters from method of moments:
        U50(x1,y1-20,3)=-1*(sqrt(6)/pi*std(Annmax(:,x1)))*log(log((50)/(50-1)))+(mean(Annmax(:,x1))-0.577216*sqrt(6)/pi*std(Annmax(:,x1)));
        %approach 4: GEV-Gumbel paramaters from maximum likelihood fitting:
        pd=fitdist(Annmax(:,x1),'ExtremeValue');
        U50(x1,y1-20,4)=pd.ParameterValues(2)*log(-1*log(0.02))+pd.ParameterValues(1);
        pci1=paramci(pd,'Parameter','mu');
        pci2=paramci(pd,'Parameter','sigma');
        U50(x1,y1-20,5)=pci2(1)*log(-1*log(0.02))+pci1(1);
        U50(x1,y1-20,6)=pci2(2)*log(-1*log(0.02))+pci1(2);
        %approach 5: Gumbel plot
        mdl=fitlm(pU',sort(Annmax(:,x1)));
        U50(x1,y1-20,7)=predict(mdl,-1*log(-1*log(1-0.02)));
        U50(x1,y1-20,8)=mdl.Rsquared.Adjusted;
        %estimates for 20-year blocks;
        for k=1:21;
            pd=fitdist(Annmax((1+(k-1):20+(k-1)),x1),'ExtremeValue');
            R50(x1,y1-20,k)=pd.ParameterValues(2)*log(-1*log(0.02))+pd.ParameterValues(1);
        end
    end
end
clearvars -except U50 R50
lat=double(ncread('ERA5_hourly_19791.nc','latitude'));
lng=double(ncread('ERA5_hourly_19791.nc','longitude'));
a1=find(lng>180);
lng(a1)=lng(a1)-360;
for i=1:681;
    U(1:1440,i,9)=lng(1:1440);
end
for j=1:1440;
    U(j,1:681,10)=lat(21:701);
end

save('PB_U50estimates.mat','U50','R50','-v7.3');
