% comp_waterlevel
%       reads in model and observed WATER LEVEL data, and returns a comparison
%       result.
%
%--------
% input
%--------
%   Observed data:
%       72**adc-a.nc    - sites 1-8 for a period of fall 2003.
%       74**adc-a.nc    - sites 1-8 for a period of spring 2004.
%
%   Model data:
%       sc_wrk9r1_his_***.nc
%
%--------
% output
%--------
%       *.mat
%
%--------
% history
%--------
%   10/2/06 Edited to run a comparison of velocities from between
%   simulations and observations. -hsk
%
%--------------------------------------------------------------------------

%--- Measured Data
%               for fall 2003
%    data filename                  site name
%========================    ======================
ddir='G:\CSULLIVAN\SOUTH_CAROLINA\DATA_DVD\DATAFILES\';

oncs{1}='7223pt-a.nc';          site{1}='site 2';  
oncs{2}='7243pt-a.nc';          site{2}='site 3';
oncs{3}='7331adc-a.nc';         site{3}='site 4';
oncs{4}='7261advs-a.nc';        site{4}='site 6a';
oncs{5}='7262advs-a.nc';        site{5}='site 6b';
oncs{6}='7281advs-a.nc';        site{6}='site 7A1';
oncs{7}='728pcvp-cal5.nc';      site{7}='site 7A2';
oncs{8}='7321adc-a.nc';         site{7}='site 7B';
oncs{9}='7311aqc-a.nc';         site{7}='site 8';

% P_1(time,depth,lat,lon) in dbar

xys=[
[-78.84902954	33.64871597];
[-78.68206787	33.77388382];
[-78.78520203	33.70080185];
[-78.76028333	33.70684814];
[-78.76028333	33.70684814];
[-78.77911667	33.72246552];
[-78.77911667	33.72246552];
[-78.78060150	33.72230148];
[-78.77200317	33.71699905]];

[xs,ys]=meshgrid(xys(:,1),xys(:,2));

% find indicies for a diagonal
indx4diag=[1:size(xs,1)*size(xs,2)];
indx4diag=reshape(indx4diag,size(xs,1),size(xs,2));
indx4diag=diag(indx4diag);

%--- Model Data
mdir='E:\HKIM\SCarolina\RomsP\large-domain\output\';
prefix='sc_wrk9r2_his_';
suffix='.nc';

%============================= MAIN =======================================

%--- Load Model Data
Gfile='E:\HKIM\SCarolina\Romsp\large-domain\grid\sc_grid_26.nc';
gnc=netcdf(Gfile,'r');
lon=gnc{'lon_rho'}(:); lat=gnc{'lat_rho'}(:); close(gnc);

% # of records each file
nrcd=72;
lon3=repmat(lon,[nrcd,1,1]); lon3=reshape(lon3,[nrcd,size(lon,1),size(lon,2)]);
lat3=repmat(lat,[nrcd,1,1]); lat3=reshape(lat3,[nrcd,size(lat,1),size(lat,2)]);

xs3=repmat(xs,[nrcd,1,1]); xs3=reshape(xs3,[nrcd,size(xs,1),size(xs,2)]);
ys3=repmat(ys,[nrcd,1,1]); ys3=reshape(ys3,[nrcd,size(ys,1),size(ys,2)]);

Zeta=[];

for N=1:23
    fname=[prefix,num2str(N,'%04d'),suffix];
    hnc=netcdf(fname,'r');
    ztmp=hnc{'zeta'}(:); close(hnc);  
    
    zi=griddata(lon3,lat3,ztmp,xs3,ys3);
    
    Zeta=[Zeta; squeeze(zi(:,indx2diag))];
end 

dbar=conv_prs_m2dbar(Zeta);

%--- Load Observed Data
nc1=netcdf([ddir,char(oncs{1})],'r');
jd2=nc1{'time'}(:)+nc1{'time2'}(:)/3600/24/1000;
P2=nc1{'P_1'}(:);
close(nc1); clear nc1

nc1=netcdf([ddir,char(oncs{2})],'r');
jd3=nc1{'time'}(:)+nc1{'time2'}(:)/3600/24/1000;
P3=nc1{'P_1'}(:);
close(nc1); clear nc1

nc1=netcdf([ddir,char(oncs{2})],'r');
jd4=nc1{'time'}(:)+nc1{'time2'}(:)/3600/24/1000;
P4=nc1{'P_1'}(:);
close(nc1); clear nc1

nc1=netcdf([ddir,char(oncs{2})],'r');
jd6a=nc1{'time'}(:)+nc1{'time2'}(:)/3600/24/1000;
P6a=nc1{'P_1'}(:);
close(nc1); clear nc1

nc1=netcdf([ddir,char(oncs{2})],'r');
jd6b=nc1{'time'}(:)+nc1{'time2'}(:)/3600/24/1000;
P6b=nc1{'P_1'}(:);
close(nc1); clear nc1

nc1=netcdf([ddir,char(oncs{2})],'r');
jd7A1=nc1{'time'}(:)+nc1{'time2'}(:)/3600/24/1000;
P7A1=nc1{'P_1'}(:);
close(nc1); clear nc1

nc1=netcdf([ddir,char(oncs{2})],'r');
jd7A2=nc1{'time'}(:)+nc1{'time2'}(:)/3600/24/1000;
P7A2=nc1{'P_1'}(:);
close(nc1); clear nc1

nc1=netcdf([ddir,char(oncs{2})],'r');
jd7B=nc1{'time'}(:)+nc1{'time2'}(:)/3600/24/1000;
P7B=nc1{'P_1'}(:);
close(nc1); clear nc1

nc1=netcdf([ddir,char(oncs{2})],'r');
jd8=nc1{'time'}(:)+nc1{'time2'}(:)/3600/24/1000;
P8=nc1{'P_1'}(:);
close(nc1); clear nc1