subroutine forcdt(u,v,a,dp,du,time,dtime,ndag,node,
     1ug,xord,ncrv,iacflg,jnc,ipess)
c
c CONTACT: ap@marc.de
c      
c* * * * *
c
c     input of time dependent forcing functions and boundary
c     conditions.
c
c     u            total displacements at a node
c     v            velocity at a node
c     a            acceleration at a node
c     dp           load increments at a node
c     du           displacement increments at a node
c     time         time
c     dtime        time increment
c     ndag         number of degrees of freedom per node
c     node         node number
c     ug           total displacements at node in global system
c     xord         orriginal coordinates
c     ncrv         number of coordinates
c     iacflg       acceleration flag - set to 1 if accelerations given
c     jnc          increment number
c
c* * * * *
c
c PURPOSE: Follow force type 2 for point loads.
c          This routine defines point loads at a node in the
c          direction along the line connecting a point in space
c          and this node in the deformed configuration.
c          A positive force means the force is directed towards
c          the point. A negative force means the force is directed
c          away from this point.
c
c NOTE: 1) Always make incrememnt 0 a null step.
c       2) This subroutine "forcdt" works in conjunction with
c          subroutines "ubginc", "uedinc" and "nodnum".
c       3) This subroutine does not work in conjunction with
c          the "auto increment" history definition option.
c
c* * * * *
c
      implicit real*8 (a-h,o-z)                                               dp
c
      include '../common/space'
      include '../common/elmcom'
      include '../common/concom'
      include '../common/dimen'
      include '../common/prepro'
      include '../common/arrays'
      include '../common/heat'
      include '../common/array2'
      include '../common/strvar'
      include '../common/blnk'
      include '../common/lass'
c
      dimension u(ndag),v(ndag),a(ndag),dp(ndag),du(ndag)
      dimension ug(1),xord(1)
c
c* * * * *
c
c local varaiables:
c xcoor1(1..3) = coordinates of 1st node in deformed configuration.
c rcoor1(1..3) = coordinates of 1st node in undeformed configuration.
c rdisp1(1..3) = total displacements of 1st node at start of increment.
c ddisp1(1..3) = displacement changes of 1st node in current increment.
c xpnt(1..3)   = coordinates of the spacial point.
c xl12(1..3)   = unit direction vector from the node to the spacial point.
c
      dimension xcoor1(6),
     *          rcoor1(6),
     *          rdisp1(6),
     *          ddisp1(6),
     *          xpnt(6),
     *          xl12(6)
c
c* * * * *
c
c save variables below in a new common block "force_com":
c f01(1..3) = force vector at the node at start of increment.
c f11(1..3) = force vector at the node in current iteration.
c ft        = magnitude of the force at the two nodes.
c             (positive means the nodes attract each other
c              negative means the nodes repel each other).
c
      common/force_com/f01(6),f11(6),ft
c
c* * * * *
c
c define a flag for additional printing in output file.
      logical debug
      data debug/.false./
c 
      if (debug) write(6,*) 'ncycle:',ncycle
c
c define the coordinates of the spacial point.
      xpnt(1) = 0d0
      xpnt(2) = 0d0
      xpnt(3) = 0d0
c
c for a 3 dimensional problem set ncrc to 3.
c for a 2 dimensional problem set ncrc to 2.
      ncrc = 3
c
c get the internal node number for the node in nin1.
      call nodnum(node,nex1,nin1,1)
c
c get the undeformed coordinates of the node in rcoor1.
      jrdpre=0
      call vecftc(rcoor1,vars(ixord),ncrdmx,ncrd,nin1,jrdpre,2,1)
c get the total displacements at the start of the increment of the
c node in rdisp1.
      jrdpre=0
      call vecftc(rdisp1,vars(idsxts),ndegmx,ndeg,nin1,jrdpre,2,5)
c get the incremental displacements in the current increment of
c the node in ddisp1.
      jrdpre=0
      call vecftc(ddisp1,vars(idsx),ndegmx,ndeg,nin1,jrdpre,2,5)
c
c compute the coordinates of the node in the current deformed
c configuration.
      if (ncycle.eq.0) then
       do i1=1,ncrc
        xcoor1(i1) = rcoor1(i1) + rdisp1(i1)
       enddo
      else
       do i1=1,ncrc
        xcoor1(i1) = rcoor1(i1) + rdisp1(i1) + ddisp1(i1)
       enddo
      endif
c
      if (debug) then
       write(6,20) node,(rcoor1(ijk),ijk=1,ncrc)
       write(6,20) node,(rdisp1(ijk),ijk=1,ncrc)
       write(6,20) node,(ddisp1(ijk),ijk=1,ncrc)
       write(6,20) node,(xcoor1(ijk),ijk=1,ncrc)
      endif
 20   format(i5,3e14.5)
c
c compute the unit direction vector from the node to the point.
      dlx = 0d0
      do i1=1,ncrc
       xl12(i1) = xpnt(i1) - xcoor1(i1)
       dlx = dlx + xl12(i1)*xl12(i1)
      enddo
      dlx = dsqrt(dlx)
      if (dlx.gt.0d0) then
       do i1=1,ncrc
        xl12(i1) = xl12(i1)/dlx
       enddo
      endif
c
c initialize the incremental load vector.
      do i1=1,ndeg
       dp(i1) = 0d0
      enddo
c define the force increments at 1st node.
      do i1=1,ncrc
       f11(i1) = ft*xl12(i1)
       dp(i1) = f11(i1) - f01(i1)
      enddo
      if (debug) write(6,20) node,(dp(ijk),ijk=1,ncrc)
c
      return
      end
      subroutine ubginc(inc,incsub)
      implicit real*8 (a-h,o-z)
c
c     user subroutine that gets called at the start of each increment.
c
      common/force_com/f01(6),f11(6),ft
c
c initialize vectors in increment 0.
      if (inc.eq.0) then
       do i1=1,3
        f01(i1) = 0d0
        f11(i1) = 0d0
       enddo
      endif
c
c define the magnitude of the force between the nodes at the start of
c a new increment.
      df = 0.1d0
      ft = dfloat(inc)*df
c
      return
      end
      subroutine uedinc(inc,incsub)
      implicit real*8 (a-h,o-z)
c
c     user subroutine that gets called at the end of each increment.
c
      common/force_com/f01(6),f11(6),ft
c
c update the force vector f01.
c the converged force vector f11 at the end of this increment
c becomes the initial vector f01 of the next increment.
      do i1=1,3
       f01(i1) = f11(i1)
      enddo
c
      return
      end
      subroutine nodnum(noid,noex,noin,ic)
c
      implicit real*8 (a-h,o-z)
c
c     ic=1   determine noex and noin for noid
c     ic=2   determine noin and noid for noex
c     ic=3   determine noid and noex for noin
c
c     noid = noid id (given by user)
c     noex = marc node number before optimize
c     noin = marc node number after optimize
c
      include '../common/dimen'
      include '../common/space'
      include '../common/arrays'
      include '../common/develp'
      include '../common/prepro'
c
      if(ic.eq.1) then
       noex=noid
       if(nnoids.ne.0) noex=ibsrch(noid,ints(inoids),numnp,1)
       noin=noex
       if(joptit.ne.0) noin=igetsh(ints(inpnum+noex-1),0)
      else if(ic.eq.2) then
       noid=noex
       if(nnoids.ne.0) noid=ints(inoids+noex-1)
       noin=noex
       if(joptit.ne.0) noin=igetsh(ints(inpnum+noex-1),0)
      else if(ic.eq.3) then
       noex=noin
       if(joptit.ne.0) noex=igetsh(ints(inpnum+noin-1),1)
       noid=noex
       if(nnoids.ne.0) noid=ints(inoids+noex-1)
      endif
c
      return
      end