Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 20011)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 20012)
@@ -3618,5 +3618,8 @@
 
 	if(computeelastic | computerigid){
+		int* indices=xNew<int>(gsize);
+		IssmDouble* values=xNew<IssmDouble>(gsize);
 		for(int i=0;i<gsize;i++){
+			indices[i]=i;
 
 			IssmDouble alpha;
@@ -3652,7 +3655,11 @@
 
 			/*Add all components to the pSgi or pSgo solution vectors:*/
-			pSgi->SetValue(i,3*rho_ice/rho_earth*area/eartharea*I*(G_rigid+G_elastic),ADD_VAL);
-			
-		}
+			values[i]=3*rho_ice/rho_earth*area/eartharea*I*(G_rigid+G_elastic);
+		}
+		pSgi->SetValues(gsize,indices,values,ADD_VAL);
+		
+		/*free ressources:*/
+		xDelete<IssmDouble>(values);
+		xDelete<int>(indices);
 	}
 	
@@ -3693,4 +3700,7 @@
 	int         M;
 
+	/*optimization:*/
+	bool store_green_functions=false;
+
 	/*ice properties: */
 	IssmDouble rho_ice,rho_water,rho_earth;
@@ -3716,4 +3726,5 @@
 	this->parameters->FindParam(&computeelastic,SealevelriseElasticEnum);
 	this->parameters->FindParam(&computeeustatic,SealevelriseEustaticEnum);
+	this->parameters->FindParam(&store_green_functions,SealevelriseStoreGreenFunctionsEnum);
 	
 	/*early return if rigid or elastic not requested:*/
@@ -3805,6 +3816,10 @@
 	}
 
+	int* indices=xNew<int>(gsize);
+	IssmDouble* values=xNewZeroInit<IssmDouble>(gsize);
+
 	for(int i=0;i<gsize;i++){
 
+		indices[i]=i; 
 		if(compute_G_elastic | compute_G_rigid){
 			IssmDouble alpha;
@@ -3839,11 +3854,19 @@
 
 		/*Add all components to the pSgi or pSgo solution vectors:*/
-		if(computerigid)pSgo->SetValue(i,3*rho_water/rho_earth*area/eartharea*I*G_rigid[i],ADD_VAL);
-		if(computeelastic)pSgo->SetValue(i,3*rho_water/rho_earth*area/eartharea*I*G_elastic[i],ADD_VAL);
-	}
+		if(computerigid)values[i]+=3*rho_water/rho_earth*area/eartharea*I*G_rigid[i];
+		if(computeelastic)values[i]+=3*rho_water/rho_earth*area/eartharea*I*G_elastic[i];
+	}
+	
+	pSgo->SetValues(gsize,indices,values,ADD_VAL);
+
+	/*free ressources:*/
+	xDelete<IssmDouble>(values);
+	xDelete<int>(indices);
 
 	/*Save G_rigid and G_elastic if computed:*/
-	if(G_rigid)this->inputs->AddInput(new DoubleArrayInput(SealevelriseGRigidEnum,G_rigid,gsize));
-	if(G_elastic)this->inputs->AddInput(new DoubleArrayInput(SealevelriseGElasticEnum,G_elastic,gsize));
+	if(store_green_functions){
+		if(G_rigid)this->inputs->AddInput(new DoubleArrayInput(SealevelriseGRigidEnum,G_rigid,gsize));
+		if(G_elastic)this->inputs->AddInput(new DoubleArrayInput(SealevelriseGElasticEnum,G_elastic,gsize));
+	}
 
 	/*Free ressources:*/
