/home/uryu/Cocal_svn_checkinout/Cocal/code/Subroutine/source_virial_gravity_WL.f90

Go to the documentation of this file.

subroutine source_virial_gravity_WL(sou_Wgra)
  use phys_constant,  only : long, pi
  use grid_parameter, only : nrg, ntg, npg
  use def_metric, only : psi, alph, tfkijkij, trk, &
  &                      bvxu, bvyu, bvzu, bvxd, bvyd, bvzd
  use def_metric_hij, only : hxxu, hxyu, hxzu, hyyu, hyzu, hzzu
  use def_ricci_tensor, only : rab
  use interface_grgrad_midpoint
  use interface_interpo_linear_type0
  use make_array_3d
  use make_array_4d
  implicit none
  integer :: ia, ib, irg, itg, ipg
  real(long) :: pi4inv
  real(long) :: alpgc, psigc, bvxgc, bvygc, bvzgc, gamu(3,3), 
               hhxxu, hhxyu, hhxzu, hhyxu, hhyyu, hhyzu, &
               hhzxu, hhzyu, hhzzu, ric(3,3), ricsc, &
               psigc2, psigc6, alp2inv, psal2inv, aijaij, trkgc, &
               dapsi(3), daalph(3), dpsi2, dalph2, bvdal, fac23
  real(long), pointer :: sou_Wgra(:,:,:)
  real(long), pointer :: dfdx(:,:,:), dfdy(:,:,:), dfdz(:,:,:)
  real(long), pointer :: grada(:,:,:,:), gradp(:,:,:,:)
!
! --- Compute source terms for virial relations.
! --- source for gravitational energies is defined on mid points.
!
  call alloc_array3d(dfdx,1,nrg,1,ntg,1,npg)
  call alloc_array3d(dfdy,1,nrg,1,ntg,1,npg)
  call alloc_array3d(dfdz,1,nrg,1,ntg,1,npg)
  call alloc_array4d(grada,1,nrg,1,ntg,1,npg,1,3)
  call alloc_array4d(gradp,1,nrg,1,ntg,1,npg,1,3)
!
  call grgrad_midpoint(psi,dfdx,dfdy,dfdz)
  gradp(1:nrg,1:ntg,1:npg,1) = dfdx(1:nrg,1:ntg,1:npg)
  gradp(1:nrg,1:ntg,1:npg,2) = dfdy(1:nrg,1:ntg,1:npg)
  gradp(1:nrg,1:ntg,1:npg,3) = dfdz(1:nrg,1:ntg,1:npg)
  call grgrad_midpoint(alph,dfdx,dfdy,dfdz)
  grada(1:nrg,1:ntg,1:npg,1) = dfdx(1:nrg,1:ntg,1:npg)
  grada(1:nrg,1:ntg,1:npg,2) = dfdy(1:nrg,1:ntg,1:npg)
  grada(1:nrg,1:ntg,1:npg,3) = dfdz(1:nrg,1:ntg,1:npg)
  pi4inv = 1.0d0/(4.0d0*pi)
!
  fac23 = 2.0d0/3.0d0
  do ipg = 1, npg
    do itg = 1, ntg
      do irg = 1, nrg
!
        call interpo_linear_type0(psigc,psi,irg,itg,ipg)
        call interpo_linear_type0(alpgc,alph,irg,itg,ipg)
        call interpo_linear_type0(bvxgc,bvxu,irg,itg,ipg)
        call interpo_linear_type0(bvygc,bvyu,irg,itg,ipg)
        call interpo_linear_type0(bvzgc,bvzu,irg,itg,ipg)
        call interpo_linear_type0(hhxxu,hxxu,irg,itg,ipg)
        call interpo_linear_type0(hhxyu,hxyu,irg,itg,ipg)
        call interpo_linear_type0(hhxzu,hxzu,irg,itg,ipg)
        call interpo_linear_type0(hhyyu,hyyu,irg,itg,ipg)
        call interpo_linear_type0(hhyzu,hyzu,irg,itg,ipg)
        call interpo_linear_type0(hhzzu,hzzu,irg,itg,ipg)
        hhyxu = hhxyu
        hhzxu = hhxzu
        hhzyu = hhyzu
        gamu(1,1) = hhxxu + 1.0d0
        gamu(1,2) = hhxyu
        gamu(1,3) = hhxzu
        gamu(2,2) = hhyyu + 1.0d0
        gamu(2,3) = hhyzu
        gamu(3,3) = hhzzu + 1.0d0
        gamu(2,1) = gamu(1,2)
        gamu(3,1) = gamu(1,3)
        gamu(3,2) = gamu(2,3)
        psigc2 = psigc**2
        psigc6 = psigc**6
        alp2inv = 1.0d0/alpgc**2
        psal2inv= (psigc/alpgc)**2
        aijaij = tfkijkij(irg,itg,ipg)
        trkgc  = trk(irg,itg,ipg)
        ric(1,1) = rab(irg,itg,ipg,1)
        ric(1,2) = rab(irg,itg,ipg,2)
        ric(1,3) = rab(irg,itg,ipg,3)
        ric(2,2) = rab(irg,itg,ipg,4)
        ric(2,3) = rab(irg,itg,ipg,5)
        ric(3,3) = rab(irg,itg,ipg,6)
        ric(2,1) = ric(1,2)
        ric(3,1) = ric(1,3)
        ric(3,2) = ric(2,3)        
        dapsi(1:3) = gradp(irg,itg,ipg,1:3)
        daalph(1:3) = grada(irg,itg,ipg,1:3)
!
        dpsi2  = 0.0d0
        dalph2 = 0.0d0
        ricsc  = 0.0d0
        do ib = 1, 3
          do ia = 1, 3
            dpsi2  = dpsi2 + gamu(ia,ib)*dapsi( ia)*dapsi( ib)
            dalph2 = dalph2+ gamu(ia,ib)*daalph(ia)*daalph(ib)
            ricsc  = ricsc + gamu(ia,ib)*ric(ia,ib)
          end do
        end do
        bvdal = bvxgc*daalph(1) +  bvygc*daalph(2) +  bvzgc*daalph(3)
!
        sou_Wgra(irg,itg,ipg) = pi4inv*(2.0d0*dpsi2 - psal2inv*dalph2 &
        &                     + 0.75d0*psigc6*(aijaij - fac23*trkgc**2) &
        &                     + psigc6*alp2inv*trkgc*bvdal + 0.25*ricsc*psigc2)
!
      end do
    end do
  end do
!
  deallocate(dfdx)
  deallocate(dfdy)
  deallocate(dfdz)
  deallocate(grada)
  deallocate(gradp)
!
end subroutine source_virial_gravity_WL

---