Index: /issm/trunk-jpl/src/c/classes/IoModel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/IoModel.cpp	(revision 27489)
+++ /issm/trunk-jpl/src/c/classes/IoModel.cpp	(revision 27490)
@@ -1124,4 +1124,69 @@
 }
 /*}}}*/
+void  IoModel::FetchData(bool** pmatrix,int* pM,int* pN,const char* data_name){/*{{{*/
+
+	/*output: */
+	int M,N;
+	IssmPDouble* matrix=NULL;
+	bool*        bool_matrix=NULL;
+	int code=0;
+
+	/*recover my_rank:*/
+	int my_rank=IssmComm::GetRank();
+
+	/*Set file pointer to beginning of the data: */
+	fid=this->SetFilePointerToData(&code,NULL,data_name);
+
+	if(code!=5 && code!=6 && code!=7)_error_("expecting a IssmDouble, integer or boolean matrix for \""<<data_name<<"\""<<" (Code is "<<code<<")");
+
+	/*Now fetch: */
+
+	/*We have to read a matrix from disk. First read the dimensions of the matrix, then the whole matrix: */
+	/*numberofelements: */
+	if(my_rank==0){
+		if(fread(&M,sizeof(int),1,fid)!=1) _error_("could not read number of rows for matrix ");
+	}
+
+	ISSM_MPI_Bcast(&M,1,ISSM_MPI_INT,0,IssmComm::GetComm());
+
+	if(my_rank==0){
+		if(fread(&N,sizeof(int),1,fid)!=1) _error_("could not read number of columns for matrix ");
+	}
+	ISSM_MPI_Bcast(&N,1,ISSM_MPI_INT,0,IssmComm::GetComm());
+
+	/*Now allocate matrix: */
+	if(M*N){
+		matrix=xNew<IssmPDouble>(M*N);
+
+		/*Read matrix on node 0, then broadcast: */
+		if(my_rank==0){
+			if(fread(matrix,M*N*sizeof(IssmPDouble),1,fid)!=1) _error_("could not read matrix ");
+		}
+
+		ISSM_MPI_Bcast(matrix,M*N,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm());
+	}
+
+	/*Now cast to bool: */
+	if(M*N){
+		bool_matrix=xNew<bool>(M*N);
+		for(int i=0;i<M;i++){
+			for(int j=0;j<N;j++){
+				bool_matrix[i*N+j]=(bool)matrix[i*N+j];
+			}
+		}
+	}
+	else{
+		bool_matrix=NULL;
+	}
+	/*Free resources:*/
+	xDelete<IssmPDouble>(matrix);
+
+	/*Assign output pointers: */
+	*pmatrix=bool_matrix;
+	if (pM)*pM=M;
+	if (pN)*pN=N;
+
+}
+/*}}}*/
 void  IoModel::FetchData(int** pmatrix,int* pM,int* pN,const char* data_name){/*{{{*/
 	int i,j;
Index: /issm/trunk-jpl/src/c/classes/IoModel.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/IoModel.h	(revision 27489)
+++ /issm/trunk-jpl/src/c/classes/IoModel.h	(revision 27490)
@@ -136,4 +136,5 @@
 		void        FetchData(char*** pstrings,int* pnumstrings,const char* data_name);
 		void        FetchData(int** pmatrix,int* pM,int* pN,const char* data_name);
+		void        FetchData(bool**  pboolmatrix,int* pM,int* pN,const char* data_name);
 		void        FetchData(IssmDouble**  pscalarmatrix,int* pM,int* pN,const char* data_name);
 #if _HAVE_AD_  && !defined(_WRAPPERS_)
Index: /issm/trunk-jpl/src/c/classes/Materials/Matlitho.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Materials/Matlitho.cpp	(revision 27489)
+++ /issm/trunk-jpl/src/c/classes/Materials/Matlitho.cpp	(revision 27490)
@@ -14,25 +14,22 @@
 /*Matlitho constructors and destructor*/
 Matlitho::Matlitho(){/*{{{*/
-	this->numlayers=0;
-	this->radius=NULL;
-	this->viscosity=NULL;
-	this->lame_lambda=NULL;
-	this->lame_mu=NULL;
-	this->burgers_viscosity=NULL;
-	this->burgers_mu=NULL;
-	this->ebm_alpha=NULL;
-	this->ebm_delta=NULL;	
-	this->ebm_taul=NULL;
-	this->ebm_tauh=NULL;
-	this->density=NULL;
-	this->rheologymodel=NULL;
-	this->issolid=NULL;
-	return;
-}
-/*}}}*/
-Matlitho::Matlitho(int matlitho_mid, IoModel* iomodel){/*{{{*/
-
-	IssmDouble* rheologymodeld=NULL;
-	IssmDouble* issolidd=NULL;
+	this->numlayers         = 0;
+	this->radius            = NULL;
+	this->viscosity         = NULL;
+	this->lame_lambda       = NULL;
+	this->lame_mu           = NULL;
+	this->burgers_viscosity = NULL;
+	this->burgers_mu        = NULL;
+	this->ebm_alpha         = NULL;
+	this->ebm_delta         = NULL;
+	this->ebm_taul          = NULL;
+	this->ebm_tauh          = NULL;
+	this->density           = NULL;
+	this->rheologymodel     = NULL;
+	this->issolid           = NULL;
+	return;
+}
+/*}}}*/
+Matlitho::Matlitho(int matlitho_mid, IoModel* iomodel, bool* issolid_in, int* rheo_in){/*{{{*/
 
 	this->mid=matlitho_mid;
@@ -72,23 +69,9 @@
 	xMemCpy<IssmDouble>(this->density, iomodel->Data("md.materials.density"),this->numlayers);
 
-	this->rheologymodel=xNew<IssmDouble>(this->numlayers);
-	xMemCpy<IssmDouble>(this->rheologymodel, iomodel->Data("md.materials.rheologymodel"),this->numlayers);
-
-	this->issolid=xNew<IssmDouble>(this->numlayers);
-	xMemCpy<IssmDouble>(this->issolid, iomodel->Data("md.materials.issolid"),this->numlayers);
-
-	/*rheologymodeld= xNew<IssmDouble>(this->numlayers);
-	this->rheologymodel=xNew<bool>(this->numlayers);
-	xMemCpy<IssmDouble>(rheologymodeld, iomodel->Data("md.materials.rheologymodel"),this->numlayers);
-	for (int i=0;i<this->numlayers;i++)this->rheologymodel[i]=reCast<bool,IssmDouble>(rheologymodeld[i]);
-
-	issolidd= xNew<IssmDouble>(this->numlayers);
+	this->rheologymodel=xNew<int>(this->numlayers);
+	xMemCpy<int>(this->rheologymodel, rheo_in, this->numlayers);
+
 	this->issolid=xNew<bool>(this->numlayers);
-	xMemCpy<IssmDouble>(issolidd, iomodel->Data("md.materials.issolid"),this->numlayers);
-	for (int i=0;i<this->numlayers;i++)this->issolid[i]=reCast<bool,IssmDouble>(issolidd[i]);*/
-
-	/*Free resources: */
-	xDelete<IssmDouble>(rheologymodeld);
-	xDelete<IssmDouble>(issolidd);
+	xMemCpy<bool>(this->issolid, issolid_in, this->numlayers);
 }
 /*}}}*/
@@ -106,6 +89,6 @@
 	xDelete<IssmDouble>(ebm_tauh);
 	xDelete<IssmDouble>(density);
-	xDelete<IssmDouble>(rheologymodel);
-	xDelete<IssmDouble>(issolid);
+	xDelete<int>(rheologymodel);
+	xDelete<bool>(issolid);
 
 	return;
@@ -141,8 +124,6 @@
 		matlitho->ebm_tauh=xNew<IssmDouble>(this->numlayers); xMemCpy<IssmDouble>(matlitho->ebm_tauh, this->ebm_tauh,this->numlayers);
 		matlitho->density=xNew<IssmDouble>(this->numlayers); xMemCpy<IssmDouble>(matlitho->density, this->density,this->numlayers);
-		matlitho->rheologymodel=xNew<IssmDouble>(this->numlayers); xMemCpy<IssmDouble>(matlitho->rheologymodel, this->rheologymodel,this->numlayers);
-		matlitho->issolid=xNew<IssmDouble>(this->numlayers); xMemCpy<IssmDouble>(matlitho->issolid, this->issolid,this->numlayers);
-		/*matlitho->rheologymodel=xNew<bool>(this->numlayers); for(int i=0;i<this->numlayers;i++)matlitho->rheologymodel[i]=this->rheologymodel[i]; 
-		matlitho->issolid=xNew<bool>(this->numlayers); for(int i=0;i<this->numlayers;i++)matlitho->issolid[i]=this->issolid[i];*/ 
+		matlitho->rheologymodel=xNew<int>(this->numlayers); xMemCpy<int>(matlitho->rheologymodel, this->rheologymodel,this->numlayers);
+		matlitho->issolid=xNew<bool>(this->numlayers); xMemCpy<bool>(matlitho->issolid, this->issolid,this->numlayers);
 	}
 
Index: /issm/trunk-jpl/src/c/classes/Materials/Matlitho.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Materials/Matlitho.h	(revision 27489)
+++ /issm/trunk-jpl/src/c/classes/Materials/Matlitho.h	(revision 27490)
@@ -28,9 +28,9 @@
 		IssmDouble*  ebm_tauh;
 		IssmDouble*  density;
-		IssmDouble*  rheologymodel;
-		IssmDouble*  issolid;
+		int*         rheologymodel;
+		bool*        issolid;
 
 		Matlitho();
-		Matlitho(int matlitho_id, IoModel* iomodel);
+		Matlitho(int matlitho_id, IoModel* iomodel, bool* issolid_in, int* rheo_in);
 		~Matlitho();
 		void SetMid(int matlitho_mid);
Index: /issm/trunk-jpl/src/c/cores/love_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/love_core.cpp	(revision 27489)
+++ /issm/trunk-jpl/src/c/cores/love_core.cpp	(revision 27490)
@@ -589,11 +589,11 @@
 	int rheo=matlitho->rheologymodel[layer_index];
 
-	if (vi!=0 && omega!=0.0){ //take into account viscosity in the rigidity if the material isn't a perfect fluid
+	if(vi!=0 && omega!=0.0){ //take into account viscosity in the rigidity if the material isn't a perfect fluid
 		doubletype ka=la0 + 2.0/3.0*mu0; //Bulk modulus
-		if (rheo==2){//EBM
+		if(rheo==2){//EBM
 			mu=muEBM<doubletype>(layer_index, omega, matlitho, femmodel);
 			la=ka-2.0/3.0*mu;
 		} 
-		else if (rheo==1){//Burgers
+		else if(rheo==1){//Burgers
 			doubletype vi2=matlitho->burgers_viscosity[layer_index];
 			doubletype mu2=matlitho->burgers_mu[layer_index];
@@ -674,5 +674,6 @@
 	doubletype frh,frhg0,fgr0,fgr,fn,rm0,rlm,flm;
 	doubletype xmin,xmax,x,dr;
-	doubletype g,ro, issolid;
+	doubletype g,ro;
+	bool       issolid;
 
 	if (pomega) { //frequency and degree dependent terms /*{{{*/
@@ -690,5 +691,5 @@
 			ro=matlitho->density[layer_index];
 			issolid=matlitho->issolid[layer_index];
-			if(issolid!=0){
+			if(issolid){
 				//GetEarthRheology<doubletype>(&la, &mu,layer_index,omega,matlitho);   
 				mu=vars->mu[layer_index*vars->nfreq+vars->ifreq];
@@ -769,5 +770,5 @@
 				g=GetGravity<doubletype>(x*ra,layer_index,femmodel,matlitho,vars);
 
-				if(issolid==1){
+				if(issolid){
 					yi_prefactor[nindex+ 1*6+3]= fn/x;                  // in dy[1*6+3]
 					yi_prefactor[nindex+ 5*6+2]= -(fgr/g0*fn)/x;        // in dy[5*6+2]
@@ -802,5 +803,5 @@
 				g=GetGravity<doubletype>(x*ra,layer_index,femmodel,matlitho,vars);
 				nindex=nsteps*36+n*36;
-				if(issolid==1){
+				if(issolid){
 					yi_prefactor[nindex+ 1*6+5]= -frhg0;       // in dy[1*6+5]
 					yi_prefactor[nindex+ 2*6+0]= -1.0/x;       // in dy[2*6+0]
@@ -825,5 +826,5 @@
 	//computes yi derivatives at r=radius[layer_index]+ n/nstep*(radius[layer_index+1]-radius[layer_index])
 
-	int issolid=matlitho->issolid[layer_index];
+	bool issolid=matlitho->issolid[layer_index];
 	int iy,id,ny, nindex, nstep, nsteps;
 	//femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum);
@@ -842,5 +843,5 @@
 			   fn=(deg*(deg+1.0));
 
-			   if(issolid==1){
+			   if(issolid){
 			   ny = 6;
 
@@ -899,5 +900,5 @@
 	nindex=nsteps*36+n*36;
 
-	if(issolid==1){
+	if(issolid){
 		ny = 6;
 	} else {
@@ -1235,7 +1236,7 @@
 
 		// Boundary Condition matrix - solid regions
-		if (matlitho->issolid[i]){
+		if(matlitho->issolid[i]){
 			one = -1.0;
-			if (i>0) if (!matlitho->issolid[i-1]) one = 1.0;
+			if(i>0) if(!matlitho->issolid[i-1]) one = 1.0;
 			for (int j=0;j<6;j++){
 				yi[(j+6*ici)+ nyi*(j+6*ici+3)] = one;
Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp	(revision 27489)
+++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp	(revision 27490)
@@ -197,7 +197,13 @@
 						break;
 					case MatlithoEnum: { /*{{{*/
-							iomodel->FetchData(13,"md.materials.radius","md.materials.viscosity","md.materials.lame_lambda","md.materials.lame_mu","md.materials.burgers_viscosity","md.materials.burgers_mu","md.materials.ebm_alpha","md.materials.ebm_delta","md.materials.ebm_taul","md.materials.ebm_tauh","md.materials.rheologymodel","md.materials.issolid","md.materials.density");
-							materials->AddObject(new Matlitho(iomodel->numberofelements+1,iomodel));
-							iomodel->DeleteData(13,"md.materials.radius","md.materials.viscosity","md.materials.lame_lambda","md.materials.lame_mu","md.materials.burgers_viscosity","md.materials.burgers_mu","md.materials.ebm_alpha","md.materials.ebm_delta","md.materials.ebm_taul","md.materials.ebm_tauh","md.materials.rheologymodel","md.materials.issolid","md.materials.density");
+							bool* issolid = NULL;
+							int*  rheologymodel = NULL;
+							iomodel->FetchData(&issolid, NULL, NULL, "md.materials.issolid");
+							iomodel->FetchData(&rheologymodel, NULL, NULL, "md.materials.rheologymodel");
+							iomodel->FetchData(11,"md.materials.radius","md.materials.viscosity","md.materials.lame_lambda","md.materials.lame_mu","md.materials.burgers_viscosity","md.materials.burgers_mu","md.materials.ebm_alpha","md.materials.ebm_delta","md.materials.ebm_taul","md.materials.ebm_tauh","md.materials.density");
+							materials->AddObject(new Matlitho(iomodel->numberofelements+1, iomodel, issolid, rheologymodel));
+							iomodel->DeleteData(11,"md.materials.radius","md.materials.viscosity","md.materials.lame_lambda","md.materials.lame_mu","md.materials.burgers_viscosity","md.materials.burgers_mu","md.materials.ebm_alpha","md.materials.ebm_delta","md.materials.ebm_taul","md.materials.ebm_tauh","md.materials.density");
+							xDelete<bool>(issolid);
+							xDelete<int>(rheologymodel);
 						}
 						/*}}}*/