function read_HK_table(filename,site,period,MODF,FillValue)

GlTb = readcell(filename,'Sheet','Global Attributes');
TsTb = readcell(filename,'Sheet','timeSeries Variable Attributes');
TsPrTb = readcell(filename,'Sheet','timeSeriesProfile Variable Attr');
TrPrTb = readcell(filename,'Sheet','trajectoryProfile Variable Attr');

n = size(TsTb,1);
n1 = size(TsTb,2);
n2 = size(TsPrTb,2);
n3 = size(TrPrTb,2);

modf_name1 = [site '_obs_timeSeries_' datestr(datenum(period(1,:)),'yyyymmdd') '_' datestr(datenum(period(2,:)),'yyyymmdd') '.nc'];
modf_name2 = [site '_obs_timeSeriesProfile_' datestr(datenum(period(1,:)),'yyyymmdd') '_' datestr(datenum(period(2,:)),'yyyymmdd') '.nc'];
modf_name3 = [site '_obs_trajectoryProfile_' datestr(datenum(period(1,:)),'yyyymmdd') '_' datestr(datenum(period(2,:)),'yyyymmdd') '.nc'];

% timeSeries file

% time dimensions
fn = fieldnames(MODF.time);
for i = 1:length(fn)
    t = MODF.time.(fn{i});
    dt = t(2)-t(1);
    nccreate(modf_name1,fn{i},'DIMENSIONS',{fn{i},length(MODF.time.(fn{i}))})
    ncwriteatt(modf_name1,fn{i},'long_name',['Valid time for observations with ' num2str(dt) '-minute cadence'])
    ncwriteatt(modf_name1,fn{i},'units','minutes since 1970-01-01 00:00:00')
    ncwriteatt(modf_name1,fn{i},'delta_t',datestr(double(dt)/60/24,'0000-00-00 HH:MM:SS'));
    ncwriteatt(modf_name1,fn{i},'standard_name','time')
    ncwriteatt(modf_name1,fn{i},'calendar','standard')
    ncwrite(modf_name1,fn{i}, MODF.time.(fn{i}))
end

% Global attributes
for i = 1:length(GlTb)
    if strcmp(GlTb{i,1},'featureType')
        ncwriteatt(modf_name1,'/',GlTb{i,1},'timeSeries');
    elseif strcmp(GlTb{i,1},'date_created')
        ncwriteatt(modf_name1,'/',GlTb{i,1},datestr(now,'yyyy-mm-dd HH:MM:SS'));
    elseif any(~ismissing(GlTb{i,2})) 
        if isa(GlTb{i,2},'datetime')
            ncwriteatt(modf_name1,'/',GlTb{i,1},datestr(GlTb{i,2},'yyyy-mm-dd HH:MM:SS'));
        else
            ncwriteatt(modf_name1,'/',GlTb{i,1},GlTb{i,2});
        end
    end
end

% single_level fixed variables
nccreate(modf_name1,'orog')
nccreate(modf_name1,'lat')
nccreate(modf_name1,'lon')

ncwriteatt(modf_name1,'orog','long_name','Surface altitude');
ncwriteatt(modf_name1,'orog','units','m');
ncwriteatt(modf_name1,'orog','standard_name','surface_altitude');
ncwriteatt(modf_name1,'orog','source','surface observation');
ncwriteatt(modf_name1,'orog','references',['https://psl.noaa.gov/arctic/observatories/' site '/']);
ncwriteatt(modf_name1,'orog','comment','psl.data@noaa.gov');
ncwrite(modf_name1,'orog',MODF.point.orog)

ncwriteatt(modf_name1,'lat','long_name','Latitude');
ncwriteatt(modf_name1,'lat','units','degrees_north');
ncwriteatt(modf_name1,'lat','standard_name','latitude')
ncwriteatt(modf_name1,'lat','source','surface observation');
ncwriteatt(modf_name1,'lat','references',['https://psl.noaa.gov/arctic/observatories/' site '/']);
ncwriteatt(modf_name1,'lat','comment','psl.data@noaa.gov');
ncwrite(modf_name1,'lat',MODF.point.lat)

ncwriteatt(modf_name1,'lon','long_name','Longitude');
ncwriteatt(modf_name1,'lon','units','degrees_east');
ncwriteatt(modf_name1,'lon','standard_name','longitude')
ncwriteatt(modf_name1,'lon','source','surface observation');
ncwriteatt(modf_name1,'lon','references',['https://psl.noaa.gov/arctic/observatories/' site '/']);
ncwriteatt(modf_name1,'lon','comment','psl.data@noaa.gov');
ncwrite(modf_name1,'lon',MODF.point.lon)

% surface variables
fn1 = fieldnames(MODF.timeSeries);

for i = 2:n1
    if any(strcmp(fn1,TsTb{1,i}))
        K = 0;
        for k =1:length(fn)
            if length(MODF.time.(fn{k})) == length(MODF.timeSeries.(TsTb{1,i}))
                K = k;
                break
            end
        end
        if K > 0
            nccreate(modf_name1,TsTb{1,i},'DIMENSIONS',{fn{K}})
            for j = 2:n
                if strcmp(TsTb{j,1},'actual_range')
                    ncwriteatt(modf_name1,TsTb{1,i},TsTb{j,1},[min(MODF.timeSeries.(TsTb{1,i})) max(MODF.timeSeries.(TsTb{1,i}))]);
                elseif strcmp(TsTb{j,1},'missing_value')
                    ncwriteatt(modf_name1,TsTb{1,i},TsTb{j,1},FillValue);
                elseif any(~ismissing(TsTb{j,i}))
                    ncwriteatt(modf_name1,TsTb{1,i},TsTb{j,1},TsTb{j,i});
                end

            end
            MODF.timeSeries.(TsTb{1,i}) = fillmissing(MODF.timeSeries.(TsTb{1,i}),'constant',FillValue);
            ncwrite(modf_name1,TsTb{1,i},MODF.timeSeries.(TsTb{1,i}));
        else
            disp(['Dimensions of ' TsTb{1,i} ' are missing.'])
        end
    else
        disp(['missing ' TsTb{1,i}])
    end
end


% timeSeriesProfile file

% time dimensions
% add a test for time vars present
for i = 1:1%length(fn)
    t = MODF.time.(fn{i});
    dt = t(2)-t(1);
    nccreate(modf_name2,fn{i},'DIMENSIONS',{fn{i},length(MODF.time.(fn{i}))})
    ncwriteatt(modf_name2,fn{i},'long_name',['Valid time for observations with ' num2str(dt) '-minute cadence'])
    ncwriteatt(modf_name2,fn{i},'units','minutes since 1970-01-01 00:00:00')
    ncwriteatt(modf_name2,fn{i},'delta_t',datestr(double(dt)/60/24,'0000-00-00 HH:MM:SS'));
    ncwriteatt(modf_name2,fn{i},'standard_name','time')
    ncwriteatt(modf_name2,fn{i},'calendar','standard')
    ncwrite(modf_name2,fn{i}, MODF.time.(fn{i}))
end

% height dimensions
% fix the attributes for heights
fnh = fieldnames(MODF.height);
h_str = '';
for i = 1:length(fnh)
    nccreate(modf_name2,fnh{i},'DIMENSIONS',{fnh{i},length(MODF.height.(fnh{i}))})
    if any(strfind(fnh{i},'depth_soilprobe'))
        ncwriteatt(modf_name2,fnh{i},'long_name','Soil observation depth')
        ncwriteatt(modf_name2,fnh{i},'units','m')
        ncwriteatt(modf_name2,fnh{i},'standard_name','depth')
    else
        if any(strfind(fnh{i},'ta'))
            h_str = 'temperature ';
        elseif any(strfind(fnh{i},'hur'))
            h_str = 'humidity ';
        elseif any(strfind(fnh{i},'wsp'))
            h_str = 'wind speed ';
        end  
        ncwriteatt(modf_name2,fnh{i},'long_name',['Tower ' h_str 'observation height'])
        ncwriteatt(modf_name2,fnh{i},'units','m')
        ncwriteatt(modf_name2,fnh{i},'standard_name','height')
    end
    ncwrite(modf_name2,fnh{i}, MODF.height.(fnh{i}))
end

% Global attributes
for i = 1:length(GlTb)
    if strcmp(GlTb{i,1},'featureType')
        ncwriteatt(modf_name2,'/',GlTb{i,1},'timeSeriesProfile');
    elseif strcmp(GlTb{i,1},'date_created')
        ncwriteatt(modf_name2,'/',GlTb{i,1},datestr(now,'yyyy-mm-dd HH:MM:SS'));
    elseif any(~ismissing(GlTb{i,2}))
        if isa(GlTb{i,2},'datetime')
            ncwriteatt(modf_name2,'/',GlTb{i,1},datestr(GlTb{i,2},'yyyy-mm-dd HH:MM:SS'));
        else
            ncwriteatt(modf_name2,'/',GlTb{i,1},GlTb{i,2});
        end
    end
end
% single_level fixed variables
nccreate(modf_name2,'orog')
nccreate(modf_name2,'lat')
nccreate(modf_name2,'lon')

ncwriteatt(modf_name2,'orog','long_name','Surface altitude');
ncwriteatt(modf_name2,'orog','units','m');
ncwriteatt(modf_name2,'orog','standard_name','surface_altitude');
ncwriteatt(modf_name2,'orog','source','surface observation');
ncwriteatt(modf_name2,'orog','references',['https://psl.noaa.gov/arctic/observatories/' site '/']);
ncwriteatt(modf_name2,'orog','comment','psl.data@noaa.gov');
ncwrite(modf_name2,'orog',MODF.point.orog)

ncwriteatt(modf_name2,'lat','long_name','Latitude');
ncwriteatt(modf_name2,'lat','units','degrees_north');
ncwriteatt(modf_name2,'lat','standard_name','latitude')
ncwriteatt(modf_name2,'lat','source','surface observation');
ncwriteatt(modf_name2,'lat','references',['https://psl.noaa.gov/arctic/observatories/' site '/']);
ncwriteatt(modf_name2,'lat','comment','psl.data@noaa.gov');
ncwrite(modf_name2,'lat',MODF.point.lat)

ncwriteatt(modf_name2,'lon','long_name','Longitude');
ncwriteatt(modf_name2,'lon','units','degrees_east');
ncwriteatt(modf_name2,'lon','standard_name','longitude')
ncwriteatt(modf_name2,'lon','source','surface observation');
ncwriteatt(modf_name2,'lon','references',['https://psl.noaa.gov/arctic/observatories/' site '/']);
ncwriteatt(modf_name2,'lon','comment','psl.data@noaa.gov');
ncwrite(modf_name2,'lon',MODF.point.lon)

% tower data
fn2 = fieldnames(MODF.timeSeriesProfile);

for i = 2:n2
    if any(strcmp(fn2,TsPrTb{1,i}))
        K1 = 0;
        K2 = 0;
        for k2 = 1:length(fn)
            if length(MODF.time.(fn{k2})) == size(MODF.timeSeriesProfile.(TsPrTb{1,i}),2)
                K2 = k2;
                break
            end
        end
        
        for k1 = 1:length(fnh)
            if length(MODF.height.(fnh{k1})) == size(MODF.timeSeriesProfile.(TsPrTb{1,i}),1)
                K1 = k1;
                break
            end
        end
        if K1 > 0 && K2 > 0
            nccreate(modf_name2,TsPrTb{1,i},'DIMENSIONS',{fnh{k1},fn{k2}})
            for j = 2:n
                if strcmp(TsPrTb{j,1},'actual_range')
                    ncwriteatt(modf_name2,TsPrTb{1,i},TsPrTb{j,1},[min(MODF.timeSeriesProfile.(TsPrTb{1,i}),[],'all') max(MODF.timeSeriesProfile.(TsPrTb{1,i}),[],'all')]);
                elseif strcmp(TsPrTb{j,1},'missing_value')
                    ncwriteatt(modf_name2,TsPrTb{1,i},TsPrTb{j,1},FillValue);
                elseif any(~ismissing(TsPrTb{j,i}))
                    ncwriteatt(modf_name2,TsPrTb{1,i},TsPrTb{j,1},TsPrTb{j,i});
                end
            end
            MODF.timeSeriesProfile.(TsPrTb{1,i}) = fillmissing(MODF.timeSeriesProfile.(TsPrTb{1,i}),'constant',FillValue);
            ncwrite(modf_name2,TsPrTb{1,i},MODF.timeSeriesProfile.(TsPrTb{1,i}));
        else
            disp(['Dimensions of ' TsPrTb{1,i} ' are missing.'])
        end
    else
        disp(['missing ' TsPrTb{1,i}])
    end
end

% trajectoryProfile file
% single_level fixed variables
nccreate(modf_name3,'orog')
nccreate(modf_name3,'lat')
nccreate(modf_name3,'lon')

% Global attributes
for i = 1:length(GlTb)
    if strcmp(GlTb{i,1},'featureType')
        ncwriteatt(modf_name3,'/',GlTb{i,1},'trajectoryProfile');
    elseif strcmp(GlTb{i,1},'date_created')
        ncwriteatt(modf_name3,'/',GlTb{i,1},datestr(now,'yyyy-mm-dd HH:MM:SS'));
    elseif any(~ismissing(GlTb{i,2}))
        if isa(GlTb{i,2},'datetime')
            ncwriteatt(modf_name3,'/',GlTb{i,1},datestr(GlTb{i,2},'yyyy-mm-dd HH:MM:SS'));
        else
            ncwriteatt(modf_name3,'/',GlTb{i,1},GlTb{i,2});
        end
    end
end

ncwriteatt(modf_name3,'orog','long_name','Surface altitude');
ncwriteatt(modf_name3,'orog','units','m');
ncwriteatt(modf_name3,'orog','standard_name','surface_altitude');
ncwriteatt(modf_name3,'orog','source','surface observation');
ncwriteatt(modf_name3,'orog','references',['https://psl.noaa.gov/arctic/observatories/' site '/']);
ncwriteatt(modf_name3,'orog','comment','psl.data@noaa.gov');
ncwrite(modf_name3,'orog',MODF.point.orog)

ncwriteatt(modf_name3,'lat','long_name','Latitude');
ncwriteatt(modf_name3,'lat','units','degrees_north');
ncwriteatt(modf_name3,'lat','standard_name','latitude')
ncwriteatt(modf_name3,'lat','source','surface observation');
ncwriteatt(modf_name3,'lat','references',['https://psl.noaa.gov/arctic/observatories/' site '/']);
ncwriteatt(modf_name3,'lat','comment','psl.data@noaa.gov');
ncwrite(modf_name3,'lat',MODF.point.lat)

ncwriteatt(modf_name3,'lon','long_name','Longitude');
ncwriteatt(modf_name3,'lon','units','degrees_east');
ncwriteatt(modf_name3,'lon','standard_name','longitude')
ncwriteatt(modf_name3,'lon','source','surface observation');
ncwriteatt(modf_name3,'lon','references',['https://psl.noaa.gov/arctic/observatories/' site '/']);
ncwriteatt(modf_name3,'lon','comment','psl.data@noaa.gov');
ncwrite(modf_name3,'lon',MODF.point.lon)

% automate this section

nccreate(modf_name3,'sample','DIMENSIONS',{'sample',length(MODF.radiosonde.time_sonde)})
nccreate(modf_name3,'source_length','DIMENSIONS',{'source_length',3})
nccreate(modf_name3,'date_length','DIMENSIONS',{'date_length',14})
nccreate(modf_name3,'time_sonde','DIMENSIONS',{'sample'})
nccreate(modf_name3,'time_release_sonde','DIMENSIONS',{'time_release_sonde',length(MODF.radiosonde.time_release_sonde)})
nccreate(modf_name3,'release_time','DIMENSIONS',{'time_release_sonde','date_length'},'Datatype','char')
nccreate(modf_name3,'count_per_sonde','DIMENSIONS',{'time_release_sonde'},'Datatype','int32')
nccreate(modf_name3,'index_sonde','DIMENSIONS',{'time_release_sonde'},'Datatype','int32')
nccreate(modf_name3,'source_sonde','DIMENSIONS',{'time_release_sonde','source_length'},'Datatype','char')

time_sonde = fillmissing(MODF.radiosonde.time_sonde,'constant',FillValue);
ncwrite(modf_name3,'time_sonde',time_sonde)
ncwriteatt(modf_name3,'time_sonde','long_name','Radiosonde valid time')
ncwriteatt(modf_name3,'time_sonde','units','seconds')
ncwriteatt(modf_name3,'time_sonde','standard_name','time')
ncwriteatt(modf_name3,'time_sonde','missing_value',FillValue)
ncwriteatt(modf_name3,'time_sonde','calendar','standard')
ncwriteatt(modf_name3,'time_sonde','comment','seconds since 1970-01-01 00:00:00')

ncwrite(modf_name3,'time_release_sonde',MODF.radiosonde.time_release_sonde);
ncwriteatt(modf_name3,'time_release_sonde','long_name','Radiosonde release time');
ncwriteatt(modf_name3,'time_release_sonde','units','minutes since 1970-01-01 00:00:00');
ncwriteatt(modf_name3,'time_release_sonde','standard_name','time')
ncwriteatt(modf_name3,'time_release_sonde','calendar','standard')

ncwrite(modf_name3,'release_time',MODF.radiosonde.release_time');
ncwriteatt(modf_name3,'release_time','long_name','Radiosonde release time in character form');
ncwriteatt(modf_name3,'release_time','cf_role','trajectory_id');

ncwriteatt(modf_name3,'sample','long_name','Index (instance dimension) of radiosonde observations in file')

ncwriteatt(modf_name3,'count_per_sonde','long_name','Number of samples for this radiosonde');
ncwriteatt(modf_name3,'count_per_sonde','sample_dimension','sample');
ncwrite(modf_name3,'count_per_sonde',MODF.radiosonde.count_per_sonde);

ncwriteatt(modf_name3,'index_sonde','long_name','Index of the first sample for this radiosonde');
ncwriteatt(modf_name3,'index_sonde','sample_dimension','sample');
ncwrite(modf_name3,'index_sonde',MODF.radiosonde.index_sonde);

ncwriteatt(modf_name3,'source_sonde','long_name','Source of radiosonde data in character form');
ncwrite(modf_name3,'source_sonde',MODF.radiosonde.source_sonde);

ncwriteatt(modf_name3,'date_length','long_name','First dimension (number of characters) of release_time')

if any(cellfun(@any,strfind(fieldnames(MODF.radiosonde),'alt_sonde')))
    nccreate(modf_name3,'lat_sonde','DIMENSIONS',{'sample'})
    nccreate(modf_name3,'lon_sonde','DIMENSIONS',{'sample'})
    nccreate(modf_name3,'alt_sonde','DIMENSIONS',{'sample'})
    
    ncwriteatt(modf_name3,'lat_sonde','long_name','Radiosonde latitude');
    ncwriteatt(modf_name3,'lat_sonde','units','degrees_north');
    ncwriteatt(modf_name3,'lat_sonde','standard_name','latitude');
    ncwriteatt(modf_name3,'lat_sonde','missing_value',FillValue);
    ncwriteatt(modf_name3,'lat_sonde','actual_range',[min(MODF.radiosonde.lat_sonde) max(MODF.radiosonde.lat_sonde)]);
    ncwriteatt(modf_name3,'lat_sonde','instrument','');
    ncwriteatt(modf_name3,'lat_sonde','source','');
    lat = fillmissing(MODF.radiosonde.lat_sonde,'constant',FillValue);
    ncwrite(modf_name3,'lat_sonde',lat);

    ncwriteatt(modf_name3,'lon_sonde','long_name','Radiosonde longitude');
    ncwriteatt(modf_name3,'lon_sonde','units','degrees_east');
    ncwriteatt(modf_name3,'lon_sonde','standard_name','longitude');
    ncwriteatt(modf_name3,'lon_sonde','missing_value',FillValue)
    ncwriteatt(modf_name3,'lon_sonde','actual_range',[min(MODF.radiosonde.lon_sonde) max(MODF.radiosonde.lon_sonde)]);
    ncwriteatt(modf_name3,'lon_sonde','instrument','');
    ncwriteatt(modf_name3,'lon_sonde','source','');
    lon = fillmissing(MODF.radiosonde.lon_sonde,'constant',FillValue);
    ncwrite(modf_name3,'lon_sonde',lon);

    ncwriteatt(modf_name3,'alt_sonde','long_name','Radiosonde altitude');
    ncwriteatt(modf_name3,'alt_sonde','units','m');
    ncwriteatt(modf_name3,'alt_sonde','standard_name','altitude');
    ncwriteatt(modf_name3,'alt_sonde','missing_value',FillValue)
    ncwriteatt(modf_name3,'alt_sonde','actual_range',[min(MODF.radiosonde.alt_sonde) max(MODF.radiosonde.alt_sonde)]);
    ncwriteatt(modf_name3,'alt_sonde','instrument','');
    ncwriteatt(modf_name3,'alt_sonde','source','');
    alt = fillmissing(MODF.radiosonde.alt_sonde,'constant',FillValue);
    ncwrite(modf_name3,'alt_sonde',alt);    
end

if any(cellfun(@any,strfind(fieldnames(MODF.radiosonde),'zg_sonde')))
    nccreate(modf_name3,'zg_sonde','DIMENSIONS',{'sample'})
    ncwriteatt(modf_name3,'zg_sonde','long_name','Radiosonde geopotential height');
    ncwriteatt(modf_name3,'zg_sonde','units','m');
    ncwriteatt(modf_name3,'zg_sonde','standard_name','geopotential_height')
    ncwriteatt(modf_name3,'zg_sonde','missing_value',FillValue)
    ncwriteatt(modf_name3,'zg_sonde','actual_range',[min(MODF.radiosonde.zg_sonde) max(MODF.radiosonde.zg_sonde)])
    ncwriteatt(modf_name3,'zg_sonde','instrument','');
    ncwriteatt(modf_name3,'zg_sonde','source','');
    zg = fillmissing(MODF.radiosonde.zg_sonde,'constant',FillValue);
    ncwrite(modf_name3,'zg_sonde',zg);
end


fn3 = fieldnames(MODF.trajectoryProfile);

for i = 2:n3
    if any(strcmp(fn3,TrPrTb{1,i}))
        nccreate(modf_name3,TrPrTb{1,i},'DIMENSIONS',{'sample'})
        for j = 2:n
            if strcmp(TrPrTb{j,1},'actual_range')
                ncwriteatt(modf_name3,TrPrTb{1,i},TrPrTb{j,1},[min(MODF.trajectoryProfile.(TrPrTb{1,i})) max(MODF.trajectoryProfile.(TrPrTb{1,i}))]);
            elseif strcmp(TrPrTb{j,1},'missing_value')
                ncwriteatt(modf_name3,TrPrTb{1,i},TrPrTb{j,1},FillValue);
            elseif any(~ismissing(TrPrTb{j,i}))
                ncwriteatt(modf_name3,TrPrTb{1,i},TrPrTb{j,1},TrPrTb{j,i});
            end
            
        end
        MODF.trajectoryProfile.(TrPrTb{1,i}) = fillmissing(MODF.trajectoryProfile.(TrPrTb{1,i}),'constant',FillValue);
        ncwrite(modf_name3,TrPrTb{1,i},MODF.trajectoryProfile.(TrPrTb{1,i}));
    else
        disp(['missing ' TrPrTb{1,i}])
    end
end