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