1022 lines
42 KiB
Matlab
1022 lines
42 KiB
Matlab
%% 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
|