anchor.diff:            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))
  +          !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, &
  
positive.code:                (Penergy-Penergy0)/Penergy0, (Kenergy-Kenergy0)/Kenergy0, (Tenergy-Tenergy0)/Tenergy0, &
                 (Availenergy-Availenergy0)/Availenergy0, RMSdiv, maxdiv, max_gradke)
  
  		  if(errormaxrel_h > 10)then
  		     !Plot fields
  		     print*, "Stopping due to large errors", errormaxrel_h
               call plotfields(ntime, time)
               exit
            end if
  
negative.code:           stop
         end if
  
         ! For tc33, h is contant, u is rotation and the f=0
         if(testcase==33)then
            fcte=0.
            fsphere=2
            !u0=200
            h_ct=((u0**2)/2._r8)*gravi
            eta_ct=(2.*u0/erad)
         else
            h_ct=(erad*omega*u0+(u0**2)/2._r8)*gravi
            eta_ct=(2.*u0/erad+2.*omega)
         end if
  
         grad_ct=u0*eta_ct
         print*, "Test case parameters"
  
         print*, "u0", u0
         print*, "h_ct", h_ct
  
  
         do i=1, mesh%nv
            !print*, i
            h%f(i)=hollgw
            !if(testcase==32)then !Only put bottom topograpy in tc32 - geostrophic balanced flow
            bt%f(i)=h0-h_ct*dsin(mesh%v(i)%lat)**2
            !elseif(testcase==33)then
            !   bt%f(i)=h0-((u0**2)/2._r8)*gravi*dsin(mesh%v(i)%lat)**2
            !end if
  