anchor.diff:      !!
       !!------------------------------------------------------------
       integer function minyw(yw,latgrid,xpos,ypos,lat)
  -        integer, intent(in)::yw,xpos,ypos
  -        real,intent(in)::latgrid(:,:),lat
  -        real::curdist,dx
  -        integer::ny
  -
  -        ny=size(latgrid,2)
  -        curdist=abs(lat-latgrid(xpos,ypos))
  -        if (ypos<ny) then
  -            dx=abs(latgrid(xpos,ypos)-latgrid(xpos,ypos+1))
  +        integer, intent(in) :: yw, xpos, ypos
  +        real,    intent(in) :: latgrid(:,:), lat
  +        real    :: curdist, dx
  +        integer :: jme
  +
  +        jme = ubound(latgrid, 2)
  +
  +        curdist = abs(lat - latgrid(xpos,ypos))
  +        if (ypos < jme) then
  +            dx = abs(latgrid(xpos, ypos) - latgrid(xpos, ypos+1))
           else
  -            dx=abs(latgrid(xpos,ypos)-latgrid(xpos,ypos-1))
  -        endif
  -        minyw=max(yw,ceiling(curdist/dx))
  -        if (minyw>yw) then
  -            if (minyw>4) then
  -                minyw=ceiling(minyw*0.8)
  +            dx = abs(latgrid(xpos,ypos)-latgrid(xpos,ypos-1))
  +        endif
  +        minyw = max(yw, ceiling(curdist/dx))
  +        if (minyw > yw) then
  +            if (minyw > 4) then
  +                minyw = ceiling(minyw*0.8)
               endif
           endif
       end function minyw
  
positive.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::nx,nz,ny
          integer:: i,j,k,l,localx,localy
          logical :: boundaries
          real::localw
          real :: local_center
  
          nx=size(fieldout,1)
          nz=size(fieldout,2)
          ny=size(fieldout,3)
  
          boundaries = .False.
          if (present(boundary_only)) boundaries = boundary_only
  
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]
  