c force_tor.f
c add the streching and bending forces to gc(i)

c force_tor(n,gc,ntime)
c gc(maxn) using to add the force_tor vector
c refurbish force_tor.f only i=1..n is used now.
c more details check HM's p52-53

	subroutine force_tor(n,gc)
	
	implicit double precision(a-h, o-z)
        parameter(maxa=500,maxn=maxa*3)
 	double precision gc(maxn),Theta(maxa,3)
        common/bf_coords/f(maxa,3),v(maxa,3),u(maxa,3)
        common/euler/alpha(maxa),beta(maxa),gama(maxa),phi(maxa)
	common/tor_param/dr,cg,phi0,dlk,phii,subseg
	common/matrix2/ r(maxa,maxa)
	
c update Theta(j)   (3.38)
	do i=1,n
	  ii=i+1
	  ii2=i-1
	  if (i.eq.1) then
	   ii2=n
 	  endif
	  if (i.eq.n) then
	   ii=1
 	  endif
	 do j=1,3
	 Theta(i,j)=phi(i)*tan(beta(i)/2)/r(i,ii)*(cos(alpha(i))*v(i,j)
     + -sin(alpha(i))*f(i,j))+phi(ii2)*tan(beta(ii2)/2)/r(i,ii)*
     + (cos(gama(ii2))*v(i,j)+sin(gama(ii2))*f(i,j))
	 enddo
	enddo

c add torsional force onto gc(maxn)
 	ii=0
	do i=1,n
	 i3=i-1
	    if (i.eq.1) then
	     i3=n
	    endif
	  do j=1,3
	    ii=ii+1 
c	    ii=3*i+j-3
	    gc(ii)=gc(ii)+cg*(Theta(i,j)-Theta(i3,j))
	   enddo
	enddo

      return
      end