c  **************************************************************
      subroutine cal_writhe(n,index,n108)

      implicit double precision(a-h, o-z)
      parameter(maxa=500,pi=3.1415926535d0,mbu=200,pi2sq=pi*pi*2.d0)
      parameter(pi2=pi*2.d0)
      integer n,i,j,index
      dimension ia1(2)
        common/matrix/x(maxa,3),bangle(maxa)
       common/matrix2/ r(maxa,maxa)
       common/figinfo/gyrat,endl2,zsum1,cmass(3),r1,shapef,writhe
       common/ischeme/icircular,iele,ihydro,iclose,iknot
     +,neuler
       common/ddx/dx(maxa),dy(maxa),dz(maxa)
       common/bf_coords/f(maxa,3),v(maxa,3),u(maxa,3)
       common/xoldtw/olddlk,xold(maxa,3),fold(maxa,3),vold(maxa,3)
     +,uold(maxa,3)   
c debug
      do i = 1, 3
         x(n+1,i) = x(1,i)
      enddo

      do i = 1, n
         dx(i)=x(i+1,1)-x(i,1)
         dy(i)=x(i+1,2)-x(i,2)
         dz(i)=x(i+1,3)-x(i,3)
      enddo

       call kwr(n,topl,jwrold,ierr,0)
       call bwr(n,2,n,beeta)

      writhe = jwrold + beeta
c      in = 1
c      ik = n
c      call calwr(n,in,ik,beta2)
c      writhe2 = beta2
      

        call cal_gyrat(n)
        call cal_twist(n,twist)

c part to check if strand-crossing in past iwr steps (usually 100)
c  if so, output the related two configurations to fort.219
c                       and the dtime of crossing to fort.218
c  and reset the old configurations, changed the random seeds, and
c  continued uncrossed configurations
      if (abs(twist+writhe-olddlk).ge.0.9) then
         write(218,*) index,twist,writhe
         call backnocross(n,index)
         call flush(218)
	 index=index-n108
c	write(*,*) 'reset at',index
      else
c         write(*,*) 'no crossing :',index,twist,writhe,olddlk
        if (iknot.eq.1) then
	       call kpoly(n,ia1,ierr)
       	       write(108,13) index, writhe,gyrat,twist,ia1(1),ia1(2)
        else
          write(108,11) index, writhe,gyrat,twist,cmass(1),cmass(2)
     +,cmass(3)
        endif
        call flush(108)

         do i=1,n
            do j=1,3
               xold(i,j)=x(i,j)
               fold(i,j)=f(i,j)
               vold(i,j)=v(i,j)
               uold(i,j)=u(i,j)
            enddo
         enddo   
      
      endif
         

      
13    format(i8,3(1x,g11.5),1x,i8,1x,i8)	    
11    format(i8,6(1x,g11.5))
12    format(i8,1x,g11.5,1x,g11.5)

      return
      end