%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Set path to include my .m files
%
path(path,'../matlab_cen');
%
% Load observational data
%
OPATH='../OBS/'
UVOF=[OPATH,'UV_INTEG.m3d']
TAUOF=[OPATH,'TAUXY.o1d']
QOF=[OPATH,'QFLUX.o1d']
%
[UVO,gridname,year,ncol,gmt]=read_obs(UVOF);
[TAUO,gridname,year,ncol,gmt]=read_obs(TAUOF);
[QO,gridname,year,ncol,gmt]=read_obs(QOF);
%
% Load model results
%
MPATH='../OBS/'
UVMF=[MPATH,'UV_INTEG-Q4.m3d']
TAUMF=[MPATH,'TAUXY.Q4.m1d']
QMF=[MPATH,'QFLUX.Q4.m1d']
%
[UVEM,gridname,year,ncol,gmt]=read_obs_trunc(UVMF);
[TAUM,gridname,year,ncol,gmt]=read_obs_trunc(TAUMF);
[QM,gridname,year,ncol,gmt]=read_obs_trunc(QMF);
%
% Remove rows of UVO which do not have UVEM counterparts
%
nmod=length(UVEM(:,1));
for i=1:length(UVEM(:,1))
   ii(i)=find(abs(UVO(:,1)-UVEM(i,1))<0.0001);
end
UVO=UVO(ii,:);
%
% Plot time-series of input data
%
tallfigure;
%
subplot(3,1,1);
hold on
hp=plot(UVO(:,1),UVO(:,4),'b--');
hp=plot(UVO(:,1),UVO(:,5),'r-');
UM=UVO(:,4)-UVEM(:,4);
VM=UVO(:,5)-UVEM(:,5);
hp=plot(UVEM(:,1),UM,'b-');
set(hp,'LineWidth',2.0)
hp=plot(UVEM(:,1),VM,'r--');
set(hp,'LineWidth',2.0)
legend('Obs U','Obs V','Mod U','Mod V');
ylabel('Velocity (m/s)');
title(['rms error=',num2str(rms(UVEM(:,4)+sqrt(-1)*UVEM(:,5))),' m/s'])
%
subplot(3,1,2);
hold on
hp=plot(TAUO(:,1),TAUO(:,2),'bo');
set(hp,'MarkerSize',4.0)
hp=plot(TAUO(:,1),TAUO(:,3),'rx');
set(hp,'MarkerSize',4.0)
hp=plot(TAUM(:,1),TAUM(:,2),'b-');
set(hp,'LineWidth',2.0)
hp=plot(TAUM(:,1),TAUM(:,3),'r--');
set(hp,'LineWidth',2.0)
legend('Obs \tau_x','Obs \tau_y','Mod \tau_x','Mod \tau_y');
ylabel('Wind Stress (Pa)');
%
subplot(3,1,3);
hold on
hp=plot(QO(:,1),QO(:,2),'bo');
set(hp,'MarkerSize',4.0)
hp=plot(QM(:,1),QM(:,2),'b-');
set(hp,'LineWidth',2.0)
legend('Obs q','Mod q');
ylabel('Heat Flux (w/m^2)');
xlabel('Time (ydgmt0)');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
