anchor.diff:                      dz = 0.5*(zpts(k-1)-zpts(k+1))
                   end if
   
  -                delta = (dx*dy*dz)**(2./3.)      
  +                delta = (dx*dy*dz)**(1./3.)      
   
   !                 if(myid==0 .and. debug>=2 .and. mod(istep,iomod)==0) then
                   if (D_sigma<-1.E-10 .or. isnan(D_sigma)) then
  
positive.code: 
      nq = 3
  
      ! let tauij = S_ij
      do m = 1, nq
          do l = 1, nq
              tauij(:,:,:,l,m) = 0.5*( duidxj(:,:,:,l,m) + duidxj(:,:,:,m,l) )
          end do
      end do
  
      pi = 4.*atan(1.)
  
      eps = 1.E-10
  
negative.code:                 if (D_sigma<0.0 .or. isnan(D_sigma)) D_sigma = 0.0
  
                  ! compute nu_sgs = (C*delta)^2 D_sigma
                  nu_sgs(i,j,k) = -(C_sigma*delta)**2*D_sigma 
  
              end do
          end do
  