function scs = read_scs_HiWinGS_2013(dfl,ddd,hhh,ta_scs_adj,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  imrain      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
18  imrspd      relative wind speed, m/s
19  imrdir      relative wind dir, +/- 180 deg from bow

input parameters: dfl = string path to scs file
                  N = number of columns of data in raw file
                  ddd = julian date
                  hhh = hour

File Format:  approx 0.1Hz
IDX     TAG    Date        Time_GMT  Dec_Lat    Dec_Lon     SPD     HDT    COG    SOG   WXTP_Ta WXTS_Ta WXTP_Pa WXTS_Pa WXTP_Ri WXTS_Ri WXTP_Rc WXTS_Rc WXTP_Dm WXTS_Dm WXTP_Sm WXTS_Sm WXTP_RH WXTS_RH WXTP_TS WXTS_TS WXTP_TD WXTS_TD SWR    LWR   ND   SBE45S  SBE48T  FLR       SSV      Depth12 Depth35 
3031875 $dsLog 2013/10/15  00:30:31  58.380150  -45.270833  000.09  052.7  333.2  1.53  7.6     7.4     1009.3  1009.2  0.0     0.0     13.07   16.27   004     359     19.6    19.4    81.7    84.4    19.53   19.26   58.99   53.99   15.652 345.6 NaN  34.9598 7.0424  00095.50  1480.41  00.0    00.00
1 2 3          4    5  6   7  8  9   10         11          12      13     14     15   16       17      18      19      20      21      22      23      24      25      26      27      28      29      30      31      32      33      34     35    36   37      38      39         40      41      42

%}

%% open file and read
disp(['Reading scs file for hour ',int2str(hhh)]);
flist = fopen(dfl,'r');
N = 43; % columns of data in raw file
temp = []; scs_raw = [];
while ~feof(flist)
    try
        temp = textscan(flist,'%2f%2f%3f %*6c %4f/%2f/%2f %2f:%2f:%2f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %*[^\n]','delimiter', ', ', 'headerlines', 2,  'emptyvalue', NaN);
        scs_raw = [scs_raw cell2mat(temp)'];
    catch
        [lnr2,cnr2] = 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});
            scs_raw = [scs_raw cell2mat(temp)'];
        end;
    end;
end;
fclose(flist);

%% define raw variables
jd_scs =  ddd + (hhh + (scs_raw(1,:)+(scs_raw(2,:)+scs_raw(3,:)/1000)/60)/60)/24;
lat =     scs_raw(10,:);        % poditive = N, negative = S
lon =     scs_raw(11,:);        % positive = E, negative = W
lon(lon<0) = lon(lon<0)+360;    % all positive to E
odec =    scs_raw(12,:)*NaN;    % ship speed through water, kts, NOT WORKING FOR HIWINGS
lrg =     scs_raw(13,:);        % heading
cogm =    scs_raw(14,:);        % course over ground, deg, C-Nav GPS receiver
sogm =    scs_raw(15,:)*0.514;	% speed over ground, m/s, C-Nav GPS receiver
WXTP_Ta = scs_raw(16,:);        % Ta, port WXT, C
WXTS_Ta = scs_raw(17,:);        % Ta, stbd WXT, C
WXTP_Pa = scs_raw(18,:);        % P, port WXT, mb
WXTS_Pa = scs_raw(19,:);        % P, stbd WXT, mb
WXTP_Ri = scs_raw(20,:);        % rain rate, port WXT, mm/h
WXTS_Ri = scs_raw(21,:);        % rain rate, stbd WXT, mm/h
WXTP_Rc = scs_raw(22,:);        % cumulative rain, port WXT, mm
WXTS_Rc = scs_raw(23,:);        % cumulative rain, stbd WXT, mm
WXTP_Dm = scs_raw(24,:);        % rwdir, port WXT, deg
WXTP_Dm(WXTP_Dm>180) = WXTP_Dm(WXTP_Dm>180)-360; % +/- 180 deg format
WXTS_Dm  =scs_raw(25,:);        % rwdir, stbd WXT, deg
WXTS_Dm(WXTS_Dm>180) = WXTS_Dm(WXTS_Dm>180)-360; % +/- 180 deg format
WXTP_Sm = scs_raw(26,:);        % rwspd, port WXT, m/s
WXTS_Sm = scs_raw(27,:);        % rwspd, stbd WXT, m/s
WXTP_RH = scs_raw(28,:);        % RH, port WXT, %
WXTS_RH = scs_raw(29,:);        % RH, stbd WXT, %
WXTP_TS = scs_raw(30,:);        % true wspd, port WXT, m/s, corrected for heading, SOG, & COG
WXTS_TS = scs_raw(31,:);        % true wspd, stbd WXT, m/s, corrected for heading, SOG, & COG
WXTP_TD = scs_raw(32,:);        % true wdir, port WXT, m/s, corrected for heading, SOG, & COG
WXTS_TD = scs_raw(33,:);        % true wdir, stbd WXT, m/s, corrected for heading, SOG, & COG
imir    = scs_raw(35,:);        % pir, IR radiation, W/m2
imsol   = scs_raw(34,:);        % psp, solar radiation, W/m2
tssm    = scs_raw(37,:);        % salinity, ppth, SBE45
tsgm    = scs_raw(38,:);        % seawater temperature, C, SBE45
flor    = scs_raw(39,:);        % fluorescence, mV, WetLabs Wet-StarWS3S-901P
imqaP = qair_p([WXTP_Ta',WXTP_RH',WXTP_Pa']); % specific humidity, port WXT, g/kg
imqaS = qair_p([WXTS_Ta',WXTS_RH',WXTS_Pa']); % specific humidity, stbd WXT, g/kg

%% check for unreasonable values
lon(lon<PosLims(1)) = NaN;
lon(lon>PosLims(2)) = NaN;
lat(lat<PosLims(3)) = NaN;
lat(lat>PosLims(4)) = NaN;
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;
imqaP(imqaP>25) = NaN; imqaS(imqaS>25) = NaN;
imqaP(imqaP<0.5) = NaN; imqaS(imqaS<0.5) = NaN;
WXTP_RH(WXTP_RH>100) = 100; WXTS_RH(WXTS_RH>100) = 100;
WXTP_RH(WXTP_RH<0) = NaN; WXTS_RH(WXTS_RH<0) = NaN;
WXTP_Ri(WXTP_Ri<0) = 0; WXTS_Ri(WXTS_Ri<0) = 0;

%% despike and replace NaNs with nearest neighbor
[latm,~] = despike2(lat);
[lonm,~] = despike2(lon);
cogm = unwrap(cogm*pi/180); % unwrap and convert to radians
[cogm,~] = despike2(cogm);
[sogm,~] = despike2(sogm);
% [tsgm,~] = despike2(tsgm);
% [tssm,~] = despike2(tssm);
imrh = (WXTP_RH + WXTS_RH)/2;
% [imrh,~] = despike2(imrh);
imta = (WXTP_Ta + WXTS_Ta)/2;
% [imta,~] = despike2(imta);
% [imsol,~] = despike2(imsol);
% [imir,~] = despike2(imir);
lrg = unwrap(lrg*pi/180); % unwrap and convert to radians
%[lrg,~] = despike2(lrg);
[WXTP_Dm,~] = despike2(WXTP_Dm);
ii = find(WXTP_Dm<=0); jj = find(WXTP_Dm>0);
imum = NaN(size(WXTP_Dm));
imum(ii) = WXTP_TS(ii); % select true wind based on rwdir
imum(jj) = WXTS_TS(jj);
% [imum,~] = despike2(imum);
imdm = NaN(size(WXTP_Dm));
imdm(ii) = WXTP_TD(ii);
imdm(jj) = WXTS_TD(jj);
imdm = unwrap(imdm*pi/180); % unwrap and convert to radians
% [imdm,~] = despike2(imdm);
imPress = (WXTP_Pa + WXTS_Pa)/2;
% [imPress,~] = despike2(imPress);
imqa = (imqaP + imqaS)/2;
% [imqa,~] = despike2(imqa);
imrspd = NaN(size(WXTP_Sm));
imrspd(ii) = WXTP_Sm(ii);
imrspd(jj) = WXTS_Sm(jj);
% [imrspd,~] = despike2(imrspd);
imrdir = NaN(size(WXTP_Dm));
imrdir(ii) = WXTP_Dm(ii);
imrdir(jj) = WXTS_Dm(jj);
% [imrdir,~] = despike2(imrdir);
imrain = (WXTP_Ri + WXTS_Ri)/2;

%% rebuild raw array
scs_raw = [jd_scs' latm' lonm' cogm' sogm' tsgm' tssm' imrh' imta' imsol' imrain' lrg' imum' imdm' imir' imPress' imqa imrspd' imrdir'];

% 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:end),jd_ref','linear','extrap');
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
scs(:,9) = scs(:,9) + ta_scs_adj;               % adjust air T if necessary

