% Utqiagvik 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_arm = [S_name ':\YOPP\MODF\Data\ARM2018\'];
path_tskin = [S_name ':\YOPP\MODF\Data\NCDC2018\'];
path_bulk = [S_name ':\YOPP\MODF\Data\'];
path_igra = [S_name ':\YOPP\MODF\Data\IGRA\'];
path_prsn = [S_name ':\YOPP\MODF\Data\SR\'];
path_soil = [S_name ':\YOPP\MODF\Data\GMD2018\'];
path_ozone = [S_name ':\YOPP\MODF\Data\Ozone\'];
period = [2018 2 1 0 0 0; 2018 3 31 0 0 0];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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]';

% ARM tower data
[time_met, tas, hurs, tdps, uas, vas] = read_twr(path_arm,period);

% ARM radiation
[time_rad, rsds, rlds, rsus, rlus, pr, ps] = read_rad(path_arm, period);

% ARM ground heat flux
[time_ghf, hfds] = read_ghf(path_arm, period);

% ARM eddy covariance fluxes
[time_ec, hfss_ec, hfls_ec] = read_ec_upd(path_arm, period);

% NCDC skin temperature
[time_tskin, tskin] = read_tskin(path_tskin, period);
% put on MODF timeline
ts = NaN(length(time05),1);
ts(ismember(datevec(time05),datevec(time_tskin),'rows')) = tskin;

% GMD ozone
[time_o, o3] = read_ozone(path_ozone, period);
% put on MODF timeline
oz = NaN(length(time01),1);
oz(ismember(datevec(time01),datevec(time_o),'rows')) = o3;

% PSD bulk turbulent fluxes
[time_bulk, z0m, z0h, tauv, tauu, hfss_bulk, hfls_bulk] = read_bulk(path_bulk, period);

% GMD soil data
[time_soil, snd_soil, tsl_soil] = read_soil(path_soil,period);
tsl = NaN(10,length(time01));
tsl(:,ismember(datevec(time01),datevec(time_soil),'rows')) = tsl_soil;
snd = NaN(length(time01),1);
snd(ismember(datevec(time01),datevec(time_soil),'rows')) = snd_soil;

% PSD snowflux data
[time_prsn, height, sr] = read_prsn(path_prsn, period);
% put on MODF timeline
for i = 1:length(time01)
    [a,b] = min(abs(time_prsn - time01(i)));
    if abs(time_prsn(b) - time01(i)) < 1/24/60
        prsn(:,i) = sr(:,b);
    else
        prsn(:,i) = NaN;
    end
end

% height vars
height_radar = height*1000;
height_tower = [2;10;20;40];

% ARM and IGRA sonde data
arm = read_sonde(path_arm,period);
igra = read_igra(path_igra,period);
% 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
x = rd_netcdf([S_name ':\YOPP\MODF\Data\ARM2018\nsatwrC1.b1.20180201.000000.cdf']);
y1 = attribute_list([S_name ':\YOPP\ECMWF\barrow_ifs-ecmwf_2018032412.nc']);

FillValue = y1.ps.varname__FillValue;
lat = x.lat;
lon = x.lon;
orog = x.alt;
clear mex

% convert time
% time01 = round((time01 - datenum([2018 1 1 0 0 0])) * 24 * 60 )/60;
% time05 = round((time05 - datenum([2018 1 1 0 0 0])) * 24 * 60 )/60;
% time10 = round((time10 - datenum([2018 1 1 0 0 0])) * 24 * 60 )/60;
% time30 = round((time30 - 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);
time05 = int32((time05 - datenum([1970 1 1 0 0 0])) * 24 * 60);
time10 = int32((time10 - datenum([1970 1 1 0 0 0])) * 24 * 60);
time30 = int32((time30 - 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);


% write MODFile
MODF.time.time01 = time01;
MODF.time.time05 = time05;
MODF.time.time10 = time10;
MODF.time.time30 = time30;

MODF.height.height_tower = height_tower;
MODF.height.height_radar = height_radar;
MODF.height.depth_soilprobe = depth;

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

MODF.timeSeries.z0m = z0m;
MODF.timeSeries.z0h = z0h;
MODF.timeSeries.ps = ps;
MODF.timeSeries.uas = uas(:,1);
MODF.timeSeries.vas = vas(:,1);
MODF.timeSeries.tas = tas(:,1);
MODF.timeSeries.tdps = tdps(:,1);
MODF.timeSeries.hurs = hurs(:,1);
MODF.timeSeries.oz = oz;
MODF.timeSeries.snd = snd;
MODF.timeSeries.ts = ts;
MODF.timeSeries.rsus = rsus;
MODF.timeSeries.rsds = rsds;
MODF.timeSeries.rlus = rlus;
MODF.timeSeries.rlds = rlds;
MODF.timeSeries.hfls_ec = hfls_ec;
MODF.timeSeries.hfss_ec = hfss_ec;
MODF.timeSeries.hfls_bulk = hfls_bulk;
MODF.timeSeries.hfss_bulk = hfss_bulk;
MODF.timeSeries.hfds = hfds;
MODF.timeSeries.tauu = tauu;
MODF.timeSeries.tauv = tauv;

MODF.timeSeriesProfile.ua = uas';
MODF.timeSeriesProfile.va = vas';
MODF.timeSeriesProfile.ta = tas';
MODF.timeSeriesProfile.tdp = tdps';
MODF.timeSeriesProfile.hur = hurs';
MODF.timeSeriesProfile.tsl = tsl;
if any(~isnan(prsn(:)))
    MODF.timeSeriesProfile.prsn = prsn;
end

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;
MODF.radiosonde.alt_sonde = alt_sonde;
MODF.radiosonde.lat_sonde = lat_sonde;
MODF.radiosonde.lon_sonde = lon_sonde;
MODF.radiosonde.zg_sonde = zg;

read_HK_table('Utqiagvik_var_att_sop1.4.xlsx','utqiagvik',period,MODF,FillValue)
