00001 subroutine excurve_WL(chgra)
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, tfkij, tfkijkij, trk, &
00006 & bvxd, bvyd, bvzd, bvxu, bvyu, bvzu
00007 use coordinate_grav_r, only : rg
00008
00009 use def_cristoffel, only : cri
00010 use def_Lie_derivatives, only : elpxx, elpxy, elpxz, elpyy, elpyz, elpzz, &
00011 & rlbxx, rlbxy, rlbxz, rlbyy, rlbyz, rlbzz
00012 use def_shift_derivatives, only : cdbvxd, cdbvyd, cdbvzd, cdivbv, &
00013 & pdbvxd, pdbvyd, pdbvzd
00014 use def_metric_hij, only : hxxd, hxyd, hxzd, hyyd, hyzd, hzzd, &
00015 & hxxu, hxyu, hxzu, hyyu, hyzu, hzzu
00016 use def_cutsw, only : cutfac
00017
00018 use interface_interpo_linear_type0
00019 use interface_grgrad_midpoint
00020 use make_array_3d
00021 use make_array_4d
00022 implicit none
00023
00024 real(8), pointer :: dbvxu(:,:,:,:), dbvyu(:,:,:,:), dbvzu(:,:,:,:)
00025 real(8), pointer :: dfdx(:,:,:), dfdy(:,:,:), dfdz(:,:,:)
00026 real(8) :: gamu(3,3)
00027 real(8) :: ainvh, bvxha, bvyha, bvzha, cutoff, diver, fa23,
00028 dbvxdx, dbvxdy, dbvxdz, &
00029 dbvydx, dbvydy, dbvydz, &
00030 dbvzdx, dbvzdy, dbvzdz, &
00031 gmxxd, gmxyd, gmxzd, gmxxu, gmxyu, gmxzu, &
00032 gmyxd, gmyyd, gmyzd, gmyxu, gmyyu, gmyzu, &
00033 gmzxd, gmzyd, gmzzd, gmzxu, gmzyu, gmzzu, &
00034 oelpxx, oelpxy, oelpxz, &
00035 oelpyy, oelpyz, oelpzz, &
00036 pdbvxdx, pdbvxdy, pdbvxdz, &
00037 pdbvydx, pdbvydy, pdbvydz, &
00038 pdbvzdx, pdbvzdy, pdbvzdz
00039 integer :: ia, ib, ic, id, info, ipg, irg, itg
00040
00041 character(len=1), intent(in) :: chgra
00042
00043 call alloc_array3d(dfdx,1,nrg,1,ntg,1,npg)
00044 call alloc_array3d(dfdy,1,nrg,1,ntg,1,npg)
00045 call alloc_array3d(dfdz,1,nrg,1,ntg,1,npg)
00046 call alloc_array4d(dbvxu,1,nrg,1,ntg,1,npg,1,3)
00047 call alloc_array4d(dbvyu,1,nrg,1,ntg,1,npg,1,3)
00048 call alloc_array4d(dbvzu,1,nrg,1,ntg,1,npg,1,3)
00049
00050
00051
00052
00053
00054 info = 0
00055
00056
00057 fa23 = 2.0d0/3.0d0
00058
00059 call grgrad_midpoint(bvxd,dfdx,dfdy,dfdz)
00060 pdbvxd(1:nrg,1:ntg,1:npg,1) = dfdx(1:nrg,1:ntg,1:npg)
00061 pdbvxd(1:nrg,1:ntg,1:npg,2) = dfdy(1:nrg,1:ntg,1:npg)
00062 pdbvxd(1:nrg,1:ntg,1:npg,3) = dfdz(1:nrg,1:ntg,1:npg)
00063
00064 call grgrad_midpoint(bvyd,dfdx,dfdy,dfdz)
00065 pdbvyd(1:nrg,1:ntg,1:npg,1) = dfdx(1:nrg,1:ntg,1:npg)
00066 pdbvyd(1:nrg,1:ntg,1:npg,2) = dfdy(1:nrg,1:ntg,1:npg)
00067 pdbvyd(1:nrg,1:ntg,1:npg,3) = dfdz(1:nrg,1:ntg,1:npg)
00068
00069 call grgrad_midpoint(bvzd,dfdx,dfdy,dfdz)
00070 pdbvzd(1:nrg,1:ntg,1:npg,1) = dfdx(1:nrg,1:ntg,1:npg)
00071 pdbvzd(1:nrg,1:ntg,1:npg,2) = dfdy(1:nrg,1:ntg,1:npg)
00072 pdbvzd(1:nrg,1:ntg,1:npg,3) = dfdz(1:nrg,1:ntg,1:npg)
00073
00074 call grgrad_midpoint(bvxu,dfdx,dfdy,dfdz)
00075 dbvxu(1:nrg,1:ntg,1:npg,1) = dfdx(1:nrg,1:ntg,1:npg)
00076 dbvxu(1:nrg,1:ntg,1:npg,2) = dfdy(1:nrg,1:ntg,1:npg)
00077 dbvxu(1:nrg,1:ntg,1:npg,3) = dfdz(1:nrg,1:ntg,1:npg)
00078
00079 call grgrad_midpoint(bvyu,dfdx,dfdy,dfdz)
00080 dbvyu(1:nrg,1:ntg,1:npg,1) = dfdx(1:nrg,1:ntg,1:npg)
00081 dbvyu(1:nrg,1:ntg,1:npg,2) = dfdy(1:nrg,1:ntg,1:npg)
00082 dbvyu(1:nrg,1:ntg,1:npg,3) = dfdz(1:nrg,1:ntg,1:npg)
00083
00084 call grgrad_midpoint(bvzu,dfdx,dfdy,dfdz)
00085 dbvzu(1:nrg,1:ntg,1:npg,1) = dfdx(1:nrg,1:ntg,1:npg)
00086 dbvzu(1:nrg,1:ntg,1:npg,2) = dfdy(1:nrg,1:ntg,1:npg)
00087 dbvzu(1:nrg,1:ntg,1:npg,3) = dfdz(1:nrg,1:ntg,1:npg)
00088
00089 do ipg = 1, npg
00090 do itg = 1, ntg
00091 do irg = 1, nrg
00092
00093 cutoff = 1.0d0
00094 if (chgra == 'c'.or.chgra == 'C'.or.chgra == 'W') then
00095 if (rg(irg) > cutfac*pi/ome) cutoff = 0.0d0
00096 end if
00097
00098 call interpo_linear_type0(bvxha,bvxd,irg,itg,ipg)
00099 call interpo_linear_type0(bvyha,bvyd,irg,itg,ipg)
00100 call interpo_linear_type0(bvzha,bvzd,irg,itg,ipg)
00101 call interpo_linear_type0(ainvh,alph,irg,itg,ipg)
00102 ainvh = 0.5d0/ainvh
00103
00104 pdbvxdx = dbvxu(irg,itg,ipg,1)
00105 pdbvxdy = dbvxu(irg,itg,ipg,2)
00106 pdbvxdz = dbvxu(irg,itg,ipg,3)
00107 pdbvydx = dbvyu(irg,itg,ipg,1)
00108 pdbvydy = dbvyu(irg,itg,ipg,2)
00109 pdbvydz = dbvyu(irg,itg,ipg,3)
00110 pdbvzdx = dbvzu(irg,itg,ipg,1)
00111 pdbvzdy = dbvzu(irg,itg,ipg,2)
00112 pdbvzdz = dbvzu(irg,itg,ipg,3)
00113
00114 cdbvxd(irg,itg,ipg,1) = pdbvxd(irg,itg,ipg,1) &
00115 & - cri(irg,itg,ipg,1,1)*bvxha &
00116 & - cri(irg,itg,ipg,2,1)*bvyha &
00117 & - cri(irg,itg,ipg,3,1)*bvzha
00118 cdbvyd(irg,itg,ipg,1) = pdbvyd(irg,itg,ipg,1) &
00119 & - cri(irg,itg,ipg,1,2)*bvxha &
00120 & - cri(irg,itg,ipg,2,2)*bvyha &
00121 & - cri(irg,itg,ipg,3,2)*bvzha
00122 cdbvzd(irg,itg,ipg,1) = pdbvzd(irg,itg,ipg,1) &
00123 & - cri(irg,itg,ipg,1,3)*bvxha &
00124 & - cri(irg,itg,ipg,2,3)*bvyha &
00125 & - cri(irg,itg,ipg,3,3)*bvzha
00126
00127 cdbvxd(irg,itg,ipg,2) = pdbvxd(irg,itg,ipg,2) &
00128 & - cri(irg,itg,ipg,1,2)*bvxha &
00129 & - cri(irg,itg,ipg,2,2)*bvyha &
00130 & - cri(irg,itg,ipg,3,2)*bvzha
00131 cdbvyd(irg,itg,ipg,2) = pdbvyd(irg,itg,ipg,2) &
00132 & - cri(irg,itg,ipg,1,4)*bvxha &
00133 & - cri(irg,itg,ipg,2,4)*bvyha &
00134 & - cri(irg,itg,ipg,3,4)*bvzha
00135 cdbvzd(irg,itg,ipg,2) = pdbvzd(irg,itg,ipg,2) &
00136 & - cri(irg,itg,ipg,1,5)*bvxha &
00137 & - cri(irg,itg,ipg,2,5)*bvyha &
00138 & - cri(irg,itg,ipg,3,5)*bvzha
00139
00140 cdbvxd(irg,itg,ipg,3) = pdbvxd(irg,itg,ipg,3) &
00141 & - cri(irg,itg,ipg,1,3)*bvxha &
00142 & - cri(irg,itg,ipg,2,3)*bvyha &
00143 & - cri(irg,itg,ipg,3,3)*bvzha
00144 cdbvyd(irg,itg,ipg,3) = pdbvyd(irg,itg,ipg,3) &
00145 & - cri(irg,itg,ipg,1,5)*bvxha &
00146 & - cri(irg,itg,ipg,2,5)*bvyha &
00147 & - cri(irg,itg,ipg,3,5)*bvzha
00148 cdbvzd(irg,itg,ipg,3) = pdbvzd(irg,itg,ipg,3) &
00149 & - cri(irg,itg,ipg,1,6)*bvxha &
00150 & - cri(irg,itg,ipg,2,6)*bvyha &
00151 & - cri(irg,itg,ipg,3,6)*bvzha
00152
00153 dbvxdx = cdbvxd(irg,itg,ipg,1)
00154 dbvydx = cdbvyd(irg,itg,ipg,1)
00155 dbvzdx = cdbvzd(irg,itg,ipg,1)
00156
00157 dbvxdy = cdbvxd(irg,itg,ipg,2)
00158 dbvydy = cdbvyd(irg,itg,ipg,2)
00159 dbvzdy = cdbvzd(irg,itg,ipg,2)
00160
00161 dbvxdz = cdbvxd(irg,itg,ipg,3)
00162 dbvydz = cdbvyd(irg,itg,ipg,3)
00163 dbvzdz = cdbvzd(irg,itg,ipg,3)
00164
00165 call interpo_linear_type0(gmxxd,hxxd,irg,itg,ipg)
00166 call interpo_linear_type0(gmxyd,hxyd,irg,itg,ipg)
00167 call interpo_linear_type0(gmxzd,hxzd,irg,itg,ipg)
00168 call interpo_linear_type0(gmyyd,hyyd,irg,itg,ipg)
00169 call interpo_linear_type0(gmyzd,hyzd,irg,itg,ipg)
00170 call interpo_linear_type0(gmzzd,hzzd,irg,itg,ipg)
00171 gmxxd = gmxxd + 1.0d0
00172 gmyyd = gmyyd + 1.0d0
00173 gmzzd = gmzzd + 1.0d0
00174 gmyxd = gmxyd
00175 gmzxd = gmxzd
00176 gmzyd = gmyzd
00177 call interpo_linear_type0(gmxxu,hxxu,irg,itg,ipg)
00178 call interpo_linear_type0(gmxyu,hxyu,irg,itg,ipg)
00179 call interpo_linear_type0(gmxzu,hxzu,irg,itg,ipg)
00180 call interpo_linear_type0(gmyyu,hyyu,irg,itg,ipg)
00181 call interpo_linear_type0(gmyzu,hyzu,irg,itg,ipg)
00182 call interpo_linear_type0(gmzzu,hzzu,irg,itg,ipg)
00183 gmyyu = gmyyu + 1.0d0
00184 gmxxu = gmxxu + 1.0d0
00185 gmzzu = gmzzu + 1.0d0
00186 gmyxu = gmxyu
00187 gmzxu = gmxzu
00188 gmzyu = gmyzu
00189
00190
00191
00192
00193 cdivbv(irg,itg,ipg) = pdbvxdx + pdbvydy + pdbvzdz
00194 diver = fa23*cdivbv(irg,itg,ipg)
00195
00196
00197
00198 oelpxx = 0.0d0
00199 oelpxy = 0.0d0
00200 oelpxz = 0.0d0
00201 oelpyy = 0.0d0
00202 oelpyz = 0.0d0
00203 oelpzz = 0.0d0
00204 if (chgra == 'h'.or.chgra == 'c'.or.chgra == 'C' &
00205 &.or.chgra == 'H'.or.chgra == 'W') then
00206 oelpxx = ainvh*ome*elpxx(irg,itg,ipg)*cutoff
00207 oelpxy = ainvh*ome*elpxy(irg,itg,ipg)*cutoff
00208 oelpxz = ainvh*ome*elpxz(irg,itg,ipg)*cutoff
00209 oelpyy = ainvh*ome*elpyy(irg,itg,ipg)*cutoff
00210 oelpyz = ainvh*ome*elpyz(irg,itg,ipg)*cutoff
00211 oelpzz = ainvh*ome*elpzz(irg,itg,ipg)*cutoff
00212 end if
00213
00214 tfkij(irg,itg,ipg,1,1) = ainvh*(2.0d0*dbvxdx-gmxxd*diver) +oelpxx
00215 tfkij(irg,itg,ipg,2,2) = ainvh*(2.0d0*dbvydy-gmyyd*diver) +oelpyy
00216 tfkij(irg,itg,ipg,3,3) = ainvh*(2.0d0*dbvzdz-gmzzd*diver) +oelpzz
00217 tfkij(irg,itg,ipg,1,2) = ainvh*(dbvydx+dbvxdy-gmxyd*diver)+oelpxy
00218 tfkij(irg,itg,ipg,1,3) = ainvh*(dbvzdx+dbvxdz-gmxzd*diver)+oelpxz
00219 tfkij(irg,itg,ipg,2,3) = ainvh*(dbvzdy+dbvydz-gmyzd*diver)+oelpyz
00220 tfkij(irg,itg,ipg,2,1) = tfkij(irg,itg,ipg,1,2)
00221 tfkij(irg,itg,ipg,3,1) = tfkij(irg,itg,ipg,1,3)
00222 tfkij(irg,itg,ipg,3,2) = tfkij(irg,itg,ipg,2,3)
00223
00224 gamu(1,1) = gmxxu
00225 gamu(1,2) = gmxyu
00226 gamu(1,3) = gmxzu
00227 gamu(2,1) = gmyxu
00228 gamu(2,2) = gmyyu
00229 gamu(2,3) = gmyzu
00230 gamu(3,1) = gmzxu
00231 gamu(3,2) = gmzyu
00232 gamu(3,3) = gmzzu
00233
00234 tfkijkij(irg,itg,ipg) = 0.0d0
00235 trk(irg,itg,ipg) = 0.0d0
00236 do id = 1, 3
00237 do ic = 1, 3
00238 do ib = 1, 3
00239 do ia = 1, 3
00240 tfkijkij(irg,itg,ipg) = tfkijkij(irg,itg,ipg) &
00241 & + gamu(ia,ic)*gamu(ib,id) &
00242 & *tfkij(irg,itg,ipg,ia,ib)*tfkij(irg,itg,ipg,ic,id)
00243 end do
00244 end do
00245 end do
00246 end do
00247
00248 if (tfkijkij(irg,itg,ipg) /= 0.) info = 1
00249
00250
00251
00252 rlbxx(irg,itg,ipg) = 2.0d0*dbvxdx
00253 rlbxy(irg,itg,ipg) = dbvydx + dbvxdy
00254 rlbxz(irg,itg,ipg) = dbvzdx + dbvxdz
00255 rlbyy(irg,itg,ipg) = 2.0d0*dbvydy
00256 rlbyz(irg,itg,ipg) = dbvzdy + dbvydz
00257 rlbzz(irg,itg,ipg) = 2.0d0*dbvzdz
00258
00259 end do
00260 end do
00261 end do
00262 if (info /= 1) write(6,*) ' ### Warning K_ij = 0 *** '
00263
00264 deallocate(dbvxu)
00265 deallocate(dbvyu)
00266 deallocate(dbvzu)
00267 deallocate(dfdx)
00268 deallocate(dfdy)
00269 deallocate(dfdz)
00270
00271 end subroutine excurve_WL
