Index: /issm/trunk/src/m/solutions/macayeal/diagnostic.m
===================================================================
--- /issm/trunk/src/m/solutions/macayeal/diagnostic.m	(revision 129)
+++ /issm/trunk/src/m/solutions/macayeal/diagnostic.m	(revision 130)
@@ -25,5 +25,4 @@
 x=md.x;
 y=md.y;
-
 index=md.elements;index=sort(index,2); %necessary 
 nods=md.numberofgrids;
@@ -47,5 +46,4 @@
 drag_type=md.drag_type;
 drag=md.drag;
-
 criterion_rel=md.eps_rel;
 criterion_abs=md.eps_abs;
@@ -53,4 +51,12 @@
 B=md.B; B_bar=(B(index(:,1))+B(index(:,2))+B(index(:,3)))/3;
 glen_coeff=md.n;
+
+%initialize velocities if any
+if (~isnan(md.vx) & ~isnan(md.vy)),
+	u=md.vx/yts; v=md.vy/yts;
+	velocity_is_present=1;
+else
+	velocity_is_present=0;
+end
 
 %average of p and q over the grids (size nel->nods)
@@ -266,8 +272,14 @@
 
 	%Compute viscosity (as in ICE and CIELO)
-	if convergence_count==1;
+	if (convergence_count==1),
 		%Initialize viscosity
-		nu_bar=viscosity(index,nel,alpha,beta,{},{},B_bar,glen_coeff);
-	elseif convergence_count==2,
+		if ~velocity_is_present;
+			nu_bar=viscosity(index,nel,alpha,beta,{},{},B_bar,glen_coeff);
+		else
+			nu_bar=viscosity(index,nel,alpha,beta,u,v,B_bar,glen_coeff);
+			nu_bar(find(nu_bar<=0))=-nu_bar(find(nu_bar<=0));
+			u_old=u; v_old=v;
+		end
+	elseif (convergence_count==2),
 		nu_bar_oldvalue=nu_bar;
 		nu_bar=viscosity(index,nel,alpha,beta,u,v,B_bar,glen_coeff);
@@ -421,5 +433,5 @@
 
 	%Test for direct shooting convergence
-	if convergence_count>1,
+	if (convergence_count>1 | velocity_is_present),
 
 		ug=[u_old;v_old];