Index: /issm/trunk-jpl/src/m/contrib/gravity/vfsa.cpp
===================================================================
--- /issm/trunk-jpl/src/m/contrib/gravity/vfsa.cpp	(revision 18434)
+++ /issm/trunk-jpl/src/m/contrib/gravity/vfsa.cpp	(revision 18435)
@@ -131,4 +131,6 @@
 double misfit(Matrix* m0,Matrix* evalid,Matrix* gobs,int dlevel,Matrix* Pobs,Matrix* xobs,Matrix* yobs,Matrix* Pp,Matrix* rho1, Matrix* rho2,int dx,int dy,int dn,int nx,int ny, int mx,int my);
 void GSLsquarefit(Matrix** pX,Matrix* A,Matrix* B);
+double signe(double a);
+void filtergrav(Matrix* A,Matrix* Ain,double ctr,double sd,int mx,int my);
 /*}}}*/
 
@@ -281,5 +283,7 @@
 	M->Echo();
 	o->Echo();*/
-
+double a;
+a=signe(-9);
+cout<<a;
 	return 0;
 }/*}}}*/
@@ -491,2 +495,81 @@
 	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){/*{{{*/
+	Matrix* m1gr=new Matrix(my,mx);
+	Matrix* m1grsm=new Matrix(my,mx);
+	Matrix* m1col=new Matrix(mx*my,1);
+	double u=0;
+	double y=0;
+	m1->MatrixEqual(m0);
+	for (int i=0;i<mx*my;i++){
+		if(icethick->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);
+			if(m1->GetValue(i,1)<=m1->GetValue(i,0)){
+				m1->SetValue(i,1,m1->GetValue(i,0)+1e-10);
+			}
+			if(m1->GetValue(i,1)>=m1->GetValue(i,0)+mmax){
+				m1->SetValue(i,1,m1->GetValue(i,0)+mmax);
+			}
+		}
+	}
+	m1->ExtractColumn(m1col,1);
+	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(icethick->GetValue(i,0)==0){
+			m1->SetValue(i,1,m1col->GetValue(i,0));
+		}
+		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);
+		}
+	}
+	delete m1gr;
+	delete m1grsm;
+	delete m1col;
+}/*}}}*/
+double signe(double a){/*{{{*/
+	if(a<0){return -1;}
+	else{return 1;}
+}/*}}}*/
+void filtergrav(Matrix* A,Matrix* Ain,double ctr,double sd,int mx,int my){/*{{{*/
+	for (int i=1;i<my-1;i++){
+		for(int j=1;j<mx-1;j++){
+			A->SetValue(i,j,(ctr*Ain->GetValue(i,j)+sd*(Ain->GetValue(i-1,j)+Ain->GetValue(i+1,j)+Ain->GetValue(i,j-1)+Ain->GetValue(i,j+1)))/(ctr+4*sd));
+		}
+	}
+}/*}}}*/