00001 subroutine sourceterm_MWtemp_WL(sou)
00002 use def_metric, only : psi
00003 use grid_parameter, only : nrg, ntg, npg
00004 use def_gamma_crist, only : gmcrix, gmcriy, gmcriz
00005 use def_metric, only : psi, alph
00006 use def_metric_hij, only : hxxu, hxyu, hxzu, hyyu, hyzu, hzzu
00007 use def_metric_hij_dirac, only : dagmabu
00008 use def_emfield, only : alva
00009 use def_emfield_derivatives, only : Lie_bAxd, Lie_bAyd, Lie_bAzd, &
00010 & Lie_bAxd_grid, Lie_bAyd_grid, Lie_bAzd_grid
00011 use def_faraday_tensor, only : fxd, fyd, fzd
00012 use make_array_4d
00013 use make_array_3d
00014 use interface_interpo_linear_type0
00015 use interface_grgrad_midpoint
00016 use interface_dadbscalar_type0
00017 use interface_dadbscalar_type3
00018 implicit none
00019
00020 real(long), pointer :: sou(:,:,:)
00021 real(long), pointer :: fnc2(:,:,:)
00022 real(long), pointer :: grad2x(:,:,:), grad2y(:,:,:), grad2z(:,:,:)
00023 real(long), pointer :: dxLbAx(:,:,:), dyLbAx(:,:,:), dzLbAx(:,:,:)
00024 real(long), pointer :: dxLbAy(:,:,:), dyLbAy(:,:,:), dzLbAy(:,:,:)
00025 real(long), pointer :: dxLbAz(:,:,:), dyLbAz(:,:,:), dzLbAz(:,:,:)
00026 real(long), pointer :: dfdx(:,:,:), dfdy(:,:,:), dfdz(:,:,:)
00027 real(long) :: hiju(3,3), dabfnc(3,3), dpafa(3,3), dLbAgc(3,3)
00028 real(long) :: psigc, alphgc, a2divp2, dxalva, dyalva, dzalva,
00029 Lie_bAxdgc, Lie_bAydgc, Lie_bAzdgc, &
00030 diracx, diracy, diracz, zfac, &
00031 hd2va, diracdva, dps2alfa, diracLieA, trdLieA, &
00032 gmxxu, gmxyu, gmxzu, gmyyu, gmyzu, gmzzu, gmyxu, gmzxu, gmzyu, &
00033 hhxxu, hhxyu, hhxzu, hhyyu, hhyzu, hhzzu, hhyxu, hhzxu, hhzyu
00034 integer :: ipg, irg, itg
00035
00036
00037
00038
00039 call alloc_array3d(dfdx,1,nrg,1,ntg,1,npg)
00040 call alloc_array3d(dfdy,1,nrg,1,ntg,1,npg)
00041 call alloc_array3d(dfdz,1,nrg,1,ntg,1,npg)
00042 call alloc_array3d(fnc2,0,nrg,0,ntg,0,npg)
00043 call alloc_array3d(grad2x,1,nrg,1,ntg,1,npg)
00044 call alloc_array3d(grad2y,1,nrg,1,ntg,1,npg)
00045 call alloc_array3d(grad2z,1,nrg,1,ntg,1,npg)
00046 call alloc_array3d(dxLbAx,1,nrg,1,ntg,1,npg)
00047 call alloc_array3d(dyLbAx,1,nrg,1,ntg,1,npg)
00048 call alloc_array3d(dzLbAx,1,nrg,1,ntg,1,npg)
00049 call alloc_array3d(dxLbAy,1,nrg,1,ntg,1,npg)
00050 call alloc_array3d(dyLbAy,1,nrg,1,ntg,1,npg)
00051 call alloc_array3d(dzLbAy,1,nrg,1,ntg,1,npg)
00052 call alloc_array3d(dxLbAz,1,nrg,1,ntg,1,npg)
00053 call alloc_array3d(dyLbAz,1,nrg,1,ntg,1,npg)
00054 call alloc_array3d(dzLbAz,1,nrg,1,ntg,1,npg)
00055 fnc2(0:nrg,0:ntg,0:npg) = psi(0:nrg,0:ntg,0:npg)**2/alph(0:nrg,0:ntg,0:npg)
00056 call grgrad_midpoint(fnc2,grad2x,grad2y,grad2z)
00057 call grgrad_midpoint(Lie_bAxd_grid,dxLbAx,dyLbAx,dzLbAx)
00058 call grgrad_midpoint(Lie_bAyd_grid,dxLbAy,dyLbAy,dzLbAy)
00059 call grgrad_midpoint(Lie_bAzd_grid,dxLbAz,dyLbAz,dzLbAz)
00060 call grgrad_midpoint(alva,dfdx,dfdy,dfdz)
00061
00062 do ipg = 1, npg
00063 do itg = 1, ntg
00064 do irg = 1, nrg
00065
00066 call interpo_linear_type0(psigc, psi,irg,itg,ipg)
00067 call interpo_linear_type0(alphgc,alph,irg,itg,ipg)
00068 call interpo_linear_type0(hhxxu,hxxu,irg,itg,ipg)
00069 call interpo_linear_type0(hhxyu,hxyu,irg,itg,ipg)
00070 call interpo_linear_type0(hhxzu,hxzu,irg,itg,ipg)
00071 call interpo_linear_type0(hhyyu,hyyu,irg,itg,ipg)
00072 call interpo_linear_type0(hhyzu,hyzu,irg,itg,ipg)
00073 call interpo_linear_type0(hhzzu,hzzu,irg,itg,ipg)
00074 hhyxu = hhxyu ; hhzxu = hhxzu ; hhzyu = hhyzu
00075 hiju(1,1) = hhxxu ; hiju(1,2) = hhxyu ; hiju(1,3) = hhxzu
00076 hiju(2,1) = hhxyu ; hiju(2,2) = hhyyu ; hiju(2,3) = hhyzu
00077 hiju(3,1) = hhxzu ; hiju(3,2) = hhyzu ; hiju(3,3) = hhzzu
00078 gmxxu = hhxxu + 1.0d0
00079 gmxyu = hhxyu
00080 gmxzu = hhxzu
00081 gmyyu = hhyyu + 1.0d0
00082 gmyzu = hhyzu
00083 gmzzu = hhzzu + 1.0d0
00084 gmyxu = hhxyu
00085 gmzxu = hhxzu
00086 gmzyu = hhyzu
00087 a2divp2 = alphgc**2/psigc**2
00088 dxalva = dfdx(irg,itg,ipg)
00089 dyalva = dfdy(irg,itg,ipg)
00090 dzalva = dfdz(irg,itg,ipg)
00091 Lie_bAxdgc = Lie_bAxd(irg,itg,ipg)
00092 Lie_bAydgc = Lie_bAyd(irg,itg,ipg)
00093 Lie_bAzdgc = Lie_bAzd(irg,itg,ipg)
00094 diracx = dagmabu(irg,itg,ipg,1)
00095 diracy = dagmabu(irg,itg,ipg,2)
00096 diracz = dagmabu(irg,itg,ipg,3)
00097
00098 dpafa(1,1) = grad2x(irg,itg,ipg)*fxd(irg,itg,ipg)
00099 dpafa(1,2) = grad2x(irg,itg,ipg)*fyd(irg,itg,ipg)
00100 dpafa(1,3) = grad2x(irg,itg,ipg)*fzd(irg,itg,ipg)
00101 dpafa(2,1) = grad2y(irg,itg,ipg)*fxd(irg,itg,ipg)
00102 dpafa(2,2) = grad2y(irg,itg,ipg)*fyd(irg,itg,ipg)
00103 dpafa(2,3) = grad2y(irg,itg,ipg)*fzd(irg,itg,ipg)
00104 dpafa(3,1) = grad2z(irg,itg,ipg)*fxd(irg,itg,ipg)
00105 dpafa(3,2) = grad2z(irg,itg,ipg)*fyd(irg,itg,ipg)
00106 dpafa(3,3) = grad2z(irg,itg,ipg)*fzd(irg,itg,ipg)
00107
00108 dLbAgc(1,1) = dxLbAx(irg,itg,ipg)
00109 dLbAgc(1,2) = dyLbAx(irg,itg,ipg)
00110 dLbAgc(1,3) = dzLbAx(irg,itg,ipg)
00111 dLbAgc(2,1) = dxLbAy(irg,itg,ipg)
00112 dLbAgc(2,2) = dyLbAy(irg,itg,ipg)
00113 dLbAgc(2,3) = dzLbAy(irg,itg,ipg)
00114 dLbAgc(3,1) = dxLbAz(irg,itg,ipg)
00115 dLbAgc(3,2) = dyLbAz(irg,itg,ipg)
00116 dLbAgc(3,3) = dzLbAz(irg,itg,ipg)
00117
00118 diracdva = diracx*dxalva + diracy*dyalva + diracz*dzalva
00119
00120
00121 call dadbscalar_type3(alva,dabfnc,irg,itg,ipg)
00122 call compute_trace(hiju,dabfnc,hd2va)
00123 call compute_trace(hiju,dpafa,dps2alfa)
00124 dps2alfa = a2divp2*dps2alfa
00125 call compute_trace(hiju,dLbAgc,trdLieA)
00126
00127 sou(irg,itg,ipg) = - hd2va - diracdva &
00128 & + dps2alfa + trdLieA
00129
00130
00131 end do
00132 end do
00133 end do
00134
00135 deallocate(dfdx)
00136 deallocate(dfdy)
00137 deallocate(dfdz)
00138 deallocate(fnc2)
00139 deallocate(grad2x)
00140 deallocate(grad2y)
00141 deallocate(grad2z)
00142 deallocate(dxLbAx)
00143 deallocate(dyLbAx)
00144 deallocate(dzLbAx)
00145 deallocate(dxLbAy)
00146 deallocate(dyLbAy)
00147 deallocate(dzLbAy)
00148 deallocate(dxLbAz)
00149 deallocate(dyLbAz)
00150 deallocate(dzLbAz)
00151
00152 end subroutine sourceterm_MWtemp_WL
