clear;clc;close all;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%setup_nctoolbox (need to use nctoolbox)
% The ERA5 TC data is included in the ty_info folder. The code for tropical cyclone identification is available from 
% https://gitlab.act.reading.ac.uk/track.
% The ERA5 climate reanalysis was generated and distributed by the European Centre for Medium-Range Weather Forecasts 
% (ECMWF; https://www.ecmwf.int/en/forecasts/datasets/ reanalysis-datasets/era5).
% For more information, please contact jianshi@hhu.edu.cn
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
era_folder='hs';
for year=1979:2022
tyinfo=load(['./ty_info/tyinfo_' num2str(year) '.txt']);
tyindex=load(['./ty_info/tyindex_' num2str(year) '.txt']);
tynum=tyindex(1,1);
lbeg=0;
lend=0;
year
for nn=1:tynum
    nn
    tyhr=tyindex(nn+1,2);
    
    tyid=tyindex(nn+1,1);
    lbeg=lend+1;
    lend=lend+tyhr;
    for i=lbeg:lend
        mon=floor(tyinfo(i,1)/10000)-floor(tyinfo(i,1)/1000000)*100;
        day=floor(tyinfo(i,1)/100)-floor(tyinfo(i,1)/10000)*100;
        itime=tyinfo(i,1)-floor(tyinfo(i,1)/100)*100;
        clon=tyinfo(i,2);
        clat=tyinfo(i,3);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        if mon<10
            filename= [ 'Hs' num2str(year) '0' num2str(mon)];
        else
            filename= [ 'Hs' num2str(year)  num2str(mon)];
        end
        if day<10
            filename =[filename '0'  num2str(day) '.grib'];
        else
            filename =[filename   num2str(day) '.grib' ];
        end
       
        grib=['./' era_folder '_t/',filename] 
        
        file_list=[grib];
        %
        nco=ncgeodataset(file_list);
        wave_list = nco.variables;
        %       
        if itime==0
            id_t=1;
        elseif itime==6
            id_t=2;
        elseif itime==12
            id_t=3;
        elseif itime==18
            id_t=4;
        end
        
        data_hs=nco.geovariable(wave_list(1));%
        data_t=nco.geovariable(wave_list(3));%
        data_lon=nco.geovariable(wave_list(6));%
        data_lat=nco.geovariable(wave_list(5));%
        hs=squeeze(data_hs.data(1,id_t,:,:));
        t=data_t.data(:);
        lat=data_lat.data(:);
        lon=data_lon.data(:);
 
        outputtywave(hs,clon,clat,lon,lat,tyid,year,mon,day,itime,era_folder)
        
    end
end
end

function []=outputtywave(hs,clon,clat,lon,lat,tyid,year,mon,day,hour,era_folder)
[m n]=size(hs);
disc=15.;
for i=1:m
    for j=1:n
        dis=((lon(j)-clon)^2+(lat(i)-clat)^2)^0.5  ;
        if dis>=disc
            hs(i,j)=NaN;
        end
        
    end
end
if mon<10
    filename= [ num2str(tyid) '_' num2str(year) '0' num2str(mon)];
else
    filename= [ num2str(tyid) '_' num2str(year)  num2str(mon)];
end
if day<10
    filename =[filename '0'  num2str(day)];
else
    filename =[filename   num2str(day)];
end
if hour<10
    filename=[filename  '0' num2str(hour)];
else
    filename=[filename  num2str(hour)];
end

str1=['./'  era_folder '_ERA5_results/' era_folder '_' filename];
save(str1,'hs') 
end

