% ttide_adcp_indiv.m  An mfile to run T_TIDE on South Carolina ADCP 
%                     velocity and pressure data.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Use of this program is described in:
%
% Armstrong, B.N., Warner, J.C., Voulgaris, G., List, J.H., Thieler, E.R., 
% Martini, M.M., Montgomery, E.T., 2011, Carolinas Coastal Change Processes 
% Project Data Report for Observations near Diamond Shoals, NC from January-April 
% 2009: U.S. Geological Survey Open-File Report 2011-1156.
%
% and
%
% Sullivan, C.M., Warner, J.C., Martini, M.A., Voulgaris, G., 
% Work, P.A., Haas, K.A., and Hanes, D.H. (2006) 
% South Carolina Coastal Erosion Study Data Report for Observations
% October 2003 - April 2004., USGS Open-File Report 2005-1429.
%
% Program written in Matlab v7.1.0 SP3
% Program ran on PC with Windows XP Professional OS.
%
% "Although this program has been used by the USGS, no warranty, 
% expressed or implied, is made by the USGS or the United States 
% Government as to the accuracy and functioning of the program 
% and related program material nor shall the fact of distribution 
% constitute any such warranty, and no responsibility is assumed 
% by the USGS in connection therewith."
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

% C. Sullivan 05/15/06, Water flow data in bins likely exposed to air has
%                       been set to the NetCDF FillValue_.  Editing code to
%                       consider the surface bin the upward-most good bin.
% C. Sullivan 01/11/05, Use data from surface, bottom, and user-selected
%                       mid-depths.  Run the analysis on each 3-month
%                       timeseries using that timeseries' start time and
%                       time step.
addpath('D:\BNA\Matlab\m_cmg\branches\barmstrong\BNA\t_tide_v1.3beta','-end');
more off

% Define data and output directories
dataDirectory = '\\IGSAGIEGGS-FILE\Stg-proc\data\HatterasNearshore10\BBV_reviewed';
plotDirectory = 'D:\BNA\data\Carolinas\webpage\dataDVD\Hatteras2010DVD\tides';

% Define adcp data filenames.
data1={'8671Aaqd-cal.nc','8671Baqd-cal.nc'};
% data2={'8731whAbeam.nc','8731whBbeam.nc','8731whCbeam.nc','8731whDbeam.nc','8731whEbeam_trm.nc'};
% data3={'8741whAbeam_trm.nc','8741whBbeam.nc','8741whCbeam.nc','8741whDbeam.nc','8741whEbeam_trm.nc'};
dataFiles={'8671'};%,'8731','8741'};



for f = 1:length(dataFiles)
    P_4_merge=[];
    uvbot_merge=[];
    uvsfc_merge=[];
    jd=[]
    eval(['no_files=length(data',num2str(f),')']);
    for no_dat=1:no_files
        
        eval(['nc=(fullfile(dataDirectory,data',num2str(f),'{no_dat}));']);
        
        %get the site, the latitude, and the depth of
        %each velocity bin
        site=dataFiles{f}(1:4);
        
        latitude=ncread(nc,'lat');
        binDepths=squeeze(ncread(nc,'depth'));
        
        %get time in julian days
        if no_dat==1
            jd=[jd double(ncread(nc,'time')) + (double(ncread(nc,'time2'))/3600/1000/24)];
        else
            jd=[jd;double(ncread(nc,'time')) + (double(ncread(nc,'time2'))/3600/1000/24)];
        end
        %get bottom pressure and convert it to height
        %with H=P/(rho*g) where rho=1025 km/m^3 and
        %g=9.8 m/s^2
        try
            P_4=squeeze(ncread(nc,'P_1294'));
        catch
            P_4=squeeze(ncread(nc,'P_1'));
        end
        P_4(P_4>=1e35)=nan; %nan bad points
        if no_dat==1
        P_4_merge=[P_4_merge P_4];
        else
        P_4_merge=[P_4_merge;P_4]; 
        end
        
        ubot = ncread(nc,'u_1205',[1 1 1 1],[1 1 1 Inf]);
        vbot = ncread(nc,'v_1206',[1 1 1 1],[1 1 1 Inf]);
        bads=find(ubot+vbot >= 1e35);
        ubot(bads) = nan;
        vbot(bads) = nan;
        
        uvbot = ubot + i*vbot;
        if no_dat==1
        uvbot_merge=[uvbot_merge squeeze(uvbot)];
        else
        uvbot_merge=[uvbot_merge;squeeze(uvbot)];            
        end
        
        dims=(size(ncread(nc,'u_1205',[1 1 1 1],[1 1 Inf 1])));
        lastGoodBin = dims(3);
        depthLastGoodBin = binDepths(lastGoodBin)
        usfc = ncread(nc,'u_1205',[1 1 lastGoodBin 1],[1 1 1 Inf]);
        vsfc = ncread(nc,'v_1206',[1 1 lastGoodBin 1],[1 1 1 Inf]);
        bads=find(usfc+vsfc >= 1e35);
        usfc(bads) = nan;
        vsfc(bads) = nan;
        
        uvsfc = usfc + i*vsfc;
        if no_dat==1
        uvsfc_merge=[uvsfc_merge squeeze(uvsfc)];
        else
        uvsfc_merge=[uvsfc_merge;squeeze(uvsfc)];
        end
    end
    Height=P_4_merge/(1025*9.8);
    display(['The length of the P_4 data set is ' num2str(length(P_4_merge))])
    display(['The length of the uvsfc data set is ' num2str(length(uvsfc_merge))])
    display(['The length of the uvbot data set is ' num2str(length(uvbot_merge))])
    
    %%%%%%run t_tide on pressure data (bottom pressure)
    if ~isempty(Height)
        startInd = find(~isnan(Height),1,'first');
        startTime = julian2datenum(jd(startInd));
        interval = mean(diff(julian2datenum(jd)))*24
        ttideOutFile = fullfile(plotDirectory,['ttide_pres_',dataFiles{f}(1:4),'.txt']);
        [nameu,fu,tidecon,xout]=t_tide(P_4_merge,'start time',startTime,...
            'interval',interval,...
            'latitude',latitude,...
            'output',ttideOutFile);
    end
    
    %%%%%%run t_tide on velocity data at the bottom
    startInd = find(~isnan(uvbot),1,'first');
    startTime = julian2datenum(jd(startInd));
    interval = mean(diff(julian2datenum(jd)))*24;
    ttideOutFile = fullfile(plotDirectory,['ttide_curr_',dataFiles{f}(1:4),'_bottom.txt']);
    [nameu,fu,tidecon,xout]=t_tide(uvbot_merge,'start time',startTime,...
        'interval',interval,...
        'latitude',latitude,...
        'output',ttideOutFile);
    
    %%%%%%run t_tide on velocity data at the surface
    %%%%%%C. Sullivan, 05/15/06: the surface is the upward-most good bin
    startInd = find(~isnan(uvsfc),1,'first');
    startTime = julian2datenum(jd(startInd));
    interval = mean(diff(julian2datenum(jd)))*24;
    ttideOutFile = fullfile(plotDirectory,['ttide_curr_',dataFiles{f}(1:4),'_surface.txt']);
    [nameu,fu,tidecon,xout]=t_tide(uvsfc_merge,'start time',startTime,...
                                         'interval',interval,...
                                         'latitude',latitude,...
                                         'output',ttideOutFile);
   
end
   
    
    
    
    