anchor.diff:  
       end function
   
  +    !>------------------------------------------------------------
  +    !! Determine if a point intersected a vertex
  +    !!
  +    !! y is the real point location
  +    !! y0 defines one end of a line
  +    !! y1 defines the other end of the line
  +    !!
  +    !! precision optionally defines the tolerance permitted to consider around the line
  +    !!
  +    !! sets top_vertex, bottom_vertex to True if the point intersected the top or bottom ends of the line (which might be the same)
  +    !! Flips the inside/outside boolean if it intersected one or both verticies
  +    !!
  +    !!------------------------------------------------------------
       subroutine check_vertex_hits(y, y0, y1, inside, top_vertex, bottom_vertex, precision)
           implicit none
           real,    intent(in)    :: y, y0, y1
  
positive.code:             enddo
          enddo
          ! then loop over all y elements on the x ends
          do k=1,ny
              do j=1,nz
                  do i=1,nx,nx-1
                      fieldout(i,j,k)=0
                      local_center = 0
                      do l=1,4
  
negative.code:     !!  loops over y,z,x but geolut is only defined over x,y (for now)
      !!
      !!------------------------------------------------------------
      subroutine geo_interp(fieldout, fieldin, geolut, boundary_only)
          implicit none
          real,                   intent(inout) :: fieldout(:,:,:)
          real,                   intent(in)    :: fieldin(:,:,:)
          type(geo_look_up_table),intent(in)    :: geolut
          logical,                intent(in),   optional :: boundary_only
  
          integer :: ims,ime, jms,jme, kms,kme
          integer:: i,j,k,l, localx,localy, nz_input
          logical :: boundaries
          real :: localw
          real :: local_center
  
          ! This can probably all be done on the 1-n assumption, but this keeps the code more flexible in the future
          ims = lbound(fieldout,1)
          ime = ubound(fieldout,1)
          jms = lbound(fieldout,2)
          jme = ubound(fieldout,2)
          kms = lbound(fieldout,3)
          kme = ubound(fieldout,3)
  
          nz_input = size(fieldin,2)
  
          boundaries = .False.
          if (present(boundary_only)) boundaries = boundary_only
  
          ! if we are only processing the boundary, then make x and y increments be the size of the array
          ! so we only hit the edges of the array
          if (boundaries) then
              call boundary_interpolate(fieldout, fieldin, geolut)
          else
          ! use the geographic lookup table generated earlier to
          ! compute a bilinear interpolation from lo to hi
          ! if we are doing the interior too, just iterate over all x and y
              do k=kms,kme
                  do j=jms,jme
                      do i=ims,ime
                          fieldout(i,j,k)=0
                          local_center = 0
                          ! This loop is used if doing a triangular bi-linear interpolation on the grid
                          ! Need to see all 4 points to compute the "local_center" value for the triangulation
                          do l=1,4
                              localx=geolut%x(l,i,k)
                              localy=geolut%y(l,i,k)
                              local_center = local_center + fieldin(localx, min(j,nz_input), localy)
                              ! This is the actual interpolation adding the value from the two raw grid points
                              if (l < 3) then
  
                                  localw=geolut%w(l,i,k)
                                  fieldout(i,j,k) = fieldout(i,j,k) + fieldin(localx, min(j,nz_input), localy) * localw
                              endif
                          enddo
                          ! This is the actual interpolation adding the value from the center average of all four grid points
                          localw = geolut%w(3,i,k)
                          fieldout(i,j,k) = fieldout(i,j,k) + local_center/4 * localw
  
                          ! This loop is used if doing a straight bi-linear interpolation on the grid
                          ! do l=1,4
                          !     localx=geolut%x(l,i,k)
                          !     localy=geolut%y(l,i,k)
  