00001 subroutine test_analytic_bns_solution_mpt(impt)
00002 use phys_constant, only : long, pi
00003 use grid_parameter, only : nrg, ntg, npg, nrf
00004 use def_metric, only : bvxd
00005 use coordinate_grav_r, only : rg
00006 use grid_points_binary_excision, only : rb
00007 implicit none
00008 integer :: irg,itg,ipg, impt
00009 real(long) :: zfac, small = 1.0d-15
00010
00011 do ipg = 0, npg
00012 do itg = 0, ntg
00013 do irg = 0, nrg
00014 if(impt.eq.1) then
00015 if(irg.lt.nrf) &
00016 & bvxd(irg,itg,ipg) = - 2.0d0*pi/3.0d0*(3.0d0-rg(irg)**2) &
00017 & -(8.0d0*pi/3.0d0)*1.0d0/rb(irg,itg,ipg)
00018 if(irg.ge.nrf) &
00019 & bvxd(irg,itg,ipg) = -(4.0d0*pi/3.0d0)*1.0d0/rg(irg) &
00020 & -(8.0d0*pi/3.0d0)*1.0d0/rb(irg,itg,ipg)
00021 end if
00022 if(impt.eq.2) then
00023 if(irg.lt.nrf) &
00024 & bvxd(irg,itg,ipg) = - 4.0d0*pi/3.0d0*(3.0d0-rg(irg)**2) &
00025 & -(4.0d0*pi/3.0d0)*1.0d0/rb(irg,itg,ipg)
00026 if(irg.ge.nrf) &
00027 & bvxd(irg,itg,ipg) = -(8.0d0*pi/3.0d0)*1.0d0/rg(irg) &
00028 & -(4.0d0*pi/3.0d0)*1.0d0/rb(irg,itg,ipg)
00029 end if
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040 end do
00041 end do
00042 end do
00043
00044 end subroutine test_analytic_bns_solution_mpt