clear all
close all
fclose all;

%% Read File
disp('Reading Wave Forcing')
radsfile = 'rads.64.nc';
ncinfo(radsfile);
X = ncread(radsfile,'x'); 
Y = ncread(radsfile,'y');
time = ncread(radsfile,'time')/60;  %ADCIRC seconds to RAS minutes
Fx = ncread(radsfile,'radstress_x')*3.2808399^2; %m^2/s to ft^2/s^2
Fy = ncread(radsfile,'radstress_y')*3.2808399^2; %m^2/s to ft^2/s^2

%% Grid
dd = 0.05; % resolution for structured grid output
ax = [ -95 -87.5  27.8   31.6]; %region to clip from data
[Xg,Yg] = meshgrid(ax(1):dd:ax(2),ax(3):dd:ax(4)); %wave forcing grid
ai = find(X>ax(1) & X<ax(2) & Y>ax(3) & Y<ax(4));
Fxg = zeros(length(Xg(:,1)),length(Xg(1,:)),length(time));
Fyg = zeros(length(Xg(:,1)),length(Xg(1,:)),length(time));

%% Interpolation
disp('Interpolating Wave Forcing')
F = scatteredInterpolant(X(ai),Y(ai),Fx(ai,1));
for f = 1:length(time)
    disp(num2str(f/length(time)))
    F.Values = Fx(ai,f);
    Fxg(:,:,f) = F(Xg,Yg); %regular grid
    F.Values = Fy(ai,f);
    Fyg(:,:,f) = F(Xg,Yg); %regular grid
end
nix = isnan(Fxg);
niy = isnan(Fyg);
Fxg(nix) = -999;
Fyg(niy) = -999;


%% White Output File
ncfileout = 'Wave_Forcing.nc'; %name of output file
delete(ncfileout)
nccreate(ncfileout,'WaveRadX','Dimensions',{'x',length(Xg(1,:)),'y',length(Xg(:,1)),'time',length(time)},'Format','classic','datatype','single')
nccreate(ncfileout,'WaveRadY','Dimensions',{'x',length(Xg(1,:)),'y',length(Xg(:,1)),'time',length(time)},'Format','classic','datatype','single')
nccreate(ncfileout,'time','Dimensions',{'time',length(time)},'Format','classic')
nccreate(ncfileout,'x','Dimensions',{'x',length(Xg(1,:))},'Format','classic')
nccreate(ncfileout,'y','Dimensions',{'y',length(Xg(:,1))},'Format','classic')
nccreate(ncfileout,'z','Dimensions',{'x',length(Xg(1,:)),'y',length(Xg(:,1))},'Format','classic')
nccreate(ncfileout,'crs','Format','classic','datatype','int32')
ncwrite(ncfileout,'WaveRadX',rot90(fliplr(Fxg)))
ncwrite(ncfileout,'WaveRadY',rot90(fliplr(Fyg)))
ncwrite(ncfileout,'time',time)
ncwrite(ncfileout,'x',ax(1):dd:ax(2))
ncwrite(ncfileout,'y',ax(3):dd:ax(4))
ncwrite(ncfileout,'z',0*Fxg(:,:,1)')

Ab = {'Conventions','CF-1.6,UGRID-0.9';
    'title','Data';
    'institution','MVN USACE';
    'source','Export NETCDF-CF_GRID from ADCIRC';
    'history','2018-04-20 13:30:52 GMT: exported from ADCIRC';
    'references','http://';
    'Metadata_Conventions','Unidata Dataset Discovery v1.0';
    'summary','Data exported from ADCIRC';
    'date_created','2020-12-20 13:30:52 GMT';
    'fews_implementation_version','2017.01';
    'fews_patch_number','73948';
    'fews_build_number','71514'};

for f = 1:12
    ncwriteatt(ncfileout,'/',char(Ab(f,1)),char(Ab(f,2)));
end

ncwriteatt(ncfileout,'time','standard_name','time')
ncwriteatt(ncfileout,'time','long_name','time')
ncwriteatt(ncfileout,'time','units','minutes since 2021-08-07 00:00:00.0 +0000') %be sure to set to ADCIRC cold start date
ncwriteatt(ncfileout,'time','axis','T')

ncwriteatt(ncfileout,'y','standard_name','latitude')
ncwriteatt(ncfileout,'y','long_name','y coordinate according to WGS 1984')
ncwriteatt(ncfileout,'y','units','degrees_north')
ncwriteatt(ncfileout,'y','axis','Y')
ncwriteatt(ncfileout,'y','_FillValue',9.96921000000000e+36)

ncwriteatt(ncfileout,'x','standard_name','longitude')
ncwriteatt(ncfileout,'x','long_name','x coordinate according to WGS 1984')
ncwriteatt(ncfileout,'x','units','degrees_east')
ncwriteatt(ncfileout,'x','axis','X')
ncwriteatt(ncfileout,'x','_FillValue',9.96921000000000e+36)

ncwriteatt(ncfileout,'z','long_name','height above mean sea level')
ncwriteatt(ncfileout,'z','units','meters')
ncwriteatt(ncfileout,'z','_FillValue',9.96921000000000e+36)

ncwriteatt(ncfileout,'crs','long_name','coordinate reference system')
ncwriteatt(ncfileout,'crs','grid_mapping_name','latitude_longitude')
ncwriteatt(ncfileout,'crs','longitude_of_prime_meridian',0)
ncwriteatt(ncfileout,'crs','semi_major_axis',6378137)
ncwriteatt(ncfileout,'crs','inverse_flattening',298.257223563000)
ncwriteatt(ncfileout,'crs','crs_wkt','GEOGCS["WGS 84", DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]], PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]], UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]], AUTHORITY["EPSG","4326"]]')
ncwriteatt(ncfileout,'crs','proj4_params','+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')
ncwriteatt(ncfileout,'crs','epsg_code','EPSG:4326')

ncwriteatt(ncfileout,'WaveRadX','long_name','Wave Radiation Stress Gradient in X')
ncwriteatt(ncfileout,'WaveRadX','units','ft^2/s^2')
ncwriteatt(ncfileout,'WaveRadX','_FillValue',-999)
ncwriteatt(ncfileout,'WaveRadX','grid_mapping','crs')

ncwriteatt(ncfileout,'WaveRadY','long_name','Wave Radiation Stress Gradient in Y')
ncwriteatt(ncfileout,'WaveRadY','units','ft^2/s^2')
ncwriteatt(ncfileout,'WaveRadY','_FillValue',-999)
ncwriteatt(ncfileout,'WaveRadY','grid_mapping','crs')

return