SUBROUTINE ALIFR6(IJ)
C     =====================
C
C     hydrostatic and radiative equilibrium quantities - 
C     derivatives of the total heating and cooling rates in the 
C     ALI points with respect to the
C     temperature, electron density, and populations
C     a variant for consistent tridiagonal operator
C
      INCLUDE 'IMPLIC.FOR'     
      INCLUDE 'BASICS.FOR'     
      INCLUDE 'ATOMIC.FOR'     
      INCLUDE 'MODELQ.FOR'     
      INCLUDE 'ALIPAR.FOR'     
      PARAMETER (T23=TWO/3.D0, T43=4.D0/3.D0)
      DIMENSION DSFP1(MLVEXP),DSFP1M(MLVEXP),DSFP1D(MLVEXP),
     *          DSFP1P(MLVEXP),DSFPMM(MLVEXP)
C
      WW=WC(IJ)
C
      DSFT1M=0.
      DSFN1M=0.
      DO II=1,NLVEXP
         DSFP1M(II)=0.
      END DO
c
c ****1. Special expressions for the first depth  -  id=1
c
      ID=1
      
      LNSKIP=.NOT.LSKIP(ID,IJ)
c
c     Basic auxliliary quantities - derivatives of the source function
c
      EMISIV=UN/EMIS1(ID)
      if(ilmcor.ne.3) then
      if(ilasct.eq.0) then
        ABST=UN/(ABSO1(ID)-ELSCAT(ID))
        S0=EMIS1(ID)*ABST
        DSFN1=S0*(DEMN1(ID)*EMISIV-(DABN1(ID)-SIGEC(IJ))*ABST)
      else
        ABST=UN/ABSO1(ID)
        S0=EMIS1(ID)*ABST
        DSFN1=S0*(DEMN1(ID)*EMISIV-DABN1(ID)*ABST)
      end if
      DSFT1=S0*(DEMT1(ID)*EMISIV-DABT1(ID)*ABST)
      DO II=1,NLVEXP
        DSFP1(II)=S0*(DEMP1(II,ID)*EMISIV-DABP1(II,ID)*ABST)
      END DO
C
c     correction for electron scattering source function contribution
c
      else
        abst=un/abso1(id)
      s0=emis1(id)*abst
      sc=elec(id)*sigeC(IJ)
      sct=sc*abst
      st=s0+sct*rad1(id)
      corr=un/(un-ali1(id)*sct)
      dsft1=corr*(s0*demt1(id)*emisiv-st*dabt1(id)*abst)
      dsfn1=corr*(s0*demn1(id)*emisiv+sigeC(IJ)*rad1(id)*abst-
     *            st*dabn1(id)*abst)
        do ii=1,nlvexp
          dsfp1(ii)=corr*(s0*demp1(ii,id)*emisiv-
     *                   st*dabp1(ii,id)*abst)
        end do
      end if
 
      EMISIP=UN/EMIS1(ID+1)
      if(ilmcor.ne.3) then
      if(ilasct.eq.0) then
        ABSTP=UN/(ABSO1(ID+1)-ELSCAT(ID+1))
        S0P=EMIS1(ID+1)*ABSTP
        DSFN1P=S0P*(DEMN1(ID+1)*EMISIP-(DABN1(ID+1)-SIGEC(IJ))*ABSTP)
      else
        ABSTP=UN/ABSO1(ID+1)
        S0P=EMIS1(ID+1)*ABSTP
        DSFN1P=S0P*(DEMN1(ID+1)*EMISIP-DABN1(ID+1)*ABSTP)
      end if
      DSFT1P=S0P*(DEMT1(ID+1)*EMISIP-DABT1(ID+1)*ABSTP)
      DO II=1,NLVEXP
        DSFP1P(II)=S0P*(DEMP1(II,ID+1)*EMISIP-DABP1(II,ID+1)*ABSTP)
      END DO
      else
      abstp=un/abso1(id+1)
      s0p=emis1(id+1)*abstp
      scp=elec(id+1)*sigeC(IJ)
      sctp=scp*abstp
      stp=s0p+sctp*rad1(id+1)
      corrp=un/(un-ali1(id+1)*sctp)
      dsft1p=corrp*(s0p*demt1(id+1)*emisip-stp*dabt1(id+1)*abstp)
      dsfn1p=corrp*(s0p*demn1(id+1)*emisip+sigeC(IJ)*rad1(id+1)*abstp-
     *            stp*dabn1(id+1)*abstp)
        do ii=1,nlvexp
          dsfp1p(ii)=corrp*(s0p*demp1(ii,id+1)*emisip-
     *                    stp*dabp1(ii,id+1)*abstp)
        end do
      end if
      
      IF(IRDER.EQ.1.OR.IRDER.EQ.3) THEN
        DSFDT(ID)=DSFT1*ALI1(ID)
        DSFDN(ID)=DSFN1*ALI1(ID)
      END IF
      IF(IRDER.GT.1) THEN
        DO II=1,NLVEXP
          DSFDP(II,ID)=DSFP1(II)*ALI1(ID)
        END DO
      END IF
 
      IF(IRDER.EQ.1.OR.IRDER.EQ.3) THEN
        DSFDTM(ID)=DSFT1M*ALIM1(ID)
        DSFDNM(ID)=DSFN1M*ALIM1(ID)
        DSFDTP(ID)=DSFT1P*ALIP1(ID)
        DSFDNP(ID)=DSFN1P*ALIP1(ID)
      END IF
      IF(IRDER.GT.1) THEN
        DO II=1,NLVEXP
          DSFDPM(II,ID)=DSFP1M(II)*ALIM1(ID)
          DSFDPP(II,ID)=DSFP1P(II)*ALIP1(ID)
        END DO
      END IF      
c
c     Hydrostatic equilibrium quantities
c
      WF=WW*FH(IJ)
      IF(LNSKIP) THEN
        FPRD(ID)=FPRD(ID)+WF*ABSO1(ID)*RAD1(ID)-
     *               WW*HEXTRD(IJ)*ABSO1(ID)
        E0=WF*RAD1(ID)
        D0=WF*ABSO1(ID)*ALI1(ID)
        HEIT(ID)=HEIT(ID)+D0*DSFT1+E0*DABT1(ID)
        HEIN(ID)=HEIN(ID)+D0*DSFN1+E0*DABN1(ID)
        DO II=1,NLVEXP
          HEIP(II,ID)=HEIP(II,ID)+D0*DSFP1(II)+E0*DABP1(II,ID)
        END DO
        IF(IFALI.GE.7) THEN
          D0P=WF*ABSO1(ID)*ALIP1(ID)
          HEITP(ID)=HEITP(ID)+D0P*DSFT1P
          HEINP(ID)=HEINP(ID)+D0P*DSFN1P
          DO II=1,NLVEXP
            HEIPP(II,ID)=HEIPP(II,ID)+D0P*DSFP1P(II)
          END DO
        END IF
      END IF
c
c     Differential equation part of radiative equilibrium
c
      FLFIX(ID)=FLFIX(ID)+WF*RAD1(ID)-WW*HEXTRD(IJ)
      IF(REDIF(ID).GT.0.) THEN
        WF=WF*ALI1(ID)
        REDT(ID)=REDT(ID)+WF*DSFT1
        REDN(ID)=REDN(ID)+WF*DSFN1
        DO II=1,NLVEXP
          REDP(II,ID)=REDP(II,ID)+WF*DSFP1(II)
        END DO
      END IF
c
C     Integral equation part of the radiative equilibrium 
C
      IF(REINT(ID).GT.0) THEN
      if(ilmcor.ne.3) then
        if(ilasct.eq.0) then
          ABST=ABSO1(ID)-ELSCAT(ID)
          WWK=WW*ABST
          FCOOLI(ID)=FCOOLI(ID)+WW*(EMIS1(ID)-ABST*RAD1(ID))
          D0=WW*(ALI1(ID)-UN)*ABST
          E0=WW*(RAD1(ID)-S0)
          REIN(ID)=REIN(ID)+D0*DSFN1+E0*(DABN1(ID)-SIGEC(IJ))
        else
          SRH=SIGE*DENS1(ID)
          ABST=ABSO1(ID)
          ABSTE=ABST-ELSCAT(ID)
          WWK=WW*ABST
          FCOOLI(ID)=FCOOLI(ID)+WW*(EMIS1(ID)-ABSTE*RAD1(ID))
          D0=WW*(ABSTE*ALI1(ID)-ABST)
          E0=WW*(RAD1(ID)-S0)
          REIN(ID)=REIN(ID)+D0*DSFN1+E0*DABN1(ID)-WW*SIGEC(IJ)*RAD1(ID)
        end if
        REIT(ID)=REIT(ID)+D0*DSFT1+E0*DABT1(ID)
        DO II=1,NLVEXP
          REIP(II,ID)=REIP(II,ID)+D0*DSFP1(II)+E0*DABP1(II,ID)
        END DO
      else
          abst=abso1(id)-elscat(id)
        d0=abst*ali1(id)
        fcooli(id)=fcooli(id)+ww*(emis1(id)-abst*rad1(id))
          rein(id)=rein(id)+ww*(d0*dsfn1+
     *                   rad1(id)*(dabn1(id)-sigeC(IJ))-demn1(id))
          do ii=1,nlvexp
            reip(ii,id)=reip(ii,id)+ww*(d0*dsfp1(ii)+
     *                     rad1(id)*dabp1(ii,id)-demp1(ii,id))
          end do
           reit(id)=reit(id)+ww*(d0*dsft1+
     *                  rad1(id)*dabt1(id)-demt1(id))
        end if
c
c       the following are needed only for tridiagonal Lambda^star
C
C       Upper sub-diagonal band
C
        WWKC=WWK*ALIP1(ID)
        CREIT(ID)=CREIT(ID)+WWKC*DSFT1P
        CREIN(ID)=CREIN(ID)+WWKC*DSFN1P
        DO II=1,NLVEXP
          CREIP(II,ID)=CREIP(II,ID)+WWKC*DSFP1P(II)
        END DO
      END IF
      
C
c ****2. loop over depths
c
      DO ID=2,ND-1
      
        LNSKIP=.NOT.LSKIP(ID,IJ)
      
        DSFTMM=DSFT1M
        DSFNMM=DSFN1M
        DO II=1,NLVEXP
          DSFPMM(II)=DSFP1M(II)
        END DO
        DSFT1M=DSFT1
        DSFN1M=DSFN1
        DO II=1,NLVEXP
          DSFP1M(II)=DSFP1(II)
        END DO
        S0=S0P
        DSFT1=DSFT1P
        DSFN1=DSFN1P
        DO II=1,NLVEXP
          DSFP1(II)=DSFP1P(II)
        END DO

        EMISIP=UN/EMIS1(ID+1)
        if(ilmcor.ne.3) then
        if(ilasct.eq.0) then
          ABSTP=UN/(ABSO1(ID+1)-ELSCAT(ID+1))
          S0P=EMIS1(ID+1)*ABSTP
          DSFN1P=S0P*(DEMN1(ID+1)*EMISIP-(DABN1(ID+1)-SIGEC(IJ))*ABSTP)
        else
          ABSTP=UN/ABSO1(ID+1)
          S0P=EMIS1(ID+1)*ABSTP
          DSFN1P=S0P*(DEMN1(ID+1)*EMISIP-DABN1(ID+1)*ABSTP)
        end if
        DSFT1P=S0P*(DEMT1(ID+1)*EMISIP-DABT1(ID+1)*ABSTP)
        DO II=1,NLVEXP
          DSFP1P(II)=S0P*(DEMP1(II,ID+1)*EMISIP-DABP1(II,ID+1)*ABSTP)
        END DO
        else
      abstp=un/abso1(id+1)
      s0p=emis1(id+1)*abstp
      scp=elec(id+1)*sigeC(IJ)
      sctp=scp*abstp
      stp=s0p+sctp*rad1(id+1)
      corrp=un/(un-ali1(id+1)*sctp)
      dsft1p=corrp*(s0p*demt1(id+1)*emisip-stp*dabt1(id+1)*abstp)
      dsfn1p=corrp*(s0p*demn1(id+1)*emisip+sigeC(IJ)*rad1(id+1)*abstp-
     *            stp*dabn1(id+1)*abstp)
        do ii=1,nlvexp
          dsfp1p(ii)=corrp*(s0p*demp1(ii,id+1)*emisip-
     *                    stp*dabp1(ii,id+1)*abstp)
        end do
        end if
      
        IF(IRDER.EQ.1.OR.IRDER.EQ.3) THEN
          DSFDT(ID)=DSFT1*ALI1(ID)
          DSFDN(ID)=DSFN1*ALI1(ID)
        END IF
        IF(IRDER.GT.1) THEN
          DO II=1,NLVEXP
            DSFDP(II,ID)=DSFP1(II)*ALI1(ID)
         END DO
        END IF

        IF(IRDER.EQ.1.OR.IRDER.EQ.3) THEN
          DSFDTM(ID)=DSFT1M*ALIM1(ID)
          DSFDNM(ID)=DSFN1M*ALIM1(ID)
          DSFDTP(ID)=DSFT1P*ALIP1(ID)
          DSFDNP(ID)=DSFN1P*ALIP1(ID)
        END IF
        IF(IRDER.GT.1) THEN
          DO II=1,NLVEXP
            DSFDPM(II,ID)=DSFP1M(II)*ALIM1(ID)
            DSFDPP(II,ID)=DSFP1P(II)*ALIP1(ID)
          END DO
        END IF
c
c       Hydrostatic equilibrium equation
c
        IF(LNSKIP) THEN
          D0=WW*FAK1(ID)
          A0=WW*FAK1(ID-1)
          FPRD(ID)=FPRD(ID)+D0*RAD1(ID)-A0*RAD1(ID-1)
          F0=D0*ALIP1(ID)
          E0=D0*ALIM1(ID)-A0*ALI1(ID-1)
          D0=D0*ALI1(ID)-A0*ALIP1(ID-1)
          HEIT(ID)=HEIT(ID)+D0*DSFT1
          HEIN(ID)=HEIN(ID)+D0*DSFN1
          HEITM(ID)=HEITM(ID)+E0*DSFT1M
          HEINM(ID)=HEINM(ID)+E0*DSFN1M
          DO II=1,NLVEXP
            HEIP(II,ID)=HEIP(II,ID)+D0*DSFP1(II)
            HEIPM(II,ID)=HEIPM(II,ID)+E0*DSFP1M(II)
          END DO
          IF(IFALI.GE.7) THEN
            HEITP(ID)=HEITP(ID)+F0*DSFT1P
            HEINP(ID)=HEINP(ID)+F0*DSFN1P
            DO II=1,NLVEXP
              HEIPP(II,ID)=HEIPP(II,ID)+F0*DSFP1P(II)
            END DO
            A0M=A0*ALIM1(ID-1)
            EHET(ID)=EHET(ID)-A0M*DSFTMM
            EHEN(ID)=EHEN(ID)-A0M*DSFNMM
            DO II=1,NLVEXP
              EHEP(II,ID)=EHEP(II,ID)-A0M*DSFPMM(II)
            END DO
          END IF
        END IF
C   
C       Differential equation part of radiative equilibrium
C
        DDT=UN/(ABSOT(ID)+ABSOT(ID-1))
        DT=DDT/DELDMZ(ID-1)
        FL=(RAD1(ID)*FAK1(ID)-RAD1(ID-1)*FAK1(ID-1))*DT
        FLFIX(ID)=FLFIX(ID)+WW*FL
        IF(REDIF(ID).GT.0) THEN
          D0=WW*FAK1(ID)*DT
          A0=WW*FAK1(ID-1)*DT
          D0M=D0*ALIM1(ID)-A0*ALI1(ID-1)
          D0P=D0*ALIP1(ID)
          D0=D0*ALI1(ID)-A0*ALIP1(ID-1)
          E0=WW*FL*DDT
          REDX(ID)=REDX(ID)+E0*ABSO1(ID)
          REDXM(ID)=REDXM(ID)+E0*ABSO1(ID-1)
          E0M=E0*DENSI(ID-1)
          E0=E0*DENSI(ID)
          REDT(ID)=REDT(ID)+D0*DSFT1-E0*DABT1(ID)
          REDTM(ID)=REDTM(ID)+D0M*DSFT1M-E0M*DABT1(ID-1)
          REDN(ID)=REDN(ID)+D0*DSFN1-E0*DABN1(ID)
          REDNM(ID)=REDNM(ID)+D0M*DSFN1M-E0M*DABN1(ID-1)
          DO II=1,NLVEXP
            REDP(II,ID)=REDP(II,ID)+D0*DSFP1(II)-E0*DABP1(II,ID)
            REDPM(II,ID)=REDPM(II,ID)+D0M*DSFP1M(II)-
     *                         E0M*DABP1(II,ID-1)
          END DO
          IF(IFALI.GE.7) THEN
            REDTP(ID)=REDTP(ID)+D0P*DSFT1M
            REDNP(ID)=REDNP(ID)+D0P*DSFN1M
            DO II=1,NLVEXP
              REDPP(II,ID)=REDPP(II,ID)+D0P*DSFP1P(II)
            END DO
            A0M=A0*ALIM1(ID-1)
            ERET(ID)=ERET(ID)-A0M*DSFTMM
            EREN(ID)=EREN(ID)-A0M*DSFNMM
            DO II=1,NLVEXP
              EREP(II,ID)=EREP(II,ID)-A0M*DSFPMM(II)
            END DO
          END IF
        END IF   
c
C       Integral equation part of the radiative equilibrium 
C
        IF(REINT(ID).GT.0) THEN
        if(ilmcor.ne.3) then
          if(ilasct.eq.0) then
            ABST=ABSO1(ID)-ELSCAT(ID)
            WWK=WW*ABST
            FCOOLI(ID)=FCOOLI(ID)+WW*(EMIS1(ID)-ABST*RAD1(ID))
            D0=WW*(ALI1(ID)-UN)*ABST
            E0=WW*(RAD1(ID)-S0)
            REIN(ID)=REIN(ID)+D0*DSFN1+E0*(DABN1(ID)-SIGEC(IJ))
          else
            SRH=SIGE*DENS1(ID)
            ABST=ABSO1(ID)
            ABSTE=ABST-ELSCAT(ID)
            WWK=WW*ABST
            FCOOLI(ID)=FCOOLI(ID)+WW*(EMIS1(ID)-ABSTE*RAD1(ID))
            D0=WW*(ABSTE*ALI1(ID)-ABST)
            E0=WW*(RAD1(ID)-S0)
            REIN(ID)=REIN(ID)+D0*DSFN1+E0*DABN1(ID)-
     *                  WW*SIGE*RAD1(ID)
          end if
          REIT(ID)=REIT(ID)+D0*DSFT1+E0*DABT1(ID)
          DO II=1,NLVEXP
            REIP(II,ID)=REIP(II,ID)+D0*DSFP1(II)+E0*DABP1(II,ID)
          END DO
        else
          abst=abso1(id)-elscat(id)
        d0=abst*ali1(id)
        fcooli(id)=fcooli(id)+ww*(emis1(id)-abst*rad1(id))
          rein(id)=rein(id)+ww*(d0*dsfn1+
     *                   rad1(id)*(dabn1(id)-sige)-demn1(id))
          do ii=1,nlvexp
            reip(ii,id)=reip(ii,id)+ww*(d0*dsfp1(ii)+
     *                     rad1(id)*dabp1(ii,id)-demp1(ii,id))
          end do
           reit(id)=reit(id)+ww*(d0*dsft1+
     *                  rad1(id)*dabt1(id)-demt1(id))
          end if
c  
c         the following are needed only for tridiagonal Lambda^star
C
C         Lower sub-diagonal band
C
          WWKA=WWK*ALIM1(ID)
          AREIT(ID)=AREIT(ID)+WWKA*DSFT1M
          AREIN(ID)=AREIN(ID)+WWKA*DSFN1M
          DO II=1,NLVEXP
            AREIP(II,ID)=AREIP(II,ID)+WWKA*DSFP1M(II)
          END DO
C
C         Upper sub-diagonal band
C
          WWKC=WWK*ALIP1(ID)
          CREIT(ID)=CREIT(ID)+WWKC*DSFT1P
          CREIN(ID)=CREIN(ID)+WWKC*DSFN1P
          DO II=1,NLVEXP
            CREIP(II,ID)=CREIP(II,ID)+WWKC*DSFP1P(II)
          END DO
        END IF

      END DO    
  
C
c ****3. deepest point - ID=ND
c
      ID=ND
      
      LNSKIP=.NOT.LSKIP(ID,IJ)
      
      DSFTMM=DSFT1M
      DSFNMM=DSFN1M
      DO II=1,NLVEXP
        DSFPMM(II)=DSFP1M(II)
      END DO
      DSFT1M=DSFT1
      DSFN1M=DSFN1
      DO II=1,NLVEXP
        DSFP1M(II)=DSFP1(II)
      END DO
      S0=S0P
      DSFT1=DSFT1P
      DSFN1=DSFN1P
      DO II=1,NLVEXP
        DSFP1(II)=DSFP1P(II)
      END DO
C
C     Improved lower boundary condition
C
      IF(IBC.GT.0.AND.IDISK.EQ.0) THEN
        DT=UN/(DELDMZ(ID-1)*(ABSOT(ID)+ABSOT(ID-1)))
        PLAD=XKFB(ID)/XKF1(ID)
        DBDT=PLAD/XKF1(ID)*HKT21(ID)*FREQ(IJ)*DT
        IF(IBC.EQ.1) THEN
          DSFT1=DSFT1+DBDT
        ELSE IF(IBC.GE.2) THEN
          PLAM=XKFB(ID-1)/XKF1(ID-1)
          TAU23=T23*DT
          TAU43=T43*DT
          D0=(PLAD*(UN+TAU43)-T43*PLAM*DT)*DT*DT
          RHD=DELDMZ(ID-1)*DENSI(ID)
          E0=D0*RHD
          DSFT1=DSFT1+DBDT*(UN+TAU23)-E0*DABT1(ID)
          DSFN1=DSFN1-E0*(DABN1(ID)+ABSO1(ID)*DENSIM(ID))
          DO II=1,NLVEXP
            DSFP1(II)=DSFP1(II)-E0*DABP1(II,ID)
          END DO
          IF(IBC.GE.3) THEN
            DBDTM=PLAM/XKF1(ID-1)*HKT21(ID-1)*FREQ(IJ)*DT
            RHD=DELDMZ(ID-1)*DENSI(ID-1)
            E0=D0*RHD
            DSFT1D=-DBDTM*DT*T23-E0*DABT1(ID-1)
            DSFN1D=-E0*(DABN1(ID-1)+ABSO1(ID-1)*DENSIM(ID-1))
            DO II=1,NLVEXP
              DSFP1D(II)=-E0*DABP1(II,ID-1)
            END DO
          END IF
        END IF
      END IF
C
      IF(IRDER.EQ.1.OR.IRDER.EQ.3) THEN
        DSFDT(ID)=DSFT1*ALI1(ID)
        DSFDN(ID)=DSFN1*ALI1(ID)
      END IF
      IF(IRDER.GT.1) THEN
        DO II=1,NLVEXP
          DSFDP(II,ID)=DSFP1(II)*ALI1(ID)
        END DO
      END IF

      IF(IRDER.EQ.1.OR.IRDER.EQ.3) THEN
        DSFDTM(ID)=DSFT1M*ALIM1(ID)
        DSFDNM(ID)=DSFN1M*ALIM1(ID)
        DSFDTP(ID)=DSFT1P*ALIP1(ID)
        DSFDNP(ID)=DSFN1P*ALIP1(ID)
      END IF
      IF(IRDER.GT.1) THEN
        DO II=1,NLVEXP
          DSFDPM(II,ID)=DSFP1M(II)*ALIM1(ID)
          DSFDPP(II,ID)=DSFP1P(II)*ALIP1(ID)
        END DO
      END IF
c
c     Hydrostatic equilibrium equation
c
      IF(LNSKIP) THEN
        D0=WW*FAK1(ID)
        A0=WW*FAK1(ID-1)
        FPRD(ID)=FPRD(ID)+D0*RAD1(ID)-A0*RAD1(ID-1)
        F0=D0*ALIP1(ID)
        E0=D0*ALIM1(ID)-A0*ALI1(ID-1)
        D0=D0*ALI1(ID)-A0*ALIP1(ID-1)
        HEIT(ID)=HEIT(ID)+D0*DSFT1
        HEIN(ID)=HEIN(ID)+D0*DSFN1
        HEITM(ID)=HEITM(ID)+E0*DSFT1M
        HEINM(ID)=HEINM(ID)+E0*DSFN1M
        DO II=1,NLVEXP
          HEIP(II,ID)=HEIP(II,ID)+D0*DSFP1(II)
          HEIPM(II,ID)=HEIPM(II,ID)+E0*DSFP1M(II)
        END DO
        IF(IFALI.GE.7) THEN
          HEITP(ID)=HEITP(ID)+F0*DSFT1P
          HEINP(ID)=HEINP(ID)+F0*DSFN1P
          DO II=1,NLVEXP
            HEIPP(II,ID)=HEIPP(II,ID)+F0*DSFP1P(II)
          END DO
          A0M=A0*ALIM1(ID-1)
          EHET(ID)=EHET(ID)-A0M*DSFTMM
          EHEN(ID)=EHEN(ID)-A0M*DSFNMM
          DO II=1,NLVEXP
            EHEP(II,ID)=EHEP(II,ID)-A0M*DSFPMM(II)
          END DO
        END IF
        IF(IBC.GE.3) THEN
          HEITM(ID)=HEITM(ID)-D0*DSFT1D
          HEINM(ID)=HEINM(ID)-D0*DSFN1D
          DO II=1,NLVEXP
            HEIPM(II,ID)=HEIPM(II,ID)-D0*DSFP1D(II)
          END DO
        END IF   
      END IF
C   
C     Differential equation part of radiative equilibrium
C
      DDT=UN/(ABSOT(ID)+ABSOT(ID-1))
      DT=DDT/DELDMZ(ID-1)
      FL=(RAD1(ID)*FAK1(ID)-RAD1(ID-1)*FAK1(ID-1))*DT
      FLFIX(ID)=FLFIX(ID)+WW*FL
      IF(REDIF(ID).GT.0) THEN
        D0=WW*FAK1(ID)*DT
        A0=WW*FAK1(ID-1)*DT
        D0M=D0*ALIM1(ID)-A0*ALI1(ID-1)
        D0P=D0*ALIP1(ID)
        D0=D0*ALI1(ID)-A0*ALIP1(ID-1)
        E0=WW*FL*DDT
        REDX(ID)=REDX(ID)+E0*ABSO1(ID)
        REDXM(ID)=REDXM(ID)+E0*ABSO1(ID-1)
        E0M=E0*DENSI(ID-1)
        E0=E0*DENSI(ID)
        REDT(ID)=REDT(ID)+D0*DSFT1-E0*DABT1(ID)
        REDTM(ID)=REDTM(ID)+D0M*DSFT1M-E0M*DABT1(ID-1)
        REDN(ID)=REDN(ID)+D0*DSFN1-E0*DABN1(ID)
        REDNM(ID)=REDNM(ID)+D0M*DSFN1M-E0M*DABN1(ID-1)
        DO II=1,NLVEXP
          REDP(II,ID)=REDP(II,ID)+D0*DSFP1(II)-E0*DABP1(II,ID)
          REDPM(II,ID)=REDPM(II,ID)+D0M*DSFP1M(II)-
     *                         E0M*DABP1(II,ID-1)
        END DO
        IF(IFALI.GE.7) THEN
          REDTP(ID)=REDTP(ID)+D0P*DSFT1M
          REDNP(ID)=REDNP(ID)+D0P*DSFN1M
          DO II=1,NLVEXP
             REDPP(II,ID)=REDPP(II,ID)+D0P*DSFP1P(II)
          END DO
          A0M=A0*ALIM1(ID-1)
          ERET(ID)=ERET(ID)-A0M*DSFTMM
          EREN(ID)=EREN(ID)-A0M*DSFNMM
          DO II=1,NLVEXP
            EREP(II,ID)=EREP(II,ID)-A0M*DSFPMM(II)
          END DO
        END IF
        IF(IBC.GE.3) THEN
          REDTM(ID)=REDTM(ID)+D0*DSFT1D
          REDNM(ID)=REDNM(ID)+D0*DSFN1D
          DO II=1,NLVEXP
            REDPM(II,ID)=REDPM(II,ID)+D0*DSFP1D(II)
          END DO
        END IF   
      END IF   
c
C     Integral equation part of the radiative equilibrium 
C
      IF(REINT(ID).GT.0) THEN
      if(ilmcor.ne.3) then
        if(ilasct.eq.0) then
          ABST=ABSO1(ID)-ELSCAT(ID)
          WWK=WW*ABST
          FCOOLI(ID)=FCOOLI(ID)+WW*(EMIS1(ID)-ABST*RAD1(ID))
          D0=WW*(ALI1(ID)-UN)*ABST
          E0=WW*(RAD1(ID)-S0)
          REIN(ID)=REIN(ID)+D0*DSFN1+E0*(DABN1(ID)-SIGEC(IJ))
        else
          SRH=SIGE*DENS1(ID)
          ABST=ABSO1(ID)
          ABSTE=ABST-ELSCAT(ID)
          WWK=WW*ABST
          FCOOLI(ID)=FCOOLI(ID)+WW*(EMIS1(ID)-ABSTE*RAD1(ID))
          D0=WW*(ABSTE*ALI1(ID)-ABST)
          E0=WW*(RAD1(ID)-S0)
          REIN(ID)=REIN(ID)+D0*DSFN1+E0*DABN1(ID)-WW*SIGEC(IJ)*RAD1(ID)
        end if
        IF(IBC.EQ.0) THEN
          REIT(ID)=REIT(ID)+D0*DSFT1+E0*DABT1(ID)
        ELSE
          REIT(ID)=REIT(ID)+D0*(DSFT1-DBDT)+E0*DABT1(ID)+
     *                  ALI1(ID)/ABST*DBDT
        END IF
        DO II=1,NLVEXP
          REIP(II,ID)=REIP(II,ID)+D0*DSFP1(II)+E0*DABP1(II,ID)
        END DO
      else
          abst=abso1(id)-elscat(id)
        d0=abst*ali1(id)
        fcooli(id)=fcooli(id)+ww*(emis1(id)-abst*rad1(id))
          rein(id)=rein(id)+ww*(d0*dsfn1+
     *                   rad1(id)*(dabn1(id)-sigeC(IJ))-demn1(id))
          do ii=1,nlvexp
            reip(ii,id)=reip(ii,id)+ww*(d0*dsfp1(ii)+
     *                     rad1(id)*dabp1(ii,id)-demp1(ii,id))
          end do
          if(ibc.eq.0) then
            reit(id)=reit(id)+ww*(d0*dsft1+
     *                  rad1(id)*dabt1(id)-demt1(id))
          else
            reit(id)=reit(id)+ww*(d0*(dsft1-dbdt)+e0*dabt1(id)+
     *                  rad1(id)*dabt1(id)-demt1(id)+
     *                  ali1(id)/abst*dbdt)
          end if
        end if
c
c       the following are needed only for tridiagonal Lambda^star
C
C       Lower sub-diagonal band
C
        WWKA=WWK*ALIM1(ID)
        AREIT(ID)=AREIT(ID)+WWKA*DSFT1M
        AREIN(ID)=AREIN(ID)+WWKA*DSFN1M
        DO II=1,NLVEXP
          AREIP(II,ID)=AREIP(II,ID)+WWKA*DSFP1M(II)
        END DO
      END IF

      RETURN
      END