00001 subroutine source_ang_mom_peos(souf)
00002 use phys_constant, only : long, pi
00003 use grid_parameter, only : nrg, ntg, npg, nrf, ntf, npf
00004 use def_matter, only : emd, rs, utf, omef
00005 use def_matter_parameter
00006 use def_metric, only : tfkijkij, psi, alph, &
00007 bvxd, bvyd, bvzd
00008 use coordinate_grav_r, only : rg
00009 use trigonometry_grav_theta, only : sinthg, costhg
00010 use trigonometry_grav_phi, only : sinphig, cosphig
00011 use def_vector_phi, only : vec_phif
00012 use make_array_3d
00013 use interface_interpo_gr2fl
00014 implicit none
00015 real(long), pointer :: souf(:,:,:)
00016 integer :: ir,it,ip
00017 real(long) :: emdw, alphw, psiw, rhow, prew, hhw, utw, omew, ene
00018 real(long) :: rjjx, rjjy, rjjz, rjjphi
00019 real(long) :: zfac, small = 1.0d-15
00020 real(long) :: vphif(1:3)
00021 real(long) :: otermx, otermy, otermz, bvxufw, bvyufw, bvzufw
00022 real(long), pointer :: alphf(:,:,:), psif(:,:,:)
00023 real(long), pointer :: bvxdf(:,:,:), bvydf(:,:,:), bvzdf(:,:,:)
00024
00025 call alloc_array3d(psif, 0, nrf, 0, ntf, 0, npf)
00026 call alloc_array3d(alphf, 0, nrf, 0, ntf, 0, npf)
00027 call alloc_array3d(bvxdf, 0, nrf, 0, ntf, 0, npf)
00028 call alloc_array3d(bvydf, 0, nrf, 0, ntf, 0, npf)
00029 call alloc_array3d(bvzdf, 0, nrf, 0, ntf, 0, npf)
00030
00031 call interpo_gr2fl(alph, alphf)
00032 call interpo_gr2fl(psi , psif )
00033 call interpo_gr2fl(bvxd, bvxdf)
00034 call interpo_gr2fl(bvyd, bvydf)
00035 call interpo_gr2fl(bvzd, bvzdf)
00036
00037 do ip = 0, npf
00038 do it = 0, ntf
00039 do ir = 0, nrf
00040 emdw = emd(ir,it,ip)
00041 if (emdw <= small) emdw = small
00042 psiw = psif(ir,it,ip)
00043 alphw = alphf(ir,it,ip)
00044 bvxufw = bvxdf(ir,it,ip)
00045 bvyufw = bvydf(ir,it,ip)
00046 bvzufw = bvzdf(ir,it,ip)
00047 call peos_q2hprho(emdw, hhw, prew, rhow, ene)
00048
00049 utw = utf(ir,it,ip)
00050 omew = omef(ir,it,ip)
00051
00052 vphif(1) = vec_phif(ir,it,ip,1)
00053 vphif(2) = vec_phif(ir,it,ip,2)
00054 vphif(3) = vec_phif(ir,it,ip,3)
00055
00056 otermx = bvxufw + omew*vphif(1)
00057 otermy = bvyufw + omew*vphif(2)
00058 otermz = bvzufw + omew*vphif(3)
00059
00060 zfac = 1.0d0
00061 if (emdw <= small) zfac = 0.0d0
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 rjjphi = rjjx*vphif(1) + rjjy*vphif(2) + rjjz*vphif(3)
00067 souf(ir,it,ip) = rjjphi*psiw**6*zfac
00068
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_ang_mom_peos
