00001 subroutine bh_boundary_psi_test_mpt(impt,char_bc,sou_surf,dsou_surf)
00002 use phys_constant, only : long
00003 use grid_parameter, only : nrg, ntg, npg, rgin
00004 use trigonometry_grav_theta, only : hsinthg
00005 use trigonometry_grav_phi, only : hcosphig
00006 use def_binary_parameter, only : sepa
00007 use grid_parameter_mpt, only : grid_param_real_
00008 implicit none
00009 character(len=1), intent(in) :: char_bc
00010 real(long), pointer :: sou_surf(:,:), dsou_surf(:,:)
00011 real(long) :: st, cp, rad1, rad2, dr2dr1, bhmass
00012 integer :: itg, ipg
00013 integer :: impt, impt_ex
00014 real(long) :: bhmass_ex, rgin_ex
00015
00016 if (impt.eq.1) impt_ex = 2
00017 if (impt.eq.2) impt_ex = 1
00018 rgin_ex = grid_param_real_(1,impt_ex)
00019 bhmass_ex = 2.0d0*rgin_ex
00020 bhmass = 2.0d0*rgin
00021
00022 do ipg = 1, npg
00023 do itg = 1, ntg
00024 st = hsinthg(itg)
00025 cp = hcosphig(ipg)
00026 rad1 = rgin
00027 rad2 = sqrt(rad1**2 - 2.0d0*rad1*sepa*st*cp + sepa**2)
00028 dr2dr1 = (rad1 - sepa*st*cp)/rad2
00029 if (char_bc.eq.'d') then
00030 sou_surf(itg,ipg) = 1.0d0 + 0.5d0*bhmass/rad1 + 0.5d0*bhmass_ex/rad2
00031 end if
00032 if (char_bc.eq.'n') then
00033 dsou_surf(itg,ipg) = - 0.5d0*bhmass/rad1**2 &
00034 & - 0.5d0*bhmass_ex/rad2**2*dr2dr1
00035 end if
00036 end do
00037 end do
00038
00039 end subroutine bh_boundary_psi_test_mpt