Index: ../trunk-jpl/src/c/analyses/SeaiceAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/SeaiceAnalysis.cpp	(revision 18500)
+++ ../trunk-jpl/src/c/analyses/SeaiceAnalysis.cpp	(revision 18501)
@@ -74,6 +74,8 @@
 
 	/*Intermediaries */
 	IssmDouble  Jdet,D_scalar,dt,thickness,concentration;
+	IssmDouble  ocean_coef,rho_ocean,ocean_lin_drag_coef,ocean_quad_drag_coef;
+	IssmDouble  vx,vy,oceanvx,oceanvy,vnorm,ocean_turning_angle;
 	IssmDouble  C[3][3];
 	IssmDouble* xyz_list = NULL;
 
@@ -89,8 +91,16 @@
 	/*Retrieve all inputs and parameters*/
 	element->GetVerticesCoordinates(&xyz_list);
 	element->FindParam(&dt,TimesteppingTimeStepEnum);
-	IssmDouble rho_i               = element->GetMaterialParameter(MaterialsRhoIceEnum);
-	Input*     thickness_input     = element->GetInput(SeaiceThicknessEnum);     _assert_(thickness_input);
+	element->FindParam(&ocean_coef,BasalforcingsOceanCoefEnum);
+	element->FindParam(&ocean_lin_drag_coef,BasalforcingsOceanLinDragCoefEnum);
+	element->FindParam(&ocean_quad_drag_coef,BasalforcingsOceanQuadDragCoefEnum);
+	element->FindParam(&ocean_turning_angle,BasalforcingsOceanTurningAngleEnum);
+	IssmDouble rho_i           = element->GetMaterialParameter(MaterialsRhoIceEnum);
+	Input*     thickness_input = element->GetInput(SeaiceThicknessEnum);       _assert_(thickness_input);
+	Input*     vx_input        = element->GetInput(VxEnum);                    _assert_(vx_input);
+	Input*     vy_input        = element->GetInput(VyEnum);                    _assert_(vy_input);
+	Input*     oceanvx_input   = element->GetInput(BasalforcingsOceanVxEnum);  _assert_(oceanvx_input);
+	Input*     oceanvy_input   = element->GetInput(BasalforcingsOceanVyEnum);  _assert_(oceanvy_input);
 
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=element->NewGauss(2);
@@ -124,7 +134,17 @@
 					B,3,numdof,0,
 					&Ke->values[0],1);
 
-		/*Create Coriolis part*/
+		/*Create Ocean forcing part*/
+		oceanvx_input->GetInputValue(&oceanvx,gauss);
+		oceanvy_input->GetInputValue(&oceanvy,gauss);
+		vx_input->GetInputValue(&vx,gauss);
+		vy_input->GetInputValue(&vy,gauss);
+		vnorm = sqrt(pow(oceanvx-vx,2) + pow(oceanvy-vy,2));
+		D_scalar = dt*concentration*ocean_coef*rho_ocean*(ocean_lin_drag_coef+ocean_quad_drag_coef*vnorm)*cos(ocean_turning_angle)*gauss->weight*Jdet;
+		TripleMultiply(M,1,numdof,1,
+					&D_scalar,1,1,0,
+					M,1,numdof,0,
+					&Ke->values[0],1);
 	}
 
 	/*Clean up and return*/
@@ -170,7 +190,6 @@
 	element->FindParam(&ocean_coef,BasalforcingsOceanCoefEnum);
 	element->FindParam(&ocean_lin_drag_coef,BasalforcingsOceanLinDragCoefEnum);
 	element->FindParam(&ocean_quad_drag_coef,BasalforcingsOceanQuadDragCoefEnum);
-	element->FindParam(&ocean_quad_drag_coef,BasalforcingsOceanQuadDragCoefEnum);
 	element->FindParam(&ocean_turning_angle,BasalforcingsOceanTurningAngleEnum);
 	rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum);
 	gravity = element->GetMaterialParameter(ConstantsGEnum);
@@ -225,7 +244,7 @@
 
 		/*Ocean forcing*/
 		vnorm = sqrt(pow(oceanvx-vx,2) + pow(oceanvy-vy,2));
-		constant_part = concentration*air_coef*rho_air*(air_lin_drag_coef+air_quad_drag_coef*vnorm);
+		constant_part = concentration*ocean_coef*rho_ocean*(ocean_lin_drag_coef+ocean_quad_drag_coef*vnorm);
 		for(int i=0;i<numnodes;i++){
 			pe->values[i*2+0] += dt*constant_part*(oceanvy-vy)*sin(ocean_turning_angle)*Jdet*gauss->weight*basis[i];
 			pe->values[i*2+1] += dt*constant_part*(vx-oceanvx)*sin(ocean_turning_angle)*Jdet*gauss->weight*basis[i];
