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