cc/*@@
c  @file      ADMConstraints.F
c  @date      Aug 98
c  @desc 
c     Calculate the ADM Constraints for output:
c    
c     Hamiltonian Constraint is:
c
c       R - K^i_j K^j_i + trK^2 - 16 Pi rho
c
c     Momentum Constraints are:
c
c       Del_j K_i^j - Del_i trK - 8 Pi j_i
c
c  @enddesc 
c@@*/

c --------------------------------------------------------------

c CCTK Header files
#include "cctk.h"
#include "declare_parameters.h"
#include "declare_arguments.h"

c Using macro definitions from Einstein
#include "../../packages/CactusEinstein/Einstein/src/Einstein.h"

c --------------------------------------------------------------

c/*@@
c  @routine    ADMConstraints
c  @date       Aug 98
c  @desc 
c    Calculate the ADM Constraints
c  @enddesc 
c  @calls     
c  @calledby   
c  @history 
c
c  @endhistory 
c@@*/

      subroutine ADMConstraints(CCTK_FARGUMENTS)
      
      USE ADM_Scalars

      implicit none

c     Declare variables in the dynamic argument list
      DECLARE_CCTK_FARGUMENTS

c     Declare parameters
      DECLARE_PARAMETERS

c     Perhaps this and others should go into cctk.h
      INTEGER CCTK_Equals



      integer :: i,j,k,i1,i2,j1,j2,k1,k2
      REAL :: dx,dy,dz

      REAL :: pi,det,uxx,uxy,uxz,uyy,uyz,uzz

c     Temporaries for the Stress-Energy tensor
      REAL :: Ttt,Ttx,Tty,Ttz,Txx,Txy,Txz,Tyy,Tyz,Tzz


#include "CalcTmunu_temps.inc"

c     Macros from Standard Einstein
#include "../../packages/CactusEinstein/Einstein/src/macro/HAMADM_declare.h"
#include "../../packages/CactusEinstein/Einstein/src/macro/MOMXADM_declare.h"
#include "../../packages/CactusEinstein/Einstein/src/macro/MOMYADM_declare.h"
#include "../../packages/CactusEinstein/Einstein/src/macro/MOMZADM_declare.h"
#include "../../packages/CactusEinstein/Einstein/src/macro/DETG_declare.h"
#include "../../packages/CactusEinstein/Einstein/src/macro/UPPERMET_declare.h"

      dx = delta_space(1)
      dy = delta_space(2)
      dz = delta_space(3)

      i1 = 2
      i2 = sh(1)-1
      j1 = 2
      j2 = sh(2)-1
      k1 = 2
      k2 = sh(2)-1
  
      do i=i1,i2
        do j=j1,j2
          do k=k1,k2
            ham(i,j,k)  = 0.0D0
            momx(i,j,k) = 0.0D0
            momy(i,j,k) = 0.0D0
            momz(i,j,k) = 0.0D0
          end do
        end do
      end do

      Pi =  ACOS(-1D0)

      do k = k1, k2
         do j = j1, j2
            do i = k1, k2
               
c              Calculate the stress energy tensor at this point
c              ------------------------------------------------

               Ttt = 0.0D0; Ttx = 0.0D0; Tty = 0.0D0; Ttz = 0.0D0
               Txx = 0.0D0; Txy = 0.0D0; Txz = 0.0D0; Tyy = 0.0D0
               Tyz = 0.0D0; Tzz = 0.0D0
               
c              This may be needed for CalcTmunu etc
#include "../../packages/CactusEinstein/Einstein/src/macro/DETG_guts.h"
               det = DETG_DETCG
#include "../../packages/CactusEinstein/Einstein/src/macro/UPPERMET_guts.h"
               uxx = UPPERMET_UXX; uxy = UPPERMET_UXY; uxz = UPPERMET_UXZ
               uyy = UPPERMET_UYY; uyz = UPPERMET_UYZ; uzz = UPPERMET_UZZ
#include "CalcTmunu.inc"
              

c              Calculate the hamiltonian constraint
c              ------------------------------------

#include "../../packages/CactusEinstein/Einstein/src/macro/HAMADM_guts.h"

               ham(i,j,k) = HAMADM_HAMADM - 16D0*Pi*Ttt/alp(i,j,k)**2

#ifdef THORN_MAHC
               if ((some_hydro .eq. MAHCHYDRO).or.
     .              (some_hydro .eq. MAHCHYDRODUST)) then
                  ham(i,j,k) = ADM_ham(i,j,k) - 16.0d0 * pi *
     .                 ( (rho(i,j,k) + rho(i,j,k)*eps(i,j,k) + press(i,j,k)) *
     .                 w_lorentz(i,j,k)**2 - press(i,j,k))
               endif
#endif


c              Calculate the Momentum constraints
c              ----------------------------------

#include "../../packages/CactusEinstein/Einstein/src/macro/MOMXADM_guts.h"
#include "../../packages/CactusEinstein/Einstein/src/macro/MOMYADM_guts.h"
#include "../../packages/CactusEinstein/Einstein/src/macro/MOMZADM_guts.h"

               momx(i,j,k) = MOMXADM_MOMXADM - 8D0*Pi*Ttx/alp(i,j,k)**2
               momy(i,j,k) = MOMYADM_MOMYADM - 8D0*Pi*Tty/alp(i,j,k)**2
               momz(i,j,k) = MOMZADM_MOMZADM - 8D0*Pi*Ttz/alp(i,j,k)**2


#ifdef THORN_MAHC
               if ((some_hydro .eq. MAHCHYDRO).or.
     .              (some_hydro .eq. MAHCHYDRODUST)) then
                  momx(i,j,k) = ADM_momx(i,j,k) -
     .                 8.0d0 * pi * sx(i,j,k) / sqrt(det)
                  momy(i,j,k) = ADM_momy(i,j,k) -
     .                 8.0d0 * pi * sy(i,j,k) / sqrt(det)
                  momz(i,j,k) = ADM_momz(i,j,k) -
     .                 8.0d0 * pi * sz(i,j,k) / sqrt(det)
               endif
#endif
               
            enddo
         enddo
      enddo

#include "../../packages/CactusEinstein/Einstein/src/macro/DETG_undefine.h"
#include "../../packages/CactusEinstein/Einstein/src/macro/HAMADM_undefine.h"
#include "../../packages/CactusEinstein/Einstein/src/macro/MOMXADM_undefine.h"
#include "../../packages/CactusEinstein/Einstein/src/macro/MOMYADM_undefine.h"
#include "../../packages/CactusEinstein/Einstein/src/macro/MOMZADM_undefine.h"
#include "../../packages/CactusEinstein/Einstein/src/macro/UPPERMET_undefine.h"


c     Synchronize and apply boundaries conditions

c      call onebound(ARENA, ADM_ham,  1d0, 1d0, 1d0)
c      call onebound(ARENA, ADM_momx,-1d0, 1d0, 1d0)
c      call onebound(ARENA, ADM_momy, 1d0,-1d0, 1d0)
c      call onebound(ARENA, ADM_momz, 1d0,-1d0, 1d0)

      if (constraint_communication.eq.1) then
         call CCTK_SyncGroup(GH,"adm::ADM_constraints")
      end if

      end subroutine ADMConstraints