Several changes:

1) Now I use the cactus routines ConfToPhys and PhysToConf instead of
doing it independently.

2) I have eliminated a goto statement that was causing the compilation
to be very slow (specially on Linux).

3) Some cleaning up.


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

1 files changed, 73 insertions, 111 deletions

diff --git a/src/AHFinder.F b/src/AHFinder.F
index f526f4f..b3f0297 100644
--- a/src/AHFinder.F
+++ b/src/AHFinder.F
@@ -34,6 +34,7 @@
       integer mtype
       integer getoutpfx
       integer istat
+      integer nhorizon
 
       CCTK_REAL atime
       CCTK_REAL rmx,r0
@@ -78,6 +79,8 @@
 !                  2.  Integral of H^2 is negative out and positive in.
 !                  3.  Integral of H^2 is negative out and in.
 !                  4.  Integral of H^2 is positive out and in.
+!
+!     nhorizon  Number of horizons to look for.
 
 
 !     *****************************************
@@ -120,9 +123,19 @@
       if (mod(int(cctk_iteration)-ahfafter,ahfso).ne.0) return
 
 
-!     *******************************************************
-!     ***   FIND OUT IF WE WANT TO FIND TRAPPED SURFACE   ***
-!     *******************************************************
+!     *******************************************
+!     ***   ADD TO COUNT OF CALLS TO FINDER   ***
+!     *******************************************
+
+!     Add to count of calls to finder.  Notice that this counter was
+!     initialized to zero, so it must be increased on first call to 1.
+
+      ahf_ncall = ahf_ncall + 1
+
+
+!     *********************************************************
+!     ***   FIND OUT IF WE WANT TO FIND A TRAPPED SURFACE   ***
+!     *********************************************************
 
       if (ahf_trapped_surface.ne.0) then
          find_trapped_surface = .true.
@@ -183,81 +196,80 @@
 
       if (conformal_state == CONFORMAL_METRIC) then
 
-         call CCTK_WARN(4,"Converting metric: conformal -> physical")
-
-         ahfconformal = .true.
-
+         ahfconformal    = .true.
          conformal_state = NOCONFORMAL_METRIC
 
-         gxx = gxx*psi**4
-         gxy = gxy*psi**4
-         gxz = gxz*psi**4
-         gyy = gyy*psi**4
-         gyz = gyz*psi**4
-         gzz = gzz*psi**4
+         call ConfToPhys(nx,ny,nz,psi,gxx,gxy,gxz,gyy,gyz,gzz)
 
       end if
 
 
-!     ***********************************************
-!     ***   CHECK IF WE WANT TO FIND 3 HORIZONS   ***
-!     ***********************************************
+!     ************************
+!     ***   FIND  {lmax}   ***
+!     ************************
 
-      if (ahf_find3.ne.0) then
-         find3 = .true.
-      else
-         find3 = .false.
+      lmax = ahf_lmax
+
+      if (lmax.ge.20) then
+         write(*,*)
+         write(*,*) 'ahf_lmax must be smaller than 20.'
+         write(*,*) 'STOPPED IN AHFinder.F'
+         write(*,*)
+         STOP
       end if
 
 
-!     *******************************************
-!     ***   ADD TO COUNT OF CALLS TO FINDER   ***
-!     *******************************************
+!     ***************************************************
+!     ***   ALLOCATE STORAGE FOR COEFFICIENT ARRAYS   ***
+!     ***************************************************
 
-!     IMPORTANT:  When looking for 3 horizons, the label 1 below
-!     marks the point where one comes back in order to look for
-!     the next horizon.
+      if (ahf_ncall.eq.1) then
 
-!     Because of this, anything that is likely to be different
-!     for the different horizons should NOT be initialized until
-!     after this point!
+         allocate(c0(0:lmax),c0_old(0:lmax))
+         allocate(cc(lmax,lmax),cc_old(lmax,lmax))
+         allocate(cs(lmax,lmax),cs_old(lmax,lmax))
 
- 1    continue
+         c0 = zero
+         cc = zero
+         cs = zero
 
-!     Say what we are doing.
+         c0_old = zero
+         cc_old = zero
+         cs_old = zero
 
-      if (find3) then
-         if (mfind.eq.0) then
-            write(*,*)
-            write(*,*) 'AHFinder: 3 HORIZONS MODE!'
-            write(*,*)
-            write(*,*) 'AHFinder: Searching for first horizon'
-         else if (mfind.eq.1) then
-            write(*,*)
-            write(*,*) 'AHFinder: Searching for second horizon'    
-         else if (mfind.eq.2) then
-            write(*,*)
-            write(*,*) 'AHFinder: Searching for third horizon' 
-         end if   
-      else
-         write(*,*)
-         write(*,*) 'AHFinder: Searching for horizon' 
       end if
 
-!     Add to count of calls to finder.  Notice that this counter was
-!     initialized to zero, so it must be increased on first call to 1.
-!     If we look for 3 horizons, this should be done before looking
-!     for the first one.
 
-      if (find3) then
-         if (mfind.eq.0) then
-            ahf_ncall = ahf_ncall + 1
-         end if
+!     ***********************************************
+!     ***   CHECK IF WE WANT TO FIND 3 HORIZONS   ***
+!     ***********************************************
+
+      if (ahf_find3.eq.0) then
+         find3 = .false.
+         nhorizon = 1
       else
-         ahf_ncall = ahf_ncall + 1
+         find3 = .true.
+         nhorizon = 3
       end if
 
 
+!     ***************************************************
+!     ***   START LOOP TO LOOK FOR SEVERAL HORIZONS   ***
+!     ***************************************************
+
+!     IMPORTANT:  This is the point where we start the loop
+!     to look for more than one horizon. Because of this,
+!     anything that is likely to be different for the different
+!     horizons should NOT be initialized until after this point!
+
+      do mfind=0,nhorizon-1
+
+!     Say what we are doing.
+
+      write(*,*)
+      write(*,*) 'AHFinder: Searching for horizon ',mfind
+
+
 !     ********************************
 !     ***   INITIALIZE  {status}   ***
 !     ********************************
@@ -562,21 +574,6 @@
       if (CCTK_Equals(ahf_octant,"high").eq.1) stepx=3
 
 
-!     ************************
-!     ***   FIND  {lmax}   ***
-!     ************************
-
-      lmax = ahf_lmax
-
-      if (lmax.ge.20) then
-         write(*,*)
-         write(*,*) 'ahf_lmax must be smaller than 20.'
-         write(*,*) 'STOPPED IN AHFinder.F'
-         write(*,*)
-         STOP
-      end if
-
-
 !     *******************************
 !     ***   FIND  {ntheta,nphi}   ***
 !     *******************************
@@ -593,27 +590,6 @@
       end if
 
 
-!     ***************************************************
-!     ***   ALLOCATE STORAGE FOR COEFFICIENT ARRAYS   ***
-!     ***************************************************
-
-      if ((ahf_ncall.eq.1).and.(mfind.eq.0)) then
-
-         allocate(c0(0:lmax),c0_old(0:lmax))
-         allocate(cc(lmax,lmax),cc_old(lmax,lmax))
-         allocate(cs(lmax,lmax),cs_old(lmax,lmax))
-
-         c0 = zero
-         cc = zero
-         cs = zero
-
-         c0_old = zero
-         cc_old = zero
-         cs_old = zero
-
-      end if
-
-
 !     ********************************************************
 !     ***   FIND TOTAL NUMBER OF TERMS IN EXPANSION {NN}   ***
 !     ********************************************************
@@ -1719,18 +1695,11 @@
       end if
 
 
-!     **********************************************
-!     ***   UPDATING COUNTER AND LOOP IF FIND3   ***
-!     **********************************************
+!     ********************************************************
+!     ***   END OF LOOP FOR LOOKING FOR SEVERAL HORIZONS   ***
+!     ********************************************************
 
-      if (find3) then
-         if (mfind.eq.2) then
-            mfind = 0
-         else
-            mfind = mfind + 1
-            goto 1
-         end if
-      end if
+      end do
 
 
 !     **************************************************
@@ -1742,16 +1711,9 @@
 
       if (ahfconformal) then
 
-         call CCTK_WARN(4,"Converting metric: physical -> conformal")
-
          conformal_state = CONFORMAL_METRIC
 
-         gxx = gxx/psi**4
-         gxy = gxy/psi**4
-         gxz = gxz/psi**4
-         gyy = gyy/psi**4
-         gyz = gyz/psi**4
-         gzz = gzz/psi**4
+         call PhysToConf(nx,ny,nz,psi,gxx,gxy,gxz,gyy,gyz,gzz)
 
       end if