Index: /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp	(revision 28261)
+++ /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp	(revision 28262)
@@ -203,18 +203,22 @@
 			iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.basin_id",BasalforcingsIsmip6BasinIdEnum);
 			iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.melt_anomaly",BasalforcingsIsmip6MeltAnomalyEnum,0.);
-			IssmDouble** array3d = NULL; int* Ms = NULL; int* Ns = NULL; int K;
-			iomodel->FetchData(&array3d,&Ms,&Ns,&K,"md.basalforcings.tf");
-			if(!array3d) _error_("md.basalforcings.tf not found in binary file");
-			for(Object* & object : elements->objects){
-				Element*       element = xDynamicCast<Element*>(object);
-				if(iomodel->domaintype!=Domain2DhorizontalEnum && !element->IsOnBase()) continue;
-				for(int kk=0;kk<K;kk++){
-					element->DatasetInputAdd(BasalforcingsIsmip6TfEnum,array3d[kk],inputs,iomodel,Ms[kk],Ns[kk],1,BasalforcingsIsmip6TfEnum,kk);
+
+			/*Deal with tf...*/
+			IssmDouble* array2d = NULL; int M,N,K; IssmDouble* temp = NULL;
+			iomodel->FetchData(&temp,&M,&K,"md.basalforcings.tf_depths"); xDelete<IssmDouble>(temp);
+			_assert_(M==1); _assert_(K>=1);
+			for(int kk=0;kk<K;kk++){
+
+				/*Fetch TF for this depth*/
+				iomodel->FetchData(&array2d, &M, &N, kk, "md.basalforcings.tf");
+				if(!array2d) _error_("md.basalforcings.tf not found in binary file");
+				for(Object* & object : elements->objects){
+					Element*  element = xDynamicCast<Element*>(object);
+					if(iomodel->domaintype!=Domain2DhorizontalEnum && !element->IsOnBase()) continue;
+					element->DatasetInputAdd(BasalforcingsIsmip6TfEnum,array2d,inputs,iomodel,M,N,1,BasalforcingsIsmip6TfEnum,kk);
 				}
-			}
-			xDelete<int>(Ms); xDelete<int>(Ns);
-			for(int i=0;i<K;i++) xDelete<IssmDouble>(array3d[i]);
-			xDelete<IssmDouble*>(array3d);
-			}
+			xDelete<IssmDouble>(array2d);
+			}
+											  }
 			break;
 		case BeckmannGoosseFloatingMeltRateEnum:
Index: /issm/trunk-jpl/src/c/classes/IoModel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/IoModel.cpp	(revision 28261)
+++ /issm/trunk-jpl/src/c/classes/IoModel.cpp	(revision 28262)
@@ -1565,4 +1565,75 @@
 }
 /*}}}*/
+void  IoModel::FetchData(IssmDouble** pmatrix,int* pM,int* pN, int layer_number,const char* data_name){/*{{{*/
+	/*Same function as above, but here we want to fetch only one specific "layer" of the dataset to avoid memory issues*/
+
+	/*output: */
+	int          M=0,N=0;
+	IssmPDouble *matrix = NULL;
+
+	/*intermediary: */
+	int  numrecords=0;
+	int  code,i;
+
+	/*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!=8)_error_("expecting a IssmDouble mat array for \""<<data_name<<"\"");
+
+	if(my_rank==0){
+		if(fread(&numrecords,sizeof(int),1,fid)!=1) _error_("could not read number of records in matrix array ");
+	}
+	ISSM_MPI_Bcast(&numrecords,1,ISSM_MPI_INT,0,IssmComm::GetComm());
+
+	if(numrecords){
+
+		/*Check consistency of layer number*/
+		_assert_(layer_number>=0);
+		if(layer_number>=numrecords) _error_("layer number of "<<data_name<<" cannot exceed "<< numrecords-1);
+
+		/*Skip until we get to the right layer*/
+		if(my_rank==0){
+			for(i=0;i<layer_number;i++){
+				if(fread(&M,sizeof(int),1,fid)!=1) _error_("could not read number of rows in " << i << "th matrix of matrix array");
+				if(fread(&N,sizeof(int),1,fid)!=1) _error_("could not read number of columns in " << i << "th matrix of matrix array");
+				fseek(fid,M*N*sizeof(double),SEEK_CUR);
+			}
+		}
+
+		/*fetch this slice!*/
+		if(my_rank==0){
+			if(fread(&M,sizeof(int),1,fid)!=1) _error_("could not read number of rows in " << i << "th matrix of matrix array");
+		}
+		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 in " << i << "th matrix of matrix array");
+		}
+		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());
+
+			*pmatrix=xNew<IssmDouble>(M*N);
+			for(int i=0;i<M*N;++i) (*pmatrix)[i]=matrix[i];
+			xDelete<IssmPDouble>(matrix);
+		}
+		else{
+			*pmatrix=NULL;
+		}
+	}
+
+	/*Assign output pointers: */
+	if(pM) *pM = M;
+	if(pN) *pN = N;
+}
+/*}}}*/
 void  IoModel::FetchData(Options* options,const char* lastnonoption){/*{{{*/
 
Index: /issm/trunk-jpl/src/c/classes/IoModel.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/IoModel.h	(revision 28261)
+++ /issm/trunk-jpl/src/c/classes/IoModel.h	(revision 28262)
@@ -144,4 +144,5 @@
 #endif
 		void        FetchData(IssmDouble*** pmatrixarray,int** pmdims,int** pndims, int* pnumrecords,const char* data_name);
+		void        FetchData(IssmDouble** pmatrix,int* pM,int* pN, int layer_number,const char* data_name);
 		void        FetchData(Options *options,const char* data_name);
 		void        FetchData(int num,...);