%------ run_Calwater_2015.m

%	Runs motion correction and flux calculations for Calwater
%   Based on LB 062014 scripts but edited & commented by BWB
%   Uses read_gprm for reformatted raw gps files - much faster
%   than reading gps files.  See comments in read_gprm.m

%   NOTE: RAW SON DATA FILES MUST BE CLEANED BEFORE THESE
%   SCRIPTS WILL RUN.  T/RH MIGHT NEED CLEANING TOO.
%   SEE COMMENTS IN read_son AND read_tairh SCRIPTS.

%% initialize
close all;
fclose all;
clear;

datasetv = 'v1';
% datasetv = 'v2';  % use this for tilt correction prior to ship
                    % speed correction

read_parameters_Calwater_2015_mac;   % set paths and other run variables
currentMname = 'Calwater_2015';
graphformat = '.png';  %select graphics format files
graphdevice = '-dpng'; %select graphic device

clearit = 0;
printit = 0;
yyyy = '2015';

decorr = 0;         % set to decorrelate winds & motion...
if decorr==1
    dcname = 'decorr';
else
    dcname = 'std';
end;

%% Run daily flux calcs
% edit start jd and end jd as desired
stjd = 24;
endjd = 24;
badSonPnts = [];       % arrays to keep track of bad sonic data
missingPnts = [];
jday_badPnts = [];

for ddd=stjd:endjd
    [m,d] = yd2md(str2num(yyyy), ddd);
    Vdate = [str2num(yyyy) m d];
    jd = sprintf('%03i',ddd);
    for hhh=12; %hhh=00:24;
        % check for data, if not continue to next iteration
        % if scs file is missing, don't run for this hour
        hr = sprintf('%02i',hhh);
        dfl = [way_raw_data_flux,yyyy(3:4),jd,'/scs0',yyyy(3:4),jd,hr,'_raw.txt'];

        if exist(dfl,'file')
            read_scs_Calwater_2015_mac;       % read scs data...
            read_met_Calwater_2015_mac;       % Read slow met file
            read_motion_Calwater_2015_mac;    % Read MotionPak file
            read_son_Calwater_2015_mac;       % Read sonic file
            read_gprm_Calwater_2015_mac;      % Read reformatted Hemisphere gps file

            % save date/time & # bad pnts from read_son
            badSonPnts = [badSonPnts; badPnts'];
            missingPnts = [missingPnts; missing'];
            jday_badPnts = [jday_badPnts; badTime'];
            clear badPnts missing badTime;

            tref = 0:1/10:35999/10;                 % 10 Hz for 1 hr
            jd_ref = ddd + hhh/24 + tref./86400;    % ref timestamp as jd

            % Get heading data and convert to radians
            % Using 1 sec ship data for slow heading right now
            % Could also use (filtered?) PSD heading...
            rdcon = pi/180;
            lrg2 = replace_NaN_nearest_neighbor(lrg);
            lrgc = cos(lrg2.*rdcon);        % convert to N
            lrgs = sin(lrg2.*rdcon);        % convert to E
            % interpolate heading to 10 Hz for motion correction
            lrgcm = interp1(jd_scs,lrgc,jd_ref);
            lrgsm = interp1(jd_scs,lrgs,jd_ref);
            % filter out NaNs in lrgsm & lrgcm following interp1
            lrgcm2 = replace_NaN_nearest_neighbor(lrgcm);
            lrgsm2 = replace_NaN_nearest_neighbor(lrgsm);
            lrgcm = lrgcm2;
            lrgsm = lrgsm2;
            clear lrgsm2 lrgcm2

            heading = mod(atan2(lrgsm,lrgcm)./rdcon + 360, 360);  % add 360 so all are positive
            shead = sin(heading./180.*pi);
            chead = cos(heading./180.*pi);
            heading = heading/180*pi;       % convert to radians

            % read licor file
            read_co2_Calwater_2015_mac;         % read licor co2 & h2o data

            % run motion corrections and raw flux calculations
            if strcmp(datasetv,'v1')
                nau_correct_Calwater_2015;     % motion corrections and EC flux calculations
                nau_process_Calwater_2015;     % inertial dissipation, compute means, run bulk model, output da file
            elseif strcmp(datasetv,'v2')
                nau_correct_Calwater_2015_v2_mac;  % version with streamlinerotation before ship speed correction
                nau_process_Calwater_2015_v2;
            else
                error('data set version unrecognized')
            end;
            close all;
        end;
    end;
end;

%% concatenate hourly files to daily and cruise total

% for ddd=stjd:endjd
%     jd=sprintf('%03i',ddd);
%     filecat([way_proc_data_flux,'procturb/da/da_',datasetv,'_',jd,'*',dcname,'.txt'],[way_proc_data_flux,'procturb/daily/',cruise,'_da_',datasetv,'_',jd,'_',dcname,'.txt']);
%     filecat([way_proc_data_flux,'procturb/cr/cr_',datasetv,'_',jd,'*',dcname,'.txt'],[way_proc_data_flux,'procturb/daily/',cruise,'_cr_',datasetv,'_',jd,'_',dcname,'.txt']);
%     filecat([way_proc_data_flux,'procturb/co/co_',datasetv,'_',jd,'*',dcname,'.txt'],[way_proc_data_flux,'procturb/daily/',cruise,'_co_',datasetv,'_',jd,'_',dcname,'.txt']);
%     filecat([way_proc_data_flux,'procturb/sp/sp_',datasetv,'_',jd,'*',dcname,'.txt'],[way_proc_data_flux,'procturb/daily/',cruise,'_sp_',datasetv,'_',jd,'_',dcname,'.txt']);
% end;
% filecat([way_proc_data_flux,'procturb/daily/',cruise,'_da_',datasetv,'_*_',dcname,'.txt'],[way_proc_data_flux,'procturb/total/',cruise,'_da_',datasetv,'_',dcname,'.txt']);
% filecat([way_proc_data_flux,'procturb/daily/',cruise,'_cr_',datasetv,'_*_',dcname,'.txt'],[way_proc_data_flux,'procturb/total/',cruise,'_cr_',datasetv,'_',dcname,'.txt']);
% filecat([way_proc_data_flux,'procturb/daily/',cruise,'_co_',datasetv,'_*_',dcname,'.txt'],[way_proc_data_flux,'procturb/total/',cruise,'_co_',datasetv,'_',dcname,'.txt']);
% filecat([way_proc_data_flux,'procturb/daily/',cruise,'_sp_',datasetv,'_*_',dcname,'.txt'],[way_proc_data_flux,'procturb/total/',cruise,'_sp_',datasetv,'_',dcname,'.txt']);


