Index: /issm/trunk-jpl/src/m/contrib/gravity/vfsa_mpi.cpp
===================================================================
--- /issm/trunk-jpl/src/m/contrib/gravity/vfsa_mpi.cpp	(revision 18490)
+++ /issm/trunk-jpl/src/m/contrib/gravity/vfsa_mpi.cpp	(revision 18491)
@@ -169,5 +169,5 @@
 	// Pp->Echo();
 
-	double  rhoi = 890;           /* ice density     */
+	double  rhoi = 917;           /* ice density     */
 	double  rhow = 1030;          /* water density   */
 	// double  rhos = 2013;		      /* sediment density */
@@ -186,5 +186,5 @@
 
 	double ctr=1;            /* parameter for filtering */
-	double sd=0.1;
+	double sd=0.5;
 
 	Matrix *xobs= new Matrix(ny,nx);
@@ -669,67 +669,90 @@
 	Matrix* m1grsm=new Matrix(my,mx);
 	Matrix* m1col=new Matrix(mx*my,1);
+	Matrix* m1gr2=new Matrix(my,mx);
+	Matrix* m1grsm2=new Matrix(my,mx);
+	Matrix* m1col2=new Matrix(mx*my,1);
+	Matrix* nptflag= new Matrix(mx*my,1);
 	double u=0;
 	double y=0;
 	m1->MatrixEqual(m0);
+	nptflag->MatrixSum(icethick,bathy);
+	/* first layer: ice */
 	for (int i=0;i<mx*my;i++){
-		if(icethick->GetValue(i,0)==0){
-			u=double (rand())/ double(RAND_MAX);
+		if(nptflag->GetValue(i,0)==0){
+			u=double(rand())/double(RAND_MAX);
 			y=signe(u-0.5)*T*(pow(1+1/T,fabs(2*u-1))-1);
 			m1->SetValue(i,1,m0->GetValue(i,1)+y*ptval);
+			m1->SetValue(i,2,m1->GetValue(i,1)+1e-10);
 			if(m1->GetValue(i,1)<=m1->GetValue(i,0)){
 				m1->SetValue(i,1,m1->GetValue(i,0)+1e-10);
+				m1->SetValue(i,2,m1->GetValue(i,1)+1e-10);
 			}
 			if(m1->GetValue(i,1)>=m1->GetValue(i,0)+mmax){
 				m1->SetValue(i,1,m1->GetValue(i,0)+mmax);
-			}
-		}
-	}
+				m1->SetValue(i,2,m1->GetValue(i,1)+1e-10);
+			}
+		}
+	}
+
 	m1->ExtractColumn(m1col,1);
 	vec2gridsimple(m1col,m1gr,mx,my);
 	filtergrav(m1grsm,m1gr,ctr,sd,mx,my);
 	reshape(m1grsm,m1col,mx,my);
+	m1->ExtractColumn(m1col2,2);
+	vec2gridsimple(m1col2,m1gr2,mx,my);
+	filtergrav(m1grsm2,m1gr2,ctr,sd,mx,my);
+	reshape(m1grsm2,m1col2,mx,my);
+
 	for (int i=0;i<mx*my;i++){
-		if(icethick->GetValue(i,0)==0){
+		if(nptflag->GetValue(i,0)==0){
 			m1->SetValue(i,1,m1col->GetValue(i,0));
+			m1->SetValue(i,2,m1col2->GetValue(i,0));
+			if(m1->GetValue(i,1)<m1->GetValue(i,0)){
+				m1->SetValue(i,1,m1->GetValue(i,0)+1e-10);
+				m1->SetValue(i,2,m1->GetValue(i,1)+1e-10);
+			}
 		}
 		else{
 			m1->SetValue(i,1,m0->GetValue(i,1));
-		}
-		if(m1->GetValue(i,1)<=m1->GetValue(i,0)){
-			m1->SetValue(i,1,m1->GetValue(i,0)+1e-10);
-		}
-	}
-
-	for (int i=0;i<mx*my;i++){
-		if(bathy->GetValue(i,0)==0){
-			u=double (rand())/ double(RAND_MAX);
-			y=signe(u-0.5)*T*(pow(1+1/T,fabs(2*u-1))-1);
-			m1->SetValue(i,2,m0->GetValue(i,2)+y*ptval);
-			if(m1->GetValue(i,2)<=m1->GetValue(i,1)){
-				m1->SetValue(i,2,m1->GetValue(i,1)+1e-10);
-			}
-			if(m1->GetValue(i,2)>=m1->GetValue(i,1)+mmax2){
-				m1->SetValue(i,2,m1->GetValue(i,1)+mmax2);
-			}
-		}
-	}
-	m1->ExtractColumn(m1col,2);
-	vec2gridsimple(m1col,m1gr,mx,my);
-	filtergrav(m1grsm,m1gr,ctr,sd,mx,my);
-	reshape(m1grsm,m1col,mx,my);
-	for (int i=0;i<mx*my;i++){
-		if(bathy->GetValue(i,0)==0){
-			m1->SetValue(i,2,m1col->GetValue(i,0));
-		}
-		else{
 			m1->SetValue(i,2,m0->GetValue(i,2));
 		}
-		if(m1->GetValue(i,2)<=m1->GetValue(i,1)){
-			m1->SetValue(i,2,m1->GetValue(i,1)+1e-10);
-		}
-	}
+	}
+
+				/* second layer: water */
+//	for (int i=0;i<mx*my;i++){
+//		if(bathy->GetValue(i,0)==0){
+//			u=double (rand())/ double(RAND_MAX);
+//			y=signe(u-0.5)*T*(pow(1+1/T,fabs(2*u-1))-1);
+//			m1->SetValue(i,2,m0->GetValue(i,2)+y*ptval);
+//			if(m1->GetValue(i,2)<=m1->GetValue(i,1)){
+//				m1->SetValue(i,2,m1->GetValue(i,1)+1e-10);
+//			}
+//			if(m1->GetValue(i,2)>=m1->GetValue(i,1)+mmax2){
+//				m1->SetValue(i,2,m1->GetValue(i,1)+mmax2);
+//			}
+//		}
+//	}
+//	m1->ExtractColumn(m1col,2);
+//	vec2gridsimple(m1col,m1gr,mx,my);
+//	filtergrav(m1grsm,m1gr,ctr,sd,mx,my);
+//	reshape(m1grsm,m1col,mx,my);
+//	for (int i=0;i<mx*my;i++){
+//		if(bathy->GetValue(i,0)==0){
+//			m1->SetValue(i,2,m1col->GetValue(i,0));
+//		}
+//		else{
+//			m1->SetValue(i,2,m0->GetValue(i,2));
+//		}
+//		if(m1->GetValue(i,2)<=m1->GetValue(i,1)){
+//			m1->SetValue(i,2,m1->GetValue(i,1)+1e-10);
+//		}
+//	}
 	delete m1gr;
 	delete m1grsm;
 	delete m1col;
+	delete m1gr2;
+	delete m1grsm2;
+	delete m1col2;
+	delete nptflag;
 }/*}}}*/
 double signe(double a){/*{{{*/