diff --git a/@dash/dash.m b/@dash/dash.m
index 93f32c52..d0e85f8b 100644
--- a/@dash/dash.m
+++ b/@dash/dash.m
@@ -4,6 +4,9 @@
     
     methods (Static)
         
+        % Analysis
+        dist = haversine(latlon1, latlon2);
+        
         % Misc
         [names, lon, lat, coord, lev, time, run, var]  = dimensionNames;
         convertToV7_3(filename);
diff --git a/@dash/haversine.m b/@dash/haversine.m
new file mode 100644
index 00000000..0253ad54
--- /dev/null
+++ b/@dash/haversine.m
@@ -0,0 +1,53 @@
+function[dist] = haversine( latlon1, latlon2 )
+%% Computes the distance between lat-lon coordinates
+%
+% dist = haversine( latlon )
+% Computes the distance between each set of latitude-longitude coordinates
+%
+% dist = haversine( latlon1, latlon2 )
+% Computes the distance between each coordinate listed in latlon1 and each
+% coordinate listed in latlon2.
+%
+% ----- Inputs -----
+%
+% latlon: A set a latitude-longitude coordinates. A numeric matrix with two
+%    columns. First column is the latitude coordinates. Second column is
+%    the longitude coordinate.
+%
+% ----- Outputs -----
+%
+% dist: The distance between the coordinates. 
+%
+%    If you provided a single input, dist is a symmetric matrix with each
+%    rows and columns corresponding to the coordinates listed in latlon.
+%
+%    If you provided two inputs, the rows of dist correspond to latlon1 and
+%    the columns correspond to latlon2.
+
+% Default for the second set of coordinates
+if ~exist('latlon2','var') || isempty(latlon2)
+    latlon2 = latlon1;
+end
+
+% Radius of the Earth
+R = 6378.137;
+
+% Convert to radians
+latlon1 = latlon1 * pi/180;
+latlon2 = latlon2 * pi/180;
+
+% Transpose latlon2 for singleton expansion
+latlon2 = latlon2';
+
+% Get the change in lat and lon
+dLat = latlon1(:,1) - latlon2(1,:);
+dLon = latlon1(:,2) - latlon2(2,:);
+
+% Get haversine function of the central angle
+a = sin(dLat/2).^2 + ( cos(latlon1(:,1)) .* cos(latlon2(1,:)) .* sin(dLon/2).^2);
+c = 2 * atan2( sqrt(a), sqrt(1-a) );
+
+% Get the distance
+dist = R * c;
+
+end
\ No newline at end of file
diff --git a/@ensembleMetadata/closestLatLon.m b/@ensembleMetadata/closestLatLon.m
new file mode 100644
index 00000000..0fcf3177
--- /dev/null
+++ b/@ensembleMetadata/closestLatLon.m
@@ -0,0 +1,54 @@
+function[rows] = closestLatLon(obj, latlon, varName, verbose)
+%% Returns the state vector rows closest to a set of lat-lon coordinates
+%
+% rows = obj.closestLatLon( coordinate )
+% Returns the state vector rows closest to the given coordinate.
+%
+% rows = obj.closestLatLon( coordinate, varName )
+% Returns the state vector rows that are closest to the given coordinate
+% for a particular variable.
+%
+% rows = obj.closestLatLon( coordinate, varName, verbose )
+% Optionally disable console messages.
+%
+% ----- Inputs -----
+%
+% coordinate: A latitude-longitude coordinate. A vector with two elements.
+%    The first element is the latitude coordinate and the second element is
+%    the longitude coordinate.
+%
+% varName: The name of a variable in the state vector ensemble. A string.
+%
+% verbose: A scalar logical that indicates whether to return console
+%    messages when determining coordinates (true -- default) or not (false)
+%
+% ----- Outputs -----
+%
+% rows: The rows of the state vector that are closest to the selected
+%    lat-lon coordinate.
+
+% Defaults
+if ~exist('varName','var') || isempty(varName)
+    varName = [];
+end
+if ~exist('verbose','var') || isempty(verbose)
+    verbose = [];
+end
+
+% Error check the user-specified coordinate
+assert( isnumeric(latlon) & isvector(latlon), 'latlon must be a numeric vector.');
+assert(numel(latlon)==2, 'latlon must have two elements.');
+dash.assertRealDefined(latlon, 'latlon');
+assert(abs(latlon(1))<=90, 'latitude coordinates must be between -90 and 90.');
+
+% Get the state vector lat-lon coordinates and get the distance to each
+svLatLon = obj.latlon(varName, verbose);
+dist = dash.haversine(svLatLon, latlon);
+
+% Find the state vector rows associated with the minimum distances
+rows = find(dist==min(dist));
+if ~isempty(varName)
+    rows = obj.findRows(varName, rows);
+end
+
+end
\ No newline at end of file
diff --git a/@ensembleMetadata/columns.m b/@ensembleMetadata/columns.m
new file mode 100644
index 00000000..20248ed8
--- /dev/null
+++ b/@ensembleMetadata/columns.m
@@ -0,0 +1,96 @@
+function[meta] = columns(obj, cols, varNames, dims, alwaysStruct)
+%% Returns the metadata at a particular column across the ensemble.
+%
+% metaStruct = obj.columns( cols )
+% Returns a structure with the metadata at the specified ensemble members
+% for each variable and dimension in the state vector ensemble.
+%
+% metaStruct = obj.columns( cols, varNames )
+% Returns metadata for the specified variables for all dimensions.
+%
+% metaStruct = obj.columns( cols, varNames, dimNames )
+% Returns metadata for the specified variables and dimensions.
+%
+% metaStruct = obj.columns( cols, varNames, dimCell )
+% Specify different dimensions for different variables.
+%
+% meta = obj.columns( cols, varName, dimName )
+% Return metadata for a single variable and single dimension directly as an
+% array.
+%
+% metaStruct = obj.columns( cols, varName, dimName, alwaysStruct )
+% Specify to always return output as a structure.
+%
+% ----- Inputs -----
+%
+% cols: The indices of the columns (ensemble members) for which to return
+%    metadata. Either a vector of linear indices, or a logical vector. If a
+%    logical vector, must have one element per ensemble member.
+%
+% varNames: A list of variables in the state vector ensemble. A string
+%    vector or cellstring vector.
+%
+% dimName: A list of dimension names. A string vector or cellstring vector.
+%
+% dimCell: A cell vector with one element per listed dimension name. Each
+%    element contains a list of dimension names for the corresponding
+%    variable.
+%
+% alwaysStruct: A scalar logical that indicates if metadata should always
+%    be returned as a structure (true) or not (false -- default)
+%
+% ----- Outputs -----
+%
+% metaStruct: A structure holding metadata at the specified column for the
+%    variables and dimensions.
+%
+% meta: An array with metadata for the specified variable and dimension.
+%    Has one row per specified ensemble member.
+
+% Error check the columns
+dash.assertVectorTypeN(cols, 'numeric', [], 'cols');
+cols = dash.checkIndices(cols, 'cols', obj.nEns, 'the number of ensemble members');
+
+% Default variable names. Error check. Indices
+if ~exist('varNames','var') || isempty(varNames)
+    varNames = obj.variableNames;
+end
+varNames = dash.assertStrList(varNames, 'varNames');
+v = dash.checkStrsInList(varNames, obj.variableNames, 'varNames', 'variables in the state vector');
+nVars = numel(v);
+
+% Default dimension names. Error check. Propagate strings into cells
+if ~exist('dims','var') || isempty(dims)
+    dims = obj.dims(v);
+elseif ~iscell(dims)
+    dims = dash.assertStrList(dims, 'dims');
+    dims = repmat({dims}, [nVars,1]);
+else
+    for d = 1:numel(dims)
+        dims{d} = dash.assertStrList(dims{d}, sprintf('Element %.f of dimCell', d));
+    end
+end
+dash.assertVectorTypeN(dims, 'cell', nVars, 'dims');
+
+% Default and error check alwaysStruct
+if ~exist('alwaysStruct','var') || isempty(alwaysStruct)
+    alwaysStruct = false;
+end
+dash.assertScalarType(alwaysStruct, 'alwaysStruct', 'logical', 'logical');
+
+% Initialize the output structure
+meta = struct();
+
+% Get the metadata for each variable
+for k = 1:nVars
+    var = varNames(k);
+    indices = repmat({cols}, [1, numel(dims{k})]);
+    meta.(var) = obj.variable(var, dims{k}, 'ensemble', indices, true);
+end
+
+% Optionally return array
+if ~alwaysStruct && nVars==1 && numel(dims{1})==1
+    meta = meta.(var).(dims{1});
+end
+
+end
\ No newline at end of file
diff --git a/@ensembleMetadata/ensembleMetadata.m b/@ensembleMetadata/ensembleMetadata.m
index e4a58166..b33d1d7b 100644
--- a/@ensembleMetadata/ensembleMetadata.m
+++ b/@ensembleMetadata/ensembleMetadata.m
@@ -163,9 +163,14 @@
         obj = rename(obj, newName);
 
         [V, meta] = regrid(obj, X, varName, dimOrder, d, keepSingletons);
-        meta = variable(obj, varName, dims, type, indices);
+        
+        meta = variable(obj, varName, dims, type, indices, alwaysStruct);
         meta = dimension(obj, dim, alwaysStruct);
-        [latlon] = coordinates(obj, dim, verbose);
+        meta = rows(obj, rows, dims, fullStruct);
+        meta = columns(obj, cols, varNames, dims, alwaysStruct);
+        
+        [latlon] = latlon(obj, varName, verbose);
+        rows = closestLatLon(obj, latlon, varName, verbose);
         rows = findRows(obj, varName, varRows);
         [nState, nEns] = sizes(obj, vars);
         
diff --git a/@ensembleMetadata/latlon.m b/@ensembleMetadata/latlon.m
new file mode 100644
index 00000000..1fa91115
--- /dev/null
+++ b/@ensembleMetadata/latlon.m
@@ -0,0 +1,205 @@
+function[latlon] = coordinates(obj, varName, verbose)
+%% Returns lat-lon coordinate metadata down the state vector. Reads metadata
+% from the "lat", "lon", and "coord" dimensions as appropriate. Returns NaN
+% coordinates for any state vector elements that use a spatial mean or
+% that lack lat-lon coordinate metadata.
+%
+% latlon = obj.coordinates
+% Returns a lat-lon coordinate for each element of a state vector.
+%
+% latlon = obj.coordinates(varName)
+% Returns the lat-lon coordinates for a particular variable.
+%
+% latlon = obj.coordinates(varName, verbose)
+% Specify whether or not to print notifications to the console. Default is
+% to display notifications.
+%
+% *** Requirements *** 
+% 1. Metadata extracted from the "lat" and "lon" dimensions must be numeric
+% and have a single column -- otherwise, NaN coordinates are returned.
+%
+% 2. Metadata extracted from the "coord" dimension must be numeric and have
+% two columns (one column for lat, and one column for lon). The method will
+% attempt to determine which column is lat and which column is lon
+% automatically. If it cannot determine which column corresponds to which
+% dimension, prints a notification to the console and uses the first column
+% as latitude coordinates.
+%
+% ----- Inputs -----
+%
+% varName: The name of a variable in a state vector ensemble. A string. Use
+%    an empty array to return coordinates for all variables.
+%
+% verbose: A scalar logical indicating whether to print notification
+%    messages to the console (true -- default) or not (false).
+%
+% ----- Outputs -----
+%
+% latlon: A numeric matrix of lat-lon coordinates. Has one row per state
+%    vector element and two columns. The first column contains latitude
+%    coordinates and the second contains longitude coordinates. The
+%    coordinates for any state vector element with a spatial mean, missing
+%    coordinate metadata, or unrecognized coordinate format will be
+%    returned as NaN.
+
+% Default and error check for notifications
+if ~exist('verbose','var') || isempty(verbose)
+    verbose = true;
+end
+dash.assertScalarType(verbose, 'verbose', 'logical', 'logical');
+notified = false;
+
+% Get the dimension names and variable(s). Preallocate the coordinates
+[~, lonName, latName, coordName] = dash.dimensionNames;
+if ~exist('varName','var') || isempty(varName)
+    v = 1:numel(obj.nEls);
+    latlon = NaN(obj.varLimit(end), 2);
+else
+    dash.assertStrFlag(varName, 'varName');
+    v = dash.checkStrsInList(varName, obj.variableNames, 'varName', 'variable in the state vector');
+    latlon = NaN(obj.nEls(v), 2);
+end
+
+% Get state metadata for the variable
+for k = 1:numel(v)
+    var = obj.variableNames(v(k));
+    meta = obj.metadata.(var).state;
+    
+    % Check whether the metadata is stored as "lat" and "lon", or "coord"
+    hasdim = false(1, 3); % lat, lon, coord
+    dims = [lonName, latName, coordName];
+    for d = 1:3
+        if ismember(dims(d), obj.dims{v(k)}) && (~isscalar(meta.(dims(d))) || ~isnan(meta.(dims(d))))
+            hasdim(d) = true;
+        end
+    end
+    
+    % Notify user if there is both lat/lon and coord metadata
+    hasdata = true;
+    if (hasdim(1)||hasdim(2)) && hasdim(3)
+        hasdata = false;
+        if verbose
+            bothMetadataTypesWarning(var, dims(find(hasdim,1)), coordName);
+            notified = true;
+        end
+    
+    % Notify if there is lat, but not lon, or vice versa
+    elseif any(hasdim(1:2)) && ~all(hasdim(1:2))
+        hasdata = false;
+        if verbose
+            missingHalfWarning(var, dims(find(hasdim,1)), dims(find(~hasdim,1)));
+            notified = true;
+        end
+        
+    % Notify if there is no metadata
+    elseif all(~hasdim)
+        hasdata = false;
+        if verbose
+            missingAllWarning(var, dims);
+            notified = true;
+        end
+    end
+    
+    % Get lat-lon metadata. Check formatting and spatial means.
+    if hasdata && ~hasdim(3)
+        latMeta = obj.variable(var, dims(2), true);
+        [hasdata, notified] = checkMetadata(latMeta, 1, var, dims(2), verbose, notified);
+        if hasdata
+            lonMeta = obj.variable(var, dims(1), true);
+            [hasdata, notified] = checkMetadata(lonMeta, 1, var, dims(1), verbose, notified);
+            varCoords = [latMeta, lonMeta];
+        end
+        
+    % Otherwise, get coordinate metadata. Check format and spatial mean
+    elseif hasdata
+        coordMeta = obj.variable(var, dims(3), true);
+        [hasdata, notified] = checkMetadata(coordMeta, 2, var, dims(3), verbose, notified);
+        
+        % Attempt to determine which column is lat and which is lon
+        if hasdata
+            over90 = any(abs(coordMeta)>90,1);
+            
+            % Notify user if using default
+            if all(over90) || all(~over90)
+                if verbose
+                    notifyDefaultColumn(var, dims(3));
+                    notified = true;
+                end
+                
+            % Otherwise, set columns and notify
+            else
+                latcol = find(~over90);
+                loncol = find(over90);
+                varCoords = coordMeta(:, [latcol, loncol]);
+                if verbose
+                    notifySelectedColumn(var, dims(3), latcol, loncol);
+                    notified = true;
+                end
+            end
+        end
+    end
+    
+    % If everything was successful, add the coordinates to the output
+    if hasdata && numel(v)==1
+        latlon = varCoords;
+    elseif hasdata
+        rows = obj.varLimit(v(k),1):obj.varLimit(v(k),2);
+        latlon(rows,:) = varCoords;
+    end
+end
+
+% Format the console after messages
+if notified
+    fprintf('\n');
+end
+
+end
+
+% Helper function to support DRY code
+function[hasdata, notified] = checkMetadata(meta, nCols, var, dim, verbose, notified)
+           
+% Check numeric with correct columns
+hasdata = true;
+if ~isnumeric(meta) || size(meta,2)~=nCols
+    hasdata = false;
+    if verbose
+        fprintf(['\nUsing NaN coordinates for "%s" because the metadata for ',...
+            'the "%s" dimension is not a numeric marix with %.f column(s).'],...
+            var, dim, nCols);
+        notified = true;
+    end
+    
+% Check for spatial mean
+elseif size(meta,3)~=1
+    hasdata = false;
+    if verbose
+        fprintf(['\nUsing NaN coordinates for "%s" because it takes a spatial ',...
+            'mean over the "%s" dimension.'], var, dim);
+        notified = true;
+    end
+end
+
+end
+
+% Long warnings
+function[] = bothMetadataTypesWarning(var, dim1, dim2)
+fprintf(['\nUsing NaN coordinates for "%s" because it has metadata for both ',...
+    'the "%s" and "%s" dimensions.'], var, dim1, dim2);
+end
+function[] = missingHalfWarning(var, dim1, dim2)
+fprintf(['\nUsing NaN coordinates for "%s" because it has metadata for the ',...
+    '"%s" dimension, but not the "%s" dimension.'], var, dim1, dim2);
+end
+function[] = missingAllWarning(var, dims)
+fprintf(['\nUsing NaN coordinates for "%s" because it does not have metadata ',...
+    'for any of the "%s", "%s", or "%s" dimensions.'], var, dims(1), dims(2), dims(3));
+end
+function[] = notifyDefaultColumn(var, dim)
+fprintf('\nCould not determine which column of the "%s" metadata for ',...
+    'variable "%s" is for latitude, and which column is for longitude.',...
+    'Using the first column as latitude coordinates.', dim, var);
+end
+function[] = notifySelectedColumn(var, dim, latcol, loncol)
+fprintf(['\nFor the "%s" metadata of variable "%s": Using column %.f for ',...
+    'latitude, and column %.f for longitude.\n'], dim, var, latcol, loncol);
+end     
\ No newline at end of file
diff --git a/@ensembleMetadata/rows.m b/@ensembleMetadata/rows.m
new file mode 100644
index 00000000..cea7b76e
--- /dev/null
+++ b/@ensembleMetadata/rows.m
@@ -0,0 +1,104 @@
+function[meta] = rows(obj, rows, dims, fullStruct)
+%% Returns the metadata at a particular row down the state vector.
+%
+% varStruct = obj.rows( rows )
+% metaStruct = obj.rows( rows )
+% Returns a structure with the metadata at the specified row for each
+% dimension. If the rows are for a single state vector variable, the fields
+% of the structure are the different dimensions. If the rows cover multiple
+% state vector variables, then the fields of the structure are the
+% variables covered by the rows. Each field contains a sub-structure with
+% the metadata for the dimensions of the variable. Each sub-structure also
+% indicates the specified rows that are in the variable.
+%
+% varStruct = obj.rows(rows, dims)
+% metaStruct = obj.rows(rows, dims)
+% meta = obj.rows(rows, dim)
+% Only returns metadata for specified dimensions. If the rows are for a
+% single state vector variable, and only a single dimension is specified,
+% then returns the metadata directly as an array.
+%
+% varStruct = obj.rows(rows, dim, fullStruct)
+% Specify to always return the full structure (in which each field
+% corresponds to a particular variable).
+%
+% varNames = obj.rows(rows, 0)
+% Return the name of the state vector variable associated with each row.
+
+% Default and error check for fullStruct
+if ~exist('fullStruct','var') || isempty(fullStruct)
+    fullStruct = false;
+end
+dash.assertScalarType(fullStruct, 'fullStruct', 'logical', 'logical');
+
+% Error check the rows
+rows = dash.checkIndices(rows, 'rows', obj.varLimit(end), 'the number of elements in the state vector');
+if ~isrow(rows)
+    rows = rows';
+end
+
+% Get the variables associated with the rows
+isVar = obj.varLimit(:,1)<=rows & obj.varLimit(:,2)>=rows;
+[whichVar, ~] = find(isVar);
+v = unique(whichVar);
+nVars = numel(v);
+
+% Default dimensions
+if ~exist('dims','var') || isempty(dims)
+    dims = unique(cat(2, obj.dims{v}));
+    
+% Only return variable names if using the 0 dimension
+elseif isequal(dims, 0)
+    meta = obj.variableNames(whichVar);
+    return;
+    
+% Error check user-specified dimensions
+else
+    allDims = unique(cat(2, obj.dims{:}));
+    dash.checkStrsInList(dims, allDims, 'dims', 'a dimension of any variable in the state vector ensemble');
+    dims = unique(string(dims));
+end
+nUserDims = numel(dims);
+
+% Initialize the output structure
+if fullStruct || nVars>1
+    meta = struct();
+end
+
+% For each variable, get the name, variable rows, and dimensions
+for k = 1:nVars
+    varName = obj.variableNames(v(k));
+    currentRows = rows(whichVar==v(k));
+    varRows = currentRows - obj.varLimit(v(k),1) + 1;
+    varDims = dims(ismember(dims, obj.dims{v(k)}));
+    nDims = numel(varDims);
+    
+    % Get the variable metadata
+    if nDims==0
+        varMeta = [];
+    else
+        indices = repmat({varRows}, [1 nDims]);
+        varMeta = obj.variable(varName, varDims, 'state', indices, true);
+    end 
+    
+    % Return the full output structure
+    if fullStruct || nVars>1
+        meta.(varName) = varMeta;
+        meta.(varName).rows = currentRows;
+        meta.(varName) = orderfields(meta.(varName), ["rows"; varDims(:)]);
+        
+    % Return a metadata structure for a single variable
+    elseif nVars==1 && nUserDims>1
+        meta = varMeta;
+        meta.stateVectorVariable = varName;
+        meta = orderfields(meta, ["stateVectorVariable"; varDims(:)]);
+        
+    % Return an array directly
+    elseif nDims==0
+        meta = varMeta;
+    else
+        meta = varMeta.(varDims(1));
+    end
+end
+
+end
\ No newline at end of file
diff --git a/@ensembleMetadata/variable.m b/@ensembleMetadata/variable.m
index efaf7441..c009f275 100644
--- a/@ensembleMetadata/variable.m
+++ b/@ensembleMetadata/variable.m
@@ -1,4 +1,4 @@
-function[meta] = variable(obj, varName, dims, type, indices)
+function[meta] = variable(obj, varName, dims, type, indices, alwaysStruct)
 %% Returns metadata down the state vector or across the ensemble for a
 % variable in a state vector ensemble.
 %
@@ -23,6 +23,9 @@
 % [...] = obj.lookup(varName, dims, type/returnState, indexCell)
 % Return the metadata for specific elements of a variable.
 %
+% [...] = obj.lookup(varName, dim, type/returnState, indices, alwaysStruct)
+% Specify if metadata should always be returned as a structure.
+%
 % ----- Inputs -----
 %
 % varName: The name of a variable in the state vector. A string.
@@ -55,6 +58,9 @@
 %    array, returns metadata at all elements down the state vector or
 %    across the ensemble, as appropriate.
 %
+% alwaysStruct: A scalar logical that indicates if metadata should always
+%    be returned as a structure (true) or not (false -- default)
+%
 % ----- Outputs -----
 %
 % meta: The metadata for a single dimension. A matrix. If the metadata is
@@ -91,8 +97,14 @@
     name = 'indexCell';
 end
 
+% Default and error check for alwaysStruct
+if ~exist('alwaysStruct','var') || isempty(alwaysStruct)
+    alwaysStruct = false;
+end
+dash.assertScalarType(alwaysStruct, 'alwaysStruct', 'logical', 'logical');
+
 % Initialize output structure
-if nDims > 1
+if nDims>1 || alwaysStruct
     meta = struct();
 end
 
@@ -133,7 +145,7 @@
     end
     
     % Return as array or structure
-    if nDims==1
+    if nDims==1 && ~alwaysStruct
         meta = dimMeta;
     else
         meta.(dims(k)) = dimMeta;
diff --git a/@stateVector/add.m b/@stateVector/add.m
index a112e029..82b8de7f 100644
--- a/@stateVector/add.m
+++ b/@stateVector/add.m
@@ -1,8 +1,17 @@
-function[obj] = add(obj, varName, file, autoCouple, overlap)
+function[obj] = add(obj, varNames, files, autoCouple, overlap)
 %% Adds a variable to a stateVector.
 %
 % obj = obj.add(varName, file)
-% Adds a variable to the state vector from a .grid file.
+% Adds a variable to the state vector and specifies the .grid file that
+% holds the data for the variable.
+%
+% obj = obj.add(varNames, file)
+% Adds multiple variables to the state vector whose data are all stored in
+% the same .grid file.
+%
+% obj = obj.add(varNames, files)
+% Adds multiple variables to the state vector and specifies the associated
+% .grid files.
 %
 % obj = obj.add(varName, file, autoCouple)
 % Specify whether the variable should be automatically coupled to other
@@ -41,28 +50,43 @@
     overlap = false;
 end
 
-% Error check, use string internally, editable
+% Ensure the vector is still editable
 obj.assertEditable;
-dash.assertScalarType(overlap, 'overlap', 'logical', 'logical');
-dash.assertScalarType(autoCouple, 'autoCouple', 'logical', 'logical');
-dash.assertStrFlag(varName, 'varName');
-varName = string(varName);
 
-% Check the name is a valid variable name and not a duplicate
-obj.checkVariableNames(varName, [], 'varName', 'add a new variable to');
-vars = obj.variableNames;
-vars(end+1) = varName;
+% Error check the variable names. Get the number of new variables
+varNames = dash.assertStrList(varNames, 'varNames');
+obj.checkVariableNames(varNames, [], 'varNames', 'add new variables to');
+nNew = numel(varNames);
+
+% Error check the other inputs.
+files = dash.assertStrList(files, 'files');
+dash.assertVectorTypeN(autoCouple, 'logical', [], 'autoCouple');
+dash.assertVectorTypeN(overlap, 'logical', [], 'overlap');
 
-% Create the new variable (error checks file).
-newVar = stateVectorVariable(varName, file);
-obj.variables = [obj.variables; newVar];
+% Check sizes and propagate scalar inputs
+fields = {files, autoCouple, overlap};
+fieldNames = {'files','autoCouple','overlap'};
+for f = 1:numel(fields)
+    nEls = numel(fields{f});
+    if ~isscalar(fields{f}) && nEls~=nNew
+        error('%s may either have 1 element or %.f elements (1 per new variable).', fieldNames{f}, nNew);
+    end
+    fields{f} = repmat(fields{f}(:), [nNew-nEls+1, 1]);
+end
+[files, autoCouple, overlap] = fields{:};
+
+% Create each new variable (also error checks the .grid files)
+for v = 1:nNew
+    newVar = stateVectorVariable(varNames(v), files(v));
+    obj.variables = [obj.variables; newVar];
+    obj.coupled(end+1, end+1) = true;
+end
 
-% Update variable coupling and overlap
-obj.overlap(end+1, 1) = overlap;
-obj.auto_Couple(end+1, 1) = autoCouple;
-obj.coupled(end+1, end+1) = true;
-if autoCouple
-    obj = obj.couple( vars(obj.auto_Couple) );
+% Update auto coupling and overlap
+obj.overlap(end+(1:nNew), 1) = overlap;
+obj.auto_Couple(end+(1:nNew), 1) = autoCouple;
+if any(autoCouple)
+    obj = obj.couple( obj.variableNames(obj.auto_Couple) );
 end
 
 end
\ No newline at end of file
diff --git a/@stateVector/checkVariableNames.m b/@stateVector/checkVariableNames.m
index 88744ff0..43bc358f 100644
--- a/@stateVector/checkVariableNames.m
+++ b/@stateVector/checkVariableNames.m
@@ -15,19 +15,20 @@
 % methodName: The name of the action being attempted.
 
 % Default index for the new names
-if isempty(newNames)
+if isempty(v)
     v = numel(obj.variables)+(1:numel(newNames));
 end
 
 % Check that the new names are valid MATLAB variable names
-if any(~isvarname(newNames))
-    bad = find(~isvarname(newNames),1);
-    str = sprintf('The value of %s (%s)', inputName, newNames);
-    if numel(newNames)>1
-        str = sprintf('Element %.f of %s (%s)', bad, inputName, newNames(bad));
+for n = 1:numel(newNames)
+    if ~isvarname(newNames(n))
+        str = sprintf('The value of %s (%s)', inputName, newNames(n));
+        if numel(newNames)>1
+            str = sprintf('Element %.f of %s (%s)', n, inputName, newNames(n));
+        end          
+        error(['%s is not a valid MATLAB variable name. Valid names must start ',...
+            'with a letter and may only include letters, numbers, and underscores.'], str);
     end
-    error(['%s is not a valid MATLAB variable name. Valid names must start ',...
-        'with a letter and may only include letters, numbers, and underscores.'], str);
 end
 
 % Combine new names with old
@@ -37,9 +38,9 @@
 % Check for duplicates
 [uniqNames, loc] = unique(names);
 nNames = numel(names);
-if nNames < numel(uniqNames)
+if nNames > numel(uniqNames)
     bad = find(~ismember(1:nNames, loc), 1);
-    error(['Cannot %s %s because there would be multiple variables named "%s".', ...
+    error(['Cannot %s %s because there would be multiple variables named "%s". ', ...
         'If you want to change existing variable names, see "stateVector.renameVariable".'],...
         methodName, obj.errorTitle, names(bad));
 end
diff --git a/@stateVectorVariable/weightedMean.m b/@stateVectorVariable/weightedMean.m
index b2fb96d2..f603d8ea 100644
--- a/@stateVectorVariable/weightedMean.m
+++ b/@stateVectorVariable/weightedMean.m
@@ -75,7 +75,7 @@
 
     % If weights is empty, this is an unweighted mean
     for k = 1:nDims
-        if isempty(obj.weightCell{d(k)})
+        if isempty(weights{k})
             obj.hasWeights(d(k)) = false;
             
         % Otherwise, error check weights and update
diff --git a/dataSource.m b/dataSource.m
index 0957dd82..c7c8a9b4 100644
--- a/dataSource.m
+++ b/dataSource.m
@@ -4,7 +4,7 @@
     % implement functionality for different types of data sourcess. (For
     % example, netCDF and .mat files and opendap files).
     
-    properties
+    properties (SetAccess = protected)
         source; % The data source. A filename or opendap url
         var; % (For hdf data sources) The name of the variable.
         dataType;  % The type of data in the file.
diff --git a/matSource.m b/matSource.m
index af829ee4..9f23db20 100644
--- a/matSource.m
+++ b/matSource.m
@@ -1,7 +1,7 @@
 classdef matSource < dataSource
     %% Reads data from a .mat file data source
     
-    properties
+    properties (SetAccess = private)
         m; % A matfile object
     end
     properties (Constant, Hidden)
diff --git a/ncSource.m b/ncSource.m
index 1885aebd..e987694f 100644
--- a/ncSource.m
+++ b/ncSource.m
@@ -2,7 +2,7 @@
     %% Used to read data from source based on a NetCDF format. Includes 
     % local NetCDF files and OPeNDAP requests.
     
-    properties
+    properties (SetAccess = protected)
         nDims; % The number of dimensions recorded in the NetCDF
     end
     
diff --git a/opendapSource.m b/opendapSource.m
index e3848830..1426910b 100644
--- a/opendapSource.m
+++ b/opendapSource.m
@@ -2,7 +2,7 @@
     %% Reads data accessed via an OPeNDAP url. Attempts to load and save 
     % the entire variable when using repeated loads to increase speed.
     
-    properties
+    properties (SetAccess = private)
         X; % The loaded and saved dataset
         attemptFullLoad; % Whether to attempt to load the entire dataset
         saved; % Whether the dataset is currently saved