00001 subroutine sourceterm_MoC_WL(souvec)
00002
00003 use phys_constant, only : pi
00004 use def_matter_parameter, only : ome
00005 use def_metric, only : psi, alph, tfkij, bvxu, bvyu, bvzu, bvxd, bvyd, bvzd
00006 use def_matter, only : emdg
00007 use grid_parameter, only : nrg, ntg, npg
00008 use coordinate_grav_r, only : rg
00009 use def_cristoffel_grid, only : cri_grid
00010 use def_gamma_crist, only : gmcrix, gmcriy, gmcriz
00011 use def_ricci_tensor, only : rab
00012 use def_formulation, only : swflu
00013 use def_shift_derivatives, only : cdbvxd, cdbvyd, cdbvzd, &
00014 & pdbvxd, pdbvyd, pdbvzd
00015 use def_Lie_derivatives, only : elpxx, elpxy, elpxz, elpyy, elpyz, elpzz
00016 use def_Lie_derivatives_grid, only : elpxx_grid, elpxy_grid, elpxz_grid, &
00017 & elpyy_grid, elpyz_grid, elpzz_grid
00018 use def_cristoffel, only : cri
00019 use def_cutsw, only : cutfac
00020 use def_metric_hij, only : hxxu, hxyu, hxzu, hyyu, hyzu, hzzu
00021 use def_metric_rotshift, only : ovxd, ovyd, ovzd
00022 use def_formulation, only : chgra
00023 use interface_interpo_linear_type0
00024 use interface_grgrad_midpoint
00025 use interface_grgrad1g_midpoint
00026 use interface_grgrad_4th_gridpoint
00027 use interface_dadbscalar_type0
00028 use interface_dadbscalar_type3
00029 use make_array_3d
00030 implicit none
00031
00032 real(8), pointer :: souvec(:,:,:,:)
00033 real(8), pointer :: grad2x(:,:,:), grad2y(:,:,:), grad2z(:,:,:)
00034 real(8), pointer :: grad4x(:,:,:), grad4y(:,:,:), grad4z(:,:,:)
00035 real(8), pointer :: fnc2(:,:,:), fnc4(:,:,:)
00036 real(8) :: dabfnc(1:3,1:3)
00037
00038 real(8), pointer :: cri1be(:,:,:), cri2be(:,:,:), cri3be(:,:,:),
00039 cri4be(:,:,:), cri5be(:,:,:), cri6be(:,:,:)
00040 real(8) :: gmxxu, gmxyu, gmxzu, gmyyu, gmyzu, gmzzu,
00041 gmyxu, gmzxu, gmzyu, &
00042 hhxxu, hhxyu, hhxzu, hhyxu, hhyyu, hhyzu, &
00043 hhzxu, hhzyu, hhzzu, &
00044 pdbvxdx, pdbvxdy, pdbvxdz, pdbvydx, pdbvydy, pdbvydz, &
00045 pdbvzdx, pdbvzdy, pdbvzdz, &
00046 bvxdc, bvxgc, bvydc, bvygc, bvzdc, bvzgc, &
00047 c11, c12, c13, c14, c15, c16, &
00048 c21, c22, c23, c24, c25, c26, &
00049 c31, c32, c33, c34, c35, c36, &
00050 gc1, gc2, gc3, gclp, gcrdb1, gcrdb2, gdcb, &
00051 gmxbvx, gmxbvy, gmxbvz, gmybvx, gmybvy, gmybvz, &
00052 gmzbvx, gmzbvy, gmzbvz, &
00053 gradfnc4, hddb, pdivlp, rax, ray, raz, tfkax, tfkay, tfkaz, &
00054 san, san2, zfac, &
00055 cutoff, afnc2inv, divlp, dxafn, &
00056 dxc1b, dxc2b, dxc3b, dxc4b, dxc5b, dxc6b, &
00057 dyc1b, dyc2b, dyc3b, dyc4b, dyc5b, dyc6b, &
00058 dzc1b, dzc2b, dzc3b, dzc4b, dzc5b, dzc6b, &
00059 dbvxdx, dbvxdy, dbvxdz, dbvydx, dbvydy, dbvydz, &
00060 dbvzdx, dbvzdy, dbvzdz, &
00061 dexxdx, dexxdy, dexxdz, dexydx, dexydy, dexydz, &
00062 dexzdx, dexzdy, dexzdz, &
00063 deyydx, deyydy, deyydz, deyzdx, deyzdy, deyzdz, &
00064 dezzdx, dezzdy, dezzdz, &
00065 dyafn, dzafn, elpxxd, elpxxm, elpxyd, elpxym, elpxzd, elpxzm, &
00066 elpyxd, elpyxm, elpyyd, elpyym, elpyzd, elpyzm, &
00067 elpzxd, elpzxm, elpzyd, elpzym, elpzzd, elpzzm, &
00068 grad(1:3), ene, fnc2gc, &
00069 bvxdgc, bvydgc, bvzdgc
00070 integer :: ipg, irg, itg, ic1, ic2, ic3, ii, ic0(1:3)
00071
00072 call alloc_array3d(fnc2,0,nrg,0,ntg,0,npg)
00073 call alloc_array3d(fnc4,0,nrg,0,ntg,0,npg)
00074 call alloc_array3d(grad2x,0,nrg,0,ntg,0,npg)
00075 call alloc_array3d(grad2y,0,nrg,0,ntg,0,npg)
00076 call alloc_array3d(grad2z,0,nrg,0,ntg,0,npg)
00077 call alloc_array3d(grad4x,0,nrg,0,ntg,0,npg)
00078 call alloc_array3d(grad4y,0,nrg,0,ntg,0,npg)
00079 call alloc_array3d(grad4z,0,nrg,0,ntg,0,npg)
00080 call alloc_array3d(cri1be,0,nrg,0,ntg,0,npg)
00081 call alloc_array3d(cri2be,0,nrg,0,ntg,0,npg)
00082 call alloc_array3d(cri3be,0,nrg,0,ntg,0,npg)
00083 call alloc_array3d(cri4be,0,nrg,0,ntg,0,npg)
00084 call alloc_array3d(cri5be,0,nrg,0,ntg,0,npg)
00085 call alloc_array3d(cri6be,0,nrg,0,ntg,0,npg)
00086
00087
00088
00089 do ipg = 0, npg
00090 do itg = 0, ntg
00091 do irg = 0, nrg
00092 hhxxu=hxxu(irg,itg,ipg)
00093 hhxyu=hxyu(irg,itg,ipg)
00094 hhxzu=hxzu(irg,itg,ipg)
00095 hhyyu=hyyu(irg,itg,ipg)
00096 hhyzu=hyzu(irg,itg,ipg)
00097 hhzzu=hzzu(irg,itg,ipg)
00098 hhyxu = hhxyu
00099 hhzxu = hhxzu
00100 hhzyu = hhyzu
00101
00102 call grgrad_4th_gridpoint(bvxd,pdbvxdx,pdbvxdy,pdbvxdz,irg,itg,ipg)
00103 call grgrad_4th_gridpoint(bvyd,pdbvydx,pdbvydy,pdbvydz,irg,itg,ipg)
00104 call grgrad_4th_gridpoint(bvzd,pdbvzdx,pdbvzdy,pdbvzdz,irg,itg,ipg)
00105
00106 fnc4(irg,itg,ipg) = hhxxu*pdbvxdx+hhxyu*pdbvxdy+hhxzu*pdbvxdz &
00107 & + hhyxu*pdbvydx+hhyyu*pdbvydy+hhyzu*pdbvydz &
00108 & + hhzxu*pdbvzdx+hhzyu*pdbvzdy+hhzzu*pdbvzdz
00109
00110 bvxdc=bvxd(irg,itg,ipg)
00111 bvydc=bvyd(irg,itg,ipg)
00112 bvzdc=bvzd(irg,itg,ipg)
00113 c11 = cri_grid(irg,itg,ipg,1,1)
00114 c12 = cri_grid(irg,itg,ipg,1,2)
00115 c13 = cri_grid(irg,itg,ipg,1,3)
00116 c14 = cri_grid(irg,itg,ipg,1,4)
00117 c15 = cri_grid(irg,itg,ipg,1,5)
00118 c16 = cri_grid(irg,itg,ipg,1,6)
00119 c21 = cri_grid(irg,itg,ipg,2,1)
00120 c22 = cri_grid(irg,itg,ipg,2,2)
00121 c23 = cri_grid(irg,itg,ipg,2,3)
00122 c24 = cri_grid(irg,itg,ipg,2,4)
00123 c25 = cri_grid(irg,itg,ipg,2,5)
00124 c26 = cri_grid(irg,itg,ipg,2,6)
00125 c31 = cri_grid(irg,itg,ipg,3,1)
00126 c32 = cri_grid(irg,itg,ipg,3,2)
00127 c33 = cri_grid(irg,itg,ipg,3,3)
00128 c34 = cri_grid(irg,itg,ipg,3,4)
00129 c35 = cri_grid(irg,itg,ipg,3,5)
00130 c36 = cri_grid(irg,itg,ipg,3,6)
00131 cri1be(irg,itg,ipg) = c11*bvxdc + c21*bvydc + c31*bvzdc
00132 cri2be(irg,itg,ipg) = c12*bvxdc + c22*bvydc + c32*bvzdc
00133 cri3be(irg,itg,ipg) = c13*bvxdc + c23*bvydc + c33*bvzdc
00134 cri4be(irg,itg,ipg) = c14*bvxdc + c24*bvydc + c34*bvzdc
00135 cri5be(irg,itg,ipg) = c15*bvxdc + c25*bvydc + c35*bvzdc
00136 cri6be(irg,itg,ipg) = c16*bvxdc + c26*bvydc + c36*bvzdc
00137
00138 end do
00139 end do
00140 end do
00141
00142 fnc2(0:nrg,0:ntg,0:npg) = psi(0:nrg,0:ntg,0:npg)**6/alph(0:nrg,0:ntg,0:npg)
00143 call grgrad_midpoint(fnc2,grad2x,grad2y,grad2z)
00144 call grgrad_midpoint(fnc4,grad4x,grad4y,grad4z)
00145
00146 san = 1.0d0/3.0d0
00147 san2= 2.0d0/3.0d0
00148
00149
00150
00151 do ii = 1, 3
00152
00153 do ipg = 1, npg
00154 do itg = 1, ntg
00155 do irg = 1, nrg
00156
00157 cutoff = 1.0d0
00158 if (chgra == 'c'.or.chgra == 'C'.or.chgra == 'W') then
00159 if (rg(irg) > cutfac*pi/ome) cutoff = 0.0d0
00160 end if
00161
00162 call interpo_linear_type0(bvxgc,bvxu,irg,itg,ipg)
00163 call interpo_linear_type0(bvygc,bvyu,irg,itg,ipg)
00164 call interpo_linear_type0(bvzgc,bvzu,irg,itg,ipg)
00165
00166 call interpo_linear_type0(hhxxu,hxxu,irg,itg,ipg)
00167 call interpo_linear_type0(hhxyu,hxyu,irg,itg,ipg)
00168 call interpo_linear_type0(hhxzu,hxzu,irg,itg,ipg)
00169 call interpo_linear_type0(hhyyu,hyyu,irg,itg,ipg)
00170 call interpo_linear_type0(hhyzu,hyzu,irg,itg,ipg)
00171 call interpo_linear_type0(hhzzu,hzzu,irg,itg,ipg)
00172 hhyxu = hhxyu
00173 hhzxu = hhxzu
00174 hhzyu = hhyzu
00175 gmxxu = hhxxu + 1.0d0
00176 gmxyu = hhxyu
00177 gmxzu = hhxzu
00178 gmyyu = hhyyu + 1.0d0
00179 gmyzu = hhyzu
00180 gmzzu = hhzzu + 1.0d0
00181 gmyxu = gmxyu
00182 gmzxu = gmxzu
00183 gmzyu = gmyzu
00184
00185 dxafn = grad2x(irg,itg,ipg)
00186 dyafn = grad2y(irg,itg,ipg)
00187 dzafn = grad2z(irg,itg,ipg)
00188 tfkax = tfkij(irg,itg,ipg,ii,1)
00189 tfkay = tfkij(irg,itg,ipg,ii,2)
00190 tfkaz = tfkij(irg,itg,ipg,ii,3)
00191
00192 dbvxdx = cdbvxd(irg,itg,ipg,1)
00193 dbvydx = cdbvyd(irg,itg,ipg,1)
00194 dbvzdx = cdbvzd(irg,itg,ipg,1)
00195 dbvxdy = cdbvxd(irg,itg,ipg,2)
00196 dbvydy = cdbvyd(irg,itg,ipg,2)
00197 dbvzdy = cdbvzd(irg,itg,ipg,2)
00198 dbvxdz = cdbvxd(irg,itg,ipg,3)
00199 dbvydz = cdbvyd(irg,itg,ipg,3)
00200 dbvzdz = cdbvzd(irg,itg,ipg,3)
00201
00202 hddb = 0.0d0
00203 gdcb = 0.0d0
00204 gcrdb1 = 0.0d0
00205 gcrdb2 = 0.0d0
00206 gradfnc4 = 0.0d0
00207 rax = 0.0d0
00208 ray = 0.0d0
00209 raz = 0.0d0
00210 divlp = 0.0d0
00211
00212 if (chgra /= 'i') then
00213
00214 if (ii == 1) gradfnc4 = grad4x(irg,itg,ipg)
00215 if (ii == 2) gradfnc4 = grad4y(irg,itg,ipg)
00216 if (ii == 3) gradfnc4 = grad4z(irg,itg,ipg)
00217
00218 c11 = cri(irg,itg,ipg,1,1)
00219 c12 = cri(irg,itg,ipg,1,2)
00220 c13 = cri(irg,itg,ipg,1,3)
00221 c14 = cri(irg,itg,ipg,1,4)
00222 c15 = cri(irg,itg,ipg,1,5)
00223 c16 = cri(irg,itg,ipg,1,6)
00224 c21 = cri(irg,itg,ipg,2,1)
00225 c22 = cri(irg,itg,ipg,2,2)
00226 c23 = cri(irg,itg,ipg,2,3)
00227 c24 = cri(irg,itg,ipg,2,4)
00228 c25 = cri(irg,itg,ipg,2,5)
00229 c26 = cri(irg,itg,ipg,2,6)
00230 c31 = cri(irg,itg,ipg,3,1)
00231 c32 = cri(irg,itg,ipg,3,2)
00232 c33 = cri(irg,itg,ipg,3,3)
00233 c34 = cri(irg,itg,ipg,3,4)
00234 c35 = cri(irg,itg,ipg,3,5)
00235 c36 = cri(irg,itg,ipg,3,6)
00236
00237 gc1 = gmcrix(irg,itg,ipg)
00238 gc2 = gmcriy(irg,itg,ipg)
00239 gc3 = gmcriz(irg,itg,ipg)
00240
00241 ic1 = 1*(ii-2)*(ii-3)/2 - 2*(ii-1)*(ii-3) + 3*(ii-1)*(ii-2)/2
00242 ic2 = 2*(ii-2)*(ii-3)/2 - 4*(ii-1)*(ii-3) + 5*(ii-1)*(ii-2)/2
00243 ic3 = 3*(ii-2)*(ii-3)/2 - 5*(ii-1)*(ii-3) + 6*(ii-1)*(ii-2)/2
00244
00245 rax = rab(irg,itg,ipg,ic1)
00246 ray = rab(irg,itg,ipg,ic2)
00247 raz = rab(irg,itg,ipg,ic3)
00248
00249 ic0(1) = ic1
00250 ic0(2) = ic2
00251 ic0(3) = ic3
00252
00253
00254
00255
00256
00257
00258 if(ii == 1) call dadbscalar_type3(bvxd,dabfnc,irg,itg,ipg)
00259 if(ii == 2) call dadbscalar_type3(bvyd,dabfnc,irg,itg,ipg)
00260 if(ii == 3) call dadbscalar_type3(bvzd,dabfnc,irg,itg,ipg)
00261 hddb = hhxxu*dabfnc(1,1) + hhxyu*dabfnc(1,2) + hhxzu*dabfnc(1,3) &
00262 & + hhyxu*dabfnc(2,1) + hhyyu*dabfnc(2,2) + hhyzu*dabfnc(2,3) &
00263 & + hhzxu*dabfnc(3,1) + hhzyu*dabfnc(3,2) + hhzzu*dabfnc(3,3)
00264
00265
00266
00267 if (ii == 1) then
00268 call grgrad1g_midpoint(cri1be,grad,irg,itg,ipg)
00269 dxc1b = grad(1)
00270 dyc1b = grad(2)
00271 dzc1b = grad(3)
00272 call grgrad1g_midpoint(cri2be,grad,irg,itg,ipg)
00273 dxc2b = grad(1)
00274 dyc2b = grad(2)
00275 dzc2b = grad(3)
00276 call grgrad1g_midpoint(cri3be,grad,irg,itg,ipg)
00277 dxc3b = grad(1)
00278 dyc3b = grad(2)
00279 dzc3b = grad(3)
00280 gdcb = gmxxu*dxc1b + gmxyu*dxc2b + gmxzu*dxc3b &
00281 & + gmyxu*dyc1b + gmyyu*dyc2b + gmyzu*dyc3b &
00282 & + gmzxu*dzc1b + gmzyu*dzc2b + gmzzu*dzc3b
00283 end if
00284 if (ii == 2) then
00285 call grgrad1g_midpoint(cri2be,grad,irg,itg,ipg)
00286 dxc2b = grad(1)
00287 dyc2b = grad(2)
00288 dzc2b = grad(3)
00289 call grgrad1g_midpoint(cri4be,grad,irg,itg,ipg)
00290 dxc4b = grad(1)
00291 dyc4b = grad(2)
00292 dzc4b = grad(3)
00293 call grgrad1g_midpoint(cri5be,grad,irg,itg,ipg)
00294 dxc5b = grad(1)
00295 dyc5b = grad(2)
00296 dzc5b = grad(3)
00297 gdcb = gmxxu*dxc2b + gmxyu*dxc4b + gmxzu*dxc5b &
00298 & + gmyxu*dyc2b + gmyyu*dyc4b + gmyzu*dyc5b &
00299 & + gmzxu*dzc2b + gmzyu*dzc4b + gmzzu*dzc5b
00300 end if
00301 if (ii == 3) then
00302 call grgrad1g_midpoint(cri3be,grad,irg,itg,ipg)
00303 dxc3b = grad(1)
00304 dyc3b = grad(2)
00305 dzc3b = grad(3)
00306 call grgrad1g_midpoint(cri5be,grad,irg,itg,ipg)
00307 dxc5b = grad(1)
00308 dyc5b = grad(2)
00309 dzc5b = grad(3)
00310 call grgrad1g_midpoint(cri6be,grad,irg,itg,ipg)
00311 dxc6b = grad(1)
00312 dyc6b = grad(2)
00313 dzc6b = grad(3)
00314 gdcb = gmxxu*dxc3b + gmxyu*dxc5b + gmxzu*dxc6b &
00315 & + gmyxu*dyc3b + gmyyu*dyc5b + gmyzu*dyc6b &
00316 & + gmzxu*dzc3b + gmzyu*dzc5b + gmzzu*dzc6b
00317 end if
00318
00319
00320
00321 if (ii == 1) gcrdb1 = gc1*dbvxdx + gc2*dbvxdy + gc3*dbvxdz
00322 if (ii == 2) gcrdb1 = gc1*dbvydx + gc2*dbvydy + gc3*dbvydz
00323 if (ii == 3) gcrdb1 = gc1*dbvzdx + gc2*dbvzdy + gc3*dbvzdz
00324
00325 gmxbvx = gmxxu*dbvxdx + gmxyu*dbvxdy + gmxzu*dbvxdz
00326 gmxbvy = gmxxu*dbvydx + gmxyu*dbvydy + gmxzu*dbvydz
00327 gmxbvz = gmxxu*dbvzdx + gmxyu*dbvzdy + gmxzu*dbvzdz
00328 gmybvx = gmyxu*dbvxdx + gmyyu*dbvxdy + gmyzu*dbvxdz
00329 gmybvy = gmyxu*dbvydx + gmyyu*dbvydy + gmyzu*dbvydz
00330 gmybvz = gmyxu*dbvzdx + gmyyu*dbvzdy + gmyzu*dbvzdz
00331 gmzbvx = gmzxu*dbvxdx + gmzyu*dbvxdy + gmzzu*dbvxdz
00332 gmzbvy = gmzxu*dbvydx + gmzyu*dbvydy + gmzzu*dbvydz
00333 gmzbvz = gmzxu*dbvzdx + gmzyu*dbvzdy + gmzzu*dbvzdz
00334
00335 if (ii == 1) &
00336 & gcrdb2 = c11*gmxbvx + c21*gmxbvy + c31*gmxbvz &
00337 & + c12*gmybvx + c22*gmybvy + c32*gmybvz &
00338 & + c13*gmzbvx + c23*gmzbvy + c33*gmzbvz
00339
00340 if (ii == 2) &
00341 & gcrdb2 = c12*gmxbvx + c22*gmxbvy + c32*gmxbvz &
00342 & + c14*gmybvx + c24*gmybvy + c34*gmybvz &
00343 & + c15*gmzbvx + c25*gmzbvy + c35*gmzbvz
00344
00345 if (ii == 3) &
00346 & gcrdb2 = c13*gmxbvx + c23*gmxbvy + c33*gmxbvz &
00347 & + c15*gmybvx + c25*gmybvy + c35*gmybvz &
00348 & + c16*gmzbvx + c26*gmzbvy + c36*gmzbvz
00349
00350
00351
00352 if (chgra == 'h'.or.chgra == 'c'.or.chgra == 'C' &
00353 &.or.chgra == 'H'.or.chgra == 'W') then
00354
00355 elpxxd = elpxx(irg,itg,ipg)
00356 elpxyd = elpxy(irg,itg,ipg)
00357 elpxzd = elpxz(irg,itg,ipg)
00358 elpyyd = elpyy(irg,itg,ipg)
00359 elpyzd = elpyz(irg,itg,ipg)
00360 elpzzd = elpzz(irg,itg,ipg)
00361 elpyxd = elpxyd
00362 elpzxd = elpxzd
00363 elpzyd = elpyzd
00364
00365 elpxxm = gmxxu*elpxxd + gmxyu*elpxyd + gmxzu*elpxzd
00366 elpxym = gmyxu*elpxxd + gmyyu*elpxyd + gmyzu*elpxzd
00367 elpxzm = gmzxu*elpxxd + gmzyu*elpxyd + gmzzu*elpxzd
00368 elpyxm = gmxxu*elpyxd + gmxyu*elpyyd + gmxzu*elpyzd
00369 elpyym = gmyxu*elpyxd + gmyyu*elpyyd + gmyzu*elpyzd
00370 elpyzm = gmzxu*elpyxd + gmzyu*elpyyd + gmzzu*elpyzd
00371 elpzxm = gmxxu*elpzxd + gmxyu*elpzyd + gmxzu*elpzzd
00372 elpzym = gmyxu*elpzxd + gmyyu*elpzyd + gmyzu*elpzzd
00373 elpzzm = gmzxu*elpzxd + gmzyu*elpzyd + gmzzu*elpzzd
00374
00375 if (ii == 1) then
00376 call grgrad1g_midpoint(elpxx_grid,grad,irg,itg,ipg)
00377 dexxdx = grad(1)
00378 dexxdy = grad(2)
00379 dexxdz = grad(3)
00380 call grgrad1g_midpoint(elpxy_grid,grad,irg,itg,ipg)
00381 dexydx = grad(1)
00382 dexydy = grad(2)
00383 dexydz = grad(3)
00384 call grgrad1g_midpoint(elpxz_grid,grad,irg,itg,ipg)
00385 dexzdx = grad(1)
00386 dexzdy = grad(2)
00387 dexzdz = grad(3)
00388 pdivlp = gmxxu*dexxdx + gmxyu*dexxdy + gmxzu*dexxdz &
00389 & + gmyxu*dexydx + gmyyu*dexydy + gmyzu*dexydz &
00390 & + gmzxu*dexzdx + gmzyu*dexzdy + gmzzu*dexzdz
00391 gclp = c11*elpxxm + c12*elpxym + c13*elpxzm &
00392 & + c21*elpyxm + c22*elpyym + c23*elpyzm &
00393 & + c31*elpzxm + c32*elpzym + c33*elpzzm
00394 divlp = pdivlp - gclp
00395 end if
00396 if (ii == 2) then
00397 call grgrad1g_midpoint(elpxy_grid,grad,irg,itg,ipg)
00398 dexydx = grad(1)
00399 dexydy = grad(2)
00400 dexydz = grad(3)
00401 call grgrad1g_midpoint(elpyy_grid,grad,irg,itg,ipg)
00402 deyydx = grad(1)
00403 deyydy = grad(2)
00404 deyydz = grad(3)
00405 call grgrad1g_midpoint(elpyz_grid,grad,irg,itg,ipg)
00406 deyzdx = grad(1)
00407 deyzdy = grad(2)
00408 deyzdz = grad(3)
00409 pdivlp = gmxxu*dexydx + gmxyu*dexydy + gmxzu*dexydz &
00410 & + gmyxu*deyydx + gmyyu*deyydy + gmyzu*deyydz &
00411 & + gmzxu*deyzdx + gmzyu*deyzdy + gmzzu*deyzdz
00412 gclp = c12*elpxxm + c14*elpxym + c15*elpxzm &
00413 & + c22*elpyxm + c24*elpyym + c25*elpyzm &
00414 & + c32*elpzxm + c34*elpzym + c35*elpzzm
00415 divlp = pdivlp - gclp
00416 end if
00417 if (ii == 3) then
00418 call grgrad1g_midpoint(elpxz_grid,grad,irg,itg,ipg)
00419 dexzdx = grad(1)
00420 dexzdy = grad(2)
00421 dexzdz = grad(3)
00422 call grgrad1g_midpoint(elpyz_grid,grad,irg,itg,ipg)
00423 deyzdx = grad(1)
00424 deyzdy = grad(2)
00425 deyzdz = grad(3)
00426 call grgrad1g_midpoint(elpzz_grid,grad,irg,itg,ipg)
00427 dezzdx = grad(1)
00428 dezzdy = grad(2)
00429 dezzdz = grad(3)
00430 pdivlp = gmxxu*dexzdx + gmxyu*dexzdy + gmxzu*dexzdz &
00431 & + gmyxu*deyzdx + gmyyu*deyzdy + gmyzu*deyzdz &
00432 & + gmzxu*dezzdx + gmzyu*dezzdy + gmzzu*dezzdz
00433 gclp = c13*elpxxm + c15*elpxym + c16*elpxzm &
00434 & + c23*elpyxm + c25*elpyym + c26*elpyzm &
00435 & + c33*elpzxm + c35*elpzym + c36*elpzzm
00436 divlp = pdivlp - gclp
00437 end if
00438
00439 end if
00440
00441 end if
00442
00443 souvec(irg,itg,ipg,ii) = &
00444 & - hddb + gdcb + gcrdb1 + gcrdb2 - san*gradfnc4 &
00445 & -(rax*bvxgc + ray*bvygc + raz*bvzgc) - ome*divlp*cutoff &
00446 & - 2.0d0*afnc2inv*( &
00447 & + hhxxu*tfkax*dxafn + hhxyu*tfkax*dyafn + hhxzu*tfkax*dzafn &
00448 & + hhyxu*tfkay*dxafn + hhyyu*tfkay*dyafn + hhyzu*tfkay*dzafn &
00449 & + hhzxu*tfkaz*dxafn + hhzyu*tfkaz*dyafn + hhzzu*tfkaz*dzafn)
00450
00451 end do
00452 end do
00453 end do
00454 end do
00455
00456
00457
00458 deallocate(fnc4)
00459 deallocate(fnc2)
00460 deallocate(grad2x)
00461 deallocate(grad2y)
00462 deallocate(grad2z)
00463 deallocate(grad4x)
00464 deallocate(grad4y)
00465 deallocate(grad4z)
00466 deallocate(cri1be)
00467 deallocate(cri2be)
00468 deallocate(cri3be)
00469 deallocate(cri4be)
00470 deallocate(cri5be)
00471 deallocate(cri6be)
00472
00473 end subroutine sourceterm_MoC_WL
