anchor.diff:  
     pure subroutine interpolmap(map, gt, lat, lng, mout)
       real, intent(in), dimension(:,:) :: map
  +    type(geotransform), intent(in) :: gt
       real(dp), intent(in) :: lat, lng
  -    type(geotransform), intent(in) :: gt
       real(dp), intent(out) :: mout
       integer :: ki, kj
       real(dp) :: ri,rj
  
positive.code:       lngm(i,j) = (lng1 + lng2) / 2 + real(j - 1, dp) / p - 180
      end forall
  
      !$omp parallel sections private(mask)
  
      !$omp section
      mask(:,:) = angdist(latm, lngm, lat1, lng1) .le. (dmax / radearth)
      call growmask(mask)
      llims(1,1) = minval(latm, mask)
  
negative.code:     llims(3,2) = minval(lngm, mask)
      ulims(3,2) = maxval(lngm, mask)
  
      mask(:,:) = angdist(latm, lngm, lat2, lng1) .le. (dmax / radearth)
      call growmask(mask)
      llims(4,1) = minval(latm, mask)
      ulims(4,1) = maxval(latm, mask)
      llims(4,2) = minval(lngm, mask)
      ulims(4,2) = maxval(lngm, mask)
  
      llat = minval(llims(:,1))
      llng = minval(llims(:,2))
  