Several changes:

1) Minor changes to messages to screen.
2) Changing counter for multiple horizons so that it starts counting
   from 1 and not from 0.
3) The most important change is that I had to reintroduce some goto
   statements that I deleted in my previous commit.  The code was not
   working properly any more.  I will think a bit more about how to
   eliminate them properly.
m


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinder/trunk@176 89daf98e-ef62-4674-b946-b8ff9de2216c

1 files changed, 45 insertions, 31 deletions

diff --git a/src/AHFinder.F b/src/AHFinder.F
index 15522c4..3b00dc6 100644
--- a/src/AHFinder.F
+++ b/src/AHFinder.F
@@ -257,12 +257,25 @@
 !     anything that is likely to be different for the different
 !     horizons should NOT be initialized until after this point!
 
-      do mfind=0,nhorizon-1
+      do mfind=1,nhorizon
 
 !     Say what we are doing.
 
-      write(*,*)
-      write(*,*) 'AHFinder: Searching for horizon ',mfind
+      if (nhorizon.gt.1) then
+         if (mfind.eq.1) then
+            write(*,*)
+            write(*,*)
+            write(*,*) 'AHFinder: MULTIPLE HORIZON MODE!'
+         end if
+         write(*,*)
+         write(*,*) 'Searching for horizon ',mfind
+      else
+         write(*,*)
+         write(*,*)
+         write(*,*) 'AHFinder: SINGLE HORIZON MODE!'
+         write(*,*)
+         write(*,*) 'Searching for horizon'
+      end if
 
 
 !     ********************************
@@ -288,7 +301,7 @@
 
       if (find3) then
             
-         if (mfind.eq.0) then
+         if (mfind.eq.1) then
 
             logf = filestr(1:nfile)//"/ahf_logfile_0"
             almf = filestr(1:nfile)//"/ahf_coeff_0.alm"
@@ -301,7 +314,7 @@
             merip1f = filestr(1:nfile)//"/ahf_meri_p1_0.tl"
             merip2f = filestr(1:nfile)//"/ahf_meri_p2_0.tl"
 
-         else if (mfind.eq.1) then
+         else if (mfind.eq.2) then
 
             logf = filestr(1:nfile)//"/ahf_logfile_1"
             almf = filestr(1:nfile)//"/ahf_coeff_1.alm"
@@ -395,11 +408,11 @@
 !     **********************************
 
       if (find3) then
-         if (mfind.eq.0) then
+         if (mfind.eq.1) then
             xc = ahf_xc_0
             yc = ahf_yc_0
             zc = ahf_zc_0
-         else if (mfind.eq.1) then   
+         else if (mfind.eq.2) then   
             xc = ahf_xc_1
             yc = ahf_yc_1
             zc = ahf_zc_1
@@ -415,24 +428,24 @@
       end if
 
 !     Check if the center of the expansion is
-!     out of the grid.
+!     out of the grid.  If it is, give up.
 
       if ((xc.gt.xmx).or.(xc.lt.xmn)) then
          write(*,*)
-         write(*,*) 'xc is out of the grid, setting it to zero.'
-         xc = 0.0D0
+         write(*,*) 'xc is out of the grid, giving up.'
+         goto 100
       end if
 
       if ((yc.gt.ymx).or.(yc.lt.ymn)) then
          write(*,*)
-         write(*,*) 'yc is out of the grid, setting it to zero.'
-         yc = 0.0D0
+         write(*,*) 'yc is out of the grid, giving up.'
+         goto 100
       end if
 
       if ((zc.gt.zmx).or.(zc.lt.zmn)) then
          write(*,*)
-         write(*,*) 'zc is out of the grid, setting it to zero.'
-         zc = 0.0D0
+         write(*,*) 'zc is out of the grid, giving up.'
+         goto 100
       end if
 
 
@@ -728,9 +741,9 @@
 !     has been forced by the parameter file.
 
       if (find3) then
-         if (mfind.eq.0) then
+         if (mfind.eq.1) then
             r0 = ahf_r0_0
-         else if (mfind.eq.1) then              
+         else if (mfind.eq.2) then              
             r0 = ahf_r0_1
          else
             r0 = ahf_r0_2
@@ -795,7 +808,7 @@
 
          if ((myproc.eq.0).and.verbose) then
             write(*,*)
-            write(*,*)  'AHFinder: Using old horizon as initial guess'
+            write(*,*)  'Using old horizon as initial guess'
          end if
 
       else if (manual_guess) then
@@ -804,7 +817,7 @@
 
          if ((myproc.eq.0).and.verbose) then
             write(*,*)
-            write(*,*)  'AHFinder: Manual guess'
+            write(*,*)  'Manual guess'
          end if
 
          c0_temp(0) = ahf_l0_guess
@@ -877,7 +890,7 @@
 
 !        Flow algorithm.
 
-         if ((ahf_ncall.eq.1).and.(mfind.eq.0)) then
+         if ((ahf_ncall.eq.1).and.(mfind.eq.1)) then
             allocate(hflow0(0:lmax),cflow0(0:lmax),nflow0(0:lmax))
             allocate(hflowc(lmax,lmax),cflowc(lmax,lmax),nflowc(lmax,lmax))
             allocate(hflows(lmax,lmax),cflows(lmax,lmax),nflows(lmax,lmax))
@@ -893,8 +906,6 @@
 !     ***   WRITE MESSAGES AND LOG FILE   ***
 !     ***************************************
 
- 5    continue
-
       status_old = status
 
       if (status) then
@@ -1043,7 +1054,7 @@
 !           M = sqrt( A / (16 pi) ) = 0.141047396 sqrt(A)
 
             write(*,*)
-            write(*,*) 'AHFinder: Surface found.'
+            write(*,*) 'Surface found, details below.'
             write(*,*)
             write(*,"(A21,ES14.6)") ' Surface area       =',intarea_h
             write(*,"(A21,ES14.6)") ' Surface mass       =',
@@ -1062,7 +1073,7 @@
 
             if (logfile) then
                write(1,*)
-               write(1,*) 'AHFinder: Surface found.'
+               write(1,*) 'Surface found, details below'
                write(1,*)
                write(1,"(A21,ES14.6)") ' Surface area       =',intarea_h
                write(1,"(A21,ES14.6)") ' Surface mass       =',
@@ -1313,12 +1324,11 @@
 !           Write messages.
 
             write(*,*)
-            write(*,*) 'AHFinder: No horizon found.'
-            write(*,*)
+            write(*,*) 'No horizon found.'
 
             if (logfile) then
                write(1,*)
-               write(1,*) 'AHFinder: No horizon found.'
+               write(1,*) 'No horizon found.'
                write(1,*)
             end if
             
@@ -1549,7 +1559,7 @@
 
       if (ahf_2Doutput.ne.0) then
          if (find3) then
-            if (mfind.eq.2) then
+            if (mfind.eq.3) then
                call CCTK_OutputVarAsByMethod(ierror,cctkGH,
      .              "ahfinder::ahfgrid3","IOFlexIO_2D","ahfgrid3")
                call CCTK_OutputVarAsByMethod(ierror,cctkGH,
@@ -1581,7 +1591,7 @@
 
          if (find3) then
 
-            if (mfind.eq.0) then
+            if (mfind.eq.1) then
 
                dhole1_xmin = xc
                dhole1_ymin = yc
@@ -1591,7 +1601,7 @@
                dhole1_ymax = yc
                dhole1_zmax = zc
 
-            else if (mfind.eq.1) then
+            else if (mfind.eq.2) then
 
                dhole2_xmin = xc
                dhole2_ymin = yc
@@ -1642,7 +1652,7 @@
 
          if (find3) then
 
-            if (mfind.eq.0) then
+            if (mfind.eq.1) then
 
 !              First hole.
 
@@ -1654,7 +1664,7 @@
                if (dhole1_ymax.lt.(yc+0.5D0*rhor)) dhole1_ymax = yc+0.5D0*rhor
                if (dhole1_zmax.lt.(zc+0.5D0*rhor)) dhole1_zmax = zc+0.5D0*rhor
 
-            else if (mfind.eq.1) then
+            else if (mfind.eq.2) then
 
 !              Second hole.
 
@@ -1699,6 +1709,10 @@
 !     ***   END OF LOOP FOR LOOKING FOR SEVERAL HORIZONS   ***
 !     ********************************************************
 
+ 100  continue
+
+      write (*,*)
+
       end do