SUBROUTINE ODF1(IMODE,IL,IU,ID,ODF)
C     ===================================
C
C     opacity distribution function for overlapping lines near the series limit
C
C     The lines converge to the edge of the (continuum) transition IL - IU,
C     IL - index of the lower level
C     IU - index of the upper level (usually the ground state of teh next ion)
C     ID - depth index
C
C     Output: ODF - opacity distribution function interpolated to the set of
C                   explicit frequencies
C
      INCLUDE 'IMPLIC.FOR'
      INCLUDE 'BASICS.FOR'
      INCLUDE 'ATOMIC.FOR'
      INCLUDE 'MODELQ.FOR'
      INCLUDE 'ODFPAR.FOR'
      PARAMETER (FRH=3.28805D15,CQT=1.284523D12)
      PARAMETER (CCOR=0.09,C00=1.25D-9,CID=0.02654,SIXTH=UN/6.)
      DIMENSION FRO(MFRO),ODF0(MFRO),ABSO0(MFRO),ODF(MFREQ),SGT(MFRO),
     *          ALAM(MFRO),
     *          FROD(MFRO),SGFR(MFRO),IODF(MFRO),IODR(MFRO),DWF(MFRO)
      SAVE FRO,SGFR,ODF0,IODF
C
      I=NQUANT(IL)
      KL=INDODF(IL)
      IELO=IEL(IL)
      N1H=NLAST(IELO)
      NQ1=NQLODF(IL)
      FRE=ENION(IL)/H
      T=TEMP(ID)
      SQT=SQRT(T)
      ANE=ELEC(ID)
      ANES=EXP(SIXTH*LOG(ANE))
      F00=C00*ANES*ANES*ANES*ANES
      DOP0=CQT*SQT
      QZ=IZ(IELO)
C
C     pseudocontinuum opacity (non-zero in all frequencies);
C     formulated through the dissolved fraction
C
      ITR=ITRA(IL,IU)
      NFR0=NFRODF(KL)
      IF(IMODE.EQ.0) THEN
         DO IJ=1,NFR0
            FRO(IJ)=FROS(IJ,KL)
            SGFR(IJ)=SIGK(FRO(IJ),ITR,1)
            ALAM(IJ)=CAS/FRO(IJ)
         END DO
      END IF
C        
C        function D(nu) - dissolved fraction
C        
c      CALL DWNFR(1,NFR0,FRE,ACOR,ANE,QZ,FRO,DWF)
      DO IJ=1,NFR0
         ABSO0(IJ)=SGFR(IJ)*DWF(IJ)
      END DO
C
C     summation over individual lines
C
      DO J=NQ1,NLMX
         XJ=J
         FXK=F00*XKIJ(KL,J)
         DOP=DOP0/WL0(KL,J)
         DBETA=WL0(KL,J)*WL0(KL,J)/CAS/FXK
         BETAD=DOP*DBETA
         FID=CID*FIJ(KL,J)*DBETA
         CALL DIVSTR(1)
         WPROB=WNHINT(J,ID)
         CALL ODFHST(NFR0,FXK,FID,WPROB,WL0(KL,J),ALAM,SGT)
         DO IJ=1,NFR0
            ABSO0(IJ)=ABSO0(IJ)+SGT(IJ)
         END DO
      END DO
C
C     opacity distribution function in the internal set of frequencies
C
      IF(IMODE.EQ.0) THEN
         ODF0(1)=ABSO0(1)
         IODF(1)=1
         DO IJ=2,NFR0
            ODF0(IJ)=ABSO0(IJ)
            IODF(IJ)=IJ
            IF(ODF0(IJ).LT.ODF0(IJ-1)) THEN
               AB=ODF0(IJ)
               IJODF=IODF(IJ)
               DO IJ0=1,IJ-1
                  IJ1=IJ-IJ0+1
                  IF(ODF0(IJ1).GE.ODF0(IJ1-1)) GO TO 71
                  ODF0(IJ1)=ODF0(IJ1-1)
                  ODF0(IJ1-1)=AB
                  IODF(IJ1)=IODF(IJ1-1)
                  IODF(IJ1-1)=IJODF
               END DO
   71          CONTINUE
            END IF
         if(odf0(ij).gt.0.001) write(6,603) ij,id,odf0(ij)
  603    format(' ij,id,odf0',2i5,1pd10.3)
         END DO
      ELSE
         ODF0(1)=ABSO0(IODF(1))
         IODR(1)=IODF(1)
         DO IJ=2,NFR0
            ODF0(IJ)=ABSO0(IODF(IJ))
            IODR(IJ)=IODF(IJ)
            IF(ODF0(IJ).LT.ODF0(IJ-1)) THEN
               AB=ODF0(IJ)
               IJODF=IODR(IJ)
               DO IJ0=1,IJ-1
                  IJ1=IJ-IJ0+1
                  IF(ODF0(IJ1).GE.ODF0(IJ1-1)) GO TO 86
                  ODF0(IJ1)=ODF0(IJ1-1)
                  ODF0(IJ1-1)=AB
                  IODR(IJ1)=IODR(IJ1-1)
                  IODR(IJ1-1)=IJODF
               END DO
   86          CONTINUE
            END IF
         if(odf0(ij).gt.0.001) write(6,603) ij,id,odf0(ij)
         END DO
         DO IJ=1,NFR0
            IODF(IJ)=IODR(IJ)
         END DO
      END IF
C
C     Reinitialization of the internal frequencies set
C
      FROD(1)=FRO(1)
      IW=IODF(1)
      IF(IW.GT.1 .AND. IW.LT.NFR0) THEN
         W1=FRO(IW-1)-FRO(IW+1)
       ELSE IF (IW.EQ.1) THEN
         W1=FRO(1)-FRO(2)
       ELSE
         W1=FRO(NFR0-1)-FRO(NFR0)
      END IF
      DO IJ=2,NFR0-1
         IW=IODF(IJ)
         IF(IW.GT.1 .AND. IW.LT.NFR0) THEN
            W2=HALF*(FRO(IW-1)-FRO(IW+1))
          ELSE IF (IW.EQ.1) THEN
            W2=HALF*(FRO(1)-FRO(2))
          ELSE
            W2=HALF*(FRO(NFR0-1)-FRO(NFR0))
         END IF
         FROD(IJ)=FROD(IJ-1)-HALF*(W1+W2)
         W1=W2
      END DO
      IJ=NFR0
      IW=IODF(IJ)
      IF(IW.GT.1 .AND. IW.LT.NFR0) THEN
         W2=FRO(IW-1)-FRO(IW+1)
       ELSE IF (IW.EQ.1) THEN
         W2=FRO(1)-FRO(2)
       ELSE
         W2=FRO(NFR0-1)-FRO(NFR0)
      END IF
      FROD(IJ)=FROD(IJ-1)-HALF*(W1+W2)
C
C     Interpolated opacity distribution function to explicit frequencies
C
      DO 150 IJ=2,NFREQ
         IF(FREQ(IJ).GT.FREQ(IJ-1)) RETURN
         ODF(IJ)=0.
         IF(FREQ(IJ).GT.FROD(1).OR.FREQ(IJ).LT.FROD(NFR0)) GO TO 150
         IF(ID.EQ.1) THEN
            IF(FREQ(IJ-1).GT.FROD(1)) I1ODF(IL)=IJ
            I2ODF(IL)=IJ
         END IF
         DO IJ1=2,NFR0
            IJ0=IJ1
            IF(FREQ(IJ).GE.FROD(IJ1)) GO TO 120
         END DO
  120    ODF(IJ)=ODF0(IJ0-1)+(ODF0(IJ0)-ODF0(IJ0-1))/
     *           (FROD(IJ0)-FROD(IJ0-1))*(FREQ(IJ)-FROD(IJ0-1))
  150 CONTINUE
      RETURN
      END