!*********************************************************************** 
!  Computes the residual of the hamiltonian and momentum constraint 
!  equations. 
!*********************************************************************** 
!  The following routine calls the hamconstraint and momconstraint routines 
!  and evaluates their norms. It is possible to dump the residual of the  
!  constraints as a diagnostic tool  at this point.  
!----------------------------------------------------------------------- 
! hamconstraint : MFH 1996 
! 
!   Routine to compute the residual of the Hamiltonian constraint 
!   Norms will be computed outside this routine. In the DAGH case  
!   the norm is calculated using the DAGH Norm routine. 
! 
subroutine h_constraint(g,k,r,gup,t,mask,hc)
implicit none
real*8,dimension(:,:,:),intent(in)  :: mask
real*8,dimension(:,:,:),intent(out) :: hc
real*8,dimension(:,:,:,:),target,intent(in)    :: g,gup,k,r
real*8,dimension(:,:,:,:),target,intent(out)   :: t
interface
   subroutine trace02sym(gup11, gup12, gup13, gup22, gup23, gup33,       &
                         t11, t12, t13, t22, t23, t33, trT)
      implicit none
      real*8,dimension(:,:,:)    :: gup11,gup12,gup13,gup22,gup23,gup33
      real*8,dimension(:,:,:)    :: t11,t12,t13,t22,t23,t33
      real*8,dimension(:,:,:)    :: trT
   end subroutine

   subroutine raise02to20s(gup11, gup12, gup13, gup22, gup23, gup33,        &
                           tin11, tin12, tin13, tin22, tin23, tin33,        &
                           tout11, tout12, tout13, tout22, tout23, tout33)
      implicit none
      real*8                        :: h, dt
      real*8,dimension(:,:,:)       :: gup11,gup12,gup13,gup22,gup23,gup33
      real*8,dimension(:,:,:)       :: tin11,tin12,tin13,tin22,tin23,tin33
      real*8,dimension(:,:,:)       :: tout11,tout12,tout13,tout22,tout23,tout33
	end subroutine

end interface

!      Compute the Ricci scalar and store in hc 
 
   call trace02sym(gup(1,:,:,:),gup(2,:,:,:),gup(3,:,:,:),gup(4,:,:,:),   &
                   gup(5,:,:,:),gup(6,:,:,:),r(1,:,:,:),r(2,:,:,:),       &
                   r(3,:,:,:),r(4,:,:,:),r(5,:,:,:),r(6,:,:,:),hc) 
 
!      Compute the trace of the extrinsic curvature and store in t temporarily 
   call trace02sym(gup(1,:,:,:),gup(2,:,:,:),gup(3,:,:,:),gup(4,:,:,:),   &
                   gup(5,:,:,:),gup(6,:,:,:),k(1,:,:,:),k(2,:,:,:),       &
                   k(3,:,:,:),k(4,:,:,:),k(5,:,:,:),k(6,:,:,:),t(1,:,:,:)) 
 
!      Add the trace of K squared to hamc 
   hc  = ( hc  + t(1,:,:,:) *t(1,:,:,:)  )  * mask  
 
!     Raise the components of the extrinsic curvature tensor 
 
   call raise02to20s(gup(1,:,:,:),gup(2,:,:,:),gup(3,:,:,:),gup(4,:,:,:),   &
                   gup(5,:,:,:),gup(6,:,:,:),k(1,:,:,:),k(2,:,:,:),         &
                   k(3,:,:,:),k(4,:,:,:),k(5,:,:,:),k(6,:,:,:),             &
                   t(1,:,:,:),t(2,:,:,:),t(3,:,:,:),t(4,:,:,:),             &
                   t(5,:,:,:),t(6,:,:,:) )

!      Contraction to get  - K^{ij}  K_{ij} : subtract from hc 
 
   hc  = hc  - ( t(1,:,:,:)*k(1,:,:,:)+t(2,:,:,:)*k(2,:,:,:)*2.0d0 +        &
                 t(3,:,:,:)*k(3,:,:,:)*2.0d0+t(4,:,:,:)*k(4,:,:,:) +        &
                 t(5,:,:,:)*k(5,:,:,:)*2.0d0+t(6,:,:,:)*k(6,:,:,:) ) * mask
 
return 
end 
!----------------------------------------------------------------------- 
! momconstraint: MFH 1996 
! 
!   Routine to compute the momentum constraint residuals.  
!   Norms will be computed outside this routine. In the DAGH case  
!   the norm is calculated using the DAGH Norm routine. 
subroutine k_constraint(g,k,gup,c,h,t,mask,m,d,nl,kc)
implicit none
integer                             :: nl
real*8,dimension(:,:,:)             :: mask
real*8,dimension(:,:,:,:)           :: g,gup,k,t,c,kc
real*8,dimension(:)              :: d
integer,dimension(:)             :: m
real*8                              :: h

interface

   subroutine raise02to20s(gup11, gup12, gup13, gup22, gup23, gup33,        &
                           tin11, tin12, tin13, tin22, tin23, tin33,        &
                           tout11, tout12, tout13, tout22, tout23, tout33)
      implicit none
      real*8                        :: h, dt
      real*8,dimension(:,:,:)       :: gup11,gup12,gup13,gup22,gup23,gup33
      real*8,dimension(:,:,:)       :: tin11,tin12,tin13,tin22,tin23,tin33
      real*8,dimension(:,:,:)       :: tout11,tout12,tout13,tout22,tout23,tout33
   end subroutine

   subroutine diff(h,f,df,dir,m,d,nl)
      implicit none
      integer                       :: nl
      real*8                        :: h
      real*8,dimension(:,:,:)    :: f, df
      real*8,dimension(:)        :: d
      integer,dimension(:)       :: m
      character                     :: dir
   end subroutine



end interface

!     Raise the components of the extrinsic curvature tensor 
 
   call raise02to20s(gup(1,:,:,:),gup(2,:,:,:),gup(3,:,:,:),gup(4,:,:,:),   &
                   gup(5,:,:,:),gup(6,:,:,:),k(1,:,:,:),k(2,:,:,:),         &
                   k(3,:,:,:),k(4,:,:,:),k(5,:,:,:),k(6,:,:,:),             &
                   t(1,:,:,:),t(2,:,:,:),t(3,:,:,:),t(4,:,:,:),             &
                   t(5,:,:,:),t(6,:,:,:) )


 
!      Compute the connection terms in the covariant derivative term 
 
!        c^{i}_{jl} K^{lj} terms 
 
   kc(1,:,:,:)=(c(1,:,:,:)*t(1,:,:,:)+c(2,:,:,:)*t(2,:,:,:)*2.0d0 +       &
                c(3,:,:,:)*t(3,:,:,:)*2.0d0+c(4,:,:,:)*t(4,:,:,:) +       &
                c(5,:,:,:)*t(5,:,:,:)*2.0d0+c(6,:,:,:)*t(6,:,:,:) )*mask

   kc(2,:,:,:)=(c(7,:,:,:)*t(1,:,:,:)+c(8,:,:,:)*t(2,:,:,:)*2.0d0 +       &
                c(9,:,:,:)*t(3,:,:,:)*2.0d0+c(10,:,:,:)*t(4,:,:,:)+      &
                c(11,:,:,:)*t(5,:,:,:)*2.0d0+c(12,:,:,:)*t(6,:,:,:) )*mask

   kc(3,:,:,:)=(c(13,:,:,:)*t(1,:,:,:)+c(14,:,:,:)*t(2,:,:,:)*2.0d0 +   &
                c(15,:,:,:)*t(3,:,:,:)*2.0d0+c(16,:,:,:)*t(4,:,:,:)+      &
                c(17,:,:,:)*t(5,:,:,:)*2.0d0+c(18,:,:,:)*t(6,:,:,:) )*mask
 
!        c^{j}_{jl} K^{il} terms 

   kc(1,:,:,:)=kc(1,:,:,:)+mask*((c(1,:,:,:)+c(8,:,:,:)+c(15,:,:,:))*      &
               t(1,:,:,:)+(c(2,:,:,:)+c(10,:,:,:)+c(17,:,:,:))*t(2,:,:,:)+ &
               (c(3,:,:,:)+c(11,:,:,:)+c(18,:,:,:))*t(3,:,:,:))

   kc(2,:,:,:)= kc(2,:,:,:)+mask*((c(1,:,:,:)+c(8,:,:,:)+c(15,:,:,:))*     &
                t(2,:,:,:)+(c(2,:,:,:)+c(10,:,:,:)+c(17,:,:,:))*t(4,:,:,:) &
                +(c(3,:,:,:)+c(11,:,:,:)+c(18,:,:,:))*t(5,:,:,:)) 
! 
   kc(3,:,:,:)= kc(3,:,:,:)+mask*((c(1,:,:,:)+c(8,:,:,:)+c(15,:,:,:))*     &
                t(3,:,:,:)+(c(2,:,:,:)+c(10,:,:,:)+c(17,:,:,:))*t(5,:,:,:) &
                +(c(3,:,:,:)+c(11,:,:,:)+c(18,:,:,:))*t(8,:,:,:)) 
! 
! 
!      Partial derivatives : K^{ij},j 

   call diff(h,t(1,:,:,:),t(7,:,:,:),'x',m,d,nl) 
   kc(1,:,:,:) = kc(1,:,:,:) + t(7,:,:,:)
   call diff(h,t(2,:,:,:),t(7,:,:,:),'y',m,d,nl) 
   kc(1,:,:,:) = kc(1,:,:,:) + t(7,:,:,:)
   call diff(h,t(3,:,:,:),t(7,:,:,:),'z',m,d,nl) 
   kc(1,:,:,:) = kc(1,:,:,:) + t(7,:,:,:)
 
   call diff(h,t(2,:,:,:),t(7,:,:,:),'x',m,d,nl) 
   kc(2,:,:,:) = kc(2,:,:,:) + t(7,:,:,:)
   call diff(h,t(4,:,:,:),t(7,:,:,:),'y',m,d,nl) 
   kc(2,:,:,:) = kc(2,:,:,:) + t(7,:,:,:)
   call diff(h,t(5,:,:,:),t(7,:,:,:),'z',m,d,nl) 
   kc(2,:,:,:) = kc(2,:,:,:) + t(7,:,:,:)
 
   call diff(h,t(3,:,:,:),t(7,:,:,:),'x',m,d,nl) 
   kc(3,:,:,:) = kc(3,:,:,:) + t(7,:,:,:)
   call diff(h,t(5,:,:,:),t(7,:,:,:),'y',m,d,nl) 
   kc(3,:,:,:) = kc(3,:,:,:) + t(7,:,:,:)
   call diff(h,t(6,:,:,:),t(7,:,:,:),'z',m,d,nl) 
   kc(3,:,:,:) = kc(3,:,:,:) + t(7,:,:,:)
 
return 
end  
!----------------------------------------------------------------------- 
