diff --git a/@msh/msh.m b/@msh/msh.m
index 0d179dc8..8a431eae 100644
--- a/@msh/msh.m
+++ b/@msh/msh.m
@@ -176,9 +176,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
@@ -187,6 +187,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
@@ -212,7 +213,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'], obj.p, varargin);
                 end
                 if ~isempty(obj.f5354)
                     writefort5354( obj.f5354, fname );
@@ -268,7 +269,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'], obj.p, varargin);
                 end
                 if any(contains(type,'5354')) && ~isempty(obj.f5354)
                     writefort5354( obj.f5354, fname );
diff --git a/@msh/private/writefort24.m b/@msh/private/writefort24.m
index 0368bd32..4e81035c 100644
--- a/@msh/private/writefort24.m
+++ b/@msh/private/writefort24.m
@@ -1,26 +1,92 @@
-function fid = writefort24( f24dat, finame )
+function fid = writefort24( f24dat, finame, points, 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';
+   warning('Using ASCII file format. Pass ''netcdf'' to write in NetCDF')
+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;
+    
+    node = length(points); 
+    tipnames = char(f24dat.tiponame); ntip = size(tipnames) ;
+    fillvalue = 0.0;
+        
+    nccreate(file,'x','Dimensions',{'node',node},'DataType','int32','DeflateLevel',dl) ;
+    ncwriteatt(file,'x','standard_name', 'longitude') ;
+    ncwriteatt(file,'x','units', 'degrees_east') ;
+    ncwriteatt(file,'x','positive', 'east') ;
+    ncwrite(file, 'x', points(:,1)) ;
 
-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,:,:));
+    nccreate(file,'y','Dimensions',{'node',node},'DataType','int32','DeflateLevel',dl) ;
+    ncwriteatt(file,'y','standard_name', 'latitude') ;
+    ncwriteatt(file,'y','units', 'degrees_north') ;
+    ncwriteatt(file,'y','positive', 'north') ;
+    ncwrite(file, 'y', points(:,2)) ;    
+    
+    nccreate(file,'constituents','Dimensions',{'num_constituents',ntip(1),'char_len',ntip(2)},'DataType','char') ;
+    ncwriteatt(file,'constituents','standard_name', 'name_of_tidal_harmonic_constituents') ;
+    ncwriteatt(file,'constituents','long_name', 'name of tidal harmonic constituents') ;
+    ncwrite(file, 'constituents', tipnames) ;
+    
+    nccreate(file,'frequency','dimensions',{'num_constituents',ntip(1)});
+    ncwriteatt(file,'frequency','standard_name','frequency_of_tidal_harmonic_constituents')
+    ncwriteatt(file,'frequency','long_name','frequency of tidal harmonic constituents')
+    ncwriteatt(file,'frequency','units','radians/second')
+    ncwrite(file, 'frequency', f24dat.omega) ;
+    
+    nccreate(file, 'sal_amplitude','Dimensions',{'num_constituents',ntip(1),'node',node},'DataType','single','DeflateLevel',dl,'FillValue',fillvalue) ;
+    ncwriteatt(file,'sal_amplitude','standard_name','amplitude_of_self_attraction_and_loading_tide_elevation')
+    ncwriteatt(file,'sal_amplitude','long_name','amplitude of self attraction and loading tide elevation')
+    ncwriteatt(file,'sal_amplitude','units','m')
+    ncwrite(file, 'sal_amplitude', squeeze(f24dat.Val(:,2,:))) ;
+    
+    nccreate(file, 'sal_phase','Dimensions',{'num_constituents',ntip(1),'node',node},'DataType','single','DeflateLevel',dl,'FillValue',fillvalue) ;
+    ncwriteatt(file,'sal_phase','long_name','phase-lag of self-attraction and loading tide elevation')
+    ncwriteatt(file,'sal_phase','standard_name','phase_lag_of_self_attraction_and_loading_tide_elevation')
+    ncwriteatt(file,'sal_phase','units','degrees with respect to GMT/UTC')
+    ncwrite(file, 'sal_phase', squeeze(f24dat.Val(:,3,:))) ;
+   
+    ncwriteatt(file,'/','title','The self-attraction and loading terms for an ADCIRC simulation');
+    ncwriteatt(file,'/','creation_date',datestr(now));
+    ncwriteatt(file,'/','source',"Made by OceanMesh2D writefort24");
+    ncwriteatt(file,'/','references',"https://github.com/CHLNDDEV/OceanMesh2D/" );
+
+else
+    error(['format = ' format ' is invalid. Choose from ascii or netcdf'])
 end
-  
-fclose(fid) ; 
 %EOF
-end
\ No newline at end of file
+end
diff --git a/README.md b/README.md
index 9e2abae2..d02d1e37 100644
--- a/README.md
+++ b/README.md
@@ -150,6 +150,7 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/)
 ### Unreleased (on current HEAD of the Projection branch)
 ## Added
 - Geoid offset nodal attribute in `Calc_f13` subroutine. https://github.com/CHLNDDEV/OceanMesh2D/pull/251
+- Support for writing Self Attraction and Loading (SAL) files in NetCDF for the ADCIRC model. https://github.com/CHLNDDEV/OceanMesh2D/pull/231
 ## Changed
 - Default mesh improvement strategy is `ds` 2.
 - Retrieve boundary indices in `msh.get_boundary_of_mesh` method. https://github.com/CHLNDDEV/OceanMesh2D/pull/259
diff --git a/utilities/Make_f24.m b/utilities/Make_f24.m
index 2a93e294..997f8279 100644
--- a/utilities/Make_f24.m
+++ b/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