% Comparison plots and adjustments to variables for the final bulk flux run

% Edit this script to make adjustments and corrections to the met data

% You don't need to run this script directly.  It is called from da_red1.m
% but, you can edit the lines below to make adjustments to the data or select
% specific variables to be used in the final COARE model run.

disp('HiWinGS_2013_fixit_da1')

% overlap indices for arrays x and MM10
xxx = 2881:7776; % x indices jd 283 through 316
mmm = 96:4991;   % MM10 indices jd 283 through 316

%% WS
figure; plot(jdy,U,'r.',jdy,U_scs,'b.',jdy,U_uh,'g.',MM10(:,1),MM10(:,9),'k.');
xlim([stdt endt]); ylim([0 30]);grid;
title('True wind speed, ship(b), PSD(r), UH(g), MM10(k)');
u_best = U_uh;                   % uh sonic generally better
u_best(jdy<260) = U(jdy<260);    % except prior to jd 260
% and a few instances when uh sonic is down but psd sonic is running
u_best(jdy>260 & isnan(u_best)) = U(jdy>260 & isnan(u_best));
% and finally, substitute HiWinGS period for MM10 result for U relative to water surface
u_best(xxx) = MM10(mmm,11);

urel_best = relspd_uh;                 % as above for rel wind spd
urel_best(jdy<260) = urel(jdy<260);
urel_best(jdy>260 & isnan(urel_best)) = urel(jdy>260 & isnan(urel_best));

%% WD
figure; plot(jdy,dir,'r.',jdy,dir_scs,'b.',jdy,dir_uh,'g.');
xlim([stdt endt]); grid;
title('True wind direction, ship(b), PSD(r), UH(g)');
dir_best = U_uh;                    % uh sonic generally better
dir_best(jdy<260) = dir(jdy<260);   % except prior to jd 260
% and a few instances when uh sonic is down but psd sonic is running
dir_best(jdy>260 & isnan(dir_best)) = dir(jdy>260 & isnan(dir_best));
dir_best(xxx) = MM10(mmm,10);

% same for relative wind direction
reldir_best = reldir_uh;                  % uh sonic generally better
reldir_best(jdy<260) = reldir(jdy<260);   % except prior to jd 260
% and a few instances when uh sonic is down but psd sonic is running
reldir_best(jdy>260 & isnan(reldir_best)) = reldir(jdy>260 & isnan(reldir_best));
reldir_best(xxx) = MM10(mmm,13);
reldir_best(reldir_best>180) = reldir_best(reldir_best>180)-360;

figure; plot(jdy,reldir_best,'b.',[jdy(1),endt],[90,90],'r-',[jdy(1),endt],[-90,-90],'r-');
xlim([stdt endt]); grid; title('Relative wind direction - best');

%% Pressure
P_im(P_im<950) = NaN;
press(press<950) = NaN; press = press + 0.12*zP; % PSD corrected to sea level from 10m
for ii=1:144:length(x)  % remove pressure spikes and fill in NaNs each day
    press(ii:ii+143) = despike2(press(ii:ii+143));
    P_im(ii:ii+143) = despike2(P_im(ii:ii+143));
end

% MM10 pressure corrected to sea level...
figure; plot(jdy,press,'r.',jdy,P_im,'b.',MM10(:,1),MM10(:,14)+0.12*15,'k.');
title('Sea Level Pressure (mb), ship(b), PSD(r), MM10(k)');
xlim([stdt endt]); grid;
P_best = P_im;  % choose avg ship bp for bulk model
P_best(xxx) = MM10(mmm,14)+0.12*15; % but use MM10 for HiWinGS

%% Tilt angle - flow distortion angle derived from the motion correction script
% this is tilt angle prior to ship speed correction
figure; plot(reldir,tiltx,'r.',reldir_uh,tilt_uh,'g.',[-90 -90],[-5 25],'b',[90 90],[-5 25],'b');
axis([-180 180 -5 25]); grid; title('Tilt Angle, degrees, red=PSD, green=UH');
clear reldir2
tilt_best = tilt_uh;                   % same as for wind speed
tilt_best(jdy<260) = tiltx(jdy<260);   % except prior to jd 260
% and a few instances when uh sonic is down but psd sonic is running
tilt_best(jdy>260 & isnan(tilt_best)) = tiltx(jdy>260 & isnan(tilt_best));

%% Tsea
tsnk(tsnk>25) = NaN; tsg(tsg>25)=NaN;
figure;plot(jdy,tsnk,'r.-',jdy,tsg,'b.',MM10(:,1),MM10(:,18),'k.'); grid;
title(' SST(degC), ship(b), PSD(r), MM10(k)');
xlim([stdt endt]);
ts_best = tsg;  % choose tsg for bulk model
ts_best(xxx) = MM10(mmm,18); % use MM10 for HiWinGS

%% RH
% psd rh had issues in response time and over saturation on HiWinGS
% possible salt contamination of the fritted metal cap???
for ii=1:144:length(x)  % remove rh spikes and fill in NaNs
    rh(ii:ii+143) = despike2(rh(ii:ii+143));
    rh_im(ii:ii+143) = despike2(rh_im(ii:ii+143));
end
ii = rh>100 & rh<102;  rh(ii)=100;
ii = rh>=102; rh(ii) = NaN;
jj = find(jdy>306.6 & jdy<306.7); rh_im(jj) = NaN;
figure; plot(jdy,rh,'r.',jdy,rh_im,'b.',MM10(:,1),MM10(:,15),'k.');
title(' RH (%), ship(b), PSD(r), MM10(k)');
xlim([stdt endt]); grid;
rh_best = rh_im;    % choose mean ship rh for bulk model
rh_best(xxx) = MM10(mmm,15); % use MM10 for HiWinGS

%% Ta
for ii=1:144:length(x)  % remove ta spikes and fill in NaNs
    ta(ii:ii+143) = despike2(ta(ii:ii+143));
    ta_im(ii:ii+143) = despike2(ta_im(ii:ii+143));
end
figure; plot(jdy,ta,'r.',jdy,ta_im,'b.',MM10(:,1),MM10(:,16),'k.'); grid;
title('Air temperature (degC), ship(b), PSD(r), MM10(k)');
xlim([stdt endt]);
ta_best = ta_im;    % choose avg ship ta for bulk model
ta_best(xxx) = MM10(mmm,16); % use MM10 for HiWinGS

%% Qa
qa = qair_p([ta,rh,press]);
qa_im = qair_p([ta_im,rh_im,P_im]);
figure; plot(jdy,qa,'r.',jdy,qa_im,'b.');
title('Air specific humidity (g/kg), blue=ship, red=PSD');
xlim([stdt endt]); ylim([0 18]); grid;
qa_best = qair_p([ta_best,rh_best,P_best]);

%% Qsea
qse = qsea_p([tsnk press]);
qs_tsg = qsea_p([tsg P_im]);
figure;plot(jdy,qse,'r.',jdy,qs_tsg,'b.');
title(' Specific humidity (g/Kg) ''near'' ocean surface, blue=ship, red=PSD');
xlim([stdt endt]); grid;
qs_best = qsea_p([ts_best,P_best]);

%% Rs
for ii=1:144:length(x)  % subtract night bias in psp each day
    % rs_im(ii+6:ii+18) is 1AM to 3AM each day - night
    rs_im(ii:ii+143) = rs_im(ii:ii+143) - nanmean1(rs_im(ii+6:ii+18));
    rs(ii:ii+143) = rs(ii:ii+143) - nanmean1(rs(ii+6:ii+18));
end
zz = find(rs<-5); rs(zz) = nanmean1(rs(zz-3:zz+3));
zz = find(rs_im<-5); rs_im(zz) = nanmean1(rs_im(zz-3:zz+3));
figure; plot(jdy,rs,'r.',jdy,rs_im,'b.',MM10(:,1),MM10(:,20),'k.');
title('PSP downward solar flux, W/m^2, ship(b), PSD(r), MM10(k)');
xlim([stdt endt]); grid;
rs_best = rs;                      % psd generally better
rs_best(jdy<249) = rs_im(jdy<249); % except prior to jd 249
rs_best(xxx) = MM10(mmm,20); % use MM10 for HiWinGS

%% Rl
ss = find(jdy==251);    % from doy 250 onward...
for ii=ss:144:length(x)  % scale ship to psd value each day
    rl_im(ii:ii+143) = rl_im(ii:ii+143) * nanmean1(rl(ii:ii+143))/nanmean1(rl_im(ii:ii+143));
end
figure; plot(jdy,rl,'r.',jdy,rl_im,'b.',MM10(:,1),MM10(:,21),'k.');
title('Downward IR flux, W/m^2, ship(b), PSD(r), MM10(k)');
xlim([stdt endt]); grid;
rl_best = rl;                      % psd generally better
rl_best(jdy<249) = rl_im(jdy<249); % except prior to jd 249
rl_best(xxx) = MM10(mmm,21); % use MM10 for HiWinGS

%% Rain
% rain difficult to interpret during snowfall
org(org<0) = 0; rain_im(rain_im<0) = 0;     % kill negative values
org(org>30) = 10; rain_im(rain_im>30) = 10; % chop down a few of the largest outliers
figure; plot(jdy,org,'r.',jdy,rain_im,'b.',MM10(:,1),MM10(:,22),'k.');
title('Rainrate, mm/hr, ship(b), PSD(r), MM10(k)');
xlim([stdt endt]); grid;
rain_best = org;    % set to psd rain
rain_best(xxx) = MM10(mmm,22); % use MM10 for HiWinGS

%% SOG/COG/Heading/Lat/Lon plots
figure;plot(jdy,sog,'r.',jdy,sog_im,'b.');title('Ship SOG (b.) & PSD SOG (r.)');
xlim([stdt endt]);ylim([0 7]); grid;
figure;plot(jdy,cog,'r.',jdy,cog_im,'b.');title('Ship COG (b.) & PSD COG (r.)');
xlim([stdt endt]); grid;
figure; plot(Lon,Lat,'r.'); title('Lat/Lon, PSD'); grid;

%% Licor fast humidity & co2
% lots of clean up required here...
Rgas = 287.1;               % Pa m-3 kg-1 K-1 for dry air
Rgas_universal = 8.314472;	% Pa m3 K-1 mol-1
Mw = 18.01528;              % molar mass of H20 g/mol
Md = 28.964;                % molar mass dry air
epsilon = Mw/Md;            % mass of water to mass of dry air = 0.622
P_lic = yxx(:,6);           % kPa
T_lic = yxx(:,7);           % deg C
sgP_lic = yxx(:,16);
AGC = yxx(:,17);

ii = find(AGC>60);  % kill all data where AGC is bad
yxx(ii,2:5) = NaN;
yxx(ii,8:end) = NaN;
cql(ii) = NaN;

co2_lic = (yxx(:,4)/1000) * Rgas_universal.*(T_lic+273.15)./(P_lic*1000) * 1e6; % mmol/m3 to ppm
ii = find(co2_lic<350 | co2_lic>450);
co2_lic(ii) = NaN;  % kill data when CO2 is bad
yxx(ii,2:5) = NaN;
yxx(ii,8:end) = NaN;
cql(ii) = NaN;

h2o_mr = (yxx(:,2)/1000).*Rgas_universal.*(T_lic+273.15)./(P_lic*1000); % mmol/m3 to mol/mol
rhoa_lic = (P_lic*1000).*(1-h2o_mr*(1-epsilon))./(Rgas.*(T_lic+273.15));% air density, kg/m3
q_lic = (yxx(:,2)*1e-3).*Mw./rhoa_lic; % specific humidity from Licor (g/kg)
ii = find(q_lic<0.1 | q_lic>20);
q_lic(ii) = NaN;
co2_lic(ii) = NaN;
yxx(ii,2:5) = NaN;  % kill data when q is bad
yxx(ii,8:end) = NaN;
cql(ii) = NaN;

sgq_lic = (yxx(:,12)*1e-3).*Mw./rhoa_lic;
ii = find(sgq_lic>3 | sgq_lic<0);
sgq_lic(ii) = NaN;  % ...when sigma q is bad
q_lic(ii) = NaN;
co2_lic(ii) = NaN;
yxx(ii,2:5) = NaN;
yxx(ii,8:end) = NaN;
cql(ii) = NaN;

sgco2_lic = (yxx(:,14)/1000) * Rgas_universal.*(T_lic+273.15)./(P_lic*1000) * 1e6;
ii = find(sgco2_lic>20);
sgco2_lic(ii) = NaN; % ...when sigma CO2 is bad
sgq_lic(ii) = NaN;
q_lic(ii) = NaN;
co2_lic(ii) = NaN;
yxx(ii,2:5) = NaN;
yxx(ii,08:end) = NaN;
cql(ii) = NaN;

% discard when number of missing data points exceeds limit
ii = find(missing_lic > missLic_lim);
sgco2_lic(ii) = NaN;
sgq_lic(ii) = NaN;
q_lic(ii) = NaN;
co2_lic(ii) = NaN;
yxx(ii,2:5) = NaN;
yxx(ii,08:end) = NaN;
cql(ii) = NaN;

wq_lic = yxx(:,8);
ii = wq_lic<-.02 | wq_lic>0.1;
wq_lic(ii) = NaN; % discard bad covariances

 % q licor plot
figure;plot(jdy,qa_best,'b.',jdy,q_lic,'r.');
title('q (g/kg), Vaisala(b.), Licor(r.)');
xlim([stdt endt]); ylim([0 15]); grid;

% rq = qa./q_lic;
% fq = 1.05*max(1,1+.89*(rq-1.044));
% ii = find(isnan(q_lic));fq(ii)=NaN;
% wq_lic = wq_lic.*fq;  % correct for separation

%% check covariances for reasonable limits
ii = find(wuj<-5 | wuj>.05);         wuj(ii) = NaN;
ii = find(wu_uh<-5 | wu_uh>.05);     wu_uh(ii) = NaN;
ii = find(wvj<-1 | wvj>1);           wvj(ii) = NaN;
ii = find(wv_uh<-1 | wv_uh>1);       wv_uh(ii) = NaN;
ii = find(wt_son<-.1 | wt_son>.25);  wt_son(ii) = NaN;
ii = find(wt_son_uh<-.1 | wt_son_uh>.25);  wt_son_uh(ii) = NaN;

% ii = find(tt>0.2); tt(ii) = NaN; % tt
% ii = find(uuj>6); uuj(ii) = NaN; % uu edson
% ii = find(vvj>6); vvj(ii) = NaN; % vv edson
% ii = find(wwj>1); wwj(ii) = NaN; % ww edson

% cu and cw
%ii = find(cu<0 | cu>1.3); cu(ii) = NaN;%  cup(ii) = NaN;
%kk = find(cw<0 | cw>1);   cw(kk) = NaN;%  cwp(kk) = NaN;
%ii = find(cu_uh<0 | cu_uh>1.3); cu_uh(ii) = NaN;%  cup_uh(ii) = NaN;
%kk = find(cw_uh<0 | cw_uh>1);   cw_uh(kk) = NaN;%  cwp_uh(kk) = NaN;

% ct & cql
ii = find(cql<1e-4 | cql>2);   cql(ii) = NaN;
ii = find(ct<1e-6 | ct>0.6);   ct(ii) = NaN;
ii = find(ct_uh<1e-6 | ct_uh>0.6);   ct_uh(ii) = NaN;

%% compute other 'best' values
sp2_best = sp2_uh;
sp2_best(isnan(sp2_best)) = sp2(isnan(sp2_best));
sp2_best(sp2_best>4) = NaN;
vv_best = vv_uh;
vv_best(isnan(vv_best)) = vvj(isnan(vv_best));
vv_best(vv_best>10) = NaN;
ww_best = ww_uh;
ww_best(isnan(ww_best)) = wwj(isnan(ww_best));
ww_best(ww_best>3) = NaN;

%% set measurement heights for 'best' met measurements
zu = 15.91;   % UH sonic, wind speed measurement height (m)
zt = 15;      % SCS air T measurement height (m)
zq = 15;      % SCS air q measurement height (m)

%% save cleaned / adjusted values back to x & yxx arrays
% substitute PSD data with 'best' values from above - used for hr ave & hr bulk
% can comment this out to preserve original x and yxx arrays in memory
x(:,2) = U_scs;
x(:,3) = ta_im;
x(:,4) = rh_im;
x(:,48) = rh_im;
x(:,5) = P_im;
x(:,6) = tsg;
x(:,7) = rs_im;
x(:,8) = rl_im;
x(:,10) = rain_im;
x(:,16) = u_best; % substitute psd Utrue with u_best
x(:,17) = ta_best;
x(:,18) = rh_best;
x(:,19) = P_best;
x(:,20) = ts_best;
x(:,21) = rs_best;
x(:,22) = rl_best;
x(:,24) = rain_best;
x(:,41) = dir_best; % substitute psd dir with u_best
x(:,42) = reldir_best;
x(:,49) = qa_best;
x(:,51) = qs_best;

x(:,117) = wuj;
x(:,119) = wvj;
x(:,121) = wt_son;
x(:,125) = uuj;
x(:,126) = vv_best;
x(:,127) = ww_best;
x(:,128) = tt;
x(:,130) = sp2_best;
x(:,185) = wu_uh;
x(:,186) = wv_uh;
x(:,187) = wt_son_uh;
x(:,214) = tilt_best;

switch IDflag
    case 'a'
        x(:,134) = cu;        % method a
        x(:,136) = cw;
        x(:,138) = ct;
        x(:,140) = cql;
        x(:,195) = cu_uh;
        x(:,197) = cw_uh;
        x(:,199) = ct_uh;
    case 'b'
        x(:,135) = cu;        % method b
        x(:,137) = cw;
        x(:,139) = ct;
        x(:,141) = cql;
        x(:,196) = cu_uh;
        x(:,198) = cw_uh;
        x(:,200) = ct_uh;
end

yxx(:,2) = q_lic;
yxx(:,4) = co2_lic;
yxx(:,8) = wq_lic;
yxx(:,12) = sgq_lic;
yxx(:,14) = sgco2_lic;
