%========================================================================
%
% AIM: Climate change impacts on reliability from both supply-demand
%     
%      !! Considering variability and size
%   
%========================================================================
clc;
load Base_Energy.mat;

% 1) Pre-stored result

RlbH_frac100_glb_demhist585 = []; RlbH_frac0_glb_demhist585   = []; 
RlbL_frac100_glb_demhist585 = []; RlbL_frac0_glb_demhist585   = []; 


% 2) Loop each model

for mdnum  =  1 : size(CMIP6_name,2) 
    
     tic
     disp(['Now Load ',CMIP6_name{mdnum}])

   % 2.1) Load lat,lon,areaa,landfraction
     lat     = CMIP6_lata{mdnum,1};
     lon     = CMIP6_lona{mdnum,1};
     LandID  = CMIP6_ID_Landa{mdnum,1};
     OceanID = CMIP6_ID_Oceana{mdnum,1};
     area    = CMIP6_areaa{mdnum,1};
     Yr_Hist = CMIP6_histtime{mdnum,2}(:,1); Yr_Hist= Yr_Hist(Yr_Hist>=1950&Yr_Hist<=2014);  
     Yr_Fut  = CMIP6_ssptime{mdnum,2}(:,1); 
   
     st_Hist = min(find(Yr_Hist==1985)); ed_Hist = max(find(Yr_Hist==2014));len_Hist = ed_Hist-st_Hist+1;
     st_Fut  = min(find(Yr_Fut==2071)); ed_Fut = max(find(Yr_Fut==2100)); len_Fut = ed_Fut-st_Fut+1;
  
   % 2.2) Load Solar power (W/m2) and Capacity Factor (CF)        
     SPfiles = dir(['/net/exo/landclim/liula/wind_solar/Global/Output/Solar_Power2/historical/',CMIP6_name{mdnum},'/*.nc']);
         % ncdisp([SPfiles(1).folder,'/',SPfiles(1).name])
     SP_Hist = ncread([SPfiles(1).folder,'/',SPfiles(1).name],'PV',[1 1 st_Hist],[inf inf len_Hist],[1 1 1]);
     CF_PV_Hist = reshape(double(SP_Hist),[size(SP_Hist,1)*size(SP_Hist,2) size(SP_Hist,3)])./PV_max;
   
     SPfiles = dir(['/net/exo/landclim/liula/wind_solar/Global/Output/Solar_Power2/ssp585/',CMIP6_name{mdnum},'/*.nc']);
     SP_Fut  = ncread([SPfiles(1).folder,'/',SPfiles(1).name],'PV',[1 1 st_Fut],[inf inf len_Fut],[1 1 1]);
     CF_PV_Fut = reshape(double(SP_Fut),[size(SP_Fut,1)*size(SP_Fut,2) size(SP_Fut,3)])./PV_max;

   % 2.3) Load Wind power (W/m2) and Capacity Factor (CF)  
     WPfiles = dir(['/net/exo/landclim/liula/wind_solar/Global/Output/Wind_Power/historical/',CMIP6_name{mdnum},'/*.nc']);
         % ncdisp([WPfiles(1).folder,'/',WPfiles(1).name])
     WP_150  = ncread([WPfiles(1).folder,'/',WPfiles(1).name],'pv_wind150',[1 1 st_Hist],[inf inf len_Hist],[1 1 1]);
     WP_150  = reshape(WP_150,[size(WP_150,1)*size(WP_150,2) size(WP_150,3)]);
     CF_Wb_Hist = WP_150./WP_max;
     
     WPfiles = dir(['/net/exo/landclim/liula/wind_solar/Global/Output/Wind_Power/ssp585/',CMIP6_name{mdnum},'/*.nc']);
     if mdnum == 13 % In EC-Earth3, WP120m or WP150m are stored independently due to size constriant
         WPfiles  = WPfiles(2);
     end
     WP_150  = ncread([WPfiles(1).folder,'/',WPfiles(1).name],'pv_wind150',[1 1 st_Fut],[inf inf len_Fut],[1 1 1]);
     WP_150  = reshape(WP_150,[size(WP_150,1)*size(WP_150,2) size(WP_150,3)]);
     CF_Wb_Fut = WP_150./WP_max; 

   % 2.4) Load Air Temperature
      filelist = dir(['/net/exo/landclim/liula/wind_solar/Global/Output/tas/historical/',CMIP6_name{mdnum},'/*.nc']);
        % ncdisp([filelist(1).folder,'/',filelist(1).name])
      A   = ncread([filelist(1).folder,'/',filelist(1).name],'tas',[1 1 st_Hist],[inf inf len_Hist],[1 1 1]);
      Ta_Hist = reshape(A,size(A,1)*size(A,2),[]) - 273.15;   
      
      filelist = dir(['/net/exo/landclim/liula/wind_solar/Global/Output/tas/ssp585/',CMIP6_name{mdnum},'/*.nc']);
      A   = ncread([filelist(1).folder,'/',filelist(1).name],'tas',[1 1 st_Fut],[inf inf len_Fut],[1 1 1]);
      Ta_Fut  = reshape(A,size(A,1)*size(A,2),[]) - 273.15;   
  
   % 2.5) Pre-store
    
        % Pre-store reliabilty/storage globally [Final output] 
        RlbH_frac100 = nan(length(area(:)),2);
        RlbH_frac0   = nan(length(area(:)),2);
        
        RlbL_frac100 = nan(length(area(:)),2);
        RlbL_frac0   = nan(length(area(:)),2);
        
  % 2.6) Loop every grid cell 
    
   % Make temporal length of hist and fut the same to ensure hist-demand == fut-demand 
    [M,I] = min([size(CF_Wb_Hist,2),size(CF_Wb_Fut,2)]);
    CF_Wb_Hist =  CF_Wb_Hist(:,1:M);
    CF_PV_Hist =  CF_PV_Hist(:,1:M);
    CF_Wb_Fut  =  CF_Wb_Hist; % CF_Wb_Fut(:,1:M);
    CF_PV_Fut  =  CF_PV_Hist; % CF_PV_Fut(:,1:M);
    Ta_Hist    =  Ta_Hist(:,1:M);
    Ta_Fut     =  Ta_Fut(:,1:M); 
     
  for k =  LandID'  % Location
   
     %  Historical Supply and Demand profile 
        CF_Wb_hist=  CF_Wb_Hist(k,:)';
        CF_PV_hist=  CF_PV_Hist(k,:)';     
        % Generate historical demand profile according to Ta inter-day variations, assuming historical time-mean demand = 1
        Ta_hist   =  Ta_Hist(k,:)';
        [Demand_High_hist,Demand_Low_hist,S_High,S_Low] = Demandfix_Ta_Liu(1,Ta_hist);
    
     %  Future supply profile
        CF_Wb_fut =  CF_Wb_Fut(k,:)';
        CF_PV_fut =  CF_PV_Fut(k,:)';
     
     %  Based on Historical system, generate future Demand profile
        % Ajust demand profile due to climate change 
        Ta_fut    =  Ta_Fut(k,:)';
        [Demand_High_fut,~] = Demandmin_Ta_Liu(Ta_fut,S_High);
        [~, Demand_Low_fut] = Demandmin_Ta_Liu(Ta_fut,S_Low);

     % !! Historical Reliabilty in prescribed storage size   
      
         % Solar
        [RlbH_hist100,NumH_Wb100,NumH_PV100,BatH_size100,rlbH_hist100day] = Power_Storage_Liu(CF_Wb_hist,CF_PV_hist,Demand_High_hist,GR_Hist,1,stg_day_Hist);
        [RlbL_hist100,NumL_Wb100,NumL_PV100,BatL_size100,rlbL_hist100day] = Power_Storage_Liu(CF_Wb_hist,CF_PV_hist,Demand_Low_hist,GR_Hist,1,stg_day_Hist);
        
         % Wind
        [RlbH_hist0,NumH_Wb0,NumH_PV0,BatH_size0,rlbH_hist0day] = Power_Storage_Liu(CF_Wb_hist,CF_PV_hist,Demand_High_hist,GR_Hist,0,stg_day_Hist);
        [RlbL_hist0,NumL_Wb0,NumL_PV0,BatL_size0,rlbL_hist0day] = Power_Storage_Liu(CF_Wb_hist,CF_PV_hist,Demand_Low_hist,GR_Hist,0,stg_day_Hist);
          
            % Average 
        RlbH_frac100(k,1)  = RlbH_hist100;
        RlbH_frac0(k,1)    = RlbH_hist0;
        
        RlbL_frac100(k,1)  = RlbL_hist100;
        RlbL_frac0(k,1)    = RlbL_hist0;
        
     % !! Future Reliabilty in prescribed storage size   
      
        % Solar 
        [RlbH_fut100,~] = Power_Storage_fixLiu(CF_Wb_fut,NumH_Wb100,CF_PV_fut,NumH_PV100,Demand_High_fut,stg_day_Hist);
        [RlbL_fut100,~] = Power_Storage_fixLiu(CF_Wb_fut,NumL_Wb100,CF_PV_fut,NumL_PV100,Demand_Low_fut,stg_day_Hist);
       
        % Wind
        [RlbH_fut0,~] = Power_Storage_fixLiu(CF_Wb_fut,NumH_Wb0,CF_PV_fut,NumH_PV0,Demand_High_fut,stg_day_Hist);
        [RlbL_fut0,~] = Power_Storage_fixLiu(CF_Wb_fut,NumL_Wb0,CF_PV_fut,NumL_PV0,Demand_Low_fut,stg_day_Hist);
        
            % Average 
        RlbH_frac100(k,2) = RlbH_fut100;
        RlbH_frac0(k,2)   = RlbH_fut0;
        
        RlbL_frac100(k,2) = RlbL_fut100;
        RlbL_frac0(k,2)   = RlbL_fut0;
                
  end
   
    % Store reliability in hist/future at Prescribed storage size
    RlbH_frac100_glb_demhist585{mdnum,1}  = RlbH_frac100; 
    RlbH_frac0_glb_demhist585{mdnum,1}    = RlbH_frac0;
    
    RlbL_frac100_glb_demhist585{mdnum,1}  = RlbL_frac100; 
    RlbL_frac0_glb_demhist585{mdnum,1}    = RlbL_frac0;
    
    toc
    
end


% %  save Result_Energy.mat RlbH_frac100_glb_demhist585  RlbH_frac0_glb_demhist585  -append;
% %  save Result_Energy.mat RlbL_frac100_glb_demhist585  RlbL_frac0_glb_demhist585  -append;


% Regrid result into 1 degree   

   RlbH_frac100_glb1deg_demhist585 = nan(size(CMIP6_name,2),180*360,2); 
   RlbH_frac0_glb1deg_demhist585   = nan(size(CMIP6_name,2),180*360,2); 
 
   RlbL_frac100_glb1deg_demhist585 = nan(size(CMIP6_name,2),180*360,2); 
   RlbL_frac0_glb1deg_demhist585   = nan(size(CMIP6_name,2),180*360,2);  

for mdnum  = 1 : size(CMIP6_name,2) 
    
   disp(['Now Load ',CMIP6_name{mdnum}])

   % Load lat,lon,areaa,landfraction
     lata  = double(CMIP6_lata{mdnum,1});
     lona  = double(CMIP6_lona{mdnum,1});
   
   % Load Reliability
     rlbH100  = RlbH_frac100_glb_demhist585{mdnum,1};
     rlbH0    = RlbH_frac0_glb_demhist585{mdnum,1};
     
     rlbL100  = RlbL_frac100_glb_demhist585{mdnum,1};
     rlbL0    = RlbL_frac0_glb_demhist585{mdnum,1};
     
   % Interpolation
      lona2 = lona; lona2(lona2<0) = lona2(lona2<0) + 360; 
      [X,Y] = meshgrid([-89.5:1:89.5],[0.5:1:359.5]);
      
      for m = 1 : 2
          
          sub   = reshape(rlbH100(:,m),size(lata,1),size(lata,2));
          sub   = [sub(fix(size(sub,1)/2)+1:end,:);sub(1:fix(size(sub,1)/2),:)];
          Vq    = rot90(griddata(lata,lona2,sub,X,Y,'linear'));
          RlbH_frac100_glb1deg_demhist585(mdnum,:,m) = Vq(:);  
          
          sub   = reshape(rlbH0(:,m),size(lata,1),size(lata,2));
          sub   = [sub(fix(size(sub,1)/2)+1:end,:);sub(1:fix(size(sub,1)/2),:)];
          Vq    = rot90(griddata(lata,lona2,sub,X,Y,'linear'));
          RlbH_frac0_glb1deg_demhist585(mdnum,:,m) = Vq(:);  
          
          sub   = reshape(rlbL100(:,m),size(lata,1),size(lata,2));
          sub   = [sub(fix(size(sub,1)/2)+1:end,:);sub(1:fix(size(sub,1)/2),:)];
          Vq    = rot90(griddata(lata,lona2,sub,X,Y,'linear'));
          RlbL_frac100_glb1deg_demhist585(mdnum,:,m) = Vq(:);  
          
          sub   = reshape(rlbL0(:,m),size(lata,1),size(lata,2));
          sub   = [sub(fix(size(sub,1)/2)+1:end,:);sub(1:fix(size(sub,1)/2),:)];
          Vq    = rot90(griddata(lata,lona2,sub,X,Y,'linear'));
          RlbL_frac0_glb1deg_demhist585(mdnum,:,m) = Vq(:);  
             
      end
      
end


% % save Result_Energy.mat RlbH_frac100_glb1deg_demhist585 RlbH_frac0_glb1deg_demhist585 -append;
% % save Result_Energy.mat RlbL_frac100_glb1deg_demhist585 RlbL_frac0_glb1deg_demhist585 -append;





