Index: ../trunk-jpl/test/NightlyRun/test804.py
===================================================================
--- ../trunk-jpl/test/NightlyRun/test804.py	(revision 18815)
+++ ../trunk-jpl/test/NightlyRun/test804.py	(revision 18816)
@@ -24,6 +24,7 @@
 md.transient.iscalving=True;
 
 md.calving.calvingrate=1000.*numpy.ones((md.mesh.numberofvertices,1))
+md.calving.meltingrate=numpy.zeros((md.mesh.numberofvertices,1))
 
 md=solve(md,TransientSolutionEnum())
 
Index: ../trunk-jpl/test/NightlyRun/test804.m
===================================================================
--- ../trunk-jpl/test/NightlyRun/test804.m	(revision 18815)
+++ ../trunk-jpl/test/NightlyRun/test804.m	(revision 18816)
@@ -14,6 +14,7 @@
 md.transient.iscalving=1;
 
 md.calving.calvingrate=1000.*ones(md.mesh.numberofvertices,1);
+md.calving.meltingrate=zeros(md.mesh.numberofvertices,1);
 
 md=solve(md,TransientSolutionEnum());
 
Index: ../trunk-jpl/test/NightlyRun/test805.py
===================================================================
--- ../trunk-jpl/test/NightlyRun/test805.py	(revision 18815)
+++ ../trunk-jpl/test/NightlyRun/test805.py	(revision 18816)
@@ -31,6 +31,7 @@
 md.transient.iscalving=True;
 
 md.calving.calvingrate=1000.*numpy.ones((md.mesh.numberofvertices,1))
+md.calving.meltingrate=numpy.zeros((md.mesh.numberofvertices,1))
 
 md=solve(md,TransientSolutionEnum())
 
Index: ../trunk-jpl/test/NightlyRun/test805.m
===================================================================
--- ../trunk-jpl/test/NightlyRun/test805.m	(revision 18815)
+++ ../trunk-jpl/test/NightlyRun/test805.m	(revision 18816)
@@ -21,6 +21,7 @@
 md.transient.iscalving=1;
 
 md.calving.calvingrate=1000.*ones(md.mesh.numberofvertices,1);
+md.calving.meltingrate=zeros(md.mesh.numberofvertices,1);
 
 md=solve(md,TransientSolutionEnum());
 
Index: ../trunk-jpl/test/NightlyRun/test806.py
===================================================================
--- ../trunk-jpl/test/NightlyRun/test806.py	(revision 18815)
+++ ../trunk-jpl/test/NightlyRun/test806.py	(revision 18816)
@@ -36,8 +36,9 @@
 md.transient.isgia=False
 md.transient.iscalving=True;
 
-md.calving=calvinglevermann();
-md.calving.coeff=4.89e13*numpy.ones((md.mesh.numberofvertices,1));
+md.calving=calvinglevermann()
+md.calving.coeff=4.89e13*numpy.ones((md.mesh.numberofvertices,1))
+md.calving.meltingrate=numpy.zeros((md.mesh.numberofvertices,1))
 
 md.transient.requested_outputs=['default','StrainRateparallel','StrainRateperpendicular','Calvingratex','Calvingratey','CalvingCalvingrate']
 
Index: ../trunk-jpl/test/NightlyRun/test806.m
===================================================================
--- ../trunk-jpl/test/NightlyRun/test806.m	(revision 18815)
+++ ../trunk-jpl/test/NightlyRun/test806.m	(revision 18816)
@@ -25,6 +25,7 @@
 
 md.calving=calvinglevermann();
 md.calving.coeff=4.89e13*ones(md.mesh.numberofvertices,1);
+md.calving.meltingrate=zeros(md.mesh.numberofvertices,1);
 
 md.transient.requested_outputs={'default','StrainRateparallel','StrainRateperpendicular','Calvingratex','Calvingratey','CalvingCalvingrate'};
 
Index: ../trunk-jpl/test/NightlyRun/IdToName.m
===================================================================
--- ../trunk-jpl/test/NightlyRun/IdToName.m	(revision 18815)
+++ ../trunk-jpl/test/NightlyRun/IdToName.m	(revision 18816)
@@ -194,7 +194,8 @@
 	case 803, name='ValleyGlacierLevelsetEnthalpyHO3d';
 	case 804, name='ValleyGlacierLevelsetCalvingSIA2d';
 	case 805, name='ValleyGlacierLevelsetEnthCalvingHO3d';
-	case 806, name='SquareShelfLevelsetCalvingSIA2dLevermann';
+	case 806, name='SquareShelfLevelsetCalvingSSA2dLevermann';
+	case 807, name='SquareShelfLevelsetMeltingSSA2d';
 	case 1101, name='ISMIPAHO';
 	case 1102, name='ISMIPAFS';
 	case 1103, name='ISMIPBHO';
Index: ../trunk-jpl/test/NightlyRun/IdToName.py
===================================================================
--- ../trunk-jpl/test/NightlyRun/IdToName.py	(revision 18815)
+++ ../trunk-jpl/test/NightlyRun/IdToName.py	(revision 18816)
@@ -179,7 +179,8 @@
 	803  : 'ValleyGlacierLevelsetEnthalpyHO3d',	
 	804  : 'ValleyGlacierLevelsetCalvingSIA2d',	
 	805  : 'ValleyGlacierLevelsetEnthCalvingHO3d',	
-	806  : 'SquareShelfLevelsetCalvingSIA2dLevermann',	
+	806  : 'SquareShelfLevelsetCalvingSSA2dLevermann',	
+	807  : 'SquareShelfLevelsetMeltingSSA2d',
 	1101 : 'ISMIPAHO',
 	1102 : 'ISMIPAFS',
 	1103 : 'ISMIPBHO',
Index: ../trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp	(revision 18815)
+++ ../trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp	(revision 18816)
@@ -47,9 +47,11 @@
 		switch(calvinglaw){
 			case DefaultCalvingEnum:
 				iomodel->FetchDataToInput(elements,CalvingCalvingrateEnum);
+				iomodel->FetchDataToInput(elements,CalvingMeltingrateEnum);
 				break;
 			case CalvingLevermannEnum:
 				iomodel->FetchDataToInput(elements,CalvinglevermannCoeffEnum);
+				iomodel->FetchDataToInput(elements,CalvinglevermannMeltingrateEnum);
 				break;
 			default:
 				_error_("Calving law "<<EnumToStringx(calvinglaw)<<" not supported yet");
@@ -124,7 +126,7 @@
 	IssmDouble Jdet, dt, D_scalar;
 	IssmDouble h,hx,hy,hz;
 	IssmDouble vel;
-	IssmDouble norm_dlsf, calvingrate;
+	IssmDouble norm_dlsf, norm_calving, calvingrate, meltingrate;
 	IssmDouble* xyz_list = NULL;
 
 	/*Get problem dimension and whether there is calving or not*/
@@ -150,6 +152,7 @@
 	IssmDouble*    v        = xNew<IssmDouble>(dim);
 	IssmDouble*    w        = xNew<IssmDouble>(dim);
 	IssmDouble*    c        = xNewZeroInit<IssmDouble>(dim);
+	IssmDouble*    m        = xNewZeroInit<IssmDouble>(dim);
 	IssmDouble*    dlsf     = xNew<IssmDouble>(dim);
 
 	/*Retrieve all inputs and parameters*/
@@ -162,6 +165,7 @@
 	Input* lsf_slopex_input   = NULL;
 	Input* lsf_slopey_input   = NULL;
 	Input* calvingrate_input  = NULL;
+	Input* meltingrate_input  = NULL;
 
 	/*Load velocities*/
 	if(domaintype==Domain2DhorizontalEnum){
@@ -185,6 +189,7 @@
 				lsf_slopex_input  = basalelement->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
 				if(dim==2) lsf_slopey_input  = basalelement->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
 				calvingrate_input = basalelement->GetInput(CalvingCalvingrateEnum);     _assert_(calvingrate_input);
+				meltingrate_input = basalelement->GetInput(CalvingMeltingrateEnum);     _assert_(meltingrate_input);
 				break;
 			case CalvingLevermannEnum:
 				if(domaintype==Domain2DhorizontalEnum){
@@ -200,6 +205,7 @@
 						calvingratey_input=basalelement->GetInput(CalvingrateyAverageEnum); _assert_(calvingratey_input);
 					}
 				}
+				meltingrate_input = basalelement->GetInput(CalvinglevermannMeltingrateEnum);     _assert_(meltingrate_input);
 				break;
 			default:
 				_error_("Calving law "<<EnumToStringx(calvinglaw)<<" not supported yet");
@@ -237,19 +243,35 @@
 					lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
 					if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],gauss);
 					calvingrate_input->GetInputValue(&calvingrate,gauss);
+					meltingrate_input->GetInputValue(&meltingrate,gauss);
 
 					norm_dlsf=0.;
 					for(i=0;i<dim;i++) norm_dlsf+=pow(dlsf[i],2);
 					norm_dlsf=sqrt(norm_dlsf);
 
 					if(norm_dlsf>1.e-10)
-					 for(i=0;i<dim;i++) c[i]=calvingrate*dlsf[i]/norm_dlsf;
+					 for(i=0;i<dim;i++){ 
+					  c[i]=calvingrate*dlsf[i]/norm_dlsf;
+					  m[i]=meltingrate*dlsf[i]/norm_dlsf;
+					 }
 					else
-					 for(i=0;i<dim;i++) c[i]=0.;
+					 for(i=0;i<dim;i++){
+						 c[i]=0.;
+						 m[i]=0.;
+					 }
+
 					break;
 				case CalvingLevermannEnum:
 					calvingratex_input->GetInputValue(&c[0],gauss); // in 3D case, add mesh velocity 
 					if(dim==2) calvingratey_input->GetInputValue(&c[1],gauss);
+					meltingrate_input->GetInputValue(&meltingrate,gauss);
+
+					norm_calving=0.;
+					for(i=0;i<dim;i++) norm_calving+=pow(c[i],2);
+					norm_calving=sqrt(norm_calving)+1.e-14;
+
+					for(i=0;i<dim;i++) m[i]=meltingrate*c[i]/norm_calving;
+					
 					break;
 				default:
 					_error_("Calving law "<<EnumToStringx(calvinglaw)<<" not supported yet");
@@ -257,7 +279,7 @@
 		}
 
 		/*Levelset speed is ice velocity - calving rate*/
-		for(i=0;i<dim;i++) w[i]=v[i]-c[i];
+		for(i=0;i<dim;i++) w[i]=v[i]-c[i]-m[i];
 
 		/*Compute D*/
 		for(row=0;row<dim;row++){
Index: ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
===================================================================
--- ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 18815)
+++ ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 18816)
@@ -212,10 +212,12 @@
 	StressIntensityFactorEnum,
 	CalvingLawEnum,
 	CalvingCalvingrateEnum,
+	CalvingMeltingrateEnum,
 	CalvingLevermannEnum,
 	DefaultCalvingEnum,
 	CalvingRequestedOutputsEnum,
 	CalvinglevermannCoeffEnum,
+	CalvinglevermannMeltingrateEnum,
 	CalvingratexEnum,
 	CalvingrateyEnum,
 	CalvingratexAverageEnum,
Index: ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
===================================================================
--- ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 18815)
+++ ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 18816)
@@ -220,10 +220,12 @@
 		case StressIntensityFactorEnum : return "StressIntensityFactor";
 		case CalvingLawEnum : return "CalvingLaw";
 		case CalvingCalvingrateEnum : return "CalvingCalvingrate";
+		case CalvingMeltingrateEnum : return "CalvingMeltingrate";
 		case CalvingLevermannEnum : return "CalvingLevermann";
 		case DefaultCalvingEnum : return "DefaultCalving";
 		case CalvingRequestedOutputsEnum : return "CalvingRequestedOutputs";
 		case CalvinglevermannCoeffEnum : return "CalvinglevermannCoeff";
+		case CalvinglevermannMeltingrateEnum : return "CalvinglevermannMeltingrate";
 		case CalvingratexEnum : return "Calvingratex";
 		case CalvingrateyEnum : return "Calvingratey";
 		case CalvingratexAverageEnum : return "CalvingratexAverage";
Index: ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
===================================================================
--- ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 18815)
+++ ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 18816)
@@ -223,10 +223,12 @@
 	      else if (strcmp(name,"StressIntensityFactor")==0) return StressIntensityFactorEnum;
 	      else if (strcmp(name,"CalvingLaw")==0) return CalvingLawEnum;
 	      else if (strcmp(name,"CalvingCalvingrate")==0) return CalvingCalvingrateEnum;
+	      else if (strcmp(name,"CalvingMeltingrate")==0) return CalvingMeltingrateEnum;
 	      else if (strcmp(name,"CalvingLevermann")==0) return CalvingLevermannEnum;
 	      else if (strcmp(name,"DefaultCalving")==0) return DefaultCalvingEnum;
 	      else if (strcmp(name,"CalvingRequestedOutputs")==0) return CalvingRequestedOutputsEnum;
 	      else if (strcmp(name,"CalvinglevermannCoeff")==0) return CalvinglevermannCoeffEnum;
+	      else if (strcmp(name,"CalvinglevermannMeltingrate")==0) return CalvinglevermannMeltingrateEnum;
 	      else if (strcmp(name,"Calvingratex")==0) return CalvingratexEnum;
 	      else if (strcmp(name,"Calvingratey")==0) return CalvingrateyEnum;
 	      else if (strcmp(name,"CalvingratexAverage")==0) return CalvingratexAverageEnum;
@@ -257,12 +259,12 @@
 	      else if (strcmp(name,"MeshVertexonbase")==0) return MeshVertexonbaseEnum;
 	      else if (strcmp(name,"MeshVertexonsurface")==0) return MeshVertexonsurfaceEnum;
 	      else if (strcmp(name,"MeshVertexonboundary")==0) return MeshVertexonboundaryEnum;
-	      else if (strcmp(name,"MeshX")==0) return MeshXEnum;
-	      else if (strcmp(name,"MeshY")==0) return MeshYEnum;
          else stage=3;
    }
    if(stage==3){
-	      if (strcmp(name,"MeshZ")==0) return MeshZEnum;
+	      if (strcmp(name,"MeshX")==0) return MeshXEnum;
+	      else if (strcmp(name,"MeshY")==0) return MeshYEnum;
+	      else if (strcmp(name,"MeshZ")==0) return MeshZEnum;
 	      else if (strcmp(name,"MeshElementtype")==0) return MeshElementtypeEnum;
 	      else if (strcmp(name,"DomainType")==0) return DomainTypeEnum;
 	      else if (strcmp(name,"DomainDimension")==0) return DomainDimensionEnum;
@@ -380,12 +382,12 @@
 	      else if (strcmp(name,"SurfaceforcingsIssmbgradients")==0) return SurfaceforcingsIssmbgradientsEnum;
 	      else if (strcmp(name,"SolutionType")==0) return SolutionTypeEnum;
 	      else if (strcmp(name,"AnalysisType")==0) return AnalysisTypeEnum;
-	      else if (strcmp(name,"ConfigurationType")==0) return ConfigurationTypeEnum;
-	      else if (strcmp(name,"AdjointBalancethicknessAnalysis")==0) return AdjointBalancethicknessAnalysisEnum;
          else stage=4;
    }
    if(stage==4){
-	      if (strcmp(name,"AdjointBalancethickness2Analysis")==0) return AdjointBalancethickness2AnalysisEnum;
+	      if (strcmp(name,"ConfigurationType")==0) return ConfigurationTypeEnum;
+	      else if (strcmp(name,"AdjointBalancethicknessAnalysis")==0) return AdjointBalancethicknessAnalysisEnum;
+	      else if (strcmp(name,"AdjointBalancethickness2Analysis")==0) return AdjointBalancethickness2AnalysisEnum;
 	      else if (strcmp(name,"AdjointHorizAnalysis")==0) return AdjointHorizAnalysisEnum;
 	      else if (strcmp(name,"AnalysisCounter")==0) return AnalysisCounterEnum;
 	      else if (strcmp(name,"DefaultAnalysis")==0) return DefaultAnalysisEnum;
@@ -503,12 +505,12 @@
 	      else if (strcmp(name,"MassconaxpbyName")==0) return MassconaxpbyNameEnum;
 	      else if (strcmp(name,"MassconaxpbyNamex")==0) return MassconaxpbyNamexEnum;
 	      else if (strcmp(name,"MassconaxpbyNamey")==0) return MassconaxpbyNameyEnum;
-	      else if (strcmp(name,"MassconaxpbyAlpha")==0) return MassconaxpbyAlphaEnum;
-	      else if (strcmp(name,"MassconaxpbyBeta")==0) return MassconaxpbyBetaEnum;
          else stage=5;
    }
    if(stage==5){
-	      if (strcmp(name,"NodeSId")==0) return NodeSIdEnum;
+	      if (strcmp(name,"MassconaxpbyAlpha")==0) return MassconaxpbyAlphaEnum;
+	      else if (strcmp(name,"MassconaxpbyBeta")==0) return MassconaxpbyBetaEnum;
+	      else if (strcmp(name,"NodeSId")==0) return NodeSIdEnum;
 	      else if (strcmp(name,"VectorParam")==0) return VectorParamEnum;
 	      else if (strcmp(name,"Riftfront")==0) return RiftfrontEnum;
 	      else if (strcmp(name,"RiftfrontType")==0) return RiftfrontTypeEnum;
@@ -626,12 +628,12 @@
 	      else if (strcmp(name,"StrainRateyy")==0) return StrainRateyyEnum;
 	      else if (strcmp(name,"StrainRateyz")==0) return StrainRateyzEnum;
 	      else if (strcmp(name,"StrainRatezz")==0) return StrainRatezzEnum;
-	      else if (strcmp(name,"Divergence")==0) return DivergenceEnum;
-	      else if (strcmp(name,"MaxDivergence")==0) return MaxDivergenceEnum;
          else stage=6;
    }
    if(stage==6){
-	      if (strcmp(name,"GiaCrossSectionShape")==0) return GiaCrossSectionShapeEnum;
+	      if (strcmp(name,"Divergence")==0) return DivergenceEnum;
+	      else if (strcmp(name,"MaxDivergence")==0) return MaxDivergenceEnum;
+	      else if (strcmp(name,"GiaCrossSectionShape")==0) return GiaCrossSectionShapeEnum;
 	      else if (strcmp(name,"GiadWdt")==0) return GiadWdtEnum;
 	      else if (strcmp(name,"GiaW")==0) return GiaWEnum;
 	      else if (strcmp(name,"P0")==0) return P0Enum;
@@ -749,12 +751,12 @@
 	      else if (strcmp(name,"Dense")==0) return DenseEnum;
 	      else if (strcmp(name,"MpiDense")==0) return MpiDenseEnum;
 	      else if (strcmp(name,"MpiSparse")==0) return MpiSparseEnum;
-	      else if (strcmp(name,"Seq")==0) return SeqEnum;
-	      else if (strcmp(name,"Mpi")==0) return MpiEnum;
          else stage=7;
    }
    if(stage==7){
-	      if (strcmp(name,"Mumps")==0) return MumpsEnum;
+	      if (strcmp(name,"Seq")==0) return SeqEnum;
+	      else if (strcmp(name,"Mpi")==0) return MpiEnum;
+	      else if (strcmp(name,"Mumps")==0) return MumpsEnum;
 	      else if (strcmp(name,"Gsl")==0) return GslEnum;
 	      else if (strcmp(name,"Option")==0) return OptionEnum;
 	      else if (strcmp(name,"GenericOption")==0) return GenericOptionEnum;
Index: ../trunk-jpl/src/m/classes/calving.py
===================================================================
--- ../trunk-jpl/src/m/classes/calving.py	(revision 18815)
+++ ../trunk-jpl/src/m/classes/calving.py	(revision 18816)
@@ -14,6 +14,7 @@
 
 	def __init__(self): # {{{
 		self.calvingrate            = float('NaN')
+		self.meltingrate            = float('NaN')
 
 		#set defaults
 		self.setdefaultparameters()
@@ -22,6 +23,7 @@
 	def __repr__(self): # {{{
 		string='   Calving parameters:'
 		string="%s\n%s"%(string,fielddisplay(self,'calvingrate','calving rate at given location [m/a]'))
+		string="%s\n%s"%(string,fielddisplay(self,'meltingrate','melting rate at given location [m/a]'))
 
 		return string
 		#}}}
@@ -36,6 +38,7 @@
 			return md
 
 		md = checkfield(md,'fieldname','calving.calvingrate','NaN',1,'size',[md.mesh.numberofvertices],'>=',0)
+		md = checkfield(md,'fieldname','calving.meltingrate','NaN',1,'size',[md.mesh.numberofvertices],'>=',0)
 		return md
 	# }}}
 	def marshall(self,md,fid):    # {{{
@@ -44,4 +47,5 @@
 
 		WriteData(fid,'enum',CalvingLawEnum(),'data',DefaultCalvingEnum(),'format','Integer');
 		WriteData(fid,'object',self,'fieldname','calvingrate','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1,'scale',1./yts)
+		WriteData(fid,'object',self,'fieldname','meltingrate','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1,'scale',1./yts)
 	# }}}
Index: ../trunk-jpl/src/m/classes/calvinglevermann.py
===================================================================
--- ../trunk-jpl/src/m/classes/calvinglevermann.py	(revision 18815)
+++ ../trunk-jpl/src/m/classes/calvinglevermann.py	(revision 18816)
@@ -13,7 +13,8 @@
 	"""
 
 	def __init__(self): # {{{
-		self.coeff = float('NaN')
+		self.coeff       = float('NaN')
+		self.meltingrate = float('NaN')
 
 		#set defaults
 		self.setdefaultparameters()
@@ -22,6 +23,7 @@
 	def __repr__(self): # {{{
 		string='   Calving Levermann parameters:'
 		string="%s\n%s"%(string,fielddisplay(self,'coeff','proportionality coefficient in Levermann model'))
+		string="%s\n%s"%(string,fielddisplay(self,'meltingrate','melting rate at given location [m/a]'))
 
 		return string
 		#}}}
@@ -37,9 +39,12 @@
 			return md
 
 		md = checkfield(md,'fieldname','calving.coeff','size',[md.mesh.numberofvertices],'>',0)
+		md = checkfield(md,'fieldname','calving.meltingrate','NaN',1,'size',[md.mesh.numberofvertices],'>=',0)
 		return md
 	# }}}
 	def marshall(self,md,fid):    # {{{
+		yts=365.*24.*3600.
 		WriteData(fid,'enum',CalvingLawEnum(),'data',CalvingLevermannEnum(),'format','Integer');
 		WriteData(fid,'enum',CalvinglevermannCoeffEnum(),'data',self.coeff,'format','DoubleMat','mattype',1)
+		WriteData(fid,'object',self,'fieldname','meltingrate','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1,'scale',1./yts)
 	# }}}
Index: ../trunk-jpl/src/m/classes/calving.m
===================================================================
--- ../trunk-jpl/src/m/classes/calving.m	(revision 18815)
+++ ../trunk-jpl/src/m/classes/calving.m	(revision 18816)
@@ -6,6 +6,7 @@
 classdef calving
 	properties (SetAccess=public) 
 		 calvingrate            = NaN;
+		 meltingrate            = NaN;
 	end
 	methods (Static)
 		function obj = loadobj(obj) % {{{
@@ -59,15 +60,18 @@
 			%Early return
 			if (solution~=TransientSolutionEnum() | md.transient.iscalving==0), return; end
 			md = checkfield(md,'fieldname','calving.calvingrate','NaN',1,'size',[md.mesh.numberofvertices 1],'>=',0);
+			md = checkfield(md,'fieldname','calving.meltingrate','NaN',1,'size',[md.mesh.numberofvertices 1],'>=',0);
 		end % }}}
 		function disp(obj) % {{{
 			disp(sprintf('   Calving parameters:'));
 			fielddisplay(obj,'calvingrate','calving rate at given location [m/a]');
+			fielddisplay(obj,'meltingrate','melting rate at given location [m/a]');
 		end % }}}
 		function marshall(obj,md,fid) % {{{
 			yts=365.0*24.0*3600.0;
 			WriteData(fid,'enum',CalvingLawEnum(),'data',DefaultCalvingEnum(),'format','Integer');
 			WriteData(fid,'object',obj,'fieldname','calvingrate','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1,'scale',1./yts);
+			WriteData(fid,'object',obj,'fieldname','meltingrate','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1,'scale',1./yts);
 		end % }}}
 	end
 end
Index: ../trunk-jpl/src/m/classes/calvinglevermann.m
===================================================================
--- ../trunk-jpl/src/m/classes/calvinglevermann.m	(revision 18815)
+++ ../trunk-jpl/src/m/classes/calvinglevermann.m	(revision 18816)
@@ -5,7 +5,8 @@
 
 classdef calvinglevermann
 	properties (SetAccess=public) 
-		 coeff = NaN;
+		 coeff       = NaN;
+		 meltingrate = NaN;
 	end
 	methods
 		function obj = calvinglevermann(varargin) % {{{
@@ -34,15 +35,19 @@
 			%Early return
 			if (solution~=TransientSolutionEnum() | md.transient.iscalving==0), return; end
 			md = checkfield(md,'fieldname','calving.coeff','>',0,'size',[md.mesh.numberofvertices 1]);
+			md = checkfield(md,'fieldname','calving.meltingrate','NaN',1,'size',[md.mesh.numberofvertices 1],'>=',0);
 		end % }}}
 		function disp(obj) % {{{
 			disp(sprintf('   Calving Levermann parameters:'));
 			fielddisplay(obj,'coeff','proportionality coefficient in Levermann model');
+			fielddisplay(obj,'meltingrate','melting rate at given location [m/a]');
 
 		end % }}}
 		function marshall(obj,md,fid) % {{{
+			yts=365.0*24.0*3600.0;
 			WriteData(fid,'enum',CalvingLawEnum(),'data',CalvingLevermannEnum(),'format','Integer');
 			WriteData(fid,'enum',CalvinglevermannCoeffEnum(),'data',obj.coeff,'format','DoubleMat','mattype',1);
+			WriteData(fid,'object',obj,'fieldname','meltingrate','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1,'scale',1./yts);
 		end % }}}
 	end
 end
