anchor.diff:  
   
       !Variables at Voronoi centers
  -    do i=1, mesh%nv
  -        lon = mesh%v(i)%lon
  -        lat = mesh%v(i)%lat
  -        call sph2cart(lon, lat, p(1), p(2), p(3))
  -        call sph2cart(0._r8, 0._r8, pc(1), pc(2), pc(3))
  -        tracer%f(i) = dexp(-5._r8*norm(p-pc)**2)
  -        htracer%f(i) = tracer%f(i)
  -        tracer_exact%f(i) = tracer%f(i)
  -        tracer_error%f(i) = tracer%f(i)
  -    end do
  +    if (order > 1) then
  +        do i=1, mesh%nv
  +            tracer%f(i) = node(i)%phi_new2
  +            htracer%f(i) = tracer%f(i)
  +            tracer_exact%f(i) = tracer%f(i)
  +            tracer_error%f(i) = tracer%f(i)
  +            !print*, node(i)%phi_new2 - dexp(-5._r8*norm(p-pc)**2)
  +        end do
  +    
  +    else
  +        do i=1, mesh%nv
  +            lon = mesh%v(i)%lon
  +            lat = mesh%v(i)%lat
  +            call sph2cart(lon, lat, p(1), p(2), p(3))
  +            call sph2cart(0._r8, 0._r8, pc(1), pc(2), pc(3))
  +            tracer%f(i) = dexp(-5._r8*norm(p-pc))**2
  +            htracer%f(i) = tracer%f(i)
  +            tracer_exact%f(i) = tracer%f(i)
  +            tracer_error%f(i) = tracer%f(i)
  +        end do
  +    endif
   
       !Calculate derived initial fields
       call tendency_moist_swm(h, u, htheta, hqv, hqc, hqr, htracer, masseq%f, momeq%f, tempeq%f, vapoureq%f, &
  
positive.code: 
        !Calculate erngies
        call calc_energies(Penergy, Kenergy, Tenergy, Availenergy)
        print '(a33, 2e16.8)','Change in mass of h*(total water):', (Twater-iniwater)/iniwater
        print*,''
        if (maxval(htracer%f)>10.d0)then
            stop
        endif
        !update fields
        u_old=u
        h_old=h
  
negative.code:     type(scalar_field), intent(in):: h  !General
  
      !Velocities (defined on edges - only normal component)
      type(scalar_field), intent(inout):: u  !General
  
      !Temperature (defined on voronoi centers) 
      type(scalar_field), intent(in):: htheta  !General
  