% Eureka MODF creator for SOP1
% change the S_name to your arctic server letter
% 
S_name = 'A'; %'X'


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% initialize your directories and the correct SOP
path = [S_name ':\YOPP\MODF\'];
path_tower = [S_name ':\eureka\surface_properties\fluxtower\towermet\products\2018\'];
path_rad1 = [S_name ':\eureka\radiometric\broadband_radiation\products\quality_controlled\2018\'];
path_rad2 = [S_name ':\eureka\surface_properties\fluxtower\towermet\products\2018\'];
path_igra = [S_name ':\YOPP\MODF\Data\IGRA\'];
period1 = [2018 2 1 0 0 0; 2018 3 31 0 0 0];
period2 = [2018 7 1 0 0 0; 2018 9 30 0 0 0];

period = period1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% time variables
time01 = [datenum(period(1,:)):1/24/60:datenum(period(2,:))+1 - 1/24/60]';
time05 = [datenum(period(1,:)):5/24/60:datenum(period(2,:))+1 - 5/24/60]';
time10 = [datenum(period(1,:)):10/24/60:datenum(period(2,:))+1 - 10/24/60]';
time30 = [datenum(period(1,:)):30/24/60:datenum(period(2,:))+1 - 30/24/60]';
time60 = [datenum(period(1,:)):1/24:datenum(period(2,:))+1 - 1/24]';
% Fluxtower data
[met, soil] = eur_twr_read(path_tower,period);

time_met = met(:,1);
time_soil = soil(:,1);

met(:,1) = [];
soil(:,1) = [];

% expand met and soil to time01
Met = NaN(length(time01),23);
Met(ismember(datevec(time01),datevec(time_met),'rows'),:) = met;

Soil = NaN(length(time01),16);
Soil(ismember(datevec(time01),datevec(time_soil),'rows'),:) = soil;

% Met data
ps = Met(:,1);
tas = Met(:,[6,4,2]);
hurs = [Met(:,7) Met(:,5) Met(:,3)]; 

uas = [-Met(:,20).*sind(Met(:,19)) -Met(:,18).*sind(Met(:,17))];
vas = [-Met(:,20).*cosd(Met(:,19)) -Met(:,18).*cosd(Met(:,17))];

% Soil data
hfds1 = Soil(:,1);
hfds2 = Soil(:,2);
hfds3 = Soil(:,3);
snd = Soil(:,4)./1000;
ts = Soil(:,5);
tsl = Soil(:,7:16);

% apply QC to Met and Soil
if datenum(period(1,:)) == datenum(period1(1,:))
    % period1
    psqc = ones(length(ps),1);
    ps = ps.*psqc;

    tasqc = ones(size(tas,1),3);
    tas = tas.*tasqc + 273.15;
    
    tslqc = ones(size(tsl,1),10);
    tsl = tsl.*tslqc + 273.15;
    
    tsqc = ones(size(ts,1),1);
    ts = ts.*tsqc + 273.15;

    hursqc = ones(size(hurs,1),3);
    hurs = hurs.*hursqc;
    
    sndqc = ones(length(snd),1);
    snd = snd.*sndqc;
else
    % period2
    psqc = ones(length(ps),1);
    ps = ps.*psqc;
    
    tasqc = ones(length(tas),3);
    tas = tas.*tasqc + 273.15;
    
    hursqc = ones(length(hurs),1);
    hurs = hurs.*hursqc;
    
    tslqc = ones(size(tsl,1),10);
    tsl = tsl.*tslqc + 273.15;
    
    tsqc = ones(size(ts,1),1);
    ts = ts.*tsqc + 273.15;
    
    uas(117700:end) = NaN;
    vas(117700:end) = NaN;
    
    snd(snd < 0.3 | snd > 0.5) = NaN;

end

% Radiation data
[rad1,~,rad2,~] = eur_rad_read(path_rad1, path_rad2, period);

time_rad1 = rad1(:,1);
rad1(:,1) = [];
time_rad2 = rad2(:,1);
rad2(:,1) = [];

Rad = NaN(length(time01),4);
Rad(ismember(datevec(time01),datevec(time_rad1),'rows'),1:4) = rad1(:,1:4);
Rad(ismember(datevec(time01),datevec(time_rad2),'rows'),5:6) = rad2(:,1:2);

rsds = Rad(:,1);
rlds = Rad(:,4);
rlus = Rad(:,5);
rsus = Rad(:,6);

rlus = NaN(length(rlus),1);

if datenum(period(1,:)) == datenum(period1(1,:))
    rsus_qc = ones(length(rsus),1);
    rsus_qc([53281:54720,57601:60480]) = NaN;
    rsus = rsus.*rsus_qc;
else
    rsus_qc = ones(length(rsus),1);
    rsus_qc([86401:89280,117698,117804,117912,117950,117966,117999,118022,118038,118083,118109,118162,118214,118238,118245,118319,118353,120800:120815]) = NaN;
    rsus = rsus.*rsus_qc;
end

% radiosondes
igra = read_igra_eureka(path_igra,period);
arm = [];
% put ARM and IGRA on the same timeline
[sonde_flag, time_rel_sonde, time_sonde, pfull, zg, ta, tdp, hur, ua, va, lat_sonde, lon_sonde, alt_sonde, release_time, count, sonde_index] = comb_sonde(arm,igra);


% constants
depth = [0.05 0.1 0.15 0.2 0.25 0.3 0.45 0.7 0.95 1.2];

% prepare to write
y1 = attribute_list([S_name ':\YOPP\ECMWF\barrow_ifs-ecmwf_2018032412.nc']);

FillValue = y1.ps.varname__FillValue;
lat = 79.9833;
lon = -85.9333;
orog = 10.04;

% convert time
% time01 = round((time01 - datenum([2018 1 1 0 0 0])) * 24 * 60 )/60;
% time60 = round((time60 - datenum([2018 1 1 0 0 0])) * 24 * 60 )/60;
% time_sonde = round((time_sonde - datenum([2018 1 1 0 0 0])) * 24 * 60 )/60;
% time_rel_sonde = round((time_rel_sonde - datenum([2018 1 1 0 0 0])) * 24 * 60 )/60;

time01 = int32((time01 - datenum([1970 1 1 0 0 0])) * 24 * 60);
%time60 = int32((time60 - datenum([1970 1 1 0 0 0])) * 24 * 60);
time_sonde = (time_sonde - datenum([1970 1 1 0 0 0])) * 24 * 60 * 60;
time_rel_sonde = int32((time_rel_sonde - datenum([1970 1 1 0 0 0])) * 24 * 60);

% height vars
height_tower_hur = [2;6;10];
height_tower_wsp = [6;11];

ps = ps*100;

MODF.time.time01 = time01;

MODF.height.height_tower_hur = height_tower_hur;
MODF.height.height_tower_wsp = height_tower_wsp;
MODF.height.depth_soilprobe = depth;

MODF.point.orog = orog;
MODF.point.lat = lat;
MODF.point.lon = lon;

MODF.timeSeries.ps = ps;
MODF.timeSeries.uas = uas(:,1);
MODF.timeSeries.vas = vas(:,1);
MODF.timeSeries.tas = tas(:,1);
MODF.timeSeries.hurs = hurs(:,1);
MODF.timeSeries.snd = snd;
MODF.timeSeries.ts = ts;
MODF.timeSeries.rsus = rsus;
MODF.timeSeries.rsds = rsds;
MODF.timeSeries.rlds = rlds;
MODF.timeSeries.hfds_1 = hfds1;
MODF.timeSeries.hfds_2 = hfds2;
MODF.timeSeries.hfds_3 = hfds3;

MODF.timeSeriesProfile.ta = tas';
MODF.timeSeriesProfile.ua = uas';
MODF.timeSeriesProfile.va = vas';
MODF.timeSeriesProfile.hur = hurs';
MODF.timeSeriesProfile.tsl = tsl';

MODF.trajectoryProfile.pa = pfull;
MODF.trajectoryProfile.ua = ua;
MODF.trajectoryProfile.va = va;
MODF.trajectoryProfile.ta = ta;
MODF.trajectoryProfile.tdp = tdp;
MODF.trajectoryProfile.hur = hur;

MODF.radiosonde.time_sonde = time_sonde;
MODF.radiosonde.time_release_sonde = time_rel_sonde;
MODF.radiosonde.release_time = release_time';
MODF.radiosonde.count_per_sonde = count;
MODF.radiosonde.index_sonde = sonde_index;
MODF.radiosonde.source_sonde = sonde_flag;

read_HK_table('Eureka_var_att_sop1.4.xlsx','eureka',period,MODF,FillValue)

