Index: /issm/trunk-jpl/src/c/analyses/SeaiceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/SeaiceAnalysis.cpp	(revision 18500)
+++ /issm/trunk-jpl/src/c/analyses/SeaiceAnalysis.cpp	(revision 18501)
@@ -75,4 +75,6 @@
 	/*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;
@@ -90,6 +92,14 @@
 	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: */
@@ -125,5 +135,15 @@
 					&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);
 	}
 
@@ -170,5 +190,4 @@
 	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);
@@ -226,5 +245,5 @@
 		/*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];