Index: issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 25766)
+++ issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 25767)
@@ -5446,6 +5446,6 @@
 
 	/*Computational flags:*/
-	bool computerigid = true;
-	bool computeelastic = true;
+	bool computerigid = false;
+	bool computeelastic = false;
 	int  horiz;
 
@@ -5469,4 +5469,8 @@
 	parameter->GetParameterValueByPointer((IssmDouble**)&G_rigid_precomputed,&M);
 
+	/*allocate indices:*/
+	indices=xNew<IssmDouble>(gsize);
+	G=xNewZeroInit<IssmDouble>(gsize);
+
 	if(computeelastic){
 		parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseGElasticEnum)); _assert_(parameter);
@@ -5478,10 +5482,5 @@
 		parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseUElasticEnum)); _assert_(parameter);
 		parameter->GetParameterValueByPointer((IssmDouble**)&U_elastic_precomputed,&M);
-	}
-
-	/*allocate indices:*/
-	indices=xNew<IssmDouble>(gsize);
-	G=xNewZeroInit<IssmDouble>(gsize);
-	if(computeelastic){
+
 		GU=xNewZeroInit<IssmDouble>(gsize);
 		if(horiz){
@@ -5549,6 +5548,6 @@
 
 	/*Add in inputs:*/
-    this->inputs->SetArrayInput(SealevelriseIndicesEnum,this->lid,indices,gsize);
-    this->inputs->SetArrayInput(SealevelriseGEnum,this->lid,G,gsize);
+	this->inputs->SetArrayInput(SealevelriseIndicesEnum,this->lid,indices,gsize);
+	this->inputs->SetArrayInput(SealevelriseGEnum,this->lid,G,gsize);
 	if(computeelastic){
 		this->inputs->SetArrayInput(SealevelriseGUEnum,this->lid,GU,gsize);
@@ -5558,4 +5557,6 @@
 		}
 	}
+	this->inputs->SetDoubleInput(AreaEnum,this->lid,area); 
+	this->AddInput(SealevelAreaEnum,&area,P0Enum);
 
 	/*Free allocations:*/
@@ -5594,5 +5595,5 @@
 	bool notfullygrounded=false;
 	bool scaleoceanarea= false;
-	bool computeelastic= true;
+	bool computerigid= false;
 	int  glfraction=1;
 
@@ -5640,10 +5641,10 @@
 	rho_ice=FindParam(MaterialsRhoIceEnum);
 	rho_water=FindParam(MaterialsRhoSeawaterEnum);
-	this->parameters->FindParam(&computeelastic,SolidearthSettingsElasticEnum);
+	this->parameters->FindParam(&computerigid,SolidearthSettingsRigidEnum);
 	this->parameters->FindParam(&scaleoceanarea,SolidearthSettingsOceanAreaScalingEnum);
 	this->parameters->FindParam(&glfraction,SolidearthSettingsGlfractionEnum);
 
 	/*retrieve precomputed G:*/
-	if(computeelastic)this->inputs->GetArrayPtr(SealevelriseGEnum,this->lid,&G,&gsize);
+	if(computerigid)this->inputs->GetArrayPtr(SealevelriseGEnum,this->lid,&G,&gsize);
 
 	/*Get area of element: precomputed in the sealevelrise_core_geometry:*/
@@ -5699,5 +5700,5 @@
 	_assert_(!xIsNan<IssmDouble>(bslrice));
 
-	if(computeelastic){
+	if(computerigid){
 		/*convert from m to kg/m^2:*/
 		I=I*rho_ice*phi;
@@ -5725,5 +5726,5 @@
 	bool notfullygrounded=false;
 	bool scaleoceanarea= false;
-	bool computeelastic= true;
+	bool computeelastic= false;
 
 	/*elastic green function:*/
@@ -5799,5 +5800,5 @@
 	IssmDouble BP;  //change in bottom pressure (Farrel and Clarke, Equ. 4)
 	IssmDouble constant;
-	bool computeelastic= true;
+	bool computeelastic= false;
 
 	/*elastic green function:*/
@@ -5901,5 +5902,5 @@
 
 	/*computational flags:*/
-	bool computeelastic= true;
+	bool computeelastic= false;
 
 	/*early return if we are not on the ocean or on an ice cap:*/
Index: issm/trunk-jpl/src/c/modules/InputUpdateFromDakotax/InputUpdateFromDakotax.cpp
===================================================================
--- issm/trunk-jpl/src/c/modules/InputUpdateFromDakotax/InputUpdateFromDakotax.cpp	(revision 25766)
+++ issm/trunk-jpl/src/c/modules/InputUpdateFromDakotax/InputUpdateFromDakotax.cpp	(revision 25767)
@@ -12,5 +12,5 @@
 #include "../InputUpdateFromConstantx/InputUpdateFromConstantx.h"
 #include "../InputUpdateFromVectorDakotax/InputUpdateFromVectorDakotax.h"
-			
+
 void  InputUpdateFromDakotax(FemModel* femmodel,double* variables,char* *variables_descriptors,int numdakotavariables){ /*{{{*/
 
@@ -34,27 +34,27 @@
 
 	/*retrieve parameters: */
-	femmodel->parameters->FindParam(&variable_partitions,&variable_partitions_num,NULL,NULL,QmuVariablePartitionsEnum); 
-	femmodel->parameters->FindParam(&variable_partitions_npart,NULL,NULL,QmuVariablePartitionsNpartEnum); 
-	femmodel->parameters->FindParam(&variable_partitions_nt,NULL,NULL,QmuVariablePartitionsNtEnum); 
-	
-
-	/*Go through all dakota descriptors, ex: "rho_ice","thermal_conductivity","thickness1","thickness2", etc ..., and 
-	 * for each descriptor, take the variable value and plug it into the inputs (more or less :)): 
-	 * We also start with distributed and standard values , as they tend to be used to pluck data from a multi-modle ensemble (mme) 
-	 * which can then be scaled. Doing the scaling first would be impractical, as the entire mme would have to be scaled, 
+	femmodel->parameters->FindParam(&variable_partitions,&variable_partitions_num,NULL,NULL,QmuVariablePartitionsEnum);
+	femmodel->parameters->FindParam(&variable_partitions_npart,NULL,NULL,QmuVariablePartitionsNpartEnum);
+	femmodel->parameters->FindParam(&variable_partitions_nt,NULL,NULL,QmuVariablePartitionsNtEnum);
+
+
+	/*Go through all dakota descriptors, ex: "rho_ice","thermal_conductivity","thickness1","thickness2", etc ..., and
+	 * for each descriptor, take the variable value and plug it into the inputs (more or less :)):
+	 * We also start with distributed and standard values , as they tend to be used to pluck data from a multi-modle ensemble (mme)
+	 * which can then be scaled. Doing the scaling first would be impractical, as the entire mme would have to be scaled,
 	 * which is a waste of time:*/
 
 	variablecount=0;
-	for(i=0;i<numdakotavariables;i++){ //these are the dakota variables, for all partitions. 
+	for(i=0;i<numdakotavariables;i++){ //these are the dakota variables, for all partitions.
 
 		descriptor=variables_descriptors[i];
-	
+
 
 		/*From descriptor, figure out if the variable is scaled, indexed, distributed or just a simple variable: */
-		if (strncmp(descriptor,"scaled_",7)==0){ 
+		if (strncmp(descriptor,"scaled_",7)==0){
 			/*we are skipping these for now.*/
 			npart=variable_partitions_npart[variablecount];
 			nt=variable_partitions_nt[variablecount];
-				
+
 			/*increment i to skip the distributed values just collected: */
 			i+=npart*nt-1; //careful, the for loop will add 1.
@@ -66,13 +66,13 @@
 			/*we are skipping these for now.*/
 		}
-		
+
 		else if (strncmp(descriptor,"distributed_",12)==0){
 			if (VerboseQmu())_printf0_("   updating variable " << descriptor << "\n");
-			
+
 			/*recover partition vector: */
 			variable_partition=variable_partitions[variablecount];
 			npart=variable_partitions_npart[variablecount];
 
-			/*Variable is distributed. Determine root name of variable (ex: distributed_DragCoefficient_1 -> DragCoefficient). 
+			/*Variable is distributed. Determine root name of variable (ex: distributed_DragCoefficient_1 -> DragCoefficient).
 			 * Allocate distributed_values and fill the distributed_values with the next npart variables: */
 
@@ -86,11 +86,11 @@
 
 			//for (int j=0;j<npart;j++)_printf_(j << ":" << distributed_values[j] << "\n");
-			
+
 			//Call specialty code:
 			InputUpdateSpecialtyCode(femmodel,distributed_values,variable_partition,npart,root);
-			
+
 			/*increment i to skip the distributed values just collected: */
 			i+=npart-1; //careful, the for loop will add 1.
-			
+
 			/*Free allocations: */
 			xDelete<double>(parameter);
@@ -104,16 +104,16 @@
 		variablecount++;
 	}
-	
+
 	variablecount=0;
 	/*now deal with scaled variabes:*/
-	for(i=0;i<numdakotavariables;i++){ //these are the dakota variables, for all partitions. 
+	for(i=0;i<numdakotavariables;i++){ //these are the dakota variables, for all partitions.
 
 		descriptor=variables_descriptors[i];
-	
+
 		/*From descriptor, figure out if the variable is scaled, indexed, distributed or just a simple variable: */
 		if (strncmp(descriptor,"scaled_",7)==0){
-		
+
 			if (VerboseQmu())_printf0_("   updating variable " << descriptor << "\n");
-		
+
 			/*recover partition vector: */
 			variable_partition=variable_partitions[variablecount];
@@ -121,5 +121,5 @@
 			nt=variable_partitions_nt[variablecount];
 
-			/* Variable is scaled, determine its root name (ex: scaled_DragCoefficient_1 -> DragCoefficient). Allocate distributed_values and fill the 
+			/* Variable is scaled, determine its root name (ex: scaled_DragCoefficient_1 -> DragCoefficient). Allocate distributed_values and fill the
 			 * distributed_values with the next npart variables coming from Dakota: */
 			memcpy(root,strstr(descriptor,"_")+1,(strlen(strstr(descriptor,"_")+1)+1)*sizeof(char));
@@ -136,5 +136,5 @@
 			/*increment i to skip the distributed values just collected: */
 			i+=npart*nt-1; //careful, the for loop will add 1.
-			
+
 			/*Free allocations: */
 			xDelete<double>(parameter);
@@ -202,44 +202,33 @@
 	transientinput2 = femmodel->inputs->GetTransientInput(DummyEnum);
 
-	for (int p=npart;p>=0;p--){
-		int pp=p; 
-		if (p==npart)pp=-1; /*so, the logic is, we want to do the -1 partition first, then 
-							 go from npart-1 to 0 in reverse order:*/
-		
-		for (int i=0;i<femmodel->elements->Size();i++){
-			int element_partition;
-
-			Tria*   element=xDynamicCast<Tria*>(femmodel->elements->GetObjectByOffset(i));
-			
-			element_partition= (int)variable_partition[element->Sid()];
-			if(element_partition!=pp)continue;
-
-			if(element_partition==-1)id=0; //grab background field
-			else id=distributed_values[element_partition]-1; //grab partition field
-
-			/*recover the right field from the mme: */
-			transientinput = datasetinput->GetTransientInputByOffset(id); _assert_(transientinput);
-
-			/*copy values from the transientinput to the final transientinput2: */
-			for (int j=0;j<N;j++){
-				TriaInput* tria_input=transientinput->GetTriaInput(j);
-				element->InputServe(tria_input);
-				if(interpolationenum==P0Enum){
-					value=tria_input->element_values[0];
-					transientinput2->AddTriaTimeInput( j,1,&(element->lid),&value,P0Enum); 
-				}
-				else if(interpolationenum==P1Enum){
-
-					/*Get values and lid list*/
-					const int   numvertices     = element->GetNumberOfVertices();
-					int        *vertexlids      = xNew<int>(numvertices);
-					int        *vertexsids      = xNew<int>(numvertices);
-
-					/*Recover vertices ids needed to initialize inputs*/
-					element->GetVerticesLidList(&vertexlids[0]);
-					element->GetVerticesSidList(&vertexsids[0]);
-					values=tria_input->element_values;
-					transientinput2->AddTriaTimeInput( j,numvertices,vertexlids,values,P1Enum); 
-				}
+	for(Object* & object : femmodel->elements->objects){
+		Tria*   element=xDynamicCast<Tria*>(object);
+
+		if((int)variable_partition[element->Sid()]==-1)id=0; //grab background field
+		else id=distributed_values[(int)variable_partition[element->Sid()]]-1; //grab partition field
+
+		/*recover the right field from the mme: */
+		transientinput = datasetinput->GetTransientInputByOffset(id); _assert_(transientinput);
+
+		/*copy values from the transientinput to the final transientinput2: */
+		for (int j=0;j<N;j++){
+			TriaInput* tria_input=transientinput->GetTriaInput(j);
+			element->InputServe(tria_input);
+			if(interpolationenum==P0Enum){
+				value=tria_input->element_values[0];
+				transientinput2->AddTriaTimeInput( j,1,&(element->lid),&value,P0Enum);
+			}
+			else if(interpolationenum==P1Enum){
+
+				/*Get values and lid list*/
+				const int   numvertices     = element->GetNumberOfVertices();
+				int        *vertexlids      = xNew<int>(numvertices);
+				int        *vertexsids      = xNew<int>(numvertices);
+
+				/*Recover vertices ids needed to initialize inputs*/
+				element->GetVerticesLidList(&vertexlids[0]);
+				element->GetVerticesSidList(&vertexsids[0]);
+				values=tria_input->element_values;
+				transientinput2->AddTriaTimeInput( j,numvertices,vertexlids,values,P1Enum);
 			}
 		}
@@ -258,10 +247,10 @@
 
 	/*Go through elements, copy input name to dummy, and scale it using the distributed_values and the partition vector:*/
-	for(int i=0;i<femmodel->elements->Size();i++){
-		Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
+	for(Object* & object : femmodel->elements->objects){
+		Element* element=xDynamicCast<Element*>(object);
 		element->InputScaleFromDakota(distributed_values,partition,npart,nt,name);
 	}
 
-	/*We created a dummy input, which was a scaled copy of the name input. Now wipe 
+	/*We created a dummy input, which was a scaled copy of the name input. Now wipe
 	 * out the name input with the new input:*/
 	femmodel->inputs->ChangeEnum(DummyEnum,name);
Index: issm/trunk-jpl/src/m/classes/model.m
===================================================================
--- issm/trunk-jpl/src/m/classes/model.m	(revision 25766)
+++ issm/trunk-jpl/src/m/classes/model.m	(revision 25767)
@@ -195,12 +195,59 @@
 			switch nargin
 				case 0
-					md=setdefaultparameters(md);
-				case 1
-					error('model constructor not supported yet');
-
+					md=setdefaultparameters(md,'earth');
 				otherwise
-					error('model constructor error message: 0 of 1 argument only in input.');
-				end
-
+					options=pairoptions(varargin{:});
+					planet=getfieldvalue(options,'planet','earth');
+					md=setdefaultparameters(md,planet);
+				end
+
+		end
+		%}}}
+		function md = setdefaultparameters(md,planet) % {{{
+
+			%initialize subclasses
+			md.mesh             = mesh2d();
+			md.mask             = mask();
+			md.constants        = constants();
+			md.geometry         = geometry();
+			md.initialization   = initialization();
+			md.smb              = SMBforcing();
+			md.basalforcings    = basalforcings();
+			md.friction         = friction();
+			md.rifts            = rifts();
+			md.solidearth       = solidearth(planet);
+			md.dsl              = dsl();
+			md.timestepping     = timestepping();
+			md.groundingline    = groundingline();
+			md.materials        = matice();
+			md.damage           = damage();
+			md.flowequation     = flowequation();
+			md.debug            = debug();
+			md.verbose          = verbose();
+			md.settings         = issmsettings();
+			md.toolkits         = toolkits();
+			md.cluster          = generic();
+			md.balancethickness = balancethickness();
+			md.stressbalance    = stressbalance();
+			md.hydrology        = hydrologyshreve();
+			md.masstransport    = masstransport();
+			md.thermal          = thermal();
+			md.steadystate      = steadystate();
+			md.transient        = transient();
+			md.levelset         = levelset();
+			md.calving          = calving();
+			md.frontalforcings  = frontalforcings();
+			md.gia              = giamme();
+			md.love             = fourierlove();
+			md.esa              = esa();
+			md.autodiff         = autodiff();
+			md.inversion        = inversion();
+			md.qmu              = qmu();
+			md.amr              = amr();
+			md.radaroverlay     = radaroverlay();
+			md.results          = struct();
+			md.outputdefinition = outputdefinition();
+			md.miscellaneous    = miscellaneous();
+			md.private          = private();
 		end
 		%}}}
@@ -1288,52 +1335,4 @@
 			end
 		end% }}}
-		function md = setdefaultparameters(md) % {{{
-
-			%initialize subclasses
-			md.mesh             = mesh2d();
-			md.mask             = mask();
-			md.constants        = constants();
-			md.geometry         = geometry();
-			md.initialization   = initialization();
-			md.smb              = SMBforcing();
-			md.basalforcings    = basalforcings();
-			md.friction         = friction();
-			md.rifts            = rifts();
-			md.solidearth       = solidearth();
-			md.dsl              = dsl();
-			md.timestepping     = timestepping();
-			md.groundingline    = groundingline();
-			md.materials        = matice();
-			md.damage           = damage();
-			md.flowequation     = flowequation();
-			md.debug            = debug();
-			md.verbose          = verbose();
-			md.settings         = issmsettings();
-			md.toolkits         = toolkits();
-			md.cluster          = generic();
-			md.balancethickness = balancethickness();
-			md.stressbalance    = stressbalance();
-			md.hydrology        = hydrologyshreve();
-			md.masstransport    = masstransport();
-			md.thermal          = thermal();
-			md.steadystate      = steadystate();
-			md.transient        = transient();
-			md.levelset         = levelset();
-			md.calving          = calving();
-			md.frontalforcings  = frontalforcings();
-			md.gia              = giamme();
-			md.love             = fourierlove();
-			md.esa              = esa();
-			md.autodiff         = autodiff();
-			md.inversion        = inversion();
-			md.qmu              = qmu();
-			md.amr              = amr();
-			md.radaroverlay     = radaroverlay();
-			md.results          = struct();
-			md.outputdefinition = outputdefinition();
-			md.miscellaneous    = miscellaneous();
-			md.private          = private();
-		end
-		%}}}
 		function md = tetras(md,varargin) % {{{
 			%TETRAS - split 3d prismatic mesh into 3 tetrahedrons
Index: issm/trunk-jpl/src/m/classes/model.py
===================================================================
--- issm/trunk-jpl/src/m/classes/model.py	(revision 25766)
+++ issm/trunk-jpl/src/m/classes/model.py	(revision 25767)
@@ -3,4 +3,5 @@
 import copy
 import sys
+from pairoptions import pairoptions
 from mesh2d import mesh2d
 from mesh3dprisms import mesh3dprisms
@@ -79,49 +80,59 @@
 class model(object):
     #properties
-    def __init__(self): #{{{
-        self.mesh = mesh2d()
-        self.mask = mask()
-        self.constants = constants()
-        self.geometry = geometry()
-        self.initialization = initialization()
-        self.smb = SMBforcing()
-        self.basalforcings = basalforcings()
-        self.friction = friction()
-        self.rifts = rifts()
-        self.solidearth = solidearth()
-        self.dsl = dsl()
-        self.timestepping = timestepping()
-        self.groundingline = groundingline()
-        self.materials = matice()
-        self.damage = damage()
-        self.flowequation = flowequation()
-        self.debug = debug()
-        self.verbose = verbose()
-        self.settings = issmsettings()
-        self.toolkits = toolkits()
-        self.cluster = generic()
-        self.balancethickness = balancethickness()
-        self.stressbalance = stressbalance()
-        self.hydrology = hydrologyshreve()
-        self.masstransport = masstransport()
-        self.thermal = thermal()
-        self.steadystate = steadystate()
-        self.transient = transient()
-        self.levelset = levelset()
-        self.calving = calving()
-        self.frontalforcings = frontalforcings()
-        self.gia = giamme()
-        self.love = fourierlove()
-        self.esa = esa()
-        self.autodiff = autodiff()
-        self.inversion = inversion()
-        self.qmu = qmu()
-        self.amr = amr()
-        self.radaroverlay = radaroverlay()
-        self.results = results()
-        self.outputdefinition = outputdefinition()
-        self.miscellaneous = miscellaneous()
-        self.private = private()
-    #}}}
+    def __init__(self, *args): #{{{
+        self.mesh = None
+        self.mask = None
+        self.constants = None
+        self.geometry = None
+        self.initialization = None
+        self.smb = None
+        self.basalforcings = None
+        self.friction = None
+        self.rifts = None
+        self.solidearth = None
+        self.dsl = None
+        self.timestepping = None
+        self.groundingline = None
+        self.materials = None
+        self.damage = None
+        self.flowequation = None
+        self.debug = None
+        self.verbose = None
+        self.settings = None
+        self.toolkits = None
+        self.cluster = None
+        self.balancethickness = None
+        self.stressbalance = None
+        self.hydrology = None
+        self.masstransport = None
+        self.thermal = None
+        self.steadystate = None
+        self.transient = None
+        self.levelset = None
+        self.calving = None
+        self.frontalforcings = None
+        self.gia = None
+        self.love = None
+        self.esa = None
+        self.autodiff = None
+        self.inversion = None
+        self.qmu = None
+        self.amr = None
+        self.radaroverlay = None
+        self.results = None
+        self.outputdefinition = None
+        self.miscellaneous = None
+        self.private = None
+
+        nargs = len(args)
+
+        if nargs == 0:
+            self.setdefaultparameters('earth')
+        else:
+            self.setdefaultparameters(args[0])
+            options = pairoptions(*args)
+            planet = options.getfieldvalue('planet', 'earth')
+            self.setdefaultparameters(planet)
+#}}}
 
     def properties(self):  # {{{
@@ -224,4 +235,50 @@
     # }}}
 
+    def setdefaultparameters(self, planet): #{{{
+        self.mesh = mesh2d()
+        self.mask = mask()
+        self.constants = constants()
+        self.geometry = geometry()
+        self.initialization = initialization()
+        self.smb = SMBforcing()
+        self.basalforcings = basalforcings()
+        self.friction = friction()
+        self.rifts = rifts()
+        self.solidearth = solidearth(planet)
+        self.dsl = dsl()
+        self.timestepping = timestepping()
+        self.groundingline = groundingline()
+        self.materials = matice()
+        self.damage = damage()
+        self.flowequation = flowequation()
+        self.debug = debug()
+        self.verbose = verbose()
+        self.settings = issmsettings()
+        self.toolkits = toolkits()
+        self.cluster = generic()
+        self.balancethickness = balancethickness()
+        self.stressbalance = stressbalance()
+        self.hydrology = hydrologyshreve()
+        self.masstransport = masstransport()
+        self.thermal = thermal()
+        self.steadystate = steadystate()
+        self.transient = transient()
+        self.levelset = levelset()
+        self.calving = calving()
+        self.frontalforcings = frontalforcings()
+        self.gia = giamme()
+        self.love = fourierlove()
+        self.esa = esa()
+        self.autodiff = autodiff()
+        self.inversion = inversion()
+        self.qmu = qmu()
+        self.amr = amr()
+        self.radaroverlay = radaroverlay()
+        self.results = results()
+        self.outputdefinition = outputdefinition()
+        self.miscellaneous = miscellaneous()
+        self.private = private()
+    #}}}
+
     def checkmessage(self, string):  # {{{
         print("model not consistent: {}".format(string))
Index: issm/trunk-jpl/src/m/classes/solidearth.m
===================================================================
--- issm/trunk-jpl/src/m/classes/solidearth.m	(revision 25766)
+++ issm/trunk-jpl/src/m/classes/solidearth.m	(revision 25767)
@@ -6,5 +6,5 @@
 classdef solidearth
 	properties (SetAccess=public) 
-		initialsealevel               = NaN;
+		initialsealevel        = NaN;
 		settings               = solidearthsettings(); 
 		external               = [];
@@ -22,26 +22,28 @@
 			switch nargin
 				case 0
-					self=setdefaultparameters(self);
+					self=setdefaultparameters(self,'earth');
+				case 1
+					self=setdefaultparameters(self,varargin{:});
 				otherwise
-					error('constructor not supported');
+					error('solidearth constructor error message: zero or one argument only!'); 
 			end
 		end % }}}
-		function self = setdefaultparameters(self) % {{{
-		
-		%output default:
-		self.requested_outputs={'default'};
+		function self = setdefaultparameters(self,planet) % {{{
 
-		%transitions should be a cell array of vectors: 
-		self.transitions={};
-		
-		%no partitions requested for barystatic contribution: 
-		self.partitionice=[];
-		self.partitionhydro=[];
+			%output default:
+			self.requested_outputs={'default'};
 
-		%no external solutions by defalt: 
-		self.external=[];
+			%transitions should be a cell array of vectors: 
+			self.transitions={};
+			
+			%no partitions requested for barystatic contribution: 
+			self.partitionice=[];
+			self.partitionhydro=[];
 
-		%earth radius
-		self.planetradius= planetradius('earth');
+			%no external solutions by defalt: 
+			self.external=[];
+
+			%earth radius
+			self.planetradius= planetradius(planet);
 		
 		end % }}}
Index: issm/trunk-jpl/src/m/classes/solidearth.py
===================================================================
--- issm/trunk-jpl/src/m/classes/solidearth.py	(revision 25766)
+++ issm/trunk-jpl/src/m/classes/solidearth.py	(revision 25767)
@@ -32,10 +32,12 @@
         self.partitionhydro     = []
 
-        nargin = len(args)
+        nargs = len(args)
 
-        if nargin == 0:
-            self.setdefaultparameters()
+        if nargs == 0:
+            self.setdefaultparameters('earth')
+        elif nargs == 1:
+            self.setdefaultparameters(args[0])
         else:
-            raise Exception('constructor not supported')
+            raise Exception('solidearth constructor error message: zero or one argument only!')
     #}}}
 
@@ -59,5 +61,5 @@
     #}}}
 
-    def setdefaultparameters(self):  # {{{
+    def setdefaultparameters(self, planet):  # {{{
         # Default output
         self.requested_outputs = ['default']
@@ -74,5 +76,5 @@
 
         # Earth radius
-        self.planetradius = planetradius('earth')
+        self.planetradius = planetradius(planet)
     #}}}
 
Index: issm/trunk-jpl/src/m/classes/solidearthsettings.m
===================================================================
--- issm/trunk-jpl/src/m/classes/solidearthsettings.m	(revision 25766)
+++ issm/trunk-jpl/src/m/classes/solidearthsettings.m	(revision 25767)
@@ -71,4 +71,9 @@
 			md = checkfield(md,'fieldname','solidearth.settings.glfraction','values',[0 1]);
 
+			%checks on computational flags
+			if self.elastic==1 & self.rigid==0,
+				error('solidearthsettings checkconsistency error message: need rigid on if elastic flag is set');
+			end
+
 			%a coupler to a planet model is provided. 
 			if self.isgrd,
Index: issm/trunk-jpl/src/m/classes/solidearthsettings.py
===================================================================
--- issm/trunk-jpl/src/m/classes/solidearthsettings.py	(revision 25766)
+++ issm/trunk-jpl/src/m/classes/solidearthsettings.py	(revision 25767)
@@ -92,4 +92,8 @@
         md = checkfield(md, 'fieldname', 'solidearth.settings.glfraction', 'values', [0, 1])
 
+        # Checks on computational flags
+        if self.elastic and not self.rigid:
+            raise Exception('solidearthsettings checkconsistency error message: need rigid on if elastic flag is set')
+
         # A coupler to planet model is provided
         if self.isgrd:
Index: issm/trunk-jpl/src/m/geometry/planetradius.py
===================================================================
--- issm/trunk-jpl/src/m/geometry/planetradius.py	(revision 25766)
+++ issm/trunk-jpl/src/m/geometry/planetradius.py	(revision 25767)
@@ -11,5 +11,7 @@
 
     if planet == 'earth':
-        radius = 6.371012 * pow(10, 6)
+        radius = 6.371012e6
+    elif planet == 'europa':
+        radius = 1.5008e6
     else:
         raise TypeError("planet type %s not supported yet!" % planet)
