Skip to content

Commit 75e2b00

Browse files
WPringleWPringle
authored
Merge pull request #231 from CHLNDDEV/sal_netcdf
Support for writing Self Attraction and Loading (SAL) files in NetCDF for the ADCIRC model
2 parents d6d17a1 + ac4c6a2 commit 75e2b00

File tree

4 files changed

+145
-64
lines changed

4 files changed

+145
-64
lines changed

‎@msh/msh.m‎

Lines changed: 5 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -176,9 +176,9 @@
176176
end
177177

178178
% write mesh to disk
179-
function write(obj,fname,type)
179+
function write(obj,fname,type,varargin)
180180
% Usage:
181-
% write(obj,fname,type)
181+
% write(obj,fname,type,varargin)
182182
%
183183
% Examples:
184184
% write(obj); % writes all available data to fort_1.xx (ADCIRC) files
@@ -187,6 +187,7 @@ function write(obj,fname,type)
187187
% write(obj,fname,'gr3'); % writes mesh data to fname.gr3 (SCHISM) file
188188
% write(obj,fname,'ww3'); % writes mesh data to fname.ww3 (WaveWatchIII) file
189189
% write(obj,fname,{'13','14'}); % writes mesh data and f13 attribute data to fname.14 and fname.13 (ADCIRC) files
190+
% write(obj,fname,'24','netcdf'); % writes fort.24 SAL data to fname.24.nc netcdf file
190191
if nargin == 1
191192
fname = 'fort_1';
192193
end
@@ -212,7 +213,7 @@ function write(obj,fname,type)
212213
writefort15( obj.f15, [fname '.15'], obj.bd );
213214
end
214215
if ~isempty(obj.f24)
215-
writefort24( obj.f24, [fname '.24'] );
216+
writefort24( obj.f24, [fname '.24'], obj.p, varargin);
216217
end
217218
if ~isempty(obj.f5354)
218219
writefort5354( obj.f5354, fname );
@@ -268,7 +269,7 @@ function write(obj,fname,type)
268269
writefort19( obj.f2001, [fname '.2001'] );
269270
end
270271
if any(contains(type,'24')) && ~isempty(obj.f24)
271-
writefort24( obj.f24, [fname '.24'] );
272+
writefort24( obj.f24, [fname '.24'], obj.p, varargin);
272273
end
273274
if any(contains(type,'5354')) && ~isempty(obj.f5354)
274275
writefort5354( obj.f5354, fname );

‎@msh/private/writefort24.m‎

Lines changed: 84 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -1,26 +1,92 @@
1-
function fid = writefort24( f24dat, finame )
1+
function fid = writefort24( f24dat, finame, points, format )
22
%
33
%
4-
if ( nargin == 1 )
5-
finputname = 'fort.24_1' ;
4+
if ( nargin == 1 )
5+
finputname = 'fort.24_1' ;
66
else
7-
finputname = finame ;
7+
finputname = finame ;
88
end
99

10-
fid = fopen(finputname,'w') ;
10+
if nargin < 3 || isempty(format)
11+
format = 'ascii';
12+
warning('Using ASCII file format. Pass ''netcdf'' to write in NetCDF')
13+
end
14+
if iscell(format)
15+
format = format{1};
16+
end
17+
18+
if strcmp(format,'ascii') % LEGACY
19+
fid = fopen(finputname,'w') ;
20+
21+
tipnames = f24dat.tiponame; ntip = length(tipnames) ;
22+
for icon = 1: ntip
23+
tipname = tipnames{icon};
24+
fprintf('Writing SAL %s data \n', char(tipname)) ;
25+
% The constituent details
26+
fprintf(fid,'%s \n',[char(tipname) ' SAL']) ;
27+
fprintf(fid,'%17.15f \n',f24dat.omega(icon)) ;
28+
fprintf(fid,'%d \n',1) ;
29+
fprintf(fid,'%s \n',char(tipname)) ;
30+
fprintf(fid,'%d \t %12.6f %12.6f \n',f24dat.Val(icon,:,:));
31+
end
32+
33+
fclose(fid) ;
34+
35+
elseif strcmp(format,'netcdf') % NETCDF
36+
% create nc file or overwrite it
37+
file = [finputname,'.nc'];
38+
if exist(file)
39+
delete(file)
40+
end
41+
% deflate level (set zero for no deflation if necessary.. using single precision so not very different in size)
42+
dl = 5;
43+
44+
node = length(points);
45+
tipnames = char(f24dat.tiponame); ntip = size(tipnames) ;
46+
fillvalue = 0.0;
47+
48+
nccreate(file,'x','Dimensions',{'node',node},'DataType','int32','DeflateLevel',dl) ;
49+
ncwriteatt(file,'x','standard_name', 'longitude') ;
50+
ncwriteatt(file,'x','units', 'degrees_east') ;
51+
ncwriteatt(file,'x','positive', 'east') ;
52+
ncwrite(file, 'x', points(:,1)) ;
1153

12-
tipnames = f24dat.tiponame; ntip = length(tipnames) ;
13-
for icon = 1: ntip
14-
tipname = tipnames{icon};
15-
fprintf('Writing SAL %s data \n', char(tipname)) ;
16-
% The constituent details
17-
fprintf(fid,'%s \n',[char(tipname) ' SAL']) ;
18-
fprintf(fid,'%17.15f \n',f24dat.omega(icon)) ;
19-
fprintf(fid,'%d \n',1) ;
20-
fprintf(fid,'%s \n',char(tipname)) ;
21-
fprintf(fid,'%d \t %12.6f %12.6f \n',f24dat.Val(icon,:,:));
54+
nccreate(file,'y','Dimensions',{'node',node},'DataType','int32','DeflateLevel',dl) ;
55+
ncwriteatt(file,'y','standard_name', 'latitude') ;
56+
ncwriteatt(file,'y','units', 'degrees_north') ;
57+
ncwriteatt(file,'y','positive', 'north') ;
58+
ncwrite(file, 'y', points(:,2)) ;
59+
60+
nccreate(file,'constituents','Dimensions',{'num_constituents',ntip(1),'char_len',ntip(2)},'DataType','char') ;
61+
ncwriteatt(file,'constituents','standard_name', 'name_of_tidal_harmonic_constituents') ;
62+
ncwriteatt(file,'constituents','long_name', 'name of tidal harmonic constituents') ;
63+
ncwrite(file, 'constituents', tipnames) ;
64+
65+
nccreate(file,'frequency','dimensions',{'num_constituents',ntip(1)});
66+
ncwriteatt(file,'frequency','standard_name','frequency_of_tidal_harmonic_constituents')
67+
ncwriteatt(file,'frequency','long_name','frequency of tidal harmonic constituents')
68+
ncwriteatt(file,'frequency','units','radians/second')
69+
ncwrite(file, 'frequency', f24dat.omega) ;
70+
71+
nccreate(file, 'sal_amplitude','Dimensions',{'num_constituents',ntip(1),'node',node},'DataType','single','DeflateLevel',dl,'FillValue',fillvalue) ;
72+
ncwriteatt(file,'sal_amplitude','standard_name','amplitude_of_self_attraction_and_loading_tide_elevation')
73+
ncwriteatt(file,'sal_amplitude','long_name','amplitude of self attraction and loading tide elevation')
74+
ncwriteatt(file,'sal_amplitude','units','m')
75+
ncwrite(file, 'sal_amplitude', squeeze(f24dat.Val(:,2,:))) ;
76+
77+
nccreate(file, 'sal_phase','Dimensions',{'num_constituents',ntip(1),'node',node},'DataType','single','DeflateLevel',dl,'FillValue',fillvalue) ;
78+
ncwriteatt(file,'sal_phase','long_name','phase-lag of self-attraction and loading tide elevation')
79+
ncwriteatt(file,'sal_phase','standard_name','phase_lag_of_self_attraction_and_loading_tide_elevation')
80+
ncwriteatt(file,'sal_phase','units','degrees with respect to GMT/UTC')
81+
ncwrite(file, 'sal_phase', squeeze(f24dat.Val(:,3,:))) ;
82+
83+
ncwriteatt(file,'/','title','The self-attraction and loading terms for an ADCIRC simulation');
84+
ncwriteatt(file,'/','creation_date',datestr(now));
85+
ncwriteatt(file,'/','source',"Made by OceanMesh2D writefort24");
86+
ncwriteatt(file,'/','references',"https://github.com/CHLNDDEV/OceanMesh2D/" );
87+
88+
else
89+
error(['format = ' format ' is invalid. Choose from ascii or netcdf'])
2290
end
23-
24-
fclose(fid) ;
2591
%EOF
26-
end
92+
end

‎README.md‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -150,6 +150,7 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/)
150150
### Unreleased (on current HEAD of the Projection branch)
151151
## Added
152152
- Geoid offset nodal attribute in `Calc_f13` subroutine. https://github.com/CHLNDDEV/OceanMesh2D/pull/251
153+
- Support for writing Self Attraction and Loading (SAL) files in NetCDF for the ADCIRC model. https://github.com/CHLNDDEV/OceanMesh2D/pull/231
153154
## Changed
154155
- Default mesh improvement strategy is `ds` 2.
155156
- Retrieve boundary indices in `msh.get_boundary_of_mesh` method. https://github.com/CHLNDDEV/OceanMesh2D/pull/259

‎utilities/Make_f24.m‎

Lines changed: 55 additions & 42 deletions
Original file line numberDiff line numberDiff line change
@@ -3,82 +3,91 @@
33
% Takes a msh object and interpolates the global SAL term into the f24
44
% struct
55
% Assumes that saldata is in the MATLAB path
6-
% The saldata required can be downloaded from:
7-
% saldata = 'FES2004' : Source at: ftp://ftp.legos.obs-mip.fr/pub/soa/...
8-
% maree/tide_model/global_solution/fes2004/
9-
%
6+
% The saldata required can be downloaded from:
7+
% saldata = 'FES2004' : Source at: ftp://ftp.legos.obs-mip.fr/pub/soa/...
8+
% maree/tide_model/global_solution/fes2004/load/
9+
%
1010
% saldata = 'FES2014' : Source at: ftp://ftp.legos.obs-mip.fr/pub/...
11-
% FES2012-project/data/LSA/FES2014/
12-
% by default saldata = 'FES2014'
11+
% FES2012-project/data/LSA/FES2014/
12+
% by default saldata = 'FES2014'
1313
%
1414
% plot_on - 1/true: to plot and print F24 values for checking
1515
% 0/false: no plotting by default
1616
%
17-
% Created by William Pringle. July 11 2018 updated to Make_f## style
17+
% Created by William Pringle. July 11 2018 updated to Make_f## style
18+
% Modified by Keith Roberts Jun 19, 2021 to specify degree format for FES
19+
% file
1820
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1921

2022
if isempty(obj.f15)
21-
error(['The msh object must have the f15 struct populated '...
22-
'with tidal potential information'])
23+
error(['The msh object must have the f15 struct populated '...
24+
'with tidal potential information'])
2325
end
2426

2527
if nargin < 2 || isempty(saldata)
26-
saldata = 'FES2014';
28+
saldata = 'FES2014';
2729
end
2830
if nargin < 3 || isempty(plot_on)
29-
plot_on = false;
30-
end
31-
32-
ll0 = obj.f15.slam(1) ;
33-
if ( ll0 < 0 )
34-
ll0 = ll0 + 360 ;
31+
plot_on = false;
3532
end
36-
lon0 = ll0*pi/180 ; lat0 = obj.f15.slam(2)*pi/180 ;
37-
38-
R = 6378206.4; % earth radius
3933

34+
% parameter for cpp conversion
35+
R = 6378206.4; % earth radius
36+
lon0 = obj.f15.slam(1) ; % central longitude
37+
lat0 = obj.f15.slam(2) ; % central latitude
4038
% Put in the tidal potential names
4139
obj.f24.tiponame = {obj.f15.tipotag.name};
4240

43-
ntip = length(obj.f24.tiponame) ;
41+
ntip = length(obj.f24.tiponame) ;
4442

4543
% choose tidal database file names and directories
4644
database = strtrim(upper(saldata)) ;
4745
direc = '';
4846

49-
% % Load tide grid data
47+
% % Load tide grid data
5048
if strcmp(database,'FES2004')
5149
tide_grid = [direc 'load.k1.nc'];
5250
tide_prefix = [direc 'load.'];
5351
tide_suffix = '.nc';
5452
lon = ncread(tide_grid,'lon');
5553
lat = ncread(tide_grid,'lat');
56-
elseif strcmp(database,'FES2014')
54+
% -180/180 degree format for 2004
55+
if (lon0 > 180); lon0 = lon0 - 360 ; end
56+
elseif strcmp(database,'FES2014')
5757
tide_grid = [direc 'K1_sal.nc'];
5858
tide_prefix = direc;
5959
tide_suffix = '_sal.nc';
6060
lon = ncread(tide_grid,'longitude');
6161
lat = ncread(tide_grid,'latitude');
6262
[lon,lat] = ndgrid(lon,flipud(lat));
63+
% 0-360 degree format for 2014
64+
if (lon0 < 0); lon0 = lon0 + 360 ; end
65+
else
66+
error(['database = ' database ' is invalid.'...
67+
' Choose from FES2004 or FES2014'])
6368
end
69+
% Convert CPP origin parameters to radians
70+
lon0 = lon0*pi/180 ; lat0 = lat0*pi/180 ;
6471

65-
% CPP Conversion of lat/lon
66-
lon = lon * pi/180; lat = lat * pi/180;
72+
% CPP Conversion of FES lat/lon
73+
lon = lon * pi/180; lat = lat * pi/180;
6774
x = R * (lon - lon0) * cos(lat0);
6875
y = R * lat;
6976

70-
% Get mesh vertices and change to FES 0 to 360 deg style
77+
% Doing the CPP conversion of the mesh
7178
VX = obj.p;
72-
VX(VX(:,1)<0,1) = VX(VX(:,1)<0,1) + 360;
73-
74-
% Doing the CPP conversion
75-
xx = VX(:,1) * pi/180; yy = VX(:,2) * pi/180;
79+
if strcmp(database,'FES2004')
80+
VX(VX(:,1)>180,1) = VX(VX(:,1)>180,1) - 360;
81+
elseif strcmp(database,'FES2014')
82+
VX(VX(:,1)<0,1) = VX(VX(:,1)<0,1) + 360;
83+
end
84+
xx = VX(:,1) * pi/180; yy = VX(:,2) * pi/180;
7685
xx = R * (xx - lon0) * cos(lat0);
7786
yy = R * yy;
7887

7988
% Now interpolate onto grid and put into fort.24 struct
8089
nnodes = length(VX) ;
81-
kvec = (1:nnodes)';
90+
kvec = (1:nnodes)';
8291
obj.f24.Val = zeros(ntip,3,nnodes) ;
8392
for icon = 1: ntip
8493
% Get the frequency
@@ -106,8 +115,8 @@
106115

107116
% Do the gridded Interpolation
108117
F = griddedInterpolant(x,y,z,'linear','none');
109-
Z = F(xx,yy);
110-
118+
Z = F(xx,yy);
119+
111120
% Convert back to amp and phase
112121
amp = abs(Z);
113122
phs = angle(Z)*180/pi;
@@ -119,20 +128,24 @@
119128

120129
% Plot interpolated results
121130
if plot_on
122-
figure(1); fastscatter(VX(:,1),VX(:,2),amp);
123-
colorbar;
124-
constituent = obj.f24.tiponame{icon};
125-
title(constituent)
126-
print(['F24_' constituent '_check'],'-dpng')
127-
end
128-
131+
figure(1); fastscatter(VX(:,1),VX(:,2),amp);
132+
colorbar;
133+
constituent = obj.f24.tiponame{icon};
134+
title(constituent)
135+
print(['F24_' constituent '_check'],'-dpng')
136+
end
137+
129138
% Put into the struct
130-
obj.f24.Val(icon,:,:) = [kvec'; amp'; phs'];
139+
obj.f24.Val(icon,:,:) = [kvec'; amp'; phs'];
140+
141+
if any(isnan(amp))
142+
warning('NaNs detected in amplitudes. Is degree format correct?')
143+
end
131144
end
132145

133146
if obj.f15.ntip ~= 2
134-
disp('Setting ntip = 2')
135-
obj.f15.ntip = 2;
147+
disp('Setting ntip = 2')
148+
obj.f15.ntip = 2;
136149
end
137150
%EOF
138151
end

0 commit comments

Comments
 (0)