Index: param.ccl
===================================================================
--- param.ccl	(revision 102)
+++ param.ccl	(working copy)
@@ -39,6 +39,10 @@
 {
 } "yes"
 
+CCTK_REAL adm_tol "Tolerance of ADM masses when give_bare_mass=no"
+{
+  (0:*) :: ""
+} 1.0e-10
 
 KEYWORD grid_setup_method "How to fill the 3D grid from the spectral grid"
 {
Index: src/TwoPunctures.c
===================================================================
--- src/TwoPunctures.c	(revision 102)
+++ src/TwoPunctures.c	(working copy)
@@ -216,7 +216,7 @@
 #endif
   static CCTK_REAL *F = NULL;
   static derivs u, v;
-  CCTK_REAL admMass, M_p, M_m, tmp_p, tmp_m, um, up, new_mass, old_mass;
+  CCTK_REAL admMass;
 
   if (! F) {
     /* Solve only when called for the first time */
@@ -251,59 +251,75 @@
       set_initial_guess(cctkGH, v);
     }
 
+    /* If bare masses are not given, iteratively solve for them given the 
+       target ADM masses target_M_plus and target_M_minus and with initial 
+       guesses given by par_m_plus and par_m_minus. */
     if(!(give_bare_mass)) {
-      * mp  = target_M_plus;
-      * mm  = target_M_minus;
-      
-      M_p= target_M_plus;
-      M_m= target_M_minus;
+      CCTK_REAL tmp_p, tmp_m, Mp_adm, Mm_adm, Mp_adm_err, Mm_adm_err, up, um;
+      char valbuf[100];
 
-      new_mass = 0;
-      old_mass = 1.;
-      while (new_mass<old_mass) {
-	old_mass = (double)*mp + (double)*mm;
-	Newton (cctkGH, nvar, n1, n2, n3, v, Newton_tol, 1);
+      CCTK_REAL M_p = target_M_plus;
+      CCTK_REAL M_m = target_M_minus;
 
-	F_of_v (cctkGH, nvar, n1, n2, n3, v, F, u);
+      CCTK_VInfo (CCTK_THORNSTRING, "Attempting to find bare masses.");
+      CCTK_VInfo (CCTK_THORNSTRING, "Target ADM masses: M_p=%g and M_m=%g",
+                  (double) M_p, (double) M_m);
+      CCTK_VInfo (CCTK_THORNSTRING, "ADM mass tolerance: %g", (double) adm_tol);
 
-	up  = PunctIntPolAtArbitPosition 
-	  (0, nvar, n1, n2, n3, v, par_b, 0., 0.);
-	um  = PunctIntPolAtArbitPosition 
-	  (0, nvar, n1, n2, n3, v, -par_b, 0., 0.);
+      /* Loop until both ADM masses are within adm_tol of their target */
+      do {
+        CCTK_VInfo (CCTK_THORNSTRING, "Bare masses: mp=%g, mm=%g",
+                    (double)*mp, (double)*mm);
+        Newton (cctkGH, nvar, n1, n2, n3, v, Newton_tol, 1);
 
-	for(int l=0; l<32;l++) {
-	  tmp_p=(4*par_b*um+ *mp +4*par_b);
-	  tmp_m=(4*par_b*up+ *mm +4*par_b);
-	  *mp = *mp - ((*mp*(up+1+M_m/tmp_p)-M_p)/
-		       (up+1+(M_m/tmp_p)-(*mp*M_m)/(pow(tmp_p,2))));
-	  *mm = *mm - ((*mm*(um+1+M_p/tmp_m)-M_m)/
-		       (um+1+(M_p/tmp_m)-(*mm*M_p)/(pow(tmp_m,2))));
-	}
-	char valbuf_p[100], valbuf_m[100];  
-	
-	sprintf (valbuf_p,"%f", (float) *mp);
-	CCTK_ParameterSet ("par_m_plus", "twopunctures", valbuf_p);
-	sprintf (valbuf_m,"%f", (float) *mm);
-	CCTK_ParameterSet ("par_m_minus", "twopunctures", valbuf_m);
-	new_mass = (double)*mp + (double)*mm;
-      }
-    } 
-    admMass = (*mp + *mm
-               - 4*par_b*PunctEvalAtArbitPosition(v.d0, 1, 0, 0, n1, n2, n3));
-    CCTK_VInfo (CCTK_THORNSTRING, "ADM mass is %g\n", (double)admMass);
+        F_of_v (cctkGH, nvar, n1, n2, n3, v, F, u);
 
+        up = PunctIntPolAtArbitPosition(0, nvar, n1, n2, n3, v, par_b, 0., 0.);
+        um = PunctIntPolAtArbitPosition(0, nvar, n1, n2, n3, v,-par_b, 0., 0.);
+
+        /* Calculate the ADM masses from to the current bare mass guess */
+        Mp_adm = (1 + up) * *mp + *mp * *mm / (4. * par_b);
+        Mm_adm = (1 + um) * *mm + *mp * *mm / (4. * par_b);
+
+        /* Check how far the current ADM masses are from the target */
+        Mp_adm_err = fabs(M_p-Mp_adm);
+        Mm_adm_err = fabs(M_m-Mm_adm);
+        CCTK_VInfo (CCTK_THORNSTRING, "ADM mass error: M_p_err=%.4g, M_m_err=%.4g",
+                    (double)Mp_adm_err, (double)Mm_adm_err);
+        
+        /* Invert the ADM mass equation and update the bare mass guess so that
+           it gives the correct target ADM masses */
+        tmp_p = -4*par_b*( 1 + um + up + um*up ) + 
+                sqrt(16*par_b*M_m*(1 + um)*(1 + up) + 
+                  pow(-M_m + M_p + 4*par_b*(1 + um)*(1 + up),2));
+        *mp = (tmp_p + M_p - M_m)/(2.*(1 + up));
+        *mm = (tmp_p - M_p + M_m)/(2.*(1 + um));
+        
+        /* Set the par_m_plus and par_m_minus parameters */
+        sprintf (valbuf,"%.17g", (double) *mp);
+        CCTK_ParameterSet ("par_m_plus", "twopunctures", valbuf);
+        
+        sprintf (valbuf,"%.17g", (double) *mm);
+        CCTK_ParameterSet ("par_m_minus", "twopunctures", valbuf);
+        
+      } while ( (Mp_adm_err > adm_tol) ||
+                (Mm_adm_err > adm_tol) );
+                
+      CCTK_VInfo (CCTK_THORNSTRING, "Found bare masses.");
+    }
+
     Newton (cctkGH, nvar, n1, n2, n3, v, Newton_tol, Newton_maxit);
  
     F_of_v (cctkGH, nvar, n1, n2, n3, v, F, u);
 
     CCTK_VInfo (CCTK_THORNSTRING,
-		  "The two puncture masses are %g and %g",
-                (double) par_m_minus, (double) par_m_plus);
+		  "The two puncture masses are mp=%.17g and mm=%.17g",
+                (double) *mp, (double) *mm);
 
     /* print out ADM mass, eq.: \Delta M_ADM=2*r*u=4*b*V for A=1,B=0,phi=0 */
     admMass = (*mp + *mm
-               - 4*par_b*PunctEvalAtArbitPosition(v.d0, 1, 0, 0, n1, n2, n3));
-    CCTK_VInfo (CCTK_THORNSTRING, "ADM mass is %g", (double)admMass);
+               - 4*par_b*PunctEvalAtArbitPosition(v.d0, 1, 0, 0, n1, n2, n3));;
+    CCTK_VInfo (CCTK_THORNSTRING, "The total ADM mass is %g", (double) admMass);
   }
 
   if (CCTK_EQUALS(grid_setup_method, "Taylor expansion"))
Index: src/Newton.c
===================================================================
--- src/Newton.c	(revision 102)
+++ src/Newton.c	(working copy)
@@ -493,7 +493,7 @@
 Newton (CCTK_POINTER_TO_CONST const cctkGH,
         int const nvar, int const n1, int const n2, int const n3,
 	derivs v,
-        CCTK_REAL const tol, int const itmax)
+        CCTK_REAL const tol, int const itmax, int const output)
 {
   DECLARE_CCTK_PARAMETERS;
 
@@ -518,8 +518,12 @@
 #pragma omp parallel for
     for (int j = 0; j < ntotal; j++)
       dv.d0[j] = 0;
-    printf ("Newton: it=%d \t |F|=%e\n", it, (double)dmax);
-    printf ("bare mass: mp=%g \t mm=%g\n", (double) par_m_plus, (double) par_m_minus);
+
+    if(verbose==1){
+      printf ("Newton: it=%d \t |F|=%e\n", it, (double)dmax);
+      printf ("bare mass: mp=%g \t mm=%g\n", (double) par_m_plus, (double) par_m_minus);
+    }
+
     fflush(stdout);
     ii =
       bicgstab (cctkGH,
@@ -536,7 +540,10 @@
       F_of_v (cctkGH, nvar, n1, n2, n3, v, F, u);
       dmax = norm_inf (F, ntotal);
   }
-  printf ("Newton: it=%d \t |F|=%e \n", it, (double)dmax);
+
+  if(verbose==1)
+    printf ("Newton: it=%d \t |F|=%e \n", it, (double)dmax);
+
   fflush(stdout);
 
   free_dvector (F, 0, ntotal - 1);