Index: /issm/trunk/src/c/objects/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Tria.cpp	(revision 2354)
+++ /issm/trunk/src/c/objects/Tria.cpp	(revision 2355)
@@ -2394,7 +2394,20 @@
 		/*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
 		for (i=0;i<numgrids;i++){
-			grade_g_gaussian[i]=
-			  -2*drag*alpha_complement*((lambda*vx+mu*vy))*Jdet*gauss_weight*l1l2l3[i]                         //standard term dJ/dki
-			  -numpar->cm_noisedmp*Jdet*gauss_weight*(dh1dh2dh3_basic[0][i]*dk[0]+dh1dh2dh3_basic[1][i]*dk[1]);  // regularization term d/dki(1/2*(dk/dx)^2)
+
+			//standard term dJ/dki
+			grade_g_gaussian[i]=-2*drag*alpha_complement*((lambda*vx+mu*vy))*Jdet*gauss_weight*l1l2l3[i];
+
+			//noise dampening d/dki(1/2*(dk/dx)^2)
+			grade_g_gaussian[i]+=-numpar->cm_noisedmp*Jdet*gauss_weight*(dh1dh2dh3_basic[0][i]*dk[0]+dh1dh2dh3_basic[1][i]*dk[1]);
+			
+			//min dampening
+			if(drag<numpar->cm_mindmp_value){ 
+				grade_g_gaussian[i]+=numpar->cm_mindmp_slope*Jdet*gauss_weight*l1l2l3[i];
+			}
+
+			//max dampening
+			if(drag<numpar->cm_maxdmp_value){ 
+				grade_g_gaussian[i]+= - numpar->cm_maxdmp_slope*Jdet*gauss_weight*l1l2l3[i];
+			}
 		}
 		
@@ -2619,5 +2632,15 @@
 
 			//Add regularization term
-			grade_g_gaussian[i]-=numpar->cm_noisedmp*Jdet*gauss_weight*(dh1dh2dh3_basic[0][i]*dk[0]+dh1dh2dh3_basic[1][i]*dk[1]);
+			grade_g_gaussian[i]+= - numpar->cm_noisedmp*Jdet*gauss_weight*(dh1dh2dh3_basic[0][i]*dk[0]+dh1dh2dh3_basic[1][i]*dk[1]);
+
+			//min dampening
+			if(drag<numpar->cm_mindmp_value){ 
+				grade_g_gaussian[i]+= numpar->cm_mindmp_slope*Jdet*gauss_weight*l1l2l3[i];
+			}
+
+			//max dampening
+			if(drag<numpar->cm_maxdmp_value){ 
+				grade_g_gaussian[i]+= - numpar->cm_maxdmp_slope*Jdet*gauss_weight*l1l2l3[i];
+			}
 		}
 
@@ -2689,4 +2712,5 @@
 	double adjy_list[numgrids];
 	double adjxadjy_list[numgrids][2];
+	double B[numgrids];
 
 	int    dofs1[1]={0};
@@ -2725,5 +2749,6 @@
 	double  thickness;
 	int     dofs[1]={0};
-	double  dk[NDOF2]; 
+	double  dB[NDOF2]; 
+	double  B_gauss;
 	
 	ParameterInputs* inputs=NULL;
@@ -2747,4 +2772,7 @@
 	if(!inputs->Recover("adjoint",&adjxadjy_list[0][0],2,dofs2,numgrids,(void**)nodes)){
 		throw ErrorException(__FUNCT__,"missing adjoint input parameter");
+	}
+	if(!inputs->Recover("B",&B[0],1,dofs1,numgrids,(void**)nodes)){
+		throw ErrorException(__FUNCT__,"parameter B not found in input");
 	}
 
@@ -2800,6 +2828,7 @@
 		GetNodalFunctionsDerivativesBasic(&dh1dh2dh3_basic[0][0],&xyz_list[0][0],gauss_l1l2l3);
 
-		/*Get k derivative: dk/dx */
-		GetParameterDerivativeValue(&dk[0], &k[0],&xyz_list[0][0], gauss_l1l2l3);
+		/*Get B derivative: dB/dx */
+		GetParameterDerivativeValue(&dB[0], &B[0],&xyz_list[0][0], gauss_l1l2l3);
+		GetParameterValue(&B_gauss, &B[0],gauss_l1l2l3);
 
 		/*Build gradje_g_gaussian vector (actually -dJ/dB): */
@@ -2809,5 +2838,16 @@
 
 			//Add regularization term
-			grade_g_gaussian[i]-=numpar->cm_noisedmp*Jdet*gauss_weight*(dh1dh2dh3_basic[0][i]*dk[0]+dh1dh2dh3_basic[1][i]*dk[1]);
+			grade_g_gaussian[i]-=numpar->cm_noisedmp*Jdet*gauss_weight*(dh1dh2dh3_basic[0][i]*dB[0]+dh1dh2dh3_basic[1][i]*dB[1]);
+
+			//min dampening
+			if(B_gauss<numpar->cm_mindmp_value){ 
+				grade_g_gaussian[i]+= numpar->cm_mindmp_slope*Jdet*gauss_weight*l1l2l3[i];
+			}
+
+			//max dampening
+			if(B_gauss<numpar->cm_maxdmp_value){ 
+				grade_g_gaussian[i]+= - numpar->cm_maxdmp_slope*Jdet*gauss_weight*l1l2l3[i];
+			}
+
 		}
 
@@ -2866,4 +2906,6 @@
 	double  gauss_weight;
 	double  gauss_l1l2l3[3];
+	double  k_gauss;
+	double  B_gauss;
 
 	/* parameters: */
@@ -2963,9 +3005,22 @@
 		#endif
 		
-		/*Add dampening term to misfit*/
+		/*Add dampening terms to misfit*/
 		if (strcmp(numpar->control_type,"drag")==0){
 			if (!shelf){
+
+				//noise dampening
 				GetParameterDerivativeValue(&dk[0], &k[0],&xyz_list[0][0], gauss_l1l2l3);
 				Jelem+=numpar->cm_noisedmp*1/2*(pow(dk[0],2)+pow(dk[1],2))*Jdet*gauss_weight;
+
+				//min dampening
+				GetParameterValue(&k_gauss, &k[0],gauss_l1l2l3);
+				if(k_gauss<numpar->cm_mindmp_value){ 
+					Jelem+=numpar->cm_mindmp_slope*k_gauss*Jdet*gauss_weight;
+				}
+
+				//max dampening
+				if(k_gauss>numpar->cm_maxdmp_value){ 
+					Jelem+=numpar->cm_maxdmp_slope*k_gauss*Jdet*gauss_weight;
+				}
 			}
 		}
@@ -2974,6 +3029,18 @@
 				throw ErrorException(__FUNCT__,"parameter B not found in input");
 			}
+			//noise dampening
 			GetParameterDerivativeValue(&dB[0], &B[0],&xyz_list[0][0], gauss_l1l2l3);
 			Jelem+=numpar->cm_noisedmp*1/2*(pow(dB[0],2)+pow(dB[1],2))*Jdet*gauss_weight;
+
+			//min dampening
+			GetParameterValue(&B_gauss, &B[0],gauss_l1l2l3);
+			if(B_gauss<numpar->cm_mindmp_value){ 
+				Jelem+=numpar->cm_mindmp_slope*B_gauss*Jdet*gauss_weight;
+			}
+
+			//max dampening
+			if(B_gauss>numpar->cm_maxdmp_value){ 
+				Jelem+=numpar->cm_maxdmp_slope*B_gauss*Jdet*gauss_weight;
+			}
 		}
 		else{