00001 subroutine excurve_WL_gridpoint
00002 use phys_constant, only : nnrg, nntg, nnpg, pi
00003 use grid_parameter, only : nrg, ntg, npg
00004 use def_matter_parameter, only : ome, ber, radi
00005 use def_metric, only : alph, bvxd, bvyd, bvzd, bvxu, bvyu, bvzu
00006 use def_metric_excurve_grid, only : tfkij_grid, tfkijkij_grid, trk_grid
00007 use coordinate_grav_r, only : rg
00008
00009 use def_cristoffel_grid, only : cri_grid
00010 use def_Lie_derivatives_grid, only : elpxx_grid, elpxy_grid, elpxz_grid, &
00011 & elpyy_grid, elpyz_grid, elpzz_grid, &
00012 & rlbxx_grid, rlbxy_grid, rlbxz_grid, &
00013 & rlbyy_grid, rlbyz_grid, rlbzz_grid
00014 use def_shift_derivatives_grid, only : cdbvxd_grid, cdbvyd_grid, &
00015 & cdbvzd_grid, cdivbv_grid, &
00016 & pdbvxd_grid, pdbvyd_grid, pdbvzd_grid
00017 use def_metric_hij, only : hxxd, hxyd, hxzd, hyyd, hyzd, hzzd, &
00018 & hxxu, hxyu, hxzu, hyyu, hyzu, hzzu
00019 use def_cutsw, only : cutfac
00020 use def_formulation, only : chgra
00021 use def_vector_x, only : vec_xg
00022
00023 use interface_grgrad_4th_gridpoint
00024 use make_array_3d
00025 use make_array_4d
00026 implicit none
00027
00028 real(8) :: dfdx, dfdy, dfdz
00029 real(8) :: gamu(3,3)
00030 real(8) :: ainvh, bvxha, bvyha, bvzha, cutoff, diver, fa23,
00031 dbvxdx, dbvxdy, dbvxdz, &
00032 dbvydx, dbvydy, dbvydz, &
00033 dbvzdx, dbvzdy, dbvzdz, &
00034 gmxxd, gmxyd, gmxzd, gmxxu, gmxyu, gmxzu, &
00035 gmyxd, gmyyd, gmyzd, gmyxu, gmyyu, gmyzu, &
00036 gmzxd, gmzyd, gmzzd, gmzxu, gmzyu, gmzzu, &
00037 oelpxx, oelpxy, oelpxz, &
00038 oelpyy, oelpyz, oelpzz, &
00039 pdbvxdx, pdbvxdy, pdbvxdz, &
00040 pdbvydx, pdbvydy, pdbvydz, &
00041 pdbvzdx, pdbvzdy, pdbvzdz
00042 integer :: ia, ib, ic, id, info, ipg, irg, itg
00043
00044
00045
00046
00047
00048 info = 0
00049
00050
00051 fa23 = 2.0d0/3.0d0
00052
00053 do ipg = 0, npg
00054 do itg = 0, ntg
00055 do irg = 0, nrg
00056
00057 call grgrad_4th_gridpoint(bvxd,dfdx,dfdy,dfdz,irg,itg,ipg)
00058 pdbvxd_grid(irg,itg,ipg,1) = dfdx
00059 pdbvxd_grid(irg,itg,ipg,2) = dfdy
00060 pdbvxd_grid(irg,itg,ipg,3) = dfdz
00061
00062 call grgrad_4th_gridpoint(bvyd,dfdx,dfdy,dfdz,irg,itg,ipg)
00063 pdbvyd_grid(irg,itg,ipg,1) = dfdx
00064 pdbvyd_grid(irg,itg,ipg,2) = dfdy
00065 pdbvyd_grid(irg,itg,ipg,3) = dfdz
00066
00067 call grgrad_4th_gridpoint(bvzd,dfdx,dfdy,dfdz,irg,itg,ipg)
00068 pdbvzd_grid(irg,itg,ipg,1) = dfdx
00069 pdbvzd_grid(irg,itg,ipg,2) = dfdy
00070 pdbvzd_grid(irg,itg,ipg,3) = dfdz
00071
00072 call grgrad_4th_gridpoint(bvxu,pdbvxdx,pdbvxdy,pdbvxdz,irg,itg,ipg)
00073 call grgrad_4th_gridpoint(bvyu,pdbvydx,pdbvydy,pdbvydz,irg,itg,ipg)
00074 call grgrad_4th_gridpoint(bvzu,pdbvzdx,pdbvzdy,pdbvzdz,irg,itg,ipg)
00075
00076 cutoff = 1.0d0
00077 if (chgra == 'c'.or.chgra == 'C'.or.chgra == 'W') then
00078 if (rg(irg) > cutfac*pi/ome) cutoff = 0.0d0
00079 end if
00080
00081 bvxha = bvxd(irg,itg,ipg)
00082 bvyha = bvyd(irg,itg,ipg)
00083 bvzha = bvzd(irg,itg,ipg)
00084 ainvh = alph(irg,itg,ipg)
00085 ainvh = 0.5d0/ainvh
00086
00087 cdbvxd_grid(irg,itg,ipg,1) = pdbvxd_grid(irg,itg,ipg,1) &
00088 & - cri_grid(irg,itg,ipg,1,1)*bvxha &
00089 & - cri_grid(irg,itg,ipg,2,1)*bvyha &
00090 & - cri_grid(irg,itg,ipg,3,1)*bvzha
00091 cdbvyd_grid(irg,itg,ipg,1) = pdbvyd_grid(irg,itg,ipg,1) &
00092 & - cri_grid(irg,itg,ipg,1,2)*bvxha &
00093 & - cri_grid(irg,itg,ipg,2,2)*bvyha &
00094 & - cri_grid(irg,itg,ipg,3,2)*bvzha
00095 cdbvzd_grid(irg,itg,ipg,1) = pdbvzd_grid(irg,itg,ipg,1) &
00096 & - cri_grid(irg,itg,ipg,1,3)*bvxha &
00097 & - cri_grid(irg,itg,ipg,2,3)*bvyha &
00098 & - cri_grid(irg,itg,ipg,3,3)*bvzha
00099
00100 cdbvxd_grid(irg,itg,ipg,2) = pdbvxd_grid(irg,itg,ipg,2) &
00101 & - cri_grid(irg,itg,ipg,1,2)*bvxha &
00102 & - cri_grid(irg,itg,ipg,2,2)*bvyha &
00103 & - cri_grid(irg,itg,ipg,3,2)*bvzha
00104 cdbvyd_grid(irg,itg,ipg,2) = pdbvyd_grid(irg,itg,ipg,2) &
00105 & - cri_grid(irg,itg,ipg,1,4)*bvxha &
00106 & - cri_grid(irg,itg,ipg,2,4)*bvyha &
00107 & - cri_grid(irg,itg,ipg,3,4)*bvzha
00108 cdbvzd_grid(irg,itg,ipg,2) = pdbvzd_grid(irg,itg,ipg,2) &
00109 & - cri_grid(irg,itg,ipg,1,5)*bvxha &
00110 & - cri_grid(irg,itg,ipg,2,5)*bvyha &
00111 & - cri_grid(irg,itg,ipg,3,5)*bvzha
00112
00113 cdbvxd_grid(irg,itg,ipg,3) = pdbvxd_grid(irg,itg,ipg,3) &
00114 & - cri_grid(irg,itg,ipg,1,3)*bvxha &
00115 & - cri_grid(irg,itg,ipg,2,3)*bvyha &
00116 & - cri_grid(irg,itg,ipg,3,3)*bvzha
00117 cdbvyd_grid(irg,itg,ipg,3) = pdbvyd_grid(irg,itg,ipg,3) &
00118 & - cri_grid(irg,itg,ipg,1,5)*bvxha &
00119 & - cri_grid(irg,itg,ipg,2,5)*bvyha &
00120 & - cri_grid(irg,itg,ipg,3,5)*bvzha
00121 cdbvzd_grid(irg,itg,ipg,3) = pdbvzd_grid(irg,itg,ipg,3) &
00122 & - cri_grid(irg,itg,ipg,1,6)*bvxha &
00123 & - cri_grid(irg,itg,ipg,2,6)*bvyha &
00124 & - cri_grid(irg,itg,ipg,3,6)*bvzha
00125
00126 dbvxdx = cdbvxd_grid(irg,itg,ipg,1)
00127 dbvydx = cdbvyd_grid(irg,itg,ipg,1)
00128 dbvzdx = cdbvzd_grid(irg,itg,ipg,1)
00129
00130 dbvxdy = cdbvxd_grid(irg,itg,ipg,2)
00131 dbvydy = cdbvyd_grid(irg,itg,ipg,2)
00132 dbvzdy = cdbvzd_grid(irg,itg,ipg,2)
00133
00134 dbvxdz = cdbvxd_grid(irg,itg,ipg,3)
00135 dbvydz = cdbvyd_grid(irg,itg,ipg,3)
00136 dbvzdz = cdbvzd_grid(irg,itg,ipg,3)
00137
00138 gmxxd = hxxd(irg,itg,ipg)
00139 gmxyd = hxyd(irg,itg,ipg)
00140 gmxzd = hxzd(irg,itg,ipg)
00141 gmyyd = hyyd(irg,itg,ipg)
00142 gmyzd = hyzd(irg,itg,ipg)
00143 gmzzd = hzzd(irg,itg,ipg)
00144 gmxxd = gmxxd + 1.0d0
00145 gmyyd = gmyyd + 1.0d0
00146 gmzzd = gmzzd + 1.0d0
00147 gmyxd = gmxyd
00148 gmzxd = gmxzd
00149 gmzyd = gmyzd
00150 gmxxu = hxxu(irg,itg,ipg)
00151 gmxyu = hxyu(irg,itg,ipg)
00152 gmxzu = hxzu(irg,itg,ipg)
00153 gmyyu = hyyu(irg,itg,ipg)
00154 gmyzu = hyzu(irg,itg,ipg)
00155 gmzzu = hzzu(irg,itg,ipg)
00156 gmyyu = gmyyu + 1.0d0
00157 gmxxu = gmxxu + 1.0d0
00158 gmzzu = gmzzu + 1.0d0
00159 gmyxu = gmxyu
00160 gmzxu = gmxzu
00161 gmzyu = gmyzu
00162
00163
00164 cdivbv_grid(irg,itg,ipg) = gmxxu*dbvxdx + gmxyu*dbvydx + gmxzu*dbvzdx &
00165 & + gmyxu*dbvxdy + gmyyu*dbvydy + gmyzu*dbvzdy &
00166 & + gmzxu*dbvxdz + gmzyu*dbvydz + gmzzu*dbvzdz
00167 diver = fa23*cdivbv_grid(irg,itg,ipg)
00168
00169 cdivbv_grid(irg,itg,ipg) = pdbvxdx + pdbvydy + pdbvzdz
00170
00171
00172
00173 oelpxx = 0.0d0
00174 oelpxy = 0.0d0
00175 oelpxz = 0.0d0
00176 oelpyy = 0.0d0
00177 oelpyz = 0.0d0
00178 oelpzz = 0.0d0
00179 if (chgra == 'h'.or.chgra == 'c'.or.chgra == 'C' &
00180 &.or.chgra == 'H'.or.chgra == 'W') then
00181 oelpxx = ainvh*ome*elpxx_grid(irg,itg,ipg)*cutoff
00182 oelpxy = ainvh*ome*elpxy_grid(irg,itg,ipg)*cutoff
00183 oelpxz = ainvh*ome*elpxz_grid(irg,itg,ipg)*cutoff
00184 oelpyy = ainvh*ome*elpyy_grid(irg,itg,ipg)*cutoff
00185 oelpyz = ainvh*ome*elpyz_grid(irg,itg,ipg)*cutoff
00186 oelpzz = ainvh*ome*elpzz_grid(irg,itg,ipg)*cutoff
00187 end if
00188
00189 tfkij_grid(irg,itg,ipg,1,1) = ainvh*(2.0d0*dbvxdx-gmxxd*diver) +oelpxx
00190 tfkij_grid(irg,itg,ipg,2,2) = ainvh*(2.0d0*dbvydy-gmyyd*diver) +oelpyy
00191 tfkij_grid(irg,itg,ipg,3,3) = ainvh*(2.0d0*dbvzdz-gmzzd*diver) +oelpzz
00192 tfkij_grid(irg,itg,ipg,1,2) = ainvh*(dbvydx+dbvxdy-gmxyd*diver)+oelpxy
00193 tfkij_grid(irg,itg,ipg,1,3) = ainvh*(dbvzdx+dbvxdz-gmxzd*diver)+oelpxz
00194 tfkij_grid(irg,itg,ipg,2,3) = ainvh*(dbvzdy+dbvydz-gmyzd*diver)+oelpyz
00195 tfkij_grid(irg,itg,ipg,2,1) = tfkij_grid(irg,itg,ipg,1,2)
00196 tfkij_grid(irg,itg,ipg,3,1) = tfkij_grid(irg,itg,ipg,1,3)
00197 tfkij_grid(irg,itg,ipg,3,2) = tfkij_grid(irg,itg,ipg,2,3)
00198
00199 gamu(1,1) = gmxxu
00200 gamu(1,2) = gmxyu
00201 gamu(1,3) = gmxzu
00202 gamu(2,1) = gmyxu
00203 gamu(2,2) = gmyyu
00204 gamu(2,3) = gmyzu
00205 gamu(3,1) = gmzxu
00206 gamu(3,2) = gmzyu
00207 gamu(3,3) = gmzzu
00208
00209 tfkijkij_grid(irg,itg,ipg) = 0.0d0
00210 trk_grid(irg,itg,ipg) = 0.0d0
00211 do id = 1, 3
00212 do ic = 1, 3
00213 do ib = 1, 3
00214 do ia = 1, 3
00215 tfkijkij_grid(irg,itg,ipg) = tfkijkij_grid(irg,itg,ipg) &
00216 & + gamu(ia,ic)*gamu(ib,id) &
00217 & * tfkij_grid(irg,itg,ipg,ia,ib) &
00218 & * tfkij_grid(irg,itg,ipg,ic,id)
00219 end do
00220 end do
00221 end do
00222 end do
00223
00224 if (tfkijkij_grid(irg,itg,ipg) /= 0.) info = 1
00225
00226
00227
00228 rlbxx_grid(irg,itg,ipg) = 2.0d0*dbvxdx
00229 rlbxy_grid(irg,itg,ipg) = dbvydx + dbvxdy
00230 rlbxz_grid(irg,itg,ipg) = dbvzdx + dbvxdz
00231 rlbyy_grid(irg,itg,ipg) = 2.0d0*dbvydy
00232 rlbyz_grid(irg,itg,ipg) = dbvzdy + dbvydz
00233 rlbzz_grid(irg,itg,ipg) = 2.0d0*dbvzdz
00234
00235 end do
00236 end do
00237 end do
00238 if (info /= 1) write(6,*) ' ### Warning K_ij = 0 *** '
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253 end subroutine excurve_WL_gridpoint
