anchor.diff:  		! Back to S-EOFs
   		if(present(pc).or.present(xeof))then
   			allocate(zeof(ns,nkeep))
  -			zeof = matmul(zff, cov(:,nt:nt-nkeep+1:-1))
  +			call sgemm('N','N',ns,nkeep,nt,1.,zff,ns, &
  +				& cov(:,nt:nt-nkeep+1:-1),nt,0.,zeof,ns)
   			deallocate(cov)
   			do i = 1, nkeep
   				zeof(:,i) = zeof(:,i) / sqrt(dot_product(ww(:), zeof(:,i)**2))
  
positive.code: 	real, allocatable :: wff(:,:), ww(:), zeof(:,:), zff(:,:)
  	real, allocatable :: eig(:)
  	integer           :: zuseteof, znkeepmax, i,j
  
  	! Setups
  	! ======
  
negative.code: 				zff(:,i) = zff(:,i) * ww(:)
  			end do
  		end if
  		! pc = matmul( transpose(zff), zeof)
  		call sgemm('T','N',nt,nkeep,ns,1.,zff,ns, &
  			& zeof,ns,0.,pc,nt)
  		do i = 1, nkeep
  			pc(:,i) = pc(:,i) / dot_product(zeof(:,i)**2, ww(:))
  		end do
  