function scs = read_scs_WHOTS_2022(dfl,ddd,hh,PosLims)
%{
    reads 1Hz scs file specified by dfl and returns 10-min table with the
    following variables:

    1   jd_ref      1 Hz timestamp
    2   lat         latitude deg, N positive
    3   lon         longitude deg, E positive
    4   sogm        gps speed over ground, m/s
    5   cogm        gps course over ground, deg
    6   lrg         heading, deg
    7   imta        Air temp, C
    8   imrh        RH humidity, g/kg
    9   imPress     air pressure, mb
    10  imqa        specific humidity, g/kg
    11  imum        true wind speed, m/s
    12  imdm        true wind direction, deg
    13  imrum       relative wind speed, m/s
    14  imrdm       relative wind direction, deg

    input parameters: dfl = string path to scs file
                      ddd = julian date
                      hhh = hour
                      PosLims = lat/lon limits for filtering position data
%}

%% open file and read
disp(['Reading scs file for hour ',int2str(hh)]);
flist = fopen(dfl,'r');
N = 16; % columns of data in raw file
tempscs = {};
scsm = [];
%%read file
while feof(flist)==0
    try
        % temp will be 1x 16 cell array, cells will be 3600 x 1 arrays
        %         07/21/2022 21:00 Sette SCS met and nav (SN: 999)
        % 0000409 2022/07/21,21:00:01,2121.9621,N,15757.9760,W,0.003,79.6,53.3,27.3,59,1018.02,73.3,20,16.4,16.4
        % 0001411 2022/07/21,21:00:02,2121.9621,N,15757.9760,W,0.003,79.6,53.3,27.3,59,1018.02,78.3,25,16.8,16.8
        tempscs = textscan(flist,'%2f%2f%3f %*f/%*f/%*f %*f:%*f:%*f %f%c %f%c %f %f %f %f %f %f %f %f %f %f %*[^\n]',...
            'delimiter',',','headerlines',1,'treatasempty',{'?????','-Unknown-'});
        tempscs{1,5} = double(tempscs{1,5}); % convert N or S to ascii code integer
        tempscs{1,7} = double(tempscs{1,7}); % convert E or W to ascii code integer
        scsm = [scsm cell2mat(tempscs)'];
    catch
        %[lnr2,~] = size(temp{1,N});
        for ll=1:N
            if length(tempscs{1,ll})~=length(tempscs{1,N})
                tempscs{1,ll}(length(tempscs{1,N})+1)=[];
            end
        end
        if ~isempty(tempscs{1,1})
            tempscs{1,5} = double(tempscs{1,5}); % convert N or S to ascii code integer
            tempscs{1,7} = double(tempscs{1,7}); % convert E or W to ascii code integer
            scsm = [scsm cell2mat(tempscs)'];
        end
    end
end
fclose(flist);

ii = isinf(scsm); scsm(ii) = NaN;

%% define raw variables
jd_scs =  ddd + (hh + (scsm(1,:)+(scsm(2,:)+scsm(3,:)/1000)/60)/60)/24;
lat =     scsm(4,:)./100;
latns =   scsm(5,:);          % 78=N, 87=W, 83=S, 69=E
lon =     scsm(6,:)./100;
lonew =   scsm(7,:);          % 78=N, 87=W, 83=S, 69=E
sogm =    scsm(8,:)*0.514;    % converted from kt to m/s
cogm =    scsm(9,:);
lrg =     scsm(10,:);         % heading
imta =    scsm(11,:);
imrh =    scsm(12,:);
imPress = scsm(13,:);         % scs,ships pressure
imqa = qair_p([imta' imrh' imPress']); % specific humidity
imdm =    scsm(14,:);         % true wind direction scs,deg
imum =    scsm(16,:)*0.514;   % true wind speed scs,m/s
imrdm =   scsm(15,:);         % rel wind direction scs,deg
imrum =   scsm(17,:)*0.514;   % rel wind speed scs,m/s

%convert lat/lon from degs.minutes to dec degrees and clean up
latmin = (lat-floor(lat))*100;
lat2 = floor(lat)+latmin/60;
lat2(latns==83) = -lat2(latns==83);
lonmin = (lon-floor(lon))*100;
lon2 = floor(lon)+lonmin/60;
lon2(lonew==87) = -lon2(lonew==87)+360;
lon2(lon2<PosLims(1)) = NaN;
lon2(lon2>PosLims(2)) = NaN;
lat2(lat2<PosLims(3)) = NaN;
lat2(lat2>PosLims(4)) = NaN;
lat = despike2(lat2); % despike & replace NaNs
lon = despike2(lon2);

%% check for unreasonable values
%imsol(imsol>2000) = NaN;
%imsol(imsol<-20) = NaN;
%imir(imir>1500) = NaN;
%imir(imir<-20) = NaN;
%tssm(tssm>40) = NaN;
%tssm(tssm<25) = NaN;
%tsgm(tsgm>50) = NaN;
%tsgm(tsgm<-2) = NaN;
cogm(cogm>360) = NaN;
cogm(cogm<0) = NaN;
sogm(sogm>20) = NaN;
imdm(imdm<0) = NaN;
imdm(imdm>360) = NaN;
imum(imum<0) = NaN;
imum(imum>50) = NaN;
imrdm(imrdm<0) = NaN;
imrdm(imrdm>360) = NaN;
imrum(imrum<0) = NaN;
imrum(imrum>50) = NaN;
imta(imta>50) = NaN;
imqa(imqa>25) = NaN;
imqa(imqa<0.5) = NaN;
imrh(imrh>100) = 100;
imrh(imrh<0) = NaN;

%% despike and replace NaNs with nearest neighbor
cogm = unwrap(cogm*pi/180); % unwrap and convert to radians
[cogm,~] = despike2(cogm);
[sogm,~] = despike2(sogm);
%[tsgm,~] = despike2(tsgm);
%[tssm,~] = despike2(tssm);
[imrh,~] = despike2(imrh);
[imta,~] = despike2(imta);
%[imsol,~] = despike2(imsol);
%[imir,~] = despike2(imir);
lrg = unwrap(lrg*pi/180); % unwrap and convert to radians
[lrg,~] = despike2(lrg);
[imum,~] = despike2(imum);
imdm = unwrap(imdm*pi/180); % unwrap and convert to radians
[imdm,~] = despike2(imdm);
imrdm = unwrap(imrdm*pi/180); % unwrap and convert to radians
[imrdm,~] = despike2(imrdm);
[imPress,~] = despike2(imPress);
[imqa,~] = despike2(imqa);

%% rebuild raw array in output format
scsm = [jd_scs',lat',lon',sogm',cogm',lrg',imta',imrh',imPress',imqa,imum',imdm',imrum',imrdm'];
[~,c]=size(scsm);

%% average to 10 min
start = (ddd+hh/24);

delta = double(10/1440);
last = start + 6*delta;
tref = start:delta:last;

scs10m = interval_avg(jd_scs',scsm(:,2:c),tref');
scs10m(:,5) = mod(scs10m(:,5)*180/pi + 360, 360);      % cogm back to degrees
scs10m(:,6) = mod(scs10m(:,6)*180/pi + 360, 360);      % lrg back to degrees
scs10m(:,12) = mod(scs10m(:,12)*180/pi + 360, 360);    % imdm back to degrees
scs10m(:,14) = mod(scs10m(:,14)*180/pi + 360, 360);    % imrdm back to degrees
ii = find(scs10m(:,14)>180); scs10m(ii,14) = scs10m(ii,14) - 360; % rwdir in +/-180° format

%% output table
vars = {'jd','lat','lon','sog','cog','hed','ta','rh','Pmb','qa','wspd','wdir','rwspd','rwdir'};
scs = array2table(scs10m,'VariableNames',vars);

end

