%% NC DUNEX Realtime script % Alexander Rey, August 12, 2019 clear all %cd /project/6008948/alex0042/RealtimeH2 addpath('NCtoolbox') addpath('m_map') setup_nctoolbox flowGrid = load('Cone7W_100.mat'); waveGrid = load('Cone7W_250.mat'); for masterTime=datetime(2019,09,02,00,00,00):hours(6):datetime(2019,09,08,00,00,00) masterDate = datetime(year(masterTime),month(masterTime),day(masterTime)); if masterTime==datetime(2019,09,02,00,00,00) FIRST_RUN = 1; else FIRST_RUN = 0; end %% Clean files from previous run delete com-* delete hot* delete rt_run.prt* delete tri-diag.* delete TMP_* delete trim-* delete wavm-* delete trih-* delete obs.loct00* delete WNDNOW delete swan.inp delete PRINT* delete hot* delete SWANOUT* delete Errfile* delete WAVE2FLOW* delete FLOW2WAVE* clear HRRR_UTM if exist(['ModelIn_4' datestr(masterTime,'YYYY-MM-DD_hh') '.mat'])==2 load(['ModelIn_4' datestr(masterTime,'YYYY-MM-DD_hh') '.mat']) waterLevels(1:37,:) = waterLevels(1:37,:); end if masterTime>datetime(2019,09,15,00,00,00) %% Import data from NOMADS %% Extratropical Surge and Tide Operational Forecast System for the Atlantic % (ESTOFS Atlantic) for water levels url = ['http://nomads.ncep.noaa.gov/dods/estofs_atl/'... datestr(masterTime,'yyyymmdd')... '/estofs_atl_conus_' datestr(masterTime,'hh') 'z']; ESTOFS = ncgeodataset(url); clear url lonESTOFS=ESTOFS{'lon'}(:); latESTOFS=ESTOFS{'lat'}(:); % Find ESTOFS index for grid points along boundry idxCount=1; for i = 2:50:length(flowGrid.data.Y)-1 % Every 5 km [gridDeg(idxCount,1),gridDeg(idxCount,2)] = utm2deg(flowGrid.data.X(1081,i),flowGrid.data.Y(1081,i),'18 S'); [idxESTOFS(idxCount,1)] = NearestValue(gridDeg(idxCount,1),latESTOFS); %Lat [idxESTOFS(idxCount,2)] = NearestValue(gridDeg(idxCount,2),lonESTOFS); %Lon wlIDX(idxCount,:) = [1082,i]; idxCount = idxCount + 1; end wlBoundCount(1) = 0; wlBoundCount(2) = idxCount -1; for i = 750:55:1080 % Every 5 km [gridDeg(idxCount,1),gridDeg(idxCount,2)] = utm2deg(flowGrid.data.X(i,1321),flowGrid.data.Y(i,1321),'18 S'); [idxESTOFS(idxCount,1)] = NearestValue(gridDeg(idxCount,1),latESTOFS); %Lat [idxESTOFS(idxCount,2)] = NearestValue(gridDeg(idxCount,2),lonESTOFS); %Lon wlIDX(idxCount,:) = [i,1322]; idxCount = idxCount + 1; end wlBoundCount(3) = idxCount -1; for i = 475:55:1080 % Every 5 km [gridDeg(idxCount,1),gridDeg(idxCount,2)] = utm2deg(flowGrid.data.X(i,1),flowGrid.data.Y(i,1),'18 S'); [idxESTOFS(idxCount,1)] = NearestValue(gridDeg(idxCount,1),latESTOFS); %Lat [idxESTOFS(idxCount,2)] = NearestValue(gridDeg(idxCount,2),lonESTOFS); %Lon wlIDX(idxCount,:) = [i,1]; idxCount = idxCount + 1; end wlBoundCount(4) = idxCount -1; clear idxCount gridDeg % Extract ESTOFS Water Levels for 36 hours at indexes waterLevelsGRID = ESTOFS{'etcwlsfc'}(2:38,min(idxESTOFS(:,1)):max(idxESTOFS(:,1)),... min(idxESTOFS(:,2)):max(idxESTOFS(:,2))); % Extract small grid % Find ESTOFS values from small grid idxCount=1; for i = 1:length(idxESTOFS) waterLevels(:,idxCount) = waterLevelsGRID(:,idxESTOFS(idxCount,1)-min(idxESTOFS(:,1))+1,... idxESTOFS(idxCount,2)-min(idxESTOFS(:,2))+1); 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 disp('Loaded ESTOFS') %% Multi-grid Wave model (WAVEWATCH III) % 4 arc minute US East Coast grid for waves url = ['http://nomads.ncep.noaa.gov/dods/wave/mww3/'... datestr(masterTime,'yyyymmdd')... '/multi_1.at_4m' datestr(masterTime,'yyyymmdd') '_' datestr(masterTime,'hh') 'z']; MWW3 = ncgeodataset(url); clear url lonMWW3=wrapTo180(MWW3{'lon'}(:)); latMWW3=MWW3{'lat'}(:); % Find MWW3 index for grid points along boundry idxCount=1; for i = 1:25:length(waveGrid.data.Y)-25 % Every ~200 km [gridDeg(idxCount,1),gridDeg(idxCount,2)] = utm2deg(waveGrid.data.X(end-2,i),waveGrid.data.Y(end-2,i),'18 S'); waveDistAlong(idxCount) = sum(sqrt(diff(waveGrid.data.X(end-2,1:i)).^2+... diff(waveGrid.data.Y(end-2,1:i)).^2)); waveUTM(idxCount,1) = waveGrid.data.X(end-2,i); waveUTM(idxCount,2) = waveGrid.data.Y(end-2,i); [idxMWW3(idxCount,1)] = NearestValue(gridDeg(idxCount,1),latMWW3); %Lat [idxMWW3(idxCount,2)] = NearestValue(gridDeg(idxCount,2),lonMWW3); %Lon idxCount = idxCount + 1; end waveBoundCount(1) = 0; waveBoundCount(2) = idxCount -1; for i = size(waveGrid.data.X,1)-25:-25:350% Every ~200 km [gridDeg(idxCount,1),gridDeg(idxCount,2)] = utm2deg(waveGrid.data.X(i,end-2),waveGrid.data.Y(i,end-2),'18 S'); waveDistAlong(idxCount) = sum(sqrt(diff(waveGrid.data.X(i:end-2,end-2)).^2+... diff(waveGrid.data.Y(i:end-2,end-2)).^2)); waveUTM(idxCount,1) = waveGrid.data.X(i,end-2); waveUTM(idxCount,2) = waveGrid.data.Y(i,end-2); [idxMWW3(idxCount,1)] = NearestValue(gridDeg(idxCount,1),latMWW3); %Lat [idxMWW3(idxCount,2)] = NearestValue(gridDeg(idxCount,2),lonMWW3); %Lon idxCount = idxCount + 1; end waveBoundCount(3) = idxCount -1; for i = size(waveGrid.data.X,1)-25:-25:240 % Every ~200 km [gridDeg(idxCount,1),gridDeg(idxCount,2)] = utm2deg(waveGrid.data.X(i,2),waveGrid.data.Y(i,2),'18 S'); waveDistAlong(idxCount) = sum(sqrt(diff(waveGrid.data.X(i:end-2,1)).^2+... diff(waveGrid.data.Y(i:end-2,1)).^2)); waveUTM(idxCount,1) = waveGrid.data.X(i,2); waveUTM(idxCount,2) = waveGrid.data.Y(i,2); [idxMWW3(idxCount,1)] = NearestValue(gridDeg(idxCount,1),latMWW3); %Lat [idxMWW3(idxCount,2)] = NearestValue(gridDeg(idxCount,2),lonMWW3); %Lon idxCount = idxCount + 1; end waveBoundCount(4) = idxCount -1; clear idxCount % Extract WW3 Hs for 36 hours at indexes wavesGRID_Hs = MWW3{'htsgwsfc'}(1:13,min(idxMWW3(:,1)):max(idxMWW3(:,1)),... min(idxMWW3(:,2)):max(idxMWW3(:,2))); % Extract small grid htsgwsfc %swell_1 % Period wavesGRID_P = MWW3{'perpwsfc'}(1:13,min(idxMWW3(:,1)):max(idxMWW3(:,1)),... % min(idxMWW3(:,2)):max(idxMWW3(:,2))); % Extract small grid perpwsfc %swper_1 % Direction wavesGRID_Dir = MWW3{'dirpwsfc'}(1:13,min(idxMWW3(:,1)):max(idxMWW3(:,1)),... min(idxMWW3(:,2)):max(idxMWW3(:,2))); % Extract small grid dirpwsfc %swdir_1 % Find WW3 values from small grid idxCount=1; for i = 1:length(idxMWW3) waves(:,idxCount,1) = wavesGRID_Hs(:,idxMWW3(idxCount,1)-min(idxMWW3(:,1))+1,... idxMWW3(idxCount,2)-min(idxMWW3(:,2))+1); %Hs waves(:,idxCount,2) = wavesGRID_P(:,idxMWW3(idxCount,1)-min(idxMWW3(:,1))+1,... idxMWW3(idxCount,2)-min(idxMWW3(:,2))+1); %Tp waves(:,idxCount,3) = wavesGRID_Dir(:,idxMWW3(idxCount,1)-min(idxMWW3(:,1))+1,... idxMWW3(idxCount,2)-min(idxMWW3(:,2))+1); %Dir idxCount = idxCount + 1; end clear idxCount disp('Loaded WW3') end if masterTime>datetime(2019,09,00,00,00,00) %% HRRR (High-Resolution Rapid Refresh) % Wind Speeds, Directions, Pressures, Precipitations % if day(masterTime)<10 % url = ['http://nomads.ncep.noaadatestr(masterTime,'yyyymmdd').gov/dods/nam/nam'... % datestr(masterTime,'yyyymmdd')... % '/nam1hr_' datestr(masterTime,'hh') 'z']; % url = ['/scratch/alex0042/HRRR/' datestr(masterTime,'yyyymmddhh') '.g2']; % untar([url '.tar'],['/home/cluster/RealtimeH/HAS011383739/nam_218_' datestr(masterTime,'yyyymmddhh')]); % else % url = ['http://nomads.ncep.noaa.gov/dods/hrrr/hrrr'... % datestr(masterTime,'yyyymmdd')... % '/hrrr_sfc.t' datestr(masterTime,'hh') 'z']; % end HRRR = ncgeodataset(['D:\Alexander\HRRR_Forecast\' ... datestr(masterTime,'yyyymmdd') ... '_' datestr(masterTime,'hh') '00_000.grb2.regrid']); % Extract Lat Lon from HRRR % lonHRRR=wrapTo180(HRRR{'lon'}(:)); % latHRRR=HRRR{'lat'}(:); m_proj('lambert conformal conic','ori',[-95 25.0],'clo',-95,'par',[25 25],'ell','sphere'); [meshX,meshY] = meshgrid(HRRR{'x'}(:),HRRR{'y'}(:)); [GridLONAll_R,GridLATAll_R] = m_xy2ll(meshX.*1000,meshY.*1000); % [latHRRR_Mesh,lonHRRR_Mesh] = meshgrid(GridLATAll_R,GridLONAll_R); %Create Mesh % Find Flow Grid IDXs with margin minGrid(1) = min(min(flowGrid.data.X))-50000; minGrid(2) = min(min(flowGrid.data.Y))-50000; maxGrid(1) = max(max(flowGrid.data.X))+50000; maxGrid(2) = max(max(flowGrid.data.Y))+50000; % [TEMPdeg(1),TEMPdeg(2)] = utm2deg(maxGrid(1),maxGrid(2),'18 T'); % maxIDX = NearestValue(TEMPdeg,latHRRR_Mesh',lonHRRR_Mesh'); % % [TEMPdeg(1),TEMPdeg(2)] = utm2deg(minGrid(1),minGrid(2),'18 T'); % minIDX = NearestValue(TEMPdeg,latHRRR_Mesh',lonHRRR_Mesh'); % clear TEMPdeg [TEMPdeg(1),TEMPdeg(2)] = utm2deg(maxGrid(1),maxGrid(2),'18 T'); maxIDX = NearestValue(TEMPdeg,GridLATAll_R,GridLONAll_R); [TEMPdeg(1),TEMPdeg(2)] = utm2deg(minGrid(1),minGrid(2),'18 T'); minIDX = NearestValue(TEMPdeg,GridLATAll_R,GridLONAll_R); clear TEMPdeg % Extraction index latIDX_HRRR = minIDX(1)-1:maxIDX(1)+1; lonIDX_HRRR = minIDX(2)-1:maxIDX(2)+1; latHRRR_Y = HRRR{'y'}(latIDX_HRRR); lonHRRR_X = HRRR{'x'}(lonIDX_HRRR); % latHRRR_Grid = HRRR{'x'}(latIDX_HRRR); % lonHRRR_Grid = HRRR{'y'}(lonIDX_HRRR); clear meshX meshY [meshX,meshY] = meshgrid(HRRR{'x'}(lonIDX_HRRR)',HRRR{'y'}(latIDX_HRRR)'); [lonHRRR_Grid,latHRRR_Grid] = m_xy2ll(meshX.*1000,meshY.*1000); clear latHRRR_Mesh minIDX maxIDX % Final inport grid for HRRR data windGridLON = lonHRRR_Grid; windGridLAT = latHRRR_Grid; % if FIRST_RUN == 1 % % Extract shoreline in model domain % [windGridLON_HR,windGridLAT_HR] = meshgrid(linspace(wrapTo180(lonHRRR_Grid(1)),wrapTo180(lonHRRR_Grid(end)),100),... % linspace(wrapTo180(latHRRR_Grid(1)),wrapTo180(latHRRR_Grid(end)),100)); % % shoreline = gshhs('gshhs_h.b',double([min(min(windGridLAT_HR)) max(max(windGridLAT_HR))]),... % double([min(min(windGridLON_HR)) max(max(windGridLON_HR))])); % clear windGridLON_HR % save('shoreline.mat','shoreline') % else load('shoreline.mat'); % end clear windU windV pressure rain temperature wind if day(masterTime)<12 for tStep = 1:37 if tStep ~=100 % Extract MWW3 Hs for 36 hours at indexes HRRR = ncgeodataset(['D:\Alexander\HRRR_Forecast\' ... datestr(masterTime,'yyyymmdd') ... '_' datestr(masterTime,'hh') '00_0' num2str(tStep-1,'%02d') '.grb2.regrid']); windU(tStep,:,:) = HRRR{'u-component_of_wind_height_above_ground'}(1,1,latIDX_HRRR,lonIDX_HRRR); %U Winds @ 10m [m/s] windV(tStep,:,:) = HRRR{'v-component_of_wind_height_above_ground'}(1,1,latIDX_HRRR,lonIDX_HRRR); %V Winds @ 10m [m/s] pressure(tStep,:,:) = HRRR{'Pressure_surface'}(1,latIDX_HRRR,lonIDX_HRRR); %Surface Pressure [pa] if tStep==1 rain(tStep,:,:) = HRRR{'Total_precipitation_surface_0_Hour_Accumulation'}(1,latIDX_HRRR,lonIDX_HRRR); %Surface precipitation [kg/m^2] elseif tStep==2 rain(tStep,:,:) = HRRR{['Total_precipitation_surface_'... num2str(tStep-1) '_Hour_Accumulation']}(1,latIDX_HRRR,lonIDX_HRRR); else rain(tStep,:,:) = HRRR{['Total_precipitation_surface_Mixed_intervals_Accumulation']}(2,latIDX_HRRR,lonIDX_HRRR); end wind(tStep,:,:) = sqrt(windU(tStep,:,:).^2+windV(tStep,:,:).^2); % Extract Temp temperature(tStep,:) = squeeze(nanmean(nanmean(HRRR{'Temperature_height_above_ground'}(1,1,latIDX_HRRR,lonIDX_HRRR)- 273.15,2),3)); %Surface Temperature [C] % rain2(tStep,:,:) = HRRR{['Precipitation_rate_surface']}(1,latIDX_HRRR,lonIDX_HRRR); else windU(tStep,:,:) = windU(tStep-1,:,:); windV(tStep,:,:) = windV(tStep-1,:,:); pressure(tStep,:,:) = pressure(tStep-1,:,:); rain(tStep,:,:) = rain(tStep-1,:,:); temperature(tStep,:) = temperature(tStep-1,:); end end else % Extract MWW3 Hs for 36 hours at indexes windU = HRRR{'ugrd10m'}(1:37,latIDX_HRRR,lonIDX_HRRR); %U Winds @ 10m [m/s] windV = HRRR{'vgrd10m'}(1:37,latIDX_HRRR,lonIDX_HRRR); %V Winds @ 10m [m/s] pressure = HRRR{'pressfc'}(1:37,latIDX_HRRR,lonIDX_HRRR); %Surface Pressure [pa] rain = HRRR{'apcpsfc'}(1:37,latIDX_HRRR,lonIDX_HRRR); %Surface precipitation [kg/m^2] wind = HRRR{'wind10m'}(1:37,latIDX_HRRR,lonIDX_HRRR); %Surface precipitation [kg/m^2] % Extract Temp temperature = squeeze(nanmean(nanmean(HRRR{'tmpsfc'}(1:37,latIDX_HRRR,lonIDX_HRRR)- 273.15,2),3)); %Surface Temperature [C] end disp('Loaded HRRR') save(['ModelIn_HRRR_' datestr(masterTime,'YYYY-MM-DD_hh') '.mat'],'-v7.3','windU','windV','pressure','rain','wind','temperature','waves','waterLevels') disp('Loaded HRRR') else load('windGrid.mat'); end load('boundData.mat') load('shoreline.mat'); %% IMPORT NCOM NCOM.Lat = ncread(['D:\Alexander\NCOM\' datestr(masterDate,'YYYYmmDD') '00_t000.nc'],'lat',[331],[240]); NCOM.Lon = ncread(['D:\Alexander\NCOM\' datestr(masterDate,'YYYYmmDD') '00_t000.nc'],'lon',91,180); NCOM.depth = ncread(['D:\Alexander\NCOM\' datestr(masterDate,'YYYYmmDD') '00_t000.nc'],'depth'); Noffset = hours(masterTime - masterDate); tStep = Noffset; tCount = 1; for Ni = masterTime:hours(3):masterTime+hours(36) noData=1; while noData==1 NCOM.water_u(:,:,:,tCount) = ncread(['D:\Alexander\NCOM\' datestr(masterDate,'YYYYmmDD') '00_t0' num2str(tStep,'%02d') '.nc'],... 'water_u',[91 331 1 1],[180 240 inf 1],[1 1 1 1]); NCOM.water_v(:,:,:,tCount) = ncread(['D:\Alexander\NCOM\' datestr(masterDate,'YYYYmmDD') '00_t0' num2str(tStep,'%02d') '.nc'],... 'water_v',[91 331 1 1],[180 240 inf 1],[1 1 1 1]); NCOM.surf_el(:,:,tCount) = ncread(['D:\Alexander\NCOM\' datestr(masterDate,'YYYYmmDD') '00_t0' num2str(tStep,'%02d') '.nc'],... 'surf_el',[91 331 1],[180 240 1],[1 1 1]); NCOM.salinity(:,:,:,tCount) = ncread(['D:\Alexander\NCOM\' datestr(masterDate,'YYYYmmDD') '00_t0' num2str(tStep,'%02d') '.nc'],... 'salinity',[91 331 1 1],[180 240 inf 1],[1 1 1 1]); if (squeeze(NCOM.water_v(100,100,1,end))==0|squeeze(NCOM.water_u(100,100,1,end))==0|squeeze(NCOM.surf_el(100,100,end))==0|squeeze(NCOM.salinity(100,100,1,end))==0)==1 noData=1; else noData=0; end end NCOM.time(tCount) = Ni; tStep=tStep+3; tCount = tCount+1; end %% NCOM Process NCOM.mag = vecmag(NCOM.water_u,NCOM.water_v); depthWeight = [diff(NCOM.depth)' 0]; for Layer = 1:40 DAVfac(:,:,Layer,:) = ones(180,240,1,13).*depthWeight(Layer); end DAVtrue = DAVfac; DAVtrue(isnan(NCOM.mag)==1) = 0; DAVfac(isnan(NCOM.mag)==1) = nan; DAVFac2 = sum(DAVtrue,3); clear DAVtrue NCOM.dav = nansum(DAVfac.*NCOM.mag,3)./DAVFac2; NCOM.davX = nansum(DAVfac.*NCOM.water_u,3)./DAVFac2; NCOM.davY = nansum(DAVfac.*NCOM.water_v,3)./DAVFac2; %% Format data for Delft3D %% Evap fileID=fopen('TemperatureRT.eva','w'); refTime=minutes(masterTime-masterDate); records=37; for i = 1 : records fprintf(fileID,'%d. 0. 0. %2.2f\n',refTime,temperature(i)); refTime = refTime + 60;% 1 hour time steps for temperature end fclose(fileID); % %% Water Level % records=37; % % boundryNames = {'OceanN','OceanS'}; % for i = 1:wlBoundCount(4) % boundryNames(i) = {['East' num2str(i)]}; % end % % % boundryTypes = {'neumann (n)','neumann (n)'}; % boundryTypes(1:wlBoundCount(4)) = {'water elevation (z)'}; % % fileID=fopen('ESTOFS4.bct','w'); % for Boundry = 1:length(boundryNames) % refTime=minutes(masterTime-masterDate); % if ismember(Boundry,wlBoundCount)==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 ''[m]''\n'... % 'parameter ''' boundryTypes{Boundry} ' end B'' unit ''[m]''\n'... % 'records-in-table %d\n'... % ]; % fprintf(fileID,header,Boundry,str2num(datestr(masterDate,'yyyymmdd')),records); % % for line = 1:size(waterLevels,1) % fprintf(fileID,'%d %1.4f %1.4f\n',refTime,waterLevels(line,Boundry),waterLevels(line,Boundry+1)); % % refTime=refTime+60; % end % end % % fclose(fileID); % % Boundry definition file % if FIRST_RUN == 1 % fileID=fopen('ESTOFS4_Bound.bnd','w'); % % for Boundry =1 : length(boundryNames) % % if Boundry == 1 % % fprintf(fileID,[boundryNames{Boundry} ' N T 750 1322 1081 1322 0.0000000e+000\n']); % % elseif Boundry == 2 % % fprintf(fileID,[boundryNames{Boundry} ' N T 475 1 1081 1 0.0000000e+000\n']); % % else % if ismember(Boundry,wlBoundCount)==1 % continue % end % % boundStart = wlIDX(Boundry,:); % boundEnd = wlIDX(Boundry+1,:); % % fprintf(fileID,[boundryNames{Boundry} ' Z T ' num2str(boundStart(1)) ' ' num2str(boundStart(2)) ' '... % num2str(boundEnd(1)) ' ' num2str(boundEnd(2)) ' 1.0000000e+004\n']); % % end % % end % fclose(fileID); % end %% NCOM ESTOFS % Currents from NCOM idxCount=1; for i = 2:50:length(flowGrid.data.Y)-1 % 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'); [idxNCOM(idxCount,1)] = NearestValue(gridDeg(idxCount,1),NCOM.Lat); %Lat [idxNCOM(idxCount,2)] = NearestValue(gridDeg(idxCount,2),wrapTo180(NCOM.Lon)); %Lon crIDX(idxCount,:) = [1081,i]; % crIDX(idxCount,:) = [2161,i]; idxCount = idxCount + 1; end crBoundCount(1) = 0; crBoundCount(2) = idxCount -1; for i = 750:55:1080 % Every 5 km % for i = round(linspace(1488,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'); [idxNCOM(idxCount,1)] = NearestValue(gridDeg(idxCount,1),NCOM.Lat); %Lat [idxNCOM(idxCount,2)] = NearestValue(gridDeg(idxCount,2),wrapTo180(NCOM.Lon)); %Lon crIDX(idxCount,:) = [i,1320]; % crIDX(idxCount,:) = [i,2641]; idxCount = idxCount + 1; end crBoundCount(3) = idxCount -1; for i = 475:55:1080 % Every 5 km % for i = round(linspace(910,2161,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 for tStep = 1:37 wlO(tStep,:) = waterLevels(tStep,:); end %% Match to grid flowBed = load('D:\Alexander\RealtimeHRRR\Dec6_100Bed.mat'); clear wlN mag DM dirGrid currents RN_1 currents_2 currents_40 currDir currDir40 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(NCOM.time)); [xOut,yOut] = meshgrid([crBoundCount(i)+1:crBoundCount(i+1)],... datenum(masterTime:hours(1):masterTime+hours(36))); 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 %% R Boundry Files records=37; 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_R2.bct','w'); for Boundry = [2:length(boundryNames)] refTime=minutes(masterTime-masterDate); 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(masterDate,'yyyymmdd')),records); for line = 1:37 fprintf(fileID,'%d %1.4f %1.4f\n',refTime,RNOut_2(line,Boundry-1),RNOut_2(line,Boundry)); refTime=refTime+60; end end fclose(fileID); %% Boundry definition file-R fileID=fopen('NCOM_R.bnd','w'); crIDX2=crIDX; crIDX2(crBoundCount(2)+1:crBoundCount(3),2) = crIDX(crBoundCount(2)+1:crBoundCount(3),2)+2; crIDX2(crBoundCount(1)+1:crBoundCount(2),1) = 1082; 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); %% Waves records = size(waves,1); fileID=fopen('WW3.bcw','w'); for i = 1:3 boundryNames(i) = {['East' num2str(i)]}; end for Boundry = 1:3 refTime=minutes(masterTime-masterDate); header=[... 'location ''East%d''\n'... 'time-function ''non-equidistant''\n'... 'reference-time %d\n'... 'time-unit ''minutes''\n'... 'interpolation ''linear''\n'... 'parameter ''time '' unit ''[min]''\n'... ]; fprintf(fileID,header,Boundry,str2num(datestr(masterTime,'yyyymmdd'))); for Segment = waveBoundCount(Boundry)+1:waveBoundCount(Boundry+1) fprintf(fileID,'parameter ''WaveHeight '' unit ''[m]''\n'); end for Segment = waveBoundCount(Boundry)+1:waveBoundCount(Boundry+1) fprintf(fileID,'parameter ''Period '' unit ''[s]''\n'); end for Segment = waveBoundCount(Boundry)+1:waveBoundCount(Boundry+1) fprintf(fileID,'parameter ''Direction '' unit ''[N^o]''\n'); end for Segment = waveBoundCount(Boundry)+1:waveBoundCount(Boundry+1) fprintf(fileID,'parameter ''DirSpreading '' unit ''[-]''\n'); end for line = 1:records fprintf(fileID,'%d ',refTime); for Segment = waveBoundCount(Boundry)+1:waveBoundCount(Boundry+1) fprintf(fileID,'%1.4f ',waves(line,Segment,1)); %Hs end fprintf(fileID,'\n '); for Segment = waveBoundCount(Boundry)+1:waveBoundCount(Boundry+1) fprintf(fileID,'%1.4f ',waves(line,Segment,2)); %Period end fprintf(fileID,'\n '); for Segment = waveBoundCount(Boundry)+1:waveBoundCount(Boundry+1) fprintf(fileID,'%1.4f ',waves(line,Segment,3)); %Direction end fprintf(fileID,'\n '); for Segment = waveBoundCount(Boundry)+1:waveBoundCount(Boundry+1) fprintf(fileID,'%1.4f ',4); %Spread end fprintf(fileID,'\n'); % fprintf(fileID,'%.2f %1.4f %1.4f %1.4f %1.4f\n',refTime,... % waves(line,Boundry,1),waves(line,Boundry,2),waves(line,Boundry,3),4); refTime=refTime+180; end end fclose(fileID); %% Modify MDF file to change reference date fileID = fopen('rt_base.mdf'); cac = textscan( fileID, '%s', 'Delimiter','\n', 'CollectOutput',true ); fclose( fileID ); fileID = fopen('rt_run.mdf', 'w' ); refTime=minutes(masterTime-masterDate); for i = 1 : 15 fprintf( fileID, '%s\n', cac{1}{i} ); end fprintf( fileID, 'Itdate = #%04d-%02d-%02d#\n',... year(masterTime),month(masterTime),day(masterTime)); for i = 17 fprintf( fileID, '%s\n', cac{1}{i} ); end fprintf( fileID,['Tstart = %#.7e \n',... 'Tstop = %#.7e \n'],refTime,refTime+2160); for i = 20:29 fprintf( fileID, '%s\n', cac{1}{i} ); end if FIRST_RUN == 0 restID = ['tri-rst.rt_run.' datestr(masterTime,'yyyymmdd.HH') '0000']; % if masterTime~=datetime(2019,09,04,06,00,00) if exist(restID)==2 eval(['movefile ' restID ' tri-rst.6h_restart']) end fprintf( fileID,['Commnt = initial conditions from restart file\n'... 'Restid = #6h_restart#\n']); delete tri-rst.rt* else fprintf( fileID,'Zeta0 = %#.7e \n',nanmean(nanmean(squeeze(waterLevels(:,:))))); end for i = 31:94 fprintf( fileID, '%s\n', cac{1}{i} ); end fprintf( fileID,[... 'Flmap = %#.7e 60 %#.7e\n',... 'Flhis = %#.7e 10 %#.7e\n',... 'Flpp = %#.7e 60 %#.7e\n'],... refTime,refTime+2160,refTime,refTime+2160,refTime,refTime+2160); for i = 98:109 fprintf( fileID, '%s\n', cac{1}{i} ); end clear cac %% Modify MDW file to add times and change reference date fileID = fopen('rt_base.mdw'); cac = textscan( fileID, '%s', 'Delimiter','\n', 'CollectOutput',true ); fclose( fileID ); fileID = fopen('rt_run.mdw', 'w' ); refTime=minutes(masterTime-masterDate); for i = 1 : 7 fprintf( fileID, '%s\n', cac{1}{i} ); end fprintf( fileID, 'ReferenceDate = %04d-%02d-%02d\n',... year(masterTime),month(masterTime),day(masterTime)); for i = 9 : 71 fprintf( fileID, '%s\n', cac{1}{i} ); end % for i = 1 : records % fprintf( fileID,'[TimePoint] \n Time = %#.7e \n',(i-1).*180); % end bounOri{1}='East'; bounOri{2}='North'; bounOri{3}='South'; bounRot{1}='counter-clockwise'; bounRot{2}='counter-clockwise'; bounRot{3}='clockwise'; for Boundry = 1 : 3 fprintf(fileID,['[Boundary]\n'... 'Name = East%d\n'... 'Definition = orientation\n'... 'Orientation = %s\n'... 'SpectrumSpec = parametric\n'... 'SpShapeType = pierson-moskowitz\n'... 'GaussSpread = 9.9999998e-003\n'... 'PeriodType = peak\n'... 'DirSpreadType = power\n'... 'PeakEnhanceFac = 3.3000000e+000\n'... 'DistanceDir = %s\n'... ],Boundry,bounOri{Boundry},bounRot{Boundry}); for Segment = waveBoundCount(Boundry)+1:waveBoundCount(Boundry+1) fprintf(fileID,'CondSpecAtDist = %d\n',waveDistAlong(Segment)) end end clear cac fclose( fileID ); %% Met Files % Create new UTM grid of the HRRR Data [UTM1D(:,1),UTM1D(:,2)] = wgs2utm(reshape(windGridLAT,1,[]),reshape(windGridLON,1,[]),... 18,'N'); HRRR_UTM(:,:,1) = double(reshape(UTM1D(:,1),size(windGridLON,1),size(windGridLON,2))); HRRR_UTM(:,:,2) = double(reshape(UTM1D(:,2),size(windGridLON,1),size(windGridLON,2))); clear UTM1D G7LR = load('C7_1400_V2_C.mat'); % Remove land rainfall % if FIRST_RUN == 1 % landmask=zeros(size(G7LR.data.X(1:end-1,1:end-1)',1),size(G7LR.data.X(1:end-1,1:end-1)',2)); % [HR_deg(:,1),HR_deg(:,2)] = utm2deg(reshape(G7LR.data.X(1:end-1,1:end-1)',1,[]),reshape(G7LR.data.Y(1:end-1,1:end-1)',1,[]),repmat('18 S',length(reshape(G7LR.data.X(1:end-1,1:end-1)',1,[])),1)); % for ShoreLevel = 1:length(shoreline) % landmask_Temp(:,:) = reshape(inpolygon(HR_deg(:,1),HR_deg(:,2),shoreline(ShoreLevel).Lat,shoreline(ShoreLevel).Lon),size(G7LR.data.X(1:end-1,1:end-1)',1),size(G7LR.data.X(1:end-1,1:end-1)',2)); % % landmask = landmask|landmask_Temp; % clear landmask_Temp % % [shoreTempX,shoreTempY,shoreTempCell] = deg2utm(shoreline(ShoreLevel).Lat,shoreline(ShoreLevel).Lon); % shoreUTM{ShoreLevel}(:,1) = shoreTempX(all((shoreTempCell==repmat('18 S',length(shoreTempCell),1))')); % shoreUTM{ShoreLevel}(:,2) = shoreTempY(all((shoreTempCell==repmat('18 S',length(shoreTempCell),1))')); % clear shoreTempX shoreTempY shoreTempCell % disp(ShoreLevel) % end % save('shorelineInfo.mat','landmask','shoreUTM'); % else load('shorelineInfo.mat'); % end clear HR_deg ShoreLevel % Zero out met array HRRR_INTER = zeros(4,37,size(G7LR.data.X(1:end-1,1:end-1),2),size(G7LR.data.X(1:end-1,1:end-1),1)); % Interpolate against flowGrid for each time step for TS = 1:37 HRRR_INTER(1,TS,:,:) = griddata(HRRR_UTM(:,:,1),HRRR_UTM(:,:,2),... double(squeeze(windU(TS,:,:))),G7LR.data.X(1:end-1,1:end-1)',G7LR.data.Y(1:end-1,1:end-1)','natural'); HRRR_INTER(2,TS,:,:) = griddata(HRRR_UTM(:,:,1),HRRR_UTM(:,:,2),... double(squeeze(windV(TS,:,:))),G7LR.data.X(1:end-1,1:end-1)',G7LR.data.Y(1:end-1,1:end-1)','natural'); rainTemp = griddata(HRRR_UTM(:,:,1),HRRR_UTM(:,:,2),... double(squeeze(rain(TS,:,:))),G7LR.data.X(1:end-1,1:end-1)',G7LR.data.Y(1:end-1,1:end-1)','natural'); rainTemp(landmask)=rainTemp(landmask)*0; % 100% infiltration HRRR_INTER(3,TS,:,:) = rainTemp; % Convert to hours clear rainTemp HRRR_INTER(4,TS,:,:) = griddata(HRRR_UTM(:,:,1),HRRR_UTM(:,:,2),... double(squeeze(pressure(TS,:,:))),G7LR.data.X(1:end-1,1:end-1)',G7LR.data.Y(1:end-1,1:end-1)','natural'); end % Set nan to Delft3D zero value HRRR_INTER(isnan(HRRR_INTER)==1)=-999.000; % Create input files m = size(HRRR_INTER,3); n = size(HRRR_INTER,4); d3ddx = 5000; d3ddy = 5000; disp(m) d3d_input_uvpFieldsCL_FUN(masterTime,m,n,... squeeze(HRRR_INTER(4,:,:,:)),squeeze(HRRR_INTER(1,:,:,:)),squeeze(HRRR_INTER(2,:,:,:)),squeeze(HRRR_INTER(3,:,:,:)),1,'hrrr',37) save(['ModelIn_4_HRRR' datestr(masterTime,'YYYY-MM-DD_hh') '.mat'],'-v7.3','windU','windV','pressure','rain','wind','temperature','waves','waterLevels') %% Run Delft3D FLOW/WAVE! % system('. /home/cluster/D7545/enviorment.sh') try system('./runDuck_Wave.sh') eval(['copyfile trim-rt_run.dat ' datestr(masterTime,'YYYY-MM-DD_hh') '_mapOut4.dat']) eval(['movefile trim-rt_run.dat /project/6008948/alex0042/RealtimeH2/' datestr(masterTime,'YYYY-MM-DD_hh') '_mapOut4.dat']) eval(['copyfile trim-rt_run.def ' datestr(masterTime,'YYYY-MM-DD_hh') '_mapOut4.def']) eval(['movefile trim-rt_run.def /project/6008948/alex0042/RealtimeH2/' datestr(masterTime,'YYYY-MM-DD_hh') '_mapOut4.def']) eval(['copyfile wavm-rt_run.dat ' datestr(masterTime,'YYYY-MM-DD_hh') '_waveOut4.dat']) eval(['movefile wavm-rt_run.dat /project/6008948/alex0042/RealtimeH2/' datestr(masterTime,'YYYY-MM-DD_hh') '_waveOut4.dat']) eval(['copyfile wavm-rt_run.def ' datestr(masterTime,'YYYY-MM-DD_hh') '_waveOut4.def']) eval(['movefile wavm-rt_run.def /project/6008948/alex0042/RealtimeH2/' datestr(masterTime,'YYYY-MM-DD_hh') '_waveOut4.def']) eval(['copyfile trih-rt_run.dat ' datestr(masterTime,'YYYY-MM-DD_hh') '_histOut4.dat']) eval(['movefile trih-rt_run.dat /project/6008948/alex0042/RealtimeH2/' datestr(masterTime,'YYYY-MM-DD_hh') '_histOut4.dat']) eval(['copyfile trih-rt_run.def ' datestr(masterTime,'YYYY-MM-DD_hh') '_histOut4.def']) eval(['movefile trih-rt_run.def /project/6008948/alex0042/RealtimeH2/' datestr(masterTime,'YYYY-MM-DD_hh') '_histOut4.def']) eval(['copyfile ModelIn_4_HRRR' datestr(masterTime,'YYYY-MM-DD_hh') '.mat /project/6008948/alex0042/RealtimeH2/ModelIn_4' datestr(masterTime,'YYYY-MM-DD_hh') '.mat']) eval(['copyfile tri-rst.rt* /project/6008948/alex0042/RealtimeH2/']) catch break end end