subroutine makebox (width, px,py,pz, intx,inty,intz, ier)
 integer ier
 real    width, px(4), py(4), pz(4), intx(1), inty(1), intz(1)

!***********************************************************
!   This subroutine creates a box to the left of a two point line on the
! surface of a sphere. This box polygon is a spherical irregular trapezoid
! or in the degenerate case a spherical triangle with great circles for
! sides and mapped in a counter clockwise order so it can be used directly
! with INSIDE to screen for points inside the box. The two point line forms
! the first side, the two adjacent sides are formed by the great circle arcs
! that bisect the angle between the line and it's adjacent lines in the
! polygon. The opposite side is then formed by constructing the parallel 
! (in a great circle sense) line which is a fixed distance from the initial
! line, but no further than the intersection of the two sides.
! The distance specified for the width of the box is an arc that is some
! fraction of the unit sphere radius and is measured from the midpoint
! of the initial line and along the perpendicular.

! On input:
!       WIDTH = width of the desired box in fraction of unit sphere radius
!       PX,PY,PZ = Arrays of length 4 containing the Cartesian
!           coordinates of the two points making the first edge
!           of the box in indexes p?(2) and p?(3) with the preceeding and
!           following points on the polygon edge in indexes p?(1) and p?(4)
!           respectively.

! p?(1) and p?(4) are altered by this subroutine.

! On output:

!       PX,PY,PZ = Arrays of length 4 containing the Cartesian
!           coordinates of the four points making the complete box
!           in indexes p?(1) and p?(4) having been modified.

!       INTX,INTY,INTZ = Cartesian coordinates of the point where the box
!            ends intersect in the same hemisphere as the original line

!       IER = Error indicator:
!             IER = 0 if no errors were encountered.
!             IER = 1 if p?(1) and p?(2) are the same point
!             IER = 2 if width is zero
!***********************************************************

integer idx
real    vp1m(3), vp1(3), vp2(3), vp2p(3), vm(3), vn(3), vf(3), vg(3)
real    vgm(3), vgp(3), vgint(3), vgn(3), vgmb(3), vgpb(3), vf1(3), vf2(3), vp1n(3), vp2n(3)
real    dint, dend, dintend

! Local parameters:
! vp1m - point preceding polygon segment initial point
! vp2p - point after polygon segment ending point

! width needs to be non zero for below methods to work
! negative width is ok
if( width .eq. 0.0 ) then
    ! zero width, just return zero width box
    px(1) = px(2)
    py(1) = py(2)
    pz(1) = pz(2)
    px(4) = px(3)
    py(4) = py(3)
    pz(4) = pz(3)
    ier = 2
    return
end if

! put two intial points into vector form (one side of box)
vp1m(1) = px(1)
vp1m(2) = py(1)
vp1m(3) = pz(1)
vp1(1) = px(2)
vp1(2) = py(2)
vp1(3) = pz(2)
vp2(1) = px(3)
vp2(2) = py(3)
vp2(3) = pz(3)
vp2p(1) = px(4)
vp2p(2) = py(4)
vp2p(3) = pz(4)

! great circle of two points
call cross_product(vp1, vp2, vg)
call unit_vector(vg, vg, ier)
if( ier .gt. 0 ) then
    ! the two points were identical
    return
end if

! midpoint is in direction of sum of the two input vectors
call add_vector(vp1, vp2, vm)
! normalize for next step
call unit_vector(vm, vm, ier)

! midpoint of second box side is fraction of great circle unit vector
! added to midpoint vector.
call smult_vector(tan(width), vg, vn)
call add_vector(vm, vn, vn)

! great circle through two midpoints
call cross_product(vm, vn, vf)

! great circle for side of box through second side midpoint
call cross_product(vf, vn, vgn)

! find ends of box
! previous side, in  counter clock wise direction is great circle bisecting
! angle between lower indexed polygon segment and this segment.
! great circle of previous segment
call cross_product(vp1m, vp1, vgm)
call unit_vector(vgm, vgm, ier)
if( ier .gt. 0 ) then
    ! the two points were identical
    vgm = vg
    ier = 0
end if

! great circle vector for bisection of angle between the two segments
! add two segment normal unit vectors to get direction
call add_vector(vg, vgm, vgmb)
! cross with first point to get great circle vector
call cross_product(vp1, vgmb, vf1)
call unit_vector(vf1, vf1, ier)

! next side, in  counter clock wise direction is great circle bisecting
! angle between higher indexed polygon segment and this segment.
! great circle of previous segment
call cross_product(vp2, vp2p, vgp)
call unit_vector(vgp, vgp, ier)
! great circle vector for bisection of angle between the two segments
! add two segment normal unit vectors to get direction
call add_vector(vg, vgp, vgpb)
! cross with first point to get great circle vector
call cross_product(vp2, vgpb, vf2)
call unit_vector(vf2, vf2, ier)

! find point of intersection between two box sides
call cross_product(vf1, vf2, vgint)
! check for correct hemishpere
if( dot_product(vgint,vp1) .lt. 0 ) then
    call cross_product(vf2, vf1, vgint)
end if
call unit_vector(vgint, vgint, ier)

! assign to return value array.
intx(1) = vgint(1)
inty(1) = vgint(2)
intz(1) = vgint(3)

! the two other corners of the box are the intersections of the 
! second side great circle with end of box great circles
call cross_product(vgn, vf1, vp1n)
call unit_vector(vp1n, vp1n, ier)
call cross_product(vgn, vf2, vp2n)
call unit_vector(vp2n, vp2n, ier)

! find distance from segment end point to point of intersection of box ends
dint = acos( dot_product(vp1,vgint))
! find distance from segment end point to box end intersection with fixed distance parallel line
dend = acos( dot_product(vp1,vp1n))
! find distance from fixed distance parallel line to point of intersection of box ends
dintend = acos( dot_product(vgint,vp1n))

! compare distances to determine if box or triangle
if( (dintend .gt. dend) .or. (dint .gt. dend) ) then
    ! intersection is on the outside of the polygon, concave outer boundary
    ! or intersection point is inside but past the end line boundary
    ! place new points in return arrays in counter clock wise order
    px(4) = vp2n(1)
    py(4) = vp2n(2)
    pz(4) = vp2n(3)
    px(1) = vp1n(1)
    py(1) = vp1n(2)
    pz(1) = vp1n(3)
else
    ! intersection is on the inside of the polygon, convex outer boundary
    ! and intersection is closer than end line boundary.
    ! place new points in return arrays in counter clock wise order
    px(4) = vgint(1)
    py(4) = vgint(2)
    pz(4) = vgint(3)
    px(1) = vgint(1)
    py(1) = vgint(2)
    pz(1) = vgint(3)
end if

! No error encountered.

ier = 0
return
end

subroutine rtrans (n, x,y,z,rlat,rlon)
integer n
real    rlat(n), rlon(n), x(n), y(n), z(n)

!   This subroutine transforms Cartesian coordinates on the unit sphere
!  back to spherical coordinates for output results. Storage for RLAT
!  and RLON may coincide with  Storage for X and Y if the latter need
!  not be saved.

! On input:

!       N = Number of nodes (points on the unit sphere)
!           whose coordinates are to be transformed.

!       X,Y,Z = Arrays of length at least N.

! The above parameters are not altered by this routine.

!       RLAT = Array of length N containing latitudinal
!              coordinates of the nodes in radians.

!       RLON = Array of length N containing longitudinal
!              coordinates of the nodes in radians.

! On output:

!       RLAT, RLON = Spherical coordinates
!               X(I)**2 + Y(I)**2 + Z(I)**2 = 1 for I = 1
!               to N.

! Modules required by TRANS:  None

! Intrinsic functions called by TRANS:  COS, SIN

integer idx, nn
real rho, sxy, phi, theta, pi

! Local parameters:

! rho = radius
! sxy = xy distance out from pole
! idx =      DO-loop index
! NN =     Local copy of N
! PHI =    Latitude
! THETA =  Longitude

pi = 2 * asin(1.0)
nn = n
do idx = 1,nn
  rho = sqrt(x(idx)*x(idx)+y(idx)*y(idx)+z(idx)*z(idx))
  sxy = sqrt(x(idx)*x(idx)+y(idx)*y(idx))
  phi = asin( z(idx) / rho )
  if( y(idx) .gt. 0.0 ) then
      theta = acos( x(idx) / sxy )
  else
      theta = - acos( x(idx) / sxy )
  end if
  rlat(idx) = phi
  rlon(idx) = theta
end do
return
end

subroutine cross_product(va, vb, vc)
real :: va(3), vb(3), vc(3)

! calculate the vector cross product of va x vb
! and returns the answer in vc 

! va, vb, vc are 3 point vectors in x, y, z order

vc(1) = va(2)*vb(3) - va(3)*vb(2)
vc(2) = va(3)*vb(1) - va(1)*vb(3)
vc(3) = va(1)*vb(2) - va(2)*vb(1)

return
end

real function norm(va)
real :: va(3)

! calculate the vector norm of va
! and returns the answer as the function value

! va is a 3 point vector in x, y, z order

norm = sqrt( dot_product( va, va ) )

return
end

subroutine unit_vector(vi, vo, ier)
integer ier
real :: vi(3), vo(3)

integer idx
real norm
real magnitude

! takes a vector vi and normalizes it to make it a unit vector
! and returns the answer in vo

! vi, vo are 3 point vectors in x, y, z order

magnitude = norm(vi)
if( magnitude .le. 0.0 ) then
    ! zero length vector on input
    ier = 1
    ! allow return of zero length vector
    magnitude = 1
else
    ier = 0
end if

do idx = 1,3
    vo(idx) = vi(idx) / magnitude
end do

return
end

subroutine add_vector(v1, v2, vsum)
real :: v1(3), v2(3), vsum(3)

integer idx

! takes two vectors, adds them together, and returns the answer in vsum
! v1, v2 and vsum are 3 point vectors in x, y, z order
do idx = 1,3
    vsum(idx) = v1(idx) + v2(idx)
end do

return
end

subroutine smult_vector(mult, vi, vo)
real :: mult, vi(3), vo(3)

integer idx

! takes a scalar value and a vector, multiplies the vector by the scalar value, and returns the answer in vo
! vi and vo are 3 point vectors in x, y, z order
do idx = 1,3
    vo(idx) = mult * vi(idx)
end do

return
end

subroutine mirror_point (pt, px,py,pz, mx,my,mz, ier)
integer ier
real    pt(3), px(4), py(4), pz(4), mx(4), my(4), mz(4)

! This subroutine takes a point, and two four point boxes. The first box
! contains the point which is to be mirrored into the second box. The
! first two points in each box is the line across which the point is to
! be mirrored. They do not need to be the same line, but should have the
! same sense as if they were. If one box is larger in any direction, the
! point will be proportionally placed. What this means is that the first
! box is counterclockwise mapped and the second box is clockwise mapped.

! On input:
!       PT = Cartesian coordinates 1,2,3 are x,y,z of the point to be mirrored
!       PX,PY,PZ = Arrays of length 4 containing the Cartesian
!           coordinates of the four points making up the box containing
!           the point to be mirrored. Indexes indexes p?(1) and p?(2) are
!           the line across which the point to be mirrored will be
!           reflected.
!       MX,MY,MZ = Arrays of length 4 containing the Cartesian
!           coordinates of the four points making up the box which will
!           contain the reflected point. Indexes indexes p?(1) and p?(2) are
!           the line across which the point to be mirrored will be
!           reflected.

! PT is altered by this subroutine.

! On output:

!       PT = Cartesian coordinates 1,2,3 are x,y,z of the point mirrored
!           into the second box

!       IER = Error indicator:
!             IER = 0 if no errors were encountered.
!             IER = 1 if p?(1) and p?(2) are the same point
!             IER = 2 if p?(3) is on the same line as p?(1) and p?(2)
!***********************************************************

integer :: idx
real :: vp1(3), vp2(3), vp3(3), vp4(3), vpi(3)
real :: vm1(3), vm2(3), vm3(3), vm4(3), vmi(3)
real :: vcs1(3), vcs2(3), vis(3), vcip(3), vcb(3), vipb(3), vcbo(3), vipbo(3)
real :: angps1, angs12, ratps1
real :: angpb, angbo, ratpb
real :: vdb(3), vdp(3)

! Local parameters:
! all coordinates are cartesian from sphere center
! vp? - vectors of the points making up the sides of the box containing the point to be mirrored
! vm? - vectors of the points making up the sides of the box to which the point will be mirrored
! vcs1 - great circle vector of the first side of the box (points 1 and 4)
! vcs2 - great circle vector of the second side of the box (points 2 and 3)
! vis - vector indicating the nearest intersection point of the two sides of the box
! vcip - great circle vector of the arc passing through the box side intersection point and the point to be mirrored
! vcb - great circle vector of the boundary arc across which the point is to be mirrored
! vipb - vector of intersection between point arc and boundary arc
! vcbo - great circle vector of the boundary arc opposite to the boundary arc across which the point is to be mirrored
! vipbo - vector of intersection between point arc and the boundary arc opposite to the boundary arc across which the point is to be mirrored
! angps1 - angle from side one to point arc along the boundary arc
! angs12 - angle from side one to side two along the boundary arc
! ratps1 - ratio of side one to point angle to side one to side two angle
! angpb - angle from boundary arc to point along point arc
! angbo - angle from boundary arc to opposite boundary along point arc
! ratpb - ratio of boundary arc to point angle to boundary arc to opposite boundary angle
! vdb - find vector in direction of boundary arc (perpendicular to great circle vector of boundary arc)
! vdp - find vector in direction of point arc (perpendicular to great circle vector of point arc)

! put four inner box points into vector form
vp1(1) = px(1)
vp1(2) = py(1)
vp1(3) = pz(1)
vp2(1) = px(2)
vp2(2) = py(2)
vp2(3) = pz(2)
vp3(1) = px(3)
vp3(2) = py(3)
vp3(3) = pz(3)
vp4(1) = px(4)
vp4(2) = py(4)
vp4(3) = pz(4)

! put four outer box points into vector form
vm1(1) = mx(1)
vm1(2) = my(1)
vm1(3) = mz(1)
vm2(1) = mx(2)
vm2(2) = my(2)
vm2(3) = mz(2)
vm3(1) = mx(3)
vm3(2) = my(3)
vm3(3) = mz(3)
vm4(1) = mx(4)
vm4(2) = my(4)
vm4(3) = mz(4)

! vcs1 - great circle vector of the first side of the box (points 1 and 4)
call cross_product(vp1, vp4, vcs1)
call unit_vector(vcs1, vcs1, ier)
if( ier .gt. 0 ) then
    ! the two points were identical
    return
end if

! vcs2 - great circle vector of the second side of the box (points 2 and 3)
call cross_product(vp2, vp3, vcs2)
call unit_vector(vcs2, vcs2, ier)
if( ier .gt. 0 ) then
    ! the two points were identical
    return
end if

! vis - vector indicating the nearest intersection point of the two sides of the box
call cross_product( vcs1, vcs2, vis )
! check vector direction
if( dot_product(vp3,vis) .lt. 0 ) then
    ! vector was in wrong direction, reverse
    call cross_product( vcs2, vcs1, vis )
end if
! unitize to get point on the surface of the sphere
call unit_vector( vis, vis, ier)

! vcip - great circle vector of the arc passing through the box side intersection point and the point to be mirrored
call cross_product(vis, pt, vcip)
call unit_vector(vcip, vcip, ier)
if( ier .gt. 0 ) then
    ! the two points were identical
    return
end if

! vcb - great circle vector of the boundary arc across which the point is to be mirrored (points 1 and 2)
call cross_product(vp1, vp2, vcb)
call unit_vector(vcb, vcb, ier)
if( ier .gt. 0 ) then
    ! the two points were identical
    return
end if

! vipb - vector of intersection between point arc and boundary arc
call cross_product( vcip, vcb, vipb )
! check vector direction
if( dot_product(vp1,vipb) .lt. 0 ) then
    ! vector was in wrong direction, reverse
    call cross_product( vcb, vcip, vipb )
end if
! unitize to get point on the surface of the sphere
call unit_vector( vipb, vipb, ier)

! vcbo - great circle vector of the boundary arc opposite to the boundary arc across which the point is to be mirrored
call cross_product(vp3, vp4, vcbo)
call unit_vector(vcbo, vcbo, ier)
if( ier .gt. 0 ) then
    ! the two points were identical, making these points also the intersection point of the two sides and the point needed for the distance calculation
    ! vipbo - vector of intersection between point arc and the boundary arc opposite to the boundary arc across which the point is to be mirrored
    vipbo = vis
else
    ! vipbo - vector of intersection between point arc and the boundary arc opposite to the boundary arc across which the point is to be mirrored
    call cross_product( vcip, vcbo, vipbo )
    ! check vector direction
    if( dot_product(vp3,vipbo) .lt. 0 ) then
        ! vector was in wrong direction, reverse
        call cross_product( vcbo, vcip, vipbo )
    end if
    ! unitize to get point on the surface of the sphere
    call unit_vector( vipbo, vipbo, ier)
end if

! angps1 - angle from side one to point arc along the boundary arc
angps1 = acos( dot_product(vp1,vipb) )

! angs12 - angle from side one to side two along the boundary arc
angs12 = acos( dot_product(vp1,vp2) )

! ratps1 - ratio of side one to point angle to side one to side two angle
ratps1 = angps1/angs12

! angpb - angle from boundary arc to point along point arc
angpb = acos( dot_product(vipb,pt) )

! angbo - angle from boundary arc to opposite boundary along point arc
angbo = acos( dot_product(vipb,vipbo) )

! ratpb - ratio of boundary arc to point angle to boundary arc to opposite boundary angle
ratpb = angpb / angbo

!write(*,*) "# ratps1, ratpb", ratps1, ratpb

! find vectors needed for second box to reverse the calculations
! vcs1 - great circle vector of the first side of the box (points 1 and 4)
call cross_product(vm1, vm4, vcs1)
call unit_vector(vcs1, vcs1, ier)
if( ier .gt. 0 ) then
    ! the two points were identical
    return
end if

! vcs2 - great circle vector of the second side of the box (points 2 and 3)
call cross_product(vm2, vm3, vcs2)
call unit_vector(vcs2, vcs2, ier)
if( ier .gt. 0 ) then
    ! the two points were identical
    return
end if

! vis - vector indicating the nearest intersection point of the two sides of the box
call cross_product( vcs1, vcs2, vis )
! check vector direction
if( dot_product(vm3,vis) .lt. 0 ) then
    ! vector was in wrong direction, reverse
    call cross_product( vcs2, vcs1, vis )
end if
! unitize to get point on the surface of the sphere
call unit_vector( vis, vis, ier)

! vcb - great circle vector of the boundary arc across which the point is to be mirrored (points 1 and 2)
call cross_product(vm1, vm2, vcb)
call unit_vector(vcb, vcb, ier)
if( ier .gt. 0 ) then
    ! the two points were identical
    return
end if

! vdb - find vector in direction of boundary arc (perpendicular to great circle vector of boundary arc)
call cross_product( vm1, vcb, vdb )
! check vector direction (direction should be toward vm2)
if( dot_product(vm2,vdb) .lt. 0 ) then
    ! vector was in wrong direction, reverse
    call cross_product( vcb, vm1, vdb )
end if
! unitize to get point on the surface of the sphere
call unit_vector( vdb, vdb, ier)

! angs12 - angle from side one to side two along the boundary arc
angs12 = acos( dot_product(vm1,vm2) )
angps1 = ratps1*angs12

! vipb - vector of intersection between point arc and boundary arc
call smult_vector( tan(angps1),vdb ,vipb )
call add_vector( vipb, vm1, vipb )
! unitize to get point on the surface of the sphere
call unit_vector( vipb, vipb, ier)

! vcip - great circle vector of the arc passing through the box side intersection point and boundary intersection of point arc
call cross_product(vis, vipb, vcip)
call unit_vector(vcip, vcip, ier)
if( ier .gt. 0 ) then
    ! the two points were identical
    return
end if

! vcbo - great circle vector of the boundary arc opposite to the boundary arc across which the point is to be mirrored
call cross_product(vm3, vm4, vcbo)
call unit_vector(vcbo, vcbo, ier)
if( ier .gt. 0 ) then
    ! the two points were identical, making these points also the intersection point of the two sides and the point needed for the distance calculation
    ! vipbo - vector of intersection between point arc and the boundary arc opposite to the boundary arc across which the point is to be mirrored
    vipbo = vis
else
    ! vipbo - vector of intersection between point arc and the boundary arc opposite to the boundary arc across which the point is to be mirrored
    call cross_product( vcip, vcbo, vipbo )
    ! check vector direction
    if( dot_product(vm3,vipbo) .lt. 0 ) then
        ! vector was in wrong direction, reverse
        call cross_product( vcbo, vcip, vipbo )
    end if
    ! unitize to get point on the surface of the sphere
    call unit_vector( vipbo, vipbo, ier)
end if

! vdp - find vector in direction of point arc (perpendicular to great circle vector of point arc)
call cross_product( vipb, vcip, vdp )
! check vector direction (direction should be toward vp2)
if( dot_product(vipbo,vdp) .lt. 0 ) then
    ! vector was in wrong direction, reverse
    call cross_product( vcip, vipb, vdp )
end if
! unitize to get point on the surface of the sphere
call unit_vector( vdp, vdp, ier)

! angbo - angle from boundary arc to opposite boundary along point arc
angbo = acos( dot_product(vipb,vipbo) )
angpb = ratpb * angbo

! find location of point
call smult_vector( tan(angpb),vdp ,pt )
call add_vector( pt, vipb, pt )
! unitize to get point on the surface of the sphere
call unit_vector( pt, pt, ier)

! no errors
ier = 0

return
end

subroutine tobound_point (pt, px,py,pz, ier)
integer ier
real    pt(3), px(4), py(4), pz(4)

! This subroutine takes a point, and a four point box which
! contains the point which is to be moved to the boundary. The
! first two points in the box are the line to which the point is to
! be moved.

! On input:
!       PT = Cartesian coordinates 1,2,3 are x,y,z of the point to be mirrored
!       PX,PY,PZ = Arrays of length 4 containing the Cartesian
!           coordinates of the four points making up the box containing
!           the point to be moved to the boundary. Indexes indexes p?(1) and p?(2) are
!           the line across which the point to be mirrored will be
!           reflected.

! PT is altered by this subroutine.

! On output:

!       PT = Cartesian coordinates 1,2,3 are x,y,z of the point mirrored
!           into the second box

!       IER = Error indicator:
!             IER = 0 if no errors were encountered.
!             IER = 1 if p?(1) and p?(2) are the same point
!***********************************************************

integer :: idx
real :: vp1(3), vp2(3), vp3(3), vp4(3), vpi(3)
real :: vcs1(3), vcs2(3), vis(3), vcip(3), vcb(3), vipb(3)

! Local parameters:
! all coordinates are cartesian from sphere center
! vp? - vectors of the points making up the sides of the box containing the point to be moved
! vcs1 - great circle vector of the first side of the box (points 1 and 4)
! vcs2 - great circle vector of the second side of the box (points 2 and 3)
! vis - vector indicating the nearest intersection point of the two sides of the box
! vcip - great circle vector of the arc passing through the box side intersection point and the point to be moved
! vcb - great circle vector of the boundary arc across which the point is to be moved
! vipb - vector of intersection between point arc and boundary arc

! put four inner box points into vector form
vp1(1) = px(1)
vp1(2) = py(1)
vp1(3) = pz(1)
vp2(1) = px(2)
vp2(2) = py(2)
vp2(3) = pz(2)
vp3(1) = px(3)
vp3(2) = py(3)
vp3(3) = pz(3)
vp4(1) = px(4)
vp4(2) = py(4)
vp4(3) = pz(4)

! vcs1 - great circle vector of the first side of the box (points 1 and 4)
call cross_product(vp1, vp4, vcs1)
call unit_vector(vcs1, vcs1, ier)
if( ier .gt. 0 ) then
    ! the two points were identical
    return
end if

! vcs2 - great circle vector of the second side of the box (points 2 and 3)
call cross_product(vp2, vp3, vcs2)
call unit_vector(vcs2, vcs2, ier)
if( ier .gt. 0 ) then
    ! the two points were identical
    return
end if

! vis - vector indicating the nearest intersection point of the two sides of the box
call cross_product( vcs1, vcs2, vis )
! check vector direction
if( dot_product(vp3,vis) .lt. 0 ) then
    ! vector was in wrong direction, reverse
    call cross_product( vcs2, vcs1, vis )
end if
! unitize to get point on the surface of the sphere
call unit_vector( vis, vis, ier)

! vcip - great circle vector of the arc passing through the box side intersection point and the point to be moved
call cross_product(vis, pt, vcip)
call unit_vector(vcip, vcip, ier)
if( ier .gt. 0 ) then
    ! the two points were identical
    return
end if

! vcb - great circle vector of the boundary arc across which the point is to be moved (points 1 and 2)
call cross_product(vp1, vp2, vcb)
call unit_vector(vcb, vcb, ier)
if( ier .gt. 0 ) then
    ! the two points were identical
    return
end if

! vipb - vector of intersection between point arc and boundary arc
call cross_product( vcip, vcb, vipb )
! check vector direction
if( dot_product(vp1,vipb) .lt. 0 ) then
    ! vector was in wrong direction, reverse
    call cross_product( vcb, vcip, vipb )
end if
! unitize to get point on the surface of the sphere
call unit_vector( vipb, vipb, ier)

! return intersected value which is on boundary
pt = vipb

! no errors
ier = 0

return
end

logical function check_cross (px,py,pz)
real    px(4), py(4), pz(4)

! This subroutine takes four points, the first two of which describe a
! great circle arc A, and points 3 and 4 which describe a great circle
! arc B. The arc is the short distance, not the long distance. The
! function returns TRUE if the two arcs have an intersection point and 
! FALSE if they do not. The intersection point in case of TRUE is not found.

! On input:
!       PX,PY,PZ = Arrays of length 4 containing the Cartesian
!           coordinates of the four points to be tested
!       These are not altered by the routine

! On return:

!       TRUE if the two arcs have an intersection point
!       FALSE if they do not
!***********************************************************

integer idx, ier
real    vp1(3), vp2(3), vp3(3), vp4(3), vf(3), vg(3), vi(3)
logical lft1, lft2

! Local parameters:

! put four intial points in to vector form
vp1(1) = px(1)
vp1(2) = py(1)
vp1(3) = pz(1)
vp2(1) = px(2)
vp2(2) = py(2)
vp2(3) = pz(2)
vp3(1) = px(3)
vp3(2) = py(3)
vp3(3) = pz(3)
vp4(1) = px(4)
vp4(2) = py(4)
vp4(3) = pz(4)

! find the great circle normal vector for the first two points
call cross_product(vp1, vp2, vg)
call unit_vector(vg, vg, ier)
if( ier .gt. 0 ) then
    ! the two points were identical
    check_cross = .FALSE.
    return
end if

! find the great circle normal vector for the last two points
call cross_product(vp3, vp4, vf)
call unit_vector(vf, vf, ier)
if( ier .gt. 0 ) then
    ! the two points were identical
    check_cross = .FALSE.
    return
end if

! check for position of first two points with relation to great circle formed by last two points
lft1 = (dot_product(vp1, vf) .gt. 0.0)
lft2 = (dot_product(vp2, vf) .gt. 0.0)
if( lft1 .eqv. lft2 ) then
    ! similar sign of the dot product means that both vectors are on the same side of the great circle.
    check_cross = .FALSE.
    return
end if

! check for position of last two points with relation to great circle formed by first two points
lft1 = (dot_product(vp3, vg) .gt. 0.0)
lft2 = (dot_product(vp4, vg) .gt. 0.0)
if( lft1 .eqv. lft2 ) then
    ! similar sign of the dot product means that both vectors are on the same side of the great circle.
    check_cross = .FALSE.
    return
end if

! if we reached here, then the two arcs must intersect
check_cross = .TRUE.
call cross_product(vf, vg, vi)
call unit_vector(vi, vi, ier)
if( ier .gt. 0 ) then
    ! the points are really colinear or the same (degenerate case)
    check_cross = .FALSE.
end if

return
end

logical function arc_longer (px,py,pz)
real    px(4), py(4), pz(4)

! This subroutine takes four points, the first two of which describe a
! great circle arc A, and points 3 and 4 which describe a great circle
! arc B. The arc is the short distance, not the long distance. The
! function returns TRUE if arc A is longer and FALSE if arc B is longer.

! On input:
!       PX,PY,PZ = Arrays of length 4 containing the Cartesian
!           coordinates of the four points to be tested
!       These are not altered by the routine

! On return:

!       TRUE if the first arc is longer
!       FALSE if the second arc is longer
!       in the case of a tie, TRUE will be returned
!***********************************************************

real    vp1(3), vp2(3), vp3(3), vp4(3)
real    dist1, dist2

! Local parameters:

! put four intial points in to vector form
vp1(1) = px(1)
vp1(2) = py(1)
vp1(3) = pz(1)
vp2(1) = px(2)
vp2(2) = py(2)
vp2(3) = pz(2)
vp3(1) = px(3)
vp3(2) = py(3)
vp3(3) = pz(3)
vp4(1) = px(4)
vp4(2) = py(4)
vp4(3) = pz(4)

! find the arc angle for the first two points
dist1 = acos( dot_product(vp1,vp2))
dist2 = acos( dot_product(vp3,vp4))

if( dist1 .ge. dist2 ) then
    arc_longer = .TRUE.
else
    arc_longer = .FALSE.
end if

return
end

subroutine outside_weights( pt, nval, x,y,z, nbnd, blist, b1, b2, i1, i2, ier, nd, dx, dy, dz)
integer :: nval, nbnd, i1, i2, ier
integer :: blist(1:nval)
real :: pt(1:3), x(1:nval), y(1:nval), z(1:nval), b1, b2
integer :: nd
real :: dx(1:nd), dy(1:nd), dz(1:nd)

! this subrountien take a point which is outside of the convex hull of a 
! triangulation, finds the nearest boundary point and it's two neighbors.
! The location of these points is used to select the two points to be used
! for averaging, and the weights calculated.

! pt - point outside of triangulation to be located and weighted (cartesian coordinates)
! nval - number of nodes in the triangulation array
! x, y, z - cartesian coordinates of triangulation nodes
! nbnd - number of boundary nodes in the triangulation
! b1, b2, - point weighting of two identified points
! i1 - index of the closest point
! i2 - index of the neighbor which will be used to weight along the boundary
! blist - order list of boundary nodes

! local variables
integer :: idx, ndx
real :: vp(3), vpn(3), vpb(3), vpa(3), vca(3), vcn(3), vcb(3), vdn(3)
integer :: in, ib, ia
real :: dist, distmin
logical :: ldirb, ldirp
real :: vp2(3), vc2(3), vcp(3), vip(3), vd12(3), vd1ip(3)

! idx - array index
! ndx - index in blist of nearest boundary node
! vp - vector form of nodal point
! vpn - vector of nearest boundary node to point
! vpb - vector of boundary node that comes before nearest point
! vpa - vector of boundary node that comes after nearest point
! in - index of nearest boundary node
! ib = index of boundary node that comes before nearest point
! ia = index of boundary node that comes after nearest point
! dist - measure of distance from point to boundary node
! distmin - measure representing the minimum distance
! ldirb - direction sense of the before boundary node with reference to the angle bisector
! ldirp - direction sense of the point with reference to the angel bisector
! vp2 - vecor of second point (selected to be before or after point depending on side)
! vc2 - great circle vector of the two boundary nodes on the same side of the angle bisector as the point
! vcp - great dircle through point and perpendicular to great circle of boundary segment (side already selected)
! vip - point of intersection where perpendicular to point hits boundary segment

! find nearest boundary node
distmin = 0.0 ! this is really the largest distance with doT_product (which give the inverse sense)
do idx = 1, nbnd
    ! place nodal point into vector
    vp(1) = x(blist(idx))
    vp(2) = y(blist(idx))
    vp(3) = z(blist(idx))
    ! find distance measure
    dist = dot_product(pt, vp)
    ! using dot product, so same inverse sense. (larger is closer)
    if( dist .gt. distmin ) then
        in = blist(idx)
        ndx = idx
        distmin = dist
    end if
end do
! place nearest boundary node into vector
vpn(1) = x(in)
vpn(2) = y(in)
vpn(3) = z(in)

! find preceeding boundary node (before)
if( ndx .eq. 1 ) then
    ib = blist(nbnd)
else
    ib = blist(ndx-1)
end if
! place into vector
vpb(1) = x(ib)
vpb(2) = y(ib)
vpb(3) = z(ib)

! find following boundary node (after)
if( ndx .eq. nbnd ) then
    ia = blist(1)
else
    ia = blist(ndx+1)
end if
vpa(1) = x(ia)
vpa(2) = y(ia)
vpa(3) = z(ia)

! find dividing line between using before and after (bisect angle made by the three points)
! great circle of previous segment
call cross_product(vpa, vpn, vca)
call unit_vector(vca, vca, ier)

! great circle of following segment
call cross_product(vpn, vpb, vcb)
call unit_vector(vcb, vcb, ier)

! add two segment normal unit vectors to get direction
call add_vector(vca, vcb, vdn)
call unit_vector(vdn, vdn, ier)

! find great circle vector of angle bisection
call cross_product(vpn, vdn, vcn)
call unit_vector(vcn, vcn, ier)

! set return values
i1 = in
! find whether point is on side of before or after point
ldirb = (dot_product(vcn, vpb) .gt. 0.0)
ldirp = (dot_product(vcn, pt) .gt. 0.0)
if( ldirb .eqv. ldirp ) then
    ! before point is on the same side
    i2 = ib
    vp2 = vpb
else
    ! after point must be on the same side (or point is on dividing line)
    i2 = ia
    vp2 = vpa
end if

! find proportional coefficients

! find great circle vector from nearest point to point 2
call cross_product(vpn, vp2, vc2)
call unit_vector(vc2, vc2, ier)

! find great dircle through point and perpendicular to great circle of boundary segment (side already selected)
call cross_product(vc2, pt, vcp)
call unit_vector(vcp, vcp, ier)

! find point of intersection where perpendicular to point hits boundary segment
call cross_product(vc2, vcp, vip)
! check that intersection is in correct hemisphere
if( dot_product( vip, vpn ) .lt. 0.0 ) then
    ! wrong hemisphere, reverse cross direction
    call cross_product(vcp, vc2, vip)
end if
call unit_vector(vip, vip, ier)

! get direction vector from nearest boundary node to intersection point of perpendicular
call cross_product(vpn, vip, vd1ip)
! get direction vector from nearest boundary node to adjacent boundary node
call cross_product(vpn, vp2, vd12)
! check if intersection point is between two boundary nodes
if( dot_product( vd1ip, vd12 ) .gt. 0.0 ) then
    ! intersection point is between two boundary nodes, find proportion
    b1 = 1.0 - acos(dot_product(vpn,vip)) / acos(dot_product(vpn,vp2))
    b2 = 1.0 - b1
else
    ! intersection point is not between two points, use nearest point only
    b1 = 1.0
    b2 = 0.0
end if

! debug
dx(1) = vip(1)
dy(1) = vip(2)
dz(1) = vip(3)

return
end