anchor.diff:      kb = kstart(ntrac)
     
       ! ===  start loop for each trajectory ===
  -    !         scrivi=.true.
       niterLoop: do 
           niter=niter+1 ! iterative step of trajectory
  -!         print *,'niter=',niter
           rg = ts 
  -!         rg=dmod(ts,1.d0) ! time interpolation constant between 0 and 1
           rr = 1.d0-rg
  -!         print *,'ts=',ts,' rg=',rg,' rr=',rr
           ! Update particle indices and locations
           x0 = x1
           y0 = y1
  
positive.code: 
          call calc_time(ds,dsmin,dt,dtmin,tss,tseas,ts,tt,dxyz,dstep,iter,rb,dsc)
  
  !         print *,'After calc_time: ds=',ds,' tss=',tss,' ts=',ts,' tt=',tt,' dt=',dt,' dtmin=',dtmin
  
          ! === calculate the new positions ===
          ! === of the trajectory           ===  
  !         if(ntrac==42) then
  !              print *,''
  !            print *,'before pos'  
  !             print '(a,f6.2,a,f6.2,a,f6.2,a,f6.2,a,f6.2,a,f6.2)','x0=',x0,' x1=',x1,&
  !                 ' y0=',y0,' y1=',y1,' z0=',z0,' z1=',z1
  !             print '(a,i3,a,i3,a,i3,a,i3,a,i3,a,i3)','ia=',ia,' ib=',ib,&
  !                 ' ja=',ja,' jb=',jb,' ka=',ka,' kb=',kb
  !         endif
  
          if(doturb==1) then
              call pos(ia,ja,ka,ib,jb,kb,x0,y0,z0,x1,y1,z1,ds,dse,dsw,dsn,dss,dsu,dsd,dsmin,dsc,&
                  ff,imt,jmt,km,rr,rb,uflux,vflux,wflux,do3d,doturb,upr)
  
negative.code:             errCode = -54
           endif
  
          ! === diffusion, which adds a random position ===
          ! === position to the new trajectory          ===
          if(doturb==2 .or. doturb==3) then
              call diffuse(x1, y1, z1, ib, jb, kb, dt,imt,jmt,km, dxv,dyu,dzt,h,ah,av,do3d,doturb)
          endif
  
          ! === Optional periodic boundary conditions ===
          ! If activated, a drifter is moved from the appropriate domain edge and 
  