Index: /issm/trunk/src/m/solutions/cielo/GradJCompute.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/GradJCompute.m	(revision 1200)
+++ /issm/trunk/src/m/solutions/cielo/GradJCompute.m	(revision 1201)
@@ -1,11 +1,18 @@
 function [u_g grad_g]=GradJCompute(m,inputs,analysis_type,sub_analysis_type);
 
+%recover parameters
+debug=m.parameters.debug;
+dim=m.parameters.dim;
+ishutter=m.parameters.ishutter;
+ismacayealpattyn=m.parameters.ismacayealpattyn;
+isstokes=m.parameters.isstokes;
+
 %Recover solution for this stiffness and right hand side: 
-displaystring(m.parameters.debug,'%s','         computing velocities...');
+displaystring(debug,'%s','         computing velocities...');
 [u_g K_ff0 K_fs0 ]=diagnostic_core_nonlinear(m,inputs,analysis_type,sub_analysis_type);
 inputs=add(inputs,'velocity',u_g,'doublevec',2,m.parameters.numberofnodes);
 
 %Buid Du, difference between observed velocity and model velocity.
-displaystring(m.parameters.debug,'%s','          computing Du...');
+displaystring(debug,'%s','          computing Du...');
 [Du_g]=Du(m.elements,m.nodes,m.loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type);
 
@@ -14,5 +21,5 @@
 
 %Solve for adjoint vector: 
-displaystring(m.parameters.debug,'%s','          computing adjoint state...');
+displaystring(debug,'%s','          computing adjoint state...');
 lambda_f=Solver(K_ff0,Du_f,[],m.parameters);
 
@@ -23,2 +30,7 @@
 %Compute gradJ 
 grad_g=Gradj(m.elements,m.nodes,m.loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type);
+
+if dim==3,
+	displaystring(debug,'%s','          extruding gradient...');
+	grad_g=FieldExtrude(m.elements,m.nodes,m.loads,m.materials,grad_g,'gradj',0);
+end