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

iii = find(MMF.yday>=275 & MMF.yday<279);

%% WS
figure; plot(jdy,U,'r.',jdy,U_scs,'b.');
xlim([stjd jdy(end)]); ylim([0 20]);grid;
title(' True wind speed - relative to earth(m/s), blue=ship, red=PSD');

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

%% Pressure
[press,~] = despike2(press); [P_im,~] = despike2(P_im);
figure; plot(jdy,press,'r.',jdy,P_im,'b.'); % imet is sfc p
title('Pressure (mb), PSD(r) and ship(b)');
xlim([stjd jdy(end)]); 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
tsnk = MMF.Tsea(iii)';
figure;plot(jdy,tsnk,'r.-',jdy,tsg,'b.'); grid;
title(' SST(degC), blue=ship tsg, red=PSD seasnake');
xlim([stjd jdy(end)]); 

%% RH
rh =((x(:,47)-3)+0.41)/0.95; % fix psd rh
figure; plot(jdy,rh,'r.',jdy,rh_im,'b.');
title(' RH (%), blue=ship, red=PSD');
xlim([stjd jdy(end)]); grid;

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

%% Qa
qa=qair_p([ta rh press]);
figure; plot(jdy,qa,'r.',jdy,qa_im,'b.');
title('Air specific humidity (g/kg), blue=ship, red=PSD');
xlim([stjd jdy(end)]); ylim([14 23]); 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=PSD');
xlim([stjd jdy(end)]); grid;

%% Rs
zz = find(rs<-5); rs(zz) = NaN; rs(zz) = nanmean1(rs(zz-3:zz+3));
figure; plot(jdy,rs,'r.',jdy,rs_im,'b.');
title('PSP downward solar flux, W/m^2, blue=ship, red=PSD');
xlim([stjd jdy(end)]); grid;

%% Rl
figure; plot(jdy,rl,'r.',jdy,rl_im,'b.');
title('Downward IR flux, W/m^2, PSD(red.), SCS(blue.)');
xlim([stjd jdy(end)]); grid;

%% Rain
figure; plot(jdy,org,'r.');
title('PSD Rainrate, mm/hr');
xlim([stjd jdy(end)]); grid;

%% SOG/COG/Heading/Lat/Lon plots

%finds weird stuff with COG SCS and replace with PSD GPS
cog_im(jdy>277 & jdy<301.8) = cog(jdy>277 & jdy<301.8);

figure;plot(jdy,sog,'r.',jdy,sog_im,'b.');title('Ship SOG (b.) & PSD SOG (r.)');
xlim([stjd jdy(end)]);ylim([0 6]); grid;
figure;plot(jdy,cog,'r.',jdy,cog_im,'b.');title('Ship COG (b.) & PSD COG (r.)');
xlim([stjd jdy(end)]);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);
wq_lic   =   yxx(:,8);

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<340 | co2_lic>375);
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<10 | q_lic>15);
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>2 | sgq_lic<5e-3);
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>15);
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;

 % q licor plot
figure;plot(jdy,qa,'b.',jdy,q_lic,'r.'); grid;
title('q (g/kg), Vaisala(b.), Licor(r.)');
xlim([stjd jdy(end)]); ylim([8 23]);

%% check covariances for reasonable limits
ii = wuj<-5 | wuj>.05;         wuj(ii) = NaN;
ii = wvj<-0.2 | wvj>0.2;       wvj(ii) = NaN;
ii = wq_lic<-.1 | wq_lic>0.2;  wq_lic(ii) = NaN;

ii = wt_son<-.01 | wt_son>.1;  wt_son(ii) = NaN;
toto = find(jdy>274 & jdy<280);
ii = wt_son(toto)<-.1 | wt_son(toto)>0.05; wt_son(toto(ii))=NaN;  %removes weird points at beginning and end

% 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 = cql<1e-4 | cql>2;   cql(ii) = NaN;
ii = 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;

x(:,134) = cu;
x(:,136) = cw;
x(:,138) = ct;
x(:,140) = cql;
% x(:,135) = cu;
% x(:,137) = cw;
% x(:,139) = ct;
% x(:,141) = cql;

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

