SUBROUTINE SABOLF(ID)
C     =====================
C
C     Saha-Boltzmann factors (SBF)
C     and "upper sums" - sum of Saha-Boltzmann factors for upper, LTE,
C     levels which are not included explicitly (USUM), and derivatives
C     wrt. temperature (T) and electron density (DUSUMN)
C
C     Input: ID  - depth index
C
      INCLUDE 'IMPLIC.FOR'
      INCLUDE 'BASICS.FOR'
      INCLUDE 'ATOMIC.FOR'
      INCLUDE 'MODELQ.FOR'
      PARAMETER (UH=1.5)
      PARAMETER (CMAX=2.154D4,CCON=2.0706D-16)
C
C     DCHI - approximate lowering of ionization potential for neutrals
C      Actual lowering is DCHI*effective charge, and is considered only
C      if IUPSUM(ION).GT.0
C
c
      if(ioptab.lt.0) return
c
      T=TEMP(ID)
      SQT=SQRT(T)
      ANE=ELEC(ID)
      STANE=SQRT(T/ANE)
      XMAX=CMAX*SQRT(STANE)
      TK=BOLK*T
      CON=CCON/T/SQT
C
C     Saha-Boltzmann factors
C
      DO ION=1,NION
         QZ=IZ(ION)
         CFN=CON/G(NNEXT(ION))
         DCH=0.
         IUPS=IUPSUM(ION)
         SSBF=0.
         DSSBFT=0.
         USUM(ION)=0.
         DUSUMT(ION)=0.
         DUSUMN(ION)=0.
         nlst=nlast(ion)
         if(ifwop(nlst).ge.0) then
             nl1up=nquant(nlst)+1
          else
             nl1up=nquant(nlst)
         end if
         DO II=NFIRST(ION),NLAST(ION)
            if(ifwop(ii).lt.0) then
               E=EH*QZ*QZ/TK
               SUM=0.
               DO J=nl1up,NLMX
                  XJ=J
                  XI=J*J
                  X=E/XI
                  FI=XI*EXP(X)*WNHINT(J,ID)
                  SUM=SUM+FI
               END DO
               g(ii)=sum*two
               gmer(imrg(ii),id)=g(ii)
            end if
            X=ENION(II)/TK
            if(x.gt.110.) x=110.
            SB=CFN*G(II)*EXP(X)
            SBF(II)=SB
            SSBF=SSBF+SB
            DSBF(II)=-(UH+ENION(II)/TK)/T
            DSSBFT=DSSBFT+SB*DSBF(II)
         END DO
C
C     Upper sums
C
         if(ifwop(nlst).ge.0) then
         IF(ION.EQ.IELHM) THEN
            USUM(ION)=0.
            DUSUMT(ION)=0.
            DUSUMN(ION)=0.
            GO TO 50
         END IF
         IF(IUPS.EQ.0) THEN
C
C     1. More exact approach - using (exact) partition functions
C
            IAT=NUMAT(IATM(NFIRST(ION)))
            XMX=XMAX*SQRT(QZ)
            CALL PARTF(IAT,IZ(ION),T,ANE,XMX,U,DUT,DUN)
            EE=ENION(NFIRST(ION))/TK
            IF(EE.GT.110.) EE=110.
            CFE=CFN*EXP(EE)
            USUM(ION)=CFE*U-SSBF
            DUSUMT(ION)=CFE*(DUT-U/T*(UH+EE))-DSSBFT
            DUSUMN(ION)=CFE*DUN
            xx=(ssbf-sbf(nfirst(ion)))/sbf(nfirst(ion))
            IF(USUM(ION).LT.0.or.ee.ge.109.or.xx.lt.1.e-7) THEN 
                USUM(ION)=0.
                DUSUMT(ION)=0.
                DUSUMN(ION)=0.
            END IF
C
C     2. Approximate approach - summation over fixed number of upper
C        levels, assumed hydrogenic (ie. their ionization energy and
C        statistical weight hydrogenic)
C
          ELSE IF(IUPS.GT.0) THEN
            SUM=0.
            DSUM=0.
            E=EH*QZ*QZ/TK
            DO J=NQUANT(NLAST(ION))+1,IUPS
               XI=J*J
               X=E/XI
               FI=XI*EXP(X)
               SUM=SUM+FI
               DSUM=DSUM-FI*(UH+X)/T
            END DO
            USUM(ION)=SUM*CON*TWO
            DUSUMT(ION)=DSUM*CON*TWO
            DUSUMN(ION)=0.
C
c        3. occupation probability form
c
          ELSE 
            SUM=0.
            DSUM=0.
            E=EH*QZ*QZ/TK
            DO J=NQUANT(NLAST(ION))+1,NLMX
               XJ=J
               XI=J*J
               X=E/XI
               FI=XI*EXP(X)*WNHINT(J,ID)
               SUM=SUM+FI
               DSUM=DSUM-FI*(UH+X)/T
            END DO
               USUM(ION)=SUM*CON*TWO
            DUSUMT(ION)=DSUM*CON*TWO
            DUSUMN(ION)=0.
         END IF
         end if
   50 CONTINUE
      END DO
      RETURN
      END