Index: /issm/trunk/src/m/kml/kml_mesh_elem.m
===================================================================
--- /issm/trunk/src/m/kml/kml_mesh_elem.m	(revision 6417)
+++ /issm/trunk/src/m/kml/kml_mesh_elem.m	(revision 6417)
@@ -0,0 +1,159 @@
+%%
+%  create kml polygons for the element mesh.
+%
+%  []=kml_mesh_elem(fid,md,params)
+%
+%  where the required input is:
+%    fid           (numeric, file ID of .kml file)
+%    md            (model, model class object)
+%
+%  the optional input is:
+%    params        (string/numeric, parameter names and values)
+%
+%  and the optional input is:
+%    data          (numeric, element or nodal results data)
+%    alt           (numeric, altitude for polygons, default 10000)
+%    cmin          (numeric, minimum of color map)
+%    cmax          (numeric, maximum of color map)
+%    cmap          (char or numeric, colormap definition)
+%
+function []=kml_mesh_elem(varargin)
+
+if ~nargin
+    help kml_mesh_elem
+    return
+end
+
+%%  process input data
+
+iarg=1;
+if (nargin >= 1)
+    fid=varargin{1};
+end
+if ~exist('fid','var') || isempty(fid) || (fid < 0)
+    error('File ID ''%d'' must be open.',fid);
+end
+
+iarg=iarg+1;
+if (nargin >= 2)
+    md=varargin{2};
+end
+if ~exist('md','var') || isempty(md) || ~isa(md,'model')
+    error(['Model ''' inputname(iarg) ''' is unrecognized class ''' class(md) '''.']);
+end
+
+%  parameters
+
+iarg=iarg+1;
+while (iarg <= nargin-1)
+    if ischar(varargin{iarg})
+        eval([varargin{iarg} '=varargin{iarg+1};']);
+        if (numel(varargin{iarg+1}) <= 20)
+            disp([varargin{iarg} '=' any2str(varargin{iarg+1}) ';']);
+        else
+            disp([varargin{iarg} '=' string_size(varargin{iarg+1}) ' ' class(varargin{iarg+1}) ';']);
+        end
+        if strcmpi(varargin{iarg},'data')
+            cdata=inputname(iarg+1);
+        end
+    else
+        error(['''' any2str(varargin{iarg}) ''' is not a parameter name.']);
+    end
+    iarg=iarg+2;
+end
+
+if exist('data','var') && ~isempty(data)
+    if     (numel(data)==md.numberofelements)
+        edata=data;
+    elseif (numel(data)==md.numberofgrids)
+        ndata=data;
+        display('Averaging nodal data to element data.');
+        edata=zeros(1,md.numberofelements);
+        for i=1:size(md.elements,1)
+            for j=1:size(md.elements,2)
+                edata(i)=edata(i)+ndata(md.elements(i,j));
+            end
+            edata(i)=edata(i)/size(md.elements,2);
+        end
+    else
+        error(['Data has incorrect number of ' num2str(numel(data)) ' values.']);
+    end
+end
+
+%  colormap command operates on a figure, so create an invisible one
+%  (could also directly call colormaps, e.g. jet(64), but risky)
+
+hfig=figure('Visible','off');
+if exist('cmap','var')
+    colormap(cmap)
+end
+cmap=colormap;
+close(hfig)
+    
+if exist('edata','var')
+    if ~exist('cmin','var')
+        cmin=min(min(edata));
+    end
+    if ~exist('cmax','var')
+        cmax=max(max(edata));
+    end
+end
+
+if ~exist('alt','var')
+    alt=10000;
+end
+
+%%  write folder for mesh
+
+fprintf(fid,'    <Folder>\n');
+if exist('cdata','var') && ~isempty(cdata)
+    fprintf(fid,'      <name>Data: %s</name>\n',cdata);
+else
+    fprintf(fid,'      <name>Mesh</name>\n');
+end
+fprintf(fid,'      <visibility>1</visibility>\n');
+fprintf(fid,'      <description>Elements=%d, Grids=%d</description>\n',...
+    md.numberofelements,md.numberofgrids);
+
+%  write each element as a polygon
+
+disp(['Writing ' num2str(size(md.elements,1)) ' tria elements as KML polygons.']);
+for i=1:size(md.elements,1)
+    fprintf(fid,'      <Placemark>\n');
+    fprintf(fid,'        <name>Element %d</name>\n',i);
+    fprintf(fid,'        <visibility>1</visibility>\n');
+    if exist('edata','var')
+        fprintf(fid,'        <description>Element data: %g</description>\n',edata(i));
+        imap = fix((edata(i)-cmin)/(cmax-cmin)*size(cmap,1))+1;
+        if     (imap >= 1) && (imap <= size(cmap,1))
+            fprintf(fid,'        <styleUrl>#MatlabColor%d</styleUrl>\n',imap);
+        elseif (edata(i) == cmax)
+            fprintf(fid,'        <styleUrl>#MatlabColor%d</styleUrl>\n',size(cmap,1));
+        else
+            fprintf(fid,'        <styleUrl>#BlackLineEmptyPoly</styleUrl>\n');
+        end
+    else
+        fprintf(fid,'        <styleUrl>#BlackLineRandomPoly</styleUrl>\n');
+    end
+    fprintf(fid,'        <Polygon>\n');
+    fprintf(fid,'          <extrude>1</extrude>\n');
+    fprintf(fid,'          <altitudeMode>relativeToGround</altitudeMode>\n');
+    fprintf(fid,'          <outerBoundaryIs>\n');
+    fprintf(fid,'            <LinearRing>\n');
+    fprintf(fid,'              <coordinates>\n');
+    for j=1:size(md.elements,2)
+        [lat(j),long(j)]=mapxy(md.x(md.elements(i,j)),md.y(md.elements(i,j)),'s');
+        fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(j),lat(j),alt);
+    end
+    fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(1),lat(1),alt);
+
+    fprintf(fid,'              </coordinates>\n');
+    fprintf(fid,'            </LinearRing>\n');
+    fprintf(fid,'          </outerBoundaryIs>\n');
+    fprintf(fid,'        </Polygon>\n');
+    fprintf(fid,'      </Placemark>\n');
+end
+fprintf(fid,'    </Folder>\n');
+
+end
+
Index: /issm/trunk/src/m/kml/kml_mesh_write.m
===================================================================
--- /issm/trunk/src/m/kml/kml_mesh_write.m	(revision 6416)
+++ /issm/trunk/src/m/kml/kml_mesh_write.m	(revision 6417)
@@ -1,3 +1,3 @@
-%
+%%
 %  write a kml file of the mesh from the model.
 %
@@ -248,1510 +248,10 @@
 fprintf(fid,'</kml>\n');
 
-fclose(fid);
-display('End of file successfully written.');
-
-end
-
-%%
-%  create kml polygons for the element mesh.
-%
-%  []=kml_mesh_elem(fid,md,params)
-%
-%  where the required input is:
-%    fid           (numeric, file ID of .kml file)
-%    md            (model, model class object)
-%
-%  the optional input is:
-%    params        (string/numeric, parameter names and values)
-%
-%  and the optional input is:
-%    data          (numeric, element or nodal results data)
-%    alt           (numeric, altitude for polygons, default 10000)
-%    cmin          (numeric, minimum of color map)
-%    cmax          (numeric, maximum of color map)
-%    cmap          (char or numeric, colormap definition)
-%
-function []=kml_mesh_elem(varargin)
-
-if ~nargin
-    help kml_mesh_elem
-    return
-end
-
-%%  process input data
-
-iarg=1;
-if (nargin >= 1)
-    fid=varargin{1};
-end
-if ~exist('fid','var') || isempty(fid) || (fid < 0)
-    error('File ID ''%d'' must be open.',fid);
-end
-
-iarg=iarg+1;
-if (nargin >= 2)
-    md=varargin{2};
-end
-if ~exist('md','var') || isempty(md) || ~isa(md,'model')
-    error(['Model ''' inputname(iarg) ''' is unrecognized class ''' class(md) '''.']);
-end
-
-%  parameters
-
-iarg=iarg+1;
-while (iarg <= nargin-1)
-    if ischar(varargin{iarg})
-        eval([varargin{iarg} '=varargin{iarg+1};']);
-        if (numel(varargin{iarg+1}) <= 20)
-            disp([varargin{iarg} '=' any2str(varargin{iarg+1}) ';']);
-        else
-            disp([varargin{iarg} '=' string_size(varargin{iarg+1}) ' ' class(varargin{iarg+1}) ';']);
-        end
-        if strcmpi(varargin{iarg},'data')
-            cdata=inputname(iarg+1);
-        end
-    else
-        error(['''' any2str(varargin{iarg}) ''' is not a parameter name.']);
-    end
-    iarg=iarg+2;
-end
-
-if exist('data','var') && ~isempty(data)
-    if     (numel(data)==md.numberofelements)
-        edata=data;
-    elseif (numel(data)==md.numberofgrids)
-        ndata=data;
-        display('Averaging nodal data to element data.');
-        edata=zeros(1,md.numberofelements);
-        for i=1:size(md.elements,1)
-            for j=1:size(md.elements,2)
-                edata(i)=edata(i)+ndata(md.elements(i,j));
-            end
-            edata(i)=edata(i)/size(md.elements,2);
-        end
-    else
-        error(['Data has incorrect number of ' num2str(numel(data)) ' values.']);
-    end
-end
-
-%  colormap command operates on a figure, so create an invisible one
-%  (could also directly call colormaps, e.g. jet(64), but risky)
-
-hfig=figure('Visible','off');
-if exist('cmap','var')
-    colormap(cmap)
-end
-cmap=colormap;
-close(hfig)
-    
-if exist('edata','var')
-    if ~exist('cmin','var')
-        cmin=min(min(edata));
-    end
-    if ~exist('cmax','var')
-        cmax=max(max(edata));
-    end
-end
-
-if ~exist('alt','var')
-    alt=10000;
-end
-
-%%  write folder for mesh
-
-fprintf(fid,'    <Folder>\n');
-if exist('cdata','var') && ~isempty(cdata)
-    fprintf(fid,'      <name>Data: %s</name>\n',cdata);
+if (fclose(fid) < 0)
+    error('File ''%s'' could not be closed.',filek2);
 else
-    fprintf(fid,'      <name>Mesh</name>\n');
-end
-fprintf(fid,'      <visibility>1</visibility>\n');
-fprintf(fid,'      <description>Elements=%d, Grids=%d</description>\n',...
-    md.numberofelements,md.numberofgrids);
-
-%  write each element as a polygon
-
-disp(['Writing ' num2str(size(md.elements,1)) ' tria elements as KML polygons.']);
-for i=1:size(md.elements,1)
-    fprintf(fid,'      <Placemark>\n');
-    fprintf(fid,'        <name>Element %d</name>\n',i);
-    fprintf(fid,'        <visibility>1</visibility>\n');
-    if exist('edata','var')
-        fprintf(fid,'        <description>Element data: %g</description>\n',edata(i));
-        imap = fix((edata(i)-cmin)/(cmax-cmin)*size(cmap,1))+1;
-        if     (imap >= 1) && (imap <= size(cmap,1))
-            fprintf(fid,'        <styleUrl>#MatlabColor%d</styleUrl>\n',imap);
-        elseif (edata(i) == cmax)
-            fprintf(fid,'        <styleUrl>#MatlabColor%d</styleUrl>\n',size(cmap,1));
-        else
-            fprintf(fid,'        <styleUrl>#BlackLineEmptyPoly</styleUrl>\n');
-        end
-    else
-        fprintf(fid,'        <styleUrl>#BlackLineRandomPoly</styleUrl>\n');
-    end
-    fprintf(fid,'        <Polygon>\n');
-    fprintf(fid,'          <extrude>1</extrude>\n');
-    fprintf(fid,'          <altitudeMode>relativeToGround</altitudeMode>\n');
-    fprintf(fid,'          <outerBoundaryIs>\n');
-    fprintf(fid,'            <LinearRing>\n');
-    fprintf(fid,'              <coordinates>\n');
-    for j=1:size(md.elements,2)
-        [lat(j),long(j)]=mapxy(md.x(md.elements(i,j)),md.y(md.elements(i,j)),'s');
-        fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(j),lat(j),alt);
-    end
-    fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(1),lat(1),alt);
-
-    fprintf(fid,'              </coordinates>\n');
-    fprintf(fid,'            </LinearRing>\n');
-    fprintf(fid,'          </outerBoundaryIs>\n');
-    fprintf(fid,'        </Polygon>\n');
-    fprintf(fid,'      </Placemark>\n');
-end
-fprintf(fid,'    </Folder>\n');
-
-end
-
-%%
-%  create kml linestrings for the flagged partition edges.
-%
-%  []=kml_part_flagedges(fid,md,params)
-%
-%  where the required input is:
-%    fid           (char, file ID of .kml file)
-%    md            (model, model class object)
-%
-%  the optional input is:
-%    params        (string/numeric, parameter names and values)
-%
-%  and the optional input is:
-%    alt           (numeric, altitude for polygons, default 10000)
-%    prtplt        (char, 'off'/'no' for partition segment plot)
-%
-function []=kml_part_flagedges(varargin)
-
-if ~nargin
-    help kml_part_flagedges
-    return
-end
-
-%%  process input data
-
-iarg=1;
-if (nargin >= 1)
-    fid=varargin{1};
-end
-if ~exist('fid','var') || isempty(fid) || (fid < 0)
-    error('File ID ''%d'' must be open.',fid);
-end
-
-iarg=iarg+1;
-if (nargin >= 2)
-    md=varargin{2};
-end
-if ~exist('md','var') || isempty(md) || ~isa(md,'model')
-    error(['Model ''' inputname(iarg) ''' is unrecognized class ''' class(md) '''.']);
-end
-
-%  parameters
-
-iarg=iarg+1;
-while (iarg <= nargin-1)
-    if ischar(varargin{iarg})
-        eval([varargin{iarg} '=varargin{iarg+1};']);
-        if (numel(varargin{iarg+1}) <= 20)
-            disp([varargin{iarg} '=' any2str(varargin{iarg+1}) ';']);
-        else
-            disp([varargin{iarg} '=' string_size(varargin{iarg+1}) ' ' class(varargin{iarg+1}) ';']);
-        end
-    else
-        error(['''' any2str(varargin{iarg}) ''' is not a parameter name.']);
-    end
-    iarg=iarg+2;
-end
-
-if ~exist('alt','var')
-    alt=10000;
-end
-
-%%  write folder for partition segments
-
-if (~exist('prtplt','var') || strncmpi(prtplt,'on' ,2) || strncmpi(prtplt,'y',1)) && ...
-    md.npart
-    [xseg,yseg]=flagedges(md.elements,md.x,md.y,md.part);
-    fprintf(fid,'    <Folder>\n');
-    fprintf(fid,'      <name>Partition Segments</name>\n');
-    fprintf(fid,'      <visibility>1</visibility>\n');
-    fprintf(fid,'      <description>Partitions=%d, Segments=%d</description>\n',...
-        md.npart,size(xseg,1));
-
-%  write each segment as a linestring
-
-    disp(['Writing ' num2str(size(xseg,1)) ' partition segments as KML linestrings.']);
-    for i=1:size(xseg,1)
-        fprintf(fid,'      <Placemark>\n');
-        fprintf(fid,'        <name>Segment %d</name>\n',i);
-        fprintf(fid,'        <visibility>1</visibility>\n');
-        fprintf(fid,'        <styleUrl>#RedLineRedPoly</styleUrl>\n');
-        fprintf(fid,'        <LineString>\n');
-        fprintf(fid,'          <extrude>1</extrude>\n');
-        fprintf(fid,'          <tessellate>1</tessellate>\n');
-        fprintf(fid,'          <altitudeMode>relativeToGround</altitudeMode>\n');
-        fprintf(fid,'          <coordinates>\n');
-        for j=1:2
-            [lat(j),long(j)]=mapxy(xseg(i,j),yseg(i,j),'s');
-            fprintf(fid,'            %0.16g,%0.16g,%0.16g\n',long(j),lat(j),alt);
-        end
-
-        fprintf(fid,'          </coordinates>\n');
-        fprintf(fid,'        </LineString>\n');
-        fprintf(fid,'      </Placemark>\n');
-    end
-    fprintf(fid,'    </Folder>\n');
-end
-
-end
-
-%%
-%  create kml linestrings for the unshared element edges.
-%
-%  []=kml_unsh_edges(fid,md,params)
-%
-%  where the required input is:
-%    fid           (char, file ID of .kml file)
-%    md            (model, model class object)
-%
-%  the optional input is:
-%    params        (string/numeric, parameter names and values)
-%
-%  and the optional input is:
-%    alt           (numeric, altitude for polygons, default 10000)
-%    prtplt        (char, 'off'/'no' for partition segment plot)
-%
-function []=kml_unsh_edges(varargin)
-
-if ~nargin
-    help kml_unsh_edges
-    return
-end
-
-%%  process input data
-
-iarg=1;
-if (nargin >= 1)
-    fid=varargin{1};
-end
-if ~exist('fid','var') || isempty(fid) || (fid < 0)
-    error('File ID ''%d'' must be open.',fid);
-end
-
-iarg=iarg+1;
-if (nargin >= 2)
-    md=varargin{2};
-end
-if ~exist('md','var') || isempty(md) || ~isa(md,'model')
-    error(['Model ''' inputname(iarg) ''' is unrecognized class ''' class(md) '''.']);
-end
-
-%  parameters
-
-iarg=iarg+1;
-while (iarg <= nargin-1)
-    if ischar(varargin{iarg})
-        eval([varargin{iarg} '=varargin{iarg+1};']);
-        if (numel(varargin{iarg+1}) <= 20)
-            disp([varargin{iarg} '=' any2str(varargin{iarg+1}) ';']);
-        else
-            disp([varargin{iarg} '=' string_size(varargin{iarg+1}) ' ' class(varargin{iarg+1}) ';']);
-        end
-    else
-        error(['''' any2str(varargin{iarg}) ''' is not a parameter name.']);
-    end
-    iarg=iarg+2;
-end
-
-if ~exist('alt','var')
-    alt=10000;
-end
-
-%%  write folder for unshared edges
-
-if (~exist('prtplt','var') || strncmpi(prtplt,'on' ,2) || strncmpi(prtplt,'y',1)) && ...
-    md.npart
-    [edgeadj]=edgeadjacency(md.elements,md.nodeconnectivity);
-    [icol,irow]=find(edgeadj'==0);
-    edgeuns=zeros(length(irow),2);
-    for i=1:length(irow)
-        edgeuns(i,1)=md.elements(irow(i),icol(i));
-        edgeuns(i,2)=md.elements(irow(i),mod(icol(i),size(md.elements,2))+1);
-    end
-    fprintf(fid,'    <Folder>\n');
-    fprintf(fid,'      <name>Unshared Edges</name>\n');
-    fprintf(fid,'      <visibility>1</visibility>\n');
-    fprintf(fid,'      <description>Partitions=%d, Edges=%d</description>\n',...
-        md.npart,size(edgeuns,1));
-
-%  write each edge as a linestring
-
-    disp(['Writing ' num2str(size(edgeuns,1)) ' unshared edges as KML linestrings.']);
-    for i=1:size(edgeuns,1)
-        fprintf(fid,'      <Placemark>\n');
-        fprintf(fid,'        <name>Edge %d</name>\n',i);
-        fprintf(fid,'        <visibility>1</visibility>\n');
-        fprintf(fid,'        <styleUrl>#RedLineRedPoly</styleUrl>\n');
-        fprintf(fid,'        <LineString>\n');
-        fprintf(fid,'          <extrude>1</extrude>\n');
-        fprintf(fid,'          <tessellate>1</tessellate>\n');
-        fprintf(fid,'          <altitudeMode>relativeToGround</altitudeMode>\n');
-        fprintf(fid,'          <coordinates>\n');
-        for j=1:2
-            [lat(j),long(j)]=mapxy(md.x(edgeuns(i,j)),md.y(edgeuns(i,j)),'s');
-            fprintf(fid,'            %0.16g,%0.16g,%0.16g\n',long(j),lat(j),alt);
-        end
-
-        fprintf(fid,'          </coordinates>\n');
-        fprintf(fid,'        </LineString>\n');
-        fprintf(fid,'      </Placemark>\n');
-    end
-    fprintf(fid,'    </Folder>\n');
-end
-
-end
-
-%%
-%  create kml polygons for the partition elements.
-%
-%  []=kml_part_elems(fid,md,params)
-%
-%  where the required input is:
-%    fid           (char, file ID of .kml file)
-%    md            (model, model class object)
-%
-%  the optional input is:
-%    params        (string/numeric, parameter names and values)
-%
-%  and the optional input is:
-%    data          (numeric, element or nodal results data)
-%    alt           (numeric, altitude for polygons, default 10000)
-%    cmin          (numeric, minimum of color map)
-%    cmax          (numeric, maximum of color map)
-%    cmap          (char or numeric, colormap definition)
-%    prtplt        (char, 'off'/'no' for partition segment plot)
-%
-function []=kml_part_elems(varargin)
-
-if ~nargin
-    help kml_part_elems
-    return
-end
-
-%%  process input data
-
-iarg=1;
-if (nargin >= 1)
-    fid=varargin{1};
-end
-if ~exist('fid','var') || isempty(fid) || (fid < 0)
-    error('File ID ''%d'' must be open.',fid);
-end
-
-iarg=iarg+1;
-if (nargin >= 2)
-    md=varargin{2};
-end
-if ~exist('md','var') || isempty(md) || ~isa(md,'model')
-    error(['Model ''' inputname(iarg) ''' is unrecognized class ''' class(md) '''.']);
-end
-
-%  parameters
-
-iarg=iarg+1;
-while (iarg <= nargin-1)
-    if ischar(varargin{iarg})
-        eval([varargin{iarg} '=varargin{iarg+1};']);
-        if (numel(varargin{iarg+1}) <= 20)
-            disp([varargin{iarg} '=' any2str(varargin{iarg+1}) ';']);
-        else
-            disp([varargin{iarg} '=' string_size(varargin{iarg+1}) ' ' class(varargin{iarg+1}) ';']);
-        end
-        if strcmpi(varargin{iarg},'data')
-            cdata=inputname(iarg+1);
-        end
-    else
-        error(['''' any2str(varargin{iarg}) ''' is not a parameter name.']);
-    end
-    iarg=iarg+2;
-end
-
-if exist('data','var') && ~isempty(data)
-    if     (numel(data)==md.numberofelements)
-        edata=data;
-    elseif (numel(data)==md.numberofgrids)
-        ndata=data;
-        display('Averaging nodal data to element data.');
-        edata=zeros(1,md.numberofelements);
-        for i=1:size(md.elements,1)
-            for j=1:size(md.elements,2)
-                edata(i)=edata(i)+ndata(md.elements(i,j));
-            end
-            edata(i)=edata(i)/size(md.elements,2);
-        end
-    else
-        error(['Data has incorrect number of ' num2str(numel(data)) ' values.']);
-    end
-end
-
-%  colormap command operates on a figure, so create an invisible one
-%  (could also directly call colormaps, e.g. jet(64), but risky)
-
-hfig=figure('Visible','off');
-if exist('cmap','var')
-    colormap(cmap)
-end
-cmap=colormap;
-close(hfig)
-    
-if exist('edata','var')
-    if ~exist('cmin','var')
-        cmin=min(min(edata));
-    end
-    if ~exist('cmax','var')
-        cmax=max(max(edata));
-    end
-end
-
-if ~exist('alt','var')
-    alt=10000;
-end
-
-%  write folder for partition elements
-
-if (~exist('prtplt','var') || strncmpi(prtplt,'on' ,2) || strncmpi(prtplt,'y',1)) && ...
-    md.npart
-    fprintf(fid,'    <Folder>\n');
-    fprintf(fid,'      <name>Partition Elements</name>\n');
-    fprintf(fid,'      <visibility>1</visibility>\n');
-    fprintf(fid,'      <description>Partitions=%d, Nodes=%d</description>\n',...
-        md.npart,md.numberofgrids);
-
-%  write each partition loop as a polygon
-
-    disp(['Writing ' num2str(md.npart) ' partitions as KML polygons.']);
-    epart=md.part(md.elements)+1;
-    if exist('ndata','var') || exist('edata','var')
-        pdata=zeros(1,md.npart);
-        pdata(:)=NaN;
-    end
-    for k=1:md.npart
-        [icol,irow]=find(epart'==k);
-        irow=unique(irow);
-        elem=md.elements(irow,:);
-        nodecon=nodeconnectivity(elem,md.numberofgrids);
-        [edgeper,elemper,iloop]=edgeperimeter(elem,nodecon);
-        iloop(end+1)=size(edgeper,1)+1;
-        if exist('ndata','var')
-            pdata(k)=ndata(find(md.part+1==k,1));
-        elseif exist('edata','var')
-            for i=1:size(epart,1)
-                if isempty(find(epart(i,:)~=k,1))
-                    pdata(k)=edata(i);
-                    break
-                end
-            end
-            if isnan(pdata(k))
-                warning('Data for Partition %d is not defined.\n',k)
-            end
-        end
-        for i=1:length(iloop)-1
-            fprintf(fid,'      <Placemark>\n');
-            if (length(iloop)-1 > 1)
-                fprintf(fid,'        <name>Partition %d, Loop %d</name>\n',k,i);
-            else
-                fprintf(fid,'        <name>Partition %d</name>\n',k);
-            end
-            fprintf(fid,'        <visibility>1</visibility>\n');
-            if exist('pdata','var')
-                fprintf(fid,'        <description>Partition data: %g</description>\n',pdata(k));
-                imap = fix((pdata(k)-cmin)/(cmax-cmin)*size(cmap,1))+1;
-                if     (imap >= 1) && (imap <= size(cmap,1))
-                    fprintf(fid,'        <styleUrl>#MatlabColor%d</styleUrl>\n',imap);
-                elseif (pdata(k) == cmax)
-                    fprintf(fid,'        <styleUrl>#MatlabColor%d</styleUrl>\n',size(cmap,1));
-                else
-                    fprintf(fid,'        <styleUrl>#BlackLineEmptyPoly</styleUrl>\n');
-                end
-            else
-                fprintf(fid,'        <styleUrl>#BlackLineRandomPoly</styleUrl>\n');
-            end
-            fprintf(fid,'        <Polygon>\n');
-            fprintf(fid,'          <extrude>1</extrude>\n');
-            fprintf(fid,'          <altitudeMode>relativeToGround</altitudeMode>\n');
-            fprintf(fid,'          <outerBoundaryIs>\n');
-            fprintf(fid,'            <LinearRing>\n');
-            fprintf(fid,'              <coordinates>\n');
-            for j=iloop(i):iloop(i+1)-1
-                [lat(j),long(j)]=mapxy(md.x(edgeper(j,1)),md.y(edgeper(j,1)),'s');
-                fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(j),lat(j),alt);
-            end
-            fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(iloop(i)),lat(iloop(i)),alt);
-
-            fprintf(fid,'              </coordinates>\n');
-            fprintf(fid,'            </LinearRing>\n');
-            fprintf(fid,'          </outerBoundaryIs>\n');
-            fprintf(fid,'        </Polygon>\n');
-            fprintf(fid,'      </Placemark>\n');
-        end
-    end
-    fprintf(fid,'    </Folder>\n');
-end
-
-end
-
-%%
-%  create kml linestrings for the partition edges.
-%
-%  []=kml_part_edges(fid,md,params)
-%
-%  where the required input is:
-%    fid           (char, file ID of .kml file)
-%    md            (model, model class object)
-%
-%  the optional input is:
-%    params        (string/numeric, parameter names and values)
-%
-%  and the optional input is:
-%    data          (numeric, element or nodal results data)
-%    alt           (numeric, altitude for polygons, default 10000)
-%    cmin          (numeric, minimum of color map)
-%    cmax          (numeric, maximum of color map)
-%    cmap          (char or numeric, colormap definition)
-%    prtplt        (char, 'off'/'no' for partition segment plot)
-%
-function []=kml_part_edges(varargin)
-
-if ~nargin
-    help kml_part_edges
-    return
-end
-
-%%  process input data
-
-iarg=1;
-if (nargin >= 1)
-    fid=varargin{1};
-end
-if ~exist('fid','var') || isempty(fid) || (fid < 0)
-    error('File ID ''%d'' must be open.',fid);
-end
-
-iarg=iarg+1;
-if (nargin >= 2)
-    md=varargin{2};
-end
-if ~exist('md','var') || isempty(md) || ~isa(md,'model')
-    error(['Model ''' inputname(iarg) ''' is unrecognized class ''' class(md) '''.']);
-end
-
-%  parameters
-
-iarg=iarg+1;
-while (iarg <= nargin-1)
-    if ischar(varargin{iarg})
-        eval([varargin{iarg} '=varargin{iarg+1};']);
-        if (numel(varargin{iarg+1}) <= 20)
-            disp([varargin{iarg} '=' any2str(varargin{iarg+1}) ';']);
-        else
-            disp([varargin{iarg} '=' string_size(varargin{iarg+1}) ' ' class(varargin{iarg+1}) ';']);
-        end
-        if strcmpi(varargin{iarg},'data')
-            cdata=inputname(iarg+1);
-        end
-    else
-        error(['''' any2str(varargin{iarg}) ''' is not a parameter name.']);
-    end
-    iarg=iarg+2;
-end
-
-if exist('data','var') && ~isempty(data)
-    if     (numel(data)==md.numberofelements)
-        edata=data;
-    elseif (numel(data)==md.numberofgrids)
-        ndata=data;
-        display('Averaging nodal data to element data.');
-        edata=zeros(1,md.numberofelements);
-        for i=1:size(md.elements,1)
-            for j=1:size(md.elements,2)
-                edata(i)=edata(i)+ndata(md.elements(i,j));
-            end
-            edata(i)=edata(i)/size(md.elements,2);
-        end
-    else
-        error(['Data has incorrect number of ' num2str(numel(data)) ' values.']);
-    end
-end
-
-%  colormap command operates on a figure, so create an invisible one
-%  (could also directly call colormaps, e.g. jet(64), but risky)
-
-hfig=figure('Visible','off');
-if exist('cmap','var')
-    colormap(cmap)
-end
-cmap=colormap;
-close(hfig)
-    
-if exist('edata','var')
-    if ~exist('cmin','var')
-        cmin=min(min(edata));
-    end
-    if ~exist('cmax','var')
-        cmax=max(max(edata));
-    end
-end
-
-if ~exist('alt','var')
-    alt=10000;
-end
-
-%%  write folder for partition edges
-
-if (~exist('prtplt','var') || strncmpi(prtplt,'on' ,2) || strncmpi(prtplt,'y',1)) && ...
-    md.npart
-    fprintf(fid,'    <Folder>\n');
-    fprintf(fid,'      <name>Partition Edges</name>\n');
-    fprintf(fid,'      <visibility>1</visibility>\n');
-    fprintf(fid,'      <description>Partitions=%d, Nodes=%d</description>\n',...
-        md.npart,md.numberofgrids);
-
-%  write each partition loop as linestrings
-
-    disp(['Writing ' num2str(md.npart) ' partitions as KML linestrings.']);
-    epart=md.part(md.elements)+1;
-    if exist('ndata','var') || exist('edata','var')
-        pdata=zeros(1,md.npart);
-        pdata(:)=NaN;
-    end
-
-%  loop over each partition
-
-    for k=1:md.npart
-        disp(['partition k=' int2str(k)])
-        
-%  for each partition, find all the included elements and determine the
-%  perimeter (including those shared by another partition)
-
-        [icol,irow]=find(epart'==k);
-        irow=unique(irow);
-        elemp=md.elements(irow,:);
-        epartp=epart(irow,:);
-        nodeconp=nodeconnectivity(elemp,md.numberofgrids);
-        [edgeadjp]=edgeadjacency(elemp,nodeconp);
-        [edgeper,elemper,iloop]=edgeperimeter(elemp,nodeconp,edgeadjp);
-        iloop(end+1)=size(edgeper,1)+1;
-        
-%  determine the data to be used for the colors (if any)
-
-        if exist('ndata','var')
-            pdata(k)=ndata(find(md.part+1==k,1));
-        elseif exist('edata','var')
-            for i=1:size(epartp,1)
-                if isempty(find(epart(i,:)~=k,1))
-                    pdata(k)=edata(irow(i));
-                    break
-                end
-            end
-            if isnan(pdata(k))
-                warning('Data for Partition %d is not defined.\n',k)
-            end
-        end
-        
-%  loop over each loop of the perimeter for the given partition
-
-        for i=1:length(iloop)-1
-            fprintf(fid,'      <Placemark>\n');
-            if (length(iloop)-1 > 1)
-                fprintf(fid,'        <name>Partition %d, Loop %d</name>\n',k,i);
-            else
-                fprintf(fid,'        <name>Partition %d</name>\n',k);
-            end
-            fprintf(fid,'        <visibility>1</visibility>\n');
-            if exist('pdata','var')
-                fprintf(fid,'        <description>Partition data: %g</description>\n',pdata(k));
-                imap = fix((pdata(k)-cmin)/(cmax-cmin)*size(cmap,1))+1;
-                if     (imap >= 1) && (imap <= size(cmap,1))
-                    fprintf(fid,'        <styleUrl>#MatlabColor%d</styleUrl>\n',imap);
-                elseif (pdata(k) == cmax)
-                    fprintf(fid,'        <styleUrl>#MatlabColor%d</styleUrl>\n',size(cmap,1));
-                else
-                    fprintf(fid,'        <styleUrl>#BlackLineEmptyPoly</styleUrl>\n');
-                end
-            else
-                fprintf(fid,'        <styleUrl>#BlackLineRandomPoly</styleUrl>\n');
-            end
-%            fprintf(fid,'        <Polygon>\n');
-%            fprintf(fid,'          <extrude>1</extrude>\n');
-%            fprintf(fid,'          <altitudeMode>relativeToGround</altitudeMode>\n');
-%            fprintf(fid,'          <outerBoundaryIs>\n');
-%            fprintf(fid,'            <LinearRing>\n');
-        fprintf(fid,'        <LineString>\n');
-        fprintf(fid,'          <extrude>1</extrude>\n');
-        fprintf(fid,'          <tessellate>1</tessellate>\n');
-        fprintf(fid,'          <altitudeMode>relativeToGround</altitudeMode>\n');
-            fprintf(fid,'              <coordinates>\n');
-            elast=0;
-            nlast=0;
-            slast=0;
-            lat=[];
-            long=[];
-            
-%  loop over the element edges on the loop of the partition
-
-            j=iloop(i);
-            while (j < iloop(i+1))
-%  find which side of element is referenced in perimeter list
-                for l=1:size(elemp,2)
-                    if ((elemp(elemper(j),l)          == edgeper(j,1)) && ...
-                        (elemp(elemper(j),mod(l,3)+1) == edgeper(j,2))) || ...
-                       ((elemp(elemper(j),l)          == edgeper(j,2)) && ...
-                        (elemp(elemper(j),mod(l,3)+1) == edgeper(j,1)))
-                        jedge=l;
-                        break
-                    end
-                end
-
-%  check if element side connects nodes in partition
-                if (epartp(elemper(j),jedge)          == k) && ...
-                   (epartp(elemper(j),mod(jedge,3)+1) == k)
-%  write out specified element side
-                    disp(['segment j=' int2str(j) ' unshared edge ' int2str(edgeper(j,1)) ' to ' int2str(edgeper(j,2)) ' on side ' int2str(jedge) ' from element ' int2str(elemper(j)) ' written.'])
-%  if first edge, write out first node
-                    if ~elast
-                        [lat(end+1),long(end+1)]=mapxy(md.x(edgeper(j,1)),md.y(edgeper(j,1)),'s');
-                        fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
-                    end
-                    [lat(end+1),long(end+1)]=mapxy(md.x(edgeper(j,2)),md.y(edgeper(j,2)),'s');
-                    fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
-                    elast=elemper(j);
-                    nlast=edgeper(j,2);
-                    slast=0;
-                    j=j+1;
-                    
-%  element not entirely within partition, so figure out boundary
-                else
-                    disp(['segment j=' int2str(j) ' from element ' int2str(elemper(j)) ' shared by other partitions.'])
-                    ielem=elemper(j);
-
-%  follow partition boundary through elements not wholly in partition
-%  (may include elements not in perimeter list)
-
-                    while 1
-%  if first edge, figure out direction from perimeter edge direction
-                        if ~nlast && ~slast
-                            nlast=find(elemp(ielem,:)==edgeper(j,1));
-                            nnext=find(elemp(ielem,:)==edgeper(j,2));
-                            if     (nlast+nnext == 3)
-                                slast=1;
-                            elseif (nlast+nnext == 5)
-                                slast=2;
-                            elseif (nlast+nnext == 4)
-                                slast=3;
-                            end
-                            if     (nnext+(6-nlast-nnext) == 3)
-                                snext=1;
-                            elseif (nnext+(6-nlast-nnext) == 5)
-                                snext=2;
-                            elseif (nnext+(6-nlast-nnext) == 4)
-                                snext=3;
-                            end
-
-%  find how many nodes of current element are in current partition
-%  (1 or 2, not 3) and handle each permutation separately
-                            ipart=find(epartp(ielem,:)==k);
-%  two nodes are in current partition, so cut off other node
-                            if (length(ipart) == 2)
-                                switch 6-sum(ipart)
-                                    case nlast
-                                        slast=6-slast-snext;
-                                        nlast=0;
-                                    case nnext
-                                        if (epartp(ielem,nnext) == k)
-                                            nlast=nnext;
-                                        end
-                                    otherwise
-                                        slast=6-slast-snext;
-                                        nlast=0;
-                                end
-%  one node is in current partition
-                            else
-%  all different, so cut through centroid
-                                if (epartp(ielem,1) ~= epartp(ielem,2)) && ...
-                                   (epartp(ielem,2) ~= epartp(ielem,3)) && ...
-                                   (epartp(ielem,3) ~= epartp(ielem,1))
-                                    switch ipart
-                                        case {nlast,nnext}
-                                            if (epartp(ielem,nnext) == k)
-                                                nlast=nnext;
-                                            end
-                                        otherwise
-                                            slast=6-slast-snext;
-                                            nlast=0;
-                                    end
-%  other two are in the same partition, so cut them off
-                                else
-                                    switch ipart
-                                        case nlast
-                                            if (epartp(ielem,nnext) == k)
-                                                nlast=nnext;
-                                            end
-                                        case nnext
-                                            slast=snext;
-                                            nlast=0;
-                                        otherwise
-                                            slast=6-slast-snext;
-                                            nlast=0;
-                                    end
-                                end
-                            end
-
-%  last edge exited last element at node
-                            if nlast
-%  write out first node of first side for half-edge to midpoint
-                                disp(['segment j=' int2str(j) ' unshared half edge from node ' int2str(elemp(ielem,nlast)) ' (node ' int2str(nlast) ') on side ' int2str(slast) ' from element ' int2str(ielem) ' written.'])
-                                [lat(end+1),long(end+1)]=mapxy((md.x(elemp(ielem,nlast))),...
-                                                               (md.y(elemp(ielem,nlast))),'s');
-                                fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
-                            end
-                            nlast=0;
-                            
-%  write out midpoint of first side
-                            [lat(end+1),long(end+1)]=mapxy((md.x(elemp(ielem,slast))...
-                                                           +md.x(elemp(ielem,mod(slast,3)+1)))/2.,...
-                                                           (md.y(elemp(ielem,slast))...
-                                                           +md.y(elemp(ielem,mod(slast,3)+1)))/2.,'s');
-                            fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
-                        end
-
-%  last edge exited last element at node
-                        if nlast
-                            if elast
-%  find where last node on previous element occurs on current element
-                                nlast=find(elemp(ielem,:)==nlast,1);
-                            end
-%  half-edge occurs on unshared side from current node (unique unless mesh
-%  is only attached at node)
-                            switch nlast
-                                case 1
-                                    if ~edgeadjp(ielem,1)
-                                        nnext=2;
-                                        slast=1;
-                                    else
-                                        nnext=3;
-                                        slast=3;
-                                    end
-                                case 2
-                                    if ~edgeadjp(ielem,2)
-                                        nnext=3;
-                                        slast=2;
-                                    else
-                                        nnext=1;
-                                        slast=1;
-                                    end
-                                case 3
-                                    if ~edgeadjp(ielem,3)
-                                        nnext=1;
-                                        slast=3;
-                                    else
-                                        nnext=2;
-                                        slast=2;
-                                    end
-                            end
-%  write out half-edge from current node to midpoint of unshared side
-                            disp(['segment j=' int2str(j) ' unshared half edge from node ' int2str(elemp(ielem,nlast)) ' (node ' int2str(nlast) ') on side ' int2str(slast) ' from element ' int2str(ielem) ' written.'])
-                            [lat(end+1),long(end+1)]=mapxy((md.x(elemp(ielem,nlast))...
-                                                           +md.x(elemp(ielem,nnext)))/2.,...
-                                                           (md.y(elemp(ielem,nlast))...
-                                                           +md.y(elemp(ielem,nnext)))/2.,'s');
-                            fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
-                            nlast=0;
-
-%  last edge exited last element at midpoint of side
-                        elseif slast
-                            if elast
-%  find where last side on previous element occurs on current element
-                                slast=find(edgeadjp(ielem,:)==elast,1);
-                            end
-                        end
-
-%  find how many nodes of current element are in current partition
-%  (1 or 2, not 3) and handle each permutation separately
-                        ipart=find(epartp(ielem,:)==k);
-                        if (length(ipart) == 2)
-%  two nodes are in current partition, so cut off other node
-                            switch 6-sum(ipart)
-                                case 1
-                                    snext=3+1-slast;
-                                case 2
-                                    snext=1+2-slast;
-                                case 3
-                                    snext=2+3-slast;
-                            end
-                        else
-                            if (epartp(ielem,1) ~= epartp(ielem,2)) && ...
-                               (epartp(ielem,2) ~= epartp(ielem,3)) && ...
-                               (epartp(ielem,3) ~= epartp(ielem,1))
-%  all different, so cut through centroid
-                                disp(['element ielem=' int2str(ielem) ' centroid written.'])
-                                [lat(end+1),long(end+1)]=mapxy(sum(md.x(elemp(ielem,:)))/3.,...
-                                                               sum(md.y(elemp(ielem,:)))/3.,'s');
-                                fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
-                            end
-%  one node is in current partition, so cut off other two nodes
-                            switch ipart
-                                case 1
-                                    snext=3+1-slast;
-                                case 2
-                                    snext=1+2-slast;
-                                case 3
-                                    snext=2+3-slast;
-                            end
-                        end
-%  write out midpoint of opposite side
-                        disp(['segment j=' int2str(j) ' internal edge from side ' int2str(slast) ' to side ' int2str(snext) ' from element ' int2str(ielem) ' written.'])
-                        [lat(end+1),long(end+1)]=mapxy((md.x(elemp(ielem,snext))...
-                                                       +md.x(elemp(ielem,mod(snext,3)+1)))/2.,...
-                                                       (md.y(elemp(ielem,snext))...
-                                                       +md.y(elemp(ielem,mod(snext,3)+1)))/2.,'s');
-                        fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
-                        elast=ielem;
-                        nlast=0;
-                        slast=snext;
-%  find adjacent element to opposite side
-                        ielem=edgeadjp(elast,slast);
-%  if opposite side is unshared, find it in edge perimeter list
-                        if ~ielem
-                            jlast=find(elemper(j:end)==elast)+j-1;
-                            j=0;
-                            for l=1:length(jlast)
-                                if ((elemp(elast,slast)          == edgeper(jlast(l),1)) && ...
-                                    (elemp(elast,mod(slast,3)+1) == edgeper(jlast(l),2))) || ...
-                                   ((elemp(elast,slast)          == edgeper(jlast(l),2)) && ...
-                                    (elemp(elast,mod(slast,3)+1) == edgeper(jlast(l),1)))
-                                    j=jlast(l);
-                                    break
-                                end
-                            end
-                            if ~j
-                                j=iloop(i+1)-1;
-                            end
-%  write out half-edge from midpoint of unshared side to node
-                            if (epartp(elast,slast) == k)
-                                nnext=slast;
-                            else
-                                nnext=mod(slast,3)+1;
-                            end
-                            disp(['segment j=' int2str(j) ' unshared half edge on side ' int2str(slast) ' to node ' int2str(elemp(elast,nnext)) ' (node ' int2str(nnext) ') from element ' int2str(elast) ' written.'])
-                            [lat(end+1),long(end+1)]=mapxy(md.x(elemp(elast,nnext)),...
-                                                           md.y(elemp(elast,nnext)),'s');
-                            fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
-                            break
-%  if not unshared, advance perimeter list and watch for end
-                        else
-                            if (elast == elemper(j))
-                                if (j+1 < iloop(i+1)) && ...
-                                   ~isempty(find(elemper(j+1:end)~=elast,1))
-                                    j=j+find(elemper(j+1:end)~=elast,1);
-                                else
-                                    break
-                                end
-                            end
-                        end
-                    end
-                    j=j+1;
-                end
-            end
-%            fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(iloop(i)),lat(iloop(i)),alt);
-
-            fprintf(fid,'              </coordinates>\n');
-        fprintf(fid,'        </LineString>\n');
-%            fprintf(fid,'            </LinearRing>\n');
-%            fprintf(fid,'          </outerBoundaryIs>\n');
-%            fprintf(fid,'        </Polygon>\n');
-            fprintf(fid,'      </Placemark>\n');
-        end
-    end
-    fprintf(fid,'    </Folder>\n');
-end
-
-end
-
-%%
-%  create kml polygons for the partitions.
-%
-%  []=kml_partitions(fid,md,params)
-%
-%  where the required input is:
-%    fid           (char, file ID of .kml file)
-%    md            (model, model class object)
-%
-%  the optional input is:
-%    params        (string/numeric, parameter names and values)
-%
-%  and the optional input is:
-%    data          (numeric, element or nodal results data)
-%    alt           (numeric, altitude for polygons, default 10000)
-%    cmin          (numeric, minimum of color map)
-%    cmax          (numeric, maximum of color map)
-%    cmap          (char or numeric, colormap definition)
-%    prtplt        (char, 'off'/'no' for partition segment plot)
-%
-function []=kml_partitions(varargin)
-
-if ~nargin
-    help kml_part_edges
-    return
-end
-
-%%  process input data
-
-iarg=1;
-if (nargin >= 1)
-    fid=varargin{1};
-end
-if ~exist('fid','var') || isempty(fid) || (fid < 0)
-    error('File ID ''%d'' must be open.',fid);
-end
-
-iarg=iarg+1;
-if (nargin >= 2)
-    md=varargin{2};
-end
-if ~exist('md','var') || isempty(md) || ~isa(md,'model')
-    error(['Model ''' inputname(iarg) ''' is unrecognized class ''' class(md) '''.']);
-end
-
-%  parameters
-
-iarg=iarg+1;
-while (iarg <= nargin-1)
-    if ischar(varargin{iarg})
-        eval([varargin{iarg} '=varargin{iarg+1};']);
-        if (numel(varargin{iarg+1}) <= 20)
-            disp([varargin{iarg} '=' any2str(varargin{iarg+1}) ';']);
-        else
-            disp([varargin{iarg} '=' string_size(varargin{iarg+1}) ' ' class(varargin{iarg+1}) ';']);
-        end
-        if strcmpi(varargin{iarg},'data')
-            cdata=inputname(iarg+1);
-        end
-    else
-        error(['''' any2str(varargin{iarg}) ''' is not a parameter name.']);
-    end
-    iarg=iarg+2;
-end
-
-if exist('data','var') && ~isempty(data)
-    if     (numel(data)==md.numberofelements)
-        edata=data;
-    elseif (numel(data)==md.numberofgrids)
-        ndata=data;
-        display('Averaging nodal data to element data.');
-        edata=zeros(1,md.numberofelements);
-        for i=1:size(md.elements,1)
-            for j=1:size(md.elements,2)
-                edata(i)=edata(i)+ndata(md.elements(i,j));
-            end
-            edata(i)=edata(i)/size(md.elements,2);
-        end
-    else
-        error(['Data has incorrect number of ' num2str(numel(data)) ' values.']);
-    end
-end
-
-%  colormap command operates on a figure, so create an invisible one
-%  (could also directly call colormaps, e.g. jet(64), but risky)
-
-hfig=figure('Visible','off');
-if exist('cmap','var')
-    colormap(cmap)
-end
-cmap=colormap;
-close(hfig)
-    
-if exist('edata','var')
-    if ~exist('cmin','var')
-        cmin=min(min(edata));
-    end
-    if ~exist('cmax','var')
-        cmax=max(max(edata));
-    end
-end
-
-if ~exist('alt','var')
-    alt=10000;
-end
-
-%%  write folder for partitions
-
-if (~exist('prtplt','var') || strncmpi(prtplt,'on' ,2) || strncmpi(prtplt,'y',1)) && ...
-    md.npart
-    fprintf(fid,'    <Folder>\n');
-    fprintf(fid,'      <name>Partitions</name>\n');
-    fprintf(fid,'      <visibility>1</visibility>\n');
-    fprintf(fid,'      <description>Partitions=%d, Nodes=%d</description>\n',...
-        md.npart,md.numberofgrids);
-
-%  write each partition loop as linestrings
-
-    disp(['Writing ' num2str(md.npart) ' partitions as KML polygons.']);
-    epart=md.part(md.elements)+1;
-    if exist('ndata','var') || exist('edata','var')
-        pdata=zeros(1,md.npart);
-        pdata(:)=NaN;
-    end
-
-%  loop over each partition
-
-    for k=1:md.npart
-        disp(['partition k=' int2str(k)])
-        
-%  for each partition, find all the included elements and determine the
-%  perimeter (including those shared by another partition)
-
-        [icol,irow]=find(epart'==k);
-        irow=unique(irow);
-        elemp=md.elements(irow,:);
-        epartp=epart(irow,:);
-        nodeconp=nodeconnectivity(elemp,md.numberofgrids);
-        [edgeadjp]=edgeadjacency(elemp,nodeconp);
-        [edgeper,elemper,iloop]=edgeperimeter(elemp,nodeconp,edgeadjp);
-        iloop(end+1)=size(edgeper,1)+1;
-        
-%  determine the data to be used for the colors (if any)
-
-        if exist('ndata','var')
-            pdata(k)=ndata(find(md.part+1==k,1));
-        elseif exist('edata','var')
-            for i=1:size(epartp,1)
-                if isempty(find(epart(i,:)~=k,1))
-                    pdata(k)=edata(irow(i));
-                    break
-                end
-            end
-            if isnan(pdata(k))
-                warning('Data for Partition %d is not defined.\n',k)
-            end
-        end
-        
-%  loop over each loop of the perimeter for the given partition
-
-        for i=1:length(iloop)-1
-            fprintf(fid,'      <Placemark>\n');
-            if (length(iloop)-1 > 1)
-                fprintf(fid,'        <name>Partition %d, Loop %d</name>\n',k,i);
-            else
-                fprintf(fid,'        <name>Partition %d</name>\n',k);
-            end
-            fprintf(fid,'        <visibility>1</visibility>\n');
-            if exist('pdata','var')
-                fprintf(fid,'        <description>Partition data: %g</description>\n',pdata(k));
-                imap = fix((pdata(k)-cmin)/(cmax-cmin)*size(cmap,1))+1;
-                if     (imap >= 1) && (imap <= size(cmap,1))
-                    fprintf(fid,'        <styleUrl>#MatlabColor%d</styleUrl>\n',imap);
-                elseif (pdata(k) == cmax)
-                    fprintf(fid,'        <styleUrl>#MatlabColor%d</styleUrl>\n',size(cmap,1));
-                else
-                    fprintf(fid,'        <styleUrl>#BlackLineEmptyPoly</styleUrl>\n');
-                end
-            else
-                fprintf(fid,'        <styleUrl>#BlackLineRandomPoly</styleUrl>\n');
-            end
-            fprintf(fid,'        <Polygon>\n');
-            fprintf(fid,'          <extrude>1</extrude>\n');
-            fprintf(fid,'          <altitudeMode>relativeToGround</altitudeMode>\n');
-            fprintf(fid,'          <outerBoundaryIs>\n');
-            fprintf(fid,'            <LinearRing>\n');
-            fprintf(fid,'              <coordinates>\n');
-            elast=0;
-            nlast=0;
-            slast=0;
-            lat=[];
-            long=[];
-            
-%  loop over the element edges on the loop of the partition
-
-            j=iloop(i);
-            while (j < iloop(i+1))
-%  find which side of element is referenced in perimeter list
-                for l=1:size(elemp,2)
-                    if ((elemp(elemper(j),l)          == edgeper(j,1)) && ...
-                        (elemp(elemper(j),mod(l,3)+1) == edgeper(j,2))) || ...
-                       ((elemp(elemper(j),l)          == edgeper(j,2)) && ...
-                        (elemp(elemper(j),mod(l,3)+1) == edgeper(j,1)))
-                        jedge=l;
-                        break
-                    end
-                end
-
-%  check if element side connects nodes in partition
-                if (epartp(elemper(j),jedge)          == k) && ...
-                   (epartp(elemper(j),mod(jedge,3)+1) == k)
-%  write out specified element side
-                    disp(['segment j=' int2str(j) ' unshared edge ' int2str(edgeper(j,1)) ' to ' int2str(edgeper(j,2)) ' on side ' int2str(jedge) ' from element ' int2str(elemper(j)) ' written.'])
-%  if first edge, write out first node
-                    if ~elast
-                        [lat(end+1),long(end+1)]=mapxy(md.x(edgeper(j,1)),md.y(edgeper(j,1)),'s');
-                        fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
-                    end
-                    [lat(end+1),long(end+1)]=mapxy(md.x(edgeper(j,2)),md.y(edgeper(j,2)),'s');
-                    fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
-                    elast=elemper(j);
-                    nlast=edgeper(j,2);
-                    slast=0;
-                    j=j+1;
-                    
-%  element not entirely within partition, so figure out boundary
-                else
-                    disp(['segment j=' int2str(j) ' from element ' int2str(elemper(j)) ' shared by other partitions.'])
-                    ielem=elemper(j);
-
-%  follow partition boundary through elements not wholly in partition
-%  (may include elements not in perimeter list)
-
-                    while 1
-%  if first edge, figure out direction from perimeter edge direction
-                        if ~nlast && ~slast
-                            nlast=find(elemp(ielem,:)==edgeper(j,1));
-                            nnext=find(elemp(ielem,:)==edgeper(j,2));
-                            if     (nlast+nnext == 3)
-                                slast=1;
-                            elseif (nlast+nnext == 5)
-                                slast=2;
-                            elseif (nlast+nnext == 4)
-                                slast=3;
-                            end
-                            if     (nnext+(6-nlast-nnext) == 3)
-                                snext=1;
-                            elseif (nnext+(6-nlast-nnext) == 5)
-                                snext=2;
-                            elseif (nnext+(6-nlast-nnext) == 4)
-                                snext=3;
-                            end
-
-%  find how many nodes of current element are in current partition
-%  (1 or 2, not 3) and handle each permutation separately
-                            ipart=find(epartp(ielem,:)==k);
-%  two nodes are in current partition, so cut off other node
-                            if (length(ipart) == 2)
-                                switch 6-sum(ipart)
-                                    case nlast
-                                        slast=6-slast-snext;
-                                        nlast=0;
-                                    case nnext
-                                        if (epartp(ielem,nnext) == k)
-                                            nlast=nnext;
-                                        end
-                                    otherwise
-                                        slast=6-slast-snext;
-                                        nlast=0;
-                                end
-%  one node is in current partition
-                            else
-%  all different, so cut through centroid
-                                if (epartp(ielem,1) ~= epartp(ielem,2)) && ...
-                                   (epartp(ielem,2) ~= epartp(ielem,3)) && ...
-                                   (epartp(ielem,3) ~= epartp(ielem,1))
-                                    switch ipart
-                                        case {nlast,nnext}
-                                            if (epartp(ielem,nnext) == k)
-                                                nlast=nnext;
-                                            end
-                                        otherwise
-                                            slast=6-slast-snext;
-                                            nlast=0;
-                                    end
-%  other two are in the same partition, so cut them off
-                                else
-                                    switch ipart
-                                        case nlast
-                                            if (epartp(ielem,nnext) == k)
-                                                nlast=nnext;
-                                            end
-                                        case nnext
-                                            slast=snext;
-                                            nlast=0;
-                                        otherwise
-                                            slast=6-slast-snext;
-                                            nlast=0;
-                                    end
-                                end
-                            end
-
-%  last edge exited last element at node
-                            if nlast
-%  write out first node of first side for half-edge to midpoint
-                                disp(['segment j=' int2str(j) ' unshared half edge from node ' int2str(elemp(ielem,nlast)) ' (node ' int2str(nlast) ') on side ' int2str(slast) ' from element ' int2str(ielem) ' written.'])
-                                [lat(end+1),long(end+1)]=mapxy((md.x(elemp(ielem,nlast))),...
-                                                               (md.y(elemp(ielem,nlast))),'s');
-                                fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
-                            end
-                            nlast=0;
-                            
-%  write out midpoint of first side
-                            [lat(end+1),long(end+1)]=mapxy((md.x(elemp(ielem,slast))...
-                                                           +md.x(elemp(ielem,mod(slast,3)+1)))/2.,...
-                                                           (md.y(elemp(ielem,slast))...
-                                                           +md.y(elemp(ielem,mod(slast,3)+1)))/2.,'s');
-                            fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
-                        end
-
-%  last edge exited last element at node
-                        if nlast
-                            if elast
-%  find where last node on previous element occurs on current element
-                                nlast=find(elemp(ielem,:)==nlast,1);
-                            end
-%  half-edge occurs on unshared side from current node (unique unless mesh
-%  is only attached at node)
-                            switch nlast
-                                case 1
-                                    if ~edgeadjp(ielem,1)
-                                        nnext=2;
-                                        slast=1;
-                                    else
-                                        nnext=3;
-                                        slast=3;
-                                    end
-                                case 2
-                                    if ~edgeadjp(ielem,2)
-                                        nnext=3;
-                                        slast=2;
-                                    else
-                                        nnext=1;
-                                        slast=1;
-                                    end
-                                case 3
-                                    if ~edgeadjp(ielem,3)
-                                        nnext=1;
-                                        slast=3;
-                                    else
-                                        nnext=2;
-                                        slast=2;
-                                    end
-                            end
-%  write out half-edge from current node to midpoint of unshared side
-                            disp(['segment j=' int2str(j) ' unshared half edge from node ' int2str(elemp(ielem,nlast)) ' (node ' int2str(nlast) ') on side ' int2str(slast) ' from element ' int2str(ielem) ' written.'])
-                            [lat(end+1),long(end+1)]=mapxy((md.x(elemp(ielem,nlast))...
-                                                           +md.x(elemp(ielem,nnext)))/2.,...
-                                                           (md.y(elemp(ielem,nlast))...
-                                                           +md.y(elemp(ielem,nnext)))/2.,'s');
-                            fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
-                            nlast=0;
-
-%  last edge exited last element at midpoint of side
-                        elseif slast
-                            if elast
-%  find where last side on previous element occurs on current element
-                                slast=find(edgeadjp(ielem,:)==elast,1);
-                            end
-                        end
-
-%  find how many nodes of current element are in current partition
-%  (1 or 2, not 3) and handle each permutation separately
-                        ipart=find(epartp(ielem,:)==k);
-                        if (length(ipart) == 2)
-%  two nodes are in current partition, so cut off other node
-                            switch 6-sum(ipart)
-                                case 1
-                                    snext=3+1-slast;
-                                case 2
-                                    snext=1+2-slast;
-                                case 3
-                                    snext=2+3-slast;
-                            end
-                        else
-                            if (epartp(ielem,1) ~= epartp(ielem,2)) && ...
-                               (epartp(ielem,2) ~= epartp(ielem,3)) && ...
-                               (epartp(ielem,3) ~= epartp(ielem,1))
-%  all different, so cut through centroid
-                                disp(['element ielem=' int2str(ielem) ' centroid written.'])
-                                [lat(end+1),long(end+1)]=mapxy(sum(md.x(elemp(ielem,:)))/3.,...
-                                                               sum(md.y(elemp(ielem,:)))/3.,'s');
-                                fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
-                            end
-%  one node is in current partition, so cut off other two nodes
-                            switch ipart
-                                case 1
-                                    snext=3+1-slast;
-                                case 2
-                                    snext=1+2-slast;
-                                case 3
-                                    snext=2+3-slast;
-                            end
-                        end
-%  write out midpoint of opposite side
-                        disp(['segment j=' int2str(j) ' internal edge from side ' int2str(slast) ' to side ' int2str(snext) ' from element ' int2str(ielem) ' written.'])
-                        [lat(end+1),long(end+1)]=mapxy((md.x(elemp(ielem,snext))...
-                                                       +md.x(elemp(ielem,mod(snext,3)+1)))/2.,...
-                                                       (md.y(elemp(ielem,snext))...
-                                                       +md.y(elemp(ielem,mod(snext,3)+1)))/2.,'s');
-                        fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
-                        elast=ielem;
-                        nlast=0;
-                        slast=snext;
-%  find adjacent element to opposite side
-                        ielem=edgeadjp(elast,slast);
-%  if opposite side is unshared, find it in edge perimeter list
-                        if ~ielem
-                            jlast=find(elemper(j:end)==elast)+j-1;
-                            j=0;
-                            for l=1:length(jlast)
-                                if ((elemp(elast,slast)          == edgeper(jlast(l),1)) && ...
-                                    (elemp(elast,mod(slast,3)+1) == edgeper(jlast(l),2))) || ...
-                                   ((elemp(elast,slast)          == edgeper(jlast(l),2)) && ...
-                                    (elemp(elast,mod(slast,3)+1) == edgeper(jlast(l),1)))
-                                    j=jlast(l);
-                                    break
-                                end
-                            end
-                            if ~j
-                                j=iloop(i+1)-1;
-                            end
-%  write out half-edge from midpoint of unshared side to node
-                            if (epartp(elast,slast) == k)
-                                nnext=slast;
-                            else
-                                nnext=mod(slast,3)+1;
-                            end
-                            disp(['segment j=' int2str(j) ' unshared half edge on side ' int2str(slast) ' to node ' int2str(elemp(elast,nnext)) ' (node ' int2str(nnext) ') from element ' int2str(elast) ' written.'])
-                            [lat(end+1),long(end+1)]=mapxy(md.x(elemp(elast,nnext)),...
-                                                           md.y(elemp(elast,nnext)),'s');
-                            fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
-                            break
-%  if not unshared, advance perimeter list and watch for end
-                        else
-                            if (elast == elemper(j))
-                                if (j+1 < iloop(i+1)) && ...
-                                   ~isempty(find(elemper(j+1:end)~=elast,1))
-                                    j=j+find(elemper(j+1:end)~=elast,1);
-                                else
-                                    break
-                                end
-                            end
-                        end
-                    end
-                    j=j+1;
-                end
-            end
-%            fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(iloop(1)),lat(iloop(1)),alt);
-
-            fprintf(fid,'              </coordinates>\n');
-            fprintf(fid,'            </LinearRing>\n');
-            fprintf(fid,'          </outerBoundaryIs>\n');
-            fprintf(fid,'        </Polygon>\n');
-            fprintf(fid,'      </Placemark>\n');
-        end
-    end
-    fprintf(fid,'    </Folder>\n');
-end
-
-end
-
+    disp(['End of file ''' filek2 ''' successfully written.']);
+end
+
+end
+
Index: /issm/trunk/src/m/kml/kml_part_edges.m
===================================================================
--- /issm/trunk/src/m/kml/kml_part_edges.m	(revision 6417)
+++ /issm/trunk/src/m/kml/kml_part_edges.m	(revision 6417)
@@ -0,0 +1,477 @@
+%%
+%  create kml linestrings for the partition edges.
+%
+%  []=kml_part_edges(fid,md,params)
+%
+%  where the required input is:
+%    fid           (char, file ID of .kml file)
+%    md            (model, model class object)
+%
+%  the optional input is:
+%    params        (string/numeric, parameter names and values)
+%
+%  and the optional input is:
+%    data          (numeric, element or nodal results data)
+%    alt           (numeric, altitude for polygons, default 10000)
+%    cmin          (numeric, minimum of color map)
+%    cmax          (numeric, maximum of color map)
+%    cmap          (char or numeric, colormap definition)
+%    prtplt        (char, 'off'/'no' for partition segment plot)
+%
+function []=kml_part_edges(varargin)
+
+if ~nargin
+    help kml_part_edges
+    return
+end
+
+%%  process input data
+
+iarg=1;
+if (nargin >= 1)
+    fid=varargin{1};
+end
+if ~exist('fid','var') || isempty(fid) || (fid < 0)
+    error('File ID ''%d'' must be open.',fid);
+end
+
+iarg=iarg+1;
+if (nargin >= 2)
+    md=varargin{2};
+end
+if ~exist('md','var') || isempty(md) || ~isa(md,'model')
+    error(['Model ''' inputname(iarg) ''' is unrecognized class ''' class(md) '''.']);
+end
+
+%  parameters
+
+iarg=iarg+1;
+while (iarg <= nargin-1)
+    if ischar(varargin{iarg})
+        eval([varargin{iarg} '=varargin{iarg+1};']);
+        if (numel(varargin{iarg+1}) <= 20)
+            disp([varargin{iarg} '=' any2str(varargin{iarg+1}) ';']);
+        else
+            disp([varargin{iarg} '=' string_size(varargin{iarg+1}) ' ' class(varargin{iarg+1}) ';']);
+        end
+        if strcmpi(varargin{iarg},'data')
+            cdata=inputname(iarg+1);
+        end
+    else
+        error(['''' any2str(varargin{iarg}) ''' is not a parameter name.']);
+    end
+    iarg=iarg+2;
+end
+
+if exist('data','var') && ~isempty(data)
+    if     (numel(data)==md.numberofelements)
+        edata=data;
+    elseif (numel(data)==md.numberofgrids)
+        ndata=data;
+        display('Averaging nodal data to element data.');
+        edata=zeros(1,md.numberofelements);
+        for i=1:size(md.elements,1)
+            for j=1:size(md.elements,2)
+                edata(i)=edata(i)+ndata(md.elements(i,j));
+            end
+            edata(i)=edata(i)/size(md.elements,2);
+        end
+    else
+        error(['Data has incorrect number of ' num2str(numel(data)) ' values.']);
+    end
+end
+
+%  colormap command operates on a figure, so create an invisible one
+%  (could also directly call colormaps, e.g. jet(64), but risky)
+
+hfig=figure('Visible','off');
+if exist('cmap','var')
+    colormap(cmap)
+end
+cmap=colormap;
+close(hfig)
+    
+if exist('edata','var')
+    if ~exist('cmin','var')
+        cmin=min(min(edata));
+    end
+    if ~exist('cmax','var')
+        cmax=max(max(edata));
+    end
+end
+
+if ~exist('alt','var')
+    alt=10000;
+end
+
+%%  write folder for partition edges
+
+if (~exist('prtplt','var') || strncmpi(prtplt,'on' ,2) || strncmpi(prtplt,'y',1)) && ...
+    md.npart
+    fprintf(fid,'    <Folder>\n');
+    fprintf(fid,'      <name>Partition Edges</name>\n');
+    fprintf(fid,'      <visibility>1</visibility>\n');
+    fprintf(fid,'      <description>Partitions=%d, Nodes=%d</description>\n',...
+        md.npart,md.numberofgrids);
+
+%  write each partition loop as linestrings
+
+    disp(['Writing ' num2str(md.npart) ' partitions as KML linestrings.']);
+    epart=md.part(md.elements)+1;
+    if exist('ndata','var') || exist('edata','var')
+        pdata=zeros(1,md.npart);
+        pdata(:)=NaN;
+    end
+
+%  loop over each partition
+
+    for k=1:md.npart
+%        disp(['partition k=' int2str(k)])
+        
+%  for each partition, find all the included elements and determine the
+%  perimeter (including those shared by another partition)
+
+        [icol,irow]=find(epart'==k);
+        irow=unique(irow);
+        elemp=md.elements(irow,:);
+        epartp=epart(irow,:);
+        nodeconp=nodeconnectivity(elemp,md.numberofgrids);
+        [edgeadjp]=edgeadjacency(elemp,nodeconp);
+        [edgeper,elemper,iloop]=edgeperimeter(elemp,nodeconp,edgeadjp);
+        iloop(end+1)=size(edgeper,1)+1;
+        
+%  determine the data to be used for the colors (if any)
+
+        if exist('ndata','var')
+            pdata(k)=ndata(find(md.part+1==k,1));
+        elseif exist('edata','var')
+            for i=1:size(epartp,1)
+                if isempty(find(epart(i,:)~=k,1))
+                    pdata(k)=edata(irow(i));
+                    break
+                end
+            end
+            if isnan(pdata(k))
+                warning('Data for Partition %d is not defined.\n',k)
+            end
+        end
+        
+%  loop over each loop of the perimeter for the given partition
+
+        for i=1:length(iloop)-1
+            fprintf(fid,'      <Placemark>\n');
+            if (length(iloop)-1 > 1)
+                fprintf(fid,'        <name>Partition %d, Loop %d</name>\n',k,i);
+            else
+                fprintf(fid,'        <name>Partition %d</name>\n',k);
+            end
+            fprintf(fid,'        <visibility>1</visibility>\n');
+            if exist('pdata','var')
+                fprintf(fid,'        <description>Partition data: %g</description>\n',pdata(k));
+                imap = fix((pdata(k)-cmin)/(cmax-cmin)*size(cmap,1))+1;
+                if     (imap >= 1) && (imap <= size(cmap,1))
+                    fprintf(fid,'        <styleUrl>#MatlabColor%d</styleUrl>\n',imap);
+                elseif (pdata(k) == cmax)
+                    fprintf(fid,'        <styleUrl>#MatlabColor%d</styleUrl>\n',size(cmap,1));
+                else
+                    fprintf(fid,'        <styleUrl>#BlackLineEmptyPoly</styleUrl>\n');
+                end
+            else
+                fprintf(fid,'        <styleUrl>#BlackLineRandomPoly</styleUrl>\n');
+            end
+%            fprintf(fid,'        <Polygon>\n');
+%            fprintf(fid,'          <extrude>1</extrude>\n');
+%            fprintf(fid,'          <altitudeMode>relativeToGround</altitudeMode>\n');
+%            fprintf(fid,'          <outerBoundaryIs>\n');
+%            fprintf(fid,'            <LinearRing>\n');
+        fprintf(fid,'        <LineString>\n');
+        fprintf(fid,'          <extrude>1</extrude>\n');
+        fprintf(fid,'          <tessellate>1</tessellate>\n');
+        fprintf(fid,'          <altitudeMode>relativeToGround</altitudeMode>\n');
+            fprintf(fid,'              <coordinates>\n');
+            elast=0;
+            nlast=0;
+            slast=0;
+            lat=[];
+            long=[];
+            
+%  loop over the element edges on the loop of the partition
+
+            j=iloop(i);
+            while (j < iloop(i+1))
+%  find which side of element is referenced in perimeter list
+                for l=1:size(elemp,2)
+                    if ((elemp(elemper(j),l)          == edgeper(j,1)) && ...
+                        (elemp(elemper(j),mod(l,3)+1) == edgeper(j,2))) || ...
+                       ((elemp(elemper(j),l)          == edgeper(j,2)) && ...
+                        (elemp(elemper(j),mod(l,3)+1) == edgeper(j,1)))
+                        jedge=l;
+                        break
+                    end
+                end
+
+%  check if element side connects nodes in partition
+                if (epartp(elemper(j),jedge)          == k) && ...
+                   (epartp(elemper(j),mod(jedge,3)+1) == k)
+%  write out specified element side
+%                    disp(['segment j=' int2str(j) ' unshared edge ' int2str(edgeper(j,1)) ' to ' int2str(edgeper(j,2)) ' on side ' int2str(jedge) ' from element ' int2str(elemper(j)) ' written.'])
+%  if first edge, write out first node
+                    if ~elast
+                        [lat(end+1),long(end+1)]=mapxy(md.x(edgeper(j,1)),md.y(edgeper(j,1)),'s');
+                        fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
+                    end
+                    [lat(end+1),long(end+1)]=mapxy(md.x(edgeper(j,2)),md.y(edgeper(j,2)),'s');
+                    fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
+                    elast=elemper(j);
+                    nlast=edgeper(j,2);
+                    slast=0;
+                    j=j+1;
+                    
+%  element not entirely within partition, so figure out boundary
+                else
+%                    disp(['segment j=' int2str(j) ' from element ' int2str(elemper(j)) ' shared by other partitions.'])
+                    ielem=elemper(j);
+
+%  follow partition boundary through elements not wholly in partition
+%  (may include elements not in perimeter list)
+
+                    while 1
+%  if first edge, figure out direction from perimeter edge direction
+                        if ~nlast && ~slast
+                            nlast=find(elemp(ielem,:)==edgeper(j,1));
+                            nnext=find(elemp(ielem,:)==edgeper(j,2));
+                            if     (nlast+nnext == 3)
+                                slast=1;
+                            elseif (nlast+nnext == 5)
+                                slast=2;
+                            elseif (nlast+nnext == 4)
+                                slast=3;
+                            end
+                            if     (nnext+(6-nlast-nnext) == 3)
+                                snext=1;
+                            elseif (nnext+(6-nlast-nnext) == 5)
+                                snext=2;
+                            elseif (nnext+(6-nlast-nnext) == 4)
+                                snext=3;
+                            end
+
+%  find how many nodes of current element are in current partition
+%  (1 or 2, not 3) and handle each permutation separately
+                            ipart=find(epartp(ielem,:)==k);
+%  two nodes are in current partition, so cut off other node
+                            if (length(ipart) == 2)
+                                switch 6-sum(ipart)
+                                    case nlast
+                                        slast=6-slast-snext;
+                                        nlast=0;
+                                    case nnext
+                                        if (epartp(ielem,nnext) == k)
+                                            nlast=nnext;
+                                        end
+                                    otherwise
+                                        slast=6-slast-snext;
+                                        nlast=0;
+                                end
+%  one node is in current partition
+                            else
+%  all different, so cut through centroid
+                                if (epartp(ielem,1) ~= epartp(ielem,2)) && ...
+                                   (epartp(ielem,2) ~= epartp(ielem,3)) && ...
+                                   (epartp(ielem,3) ~= epartp(ielem,1))
+                                    switch ipart
+                                        case {nlast,nnext}
+                                            if (epartp(ielem,nnext) == k)
+                                                nlast=nnext;
+                                            end
+                                        otherwise
+                                            slast=6-slast-snext;
+                                            nlast=0;
+                                    end
+%  other two are in the same partition, so cut them off
+                                else
+                                    switch ipart
+                                        case nlast
+                                            if (epartp(ielem,nnext) == k)
+                                                nlast=nnext;
+                                            end
+                                        case nnext
+                                            slast=snext;
+                                            nlast=0;
+                                        otherwise
+                                            slast=6-slast-snext;
+                                            nlast=0;
+                                    end
+                                end
+                            end
+
+%  last edge exited last element at node
+                            if nlast
+%  write out first node of first side for half-edge to midpoint
+%                                disp(['segment j=' int2str(j) ' unshared half edge from node ' int2str(elemp(ielem,nlast)) ' (node ' int2str(nlast) ') on side ' int2str(slast) ' from element ' int2str(ielem) ' written.'])
+                                [lat(end+1),long(end+1)]=mapxy((md.x(elemp(ielem,nlast))),...
+                                                               (md.y(elemp(ielem,nlast))),'s');
+                                fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
+                            end
+                            nlast=0;
+                            
+%  write out midpoint of first side
+                            [lat(end+1),long(end+1)]=mapxy((md.x(elemp(ielem,slast))...
+                                                           +md.x(elemp(ielem,mod(slast,3)+1)))/2.,...
+                                                           (md.y(elemp(ielem,slast))...
+                                                           +md.y(elemp(ielem,mod(slast,3)+1)))/2.,'s');
+                            fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
+                        end
+
+%  last edge exited last element at node
+                        if nlast
+                            if elast
+%  find where last node on previous element occurs on current element
+                                nlast=find(elemp(ielem,:)==nlast,1);
+                            end
+%  half-edge occurs on unshared side from current node (unique unless mesh
+%  is only attached at node)
+                            switch nlast
+                                case 1
+                                    if ~edgeadjp(ielem,1)
+                                        nnext=2;
+                                        slast=1;
+                                    else
+                                        nnext=3;
+                                        slast=3;
+                                    end
+                                case 2
+                                    if ~edgeadjp(ielem,2)
+                                        nnext=3;
+                                        slast=2;
+                                    else
+                                        nnext=1;
+                                        slast=1;
+                                    end
+                                case 3
+                                    if ~edgeadjp(ielem,3)
+                                        nnext=1;
+                                        slast=3;
+                                    else
+                                        nnext=2;
+                                        slast=2;
+                                    end
+                            end
+%  write out half-edge from current node to midpoint of unshared side
+%                            disp(['segment j=' int2str(j) ' unshared half edge from node ' int2str(elemp(ielem,nlast)) ' (node ' int2str(nlast) ') on side ' int2str(slast) ' from element ' int2str(ielem) ' written.'])
+                            [lat(end+1),long(end+1)]=mapxy((md.x(elemp(ielem,nlast))...
+                                                           +md.x(elemp(ielem,nnext)))/2.,...
+                                                           (md.y(elemp(ielem,nlast))...
+                                                           +md.y(elemp(ielem,nnext)))/2.,'s');
+                            fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
+                            nlast=0;
+
+%  last edge exited last element at midpoint of side
+                        elseif slast
+                            if elast
+%  find where last side on previous element occurs on current element
+                                slast=find(edgeadjp(ielem,:)==elast,1);
+                            end
+                        end
+
+%  find how many nodes of current element are in current partition
+%  (1 or 2, not 3) and handle each permutation separately
+                        ipart=find(epartp(ielem,:)==k);
+                        if (length(ipart) == 2)
+%  two nodes are in current partition, so cut off other node
+                            switch 6-sum(ipart)
+                                case 1
+                                    snext=3+1-slast;
+                                case 2
+                                    snext=1+2-slast;
+                                case 3
+                                    snext=2+3-slast;
+                            end
+                        else
+                            if (epartp(ielem,1) ~= epartp(ielem,2)) && ...
+                               (epartp(ielem,2) ~= epartp(ielem,3)) && ...
+                               (epartp(ielem,3) ~= epartp(ielem,1))
+%  all different, so cut through centroid
+%                                disp(['element ielem=' int2str(ielem) ' centroid written.'])
+                                [lat(end+1),long(end+1)]=mapxy(sum(md.x(elemp(ielem,:)))/3.,...
+                                                               sum(md.y(elemp(ielem,:)))/3.,'s');
+                                fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
+                            end
+%  one node is in current partition, so cut off other two nodes
+                            switch ipart
+                                case 1
+                                    snext=3+1-slast;
+                                case 2
+                                    snext=1+2-slast;
+                                case 3
+                                    snext=2+3-slast;
+                            end
+                        end
+%  write out midpoint of opposite side
+%                        disp(['segment j=' int2str(j) ' internal edge from side ' int2str(slast) ' to side ' int2str(snext) ' from element ' int2str(ielem) ' written.'])
+                        [lat(end+1),long(end+1)]=mapxy((md.x(elemp(ielem,snext))...
+                                                       +md.x(elemp(ielem,mod(snext,3)+1)))/2.,...
+                                                       (md.y(elemp(ielem,snext))...
+                                                       +md.y(elemp(ielem,mod(snext,3)+1)))/2.,'s');
+                        fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
+                        elast=ielem;
+                        nlast=0;
+                        slast=snext;
+%  find adjacent element to opposite side
+                        ielem=edgeadjp(elast,slast);
+%  if opposite side is unshared, find it in edge perimeter list
+                        if ~ielem
+                            jlast=find(elemper(j:end)==elast)+j-1;
+                            j=0;
+                            for l=1:length(jlast)
+                                if ((elemp(elast,slast)          == edgeper(jlast(l),1)) && ...
+                                    (elemp(elast,mod(slast,3)+1) == edgeper(jlast(l),2))) || ...
+                                   ((elemp(elast,slast)          == edgeper(jlast(l),2)) && ...
+                                    (elemp(elast,mod(slast,3)+1) == edgeper(jlast(l),1)))
+                                    j=jlast(l);
+                                    break
+                                end
+                            end
+                            if ~j
+                                j=iloop(i+1)-1;
+                            end
+%  write out half-edge from midpoint of unshared side to node
+                            if (epartp(elast,slast) == k)
+                                nnext=slast;
+                            else
+                                nnext=mod(slast,3)+1;
+                            end
+%                            disp(['segment j=' int2str(j) ' unshared half edge on side ' int2str(slast) ' to node ' int2str(elemp(elast,nnext)) ' (node ' int2str(nnext) ') from element ' int2str(elast) ' written.'])
+                            [lat(end+1),long(end+1)]=mapxy(md.x(elemp(elast,nnext)),...
+                                                           md.y(elemp(elast,nnext)),'s');
+                            fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
+                            break
+%  if not unshared, advance perimeter list and watch for end
+                        else
+                            if (elast == elemper(j))
+                                if (j+1 < iloop(i+1)) && ...
+                                   ~isempty(find(elemper(j+1:end)~=elast,1))
+                                    j=j+find(elemper(j+1:end)~=elast,1);
+                                else
+                                    break
+                                end
+                            end
+                        end
+                    end
+                    j=j+1;
+                end
+            end
+%            fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(iloop(i)),lat(iloop(i)),alt);
+
+            fprintf(fid,'              </coordinates>\n');
+        fprintf(fid,'        </LineString>\n');
+%            fprintf(fid,'            </LinearRing>\n');
+%            fprintf(fid,'          </outerBoundaryIs>\n');
+%            fprintf(fid,'        </Polygon>\n');
+            fprintf(fid,'      </Placemark>\n');
+        end
+    end
+    fprintf(fid,'    </Folder>\n');
+end
+
+end
+
Index: /issm/trunk/src/m/kml/kml_part_elems.m
===================================================================
--- /issm/trunk/src/m/kml/kml_part_elems.m	(revision 6417)
+++ /issm/trunk/src/m/kml/kml_part_elems.m	(revision 6417)
@@ -0,0 +1,189 @@
+%%
+%  create kml polygons for the partition elements.
+%
+%  []=kml_part_elems(fid,md,params)
+%
+%  where the required input is:
+%    fid           (char, file ID of .kml file)
+%    md            (model, model class object)
+%
+%  the optional input is:
+%    params        (string/numeric, parameter names and values)
+%
+%  and the optional input is:
+%    data          (numeric, element or nodal results data)
+%    alt           (numeric, altitude for polygons, default 10000)
+%    cmin          (numeric, minimum of color map)
+%    cmax          (numeric, maximum of color map)
+%    cmap          (char or numeric, colormap definition)
+%    prtplt        (char, 'off'/'no' for partition segment plot)
+%
+function []=kml_part_elems(varargin)
+
+if ~nargin
+    help kml_part_elems
+    return
+end
+
+%%  process input data
+
+iarg=1;
+if (nargin >= 1)
+    fid=varargin{1};
+end
+if ~exist('fid','var') || isempty(fid) || (fid < 0)
+    error('File ID ''%d'' must be open.',fid);
+end
+
+iarg=iarg+1;
+if (nargin >= 2)
+    md=varargin{2};
+end
+if ~exist('md','var') || isempty(md) || ~isa(md,'model')
+    error(['Model ''' inputname(iarg) ''' is unrecognized class ''' class(md) '''.']);
+end
+
+%  parameters
+
+iarg=iarg+1;
+while (iarg <= nargin-1)
+    if ischar(varargin{iarg})
+        eval([varargin{iarg} '=varargin{iarg+1};']);
+        if (numel(varargin{iarg+1}) <= 20)
+            disp([varargin{iarg} '=' any2str(varargin{iarg+1}) ';']);
+        else
+            disp([varargin{iarg} '=' string_size(varargin{iarg+1}) ' ' class(varargin{iarg+1}) ';']);
+        end
+        if strcmpi(varargin{iarg},'data')
+            cdata=inputname(iarg+1);
+        end
+    else
+        error(['''' any2str(varargin{iarg}) ''' is not a parameter name.']);
+    end
+    iarg=iarg+2;
+end
+
+if exist('data','var') && ~isempty(data)
+    if     (numel(data)==md.numberofelements)
+        edata=data;
+    elseif (numel(data)==md.numberofgrids)
+        ndata=data;
+        display('Averaging nodal data to element data.');
+        edata=zeros(1,md.numberofelements);
+        for i=1:size(md.elements,1)
+            for j=1:size(md.elements,2)
+                edata(i)=edata(i)+ndata(md.elements(i,j));
+            end
+            edata(i)=edata(i)/size(md.elements,2);
+        end
+    else
+        error(['Data has incorrect number of ' num2str(numel(data)) ' values.']);
+    end
+end
+
+%  colormap command operates on a figure, so create an invisible one
+%  (could also directly call colormaps, e.g. jet(64), but risky)
+
+hfig=figure('Visible','off');
+if exist('cmap','var')
+    colormap(cmap)
+end
+cmap=colormap;
+close(hfig)
+    
+if exist('edata','var')
+    if ~exist('cmin','var')
+        cmin=min(min(edata));
+    end
+    if ~exist('cmax','var')
+        cmax=max(max(edata));
+    end
+end
+
+if ~exist('alt','var')
+    alt=10000;
+end
+
+%  write folder for partition elements
+
+if (~exist('prtplt','var') || strncmpi(prtplt,'on' ,2) || strncmpi(prtplt,'y',1)) && ...
+    md.npart
+    fprintf(fid,'    <Folder>\n');
+    fprintf(fid,'      <name>Partition Elements</name>\n');
+    fprintf(fid,'      <visibility>1</visibility>\n');
+    fprintf(fid,'      <description>Partitions=%d, Nodes=%d</description>\n',...
+        md.npart,md.numberofgrids);
+
+%  write each partition loop as a polygon
+
+    disp(['Writing ' num2str(md.npart) ' partitions as KML polygons.']);
+    epart=md.part(md.elements)+1;
+    if exist('ndata','var') || exist('edata','var')
+        pdata=zeros(1,md.npart);
+        pdata(:)=NaN;
+    end
+    for k=1:md.npart
+        [icol,irow]=find(epart'==k);
+        irow=unique(irow);
+        elem=md.elements(irow,:);
+        nodecon=nodeconnectivity(elem,md.numberofgrids);
+        [edgeper,elemper,iloop]=edgeperimeter(elem,nodecon);
+        iloop(end+1)=size(edgeper,1)+1;
+        if exist('ndata','var')
+            pdata(k)=ndata(find(md.part+1==k,1));
+        elseif exist('edata','var')
+            for i=1:size(epart,1)
+                if isempty(find(epart(i,:)~=k,1))
+                    pdata(k)=edata(i);
+                    break
+                end
+            end
+            if isnan(pdata(k))
+                warning('Data for Partition %d is not defined.\n',k)
+            end
+        end
+        for i=1:length(iloop)-1
+            fprintf(fid,'      <Placemark>\n');
+            if (length(iloop)-1 > 1)
+                fprintf(fid,'        <name>Partition %d, Loop %d</name>\n',k,i);
+            else
+                fprintf(fid,'        <name>Partition %d</name>\n',k);
+            end
+            fprintf(fid,'        <visibility>1</visibility>\n');
+            if exist('pdata','var')
+                fprintf(fid,'        <description>Partition data: %g</description>\n',pdata(k));
+                imap = fix((pdata(k)-cmin)/(cmax-cmin)*size(cmap,1))+1;
+                if     (imap >= 1) && (imap <= size(cmap,1))
+                    fprintf(fid,'        <styleUrl>#MatlabColor%d</styleUrl>\n',imap);
+                elseif (pdata(k) == cmax)
+                    fprintf(fid,'        <styleUrl>#MatlabColor%d</styleUrl>\n',size(cmap,1));
+                else
+                    fprintf(fid,'        <styleUrl>#BlackLineEmptyPoly</styleUrl>\n');
+                end
+            else
+                fprintf(fid,'        <styleUrl>#BlackLineRandomPoly</styleUrl>\n');
+            end
+            fprintf(fid,'        <Polygon>\n');
+            fprintf(fid,'          <extrude>1</extrude>\n');
+            fprintf(fid,'          <altitudeMode>relativeToGround</altitudeMode>\n');
+            fprintf(fid,'          <outerBoundaryIs>\n');
+            fprintf(fid,'            <LinearRing>\n');
+            fprintf(fid,'              <coordinates>\n');
+            for j=iloop(i):iloop(i+1)-1
+                [lat(j),long(j)]=mapxy(md.x(edgeper(j,1)),md.y(edgeper(j,1)),'s');
+                fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(j),lat(j),alt);
+            end
+            fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(iloop(i)),lat(iloop(i)),alt);
+
+            fprintf(fid,'              </coordinates>\n');
+            fprintf(fid,'            </LinearRing>\n');
+            fprintf(fid,'          </outerBoundaryIs>\n');
+            fprintf(fid,'        </Polygon>\n');
+            fprintf(fid,'      </Placemark>\n');
+        end
+    end
+    fprintf(fid,'    </Folder>\n');
+end
+
+end
+
Index: /issm/trunk/src/m/kml/kml_part_flagedges.m
===================================================================
--- /issm/trunk/src/m/kml/kml_part_flagedges.m	(revision 6417)
+++ /issm/trunk/src/m/kml/kml_part_flagedges.m	(revision 6417)
@@ -0,0 +1,100 @@
+%%
+%  create kml linestrings for the flagged partition edges.
+%
+%  []=kml_part_flagedges(fid,md,params)
+%
+%  where the required input is:
+%    fid           (char, file ID of .kml file)
+%    md            (model, model class object)
+%
+%  the optional input is:
+%    params        (string/numeric, parameter names and values)
+%
+%  and the optional input is:
+%    alt           (numeric, altitude for polygons, default 10000)
+%    prtplt        (char, 'off'/'no' for partition segment plot)
+%
+function []=kml_part_flagedges(varargin)
+
+if ~nargin
+    help kml_part_flagedges
+    return
+end
+
+%%  process input data
+
+iarg=1;
+if (nargin >= 1)
+    fid=varargin{1};
+end
+if ~exist('fid','var') || isempty(fid) || (fid < 0)
+    error('File ID ''%d'' must be open.',fid);
+end
+
+iarg=iarg+1;
+if (nargin >= 2)
+    md=varargin{2};
+end
+if ~exist('md','var') || isempty(md) || ~isa(md,'model')
+    error(['Model ''' inputname(iarg) ''' is unrecognized class ''' class(md) '''.']);
+end
+
+%  parameters
+
+iarg=iarg+1;
+while (iarg <= nargin-1)
+    if ischar(varargin{iarg})
+        eval([varargin{iarg} '=varargin{iarg+1};']);
+        if (numel(varargin{iarg+1}) <= 20)
+            disp([varargin{iarg} '=' any2str(varargin{iarg+1}) ';']);
+        else
+            disp([varargin{iarg} '=' string_size(varargin{iarg+1}) ' ' class(varargin{iarg+1}) ';']);
+        end
+    else
+        error(['''' any2str(varargin{iarg}) ''' is not a parameter name.']);
+    end
+    iarg=iarg+2;
+end
+
+if ~exist('alt','var')
+    alt=10000;
+end
+
+%%  write folder for partition segments
+
+if (~exist('prtplt','var') || strncmpi(prtplt,'on' ,2) || strncmpi(prtplt,'y',1)) && ...
+    md.npart
+    [xseg,yseg]=flagedges(md.elements,md.x,md.y,md.part);
+    fprintf(fid,'    <Folder>\n');
+    fprintf(fid,'      <name>Partition Segments</name>\n');
+    fprintf(fid,'      <visibility>1</visibility>\n');
+    fprintf(fid,'      <description>Partitions=%d, Segments=%d</description>\n',...
+        md.npart,size(xseg,1));
+
+%  write each segment as a linestring
+
+    disp(['Writing ' num2str(size(xseg,1)) ' partition segments as KML linestrings.']);
+    for i=1:size(xseg,1)
+        fprintf(fid,'      <Placemark>\n');
+        fprintf(fid,'        <name>Segment %d</name>\n',i);
+        fprintf(fid,'        <visibility>1</visibility>\n');
+        fprintf(fid,'        <styleUrl>#RedLineRedPoly</styleUrl>\n');
+        fprintf(fid,'        <LineString>\n');
+        fprintf(fid,'          <extrude>1</extrude>\n');
+        fprintf(fid,'          <tessellate>1</tessellate>\n');
+        fprintf(fid,'          <altitudeMode>relativeToGround</altitudeMode>\n');
+        fprintf(fid,'          <coordinates>\n');
+        for j=1:2
+            [lat(j),long(j)]=mapxy(xseg(i,j),yseg(i,j),'s');
+            fprintf(fid,'            %0.16g,%0.16g,%0.16g\n',long(j),lat(j),alt);
+        end
+
+        fprintf(fid,'          </coordinates>\n');
+        fprintf(fid,'        </LineString>\n');
+        fprintf(fid,'      </Placemark>\n');
+    end
+    fprintf(fid,'    </Folder>\n');
+end
+
+end
+
Index: /issm/trunk/src/m/kml/kml_partitions.m
===================================================================
--- /issm/trunk/src/m/kml/kml_partitions.m	(revision 6417)
+++ /issm/trunk/src/m/kml/kml_partitions.m	(revision 6417)
@@ -0,0 +1,472 @@
+%%
+%  create kml polygons for the partitions.
+%
+%  []=kml_partitions(fid,md,params)
+%
+%  where the required input is:
+%    fid           (char, file ID of .kml file)
+%    md            (model, model class object)
+%
+%  the optional input is:
+%    params        (string/numeric, parameter names and values)
+%
+%  and the optional input is:
+%    data          (numeric, element or nodal results data)
+%    alt           (numeric, altitude for polygons, default 10000)
+%    cmin          (numeric, minimum of color map)
+%    cmax          (numeric, maximum of color map)
+%    cmap          (char or numeric, colormap definition)
+%    prtplt        (char, 'off'/'no' for partition segment plot)
+%
+function []=kml_partitions(varargin)
+
+if ~nargin
+    help kml_part_edges
+    return
+end
+
+%%  process input data
+
+iarg=1;
+if (nargin >= 1)
+    fid=varargin{1};
+end
+if ~exist('fid','var') || isempty(fid) || (fid < 0)
+    error('File ID ''%d'' must be open.',fid);
+end
+
+iarg=iarg+1;
+if (nargin >= 2)
+    md=varargin{2};
+end
+if ~exist('md','var') || isempty(md) || ~isa(md,'model')
+    error(['Model ''' inputname(iarg) ''' is unrecognized class ''' class(md) '''.']);
+end
+
+%  parameters
+
+iarg=iarg+1;
+while (iarg <= nargin-1)
+    if ischar(varargin{iarg})
+        eval([varargin{iarg} '=varargin{iarg+1};']);
+        if (numel(varargin{iarg+1}) <= 20)
+            disp([varargin{iarg} '=' any2str(varargin{iarg+1}) ';']);
+        else
+            disp([varargin{iarg} '=' string_size(varargin{iarg+1}) ' ' class(varargin{iarg+1}) ';']);
+        end
+        if strcmpi(varargin{iarg},'data')
+            cdata=inputname(iarg+1);
+        end
+    else
+        error(['''' any2str(varargin{iarg}) ''' is not a parameter name.']);
+    end
+    iarg=iarg+2;
+end
+
+if exist('data','var') && ~isempty(data)
+    if     (numel(data)==md.numberofelements)
+        edata=data;
+    elseif (numel(data)==md.numberofgrids)
+        ndata=data;
+        display('Averaging nodal data to element data.');
+        edata=zeros(1,md.numberofelements);
+        for i=1:size(md.elements,1)
+            for j=1:size(md.elements,2)
+                edata(i)=edata(i)+ndata(md.elements(i,j));
+            end
+            edata(i)=edata(i)/size(md.elements,2);
+        end
+    else
+        error(['Data has incorrect number of ' num2str(numel(data)) ' values.']);
+    end
+end
+
+%  colormap command operates on a figure, so create an invisible one
+%  (could also directly call colormaps, e.g. jet(64), but risky)
+
+hfig=figure('Visible','off');
+if exist('cmap','var')
+    colormap(cmap)
+end
+cmap=colormap;
+close(hfig)
+    
+if exist('edata','var')
+    if ~exist('cmin','var')
+        cmin=min(min(edata));
+    end
+    if ~exist('cmax','var')
+        cmax=max(max(edata));
+    end
+end
+
+if ~exist('alt','var')
+    alt=10000;
+end
+
+%%  write folder for partitions
+
+if (~exist('prtplt','var') || strncmpi(prtplt,'on' ,2) || strncmpi(prtplt,'y',1)) && ...
+    md.npart
+    fprintf(fid,'    <Folder>\n');
+    fprintf(fid,'      <name>Partitions</name>\n');
+    fprintf(fid,'      <visibility>1</visibility>\n');
+    fprintf(fid,'      <description>Partitions=%d, Nodes=%d</description>\n',...
+        md.npart,md.numberofgrids);
+
+%  write each partition loop as linestrings
+
+    disp(['Writing ' num2str(md.npart) ' partitions as KML polygons.']);
+    epart=md.part(md.elements)+1;
+    if exist('ndata','var') || exist('edata','var')
+        pdata=zeros(1,md.npart);
+        pdata(:)=NaN;
+    end
+
+%  loop over each partition
+
+    for k=1:md.npart
+%        disp(['partition k=' int2str(k)])
+        
+%  for each partition, find all the included elements and determine the
+%  perimeter (including those shared by another partition)
+
+        [icol,irow]=find(epart'==k);
+        irow=unique(irow);
+        elemp=md.elements(irow,:);
+        epartp=epart(irow,:);
+        nodeconp=nodeconnectivity(elemp,md.numberofgrids);
+        [edgeadjp]=edgeadjacency(elemp,nodeconp);
+        [edgeper,elemper,iloop]=edgeperimeter(elemp,nodeconp,edgeadjp);
+        iloop(end+1)=size(edgeper,1)+1;
+        
+%  determine the data to be used for the colors (if any)
+
+        if exist('ndata','var')
+            pdata(k)=ndata(find(md.part+1==k,1));
+        elseif exist('edata','var')
+            for i=1:size(epartp,1)
+                if isempty(find(epart(i,:)~=k,1))
+                    pdata(k)=edata(irow(i));
+                    break
+                end
+            end
+            if isnan(pdata(k))
+                warning('Data for Partition %d is not defined.\n',k)
+            end
+        end
+        
+%  loop over each loop of the perimeter for the given partition
+
+        for i=1:length(iloop)-1
+            fprintf(fid,'      <Placemark>\n');
+            if (length(iloop)-1 > 1)
+                fprintf(fid,'        <name>Partition %d, Loop %d</name>\n',k,i);
+            else
+                fprintf(fid,'        <name>Partition %d</name>\n',k);
+            end
+            fprintf(fid,'        <visibility>1</visibility>\n');
+            if exist('pdata','var')
+                fprintf(fid,'        <description>Partition data: %g</description>\n',pdata(k));
+                imap = fix((pdata(k)-cmin)/(cmax-cmin)*size(cmap,1))+1;
+                if     (imap >= 1) && (imap <= size(cmap,1))
+                    fprintf(fid,'        <styleUrl>#MatlabColor%d</styleUrl>\n',imap);
+                elseif (pdata(k) == cmax)
+                    fprintf(fid,'        <styleUrl>#MatlabColor%d</styleUrl>\n',size(cmap,1));
+                else
+                    fprintf(fid,'        <styleUrl>#BlackLineEmptyPoly</styleUrl>\n');
+                end
+            else
+                fprintf(fid,'        <styleUrl>#BlackLineRandomPoly</styleUrl>\n');
+            end
+            fprintf(fid,'        <Polygon>\n');
+            fprintf(fid,'          <extrude>1</extrude>\n');
+            fprintf(fid,'          <altitudeMode>relativeToGround</altitudeMode>\n');
+            fprintf(fid,'          <outerBoundaryIs>\n');
+            fprintf(fid,'            <LinearRing>\n');
+            fprintf(fid,'              <coordinates>\n');
+            elast=0;
+            nlast=0;
+            slast=0;
+            lat=[];
+            long=[];
+            
+%  loop over the element edges on the loop of the partition
+
+            j=iloop(i);
+            while (j < iloop(i+1))
+%  find which side of element is referenced in perimeter list
+                for l=1:size(elemp,2)
+                    if ((elemp(elemper(j),l)          == edgeper(j,1)) && ...
+                        (elemp(elemper(j),mod(l,3)+1) == edgeper(j,2))) || ...
+                       ((elemp(elemper(j),l)          == edgeper(j,2)) && ...
+                        (elemp(elemper(j),mod(l,3)+1) == edgeper(j,1)))
+                        jedge=l;
+                        break
+                    end
+                end
+
+%  check if element side connects nodes in partition
+                if (epartp(elemper(j),jedge)          == k) && ...
+                   (epartp(elemper(j),mod(jedge,3)+1) == k)
+%  write out specified element side
+%                    disp(['segment j=' int2str(j) ' unshared edge ' int2str(edgeper(j,1)) ' to ' int2str(edgeper(j,2)) ' on side ' int2str(jedge) ' from element ' int2str(elemper(j)) ' written.'])
+%  if first edge, write out first node
+                    if ~elast
+                        [lat(end+1),long(end+1)]=mapxy(md.x(edgeper(j,1)),md.y(edgeper(j,1)),'s');
+                        fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
+                    end
+                    [lat(end+1),long(end+1)]=mapxy(md.x(edgeper(j,2)),md.y(edgeper(j,2)),'s');
+                    fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
+                    elast=elemper(j);
+                    nlast=edgeper(j,2);
+                    slast=0;
+                    j=j+1;
+                    
+%  element not entirely within partition, so figure out boundary
+                else
+%                    disp(['segment j=' int2str(j) ' from element ' int2str(elemper(j)) ' shared by other partitions.'])
+                    ielem=elemper(j);
+
+%  follow partition boundary through elements not wholly in partition
+%  (may include elements not in perimeter list)
+
+                    while 1
+%  if first edge, figure out direction from perimeter edge direction
+                        if ~nlast && ~slast
+                            nlast=find(elemp(ielem,:)==edgeper(j,1));
+                            nnext=find(elemp(ielem,:)==edgeper(j,2));
+                            if     (nlast+nnext == 3)
+                                slast=1;
+                            elseif (nlast+nnext == 5)
+                                slast=2;
+                            elseif (nlast+nnext == 4)
+                                slast=3;
+                            end
+                            if     (nnext+(6-nlast-nnext) == 3)
+                                snext=1;
+                            elseif (nnext+(6-nlast-nnext) == 5)
+                                snext=2;
+                            elseif (nnext+(6-nlast-nnext) == 4)
+                                snext=3;
+                            end
+
+%  find how many nodes of current element are in current partition
+%  (1 or 2, not 3) and handle each permutation separately
+                            ipart=find(epartp(ielem,:)==k);
+%  two nodes are in current partition, so cut off other node
+                            if (length(ipart) == 2)
+                                switch 6-sum(ipart)
+                                    case nlast
+                                        slast=6-slast-snext;
+                                        nlast=0;
+                                    case nnext
+                                        if (epartp(ielem,nnext) == k)
+                                            nlast=nnext;
+                                        end
+                                    otherwise
+                                        slast=6-slast-snext;
+                                        nlast=0;
+                                end
+%  one node is in current partition
+                            else
+%  all different, so cut through centroid
+                                if (epartp(ielem,1) ~= epartp(ielem,2)) && ...
+                                   (epartp(ielem,2) ~= epartp(ielem,3)) && ...
+                                   (epartp(ielem,3) ~= epartp(ielem,1))
+                                    switch ipart
+                                        case {nlast,nnext}
+                                            if (epartp(ielem,nnext) == k)
+                                                nlast=nnext;
+                                            end
+                                        otherwise
+                                            slast=6-slast-snext;
+                                            nlast=0;
+                                    end
+%  other two are in the same partition, so cut them off
+                                else
+                                    switch ipart
+                                        case nlast
+                                            if (epartp(ielem,nnext) == k)
+                                                nlast=nnext;
+                                            end
+                                        case nnext
+                                            slast=snext;
+                                            nlast=0;
+                                        otherwise
+                                            slast=6-slast-snext;
+                                            nlast=0;
+                                    end
+                                end
+                            end
+
+%  last edge exited last element at node
+                            if nlast
+%  write out first node of first side for half-edge to midpoint
+%                                disp(['segment j=' int2str(j) ' unshared half edge from node ' int2str(elemp(ielem,nlast)) ' (node ' int2str(nlast) ') on side ' int2str(slast) ' from element ' int2str(ielem) ' written.'])
+                                [lat(end+1),long(end+1)]=mapxy((md.x(elemp(ielem,nlast))),...
+                                                               (md.y(elemp(ielem,nlast))),'s');
+                                fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
+                            end
+                            nlast=0;
+                            
+%  write out midpoint of first side
+                            [lat(end+1),long(end+1)]=mapxy((md.x(elemp(ielem,slast))...
+                                                           +md.x(elemp(ielem,mod(slast,3)+1)))/2.,...
+                                                           (md.y(elemp(ielem,slast))...
+                                                           +md.y(elemp(ielem,mod(slast,3)+1)))/2.,'s');
+                            fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
+                        end
+
+%  last edge exited last element at node
+                        if nlast
+                            if elast
+%  find where last node on previous element occurs on current element
+                                nlast=find(elemp(ielem,:)==nlast,1);
+                            end
+%  half-edge occurs on unshared side from current node (unique unless mesh
+%  is only attached at node)
+                            switch nlast
+                                case 1
+                                    if ~edgeadjp(ielem,1)
+                                        nnext=2;
+                                        slast=1;
+                                    else
+                                        nnext=3;
+                                        slast=3;
+                                    end
+                                case 2
+                                    if ~edgeadjp(ielem,2)
+                                        nnext=3;
+                                        slast=2;
+                                    else
+                                        nnext=1;
+                                        slast=1;
+                                    end
+                                case 3
+                                    if ~edgeadjp(ielem,3)
+                                        nnext=1;
+                                        slast=3;
+                                    else
+                                        nnext=2;
+                                        slast=2;
+                                    end
+                            end
+%  write out half-edge from current node to midpoint of unshared side
+%                            disp(['segment j=' int2str(j) ' unshared half edge from node ' int2str(elemp(ielem,nlast)) ' (node ' int2str(nlast) ') on side ' int2str(slast) ' from element ' int2str(ielem) ' written.'])
+                            [lat(end+1),long(end+1)]=mapxy((md.x(elemp(ielem,nlast))...
+                                                           +md.x(elemp(ielem,nnext)))/2.,...
+                                                           (md.y(elemp(ielem,nlast))...
+                                                           +md.y(elemp(ielem,nnext)))/2.,'s');
+                            fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
+                            nlast=0;
+
+%  last edge exited last element at midpoint of side
+                        elseif slast
+                            if elast
+%  find where last side on previous element occurs on current element
+                                slast=find(edgeadjp(ielem,:)==elast,1);
+                            end
+                        end
+
+%  find how many nodes of current element are in current partition
+%  (1 or 2, not 3) and handle each permutation separately
+                        ipart=find(epartp(ielem,:)==k);
+                        if (length(ipart) == 2)
+%  two nodes are in current partition, so cut off other node
+                            switch 6-sum(ipart)
+                                case 1
+                                    snext=3+1-slast;
+                                case 2
+                                    snext=1+2-slast;
+                                case 3
+                                    snext=2+3-slast;
+                            end
+                        else
+                            if (epartp(ielem,1) ~= epartp(ielem,2)) && ...
+                               (epartp(ielem,2) ~= epartp(ielem,3)) && ...
+                               (epartp(ielem,3) ~= epartp(ielem,1))
+%  all different, so cut through centroid
+%                                disp(['element ielem=' int2str(ielem) ' centroid written.'])
+                                [lat(end+1),long(end+1)]=mapxy(sum(md.x(elemp(ielem,:)))/3.,...
+                                                               sum(md.y(elemp(ielem,:)))/3.,'s');
+                                fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
+                            end
+%  one node is in current partition, so cut off other two nodes
+                            switch ipart
+                                case 1
+                                    snext=3+1-slast;
+                                case 2
+                                    snext=1+2-slast;
+                                case 3
+                                    snext=2+3-slast;
+                            end
+                        end
+%  write out midpoint of opposite side
+%                        disp(['segment j=' int2str(j) ' internal edge from side ' int2str(slast) ' to side ' int2str(snext) ' from element ' int2str(ielem) ' written.'])
+                        [lat(end+1),long(end+1)]=mapxy((md.x(elemp(ielem,snext))...
+                                                       +md.x(elemp(ielem,mod(snext,3)+1)))/2.,...
+                                                       (md.y(elemp(ielem,snext))...
+                                                       +md.y(elemp(ielem,mod(snext,3)+1)))/2.,'s');
+                        fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
+                        elast=ielem;
+                        nlast=0;
+                        slast=snext;
+%  find adjacent element to opposite side
+                        ielem=edgeadjp(elast,slast);
+%  if opposite side is unshared, find it in edge perimeter list
+                        if ~ielem
+                            jlast=find(elemper(j:end)==elast)+j-1;
+                            j=0;
+                            for l=1:length(jlast)
+                                if ((elemp(elast,slast)          == edgeper(jlast(l),1)) && ...
+                                    (elemp(elast,mod(slast,3)+1) == edgeper(jlast(l),2))) || ...
+                                   ((elemp(elast,slast)          == edgeper(jlast(l),2)) && ...
+                                    (elemp(elast,mod(slast,3)+1) == edgeper(jlast(l),1)))
+                                    j=jlast(l);
+                                    break
+                                end
+                            end
+                            if ~j
+                                j=iloop(i+1)-1;
+                            end
+%  write out half-edge from midpoint of unshared side to node
+                            if (epartp(elast,slast) == k)
+                                nnext=slast;
+                            else
+                                nnext=mod(slast,3)+1;
+                            end
+%                            disp(['segment j=' int2str(j) ' unshared half edge on side ' int2str(slast) ' to node ' int2str(elemp(elast,nnext)) ' (node ' int2str(nnext) ') from element ' int2str(elast) ' written.'])
+                            [lat(end+1),long(end+1)]=mapxy(md.x(elemp(elast,nnext)),...
+                                                           md.y(elemp(elast,nnext)),'s');
+                            fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
+                            break
+%  if not unshared, advance perimeter list and watch for end
+                        else
+                            if (elast == elemper(j))
+                                if (j+1 < iloop(i+1)) && ...
+                                   ~isempty(find(elemper(j+1:end)~=elast,1))
+                                    j=j+find(elemper(j+1:end)~=elast,1);
+                                else
+                                    break
+                                end
+                            end
+                        end
+                    end
+                    j=j+1;
+                end
+            end
+%            fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(iloop(1)),lat(iloop(1)),alt);
+
+            fprintf(fid,'              </coordinates>\n');
+            fprintf(fid,'            </LinearRing>\n');
+            fprintf(fid,'          </outerBoundaryIs>\n');
+            fprintf(fid,'        </Polygon>\n');
+            fprintf(fid,'      </Placemark>\n');
+        end
+    end
+    fprintf(fid,'    </Folder>\n');
+end
+
+end
+
Index: /issm/trunk/src/m/kml/kml_unsh_edges.m
===================================================================
--- /issm/trunk/src/m/kml/kml_unsh_edges.m	(revision 6417)
+++ /issm/trunk/src/m/kml/kml_unsh_edges.m	(revision 6417)
@@ -0,0 +1,106 @@
+%%
+%  create kml linestrings for the unshared element edges.
+%
+%  []=kml_unsh_edges(fid,md,params)
+%
+%  where the required input is:
+%    fid           (char, file ID of .kml file)
+%    md            (model, model class object)
+%
+%  the optional input is:
+%    params        (string/numeric, parameter names and values)
+%
+%  and the optional input is:
+%    alt           (numeric, altitude for polygons, default 10000)
+%    prtplt        (char, 'off'/'no' for partition segment plot)
+%
+function []=kml_unsh_edges(varargin)
+
+if ~nargin
+    help kml_unsh_edges
+    return
+end
+
+%%  process input data
+
+iarg=1;
+if (nargin >= 1)
+    fid=varargin{1};
+end
+if ~exist('fid','var') || isempty(fid) || (fid < 0)
+    error('File ID ''%d'' must be open.',fid);
+end
+
+iarg=iarg+1;
+if (nargin >= 2)
+    md=varargin{2};
+end
+if ~exist('md','var') || isempty(md) || ~isa(md,'model')
+    error(['Model ''' inputname(iarg) ''' is unrecognized class ''' class(md) '''.']);
+end
+
+%  parameters
+
+iarg=iarg+1;
+while (iarg <= nargin-1)
+    if ischar(varargin{iarg})
+        eval([varargin{iarg} '=varargin{iarg+1};']);
+        if (numel(varargin{iarg+1}) <= 20)
+            disp([varargin{iarg} '=' any2str(varargin{iarg+1}) ';']);
+        else
+            disp([varargin{iarg} '=' string_size(varargin{iarg+1}) ' ' class(varargin{iarg+1}) ';']);
+        end
+    else
+        error(['''' any2str(varargin{iarg}) ''' is not a parameter name.']);
+    end
+    iarg=iarg+2;
+end
+
+if ~exist('alt','var')
+    alt=10000;
+end
+
+%%  write folder for unshared edges
+
+if (~exist('prtplt','var') || strncmpi(prtplt,'on' ,2) || strncmpi(prtplt,'y',1)) && ...
+    md.npart
+    [edgeadj]=edgeadjacency(md.elements,md.nodeconnectivity);
+    [icol,irow]=find(edgeadj'==0);
+    edgeuns=zeros(length(irow),2);
+    for i=1:length(irow)
+        edgeuns(i,1)=md.elements(irow(i),icol(i));
+        edgeuns(i,2)=md.elements(irow(i),mod(icol(i),size(md.elements,2))+1);
+    end
+    fprintf(fid,'    <Folder>\n');
+    fprintf(fid,'      <name>Unshared Edges</name>\n');
+    fprintf(fid,'      <visibility>1</visibility>\n');
+    fprintf(fid,'      <description>Partitions=%d, Edges=%d</description>\n',...
+        md.npart,size(edgeuns,1));
+
+%  write each edge as a linestring
+
+    disp(['Writing ' num2str(size(edgeuns,1)) ' unshared edges as KML linestrings.']);
+    for i=1:size(edgeuns,1)
+        fprintf(fid,'      <Placemark>\n');
+        fprintf(fid,'        <name>Edge %d</name>\n',i);
+        fprintf(fid,'        <visibility>1</visibility>\n');
+        fprintf(fid,'        <styleUrl>#RedLineRedPoly</styleUrl>\n');
+        fprintf(fid,'        <LineString>\n');
+        fprintf(fid,'          <extrude>1</extrude>\n');
+        fprintf(fid,'          <tessellate>1</tessellate>\n');
+        fprintf(fid,'          <altitudeMode>relativeToGround</altitudeMode>\n');
+        fprintf(fid,'          <coordinates>\n');
+        for j=1:2
+            [lat(j),long(j)]=mapxy(md.x(edgeuns(i,j)),md.y(edgeuns(i,j)),'s');
+            fprintf(fid,'            %0.16g,%0.16g,%0.16g\n',long(j),lat(j),alt);
+        end
+
+        fprintf(fid,'          </coordinates>\n');
+        fprintf(fid,'        </LineString>\n');
+        fprintf(fid,'      </Placemark>\n');
+    end
+    fprintf(fid,'    </Folder>\n');
+end
+
+end
+