% ttide_aqdp_comb.m An mfile to run T_TIDE on South Carolina Aquadopp % velocity and pressure data. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Use of this program is described in: % % 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/31/05, This version use velocity data from surface, % bottom, and user-selected mid-depths. The pressure % used is the bottom pressure (P_4). The harmonic % analysis is run using the entire timeseries % thus combining the data from deployments 1 and 2. more off % Define data and output directories dataDirectory = 'E:\CSULLIVAN\SOUTH_CAROLINA\DATA_DVD\DATAFILES'; plotDirectory = 'E:\CSULLIVAN\SOUTH_CAROLINA\PLOTS\TTIDE'; % Define ADCP data filenames. dataFiles={'7311aqc-a.nc','7511aqc-a.nc'}; % The following loop loads the entire timeseries of % pressure and velocities at site 2 and runs harmonic % analysis on the data. for f = 2:2:length(dataFiles) nc1=netcdf(fullfile(dataDirectory,dataFiles{f-1}), 'nowrite'); nc2=netcdf(fullfile(dataDirectory,dataFiles{f}), 'nowrite'); site=nc1.DESCRIPT(6); latitude=nc1.latitude(:); binDepths1=nc1{'depth'}(:); binDepths2=nc2{'depth'}(:); %time in julian days jd1=nc1{'time'}(:) + (nc1{'time2'}(:)/3600/1000/24); jd2=nc2{'time'}(:) + (nc2{'time2'}(:)/3600/1000/24); rawTime=[jd1;jd2]; %time interval assuming dt during deployment 1 %is the same as dt during deployment 2 interval1=mean(diff(jd1)); interval2=mean(diff(jd2)); interval=interval1; %will interpolate data onto an array that is time %in julian days. interpTime=[rawTime(1):interval:rawTime(end)]; %%%%%%run t_tide on pressure data (bottom pressure) %bottom pressure converted to height using H=P/(rho*g) %where rho=1025 km/m^3 and g=9.8 m/s^2 P_41=nc1{'P_4'}(:); P_41(P_41>=1e35)=nan; Height1=P_41/(1025*9.8); P_42=nc2{'P_4'}(:); P_42(P_42>=1e35)=nan; Height2=P_42/(1025*9.8); %combine deployments rawHeight=[Height1;Height2]; %interpolate data [interpHeight]=smart_interp(rawTime,rawHeight,interpTime,3); startInd = find(~isnan(interpHeight),1,'first'); startTime = julian2datenum(interpTime(startInd)); interval = mean(diff(julian2datenum(interpTime)))*24; ttideOutFile = fullfile(plotDirectory,['ttide_pres_site',site,'.txt']); [nameu,fu,tidecon,xout]=t_tide(interpHeight,'start time',startTime,... 'interval',interval,... 'latitude',latitude,... 'output',ttideOutFile); %%%%%%run t_tide on velocity data at the bottom ubot1 = [nc1{'u_1205'}(:,1)]; vbot1 = [nc1{'v_1206'}(:,1)]; bads=find(ubot1+vbot1 >= 1e35); ubot1(bads) = nan; vbot1(bads) = nan; ubot2 = [nc2{'u_1205'}(:,1)]; vbot2 = [nc2{'v_1206'}(:,1)]; bads=find(ubot2+vbot2 >= 1e35); ubot2(bads) = nan; vbot2(bads) = nan; raw_uvbot=[ubot1 + i*vbot1; ubot2 + i*vbot2]; %interpolate data [interp_uvbot]=smart_interp(rawTime,raw_uvbot,interpTime,3); startInd = find(~isnan(interp_uvbot),1,'first'); startTime = julian2datenum(interpTime(startInd)); interval = mean(diff(julian2datenum(interpTime)))*24; ttideOutFile = fullfile(plotDirectory,['ttide_curr_site',site,'_bottom.txt']); [nameu,fu,tidecon,xout]=t_tide(interp_uvbot,'start time',startTime,... 'interval',interval,... 'latitude',latitude,... 'output',ttideOutFile); %%%%%%run t_tide on velocity data at the surface %%%%%%C. Sullivan, 05/15/06: don't ttide the combined at the surface %%%%%%b/c the surface (the upward-most good bins) for deployments 1 %%%%%%and 2 are > 0.5 m apart % usfc1 = [nc1{'u_1205'}(:,end)]; % vsfc1 = [nc1{'v_1206'}(:,end)]; % bads=find(usfc1+vsfc1 >= 1e35); % usfc1(bads) = nan; % vsfc1(bads) = nan; % % %addtional bad data at beginning of 1st deply. timeseries % usfc1=usfc1(200:end); % vsfc1=vsfc1(200:end); % % usfc2 = [nc2{'u_1205'}(:,1)]; % vsfc2 = [nc2{'v_1206'}(:,1)]; % bads=find(usfc2+vsfc2 >= 1e35); % usfc2(bads) = nan; % vsfc2(bads) = nan; % % raw_uvsfc=[usfc1 + i*vsfc1; usfc2 + i*vsfc2]; % % %interpolate data % [interp_uvsfc]=smart_interp(rawTime,raw_uvsfc,interpTime,3); % % startInd = find(~isnan(interp_uvsfc),1,'first'); % startTime = julian2datenum(interpTime(startInd)); % interval = mean(diff(julian2datenum(interpTime)))*24; % ttideOutFile = fullfile(plotDirectory,['ttide_curr_site',site,'_surface.txt']); % [nameu,fu,tidecon,xout]=t_tide(interp_uvsfc,'start time',startTime,... % 'interval',interval,... % 'latitude',latitude,... % 'output',ttideOutFile); %%%%%%run t_tide on velocity data at 4 mbs and 6 mbs targetDepths=[4,6]; %mbs for b=1:size(targetDepths,2) [index_num,dist_to_hit] = value2Index(binDepths1,targetDepths(b),0.5); u1 = [nc1{'u_1205'}(:,index_num)]; v1 = [nc1{'v_1206'}(:,index_num)]; bads=find(u1+v1 >= 1e35); u1(bads) = nan; v1(bads) = nan; [index_num,dist_to_hit] = value2Index(binDepths2,targetDepths(b),0.5); u2 = [nc2{'u_1205'}(:,index_num)]; v2 = [nc2{'v_1206'}(:,index_num)]; bads=find(u2+v2 >= 1e35); u2(bads) = nan; v2(bads) = nan; raw_uv=[u1 + i*v1; u2 + i*v2]; %interpolate data [interp_uv]=smart_interp(rawTime,raw_uv,interpTime,3); startInd = find(~isnan(interp_uv),1,'first'); startTime = julian2datenum(interpTime(startInd)); interval = mean(diff(julian2datenum(interpTime)))*24; ttideOutFile = fullfile(plotDirectory,['ttide_curr_site',site,'_',num2str(targetDepths(b)),'mbs.txt']); [nameu,fu,tidecon,xout]=t_tide(interp_uv,'start time',startTime,... 'interval',interval,... 'latitude',latitude,... 'output',ttideOutFile); end nc1=close(nc1); nc2=close(nc2); end