anchor.diff:      end if
   
       if(k==0 .or. mod(k,nprints*40)==0 )then
  -        buffer="        n    mvdist  tcase       k       nt    dt(s)  "//&
  -         "    time(dys) "//&
  -         "     mass "//&
  -         "        Penergy          Kenergy       Tenergy        Availenergy       RMSdiv"
  +       buffer="        n    mvdist  tcase       k       nt    dt(s)  "//&
  +            "    time(dys) "//&
  +            "     mass "//&
  +            "        Penergy          Kenergy       Tenergy        Availenergy       RMSdiv"
          write(*, '(a)') trim(buffer)
       end if
   
  
positive.code:           maxdiv=maxval(abs(divu%f))
            max_gradke=0.0_r8
      	  do l=1, mesh%ne
           	!print*, l
  	         !Calculate gradient of Ke on edges
  	         if(useStagHC)then
  	            signcor=dsign( 1._r8, dot_product(mesh%edhx(l)%nr, &
  	                 mesh%v(mesh%edhx(l)%sh(2))%p-mesh%v(mesh%edhx(l)%sh(1))%p ))
  	            gradke_tmp=signcor*(Kin_energy%f(mesh%edhx(l)%sh(2))-Kin_energy%f(mesh%edhx(l)%sh(1)))
  	         elseif(useStagHTC)then
  	          ! Obs: the tangent of a triangle edge (which is used as normal
  	          !   of the voronoi cell edges) is always defined such
  	          !   that it point from its vertex 1 to its vertex 2
  	          !   Therefore, the gradient needs no correction term (n_ei)
  	            gradke_tmp=(Kin_energy%f(mesh%ed(l)%v(2))-Kin_energy%f(mesh%ed(l)%v(1)))
  	          end if
  	          max_gradke=max(max_gradke, abs(gradke_tmp/mesh%ed(l)%leng/erad))
         	  end do
  
            call write_evol_error(k, time, (Tmass-inimass)/inimass, errormaxrel_h, error2_h, &
  
negative.code:     !File for errors
      filename=trim(datadir)//trim(swmname)//"_mass_energy_"//trim(mesh%name)//".txt"
      buffer="        n        mvdist  tcase       k       nt        dt(s)       "//&
           "          cfl                 time(dys)             mass               Penergy    "//&
           "          Kenergy            Tenergy            Availenergy           RMSdiv"
  
      inquire(file=filename, opened=iopen)
      if(iopen)then
  