00001 subroutine grgrad_midpoint_r2nd(fnc,dfdx,dfdy,dfdz)
00002 use phys_constant, only : long
00003 use grid_parameter, only : nrg, ntg, npg
00004 use coordinate_grav_r, only : hrginv, drginv
00005 use coordinate_grav_theta, only : dthginv
00006 use coordinate_grav_phi, only : dphiginv
00007 use trigonometry_grav_theta, only : hsinthg, hcosthg, hcosecthg
00008 use trigonometry_grav_phi, only : hsinphig, hcosphig
00009 implicit none
00010 real(long), pointer :: fnc(:,:,:)
00011 real(long), pointer :: dfdx(:,:,:), dfdy(:,:,:), dfdz(:,:,:)
00012 real(long) :: gr1, gr2, gr3, dfncdr, dfncdth, dfncdphi
00013 integer :: irg, itg, ipg
00014
00015
00016
00017
00018
00019
00020 do irg = 1, nrg
00021 do itg = 1, ntg
00022 do ipg = 1, npg
00023 dfncdr = 0.25d0 &
00024 & *(fnc(irg,itg ,ipg ) - fnc(irg-1,itg ,ipg ) &
00025 & + fnc(irg,itg-1,ipg ) - fnc(irg-1,itg-1,ipg ) &
00026 & + fnc(irg,itg ,ipg-1) - fnc(irg-1,itg ,ipg-1) &
00027 & + fnc(irg,itg-1,ipg-1) - fnc(irg-1,itg-1,ipg-1))*drginv(irg)
00028 dfncdth = 0.25d0 &
00029 & *(fnc(irg ,itg,ipg ) - fnc(irg ,itg-1,ipg ) &
00030 & + fnc(irg-1,itg,ipg ) - fnc(irg-1,itg-1,ipg ) &
00031 & + fnc(irg ,itg,ipg-1) - fnc(irg ,itg-1,ipg-1) &
00032 & + fnc(irg-1,itg,ipg-1) - fnc(irg-1,itg-1,ipg-1))*dthginv
00033 dfncdphi = 0.25d0 &
00034 & *(fnc(irg ,itg ,ipg) - fnc(irg ,itg ,ipg-1) &
00035 & + fnc(irg-1,itg ,ipg) - fnc(irg-1,itg ,ipg-1) &
00036 & + fnc(irg ,itg-1,ipg) - fnc(irg ,itg-1,ipg-1) &
00037 & + fnc(irg-1,itg-1,ipg) - fnc(irg-1,itg-1,ipg-1))*dphiginv
00038
00039
00040
00041 gr1 = dfncdr
00042 gr2 = dfncdth*hrginv(irg)
00043 gr3 = dfncdphi*hrginv(irg)*hcosecthg(itg)
00044 dfdx(irg,itg,ipg) = gr1 * hsinthg(itg) * hcosphig(ipg) &
00045 & + gr2 * hcosthg(itg) * hcosphig(ipg) &
00046 & - gr3 * hsinphig(ipg)
00047 dfdy(irg,itg,ipg) = gr1 * hsinthg(itg) * hsinphig(ipg) &
00048 & + gr2 * hcosthg(itg) * hsinphig(ipg) &
00049 & + gr3 * hcosphig(ipg)
00050 dfdz(irg,itg,ipg) = gr1 * hcosthg(itg) &
00051 & - gr2 * hsinthg(itg)
00052 end do
00053 end do
00054 end do
00055 end subroutine grgrad_midpoint_r2nd