c ********************************************************************
c  brownian dynamics algorithm. taking into account the hydrodynamics.
c  revised second-order BD (less update D).
c     hyd0   -- modifed scheme  (use old D)
c ********************************************************************

      subroutine bd_hy0(n,ntime) 

      implicit double precision(a-h, o-z)
      parameter(maxa=500,maxn=maxa*3,bolz=8.31d-3/4.18)
      integer nout, nn,ntime,n3,n
      double precision gc(maxn),xold(maxa,3),gcnew(maxn),fcnew,
     +   y(maxn),rran(maxn),fold(maxa,3),vold(maxa,3),uold(maxa,3)
     +	 ,ftor(maxa),phiold(maxa),sum(maxn)
      double precision drr(maxa,3)
      common/matrix/x(maxa,3),bangle(maxa),rz(maxa)
      common/param/temp,bendc,r0,beadr,ss1,interval,dtime,gamma,
     +amass,am(maxn),coef,dtime2
      common/diffu_tensor/diffu(maxn,maxn),w(maxn,maxn),varr(maxn)
      common/myrandom/amean,var(maxn),fvar(maxn)
      common/euler/alpha(maxa),beta(maxa),gama(maxa),phi(maxa)
      common/bf_coords/f(maxa,3),v(maxa,3),u(maxa,3)
      common/ischeme/icircular,iele,ihydro,iclose,iknot
     +,neuler

       n3=n*3
c       call update_rij(n)

       do i = 1, n
         xold(i,1) = x(i,1)
         xold(i,2) = x(i,2)
         xold(i,3) = x(i,3)
         if (icircular .eq. 1) then
            do j = 1, 3
               fold(i,j) = f(i,j) 
               uold(i,j) = u(i,j) 
               vold(i,j) = v(i,j)
            enddo
            phiold(i) = phi(i)
         endif
       enddo

      call force4(n,gc,ntime)
      if ((icircular .eq. 1).and.(iclose.eq.1)) then
          iout = 0
	  call torsion_bd(n,ftor,phiold,iout)
      endif

      call rannv2(n3,y,amean,varr)

      do i = 1, n3
         rran(i) = 0.d0
         do j = 1, i
            rran(i) = rran(i) + w(i,j)*y(j)
         enddo
      enddo

      do i = 1, n3
         sum(i) = 0.d0
      enddo
      do i=1, n3
         do j = 1, n3
            if (i .ge. j) then
               sum(i) = sum(i) + diffu(j,i)*gc(j)
            else
               sum(i) = sum(i) + diffu(i,j)*gc(j)
            endif
         enddo
      enddo
c in my model  i=3*id+jd-3    jd=1:3  id=1:n
        ii=1
        do i=1,n
           do j=1,3
           x(i,j)=x(i,j)+dtime*sum(ii)+rran(ii)
           ii=ii+1
           enddo
        enddo

	call update_rij(n)	
        call update_vector(n)

      if ((icircular .eq. 1).and.(iclose.eq.1)) then
         call update_euler(n)
      endif

C Go further to use second-order BD.
c      update dphiz(rz),f,v by ftor 

 
      call force4(n,gcnew,ntime)

      if ((icircular .eq. 1).and.(iclose.eq.1)) then
         do i = 1, n
            do j = 1, 3
               f(i,j) = fold(i,j)
               v(i,j) = vold(i,j)
               u(i,j) = uold(i,j)
            enddo
         enddo
         iout = 1
         call torsion_bd(n,ftor,phiold,iout)
      endif

      do i = 1, n3
         sum(i) = 0.d0
      enddo
      do i=1, n3
         do j = 1, n3
            if (i .ge. j) then
               sum(i) = sum(i) +diffu(j,i)*(gc(j)+gcnew(j))
            else
               sum(i) = sum(i) +diffu(i,j)*(gc(j)+gcnew(j))
            endif
         enddo
      enddo
c in my model  i=3*id+jd-3    jd=1:3  id=1:n
        ii=1
        do i=1,n
           do j=1,3
           x(i,j)=xold(i,j)+0.5*dtime*sum(ii)+rran(ii)
           ii=ii+1
           enddo
        enddo

c	do i=1,n
c           write(99,700)i, drr(i,1),drr(i,2),drr(i,3)
c        enddo
c

      call update_rij(n)
      call update_vector(n)
      if ((icircular .eq. 1).and.(iclose.eq.1)) then
         call update_euler(n)
      endif
700	format(i8,3(1x,g13.6))
      

      return 
      end