diff --git a/@edgefx/edgefx.m b/@edgefx/edgefx.m
index 16aabb8f..e9452db5 100644
--- a/@edgefx/edgefx.m
+++ b/@edgefx/edgefx.m
@@ -399,6 +399,11 @@
             % interpolate DEM's bathy linearly onto our edgefunction grid.
             [xg,yg] = CreateStructGrid(obj);
             tmpz    = feat.Fb(xg,yg);
+            if ~isempty(feat.Fb2)
+               tmpz2 = feat.Fb2(xg,yg); 
+               tmpz(isnan(tmpz)) = tmpz2(isnan(tmpz));
+               clear tmpz2
+            end
             % Initialise wld
             obj.wld = NaN([obj.nx,obj.ny]);
             grav    = 9.807;
@@ -433,6 +438,11 @@
             [xg,yg] = CreateStructGrid(obj);
             
             tmpz    = feat.Fb(xg,yg);
+            if ~isempty(feat.Fb2)
+               tmpz2 = feat.Fb2(xg,yg); 
+               tmpz(isnan(tmpz)) = tmpz2(isnan(tmpz));
+               clear tmpz2
+            end
             tmpz(tmpz > 50) = 50; % ensure no larger than 50 m above land
             % use a harvestine assumption
             dx = obj.h0*cosd(min(yg(1,:),85)); % for gradient function
@@ -465,12 +475,12 @@
             else
                 filtit = 0;
                 for lambda = obj.fl'
-                    if all(lambda ~= 0)
-                        % do a bandpass filter
-                        tmpz_ft  = filt2(tmpz,dy,lambda,'bp') ;
-                    elseif lambda(2) == 0
+                    if length(lambda) == 1 || lambda(2) == 0
                         % do a low pass filter
                         tmpz_ft  = filt2(tmpz,dy,lambda(1),'lp') ;
+                    elseif all(lambda ~= 0)
+                        % do a bandpass filter
+                        tmpz_ft  = filt2(tmpz,dy,lambda,'bp') ;
                     else
                         % highpass filter not recommended
                         warning(['Highpass filter on bathymetry in slope' ...
@@ -636,7 +646,7 @@
                 tempbb{jj} = feat.boubox;
             end
             
-            [xg,yg]=CreateStructGrid(obj);
+            [xg,yg] = CreateStructGrid(obj);
             
             % STEP 2: For each channel point, set the resolution around a
             % neighborhood of points as |h|/ch where ch is non-dim # between 0.5-1.5
@@ -670,7 +680,13 @@
                 end
             end
             toc 
-            dp = max(1,-feat.Fb(xg(tidx),yg(tidx)));
+            tmpz = feat.Fb(xg(tidx),yg(tidx));
+            if ~isempty(feat.Fb2)
+               tmpz2 = feat.Fb2(xg(tidx),yg(tidx)); 
+               tmpz(isnan(tmpz)) = tmpz2(isnan(tmpz));
+               clear tmpz2
+            end
+            dp = max(1,-tmpz);
             obj.chd(tidx) = dp/obj.ch;
             obj.chd(obj.chd < obj.min_el_ch) = obj.min_el_ch;
             clearvars Fb pts dp radii tempbb xg yg
@@ -804,6 +820,17 @@ function plot(obj,type)
                 hh_m(nearshore & hh_m > obj.max_el_ns) = obj.max_el_ns;
             end
             hh_m(hh_m < obj.h0 ) = obj.h0;
+            
+            % Compute tmpz if necessary
+            if length(obj.max_el) > 1 || length(obj.g) > 1 || obj.dt(1) >= 0
+                tmpz    = feat.Fb(xg,yg);
+                if ~isempty(feat.Fb2)
+                   tmpz2 = feat.Fb2(xg,yg); 
+                   tmpz(isnan(tmpz)) = tmpz2(isnan(tmpz));
+                   clear tmpz2
+                end
+            end
+            
             for param = obj.max_el'
                 if numel(param)==1 && param~=0
                     mx   = obj.max_el(1);
@@ -815,8 +842,8 @@ function plot(obj,type)
                     dp1 = param(2);
                     dp2 = param(3);
                     
-                    limidx = (feat.Fb(xg,yg) < dp1 & ...
-                        feat.Fb(xg,yg) > dp2) & (hh_m > mx | isnan(hh_m));
+                    limidx = (tmpz < dp1 & tmpz > dp2) & ...
+                             (hh_m > mx | isnan(hh_m));
                     
                     hh_m( limidx ) = mx;
                 end
@@ -858,8 +885,7 @@ function plot(obj,type)
                     dp1  = param(2);
                     dp2  = param(3);
                     
-                    limidx = (feat.Fb(xg,yg) < dp1 & ...
-                        feat.Fb(xg,yg) > dp2) ;
+                    limidx = (tmpz < dp1 & tmpz > dp2) ;
                     
                     dmy( limidx ) = lim;
                 end
@@ -895,7 +921,6 @@ function plot(obj,type)
                 end
                 % Estimate wavespeed in ocean (second term represents orbital
                 % velocity at 0 degree phase for 1-m amp. wave).
-                tmpz    = feat.Fb(xg,yg);
                 grav = 9.807;
                 % Limit the minimum depth to 1 m (ignore overland)
                 tmpz(tmpz > - 1) = -1;
diff --git a/@geodata/geodata.m b/@geodata/geodata.m
index 9de8994d..a4ca80be 100644
--- a/@geodata/geodata.m
+++ b/@geodata/geodata.m
@@ -306,6 +306,8 @@
             
             obj = ParseDEM(obj) ;
             
+            
+            
         end
         
         function obj = ParseShoreline(obj)
@@ -556,9 +558,14 @@
                     clear x y demz
                 else
                     % main interpolant
-                    obj.Fb   = griddedInterpolant({x,y},demz,...
-                        'linear','nearest');
-                    obj.x0y0 = [x(1),y(1)];
+                    if obj.BACKUPdemfile~=0
+                        obj.Fb   = griddedInterpolant({x,y},demz,...
+                            'linear','none'); % no extrapolation (so use Fb2)
+                    else
+                        obj.Fb   = griddedInterpolant({x,y},demz,...
+                            'linear','nearest');
+                        obj.x0y0 = [x(1),y(1)];
+                    end
                     % clear data from memory
                     clear x y demz
                 end
diff --git a/@msh/msh.m b/@msh/msh.m
index 356eeb1e..1fbb9b71 100644
--- a/@msh/msh.m
+++ b/@msh/msh.m
@@ -871,7 +871,8 @@ function write(obj,fname,type)
             %                         relevant for CA interpolation method).
             %                         default value N=1.
             %
-            %        nan - 'fill' to fill in any NaNs appearing in bathy
+            %        nan - 'fill' to fill in any NaNs everywhere
+            %            -'fillinside' to fill NaNs only in DEM extents
             %
             %   mindepth - ensure the minimum depth is bounded in the
             %                         interpolated region
diff --git a/@msh/private/GridData.m b/@msh/private/GridData.m
index 0a673f05..21a58b74 100644
--- a/@msh/private/GridData.m
+++ b/@msh/private/GridData.m
@@ -30,7 +30,8 @@
 %                         relevant for CA interpolation method). 
 %                         default value N=1. 
 %
-%        nan (optional) - 'fill' to fill in any NaNs appearing in bathy
+%        nan (optional) - 'fill' to fill in any NaNs everywhere
+%                        -'fillinside' to fill NaNs only in DEM extents
 %
 %   mindepth (optional) - ensure the minimum depth is bounded in the 
 %                         interpolated region 
diff --git a/utilities/Calc_f13.m b/utilities/Calc_f13.m
index 5bc01a54..36d658ec 100644
--- a/utilities/Calc_f13.m
+++ b/utilities/Calc_f13.m
@@ -11,6 +11,7 @@
 %              'Mn' ('mannings_n_at_sea_floor')
 %              'Ss' ('surface_submergence_state')
 %              'Re' ('initial_river_elevation')
+%              'Ad' ('advection_state')
 %
 %            3) then either: 
 %              'inpoly' followed by...
@@ -47,6 +48,9 @@
 elseif strcmpi(attribute,'Re')
     attrname = 'initial_river_elevation';
     default_val = 0;
+elseif strcmpi(attribute,'Ad')
+    attrname = 'advection_state';
+    default_val = -999;
 else
     error(['Attribute ' attribute ' not currently supported'])
 end