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