/*@@
   @file      GRHydro_Tmunu.F90
   @date      Aug 30, 2010
   @author    Joshua Faber, Scott Noble, Bruno Mundim, Ian Hawke
   @histpry
    Apr. 2009: Luca Baiotti copied and adapted for the Tmunu-thorn mechanism the original include file
   @desc 
     The calculation of the stress energy tensor.
     The version used here was worked out by Miguel Alcubierre. I
     think it was an extension of the routine from GR3D, written
     by Mark Miller.
     C version added by Ian Hawke.

     Lower components of the stress-energy tensor obtained from
     the hydro variables.  The components are given by:

     T      =  (rho h +b^2) u   u    +  (P+b^2/2)  g      - b   b
      mu nu                   mu  nu                mu nu     mu nu 

     where rho is the energy density of the fluid, h the enthalpy
     and P the pressure.  The enthalpy is given in terms of the
     basic variables as:

     h  =  1  +  e  +  P/rho

     with e the internal energy (eps here).

     In the expresion for T_{mu,nu} we also have the four-velocity
     of the fluid given by (v_i is the 3-velocity field):

                                 i
     u  =  W ( - alpha +  v  beta  )
      0                    i

     u  =  W v
      i       i
                                                i  -1/2
     with W the Lorentz factor:   W = ( 1 -  v v  )
                                              i

     and where alpha and beta are the lapse and shift vector.

     Finally, the 4 metric is given by

                   2             i
     g   =  - alpha  + beta  beta
      00                   i

     g   =  beta
      0i        i


     g   =  gamma      (the spatial metric)
      ij        ij


   @enddesc 
 @@*/
#include "cctk.h"
#include "cctk_Arguments.h"
#include "cctk_Parameters.h"
#include "SpaceMask.h"

#define velx(i,j,k) vel(i,j,k,1)
#define vely(i,j,k) vel(i,j,k,2)
#define velz(i,j,k) vel(i,j,k,3)
#define Bvecx(i,j,k) Bvec(i,j,k,1)
#define Bvecy(i,j,k) Bvec(i,j,k,2)
#define Bvecz(i,j,k) Bvec(i,j,k,3)

 subroutine GRHydro_TmunuM(CCTK_ARGUMENTS)

   implicit none

   DECLARE_CCTK_ARGUMENTS

   CCTK_REAL :: velxlow, velylow, velzlow
   CCTK_REAL :: betaxlow, betaylow, betazlow, beta2
   CCTK_REAL :: Bvecxlow,Bvecylow,Bveczlow
   CCTK_REAL :: bdotv,b2,bxlow,bylow,bzlow,btlow,dum
   CCTK_REAL :: utlow,rhohstarw2,pstar
   CCTK_REAL :: bdotbeta,vdotbeta
   CCTK_INT :: i,j,k


   !$OMP PARALLEL DO PRIVATE(i,j,velxlow, velylow, velzlow,&
   !$OMP Bvecxlow,Bvecylow,Bveczlow, bdotv,dum,b2,bxlow,bylow,bzlow,&
   !$OMP betaxlow, betaylow, betazlow, beta2, bdotbeta,vdotbeta,utlow, btlow,&
   !$OMP rhohstarw2,pstar)

   do k = 1, cctk_lsh(3)
     do j = 1, cctk_lsh(2)
       do i = 1, cctk_lsh(1)

          call calc_vlow_blow(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),&
               gyy(i,j,k),gyz(i,j,k),gzz(i,j,k), &
               velx(i,j,k),vely(i,j,k),velz(i,j,k),Bvecx(i,j,k),Bvecy(i,j,k),Bvecz(i,j,k), &
               velxlow,velylow,velzlow,Bvecxlow,Bvecylow,Bveczlow, &
               bdotv,dum,dum,b2,bxlow,bylow,bzlow)
          
!!$       Calculate lower components and square of shift vector.
          
          if (shift_state .ne. 0) then
             
             betaxlow = gxx(i,j,k)*betax(i,j,k) + gxy(i,j,k)*betay(i,j,k) + gxz(i,j,k)*betaz(i,j,k)
             betaylow = gxy(i,j,k)*betax(i,j,k) + gyy(i,j,k)*betay(i,j,k) + gyz(i,j,k)*betaz(i,j,k)
             betazlow = gxz(i,j,k)*betax(i,j,k) + gyz(i,j,k)*betay(i,j,k) + gzz(i,j,k)*betaz(i,j,k)
             beta2 = betax(i,j,k)*betaxlow + betay(i,j,k)*betaylow + betaz(i,j,k)*betazlow 
            
             bdotbeta = betaxlow*Bvecx(i,j,k)+betaylow*Bvecy(i,j,k)+betazlow*Bvecz(i,j,k)
             vdotbeta = betaxlow*velx(i,j,k)+betaylow*vely(i,j,k)+betazlow*velz(i,j,k)

!!$   u0 low is missing the w_lorentz factor (see below)!!
             utlow = -1.d0*alp(i,j,k) + vdotbeta

             btlow = -1.0d0*w_lorentz(i,j,k)*alp(i,j,k)*bdotv + &
                  bdotbeta/w_lorentz(i,j,k) + w_lorentz(i,j,k)*bdotv*vdotbeta

            
         else
            
            betaxlow = 0.0D0
            betaylow = 0.0D0
            betazlow = 0.0D0            
            beta2 = 0.0D0

!!$   u0 low is missing the w_lorentz factor (see below)!!
            utlow = -1.0*alp(i,j,k)
            btlow = utlow*w_lorentz(i,j,k)*bdotv
            
         end if

!!$       Calculate the specific relativistic enthalpy times rho + the mag. field contribution times the
!!$       square of the lorentz factor.

         rhohstarw2 = w_lorentz(i,j,k)**2*(rho(i,j,k)*(1.0D0 + eps(i,j,k)) + press(i,j,k) + b2)
         pstar =  press(i,j,k)+b2/2.d0

!!$       Calculate lower components of 4-velocity (without the Lorent factor).
!!$         uxlow = velxlow
!!$         uylow = velylow
!!$         uzlow = velzlow

!!$       Calculate Tmunu (the lower components!).

         eTtt(i,j,k) = eTtt(i,j,k) + rhohstarw2*utlow**2 + pstar*(beta2 - alp(i,j,k)**2) - btlow**2

         eTtx(i,j,k) = eTtx(i,j,k) + rhohstarw2*utlow*velxlow + pstar*betaxlow - btlow*bxlow
         eTty(i,j,k) = eTty(i,j,k) + rhohstarw2*utlow*velylow + pstar*betaylow - btlow*bylow
         eTtz(i,j,k) = eTtz(i,j,k) + rhohstarw2*utlow*velzlow + pstar*betazlow - btlow*bzlow

         eTxx(i,j,k) = eTxx(i,j,k) + rhohstarw2*velxlow**2 + pstar*gxx(i,j,k) - bxlow**2
         eTyy(i,j,k) = eTyy(i,j,k) + rhohstarw2*velylow**2 + pstar*gyy(i,j,k) - bylow**2
         eTzz(i,j,k) = eTzz(i,j,k) + rhohstarw2*velzlow**2 + pstar*gzz(i,j,k) - bzlow**2
         
         eTxy(i,j,k) = eTxy(i,j,k) + rhohstarw2*velxlow*velylow + pstar*gxy(i,j,k) - bxlow*bylow
         eTxz(i,j,k) = eTxz(i,j,k) + rhohstarw2*velxlow*velzlow + pstar*gxz(i,j,k) - bxlow*bzlow
         eTyz(i,j,k) = eTyz(i,j,k) + rhohstarw2*velylow*velzlow + pstar*gyz(i,j,k) - bylow*bzlow

       end do
     end do
   end do
   !$OMP END PARALLEL DO

   return
 
 end subroutine GRHydro_TmunuM