SUBROUTINE LYMLIN(IJ)
C     =====================
C
C     opacity and emissibvity in first 30 Lyman lines
C
      INCLUDE 'IMPLIC.FOR'
      INCLUDE 'BASICS.FOR'
      INCLUDE 'ATOMIC.FOR'
      INCLUDE 'MODELQ.FOR'
C
      PARAMETER (SIXTH=1./6.,TTW=2./3.,
     *           OS0=0.02654,CPP=4.1412E-16,CPJ=157803.,
     *           C00=1.25E-9,C18=2.997925E18)
      PARAMETER (MLEVL=30)
C
      DIMENSION ABLYM(MDEPTH),EMLYM(MDEPTH),
     *          XKIJL(MLEVL),FIJL(MLEVL),WL0L(MLEVL),FR0L(MLEVL),
     *          F00(MDEPTH),DOP0(MDEPTH),PJ(MLEVL,MDEPTH),
     *          ABTR(MLEVL,MDEPTH),EMTR(MLEVL,MDEPTH),FID(MLEVL,MDEPTH),
     *          AD0(MLEVL,MDEPTH),DIV0(MLEVL,MDEPTH),
     *          DBET0(MLEVL,MDEPTH),BETAD0(MLEVL,MDEPTH)
C
      DATA INIP/1/
C
      IF(IELH.EQ.0) RETURN
C
      FR=FREQ(IJ)
      IF(FR.GT.3.28805E15.OR.FR.LT.1.5E15) RETURN
C
      IF(INIP.EQ.1) THEN
         N0H=NFIRST(IELH)
         N1H=NLAST(IELH)
         NKH=NNEXT(IELH)
         NLH=N1H-N0H
         XII=1.
         DO J=2,30 
            CALL STARK0(1,J,1,XKIJ0,WL00,FIJ0)
            XKIJL(J)=XKIJ0
            FIJL(J)=FIJ0
            WL0L(J)=WL00
            FR0L(J)=C18/WL0L(J)
         END DO
         DO ID=1,ND
            T=TEMP(ID)
            T1=1./T
            SQT=SQRT(T)
            ANE=ELEC(ID)
            ANP=POPUL(NKH,ID)
            F00(ID)=C00*ANE**TTW
            DOP0(ID)=1.E8*SQRT(1.65E8*T+VTURB(ID))
C           
            P0=CPP*ANE*ANP*T1/SQT
            P1=POPUL(N0H,ID)
            DO J=2,30
               X=J*J
               JJ=N0H+J-1
               XJJ=1./JJ
               IF(J.LE.NLH) THEN
                  PJ(J,ID)=POPUL(JJ,ID)
                ELSE
                  PJ(J,ID)=P0*EXP(CPJ/(X*T))*X*WNHINT(J,ID)
               END IF
               ABTR(J,ID)=P1*WNHINT(J,ID)
               EMTR(J,ID)=PJ(J,ID)*XJJ*EXP(CPJ*(XII-XJJ)*T1)
C              
               FXK=F00(ID)*XKIJL(J)
               DBETA=WL0L(J)*WL0L(J)/(C18*FXK)
               FID(J,ID)=OS0*FIJL(J)*DBETA
               DOP=DOP0(ID)/WL0L(J)
               BETAD=DOP*DBETA
               CALL DIVSTR(1)
               AD0(J,ID)=ADH
               DIV0(J,ID)=DIVH
               DBET0(J,ID)=DBETA
               BETAD0(J,ID)=BETAD
           END DO
        END DO
        INIP=0
      END IF
C
      FR=FREQ(IJ)
      WL=C18/FR
      F15=FR*1.E-15
      DO ID=1,ND
         ABLYM(ID)=0.
         EMLYM(ID)=0.
         DO J=2,30
            BETA=DBET0(J,ID)*ABS(FR-FR0L(J))
            BETAD=BETAD0(J,ID)
            ADH=AD0(J,ID)
            DIVH=DIV0(J,ID)
            SG=STARKA(BETA,TWO)*FID(J,ID)
            ABLYM(ID)=ABLYM(ID)+SG*ABTR(J,ID)
            EMLYM(ID)=EMLYM(ID)+SG*EMTR(J,ID)
c      if(wl.gt.1120.0.and.wl.lt.1120.3.and.id.eq.50)
c    *  write(6,600) j,fr,fr0l(j),abs(fr-fr0l(j)),
c    *  dbet0(j,id),beta,fid(j,id),sg,abtr(j,id),
c    *  emtr(j,id),sg*abtr(j,id),ablym(id)
c      if(wl.gt.1120.0.and.wl.lt.1120.3.and.id.eq.50)
c    *  write(6,600) j,fid(j,id),dbet0(j,id),wl0l(j),xkijl(j),
c    *  dop0(id),fijl(j) 
c 600 format('lymsig',i4,1p11e11.3)
         END DO
         XKT=EXP(-4.79928e-11*FR/TEMP(ID))
         XKB=XKT*1.4743E-2*F15*F15*F15
         ABLYM(ID)=ABLYM(ID)-XKT*EMLYM(ID)
         EMLYM(ID)=XKB*EMLYM(ID)
         ABSO1(ID)=ABSO1(ID)+ABLYM(ID)
         EMIS1(ID)=EMIS1(ID)+EMLYM(ID)
c       if(wl.gt.1120.0.and.wl.lt.1120.3)
c    *  write(6,601) ij,wl,ablym(50),emlym(50),xkt,xkb
c 601 format('lymlin',i6,f10.3,1p4e11.3)
      END DO
      RETURN
      END