%% Script to create inputs for Duck Delft3D %% Import flowGrid = load('D:\Fateme_Alex\postRap\Hermine\C7HR100.mat'); % flowGrid = load('D:\Alexander\RT_25m\C7_25.mat'); % flowGrid = load('D:\NCBathy\Cone7\C9S_50.mat'); % waveGrid = load('D:\Alexander\RealtimeHRRR.mat'); % waveGrid = load('D:\NCBathy\Cone7\C9S_100.mat'); load('D:\Temp\NC_Data.mat') load('measuredNOAA.mat'); load('D:\Fateme_Alex\postRap\Hermine\NCOM2016Final.mat') clear NC_Data_PAR NC_Data_NOAA NC_Data_WQ %% Process % axisLimits=[datetime(2016,10,07,00,00,00) datetime(2016,10,12)]; % axisLimits=[datetime(2016,08,27,00,00,00) datetime(2016,09,05)]; % axisLimits=[datetime(2016,08,15,00,00,00) datetime(2016,11,01)]; axisLimits=[datetime(2019,09,02) datetime(2019,09,09,18,00,00)]; dateStart = axisLimits(1); dateEnd = axisLimits(2); delfTime = axisLimits(1):minutes(10):axisLimits(2); clear windIN windOut rainOut tempOut tideOut wlOut windOut_p_20 windOut_m_20 waveOut %% Process Wind for Site=[1 3 4 5] windIN(1,Site,:) = interp1(datenum(datetime(NC_Data_Wind{Site}(:,4),'ConvertFrom','posixtime')),... NC_Data_Wind{Site}(:,5),... datenum(axisLimits(1):minutes(10):axisLimits(2))); windIN(2,Site,:) = interp1(datenum(datetime(NC_Data_Wind{Site}(:,4),'ConvertFrom','posixtime')),... NC_Data_Wind{Site}(:,9),... datenum(axisLimits(1):minutes(10):axisLimits(2)),'nearest'); end windIN(windIN == 0) = NaN; delfTime = axisLimits(1):minutes(10):axisLimits(2); windOut(1,:) = nanmean(windIN(1,:,:),2); windOut(2,:) = wrapTo360(rad2deg(circ_mean(deg2rad(windIN(2,:,:)), [], 2))); windOut_p_20(1,:) = nanmean(windIN(1,:,:),2); windOut_p_20(2,:) = wrapTo360(rad2deg(circ_mean(deg2rad(windIN(2,:,:)+20), [], 2))); windOut_m_20(1,:) = nanmean(windIN(1,:,:),2); windOut_m_20(2,:) = wrapTo360(rad2deg(circ_mean(deg2rad(windIN(2,:,:)-20), [], 2))); %% Process Pier Wind [UA,~,idx] = unique(NC_Data_pierWind(:,4)); windSpeedU = [accumarray(idx,NC_Data_pierWind(:,5),[],@mean)]; windDirU = [accumarray(idx,NC_Data_pierWind(:,9),[],@mean)]; windIN_P(1,:) = interp1(datenum(datetime(UA,'ConvertFrom','posixtime')),... windSpeedU,... datenum(axisLimits(1):minutes(10):axisLimits(2)),'nearest'); windIN_P(2,:) = interp1(datenum(datetime(UA,'ConvertFrom','posixtime')),... windDirU,... datenum(axisLimits(1):minutes(10):axisLimits(2)),'nearest'); delfTime = axisLimits(1):minutes(10):axisLimits(2); windOut_P(1,:) = windIN_P(1,:); windOut_P(2,:) = windIN_P(2,:); clear UA windSpeedU windDirU windIN %% Process Rain % Rain is given in total per 10 minutes, so multiplied by 6 to become % hourly rainOut = interp1(datenum(datetime(NC_Data_Rain(:,3),'ConvertFrom','posixtime')),... NC_Data_Rain(:,5),... datenum(axisLimits(1):minutes(10):axisLimits(2))).*6; %% Process Temperature tempOut = interp1(datenum(datetime(NC_Data_Temp(:,1),'ConvertFrom','posixtime')),... NC_Data_Temp(:,4),... datenum(axisLimits(1):minutes(10):axisLimits(2))); %% Process Tides stations{6} = 8651370; %Dck stations{9} = 8656483; %Beau % monthDay = [31,30,31]; for Stat=[6 9] for Month = 1:3 if Month==1 dateStartM=dateStart +((Month-1).*days(31)); dateEndM=dateStart +((Month).*days(31)); else dateStartM=dateEndM + days(1); dateEndM=dateStart +((Month).*days(31)); end url = ['https://tidesandcurrents.noaa.gov/api/datagetter?product=water_level&application=NOS.COOPS.TAC.WL&'... 'begin_date=' datestr(dateStartM,'yyyymmdd') '&end_date=' datestr(dateEndM,'yyyymmdd') '&datum=NAVD' '&station=' num2str(stations{Stat})... '&time_zone=GMT&units=metric&format=csv']; wlIN = webread(url); if Month==1 measuredWL_D = datenum(wlIN.DateTime); measuredWL_W = wlIN.WaterLevel; else measuredWL_D = vertcat(measuredWL_D,datenum(wlIN.DateTime)); measuredWL_W = vertcat(measuredWL_W,wlIN.WaterLevel); end clear windIN wlIN end switch Stat case 6 tideOutN = interp1(measuredWL_D,... measuredWL_W,... datenum(axisLimits(1):minutes(10):axisLimits(2))); case 9 tideOutS = interp1(measuredWL_D,... measuredWL_W,... datenum(axisLimits(1):minutes(10):axisLimits(2))); end end %% Process Waves waveOut(:,1) = interp1(datenum(datetime(NC_Data_WAVE{3}(:,5),'ConvertFrom','posixtime')),... NC_Data_WAVE{3}(:,6),... datenum(axisLimits(1):minutes(10):axisLimits(2))); waveOut(:,2) = interp1(datenum(datetime(NC_Data_WAVE{3}(:,5),'ConvertFrom','posixtime')),... 1./NC_Data_WAVE{3}(:,7),... datenum(axisLimits(1):minutes(10):axisLimits(2))); waveOut(:,3) = windOut(2,:)'; waveOut(:,4) = ones(length(waveOut),1).*15; %% Import NDBC Waves % py_addpath('D:\MATLAB4\Downloads\wave_buoys_exploit-master') % py.importlib.import_module('utils_xuan') % py_addpath('D:\MATLAB4\Downloads\wave_buys_exploit-master\buoy_original_wanghe_source_code_buoystromassociator_31mai2018\buoy\utilos') % py.importlib.import_module('wav_utils') py_addpath('D:\MATLAB4\Downloads\DirectionalSpectra-master') mod = py.importlib.import_module('getmem'); addpath('Python') freqNDBC = [0.02 0.0325 0.0375 0.0425 0.0475 0.0525 0.0575 0.0625 0.0675 0.0725 0.0775 0.0825 0.0875 0.0925 0.1 0.11 0.12 0.13 0.14 0.15 0.16 0.17 0.18 0.19 0.2 0.21 0.22 0.23 0.24 0.25 0.26 0.27 0.28 0.29 0.3 0.31 0.32 0.33 0.34 0.35 0.365 0.385 0.405 0.425 0.445 0.465 0.485]; StatCount=1; clear buoyIn dirSpec tmpSpc % datesWW3=dateStart:minutes(20):dateEnd; for Station = [217 41025 44014 147] VarCount=1; if ismember(Station,[44014,41025])==1 for Var = [23 4 9 10 11] [buoyIn{StatCount,1}(:,:,VarCount),buoyIn{StatCount,2}]... = NDBCsave(['D:\WW3\' num2str(Station) alphabet(Var) '2016.txt']); % [buoyIn{StatCount,1}(:,:,VarCount),buoyIn{StatCount,2}]... % = NDBCsave(['D:\WW3\' num2str(Station) alphabet(Var) '2019.txt']); VarCount=VarCount+1; end else url = ['http://thredds.cdip.ucsd.edu/thredds/dodsC/cdip/archive/' num2str(Station) 'p1/' num2str(Station) 'p1_historic.nc']; % url = ['http://thredds.cdip.ucsd.edu/thredds/dodsC/cdip/realtime/' num2str(Station) 'p1_rt.nc']; waveTime = datetime(ncread(url,'waveTime'),'ConvertFrom','POSIX'); cDipIDX(1) = NearestValue(datenum(dateStart),datenum(waveTime)); cDipIDX(2) = NearestValue(datenum(dateEnd),datenum(waveTime)); freqCdip = ncread(url,'waveFrequency'); buoyIn{StatCount,1}(:,:,1) = ncread(url,'waveEnergyDensity',[1 cDipIDX(1)],[64 cDipIDX(2)-cDipIDX(1)])'; buoyIn{StatCount,1}(:,:,2) = ncread(url,'waveA1Value',[1 cDipIDX(1)],[64 cDipIDX(2)-cDipIDX(1)])'; buoyIn{StatCount,1}(:,:,3) = ncread(url,'waveA2Value',[1 cDipIDX(1)],[64 cDipIDX(2)-cDipIDX(1)])'; buoyIn{StatCount,1}(:,:,4) = ncread(url,'waveB1Value',[1 cDipIDX(1)],[64 cDipIDX(2)-cDipIDX(1)])'; buoyIn{StatCount,1}(:,:,5) = ncread(url,'waveB2Value',[1 cDipIDX(1)],[64 cDipIDX(2)-cDipIDX(1)])'; buoyIn{StatCount,2} = datetime(ncread(url,'waveTime',[cDipIDX(1)],[cDipIDX(2)-cDipIDX(1)]),'ConvertFrom','POSIX'); end DirCount=1; directions = [5:5:360]; directions2 = [1:1:360]; directionsR = deg2rad(directions); clear tmpSpc tmpSpc if ismember(Station,[44014,41025])==1 for Direction = directions tmpSpc(:,:,DirCount) = buoyIn{StatCount,1}(:,:,1).*... (1/180).*(0.5+0.01.*buoyIn{StatCount,1}(:,:,4).*cosd(1.*(Direction-buoyIn{StatCount,1}(:,:,2)))+... (0.01.*buoyIn{StatCount,1}(:,:,5).*cosd(2.*(Direction-buoyIn{StatCount,1}(:,:,3))))); %https://www.ndbc.noaa.gov/measdes.shtml DirCount=DirCount+1; end else tmpSpc = zeros(length(buoyIn{StatCount,2}),length(freqCdip),length(directions2)); for Time = 1:length(buoyIn{StatCount,2}) pySpc = mod.GetMem(py.numpy.array(squeeze(buoyIn{StatCount}(Time,:,2))),... py.numpy.array(squeeze(buoyIn{StatCount}(Time,:,4))),py.numpy.array(squeeze(buoyIn{StatCount}(Time,:,3))),... py.numpy.array(squeeze(buoyIn{StatCount}(Time,:,5)))); tmpSpc(Time,:,:) = np2mat(pySpc).*repmat(squeeze(buoyIn{StatCount}(Time,:,1))',1,size(np2mat(pySpc),2)); disp(Time) end clear pySpc end if ismember(Station,[44014,41025])~=1 [meshIN1,meshIN2] = meshgrid(freqCdip,directions2); [meshOUT1,meshOUT2] = meshgrid(freqNDBC',directions); meshIN1=double(meshIN1); dirSpec{StatCount,1} = zeros(length(buoyIn{StatCount,2}),length(freqNDBC),length(directions)); for Time = 1:length(buoyIn{StatCount,2}) dirSpec{StatCount,1}(Time,:,:) = interp2(meshIN1,meshIN2,squeeze(tmpSpc(Time,:,:))',... meshOUT1,meshOUT2)'; end dirSpec{StatCount,1}(isnan(dirSpec{StatCount,1})==1)=0; dirSpec{StatCount,2}=buoyIn{StatCount,2}; else dirSpec{StatCount,1}=tmpSpc; dirSpec{StatCount,2}=buoyIn{StatCount,2}; end clear tmpSpc StatCount=StatCount+1; end clear Var Direction DirCount StatCount VarCount %% Processs WW3 waves % [meshDateQ,mesgFreqQ,meshDirQ] = meshgrid(freq,datenum(datesWW3),directions); clear swanIN swanINTER % for StatCount=1:2 dateCount=1; % datesWW3=datetime(2016,09,02):minutes(20):datetime(2016,09,05); datesWW3=datetime(2016,10,06):minutes(20):datetime(2016,10,11); % datesWW3=datetime(2019,09,02):minutes(60):datetime(2019,09,05); swanINTER=zeros(length(datesWW3),length(freqNDBC),length(directions),7); for dateWW3 = datesWW3 [idx(1),waveDif(1)] = NearestValue(dateWW3,dirSpec{1,2}); [idx(2),waveDif(2)] = NearestValue(dateWW3,dirSpec{2,2}); [idx(3),waveDif(3)] = NearestValue(dateWW3,dirSpec{3,2}); [idx(4),waveDif(4)] = NearestValue(dateWW3,dirSpec{4,2}); if max(waveDif)>hours(1) disp(['noData for ' datestr(dateWW3)]) end % [meshDate,mesgFreq,meshDir] = meshgrid(freq,datenum(dirSpec{StatCount,2}),directions); % swanIN(:,:,:,StatCount) = interp3(meshDate,mesgFreq,meshDir,... % dirSpec{StatCount,1},meshDateQ,mesgFreqQ,meshDirQ); % swanIN(dateCount,:,:,StatCount) = dirSpec{StatCount,1}(idx,:,:); % [Xq,Yq,Zq] = meshgrid( 1:size(dirSpec{StatCount,1},3),1:size(dirSpec{StatCount,1},2),... % linspace(1,2,size(flowGrid.data.Y,2)./2)); % swanINTER(dateCount,:,:,:) = interp3(... % cat(3,squeeze(dirSpec{1,1}(idx(1),:,:)),squeeze(dirSpec{2,1}(idx(2),:,:))),... % Xq,Yq,Zq); % swanINTER(dateCount,:,:,:) = linspaceNDim(squeeze(dirSpec{1,1}(idx(1),:,:)),... % squeeze(dirSpec{2,1}(idx(2),:,:)),floor(size(flowGrid.data.Y,2)./25)); swanINTER(dateCount,:,:,1:3) = linspaceNDim(squeeze(dirSpec{1,1}(idx(1),:,:)),... squeeze(dirSpec{2,1}(idx(2),:,:)),3); swanINTER(dateCount,:,:,3:7) = linspaceNDim(squeeze(dirSpec{2,1}(idx(2),:,:)),... squeeze(dirSpec{3,1}(idx(3),:,:)),5); swanINTER(dateCount,:,:,7:8) = linspaceNDim(squeeze(dirSpec{3,1}(idx(3),:,:)),... squeeze(dirSpec{4,1}(idx(4),:,:)),2); dateCount=dateCount+1; disp(dateWW3) end % end % xPts = [1300 2161 2161 1650]; % yPts = [2 2 2641 2641]; % xPts = [650 1081 1081 825]; % yPts = [2 2 1321 1321]; % xPts = [1300 1749lear akll 1749 1650]; % yPts = [2 2 2641 2641]; xPts = [490 874 874 740]; yPts = [2 2 1321 1321]; xPtsSWAN(1:3) = floor(linspace(xPts(1),xPts(2),3)); xPtsSWAN(3:7) = floor(linspace(xPts(2),xPts(3),5)); xPtsSWAN(7:8) = floor(linspace(xPts(3),xPts(4),2)); yPtsSWAN(1:3) = floor(linspace(yPts(1),yPts(2),3)); yPtsSWAN(3:7) = floor(linspace(yPts(2),yPts(3),5)); yPtsSWAN(7:8) = floor(linspace(yPts(3),yPts(4),2)); for i=1:8 xPtsSWANutm(i) = waveGrid.data.X(xPtsSWAN(i),yPtsSWAN(i)); yPtsSWANutm(i) = waveGrid.data.Y(xPtsSWAN(i),yPtsSWAN(i)); end % dwr_2d_swan(freq,directions,swanIN(:,:,:,1),'ObsN',datesWW3,[5.1431e+05],[4.0512e+06],1,'NDBC_N.bcw') dwr_2d_swanV2(freqNDBC,directions,swanINTER,'NDBC',datesWW3,... [xPtsSWANutm],[yPtsSWANutm],1,'NDBC_2016Oct_C9S.bcw') % dwr_2d_swan(freq,directions,swanIN(:,:,:,2),'ObsS',datesWW3,[4.6689e+05],[3.8759e+06],1,'NDBC_S.bcw') %% Process WL for Site=[1 2 3 5] wlIN(Site,:) = interp1(datenum(datetime(NC_Data_WL{Site}(:,4),'ConvertFrom','posixtime')),... NC_Data_WL{Site}(:,5),... datenum(axisLimits(1):minutes(10):axisLimits(2))); end wlIN(wlIN == 0) = NaN; delfTime = axisLimits(1):minutes(10):axisLimits(2); wlOut = nanmean(wlIN(:,:)); %% load Measured NOAA WL load('measuredNOAA.mat') tideOutN = interp1(measuredWL{6}(:,1),... measuredWL{6}(:,2),... datenum(axisLimits(1):minutes(10):axisLimits(2))); tideOutS = interp1(measuredWL{9}(:,1),... measuredWL{9}(:,2),... datenum(axisLimits(1):minutes(10):axisLimits(2))); %% Process Stream % Designate Input Coordinates % streamDelftLocation=[325409 4034751;... %1 % 307140 3937494;... %2 % 296176 3903604;... %3 % 315241 3886276;... %4 % 388920 4024591;... %5 % 317623 3970513]; %6 % streamDelftNum = [5 1 6 2 3 4 1 1 1]; streamDelftLocation=[325409 4034751;... %1 313151 3935550;... %2 296176 3903604;... %3 315241 3886276;... %4 388920 4024591;... %5 322564 3971902]; %6 streamDelftNum = [5 1 6 2 3 4 1 1 1]; for Stream=1:9 streamIN(Stream,:) = interp1(measuredStream{Stream}(:,1),... measuredStream{Stream}(:,2),... datenum(axisLimits(1):minutes(10):axisLimits(2))); end for Stream=1:6 streamOut(:,Stream) = nansum(streamIN(streamDelftNum==Stream,:),1); end streamOut(streamOut<0)=0; %% Load Grid and find stream locations grid = load('D:\Fateme_Alex\postRap\Hermine\C7HR100.mat'); for Stream=1:6 [streamIDX(Stream,:)] = NearestValue(streamDelftLocation(Stream,:),grid.data.X,grid.data.Y); end %% Wind fileID=fopen('D:\Fateme_Alex\CS2019\windOut_m_20.wnd','w'); refTime=0; for line=1:length(delfTime) fprintf(fileID,'%d. %2.1f %3.0f\n',refTime,windOut_m_20(1,line),windOut_m_20(2,line)); refTime=refTime+10; end fclose(fileID) %% Wind P fileID=fopen('D:\Fateme_Alex\CS2019\WindPH.wnd','w'); refTime=0; for line=1:length(delfTime) fprintf(fileID,'%d. %2.1f %3.0f\n',refTime,windOut_P(1,line),windOut_P(2,line)); refTime=refTime+10; end fclose(fileID) %% Rain fileID=fopen('D:\Fateme_Alex\CS2019\RainFull.eva','w'); refTime=0; for line=1:length(delfTime) fprintf(fileID,'%d. %2.2f -999.0 %2.2f\n',refTime,rainOut(line),tempOut(line)); refTime=refTime+10; end fclose(fileID) %% Streams fileID=fopen('D:\Fateme_Alex\CS2019\Stream_100.src','w'); for line=1:length(streamIDX) fprintf(fileID,'Stream_%d Y %d %d 0 N\n',line,streamIDX(line,1),streamIDX(line,2)); end fclose(fileID) fileID=fopen('D:\Fateme_Alex\CS2019\Stream_100.dis','w'); refTime=0; records=length(delfTime); for Stream=1:length(streamIDX) refTime=0; header=[... 'table-name ''Discharge : %d''\n'... 'contents ''regular ''\n'... 'location ''Stream_%d ''\n'... 'time-function ''non-equidistant''\n'... 'reference-time %d\n'... 'time-unit ''minutes''\n'... 'interpolation ''linear''\n'... 'parameter ''time '' unit ''[min]''\n'... 'parameter ''flux/discharge rate'' unit ''[m3/s]''\n'... 'records-in-table %d\n'... ]; fprintf(fileID,header,Stream,Stream,str2num(datestr(axisLimits(1),'yyyymmdd')),records) for line=1:length(delfTime) fprintf(fileID,'%d %1.7f\n',refTime,streamOut(line,Stream)); refTime=refTime+10; end end fclose(fileID) %% Obs %% Create Observation File load('CSpoints.mat'); locations = [36.606,-74.840;... 36.260,-75.594;... 36.001,-75.421;... 35.750,-75.330;... 35.025,-75.363; 36.189244, -75.739155;... 36.3205,-75.872;... 36.1656,-75.8155;... 36.0866, -75.76833;... 35.8964389,-75.6220889;... 34.7113528,-76.7368333;... 36.20002, -75.7151]; locNames = {'Virginia Beach Wave',... 'Duck 26 Wave',... 'Nages Head Wave',... 'Oregon Inlet Wave',... 'Diamond Shoals Wave',... 'Duck 11 AWAC',... 'CS North',... 'CS South',... 'USGS CURRITUCK SOUND NR POINT HARBOR',... 'USGS ROANOKE SOUND AT POND ISLAND',... 'USGS BOGUE SOUND AT ATLANTIC BEACH',... 'Duck 17 Wave'}; % Realtime Points % Create Observation File locations = [36.606,-74.840;... 36.260,-75.594;... 36.001,-75.421;... 35.750,-75.330;... 35.025,-75.363; 36.189244, -75.739155;... 36.3205,-75.872;... 36.1656,-75.8155;... 36.0866, -75.76833;... 35.8964389,-75.6220889;... 34.7113528,-76.7368333;... 35.96176,-75.64163;... 36.04318,-75.68936;... 36.12932,-75.74501;... 36.20002, -75.7151]; locNames = {'Virginia Beach Wave',... 'Duck 26 Wave',... 'Nages Head Wave',... 'Oregon Inlet Wave',... 'Diamond Shoals Wave',... 'Duck 11 AWAC',... 'CS North',... 'CS South',... 'USGS CURRITUCK SOUND NR POINT HARBOR',... 'USGS ROANOKE SOUND AT POND ISLAND',... 'USGS BOGUE SOUND AT ATLANTIC BEACH',... 'Villa Dunes SW Dock Shoreline',... 'Hayman St Dock SW Shoreline',... 'Gunther Telephone Pole'... 'Duck 17 Wave'}; % if FIRST_RUN == 0 for i = 1:length(locations) pointsCS{i+11,1} = locNames{i}; [pointsCS{i+11,2},pointsCS{i+11,3}] = wgs2utm(locations(i,1),locations(i,2),18,'N'); end pointsCS{27,1} = 'Beaufort 2'; pointsCS{27,2} = 345563; pointsCS{27,3} = 3842033; pointsCS{27,1} = 'Beaufort 3'; pointsCS{27,2} = 346832; pointsCS{27,3} = 3838721; pointsCS{8,3} = pointsCS{8,3} + 2500; pointsCS{7,2} = pointsCS{7,2} - 500; fileID=fopen('obs_f_RT.obs','w'); for line=1:length(pointsCS) obsIDX = NearestValue([pointsCS{line,2:3}],flowGrid.data.X,flowGrid.data.Y); if length(pointsCS{line,1})>20 fprintf(fileID,[pointsCS{line,1}(1:20) ' %d %d\n'],obsIDX(1),obsIDX(2)); else fprintf(fileID,[pointsCS{line,1}(1:end) ' %d %d\n'],obsIDX(1),obsIDX(2)); end end fclose(fileID); fileID=fopen('obs_w_RT.loc','w'); for i = 1:length(pointsCS) fprintf(fileID,' %1.7f %1.7f\n',pointsCS{i,2},pointsCS{i,3}); end fclose(fileID) %end %% Currents from NCOM idxCount=1; for i = round(linspace(2,length(flowGrid.data.Y)-1,53)) % Every 5 km % [gridDeg(idxCount,1),gridDeg(idxCount,2)] = utm2deg(flowGrid.data.X(1081,i),flowGrid.data.Y(1081,i),'18 S'); % [gridDeg(idxCount,1),gridDeg(idxCount,2)] = utm2deg(flowGrid.data.X(2161,i),flowGrid.data.Y(2161,i),'18 S'); [gridDeg(idxCount,1),gridDeg(idxCount,2)] = utm2deg(flowGrid.data.X(4322,i),flowGrid.data.Y(4322,i),'18 S'); [idxNCOM(idxCount,1)] = NearestValue(gridDeg(idxCount,1),NCOM.Lat); %Lat [idxNCOM(idxCount,2)] = NearestValue(gridDeg(idxCount,2),wrapTo180(NCOM.Lon)); %Lon % crIDX(idxCount,:) = [1082,i]; % crIDX(idxCount,:) = [2161,i]; crIDX(idxCount,:) = [4322,i]; idxCount = idxCount + 1; end crBoundCount(1) = 0; crBoundCount(2) = idxCount -1; % for i = round(linspace(744,1080,14)) % Every 5 km % for i = round(linspace(1488,2161,14)) % Every 5 km for i = round(linspace(2976,2161,14)) % Every 5 km % [gridDeg(idxCount,1),gridDeg(idxCount,2)] = utm2deg(flowGrid.data.X(i,1321),flowGrid.data.Y(i,1321),'18 S'); % [gridDeg(idxCount,1),gridDeg(idxCount,2)] = utm2deg(flowGrid.data.X(i,2641),flowGrid.data.Y(i,2641),'18 S'); [gridDeg(idxCount,1),gridDeg(idxCount,2)] = utm2deg(flowGrid.data.X(i,5282),flowGrid.data.Y(i,5282),'18 S'); [idxNCOM(idxCount,1)] = NearestValue(gridDeg(idxCount,1),NCOM.Lat); %Lat [idxNCOM(idxCount,2)] = NearestValue(gridDeg(idxCount,2),wrapTo180(NCOM.Lon)); %Lon % crIDX(idxCount,:) = [i,1321]; % crIDX(idxCount,:) = [i,2641]; crIDX(idxCount,:) = [i,5282]; idxCount = idxCount + 1; end crBoundCount(3) = idxCount -1; % for i = round(linspace(455,1080,26)) % Every 5 km % for i = round(linspace(910,2161,26)) % Every 5 km for i = round(linspace(1820,4322,26)) % Every 5 km % [gridDeg(idxCount,1),gridDeg(idxCount,2)] = utm2deg(flowGrid.data.X(i,1),flowGrid.data.Y(i,1),'18 S'); [gridDeg(idxCount,1),gridDeg(idxCount,2)] = utm2deg(flowGrid.data.X(i,1),flowGrid.data.Y(i,1),'18 S'); [idxNCOM(idxCount,1)] = NearestValue(gridDeg(idxCount,1),NCOM.Lat); %Lat [idxNCOM(idxCount,2)] = NearestValue(gridDeg(idxCount,2),wrapTo180(NCOM.Lon)); %Lon % crIDX(idxCount,:) = [i,1]; crIDX(idxCount,:) = [i,1]; idxCount = idxCount + 1; end crBoundCount(4) = idxCount -1; clear idxCount gridDeg %% Interpolate Water Level for tStep = 1:length(tideOutN) wlO(tStep,1:crBoundCount(2)) = linspace(tideOutS(tStep),tideOutN(tStep),crBoundCount(2)); wlO(tStep,crBoundCount(2)+1:crBoundCount(3)) = linspace(tideOutN(tStep),tideOutN(tStep),crBoundCount(3)-crBoundCount(2)); wlO(tStep,crBoundCount(3)+1:crBoundCount(4)) = linspace(tideOutS(tStep),tideOutS(tStep),crBoundCount(4)-crBoundCount(3)); end %% Match to grid flowBed = load('D:\Fateme_Alex\postRap\Hermine\C7HR50_Bed.mat'); idxCount=1; for i = 1:length(idxNCOM) currX = squeeze(NCOM.davX(idxNCOM(idxCount,2),idxNCOM(idxCount,1),1,:)); currY = squeeze(NCOM.davY(idxNCOM(idxCount,2),idxNCOM(idxCount,1),1,:)); currX40 = squeeze(NCOM.water_u(idxNCOM(idxCount,2),idxNCOM(idxCount,1),:,:)); currY40 = squeeze(NCOM.water_v(idxNCOM(idxCount,2),idxNCOM(idxCount,1),:,:)); currDir(i,:) = wrapTo360(atan2d(-currY,currX)+90); currDir40(i,:,:) = wrapTo360(atan2d(-currY40,currX40)+90); xGridDir(i,:) = ([flowGrid.data.X(crIDX(idxCount,1)-2,crIDX(idxCount,2)-0) flowGrid.data.X(crIDX(idxCount,1)-1,crIDX(idxCount,2)-0)]); yGridDir(i,:) = ([flowGrid.data.Y(crIDX(idxCount,1)-2,crIDX(idxCount,2)-0) flowGrid.data.Y(crIDX(idxCount,1)-1,crIDX(idxCount,2)-0)]); wlN(:,idxCount) = squeeze(NCOM.surf_el(idxNCOM(idxCount,2),idxNCOM(idxCount,1),:))+0.112; mag(:,idxCount) = vecmag(currX,currY); DM(idxCount) = -squeeze(flowBed.data.Val(crIDX(idxCount,1)-1,crIDX(idxCount,2))); if i>crBoundCount(2) && i<=crBoundCount(3) dirGrid(i) = wrapTo360(atan2d(-diff(yGridDir(i,:)),diff(xGridDir(i,:)))+90); currents(:,idxCount) = sind(dirGrid(i)-currDir(i,:))'.*vecmag(currX,currY); RN_1(:,idxCount) = currents(:,idxCount) - wlN(:,idxCount).*sqrt(9.81./DM(idxCount)); currents_2(:,idxCount) = currY; for j=1:40 currents_40(:,j,idxCount) = squeeze(sind(dirGrid(i)-squeeze(currDir40(i,j,:))).*vecmag(currX40(j,:),currY40(j,:))'); end elseif i>crBoundCount(3) dirGrid(i) = wrapTo360(atan2d(-diff(yGridDir(i,:)),diff(xGridDir(i,:)))+90); currents(:,idxCount) = sind(dirGrid(i)-currDir(i,:))'.*vecmag(currX,currY); RN_1(:,idxCount) = currents(:,idxCount) + wlN(:,idxCount).*sqrt(9.81./DM(idxCount)); currents_2(:,idxCount) = currY; for j=1:40 currents_40(:,j,idxCount) = squeeze(sind(dirGrid(i)-squeeze(currDir40(i,j,:))).*vecmag(currX40(j,:),currY40(j,:))'); end else dirGrid(i) = wrapTo360(atan2d(-diff(yGridDir(i,:)),diff(xGridDir(i,:)))); currents(:,idxCount) = -sind(dirGrid(i)-currDir(i,:))'.*vecmag(currX,currY); RN_1(:,idxCount) = currents(:,idxCount) - wlN(:,idxCount).*sqrt(9.81./DM(idxCount)); currents_2(:,idxCount) = currX; for j=1:40 currents_40(:,j,idxCount) = squeeze(-sind(dirGrid(i)-squeeze(currDir40(i,j,:))).*vecmag(currX40(j,:),currY40(j,:))'); end end % disp(idxCount) idxCount = idxCount + 1; end % fprintf(fileID,[boundryNames{Boundry} ' N T 750 1322 1081 1322 0.0000000e+000\n']); % fprintf(fileID,[boundryNames{Boundry} ' N T 475 1 1081 1 0.0000000e+000\n']); clear idxCount %% interpolate Currents clear currentsOut for i = 1:3 [xIN,yIN] = meshgrid([crBoundCount(i)+1:crBoundCount(i+1)],... datenum(hours(NCOM.time)+datetime(2013,04,05))); [xOut,yOut] = meshgrid([crBoundCount(i)+1:crBoundCount(i+1)],... datenum(axisLimits(1):minutes(10):axisLimits(2))); currentsOut(:,crBoundCount(i)+1:crBoundCount(i+1)) = ... interp2(xIN,yIN,... currents(:,crBoundCount(i)+1:crBoundCount(i+1)),... xOut,yOut); currents2Out(:,crBoundCount(i)+1:crBoundCount(i+1)) = ... interp2(xIN,yIN,... currents_2(:,crBoundCount(i)+1:crBoundCount(i+1)),... xOut,yOut); wlNOut(:,crBoundCount(i)+1:crBoundCount(i+1)) = ... interp2(xIN,yIN,... wlN(:,crBoundCount(i)+1:crBoundCount(i+1)),... xOut,yOut); RNOut_1(:,crBoundCount(i)+1:crBoundCount(i+1)) = ... interp2(xIN,yIN,... RN_1(:,crBoundCount(i)+1:crBoundCount(i+1)),... xOut,yOut); for j = 1:length(idxNCOM) currents40Out(:,crBoundCount(i)+1:crBoundCount(i+1)) = ... interp2(xIN,yIN,... squeeze(currents_40(:,j,crBoundCount(i)+1:crBoundCount(i+1)),... xOut,yOut)); end end idxCount=1; for i = 1:length(idxNCOM) if i>crBoundCount(2) && i<=crBoundCount(3) RNOut_2(:,i) = currentsOut(:,idxCount) - wlO(:,idxCount).*sqrt(9.81./DM(idxCount)); elseif i>crBoundCount(3) RNOut_2(:,i) = currentsOut(:,idxCount) + wlO(:,idxCount).*sqrt(9.81./DM(idxCount)); else RNOut_2(:,i) = currentsOut(:,idxCount) - wlO(:,idxCount).*sqrt(9.81./DM(idxCount)); end idxCount = idxCount+1; end %% Combined Boundry Files records=size(tideOutS,2); boundryNames(1) = {'OceanE'}; for i = 2:crBoundCount(4) boundryNames(i) = {['NCOM' num2str(i)]}; end boundryTypes(1) = {'water elevation (z)'}; boundryTypes(2:crBoundCount(4)+1) = {'current (c)'}; boundryUnits(1) = {'[m]'}; boundryUnits(2:crBoundCount(4)+1) = {'[m/s]'}; fileID=fopen('NCOM_Tides2.bct','w'); for Boundry = [1 crBoundCount(2)+1:length(boundryNames)] refTime=0; if Boundry==1 elseif ismember(Boundry-1,crBoundCount)==1 continue end header=[... 'table-name ''Boundary Section : %d''\n'... 'contents ''Uniform ''\n'... 'location ''' boundryNames{Boundry} ' ''\n'... 'time-function ''non-equidistant''\n'... 'reference-time %d\n'... 'time-unit ''minutes''\n'... 'interpolation ''linear''\n'... 'parameter ''time '' unit ''[min]''\n'... 'parameter ''' boundryTypes{Boundry} ' end A'' unit ''' boundryUnits{Boundry} '''\n'... 'parameter ''' boundryTypes{Boundry} ' end B'' unit ''' boundryUnits{Boundry} '''\n'... 'records-in-table %d\n'... ]; fprintf(fileID,header,Boundry,str2num(datestr(dateStart,'yyyymmdd')),records); if Boundry<2 for line = 1:size(tideOutS,2) fprintf(fileID,'%d %1.4f %1.4f\n',refTime,tideOutS(line),tideOutN(line)); refTime=refTime+10; end else for line = 1:size(tideOutS,2) fprintf(fileID,'%d %1.4f %1.4f\n',refTime,currents2Out(line,Boundry-1),currents2Out(line,Boundry)); refTime=refTime+10; end end end fclose(fileID); %% R Boundry Files records=size(tideOutS,2); for i = 1:crBoundCount(4) boundryNames(i) = {['NCOM' num2str(i)]}; end % boundryTypes(1) = {'water elevation (z)'}; boundryTypes(1:crBoundCount(4)) = {'riemann (r)'}; boundryUnits(1:crBoundCount(4)) = {'[m/s]'}; fileID=fopen('NCOM_R1.bct','w'); for Boundry = [2:length(boundryNames)] refTime=0; if Boundry==1 elseif ismember(Boundry-1,crBoundCount)==1 continue end header=[... 'table-name ''Boundary Section : %d''\n'... 'contents ''Uniform ''\n'... 'location ''' boundryNames{Boundry} ' ''\n'... 'time-function ''non-equidistant''\n'... 'reference-time %d\n'... 'time-unit ''minutes''\n'... 'interpolation ''linear''\n'... 'parameter ''time '' unit ''[min]''\n'... 'parameter ''' boundryTypes{Boundry} ' end A'' unit ''' boundryUnits{Boundry} '''\n'... 'parameter ''' boundryTypes{Boundry} ' end B'' unit ''' boundryUnits{Boundry} '''\n'... 'records-in-table %d\n'... ]; fprintf(fileID,header,Boundry,str2num(datestr(dateStart,'yyyymmdd')),records); for line = 1:size(tideOutS,2) fprintf(fileID,'%d %1.4f %1.4f\n',refTime,RNOut_1(line,Boundry-1),RNOut_1(line,Boundry)); refTime=refTime+10; end end fclose(fileID); %% Boundry definition file fileID=fopen('NCOM_25m.bnd','w'); crIDX2=crIDX; crIDX2(crBoundCount(2)+1:crBoundCount(3),2) = crIDX(crBoundCount(2)+1:crBoundCount(3),2)+1; for Boundry = [1 crBoundCount(2)+1:length(boundryNames)] if Boundry == 1 fprintf(fileID,[boundryNames{Boundry} ' Z T 1082 1320 1082 2 1.0000000e+004\n']); elseif ismember(Boundry-1,crBoundCount)==1 continue else boundStart = crIDX2(Boundry-1,:); boundEnd = crIDX2(Boundry,:); fprintf(fileID,[boundryNames{Boundry} ' C T ' num2str(boundStart(1)) ' ' num2str(boundStart(2)) ' '... num2str(boundEnd(1)) ' ' num2str(boundEnd(2)) ' 1.0000000e+004 Uniform \n']); end end fclose(fileID); %% Boundry definition file-R fileID=fopen('NCOM_R25.bnd','w'); crIDX2=crIDX; % crIDX2(crBoundCount(2)+1:crBoundCount(3),2) = crIDX(crBoundCount(2)+1:crBoundCount(3),2)+1; for Boundry = [2:length(boundryNames)] if ismember(Boundry-1,crBoundCount)==1 continue else boundStart = crIDX2(Boundry-1,:); boundEnd = crIDX2(Boundry,:); fprintf(fileID,[boundryNames{Boundry} ' R T ' num2str(boundStart(1)) ' ' num2str(boundStart(2)) ' '... num2str(boundEnd(1)) ' ' num2str(boundEnd(2)) ' 0.0000000e+000 Uniform \n']); end end fclose(fileID); %% Tides refTime=0; records=length(delfTime); boundryNames = {'East','OceanN','OceanS'}; tideOutNSmooth=smooth(tideOutN,10); tideOutSSmooth=smooth(tideOutS,10); fileID=fopen('D:\Fateme_Alex\CSBoundryTS_SmoothNS_NOAA_NS_100.bct','w'); for Boundry =[2 1 3] refTime=0; header=[... 'table-name ''Boundary Section : %d''\n'... 'contents ''Uniform ''\n'... 'location ''' boundryNames{Boundry} ' ''\n'... 'time-function ''non-equidistant''\n'... 'reference-time %d\n'... 'time-unit ''minutes''\n'... 'interpolation ''linear''\n'... 'parameter ''time '' unit ''[min]''\n'... 'parameter ''water elevation (z) end A'' unit ''[m]''\n'... 'parameter ''water elevation (z) end B'' unit ''[m]''\n'... 'records-in-table %d\n'... ]; fprintf(fileID,header,Boundry,str2num(datestr(axisLimits(1),'yyyymmdd')),records) for line=1:length(delfTime) switch Boundry case 1 fprintf(fileID,'%d %1.4f %1.4f\n',refTime,tideOutN(line),tideOutS(line)); case 2 fprintf(fileID,'%d %1.4f %1.4f\n',refTime,tideOutN(line),tideOutN(line)); case 3 fprintf(fileID,'%d %1.4f %1.4f\n',refTime,tideOutS(line),tideOutS(line)); end refTime=refTime+10; end end fclose(fileID) refTime=0; records=length(delfTime); boundryNames = {'East','OceanN','OceanS'}; boundryTypes = {'water elevation (z)','neumann (n)','neumann (n)'}; tideOutNSmooth=zeros(length(tideOutN),3); tideOutSSmooth=zeros(length(tideOutN),3); tideOutNSmooth(:,1)=smooth(tideOutN,10); tideOutSSmooth(:,1)=smooth(tideOutS,10); fileID=fopen('D:\Fateme_Alex\CSBoundryTS_SmoothNS_NOAA_NS_Neu_100.bct','w'); crBoundCount=1; for Boundry =[2 1 3] refTime=0; header=[... 'table-name ''Boundary Section : %d''\n'... 'contents ''Uniform ''\n'... 'location ''' boundryNames{Boundry} ' ''\n'... 'time-function ''non-equidistant''\n'... 'reference-time %d\n'... 'time-unit ''minutes''\n'... 'interpolation ''linear''\n'... 'parameter ''time '' unit ''[min]''\n'... 'parameter ''' boundryTypes{Boundry} ' end A'' unit ''[m]''\n'... 'parameter ''' boundryTypes{Boundry} ' end B'' unit ''[m]''\n'... 'records-in-table %d\n'... ]; fprintf(fileID,header,crBoundCount,str2num(datestr(axisLimits(1),'yyyymmdd')),records) for line=1:length(delfTime) switch Boundry case 1 fprintf(fileID,'%d %1.4f %1.4f\n',refTime,tideOutN(line),tideOutS(line)); case 2 fprintf(fileID,'%d %1.4f %1.4f\n',refTime,0,0); case 3 fprintf(fileID,'%d %1.4f %1.4f\n',refTime,0,0); end refTime=refTime+10; end crBoundCount=crBoundCount+1; end fclose(fileID) % % refTime=0; % records=length(delfTime); % boundryNames = {'North','South'}; % boundryTypes = {'water elevation (z)','water elevation (z)'}; % % tideOutSmooth(:,1)=smooth(wlIN(3,:),25); % tideOutSmooth(:,2)=smooth(wlIN(1,:),25); % % fileID=fopen('D:\Fateme_Alex\CS2019\CS_Small_Bound.bct','w'); % % for Boundry =1:2 % refTime=0; % % header=[... % 'table-name ''Boundary Section : %d''\n'... % 'contents ''Uniform ''\n'... % 'location ''' boundryNames{Boundry} ' ''\n'... % 'time-function ''non-equidistant''\n'... % 'reference-time %d\n'... % 'time-unit ''minutes''\n'... % 'interpolation ''linear''\n'... % 'parameter ''time '' unit ''[min]''\n'... % 'parameter ''' boundryTypes{Boundry} ' end A'' unit ''[m]''\n'... % 'parameter ''' boundryTypes{Boundry} ' end B'' unit ''[m]''\n'... % 'records-in-table %d\n'... % ]; % % % fprintf(fileID,header,Boundry,str2num(datestr(axisLimits(1),'yyyymmdd')),records) % % for line=1:length(delfTime) % fprintf(fileID,'%d %1.4f %1.4f\n',refTime,tideOutSmooth(line,Boundry),tideOutSmooth(line,Boundry)); % refTime=refTime+10; % end % % % end % fclose(fileID) %% Waves refTime=0; records=length(delfTime); boundryNames = {'NorthWave'}; fileID=fopen('D:\Fateme_Alex\CS2019\CS_Wave_Bound.bcw','w'); for Boundry =1 refTime=0; header=[... 'location ''' boundryNames{Boundry} ' ''\n'... 'time-function ''non-equidistant''\n'... 'reference-time %d\n'... 'time-unit ''minutes''\n'... 'interpolation ''linear''\n'... 'parameter ''time'' unit ''[min]''\n'... 'parameter ''WaveHeight'' unit ''[m]''\n'... 'parameter ''Period'' unit ''[s]''\n'... 'parameter ''Direction'' unit ''[N^o]''\n'... 'parameter ''DirSpreading'' unit ''[degrees]''\n'... ]; fprintf(fileID,header,str2num(datestr(axisLimits(1),'yyyymmdd'))) for line=1:length(delfTime) if max(isnan(waveOut(line,:)))==0 fprintf(fileID,'%d %1.2f %1.2f %3.0f %d\n',refTime,waveOut(line,1),waveOut(line,2),waveOut(line,3),waveOut(line,4)); end refTime=refTime+10; end end fclose(fileID)