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

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

% Adjustments to variables used for COARE input will be used in the model run
% but mean met values in the 10-min and 1-hr flux files will be the unadjusted
% original unless the adjustments are also copied back into the x array.

% 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('CAPRICORN_2016_fixit_da1')

%% WS
figure; plot(jdy,U,'r.',jdy,U_scs,'b.');
xlim([stdt endt]); grid;
title(' True wind speed - relative to earth(m/s), blue=ship, red=NOAA');

%% WD
figure; plot(jdy,dir,'r.',jdy,dir_scs,'b.');
xlim([stdt endt]); grid;
title(' True wind direction - relative to earth (deg), blue=ship, red=NOAA');

%% Pressure
figure; plot(jdy,press,'r.',jdy,P_im,'b.'); % imet is sfc p
title('Pressure (mb), NOAA(r) and ship(b)');
xlim([stdt endt]); grid;

%% Tilt angle - flow distortion angle derived from the motion correction script
% this is tilt angle prior to ship speed correction
figure; plot(reldir,tiltx,'b.',[-90 -90],[-5 25],'r',[90 90],[-5 25],'r');
axis([-180 180 -5 25]); grid; title('Tilt Angle, degrees');
clear reldir2

%% Tsea
% eliminate periods when snake is out of water...
bbb = find(jdy<75 | (jdy>75.8&jdy<76.6) | (jdy>78.83&jdy<79.26) | ...
    (jdy>79.825&jdy<80.28) | (jdy>82.040&jdy<82.93) | (jdy>85.84&jdy<88.345) | ...
    (jdy>88.87&jdy<90.06) | (jdy>96.4&jdy<97.2) | (jdy>101.86&jdy<102.12) | ...
    (jdy>104.888));
tsnk(bbb) = NaN;
figure;plot(jdy,tsg,'b.',jdy,tsnk,'r-'); grid;
title(' SST(degC), blue=ship tsg, red=NOAA seasnake');
xlim([stdt endt]);

%% RH
rh(rh>100) = 100; rh_im(rh_im>100)  = 100; % fix for foggy conditions
figure; plot(jdy,rh,'r.',jdy,rh_im,'b.');
title(' RH (%), blue=ship, red=NOAA');
xlim([stdt endt]); ylim([-2,102]);grid;

%% Ta
figure; plot(jdy,ta,'r.',jdy,ta_im,'b.');
title('Air temperature (degC), blue=ship, red=NOAA');
xlim([stdt endt]); grid;

%% Qa
figure; plot(jdy,qa,'r.',jdy,qa_im,'b.');
title('Air specific humidity (g/kg), blue=ship, red=NOAA');
xlim([stdt endt]); grid;

%% 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=NOAA');
xlim([stdt endt]); grid;

%% Rs
zz = find(rs<-5); rs(zz) = NaN; rs(zz) = nanmean1(rs(zz-3:zz+3));
zz = find(jdy>83.34 & jdy<83.85);
rs_im(zz) = rs(zz); rsmax_im(zz) = rs(zz);
figure; plot(jdy,rs,'r.',jdy,rs_im,'b.'); %,jdy,rsmax_im,'k--');
title('PSP downward solar flux, W/m^2, blue=ship, red=NOAA');
xlim([stdt endt]); grid;

%% Rl
zz = find(jdy>83.32 & jdy<84.1);
rl_im(zz) = NaN;
figure; plot(jdy,rl,'r.',jdy,rl_im,'b.');
title('Downward IR flux, W/m^2, NOAA(red.), Ship(blue.)');
xlim([stdt endt]); grid;

%% Rain
figure; plot(jdy,org,'r.',jdy,rain_im,'b.');
title('Rainrate, mm/hr'); legend('NOAA Org','Ship RM Young','location','northeast');
xlim([stdt endt]); grid;

%% SOG/COG/Heading/Lat/Lon plots
figure;plot(jdy,sog,'r.',jdy,sog_im,'b.');title('Ship SOG (b.) & NOAA SOG (r.)');
xlim([stdt endt]);ylim([0 6]); grid;
figure;plot(jdy,cog,'r.',jdy,cog_im,'b.');title('Ship COG (b.) & NOAA COG (r.)');
xlim([stdt endt]);grid;
figure; plot(Lon,Lat,'r.'); title('Lat/Lon, NOAA'); 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
co2_lic = (yxx(:,4)/1000) * Rgas_universal.*(ta+273.15)./(press*100) * 1e6; % using psd ta and press
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
h2o_mr = (yxx(:,2)/1000).*Rgas_universal.*(ta+273.15)./(press*100); % using psd ta and press
rhoa_lic = (press*100).*(1-h2o_mr*(1-epsilon))./(Rgas.*(ta+273.15));
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;
sgco2_lic = (yxx(:,14)/1000) * Rgas_universal.*(ta+273.15)./(press*100) * 1e6; % using psd ta and press
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;

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,'b.',jdy,q_lic,'r.'); grid;
title('q (g/kg), Vaisala(b.), Licor(r.)');
xlim([stdt endt]);

%% check covariances for reasonable limits
ii = wuj<-5 | wuj>.05;         wuj(ii) = NaN;
ii = wvj<-1 | wvj>1;           wvj(ii) = NaN;
ii = wt_son<-.1 | wt_son>.25;  wt_son(ii) = NaN;
ii = wq_lic<-.02 | wq_lic>0.1; wq_lic(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
cup = cu;
cwp = cw;
ii = find(cup<0 |cup >1.3); cup(ii) = NaN;  cu(ii) = NaN;
kk = find(cwp<0 |cwp >1);   cwp(kk) = NaN;  cw(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;

% Ts noise level
ii = find(Ts_noise > 0.001);    ct(ii) = NaN;  wt_son(ii) = NaN;

% 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

%% save changes to arrays in memory
% can comment this out to preserve original x and yxx arrays in memory
% x(:,2) = U_scs;
% x(:,3) = ta_im;
% x(:,5) = P_im;
% x(:,6) = tsg;
x(:,7) = rs_im;
% x(:,8) = rl_im;
% x(:,16) = U;
% x(:,17) = ta;
% x(:,18) = rh;
% x(:,19) = press;
x(:,20) = tsnk;
x(:,21) = rs;
% x(:,22) = rl;
% x(:,24) = org;
% x(:,41) = dir;
% x(:,44) = tiltx;


x(:,117) = wuj;
x(:,119) = wvj;
x(:,121) = wt_son;

x(:,128) = tt;
x(:,125) = uuj;
x(:,126) = vvj;
x(:,127) = wwj;

switch IDflag
    case 'a'
        x(:,134) = cu;      % ID METHOD A
        x(:,136) = cw;
        x(:,138) = ct;
        x(:,140) = cql;
    case 'b'
        x(:,135) = cu;      % ID METHOD B
        x(:,137) = cw;
        x(:,139) = ct;
        x(:,141) = cql;
end;

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

