00001 subroutine source_ang_mom_smarr(sous)
00002   use phys_constant, only : long
00003   use grid_parameter, only : nrg, ntg, npg
00004   use make_array_2d
00005   use def_metric
00006   use def_metric_excurve_grid, only : tfkij_grid
00007   use def_bh_parameter, only : ome_bh
00008   use coordinate_grav_r, only : rg
00009   use interface_interpo_linear_type0_2Dsurf
00010   use def_vector_bh, only : hvec_bh_cbh_xg, hvec_bh_cm_phig
00011   use interface_sourceterm_surface_int
00012   implicit none
00013   real(long), pointer :: sous(:,:), sous1(:,:), sous2(:,:), psi_bh(:,:)
00014   real(long), pointer :: bvxd_bh(:,:), bvyd_bh(:,:), bvzd_bh(:,:)
00015   real(long), pointer :: sou_bhsurf(:,:), dsou_bhsurf(:,:)
00016   real(long) :: ni, Aij_surf, val, psi2, psi6, work(2,2), bv(3), beta, vphi_cm
00017   integer    :: irg, itg, ipg, ia, ib
00018 
00019   call alloc_array2d(psi_bh,0,ntg,0,npg)
00020   call alloc_array2d(bvxd_bh,0,ntg,0,npg)
00021   call alloc_array2d(bvyd_bh,0,ntg,0,npg)
00022   call alloc_array2d(bvzd_bh,0,ntg,0,npg)
00023   call alloc_array2d(sous1, 0,ntg,0,npg)
00024   call alloc_array2d(sous2, 0,ntg,0,npg)
00025   call alloc_array2d(sou_bhsurf ,0,ntg,0,npg)
00026   call alloc_array2d(dsou_bhsurf,0,ntg,0,npg)
00027 
00028   psi_bh(0:ntg,0:npg)  = psi(0,0:ntg,0:npg)
00029   bvxd_bh(0:ntg,0:npg) = bvxd(0,0:ntg,0:npg)
00030   bvyd_bh(0:ntg,0:npg) = bvyd(0,0:ntg,0:npg)
00031   bvzd_bh(0:ntg,0:npg) = bvzd(0,0:ntg,0:npg)
00032 
00033   call sourceterm_surface_int(alph,0,sou_bhsurf,dsou_bhsurf)
00034 
00035   do ipg = 1, npg
00036     do itg = 1, ntg
00037       call interpo_linear_type0_2Dsurf(val,psi_bh,itg,ipg)
00038       psi2 = val**2
00039       psi6 = val**6
00040       sous1(itg,ipg) = -psi2*dsou_bhsurf(itg,ipg)
00041 
00042       call interpo_linear_type0_2Dsurf(bv(1),bvxd_bh,itg,ipg)
00043       call interpo_linear_type0_2Dsurf(bv(2),bvyd_bh,itg,ipg)
00044       call interpo_linear_type0_2Dsurf(bv(3),bvzd_bh,itg,ipg)
00045       sous2(itg,ipg)=0.0d0
00046       do ib = 1, 3
00047         do ia = 1, 3
00048           work(1:2,1:2) = tfkij_grid(0,itg-1:itg,ipg-1:ipg,ia,ib)
00049           call interpo_linear1p_type0_2Dsurf(Aij_surf,work)
00050 
00051 
00052           ni = -hvec_bh_cbh_xg(itg,ipg,ib)/rg(0)
00053           vphi_cm = hvec_bh_cm_phig(itg,ipg,ia)
00054           beta = bv(ia) + ome_bh*vphi_cm
00055 
00056           sous2(itg,ipg) = sous2(itg,ipg) + Aij_surf*beta*ni
00057         end do
00058       end do
00059       sous2(itg,ipg) = sous2(itg,ipg)*psi6
00060 
00061       sous(itg,ipg)  = sous1(itg,ipg) - sous2(itg,ipg)
00062     end do
00063   end do
00064 
00065   deallocate(psi_bh)
00066   deallocate(bvxd_bh)
00067   deallocate(bvyd_bh)
00068   deallocate(bvzd_bh)
00069   deallocate(sous1)
00070   deallocate(sous2)
00071   deallocate(sou_bhsurf)
00072   deallocate(dsou_bhsurf)
00073 end subroutine source_ang_mom_smarr
00074