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