00001 subroutine bh_boundary_CF(sou_surf,dsou_surf,char_mp)
00002 use phys_constant, only : long
00003 use grid_parameter, only : nrg, ntg, npg, rgin
00004 use def_metric, only : psi
00005 use def_bh_parameter, only : ome_bh, spin_bh, alph_bh, psi_bh, bh_bctype
00006 use def_binary_parameter, only : sepa
00007 use def_metric_excurve_grid, only : tfkij_grid
00008 use def_vector_bh, only : hvec_bh_cbh_xg,hvec_bh_cm_phig,hvec_bh_cbh_spin
00009 use make_array_2d
00010 use interface_grdr_gridpoint_type0_nosym
00011 use interface_interpo_linear_type0_2Dsurf
00012 implicit none
00013 real(long), pointer :: sou_surf(:,:), dsou_surf(:,:)
00014 real(long), pointer :: psi_bhsurf(:,:), dpsi_bhsurf(:,:)
00015 real(long) :: work(2,2), deriv, dpsidr
00016 real(long) :: Aij_surf, Aij_sisj, si, sj
00017 real(long) :: psigc, psigc3, psigcinv2
00018 real(long) :: vphi_cm, vphi_cbh
00019 integer :: itg, ipg, ii, jj
00020 character(len=4), intent(in) :: char_mp
00021
00022
00023
00024 if (bh_bctype.eq.'TU') then
00025 if (char_mp.eq.'psi ') then
00026 sou_surf(1:ntg,1:npg) = psi_bh
00027 dsou_surf(1:ntg,1:npg)= 0.0d0
00028 end if
00029 if (char_mp.eq.'alps') then
00030 sou_surf(1:ntg,1:npg) = alph_bh*psi_bh
00031 dsou_surf(1:ntg,1:npg)= 0.0d0
00032 end if
00033 if (char_mp.eq.'bvxd'.or.char_mp.eq.'bvyd'.or.char_mp.eq.'bvzd') then
00034 if (char_mp.eq.'bvxd') ii = 1
00035 if (char_mp.eq.'bvyd') ii = 2
00036 if (char_mp.eq.'bvzd') ii = 3
00037 do ipg = 1, npg
00038 do itg = 1, ntg
00039 vphi_cm = hvec_bh_cm_phig(itg,ipg,ii)
00040 sou_surf(itg,ipg) = - ome_bh*vphi_cm
00041 dsou_surf(itg,ipg) = 0.0d0
00042 end do
00043 end do
00044 end if
00045 end if
00046
00047
00048
00049 if (bh_bctype.eq.'AH') then
00050
00051 call alloc_array2d(psi_bhsurf, 0, ntg, 0, npg)
00052 call alloc_array2d(dpsi_bhsurf, 0, ntg, 0, npg)
00053
00054 psi_bhsurf(0:ntg,0:npg) = psi(0,0:ntg,0:npg)
00055
00056 if (char_mp.eq.'psi '.or.char_mp.eq.'alps') then
00057 do ipg = 0, npg
00058 do itg = 0, ntg
00059 call grdr_gridpoint_type0_nosym(psi,deriv,0,itg,ipg)
00060 dpsi_bhsurf(itg,ipg) = deriv
00061 end do
00062 end do
00063
00064 do ipg = 1, npg
00065 do itg = 1, ntg
00066 call interpo_linear_type0_2Dsurf(dpsidr,dpsi_bhsurf,itg,ipg)
00067 call interpo_linear_type0_2Dsurf(psigc,psi_bhsurf,itg,ipg)
00068 psigc3 = psigc**3
00069 Aij_sisj = 0.0d0
00070 do ii = 1, 3
00071 do jj = 1, 3
00072 work(1:2,1:2) = tfkij_grid(0,itg-1:itg,ipg-1:ipg,ii,jj)
00073 call interpo_linear1p_type0_2Dsurf(Aij_surf,work)
00074 si = hvec_bh_cbh_xg(itg,ipg,ii)/rgin
00075 sj = hvec_bh_cbh_xg(itg,ipg,jj)/rgin
00076 Aij_sisj = Aij_sisj + Aij_surf*si*sj
00077 end do
00078 end do
00079 sou_surf(itg,ipg) = 2.0d0*rgin*(-dpsidr - 0.25d0*psigc3*Aij_sisj)
00080 dsou_surf(itg,ipg)= - 0.5d0*psigc/rgin - 0.25d0*psigc3*Aij_sisj
00081 end do
00082 end do
00083 if (char_mp.eq.'alps') then
00084 sou_surf(1:ntg,1:npg) = alph_bh*sou_surf(1:ntg,1:npg)
00085 dsou_surf(1:ntg,1:npg)= alph_bh*dsou_surf(1:ntg,1:npg)
00086 end if
00087 end if
00088
00089 if (char_mp.eq.'bvxd'.or.char_mp.eq.'bvyd'.or.char_mp.eq.'bvzd') then
00090 if (char_mp.eq.'bvxd') ii = 1
00091 if (char_mp.eq.'bvyd') ii = 2
00092 if (char_mp.eq.'bvzd') ii = 3
00093
00094 do ipg = 1, npg
00095 do itg = 1, ntg
00096 call interpo_linear_type0_2Dsurf(psigc,psi_bhsurf,itg,ipg)
00097 psigcinv2 = 1.0d0/psigc**2
00098 si = hvec_bh_cbh_xg(itg,ipg,ii)/rgin
00099 vphi_cm = hvec_bh_cm_phig(itg,ipg,ii)
00100 vphi_cbh = hvec_bh_cbh_spin(itg,ipg,ii)
00101 sou_surf(itg,ipg) = alph_bh*psigcinv2*si &
00102 & - ome_bh*vphi_cm - spin_bh*vphi_cbh
00103 end do
00104 end do
00105 end if
00106
00107 deallocate(psi_bhsurf)
00108 deallocate(dpsi_bhsurf)
00109
00110 end if
00111
00112 end subroutine bh_boundary_CF
