anchor.diff:      use interpack, only: &
       plot_scalarfield, &
       plot_cart_vectorfield, &
  +    vector_reconstruct, &
       aplyr, &
       constr
   
  
positive.code:       theta_sp = -40._r8/300._r8 
        theta_eq =  30._r8/300._r8 
        theta_np = -20._r8/300._r8 
        if(testcase==7)then
          u00 = 80.0
          lat0 = pi/7.0
          lat1 = pi/2.0 - lat0
        else
          if(testcase==8)then ! Jet in Southern Hemisphere
             u00 = 80.0
             lat0 = -5.d0*deg2rad
             lat1 = -45.d0*deg2rad
          end if
        end if
  	
        en = exp(-4/(lat1 - lat0)**2)
        umen = u00/en
        totvol = 0.0D0
  
negative.code:           p1=node(i)%edge(n)%xyz2(1,:) 
            p2=node(i)%edge(n)%xyz2(2,:) 
            do s=1,diml
               ! Reconstruct the vector field at quadrature points
               urecon = vector_reconstruct (node(i)%G(s)%lpg(n,1:3), u, mesh, reconadvmtd)
  
               !Calculando o produto interno entre o vetor velocidade e o vetor normal unitario a edge 
               dot=dot_product(urecon, node(i)%G(s)%lvn(n,1:3))
               !dot=dot_product(velocity(node(i)%G(s)%lpg(n,1:3), time),node(i)%G(s)%lvn(n,1:3))
  
               !Determinando os coeficientes da matriz de rotação para o node i 
               xx=mesh%v(i)%p(1)
               yy=mesh%v(i)%p(2)
  