SUBROUTINE NLTE(IL,ILW,IUN,GI,GJ)
C     ===========================================
C
C     Control procedure for the NLTE option
C
      INCLUDE 'PARAMS.FOR'
      INCLUDE 'MODELP.FOR'
      INCLUDE 'LINDAT.FOR'
      COMMON/NLTPOP/PNLT(MATOM,MION,MDEPTH)
      PARAMETER (UN     = 1., 
     *           C3     = 1.4387886,
     *           XET    = 8067.6,
     *           XET3   = XET*C3)
C
C     CALCULATION OF THE 
C     CENTRAL OPACITY (ABCENT) AND THE LINE SOURCE FUNCTION  (SLIN)
C
      if(gi.le.0..or.gj.le.0.) return
      ILNLT=INDNLT(IL)
      IF(ILNLT.LE.0) RETURN
      IAT=INDAT(IL)/100
      ION=MOD(INDAT(IL),100)
      EGF=EXP(GF0(IL))
      BNU=BN*(FREQ0(IL)*1.E-15)**3
      DP0=3.33564E-11*FREQ0(IL)
      DP1=1.651E8/AMAS(IAT)
      IF(ILW.LE.0) GO TO 100
C
C     line is a transition between explicit levels of the
C     input model
C
      NKI=NNEXT(IEL(ILW))
      DO 60 ID=1,ND
         T=TEMP(ID)
         COR=1.
         PP=PNLT(IAT,ION,ID)
         IF(ILW.GT.0) THEN
            PI=POPUL(ILW,ID)/G(ILW)
          ELSE
            PI=PP*EXP((ENEV(IAT,ION)*XET3-EXCL0(IL))/T)
         END IF
         IF(IUN.GT.0) THEN
            PJ=POPUL(IUN,ID)/G(IUN)
            cor=(excu0(il)-excl0(il)+
     *         (enion(iun)-enion(ilw))/1.38054e-16)/t
            cor=exp(cor)
          ELSE
            PJ=PP*EXP((ENEV(IAT,ION)*XET3-EXCU0(IL))/T)
         END IF
         if(pj.gt.0.) then
            X=PI/PJ*cor
          else
            x=un
         end if
         IF(X.EQ.UN) X=EXP(4.79928E-11*FREQ0(IL)/T)
         DOP=DP0*SQRT(DP1*T+VTURB(ID))
         SLIN(ILNLT,ID)=BNU/(X-UN)
         if(pi.gt.0.) ABCENT(ILNLT,ID)=PI*(UN-UN/X)*EGF/DOP
   60 CONTINUE
      RETURN
C
C     Approximate NLTE for resonance lines - second order escape
C     probablity theory form of the source function
C
C     Optical depth scale
C
  100 CONTINUE
      ALMIL=2.997925E17/FREQ0(IL)
      HKF=HK*FREQ0(IL)
      DO 110 ID=1,ND
         T=TEMP(ID)
         DOP=DP0*SQRT(DP1*T+VTURB(ID))
         X=EXP(HKF/T)
         ABCENT(ILNLT,ID)=EGF*EXP(-EXCL0(IL)/T)*RRR(ID,ION,IAT)/
     *                    DOP*(1.-1./X)
         AB=ABSTD(ID)+ABCENT(ILNLT,ID)*1.77245
         if(ifwin.gt.0) 
     *   AB=ABSTDW(IJCONT(IL),ID)+ABCENT(ILNLT,ID)*1.77245
         IF(ID.EQ.1) THEN
            ABM=AB/DENS(1)
            TAU=0.5*DM(1)*ABM
          ELSE
            AB0=AB/DENS(ID)
            TAU=TAU+0.5*(DM(ID)-DM(ID-1))*(AB0+ABM)
            ABM=AB0
         END IF
C
C        approximate epsilon after Kastner
C
         E=EPS(T,ELEC(ID),ALMIL,ION,IUN)
         XK2=XK2DOP(TAU)
         SLIN(ILNLT,ID)=SQRT(E/(E+(1.-E)*XK2))*BNU/(X-1.)
  110 CONTINUE
      RETURN
      END