anchor.diff:    end subroutine integrate
   
     double precision function folding(ip,f,nt)
  +    use dbsym
       integer, intent(in) :: ip,nt
       interface
          function f(x,i)
  
positive.code:   end function sigmaij
  
    subroutine inomp(i,nt)
      integer, intent(in) :: i, nt
      call integrate(nstroke_coordinates(i,:),node_coordinates(i,:),i,nt)
      intphi_over(i) = folding(i,f,nt)
    end subroutine inomp
  
negative.code:     do iz=1,Nz
         do ik=1,Nk
            tmp = f(nodes(nthread,iz,ik,:),ip)
            gtmp = gtmp + realpart(jacobian(nthread,iz,ik))*tmp
         end do
      end do
      folding = gtmp
  