00001 subroutine source_ang_mom(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, rs, utf
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
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 rhow = emdw**pinx
00048 prew = rhow*emdw
00049 hhw = 1.0d0 + (pinx+1.0d0)*emdw
00050
00051
00052 utw = utf(ir,it,ip)
00053
00054 vphif(1) = vec_phif(ir,it,ip,1)
00055 vphif(2) = vec_phif(ir,it,ip,2)
00056 vphif(3) = vec_phif(ir,it,ip,3)
00057
00058 otermx = bvxufw + ome*vphif(1)
00059 otermy = bvyufw + ome*vphif(2)
00060 otermz = bvzufw + ome*vphif(3)
00061
00062 zfac = 1.0d0
00063 if (emdw <= small) zfac = 0.0d0
00064 rjjx = hhw*rhow*alphw*utw**2*psiw**4*otermx
00065 rjjy = hhw*rhow*alphw*utw**2*psiw**4*otermy
00066 rjjz = hhw*rhow*alphw*utw**2*psiw**4*otermz
00067
00068 rjjphi = rjjx*vphif(1) + rjjy*vphif(2) + rjjz*vphif(3)
00069 souf(ir,it,ip) = rjjphi*psiw**6*zfac
00070
00071 end do
00072 end do
00073 end do
00074
00075 deallocate(alphf)
00076 deallocate(psif)
00077 deallocate(bvxdf)
00078 deallocate(bvydf)
00079 deallocate(bvzdf)
00080 end subroutine source_ang_mom
