function scs = read_scs_DYNAMO_2011(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  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
                  N = number of columns of data in raw file
                  ddd = julian date
                  hhh = hour
%}

%% open file and read
disp(['Reading scs file for hour ',int2str(hhh)]);
flist = fopen(dfl,'r');
N=63; % columns of data in raw file, SCS is 15s for DYNAMO2011 (Revelle)
temp = {}; scs_raw = [];
while ~feof(flist)
    try
        temp = textscan(flist,'%f %*s %f %f %f %*s %f %*s %f %*s %f %*s %f %*s %f %*s %f %*s %f %*s %f %*s %f %*s %f %*s %f %*s %f %*s %f %*s %f %*s %f %*s %f %*s %f %*s %f %*s %f %*s %f %*s %f %*s %f %*s %f %*s %f %*s %f %*s %f %*s %f %*s %f %*s %f %*s %f %*s %f %*s %f %*s %f %*s %f %*s %f %*s %f %*s %f %*s %f %*s %f %*s %f %*s %f %*s %f %*s %f %*s %f %*s %f %*s %f %*s %f %*s %f %*s %f %*s %f %*s %f %*s %f %*s %f %*s %f %*s %f %*s %f %*s %f %*s %f %*s %f  %*[^\n]','delimiter', ',', 'headerlines', 1,  '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;
scs_raw(2,:)=[];
fclose(flist);

%% define raw variables

gpshour = fix(scs_raw(2,:)/1e4);
gpsmin = fix(scs_raw(1,:)/1e5);
gpssec = (scs_raw(1,:) - gpsmin*1e5)/1e3;
gpstim = gpshour + gpsmin/60 + gpssec/60/60;      %Decimal time in hours
jd_scs = gpstim/24 + ddd;

lat =     scs_raw(56,:);
lon =     scs_raw(57,:);
cogm =    scs_raw(59,:);
sogm =    scs_raw(60,:)*0.514;   % converted from kt to m/s
tsgm =    scs_raw(24,:);
tssm =    scs_raw(26,:);
imrh =    scs_raw(7,:);
imta =    scs_raw(3,:);
imsol =   scs_raw(14,:);
orgm =    scs_raw(10,:);
lrg =     scs_raw(50,:);
%odec =    scs_raw(15,:)*.514; % not working for DYNAMO
imum =    scs_raw(18,:);
imdm =    scs_raw(19,:);
imir  =   scs_raw(13,:);
imPress = scs_raw(6,:);       % ships pressure
imqa = qair_p([imta' imrh' imPress']); % specific humidity

% clean up lat/lon
lon(lon<PosLims(1)) = NaN;
lon(lon>PosLims(2)) = NaN;
lat(lat<PosLims(3)) = NaN;
lat(lat>PosLims(4)) = NaN;
lat = despike2(lat); % despike & replace NaNs
lon = despike2(lon);

%% 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','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

