*******************************************************************************
C                          WINDNUM_GRSOL.F
C Calculate the winding number for a (quasi) 2D system
*******************************************************************************
      subroutine windnum(x,y,image,xwind,ywind)
      include 'pimc_sqr3sc.dcl'
      real x(2*m,nat),y(2*m,nat)
      integer image(2*m,nat)

      xwind=0.0
      ywind=0.0
      do i=1,nat
         do j=1,m
            sgn=0.
            xd=x(j+1,i)-x(j,i)
            axd=abs(xd)
            if(axd.ne.0.) sgn=xd/axd
            if(axd.gt.xdim/2.) then
               sgn=-sgn
               xwind=xwind+sgn
            end if
            
            sgn=0.0
            yd=y(j+1,i)-y(j,i)
            ayd=abs(yd)
            if(ayd.ne.0.) sgn=yd/ayd
            if(ayd.gt.ydim/2.) then 
               sgn=-sgn
               ywind=ywind+sgn
            end if
         end do

         if(abs(xwind).lt.5.e-6) xwind=0.
         if(abs(ywind).lt.5.e-6) ywind=0.
      end do
      end

