/home/uryu/Cocal_svn_checkinout/Cocal/code/Subroutine/sourceterm_MoC_CF_with_divshift.f90

Go to the documentation of this file.

subroutine sourceterm_MoC_CF_with_divshift(souvec)
  use phys_constant, only : long
  use grid_parameter, only : nrg, ntg, npg
  use def_metric
  use def_metric_excurve_grid, only : trk_grid
  use def_bh_parameter, only : ome_bh
  use def_vector_x, only : hvec_xg
  use make_array_3d
  use make_array_4d
  use interface_interpo_linear_type0
  use interface_grgrad_midpoint
  use interface_grgrad_midpoint_type0
  use interface_grgrad_gridpoint_type0
  use interface_dadbscalar_3rd2nd_type0
  use def_vector_phi, only : vec_phig
  implicit none
  real(long), pointer :: souvec(:,:,:,:), ddivBeta(:,:,:,:)
  real(long), pointer :: fnc2(:,:,:), divBeta(:,:,:), gradtrK(:,:,:,:)
  real(long), pointer :: grad2x(:,:,:), grad2y(:,:,:), grad2z(:,:,:)
  real(long) :: san, san2, san4
  real(long) :: alpgc, psigc, fnc2gc, afnc2inv, dxafn, dyafn, dzafn
  real(long) :: ovxu, ovyu, ovzu
  real(long) :: tfkax, tfkay, tfkaz, dtrK
  real(long) :: dbxdx,dbxdy,dbxdz, dbydx,dbydy,dbydz, dbzdx,dbzdy,dbzdz 
  real(long) :: d2bvx(3,3), d2bvy(3,3), d2bvz(3,3)
  integer :: ii, irg, itg, ipg, impt
!
  call alloc_array3d(fnc2,0,nrg,0,ntg,0,npg)
  call alloc_array3d(grad2x,1,nrg,1,ntg,1,npg)
  call alloc_array3d(grad2y,1,nrg,1,ntg,1,npg)
  call alloc_array3d(grad2z,1,nrg,1,ntg,1,npg)
  call alloc_array3d(divBeta,0,nrg,0,ntg,0,npg)
  call alloc_array4d(ddivBeta,1,nrg,1,ntg,1,npg,1,3)
  call alloc_array4d(gradtrK,1,nrg,1,ntg,1,npg,1,3)
!
  san = 1.0d0/3.0d0
  san2= 2.0d0/3.0d0
  san4= 4.0d0/3.0d0
!
! --- Source terms of Momentum constraint 
! --- for computing shift.  
!
  call grgrad_midpoint(trk_grid,grad2x,grad2y,grad2z)
  gradtrK(1:nrg,1:ntg,1:npg,1) = grad2x(1:nrg,1:ntg,1:npg)
  gradtrK(1:nrg,1:ntg,1:npg,2) = grad2y(1:nrg,1:ntg,1:npg)
  gradtrK(1:nrg,1:ntg,1:npg,3) = grad2z(1:nrg,1:ntg,1:npg)
!
  do ipg = 0, npg
    do itg = 0, ntg
      do irg = 0, nrg
! value inside of the companion BH causes a problem.
! value inside of the companion BH is set to zero in interpo_gr2gr_4th(_mpt)
        psigc = psi(irg,itg,ipg)
        alpgc = alph(irg,itg,ipg)
        if (alpgc.le.0.0d0) alpgc = 1.0d-20
        fnc2(irg,itg,ipg) = psi(irg,itg,ipg)**6/alpgc
        call grgrad_gridpoint_type0(bvxd,dbxdx,dbxdy,dbxdz,irg,itg,ipg)
        call grgrad_gridpoint_type0(bvyd,dbydx,dbydy,dbydz,irg,itg,ipg)
        call grgrad_gridpoint_type0(bvzd,dbzdx,dbzdy,dbzdz,irg,itg,ipg)
        divBeta(irg,itg,ipg) = dbxdx + dbydy + dbzdz
!!!!        call grgrad_gridpoint_type0(psi,dxafn,dyafn,dzafn,irg,itg,ipg)
!!!!        ovxu = bvxd(irg,itg,ipg) + ome_bh*vec_phig(irg,itg,ipg,1)
!!!!        ovyu = bvyd(irg,itg,ipg) + ome_bh*vec_phig(irg,itg,ipg,2)
!!!!        ovzu = bvzd(irg,itg,ipg) + ome_bh*vec_phig(irg,itg,ipg,3)
!!!!        divBeta(irg,itg,ipg) = - 6.0d0/psigc &
!!!!        &                      * (ovxu*dxafn + ovyu*dyafn + ovzu*dzafn)
      end do
    end do
  end do
!
  call grgrad_midpoint(fnc2,grad2x,grad2y,grad2z)
!
  do ipg = 1, npg
    do itg = 1, ntg
      do irg = 1, nrg
        call grgrad_midpoint_type0(divBeta,dxafn,dyafn,dzafn,irg,itg,ipg)
        ddivBeta(irg,itg,ipg,1) = dxafn        
        ddivBeta(irg,itg,ipg,2) = dyafn
        ddivBeta(irg,itg,ipg,3) = dzafn
!!        call dadbscalar_3rd2nd_type0(bvxd,d2bvx,irg,itg,ipg,'bh')
!!        call dadbscalar_3rd2nd_type0(bvyd,d2bvy,irg,itg,ipg,'bh')
!!        call dadbscalar_3rd2nd_type0(bvzd,d2bvz,irg,itg,ipg,'bh')
!!        ddivBeta(irg,itg,ipg,1) = d2bvx(1,1) + d2bvy(1,2) + d2bvz(1,3)
!!        ddivBeta(irg,itg,ipg,2) = d2bvx(2,1) + d2bvy(2,2) + d2bvz(2,3)
!!        ddivBeta(irg,itg,ipg,3) = d2bvx(3,1) + d2bvy(3,2) + d2bvz(3,3)
      end do
    end do
  end do
!
  do ii = 1, 3
    do ipg = 1, npg
      do itg = 1, ntg
        do irg = 1, nrg
          call interpo_linear_type0(alpgc,alph,irg,itg,ipg)
          call interpo_linear_type0(fnc2gc,fnc2,irg,itg,ipg)
          afnc2inv = alpgc/fnc2gc
!          call interpo_linear_type0(psigc,psi,irg,itg,ipg)
!          if (psigc.le.0.0d0) psigc = 3.0
!          afnc2inv = alpgc**2/psigc**6
          dxafn = grad2x(irg,itg,ipg)
          dyafn = grad2y(irg,itg,ipg)
          dzafn = grad2z(irg,itg,ipg)
          tfkax  = tfkij(irg,itg,ipg,ii,1)
          tfkay  = tfkij(irg,itg,ipg,ii,2)
          tfkaz  = tfkij(irg,itg,ipg,ii,3)
          dtrK = gradtrK(irg,itg,ipg,ii)
!
          souvec(irg,itg,ipg,ii) = - 2.0d0*afnc2inv &
                 &   *(tfkax*dxafn + tfkay*dyafn + tfkaz*dzafn) &
                 &   + san4*alpgc*dtrK &
                 &   - san*ddivBeta(irg,itg,ipg,ii)
!
!if (irg == 2.and.itg == ntg/2-2.and.ipg == 2) then
!write(6,'(1p,3e16.8)') hvec_xg(irg,itg,ipg,1), hvec_xg(irg,itg,ipg,2), &
!&                      hvec_xg(irg,itg,ipg,3)
!write(6,'(1p,3e16.8)') san4*alpgc*dtrK, ddivBeta(irg,itg,ipg,ii)
!end if
!
        end do
      end do
    end do
  end do     
!
  deallocate(fnc2)
  deallocate(grad2x)
  deallocate(grad2y)
  deallocate(grad2z)
  deallocate(divBeta)
  deallocate(ddivBeta)
  deallocate(gradtrK)
!
end subroutine sourceterm_MoC_CF_with_divshift

---