00001 subroutine calc_surface_normal_midpoint(rs,rnx,rny,rnz,it,ip)
00002 use trigonometry_grav_theta, only : hsinthg, hcosthg, hcosecthg
00003 use trigonometry_grav_phi, only : hsinphig, hcosphig
00004 use interface_flgrad_midpoint_type0_2Dsurf_sph
00005 use interface_interpo_linear_type0_2Dsurf
00006 implicit none
00007 real(long), pointer :: rs(:,:)
00008 real(long), intent(out) :: rnx, rny, rnz
00009 real(long) :: rnr, rnth, rnphi, dtrs, dprs, hrs, hrsinv
00010 real(long) :: erx, ery, erz, ethx, ethy, ethz, ephix, ephiy, ephiz
00011 integer :: ir, it, ip
00012
00013 call flgrad_midpoint_type0_2Dsurf_sph(rs,dtrs,dprs,it,ip)
00014 call interpo_linear_type0_2Dsurf(hrs,rs,it,ip)
00015 hrsinv = 1.0d0/hrs
00016
00017 rnr = 1.0d0
00018 rnth = - hrsinv*dtrs
00019 rnphi = - hrsinv*hcosecthg(it)*dprs
00020 erx = hsinthg(it)*hcosphig(ip)
00021 ery = hsinthg(it)*hsinphig(ip)
00022 erz = hcosthg(it)
00023 ethx = hcosthg(it)*hcosphig(ip)
00024 ethy = hcosthg(it)*hsinphig(ip)
00025 ethz = - hsinthg(it)
00026 ephix = - hsinphig(ip)
00027 ephiy = hcosphig(ip)
00028 ephiz = 0.0d0
00029 rnx = rnr*erx + rnth*ethx + rnphi*ephix
00030 rny = rnr*ery + rnth*ethy + rnphi*ephiy
00031 rnz = rnr*erz + rnth*ethz + rnphi*ephiz
00032
00033 end subroutine calc_surface_normal_midpoint
