SUBROUTINE ACCEL2
C     =================
C
C     Acceleration of convergence (from Auer 1987, in Numerical
C     Radiative Transfer p. 101)
C
      INCLUDE 'IMPLIC.FOR'
      INCLUDE 'BASICS.FOR'
      INCLUDE 'ITERAT.FOR'
      INCLUDE 'MODELQ.FOR'
C
      IF(NITER.LT.IACC .OR. ITER.LT.IACC0) RETURN
      ipng=1
      if(iacd.gt.0) ipng=mod((iter-iacc),iacd)
      if(.not.lac2) then
      IPT=MOD(ITER,3)
      IPT0=MOD(IACC,3)
      IPT1=MOD((IACC+1),3)
      IPT2=MOD((IACC+2),3)
      IF(ITER.EQ.IACC0) THEN
         DO ID=1,ND
            DO IX=1,NN
               PSY3(IX,ID)=PSY0(IX,ID)
            END DO
         END DO
       ELSE IF(IPT.EQ.IPT1) THEN
         DO ID=1,ND
            DO IX=1,NN
               PSY2(IX,ID)=PSY0(IX,ID)
            END DO
         END DO
       ELSE IF(IPT.EQ.IPT2) THEN
         DO ID=1,ND
            DO IX=1,NN
               PSY1(IX,ID)=PSY0(IX,ID)
            END DO
         END DO
      END IF
      else if (ipng.ne.0) then
         DO ID=1,ND
            DO IX=1,NN
               PSY3(IX,ID)=PSY2(IX,ID)
            END DO
         END DO
         DO ID=1,ND
            DO IX=1,NN
               PSY2(IX,ID)=PSY1(IX,ID)
            END DO
         END DO
         DO ID=1,ND
            DO IX=1,NN
               PSY1(IX,ID)=PSY0(IX,ID)
            END DO
         END DO
         RETURN
      end if

      IF(ITER.LT.IACC) RETURN
C
      A1=0.
      B1=0.
      B2=0.
      C1=0.
      C2=0.
      DO IX=1,NN
         IF(LSNG(IX)) THEN
         DO ID=1,ND
            WT=0.
            IF(PSY0(IX,ID).NE.0.) WT=1./ABS(PSY0(IX,ID))
            D0=PSY0(IX,ID)-PSY1(IX,ID)
            D1=D0-PSY1(IX,ID)+PSY2(IX,ID)
            D2=D0-PSY2(IX,ID)+PSY3(IX,ID)
            A1=A1+WT*D1*D1
            B1=B1+WT*D1*D2
            B2=B2+WT*D2*D2
            C1=C1+WT*D0*D1
            C2=C2+WT*D0*D2
         END DO
         END IF
      END DO
      AB=B2*A1-B1*B1
      IF(AB.EQ.0.) THEN
         WRITE(6,601) ITER,AB
         WRITE(10,601) ITER,AB
         IACC=IACC+IACD
         IACC0=IACC-3
         RETURN
      ENDIF
      A=(B2*C1-B1*C2)/AB
      B=(A1*C2-B1*C1)/AB
C
      DO ID=1,ND
         DO IX=1,NN
            PSY0(IX,ID)=(1.-A-B)*PSY0(IX,ID)+A*PSY1(IX,ID)+
     *                   B*PSY2(IX,ID)
         END DO
      END DO
      WRITE(6,600) ITER
      WRITE(10,600) ITER
      LAC2=.TRUE.
      LRES2=.FALSE.
c
c     call RESOLV after evaluating the accelerated estimate
c
      CALL RESOLV
      LRES2=.TRUE.
      RETURN
  600 FORMAT(' **** ACCEL2, ITER=',I4)
  601 FORMAT(' **** ACCEL2, ITER=',I4,'   AB = ',F7.3)
      END