anchor.diff:  
           w3 = 1 - w1 - w2
   
  -        if (minval([w1,w2,w3]) < -1e-3) then
  +        if (minval([w1,w2,w3]) < -1e-4) then
               write(*,*) "ERROR: Point not located in bounding triangle"
               write(*,*) xi, yi
               write(*,*) "Triangle vertices"
  
positive.code:         ! locals
          real::mindist, curdist,dx,dy,xsign,ysign,x,y    ! temporary variables to use during the search
          integer::nx,ny,xc,yc,xw,yw,iterations,xstep,ystep   ! more temporary variables
  
          nx = size(lo%lat, 1)
          ny = size(lo%lat, 2)
          ! calcualte dx/dy at the middle of the grid
          dx = lo%lon(2,1) - lo%lon(1,1)
          dy = lo%lat(1,2) - lo%lat(1,1)
  
          ! Handle weird edge case for WRF ideal simulations
          if (dx == 0) then
  
negative.code:                 if (.not.test_triangle(lat, lon, lo, find_surrounding, nx, ny)) then
                      find_surrounding%x = [pos%x, pos%x,   pos%x+dx, pos%x+dx ]
                      find_surrounding%y = [pos%y, pos%y+dy, pos%y,   pos%y+dy]
  
                      if (.not.test_triangle(lat, lon, lo, find_surrounding, nx, ny)) then
                          find_surrounding%x = [pos%x+dx, pos%x,   pos%x+dx, pos%x]
                          find_surrounding%y = [pos%y+dy, pos%y+dy, pos%y,   pos%y]
  
                          if (.not.test_triangle(lat, lon, lo, find_surrounding, nx, ny)) then
                              find_surrounding%x = [pos%x+dx, pos%x+dx,   pos%x, pos%x]
                              find_surrounding%y = [pos%y+dy, pos%y,   pos%y+dy, pos%y]
  