%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%% Subroutine used in ACE_spatial.m that calculates ACE, TCI and Density of TC (all or major hurricanes only)
%% over the temporal period defined by [F_year - L_year]
%%      the spatial period defined  by [lon_min,lon_max,lat_min,lat_max]
%% using the best tracks data : WS_EP,LAT_EP,LON_EP,HOUR_EP,DAY_EP,MONTH_EP,YEAR_EP
%%
%%  As an integrated metric of TC activity, we use the Accumulated Cyclone Energy (ACE) index, defined as the sum of the squares of the maximum wind speeds over all six-hour periods during which the maximum wind speed is over 35kt for every storms [Bell et al., 2000]. We also considered the TC Intensity (TCI, i.e. the simple maximum wind-speed average) since it relates more to oceanic conditions than ACE, which encompasses atmospheric effects on storm intensification as well [Lin and Chan, 2015]. 
%%
%%
%%
%%
%% Bell, G. D., M. S. Halpert, R. C. Schnell, R. W. Higgins, J. Lawrimore, V. E. Kousky, R. Tinker, W. Thiaw, M. Chelliah, and A. Artusa, Climate %% assessment for 1999. Bull. Amer. Meteor. Soc., 81, S1-S50, (2000).
%%
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
function [ACEd2,INTENSITYd2,DENSITYalld2,ACEm2,INTENSITYm2,DENSITYallm2,DENSITYmajorm2,ACEm_anom2,INTENSITYm_anom2,DENSITYallm_anom2,DENSITYmajorm_anom2,Seasonal_ACEm,Seasonal_INTENSITYm,Seasonal_DENSITYallm,Seasonal_DENSITYmajorm] = calc_ACE_box_w_daily(F_year,L_year,lon_min,lon_max,lat_min,lat_max,WS_EP,LAT_EP,LON_EP,HOUR_EP,DAY_EP,MONTH_EP,YEAR_EP);




		disp(' ==================================================')
		disp('!                                                  !')		
		disp('! Compute ACE & Intensity                          !')
		disp('!                                                  !')
		disp('!         For the Period:                          !')
		disp('!         ---------------                          !')
		disp(['!          ',num2str(F_year),' - ',num2str(L_year),'                             !'])
		disp('!                                                  !')	
		disp('!         For the Region:                          !')
		disp('!         ---------------                          !')
		disp(['!          ',num2str(lon_min),' - ',num2str(lon_max),'E                       !'])
		disp(['!            ',num2str(lat_min),' - ',num2str(lat_max),'N                         !'])
		disp('!                                                  !')
		disp(' ==================================================')
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp(' ==> Extracting spatial region')
disp('   ')
  WS_EP(find(LAT_EP < lat_min))=0;% 
  WS_EP(find(LAT_EP > lat_max))=0;% Keep only storms occurring 
  WS_EP(find(LON_EP < lon_min))=0;% in the defined region
  WS_EP(find(LON_EP > lon_max))=0;% 
%
disp(' ==> Computing ACE, Intensity & Density')
disp('   ')
%
sq_WS=(WS_EP).^2;
for i=1:1:length(HOUR_EP(1,:))
		% time interval is 6-hour
		ACE(i)=(sum(sq_WS(find(mod(HOUR_EP(:,i),6) == 0),i)))/10000; % in [10^-4 Knots^2]
		INTENSITY(i)=(max(WS_EP(:,i)));                              % in [Knots]
end
%
TC=[];
TC_intensity=[];
TC_density=[];
%
for i=1:1:length(MONTH_EP(1,:));   
       TC=[TC;YEAR_EP(i) MONTH_EP(1,i) DAY_EP(1,i) ACE(i)];
       TC_intensity=[TC_intensity;YEAR_EP(i) MONTH_EP(1,i) DAY_EP(1,i) INTENSITY(i)];
       TC_density=[TC_density;YEAR_EP(i) MONTH_EP(1,i) DAY_EP(1,i) max(WS_EP(:,i)) 1];
end
%
% Create Monthly time series
% ==========================
%
for i=F_year:1:L_year;
   for m=1:1:12;
    ACEm(i-(F_year - 1),m)	  =	sum(TC(find(TC(:,1) == i & TC(:,2) == m),4));  % Cumulative sum every 6 hours of square wind speed for each month of each year of the defined period
    INTENSITYm(i-(F_year - 1),m)    =	nanmean(TC_intensity(find(TC_intensity(:,1) == i & TC_intensity(:,2) == m),4)); % Average of all maximum wind speed for each month of each year of the defined period
    DENSITYallm(i-(F_year - 1),m)   =	sum(TC_density(find(TC_density(:,1) == i & TC_density(:,2) == m & TC_density(:,4) > 32),5)); % Count of all storms above 32 knots (Tropical Depressions) for each month of each year of the defined period
    DENSITYmajorm(i-(F_year - 1),m) =	sum(TC_density(find(TC_density(:,1) == i & TC_density(:,2) == m & TC_density(:,4) > 94),5)); % Count of major hurricanes, i.e; above 94 knots (Category 3) for each month of each year of the defined period
   	for d=1:1:31; 
    ACEd(i-(F_year - 1),m,d)		=	sum(TC(find(TC(:,1) == i & TC(:,2) == m & TC(:,3) == d),4));
    INTENSITYd(i-(F_year - 1),m,d)	=	nanmean(TC_intensity(find(TC_intensity(:,1) == i & TC_intensity(:,2) == m & TC_intensity(:,3) == d),4));
    DENSITYalld(i-(F_year - 1),m,d)   =	sum(TC_density(find(TC_density(:,1) == i & TC_density(:,2) == m  & TC_intensity(:,3) == d & TC_density(:,4) > 32),5)); % Count of all storms above 32 knots (Tropical Depressions) for each month of each year of the defined period
	
   	end
	
    end
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% correct: Months 30/31 days and leap year for daily ACE/Intensity
disp(' ==> Correcting time')
disp('   ')
for i=F_year:1:L_year;
			
			if isleap(i)==0
			ACEd(i-(F_year - 1),2,29)=NaN;
			ACEd(i-(F_year - 1),2,30)=NaN;
			ACEd(i-(F_year - 1),2,31)=NaN;
			INTENSITYd(i-(F_year - 1),2,29)=NaN;
			INTENSITYd(i-(F_year - 1),2,30)=NaN;
			INTENSITYd(i-(F_year - 1),2,31)=NaN;
			DENSITYalld(i-(F_year - 1),2,29)=NaN;
			DENSITYalld(i-(F_year - 1),2,30)=NaN;
			DENSITYalld(i-(F_year - 1),2,31)=NaN;						
			elseif isleap(i)==1
			ACEd(i-(F_year - 1),2,30)=NaN;
			ACEd(i-(F_year - 1),2,31)=NaN;
			INTENSITYd(i-(F_year - 1),2,30)=NaN;
			INTENSITYd(i-(F_year - 1),2,31)=NaN;
			DENSITYalld(i-(F_year - 1),2,30)=NaN;
			DENSITYalld(i-(F_year - 1),2,31)=NaN;				
			end

end

ACEd(:,4,31)=NaN;
INTENSITYd(:,4,31)=NaN;
DENSITYalld(:,4,31)=NaN;

ACEd(:,6,31)=NaN;
INTENSITYd(:,6,31)=NaN;
DENSITYalld(:,6,31)=NaN;

ACEd(:,9,31)=NaN;
INTENSITYd(:,9,31)=NaN;
DENSITYalld(:,9,31)=NaN;

ACEd(:,11,31)=NaN;
INTENSITYd(:,11,31)=NaN;
DENSITYalld(:,11,31)=NaN;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate Anomalies
% ===================
%
if L_year - F_year >= 1
disp(' ==> Computing anomalies')
disp('   ')
		for i=F_year:1:L_year;
		   for m=1:1:12;
		   ACEm_anom(i-(F_year - 1),m) = ACEm(i-(F_year - 1),m) - nanmean(ACEm(:,m),1);
		   INTENSITYm_anom(i-(F_year - 1),m) = INTENSITYm(i-(F_year - 1),m) - nanmean(INTENSITYm(:,m),1);   
		   DENSITYallm_anom(i-(F_year - 1),m) = DENSITYallm(i-(F_year - 1),m) - nanmean(DENSITYallm(:,m),1);   
		   DENSITYmajorm_anom(i-(F_year - 1),m) = DENSITYmajorm(i-(F_year - 1),m) - nanmean(DENSITYmajorm(:,m),1);   
		   end
		end

end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute Seasonal cycle
% ======================
%
if L_year - F_year >= 1
disp(' ==> Computing seasonal cycle')
disp('   ')

		for m=1:1:12;
				   Seasonal_ACEm(m) 		= nanmean(ACEm(:,m),1);
				   Seasonal_INTENSITYm(m) 	= nanmean(INTENSITYm(:,m),1); 
				   Seasonal_DENSITYallm(m) 	= nanmean(DENSITYallm(:,m),1); 
				   Seasonal_DENSITYmajorm(m)= nanmean(DENSITYmajorm(:,m),1); 
				     
		end
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Make Timeseries vector for ACE/Intensity/density values and anomalies
% =====================================================================
%
cptm=0;
cpt=0;
for i=F_year:1:L_year;
   for m=1:1:12;
   cptm=cptm+1;
   ACEm2(cptm)=ACEm(i-(F_year - 1),m);
   INTENSITYm2(cptm)=INTENSITYm(i-(F_year - 1),m);
   DENSITYallm2(cptm)=DENSITYallm(i-(F_year - 1),m);
   DENSITYmajorm2(cptm)=DENSITYmajorm(i-(F_year - 1),m);
   
  	for d=1:1:31;
	cpt=cpt+1;
	ACEd2(cpt)=ACEd(i-(F_year - 1),m,d);
	INTENSITYd2(cpt)=INTENSITYd(i-(F_year - 1),m,d);		
	DENSITYalld2(cpt)=DENSITYalld(i-(F_year - 1),m,d);		
  	end
	
   end
end

ACEd2(isnan(ACEd2))=[];
INTENSITYd2(isnan(INTENSITYd2))=[];
DENSITYalld2(isnan(DENSITYalld2))=[];


if L_year - F_year >= 1

cptm=0;
cpt=0;
for i=F_year:1:L_year;
   for m=1:1:12;
   cptm=cptm+1;
   ACEm_anom2(cptm)=ACEm_anom(i-(F_year - 1),m);
   INTENSITYm_anom2(cptm)=INTENSITYm_anom(i-(F_year - 1),m);   
   DENSITYallm_anom2(cptm)=DENSITYallm_anom(i-(F_year - 1),m);  
   DENSITYmajorm_anom2(cptm)=DENSITYmajorm_anom(i-(F_year - 1),m); 
   end   
end

end