00001 subroutine source_komar_mass_compact(souf)
00002   use phys_constant, only  : long, pi
00003   use grid_parameter, only : nrf, ntf, npf
00004   use def_matter
00005   use def_matter_parameter, only  :   radi, pinx, ber, ome
00006   use def_metric, only  :   psi, alph, bvxd, bvyd, bvzd
00007   use coordinate_grav_r, only       : rg
00008   use trigonometry_grav_theta, only : sinthg, costhg
00009   use trigonometry_grav_phi, only   : sinphig, cosphig
00010   use def_vector_phi, only : vec_phif
00011   use make_array_3d
00012   use interface_interpo_gr2fl
00013   implicit none
00014   real(long), pointer :: souf(:,:,:)
00015   real(long), pointer :: alphf(:,:,:),  psif(:,:,:) 
00016   real(long), pointer :: bvxdf(:,:,:), bvydf(:,:,:), bvzdf(:,:,:)
00017   real(long)  ::   emdw, alphw, psiw, rhow, prew, hhw, utw, rhoHw, esseS
00018   real(long)  ::   rjjx, rjjy, rjjz, rjjbeta
00019   real(long)  ::   vphif(1:3)
00020   real(long)  ::   otermx, otermy, otermz, bvxufw, bvyufw, bvzufw
00021   integer     ::   irf,itf,ipf
00022   real(long)  ::   zfac, small = 1.0d-15
00023 
00024   call alloc_array3d(psif, 0, nrf, 0, ntf, 0, npf)
00025   call alloc_array3d(alphf, 0, nrf, 0, ntf, 0, npf)
00026   call alloc_array3d(bvxdf, 0, nrf, 0, ntf, 0, npf)
00027   call alloc_array3d(bvydf, 0, nrf, 0, ntf, 0, npf)
00028   call alloc_array3d(bvzdf, 0, nrf, 0, ntf, 0, npf)
00029 
00030   call interpo_gr2fl(alph, alphf)
00031   call interpo_gr2fl(psi, psif)
00032   call interpo_gr2fl(bvxd, bvxdf)
00033   call interpo_gr2fl(bvyd, bvydf)
00034   call interpo_gr2fl(bvzd, bvzdf)
00035 
00036   do ipf = 0, npf
00037     do itf = 0, ntf
00038       do irf = 0, nrf
00039         emdw = emd(irf,itf,ipf)
00040         if (emdw <= small) emdw = small
00041         psiw = psif(irf,itf,ipf)
00042         alphw = alphf(irf,itf,ipf)
00043         bvxufw = bvxdf(irf,itf,ipf)
00044         bvyufw = bvydf(irf,itf,ipf)
00045         bvzufw = bvzdf(irf,itf,ipf)
00046         rhow = emdw**pinx
00047         prew = rhow*emdw
00048         hhw  = 1.0d0 + (pinx+1.0d0)*emdw
00049         utw = hhw/ber
00050 
00051         vphif(1) = vec_phif(irf,itf,ipf,1)
00052         vphif(2) = vec_phif(irf,itf,ipf,2)
00053         vphif(3) = vec_phif(irf,itf,ipf,3)
00054 
00055         otermx = bvxufw + ome*vphif(1)
00056         otermy = bvyufw + ome*vphif(2)
00057         otermz = bvzufw + ome*vphif(3)
00058 
00059         rhoHw = hhw*rhow*(alphw*utw)**2 - prew
00060         esseS = -hhw*rhow + 4.0d0*prew + rhoHw
00061 
00062         rjjx = hhw*rhow*alphw*utw**2*psiw**4*otermx
00063         rjjy = hhw*rhow*alphw*utw**2*psiw**4*otermy
00064         rjjz = hhw*rhow*alphw*utw**2*psiw**4*otermz
00065 
00066         rjjbeta = rjjx*bvxufw + rjjy*bvyufw + rjjz*bvzufw
00067 
00068         souf(irf,itf,ipf) = (alphw*(esseS+rhoHw) - 2.0d0*rjjbeta)*psiw**6
00069       end do
00070     end do
00071   end do
00072 
00073   deallocate(alphf)
00074   deallocate(psif)
00075   deallocate(bvxdf)
00076   deallocate(bvydf)
00077   deallocate(bvzdf)
00078 end subroutine source_komar_mass_compact