Bug fix/interp and cleaning

@msh/msh.m

@@ -1106,6 +1106,14 @@ function plotter(cmap,round_dec,yylabel,apply_pivot)
                     obj.renum = varargin{ii + 1};
                 end
             end
+            if ~isempty(pfixV)
+                if opt.ds == 2
+                    warning(['Fixed points is not compatible with ' ...
+                        'hill-climbing smoothing (ds = 2), changing ' ...
+                        'to implicit smoothing (ds = 1)'])
+                    opt.ds = 1;
+                end
+            end
 
             % display options
             disp('INFO: The following cleaning options have been enabled..')
@@ -1335,7 +1343,8 @@ function plotter(cmap,round_dec,yylabel,apply_pivot)
             if nargin < 2
                 error('Needs type: one of "auto", "outer", "inner", "delete", or "weirs"')
             end
-            if iscell(varargin{1})
+
+            if ~isempty(varargin) && iscell(varargin{1})
                 varargin = varargin{1};
             end
             L = 1e3;
@@ -1447,7 +1456,7 @@ function plotter(cmap,round_dec,yylabel,apply_pivot)
                             else
                                % set distance to be larger than dist_lim  
                                % everywhere when no land exists
-                               ldst = eb_mid(:,1)*0 + 2*dist_lim
+                               ldst = eb_mid(:,1)*0 + 2*dist_lim;
                             end    
                             eb_class = ldst > dist_lim;
                             if strcmp(classifier,'both')

@msh/private/GridData.m

@@ -389,8 +389,7 @@
         if nanfill
             if sum(isnan(b)) > 0
                 localcoord = obj.p(K,:);
-                [KI, ~] = ourKNNsearch(localcoord(~isnan(b),:),',localcoord(isnan(b)',1); 
-                %KI = knnsearch(localcoord(~isnan(b),:),localcoord(isnan(b),:));
+                [KI, ~] = ourKNNsearch(localcoord(~isnan(b),:)',localcoord(isnan(b),:)',1); 
                 bb = b(~isnan(b));
                 b(isnan(b)) = bb(KI); clear bb localcoord
             end
@@ -420,14 +419,12 @@
             if ~isempty(find(isnan(bx),1))
                 localcoord = obj.p(K,:);
                 [KI, ~] = ourKNNsearch(localcoord(~isnan(bx),:)',localcoord(isnan(bx),:)',1); 
-                %KI = knnsearch(localcoord(~isnan(bx),:),localcoord(isnan(bx),:));
                 bb = bx(~isnan(bx));
                 bx(isnan(bx)) = bb(KI); clear bb localcoord
             end
             if ~isempty(find(isnan(by),1))
                 localcoord = obj.p(K,:);
                 [KI, ~] = ourKNNsearch(localcoord(~isnan(by),:)',localcoord(isnan(by),:)',1); 
-                %KI = knnsearch(localcoord(~isnan(by),:),localcoord(isnan(by),:));
                 bb = by(~isnan(by));
                 by(isnan(by)) = bb(KI); clear bb localcoord
             end

@msh/private/writefort15.m

@@ -88,17 +88,32 @@
 fprintf( fid, '%g  \t \t ! REFTIM \n', f15dat.reftim ) ; 
 
 % WTIMINC
-if f15dat.nws == 8
-  fprintf( fid, '%d %d %d %d %d %g', f15dat.wtimnc ) ;
-  fprintf( fid, '  \t ! YYYY MM DD HH24 StormNumber BLAdj \n' ) ;  
-elseif f15dat.nws >= 19
-  fprintf( fid, '%d %d %d %d %d %g %d', f15dat.wtimnc ) ;
-  fprintf( fid, '  \t ! YYYY MM DD HH24 StormNumber BLAdj geofactor \n' ) ;  
-elseif f15dat.nws > 0
-  fprintf( fid, '%d  ', f15dat.wtimnc ) ;
-  fprintf( fid, '  \t ! WTMINC \n' ) ;       
-end
-
+if f15dat.nws > 0
+   nws = f15dat.nws;
+   wtimnc = f15dat.wtimnc;
+   rstimnc_fmt = ''; rstimnc_msg = '';
+   if nws > 300
+      wtimnc = [f15dat.wtimnc f15dat.rstimnc];
+      nws = f15dat.nws - 300;
+      rstimnc_fmt = ' %d'; 
+      rstimnc_msg = ' RSTIMINC';
+   end
+   if nws == 8
+      wtimnc_fmt = '%d %d %d %d %d %g';
+      wtimnc_msg = 'YYYY MM DD HH24 StormNumber BLAdj';
+   elseif nws >= 19
+      wtimnc_fmt = '%d %d %d %d %d %g %d';
+      wtimnc_msg = 'YYYY MM DD HH24 StormNumber BLAdj geofactor';
+   else
+      wtimnc_fmt = '%d  ';
+      wtimnc_msg = 'WTIMINC';
+   end
+   wtimnc_fmt = [wtimnc_fmt rstimnc_fmt]; 
+   wtimnc_msg = [wtimnc_msg rstimnc_msg]; 
+   fprintf( fid, wtimnc_fmt, wtimnc ) ;
+   fprintf( fid, ['  \t ! ' wtimnc_msg ' \n'] ) ; 
+end  
+
 % RNDY
 fprintf( fid, '%g   \t \t ! RNDY \n', f15dat.rndy ) ;
 

README.md

@@ -157,18 +157,27 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/)
 
 ### Unreleased (on current HEAD of the Projection branch)
 ## Added
+- `namelist` and `RSTIMNC` input arguments for `Make_f15.m` fort.15 generator. updated the help message for all input argumebts to `Make_f15`. https://github.com/CHLNDDEV/OceanMesh2D/pull/283 
+- New stability namelist options to `Make_f15.m` fort.15 generator. https://github.com/CHLNDDEV/OceanMesh2D/pull/283 
 - Optionality for mesh2dgen to choose the method of mesh generation `kind`, and the maximum iteration count `iter`. https://github.com/CHLNDDEV/OceanMesh2D/pull/272
 - Option `improve_with_reduced_quality` to `meshgen` for allowing for mesh improvements even when quality is decreasing or large number of nodes eliminated, which sometimes is necessary to force the advancement in mesh quality.
 - Option `delaunay_elim_on_exit` to `meshgen` to skip the last call to `delaunay_elim` to potentially avoid deleting boundary elements.
 - 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
+- removed namelists from default setup in `Make_f15.m` fort.15 generator to be invoked by user as an input argument instead. https://github.com/CHLNDDEV/OceanMesh2D/pull/283 
+- Use implicit smoother (`ds=1`) in `msh.clean` when fix points are present. https://github.com/CHLNDDEV/OceanMesh2D/pull/283
 - Default filename for the dynamicWaterLevelCorrection is now `null` so that it is not evoked by default. https://github.com/CHLNDDEV/OceanMesh2D/pull/272
 - Default mesh improvement strategy is `ds` 2.
 - Retrieve boundary indices in `msh.get_boundary_of_mesh` method. https://github.com/CHLNDDEV/OceanMesh2D/pull/259
 - `msh.offset63` struct and associated write/make routines for dynamicwaterlevel offset functionality. https://github.com/CHLNDDEV/OceanMesh2D/pull/259
 - dynamicWaterLevelCorrection to fort.15 namelist, and PRBCKGRND option to met fort.15 namelist. https://github.com/CHLNDDEV/OceanMesh2D/pull/261
 ## Fixed
+- writing out the wave coupling timestep `RSTIMNC` on the `WTIMNC` line of fort.15 when `NWS > 300`. https://github.com/CHLNDDEV/OceanMesh2D/pull/283 
+- `make_bc` call with empty varargin, e.g., when calling inner. https://github.com/CHLNDDEV/OceanMesh2D/pull/283
+- `Make_offset63` call with constant offset (length 1 `time_vector`). https://github.com/CHLNDDEV/OceanMesh2D/pull/283
+- `ourKNNsearch` call in `nanfill` option of `Griddata` (`msh.interp`). https://github.com/CHLNDDEV/OceanMesh2D/pull/283
+- `m_map` link in `setup.sh` Issue #277. https://github.com/CHLNDDEV/OceanMesh2D/pull/283
 - Updated `Calc_f13.m` to avoid an "Unrecognized variable" error by ensuring "broken" is always defined. https://github.com/CHLNDDEV/OceanMesh2D/pull/282
 - Fixed test for likely geographic coordinates in `Make_f15.m`. https://github.com/CHLNDDEV/OceanMesh2D/pull/282
 - updated `Gridded_to_Mesh_SeaBed_DepthAveraged.m` to fix the infinite loop in using `Cal_IT_Fric.m` by filling in the NaNs at greater depths with values from above. https://github.com/CHLNDDEV/OceanMesh2D/pull/280

setup.sh

@@ -17,7 +17,7 @@ if $m_map; then
      echo "m_map directory already exists"
   else
      # download m_map archive and unzip
-     wget "http://www.eos.ubc.ca/%7Erich/m_map1.4.zip" -O m_map.zip
+     wget "https://www.eoas.ubc.ca/~rich/m_map1.4.zip" -O m_map.zip
      unzip m_map.zip
      rm m_map.zip
   fi

utilities/Make_f15.m

@@ -43,6 +43,24 @@
 %               ...('sta_database',{'CO-OPS',1,'NDBC',2,'NDBC',3}) to output 
 %               the CO-OPS stations for elevation and NDBC stations for velocity and 
 %               meteorological recording
+% 
+%               'tau0' : specify the GWCE control parameter: https://wiki.adcirc.org/TAU0 
+%
+%               'tau0minmax' : specify the min/max of GWCE control parameter 
+%                (must be supplied when tau0 == -5) : https://wiki.adcirc.org/TAU0 
+%
+%               'NWS' : specify the wind control parameter: https://wiki.adcirc.org/NWS 
+%
+%               'WTIMNC' : specify the wind timestep (must be supplied of NWS > 0): 
+%                https://wiki.adcirc.org/WTIMINC
+%                https://wiki.adcirc.org/Supplemental_meteorological/wave/ice_parameters 
+%
+%               'RSTIMNC' : specify the SWAN wave timestep (must be supplied of NWS > 300): 
+%                https://wiki.adcirc.org/Supplemental_meteorological/wave/ice_parameters
+%    
+%               'namelist' : cell/string list of namelists to add from the following:
+%                {'met', 'dynamicWaterLevelCorrection', 'limit', 'velwd'}
+%                https://wiki.adcirc.org/Fort.15_file_format#Namelists
 %
 %  Outputs:  1) msh class obj with f15 struct populated
 %
@@ -94,6 +112,7 @@
     % NWS
     f15dat.nws = 0 ;
     f15dat.wtimnc = 0;
+    f15dat.rstimnc = 0;
     % NRAMP
     f15dat.nramp = 0 ;
     % G
@@ -185,34 +204,6 @@
     f15dat.extraline(9).msg = 'name@instit.edu';
     f15dat.extraline(10).msg = '';
 
-    % control lists
-    % met
-    f15dat.controllist(1).type = 'met';
-    f15dat.controllist(1).var(1).name = 'WindDragLimit';
-    f15dat.controllist(1).var(1).val = 0.0025;
-    f15dat.controllist(1).var(2).name = 'PRBCKGRND';
-    f15dat.controllist(1).var(2).val = 1013;
-    f15dat.controllist(1).var(3).name = 'DragLawString';
-    f15dat.controllist(1).var(3).val = 'default'; 
-    f15dat.controllist(1).var(4).name = 'outputWindDrag';
-    f15dat.controllist(1).var(4).val = 'F'; 
-    f15dat.controllist(1).var(5).name = 'invertedBarometerOnElevationBoundary';
-    f15dat.controllist(1).var(5).val = 'F'; 
-    % dynamicwaterlevelcorrection control
-    f15dat.controllist(2).type = 'dynamicWaterLevelCorrection';
-    f15dat.controllist(2).var(1).name = [f15dat.controllist(2).type 'FileName'];
-    f15dat.controllist(2).var(1).val = 'null';
-    f15dat.controllist(2).var(2).name = [f15dat.controllist(2).type 'Multiplier'];
-    f15dat.controllist(2).var(2).val = 1.0; 
-    f15dat.controllist(2).var(3).name = [f15dat.controllist(2).type 'RampStart'];
-    f15dat.controllist(2).var(3).val = 0.0; 
-    f15dat.controllist(2).var(4).name = [f15dat.controllist(2).type 'RampEnd'];
-    f15dat.controllist(2).var(4).val = 0.0; 
-    f15dat.controllist(2).var(5).name = [f15dat.controllist(2).type 'RampReferenceTime'];
-    f15dat.controllist(2).var(5).val = 'coldstart'; 
-    f15dat.controllist(2).var(6).name = [f15dat.controllist(2).type 'SkipSnaps'];
-    f15dat.controllist(2).var(6).val = 0; 
-
     % Put into the msh class
     obj.f15 = f15dat;
 end
@@ -253,9 +244,10 @@
 const = [];
 tidal_database = [];
 sta_database = [];
+namelists = [];
 if ~isempty(varargin)
     names = {'const','tidal_database','NWS','WTIMNC','tau0','tau0minmax',...
-             'sta_database'};
+             'sta_database','RSTIMNC','namelist'};
     for ii = 1:length(names)
         ind = find(~cellfun(@isempty,strfind(varargin(1:2:end),names{ii})));
         if ~isempty(ind)
@@ -273,6 +265,10 @@
                 obj.f15.tau0minmax = varargin{ind*2};
             elseif ii == 7
                 sta_database = varargin{ind*2};
+            elseif ii == 8
+                obj.f15.rstimnc = varargin{ind*2};
+            elseif ii == 9
+                namelists = varargin{ind*2};
             end
         end    
     end
@@ -284,6 +280,10 @@
     if obj.f15.wtimnc == 0
         error(['Must specify WTIMNC option with NWS = ' num2str(obj.f15.nws)])
     end
+
+    if obj.f15.nws > 300 && obj.f15.rstimnc == 0
+        error(['Must specify RSTIMNC option with NWS = ' num2str(obj.f15.nws)])
+    end
 end
 
 % Tau0FullDomainMin, Tau0FullDomainMax 
@@ -356,15 +356,58 @@
     obj = tidal_stations_parser(obj,string(sta_database(1:2:end-1)),...
                                 sta_database(2:2:end));
 end
-% 
-% % NOUTGC
-% if ( f15dat.im == 10 )
-%     f15dat.outgc = readlinevec( fid ) ; 
-% end
-% 
-% % NOUTGW
-% if ( f15dat.nws ~= 0 ) 
-%     f15dat.outgw = readlinevec( fid ) ;
-% end
+
+% adding the namelist controls
+if ~isempty(namelists)
+    ci = 0;
+    if find(contains(namelists,'met'),1)
+       % met
+       ci = ci + 1;
+       obj.f15.controllist(ci).type = 'met';
+       obj.f15.controllist(ci).var(1).name = 'WindDragLimit';
+       obj.f15.controllist(ci).var(1).val = 0.0025;
+       obj.f15.controllist(ci).var(2).name = 'PRBCKGRND';
+       obj.f15.controllist(ci).var(2).val = 1013;
+       obj.f15.controllist(ci).var(3).name = 'DragLawString';
+       obj.f15.controllist(ci).var(3).val = 'default'; 
+       obj.f15.controllist(ci).var(4).name = 'outputWindDrag';
+       obj.f15.controllist(ci).var(4).val = 'F'; 
+       obj.f15.controllist(ci).var(5).name = 'invertedBarometerOnElevationBoundary';
+       obj.f15.controllist(ci).var(5).val = 'F'; 
+    end
+    if find(contains(namelists,'dynamicWaterLevelCorrection','IgnoreCase',true),1)
+       % dynamicwaterlevelcorrection control
+       ci = ci + 1;
+       obj.f15.controllist(ci).type = 'dynamicWaterLevelCorrection';
+       obj.f15.controllist(ci).var(1).name = [obj.f15.controllist(ci).type 'FileName'];
+       obj.f15.controllist(ci).var(1).val = 'null';
+       obj.f15.controllist(ci).var(2).name = [obj.f15.controllist(ci).type 'Multiplier'];
+       obj.f15.controllist(ci).var(2).val = 1.0; 
+       obj.f15.controllist(ci).var(3).name = [obj.f15.controllist(ci).type 'RampStart'];
+       obj.f15.controllist(ci).var(3).val = 0.0; 
+       obj.f15.controllist(ci).var(4).name = [obj.f15.controllist(ci).type 'RampEnd'];
+       obj.f15.controllist(ci).var(4).val = 0.0; 
+       obj.f15.controllist(ci).var(5).name = [obj.f15.controllist(ci).type 'RampReferenceTime'];
+       obj.f15.controllist(ci).var(5).val = 'coldstart'; 
+       obj.f15.controllist(ci).var(6).name = [obj.f15.controllist(ci).type 'SkipSnaps'];
+       obj.f15.controllist(ci).var(6).val = 0;
+    end 
+    if find(contains(namelists,'limit'),1)
+       ci = ci + 1;
+       % limit control
+       obj.f15.controllist(ci).type = 'limit';
+       obj.f15.controllist(ci).var(1).name = 'slim';
+       obj.f15.controllist(ci).var(1).val = 4e-4;
+       obj.f15.controllist(ci).var(2).name = 'windlim';
+       obj.f15.controllist(ci).var(2).val = 'T';
+    end
+    if find(contains(namelists,'velwd'),1)
+       ci = ci + 1;
+       % wetdry velocity control
+       obj.f15.controllist(ci).type = 'velwd';
+       obj.f15.controllist(ci).var(1).name = 'directvelWD';
+       obj.f15.controllist(ci).var(1).val = 'T';
+    end
+end
 
 end

utilities/Make_offset63.m

@@ -32,6 +32,9 @@
                               seconds(time_vector(end) - time_vector(1))/ ...
                               max(length(time_vector)-1,1)...
                              );
+if isnan(obj.offset63.time_interval)
+    obj.offset63.time_interval = 999*3600*24;
+end 
 obj.offset63.default_val = 0;
 obj.offset63.offset_nodes = offset_nodes;
 obj.offset63.offset_values = offset_values;