Support for writing Self Attraction and Loading (SAL) files in NetCDF for the ADCIRC model

@msh/msh.m

@@ -175,9 +175,9 @@
         end
 
         % write mesh to disk
-        function write(obj,fname,type)
+        function write(obj,fname,type,varargin)
             % Usage:
-            % write(obj,fname,type)
+            % write(obj,fname,type,varargin)
             %
             % Examples:
             % write(obj);       % writes all available data to fort_1.xx (ADCIRC) files
@@ -186,6 +186,7 @@ function write(obj,fname,type)
             % write(obj,fname,'gr3'); % writes mesh data to fname.gr3 (SCHISM) file
             % write(obj,fname,'ww3'); % writes mesh data to fname.ww3 (WaveWatchIII) file
             % write(obj,fname,{'13','14'}); % writes mesh data and f13 attribute data to fname.14 and fname.13 (ADCIRC) files
+            % write(obj,fname,'24','netcdf'); % writes fort.24 SAL data to fname.24.nc netcdf file
             if nargin == 1
                 fname = 'fort_1';
             end
@@ -211,7 +212,7 @@ function write(obj,fname,type)
                     writefort15( obj.f15, [fname '.15'], obj.bd );
                 end
                 if ~isempty(obj.f24)
-                    writefort24( obj.f24, [fname '.24'] );
+                    writefort24( obj.f24, [fname '.24'], varargin);
                 end
                 if ~isempty(obj.f5354)
                     writefort5354( obj.f5354, fname );
@@ -264,7 +265,7 @@ function write(obj,fname,type)
                     writefort19( obj.f2001, [fname '.2001'] );
                 end
                 if any(contains(type,'24')) && ~isempty(obj.f24)
-                    writefort24( obj.f24, [fname '.24'] );
+                    writefort24( obj.f24, [fname '.24'], varargin);
                 end
                 if any(contains(type,'5354')) && ~isempty(obj.f5354)
                     writefort5354( obj.f5354, fname );

@msh/private/writefort24.m

@@ -1,26 +1,73 @@
-function fid = writefort24( f24dat, finame )
+function fid = writefort24( f24dat, finame, format )
 %
 %
-if ( nargin == 1 ) 
-	finputname = 'fort.24_1' ;
+if ( nargin == 1 )
+    finputname = 'fort.24_1' ;
 else
-	finputname = finame ;   
+    finputname = finame ;
 end
 
-fid = fopen(finputname,'w') ;
+if nargin < 3 || isempty(format) 
+   format = 'ascii';
+end
+if iscell(format)
+   format = format{1};
+end
+
+if strcmp(format,'ascii') % LEGACY
+    fid = fopen(finputname,'w') ;
+
+    tipnames = f24dat.tiponame; ntip = length(tipnames) ;
+    for icon = 1: ntip
+        tipname = tipnames{icon};
+        fprintf('Writing SAL %s data \n', char(tipname)) ;
+        % The constituent details
+        fprintf(fid,'%s \n',[char(tipname) ' SAL']) ;
+        fprintf(fid,'%17.15f \n',f24dat.omega(icon)) ;
+        fprintf(fid,'%d \n',1) ;
+        fprintf(fid,'%s \n',char(tipname)) ;
+        fprintf(fid,'%d \t %12.6f  %12.6f \n',f24dat.Val(icon,:,:));
+    end
+
+    fclose(fid) ;
+
+elseif strcmp(format,'netcdf') % NETCDF
+    % create nc file or overwrite it
+    file = [finputname,'.nc'];
+    if exist(file)
+        delete(file)
+    end
+    % deflate level (set zero for no deflation if necessary.. using single precision so not very different in size)
+    dl = 5;
+    % SAL indices
+    indices = squeeze(f24dat.Val(1,1,:)); 
+    sal_nodes = length(indices);
+    nccreate(file,'mesh_indices','Dimensions',{'sal_nodes',sal_nodes},'DataType','int32','DeflateLevel',dl) ;
+    ncwriteatt(file,'mesh_indices','long_name','index of SAL values into mesh coordinates')
+    ncwrite(file, 'mesh_indices', indices) ;
+    % SAL name
+    tipnames = char(f24dat.tiponame); ntip = size(tipnames) ;
+    nccreate(file,'const','dimensions',{'num_const',ntip(1),'char_len',ntip(2)},'DataType','char') ;
+    ncwriteatt(file,'const','long_name','names of the tidal harmonic constituents')
+    ncwrite(file, 'const', tipnames) ;
+    % SAL frequencies
+    nccreate(file,'frequency','dimensions',{'num_const',ntip(1)});
+    ncwriteatt(file,'frequency','long_name','frequency of harmonic constituents')
+    ncwriteatt(file,'frequency','units','rad/s')
+    ncwrite(file, 'frequency', f24dat.omega) ;
+    % SAL amplitudes
+    nccreate(file, 'salamp','Dimensions',{'num_const',ntip(1),'sal_nodes',sal_nodes},'DataType','single','DeflateLevel',dl) ;
+    ncwriteatt(file,'salamp','long_name','amplitude of self-attraction and loading tide elevation')
+    ncwriteatt(file,'salamp','units','m')
+    ncwrite(file, 'salamp', squeeze(f24dat.Val(:,2,:))) ;
+    % SAL phase lags
+    nccreate(file, 'salphs','Dimensions',{'num_const',ntip(1),'sal_nodes',sal_nodes},'DataType','single','DeflateLevel',dl) ;
+    ncwriteatt(file,'salphs','long_name','phase-lag of self-attraction and loading tide elevation')
+    ncwriteatt(file,'salphs','units','degrees with respect to GMT/UTC')
+    ncwrite(file, 'salphs', squeeze(f24dat.Val(:,3,:))) ;
 
-tipnames = f24dat.tiponame; ntip = length(tipnames) ; 
-for icon = 1: ntip
-    tipname = tipnames{icon};
-    fprintf('Writing SAL %s data \n', char(tipname)) ; 
-    % The constituent details
-    fprintf(fid,'%s \n',[char(tipname) ' SAL']) ;
-    fprintf(fid,'%17.15f \n',f24dat.omega(icon)) ;
-    fprintf(fid,'%d \n',1) ;  
-    fprintf(fid,'%s \n',char(tipname)) ;
-    fprintf(fid,'%d \t %12.6f  %12.6f \n',f24dat.Val(icon,:,:));
+else
+    error(['format = ' format ' is invalid. Choose from ascii or netcdf'])
 end
-
-fclose(fid) ; 
 %EOF
-end
+end

utilities/Make_f24.m

@@ -3,82 +3,91 @@
 % Takes a msh object and interpolates the global SAL term into the f24
 % struct
 % Assumes that saldata is in the MATLAB path
-% The saldata required can be downloaded from:                                                          
-% saldata = 'FES2004' : Source at: ftp://ftp.legos.obs-mip.fr/pub/soa/...  
-%                                 maree/tide_model/global_solution/fes2004/           
-%                                                                         
+% The saldata required can be downloaded from:
+% saldata = 'FES2004' : Source at: ftp://ftp.legos.obs-mip.fr/pub/soa/...
+%                                 maree/tide_model/global_solution/fes2004/load/
+%
 % saldata = 'FES2014' : Source at: ftp://ftp.legos.obs-mip.fr/pub/...
-%                                  FES2012-project/data/LSA/FES2014/             
-% by default saldata = 'FES2014' 
+%                                  FES2012-project/data/LSA/FES2014/
+% by default saldata = 'FES2014'
 %
 % plot_on -  1/true: to plot and print F24 values for checking
 %            0/false: no plotting by default
 %
-% Created by William Pringle. July 11 2018 updated to Make_f## style             
+% Created by William Pringle. July 11 2018 updated to Make_f## style
+% Modified by Keith Roberts Jun 19, 2021 to specify degree format for FES
+% file
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 if isempty(obj.f15)
-   error(['The msh object must have the f15 struct populated '...
-          'with tidal potential information']) 
+    error(['The msh object must have the f15 struct populated '...
+        'with tidal potential information'])
 end
 
 if nargin < 2 || isempty(saldata)
-   saldata = 'FES2014'; 
+    saldata = 'FES2014';
 end
 if nargin < 3 || isempty(plot_on)
-   plot_on = false; 
-end
-
-ll0 = obj.f15.slam(1) ;
-if ( ll0 < 0 ) 
-    ll0 = ll0 + 360 ; 
+    plot_on = false;
 end
-lon0 = ll0*pi/180 ; lat0 = obj.f15.slam(2)*pi/180 ; 
-
-R = 6378206.4; % earth radius  
 
+% parameter for cpp conversion
+R = 6378206.4; % earth radius
+lon0 = obj.f15.slam(1) ; % central longitude
+lat0 = obj.f15.slam(2) ; % central latitude
 % Put in the tidal potential names
 obj.f24.tiponame = {obj.f15.tipotag.name};
 
-ntip = length(obj.f24.tiponame) ;  
+ntip = length(obj.f24.tiponame) ;
 
 % choose tidal database file names and directories
 database = strtrim(upper(saldata)) ;
 direc    = '';
 
-% % Load tide grid data 
+% % Load tide grid data
 if strcmp(database,'FES2004')
     tide_grid     = [direc 'load.k1.nc'];
     tide_prefix   = [direc 'load.'];
     tide_suffix   = '.nc';
     lon = ncread(tide_grid,'lon');
     lat = ncread(tide_grid,'lat');
-elseif  strcmp(database,'FES2014')
+    % -180/180 degree format for 2004
+    if (lon0 > 180); lon0 = lon0 - 360 ; end
+elseif strcmp(database,'FES2014')
     tide_grid     = [direc  'K1_sal.nc'];
     tide_prefix   = direc;
     tide_suffix   = '_sal.nc';
     lon = ncread(tide_grid,'longitude');
     lat = ncread(tide_grid,'latitude');
     [lon,lat] = ndgrid(lon,flipud(lat));
+    % 0-360 degree format for 2014
+    if (lon0 < 0); lon0 = lon0 + 360 ; end
+else
+    error(['database = ' database ' is invalid.'...
+          ' Choose from FES2004 or FES2014'])
 end
+% Convert CPP origin parameters to radians
+lon0 = lon0*pi/180 ; lat0 = lat0*pi/180 ;
 
-% CPP Conversion of lat/lon
-lon = lon * pi/180; lat = lat * pi/180; 
+% CPP Conversion of FES lat/lon
+lon = lon * pi/180; lat = lat * pi/180;
 x = R * (lon - lon0) * cos(lat0);
 y = R * lat;
 
-% Get mesh vertices and change to FES 0 to 360 deg style
+% Doing the CPP conversion of the mesh
 VX = obj.p;
-VX(VX(:,1)<0,1) = VX(VX(:,1)<0,1) + 360; 
-
-% Doing the CPP conversion
-xx = VX(:,1) * pi/180; yy = VX(:,2) * pi/180; 
+if strcmp(database,'FES2004')
+   VX(VX(:,1)>180,1) = VX(VX(:,1)>180,1) - 360;
+elseif strcmp(database,'FES2014')
+   VX(VX(:,1)<0,1) = VX(VX(:,1)<0,1) + 360;
+end
+xx = VX(:,1) * pi/180; yy = VX(:,2) * pi/180;
 xx = R * (xx - lon0) * cos(lat0);
 yy = R * yy;
 
 % Now interpolate onto grid and put into fort.24 struct
 nnodes = length(VX) ;
-kvec = (1:nnodes)'; 
+kvec = (1:nnodes)';
 obj.f24.Val = zeros(ntip,3,nnodes) ;
 for icon = 1: ntip
     % Get the frequency
@@ -106,8 +115,8 @@
 
     % Do the gridded Interpolation
     F = griddedInterpolant(x,y,z,'linear','none');
-    Z = F(xx,yy);  
-
+    Z = F(xx,yy);
+    
     % Convert back to amp and phase
     amp = abs(Z);
     phs = angle(Z)*180/pi;
@@ -119,20 +128,24 @@
 
     % Plot interpolated results
     if plot_on
-       figure(1); fastscatter(VX(:,1),VX(:,2),amp); 
-       colorbar;
-       constituent = obj.f24.tiponame{icon};
-       title(constituent)
-       print(['F24_' constituent '_check'],'-dpng')
-    end    
-
+        figure(1); fastscatter(VX(:,1),VX(:,2),amp);
+        colorbar;
+        constituent = obj.f24.tiponame{icon};
+        title(constituent)
+        print(['F24_' constituent '_check'],'-dpng')
+    end
+    
     % Put into the struct
-    obj.f24.Val(icon,:,:) = [kvec'; amp'; phs']; 
+    obj.f24.Val(icon,:,:) = [kvec'; amp'; phs'];
+
+    if any(isnan(amp))
+        warning('NaNs detected in amplitudes. Is degree format correct?')
+    end
 end
 
 if obj.f15.ntip ~= 2
-   disp('Setting ntip = 2')
-   obj.f15.ntip = 2;
+    disp('Setting ntip = 2')
+    obj.f15.ntip = 2;
 end
 %EOF
 end