SUBROUTINE RADPRE
C     =================
C
C     radiative acceleration
C
C      exclude automatically the strongest lines of total
C       radiation pressure (Keyword XGRAD)
C
C      depth-dependent criterion
C
C      automatic explicit frequencies if XGRAD>=0
C
      INCLUDE 'IMPLIC.FOR'
      INCLUDE 'BASICS.FOR'
      INCLUDE 'ATOMIC.FOR'
      INCLUDE 'MODELQ.FOR'
      INCLUDE 'ODFPAR.FOR'
      INCLUDE 'ALIPAR.FOR'
      DIMENSION GRADA(MDEPTH),PRID(MDEPTH),GRADI(MFREQ)
      DIMENSION PGT(MDEPTH),GGRT(MDEPTH),XGRD(MDEPTH),XGRD0(10)
      DIMENSION XGRD1(20),XGRD2(20)
      DIMENSION IIGR(MFREQ)
      PARAMETER(PGRD=4.1916825D-10)
      DATA XGRD0/0.1,0.3,0.5,0.7,0.9,0.92,0.94,0.96,0.98,0.99/
      DATA XGRD1/0.1,0.2,0.3,0.4,0.5,0.6,0.65,0.7,0.75,0.8,
     &           0.85,0.9,0.92,0.94,0.96,0.98,0.99,0.99,0.99,0.99/
      DATA XGRD2/0.1,0.2,0.3,0.4,0.45,0.5,0.55,0.6,0.65,0.7,
     &           0.75,0.8,0.84,0.87,0.9,0.93,0.95,0.97,0.98,0.99/
      
      IF(XGRAD.EQ.0) THEN
        DO ID=1,10
        XGRD(ID)=XGRD0(ID)
      END DO
      DO ID=11,ND
        XGRD(ID)=XGRD(ID-1)
      END DO
      ELSE IF(XGRAD.EQ.-1.) THEN
        DO ID=1,20
        XGRD(ID)=XGRD1(ID)
      END DO
      DO ID=21,ND
        XGRD(ID)=XGRD(ID-1)
      END DO
      ELSE IF(XGRAD.EQ.-2.) THEN
        DO ID=1,20
          XGRD(ID)=XGRD2(ID)
        END DO
        DO ID=21,ND
          XGRD(ID)=XGRD(ID-1)
        END DO
      ELSE
        DO ID=1,ND
        XGRD(ID)=XGRAD
      END DO
      END IF

C     Acceleration due to gas and turbulent pressure

      DO ID=1,ND
        PGAS=(DENS(ID)/WMM(ID)+ELEC(ID))*BOLK*TEMP(ID)
        PGT(ID)=PGAS+HALF*DENS(ID)*VTURB(ID)*VTURB(ID)
      END DO
      DO ID=2,ND
        GGRT(ID)=(PGT(ID)-PGT(ID-1))/(DM(ID)-DM(ID-1))
      END DO
      GGRT(1)=GGRT(2)

C     Compute total radiative acceleration at every depth points
 
      DO ID=1,ND
        GRAD(ID)=0.
      GRADA(ID)=0.
        PRADT(ID)=0.
      END DO
      PRD0=0.
      DO ID=2,ND
      PRID(ID)=PGRD/(DM(ID)-DM(ID-1))
      END DO
      PGRD1=PGRD/DENS(1)
      DO IJ=1,NFREQ
        CALL OPACF1(IJ)
        CALL RTEFR1(IJ)
      FLUXW=W(IJ)*(RAD1(1)*FH(IJ)-HEXTRD(IJ))
      GRADF(1,IJ)=FLUXW*ABSO1(1)*PGRD1
        GRADA(1)=GRADA(1)+GRADF(1,IJ)
      DO ID=2,ND
          FRD=FAK1(ID)*RAD1(ID)-FAK1(ID-1)*RAD1(ID-1)
          GRADF(ID,IJ)=W(IJ)*FRD*PRID(ID)
          GRADA(ID)=GRADA(ID)+GRADF(ID,IJ)
      END DO
      END DO
      DO ID=1,ND
        GGRT(ID)=GRADA(ID)+GGRT(ID)
      END DO
C
C     radiation pressure
C
      DO ID=1,ND
         PRADT(ID)=PRADT(ID)*PCK
      END DO
      PRD0=PRD0/DENS1(1)*DM(1)*PCK
C
C     Depth-dependent rejection: set up LSKIP(ID,IJ)
C      only if XGRAD<=0

      NFE=0
      DO ID=1,ND
        XGR0=GRAV*ABS(XGRD(ID))
        DO IJ=1,NFREQ
        GRADI(IJ)=GRADF(ID,IJ)
          LSKIP(ID,IJ)=.FALSE.
        END DO
        if(ifprad.eq.0) then
            do ij=1,nfreq
               lskip(id,ij)=.true.
            end do
        end if
        CALL INDEXX(NFREQ,GRADI,IIGR)
      GRAD(ID)=GRADA(ID)
      GGRT0=GGRT(ID)
      
      IF(XGRAD.GT.0. .AND. ID.GT.1) THEN
        DO IJ=1,NFREQ
          LSKIP(ID,IJ)=LSKIP(1,IJ)
          IF(LSKIP(ID,IJ)) THEN
            GGRT0=GGRT0-GRADI(IJ)
            GRAD(ID)=GRAD(ID)-GRADI(IJ)
          END IF
        END DO
      GO TO 110
        END IF
      
      IF(ID.GE.ND-1) GO TO 110

        IJR=NFREQ
        NSK=0
        DO WHILE(GRAD(ID).GT.XGR0 .AND. IJR.GT.0)
          IJ=IIGR(IJR)
        IF(IJLIN(IJ).EQ.0 .AND. NLINES(IJ).EQ.0) GO TO 99
          LSKIP(ID,IJ)=.TRUE.
        GGRT0=GGRT0-GRADI(IJ)
        GRAD(ID)=GRAD(ID)-GRADI(IJ)
        NSK=NSK+1
        IF(XGRD(ID).GE.0.) GO TO 99
          IF(NFE.LT.10) THEN
            IF(ISPODF.EQ.0) THEN
              ITR=IJLIN(IJ)
            IF(ITR.EQ.0) GO TO 99 
              INDXPA=IABS(INDEXP(ITR))
              DX=(FREQ(IJ)-FREQ(IJ+1))*0.25
              DZ=ABS(FREQ(IJ)-FR0(ITR))
              IF(DZ.LT.DX .AND. INDXPA.EQ.1) THEN
                IF(INDEXP(ITR).LT.0) THEN
                  INDEXP(ITR)=-9
                ELSE
                  INDEXP(ITR)=9
                END IF
                IF(.NOT.LEXP(ITR)) THEN
                  LEXP(ITR)=.TRUE.
                  NFREQE=NFREQE+1
                  if(nfreqe.gt.mfrex) 
     *            CALL QUIT('nfreqe.gt.mfrex',nfreqe,mfrex)
                  NN=NN+1
                  IJALI(IJ)=0
                  IJEX(IJ)=NFREQE
                  IJFR(NFREQE)=IJ
                  IJX(IJ)=1
                  WC(IJ)=0.
                  WRITE(10,612) FREQ(IJ),ITR,IJ,NFREQE
                  NFE=NFE+1
                END IF
              END IF
            ELSE
              DO 100 ILINT=1,NLINES(IJ)
                ITR=ITRLIN(ILINT,IJ)
                INDXPA=IABS(INDEXP(ITR))
                DX=(FREQ(IJ)-FREQ(IJ+1))*0.25
                DZ=ABS(FREQ(IJ)-FR0(ITR))
                IF(DZ.GT.DX .OR. INDXPA.NE.1) GOTO 100
                IF(INDEXP(ITR).LT.0) THEN
                  INDEXP(ITR)=-9
                ELSE
                  INDEXP(ITR)=9
                END IF
                IF(.NOT.LEXP(ITR)) THEN
                  LEXP(ITR)=.TRUE.
                  NFREQE=NFREQE+1
                  if(nfreqe.gt.mfrex) 
     *            CALL QUIT('nfreqe.gt.mfrex',nfreqe,mfrex)
                  NN=NN+1
                  IJALI(IJ)=0
                  IJEX(IJ)=NFREQE
                  IJFR(NFREQE)=IJ
                  IJX(IJ)=1
                  WC(IJ)=0.
                  WRITE(10,612) FREQ(IJ),ITR,IJ,NFREQE
                  NFE=NFE+1
                END IF
  100         CONTINUE
            END IF
          END IF
   99     IJR=IJR-1
        END DO
      
  110   IF(ID.EQ.1) THEN
        TAUR=HALF*DEDM1*ABROSD(ID)*DENS(ID)
      ELSE
          DTAUR=DELDM(ID-1)*(ABROSD(ID)+ABROSD(ID-1))
          TAUR=TAUR+DTAUR
      END IF
      
      rgrt=ggrt(id)/ggrt0
      if(rgrt.gt.0.) then
        rgrt=dlog10(rgrt)
      else
        rgrt=-9.
      end if
      
      END DO
  612 FORMAT(' AUTOMATIC EXPLICIT FREQ.  ',1PE12.6,3I8)

      RETURN
      END