Index: /issm/trunk-jpl/src/m/boundaryconditions/SetIceSheetBC.m
===================================================================
--- /issm/trunk-jpl/src/m/boundaryconditions/SetIceSheetBC.m	(revision 17070)
+++ /issm/trunk-jpl/src/m/boundaryconditions/SetIceSheetBC.m	(revision 17071)
@@ -29,22 +29,13 @@
 %No ice front: do nothing
 
-%Create zeros basal melting rate and surface mass balance if not specified
-if (isnan(md.surfaceforcings.precipitation) & (md.surfaceforcings.ispdd==1)),
-	md.surfaceforcings.precipitation=zeros(md.mesh.numberofvertices,1);
-	disp('      no surfaceforcings.precipitation specified: values set as zero');
-end
-if (isnan(md.surfaceforcings.mass_balance) & (md.surfaceforcings.ispdd==0)),
-	md.surfaceforcings.mass_balance=zeros(md.mesh.numberofvertices,1);
-	disp('      no surfaceforcings.mass_balance specified: values set as zero');
-end
-if isnan(md.basalforcings.melting_rate),
-	md.basalforcings.melting_rate=zeros(md.mesh.numberofvertices,1);
-	disp('      no basalforcings.melting_rate specified: values set as zero');
-end
+%Initialize surface and basal forcings
+md.surfaceforcings = initialize(md.surfaceforcings,md);
+md.basalforcings   = initialize(md.basalforcings,md);
+
+%Deal with other boundary conditions
 if isnan(md.balancethickness.thickening_rate),
 	md.balancethickness.thickening_rate=zeros(md.mesh.numberofvertices,1);
 	disp('      no balancethickness.thickening_rate specified: values set as zero');
 end
-
 md.masstransport.spcthickness=NaN*ones(md.mesh.numberofvertices,1);
 md.balancethickness.spcthickness=NaN*ones(md.mesh.numberofvertices,1);
Index: /issm/trunk-jpl/src/m/boundaryconditions/SetIceSheetBC.py
===================================================================
--- /issm/trunk-jpl/src/m/boundaryconditions/SetIceSheetBC.py	(revision 17070)
+++ /issm/trunk-jpl/src/m/boundaryconditions/SetIceSheetBC.py	(revision 17071)
@@ -34,18 +34,12 @@
 	#No ice front -> do nothing
 
-	#Create zeros basal melting rate and surface mass balance if not specified
-	if numpy.all(numpy.isnan(md.surfaceforcings.precipitation)) and (md.surfaceforcings.ispdd==1):
-		md.surfaceforcings.precipitation=numpy.zeros((md.mesh.numberofvertices,1))
-		print "      no surfaceforcings.precipitation specified: values set as zero"
-	if numpy.all(numpy.isnan(md.surfaceforcings.mass_balance)) and (md.surfaceforcings.ispdd==0):
-		md.surfaceforcings.mass_balance=numpy.zeros((md.mesh.numberofvertices,1))
-		print "      no surfaceforcings.mass_balance specified: values set as zero"
-	if numpy.all(numpy.isnan(md.basalforcings.melting_rate)):
-		md.basalforcings.melting_rate=numpy.zeros((md.mesh.numberofvertices,1))
-		print "      no basalforcings.melting_rate specified: values set as zero"
+	#Create zeros basalforcings and surfaceforcings
+	md.surfaceforcings.initialize(md)
+	md.basalforcings.initialize(md)
+
+	#Deal with other boundary conditions
 	if numpy.all(numpy.isnan(md.balancethickness.thickening_rate)):
 		md.balancethickness.thickening_rate=numpy.zeros((md.mesh.numberofvertices,1))
 		print "      no balancethickness.thickening_rate specified: values set as zero"
-
 	md.masstransport.spcthickness=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
 	md.balancethickness.spcthickness=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
Index: /issm/trunk-jpl/src/m/boundaryconditions/SetIceShelfBC.m
===================================================================
--- /issm/trunk-jpl/src/m/boundaryconditions/SetIceShelfBC.m	(revision 17070)
+++ /issm/trunk-jpl/src/m/boundaryconditions/SetIceShelfBC.m	(revision 17071)
@@ -62,22 +62,13 @@
 end
 
-%Create zeros basalforcings and surfaceforcings
-if (isnan(md.surfaceforcings.precipitation) & (md.surfaceforcings.ispdd==1)),
-	md.surfaceforcings.precipitation=zeros(md.mesh.numberofvertices,1);
-	disp('      no surfaceforcings.precipitation specified: values set as zero');
-end
-if (isnan(md.surfaceforcings.mass_balance) & (md.surfaceforcings.ispdd==0)),
-	        md.surfaceforcings.mass_balance=zeros(md.mesh.numberofvertices,1);
-		disp('      no surfaceforcings.mass_balance specified: values set as zero');
-end
-if isnan(md.basalforcings.melting_rate),
-	md.basalforcings.melting_rate=zeros(md.mesh.numberofvertices,1);
-	disp('      no basalforcings.melting_rate specified: values set as zero');
-end
+%Initialize surface and basal forcings
+md.surfaceforcings = initialize(md.surfaceforcings,md);
+md.basalforcings   = initialize(md.basalforcings,md);
+
+%Deal with other boundary conditions
 if isnan(md.balancethickness.thickening_rate),
 	md.balancethickness.thickening_rate=zeros(md.mesh.numberofvertices,1);
 	disp('      no balancethickness.thickening_rate specified: values set as zero');
 end
-
 md.masstransport.spcthickness=NaN*ones(md.mesh.numberofvertices,1);
 md.balancethickness.spcthickness=NaN*ones(md.mesh.numberofvertices,1);
Index: /issm/trunk-jpl/src/m/boundaryconditions/SetIceShelfBC.py
===================================================================
--- /issm/trunk-jpl/src/m/boundaryconditions/SetIceShelfBC.py	(revision 17070)
+++ /issm/trunk-jpl/src/m/boundaryconditions/SetIceShelfBC.py	(revision 17071)
@@ -77,17 +77,11 @@
 
 	#Create zeros basalforcings and surfaceforcings
-	if numpy.all(numpy.isnan(md.surfaceforcings.precipitation)) and (md.surfaceforcings.ispdd==1):
-		md.surfaceforcings.precipitation=numpy.zeros((md.mesh.numberofvertices,1))
-		print "      no surfaceforcings.precipitation specified: values set as zero"
-	if numpy.all(numpy.isnan(md.surfaceforcings.mass_balance)) and (md.surfaceforcings.ispdd==0):
-		md.surfaceforcings.mass_balance=numpy.zeros((md.mesh.numberofvertices,1))
-		print "      no surfaceforcings.mass_balance specified: values set as zero"
-	if numpy.all(numpy.isnan(md.basalforcings.melting_rate)):
-		md.basalforcings.melting_rate=numpy.zeros((md.mesh.numberofvertices,1))
-		print "      no basalforcings.melting_rate specified: values set as zero"
+	md.surfaceforcings.initialize(md)
+	md.basalforcings.initialize(md)
+
+	#Deal with other boundary conditions
 	if numpy.all(numpy.isnan(md.balancethickness.thickening_rate)):
 		md.balancethickness.thickening_rate=numpy.zeros((md.mesh.numberofvertices,1))
 		print "      no balancethickness.thickening_rate specified: values set as zero"
-
 	md.masstransport.spcthickness=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
 	md.balancethickness.spcthickness=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
Index: /issm/trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.m
===================================================================
--- /issm/trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.m	(revision 17070)
+++ /issm/trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.m	(revision 17071)
@@ -78,17 +78,9 @@
 md.hydrology.spcwatercolumn(pos,1)=1;
 
-%Create zeros basalforcings and surfaceforcings
-if (isnan(md.surfaceforcings.precipitation) & (md.surfaceforcings.ispdd==1)),
-	md.surfaceforcings.precipitation=zeros(md.mesh.numberofvertices,1);
-	disp('      no surfaceforcings.precipitation specified: values set as zero');
-end
-if (isnan(md.surfaceforcings.mass_balance) & (md.surfaceforcings.ispdd==0)),
-	md.surfaceforcings.mass_balance=zeros(md.mesh.numberofvertices,1);
-	disp('      no surfaceforcings.mass_balance specified: values set as zero');
-end
-if isnan(md.basalforcings.melting_rate),
-	md.basalforcings.melting_rate=zeros(md.mesh.numberofvertices,1);
-	disp('      no basalforcings.melting_rate specified: values set as zero');
-end
+%Initialize surface and basal forcings
+md.surfaceforcings = initialize(md.surfaceforcings,md);
+md.basalforcings   = initialize(md.basalforcings,md);
+
+%Deal with other boundary conditions
 if isnan(md.balancethickness.thickening_rate),
 	md.balancethickness.thickening_rate=zeros(md.mesh.numberofvertices,1);
Index: /issm/trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.py
===================================================================
--- /issm/trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.py	(revision 17070)
+++ /issm/trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.py	(revision 17071)
@@ -82,13 +82,8 @@
 
 	#Create zeros basalforcings and surfaceforcings
-	if numpy.all(numpy.isnan(md.surfaceforcings.precipitation)) and (md.surfaceforcings.ispdd==1):
-		md.surfaceforcings.precipitation=numpy.zeros((md.mesh.numberofvertices,1))
-		print "      no surfaceforcings.precipitation specified: values set as zero"
-	if numpy.all(numpy.isnan(md.surfaceforcings.mass_balance)) and (md.surfaceforcings.ispdd==0):
-		md.surfaceforcings.mass_balance=numpy.zeros((md.mesh.numberofvertices,1))
-		print "      no surfaceforcings.mass_balance specified: values set as zero"
-	if numpy.all(numpy.isnan(md.basalforcings.melting_rate)):
-		md.basalforcings.melting_rate=numpy.zeros((md.mesh.numberofvertices,1))
-		print "      no basalforcings.melting_rate specified: values set as zero"
+	md.surfaceforcings.initialize(md)
+	md.basalforcings.initialize(md)
+
+	#Deal with other boundary conditions
 	if numpy.all(numpy.isnan(md.balancethickness.thickening_rate)):
 		md.balancethickness.thickening_rate=numpy.zeros((md.mesh.numberofvertices,1))
Index: /issm/trunk-jpl/src/m/classes/basalforcings.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/basalforcings.m	(revision 17070)
+++ /issm/trunk-jpl/src/m/classes/basalforcings.m	(revision 17071)
@@ -18,4 +18,12 @@
 					error('constructor not supported');
 			end
+		end % }}}
+		function self = initialize(self,md) % {{{
+
+			if isnan(self.melting_rate),
+				self.melting_rate=zeros(md.mesh.numberofvertices,1);
+				disp('      no basalforcings.melting_rate specified: values set as zero');
+			end
+
 		end % }}}
 		function obj = setdefaultparameters(obj) % {{{
Index: /issm/trunk-jpl/src/m/classes/basalforcings.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/basalforcings.py	(revision 17070)
+++ /issm/trunk-jpl/src/m/classes/basalforcings.py	(revision 17071)
@@ -30,4 +30,12 @@
 		return string
 		#}}}
+	def initialize(self,md): # {{{
+
+		if numpy.all(numpy.isnan(self.melting_rate)):
+			self.melting_rate=numpy.zeros((md.mesh.numberofvertices,1))
+			print "      no basalforcings.melting_rate specified: values set as zero"
+
+		return self
+	#}}}
 	def setdefaultparameters(self): # {{{
 		return self
Index: /issm/trunk-jpl/src/m/classes/surfaceforcings.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/surfaceforcings.m	(revision 17070)
+++ /issm/trunk-jpl/src/m/classes/surfaceforcings.m	(revision 17071)
@@ -32,4 +32,16 @@
 					error('constructor not supported');
 			end
+		end % }}}
+		function self = initialize(self,md) % {{{
+
+			if (isnan(self.precipitation) & (self.ispdd==1)),
+				self.precipitation=zeros(md.mesh.numberofvertices,1);
+				disp('      no surfaceforcings.precipitation specified: values set as zero');
+			end
+			if (isnan(self.mass_balance) & (self.ispdd==0)),
+				self.mass_balance=zeros(md.mesh.numberofvertices,1);
+				disp('      no surfaceforcings.mass_balance specified: values set as zero');
+			end
+
 		end % }}}
 		function obj = setdefaultparameters(obj) % {{{
Index: /issm/trunk-jpl/src/m/classes/surfaceforcings.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/surfaceforcings.py	(revision 17070)
+++ /issm/trunk-jpl/src/m/classes/surfaceforcings.py	(revision 17071)
@@ -59,4 +59,15 @@
 		return string
 		#}}}
+	def initialize(self,md): # {{{
+
+		if numpy.all(numpy.isnan(self.precipitation)) and (self.ispdd==1):
+			self.precipitation=numpy.zeros((md.mesh.numberofvertices,1))
+			print "      no surfaceforcings.precipitation specified: values set as zero"
+		if numpy.all(numpy.isnan(self.mass_balance)) and (self.ispdd==0):
+			self.mass_balance=numpy.zeros((md.mesh.numberofvertices,1))
+			print "      no surfaceforcings.mass_balance specified: values set as zero"
+
+		return self
+	#}}}
 	def setdefaultparameters(self): # {{{