Index: /issm/trunk-jpl/src/c/classes/Cfsurfacesquaretransient.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Cfsurfacesquaretransient.cpp	(revision 27720)
+++ /issm/trunk-jpl/src/c/classes/Cfsurfacesquaretransient.cpp	(revision 27721)
@@ -32,7 +32,8 @@
 	this->datatimes      = NULL;
 	this->passedflags    = NULL;
-}
-/*}}}*/
-Cfsurfacesquaretransient::Cfsurfacesquaretransient(char* in_name, int in_definitionenum, int in_model_enum, int in_num_datatimes, IssmDouble* in_datatimes, bool* in_passedflags){/*{{{*/
+	this->J              = 0.;
+}
+/*}}}*/
+Cfsurfacesquaretransient::Cfsurfacesquaretransient(char* in_name, int in_definitionenum, int in_model_enum, int in_num_datatimes, IssmDouble* in_datatimes, bool* in_passedflags, IssmDouble in_J){/*{{{*/
 
 	this->definitionenum=in_definitionenum;
@@ -50,4 +51,6 @@
 	xMemCpy<IssmDouble>(this->datatimes,in_datatimes,this->num_datatimes);
 	xMemCpy<bool>(this->passedflags,in_passedflags,this->num_datatimes);
+
+	this->J = in_J;
 }
 /*}}}*/
@@ -70,4 +73,5 @@
 	/*initialize passedtimes to false*/
 	for(int i=0;i<this->num_datatimes;i++) this->passedflags[i]= false;
+	this->J = 0;
 }
 /*}}}*/
@@ -81,5 +85,5 @@
 /*Object virtual function resolutoin: */
 Object* Cfsurfacesquaretransient::copy() {/*{{{*/
-	Cfsurfacesquaretransient* output = new Cfsurfacesquaretransient(this->name,this->definitionenum, this->model_enum, this->num_datatimes, this->datatimes,this->passedflags);
+	Cfsurfacesquaretransient* output = new Cfsurfacesquaretransient(this->name,this->definitionenum, this->model_enum, this->num_datatimes, this->datatimes,this->passedflags, this->J);
 	return (Object*)output;
 }
@@ -110,4 +114,5 @@
 	marshallhandle->call(this->datatimes,this->num_datatimes);
 	marshallhandle->call(this->passedflags,this->num_datatimes);
+	marshallhandle->call(this->J);
 } 
 /*}}}*/
@@ -130,10 +135,103 @@
 /*}}}*/
 IssmDouble Cfsurfacesquaretransient::Response(FemModel* femmodel){/*{{{*/
-	 _error_("Not implemented yet");
-
-}
-/*}}}*/
+
+	/*recover model time parameters: */
+	IssmDouble time;
+	femmodel->parameters->FindParam(&time,TimeEnum);
+
+	/*Find closest datatime that is less than time*/
+	int pos=-1;
+	for(int i=0;i<this->num_datatimes;i++){
+		if(this->datatimes[i]<=time){
+			pos = i;
+		}
+		else{
+			break;
+		}
+	}
+
+	/*if pos=-1, time is earlier than the first data observation in this dataset*/
+	if(pos==-1){
+		_assert_(this->J==0.);
+		return 0.;
+	}
+
+	/*Check that we have not yet calculated this cost function*/
+	if(this->passedflags[pos]){
+		return 0.;
+		//return this->J; //FIXME
+	}
+
+	/*Calculate cost function for this time slice*/
+	IssmDouble J_part=0.;
+	for(Object* & object : femmodel->elements->objects){
+		Element* element=xDynamicCast<Element*>(object);
+		J_part+=this->Cfsurfacesquaretransient_Calculation(element,model_enum);
+	}
+
+	/*Sum across partition*/
+	IssmDouble J_sum;
+	ISSM_MPI_Allreduce( (void*)&J_part,(void*)&J_sum,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
+	ISSM_MPI_Bcast(&J_sum,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+
+	/*Record this cost function so that we do not recalculate it later*/
+	this->passedflags[pos]= true;
+	this->J += J_sum;
+
+	/*Return full cost function this far*/
+	//return this->J; //FIXME
+	return J_sum;
+}/*}}}*/
 IssmDouble Cfsurfacesquaretransient::Cfsurfacesquaretransient_Calculation(Element* element, int model_enum){/*{{{*/
 
-	_error_("not implemented");
+	IssmDouble Jelem=0.;
+	IssmDouble misfit,Jdet;
+	IssmDouble model,obs,weight;
+	IssmDouble* xyz_list = NULL;
+
+	/*Get basal element*/
+	if(!element->IsOnSurface()) return 0.;
+
+	/*If on water, return 0: */
+	if(!element->IsIceInElement()) return 0.;
+
+	/*Spawn surface element*/
+	Element* topelement = element->SpawnTopElement();
+
+	/* Get node coordinates*/
+	topelement->GetVerticesCoordinates(&xyz_list);
+
+	/*Retrieve all inputs we will be needing: */
+	DatasetInput *datasetinput = topelement->GetDatasetInput(definitionenum); _assert_(datasetinput);
+	Input        *model_input  = topelement->GetInput(model_enum);            _assert_(model_input);
+
+	/* Start  looping on the number of gaussian points: */
+	Gauss* gauss=topelement->NewGauss(2);
+	while(gauss->next()){
+
+		/* Get Jacobian determinant: */
+		topelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
+
+		/*Get all parameters at gaussian point*/
+		datasetinput->GetInputValue(&weight,gauss,WeightsSurfaceObservationEnum);
+		datasetinput->GetInputValue(&obs,gauss,SurfaceObservationEnum);
+		model_input->GetInputValue(&model,gauss);
+
+		/*Compute Misfit
+		 *     
+		 *       1  [           2 ]
+		 *  J = --- | (x - x   )  |
+		 *       2  [       obs   ]
+		 **/
+		misfit=0.5*(model-obs)*(model-obs);
+
+		/*Add to cost function*/
+		Jelem+=misfit*weight*Jdet*gauss->weight;
+	}
+
+	/*clean up and Return: */
+	if(topelement->IsSpawnedElement()){topelement->DeleteMaterials(); delete topelement;};
+	xDelete<IssmDouble>(xyz_list);
+	delete gauss;
+	return Jelem;
 }/*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Cfsurfacesquaretransient.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Cfsurfacesquaretransient.h	(revision 27720)
+++ /issm/trunk-jpl/src/c/classes/Cfsurfacesquaretransient.h	(revision 27721)
@@ -22,9 +22,10 @@
 		IssmDouble *datatimes;
 		bool       *passedflags;
+		IssmDouble  J;
 
 		/*Cfsurfacesquaretransient constructors, destructors :*/
 		Cfsurfacesquaretransient();
 		Cfsurfacesquaretransient(char* in_name, int in_definitionenum, int in_model_enum,int num_datatimes, IssmDouble* in_datatime);
-		Cfsurfacesquaretransient(char* in_name, int in_definitionenum, int in_model_enum,int num_datatimes, IssmDouble* in_datatime, bool* in_timepassedflag);
+		Cfsurfacesquaretransient(char* in_name, int in_definitionenum, int in_model_enum,int num_datatimes, IssmDouble* in_datatime, bool* in_timepassedflag, IssmDouble in_J);
 		~Cfsurfacesquaretransient();