Index: /issm/trunk-jpl/src/m/contrib/gravity/vfsa_mpi.cpp
===================================================================
--- /issm/trunk-jpl/src/m/contrib/gravity/vfsa_mpi.cpp	(revision 18549)
+++ /issm/trunk-jpl/src/m/contrib/gravity/vfsa_mpi.cpp	(revision 18550)
@@ -133,5 +133,5 @@
 double signe(double a);
 void filtergrav(Matrix* A,Matrix* Ain,double ctr,double sd,int mx,int my);
-void newmodelgen(Matrix* m0,Matrix* m1,Matrix* bathy,Matrix* icethick,int mx,int my,double T,double ptval,double mmax,double mmax2,double ctr,double sd);
+void newmodelgen(Matrix* m0,Matrix* m1,Matrix* bathy,Matrix* icethick,int mx,int my,double T,double ptval,double mmax,double mmax2,double ctr,double sd, Matrix *landmask);
 double coolshed(double T0,double k,double c,double D);
 /*}}}*/
@@ -165,9 +165,9 @@
 
 
-	Matrix *Pp=new Matrix(mx*my,2); /* data positions */
+	Matrix *Pp=new Matrix(mx*my,2); /* prisms positions */
 	makep(Pp,mx,my,dx,dy);
 	// Pp->Echo();
 
-	double  rhoi = 917;           /* ice density     */
+	double  rhoi = 890;           /* ice density     */
 	double  rhow = 1030;          /* water density   */
 	// double  rhos = 2013;		      /* sediment density */
@@ -185,4 +185,6 @@
 	rho2->SetValue(0,1,rhow-rhoc);
 
+	double dlevel=860;         /* level of data acquisition */
+
 	double ctr=1;            /* parameter for filtering */
 	double sd=0.5;
@@ -195,5 +197,4 @@
 	//	yobs->Echo();
 
-	double dlevel= 860;                /* mean level of data acquisition */
 
 	double mmax  = 1000;               /* max value for layer interfaces */
@@ -223,9 +224,19 @@
 	/* load the data {{{*/
 
+	/*landmask */
+
+	ifstream file("landmask.txt");
+	Matrix * landmask= new Matrix(nx*ny,1);
+	double inputnumber;
+	for(int i=0;i<ny*nx; i++){ 
+		file >> inputnumber;
+		landmask->SetValue(i,0,inputnumber);
+	}
+	file.close();
+
 	/* Observed gravity anomaly */
 
 	ifstream file1("store_fa500_36s.txt");
 	Matrix * gobs= new Matrix(nx*ny,1);
-	double inputnumber;
 	for(int i=0;i<ny*nx; i++){ 
 		file1 >> inputnumber;
@@ -261,6 +272,6 @@
 
 	ifstream file4("store_flag_eval500.txt");
-	Matrix * evalid= new Matrix(mx*my,1);
-	for(int s=0;s<mx*my; s++){ 
+	Matrix * evalid= new Matrix(nx*ny,1);
+	for(int s=0;s<nx*ny; s++){ 
 		file4 >> inputnumber;
 		evalid->SetValue(s,0,inputnumber);
@@ -346,5 +357,5 @@
 		/* update model and calculate energy */
 
-		newmodelgen(m_old,m_new,bathy,icethick,mx,my,Tm,ptval,mmax,mmax2,ctr,sd);  /* new model */
+		newmodelgen(m_old,m_new,bathy,icethick,mx,my,Tm,ptval,mmax,mmax2,ctr,sd, landmask);  /* new model */
 		E_new=misfit(m_new,evalid,gobs,dlevel,Pobs,xobs,yobs,Pp,rho1,rho2,dx,dy,dn,nx,ny,mx,my, my_rank, num_procs); /* new energy */
 		dE=E_new-E_old;                                        /* energy difference */
@@ -508,4 +519,6 @@
 					m2->SetValue(i,j,i*1e-10);
 				}
+				m1->SetValue(i,j,m1->GetValue(i,j));
+				m2->SetValue(i,j,m2->GetValue(i,j));
 			}
 			else{
@@ -619,7 +632,7 @@
 	Matrix* g=new Matrix(nx*ny,1);
 	Matrix* gcalgr=new Matrix(ny,nx);
-	Matrix* gcalvec=new Matrix(mx*my,1);
-	Matrix* df=new Matrix(mx*my,1);
-	Matrix* G=new Matrix(mx*my,3);
+	Matrix* gcalvec=new Matrix(nx*ny,1);
+	Matrix* df=new Matrix(nx*ny,1);
+	Matrix* G=new Matrix(nx*ny,3);
 	double a=0;
 	double b=0;
@@ -630,6 +643,6 @@
 	g->MatrixSum(g1,g2);
 	vec2gridsimple(g,gcalgr,nx,ny);
-	reshape(gcalgr,gcalvec,mx,my);
-	for (int i=0;i<mx*my;i++){
+	reshape(gcalgr,gcalvec,nx,ny);
+	for (int i=0;i<nx*ny;i++){
 		df->SetValue(i,0,evalid->GetValue(i,0)*(gobs->GetValue(i,0)-gcalvec->GetValue(i,0)));
 		G->SetValue(i,0,evalid->GetValue(i,0)*Pobs->GetValue(i,0));
@@ -640,11 +653,11 @@
 	GSLsquarefit(&M,G,df);
 
-	for (int i=0;i<my;i++){
-		for(int j=0;j<mx;j++){
+	for (int i=0;i<ny;i++){
+		for(int j=0;j<nx;j++){
 			gcalgr->SetValue(i,j,gcalgr->GetValue(i,j)+xobs->GetValue(i,j)*M->GetValue(0,0)+yobs->GetValue(i,j)*M->GetValue(1,0)+M->GetValue(2,0));
 		}
 	}
-	reshape(gcalgr,g,mx,my);
-	for (int i=0;i<mx*my;i++){
+	reshape(gcalgr,g,nx,ny);
+	for (int i=0;i<nx*ny;i++){
 		a=a+fabs(evalid->GetValue(i,0)*(gobs->GetValue(i,0)-g->GetValue(i,0)));
 		b=b+fabs(evalid->GetValue(i,0)*(gobs->GetValue(i,0)+g->GetValue(i,0)));
@@ -665,5 +678,5 @@
 	return e;
 }/*}}}*/
-void newmodelgen(Matrix* m0,Matrix* m1,Matrix* bathy,Matrix* icethick,int mx,int my,double T,double ptval,double mmax,double mmax2,double ctr,double sd){/*{{{*/
+void newmodelgen(Matrix* m0,Matrix* m1,Matrix* bathy,Matrix* icethick,int mx,int my,double T,double ptval,double mmax,double mmax2,double ctr,double sd, Matrix *landmask){/*{{{*/
 	Matrix* m1gr=new Matrix(my,mx);
 	Matrix* m1grsm=new Matrix(my,mx);
@@ -679,16 +692,31 @@
 	/* first layer: ice */
 	for (int i=0;i<mx*my;i++){
-		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);
+		if(landmask->GetValue(i,0)==2){
+			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)+mmax){
-				m1->SetValue(i,1,m1->GetValue(i,0)+mmax);
-				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);
+				}
+			}
+		}
+		else if(landmask->GetValue(i,0)==0){
+			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,2,m0->GetValue(i,2)+y*ptval);
+				if(m1->GetValue(i,2)<=m1->GetValue(i,0)){
+					m1->SetValue(i,2,m1->GetValue(i,0)+1e-10);
+				}
+				if(m1->GetValue(i,2)>=m1->GetValue(i,0)+mmax2){
+					m1->SetValue(i,2,m1->GetValue(i,0)+mmax2);
+				}
 			}
 		}
@@ -705,15 +733,49 @@
 
 	for (int i=0;i<mx*my;i++){
-		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));
-			m1->SetValue(i,2,m0->GetValue(i,2));
+		if(landmask->GetValue(i,0)==2){
+			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);
+				}
+				if(fabs(m1->GetValue(i,2)-m1->GetValue(i,1))>1){
+					m1->SetValue(i,2,m1->GetValue(i,1)+1e-10);
+				}
+			}
+			else{
+				m1->SetValue(i,1,m0->GetValue(i,1));
+				m1->SetValue(i,2,m0->GetValue(i,2));
+			}
+		}
+		else if(landmask->GetValue(i,0)==0){
+			if(nptflag->GetValue(i,0)==0){
+				m1->SetValue(i,2,m1col2->GetValue(i,0));
+				if(m1->GetValue(i,2)<=m1->GetValue(i,0)){
+					m1->SetValue(i,2,m1->GetValue(i,0)+1e-10);
+				}
+				if(fabs(m1->GetValue(i,0)-m1->GetValue(i,1))>1){
+					m1->SetValue(i,1,m1->GetValue(i,0)+1e-10);
+				}
+			}
+			else{
+				m1->SetValue(i,1,m0->GetValue(i,1));
+				m1->SetValue(i,2,m0->GetValue(i,2));
+			}
+		}
+		else {
+			if(nptflag->GetValue(i,0)==0){
+				if(fabs(m1->GetValue(i,0)-m1->GetValue(i,1))>1){
+					m1->SetValue(i,1,m1->GetValue(i,0));
+				}
+				if(fabs(m1->GetValue(i,0)-m1->GetValue(i,2))>1){
+					m1->SetValue(i,2,m1->GetValue(i,0));
+				}
+			}
+			else{
+				m1->SetValue(i,1,m0->GetValue(i,1));
+				m1->SetValue(i,2,m0->GetValue(i,2));
+			}
 		}
 	}