anchor.diff:      divu%name="divu"
       allocate(divu%f(1:divu%n), stat=ist)
   
  +    !Gradient of divergence
  +    grad_ed_div%n=mesh%ne
  +    grad_ed_div%pos=edpos
  +    allocate(grad_ed_div%f(1:grad_ed_div%n), stat=ist)
  +    grad_ed_div%f = 0._r8
  +
  +    !Gradient of vorticity
  +    grad_ed_vort%n=mesh%ne
  +    grad_ed_vort%pos=edpos
  +    allocate(grad_ed_vort%f(1:grad_ed_vort%n), stat=ist)
  +    grad_ed_vort%f = 0._r8
  +
  +    !Laplacian at edges
       lapu%n=mesh%ne
       lapu%pos=edpos
       lapu%name="lapu"
  
positive.code:   type(scalar_field):: q_ed  !pv on edges
    type(scalar_field):: q_ed_exact  !pv on edges
    type(scalar_field):: q_ed_error  !pv on edges
  
  
    type(scalar_field):: q_grad_ed  !grad pv on edges
  
negative.code:     !$OMP SHARED(u_error, u_exact, uh, uhq_perp) &
      !$OMP SHARED(v_hx, v_ed, vhq_hx) &
      !$OMP SHARED(h, h_old, h_0, h_error, h_exact, h_ed, h_rhb) &
      !$OMP SHARED(h_tr, zeta, eta, eta_ed, q_tr, q_ed) &
      !$OMP SHARED(ke_hx, ke_tr, ghbK, grad_ghbK, hbt, bt) &
      !$OMP SHARED(grad_ed_div, grad_ed_vort) &
      !$OMP SHARED(divuh, q_tr_exact, q_tr_error, q_grad_ed, divu, lapu ) &
      !$OMP SHARED(massf0, massf1, massf2, massf3) &
      !$OMP SHARED(momf0, momf1, momf2, momf3)
  