SUBROUTINE LUCY
C     ===============
C
C     NLTE Lucy-Unsold temperature correction scheme, following
C       Werner & Dreizler
C
      INCLUDE 'IMPLIC.FOR'
      INCLUDE 'BASICS.FOR'
      INCLUDE 'ATOMIC.FOR'
      INCLUDE 'MODELQ.FOR'
      INCLUDE 'ODFPAR.FOR'
      INCLUDE 'ITERAT.FOR'
      INCLUDE 'ALIPAR.FOR'
      INCLUDE 'ARRAY1.FOR'
      DIMENSION HEAT(MDEPTH),ABSZ(MDEPTH),ABSP(MDEPTH),ABSH(MDEPTH),
     *    TOTJ(MDEPTH),TOTH(MDEPTH),TOTB(MDEPTH),DELH(MDEPTH),
     *    EDDF(MDEPTH),DELTAT(MDEPTH),dt1(mdepth),dt2(mdepth),
     *    ANTC(MDEPTH),XE(MDEPTH),heab(mdepth),tau(mdepth),
     *    TEM0(MDEPTH),TEM1(MDEPTH),TEM2(MDEPTH),TEM3(MDEPTH)
      DIMENSION COL(MTRANS),CLOC(MTRANS)
      PARAMETER (THIRD=1./3.)
      DATA ILINIT /0/
C
c
      IF(ILINIT.EQ.0) THEN
      ILINIT=1
      do itr=1,ntrans
         iluctr(itr)=0
      end do
      ntrl=0
      read(2,*,end=5,err=5) ntrl
    5 continue
      if(ntrl.gt.0) then
      do i=1,ntrl
         read(2,*,end=6,err=6) itrl
       iluctr(itrl)=1
      end do
      end if
    6 continue
      end if
c
c     IF(MOD(ILAM,4).NE.3) RETURN
      if(itlucy.le.0) return
      LAC2T=.FALSE.
      iacc0t=IACLT-3
C
      ilucy=1
   10 continue
      CALL OPAINI(0)
C
      DO ID=1,ND
        HEAT(ID)=0.
        heab(id)=0.
        ABSZ(ID)=0.
        ABSP(ID)=0.
        ABSH(ID)=0.
        TOTJ(ID)=0.
        TOTH(ID)=0.
        TOTB(ID)=0.
        EDDF(ID)=0.
      END DO
      EDDH=0.
C
      DO IJ=1,NFREQ
        W0=W(IJ)
        FR=FREQ(IJ)
        FR15=FR*1.D-15
        BNU=BN*FR15*FR15*FR15
      CALL OPACFL(IJ)
      CALL RTEFR1(IJ)
      tau(1)=abso1(1)/dens(1)*dm(1)
      DO ID=1,ND
        if(id.ge.2) tau(id)=tau(id-1)+half*(abso1(id)/dens(id)+
     *                abso1(id-1)/dens(id-1))*(dm(id)-dm(id-1))
          PLAND=BNU/(EXP(HK*FR/TEMP(ID))-UN)
        ABSOT0=ABSO1L(ID)-SCAT1(ID)
        HEAT(ID)=HEAT(ID)+W0*(ABSOT0*RAD1(ID)-EMIS1L(ID))
        HEAB(ID)=HEAB(ID)+W0*(ABSOT0*RAD1(ID)-EMIS1L(ID))/
     *             (ABSOT0*PLAND)
        ABSP(ID)=ABSP(ID)+W0*ABSOT0*PLAND
        TOTJ(ID)=TOTJ(ID)+W0*RAD1(ID)
        TOTB(ID)=TOTB(ID)+W0*PLAND
        EDDF(ID)=EDDF(ID)+W0*FAK1(ID)*RAD1(ID)
      END DO
C
      ID=1
      EDDH=EDDH+W0*RAD1(ID)*FH(IJ)
      TOTH(ID)=TOTH(ID)+W0*(RAD1(ID)*FH(IJ)-HEXTRD(IJ))
      ABSOT0=ABSO1(ID)
      ABSH(ID)=ABSH(ID)+W0*ABSOT0*(RAD1(ID)*FH(IJ)-HEXTRD(IJ))
      DO ID=2,ND
        ABSOT0=HALF*(ABSO1(ID)+ABSO1(ID-1))
        DTM=DELDM(ID-1)*(ABSO1(ID)*DENS1(ID)+ABSO1(ID-1)*DENS1(ID-1))
        FLUZ=(RAD1(ID)*FAK1(ID)-RAD1(ID-1)*FAK1(ID-1))/DTM
        TOTH(ID)=TOTH(ID)+W0*FLUZ
        ABSH(ID)=ABSH(ID)+W0*ABSOT0*FLUZ
      END DO
      END DO
C
      EDDH=EDDH/TOTJ(1)
      TEF4=SIG4P*TEFF**4
      TOTH(ND)=TEF4
      HEAT(ND)=0.
      HEAT(ND-1)=0.
      DO ID=1,ND
        ABSZ(ID)=(HEAT(ID)+ABSP(ID))/TOTJ(ID)
      ABSP(ID)=ABSP(ID)/TOTB(ID)
      ABSH(ID)=ABSH(ID)/TOTH(ID)
      EDDF(ID)=EDDF(ID)/TOTJ(ID)
      DELH(ID)=-TOTH(ID)+TEF4
      END DO
      DHHMX1=0.
      DO ID=1,ND-1
        DHH=ABS(DELH(ID)/TOTH(ID))
        IF(DHH.GT.DHHMX1) DHHMX1=DHH
      END DO
C
      ID=1
      TP3=TEMP1(ID)*TEMP1(ID)*TEMP1(ID)
      XX=EDDF(ID)/EDDH*DELH(ID)
      XX1=XX
      DT1(id)=HEAT(ID)/16./SIG4P*TP3/ABSP(ID)
      DT2(id)=ABSZ(ID)/EDDF(ID)*XX/16./SIG4P*TP3/ABSP(ID)
      DELTAT(ID)=dt1(id)+dt2(id)
      DO ID=2,ND
        TP3=TEMP1(ID)*TEMP1(ID)*TEMP1(ID)
        XX=XX+DELDM(ID)*(ABSH(ID)*DENS1(ID)*DELH(ID)+
     *                   ABSH(ID+1)*DENS1(ID+1)*DELH(ID+1))
        DT1(id)=HEAT(ID)/16./SIG4P*TP3/ABSP(ID)
        DT2(id)=ABSZ(ID)/EDDF(ID)*XX/16./SIG4P*TP3/ABSP(ID)
        DELTAT(ID)=dt1(id)+dt2(id)
      END DO
C
C     New temperature
C
      DO ID=1,ND
         TEMP(ID)=TEMP(ID)+DELTAT(ID)
         tem0(id)=temp(id)
         AOLD=DENS(ID)/WMM(ID)+ELEC(ID)
         XE(ID)=UN-ELEC(ID)/AOLD
      END DO
c
c     acceleration
c
      if(itlucy.lt.IACLT .or. ilucy.lt.iacc0t) go to 20
      ipng=1
      if(IACLDT.gt.0) ipng=mod((ilucy-IACLT),IACLDT)
      if(.not.lac2t) then
      IPT=MOD(ILUCY,3)
      IPT0=MOD(IACLT,3)
      IPT1=MOD((IACLT+1),3)
      IPT2=MOD((IACLT+2),3)
      IF(ILUCY.EQ.IACC0T) THEN
         DO ID=1,ND
            TEM3(ID)=TEM0(ID)
         END DO
       ELSE IF(IPT.EQ.IPT1) THEN
         DO ID=1,ND
            TEM2(ID)=TEM0(ID)
         END DO
       ELSE IF(IPT.EQ.IPT2) THEN
         DO ID=1,ND
            TEM1(ID)=TEM0(ID)
         END DO
      ENDIF
      else if (ipng.ne.0) then
         DO ID=1,ND
            TEM3(ID)=TEM2(ID)
         END DO
         DO ID=1,ND
            TEM2(ID)=TEM1(ID)
         END DO
         DO ID=1,ND
            TEM1(ID)=TEM0(ID)
         END DO
         GO TO 20
      end if

      IF(ILUCY.LT.IACLT) go to 20
C
      A1=0.
      B1=0.
      B2=0.
      C1=0.
      C2=0.
      DO ID=1,ND
       WT=0.
         IF(TEM0(ID).NE.0.) WT=1./ABS(TEM0(ID))
         D0=TEM0(ID)-TEM1(ID)
         D1=D0-TEM1(ID)+TEM2(ID)
         D2=D0-TEM2(ID)+TEM3(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
      AB=B2*A1-B1*B1
      IF(AB.EQ.0.) THEN
         WRITE(6,601) ILUCY,AB
         IACLT=IACLT+IACLDT
         IACC0T=IACLT-3
         GO TO 20
      ENDIF
      A0=(B2*C1-B1*C2)/AB
      B0=(A1*C2-B1*C1)/AB
      DO ID=1,ND
         TEM0(ID)=(1.-A0-B0)*TEM0(ID)+A0*TEM1(ID)+
     *                   B0*TEM2(ID)
         TEMP(ID)=TEM0(ID)
      END DO
      LAC2T=.TRUE.
  601 FORMAT(/,' **** ACCELT, ITER=',I4,'   AB = ',F7.3,/)
   20 CONTINUE

      CALL TDPINI
C
C     Integrate hydrostatic equilibrium      
C
      IF(IHECOR.GE.1) THEN
         ID=1
         PTUR=HALF*VTURB(ID)*VTURB(ID)*DENS(ID)
       ANTC(ID)=(DM(ID)*GRAV-PRD0-PTUR)/BOLK/TEMP(ID)
       IF(ANTC(ID).LE.0) ANTC(ID)=DENS(ID)/WMM(ID)+ELEC(ID)
         DO ID=2,ND
            PTUR=HALF*VTURB(ID)*VTURB(ID)*DENS(ID)
            PTURM=HALF*VTURB(ID-1)*VTURB(ID-1)*DENS(ID-1)
            ANTC(ID)=(GRAV*(DM(ID)-DM(ID-1))+
     *                BOLK*TEMP(ID-1)*ANTC(ID-1)-
     *                PRADT(ID)+PRADT(ID-1)-PTUR+PTURM)/
     *                BOLK/TEMP(ID)
         END DO
         DO ID=1,ND
            ELEC(ID)=(UN-XE(ID))*ANTC(ID)
            DENS(ID)=WMM(ID)*(ANTC(ID)-ELEC(ID))
          dens1(id)=un/dens(id)
         END DO
      END IF
C
C     Other depth-dependent quantities
C
      DO ID=1,ND
         PGS(ID)=(DENS(ID)/WMM(ID)+ELEC(ID))*BOLK*TEMP(ID)
         CALL WNSTOR(ID)
         IF(LTE) GO TO 60
         CALL SABOLF(ID)
         CALL COLIS(ID,TEMP(ID),COL,CLOC)
         DO I=1,NTRANS
            COLRAT(I,ID)=COL(I)
            COLTAR(I,ID)=CLOC(I)
         END DO
   60 CONTINUE
      END DO
C
C     new populations
C
      DO ID=1,ND
         CALL STEQEQ(ID,POP,1)
         IF(.NOT.LCHC.and.iter.lt.ielcor) CALL ELCOR(ID)
      END DO
C
      CALL CONCOR
      CALL ODFMER
      ilucy=ilucy+1
      if(ilucy.le.itlucy) go to 10
C
      RETURN
      END