% Organize wind data for WAP from ERA5
% August 19, 2022
% February 2, 2024
% Version 4: also calculates wind direction from u,v
%
datafolder = '/Users/jsturner/satellite_data/WAP_Wind/ERA5/';
dir1 = dir(datafolder);
lf=length(dir1);
filelist={};
for i=1:lf
    nf=dir1(i).name; tf=contains(nf,'.nc');  % is it an nc file
    if tf==1 % if this IS an nc file
        filelist=[filelist;nf];
    end
end; clearvars lf nf tf dir1 i

% Preallocate 3D array for wind speed and 1D array for date
date_all = [];
windspeed_daily = [];
winddir_daily = [];
u_day = [];
v_day = [];
for k = 1:length(filelist)
    disp('  ')
    disp(['File ' num2str(k) ' of ' num2str(length(filelist))]);
    % Load hourly wind data and time vector
    file1 = [datafolder char(filelist{k})];
    % ncdisp(file1)
    time_hours = ncread(file1,'time');
    u10 = ncread(file1,'u10');
    v10 = ncread(file1,'v10');
    % Convert hours to matlab date format
    date_mat = nan(length(time_hours),1);
    for i = 1:length(time_hours)
        hour = time_hours(i);
        date_mat(i) = addtodate(datenum(1900,01,01,0,0,0),hour,'hour');
    end; clearvars hour i

    % Calculate hourly windspeed, before averaging
    disp('Calculating hourly windspeed...');
    windspeed_hourly = nan(size(v10,1),size(v10,2),length(date_mat));
    for i = 1:length(date_mat)
        windspeed_hourly(:,:,i) = sqrt( (u10(:,:,i).^2) + (v10(:,:,i).^2) );
    end
    
    % Average into daily wind U10 and V10 for wind direction
    disp('Averaging hourly u,v to daily u,v...');
    time_days = (date_mat(1):1:date_mat(end))';
    u10_day = nan(size(u10,1),size(u10,2),length(time_days));
    v10_day = nan(size(u10,1),size(u10,2),length(time_days));
    for i = 1:length(time_days)
        day1 = time_days(i);
        day2 = time_days(i)+1;
        ind_thisday = find(date_mat >= day1 & date_mat < day2);
        u10_day (:,:,i) = nanmean(u10(:,:,ind_thisday),3);
        v10_day (:,:,i) = nanmean(v10(:,:,ind_thisday),3);
    end
    
    % Calculating daily wind direction from daily u,v
    disp('Calculating daily wind direction from daily u,v...');
    winddir_day = nan(size(v10,1),size(v10,2),length(time_days));
    Angles360 = @(a) rem(360+a, 360);
    for i = 1:length(time_days)
        a = atan2d(u10_day(:,:,i),v10_day(:,:,i));
        wdir = Angles360(a);
        winddir_day(:,:,i) = wdir;
    end
    
    % Example: file1 record starts 1-Jan-1997 and ends 31-Dec-2001
    disp('Averaging hourly windspeed to daily windspeed');
    windspeed_day = nan(size(u10,1),size(u10,2),length(time_days));
    for i = 1:length(time_days)
        day1 = time_days(i);
        day2 = time_days(i)+1;
        ind_thisday = find(date_mat >= day1 & date_mat < day2);
        windspeed_day(:,:,i) = nanmean(windspeed_hourly(:,:,ind_thisday),3);
    end

    date_all = [date_all; time_days];
    windspeed_daily = cat(3,windspeed_daily,windspeed_day);
    winddir_daily = cat(3,winddir_daily,winddir_day);
    u_day = cat(3,u_day,u10_day);
    v_day = cat(3,v_day,v10_day);
end; clearvars k windspeed_thisfile time_days day1 day2 time_hours windspeed_hourly

% Load latitude and longitude from last file
lat0 = ncread(file1,'latitude');
lon0 = ncread(file1,'longitude');
lon = repmat(lon0,1,size(windspeed_daily,2));
lat = repmat(lat0,1,size(windspeed_daily,1))';

clearvars a wdir ans directory file1 lat0 lon0 folder windspeed i u10 v10 date_mat
% metadata = ['File created by Jessie Turner in MATLAB from ERA5 hourly data using script "load_windspeed_daily_from_hourly_v3.m" on ' datestr(now)];
% 
% % Save
metadata =['This file was created in MATLAB from ERA5 hourly wind data (u10,v10) by Jessie Turner on ' datestr(now,'yyyy-mm-dd')];
save('daily_windspeed_era5_v2.mat','date_all','winddir_daily','windspeed_daily','u_day','v_day','lon','lat','metadata'); % about 290 MB
% 


