00001 subroutine sourceterm_trG_WL(sou)
00002 use def_metric, only : psi, alps, tfkijkij
00003 use grid_parameter, only : nrg, ntg, npg
00004 use def_gamma_crist, only : gmcrix, gmcriy, gmcriz
00005 use def_ricci_tensor, only : rab
00006 use def_metric_hij, only : hxxu, hxyu, hxzu, hyyu, hyzu, hzzu
00007
00008 use interface_interpo_linear_type0
00009 use interface_grgrad_midpoint
00010 use interface_dadbscalar_type0
00011 use interface_dadbscalar_type3
00012 use make_array_4d
00013 use make_array_3d
00014 implicit none
00015
00016 real(8), pointer :: sou(:,:,:)
00017 real(8), pointer :: grad(:,:,:,:), dfdx(:,:,:), dfdy(:,:,:), dfdz(:,:,:)
00018 real(8) :: dabfnc(3,3)
00019 real(8) :: alpsigc, gcdp, hd2p, rics,
00020 rxx, rxy, rxz, ryx, ryy, ryz, rzx, rzy, rzz, &
00021 gmxxu, gmxyu, gmxzu, gmyyu, gmyzu, gmzzu, gmyxu, gmzxu, gmzyu, &
00022 hhxxu, hhxyu, hhxzu, hhyyu, hhyzu, hhzzu, hhyxu, hhzxu, hhzyu
00023 integer :: ipg, irg, itg
00024
00025 call alloc_array4d(grad,1,nrg,1,ntg,1,npg,1,3)
00026 call alloc_array3d(dfdx,1,nrg,1,ntg,1,npg)
00027 call alloc_array3d(dfdy,1,nrg,1,ntg,1,npg)
00028 call alloc_array3d(dfdz,1,nrg,1,ntg,1,npg)
00029
00030
00031
00032 call grgrad_midpoint(alps,dfdx,dfdy,dfdz)
00033 grad(1:nrg,1:ntg,1:npg,1) = dfdx(1:nrg,1:ntg,1:npg)
00034 grad(1:nrg,1:ntg,1:npg,2) = dfdy(1:nrg,1:ntg,1:npg)
00035 grad(1:nrg,1:ntg,1:npg,3) = dfdz(1:nrg,1:ntg,1:npg)
00036
00037 do ipg = 1, npg
00038 do itg = 1, ntg
00039 do irg = 1, nrg
00040
00041 call interpo_linear_type0(alpsigc,alps,irg,itg,ipg)
00042
00043 call interpo_linear_type0(hhxxu,hxxu,irg,itg,ipg)
00044 call interpo_linear_type0(hhxyu,hxyu,irg,itg,ipg)
00045 call interpo_linear_type0(hhxzu,hxzu,irg,itg,ipg)
00046 call interpo_linear_type0(hhyyu,hyyu,irg,itg,ipg)
00047 call interpo_linear_type0(hhyzu,hyzu,irg,itg,ipg)
00048 call interpo_linear_type0(hhzzu,hzzu,irg,itg,ipg)
00049 hhyxu = hhxyu
00050 hhzxu = hhxzu
00051 hhzyu = hhyzu
00052 gmxxu = hhxxu + 1.0d0
00053 gmxyu = hhxyu
00054 gmxzu = hhxzu
00055 gmyyu = hhyyu + 1.0d0
00056 gmyzu = hhyzu
00057 gmzzu = hhzzu + 1.0d0
00058 gmyxu = hhxyu
00059 gmzxu = hhxzu
00060 gmzyu = hhyzu
00061
00062 rxx = rab(irg,itg,ipg,1)
00063 rxy = rab(irg,itg,ipg,2)
00064 rxz = rab(irg,itg,ipg,3)
00065 ryy = rab(irg,itg,ipg,4)
00066 ryz = rab(irg,itg,ipg,5)
00067 rzz = rab(irg,itg,ipg,6)
00068 ryx = rxy
00069 rzx = rxz
00070 rzy = ryz
00071
00072 rics = 0.0d0
00073 gcdp = 0.0d0
00074 hd2p = 0.0d0
00075
00076 rics = gmxxu*rxx + gmxyu*rxy + gmxzu*rxz &
00077 & + gmyxu*ryx + gmyyu*ryy + gmyzu*ryz &
00078 & + gmzxu*rzx + gmzyu*rzy + gmzzu*rzz
00079
00080 gcdp = gmcrix(irg,itg,ipg)*grad(irg,itg,ipg,1) &
00081 & + gmcriy(irg,itg,ipg)*grad(irg,itg,ipg,2) &
00082 & + gmcriz(irg,itg,ipg)*grad(irg,itg,ipg,3)
00083
00084
00085 call dadbscalar_type3(alps,dabfnc,irg,itg,ipg)
00086 hd2p = hhxxu* dabfnc(1,1) &
00087 & + hhxyu*(dabfnc(1,2) + dabfnc(2,1)) &
00088 & + hhxzu*(dabfnc(1,3) + dabfnc(3,1)) &
00089 & + hhyyu* dabfnc(2,2) &
00090 & + hhyzu*(dabfnc(2,3) + dabfnc(3,2)) &
00091 & + hhzzu* dabfnc(3,3)
00092
00093 sou(irg,itg,ipg) = - hd2p + gcdp + 0.125d0*alpsigc*rics
00094
00095 end do
00096 end do
00097 end do
00098
00099 deallocate(grad)
00100 deallocate(dfdx)
00101 deallocate(dfdy)
00102 deallocate(dfdz)
00103
00104 end subroutine sourceterm_trG_WL
00105
