      SUBROUTINE BBKICK(T1,IS,SS)
      IMPLICIT NONE
      INTEGER IS
      REAL*8 T1,SS(2)
      INCLUDE 'ctrlcm.h'
      INCLUDE 'beamcm.h'
      INCLUDE 'bbcom.h'
      INCLUDE 'cnstcm.h'
      INTEGER L,N,I,NDIV,IDIV,KIN
      REAL*8 XYS(3),XYS0(3),EP0(0:3),EP1(0:3),FXYS(3),FXYS0(3),
     %   DT,DT1,DS,SGN,EMASS,EMASS2,ERXY(2),ER,
     %   BBFLD(3,2),BBFLD0(3,2),BBFL1(3,2),EXTFXY(3),FF(3),WG
      INTEGER I2(3)/2,3,1/,I3(3)/3,1,2/
      REAL*8 ERMIN/1D-3/
C      ERMIN: maximum error in x and y in one time step in units
C             of mesh size.
C
      integer ncal/0/
      CALL CPUTIM('BBKICK',1)
      ncal=ncal+1
      DS=SS(2)-SS(1)
      EMASS=MASS(2)
      EMASS2=EMASS**2
      DO 400 N=1,NP
        IF(ISBIN(N).NE.IS) GOTO 400
        IF(LBBFIN(N).NE.0) GOTO 400
        IF(EP(3,N).GE.0) THEN
          L=1
          SGN=1
        ELSE
          L=2
          SGN=-1
        ENDIF
        IF(BBFLG(3-L).EQ.0) GOTO 400
        KIN=KIND(N)
        IF(KIN.EQ.1) THEN
          IF(LPAIR.LE.0) GOTO 400
C    Beam field on photon for coherent-pair calculation
          DT=T1-TXYS(0,N)
          DO 160 I=1,3
            XYS(I)=TXYS(I,N)+EP(I,N)/EP(0,N)*DT
 160      CONTINUE
          CALL BBKICK0(L,DS,1,TXYS(1,N),EP(0,N),BBFLD0,FF,NXY,BBQ)
          CALL BBKICK0(L,DS,1,XYS,EP(0,N),BBFLD,FF,NXY,BBQ)
          DO 170 I=1,3
            FLD(I,1,N)=FLD(I,1,N)+0.5D0*(BBFLD0(I,1)+BBFLD(I,1))
            FLD(I,2,N)=FLD(I,2,N)+0.5D0*(BBFLD0(I,2)+BBFLD(I,2))
 170      CONTINUE
          GOTO 390
        ENDIF
        DT=T1-TXYS(0,N)
        DO 180 I=1,3
          XYS0(I)=TXYS(I,N)
 180    CONTINUE
        DO 190 I=0,3
          EP0(I)=EP(I,N)
 190    CONTINUE
        DO 200 I=1,3
          EXTFXY(I)=CHARGE(KIN)*(FLD(I,1,N)
     %      +(EP0(I2(I))*FLD(I3(I),2,N)
     %       -EP0(I3(I))*FLD(I2(I),2,N))/EP0(0))
 200    CONTINUE
C-- Determine time step size
        CALL BBKICK0(L,DS,KIN,XYS0,EP0,BBFLD0,FF,NXY,BBQ)
        DO 210 I=1,3
          FXYS0(I)=EXTFXY(I)+FF(I)
          EP1(I)=EP0(I)+FXYS0(I)*DT
 210    CONTINUE
        EP1(0)=SQRT(EMASS2+EP1(1)**2+EP1(2)**2+EP1(3)**2)
        DO 215 I=1,3
          XYS(I)=XYS0(I)+(EP0(I)/EP0(0)+EP1(I)/EP1(0))/2*DT
 215    CONTINUE
        CALL BBKICK0(L,DS,KIN,XYS,EP1,BBFLD,FF,NXY,BBQ)
        DO 220 I=1,3
          FXYS(I)=EXTFXY(I)+FF(I)
 220    CONTINUE
        DO 225 I=1,2
          ERXY(I)=((FXYS(I)-FXYS0(I))/EP1(0)*DT**2/4)/BBDXY(I,L)
 225    CONTINUE
        ER=MAX(ABS(ERXY(1)),ABS(ERXY(2)))
C          Error in one step is proportional to 1/NDIV**3
C          Total error is to 1/NDIV**2
        NDIV=MAX(1,NINT(SQRT(ER/ERMIN)))
        IF(NDIV.EQ.1) THEN
C-- One step is enough
          DO 230 I=1,3
            BBFL1(I,1)=0.5D0*(BBFLD0(I,1)+BBFLD(I,1))
            BBFL1(I,2)=0.5D0*(BBFLD0(I,2)+BBFLD(I,2))
            EP1(I)=EP0(I)+(FXYS0(I)+FXYS(I))/2*DT
 230      CONTINUE
          EP1(0)=SQRT(EMASS2+EP1(1)**2+EP1(2)**2+EP1(3)**2)
          DO 235 I=1,3
            XYS(I)=XYS0(I)+(EP0(I)/EP0(0)+EP1(I)/EP1(0))/2*DT
 235      CONTINUE
        ELSE
C-- NDIV steps
          DT1=DT/NDIV
          DO 300 IDIV=1,NDIV
            WG=1
            IF(IDIV.EQ.NDIV) WG=0.5D0
            IF(NDIV.NE.1) THEN
              DO 240 I=0,3
                EP0(I)=EP1(I)
 240          CONTINUE
            ENDIF
            DO 245 I=1,2
              XYS(I)=XYS0(I)+EP0(I)/EP0(0)*DT1+FXYS0(I)/EP0(0)*DT1**2/2
 245        CONTINUE
            CALL BBKICK0(L,DS,KIN,XYS,EP0,BBFLD,FF,NXY,BBQ)
            DO 250 I=1,3
              FXYS(I)=EXTFXY(I)+FF(I)
              BBFL1(I,1)=BBFL1(I,1)+WG*BBFLD(I,1)
              BBFL1(I,2)=BBFL1(I,2)+WG*BBFLD(I,2)
 250        CONTINUE
            DO 260 I=1,3
              EP1(I)=EP0(I)+(FXYS0(I)+FXYS(I))/2*DT1
 260        CONTINUE
            EP1(0)=SQRT(EMASS2+EP1(1)**2+EP1(2)**2+EP1(3)**2)
            DO 265 I=1,3
              XYS(I)=XYS0(I)+(EP0(I)/EP0(0)+EP1(I)/EP1(0))/2*DT1
              XYS0(I)=XYS(I)
              FXYS0(I)=FXYS(I)
 265        CONTINUE
 300      CONTINUE
        ENDIF
        DO 320 I=1,3
          FLD(I,1,N)=FLD(I,1,N)+BBFL1(I,1)/NDIV
          FLD(I,2,N)=FLD(I,2,N)+BBFL1(I,2)/NDIV
 320    CONTINUE
        IF(ISPIN.GE.1) THEN
          CALL BMTEQ(EP(0,N),CHARGE(KIN),MASS(KIN),ANOM(KIN),
     %       FLD(1,1,N),DT,SPIN(1,N))
        ENDIF
C-- Store new variables. 
C    Extrapolate back to the initial time by a straight line.
C    Effects of the external field have already taken into account
C    so that only straight line extrapolation is enough at ENDPUSH
C    for particles with LBBFIN=1.
        DO 340 I=0,3
          EP(I,N)=EP1(I)
 340    CONTINUE
        DO 360 I=1,3
          TXYS(I,N)=XYS(I)-DT*EP(I,N)/EP(0,N)
 360    CONTINUE
C--
 390    LBBFIN(N)=1
 400  CONTINUE
      CALL CPUTIM('BBKICK',2)
      RETURN
      END


