anchor.diff:      real (kind=dp_t), intent(in) :: t
       integer, intent(in) :: n
   
  -    real (kind=dp_t) :: mach_max
  +    real (kind=dp_t) :: mach_max, mach_L2, M
   
       real (kind=dp_t) :: e, rho, p
   
  
positive.code:        open(newunit=lun, file=trim(outfile), status="old", position="append")
      endif
  
  100 format(1x, g20.10, 1x, g20.10)
      write(unit=lun,fmt=100) t, mach_max
  
      close (unit=lun)
  
  
negative.code:        rho = U%data(i,iudens)
         call eos(eos_input_e, p, e, rho)
  
         M = abs(U%data(i,iumomx)/U%data(i,iudens))/sqrt(gamma*p/rho)
  
         mach_max = max(mach_max, M)
         mach_L2 = mach_L2 + M**2
      enddo
  
      mach_L2 = sqrt(U%grid%dx*mach_L2)
  
      outfile = trim(problem_name) // "_diag.out"
  
      if (n == 1) then
  