00001 subroutine source_quad_pole_peos(souf5d, iquad)
00002 use phys_constant, only : long
00003 use grid_parameter, only : nrf, ntf, npf
00004 use def_matter, only : emd, rs
00005 use def_matter_parameter, only : ber
00006 use def_metric, only : psi, alph
00007 use coordinate_grav_r, only : rg
00008 use trigonometry_grav_phi, only : sinphig, cosphig
00009 use trigonometry_grav_theta, only : sinthg, costhg
00010 use def_vector_x, only : vec_xf
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 :: souf5d(:,:,:,:,:)
00016 real(long), pointer :: alphf(:,:,:), psif(:,:,:)
00017 real(long) :: emdw, alphw, psiw, rhow, hhw, utw
00018 real(long) :: small = 1.0d-15
00019 real(long) :: prew, ene
00020 real(long) :: xx(1:3), rr, ox(1:3), px(1:3)
00021 integer :: ir, it, ip, i, j, iquad
00022
00023 call alloc_array3d(psif, 0, nrf, 0, ntf, 0, npf)
00024 call alloc_array3d(alphf, 0, nrf, 0, ntf, 0, npf)
00025 call interpo_gr2fl(alph, alphf)
00026 call interpo_gr2fl(psi, psif)
00027
00028 do ip = 0, npf
00029 do it = 0, ntf
00030 do ir = 0, nrf
00031 emdw = emd(ir,it,ip)
00032 if (emdw <= small) emdw = small
00033 psiw = psif(ir,it,ip)
00034 alphw = alphf(ir,it,ip)
00035 call peos_q2hprho(emdw, hhw, prew, rhow, ene)
00036 utw = hhw/ber
00037
00038 xx(1) = vec_xf(ir,it,ip,1)
00039 xx(2) = vec_xf(ir,it,ip,2)
00040 xx(3) = vec_xf(ir,it,ip,3)
00041 rr = xx(1)**2+xx(2)**2+xx(3)**2
00042 ox(1:2) = xx(1:2)
00043 ox(3) = 0.0d0
00044 px(1) = vec_phif(ir,it,ip,1)
00045 px(2) = vec_phif(ir,it,ip,2)
00046 px(3) = vec_phif(ir,it,ip,3)
00047
00048 do i=1, 3
00049 do j=1, 3
00050 if (iquad == 0 .or. iquad == 1) &
00051 & souf5d(ir,it,ip,i,j) = xx(i)*xx(j)
00052 if (iquad == 1 .and. i == j) &
00053 & souf5d(ir,it,ip,i,j) = xx(i)*xx(j) - rr/3.0d0
00054 if (iquad == 2) &
00055 & souf5d(ir,it,ip,i,j) = px(i)*xx(j) + xx(i)*px(j)
00056 if (iquad == 3) &
00057 & souf5d(ir,it,ip,i,j) = &
00058 & ox(i)*xx(j) + xx(i)*ox(j) - 2.0d0*px(i)*px(j)
00059 if (iquad == 4) &
00060 & souf5d(ir,it,ip,i,j) = &
00061 & px(i)*xx(j) + 3.0d0*ox(i)*px(j) + 3.0d0*px(i)*ox(j) + xx(i)*px(j)
00062 souf5d(ir,it,ip,i,j) = souf5d(ir,it,ip,i,j)*rhow*alphw*utw*psiw**6
00063 end do
00064 end do
00065
00066 end do
00067 end do
00068 end do
00069
00070 deallocate(alphf)
00071 deallocate(psif)
00072 end subroutine source_quad_pole_peos
