Index: /issm/trunk/src/c/modules/MeshProfileIntersectionx/ElementSegment.cpp
===================================================================
--- /issm/trunk/src/c/modules/MeshProfileIntersectionx/ElementSegment.cpp	(revision 4993)
+++ /issm/trunk/src/c/modules/MeshProfileIntersectionx/ElementSegment.cpp	(revision 4994)
@@ -85,3 +85,9 @@
 		segments_dataset->AddObject(new  Segment(el+1,xfinal[0],yfinal[0],xfinal[1],yfinal[1]));
 	}
+	else{
+		/*No interesections, but the segment might be entirely inside this triangle!: */
+		if ( (NodeInElement(xgrids,ygrids,xsegment[0],ysegment[0])) && (NodeInElement(xgrids,ygrids,xsegment[1],ysegment[1])) ){
+			segments_dataset->AddObject(new  Segment(el+1,xsegment[0],ysegment[0],xsegment[1],ysegment[1]));
+		}
+	}
 }
Index: /issm/trunk/src/c/modules/ModelProcessorx/Qmu/CreateParametersQmu.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/Qmu/CreateParametersQmu.cpp	(revision 4993)
+++ /issm/trunk/src/c/modules/ModelProcessorx/Qmu/CreateParametersQmu.cpp	(revision 4994)
@@ -151,6 +151,7 @@
 			descriptor=variabledescriptors[i];
 
-			if ((strcmp(descriptor,"Thickness")==0) ||
-				(strcmp(descriptor,"DragCoefficient")     ==0)
+			if ((strcmp(descriptor,"Thickness")==0) || (strcmp(descriptor,"DragCoefficient")     ==0) || (strcmp(descriptor,"Surface")     ==0) || (strcmp(descriptor,"Bed")     ==0)
+
+
 				){
 
Index: /issm/trunk/src/m/classes/public/plot/plot_manager.m
===================================================================
--- /issm/trunk/src/m/classes/public/plot/plot_manager.m	(revision 4993)
+++ /issm/trunk/src/m/classes/public/plot/plot_manager.m	(revision 4994)
@@ -45,4 +45,8 @@
 			plot_gridnumbering(md,options,subplotwidth,i);
 			return;
+		case 'qmu_mass_flux_segments',
+			plot_qmu_mass_flux_segments(md,options,nlines,ncols,i);
+			return;
+
 		case 'highlightgrids',
 			plot_highlightgrids(md,options,subplotwidth,i);
Index: /issm/trunk/src/m/classes/public/plot/plot_qmu_mass_flux_segments.m
===================================================================
--- /issm/trunk/src/m/classes/public/plot/plot_qmu_mass_flux_segments.m	(revision 4994)
+++ /issm/trunk/src/m/classes/public/plot/plot_qmu_mass_flux_segments.m	(revision 4994)
@@ -0,0 +1,49 @@
+function plot_qmu_mass_flux_segments(md,options,nlines,ncols,i);
+%PLOT_QMU_MASS_FLUX_SEGMENTS - plot segments from the qmu analysis of mass fluxes
+%
+%   Usage:
+%      plot_qmu_mass_flux_segments(md,options,nlines,ncols,i);
+%
+
+subplot(nlines,ncols,i); 
+
+%process mesh and data
+[x y z elements is2d]=processmesh(md,[],options);
+
+allsegments=md.qmu_mass_flux_segments;
+
+if (md.dim==2),
+
+	%recover segments
+	hold on
+	for i=1:length(allsegments),
+		segments=allsegments{i};
+
+		%plot semgnets
+		for j=1:length(segments),
+			plot([segments(j,1) segments(j,3)],[segments(j,2) segments(j,4)]);
+		end
+		text(segments(j,1),segments(j,2),['Profile #' num2str(i)]);
+
+		%plot normals
+		
+		for j=1:length(segments),
+			xstart=mean([segments(j,1) segments(j,3)]);
+			ystart=mean([segments(j,2) segments(j,4)]);
+			length1=sqrt((segments(j,1)-segments(j,3)).^2 + (segments(j,2)-segments(j,4)).^2);
+			normal(:,1)=cos(atan2(segments(j,1)-segments(j,3) , segments(j,4)-segments(j,2)));
+			normal(:,2)=sin(atan2(segments(j,1)-segments(j,3) , segments(j,4)-segments(j,2)));
+			xend=xstart+length1.*normal(:,1);
+			yend=ystart+length1.*normal(:,2);
+			plot([xstart xend],[ystart yend],'r-');
+		end
+
+	end
+else
+	error('plot_qmu_mass_flux_segments: 3d plot of segments not supported yet!');
+end
+
+%apply options
+options=addfielddefault(options,'title','Mass Flux segments and normals');
+options=addfielddefault(options,'colorbar',0);
+applyoptions(md,[],options);
Index: /issm/trunk/src/m/qmu/setupdesign/QmuSetupDesign.m
===================================================================
--- /issm/trunk/src/m/qmu/setupdesign/QmuSetupDesign.m	(revision 4993)
+++ /issm/trunk/src/m/qmu/setupdesign/QmuSetupDesign.m	(revision 4994)
@@ -29,4 +29,12 @@
 	dvar=setupthickness(dvar,variables,params,varargin{:});
 
+elseif strcmpi(descriptor,'Surface')
+
+	dvar=setupsurface(dvar,variables,params,varargin{:});
+
+elseif strcmpi(descriptor,'Bed')
+
+	dvar=setupbed(dvar,variables,params,varargin{:});
+
 elseif strcmpi(descriptor,'DragCoefficient')
 
Index: /issm/trunk/src/m/qmu/setupdesign/setupbed.m
===================================================================
--- /issm/trunk/src/m/qmu/setupdesign/setupbed.m	(revision 4994)
+++ /issm/trunk/src/m/qmu/setupdesign/setupbed.m	(revision 4994)
@@ -0,0 +1,27 @@
+function dvar=setupbed(dvar,variables,params,varargin)
+
+for i=1:length(varargin)
+    if strcmp(class(varargin{i}),'model')
+        md=varargin{i};
+        break;
+    end
+end
+%capture stddev:
+stddev=variables.stddev;
+
+if length(stddev)>md.npart,
+	error('setupbed error message: stddev should be either a scalar or a ''npart'' length vector');
+end
+
+%ok, dealing with semi-discrete distributed variable. Distribute according to how many 
+%partitions we want
+
+for j=1:md.npart
+	dvar(end+1)           =variables;
+	dvar(end  ).descriptor=sprintf('%s%d',variables.descriptor,j);
+	if length(stddev)>1,
+		dvar(end  ).stddev=stddev(j);
+	end
+end
+	
+end
Index: /issm/trunk/src/m/qmu/setupdesign/setupsurface.m
===================================================================
--- /issm/trunk/src/m/qmu/setupdesign/setupsurface.m	(revision 4994)
+++ /issm/trunk/src/m/qmu/setupdesign/setupsurface.m	(revision 4994)
@@ -0,0 +1,27 @@
+function dvar=setupsurface(dvar,variables,params,varargin)
+
+for i=1:length(varargin)
+    if strcmp(class(varargin{i}),'model')
+        md=varargin{i};
+        break;
+    end
+end
+%capture stddev:
+stddev=variables.stddev;
+
+if length(stddev)>md.npart,
+	error('setupsurface error message: stddev should be either a scalar or a ''npart'' length vector');
+end
+
+%ok, dealing with semi-discrete distributed variable. Distribute according to how many 
+%partitions we want
+
+for j=1:md.npart
+	dvar(end+1)           =variables;
+	dvar(end  ).descriptor=sprintf('%s%d',variables.descriptor,j);
+	if length(stddev)>1,
+		dvar(end  ).stddev=stddev(j);
+	end
+end
+	
+end
Index: /issm/trunk/src/m/utils/Exp/expcreatecircle.m
===================================================================
--- /issm/trunk/src/m/utils/Exp/expcreatecircle.m	(revision 4993)
+++ /issm/trunk/src/m/utils/Exp/expcreatecircle.m	(revision 4994)
@@ -11,8 +11,9 @@
 
 %Calculate the cartesians coordinates of the points
-x_list=ones(numberofgrids,1);
-y_list=ones(numberofgrids,1);
+x_list=ones(numberofgrids+1,1);
+y_list=ones(numberofgrids+1,1);
 
 theta=(0:2*pi/numberofgrids:2*pi*(1-1/numberofgrids))';
+theta=[theta;0];
 
 x_list=radius*x_list.*cos(theta);