      PROGRAM MDRGF
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C  EVALUATES BURGESS DR GENERAL FORMULA I OR II.
C
C  Z=RESIDUAL CHARGE OF THE TARGET ION, .LT.0 TO FLAG H-LIKE, ELSE .GT.0
C  DE=CORE EXCITATION ENERGY (RYD).
C  FUP=UPWARD ABSORPTION OSCILLATOR STRENGTH.
C  TMIN,TMAX=TEMPERATURE RANGE IN UNITS 10**6 KELVIN
C            IF TMIN.GT.TMAX THEN STANDARD ADAS TEMPS USED.
C  NTEMP=NO. OF TEMPS REQUIRED, >1 .
C  NMIN .GT. NMAX .OR. NMIN .LT. 1 , GF I
C  .ELSE. SUM OVER N FROM NMIN TO NMAX FOR GF II.
C
C  MOST ALTERNATIVE INPUTS TO FUP EXIST:
C  S=LINE STRENGTH
C  AR=RADIATIVE RATE
C  GF=SYMMETRIC F-VALUE
C  FDOWN=DOWNWARD EMISSION OSCILLATOR STRENGTH
C
C  THESE REQUIRE STAT WEIGHT INFO
C  GLOW=LOWER STAT WEIGHT
C  GUP=UPPER STAT WEIGHT
C
C  UNITS=IONIZATION POTENTIAL OF HYDROGEN. 
C        DEFAULT=1.0, ASSUMES DE IN RYDBERGS.
C
C
      PARAMETER (NDIM1=501)
      PARAMETER (NDIM28=19)
C
      LOGICAL BGFI,BHLIKE
      CHARACTER TEXT*80
C
      DIMENSION ALFS(NDIM1),TEMP(NDIM1)
      DIMENSION THTIC(NDIM28)
C
      DATA
     XTHTIC/1.0E1,2.0E1,5.0E1,1.0E2,2.0E2,5.0E2,1.0E3,2.0E3,5.0E3,1.0E4
     X     ,2.0E4,5.0E4,1.0E5,2.0E5,5.0E5,1.0E6,2.0E6,5.0E6,1.0E7/
C
      NAMELIST/SGF/ Z,DE,FUP,TMIN,TMAX,NTEMP,NMIN,NMAX
     X,FDOWN,GF,S,AR,GLOW,GUP,UNITS,TEXT
C
      OPEN(UNIT=5,FILE='drgfin',STATUS='UNKNOWN')
      OPEN(UNIT=6,FILE='drgfout',STATUS='UNKNOWN')
C
      DO 1 I=1,100
      ALFS(I)=0.0
   1  CONTINUE
C
   2  Z=0.0
      TMIN=100.
      TMAX=10.
      NTEMP=31
      NMIN=-1
C
      FUP=0.0
      FDOWN=0.0
      GF=0.0
      S=0.0
      AR=0.0
      GLOW=0.0
      GUP=0.0
C
      TEXT='  '
      UNITS=1.0                                      !RYDBERGS
C
      READ (5,SGF,END=4)
C
      IF(Z.EQ.0.0)GO TO 4                            !EXIT
C
      IF(DE.EQ.0.0)STOP '***ERROR: DE IS ZERO!'
      DE=DE/UNITS
      IF(DE.LT.0.0)DE=-DE
C
      IF(FUP.EQ.0.0)THEN                           !LOOK FOR OTHER INPUT
        IF(S.NE.0.0)THEN                                   !CONVERT TO FUP
          IF(GLOW.EQ.0.0)STOP'S INPUT REQUIRES LOWER STAT WEIGHT (GLOW)'
          FUP=DE*S/(3*GLOW)
        ENDIF
        IF(GF.NE.0.0)THEN
          IF(GLOW.EQ.0.0)STOP'GF INPUT REQUIRES LOWER STAT WEIGHT(GLOW)'
          FUP=GF/GLOW
        ENDIF
        IF(FDOWN.NE.0.0)THEN
          IF(GLOW*GUP.EQ.0.0)STOP
     X    'FDOWN INPUT REQUIRES LOWER AND UPPER STAT WEIGHTS(GLOW,GUP)'
          FUP=FDOWN*GUP/GLOW
        ENDIF
        IF(AR.NE.0.0)THEN
          IF(GLOW*GUP.EQ.0.0)STOP
     X    'AR INPUT REQUIRES LOWER AND UPPER STAT WEIGHTS(GLOW,GUP)'
          FUP=1.245E-10*AR*GUP/(GLOW*DE**3)
        ENDIF
        IF(FUP.EQ.0.0)STOP'***ERROR: FUP IS ZERO!'
      ENDIF
C
      WRITE(6,104)Z,DE,FUP,NTEMP,NMIN,NMAX
C
      IF(NTEMP.GT.NDIM1) STOP 'INCREASE NDIM1 TO NTEMP'
      BGFI=NMIN.GT.NMAX.OR.NMIN.LT.1
      BHLIKE=Z.LT.0.0
      IF(BHLIKE)Z=-Z
      IF(NTEMP.LT.2)NTEMP=2
C
      ZP=Z+1.0
      A0=1.0+0.015*Z**3/ZP**2
      ET=0.158*DE/A0
      IF(.NOT.BGFI)GO TO 9
      EPS=DE/ZP**2
      X=ZP*EPS
      AX=SQRT(X)/(1.0+0.105*X+0.015*X**2)
      BZ=SQRT(Z)*SQRT(ZP)*ZP**2/SQRT(Z**2+13.4)
      HZ=1.0
      IF(BHLIKE)THEN
        Z11=Z-11
        T=1.0+4.5E-5*Z11*Z11*Z11
        HZ=0.84+0.5/ZP**2+0.03*Z11/T
        WRITE(6,105)HZ
      ENDIF
C
   9  IF(TMIN.LE.TMAX)THEN
        IF(TMIN.LT.1.0E-4)TMIN=1.0E-4
        T1=LOG10(TMIN)
        IF(TMAX.LT.1.0E-4)TMAX=1.0E-4
        IF(TMIN.EQ.TMAX)TMAX=10.0*TMIN
        T2=LOG10(TMAX)
        TN=NTEMP-1
        DT=T2-T1
        DT=DT/TN
        DO I=1,NTEMP
          TI=I-1
          TEMP(I)=TI*DT+T1
          TEMP(I)=10.0**TEMP(I)
        ENDDO
      ELSE
        NTEMP=NDIM28
        DO I=1,NTEMP
          TEMP(I)=THTIC(I)*Z**2/1.E6
        ENDDO
      ENDIF
      WRITE(6,101)
C
      IF(BGFI)GO TO 8
      T0=2.0E6*Z*Z/DE**3
      ALFN=0.0
      DO 7 N=NMIN,NMAX
        T3=N**3
        T5=N**5
        T2=1.0+T0/T5
        T2=1.0/T2
        T2=T2/T3
        ALFN=ALFN+T2
   7  CONTINUE
C
   8  DO 3 I=1,NTEMP
        IF(.NOT.BGFI)ALF=ALFN*8.0E-9*Z*Z*FUP/DE
        IF(BGFI)ALF=3.0E-12*BZ*FUP*AX*HZ
        ALF=ALF*EXP(-ET/TEMP(I))/(TEMP(I)*SQRT(TEMP(I)))
        ALFS(I)=ALFS(I)+ALF
        T6=TEMP(I)*1.E6
        WRITE(6,102)ALOG10(T6),T6,ALF
   3  CONTINUE
      NTEMP0=NTEMP
      GO TO 2
   4  WRITE(6,106)TEXT
      WRITE(6,103)
      DO 5 I=1,NTEMP0
        T6=TEMP(I)*1.E6
        WRITE(6,102)ALOG10(T6),T6,ALFS(I)
   5  CONTINUE
      STOP
C
 101  FORMAT(/'LOG(T)',5X,'TEMP',10X,'ALF GF')
 102  FORMAT(F5.2,1PE12.2,E15.2)
 103  FORMAT(/'LOG(T)',5X,'TEMP',8X,'SUM ALF GF')
 104  FORMAT(//' Z=',F4.1,5X,'DE=',F10.4,' RYD',5X,'FUP=',F10.4,5X
     X,'NTEMP=',I4,5X,'NMIN=',I4,5X,'NMAX=',I4)
 105  FORMAT(/'H-LIKE CORRECTION FACTOR: ',F7.3)
 106  FORMAT(///'FINAL SUMMED OUTPUT FOR CASE: ',A80)
C
      END
