source_komar_mass_compact_WL.f90

Go to the documentation of this file.

subroutine source_komar_mass_compact_WL(souf)
  use phys_constant, only  : long, pi
  use grid_parameter, only : nrf, ntf, npf
  use coordinate_grav_r, only       : rg
  use trigonometry_grav_theta, only : sinthg, costhg
  use trigonometry_grav_phi, only   : sinphig, cosphig
  use def_matter
  use def_matter_parameter, only  :   radi, ber, ome
  use def_metric, only : psi, alph, bvxd, bvyd, bvzd, bvxu, bvyu, bvzu
  use def_metric_hij, only : hxxd, hxyd, hxzd, hyyd, hyzd, hzzd
  use def_vector_phi, only : vec_phif
  use make_array_3d
  use interface_interpo_gr2fl
  implicit none
  real(long), pointer :: souf(:,:,:)
  real(long), pointer :: alphf(:,:,:),  psif(:,:,:) 
  real(long), pointer :: bvxdf(:,:,:), bvydf(:,:,:), bvzdf(:,:,:)
  real(long), pointer :: bvxuf(:,:,:), bvyuf(:,:,:), bvzuf(:,:,:)
  real(long), pointer :: hxxdf(:,:,:), hxydf(:,:,:), hxzdf(:,:,:), 
                        hyydf(:,:,:), hyzdf(:,:,:), hzzdf(:,:,:)
  real(long) :: emdw, alphw, psiw, rhow, prew, hhw, utw, rhoHw, esseS
  real(long) :: rjjx, rjjy, rjjz, rjjbeta, ene
  real(long) :: vphif(1:3)
  real(long) :: otermx, otermy, otermz
  real(long) :: bvxdfw, bvydfw, bvzdfw, bvxufw, bvyufw, bvzufw
  integer    :: irf, itf, ipf
  real(long) :: zfac, small = 1.0d-15
  real(long) :: hhxxdf, hhxydf, hhxzdf, hhyxdf, hhyydf, hhyzdf, 
               hhzxdf, hhzydf, hhzzdf
!
  call alloc_array3d(psif, 0, nrf, 0, ntf, 0, npf)
  call alloc_array3d(alphf, 0, nrf, 0, ntf, 0, npf)
  call alloc_array3d(bvxdf, 0, nrf, 0, ntf, 0, npf)
  call alloc_array3d(bvydf, 0, nrf, 0, ntf, 0, npf)
  call alloc_array3d(bvzdf, 0, nrf, 0, ntf, 0, npf)
  call alloc_array3d(bvxuf, 0, nrf, 0, ntf, 0, npf)
  call alloc_array3d(bvyuf, 0, nrf, 0, ntf, 0, npf)
  call alloc_array3d(bvzuf, 0, nrf, 0, ntf, 0, npf)
!
  call alloc_array3d(hxxdf, 0, nrf, 0, ntf, 0, npf)
  call alloc_array3d(hxydf, 0, nrf, 0, ntf, 0, npf)
  call alloc_array3d(hxzdf, 0, nrf, 0, ntf, 0, npf)
  call alloc_array3d(hyydf, 0, nrf, 0, ntf, 0, npf)
  call alloc_array3d(hyzdf, 0, nrf, 0, ntf, 0, npf)
  call alloc_array3d(hzzdf, 0, nrf, 0, ntf, 0, npf)
!
  call interpo_gr2fl(alph, alphf)
  call interpo_gr2fl(psi, psif)
  call interpo_gr2fl(bvxd, bvxdf)
  call interpo_gr2fl(bvyd, bvydf)
  call interpo_gr2fl(bvzd, bvzdf)
  call interpo_gr2fl(bvxu, bvxuf)
  call interpo_gr2fl(bvyu, bvyuf)
  call interpo_gr2fl(bvzu, bvzuf)
!
  call interpo_gr2fl(hxxd, hxxdf)
  call interpo_gr2fl(hxyd, hxydf)
  call interpo_gr2fl(hxzd, hxzdf)
  call interpo_gr2fl(hyyd, hyydf)
  call interpo_gr2fl(hyzd, hyzdf)
  call interpo_gr2fl(hzzd, hzzdf)
!
  do ipf = 0, npf
    do itf = 0, ntf
      do irf = 0, nrf
        emdw = emd(irf,itf,ipf)
        if (emdw <= small) emdw = small
        psiw = psif(irf,itf,ipf)
        alphw = alphf(irf,itf,ipf)
        bvxufw = bvxuf(irf,itf,ipf)
        bvyufw = bvyuf(irf,itf,ipf)
        bvzufw = bvzuf(irf,itf,ipf)
        bvxdfw = bvxdf(irf,itf,ipf)
        bvydfw = bvydf(irf,itf,ipf)
        bvzdfw = bvzdf(irf,itf,ipf)
        call peos_q2hprho(emdw, hhw, prew, rhow, ene)
        utw = hhw/ber
!
        vphif(1) = vec_phif(irf,itf,ipf,1)
        vphif(2) = vec_phif(irf,itf,ipf,2)
        vphif(3) = vec_phif(irf,itf,ipf,3)
!
        hhxxdf = hxxdf(irf,itf,ipf)
        hhxydf = hxydf(irf,itf,ipf)
        hhxzdf = hxzdf(irf,itf,ipf)
        hhyydf = hyydf(irf,itf,ipf)
        hhyzdf = hyzdf(irf,itf,ipf)
        hhzzdf = hzzdf(irf,itf,ipf)
        hhyxdf = hhxydf
        hhzxdf = hhxzdf
        hhzydf = hhyzdf
!        
        otermx = bvxdfw + ome*vphif(1) &
             &   + hhxxdf*ome*vphif(1) &
             &   + hhxydf*ome*vphif(2) &
             &   + hhxzdf*ome*vphif(3)
        otermy = bvydfw + ome*vphif(2) &
             &   + hhyxdf*ome*vphif(1) &
             &   + hhyydf*ome*vphif(2) &
             &   + hhyzdf*ome*vphif(3)
        otermz = bvzdfw + ome*vphif(3) &
             &   + hhzxdf*ome*vphif(1) &
             &   + hhzydf*ome*vphif(2) &
             &   + hhzzdf*ome*vphif(3)
!
        rhoHw = hhw*rhow*(alphw*utw)**2 - prew
        esseS = -hhw*rhow + 4.0d0*prew + rhoHw
! 
        rjjx = hhw*rhow*alphw*utw**2*psiw**4*otermx
        rjjy = hhw*rhow*alphw*utw**2*psiw**4*otermy
        rjjz = hhw*rhow*alphw*utw**2*psiw**4*otermz
!
        rjjbeta = rjjx*bvxufw + rjjy*bvyufw + rjjz*bvzufw
!
        souf(irf,itf,ipf) = (alphw*(esseS+rhoHw) - 2.0d0*rjjbeta)*psiw**6
      end do
    end do
  end do
!
  deallocate(psif)
  deallocate(alphf)
  deallocate(bvxdf)
  deallocate(bvydf)
  deallocate(bvzdf)
  deallocate(bvxuf)
  deallocate(bvyuf)
  deallocate(bvzuf)
!
  deallocate(hxxdf)
  deallocate(hxydf)
  deallocate(hxzdf)
  deallocate(hyydf)
  deallocate(hyzdf)
  deallocate(hzzdf)
!
end subroutine source_komar_mass_compact_WL