/*@@
  @file      ADM_Lapse.F
  @date      Aug 1998
  @desc 
     Set up lapse evolve for ADM routines
  
     We will evolve the equation 

       d_t alpha = ADM_gauge * K + flatNablaAlpCoeff D_flat^2 alpha
                 + nablaAlpCoeff D_g^2 alpha

  @enddesc 
@@*/

#include "cctk.h"
#include "cctk_Parameters.h"
#include "cctk_Arguments.h"


/*@@
  @routine    ADM_Lapse
  @date       Aug 1998
  @desc 
     Set up lapse evolve for ADM routines

     The dispersion coefficients are calculated once in ADM_Preloop.F
  @enddesc 
  @calls     
  @calledby  ADM_CCTK_REALLeap.F, ADM_StagLeap.F 
  @history 
   @hauthor Gerd Lanfermann
   @hdate July 99
   @hdesc
     For the slicing, the routine have been changed to use the slicing handles.
     Read ./Einstein/Slicing.c for the details. There is a grid scalar
     called "active_slicing_handle" which defines what slicing to use in the
     next iteration(it's a number). Every slicing checks against this number 
     by calling GetSlicingHandle, if it matches, it does the slicing.

     This basically substituted the prev. CCTK_Equals(slicing,"  ") comparison.
  @endhistory 
@@*/

      subroutine ADM_Lapse(CCTK_ARGUMENTS)

      USE ADM_Scalars

      implicit none

c     Declare variables in the dynamic argument list
      DECLARE_CCTK_ARGUMENTS

c     Declare parameters
      DECLARE_CCTK_PARAMETERS

c     Declare CCTK functions and others used here
      DECLARE_CCTK_FUNCTIONS
      INTEGER Einstein_GetSlicingHandle

c     Diffusion coefficients for lapse evolution
      CCTK_REAL :: alp_Geps,alp_nablaeps

      CHARACTER*200 :: fortran_string
      CCTK_INT      :: fortran_string_len


c     GEODESIC SLICING
c     ----------------

      if (active_slicing_handle .eq. Einstein_GetSlicingHandle("geodesic"))  then

	 alp = 1D0
         ADM_gauge = 10000000.0D0

c     STATIC SLICING
c     --------------

      else if (active_slicing_handle .eq. Einstein_GetSlicingHandle("static")) then

         ADM_gauge = 10000000.0D0

c     1+LOG SLICING
c     -------------

      else if (active_slicing_handle .eq. Einstein_GetSlicingHandle("1+log")) then

         alp_Geps     = ADM_flatNablaAlpCoeff
         alp_nablaeps = ADM_nablaAlpCoeff

c        In the NCSA style

         if (CCTK_EQUALS(slicing_flavour,'kleing')) then

            if (store_initial_lapse == 1) then
               ADM_gauge = -2d0*alp*ADM_alp0
            else
               ADM_gauge = -2D0*alp
	    end if

         else

            ADM_gauge = -harmonic_f*alp

         endif

         if (alp_nablaeps .lt. 1e-10) then
            alp_nablaeps = 0D0  ! This is for performance
         end if


c     BONA-MASSO TYPE HARMONIC SLICINGS
c     ---------------------------------

      else if (CCTK_EQUALS(slicing,'harmonic')) then

       
         ADM_gauge = -alp**2*harmonic_f


c     NCSA TYPE SLICING WITH g(gamma)=sqrt(det g) 
c     -----------------
c     This is just plain "harmonic"

      else if (active_slicing_handle .eq. Einstein_GetSlicingHandle("detg")) then

         ADM_gauge = -alp**2


c     EXTERNAL SLICING

      else if (active_slicing_handle .eq. Einstein_GetSlicingHandle("external")) then

         ADM_gauge = 100000000.0

c     MAXIMAL SLICING

      else if (active_slicing_handle .eq. Einstein_GetSlicingHandle("maximal")) then

         ADM_gauge = 100000000.0

c     and the final sanity check if whatever is left is registered at all!
      else
c        Get a fortran string for the slicing
         call CCTK_FortranString(fortran_string_len,slicing,fortran_string)
         if (Einstein_GetSlicingHandle(fortran_string(1:fortran_string_len)).lt.0) then
            call CCTK_WARN(1,"thorn ADM: slicing in parfile not registered by ADM")
         else
            call CCTK_WARN(1,"thorn ADM: slicing registered but not coded!") 
         end if
      end if

      return
      end



      subroutine ADM_alp0_init(CCTK_ARGUMENTS)

      implicit none

      DECLARE_CCTK_ARGUMENTS
      DECLARE_CCTK_PARAMETERS

      ADM_alp0 = alp

      return
      end