00001 subroutine sourceterm_MoC_peos(souvec,sou)
00002 use phys_constant, only : long, pi
00003 use grid_parameter, only : nrg, ntg, npg
00004 use coordinate_grav_r, only : hrg, rg
00005 use trigonometry_grav_theta, only : hsinthg, hcosthg
00006 use trigonometry_grav_phi, only : hsinphig, hcosphig
00007 use def_metric, only : psi, alph, tfkijkij, bvxd, bvyd, bvzd, tfkij
00008 use def_matter, only : emdg
00009 use def_matter_parameter, only : ome, ber, radi
00010 use def_vector_x, only : hvec_xg
00011 use def_vector_phi, only : hvec_phig
00012 use make_array_3d
00013 use interface_grgrad_midpoint
00014 use interface_interpo_linear_type0
00015 implicit none
00016 real(long), pointer :: sou(:,:,:)
00017 real(long), pointer :: souvec(:,:,:,:)
00018 real(long), pointer :: fnc2(:,:,:)
00019 real(long), pointer :: grad2x(:,:,:), grad2y(:,:,:), grad2z(:,:,:)
00020 real(long) :: vphig(3), xxx, yyy, zzz, san, san2
00021 real(long) :: emdgc, rhogc, pregc, hhgc, utgc, oterm, zfac, rjj
00022 real(long) :: psigc, alpgc, fnc2gc, afnc2inv, dxafn, dyafn, dzafn
00023 real(long) :: bvxdgc, bvydgc, bvzdgc, tfkax, tfkay, tfkaz, ene
00024 integer :: ii, irg, itg, ipg
00025
00026 call alloc_array3d(fnc2,0,nrg,0,ntg,0,npg)
00027 call alloc_array3d(grad2x,1,nrg,1,ntg,1,npg)
00028 call alloc_array3d(grad2y,1,nrg,1,ntg,1,npg)
00029 call alloc_array3d(grad2z,1,nrg,1,ntg,1,npg)
00030
00031 san = 1.0d0/3.0d0
00032 san2= 2.0d0/3.0d0
00033
00034
00035
00036
00037 do ipg = 0, npg
00038 do itg = 0, ntg
00039 do irg = 0, nrg
00040 fnc2(irg,itg,ipg) = psi(irg,itg,ipg)**6/alph(irg,itg,ipg)
00041 end do
00042 end do
00043 end do
00044
00045 call grgrad_midpoint(fnc2,grad2x,grad2y,grad2z)
00046
00047 do ii = 1, 3
00048 do ipg = 1, npg
00049 do itg = 1, ntg
00050 do irg = 1, nrg
00051 call interpo_linear_type0(emdgc,emdg,irg,itg,ipg)
00052 call interpo_linear_type0(psigc,psi,irg,itg,ipg)
00053 call interpo_linear_type0(alpgc,alph,irg,itg,ipg)
00054 call interpo_linear_type0(fnc2gc,fnc2,irg,itg,ipg)
00055 call interpo_linear_type0(bvxdgc,bvxd,irg,itg,ipg)
00056 call interpo_linear_type0(bvydgc,bvyd,irg,itg,ipg)
00057 call interpo_linear_type0(bvzdgc,bvzd,irg,itg,ipg)
00058 afnc2inv = alpgc/fnc2gc
00059 dxafn = grad2x(irg,itg,ipg)
00060 dyafn = grad2y(irg,itg,ipg)
00061 dzafn = grad2z(irg,itg,ipg)
00062 tfkax = tfkij(irg,itg,ipg,ii,1)
00063 tfkay = tfkij(irg,itg,ipg,ii,2)
00064 tfkaz = tfkij(irg,itg,ipg,ii,3)
00065
00066 zfac = 1.0d0
00067 if (emdgc <= 1.0d-15) then
00068 emdgc = 1.0d-15
00069 zfac = 0.0d0
00070 end if
00071 call peos_q2hprho(emdgc, hhgc, pregc, rhogc, ene)
00072 utgc = hhgc/ber
00073 vphig(1) = hvec_phig(irg,itg,ipg,1)
00074 vphig(2) = hvec_phig(irg,itg,ipg,2)
00075 vphig(3) = hvec_phig(irg,itg,ipg,3)
00076 oterm = 0.0d0
00077 if (ii == 1) oterm = bvxdgc + ome*vphig(1)
00078 if (ii == 2) oterm = bvydgc + ome*vphig(2)
00079 if (ii == 3) oterm = bvzdgc + ome*vphig(3)
00080 rjj = hhgc*rhogc*alpgc*utgc**2*psigc**4*oterm
00081
00082 souvec(irg,itg,ipg,ii) = - 2.0d0*afnc2inv &
00083 & *(tfkax*dxafn + tfkay*dyafn + tfkaz*dzafn) &
00084 & + radi**2*16.0d0*pi*alpgc*rjj*zfac
00085 end do
00086 end do
00087 end do
00088 end do
00089
00090
00091
00092 do ipg = 1, npg
00093 do itg = 1, ntg
00094 do irg = 1, nrg
00095 xxx = hvec_xg(irg,itg,ipg,1)
00096 yyy = hvec_xg(irg,itg,ipg,2)
00097 zzz = hvec_xg(irg,itg,ipg,3)
00098 sou(irg,itg,ipg) = xxx*souvec(irg,itg,ipg,1) &
00099 & + yyy*souvec(irg,itg,ipg,2) &
00100 & + zzz*souvec(irg,itg,ipg,3)
00101 end do
00102 end do
00103 end do
00104
00105 deallocate(fnc2)
00106 deallocate(grad2x)
00107 deallocate(grad2y)
00108 deallocate(grad2z)
00109 end subroutine sourceterm_MoC_peos
