function scs = read_scs_STRATUS_2017(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   cogm        gps course over ground, deg
5   sogm        gps speed over ground, m/s
6   tsgm        tsg sst from ship intake
7   tssm        tsg salinity
8   imrh        RH
9   imta        Air temp, C
10  imsol       psp radiation, W/m2
11  orgm        rain, mm/hr
12  lrg         laser ring gyro heading, deg
13  imum        true wind speed, m/s
14  imdm        true wind directin, deg
15  imir        pir radiation, W/m2
16  imPress     air pressure, mb
17  imqa        air specific humidity, g/kg

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=19; % columns of data in raw file
temp = []; scs_raw = [];
while ~feof(flist)
    try
        temp = textscan(flist,'%2f%2f%3f %f%c %f%c %f %f %f %f %f %f %f %f %f %f %f %f %*[^\n]','delimiter',',','headerlines',1,'emptyvalue',NaN);
        temp{1,5} = double(temp{1,5}); % convert N or S to ascii code integer
        temp{1,7} = double(temp{1,7}); % convert E or W to ascii code integer
        scs_raw = [scs_raw cell2mat(temp)'];
    catch
        [lnr2,~] = size(temp{1,N});
        for ll=1:N
            if length(temp{1,ll})~=lnr2
                temp{1,ll}(lnr2+1)=[];
            end
        end
        if ~isempty(temp{1,1})
            temp{1,5} = double(temp{1,5}); % convert N or S to ascii code integer
            temp{1,7} = double(temp{1,7}); % convert E or W to ascii code integer
            scs_raw = [scs_raw cell2mat(temp)'];
        end
    end
end
fclose(flist);

%% define raw variables

[~,cols] = size(scs_raw);
jd_scs =  ddd + (hhh + (scs_raw(1,:)+(scs_raw(2,:)+scs_raw(3,:)/1000)/60)/60)/24;
lat =     scs_raw(4,:)./100;
latns =   scs_raw(5,:);         % 78=N, 87=W, 83=S, 69=E
lon =     scs_raw(6,:)./100;
lonew =   scs_raw(7,:);         % 78=N, 87=W, 83=S, 69=E
cogm =    scs_raw(8,:);
sogm =    scs_raw(9,:)*0.514;   % converted from kt to m/s
tsgm =    scs_raw(10,:);
tssm =    scs_raw(11,:);
imrh =    scs_raw(12,:);
imta =    scs_raw(13,:);
imsol =   scs_raw(14,:);
% orgm =    scs_raw(15,:);  % not available for STRATUS
orgm = NaN(1,cols);
lrg =     scs_raw(15,:);
%odec =    stx(16,:)*.514; % not available for STRATUS
imum =    scs_raw(16,:)*.514;  % converted from kt to m/s
imdm =    scs_raw(17,:);
imir  =   scs_raw(18,:);
imPress = scs_raw(19,:);       % ships pressure
imqa = qair_p([imta' imrh' imPress']); % specific humidity

%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
scs_raw = [jd_scs' lat' lon' cogm' sogm' tsgm' tssm' imrh' imta' imsol' orgm' lrg' imum' imdm' imir' imPress' imqa];

% interpolate to exactly 1 Hz
start = (ddd+hhh/24);
delta = double(1.0/86400);
last = start + 3599*delta;
jd_ref = [start:delta:last];
scs = interp1(jd_scs',scs_raw(:,2:17),jd_ref','linear');
scs = [jd_ref' scs];
scs(:,4) = mod(scs(:,4)*180/pi + 360, 360);     % cogm back to degrees
scs(:,12) = mod(scs(:,12)*180/pi + 360, 360);   % lrg back to degrees
scs(:,14) = mod(scs(:,14)*180/pi + 360, 360);   % imdm back to degrees

