anchor.diff:                  ktype, skmean,  koption, itrend, nxdis,nydis,nzdis,xdb,ydb,zdb, &
                   bv,ndmax, ndmin, xa,ya,za,vra, vea, na, &
                   ndjak, xlj, ylj, zlj, vrlj, extvj, &
  -                UNEST)
  +                UNEST, idrift )
       !-----------------------------------------------------------------------
   
       !                Krige a 3-D Grid of Rectangular Blocks
  
positive.code:     ! TODO: consider rewriting this as a single block estimate and iterate in Python 
      !       use Cython, Numba and parallel computing as way to speed the process 
  
      !for safety reason we don't want undeclared variables
      IMPLICIT NONE 
  
  
      ! in 
      ! variogram
      integer, intent(in)  :: nst
      integer, intent(in), dimension(1)     :: it
      real*8, intent(in), dimension(1)      :: c0
      real*8, intent(in), dimension(nst) :: cc,aa, aa1, aa2,ang1,ang2,ang3
      ! search
      real*8, intent(in) :: radius, radius1, radius2, sang1, sang2, sang3
      integer, intent(in) :: MAXSBX, MAXSBY, MAXSBZ
      integer, intent(in) :: ndmax, ndmin
      ! grid 
      integer, intent(in) :: nx,ny,nz
      real*8 , intent(in) :: xmn,ymn,zmn,xsiz,ysiz,zsiz
      real*8 , intent(in), dimension(nx*ny*nz) :: extve
      !Jacknife
      integer, intent(in)               :: ndjak
      real*8, intent(in), dimension(nd) :: xlj, ylj, zlj, vrlj, extvj
      ! data 
      integer, intent(in)               :: nd
      real*8, intent(in), dimension(nd) :: x,y,z,vr, ve                   !vr->variable, ve->external drift
      real*8, intent(in) :: tmin, tmax
      !kriging parameters     
      integer, intent(in) :: ktype, koption, noct, itrend
      integer, intent(inout) :: nxdis,nydis,nzdis                      
      real*8, intent(in) :: UNEST
      real*8, intent(inout) :: skmean
      
  
      ! out
      ! kriging parameters
  
negative.code:     real*8 :: a,b,c,d,e,f,g,h
      dimension a(*),b(*),c(*),d(*),e(*),f(*),g(*),h(*)
  
      real*8 :: ta,tb,tc,td,te,tf,tg,th
      real*8 :: xa,xb,xc,xd,xe,xf,xg,xh
  
      ! The dimensions for lt and ut have to be at least log (base 2) n
  
      integer ::   lt(64),ut(64),i,j,k,m,p,q
   
  
      ! Initialize:
  
  