anchor.diff:    end subroutine laplacian_truncation
   
   
  +  subroutine poisson_error(mesh)
  +
  +    !Mesh
  +    type(grid_structure) :: mesh
  +
  +
  +    !Scalar fields
  +    type(scalar_field):: func
  +    type(scalar_field):: lap_ex
  +    type(scalar_field):: lap
  +    type(scalar_field):: error
  +    type(scalar_field):: ergrad
  +    type(scalar_field):: eddif
  +
  +    !Counters
  +    integer:: i
  +    integer:: j
  +    integer:: k
  +
  +    !Errors
  +    logical:: ifile
  +    integer:: iunit
  +    real (r8)::  error2
  +    real (r8):: errormax
  +    real (r8):: errorgrad2
  +    real (r8):: errorgradsup
  +    real (r8):: maxeddif
  +    character (len=256):: filename
  +
  +    !Aux
  +    real(r8):: p(1:3)
  +    real(r8):: gradex
  +    real(r8):: gradest
  +    real(r8):: len
  +
  +
  +    print*
  +    print*,"Laplacian Truncation Analysis "
  +    print*
  +
  +    !-------------------------------------------
  +    !  Read parameters from file "simul.par"
  +    !------------------------------------------
  +
  +    call getsimulpars(mesh)
  +
  +
  +    print*,"Setting up variables ..."
  +
  +    !Scalars on hexagon centers (nodes) - Exact
  +    lap_ex%pos=0
  +    lap_ex%n=mesh%nv
  +    lap_ex%name="lap_ex"
  +    allocate(lap_ex%f(1:lap_ex%n))
  +    do i=1, mesh%nv
  +      !Laplacian on node
  +      p=mesh%v(i)%p
  +      !lap_ex%f(i)=lap_exact(p)
  +
  +      !Laplacian on barycenter
  +      !p=mesh%hx(i)%b%p
  +      lap_ex%f(i)=lap_exact(p)
  +       !print "(i8, 4f16.8)", i, lap_ex%f(i), p
  +    end do
  +
  +    !Edge displacement
  +    eddif%n=mesh%nv
  +    eddif%name="edgedif"
  +    eddif%pos=0
  +    allocate(eddif%f(1:mesh%nv))
  +    do i=1, mesh%nv
  +      eddif%f(i)=0._r8
  +      do j=1,mesh%v(i)%nnb
  +        k=mesh%v(i)%ed(j)
  +        eddif%f(i)=max(eddif%f(i), (arclen(mesh%ed(k)%c%p, mesh%edhx(k)%c%p)/ &
  +          mesh%edhx(k)%leng))
  +      end do
  +       !eddif%f(i)=eddif%f(i)/mesh%v(i)%nnb
  +    end do
  +
  +    !Edge intersection maximum difference
  +    maxeddif=maxval(abs(eddif%f(1:eddif%n)))
  +
  +    !Gradient error
  +    ergrad%n=mesh%nv
  +    ergrad%name="ergrad"
  +    ergrad%pos=0
  +    allocate(ergrad%f(1:mesh%nv))
  +    do i=1, mesh%nv
  +      ergrad%f(i)=0._r8
  +    end do
  +
  +    !Laplacian error
  +    error%n=mesh%nv
  +    error%name="erlap"
  +    error%pos=0
  +    allocate(error%f(1:mesh%nv))
  +
  +    !Test function - Scalars on hexagon centers (nodes)
  +    func%pos=0
  +    func%n=mesh%nv
  +    func%name="func"
  +    allocate(func%f(1:func%n))
  +    do i=1, mesh%nv
  +      !Function on node
  +      p=mesh%v(i)%p
  +      func%f(i)=f(p)
  +    end do
  +
  +    !Calculate Numerical Laplacian
  +    lap%pos=0
  +    lap%n=mesh%nv
  +    lap%name="lap_est"
  +    allocate(lap%f(1:lap%n))
  +    do i=1, mesh%nv
  +      lap%f(i)=0._r8
  +      ergrad%f(i)=0._r8
  +      do j=1, mesh%v(i)%nnb
  +        !Edge index
  +        k=mesh%v(i)%ed(j)
  +        !Hexagonal edge midpoint
  +        p=mesh%edhx(k)%c%p
  +        !ExactGrad=ExactGradVector*NormalVectorEdge*CorrectionOutHx
  +        gradex=dot_product(df(p),mesh%edhx(k)%nr)*mesh%hx(i)%nr(j)
  +        !Estimated Gradient Normal component
  +        gradest=func%f(mesh%v(i)%nb(j))-func%f(i)
  +        len=mesh%ed(k)%leng
  +        !arclen(p, mesh%v(i)%p) + &
  +        !arclen(p, mesh%v(mesh%v(i)%nb(j))%p)
  +        gradest=gradest/ len !mesh%ed(k)%leng
  +        !Maximum gradient error for cell
  +        ergrad%f(i)=max(abs(gradest-gradex), ergrad%f(i))
  +        !Updade Laplacian
  +        lap%f(i)=lap%f(i)+gradest*mesh%edhx(k)%leng
  +      end do
  +      lap%f(i)=lap%f(i)/mesh%hx(i)%areag
  +      !Error in laplacian
  +      error%f(i)=lap_ex%f(i)-lap%f(i)
  +       !print"(i4, 3f16.8)", i, lap%f(i), lap_ex%f(i), error%f(i)
  +    end do
  +
  +    !Global Errors
  +
  +    error2=error_norm_2(lap%f, lap_ex%f, error%n)
  +    errormax=error_norm_max(lap%f, lap_ex%f, error%n)
  +    print*, "Error Lap (max, L2): ", errormax, error2
  +
  +    !Global Errors for gradients
  +    errorgrad2=dsqrt(dot_product(ergrad%f,ergrad%f)/ergrad%n)
  +    errorgradsup=maxval(abs(ergrad%f(1:ergrad%n)))
  +    print*, "Error GRAD (max, L2): ", errorgradsup, errorgrad2
  +    print*
  +
  +    ! Save error estimates
  +    !-------------------------------------------------------
  +
  +    !For hexagonal methods
  +    filename=trim(datadir)//"lap_errors_"//trim(mesh%kind)// &
  +      "_"//trim(mesh%pos)//"_"//trim(mesh%optm)//".txt"
  +    call getunit(iunit)
  +
  +    inquire(file=filename, exist=ifile)
  +    if(ifile)then
  +      open(iunit,file=filename, status='old', position='append')
  +    else
  +      open(iunit,file=filename, status='replace')
  +      write(iunit, '(a)') &
  +        " n    distance stag testfunc errormax error2 errorgradsup errorgrad2 maxeddif"
  +    end if
  +
  +    write(iunit, '(i8, f18.8, i8, 5f22.12)') mesh%nv, mesh%meanvdist*rad2deg, &
  +      testfunc, errormax, error2, errorgradsup, errorgrad2, maxeddif
  +
  +    close(iunit)
  +
  +    ! Plot fields
  +    !-------------------------
  +    if(plots) then
  +      print*, "Plotting variables ... "
  +
  +      lap_ex%name=trim(simulname)//"_exact"
  +      call plot_scalarfield(lap_ex, mesh)
  +
  +      lap%name=trim(simulname)//"_est"
  +      call plot_scalarfield(lap, mesh)
  +
  +      error%name=trim(simulname)//"_error"
  +      call plot_scalarfield(error, mesh)
  +
  +      eddif%name=trim(simulname)//"_eddif"
  +      call plot_scalarfield(eddif, mesh)
  +
  +      ergrad%name=trim(simulname)//"_ergrad"
  +      call plot_scalarfield(ergrad, mesh)
  +
  +    end if
  +
  +    return
  +
  +  end subroutine poisson_error
   
     subroutine relaxac(mesh)
       !--------------------------------------------------------
  
positive.code: 
    !Global constants
    use constants, only: &
         datadir, &
         deg2rad, &
         eps, &
         erad, &
         i4, &
         pardir, &
         pi, &
         r8, &
         rad2deg
  
  
    !Data structures
    use datastruct, only:  &
         grid_structure, &
         scalar_field, &
         vector, &
         vector_field_cart
  
    !Spherical mesh routines
    use smeshpack, only: &
         alignind, &
         alignindlimit, &
         arcdistll, &
         arcintersec, &
         arclen, &
         cart2sph, &
         choleskydecomp, &
         choleskysolve, &
         convert_vec_sph2cart, &
         convijll, &
         convllij, &
         distortionhx, &
         error_norm_2, &
         error_norm_max, &
         gcircarcintersec, &
         getnearnode, &
         getnearnodes, &
         gettr, &
         getunit, &
         hxtr_intersec_areas, &
         modint, &
         norm, &
         ortogonalarc, &
         proj_vec_sphere, &
         sph2cart, &
         sphpolarea, &
         sqtriintersec, &
         vorbarycenter
  
  
  
    !Interpolation pack
    use interpack, only: &
         plot_scalarfield, &
         plot_cart_vectorfield, &
         plot_grad_vectorfield, &
         plot_scalarfield_sphgrid, &
         plot_vectorfield_sphgrid
  
    implicit none	
  
  
negative.code: 
      !Set ploting parameter
      if(iplots==1)then
        plots=.true.
      else
        plots=.false.
      end if
  
  
      !Set discretization position (centroid?)
      if(idiscretcentroid==1)then
        discretcentroid=.true.
      else
        discretcentroid=.false.
      end if
  
  
  