Index: ../trunk-jpl/src/c/analyses/SeaiceAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/SeaiceAnalysis.cpp	(revision 18495)
+++ ../trunk-jpl/src/c/analyses/SeaiceAnalysis.cpp	(revision 18496)
@@ -38,6 +38,7 @@
 	iomodel->FetchDataToInput(elements,VyStarEnum);
 	iomodel->FetchDataToInput(elements,VxEnum);
 	iomodel->FetchDataToInput(elements,VyEnum);
+	iomodel->FetchDataToInput(elements,SeaiceCoriolisFactorEnum);
 }/*}}}*/
 void SeaiceAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
 
@@ -72,7 +73,7 @@
 	if(concentration_input->Max()==0.) return NULL;
 
 	/*Intermediaries */
-	IssmDouble  Jdet,D_scalar,dt;
+	IssmDouble  Jdet,D_scalar,dt,thickness,concentration;
 	IssmDouble  C[3][3];
 	IssmDouble* xyz_list = NULL;
 
@@ -88,17 +89,23 @@
 	/*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);
 
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=element->NewGauss(2);
 	for(int ig=gauss->begin();ig<gauss->end();ig++){
 
 		gauss->GaussPoint(ig);
+
+		/*Get Jacobian Determinant and thickness*/
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+		thickness_input->GetInputValue(&thickness,gauss);
+		concentration_input->GetInputValue(&concentration,gauss);
 
 		/*Create first part of the stiffness matrix*/
 		this->GetM(M,element,xyz_list,gauss);
-		D_scalar=gauss->weight*Jdet;
+		D_scalar=rho_i*thickness*gauss->weight*Jdet;
 		TripleMultiply(M,1,numdof,1,
 					&D_scalar,1,1,0,
 					M,1,numdof,0,
@@ -116,6 +123,8 @@
 					&C[0][0],3,3,0,
 					B,3,numdof,0,
 					&Ke->values[0],1);
+
+		/*Create Coriolis part*/
 	}
 
 	/*Clean up and return*/
@@ -132,9 +141,13 @@
 	if(concentration_input->Max()==0.) return NULL;
 
 	/*Intermediaries */
-	IssmDouble  concentration,thickness,dt,Jdet,vx,vy;
-	IssmDouble  rho_air,air_coef,air_lin_drag_coef,air_quad_drag_coef;
-	IssmDouble  windvx,windvy,wind_vnorm;
+	IssmDouble  air_coef,ocean_coef,constant_part;
+	IssmDouble  rho_ice,rho_air,rho_ocean,gravity;
+	IssmDouble  vx,vy,vxstar,vystar,windvx,windvy,oceanvx,oceanvy,vnorm;
+	IssmDouble  air_lin_drag_coef,air_quad_drag_coef;
+	IssmDouble  ocean_lin_drag_coef,ocean_quad_drag_coef;
+	IssmDouble  concentration,thickness,coriolis_factor,dt,Jdet;
+	IssmDouble  ocean_turning_angle,dssh[2];
 	IssmDouble* xyz_list = NULL;
 
 	/*Fetch number of nodes and dof for this finite element*/
@@ -152,11 +165,25 @@
 	element->FindParam(&air_coef,SurfaceforcingsAirCoefEnum);
 	element->FindParam(&air_lin_drag_coef,SurfaceforcingsAirLinDragCoefEnum);
 	element->FindParam(&air_quad_drag_coef,SurfaceforcingsAirQuadDragCoefEnum);
+	element->FindParam(&rho_ocean,BasalforcingsRhoOceanEnum);
+	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);
+	Input* thickness_input     = element->GetInput(SeaiceThicknessEnum);       _assert_(thickness_input);
+	Input* coriolis_fact_input = element->GetInput(SeaiceCoriolisFactorEnum);  _assert_(coriolis_fact_input);
 	Input* vx_input            = element->GetInput(VxEnum);                    _assert_(vx_input);
 	Input* vy_input            = element->GetInput(VyEnum);                    _assert_(vy_input);
+	Input* vxstar_input        = element->GetInput(VxStarEnum);                _assert_(vxstar_input);
+	Input* vystar_input        = element->GetInput(VyStarEnum);                _assert_(vystar_input);
 	Input* windvx_input        = element->GetInput(SurfaceforcingsWindVxEnum); _assert_(windvx_input);
 	Input* windvy_input        = element->GetInput(SurfaceforcingsWindVyEnum); _assert_(windvy_input);
-	Input* thickness_input     = element->GetInput(SeaiceThicknessEnum);       _assert_(thickness_input);
+	Input* oceanvx_input       = element->GetInput(BasalforcingsOceanVxEnum);  _assert_(oceanvx_input);
+	Input* oceanvy_input       = element->GetInput(BasalforcingsOceanVyEnum);  _assert_(oceanvy_input);
+	Input* oceanssh_input      = element->GetInput(BasalforcingsOceanSshEnum); _assert_(oceanssh_input);
 
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=element->NewGauss(2);
@@ -167,20 +194,53 @@
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
 		element->NodalFunctions(basis, gauss);
 
+		/*Get all inputs on gauss point*/
 		thickness_input->GetInputValue(&thickness,gauss);
 		concentration_input->GetInputValue(&concentration,gauss);
+		coriolis_fact_input->GetInputValue(&coriolis_factor,gauss);
 		vx_input->GetInputValue(&vx,gauss);
 		vy_input->GetInputValue(&vy,gauss);
+		vxstar_input->GetInputValue(&vxstar,gauss);
+		vystar_input->GetInputValue(&vystar,gauss);
 		windvx_input->GetInputValue(&windvx,gauss);
 		windvy_input->GetInputValue(&windvy,gauss);
+		oceanvx_input->GetInputValue(&oceanvx,gauss);
+		oceanvy_input->GetInputValue(&oceanvy,gauss);
+		oceanssh_input->GetInputDerivativeValue(&dssh[2],xyz_list,gauss);
 
+		/*Previous speed (finite differencing)*/
+		for(int i=0;i<numnodes;i++){
+			pe->values[i*2+0]+=Jdet*gauss->weight*rho_ice*thickness*vx*basis[i];
+			pe->values[i*2+1]+=Jdet*gauss->weight*rho_ice*thickness*vy*basis[i];
+		}
+
 		/*Atmespheric forcing*/
-		wind_vnorm = sqrt(windvx*windvx + windvy*windvy);
+		vnorm    = sqrt(windvx*windvx + windvy*windvy);
+		constant_part = concentration*air_coef*rho_air*(air_lin_drag_coef+air_quad_drag_coef*vnorm);
+		for(int i=0;i<numnodes;i++){
+			pe->values[i*2+0] += dt*constant_part*windvx*Jdet*gauss->weight*basis[i];
+			pe->values[i*2+1] += dt*constant_part*windvy*Jdet*gauss->weight*basis[i];
+		}
 
+		/*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);
 		for(int i=0;i<numnodes;i++){
-			pe->values[i*2+0] += dt*concentration*air_coef*rho_air*(air_lin_drag_coef+air_quad_drag_coef*wind_vnorm)*windvx*Jdet*gauss->weight*basis[i];
-			pe->values[i*2+1] += dt*concentration*air_coef*rho_air*(air_lin_drag_coef+air_quad_drag_coef*wind_vnorm)*windvy*Jdet*gauss->weight*basis[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];
+
+			pe->values[i*2+0] += dt*constant_part*oceanvx*cos(ocean_turning_angle)*Jdet*gauss->weight*basis[i];
+			pe->values[i*2+1] += dt*constant_part*oceanvy*cos(ocean_turning_angle)*Jdet*gauss->weight*basis[i];
 		}
+
+		/*Coriolis forces (use ustar)*/
+		for(int i=0;i<numnodes;i++){
+			pe->values[i*2+0] += dt*(-rho_ice*thickness*coriolis_factor*vystar)*Jdet*gauss->weight*basis[i];
+			pe->values[i*2+1] += dt*(+rho_ice*thickness*coriolis_factor*vxstar)*Jdet*gauss->weight*basis[i];
+
+			pe->values[i*2+0] += dt*(-rho_ice*thickness*gravity*dssh[0])*Jdet*gauss->weight*basis[i];
+			pe->values[i*2+1] += dt*(-rho_ice*thickness*gravity*dssh[1])*Jdet*gauss->weight*basis[i];
+		}
 	}
 
 	/*Clean up and return*/
