function scs = auto_read_scs_WHOTS_2019(dfl,ddd,hhh,PosLims)
%{
reads scs file specified by dfl and returns array with the
following columns:

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   lrg        heading, deg
6   imta        Air temp, C
7   imrh        RH humidity, g/kg
8   imPress     air pressure, mb
9   imdm        true wind direction, deg
10  rmdm        relative wind direction, deg
11  imum        true wind speed, m/s
12  rmum        relative wind speed, m/s
   not collected cogm        course over ground, 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(hhh)]);
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
        %10/05/2019 22:00 Generic with ASCII Output (SN: 123456)
        %0000290 2019/10/05,22:00:00,2237.2889N,15802.0438W,1.072,$HEHDT,062.9,T*22,026.3,069.8,1015.44,54.41,353.0,6.67,007.7
        %0001304 2019/10/05,22:00:01,2237.2888N,15802.0435W,1.053,$HEHDT,062.4,T*2F,026.3,069.8,1015.47,52.56,353.0,6.71,007.7
        
        tempscs = textscan(flist,'%2f%2f%3f %*f/%*f/%*f %*f:%*f:%*f %f%c %f%c %f %*6c %f %*4c %f %f %f %f %f %f %f %*[^\n]','delimiter',',','headerlines',1,'emptyvalue',NaN);
        
       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);

%% define raw variables
jd_scs =  ddd + (hhh + (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
lrg =     scsm(9,:);  %Heading
imta =    scsm(10,:);
imrh =    scsm(11,:);
imPress = scsm(12,:);       % scs,ships pressure
imqa = qair_p([imta' imrh' imPress']); % specific humidity
imdm =    scsm(13,:); %true wind direction scs,deg
imum =    scsm(15,:)*0.514;  % true wind speed scs,m/s

%cogm =    scsm(8,:);
%tsgm =    scsm(10,:);
%tssm =    scsm(11,:);
%imsol =   scsm(14,:);
% orgm =    scsm(15,:);  % not available for STRATUS
%orgm = NaN(1,cols);
%odec =    stx(16,:)*.514; % not available for STRATUS
%imir  =   scsm(18,:);


%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;
%imdm(imdm>360) = NaN;
%imdm(imdm<0) = NaN;
%imum(imum<0) = NaN;
%imum(imum>50) = NaN;
sogm(sogm>20) = NaN;
imta(imta>50) = NaN;
%orgm(orgm>100) = 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);
[imPress,~] = despike2(imPress);
[imqa,~] = despike2(imqa);

%% rebuild raw array
scsm = [jd_scs' lat' lon' sogm' imrh' imta' lrg' imum' imdm' imPress' imqa];
[l,c]=size(scsm);

% interpolate to exactly 1 Hz
% start = (ddd+hhh/24);
% delta = double(1.0/1440);
% last = start + 3599*delta;
% jd_ref = [start:delta:last];
% scs = interp1(jd_scs',scsm(:,2:c),jd_ref','linear');
% scs = [jd_ref' scs];

%% average to 10 min
start = (ddd+hhh/24);
% delta = double(1/1440);
% last = start + 60*delta;
delta = double(10/1440);
last = start + 6*delta;
tref = start:delta:last;
% scs = interp1(jd_scs',scsm(:,2:N),tref','linear');
% scs = [tref' scs];
scs = interval_avg(jd_scs',scsm(:,2:c),tref');
%scs(:,4) = mod(scs(:,4)*180/pi + 360, 360);     % cogm back to degrees
scs(:,7) = mod(scs(:,7)*180/pi + 360, 360);   % lrg back to degrees
scs(:,9) = mod(scs(:,9)*180/pi + 360, 360);   % imdm back to degrees

