SUBROUTINE ACCELP
C     =================
C
C     Acceleration of convergence for populations
C     (from Auer 1987, in Numerical Radiative Transfer p. 101)
C
      INCLUDE 'IMPLIC.FOR'
      INCLUDE 'BASICS.FOR'
      INCLUDE 'MODELQ.FOR'
      INCLUDE 'ITERAT.FOR'
      COMMON/POPULS/POPUL1(MLEVEL,MDEPTH),
     *              POPUL2(MLEVEL,MDEPTH),POPUL3(MLEVEL,MDEPTH)
C
      IF(NLAMBD.LT.IACPP.OR. ILAM.LT.IACC0P) RETURN
      ipng=1
      if(iacdp.gt.0) ipng=mod((ILAM-IACPP),IACDP)
      if(.not.lac2p) then
      IPT=MOD(ILAM,3)
      IPT0=MOD(IACPP,3)
      IPT1=MOD((IACPP+1),3)
      IPT2=MOD((IACPP+2),3)
      IF(ILAM.EQ.IACC0P) THEN
         DO ID=1,ND
            DO IX=1,NLEVEL
               POPUL3(IX,ID)=POPUL(IX,ID)
            END DO
         END DO
       ELSE IF(IPT.EQ.IPT1) THEN
         DO ID=1,ND
            DO IX=1,NLEVEL
               POPUL2(IX,ID)=POPUL(IX,ID)
            END DO
         END DO
       ELSE IF(IPT.EQ.IPT2) THEN
         DO ID=1,ND
            DO IX=1,NLEVEL
               POPUL1(IX,ID)=POPUL(IX,ID)
            END DO
         END DO
      END IF
      else if (ipng.ne.0) then
         DO ID=1,ND
            DO IX=1,NLEVEL
               POPUL3(IX,ID)=POPUL2(IX,ID)
            END DO
         END DO
         DO ID=1,ND
            DO IX=1,NLEVEL
               POPUL2(IX,ID)=POPUL1(IX,ID)
            END DO
         END DO
         DO ID=1,ND
            DO IX=1,NLEVEL
               POPUL1(IX,ID)=POPUL(IX,ID)
            END DO
         END DO
         RETURN
      end if

      IF(ILAM.LT.IACPP) RETURN
C
      A1=0.
      B1=0.
      B2=0.
      C1=0.
      C2=0.
      DO ID=1,ND
         DO IX=1,NLEVEL
            IF(POPUL(IX,ID).NE.0.) WT=1./ABS(POPUL(IX,ID))
            D0=POPUL(IX,ID)-POPUL1(IX,ID)
            D1=D0-POPUL1(IX,ID)+POPUL2(IX,ID)
            D2=D0-POPUL2(IX,ID)+POPUL3(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 DO
      AB=B2*A1-B1*B1
      IF(AB.EQ.0.) THEN
         WRITE(6,601) ILAM,AB
         WRITE(10,601) ILAM,AB
         IACPP=IACPP+IACDP
         IACC0P=IACPP-3
         RETURN
      ENDIF
      A=(B2*C1-B1*C2)/AB
      B=(A1*C2-B1*C1)/AB
C
      DO ID=1,ND
         DO IX=1,NLEVEL
            POPUL(IX,ID)=(1.-A-B)*POPUL(IX,ID)+A*POPUL1(IX,ID)+
     *                   B*POPUL2(IX,ID)
         END DO
      END DO
      WRITE(6,600) ILAM
      WRITE(10,600) ILAM
      LAC2P=.TRUE.
  600 FORMAT(' **** ACCELP, ITER=',I4)
  601 FORMAT(' **** ACCELP, ITER=',I4,'   AB = ',F7.3)
      RETURN
      END