anchor.diff:      C_sigma = 1.35 ! parameter from Nicoud et al.
   
       tauij = 0.0
  +    nu_sgs = 0.0
   
       do k = 1, nzp
           do j = 1, nypl
  
positive.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
  
negative.code:     ! tau_ij = -2*nu*S_ij
      do m = 1, nq
          do l = 1, nq
              tauij(:,:,:,l,m) = 2.0*nu_sgs*0.5*( duidxj(:,:,:,l,m) + duidxj(:,:,:,m,l) )
          end do
      end do
  
  