SUBROUTINE RATSP1
C     =================
C
C     Evaluation of "preconditioned" radiative rates
C
      INCLUDE 'IMPLIC.FOR'
      INCLUDE 'BASICS.FOR'
      INCLUDE 'ATOMIC.FOR'
      INCLUDE 'MODELQ.FOR'
      INCLUDE 'ODFPAR.FOR'
      INCLUDE 'ALIPAR.FOR'
      INCLUDE 'ARRAY1.FOR'
      INCLUDE 'ITERAT.FOR'
      PARAMETER(PGRD=4.1916825D-10)
      DIMENSION EHK(MDEPTH)
C     DIMENSION EHKL(MFREQL)
C
C     zero the rates
C
      DO ID=1,ND
         PRADT(ID)=0.
         PRADA(ID)=0.
         FLRD(ID)=0.
         DO ITR=1,NTRANS
            RRU(ITR,ID)=0.
            RRD(ITR,ID)=0.
         END DO
         DO ITRP=1,NTRPRD
            PJBAR(ITRP,ID)=0.
         END DO
      END DO
      PRD0=0.
C
      LROSS=NDRE.LE.0.AND.ITER.EQ.1.OR.LFIN
      IF(HMIX0.GT.0.) LROSS=.TRUE.
      IF(LROSS) THEN
         DO ID=1,ND
            ABROSD(ID)=0.
            SUMDPL(ID)=0.
         END DO
      END IF
C
      DO 500 IJ=1,NFREQ
         IF(IJX(IJ).EQ.-1) GO TO 500
         FR=FREQ(IJ)
         W0=W0E(IJ)
         WW=W(IJ)
         CALL OPACF1(IJ)
         CALL RTEFR1(IJ)
         IF(LROSS) CALL ROSSTD(IJ)
         FLUXW=W(IJ)*RAD1(1)*FH(IJ)
         GRADF(1,IJ)=FLUXW*ABSO1(1)/DENS(1)*PGRD
         FLRD(1)=FLRD(1)+WW*FH(IJ)*RAD1(1)-WW*HEXTRD(IJ)
         DO ID=2,ND
            DT=UN/(ABSOT(ID)+ABSOT(ID-1))/DELDMZ(ID-1)
            FL=(RAD1(ID)*FAK1(ID)-RAD1(ID-1)*FAK1(ID-1))*DT
            FLRD(ID)=FLRD(ID)+WW*FL
         END DO
         if(ioptab.lt.0) go to 500
C
C        ---------------------
C        Continuum transitions
C        ---------------------
C
         DO ID=1,ND
            EHK(ID)=EXP(-HKT1(ID)*FR)
            if(ilpsct.eq.0) then
               ALAB(ID)=ALI1(ID)/(ABSO1(ID)-ELSCAT(ID))
            else
               ALAB(ID)=ALI1(ID)/ABSO1(ID)
            end if
            DO 10 IBFT=1,NTRANC
               ITR=ITRBF(IBFT)
               SG=CROSS(IBFT,IJ)
               IF(SG.LE.0.) GO TO 10
               II=ILOW(ITR)
               JJ=IUP(ITR)
               IF(IPZERO(II,ID).NE.0.OR.IPZERO(JJ,ID).NE.0) GO TO 10
               JC=ITRA(JJ,II)
               IF(IFWOP(II).GE.0) THEN
                 ICDW=MCDW(ITR)
                 IF(ICDW.GE.1) SG=SG*DWF1(ICDW,ID)
               ELSE
                 IMER=IMRG(II)
                 SG=SGMG(IMER,ID)
               ENDIF
               SGW0=SG*W0
               RLAM=SG*ALAB(ID)
               ELIN=EMTRA(ITR,ID)*EHK(ID)
               RADRES=RAD1(ID)-RLAM*BNUE(IJ)*ELIN
               BNURES=BNUE(IJ)*(UN-RLAM*(ABTRA(ITR,ID)-ELIN))
               RRU(ITR,ID)=RRU(ITR,ID)+SGW0*RADRES
               RRD(ITR,ID)=RRD(ITR,ID)+SGW0*(RADRES+BNURES)*EHK(ID)
   10       CONTINUE
         END DO
C            
C        ----------------
C        Line transitions 
C        ----------------
C             
         IF(ISPODF.EQ.0) THEN
         IF(IJLIN(IJ).GT.0) THEN
C
C        the "primary" line at the given frequency
C
            ITR=IJLIN(IJ)
            II=ILOW(ITR)
            JJ=IUP(ITR)
            DO 50 ID=1,ND
               IF(IPZERO(II,ID).NE.0.OR.IPZERO(JJ,ID).NE.0) GO TO 50
               SG=PRFLIN(ID,IJ)
               SGW0=SG*W0
               RLAM=SG*ALAB(ID)
               ELIN=EMTRA(ITR,ID)*EHK(ID)
               RADRES=RAD1(ID)-RLAM*BNUE(IJ)*ELIN
               BNURES=BNUE(IJ)*(UN-RLAM*(ABTRA(ITR,ID)-ELIN))
               RRU(ITR,ID)=RRU(ITR,ID)+SGW0*RADRES
               RRD(ITR,ID)=RRD(ITR,ID)+SGW0*(RADRES+BNURES)*EHK(ID)
   50       CONTINUE
c
            itrprd=iprd(itr)
            if(itrprd.gt.0) then
            s=un/(0.02654*osc0(itr))
            do id=1,nd
               sg=prflin(id,ij)*s
               pjbar(itrprd,id)=pjbar(itrprd,id)+sg*w(ij)*rad1(id)
            end do
            end if 
c
         ENDIF
         IF(NLINES(IJ).LE.0) GO TO 100
C
C        the "overlapping" lines at the given frequency
C
         DO 90 ILINT=1,NLINES(IJ)
            ITR=ITRLIN(ILINT,IJ)
            if(linexp(itr)) go to 90
            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
            II=ILOW(ITR)
            JJ=IUP(ITR)
            DO 80 ID=1,ND
               IF(IPZERO(II,ID).NE.0.OR.IPZERO(JJ,ID).NE.0) GO TO 80
               SG=A1*PRFLIN(ID,IJ1)+A2*PRFLIN(ID,IJ0)
               SGW0=SG*W0
               RLAM=SG*ALAB(ID)
               ELIN=EMTRA(ITR,ID)*EHK(ID)
               RADRES=RAD1(ID)-RLAM*BNUE(IJ)*ELIN
               BNURES=BNUE(IJ)*(UN-RLAM*(ABTRA(ITR,ID)-ELIN))
               RRU(ITR,ID)=RRU(ITR,ID)+SGW0*RADRES
               RRD(ITR,ID)=RRD(ITR,ID)+SGW0*(RADRES+BNURES)*EHK(ID)
   80       CONTINUE
c
            itrprd=iprd(itr)
            if(itrprd.gt.0) then
               s=un/(0.02654*osc0(itr))
               do id=1,nd
                  SG=(A1*PRFLIN(ID,IJ1)+A2*PRFLIN(ID,IJ0))*s
                  pjbar(itrprd,id)=pjbar(itrprd,id)+sg*w(ij)*rad1(id)
               end do
            end if 
c
   90    CONTINUE
  100 CONTINUE
C
C     Opacity sampling option
C
      ELSE
        IF(NLINES(IJ).LE.0) GO TO 200
        DO 190 ILINT=1,NLINES(IJ)
           ITR=ITRLIN(ILINT,IJ)
           KJ=IJ-IFR0(ITR)+KFR0(ITR)
           INDXPA=IABS(INDEXP(ITR))
           II=ILOW(ITR)
           JJ=IUP(ITR)
           IF(INDXPA.NE.3 .AND. INDXPA.NE.4) THEN
              DO 150 ID=1,ND
                 IF(IPZERO(II,ID).NE.0.OR.IPZERO(JJ,ID).NE.0) GO TO 150
                 SG=PRFLIN(ID,KJ)
                 SGW0=SG*W0
                 RLAM=SG*ALAB(ID)
                 ELIN=EMTRA(ITR,ID)*EHK(ID)
                 RADRES=RAD1(ID)-RLAM*BNUE(IJ)*ELIN
                 BNURES=BNUE(IJ)*(UN-RLAM*(ABTRA(ITR,ID)-ELIN))
                 RRU(ITR,ID)=RRU(ITR,ID)+SGW0*RADRES
                 RRD(ITR,ID)=RRD(ITR,ID)+SGW0*(RADRES+BNURES)*EHK(ID)
  150         CONTINUE
           ELSE
              DO 160 ID=1,ND
                 IF(IPZERO(II,ID).NE.0.OR.IPZERO(JJ,ID).NE.0) GO TO 160
                 KJD=JIDI(ID)
                 SG=EXP(XJID(ID)*SIGFE(KJD,KJ)+
     *              (UN-XJID(ID))*SIGFE(KJD+1,KJ))
                 SGW0=SG*W0
                 RLAM=SG*ALAB(ID)
                 ELIN=EMTRA(ITR,ID)*EHK(ID)
                 RADRES=RAD1(ID)-RLAM*BNUE(IJ)*ELIN
                 BNURES=BNUE(IJ)*(UN-RLAM*(ABTRA(ITR,ID)-ELIN))
                 RRU(ITR,ID)=RRU(ITR,ID)+SGW0*RADRES
                 RRD(ITR,ID)=RRD(ITR,ID)+SGW0*(RADRES+BNURES)*EHK(ID)
  160         CONTINUE
           END IF
  190   CONTINUE
  200 CONTINUE
      END IF
  500 CONTINUE
C
      DO ID=1,ND
         IF(CRSW(ID).NE.UN) THEN
            DO ITR=1,NTRANS
               RRU(ITR,ID)=RRU(ITR,ID)*CRSW(ID)
               RRD(ITR,ID)=RRD(ITR,ID)*CRSW(ID)
            END DO
         END IF
      END DO
C
C     radiation pressure
C
      PRDX=1.
      DO ID=1,ND
         PRADT(ID)=PRADT(ID)*PCK
         PRADA(ID)=PRADA(ID)*PCK
         if(prada(id).gt.0.) PRDR=PRADT(ID)/PRADA(ID)
         IF(PRDR.LT.PRDX) PRDX=PRDR
      END DO
      PRD0=PRD0/DENS1(1)*DM(1)*PCK
      IF(LFIN) WRITE(10,1100) PRDX,ITER
 1100 FORMAT(' PRAD MIN RATIO ',F10.6,I4)
C
C     Rosseland mean opacity
C
      IF(LROSS) THEN
      DO ID=1,ND
         ABROSD(ID)=SUMDPL(ID)/(ABROSD(ID)*DENS(ID))
      END DO
      END IF
            
      RETURN
      END