2334 lines
108 KiB
Matlab
2334 lines
108 KiB
Matlab
%% NC DUNEX Realtime script
|
|
% Alexander Rey, August 12, 2019
|
|
% Moved to COASTLINES Feb 29, 2020
|
|
|
|
clear all
|
|
cd D:\DUNEX_RT\Operation
|
|
echo off;
|
|
addpath('NCtoolbox')
|
|
addpath('m_map')
|
|
setup_nctoolbox
|
|
flowGrid = load('Cone7W_100.mat');
|
|
waveGrid = load('Cone7W_250.mat');
|
|
load('LMSL_RT.mat')
|
|
addpath('D:\DUNEX_RT\Operation\abogaard-matlab2gsheets-6d76271')
|
|
setenv('PATH', [getenv('PATH') ';D:\DUNEX_RT\Windows_64']);
|
|
|
|
%% Details for e-mail server
|
|
setpref('Internet','SMTP_Server','smtp.sendgrid.net');
|
|
setpref('Internet','E_mail','mail@alexanderrey.ca');
|
|
setpref('Internet','SMTP_Username','apikey');
|
|
setpref('Internet','SMTP_Password','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');
|
|
props = java.lang.System.getProperties;
|
|
props.setProperty('mail.smtp.auth','true');
|
|
props.setProperty('mail.smtp.socketFactory.class', 'javax.net.ssl.SSLSocketFactory');
|
|
props.setProperty('mail.smtp.socketFactory.port','465');
|
|
props.setProperty('mail.smtp.starttls.enable','true');
|
|
|
|
|
|
%% Master Model Time Step
|
|
% Every 6 hours
|
|
|
|
nistTime = datetime(datevec(now())) - tzoffset(datetime(datevec(now()),'TimeZone','local'));
|
|
masterTime = dateshift(nistTime,'start','hour')-hours(2) % Assume script will run at 5 and a half hours after model start
|
|
clear nistTime
|
|
masterDate = datetime(year(masterTime),month(masterTime),day(masterTime));
|
|
|
|
if exist(['tri-rst.rt_run.' datestr(masterTime,'yyyymmdd.HH') '0000'])==2
|
|
FIRST_RUN = 0;
|
|
else
|
|
FIRST_RUN = 1;
|
|
end
|
|
|
|
%% Import
|
|
try
|
|
Import=0;
|
|
|
|
|
|
%% Clean files from previous run
|
|
system('taskkill /IM "d_hydro.exe" /F')
|
|
system('taskkill /IM "wave.exe" /F')
|
|
|
|
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
|
|
while Import==0
|
|
try
|
|
%% Import data from NOMADS
|
|
%% Extratropical Surge and Tide Operational Forecast System for the Atlantic
|
|
% (ESTOFS Atlantic) for water levels
|
|
options = weboptions;
|
|
options.Timeout = 90;
|
|
try
|
|
% url = ['https://nomads.ncep.noaa.gov/pub/data/nccf/com/estofs/prod/estofs.'...
|
|
% datestr(masterTime-hours(6),'yyyymmdd')...
|
|
% '/estofs.t' datestr(masterTime-hours(6),'hh') 'z.conus.east.f000.grib2'];
|
|
|
|
url = ['https://noaa-gestofs-pds.s3.amazonaws.com/estofs.'...
|
|
datestr(masterTime-hours(6),'yyyymmdd')...
|
|
'/estofs.t' datestr(masterTime-hours(6),'hh') 'z.conus.east.f000.grib2'];
|
|
|
|
websave('ESTOFS.grib',url,options)
|
|
ESTOFS = ncgeodataset('ESTOFS.grib');
|
|
|
|
m_proj('lambert conformal conic','ori',[-95 25.0],'clo',-95,'par',[25 25],'ell','sphere');
|
|
[meshX,meshY] = meshgrid(ESTOFS{'x'}(:),ESTOFS{'y'}(:));
|
|
[lonESTOFS,latESTOFS] = m_xy2ll(meshX.*1000,meshY.*1000);
|
|
clear meshX meshY
|
|
ESTstart = 7;
|
|
ESTrun = 6;
|
|
catch
|
|
% url = ['http://nomads.ncep.noaa.gov/dods/estofs_atl/'...
|
|
% datestr(masterTime-days(1),'yyyymmdd')...
|
|
% '/estofs_atl_conus_' datestr(masterTime-days(1),'hh') 'z'];
|
|
% ESTOFS = ncgeodataset(url);
|
|
% lonESTOFS=ESTOFS{'lon'}(:);
|
|
% latESTOFS=ESTOFS{'lat'}(:);
|
|
% ESTstart = 24+ESTstart;
|
|
% url = ['https://nomads.ncep.noaa.gov/pub/data/nccf/com/estofs/prod/estofs.'...
|
|
% datestr(masterTime-days(1),'yyyymmdd')...
|
|
% '/estofs.t' datestr(masterTime-days(1),'hh') 'z.conus.east.f000.grib2'];
|
|
url = ['https://noaa-gestofs-pds.s3.amazonaws.com/estofs.'...
|
|
datestr(masterTime-days(1),'yyyymmdd')...
|
|
'/estofs.t' datestr(masterTime-days(1),'hh') 'z.conus.east.f000.grib2'];
|
|
|
|
websave('ESTOFS.grib',url,options)
|
|
ESTOFS = ncgeodataset('ESTOFS.grib');
|
|
|
|
m_proj('lambert conformal conic','ori',[-95 25.0],'clo',-95,'par',[25 25],'ell','sphere');
|
|
[meshX,meshY] = meshgrid(ESTOFS{'x'}(:),ESTOFS{'y'}(:));
|
|
[lonESTOFS,latESTOFS] = m_xy2ll(meshX.*1000,meshY.*1000);
|
|
clear meshX meshY
|
|
ESTstart = 25;
|
|
ESTrun = 24;
|
|
end
|
|
|
|
|
|
clear url
|
|
|
|
% Find ESTOFS index for grid points along boundry
|
|
idxCount=1;
|
|
for i = 2:10: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,:)] = ...
|
|
NearestValue([gridDeg(idxCount,:)],latESTOFS,lonESTOFS); %Lat
|
|
wlIDX(idxCount,:) = [1082,i];
|
|
|
|
idxCount = idxCount + 1;
|
|
end
|
|
|
|
wlBoundCount(1) = 0;
|
|
wlBoundCount(2) = idxCount -1;
|
|
|
|
|
|
for i = 750:11: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
|
|
[idxESTOFS(idxCount,:)] = ...
|
|
NearestValue([gridDeg(idxCount,:)],latESTOFS,lonESTOFS);
|
|
wlIDX(idxCount,:) = [i,1322];
|
|
idxCount = idxCount + 1;
|
|
end
|
|
|
|
wlBoundCount(3) = idxCount -1;
|
|
|
|
for i = 475:5: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
|
|
[idxESTOFS(idxCount,:)] = ...
|
|
NearestValue([gridDeg(idxCount,:)],latESTOFS,lonESTOFS);
|
|
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{'Ocean_Surface_Elevation_Relative_to_Geoid_surface'}(ESTstart:ESTstart+36,min(idxESTOFS(:,1)):max(idxESTOFS(:,1)),...
|
|
% min(idxESTOFS(:,2)):max(idxESTOFS(:,2))); % Extract small grid
|
|
|
|
% Find ESTOFS values from small grid
|
|
idxCount=1;
|
|
errorET = 0;
|
|
for i = ESTstart:ESTstart+96
|
|
if errorET<10
|
|
try
|
|
% url = ['https://nomads.ncep.noaa.gov/pub/data/nccf/com/estofs/prod/estofs.'...
|
|
% datestr(masterTime-hours(ESTrun),'yyyymmdd')...
|
|
% '/estofs.t' datestr(masterTime-hours(ESTrun),'hh') 'z.conus.east.f'...
|
|
% num2str(i,'%03d') '.grib2'];
|
|
|
|
url = ['https://noaa-gestofs-pds.s3.amazonaws.com/estofs.'...
|
|
datestr(masterTime-hours(ESTrun),'yyyymmdd')...
|
|
'/estofs.t' datestr(masterTime-hours(ESTrun),'hh') 'z.conus.east.f'...
|
|
num2str(i,'%03d') '.grib2'];
|
|
websave(['D:\DUNEX_RT\Archive\ESTOFS\ESTOFS_' ...
|
|
datestr(masterTime,'YYYY-mm-DD_HH') '_f' num2str(idxCount-1,'%03d') '.grib'],url,options)
|
|
|
|
if i<=ESTstart+48
|
|
ESTOFS = ncgeodataset(['D:\DUNEX_RT\Archive\ESTOFS\ESTOFS_' ...
|
|
datestr(masterTime,'YYYY-mm-DD_HH') '_f' num2str(idxCount-1,'%03d') '.grib']);
|
|
|
|
waterLevelsGRID(idxCount,:,:) = ESTOFS{'Ocean_Surface_Elevation_Relative_to_Geoid_surface'}(1,min(idxESTOFS(:,1)):max(idxESTOFS(:,1)),...
|
|
min(idxESTOFS(:,2)):max(idxESTOFS(:,2))); % Extract small grid
|
|
end
|
|
idxCount = idxCount + 1;
|
|
catch
|
|
errorET = errorET+1;
|
|
end
|
|
else
|
|
error('ESTOFS ERROR')
|
|
end
|
|
|
|
end
|
|
|
|
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
|
|
|
|
|
|
% Save ESTOFS Points
|
|
% url = ['https://nomads.ncep.noaa.gov/pub/data/nccf/com/estofs/prod/estofs.'...
|
|
% datestr(masterTime-hours(ESTrun),'yyyymmdd')...
|
|
% '/estofs.t' datestr(masterTime-hours(ESTrun),'hh') 'z.points.cwl.nc'];
|
|
url = ['https://noaa-gestofs-pds.s3.amazonaws.com/estofs.'...
|
|
datestr(masterTime-hours(ESTrun),'yyyymmdd')...
|
|
'/estofs.t' datestr(masterTime-hours(ESTrun),'hh') 'z.points.cwl.nc'];
|
|
|
|
websave(['D:\DUNEX_RT\Archive\ESTOFS\ESTOFS_' ...
|
|
datestr(masterTime,'YYYY-mm-DD_HH') '_points.nc'],url,options)
|
|
% 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 = ['https://nomads.ncep.noaa.gov/dods/wave/gfswave/'...
|
|
datestr(masterTime-hours(6),'yyyymmdd')...
|
|
'/gfswave.atlocn.0p16_' datestr(masterTime-hours(6),'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
|
|
% 0 132 163 285
|
|
|
|
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'}(7:23,min(idxMWW3(:,1)):max(idxMWW3(:,1)),...
|
|
min(idxMWW3(:,2)):max(idxMWW3(:,2))); % Extract small grid htsgwsfc
|
|
%swell_1
|
|
% Period
|
|
wavesGRID_P = MWW3{'perpwsfc'}(7:23,min(idxMWW3(:,1)):max(idxMWW3(:,1)),... %
|
|
min(idxMWW3(:,2)):max(idxMWW3(:,2))); % Extract small grid perpwsfc
|
|
%swper_1
|
|
% Direction
|
|
wavesGRID_Dir = MWW3{'dirpwsfc'}(7:23,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')
|
|
|
|
|
|
%% HRRR (High-Resolution Rapid Refresh)
|
|
% Wind Speeds, Directions, Pressures, Precipitations
|
|
% url = ['http://nomads.ncep.noaa.gov/dods/hrrr/hrrr'...
|
|
% datestr(masterTime,'yyyymmdd')...
|
|
% '/hrrr_sfc.t' datestr(masterTime,'hh') 'z'];
|
|
% delete(gcp('nocreate'))
|
|
% parpool(8);
|
|
|
|
for ForcastHour = 0:48
|
|
% url = ['http://nomads.ncep.noaa.gov/pub/data/nccf/com/hrrr/prod/'...
|
|
% 'hrrr.' datestr(masterTime,'yyyymmdd')...
|
|
% '/conus/hrrr.t' num2str(hour(masterTime),'%02d') 'z.'...
|
|
% 'wrfsfcf' num2str(ForcastHour,'%02d') '.grib2'];
|
|
|
|
url = ['https://nomads.ncep.noaa.gov/cgi-bin/filter_hrrr_2d.pl?file=hrrr.'...
|
|
't' num2str(hour(masterTime),'%02d') 'z.wrfsfc'...
|
|
'f' num2str(ForcastHour,'%02d') '.grib2&'...
|
|
'lev_10_m_above_ground=on&lev_surface=on&var_APCP=on&var_PRES=on&var_UGRD=on&var_VGRD=on&var_TMP=on&subregion=&leftlon=-78&rightlon=-74&toplat=38&bottomlat=32.5&dir=%2F'...
|
|
'hrrr.' datestr(masterTime,'yyyymmdd') '%2Fconus'];
|
|
|
|
websave(['D:\DUNEX_RT\Operation\HRRR\hrrrT' num2str(ForcastHour,'%02d') '.grib2'],url,options)
|
|
|
|
[pCheck,gridParams] = system(['perl "D:\DUNEX_RT\Windows_64\grid_defn.pl" '...
|
|
'D:\DUNEX_RT\Operation\HRRR\hrrrT' num2str(ForcastHour,'%02d') '.grib2']);
|
|
[wCheck,gridParams] = system(['wgrib2 D:\DUNEX_RT\Operation\HRRR\hrrrT' ...
|
|
num2str(ForcastHour,'%02d') '.grib2 -new_grid_winds earth -new_grid ' ...
|
|
gridParams ' D:\DUNEX_RT\Operation\HRRR\hrrrT'...
|
|
num2str(ForcastHour,'%02d') '.grib2.regrid']);
|
|
end
|
|
% delete(gcp('nocreate'))
|
|
|
|
HRRR = ncgeodataset(['D:\DUNEX_RT\Operation\HRRR\hrrrT' num2str(0,'%02d') '.grib2.regrid']);
|
|
clear url
|
|
|
|
m_proj('lambert conformal conic','ori',[262.5 38.5],'clo',262.5,'par',[38.5 38.5],'ell','sphere');%,'lat',[20 55],'lon',[-60 -130]
|
|
[meshX,meshY] = meshgrid(HRRR{'x'}(:),HRRR{'y'}(:));
|
|
[windGridLON,windGridLAT] = m_xy2ll(meshX.*1000,meshY.*1000);
|
|
windGridLON=wrapTo180(windGridLON);
|
|
clear meshX meshY
|
|
|
|
% Extract MWW3 Hs for 48 hours at indexes
|
|
for ForcastHour = 0:48
|
|
|
|
HRRR = ncgeodataset(['D:\DUNEX_RT\Operation\HRRR\hrrrT' num2str(ForcastHour,'%02d') '.grib2.regrid']);
|
|
windU(ForcastHour+1,:,:) = HRRR{'u-component_of_wind_height_above_ground'}(:,:,:); %U Winds @ 10m [m/s]
|
|
windV(ForcastHour+1,:,:) = HRRR{'v-component_of_wind_height_above_ground'}(:,:,:); %V Winds @ 10m [m/s]
|
|
|
|
HRRR = ncgeodataset(['D:\DUNEX_RT\Operation\HRRR\hrrrT' num2str(ForcastHour,'%02d') '.grib2']);
|
|
pressure(ForcastHour+1,:,:) = HRRR{'Pressure_surface'}(:,:,:); %Surface Pressure [pa]
|
|
|
|
if ForcastHour==0
|
|
rain(ForcastHour+1,:,:) = HRRR{'Total_precipitation_surface_0_Hour_Accumulation'}(1,:,:); %Surface precipitation [kg/m^2]
|
|
elseif ForcastHour==1
|
|
rain(ForcastHour+1,:,:) = HRRR{'Total_precipitation_surface_1_Hour_Accumulation'}(1,:,:); %Surface precipitation [kg/m^2]
|
|
else
|
|
rain(ForcastHour+1,:,:) = HRRR{'Total_precipitation_surface_Mixed_intervals_Accumulation'}(2,:,:); %Surface precipitation [kg/m^2]
|
|
end
|
|
temperature(ForcastHour+1) = squeeze(nanmean(nanmean(HRRR{'Temperature_surface'}(:,:,:)- 273.15,2),3)); %Surface Temperature [C]
|
|
end
|
|
% Extract Temp
|
|
|
|
disp('Loaded HRRR')
|
|
|
|
|
|
%% IMPORT NCOM
|
|
|
|
try
|
|
nComDate=masterDate;
|
|
% ncomeURL = ['http://nomads.ncep.noaa.gov:80/dods/ncom/ncom' datestr(nComDate,'YYYYmmDD') '/ncom_useast_'];
|
|
% ncomeURL = ['https://www.ncei.noaa.gov/thredds-coastal/dodsC/us_east/us_east_20201218_to_current/20201220/coamps_ncom_useast_u_1_2020122000_00000000.nc'];
|
|
% ncomeURL = ['https://www.ncei.noaa.gov/thredds-coastal/dodsC/us_east/us_east_20201218_to_current/' datestr(nComDate,'YYYYmmDD')...
|
|
% '/coamps_ncom_useast_u_1_'];
|
|
|
|
% mget(ncom_ftp,['grids/operational/NCOM/regional/GRIB2/ncom_useast_' datestr(nComDate,'YYYYmmDD') '00.grb2'],['D:\DUNEX_RT\Operation\NCOM\']);
|
|
ncomURL = ['https://nomads.ncep.noaa.gov/pub/data/nccf/com/ncom/prod/ncom.' datestr(nComDate,'YYYYmmDD') '/' 'useast_u_ocn_ncout_grid1_'...
|
|
datestr(nComDate,'YYYYmmDD') '00_t0025-0048.tgz'];
|
|
|
|
ncom_downloadPath = ['D:\DUNEX_RT\Operation\NCOM\' datestr(masterDate,'YYYYmmDD') '\ncom_useast_' datestr(nComDate,'YYYYmmDD') '0049-0072.tgz'];
|
|
if exist(ncom_downloadPath) ~= 2
|
|
mkdir(['D:\DUNEX_RT\Operation\NCOM\' datestr(masterDate,'YYYYmmDD')])
|
|
|
|
ncom_downloadPath = ['D:\DUNEX_RT\Operation\NCOM\' datestr(masterDate,'YYYYmmDD') '\ncom_useast_' datestr(nComDate,'YYYYmmDD') '0025-0048.tgz'];
|
|
websave(ncom_downloadPath,ncomURL,options)
|
|
untar(ncom_downloadPath,['D:\DUNEX_RT\Operation\NCOM\' datestr(masterDate,'YYYYmmDD')])
|
|
|
|
ncomURL = ['https://nomads.ncep.noaa.gov/pub/data/nccf/com/ncom/prod/ncom.' datestr(nComDate,'YYYYmmDD') '/' 'useast_u_ocn_ncout_grid1_'...
|
|
datestr(nComDate,'YYYYmmDD') '00_t0000-0024.tgz'];
|
|
ncom_downloadPath = ['D:\DUNEX_RT\Operation\NCOM\' datestr(masterDate,'YYYYmmDD') '\ncom_useast_' datestr(nComDate,'YYYYmmDD') '0000-0024.tgz'];
|
|
websave(ncom_downloadPath,ncomURL,options)
|
|
untar(ncom_downloadPath,['D:\DUNEX_RT\Operation\NCOM\' datestr(masterDate,'YYYYmmDD')])
|
|
|
|
ncomURL = ['https://nomads.ncep.noaa.gov/pub/data/nccf/com/ncom/prod/ncom.' datestr(nComDate,'YYYYmmDD') '/' 'useast_u_ocn_ncout_grid1_'...
|
|
datestr(nComDate,'YYYYmmDD') '00_t0049-0072.tgz'];
|
|
ncom_downloadPath = ['D:\DUNEX_RT\Operation\NCOM\' datestr(masterDate,'YYYYmmDD') '\ncom_useast_' datestr(nComDate,'YYYYmmDD') '0049-0072.tgz'];
|
|
websave(ncom_downloadPath,ncomURL,options)
|
|
untar(ncom_downloadPath,['D:\DUNEX_RT\Operation\NCOM\' datestr(masterDate,'YYYYmmDD')])
|
|
|
|
ncomURL = ['https://nomads.ncep.noaa.gov/pub/data/nccf/com/ncom/prod/ncom.' datestr(nComDate,'YYYYmmDD') '/' 'useast_u_ocn_ncout_grid1_'...
|
|
datestr(nComDate,'YYYYmmDD') '00_t0073-0096.tgz'];
|
|
ncom_downloadPath = ['D:\DUNEX_RT\Operation\NCOM\' datestr(masterDate,'YYYYmmDD') '\ncom_useast_' datestr(nComDate,'YYYYmmDD') '0073-0096.tgz'];
|
|
websave(ncom_downloadPath,ncomURL,options)
|
|
untar(ncom_downloadPath,['D:\DUNEX_RT\Operation\NCOM\' datestr(masterDate,'YYYYmmDD')])
|
|
end
|
|
|
|
if exist(['D:\DUNEX_RT\Operation\NCOM\' datestr(masterDate-days(2),'YYYYmmDD')]) == 7
|
|
rmdir(['D:\DUNEX_RT\Operation\NCOM\' datestr(masterDate-days(2),'YYYYmmDD')], 's')
|
|
end
|
|
ncom_filePath = ['D:\DUNEX_RT\Operation\NCOM\' datestr(masterDate,'YYYYmmDD') '\coamps_ncom_useast_u_1_' datestr(nComDate,'YYYYmmDD') '00_00'];
|
|
% NCOM.Lat = ncread([ncomeURL datestr(nComDate,'YYYYmmDD')],'lat',[331],[240]);
|
|
% NCOM.Lat = ncread([ncomeURL datestr(nComDate,'YYYYmmDD') '00_00000000.nc'],'lat',[331],[240]);
|
|
NCOM.Lat = ncread([ncom_filePath '48' '0000.nc'],'lat',[331],[240]);
|
|
NCOM.Lon = ncread([ncom_filePath '48' '0000.nc'],'lon',91,180);
|
|
% NCOM.Lat = ncom_grib{'lat',[331],[240]};
|
|
% NCOM.Lon = ncom_grib{'lon',91,180};
|
|
catch
|
|
try
|
|
if hour(masterTime) ~= 18
|
|
nComDate=masterDate-days(1);
|
|
else
|
|
nComDate=masterDate;
|
|
end
|
|
% ncomeURL = ['http://nomads.ncep.noaa.gov:80/dods/ncom/ncom' datestr(nComDate,'YYYYmmDD') '/ncom_useast_'];
|
|
% ncomeURL = ['https://www.ncei.noaa.gov/thredds-coastal/dodsC/us_east/us_east_20201218_to_current/' datestr(nComDate,'YYYYmmDD')...
|
|
% '/coamps_ncom_useast_u_1_'];
|
|
% mget(ncom_ftp,['grids/operational/NCOM/regional/GRIB2/ncom_useast_' datestr(nComDate,'YYYYmmDD') '00.grb2'],['D:\DUNEX_RT\Operation\NCOM\']);
|
|
% ncom_downloadPath = ['D:\DUNEX_RT\Operation\NCOM\grids\operational\NCOM\regional\GRIB2\ncom_useast_' datestr(nComDate,'YYYYmmDD') '00.grb2'];
|
|
% ncom_grib = ncgeodataset(ncom_downloadPath);
|
|
ncom_downloadPath = ['D:\DUNEX_RT\Operation\NCOM\' datestr(masterDate,'YYYYmmDD') '\ncom_useast_' datestr(nComDate,'YYYYmmDD') '0049-0072.tgz'];
|
|
if exist(ncom_downloadPath) ~= 2
|
|
ncomURL = ['https://nomads.ncep.noaa.gov/pub/data/nccf/com/ncom/prod/ncom.' datestr(nComDate,'YYYYmmDD') '/' 'useast_u_ocn_ncout_grid1_'...
|
|
datestr(nComDate,'YYYYmmDD') '00_t0025-0048.tgz'];
|
|
ncom_downloadPath = ['D:\DUNEX_RT\Operation\NCOM\' datestr(masterDate,'YYYYmmDD') '\ncom_useast_' datestr(nComDate,'YYYYmmDD') '00.tgz'];
|
|
websave(ncom_downloadPath,ncomURL,options)
|
|
untar(ncom_downloadPath,['D:\DUNEX_RT\Operation\NCOM\' datestr(masterDate,'YYYYmmDD')])
|
|
|
|
ncomURL = ['https://nomads.ncep.noaa.gov/pub/data/nccf/com/ncom/prod/ncom.' datestr(nComDate,'YYYYmmDD') '/' 'useast_u_ocn_ncout_grid1_'...
|
|
datestr(nComDate,'YYYYmmDD') '00_t0000-0024.tgz'];
|
|
ncom_downloadPath = ['D:\DUNEX_RT\Operation\NCOM\' datestr(masterDate,'YYYYmmDD') '\ncom_useast_' datestr(nComDate,'YYYYmmDD') '0000-0024.tgz'];
|
|
websave(ncom_downloadPath,ncomURL,options)
|
|
untar(ncom_downloadPath,['D:\DUNEX_RT\Operation\NCOM\' datestr(masterDate,'YYYYmmDD')])
|
|
|
|
ncomURL = ['https://nomads.ncep.noaa.gov/pub/data/nccf/com/ncom/prod/ncom.' datestr(nComDate,'YYYYmmDD') '/' 'useast_u_ocn_ncout_grid1_'...
|
|
datestr(nComDate,'YYYYmmDD') '00_t0049-0072.tgz'];
|
|
ncom_downloadPath = ['D:\DUNEX_RT\Operation\NCOM\' datestr(masterDate,'YYYYmmDD') '\ncom_useast_' datestr(nComDate,'YYYYmmDD') '0049-0072.tgz'];
|
|
websave(ncom_downloadPath,ncomURL,options)
|
|
untar(ncom_downloadPath,['D:\DUNEX_RT\Operation\NCOM\' datestr(masterDate,'YYYYmmDD')])
|
|
|
|
ncomURL = ['https://nomads.ncep.noaa.gov/pub/data/nccf/com/ncom/prod/ncom.' datestr(nComDate,'YYYYmmDD') '/' 'useast_u_ocn_ncout_grid1_'...
|
|
datestr(nComDate,'YYYYmmDD') '00_t0073-0096.tgz'];
|
|
ncom_downloadPath = ['D:\DUNEX_RT\Operation\NCOM\' datestr(masterDate,'YYYYmmDD') '\ncom_useast_' datestr(nComDate,'YYYYmmDD') '0073-0096.tgz'];
|
|
websave(ncom_downloadPath,ncomURL,options)
|
|
untar(ncom_downloadPath,['D:\DUNEX_RT\Operation\NCOM\' datestr(masterDate,'YYYYmmDD')])
|
|
|
|
end
|
|
|
|
if exist(['D:\DUNEX_RT\Operation\NCOM\' datestr(masterDate-days(2),'YYYYmmDD')]) == 7
|
|
rmdir(['D:\DUNEX_RT\Operation\NCOM\' datestr(masterDate-days(2),'YYYYmmDD')], 's')
|
|
end
|
|
ncom_filePath = ['D:\DUNEX_RT\Operation\NCOM\' datestr(masterDate,'YYYYmmDD') '\coamps_ncom_useast_u_1_' datestr(nComDate,'YYYYmmDD') '00_00'];
|
|
% NCOM.Lat = ncread([ncomeURL datestr(nComDate,'YYYYmmDD')],'lat',[331],[240]);
|
|
% NCOM.Lat = ncread([ncomeURL datestr(nComDate,'YYYYmmDD') '00_00000000.nc'],'lat',[331],[240]);
|
|
NCOM.Lat = ncread([ncom_filePath '48' '0000.nc'],'lat',[331],[240]);
|
|
NCOM.Lon = ncread([ncom_filePath '48' '0000.nc'],'lon',91,180);
|
|
% NCOM.Lat = ncom_grib{'lat',[331],[240]};
|
|
% NCOM.Lon = ncom_grib{'lon',91,180};
|
|
|
|
catch
|
|
load(['ModelOut_03-Mar-2020_18Z'],'NCOM');
|
|
fid = fopen('errorFile.log','a+');
|
|
% write the error to file
|
|
% first line: message
|
|
fprintf(fid,'%s : %s\n',datestr(masterTime),'NCOM Fallback');
|
|
sendmail('alexander.rey@queensu.ca',[datestr(masterTime) ': NCOM Fallback']) ;
|
|
|
|
% close file
|
|
fclose(fid)
|
|
end
|
|
end
|
|
|
|
if FIRST_RUN == 1
|
|
NCOM.depth = ncread('D:\Alexander\NCOM\2019090500_t000.nc','depth');
|
|
end
|
|
|
|
Noffset = hours(masterTime - nComDate);
|
|
tStep = Noffset;
|
|
tCount = 1;
|
|
|
|
try
|
|
for Ni = masterTime:hours(3):masterTime+hours(48)
|
|
noData=1;
|
|
while noData==1
|
|
try
|
|
% NCOMstep = tStep./3+1;
|
|
NCOMstep = 1;
|
|
|
|
% NCOM.water_u(:,:,:,tCount) = ncread([ncomeURL datestr(masterDate,'YYYYmmDD')],...
|
|
% 'water_u',[NCOMstep 1 331 91],[1 40 240 180],[1 1 1 1]);
|
|
% NCOM.water_u(:,:,:,tCount) = ncread([ncomeURL datestr(nComDate,'YYYYmmDD') '00_00' num2str(tStep) '0000.nc'],...
|
|
% 'water_u',[91 331 1 NCOMstep],[180 240 inf 1],[1 1 1 1]);
|
|
% NCOM.water_v(:,:,:,tCount) = ncread([ncomeURL datestr(nComDate,'YYYYmmDD') '00_00' num2str(tStep) '0000.nc'],...
|
|
% 'water_v',[91 331 1 NCOMstep],[180 240 inf 1],[1 1 1 1]);
|
|
% NCOM.surf_el(:,:,tCount) = ncread([ncomeURL datestr(nComDate,'YYYYmmDD') '00_00' num2str(tStep) '0000.nc'],...
|
|
% 'surf_el',[91 331 NCOMstep],[180 240 1],[1 1 1]);
|
|
% NCOM.salinity(:,:,:,tCount) = ncread([ncomeURL datestr(nComDate,'YYYYmmDD') '00_00' num2str(tStep) '0000.nc'],...
|
|
% 'salinity',[91 331 1 NCOMstep],[180 240 inf 1],[1 1 1 1]);
|
|
NCOM.water_u(:,:,:,tCount) = ncread([ncom_filePath num2str(tStep,'%02d') '0000.nc'],...
|
|
'water_u',[91 331 1 1],[180 240 inf 1],[1 1 1 1]);
|
|
NCOM.water_v(:,:,:,tCount) = ncread([ncom_filePath num2str(tStep,'%02d') '0000.nc'],...
|
|
'water_v',[91 331 1 1],[180 240 inf 1],[1 1 1 1]);
|
|
NCOM.surf_el(:,:,tCount) = ncread([ncom_filePath num2str(tStep,'%02d') '0000.nc'],...
|
|
'surf_el',[91 331 1],[180 240 1],[1 1 1]);
|
|
NCOM.salinity(:,:,:,tCount) =ncread([ncom_filePath num2str(tStep,'%02d') '0000.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
|
|
catch
|
|
noData=1;
|
|
end
|
|
disp(Ni)
|
|
end
|
|
NCOM.time(tCount) = Ni;
|
|
tStep=tStep+3;
|
|
tCount = tCount+1;
|
|
end
|
|
catch
|
|
load(['ModelOut_03-Mar-2020_18Z'],'NCOM');
|
|
fid = fopen('errorFile.log','a+');
|
|
% write the error to file
|
|
% first line: message
|
|
fprintf(fid,'%s : %s\n',datestr(masterTime),'NCOM Fallback');
|
|
sendmail('alexander.rey@queensu.ca',[datestr(masterTime) ': NCOM Fallback']) ;
|
|
|
|
% close file
|
|
fclose(fid)
|
|
end
|
|
|
|
|
|
|
|
%% NCOM Process
|
|
NCOM.mag = vecmag(NCOM.water_u,NCOM.water_v);
|
|
|
|
if FIRST_RUN == 1
|
|
depthWeight = [diff(NCOM.depth)' 0];
|
|
|
|
for Layer = 1:40
|
|
DAVfac(:,:,Layer,:) = ones(180,240,1,17).*depthWeight(Layer);
|
|
end
|
|
DAVtrue = DAVfac;
|
|
N_depth = NCOM.depth;
|
|
save('NCOMstuff.mat','DAVtrue','DAVfac','N_depth')
|
|
else
|
|
load('NCOMstuff.mat')
|
|
NCOM.depth = N_depth;
|
|
clear N_depth;
|
|
end
|
|
DAVtrue(isnan(NCOM.mag(:,:,1:40,:))==1) = 0;
|
|
DAVfac(isnan(NCOM.mag(:,:,1:40,:))==1) = nan;
|
|
DAVFac2 = sum(DAVtrue,3);
|
|
clear DAVtrue
|
|
NCOM.dav = nansum(DAVfac.*NCOM.mag(:,:,1:40,:),3)./DAVFac2;
|
|
NCOM.davX = nansum(DAVfac.*NCOM.water_u(:,:,1:40,:),3)./DAVFac2;
|
|
NCOM.davY = nansum(DAVfac.*NCOM.water_v(:,:,1:40,:),3)./DAVFac2;
|
|
|
|
Import=1
|
|
catch err
|
|
|
|
clear waterLevels
|
|
|
|
%open file
|
|
fid = fopen('errorFile.log','a+');
|
|
% write the error to file
|
|
% first line: message
|
|
fprintf(fid,'Line %d on %s : %s\n',err.stack.line,datestr(masterTime),err.message);
|
|
sendmail('alexander.rey@queensu.ca',['DUNEX-RT: ' datestr(masterTime) err.message]) ;
|
|
|
|
% close file
|
|
fclose(fid)
|
|
disp('FAILED TO IMPORT FROM NOMADS')
|
|
pause(600)
|
|
end
|
|
end
|
|
|
|
|
|
load('shoreline.mat');
|
|
|
|
%% Format data for Delft3D
|
|
%% Evap
|
|
fileID=fopen('TemperatureRT.eva','w');
|
|
refTime=minutes(masterTime-masterDate);
|
|
records=49;
|
|
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);
|
|
|
|
%% NCOM ESTOFS
|
|
% Currents from NCOM
|
|
|
|
idxCount=1;
|
|
for i = 2:10: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:11: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:5: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
|
|
|
|
%% 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
|
|
|
|
vDatOrder = [crBoundCount(3)+1:1:crBoundCount(4) ...
|
|
crBoundCount(1)+1:1:crBoundCount(2) ...
|
|
crBoundCount(3):-1:crBoundCount(2)+1];
|
|
|
|
crBoundCountOld = [0 27 34 46];
|
|
vDatOrderOld = [crBoundCountOld(3)+1:1:crBoundCountOld(4) ...
|
|
crBoundCountOld(1)+1:1:crBoundCountOld(2) ...
|
|
crBoundCountOld(3):-1:crBoundCountOld(2)+1];
|
|
|
|
vDatNavD = LMSL_RT(:,3);
|
|
vDatNavD(vDatNavD==-999999)=nan;
|
|
|
|
% interp1(vDatNavD(vDatOrder),1:46)
|
|
vDatNavDInter(vDatOrder) = naninterp(interp1(...
|
|
linspace(1,length(vDatOrder),length(vDatOrderOld)),...
|
|
vDatNavD(vDatOrderOld),1:length(vDatOrder)));
|
|
for i = 1:length(idxNCOM)
|
|
wlO(:,i) = interp1(masterTime:hours(1):masterTime+hours(48),waterLevels(:,i)+0,masterTime:minutes(10):masterTime+hours(48),'pchip');
|
|
end
|
|
|
|
%% Match to grid
|
|
flowBed = load('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),:))+vDatNavDInter(i);
|
|
mag(:,idxCount) = vecmag(currX,currY);
|
|
|
|
% if idxCount==1
|
|
% DM(idxCount) = -squeeze(flowBed.data.Val(crIDX(idxCount+1,1)-1,crIDX(idxCount+1,2)));
|
|
% else
|
|
DM(idxCount) = -squeeze(flowBed.data.Val(crIDX(idxCount,1)-1,crIDX(idxCount,2)));
|
|
% end
|
|
|
|
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);
|
|
% currents(:,idxCount) = sind(dirGrid(i)-0)'.*vecmag(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:minutes(10):masterTime+hours(48)));
|
|
|
|
|
|
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) = 1.0.*currentsOut(:,idxCount) - wlO(:,idxCount).*sqrt(9.81./DM(idxCount));
|
|
% RNOut_2(:,i) = 1.0.*currentsOut(:,idxCount) - 0.*sqrt(9.81./DM(idxCount));
|
|
|
|
elseif i>crBoundCount(3)
|
|
if i<230
|
|
RNOut_2(:,i) = 1.0.*currentsOut(:,idxCount) + wlO(:,idxCount).*sqrt(9.81./DM(idxCount));
|
|
% RNOut_2(:,i) = 1.0.*currentsOut(:,idxCount) + 1.*sqrt(9.81./DM(idxCount));
|
|
else
|
|
RNOut_2(:,i) = 0.0.*currentsOut(:,idxCount) + wlO(:,idxCount).*sqrt(9.81./DM(idxCount));
|
|
% RNOut_2(:,i) = 1.0.*currentsOut(:,idxCount) + 0.*sqrt(9.81./DM(idxCount));
|
|
end
|
|
else
|
|
RNOut_2(:,i) = 0.0.*currentsOut(:,idxCount) - wlO(:,idxCount).*sqrt(9.81./DM(idxCount));
|
|
% RNOut_2(:,i) = -1.0.*currentsOut(:,idxCount) - 0.*sqrt(9.81./DM(idxCount));
|
|
end
|
|
|
|
|
|
|
|
idxCount = idxCount+1;
|
|
end
|
|
|
|
%% R Boundry Files
|
|
records=size(RNOut_2,1);
|
|
|
|
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:records
|
|
fprintf(fileID,'%d %1.4f %1.4f\n',refTime,RNOut_2(line,Boundry-1),RNOut_2(line,Boundry));
|
|
|
|
refTime=refTime+10;
|
|
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+2880);
|
|
|
|
|
|
|
|
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 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,['Commnt = initial conditions from map file\n'...
|
|
% 'Restid = #Jan16_06#\n']);
|
|
|
|
% fprintf( fileID,'Zeta0 = %#.7e \n',nanmean(nanmean(squeeze(wlO(1,:)))));
|
|
fprintf( fileID,'Zeta0 = %#.7e \n',0);
|
|
|
|
end
|
|
|
|
for i = 31:53
|
|
fprintf( fileID, '%s\n', cac{1}{i} );
|
|
end
|
|
|
|
% if FIRST_RUN==1
|
|
% fprintf( fileID, 'Vicouv = 1.0000000e+002 \n');
|
|
% else
|
|
fprintf( fileID, '%s\n', cac{1}{54} );
|
|
% end
|
|
|
|
for i = 55:64
|
|
fprintf( fileID, '%s\n', cac{1}{i} );
|
|
end
|
|
|
|
if FIRST_RUN==1
|
|
fprintf( fileID, 'Tlfsmo = 0.0000000e+001 \n');
|
|
else
|
|
fprintf( fileID, 'Tlfsmo = 0.0000000e+000 \n');
|
|
end
|
|
|
|
for i = 66: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+2880,refTime,refTime+2880,refTime,refTime+2880);
|
|
|
|
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';
|
|
bounDef{1}='grid-coordinates';
|
|
bounDef{2}='orientation';
|
|
bounDef{3}='orientation';
|
|
|
|
|
|
for Boundry = 1 : 3
|
|
|
|
if Boundry ==1
|
|
fprintf(fileID,['[Boundary]\n'...
|
|
'Name = East%d\n'...
|
|
'Definition = %s\n'...
|
|
'StartCoordM = %d\n'...
|
|
'EndCoordM = %d\n'...
|
|
'StartCoordN = %d\n'...
|
|
'EndCoordN = %d\n'...
|
|
'SpectrumSpec = parametric\n'...
|
|
'SpShapeType = pierson-moskowitz\n'...
|
|
'GaussSpread = 9.9999998e-003\n'...
|
|
'PeriodType = peak\n'...
|
|
'DirSpreadType = power\n'...
|
|
'PeakEnhanceFac = 3.3000000e+000\n'...
|
|
],Boundry,bounDef{Boundry},432,432,2,528);
|
|
else
|
|
fprintf(fileID,['[Boundary]\n'...
|
|
'Name = East%d\n'...
|
|
'Definition = %s\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,bounDef{Boundry},bounOri{Boundry},bounRot{Boundry});
|
|
end
|
|
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');
|
|
G7bath = load('C7_1400_V2_D.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
|
|
landmask(G7bath.data.Val'<0)=1;
|
|
landmask=logical(landmask);
|
|
|
|
if masterTime>=datetime(2019,09,05,00,00,00)
|
|
runSteps = 49;
|
|
else
|
|
runSteps = 8;
|
|
end
|
|
|
|
% Zero out met array
|
|
HRRR_INTER = zeros(4,runSteps,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:runSteps
|
|
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',runSteps)
|
|
|
|
save(['ModelIn_6_HRRR' datestr(masterTime,'YYYY-mm-DD_HH') '.mat'],'-v7.3','windU','windV','pressure','rain','temperature','waves','waterLevels','NCOM')
|
|
|
|
%% Run Delft3D FLOW/WAVE!
|
|
% system('. /home/cluster/D7545/enviorment.sh'
|
|
|
|
system('runcoast_flow2d3d_parallel_wave.bat')
|
|
if exist('wavm-rt_run.dat')~=2
|
|
fid = fopen('errorFile.log','a+');
|
|
fprintf(fid,'Delft Batch File Fallback\n');
|
|
|
|
% close file
|
|
fclose(fid)
|
|
disp('Run Fallback')
|
|
pause(60)
|
|
|
|
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*
|
|
delete swn-diag.*
|
|
|
|
system('runcoast_flow2d3d_parallel_wave.bat')
|
|
end
|
|
|
|
%% Move things
|
|
pause(60)
|
|
|
|
eval(['copyfile trih-rt_run.dat D:\DUNEX_RT\Archive\' datestr(masterTime,'YYYY-mm-DD_HH') '_histOut7.dat'])
|
|
eval(['copyfile trih-rt_run.def D:\DUNEX_RT\Archive\' datestr(masterTime,'YYYY-mm-DD_HH') '_histOut7.def'])
|
|
|
|
eval(['movefile ModelIn_6_HRRR' datestr(masterTime,'YYYY-mm-DD_HH') '.mat D:\DUNEX_RT\Archive\ModelIn_7' datestr(masterTime,'YYYY-mm-DD_HH') '.mat'])
|
|
% eval(['copyfile tri-rst.rt* D:\Alexander\RealtimeHRRR\'])
|
|
|
|
|
|
eval(['copyfile hrrr.amu D:\DUNEX_RT\Archive\' datestr(masterTime,'YYYY-mm-DD_HH') '_hrrr7.amu'])
|
|
eval(['copyfile hrrr.amv D:\DUNEX_RT\Archive\' datestr(masterTime,'YYYY-mm-DD_HH') '_hrrr7.amv'])
|
|
eval(['copyfile hrrr.ampr D:\DUNEX_RT\Archive\' datestr(masterTime,'YYYY-mm-DD_HH') '_hrrr7.ampr'])
|
|
eval(['copyfile hrrr.amp D:\DUNEX_RT\Archive\' datestr(masterTime,'YYYY-mm-DD_HH') '_hrrr7.amp'])
|
|
eval(['copyfile trim-rt_run.dat D:\DUNEX_RT\Archive\' datestr(masterTime,'YYYY-mm-DD_HH') '_mapOut7.dat'])
|
|
eval(['copyfile trim-rt_run.def D:\DUNEX_RT\Archive\' datestr(masterTime,'YYYY-mm-DD_HH') '_mapOut7.def'])
|
|
eval(['copyfile wavm-rt_run.dat D:\DUNEX_RT\Archive\' datestr(masterTime,'YYYY-mm-DD_HH') '_waveOut7.dat'])
|
|
eval(['copyfile wavm-rt_run.def D:\DUNEX_RT\Archive\' datestr(masterTime,'YYYY-mm-DD_HH') '_waveOut7.def'])
|
|
|
|
|
|
|
|
%% Clear things not needed for plotting
|
|
clear HRRR_INTER flowGrid lonHRRR_Mesh windU windV rain pressure waveGrid
|
|
|
|
|
|
delete(gcp('nocreate'))
|
|
parpool(3);
|
|
%% Import Model Results
|
|
addpath(genpath('matlab/applications/delft3d_matlab'),'-begin')
|
|
|
|
modelWater = qpread(qpfopen('trim-rt_run.dat'),1,'water level','griddata',0,0,0);
|
|
flowGrid = load('Cone7W_100.mat');
|
|
waveGrid = load('Cone7W_250.mat');
|
|
flowBed = load('Dec6_100Bed.mat');
|
|
waveBed = load('Dec6_250Bed.mat');
|
|
|
|
%% Create GeoTIFF for water levels
|
|
clear X Y Val geoDat geoMeshX geoMeshY R xWorldLimits yWorldLimits rasterSize
|
|
|
|
% timeIDX = [1 4 7 13 25 37];
|
|
timeIDX = [1 7 13 25 37 49];
|
|
|
|
yMinOffset = 75;
|
|
yMaxOffset = 25;
|
|
|
|
parfor Time = 1:6
|
|
% X{Time} = modelWater.X;
|
|
% Y{Time} = modelWater.Y;
|
|
X{Time} = flowGrid.data.X;
|
|
Y{Time} = flowGrid.data.Y;
|
|
|
|
modelPlot_1 = squeeze(modelWater.Val(timeIDX(Time),1:end-1,1:end-1));
|
|
|
|
xWorldLimits{Time} = [min(min(X{Time}(1:end,1:end))) max(max(X{Time}(1:end,1:end)))];
|
|
yWorldLimits{Time} = [min(min(Y{Time}(1:end,1:end))) max(max(Y{Time}(1:end,1:end)))];
|
|
% rasterSize = [(size(data.X,1)-2).*2 (size(data.X,2)-2).*2];
|
|
rasterSize{Time} = [length(xWorldLimits{Time}(1):100:xWorldLimits{Time}(2)) length(yWorldLimits{Time}(1):100:yWorldLimits{Time}(2))];
|
|
|
|
R{Time} = maprefpostings(xWorldLimits{Time},yWorldLimits{Time},[rasterSize{Time}(2) rasterSize{Time}(1)])
|
|
[geoMeshX{Time},geoMeshY{Time}] = meshgrid(R{Time}.XWorldLimits(1):100:R{Time}.XWorldLimits(2),...
|
|
R{Time}.YWorldLimits(1):100:R{Time}.YWorldLimits(2));
|
|
|
|
modelMask = squeeze(flowBed.data.Val>1);
|
|
modelPlot_1(modelMask)=nan;
|
|
modelPlot_1(400:520,525:620)=nan;
|
|
modelPlot_1(:,1:yMinOffset)=nan;
|
|
modelPlot_1(:,size(modelMask,2)-yMaxOffset:size(modelMask,2))=nan;
|
|
|
|
geoDat{Time} = -griddata(reshape(X{Time}(1:end-1,1:end-1),1,[]),...
|
|
reshape(Y{Time}(1:end-1,1:end-1),1,[]),...
|
|
reshape(modelPlot_1,1,[]),...
|
|
geoMeshX{Time},geoMeshY{Time},'natural');
|
|
|
|
geoDat{Time}(isnan(geoDat{Time})==1) = -999;
|
|
|
|
geotiffwrite(['WL_' num2str(timeIDX(Time)) '.tif'],geoDat{Time},...
|
|
R{Time},'CoordRefSysCode',32618)
|
|
disp(Time)
|
|
% clear X Y Val geoDat geoMeshX geoMeshY
|
|
end
|
|
|
|
clear modelWater
|
|
|
|
modelWaterDAV = qpread(qpfopen('trim-rt_run.dat'),1,'depth averaged velocity','griddata',0,0,0);
|
|
|
|
%% Create GeoTIFF for DAVs
|
|
clear X Y Val geoDat geoMeshX geoMeshY R xWorldLimits yWorldLimits rasterSize geoDat modelPlot_1
|
|
|
|
% timeIDX = [1 4 7 13 25 37];
|
|
timeIDX = [1 7 13 25 37 49];
|
|
if FIRST_RUN == 1
|
|
Times = 2:6;
|
|
else
|
|
Times = 1:6;
|
|
end
|
|
|
|
|
|
parfor Time = Times
|
|
X{Time} = flowGrid.data.X;
|
|
Y{Time} = flowGrid.data.Y;
|
|
modelPlot_1 = sqrt(squeeze(modelWaterDAV.XComp(timeIDX(Time),1:end-1,1:end-1)).^2+...
|
|
squeeze(modelWaterDAV.YComp(timeIDX(Time),1:end-1,1:end-1)).^2);
|
|
|
|
xWorldLimits{Time} = [min(min(X{Time}(1:end,1:end))) max(max(X{Time}(1:end,1:end)))];
|
|
yWorldLimits{Time} = [min(min(Y{Time}(1:end,1:end))) max(max(Y{Time}(1:end,1:end)))];
|
|
% rasterSize = [(size(data.X,1)-2).*2 (size(data.X,2)-2).*2];
|
|
rasterSize{Time} = [length(xWorldLimits{Time}(1):100:xWorldLimits{Time}(2)) length(yWorldLimits{Time}(1):100:yWorldLimits{Time}(2))];
|
|
|
|
R{Time} = maprefpostings(xWorldLimits{Time},yWorldLimits{Time},[rasterSize{Time}(2) rasterSize{Time}(1)])
|
|
[geoMeshX{Time},geoMeshY{Time}] = meshgrid(R{Time}.XWorldLimits(1):100:R{Time}.XWorldLimits(2),...
|
|
R{Time}.YWorldLimits(1):100:R{Time}.YWorldLimits(2));
|
|
|
|
modelPlot_1(modelPlot_1==0)=nan;
|
|
|
|
modelMask = squeeze(flowBed.data.Val>1);
|
|
modelPlot_1(modelMask)=nan;
|
|
modelPlot_1(400:520,525:620)=nan;
|
|
modelPlot_1(:,1:yMinOffset)=nan;
|
|
modelPlot_1(:,size(modelMask,2)-yMaxOffset:size(modelMask,2))=nan;
|
|
|
|
|
|
geoDat{Time} = -griddata(reshape(X{Time}(1:end-1,1:end-1),1,[]),...
|
|
reshape(Y{Time}(1:end-1,1:end-1),1,[]),...
|
|
reshape(modelPlot_1,1,[]),...
|
|
geoMeshX{Time},geoMeshY{Time},'natural');
|
|
geoDat{Time}(isnan(geoDat{Time})==1) = -999;
|
|
|
|
geotiffwrite(['DAV_' num2str(timeIDX(Time)) '.tif'],geoDat{Time},...
|
|
R{Time},'CoordRefSysCode',32618)
|
|
disp(Time)
|
|
% clear X Y Val geoDat geoMeshX geoMeshY
|
|
end
|
|
|
|
%% Quiver for DAV
|
|
clear X Y ValX ValY XXd YYd geoS XdatE YdatE Xdat Ydat Udat Vdat geoDat
|
|
vecScale(4) = 100000;
|
|
vecScale(3) = 50000;
|
|
vecScale(2) = 2500;
|
|
vecScale(1) = 1000;
|
|
|
|
for Time = Times
|
|
|
|
X{Time} = flowGrid.data.X(2:end,2:end);
|
|
Y{Time} = flowGrid.data.Y(2:end,2:end);
|
|
modelPlot_1 = squeeze(modelWaterDAV.XComp(timeIDX(Time),2:end,2:end));
|
|
modelPlot_2 = squeeze(modelWaterDAV.YComp(timeIDX(Time),2:end,2:end));
|
|
|
|
modelMask = squeeze(flowBed.data.Val>1);
|
|
modelPlot_1(modelMask)=nan;
|
|
modelPlot_1(400:520,525:620)=nan;
|
|
modelPlot_2(modelMask)=nan;
|
|
modelPlot_2(400:520,525:620)=nan;
|
|
|
|
modelPlot_1(:,1:yMinOffset)=nan;
|
|
modelPlot_1(:,size(modelMask,2)-yMaxOffset:size(modelMask,2))=nan;
|
|
modelPlot_2(:,1:yMinOffset)=nan;
|
|
modelPlot_2(:,size(modelMask,2)-yMaxOffset:size(modelMask,2))=nan;
|
|
|
|
|
|
qSpace = [4 16 50 100]
|
|
|
|
Count{Time}=1;
|
|
|
|
for Detail = 1:4
|
|
q{Time} = quiver(X{Time}(1:qSpace(Detail):end-0,1:qSpace(Detail):end-0),Y{Time}(1:qSpace(Detail):end-0,1:qSpace(Detail):end-0),...
|
|
squeeze(modelPlot_1(1:qSpace(Detail):end-0,1:qSpace(Detail):end-0)),squeeze(modelPlot_2(1:qSpace(Detail):end-0,1:qSpace(Detail):end-0)),'color','k');
|
|
|
|
Xdat{Time} = reshape(q{Time}.XData,1,[]);
|
|
Ydat{Time} = reshape(q{Time}.YData,1,[]);
|
|
|
|
Udat{Time} = reshape(q{Time}.UData,1,[]);
|
|
Vdat{Time} = reshape(q{Time}.VData,1,[]);
|
|
|
|
for i = 1:length(Xdat{Time})
|
|
if ((Udat{Time}(i))==0 && (Vdat{Time}(i))==0) ||...
|
|
(isnan(Xdat{Time}(i))==1 || isnan(Ydat{Time}(i))==1) ||...
|
|
((isnan(Udat{Time}(i)))==1 && isnan(Vdat{Time}(i))==1)
|
|
continue
|
|
end
|
|
|
|
geoS{Time}(Count{Time}).Geometry = 'Line';
|
|
|
|
geoS{Time}(Count{Time}).Detail = Detail;
|
|
XdatE{Time} = Xdat{Time}(i)+Udat{Time}(i).*vecScale(Detail);
|
|
YdatE{Time} = Ydat{Time}(i)+Vdat{Time}(i).*vecScale(Detail);
|
|
[XXd{Time}(1),YYd{Time}(1)] = utm2deg(Xdat{Time}(i),Ydat{Time}(i),'18 S');
|
|
[XXd{Time}(2),YYd{Time}(2)] = utm2deg(XdatE{Time}',YdatE{Time}','18 S');
|
|
geoS{Time}(Count{Time}).BoundingBox = [min(XXd{Time}) min(YYd{Time});...
|
|
max(XXd{Time}) max(YYd{Time})];
|
|
geoS{Time}(Count{Time}).Lat = [XXd{Time}(1) XXd{Time}(2)];
|
|
geoS{Time}(Count{Time}).Lon = [YYd{Time}(1) YYd{Time}(2)];
|
|
|
|
|
|
Count{Time}=Count{Time}+1;
|
|
end
|
|
close all
|
|
|
|
Xdat{Time}=[];
|
|
Ydat{Time}=[];
|
|
Udat{Time}=[];
|
|
Vdat{Time}=[];
|
|
q{Time}=[];
|
|
end
|
|
disp(Time)
|
|
|
|
shapewrite(geoS{Time}(1:end-1),['DAV_' num2str(timeIDX(Time)) '_D.shp'])
|
|
eval(['copyfile vector.prj DAV_' num2str(timeIDX(Time)) '_D.prj'])
|
|
|
|
end
|
|
|
|
clear modelWaterDAV
|
|
|
|
%% Create GeoTIFF for HS
|
|
|
|
modelWaves = qpread(qpfopen('wavm-rt_run.dat'),1,'hsig wave height','griddata',0,0,0);
|
|
modelWavesDir = qpread(qpfopen('wavm-rt_run.dat'),1,'hsig wave vector (mean direction)','griddata',0,0,0);
|
|
|
|
clear X Y Val geoDat geoMeshX geoMeshY R xWorldLimits yWorldLimits rasterSize GeoS
|
|
|
|
% timeIDX = [1 4 7 13 25 37];
|
|
timeIDX = [1 7 13 25 37 49];
|
|
parfor Time = 1:6
|
|
X{Time} = waveGrid.data.X(1:end-1,1:end-1);
|
|
Y{Time} = waveGrid.data.Y(1:end-1,1:end-1);
|
|
modelPlot_1 = (squeeze(modelWaves.Val(timeIDX(Time),:,:)));
|
|
|
|
xWorldLimits{Time} = [min(min(X{Time}(1:end,1:end))) max(max(X{Time}(1:end,1:end)))];
|
|
yWorldLimits{Time} = [min(min(Y{Time}(1:end,1:end))) max(max(Y{Time}(1:end,1:end)))];
|
|
% rasterSize = [(size(data.X,1)-2).*2 (size(data.X,2)-2).*2];
|
|
rasterSize{Time} = [length(xWorldLimits{Time}(1):250:xWorldLimits{Time}(2)) length(yWorldLimits{Time}(1):250:yWorldLimits{Time}(2))];
|
|
|
|
R{Time} = maprefpostings(xWorldLimits{Time},yWorldLimits{Time},[rasterSize{Time}(2) rasterSize{Time}(1)])
|
|
[geoMeshX{Time},geoMeshY{Time}] = meshgrid(R{Time}.XWorldLimits(1):250:R{Time}.XWorldLimits(2),...
|
|
R{Time}.YWorldLimits(1):250:R{Time}.YWorldLimits(2));
|
|
|
|
modelPlot_1(modelPlot_1==0)=nan;
|
|
|
|
modelMask = squeeze(waveBed.data.Val<-1);
|
|
modelPlot_1(modelMask)=nan;
|
|
modelPlot_1(155:205,210:265)=nan;
|
|
|
|
geoDat{Time} = -griddata(reshape(X{Time}(1:end-0,1:end-0),1,[]),...
|
|
reshape(Y{Time}(1:end-0,1:end-0),1,[]),...
|
|
reshape(modelPlot_1,1,[]),...
|
|
geoMeshX{Time},geoMeshY{Time},'natural');
|
|
geoDat{Time}(isnan(geoDat{Time})==1) = -999;
|
|
|
|
geotiffwrite(['HS_' num2str(timeIDX(Time)) '.tif'],geoDat{Time},...
|
|
R{Time},'CoordRefSysCode',32618)
|
|
disp(Time)
|
|
% clear X Y Val geoDat geoMeshX geoMeshY
|
|
end
|
|
|
|
|
|
|
|
%% Quiver for HS
|
|
clear X Y ValX ValY XXd YYd geoS XdatE YdatE Xdat Ydat Udat Vdat geoDat
|
|
|
|
vecScale(4) = 15000;
|
|
vecScale(3) = 5000;
|
|
vecScale(2) = 1000;
|
|
vecScale(1) = 500;
|
|
|
|
for Time = 1:6
|
|
|
|
X{Time} = waveGrid.data.X(1:end-1,1:end-1);
|
|
Y{Time} = waveGrid.data.Y(1:end-1,1:end-1);
|
|
modelPlot_1 = squeeze(modelWavesDir.XComp(timeIDX(Time),:,:));
|
|
modelPlot_2 = squeeze(modelWavesDir.YComp(timeIDX(Time),:,:));
|
|
|
|
modelMask = squeeze(waveBed.data.Val<-1);
|
|
modelPlot_1(modelMask)=nan;
|
|
modelPlot_1(155:205,210:265)=nan;
|
|
modelPlot_2(modelMask)=nan;
|
|
modelPlot_2(155:205,210:265)=nan;
|
|
|
|
|
|
qSpace = [2 10 20 40]
|
|
Count{Time}=1;
|
|
|
|
for Detail = 1:4
|
|
q{Time} = quiver(X{Time}(1:qSpace(Detail):end-0,1:qSpace(Detail):end-0),Y{Time}(1:qSpace(Detail):end-0,1:qSpace(Detail):end-0),...
|
|
squeeze(modelPlot_1(1:qSpace(Detail):end-0,1:qSpace(Detail):end-0)),squeeze(modelPlot_2(1:qSpace(Detail):end-0,1:qSpace(Detail):end-0)),'color','k');
|
|
|
|
Xdat{Time} = reshape(q{Time}.XData,1,[]);
|
|
Ydat{Time} = reshape(q{Time}.YData,1,[]);
|
|
|
|
Udat{Time} = reshape(q{Time}.UData,1,[]);
|
|
Vdat{Time} = reshape(q{Time}.VData,1,[]);
|
|
|
|
|
|
for i = 1:length(Xdat{Time})
|
|
if (Udat{Time}(i))==0 && (Vdat{Time}(i))==0
|
|
continue
|
|
elseif isnan(Xdat{Time}(i))==1 || isnan(Ydat{Time}(i))==1
|
|
continue
|
|
elseif isnan(Udat{Time}(i))==1 && isnan(Vdat{Time}(i))==1
|
|
continue
|
|
end
|
|
|
|
geoS{Time}(Count{Time}).Geometry = 'Line';
|
|
|
|
geoS{Time}(Count{Time}).Detail = Detail;
|
|
XdatE{Time} = Xdat{Time}(i)+Udat{Time}(i).*vecScale(Detail);
|
|
YdatE{Time} = Ydat{Time}(i)+Vdat{Time}(i).*vecScale(Detail);
|
|
[XXd{Time}(1),YYd{Time}(1)] = utm2deg(Xdat{Time}(i),Ydat{Time}(i),'18 S');
|
|
[XXd{Time}(2),YYd{Time}(2)] = utm2deg(XdatE{Time}',YdatE{Time}','18 S');
|
|
geoS{Time}(Count{Time}).BoundingBox = [min(XXd{Time}) min(YYd{Time});...
|
|
max(XXd{Time}) max(YYd{Time})];
|
|
geoS{Time}(Count{Time}).Lat = [XXd{Time}(1) XXd{Time}(2)];
|
|
geoS{Time}(Count{Time}).Lon = [YYd{Time}(1) YYd{Time}(2)];
|
|
|
|
|
|
Count{Time}=Count{Time}+1;
|
|
end
|
|
close all
|
|
end
|
|
disp(Time)
|
|
|
|
shapewrite(geoS{Time}(1:end-1),['HS_' num2str(timeIDX(Time)) '_D.shp'])
|
|
eval(['copyfile vector.prj HS_' num2str(timeIDX(Time)) '_D.prj'])
|
|
|
|
end
|
|
|
|
clear modelWaves
|
|
clear modelWavesDIR
|
|
|
|
%% Create GeoTIFF for Wind
|
|
modelWaterWind = qpread(qpfopen('trim-rt_run.dat'),1,'wind speed','griddata',0,0,0);
|
|
|
|
|
|
clear X Y Val geoDat geoMeshX geoMeshY R xWorldLimits yWorldLimits rasterSize geoS ValX ValY
|
|
|
|
% timeIDX = [1 4 7 13 25 37];
|
|
timeIDX = [1 7 13 25 37 49];
|
|
parfor Time = 1:6
|
|
X{Time} = flowGrid.data.X;
|
|
Y{Time} = flowGrid.data.Y;
|
|
Val{Time} = sqrt(squeeze(modelWaterWind.XComp(timeIDX(Time),:,:)).^2+...
|
|
squeeze(modelWaterWind.YComp(timeIDX(Time),:,:)).^2);
|
|
|
|
xWorldLimits{Time} = [min(min(X{Time}(1:end,1:end))) max(max(X{Time}(1:end,1:end)))];
|
|
yWorldLimits{Time} = [min(min(Y{Time}(1:end,1:end))) max(max(Y{Time}(1:end,1:end)))];
|
|
% rasterSize = [(size(data.X,1)-2).*2 (size(data.X,2)-2).*2];
|
|
rasterSize{Time} = [length(xWorldLimits{Time}(1):1000:xWorldLimits{Time}(2)) length(yWorldLimits{Time}(1):1000:yWorldLimits{Time}(2))];
|
|
|
|
R{Time} = maprefpostings(xWorldLimits{Time},yWorldLimits{Time},[rasterSize{Time}(2) rasterSize{Time}(1)])
|
|
[geoMeshX{Time},geoMeshY{Time}] = meshgrid(R{Time}.XWorldLimits(1):1000:R{Time}.XWorldLimits(2),...
|
|
R{Time}.YWorldLimits(1):1000:R{Time}.YWorldLimits(2));
|
|
|
|
Val{Time}(Val{Time}==0)=nan;
|
|
geoDat{Time} = -griddata(reshape(X{Time}(2:end-2,2:end-2),1,[]),...
|
|
reshape(Y{Time}(2:end-2,2:end-2),1,[]),...
|
|
reshape(Val{Time}(2:end-2,2:end-2),1,[]),...
|
|
geoMeshX{Time},geoMeshY{Time},'natural');
|
|
geoDat{Time}(isnan(geoDat{Time})==1) = -999;
|
|
|
|
geotiffwrite(['Uw_' num2str(timeIDX(Time)) '.tif'],geoDat{Time},...
|
|
R{Time},'CoordRefSysCode',32618)
|
|
disp(Time)
|
|
% clear X Y Val geoDat geoMeshX geoMeshY
|
|
end
|
|
|
|
|
|
|
|
%% Quiver for Wind
|
|
clear X Y ValX ValY XXd YYd geoS XdatE YdatE Xdat Ydat Udat Vdat geoDat
|
|
|
|
|
|
vecScale(4) = 3000;
|
|
vecScale(3) = 500;
|
|
vecScale(2) = 100;
|
|
vecScale(1) = 40;
|
|
|
|
for Time = 1:6
|
|
|
|
X{Time} = flowGrid.data.X;
|
|
Y{Time} = flowGrid.data.Y;
|
|
ValX{Time} = squeeze(modelWaterWind.XComp(timeIDX(Time),:,:));
|
|
ValY{Time} = squeeze(modelWaterWind.YComp(timeIDX(Time),:,:));
|
|
|
|
|
|
qSpace = [5 25 50 100]
|
|
Count{Time}=1;
|
|
|
|
for Detail = 1:4
|
|
q{Time} = quiver(X{Time}(2:qSpace(Detail):end-0,2:qSpace(Detail):end-0),Y{Time}(2:qSpace(Detail):end-0,2:qSpace(Detail):end-0),...
|
|
squeeze(ValX{Time}(2:qSpace(Detail):end-0,2:qSpace(Detail):end-0)),squeeze(ValY{Time}(2:qSpace(Detail):end-0,2:qSpace(Detail):end-0)),'color','k');
|
|
|
|
Xdat{Time} = reshape(q{Time}.XData,1,[]);
|
|
Ydat{Time} = reshape(q{Time}.YData,1,[]);
|
|
|
|
Udat{Time} = reshape(q{Time}.UData,1,[]);
|
|
Vdat{Time} = reshape(q{Time}.VData,1,[]);
|
|
|
|
|
|
for i = 1:length(Xdat{Time})
|
|
if (Udat{Time}(i))==0 && (Vdat{Time}(i))==0
|
|
continue
|
|
elseif isnan(Xdat{Time}(i))==1 || isnan(Ydat{Time}(i))==1
|
|
continue
|
|
end
|
|
|
|
geoS{Time}(Count{Time}).Geometry = 'Line';
|
|
|
|
geoS{Time}(Count{Time}).Detail = Detail;
|
|
XdatE{Time} = Xdat{Time}(i)+Udat{Time}(i).*vecScale(Detail);
|
|
YdatE{Time} = Ydat{Time}(i)+Vdat{Time}(i).*vecScale(Detail);
|
|
[XXd{Time}(1),YYd{Time}(1)] = utm2deg(Xdat{Time}(i),Ydat{Time}(i),'18 S');
|
|
[XXd{Time}(2),YYd{Time}(2)] = utm2deg(XdatE{Time}',YdatE{Time}','18 S');
|
|
geoS{Time}(Count{Time}).BoundingBox = [min(XXd{Time}) min(YYd{Time});...
|
|
max(XXd{Time}) max(YYd{Time})];
|
|
geoS{Time}(Count{Time}).Lat = [XXd{Time}(1) XXd{Time}(2)];
|
|
geoS{Time}(Count{Time}).Lon = [YYd{Time}(1) YYd{Time}(2)];
|
|
|
|
|
|
Count{Time}=Count{Time}+1;
|
|
end
|
|
close all
|
|
end
|
|
disp(Time)
|
|
|
|
shapewrite(geoS{Time}(1:end-1),['Uw_' num2str(timeIDX(Time)) '_D.shp'])
|
|
eval(['copyfile vector.prj Uw_' num2str(timeIDX(Time)) '_D.prj'])
|
|
|
|
end
|
|
|
|
clear modelWaterWind
|
|
|
|
%% Create GeoTIFF for Prep
|
|
|
|
modelWaterPrep = qpread(qpfopen('trim-rt_run.dat'),1,'precipitation rate','griddata',0,0,0);
|
|
|
|
clear X Y Val geoDat geoMeshX geoMeshY R xWorldLimits yWorldLimits rasterSize geoS
|
|
|
|
% timeIDX = [1 4 7 13 25 37];
|
|
timeIDX = [1 7 13 25 37 49];
|
|
parfor Time = 1:6
|
|
X{Time} = flowGrid.data.X;
|
|
Y{Time} = flowGrid.data.Y;
|
|
Val{Time} = (squeeze(modelWaterPrep.Val(timeIDX(Time),:,:)));
|
|
|
|
xWorldLimits{Time} = [min(min(X{Time}(1:end,1:end))) max(max(X{Time}(1:end,1:end)))];
|
|
yWorldLimits{Time} = [min(min(Y{Time}(1:end,1:end))) max(max(Y{Time}(1:end,1:end)))];
|
|
% rasterSize = [(size(data.X,1)-2).*2 (size(data.X,2)-2).*2];
|
|
rasterSize{Time} = [length(xWorldLimits{Time}(1):10000:xWorldLimits{Time}(2)) length(yWorldLimits{Time}(1):10000:yWorldLimits{Time}(2))];
|
|
|
|
R{Time} = maprefpostings(xWorldLimits{Time},yWorldLimits{Time},[rasterSize{Time}(2) rasterSize{Time}(1)])
|
|
[geoMeshX{Time},geoMeshY{Time}] = meshgrid(R{Time}.XWorldLimits(1):10000:R{Time}.XWorldLimits(2),...
|
|
R{Time}.YWorldLimits(1):10000:R{Time}.YWorldLimits(2));
|
|
|
|
geoDat{Time} = -griddata(reshape(X{Time}(2:end-2,2:end-2),1,[]),...
|
|
reshape(Y{Time}(2:end-2,2:end-2),1,[]),...
|
|
reshape(Val{Time}(2:end-2,2:end-2),1,[]),...
|
|
geoMeshX{Time},geoMeshY{Time},'natural');
|
|
|
|
geotiffwrite(['Prep_' num2str(timeIDX(Time)) '.tif'],geoDat{Time},...
|
|
R{Time},'CoordRefSysCode',32618)
|
|
disp(Time)
|
|
% clear X Y Val geoDat geoMeshX geoMeshY
|
|
end
|
|
|
|
clear modelWaterPrep
|
|
|
|
%% Export to GEOSERVER
|
|
copyfile *.tif \\geoserver.engineering.queensu.ca\TIFFS\
|
|
copyfile *.dbf \\geoserver.engineering.queensu.ca\TIFFS\
|
|
copyfile *.shx \\geoserver.engineering.queensu.ca\TIFFS\
|
|
copyfile *.shp \\geoserver.engineering.queensu.ca\TIFFS\
|
|
copyfile *.prj \\geoserver.engineering.queensu.ca\TIFFS\
|
|
|
|
|
|
%% Import Model forecast
|
|
modelWaterHist = qpread(qpfopen('trih-rt_run.dat'),1,'water level','data',0,0);
|
|
modelCurHist = qpread(qpfopen('trih-rt_run.dat'),1,'depth averaged velocity','data',0,0);
|
|
|
|
%% Validation import ESTOFS
|
|
load('DUNEX_Nov20.mat')
|
|
clear ESTOFS lonESTOFS lonESTOFS
|
|
|
|
wLIDX=[6 7 8 38];
|
|
wlGrid = [5000 -4000 0 -2000];
|
|
wlGridY = [0 0 2000 5000]; %7500
|
|
ESTOFStimeOffset = [0 12 24 48];
|
|
ESTOFSrunLength = [48 60 72 96];
|
|
|
|
for e = [1 4]
|
|
try
|
|
clear meshX meshY
|
|
m_proj('lambert conformal conic','ori',[-95 25.0],'clo',-95,'par',[25 25],'ell','sphere');
|
|
ESTOFS = ncgeodataset(['D:\DUNEX_RT\Archive\ESTOFS\ESTOFS_' ...
|
|
datestr(masterTime-hours(ESTOFStimeOffset(e)),'YYYY-mm-DD_HH') '_f' num2str(0,'%03d') '.grib']);
|
|
|
|
[meshX,meshY] = meshgrid(ESTOFS{'x'}(:),ESTOFS{'y'}(:));
|
|
[lonESTOFS,latESTOFS] = m_xy2ll(meshX.*1000,meshY.*1000);
|
|
clear meshX meshY
|
|
ESTOFS_Val{e} = 1;
|
|
catch
|
|
ESTOFS_Val{e} = 0;
|
|
end
|
|
end
|
|
clear waterLevelsESTOFSINTER waterLevelsESTOFS
|
|
for e = [1 4]
|
|
try
|
|
for Stat = 1:4
|
|
if ESTOFS_Val{e} == 0
|
|
estofsWL(:,Stat,e) = zeros(49,1);
|
|
else
|
|
[wlDeg(2),wlDeg(1)] = utm2deg(pointsCS{wLIDX(Stat),2}+wlGrid(Stat),pointsCS{wLIDX(Stat),3}+wlGridY(Stat),'18 S');
|
|
[estofIDX] = NearestValue(wlDeg,lonESTOFS,latESTOFS); %Lat
|
|
|
|
for i = 1:ESTOFSrunLength(e)+1
|
|
ESTOFS = ncgeodataset(['D:\DUNEX_RT\Archive\ESTOFS\ESTOFS_' ...
|
|
datestr(masterTime-hours(ESTOFStimeOffset(e)),'YYYY-mm-DD_HH') '_f' num2str(i-1,'%03d') '.grib']);
|
|
|
|
waterLevelsESTOFS{e}(i,Stat) = squeeze(ESTOFS{'Ocean_Surface_Elevation_Relative_to_Geoid_surface'}(1,estofIDX(1),estofIDX(2)));
|
|
end
|
|
waterLevelsESTOFSINTER{e}(:,Stat) = interp1(masterTime-hours(ESTOFStimeOffset(e)):minutes(60):masterTime+hours(48),...
|
|
waterLevelsESTOFS{e}(:,Stat),masterTime-hours(ESTOFStimeOffset(e)):minutes(15):masterTime+hours(48),'spline');
|
|
end
|
|
end
|
|
catch
|
|
waterLevelsESTOFSINTER{e}(:,Stat) = zeros(1,385);
|
|
end
|
|
end
|
|
|
|
|
|
%% Validation import WL
|
|
% Duck
|
|
% Oregon Inlet
|
|
% Hatteras
|
|
% Beaufort
|
|
clear stations statWL
|
|
options = weboptions;
|
|
options.Timeout = 30;
|
|
|
|
stations = {8651370,8652587,8654467,8656483};
|
|
dateStartM = masterTime-hours(48);
|
|
dateEndM = masterTime+hours(5);
|
|
|
|
for Stat = 1:4
|
|
try
|
|
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,options);
|
|
|
|
statWL(:,Stat) = interp1(wlIN.DateTime,wlIN.WaterLevel,dateStartM:minutes(15):dateEndM);
|
|
catch err
|
|
%open file
|
|
fid = fopen('errorFile.log','a+');
|
|
% write the error to file
|
|
% first line: message
|
|
fprintf(fid,'Line %d on %s : %s\n',err.stack.line,datestr(masterTime),err.message);
|
|
sendmail('alexander.rey@queensu.ca',['Ext-DUNEX-RT: ' datestr(masterTime) err.message]) ;
|
|
|
|
% close file
|
|
fclose(fid)
|
|
|
|
statWL(:,Stat) = 0;
|
|
end
|
|
|
|
|
|
clear wlIN
|
|
|
|
end
|
|
%% Validation import Currents
|
|
% Duck
|
|
clear statCur
|
|
for Stat = 1
|
|
% url = 'https://chlthredds.erdc.dren.mil/thredds/dodsC/frf/oceanography/currents/awac-11m/awac-11m.ncml';
|
|
% curTime = ncread(url,'time');
|
|
% curVal = ncread(url,'currentSpeed');
|
|
try
|
|
if day(masterTime)<2
|
|
websave('statCNC.nc',['https://chlthredds.erdc.dren.mil/thredds/fileServer/frf/oceanography/currents/awac-11m/' datestr(masterTime-day(3),'YYYY') '/FRF-ocean_currents_awac-11m_' ...
|
|
datestr(masterTime-days(3),'YYYYmm') '.nc'],options);
|
|
curTime = ncread('statCNC.nc','time');
|
|
curVal = ncread('statCNC.nc','currentSpeed');
|
|
|
|
websave('statCNC.nc',['https://chlthredds.erdc.dren.mil/thredds/fileServer/frf/oceanography/currents/awac-11m/' datestr(masterTime,'YYYY') '/FRF-ocean_currents_awac-11m_' ...
|
|
datestr(masterTime,'YYYYmm') '.nc'],options);
|
|
curTime = vertcat(curTime,ncread('statCNC.nc','time'));
|
|
curVal = vertcat(curVal,ncread('statCNC.nc','currentSpeed'));
|
|
|
|
else
|
|
websave('statCNC.nc',['https://chlthredds.erdc.dren.mil/thredds/fileServer/frf/oceanography/currents/awac-11m/' datestr(masterTime,'YYYY') '/FRF-ocean_currents_awac-11m_' ...
|
|
datestr(masterTime,'YYYYmm') '.nc'],options);
|
|
|
|
curTime = ncread('statCNC.nc','time');
|
|
curVal = ncread('statCNC.nc','currentSpeed');
|
|
end
|
|
|
|
statCur(:) = interp1(datetime(curTime,'convertFrom','POSIX'),curVal,dateStartM:minutes(15):dateEndM);
|
|
clear wlIN
|
|
catch err
|
|
%open file
|
|
fid = fopen('errorFile.log','a+');
|
|
% write the error to file
|
|
% first line: message
|
|
fprintf(fid,'Line %d on %s : %s\n',err.stack.line,datestr(masterTime),err.message);
|
|
sendmail('alexander.rey@queensu.ca',['Ext-DUNEX-RT: ' datestr(masterTime) err.message]) ;
|
|
|
|
% close file
|
|
fclose(fid)
|
|
|
|
statCur(:) = zeros(length(dateStartM:minutes(10):dateEndM),1);
|
|
end
|
|
end
|
|
|
|
%% Validation import Waves
|
|
% Virginia Beach
|
|
% Duck 26
|
|
% Nages Head
|
|
% Oregon Inlet
|
|
% Diamond Shoals
|
|
clear stations
|
|
|
|
stations = {44014,44100,44086,44095,41025};
|
|
dateStartM = masterTime-hours(48);
|
|
dateEndM = masterTime+hours(5);
|
|
|
|
clear statWave
|
|
|
|
for Stat = 1:5
|
|
try
|
|
if ismember(Stat,[1,5]) == 1
|
|
url = ['https://www.ndbc.noaa.gov/data/5day2/' num2str(stations{Stat}) '_5day.txt'];
|
|
|
|
websave('NDBCdata.txt',url,options);
|
|
opts = delimitedTextImportOptions("NumVariables", 39);
|
|
|
|
% Specify range and delimiter
|
|
opts.DataLines = [3, Inf];
|
|
opts.Delimiter = " ";
|
|
|
|
% Specify column names and types
|
|
opts.VariableNames = ["YY", "MM", "DD", "hh", "mm", "WDIR", "WSPD", "GST", "WVHT", "DPD", "APD", "MWD", "PRES", "ATMP", "WTMP", "DEWP", "VIS", "PTDY", "TIDE", "Var20", "Var21", "Var22", "Var23", "Var24", "Var25", "Var26", "Var27", "Var28", "Var29", "Var30", "Var31", "Var32", "Var33", "Var34", "Var35", "Var36", "Var37", "Var38", "Var39"];
|
|
opts.SelectedVariableNames = ["YY", "MM", "DD", "hh", "mm", "WDIR", "WSPD", "GST", "WVHT", "DPD", "APD", "MWD", "PRES", "ATMP", "WTMP", "DEWP", "VIS", "PTDY", "TIDE"];
|
|
opts.VariableTypes = ["double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "string", "string", "string", "string", "string", "string", "string", "string", "string", "string", "string", "string", "string", "string", "string", "string", "string", "string", "string", "string"];
|
|
opts = setvaropts(opts, [20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39], "WhitespaceRule", "preserve");
|
|
opts = setvaropts(opts, 15, "TrimNonNumeric", true);
|
|
opts = setvaropts(opts, 15, "ThousandsSeparator", ",");
|
|
opts = setvaropts(opts, [20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39], "EmptyFieldRule", "auto");
|
|
opts.ExtraColumnsRule = "ignore";
|
|
opts.EmptyLineRule = "read";
|
|
opts.ConsecutiveDelimitersRule = "join";
|
|
opts.LeadingDelimitersRule = "ignore";
|
|
|
|
else
|
|
url = ['https://www.ndbc.noaa.gov/data/realtime2/' num2str(stations{Stat}) '.txt'];
|
|
websave('NDBCdata.txt',url);
|
|
|
|
opts = delimitedTextImportOptions("NumVariables", 52);
|
|
|
|
% Specify range and delimiter
|
|
opts.DataLines = [3, Inf];
|
|
opts.Delimiter = " ";
|
|
|
|
% Specify column names and types
|
|
opts.VariableNames = ["YY", "MM", "DD", "hh", "mm", "WDIR", "WSPD", "GST", "WVHT", "DPD", "APD", "MWD", "PRES", "ATMP", "WTMP", "DEWP", "VIS", "PTDY", "TIDE", "Var20", "Var21", "Var22", "Var23", "Var24", "Var25", "Var26", "Var27", "Var28", "Var29", "Var30", "Var31", "Var32", "Var33", "Var34", "Var35", "Var36", "Var37", "Var38", "Var39", "Var40", "Var41", "Var42", "Var43", "Var44", "Var45", "Var46", "Var47", "Var48", "Var49", "Var50", "Var51", "Var52"];
|
|
opts.SelectedVariableNames = ["YY", "MM", "DD", "hh", "mm", "WDIR", "WSPD", "GST", "WVHT", "DPD", "APD", "MWD", "PRES", "ATMP", "WTMP", "DEWP", "VIS", "PTDY", "TIDE"];
|
|
opts.VariableTypes = ["double", "double", "double", "double", "double", "categorical", "categorical", "categorical", "double", "double", "double", "double", "categorical", "categorical", "double", "double", "categorical", "categorical", "categorical", "string", "string", "string", "string", "string", "string", "string", "string", "string", "string", "string", "string", "string", "string", "string", "string", "string", "string", "string", "string", "string", "string", "string", "string", "string", "string", "string", "string", "string", "string", "string", "string", "string"];
|
|
opts = setvaropts(opts, [20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52], "WhitespaceRule", "preserve");
|
|
opts = setvaropts(opts, [9, 10, 11, 12, 15], "TrimNonNumeric", true);
|
|
opts = setvaropts(opts, [9, 10, 11, 12, 15], "ThousandsSeparator", ",");
|
|
opts = setvaropts(opts, [6, 7, 8, 13, 14, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52], "EmptyFieldRule", "auto");
|
|
opts.ExtraColumnsRule = "ignore";
|
|
opts.EmptyLineRule = "read";
|
|
opts.ConsecutiveDelimitersRule = "join";
|
|
opts.LeadingDelimitersRule = "ignore";
|
|
|
|
end
|
|
|
|
% Import the data
|
|
waveIN = readtable("NDBCdata.txt", opts);
|
|
waveNAN = isnan(waveIN.WVHT)==0;
|
|
% Hs
|
|
statWave(:,1,Stat) = interp1(datetime(waveIN.YY(waveNAN),(waveIN.MM(waveNAN)),(waveIN.DD(waveNAN)),...
|
|
(waveIN.hh(waveNAN)),(waveIN.mm(waveNAN)),zeros(length(waveIN.YY(waveNAN)),1)),...
|
|
(waveIN.WVHT(waveNAN)),dateStartM:minutes(15):dateEndM);
|
|
% DPD
|
|
% statWave(:,2,Stat) = interp1(datetime(waveIN.YY(waveNAN),(waveIN.MM(waveNAN)),(waveIN.DD(waveNAN)),...
|
|
% (waveIN.hh(waveNAN)),(waveIN.mm(waveNAN)),zeros(length(waveIN.YY(waveNAN)),1)),...
|
|
% (waveIN.DPD(waveNAN)),dateStartM:minutes(10):dateEndM);
|
|
% APD
|
|
statWave(:,3,Stat) = interp1(datetime(waveIN.YY(waveNAN),(waveIN.MM(waveNAN)),(waveIN.DD(waveNAN)),...
|
|
(waveIN.hh(waveNAN)),(waveIN.mm(waveNAN)),zeros(length(waveIN.YY(waveNAN)),1)),...
|
|
(waveIN.APD(waveNAN)),dateStartM:minutes(15):dateEndM);
|
|
%MWD
|
|
statWave(:,4,Stat) = interp1(datetime(waveIN.YY(waveNAN),(waveIN.MM(waveNAN)),(waveIN.DD(waveNAN)),...
|
|
(waveIN.hh(waveNAN)),(waveIN.mm(waveNAN)),zeros(length(waveIN.YY(waveNAN)),1)),...
|
|
(waveIN.MWD(waveNAN)),dateStartM:minutes(15):dateEndM);
|
|
clear waveIN opts
|
|
catch err
|
|
%open file
|
|
fid = fopen('errorFile.log','a+');
|
|
% write the error to file
|
|
% first line: message
|
|
%fprintf(fid,'Line %d on %s : %s\n',err.stack.line,datestr(masterTime),err.message);
|
|
%sendmail('alexander.rey@queensu.ca',['Ext-DUNEX-RT: ' datestr(masterTime) err.message]);
|
|
|
|
% close file
|
|
fclose(fid)
|
|
statWave(:,:,Stat) = zeros(size(dateStartM:minutes(10):dateEndM,2),4,1);
|
|
end
|
|
end
|
|
statWave(isnan(statWave))=0;
|
|
%% Import Model Waves
|
|
|
|
for i = 1:49
|
|
startRow = 8;
|
|
formatSpec = '%13f%14f%14f%14f%14f%14f%14f%14f%14f%14f%14f%14f%f%[^\n\r]';
|
|
|
|
fileID = fopen(['Dec28.loct00000' num2str(i,'%02d') '.tab'],'r');
|
|
|
|
% This call is based on the structure of the file used to generate this
|
|
% code. If an error occurs for a different file, try regenerating the code
|
|
% from the Import Tool.
|
|
dataArray = textscan(fileID, formatSpec, 'Delimiter', '', 'WhiteSpace', '', 'TextType', 'string', 'HeaderLines' ,startRow-1, 'ReturnOnError', false, 'EndOfLine', '\r\n');
|
|
|
|
fclose(fileID);
|
|
modelWaveIN= table(dataArray{1:end-1}, 'VariableNames', {'Xp','Yp','Depth','Hsig','Dir','RTpea','kTm01','Dspr','Ubot','XWindv','YWindv','XVel','YVel'});
|
|
|
|
modelWave(:,1,i) = modelWaveIN.Hsig;
|
|
modelWave(:,3,i) = modelWaveIN.RTpea;
|
|
modelWave(:,4,i) = modelWaveIN.Dir;
|
|
|
|
|
|
clearvars filename startRow formatSpec fileID dataArray ans modelWaveIN;
|
|
|
|
end
|
|
|
|
%% Import Model WL
|
|
modelWL = modelWaterHist.Val;
|
|
modelCur = sqrt(modelCurHist.XComp.^2+modelCurHist.YComp.^2);
|
|
|
|
for idx = 1:size(modelWL,2)
|
|
modelWLInter(:,idx) = interp1(modelWaterHist.Time,modelWL(:,idx),datenum(masterTime:minutes(15):masterTime+hours(48)));
|
|
end
|
|
for idx = 1:size(modelCur,2)
|
|
modelCurInter(:,idx) = sqrt(interp1(modelCurHist.Time,modelCurHist.XComp(:,idx),datenum(masterTime:minutes(15):masterTime+hours(48))).^2+...
|
|
interp1(modelCurHist.Time,modelCurHist.YComp(:,idx),datenum(masterTime:minutes(15):masterTime+hours(48))).^2);
|
|
end
|
|
|
|
save(['ModelOut_' datestr(masterTime,'dd-mmm-yyyy_HH') 'Z.mat'],'modelWL','modelWave','statWL','statWave','modelCur','statCur','modelWLInter','modelCurInter');
|
|
|
|
%% Import Historic Model WL
|
|
if masterTime>=datetime(2021,07,08)
|
|
try
|
|
H12 = load(['ModelOut_' datestr(masterTime-hours(12),'dd-mmm-yyyy_HH') 'Z.mat']);
|
|
catch
|
|
end
|
|
try
|
|
H24 = load(['ModelOut_' datestr(masterTime-hours(24),'dd-mmm-yyyy_HH') 'Z.mat']);
|
|
catch
|
|
end
|
|
try
|
|
H48 = load(['ModelOut_' datestr(masterTime-hours(48),'dd-mmm-yyyy_HH') 'Z.mat']);
|
|
catch
|
|
end
|
|
end
|
|
%% Save html for Currents
|
|
googleDates = days(masterTime-hours(48)+minutes((0:384)*15)-datetime(1899,12,30,00,00,00));
|
|
% mat2sheets('1K6toMh7L8nL3rtoPLJxD1Ebfqaty0wkPJSUGXNAbhoU','1763350398',[2,1],googleDates')
|
|
googleDatesCell(2:length(googleDates)+1) = num2cell(googleDates);
|
|
googleDatesCell{1} = 'Time';
|
|
locLine = 1:5:67;
|
|
locationsWave ={'FRF 11 m AWAC '};
|
|
|
|
% % for i = 1:433
|
|
% % for j = 1:100
|
|
% % modelDavDel(i,j)={''};
|
|
% % end
|
|
% % end
|
|
% % mat2sheets('1K6toMh7L8nL3rtoPLJxD1Ebfqaty0wkPJSUGXNAbhoU','1763350398',[2,2],modelDavDel)
|
|
|
|
davIDX=[17 161:172];
|
|
|
|
clear modelCurPlot curValLoc
|
|
for Loc = 1:13
|
|
for Line=1:5
|
|
switch Line
|
|
case 1
|
|
curValLoc{Line} = ['Forecast'];
|
|
if exist('modelCur','var') ==0
|
|
modelCurPlot(1:385,locLine(Loc)+0) = num2cell(zeros(385,1));
|
|
continue
|
|
end
|
|
modelCurPlot(1:192,locLine(Loc)+0) = {''};
|
|
modelCurPlot(193:1:385,locLine(Loc)+0) = num2cell(modelCurInter(:,davIDX(Loc)));
|
|
case 2
|
|
curValLoc{Line} = ['-12 Hour: ' datestr(masterTime-hours(12),'mm/dd HH') 'Z Run'];
|
|
if exist('H12','var') ==0
|
|
modelCurPlot(1:385,locLine(Loc)+1) = num2cell(zeros(385,1));
|
|
continue
|
|
end
|
|
modelCurPlot(1:143,locLine(Loc)+1) = {''};
|
|
modelCurPlot(337:385,locLine(Loc)+1) = {''};
|
|
modelCurPlot(144:1:336,locLine(Loc)+1) = num2cell(H12.modelCurInter(:,davIDX(Loc)));
|
|
|
|
case 3
|
|
curValLoc{Line} = ['-24 Hour: ' datestr(masterTime-hours(24),'mm/dd HH') 'Z Run'];
|
|
if exist('H24','var') ==0
|
|
modelCurPlot(1:385,locLine(Loc)+2) = num2cell(zeros(385,1));
|
|
continue
|
|
end
|
|
modelCurPlot(1:96,locLine(Loc)+2) = {''};
|
|
modelCurPlot(290:385,locLine(Loc)+2) = {''};
|
|
modelCurPlot(97:1:289,locLine(Loc)+2) = num2cell(H24.modelCurInter(:,davIDX(Loc)));
|
|
|
|
case 4
|
|
curValLoc{Line} = ['-48 Hour: ' datestr(masterTime-hours(48),'mm/dd HH') 'Z Run'];
|
|
if exist('H48','var') ==0
|
|
modelCurPlot(1:385,locLine(Loc)+3) = num2cell(zeros(385,1));
|
|
continue
|
|
end
|
|
modelCurPlot(193:385,locLine(Loc)+3) = {''};
|
|
modelCurPlot(1:1:192,locLine(Loc)+3) = num2cell(H48.modelCurInter(:,davIDX(Loc)));
|
|
case 5
|
|
curValLoc{Line} = ['Observed Velocity'];
|
|
if Loc> 1 || isnan(statCur(1,Loc))
|
|
modelCurPlot(1:385,locLine(Loc)+4) = num2cell(zeros(385,1));
|
|
continue
|
|
end
|
|
modelCurPlot(214:385,locLine(Loc)+4) = {''};
|
|
modelCurPlot(1:213,locLine(Loc)+4) = num2cell(statCur);
|
|
|
|
end
|
|
|
|
end
|
|
end
|
|
|
|
% mat2sheets('1K6toMh7L8nL3rtoPLJxD1Ebfqaty0wkPJSUGXNAbhoU','1763350398',...
|
|
% [1,locLine(Loc)],modelCurPlot)
|
|
%
|
|
% mat2sheets('1K6toMh7L8nL3rtoPLJxD1Ebfqaty0wkPJSUGXNAbhoU','1763350398',...
|
|
% [1,locLine(Loc)],curValLoc)
|
|
curValLoc = repmat(curValLoc,1,Loc);
|
|
|
|
mat2sheets('1K6toMh7L8nL3rtoPLJxD1Ebfqaty0wkPJSUGXNAbhoU','1763350398',...
|
|
[1,1],horzcat(googleDatesCell',vertcat(curValLoc(1:size(modelCurPlot,2)),modelCurPlot)))
|
|
%% Save html for WL
|
|
googleDates = days(masterTime-hours(48)+minutes((0:384)*15)-datetime(1899,12,30,00,00,00));
|
|
% mat2sheets('1K6toMh7L8nL3rtoPLJxD1Ebfqaty0wkPJSUGXNAbhoU','1102090436',[2,1],googleDates')
|
|
googleDatesCell(2:length(googleDates)+1) = num2cell(googleDates);
|
|
googleDatesCell{1} = 'Time';
|
|
|
|
% locLine = [2 7 12 17 22 27];
|
|
% locLine = [2 9 16 23 30 37];
|
|
locLine = 1:7:114;
|
|
|
|
locationsWl ={'FRF Pier',...
|
|
'Oregon Inlet',...
|
|
'Hatteras',...
|
|
'Beaufort'};
|
|
|
|
% % for i = 1:433
|
|
% % for j = 1:121
|
|
% % modelWlDel(i,j)={''};
|
|
% % end
|
|
% % end
|
|
% % mat2sheets('1K6toMh7L8nL3rtoPLJxD1Ebfqaty0wkPJSUGXNAbhoU','1102090436',[2,2],modelWlDel)
|
|
|
|
wLIDX=[6 7 8 38 161:172];
|
|
clear modelWlPlot wlValLoc
|
|
for Loc = 1:length(wLIDX)
|
|
% modelwlPlot = cell(432,4);
|
|
for Line=[1:5 6 9]
|
|
switch Line
|
|
case 1
|
|
wlValLoc{Line} = ['Forecast'];
|
|
if exist('modelCur','var') ==0
|
|
modelwlPlot(1:385,locLine(Loc)+0) = num2cell(zeros(385,1));
|
|
continue
|
|
end
|
|
modelwlPlot(1:192,locLine(Loc)+0) = {''};
|
|
modelwlPlot(193:1:385,locLine(Loc)+0) = num2cell(modelWLInter(:,wLIDX(Loc)));
|
|
case 2
|
|
wlValLoc{Line} = ['-12 Hour: ' datestr(masterTime-hours(12),'mm/dd HH') 'Z Run'];
|
|
if exist('H12','var') ==0
|
|
modelwlPlot(1:385,locLine(Loc)+1) = num2cell(zeros(385,1));
|
|
continue
|
|
end
|
|
modelwlPlot(1:143,locLine(Loc)+1) = {''};
|
|
modelwlPlot(337:385,locLine(Loc)+1) = {''};
|
|
modelwlPlot(144:1:336,locLine(Loc)+1) = num2cell(H12.modelWLInter(:,wLIDX(Loc)));
|
|
|
|
case 3
|
|
wlValLoc{Line} = ['-24 Hour: ' datestr(masterTime-hours(24),'mm/dd HH') 'Z Run'];
|
|
if exist('H24','var') ==0
|
|
modelwlPlot(1:385,locLine(Loc)+2) = num2cell(zeros(385,1));
|
|
continue
|
|
end
|
|
modelwlPlot(1:96,locLine(Loc)+2) = {''};
|
|
modelwlPlot(290:385,locLine(Loc)+2) = {''};
|
|
modelwlPlot(97:1:289,locLine(Loc)+2) = num2cell(H24.modelWLInter(:,wLIDX(Loc)));
|
|
|
|
case 4
|
|
wlValLoc{Line} = ['-48 Hour: ' datestr(masterTime-hours(48),'mm/dd HH') 'Z Run'];
|
|
if exist('H48','var') ==0
|
|
modelwlPlot(1:385,locLine(Loc)+3) = num2cell(zeros(385,1));
|
|
continue
|
|
end
|
|
modelwlPlot(218:385,locLine(Loc)+3) = {''};
|
|
modelwlPlot(1:1:217,locLine(Loc)+3) = num2cell(H48.modelWLInter(:,wLIDX(Loc)));
|
|
case 5
|
|
wlValLoc{Line} = ['Observed Water Level'];
|
|
if Loc>4 || isnan(statWL(1,Loc))
|
|
modelwlPlot(1:385,locLine(Loc)+4) = num2cell(zeros(385,1));
|
|
continue
|
|
end
|
|
modelwlPlot(214:385,locLine(Loc)+4) = {''};
|
|
modelwlPlot(1:213,locLine(Loc)+4) = num2cell(statWL(:,Loc));
|
|
case 6
|
|
wlValLoc{Line} = ['ESTOFS Forecast'];
|
|
if Loc>4 || ESTOFS_Val{1} == 0
|
|
modelwlPlot(1:385,locLine(Loc)+5) = num2cell(zeros(385,1));
|
|
continue
|
|
end
|
|
modelwlPlot(1:192,locLine(Loc)+5) = {''};
|
|
modelwlPlot(193:385,locLine(Loc)+5) = num2cell(waterLevelsESTOFSINTER{1}(:,Loc));
|
|
case 9
|
|
wlValLoc{7} = ['-48 Hour ESTOFS'];
|
|
if Loc>4 || ESTOFS_Val{4} == 0 || size(ESTOFS_Val{4},2) == 1
|
|
modelwlPlot(1:385,locLine(Loc)+6) = num2cell(zeros(385,1));
|
|
continue
|
|
end
|
|
modelwlPlot(1:1:385,locLine(Loc)+6) = num2cell(waterLevelsESTOFSINTER{4}(:,Loc));
|
|
end
|
|
end
|
|
disp(Loc)
|
|
|
|
end
|
|
|
|
% mat2sheets('1K6toMh7L8nL3rtoPLJxD1Ebfqaty0wkPJSUGXNAbhoU','1102090436',...
|
|
% [2,locLine(Loc)],modelwlPlot)
|
|
%
|
|
% mat2sheets('1K6toMh7L8nL3rtoPLJxD1Ebfqaty0wkPJSUGXNAbhoU','1102090436',...
|
|
% [1,locLine(Loc)],wlValLoc)
|
|
|
|
wlValLoc = repmat(wlValLoc,1,Loc);
|
|
|
|
mat2sheets('1K6toMh7L8nL3rtoPLJxD1Ebfqaty0wkPJSUGXNAbhoU','1102090436',...
|
|
[1,1],horzcat(googleDatesCell',vertcat(wlValLoc(1:size(modelwlPlot,2)),modelwlPlot)))
|
|
|
|
%% Save html for Waves
|
|
googleDates = days(masterTime-hours(48)+minutes((0:384)*15)-datetime(1899,12,30,00,00,00));
|
|
googleDatesCell(2:length(googleDates)+1) = num2cell(googleDates);
|
|
googleDatesCell{1} = 'Time';
|
|
|
|
% mat2sheets('1K6toMh7L8nL3rtoPLJxD1Ebfqaty0wkPJSUGXNAbhoU','319734265',[2,1],googleDates')
|
|
|
|
% locLine = [2 7 12 17 22 27];
|
|
locLine = 1:5:114;
|
|
|
|
locationsWave ={'Virginia Beach',...
|
|
'FRF 26 m',...
|
|
'Nags Head',...
|
|
'Oregon Inlet',...
|
|
'Diamond Shoals'};
|
|
|
|
clear waveValLoc modelwavePlot
|
|
|
|
for i = 1:385
|
|
for j = 1:85
|
|
modelwavePlot(i,j)={''};
|
|
end
|
|
end
|
|
% % mat2sheets('1K6toMh7L8nL3rtoPLJxD1Ebfqaty0wkPJSUGXNAbhoU','319734265',[2,2],modelWaveDel)
|
|
|
|
waveIDX = [12 13 14 15 16 161:172];
|
|
for Loc = 1:length(waveIDX)
|
|
% modelwavePlot = cell(432,5);
|
|
for Line=1:5
|
|
switch Line
|
|
case 1
|
|
waveValLoc{Line} = ['Forecast'];
|
|
if exist('modelWave','var') ==0
|
|
modelwavePlot(1:385,locLine(Loc)+0) = num2cell(zeros(433,1));
|
|
continue
|
|
end
|
|
modelwavePlot(1:191,locLine(Loc)+0) = {''};
|
|
modelwavePlot(192,locLine(Loc)+0) = {0};
|
|
modelwavePlot(193:4:385,locLine(Loc)+0) = num2cell(modelWave(waveIDX(Loc),1,:));
|
|
case 2
|
|
waveValLoc{Line} = ['-12 Hour: ' datestr(masterTime-hours(12),'mm/dd HH') 'Z Run'];
|
|
if exist('H12','var') ==0 || size(H12.modelWave,1)<waveIDX(Loc)
|
|
modelwavePlot(1:385,locLine(Loc)+1) = num2cell(zeros(385,1));
|
|
continue
|
|
end
|
|
modelwavePlot(1:142,locLine(Loc)+1) = {''};
|
|
modelwavePlot(338:385,locLine(Loc)+1) = {''};
|
|
modelwavePlot(143,locLine(Loc)+1) = {0};
|
|
modelwavePlot(337,locLine(Loc)+1) = {0};
|
|
modelwavePlot(144:4:336,locLine(Loc)+1) = num2cell(H12.modelWave(waveIDX(Loc),1,:));
|
|
|
|
case 3
|
|
waveValLoc{Line} = ['-24 Hour: ' datestr(masterTime-hours(24),'mm/dd HH') 'Z Run'];
|
|
if exist('H24','var') ==0 || size(H24.modelWave,1)<waveIDX(Loc)
|
|
modelwavePlot(1:385,locLine(Loc)+2) = num2cell(zeros(385,1));
|
|
continue
|
|
end
|
|
modelwavePlot(1:95,locLine(Loc)+2) = {''};
|
|
modelwavePlot(291:385,locLine(Loc)+2) = {''};
|
|
modelwavePlot(96,locLine(Loc)+2) = {0};
|
|
modelwavePlot(290,locLine(Loc)+2) = {0};
|
|
modelwavePlot(97:4:289,locLine(Loc)+2) = num2cell(H24.modelWave(waveIDX(Loc),1,:));
|
|
|
|
case 4
|
|
waveValLoc{Line} = ['-48 Hour: ' datestr(masterTime-hours(36),'mm/dd HH') 'Z Run'];
|
|
if exist('H48','var') ==0 || size(H36.modelWave,1)<waveIDX(Loc)
|
|
modelwavePlot(1:385,locLine(Loc)+3) = num2cell(zeros(385,1));
|
|
continue
|
|
end
|
|
modelwavePlot(218:385,locLine(Loc)+3) = {''};
|
|
modelwavePlot(218,locLine(Loc)+3) = {0};
|
|
modelwavePlot(1:6:217,locLine(Loc)+3) = num2cell(H36.modelWave(waveIDX(Loc),1,:));
|
|
case 5
|
|
waveValLoc{Line} = ['Observed Wave'];
|
|
if Loc>5 || isnan(statWave(1,1,Loc))
|
|
modelwavePlot(1:385,locLine(Loc)+4) = num2cell(zeros(385,1));
|
|
continue
|
|
end
|
|
modelwavePlot(215:385,locLine(Loc)+4) = {''};
|
|
modelwavePlot(214,locLine(Loc)+4) = {0};
|
|
modelwavePlot(1:1:213,locLine(Loc)+4) = num2cell(statWave(:,1,Loc));
|
|
|
|
end
|
|
|
|
end
|
|
% mat2sheets('1K6toMh7L8nL3rtoPLJxD1Ebfqaty0wkPJSUGXNAbhoU','319734265',...
|
|
% [2,locLine(Loc)],modelwavePlot)
|
|
%
|
|
% mat2sheets('1K6toMh7L8nL3rtoPLJxD1Ebfqaty0wkPJSUGXNAbhoU','319734265',...
|
|
% [1,locLine(Loc)],waveValLoc)
|
|
end
|
|
|
|
waveValLoc = repmat(waveValLoc,1,Loc);
|
|
|
|
mat2sheets('1K6toMh7L8nL3rtoPLJxD1Ebfqaty0wkPJSUGXNAbhoU','319734265',...
|
|
[1,1],horzcat(googleDatesCell',vertcat(waveValLoc(1:size(modelwavePlot,2)),modelwavePlot)))
|
|
|
|
|
|
|
|
|
|
%% Import and modify times
|
|
% modelHour = [0 3 6 12 24 48];
|
|
modelHour = [0 6 12 24 36 48];
|
|
for i = 1:6
|
|
dateText(i) = {[datestr(masterTime+hours(modelHour(i)),'mmm dd, HH:MM') ' UTC (Hour '...
|
|
num2str(modelHour(i),'%d') ')']};
|
|
end
|
|
mat2sheets('1K6toMh7L8nL3rtoPLJxD1Ebfqaty0wkPJSUGXNAbhoU','0',[1,1],...
|
|
dateText)
|
|
|
|
catch err
|
|
%open file
|
|
fid = fopen('errorFile.log','a+');
|
|
% write the error to file
|
|
% first line: message
|
|
fprintf(fid,'Line %d on %s : %s\n',err.stack.line,datestr(masterTime),err.message);
|
|
sendmail('alexander.rey@queensu.ca',['ALL-DUNEX-RT: ' datestr(masterTime) err.message]) ;
|
|
|
|
% close file
|
|
fclose(fid)
|
|
|
|
end
|
|
|
|
%%
|
|
% quit
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|