00001 subroutine source_komar_mass(soug,souf)
00002   use phys_constant, only  :   long, pi
00003   use grid_parameter, only  :   nrg, ntg, npg, nrf, ntf, npf
00004   use def_matter, only  :   emdg, emd
00005   use def_matter_parameter, only  :   radi, pinx, ber
00006   use def_metric, only  :   tfkijkij,  psi, alph
00007   use make_array_3d
00008   use interface_interpo_gr2fl
00009   use interface_interpo_linear_type0
00010   implicit none
00011   real(long),pointer ::   soug(:,:,:), souf(:,:,:)
00012   real(long), pointer ::  alphf(:,:,:), psif(:,:,:) 
00013   real(long)  ::   psiwm6
00014   real(long)  ::   emdw, alphw, psiw, rhow, prew, hhw, utw, rhoHw, esseS
00015   real(long)  ::   zfac, small = 1.0d-15
00016   real(long)  ::   otermx, otermy, otermz
00017   integer     ::   irg,itg,ipg,irf,itf,ipf
00018 
00019   call alloc_array3d(psif, 0, nrf, 0, ntf, 0, npf)
00020   call alloc_array3d(alphf, 0, nrf, 0, ntf, 0, npf)
00021   call interpo_gr2fl(alph, alphf)
00022   call interpo_gr2fl(psi, psif)
00023 
00024   do ipg = 1, npg
00025     do itg = 1, ntg
00026       do irg = 1, nrg
00027         call interpo_linear_type0(psiw,psi,irg,itg,ipg)
00028         call interpo_linear_type0(alphw,alph,irg,itg,ipg)
00029         psiwm6 = psiw**6
00030         soug(irg,itg,ipg) = alphw*psiwm6*tfkijkij(irg,itg,ipg)
00031       end do
00032     end do
00033   end do
00034 
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         rhow = emdw**pinx
00044         prew = rhow*emdw
00045         hhw  = 1.0d0 + (pinx+1.0d0)*emdw
00046         utw = hhw/ber
00047 
00048         rhoHw = hhw*rhow*(alphw*utw)**2 - prew
00049         esseS = -hhw*rhow + 4.0d0*prew + rhoHw
00050 
00051         souf(irf,itf,ipf) = 4.0d0*pi*alphw*psiw**6*(esseS+rhoHw)
00052       end do
00053     end do
00054   end do
00055 
00056   deallocate(alphf)
00057   deallocate(psif)
00058 end subroutine source_komar_mass