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, utf
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
00047 utw = utf(irf,itf,ipf)
00048
00049 rhoHw = hhw*rhow*(alphw*utw)**2 - prew
00050 esseS = -hhw*rhow + 4.0d0*prew + rhoHw
00051
00052 souf(irf,itf,ipf) = 4.0d0*pi*alphw*psiw**6*(esseS+rhoHw)
00053 end do
00054 end do
00055 end do
00056
00057 deallocate(alphf)
00058 deallocate(psif)
00059 end subroutine source_komar_mass