subroutine fourstep(t,dt,fm,lm,ifirst,ilast,
     &  xold, xnew, km, km1, km2, km3, Gforce)
      IMPLICIT NONE
c     ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      include "./linearad.i"
c     ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      integer fm,lm,ifirst,ilast
      double precision t,dt
      double precision xold(fm:lm,1:4), xnew(fm:lm,1:4),
     &                 km(fm:lm,1:4), km1(fm:lm,1:4),
     &                 km2(fm:lm,1:4), km3(fm:lm, 1:4),
     &                 Gforce(fm:lm,1:2)
      integer ie, je
      double precision force(1:4), force1, force2
c     ******************************************************************
c     multiply fluxes by dt (ie, compute time integral over cell side)
c     ******************************************************************

      do ie=ifirst,ilast
         do je=1,4
            xnew(ie,je)=xold(ie,je)
         enddo
      enddo
      
      do ie=ifirst,ilast
         call force_law(fm, lm, ifirst, ilast, ie, xold, force,
     &                  force1, force2)
         Gforce(ie,1) = force1
         Gforce(ie,2) = force2

         do je=1,4
            km(ie,je) = force(je)
            xnew(ie,je) = xnew(ie,je) 
     &                   + dt/24.d0*( 55.d0*km(ie,je)
     &                               -59.d0*km1(ie,je)
     &                               +37.d0*km2(ie,je)
     &                               - 9.d0*km3(ie,je) )
         enddo
      enddo

      return
      end

c*************************************************************************
      subroutine force_law(fm, lm, ifirst, ilast, ie, x, force,
     &                     force1, force2)
      IMPLICIT NONE
c     ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      include "./linearad.i"
c     ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      integer fm,lm,ifirst,ilast, ie
      double precision x(fm:lm,1:4), force(1:4)
      double precision force1, force2
      integer je
      double precision vlen, xlen
      double precision cla,clr
      double precision temp1, temp2
c     ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      cla=ca/la
      clr=cr/lr

      vlen=(x(ie,3)**2+x(ie,4)**2)**0.5d0
      force(1)=0.d0
      force(2)=0.d0
      force(3)=0.d0
      force(4)=0.d0

      force1 = 0.d0
      force2 = 0.d0

      do je=ifirst,ilast
         xlen=( (x(ie,1)-x(je,1))**2
     &         +(x(ie,2)-x(je,2))**2 )**0.5d0
         if (xlen.ne.0.d0) then
            temp1 = -cla*(x(ie,1)-x(je,1))/xlen*exp(-xlen/la)
     &              +clr*(x(ie,1)-x(je,1))/xlen*exp(-xlen/lr)
            temp2 = -cla*(x(ie,2)-x(je,2))/xlen*exp(-xlen/la)
     &              +clr*(x(ie,2)-x(je,2))/xlen*exp(-xlen/lr)
            force(3) = force(3) + temp1
            force(4) = force(4) + temp2

            force1 = force1 + temp1*xlen
            force2 = force2 + temp2*xlen
         endif
      enddo
      
      force(1)=x(ie,3)
      force(2)=x(ie,4)
      force(3)=( alpha*x(ie,3)
     &          -beta*x(ie,3)*vlen*vlen
     &          +force(3) )
      force(4)=( alpha*x(ie,4)
     &          -beta*x(ie,4)*vlen*vlen
     &          +force(4) )
      
      return
      end

c********************************************************************