%------ read_hp_HiWinGS_2013.m
%	reads slow bulk met data files from PSD data system

%	These variables must exist in the current workspace
%		hhh: numerical value of current hour (00-23)
%		ddd: numerical value of current day of year (0-365)
%		yyyy: string variable of current year (e.g.'2014')
%		way_raw_data_flux: string path to data folder
%       clearit: sets clear flag
%		
%	These arrays & variables output:
%		press, psp, pir, org, Tsea, T, RH, Td, jd_trh
%		Rhvais, Tvais, jd_pc
%       

%% initialize
disp('read_hp_hiwings2013.m')
fclose('all');

g=0.0098;					% adiabatic lapse rate
sig_sb=5.67e-8;				% Stefan Boltzmann constant
jd=sprintf('%03i',ddd);     % day of year string
hr=sprintf('%02i',hhh);     % hour string

% Full path to file on Windows
% e=[way_raw_data_flux, yyyy(3:4),jd,'\met0' yyyy(3:4) jd hr '_raw.txt'];

% Full path to file on Mac
dfl=[way_raw_data_flux, yyyy(3:4),jd,'/met0' yyyy(3:4) jd hr '_raw.txt'];

disp(['Reading means file from hour ',hr]);
flist=fopen(dfl,'r');

% Declare Arrays
nr1=61;                 % number of variables for HiWinGS
jd_pc=ones(1,60)*NaN;   % empty array for 1 min time stamps

%% load raw text files
if flist>0,             % if the file exists, 
    stx=[];             % empty array for output
    temp={};            % empty cell array for textscan
    while ~feof(flist)  % read file, nr1 variables       
        try
            % temp is 1x(nr1) cell array containing ~60x1 arrays of data
            % Format string must have total of nr1 fields
            %                       1                               10                            20                            30                            40                            50                            60
            temp =textscan(flist,'%2f%2f%3f %2f/%2f/%4f %2f:%2f:%2f %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %*[^\n]','delimiter',', ','headerlines',1,'emptyvalue',NaN);
            % stx is 60x9 array of time stamp data from file
            stx=[stx temp{1} temp{2} temp{3} temp{4} temp{5} temp{6} temp{7} temp{8} temp{9}];
            % For remaining cells, scan for NaNs and blanks, format as
            % numbers, and append to stx
            for ii=10:length(temp)  % Look for NaNs and blanks in data fields
                temp{ii}(ismember(temp{ii},'NAN'))=cellstr('NaN');
                temp{ii}(ismember(temp{ii},''))=cellstr('NaN');
                stx=[stx str2num(char(temp{ii}))];  % 60xnr1 array
            end;
%           lud=find((ismember(temp{ii},''))); %sometimes a comma is
%           missing thus offsetting the reading and empty value are
%           found which crashes the code - LB - 10192013
%           will edit manually the file for now...
        catch
            % length of last cell reflects missing values from any/all fields
            [lnr1,cnr1]=size(temp{1,nr1});
            if lnr1==0 || cnr1==0; continue; end; % read next batch of lines if no data
            for ii=1:nr1
                if length(temp{1,ii})~=lnr1
                    temp{1,ii}(lnr1+1)=[];  %truncate length of all cells to length(temp{1,nr1})
                end;
            end;
            if isempty(temp{1,1});
                temp{1,nr1}=temp{1,nr1}';
            end;
            stx = [stx cell2mat(temp)'];      % need to handle stx as above here???
        end;   %end try/catch
    end; %end ~feof(flist)
    stx = stx';
      % PSD timestamp marks the end of each 1-min interval
%     jd_pc=ddd+(hhh+(stx(1,:)+(stx(2,:)+stx(3,:)/1000)/60)/60)/24;
    % Datalogger timestamp marks the start of each interval, which conforms
    % to the assumed convention.
    jd_pc = ddd + stx(7,:)/24 + stx(8,:)/1440 + stx(9,:)/86400;  % timestamp from met file
    fclose(flist);
else  % no data default
    disp(['Reading means file from hour ',hr]);
    disp(sprintf('\b<--File not found.'))
    stx=ones(nr1,60)*NaN;
    jd_pc=ddd + (hhh+(0:59)/60)/24;
end;

%% adjustments & postcalibrations for flux observations (done on ship)
tsnake_adjustment=0; %degC;
ta_etl_adjustment=0;
rh_etl_adjustment=0;
bp_etl_adjustment=0;

% TRH=load('C:\data\CRUISES\HiWinGS_2013\Knorr\flux\Processed\cleanedTRH');

t_case1_adjustment=0;
t_dome1_adjustment=0;
t_case2_adjustment=0;
t_dome2_adjustment=0;

press=stx(35,:)+bp_etl_adjustment;  %pressure (mbar)

psp1=-stx(25,:);                    %Eppley PSP bridge
psp2=-stx(27,:);                    %K&Z PSP bridge

org=stx(30,:);                      %mm/hr

Tsea=stx(34,:)+tsnake_adjustment;

Tc1=stx(13,:) +t_case1_adjustment;  %case temperature for PIR1 Eppley
Td1=stx(15,:) +t_dome1_adjustment;  %dome temperature for PIR1 Eppley
Tc2=stx(17,:) +t_case1_adjustment;  %case temperature for PIR2 Eppley
Td2=stx(19,:) +t_dome1_adjustment;  %dome temperature for PIR2 Eppley

therm1=-stx(21,:);                  % thermopile voltages
therm2=-stx(23,:);

%% read in vaisala T/RH
read_tairh_HiWinGS_2013;    % reads in 1Hz arrays: RH, T, Td, jd_trh

Rhvais=intavg(jd_trh,RH,jd_pc);
Rhvais=Rhvais';
Tvais=intavg(jd_trh,T,jd_pc);
Tvais=Tvais';

%% compute radiation
C2K=273.15;  
pir1=therm1+sig_sb*(Tc1+C2K).^4-3*sig_sb*((Td1+C2K).^4-(Tc1+C2K).^4); %adjust B coeff to 3 for leg1
pir2=therm2+sig_sb*(Tc2+C2K).^4-4*sig_sb*((Td2+C2K).^4-(Tc2+C2K).^4);

% PIR1 had some solar contamination or overheating. Will use PIR2 only
pir=pir2;
% pir=(pir1+pir2)/2;
psp=(psp1+psp2)/2;

%% cleanup
if length(Rhvais)<60,
   Rhvais=[Rhvais (1:(60-length(Rhvais))).*0+NaN];
end;
if length(pir)<60,
   pir=[pir (1:(60-length(pir))).*0+NaN];
end;
if length(psp)<60,
   psp=[psp (1:(60-length(psp))).*0+NaN];
end;
if length(org)<60,
   org=[org (1:(60-length(org))).*0+NaN];
end;
if length(Tvais)<60,
   Tvais=[Tvais (1:(60-length(Tvais))).*0+NaN];
end;
if length(Tsea)<60,
   Tsea=[Tsea (1:(60-length(Tsea))).*0+NaN];
end;
if length(press)<60,
   press=[press (1:(60-length(press))).*0+NaN];
end;

if length(Rhvais)>60,
   Rhvais=Rhvais(1:60);
end;
if length(pir)>60,
   pir=pir(1:60);
end;
if length(psp)>60,
   psp=psp(1:60);
end;
if length(org)>60,
   org=org(1:60);
end;
if length(Tvais)>60,
   Tvais=Tvais(1:60);
end;
if length(Tsea)>60,
   Tsea=Tsea(1:60);
end;
if length(press)>60,
   press=press(1:60);
end;

if clearit
    clear psp1 psp2 pir2 pir1 stx C2K dfl
    clear g flist ans therm1 therm2 st2 sig_sb Tc1 Tc2 Td1 Td2 nr1
    clear tsnake_adjustment ta_etl_adjustment rh_etl_adjustment bp_etl_adjustment
    clear t_case1_adjustment t_dome1_adjustment t_case2_adjustment t_dome2_adjustment
end;

