Index: /issm/trunk-jpl/src/m/classes/friction.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/friction.py	(revision 19706)
+++ /issm/trunk-jpl/src/m/classes/friction.py	(revision 19707)
@@ -10,5 +10,5 @@
 
 	   Usage:
-	      friction=friction();
+	      friction=friction()
 	"""
 
Index: /issm/trunk-jpl/src/m/classes/frictioncoulomb.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/frictioncoulomb.m	(revision 19706)
+++ /issm/trunk-jpl/src/m/classes/frictioncoulomb.m	(revision 19707)
@@ -48,12 +48,11 @@
 		end % }}}
 		function disp(self) % {{{
-			disp(sprintf('Basal shear stress parameters: Sigma_b = min( coefficient^2 * Neff ^r * |u_b|^(s-1) * u_b\n, coefficientcoulomb^2 * rho_i * g * (h-h_f) (effective stress Neff=rho_ice*g*thickness+rho_water*g*bed, r=q/p and s=1/p, floatation thickness h_f=max(0,-rho_fw / rho_i * bed))'));
-			fielddisplay(self,'coefficient','friction coefficient [SI]');
-			fielddisplay(self,'coefficientcoulomb','coulomb friction coefficient [SI]');
+			disp(sprintf('Basal shear stress parameters: Sigma_b = min( coefficient^2 * Neff ^r * |u_b|^(s-1) * u_b\n, coefficientcoulomb^2 * rho_i * g * (h-h_f)) (effective stress Neff=rho_ice*g*thickness+rho_water*g*bed, r=q/p and s=1/p, floatation thickness h_f=max(0,-rho_sw / rho_i * bed))'));
+			fielddisplay(self,'coefficient','power law (Weertman) friction coefficient [SI]');
+			fielddisplay(self,'coefficientcoulomb','Coulomb friction coefficient [SI]');
 			fielddisplay(self,'p','p exponent');
 			fielddisplay(self,'q','q exponent');
 		end % }}}
 		function marshall(self,md,fid) % {{{
-			yts=365.0*24.0*3600.0;
 
 			WriteData(fid,'enum',FrictionLawEnum,'data',7,'format','Integer');
Index: /issm/trunk-jpl/src/m/classes/frictioncoulomb.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/frictioncoulomb.py	(revision 19707)
+++ /issm/trunk-jpl/src/m/classes/frictioncoulomb.py	(revision 19707)
@@ -0,0 +1,63 @@
+from fielddisplay import fielddisplay
+from project3d import project3d
+from EnumDefinitions import *
+from checkfield import checkfield
+from WriteData import WriteData
+
+class frictioncoulomb(object):
+    """
+    FRICTIONCOULOMB class definition
+
+    Usage:
+        frictioncoulomb=frictioncoulomb()
+    """
+
+    def __init__(self): # {{{
+        self.coefficient = float('NaN')
+        self.coefficientcoulomb = float('NaN')
+        self.p = float('NaN')
+	self.q = float('NaN')
+
+	#set defaults
+	self.setdefaultparameters()
+
+    #}}}
+    def __repr__(self): # {{{
+	string="Basal shear stress parameters: Sigma_b = min(coefficient^2 * Neff ^r * |u_b|^(s-1) * u_b,\n coefficientcoulomb^2 * rho_i * g * (h-h_f)), (effective stress Neff=rho_ice*g*thickness+rho_water*g*bed, r=q/p and s=1/p)."
+
+	string="%s\n%s"%(string,fielddisplay(self,"coefficient","power law (Weertman) friction coefficient [SI]"))
+	string="%s\n%s"%(string,fielddisplay(self,"coefficientcoulomb","Coulomb friction coefficient [SI]"))
+	string="%s\n%s"%(string,fielddisplay(self,"p","p exponent"))
+	string="%s\n%s"%(string,fielddisplay(self,"q","q exponent"))
+	return string
+    #}}}
+    def extrude(self,md): # {{{
+	self.coefficient=project3d(md,'vector',self.coefficient,'type','node','layer',1)
+	self.coefficientcoulomb=project3d(md,'vector',self.coefficientcoulomb,'type','node','layer',1)
+	self.p=project3d(md,'vector',self.p,'type','element')
+	self.q=project3d(md,'vector',self.q,'type','element')
+	return self
+    #}}}
+    def setdefaultparameters(self): # {{{
+	return self
+    #}}}
+    def checkconsistency(self,md,solution,analyses):    # {{{
+
+	#Early return
+	if StressbalanceAnalysisEnum() not in analyses and ThermalAnalysisEnum() not in analyses:
+	    return md
+
+	md = checkfield(md,'fieldname','friction.coefficient','timeseries',1,'NaN',1)
+	md = checkfield(md,'fieldname','friction.coefficientcoulomb','timeseries',1,'NaN',1)
+	md = checkfield(md,'fieldname','friction.q','NaN',1,'size',[md.mesh.numberofelements])
+	md = checkfield(md,'fieldname','friction.p','NaN',1,'size',[md.mesh.numberofelements])
+
+	return md
+    # }}}
+    def marshall(self,md,fid):    # {{{
+	WriteData(fid,'enum',FrictionLawEnum(),'data',1,'format','Integer')
+	WriteData(fid,'object',self,'fieldname','coefficient','format','DoubleMat','mattype',1)
+	WriteData(fid,'object',self,'fieldname','coefficientcoulomb','format','DoubleMat','mattype',1)
+	WriteData(fid,'object',self,'fieldname','p','format','DoubleMat','mattype',2)
+	WriteData(fid,'object',self,'fieldname','q','format','DoubleMat','mattype',2)
+    # }}}
Index: /issm/trunk-jpl/src/m/classes/mismipbasalforcings.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/mismipbasalforcings.m	(revision 19706)
+++ /issm/trunk-jpl/src/m/classes/mismipbasalforcings.m	(revision 19707)
@@ -72,5 +72,5 @@
 		end % }}}
 		function disp(self) % {{{
-			disp(sprintf('   mismip basal forcings parameters:'));
+			disp(sprintf('   MISMIP+ basal melt parameterization:'));
 
 			fielddisplay(self,'groundedice_melting_rate','basal melting rate (positive if melting) [m/yr]');
Index: /issm/trunk-jpl/src/m/classes/mismipbasalforcings.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/mismipbasalforcings.py	(revision 19707)
+++ /issm/trunk-jpl/src/m/classes/mismipbasalforcings.py	(revision 19707)
@@ -0,0 +1,98 @@
+from fielddisplay import fielddisplay
+from EnumDefinitions import *
+from checkfield import checkfield
+from WriteData import WriteData
+import numpy
+
+class mismipbasalforcings(object):
+    """
+    MISMIP Basal Forcings class definition
+
+        Usage:
+	    mismipbasalforcings=mismipbasalforcings()
+    """
+
+    def __init__(self,md): # {{{
+
+        self.groundedice_melting_rate = float('NaN')
+        self.meltrate_factor = float('NaN')
+        self.threshold_thickness = float('NaN')
+        self.upperdepth_melt = float('NaN')
+        self.geothermalflux = float('NaN')
+
+        if numpy.all(numpy.isnan(self.groundedice_melting_rate)):
+            self.groundedice_melting_rate=numpy.zeros(md.mesh.numberofvertices)
+            print ' no basalforcings.groundedice_melting_rate specified: values set as zero'
+
+	self.setdefaultparameters()
+
+    #}}}
+    def __repr__(self): # {{{
+        string=" MISMIP+ basal melt parameterization\n"
+        string="%s\n%s"%(string,fielddisplay(self,"groundedice_melting_rate","basal melting rate (positive if melting) [m/yr]"))
+        string="%s\n%s"%(string,fielddisplay(self,"meltrate_factor","Melt-rate rate factor [1/yr]"))
+        string="%s\n%s"%(string,fielddisplay(self,"threshold_thickness","Threshold thickness for saturation of basal melting [m]"))
+        string="%s\n%s"%(string,fielddisplay(self,"upperdepth_melt","Depth above which melt rate is zero [m]"))
+        string="%s\n%s"%(string,fielddisplay(self,"geothermalflux","Geothermal heat flux [W/m^2]"))
+
+	return string
+    #}}}
+    def extrude(self,md): # {{{
+	self.coefficient=project3d(md,'vector',self.coefficient,'type','node','layer',1)
+	self.p=project3d(md,'vector',self.p,'type','element')
+	self.q=project3d(md,'vector',self.q,'type','element')
+	return self
+    #}}}
+    def setdefaultparameters(self): # {{{
+
+        # default values for melting parameterization
+        self.meltrate_factor = 0.2
+        self.threshold_thickness = 75.
+        self.upperdepth_melt = -100.
+
+	return self
+    #}}}
+    def checkconsistency(self,md,solution,analyses):    # {{{
+
+	#Early return
+        if MasstransportAnalysisEnum() in analyses and not (solution==TransientSolutionEnum() and md.transient.ismasstransport==0):
+
+	    md = checkfield(md,'fieldname','basalforcings.groundedice_melting_rate','NaN',1,'timeseries',1)
+	    md = checkfield(md,'fieldname','basalforcings.meltrate_factor','>=',0,'numel',1)
+	    md = checkfield(md,'fieldname','basalforcings.threshold_thickness','>=',0,'numel',1)
+	    md = checkfield(md,'fieldname','basalforcings.upperdepth_melt','<=',0,'numel',1)
+
+        if BalancethicknessAnalysisEnum() in analyses:
+
+	    md = checkfield(md,'fieldname','basalforcings.groundedice_melting_rate','NaN',1,'size',[md.mesh.numberofvertices])
+	    md = checkfield(md,'fieldname','basalforcings.meltrate_factor','>=',0,'numel',1)
+	    md = checkfield(md,'fieldname','basalforcings.threshold_thickness','>=',0,'numel',1)
+	    md = checkfield(md,'fieldname','basalforcings.upperdepth_melt','<=',0,'numel',1)
+
+        if ThermalAnalysisEnum() in analyses and not (solution==TransientSolutionEnum() and md.transient.isthermal==0):
+
+	    md = checkfield(md,'fieldname','basalforcings.groundedice_melting_rate','NaN',1,'timeseries',1)
+	    md = checkfield(md,'fieldname','basalforcings.meltrate_factor','>=',0,'numel',1)
+	    md = checkfield(md,'fieldname','basalforcings.threshold_thickness','>=',0,'numel',1)
+	    md = checkfield(md,'fieldname','basalforcings.upperdepth_melt','<=',0,'numel',1)
+	    md = checkfield(md,'fieldname','basalforcings.geothermal_flux','NaN',1,'timeseries',1,'>=',0)
+	return md
+    # }}}
+    def marshall(self,md,fid):    # {{{
+
+        yts=md.constants.yts
+        if yts!=365.2422*24.*3600.:
+            print 'WARNING: value of yts for MISMIP+ runs different from ISSM default!'
+
+        floatingice_melting_rate = numpy.zeros((md.mesh.numberofvertices,1))
+        floatingice_melting_rate = md.basalforcings.meltrate_factor*numpy.tanh((md.geometry.base-md.geometry.bed)/md.basalforcings.threshold_thickness)*numpy.amax(md.basalforcings.upperdepth_melt-md.geometry.base,0)
+
+	WriteData(fid,'enum',BasalforcingsEnum(),'data',MismipFloatingMeltRateEnum(),'format','Integer')
+	WriteData(fid,'data',floatingice_melting_rate,'format','DoubleMat','enum',BasalforcingsFloatingiceMeltingRateEnum(),'mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1)
+	WriteData(fid,'object',self,'fieldname','groundedice_melting_rate','format','DoubleMat','enum',BasalforcingsGroundediceMeltingRateEnum(),'mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1)
+	WriteData(fid,'object',self,'fieldname','geothermal_flux','enum',BasalforcingsGeothermalfluxEnum(),'format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1)
+	WriteData(fid,'object',self,'fieldname','meltrate_factor','format','Double','enum',BasalforcingsMeltrateFactorEnum(),'scale',1./yts)
+	WriteData(fid,'object',self,'fieldname','threshold_thickness','format','Double','enum',BasalforcingsThresholdThicknessEnum())
+	WriteData(fid,'object',self,'fieldname','upperdepth_melt','format','Double','enum',BasalforcingsUpperdepthMeltEnum())
+
+    # }}}
Index: /issm/trunk-jpl/src/m/classes/timestepping.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/timestepping.py	(revision 19706)
+++ /issm/trunk-jpl/src/m/classes/timestepping.py	(revision 19707)
@@ -66,5 +66,5 @@
 	def marshall(self,md,fid):    # {{{
 
-		yts=365.0*24.0*3600.0
+                yts=md.constants.yts
 
 		WriteData(fid,'object',self,'fieldname','start_time','format','Double','scale',yts)