anchor.diff:          iam = ia-1
           ja = jb
           ka = kb
  -!       print '(a,f8.1,a,f8.1,a,f12.2)','uflux(ia,ja,ka,1)=',uflux(ia,ja,ka,1),' uflux(ia ,ja,ka,2)=',uflux(ia ,ja,ka,2),' dxyz=',dxyz
  -!       print '(a,f8.1,a,f8.1,a,f12.2)','vflux(ia,ja,ka,1)=',vflux(ia,ja,ka,1),' vflux(ia ,ja,ka,2)=',vflux(ia ,ja,ka,2),' dxyz=',dxyz
  -!       print *,'ia=',ia,' ja=',ja,' ka=',ka
  -!       print *,'iter=',iter,'ntractot=',ntractot
  -!         ! start track off with no error (maybe export these later to keep track of)
  +        ! start track off with no error (maybe export these later to keep track of)
           errCode = 0
   
  -! print *,'ib=',ib,' ia=',ia,' jb=',jb,' ja=',ja
  -
           call calc_dxyz(ib,jb,kb,rr,imt,jmt,km,dzt,dxdy,dxyz)
  -        ! I am putting the error checks directly in the loop for convenience
  -        !         call errorCheck('dxyzError',errCode,dxyz,flag)
  +
           ! Check the grid box volume
           if(dxyz == 0.d0) then
               print *,'====================================='
  
positive.code:             stop
          end if
  
  !         call errorCheck('infLoopError'  ,errCode)
  
          if(niter.gt.30000) then ! break infinite loops
  !         if(niter-nrj(ntrac,4).gt.30000) then ! break infinite loops
  !             nloop=nloop+1             
              print *,'====================================='
              print *,'Warning: Particle in infinite loop '
              print *,'ntrac:',ntrac
              print *,'niter:',niter!,'nrj:',nrj(ntrac,4)
  !             print *,'dxdy:',dxdy(ib,jb),'dxyz:',dxyz
  !             print *,'kmt:',kmt(ia-1,ja-1),'dz(k):',dz(ka-1)
              print *,'ia=',ia,' ib=',ib,' ja=',ja,' jb=',jb, & 
                 ' ka=',ka,' kb=',kb
              print *,'x1=',x1,' x0=',x0,' y1=',y1,' y0=',y0, & 
  
negative.code:                 ff,imt,jmt,km,rr,rb,uflux,vflux,wflux,do3d,doturb)
          endif
  
          ! === make sure that trajectory ===
          ! === is inside ib,jb,kb box    ===
          if(x1.ne.dble(idint(x1))) ib=idint(x1)+1 ! index for correct cell?
          if(y1.ne.dble(idint(y1))) jb=idint(y1)+1 ! index for correct cell?
  
          if(ia>imt .or. ib>imt .or. ja>jmt .or. jb>jmt &
               .or. ia<1 .or. ib<1 .or. ja<1 .or. jb<1) then
              print *,'====================================='
              print *,'Warning: Trajectory leaving model area'
              print *,'-------------------------------------'
              print *,'-------------------------------------'
              print *,'The trajectory is killed'
              print *,'====================================='
              errCode = -50
              flag(ntrac) = 1
          endif
          if (errCode.ne.0) cycle ntracLoop
  
          ! if trajectory above sea level,
          ! then put back in the middle of shallowest layer (evaporation)
          if( z1.ge.dble(km) ) then
              z1=dble(km)-0.5d0
              kb=km
              errCode = -50
          endif
  
          ! sets the right level for the corresponding trajectory depth
          if(z1.ne.dble(idint(z1))) then
              kb=idint(z1)+1
              if(kb == km+1) kb=km  ! (should perhaps be removed)
              errCode = -52
          endif
  
          if(x1 == dble(idint(x1)) .and. y1 == dble(idint(y1)))  then
              print *,'corner problem'
              ! corner problems may be solved the following way 
              ! but should really not happen at all
              if(ds == dse .or. ds == dsw) then
  