00001 subroutine source_virial_gravity_CF(sou_Wgra)
00002 use phys_constant, only : long, pi
00003 use grid_parameter, only : nrg, ntg, npg
00004 use def_metric, only : psi, alph, tfkijkij, trk, &
00005 & bvxu, bvyu, bvzu, bvxd, bvyd, bvzd
00006
00007 use interface_grgrad_midpoint
00008 use interface_interpo_linear_type0
00009 use make_array_3d
00010 use make_array_4d
00011 implicit none
00012 integer :: ia, ib, irg, itg, ipg
00013 real(long) :: psifc, psifc6, alpfc, pi4inv
00014 real(long) :: alpgc, psigc, bvxgc, bvygc, bvzgc,
00015 psigc2, psigc6, alp2inv, psal2inv, aijaij, trkgc, &
00016 dapsi(3), daalph(3), dpsi2, dalph2, bvdal
00017 real(long), pointer :: sou_Wgra(:,:,:)
00018 real(long), pointer :: dfdx(:,:,:), dfdy(:,:,:), dfdz(:,:,:)
00019 real(long), pointer :: grada(:,:,:,:), gradp(:,:,:,:)
00020
00021
00022
00023
00024 call alloc_array3d(dfdx,1,nrg,1,ntg,1,npg)
00025 call alloc_array3d(dfdy,1,nrg,1,ntg,1,npg)
00026 call alloc_array3d(dfdz,1,nrg,1,ntg,1,npg)
00027 call alloc_array4d(grada,1,nrg,1,ntg,1,npg,1,3)
00028 call alloc_array4d(gradp,1,nrg,1,ntg,1,npg,1,3)
00029
00030 call grgrad_midpoint(psi,dfdx,dfdy,dfdz)
00031 gradp(1:nrg,1:ntg,1:npg,1) = dfdx(1:nrg,1:ntg,1:npg)
00032 gradp(1:nrg,1:ntg,1:npg,2) = dfdy(1:nrg,1:ntg,1:npg)
00033 gradp(1:nrg,1:ntg,1:npg,3) = dfdz(1:nrg,1:ntg,1:npg)
00034 call grgrad_midpoint(alph,dfdx,dfdy,dfdz)
00035 grada(1:nrg,1:ntg,1:npg,1) = dfdx(1:nrg,1:ntg,1:npg)
00036 grada(1:nrg,1:ntg,1:npg,2) = dfdy(1:nrg,1:ntg,1:npg)
00037 grada(1:nrg,1:ntg,1:npg,3) = dfdz(1:nrg,1:ntg,1:npg)
00038 pi4inv = 1.0d0/(4.0d0*pi)
00039
00040 do ipg = 1, npg
00041 do itg = 1, ntg
00042 do irg = 1, nrg
00043
00044 call interpo_linear_type0(psigc,psi,irg,itg,ipg)
00045 call interpo_linear_type0(alpgc,alph,irg,itg,ipg)
00046 call interpo_linear_type0(bvxgc,bvxu,irg,itg,ipg)
00047 call interpo_linear_type0(bvygc,bvyu,irg,itg,ipg)
00048 call interpo_linear_type0(bvzgc,bvzu,irg,itg,ipg)
00049
00050 psigc2 = psigc**2
00051 psigc6 = psigc**6
00052 alp2inv = 1.0d0/alpgc**2
00053 psal2inv= (psigc/alpgc)**2
00054 aijaij = tfkijkij(irg,itg,ipg)
00055 trkgc = trk(irg,itg,ipg)
00056 dapsi(1:3) = gradp(irg,itg,ipg,1:3)
00057 daalph(1:3) = grada(irg,itg,ipg,1:3)
00058
00059 dpsi2 = 0.0d0
00060 dalph2 = 0.0d0
00061 dpsi2 = dapsi(1)**2 + dapsi(2)**2 + dapsi(3)**2
00062 dalph2 = daalph(1)**2 + daalph(2)**2 + daalph(3)**2
00063 bvdal = bvxgc*daalph(1) + bvygc*daalph(2) + bvzgc*daalph(3)
00064
00065 sou_Wgra(irg,itg,ipg) = pi4inv*(2.0d0*dpsi2 - psal2inv*dalph2 &
00066 & + 0.75d0*psigc6*(aijaij - 2.0d0*trkgc/3.0d0) &
00067 & + psigc6*alp2inv*trkgc*bvdal)
00068
00069 end do
00070 end do
00071 end do
00072
00073 deallocate(dfdx)
00074 deallocate(dfdy)
00075 deallocate(dfdz)
00076 deallocate(grada)
00077 deallocate(gradp)
00078
00079 end subroutine source_virial_gravity_CF
