00001 subroutine sourceterm_2pot_test(sou)
00002   use phys_constant, only : long, pi
00003   use grid_parameter, only : nrg, ntg, npg
00004   use def_metric, only : psi, alph
00005 
00006 
00007   use interface_interpo_linear_type0
00008   use interface_grgrad_midpoint_r3rd_type0
00009   use interface_grgrad_midpoint_r4th_type0
00010   use make_array_3d
00011   implicit none
00012   real(long), pointer :: sou(:,:,:)
00013 
00014   real(long) :: psigc, alpgc, alpsigc, psiinv
00015   real(long) :: dxpsi,dypsi,dzpsi,dxalph,dyalph,dzalph
00016   integer    :: irg, itg, ipg
00017 
00018   do ipg = 1, npg
00019     do itg = 1, ntg
00020       do irg = 1, nrg
00021         call interpo_linear_type0(psigc,psi,irg,itg,ipg)
00022         call grgrad_midpoint_r3rd_type0(psi,dxpsi,dypsi,dzpsi,irg,itg,ipg,'bh')
00023         call grgrad_midpoint_r3rd_type0(alph,dxalph,dyalph,dzalph,&
00024                                                           &   irg,itg,ipg,'bh')
00025 
00026 
00027 
00028         psiinv = 1.0d0/psigc
00029         sou(irg,itg,ipg) = -2.0d0*psiinv & 
00030         &                *(dxpsi*dxalph + dypsi*dyalph + dzpsi*dzalph)
00031       end do
00032     end do
00033   end do
00034 
00035 end subroutine sourceterm_2pot_test