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

Go to the documentation of this file.

subroutine source_virial_WL(sou_Tkin,sou_Pint,sou_Wgra)
  use phys_constant,  only : long, pi
  use grid_parameter, only : nrg, ntg, npg, nrf, ntf, npf
  use def_metric, only : psi, alph, tfkijkij, &
  &                      bvxu, bvyu, bvzu, bvxd, bvyd, bvzd
  use def_metric_hij, only : hxxu, hxyu, hxzu, hyyu, hyzu, hzzu
  use def_ricci_tensor, only : rab
  use def_metric_on_SFC_WL, only : psif, alphf
  use def_matter,     only : emd, utf
  use interface_grgrad_midpoint
  use interface_interpo_linear_type0
  use make_array_3d
  use make_array_4d
  implicit none
  integer :: irf, itf, ipf, ia, ib, irg, itg, ipg
  real(long) :: psifc, psifc6, alpfc, pi4inv
  real(long) :: emdfc, hhfc, prefc, rhofc, enefc, utfc
  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
  real(long), pointer :: sou_Tkin(:,:,:), sou_Pint(:,:,:), sou_Wgra(:,:,:)
  real(long), pointer :: dfdx(:,:,:), dfdy(:,:,:), dfdz(:,:,:)
  real(long), pointer :: grada(:,:,:,:), gradp(:,:,:,:)
!
! --- Compute source terms for virial relations.
! --- sources for kinetic and internal energies are defined on SCF grid points.
! --- sources for EM, and gravitational energies are 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)
!
  do ipf = 0, npf
    do itf = 0, ntf
      do irf = 0, nrf
!
        emdfc = emd(irf,itf,ipf)
        call peos_q2hprho(emdfc, hhfc, prefc, rhofc, enefc)
        utfc = utf(irf,itf,ipf)
!
        psifc = psif(irf,itf,ipf)
        alpfc = alphf(irf,itf,ipf)
        psifc6 = psifc**6
!
        sou_Tkin(irf,itf,ipf)= 0.5d0*hhfc*rhofc*((alpfc*utfc)**2-1.0d0)*psifc6
        sou_Pint(irf,itf,ipf)= prefc*psifc6
!
      end do
    end do
  end do
!
  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)
!
  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  = 0.0d0
        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 - 2.0d0*trkgc/3.0d0) &
        &                     + 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_WL

---