00001 subroutine excurve_CF_gridpoint_bhex
00002 use phys_constant, only : long
00003 use grid_parameter, only : nrg, ntg, npg
00004 use def_metric, only : alph, bvxd, bvyd, bvzd
00005 use def_metric_excurve_grid, only : tfkij_grid, tfkijkij_grid
00006 use interface_grgrad_4th_gridpoint_bhex
00007 implicit none
00008 integer :: ia, ib, info, irg, itg, ipg
00009 real(long) :: ainvh, diver, cdivbv, fa23,
00010 dbvxdx, dbvxdy, dbvxdz, &
00011 dbvydx, dbvydy, dbvydz, &
00012 dbvzdx, dbvzdy, dbvzdz
00013
00014
00015
00016
00017 fa23 = 2.0d0/3.0d0
00018
00019 do ipg = 0, npg
00020 do itg = 0, ntg
00021 do irg = 0, nrg
00022
00023 call grgrad_4th_gridpoint_bhex(bvxd,dbvxdx,dbvxdy,dbvxdz,irg,itg,ipg)
00024 call grgrad_4th_gridpoint_bhex(bvyd,dbvydx,dbvydy,dbvydz,irg,itg,ipg)
00025 call grgrad_4th_gridpoint_bhex(bvzd,dbvzdx,dbvzdy,dbvzdz,irg,itg,ipg)
00026
00027 ainvh = 0.5d0/alph(irg,itg,ipg)
00028 cdivbv = dbvxdx + dbvydy + dbvzdz
00029 diver = fa23*cdivbv
00030
00031 tfkij_grid(irg,itg,ipg,1,1) = ainvh*(2.0d0*dbvxdx - diver)
00032 tfkij_grid(irg,itg,ipg,2,2) = ainvh*(2.0d0*dbvydy - diver)
00033 tfkij_grid(irg,itg,ipg,3,3) = ainvh*(2.0d0*dbvzdz - diver)
00034 tfkij_grid(irg,itg,ipg,1,2) = ainvh*(dbvydx + dbvxdy)
00035 tfkij_grid(irg,itg,ipg,1,3) = ainvh*(dbvzdx + dbvxdz)
00036 tfkij_grid(irg,itg,ipg,2,3) = ainvh*(dbvzdy + dbvydz)
00037 tfkij_grid(irg,itg,ipg,2,1) = tfkij_grid(irg,itg,ipg,1,2)
00038 tfkij_grid(irg,itg,ipg,3,1) = tfkij_grid(irg,itg,ipg,1,3)
00039 tfkij_grid(irg,itg,ipg,3,2) = tfkij_grid(irg,itg,ipg,2,3)
00040
00041 tfkijkij_grid(irg,itg,ipg) = 0.0d0
00042 do ib = 1, 3
00043 do ia = 1, 3
00044 tfkijkij_grid(irg,itg,ipg) = tfkijkij_grid(irg,itg,ipg) &
00045 & +tfkij_grid(irg,itg,ipg,ia,ib)*tfkij_grid(irg,itg,ipg,ia,ib)
00046 end do
00047 end do
00048
00049 if (tfkijkij_grid(irg,itg,ipg) /= 0.) info = 1
00050
00051 end do
00052 end do
00053 end do
00054 if (info /= 1) write(6,*) ' ### Warning K_ij = 0 *** '
00055
00056 end subroutine excurve_CF_gridpoint_bhex
