c/*@@
c  @file      AHFinder_initguess.F
c  @date      July 1999
c  @author    Miguel Alcubierre
c  @desc 
c             Find initial guess for minimization algorithm
c  @enddesc 
c@@*/

c Note that including cactus.h will also include AHFinder.h
!#include "cactus.h"
#include "cctk.h"
#include "cctk_parameters.h"
#include "cctk_arguments.h"

      subroutine AHFinder_initguess(CCTK_FARGUMENTS,rmx)

      use AHFinder_dat

      implicit none

      DECLARE_CCTK_FARGUMENTS
      DECLARE_CCTK_PARAMETERS

      integer i,j
      CCTK_INT nn0,nn2,min_i,min_j

      CCTK_REAL zero,one,cc0,cc2,rmx,aux,min_aa
      CCTK_REAL, allocatable, dimension(:,:) :: aa,a0,a2
      
!     Description of variables
!
!     i,j       Counters.     
!
!     aa        Area array for initial guess.
!     a0        c0(0) array for initial guess.
!     a2        c0(2) array for initial guess.
!
!     rmx       Radius of largest sphere that fits in the grid.
!
!     nn0       Number of subdivisions of c0(0) for initial guess.
!     nn2       Number of subdivisions of c0(2) for initial guess.



!     **********************
!     ***   INITIALIZE   ***
!     **********************

!     define numbers

      zero = 0.0D0
      one  = 1.0D0
      aux = rmx

!     Initialize arrays.

      c0 = zero
      cc = zero
      cs = zero
      cc0 = zero
      cc2 = zero

!     Get {nn0,nn2}.

!      nn0 = geti("ahf_nn0")
!      nn2 = geti("ahf_nn2")
      nn0 = ahf_nn0
      nn2 = ahf_nn0

      if (cartoon) nn2 = 0

!     Allocate storage.

      allocate(aa(0:nn0,0:nn2))
      allocate(a0(0:nn0,0:nn2))
      allocate(a2(0:nn0,0:nn2))

      aa = zero
      a0 = zero
      a2 = zero

!     ****************************
!     ***  FIND INITIAL GUESS  ***
!     ****************************

!     IMPORTANT:  In order to choose a good initial guess,
!     I map a section of the space {c0(0),c0(2)} looking
!     for a local minimum.  This is VERY expensive, so the
!     parameter space is not very well convered.

      if (veryver.and.(myproc.eq.0)) then
         write(*,*)
         write(*,*) 'AHFinder: Choosing initial guess.'
         write(*,*) 'Be patient, this will take a while.'
      end if


!     lmax >= 2 and good guess.

      if ((lmax.ge.2).and.(.not.sloppy)) then

         do j=0,nn2

            c0(2) = (0.5D0*dble(j)/dble(nn2) - 0.25D0)*aux

            call AHFinder_fun(CCTK_FARGUMENTS)
            call AHFinder_exp(CCTK_FARGUMENTS)

            do i=0,nn0

               c0(0) = 4.0D0*dx + dble(i)/dble(nn0)*(aux - 4.0D0*dx)

               a0(i,j) = c0(0)
               a2(i,j) = c0(2)

               call AHFinder_int(CCTK_FARGUMENTS)

               if (minarea) then
                  aa(i,j) = intarea
               else
                  aa(i,j) = intexp2
               end if
               if(guessverbose) then
                  write(*,"(A11,2ES14.6)") 'a0 and a2: ',a0(i,j),a2(i,j)
                  write(*,"(A8,ES14.6)") '   H  : ',intexp
                  write(*,"(A8,ES14.6)") '   H^2: ',intexp2
                  write(*,"(A13,ES14.6)") '   H/area  : ',intexp/intarea
                  write(*,"(A13,ES14.6)") '   H^2/area: ',intexp2/intarea
               end if

            end do
         end do

         if ((.not.inner).or.(minarea)) then

            min_i = 1
            min_j = 1
            min_aa = 1.0d15
            do i=1,nn0-1
               do j=1,nn2-1
                  if ((aa(i,j).gt.zero).and.
     .                (aa(i,j).lt.aa(i+1,j)).and.
     .                (aa(i,j).lt.aa(i-1,j)).and.
     .                (aa(i,j).lt.aa(i,j+1)).and.
     .                (aa(i,j).lt.aa(i,j-1))) then
                     cc0 = a0(i,j)
                     cc2 = a2(i,j)
                  endif
                  if (aa(i,j) .le. min_aa) then
                     min_i = i
                     min_j = j
                     min_aa = aa(i,j)
                  endif
               end do
            end do

         else

            min_i = nn0-1
            min_j = 1
            min_aa = 1.0d15
            do i=nn0-1,1,-1
               do j=1,nn2-1
                  if ((aa(i,j).gt.zero).and.
     .                (aa(i,j).lt.aa(i+1,j)).and.
     .                (aa(i,j).lt.aa(i-1,j)).and.
     .                (aa(i,j).lt.aa(i,j+1)).and.
     .                (aa(i,j).lt.aa(i,j-1))) then
                     cc0 = a0(i,j)
                     cc2 = a2(i,j)
                  endif
                  if (aa(i,j) .le. min_aa) then
                     min_i = i
                     min_j = j
                     min_aa = aa(i,j)
                  endif
               end do
            end do

         end if
         if (guess_absmin) then
            cc0 = a0(min_i,min_j)
            cc2 = a2(min_i,min_j)
         endif

!        Set up initial guess.

         c0(0) = cc0
         c0(2) = cc2

!     lmax < 2 or sloppy guess.

      else

         call AHFinder_fun(CCTK_FARGUMENTS)
         call AHFinder_exp(CCTK_FARGUMENTS)

         do i=0,nn0

            c0(0) = 4.0D0*dx + dble(i)/dble(nn0)*(aux - 4.0D0*dx)

            a0(i,0) = c0(0)

            call AHFinder_int(CCTK_FARGUMENTS)

            if (minarea) then
               aa(i,0) = intarea
            else
               aa(i,0) = intexp2
            end if
            if(guessverbose) then
               write(*,"(A11,2ES14.6)") 'a0 and a2: ',a0(i,0),a2(i,0)
               write(*,"(A8,ES14.6)") '   H  : ',intexp
               write(*,"(A8,ES14.6)") '   H^2: ',intexp2
               write(*,"(A13,ES14.6)") '   H/area  : ',intexp/intarea
               write(*,"(A13,ES14.6)") '   H^2/area: ',intexp2/intarea
            end if

         end do

         if ((.not.inner).or.(minarea)) then

            do i=1,nn0-1
               if ((aa(i,0).gt.zero).and.
     .             (aa(i,0).lt.0.9D0*aa(i+1,0)).and.
     .             (aa(i,0).lt.0.9D0*aa(i-1,0))) then
                  cc0 = a0(i,0)
               end if
            end do

         else

            do i=nn0-1,1,-1
               if ((aa(i,0).gt.zero).and.
     .             (aa(i,0).lt.0.9D0*aa(i+1,0)).and.
     .             (aa(i,0).lt.0.9D0*aa(i-1,0))) then
                  cc0 = a0(i,0)
               end if
            end do

         end if

!        Set up initial guess.

         c0(0) = cc0

      end if

!     If no local minimum was found, start with a large
!     sphere anyway.

      if (c0(0).eq.zero) then
         c0(0) = 0.8D0*rmx
         c0(2) = zero
      end if


!     ***************
!     ***   END   ***
!     ***************

!     Print message.

      if (veryver.and.(myproc.eq.0)) then
         write(*,*)
         write(*,*) 'Initial guess found, pinning down the horizon.'
      end if

!     Deallocate storage.

      deallocate(aa,a0,a2)


      end subroutine AHFinder_initguess