Index: /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp	(revision 27202)
+++ /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp	(revision 27203)
@@ -225,4 +225,5 @@
 			parameters->AddObject(iomodel->CopyConstantObject("md.calving.xoffset",CalvingXoffsetEnum));
 			parameters->AddObject(iomodel->CopyConstantObject("md.calving.yoffset",CalvingYoffsetEnum));
+			parameters->AddObject(iomodel->CopyConstantObject("md.calving.vel_threshold",CalvingVelThresholdEnum));
 			break;
 		default:
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 27202)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 27203)
@@ -844,5 +844,5 @@
 	IssmDouble  arate, rho_ice, rho_water, thickness;
 	int			use_parameter=-1;
-	IssmDouble  gamma, theta, alpha, xoffset, yoffset;
+	IssmDouble  gamma, theta, alpha, xoffset, yoffset, vel_threshold, vrate;
 
 	/* Get node coordinates and dof list: */
@@ -870,4 +870,5 @@
 	this->FindParam(&xoffset, CalvingXoffsetEnum);
 	this->FindParam(&yoffset, CalvingYoffsetEnum);
+	this->FindParam(&vel_threshold, CalvingVelThresholdEnum);
 
 	/* Start looping on the number of vertices: */
@@ -883,7 +884,10 @@
 		gr_input->GetInputValue(&groundedice,gauss);
 		bs_input->GetInputValue(&bed,gauss);
+		H_input->GetInputValue(&thickness,gauss);
 		vel=sqrt(vx*vx+vy*vy)+1.e-14;
 		sl_input->GetInputValue(&sealevel,gauss);
 		arate_input->GetInputValue(&arate,gauss);
+		vrate = 1.0;
+		if (vel < vel_threshold) vrate = vel / vel_threshold;
 
 		/*Compute strain rate and viscosity: */
@@ -912,4 +916,8 @@
 				gamma = yoffset -  0.5*theta*tanh(alpha*(bed+xoffset));
 				break;
+			case 2:
+				/* 2 tanh(thicknes): f(x)=y_{o}-\frac{\theta}{2}\tanh(\alpha(x+x_{o})) */
+				gamma = yoffset -  0.5*theta*tanh(alpha*(-thickness+xoffset));
+				break;
 			case -1:
 				/* nothing, just the arate*/
@@ -926,5 +934,5 @@
 
 		/*-------------------------------------------*/
-		calvingrate[iv] = arate*gamma;
+		calvingrate[iv] = arate*gamma*vrate;
 	}
 
Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 27202)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 27203)
@@ -116,4 +116,5 @@
 	CalvingXoffsetEnum,
 	CalvingYoffsetEnum,
+	CalvingVelThresholdEnum,
 	ConfigurationTypeEnum,
 	ConstantsGEnum,
Index: /issm/trunk-jpl/src/m/classes/calvingparameterization.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/calvingparameterization.m	(revision 27202)
+++ /issm/trunk-jpl/src/m/classes/calvingparameterization.m	(revision 27203)
@@ -12,4 +12,5 @@
 		xoffset = 0;
 		yoffset = 0;
+		vel_threshold = 6000;
 	end
 	methods
@@ -51,4 +52,6 @@
 			% offset in y-axis 
 			self.yoffset = 0;
+			% velocity threshold to reduce calving rate
+			vel_threshold = 6000; % m/a
 		end % }}}
 		function md = checkconsistency(self,md,solution,analyses) % {{{
@@ -57,18 +60,20 @@
 
 			md = checkfield(md,'fieldname','calving.min_thickness','>=',0,'NaN',1,'Inf',1,'numel',1);
-			md = checkfield(md,'fieldname','calving.use_param','values',[-1, 0, 1]);
+			md = checkfield(md,'fieldname','calving.use_param','values',[-1, 0, 1, 2]);
 			md = checkfield(md,'fieldname','calving.theta','NaN',1,'Inf',1,'numel',1);
 			md = checkfield(md,'fieldname','calving.alpha','NaN',1,'Inf',1,'numel',1);
 			md = checkfield(md,'fieldname','calving.xoffset','NaN',1,'Inf',1,'numel',1);
 			md = checkfield(md,'fieldname','calving.yoffset','NaN',1,'Inf',1,'numel',1);
+			md = checkfield(md,'fieldname','calving.vel_threshold','NaN',1,'Inf',1,'numel',1);
 		end % }}}
 		function disp(self) % {{{
 			disp(sprintf('   Calving test parameters:'));
 			fielddisplay(self,'min_thickness','minimum thickness below which no ice is allowed [m]');
-			fielddisplay(self,'use_param','-1 - just use frontal ablation rate, 0 - f(x) = y_{o} + \alpha (x+x_{o}), 1 - f(x)=y_{o}-\frac{\theta}{2}\tanh(\alpha(x+x_{o}))');
+			fielddisplay(self,'use_param','-1 - just use frontal ablation rate, 0 - f(x) = y_{o} + \alpha (x+x_{o}), 1 - f(x)=y_{o}-\frac{\theta}{2}\tanh(\alpha(x+x_{o})), 2 - tanh(thickness)');
 			fielddisplay(self,'theta','the amplifier');
 			fielddisplay(self,'alpha','the slope');
 			fielddisplay(self,'xoffset','offset in x-axis');
 			fielddisplay(self,'yoffset','offset in y-axis');
+			fielddisplay(self,'vel_threshold','threshold of ice velocity to reduce the calving rate');
 		end % }}}
 		function marshall(self,prefix,md,fid) % {{{
@@ -81,4 +86,5 @@
 			WriteData(fid,prefix,'object',self,'fieldname','xoffset','format','Double');
 			WriteData(fid,prefix,'object',self,'fieldname','yoffset','format','Double');
+			WriteData(fid,prefix,'object',self,'fieldname','vel_threshold','format','Double','scale', 1./yts);
 		end % }}}
 	end