Index: /issm/trunk-jpl/src/c/Makefile.am
===================================================================
--- /issm/trunk-jpl/src/c/Makefile.am	(revision 23945)
+++ /issm/trunk-jpl/src/c/Makefile.am	(revision 23946)
@@ -117,4 +117,5 @@
 					./classes/Constraints/SpcStatic.cpp\
 					./classes/Constraints/SpcDynamic.cpp\
+					./classes/Loads/Channel.cpp\
 					./classes/Loads/Loads.cpp\
 					./classes/Loads/Penpair.cpp\
Index: /issm/trunk-jpl/src/c/classes/Loads/Channel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Loads/Channel.cpp	(revision 23946)
+++ /issm/trunk-jpl/src/c/classes/Loads/Channel.cpp	(revision 23946)
@@ -0,0 +1,373 @@
+/*!\file Channel.c
+ * \brief: implementation of the Channel object
+ */
+
+/*Headers:*/
+/*{{{*/
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+#include "shared/shared.h"
+#include "../classes.h"
+/*}}}*/	
+#define NUMNODES 2
+
+/*Channel constructors and destructor*/
+Channel::Channel(){/*{{{*/
+	this->parameters = NULL;
+	this->helement   = NULL;
+	this->element    = NULL;
+	this->hnodes     = NULL;
+	this->hvertices  = NULL;
+	this->nodes      = NULL;
+}
+/*}}}*/
+Channel::Channel(int channel_id,int i,int index,IoModel* iomodel){/*{{{*/
+
+	/* Intermediary */
+	int  j;
+	int  pos1,pos2,pos3,pos4;
+	int  num_nodes;
+
+	/*channel constructor data: */
+	int   channel_elem_ids[2];
+	int   channel_vertex_ids[2];
+	int   channel_node_ids[4];
+	int   channel_type;
+
+	/*Get edge*/
+	int i1 = iomodel->faces[4*index+0];
+	int i2 = iomodel->faces[4*index+1];
+	int e1 = iomodel->faces[4*index+2];
+	int e2 = iomodel->faces[4*index+3];
+
+	/*First, see wether this is an internal or boundary edge (if e2=-1)*/
+	if(e2==-1){
+		/* Boundary edge, only one element */
+		num_nodes=2;
+		channel_type=BoundaryEnum;
+		channel_elem_ids[0]=e1;
+	}
+	else{
+		/* internal edge: connected to 2 elements */
+		 num_nodes=4;
+		channel_type=InternalEnum;
+		channel_elem_ids[0]=e1;
+		channel_elem_ids[1]=e2;
+	}
+
+	/*1: Get vertices ids*/
+	channel_vertex_ids[0]=i1;
+	channel_vertex_ids[1]=i2;
+
+	/*2: Get node ids*/
+	if (channel_type==InternalEnum){
+
+		/*Now, we must get the nodes of the 4 nodes located on the edge*/
+
+		/*2: Get the column where these ids are located in the index*/
+		pos1=pos2=pos3=pos4=UNDEF;
+		for(j=0;j<3;j++){
+			if(iomodel->elements[3*(e1-1)+j]==i1) pos1=j+1;
+			if(iomodel->elements[3*(e1-1)+j]==i2) pos2=j+1;
+			if(iomodel->elements[3*(e2-1)+j]==i1) pos3=j+1;
+			if(iomodel->elements[3*(e2-1)+j]==i2) pos4=j+1;
+		}
+		_assert_(pos1!=UNDEF && pos2!=UNDEF && pos3!=UNDEF && pos4!=UNDEF);
+
+		/*3: We have the id of the elements and the position of the vertices in the index
+		 * we can compute their dofs!*/
+		channel_node_ids[0]=3*(e1-1)+pos1;
+		channel_node_ids[1]=3*(e1-1)+pos2;
+		channel_node_ids[2]=3*(e2-1)+pos3;
+		channel_node_ids[3]=3*(e2-1)+pos4;
+	}
+	else{
+
+		/*2: Get the column where these ids are located in the index*/
+		pos1=pos2=UNDEF;
+		for(j=0;j<3;j++){
+			if(iomodel->elements[3*(e1-1)+j]==i1) pos1=j+1;
+			if(iomodel->elements[3*(e1-1)+j]==i2) pos2=j+1;
+		}
+		_assert_(pos1!=UNDEF && pos2!=UNDEF);
+
+		/*3: We have the id of the elements and the position of the vertices in the index
+		 * we can compute their dofs!*/
+		channel_node_ids[0]=3*(e1-1)+pos1;
+		channel_node_ids[1]=3*(e1-1)+pos2;
+	}
+
+	/*Ok, we have everything to build the object: */
+	this->id=channel_id;
+
+	/*Hooks: */
+	this->hnodes    =new Hook(channel_node_ids,num_nodes);
+	this->hvertices =new Hook(&channel_vertex_ids[0],2);
+	this->helement  =new Hook(channel_elem_ids,1); // take only the first element for now
+
+	//this->parameters: we still can't point to it, it may not even exist. Configure will handle this.
+	this->parameters=NULL;
+	this->element=NULL;
+	this->nodes=NULL;
+}
+/*}}}*/
+Channel::~Channel(){/*{{{*/
+	this->parameters=NULL;
+	delete helement;
+	delete hnodes;
+	delete hvertices;
+}
+/*}}}*/
+
+/*Object virtual functions definitions:*/
+Object* Channel::copy() {/*{{{*/
+
+	Channel* channel=NULL;
+
+	channel=new Channel();
+
+	/*copy fields: */
+	channel->id=this->id;
+
+	/*point parameters: */
+	channel->parameters=this->parameters;
+
+	/*now deal with hooks and objects: */
+	channel->hnodes    = (Hook*)this->hnodes->copy();
+	channel->hvertices = (Hook*)this->hvertices->copy();
+	channel->helement  = (Hook*)this->helement->copy();
+
+	/*corresponding fields*/
+	channel->nodes    = (Node**)channel->hnodes->deliverp();
+	channel->vertices = (Vertex**)channel->hvertices->deliverp();
+	channel->element  = (Element*)channel->helement->delivers();
+
+	return channel;
+}
+/*}}}*/
+void    Channel::DeepEcho(void){/*{{{*/
+
+	_printf_("Channel:\n");
+	_printf_("   id: " << id << "\n");
+	hnodes->DeepEcho();
+	hvertices->DeepEcho();
+	helement->DeepEcho();
+	_printf_("   parameters\n");
+	if(parameters)
+	 parameters->DeepEcho();
+	else
+	 _printf_("      NULL\n");
+}		
+/*}}}*/
+void    Channel::Echo(void){/*{{{*/
+	_printf_("Channel:\n");
+	_printf_("   id: " << id << "\n");
+	hnodes->Echo();
+	hvertices->Echo();
+	helement->Echo();
+	_printf_("   parameters: " << parameters << "\n");
+}
+/*}}}*/
+int     Channel::Id(void){/*{{{*/
+	return id;
+}
+/*}}}*/
+void    Channel::Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){ /*{{{*/
+
+	_assert_(this);
+
+	/*ok, marshall operations: */
+	MARSHALLING_ENUM(ChannelEnum);
+	MARSHALLING(id);
+
+	if(marshall_direction==MARSHALLING_BACKWARD){
+		this->hnodes      = new Hook();
+		this->hvertices   = new Hook();
+		this->helement    = new Hook();
+	}
+
+	this->hnodes->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
+	this->helement->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
+	this->hvertices->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
+
+	/*corresponding fields*/
+	nodes    =(Node**)this->hnodes->deliverp();
+	vertices =(Vertex**)this->hvertices->deliverp();
+	element  =(Element*)this->helement->delivers();
+
+}
+/*}}}*/
+int     Channel::ObjectEnum(void){/*{{{*/
+
+	return ChannelEnum;
+
+}
+/*}}}*/
+
+/*Load virtual functions definitions:*/
+void  Channel::Configure(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){/*{{{*/
+
+	/*Take care of hooking up all objects for this element, ie links the objects in the hooks to their respective 
+	 * datasets, using internal ids and offsets hidden in hooks: */
+	hnodes->configure((DataSet*)nodesin);
+	hvertices->configure((DataSet*)verticesin);
+	helement->configure((DataSet*)elementsin);
+
+	/*Initialize hooked fields*/
+	this->nodes    = (Node**)hnodes->deliverp();
+	this->vertices = (Vertex**)hvertices->deliverp();
+	this->element  = (Element*)helement->delivers();
+
+	/*point parameters to real dataset: */
+	this->parameters=parametersin;
+}
+/*}}}*/
+void  Channel::CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs){/*{{{*/
+
+	/*recover some parameters*/
+	ElementMatrix* Ke=NULL;
+	int analysis_type;
+	this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
+
+	if(analysis_type != HydrologyGlaDSAnalysisEnum) _error_("Analysis not supported!!");
+	_error_("STOP");
+
+	/*Add to global matrix*/
+	if(Ke){
+		Ke->AddToGlobal(Kff,Kfs);
+		delete Ke;
+	}
+
+}
+/*}}}*/
+void  Channel::CreatePVector(Vector<IssmDouble>* pf){/*{{{*/
+
+	/*recover some parameters*/
+	ElementVector* pe=NULL;
+	int analysis_type;
+	this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
+	if(analysis_type != HydrologyGlaDSAnalysisEnum) _error_("Analysis not supported!!");
+	_error_("STOP");
+
+	/*Add to global matrix*/
+	if(pe){
+		pe->AddToGlobal(pf);
+		delete pe;
+	}
+
+}
+/*}}}*/
+void  Channel::GetNodesLidList(int* lidlist){/*{{{*/
+
+	_assert_(lidlist);
+	_assert_(nodes);
+
+	for(int i=0;i<NUMNODES;i++) lidlist[i]=nodes[i]->Lid();
+}
+/*}}}*/
+void  Channel::GetNodesSidList(int* sidlist){/*{{{*/
+
+	_assert_(sidlist);
+	_assert_(nodes);
+
+	for(int i=0;i<NUMNODES;i++) sidlist[i]=nodes[i]->Sid();
+}
+/*}}}*/
+int   Channel::GetNumberOfNodes(void){/*{{{*/
+	return NUMNODES;
+}
+/*}}}*/
+bool  Channel::IsPenalty(void){/*{{{*/
+	return false;
+}
+/*}}}*/
+void  Channel::PenaltyCreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs,IssmDouble kmax){/*{{{*/
+
+	/*No stiffness loads applied, do nothing: */
+	return;
+
+}
+/*}}}*/
+void  Channel::PenaltyCreatePVector(Vector<IssmDouble>* pf,IssmDouble kmax){/*{{{*/
+
+	/*No penalty loads applied, do nothing: */
+	return;
+
+}
+/*}}}*/
+void  Channel::ResetHooks(){/*{{{*/
+
+	this->nodes=NULL;
+	this->vertices=NULL;
+	this->element=NULL;
+	this->parameters=NULL;
+
+	/*Get Element type*/
+	this->hnodes->reset();
+	this->hvertices->reset();
+	this->helement->reset();
+
+}
+/*}}}*/
+void  Channel::SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){/*{{{*/
+
+}
+/*}}}*/
+void  Channel::SetwiseNodeConnectivity(int* pd_nz,int* po_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum){/*{{{*/
+
+	/*Output */
+	int d_nz = 0;
+	int o_nz = 0;
+
+	/*Loop over all nodes*/
+	for(int i=0;i<this->GetNumberOfNodes();i++){
+
+		if(!flags[this->nodes[i]->Lid()]){
+
+			/*flag current node so that no other element processes it*/
+			flags[this->nodes[i]->Lid()]=true;
+
+			int counter=0;
+			while(flagsindices[counter]>=0) counter++;
+			flagsindices[counter]=this->nodes[i]->Lid();
+
+			/*if node is clone, we have an off-diagonal non-zero, else it is a diagonal non-zero*/
+			switch(set2_enum){
+				case FsetEnum:
+					if(nodes[i]->fsize){
+						if(this->nodes[i]->IsClone())
+						 o_nz += 1;
+						else
+						 d_nz += 1;
+					}
+					break;
+				case GsetEnum:
+					if(nodes[i]->gsize){
+						if(this->nodes[i]->IsClone())
+						 o_nz += 1;
+						else
+						 d_nz += 1;
+					}
+					break;
+				case SsetEnum:
+					if(nodes[i]->ssize){
+						if(this->nodes[i]->IsClone())
+						 o_nz += 1;
+						else
+						 d_nz += 1;
+					}
+					break;
+				default: _error_("not supported");
+			}
+		}
+	}
+
+	/*Assign output pointers: */
+	*pd_nz=d_nz;
+	*po_nz=o_nz;
+}
+/*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Loads/Channel.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Loads/Channel.h	(revision 23946)
+++ /issm/trunk-jpl/src/c/classes/Loads/Channel.h	(revision 23946)
@@ -0,0 +1,78 @@
+/*!\file Channel.h
+ * \brief: header file for icefront object
+ */
+
+#ifndef _CHANNEL_H_
+#define _CHANNEL_H_
+
+/*Headers:*/
+#include "./Load.h"
+class Hook;
+class Parameters;
+class IoModel;
+class Element;
+class Vertex;
+class ElementMatrix;
+class ElementVector;
+
+class Channel: public Load {
+
+	public: 
+		int id;
+
+		/*Hooks*/
+		Hook *helement;
+		Hook *hnodes;
+		Hook *hvertices;
+
+		/*Corresponding fields*/
+		Element     *element;
+		Vertex     **vertices;
+		Node       **nodes;
+		Parameters  *parameters;
+
+		/*Channel constructors,destructors {{{*/
+		Channel();
+		Channel(int numericalflux_id,int i,int index,IoModel* iomodel);
+		~Channel();
+		/*}}}*/
+		/*Object virtual functions definitions:{{{ */
+		Object *copy();
+		void    DeepEcho();
+		void    Echo();
+		int     Id();
+		void Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction);
+		int     ObjectEnum();
+		/*}}}*/
+		/*Update virtual functions resolution: {{{*/
+		void InputUpdateFromConstant(IssmDouble constant, int name){/*Do nothing*/};
+		void InputUpdateFromConstant(int constant, int name){/*Do nothing*/};
+		void InputUpdateFromConstant(bool constant, int name){_error_("Not implemented yet!");}
+		void InputUpdateFromIoModel(int index, IoModel* iomodel){_error_("not implemented yet");};
+		void InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols, int name, int type){/*Do nothing*/}
+		void InputUpdateFromVector(IssmDouble* vector, int name, int type){/*Do nothing*/}
+		void InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){/*Do nothing*/}
+		/*}}}*/
+		/*Load virtual functions definitions: {{{*/
+		void Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters);
+		void CreateJacobianMatrix(Matrix<IssmDouble>* Jff){_error_("Not implemented yet");};
+		void CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs);
+		void CreatePVector(Vector<IssmDouble>* pf);
+		void GetNodesLidList(int* lidlist);
+		void GetNodesSidList(int* sidlist);
+		int  GetNumberOfNodes(void);
+		bool IsPenalty(void);
+		void PenaltyCreateJacobianMatrix(Matrix<IssmDouble>* Jff,IssmDouble kmax){_error_("Not implemented yet");};
+		void PenaltyCreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* kfs, IssmDouble kmax);
+		void PenaltyCreatePVector(Vector<IssmDouble>* pf, IssmDouble kmax);
+		void ResetHooks();
+		void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters);
+		void SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum);
+		/*}}}*/
+		/*Channel management:{{{*/
+		void           GetNormal(IssmDouble* normal,IssmDouble xyz_list[4][3]);
+		/*}}}*/
+
+};
+
+#endif  /* _NUMERICALFLUX_H_ */
Index: /issm/trunk-jpl/src/c/classes/classes.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/classes.h	(revision 23945)
+++ /issm/trunk-jpl/src/c/classes/classes.h	(revision 23946)
@@ -35,4 +35,5 @@
 
 /*Loads: */
+#include "./Loads/Channel.h"
 #include "./Loads/Loads.h"
 #include "./Loads/Load.h"
Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 23945)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 23946)
@@ -928,4 +928,5 @@
 	CfdragcoeffabsgradEnum,
 	ClosedEnum,
+	ChannelEnum,
 	ColinearEnum,
 	ConstraintsEnum,
Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 23945)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 23946)
@@ -932,4 +932,5 @@
 		case CfdragcoeffabsgradEnum : return "Cfdragcoeffabsgrad";
 		case ClosedEnum : return "Closed";
+		case ChannelEnum : return "Channel";
 		case ColinearEnum : return "Colinear";
 		case ConstraintsEnum : return "Constraints";
Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 23945)
+++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 23946)
@@ -953,4 +953,5 @@
 	      else if (strcmp(name,"Cfdragcoeffabsgrad")==0) return CfdragcoeffabsgradEnum;
 	      else if (strcmp(name,"Closed")==0) return ClosedEnum;
+	      else if (strcmp(name,"Channel")==0) return ChannelEnum;
 	      else if (strcmp(name,"Colinear")==0) return ColinearEnum;
 	      else if (strcmp(name,"Constraints")==0) return ConstraintsEnum;
@@ -997,9 +998,9 @@
 	      else if (strcmp(name,"ExtrapolationAnalysis")==0) return ExtrapolationAnalysisEnum;
 	      else if (strcmp(name,"ExtrudeFromBaseAnalysis")==0) return ExtrudeFromBaseAnalysisEnum;
-	      else if (strcmp(name,"ExtrudeFromTopAnalysis")==0) return ExtrudeFromTopAnalysisEnum;
          else stage=9;
    }
    if(stage==9){
-	      if (strcmp(name,"FemModel")==0) return FemModelEnum;
+	      if (strcmp(name,"ExtrudeFromTopAnalysis")==0) return ExtrudeFromTopAnalysisEnum;
+	      else if (strcmp(name,"FemModel")==0) return FemModelEnum;
 	      else if (strcmp(name,"FileParam")==0) return FileParamEnum;
 	      else if (strcmp(name,"FixedTimestepping")==0) return FixedTimesteppingEnum;
@@ -1120,9 +1121,9 @@
 	      else if (strcmp(name,"MeshX")==0) return MeshXEnum;
 	      else if (strcmp(name,"MeshY")==0) return MeshYEnum;
-	      else if (strcmp(name,"MINIcondensed")==0) return MINIcondensedEnum;
          else stage=10;
    }
    if(stage==10){
-	      if (strcmp(name,"MINI")==0) return MINIEnum;
+	      if (strcmp(name,"MINIcondensed")==0) return MINIcondensedEnum;
+	      else if (strcmp(name,"MINI")==0) return MINIEnum;
 	      else if (strcmp(name,"MinVel")==0) return MinVelEnum;
 	      else if (strcmp(name,"MinVx")==0) return MinVxEnum;
@@ -1243,9 +1244,9 @@
 	      else if (strcmp(name,"TotalSmb")==0) return TotalSmbEnum;
 	      else if (strcmp(name,"TotalSmbScaled")==0) return TotalSmbScaledEnum;
-	      else if (strcmp(name,"TransientArrayParam")==0) return TransientArrayParamEnum;
          else stage=11;
    }
    if(stage==11){
-	      if (strcmp(name,"TransientInput")==0) return TransientInputEnum;
+	      if (strcmp(name,"TransientArrayParam")==0) return TransientArrayParamEnum;
+	      else if (strcmp(name,"TransientInput")==0) return TransientInputEnum;
 	      else if (strcmp(name,"TransientParam")==0) return TransientParamEnum;
 	      else if (strcmp(name,"TransientSolution")==0) return TransientSolutionEnum;