00001 subroutine grgrad_4th(fnc,dfdx,dfdy,dfdz,irg,itg,ipg)
00002 use phys_constant, only : long, pi
00003 use grid_parameter, only : nrg, ntg, npg
00004 use coordinate_grav_r, only : rginv
00005 use coordinate_grav_theta, only : dthginv
00006 use coordinate_grav_phi, only : dphiginv
00007 use coordinate_grav_extended
00008 use trigonometry_grav_theta, only : sinthg, costhg, cosecthg
00009 use trigonometry_grav_phi, only : sinphig, cosphig
00010 implicit none
00011 real(long), pointer :: fnc(:,:,:)
00012 real(long), intent(out) :: dfdx, dfdy, dfdz
00013 real(long) :: gr1, gr2, gr3, rv, tv, pv
00014 real(long) :: r5(5), th5(5), phi5(5), fr5(5), ft5(5), fp5(5)
00015 integer :: irg, itg, ipg, irgex, itgex, ipgex
00016 integer :: ir0, it0, ip0, irg0 , itg0 , ipg0, ii
00017 real(long), external :: dfdx_4th
00018
00019
00020
00021
00022
00023
00024 ir0 = min0(irg-2,nrg-4)
00025 it0 = itg-2
00026 ip0 = ipg-2
00027
00028 do ii = 1, 5
00029 irg0 = ir0 + ii - 1
00030 itg0 = it0 + ii - 1
00031 ipg0 = ip0 + ii - 1
00032
00033 r5(ii) = rgex(irg0)
00034 th5(ii) = thgex(itg0)
00035 phi5(ii) = phigex(ipg0)
00036
00037 if (irg.eq.0) then
00038 irgex = irgex_r(irg0)
00039 itgex = ntg/2
00040 ipgex = ipgex_r(0,irg0)
00041 fr5(ii) = fnc(irgex,itgex,ipgex)
00042 irgex = irgex_r(irg0)
00043 itgex = ntg/2
00044 ipgex = ipgex_r(npg/4,irg0)
00045 ft5(ii) = fnc(irgex,itgex,ipgex)
00046 irgex = irgex_r(irg0)
00047 itgex = itgex_r(0,irg0)
00048 ipgex = 0
00049 fp5(ii) = fnc(irgex,itgex,ipgex)
00050 else
00051 if (itg.eq.0.or.itg.eq.ntg) then
00052 irgex = irg
00053 itgex = itgex_th(itg0)
00054 ipgex = ipgex_th(0,itg0)
00055 fr5(ii) = fnc(irgex,itgex,ipgex)
00056 irgex = irg
00057 itgex = itgex_th(itg0)
00058 ipgex = ipgex_th(npg/4,itg0)
00059 ft5(ii) = fnc(irgex,itgex,ipgex)
00060 irgex = irgex_r(irg0)
00061 itgex = itgex_r(itg,irg0)
00062 ipgex = 0
00063 fp5(ii) = fnc(irgex,itgex,ipgex)
00064 else
00065 irgex = irgex_r(irg0)
00066 itgex = itgex_r(itg,irg0)
00067 ipgex = ipgex_r(ipg,irg0)
00068 fr5(ii) = fnc(irgex,itgex,ipgex)
00069 irgex = irg
00070 itgex = itgex_th(itg0)
00071 ipgex = ipgex_th(ipg,itg0)
00072 ft5(ii) = fnc(irgex,itgex,ipgex)
00073 irgex = irg
00074 itgex = itg
00075 ipgex = ipgex_phi(ipg0)
00076 fp5(ii) = fnc(irgex,itgex,ipgex)
00077 end if
00078 end if
00079 end do
00080
00081
00082
00083 rv = rg(irg)
00084 tv = thg(itg)
00085 pv = phig(ipg)
00086
00087 if (irg.eq.0) then
00088 dfdx = dfdx_4th(r5,fr5,rv)
00089 dfdy = dfdx_4th(r5,ft5,rv)
00090 dfdz = dfdx_4th(r5,fp5,rv)
00091 else
00092 if (itg.eq.0) then
00093 dfdx = dfdx_4th(th5,fr5,tv)*rginv(irg)
00094 dfdy = dfdx_4th(th5,ft5,tv)*rginv(irg)
00095 dfdz = dfdx_4th(r5,fp5,rv)
00096 else if (itg.eq.ntg) then
00097 dfdx = - dfdx_4th(th5,fr5,tv)*rginv(irg)
00098 dfdy = - dfdx_4th(th5,ft5,tv)*rginv(irg)
00099 dfdz = - dfdx_4th(r5,fp5,rv)
00100 else
00101 gr1 = dfdx_4th(r5,fr5,rv)
00102 gr2 = dfdx_4th(th5,ft5,tv)*rginv(irg)
00103 gr3 = dfdx_4th(phi5,fp5,pv)*rginv(irg)*cosecthg(itg)
00104 dfdx = gr1*sinthg(itg)*cosphig(ipg) &
00105 & + gr2*costhg(itg)*cosphig(ipg) &
00106 & - gr3*sinphig(ipg)
00107 dfdy = gr1*sinthg(itg)*sinphig(ipg) &
00108 & + gr2*costhg(itg)*sinphig(ipg) &
00109 & + gr3*cosphig(ipg)
00110 dfdz = gr1*costhg(itg) &
00111 & - gr2*sinthg(itg)
00112 end if
00113 end if
00114 end subroutine grgrad_4th
