anchor.diff:      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)
  +      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
   
  
positive.code: 
      !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
  
  
negative.code: 
      !Laplacian
      select case(testfunc)
        case(4) !Smooth Trigonometric
          lap_exact=dcos(lat) * (-5.0_r8 * dcos(lon) ** 2 + 7.0_r8 &
            - 12.0_r8 * dcos(lat)** 2 * dsin(lon) ** 2)
        case(7) !Constant 1
          lap_exact=0
        case (8) !Oscilat trig (beta of HR95)
          m=1
          n=1
          lap_exact = -dcos(m*lon)*dcos(n*lat)**2 * (m**2*dcos(n * lat) **2 &
            - 4._r8 * dsin(lat) * dcos(n*lat)*dsin(n*lat)*n*dcos(lat) &
            - 12._r8 * dcos(lat)** 2 * n** 2 + 16._r8*dcos(lat)**2 *dcos(n* lat)**2* n** 2) &
            / dcos(lat) ** 2
           !lap_exact= (dcos(lat)**2) * (dsin(lon) **5 )* &
           !    (33._r8 * dcos(lon) ** 2 + 9.0_r8 &
           !    - 20.0_r8 * dcos(lat)** 2 * dsin(lon) ** 2)
        case default
          print*, "Laplacian (lap_exact) error: unknown scalar field (testfunc) :", testfunc
          stop
      end select
  
      return
  