anchor.diff:      ! tau_ij = -2*nu*S_ij
       do m = 1, nq
           do l = 1, nq
  -            tauij(:,:,:,l,m) = -2.0*nu_sgs*tauij(:,:,:,l,m)
  +            tauij(:,:,:,l,m) = 2.0*nu_sgs*0.5*( duidxj(:,:,:,l,m) + duidxj(:,:,:,m,l) )
           end do
       end do
   
  
positive.code:                 ctot(i,j,k) = ctot(i,j,k) / cd(i,j,k)
                  ! clip negative values
                  if ( ctot(i,j,k)<0.0 .or. isnan( ctot(i,j,k) ) ) ctot(i,j,k) = 0.0 
                  if ( ctot(i,j,k)>0.2 ) ctot(i,j,k) = 0.2 ! clip values above 0.2
              end do
          end do
      end do
  
negative.code: 
      nq = 3
  
      pi = 4.*atan(1.)
  
      eps = 1.E-10
  