00001 subroutine source_komar_mass_compact_peos(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, 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, omew
00018 real(long) :: rjjx, rjjy, rjjz, rjjbeta, ene, rhoHw, esseS
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 call peos_q2hprho(emdw, hhw, prew, rhow, ene)
00047 utw = utf(irf,itf,ipf)
00048 omew = omef(irf,itf,ipf)
00049
00050 vphif(1) = vec_phif(irf,itf,ipf,1)
00051 vphif(2) = vec_phif(irf,itf,ipf,2)
00052 vphif(3) = vec_phif(irf,itf,ipf,3)
00053
00054 otermx = bvxufw + omew*vphif(1)
00055 otermy = bvyufw + omew*vphif(2)
00056 otermz = bvzufw + omew*vphif(3)
00057
00058 rhoHw = hhw*rhow*(alphw*utw)**2 - prew
00059 esseS = -hhw*rhow + 4.0d0*prew + rhoHw
00060
00061 rjjx = hhw*rhow*alphw*utw**2*psiw**4*otermx
00062 rjjy = hhw*rhow*alphw*utw**2*psiw**4*otermy
00063 rjjz = hhw*rhow*alphw*utw**2*psiw**4*otermz
00064
00065 rjjbeta = rjjx*bvxufw + rjjy*bvyufw + rjjz*bvzufw
00066
00067 souf(irf,itf,ipf) = (alphw*(esseS+rhoHw) - 2.0d0*rjjbeta)*psiw**6
00068 end do
00069 end do
00070 end do
00071
00072 deallocate(alphf)
00073 deallocate(psif)
00074 deallocate(bvxdf)
00075 deallocate(bvydf)
00076 deallocate(bvzdf)
00077 end subroutine source_komar_mass_compact_peos
