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 'PARAMS.FOR'
      INCLUDE 'MODELP.FOR'
      PARAMETER (UH=1.5)
      PARAMETER (CMAX=2.154D4,CCON=2.0706D-16,TWO=2.D0)
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
      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 50 ION=1,NION
         QZ=IZ(ION)
         CFN=CON/G(NNEXT(ION))
         DCH=0.
         IUPS=IUPSUM(ION)
         SSBF=0.
         USUM(ION)=0.
         nlst=nlast(ion)
         if(ifwop(nlst).ge.0) then
             nl1up=nquant(nlst)+1
          else
             nl1up=nquant(nlst)
         end if
         DO 10 II=NFIRST(ION),NLAST(ION)
            if(ifwop(ii).lt.0) then
               E=EH*QZ*QZ/TK
               SUM=0.
               DO 5 J=nl1up,NLMX
                  XJ=J
                  XI=J*J
                  X=E/XI
                  FI=XI*EXP(X)*WNHINT(J,ID)
                  SUM=SUM+FI
    5          CONTINUE
               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
   10    CONTINUE
C
C     Upper sums
C
         if(ifwop(nlst).lt.0) go to 50
         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)
         EE=ENION(NFIRST(ION))/TK
         if(ee.gt.110.) ee=110.
         CFE=CFN*EXP(EE)
         USUM(ION)=CFE*U-SSBF
         xx=(ssbf-sbf(nfirst(ion)))/sbf(nfirst(ion))
         IF(USUM(ION).LT.0.or.ee.ge.109.or.xx.lt.1.e-7) USUM(ION)=0.
         IF(USUM(ION).LT.0.) USUM(ION)=0.
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 30 J=NQUANT(NLAST(ION))+1,IUPS
            XI=J*J
            X=E/XI
            FI=XI*EXP(X)
            SUM=SUM+FI
   30    CONTINUE
         USUM(ION)=SUM*CON*TWO
C
c        3. occupation probability form
c
         else 
         SUM=0.
         DSUM=0.
         E=EH*QZ*QZ/TK
         DO 40 J=NQUANT(NLAST(ION))+1,NLMX
            XJ=J
            XI=J*J
            X=E/XI
            FI=XI*EXP(X)*WNHINT(J,ID)
            SUM=SUM+FI
   40    CONTINUE
         USUM(ION)=SUM*CON*TWO
         end if
   50 CONTINUE
      RETURN
      END