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