SUBROUTINE PRD(IJ)
c     ==================
c     
c     modification of the line emission coefficient 
c     and the scattering coefficient in the case of PRD
c
      INCLUDE 'IMPLIC.FOR'
      INCLUDE 'BASICS.FOR'
      INCLUDE 'ATOMIC.FOR'
      INCLUDE 'MODELQ.FOR'
      INCLUDE 'ITERAT.FOR'
      parameter(a21=4.699e8,pi2=6.28318531,gr=2.*4.8e-8)
c
      if(ij.gt.0) then
c     if(ilam.le.1) return
      FR=FREQ(IJ)
      IF(ISPODF.EQ.0) THEN
      IF(IJLIN(IJ).GT.0) THEN
C
C     the "primary" line at the given frequency
C
         ITR=IJLIN(IJ)
         ITRPRD=IPRD(ITR)
         IF(ITRPRD.GT.0) THEN
         DFR=ABS(FREQ(IJ)-FR0(ITR))
         if(ilow(itr).eq.nfirst(ielh)) then
            omeg=dfr*pi2
            gra=a21+gr*popul(nfirst(ielh),id)
            do id=1,nd
               coher(itrprd,id)=a21/
     *             (gra+gami(2,'elec',omeg,temp(id),elec(id)))
            end do
         end if
         DO ID=1,ND
            SG=PRFLIN(ID,IJ)
            IF(DFR/DOPTR(ITRPRD,ID).LE.XPDIV) SG=0.
            SCALIN=SG*ABTRA(ITR,ID)*COHER(ITRPRD,ID)
            SCAT1(ID)=SCAT1(ID)+SCALIN
            scem=sg*emtra(itr,id)*coher(itrprd,id)*xkfb(id)
c           EMIS1(ID)=EMIS1(ID)-SCALIN*RJBAR(ITRPRD,ID)
            EMIS1(ID)=EMIS1(ID)-SCEM
         END DO 
         END IF
      END IF
      IF(NLINES(IJ).GT.0) THEN
C
C     the "overlapping" lines at the given frequency
C
      DO 100 ILINT=1,NLINES(IJ)
         ITR=ITRLIN(ILINT,IJ)
         ITRPRD=IPRD(ITR)
         IF(ITRPRD.EQ.0) GO TO 100
         IJ0=IFR0(ITR)
         DO IJT=IJ0,IFR1(ITR)
            IF(FREQ(IJT).LE.FR) THEN
               IJ0=IJT
               GO TO 70
            END IF
         END DO
   70    IJ1=IJ0-1
         A1=(FR-FREQ(IJ0))/(FREQ(IJ1)-FREQ(IJ0))
         A2=UN-A1
         DFR=ABS(FREQ(IJ)-FR0(ITR))
         if(ilow(itr).eq.nfirst(ielh)) then
            omeg=dfr*pi2
            gra=a21+gr*popul(nfirst(ielh),id)
            do id=1,nd
               coher(itrprd,id)=a21/
     *              (gra+gami(2,'elec',omeg,temp(id),elec(id)))
            end do
         end if
         DO ID=1,ND
            SG=A1*PRFLIN(ID,IJ1)+A2*PRFLIN(ID,IJ0)
            IF(DFR/DOPTR(ITRPRD,ID).LE.XPDIV) SG=0.
            SCALIN=SG*ABTRA(ITR,ID)*COHER(ITRPRD,ID)
            scem=sg*emtra(itr,id)*coher(itrprd,id)*xkfb(id)
            SCAT1(ID)=SCAT1(ID)+SCALIN
c           EMIS1(ID)=EMIS1(ID)-SCALIN*RJBAR(ITRPRD,ID)
            EMIS1(ID)=EMIS1(ID)-SCEM
         END DO 
  100 CONTINUE
      END IF
C
C     Opacity sampling option
C
      ELSE
      IF(NLINES(IJ).GT.0) THEN
      DO 300 ILINT=1,NLINES(IJ)
         ITR=ITRLIN(ILINT,IJ)
         ITRPRD=IPRD(ITR)
         IF(ITRPRD.EQ.0) GO TO 300
         KJ=IJ-IFR0(ITR)+KFR0(ITR)
         INDXPA=IABS(INDEXP(ITR))
         IF(INDXPA.NE.3 .AND. INDXPA.NE.4) THEN
         DFR=ABS(FREQ(IJ)-FR0(ITR))
         if(ilow(itr).eq.nfirst(ielh)) then
            omeg=dfr*pi2
            gra=a21+gr*popul(nfirst(ielh),id)
            do id=1,nd
               coher(itrprd,id)=a21/
     *         (gra+gami(2,'elec',omeg,temp(id),elec(id)))
            end do
         end if
         DO ID=1,ND
            SG=PRFLIN(ID,KJ)
            IF(DFR/DOPTR(ITRPRD,ID).LE.XPDIV) SG=0.
            SCALIN=SG*ABTRA(ITR,ID)*COHER(ITRPRD,ID)
            SCAT1(ID)=SCAT1(ID)+SCALIN
            EMIS1(ID)=EMIS1(ID)-SCEM
         END DO
         END IF
  300 CONTINUE
      END IF
      END IF
      RETURN
c


      end if
c
      do itrp=1,ntrprd
         itr=itrtot(itrp)
         aji=osc0(itr)*g(ilow(itr))/g(iup(itr))*7.42163e-22*
     *       fr0(itr)**2
         omeg=0.
         do id=1,nd
            t=temp(id)
            ane=elec(id)
            call dopgam(itr,id,t,dop,agam) 
            doptr(itrp,id)=dop
            coher(itrp,id)=0.99
            if(agam.gt.0.) coher(itrp,id)=aji/(12.5664*dop*agam)
            if(coher(itrp,id).gt.0.999) coher(itrp,id)=0.999
c            
c           special expression for Lyman-alpha
c
            coher(itrp,id)=aji/(aji+9.8e-8*popul(nfirst(ielh),id)+
     *                     0.667*(gami(2,'iont',omeg,t,ane)+
     *                     gami(2,'elec',omeg,t,ane)))
            rjbar(itrp,id)=pjbar(itrp,id)
         end do
      end do
      return
      END