Index: /issm/trunk-jpl/src/m/Makefile.am
===================================================================
--- /issm/trunk-jpl/src/m/Makefile.am	(revision 22203)
+++ /issm/trunk-jpl/src/m/Makefile.am	(revision 22204)
@@ -65,4 +65,5 @@
 					${ISSM_DIR}/src/m/mesh/*.py \
 					${ISSM_DIR}/src/m/mesh/rifts/*.py \
+					${ISSM_DIR}/src/m/mesh/planet/gmsh/*.py \
 					${ISSM_DIR}/src/m/miscellaneous/*.py \
 					${ISSM_DIR}/src/m/modules/*.py \
Index: /issm/trunk-jpl/src/m/classes/frictionwaterlayer.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/frictionwaterlayer.py	(revision 22204)
+++ /issm/trunk-jpl/src/m/classes/frictionwaterlayer.py	(revision 22204)
@@ -0,0 +1,67 @@
+import numpy as np
+from project3d import project3d
+from fielddisplay import fielddisplay
+from checkfield import checkfield
+from WriteData import WriteData
+
+
+class frictionwaterlayer(object):
+	"""
+	frictionwaterlayer class definition
+
+	Usage:
+		friction=frictionwaterlayer(md);
+	"""
+	def __init__(self,md): # {{{
+		self.coefficient  		  	= md.friction.coefficient
+		self.f    				= float('NaN')
+		self.p    				= md.friction.p
+		self.q    				= md.friction.q
+		self.water_layer    			= float('NaN')
+	#}}}
+		
+	def checkconsistency(self,md,solution,analyses): #{{{ 
+		#Early return
+		if ('StressbalanceAnalysis' not in analyses) and ('ThermalAnalysis' not in analyses):
+			return
+
+		md = checkfield(md,'fieldname','friction.coefficient','timeseries',1,'NaN',1,'Inf',1)
+		md = checkfield(md,'fieldname','friction.f','size',[1,1],'NaN',1,'Inf',1)
+		md = checkfield(md,'fieldname','friction.q','NaN',1,'Inf',1,'size',[md.mesh.numberofelements])
+		md = checkfield(md,'fieldname','friction.p','NaN',1,'Inf',1,'size',[md.mesh.numberofelements])
+		md = checkfield(md,'fieldname','thermal.spctemperature','Inf',1,'timeseries',1,'>=',0.)
+	# }}}
+
+	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')
+		self.water_layer=project3d(md,'vector',self.water_layer,'type','node','layer',1)
+		return self	
+	 # }}} 
+
+
+	def __repr__(self): # {{{
+		string='Basal shear stress parameters: tau_b = coefficient^2 * Neff ^r * |u_b|^(s-1) * u_b * 1/f(T)\n(effective stress Neff=rho_ice*g*thickness+rho_water*g*(bed+water_layer), r=q/p and s=1/p)'
+		string="%s\n%s"%(string,fielddisplay(self,'coefficient','frictiontemp coefficient [SI]'))
+		string="%s\n%s"%(string,fielddisplay(self,'f','f variable for effective pressure'))
+		string="%s\n%s"%(string,fielddisplay(self,'p','p exponent'))
+		string="%s\n%s"%(string,fielddisplay(self,'q','q exponent'))
+		string="%s\n%s"%(string,fielddisplay(self,'water_layer','water thickness at the base of the ice (m)'))
+
+		return string
+	#}}}
+
+
+	def marshall(self,prefix,md,fid): #{{{
+		yts=md.constants.yts
+		
+		WriteData(fid,prefix,'name','md.friction.law','data',5,'format','Integer')
+		WriteData(fid,prefix,'class','friction','object',self,'fieldname','coefficient','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts)
+		WriteData(fid,prefix,'class','friction','object',self,'fieldname','f','format','Double')
+		WriteData(fid,prefix,'class','friction','object',self,'fieldname','p','format','DoubleMat','mattype',2)
+		WriteData(fid,prefix,'class','friction','object',self,'fieldname','q','format','DoubleMat','mattype',2)
+		WriteData(fid,prefix,'class','friction','object',self,'fieldname','water_layer','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts)
+	#}}}
+
+
Index: /issm/trunk-jpl/src/m/solve/WriteData.py
===================================================================
--- /issm/trunk-jpl/src/m/solve/WriteData.py	(revision 22203)
+++ /issm/trunk-jpl/src/m/solve/WriteData.py	(revision 22204)
@@ -44,15 +44,14 @@
 		data=np.array(data)
 		scale = options.getfieldvalue('scale')
-		if np.size(data) > 1 :
-			if np.size(data,0)==timeserieslength:
-				data[0:-1,:] = scale*data[0:-1,:]
-			else:
-				data  = scale*data
+		if np.size(data) > 1 and np.ndim(data) > 1 and np.size(data,0)==timeserieslength:
+			data[0:-1,:] = scale*data[0:-1,:]
 		else:
 			data  = scale*data
-	if np.size(data) > 1 :
-		if np.size(data,0)==timeserieslength:
-			yts = options.getfieldvalue('yts')
-			data[-1,:] = yts*data[-1,:]
+	if np.size(data) > 1 and np.size(data,0)==timeserieslength:
+		yts = options.getfieldvalue('yts')
+		if np.ndim(data) > 1:
+			data[-1,:] = data[-1,:]*yts
+		else:
+			data[-1] = data[-1]*yts
 
 	#Step 1: write the enum to identify this record uniquely